* This file is part of the GROMACS molecular simulation package.
*
* Copyright (c) 2014,2015,2016,2017,2018 by the GROMACS development team.
- * Copyright (c) 2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2019,2020,2021, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
#include "gromacs/listed_forces/bonded.h"
#include "gromacs/math/functions.h"
#include "gromacs/math/vec.h"
+#include "gromacs/mdtypes/forcerec.h"
#include "gromacs/mdtypes/group.h"
+#include "gromacs/mdtypes/interaction_const.h"
#include "gromacs/mdtypes/md_enums.h"
#include "gromacs/mdtypes/nblist.h"
#include "gromacs/mdtypes/simulation_workload.h"
#include "gromacs/pbcutil/ishift.h"
-#include "gromacs/pbcutil/mshift.h"
#include "gromacs/pbcutil/pbc.h"
#include "gromacs/pbcutil/pbc_simd.h"
#include "gromacs/simd/simd.h"
"IMPORTANT: This should not happen in a stable simulation, so there is\n"
"probably something wrong with your system. Only change the table-extension\n"
"distance in the mdp file if you are really sure that is the reason.\n",
- glatnr(global_atom_index, ai), glatnr(global_atom_index, aj), r, rlimit);
+ glatnr(global_atom_index, ai),
+ glatnr(global_atom_index, aj),
+ r,
+ rlimit);
if (debug)
{
fprintf(debug,
"%8f %8f %8f\n%8f %8f %8f\n1-4 (%d,%d) interaction not within cut-off! r=%g. "
"Ignored\n",
- x[ai][XX], x[ai][YY], x[ai][ZZ], x[aj][XX], x[aj][YY], x[aj][ZZ],
- glatnr(global_atom_index, ai), glatnr(global_atom_index, aj), r);
+ x[ai][XX],
+ x[ai][YY],
+ x[ai][ZZ],
+ x[aj][XX],
+ x[aj][YY],
+ x[aj][ZZ],
+ glatnr(global_atom_index, ai),
+ glatnr(global_atom_index, aj),
+ r);
}
}
return fscal;
}
+static inline real sixthRoot(const real r)
+{
+ return gmx::invsqrt(std::cbrt(r));
+}
+
/*! \brief Compute the energy and force for a single pair interaction under FEP */
-static real free_energy_evaluate_single(real r2,
- real sc_r_power,
- real alpha_coul,
- real alpha_vdw,
- real tabscale,
- const real* vftab,
- real tableStride,
- real qqA,
- real c6A,
- real c12A,
- real qqB,
- real c6B,
- real c12B,
- const real LFC[2],
- const real LFV[2],
- const real DLF[2],
- const real lfac_coul[2],
- const real lfac_vdw[2],
- const real dlfac_coul[2],
- const real dlfac_vdw[2],
- real sigma6_def,
- real sigma6_min,
- real sigma2_def,
- real sigma2_min,
- real* velectot,
- real* vvdwtot,
- real* dvdl)
+template<KernelSoftcoreType softcoreType>
+static real free_energy_evaluate_single(real r2,
+ real rCoulCutoff,
+ const interaction_const_t::SoftCoreParameters& scParams,
+ real tabscale,
+ const real* vftab,
+ real tableStride,
+ real qqA,
+ real c6A,
+ real c12A,
+ real qqB,
+ real c6B,
+ real c12B,
+ real facel,
+ const real LFC[2],
+ const real LFV[2],
+ const real DLF[2],
+ const real lfac_coul[2],
+ const real lfac_vdw[2],
+ const real dlfac_coul[2],
+ const real dlfac_vdw[2],
+ real* velectot,
+ real* vvdwtot,
+ real* dvdl)
{
- real rp, rpm2, rtab, eps, eps2, Y, F, Geps, Heps2, Fp, VV, FF, fscal;
- real qq[2], c6[2], c12[2], sigma6[2], sigma2[2], sigma_pow[2];
- real alpha_coul_eff, alpha_vdw_eff, dvdl_coul, dvdl_vdw;
+ real rtab, eps, eps2, Y, F, Geps, Heps2, Fp, VV, FF, fscal;
+ real qq[2], c6[2], c12[2], sigma6[2], sigma_pow[2];
+ real alpha_coul_eff, alpha_vdw_eff, dvdl_coul_sum, dvdl_vdw_sum;
real rpinv, r_coul, r_vdw, velecsum, vvdwsum;
real fscal_vdw[2], fscal_elec[2];
real velec[2], vvdw[2];
+ real dvdl_elec[2], dvdl_vdw[2];
+ real gapsysScaleLinpointCoul, gapsysScaleLinpointVdW, gapsysSigma6VdW[2];
+ real rQ, rLJ;
+ real scaleDvdlRCoul;
int i, ntab;
- const real half = 0.5;
- const real minusOne = -1.0;
- const real one = 1.0;
- const real two = 2.0;
- const real six = 6.0;
- const real fourtyeight = 48.0;
+ const real half = 0.5_real;
+ const real one = 1.0_real;
+ const real two = 2.0_real;
qq[0] = qqA;
qq[1] = qqB;
c12[0] = c12A;
c12[1] = c12B;
- if (sc_r_power == six)
- {
- rpm2 = r2 * r2; /* r4 */
- rp = rpm2 * r2; /* r6 */
- }
- else if (sc_r_power == fourtyeight)
- {
- rp = r2 * r2 * r2; /* r6 */
- rp = rp * rp; /* r12 */
- rp = rp * rp; /* r24 */
- rp = rp * rp; /* r48 */
- rpm2 = rp / r2; /* r46 */
- }
- else
- {
- rp = std::pow(r2, half * sc_r_power); /* not currently supported as input, but can handle it */
- rpm2 = rp / r2;
- }
+ const real rpm2 = r2 * r2; /* r4 */
+ const real rp = rpm2 * r2; /* r6 */
+ const real r = sqrt(r2); /* r1 */
/* Loop over state A(0) and B(1) */
for (i = 0; i < 2; i++)
{
- if ((c6[i] > 0) && (c12[i] > 0))
+ if constexpr (softcoreType == KernelSoftcoreType::Beutler)
{
- /* The c6 & c12 coefficients now contain the constants 6.0 and 12.0, respectively.
- * Correct for this by multiplying with (1/12.0)/(1/6.0)=6.0/12.0=0.5.
- */
- sigma6[i] = half * c12[i] / c6[i];
- sigma2[i] = std::cbrt(half * c12[i] / c6[i]);
- /* should be able to get rid of this ^^^ internal pow call eventually. Will require agreement on
- what data to store externally. Can't be fixed without larger scale changes, so not 5.0 */
- if (sigma6[i] < sigma6_min) /* for disappearing coul and vdw with soft core at the same time */
+ if ((c6[i] > 0) && (c12[i] > 0))
+ {
+ /* The c6 & c12 coefficients now contain the constants 6.0 and 12.0, respectively.
+ * Correct for this by multiplying with (1/12.0)/(1/6.0)=6.0/12.0=0.5.
+ */
+ sigma6[i] = half * c12[i] / c6[i];
+ if (sigma6[i] < scParams.sigma6Minimum) /* for disappearing coul and vdw with soft core at the same time */
+ {
+ sigma6[i] = scParams.sigma6Minimum;
+ }
+ }
+ else
{
- sigma6[i] = sigma6_min;
- sigma2[i] = sigma2_min;
+ sigma6[i] = scParams.sigma6WithInvalidSigma;
}
- }
- else
- {
- sigma6[i] = sigma6_def;
- sigma2[i] = sigma2_def;
- }
- if (sc_r_power == six)
- {
sigma_pow[i] = sigma6[i];
}
- else if (sc_r_power == fourtyeight)
+ // NOLINTNEXTLINE(readability-misleading-indentation) remove when clang-tidy-13 is required
+ if constexpr (softcoreType == KernelSoftcoreType::Gapsys)
{
- sigma_pow[i] = sigma6[i] * sigma6[i]; /* sigma^12 */
- sigma_pow[i] = sigma_pow[i] * sigma_pow[i]; /* sigma^24 */
- sigma_pow[i] = sigma_pow[i] * sigma_pow[i]; /* sigma^48 */
- }
- else
- { /* not really supported as input, but in here for testing the general case*/
- sigma_pow[i] = std::pow(sigma2[i], sc_r_power / 2);
+ if ((c6[i] > 0) && (c12[i] > 0))
+ {
+ /* The c6 & c12 coefficients now contain the constants 6.0 and 12.0, respectively.
+ * Correct for this by multiplying with (1/12.0)/(1/6.0)=6.0/12.0=0.5.
+ */
+ gapsysSigma6VdW[i] = half * c12[i] / c6[i];
+ }
+ else
+ {
+ gapsysSigma6VdW[i] = scParams.gapsysSigma6VdW;
+ }
}
}
- /* only use softcore if one of the states has a zero endstate - softcore is for avoiding infinities!*/
- if ((c12[0] > 0) && (c12[1] > 0))
+ // NOLINTNEXTLINE(readability-misleading-indentation) remove when clang-tidy-13 is required
+ if constexpr (softcoreType == KernelSoftcoreType::Beutler)
{
- alpha_vdw_eff = 0;
- alpha_coul_eff = 0;
+ /* only use softcore if one of the states has a zero endstate - softcore is for avoiding infinities!*/
+ if ((c12[0] > 0) && (c12[1] > 0))
+ {
+ alpha_vdw_eff = 0;
+ alpha_coul_eff = 0;
+ }
+ else
+ {
+ alpha_vdw_eff = scParams.alphaVdw;
+ alpha_coul_eff = scParams.alphaCoulomb;
+ }
}
- else
+ // NOLINTNEXTLINE(readability-misleading-indentation) remove when clang-tidy-13 is required
+ if constexpr (softcoreType == KernelSoftcoreType::Gapsys)
{
- alpha_vdw_eff = alpha_vdw;
- alpha_coul_eff = alpha_coul;
+ /* only use softcore if one of the states has a zero endstate - softcore is for avoiding infinities!*/
+ if ((c12[0] > 0) && (c12[1] > 0))
+ {
+ gapsysScaleLinpointVdW = 0;
+ gapsysScaleLinpointCoul = 0;
+ }
+ else
+ {
+ gapsysScaleLinpointVdW = scParams.gapsysScaleLinpointVdW;
+ gapsysScaleLinpointCoul = scParams.gapsysScaleLinpointCoul;
+ }
}
/* Loop over A and B states again */
+ // NOLINTNEXTLINE(readability-misleading-indentation) remove when clang-tidy-13 is required
for (i = 0; i < 2; i++)
{
fscal_elec[i] = 0;
fscal_vdw[i] = 0;
velec[i] = 0;
vvdw[i] = 0;
+ dvdl_elec[i] = 0;
+ dvdl_vdw[i] = 0;
+ rQ = 0;
+ rLJ = 0;
/* Only spend time on A or B state if it is non-zero */
if ((qq[i] != 0) || (c6[i] != 0) || (c12[i] != 0))
{
/* Coulomb */
- rpinv = one / (alpha_coul_eff * lfac_coul[i] * sigma_pow[i] + rp);
- r_coul = std::pow(rpinv, minusOne / sc_r_power);
-
- /* Electrostatics table lookup data */
- rtab = r_coul * tabscale;
- ntab = static_cast<int>(rtab);
- eps = rtab - ntab;
- eps2 = eps * eps;
- ntab = static_cast<int>(tableStride * ntab);
- /* Electrostatics */
- Y = vftab[ntab];
- F = vftab[ntab + 1];
- Geps = eps * vftab[ntab + 2];
- Heps2 = eps2 * vftab[ntab + 3];
- Fp = F + Geps + Heps2;
- VV = Y + eps * Fp;
- FF = Fp + Geps + two * Heps2;
- velec[i] = qq[i] * VV;
- fscal_elec[i] = -qq[i] * FF * r_coul * rpinv * tabscale;
+ if constexpr (softcoreType == KernelSoftcoreType::Beutler)
+ {
+ rpinv = one / (alpha_coul_eff * lfac_coul[i] * sigma_pow[i] + rp);
+ r_coul = sixthRoot(rpinv);
+ }
+ else
+ {
+ rpinv = one / rp;
+ r_coul = r;
+ }
+
+ // NOLINTNEXTLINE(readability-misleading-indentation) remove when clang-tidy-13 is required
+ if constexpr (softcoreType == KernelSoftcoreType::Gapsys)
+ {
+ if ((facel != 0) && (LFC[i] < 1))
+ {
+ rQ = gmx::sixthroot(one - LFC[i]) * (one + std::fabs(qq[i] / facel));
+ rQ *= gapsysScaleLinpointCoul;
+ }
+ else
+ {
+ rQ = 0;
+ }
+ scaleDvdlRCoul = 1;
+ if (rQ > rCoulCutoff)
+ {
+ rQ = rCoulCutoff;
+ scaleDvdlRCoul = 0;
+ }
+ }
+
+ // NOLINTNEXTLINE(readability-misleading-indentation) remove when clang-tidy-13 is required
+ if ((softcoreType == KernelSoftcoreType::Gapsys) && (r < rQ))
+ {
+ real rInvQ = one / rQ;
+ real constFac = qq[i] * rInvQ;
+ real linFac = constFac * r * rInvQ;
+ real quadrFac = linFac * r * rInvQ;
+
+ /* Computing Coulomb force and potential energy */
+ fscal_elec[i] = -2 * quadrFac + 3 * linFac;
+ fscal_elec[i] *= rpinv;
+
+ velec[i] = quadrFac - 3 * (linFac - constFac);
+
+ dvdl_elec[i] += scaleDvdlRCoul * DLF[i] * half * (LFC[i] / (1 - LFC[i]))
+ * (quadrFac - 2 * linFac + constFac);
+ }
+ else // Beutler, resp. hardcore
+ {
+ /* Electrostatics table lookup data */
+ rtab = r_coul * tabscale;
+ ntab = static_cast<int>(rtab);
+ eps = rtab - ntab;
+ eps2 = eps * eps;
+ ntab = static_cast<int>(tableStride * ntab);
+ /* Electrostatics */
+ Y = vftab[ntab];
+ F = vftab[ntab + 1];
+ Geps = eps * vftab[ntab + 2];
+ Heps2 = eps2 * vftab[ntab + 3];
+ Fp = F + Geps + Heps2;
+ VV = Y + eps * Fp;
+ FF = Fp + Geps + two * Heps2;
+ velec[i] = qq[i] * VV;
+ fscal_elec[i] = -qq[i] * FF * r_coul * rpinv * tabscale;
+ }
/* Vdw */
- rpinv = one / (alpha_vdw_eff * lfac_vdw[i] * sigma_pow[i] + rp);
- r_vdw = std::pow(rpinv, minusOne / sc_r_power);
- /* Vdw table lookup data */
- rtab = r_vdw * tabscale;
- ntab = static_cast<int>(rtab);
- eps = rtab - ntab;
- eps2 = eps * eps;
- ntab = 12 * ntab;
- /* Dispersion */
- Y = vftab[ntab + 4];
- F = vftab[ntab + 5];
- Geps = eps * vftab[ntab + 6];
- Heps2 = eps2 * vftab[ntab + 7];
- Fp = F + Geps + Heps2;
- VV = Y + eps * Fp;
- FF = Fp + Geps + two * Heps2;
- vvdw[i] = c6[i] * VV;
- fscal_vdw[i] = -c6[i] * FF;
-
- /* Repulsion */
- Y = vftab[ntab + 8];
- F = vftab[ntab + 9];
- Geps = eps * vftab[ntab + 10];
- Heps2 = eps2 * vftab[ntab + 11];
- Fp = F + Geps + Heps2;
- VV = Y + eps * Fp;
- FF = Fp + Geps + two * Heps2;
- vvdw[i] += c12[i] * VV;
- fscal_vdw[i] -= c12[i] * FF;
- fscal_vdw[i] *= r_vdw * rpinv * tabscale;
+ if constexpr (softcoreType == KernelSoftcoreType::Beutler)
+ {
+ rpinv = one / (alpha_vdw_eff * lfac_vdw[i] * sigma_pow[i] + rp);
+ r_vdw = sixthRoot(rpinv);
+ }
+ else
+ {
+ rpinv = one / rp;
+ r_vdw = r;
+ }
+
+ if constexpr (softcoreType == KernelSoftcoreType::Gapsys)
+ {
+ constexpr real c_twentySixSeventh = 26.0_real / 7.0_real;
+
+ if (LFV[i] < 1)
+ {
+
+ rLJ = gmx::sixthroot(c_twentySixSeventh * gapsysSigma6VdW[i] * (one - LFV[i]));
+ rLJ *= gapsysScaleLinpointVdW;
+ }
+ else
+ {
+ rLJ = 0;
+ }
+ }
+
+ // NOLINTNEXTLINE(readability-misleading-indentation) remove when clang-tidy-13 is required
+ if ((softcoreType == KernelSoftcoreType::Gapsys) && (r < rLJ))
+ {
+ // scaled values for c6 and c12
+ real c6s, c12s;
+ c6s = c6[i] / 6.0_real;
+ c12s = c12[i] / 12.0_real;
+
+ /* Temporary variables for inverted values */
+ real rInvLJ = one / rLJ;
+ real rInv14, rInv13, rInv12;
+ real rInv8, rInv7, rInv6;
+ rInv6 = rInvLJ * rInvLJ * rInvLJ;
+ rInv6 *= rInv6;
+ rInv7 = rInv6 * rInvLJ;
+ rInv8 = rInv7 * rInvLJ;
+ rInv14 = c12s * rInv7 * rInv7 * r2;
+ rInv13 = c12s * rInv7 * rInv6 * r;
+ rInv12 = c12s * rInv6 * rInv6;
+ rInv8 *= c6s * r2;
+ rInv7 *= c6s * r;
+ rInv6 *= c6s;
+
+ /* Temporary variables for A and B */
+ real quadrFac, linearFac, constFac;
+ quadrFac = 156 * rInv14 - 42 * rInv8;
+ linearFac = 168 * rInv13 - 48 * rInv7;
+ constFac = 91 * rInv12 - 28 * rInv6;
+
+ /* Computing LJ force and potential energy*/
+ fscal_vdw[i] = -quadrFac + linearFac;
+ fscal_vdw[i] *= rpinv;
+
+ vvdw[i] = 0.5_real * quadrFac - linearFac + constFac;
+
+ dvdl_vdw[i] += DLF[i] * 28 * (LFV[i] / (one - LFV[i]))
+ * ((6.5_real * rInv14 - rInv8) - (13 * rInv13 - 2 * rInv7)
+ + (6.5_real * rInv12 - rInv6));
+ }
+ else // Beutler, resp. hardcore
+ {
+ /* Vdw table lookup data */
+ rtab = r_vdw * tabscale;
+ ntab = static_cast<int>(rtab);
+ eps = rtab - ntab;
+ eps2 = eps * eps;
+ ntab = 12 * ntab;
+ /* Dispersion */
+ Y = vftab[ntab + 4];
+ F = vftab[ntab + 5];
+ Geps = eps * vftab[ntab + 6];
+ Heps2 = eps2 * vftab[ntab + 7];
+ Fp = F + Geps + Heps2;
+ VV = Y + eps * Fp;
+ FF = Fp + Geps + two * Heps2;
+ vvdw[i] = c6[i] * VV;
+ fscal_vdw[i] = -c6[i] * FF;
+
+ /* Repulsion */
+ Y = vftab[ntab + 8];
+ F = vftab[ntab + 9];
+ Geps = eps * vftab[ntab + 10];
+ Heps2 = eps2 * vftab[ntab + 11];
+ Fp = F + Geps + Heps2;
+ VV = Y + eps * Fp;
+ FF = Fp + Geps + two * Heps2;
+ vvdw[i] += c12[i] * VV;
+ fscal_vdw[i] -= c12[i] * FF;
+ fscal_vdw[i] *= r_vdw * rpinv * tabscale;
+ }
}
}
/* Now we have velec[i], vvdw[i], and fscal[i] for both states */
/* Assemble A and B states */
- velecsum = 0;
- vvdwsum = 0;
- dvdl_coul = 0;
- dvdl_vdw = 0;
- fscal = 0;
+ velecsum = 0;
+ vvdwsum = 0;
+ dvdl_coul_sum = 0;
+ dvdl_vdw_sum = 0;
+ fscal = 0;
for (i = 0; i < 2; i++)
{
velecsum += LFC[i] * velec[i];
fscal += (LFC[i] * fscal_elec[i] + LFV[i] * fscal_vdw[i]) * rpm2;
- dvdl_coul += velec[i] * DLF[i]
- + LFC[i] * alpha_coul_eff * dlfac_coul[i] * fscal_elec[i] * sigma_pow[i];
- dvdl_vdw += vvdw[i] * DLF[i]
- + LFV[i] * alpha_vdw_eff * dlfac_vdw[i] * fscal_vdw[i] * sigma_pow[i];
+ if constexpr (softcoreType == KernelSoftcoreType::Gapsys)
+ {
+ dvdl_coul_sum += dvdl_elec[i];
+ dvdl_vdw_sum += dvdl_vdw[i];
+ }
+ // NOLINTNEXTLINE(readability-misleading-indentation) remove when clang-tidy-13 is required
+ dvdl_coul_sum += velec[i] * DLF[i];
+ dvdl_vdw_sum += vvdw[i] * DLF[i];
+ if constexpr (softcoreType == KernelSoftcoreType::Beutler)
+ {
+ dvdl_coul_sum += LFC[i] * alpha_coul_eff * dlfac_coul[i] * fscal_elec[i] * sigma_pow[i];
+ dvdl_vdw_sum += LFV[i] * alpha_vdw_eff * dlfac_vdw[i] * fscal_vdw[i] * sigma_pow[i];
+ }
}
- dvdl[efptCOUL] += dvdl_coul;
- dvdl[efptVDW] += dvdl_vdw;
+ dvdl[static_cast<int>(FreeEnergyPerturbationCouplingType::Coul)] += dvdl_coul_sum;
+ dvdl[static_cast<int>(FreeEnergyPerturbationCouplingType::Vdw)] += dvdl_vdw_sum;
*velectot = velecsum;
*vvdwtot = vvdwsum;
/*! \brief Calculate pair interactions, supports all types and conditions. */
template<BondedKernelFlavor flavor>
-static real do_pairs_general(int ftype,
- int nbonds,
- const t_iatom iatoms[],
- const t_iparams iparams[],
- const rvec x[],
- rvec4 f[],
- rvec fshift[],
- const struct t_pbc* pbc,
- const struct t_graph* g,
- const real* lambda,
- real* dvdl,
- const t_mdatoms* md,
- const t_forcerec* fr,
- gmx_grppairener_t* grppener,
- int* global_atom_index)
+static real do_pairs_general(int ftype,
+ int nbonds,
+ const t_iatom iatoms[],
+ const t_iparams iparams[],
+ const rvec x[],
+ rvec4 f[],
+ rvec fshift[],
+ const struct t_pbc* pbc,
+ const real* lambda,
+ real* dvdl,
+ gmx::ArrayRef<real> chargeA,
+ gmx::ArrayRef<real> chargeB,
+ gmx::ArrayRef<bool> atomIsPerturbed,
+ gmx::ArrayRef<unsigned short> cENER,
+ int numEnergyGroups,
+ const t_forcerec* fr,
+ gmx_grppairener_t* grppener,
+ int* global_atom_index)
{
real qq, c6, c12;
rvec dx;
- ivec dt;
int i, itype, ai, aj, gid;
int fshift_index;
real r2;
real* energygrp_vdw;
static gmx_bool warned_rlimit = FALSE;
/* Free energy stuff */
- gmx_bool bFreeEnergy;
- real LFC[2], LFV[2], DLF[2], lfac_coul[2], lfac_vdw[2], dlfac_coul[2], dlfac_vdw[2];
- real qqB, c6B, c12B, sigma2_def, sigma2_min;
+ gmx_bool bFreeEnergy;
+ real LFC[2], LFV[2], DLF[2], lfac_coul[2], lfac_vdw[2], dlfac_coul[2], dlfac_vdw[2];
+ real qqB, c6B, c12B;
+ const real oneSixth = 1.0_real / 6.0_real;
switch (ftype)
{
case F_LJ14:
case F_LJC14_Q:
- energygrp_elec = grppener->ener[egCOUL14].data();
- energygrp_vdw = grppener->ener[egLJ14].data();
+ energygrp_elec = grppener->energyGroupPairTerms[NonBondedEnergyTerms::Coulomb14].data();
+ energygrp_vdw = grppener->energyGroupPairTerms[NonBondedEnergyTerms::LJ14].data();
break;
case F_LJC_PAIRS_NB:
- energygrp_elec = grppener->ener[egCOULSR].data();
- energygrp_vdw = grppener->ener[egLJSR].data();
+ energygrp_elec = grppener->energyGroupPairTerms[NonBondedEnergyTerms::CoulombSR].data();
+ energygrp_vdw = grppener->energyGroupPairTerms[NonBondedEnergyTerms::LJSR].data();
break;
default:
energygrp_elec = nullptr; /* Keep compiler happy */
gmx_fatal(FARGS, "Unknown function type %d in do_nonbonded14", ftype);
}
- if (fr->efep != efepNO)
+ if (fr->efep != FreeEnergyPerturbationType::No)
{
/* Lambda factor for state A=1-lambda and B=lambda */
- LFC[0] = 1.0 - lambda[efptCOUL];
- LFV[0] = 1.0 - lambda[efptVDW];
- LFC[1] = lambda[efptCOUL];
- LFV[1] = lambda[efptVDW];
+ LFC[0] = 1.0 - lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Coul)];
+ LFV[0] = 1.0 - lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Vdw)];
+ LFC[1] = lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Coul)];
+ LFV[1] = lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Vdw)];
/*derivative of the lambda factor for state A and B */
DLF[0] = -1;
DLF[1] = 1;
- /* precalculate */
- sigma2_def = std::cbrt(fr->sc_sigma6_def);
- sigma2_min = std::cbrt(fr->sc_sigma6_min);
+ GMX_ASSERT(fr->ic->softCoreParameters, "We need soft-core parameters");
+ const auto& scParams = *fr->ic->softCoreParameters;
for (i = 0; i < 2; i++)
{
- lfac_coul[i] = (fr->sc_power == 2 ? (1 - LFC[i]) * (1 - LFC[i]) : (1 - LFC[i]));
- dlfac_coul[i] =
- DLF[i] * fr->sc_power / fr->sc_r_power * (fr->sc_power == 2 ? (1 - LFC[i]) : 1);
- lfac_vdw[i] = (fr->sc_power == 2 ? (1 - LFV[i]) * (1 - LFV[i]) : (1 - LFV[i]));
- dlfac_vdw[i] =
- DLF[i] * fr->sc_power / fr->sc_r_power * (fr->sc_power == 2 ? (1 - LFV[i]) : 1);
+ lfac_coul[i] = (scParams.lambdaPower == 2 ? (1 - LFC[i]) * (1 - LFC[i]) : (1 - LFC[i]));
+ dlfac_coul[i] = DLF[i] * scParams.lambdaPower * oneSixth
+ * (scParams.lambdaPower == 2 ? (1 - LFC[i]) : 1);
+ lfac_vdw[i] = (scParams.lambdaPower == 2 ? (1 - LFV[i]) * (1 - LFV[i]) : (1 - LFV[i]));
+ dlfac_vdw[i] = DLF[i] * scParams.lambdaPower * oneSixth
+ * (scParams.lambdaPower == 2 ? (1 - LFV[i]) : 1);
}
}
- else
- {
- sigma2_min = sigma2_def = 0;
- }
/* TODO This code depends on the logic in tables.c that constructs
the table layout, which should be made explicit in future
etiNR == 3,
"Pair-interaction code that uses GROMACS interaction tables supports exactly 3 tables");
GMX_ASSERT(
- fr->pairsTable->interaction == GMX_TABLE_INTERACTION_ELEC_VDWREP_VDWDISP,
+ fr->pairsTable->interaction == TableInteraction::ElectrostaticVdwRepulsionVdwDispersion,
"Pair interaction kernels need a table with Coulomb, repulsion and dispersion entries");
const real epsfac = fr->ic->epsfac;
itype = iatoms[i++];
ai = iatoms[i++];
aj = iatoms[i++];
- gid = GID(md->cENER[ai], md->cENER[aj], md->nenergrp);
+ gid = GID(cENER[ai], cENER[aj], numEnergyGroups);
/* Get parameters */
switch (ftype)
{
case F_LJ14:
- bFreeEnergy = (fr->efep != efepNO
- && (((md->nPerturbed != 0) && (md->bPerturbed[ai] || md->bPerturbed[aj]))
- || iparams[itype].lj14.c6A != iparams[itype].lj14.c6B
- || iparams[itype].lj14.c12A != iparams[itype].lj14.c12B));
- qq = md->chargeA[ai] * md->chargeA[aj] * epsfac * fr->fudgeQQ;
- c6 = iparams[itype].lj14.c6A;
- c12 = iparams[itype].lj14.c12A;
+ bFreeEnergy =
+ (fr->efep != FreeEnergyPerturbationType::No
+ && ((!atomIsPerturbed.empty() && (atomIsPerturbed[ai] || atomIsPerturbed[aj]))
+ || iparams[itype].lj14.c6A != iparams[itype].lj14.c6B
+ || iparams[itype].lj14.c12A != iparams[itype].lj14.c12B));
+ qq = chargeA[ai] * chargeA[aj] * epsfac * fr->fudgeQQ;
+ c6 = iparams[itype].lj14.c6A;
+ c12 = iparams[itype].lj14.c12A;
break;
case F_LJC14_Q:
qq = iparams[itype].ljc14.qi * iparams[itype].ljc14.qj * epsfac
}
else
{
- fshift_index = CENTRAL;
+ fshift_index = c_centralShiftIndex;
rvec_sub(x[ai], x[aj], dx);
}
r2 = norm2(dx);
- if (r2 >= fr->pairsTable->r * fr->pairsTable->r)
+ if (r2 >= fr->pairsTable->interactionRange * fr->pairsTable->interactionRange)
{
/* This check isn't race free. But it doesn't matter because if a race occurs the only
* disadvantage is that the warning is printed twice */
if (!warned_rlimit)
{
- warning_rlimit(x, ai, aj, global_atom_index, sqrt(r2), fr->pairsTable->r);
+ warning_rlimit(x, ai, aj, global_atom_index, sqrt(r2), fr->pairsTable->interactionRange);
warned_rlimit = TRUE;
}
continue;
if (bFreeEnergy)
{
/* Currently free energy is only supported for F_LJ14, so no need to check for that if we got here */
- qqB = md->chargeB[ai] * md->chargeB[aj] * epsfac * fr->fudgeQQ;
+ qqB = chargeB[ai] * chargeB[aj] * epsfac * fr->fudgeQQ;
c6B = iparams[itype].lj14.c6B * 6.0;
c12B = iparams[itype].lj14.c12B * 12.0;
- fscal = free_energy_evaluate_single(
- r2, fr->sc_r_power, fr->sc_alphacoul, fr->sc_alphavdw, fr->pairsTable->scale,
- fr->pairsTable->data, fr->pairsTable->stride, qq, c6, c12, qqB, c6B, c12B, LFC,
- LFV, DLF, lfac_coul, lfac_vdw, dlfac_coul, dlfac_vdw, fr->sc_sigma6_def,
- fr->sc_sigma6_min, sigma2_def, sigma2_min, &velec, &vvdw, dvdl);
+ const auto& scParams = *fr->ic->softCoreParameters;
+ if (scParams.softcoreType == SoftcoreType::Beutler)
+ {
+ if (scParams.alphaCoulomb == 0 && scParams.alphaVdw == 0)
+ {
+ fscal = free_energy_evaluate_single<KernelSoftcoreType::None>(
+ r2,
+ fr->ic->rcoulomb,
+ *fr->ic->softCoreParameters,
+ fr->pairsTable->scale,
+ fr->pairsTable->data.data(),
+ fr->pairsTable->stride,
+ qq,
+ c6,
+ c12,
+ qqB,
+ c6B,
+ c12B,
+ epsfac,
+ LFC,
+ LFV,
+ DLF,
+ lfac_coul,
+ lfac_vdw,
+ dlfac_coul,
+ dlfac_vdw,
+ &velec,
+ &vvdw,
+ dvdl);
+ }
+ else
+ {
+ fscal = free_energy_evaluate_single<KernelSoftcoreType::Beutler>(
+ r2,
+ fr->ic->rcoulomb,
+ *fr->ic->softCoreParameters,
+ fr->pairsTable->scale,
+ fr->pairsTable->data.data(),
+ fr->pairsTable->stride,
+ qq,
+ c6,
+ c12,
+ qqB,
+ c6B,
+ c12B,
+ epsfac,
+ LFC,
+ LFV,
+ DLF,
+ lfac_coul,
+ lfac_vdw,
+ dlfac_coul,
+ dlfac_vdw,
+ &velec,
+ &vvdw,
+ dvdl);
+ }
+ }
+ else // Gapsys
+ {
+ if (scParams.gapsysScaleLinpointCoul == 0 && scParams.gapsysScaleLinpointVdW == 0)
+ {
+ fscal = free_energy_evaluate_single<KernelSoftcoreType::None>(
+ r2,
+ fr->ic->rcoulomb,
+ *fr->ic->softCoreParameters,
+ fr->pairsTable->scale,
+ fr->pairsTable->data.data(),
+ fr->pairsTable->stride,
+ qq,
+ c6,
+ c12,
+ qqB,
+ c6B,
+ c12B,
+ epsfac,
+ LFC,
+ LFV,
+ DLF,
+ lfac_coul,
+ lfac_vdw,
+ dlfac_coul,
+ dlfac_vdw,
+ &velec,
+ &vvdw,
+ dvdl);
+ }
+ else
+ {
+ fscal = free_energy_evaluate_single<KernelSoftcoreType::Gapsys>(
+ r2,
+ fr->ic->rcoulomb,
+ *fr->ic->softCoreParameters,
+ fr->pairsTable->scale,
+ fr->pairsTable->data.data(),
+ fr->pairsTable->stride,
+ qq,
+ c6,
+ c12,
+ qqB,
+ c6B,
+ c12B,
+ epsfac,
+ LFC,
+ LFV,
+ DLF,
+ lfac_coul,
+ lfac_vdw,
+ dlfac_coul,
+ dlfac_vdw,
+ &velec,
+ &vvdw,
+ dvdl);
+ }
+ }
}
else
{
/* Evaluate tabulated interaction without free energy */
- fscal = evaluate_single(r2, fr->pairsTable->scale, fr->pairsTable->data,
- fr->pairsTable->stride, qq, c6, c12, &velec, &vvdw);
+ fscal = evaluate_single(r2,
+ fr->pairsTable->scale,
+ fr->pairsTable->data.data(),
+ fr->pairsTable->stride,
+ qq,
+ c6,
+ c12,
+ &velec,
+ &vvdw);
}
energygrp_elec[gid] += velec;
if (computeVirial(flavor))
{
- if (g)
- {
- /* Correct the shift forces using the graph */
- ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, aj), dt);
- fshift_index = IVEC2IS(dt);
- }
- if (fshift_index != CENTRAL)
+ if (fshift_index != c_centralShiftIndex)
{
rvec_inc(fshift[fshift_index], dx);
- rvec_dec(fshift[CENTRAL], dx);
+ rvec_dec(fshift[c_centralShiftIndex], dx);
}
}
}
* This function is templated for real/SimdReal and for optimization.
*/
template<typename T, int pack_size, typename pbc_type>
-static void do_pairs_simple(int nbonds,
- const t_iatom iatoms[],
- const t_iparams iparams[],
- const rvec x[],
- rvec4 f[],
- const pbc_type pbc,
- const t_mdatoms* md,
- const real scale_factor)
+static void do_pairs_simple(int nbonds,
+ const t_iatom iatoms[],
+ const t_iparams iparams[],
+ const rvec x[],
+ rvec4 f[],
+ const pbc_type pbc,
+ gmx::ArrayRef<real> charge,
+ const real scale_factor)
{
const int nfa1 = 1 + 2;
{
coeff[0 * pack_size + s] = iparams[itype].lj14.c6A;
coeff[1 * pack_size + s] = iparams[itype].lj14.c12A;
- coeff[2 * pack_size + s] = md->chargeA[ai[s]] * md->chargeA[aj[s]];
+ coeff[2 * pack_size + s] = charge[ai[s]] * charge[aj[s]];
/* Avoid indexing the iatoms array out of bounds.
* We pad the coordinate indices with the last atom pair.
}
/*! \brief Calculate all listed pair interactions */
-void do_pairs(int ftype,
- int nbonds,
- const t_iatom iatoms[],
- const t_iparams iparams[],
- const rvec x[],
- rvec4 f[],
- rvec fshift[],
- const struct t_pbc* pbc,
- const struct t_graph* g,
- const real* lambda,
- real* dvdl,
- const t_mdatoms* md,
- const t_forcerec* fr,
- const bool havePerturbedInteractions,
- const gmx::StepWorkload& stepWork,
- gmx_grppairener_t* grppener,
- int* global_atom_index)
+void do_pairs(int ftype,
+ int nbonds,
+ const t_iatom iatoms[],
+ const t_iparams iparams[],
+ const rvec x[],
+ rvec4 f[],
+ rvec fshift[],
+ const struct t_pbc* pbc,
+ const real* lambda,
+ real* dvdl,
+ gmx::ArrayRef<real> chargeA,
+ gmx::ArrayRef<real> chargeB,
+ gmx::ArrayRef<bool> atomIsPerturbed,
+ gmx::ArrayRef<unsigned short> cENER,
+ const int numEnergyGroups,
+ const t_forcerec* fr,
+ const bool havePerturbedInteractions,
+ const gmx::StepWorkload& stepWork,
+ gmx_grppairener_t* grppener,
+ int* global_atom_index)
{
- if (ftype == F_LJ14 && fr->ic->vdwtype != evdwUSER && !EEL_USER(fr->ic->eeltype)
+ if (ftype == F_LJ14 && fr->ic->vdwtype != VanDerWaalsType::User && !EEL_USER(fr->ic->eeltype)
&& !havePerturbedInteractions && (!stepWork.computeVirial && !stepWork.computeEnergy))
{
/* We use a fast code-path for plain LJ 1-4 without FEP.
set_pbc_simd(pbc, pbc_simd);
do_pairs_simple<SimdReal, GMX_SIMD_REAL_WIDTH, const real*>(
- nbonds, iatoms, iparams, x, f, pbc_simd, md, fr->ic->epsfac * fr->fudgeQQ);
+ nbonds, iatoms, iparams, x, f, pbc_simd, chargeA, fr->ic->epsfac * fr->fudgeQQ);
}
else
#endif
pbc_nonnull = &pbc_no;
}
- do_pairs_simple<real, 1, const t_pbc*>(nbonds, iatoms, iparams, x, f, pbc_nonnull, md,
- fr->ic->epsfac * fr->fudgeQQ);
+ do_pairs_simple<real, 1, const t_pbc*>(
+ nbonds, iatoms, iparams, x, f, pbc_nonnull, chargeA, fr->ic->epsfac * fr->fudgeQQ);
}
}
else if (stepWork.computeVirial)
{
- do_pairs_general<BondedKernelFlavor::ForcesAndVirialAndEnergy>(
- ftype, nbonds, iatoms, iparams, x, f, fshift, pbc, g, lambda, dvdl, md, fr,
- grppener, global_atom_index);
+ do_pairs_general<BondedKernelFlavor::ForcesAndVirialAndEnergy>(ftype,
+ nbonds,
+ iatoms,
+ iparams,
+ x,
+ f,
+ fshift,
+ pbc,
+ lambda,
+ dvdl,
+ chargeA,
+ chargeB,
+ atomIsPerturbed,
+ cENER,
+ numEnergyGroups,
+ fr,
+ grppener,
+ global_atom_index);
}
else
{
- do_pairs_general<BondedKernelFlavor::ForcesAndEnergy>(ftype, nbonds, iatoms, iparams, x, f,
- fshift, pbc, g, lambda, dvdl, md, fr,
- grppener, global_atom_index);
+ do_pairs_general<BondedKernelFlavor::ForcesAndEnergy>(ftype,
+ nbonds,
+ iatoms,
+ iparams,
+ x,
+ f,
+ fshift,
+ pbc,
+ lambda,
+ dvdl,
+ chargeA,
+ chargeB,
+ atomIsPerturbed,
+ cENER,
+ numEnergyGroups,
+ fr,
+ grppener,
+ global_atom_index);
}
}
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff