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