2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2012,2013,2014,2015,2016,2017,2018, 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.
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.
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.
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.
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.
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.
38 const nbnxn_ci_t *nbln;
39 const nbnxn_cj_t *l_cj;
46 gmx_bool do_LJ, half_LJ, do_coul;
47 int cjind0, cjind1, cjind;
51 int egps_ishift, egps_imask;
52 int egps_jshift, egps_jmask, egps_jstride;
54 real *vvdwtp[UNROLLI];
61 SimdReal ix_S0, iy_S0, iz_S0;
62 SimdReal ix_S2, iy_S2, iz_S2;
63 SimdReal fix_S0, fiy_S0, fiz_S0;
64 SimdReal fix_S2, fiy_S2, fiz_S2;
66 SimdReal diagonal_jmi_S;
67 #if UNROLLI == UNROLLJ
68 SimdBool diagonal_mask_S0, diagonal_mask_S2;
70 SimdBool diagonal_mask0_S0, diagonal_mask0_S2;
71 SimdBool diagonal_mask1_S0, diagonal_mask1_S2;
74 unsigned *exclusion_filter;
75 SimdBitMask filter_S0, filter_S2;
80 SimdReal iq_S0 = setZero();
81 SimdReal iq_S2 = setZero();
86 SimdReal hrc_3_S, moh_rc_S;
91 /* Coulomb table variables */
93 const real *tab_coul_F;
94 #if defined CALC_ENERGIES && !defined TAB_FDV0
95 const real *tab_coul_V;
103 #ifdef CALC_COUL_EWALD
104 SimdReal beta2_S, beta_S;
107 #if defined CALC_ENERGIES && (defined CALC_COUL_EWALD || defined CALC_COUL_TAB)
111 #if defined LJ_CUT && defined CALC_ENERGIES
112 SimdReal p6_cpot_S, p12_cpot_S;
116 SimdReal swV3_S, swV4_S, swV5_S;
117 SimdReal swF2_S, swF3_S, swF4_S;
119 #ifdef LJ_FORCE_SWITCH
121 SimdReal p6_fc2_S, p6_fc3_S;
122 SimdReal p12_fc2_S, p12_fc3_S;
124 SimdReal p6_vc3_S, p6_vc4_S;
125 SimdReal p12_vc3_S, p12_vc4_S;
126 SimdReal p6_6cpot_S, p12_12cpot_S;
130 real lj_ewaldcoeff2, lj_ewaldcoeff6_6;
131 SimdReal mone_S, half_S, lje_c2_S, lje_c6_6_S;
137 SimdReal hsig_i_S0, seps_i_S0;
138 SimdReal hsig_i_S2, seps_i_S2;
141 alignas(GMX_SIMD_ALIGNMENT) real pvdw_c6[2*UNROLLI*UNROLLJ];
142 real *pvdw_c12 = pvdw_c6 + UNROLLI*UNROLLJ;
145 #if defined LJ_COMB_GEOM || defined LJ_EWALD_GEOM
148 #endif /* LJ_COMB_LB */
152 #ifdef VDW_CUTOFF_CHECK
162 #if defined LJ_COMB_GEOM || defined LJ_COMB_LB || defined LJ_EWALD_GEOM
165 #if !(defined LJ_COMB_GEOM || defined LJ_COMB_LB || defined FIX_LJ_C)
166 /* No combination rule used */
167 real *nbfp_ptr = nbat->nbfp_aligned;
168 const int *type = nbat->type;
171 /* Load j-i for the first i */
172 diagonal_jmi_S = load<SimdReal>(nbat->simd_2xnn_diagonal_j_minus_i);
173 /* Generate all the diagonal masks as comparison results */
174 #if UNROLLI == UNROLLJ
175 diagonal_mask_S0 = (zero_S < diagonal_jmi_S);
176 diagonal_jmi_S = diagonal_jmi_S - one_S;
177 diagonal_jmi_S = diagonal_jmi_S - one_S;
178 diagonal_mask_S2 = (zero_S < diagonal_jmi_S);
180 #if 2*UNROLLI == UNROLLJ
181 diagonal_mask0_S0 = (zero_S < diagonal_jmi_S);
182 diagonal_jmi_S = diagonal_jmi_S - one_S;
183 diagonal_jmi_S = diagonal_jmi_S - one_S;
184 diagonal_mask0_S2 = (zero_S < diagonal_jmi_S);
185 diagonal_jmi_S = diagonal_jmi_S - one_S;
186 diagonal_jmi_S = diagonal_jmi_S - one_S;
187 diagonal_mask1_S0 = (zero_S < diagonal_jmi_S);
188 diagonal_jmi_S = diagonal_jmi_S - one_S;
189 diagonal_jmi_S = diagonal_jmi_S - one_S;
190 diagonal_mask1_S2 = (zero_S < diagonal_jmi_S);
194 /* Load masks for topology exclusion masking. filter_stride is
195 static const, so the conditional will be optimized away. */
196 #if GMX_DOUBLE && !GMX_SIMD_HAVE_INT32_LOGICAL
197 exclusion_filter = nbat->simd_exclusion_filter64;
199 exclusion_filter = nbat->simd_exclusion_filter;
202 /* Here we cast the exclusion filters from unsigned * to int * or real *.
203 * Since we only check bits, the actual value they represent does not
204 * matter, as long as both filter and mask data are treated the same way.
206 #if GMX_SIMD_HAVE_INT32_LOGICAL
207 filter_S0 = load<SimdBitMask>(reinterpret_cast<const int *>(exclusion_filter + 0*UNROLLJ));
208 filter_S2 = load<SimdBitMask>(reinterpret_cast<const int *>(exclusion_filter + 2*UNROLLJ));
210 filter_S0 = load<SimdBitMask>(reinterpret_cast<const real *>(exclusion_filter + 0*UNROLLJ));
211 filter_S2 = load<SimdBitMask>(reinterpret_cast<const real *>(exclusion_filter + 2*UNROLLJ));
215 /* Reaction-field constants */
216 mrc_3_S = SimdReal(-2*ic->k_rf);
218 hrc_3_S = SimdReal(ic->k_rf);
219 moh_rc_S = SimdReal(-ic->c_rf);
225 invtsp_S = SimdReal(ic->tabq_scale);
227 mhalfsp_S = SimdReal(-0.5/ic->tabq_scale);
231 tab_coul_F = ic->tabq_coul_FDV0;
233 tab_coul_F = ic->tabq_coul_F;
235 tab_coul_V = ic->tabq_coul_V;
238 #endif /* CALC_COUL_TAB */
240 #ifdef CALC_COUL_EWALD
241 beta2_S = SimdReal(ic->ewaldcoeff_q*ic->ewaldcoeff_q);
242 beta_S = SimdReal(ic->ewaldcoeff_q);
245 #if (defined CALC_COUL_TAB || defined CALC_COUL_EWALD) && defined CALC_ENERGIES
246 sh_ewald_S = SimdReal(ic->sh_ewald);
249 /* LJ function constants */
250 #if defined CALC_ENERGIES || defined LJ_POT_SWITCH
251 SimdReal sixth_S = SimdReal(1.0/6.0);
252 SimdReal twelveth_S = SimdReal(1.0/12.0);
255 #if defined LJ_CUT && defined CALC_ENERGIES
256 /* We shift the potential by cpot, which can be zero */
257 p6_cpot_S = SimdReal(ic->dispersion_shift.cpot);
258 p12_cpot_S = SimdReal(ic->repulsion_shift.cpot);
261 rswitch_S = SimdReal(ic->rvdw_switch);
262 swV3_S = SimdReal(ic->vdw_switch.c3);
263 swV4_S = SimdReal(ic->vdw_switch.c4);
264 swV5_S = SimdReal(ic->vdw_switch.c5);
265 swF2_S = SimdReal(3*ic->vdw_switch.c3);
266 swF3_S = SimdReal(4*ic->vdw_switch.c4);
267 swF4_S = SimdReal(5*ic->vdw_switch.c5);
269 #ifdef LJ_FORCE_SWITCH
270 rswitch_S = SimdReal(ic->rvdw_switch);
271 p6_fc2_S = SimdReal(ic->dispersion_shift.c2);
272 p6_fc3_S = SimdReal(ic->dispersion_shift.c3);
273 p12_fc2_S = SimdReal(ic->repulsion_shift.c2);
274 p12_fc3_S = SimdReal(ic->repulsion_shift.c3);
277 SimdReal mthird_S = SimdReal(-1.0/3.0);
278 SimdReal mfourth_S = SimdReal(-1.0/4.0);
280 p6_vc3_S = mthird_S * p6_fc2_S;
281 p6_vc4_S = mfourth_S * p6_fc3_S;
282 p6_6cpot_S = SimdReal(ic->dispersion_shift.cpot/6);
283 p12_vc3_S = mthird_S * p12_fc2_S;
284 p12_vc4_S = mfourth_S * p12_fc3_S;
285 p12_12cpot_S = SimdReal(ic->repulsion_shift.cpot/12);
290 mone_S = SimdReal(-1.0);
291 half_S = SimdReal(0.5);
292 lj_ewaldcoeff2 = ic->ewaldcoeff_lj*ic->ewaldcoeff_lj;
293 lj_ewaldcoeff6_6 = lj_ewaldcoeff2*lj_ewaldcoeff2*lj_ewaldcoeff2/6;
294 lje_c2_S = SimdReal(lj_ewaldcoeff2);
295 lje_c6_6_S = SimdReal(lj_ewaldcoeff6_6);
297 /* Determine the grid potential at the cut-off */
298 SimdReal lje_vc_S = SimdReal(ic->sh_lj_ewald);
302 /* The kernel either supports rcoulomb = rvdw or rcoulomb >= rvdw */
303 rc2_S = SimdReal(ic->rcoulomb*ic->rcoulomb);
304 #ifdef VDW_CUTOFF_CHECK
305 rcvdw2_S = SimdReal(ic->rvdw*ic->rvdw);
308 minRsq_S = SimdReal(NBNXN_MIN_RSQ);
312 shiftvec = shift_vec[0];
317 for (jp = 0; jp < UNROLLJ; jp++)
319 pvdw_c6 [0*UNROLLJ+jp] = nbat->nbfp[0*2];
320 pvdw_c6 [1*UNROLLJ+jp] = nbat->nbfp[0*2];
321 pvdw_c6 [2*UNROLLJ+jp] = nbat->nbfp[0*2];
322 pvdw_c6 [3*UNROLLJ+jp] = nbat->nbfp[0*2];
324 pvdw_c12[0*UNROLLJ+jp] = nbat->nbfp[0*2+1];
325 pvdw_c12[1*UNROLLJ+jp] = nbat->nbfp[0*2+1];
326 pvdw_c12[2*UNROLLJ+jp] = nbat->nbfp[0*2+1];
327 pvdw_c12[3*UNROLLJ+jp] = nbat->nbfp[0*2+1];
329 SimdReal c6_S0 = load<SimdReal>(pvdw_c6 +0*UNROLLJ);
330 SimdReal c6_S1 = load<SimdReal>(pvdw_c6 +1*UNROLLJ);
331 SimdReal c6_S2 = load<SimdReal>(pvdw_c6 +2*UNROLLJ);
332 SimdReal c6_S3 = load<SimdReal>(pvdw_c6 +3*UNROLLJ);
334 SimdReal c12_S0 = load<SimdReal>(pvdw_c12+0*UNROLLJ);
335 SimdReal c12_S1 = load<SimdReal>(pvdw_c12+1*UNROLLJ);
336 SimdReal c12_S2 = load<SimdReal>(pvdw_c12+2*UNROLLJ);
337 SimdReal c12_S3 = load<SimdReal>(pvdw_c12+3*UNROLLJ);
338 #endif /* FIX_LJ_C */
341 egps_ishift = nbat->neg_2log;
342 egps_imask = (1<<egps_ishift) - 1;
343 egps_jshift = 2*nbat->neg_2log;
344 egps_jmask = (1<<egps_jshift) - 1;
345 egps_jstride = (UNROLLJ>>1)*UNROLLJ;
346 /* Major division is over i-particle energy groups, determine the stride */
347 Vstride_i = nbat->nenergrp*(1<<nbat->neg_2log)*egps_jstride;
353 for (n = 0; n < nbl->nci; n++)
357 ish = (nbln->shift & NBNXN_CI_SHIFT);
359 cjind0 = nbln->cj_ind_start;
360 cjind1 = nbln->cj_ind_end;
362 ci_sh = (ish == CENTRAL ? ci : -1);
364 shX_S = SimdReal(shiftvec[ish3]);
365 shY_S = SimdReal(shiftvec[ish3+1]);
366 shZ_S = SimdReal(shiftvec[ish3+2]);
371 #if defined LJ_COMB_LB || defined LJ_COMB_GEOM || defined LJ_EWALD_GEOM
375 int sci = (ci>>1)*STRIDE;
376 int scix = sci*DIM + (ci & 1)*(STRIDE>>1);
377 #if defined LJ_COMB_LB || defined LJ_COMB_GEOM || defined LJ_EWALD_GEOM
378 int sci2 = sci*2 + (ci & 1)*(STRIDE>>1);
380 sci += (ci & 1)*(STRIDE>>1);
383 /* We have 5 LJ/C combinations, but use only three inner loops,
384 * as the other combinations are unlikely and/or not much faster:
385 * inner half-LJ + C for half-LJ + C / no-LJ + C
386 * inner LJ + C for full-LJ + C
387 * inner LJ for full-LJ + no-C / half-LJ + no-C
389 do_LJ = (nbln->shift & NBNXN_CI_DO_LJ(0));
390 do_coul = (nbln->shift & NBNXN_CI_DO_COUL(0));
391 half_LJ = ((nbln->shift & NBNXN_CI_HALF_LJ(0)) || !do_LJ) && do_coul;
394 egps_i = nbat->energrp[ci];
398 for (ia = 0; ia < UNROLLI; ia++)
400 egp_ia = (egps_i >> (ia*egps_ishift)) & egps_imask;
401 vvdwtp[ia] = Vvdw + egp_ia*Vstride_i;
402 vctp[ia] = Vc + egp_ia*Vstride_i;
409 gmx_bool do_self = TRUE;
411 gmx_bool do_self = do_coul;
414 if (do_self && l_cj[nbln->cj_ind_start].cj == ci_sh)
417 if (do_self && l_cj[nbln->cj_ind_start].cj == (ci_sh>>1))
426 Vc_sub_self = 0.5*ic->c_rf;
430 Vc_sub_self = 0.5*tab_coul_F[2];
432 Vc_sub_self = 0.5*tab_coul_V[0];
435 #ifdef CALC_COUL_EWALD
437 Vc_sub_self = 0.5*ic->ewaldcoeff_q*M_2_SQRTPI;
440 for (ia = 0; ia < UNROLLI; ia++)
446 vctp[ia][((egps_i>>(ia*egps_ishift)) & egps_imask)*egps_jstride]
450 -= facel*qi*qi*Vc_sub_self;
458 for (ia = 0; ia < UNROLLI; ia++)
462 c6_i = nbat->nbfp[nbat->type[sci+ia]*(nbat->ntype + 1)*2]/6;
464 vvdwtp[ia][((egps_i>>(ia*egps_ishift)) & egps_imask)*egps_jstride]
468 += 0.5*c6_i*lj_ewaldcoeff6_6;
471 #endif /* LJ_EWALD */
475 /* Load i atom data */
476 int sciy = scix + STRIDE;
477 int sciz = sciy + STRIDE;
478 ix_S0 = loadU1DualHsimd(x+scix);
479 ix_S2 = loadU1DualHsimd(x+scix+2);
480 iy_S0 = loadU1DualHsimd(x+sciy);
481 iy_S2 = loadU1DualHsimd(x+sciy+2);
482 iz_S0 = loadU1DualHsimd(x+sciz);
483 iz_S2 = loadU1DualHsimd(x+sciz+2);
484 ix_S0 = ix_S0 + shX_S;
485 ix_S2 = ix_S2 + shX_S;
486 iy_S0 = iy_S0 + shY_S;
487 iy_S2 = iy_S2 + shY_S;
488 iz_S0 = iz_S0 + shZ_S;
489 iz_S2 = iz_S2 + shZ_S;
495 facel_S = SimdReal(facel);
497 iq_S0 = loadU1DualHsimd(q+sci);
498 iq_S2 = loadU1DualHsimd(q+sci+2);
499 iq_S0 = facel_S * iq_S0;
500 iq_S2 = facel_S * iq_S2;
504 hsig_i_S0 = loadU1DualHsimd(ljc+sci2);
505 hsig_i_S2 = loadU1DualHsimd(ljc+sci2+2);
506 seps_i_S0 = loadU1DualHsimd(ljc+sci2+STRIDE);
507 seps_i_S2 = loadU1DualHsimd(ljc+sci2+STRIDE+2);
510 SimdReal c6s_S0, c12s_S0;
511 SimdReal c6s_S2, c12s_S2;
513 c6s_S0 = loadU1DualHsimd(ljc+sci2);
517 c6s_S2 = loadU1DualHsimd(ljc+sci2+2);
519 c12s_S0 = loadU1DualHsimd(ljc+sci2+STRIDE);
522 c12s_S2 = loadU1DualHsimd(ljc+sci2+STRIDE+2);
524 #elif !defined LJ_COMB_LB && !defined FIX_LJ_C
525 const real *nbfp0 = nbfp_ptr + type[sci ]*nbat->ntype*c_simdBestPairAlignment;
526 const real *nbfp1 = nbfp_ptr + type[sci+1]*nbat->ntype*c_simdBestPairAlignment;
527 const real *nbfp2 = NULL, *nbfp3 = NULL;
530 nbfp2 = nbfp_ptr + type[sci+2]*nbat->ntype*c_simdBestPairAlignment;
531 nbfp3 = nbfp_ptr + type[sci+3]*nbat->ntype*c_simdBestPairAlignment;
536 /* We need the geometrically combined C6 for the PME grid correction */
537 SimdReal c6s_S0, c6s_S2;
538 c6s_S0 = loadU1DualHsimd(ljc+sci2);
541 c6s_S2 = loadU1DualHsimd(ljc+sci2+2);
545 /* Zero the potential energy for this list */
547 SimdReal Vvdwtot_S = setZero();
548 SimdReal vctot_S = setZero();
551 /* Clear i atom forces */
561 /* Currently all kernels use (at least half) LJ */
565 /* Coulomb: all i-atoms, LJ: first half i-atoms */
569 while (cjind < cjind1 && nbl->cj[cjind].excl != NBNXN_INTERACTION_MASK_ALL)
571 #include "gromacs/mdlib/nbnxn_kernels/simd_2xnn/nbnxn_kernel_simd_2xnn_inner.h"
575 for (; (cjind < cjind1); cjind++)
577 #include "gromacs/mdlib/nbnxn_kernels/simd_2xnn/nbnxn_kernel_simd_2xnn_inner.h"
584 /* Coulomb: all i-atoms, LJ: all i-atoms */
587 while (cjind < cjind1 && nbl->cj[cjind].excl != NBNXN_INTERACTION_MASK_ALL)
589 #include "gromacs/mdlib/nbnxn_kernels/simd_2xnn/nbnxn_kernel_simd_2xnn_inner.h"
593 for (; (cjind < cjind1); cjind++)
595 #include "gromacs/mdlib/nbnxn_kernels/simd_2xnn/nbnxn_kernel_simd_2xnn_inner.h"
601 /* Coulomb: none, LJ: all i-atoms */
603 while (cjind < cjind1 && nbl->cj[cjind].excl != NBNXN_INTERACTION_MASK_ALL)
605 #include "gromacs/mdlib/nbnxn_kernels/simd_2xnn/nbnxn_kernel_simd_2xnn_inner.h"
609 for (; (cjind < cjind1); cjind++)
611 #include "gromacs/mdlib/nbnxn_kernels/simd_2xnn/nbnxn_kernel_simd_2xnn_inner.h"
615 ninner += cjind1 - cjind0;
617 /* Add accumulated i-forces to the force array */
618 real fShiftX = reduceIncr4ReturnSumHsimd(f+scix, fix_S0, fix_S2);
619 real fShiftY = reduceIncr4ReturnSumHsimd(f+sciy, fiy_S0, fiy_S2);
620 real fShiftZ = reduceIncr4ReturnSumHsimd(f+sciz, fiz_S0, fiz_S2);
622 #ifdef CALC_SHIFTFORCES
623 fshift[ish3+0] += fShiftX;
624 fshift[ish3+1] += fShiftY;
625 fshift[ish3+2] += fShiftZ;
631 *Vc += reduce(vctot_S);
633 *Vvdw += reduce(Vvdwtot_S);
636 /* Outer loop uses 6 flops/iteration */
640 printf("atom pairs %d\n", npair);