2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2014,2015,2016,2017,2018,2019, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
37 * \brief This file defines functions for "pair" interactions
38 * (i.e. listed non-bonded interactions, e.g. 1-4 interactions)
40 * \author Mark Abraham <mark.j.abraham@gmail.com>
42 * \ingroup module_listed_forces
50 #include "gromacs/listed_forces/bonded.h"
51 #include "gromacs/math/functions.h"
52 #include "gromacs/math/vec.h"
53 #include "gromacs/mdtypes/group.h"
54 #include "gromacs/mdtypes/md_enums.h"
55 #include "gromacs/mdtypes/nblist.h"
56 #include "gromacs/mdtypes/simulation_workload.h"
57 #include "gromacs/pbcutil/ishift.h"
58 #include "gromacs/pbcutil/mshift.h"
59 #include "gromacs/pbcutil/pbc.h"
60 #include "gromacs/pbcutil/pbc_simd.h"
61 #include "gromacs/simd/simd.h"
62 #include "gromacs/simd/simd_math.h"
63 #include "gromacs/simd/vector_operations.h"
64 #include "gromacs/tables/forcetable.h"
65 #include "gromacs/utility/basedefinitions.h"
66 #include "gromacs/utility/fatalerror.h"
67 #include "gromacs/utility/gmxassert.h"
69 #include "listed_internal.h"
71 using namespace gmx; // TODO: Remove when this file is moved into gmx namespace
73 /*! \brief Issue a warning if a listed interaction is beyond a table limit */
74 static void warning_rlimit(const rvec* x, int ai, int aj, int* global_atom_index, real r, real rlimit)
77 "Listed nonbonded interaction between particles %d and %d\n"
78 "at distance %.3f which is larger than the table limit %.3f nm.\n\n"
79 "This is likely either a 1,4 interaction, or a listed interaction inside\n"
80 "a smaller molecule you are decoupling during a free energy calculation.\n"
81 "Since interactions at distances beyond the table cannot be computed,\n"
82 "they are skipped until they are inside the table limit again. You will\n"
83 "only see this message once, even if it occurs for several interactions.\n\n"
84 "IMPORTANT: This should not happen in a stable simulation, so there is\n"
85 "probably something wrong with your system. Only change the table-extension\n"
86 "distance in the mdp file if you are really sure that is the reason.\n",
87 glatnr(global_atom_index, ai), glatnr(global_atom_index, aj), r, rlimit);
92 "%8f %8f %8f\n%8f %8f %8f\n1-4 (%d,%d) interaction not within cut-off! r=%g. "
94 x[ai][XX], x[ai][YY], x[ai][ZZ], x[aj][XX], x[aj][YY], x[aj][ZZ],
95 glatnr(global_atom_index, ai), glatnr(global_atom_index, aj), r);
99 /*! \brief Compute the energy and force for a single pair interaction */
100 static real evaluate_single(real r2,
110 real rinv, r, rtab, eps, eps2, Y, F, Geps, Heps2, Fp, VVe, FFe, VVd, FFd, VVr, FFr, fscal;
113 /* Do the tabulated interactions - first table lookup */
114 rinv = gmx::invsqrt(r2);
117 ntab = static_cast<int>(rtab);
120 ntab = static_cast<int>(tableStride * ntab);
124 Geps = eps * vftab[ntab + 2];
125 Heps2 = eps2 * vftab[ntab + 3];
126 Fp = F + Geps + Heps2;
128 FFe = Fp + Geps + 2.0 * Heps2;
132 Geps = eps * vftab[ntab + 6];
133 Heps2 = eps2 * vftab[ntab + 7];
134 Fp = F + Geps + Heps2;
136 FFd = Fp + Geps + 2.0 * Heps2;
140 Geps = eps * vftab[ntab + 10];
141 Heps2 = eps2 * vftab[ntab + 11];
142 Fp = F + Geps + Heps2;
144 FFr = Fp + Geps + 2.0 * Heps2;
147 *vvdw = c6 * VVd + c12 * VVr;
149 fscal = -(qq * FFe + c6 * FFd + c12 * FFr) * tabscale * rinv;
154 /*! \brief Compute the energy and force for a single pair interaction under FEP */
155 static real free_energy_evaluate_single(real r2,
171 const real lfac_coul[2],
172 const real lfac_vdw[2],
173 const real dlfac_coul[2],
174 const real dlfac_vdw[2],
183 real rp, rpm2, rtab, eps, eps2, Y, F, Geps, Heps2, Fp, VV, FF, fscal;
184 real qq[2], c6[2], c12[2], sigma6[2], sigma2[2], sigma_pow[2];
185 real alpha_coul_eff, alpha_vdw_eff, dvdl_coul, dvdl_vdw;
186 real rpinv, r_coul, r_vdw, velecsum, vvdwsum;
187 real fscal_vdw[2], fscal_elec[2];
188 real velec[2], vvdw[2];
190 const real half = 0.5;
191 const real minusOne = -1.0;
192 const real one = 1.0;
193 const real two = 2.0;
194 const real six = 6.0;
195 const real fourtyeight = 48.0;
204 if (sc_r_power == six)
206 rpm2 = r2 * r2; /* r4 */
207 rp = rpm2 * r2; /* r6 */
209 else if (sc_r_power == fourtyeight)
211 rp = r2 * r2 * r2; /* r6 */
212 rp = rp * rp; /* r12 */
213 rp = rp * rp; /* r24 */
214 rp = rp * rp; /* r48 */
215 rpm2 = rp / r2; /* r46 */
219 rp = std::pow(r2, half * sc_r_power); /* not currently supported as input, but can handle it */
223 /* Loop over state A(0) and B(1) */
224 for (i = 0; i < 2; i++)
226 if ((c6[i] > 0) && (c12[i] > 0))
228 /* The c6 & c12 coefficients now contain the constants 6.0 and 12.0, respectively.
229 * Correct for this by multiplying with (1/12.0)/(1/6.0)=6.0/12.0=0.5.
231 sigma6[i] = half * c12[i] / c6[i];
232 sigma2[i] = std::cbrt(half * c12[i] / c6[i]);
233 /* should be able to get rid of this ^^^ internal pow call eventually. Will require agreement on
234 what data to store externally. Can't be fixed without larger scale changes, so not 5.0 */
235 if (sigma6[i] < sigma6_min) /* for disappearing coul and vdw with soft core at the same time */
237 sigma6[i] = sigma6_min;
238 sigma2[i] = sigma2_min;
243 sigma6[i] = sigma6_def;
244 sigma2[i] = sigma2_def;
246 if (sc_r_power == six)
248 sigma_pow[i] = sigma6[i];
250 else if (sc_r_power == fourtyeight)
252 sigma_pow[i] = sigma6[i] * sigma6[i]; /* sigma^12 */
253 sigma_pow[i] = sigma_pow[i] * sigma_pow[i]; /* sigma^24 */
254 sigma_pow[i] = sigma_pow[i] * sigma_pow[i]; /* sigma^48 */
257 { /* not really supported as input, but in here for testing the general case*/
258 sigma_pow[i] = std::pow(sigma2[i], sc_r_power / 2);
262 /* only use softcore if one of the states has a zero endstate - softcore is for avoiding infinities!*/
263 if ((c12[0] > 0) && (c12[1] > 0))
270 alpha_vdw_eff = alpha_vdw;
271 alpha_coul_eff = alpha_coul;
274 /* Loop over A and B states again */
275 for (i = 0; i < 2; i++)
282 /* Only spend time on A or B state if it is non-zero */
283 if ((qq[i] != 0) || (c6[i] != 0) || (c12[i] != 0))
286 rpinv = one / (alpha_coul_eff * lfac_coul[i] * sigma_pow[i] + rp);
287 r_coul = std::pow(rpinv, minusOne / sc_r_power);
289 /* Electrostatics table lookup data */
290 rtab = r_coul * tabscale;
291 ntab = static_cast<int>(rtab);
294 ntab = static_cast<int>(tableStride * ntab);
298 Geps = eps * vftab[ntab + 2];
299 Heps2 = eps2 * vftab[ntab + 3];
300 Fp = F + Geps + Heps2;
302 FF = Fp + Geps + two * Heps2;
303 velec[i] = qq[i] * VV;
304 fscal_elec[i] = -qq[i] * FF * r_coul * rpinv * tabscale;
307 rpinv = one / (alpha_vdw_eff * lfac_vdw[i] * sigma_pow[i] + rp);
308 r_vdw = std::pow(rpinv, minusOne / sc_r_power);
309 /* Vdw table lookup data */
310 rtab = r_vdw * tabscale;
311 ntab = static_cast<int>(rtab);
318 Geps = eps * vftab[ntab + 6];
319 Heps2 = eps2 * vftab[ntab + 7];
320 Fp = F + Geps + Heps2;
322 FF = Fp + Geps + two * Heps2;
323 vvdw[i] = c6[i] * VV;
324 fscal_vdw[i] = -c6[i] * FF;
329 Geps = eps * vftab[ntab + 10];
330 Heps2 = eps2 * vftab[ntab + 11];
331 Fp = F + Geps + Heps2;
333 FF = Fp + Geps + two * Heps2;
334 vvdw[i] += c12[i] * VV;
335 fscal_vdw[i] -= c12[i] * FF;
336 fscal_vdw[i] *= r_vdw * rpinv * tabscale;
339 /* Now we have velec[i], vvdw[i], and fscal[i] for both states */
340 /* Assemble A and B states */
346 for (i = 0; i < 2; i++)
348 velecsum += LFC[i] * velec[i];
349 vvdwsum += LFV[i] * vvdw[i];
351 fscal += (LFC[i] * fscal_elec[i] + LFV[i] * fscal_vdw[i]) * rpm2;
353 dvdl_coul += velec[i] * DLF[i]
354 + LFC[i] * alpha_coul_eff * dlfac_coul[i] * fscal_elec[i] * sigma_pow[i];
355 dvdl_vdw += vvdw[i] * DLF[i]
356 + LFV[i] * alpha_vdw_eff * dlfac_vdw[i] * fscal_vdw[i] * sigma_pow[i];
359 dvdl[efptCOUL] += dvdl_coul;
360 dvdl[efptVDW] += dvdl_vdw;
362 *velectot = velecsum;
368 /*! \brief Calculate pair interactions, supports all types and conditions. */
369 template<BondedKernelFlavor flavor>
370 static real do_pairs_general(int ftype,
372 const t_iatom iatoms[],
373 const t_iparams iparams[],
377 const struct t_pbc* pbc,
378 const struct t_graph* g,
382 const t_forcerec* fr,
383 gmx_grppairener_t* grppener,
384 int* global_atom_index)
389 int i, itype, ai, aj, gid;
392 real fscal, velec, vvdw;
393 real* energygrp_elec;
395 static gmx_bool warned_rlimit = FALSE;
396 /* Free energy stuff */
397 gmx_bool bFreeEnergy;
398 real LFC[2], LFV[2], DLF[2], lfac_coul[2], lfac_vdw[2], dlfac_coul[2], dlfac_vdw[2];
399 real qqB, c6B, c12B, sigma2_def, sigma2_min;
405 energygrp_elec = grppener->ener[egCOUL14].data();
406 energygrp_vdw = grppener->ener[egLJ14].data();
409 energygrp_elec = grppener->ener[egCOULSR].data();
410 energygrp_vdw = grppener->ener[egLJSR].data();
413 energygrp_elec = nullptr; /* Keep compiler happy */
414 energygrp_vdw = nullptr; /* Keep compiler happy */
415 gmx_fatal(FARGS, "Unknown function type %d in do_nonbonded14", ftype);
418 if (fr->efep != efepNO)
420 /* Lambda factor for state A=1-lambda and B=lambda */
421 LFC[0] = 1.0 - lambda[efptCOUL];
422 LFV[0] = 1.0 - lambda[efptVDW];
423 LFC[1] = lambda[efptCOUL];
424 LFV[1] = lambda[efptVDW];
426 /*derivative of the lambda factor for state A and B */
431 sigma2_def = std::cbrt(fr->sc_sigma6_def);
432 sigma2_min = std::cbrt(fr->sc_sigma6_min);
434 for (i = 0; i < 2; i++)
436 lfac_coul[i] = (fr->sc_power == 2 ? (1 - LFC[i]) * (1 - LFC[i]) : (1 - LFC[i]));
438 DLF[i] * fr->sc_power / fr->sc_r_power * (fr->sc_power == 2 ? (1 - LFC[i]) : 1);
439 lfac_vdw[i] = (fr->sc_power == 2 ? (1 - LFV[i]) * (1 - LFV[i]) : (1 - LFV[i]));
441 DLF[i] * fr->sc_power / fr->sc_r_power * (fr->sc_power == 2 ? (1 - LFV[i]) : 1);
446 sigma2_min = sigma2_def = 0;
449 /* TODO This code depends on the logic in tables.c that constructs
450 the table layout, which should be made explicit in future
454 "Pair-interaction code that uses GROMACS interaction tables supports exactly 3 tables");
456 fr->pairsTable->interaction == GMX_TABLE_INTERACTION_ELEC_VDWREP_VDWDISP,
457 "Pair interaction kernels need a table with Coulomb, repulsion and dispersion entries");
459 const real epsfac = fr->ic->epsfac;
462 for (i = 0; (i < nbonds);)
467 gid = GID(md->cENER[ai], md->cENER[aj], md->nenergrp);
473 bFreeEnergy = (fr->efep != efepNO
474 && (((md->nPerturbed != 0) && (md->bPerturbed[ai] || md->bPerturbed[aj]))
475 || iparams[itype].lj14.c6A != iparams[itype].lj14.c6B
476 || iparams[itype].lj14.c12A != iparams[itype].lj14.c12B));
477 qq = md->chargeA[ai] * md->chargeA[aj] * epsfac * fr->fudgeQQ;
478 c6 = iparams[itype].lj14.c6A;
479 c12 = iparams[itype].lj14.c12A;
482 qq = iparams[itype].ljc14.qi * iparams[itype].ljc14.qj * epsfac
483 * iparams[itype].ljc14.fqq;
484 c6 = iparams[itype].ljc14.c6;
485 c12 = iparams[itype].ljc14.c12;
488 qq = iparams[itype].ljcnb.qi * iparams[itype].ljcnb.qj * epsfac;
489 c6 = iparams[itype].ljcnb.c6;
490 c12 = iparams[itype].ljcnb.c12;
493 /* Cannot happen since we called gmx_fatal() above in this case */
494 qq = c6 = c12 = 0; /* Keep compiler happy */
498 /* To save flops in the optimized kernels, c6/c12 have 6.0/12.0 derivative prefactors
499 * included in the general nfbp array now. This means the tables are scaled down by the
500 * same factor, so when we use the original c6/c12 parameters from iparams[] they must
506 /* Do we need to apply full periodic boundary conditions? */
509 fshift_index = pbc_dx_aiuc(pbc, x[ai], x[aj], dx);
513 fshift_index = CENTRAL;
514 rvec_sub(x[ai], x[aj], dx);
518 if (r2 >= fr->pairsTable->r * fr->pairsTable->r)
520 /* This check isn't race free. But it doesn't matter because if a race occurs the only
521 * disadvantage is that the warning is printed twice */
524 warning_rlimit(x, ai, aj, global_atom_index, sqrt(r2), fr->pairsTable->r);
525 warned_rlimit = TRUE;
532 /* Currently free energy is only supported for F_LJ14, so no need to check for that if we got here */
533 qqB = md->chargeB[ai] * md->chargeB[aj] * epsfac * fr->fudgeQQ;
534 c6B = iparams[itype].lj14.c6B * 6.0;
535 c12B = iparams[itype].lj14.c12B * 12.0;
537 fscal = free_energy_evaluate_single(
538 r2, fr->sc_r_power, fr->sc_alphacoul, fr->sc_alphavdw, fr->pairsTable->scale,
539 fr->pairsTable->data, fr->pairsTable->stride, qq, c6, c12, qqB, c6B, c12B, LFC,
540 LFV, DLF, lfac_coul, lfac_vdw, dlfac_coul, dlfac_vdw, fr->sc_sigma6_def,
541 fr->sc_sigma6_min, sigma2_def, sigma2_min, &velec, &vvdw, dvdl);
545 /* Evaluate tabulated interaction without free energy */
546 fscal = evaluate_single(r2, fr->pairsTable->scale, fr->pairsTable->data,
547 fr->pairsTable->stride, qq, c6, c12, &velec, &vvdw);
550 energygrp_elec[gid] += velec;
551 energygrp_vdw[gid] += vvdw;
552 svmul(fscal, dx, dx);
558 if (computeVirial(flavor))
562 /* Correct the shift forces using the graph */
563 ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, aj), dt);
564 fshift_index = IVEC2IS(dt);
566 if (fshift_index != CENTRAL)
568 rvec_inc(fshift[fshift_index], dx);
569 rvec_dec(fshift[CENTRAL], dx);
576 /*! \brief Calculate pairs, only for plain-LJ + plain Coulomb normal type.
578 * This function is templated for real/SimdReal and for optimization.
580 template<typename T, int pack_size, typename pbc_type>
581 static void do_pairs_simple(int nbonds,
582 const t_iatom iatoms[],
583 const t_iparams iparams[],
588 const real scale_factor)
590 const int nfa1 = 1 + 2;
596 #if GMX_SIMD_HAVE_REAL
597 alignas(GMX_SIMD_ALIGNMENT) std::int32_t ai[pack_size];
598 alignas(GMX_SIMD_ALIGNMENT) std::int32_t aj[pack_size];
599 alignas(GMX_SIMD_ALIGNMENT) real coeff[3 * pack_size];
601 std::int32_t ai[pack_size];
602 std::int32_t aj[pack_size];
603 real coeff[3 * pack_size];
606 /* nbonds is #pairs*nfa1, here we step pack_size pairs */
607 for (int i = 0; i < nbonds; i += pack_size * nfa1)
609 /* Collect atoms for pack_size pairs.
610 * iu indexes into iatoms, we should not let iu go beyond nbonds.
613 for (int s = 0; s < pack_size; s++)
615 int itype = iatoms[iu];
616 ai[s] = iatoms[iu + 1];
617 aj[s] = iatoms[iu + 2];
619 if (i + s * nfa1 < nbonds)
621 coeff[0 * pack_size + s] = iparams[itype].lj14.c6A;
622 coeff[1 * pack_size + s] = iparams[itype].lj14.c12A;
623 coeff[2 * pack_size + s] = md->chargeA[ai[s]] * md->chargeA[aj[s]];
625 /* Avoid indexing the iatoms array out of bounds.
626 * We pad the coordinate indices with the last atom pair.
628 if (iu + nfa1 < nbonds)
635 /* Pad the coefficient arrays with zeros to get zero forces */
636 coeff[0 * pack_size + s] = 0;
637 coeff[1 * pack_size + s] = 0;
638 coeff[2 * pack_size + s] = 0;
642 /* Load the coordinates */
644 gatherLoadUTranspose<3>(reinterpret_cast<const real*>(x), ai, &xi[XX], &xi[YY], &xi[ZZ]);
645 gatherLoadUTranspose<3>(reinterpret_cast<const real*>(x), aj, &xj[XX], &xj[YY], &xj[ZZ]);
647 T c6 = load<T>(coeff + 0 * pack_size);
648 T c12 = load<T>(coeff + 1 * pack_size);
649 T qq = load<T>(coeff + 2 * pack_size);
651 /* We could save these operations by storing 6*C6,12*C12 */
656 pbc_dx_aiuc(pbc, xi, xj, dr);
658 T rsq = dr[XX] * dr[XX] + dr[YY] * dr[YY] + dr[ZZ] * dr[ZZ];
659 T rinv = gmx::invsqrt(rsq);
660 T rinv2 = rinv * rinv;
661 T rinv6 = rinv2 * rinv2 * rinv2;
663 /* Calculate the Coulomb force * r */
664 T cfr = ef * qq * rinv;
666 /* Calculate the LJ force * r and add it to the Coulomb part */
667 T fr = gmx::fma(fms(c12, rinv6, c6), rinv6, cfr);
669 T finvr = fr * rinv2;
670 T fx = finvr * dr[XX];
671 T fy = finvr * dr[YY];
672 T fz = finvr * dr[ZZ];
674 /* Add the pair forces to the force array.
675 * Note that here we might add multiple force components for some atoms
676 * due to the SIMD padding. But the extra force components are zero.
678 transposeScatterIncrU<4>(reinterpret_cast<real*>(f), ai, fx, fy, fz);
679 transposeScatterDecrU<4>(reinterpret_cast<real*>(f), aj, fx, fy, fz);
683 /*! \brief Calculate all listed pair interactions */
684 void do_pairs(int ftype,
686 const t_iatom iatoms[],
687 const t_iparams iparams[],
691 const struct t_pbc* pbc,
692 const struct t_graph* g,
696 const t_forcerec* fr,
697 const bool havePerturbedInteractions,
698 const gmx::StepWorkload& stepWork,
699 gmx_grppairener_t* grppener,
700 int* global_atom_index)
702 if (ftype == F_LJ14 && fr->ic->vdwtype != evdwUSER && !EEL_USER(fr->ic->eeltype)
703 && !havePerturbedInteractions && (!stepWork.computeVirial && !stepWork.computeEnergy))
705 /* We use a fast code-path for plain LJ 1-4 without FEP.
707 * TODO: Add support for energies (straightforward) and virial
708 * in the SIMD template. For the virial it's inconvenient to store
709 * the force sums for the shifts and we should directly calculate
710 * and sum the virial for the shifts. But we should do this
711 * at once for the angles and dihedrals as well.
713 #if GMX_SIMD_HAVE_REAL
714 if (fr->use_simd_kernels)
716 alignas(GMX_SIMD_ALIGNMENT) real pbc_simd[9 * GMX_SIMD_REAL_WIDTH];
717 set_pbc_simd(pbc, pbc_simd);
719 do_pairs_simple<SimdReal, GMX_SIMD_REAL_WIDTH, const real*>(
720 nbonds, iatoms, iparams, x, f, pbc_simd, md, fr->ic->epsfac * fr->fudgeQQ);
725 /* This construct is needed because pbc_dx_aiuc doesn't accept pbc=NULL */
727 const t_pbc* pbc_nonnull;
735 set_pbc(&pbc_no, epbcNONE, nullptr);
736 pbc_nonnull = &pbc_no;
739 do_pairs_simple<real, 1, const t_pbc*>(nbonds, iatoms, iparams, x, f, pbc_nonnull, md,
740 fr->ic->epsfac * fr->fudgeQQ);
743 else if (stepWork.computeVirial)
745 do_pairs_general<BondedKernelFlavor::ForcesAndVirialAndEnergy>(
746 ftype, nbonds, iatoms, iparams, x, f, fshift, pbc, g, lambda, dvdl, md, fr,
747 grppener, global_atom_index);
751 do_pairs_general<BondedKernelFlavor::ForcesAndEnergy>(ftype, nbonds, iatoms, iparams, x, f,
752 fshift, pbc, g, lambda, dvdl, md, fr,
753 grppener, global_atom_index);