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,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.
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/pbc.h"
63 #include "gromacs/pbcutil/pbc_simd.h"
64 #include "gromacs/simd/simd.h"
65 #include "gromacs/simd/simd_math.h"
66 #include "gromacs/simd/vector_operations.h"
67 #include "gromacs/tables/forcetable.h"
68 #include "gromacs/utility/basedefinitions.h"
69 #include "gromacs/utility/fatalerror.h"
70 #include "gromacs/utility/gmxassert.h"
72 #include "listed_internal.h"
74 using namespace gmx; // TODO: Remove when this file is moved into gmx namespace
76 /*! \brief Issue a warning if a listed interaction is beyond a table limit */
77 static void warning_rlimit(const rvec* x, int ai, int aj, int* global_atom_index, real r, real rlimit)
80 "Listed nonbonded interaction between particles %d and %d\n"
81 "at distance %.3f which is larger than the table limit %.3f nm.\n\n"
82 "This is likely either a 1,4 interaction, or a listed interaction inside\n"
83 "a smaller molecule you are decoupling during a free energy calculation.\n"
84 "Since interactions at distances beyond the table cannot be computed,\n"
85 "they are skipped until they are inside the table limit again. You will\n"
86 "only see this message once, even if it occurs for several interactions.\n\n"
87 "IMPORTANT: This should not happen in a stable simulation, so there is\n"
88 "probably something wrong with your system. Only change the table-extension\n"
89 "distance in the mdp file if you are really sure that is the reason.\n",
90 glatnr(global_atom_index, ai),
91 glatnr(global_atom_index, aj),
98 "%8f %8f %8f\n%8f %8f %8f\n1-4 (%d,%d) interaction not within cut-off! r=%g. "
106 glatnr(global_atom_index, ai),
107 glatnr(global_atom_index, aj),
112 /*! \brief Compute the energy and force for a single pair interaction */
113 static real evaluate_single(real r2,
123 real rinv, r, rtab, eps, eps2, Y, F, Geps, Heps2, Fp, VVe, FFe, VVd, FFd, VVr, FFr, fscal;
126 /* Do the tabulated interactions - first table lookup */
127 rinv = gmx::invsqrt(r2);
130 ntab = static_cast<int>(rtab);
133 ntab = static_cast<int>(tableStride * ntab);
137 Geps = eps * vftab[ntab + 2];
138 Heps2 = eps2 * vftab[ntab + 3];
139 Fp = F + Geps + Heps2;
141 FFe = Fp + Geps + 2.0 * Heps2;
145 Geps = eps * vftab[ntab + 6];
146 Heps2 = eps2 * vftab[ntab + 7];
147 Fp = F + Geps + Heps2;
149 FFd = Fp + Geps + 2.0 * Heps2;
153 Geps = eps * vftab[ntab + 10];
154 Heps2 = eps2 * vftab[ntab + 11];
155 Fp = F + Geps + Heps2;
157 FFr = Fp + Geps + 2.0 * Heps2;
160 *vvdw = c6 * VVd + c12 * VVr;
162 fscal = -(qq * FFe + c6 * FFd + c12 * FFr) * tabscale * rinv;
167 static inline real sixthRoot(const real r)
169 return gmx::invsqrt(std::cbrt(r));
172 /*! \brief Compute the energy and force for a single pair interaction under FEP */
173 static real free_energy_evaluate_single(real r2,
174 const interaction_const_t::SoftCoreParameters& scParams,
187 const real lfac_coul[2],
188 const real lfac_vdw[2],
189 const real dlfac_coul[2],
190 const real dlfac_vdw[2],
195 real rtab, eps, eps2, Y, F, Geps, Heps2, Fp, VV, FF, fscal;
196 real qq[2], c6[2], c12[2], sigma6[2], sigma_pow[2];
197 real alpha_coul_eff, alpha_vdw_eff, dvdl_coul, dvdl_vdw;
198 real rpinv, r_coul, r_vdw, velecsum, vvdwsum;
199 real fscal_vdw[2], fscal_elec[2];
200 real velec[2], vvdw[2];
202 const real half = 0.5_real;
203 const real one = 1.0_real;
204 const real two = 2.0_real;
213 const real rpm2 = r2 * r2; /* r4 */
214 const real rp = rpm2 * r2; /* r6 */
216 /* Loop over state A(0) and B(1) */
217 for (i = 0; i < 2; i++)
219 if ((c6[i] > 0) && (c12[i] > 0))
221 /* The c6 & c12 coefficients now contain the constants 6.0 and 12.0, respectively.
222 * Correct for this by multiplying with (1/12.0)/(1/6.0)=6.0/12.0=0.5.
224 sigma6[i] = half * c12[i] / c6[i];
225 if (sigma6[i] < scParams.sigma6Minimum) /* for disappearing coul and vdw with soft core at the same time */
227 sigma6[i] = scParams.sigma6Minimum;
232 sigma6[i] = scParams.sigma6WithInvalidSigma;
234 sigma_pow[i] = sigma6[i];
237 /* only use softcore if one of the states has a zero endstate - softcore is for avoiding infinities!*/
238 if ((c12[0] > 0) && (c12[1] > 0))
245 alpha_vdw_eff = scParams.alphaVdw;
246 alpha_coul_eff = scParams.alphaCoulomb;
249 /* Loop over A and B states again */
250 for (i = 0; i < 2; i++)
257 /* Only spend time on A or B state if it is non-zero */
258 if ((qq[i] != 0) || (c6[i] != 0) || (c12[i] != 0))
261 rpinv = one / (alpha_coul_eff * lfac_coul[i] * sigma_pow[i] + rp);
262 r_coul = sixthRoot(rpinv);
264 /* Electrostatics table lookup data */
265 rtab = r_coul * tabscale;
266 ntab = static_cast<int>(rtab);
269 ntab = static_cast<int>(tableStride * ntab);
273 Geps = eps * vftab[ntab + 2];
274 Heps2 = eps2 * vftab[ntab + 3];
275 Fp = F + Geps + Heps2;
277 FF = Fp + Geps + two * Heps2;
278 velec[i] = qq[i] * VV;
279 fscal_elec[i] = -qq[i] * FF * r_coul * rpinv * tabscale;
282 rpinv = one / (alpha_vdw_eff * lfac_vdw[i] * sigma_pow[i] + rp);
283 r_vdw = sixthRoot(rpinv);
284 /* Vdw table lookup data */
285 rtab = r_vdw * tabscale;
286 ntab = static_cast<int>(rtab);
293 Geps = eps * vftab[ntab + 6];
294 Heps2 = eps2 * vftab[ntab + 7];
295 Fp = F + Geps + Heps2;
297 FF = Fp + Geps + two * Heps2;
298 vvdw[i] = c6[i] * VV;
299 fscal_vdw[i] = -c6[i] * FF;
304 Geps = eps * vftab[ntab + 10];
305 Heps2 = eps2 * vftab[ntab + 11];
306 Fp = F + Geps + Heps2;
308 FF = Fp + Geps + two * Heps2;
309 vvdw[i] += c12[i] * VV;
310 fscal_vdw[i] -= c12[i] * FF;
311 fscal_vdw[i] *= r_vdw * rpinv * tabscale;
314 /* Now we have velec[i], vvdw[i], and fscal[i] for both states */
315 /* Assemble A and B states */
321 for (i = 0; i < 2; i++)
323 velecsum += LFC[i] * velec[i];
324 vvdwsum += LFV[i] * vvdw[i];
326 fscal += (LFC[i] * fscal_elec[i] + LFV[i] * fscal_vdw[i]) * rpm2;
328 dvdl_coul += velec[i] * DLF[i]
329 + LFC[i] * alpha_coul_eff * dlfac_coul[i] * fscal_elec[i] * sigma_pow[i];
330 dvdl_vdw += vvdw[i] * DLF[i]
331 + LFV[i] * alpha_vdw_eff * dlfac_vdw[i] * fscal_vdw[i] * sigma_pow[i];
334 dvdl[static_cast<int>(FreeEnergyPerturbationCouplingType::Coul)] += dvdl_coul;
335 dvdl[static_cast<int>(FreeEnergyPerturbationCouplingType::Vdw)] += dvdl_vdw;
337 *velectot = velecsum;
343 /*! \brief Calculate pair interactions, supports all types and conditions. */
344 template<BondedKernelFlavor flavor>
345 static real do_pairs_general(int ftype,
347 const t_iatom iatoms[],
348 const t_iparams iparams[],
352 const struct t_pbc* pbc,
356 const t_forcerec* fr,
357 gmx_grppairener_t* grppener,
358 int* global_atom_index)
362 int i, itype, ai, aj, gid;
365 real fscal, velec, vvdw;
366 real* energygrp_elec;
368 static gmx_bool warned_rlimit = FALSE;
369 /* Free energy stuff */
370 gmx_bool bFreeEnergy;
371 real LFC[2], LFV[2], DLF[2], lfac_coul[2], lfac_vdw[2], dlfac_coul[2], dlfac_vdw[2];
373 const real oneSixth = 1.0_real / 6.0_real;
379 energygrp_elec = grppener->energyGroupPairTerms[NonBondedEnergyTerms::Coulomb14].data();
380 energygrp_vdw = grppener->energyGroupPairTerms[NonBondedEnergyTerms::LJ14].data();
383 energygrp_elec = grppener->energyGroupPairTerms[NonBondedEnergyTerms::CoulombSR].data();
384 energygrp_vdw = grppener->energyGroupPairTerms[NonBondedEnergyTerms::LJSR].data();
387 energygrp_elec = nullptr; /* Keep compiler happy */
388 energygrp_vdw = nullptr; /* Keep compiler happy */
389 gmx_fatal(FARGS, "Unknown function type %d in do_nonbonded14", ftype);
392 if (fr->efep != FreeEnergyPerturbationType::No)
394 /* Lambda factor for state A=1-lambda and B=lambda */
395 LFC[0] = 1.0 - lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Coul)];
396 LFV[0] = 1.0 - lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Vdw)];
397 LFC[1] = lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Coul)];
398 LFV[1] = lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Vdw)];
400 /*derivative of the lambda factor for state A and B */
404 GMX_ASSERT(fr->ic->softCoreParameters, "We need soft-core parameters");
405 const auto& scParams = *fr->ic->softCoreParameters;
407 for (i = 0; i < 2; i++)
409 lfac_coul[i] = (scParams.lambdaPower == 2 ? (1 - LFC[i]) * (1 - LFC[i]) : (1 - LFC[i]));
410 dlfac_coul[i] = DLF[i] * scParams.lambdaPower * oneSixth
411 * (scParams.lambdaPower == 2 ? (1 - LFC[i]) : 1);
412 lfac_vdw[i] = (scParams.lambdaPower == 2 ? (1 - LFV[i]) * (1 - LFV[i]) : (1 - LFV[i]));
413 dlfac_vdw[i] = DLF[i] * scParams.lambdaPower * oneSixth
414 * (scParams.lambdaPower == 2 ? (1 - LFV[i]) : 1);
418 /* TODO This code depends on the logic in tables.c that constructs
419 the table layout, which should be made explicit in future
423 "Pair-interaction code that uses GROMACS interaction tables supports exactly 3 tables");
425 fr->pairsTable->interaction == GMX_TABLE_INTERACTION_ELEC_VDWREP_VDWDISP,
426 "Pair interaction kernels need a table with Coulomb, repulsion and dispersion entries");
428 const real epsfac = fr->ic->epsfac;
431 for (i = 0; (i < nbonds);)
436 gid = GID(md->cENER[ai], md->cENER[aj], md->nenergrp);
442 bFreeEnergy = (fr->efep != FreeEnergyPerturbationType::No
443 && (((md->nPerturbed != 0) && (md->bPerturbed[ai] || md->bPerturbed[aj]))
444 || iparams[itype].lj14.c6A != iparams[itype].lj14.c6B
445 || iparams[itype].lj14.c12A != iparams[itype].lj14.c12B));
446 qq = md->chargeA[ai] * md->chargeA[aj] * epsfac * fr->fudgeQQ;
447 c6 = iparams[itype].lj14.c6A;
448 c12 = iparams[itype].lj14.c12A;
451 qq = iparams[itype].ljc14.qi * iparams[itype].ljc14.qj * epsfac
452 * iparams[itype].ljc14.fqq;
453 c6 = iparams[itype].ljc14.c6;
454 c12 = iparams[itype].ljc14.c12;
457 qq = iparams[itype].ljcnb.qi * iparams[itype].ljcnb.qj * epsfac;
458 c6 = iparams[itype].ljcnb.c6;
459 c12 = iparams[itype].ljcnb.c12;
462 /* Cannot happen since we called gmx_fatal() above in this case */
463 qq = c6 = c12 = 0; /* Keep compiler happy */
467 /* To save flops in the optimized kernels, c6/c12 have 6.0/12.0 derivative prefactors
468 * included in the general nfbp array now. This means the tables are scaled down by the
469 * same factor, so when we use the original c6/c12 parameters from iparams[] they must
475 /* Do we need to apply full periodic boundary conditions? */
478 fshift_index = pbc_dx_aiuc(pbc, x[ai], x[aj], dx);
482 fshift_index = CENTRAL;
483 rvec_sub(x[ai], x[aj], dx);
487 if (r2 >= fr->pairsTable->r * fr->pairsTable->r)
489 /* This check isn't race free. But it doesn't matter because if a race occurs the only
490 * disadvantage is that the warning is printed twice */
493 warning_rlimit(x, ai, aj, global_atom_index, sqrt(r2), fr->pairsTable->r);
494 warned_rlimit = TRUE;
501 /* Currently free energy is only supported for F_LJ14, so no need to check for that if we got here */
502 qqB = md->chargeB[ai] * md->chargeB[aj] * epsfac * fr->fudgeQQ;
503 c6B = iparams[itype].lj14.c6B * 6.0;
504 c12B = iparams[itype].lj14.c12B * 12.0;
506 fscal = free_energy_evaluate_single(r2,
507 *fr->ic->softCoreParameters,
508 fr->pairsTable->scale,
509 fr->pairsTable->data.data(),
510 fr->pairsTable->stride,
530 /* Evaluate tabulated interaction without free energy */
531 fscal = evaluate_single(r2,
532 fr->pairsTable->scale,
533 fr->pairsTable->data.data(),
534 fr->pairsTable->stride,
542 energygrp_elec[gid] += velec;
543 energygrp_vdw[gid] += vvdw;
544 svmul(fscal, dx, dx);
550 if (computeVirial(flavor))
552 if (fshift_index != CENTRAL)
554 rvec_inc(fshift[fshift_index], dx);
555 rvec_dec(fshift[CENTRAL], dx);
562 /*! \brief Calculate pairs, only for plain-LJ + plain Coulomb normal type.
564 * This function is templated for real/SimdReal and for optimization.
566 template<typename T, int pack_size, typename pbc_type>
567 static void do_pairs_simple(int nbonds,
568 const t_iatom iatoms[],
569 const t_iparams iparams[],
574 const real scale_factor)
576 const int nfa1 = 1 + 2;
582 #if GMX_SIMD_HAVE_REAL
583 alignas(GMX_SIMD_ALIGNMENT) std::int32_t ai[pack_size];
584 alignas(GMX_SIMD_ALIGNMENT) std::int32_t aj[pack_size];
585 alignas(GMX_SIMD_ALIGNMENT) real coeff[3 * pack_size];
587 std::int32_t ai[pack_size];
588 std::int32_t aj[pack_size];
589 real coeff[3 * pack_size];
592 /* nbonds is #pairs*nfa1, here we step pack_size pairs */
593 for (int i = 0; i < nbonds; i += pack_size * nfa1)
595 /* Collect atoms for pack_size pairs.
596 * iu indexes into iatoms, we should not let iu go beyond nbonds.
599 for (int s = 0; s < pack_size; s++)
601 int itype = iatoms[iu];
602 ai[s] = iatoms[iu + 1];
603 aj[s] = iatoms[iu + 2];
605 if (i + s * nfa1 < nbonds)
607 coeff[0 * pack_size + s] = iparams[itype].lj14.c6A;
608 coeff[1 * pack_size + s] = iparams[itype].lj14.c12A;
609 coeff[2 * pack_size + s] = md->chargeA[ai[s]] * md->chargeA[aj[s]];
611 /* Avoid indexing the iatoms array out of bounds.
612 * We pad the coordinate indices with the last atom pair.
614 if (iu + nfa1 < nbonds)
621 /* Pad the coefficient arrays with zeros to get zero forces */
622 coeff[0 * pack_size + s] = 0;
623 coeff[1 * pack_size + s] = 0;
624 coeff[2 * pack_size + s] = 0;
628 /* Load the coordinates */
630 gatherLoadUTranspose<3>(reinterpret_cast<const real*>(x), ai, &xi[XX], &xi[YY], &xi[ZZ]);
631 gatherLoadUTranspose<3>(reinterpret_cast<const real*>(x), aj, &xj[XX], &xj[YY], &xj[ZZ]);
633 T c6 = load<T>(coeff + 0 * pack_size);
634 T c12 = load<T>(coeff + 1 * pack_size);
635 T qq = load<T>(coeff + 2 * pack_size);
637 /* We could save these operations by storing 6*C6,12*C12 */
642 pbc_dx_aiuc(pbc, xi, xj, dr);
644 T rsq = dr[XX] * dr[XX] + dr[YY] * dr[YY] + dr[ZZ] * dr[ZZ];
645 T rinv = gmx::invsqrt(rsq);
646 T rinv2 = rinv * rinv;
647 T rinv6 = rinv2 * rinv2 * rinv2;
649 /* Calculate the Coulomb force * r */
650 T cfr = ef * qq * rinv;
652 /* Calculate the LJ force * r and add it to the Coulomb part */
653 T fr = gmx::fma(fms(c12, rinv6, c6), rinv6, cfr);
655 T finvr = fr * rinv2;
656 T fx = finvr * dr[XX];
657 T fy = finvr * dr[YY];
658 T fz = finvr * dr[ZZ];
660 /* Add the pair forces to the force array.
661 * Note that here we might add multiple force components for some atoms
662 * due to the SIMD padding. But the extra force components are zero.
664 transposeScatterIncrU<4>(reinterpret_cast<real*>(f), ai, fx, fy, fz);
665 transposeScatterDecrU<4>(reinterpret_cast<real*>(f), aj, fx, fy, fz);
669 /*! \brief Calculate all listed pair interactions */
670 void do_pairs(int ftype,
672 const t_iatom iatoms[],
673 const t_iparams iparams[],
677 const struct t_pbc* pbc,
681 const t_forcerec* fr,
682 const bool havePerturbedInteractions,
683 const gmx::StepWorkload& stepWork,
684 gmx_grppairener_t* grppener,
685 int* global_atom_index)
687 if (ftype == F_LJ14 && fr->ic->vdwtype != VanDerWaalsType::User && !EEL_USER(fr->ic->eeltype)
688 && !havePerturbedInteractions && (!stepWork.computeVirial && !stepWork.computeEnergy))
690 /* We use a fast code-path for plain LJ 1-4 without FEP.
692 * TODO: Add support for energies (straightforward) and virial
693 * in the SIMD template. For the virial it's inconvenient to store
694 * the force sums for the shifts and we should directly calculate
695 * and sum the virial for the shifts. But we should do this
696 * at once for the angles and dihedrals as well.
698 #if GMX_SIMD_HAVE_REAL
699 if (fr->use_simd_kernels)
701 alignas(GMX_SIMD_ALIGNMENT) real pbc_simd[9 * GMX_SIMD_REAL_WIDTH];
702 set_pbc_simd(pbc, pbc_simd);
704 do_pairs_simple<SimdReal, GMX_SIMD_REAL_WIDTH, const real*>(
705 nbonds, iatoms, iparams, x, f, pbc_simd, md, fr->ic->epsfac * fr->fudgeQQ);
710 /* This construct is needed because pbc_dx_aiuc doesn't accept pbc=NULL */
712 const t_pbc* pbc_nonnull;
720 set_pbc(&pbc_no, PbcType::No, nullptr);
721 pbc_nonnull = &pbc_no;
724 do_pairs_simple<real, 1, const t_pbc*>(
725 nbonds, iatoms, iparams, x, f, pbc_nonnull, md, fr->ic->epsfac * fr->fudgeQQ);
728 else if (stepWork.computeVirial)
730 do_pairs_general<BondedKernelFlavor::ForcesAndVirialAndEnergy>(
731 ftype, nbonds, iatoms, iparams, x, f, fshift, pbc, lambda, dvdl, md, fr, grppener, global_atom_index);
735 do_pairs_general<BondedKernelFlavor::ForcesAndEnergy>(
736 ftype, nbonds, iatoms, iparams, x, f, fshift, pbc, lambda, dvdl, md, fr, grppener, global_atom_index);