From 24a1974c10922a655213bfe33790f2fd179f9aff Mon Sep 17 00:00:00 2001
From: Sebastian Kehl <sebastian.kehl@mpcdf.mpg.de>
Date: Wed, 29 Sep 2021 08:43:00 +0200
Subject: [PATCH] Guard for epsfac equals zero in sc-gapsys.

---
 src/gromacs/gmxlib/nonbonded/nb_softcore.h |  5 ++---
 src/gromacs/listed_forces/pairs.cpp        | 23 ++++++++++++++++++----
 2 files changed, 21 insertions(+), 7 deletions(-)

diff --git a/src/gromacs/gmxlib/nonbonded/nb_softcore.h b/src/gromacs/gmxlib/nonbonded/nb_softcore.h
index f06c48b3fc..a2615ede4a 100644
--- 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 2d11a3cf71..12b22f51c3 100644
--- 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))
-- 
2.22.0