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,2018,2019, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include "nb_free_energy.h"
  40
  41 #include <cmath>
  42
  43 #include <algorithm>
  44
  45 #include "gromacs/gmxlib/nrnb.h"
  46 #include "gromacs/gmxlib/nonbonded/nb_kernel.h"
  47 #include "gromacs/gmxlib/nonbonded/nonbonded.h"
  48 #include "gromacs/math/functions.h"
  49 #include "gromacs/math/vec.h"
  50 #include "gromacs/mdtypes/forcerec.h"
  51 #include "gromacs/mdtypes/md_enums.h"
  52 #include "gromacs/utility/fatalerror.h"
  53
  54
  55 //! Enum for templating the soft-core treatment in the kernel
  56 enum class SoftCoreTreatment
  57 {
  58     None,    //!< No soft-core
  59     RPower6, //!< Soft-core with r-power = 6
  60     RPower48 //!< Soft-core with r-power = 48
  61 };
  62
  63 //! Most treatments are fine with float in mixed-precision mode.
  64 template <SoftCoreTreatment softCoreTreatment>
  65 struct SoftCoreReal
  66 {
  67     //! Real type for soft-core calculations
  68     using Real = real;
  69 };
  70
  71 //! This treatment requires double precision for some computations.
  72 template <>
  73 struct SoftCoreReal<SoftCoreTreatment::RPower48>
  74 {
  75     //! Real type for soft-core calculations
  76     using Real = double;
  77 };
  78
  79 //! Templated free-energy non-bonded kernel
  80 template<SoftCoreTreatment softCoreTreatment, bool scLambdasOrAlphasDiffer>
  81 static void
  82 nb_free_energy_kernel(const t_nblist * gmx_restrict    nlist,
  83                       rvec * gmx_restrict              xx,
  84                       rvec * gmx_restrict              ff,
  85                       t_forcerec * gmx_restrict        fr,
  86                       const t_mdatoms * gmx_restrict   mdatoms,
  87                       nb_kernel_data_t * gmx_restrict  kernel_data,
  88                       t_nrnb * gmx_restrict            nrnb)
  89 {
  90     using SCReal = typename SoftCoreReal<softCoreTreatment>::Real;
  91
  92     constexpr bool useSoftCore = (softCoreTreatment != SoftCoreTreatment::None);
  93
  94 #define  STATE_A  0
  95 #define  STATE_B  1
  96 #define  NSTATES  2
  97     int           i, n, ii, is3, ii3, k, nj0, nj1, jnr, j3, ggid;
  98     real          shX, shY, shZ;
  99     real          tx, ty, tz, Fscal;
 100     SCReal        FscalC[NSTATES], FscalV[NSTATES];  /* Needs double for sc_power==48 */
 101     SCReal        Vcoul[NSTATES], Vvdw[NSTATES];     /* Needs double for sc_power==48 */
 102     real          rinv6, r;
 103     real          iqA, iqB;
 104     real          qq[NSTATES], vctot;
 105     int           ntiA, ntiB, tj[NSTATES];
 106     real          Vvdw6, Vvdw12, vvtot;
 107     real          ix, iy, iz, fix, fiy, fiz;
 108     real          dx, dy, dz, rsq, rinv;
 109     real          c6[NSTATES], c12[NSTATES], c6grid;
 110     real          LFC[NSTATES], LFV[NSTATES], DLF[NSTATES];
 111     SCReal        dvdl_coul, dvdl_vdw;
 112     real          lfac_coul[NSTATES], dlfac_coul[NSTATES], lfac_vdw[NSTATES], dlfac_vdw[NSTATES];
 113     real          sigma6[NSTATES], alpha_vdw_eff, alpha_coul_eff;
 114     SCReal        rp, rpm2, rC, rV, rinvC, rpinvC, rinvV, rpinvV; /* Needs double for sc_power==48 */
 115     real          sigma_pow[NSTATES];
 116     real          VV, FF;
 117     int           icoul, ivdw;
 118     int           nri;
 119     const int *   iinr;
 120     const int *   jindex;
 121     const int *   jjnr;
 122     const int *   shift;
 123     const int *   gid;
 124     const int *   typeA;
 125     const int *   typeB;
 126     int           ntype;
 127     const real *  shiftvec;
 128     real *        fshift;
 129     const real *  x;
 130     real *        f;
 131     const real *  chargeA;
 132     const real *  chargeB;
 133     real          sigma6_min, sigma6_def, lam_power;
 134     real          alpha_coul, alpha_vdw, lambda_coul, lambda_vdw;
 135     const real *  nbfp, *nbfp_grid;
 136     real *        dvdl;
 137     real *        Vv;
 138     real *        Vc;
 139     gmx_bool      bDoForces, bDoShiftForces, bDoPotential;
 140     gmx_bool      bEwald, bEwaldLJ;
 141     real          rcutoff_max2;
 142     const real *  tab_ewald_F_lj = nullptr;
 143     const real *  tab_ewald_V_lj = nullptr;
 144     real          d, d2, sw, dsw;
 145     real          vdw_swV3, vdw_swV4, vdw_swV5, vdw_swF2, vdw_swF3, vdw_swF4;
 146     gmx_bool      bComputeVdwInteraction, bComputeElecInteraction;
 147     const real *  ewtab = nullptr;
 148     int           ewitab;
 149     real          ewrt, eweps, ewtabscale = 0, ewtabhalfspace = 0, sh_ewald = 0;
 150
 151     const real    onetwelfth  = 1.0/12.0;
 152     const real    onesixth    = 1.0/6.0;
 153     const real    zero        = 0.0;
 154     const real    half        = 0.5;
 155     const real    one         = 1.0;
 156     const real    two         = 2.0;
 157     const real    six         = 6.0;
 158
 159     /* Extract pointer to non-bonded interaction constants */
 160     const interaction_const_t *ic = fr->ic;
 161
 162     x                   = xx[0];
 163     f                   = ff[0];
 164
 165     fshift              = fr->fshift[0];
 166
 167     nri                 = nlist->nri;
 168     iinr                = nlist->iinr;
 169     jindex              = nlist->jindex;
 170     jjnr                = nlist->jjnr;
 171     shift               = nlist->shift;
 172     gid                 = nlist->gid;
 173
 174     shiftvec            = fr->shift_vec[0];
 175     chargeA             = mdatoms->chargeA;
 176     chargeB             = mdatoms->chargeB;
 177     Vc                  = kernel_data->energygrp_elec;
 178     typeA               = mdatoms->typeA;
 179     typeB               = mdatoms->typeB;
 180     ntype               = fr->ntype;
 181     nbfp                = fr->nbfp;
 182     nbfp_grid           = fr->ljpme_c6grid;
 183     Vv                  = kernel_data->energygrp_vdw;
 184     lambda_coul         = kernel_data->lambda[efptCOUL];
 185     lambda_vdw          = kernel_data->lambda[efptVDW];
 186     dvdl                = kernel_data->dvdl;
 187     alpha_coul          = fr->sc_alphacoul;
 188     alpha_vdw           = fr->sc_alphavdw;
 189     lam_power           = fr->sc_power;
 190     sigma6_def          = fr->sc_sigma6_def;
 191     sigma6_min          = fr->sc_sigma6_min;
 192     bDoForces           = ((kernel_data->flags & GMX_NONBONDED_DO_FORCE) != 0);
 193     bDoShiftForces      = ((kernel_data->flags & GMX_NONBONDED_DO_SHIFTFORCE) != 0);
 194     bDoPotential        = ((kernel_data->flags & GMX_NONBONDED_DO_POTENTIAL) != 0);
 195
 196     // Extract data from interaction_const_t
 197     const real facel           = ic->epsfac;
 198     const real rcoulomb        = ic->rcoulomb;
 199     const real krf             = ic->k_rf;
 200     const real crf             = ic->c_rf;
 201     const real sh_lj_ewald     = ic->sh_lj_ewald;
 202     const real rvdw            = ic->rvdw;
 203     const real dispersionShift = ic->dispersion_shift.cpot;
 204     const real repulsionShift  = ic->repulsion_shift.cpot;
 205
 206     // Note that the nbnxm kernels do not support Coulomb potential switching at all
 207     GMX_ASSERT(ic->coulomb_modifier != eintmodPOTSWITCH,
 208                "Potential switching is not supported for Coulomb with FEP");
 209
 210     if (ic->vdw_modifier == eintmodPOTSWITCH)
 211     {
 212         d               = ic->rvdw - ic->rvdw_switch;
 213         vdw_swV3        = -10.0/(d*d*d);
 214         vdw_swV4        =  15.0/(d*d*d*d);
 215         vdw_swV5        =  -6.0/(d*d*d*d*d);
 216         vdw_swF2        = -30.0/(d*d*d);
 217         vdw_swF3        =  60.0/(d*d*d*d);
 218         vdw_swF4        = -30.0/(d*d*d*d*d);
 219     }
 220     else
 221     {
 222         /* Avoid warnings from stupid compilers (looking at you, Clang!) */
 223         vdw_swV3 = vdw_swV4 = vdw_swV5 = vdw_swF2 = vdw_swF3 = vdw_swF4 = 0.0;
 224     }
 225
 226     if (EVDW_PME(ic->vdwtype))
 227     {
 228         ivdw         = GMX_NBKERNEL_VDW_LJEWALD;
 229     }
 230     else
 231     {
 232         ivdw         = GMX_NBKERNEL_VDW_LENNARDJONES;
 233     }
 234
 235     if (ic->eeltype == eelCUT || EEL_RF(ic->eeltype))
 236     {
 237         icoul        = GMX_NBKERNEL_ELEC_REACTIONFIELD;
 238     }
 239     else if (EEL_PME_EWALD(ic->eeltype))
 240     {
 241         icoul        = GMX_NBKERNEL_ELEC_EWALD;
 242     }
 243     else
 244     {
 245         gmx_incons("Unsupported eeltype with Verlet and free-energy");
 246     }
 247
 248     rcutoff_max2    = std::max(ic->rcoulomb, ic->rvdw);
 249     rcutoff_max2    = rcutoff_max2*rcutoff_max2;
 250
 251     bEwald          = (icoul == GMX_NBKERNEL_ELEC_EWALD);
 252     bEwaldLJ        = (ivdw == GMX_NBKERNEL_VDW_LJEWALD);
 253
 254     if (bEwald || bEwaldLJ)
 255     {
 256         sh_ewald       = ic->sh_ewald;
 257         ewtab          = ic->tabq_coul_FDV0;
 258         ewtabscale     = ic->tabq_scale;
 259         ewtabhalfspace = half/ewtabscale;
 260         tab_ewald_F_lj = ic->tabq_vdw_F;
 261         tab_ewald_V_lj = ic->tabq_vdw_V;
 262     }
 263
 264     /* For Ewald/PME interactions we cannot easily apply the soft-core component to
 265      * reciprocal space. When we use non-switched Ewald interactions, we
 266      * assume the soft-coring does not significantly affect the grid contribution
 267      * and apply the soft-core only to the full 1/r (- shift) pair contribution.
 268      *
 269      * However, we cannot use this approach for switch-modified since we would then
 270      * effectively end up evaluating a significantly different interaction here compared to the
 271      * normal (non-free-energy) kernels, either by applying a cutoff at a different
 272      * position than what the user requested, or by switching different
 273      * things (1/r rather than short-range Ewald). For these settings, we just
 274      * use the traditional short-range Ewald interaction in that case.
 275      */
 276     GMX_RELEASE_ASSERT(!(bEwald   && ic->coulomb_modifier == eintmodPOTSWITCH) &&
 277                        !(bEwaldLJ && ic->vdw_modifier   == eintmodPOTSWITCH),
 278                        "Can not apply soft-core to switched Ewald potentials");
 279
 280     dvdl_coul  = 0;
 281     dvdl_vdw   = 0;
 282
 283     /* Lambda factor for state A, 1-lambda*/
 284     LFC[STATE_A] = one - lambda_coul;
 285     LFV[STATE_A] = one - lambda_vdw;
 286
 287     /* Lambda factor for state B, lambda*/
 288     LFC[STATE_B] = lambda_coul;
 289     LFV[STATE_B] = lambda_vdw;
 290
 291     /*derivative of the lambda factor for state A and B */
 292     DLF[STATE_A] = -1;
 293     DLF[STATE_B] = 1;
 294
 295     constexpr real sc_r_power = (softCoreTreatment == SoftCoreTreatment::RPower48 ? 48.0_real : 6.0_real);
 296     for (i = 0; i < NSTATES; i++)
 297     {
 298         lfac_coul[i]  = (lam_power == 2 ? (1-LFC[i])*(1-LFC[i]) : (1-LFC[i]));
 299         dlfac_coul[i] = DLF[i]*lam_power/sc_r_power*(lam_power == 2 ? (1-LFC[i]) : 1);
 300         lfac_vdw[i]   = (lam_power == 2 ? (1-LFV[i])*(1-LFV[i]) : (1-LFV[i]));
 301         dlfac_vdw[i]  = DLF[i]*lam_power/sc_r_power*(lam_power == 2 ? (1-LFV[i]) : 1);
 302     }
 303
 304     for (n = 0; (n < nri); n++)
 305     {
 306         int npair_within_cutoff;
 307
 308         npair_within_cutoff = 0;
 309
 310         is3              = 3*shift[n];
 311         shX              = shiftvec[is3];
 312         shY              = shiftvec[is3+1];
 313         shZ              = shiftvec[is3+2];
 314         nj0              = jindex[n];
 315         nj1              = jindex[n+1];
 316         ii               = iinr[n];
 317         ii3              = 3*ii;
 318         ix               = shX + x[ii3+0];
 319         iy               = shY + x[ii3+1];
 320         iz               = shZ + x[ii3+2];
 321         iqA              = facel*chargeA[ii];
 322         iqB              = facel*chargeB[ii];
 323         ntiA             = 2*ntype*typeA[ii];
 324         ntiB             = 2*ntype*typeB[ii];
 325         vctot            = 0;
 326         vvtot            = 0;
 327         fix              = 0;
 328         fiy              = 0;
 329         fiz              = 0;
 330
 331         for (k = nj0; (k < nj1); k++)
 332         {
 333             jnr              = jjnr[k];
 334             j3               = 3*jnr;
 335             dx               = ix - x[j3];
 336             dy               = iy - x[j3+1];
 337             dz               = iz - x[j3+2];
 338             rsq              = dx*dx + dy*dy + dz*dz;
 339
 340             if (rsq >= rcutoff_max2)
 341             {
 342                 /* We save significant time by skipping all code below.
 343                  * Note that with soft-core interactions, the actual cut-off
 344                  * check might be different. But since the soft-core distance
 345                  * is always larger than r, checking on r here is safe.
 346                  */
 347                 continue;
 348             }
 349             npair_within_cutoff++;
 350
 351             if (rsq > 0)
 352             {
 353                 /* Note that unlike in the nbnxn kernels, we do not need
 354                  * to clamp the value of rsq before taking the invsqrt
 355                  * to avoid NaN in the LJ calculation, since here we do
 356                  * not calculate LJ interactions when C6 and C12 are zero.
 357                  */
 358
 359                 rinv         = gmx::invsqrt(rsq);
 360                 r            = rsq*rinv;
 361             }
 362             else
 363             {
 364                 /* The force at r=0 is zero, because of symmetry.
 365                  * But note that the potential is in general non-zero,
 366                  * since the soft-cored r will be non-zero.
 367                  */
 368                 rinv         = 0;
 369                 r            = 0;
 370             }
 371
 372             if (softCoreTreatment == SoftCoreTreatment::None)
 373             {
 374                 /* The soft-core power p will not affect the results
 375                  * with not using soft-core, so we use power of 0 which gives
 376                  * the simplest math and cheapest code.
 377                  */
 378                 rpm2         = rinv*rinv;
 379                 rp           = 1;
 380             }
 381             if (softCoreTreatment == SoftCoreTreatment::RPower6)
 382             {
 383                 rpm2         = rsq*rsq;  /* r4 */
 384                 rp           = rpm2*rsq; /* r6 */
 385             }
 386             if (softCoreTreatment == SoftCoreTreatment::RPower48)
 387             {
 388                 rp           = rsq*rsq*rsq; /* r6 */
 389                 rp           = rp*rp;       /* r12 */
 390                 rp           = rp*rp;       /* r24 */
 391                 rp           = rp*rp;       /* r48 */
 392                 rpm2         = rp/rsq;      /* r46 */
 393             }
 394
 395             Fscal = 0;
 396
 397             qq[STATE_A]      = iqA*chargeA[jnr];
 398             qq[STATE_B]      = iqB*chargeB[jnr];
 399
 400             tj[STATE_A]      = ntiA+2*typeA[jnr];
 401             tj[STATE_B]      = ntiB+2*typeB[jnr];
 402
 403             if (nlist->excl_fep == nullptr || nlist->excl_fep[k])
 404             {
 405                 c6[STATE_A]      = nbfp[tj[STATE_A]];
 406                 c6[STATE_B]      = nbfp[tj[STATE_B]];
 407
 408                 for (i = 0; i < NSTATES; i++)
 409                 {
 410                     c12[i]             = nbfp[tj[i]+1];
 411                     if (useSoftCore)
 412                     {
 413                         if ((c6[i] > 0) && (c12[i] > 0))
 414                         {
 415                             /* c12 is stored scaled with 12.0 and c6 is scaled with 6.0 - correct for this */
 416                             sigma6[i]       = half*c12[i]/c6[i];
 417                             if (sigma6[i] < sigma6_min)   /* for disappearing coul and vdw with soft core at the same time */
 418                             {
 419                                 sigma6[i] = sigma6_min;
 420                             }
 421                         }
 422                         else
 423                         {
 424                             sigma6[i]       = sigma6_def;
 425                         }
 426                         if (softCoreTreatment == SoftCoreTreatment::RPower6)
 427                         {
 428                             sigma_pow[i]    = sigma6[i];
 429                         }
 430                         else
 431                         {
 432                             sigma_pow[i]    = sigma6[i]*sigma6[i];       /* sigma^12 */
 433                             sigma_pow[i]    = sigma_pow[i]*sigma_pow[i]; /* sigma^24 */
 434                             sigma_pow[i]    = sigma_pow[i]*sigma_pow[i]; /* sigma^48 */
 435                         }
 436                     }
 437                 }
 438
 439                 if (useSoftCore)
 440                 {
 441                     /* only use softcore if one of the states has a zero endstate - softcore is for avoiding infinities!*/
 442                     if ((c12[STATE_A] > 0) && (c12[STATE_B] > 0))
 443                     {
 444                         alpha_vdw_eff    = 0;
 445                         alpha_coul_eff   = 0;
 446                     }
 447                     else
 448                     {
 449                         alpha_vdw_eff    = alpha_vdw;
 450                         alpha_coul_eff   = alpha_coul;
 451                     }
 452                 }
 453
 454                 for (i = 0; i < NSTATES; i++)
 455                 {
 456                     FscalC[i]    = 0;
 457                     FscalV[i]    = 0;
 458                     Vcoul[i]     = 0;
 459                     Vvdw[i]      = 0;
 460
 461                     /* Only spend time on A or B state if it is non-zero */
 462                     if ( (qq[i] != 0) || (c6[i] != 0) || (c12[i] != 0) )
 463                     {
 464                         /* this section has to be inside the loop because of the dependence on sigma_pow */
 465                         if (useSoftCore)
 466                         {
 467                             rpinvC     = one/(alpha_coul_eff*lfac_coul[i]*sigma_pow[i]+rp);
 468                             rinvC      = std::pow(rpinvC, one/sc_r_power);
 469                             rC         = one/rinvC;
 470
 471                             if (scLambdasOrAlphasDiffer)
 472                             {
 473                                 rpinvV = one/(alpha_vdw_eff*lfac_vdw[i]*sigma_pow[i]+rp);
 474                                 rinvV  = std::pow(rpinvV, one/sc_r_power);
 475                                 rV     = one/rinvV;
 476                             }
 477                             else
 478                             {
 479                                 /* We can avoid one expensive pow and one / operation */
 480                                 rpinvV = rpinvC;
 481                                 rinvV  = rinvC;
 482                                 rV     = rC;
 483                             }
 484                         }
 485                         else
 486                         {
 487                             rpinvC     = 1;
 488                             rinvC      = rinv;
 489                             rC         = r;
 490
 491                             rpinvV     = 1;
 492                             rinvV      = rinv;
 493                             rV         = r;
 494                         }
 495
 496                         /* Only process the coulomb interactions if we have charges,
 497                          * and if we either include all entries in the list (no cutoff
 498                          * used in the kernel), or if we are within the cutoff.
 499                          */
 500                         bComputeElecInteraction =
 501                             ( bEwald && r < rcoulomb) ||
 502                             (!bEwald && rC < rcoulomb);
 503
 504                         if ( (qq[i] != 0) && bComputeElecInteraction)
 505                         {
 506                             if (bEwald)
 507                             {
 508                                 /* Ewald FEP is done only on the 1/r part */
 509                                 Vcoul[i]   = qq[i]*(rinvC-sh_ewald);
 510                                 FscalC[i]  = qq[i]*rinvC;
 511                             }
 512                             else
 513                             {
 514                                 /* reaction-field */
 515                                 Vcoul[i]   = qq[i]*(rinvC + krf*rC*rC-crf);
 516                                 FscalC[i]  = qq[i]*(rinvC - two*krf*rC*rC);
 517                             }
 518                         }
 519
 520                         /* Only process the VDW interactions if we have
 521                          * some non-zero parameters, and if we either
 522                          * include all entries in the list (no cutoff used
 523                          * in the kernel), or if we are within the cutoff.
 524                          */
 525                         bComputeVdwInteraction =
 526                             ( bEwaldLJ && r < rvdw) ||
 527                             (!bEwaldLJ && rV < rvdw);
 528                         if ((c6[i] != 0 || c12[i] != 0) && bComputeVdwInteraction)
 529                         {
 530                             /* cutoff LJ, also handles part of Ewald LJ */
 531                             if (softCoreTreatment == SoftCoreTreatment::RPower6)
 532                             {
 533                                 rinv6        = rpinvV;
 534                             }
 535                             else
 536                             {
 537                                 rinv6        = rinvV*rinvV;
 538                                 rinv6        = rinv6*rinv6*rinv6;
 539                             }
 540                             Vvdw6            = c6[i]*rinv6;
 541                             Vvdw12           = c12[i]*rinv6*rinv6;
 542
 543                             Vvdw[i]          = ( (Vvdw12 + c12[i]*repulsionShift)*onetwelfth
 544                                                  - (Vvdw6 + c6[i]*dispersionShift)*onesixth);
 545                             FscalV[i]        = Vvdw12 - Vvdw6;
 546
 547                             if (bEwaldLJ)
 548                             {
 549                                 /* Subtract the grid potential at the cut-off */
 550                                 c6grid      = nbfp_grid[tj[i]];
 551                                 Vvdw[i]    += c6grid*sh_lj_ewald*onesixth;
 552                             }
 553
 554                             if (ic->vdw_modifier == eintmodPOTSWITCH)
 555                             {
 556                                 d                = rV - ic->rvdw_switch;
 557                                 d                = (d > zero) ? d : zero;
 558                                 d2               = d*d;
 559                                 sw               = one+d2*d*(vdw_swV3+d*(vdw_swV4+d*vdw_swV5));
 560                                 dsw              = d2*(vdw_swF2+d*(vdw_swF3+d*vdw_swF4));
 561
 562                                 FscalV[i]        = FscalV[i]*sw - rV*Vvdw[i]*dsw;
 563                                 Vvdw[i]         *= sw;
 564
 565                                 FscalV[i]  = (rV < rvdw) ? FscalV[i] : zero;
 566                                 Vvdw[i]    = (rV < rvdw) ? Vvdw[i] : zero;
 567                             }
 568                         }
 569
 570                         /* FscalC (and FscalV) now contain: dV/drC * rC
 571                          * Now we multiply by rC^-p, so it will be: dV/drC * rC^1-p
 572                          * Further down we first multiply by r^p-2 and then by
 573                          * the vector r, which in total gives: dV/drC * (r/rC)^1-p
 574                          */
 575                         FscalC[i] *= rpinvC;
 576                         FscalV[i] *= rpinvV;
 577                     }
 578                 }
 579
 580                 /* Assemble A and B states */
 581                 for (i = 0; i < NSTATES; i++)
 582                 {
 583                     vctot         += LFC[i]*Vcoul[i];
 584                     vvtot         += LFV[i]*Vvdw[i];
 585
 586                     Fscal     += LFC[i]*FscalC[i]*rpm2;
 587                     Fscal     += LFV[i]*FscalV[i]*rpm2;
 588
 589                     if (useSoftCore)
 590                     {
 591                         dvdl_coul += Vcoul[i]*DLF[i] + LFC[i]*alpha_coul_eff*dlfac_coul[i]*FscalC[i]*sigma_pow[i];
 592                         dvdl_vdw  += Vvdw[i]*DLF[i] + LFV[i]*alpha_vdw_eff*dlfac_vdw[i]*FscalV[i]*sigma_pow[i];
 593                     }
 594                     else
 595                     {
 596                         dvdl_coul += Vcoul[i]*DLF[i];
 597                         dvdl_vdw  += Vvdw[i]*DLF[i];
 598                     }
 599                 }
 600             }
 601             else if (icoul == GMX_NBKERNEL_ELEC_REACTIONFIELD)
 602             {
 603                 /* For excluded pairs, which are only in this pair list when
 604                  * using the Verlet scheme, we don't use soft-core.
 605                  * The group scheme also doesn't soft-core for these.
 606                  * As there is no singularity, there is no need for soft-core.
 607                  */
 608                 VV = krf*rsq - crf;
 609                 FF = -two*krf;
 610
 611                 if (ii == jnr)
 612                 {
 613                     VV *= half;
 614                 }
 615
 616                 for (i = 0; i < NSTATES; i++)
 617                 {
 618                     vctot      += LFC[i]*qq[i]*VV;
 619                     Fscal      += LFC[i]*qq[i]*FF;
 620                     dvdl_coul  += DLF[i]*qq[i]*VV;
 621                 }
 622             }
 623
 624             if (bEwald && r < rcoulomb)
 625             {
 626                 /* See comment in the preamble. When using Ewald interactions
 627                  * (unless we use a switch modifier) we subtract the reciprocal-space
 628                  * Ewald component here which made it possible to apply the free
 629                  * energy interaction to 1/r (vanilla coulomb short-range part)
 630                  * above. This gets us closer to the ideal case of applying
 631                  * the softcore to the entire electrostatic interaction,
 632                  * including the reciprocal-space component.
 633                  */
 634                 real v_lr, f_lr;
 635
 636                 ewrt      = r*ewtabscale;
 637                 ewitab    = static_cast<int>(ewrt);
 638                 eweps     = ewrt-ewitab;
 639                 ewitab    = 4*ewitab;
 640                 f_lr      = ewtab[ewitab]+eweps*ewtab[ewitab+1];
 641                 v_lr      = (ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+f_lr));
 642                 f_lr     *= rinv;
 643
 644                 /* Note that any possible Ewald shift has already been applied in
 645                  * the normal interaction part above.
 646                  */
 647
 648                 if (ii == jnr)
 649                 {
 650                     /* If we get here, the i particle (ii) has itself (jnr)
 651                      * in its neighborlist. This can only happen with the Verlet
 652                      * scheme, and corresponds to a self-interaction that will
 653                      * occur twice. Scale it down by 50% to only include it once.
 654                      */
 655                     v_lr *= half;
 656                 }
 657
 658                 for (i = 0; i < NSTATES; i++)
 659                 {
 660                     vctot      -= LFC[i]*qq[i]*v_lr;
 661                     Fscal      -= LFC[i]*qq[i]*f_lr;
 662                     dvdl_coul  -= (DLF[i]*qq[i])*v_lr;
 663                 }
 664             }
 665
 666             if (bEwaldLJ && r < rvdw)
 667             {
 668                 /* See comment in the preamble. When using LJ-Ewald interactions
 669                  * (unless we use a switch modifier) we subtract the reciprocal-space
 670                  * Ewald component here which made it possible to apply the free
 671                  * energy interaction to r^-6 (vanilla LJ6 short-range part)
 672                  * above. This gets us closer to the ideal case of applying
 673                  * the softcore to the entire VdW interaction,
 674                  * including the reciprocal-space component.
 675                  */
 676                 /* We could also use the analytical form here
 677                  * iso a table, but that can cause issues for
 678                  * r close to 0 for non-interacting pairs.
 679                  */
 680                 real rs, frac, f_lr;
 681                 int  ri;
 682
 683                 rs     = rsq*rinv*ewtabscale;
 684                 ri     = static_cast<int>(rs);
 685                 frac   = rs - ri;
 686                 f_lr   = (1 - frac)*tab_ewald_F_lj[ri] + frac*tab_ewald_F_lj[ri+1];
 687                 /* TODO: Currently the Ewald LJ table does not contain
 688                  * the factor 1/6, we should add this.
 689                  */
 690                 FF     = f_lr*rinv/six;
 691                 VV     = (tab_ewald_V_lj[ri] - ewtabhalfspace*frac*(tab_ewald_F_lj[ri] + f_lr))/six;
 692
 693                 if (ii == jnr)
 694                 {
 695                     /* If we get here, the i particle (ii) has itself (jnr)
 696                      * in its neighborlist. This can only happen with the Verlet
 697                      * scheme, and corresponds to a self-interaction that will
 698                      * occur twice. Scale it down by 50% to only include it once.
 699                      */
 700                     VV *= half;
 701                 }
 702
 703                 for (i = 0; i < NSTATES; i++)
 704                 {
 705                     c6grid      = nbfp_grid[tj[i]];
 706                     vvtot      += LFV[i]*c6grid*VV;
 707                     Fscal      += LFV[i]*c6grid*FF;
 708                     dvdl_vdw   += (DLF[i]*c6grid)*VV;
 709                 }
 710             }
 711
 712             if (bDoForces)
 713             {
 714                 tx         = Fscal*dx;
 715                 ty         = Fscal*dy;
 716                 tz         = Fscal*dz;
 717                 fix        = fix + tx;
 718                 fiy        = fiy + ty;
 719                 fiz        = fiz + tz;
 720                 /* OpenMP atomics are expensive, but this kernels is also
 721                  * expensive, so we can take this hit, instead of using
 722                  * thread-local output buffers and extra reduction.
 723                  *
 724                  * All the OpenMP regions in this file are trivial and should
 725                  * not throw, so no need for try/catch.
 726                  */
 727 #pragma omp atomic
 728                 f[j3]     -= tx;
 729 #pragma omp atomic
 730                 f[j3+1]   -= ty;
 731 #pragma omp atomic
 732                 f[j3+2]   -= tz;
 733             }
 734         }
 735
 736         /* The atomics below are expensive with many OpenMP threads.
 737          * Here unperturbed i-particles will usually only have a few
 738          * (perturbed) j-particles in the list. Thus with a buffered list
 739          * we can skip a significant number of i-reductions with a check.
 740          */
 741         if (npair_within_cutoff > 0)
 742         {
 743             if (bDoForces)
 744             {
 745 #pragma omp atomic
 746                 f[ii3]        += fix;
 747 #pragma omp atomic
 748                 f[ii3+1]      += fiy;
 749 #pragma omp atomic
 750                 f[ii3+2]      += fiz;
 751             }
 752             if (bDoShiftForces)
 753             {
 754 #pragma omp atomic
 755                 fshift[is3]   += fix;
 756 #pragma omp atomic
 757                 fshift[is3+1] += fiy;
 758 #pragma omp atomic
 759                 fshift[is3+2] += fiz;
 760             }
 761             if (bDoPotential)
 762             {
 763                 ggid               = gid[n];
 764 #pragma omp atomic
 765                 Vc[ggid]          += vctot;
 766 #pragma omp atomic
 767                 Vv[ggid]          += vvtot;
 768             }
 769         }
 770     }
 771
 772 #pragma omp atomic
 773     dvdl[efptCOUL]     += dvdl_coul;
 774  #pragma omp atomic
 775     dvdl[efptVDW]      += dvdl_vdw;
 776
 777     /* Estimate flops, average for free energy stuff:
 778      * 12  flops per outer iteration
 779      * 150 flops per inner iteration
 780      */
 781 #pragma omp atomic
 782     inc_nrnb(nrnb, eNR_NBKERNEL_FREE_ENERGY, nlist->nri*12 + nlist->jindex[n]*150);
 783 }
 784
 785 void gmx_nb_free_energy_kernel(const t_nblist   *nlist,
 786                                rvec             *xx,
 787                                rvec             *ff,
 788                                t_forcerec       *fr,
 789                                const t_mdatoms  *mdatoms,
 790                                nb_kernel_data_t *kernel_data,
 791                                t_nrnb           *nrnb)
 792 {
 793     if (fr->sc_alphacoul == 0 && fr->sc_alphavdw == 0)
 794     {
 795         nb_free_energy_kernel<SoftCoreTreatment::None, false>(nlist, xx, ff, fr, mdatoms, kernel_data, nrnb);
 796     }
 797     else if (fr->sc_r_power == 6.0_real)
 798     {
 799         if (kernel_data->lambda[efptCOUL] == kernel_data->lambda[efptVDW] &&
 800             fr->sc_alphacoul == fr->sc_alphavdw)
 801         {
 802             nb_free_energy_kernel<SoftCoreTreatment::RPower6, false>(nlist, xx, ff, fr, mdatoms, kernel_data, nrnb);
 803         }
 804         else
 805         {
 806             nb_free_energy_kernel<SoftCoreTreatment::RPower6, true>(nlist, xx, ff, fr, mdatoms, kernel_data, nrnb);
 807         }
 808     }
 809     else if (fr->sc_r_power == 48.0_real)
 810     {
 811         nb_free_energy_kernel<SoftCoreTreatment::RPower48, true>(nlist, xx, ff, fr, mdatoms, kernel_data, nrnb);
 812     }
 813     else
 814     {
 815         GMX_RELEASE_ASSERT(false, "Unsupported soft-core r-power");
 816     }
 817 }