Implemented LJ-PME nbnxn kernels
[alexxy/gromacs.git] / src / gromacs / mdlib / nbnxn_kernels / nbnxn_kernel_ref_inner.h
1 /*
2  * This file is part of the GROMACS molecular simulation package.
3  *
4  * Copyright (c) 2012,2013,2014, by the GROMACS development team, led by
5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6  * and including many others, as listed in the AUTHORS file in the
7  * top-level source directory and at http://www.gromacs.org.
8  *
9  * GROMACS is free software; you can redistribute it and/or
10  * modify it under the terms of the GNU Lesser General Public License
11  * as published by the Free Software Foundation; either version 2.1
12  * of the License, or (at your option) any later version.
13  *
14  * GROMACS is distributed in the hope that it will be useful,
15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
17  * Lesser General Public License for more details.
18  *
19  * You should have received a copy of the GNU Lesser General Public
20  * License along with GROMACS; if not, see
21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
23  *
24  * If you want to redistribute modifications to GROMACS, please
25  * consider that scientific software is very special. Version
26  * control is crucial - bugs must be traceable. We will be happy to
27  * consider code for inclusion in the official distribution, but
28  * derived work must not be called official GROMACS. Details are found
29  * in the README & COPYING files - if they are missing, get the
30  * official version at http://www.gromacs.org.
31  *
32  * To help us fund GROMACS development, we humbly ask that you cite
33  * the research papers on the package. Check out http://www.gromacs.org.
34  */
35
36 /* When calculating RF or Ewald interactions we calculate the electrostatic
37  * forces and energies on excluded atom pairs here in the non-bonded loops.
38  */
39 #if defined CHECK_EXCLS && (defined CALC_COULOMB || defined LJ_EWALD)
40 #define EXCL_FORCES
41 #endif
42
43 {
44     int cj;
45 #ifdef ENERGY_GROUPS
46     int egp_cj;
47 #endif
48     int i;
49
50     cj = l_cj[cjind].cj;
51
52 #ifdef ENERGY_GROUPS
53     egp_cj = nbat->energrp[cj];
54 #endif
55     for (i = 0; i < UNROLLI; i++)
56     {
57         int ai;
58         int type_i_off;
59         int j;
60
61         ai = ci*UNROLLI + i;
62
63         type_i_off = type[ai]*ntype2;
64
65         for (j = 0; j < UNROLLJ; j++)
66         {
67             int  aj;
68             real dx, dy, dz;
69             real rsq, rinv;
70             real rinvsq, rinvsix;
71             real c6, c12;
72             real FrLJ6 = 0, FrLJ12 = 0, frLJ = 0, VLJ = 0;
73 #if defined LJ_FORCE_SWITCH || defined LJ_POT_SWITCH
74             real r, rsw;
75 #endif
76
77 #ifdef CALC_COULOMB
78             real qq;
79             real fcoul;
80 #ifdef CALC_COUL_TAB
81             real rs, frac;
82             int  ri;
83             real fexcl;
84 #endif
85 #ifdef CALC_ENERGIES
86             real vcoul;
87 #endif
88 #endif
89             real fscal;
90             real fx, fy, fz;
91
92             /* A multiply mask used to zero an interaction
93              * when either the distance cutoff is exceeded, or
94              * (if appropriate) the i and j indices are
95              * unsuitable for this kind of inner loop. */
96             real skipmask;
97
98 #ifdef CHECK_EXCLS
99             /* A multiply mask used to zero an interaction
100              * when that interaction should be excluded
101              * (e.g. because of bonding). */
102             int interact;
103
104             interact = ((l_cj[cjind].excl>>(i*UNROLLI + j)) & 1);
105 #ifndef EXCL_FORCES
106             skipmask = interact;
107 #else
108             skipmask = !(cj == ci_sh && j <= i);
109 #endif
110 #else
111 #define interact 1.0
112             skipmask = 1.0;
113 #endif
114
115             aj = cj*UNROLLJ + j;
116
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];
120
121             rsq = dx*dx + dy*dy + dz*dz;
122
123             /* Prepare to enforce the cut-off. */
124             skipmask = (rsq >= rcut2) ? 0 : skipmask;
125             /* 9 flops for r^2 + cut-off check */
126
127 #ifdef CHECK_EXCLS
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.
131              */
132             rsq += (1 - interact)*NBNXN_AVOID_SING_R2_INC;
133 #endif
134
135 #ifdef COUNT_PAIRS
136             npair++;
137 #endif
138
139             rinv = gmx_invsqrt(rsq);
140             /* 5 flops for invsqrt */
141
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;
147
148             rinvsq  = rinv*rinv;
149
150 #ifdef HALF_LJ
151             if (i < UNROLLI/2)
152 #endif
153             {
154                 c6      = nbfp[type_i_off+type[aj]*2  ];
155                 c12     = nbfp[type_i_off+type[aj]*2+1];
156
157 #if defined LJ_CUT || defined LJ_FORCE_SWITCH || defined LJ_POT_SWITCH
158                 rinvsix = interact*rinvsq*rinvsq*rinvsq;
159                 FrLJ6   = c6*rinvsix;
160                 FrLJ12  = c12*rinvsix*rinvsix;
161                 frLJ    = FrLJ12 - FrLJ6;
162                 /* 7 flops for r^-2 + LJ force */
163 #if defined CALC_ENERGIES || defined LJ_POT_SWITCH
164                 VLJ     = (FrLJ12 + c12*ic->repulsion_shift.cpot)/12 -
165                     (FrLJ6 + c6*ic->dispersion_shift.cpot)/6;
166                 /* 7 flops for LJ energy */
167 #endif
168 #endif
169
170 #if defined LJ_FORCE_SWITCH || defined LJ_POT_SWITCH
171                 /* Force or potential switching from ic->rvdw_switch */
172                 r       = rsq*rinv;
173                 rsw     = r - ic->rvdw_switch;
174                 rsw     = (rsw >= 0.0 ? rsw : 0.0);
175 #endif
176 #ifdef LJ_FORCE_SWITCH
177                 frLJ   +=
178                     -c6*(ic->dispersion_shift.c2 + ic->dispersion_shift.c3*rsw)*rsw*rsw*r
179                     + c12*(ic->repulsion_shift.c2 + ic->repulsion_shift.c3*rsw)*rsw*rsw*r;
180 #if defined CALC_ENERGIES
181                 VLJ    +=
182                     -c6*(-ic->dispersion_shift.c2/3 - ic->dispersion_shift.c3/4*rsw)*rsw*rsw*rsw
183                     + c12*(-ic->repulsion_shift.c2/3 - ic->repulsion_shift.c3/4*rsw)*rsw*rsw*rsw;
184 #endif
185 #endif
186
187 #if defined CALC_ENERGIES || defined LJ_POT_SWITCH
188                 /* Masking should be done after force switching,
189                  * but before potential switching.
190                  */
191                 /* Need to zero the interaction if there should be exclusion. */
192                 VLJ     = VLJ * interact;
193 #endif
194
195 #ifdef LJ_POT_SWITCH
196                 {
197                     real sw, dsw;
198
199                     sw    = 1.0 + (swV3 + (swV4+ swV5*rsw)*rsw)*rsw*rsw*rsw;
200                     dsw   = (swF2 + (swF3 + swF4*rsw)*rsw)*rsw*rsw;
201
202                     frLJ  = frLJ*sw - r*VLJ*dsw;
203                     VLJ  *= sw;
204                 }
205 #endif
206
207 #ifdef LJ_EWALD
208                 {
209                     real c6grid, rinvsix_nm, cr2, expmcr2, poly, sh_mask;
210
211 #ifdef LJ_EWALD_COMB_GEOM
212                     c6grid       = ljc[type[ai]*2]*ljc[type[aj]*2];
213 #elif defined LJ_EWALD_COMB_LB
214                     {
215                         real sigma, sigma2, epsilon;
216
217                         /* These sigma and epsilon are scaled to give 6*C6 */
218                         sigma   = ljc[type[ai]*2] + ljc[type[aj]*2];
219                         epsilon = ljc[type[ai]*2+1]*ljc[type[aj]*2+1];
220
221                         sigma2  = sigma*sigma;
222                         c6grid  = epsilon*sigma2*sigma2*sigma2;
223                     }
224 #else
225 #error "No LJ Ewald combination rule defined"
226 #endif
227
228 #ifdef CHECK_EXCLS
229                     /* Recalculate rinvsix without exclusion mask */
230                     rinvsix_nm   = rinvsq*rinvsq*rinvsq;
231 #else
232                     rinvsix_nm   = rinvsix;
233 #endif
234                     cr2          = lje_coeff2*rsq;
235 #ifdef GMX_DOUBLE
236                     expmcr2      = exp(-cr2);
237 #else
238                     expmcr2      = expf(-cr2);
239 #endif
240                     poly         = 1 + cr2 + 0.5*cr2*cr2;
241
242                     /* Subtract the grid force from the total LJ force */
243                     frLJ        += c6grid*(rinvsix_nm - expmcr2*(rinvsix_nm*poly + lje_coeff6_6));
244 #ifdef CALC_ENERGIES
245                     /* Shift should only be applied to real LJ pairs */
246                     sh_mask      = lje_vc*interact;
247
248                     VLJ         += c6grid/6*(rinvsix_nm*(1 - expmcr2*poly) + sh_mask);
249 #endif
250                 }
251 #endif          /* LJ_EWALD */
252
253 #ifdef VDW_CUTOFF_CHECK
254                 /* Mask for VdW cut-off shorter than Coulomb cut-off */
255                 {
256                     real skipmask_rvdw;
257
258                     skipmask_rvdw = (rsq < rvdw2);
259                     frLJ         *= skipmask_rvdw;
260 #ifdef CALC_ENERGIES
261                     VLJ    *= skipmask_rvdw;
262 #endif
263                 }
264 #else
265 #if defined CALC_ENERGIES
266                 /* Need to zero the interaction if r >= rcut */
267                 VLJ     = VLJ * skipmask;
268                 /* 1 more flop for LJ energy */
269 #endif
270 #endif          /* VDW_CUTOFF_CHECK */
271
272
273 #ifdef CALC_ENERGIES
274 #ifdef ENERGY_GROUPS
275                 Vvdw[egp_sh_i[i]+((egp_cj>>(nbat->neg_2log*j)) & egp_mask)] += VLJ;
276 #else
277                 Vvdw_ci += VLJ;
278                 /* 1 flop for LJ energy addition */
279 #endif
280 #endif
281             }
282
283 #ifdef CALC_COULOMB
284             /* Enforce the cut-off and perhaps exclusions. In
285              * those cases, rinv is zero because of skipmask,
286              * but fcoul and vcoul will later be non-zero (in
287              * both RF and table cases) because of the
288              * contributions that do not depend on rinv. These
289              * contributions cannot be allowed to accumulate
290              * to the force and potential, and the easiest way
291              * to do this is to zero the charges in
292              * advance. */
293             qq = skipmask * qi[i] * q[aj];
294
295 #ifdef CALC_COUL_RF
296             fcoul  = qq*(interact*rinv*rinvsq - k_rf2);
297             /* 4 flops for RF force */
298 #ifdef CALC_ENERGIES
299             vcoul  = qq*(interact*rinv + k_rf*rsq - c_rf);
300             /* 4 flops for RF energy */
301 #endif
302 #endif
303
304 #ifdef CALC_COUL_TAB
305             rs     = rsq*rinv*ic->tabq_scale;
306             ri     = (int)rs;
307             frac   = rs - ri;
308 #ifndef GMX_DOUBLE
309             /* fexcl = F_i + frac * (F_(i+1)-F_i) */
310             fexcl  = tab_coul_FDV0[ri*4] + frac*tab_coul_FDV0[ri*4+1];
311 #else
312             /* fexcl = (1-frac) * F_i + frac * F_(i+1) */
313             fexcl  = (1 - frac)*tab_coul_F[ri] + frac*tab_coul_F[ri+1];
314 #endif
315             fcoul  = interact*rinvsq - fexcl;
316             /* 7 flops for float 1/r-table force */
317 #ifdef CALC_ENERGIES
318 #ifndef GMX_DOUBLE
319             vcoul  = qq*(interact*(rinv - ic->sh_ewald)
320                          -(tab_coul_FDV0[ri*4+2]
321                            -halfsp*frac*(tab_coul_FDV0[ri*4] + fexcl)));
322             /* 7 flops for float 1/r-table energy (8 with excls) */
323 #else
324             vcoul  = qq*(interact*(rinv - ic->sh_ewald)
325                          -(tab_coul_V[ri]
326                            -halfsp*frac*(tab_coul_F[ri] + fexcl)));
327 #endif
328 #endif
329             fcoul *= qq*rinv;
330 #endif
331
332 #ifdef CALC_ENERGIES
333 #ifdef ENERGY_GROUPS
334             Vc[egp_sh_i[i]+((egp_cj>>(nbat->neg_2log*j)) & egp_mask)] += vcoul;
335 #else
336             Vc_ci += vcoul;
337             /* 1 flop for Coulomb energy addition */
338 #endif
339 #endif
340 #endif
341
342 #ifdef CALC_COULOMB
343 #ifdef HALF_LJ
344             if (i < UNROLLI/2)
345 #endif
346             {
347                 fscal = frLJ*rinvsq + fcoul;
348                 /* 2 flops for scalar LJ+Coulomb force */
349             }
350 #ifdef HALF_LJ
351             else
352             {
353                 fscal = fcoul;
354             }
355 #endif
356 #else
357             fscal = frLJ*rinvsq;
358 #endif
359             fx = fscal*dx;
360             fy = fscal*dy;
361             fz = fscal*dz;
362
363             /* Increment i-atom force */
364             fi[i*FI_STRIDE+XX] += fx;
365             fi[i*FI_STRIDE+YY] += fy;
366             fi[i*FI_STRIDE+ZZ] += fz;
367             /* Decrement j-atom force */
368             f[aj*F_STRIDE+XX]  -= fx;
369             f[aj*F_STRIDE+YY]  -= fy;
370             f[aj*F_STRIDE+ZZ]  -= fz;
371             /* 9 flops for force addition */
372         }
373     }
374 }
375
376 #undef interact
377 #undef EXCL_FORCES