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