Guard for epsfac equals zero in sc-gapsys.

diff --git a/src/gromacs/gmxlib/nonbonded/nb_softcore.h b/src/gromacs/gmxlib/nonbonded/nb_softcore.h
index f06c48b3fc7a68775ec192aa96ba45b6055b83a1..a2615ede4a7aedd78199d58b71bf94aa4de1bbda 100644 (file)
--- a/src/gromacs/gmxlib/nonbonded/nb_softcore.h
+++ b/src/gromacs/gmxlib/nonbonded/nb_softcore.h
@@ -84,7 +84,7 @@ static inline void reactionFieldQuadraticPotential(const RealType qq,
                                                    BoolType       mask)
 {
     /* check if we have to use the hardcore values */
-    BoolType computeValues = mask && (lambdaFac < 1 && 0 < alphaEff);
+    BoolType computeValues = mask && (lambdaFac < 1 && 0 < alphaEff && facel != 0);
     if (gmx::anyTrue(computeValues))
     {
         RealType lambdaFacRev = gmx::selectByMask(1 - lambdaFac, computeValues);
@@ -140,9 +140,8 @@ static inline void ewaldQuadraticPotential(const RealType qq,
                                            RealType*      dvdl,
                                            BoolType       mask)
 {
-
     /* check if we have to use the hardcore values */
-    BoolType computeValues = mask && (lambdaFac < 1 && 0 < alphaEff);
+    BoolType computeValues = mask && (lambdaFac < 1 && 0 < alphaEff && facel != 0);
     if (gmx::anyTrue(computeValues))
     {
         RealType lambdaFacRev = gmx::selectByMask(1 - lambdaFac, computeValues);

diff --git a/src/gromacs/listed_forces/pairs.cpp b/src/gromacs/listed_forces/pairs.cpp
index 2d11a3cf71eafeffdefe708cde6d765e00bf43a2..12b22f51c394f1c6d027ca1199bd3e0a714ac66c 100644 (file)
--- a/src/gromacs/listed_forces/pairs.cpp
+++ b/src/gromacs/listed_forces/pairs.cpp
@@ -316,8 +316,15 @@ static real free_energy_evaluate_single(real
 
             if constexpr (softcoreType == KernelSoftcoreType::Gapsys)
             {
-                rQ = gmx::sixthroot(one - LFC[i]) * (one + std::fabs(qq[i] / facel));
-                rQ *= scaleLinpointCoulGapsys;
+                if ((facel != 0) && (LFC[i] < 1))
+                {
+                    rQ = gmx::sixthroot(one - LFC[i]) * (one + std::fabs(qq[i] / facel));
+                    rQ *= scaleLinpointCoulGapsys;
+                }
+                else
+                {
+                    rQ = 0;
+                }
                 scaleDvdlRCoul = 1;
                 if (rQ > rCoulCutoff)
                 {
@@ -378,8 +385,16 @@ static real free_energy_evaluate_single(real
             {
                 constexpr real c_twentySixSeventh = 26.0_real / 7.0_real;
 
-                rLJ = gmx::sixthroot(c_twentySixSeventh * sigma6VdWGapsys[i] * (one - LFV[i]));
-                rLJ *= scaleLinpointVdWGapsys;
+                if (LFV[i] < 1)
+                {
+
+                    rLJ = gmx::sixthroot(c_twentySixSeventh * sigma6VdWGapsys[i] * (one - LFV[i]));
+                    rLJ *= scaleLinpointVdWGapsys;
+                }
+                else
+                {
+                    rLJ = 0;
+                }
             }
 
             if ((softcoreType == KernelSoftcoreType::Gapsys) && (r < rLJ))