Minor code clean-up in the NB FEP kernel.
[alexxy/gromacs.git] / src / gromacs / gmxlib / nonbonded / nb_free_energy.cpp
index c8e46d15d3a5cb381c876ddc79941fc6fabf9a39..214e41a9d93de0b01c29b8a0db279bbb4b6e5e79 100644 (file)
@@ -63,6 +63,7 @@
 #include "gromacs/utility/fatalerror.h"
 #include "gromacs/utility/arrayref.h"
 
+#include "nb_softcore.h"
 
 //! Scalar (non-SIMD) data types.
 struct ScalarDataTypes
@@ -257,7 +258,7 @@ static inline RealType potSwitchPotentialMod(const RealType potentialInp, const
 
 
 //! Templated free-energy non-bonded kernel
-template<typename DataTypes, bool useSoftCore, bool scLambdasOrAlphasDiffer, bool vdwInteractionTypeIsEwald, bool elecInteractionTypeIsEwald, bool vdwModifierIsPotSwitch, bool computeForces>
+template<typename DataTypes, KernelSoftcoreType softcoreType, bool scLambdasOrAlphasDiffer, bool vdwInteractionTypeIsEwald, bool elecInteractionTypeIsEwald, bool vdwModifierIsPotSwitch, bool computeForces>
 static void nb_free_energy_kernel(const t_nblist&                                  nlist,
                                   const gmx::ArrayRefWithPadding<const gmx::RVec>& coords,
                                   const int                                        ntype,
@@ -303,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 gapsysScaleLinpointCoul = scParams.gapsysScaleLinpointCoul;
+    const real gmx_unused gapsysScaleLinpointVdW  = scParams.gapsysScaleLinpointVdW;
+    const real gmx_unused gapsysSigma6VdW         = scParams.gapsysSigma6VdW;
+
     const bool gmx_unused doShiftForces = ((flags & GMX_NONBONDED_DO_SHIFTFORCE) != 0);
     const bool            doPotential   = ((flags & GMX_NONBONDED_DO_POTENTIAL) != 0);
 
@@ -364,8 +372,9 @@ static void nb_free_energy_kernel(const t_nblist&
         icoul = NbkernelElecType::None;
     }
 
-    real rcutoff_max2 = std::max(ic.rcoulomb, ic.rvdw);
-    rcutoff_max2      = rcutoff_max2 * rcutoff_max2;
+    real rcutoff_max2                 = std::max(ic.rcoulomb, ic.rvdw);
+    rcutoff_max2                      = rcutoff_max2 * rcutoff_max2;
+    const real gmx_unused rCutoffCoul = ic.rcoulomb;
 
     real gmx_unused sh_ewald = 0;
     if constexpr (elecInteractionTypeIsEwald || vdwInteractionTypeIsEwald)
@@ -392,8 +401,8 @@ static void nb_free_energy_kernel(const t_nblist&
     const RealType            maxRInvSix(c_maxRInvSix);
     const RealType gmx_unused floatMin(GMX_FLOAT_MIN);
 
-    RealType dvdlCoul(zero);
-    RealType dvdlVdw(zero);
+    RealType dvdlCoul{ zero };
+    RealType dvdlVdw{ zero };
 
     /* Lambda factor for state A, 1-lambda*/
     real LFC[NSTATES], LFV[NSTATES];
@@ -470,10 +479,10 @@ static void nb_free_energy_kernel(const t_nblist&
         int            preloadIi[DataTypes::simdRealWidth];
         int gmx_unused preloadIs[DataTypes::simdRealWidth];
 #endif
-        for (int s = 0; s < DataTypes::simdRealWidth; s++)
+        for (int i = 0; i < DataTypes::simdRealWidth; i++)
         {
-            preloadIi[s] = ii;
-            preloadIs[s] = shift[n];
+            preloadIi[i] = ii;
+            preloadIs[i] = shift[n];
         }
         IntType ii_s = gmx::load<IntType>(preloadIi);
 
@@ -490,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 preloadGapsysScaleLinpointVdW[DataTypes::simdRealWidth];
+            alignas(GMX_SIMD_ALIGNMENT)
+                    real gmx_unused preloadGapsysScaleLinpointCoul[DataTypes::simdRealWidth];
+            alignas(GMX_SIMD_ALIGNMENT)
+                    real gmx_unused preloadGapsysSigma6VdW[NSTATES][DataTypes::simdRealWidth];
             alignas(GMX_SIMD_ALIGNMENT) real preloadLjPmeC6Grid[NSTATES][DataTypes::simdRealWidth];
 #else
             real            preloadPairIsValid[DataTypes::simdRealWidth];
@@ -500,84 +515,120 @@ 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 preloadGapsysScaleLinpointVdW[DataTypes::simdRealWidth];
+            real gmx_unused preloadGapsysScaleLinpointCoul[DataTypes::simdRealWidth];
+            real gmx_unused preloadGapsysSigma6VdW[NSTATES][DataTypes::simdRealWidth];
             real            preloadLjPmeC6Grid[NSTATES][DataTypes::simdRealWidth];
 #endif
-            for (int s = 0; s < DataTypes::simdRealWidth; s++)
+            for (int j = 0; j < DataTypes::simdRealWidth; j++)
             {
-                if (k + s < nj1)
+                if (k + j < nj1)
                 {
-                    preloadPairIsValid[s] = true;
+                    preloadPairIsValid[j] = true;
                     /* Check if this pair on the exclusions list.*/
-                    preloadPairIncluded[s]  = (nlist.excl_fep.empty() || nlist.excl_fep[k + s]);
-                    const int jnr           = jjnr[k + s];
-                    preloadJnr[s]           = jnr;
-                    typeIndices[STATE_A][s] = ntiA + typeA[jnr];
-                    typeIndices[STATE_B][s] = ntiB + typeB[jnr];
-                    preloadQq[STATE_A][s]   = iqA * chargeA[jnr];
-                    preloadQq[STATE_B][s]   = iqB * chargeB[jnr];
+                    preloadPairIncluded[j]  = (nlist.excl_fep.empty() || nlist.excl_fep[k + j]);
+                    const int jnr           = jjnr[k + j];
+                    preloadJnr[j]           = jnr;
+                    typeIndices[STATE_A][j] = ntiA + typeA[jnr];
+                    typeIndices[STATE_B][j] = ntiB + typeB[jnr];
+                    preloadQq[STATE_A][j]   = iqA * chargeA[jnr];
+                    preloadQq[STATE_B][j]   = iqB * chargeB[jnr];
 
                     for (int i = 0; i < NSTATES; i++)
                     {
                         if constexpr (vdwInteractionTypeIsEwald)
                         {
-                            preloadLjPmeC6Grid[i][s] = nbfp_grid[2 * typeIndices[i][s]];
+                            preloadLjPmeC6Grid[i][j] = nbfp_grid[2 * typeIndices[i][j]];
                         }
                         else
                         {
-                            preloadLjPmeC6Grid[i][s] = 0;
+                            preloadLjPmeC6Grid[i][j] = 0;
                         }
-                        if constexpr (useSoftCore)
+                        if constexpr (softcoreType == KernelSoftcoreType::Beutler)
                         {
-                            const real c6  = nbfp[2 * typeIndices[i][s]];
-                            const real c12 = nbfp[2 * typeIndices[i][s] + 1];
+                            const real c6  = nbfp[2 * typeIndices[i][j]];
+                            const real c12 = nbfp[2 * typeIndices[i][j] + 1];
                             if (c6 > 0 && c12 > 0)
                             {
                                 /* c12 is stored scaled with 12.0 and c6 is scaled with 6.0 - correct for this */
-                                preloadSigma6[i][s] = 0.5_real * c12 / c6;
-                                if (preloadSigma6[i][s]
+                                preloadSigma6[i][j] = 0.5_real * c12 / c6;
+                                if (preloadSigma6[i][j]
                                     < sigma6_min) /* for disappearing coul and vdw with soft core at the same time */
                                 {
-                                    preloadSigma6[i][s] = sigma6_min;
+                                    preloadSigma6[i][j] = sigma6_min;
                                 }
                             }
                             else
                             {
-                                preloadSigma6[i][s] = sigma6_def;
+                                preloadSigma6[i][j] = sigma6_def;
                             }
                         }
+                        if constexpr (softcoreType == KernelSoftcoreType::Gapsys)
+                        {
+                            const real c6  = nbfp[2 * typeIndices[i][j]];
+                            const real c12 = nbfp[2 * typeIndices[i][j] + 1];
+                            if (c6 > 0 && c12 > 0)
+                            {
+                                /* c12 is stored scaled with 12.0 and c6 is scaled with 6.0 - correct for this */
+                                preloadGapsysSigma6VdW[i][j] = 0.5_real * c12 / c6;
+                            }
+                            else
+                            {
+                                preloadGapsysSigma6VdW[i][j] = gapsysSigma6VdW;
+                            }
+                        }
+                    }
+                    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][j] + 1];
+                        const real c12B = nbfp[2 * typeIndices[STATE_B][j] + 1];
+                        if (c12A > 0 && c12B > 0)
+                        {
+                            preloadAlphaVdwEff[j]  = 0;
+                            preloadAlphaCoulEff[j] = 0;
+                        }
+                        else
+                        {
+                            preloadAlphaVdwEff[j]  = alpha_vdw;
+                            preloadAlphaCoulEff[j] = alpha_coul;
+                        }
                     }
-                    if constexpr (useSoftCore)
+                    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];
+                        const real c12A = nbfp[2 * typeIndices[STATE_A][j] + 1];
+                        const real c12B = nbfp[2 * typeIndices[STATE_B][j] + 1];
                         if (c12A > 0 && c12B > 0)
                         {
-                            preloadAlphaVdwEff[s]  = 0;
-                            preloadAlphaCoulEff[s] = 0;
+                            preloadGapsysScaleLinpointVdW[j]  = 0;
+                            preloadGapsysScaleLinpointCoul[j] = 0;
                         }
                         else
                         {
-                            preloadAlphaVdwEff[s]  = alpha_vdw;
-                            preloadAlphaCoulEff[s] = alpha_coul;
+                            preloadGapsysScaleLinpointVdW[j]  = gapsysScaleLinpointVdW;
+                            preloadGapsysScaleLinpointCoul[j] = gapsysScaleLinpointCoul;
                         }
                     }
                 }
                 else
                 {
-                    preloadJnr[s]          = jjnr[k];
-                    preloadPairIsValid[s]  = false;
-                    preloadPairIncluded[s] = false;
-                    preloadAlphaVdwEff[s]  = 0;
-                    preloadAlphaCoulEff[s] = 0;
-
+                    preloadJnr[j]                     = jjnr[k];
+                    preloadPairIsValid[j]             = false;
+                    preloadPairIncluded[j]            = false;
+                    preloadAlphaVdwEff[j]             = 0;
+                    preloadAlphaCoulEff[j]            = 0;
+                    preloadGapsysScaleLinpointVdW[j]  = 0;
+                    preloadGapsysScaleLinpointCoul[j] = 0;
+
+                    typeIndices[STATE_A][j] = ntiA + typeA[jjnr[k]];
+                    typeIndices[STATE_B][j] = ntiB + typeB[jjnr[k]];
                     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;
+                        preloadLjPmeC6Grid[i][j]     = 0;
+                        preloadQq[i][j]              = 0;
+                        preloadSigma6[i][j]          = 0;
+                        preloadGapsysSigma6VdW[i][j] = 0;
                     }
                 }
             }
