Fix use of simd-masks in gapsys softcore.

[alexxy/gromacs.git] / src / gromacs / gmxlib / nonbonded / nb_softcore.h

diff --git a/src/gromacs/gmxlib/nonbonded/nb_softcore.h b/src/gromacs/gmxlib/nonbonded/nb_softcore.h
index a2615ede4a7aedd78199d58b71bf94aa4de1bbda..de1852d3b718014fe904d5d09d70f02f822d4bcc 100644 (file)
--- a/src/gromacs/gmxlib/nonbonded/nb_softcore.h
+++ b/src/gromacs/gmxlib/nonbonded/nb_softcore.h
@@ -43,23 +43,26 @@
 #include "gromacs/simd/simd_math.h"
 
 /* linearized electrostatics */
-template<class RealType, class BoolType>
+template<bool computeForces, class RealType, class BoolType>
 static inline void quadraticApproximationCoulomb(const RealType qq,
                                                  const RealType rInvQ,
                                                  const RealType r,
                                                  const real     lambdaFac,
                                                  const real     dLambdaFac,
-                                                 RealType*      force,
-                                                 RealType*      potential,
-                                                 RealType*      dvdl,
-                                                 BoolType       dvdlMask)
+                                                 RealType gmx_unused* force,
+                                                 RealType*            potential,
+                                                 RealType*            dvdl,
+                                                 BoolType             dvdlMask)
 {
     RealType constFac = qq * rInvQ;
     RealType linFac   = constFac * r * rInvQ;
     RealType quadrFac = linFac * r * rInvQ;
 
     /* Computing Coulomb force and potential energy */
-    *force = -2 * quadrFac + 3 * linFac;
+    if constexpr (computeForces)
+    {
+        *force = -2 * quadrFac + 3 * linFac;
+    }
 
     *potential = quadrFac - 3 * (linFac - constFac);
 
@@ -68,7 +71,7 @@ static inline void quadraticApproximationCoulomb(const RealType qq,
 }
 
 /* reaction-field linearized electrostatics */
-template<class RealType, class BoolType>
+template<bool computeForces, class RealType, class BoolType>
 static inline void reactionFieldQuadraticPotential(const RealType qq,
                                                    const real     facel,
                                                    const RealType r,
@@ -78,13 +81,16 @@ static inline void reactionFieldQuadraticPotential(const RealType qq,
                                                    const RealType alphaEff,
                                                    const real     krf,
                                                    const real     potentialShift,
-                                                   RealType*      force,
-                                                   RealType*      potential,
-                                                   RealType*      dvdl,
-                                                   BoolType       mask)
+                                                   RealType gmx_unused* force,
+                                                   RealType*            potential,
+                                                   RealType*            dvdl,
+                                                   BoolType             mask)
 {
+    RealType one(1);
+    RealType zero(0);
+
     /* check if we have to use the hardcore values */
-    BoolType computeValues = mask && (lambdaFac < 1 && 0 < alphaEff && facel != 0);
+    BoolType computeValues = mask && (lambdaFac < one && zero < alphaEff && facel != zero);
     if (gmx::anyTrue(computeValues))
     {
         RealType lambdaFacRev = gmx::selectByMask(1 - lambdaFac, computeValues);
@@ -110,7 +116,7 @@ static inline void reactionFieldQuadraticPotential(const RealType qq,
             RealType forceQuad(0);
             RealType potentialQuad(0);
             RealType dvdlQuad(0);
-            quadraticApproximationCoulomb(
+            quadraticApproximationCoulomb<computeForces>(
                     qq, rInvQ, r, lambdaFac, dLambdaFac, &forceQuad, &potentialQuad, &dvdlQuad, computeValues);
 
             // rf modification
@@ -118,7 +124,10 @@ static inline void reactionFieldQuadraticPotential(const RealType qq,
             potentialQuad = potentialQuad + qq * (krf * r * r - potentialShift);
 
             // update
-            *force     = gmx::blend(*force, forceQuad, computeValues);
+            if constexpr (computeForces)
+            {
+                *force = gmx::blend(*force, forceQuad, computeValues);
+            }
             *potential = gmx::blend(*potential, potentialQuad, computeValues);
             *dvdl      = *dvdl + gmx::selectByMask(dvdlQuad, computeValues && withinCutoff);
         }
@@ -126,7 +135,7 @@ static inline void reactionFieldQuadraticPotential(const RealType qq,
 }
 
 /* ewald linearized electrostatics */
-template<class RealType, class BoolType>
+template<bool computeForces, class RealType, class BoolType>
 static inline void ewaldQuadraticPotential(const RealType qq,
                                            const real     facel,
                                            const RealType r,
@@ -135,13 +144,16 @@ static inline void ewaldQuadraticPotential(const RealType qq,
                                            const real     dLambdaFac,
                                            const RealType alphaEff,
                                            const real     potentialShift,
-                                           RealType*      force,
-                                           RealType*      potential,
-                                           RealType*      dvdl,
-                                           BoolType       mask)
+                                           RealType gmx_unused* force,
+                                           RealType*            potential,
+                                           RealType*            dvdl,
+                                           BoolType             mask)
 {
+    RealType one(1);
+    RealType zero(0);
+
     /* check if we have to use the hardcore values */
-    BoolType computeValues = mask && (lambdaFac < 1 && 0 < alphaEff && facel != 0);
+    BoolType computeValues = mask && (lambdaFac < one && zero < alphaEff && facel != zero);
     if (gmx::anyTrue(computeValues))
     {
         RealType lambdaFacRev = gmx::selectByMask(1 - lambdaFac, computeValues);
@@ -167,14 +179,17 @@ static inline void ewaldQuadraticPotential(const RealType qq,
             RealType forceQuad(0);
             RealType potentialQuad(0);
             RealType dvdlQuad(0);
-            quadraticApproximationCoulomb(
+            quadraticApproximationCoulomb<computeForces>(
                     qq, rInvQ, r, lambdaFac, dLambdaFac, &forceQuad, &potentialQuad, &dvdlQuad, computeValues);
 
             // ewald modification
             potentialQuad = potentialQuad - qq * potentialShift;
 
             // update
-            *force     = gmx::blend(*force, forceQuad, computeValues);
+            if constexpr (computeForces)
+            {
+                *force = gmx::blend(*force, forceQuad, computeValues);
+            }
             *potential = gmx::blend(*potential, potentialQuad, computeValues);
             *dvdl      = *dvdl + gmx::selectByMask(dvdlQuad, computeValues && withinCutoff);
         }
@@ -182,7 +197,7 @@ static inline void ewaldQuadraticPotential(const RealType qq,
 }
 
 /* cutoff LJ with quadratic appximation of lj-potential */
-template<class RealType, class BoolType>
+template<bool computeForces, class RealType, class BoolType>
 static inline void lennardJonesQuadraticPotential(const RealType c6,
                                                   const RealType c12,
                                                   const RealType r,
@@ -193,18 +208,21 @@ static inline void lennardJonesQuadraticPotential(const RealType c6,
                                                   const RealType alphaEff,
                                                   const real     repulsionShift,
                                                   const real     dispersionShift,
-                                                  RealType*      force,
-                                                  RealType*      potential,
-                                                  RealType*      dvdl,
-                                                  BoolType       mask)
+                                                  RealType gmx_unused* force,
+                                                  RealType*            potential,
+                                                  RealType*            dvdl,
+                                                  BoolType             mask)
 {
     constexpr real c_twentySixSeventh = 26.0_real / 7.0_real;
     constexpr real c_oneSixth         = 1.0_real / 6.0_real;
     constexpr real c_oneTwelth        = 1.0_real / 12.0_real;
     constexpr real c_half             = 1.0_real / 2.0_real;
 
+    RealType one(1);
+    RealType zero(0);
+
     /* check if we have to use the hardcore values */
-    BoolType computeValues = mask && (lambdaFac < 1 && 0 < alphaEff);
+    BoolType computeValues = mask && (lambdaFac < one && zero < alphaEff);
     if (gmx::anyTrue(computeValues))
     {
         RealType lambdaFacRev    = gmx::selectByMask(1 - lambdaFac, computeValues);
@@ -243,16 +261,19 @@ static inline void lennardJonesQuadraticPotential(const RealType c6,
             constFac  = 91 * rInv12C - 28 * rInv6C;
 
             /* Computing LJ force and potential energy */
-            RealType forceQuad     = -quadrFac + linearFac;
-            RealType potentialQuad = c_half * quadrFac - linearFac + constFac;
-            RealType dvdlQuad      = dLambdaFac * 28 * (lambdaFac * lambdaFacRevInv)
+            RealType gmx_unused forceQuad     = -quadrFac + linearFac;
+            RealType            potentialQuad = c_half * quadrFac - linearFac + constFac;
+            RealType            dvdlQuad      = dLambdaFac * 28 * (lambdaFac * lambdaFacRevInv)
                                 * ((6.5_real * rInv14C - rInv8C) - (13 * rInv13C - 2 * rInv7C)
                                    + (6.5_real * rInv12C - rInv6C));
 
             potentialQuad = potentialQuad
                             + gmx::selectByMask(((c12s * repulsionShift) - (c6s * dispersionShift)),
                                                 computeValues);
-            *force     = gmx::blend(*force, forceQuad, computeValues);
+            if constexpr (computeForces)
+            {
+                *force = gmx::blend(*force, forceQuad, computeValues);
+            }
             *potential = gmx::blend(*potential, potentialQuad, computeValues);
             *dvdl      = *dvdl + gmx::selectByMask(dvdlQuad, computeValues);
         }