2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
37 /* This file is completely threadsafe - keep it that way! */
42 #include "gromacs/math/vectypes.h"
43 #include "gromacs/mdlib/force.h"
44 #include "gromacs/mdlib/gmx_omp_nthreads.h"
45 #include "gromacs/pbcutil/ishift.h"
46 #include "gromacs/pbcutil/mshift.h"
47 #include "gromacs/pbcutil/pbc.h"
59 static void upd_vir(rvec vir
, real dvx
, real dvy
, real dvz
)
66 void calc_vir(int nxf
, rvec x
[], rvec f
[], tensor vir
,
67 gmx_bool bScrewPBC
, matrix box
)
70 double dvxx
= 0, dvxy
= 0, dvxz
= 0, dvyx
= 0, dvyy
= 0, dvyz
= 0, dvzx
= 0, dvzy
= 0, dvzz
= 0;
72 #pragma omp parallel for num_threads(gmx_omp_nthreads_get(emntDefault)) \
74 reduction(+: dvxx, dvxy, dvxz, dvyx, dvyy, dvyz, dvzx, dvzy, dvzz)
75 for (i
= 0; i
< nxf
; i
++)
77 dvxx
+= x
[i
][XX
]*f
[i
][XX
];
78 dvxy
+= x
[i
][XX
]*f
[i
][YY
];
79 dvxz
+= x
[i
][XX
]*f
[i
][ZZ
];
81 dvyx
+= x
[i
][YY
]*f
[i
][XX
];
82 dvyy
+= x
[i
][YY
]*f
[i
][YY
];
83 dvyz
+= x
[i
][YY
]*f
[i
][ZZ
];
85 dvzx
+= x
[i
][ZZ
]*f
[i
][XX
];
86 dvzy
+= x
[i
][ZZ
]*f
[i
][YY
];
87 dvzz
+= x
[i
][ZZ
]*f
[i
][ZZ
];
92 /* We should correct all odd x-shifts, but the range of isx is -2 to 2 */
93 if (isx
== 1 || isx
== -1)
95 dvyy
+= box
[YY
][YY
]*f
[i
][YY
];
96 dvyz
+= box
[YY
][YY
]*f
[i
][ZZ
];
98 dvzy
+= box
[ZZ
][ZZ
]*f
[i
][YY
];
99 dvzz
+= box
[ZZ
][ZZ
]*f
[i
][ZZ
];
104 upd_vir(vir
[XX
], dvxx
, dvxy
, dvxz
);
105 upd_vir(vir
[YY
], dvyx
, dvyy
, dvyz
);
106 upd_vir(vir
[ZZ
], dvzx
, dvzy
, dvzz
);
110 static void lo_fcv(int i0
, int i1
,
111 real x
[], real f
[], tensor vir
,
112 int is
[], real box
[], gmx_bool bTriclinic
)
114 int i
, i3
, tx
, ty
, tz
;
116 real dvxx
= 0, dvxy
= 0, dvxz
= 0, dvyx
= 0, dvyy
= 0, dvyz
= 0, dvzx
= 0, dvzy
= 0, dvzz
= 0;
120 for (i
= i0
; (i
< i1
); i
++)
127 xx
= x
[i3
+XX
]-tx
*box
[XXXX
]-ty
*box
[YYXX
]-tz
*box
[ZZXX
];
132 yy
= x
[i3
+YY
]-ty
*box
[YYYY
]-tz
*box
[ZZYY
];
137 zz
= x
[i3
+ZZ
]-tz
*box
[ZZZZ
];
145 for (i
= i0
; (i
< i1
); i
++)
152 xx
= x
[i3
+XX
]-tx
*box
[XXXX
];
157 yy
= x
[i3
+YY
]-ty
*box
[YYYY
];
162 zz
= x
[i3
+ZZ
]-tz
*box
[ZZZZ
];
169 upd_vir(vir
[XX
], dvxx
, dvxy
, dvxz
);
170 upd_vir(vir
[YY
], dvyx
, dvyy
, dvyz
);
171 upd_vir(vir
[ZZ
], dvzx
, dvzy
, dvzz
);
174 void f_calc_vir(int i0
, int i1
, rvec x
[], rvec f
[], tensor vir
,
175 t_graph
*g
, matrix box
)
179 if (g
&& (g
->nnodes
> 0))
181 /* Calculate virial for bonded forces only when they belong to
184 start
= std::max(i0
, g
->at_start
);
185 end
= std::min(i1
, g
->at_end
);
186 lo_fcv(start
, end
, x
[0], f
[0], vir
, g
->ishift
[0], box
[0], TRICLINIC(box
));
188 /* If not all atoms are bonded, calculate their virial contribution
189 * anyway, without shifting back their coordinates.
190 * Note the nifty pointer arithmetic...
194 calc_vir(start
-i0
, x
+ i0
, f
+ i0
, vir
, FALSE
, box
);
198 calc_vir(i1
-end
, x
+ end
, f
+ end
, vir
, FALSE
, box
);
203 calc_vir(i1
-i0
, x
+ i0
, f
+ i0
, vir
, FALSE
, box
);