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