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/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), glatnr(global_atom_index, aj), r, rlimit);
95 "%8f %8f %8f\n%8f %8f %8f\n1-4 (%d,%d) interaction not within cut-off! r=%g. "
97 x[ai][XX], x[ai][YY], x[ai][ZZ], x[aj][XX], x[aj][YY], x[aj][ZZ],
98 glatnr(global_atom_index, ai), glatnr(global_atom_index, aj), r);
102 /*! \brief Compute the energy and force for a single pair interaction */
103 static real evaluate_single(real r2,
113 real rinv, r, rtab, eps, eps2, Y, F, Geps, Heps2, Fp, VVe, FFe, VVd, FFd, VVr, FFr, fscal;
116 /* Do the tabulated interactions - first table lookup */
117 rinv = gmx::invsqrt(r2);
120 ntab = static_cast<int>(rtab);
123 ntab = static_cast<int>(tableStride * ntab);
127 Geps = eps * vftab[ntab + 2];
128 Heps2 = eps2 * vftab[ntab + 3];
129 Fp = F + Geps + Heps2;
131 FFe = Fp + Geps + 2.0 * Heps2;
135 Geps = eps * vftab[ntab + 6];
136 Heps2 = eps2 * vftab[ntab + 7];
137 Fp = F + Geps + Heps2;
139 FFd = Fp + Geps + 2.0 * Heps2;
143 Geps = eps * vftab[ntab + 10];
144 Heps2 = eps2 * vftab[ntab + 11];
145 Fp = F + Geps + Heps2;
147 FFr = Fp + Geps + 2.0 * Heps2;
150 *vvdw = c6 * VVd + c12 * VVr;
152 fscal = -(qq * FFe + c6 * FFd + c12 * FFr) * tabscale * rinv;
157 static inline real sixthRoot(const real r)
159 return gmx::invsqrt(std::cbrt(r));
162 /*! \brief Compute the energy and force for a single pair interaction under FEP */
163 static real free_energy_evaluate_single(real r2,
164 const interaction_const_t::SoftCoreParameters& scParams,
177 const real lfac_coul[2],
178 const real lfac_vdw[2],
179 const real dlfac_coul[2],
180 const real dlfac_vdw[2],
185 real rtab, eps, eps2, Y, F, Geps, Heps2, Fp, VV, FF, fscal;
186 real qq[2], c6[2], c12[2], sigma6[2], sigma_pow[2];
187 real alpha_coul_eff, alpha_vdw_eff, dvdl_coul, dvdl_vdw;
188 real rpinv, r_coul, r_vdw, velecsum, vvdwsum;
189 real fscal_vdw[2], fscal_elec[2];
190 real velec[2], vvdw[2];
192 const real half = 0.5_real;
193 const real one = 1.0_real;
194 const real two = 2.0_real;
203 const real rpm2 = r2 * r2; /* r4 */
204 const real rp = rpm2 * r2; /* r6 */
206 /* Loop over state A(0) and B(1) */
207 for (i = 0; i < 2; i++)
209 if ((c6[i] > 0) && (c12[i] > 0))
211 /* The c6 & c12 coefficients now contain the constants 6.0 and 12.0, respectively.
212 * Correct for this by multiplying with (1/12.0)/(1/6.0)=6.0/12.0=0.5.
214 sigma6[i] = half * c12[i] / c6[i];
215 if (sigma6[i] < scParams.sigma6Minimum) /* for disappearing coul and vdw with soft core at the same time */
217 sigma6[i] = scParams.sigma6Minimum;
222 sigma6[i] = scParams.sigma6WithInvalidSigma;
224 sigma_pow[i] = sigma6[i];
227 /* only use softcore if one of the states has a zero endstate - softcore is for avoiding infinities!*/
228 if ((c12[0] > 0) && (c12[1] > 0))
235 alpha_vdw_eff = scParams.alphaVdw;
236 alpha_coul_eff = scParams.alphaCoulomb;
239 /* Loop over A and B states again */
240 for (i = 0; i < 2; i++)
247 /* Only spend time on A or B state if it is non-zero */
248 if ((qq[i] != 0) || (c6[i] != 0) || (c12[i] != 0))
251 rpinv = one / (alpha_coul_eff * lfac_coul[i] * sigma_pow[i] + rp);
252 r_coul = sixthRoot(rpinv);
254 /* Electrostatics table lookup data */
255 rtab = r_coul * tabscale;
256 ntab = static_cast<int>(rtab);
259 ntab = static_cast<int>(tableStride * ntab);
263 Geps = eps * vftab[ntab + 2];
264 Heps2 = eps2 * vftab[ntab + 3];
265 Fp = F + Geps + Heps2;
267 FF = Fp + Geps + two * Heps2;
268 velec[i] = qq[i] * VV;
269 fscal_elec[i] = -qq[i] * FF * r_coul * rpinv * tabscale;
272 rpinv = one / (alpha_vdw_eff * lfac_vdw[i] * sigma_pow[i] + rp);
273 r_vdw = sixthRoot(rpinv);
274 /* Vdw table lookup data */
275 rtab = r_vdw * tabscale;
276 ntab = static_cast<int>(rtab);
283 Geps = eps * vftab[ntab + 6];
284 Heps2 = eps2 * vftab[ntab + 7];
285 Fp = F + Geps + Heps2;
287 FF = Fp + Geps + two * Heps2;
288 vvdw[i] = c6[i] * VV;
289 fscal_vdw[i] = -c6[i] * FF;
294 Geps = eps * vftab[ntab + 10];
295 Heps2 = eps2 * vftab[ntab + 11];
296 Fp = F + Geps + Heps2;
298 FF = Fp + Geps + two * Heps2;
299 vvdw[i] += c12[i] * VV;
300 fscal_vdw[i] -= c12[i] * FF;
301 fscal_vdw[i] *= r_vdw * rpinv * tabscale;
304 /* Now we have velec[i], vvdw[i], and fscal[i] for both states */
305 /* Assemble A and B states */
311 for (i = 0; i < 2; i++)
313 velecsum += LFC[i] * velec[i];
314 vvdwsum += LFV[i] * vvdw[i];
316 fscal += (LFC[i] * fscal_elec[i] + LFV[i] * fscal_vdw[i]) * rpm2;
318 dvdl_coul += velec[i] * DLF[i]
319 + LFC[i] * alpha_coul_eff * dlfac_coul[i] * fscal_elec[i] * sigma_pow[i];
320 dvdl_vdw += vvdw[i] * DLF[i]
321 + LFV[i] * alpha_vdw_eff * dlfac_vdw[i] * fscal_vdw[i] * sigma_pow[i];
324 dvdl[efptCOUL] += dvdl_coul;
325 dvdl[efptVDW] += dvdl_vdw;
327 *velectot = velecsum;
333 /*! \brief Calculate pair interactions, supports all types and conditions. */
334 template<BondedKernelFlavor flavor>
335 static real do_pairs_general(int ftype,
337 const t_iatom iatoms[],
338 const t_iparams iparams[],
342 const struct t_pbc* pbc,
346 const t_forcerec* fr,
347 gmx_grppairener_t* grppener,
348 int* global_atom_index)
352 int i, itype, ai, aj, gid;
355 real fscal, velec, vvdw;
356 real* energygrp_elec;
358 static gmx_bool warned_rlimit = FALSE;
359 /* Free energy stuff */
360 gmx_bool bFreeEnergy;
361 real LFC[2], LFV[2], DLF[2], lfac_coul[2], lfac_vdw[2], dlfac_coul[2], dlfac_vdw[2];
363 const real oneSixth = 1.0_real / 6.0_real;
369 energygrp_elec = grppener->ener[egCOUL14].data();
370 energygrp_vdw = grppener->ener[egLJ14].data();
373 energygrp_elec = grppener->ener[egCOULSR].data();
374 energygrp_vdw = grppener->ener[egLJSR].data();
377 energygrp_elec = nullptr; /* Keep compiler happy */
378 energygrp_vdw = nullptr; /* Keep compiler happy */
379 gmx_fatal(FARGS, "Unknown function type %d in do_nonbonded14", ftype);
382 if (fr->efep != efepNO)
384 /* Lambda factor for state A=1-lambda and B=lambda */
385 LFC[0] = 1.0 - lambda[efptCOUL];
386 LFV[0] = 1.0 - lambda[efptVDW];
387 LFC[1] = lambda[efptCOUL];
388 LFV[1] = lambda[efptVDW];
390 /*derivative of the lambda factor for state A and B */
394 GMX_ASSERT(fr->ic->softCoreParameters, "We need soft-core parameters");
395 const auto& scParams = *fr->ic->softCoreParameters;
397 for (i = 0; i < 2; i++)
399 lfac_coul[i] = (scParams.lambdaPower == 2 ? (1 - LFC[i]) * (1 - LFC[i]) : (1 - LFC[i]));
400 dlfac_coul[i] = DLF[i] * scParams.lambdaPower * oneSixth
401 * (scParams.lambdaPower == 2 ? (1 - LFC[i]) : 1);
402 lfac_vdw[i] = (scParams.lambdaPower == 2 ? (1 - LFV[i]) * (1 - LFV[i]) : (1 - LFV[i]));
403 dlfac_vdw[i] = DLF[i] * scParams.lambdaPower * oneSixth
404 * (scParams.lambdaPower == 2 ? (1 - LFV[i]) : 1);
408 /* TODO This code depends on the logic in tables.c that constructs
409 the table layout, which should be made explicit in future
413 "Pair-interaction code that uses GROMACS interaction tables supports exactly 3 tables");
415 fr->pairsTable->interaction == GMX_TABLE_INTERACTION_ELEC_VDWREP_VDWDISP,
416 "Pair interaction kernels need a table with Coulomb, repulsion and dispersion entries");
418 const real epsfac = fr->ic->epsfac;
421 for (i = 0; (i < nbonds);)
426 gid = GID(md->cENER[ai], md->cENER[aj], md->nenergrp);
432 bFreeEnergy = (fr->efep != efepNO
433 && (((md->nPerturbed != 0) && (md->bPerturbed[ai] || md->bPerturbed[aj]))
434 || iparams[itype].lj14.c6A != iparams[itype].lj14.c6B
435 || iparams[itype].lj14.c12A != iparams[itype].lj14.c12B));
436 qq = md->chargeA[ai] * md->chargeA[aj] * epsfac * fr->fudgeQQ;
437 c6 = iparams[itype].lj14.c6A;
438 c12 = iparams[itype].lj14.c12A;
441 qq = iparams[itype].ljc14.qi * iparams[itype].ljc14.qj * epsfac
442 * iparams[itype].ljc14.fqq;
443 c6 = iparams[itype].ljc14.c6;
444 c12 = iparams[itype].ljc14.c12;
447 qq = iparams[itype].ljcnb.qi * iparams[itype].ljcnb.qj * epsfac;
448 c6 = iparams[itype].ljcnb.c6;
449 c12 = iparams[itype].ljcnb.c12;
452 /* Cannot happen since we called gmx_fatal() above in this case */
453 qq = c6 = c12 = 0; /* Keep compiler happy */
457 /* To save flops in the optimized kernels, c6/c12 have 6.0/12.0 derivative prefactors
458 * included in the general nfbp array now. This means the tables are scaled down by the
459 * same factor, so when we use the original c6/c12 parameters from iparams[] they must
465 /* Do we need to apply full periodic boundary conditions? */
468 fshift_index = pbc_dx_aiuc(pbc, x[ai], x[aj], dx);
472 fshift_index = CENTRAL;
473 rvec_sub(x[ai], x[aj], dx);
477 if (r2 >= fr->pairsTable->r * fr->pairsTable->r)
479 /* This check isn't race free. But it doesn't matter because if a race occurs the only
480 * disadvantage is that the warning is printed twice */
483 warning_rlimit(x, ai, aj, global_atom_index, sqrt(r2), fr->pairsTable->r);
484 warned_rlimit = TRUE;
491 /* Currently free energy is only supported for F_LJ14, so no need to check for that if we got here */
492 qqB = md->chargeB[ai] * md->chargeB[aj] * epsfac * fr->fudgeQQ;
493 c6B = iparams[itype].lj14.c6B * 6.0;
494 c12B = iparams[itype].lj14.c12B * 12.0;
496 fscal = free_energy_evaluate_single(
497 r2, *fr->ic->softCoreParameters, fr->pairsTable->scale, fr->pairsTable->data.data(),
498 fr->pairsTable->stride, qq, c6, c12, qqB, c6B, c12B, LFC, LFV, DLF, lfac_coul,
499 lfac_vdw, dlfac_coul, dlfac_vdw, &velec, &vvdw, dvdl);
503 /* Evaluate tabulated interaction without free energy */
504 fscal = evaluate_single(r2, fr->pairsTable->scale, fr->pairsTable->data.data(),
505 fr->pairsTable->stride, qq, c6, c12, &velec, &vvdw);
508 energygrp_elec[gid] += velec;
509 energygrp_vdw[gid] += vvdw;
510 svmul(fscal, dx, dx);
516 if (computeVirial(flavor))
518 if (fshift_index != CENTRAL)
520 rvec_inc(fshift[fshift_index], dx);
521 rvec_dec(fshift[CENTRAL], dx);
528 /*! \brief Calculate pairs, only for plain-LJ + plain Coulomb normal type.
530 * This function is templated for real/SimdReal and for optimization.
532 template<typename T, int pack_size, typename pbc_type>
533 static void do_pairs_simple(int nbonds,
534 const t_iatom iatoms[],
535 const t_iparams iparams[],
540 const real scale_factor)
542 const int nfa1 = 1 + 2;
548 #if GMX_SIMD_HAVE_REAL
549 alignas(GMX_SIMD_ALIGNMENT) std::int32_t ai[pack_size];
550 alignas(GMX_SIMD_ALIGNMENT) std::int32_t aj[pack_size];
551 alignas(GMX_SIMD_ALIGNMENT) real coeff[3 * pack_size];
553 std::int32_t ai[pack_size];
554 std::int32_t aj[pack_size];
555 real coeff[3 * pack_size];
558 /* nbonds is #pairs*nfa1, here we step pack_size pairs */
559 for (int i = 0; i < nbonds; i += pack_size * nfa1)
561 /* Collect atoms for pack_size pairs.
562 * iu indexes into iatoms, we should not let iu go beyond nbonds.
565 for (int s = 0; s < pack_size; s++)
567 int itype = iatoms[iu];
568 ai[s] = iatoms[iu + 1];
569 aj[s] = iatoms[iu + 2];
571 if (i + s * nfa1 < nbonds)
573 coeff[0 * pack_size + s] = iparams[itype].lj14.c6A;
574 coeff[1 * pack_size + s] = iparams[itype].lj14.c12A;
575 coeff[2 * pack_size + s] = md->chargeA[ai[s]] * md->chargeA[aj[s]];
577 /* Avoid indexing the iatoms array out of bounds.
578 * We pad the coordinate indices with the last atom pair.
580 if (iu + nfa1 < nbonds)
587 /* Pad the coefficient arrays with zeros to get zero forces */
588 coeff[0 * pack_size + s] = 0;
589 coeff[1 * pack_size + s] = 0;
590 coeff[2 * pack_size + s] = 0;
594 /* Load the coordinates */
596 gatherLoadUTranspose<3>(reinterpret_cast<const real*>(x), ai, &xi[XX], &xi[YY], &xi[ZZ]);
597 gatherLoadUTranspose<3>(reinterpret_cast<const real*>(x), aj, &xj[XX], &xj[YY], &xj[ZZ]);
599 T c6 = load<T>(coeff + 0 * pack_size);
600 T c12 = load<T>(coeff + 1 * pack_size);
601 T qq = load<T>(coeff + 2 * pack_size);
603 /* We could save these operations by storing 6*C6,12*C12 */
608 pbc_dx_aiuc(pbc, xi, xj, dr);
610 T rsq = dr[XX] * dr[XX] + dr[YY] * dr[YY] + dr[ZZ] * dr[ZZ];
611 T rinv = gmx::invsqrt(rsq);
612 T rinv2 = rinv * rinv;
613 T rinv6 = rinv2 * rinv2 * rinv2;
615 /* Calculate the Coulomb force * r */
616 T cfr = ef * qq * rinv;
618 /* Calculate the LJ force * r and add it to the Coulomb part */
619 T fr = gmx::fma(fms(c12, rinv6, c6), rinv6, cfr);
621 T finvr = fr * rinv2;
622 T fx = finvr * dr[XX];
623 T fy = finvr * dr[YY];
624 T fz = finvr * dr[ZZ];
626 /* Add the pair forces to the force array.
627 * Note that here we might add multiple force components for some atoms
628 * due to the SIMD padding. But the extra force components are zero.
630 transposeScatterIncrU<4>(reinterpret_cast<real*>(f), ai, fx, fy, fz);
631 transposeScatterDecrU<4>(reinterpret_cast<real*>(f), aj, fx, fy, fz);
635 /*! \brief Calculate all listed pair interactions */
636 void do_pairs(int ftype,
638 const t_iatom iatoms[],
639 const t_iparams iparams[],
643 const struct t_pbc* pbc,
647 const t_forcerec* fr,
648 const bool havePerturbedInteractions,
649 const gmx::StepWorkload& stepWork,
650 gmx_grppairener_t* grppener,
651 int* global_atom_index)
653 if (ftype == F_LJ14 && fr->ic->vdwtype != evdwUSER && !EEL_USER(fr->ic->eeltype)
654 && !havePerturbedInteractions && (!stepWork.computeVirial && !stepWork.computeEnergy))
656 /* We use a fast code-path for plain LJ 1-4 without FEP.
658 * TODO: Add support for energies (straightforward) and virial
659 * in the SIMD template. For the virial it's inconvenient to store
660 * the force sums for the shifts and we should directly calculate
661 * and sum the virial for the shifts. But we should do this
662 * at once for the angles and dihedrals as well.
664 #if GMX_SIMD_HAVE_REAL
665 if (fr->use_simd_kernels)
667 alignas(GMX_SIMD_ALIGNMENT) real pbc_simd[9 * GMX_SIMD_REAL_WIDTH];
668 set_pbc_simd(pbc, pbc_simd);
670 do_pairs_simple<SimdReal, GMX_SIMD_REAL_WIDTH, const real*>(
671 nbonds, iatoms, iparams, x, f, pbc_simd, md, fr->ic->epsfac * fr->fudgeQQ);
676 /* This construct is needed because pbc_dx_aiuc doesn't accept pbc=NULL */
678 const t_pbc* pbc_nonnull;
686 set_pbc(&pbc_no, PbcType::No, nullptr);
687 pbc_nonnull = &pbc_no;
690 do_pairs_simple<real, 1, const t_pbc*>(nbonds, iatoms, iparams, x, f, pbc_nonnull, md,
691 fr->ic->epsfac * fr->fudgeQQ);
694 else if (stepWork.computeVirial)
696 do_pairs_general<BondedKernelFlavor::ForcesAndVirialAndEnergy>(
697 ftype, nbonds, iatoms, iparams, x, f, fshift, pbc, lambda, dvdl, md, fr, grppener,
702 do_pairs_general<BondedKernelFlavor::ForcesAndEnergy>(ftype, nbonds, iatoms, iparams, x, f,
703 fshift, pbc, lambda, dvdl, md, fr,
704 grppener, global_atom_index);