1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This source code is part of
8 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
9 * Copyright (c) 2001-2009, The GROMACS Development Team
11 * Gromacs is a library for molecular simulation and trajectory analysis,
12 * written by Erik Lindahl, David van der Spoel, Berk Hess, and others - for
13 * a full list of developers and information, check out http://www.gromacs.org
15 * This program is free software; you can redistribute it and/or modify it under
16 * the terms of the GNU Lesser General Public License as published by the Free
17 * Software Foundation; either version 2 of the License, or (at your option) any
19 * As a special exception, you may use this file as part of a free software
20 * library without restriction. Specifically, if other files instantiate
21 * templates or use macros or inline functions from this file, or you compile
22 * this file and link it with other files to produce an executable, this
23 * file does not by itself cause the resulting executable to be covered by
24 * the GNU Lesser General Public License.
26 * In plain-speak: do not worry about classes/macros/templates either - only
27 * changes to the library have to be LGPL, not an application linking with it.
29 * To help fund GROMACS development, we humbly ask that you cite
30 * the papers people have written on it - you can find them on the website!
33 /* When calculating RF or Ewald interactions we calculate the electrostatic
34 * forces and energies on excluded atom pairs here in the non-bonded loops.
36 #if defined CHECK_EXCLS && defined CALC_COULOMB
50 egp_cj
= nbat
->energrp
[cj
];
52 for(i
=0; i
<UNROLLI
; i
++)
60 type_i_off
= type
[ai
]*ntype2
;
62 for(j
=0; j
<UNROLLJ
; j
++)
69 real FrLJ6
=0,FrLJ12
=0,VLJ
=0;
85 /* A multiply mask used to zero an interaction
86 * when either the distance cutoff is exceeded, or
87 * (if appropriate) the i and j indices are
88 * unsuitable for this kind of inner loop. */
90 #ifdef VDW_CUTOFF_CHECK
94 /* A multiply mask used to zero an interaction
95 * when that interaction should be excluded
96 * (e.g. because of bonding). */
99 interact
= ((l_cj
[cjind
].excl
>>(i
*UNROLLI
+ j
)) & 1);
103 skipmask
= !(cj
== ci_sh
&& j
<= i
);
112 dx
= xi
[i
*XI_STRIDE
+XX
] - x
[aj
*X_STRIDE
+XX
];
113 dy
= xi
[i
*XI_STRIDE
+YY
] - x
[aj
*X_STRIDE
+YY
];
114 dz
= xi
[i
*XI_STRIDE
+ZZ
] - x
[aj
*X_STRIDE
+ZZ
];
116 rsq
= dx
*dx
+ dy
*dy
+ dz
*dz
;
118 /* Prepare to enforce the cut-off. */
119 skipmask
= (rsq
>= rcut2
) ? 0 : skipmask
;
120 /* 9 flops for r^2 + cut-off check */
123 /* Excluded atoms are allowed to be on top of each other.
124 * To avoid overflow of rinv, rinvsq and rinvsix
125 * we add a small number to rsq for excluded pairs only.
127 rsq
+= (1 - interact
)*NBNXN_AVOID_SING_R2_INC
;
134 rinv
= gmx_invsqrt(rsq
);
135 /* 5 flops for invsqrt */
137 /* Partially enforce the cut-off (and perhaps
138 * exclusions) to avoid possible overflow of
139 * rinvsix when computing LJ, and/or overflowing
140 * the Coulomb table during lookup. */
141 rinv
= rinv
* skipmask
;
149 rinvsix
= interact
*rinvsq
*rinvsq
*rinvsq
;
151 #ifdef VDW_CUTOFF_CHECK
152 skipmask_rvdw
= (rsq
< rvdw2
);
153 rinvsix
*= skipmask_rvdw
;
156 c6
= nbfp
[type_i_off
+type
[aj
]*2 ];
157 c12
= nbfp
[type_i_off
+type
[aj
]*2+1];
159 FrLJ12
= c12
*rinvsix
*rinvsix
;
160 /* 6 flops for r^-2 + LJ force */
162 VLJ
= (FrLJ12
- c12
*sh_invrc6
*sh_invrc6
)/12 -
163 (FrLJ6
- c6
*sh_invrc6
)/6;
164 /* Need to zero the interaction if r >= rcut
165 * or there should be exclusion. */
166 VLJ
= VLJ
* skipmask
* interact
;
167 /* 9 flops for LJ energy */
168 #ifdef VDW_CUTOFF_CHECK
169 VLJ
*= skipmask_rvdw
;
172 Vvdw
[egp_sh_i
[i
]+((egp_cj
>>(nbat
->neg_2log
*j
)) & egp_mask
)] += VLJ
;
175 /* 1 flop for LJ energy addition */
181 /* Enforce the cut-off and perhaps exclusions. In
182 * those cases, rinv is zero because of skipmask,
183 * but fcoul and vcoul will later be non-zero (in
184 * both RF and table cases) because of the
185 * contributions that do not depend on rinv. These
186 * contributions cannot be allowed to accumulate
187 * to the force and potential, and the easiest way
188 * to do this is to zero the charges in
190 qq
= skipmask
* qi
[i
] * q
[aj
];
193 fcoul
= qq
*(interact
*rinv
*rinvsq
- k_rf2
);
194 /* 4 flops for RF force */
196 vcoul
= qq
*(interact
*rinv
+ k_rf
*rsq
- c_rf
);
197 /* 4 flops for RF energy */
202 rs
= rsq
*rinv
*ic
->tabq_scale
;
206 /* fexcl = F_i + frac * (F_(i+1)-F_i) */
207 fexcl
= tab_coul_FDV0
[ri
*4] + frac
*tab_coul_FDV0
[ri
*4+1];
209 /* fexcl = (1-frac) * F_i + frac * F_(i+1) */
210 fexcl
= (1 - frac
)*tab_coul_F
[ri
] + frac
*tab_coul_F
[ri
+1];
212 fcoul
= interact
*rinvsq
- fexcl
;
213 /* 7 flops for float 1/r-table force */
216 vcoul
= qq
*(interact
*(rinv
- ic
->sh_ewald
)
217 -(tab_coul_FDV0
[ri
*4+2]
218 -halfsp
*frac
*(tab_coul_FDV0
[ri
*4] + fexcl
)));
219 /* 7 flops for float 1/r-table energy (8 with excls) */
221 vcoul
= qq
*(interact
*(rinv
- ic
->sh_ewald
)
223 -halfsp
*frac
*(tab_coul_F
[ri
] + fexcl
)));
231 Vc
[egp_sh_i
[i
]+((egp_cj
>>(nbat
->neg_2log
*j
)) & egp_mask
)] += vcoul
;
234 /* 1 flop for Coulomb energy addition */
244 fscal
= (FrLJ12
- FrLJ6
)*rinvsq
+ fcoul
;
245 /* 3 flops for scalar LJ+Coulomb force */
254 fscal
= (FrLJ12
- FrLJ6
)*rinvsq
;
260 /* Increment i-atom force */
261 fi
[i
*FI_STRIDE
+XX
] += fx
;
262 fi
[i
*FI_STRIDE
+YY
] += fy
;
263 fi
[i
*FI_STRIDE
+ZZ
] += fz
;
264 /* Decrement j-atom force */
265 f
[aj
*F_STRIDE
+XX
] -= fx
;
266 f
[aj
*F_STRIDE
+YY
] -= fy
;
267 f
[aj
*F_STRIDE
+ZZ
] -= fz
;
268 /* 9 flops for force addition */