From 7064faa9a2d967a34df366ec97b304bd79ef8b32 Mon Sep 17 00:00:00 2001 From: Magnus Lundborg Date: Thu, 7 Nov 2019 14:47:38 +0100 Subject: [PATCH] Removed SC RPower 48 from FEP kernel With RPower 48 deprecated the FEP kernel has been simplified quite a lot. Fixes #3253. Change-Id: Ic18d9f938679b2a7499fbcc31bfc9d28e89fdb98 --- docs/user-guide/mdp-options.rst | 4 - src/gromacs/fileio/tpxio.cpp | 4 + .../gmxlib/nonbonded/nb_free_energy.cpp | 203 +++++------------- src/gromacs/gmxpreprocess/readir.cpp | 17 +- src/gromacs/listed_forces/pairs.cpp | 25 +-- 5 files changed, 67 insertions(+), 186 deletions(-) diff --git a/docs/user-guide/mdp-options.rst b/docs/user-guide/mdp-options.rst index 35d783d61c..282fc9abd1 100644 --- a/docs/user-guide/mdp-options.rst +++ b/docs/user-guide/mdp-options.rst @@ -2457,10 +2457,6 @@ Free energy calculations (6) power 6 for the radial term in the soft-core equation. - (48) - (deprecated) power 48 for the radial term in the soft-core equation. - Note that sc-alpha should generally be much lower (between 0.001 and 0.003). - .. mdp:: sc-coul (no) diff --git a/src/gromacs/fileio/tpxio.cpp b/src/gromacs/fileio/tpxio.cpp index c39d7f0edb..301d07bb57 100644 --- a/src/gromacs/fileio/tpxio.cpp +++ b/src/gromacs/fileio/tpxio.cpp @@ -562,6 +562,10 @@ static void do_fepvals(gmx::ISerializer* serializer, t_lambda* fepvals, int file { fepvals->sc_r_power = 6.0; } + if (fepvals->sc_r_power != 6.0) + { + gmx_fatal(FARGS, "sc-r-power=48 is no longer supported"); + } serializer->doReal(&fepvals->sc_sigma); if (serializer->reading()) { diff --git a/src/gromacs/gmxlib/nonbonded/nb_free_energy.cpp b/src/gromacs/gmxlib/nonbonded/nb_free_energy.cpp index d3d5ee60ba..ecd6aa7917 100644 --- a/src/gromacs/gmxlib/nonbonded/nb_free_energy.cpp +++ b/src/gromacs/gmxlib/nonbonded/nb_free_energy.cpp @@ -54,85 +54,17 @@ #include "gromacs/utility/fatalerror.h" -//! Enum for templating the soft-core treatment in the kernel -enum class SoftCoreTreatment -{ - None, //!< No soft-core - RPower6, //!< Soft-core with r-power = 6 - RPower48 //!< Soft-core with r-power = 48 -}; - -//! Most treatments are fine with float in mixed-precision mode. -template -struct SoftCoreReal -{ - //! Real type for soft-core calculations - using Real = real; -}; - -//! This treatment requires double precision for some computations. -template<> -struct SoftCoreReal -{ - //! Real type for soft-core calculations - using Real = double; -}; - //! Computes r^(1/p) and 1/r^(1/p) for the standard p=6 -template static inline void pthRoot(const real r, real* pthRoot, real* invPthRoot) { *invPthRoot = gmx::invsqrt(std::cbrt(r)); *pthRoot = 1 / (*invPthRoot); } -// We need a double version to make the specialization below work -#if !GMX_DOUBLE -//! Computes r^(1/p) and 1/r^(1/p) for the standard p=6 -template -static inline void pthRoot(const double r, real* pthRoot, double* invPthRoot) -{ - *invPthRoot = gmx::invsqrt(std::cbrt(r)); - *pthRoot = 1 / (*invPthRoot); -} -#endif - -//! Computes r^(1/p) and 1/r^(1/p) for p=48 -template<> -inline void pthRoot(const double r, real* pthRoot, double* invPthRoot) -{ - *pthRoot = std::pow(r, 1.0 / 48.0); - *invPthRoot = 1 / (*pthRoot); -} - -template -static inline real calculateSigmaPow(const real sigma6) -{ - if (softCoreTreatment == SoftCoreTreatment::RPower6) - { - return sigma6; - } - else - { - real sigmaPow = sigma6 * sigma6; /* sigma^12 */ - sigmaPow = sigmaPow * sigmaPow; /* sigma^24 */ - sigmaPow = sigmaPow * sigmaPow; /* sigma^48 */ - return (sigmaPow); - } -} - -template -static inline real calculateRinv6(const SCReal rinvV) +static inline real calculateRinv6(const real rinvV) { - if (softCoreTreatment == SoftCoreTreatment::RPower6) - { - return rinvV; - } - else - { - real rinv6 = rinvV * rinvV; - return (rinv6 * rinv6 * rinv6); - } + real rinv6 = rinvV * rinvV; + return (rinv6 * rinv6 * rinv6); } static inline real calculateVdw6(const real c6, const real rinv6) @@ -146,15 +78,12 @@ static inline real calculateVdw12(const real c12, const real rinv6) } /* reaction-field electrostatics */ -template -static inline SCReal -reactionFieldScalarForce(const real qq, const real rinv, const SCReal r, const real krf, const real two) +static inline real +reactionFieldScalarForce(const real qq, const real rinv, const real r, const real krf, const real two) { return (qq * (rinv - two * krf * r * r)); } -template -static inline real -reactionFieldPotential(const real qq, const real rinv, const SCReal r, const real krf, const real potentialShift) +static inline real reactionFieldPotential(const real qq, const real rinv, const real r, const real krf, const real potentialShift) { return (qq * (rinv + krf * r * r - potentialShift)); } @@ -193,25 +122,23 @@ static inline real ewaldLennardJonesGridSubtract(const real c6grid, const real p } /* LJ Potential switch */ -template -static inline SCReal potSwitchScalarForceMod(const SCReal fScalarInp, - const real potential, - const real sw, - const SCReal r, - const real rVdw, - const real dsw, - const real zero) +static inline real potSwitchScalarForceMod(const real fScalarInp, + const real potential, + const real sw, + const real r, + const real rVdw, + const real dsw, + const real zero) { if (r < rVdw) { - SCReal fScalar = fScalarInp * sw - r * potential * dsw; + real fScalar = fScalarInp * sw - r * potential * dsw; return (fScalar); } return (zero); } -template static inline real -potSwitchPotentialMod(const real potentialInp, const real sw, const SCReal r, const real rVdw, const real zero) +potSwitchPotentialMod(const real potentialInp, const real sw, const real r, const real rVdw, const real zero) { if (r < rVdw) { @@ -223,7 +150,7 @@ potSwitchPotentialMod(const real potentialInp, const real sw, const SCReal r, co //! Templated free-energy non-bonded kernel -template +template static void nb_free_energy_kernel(const t_nblist* gmx_restrict nlist, rvec* gmx_restrict xx, gmx::ForceWithShiftForces* forceWithShiftForces, @@ -232,10 +159,6 @@ static void nb_free_energy_kernel(const t_nblist* gmx_restrict nlist, nb_kernel_data_t* gmx_restrict kernel_data, t_nrnb* gmx_restrict nrnb) { - using SCReal = typename SoftCoreReal::Real; - - constexpr bool useSoftCore = (softCoreTreatment != SoftCoreTreatment::None); - #define STATE_A 0 #define STATE_B 1 #define NSTATES 2 @@ -357,8 +280,8 @@ static void nb_free_energy_kernel(const t_nblist* gmx_restrict nlist, GMX_RELEASE_ASSERT(!(vdwInteractionTypeIsEwald && vdwModifierIsPotSwitch), "Can not apply soft-core to switched Ewald potentials"); - SCReal dvdl_coul = 0; /* Needs double for sc_power==48 */ - SCReal dvdl_vdw = 0; /* Needs double for sc_power==48 */ + real dvdl_coul = 0; + real dvdl_vdw = 0; /* Lambda factor for state A, 1-lambda*/ real LFC[NSTATES], LFV[NSTATES]; @@ -375,7 +298,7 @@ static void nb_free_energy_kernel(const t_nblist* gmx_restrict nlist, DLF[STATE_B] = 1; real lfac_coul[NSTATES], dlfac_coul[NSTATES], lfac_vdw[NSTATES], dlfac_vdw[NSTATES]; - constexpr real sc_r_power = (softCoreTreatment == SoftCoreTreatment::RPower48 ? 48.0_real : 6.0_real); + constexpr real sc_r_power = 6.0_real; for (int i = 0; i < NSTATES; i++) { lfac_coul[i] = (lam_power == 2 ? (1 - LFC[i]) * (1 - LFC[i]) : (1 - LFC[i])); @@ -421,12 +344,12 @@ static void nb_free_energy_kernel(const t_nblist* gmx_restrict nlist, const int j3 = 3 * jnr; real c6[NSTATES], c12[NSTATES], qq[NSTATES], Vcoul[NSTATES], Vvdw[NSTATES]; real r, rinv, rp, rpm2; - real alpha_vdw_eff, alpha_coul_eff, sigma_pow[NSTATES]; + real alpha_vdw_eff, alpha_coul_eff, sigma6[NSTATES]; const real dx = ix - x[j3]; const real dy = iy - x[j3 + 1]; const real dz = iz - x[j3 + 2]; const real rsq = dx * dx + dy * dy + dz * dz; - SCReal FscalC[NSTATES], FscalV[NSTATES]; /* Needs double for sc_power==48 */ + real FscalC[NSTATES], FscalV[NSTATES]; if (rsq >= rcutoff_max2) { @@ -460,7 +383,12 @@ static void nb_free_energy_kernel(const t_nblist* gmx_restrict nlist, r = 0; } - if (softCoreTreatment == SoftCoreTreatment::None) + if (useSoftCore) + { + rpm2 = rsq * rsq; /* r4 */ + rp = rpm2 * rsq; /* r6 */ + } + else { /* The soft-core power p will not affect the results * with not using soft-core, so we use power of 0 which gives @@ -469,19 +397,6 @@ static void nb_free_energy_kernel(const t_nblist* gmx_restrict nlist, rpm2 = rinv * rinv; rp = 1; } - if (softCoreTreatment == SoftCoreTreatment::RPower6) - { - rpm2 = rsq * rsq; /* r4 */ - rp = rpm2 * rsq; /* r6 */ - } - if (softCoreTreatment == SoftCoreTreatment::RPower48) - { - rp = rsq * rsq * rsq; /* r6 */ - rp = rp * rp; /* r12 */ - rp = rp * rp; /* r24 */ - rp = rp * rp; /* r48 */ - rpm2 = rp / rsq; /* r46 */ - } real Fscal = 0; @@ -501,7 +416,6 @@ static void nb_free_energy_kernel(const t_nblist* gmx_restrict nlist, c12[i] = nbfp[tj[i] + 1]; if (useSoftCore) { - real sigma6[NSTATES]; if ((c6[i] > 0) && (c12[i] > 0)) { /* c12 is stored scaled with 12.0 and c6 is scaled with 6.0 - correct for this */ @@ -515,7 +429,6 @@ static void nb_free_energy_kernel(const t_nblist* gmx_restrict nlist, { sigma6[i] = sigma6_def; } - sigma_pow[i] = calculateSigmaPow(sigma6[i]); } } @@ -541,21 +454,20 @@ static void nb_free_energy_kernel(const t_nblist* gmx_restrict nlist, Vcoul[i] = 0; Vvdw[i] = 0; - real rinvC, rinvV; - SCReal rC, rV, rpinvC, rpinvV; /* Needs double for sc_power==48 */ + real rinvC, rinvV, rC, rV, rpinvC, rpinvV; /* Only spend time on A or B state if it is non-zero */ if ((qq[i] != 0) || (c6[i] != 0) || (c12[i] != 0)) { - /* this section has to be inside the loop because of the dependence on sigma_pow */ + /* this section has to be inside the loop because of the dependence on sigma6 */ if (useSoftCore) { - rpinvC = one / (alpha_coul_eff * lfac_coul[i] * sigma_pow[i] + rp); - pthRoot(rpinvC, &rinvC, &rC); + rpinvC = one / (alpha_coul_eff * lfac_coul[i] * sigma6[i] + rp); + pthRoot(rpinvC, &rinvC, &rC); if (scLambdasOrAlphasDiffer) { - rpinvV = one / (alpha_vdw_eff * lfac_vdw[i] * sigma_pow[i] + rp); - pthRoot(rpinvV, &rinvV, &rV); + rpinvV = one / (alpha_vdw_eff * lfac_vdw[i] * sigma6[i] + rp); + pthRoot(rpinvV, &rinvV, &rV); } else { @@ -607,13 +519,13 @@ static void nb_free_energy_kernel(const t_nblist* gmx_restrict nlist, if ((c6[i] != 0 || c12[i] != 0) && computeVdwInteraction) { real rinv6; - if (softCoreTreatment == SoftCoreTreatment::RPower6) + if (useSoftCore) { - rinv6 = calculateRinv6(rpinvV); + rinv6 = rpinvV; } else { - rinv6 = calculateRinv6(rinvV); + rinv6 = calculateRinv6(rinvV); } real Vvdw6 = calculateVdw6(c6[i], rinv6); real Vvdw12 = calculateVdw12(c12[i], rinv6); @@ -664,11 +576,10 @@ static void nb_free_energy_kernel(const t_nblist* gmx_restrict nlist, if (useSoftCore) { - dvdl_coul += - Vcoul[i] * DLF[i] - + LFC[i] * alpha_coul_eff * dlfac_coul[i] * FscalC[i] * sigma_pow[i]; + dvdl_coul += Vcoul[i] * DLF[i] + + LFC[i] * alpha_coul_eff * dlfac_coul[i] * FscalC[i] * sigma6[i]; dvdl_vdw += Vvdw[i] * DLF[i] - + LFV[i] * alpha_vdw_eff * dlfac_vdw[i] * FscalV[i] * sigma_pow[i]; + + LFV[i] * alpha_vdw_eff * dlfac_vdw[i] * FscalV[i] * sigma6[i]; } else { @@ -866,55 +777,55 @@ typedef void (*KernelFunction)(const t_nblist* gmx_restrict nlist, nb_kernel_data_t* gmx_restrict kernel_data, t_nrnb* gmx_restrict nrnb); -template +template static KernelFunction dispatchKernelOnVdwModifier(const bool vdwModifierIsPotSwitch) { if (vdwModifierIsPotSwitch) { - return (nb_free_energy_kernel); } else { - return (nb_free_energy_kernel); } } -template +template static KernelFunction dispatchKernelOnElecInteractionType(const bool elecInteractionTypeIsEwald, const bool vdwModifierIsPotSwitch) { if (elecInteractionTypeIsEwald) { - return (dispatchKernelOnVdwModifier( + return (dispatchKernelOnVdwModifier( vdwModifierIsPotSwitch)); } else { - return (dispatchKernelOnVdwModifier( + return (dispatchKernelOnVdwModifier( vdwModifierIsPotSwitch)); } } -template +template static KernelFunction dispatchKernelOnVdwInteractionType(const bool vdwInteractionTypeIsEwald, const bool elecInteractionTypeIsEwald, const bool vdwModifierIsPotSwitch) { if (vdwInteractionTypeIsEwald) { - return (dispatchKernelOnElecInteractionType( + return (dispatchKernelOnElecInteractionType( elecInteractionTypeIsEwald, vdwModifierIsPotSwitch)); } else { - return (dispatchKernelOnElecInteractionType( + return (dispatchKernelOnElecInteractionType( elecInteractionTypeIsEwald, vdwModifierIsPotSwitch)); } } -template +template static KernelFunction dispatchKernelOnScLambdasOrAlphasDifference(const bool scLambdasOrAlphasDiffer, const bool vdwInteractionTypeIsEwald, const bool elecInteractionTypeIsEwald, @@ -922,12 +833,12 @@ static KernelFunction dispatchKernelOnScLambdasOrAlphasDifference(const bool scL { if (scLambdasOrAlphasDiffer) { - return (dispatchKernelOnVdwInteractionType( + return (dispatchKernelOnVdwInteractionType( vdwInteractionTypeIsEwald, elecInteractionTypeIsEwald, vdwModifierIsPotSwitch)); } else { - return (dispatchKernelOnVdwInteractionType( + return (dispatchKernelOnVdwInteractionType( vdwInteractionTypeIsEwald, elecInteractionTypeIsEwald, vdwModifierIsPotSwitch)); } } @@ -940,19 +851,13 @@ static KernelFunction dispatchKernel(const bool scLambdasOrAlphasDiffer, { if (fr->sc_alphacoul == 0 && fr->sc_alphavdw == 0) { - return (dispatchKernelOnScLambdasOrAlphasDifference( - scLambdasOrAlphasDiffer, vdwInteractionTypeIsEwald, elecInteractionTypeIsEwald, - vdwModifierIsPotSwitch)); - } - else if (fr->sc_r_power == 6.0_real) - { - return (dispatchKernelOnScLambdasOrAlphasDifference( + return (dispatchKernelOnScLambdasOrAlphasDifference( scLambdasOrAlphasDiffer, vdwInteractionTypeIsEwald, elecInteractionTypeIsEwald, vdwModifierIsPotSwitch)); } else { - return (dispatchKernelOnScLambdasOrAlphasDifference( + return (dispatchKernelOnScLambdasOrAlphasDifference( scLambdasOrAlphasDiffer, vdwInteractionTypeIsEwald, elecInteractionTypeIsEwald, vdwModifierIsPotSwitch)); } @@ -979,7 +884,7 @@ void gmx_nb_free_energy_kernel(const t_nblist* nlist, { scLambdasOrAlphasDiffer = false; } - else if (fr->sc_r_power == 6.0_real || fr->sc_r_power == 48.0_real) + else if (fr->sc_r_power == 6.0_real) { if (kernel_data->lambda[efptCOUL] == kernel_data->lambda[efptVDW] && fr->sc_alphacoul == fr->sc_alphavdw) { diff --git a/src/gromacs/gmxpreprocess/readir.cpp b/src/gromacs/gmxpreprocess/readir.cpp index 657fcea1e6..f4699da88f 100644 --- a/src/gromacs/gmxpreprocess/readir.cpp +++ b/src/gromacs/gmxpreprocess/readir.cpp @@ -667,9 +667,11 @@ void check_ir(const char* mdparin, sprintf(err_buf, "The soft-core power is %d and can only be 1 or 2", fep->sc_power); CHECK(fep->sc_alpha != 0 && fep->sc_power != 1 && fep->sc_power != 2); - sprintf(err_buf, "The soft-core sc-r-power is %d and can only be 6 or 48", + sprintf(err_buf, + "The soft-core sc-r-power is %d and can only be 6. (sc-r-power 48 is no longer " + "supported.)", static_cast(fep->sc_r_power)); - CHECK(fep->sc_alpha != 0 && fep->sc_r_power != 6.0 && fep->sc_r_power != 48.0); + CHECK(fep->sc_alpha != 0 && fep->sc_r_power != 6.0); sprintf(err_buf, "Can't use positive delta-lambda (%g) if initial state/lambda does not start at " @@ -1566,17 +1568,6 @@ static void do_fep_params(t_inputrec* ir, char fep_lambda[][STRLEN], char weight } - /* make it easier if sc_r_power = 48 by increasing it to the 4th power, to be in the right scale. */ - if (fep->sc_r_power == 48) - { - if (fep->sc_alpha > 0.1) - { - gmx_fatal(FARGS, - "sc_alpha (%f) for sc_r_power = 48 should usually be between 0.001 and 0.004", - fep->sc_alpha); - } - } - /* now read in the weights */ parse_n_real(weights, &nweights, &(expand->init_lambda_weights), wi); if (nweights == 0) diff --git a/src/gromacs/listed_forces/pairs.cpp b/src/gromacs/listed_forces/pairs.cpp index c3fee5885f..bffbcb42ef 100644 --- a/src/gromacs/listed_forces/pairs.cpp +++ b/src/gromacs/listed_forces/pairs.cpp @@ -188,12 +188,11 @@ static real free_energy_evaluate_single(real r2, real fscal_vdw[2], fscal_elec[2]; real velec[2], vvdw[2]; 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; + const real minusOne = -1.0; + const real one = 1.0; + const real two = 2.0; + const real six = 6.0; qq[0] = qqA; qq[1] = qqB; @@ -207,14 +206,6 @@ static real free_energy_evaluate_single(real r2, 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 */ @@ -248,12 +239,6 @@ static real free_energy_evaluate_single(real r2, { sigma_pow[i] = sigma6[i]; } - else if (sc_r_power == fourtyeight) - { - 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); -- 2.22.0