#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);
- RealType lambdaFacRevInv = gmx::maskzInv(1.0 - lambdaFac, dvdlMask);
+ RealType lambdaFacRevInv = gmx::maskzInv(1 - lambdaFac, dvdlMask);
*dvdl = dLambdaFac * 0.5_real * (lambdaFac * lambdaFacRevInv) * (quadrFac - 2 * linFac + constFac);
}
/* 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,
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);
+ BoolType computeValues = mask && (lambdaFac < one && zero < alphaEff && facel != zero);
if (gmx::anyTrue(computeValues))
{
RealType lambdaFacRev = gmx::selectByMask(1 - lambdaFac, computeValues);
RealType forceQuad(0);
RealType potentialQuad(0);
RealType dvdlQuad(0);
- quadraticApproximationCoulomb(
+ quadraticApproximationCoulomb<computeForces>(
qq, rInvQ, r, lambdaFac, dLambdaFac, &forceQuad, &potentialQuad, &dvdlQuad, computeValues);
// rf modification
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);
}
}
/* 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,
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);
+ BoolType computeValues = mask && (lambdaFac < one && zero < alphaEff && facel != zero);
if (gmx::anyTrue(computeValues))
{
RealType lambdaFacRev = gmx::selectByMask(1 - lambdaFac, computeValues);
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);
}
}
/* 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,
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);
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)
- * ((6.5_real * rInv14C - rInv8C) - (13 * rInv13C - 2. * rInv7C)
+ 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);
+ + gmx::selectByMask(((c12s * repulsionShift) - (c6s * dispersionShift)),
+ computeValues);
+ if constexpr (computeForces)
+ {
+ *force = gmx::blend(*force, forceQuad, computeValues);
+ }
*potential = gmx::blend(*potential, potentialQuad, computeValues);
*dvdl = *dvdl + gmx::selectByMask(dvdlQuad, computeValues);
}