}
/*! \brief Compute the energy and force for a single pair interaction under FEP */
-template<SoftcoreType softcoreType>
+template<KernelSoftcoreType softcoreType>
static real free_energy_evaluate_single(real r2,
real rCoulCutoff,
const interaction_const_t::SoftCoreParameters& scParams,
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;
/* Loop over state A(0) and B(1) */
for (i = 0; i < 2; i++)
{
- if (softcoreType == SoftcoreType::Beutler || softcoreType == SoftcoreType::Gapsys)
+ if constexpr (softcoreType == KernelSoftcoreType::Beutler)
{
if ((c6[i] > 0) && (c12[i] > 0))
{
* 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) && softcoreType == SoftcoreType::Beutler) /* for disappearing coul and vdw with soft core at the same time */
+ if (sigma6[i] < scParams.sigma6Minimum) /* for disappearing coul and vdw with soft core at the same time */
{
sigma6[i] = scParams.sigma6Minimum;
}
}
sigma_pow[i] = sigma6[i];
}
+ // NOLINTNEXTLINE(readability-misleading-indentation) remove when clang-tidy-13 is required
+ if constexpr (softcoreType == KernelSoftcoreType::Gapsys)
+ {
+ 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;
+ }
+ }
}
- if (softcoreType == SoftcoreType::Beutler || softcoreType == SoftcoreType::Gapsys)
+ // NOLINTNEXTLINE(readability-misleading-indentation) remove when clang-tidy-13 is required
+ if constexpr (softcoreType == KernelSoftcoreType::Beutler)
{
/* only use softcore if one of the states has a zero endstate - softcore is for avoiding infinities!*/
if ((c12[0] > 0) && (c12[1] > 0))
}
else
{
- if (softcoreType == SoftcoreType::Beutler)
- {
- alpha_vdw_eff = scParams.alphaVdw;
- alpha_coul_eff = scParams.alphaCoulomb;
- }
- else if (softcoreType == SoftcoreType::Gapsys)
- {
- alpha_vdw_eff = scParams.alphaVdw;
- alpha_coul_eff = gmx::sixthroot(scParams.sigma6WithInvalidSigma);
- }
+ alpha_vdw_eff = scParams.alphaVdw;
+ alpha_coul_eff = scParams.alphaCoulomb;
+ }
+ }
+ // NOLINTNEXTLINE(readability-misleading-indentation) remove when clang-tidy-13 is required
+ if constexpr (softcoreType == KernelSoftcoreType::Gapsys)
+ {
+ /* 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;
if ((qq[i] != 0) || (c6[i] != 0) || (c12[i] != 0))
{
/* Coulomb */
- if (softcoreType == SoftcoreType::Beutler)
+ if constexpr (softcoreType == KernelSoftcoreType::Beutler)
{
rpinv = one / (alpha_coul_eff * lfac_coul[i] * sigma_pow[i] + rp);
r_coul = sixthRoot(rpinv);
r_coul = r;
}
- if (softcoreType == SoftcoreType::Gapsys)
+ // NOLINTNEXTLINE(readability-misleading-indentation) remove when clang-tidy-13 is required
+ if constexpr (softcoreType == KernelSoftcoreType::Gapsys)
{
- rQ = gmx::sixthroot(1. - LFC[i]) * (1. + std::fabs(qq[i] / facel));
- rQ *= alpha_coul_eff;
- scaleDvdlRCoul = 1.0;
+ 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.0;
+ scaleDvdlRCoul = 0;
}
}
- if ((softcoreType == SoftcoreType::Gapsys) && (r < rQ))
+ // NOLINTNEXTLINE(readability-misleading-indentation) remove when clang-tidy-13 is required
+ if ((softcoreType == KernelSoftcoreType::Gapsys) && (r < rQ))
{
- real rInvQ = 1.0 / 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] = -2 * quadrFac + 3 * linFac;
fscal_elec[i] *= rpinv;
- velec[i] = quadrFac - 3. * (linFac - constFac);
+ velec[i] = quadrFac - 3 * (linFac - constFac);
- dvdl_elec[i] += scaleDvdlRCoul * DLF[i] * 0.5 * (LFC[i] / (1. - LFC[i]))
- * (quadrFac - 2. * linFac + constFac);
+ dvdl_elec[i] += scaleDvdlRCoul * DLF[i] * half * (LFC[i] / (1 - LFC[i]))
+ * (quadrFac - 2 * linFac + constFac);
}
else // Beutler, resp. hardcore
{
}
/* Vdw */
- if (softcoreType == SoftcoreType::Beutler)
+ if constexpr (softcoreType == KernelSoftcoreType::Beutler)
{
rpinv = one / (alpha_vdw_eff * lfac_vdw[i] * sigma_pow[i] + rp);
r_vdw = sixthRoot(rpinv);
r_vdw = r;
}
- if (softcoreType == SoftcoreType::Gapsys)
+ if constexpr (softcoreType == KernelSoftcoreType::Gapsys)
{
- constexpr real c_twentySixSeventh = 26.0 / 7.0;
+ constexpr real c_twentySixSeventh = 26.0_real / 7.0_real;
- rLJ = gmx::sixthroot(c_twentySixSeventh * sigma6[i] * (1. - LFV[i]));
- rLJ *= alpha_vdw_eff;
+ if (LFV[i] < 1)
+ {
+
+ rLJ = gmx::sixthroot(c_twentySixSeventh * gapsysSigma6VdW[i] * (one - LFV[i]));
+ rLJ *= gapsysScaleLinpointVdW;
+ }
+ else
+ {
+ rLJ = 0;
+ }
}
- if ((softcoreType == SoftcoreType::Gapsys) && (r < rLJ))
+ // 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;
- c12s = c12[i] / 12.0;
+ c6s = c6[i] / 6.0_real;
+ c12s = c12[i] / 12.0_real;
/* Temporary variables for inverted values */
- real rInvLJ = 1.0 / rLJ;
+ real rInvLJ = one / rLJ;
real rInv14, rInv13, rInv12;
real rInv8, rInv7, rInv6;
rInv6 = rInvLJ * rInvLJ * rInvLJ;
/* 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;
+ 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] = -quadrFac + linearFac;
fscal_vdw[i] *= rpinv;
- vvdw[i] = 0.5 * quadrFac - linearFac + constFac;
+ vvdw[i] = 0.5_real * quadrFac - linearFac + constFac;
- dvdl_vdw[i] += DLF[i] * 28. * (LFV[i] / (1. - LFV[i]))
- * ((6.5 * rInv14 - rInv8) - (13. * rInv13 - 2. * rInv7)
- + (6.5 * rInv12 - rInv6));
+ 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
{
fscal += (LFC[i] * fscal_elec[i] + LFV[i] * fscal_vdw[i]) * rpm2;
- if (softcoreType == SoftcoreType::Gapsys)
+ 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 (softcoreType == SoftcoreType::Beutler)
+ 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];
c6B = iparams[itype].lj14.c6B * 6.0;
c12B = iparams[itype].lj14.c12B * 12.0;
- if (fr->ic->softCoreParameters->softcoreType == SoftcoreType::Gapsys)
+ const auto& scParams = *fr->ic->softCoreParameters;
+ if (scParams.softcoreType == SoftcoreType::Beutler)
{
- fscal = free_energy_evaluate_single<SoftcoreType::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);
+ 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 if (fr->ic->softCoreParameters->softcoreType == SoftcoreType::Beutler)
+ else // Gapsys
{
- fscal = free_energy_evaluate_single<SoftcoreType::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
- {
- fscal = free_energy_evaluate_single<SoftcoreType::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);
+ 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
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff