2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2012,2013,2014,2015,2016 by the GROMACS development team.
5 * Copyright (c) 2017,2018,2019,2020,2021, by the GROMACS development team, led by
6 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
7 * and including many others, as listed in the AUTHORS file in the
8 * top-level source directory and at http://www.gromacs.org.
10 * GROMACS is free software; you can redistribute it and/or
11 * modify it under the terms of the GNU Lesser General Public License
12 * as published by the Free Software Foundation; either version 2.1
13 * of the License, or (at your option) any later version.
15 * GROMACS is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 * Lesser General Public License for more details.
20 * You should have received a copy of the GNU Lesser General Public
21 * License along with GROMACS; if not, see
22 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
23 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
25 * If you want to redistribute modifications to GROMACS, please
26 * consider that scientific software is very special. Version
27 * control is crucial - bugs must be traceable. We will be happy to
28 * consider code for inclusion in the official distribution, but
29 * derived work must not be called official GROMACS. Details are found
30 * in the README & COPYING files - if they are missing, get the
31 * official version at http://www.gromacs.org.
33 * To help us fund GROMACS development, we humbly ask that you cite
34 * the research papers on the package. Check out http://www.gromacs.org.
41 /* Unpack pointers for output */
42 real* f = out->f.data();
43 real* fshift = out->fshift.data();
46 real* Vvdw = out->VSvdw.data();
47 real* Vc = out->VSc.data();
49 real* Vvdw = out->Vvdw.data();
50 real* Vc = out->Vc.data();
57 SimdReal ix_S0, iy_S0, iz_S0;
58 SimdReal ix_S2, iy_S2, iz_S2;
59 SimdReal fix_S0, fiy_S0, fiz_S0;
60 SimdReal fix_S2, fiy_S2, fiz_S2;
62 SimdReal diagonal_jmi_S;
63 #if UNROLLI == UNROLLJ
64 SimdBool diagonal_mask_S0, diagonal_mask_S2;
66 SimdBool diagonal_mask0_S0, diagonal_mask0_S2;
67 SimdBool diagonal_mask1_S0, diagonal_mask1_S2;
70 SimdBitMask filter_S0, filter_S2;
75 SimdReal iq_S0 = setZero();
76 SimdReal iq_S2 = setZero();
81 SimdReal hrc_3_S, moh_rc_S;
86 /* Coulomb table variables */
94 #ifdef CALC_COUL_EWALD
95 SimdReal beta2_S, beta_S;
98 #if defined CALC_ENERGIES && (defined CALC_COUL_EWALD || defined CALC_COUL_TAB)
102 #if defined LJ_CUT && defined CALC_ENERGIES
103 SimdReal p6_cpot_S, p12_cpot_S;
107 SimdReal swV3_S, swV4_S, swV5_S;
108 SimdReal swF2_S, swF3_S, swF4_S;
110 #ifdef LJ_FORCE_SWITCH
112 SimdReal p6_fc2_S, p6_fc3_S;
113 SimdReal p12_fc2_S, p12_fc3_S;
114 # ifdef CALC_ENERGIES
115 SimdReal p6_vc3_S, p6_vc4_S;
116 SimdReal p12_vc3_S, p12_vc4_S;
117 SimdReal p6_6cpot_S, p12_12cpot_S;
121 SimdReal half_S, lje_c2_S, lje_c6_6_S;
125 SimdReal hsig_i_S0, seps_i_S0;
126 SimdReal hsig_i_S2, seps_i_S2;
129 alignas(GMX_SIMD_ALIGNMENT) real pvdw_c6[2 * UNROLLI * UNROLLJ];
130 real* pvdw_c12 = pvdw_c6 + UNROLLI * UNROLLJ;
132 #endif /* LJ_COMB_LB */
136 #ifdef VDW_CUTOFF_CHECK
144 const nbnxn_atomdata_t::Params& nbatParams = nbat->params();
146 #if defined LJ_COMB_GEOM || defined LJ_COMB_LB || defined LJ_EWALD_GEOM
147 const real* gmx_restrict ljc = nbatParams.lj_comb.data();
149 #if !(defined LJ_COMB_GEOM || defined LJ_COMB_LB || defined FIX_LJ_C)
150 /* No combination rule used */
151 const real* gmx_restrict nbfp_ptr = nbatParams.nbfp_aligned.data();
152 const int* gmx_restrict type = nbatParams.type.data();
155 /* Load j-i for the first i */
156 diagonal_jmi_S = load<SimdReal>(nbat->simdMasks.diagonal_2xnn_j_minus_i.data());
157 /* Generate all the diagonal masks as comparison results */
158 #if UNROLLI == UNROLLJ
159 diagonal_mask_S0 = (zero_S < diagonal_jmi_S);
160 diagonal_jmi_S = diagonal_jmi_S - one_S;
161 diagonal_jmi_S = diagonal_jmi_S - one_S;
162 diagonal_mask_S2 = (zero_S < diagonal_jmi_S);
164 # if 2 * UNROLLI == UNROLLJ
165 diagonal_mask0_S0 = (zero_S < diagonal_jmi_S);
166 diagonal_jmi_S = diagonal_jmi_S - one_S;
167 diagonal_jmi_S = diagonal_jmi_S - one_S;
168 diagonal_mask0_S2 = (zero_S < diagonal_jmi_S);
169 diagonal_jmi_S = diagonal_jmi_S - one_S;
170 diagonal_jmi_S = diagonal_jmi_S - one_S;
171 diagonal_mask1_S0 = (zero_S < diagonal_jmi_S);
172 diagonal_jmi_S = diagonal_jmi_S - one_S;
173 diagonal_jmi_S = diagonal_jmi_S - one_S;
174 diagonal_mask1_S2 = (zero_S < diagonal_jmi_S);
178 /* Load masks for topology exclusion masking. filter_stride is
179 static const, so the conditional will be optimized away. */
180 #if GMX_DOUBLE && !GMX_SIMD_HAVE_INT32_LOGICAL
181 const std::uint64_t* gmx_restrict exclusion_filter = nbat->simdMasks.exclusion_filter64.data();
183 const std::uint32_t* gmx_restrict exclusion_filter = nbat->simdMasks.exclusion_filter.data();
186 /* Here we cast the exclusion filters from unsigned * to int * or real *.
187 * Since we only check bits, the actual value they represent does not
188 * matter, as long as both filter and mask data are treated the same way.
190 #if GMX_SIMD_HAVE_INT32_LOGICAL
191 filter_S0 = load<SimdBitMask>(reinterpret_cast<const int*>(exclusion_filter + 0 * UNROLLJ));
192 filter_S2 = load<SimdBitMask>(reinterpret_cast<const int*>(exclusion_filter + 2 * UNROLLJ));
194 filter_S0 = load<SimdBitMask>(reinterpret_cast<const real*>(exclusion_filter + 0 * UNROLLJ));
195 filter_S2 = load<SimdBitMask>(reinterpret_cast<const real*>(exclusion_filter + 2 * UNROLLJ));
199 /* Reaction-field constants */
200 mrc_3_S = SimdReal(-2 * ic->reactionFieldCoefficient);
201 # ifdef CALC_ENERGIES
202 hrc_3_S = SimdReal(ic->reactionFieldCoefficient);
203 moh_rc_S = SimdReal(-ic->reactionFieldShift);
209 invtsp_S = SimdReal(ic->coulombEwaldTables->scale);
210 # ifdef CALC_ENERGIES
211 mhalfsp_S = SimdReal(-0.5_real / ic->coulombEwaldTables->scale);
215 const real* tab_coul_F = ic->coulombEwaldTables->tableFDV0.data();
217 const real* tab_coul_F = ic->coulombEwaldTables->tableF.data();
218 # ifdef CALC_ENERGIES
219 const real* tab_coul_V = ic->coulombEwaldTables->tableV.data();
222 #endif /* CALC_COUL_TAB */
224 #ifdef CALC_COUL_EWALD
225 beta2_S = SimdReal(ic->ewaldcoeff_q * ic->ewaldcoeff_q);
226 beta_S = SimdReal(ic->ewaldcoeff_q);
229 #if (defined CALC_COUL_TAB || defined CALC_COUL_EWALD) && defined CALC_ENERGIES
230 sh_ewald_S = SimdReal(ic->sh_ewald);
233 /* LJ function constants */
234 #if defined CALC_ENERGIES || defined LJ_POT_SWITCH
235 SimdReal sixth_S = SimdReal(1.0 / 6.0);
236 SimdReal twelveth_S = SimdReal(1.0 / 12.0);
239 #if defined LJ_CUT && defined CALC_ENERGIES
240 /* We shift the potential by cpot, which can be zero */
241 p6_cpot_S = SimdReal(ic->dispersion_shift.cpot);
242 p12_cpot_S = SimdReal(ic->repulsion_shift.cpot);
245 rswitch_S = SimdReal(ic->rvdw_switch);
246 swV3_S = SimdReal(ic->vdw_switch.c3);
247 swV4_S = SimdReal(ic->vdw_switch.c4);
248 swV5_S = SimdReal(ic->vdw_switch.c5);
249 swF2_S = SimdReal(3 * ic->vdw_switch.c3);
250 swF3_S = SimdReal(4 * ic->vdw_switch.c4);
251 swF4_S = SimdReal(5 * ic->vdw_switch.c5);
253 #ifdef LJ_FORCE_SWITCH
254 rswitch_S = SimdReal(ic->rvdw_switch);
255 p6_fc2_S = SimdReal(ic->dispersion_shift.c2);
256 p6_fc3_S = SimdReal(ic->dispersion_shift.c3);
257 p12_fc2_S = SimdReal(ic->repulsion_shift.c2);
258 p12_fc3_S = SimdReal(ic->repulsion_shift.c3);
259 # ifdef CALC_ENERGIES
261 SimdReal mthird_S = SimdReal(-1.0 / 3.0);
262 SimdReal mfourth_S = SimdReal(-1.0 / 4.0);
264 p6_vc3_S = mthird_S * p6_fc2_S;
265 p6_vc4_S = mfourth_S * p6_fc3_S;
266 p6_6cpot_S = SimdReal(ic->dispersion_shift.cpot / 6);
267 p12_vc3_S = mthird_S * p12_fc2_S;
268 p12_vc4_S = mfourth_S * p12_fc3_S;
269 p12_12cpot_S = SimdReal(ic->repulsion_shift.cpot / 12);
274 half_S = SimdReal(0.5);
275 const real lj_ewaldcoeff2 = ic->ewaldcoeff_lj * ic->ewaldcoeff_lj;
276 const real lj_ewaldcoeff6_6 = lj_ewaldcoeff2 * lj_ewaldcoeff2 * lj_ewaldcoeff2 / 6;
277 lje_c2_S = SimdReal(lj_ewaldcoeff2);
278 lje_c6_6_S = SimdReal(lj_ewaldcoeff6_6);
279 # ifdef CALC_ENERGIES
280 /* Determine the grid potential at the cut-off */
281 SimdReal lje_vc_S = SimdReal(ic->sh_lj_ewald);
285 /* The kernel either supports rcoulomb = rvdw or rcoulomb >= rvdw */
286 rc2_S = SimdReal(ic->rcoulomb * ic->rcoulomb);
287 #ifdef VDW_CUTOFF_CHECK
288 rcvdw2_S = SimdReal(ic->rvdw * ic->rvdw);
291 minRsq_S = SimdReal(c_nbnxnMinDistanceSquared);
293 const real* gmx_restrict q = nbatParams.q.data();
294 const real facel = ic->epsfac;
295 const real* gmx_restrict shiftvec = shift_vec[0];
296 const real* gmx_restrict x = nbat->x().data();
300 for (jp = 0; jp < UNROLLJ; jp++)
302 pvdw_c6[0 * UNROLLJ + jp] = nbat->nbfp[0 * 2];
303 pvdw_c6[1 * UNROLLJ + jp] = nbat->nbfp[0 * 2];
304 pvdw_c6[2 * UNROLLJ + jp] = nbat->nbfp[0 * 2];
305 pvdw_c6[3 * UNROLLJ + jp] = nbat->nbfp[0 * 2];
307 pvdw_c12[0 * UNROLLJ + jp] = nbat->nbfp[0 * 2 + 1];
308 pvdw_c12[1 * UNROLLJ + jp] = nbat->nbfp[0 * 2 + 1];
309 pvdw_c12[2 * UNROLLJ + jp] = nbat->nbfp[0 * 2 + 1];
310 pvdw_c12[3 * UNROLLJ + jp] = nbat->nbfp[0 * 2 + 1];
312 SimdReal c6_S0 = load<SimdReal>(pvdw_c6 + 0 * UNROLLJ);
313 SimdReal c6_S1 = load<SimdReal>(pvdw_c6 + 1 * UNROLLJ);
314 SimdReal c6_S2 = load<SimdReal>(pvdw_c6 + 2 * UNROLLJ);
315 SimdReal c6_S3 = load<SimdReal>(pvdw_c6 + 3 * UNROLLJ);
317 SimdReal c12_S0 = load<SimdReal>(pvdw_c12 + 0 * UNROLLJ);
318 SimdReal c12_S1 = load<SimdReal>(pvdw_c12 + 1 * UNROLLJ);
319 SimdReal c12_S2 = load<SimdReal>(pvdw_c12 + 2 * UNROLLJ);
320 SimdReal c12_S3 = load<SimdReal>(pvdw_c12 + 3 * UNROLLJ);
321 #endif /* FIX_LJ_C */
324 const int egps_ishift = nbatParams.neg_2log;
325 const int egps_imask = (1 << egps_ishift) - 1;
326 const int egps_jshift = 2 * nbatParams.neg_2log;
327 const int egps_jmask = (1 << egps_jshift) - 1;
328 const int egps_jstride = (UNROLLJ >> 1) * UNROLLJ;
329 /* Major division is over i-particle energy groups, determine the stride */
330 const int Vstride_i = nbatParams.nenergrp * (1 << nbatParams.neg_2log) * egps_jstride;
333 const nbnxn_cj_t* l_cj = nbl->cj.data();
336 for (const nbnxn_ci_t& ciEntry : nbl->ci)
338 const int ish = (ciEntry.shift & NBNXN_CI_SHIFT);
339 const int ish3 = ish * 3;
340 const int cjind0 = ciEntry.cj_ind_start;
341 const int cjind1 = ciEntry.cj_ind_end;
342 const int ci = ciEntry.ci;
343 const int ci_sh = (ish == gmx::c_centralShiftIndex ? ci : -1);
345 shX_S = SimdReal(shiftvec[ish3]);
346 shY_S = SimdReal(shiftvec[ish3 + 1]);
347 shZ_S = SimdReal(shiftvec[ish3 + 2]);
350 int sci = ci * STRIDE;
351 int scix = sci * DIM;
352 # if defined LJ_COMB_LB || defined LJ_COMB_GEOM || defined LJ_EWALD_GEOM
356 int sci = (ci >> 1) * STRIDE;
357 int scix = sci * DIM + (ci & 1) * (STRIDE >> 1);
358 # if defined LJ_COMB_LB || defined LJ_COMB_GEOM || defined LJ_EWALD_GEOM
359 int sci2 = sci * 2 + (ci & 1) * (STRIDE >> 1);
361 sci += (ci & 1) * (STRIDE >> 1);
364 /* We have 5 LJ/C combinations, but use only three inner loops,
365 * as the other combinations are unlikely and/or not much faster:
366 * inner half-LJ + C for half-LJ + C / no-LJ + C
367 * inner LJ + C for full-LJ + C
368 * inner LJ for full-LJ + no-C / half-LJ + no-C
370 const bool do_LJ = ((ciEntry.shift & NBNXN_CI_DO_LJ(0)) != 0);
371 const bool do_coul = ((ciEntry.shift & NBNXN_CI_DO_COUL(0)) != 0);
372 const bool half_LJ = (((ciEntry.shift & NBNXN_CI_HALF_LJ(0)) != 0) || !do_LJ) && do_coul;
375 const int egps_i = nbatParams.energrp[ci];
376 real* vvdwtp[UNROLLI];
379 for (int ia = 0; ia < UNROLLI; ia++)
381 int egp_ia = (egps_i >> (ia * egps_ishift)) & egps_imask;
382 vvdwtp[ia] = Vvdw + egp_ia * Vstride_i;
383 vctp[ia] = Vc + egp_ia * Vstride_i;
389 # ifdef LJ_EWALD_GEOM
390 gmx_bool do_self = TRUE;
392 gmx_bool do_self = do_coul;
395 if (do_self && l_cj[ciEntry.cj_ind_start].cj == ci_sh)
398 if (do_self && l_cj[ciEntry.cj_ind_start].cj == (ci_sh >> 1))
404 const real Vc_sub_self = 0.5 * ic->reactionFieldShift;
406 # ifdef CALC_COUL_TAB
408 const real Vc_sub_self = 0.5 * tab_coul_F[2];
410 const real Vc_sub_self = 0.5 * tab_coul_V[0];
413 # ifdef CALC_COUL_EWALD
415 const real Vc_sub_self = 0.5 * ic->ewaldcoeff_q * M_2_SQRTPI;
418 for (int ia = 0; ia < UNROLLI; ia++)
420 const real qi = q[sci + ia];
421 # ifdef ENERGY_GROUPS
422 vctp[ia][((egps_i >> (ia * egps_ishift)) & egps_imask) * egps_jstride]
426 -= facel * qi * qi * Vc_sub_self;
430 # ifdef LJ_EWALD_GEOM
432 for (int ia = 0; ia < UNROLLI; ia++)
435 nbatParams.nbfp[nbatParams.type[sci + ia] * (nbatParams.numTypes + 1) * 2]
437 # ifdef ENERGY_GROUPS
438 vvdwtp[ia][((egps_i >> (ia * egps_ishift)) & egps_imask) * egps_jstride]
442 += 0.5 * c6_i * lj_ewaldcoeff6_6;
445 # endif /* LJ_EWALD */
449 /* Load i atom data */
450 int sciy = scix + STRIDE;
451 int sciz = sciy + STRIDE;
452 ix_S0 = loadU1DualHsimd(x + scix);
453 ix_S2 = loadU1DualHsimd(x + scix + 2);
454 iy_S0 = loadU1DualHsimd(x + sciy);
455 iy_S2 = loadU1DualHsimd(x + sciy + 2);
456 iz_S0 = loadU1DualHsimd(x + sciz);
457 iz_S2 = loadU1DualHsimd(x + sciz + 2);
458 ix_S0 = ix_S0 + shX_S;
459 ix_S2 = ix_S2 + shX_S;
460 iy_S0 = iy_S0 + shY_S;
461 iy_S2 = iy_S2 + shY_S;
462 iz_S0 = iz_S0 + shZ_S;
463 iz_S2 = iz_S2 + shZ_S;
469 facel_S = SimdReal(facel);
471 iq_S0 = loadU1DualHsimd(q + sci);
472 iq_S2 = loadU1DualHsimd(q + sci + 2);
473 iq_S0 = facel_S * iq_S0;
474 iq_S2 = facel_S * iq_S2;
478 hsig_i_S0 = loadU1DualHsimd(ljc + sci2);
479 hsig_i_S2 = loadU1DualHsimd(ljc + sci2 + 2);
480 seps_i_S0 = loadU1DualHsimd(ljc + sci2 + STRIDE);
481 seps_i_S2 = loadU1DualHsimd(ljc + sci2 + STRIDE + 2);
484 SimdReal c6s_S0, c12s_S0;
485 SimdReal c6s_S2, c12s_S2;
487 c6s_S0 = loadU1DualHsimd(ljc + sci2);
491 c6s_S2 = loadU1DualHsimd(ljc + sci2 + 2);
493 c12s_S0 = loadU1DualHsimd(ljc + sci2 + STRIDE);
496 c12s_S2 = loadU1DualHsimd(ljc + sci2 + STRIDE + 2);
498 # elif !defined LJ_COMB_LB && !defined FIX_LJ_C
499 const int numTypes = nbatParams.numTypes;
500 const real* nbfp0 = nbfp_ptr + type[sci] * numTypes * c_simdBestPairAlignment;
501 const real* nbfp1 = nbfp_ptr + type[sci + 1] * numTypes * c_simdBestPairAlignment;
502 const real *nbfp2 = nullptr, *nbfp3 = nullptr;
505 nbfp2 = nbfp_ptr + type[sci + 2] * numTypes * c_simdBestPairAlignment;
506 nbfp3 = nbfp_ptr + type[sci + 3] * numTypes * c_simdBestPairAlignment;
511 /* We need the geometrically combined C6 for the PME grid correction */
512 SimdReal c6s_S0, c6s_S2;
513 c6s_S0 = loadU1DualHsimd(ljc + sci2);
516 c6s_S2 = loadU1DualHsimd(ljc + sci2 + 2);
520 /* Zero the potential energy for this list */
522 SimdReal Vvdwtot_S = setZero();
523 SimdReal vctot_S = setZero();
526 /* Clear i atom forces */
536 /* Currently all kernels use (at least half) LJ */
540 /* Coulomb: all i-atoms, LJ: first half i-atoms */
544 while (cjind < cjind1 && nbl->cj[cjind].excl != NBNXN_INTERACTION_MASK_ALL)
546 #include "kernel_inner.h"
550 for (; (cjind < cjind1); cjind++)
552 #include "kernel_inner.h"
559 /* Coulomb: all i-atoms, LJ: all i-atoms */
562 while (cjind < cjind1 && nbl->cj[cjind].excl != NBNXN_INTERACTION_MASK_ALL)
564 #include "kernel_inner.h"
568 for (; (cjind < cjind1); cjind++)
570 #include "kernel_inner.h"
576 /* Coulomb: none, LJ: all i-atoms */
578 while (cjind < cjind1 && nbl->cj[cjind].excl != NBNXN_INTERACTION_MASK_ALL)
580 #include "kernel_inner.h"
584 for (; (cjind < cjind1); cjind++)
586 #include "kernel_inner.h"
590 ninner += cjind1 - cjind0;
592 /* Add accumulated i-forces to the force array */
593 real fShiftX = reduceIncr4ReturnSumHsimd(f + scix, fix_S0, fix_S2);
594 real fShiftY = reduceIncr4ReturnSumHsimd(f + sciy, fiy_S0, fiy_S2);
595 real fShiftZ = reduceIncr4ReturnSumHsimd(f + sciz, fiz_S0, fiz_S2);
597 #ifdef CALC_SHIFTFORCES
598 fshift[ish3 + 0] += fShiftX;
599 fshift[ish3 + 1] += fShiftY;
600 fshift[ish3 + 2] += fShiftZ;
606 *Vc += reduce(vctot_S);
608 *Vvdw += reduce(Vvdwtot_S);
611 /* Outer loop uses 6 flops/iteration */
615 printf("atom pairs %d\n", npair);