2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2014,2015,2016,2017,2018 by the GROMACS development team.
5 * Copyright (c) 2019,2020, 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.
38 * \brief This file defines functions for "pair" interactions
39 * (i.e. listed non-bonded interactions, e.g. 1-4 interactions)
41 * \author Mark Abraham <mark.j.abraham@gmail.com>
43 * \ingroup module_listed_forces
51 #include "gromacs/listed_forces/bonded.h"
52 #include "gromacs/math/functions.h"
53 #include "gromacs/math/vec.h"
54 #include "gromacs/mdtypes/forcerec.h"
55 #include "gromacs/mdtypes/group.h"
56 #include "gromacs/mdtypes/interaction_const.h"
57 #include "gromacs/mdtypes/md_enums.h"
58 #include "gromacs/mdtypes/mdatom.h"
59 #include "gromacs/mdtypes/nblist.h"
60 #include "gromacs/mdtypes/simulation_workload.h"
61 #include "gromacs/pbcutil/ishift.h"
62 #include "gromacs/pbcutil/mshift.h"
63 #include "gromacs/pbcutil/pbc.h"
64 #include "gromacs/pbcutil/pbc_simd.h"
65 #include "gromacs/simd/simd.h"
66 #include "gromacs/simd/simd_math.h"
67 #include "gromacs/simd/vector_operations.h"
68 #include "gromacs/tables/forcetable.h"
69 #include "gromacs/utility/basedefinitions.h"
70 #include "gromacs/utility/fatalerror.h"
71 #include "gromacs/utility/gmxassert.h"
73 #include "listed_internal.h"
75 using namespace gmx; // TODO: Remove when this file is moved into gmx namespace
77 /*! \brief Issue a warning if a listed interaction is beyond a table limit */
78 static void warning_rlimit(const rvec* x, int ai, int aj, int* global_atom_index, real r, real rlimit)
81 "Listed nonbonded interaction between particles %d and %d\n"
82 "at distance %.3f which is larger than the table limit %.3f nm.\n\n"
83 "This is likely either a 1,4 interaction, or a listed interaction inside\n"
84 "a smaller molecule you are decoupling during a free energy calculation.\n"
85 "Since interactions at distances beyond the table cannot be computed,\n"
86 "they are skipped until they are inside the table limit again. You will\n"
87 "only see this message once, even if it occurs for several interactions.\n\n"
88 "IMPORTANT: This should not happen in a stable simulation, so there is\n"
89 "probably something wrong with your system. Only change the table-extension\n"
90 "distance in the mdp file if you are really sure that is the reason.\n",
91 glatnr(global_atom_index, ai), glatnr(global_atom_index, aj), r, rlimit);
96 "%8f %8f %8f\n%8f %8f %8f\n1-4 (%d,%d) interaction not within cut-off! r=%g. "
98 x[ai][XX], x[ai][YY], x[ai][ZZ], x[aj][XX], x[aj][YY], x[aj][ZZ],
99 glatnr(global_atom_index, ai), glatnr(global_atom_index, aj), r);
103 /*! \brief Compute the energy and force for a single pair interaction */
104 static real evaluate_single(real r2,
114 real rinv, r, rtab, eps, eps2, Y, F, Geps, Heps2, Fp, VVe, FFe, VVd, FFd, VVr, FFr, fscal;
117 /* Do the tabulated interactions - first table lookup */
118 rinv = gmx::invsqrt(r2);
121 ntab = static_cast<int>(rtab);
124 ntab = static_cast<int>(tableStride * ntab);
128 Geps = eps * vftab[ntab + 2];
129 Heps2 = eps2 * vftab[ntab + 3];
130 Fp = F + Geps + Heps2;
132 FFe = Fp + Geps + 2.0 * Heps2;
136 Geps = eps * vftab[ntab + 6];
137 Heps2 = eps2 * vftab[ntab + 7];
138 Fp = F + Geps + Heps2;
140 FFd = Fp + Geps + 2.0 * Heps2;
144 Geps = eps * vftab[ntab + 10];
145 Heps2 = eps2 * vftab[ntab + 11];
146 Fp = F + Geps + Heps2;
148 FFr = Fp + Geps + 2.0 * Heps2;
151 *vvdw = c6 * VVd + c12 * VVr;
153 fscal = -(qq * FFe + c6 * FFd + c12 * FFr) * tabscale * rinv;
158 /*! \brief Compute the energy and force for a single pair interaction under FEP */
159 static real free_energy_evaluate_single(real r2,
175 const real lfac_coul[2],
176 const real lfac_vdw[2],
177 const real dlfac_coul[2],
178 const real dlfac_vdw[2],
187 real rp, rpm2, rtab, eps, eps2, Y, F, Geps, Heps2, Fp, VV, FF, fscal;
188 real qq[2], c6[2], c12[2], sigma6[2], sigma2[2], sigma_pow[2];
189 real alpha_coul_eff, alpha_vdw_eff, dvdl_coul, dvdl_vdw;
190 real rpinv, r_coul, r_vdw, velecsum, vvdwsum;
191 real fscal_vdw[2], fscal_elec[2];
192 real velec[2], vvdw[2];
194 const real half = 0.5;
195 const real minusOne = -1.0;
196 const real one = 1.0;
197 const real two = 2.0;
198 const real six = 6.0;
207 if (sc_r_power == six)
209 rpm2 = r2 * r2; /* r4 */
210 rp = rpm2 * r2; /* r6 */
214 rp = std::pow(r2, half * sc_r_power); /* not currently supported as input, but can handle it */
218 /* Loop over state A(0) and B(1) */
219 for (i = 0; i < 2; i++)
221 if ((c6[i] > 0) && (c12[i] > 0))
223 /* The c6 & c12 coefficients now contain the constants 6.0 and 12.0, respectively.
224 * Correct for this by multiplying with (1/12.0)/(1/6.0)=6.0/12.0=0.5.
226 sigma6[i] = half * c12[i] / c6[i];
227 sigma2[i] = std::cbrt(half * c12[i] / c6[i]);
228 /* should be able to get rid of this ^^^ internal pow call eventually. Will require agreement on
229 what data to store externally. Can't be fixed without larger scale changes, so not 5.0 */
230 if (sigma6[i] < sigma6_min) /* for disappearing coul and vdw with soft core at the same time */
232 sigma6[i] = sigma6_min;
233 sigma2[i] = sigma2_min;
238 sigma6[i] = sigma6_def;
239 sigma2[i] = sigma2_def;
241 if (sc_r_power == six)
243 sigma_pow[i] = sigma6[i];
246 { /* not really supported as input, but in here for testing the general case*/
247 sigma_pow[i] = std::pow(sigma2[i], sc_r_power / 2);
251 /* only use softcore if one of the states has a zero endstate - softcore is for avoiding infinities!*/
252 if ((c12[0] > 0) && (c12[1] > 0))
259 alpha_vdw_eff = alpha_vdw;
260 alpha_coul_eff = alpha_coul;
263 /* Loop over A and B states again */
264 for (i = 0; i < 2; i++)
271 /* Only spend time on A or B state if it is non-zero */
272 if ((qq[i] != 0) || (c6[i] != 0) || (c12[i] != 0))
275 rpinv = one / (alpha_coul_eff * lfac_coul[i] * sigma_pow[i] + rp);
276 r_coul = std::pow(rpinv, minusOne / sc_r_power);
278 /* Electrostatics table lookup data */
279 rtab = r_coul * tabscale;
280 ntab = static_cast<int>(rtab);
283 ntab = static_cast<int>(tableStride * ntab);
287 Geps = eps * vftab[ntab + 2];
288 Heps2 = eps2 * vftab[ntab + 3];
289 Fp = F + Geps + Heps2;
291 FF = Fp + Geps + two * Heps2;
292 velec[i] = qq[i] * VV;
293 fscal_elec[i] = -qq[i] * FF * r_coul * rpinv * tabscale;
296 rpinv = one / (alpha_vdw_eff * lfac_vdw[i] * sigma_pow[i] + rp);
297 r_vdw = std::pow(rpinv, minusOne / sc_r_power);
298 /* Vdw table lookup data */
299 rtab = r_vdw * tabscale;
300 ntab = static_cast<int>(rtab);
307 Geps = eps * vftab[ntab + 6];
308 Heps2 = eps2 * vftab[ntab + 7];
309 Fp = F + Geps + Heps2;
311 FF = Fp + Geps + two * Heps2;
312 vvdw[i] = c6[i] * VV;
313 fscal_vdw[i] = -c6[i] * FF;
318 Geps = eps * vftab[ntab + 10];
319 Heps2 = eps2 * vftab[ntab + 11];
320 Fp = F + Geps + Heps2;
322 FF = Fp + Geps + two * Heps2;
323 vvdw[i] += c12[i] * VV;
324 fscal_vdw[i] -= c12[i] * FF;
325 fscal_vdw[i] *= r_vdw * rpinv * tabscale;
328 /* Now we have velec[i], vvdw[i], and fscal[i] for both states */
329 /* Assemble A and B states */
335 for (i = 0; i < 2; i++)
337 velecsum += LFC[i] * velec[i];
338 vvdwsum += LFV[i] * vvdw[i];
340 fscal += (LFC[i] * fscal_elec[i] + LFV[i] * fscal_vdw[i]) * rpm2;
342 dvdl_coul += velec[i] * DLF[i]
343 + LFC[i] * alpha_coul_eff * dlfac_coul[i] * fscal_elec[i] * sigma_pow[i];
344 dvdl_vdw += vvdw[i] * DLF[i]
345 + LFV[i] * alpha_vdw_eff * dlfac_vdw[i] * fscal_vdw[i] * sigma_pow[i];
348 dvdl[efptCOUL] += dvdl_coul;
349 dvdl[efptVDW] += dvdl_vdw;
351 *velectot = velecsum;
357 /*! \brief Calculate pair interactions, supports all types and conditions. */
358 template<BondedKernelFlavor flavor>
359 static real do_pairs_general(int ftype,
361 const t_iatom iatoms[],
362 const t_iparams iparams[],
366 const struct t_pbc* pbc,
367 const struct t_graph* g,
371 const t_forcerec* fr,
372 gmx_grppairener_t* grppener,
373 int* global_atom_index)
378 int i, itype, ai, aj, gid;
381 real fscal, velec, vvdw;
382 real* energygrp_elec;
384 static gmx_bool warned_rlimit = FALSE;
385 /* Free energy stuff */
386 gmx_bool bFreeEnergy;
387 real LFC[2], LFV[2], DLF[2], lfac_coul[2], lfac_vdw[2], dlfac_coul[2], dlfac_vdw[2];
388 real qqB, c6B, c12B, sigma2_def, sigma2_min;
394 energygrp_elec = grppener->ener[egCOUL14].data();
395 energygrp_vdw = grppener->ener[egLJ14].data();
398 energygrp_elec = grppener->ener[egCOULSR].data();
399 energygrp_vdw = grppener->ener[egLJSR].data();
402 energygrp_elec = nullptr; /* Keep compiler happy */
403 energygrp_vdw = nullptr; /* Keep compiler happy */
404 gmx_fatal(FARGS, "Unknown function type %d in do_nonbonded14", ftype);
407 if (fr->efep != efepNO)
409 /* Lambda factor for state A=1-lambda and B=lambda */
410 LFC[0] = 1.0 - lambda[efptCOUL];
411 LFV[0] = 1.0 - lambda[efptVDW];
412 LFC[1] = lambda[efptCOUL];
413 LFV[1] = lambda[efptVDW];
415 /*derivative of the lambda factor for state A and B */
420 sigma2_def = std::cbrt(fr->sc_sigma6_def);
421 sigma2_min = std::cbrt(fr->sc_sigma6_min);
423 for (i = 0; i < 2; i++)
425 lfac_coul[i] = (fr->sc_power == 2 ? (1 - LFC[i]) * (1 - LFC[i]) : (1 - LFC[i]));
427 DLF[i] * fr->sc_power / fr->sc_r_power * (fr->sc_power == 2 ? (1 - LFC[i]) : 1);
428 lfac_vdw[i] = (fr->sc_power == 2 ? (1 - LFV[i]) * (1 - LFV[i]) : (1 - LFV[i]));
430 DLF[i] * fr->sc_power / fr->sc_r_power * (fr->sc_power == 2 ? (1 - LFV[i]) : 1);
435 sigma2_min = sigma2_def = 0;
438 /* TODO This code depends on the logic in tables.c that constructs
439 the table layout, which should be made explicit in future
443 "Pair-interaction code that uses GROMACS interaction tables supports exactly 3 tables");
445 fr->pairsTable->interaction == GMX_TABLE_INTERACTION_ELEC_VDWREP_VDWDISP,
446 "Pair interaction kernels need a table with Coulomb, repulsion and dispersion entries");
448 const real epsfac = fr->ic->epsfac;
451 for (i = 0; (i < nbonds);)
456 gid = GID(md->cENER[ai], md->cENER[aj], md->nenergrp);
462 bFreeEnergy = (fr->efep != efepNO
463 && (((md->nPerturbed != 0) && (md->bPerturbed[ai] || md->bPerturbed[aj]))
464 || iparams[itype].lj14.c6A != iparams[itype].lj14.c6B
465 || iparams[itype].lj14.c12A != iparams[itype].lj14.c12B));
466 qq = md->chargeA[ai] * md->chargeA[aj] * epsfac * fr->fudgeQQ;
467 c6 = iparams[itype].lj14.c6A;
468 c12 = iparams[itype].lj14.c12A;
471 qq = iparams[itype].ljc14.qi * iparams[itype].ljc14.qj * epsfac
472 * iparams[itype].ljc14.fqq;
473 c6 = iparams[itype].ljc14.c6;
474 c12 = iparams[itype].ljc14.c12;
477 qq = iparams[itype].ljcnb.qi * iparams[itype].ljcnb.qj * epsfac;
478 c6 = iparams[itype].ljcnb.c6;
479 c12 = iparams[itype].ljcnb.c12;
482 /* Cannot happen since we called gmx_fatal() above in this case */
483 qq = c6 = c12 = 0; /* Keep compiler happy */
487 /* To save flops in the optimized kernels, c6/c12 have 6.0/12.0 derivative prefactors
488 * included in the general nfbp array now. This means the tables are scaled down by the
489 * same factor, so when we use the original c6/c12 parameters from iparams[] they must
495 /* Do we need to apply full periodic boundary conditions? */
498 fshift_index = pbc_dx_aiuc(pbc, x[ai], x[aj], dx);
502 fshift_index = CENTRAL;
503 rvec_sub(x[ai], x[aj], dx);
507 if (r2 >= fr->pairsTable->r * fr->pairsTable->r)
509 /* This check isn't race free. But it doesn't matter because if a race occurs the only
510 * disadvantage is that the warning is printed twice */
513 warning_rlimit(x, ai, aj, global_atom_index, sqrt(r2), fr->pairsTable->r);
514 warned_rlimit = TRUE;
521 /* Currently free energy is only supported for F_LJ14, so no need to check for that if we got here */
522 qqB = md->chargeB[ai] * md->chargeB[aj] * epsfac * fr->fudgeQQ;
523 c6B = iparams[itype].lj14.c6B * 6.0;
524 c12B = iparams[itype].lj14.c12B * 12.0;
526 fscal = free_energy_evaluate_single(
527 r2, fr->sc_r_power, fr->sc_alphacoul, fr->sc_alphavdw, fr->pairsTable->scale,
528 fr->pairsTable->data, fr->pairsTable->stride, qq, c6, c12, qqB, c6B, c12B, LFC,
529 LFV, DLF, lfac_coul, lfac_vdw, dlfac_coul, dlfac_vdw, fr->sc_sigma6_def,
530 fr->sc_sigma6_min, sigma2_def, sigma2_min, &velec, &vvdw, dvdl);
534 /* Evaluate tabulated interaction without free energy */
535 fscal = evaluate_single(r2, fr->pairsTable->scale, fr->pairsTable->data,
536 fr->pairsTable->stride, qq, c6, c12, &velec, &vvdw);
539 energygrp_elec[gid] += velec;
540 energygrp_vdw[gid] += vvdw;
541 svmul(fscal, dx, dx);
547 if (computeVirial(flavor))
551 /* Correct the shift forces using the graph */
552 ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, aj), dt);
553 fshift_index = IVEC2IS(dt);
555 if (fshift_index != CENTRAL)
557 rvec_inc(fshift[fshift_index], dx);
558 rvec_dec(fshift[CENTRAL], dx);
565 /*! \brief Calculate pairs, only for plain-LJ + plain Coulomb normal type.
567 * This function is templated for real/SimdReal and for optimization.
569 template<typename T, int pack_size, typename pbc_type>
570 static void do_pairs_simple(int nbonds,
571 const t_iatom iatoms[],
572 const t_iparams iparams[],
577 const real scale_factor)
579 const int nfa1 = 1 + 2;
585 #if GMX_SIMD_HAVE_REAL
586 alignas(GMX_SIMD_ALIGNMENT) std::int32_t ai[pack_size];
587 alignas(GMX_SIMD_ALIGNMENT) std::int32_t aj[pack_size];
588 alignas(GMX_SIMD_ALIGNMENT) real coeff[3 * pack_size];
590 std::int32_t ai[pack_size];
591 std::int32_t aj[pack_size];
592 real coeff[3 * pack_size];
595 /* nbonds is #pairs*nfa1, here we step pack_size pairs */
596 for (int i = 0; i < nbonds; i += pack_size * nfa1)
598 /* Collect atoms for pack_size pairs.
599 * iu indexes into iatoms, we should not let iu go beyond nbonds.
602 for (int s = 0; s < pack_size; s++)
604 int itype = iatoms[iu];
605 ai[s] = iatoms[iu + 1];
606 aj[s] = iatoms[iu + 2];
608 if (i + s * nfa1 < nbonds)
610 coeff[0 * pack_size + s] = iparams[itype].lj14.c6A;
611 coeff[1 * pack_size + s] = iparams[itype].lj14.c12A;
612 coeff[2 * pack_size + s] = md->chargeA[ai[s]] * md->chargeA[aj[s]];
614 /* Avoid indexing the iatoms array out of bounds.
615 * We pad the coordinate indices with the last atom pair.
617 if (iu + nfa1 < nbonds)
624 /* Pad the coefficient arrays with zeros to get zero forces */
625 coeff[0 * pack_size + s] = 0;
626 coeff[1 * pack_size + s] = 0;
627 coeff[2 * pack_size + s] = 0;
631 /* Load the coordinates */
633 gatherLoadUTranspose<3>(reinterpret_cast<const real*>(x), ai, &xi[XX], &xi[YY], &xi[ZZ]);
634 gatherLoadUTranspose<3>(reinterpret_cast<const real*>(x), aj, &xj[XX], &xj[YY], &xj[ZZ]);
636 T c6 = load<T>(coeff + 0 * pack_size);
637 T c12 = load<T>(coeff + 1 * pack_size);
638 T qq = load<T>(coeff + 2 * pack_size);
640 /* We could save these operations by storing 6*C6,12*C12 */
645 pbc_dx_aiuc(pbc, xi, xj, dr);
647 T rsq = dr[XX] * dr[XX] + dr[YY] * dr[YY] + dr[ZZ] * dr[ZZ];
648 T rinv = gmx::invsqrt(rsq);
649 T rinv2 = rinv * rinv;
650 T rinv6 = rinv2 * rinv2 * rinv2;
652 /* Calculate the Coulomb force * r */
653 T cfr = ef * qq * rinv;
655 /* Calculate the LJ force * r and add it to the Coulomb part */
656 T fr = gmx::fma(fms(c12, rinv6, c6), rinv6, cfr);
658 T finvr = fr * rinv2;
659 T fx = finvr * dr[XX];
660 T fy = finvr * dr[YY];
661 T fz = finvr * dr[ZZ];
663 /* Add the pair forces to the force array.
664 * Note that here we might add multiple force components for some atoms
665 * due to the SIMD padding. But the extra force components are zero.
667 transposeScatterIncrU<4>(reinterpret_cast<real*>(f), ai, fx, fy, fz);
668 transposeScatterDecrU<4>(reinterpret_cast<real*>(f), aj, fx, fy, fz);
672 /*! \brief Calculate all listed pair interactions */
673 void do_pairs(int ftype,
675 const t_iatom iatoms[],
676 const t_iparams iparams[],
680 const struct t_pbc* pbc,
681 const struct t_graph* g,
685 const t_forcerec* fr,
686 const bool havePerturbedInteractions,
687 const gmx::StepWorkload& stepWork,
688 gmx_grppairener_t* grppener,
689 int* global_atom_index)
691 if (ftype == F_LJ14 && fr->ic->vdwtype != evdwUSER && !EEL_USER(fr->ic->eeltype)
692 && !havePerturbedInteractions && (!stepWork.computeVirial && !stepWork.computeEnergy))
694 /* We use a fast code-path for plain LJ 1-4 without FEP.
696 * TODO: Add support for energies (straightforward) and virial
697 * in the SIMD template. For the virial it's inconvenient to store
698 * the force sums for the shifts and we should directly calculate
699 * and sum the virial for the shifts. But we should do this
700 * at once for the angles and dihedrals as well.
702 #if GMX_SIMD_HAVE_REAL
703 if (fr->use_simd_kernels)
705 alignas(GMX_SIMD_ALIGNMENT) real pbc_simd[9 * GMX_SIMD_REAL_WIDTH];
706 set_pbc_simd(pbc, pbc_simd);
708 do_pairs_simple<SimdReal, GMX_SIMD_REAL_WIDTH, const real*>(
709 nbonds, iatoms, iparams, x, f, pbc_simd, md, fr->ic->epsfac * fr->fudgeQQ);
714 /* This construct is needed because pbc_dx_aiuc doesn't accept pbc=NULL */
716 const t_pbc* pbc_nonnull;
724 set_pbc(&pbc_no, PbcType::No, nullptr);
725 pbc_nonnull = &pbc_no;
728 do_pairs_simple<real, 1, const t_pbc*>(nbonds, iatoms, iparams, x, f, pbc_nonnull, md,
729 fr->ic->epsfac * fr->fudgeQQ);
732 else if (stepWork.computeVirial)
734 do_pairs_general<BondedKernelFlavor::ForcesAndVirialAndEnergy>(
735 ftype, nbonds, iatoms, iparams, x, f, fshift, pbc, g, lambda, dvdl, md, fr,
736 grppener, global_atom_index);
740 do_pairs_general<BondedKernelFlavor::ForcesAndEnergy>(ftype, nbonds, iatoms, iparams, x, f,
741 fshift, pbc, g, lambda, dvdl, md, fr,
742 grppener, global_atom_index);