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/nblist.h"
59 #include "gromacs/mdtypes/simulation_workload.h"
60 #include "gromacs/pbcutil/ishift.h"
61 #include "gromacs/pbcutil/pbc.h"
62 #include "gromacs/pbcutil/pbc_simd.h"
63 #include "gromacs/simd/simd.h"
64 #include "gromacs/simd/simd_math.h"
65 #include "gromacs/simd/vector_operations.h"
66 #include "gromacs/tables/forcetable.h"
67 #include "gromacs/utility/basedefinitions.h"
68 #include "gromacs/utility/fatalerror.h"
69 #include "gromacs/utility/gmxassert.h"
71 #include "listed_internal.h"
73 using namespace gmx; // TODO: Remove when this file is moved into gmx namespace
75 /*! \brief Issue a warning if a listed interaction is beyond a table limit */
76 static void warning_rlimit(const rvec* x, int ai, int aj, int* global_atom_index, real r, real rlimit)
79 "Listed nonbonded interaction between particles %d and %d\n"
80 "at distance %.3f which is larger than the table limit %.3f nm.\n\n"
81 "This is likely either a 1,4 interaction, or a listed interaction inside\n"
82 "a smaller molecule you are decoupling during a free energy calculation.\n"
83 "Since interactions at distances beyond the table cannot be computed,\n"
84 "they are skipped until they are inside the table limit again. You will\n"
85 "only see this message once, even if it occurs for several interactions.\n\n"
86 "IMPORTANT: This should not happen in a stable simulation, so there is\n"
87 "probably something wrong with your system. Only change the table-extension\n"
88 "distance in the mdp file if you are really sure that is the reason.\n",
89 glatnr(global_atom_index, ai),
90 glatnr(global_atom_index, aj),
97 "%8f %8f %8f\n%8f %8f %8f\n1-4 (%d,%d) interaction not within cut-off! r=%g. "
105 glatnr(global_atom_index, ai),
106 glatnr(global_atom_index, aj),
111 /*! \brief Compute the energy and force for a single pair interaction */
112 static real evaluate_single(real r2,
122 real rinv, r, rtab, eps, eps2, Y, F, Geps, Heps2, Fp, VVe, FFe, VVd, FFd, VVr, FFr, fscal;
125 /* Do the tabulated interactions - first table lookup */
126 rinv = gmx::invsqrt(r2);
129 ntab = static_cast<int>(rtab);
132 ntab = static_cast<int>(tableStride * ntab);
136 Geps = eps * vftab[ntab + 2];
137 Heps2 = eps2 * vftab[ntab + 3];
138 Fp = F + Geps + Heps2;
140 FFe = Fp + Geps + 2.0 * Heps2;
144 Geps = eps * vftab[ntab + 6];
145 Heps2 = eps2 * vftab[ntab + 7];
146 Fp = F + Geps + Heps2;
148 FFd = Fp + Geps + 2.0 * Heps2;
152 Geps = eps * vftab[ntab + 10];
153 Heps2 = eps2 * vftab[ntab + 11];
154 Fp = F + Geps + Heps2;
156 FFr = Fp + Geps + 2.0 * Heps2;
159 *vvdw = c6 * VVd + c12 * VVr;
161 fscal = -(qq * FFe + c6 * FFd + c12 * FFr) * tabscale * rinv;
166 static inline real sixthRoot(const real r)
168 return gmx::invsqrt(std::cbrt(r));
171 /*! \brief Compute the energy and force for a single pair interaction under FEP */
172 static real free_energy_evaluate_single(real r2,
173 const interaction_const_t::SoftCoreParameters& scParams,
186 const real lfac_coul[2],
187 const real lfac_vdw[2],
188 const real dlfac_coul[2],
189 const real dlfac_vdw[2],
194 real rtab, eps, eps2, Y, F, Geps, Heps2, Fp, VV, FF, fscal;
195 real qq[2], c6[2], c12[2], sigma6[2], sigma_pow[2];
196 real alpha_coul_eff, alpha_vdw_eff, dvdl_coul, dvdl_vdw;
197 real rpinv, r_coul, r_vdw, velecsum, vvdwsum;
198 real fscal_vdw[2], fscal_elec[2];
199 real velec[2], vvdw[2];
201 const real half = 0.5_real;
202 const real one = 1.0_real;
203 const real two = 2.0_real;
212 const real rpm2 = r2 * r2; /* r4 */
213 const real rp = rpm2 * r2; /* r6 */
215 /* Loop over state A(0) and B(1) */
216 for (i = 0; i < 2; i++)
218 if ((c6[i] > 0) && (c12[i] > 0))
220 /* The c6 & c12 coefficients now contain the constants 6.0 and 12.0, respectively.
221 * Correct for this by multiplying with (1/12.0)/(1/6.0)=6.0/12.0=0.5.
223 sigma6[i] = half * c12[i] / c6[i];
224 if (sigma6[i] < scParams.sigma6Minimum) /* for disappearing coul and vdw with soft core at the same time */
226 sigma6[i] = scParams.sigma6Minimum;
231 sigma6[i] = scParams.sigma6WithInvalidSigma;
233 sigma_pow[i] = sigma6[i];
236 /* only use softcore if one of the states has a zero endstate - softcore is for avoiding infinities!*/
237 if ((c12[0] > 0) && (c12[1] > 0))
244 alpha_vdw_eff = scParams.alphaVdw;
245 alpha_coul_eff = scParams.alphaCoulomb;
248 /* Loop over A and B states again */
249 for (i = 0; i < 2; i++)
256 /* Only spend time on A or B state if it is non-zero */
257 if ((qq[i] != 0) || (c6[i] != 0) || (c12[i] != 0))
260 rpinv = one / (alpha_coul_eff * lfac_coul[i] * sigma_pow[i] + rp);
261 r_coul = sixthRoot(rpinv);
263 /* Electrostatics table lookup data */
264 rtab = r_coul * tabscale;
265 ntab = static_cast<int>(rtab);
268 ntab = static_cast<int>(tableStride * ntab);
272 Geps = eps * vftab[ntab + 2];
273 Heps2 = eps2 * vftab[ntab + 3];
274 Fp = F + Geps + Heps2;
276 FF = Fp + Geps + two * Heps2;
277 velec[i] = qq[i] * VV;
278 fscal_elec[i] = -qq[i] * FF * r_coul * rpinv * tabscale;
281 rpinv = one / (alpha_vdw_eff * lfac_vdw[i] * sigma_pow[i] + rp);
282 r_vdw = sixthRoot(rpinv);
283 /* Vdw table lookup data */
284 rtab = r_vdw * tabscale;
285 ntab = static_cast<int>(rtab);
292 Geps = eps * vftab[ntab + 6];
293 Heps2 = eps2 * vftab[ntab + 7];
294 Fp = F + Geps + Heps2;
296 FF = Fp + Geps + two * Heps2;
297 vvdw[i] = c6[i] * VV;
298 fscal_vdw[i] = -c6[i] * FF;
303 Geps = eps * vftab[ntab + 10];
304 Heps2 = eps2 * vftab[ntab + 11];
305 Fp = F + Geps + Heps2;
307 FF = Fp + Geps + two * Heps2;
308 vvdw[i] += c12[i] * VV;
309 fscal_vdw[i] -= c12[i] * FF;
310 fscal_vdw[i] *= r_vdw * rpinv * tabscale;
313 /* Now we have velec[i], vvdw[i], and fscal[i] for both states */
314 /* Assemble A and B states */
320 for (i = 0; i < 2; i++)
322 velecsum += LFC[i] * velec[i];
323 vvdwsum += LFV[i] * vvdw[i];
325 fscal += (LFC[i] * fscal_elec[i] + LFV[i] * fscal_vdw[i]) * rpm2;
327 dvdl_coul += velec[i] * DLF[i]
328 + LFC[i] * alpha_coul_eff * dlfac_coul[i] * fscal_elec[i] * sigma_pow[i];
329 dvdl_vdw += vvdw[i] * DLF[i]
330 + LFV[i] * alpha_vdw_eff * dlfac_vdw[i] * fscal_vdw[i] * sigma_pow[i];
333 dvdl[static_cast<int>(FreeEnergyPerturbationCouplingType::Coul)] += dvdl_coul;
334 dvdl[static_cast<int>(FreeEnergyPerturbationCouplingType::Vdw)] += dvdl_vdw;
336 *velectot = velecsum;
342 /*! \brief Calculate pair interactions, supports all types and conditions. */
343 template<BondedKernelFlavor flavor>
344 static real do_pairs_general(int ftype,
346 const t_iatom iatoms[],
347 const t_iparams iparams[],
351 const struct t_pbc* pbc,
354 gmx::ArrayRef<real> chargeA,
355 gmx::ArrayRef<real> chargeB,
356 gmx::ArrayRef<bool> atomIsPerturbed,
357 gmx::ArrayRef<unsigned short> cENER,
359 const t_forcerec* fr,
360 gmx_grppairener_t* grppener,
361 int* global_atom_index)
365 int i, itype, ai, aj, gid;
368 real fscal, velec, vvdw;
369 real* energygrp_elec;
371 static gmx_bool warned_rlimit = FALSE;
372 /* Free energy stuff */
373 gmx_bool bFreeEnergy;
374 real LFC[2], LFV[2], DLF[2], lfac_coul[2], lfac_vdw[2], dlfac_coul[2], dlfac_vdw[2];
376 const real oneSixth = 1.0_real / 6.0_real;
382 energygrp_elec = grppener->energyGroupPairTerms[NonBondedEnergyTerms::Coulomb14].data();
383 energygrp_vdw = grppener->energyGroupPairTerms[NonBondedEnergyTerms::LJ14].data();
386 energygrp_elec = grppener->energyGroupPairTerms[NonBondedEnergyTerms::CoulombSR].data();
387 energygrp_vdw = grppener->energyGroupPairTerms[NonBondedEnergyTerms::LJSR].data();
390 energygrp_elec = nullptr; /* Keep compiler happy */
391 energygrp_vdw = nullptr; /* Keep compiler happy */
392 gmx_fatal(FARGS, "Unknown function type %d in do_nonbonded14", ftype);
395 if (fr->efep != FreeEnergyPerturbationType::No)
397 /* Lambda factor for state A=1-lambda and B=lambda */
398 LFC[0] = 1.0 - lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Coul)];
399 LFV[0] = 1.0 - lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Vdw)];
400 LFC[1] = lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Coul)];
401 LFV[1] = lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Vdw)];
403 /*derivative of the lambda factor for state A and B */
407 GMX_ASSERT(fr->ic->softCoreParameters, "We need soft-core parameters");
408 const auto& scParams = *fr->ic->softCoreParameters;
410 for (i = 0; i < 2; i++)
412 lfac_coul[i] = (scParams.lambdaPower == 2 ? (1 - LFC[i]) * (1 - LFC[i]) : (1 - LFC[i]));
413 dlfac_coul[i] = DLF[i] * scParams.lambdaPower * oneSixth
414 * (scParams.lambdaPower == 2 ? (1 - LFC[i]) : 1);
415 lfac_vdw[i] = (scParams.lambdaPower == 2 ? (1 - LFV[i]) * (1 - LFV[i]) : (1 - LFV[i]));
416 dlfac_vdw[i] = DLF[i] * scParams.lambdaPower * oneSixth
417 * (scParams.lambdaPower == 2 ? (1 - LFV[i]) : 1);
421 /* TODO This code depends on the logic in tables.c that constructs
422 the table layout, which should be made explicit in future
426 "Pair-interaction code that uses GROMACS interaction tables supports exactly 3 tables");
428 fr->pairsTable->interaction == TableInteraction::ElectrostaticVdwRepulsionVdwDispersion,
429 "Pair interaction kernels need a table with Coulomb, repulsion and dispersion entries");
431 const real epsfac = fr->ic->epsfac;
434 for (i = 0; (i < nbonds);)
439 gid = GID(cENER[ai], cENER[aj], numEnergyGroups);
446 (fr->efep != FreeEnergyPerturbationType::No
447 && ((!atomIsPerturbed.empty() && (atomIsPerturbed[ai] || atomIsPerturbed[aj]))
448 || iparams[itype].lj14.c6A != iparams[itype].lj14.c6B
449 || iparams[itype].lj14.c12A != iparams[itype].lj14.c12B));
450 qq = chargeA[ai] * chargeA[aj] * epsfac * fr->fudgeQQ;
451 c6 = iparams[itype].lj14.c6A;
452 c12 = iparams[itype].lj14.c12A;
455 qq = iparams[itype].ljc14.qi * iparams[itype].ljc14.qj * epsfac
456 * iparams[itype].ljc14.fqq;
457 c6 = iparams[itype].ljc14.c6;
458 c12 = iparams[itype].ljc14.c12;
461 qq = iparams[itype].ljcnb.qi * iparams[itype].ljcnb.qj * epsfac;
462 c6 = iparams[itype].ljcnb.c6;
463 c12 = iparams[itype].ljcnb.c12;
466 /* Cannot happen since we called gmx_fatal() above in this case */
467 qq = c6 = c12 = 0; /* Keep compiler happy */
471 /* To save flops in the optimized kernels, c6/c12 have 6.0/12.0 derivative prefactors
472 * included in the general nfbp array now. This means the tables are scaled down by the
473 * same factor, so when we use the original c6/c12 parameters from iparams[] they must
479 /* Do we need to apply full periodic boundary conditions? */
482 fshift_index = pbc_dx_aiuc(pbc, x[ai], x[aj], dx);
486 fshift_index = c_centralShiftIndex;
487 rvec_sub(x[ai], x[aj], dx);
491 if (r2 >= fr->pairsTable->interactionRange * fr->pairsTable->interactionRange)
493 /* This check isn't race free. But it doesn't matter because if a race occurs the only
494 * disadvantage is that the warning is printed twice */
497 warning_rlimit(x, ai, aj, global_atom_index, sqrt(r2), fr->pairsTable->interactionRange);
498 warned_rlimit = TRUE;
505 /* Currently free energy is only supported for F_LJ14, so no need to check for that if we got here */
506 qqB = chargeB[ai] * chargeB[aj] * epsfac * fr->fudgeQQ;
507 c6B = iparams[itype].lj14.c6B * 6.0;
508 c12B = iparams[itype].lj14.c12B * 12.0;
510 fscal = free_energy_evaluate_single(r2,
511 *fr->ic->softCoreParameters,
512 fr->pairsTable->scale,
513 fr->pairsTable->data.data(),
514 fr->pairsTable->stride,
534 /* Evaluate tabulated interaction without free energy */
535 fscal = evaluate_single(r2,
536 fr->pairsTable->scale,
537 fr->pairsTable->data.data(),
538 fr->pairsTable->stride,
546 energygrp_elec[gid] += velec;
547 energygrp_vdw[gid] += vvdw;
548 svmul(fscal, dx, dx);
554 if (computeVirial(flavor))
556 if (fshift_index != c_centralShiftIndex)
558 rvec_inc(fshift[fshift_index], dx);
559 rvec_dec(fshift[c_centralShiftIndex], dx);
566 /*! \brief Calculate pairs, only for plain-LJ + plain Coulomb normal type.
568 * This function is templated for real/SimdReal and for optimization.
570 template<typename T, int pack_size, typename pbc_type>
571 static void do_pairs_simple(int nbonds,
572 const t_iatom iatoms[],
573 const t_iparams iparams[],
577 gmx::ArrayRef<real> charge,
578 const real scale_factor)
580 const int nfa1 = 1 + 2;
586 #if GMX_SIMD_HAVE_REAL
587 alignas(GMX_SIMD_ALIGNMENT) std::int32_t ai[pack_size];
588 alignas(GMX_SIMD_ALIGNMENT) std::int32_t aj[pack_size];
589 alignas(GMX_SIMD_ALIGNMENT) real coeff[3 * pack_size];
591 std::int32_t ai[pack_size];
592 std::int32_t aj[pack_size];
593 real coeff[3 * pack_size];
596 /* nbonds is #pairs*nfa1, here we step pack_size pairs */
597 for (int i = 0; i < nbonds; i += pack_size * nfa1)
599 /* Collect atoms for pack_size pairs.
600 * iu indexes into iatoms, we should not let iu go beyond nbonds.
603 for (int s = 0; s < pack_size; s++)
605 int itype = iatoms[iu];
606 ai[s] = iatoms[iu + 1];
607 aj[s] = iatoms[iu + 2];
609 if (i + s * nfa1 < nbonds)
611 coeff[0 * pack_size + s] = iparams[itype].lj14.c6A;
612 coeff[1 * pack_size + s] = iparams[itype].lj14.c12A;
613 coeff[2 * pack_size + s] = charge[ai[s]] * charge[aj[s]];
615 /* Avoid indexing the iatoms array out of bounds.
616 * We pad the coordinate indices with the last atom pair.
618 if (iu + nfa1 < nbonds)
625 /* Pad the coefficient arrays with zeros to get zero forces */
626 coeff[0 * pack_size + s] = 0;
627 coeff[1 * pack_size + s] = 0;
628 coeff[2 * pack_size + s] = 0;
632 /* Load the coordinates */
634 gatherLoadUTranspose<3>(reinterpret_cast<const real*>(x), ai, &xi[XX], &xi[YY], &xi[ZZ]);
635 gatherLoadUTranspose<3>(reinterpret_cast<const real*>(x), aj, &xj[XX], &xj[YY], &xj[ZZ]);
637 T c6 = load<T>(coeff + 0 * pack_size);
638 T c12 = load<T>(coeff + 1 * pack_size);
639 T qq = load<T>(coeff + 2 * pack_size);
641 /* We could save these operations by storing 6*C6,12*C12 */
646 pbc_dx_aiuc(pbc, xi, xj, dr);
648 T rsq = dr[XX] * dr[XX] + dr[YY] * dr[YY] + dr[ZZ] * dr[ZZ];
649 T rinv = gmx::invsqrt(rsq);
650 T rinv2 = rinv * rinv;
651 T rinv6 = rinv2 * rinv2 * rinv2;
653 /* Calculate the Coulomb force * r */
654 T cfr = ef * qq * rinv;
656 /* Calculate the LJ force * r and add it to the Coulomb part */
657 T fr = gmx::fma(fms(c12, rinv6, c6), rinv6, cfr);
659 T finvr = fr * rinv2;
660 T fx = finvr * dr[XX];
661 T fy = finvr * dr[YY];
662 T fz = finvr * dr[ZZ];
664 /* Add the pair forces to the force array.
665 * Note that here we might add multiple force components for some atoms
666 * due to the SIMD padding. But the extra force components are zero.
668 transposeScatterIncrU<4>(reinterpret_cast<real*>(f), ai, fx, fy, fz);
669 transposeScatterDecrU<4>(reinterpret_cast<real*>(f), aj, fx, fy, fz);
673 /*! \brief Calculate all listed pair interactions */
674 void do_pairs(int ftype,
676 const t_iatom iatoms[],
677 const t_iparams iparams[],
681 const struct t_pbc* pbc,
684 gmx::ArrayRef<real> chargeA,
685 gmx::ArrayRef<real> chargeB,
686 gmx::ArrayRef<bool> atomIsPerturbed,
687 gmx::ArrayRef<unsigned short> cENER,
688 const int numEnergyGroups,
689 const t_forcerec* fr,
690 const bool havePerturbedInteractions,
691 const gmx::StepWorkload& stepWork,
692 gmx_grppairener_t* grppener,
693 int* global_atom_index)
695 if (ftype == F_LJ14 && fr->ic->vdwtype != VanDerWaalsType::User && !EEL_USER(fr->ic->eeltype)
696 && !havePerturbedInteractions && (!stepWork.computeVirial && !stepWork.computeEnergy))
698 /* We use a fast code-path for plain LJ 1-4 without FEP.
700 * TODO: Add support for energies (straightforward) and virial
701 * in the SIMD template. For the virial it's inconvenient to store
702 * the force sums for the shifts and we should directly calculate
703 * and sum the virial for the shifts. But we should do this
704 * at once for the angles and dihedrals as well.
706 #if GMX_SIMD_HAVE_REAL
707 if (fr->use_simd_kernels)
709 alignas(GMX_SIMD_ALIGNMENT) real pbc_simd[9 * GMX_SIMD_REAL_WIDTH];
710 set_pbc_simd(pbc, pbc_simd);
712 do_pairs_simple<SimdReal, GMX_SIMD_REAL_WIDTH, const real*>(
713 nbonds, iatoms, iparams, x, f, pbc_simd, chargeA, fr->ic->epsfac * fr->fudgeQQ);
718 /* This construct is needed because pbc_dx_aiuc doesn't accept pbc=NULL */
720 const t_pbc* pbc_nonnull;
728 set_pbc(&pbc_no, PbcType::No, nullptr);
729 pbc_nonnull = &pbc_no;
732 do_pairs_simple<real, 1, const t_pbc*>(
733 nbonds, iatoms, iparams, x, f, pbc_nonnull, chargeA, fr->ic->epsfac * fr->fudgeQQ);
736 else if (stepWork.computeVirial)
738 do_pairs_general<BondedKernelFlavor::ForcesAndVirialAndEnergy>(ftype,
759 do_pairs_general<BondedKernelFlavor::ForcesAndEnergy>(ftype,