From b72a2df243c36d5022783844d027e1ae43a4525b Mon Sep 17 00:00:00 2001 From: Magnus Lundborg Date: Wed, 19 Jun 2019 16:13:33 +0200 Subject: [PATCH] Inline functions of calculations of FE interaction types. The functions are so far only templated on the real requirements of the softcore that is used. This is one step towards templating the calculations for SIMD. Refs #2997. Change-Id: I3fd119dce30f95eba4f8cd6f139f99260acf0e22 --- .../gmxlib/nonbonded/nb_free_energy.cpp | 180 +++++++++++++----- 1 file changed, 137 insertions(+), 43 deletions(-) diff --git a/src/gromacs/gmxlib/nonbonded/nb_free_energy.cpp b/src/gromacs/gmxlib/nonbonded/nb_free_energy.cpp index 8303ae6007..cd821ef414 100644 --- a/src/gromacs/gmxlib/nonbonded/nb_free_energy.cpp +++ b/src/gromacs/gmxlib/nonbonded/nb_free_energy.cpp @@ -110,6 +110,119 @@ inline void pthRoot(const double r, *invPthRoot = 1/(*pthRoot); } +template +static inline real calculateSigmaPow(const real sigma6) +{ + if (softCoreTreatment == SoftCoreTreatment::RPower6) + { + return sigma6; + } + else + { + real sigmaPow = sigma6 * sigma6; /* sigma^12 */ + sigmaPow = sigmaPow * sigmaPow; /* sigma^24 */ + sigmaPow = sigmaPow * sigmaPow; /* sigma^48 */ + return (sigmaPow); + } +} + +template +static inline real calculateRinv6(const SCReal rinvV) +{ + if (softCoreTreatment == SoftCoreTreatment::RPower6) + { + return rinvV; + } + else + { + real rinv6 = rinvV * rinvV; + return (rinv6 * rinv6 * rinv6); + } +} + +static inline real calculateVdw6(const real c6, const real rinv6) +{ + return (c6 * rinv6); +} + +static inline real calculateVdw12(const real c12, const real rinv6) +{ + return (c12 * rinv6 * rinv6); +} + +/* reaction-field electrostatics */ +template +static inline SCReal +reactionFieldScalarForce(const real qq, const real rinv, const SCReal r, const real krf, + const real two) +{ + return (qq*(rinv - two*krf*r*r)); +} +template +static inline real +reactionFieldPotential(const real qq, const real rinv, const SCReal r, const real krf, + const real potentialShift) +{ + return (qq*(rinv + krf*r*r-potentialShift)); +} + +/* Ewald electrostatics */ +static inline real +ewaldScalarForce(const real coulomb, const real rinv) +{ + return (coulomb*rinv); +} +static inline real +ewaldPotential(const real coulomb, const real rinv, const real potentialShift) +{ + return (coulomb*(rinv-potentialShift)); +} + +/* cutoff LJ */ +static inline real +lennardJonesScalarForce(const real v6, const real v12) +{ + return (v12 - v6); +} +static inline real +lennardJonesPotential(const real v6, const real v12, const real c6, const real c12, const real repulsionShift, + const real dispersionShift, const real onesixth, const real onetwelfth) +{ + return ((v12 + c12*repulsionShift)*onetwelfth - (v6 + c6*dispersionShift)*onesixth); +} + +/* Ewald LJ */ +static inline real +ewaldLennardJonesGridSubtract(const real c6grid, const real potentialShift, const real onesixth) +{ + return (c6grid * potentialShift * onesixth); +} + +/* LJ Potential switch */ +template +static inline SCReal +potSwitchScalarForceMod(const SCReal fScalarInp, const real potential, const real sw, const SCReal r, const real rVdw, const real dsw, const real zero) +{ + if (r < rVdw) + { + SCReal fScalar = fScalarInp * sw - r * potential * dsw; + return (fScalar); + } + return (zero); +} +template +static inline real +potSwitchPotentialMod(const real potentialInp, const real sw, const SCReal r, const real rVdw, const real zero) +{ + if (r < rVdw) + { + real potential = potentialInp * sw; + return (potential); + } + return (zero); +} + + //! Templated free-energy non-bonded kernel template static void @@ -133,13 +246,12 @@ nb_free_energy_kernel(const t_nblist * gmx_restrict nlist, real tx, ty, tz, Fscal; SCReal FscalC[NSTATES], FscalV[NSTATES]; /* Needs double for sc_power==48 */ real Vcoul[NSTATES], Vvdw[NSTATES]; - real rinv6, r; real iqA, iqB; real qq[NSTATES], vctot; int ntiA, ntiB, tj[NSTATES]; - real Vvdw6, Vvdw12, vvtot; + real vvtot; real ix, iy, iz, fix, fiy, fiz; - real dx, dy, dz, rsq, rinv; + real dx, dy, dz, r, rsq, rinv; real c6[NSTATES], c12[NSTATES], c6grid; real LFC[NSTATES], LFV[NSTATES], DLF[NSTATES]; SCReal dvdl_coul, dvdl_vdw; @@ -173,7 +285,6 @@ nb_free_energy_kernel(const t_nblist * gmx_restrict nlist, real rcutoff_max2; const real * tab_ewald_F_lj = nullptr; const real * tab_ewald_V_lj = nullptr; - real d, d2, sw, dsw; real vdw_swV3, vdw_swV4, vdw_swV5, vdw_swF2, vdw_swF3, vdw_swF4; gmx_bool bComputeVdwInteraction, bComputeElecInteraction; const real * ewtab = nullptr; @@ -243,7 +354,7 @@ nb_free_energy_kernel(const t_nblist * gmx_restrict nlist, if (ic->vdw_modifier == eintmodPOTSWITCH) { - d = ic->rvdw - ic->rvdw_switch; + real d = ic->rvdw - ic->rvdw_switch; vdw_swV3 = -10.0/(d*d*d); vdw_swV4 = 15.0/(d*d*d*d); vdw_swV5 = -6.0/(d*d*d*d*d); @@ -458,16 +569,7 @@ nb_free_energy_kernel(const t_nblist * gmx_restrict nlist, { sigma6[i] = sigma6_def; } - if (softCoreTreatment == SoftCoreTreatment::RPower6) - { - sigma_pow[i] = sigma6[i]; - } - else - { - sigma_pow[i] = sigma6[i]*sigma6[i]; /* sigma^12 */ - sigma_pow[i] = sigma_pow[i]*sigma_pow[i]; /* sigma^24 */ - sigma_pow[i] = sigma_pow[i]*sigma_pow[i]; /* sigma^48 */ - } + sigma_pow[i] = calculateSigmaPow(sigma6[i]); } } @@ -538,15 +640,13 @@ nb_free_energy_kernel(const t_nblist * gmx_restrict nlist, { if (bEwald) { - /* Ewald FEP is done only on the 1/r part */ - Vcoul[i] = qq[i]*(rinvC-sh_ewald); - FscalC[i] = qq[i]*rinvC; + Vcoul[i] = ewaldPotential(qq[i], rinvC, sh_ewald); + FscalC[i] = ewaldScalarForce(qq[i], rinvC); } else { - /* reaction-field */ - Vcoul[i] = qq[i]*(rinvC + krf*rC*rC-crf); - FscalC[i] = qq[i]*(rinvC - two*krf*rC*rC); + Vcoul[i] = reactionFieldPotential(qq[i], rinvC, rC, krf, crf); + FscalC[i] = reactionFieldScalarForce(qq[i], rinvC, rC, krf, two); } } @@ -560,43 +660,37 @@ nb_free_energy_kernel(const t_nblist * gmx_restrict nlist, (!bEwaldLJ && rV < rvdw); if ((c6[i] != 0 || c12[i] != 0) && bComputeVdwInteraction) { - /* cutoff LJ, also handles part of Ewald LJ */ + real rinv6; if (softCoreTreatment == SoftCoreTreatment::RPower6) { - rinv6 = rpinvV; + rinv6 = calculateRinv6(rpinvV); } else { - rinv6 = rinvV*rinvV; - rinv6 = rinv6*rinv6*rinv6; + rinv6 = calculateRinv6(rinvV); } - Vvdw6 = c6[i]*rinv6; - Vvdw12 = c12[i]*rinv6*rinv6; + real Vvdw6 = calculateVdw6(c6[i], rinv6); + real Vvdw12 = calculateVdw12(c12[i], rinv6); - Vvdw[i] = ( (Vvdw12 + c12[i]*repulsionShift)*onetwelfth - - (Vvdw6 + c6[i]*dispersionShift)*onesixth); - FscalV[i] = Vvdw12 - Vvdw6; + Vvdw[i] = lennardJonesPotential(Vvdw6, Vvdw12, c6[i], c12[i], repulsionShift, dispersionShift, onesixth, onetwelfth); + FscalV[i] = lennardJonesScalarForce(Vvdw6, Vvdw12); if (bEwaldLJ) { /* Subtract the grid potential at the cut-off */ - c6grid = nbfp_grid[tj[i]]; - Vvdw[i] += c6grid*sh_lj_ewald*onesixth; + Vvdw[i] += ewaldLennardJonesGridSubtract(nbfp_grid[tj[i]], sh_lj_ewald, onesixth); } if (ic->vdw_modifier == eintmodPOTSWITCH) { - d = rV - ic->rvdw_switch; - d = (d > zero) ? d : zero; - d2 = d*d; - sw = one+d2*d*(vdw_swV3+d*(vdw_swV4+d*vdw_swV5)); - dsw = d2*(vdw_swF2+d*(vdw_swF3+d*vdw_swF4)); - - FscalV[i] = FscalV[i]*sw - rV*Vvdw[i]*dsw; - Vvdw[i] *= sw; - - FscalV[i] = (rV < rvdw) ? FscalV[i] : zero; - Vvdw[i] = (rV < rvdw) ? Vvdw[i] : zero; + real d = rV - ic->rvdw_switch; + d = (d > zero) ? d : zero; + real d2 = d*d; + real sw = one+d2*d*(vdw_swV3+d*(vdw_swV4+d*vdw_swV5)); + real dsw = d2*(vdw_swF2+d*(vdw_swF3+d*vdw_swF4)); + + FscalV[i] = potSwitchScalarForceMod(FscalV[i], Vvdw[i], sw, rV, rvdw, dsw, zero); + Vvdw[i] = potSwitchPotentialMod(Vvdw[i], sw, rV, rvdw, zero); } } -- 2.22.0