@@ -629,21 +680,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 gapsysScaleLinpointVdWEff;
+            RealType gmx_unused gapsysScaleLinpointCoulEff;
+            RealType gmx_unused gapsysSigma6VdWEff[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 (useSoftCore)
+                if constexpr (softcoreType == KernelSoftcoreType::Beutler)
                 {
                     sigma6[i] = gmx::load<RealType>(preloadSigma6[i]);
                 }
+                if constexpr (softcoreType == KernelSoftcoreType::Gapsys)
+                {
+                    gapsysSigma6VdWEff[i] = gmx::load<RealType>(preloadGapsysSigma6VdW[i]);
+                }
             }
-            if constexpr (useSoftCore)
+            if constexpr (softcoreType == KernelSoftcoreType::Beutler)
             {
                 alphaVdwEff  = gmx::load<RealType>(preloadAlphaVdwEff);
                 alphaCoulEff = gmx::load<RealType>(preloadAlphaCoulEff);
             }
+            if constexpr (softcoreType == KernelSoftcoreType::Gapsys)
+            {
+                gapsysScaleLinpointVdWEff  = gmx::load<RealType>(preloadGapsysScaleLinpointVdW);
+                gapsysScaleLinpointCoulEff = gmx::load<RealType>(preloadGapsysScaleLinpointCoul);
+            }
 
             // Avoid overflow of r^-12 at distances near zero
             rSq  = gmx::max(rSq, minDistanceSquared);
