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/group.h"
  55 #include "gromacs/mdtypes/md_enums.h"
  56 #include "gromacs/mdtypes/nblist.h"
  57 #include "gromacs/mdtypes/simulation_workload.h"
  58 #include "gromacs/pbcutil/ishift.h"
  59 #include "gromacs/pbcutil/mshift.h"
  60 #include "gromacs/pbcutil/pbc.h"
  61 #include "gromacs/pbcutil/pbc_simd.h"
  62 #include "gromacs/simd/simd.h"
  63 #include "gromacs/simd/simd_math.h"
  64 #include "gromacs/simd/vector_operations.h"
  65 #include "gromacs/tables/forcetable.h"
  66 #include "gromacs/utility/basedefinitions.h"
  67 #include "gromacs/utility/fatalerror.h"
  68 #include "gromacs/utility/gmxassert.h"
  69
  70 #include "listed_internal.h"
  71
  72 using namespace gmx; // TODO: Remove when this file is moved into gmx namespace
  73
  74 /*! \brief Issue a warning if a listed interaction is beyond a table limit */
  75 static void warning_rlimit(const rvec* x, int ai, int aj, int* global_atom_index, real r, real rlimit)
  76 {
  77     gmx_warning(
  78             "Listed nonbonded interaction between particles %d and %d\n"
  79             "at distance %.3f which is larger than the table limit %.3f nm.\n\n"
  80             "This is likely either a 1,4 interaction, or a listed interaction inside\n"
  81             "a smaller molecule you are decoupling during a free energy calculation.\n"
  82             "Since interactions at distances beyond the table cannot be computed,\n"
  83             "they are skipped until they are inside the table limit again. You will\n"
  84             "only see this message once, even if it occurs for several interactions.\n\n"
  85             "IMPORTANT: This should not happen in a stable simulation, so there is\n"
  86             "probably something wrong with your system. Only change the table-extension\n"
  87             "distance in the mdp file if you are really sure that is the reason.\n",
  88             glatnr(global_atom_index, ai), glatnr(global_atom_index, aj), r, rlimit);
  89
  90     if (debug)
  91     {
  92         fprintf(debug,
  93                 "%8f %8f %8f\n%8f %8f %8f\n1-4 (%d,%d) interaction not within cut-off! r=%g. "
  94                 "Ignored\n",
  95                 x[ai][XX], x[ai][YY], x[ai][ZZ], x[aj][XX], x[aj][YY], x[aj][ZZ],
  96                 glatnr(global_atom_index, ai), glatnr(global_atom_index, aj), r);
  97     }
  98 }
  99
 100 /*! \brief Compute the energy and force for a single pair interaction */
 101 static real evaluate_single(real        r2,
 102                             real        tabscale,
 103                             const real* vftab,
 104                             real        tableStride,
 105                             real        qq,
 106                             real        c6,
 107                             real        c12,
 108                             real*       velec,
 109                             real*       vvdw)
 110 {
 111     real rinv, r, rtab, eps, eps2, Y, F, Geps, Heps2, Fp, VVe, FFe, VVd, FFd, VVr, FFr, fscal;
 112     int  ntab;
 113
 114     /* Do the tabulated interactions - first table lookup */
 115     rinv = gmx::invsqrt(r2);
 116     r    = r2 * rinv;
 117     rtab = r * tabscale;
 118     ntab = static_cast<int>(rtab);
 119     eps  = rtab - ntab;
 120     eps2 = eps * eps;
 121     ntab = static_cast<int>(tableStride * ntab);
 122     /* Electrostatics */
 123     Y     = vftab[ntab];
 124     F     = vftab[ntab + 1];
 125     Geps  = eps * vftab[ntab + 2];
 126     Heps2 = eps2 * vftab[ntab + 3];
 127     Fp    = F + Geps + Heps2;
 128     VVe   = Y + eps * Fp;
 129     FFe   = Fp + Geps + 2.0 * Heps2;
 130     /* Dispersion */
 131     Y     = vftab[ntab + 4];
 132     F     = vftab[ntab + 5];
 133     Geps  = eps * vftab[ntab + 6];
 134     Heps2 = eps2 * vftab[ntab + 7];
 135     Fp    = F + Geps + Heps2;
 136     VVd   = Y + eps * Fp;
 137     FFd   = Fp + Geps + 2.0 * Heps2;
 138     /* Repulsion */
 139     Y     = vftab[ntab + 8];
 140     F     = vftab[ntab + 9];
 141     Geps  = eps * vftab[ntab + 10];
 142     Heps2 = eps2 * vftab[ntab + 11];
 143     Fp    = F + Geps + Heps2;
 144     VVr   = Y + eps * Fp;
 145     FFr   = Fp + Geps + 2.0 * Heps2;
 146
 147     *velec = qq * VVe;
 148     *vvdw  = c6 * VVd + c12 * VVr;
 149
 150     fscal = -(qq * FFe + c6 * FFd + c12 * FFr) * tabscale * rinv;
 151
 152     return fscal;
 153 }
 154
 155 /*! \brief Compute the energy and force for a single pair interaction under FEP */
 156 static real free_energy_evaluate_single(real        r2,
 157                                         real        sc_r_power,
 158                                         real        alpha_coul,
 159                                         real        alpha_vdw,
 160                                         real        tabscale,
 161                                         const real* vftab,
 162                                         real        tableStride,
 163                                         real        qqA,
 164                                         real        c6A,
 165                                         real        c12A,
 166                                         real        qqB,
 167                                         real        c6B,
 168                                         real        c12B,
 169                                         const real  LFC[2],
 170                                         const real  LFV[2],
 171                                         const real  DLF[2],
 172                                         const real  lfac_coul[2],
 173                                         const real  lfac_vdw[2],
 174                                         const real  dlfac_coul[2],
 175                                         const real  dlfac_vdw[2],
 176                                         real        sigma6_def,
 177                                         real        sigma6_min,
 178                                         real        sigma2_def,
 179                                         real        sigma2_min,
 180                                         real*       velectot,
 181                                         real*       vvdwtot,
 182                                         real*       dvdl)
 183 {
 184     real       rp, rpm2, rtab, eps, eps2, Y, F, Geps, Heps2, Fp, VV, FF, fscal;
 185     real       qq[2], c6[2], c12[2], sigma6[2], sigma2[2], sigma_pow[2];
 186     real       alpha_coul_eff, alpha_vdw_eff, dvdl_coul, dvdl_vdw;
 187     real       rpinv, r_coul, r_vdw, velecsum, vvdwsum;
 188     real       fscal_vdw[2], fscal_elec[2];
 189     real       velec[2], vvdw[2];
 190     int        i, ntab;
 191     const real half     = 0.5;
 192     const real minusOne = -1.0;
 193     const real one      = 1.0;
 194     const real two      = 2.0;
 195     const real six      = 6.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
 210     {
 211         rp = std::pow(r2, half * sc_r_power); /* not currently supported as input, but can handle it */
 212         rpm2 = rp / r2;
 213     }
 214
 215     /* Loop over state A(0) and B(1) */
 216     for (i = 0; i < 2; i++)
 217     {
 218         if ((c6[i] > 0) && (c12[i] > 0))
 219         {
 220             /* The c6 & c12 coefficients now contain the constants 6.0 and 12.0, respectively.
 221              * Correct for this by multiplying with (1/12.0)/(1/6.0)=6.0/12.0=0.5.
 222              */
 223             sigma6[i] = half * c12[i] / c6[i];
 224             sigma2[i] = std::cbrt(half * c12[i] / c6[i]);
 225             /* should be able to get rid of this ^^^ internal pow call eventually.  Will require agreement on
 226                what data to store externally.  Can't be fixed without larger scale changes, so not 5.0 */
 227             if (sigma6[i] < sigma6_min) /* for disappearing coul and vdw with soft core at the same time */
 228             {
 229                 sigma6[i] = sigma6_min;
 230                 sigma2[i] = sigma2_min;
 231             }
 232         }
 233         else
 234         {
 235             sigma6[i] = sigma6_def;
 236             sigma2[i] = sigma2_def;
 237         }
 238         if (sc_r_power == six)
 239         {
 240             sigma_pow[i] = sigma6[i];
 241         }
 242         else
 243         { /* not really supported as input, but in here for testing the general case*/
 244             sigma_pow[i] = std::pow(sigma2[i], sc_r_power / 2);
 245         }
 246     }
 247
 248     /* only use softcore if one of the states has a zero endstate - softcore is for avoiding infinities!*/
 249     if ((c12[0] > 0) && (c12[1] > 0))
 250     {
 251         alpha_vdw_eff  = 0;
 252         alpha_coul_eff = 0;
 253     }
 254     else
 255     {
 256         alpha_vdw_eff  = alpha_vdw;
 257         alpha_coul_eff = alpha_coul;
 258     }
 259
 260     /* Loop over A and B states again */
 261     for (i = 0; i < 2; i++)
 262     {
 263         fscal_elec[i] = 0;
 264         fscal_vdw[i]  = 0;
 265         velec[i]      = 0;
 266         vvdw[i]       = 0;
 267
 268         /* Only spend time on A or B state if it is non-zero */
 269         if ((qq[i] != 0) || (c6[i] != 0) || (c12[i] != 0))
 270         {
 271             /* Coulomb */
 272             rpinv  = one / (alpha_coul_eff * lfac_coul[i] * sigma_pow[i] + rp);
 273             r_coul = std::pow(rpinv, minusOne / sc_r_power);
 274
 275             /* Electrostatics table lookup data */
 276             rtab = r_coul * tabscale;
 277             ntab = static_cast<int>(rtab);
 278             eps  = rtab - ntab;
 279             eps2 = eps * eps;
 280             ntab = static_cast<int>(tableStride * ntab);
 281             /* Electrostatics */
 282             Y             = vftab[ntab];
 283             F             = vftab[ntab + 1];
 284             Geps          = eps * vftab[ntab + 2];
 285             Heps2         = eps2 * vftab[ntab + 3];
 286             Fp            = F + Geps + Heps2;
 287             VV            = Y + eps * Fp;
 288             FF            = Fp + Geps + two * Heps2;
 289             velec[i]      = qq[i] * VV;
 290             fscal_elec[i] = -qq[i] * FF * r_coul * rpinv * tabscale;
 291
 292             /* Vdw */
 293             rpinv = one / (alpha_vdw_eff * lfac_vdw[i] * sigma_pow[i] + rp);
 294             r_vdw = std::pow(rpinv, minusOne / sc_r_power);
 295             /* Vdw table lookup data */
 296             rtab = r_vdw * tabscale;
 297             ntab = static_cast<int>(rtab);
 298             eps  = rtab - ntab;
 299             eps2 = eps * eps;
 300             ntab = 12 * ntab;
 301             /* Dispersion */
 302             Y            = vftab[ntab + 4];
 303             F            = vftab[ntab + 5];
 304             Geps         = eps * vftab[ntab + 6];
 305             Heps2        = eps2 * vftab[ntab + 7];
 306             Fp           = F + Geps + Heps2;
 307             VV           = Y + eps * Fp;
 308             FF           = Fp + Geps + two * Heps2;
 309             vvdw[i]      = c6[i] * VV;
 310             fscal_vdw[i] = -c6[i] * FF;
 311
 312             /* Repulsion */
 313             Y     = vftab[ntab + 8];
 314             F     = vftab[ntab + 9];
 315             Geps  = eps * vftab[ntab + 10];
 316             Heps2 = eps2 * vftab[ntab + 11];
 317             Fp    = F + Geps + Heps2;
 318             VV    = Y + eps * Fp;
 319             FF    = Fp + Geps + two * Heps2;
 320             vvdw[i] += c12[i] * VV;
 321             fscal_vdw[i] -= c12[i] * FF;
 322             fscal_vdw[i] *= r_vdw * rpinv * tabscale;
 323         }
 324     }
 325     /* Now we have velec[i], vvdw[i], and fscal[i] for both states */
 326     /* Assemble A and B states */
 327     velecsum  = 0;
 328     vvdwsum   = 0;
 329     dvdl_coul = 0;
 330     dvdl_vdw  = 0;
 331     fscal     = 0;
 332     for (i = 0; i < 2; i++)
 333     {
 334         velecsum += LFC[i] * velec[i];
 335         vvdwsum += LFV[i] * vvdw[i];
 336
 337         fscal += (LFC[i] * fscal_elec[i] + LFV[i] * fscal_vdw[i]) * rpm2;
 338
 339         dvdl_coul += velec[i] * DLF[i]
 340                      + LFC[i] * alpha_coul_eff * dlfac_coul[i] * fscal_elec[i] * sigma_pow[i];
 341         dvdl_vdw += vvdw[i] * DLF[i]
 342                     + LFV[i] * alpha_vdw_eff * dlfac_vdw[i] * fscal_vdw[i] * sigma_pow[i];
 343     }
 344
 345     dvdl[efptCOUL] += dvdl_coul;
 346     dvdl[efptVDW] += dvdl_vdw;
 347
 348     *velectot = velecsum;
 349     *vvdwtot  = vvdwsum;
 350
 351     return fscal;
 352 }
 353
 354 /*! \brief Calculate pair interactions, supports all types and conditions. */
 355 template<BondedKernelFlavor flavor>
 356 static real do_pairs_general(int                   ftype,
 357                              int                   nbonds,
 358                              const t_iatom         iatoms[],
 359                              const t_iparams       iparams[],
 360                              const rvec            x[],
 361                              rvec4                 f[],
 362                              rvec                  fshift[],
 363                              const struct t_pbc*   pbc,
 364                              const struct t_graph* g,
 365                              const real*           lambda,
 366                              real*                 dvdl,
 367                              const t_mdatoms*      md,
 368                              const t_forcerec*     fr,
 369                              gmx_grppairener_t*    grppener,
 370                              int*                  global_atom_index)
 371 {
 372     real            qq, c6, c12;
 373     rvec            dx;
 374     ivec            dt;
 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, fr->pairsTable->stride, qq, c6, c12, qqB, c6B, c12B, LFC,
 526                     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,
 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 (g)
 547             {
 548                 /* Correct the shift forces using the graph */
 549                 ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, aj), dt);
 550                 fshift_index = IVEC2IS(dt);
 551             }
 552             if (fshift_index != CENTRAL)
 553             {
 554                 rvec_inc(fshift[fshift_index], dx);
 555                 rvec_dec(fshift[CENTRAL], dx);
 556             }
 557         }
 558     }
 559     return 0.0;
 560 }
 561
 562 /*! \brief Calculate pairs, only for plain-LJ + plain Coulomb normal type.
 563  *
 564  * This function is templated for real/SimdReal and for optimization.
 565  */
 566 template<typename T, int pack_size, typename pbc_type>
 567 static void do_pairs_simple(int              nbonds,
 568                             const t_iatom    iatoms[],
 569                             const t_iparams  iparams[],
 570                             const rvec       x[],
 571                             rvec4            f[],
 572                             const pbc_type   pbc,
 573                             const t_mdatoms* md,
 574                             const real       scale_factor)
 575 {
 576     const int nfa1 = 1 + 2;
 577
 578     T six(6);
 579     T twelve(12);
 580     T ef(scale_factor);
 581
 582 #if GMX_SIMD_HAVE_REAL
 583     alignas(GMX_SIMD_ALIGNMENT) std::int32_t ai[pack_size];
 584     alignas(GMX_SIMD_ALIGNMENT) std::int32_t aj[pack_size];
 585     alignas(GMX_SIMD_ALIGNMENT) real         coeff[3 * pack_size];
 586 #else
 587     std::int32_t ai[pack_size];
 588     std::int32_t aj[pack_size];
 589     real         coeff[3 * pack_size];
 590 #endif
 591
 592     /* nbonds is #pairs*nfa1, here we step pack_size pairs */
 593     for (int i = 0; i < nbonds; i += pack_size * nfa1)
 594     {
 595         /* Collect atoms for pack_size pairs.
 596          * iu indexes into iatoms, we should not let iu go beyond nbonds.
 597          */
 598         int iu = i;
 599         for (int s = 0; s < pack_size; s++)
 600         {
 601             int itype = iatoms[iu];
 602             ai[s]     = iatoms[iu + 1];
 603             aj[s]     = iatoms[iu + 2];
 604
 605             if (i + s * nfa1 < nbonds)
 606             {
 607                 coeff[0 * pack_size + s] = iparams[itype].lj14.c6A;
 608                 coeff[1 * pack_size + s] = iparams[itype].lj14.c12A;
 609                 coeff[2 * pack_size + s] = md->chargeA[ai[s]] * md->chargeA[aj[s]];
 610
 611                 /* Avoid indexing the iatoms array out of bounds.
 612                  * We pad the coordinate indices with the last atom pair.
 613                  */
 614                 if (iu + nfa1 < nbonds)
 615                 {
 616                     iu += nfa1;
 617                 }
 618             }
 619             else
 620             {
 621                 /* Pad the coefficient arrays with zeros to get zero forces */
 622                 coeff[0 * pack_size + s] = 0;
 623                 coeff[1 * pack_size + s] = 0;
 624                 coeff[2 * pack_size + s] = 0;
 625             }
 626         }
 627
 628         /* Load the coordinates */
 629         T xi[DIM], xj[DIM];
 630         gatherLoadUTranspose<3>(reinterpret_cast<const real*>(x), ai, &xi[XX], &xi[YY], &xi[ZZ]);
 631         gatherLoadUTranspose<3>(reinterpret_cast<const real*>(x), aj, &xj[XX], &xj[YY], &xj[ZZ]);
 632
 633         T c6  = load<T>(coeff + 0 * pack_size);
 634         T c12 = load<T>(coeff + 1 * pack_size);
 635         T qq  = load<T>(coeff + 2 * pack_size);
 636
 637         /* We could save these operations by storing 6*C6,12*C12 */
 638         c6  = six * c6;
 639         c12 = twelve * c12;
 640
 641         T dr[DIM];
 642         pbc_dx_aiuc(pbc, xi, xj, dr);
 643
 644         T rsq   = dr[XX] * dr[XX] + dr[YY] * dr[YY] + dr[ZZ] * dr[ZZ];
 645         T rinv  = gmx::invsqrt(rsq);
 646         T rinv2 = rinv * rinv;
 647         T rinv6 = rinv2 * rinv2 * rinv2;
 648
 649         /* Calculate the Coulomb force * r */
 650         T cfr = ef * qq * rinv;
 651
 652         /* Calculate the LJ force * r and add it to the Coulomb part */
 653         T fr = gmx::fma(fms(c12, rinv6, c6), rinv6, cfr);
 654
 655         T finvr = fr * rinv2;
 656         T fx    = finvr * dr[XX];
 657         T fy    = finvr * dr[YY];
 658         T fz    = finvr * dr[ZZ];
 659
 660         /* Add the pair forces to the force array.
 661          * Note that here we might add multiple force components for some atoms
 662          * due to the SIMD padding. But the extra force components are zero.
 663          */
 664         transposeScatterIncrU<4>(reinterpret_cast<real*>(f), ai, fx, fy, fz);
 665         transposeScatterDecrU<4>(reinterpret_cast<real*>(f), aj, fx, fy, fz);
 666     }
 667 }
 668
 669 /*! \brief Calculate all listed pair interactions */
 670 void do_pairs(int                      ftype,
 671               int                      nbonds,
 672               const t_iatom            iatoms[],
 673               const t_iparams          iparams[],
 674               const rvec               x[],
 675               rvec4                    f[],
 676               rvec                     fshift[],
 677               const struct t_pbc*      pbc,
 678               const struct t_graph*    g,
 679               const real*              lambda,
 680               real*                    dvdl,
 681               const t_mdatoms*         md,
 682               const t_forcerec*        fr,
 683               const bool               havePerturbedInteractions,
 684               const gmx::StepWorkload& stepWork,
 685               gmx_grppairener_t*       grppener,
 686               int*                     global_atom_index)
 687 {
 688     if (ftype == F_LJ14 && fr->ic->vdwtype != evdwUSER && !EEL_USER(fr->ic->eeltype)
 689         && !havePerturbedInteractions && (!stepWork.computeVirial && !stepWork.computeEnergy))
 690     {
 691         /* We use a fast code-path for plain LJ 1-4 without FEP.
 692          *
 693          * TODO: Add support for energies (straightforward) and virial
 694          * in the SIMD template. For the virial it's inconvenient to store
 695          * the force sums for the shifts and we should directly calculate
 696          * and sum the virial for the shifts. But we should do this
 697          * at once for the angles and dihedrals as well.
 698          */
 699 #if GMX_SIMD_HAVE_REAL
 700         if (fr->use_simd_kernels)
 701         {
 702             alignas(GMX_SIMD_ALIGNMENT) real pbc_simd[9 * GMX_SIMD_REAL_WIDTH];
 703             set_pbc_simd(pbc, pbc_simd);
 704
 705             do_pairs_simple<SimdReal, GMX_SIMD_REAL_WIDTH, const real*>(
 706                     nbonds, iatoms, iparams, x, f, pbc_simd, md, fr->ic->epsfac * fr->fudgeQQ);
 707         }
 708         else
 709 #endif
 710         {
 711             /* This construct is needed because pbc_dx_aiuc doesn't accept pbc=NULL */
 712             t_pbc        pbc_no;
 713             const t_pbc* pbc_nonnull;
 714
 715             if (pbc != nullptr)
 716             {
 717                 pbc_nonnull = pbc;
 718             }
 719             else
 720             {
 721                 set_pbc(&pbc_no, PbcType::No, nullptr);
 722                 pbc_nonnull = &pbc_no;
 723             }
 724
 725             do_pairs_simple<real, 1, const t_pbc*>(nbonds, iatoms, iparams, x, f, pbc_nonnull, md,
 726                                                    fr->ic->epsfac * fr->fudgeQQ);
 727         }
 728     }
 729     else if (stepWork.computeVirial)
 730     {
 731         do_pairs_general<BondedKernelFlavor::ForcesAndVirialAndEnergy>(
 732                 ftype, nbonds, iatoms, iparams, x, f, fshift, pbc, g, lambda, dvdl, md, fr,
 733                 grppener, global_atom_index);
 734     }
 735     else
 736     {
 737         do_pairs_general<BondedKernelFlavor::ForcesAndEnergy>(ftype, nbonds, iatoms, iparams, x, f,
 738                                                               fshift, pbc, g, lambda, dvdl, md, fr,
 739                                                               grppener, global_atom_index);
 740     }
 741 }