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 upd_vir(rvec vir,real dvx,real dvy,real dvz)
53 void calc_vir(FILE *log,int nxf,rvec x[],rvec f[],tensor vir,
54 gmx_bool bScrewPBC,matrix box)
57 double dvxx=0,dvxy=0,dvxz=0,dvyx=0,dvyy=0,dvyz=0,dvzx=0,dvzy=0,dvzz=0;
59 for(i=0; (i<nxf); i++) {
60 dvxx+=x[i][XX]*f[i][XX];
61 dvxy+=x[i][XX]*f[i][YY];
62 dvxz+=x[i][XX]*f[i][ZZ];
64 dvyx+=x[i][YY]*f[i][XX];
65 dvyy+=x[i][YY]*f[i][YY];
66 dvyz+=x[i][YY]*f[i][ZZ];
68 dvzx+=x[i][ZZ]*f[i][XX];
69 dvzy+=x[i][ZZ]*f[i][YY];
70 dvzz+=x[i][ZZ]*f[i][ZZ];
74 /* We should correct all odd x-shifts, but the range of isx is -2 to 2 */
75 if (isx == 1 || isx == -1) {
76 dvyy += box[YY][YY]*f[i][YY];
77 dvyz += box[YY][YY]*f[i][ZZ];
79 dvzy += box[ZZ][ZZ]*f[i][YY];
80 dvzz += box[ZZ][ZZ]*f[i][ZZ];
85 upd_vir(vir[XX],dvxx,dvxy,dvxz);
86 upd_vir(vir[YY],dvyx,dvyy,dvyz);
87 upd_vir(vir[ZZ],dvzx,dvzy,dvzz);
91 static void lo_fcv(int i0,int i1,
92 real x[],real f[],tensor vir,
93 int is[],real box[], gmx_bool bTriclinic)
97 real dvxx=0,dvxy=0,dvxz=0,dvyx=0,dvyy=0,dvyz=0,dvzx=0,dvzy=0,dvzz=0;
100 for(i=i0; (i<i1); i++) {
106 xx=x[i3+XX]-tx*box[XXXX]-ty*box[YYXX]-tz*box[ZZXX];
111 yy=x[i3+YY]-ty*box[YYYY]-tz*box[ZZYY];
116 zz=x[i3+ZZ]-tz*box[ZZZZ];
122 for(i=i0; (i<i1); i++) {
128 xx=x[i3+XX]-tx*box[XXXX];
133 yy=x[i3+YY]-ty*box[YYYY];
138 zz=x[i3+ZZ]-tz*box[ZZZZ];
145 upd_vir(vir[XX],dvxx,dvxy,dvxz);
146 upd_vir(vir[YY],dvyx,dvyy,dvyz);
147 upd_vir(vir[ZZ],dvzx,dvzy,dvzz);
150 static void lo_fcv2(int i0,int i1,
151 rvec x[],rvec f[],tensor vir,
152 ivec is[],matrix box, gmx_bool bTriclinic)
156 real dvxx=0,dvxy=0,dvxz=0,dvyx=0,dvyy=0,dvyz=0,dvzx=0,dvzy=0,dvzz=0;
159 for(i=i0,gg=0; (i<i1); i++,gg++) {
164 xx=x[i][XX]-tx*box[XX][XX]-ty*box[YY][XX]-tz*box[ZZ][XX];
169 yy=x[i][YY]-ty*box[YY][YY]-tz*box[ZZ][YY];
174 zz=x[i][ZZ]-tz*box[ZZ][ZZ];
180 for(i=i0,gg=0; (i<i1); i++,gg++) {
185 xx=x[i][XX]-tx*box[XX][XX];
190 yy=x[i][YY]-ty*box[YY][YY];
195 zz=x[i][ZZ]-tz*box[ZZ][ZZ];
202 upd_vir(vir[XX],dvxx,dvxy,dvxz);
203 upd_vir(vir[YY],dvyx,dvyy,dvyz);
204 upd_vir(vir[ZZ],dvzx,dvzy,dvzz);
207 void f_calc_vir(FILE *log,int i0,int i1,rvec x[],rvec f[],tensor vir,
208 t_graph *g,matrix box)
212 if (g && (g->nnodes > 0)) {
213 /* Calculate virial for bonded forces only when they belong to
216 start = max(i0,g->at_start);
217 end = min(i1,g->at_end);
219 lo_fcv2(start,end,x,f,vir,g->ishift,box,TRICLINIC(box));
221 lo_fcv(start,end,x[0],f[0],vir,g->ishift[0],box[0],TRICLINIC(box));
224 /* If not all atoms are bonded, calculate their virial contribution
225 * anyway, without shifting back their coordinates.
226 * Note the nifty pointer arithmetic...
229 calc_vir(log,start-i0,x + i0,f + i0,vir,FALSE,box);
231 calc_vir(log,i1-end,x + end,f + end,vir,FALSE,box);
234 calc_vir(log,i1-i0,x + i0,f + i0,vir,FALSE,box);