From e6ef503126b61c7c5b06a0f6966484a85512cef1 Mon Sep 17 00:00:00 2001 From: Sebastian Kehl Date: Thu, 23 Sep 2021 15:31:11 +0200 Subject: [PATCH] Introduce explicit parameters for gapsys-sc. --- src/gromacs/fileio/tpxio.cpp | 8 +- .../gmxlib/nonbonded/nb_free_energy.cpp | 171 ++++++++---- src/gromacs/gmxpreprocess/readir.cpp | 29 +- ...arametersWithValuesIncludingAssignment.xml | 3 + .../GetIrTest_AcceptsElectricField.xml | 3 + ...IrTest_AcceptsElectricFieldOscillating.xml | 3 + .../GetIrTest_AcceptsElectricFieldPulsed.xml | 3 + .../refdata/GetIrTest_AcceptsEmptyLines.xml | 3 + .../GetIrTest_AcceptsImplicitSolventNo.xml | 3 + .../GetIrTest_AcceptsKeyWithoutValue.xml | 3 + .../tests/refdata/GetIrTest_AcceptsMimic.xml | 3 + ...IrTest_HandlesDifferentKindsOfMdpLines.xml | 3 + src/gromacs/listed_forces/pairs.cpp | 257 +++++++++++------- src/gromacs/mdtypes/inputrec.cpp | 20 +- src/gromacs/mdtypes/inputrec.h | 6 + src/gromacs/mdtypes/interaction_const.cpp | 5 +- src/gromacs/mdtypes/interaction_const.h | 6 + 17 files changed, 359 insertions(+), 170 deletions(-) diff --git a/src/gromacs/fileio/tpxio.cpp b/src/gromacs/fileio/tpxio.cpp index 33e63d4f11..39943911fa 100644 --- a/src/gromacs/fileio/tpxio.cpp +++ b/src/gromacs/fileio/tpxio.cpp @@ -626,10 +626,16 @@ static void do_fepvals(gmx::ISerializer* serializer, t_lambda* fepvals, int file if (file_version >= tpxv_SoftcoreGapsys) { serializer->doInt(reinterpret_cast(&fepvals->softcoreFunction)); + serializer->doReal(&fepvals->scScaleLinpointLJGapsys); + serializer->doReal(&fepvals->scScaleLinpointQGapsys); + serializer->doReal(&fepvals->scSigmaLJGapsys); } else { - fepvals->softcoreFunction = SoftcoreType::Beutler; + fepvals->softcoreFunction = SoftcoreType::Beutler; + fepvals->scScaleLinpointLJGapsys = 0.85; + fepvals->scScaleLinpointQGapsys = 0.3; + fepvals->scSigmaLJGapsys = 0.3; } /* handle lambda_neighbors */ diff --git a/src/gromacs/gmxlib/nonbonded/nb_free_energy.cpp b/src/gromacs/gmxlib/nonbonded/nb_free_energy.cpp index 33e485d191..41bcedc058 100644 --- a/src/gromacs/gmxlib/nonbonded/nb_free_energy.cpp +++ b/src/gromacs/gmxlib/nonbonded/nb_free_energy.cpp @@ -304,14 +304,21 @@ static void nb_free_energy_kernel(const t_nblist& gmx::ArrayRef shift = nlist.shift; gmx::ArrayRef gid = nlist.gid; - const real lambda_coul = lambda[static_cast(FreeEnergyPerturbationCouplingType::Coul)]; - const real lambda_vdw = lambda[static_cast(FreeEnergyPerturbationCouplingType::Vdw)]; - const auto& scParams = *ic.softCoreParameters; - const real gmx_unused alpha_coul = scParams.alphaCoulomb; - const real gmx_unused alpha_vdw = scParams.alphaVdw; - const real lam_power = scParams.lambdaPower; - const real gmx_unused sigma6_def = scParams.sigma6WithInvalidSigma; - const real gmx_unused sigma6_min = scParams.sigma6Minimum; + const real lambda_coul = lambda[static_cast(FreeEnergyPerturbationCouplingType::Coul)]; + const real lambda_vdw = lambda[static_cast(FreeEnergyPerturbationCouplingType::Vdw)]; + + // Extract softcore parameters + const auto& scParams = *ic.softCoreParameters; + const real lam_power = scParams.lambdaPower; + const real gmx_unused alpha_coul = scParams.alphaCoulomb; + const real gmx_unused alpha_vdw = scParams.alphaVdw; + const real gmx_unused sigma6_def = scParams.sigma6WithInvalidSigma; + const real gmx_unused sigma6_min = scParams.sigma6Minimum; + + const real gmx_unused scaleLinpointCoulGapsys = scParams.scaleLinpointCoulGapsys; + const real gmx_unused scaleLinpointVdWGapsys = scParams.scaleLinpointVdWGapsys; + const real gmx_unused sigma6VdWGapsys = scParams.sigma6VdWGapsys; + const bool gmx_unused doShiftForces = ((flags & GMX_NONBONDED_DO_SHIFTFORCE) != 0); const bool doPotential = ((flags & GMX_NONBONDED_DO_POTENTIAL) != 0); @@ -492,6 +499,12 @@ static void nb_free_energy_kernel(const t_nblist& alignas(GMX_SIMD_ALIGNMENT) real gmx_unused preloadSigma6[NSTATES][DataTypes::simdRealWidth]; alignas(GMX_SIMD_ALIGNMENT) real gmx_unused preloadAlphaVdwEff[DataTypes::simdRealWidth]; alignas(GMX_SIMD_ALIGNMENT) real gmx_unused preloadAlphaCoulEff[DataTypes::simdRealWidth]; + alignas(GMX_SIMD_ALIGNMENT) + real gmx_unused preloadScaleLinpointVdWGapsys[DataTypes::simdRealWidth]; + alignas(GMX_SIMD_ALIGNMENT) + real gmx_unused preloadScaleLinpointCoulGapsys[DataTypes::simdRealWidth]; + alignas(GMX_SIMD_ALIGNMENT) + real gmx_unused preloadSigma6VdWGapsys[NSTATES][DataTypes::simdRealWidth]; alignas(GMX_SIMD_ALIGNMENT) real preloadLjPmeC6Grid[NSTATES][DataTypes::simdRealWidth]; #else real preloadPairIsValid[DataTypes::simdRealWidth]; @@ -502,6 +515,9 @@ static void nb_free_energy_kernel(const t_nblist& real gmx_unused preloadSigma6[NSTATES][DataTypes::simdRealWidth]; real gmx_unused preloadAlphaVdwEff[DataTypes::simdRealWidth]; real gmx_unused preloadAlphaCoulEff[DataTypes::simdRealWidth]; + real gmx_unused preloadScaleLinpointVdWGapsys[DataTypes::simdRealWidth]; + real gmx_unused preloadScaleLinpointCoulGapsys[DataTypes::simdRealWidth]; + real gmx_unused preloadSigma6VdWGapsys[NSTATES][DataTypes::simdRealWidth]; real preloadLjPmeC6Grid[NSTATES][DataTypes::simdRealWidth]; #endif for (int s = 0; s < DataTypes::simdRealWidth; s++) @@ -528,8 +544,7 @@ static void nb_free_energy_kernel(const t_nblist& { preloadLjPmeC6Grid[i][s] = 0; } - if constexpr (softcoreType == KernelSoftcoreType::Beutler - || softcoreType == KernelSoftcoreType::Gapsys) + if constexpr (softcoreType == KernelSoftcoreType::Beutler) { const real c6 = nbfp[2 * typeIndices[i][s]]; const real c12 = nbfp[2 * typeIndices[i][s] + 1]; @@ -548,9 +563,26 @@ static void nb_free_energy_kernel(const t_nblist& preloadSigma6[i][s] = sigma6_def; } } + if constexpr (softcoreType == KernelSoftcoreType::Gapsys) + { + const real c6 = nbfp[2 * typeIndices[i][s]]; + const real c12 = nbfp[2 * typeIndices[i][s] + 1]; + if (c6 > 0 && c12 > 0) + { + /* c12 is stored scaled with 12.0 and c6 is scaled with 6.0 - correct for this */ + preloadSigma6VdWGapsys[i][s] = 0.5_real * c12 / c6; + if (preloadSigma6VdWGapsys[i][s] < sigma6VdWGapsys) + { + preloadSigma6VdWGapsys[i][s] = sigma6VdWGapsys; + } + } + else + { + preloadSigma6VdWGapsys[i][s] = sigma6VdWGapsys; + } + } } - if constexpr (softcoreType == KernelSoftcoreType::Beutler - || softcoreType == KernelSoftcoreType::Gapsys) + if constexpr (softcoreType == KernelSoftcoreType::Beutler) { /* only use softcore if one of the states has a zero endstate - softcore is for avoiding infinities!*/ const real c12A = nbfp[2 * typeIndices[STATE_A][s] + 1]; @@ -562,34 +594,45 @@ static void nb_free_energy_kernel(const t_nblist& } else { - if constexpr (softcoreType == KernelSoftcoreType::Beutler) - { - preloadAlphaVdwEff[s] = alpha_vdw; - preloadAlphaCoulEff[s] = alpha_coul; - } - else if constexpr (softcoreType == KernelSoftcoreType::Gapsys) - { - preloadAlphaVdwEff[s] = alpha_vdw; - preloadAlphaCoulEff[s] = gmx::sixthroot(sigma6_def); - } + preloadAlphaVdwEff[s] = alpha_vdw; + preloadAlphaCoulEff[s] = alpha_coul; + } + } + if constexpr (softcoreType == KernelSoftcoreType::Gapsys) + { + /* only use softcore if one of the states has a zero endstate - softcore is for avoiding infinities!*/ + const real c12A = nbfp[2 * typeIndices[STATE_A][s] + 1]; + const real c12B = nbfp[2 * typeIndices[STATE_B][s] + 1]; + if (c12A > 0 && c12B > 0) + { + preloadScaleLinpointVdWGapsys[s] = 0; + preloadScaleLinpointCoulGapsys[s] = 0; + } + else + { + preloadScaleLinpointVdWGapsys[s] = scaleLinpointVdWGapsys; + preloadScaleLinpointCoulGapsys[s] = scaleLinpointCoulGapsys; } } } else { - preloadJnr[s] = jjnr[k]; - preloadPairIsValid[s] = false; - preloadPairIncluded[s] = false; - preloadAlphaVdwEff[s] = 0; - preloadAlphaCoulEff[s] = 0; + preloadJnr[s] = jjnr[k]; + preloadPairIsValid[s] = false; + preloadPairIncluded[s] = false; + preloadAlphaVdwEff[s] = 0; + preloadAlphaCoulEff[s] = 0; + preloadScaleLinpointVdWGapsys[s] = 0; + preloadScaleLinpointCoulGapsys[s] = 0; for (int i = 0; i < NSTATES; i++) { - typeIndices[STATE_A][s] = ntiA + typeA[jjnr[k]]; - typeIndices[STATE_B][s] = ntiB + typeB[jjnr[k]]; - preloadLjPmeC6Grid[i][s] = 0; - preloadQq[i][s] = 0; - preloadSigma6[i][s] = 0; + typeIndices[STATE_A][s] = ntiA + typeA[jjnr[k]]; + typeIndices[STATE_B][s] = ntiB + typeB[jjnr[k]]; + preloadLjPmeC6Grid[i][s] = 0; + preloadQq[i][s] = 0; + preloadSigma6[i][s] = 0; + preloadSigma6VdWGapsys[i][s] = 0; } } } @@ -641,23 +684,33 @@ static void nb_free_energy_kernel(const t_nblist& RealType gmx_unused ljPmeC6Grid[NSTATES]; RealType gmx_unused alphaVdwEff; RealType gmx_unused alphaCoulEff; + RealType gmx_unused scaleLinpointVdWGapsysEff; + RealType gmx_unused scaleLinpointCoulGapsysEff; + RealType gmx_unused sigmaVdWGapsysEff[NSTATES]; for (int i = 0; i < NSTATES; i++) { gmx::gatherLoadTranspose<2>(nbfp.data(), typeIndices[i], &c6[i], &c12[i]); qq[i] = gmx::load(preloadQq[i]); ljPmeC6Grid[i] = gmx::load(preloadLjPmeC6Grid[i]); - if constexpr (softcoreType == KernelSoftcoreType::Beutler - || softcoreType == KernelSoftcoreType::Gapsys) + if constexpr (softcoreType == KernelSoftcoreType::Beutler) { sigma6[i] = gmx::load(preloadSigma6[i]); } + if constexpr (softcoreType == KernelSoftcoreType::Gapsys) + { + sigmaVdWGapsysEff[i] = gmx::load(preloadSigma6VdWGapsys[i]); + } } - if constexpr (softcoreType == KernelSoftcoreType::Beutler - || softcoreType == KernelSoftcoreType::Gapsys) + if constexpr (softcoreType == KernelSoftcoreType::Beutler) { alphaVdwEff = gmx::load(preloadAlphaVdwEff); alphaCoulEff = gmx::load(preloadAlphaCoulEff); } + if constexpr (softcoreType == KernelSoftcoreType::Gapsys) + { + scaleLinpointVdWGapsysEff = gmx::load(preloadScaleLinpointVdWGapsys); + scaleLinpointCoulGapsysEff = gmx::load(preloadScaleLinpointCoulGapsys); + } // Avoid overflow of r^-12 at distances near zero rSq = gmx::max(rSq, minDistanceSquared); @@ -765,7 +818,7 @@ static void nb_free_energy_kernel(const t_nblist& rCutoffCoul, LFC[i], DLF[i], - alphaCoulEff, + scaleLinpointCoulGapsysEff, sh_ewald, &fScalC[i], &vCoul[i], @@ -789,7 +842,7 @@ static void nb_free_energy_kernel(const t_nblist& rCutoffCoul, LFC[i], DLF[i], - alphaCoulEff, + scaleLinpointCoulGapsysEff, krf, crf, &fScalC[i], @@ -856,8 +909,8 @@ static void nb_free_energy_kernel(const t_nblist& rSq, LFV[i], DLF[i], - sigma6[i], - alphaVdwEff, + sigmaVdWGapsysEff[i], + scaleLinpointVdWGapsysEff, repulsionShift, dispersionShift, &fScalV[i], @@ -1242,21 +1295,12 @@ static KernelFunction dispatchKernel(const bool scLambdasOrAlpha const bool useSimd, const interaction_const_t& ic) { - if (ic.softCoreParameters->alphaCoulomb == 0 && ic.softCoreParameters->alphaVdw == 0) - { - return (dispatchKernelOnScLambdasOrAlphasDifference( - scLambdasOrAlphasDiffer, - vdwInteractionTypeIsEwald, - elecInteractionTypeIsEwald, - vdwModifierIsPotSwitch, - computeForces, - useSimd)); - } - else + const auto& scParams = *ic.softCoreParameters; + if (scParams.softcoreType == SoftcoreType::Beutler) { - if (ic.softCoreParameters->softcoreType == SoftcoreType::Beutler) + if (scParams.alphaCoulomb == 0 && scParams.alphaVdw == 0) { - return (dispatchKernelOnScLambdasOrAlphasDifference( + return (dispatchKernelOnScLambdasOrAlphasDifference( scLambdasOrAlphasDiffer, vdwInteractionTypeIsEwald, elecInteractionTypeIsEwald, @@ -1264,9 +1308,19 @@ static KernelFunction dispatchKernel(const bool scLambdasOrAlpha computeForces, useSimd)); } - else + return (dispatchKernelOnScLambdasOrAlphasDifference( + scLambdasOrAlphasDiffer, + vdwInteractionTypeIsEwald, + elecInteractionTypeIsEwald, + vdwModifierIsPotSwitch, + computeForces, + useSimd)); + } + else // Gapsys + { + if (scParams.scaleLinpointCoulGapsys == 0 && scParams.scaleLinpointVdWGapsys == 0) { - return (dispatchKernelOnScLambdasOrAlphasDifference( + return (dispatchKernelOnScLambdasOrAlphasDifference( scLambdasOrAlphasDiffer, vdwInteractionTypeIsEwald, elecInteractionTypeIsEwald, @@ -1274,6 +1328,13 @@ static KernelFunction dispatchKernel(const bool scLambdasOrAlpha computeForces, useSimd)); } + return (dispatchKernelOnScLambdasOrAlphasDifference( + scLambdasOrAlphasDiffer, + vdwInteractionTypeIsEwald, + elecInteractionTypeIsEwald, + vdwModifierIsPotSwitch, + computeForces, + useSimd)); } } diff --git a/src/gromacs/gmxpreprocess/readir.cpp b/src/gromacs/gmxpreprocess/readir.cpp index 99e02def68..c5873ab705 100644 --- a/src/gromacs/gmxpreprocess/readir.cpp +++ b/src/gromacs/gmxpreprocess/readir.cpp @@ -2322,19 +2322,22 @@ void get_ir(const char* mdparin, setStringEntry(&inp, "temperature-lambdas", ""); fep->lambda_neighbors = get_eint(&inp, "calc-lambda-neighbors", 1, wi); setStringEntry(&inp, "init-lambda-weights", inputrecStrings->lambda_weights, nullptr); - fep->edHdLPrintEnergy = getEnum(&inp, "dhdl-print-energy", wi); - fep->softcoreFunction = getEnum(&inp, "sc-function", wi); - fep->sc_alpha = get_ereal(&inp, "sc-alpha", 0.0, wi); - fep->sc_power = get_eint(&inp, "sc-power", 1, wi); - fep->sc_r_power = get_ereal(&inp, "sc-r-power", 6.0, wi); - fep->sc_sigma = get_ereal(&inp, "sc-sigma", 0.3, wi); - fep->bScCoul = (getEnum(&inp, "sc-coul", wi) != Boolean::No); - fep->dh_hist_size = get_eint(&inp, "dh_hist_size", 0, wi); - fep->dh_hist_spacing = get_ereal(&inp, "dh_hist_spacing", 0.1, wi); - fep->separate_dhdl_file = getEnum(&inp, "separate-dhdl-file", wi); - fep->dhdl_derivatives = getEnum(&inp, "dhdl-derivatives", wi); - fep->dh_hist_size = get_eint(&inp, "dh_hist_size", 0, wi); - fep->dh_hist_spacing = get_ereal(&inp, "dh_hist_spacing", 0.1, wi); + fep->edHdLPrintEnergy = getEnum(&inp, "dhdl-print-energy", wi); + fep->softcoreFunction = getEnum(&inp, "sc-function", wi); + fep->sc_alpha = get_ereal(&inp, "sc-alpha", 0.0, wi); + fep->sc_power = get_eint(&inp, "sc-power", 1, wi); + fep->sc_r_power = get_ereal(&inp, "sc-r-power", 6.0, wi); + fep->sc_sigma = get_ereal(&inp, "sc-sigma", 0.3, wi); + fep->bScCoul = (getEnum(&inp, "sc-coul", wi) != Boolean::No); + fep->scScaleLinpointLJGapsys = get_ereal(&inp, "sc-scale-linpoint-LJ-gapsys", 0.85, wi); + fep->scScaleLinpointQGapsys = get_ereal(&inp, "sc-scale-linpoint-Q-gapsys", 0.3, wi); + fep->scSigmaLJGapsys = get_ereal(&inp, "sc-sigma-LJ-gapsys", 0.3, wi); + fep->dh_hist_size = get_eint(&inp, "dh_hist_size", 0, wi); + fep->dh_hist_spacing = get_ereal(&inp, "dh_hist_spacing", 0.1, wi); + fep->separate_dhdl_file = getEnum(&inp, "separate-dhdl-file", wi); + fep->dhdl_derivatives = getEnum(&inp, "dhdl-derivatives", wi); + fep->dh_hist_size = get_eint(&inp, "dh_hist_size", 0, wi); + fep->dh_hist_spacing = get_ereal(&inp, "dh_hist_spacing", 0.1, wi); /* Non-equilibrium MD stuff */ printStringNewline(&inp, "Non-equilibrium MD stuff"); diff --git a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsDefineParametersWithValuesIncludingAssignment.xml b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsDefineParametersWithValuesIncludingAssignment.xml index 1c2029b4f2..ed5a7eb1fc 100644 --- a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsDefineParametersWithValuesIncludingAssignment.xml +++ b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsDefineParametersWithValuesIncludingAssignment.xml @@ -255,6 +255,9 @@ sc-power = 1 sc-r-power = 6 sc-sigma = 0.3 sc-coul = no +sc-scale-linpoint-LJ-gapsys = 0.85 +sc-scale-linpoint-Q-gapsys = 0.3 +sc-sigma-LJ-gapsys = 0.3 separate-dhdl-file = yes dhdl-derivatives = yes dh_hist_size = 0 diff --git a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsElectricField.xml b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsElectricField.xml index 21fe54ad8a..34dea51cca 100644 --- a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsElectricField.xml +++ b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsElectricField.xml @@ -255,6 +255,9 @@ sc-power = 1 sc-r-power = 6 sc-sigma = 0.3 sc-coul = no +sc-scale-linpoint-LJ-gapsys = 0.85 +sc-scale-linpoint-Q-gapsys = 0.3 +sc-sigma-LJ-gapsys = 0.3 separate-dhdl-file = yes dhdl-derivatives = yes dh_hist_size = 0 diff --git a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsElectricFieldOscillating.xml b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsElectricFieldOscillating.xml index e8c72614b2..9fc0ed9f50 100644 --- a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsElectricFieldOscillating.xml +++ b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsElectricFieldOscillating.xml @@ -255,6 +255,9 @@ sc-power = 1 sc-r-power = 6 sc-sigma = 0.3 sc-coul = no +sc-scale-linpoint-LJ-gapsys = 0.85 +sc-scale-linpoint-Q-gapsys = 0.3 +sc-sigma-LJ-gapsys = 0.3 separate-dhdl-file = yes dhdl-derivatives = yes dh_hist_size = 0 diff --git a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsElectricFieldPulsed.xml b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsElectricFieldPulsed.xml index 96b74dfec2..902a7b45fe 100644 --- a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsElectricFieldPulsed.xml +++ b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsElectricFieldPulsed.xml @@ -255,6 +255,9 @@ sc-power = 1 sc-r-power = 6 sc-sigma = 0.3 sc-coul = no +sc-scale-linpoint-LJ-gapsys = 0.85 +sc-scale-linpoint-Q-gapsys = 0.3 +sc-sigma-LJ-gapsys = 0.3 separate-dhdl-file = yes dhdl-derivatives = yes dh_hist_size = 0 diff --git a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsEmptyLines.xml b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsEmptyLines.xml index 343af36cad..5ca2760ed3 100644 --- a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsEmptyLines.xml +++ b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsEmptyLines.xml @@ -255,6 +255,9 @@ sc-power = 1 sc-r-power = 6 sc-sigma = 0.3 sc-coul = no +sc-scale-linpoint-LJ-gapsys = 0.85 +sc-scale-linpoint-Q-gapsys = 0.3 +sc-sigma-LJ-gapsys = 0.3 separate-dhdl-file = yes dhdl-derivatives = yes dh_hist_size = 0 diff --git a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsImplicitSolventNo.xml b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsImplicitSolventNo.xml index 343af36cad..5ca2760ed3 100644 --- a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsImplicitSolventNo.xml +++ b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsImplicitSolventNo.xml @@ -255,6 +255,9 @@ sc-power = 1 sc-r-power = 6 sc-sigma = 0.3 sc-coul = no +sc-scale-linpoint-LJ-gapsys = 0.85 +sc-scale-linpoint-Q-gapsys = 0.3 +sc-sigma-LJ-gapsys = 0.3 separate-dhdl-file = yes dhdl-derivatives = yes dh_hist_size = 0 diff --git a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsKeyWithoutValue.xml b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsKeyWithoutValue.xml index 343af36cad..5ca2760ed3 100644 --- a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsKeyWithoutValue.xml +++ b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsKeyWithoutValue.xml @@ -255,6 +255,9 @@ sc-power = 1 sc-r-power = 6 sc-sigma = 0.3 sc-coul = no +sc-scale-linpoint-LJ-gapsys = 0.85 +sc-scale-linpoint-Q-gapsys = 0.3 +sc-sigma-LJ-gapsys = 0.3 separate-dhdl-file = yes dhdl-derivatives = yes dh_hist_size = 0 diff --git a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsMimic.xml b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsMimic.xml index 318e2b498f..d8cff3f545 100644 --- a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsMimic.xml +++ b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsMimic.xml @@ -255,6 +255,9 @@ sc-power = 1 sc-r-power = 6 sc-sigma = 0.3 sc-coul = no +sc-scale-linpoint-LJ-gapsys = 0.85 +sc-scale-linpoint-Q-gapsys = 0.3 +sc-sigma-LJ-gapsys = 0.3 separate-dhdl-file = yes dhdl-derivatives = yes dh_hist_size = 0 diff --git a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_HandlesDifferentKindsOfMdpLines.xml b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_HandlesDifferentKindsOfMdpLines.xml index 1e25081dd8..c7b68c0a8b 100644 --- a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_HandlesDifferentKindsOfMdpLines.xml +++ b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_HandlesDifferentKindsOfMdpLines.xml @@ -255,6 +255,9 @@ sc-power = 1 sc-r-power = 6 sc-sigma = 0.3 sc-coul = no +sc-scale-linpoint-LJ-gapsys = 0.85 +sc-scale-linpoint-Q-gapsys = 0.3 +sc-sigma-LJ-gapsys = 0.3 separate-dhdl-file = yes dhdl-derivatives = yes dh_hist_size = 0 diff --git a/src/gromacs/listed_forces/pairs.cpp b/src/gromacs/listed_forces/pairs.cpp index e8e4a02ced..c97718459e 100644 --- a/src/gromacs/listed_forces/pairs.cpp +++ b/src/gromacs/listed_forces/pairs.cpp @@ -201,6 +201,7 @@ static real free_energy_evaluate_single(real real fscal_vdw[2], fscal_elec[2]; real velec[2], vvdw[2]; real dvdl_elec[2], dvdl_vdw[2]; + real scaleLinpointCoulGapsys, scaleLinpointVdWGapsys, sigma6VdWGapsys[2]; real rQ, rLJ; real scaleDvdlRCoul; int i, ntab; @@ -222,7 +223,7 @@ static real free_energy_evaluate_single(real /* Loop over state A(0) and B(1) */ for (i = 0; i < 2; i++) { - if (softcoreType == KernelSoftcoreType::Beutler || softcoreType == KernelSoftcoreType::Gapsys) + if constexpr (softcoreType == KernelSoftcoreType::Beutler) { if ((c6[i] > 0) && (c12[i] > 0)) { @@ -230,7 +231,7 @@ static real free_energy_evaluate_single(real * 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 (softcoreType == KernelSoftcoreType::Beutler && (sigma6[i] < scParams.sigma6Minimum)) /* 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; } @@ -241,9 +242,28 @@ static real free_energy_evaluate_single(real } sigma_pow[i] = sigma6[i]; } + 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. + */ + sigma6VdWGapsys[i] = half * c12[i] / c6[i]; + if (sigma6VdWGapsys[i] + < scParams.sigma6VdWGapsys) /* for disappearing coul and vdw with soft core at the same time */ + { + sigma6VdWGapsys[i] = scParams.sigma6VdWGapsys; + } + } + else + { + sigma6VdWGapsys[i] = scParams.sigma6VdWGapsys; + } + } } - if (softcoreType == KernelSoftcoreType::Beutler || softcoreType == KernelSoftcoreType::Gapsys) + 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)) @@ -253,16 +273,22 @@ static real free_energy_evaluate_single(real } else { - if (softcoreType == KernelSoftcoreType::Beutler) - { - alpha_vdw_eff = scParams.alphaVdw; - alpha_coul_eff = scParams.alphaCoulomb; - } - else if (softcoreType == KernelSoftcoreType::Gapsys) - { - alpha_vdw_eff = scParams.alphaVdw; - alpha_coul_eff = gmx::sixthroot(scParams.sigma6WithInvalidSigma); - } + alpha_vdw_eff = scParams.alphaVdw; + alpha_coul_eff = scParams.alphaCoulomb; + } + } + 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)) + { + scaleLinpointVdWGapsys = 0; + scaleLinpointCoulGapsys = 0; + } + else + { + scaleLinpointVdWGapsys = scParams.scaleLinpointVdWGapsys; + scaleLinpointCoulGapsys = scParams.scaleLinpointCoulGapsys; } } @@ -282,7 +308,7 @@ static real free_energy_evaluate_single(real if ((qq[i] != 0) || (c6[i] != 0) || (c12[i] != 0)) { /* Coulomb */ - if (softcoreType == KernelSoftcoreType::Beutler) + if constexpr (softcoreType == KernelSoftcoreType::Beutler) { rpinv = one / (alpha_coul_eff * lfac_coul[i] * sigma_pow[i] + rp); r_coul = sixthRoot(rpinv); @@ -293,10 +319,10 @@ static real free_energy_evaluate_single(real r_coul = r; } - if (softcoreType == KernelSoftcoreType::Gapsys) + if constexpr (softcoreType == KernelSoftcoreType::Gapsys) { rQ = gmx::sixthroot(one - LFC[i]) * (one + std::fabs(qq[i] / facel)); - rQ *= alpha_coul_eff; + rQ *= scaleLinpointCoulGapsys; scaleDvdlRCoul = 1; if (rQ > rCoulCutoff) { @@ -342,7 +368,7 @@ static real free_energy_evaluate_single(real } /* Vdw */ - if (softcoreType == KernelSoftcoreType::Beutler) + if constexpr (softcoreType == KernelSoftcoreType::Beutler) { rpinv = one / (alpha_vdw_eff * lfac_vdw[i] * sigma_pow[i] + rp); r_vdw = sixthRoot(rpinv); @@ -353,12 +379,12 @@ static real free_energy_evaluate_single(real r_vdw = r; } - if (softcoreType == KernelSoftcoreType::Gapsys) + if constexpr (softcoreType == KernelSoftcoreType::Gapsys) { constexpr real c_twentySixSeventh = 26.0_real / 7.0_real; - rLJ = gmx::sixthroot(c_twentySixSeventh * sigma6[i] * (one - LFV[i])); - rLJ *= alpha_vdw_eff; + rLJ = gmx::sixthroot(c_twentySixSeventh * sigma6VdWGapsys[i] * (one - LFV[i])); + rLJ *= scaleLinpointVdWGapsys; } if ((softcoreType == KernelSoftcoreType::Gapsys) && (r < rLJ)) @@ -446,14 +472,14 @@ static real free_energy_evaluate_single(real fscal += (LFC[i] * fscal_elec[i] + LFV[i] * fscal_vdw[i]) * rpm2; - if (softcoreType == KernelSoftcoreType::Gapsys) + if constexpr (softcoreType == KernelSoftcoreType::Gapsys) { dvdl_coul_sum += dvdl_elec[i]; dvdl_vdw_sum += dvdl_vdw[i]; } dvdl_coul_sum += velec[i] * DLF[i]; dvdl_vdw_sum += vvdw[i] * DLF[i]; - if (softcoreType == KernelSoftcoreType::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]; @@ -638,86 +664,119 @@ static real do_pairs_general(int ftype, c12B = iparams[itype].lj14.c12B * 12.0; const auto& scParams = *fr->ic->softCoreParameters; - if (scParams.softcoreType == SoftcoreType::Gapsys) - { - fscal = free_energy_evaluate_single( - 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 (scParams.softcoreType == SoftcoreType::Beutler) + if (scParams.softcoreType == SoftcoreType::Beutler) { - fscal = free_energy_evaluate_single( - 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( + 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( + 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 + else // Gapsys { - fscal = free_energy_evaluate_single( - 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.scaleLinpointCoulGapsys == 0 && scParams.scaleLinpointVdWGapsys == 0) + { + fscal = free_energy_evaluate_single( + 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( + 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 diff --git a/src/gromacs/mdtypes/inputrec.cpp b/src/gromacs/mdtypes/inputrec.cpp index 87f03085e8..97cd3ab741 100644 --- a/src/gromacs/mdtypes/inputrec.cpp +++ b/src/gromacs/mdtypes/inputrec.cpp @@ -601,6 +601,9 @@ static void pr_fepvals(FILE* fp, int indent, const t_lambda* fep, gmx_bool bMDPf PS("separate-dhdl-file", enumValueToString(fep->separate_dhdl_file)); PS("dhdl-derivatives", enumValueToString(fep->dhdl_derivatives)); PS("sc-function", enumValueToString(fep->softcoreFunction)); + PR("sc-scale-linpoint-LJ-gapsys", fep->scScaleLinpointLJGapsys); + PR("sc-scale-linpoint-Q-gapsys", fep->scScaleLinpointQGapsys); + PR("sc-sigma-LJ-gapsys", fep->scSigmaLJGapsys); }; static void pr_pull(FILE* fp, int indent, const pull_params_t& pull) @@ -1347,9 +1350,24 @@ static void cmp_fepvals(FILE* fp, const t_lambda* fep1, const t_lambda* fep2, re cmp_bool(fp, "inputrec->fepvals->bScCoul", -1, fep1->bScCoul, fep1->bScCoul); cmpEnum(fp, "inputrec->separate_dhdl_file", fep1->separate_dhdl_file, fep2->separate_dhdl_file); cmpEnum(fp, "inputrec->dhdl_derivatives", fep1->dhdl_derivatives, fep2->dhdl_derivatives); - cmpEnum(fp, "inputrec->softcoreFunction", fep1->softcoreFunction, fep2->softcoreFunction); cmp_int(fp, "inputrec->dh_hist_size", -1, fep1->dh_hist_size, fep2->dh_hist_size); cmp_double(fp, "inputrec->dh_hist_spacing", -1, fep1->dh_hist_spacing, fep2->dh_hist_spacing, ftol, abstol); + cmpEnum(fp, "inputrec->fepvals->softcoreFunction", fep1->softcoreFunction, fep2->softcoreFunction); + cmp_real(fp, + "inputrec->fepvals->scScaleLinpointLJGapsys", + -1, + fep1->scScaleLinpointLJGapsys, + fep2->scScaleLinpointLJGapsys, + ftol, + abstol); + cmp_real(fp, + "inputrec->fepvals->scScaleLinpointQGapsys", + -1, + fep1->scScaleLinpointQGapsys, + fep2->scScaleLinpointQGapsys, + ftol, + abstol); + cmp_real(fp, "inputrec->fepvals->scSigmaLJGapsys", -1, fep1->scSigmaLJGapsys, fep2->scSigmaLJGapsys, ftol, abstol); } void cmp_inputrec(FILE* fp, const t_inputrec* ir1, const t_inputrec* ir2, real ftol, real abstol) diff --git a/src/gromacs/mdtypes/inputrec.h b/src/gromacs/mdtypes/inputrec.h index aa40d516f0..404d486821 100644 --- a/src/gromacs/mdtypes/inputrec.h +++ b/src/gromacs/mdtypes/inputrec.h @@ -145,6 +145,12 @@ struct t_lambda bool bScCoul = false; //! The specific soft-core function to use SoftcoreType softcoreFunction = SoftcoreType::Beutler; + //! scale for the linearization point for the vdw interaction with gapsys soft-core + real scScaleLinpointLJGapsys = 0.85; + //! scale for the linearization point for the coulomb interaction with gapsys soft-core + real scScaleLinpointQGapsys = 0.3; + //! lower bound for c12/c6 in gapsys soft-core + real scSigmaLJGapsys = 0.3; //! Whether to print the dvdl term associated with this term; if it is not specified as separate, it is lumped with the FEP term gmx::EnumerationArray separate_dvdl; //! Whether to write a separate dhdl.xvg file note: NOT a gmx_bool, but an enum diff --git a/src/gromacs/mdtypes/interaction_const.cpp b/src/gromacs/mdtypes/interaction_const.cpp index 7463d3222b..0db07b5eb9 100644 --- a/src/gromacs/mdtypes/interaction_const.cpp +++ b/src/gromacs/mdtypes/interaction_const.cpp @@ -54,7 +54,10 @@ interaction_const_t::SoftCoreParameters::SoftCoreParameters(const t_lambda& fepv lambdaPower(fepvals.sc_power), sigma6WithInvalidSigma(gmx::power6(fepvals.sc_sigma)), sigma6Minimum(fepvals.bScCoul ? gmx::power6(fepvals.sc_sigma_min) : 0), - softcoreType(fepvals.softcoreFunction) + softcoreType(fepvals.softcoreFunction), + scaleLinpointVdWGapsys(fepvals.scScaleLinpointLJGapsys), + scaleLinpointCoulGapsys(fepvals.scScaleLinpointQGapsys), + sigma6VdWGapsys(gmx::power6(fepvals.scSigmaLJGapsys)) { // This is checked during tpr reading, so we can assert here GMX_RELEASE_ASSERT(fepvals.sc_r_power == 6.0, "We only support soft-core r-power 6"); diff --git a/src/gromacs/mdtypes/interaction_const.h b/src/gromacs/mdtypes/interaction_const.h index 29c3255011..035cda1dbe 100644 --- a/src/gromacs/mdtypes/interaction_const.h +++ b/src/gromacs/mdtypes/interaction_const.h @@ -129,6 +129,12 @@ struct interaction_const_t real sigma6Minimum; // soft-core function SoftcoreType softcoreType; + // (gapsys sc) linearization point scaling for vdW interactions + real scaleLinpointVdWGapsys; + // (gapsys sc) linearization point scaling for Coulomb interactions + real scaleLinpointCoulGapsys; + // (gapsys sc) lower bound/replacement for c12/c6 in vdw interactions + real sigma6VdWGapsys; }; /* VdW */ -- 2.22.0