Introduce explicit parameters for gapsys-sc.

diff --git a/src/gromacs/fileio/tpxio.cpp b/src/gromacs/fileio/tpxio.cpp
index 33e63d4f11b5029b2fb7b6f4dc38fb448af4e6a9..39943911faed17b4f26efa51f05c207ed1b6b601 100644 (file)
--- 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<int*>(&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 33e485d191f139b7152b1d28cf80e66f5d1520b5..41bcedc0585e0de8d7559f0e9d256797d501792b 100644 (file)
--- 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<const int> shift  = nlist.shift;
     gmx::ArrayRef<const int> gid    = nlist.gid;
 
-    const real  lambda_coul = lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Coul)];
-    const real  lambda_vdw  = lambda[static_cast<int>(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<int>(FreeEnergyPerturbationCouplingType::Coul)];
+    const real lambda_vdw  = lambda[static_cast<int>(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<RealType>(preloadQq[i]);
                 ljPmeC6Grid[i] = gmx::load<RealType>(preloadLjPmeC6Grid[i]);
-                if constexpr (softcoreType == KernelSoftcoreType::Beutler
-                              || softcoreType == KernelSoftcoreType::Gapsys)
+                if constexpr (softcoreType == KernelSoftcoreType::Beutler)
                 {
                     sigma6[i] = gmx::load<RealType>(preloadSigma6[i]);
                 }
+                if constexpr (softcoreType == KernelSoftcoreType::Gapsys)
+                {
+                    sigmaVdWGapsysEff[i] = gmx::load<RealType>(preloadSigma6VdWGapsys[i]);
+                }
             }
-            if constexpr (softcoreType == KernelSoftcoreType::Beutler
-                          || softcoreType == KernelSoftcoreType::Gapsys)
+            if constexpr (softcoreType == KernelSoftcoreType::Beutler)
             {
                 alphaVdwEff  = gmx::load<RealType>(preloadAlphaVdwEff);
                 alphaCoulEff = gmx::load<RealType>(preloadAlphaCoulEff);
             }
+            if constexpr (softcoreType == KernelSoftcoreType::Gapsys)
+            {
+                scaleLinpointVdWGapsysEff  = gmx::load<RealType>(preloadScaleLinpointVdWGapsys);
+                scaleLinpointCoulGapsysEff = gmx::load<RealType>(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<KernelSoftcoreType::None>(
-                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<KernelSoftcoreType::Beutler>(
+            return (dispatchKernelOnScLambdasOrAlphasDifference<KernelSoftcoreType::None>(
                     scLambdasOrAlphasDiffer,
                     vdwInteractionTypeIsEwald,
                     elecInteractionTypeIsEwald,
@@ -1264,9 +1308,19 @@ static KernelFunction dispatchKernel(const bool                 scLambdasOrAlpha
                     computeForces,
                     useSimd));
         }
-        else
+        return (dispatchKernelOnScLambdasOrAlphasDifference<KernelSoftcoreType::Beutler>(
+                scLambdasOrAlphasDiffer,
+                vdwInteractionTypeIsEwald,
+                elecInteractionTypeIsEwald,
+                vdwModifierIsPotSwitch,
+                computeForces,
+                useSimd));
+    }
+    else // Gapsys
+    {
+        if (scParams.scaleLinpointCoulGapsys == 0 && scParams.scaleLinpointVdWGapsys == 0)
         {
-            return (dispatchKernelOnScLambdasOrAlphasDifference<KernelSoftcoreType::Gapsys>(
+            return (dispatchKernelOnScLambdasOrAlphasDifference<KernelSoftcoreType::None>(
                     scLambdasOrAlphasDiffer,
                     vdwInteractionTypeIsEwald,
                     elecInteractionTypeIsEwald,
@@ -1274,6 +1328,13 @@ static KernelFunction dispatchKernel(const bool                 scLambdasOrAlpha
                     computeForces,
                     useSimd));
         }
+        return (dispatchKernelOnScLambdasOrAlphasDifference<KernelSoftcoreType::Gapsys>(
+                scLambdasOrAlphasDiffer,
+                vdwInteractionTypeIsEwald,
+                elecInteractionTypeIsEwald,
+                vdwModifierIsPotSwitch,
+                computeForces,
+                useSimd));
     }
 }
 

diff --git a/src/gromacs/gmxpreprocess/readir.cpp b/src/gromacs/gmxpreprocess/readir.cpp
index 99e02def68acd050bc6241b535c7534e8a01bead..c5873ab705eb2e5ba15787c03aa2cb3f0e925546 100644 (file)
--- 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<FreeEnergyPrintEnergy>(&inp, "dhdl-print-energy", wi);
-    fep->softcoreFunction   = getEnum<SoftcoreType>(&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<Boolean>(&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<SeparateDhdlFile>(&inp, "separate-dhdl-file", wi);
-    fep->dhdl_derivatives   = getEnum<DhDlDerivativeCalculation>(&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<FreeEnergyPrintEnergy>(&inp, "dhdl-print-energy", wi);
+    fep->softcoreFunction        = getEnum<SoftcoreType>(&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<Boolean>(&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<SeparateDhdlFile>(&inp, "separate-dhdl-file", wi);
+    fep->dhdl_derivatives        = getEnum<DhDlDerivativeCalculation>(&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 1c2029b4f2aea7a4d16661e8da06f890a9e07acc..ed5a7eb1fce3f5a9f514ce6c92e74320f5c5b5d3 100644 (file)
--- 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 21fe54ad8a10ca0fcfa32d62ebe502c3d7c1b518..34dea51cca70bce4d72ea8c12a0cc963b2c9cc39 100644 (file)
--- 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 e8c72614b2a58713974f759b5f3c242686e05039..9fc0ed9f5077eab6418b791773af75b9f9e1b094 100644 (file)
--- 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 96b74dfec237da7d6bfbfb963ba2f0b71b9af697..902a7b45fe7d5f40914d1e2b023d84224cf24fcb 100644 (file)
--- 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 343af36cad09a179b24be8b4f2e3de5bc3db22bc..5ca2760ed362ad5cd16ca2434ee2bd4e2f39d9c0 100644 (file)
--- 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 343af36cad09a179b24be8b4f2e3de5bc3db22bc..5ca2760ed362ad5cd16ca2434ee2bd4e2f39d9c0 100644 (file)
--- 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 343af36cad09a179b24be8b4f2e3de5bc3db22bc..5ca2760ed362ad5cd16ca2434ee2bd4e2f39d9c0 100644 (file)
--- 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 318e2b498f7c1eeba84a19248c13ff8ba1f7bbed..d8cff3f54527902429175e5303e882a8ac224d58 100644 (file)
--- 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 1e25081dd86dfc6189f0ccbfe09f19df28c9c6eb..c7b68c0a8bda561f02d3ee62c3a2df5397689232 100644 (file)
--- 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 e8e4a02ced7f4d5d1b9b6f3c968a427dd9341703..c97718459e7bab513f94d8dcb7811e1437c5d406 100644 (file)
--- 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<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 if (scParams.softcoreType == SoftcoreType::Beutler)
+            if (scParams.softcoreType == SoftcoreType::Beutler)
             {
-                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);
+                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
+            else // Gapsys
             {
-                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);
+                if (scParams.scaleLinpointCoulGapsys == 0 && scParams.scaleLinpointVdWGapsys == 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

diff --git a/src/gromacs/mdtypes/inputrec.cpp b/src/gromacs/mdtypes/inputrec.cpp
index 87f03085e8cc40df25e5676eea0ead8b902972d5..97cd3ab7411f315c2048a23bfbb35395eefe6994 100644 (file)
--- 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 aa40d516f0845866e82303c5193238e30a1ddaca..404d486821de04b0195bae87c7f834cd1d7d12ae 100644 (file)
--- 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<FreeEnergyPerturbationCouplingType, bool> 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 7463d3222bb130c2ecd74e053f23dab2dac9aea1..0db07b5eb929659f07b559f7620910ba3811e650 100644 (file)
--- 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 29c32550111184bc98247f3478c4e9c952531c4a..035cda1dbe8371ec5f75a182abeb390954f1a1af 100644 (file)
--- 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 */