@@ -651,7 +714,7 @@ static void nb_free_energy_kernel(const t_nblist&
             r    = rSq * rInv;
 
             RealType gmx_unused rp, rpm2;
-            if constexpr (useSoftCore)
+            if constexpr (softcoreType == KernelSoftcoreType::Beutler)
             {
                 rpm2 = rSq * rSq;  /* r4 */
                 rp   = rpm2 * rSq; /* r6 */
@@ -689,7 +752,7 @@ static void nb_free_energy_kernel(const t_nblist&
                                              && bPairIncluded && withinCutoffMask);
                     if (gmx::anyTrue(nonZeroState))
                     {
-                        if constexpr (useSoftCore)
+                        if constexpr (softcoreType == KernelSoftcoreType::Beutler)
                         {
                             RealType divisor = (alphaCoulEff * lFacCoul[i] * sigma6[i] + rp);
                             rPInvC           = gmx::inv(divisor);
@@ -742,6 +805,22 @@ static void nb_free_energy_kernel(const t_nblist&
                                 {
                                     fScalC[i] = ewaldScalarForce(qq[i], rInvC);
                                 }
+
+                                if constexpr (softcoreType == KernelSoftcoreType::Gapsys)
+                                {
+                                    ewaldQuadraticPotential<computeForces>(qq[i],
+                                                                           facel,
+                                                                           rC,
+                                                                           rCutoffCoul,
+                                                                           LFC[i],
+                                                                           DLF[i],
+                                                                           gapsysScaleLinpointCoulEff,
+                                                                           sh_ewald,
+                                                                           &fScalC[i],
+                                                                           &vCoul[i],
+                                                                           &dvdlCoul,
+                                                                           computeElecInteraction);
+                                }
                             }
                             else
                             {
@@ -750,6 +829,24 @@ static void nb_free_energy_kernel(const t_nblist&
                                 {
                                     fScalC[i] = reactionFieldScalarForce(qq[i], rInvC, rC, krf, two);
                                 }
+
+                                if constexpr (softcoreType == KernelSoftcoreType::Gapsys)
+                                {
+                                    reactionFieldQuadraticPotential<computeForces>(
+                                            qq[i],
+                                            facel,
+                                            rC,
+                                            rCutoffCoul,
+                                            LFC[i],
+                                            DLF[i],
+                                            gapsysScaleLinpointCoulEff,
+                                            krf,
+                                            crf,
+                                            &fScalC[i],
+                                            &vCoul[i],
+                                            &dvdlCoul,
+                                            computeElecInteraction);
+                                }
                             }
 
                             vCoul[i] = gmx::selectByMask(vCoul[i], computeElecInteraction);
@@ -777,7 +874,7 @@ static void nb_free_energy_kernel(const t_nblist&
                         if (gmx::anyTrue(computeVdwInteraction))
                         {
                             RealType rInv6;
-                            if constexpr (useSoftCore)
+                            if constexpr (softcoreType == KernelSoftcoreType::Beutler)
                             {
                                 rInv6 = rPInvV;
                             }
@@ -801,6 +898,24 @@ static void nb_free_energy_kernel(const t_nblist&
                                 fScalV[i] = lennardJonesScalarForce(vVdw6, vVdw12);
                             }
 
+                            if constexpr (softcoreType == KernelSoftcoreType::Gapsys)
+                            {
+                                lennardJonesQuadraticPotential<computeForces>(c6[i],
+                                                                              c12[i],
+                                                                              r,
+                                                                              rSq,
+                                                                              LFV[i],
+                                                                              DLF[i],
+                                                                              gapsysSigma6VdWEff[i],
+                                                                              gapsysScaleLinpointVdWEff,
+                                                                              repulsionShift,
+                                                                              dispersionShift,
+                                                                              &fScalV[i],
+                                                                              &vVdw[i],
+                                                                              &dvdlVdw,
+                                                                              computeVdwInteraction);
+                            }
+
                             if constexpr (vdwInteractionTypeIsEwald)
                             {
                                 /* Subtract the grid potential at the cut-off */
@@ -864,7 +979,7 @@ static void nb_free_energy_kernel(const t_nblist&
                             fScal = fScal + LFV[i] * fScalV[i] * rpm2;
                         }
 
-                        if constexpr (useSoftCore)
+                        if constexpr (softcoreType == KernelSoftcoreType::Beutler)
                         {
                             dvdlCoul = dvdlCoul + vCoul[i] * DLF[i]
                                        + LFC[i] * alphaCoulEff * dlFacCoul[i] * fScalC[i] * sigma6[i];
@@ -1063,56 +1178,56 @@ typedef void (*KernelFunction)(const t_nblist&
                                gmx::ArrayRef<real> threadVv,
                                gmx::ArrayRef<real> threadDvdl);
 
-template<bool useSoftCore, bool scLambdasOrAlphasDiffer, bool vdwInteractionTypeIsEwald, bool elecInteractionTypeIsEwald, bool vdwModifierIsPotSwitch, bool computeForces>
+template<KernelSoftcoreType softcoreType, bool scLambdasOrAlphasDiffer, bool vdwInteractionTypeIsEwald, bool elecInteractionTypeIsEwald, bool vdwModifierIsPotSwitch, bool computeForces>
 static KernelFunction dispatchKernelOnUseSimd(const bool useSimd)
 {
     if (useSimd)
     {
 #if GMX_SIMD_HAVE_REAL && GMX_SIMD_HAVE_INT32_ARITHMETICS && GMX_USE_SIMD_KERNELS
-        return (nb_free_energy_kernel<SimdDataTypes, useSoftCore, scLambdasOrAlphasDiffer, vdwInteractionTypeIsEwald, elecInteractionTypeIsEwald, vdwModifierIsPotSwitch, computeForces>);
+        return (nb_free_energy_kernel<SimdDataTypes, softcoreType, scLambdasOrAlphasDiffer, vdwInteractionTypeIsEwald, elecInteractionTypeIsEwald, vdwModifierIsPotSwitch, computeForces>);
 #else
-        return (nb_free_energy_kernel<ScalarDataTypes, useSoftCore, scLambdasOrAlphasDiffer, vdwInteractionTypeIsEwald, elecInteractionTypeIsEwald, vdwModifierIsPotSwitch, computeForces>);
+        return (nb_free_energy_kernel<ScalarDataTypes, softcoreType, scLambdasOrAlphasDiffer, vdwInteractionTypeIsEwald, elecInteractionTypeIsEwald, vdwModifierIsPotSwitch, computeForces>);
 #endif
     }
     else
     {
-        return (nb_free_energy_kernel<ScalarDataTypes, useSoftCore, scLambdasOrAlphasDiffer, vdwInteractionTypeIsEwald, elecInteractionTypeIsEwald, vdwModifierIsPotSwitch, computeForces>);
+        return (nb_free_energy_kernel<ScalarDataTypes, softcoreType, scLambdasOrAlphasDiffer, vdwInteractionTypeIsEwald, elecInteractionTypeIsEwald, vdwModifierIsPotSwitch, computeForces>);
     }
 }
 
-template<bool useSoftCore, bool scLambdasOrAlphasDiffer, bool vdwInteractionTypeIsEwald, bool elecInteractionTypeIsEwald, bool vdwModifierIsPotSwitch>
+template<KernelSoftcoreType softcoreType, bool scLambdasOrAlphasDiffer, bool vdwInteractionTypeIsEwald, bool elecInteractionTypeIsEwald, bool vdwModifierIsPotSwitch>
 static KernelFunction dispatchKernelOnComputeForces(const bool computeForces, const bool useSimd)
 {
     if (computeForces)
     {
-        return (dispatchKernelOnUseSimd<useSoftCore, scLambdasOrAlphasDiffer, vdwInteractionTypeIsEwald, elecInteractionTypeIsEwald, vdwModifierIsPotSwitch, true>(
+        return (dispatchKernelOnUseSimd<softcoreType, scLambdasOrAlphasDiffer, vdwInteractionTypeIsEwald, elecInteractionTypeIsEwald, vdwModifierIsPotSwitch, true>(
                 useSimd));
     }
     else
     {
-        return (dispatchKernelOnUseSimd<useSoftCore, scLambdasOrAlphasDiffer, vdwInteractionTypeIsEwald, elecInteractionTypeIsEwald, vdwModifierIsPotSwitch, false>(
+        return (dispatchKernelOnUseSimd<softcoreType, scLambdasOrAlphasDiffer, vdwInteractionTypeIsEwald, elecInteractionTypeIsEwald, vdwModifierIsPotSwitch, false>(
                 useSimd));
     }
 }
 
-template<bool useSoftCore, bool scLambdasOrAlphasDiffer, bool vdwInteractionTypeIsEwald, bool elecInteractionTypeIsEwald>
+template<KernelSoftcoreType softcoreType, bool scLambdasOrAlphasDiffer, bool vdwInteractionTypeIsEwald, bool elecInteractionTypeIsEwald>
 static KernelFunction dispatchKernelOnVdwModifier(const bool vdwModifierIsPotSwitch,
                                                   const bool computeForces,
                                                   const bool useSimd)
 {
     if (vdwModifierIsPotSwitch)
     {
-        return (dispatchKernelOnComputeForces<useSoftCore, scLambdasOrAlphasDiffer, vdwInteractionTypeIsEwald, elecInteractionTypeIsEwald, true>(
+        return (dispatchKernelOnComputeForces<softcoreType, scLambdasOrAlphasDiffer, vdwInteractionTypeIsEwald, elecInteractionTypeIsEwald, true>(
                 computeForces, useSimd));
     }
     else
     {
-        return (dispatchKernelOnComputeForces<useSoftCore, scLambdasOrAlphasDiffer, vdwInteractionTypeIsEwald, elecInteractionTypeIsEwald, false>(
+        return (dispatchKernelOnComputeForces<softcoreType, scLambdasOrAlphasDiffer, vdwInteractionTypeIsEwald, elecInteractionTypeIsEwald, false>(
                 computeForces, useSimd));
     }
 }
 
-template<bool useSoftCore, bool scLambdasOrAlphasDiffer, bool vdwInteractionTypeIsEwald>
+template<KernelSoftcoreType softcoreType, bool scLambdasOrAlphasDiffer, bool vdwInteractionTypeIsEwald>
 static KernelFunction dispatchKernelOnElecInteractionType(const bool elecInteractionTypeIsEwald,
                                                           const bool vdwModifierIsPotSwitch,
                                                           const bool computeForces,
@@ -1120,17 +1235,17 @@ static KernelFunction dispatchKernelOnElecInteractionType(const bool elecInterac
 {
     if (elecInteractionTypeIsEwald)
     {
-        return (dispatchKernelOnVdwModifier<useSoftCore, scLambdasOrAlphasDiffer, vdwInteractionTypeIsEwald, true>(
+        return (dispatchKernelOnVdwModifier<softcoreType, scLambdasOrAlphasDiffer, vdwInteractionTypeIsEwald, true>(
                 vdwModifierIsPotSwitch, computeForces, useSimd));
     }
     else
     {
-        return (dispatchKernelOnVdwModifier<useSoftCore, scLambdasOrAlphasDiffer, vdwInteractionTypeIsEwald, false>(
+        return (dispatchKernelOnVdwModifier<softcoreType, scLambdasOrAlphasDiffer, vdwInteractionTypeIsEwald, false>(
                 vdwModifierIsPotSwitch, computeForces, useSimd));
     }
 }
 
-template<bool useSoftCore, bool scLambdasOrAlphasDiffer>
+template<KernelSoftcoreType softcoreType, bool scLambdasOrAlphasDiffer>
 static KernelFunction dispatchKernelOnVdwInteractionType(const bool vdwInteractionTypeIsEwald,
                                                          const bool elecInteractionTypeIsEwald,
                                                          const bool vdwModifierIsPotSwitch,
@@ -1139,17 +1254,17 @@ static KernelFunction dispatchKernelOnVdwInteractionType(const bool vdwInteracti
 {
     if (vdwInteractionTypeIsEwald)
     {
-        return (dispatchKernelOnElecInteractionType<useSoftCore, scLambdasOrAlphasDiffer, true>(
+        return (dispatchKernelOnElecInteractionType<softcoreType, scLambdasOrAlphasDiffer, true>(
                 elecInteractionTypeIsEwald, vdwModifierIsPotSwitch, computeForces, useSimd));
     }
     else
     {
-        return (dispatchKernelOnElecInteractionType<useSoftCore, scLambdasOrAlphasDiffer, false>(
+        return (dispatchKernelOnElecInteractionType<softcoreType, scLambdasOrAlphasDiffer, false>(
                 elecInteractionTypeIsEwald, vdwModifierIsPotSwitch, computeForces, useSimd));
     }
 }
 
-template<bool useSoftCore>
+template<KernelSoftcoreType softcoreType>
 static KernelFunction dispatchKernelOnScLambdasOrAlphasDifference(const bool scLambdasOrAlphasDiffer,
                                                                   const bool vdwInteractionTypeIsEwald,
                                                                   const bool elecInteractionTypeIsEwald,
@@ -1159,12 +1274,12 @@ static KernelFunction dispatchKernelOnScLambdasOrAlphasDifference(const bool scL
 {
     if (scLambdasOrAlphasDiffer)
     {
-        return (dispatchKernelOnVdwInteractionType<useSoftCore, true>(
+        return (dispatchKernelOnVdwInteractionType<softcoreType, true>(
                 vdwInteractionTypeIsEwald, elecInteractionTypeIsEwald, vdwModifierIsPotSwitch, computeForces, useSimd));
     }
     else
     {
-        return (dispatchKernelOnVdwInteractionType<useSoftCore, false>(
+        return (dispatchKernelOnVdwInteractionType<softcoreType, false>(
                 vdwInteractionTypeIsEwald, elecInteractionTypeIsEwald, vdwModifierIsPotSwitch, computeForces, useSimd));
     }
 }
@@ -1177,23 +1292,46 @@ static KernelFunction dispatchKernel(const bool                 scLambdasOrAlpha
                                      const bool                 useSimd,
                                      const interaction_const_t& ic)
 {
-    if (ic.softCoreParameters->alphaCoulomb == 0 && ic.softCoreParameters->alphaVdw == 0)
+    const auto& scParams = *ic.softCoreParameters;
+    if (scParams.softcoreType == SoftcoreType::Beutler)
     {
-        return (dispatchKernelOnScLambdasOrAlphasDifference<false>(scLambdasOrAlphasDiffer,
-                                                                   vdwInteractionTypeIsEwald,
-                                                                   elecInteractionTypeIsEwald,
-                                                                   vdwModifierIsPotSwitch,
-                                                                   computeForces,
-                                                                   useSimd));
+        if (scParams.alphaCoulomb == 0 && scParams.alphaVdw == 0)
+        {
+            return (dispatchKernelOnScLambdasOrAlphasDifference<KernelSoftcoreType::None>(
+                    scLambdasOrAlphasDiffer,
+                    vdwInteractionTypeIsEwald,
+                    elecInteractionTypeIsEwald,
+                    vdwModifierIsPotSwitch,
+                    computeForces,
+                    useSimd));
+        }
+        return (dispatchKernelOnScLambdasOrAlphasDifference<KernelSoftcoreType::Beutler>(
+                scLambdasOrAlphasDiffer,
+                vdwInteractionTypeIsEwald,
+                elecInteractionTypeIsEwald,
+                vdwModifierIsPotSwitch,
+                computeForces,
+                useSimd));
     }
-    else
+    else // Gapsys
     {
-        return (dispatchKernelOnScLambdasOrAlphasDifference<true>(scLambdasOrAlphasDiffer,
-                                                                  vdwInteractionTypeIsEwald,
-                                                                  elecInteractionTypeIsEwald,
-                                                                  vdwModifierIsPotSwitch,
-                                                                  computeForces,
-                                                                  useSimd));
+        if (scParams.gapsysScaleLinpointCoul == 0 && scParams.gapsysScaleLinpointVdW == 0)
+        {
+            return (dispatchKernelOnScLambdasOrAlphasDifference<KernelSoftcoreType::None>(
+                    scLambdasOrAlphasDiffer,
+                    vdwInteractionTypeIsEwald,
+                    elecInteractionTypeIsEwald,
+                    vdwModifierIsPotSwitch,
+                    computeForces,
+                    useSimd));
+        }
+        return (dispatchKernelOnScLambdasOrAlphasDifference<KernelSoftcoreType::Gapsys>(
+                scLambdasOrAlphasDiffer,
+                vdwInteractionTypeIsEwald,
+                elecInteractionTypeIsEwald,
+                vdwModifierIsPotSwitch,
+                computeForces,
+                useSimd));
     }
 }