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 template<KernelSoftcoreType softcoreType>
173 static real free_energy_evaluate_single(real r2,
175 const interaction_const_t::SoftCoreParameters& scParams,
189 const real lfac_coul[2],
190 const real lfac_vdw[2],
191 const real dlfac_coul[2],
192 const real dlfac_vdw[2],
197 real rtab, eps, eps2, Y, F, Geps, Heps2, Fp, VV, FF, fscal;
198 real qq[2], c6[2], c12[2], sigma6[2], sigma_pow[2];
199 real alpha_coul_eff, alpha_vdw_eff, dvdl_coul_sum, dvdl_vdw_sum;
200 real rpinv, r_coul, r_vdw, velecsum, vvdwsum;
201 real fscal_vdw[2], fscal_elec[2];
202 real velec[2], vvdw[2];
203 real dvdl_elec[2], dvdl_vdw[2];
207 const real half = 0.5_real;
208 const real one = 1.0_real;
209 const real two = 2.0_real;
218 const real rpm2 = r2 * r2; /* r4 */
219 const real rp = rpm2 * r2; /* r6 */
220 const real r = sqrt(r2); /* r1 */
222 /* Loop over state A(0) and B(1) */
223 for (i = 0; i < 2; i++)
225 if (softcoreType == KernelSoftcoreType::Beutler || softcoreType == KernelSoftcoreType::Gapsys)
227 if ((c6[i] > 0) && (c12[i] > 0))
229 /* The c6 & c12 coefficients now contain the constants 6.0 and 12.0, respectively.
230 * Correct for this by multiplying with (1/12.0)/(1/6.0)=6.0/12.0=0.5.
232 sigma6[i] = half * c12[i] / c6[i];
233 if (softcoreType == KernelSoftcoreType::Beutler && (sigma6[i] < scParams.sigma6Minimum)) /* for disappearing coul and vdw with soft core at the same time */
235 sigma6[i] = scParams.sigma6Minimum;
240 sigma6[i] = scParams.sigma6WithInvalidSigma;
242 sigma_pow[i] = sigma6[i];
246 if (softcoreType == KernelSoftcoreType::Beutler || softcoreType == KernelSoftcoreType::Gapsys)
248 /* only use softcore if one of the states has a zero endstate - softcore is for avoiding infinities!*/
249 if ((c12[0] > 0) && (c12[1] > 0))
256 if (softcoreType == KernelSoftcoreType::Beutler)
258 alpha_vdw_eff = scParams.alphaVdw;
259 alpha_coul_eff = scParams.alphaCoulomb;
261 else if (softcoreType == KernelSoftcoreType::Gapsys)
263 alpha_vdw_eff = scParams.alphaVdw;
264 alpha_coul_eff = gmx::sixthroot(scParams.sigma6WithInvalidSigma);
269 /* Loop over A and B states again */
270 for (i = 0; i < 2; i++)
281 /* Only spend time on A or B state if it is non-zero */
282 if ((qq[i] != 0) || (c6[i] != 0) || (c12[i] != 0))
285 if (softcoreType == KernelSoftcoreType::Beutler)
287 rpinv = one / (alpha_coul_eff * lfac_coul[i] * sigma_pow[i] + rp);
288 r_coul = sixthRoot(rpinv);
296 if (softcoreType == KernelSoftcoreType::Gapsys)
298 rQ = gmx::sixthroot(one - LFC[i]) * (one + std::fabs(qq[i] / facel));
299 rQ *= alpha_coul_eff;
301 if (rQ > rCoulCutoff)
308 if ((softcoreType == KernelSoftcoreType::Gapsys) && (r < rQ))
310 real rInvQ = one / rQ;
311 real constFac = qq[i] * rInvQ;
312 real linFac = constFac * r * rInvQ;
313 real quadrFac = linFac * r * rInvQ;
315 /* Computing Coulomb force and potential energy */
316 fscal_elec[i] = 2 * quadrFac - 3 * linFac;
317 fscal_elec[i] *= rpinv;
319 velec[i] = quadrFac - 3 * (linFac - constFac);
321 dvdl_elec[i] += scaleDvdlRCoul * DLF[i] * half * (LFC[i] / (1 - LFC[i]))
322 * (quadrFac - 2 * linFac + constFac);
324 else // Beutler, resp. hardcore
326 /* Electrostatics table lookup data */
327 rtab = r_coul * tabscale;
328 ntab = static_cast<int>(rtab);
331 ntab = static_cast<int>(tableStride * ntab);
335 Geps = eps * vftab[ntab + 2];
336 Heps2 = eps2 * vftab[ntab + 3];
337 Fp = F + Geps + Heps2;
339 FF = Fp + Geps + two * Heps2;
340 velec[i] = qq[i] * VV;
341 fscal_elec[i] = -qq[i] * FF * r_coul * rpinv * tabscale;
345 if (softcoreType == KernelSoftcoreType::Beutler)
347 rpinv = one / (alpha_vdw_eff * lfac_vdw[i] * sigma_pow[i] + rp);
348 r_vdw = sixthRoot(rpinv);
356 if (softcoreType == KernelSoftcoreType::Gapsys)
358 constexpr real c_twentySixSeventh = 26.0_real / 7.0_real;
360 rLJ = gmx::sixthroot(c_twentySixSeventh * sigma6[i] * (one - LFV[i]));
361 rLJ *= alpha_vdw_eff;
364 if ((softcoreType == KernelSoftcoreType::Gapsys) && (r < rLJ))
366 // scaled values for c6 and c12
368 c6s = c6[i] / 6.0_real;
369 c12s = c12[i] / 12.0_real;
371 /* Temporary variables for inverted values */
372 real rInvLJ = one / rLJ;
373 real rInv14, rInv13, rInv12;
374 real rInv8, rInv7, rInv6;
375 rInv6 = rInvLJ * rInvLJ * rInvLJ;
377 rInv7 = rInv6 * rInvLJ;
378 rInv8 = rInv7 * rInvLJ;
379 rInv14 = c12s * rInv7 * rInv7 * r2;
380 rInv13 = c12s * rInv7 * rInv6 * r;
381 rInv12 = c12s * rInv6 * rInv6;
386 /* Temporary variables for A and B */
387 real quadrFac, linearFac, constFac;
388 quadrFac = 156 * rInv14 - 42 * rInv8;
389 linearFac = 168 * rInv13 - 48 * rInv7;
390 constFac = 91 * rInv12 - 28 * rInv6;
392 /* Computing LJ force and potential energy*/
393 fscal_vdw[i] = quadrFac - linearFac;
394 fscal_vdw[i] *= rpinv;
396 vvdw[i] = 0.5 * quadrFac - linearFac + constFac;
398 dvdl_vdw[i] += DLF[i] * 28 * (LFV[i] / (one - LFV[i]))
399 * ((6.5_real * rInv14 - rInv8) - (13 * rInv13 - 2 * rInv7)
400 + (6.5_real * rInv12 - rInv6));
402 else // Beutler, resp. hardcore
404 /* Vdw table lookup data */
405 rtab = r_vdw * tabscale;
406 ntab = static_cast<int>(rtab);
413 Geps = eps * vftab[ntab + 6];
414 Heps2 = eps2 * vftab[ntab + 7];
415 Fp = F + Geps + Heps2;
417 FF = Fp + Geps + two * Heps2;
418 vvdw[i] = c6[i] * VV;
419 fscal_vdw[i] = -c6[i] * FF;
424 Geps = eps * vftab[ntab + 10];
425 Heps2 = eps2 * vftab[ntab + 11];
426 Fp = F + Geps + Heps2;
428 FF = Fp + Geps + two * Heps2;
429 vvdw[i] += c12[i] * VV;
430 fscal_vdw[i] -= c12[i] * FF;
431 fscal_vdw[i] *= r_vdw * rpinv * tabscale;
435 /* Now we have velec[i], vvdw[i], and fscal[i] for both states */
436 /* Assemble A and B states */
442 for (i = 0; i < 2; i++)
444 velecsum += LFC[i] * velec[i];
445 vvdwsum += LFV[i] * vvdw[i];
447 fscal += (LFC[i] * fscal_elec[i] + LFV[i] * fscal_vdw[i]) * rpm2;
449 if (softcoreType == KernelSoftcoreType::Gapsys)
451 dvdl_coul_sum += dvdl_elec[i];
452 dvdl_vdw_sum += dvdl_vdw[i];
454 dvdl_coul_sum += velec[i] * DLF[i];
455 dvdl_vdw_sum += vvdw[i] * DLF[i];
456 if (softcoreType == KernelSoftcoreType::Beutler)
458 dvdl_coul_sum += LFC[i] * alpha_coul_eff * dlfac_coul[i] * fscal_elec[i] * sigma_pow[i];
459 dvdl_vdw_sum += LFV[i] * alpha_vdw_eff * dlfac_vdw[i] * fscal_vdw[i] * sigma_pow[i];
463 dvdl[static_cast<int>(FreeEnergyPerturbationCouplingType::Coul)] += dvdl_coul_sum;
464 dvdl[static_cast<int>(FreeEnergyPerturbationCouplingType::Vdw)] += dvdl_vdw_sum;
466 *velectot = velecsum;
472 /*! \brief Calculate pair interactions, supports all types and conditions. */
473 template<BondedKernelFlavor flavor>
474 static real do_pairs_general(int ftype,
476 const t_iatom iatoms[],
477 const t_iparams iparams[],
481 const struct t_pbc* pbc,
484 gmx::ArrayRef<real> chargeA,
485 gmx::ArrayRef<real> chargeB,
486 gmx::ArrayRef<bool> atomIsPerturbed,
487 gmx::ArrayRef<unsigned short> cENER,
489 const t_forcerec* fr,
490 gmx_grppairener_t* grppener,
491 int* global_atom_index)
495 int i, itype, ai, aj, gid;
498 real fscal, velec, vvdw;
499 real* energygrp_elec;
501 static gmx_bool warned_rlimit = FALSE;
502 /* Free energy stuff */
503 gmx_bool bFreeEnergy;
504 real LFC[2], LFV[2], DLF[2], lfac_coul[2], lfac_vdw[2], dlfac_coul[2], dlfac_vdw[2];
506 const real oneSixth = 1.0_real / 6.0_real;
512 energygrp_elec = grppener->energyGroupPairTerms[NonBondedEnergyTerms::Coulomb14].data();
513 energygrp_vdw = grppener->energyGroupPairTerms[NonBondedEnergyTerms::LJ14].data();
516 energygrp_elec = grppener->energyGroupPairTerms[NonBondedEnergyTerms::CoulombSR].data();
517 energygrp_vdw = grppener->energyGroupPairTerms[NonBondedEnergyTerms::LJSR].data();
520 energygrp_elec = nullptr; /* Keep compiler happy */
521 energygrp_vdw = nullptr; /* Keep compiler happy */
522 gmx_fatal(FARGS, "Unknown function type %d in do_nonbonded14", ftype);
525 if (fr->efep != FreeEnergyPerturbationType::No)
527 /* Lambda factor for state A=1-lambda and B=lambda */
528 LFC[0] = 1.0 - lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Coul)];
529 LFV[0] = 1.0 - lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Vdw)];
530 LFC[1] = lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Coul)];
531 LFV[1] = lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Vdw)];
533 /*derivative of the lambda factor for state A and B */
537 GMX_ASSERT(fr->ic->softCoreParameters, "We need soft-core parameters");
538 const auto& scParams = *fr->ic->softCoreParameters;
540 for (i = 0; i < 2; i++)
542 lfac_coul[i] = (scParams.lambdaPower == 2 ? (1 - LFC[i]) * (1 - LFC[i]) : (1 - LFC[i]));
543 dlfac_coul[i] = DLF[i] * scParams.lambdaPower * oneSixth
544 * (scParams.lambdaPower == 2 ? (1 - LFC[i]) : 1);
545 lfac_vdw[i] = (scParams.lambdaPower == 2 ? (1 - LFV[i]) * (1 - LFV[i]) : (1 - LFV[i]));
546 dlfac_vdw[i] = DLF[i] * scParams.lambdaPower * oneSixth
547 * (scParams.lambdaPower == 2 ? (1 - LFV[i]) : 1);
551 /* TODO This code depends on the logic in tables.c that constructs
552 the table layout, which should be made explicit in future
556 "Pair-interaction code that uses GROMACS interaction tables supports exactly 3 tables");
558 fr->pairsTable->interaction == TableInteraction::ElectrostaticVdwRepulsionVdwDispersion,
559 "Pair interaction kernels need a table with Coulomb, repulsion and dispersion entries");
561 const real epsfac = fr->ic->epsfac;
564 for (i = 0; (i < nbonds);)
569 gid = GID(cENER[ai], cENER[aj], numEnergyGroups);
576 (fr->efep != FreeEnergyPerturbationType::No
577 && ((!atomIsPerturbed.empty() && (atomIsPerturbed[ai] || atomIsPerturbed[aj]))
578 || iparams[itype].lj14.c6A != iparams[itype].lj14.c6B
579 || iparams[itype].lj14.c12A != iparams[itype].lj14.c12B));
580 qq = chargeA[ai] * chargeA[aj] * epsfac * fr->fudgeQQ;
581 c6 = iparams[itype].lj14.c6A;
582 c12 = iparams[itype].lj14.c12A;
585 qq = iparams[itype].ljc14.qi * iparams[itype].ljc14.qj * epsfac
586 * iparams[itype].ljc14.fqq;
587 c6 = iparams[itype].ljc14.c6;
588 c12 = iparams[itype].ljc14.c12;
591 qq = iparams[itype].ljcnb.qi * iparams[itype].ljcnb.qj * epsfac;
592 c6 = iparams[itype].ljcnb.c6;
593 c12 = iparams[itype].ljcnb.c12;
596 /* Cannot happen since we called gmx_fatal() above in this case */
597 qq = c6 = c12 = 0; /* Keep compiler happy */
601 /* To save flops in the optimized kernels, c6/c12 have 6.0/12.0 derivative prefactors
602 * included in the general nfbp array now. This means the tables are scaled down by the
603 * same factor, so when we use the original c6/c12 parameters from iparams[] they must
609 /* Do we need to apply full periodic boundary conditions? */
612 fshift_index = pbc_dx_aiuc(pbc, x[ai], x[aj], dx);
616 fshift_index = c_centralShiftIndex;
617 rvec_sub(x[ai], x[aj], dx);
621 if (r2 >= fr->pairsTable->interactionRange * fr->pairsTable->interactionRange)
623 /* This check isn't race free. But it doesn't matter because if a race occurs the only
624 * disadvantage is that the warning is printed twice */
627 warning_rlimit(x, ai, aj, global_atom_index, sqrt(r2), fr->pairsTable->interactionRange);
628 warned_rlimit = TRUE;
635 /* Currently free energy is only supported for F_LJ14, so no need to check for that if we got here */
636 qqB = chargeB[ai] * chargeB[aj] * epsfac * fr->fudgeQQ;
637 c6B = iparams[itype].lj14.c6B * 6.0;
638 c12B = iparams[itype].lj14.c12B * 12.0;
640 const auto& scParams = *fr->ic->softCoreParameters;
641 if (scParams.softcoreType == SoftcoreType::Gapsys)
643 fscal = free_energy_evaluate_single<KernelSoftcoreType::Gapsys>(
646 *fr->ic->softCoreParameters,
647 fr->pairsTable->scale,
648 fr->pairsTable->data.data(),
649 fr->pairsTable->stride,
668 else if (scParams.softcoreType == SoftcoreType::Beutler)
670 fscal = free_energy_evaluate_single<KernelSoftcoreType::Beutler>(
673 *fr->ic->softCoreParameters,
674 fr->pairsTable->scale,
675 fr->pairsTable->data.data(),
676 fr->pairsTable->stride,
697 fscal = free_energy_evaluate_single<KernelSoftcoreType::None>(
700 *fr->ic->softCoreParameters,
701 fr->pairsTable->scale,
702 fr->pairsTable->data.data(),
703 fr->pairsTable->stride,
725 /* Evaluate tabulated interaction without free energy */
726 fscal = evaluate_single(r2,
727 fr->pairsTable->scale,
728 fr->pairsTable->data.data(),
729 fr->pairsTable->stride,
737 energygrp_elec[gid] += velec;
738 energygrp_vdw[gid] += vvdw;
739 svmul(fscal, dx, dx);
745 if (computeVirial(flavor))
747 if (fshift_index != c_centralShiftIndex)
749 rvec_inc(fshift[fshift_index], dx);
750 rvec_dec(fshift[c_centralShiftIndex], dx);
757 /*! \brief Calculate pairs, only for plain-LJ + plain Coulomb normal type.
759 * This function is templated for real/SimdReal and for optimization.
761 template<typename T, int pack_size, typename pbc_type>
762 static void do_pairs_simple(int nbonds,
763 const t_iatom iatoms[],
764 const t_iparams iparams[],
768 gmx::ArrayRef<real> charge,
769 const real scale_factor)
771 const int nfa1 = 1 + 2;
777 #if GMX_SIMD_HAVE_REAL
778 alignas(GMX_SIMD_ALIGNMENT) std::int32_t ai[pack_size];
779 alignas(GMX_SIMD_ALIGNMENT) std::int32_t aj[pack_size];
780 alignas(GMX_SIMD_ALIGNMENT) real coeff[3 * pack_size];
782 std::int32_t ai[pack_size];
783 std::int32_t aj[pack_size];
784 real coeff[3 * pack_size];
787 /* nbonds is #pairs*nfa1, here we step pack_size pairs */
788 for (int i = 0; i < nbonds; i += pack_size * nfa1)
790 /* Collect atoms for pack_size pairs.
791 * iu indexes into iatoms, we should not let iu go beyond nbonds.
794 for (int s = 0; s < pack_size; s++)
796 int itype = iatoms[iu];
797 ai[s] = iatoms[iu + 1];
798 aj[s] = iatoms[iu + 2];
800 if (i + s * nfa1 < nbonds)
802 coeff[0 * pack_size + s] = iparams[itype].lj14.c6A;
803 coeff[1 * pack_size + s] = iparams[itype].lj14.c12A;
804 coeff[2 * pack_size + s] = charge[ai[s]] * charge[aj[s]];
806 /* Avoid indexing the iatoms array out of bounds.
807 * We pad the coordinate indices with the last atom pair.
809 if (iu + nfa1 < nbonds)
816 /* Pad the coefficient arrays with zeros to get zero forces */
817 coeff[0 * pack_size + s] = 0;
818 coeff[1 * pack_size + s] = 0;
819 coeff[2 * pack_size + s] = 0;
823 /* Load the coordinates */
825 gatherLoadUTranspose<3>(reinterpret_cast<const real*>(x), ai, &xi[XX], &xi[YY], &xi[ZZ]);
826 gatherLoadUTranspose<3>(reinterpret_cast<const real*>(x), aj, &xj[XX], &xj[YY], &xj[ZZ]);
828 T c6 = load<T>(coeff + 0 * pack_size);
829 T c12 = load<T>(coeff + 1 * pack_size);
830 T qq = load<T>(coeff + 2 * pack_size);
832 /* We could save these operations by storing 6*C6,12*C12 */
837 pbc_dx_aiuc(pbc, xi, xj, dr);
839 T rsq = dr[XX] * dr[XX] + dr[YY] * dr[YY] + dr[ZZ] * dr[ZZ];
840 T rinv = gmx::invsqrt(rsq);
841 T rinv2 = rinv * rinv;
842 T rinv6 = rinv2 * rinv2 * rinv2;
844 /* Calculate the Coulomb force * r */
845 T cfr = ef * qq * rinv;
847 /* Calculate the LJ force * r and add it to the Coulomb part */
848 T fr = gmx::fma(fms(c12, rinv6, c6), rinv6, cfr);
850 T finvr = fr * rinv2;
851 T fx = finvr * dr[XX];
852 T fy = finvr * dr[YY];
853 T fz = finvr * dr[ZZ];
855 /* Add the pair forces to the force array.
856 * Note that here we might add multiple force components for some atoms
857 * due to the SIMD padding. But the extra force components are zero.
859 transposeScatterIncrU<4>(reinterpret_cast<real*>(f), ai, fx, fy, fz);
860 transposeScatterDecrU<4>(reinterpret_cast<real*>(f), aj, fx, fy, fz);
864 /*! \brief Calculate all listed pair interactions */
865 void do_pairs(int ftype,
867 const t_iatom iatoms[],
868 const t_iparams iparams[],
872 const struct t_pbc* pbc,
875 gmx::ArrayRef<real> chargeA,
876 gmx::ArrayRef<real> chargeB,
877 gmx::ArrayRef<bool> atomIsPerturbed,
878 gmx::ArrayRef<unsigned short> cENER,
879 const int numEnergyGroups,
880 const t_forcerec* fr,
881 const bool havePerturbedInteractions,
882 const gmx::StepWorkload& stepWork,
883 gmx_grppairener_t* grppener,
884 int* global_atom_index)
886 if (ftype == F_LJ14 && fr->ic->vdwtype != VanDerWaalsType::User && !EEL_USER(fr->ic->eeltype)
887 && !havePerturbedInteractions && (!stepWork.computeVirial && !stepWork.computeEnergy))
889 /* We use a fast code-path for plain LJ 1-4 without FEP.
891 * TODO: Add support for energies (straightforward) and virial
892 * in the SIMD template. For the virial it's inconvenient to store
893 * the force sums for the shifts and we should directly calculate
894 * and sum the virial for the shifts. But we should do this
895 * at once for the angles and dihedrals as well.
897 #if GMX_SIMD_HAVE_REAL
898 if (fr->use_simd_kernels)
900 alignas(GMX_SIMD_ALIGNMENT) real pbc_simd[9 * GMX_SIMD_REAL_WIDTH];
901 set_pbc_simd(pbc, pbc_simd);
903 do_pairs_simple<SimdReal, GMX_SIMD_REAL_WIDTH, const real*>(
904 nbonds, iatoms, iparams, x, f, pbc_simd, chargeA, fr->ic->epsfac * fr->fudgeQQ);
909 /* This construct is needed because pbc_dx_aiuc doesn't accept pbc=NULL */
911 const t_pbc* pbc_nonnull;
919 set_pbc(&pbc_no, PbcType::No, nullptr);
920 pbc_nonnull = &pbc_no;
923 do_pairs_simple<real, 1, const t_pbc*>(
924 nbonds, iatoms, iparams, x, f, pbc_nonnull, chargeA, fr->ic->epsfac * fr->fudgeQQ);
927 else if (stepWork.computeVirial)
929 do_pairs_general<BondedKernelFlavor::ForcesAndVirialAndEnergy>(ftype,
950 do_pairs_general<BondedKernelFlavor::ForcesAndEnergy>(ftype,