src/gromacs/gmxlib/nonbonded/nb_free_energy.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include "nb_free_energy.h"
  41
  42 #include "config.h"
  43
  44 #include <cmath>
  45
  46 #include <algorithm>
  47
  48 #include "gromacs/gmxlib/nrnb.h"
  49 #include "gromacs/gmxlib/nonbonded/nb_kernel.h"
  50 #include "gromacs/gmxlib/nonbonded/nonbonded.h"
  51 #include "gromacs/math/functions.h"
  52 #include "gromacs/math/vec.h"
  53 #include "gromacs/mdtypes/forceoutput.h"
  54 #include "gromacs/mdtypes/forcerec.h"
  55 #include "gromacs/mdtypes/interaction_const.h"
  56 #include "gromacs/mdtypes/md_enums.h"
  57 #include "gromacs/mdtypes/mdatom.h"
  58 #include "gromacs/simd/simd.h"
  59 #include "gromacs/utility/fatalerror.h"
  60
  61
  62 //! Scalar (non-SIMD) data types.
  63 struct ScalarDataTypes
  64 {
  65     using RealType                     = real; //!< The data type to use as real.
  66     using IntType                      = int;  //!< The data type to use as int.
  67     static constexpr int simdRealWidth = 1;    //!< The width of the RealType.
  68     static constexpr int simdIntWidth  = 1;    //!< The width of the IntType.
  69 };
  70
  71 #if GMX_SIMD_HAVE_REAL && GMX_SIMD_HAVE_INT32_ARITHMETICS
  72 //! SIMD data types.
  73 struct SimdDataTypes
  74 {
  75     using RealType                     = gmx::SimdReal;         //!< The data type to use as real.
  76     using IntType                      = gmx::SimdInt32;        //!< The data type to use as int.
  77     static constexpr int simdRealWidth = GMX_SIMD_REAL_WIDTH;   //!< The width of the RealType.
  78     static constexpr int simdIntWidth  = GMX_SIMD_FINT32_WIDTH; //!< The width of the IntType.
  79 };
  80 #endif
  81
  82 //! Computes r^(1/p) and 1/r^(1/p) for the standard p=6
  83 template<class RealType>
  84 static inline void pthRoot(const RealType r, RealType* pthRoot, RealType* invPthRoot)
  85 {
  86     *invPthRoot = gmx::invsqrt(std::cbrt(r));
  87     *pthRoot    = 1 / (*invPthRoot);
  88 }
  89
  90 template<class RealType>
  91 static inline RealType calculateRinv6(const RealType rInvV)
  92 {
  93     RealType rInv6 = rInvV * rInvV;
  94     return (rInv6 * rInv6 * rInv6);
  95 }
  96
  97 template<class RealType>
  98 static inline RealType calculateVdw6(const RealType c6, const RealType rInv6)
  99 {
 100     return (c6 * rInv6);
 101 }
 102
 103 template<class RealType>
 104 static inline RealType calculateVdw12(const RealType c12, const RealType rInv6)
 105 {
 106     return (c12 * rInv6 * rInv6);
 107 }
 108
 109 /* reaction-field electrostatics */
 110 template<class RealType>
 111 static inline RealType reactionFieldScalarForce(const RealType qq,
 112                                                 const RealType rInv,
 113                                                 const RealType r,
 114                                                 const real     krf,
 115                                                 const real     two)
 116 {
 117     return (qq * (rInv - two * krf * r * r));
 118 }
 119 template<class RealType>
 120 static inline RealType reactionFieldPotential(const RealType qq,
 121                                               const RealType rInv,
 122                                               const RealType r,
 123                                               const real     krf,
 124                                               const real     potentialShift)
 125 {
 126     return (qq * (rInv + krf * r * r - potentialShift));
 127 }
 128
 129 /* Ewald electrostatics */
 130 template<class RealType>
 131 static inline RealType ewaldScalarForce(const RealType coulomb, const RealType rInv)
 132 {
 133     return (coulomb * rInv);
 134 }
 135 template<class RealType>
 136 static inline RealType ewaldPotential(const RealType coulomb, const RealType rInv, const real potentialShift)
 137 {
 138     return (coulomb * (rInv - potentialShift));
 139 }
 140
 141 /* cutoff LJ */
 142 template<class RealType>
 143 static inline RealType lennardJonesScalarForce(const RealType v6, const RealType v12)
 144 {
 145     return (v12 - v6);
 146 }
 147 template<class RealType>
 148 static inline RealType lennardJonesPotential(const RealType v6,
 149                                              const RealType v12,
 150                                              const RealType c6,
 151                                              const RealType c12,
 152                                              const real     repulsionShift,
 153                                              const real     dispersionShift,
 154                                              const real     oneSixth,
 155                                              const real     oneTwelfth)
 156 {
 157     return ((v12 + c12 * repulsionShift) * oneTwelfth - (v6 + c6 * dispersionShift) * oneSixth);
 158 }
 159
 160 /* Ewald LJ */
 161 static inline real ewaldLennardJonesGridSubtract(const real c6grid, const real potentialShift, const real oneSixth)
 162 {
 163     return (c6grid * potentialShift * oneSixth);
 164 }
 165
 166 /* LJ Potential switch */
 167 template<class RealType>
 168 static inline RealType potSwitchScalarForceMod(const RealType fScalarInp,
 169                                                const RealType potential,
 170                                                const RealType sw,
 171                                                const RealType r,
 172                                                const RealType rVdw,
 173                                                const RealType dsw,
 174                                                const real     zero)
 175 {
 176     if (r < rVdw)
 177     {
 178         real fScalar = fScalarInp * sw - r * potential * dsw;
 179         return (fScalar);
 180     }
 181     return (zero);
 182 }
 183 template<class RealType>
 184 static inline RealType potSwitchPotentialMod(const RealType potentialInp,
 185                                              const RealType sw,
 186                                              const RealType r,
 187                                              const RealType rVdw,
 188                                              const real     zero)
 189 {
 190     if (r < rVdw)
 191     {
 192         real potential = potentialInp * sw;
 193         return (potential);
 194     }
 195     return (zero);
 196 }
 197
 198
 199 //! Templated free-energy non-bonded kernel
 200 template<typename DataTypes, bool useSoftCore, bool scLambdasOrAlphasDiffer, bool vdwInteractionTypeIsEwald, bool elecInteractionTypeIsEwald, bool vdwModifierIsPotSwitch>
 201 static void nb_free_energy_kernel(const t_nblist* gmx_restrict nlist,
 202                                   rvec* gmx_restrict         xx,
 203                                   gmx::ForceWithShiftForces* forceWithShiftForces,
 204                                   const t_forcerec* gmx_restrict fr,
 205                                   const t_mdatoms* gmx_restrict mdatoms,
 206                                   nb_kernel_data_t* gmx_restrict kernel_data,
 207                                   t_nrnb* gmx_restrict nrnb)
 208 {
 209 #define STATE_A 0
 210 #define STATE_B 1
 211 #define NSTATES 2
 212
 213     using RealType = typename DataTypes::RealType;
 214     using IntType  = typename DataTypes::IntType;
 215
 216     /* FIXME: How should these be handled with SIMD? */
 217     constexpr real oneTwelfth = 1.0 / 12.0;
 218     constexpr real oneSixth   = 1.0 / 6.0;
 219     constexpr real zero       = 0.0;
 220     constexpr real half       = 0.5;
 221     constexpr real one        = 1.0;
 222     constexpr real two        = 2.0;
 223     constexpr real six        = 6.0;
 224
 225     /* Extract pointer to non-bonded interaction constants */
 226     const interaction_const_t* ic = fr->ic;
 227
 228     // Extract pair list data
 229     const int  nri    = nlist->nri;
 230     const int* iinr   = nlist->iinr;
 231     const int* jindex = nlist->jindex;
 232     const int* jjnr   = nlist->jjnr;
 233     const int* shift  = nlist->shift;
 234     const int* gid    = nlist->gid;
 235
 236     const real* shiftvec      = fr->shift_vec[0];
 237     const real* chargeA       = mdatoms->chargeA;
 238     const real* chargeB       = mdatoms->chargeB;
 239     real*       Vc            = kernel_data->energygrp_elec;
 240     const int*  typeA         = mdatoms->typeA;
 241     const int*  typeB         = mdatoms->typeB;
 242     const int   ntype         = fr->ntype;
 243     const real* nbfp          = fr->nbfp.data();
 244     const real* nbfp_grid     = fr->ljpme_c6grid;
 245     real*       Vv            = kernel_data->energygrp_vdw;
 246     const real  lambda_coul   = kernel_data->lambda[efptCOUL];
 247     const real  lambda_vdw    = kernel_data->lambda[efptVDW];
 248     real*       dvdl          = kernel_data->dvdl;
 249     const auto& scParams      = *ic->softCoreParameters;
 250     const real  alpha_coul    = scParams.alphaCoulomb;
 251     const real  alpha_vdw     = scParams.alphaVdw;
 252     const real  lam_power     = scParams.lambdaPower;
 253     const real  sigma6_def    = scParams.sigma6WithInvalidSigma;
 254     const real  sigma6_min    = scParams.sigma6Minimum;
 255     const bool  doForces      = ((kernel_data->flags & GMX_NONBONDED_DO_FORCE) != 0);
 256     const bool  doShiftForces = ((kernel_data->flags & GMX_NONBONDED_DO_SHIFTFORCE) != 0);
 257     const bool  doPotential   = ((kernel_data->flags & GMX_NONBONDED_DO_POTENTIAL) != 0);
 258
 259     // Extract data from interaction_const_t
 260     const real facel           = ic->epsfac;
 261     const real rCoulomb        = ic->rcoulomb;
 262     const real krf             = ic->k_rf;
 263     const real crf             = ic->c_rf;
 264     const real shLjEwald       = ic->sh_lj_ewald;
 265     const real rVdw            = ic->rvdw;
 266     const real dispersionShift = ic->dispersion_shift.cpot;
 267     const real repulsionShift  = ic->repulsion_shift.cpot;
 268
 269     // Note that the nbnxm kernels do not support Coulomb potential switching at all
 270     GMX_ASSERT(ic->coulomb_modifier != eintmodPOTSWITCH,
 271                "Potential switching is not supported for Coulomb with FEP");
 272
 273     real vdw_swV3, vdw_swV4, vdw_swV5, vdw_swF2, vdw_swF3, vdw_swF4;
 274     if (vdwModifierIsPotSwitch)
 275     {
 276         const real d = ic->rvdw - ic->rvdw_switch;
 277         vdw_swV3     = -10.0 / (d * d * d);
 278         vdw_swV4     = 15.0 / (d * d * d * d);
 279         vdw_swV5     = -6.0 / (d * d * d * d * d);
 280         vdw_swF2     = -30.0 / (d * d * d);
 281         vdw_swF3     = 60.0 / (d * d * d * d);
 282         vdw_swF4     = -30.0 / (d * d * d * d * d);
 283     }
 284     else
 285     {
 286         /* Avoid warnings from stupid compilers (looking at you, Clang!) */
 287         vdw_swV3 = vdw_swV4 = vdw_swV5 = vdw_swF2 = vdw_swF3 = vdw_swF4 = 0.0;
 288     }
 289
 290     int icoul;
 291     if (ic->eeltype == eelCUT || EEL_RF(ic->eeltype))
 292     {
 293         icoul = GMX_NBKERNEL_ELEC_REACTIONFIELD;
 294     }
 295     else
 296     {
 297         icoul = GMX_NBKERNEL_ELEC_NONE;
 298     }
 299
 300     real rcutoff_max2 = std::max(ic->rcoulomb, ic->rvdw);
 301     rcutoff_max2      = rcutoff_max2 * rcutoff_max2;
 302
 303     const real* tab_ewald_F_lj           = nullptr;
 304     const real* tab_ewald_V_lj           = nullptr;
 305     const real* ewtab                    = nullptr;
 306     real        coulombTableScale        = 0;
 307     real        coulombTableScaleInvHalf = 0;
 308     real        vdwTableScale            = 0;
 309     real        vdwTableScaleInvHalf     = 0;
 310     real        sh_ewald                 = 0;
 311     if (elecInteractionTypeIsEwald || vdwInteractionTypeIsEwald)
 312     {
 313         sh_ewald = ic->sh_ewald;
 314     }
 315     if (elecInteractionTypeIsEwald)
 316     {
 317         const auto& coulombTables = *ic->coulombEwaldTables;
 318         ewtab                     = coulombTables.tableFDV0.data();
 319         coulombTableScale         = coulombTables.scale;
 320         coulombTableScaleInvHalf  = half / coulombTableScale;
 321     }
 322     if (vdwInteractionTypeIsEwald)
 323     {
 324         const auto& vdwTables = *ic->vdwEwaldTables;
 325         tab_ewald_F_lj        = vdwTables.tableF.data();
 326         tab_ewald_V_lj        = vdwTables.tableV.data();
 327         vdwTableScale         = vdwTables.scale;
 328         vdwTableScaleInvHalf  = half / vdwTableScale;
 329     }
 330
 331     /* For Ewald/PME interactions we cannot easily apply the soft-core component to
 332      * reciprocal space. When we use non-switched Ewald interactions, we
 333      * assume the soft-coring does not significantly affect the grid contribution
 334      * and apply the soft-core only to the full 1/r (- shift) pair contribution.
 335      *
 336      * However, we cannot use this approach for switch-modified since we would then
 337      * effectively end up evaluating a significantly different interaction here compared to the
 338      * normal (non-free-energy) kernels, either by applying a cutoff at a different
 339      * position than what the user requested, or by switching different
 340      * things (1/r rather than short-range Ewald). For these settings, we just
 341      * use the traditional short-range Ewald interaction in that case.
 342      */
 343     GMX_RELEASE_ASSERT(!(vdwInteractionTypeIsEwald && vdwModifierIsPotSwitch),
 344                        "Can not apply soft-core to switched Ewald potentials");
 345
 346     real dvdlCoul = 0;
 347     real dvdlVdw  = 0;
 348
 349     /* Lambda factor for state A, 1-lambda*/
 350     real LFC[NSTATES], LFV[NSTATES];
 351     LFC[STATE_A] = one - lambda_coul;
 352     LFV[STATE_A] = one - lambda_vdw;
 353
 354     /* Lambda factor for state B, lambda*/
 355     LFC[STATE_B] = lambda_coul;
 356     LFV[STATE_B] = lambda_vdw;
 357
 358     /*derivative of the lambda factor for state A and B */
 359     real DLF[NSTATES];
 360     DLF[STATE_A] = -1;
 361     DLF[STATE_B] = 1;
 362
 363     real           lFacCoul[NSTATES], dlFacCoul[NSTATES], lFacVdw[NSTATES], dlFacVdw[NSTATES];
 364     constexpr real sc_r_power = 6.0_real;
 365     for (int i = 0; i < NSTATES; i++)
 366     {
 367         lFacCoul[i]  = (lam_power == 2 ? (1 - LFC[i]) * (1 - LFC[i]) : (1 - LFC[i]));
 368         dlFacCoul[i] = DLF[i] * lam_power / sc_r_power * (lam_power == 2 ? (1 - LFC[i]) : 1);
 369         lFacVdw[i]   = (lam_power == 2 ? (1 - LFV[i]) * (1 - LFV[i]) : (1 - LFV[i]));
 370         dlFacVdw[i]  = DLF[i] * lam_power / sc_r_power * (lam_power == 2 ? (1 - LFV[i]) : 1);
 371     }
 372
 373     // TODO: We should get rid of using pointers to real
 374     const real* x             = xx[0];
 375     real* gmx_restrict f      = &(forceWithShiftForces->force()[0][0]);
 376     real* gmx_restrict fshift = &(forceWithShiftForces->shiftForces()[0][0]);
 377
 378     const real rlistSquared = gmx::square(fr->rlist);
 379
 380     int numExcludedPairsBeyondRlist = 0;
 381
 382     for (int n = 0; n < nri; n++)
 383     {
 384         int npair_within_cutoff = 0;
 385
 386         const int  is3   = 3 * shift[n];
 387         const real shX   = shiftvec[is3];
 388         const real shY   = shiftvec[is3 + 1];
 389         const real shZ   = shiftvec[is3 + 2];
 390         const int  nj0   = jindex[n];
 391         const int  nj1   = jindex[n + 1];
 392         const int  ii    = iinr[n];
 393         const int  ii3   = 3 * ii;
 394         const real ix    = shX + x[ii3 + 0];
 395         const real iy    = shY + x[ii3 + 1];
 396         const real iz    = shZ + x[ii3 + 2];
 397         const real iqA   = facel * chargeA[ii];
 398         const real iqB   = facel * chargeB[ii];
 399         const int  ntiA  = 2 * ntype * typeA[ii];
 400         const int  ntiB  = 2 * ntype * typeB[ii];
 401         real       vCTot = 0;
 402         real       vVTot = 0;
 403         real       fIX   = 0;
 404         real       fIY   = 0;
 405         real       fIZ   = 0;
 406
 407         for (int k = nj0; k < nj1; k++)
 408         {
 409             int            tj[NSTATES];
 410             const int      jnr = jjnr[k];
 411             const int      j3  = 3 * jnr;
 412             RealType       c6[NSTATES], c12[NSTATES], qq[NSTATES], vCoul[NSTATES], vVdw[NSTATES];
 413             RealType       r, rInv, rp, rpm2;
 414             RealType       alphaVdwEff, alphaCoulEff, sigma6[NSTATES];
 415             const RealType dX  = ix - x[j3];
 416             const RealType dY  = iy - x[j3 + 1];
 417             const RealType dZ  = iz - x[j3 + 2];
 418             const RealType rSq = dX * dX + dY * dY + dZ * dZ;
 419             RealType       fScalC[NSTATES], fScalV[NSTATES];
 420             /* Check if this pair on the exlusions list.*/
 421             const bool bPairIncluded = nlist->excl_fep == nullptr || nlist->excl_fep[k];
 422
 423             if (rSq >= rcutoff_max2 && bPairIncluded)
 424             {
 425                 /* We save significant time by skipping all code below.
 426                  * Note that with soft-core interactions, the actual cut-off
 427                  * check might be different. But since the soft-core distance
 428                  * is always larger than r, checking on r here is safe.
 429                  * Exclusions outside the cutoff can not be skipped as
 430                  * when using Ewald: the reciprocal-space
 431                  * Ewald component still needs to be subtracted.
 432                  */
 433
 434                 continue;
 435             }
 436             npair_within_cutoff++;
 437
 438             if (rSq > rlistSquared)
 439             {
 440                 numExcludedPairsBeyondRlist++;
 441             }
 442
 443             if (rSq > 0)
 444             {
 445                 /* Note that unlike in the nbnxn kernels, we do not need
 446                  * to clamp the value of rSq before taking the invsqrt
 447                  * to avoid NaN in the LJ calculation, since here we do
 448                  * not calculate LJ interactions when C6 and C12 are zero.
 449                  */
 450
 451                 rInv = gmx::invsqrt(rSq);
 452                 r    = rSq * rInv;
 453             }
 454             else
 455             {
 456                 /* The force at r=0 is zero, because of symmetry.
 457                  * But note that the potential is in general non-zero,
 458                  * since the soft-cored r will be non-zero.
 459                  */
 460                 rInv = 0;
 461                 r    = 0;
 462             }
 463
 464             if (useSoftCore)
 465             {
 466                 rpm2 = rSq * rSq;  /* r4 */
 467                 rp   = rpm2 * rSq; /* r6 */
 468             }
 469             else
 470             {
 471                 /* The soft-core power p will not affect the results
 472                  * with not using soft-core, so we use power of 0 which gives
 473                  * the simplest math and cheapest code.
 474                  */
 475                 rpm2 = rInv * rInv;
 476                 rp   = 1;
 477             }
 478
 479             RealType fScal = 0;
 480
 481             qq[STATE_A] = iqA * chargeA[jnr];
 482             qq[STATE_B] = iqB * chargeB[jnr];
 483
 484             tj[STATE_A] = ntiA + 2 * typeA[jnr];
 485             tj[STATE_B] = ntiB + 2 * typeB[jnr];
 486
 487             if (bPairIncluded)
 488             {
 489                 c6[STATE_A] = nbfp[tj[STATE_A]];
 490                 c6[STATE_B] = nbfp[tj[STATE_B]];
 491
 492                 for (int i = 0; i < NSTATES; i++)
 493                 {
 494                     c12[i] = nbfp[tj[i] + 1];
 495                     if (useSoftCore)
 496                     {
 497                         if ((c6[i] > 0) && (c12[i] > 0))
 498                         {
 499                             /* c12 is stored scaled with 12.0 and c6 is scaled with 6.0 - correct for this */
 500                             sigma6[i] = half * c12[i] / c6[i];
 501                             if (sigma6[i] < sigma6_min) /* for disappearing coul and vdw with soft core at the same time */
 502                             {
 503                                 sigma6[i] = sigma6_min;
 504                             }
 505                         }
 506                         else
 507                         {
 508                             sigma6[i] = sigma6_def;
 509                         }
 510                     }
 511                 }
 512
 513                 if (useSoftCore)
 514                 {
 515                     /* only use softcore if one of the states has a zero endstate - softcore is for avoiding infinities!*/
 516                     if ((c12[STATE_A] > 0) && (c12[STATE_B] > 0))
 517                     {
 518                         alphaVdwEff  = 0;
 519                         alphaCoulEff = 0;
 520                     }
 521                     else
 522                     {
 523                         alphaVdwEff  = alpha_vdw;
 524                         alphaCoulEff = alpha_coul;
 525                     }
 526                 }
 527
 528                 for (int i = 0; i < NSTATES; i++)
 529                 {
 530                     fScalC[i] = 0;
 531                     fScalV[i] = 0;
 532                     vCoul[i]  = 0;
 533                     vVdw[i]   = 0;
 534
 535                     RealType rInvC, rInvV, rC, rV, rPInvC, rPInvV;
 536
 537                     /* Only spend time on A or B state if it is non-zero */
 538                     if ((qq[i] != 0) || (c6[i] != 0) || (c12[i] != 0))
 539                     {
 540                         /* this section has to be inside the loop because of the dependence on sigma6 */
 541                         if (useSoftCore)
 542                         {
 543                             rPInvC = one / (alphaCoulEff * lFacCoul[i] * sigma6[i] + rp);
 544                             pthRoot(rPInvC, &rInvC, &rC);
 545                             if (scLambdasOrAlphasDiffer)
 546                             {
 547                                 rPInvV = one / (alphaVdwEff * lFacVdw[i] * sigma6[i] + rp);
 548                                 pthRoot(rPInvV, &rInvV, &rV);
 549                             }
 550                             else
 551                             {
 552                                 /* We can avoid one expensive pow and one / operation */
 553                                 rPInvV = rPInvC;
 554                                 rInvV  = rInvC;
 555                                 rV     = rC;
 556                             }
 557                         }
 558                         else
 559                         {
 560                             rPInvC = 1;
 561                             rInvC  = rInv;
 562                             rC     = r;
 563
 564                             rPInvV = 1;
 565                             rInvV  = rInv;
 566                             rV     = r;
 567                         }
 568
 569                         /* Only process the coulomb interactions if we have charges,
 570                          * and if we either include all entries in the list (no cutoff
 571                          * used in the kernel), or if we are within the cutoff.
 572                          */
 573                         bool computeElecInteraction = (elecInteractionTypeIsEwald && r < rCoulomb)
 574                                                       || (!elecInteractionTypeIsEwald && rC < rCoulomb);
 575
 576                         if ((qq[i] != 0) && computeElecInteraction)
 577                         {
 578                             if (elecInteractionTypeIsEwald)
 579                             {
 580                                 vCoul[i]  = ewaldPotential(qq[i], rInvC, sh_ewald);
 581                                 fScalC[i] = ewaldScalarForce(qq[i], rInvC);
 582                             }
 583                             else
 584                             {
 585                                 vCoul[i]  = reactionFieldPotential(qq[i], rInvC, rC, krf, crf);
 586                                 fScalC[i] = reactionFieldScalarForce(qq[i], rInvC, rC, krf, two);
 587                             }
 588                         }
 589
 590                         /* Only process the VDW interactions if we have
 591                          * some non-zero parameters, and if we either
 592                          * include all entries in the list (no cutoff used
 593                          * in the kernel), or if we are within the cutoff.
 594                          */
 595                         bool computeVdwInteraction = (vdwInteractionTypeIsEwald && r < rVdw)
 596                                                      || (!vdwInteractionTypeIsEwald && rV < rVdw);
 597                         if ((c6[i] != 0 || c12[i] != 0) && computeVdwInteraction)
 598                         {
 599                             RealType rInv6;
 600                             if (useSoftCore)
 601                             {
 602                                 rInv6 = rPInvV;
 603                             }
 604                             else
 605                             {
 606                                 rInv6 = calculateRinv6(rInvV);
 607                             }
 608                             RealType vVdw6  = calculateVdw6(c6[i], rInv6);
 609                             RealType vVdw12 = calculateVdw12(c12[i], rInv6);
 610
 611                             vVdw[i] = lennardJonesPotential(
 612                                     vVdw6, vVdw12, c6[i], c12[i], repulsionShift, dispersionShift, oneSixth, oneTwelfth);
 613                             fScalV[i] = lennardJonesScalarForce(vVdw6, vVdw12);
 614
 615                             if (vdwInteractionTypeIsEwald)
 616                             {
 617                                 /* Subtract the grid potential at the cut-off */
 618                                 vVdw[i] += ewaldLennardJonesGridSubtract(
 619                                         nbfp_grid[tj[i]], shLjEwald, oneSixth);
 620                             }
 621
 622                             if (vdwModifierIsPotSwitch)
 623                             {
 624                                 RealType d        = rV - ic->rvdw_switch;
 625                                 d                 = (d > zero) ? d : zero;
 626                                 const RealType d2 = d * d;
 627                                 const RealType sw =
 628                                         one + d2 * d * (vdw_swV3 + d * (vdw_swV4 + d * vdw_swV5));
 629                                 const RealType dsw = d2 * (vdw_swF2 + d * (vdw_swF3 + d * vdw_swF4));
 630
 631                                 fScalV[i] = potSwitchScalarForceMod(
 632                                         fScalV[i], vVdw[i], sw, rV, rVdw, dsw, zero);
 633                                 vVdw[i] = potSwitchPotentialMod(vVdw[i], sw, rV, rVdw, zero);
 634                             }
 635                         }
 636
 637                         /* fScalC (and fScalV) now contain: dV/drC * rC
 638                          * Now we multiply by rC^-p, so it will be: dV/drC * rC^1-p
 639                          * Further down we first multiply by r^p-2 and then by
 640                          * the vector r, which in total gives: dV/drC * (r/rC)^1-p
 641                          */
 642                         fScalC[i] *= rPInvC;
 643                         fScalV[i] *= rPInvV;
 644                     }
 645                 } // end for (int i = 0; i < NSTATES; i++)
 646
 647                 /* Assemble A and B states */
 648                 for (int i = 0; i < NSTATES; i++)
 649                 {
 650                     vCTot += LFC[i] * vCoul[i];
 651                     vVTot += LFV[i] * vVdw[i];
 652
 653                     fScal += LFC[i] * fScalC[i] * rpm2;
 654                     fScal += LFV[i] * fScalV[i] * rpm2;
 655
 656                     if (useSoftCore)
 657                     {
 658                         dvdlCoul += vCoul[i] * DLF[i]
 659                                     + LFC[i] * alphaCoulEff * dlFacCoul[i] * fScalC[i] * sigma6[i];
 660                         dvdlVdw += vVdw[i] * DLF[i]
 661                                    + LFV[i] * alphaVdwEff * dlFacVdw[i] * fScalV[i] * sigma6[i];
 662                     }
 663                     else
 664                     {
 665                         dvdlCoul += vCoul[i] * DLF[i];
 666                         dvdlVdw += vVdw[i] * DLF[i];
 667                     }
 668                 }
 669             } // end if (bPairIncluded)
 670             else if (icoul == GMX_NBKERNEL_ELEC_REACTIONFIELD)
 671             {
 672                 /* For excluded pairs, which are only in this pair list when
 673                  * using the Verlet scheme, we don't use soft-core.
 674                  * As there is no singularity, there is no need for soft-core.
 675                  */
 676                 const real FF = -two * krf;
 677                 RealType   VV = krf * rSq - crf;
 678
 679                 if (ii == jnr)
 680                 {
 681                     VV *= half;
 682                 }
 683
 684                 for (int i = 0; i < NSTATES; i++)
 685                 {
 686                     vCTot += LFC[i] * qq[i] * VV;
 687                     fScal += LFC[i] * qq[i] * FF;
 688                     dvdlCoul += DLF[i] * qq[i] * VV;
 689                 }
 690             }
 691
 692             if (elecInteractionTypeIsEwald && (r < rCoulomb || !bPairIncluded))
 693             {
 694                 /* See comment in the preamble. When using Ewald interactions
 695                  * (unless we use a switch modifier) we subtract the reciprocal-space
 696                  * Ewald component here which made it possible to apply the free
 697                  * energy interaction to 1/r (vanilla coulomb short-range part)
 698                  * above. This gets us closer to the ideal case of applying
 699                  * the softcore to the entire electrostatic interaction,
 700                  * including the reciprocal-space component.
 701                  */
 702                 real v_lr, f_lr;
 703
 704                 const RealType ewrt   = r * coulombTableScale;
 705                 IntType        ewitab = static_cast<IntType>(ewrt);
 706                 const RealType eweps  = ewrt - ewitab;
 707                 ewitab                = 4 * ewitab;
 708                 f_lr                  = ewtab[ewitab] + eweps * ewtab[ewitab + 1];
 709                 v_lr = (ewtab[ewitab + 2] - coulombTableScaleInvHalf * eweps * (ewtab[ewitab] + f_lr));
 710                 f_lr *= rInv;
 711
 712                 /* Note that any possible Ewald shift has already been applied in
 713                  * the normal interaction part above.
 714                  */
 715
 716                 if (ii == jnr)
 717                 {
 718                     /* If we get here, the i particle (ii) has itself (jnr)
 719                      * in its neighborlist. This can only happen with the Verlet
 720                      * scheme, and corresponds to a self-interaction that will
 721                      * occur twice. Scale it down by 50% to only include it once.
 722                      */
 723                     v_lr *= half;
 724                 }
 725
 726                 for (int i = 0; i < NSTATES; i++)
 727                 {
 728                     vCTot -= LFC[i] * qq[i] * v_lr;
 729                     fScal -= LFC[i] * qq[i] * f_lr;
 730                     dvdlCoul -= (DLF[i] * qq[i]) * v_lr;
 731                 }
 732             }
 733
 734             if (vdwInteractionTypeIsEwald && (r < rVdw || !bPairIncluded))
 735             {
 736                 /* See comment in the preamble. When using LJ-Ewald interactions
 737                  * (unless we use a switch modifier) we subtract the reciprocal-space
 738                  * Ewald component here which made it possible to apply the free
 739                  * energy interaction to r^-6 (vanilla LJ6 short-range part)
 740                  * above. This gets us closer to the ideal case of applying
 741                  * the softcore to the entire VdW interaction,
 742                  * including the reciprocal-space component.
 743                  */
 744                 /* We could also use the analytical form here
 745                  * iso a table, but that can cause issues for
 746                  * r close to 0 for non-interacting pairs.
 747                  */
 748
 749                 const RealType rs   = rSq * rInv * vdwTableScale;
 750                 const IntType  ri   = static_cast<IntType>(rs);
 751                 const RealType frac = rs - ri;
 752                 const RealType f_lr = (1 - frac) * tab_ewald_F_lj[ri] + frac * tab_ewald_F_lj[ri + 1];
 753                 /* TODO: Currently the Ewald LJ table does not contain
 754                  * the factor 1/6, we should add this.
 755                  */
 756                 const RealType FF = f_lr * rInv / six;
 757                 RealType       VV =
 758                         (tab_ewald_V_lj[ri] - vdwTableScaleInvHalf * frac * (tab_ewald_F_lj[ri] + f_lr))
 759                         / six;
 760
 761                 if (ii == jnr)
 762                 {
 763                     /* If we get here, the i particle (ii) has itself (jnr)
 764                      * in its neighborlist. This can only happen with the Verlet
 765                      * scheme, and corresponds to a self-interaction that will
 766                      * occur twice. Scale it down by 50% to only include it once.
 767                      */
 768                     VV *= half;
 769                 }
 770
 771                 for (int i = 0; i < NSTATES; i++)
 772                 {
 773                     const real c6grid = nbfp_grid[tj[i]];
 774                     vVTot += LFV[i] * c6grid * VV;
 775                     fScal += LFV[i] * c6grid * FF;
 776                     dvdlVdw += (DLF[i] * c6grid) * VV;
 777                 }
 778             }
 779
 780             if (doForces)
 781             {
 782                 const real tX = fScal * dX;
 783                 const real tY = fScal * dY;
 784                 const real tZ = fScal * dZ;
 785                 fIX           = fIX + tX;
 786                 fIY           = fIY + tY;
 787                 fIZ           = fIZ + tZ;
 788                 /* OpenMP atomics are expensive, but this kernels is also
 789                  * expensive, so we can take this hit, instead of using
 790                  * thread-local output buffers and extra reduction.
 791                  *
 792                  * All the OpenMP regions in this file are trivial and should
 793                  * not throw, so no need for try/catch.
 794                  */
 795 #pragma omp atomic
 796                 f[j3] -= tX;
 797 #pragma omp atomic
 798                 f[j3 + 1] -= tY;
 799 #pragma omp atomic
 800                 f[j3 + 2] -= tZ;
 801             }
 802         } // end for (int k = nj0; k < nj1; k++)
 803
 804         /* The atomics below are expensive with many OpenMP threads.
 805          * Here unperturbed i-particles will usually only have a few
 806          * (perturbed) j-particles in the list. Thus with a buffered list
 807          * we can skip a significant number of i-reductions with a check.
 808          */
 809         if (npair_within_cutoff > 0)
 810         {
 811             if (doForces)
 812             {
 813 #pragma omp atomic
 814                 f[ii3] += fIX;
 815 #pragma omp atomic
 816                 f[ii3 + 1] += fIY;
 817 #pragma omp atomic
 818                 f[ii3 + 2] += fIZ;
 819             }
 820             if (doShiftForces)
 821             {
 822 #pragma omp atomic
 823                 fshift[is3] += fIX;
 824 #pragma omp atomic
 825                 fshift[is3 + 1] += fIY;
 826 #pragma omp atomic
 827                 fshift[is3 + 2] += fIZ;
 828             }
 829             if (doPotential)
 830             {
 831                 int ggid = gid[n];
 832 #pragma omp atomic
 833                 Vc[ggid] += vCTot;
 834 #pragma omp atomic
 835                 Vv[ggid] += vVTot;
 836             }
 837         }
 838     } // end for (int n = 0; n < nri; n++)
 839
 840 #pragma omp atomic
 841     dvdl[efptCOUL] += dvdlCoul;
 842 #pragma omp atomic
 843     dvdl[efptVDW] += dvdlVdw;
 844
 845     /* Estimate flops, average for free energy stuff:
 846      * 12  flops per outer iteration
 847      * 150 flops per inner iteration
 848      */
 849 #pragma omp atomic
 850     inc_nrnb(nrnb, eNR_NBKERNEL_FREE_ENERGY, nlist->nri * 12 + nlist->jindex[nri] * 150);
 851
 852     if (numExcludedPairsBeyondRlist > 0)
 853     {
 854         gmx_fatal(FARGS,
 855                   "There are %d perturbed non-bonded pair interactions beyond the pair-list cutoff "
 856                   "of %g nm, which is not supported. This can happen because the system is "
 857                   "unstable or because intra-molecular interactions at long distances are "
 858                   "excluded. If the "
 859                   "latter is the case, you can try to increase nstlist or rlist to avoid this."
 860                   "The error is likely triggered by the use of couple-intramol=no "
 861                   "and the maximal distance in the decoupled molecule exceeding rlist.",
 862                   numExcludedPairsBeyondRlist,
 863                   fr->rlist);
 864     }
 865 }
 866
 867 typedef void (*KernelFunction)(const t_nblist* gmx_restrict nlist,
 868                                rvec* gmx_restrict         xx,
 869                                gmx::ForceWithShiftForces* forceWithShiftForces,
 870                                const t_forcerec* gmx_restrict fr,
 871                                const t_mdatoms* gmx_restrict mdatoms,
 872                                nb_kernel_data_t* gmx_restrict kernel_data,
 873                                t_nrnb* gmx_restrict nrnb);
 874
 875 template<bool useSoftCore, bool scLambdasOrAlphasDiffer, bool vdwInteractionTypeIsEwald, bool elecInteractionTypeIsEwald, bool vdwModifierIsPotSwitch>
 876 static KernelFunction dispatchKernelOnUseSimd(const bool useSimd)
 877 {
 878     if (useSimd)
 879     {
 880 #if GMX_SIMD_HAVE_REAL && GMX_SIMD_HAVE_INT32_ARITHMETICS && GMX_USE_SIMD_KERNELS
 881         /* FIXME: Here SimdDataTypes should be used to enable SIMD. So far, the code in nb_free_energy_kernel is not adapted to SIMD */
 882         return (nb_free_energy_kernel<ScalarDataTypes, useSoftCore, scLambdasOrAlphasDiffer, vdwInteractionTypeIsEwald, elecInteractionTypeIsEwald, vdwModifierIsPotSwitch>);
 883 #else
 884         return (nb_free_energy_kernel<ScalarDataTypes, useSoftCore, scLambdasOrAlphasDiffer, vdwInteractionTypeIsEwald, elecInteractionTypeIsEwald, vdwModifierIsPotSwitch>);
 885 #endif
 886     }
 887     else
 888     {
 889         return (nb_free_energy_kernel<ScalarDataTypes, useSoftCore, scLambdasOrAlphasDiffer, vdwInteractionTypeIsEwald, elecInteractionTypeIsEwald, vdwModifierIsPotSwitch>);
 890     }
 891 }
 892
 893 template<bool useSoftCore, bool scLambdasOrAlphasDiffer, bool vdwInteractionTypeIsEwald, bool elecInteractionTypeIsEwald>
 894 static KernelFunction dispatchKernelOnVdwModifier(const bool vdwModifierIsPotSwitch, const bool useSimd)
 895 {
 896     if (vdwModifierIsPotSwitch)
 897     {
 898         return (dispatchKernelOnUseSimd<useSoftCore, scLambdasOrAlphasDiffer, vdwInteractionTypeIsEwald, elecInteractionTypeIsEwald, true>(
 899                 useSimd));
 900     }
 901     else
 902     {
 903         return (dispatchKernelOnUseSimd<useSoftCore, scLambdasOrAlphasDiffer, vdwInteractionTypeIsEwald, elecInteractionTypeIsEwald, false>(
 904                 useSimd));
 905     }
 906 }
 907
 908 template<bool useSoftCore, bool scLambdasOrAlphasDiffer, bool vdwInteractionTypeIsEwald>
 909 static KernelFunction dispatchKernelOnElecInteractionType(const bool elecInteractionTypeIsEwald,
 910                                                           const bool vdwModifierIsPotSwitch,
 911                                                           const bool useSimd)
 912 {
 913     if (elecInteractionTypeIsEwald)
 914     {
 915         return (dispatchKernelOnVdwModifier<useSoftCore, scLambdasOrAlphasDiffer, vdwInteractionTypeIsEwald, true>(
 916                 vdwModifierIsPotSwitch, useSimd));
 917     }
 918     else
 919     {
 920         return (dispatchKernelOnVdwModifier<useSoftCore, scLambdasOrAlphasDiffer, vdwInteractionTypeIsEwald, false>(
 921                 vdwModifierIsPotSwitch, useSimd));
 922     }
 923 }
 924
 925 template<bool useSoftCore, bool scLambdasOrAlphasDiffer>
 926 static KernelFunction dispatchKernelOnVdwInteractionType(const bool vdwInteractionTypeIsEwald,
 927                                                          const bool elecInteractionTypeIsEwald,
 928                                                          const bool vdwModifierIsPotSwitch,
 929                                                          const bool useSimd)
 930 {
 931     if (vdwInteractionTypeIsEwald)
 932     {
 933         return (dispatchKernelOnElecInteractionType<useSoftCore, scLambdasOrAlphasDiffer, true>(
 934                 elecInteractionTypeIsEwald, vdwModifierIsPotSwitch, useSimd));
 935     }
 936     else
 937     {
 938         return (dispatchKernelOnElecInteractionType<useSoftCore, scLambdasOrAlphasDiffer, false>(
 939                 elecInteractionTypeIsEwald, vdwModifierIsPotSwitch, useSimd));
 940     }
 941 }
 942
 943 template<bool useSoftCore>
 944 static KernelFunction dispatchKernelOnScLambdasOrAlphasDifference(const bool scLambdasOrAlphasDiffer,
 945                                                                   const bool vdwInteractionTypeIsEwald,
 946                                                                   const bool elecInteractionTypeIsEwald,
 947                                                                   const bool vdwModifierIsPotSwitch,
 948                                                                   const bool useSimd)
 949 {
 950     if (scLambdasOrAlphasDiffer)
 951     {
 952         return (dispatchKernelOnVdwInteractionType<useSoftCore, true>(
 953                 vdwInteractionTypeIsEwald, elecInteractionTypeIsEwald, vdwModifierIsPotSwitch, useSimd));
 954     }
 955     else
 956     {
 957         return (dispatchKernelOnVdwInteractionType<useSoftCore, false>(
 958                 vdwInteractionTypeIsEwald, elecInteractionTypeIsEwald, vdwModifierIsPotSwitch, useSimd));
 959     }
 960 }
 961
 962 static KernelFunction dispatchKernel(const bool                 scLambdasOrAlphasDiffer,
 963                                      const bool                 vdwInteractionTypeIsEwald,
 964                                      const bool                 elecInteractionTypeIsEwald,
 965                                      const bool                 vdwModifierIsPotSwitch,
 966                                      const bool                 useSimd,
 967                                      const interaction_const_t& ic)
 968 {
 969     if (ic.softCoreParameters->alphaCoulomb == 0 && ic.softCoreParameters->alphaVdw == 0)
 970     {
 971         return (dispatchKernelOnScLambdasOrAlphasDifference<false>(scLambdasOrAlphasDiffer,
 972                                                                    vdwInteractionTypeIsEwald,
 973                                                                    elecInteractionTypeIsEwald,
 974                                                                    vdwModifierIsPotSwitch,
 975                                                                    useSimd));
 976     }
 977     else
 978     {
 979         return (dispatchKernelOnScLambdasOrAlphasDifference<true>(scLambdasOrAlphasDiffer,
 980                                                                   vdwInteractionTypeIsEwald,
 981                                                                   elecInteractionTypeIsEwald,
 982                                                                   vdwModifierIsPotSwitch,
 983                                                                   useSimd));
 984     }
 985 }
 986
 987
 988 void gmx_nb_free_energy_kernel(const t_nblist*            nlist,
 989                                rvec*                      xx,
 990                                gmx::ForceWithShiftForces* ff,
 991                                const t_forcerec*          fr,
 992                                const t_mdatoms*           mdatoms,
 993                                nb_kernel_data_t*          kernel_data,
 994                                t_nrnb*                    nrnb)
 995 {
 996     const interaction_const_t& ic = *fr->ic;
 997     GMX_ASSERT(EEL_PME_EWALD(ic.eeltype) || ic.eeltype == eelCUT || EEL_RF(ic.eeltype),
 998                "Unsupported eeltype with free energy");
 999     GMX_ASSERT(ic.softCoreParameters, "We need soft-core parameters");
1000
1001     const auto& scParams                   = *ic.softCoreParameters;
1002     const bool  vdwInteractionTypeIsEwald  = (EVDW_PME(ic.vdwtype));
1003     const bool  elecInteractionTypeIsEwald = (EEL_PME_EWALD(ic.eeltype));
1004     const bool  vdwModifierIsPotSwitch     = (ic.vdw_modifier == eintmodPOTSWITCH);
1005     bool        scLambdasOrAlphasDiffer    = true;
1006     const bool  useSimd                    = fr->use_simd_kernels;
1007
1008     if (scParams.alphaCoulomb == 0 && scParams.alphaVdw == 0)
1009     {
1010         scLambdasOrAlphasDiffer = false;
1011     }
1012     else
1013     {
1014         if (kernel_data->lambda[efptCOUL] == kernel_data->lambda[efptVDW]
1015             && scParams.alphaCoulomb == scParams.alphaVdw)
1016         {
1017             scLambdasOrAlphasDiffer = false;
1018         }
1019     }
1020
1021     KernelFunction kernelFunc;
1022     kernelFunc = dispatchKernel(scLambdasOrAlphasDiffer,
1023                                 vdwInteractionTypeIsEwald,
1024                                 elecInteractionTypeIsEwald,
1025                                 vdwModifierIsPotSwitch,
1026                                 useSimd,
1027                                 ic);
1028     kernelFunc(nlist, xx, ff, fr, mdatoms, kernel_data, nrnb);
1029 }