src/gromacs/listed_forces/pairs.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2014,2015,2016,2017,2018,2019, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 /*! \internal \file
  36  *
  37  * \brief This file defines functions for "pair" interactions
  38  * (i.e. listed non-bonded interactions, e.g. 1-4 interactions)
  39  *
  40  * \author Mark Abraham <mark.j.abraham@gmail.com>
  41  *
  42  * \ingroup module_listed_forces
  43  */
  44 #include "gmxpre.h"
  45
  46 #include "pairs.h"
  47
  48 #include <cmath>
  49
  50 #include "gromacs/listed_forces/bonded.h"
  51 #include "gromacs/math/functions.h"
  52 #include "gromacs/math/vec.h"
  53 #include "gromacs/mdtypes/group.h"
  54 #include "gromacs/mdtypes/md_enums.h"
  55 #include "gromacs/mdtypes/nblist.h"
  56 #include "gromacs/mdtypes/simulation_workload.h"
  57 #include "gromacs/pbcutil/ishift.h"
  58 #include "gromacs/pbcutil/mshift.h"
  59 #include "gromacs/pbcutil/pbc.h"
  60 #include "gromacs/pbcutil/pbc_simd.h"
  61 #include "gromacs/simd/simd.h"
  62 #include "gromacs/simd/simd_math.h"
  63 #include "gromacs/simd/vector_operations.h"
  64 #include "gromacs/tables/forcetable.h"
  65 #include "gromacs/utility/basedefinitions.h"
  66 #include "gromacs/utility/fatalerror.h"
  67 #include "gromacs/utility/gmxassert.h"
  68
  69 #include "listed_internal.h"
  70
  71 using namespace gmx; // TODO: Remove when this file is moved into gmx namespace
  72
  73 /*! \brief Issue a warning if a listed interaction is beyond a table limit */
  74 static void warning_rlimit(const rvec* x, int ai, int aj, int* global_atom_index, real r, real rlimit)
  75 {
  76     gmx_warning(
  77             "Listed nonbonded interaction between particles %d and %d\n"
  78             "at distance %.3f which is larger than the table limit %.3f nm.\n\n"
  79             "This is likely either a 1,4 interaction, or a listed interaction inside\n"
  80             "a smaller molecule you are decoupling during a free energy calculation.\n"
  81             "Since interactions at distances beyond the table cannot be computed,\n"
  82             "they are skipped until they are inside the table limit again. You will\n"
  83             "only see this message once, even if it occurs for several interactions.\n\n"
  84             "IMPORTANT: This should not happen in a stable simulation, so there is\n"
  85             "probably something wrong with your system. Only change the table-extension\n"
  86             "distance in the mdp file if you are really sure that is the reason.\n",
  87             glatnr(global_atom_index, ai), glatnr(global_atom_index, aj), r, rlimit);
  88
  89     if (debug)
  90     {
  91         fprintf(debug,
  92                 "%8f %8f %8f\n%8f %8f %8f\n1-4 (%d,%d) interaction not within cut-off! r=%g. "
  93                 "Ignored\n",
  94                 x[ai][XX], x[ai][YY], x[ai][ZZ], x[aj][XX], x[aj][YY], x[aj][ZZ],
  95                 glatnr(global_atom_index, ai), glatnr(global_atom_index, aj), r);
  96     }
  97 }
  98
  99 /*! \brief Compute the energy and force for a single pair interaction */
 100 static real evaluate_single(real        r2,
 101                             real        tabscale,
 102                             const real* vftab,
 103                             real        tableStride,
 104                             real        qq,
 105                             real        c6,
 106                             real        c12,
 107                             real*       velec,
 108                             real*       vvdw)
 109 {
 110     real rinv, r, rtab, eps, eps2, Y, F, Geps, Heps2, Fp, VVe, FFe, VVd, FFd, VVr, FFr, fscal;
 111     int  ntab;
 112
 113     /* Do the tabulated interactions - first table lookup */
 114     rinv = gmx::invsqrt(r2);
 115     r    = r2 * rinv;
 116     rtab = r * tabscale;
 117     ntab = static_cast<int>(rtab);
 118     eps  = rtab - ntab;
 119     eps2 = eps * eps;
 120     ntab = static_cast<int>(tableStride * ntab);
 121     /* Electrostatics */
 122     Y     = vftab[ntab];
 123     F     = vftab[ntab + 1];
 124     Geps  = eps * vftab[ntab + 2];
 125     Heps2 = eps2 * vftab[ntab + 3];
 126     Fp    = F + Geps + Heps2;
 127     VVe   = Y + eps * Fp;
 128     FFe   = Fp + Geps + 2.0 * Heps2;
 129     /* Dispersion */
 130     Y     = vftab[ntab + 4];
 131     F     = vftab[ntab + 5];
 132     Geps  = eps * vftab[ntab + 6];
 133     Heps2 = eps2 * vftab[ntab + 7];
 134     Fp    = F + Geps + Heps2;
 135     VVd   = Y + eps * Fp;
 136     FFd   = Fp + Geps + 2.0 * Heps2;
 137     /* Repulsion */
 138     Y     = vftab[ntab + 8];
 139     F     = vftab[ntab + 9];
 140     Geps  = eps * vftab[ntab + 10];
 141     Heps2 = eps2 * vftab[ntab + 11];
 142     Fp    = F + Geps + Heps2;
 143     VVr   = Y + eps * Fp;
 144     FFr   = Fp + Geps + 2.0 * Heps2;
 145
 146     *velec = qq * VVe;
 147     *vvdw  = c6 * VVd + c12 * VVr;
 148
 149     fscal = -(qq * FFe + c6 * FFd + c12 * FFr) * tabscale * rinv;
 150
 151     return fscal;
 152 }
 153
 154 /*! \brief Compute the energy and force for a single pair interaction under FEP */
 155 static real free_energy_evaluate_single(real        r2,
 156                                         real        sc_r_power,
 157                                         real        alpha_coul,
 158                                         real        alpha_vdw,
 159                                         real        tabscale,
 160                                         const real* vftab,
 161                                         real        tableStride,
 162                                         real        qqA,
 163                                         real        c6A,
 164                                         real        c12A,
 165                                         real        qqB,
 166                                         real        c6B,
 167                                         real        c12B,
 168                                         const real  LFC[2],
 169                                         const real  LFV[2],
 170                                         const real  DLF[2],
 171                                         const real  lfac_coul[2],
 172                                         const real  lfac_vdw[2],
 173                                         const real  dlfac_coul[2],
 174                                         const real  dlfac_vdw[2],
 175                                         real        sigma6_def,
 176                                         real        sigma6_min,
 177                                         real        sigma2_def,
 178                                         real        sigma2_min,
 179                                         real*       velectot,
 180                                         real*       vvdwtot,
 181                                         real*       dvdl)
 182 {
 183     real       rp, rpm2, rtab, eps, eps2, Y, F, Geps, Heps2, Fp, VV, FF, fscal;
 184     real       qq[2], c6[2], c12[2], sigma6[2], sigma2[2], sigma_pow[2];
 185     real       alpha_coul_eff, alpha_vdw_eff, dvdl_coul, dvdl_vdw;
 186     real       rpinv, r_coul, r_vdw, velecsum, vvdwsum;
 187     real       fscal_vdw[2], fscal_elec[2];
 188     real       velec[2], vvdw[2];
 189     int        i, ntab;
 190     const real half        = 0.5;
 191     const real minusOne    = -1.0;
 192     const real one         = 1.0;
 193     const real two         = 2.0;
 194     const real six         = 6.0;
 195     const real fourtyeight = 48.0;
 196
 197     qq[0]  = qqA;
 198     qq[1]  = qqB;
 199     c6[0]  = c6A;
 200     c6[1]  = c6B;
 201     c12[0] = c12A;
 202     c12[1] = c12B;
 203
 204     if (sc_r_power == six)
 205     {
 206         rpm2 = r2 * r2;   /* r4 */
 207         rp   = rpm2 * r2; /* r6 */
 208     }
 209     else if (sc_r_power == fourtyeight)
 210     {
 211         rp   = r2 * r2 * r2; /* r6 */
 212         rp   = rp * rp;      /* r12 */
 213         rp   = rp * rp;      /* r24 */
 214         rp   = rp * rp;      /* r48 */
 215         rpm2 = rp / r2;      /* r46 */
 216     }
 217     else
 218     {
 219         rp = std::pow(r2, half * sc_r_power); /* not currently supported as input, but can handle it */
 220         rpm2 = rp / r2;
 221     }
 222
 223     /* Loop over state A(0) and B(1) */
 224     for (i = 0; i < 2; i++)
 225     {
 226         if ((c6[i] > 0) && (c12[i] > 0))
 227         {
 228             /* The c6 & c12 coefficients now contain the constants 6.0 and 12.0, respectively.
 229              * Correct for this by multiplying with (1/12.0)/(1/6.0)=6.0/12.0=0.5.
 230              */
 231             sigma6[i] = half * c12[i] / c6[i];
 232             sigma2[i] = std::cbrt(half * c12[i] / c6[i]);
 233             /* should be able to get rid of this ^^^ internal pow call eventually.  Will require agreement on
 234                what data to store externally.  Can't be fixed without larger scale changes, so not 5.0 */
 235             if (sigma6[i] < sigma6_min) /* for disappearing coul and vdw with soft core at the same time */
 236             {
 237                 sigma6[i] = sigma6_min;
 238                 sigma2[i] = sigma2_min;
 239             }
 240         }
 241         else
 242         {
 243             sigma6[i] = sigma6_def;
 244             sigma2[i] = sigma2_def;
 245         }
 246         if (sc_r_power == six)
 247         {
 248             sigma_pow[i] = sigma6[i];
 249         }
 250         else if (sc_r_power == fourtyeight)
 251         {
 252             sigma_pow[i] = sigma6[i] * sigma6[i];       /* sigma^12 */
 253             sigma_pow[i] = sigma_pow[i] * sigma_pow[i]; /* sigma^24 */
 254             sigma_pow[i] = sigma_pow[i] * sigma_pow[i]; /* sigma^48 */
 255         }
 256         else
 257         { /* not really supported as input, but in here for testing the general case*/
 258             sigma_pow[i] = std::pow(sigma2[i], sc_r_power / 2);
 259         }
 260     }
 261
 262     /* only use softcore if one of the states has a zero endstate - softcore is for avoiding infinities!*/
 263     if ((c12[0] > 0) && (c12[1] > 0))
 264     {
 265         alpha_vdw_eff  = 0;
 266         alpha_coul_eff = 0;
 267     }
 268     else
 269     {
 270         alpha_vdw_eff  = alpha_vdw;
 271         alpha_coul_eff = alpha_coul;
 272     }
 273
 274     /* Loop over A and B states again */
 275     for (i = 0; i < 2; i++)
 276     {
 277         fscal_elec[i] = 0;
 278         fscal_vdw[i]  = 0;
 279         velec[i]      = 0;
 280         vvdw[i]       = 0;
 281
 282         /* Only spend time on A or B state if it is non-zero */
 283         if ((qq[i] != 0) || (c6[i] != 0) || (c12[i] != 0))
 284         {
 285             /* Coulomb */
 286             rpinv  = one / (alpha_coul_eff * lfac_coul[i] * sigma_pow[i] + rp);
 287             r_coul = std::pow(rpinv, minusOne / sc_r_power);
 288
 289             /* Electrostatics table lookup data */
 290             rtab = r_coul * tabscale;
 291             ntab = static_cast<int>(rtab);
 292             eps  = rtab - ntab;
 293             eps2 = eps * eps;
 294             ntab = static_cast<int>(tableStride * ntab);
 295             /* Electrostatics */
 296             Y             = vftab[ntab];
 297             F             = vftab[ntab + 1];
 298             Geps          = eps * vftab[ntab + 2];
 299             Heps2         = eps2 * vftab[ntab + 3];
 300             Fp            = F + Geps + Heps2;
 301             VV            = Y + eps * Fp;
 302             FF            = Fp + Geps + two * Heps2;
 303             velec[i]      = qq[i] * VV;
 304             fscal_elec[i] = -qq[i] * FF * r_coul * rpinv * tabscale;
 305
 306             /* Vdw */
 307             rpinv = one / (alpha_vdw_eff * lfac_vdw[i] * sigma_pow[i] + rp);
 308             r_vdw = std::pow(rpinv, minusOne / sc_r_power);
 309             /* Vdw table lookup data */
 310             rtab = r_vdw * tabscale;
 311             ntab = static_cast<int>(rtab);
 312             eps  = rtab - ntab;
 313             eps2 = eps * eps;
 314             ntab = 12 * ntab;
 315             /* Dispersion */
 316             Y            = vftab[ntab + 4];
 317             F            = vftab[ntab + 5];
 318             Geps         = eps * vftab[ntab + 6];
 319             Heps2        = eps2 * vftab[ntab + 7];
 320             Fp           = F + Geps + Heps2;
 321             VV           = Y + eps * Fp;
 322             FF           = Fp + Geps + two * Heps2;
 323             vvdw[i]      = c6[i] * VV;
 324             fscal_vdw[i] = -c6[i] * FF;
 325
 326             /* Repulsion */
 327             Y     = vftab[ntab + 8];
 328             F     = vftab[ntab + 9];
 329             Geps  = eps * vftab[ntab + 10];
 330             Heps2 = eps2 * vftab[ntab + 11];
 331             Fp    = F + Geps + Heps2;
 332             VV    = Y + eps * Fp;
 333             FF    = Fp + Geps + two * Heps2;
 334             vvdw[i] += c12[i] * VV;
 335             fscal_vdw[i] -= c12[i] * FF;
 336             fscal_vdw[i] *= r_vdw * rpinv * tabscale;
 337         }
 338     }
 339     /* Now we have velec[i], vvdw[i], and fscal[i] for both states */
 340     /* Assemble A and B states */
 341     velecsum  = 0;
 342     vvdwsum   = 0;
 343     dvdl_coul = 0;
 344     dvdl_vdw  = 0;
 345     fscal     = 0;
 346     for (i = 0; i < 2; i++)
 347     {
 348         velecsum += LFC[i] * velec[i];
 349         vvdwsum += LFV[i] * vvdw[i];
 350
 351         fscal += (LFC[i] * fscal_elec[i] + LFV[i] * fscal_vdw[i]) * rpm2;
 352
 353         dvdl_coul += velec[i] * DLF[i]
 354                      + LFC[i] * alpha_coul_eff * dlfac_coul[i] * fscal_elec[i] * sigma_pow[i];
 355         dvdl_vdw += vvdw[i] * DLF[i]
 356                     + LFV[i] * alpha_vdw_eff * dlfac_vdw[i] * fscal_vdw[i] * sigma_pow[i];
 357     }
 358
 359     dvdl[efptCOUL] += dvdl_coul;
 360     dvdl[efptVDW] += dvdl_vdw;
 361
 362     *velectot = velecsum;
 363     *vvdwtot  = vvdwsum;
 364
 365     return fscal;
 366 }
 367
 368 /*! \brief Calculate pair interactions, supports all types and conditions. */
 369 template<BondedKernelFlavor flavor>
 370 static real do_pairs_general(int                   ftype,
 371                              int                   nbonds,
 372                              const t_iatom         iatoms[],
 373                              const t_iparams       iparams[],
 374                              const rvec            x[],
 375                              rvec4                 f[],
 376                              rvec                  fshift[],
 377                              const struct t_pbc*   pbc,
 378                              const struct t_graph* g,
 379                              const real*           lambda,
 380                              real*                 dvdl,
 381                              const t_mdatoms*      md,
 382                              const t_forcerec*     fr,
 383                              gmx_grppairener_t*    grppener,
 384                              int*                  global_atom_index)
 385 {
 386     real            qq, c6, c12;
 387     rvec            dx;
 388     ivec            dt;
 389     int             i, itype, ai, aj, gid;
 390     int             fshift_index;
 391     real            r2;
 392     real            fscal, velec, vvdw;
 393     real*           energygrp_elec;
 394     real*           energygrp_vdw;
 395     static gmx_bool warned_rlimit = FALSE;
 396     /* Free energy stuff */
 397     gmx_bool bFreeEnergy;
 398     real     LFC[2], LFV[2], DLF[2], lfac_coul[2], lfac_vdw[2], dlfac_coul[2], dlfac_vdw[2];
 399     real     qqB, c6B, c12B, sigma2_def, sigma2_min;
 400
 401     switch (ftype)
 402     {
 403         case F_LJ14:
 404         case F_LJC14_Q:
 405             energygrp_elec = grppener->ener[egCOUL14].data();
 406             energygrp_vdw  = grppener->ener[egLJ14].data();
 407             break;
 408         case F_LJC_PAIRS_NB:
 409             energygrp_elec = grppener->ener[egCOULSR].data();
 410             energygrp_vdw  = grppener->ener[egLJSR].data();
 411             break;
 412         default:
 413             energygrp_elec = nullptr; /* Keep compiler happy */
 414             energygrp_vdw  = nullptr; /* Keep compiler happy */
 415             gmx_fatal(FARGS, "Unknown function type %d in do_nonbonded14", ftype);
 416     }
 417
 418     if (fr->efep != efepNO)
 419     {
 420         /* Lambda factor for state A=1-lambda and B=lambda */
 421         LFC[0] = 1.0 - lambda[efptCOUL];
 422         LFV[0] = 1.0 - lambda[efptVDW];
 423         LFC[1] = lambda[efptCOUL];
 424         LFV[1] = lambda[efptVDW];
 425
 426         /*derivative of the lambda factor for state A and B */
 427         DLF[0] = -1;
 428         DLF[1] = 1;
 429
 430         /* precalculate */
 431         sigma2_def = std::cbrt(fr->sc_sigma6_def);
 432         sigma2_min = std::cbrt(fr->sc_sigma6_min);
 433
 434         for (i = 0; i < 2; i++)
 435         {
 436             lfac_coul[i] = (fr->sc_power == 2 ? (1 - LFC[i]) * (1 - LFC[i]) : (1 - LFC[i]));
 437             dlfac_coul[i] =
 438                     DLF[i] * fr->sc_power / fr->sc_r_power * (fr->sc_power == 2 ? (1 - LFC[i]) : 1);
 439             lfac_vdw[i] = (fr->sc_power == 2 ? (1 - LFV[i]) * (1 - LFV[i]) : (1 - LFV[i]));
 440             dlfac_vdw[i] =
 441                     DLF[i] * fr->sc_power / fr->sc_r_power * (fr->sc_power == 2 ? (1 - LFV[i]) : 1);
 442         }
 443     }
 444     else
 445     {
 446         sigma2_min = sigma2_def = 0;
 447     }
 448
 449     /* TODO This code depends on the logic in tables.c that constructs
 450        the table layout, which should be made explicit in future
 451        cleanup. */
 452     GMX_ASSERT(
 453             etiNR == 3,
 454             "Pair-interaction code that uses GROMACS interaction tables supports exactly 3 tables");
 455     GMX_ASSERT(
 456             fr->pairsTable->interaction == GMX_TABLE_INTERACTION_ELEC_VDWREP_VDWDISP,
 457             "Pair interaction kernels need a table with Coulomb, repulsion and dispersion entries");
 458
 459     const real epsfac = fr->ic->epsfac;
 460
 461     bFreeEnergy = FALSE;
 462     for (i = 0; (i < nbonds);)
 463     {
 464         itype = iatoms[i++];
 465         ai    = iatoms[i++];
 466         aj    = iatoms[i++];
 467         gid   = GID(md->cENER[ai], md->cENER[aj], md->nenergrp);
 468
 469         /* Get parameters */
 470         switch (ftype)
 471         {
 472             case F_LJ14:
 473                 bFreeEnergy = (fr->efep != efepNO
 474                                && (((md->nPerturbed != 0) && (md->bPerturbed[ai] || md->bPerturbed[aj]))
 475                                    || iparams[itype].lj14.c6A != iparams[itype].lj14.c6B
 476                                    || iparams[itype].lj14.c12A != iparams[itype].lj14.c12B));
 477                 qq          = md->chargeA[ai] * md->chargeA[aj] * epsfac * fr->fudgeQQ;
 478                 c6          = iparams[itype].lj14.c6A;
 479                 c12         = iparams[itype].lj14.c12A;
 480                 break;
 481             case F_LJC14_Q:
 482                 qq = iparams[itype].ljc14.qi * iparams[itype].ljc14.qj * epsfac
 483                      * iparams[itype].ljc14.fqq;
 484                 c6  = iparams[itype].ljc14.c6;
 485                 c12 = iparams[itype].ljc14.c12;
 486                 break;
 487             case F_LJC_PAIRS_NB:
 488                 qq  = iparams[itype].ljcnb.qi * iparams[itype].ljcnb.qj * epsfac;
 489                 c6  = iparams[itype].ljcnb.c6;
 490                 c12 = iparams[itype].ljcnb.c12;
 491                 break;
 492             default:
 493                 /* Cannot happen since we called gmx_fatal() above in this case */
 494                 qq = c6 = c12 = 0; /* Keep compiler happy */
 495                 break;
 496         }
 497
 498         /* To save flops in the optimized kernels, c6/c12 have 6.0/12.0 derivative prefactors
 499          * included in the general nfbp array now. This means the tables are scaled down by the
 500          * same factor, so when we use the original c6/c12 parameters from iparams[] they must
 501          * be scaled up.
 502          */
 503         c6 *= 6.0;
 504         c12 *= 12.0;
 505
 506         /* Do we need to apply full periodic boundary conditions? */
 507         if (fr->bMolPBC)
 508         {
 509             fshift_index = pbc_dx_aiuc(pbc, x[ai], x[aj], dx);
 510         }
 511         else
 512         {
 513             fshift_index = CENTRAL;
 514             rvec_sub(x[ai], x[aj], dx);
 515         }
 516         r2 = norm2(dx);
 517
 518         if (r2 >= fr->pairsTable->r * fr->pairsTable->r)
 519         {
 520             /* This check isn't race free. But it doesn't matter because if a race occurs the only
 521              * disadvantage is that the warning is printed twice */
 522             if (!warned_rlimit)
 523             {
 524                 warning_rlimit(x, ai, aj, global_atom_index, sqrt(r2), fr->pairsTable->r);
 525                 warned_rlimit = TRUE;
 526             }
 527             continue;
 528         }
 529
 530         if (bFreeEnergy)
 531         {
 532             /* Currently free energy is only supported for F_LJ14, so no need to check for that if we got here */
 533             qqB  = md->chargeB[ai] * md->chargeB[aj] * epsfac * fr->fudgeQQ;
 534             c6B  = iparams[itype].lj14.c6B * 6.0;
 535             c12B = iparams[itype].lj14.c12B * 12.0;
 536
 537             fscal = free_energy_evaluate_single(
 538                     r2, fr->sc_r_power, fr->sc_alphacoul, fr->sc_alphavdw, fr->pairsTable->scale,
 539                     fr->pairsTable->data, fr->pairsTable->stride, qq, c6, c12, qqB, c6B, c12B, LFC,
 540                     LFV, DLF, lfac_coul, lfac_vdw, dlfac_coul, dlfac_vdw, fr->sc_sigma6_def,
 541                     fr->sc_sigma6_min, sigma2_def, sigma2_min, &velec, &vvdw, dvdl);
 542         }
 543         else
 544         {
 545             /* Evaluate tabulated interaction without free energy */
 546             fscal = evaluate_single(r2, fr->pairsTable->scale, fr->pairsTable->data,
 547                                     fr->pairsTable->stride, qq, c6, c12, &velec, &vvdw);
 548         }
 549
 550         energygrp_elec[gid] += velec;
 551         energygrp_vdw[gid] += vvdw;
 552         svmul(fscal, dx, dx);
 553
 554         /* Add the forces */
 555         rvec_inc(f[ai], dx);
 556         rvec_dec(f[aj], dx);
 557
 558         if (computeVirial(flavor))
 559         {
 560             if (g)
 561             {
 562                 /* Correct the shift forces using the graph */
 563                 ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, aj), dt);
 564                 fshift_index = IVEC2IS(dt);
 565             }
 566             if (fshift_index != CENTRAL)
 567             {
 568                 rvec_inc(fshift[fshift_index], dx);
 569                 rvec_dec(fshift[CENTRAL], dx);
 570             }
 571         }
 572     }
 573     return 0.0;
 574 }
 575
 576 /*! \brief Calculate pairs, only for plain-LJ + plain Coulomb normal type.
 577  *
 578  * This function is templated for real/SimdReal and for optimization.
 579  */
 580 template<typename T, int pack_size, typename pbc_type>
 581 static void do_pairs_simple(int              nbonds,
 582                             const t_iatom    iatoms[],
 583                             const t_iparams  iparams[],
 584                             const rvec       x[],
 585                             rvec4            f[],
 586                             const pbc_type   pbc,
 587                             const t_mdatoms* md,
 588                             const real       scale_factor)
 589 {
 590     const int nfa1 = 1 + 2;
 591
 592     T six(6);
 593     T twelve(12);
 594     T ef(scale_factor);
 595
 596 #if GMX_SIMD_HAVE_REAL
 597     alignas(GMX_SIMD_ALIGNMENT) std::int32_t ai[pack_size];
 598     alignas(GMX_SIMD_ALIGNMENT) std::int32_t aj[pack_size];
 599     alignas(GMX_SIMD_ALIGNMENT) real         coeff[3 * pack_size];
 600 #else
 601     std::int32_t ai[pack_size];
 602     std::int32_t aj[pack_size];
 603     real         coeff[3 * pack_size];
 604 #endif
 605
 606     /* nbonds is #pairs*nfa1, here we step pack_size pairs */
 607     for (int i = 0; i < nbonds; i += pack_size * nfa1)
 608     {
 609         /* Collect atoms for pack_size pairs.
 610          * iu indexes into iatoms, we should not let iu go beyond nbonds.
 611          */
 612         int iu = i;
 613         for (int s = 0; s < pack_size; s++)
 614         {
 615             int itype = iatoms[iu];
 616             ai[s]     = iatoms[iu + 1];
 617             aj[s]     = iatoms[iu + 2];
 618
 619             if (i + s * nfa1 < nbonds)
 620             {
 621                 coeff[0 * pack_size + s] = iparams[itype].lj14.c6A;
 622                 coeff[1 * pack_size + s] = iparams[itype].lj14.c12A;
 623                 coeff[2 * pack_size + s] = md->chargeA[ai[s]] * md->chargeA[aj[s]];
 624
 625                 /* Avoid indexing the iatoms array out of bounds.
 626                  * We pad the coordinate indices with the last atom pair.
 627                  */
 628                 if (iu + nfa1 < nbonds)
 629                 {
 630                     iu += nfa1;
 631                 }
 632             }
 633             else
 634             {
 635                 /* Pad the coefficient arrays with zeros to get zero forces */
 636                 coeff[0 * pack_size + s] = 0;
 637                 coeff[1 * pack_size + s] = 0;
 638                 coeff[2 * pack_size + s] = 0;
 639             }
 640         }
 641
 642         /* Load the coordinates */
 643         T xi[DIM], xj[DIM];
 644         gatherLoadUTranspose<3>(reinterpret_cast<const real*>(x), ai, &xi[XX], &xi[YY], &xi[ZZ]);
 645         gatherLoadUTranspose<3>(reinterpret_cast<const real*>(x), aj, &xj[XX], &xj[YY], &xj[ZZ]);
 646
 647         T c6  = load<T>(coeff + 0 * pack_size);
 648         T c12 = load<T>(coeff + 1 * pack_size);
 649         T qq  = load<T>(coeff + 2 * pack_size);
 650
 651         /* We could save these operations by storing 6*C6,12*C12 */
 652         c6  = six * c6;
 653         c12 = twelve * c12;
 654
 655         T dr[DIM];
 656         pbc_dx_aiuc(pbc, xi, xj, dr);
 657
 658         T rsq   = dr[XX] * dr[XX] + dr[YY] * dr[YY] + dr[ZZ] * dr[ZZ];
 659         T rinv  = gmx::invsqrt(rsq);
 660         T rinv2 = rinv * rinv;
 661         T rinv6 = rinv2 * rinv2 * rinv2;
 662
 663         /* Calculate the Coulomb force * r */
 664         T cfr = ef * qq * rinv;
 665
 666         /* Calculate the LJ force * r and add it to the Coulomb part */
 667         T fr = gmx::fma(fms(c12, rinv6, c6), rinv6, cfr);
 668
 669         T finvr = fr * rinv2;
 670         T fx    = finvr * dr[XX];
 671         T fy    = finvr * dr[YY];
 672         T fz    = finvr * dr[ZZ];
 673
 674         /* Add the pair forces to the force array.
 675          * Note that here we might add multiple force components for some atoms
 676          * due to the SIMD padding. But the extra force components are zero.
 677          */
 678         transposeScatterIncrU<4>(reinterpret_cast<real*>(f), ai, fx, fy, fz);
 679         transposeScatterDecrU<4>(reinterpret_cast<real*>(f), aj, fx, fy, fz);
 680     }
 681 }
 682
 683 /*! \brief Calculate all listed pair interactions */
 684 void do_pairs(int                      ftype,
 685               int                      nbonds,
 686               const t_iatom            iatoms[],
 687               const t_iparams          iparams[],
 688               const rvec               x[],
 689               rvec4                    f[],
 690               rvec                     fshift[],
 691               const struct t_pbc*      pbc,
 692               const struct t_graph*    g,
 693               const real*              lambda,
 694               real*                    dvdl,
 695               const t_mdatoms*         md,
 696               const t_forcerec*        fr,
 697               const bool               havePerturbedInteractions,
 698               const gmx::StepWorkload& stepWork,
 699               gmx_grppairener_t*       grppener,
 700               int*                     global_atom_index)
 701 {
 702     if (ftype == F_LJ14 && fr->ic->vdwtype != evdwUSER && !EEL_USER(fr->ic->eeltype)
 703         && !havePerturbedInteractions && (!stepWork.computeVirial && !stepWork.computeEnergy))
 704     {
 705         /* We use a fast code-path for plain LJ 1-4 without FEP.
 706          *
 707          * TODO: Add support for energies (straightforward) and virial
 708          * in the SIMD template. For the virial it's inconvenient to store
 709          * the force sums for the shifts and we should directly calculate
 710          * and sum the virial for the shifts. But we should do this
 711          * at once for the angles and dihedrals as well.
 712          */
 713 #if GMX_SIMD_HAVE_REAL
 714         if (fr->use_simd_kernels)
 715         {
 716             alignas(GMX_SIMD_ALIGNMENT) real pbc_simd[9 * GMX_SIMD_REAL_WIDTH];
 717             set_pbc_simd(pbc, pbc_simd);
 718
 719             do_pairs_simple<SimdReal, GMX_SIMD_REAL_WIDTH, const real*>(
 720                     nbonds, iatoms, iparams, x, f, pbc_simd, md, fr->ic->epsfac * fr->fudgeQQ);
 721         }
 722         else
 723 #endif
 724         {
 725             /* This construct is needed because pbc_dx_aiuc doesn't accept pbc=NULL */
 726             t_pbc        pbc_no;
 727             const t_pbc* pbc_nonnull;
 728
 729             if (pbc != nullptr)
 730             {
 731                 pbc_nonnull = pbc;
 732             }
 733             else
 734             {
 735                 set_pbc(&pbc_no, epbcNONE, nullptr);
 736                 pbc_nonnull = &pbc_no;
 737             }
 738
 739             do_pairs_simple<real, 1, const t_pbc*>(nbonds, iatoms, iparams, x, f, pbc_nonnull, md,
 740                                                    fr->ic->epsfac * fr->fudgeQQ);
 741         }
 742     }
 743     else if (stepWork.computeVirial)
 744     {
 745         do_pairs_general<BondedKernelFlavor::ForcesAndVirialAndEnergy>(
 746                 ftype, nbonds, iatoms, iparams, x, f, fshift, pbc, g, lambda, dvdl, md, fr,
 747                 grppener, global_atom_index);
 748     }
 749     else
 750     {
 751         do_pairs_general<BondedKernelFlavor::ForcesAndEnergy>(ftype, nbonds, iatoms, iparams, x, f,
 752                                                               fshift, pbc, g, lambda, dvdl, md, fr,
 753                                                               grppener, global_atom_index);
 754     }
 755 }