3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * GROwing Monsters And Cloning Shrimps
35 /* This file is completely threadsafe - keep it that way! */
46 static void dprod1(rvec vir,real x,rvec f)
53 static void upd_vir(rvec vir,real dvx,real dvy,real dvz)
60 void calc_vir(FILE *log,int nxf,rvec x[],rvec f[],tensor vir,
61 bool bScrewPBC,matrix box)
64 double dvxx=0,dvxy=0,dvxz=0,dvyx=0,dvyy=0,dvyz=0,dvzx=0,dvzy=0,dvzz=0;
66 for(i=0; (i<nxf); i++) {
67 dvxx+=x[i][XX]*f[i][XX];
68 dvxy+=x[i][XX]*f[i][YY];
69 dvxz+=x[i][XX]*f[i][ZZ];
71 dvyx+=x[i][YY]*f[i][XX];
72 dvyy+=x[i][YY]*f[i][YY];
73 dvyz+=x[i][YY]*f[i][ZZ];
75 dvzx+=x[i][ZZ]*f[i][XX];
76 dvzy+=x[i][ZZ]*f[i][YY];
77 dvzz+=x[i][ZZ]*f[i][ZZ];
81 /* We should correct all odd x-shifts, but the range of isx is -2 to 2 */
82 if (isx == 1 || isx == -1) {
83 dvyy += box[YY][YY]*f[i][YY];
84 dvyz += box[YY][YY]*f[i][ZZ];
86 dvzy += box[ZZ][ZZ]*f[i][YY];
87 dvzz += box[ZZ][ZZ]*f[i][ZZ];
92 upd_vir(vir[XX],dvxx,dvxy,dvxz);
93 upd_vir(vir[YY],dvyx,dvyy,dvyz);
94 upd_vir(vir[ZZ],dvzx,dvzy,dvzz);
98 static void lo_fcv(int i0,int i1,int g0,
99 real x[],real f[],tensor vir,
100 int is[],real box[], bool bTriclinic)
102 int i,i3,gg,g3,tx,ty,tz;
104 real dvxx=0,dvxy=0,dvxz=0,dvyx=0,dvyy=0,dvyz=0,dvzx=0,dvzy=0,dvzz=0;
107 for(i=i0,gg=g0; (i<i1); i++,gg++) {
114 xx=x[i3+XX]-tx*box[XXXX]-ty*box[YYXX]-tz*box[ZZXX];
119 yy=x[i3+YY]-ty*box[YYYY]-tz*box[ZZYY];
124 zz=x[i3+ZZ]-tz*box[ZZZZ];
130 for(i=i0,gg=g0; (i<i1); i++,gg++) {
137 xx=x[i3+XX]-tx*box[XXXX];
142 yy=x[i3+YY]-ty*box[YYYY];
147 zz=x[i3+ZZ]-tz*box[ZZZZ];
154 upd_vir(vir[XX],dvxx,dvxy,dvxz);
155 upd_vir(vir[YY],dvyx,dvyy,dvyz);
156 upd_vir(vir[ZZ],dvzx,dvzy,dvzz);
159 static void lo_fcv2(int i0,int i1,
160 rvec x[],rvec f[],tensor vir,
161 ivec is[],matrix box, bool bTriclinic)
165 real dvxx=0,dvxy=0,dvxz=0,dvyx=0,dvyy=0,dvyz=0,dvzx=0,dvzy=0,dvzz=0;
168 for(i=i0,gg=0; (i<i1); i++,gg++) {
173 xx=x[i][XX]-tx*box[XX][XX]-ty*box[YY][XX]-tz*box[ZZ][XX];
178 yy=x[i][YY]-ty*box[YY][YY]-tz*box[ZZ][YY];
183 zz=x[i][ZZ]-tz*box[ZZ][ZZ];
189 for(i=i0,gg=0; (i<i1); i++,gg++) {
194 xx=x[i][XX]-tx*box[XX][XX];
199 yy=x[i][YY]-ty*box[YY][YY];
204 zz=x[i][ZZ]-tz*box[ZZ][ZZ];
211 upd_vir(vir[XX],dvxx,dvxy,dvxz);
212 upd_vir(vir[YY],dvyx,dvyy,dvyz);
213 upd_vir(vir[ZZ],dvzx,dvzy,dvzz);
216 void f_calc_vir(FILE *log,int i0,int i1,rvec x[],rvec f[],tensor vir,
217 t_graph *g,matrix box)
221 if (g && (g->nnodes > 0)) {
222 /* Calculate virial for bonded forces only when they belong to
225 start = max(i0,g->start);
226 end = min(i1,g->end+1);
228 lo_fcv2(start,end,x,f,vir,g->ishift,box,TRICLINIC(box));
230 lo_fcv(start,end,0,x[0],f[0],vir,g->ishift[0],box[0],TRICLINIC(box));
233 /* If not all atoms are bonded, calculate their virial contribution
234 * anyway, without shifting back their coordinates.
235 * Note the nifty pointer arithmetic...
238 calc_vir(log,start-i0,x + i0,f + i0,vir,FALSE,box);
240 calc_vir(log,i1-end,x + end,f + end,vir,FALSE,box);
243 calc_vir(log,i1-i0,x + i0,f + i0,vir,FALSE,box);