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/math/functions.h"
51 #include "gromacs/math/vec.h"
52 #include "gromacs/mdtypes/group.h"
53 #include "gromacs/mdtypes/md_enums.h"
54 #include "gromacs/mdtypes/nblist.h"
55 #include "gromacs/pbcutil/ishift.h"
56 #include "gromacs/pbcutil/mshift.h"
57 #include "gromacs/pbcutil/pbc.h"
58 #include "gromacs/pbcutil/pbc_simd.h"
59 #include "gromacs/simd/simd.h"
60 #include "gromacs/simd/simd_math.h"
61 #include "gromacs/simd/vector_operations.h"
62 #include "gromacs/tables/forcetable.h"
63 #include "gromacs/utility/basedefinitions.h"
64 #include "gromacs/utility/fatalerror.h"
65 #include "gromacs/utility/gmxassert.h"
67 #include "listed_internal.h"
69 using namespace gmx; // TODO: Remove when this file is moved into gmx namespace
71 /*! \brief Issue a warning if a listed interaction is beyond a table limit */
73 warning_rlimit(const rvec *x, int ai, int aj, int * global_atom_index, real r, real rlimit)
75 gmx_warning("Listed nonbonded interaction between particles %d and %d\n"
76 "at distance %.3f which is larger than the table limit %.3f nm.\n\n"
77 "This is likely either a 1,4 interaction, or a listed interaction inside\n"
78 "a smaller molecule you are decoupling during a free energy calculation.\n"
79 "Since interactions at distances beyond the table cannot be computed,\n"
80 "they are skipped until they are inside the table limit again. You will\n"
81 "only see this message once, even if it occurs for several interactions.\n\n"
82 "IMPORTANT: This should not happen in a stable simulation, so there is\n"
83 "probably something wrong with your system. Only change the table-extension\n"
84 "distance in the mdp file if you are really sure that is the reason.\n",
85 glatnr(global_atom_index, ai), glatnr(global_atom_index, aj), r, rlimit);
90 "%8f %8f %8f\n%8f %8f %8f\n1-4 (%d,%d) interaction not within cut-off! r=%g. Ignored\n",
91 x[ai][XX], x[ai][YY], x[ai][ZZ], x[aj][XX], x[aj][YY], x[aj][ZZ],
92 glatnr(global_atom_index, ai), glatnr(global_atom_index, aj), r);
96 /*! \brief Compute the energy and force for a single pair interaction */
98 evaluate_single(real r2, real tabscale, const real *vftab, real tableStride,
99 real qq, real c6, real c12, real *velec, real *vvdw)
101 real rinv, r, rtab, eps, eps2, Y, F, Geps, Heps2, Fp, VVe, FFe, VVd, FFd, VVr, FFr, fscal;
104 /* Do the tabulated interactions - first table lookup */
105 rinv = gmx::invsqrt(r2);
108 ntab = static_cast<int>(rtab);
111 ntab = static_cast<int>(tableStride*ntab);
115 Geps = eps*vftab[ntab+2];
116 Heps2 = eps2*vftab[ntab+3];
119 FFe = Fp+Geps+2.0*Heps2;
123 Geps = eps*vftab[ntab+6];
124 Heps2 = eps2*vftab[ntab+7];
127 FFd = Fp+Geps+2.0*Heps2;
131 Geps = eps*vftab[ntab+10];
132 Heps2 = eps2*vftab[ntab+11];
135 FFr = Fp+Geps+2.0*Heps2;
138 *vvdw = c6*VVd+c12*VVr;
140 fscal = -(qq*FFe+c6*FFd+c12*FFr)*tabscale*rinv;
145 /*! \brief Compute the energy and force for a single pair interaction under FEP */
147 free_energy_evaluate_single(real r2, real sc_r_power, real alpha_coul,
148 real alpha_vdw, real tabscale, const real *vftab, real tableStride,
149 real qqA, real c6A, real c12A, real qqB,
150 real c6B, real c12B, const real LFC[2], const real LFV[2], const real DLF[2],
151 const real lfac_coul[2], const real lfac_vdw[2], const real dlfac_coul[2],
152 const real dlfac_vdw[2], real sigma6_def, real sigma6_min,
153 real sigma2_def, real sigma2_min,
154 real *velectot, real *vvdwtot, real *dvdl)
156 real rp, rpm2, rtab, eps, eps2, Y, F, Geps, Heps2, Fp, VV, FF, fscal;
157 real qq[2], c6[2], c12[2], sigma6[2], sigma2[2], sigma_pow[2];
158 real alpha_coul_eff, alpha_vdw_eff, dvdl_coul, dvdl_vdw;
159 real rpinv, r_coul, r_vdw, velecsum, vvdwsum;
160 real fscal_vdw[2], fscal_elec[2];
161 real velec[2], vvdw[2];
163 const real half = 0.5;
164 const real minusOne = -1.0;
165 const real one = 1.0;
166 const real two = 2.0;
167 const real six = 6.0;
168 const real fourtyeight = 48.0;
177 if (sc_r_power == six)
179 rpm2 = r2*r2; /* r4 */
180 rp = rpm2*r2; /* r6 */
182 else if (sc_r_power == fourtyeight)
184 rp = r2*r2*r2; /* r6 */
185 rp = rp*rp; /* r12 */
186 rp = rp*rp; /* r24 */
187 rp = rp*rp; /* r48 */
188 rpm2 = rp/r2; /* r46 */
192 rp = std::pow(r2, half * sc_r_power); /* not currently supported as input, but can handle it */
196 /* Loop over state A(0) and B(1) */
197 for (i = 0; i < 2; i++)
199 if ((c6[i] > 0) && (c12[i] > 0))
201 /* The c6 & c12 coefficients now contain the constants 6.0 and 12.0, respectively.
202 * Correct for this by multiplying with (1/12.0)/(1/6.0)=6.0/12.0=0.5.
204 sigma6[i] = half*c12[i]/c6[i];
205 sigma2[i] = std::cbrt(half*c12[i]/c6[i]);
206 /* should be able to get rid of this ^^^ internal pow call eventually. Will require agreement on
207 what data to store externally. Can't be fixed without larger scale changes, so not 5.0 */
208 if (sigma6[i] < sigma6_min) /* for disappearing coul and vdw with soft core at the same time */
210 sigma6[i] = sigma6_min;
211 sigma2[i] = sigma2_min;
216 sigma6[i] = sigma6_def;
217 sigma2[i] = sigma2_def;
219 if (sc_r_power == six)
221 sigma_pow[i] = sigma6[i];
223 else if (sc_r_power == fourtyeight)
225 sigma_pow[i] = sigma6[i]*sigma6[i]; /* sigma^12 */
226 sigma_pow[i] = sigma_pow[i]*sigma_pow[i]; /* sigma^24 */
227 sigma_pow[i] = sigma_pow[i]*sigma_pow[i]; /* sigma^48 */
230 { /* not really supported as input, but in here for testing the general case*/
231 sigma_pow[i] = std::pow(sigma2[i], sc_r_power/2);
235 /* only use softcore if one of the states has a zero endstate - softcore is for avoiding infinities!*/
236 if ((c12[0] > 0) && (c12[1] > 0))
243 alpha_vdw_eff = alpha_vdw;
244 alpha_coul_eff = alpha_coul;
247 /* Loop over A and B states again */
248 for (i = 0; i < 2; i++)
255 /* Only spend time on A or B state if it is non-zero */
256 if ( (qq[i] != 0) || (c6[i] != 0) || (c12[i] != 0) )
259 rpinv = one/(alpha_coul_eff*lfac_coul[i]*sigma_pow[i]+rp);
260 r_coul = std::pow(rpinv, minusOne / sc_r_power);
262 /* Electrostatics table lookup data */
263 rtab = r_coul*tabscale;
264 ntab = static_cast<int>(rtab);
267 ntab = static_cast<int>(tableStride*ntab);
271 Geps = eps*vftab[ntab+2];
272 Heps2 = eps2*vftab[ntab+3];
275 FF = Fp+Geps+two*Heps2;
277 fscal_elec[i] = -qq[i]*FF*r_coul*rpinv*tabscale;
280 rpinv = one/(alpha_vdw_eff*lfac_vdw[i]*sigma_pow[i]+rp);
281 r_vdw = std::pow(rpinv, minusOne / sc_r_power);
282 /* Vdw table lookup data */
283 rtab = r_vdw*tabscale;
284 ntab = static_cast<int>(rtab);
291 Geps = eps*vftab[ntab+6];
292 Heps2 = eps2*vftab[ntab+7];
295 FF = Fp+Geps+two*Heps2;
297 fscal_vdw[i] = -c6[i]*FF;
302 Geps = eps*vftab[ntab+10];
303 Heps2 = eps2*vftab[ntab+11];
306 FF = Fp+Geps+two*Heps2;
307 vvdw[i] += c12[i]*VV;
308 fscal_vdw[i] -= c12[i]*FF;
309 fscal_vdw[i] *= r_vdw*rpinv*tabscale;
312 /* Now we have velec[i], vvdw[i], and fscal[i] for both states */
313 /* Assemble A and B states */
319 for (i = 0; i < 2; i++)
321 velecsum += LFC[i]*velec[i];
322 vvdwsum += LFV[i]*vvdw[i];
324 fscal += (LFC[i]*fscal_elec[i]+LFV[i]*fscal_vdw[i])*rpm2;
326 dvdl_coul += velec[i]*DLF[i] + LFC[i]*alpha_coul_eff*dlfac_coul[i]*fscal_elec[i]*sigma_pow[i];
327 dvdl_vdw += vvdw[i]*DLF[i] + LFV[i]*alpha_vdw_eff*dlfac_vdw[i]*fscal_vdw[i]*sigma_pow[i];
330 dvdl[efptCOUL] += dvdl_coul;
331 dvdl[efptVDW] += dvdl_vdw;
333 *velectot = velecsum;
339 /*! \brief Calculate pair interactions, supports all types and conditions. */
341 do_pairs_general(int ftype, int nbonds,
342 const t_iatom iatoms[], const t_iparams iparams[],
343 const rvec x[], rvec4 f[], rvec fshift[],
344 const struct t_pbc *pbc, const struct t_graph *g,
345 const real *lambda, real *dvdl,
347 const t_forcerec *fr, gmx_grppairener_t *grppener,
348 int *global_atom_index)
353 int i, itype, ai, aj, gid;
356 real fscal, velec, vvdw;
357 real * energygrp_elec;
358 real * energygrp_vdw;
359 static gmx_bool warned_rlimit = FALSE;
360 /* Free energy stuff */
361 gmx_bool bFreeEnergy;
362 real LFC[2], LFV[2], DLF[2], lfac_coul[2], lfac_vdw[2], dlfac_coul[2], dlfac_vdw[2];
363 real qqB, c6B, c12B, sigma2_def, sigma2_min;
369 energygrp_elec = grppener->ener[egCOUL14];
370 energygrp_vdw = grppener->ener[egLJ14];
373 energygrp_elec = grppener->ener[egCOULSR];
374 energygrp_vdw = grppener->ener[egLJSR];
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 */
395 sigma2_def = std::cbrt(fr->sc_sigma6_def);
396 sigma2_min = std::cbrt(fr->sc_sigma6_min);
398 for (i = 0; i < 2; i++)
400 lfac_coul[i] = (fr->sc_power == 2 ? (1-LFC[i])*(1-LFC[i]) : (1-LFC[i]));
401 dlfac_coul[i] = DLF[i]*fr->sc_power/fr->sc_r_power*(fr->sc_power == 2 ? (1-LFC[i]) : 1);
402 lfac_vdw[i] = (fr->sc_power == 2 ? (1-LFV[i])*(1-LFV[i]) : (1-LFV[i]));
403 dlfac_vdw[i] = DLF[i]*fr->sc_power/fr->sc_r_power*(fr->sc_power == 2 ? (1-LFV[i]) : 1);
408 sigma2_min = sigma2_def = 0;
411 /* TODO This code depends on the logic in tables.c that constructs
412 the table layout, which should be made explicit in future
414 GMX_ASSERT(etiNR == 3, "Pair-interaction code that uses GROMACS interaction tables supports exactly 3 tables");
415 GMX_ASSERT(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);
433 (fr->efep != efepNO &&
434 (((md->nPerturbed != 0) && (md->bPerturbed[ai] || md->bPerturbed[aj])) ||
435 iparams[itype].lj14.c6A != iparams[itype].lj14.c6B ||
436 iparams[itype].lj14.c12A != iparams[itype].lj14.c12B));
437 qq = md->chargeA[ai]*md->chargeA[aj]*epsfac*fr->fudgeQQ;
438 c6 = iparams[itype].lj14.c6A;
439 c12 = iparams[itype].lj14.c12A;
442 qq = iparams[itype].ljc14.qi*iparams[itype].ljc14.qj*epsfac*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(r2, fr->sc_r_power, fr->sc_alphacoul, fr->sc_alphavdw,
497 fr->pairsTable->scale, fr->pairsTable->data, fr->pairsTable->stride,
498 qq, c6, c12, qqB, c6B, c12B,
499 LFC, LFV, DLF, lfac_coul, lfac_vdw, dlfac_coul, dlfac_vdw,
500 fr->sc_sigma6_def, fr->sc_sigma6_min, sigma2_def, sigma2_min, &velec, &vvdw, dvdl);
504 /* Evaluate tabulated interaction without free energy */
505 fscal = evaluate_single(r2, fr->pairsTable->scale, fr->pairsTable->data, fr->pairsTable->stride,
506 qq, c6, c12, &velec, &vvdw);
509 energygrp_elec[gid] += velec;
510 energygrp_vdw[gid] += vvdw;
511 svmul(fscal, dx, dx);
519 /* Correct the shift forces using the graph */
520 ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, aj), dt);
521 fshift_index = IVEC2IS(dt);
523 if (fshift_index != CENTRAL)
525 rvec_inc(fshift[fshift_index], dx);
526 rvec_dec(fshift[CENTRAL], dx);
532 /*! \brief Calculate pairs, only for plain-LJ + plain Coulomb normal type.
534 * This function is templated for real/SimdReal and for optimization.
536 template<typename T, int pack_size,
539 do_pairs_simple(int nbonds,
540 const t_iatom iatoms[], const t_iparams iparams[],
541 const rvec x[], rvec4 f[],
544 const real scale_factor)
546 const int nfa1 = 1 + 2;
552 #if GMX_SIMD_HAVE_REAL
553 alignas(GMX_SIMD_ALIGNMENT) std::int32_t ai[pack_size];
554 alignas(GMX_SIMD_ALIGNMENT) std::int32_t aj[pack_size];
555 alignas(GMX_SIMD_ALIGNMENT) real coeff[3*pack_size];
557 std::int32_t ai[pack_size];
558 std::int32_t aj[pack_size];
559 real coeff[3*pack_size];
562 /* nbonds is #pairs*nfa1, here we step pack_size pairs */
563 for (int i = 0; i < nbonds; i += pack_size*nfa1)
565 /* Collect atoms for pack_size pairs.
566 * iu indexes into iatoms, we should not let iu go beyond nbonds.
569 for (int s = 0; s < pack_size; s++)
571 int itype = iatoms[iu];
572 ai[s] = iatoms[iu + 1];
573 aj[s] = iatoms[iu + 2];
575 if (i + s*nfa1 < nbonds)
577 coeff[0*pack_size + s] = iparams[itype].lj14.c6A;
578 coeff[1*pack_size + s] = iparams[itype].lj14.c12A;
579 coeff[2*pack_size + s] = md->chargeA[ai[s]]*md->chargeA[aj[s]];
581 /* Avoid indexing the iatoms array out of bounds.
582 * We pad the coordinate indices with the last atom pair.
584 if (iu + nfa1 < nbonds)
591 /* Pad the coefficient arrays with zeros to get zero forces */
592 coeff[0*pack_size + s] = 0;
593 coeff[1*pack_size + s] = 0;
594 coeff[2*pack_size + s] = 0;
598 /* Load the coordinates */
600 gatherLoadUTranspose<3>(reinterpret_cast<const real *>(x), ai, &xi[XX], &xi[YY], &xi[ZZ]);
601 gatherLoadUTranspose<3>(reinterpret_cast<const real *>(x), aj, &xj[XX], &xj[YY], &xj[ZZ]);
603 T c6 = load<T>(coeff + 0*pack_size);
604 T c12 = load<T>(coeff + 1*pack_size);
605 T qq = load<T>(coeff + 2*pack_size);
607 /* We could save these operations by storing 6*C6,12*C12 */
612 pbc_dx_aiuc(pbc, xi, xj, dr);
614 T rsq = dr[XX]*dr[XX] + dr[YY]*dr[YY] + dr[ZZ]*dr[ZZ];
615 T rinv = gmx::invsqrt(rsq);
617 T rinv6 = rinv2*rinv2*rinv2;
619 /* Calculate the Coulomb force * r */
622 /* Calculate the LJ force * r and add it to the Coulomb part */
623 T fr = gmx::fma(fms(c12, rinv6, c6), rinv6, cfr);
630 /* Add the pair forces to the force array.
631 * Note that here we might add multiple force components for some atoms
632 * due to the SIMD padding. But the extra force components are zero.
634 transposeScatterIncrU<4>(reinterpret_cast<real *>(f), ai, fx, fy, fz);
635 transposeScatterDecrU<4>(reinterpret_cast<real *>(f), aj, fx, fy, fz);
639 /*! \brief Calculate all listed pair interactions */
641 do_pairs(int ftype, int nbonds,
642 const t_iatom iatoms[], const t_iparams iparams[],
643 const rvec x[], rvec4 f[], rvec fshift[],
644 const struct t_pbc *pbc, const struct t_graph *g,
645 const real *lambda, real *dvdl,
647 const t_forcerec *fr,
648 const gmx_bool computeForcesOnly,
649 gmx_grppairener_t *grppener,
650 int *global_atom_index)
652 if (ftype == F_LJ14 &&
653 fr->ic->vdwtype != evdwUSER && !EEL_USER(fr->ic->eeltype) &&
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.
665 alignas(GMX_SIMD_ALIGNMENT) real pbc_simd[9*GMX_SIMD_REAL_WIDTH];
666 set_pbc_simd(pbc, pbc_simd);
668 do_pairs_simple<SimdReal, GMX_SIMD_REAL_WIDTH,
669 const real *>(nbonds, iatoms, iparams,
671 md, fr->ic->epsfac*fr->fudgeQQ);
673 /* This construct is needed because pbc_dx_aiuc doesn't accept pbc=NULL */
675 const t_pbc *pbc_nonnull;
683 set_pbc(&pbc_no, epbcNONE, nullptr);
684 pbc_nonnull = &pbc_no;
687 do_pairs_simple<real, 1,
688 const t_pbc *>(nbonds, iatoms, iparams,
690 md, fr->ic->epsfac*fr->fudgeQQ);
695 do_pairs_general(ftype, nbonds, iatoms, iparams,
696 x, f, fshift, pbc, g,
698 md, fr, grppener, global_atom_index);