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-2009, The GROMACS Development Team
6 * Copyright (c) 2012,2013, by the GROMACS development team, led by
7 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
8 * others, as listed in the AUTHORS file in the top-level source
9 * 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.
38 /* When calculating RF or Ewald interactions we calculate the electrostatic
39 * forces and energies on excluded atom pairs here in the non-bonded loops.
41 #if defined CHECK_EXCLS && defined CALC_COULOMB
55 egp_cj = nbat->energrp[cj];
57 for(i=0; i<UNROLLI; i++)
65 type_i_off = type[ai]*ntype2;
67 for(j=0; j<UNROLLJ; j++)
74 real FrLJ6=0,FrLJ12=0,VLJ=0;
90 /* A multiply mask used to zero an interaction
91 * when either the distance cutoff is exceeded, or
92 * (if appropriate) the i and j indices are
93 * unsuitable for this kind of inner loop. */
95 #ifdef VDW_CUTOFF_CHECK
99 /* A multiply mask used to zero an interaction
100 * when that interaction should be excluded
101 * (e.g. because of bonding). */
104 interact = ((l_cj[cjind].excl>>(i*UNROLLI + j)) & 1);
108 skipmask = !(cj == ci_sh && j <= i);
117 dx = xi[i*XI_STRIDE+XX] - x[aj*X_STRIDE+XX];
118 dy = xi[i*XI_STRIDE+YY] - x[aj*X_STRIDE+YY];
119 dz = xi[i*XI_STRIDE+ZZ] - x[aj*X_STRIDE+ZZ];
121 rsq = dx*dx + dy*dy + dz*dz;
123 /* Prepare to enforce the cut-off. */
124 skipmask = (rsq >= rcut2) ? 0 : skipmask;
125 /* 9 flops for r^2 + cut-off check */
128 /* Excluded atoms are allowed to be on top of each other.
129 * To avoid overflow of rinv, rinvsq and rinvsix
130 * we add a small number to rsq for excluded pairs only.
132 rsq += (1 - interact)*NBNXN_AVOID_SING_R2_INC;
139 rinv = gmx_invsqrt(rsq);
140 /* 5 flops for invsqrt */
142 /* Partially enforce the cut-off (and perhaps
143 * exclusions) to avoid possible overflow of
144 * rinvsix when computing LJ, and/or overflowing
145 * the Coulomb table during lookup. */
146 rinv = rinv * skipmask;
154 rinvsix = interact*rinvsq*rinvsq*rinvsq;
156 #ifdef VDW_CUTOFF_CHECK
157 skipmask_rvdw = (rsq < rvdw2);
158 rinvsix *= skipmask_rvdw;
161 c6 = nbfp[type_i_off+type[aj]*2 ];
162 c12 = nbfp[type_i_off+type[aj]*2+1];
164 FrLJ12 = c12*rinvsix*rinvsix;
165 /* 6 flops for r^-2 + LJ force */
167 VLJ = (FrLJ12 - c12*sh_invrc6*sh_invrc6)/12 -
168 (FrLJ6 - c6*sh_invrc6)/6;
169 /* Need to zero the interaction if r >= rcut
170 * or there should be exclusion. */
171 VLJ = VLJ * skipmask * interact;
172 /* 9 flops for LJ energy */
173 #ifdef VDW_CUTOFF_CHECK
174 VLJ *= skipmask_rvdw;
177 Vvdw[egp_sh_i[i]+((egp_cj>>(nbat->neg_2log*j)) & egp_mask)] += VLJ;
180 /* 1 flop for LJ energy addition */
186 /* Enforce the cut-off and perhaps exclusions. In
187 * those cases, rinv is zero because of skipmask,
188 * but fcoul and vcoul will later be non-zero (in
189 * both RF and table cases) because of the
190 * contributions that do not depend on rinv. These
191 * contributions cannot be allowed to accumulate
192 * to the force and potential, and the easiest way
193 * to do this is to zero the charges in
195 qq = skipmask * qi[i] * q[aj];
198 fcoul = qq*(interact*rinv*rinvsq - k_rf2);
199 /* 4 flops for RF force */
201 vcoul = qq*(interact*rinv + k_rf*rsq - c_rf);
202 /* 4 flops for RF energy */
207 rs = rsq*rinv*ic->tabq_scale;
211 /* fexcl = F_i + frac * (F_(i+1)-F_i) */
212 fexcl = tab_coul_FDV0[ri*4] + frac*tab_coul_FDV0[ri*4+1];
214 /* fexcl = (1-frac) * F_i + frac * F_(i+1) */
215 fexcl = (1 - frac)*tab_coul_F[ri] + frac*tab_coul_F[ri+1];
217 fcoul = interact*rinvsq - fexcl;
218 /* 7 flops for float 1/r-table force */
221 vcoul = qq*(interact*(rinv - ic->sh_ewald)
222 -(tab_coul_FDV0[ri*4+2]
223 -halfsp*frac*(tab_coul_FDV0[ri*4] + fexcl)));
224 /* 7 flops for float 1/r-table energy (8 with excls) */
226 vcoul = qq*(interact*(rinv - ic->sh_ewald)
228 -halfsp*frac*(tab_coul_F[ri] + fexcl)));
236 Vc[egp_sh_i[i]+((egp_cj>>(nbat->neg_2log*j)) & egp_mask)] += vcoul;
239 /* 1 flop for Coulomb energy addition */
249 fscal = (FrLJ12 - FrLJ6)*rinvsq + fcoul;
250 /* 3 flops for scalar LJ+Coulomb force */
259 fscal = (FrLJ12 - FrLJ6)*rinvsq;
265 /* Increment i-atom force */
266 fi[i*FI_STRIDE+XX] += fx;
267 fi[i*FI_STRIDE+YY] += fy;
268 fi[i*FI_STRIDE+ZZ] += fz;
269 /* Decrement j-atom force */
270 f[aj*F_STRIDE+XX] -= fx;
271 f[aj*F_STRIDE+YY] -= fy;
272 f[aj*F_STRIDE+ZZ] -= fz;
273 /* 9 flops for force addition */