src/gromacs/listed_forces/bonded.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 /*! \internal \file
  39  *
  40  * \brief This file defines low-level functions necessary for
  41  * computing energies and forces for listed interactions.
  42  *
  43  * \author Mark Abraham <mark.j.abraham@gmail.com>
  44  *
  45  * \ingroup module_listed_forces
  46  */
  47 #include "gmxpre.h"
  48
  49 #include "bonded.h"
  50
  51 #include "config.h"
  52
  53 #include <cassert>
  54 #include <cmath>
  55
  56 #include <algorithm>
  57
  58 #include "gromacs/gmxlib/nrnb.h"
  59 #include "gromacs/listed_forces/disre.h"
  60 #include "gromacs/listed_forces/orires.h"
  61 #include "gromacs/math/functions.h"
  62 #include "gromacs/math/units.h"
  63 #include "gromacs/math/utilities.h"
  64 #include "gromacs/math/vec.h"
  65 #include "gromacs/mdtypes/fcdata.h"
  66 #include "gromacs/mdtypes/mdatom.h"
  67 #include "gromacs/pbcutil/ishift.h"
  68 #include "gromacs/pbcutil/pbc.h"
  69 #include "gromacs/pbcutil/pbc_simd.h"
  70 #include "gromacs/simd/simd.h"
  71 #include "gromacs/simd/simd_math.h"
  72 #include "gromacs/simd/vector_operations.h"
  73 #include "gromacs/utility/basedefinitions.h"
  74 #include "gromacs/utility/enumerationhelpers.h"
  75 #include "gromacs/utility/fatalerror.h"
  76 #include "gromacs/utility/real.h"
  77 #include "gromacs/utility/smalloc.h"
  78
  79 #include "listed_internal.h"
  80 #include "restcbt.h"
  81
  82 using namespace gmx; // TODO: Remove when this file is moved into gmx namespace
  83
  84 const EnumerationArray<BondedKernelFlavor, std::string> c_bondedKernelFlavorStrings = {
  85     "forces, using SIMD when available", "forces, not using SIMD",
  86     "forces, virial, and energy (ie. not using SIMD)", "forces and energy (ie. not using SIMD)"
  87 };
  88 namespace
  89 {
  90
  91 //! Type of CPU function to compute a bonded interaction.
  92 using BondedFunction = real (*)(int              nbonds,
  93                                 const t_iatom    iatoms[],
  94                                 const t_iparams  iparams[],
  95                                 const rvec       x[],
  96                                 rvec4            f[],
  97                                 rvec             fshift[],
  98                                 const t_pbc*     pbc,
  99                                 real             lambda,
 100                                 real*            dvdlambda,
 101                                 const t_mdatoms* md,
 102                                 t_fcdata*        fcd,
 103                                 int*             ddgatindex);
 104
 105 /*! \brief Mysterious CMAP coefficient matrix */
 106 const int cmap_coeff_matrix[] = {
 107     1,  0,  -3, 2,  0,  0, 0,  0,  -3, 0,  9,  -6, 2, 0,  -6, 4,  0,  0,  0, 0,  0, 0, 0,  0,
 108     3,  0,  -9, 6,  -2, 0, 6,  -4, 0,  0,  0,  0,  0, 0,  0,  0,  0,  0,  9, -6, 0, 0, -6, 4,
 109     0,  0,  3,  -2, 0,  0, 0,  0,  0,  0,  -9, 6,  0, 0,  6,  -4, 0,  0,  0, 0,  1, 0, -3, 2,
 110     -2, 0,  6,  -4, 1,  0, -3, 2,  0,  0,  0,  0,  0, 0,  0,  0,  -1, 0,  3, -2, 1, 0, -3, 2,
 111     0,  0,  0,  0,  0,  0, 0,  0,  0,  0,  -3, 2,  0, 0,  3,  -2, 0,  0,  0, 0,  0, 0, 3,  -2,
 112     0,  0,  -6, 4,  0,  0, 3,  -2, 0,  1,  -2, 1,  0, 0,  0,  0,  0,  -3, 6, -3, 0, 2, -4, 2,
 113     0,  0,  0,  0,  0,  0, 0,  0,  0,  3,  -6, 3,  0, -2, 4,  -2, 0,  0,  0, 0,  0, 0, 0,  0,
 114     0,  0,  -3, 3,  0,  0, 2,  -2, 0,  0,  -1, 1,  0, 0,  0,  0,  0,  0,  3, -3, 0, 0, -2, 2,
 115     0,  0,  0,  0,  0,  1, -2, 1,  0,  -2, 4,  -2, 0, 1,  -2, 1,  0,  0,  0, 0,  0, 0, 0,  0,
 116     0,  -1, 2,  -1, 0,  1, -2, 1,  0,  0,  0,  0,  0, 0,  0,  0,  0,  0,  1, -1, 0, 0, -1, 1,
 117     0,  0,  0,  0,  0,  0, -1, 1,  0,  0,  2,  -2, 0, 0,  -1, 1
 118 };
 119
 120
 121 /*! \brief Compute dx = xi - xj, modulo PBC if non-NULL
 122  *
 123  * \todo This kind of code appears in many places. Consolidate it */
 124 int pbc_rvec_sub(const t_pbc* pbc, const rvec xi, const rvec xj, rvec dx)
 125 {
 126     if (pbc)
 127     {
 128         return pbc_dx_aiuc(pbc, xi, xj, dx);
 129     }
 130     else
 131     {
 132         rvec_sub(xi, xj, dx);
 133         return CENTRAL;
 134     }
 135 }
 136
 137 } // namespace
 138
 139 //! \cond
 140 real bond_angle(const rvec   xi,
 141                 const rvec   xj,
 142                 const rvec   xk,
 143                 const t_pbc* pbc,
 144                 rvec         r_ij,
 145                 rvec         r_kj,
 146                 real*        costh,
 147                 int*         t1,
 148                 int*         t2)
 149 /* Return value is the angle between the bonds i-j and j-k */
 150 {
 151     /* 41 FLOPS */
 152     real th;
 153
 154     *t1 = pbc_rvec_sub(pbc, xi, xj, r_ij); /*  3                */
 155     *t2 = pbc_rvec_sub(pbc, xk, xj, r_kj); /*  3                */
 156
 157     *costh = cos_angle(r_ij, r_kj); /* 25               */
 158     th     = std::acos(*costh);     /* 10               */
 159     /* 41 TOTAL */
 160     return th;
 161 }
 162
 163 real dih_angle(const rvec   xi,
 164                const rvec   xj,
 165                const rvec   xk,
 166                const rvec   xl,
 167                const t_pbc* pbc,
 168                rvec         r_ij,
 169                rvec         r_kj,
 170                rvec         r_kl,
 171                rvec         m,
 172                rvec         n,
 173                int*         t1,
 174                int*         t2,
 175                int*         t3)
 176 {
 177     *t1 = pbc_rvec_sub(pbc, xi, xj, r_ij); /*  3        */
 178     *t2 = pbc_rvec_sub(pbc, xk, xj, r_kj); /*  3                */
 179     *t3 = pbc_rvec_sub(pbc, xk, xl, r_kl); /*  3                */
 180
 181     cprod(r_ij, r_kj, m);        /*  9        */
 182     cprod(r_kj, r_kl, n);        /*  9          */
 183     real phi  = gmx_angle(m, n); /* 49 (assuming 25 for atan2) */
 184     real ipr  = iprod(r_ij, n);  /*  5        */
 185     real sign = (ipr < 0.0) ? -1.0 : 1.0;
 186     phi       = sign * phi; /*  1               */
 187     /* 82 TOTAL */
 188     return phi;
 189 }
 190 //! \endcond
 191
 192 void make_dp_periodic(real* dp) /* 1 flop? */
 193 {
 194     /* dp cannot be outside (-pi,pi) */
 195     if (*dp >= M_PI)
 196     {
 197         *dp -= 2 * M_PI;
 198     }
 199     else if (*dp < -M_PI)
 200     {
 201         *dp += 2 * M_PI;
 202     }
 203 }
 204
 205 namespace
 206 {
 207
 208 /*! \brief Spreads and accumulates the bonded forces to the two atoms and adds the virial contribution when needed
 209  *
 210  * \p shiftIndex is used as the periodic shift.
 211  */
 212 template<BondedKernelFlavor flavor>
 213 inline void spreadBondForces(const real bondForce,
 214                              const rvec dx,
 215                              const int  ai,
 216                              const int  aj,
 217                              rvec4*     f,
 218                              int        shiftIndex,
 219                              rvec*      fshift)
 220 {
 221     for (int m = 0; m < DIM; m++) /*  15          */
 222     {
 223         const real fij = bondForce * dx[m];
 224         f[ai][m] += fij;
 225         f[aj][m] -= fij;
 226         if (computeVirial(flavor))
 227         {
 228             fshift[shiftIndex][m] += fij;
 229             fshift[CENTRAL][m] -= fij;
 230         }
 231     }
 232 }
 233
 234 /*! \brief Morse potential bond
 235  *
 236  * By Frank Everdij. Three parameters needed:
 237  *
 238  * b0 = equilibrium distance in nm
 239  * be = beta in nm^-1 (actually, it's nu_e*Sqrt(2*pi*pi*mu/D_e))
 240  * cb = well depth in kJ/mol
 241  *
 242  * Note: the potential is referenced to be +cb at infinite separation
 243  *       and zero at the equilibrium distance!
 244  */
 245 template<BondedKernelFlavor flavor>
 246 real morse_bonds(int             nbonds,
 247                  const t_iatom   forceatoms[],
 248                  const t_iparams forceparams[],
 249                  const rvec      x[],
 250                  rvec4           f[],
 251                  rvec            fshift[],
 252                  const t_pbc*    pbc,
 253                  real            lambda,
 254                  real*           dvdlambda,
 255                  const t_mdatoms gmx_unused* md,
 256                  t_fcdata gmx_unused* fcd,
 257                  int gmx_unused* global_atom_index)
 258 {
 259     const real one = 1.0;
 260     const real two = 2.0;
 261     real       dr, dr2, temp, omtemp, cbomtemp, fbond, vbond, vtot;
 262     real       b0, be, cb, b0A, beA, cbA, b0B, beB, cbB, L1;
 263     rvec       dx;
 264     int        i, ki, type, ai, aj;
 265
 266     vtot = 0.0;
 267     for (i = 0; (i < nbonds);)
 268     {
 269         type = forceatoms[i++];
 270         ai   = forceatoms[i++];
 271         aj   = forceatoms[i++];
 272
 273         b0A = forceparams[type].morse.b0A;
 274         beA = forceparams[type].morse.betaA;
 275         cbA = forceparams[type].morse.cbA;
 276
 277         b0B = forceparams[type].morse.b0B;
 278         beB = forceparams[type].morse.betaB;
 279         cbB = forceparams[type].morse.cbB;
 280
 281         L1 = one - lambda;            /* 1 */
 282         b0 = L1 * b0A + lambda * b0B; /* 3 */
 283         be = L1 * beA + lambda * beB; /* 3 */
 284         cb = L1 * cbA + lambda * cbB; /* 3 */
 285
 286         ki   = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /*   3          */
 287         dr2  = iprod(dx, dx);                       /*   5          */
 288         dr   = dr2 * gmx::invsqrt(dr2);             /*  10          */
 289         temp = std::exp(-be * (dr - b0));           /*  12          */
 290
 291         if (temp == one)
 292         {
 293             /* bonds are constrainted. This may _not_ include bond constraints if they are lambda dependent */
 294             *dvdlambda += cbB - cbA;
 295             continue;
 296         }
 297
 298         omtemp   = one - temp;                                      /*   1          */
 299         cbomtemp = cb * omtemp;                                     /*   1          */
 300         vbond    = cbomtemp * omtemp;                               /*   1          */
 301         fbond    = -two * be * temp * cbomtemp * gmx::invsqrt(dr2); /*   9          */
 302         vtot += vbond;                                              /*   1          */
 303
 304         *dvdlambda += (cbB - cbA) * omtemp * omtemp
 305                       - (2 - 2 * omtemp) * omtemp * cb * ((b0B - b0A) * be - (beB - beA) * (dr - b0)); /* 15 */
 306
 307         spreadBondForces<flavor>(fbond, dx, ai, aj, f, ki, fshift); /* 15 */
 308     }                                                               /*  83 TOTAL    */
 309     return vtot;
 310 }
 311
 312 //! \cond
 313 template<BondedKernelFlavor flavor>
 314 real cubic_bonds(int             nbonds,
 315                  const t_iatom   forceatoms[],
 316                  const t_iparams forceparams[],
 317                  const rvec      x[],
 318                  rvec4           f[],
 319                  rvec            fshift[],
 320                  const t_pbc*    pbc,
 321                  real gmx_unused lambda,
 322                  real gmx_unused* dvdlambda,
 323                  const t_mdatoms gmx_unused* md,
 324                  t_fcdata gmx_unused* fcd,
 325                  int gmx_unused* global_atom_index)
 326 {
 327     const real three = 3.0;
 328     const real two   = 2.0;
 329     real       kb, b0, kcub;
 330     real       dr, dr2, dist, kdist, kdist2, fbond, vbond, vtot;
 331     rvec       dx;
 332     int        i, ki, type, ai, aj;
 333
 334     vtot = 0.0;
 335     for (i = 0; (i < nbonds);)
 336     {
 337         type = forceatoms[i++];
 338         ai   = forceatoms[i++];
 339         aj   = forceatoms[i++];
 340
 341         b0   = forceparams[type].cubic.b0;
 342         kb   = forceparams[type].cubic.kb;
 343         kcub = forceparams[type].cubic.kcub;
 344
 345         ki  = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /*   3          */
 346         dr2 = iprod(dx, dx);                       /*   5          */
 347
 348         if (dr2 == 0.0)
 349         {
 350             continue;
 351         }
 352
 353         dr     = dr2 * gmx::invsqrt(dr2); /*  10          */
 354         dist   = dr - b0;
 355         kdist  = kb * dist;
 356         kdist2 = kdist * dist;
 357
 358         vbond = kdist2 + kcub * kdist2 * dist;
 359         fbond = -(two * kdist + three * kdist2 * kcub) / dr;
 360
 361         vtot += vbond; /* 21 */
 362
 363         spreadBondForces<flavor>(fbond, dx, ai, aj, f, ki, fshift); /* 15 */
 364     }                                                               /*  54 TOTAL    */
 365     return vtot;
 366 }
 367
 368 template<BondedKernelFlavor flavor>
 369 real FENE_bonds(int             nbonds,
 370                 const t_iatom   forceatoms[],
 371                 const t_iparams forceparams[],
 372                 const rvec      x[],
 373                 rvec4           f[],
 374                 rvec            fshift[],
 375                 const t_pbc*    pbc,
 376                 real gmx_unused lambda,
 377                 real gmx_unused* dvdlambda,
 378                 const t_mdatoms gmx_unused* md,
 379                 t_fcdata gmx_unused* fcd,
 380                 int*                 global_atom_index)
 381 {
 382     const real half = 0.5;
 383     const real one  = 1.0;
 384     real       bm, kb;
 385     real       dr2, bm2, omdr2obm2, fbond, vbond, vtot;
 386     rvec       dx;
 387     int        i, ki, type, ai, aj;
 388
 389     vtot = 0.0;
 390     for (i = 0; (i < nbonds);)
 391     {
 392         type = forceatoms[i++];
 393         ai   = forceatoms[i++];
 394         aj   = forceatoms[i++];
 395
 396         bm = forceparams[type].fene.bm;
 397         kb = forceparams[type].fene.kb;
 398
 399         ki  = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /*   3          */
 400         dr2 = iprod(dx, dx);                       /*   5          */
 401
 402         if (dr2 == 0.0)
 403         {
 404             continue;
 405         }
 406
 407         bm2 = bm * bm;
 408
 409         if (dr2 >= bm2)
 410         {
 411             gmx_fatal(FARGS, "r^2 (%f) >= bm^2 (%f) in FENE bond between atoms %d and %d", dr2, bm2,
 412                       glatnr(global_atom_index, ai), glatnr(global_atom_index, aj));
 413         }
 414
 415         omdr2obm2 = one - dr2 / bm2;
 416
 417         vbond = -half * kb * bm2 * std::log(omdr2obm2);
 418         fbond = -kb / omdr2obm2;
 419
 420         vtot += vbond; /* 35 */
 421
 422         spreadBondForces<flavor>(fbond, dx, ai, aj, f, ki, fshift); /* 15 */
 423     }                                                               /*  58 TOTAL    */
 424     return vtot;
 425 }
 426
 427 real harmonic(real kA, real kB, real xA, real xB, real x, real lambda, real* V, real* F)
 428 {
 429     const real half = 0.5;
 430     real       L1, kk, x0, dx, dx2;
 431     real       v, f, dvdlambda;
 432
 433     L1 = 1.0 - lambda;
 434     kk = L1 * kA + lambda * kB;
 435     x0 = L1 * xA + lambda * xB;
 436
 437     dx  = x - x0;
 438     dx2 = dx * dx;
 439
 440     f         = -kk * dx;
 441     v         = half * kk * dx2;
 442     dvdlambda = half * (kB - kA) * dx2 + (xA - xB) * kk * dx;
 443
 444     *F = f;
 445     *V = v;
 446
 447     return dvdlambda;
 448
 449     /* That was 19 flops */
 450 }
 451
 452 template<BondedKernelFlavor flavor>
 453 std::enable_if_t<flavor != BondedKernelFlavor::ForcesSimdWhenAvailable || !GMX_SIMD_HAVE_REAL, real>
 454 bonds(int             nbonds,
 455       const t_iatom   forceatoms[],
 456       const t_iparams forceparams[],
 457       const rvec      x[],
 458       rvec4           f[],
 459       rvec            fshift[],
 460       const t_pbc*    pbc,
 461       real            lambda,
 462       real*           dvdlambda,
 463       const t_mdatoms gmx_unused* md,
 464       t_fcdata gmx_unused* fcd,
 465       int gmx_unused* global_atom_index)
 466 {
 467     int  i, ki, ai, aj, type;
 468     real dr, dr2, fbond, vbond, vtot;
 469     rvec dx;
 470
 471     vtot = 0.0;
 472     for (i = 0; (i < nbonds);)
 473     {
 474         type = forceatoms[i++];
 475         ai   = forceatoms[i++];
 476         aj   = forceatoms[i++];
 477
 478         ki  = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /*   3      */
 479         dr2 = iprod(dx, dx);                       /*   5               */
 480         dr  = std::sqrt(dr2);                      /*  10               */
 481
 482         *dvdlambda += harmonic(forceparams[type].harmonic.krA, forceparams[type].harmonic.krB,
 483                                forceparams[type].harmonic.rA, forceparams[type].harmonic.rB, dr,
 484                                lambda, &vbond, &fbond); /*  19  */
 485
 486         if (dr2 == 0.0)
 487         {
 488             continue;
 489         }
 490
 491
 492         vtot += vbond;              /* 1*/
 493         fbond *= gmx::invsqrt(dr2); /*   6              */
 494
 495         spreadBondForces<flavor>(fbond, dx, ai, aj, f, ki, fshift); /* 15 */
 496     }                                                               /* 59 TOTAL */
 497     return vtot;
 498 }
 499
 500 #if GMX_SIMD_HAVE_REAL
 501
 502 /*! \brief Computes forces (only) for harmonic bonds using SIMD intrinsics
 503  *
 504  * As plain-C bonds(), but using SIMD to calculate many bonds at once.
 505  * This routines does not calculate energies and shift forces.
 506  */
 507 template<BondedKernelFlavor flavor>
 508 std::enable_if_t<flavor == BondedKernelFlavor::ForcesSimdWhenAvailable, real>
 509 bonds(int             nbonds,
 510       const t_iatom   forceatoms[],
 511       const t_iparams forceparams[],
 512       const rvec      x[],
 513       rvec4           f[],
 514       rvec gmx_unused fshift[],
 515       const t_pbc*    pbc,
 516       real gmx_unused lambda,
 517       real gmx_unused* dvdlambda,
 518       const t_mdatoms gmx_unused* md,
 519       t_fcdata gmx_unused* fcd,
 520       int gmx_unused* global_atom_index)
 521 {
 522     constexpr int                            nfa1 = 3;
 523     alignas(GMX_SIMD_ALIGNMENT) std::int32_t ai[GMX_SIMD_REAL_WIDTH];
 524     alignas(GMX_SIMD_ALIGNMENT) std::int32_t aj[GMX_SIMD_REAL_WIDTH];
 525     alignas(GMX_SIMD_ALIGNMENT) real         coeff[2 * GMX_SIMD_REAL_WIDTH];
 526
 527     alignas(GMX_SIMD_ALIGNMENT) real pbc_simd[9 * GMX_SIMD_REAL_WIDTH];
 528
 529     set_pbc_simd(pbc, pbc_simd);
 530
 531     const SimdReal real_eps = SimdReal(GMX_REAL_EPS);
 532
 533     /* nbonds is the number of bonds times nfa1, here we step GMX_SIMD_REAL_WIDTH angles */
 534     for (int i = 0; i < nbonds; i += GMX_SIMD_REAL_WIDTH * nfa1)
 535     {
 536         /* Collect atoms for GMX_SIMD_REAL_WIDTH angles.
 537          * iu indexes into forceatoms, we should not let iu go beyond nbonds.
 538          */
 539         int iu = i;
 540         for (int s = 0; s < GMX_SIMD_REAL_WIDTH; s++)
 541         {
 542             const int type = forceatoms[iu];
 543             ai[s]          = forceatoms[iu + 1];
 544             aj[s]          = forceatoms[iu + 2];
 545
 546             /* At the end fill the arrays with the last atoms and 0 params */
 547             if (i + s * nfa1 < nbonds)
 548             {
 549                 coeff[s]                       = forceparams[type].harmonic.krA;
 550                 coeff[GMX_SIMD_REAL_WIDTH + s] = forceparams[type].harmonic.rA;
 551
 552                 if (iu + nfa1 < nbonds)
 553                 {
 554                     iu += nfa1;
 555                 }
 556             }
 557             else
 558             {
 559                 coeff[s]                       = 0;
 560                 coeff[GMX_SIMD_REAL_WIDTH + s] = 0;
 561             }
 562         }
 563
 564         /* Store the non PBC corrected distances packed and aligned */
 565         SimdReal xi, yi, zi;
 566         SimdReal xj, yj, zj;
 567         gatherLoadUTranspose<3>(reinterpret_cast<const real*>(x), ai, &xi, &yi, &zi);
 568         gatherLoadUTranspose<3>(reinterpret_cast<const real*>(x), aj, &xj, &yj, &zj);
 569         SimdReal rij_x = xi - xj;
 570         SimdReal rij_y = yi - yj;
 571         SimdReal rij_z = zi - zj;
 572
 573         pbc_correct_dx_simd(&rij_x, &rij_y, &rij_z, pbc_simd);
 574
 575         const SimdReal dist2 = rij_x * rij_x + rij_y * rij_y + rij_z * rij_z;
 576         // Here we avoid sqrt(0), the force will be zero because rij=0
 577         const SimdReal invDist = invsqrt(max(dist2, real_eps));
 578
 579         const SimdReal k  = load<SimdReal>(coeff);
 580         const SimdReal r0 = load<SimdReal>(coeff + GMX_SIMD_REAL_WIDTH);
 581
 582         // Compute the force divided by the distance
 583         const SimdReal forceOverR = -k * (dist2 * invDist - r0) * invDist;
 584
 585         const SimdReal f_x = forceOverR * rij_x;
 586         const SimdReal f_y = forceOverR * rij_y;
 587         const SimdReal f_z = forceOverR * rij_z;
 588
 589         transposeScatterIncrU<4>(reinterpret_cast<real*>(f), ai, f_x, f_y, f_z);
 590         transposeScatterDecrU<4>(reinterpret_cast<real*>(f), aj, f_x, f_y, f_z);
 591     }
 592
 593     return 0;
 594 }
 595
 596 #endif // GMX_SIMD_HAVE_REAL
 597
 598 template<BondedKernelFlavor flavor>
 599 real restraint_bonds(int             nbonds,
 600                      const t_iatom   forceatoms[],
 601                      const t_iparams forceparams[],
 602                      const rvec      x[],
 603                      rvec4           f[],
 604                      rvec            fshift[],
 605                      const t_pbc*    pbc,
 606                      real            lambda,
 607                      real*           dvdlambda,
 608                      const t_mdatoms gmx_unused* md,
 609                      t_fcdata gmx_unused* fcd,
 610                      int gmx_unused* global_atom_index)
 611 {
 612     int  i, ki, ai, aj, type;
 613     real dr, dr2, fbond, vbond, vtot;
 614     real L1;
 615     real low, dlow, up1, dup1, up2, dup2, k, dk;
 616     real drh, drh2;
 617     rvec dx;
 618
 619     L1 = 1.0 - lambda;
 620
 621     vtot = 0.0;
 622     for (i = 0; (i < nbonds);)
 623     {
 624         type = forceatoms[i++];
 625         ai   = forceatoms[i++];
 626         aj   = forceatoms[i++];
 627
 628         ki  = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /*   3      */
 629         dr2 = iprod(dx, dx);                       /*   5               */
 630         dr  = dr2 * gmx::invsqrt(dr2);             /*  10               */
 631
 632         low  = L1 * forceparams[type].restraint.lowA + lambda * forceparams[type].restraint.lowB;
 633         dlow = -forceparams[type].restraint.lowA + forceparams[type].restraint.lowB;
 634         up1  = L1 * forceparams[type].restraint.up1A + lambda * forceparams[type].restraint.up1B;
 635         dup1 = -forceparams[type].restraint.up1A + forceparams[type].restraint.up1B;
 636         up2  = L1 * forceparams[type].restraint.up2A + lambda * forceparams[type].restraint.up2B;
 637         dup2 = -forceparams[type].restraint.up2A + forceparams[type].restraint.up2B;
 638         k    = L1 * forceparams[type].restraint.kA + lambda * forceparams[type].restraint.kB;
 639         dk   = -forceparams[type].restraint.kA + forceparams[type].restraint.kB;
 640         /* 24 */
 641
 642         if (dr < low)
 643         {
 644             drh   = dr - low;
 645             drh2  = drh * drh;
 646             vbond = 0.5 * k * drh2;
 647             fbond = -k * drh;
 648             *dvdlambda += 0.5 * dk * drh2 - k * dlow * drh;
 649         } /* 11 */
 650         else if (dr <= up1)
 651         {
 652             vbond = 0;
 653             fbond = 0;
 654         }
 655         else if (dr <= up2)
 656         {
 657             drh   = dr - up1;
 658             drh2  = drh * drh;
 659             vbond = 0.5 * k * drh2;
 660             fbond = -k * drh;
 661             *dvdlambda += 0.5 * dk * drh2 - k * dup1 * drh;
 662         } /* 11 */
 663         else
 664         {
 665             drh   = dr - up2;
 666             vbond = k * (up2 - up1) * (0.5 * (up2 - up1) + drh);
 667             fbond = -k * (up2 - up1);
 668             *dvdlambda += dk * (up2 - up1) * (0.5 * (up2 - up1) + drh)
 669                           + k * (dup2 - dup1) * (up2 - up1 + drh) - k * (up2 - up1) * dup2;
 670         }
 671
 672         if (dr2 == 0.0)
 673         {
 674             continue;
 675         }
 676
 677         vtot += vbond;              /* 1*/
 678         fbond *= gmx::invsqrt(dr2); /*   6              */
 679
 680         spreadBondForces<flavor>(fbond, dx, ai, aj, f, ki, fshift); /* 15 */
 681     }                                                               /* 59 TOTAL */
 682
 683     return vtot;
 684 }
 685
 686 template<BondedKernelFlavor flavor>
 687 real polarize(int              nbonds,
 688               const t_iatom    forceatoms[],
 689               const t_iparams  forceparams[],
 690               const rvec       x[],
 691               rvec4            f[],
 692               rvec             fshift[],
 693               const t_pbc*     pbc,
 694               real             lambda,
 695               real*            dvdlambda,
 696               const t_mdatoms* md,
 697               t_fcdata gmx_unused* fcd,
 698               int gmx_unused* global_atom_index)
 699 {
 700     int  i, ki, ai, aj, type;
 701     real dr, dr2, fbond, vbond, vtot, ksh;
 702     rvec dx;
 703
 704     vtot = 0.0;
 705     for (i = 0; (i < nbonds);)
 706     {
 707         type = forceatoms[i++];
 708         ai   = forceatoms[i++];
 709         aj   = forceatoms[i++];
 710         ksh  = gmx::square(md->chargeA[aj]) * ONE_4PI_EPS0 / forceparams[type].polarize.alpha;
 711
 712         ki  = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /*   3      */
 713         dr2 = iprod(dx, dx);                       /*   5               */
 714         dr  = std::sqrt(dr2);                      /*  10               */
 715
 716         *dvdlambda += harmonic(ksh, ksh, 0, 0, dr, lambda, &vbond, &fbond); /*  19  */
 717
 718         if (dr2 == 0.0)
 719         {
 720             continue;
 721         }
 722
 723         vtot += vbond;              /* 1*/
 724         fbond *= gmx::invsqrt(dr2); /*   6              */
 725
 726         spreadBondForces<flavor>(fbond, dx, ai, aj, f, ki, fshift); /* 15 */
 727     }                                                               /* 59 TOTAL */
 728     return vtot;
 729 }
 730
 731 template<BondedKernelFlavor flavor>
 732 real anharm_polarize(int              nbonds,
 733                      const t_iatom    forceatoms[],
 734                      const t_iparams  forceparams[],
 735                      const rvec       x[],
 736                      rvec4            f[],
 737                      rvec             fshift[],
 738                      const t_pbc*     pbc,
 739                      real             lambda,
 740                      real*            dvdlambda,
 741                      const t_mdatoms* md,
 742                      t_fcdata gmx_unused* fcd,
 743                      int gmx_unused* global_atom_index)
 744 {
 745     int  i, ki, ai, aj, type;
 746     real dr, dr2, fbond, vbond, vtot, ksh, khyp, drcut, ddr, ddr3;
 747     rvec dx;
 748
 749     vtot = 0.0;
 750     for (i = 0; (i < nbonds);)
 751     {
 752         type = forceatoms[i++];
 753         ai   = forceatoms[i++];
 754         aj   = forceatoms[i++];
 755         ksh = gmx::square(md->chargeA[aj]) * ONE_4PI_EPS0 / forceparams[type].anharm_polarize.alpha; /* 7*/
 756         khyp  = forceparams[type].anharm_polarize.khyp;
 757         drcut = forceparams[type].anharm_polarize.drcut;
 758
 759         ki  = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /*   3      */
 760         dr2 = iprod(dx, dx);                       /*   5               */
 761         dr  = dr2 * gmx::invsqrt(dr2);             /*  10               */
 762
 763         *dvdlambda += harmonic(ksh, ksh, 0, 0, dr, lambda, &vbond, &fbond); /*  19  */
 764
 765         if (dr2 == 0.0)
 766         {
 767             continue;
 768         }
 769
 770         if (dr > drcut)
 771         {
 772             ddr  = dr - drcut;
 773             ddr3 = ddr * ddr * ddr;
 774             vbond += khyp * ddr * ddr3;
 775             fbond -= 4 * khyp * ddr3;
 776         }
 777         fbond *= gmx::invsqrt(dr2); /*   6              */
 778         vtot += vbond;              /* 1*/
 779
 780         spreadBondForces<flavor>(fbond, dx, ai, aj, f, ki, fshift); /* 15 */
 781     }                                                               /* 72 TOTAL */
 782     return vtot;
 783 }
 784
 785 template<BondedKernelFlavor flavor>
 786 real water_pol(int             nbonds,
 787                const t_iatom   forceatoms[],
 788                const t_iparams forceparams[],
 789                const rvec      x[],
 790                rvec4           f[],
 791                rvec gmx_unused fshift[],
 792                const t_pbc gmx_unused* pbc,
 793                real gmx_unused lambda,
 794                real gmx_unused* dvdlambda,
 795                const t_mdatoms gmx_unused* md,
 796                t_fcdata gmx_unused* fcd,
 797                int gmx_unused* global_atom_index)
 798 {
 799     /* This routine implements anisotropic polarizibility for water, through
 800      * a shell connected to a dummy with spring constant that differ in the
 801      * three spatial dimensions in the molecular frame.
 802      */
 803     int  i, m, aO, aH1, aH2, aD, aS, type, type0, ki;
 804     rvec dOH1, dOH2, dHH, dOD, dDS, nW, kk, dx, kdx, proj;
 805     real vtot, fij, r_HH, r_OD, r_nW, tx, ty, tz, qS;
 806
 807     vtot = 0.0;
 808     if (nbonds > 0)
 809     {
 810         type0  = forceatoms[0];
 811         aS     = forceatoms[5];
 812         qS     = md->chargeA[aS];
 813         kk[XX] = gmx::square(qS) * ONE_4PI_EPS0 / forceparams[type0].wpol.al_x;
 814         kk[YY] = gmx::square(qS) * ONE_4PI_EPS0 / forceparams[type0].wpol.al_y;
 815         kk[ZZ] = gmx::square(qS) * ONE_4PI_EPS0 / forceparams[type0].wpol.al_z;
 816         r_HH   = 1.0 / forceparams[type0].wpol.rHH;
 817         for (i = 0; (i < nbonds); i += 6)
 818         {
 819             type = forceatoms[i];
 820             if (type != type0)
 821             {
 822                 gmx_fatal(FARGS, "Sorry, type = %d, type0 = %d, file = %s, line = %d", type, type0,
 823                           __FILE__, __LINE__);
 824             }
 825             aO  = forceatoms[i + 1];
 826             aH1 = forceatoms[i + 2];
 827             aH2 = forceatoms[i + 3];
 828             aD  = forceatoms[i + 4];
 829             aS  = forceatoms[i + 5];
 830
 831             /* Compute vectors describing the water frame */
 832             pbc_rvec_sub(pbc, x[aH1], x[aO], dOH1);
 833             pbc_rvec_sub(pbc, x[aH2], x[aO], dOH2);
 834             pbc_rvec_sub(pbc, x[aH2], x[aH1], dHH);
 835             pbc_rvec_sub(pbc, x[aD], x[aO], dOD);
 836             ki = pbc_rvec_sub(pbc, x[aS], x[aD], dDS);
 837             cprod(dOH1, dOH2, nW);
 838
 839             /* Compute inverse length of normal vector
 840              * (this one could be precomputed, but I'm too lazy now)
 841              */
 842             r_nW = gmx::invsqrt(iprod(nW, nW));
 843             /* This is for precision, but does not make a big difference,
 844              * it can go later.
 845              */
 846             r_OD = gmx::invsqrt(iprod(dOD, dOD));
 847
 848             /* Normalize the vectors in the water frame */
 849             svmul(r_nW, nW, nW);
 850             svmul(r_HH, dHH, dHH);
 851             svmul(r_OD, dOD, dOD);
 852
 853             /* Compute displacement of shell along components of the vector */
 854             dx[ZZ] = iprod(dDS, dOD);
 855             /* Compute projection on the XY plane: dDS - dx[ZZ]*dOD */
 856             for (m = 0; (m < DIM); m++)
 857             {
 858                 proj[m] = dDS[m] - dx[ZZ] * dOD[m];
 859             }
 860
 861             /*dx[XX] = iprod(dDS,nW);
 862                dx[YY] = iprod(dDS,dHH);*/
 863             dx[XX] = iprod(proj, nW);
 864             for (m = 0; (m < DIM); m++)
 865             {
 866                 proj[m] -= dx[XX] * nW[m];
 867             }
 868             dx[YY] = iprod(proj, dHH);
 869             /* Now compute the forces and energy */
 870             kdx[XX] = kk[XX] * dx[XX];
 871             kdx[YY] = kk[YY] * dx[YY];
 872             kdx[ZZ] = kk[ZZ] * dx[ZZ];
 873             vtot += iprod(dx, kdx);
 874
 875             for (m = 0; (m < DIM); m++)
 876             {
 877                 /* This is a tensor operation but written out for speed */
 878                 tx  = nW[m] * kdx[XX];
 879                 ty  = dHH[m] * kdx[YY];
 880                 tz  = dOD[m] * kdx[ZZ];
 881                 fij = -tx - ty - tz;
 882                 f[aS][m] += fij;
 883                 f[aD][m] -= fij;
 884                 if (computeVirial(flavor))
 885                 {
 886                     fshift[ki][m] += fij;
 887                     fshift[CENTRAL][m] -= fij;
 888                 }
 889             }
 890         }
 891     }
 892     return 0.5 * vtot;
 893 }
 894
 895 template<BondedKernelFlavor flavor>
 896 static real
 897 do_1_thole(const rvec xi, const rvec xj, rvec fi, rvec fj, const t_pbc* pbc, real qq, rvec fshift[], real afac)
 898 {
 899     rvec r12;
 900     real r12sq, r12_1, r12bar, v0, v1, fscal, ebar, fff;
 901     int  m, t;
 902
 903     t = pbc_rvec_sub(pbc, xi, xj, r12); /*  3 */
 904
 905     r12sq  = iprod(r12, r12);                                                     /*  5 */
 906     r12_1  = gmx::invsqrt(r12sq);                                                 /*  5 */
 907     r12bar = afac / r12_1;                                                        /*  5 */
 908     v0     = qq * ONE_4PI_EPS0 * r12_1;                                           /*  2 */
 909     ebar   = std::exp(-r12bar);                                                   /*  5 */
 910     v1     = (1 - (1 + 0.5 * r12bar) * ebar);                                     /*  4 */
 911     fscal  = ((v0 * r12_1) * v1 - v0 * 0.5 * afac * ebar * (r12bar + 1)) * r12_1; /* 9 */
 912
 913     for (m = 0; (m < DIM); m++)
 914     {
 915         fff = fscal * r12[m];
 916         fi[m] += fff;
 917         fj[m] -= fff;
 918         if (computeVirial(flavor))
 919         {
 920             fshift[t][m] += fff;
 921             fshift[CENTRAL][m] -= fff;
 922         }
 923     } /* 15 */
 924
 925     return v0 * v1; /* 1 */
 926     /* 54 */
 927 }
 928
 929 template<BondedKernelFlavor flavor>
 930 real thole_pol(int             nbonds,
 931                const t_iatom   forceatoms[],
 932                const t_iparams forceparams[],
 933                const rvec      x[],
 934                rvec4           f[],
 935                rvec            fshift[],
 936                const t_pbc*    pbc,
 937                real gmx_unused lambda,
 938                real gmx_unused* dvdlambda,
 939                const t_mdatoms* md,
 940                t_fcdata gmx_unused* fcd,
 941                int gmx_unused* global_atom_index)
 942 {
 943     /* Interaction between two pairs of particles with opposite charge */
 944     int  i, type, a1, da1, a2, da2;
 945     real q1, q2, qq, a, al1, al2, afac;
 946     real V = 0;
 947
 948     for (i = 0; (i < nbonds);)
 949     {
 950         type = forceatoms[i++];
 951         a1   = forceatoms[i++];
 952         da1  = forceatoms[i++];
 953         a2   = forceatoms[i++];
 954         da2  = forceatoms[i++];
 955         q1   = md->chargeA[da1];
 956         q2   = md->chargeA[da2];
 957         a    = forceparams[type].thole.a;
 958         al1  = forceparams[type].thole.alpha1;
 959         al2  = forceparams[type].thole.alpha2;
 960         qq   = q1 * q2;
 961         afac = a * gmx::invsixthroot(al1 * al2);
 962         V += do_1_thole<flavor>(x[a1], x[a2], f[a1], f[a2], pbc, qq, fshift, afac);
 963         V += do_1_thole<flavor>(x[da1], x[a2], f[da1], f[a2], pbc, -qq, fshift, afac);
 964         V += do_1_thole<flavor>(x[a1], x[da2], f[a1], f[da2], pbc, -qq, fshift, afac);
 965         V += do_1_thole<flavor>(x[da1], x[da2], f[da1], f[da2], pbc, qq, fshift, afac);
 966     }
 967     /* 290 flops */
 968     return V;
 969 }
 970
 971 // Avoid gcc 386 -O3 code generation bug in this function (see Issue
 972 // #3205 for more information)
 973 #if defined(__GNUC__) && defined(__i386__) && defined(__OPTIMIZE__)
 974 #    pragma GCC push_options
 975 #    pragma GCC optimize("O1")
 976 #    define avoid_gcc_i386_o3_code_generation_bug
 977 #endif
 978
 979 template<BondedKernelFlavor flavor>
 980 std::enable_if_t<flavor != BondedKernelFlavor::ForcesSimdWhenAvailable || !GMX_SIMD_HAVE_REAL, real>
 981 angles(int             nbonds,
 982        const t_iatom   forceatoms[],
 983        const t_iparams forceparams[],
 984        const rvec      x[],
 985        rvec4           f[],
 986        rvec            fshift[],
 987        const t_pbc*    pbc,
 988        real            lambda,
 989        real*           dvdlambda,
 990        const t_mdatoms gmx_unused* md,
 991        t_fcdata gmx_unused* fcd,
 992        int gmx_unused* global_atom_index)
 993 {
 994     int  i, ai, aj, ak, t1, t2, type;
 995     rvec r_ij, r_kj;
 996     real cos_theta, cos_theta2, theta, dVdt, va, vtot;
 997
 998     vtot = 0.0;
 999     for (i = 0; i < nbonds;)
1000     {
1001         type = forceatoms[i++];
1002         ai   = forceatoms[i++];
1003         aj   = forceatoms[i++];
1004         ak   = forceatoms[i++];
1005
1006         theta = bond_angle(x[ai], x[aj], x[ak], pbc, r_ij, r_kj, &cos_theta, &t1, &t2); /*  41 */
1007
1008         *dvdlambda += harmonic(forceparams[type].harmonic.krA, forceparams[type].harmonic.krB,
1009                                forceparams[type].harmonic.rA * DEG2RAD,
1010                                forceparams[type].harmonic.rB * DEG2RAD, theta, lambda, &va, &dVdt); /*  21  */
1011         vtot += va;
1012
1013         cos_theta2 = gmx::square(cos_theta);
1014         if (cos_theta2 < 1)
1015         {
1016             int  m;
1017             real st, sth;
1018             real cik, cii, ckk;
1019             real nrkj2, nrij2;
1020             real nrkj_1, nrij_1;
1021             rvec f_i, f_j, f_k;
1022
1023             st    = dVdt * gmx::invsqrt(1 - cos_theta2); /*  12         */
1024             sth   = st * cos_theta;                      /*   1         */
1025             nrij2 = iprod(r_ij, r_ij);                   /*   5         */
1026             nrkj2 = iprod(r_kj, r_kj);                   /*   5         */
1027
1028             nrij_1 = gmx::invsqrt(nrij2); /*  10                */
1029             nrkj_1 = gmx::invsqrt(nrkj2); /*  10                */
1030
1031             cik = st * nrij_1 * nrkj_1;  /*   2         */
1032             cii = sth * nrij_1 * nrij_1; /*   2         */
1033             ckk = sth * nrkj_1 * nrkj_1; /*   2         */
1034
1035             for (m = 0; m < DIM; m++)
1036             { /*  39            */
1037                 f_i[m] = -(cik * r_kj[m] - cii * r_ij[m]);
1038                 f_k[m] = -(cik * r_ij[m] - ckk * r_kj[m]);
1039                 f_j[m] = -f_i[m] - f_k[m];
1040                 f[ai][m] += f_i[m];
1041                 f[aj][m] += f_j[m];
1042                 f[ak][m] += f_k[m];
1043             }
1044             if (computeVirial(flavor))
1045             {
1046                 rvec_inc(fshift[t1], f_i);
1047                 rvec_inc(fshift[CENTRAL], f_j);
1048                 rvec_inc(fshift[t2], f_k);
1049             }
1050         } /* 161 TOTAL  */
1051     }
1052
1053     return vtot;
1054 }
1055
1056 #ifdef avoid_gcc_i386_o3_code_generation_bug
1057 #    pragma GCC pop_options
1058 #    undef avoid_gcc_i386_o3_code_generation_bug
1059 #endif
1060
1061 #if GMX_SIMD_HAVE_REAL
1062
1063 /* As angles, but using SIMD to calculate many angles at once.
1064  * This routines does not calculate energies and shift forces.
1065  */
1066 template<BondedKernelFlavor flavor>
1067 std::enable_if_t<flavor == BondedKernelFlavor::ForcesSimdWhenAvailable, real>
1068 angles(int             nbonds,
1069        const t_iatom   forceatoms[],
1070        const t_iparams forceparams[],
1071        const rvec      x[],
1072        rvec4           f[],
1073        rvec gmx_unused fshift[],
1074        const t_pbc*    pbc,
1075        real gmx_unused lambda,
1076        real gmx_unused* dvdlambda,
1077        const t_mdatoms gmx_unused* md,
1078        t_fcdata gmx_unused* fcd,
1079        int gmx_unused* global_atom_index)
1080 {
1081     const int                                nfa1 = 4;
1082     int                                      i, iu, s;
1083     int                                      type;
1084     alignas(GMX_SIMD_ALIGNMENT) std::int32_t ai[GMX_SIMD_REAL_WIDTH];
1085     alignas(GMX_SIMD_ALIGNMENT) std::int32_t aj[GMX_SIMD_REAL_WIDTH];
1086     alignas(GMX_SIMD_ALIGNMENT) std::int32_t ak[GMX_SIMD_REAL_WIDTH];
1087     alignas(GMX_SIMD_ALIGNMENT) real         coeff[2 * GMX_SIMD_REAL_WIDTH];
1088     SimdReal                                 deg2rad_S(DEG2RAD);
1089     SimdReal                                 xi_S, yi_S, zi_S;
1090     SimdReal                                 xj_S, yj_S, zj_S;
1091     SimdReal                                 xk_S, yk_S, zk_S;
1092     SimdReal                                 k_S, theta0_S;
1093     SimdReal                                 rijx_S, rijy_S, rijz_S;
1094     SimdReal                                 rkjx_S, rkjy_S, rkjz_S;
1095     SimdReal                                 one_S(1.0);
1096     SimdReal one_min_eps_S(1.0_real - GMX_REAL_EPS); // Largest number < 1
1097
1098     SimdReal                         rij_rkj_S;
1099     SimdReal                         nrij2_S, nrij_1_S;
1100     SimdReal                         nrkj2_S, nrkj_1_S;
1101     SimdReal                         cos_S, invsin_S;
1102     SimdReal                         cos2_S;
1103     SimdReal                         theta_S;
1104     SimdReal                         st_S, sth_S;
1105     SimdReal                         cik_S, cii_S, ckk_S;
1106     SimdReal                         f_ix_S, f_iy_S, f_iz_S;
1107     SimdReal                         f_kx_S, f_ky_S, f_kz_S;
1108     alignas(GMX_SIMD_ALIGNMENT) real pbc_simd[9 * GMX_SIMD_REAL_WIDTH];
1109
1110     set_pbc_simd(pbc, pbc_simd);
1111
1112     /* nbonds is the number of angles times nfa1, here we step GMX_SIMD_REAL_WIDTH angles */
1113     for (i = 0; (i < nbonds); i += GMX_SIMD_REAL_WIDTH * nfa1)
1114     {
1115         /* Collect atoms for GMX_SIMD_REAL_WIDTH angles.
1116          * iu indexes into forceatoms, we should not let iu go beyond nbonds.
1117          */
1118         iu = i;
1119         for (s = 0; s < GMX_SIMD_REAL_WIDTH; s++)
1120         {
1121             type  = forceatoms[iu];
1122             ai[s] = forceatoms[iu + 1];
1123             aj[s] = forceatoms[iu + 2];
1124             ak[s] = forceatoms[iu + 3];
1125
1126             /* At the end fill the arrays with the last atoms and 0 params */
1127             if (i + s * nfa1 < nbonds)
1128             {
1129                 coeff[s]                       = forceparams[type].harmonic.krA;
1130                 coeff[GMX_SIMD_REAL_WIDTH + s] = forceparams[type].harmonic.rA;
1131
1132                 if (iu + nfa1 < nbonds)
1133                 {
1134                     iu += nfa1;
1135                 }
1136             }
1137             else
1138             {
1139                 coeff[s]                       = 0;
1140                 coeff[GMX_SIMD_REAL_WIDTH + s] = 0;
1141             }
1142         }
1143
1144         /* Store the non PBC corrected distances packed and aligned */
1145         gatherLoadUTranspose<3>(reinterpret_cast<const real*>(x), ai, &xi_S, &yi_S, &zi_S);
1146         gatherLoadUTranspose<3>(reinterpret_cast<const real*>(x), aj, &xj_S, &yj_S, &zj_S);
1147         gatherLoadUTranspose<3>(reinterpret_cast<const real*>(x), ak, &xk_S, &yk_S, &zk_S);
1148         rijx_S = xi_S - xj_S;
1149         rijy_S = yi_S - yj_S;
1150         rijz_S = zi_S - zj_S;
1151         rkjx_S = xk_S - xj_S;
1152         rkjy_S = yk_S - yj_S;
1153         rkjz_S = zk_S - zj_S;
1154
1155         k_S      = load<SimdReal>(coeff);
1156         theta0_S = load<SimdReal>(coeff + GMX_SIMD_REAL_WIDTH) * deg2rad_S;
1157
1158         pbc_correct_dx_simd(&rijx_S, &rijy_S, &rijz_S, pbc_simd);
1159         pbc_correct_dx_simd(&rkjx_S, &rkjy_S, &rkjz_S, pbc_simd);
1160
1161         rij_rkj_S = iprod(rijx_S, rijy_S, rijz_S, rkjx_S, rkjy_S, rkjz_S);
1162
1163         nrij2_S = norm2(rijx_S, rijy_S, rijz_S);
1164         nrkj2_S = norm2(rkjx_S, rkjy_S, rkjz_S);
1165
1166         nrij_1_S = invsqrt(nrij2_S);
1167         nrkj_1_S = invsqrt(nrkj2_S);
1168
1169         cos_S = rij_rkj_S * nrij_1_S * nrkj_1_S;
1170
1171         /* We compute cos^2 using a division instead of squaring cos_S,
1172          * as we would then get additional error contributions from
1173          * the two invsqrt operations, which get amplified by
1174          * the 1/sqrt(1-cos^2) for cos values close to 1.
1175          *
1176          * Note that the non-SIMD version of angles() avoids this issue
1177          * by computing the cosine using double precision.
1178          */
1179         cos2_S = rij_rkj_S * rij_rkj_S / (nrij2_S * nrkj2_S);
1180
1181         /* To allow for 0 and 180 degrees, we need to avoid issues when
1182          * the cosine is close to -1 and 1. For acos() the argument needs
1183          * to be in that range. For cos^2 we take the min of cos^2 and 1 - 1bit,
1184          * so we can safely compute 1/sin as 1/sqrt(1 - cos^2).
1185          */
1186         cos_S  = max(cos_S, -one_S);
1187         cos_S  = min(cos_S, one_S);
1188         cos2_S = min(cos2_S, one_min_eps_S);
1189
1190         theta_S = acos(cos_S);
1191
1192         invsin_S = invsqrt(one_S - cos2_S);
1193
1194         st_S  = k_S * (theta0_S - theta_S) * invsin_S;
1195         sth_S = st_S * cos_S;
1196
1197         cik_S = st_S * nrij_1_S * nrkj_1_S;
1198         cii_S = sth_S * nrij_1_S * nrij_1_S;
1199         ckk_S = sth_S * nrkj_1_S * nrkj_1_S;
1200
1201         f_ix_S = cii_S * rijx_S;
1202         f_ix_S = fnma(cik_S, rkjx_S, f_ix_S);
1203         f_iy_S = cii_S * rijy_S;
1204         f_iy_S = fnma(cik_S, rkjy_S, f_iy_S);
1205         f_iz_S = cii_S * rijz_S;
1206         f_iz_S = fnma(cik_S, rkjz_S, f_iz_S);
1207         f_kx_S = ckk_S * rkjx_S;
1208         f_kx_S = fnma(cik_S, rijx_S, f_kx_S);
1209         f_ky_S = ckk_S * rkjy_S;
1210         f_ky_S = fnma(cik_S, rijy_S, f_ky_S);
1211         f_kz_S = ckk_S * rkjz_S;
1212         f_kz_S = fnma(cik_S, rijz_S, f_kz_S);
1213
1214         transposeScatterIncrU<4>(reinterpret_cast<real*>(f), ai, f_ix_S, f_iy_S, f_iz_S);
1215         transposeScatterDecrU<4>(reinterpret_cast<real*>(f), aj, f_ix_S + f_kx_S, f_iy_S + f_ky_S,
1216                                  f_iz_S + f_kz_S);
1217         transposeScatterIncrU<4>(reinterpret_cast<real*>(f), ak, f_kx_S, f_ky_S, f_kz_S);
1218     }
1219
1220     return 0;
1221 }
1222
1223 #endif // GMX_SIMD_HAVE_REAL
1224
1225 template<BondedKernelFlavor flavor>
1226 real linear_angles(int             nbonds,
1227                    const t_iatom   forceatoms[],
1228                    const t_iparams forceparams[],
1229                    const rvec      x[],
1230                    rvec4           f[],
1231                    rvec            fshift[],
1232                    const t_pbc*    pbc,
1233                    real            lambda,
1234                    real*           dvdlambda,
1235                    const t_mdatoms gmx_unused* md,
1236                    t_fcdata gmx_unused* fcd,
1237                    int gmx_unused* global_atom_index)
1238 {
1239     int  i, m, ai, aj, ak, t1, t2, type;
1240     rvec f_i, f_j, f_k;
1241     real L1, kA, kB, aA, aB, dr, dr2, va, vtot, a, b, klin;
1242     rvec r_ij, r_kj, r_ik, dx;
1243
1244     L1   = 1 - lambda;
1245     vtot = 0.0;
1246     for (i = 0; (i < nbonds);)
1247     {
1248         type = forceatoms[i++];
1249         ai   = forceatoms[i++];
1250         aj   = forceatoms[i++];
1251         ak   = forceatoms[i++];
1252
1253         kA   = forceparams[type].linangle.klinA;
1254         kB   = forceparams[type].linangle.klinB;
1255         klin = L1 * kA + lambda * kB;
1256
1257         aA = forceparams[type].linangle.aA;
1258         aB = forceparams[type].linangle.aB;
1259         a  = L1 * aA + lambda * aB;
1260         b  = 1 - a;
1261
1262         t1 = pbc_rvec_sub(pbc, x[ai], x[aj], r_ij);
1263         t2 = pbc_rvec_sub(pbc, x[ak], x[aj], r_kj);
1264         rvec_sub(r_ij, r_kj, r_ik);
1265
1266         dr2 = 0;
1267         for (m = 0; (m < DIM); m++)
1268         {
1269             dr = -a * r_ij[m] - b * r_kj[m];
1270             dr2 += dr * dr;
1271             dx[m]  = dr;
1272             f_i[m] = a * klin * dr;
1273             f_k[m] = b * klin * dr;
1274             f_j[m] = -(f_i[m] + f_k[m]);
1275             f[ai][m] += f_i[m];
1276             f[aj][m] += f_j[m];
1277             f[ak][m] += f_k[m];
1278         }
1279         va = 0.5 * klin * dr2;
1280         *dvdlambda += 0.5 * (kB - kA) * dr2 + klin * (aB - aA) * iprod(dx, r_ik);
1281
1282         vtot += va;
1283
1284         if (computeVirial(flavor))
1285         {
1286             rvec_inc(fshift[t1], f_i);
1287             rvec_inc(fshift[CENTRAL], f_j);
1288             rvec_inc(fshift[t2], f_k);
1289         }
1290     } /* 57 TOTAL       */
1291     return vtot;
1292 }
1293
1294 template<BondedKernelFlavor flavor>
1295 std::enable_if_t<flavor != BondedKernelFlavor::ForcesSimdWhenAvailable || !GMX_SIMD_HAVE_REAL, real>
1296 urey_bradley(int             nbonds,
1297              const t_iatom   forceatoms[],
1298              const t_iparams forceparams[],
1299              const rvec      x[],
1300              rvec4           f[],
1301              rvec            fshift[],
1302              const t_pbc*    pbc,
1303              real            lambda,
1304              real*           dvdlambda,
1305              const t_mdatoms gmx_unused* md,
1306              t_fcdata gmx_unused* fcd,
1307              int gmx_unused* global_atom_index)
1308 {
1309     int  i, m, ai, aj, ak, t1, t2, type, ki;
1310     rvec r_ij, r_kj, r_ik;
1311     real cos_theta, cos_theta2, theta;
1312     real dVdt, va, vtot, dr, dr2, vbond, fbond, fik;
1313     real kthA, th0A, kUBA, r13A, kthB, th0B, kUBB, r13B;
1314
1315     vtot = 0.0;
1316     for (i = 0; (i < nbonds);)
1317     {
1318         type = forceatoms[i++];
1319         ai   = forceatoms[i++];
1320         aj   = forceatoms[i++];
1321         ak   = forceatoms[i++];
1322         th0A = forceparams[type].u_b.thetaA * DEG2RAD;
1323         kthA = forceparams[type].u_b.kthetaA;
1324         r13A = forceparams[type].u_b.r13A;
1325         kUBA = forceparams[type].u_b.kUBA;
1326         th0B = forceparams[type].u_b.thetaB * DEG2RAD;
1327         kthB = forceparams[type].u_b.kthetaB;
1328         r13B = forceparams[type].u_b.r13B;
1329         kUBB = forceparams[type].u_b.kUBB;
1330
1331         theta = bond_angle(x[ai], x[aj], x[ak], pbc, r_ij, r_kj, &cos_theta, &t1, &t2); /*  41 */
1332
1333         *dvdlambda += harmonic(kthA, kthB, th0A, th0B, theta, lambda, &va, &dVdt); /*  21  */
1334         vtot += va;
1335
1336         ki  = pbc_rvec_sub(pbc, x[ai], x[ak], r_ik); /*   3      */
1337         dr2 = iprod(r_ik, r_ik);                     /*   5             */
1338         dr  = dr2 * gmx::invsqrt(dr2);               /*  10             */
1339
1340         *dvdlambda += harmonic(kUBA, kUBB, r13A, r13B, dr, lambda, &vbond, &fbond); /*  19  */
1341
1342         cos_theta2 = gmx::square(cos_theta); /*   1             */
1343         if (cos_theta2 < 1)
1344         {
1345             real st, sth;
1346             real cik, cii, ckk;
1347             real nrkj2, nrij2;
1348             rvec f_i, f_j, f_k;
1349
1350             st    = dVdt * gmx::invsqrt(1 - cos_theta2); /*  12         */
1351             sth   = st * cos_theta;                      /*   1         */
1352             nrkj2 = iprod(r_kj, r_kj);                   /*   5         */
1353             nrij2 = iprod(r_ij, r_ij);
1354
1355             cik = st * gmx::invsqrt(nrkj2 * nrij2); /*  12              */
1356             cii = sth / nrij2;                      /*  10              */
1357             ckk = sth / nrkj2;                      /*  10              */
1358
1359             for (m = 0; (m < DIM); m++) /*  39          */
1360             {
1361                 f_i[m] = -(cik * r_kj[m] - cii * r_ij[m]);
1362                 f_k[m] = -(cik * r_ij[m] - ckk * r_kj[m]);
1363                 f_j[m] = -f_i[m] - f_k[m];
1364                 f[ai][m] += f_i[m];
1365                 f[aj][m] += f_j[m];
1366                 f[ak][m] += f_k[m];
1367             }
1368             if (computeVirial(flavor))
1369             {
1370                 rvec_inc(fshift[t1], f_i);
1371                 rvec_inc(fshift[CENTRAL], f_j);
1372                 rvec_inc(fshift[t2], f_k);
1373             }
1374         } /* 161 TOTAL  */
1375         /* Time for the bond calculations */
1376         if (dr2 == 0.0)
1377         {
1378             continue;
1379         }
1380
1381         vtot += vbond;              /* 1*/
1382         fbond *= gmx::invsqrt(dr2); /*   6              */
1383
1384         for (m = 0; (m < DIM); m++) /*  15              */
1385         {
1386             fik = fbond * r_ik[m];
1387             f[ai][m] += fik;
1388             f[ak][m] -= fik;
1389             if (computeVirial(flavor))
1390             {
1391                 fshift[ki][m] += fik;
1392                 fshift[CENTRAL][m] -= fik;
1393             }
1394         }
1395     }
1396     return vtot;
1397 }
1398
1399 #if GMX_SIMD_HAVE_REAL
1400
1401 /* As urey_bradley, but using SIMD to calculate many potentials at once.
1402  * This routines does not calculate energies and shift forces.
1403  */
1404 template<BondedKernelFlavor flavor>
1405 std::enable_if_t<flavor == BondedKernelFlavor::ForcesSimdWhenAvailable, real>
1406 urey_bradley(int             nbonds,
1407              const t_iatom   forceatoms[],
1408              const t_iparams forceparams[],
1409              const rvec      x[],
1410              rvec4           f[],
1411              rvec gmx_unused fshift[],
1412              const t_pbc*    pbc,
1413              real gmx_unused lambda,
1414              real gmx_unused* dvdlambda,
1415              const t_mdatoms gmx_unused* md,
1416              t_fcdata gmx_unused* fcd,
1417              int gmx_unused* global_atom_index)
1418 {
1419     constexpr int                            nfa1 = 4;
1420     alignas(GMX_SIMD_ALIGNMENT) std::int32_t ai[GMX_SIMD_REAL_WIDTH];
1421     alignas(GMX_SIMD_ALIGNMENT) std::int32_t aj[GMX_SIMD_REAL_WIDTH];
1422     alignas(GMX_SIMD_ALIGNMENT) std::int32_t ak[GMX_SIMD_REAL_WIDTH];
1423     alignas(GMX_SIMD_ALIGNMENT) real         coeff[4 * GMX_SIMD_REAL_WIDTH];
1424     alignas(GMX_SIMD_ALIGNMENT) real         pbc_simd[9 * GMX_SIMD_REAL_WIDTH];
1425
1426     set_pbc_simd(pbc, pbc_simd);
1427
1428     /* nbonds is the number of angles times nfa1, here we step GMX_SIMD_REAL_WIDTH angles */
1429     for (int i = 0; i < nbonds; i += GMX_SIMD_REAL_WIDTH * nfa1)
1430     {
1431         /* Collect atoms for GMX_SIMD_REAL_WIDTH angles.
1432          * iu indexes into forceatoms, we should not let iu go beyond nbonds.
1433          */
1434         int iu = i;
1435         for (int s = 0; s < GMX_SIMD_REAL_WIDTH; s++)
1436         {
1437             const int type = forceatoms[iu];
1438             ai[s]          = forceatoms[iu + 1];
1439             aj[s]          = forceatoms[iu + 2];
1440             ak[s]          = forceatoms[iu + 3];
1441
1442             /* At the end fill the arrays with the last atoms and 0 params */
1443             if (i + s * nfa1 < nbonds)
1444             {
1445                 coeff[s]                           = forceparams[type].u_b.kthetaA;
1446                 coeff[GMX_SIMD_REAL_WIDTH + s]     = forceparams[type].u_b.thetaA;
1447                 coeff[GMX_SIMD_REAL_WIDTH * 2 + s] = forceparams[type].u_b.kUBA;
1448                 coeff[GMX_SIMD_REAL_WIDTH * 3 + s] = forceparams[type].u_b.r13A;
1449
1450                 if (iu + nfa1 < nbonds)
1451                 {
1452                     iu += nfa1;
1453                 }
1454             }
1455             else
1456             {
1457                 coeff[s]                           = 0;
1458                 coeff[GMX_SIMD_REAL_WIDTH + s]     = 0;
1459                 coeff[GMX_SIMD_REAL_WIDTH * 2 + s] = 0;
1460                 coeff[GMX_SIMD_REAL_WIDTH * 3 + s] = 0;
1461             }
1462         }
1463
1464         SimdReal xi_S, yi_S, zi_S;
1465         SimdReal xj_S, yj_S, zj_S;
1466         SimdReal xk_S, yk_S, zk_S;
1467
1468         /* Store the non PBC corrected distances packed and aligned */
1469         gatherLoadUTranspose<3>(reinterpret_cast<const real*>(x), ai, &xi_S, &yi_S, &zi_S);
1470         gatherLoadUTranspose<3>(reinterpret_cast<const real*>(x), aj, &xj_S, &yj_S, &zj_S);
1471         gatherLoadUTranspose<3>(reinterpret_cast<const real*>(x), ak, &xk_S, &yk_S, &zk_S);
1472         SimdReal rijx_S = xi_S - xj_S;
1473         SimdReal rijy_S = yi_S - yj_S;
1474         SimdReal rijz_S = zi_S - zj_S;
1475         SimdReal rkjx_S = xk_S - xj_S;
1476         SimdReal rkjy_S = yk_S - yj_S;
1477         SimdReal rkjz_S = zk_S - zj_S;
1478         SimdReal rikx_S = xi_S - xk_S;
1479         SimdReal riky_S = yi_S - yk_S;
1480         SimdReal rikz_S = zi_S - zk_S;
1481
1482         const SimdReal ktheta_S = load<SimdReal>(coeff);
1483         const SimdReal theta0_S = load<SimdReal>(coeff + GMX_SIMD_REAL_WIDTH) * DEG2RAD;
1484         const SimdReal kUB_S    = load<SimdReal>(coeff + 2 * GMX_SIMD_REAL_WIDTH);
1485         const SimdReal r13_S    = load<SimdReal>(coeff + 3 * GMX_SIMD_REAL_WIDTH);
1486
1487         pbc_correct_dx_simd(&rijx_S, &rijy_S, &rijz_S, pbc_simd);
1488         pbc_correct_dx_simd(&rkjx_S, &rkjy_S, &rkjz_S, pbc_simd);
1489         pbc_correct_dx_simd(&rikx_S, &riky_S, &rikz_S, pbc_simd);
1490
1491         const SimdReal rij_rkj_S = iprod(rijx_S, rijy_S, rijz_S, rkjx_S, rkjy_S, rkjz_S);
1492
1493         const SimdReal dr2_S = iprod(rikx_S, riky_S, rikz_S, rikx_S, riky_S, rikz_S);
1494
1495         const SimdReal nrij2_S = norm2(rijx_S, rijy_S, rijz_S);
1496         const SimdReal nrkj2_S = norm2(rkjx_S, rkjy_S, rkjz_S);
1497
1498         const SimdReal nrij_1_S = invsqrt(nrij2_S);
1499         const SimdReal nrkj_1_S = invsqrt(nrkj2_S);
1500         const SimdReal invdr2_S = invsqrt(dr2_S);
1501         const SimdReal dr_S     = dr2_S * invdr2_S;
1502
1503         constexpr real min_one_plus_eps = -1.0 + 2.0 * GMX_REAL_EPS; // Smallest number > -1
1504
1505         /* To allow for 180 degrees, we take the max of cos and -1 + 1bit,
1506          * so we can safely get the 1/sin from 1/sqrt(1 - cos^2).
1507          * This also ensures that rounding errors would cause the argument
1508          * of simdAcos to be < -1.
1509          * Note that we do not take precautions for cos(0)=1, so the bonds
1510          * in an angle should not align at an angle of 0 degrees.
1511          */
1512         const SimdReal cos_S = max(rij_rkj_S * nrij_1_S * nrkj_1_S, min_one_plus_eps);
1513
1514         const SimdReal theta_S  = acos(cos_S);
1515         const SimdReal invsin_S = invsqrt(1.0 - cos_S * cos_S);
1516         const SimdReal st_S     = ktheta_S * (theta0_S - theta_S) * invsin_S;
1517         const SimdReal sth_S    = st_S * cos_S;
1518
1519         const SimdReal cik_S = st_S * nrij_1_S * nrkj_1_S;
1520         const SimdReal cii_S = sth_S * nrij_1_S * nrij_1_S;
1521         const SimdReal ckk_S = sth_S * nrkj_1_S * nrkj_1_S;
1522
1523         const SimdReal sUB_S = kUB_S * (r13_S - dr_S) * invdr2_S;
1524
1525         const SimdReal f_ikx_S = sUB_S * rikx_S;
1526         const SimdReal f_iky_S = sUB_S * riky_S;
1527         const SimdReal f_ikz_S = sUB_S * rikz_S;
1528
1529         const SimdReal f_ix_S = fnma(cik_S, rkjx_S, cii_S * rijx_S) + f_ikx_S;
1530         const SimdReal f_iy_S = fnma(cik_S, rkjy_S, cii_S * rijy_S) + f_iky_S;
1531         const SimdReal f_iz_S = fnma(cik_S, rkjz_S, cii_S * rijz_S) + f_ikz_S;
1532         const SimdReal f_kx_S = fnma(cik_S, rijx_S, ckk_S * rkjx_S) - f_ikx_S;
1533         const SimdReal f_ky_S = fnma(cik_S, rijy_S, ckk_S * rkjy_S) - f_iky_S;
1534         const SimdReal f_kz_S = fnma(cik_S, rijz_S, ckk_S * rkjz_S) - f_ikz_S;
1535
1536         transposeScatterIncrU<4>(reinterpret_cast<real*>(f), ai, f_ix_S, f_iy_S, f_iz_S);
1537         transposeScatterDecrU<4>(reinterpret_cast<real*>(f), aj, f_ix_S + f_kx_S, f_iy_S + f_ky_S,
1538                                  f_iz_S + f_kz_S);
1539         transposeScatterIncrU<4>(reinterpret_cast<real*>(f), ak, f_kx_S, f_ky_S, f_kz_S);
1540     }
1541
1542     return 0;
1543 }
1544
1545 #endif // GMX_SIMD_HAVE_REAL
1546
1547 template<BondedKernelFlavor flavor>
1548 real quartic_angles(int             nbonds,
1549                     const t_iatom   forceatoms[],
1550                     const t_iparams forceparams[],
1551                     const rvec      x[],
1552                     rvec4           f[],
1553                     rvec            fshift[],
1554                     const t_pbc*    pbc,
1555                     real gmx_unused lambda,
1556                     real gmx_unused* dvdlambda,
1557                     const t_mdatoms gmx_unused* md,
1558                     t_fcdata gmx_unused* fcd,
1559                     int gmx_unused* global_atom_index)
1560 {
1561     int  i, j, ai, aj, ak, t1, t2, type;
1562     rvec r_ij, r_kj;
1563     real cos_theta, cos_theta2, theta, dt, dVdt, va, dtp, c, vtot;
1564
1565     vtot = 0.0;
1566     for (i = 0; (i < nbonds);)
1567     {
1568         type = forceatoms[i++];
1569         ai   = forceatoms[i++];
1570         aj   = forceatoms[i++];
1571         ak   = forceatoms[i++];
1572
1573         theta = bond_angle(x[ai], x[aj], x[ak], pbc, r_ij, r_kj, &cos_theta, &t1, &t2); /*  41 */
1574
1575         dt = theta - forceparams[type].qangle.theta * DEG2RAD; /* 2          */
1576
1577         dVdt = 0;
1578         va   = forceparams[type].qangle.c[0];
1579         dtp  = 1.0;
1580         for (j = 1; j <= 4; j++)
1581         {
1582             c = forceparams[type].qangle.c[j];
1583             dVdt -= j * c * dtp;
1584             dtp *= dt;
1585             va += c * dtp;
1586         }
1587         /* 20 */
1588
1589         vtot += va;
1590
1591         cos_theta2 = gmx::square(cos_theta); /*   1             */
1592         if (cos_theta2 < 1)
1593         {
1594             int  m;
1595             real st, sth;
1596             real cik, cii, ckk;
1597             real nrkj2, nrij2;
1598             rvec f_i, f_j, f_k;
1599
1600             st    = dVdt * gmx::invsqrt(1 - cos_theta2); /*  12         */
1601             sth   = st * cos_theta;                      /*   1         */
1602             nrkj2 = iprod(r_kj, r_kj);                   /*   5         */
1603             nrij2 = iprod(r_ij, r_ij);
1604
1605             cik = st * gmx::invsqrt(nrkj2 * nrij2); /*  12              */
1606             cii = sth / nrij2;                      /*  10              */
1607             ckk = sth / nrkj2;                      /*  10              */
1608
1609             for (m = 0; (m < DIM); m++) /*  39          */
1610             {
1611                 f_i[m] = -(cik * r_kj[m] - cii * r_ij[m]);
1612                 f_k[m] = -(cik * r_ij[m] - ckk * r_kj[m]);
1613                 f_j[m] = -f_i[m] - f_k[m];
1614                 f[ai][m] += f_i[m];
1615                 f[aj][m] += f_j[m];
1616                 f[ak][m] += f_k[m];
1617             }
1618
1619             if (computeVirial(flavor))
1620             {
1621                 rvec_inc(fshift[t1], f_i);
1622                 rvec_inc(fshift[CENTRAL], f_j);
1623                 rvec_inc(fshift[t2], f_k);
1624             }
1625         } /* 153 TOTAL  */
1626     }
1627     return vtot;
1628 }
1629
1630
1631 #if GMX_SIMD_HAVE_REAL
1632
1633 /* As dih_angle above, but calculates 4 dihedral angles at once using SIMD,
1634  * also calculates the pre-factor required for the dihedral force update.
1635  * Note that bv and buf should be register aligned.
1636  */
1637 inline void dih_angle_simd(const rvec* x,
1638                            const int*  ai,
1639                            const int*  aj,
1640                            const int*  ak,
1641                            const int*  al,
1642                            const real* pbc_simd,
1643                            SimdReal*   phi_S,
1644                            SimdReal*   mx_S,
1645                            SimdReal*   my_S,
1646                            SimdReal*   mz_S,
1647                            SimdReal*   nx_S,
1648                            SimdReal*   ny_S,
1649                            SimdReal*   nz_S,
1650                            SimdReal*   nrkj_m2_S,
1651                            SimdReal*   nrkj_n2_S,
1652                            SimdReal*   p_S,
1653                            SimdReal*   q_S)
1654 {
1655     SimdReal xi_S, yi_S, zi_S;
1656     SimdReal xj_S, yj_S, zj_S;
1657     SimdReal xk_S, yk_S, zk_S;
1658     SimdReal xl_S, yl_S, zl_S;
1659     SimdReal rijx_S, rijy_S, rijz_S;
1660     SimdReal rkjx_S, rkjy_S, rkjz_S;
1661     SimdReal rklx_S, rkly_S, rklz_S;
1662     SimdReal cx_S, cy_S, cz_S;
1663     SimdReal cn_S;
1664     SimdReal s_S;
1665     SimdReal ipr_S;
1666     SimdReal iprm_S, iprn_S;
1667     SimdReal nrkj2_S, nrkj_1_S, nrkj_2_S, nrkj_S;
1668     SimdReal toler_S;
1669     SimdReal nrkj2_min_S;
1670     SimdReal real_eps_S;
1671
1672     /* Used to avoid division by zero.
1673      * We take into acount that we multiply the result by real_eps_S.
1674      */
1675     nrkj2_min_S = SimdReal(GMX_REAL_MIN / (2 * GMX_REAL_EPS));
1676
1677     /* The value of the last significant bit (GMX_REAL_EPS is half of that) */
1678     real_eps_S = SimdReal(2 * GMX_REAL_EPS);
1679
1680     /* Store the non PBC corrected distances packed and aligned */
1681     gatherLoadUTranspose<3>(reinterpret_cast<const real*>(x), ai, &xi_S, &yi_S, &zi_S);
1682     gatherLoadUTranspose<3>(reinterpret_cast<const real*>(x), aj, &xj_S, &yj_S, &zj_S);
1683     gatherLoadUTranspose<3>(reinterpret_cast<const real*>(x), ak, &xk_S, &yk_S, &zk_S);
1684     gatherLoadUTranspose<3>(reinterpret_cast<const real*>(x), al, &xl_S, &yl_S, &zl_S);
1685     rijx_S = xi_S - xj_S;
1686     rijy_S = yi_S - yj_S;
1687     rijz_S = zi_S - zj_S;
1688     rkjx_S = xk_S - xj_S;
1689     rkjy_S = yk_S - yj_S;
1690     rkjz_S = zk_S - zj_S;
1691     rklx_S = xk_S - xl_S;
1692     rkly_S = yk_S - yl_S;
1693     rklz_S = zk_S - zl_S;
1694
1695     pbc_correct_dx_simd(&rijx_S, &rijy_S, &rijz_S, pbc_simd);
1696     pbc_correct_dx_simd(&rkjx_S, &rkjy_S, &rkjz_S, pbc_simd);
1697     pbc_correct_dx_simd(&rklx_S, &rkly_S, &rklz_S, pbc_simd);
1698
1699     cprod(rijx_S, rijy_S, rijz_S, rkjx_S, rkjy_S, rkjz_S, mx_S, my_S, mz_S);
1700
1701     cprod(rkjx_S, rkjy_S, rkjz_S, rklx_S, rkly_S, rklz_S, nx_S, ny_S, nz_S);
1702
1703     cprod(*mx_S, *my_S, *mz_S, *nx_S, *ny_S, *nz_S, &cx_S, &cy_S, &cz_S);
1704
1705     cn_S = sqrt(norm2(cx_S, cy_S, cz_S));
1706
1707     s_S = iprod(*mx_S, *my_S, *mz_S, *nx_S, *ny_S, *nz_S);
1708
1709     /* Determine the dihedral angle, the sign might need correction */
1710     *phi_S = atan2(cn_S, s_S);
1711
1712     ipr_S = iprod(rijx_S, rijy_S, rijz_S, *nx_S, *ny_S, *nz_S);
1713
1714     iprm_S = norm2(*mx_S, *my_S, *mz_S);
1715     iprn_S = norm2(*nx_S, *ny_S, *nz_S);
1716
1717     nrkj2_S = norm2(rkjx_S, rkjy_S, rkjz_S);
1718
1719     /* Avoid division by zero. When zero, the result is multiplied by 0
1720      * anyhow, so the 3 max below do not affect the final result.
1721      */
1722     nrkj2_S  = max(nrkj2_S, nrkj2_min_S);
1723     nrkj_1_S = invsqrt(nrkj2_S);
1724     nrkj_2_S = nrkj_1_S * nrkj_1_S;
1725     nrkj_S   = nrkj2_S * nrkj_1_S;
1726
1727     toler_S = nrkj2_S * real_eps_S;
1728
1729     /* Here the plain-C code uses a conditional, but we can't do that in SIMD.
1730      * So we take a max with the tolerance instead. Since we multiply with
1731      * m or n later, the max does not affect the results.
1732      */
1733     iprm_S     = max(iprm_S, toler_S);
1734     iprn_S     = max(iprn_S, toler_S);
1735     *nrkj_m2_S = nrkj_S * inv(iprm_S);
1736     *nrkj_n2_S = nrkj_S * inv(iprn_S);
1737
1738     /* Set sign of phi_S with the sign of ipr_S; phi_S is currently positive */
1739     *phi_S = copysign(*phi_S, ipr_S);
1740     *p_S   = iprod(rijx_S, rijy_S, rijz_S, rkjx_S, rkjy_S, rkjz_S);
1741     *p_S   = *p_S * nrkj_2_S;
1742
1743     *q_S = iprod(rklx_S, rkly_S, rklz_S, rkjx_S, rkjy_S, rkjz_S);
1744     *q_S = *q_S * nrkj_2_S;
1745 }
1746
1747 #endif // GMX_SIMD_HAVE_REAL
1748
1749 } // namespace
1750
1751 template<BondedKernelFlavor flavor>
1752 void do_dih_fup(int          i,
1753                 int          j,
1754                 int          k,
1755                 int          l,
1756                 real         ddphi,
1757                 rvec         r_ij,
1758                 rvec         r_kj,
1759                 rvec         r_kl,
1760                 rvec         m,
1761                 rvec         n,
1762                 rvec4        f[],
1763                 rvec         fshift[],
1764                 const t_pbc* pbc,
1765                 const rvec   x[],
1766                 int          t1,
1767                 int          t2,
1768                 int          t3)
1769 {
1770     /* 143 FLOPS */
1771     rvec f_i, f_j, f_k, f_l;
1772     rvec uvec, vvec, svec, dx_jl;
1773     real iprm, iprn, nrkj, nrkj2, nrkj_1, nrkj_2;
1774     real a, b, p, q, toler;
1775
1776     iprm  = iprod(m, m);       /*  5    */
1777     iprn  = iprod(n, n);       /*  5    */
1778     nrkj2 = iprod(r_kj, r_kj); /*  5    */
1779     toler = nrkj2 * GMX_REAL_EPS;
1780     if ((iprm > toler) && (iprn > toler))
1781     {
1782         nrkj_1 = gmx::invsqrt(nrkj2);  /* 10    */
1783         nrkj_2 = nrkj_1 * nrkj_1;      /*  1    */
1784         nrkj   = nrkj2 * nrkj_1;       /*  1    */
1785         a      = -ddphi * nrkj / iprm; /* 11    */
1786         svmul(a, m, f_i);              /*  3    */
1787         b = ddphi * nrkj / iprn;       /* 11    */
1788         svmul(b, n, f_l);              /*  3  */
1789         p = iprod(r_ij, r_kj);         /*  5    */
1790         p *= nrkj_2;                   /*  1    */
1791         q = iprod(r_kl, r_kj);         /*  5    */
1792         q *= nrkj_2;                   /*  1    */
1793         svmul(p, f_i, uvec);           /*  3    */
1794         svmul(q, f_l, vvec);           /*  3    */
1795         rvec_sub(uvec, vvec, svec);    /*  3    */
1796         rvec_sub(f_i, svec, f_j);      /*  3    */
1797         rvec_add(f_l, svec, f_k);      /*  3    */
1798         rvec_inc(f[i], f_i);           /*  3    */
1799         rvec_dec(f[j], f_j);           /*  3    */
1800         rvec_dec(f[k], f_k);           /*  3    */
1801         rvec_inc(f[l], f_l);           /*  3    */
1802
1803         if (computeVirial(flavor))
1804         {
1805             if (pbc)
1806             {
1807                 t3 = pbc_rvec_sub(pbc, x[l], x[j], dx_jl);
1808             }
1809             else
1810             {
1811                 t3 = CENTRAL;
1812             }
1813
1814             rvec_inc(fshift[t1], f_i);
1815             rvec_dec(fshift[CENTRAL], f_j);
1816             rvec_dec(fshift[t2], f_k);
1817             rvec_inc(fshift[t3], f_l);
1818         }
1819     }
1820     /* 112 TOTAL    */
1821 }
1822
1823 namespace
1824 {
1825
1826 #if GMX_SIMD_HAVE_REAL
1827 /* As do_dih_fup_noshiftf above, but with SIMD and pre-calculated pre-factors */
1828 inline void gmx_simdcall do_dih_fup_noshiftf_simd(const int* ai,
1829                                                   const int* aj,
1830                                                   const int* ak,
1831                                                   const int* al,
1832                                                   SimdReal   p,
1833                                                   SimdReal   q,
1834                                                   SimdReal   f_i_x,
1835                                                   SimdReal   f_i_y,
1836                                                   SimdReal   f_i_z,
1837                                                   SimdReal   mf_l_x,
1838                                                   SimdReal   mf_l_y,
1839                                                   SimdReal   mf_l_z,
1840                                                   rvec4      f[])
1841 {
1842     SimdReal sx    = p * f_i_x + q * mf_l_x;
1843     SimdReal sy    = p * f_i_y + q * mf_l_y;
1844     SimdReal sz    = p * f_i_z + q * mf_l_z;
1845     SimdReal f_j_x = f_i_x - sx;
1846     SimdReal f_j_y = f_i_y - sy;
1847     SimdReal f_j_z = f_i_z - sz;
1848     SimdReal f_k_x = mf_l_x - sx;
1849     SimdReal f_k_y = mf_l_y - sy;
1850     SimdReal f_k_z = mf_l_z - sz;
1851     transposeScatterIncrU<4>(reinterpret_cast<real*>(f), ai, f_i_x, f_i_y, f_i_z);
1852     transposeScatterDecrU<4>(reinterpret_cast<real*>(f), aj, f_j_x, f_j_y, f_j_z);
1853     transposeScatterIncrU<4>(reinterpret_cast<real*>(f), ak, f_k_x, f_k_y, f_k_z);
1854     transposeScatterDecrU<4>(reinterpret_cast<real*>(f), al, mf_l_x, mf_l_y, mf_l_z);
1855 }
1856 #endif // GMX_SIMD_HAVE_REAL
1857
1858 /*! \brief Computes and returns the proper dihedral force
1859  *
1860  * With the appropriate kernel flavor, also computes and accumulates
1861  * the energy and dV/dlambda.
1862  */
1863 template<BondedKernelFlavor flavor>
1864 real dopdihs(real cpA, real cpB, real phiA, real phiB, int mult, real phi, real lambda, real* V, real* dvdlambda)
1865 {
1866     const real L1   = 1.0 - lambda;
1867     const real ph0  = (L1 * phiA + lambda * phiB) * DEG2RAD;
1868     const real dph0 = (phiB - phiA) * DEG2RAD;
1869     const real cp   = L1 * cpA + lambda * cpB;
1870
1871     const real mdphi = mult * phi - ph0;
1872     const real sdphi = std::sin(mdphi);
1873     const real ddphi = -cp * mult * sdphi;
1874     if (computeEnergy(flavor))
1875     {
1876         const real v1 = 1 + std::cos(mdphi);
1877         *V += cp * v1;
1878
1879         *dvdlambda += (cpB - cpA) * v1 + cp * dph0 * sdphi;
1880     }
1881
1882     return ddphi;
1883     /* That was 40 flops */
1884 }
1885
1886 /*! \brief Similar to \p dopdihs(), except for a minus sign and a different treatment of mult/phi0 */
1887 real dopdihs_min(real cpA, real cpB, real phiA, real phiB, int mult, real phi, real lambda, real* V, real* F)
1888 {
1889     real v, dvdlambda, mdphi, v1, sdphi, ddphi;
1890     real L1   = 1.0 - lambda;
1891     real ph0  = (L1 * phiA + lambda * phiB) * DEG2RAD;
1892     real dph0 = (phiB - phiA) * DEG2RAD;
1893     real cp   = L1 * cpA + lambda * cpB;
1894
1895     mdphi = mult * (phi - ph0);
1896     sdphi = std::sin(mdphi);
1897     ddphi = cp * mult * sdphi;
1898     v1    = 1.0 - std::cos(mdphi);
1899     v     = cp * v1;
1900
1901     dvdlambda = (cpB - cpA) * v1 + cp * dph0 * sdphi;
1902
1903     *V = v;
1904     *F = ddphi;
1905
1906     return dvdlambda;
1907
1908     /* That was 40 flops */
1909 }
1910
1911 template<BondedKernelFlavor flavor>
1912 std::enable_if_t<flavor != BondedKernelFlavor::ForcesSimdWhenAvailable || !GMX_SIMD_HAVE_REAL, real>
1913 pdihs(int             nbonds,
1914       const t_iatom   forceatoms[],
1915       const t_iparams forceparams[],
1916       const rvec      x[],
1917       rvec4           f[],
1918       rvec            fshift[],
1919       const t_pbc*    pbc,
1920       real            lambda,
1921       real*           dvdlambda,
1922       const t_mdatoms gmx_unused* md,
1923       t_fcdata gmx_unused* fcd,
1924       int gmx_unused* global_atom_index)
1925 {
1926     int  t1, t2, t3;
1927     rvec r_ij, r_kj, r_kl, m, n;
1928
1929     real vtot = 0.0;
1930
1931     for (int i = 0; i < nbonds;)
1932     {
1933         const int ai = forceatoms[i + 1];
1934         const int aj = forceatoms[i + 2];
1935         const int ak = forceatoms[i + 3];
1936         const int al = forceatoms[i + 4];
1937
1938         const real phi = dih_angle(x[ai], x[aj], x[ak], x[al], pbc, r_ij, r_kj, r_kl, m, n, &t1,
1939                                    &t2, &t3); /*  84      */
1940
1941         /* Loop over dihedrals working on the same atoms,
1942          * so we avoid recalculating angles and distributing forces.
1943          */
1944         real ddphi_tot = 0;
1945         do
1946         {
1947             const int type = forceatoms[i];
1948             ddphi_tot += dopdihs<flavor>(forceparams[type].pdihs.cpA, forceparams[type].pdihs.cpB,
1949                                          forceparams[type].pdihs.phiA, forceparams[type].pdihs.phiB,
1950                                          forceparams[type].pdihs.mult, phi, lambda, &vtot, dvdlambda);
1951
1952             i += 5;
1953         } while (i < nbonds && forceatoms[i + 1] == ai && forceatoms[i + 2] == aj
1954                  && forceatoms[i + 3] == ak && forceatoms[i + 4] == al);
1955
1956         do_dih_fup<flavor>(ai, aj, ak, al, ddphi_tot, r_ij, r_kj, r_kl, m, n, f, fshift, pbc, x, t1,
1957                            t2, t3); /* 112              */
1958     }                               /* 223 TOTAL  */
1959
1960     return vtot;
1961 }
1962
1963 #if GMX_SIMD_HAVE_REAL
1964
1965 /* As pdihs above, but using SIMD to calculate multiple dihedrals at once */
1966 template<BondedKernelFlavor flavor>
1967 std::enable_if_t<flavor == BondedKernelFlavor::ForcesSimdWhenAvailable, real>
1968 pdihs(int             nbonds,
1969       const t_iatom   forceatoms[],
1970       const t_iparams forceparams[],
1971       const rvec      x[],
1972       rvec4           f[],
1973       rvec gmx_unused fshift[],
1974       const t_pbc*    pbc,
1975       real gmx_unused lambda,
1976       real gmx_unused* dvdlambda,
1977       const t_mdatoms gmx_unused* md,
1978       t_fcdata gmx_unused* fcd,
1979       int gmx_unused* global_atom_index)
1980 {
1981     const int                                nfa1 = 5;
1982     int                                      i, iu, s;
1983     int                                      type;
1984     alignas(GMX_SIMD_ALIGNMENT) std::int32_t ai[GMX_SIMD_REAL_WIDTH];
1985     alignas(GMX_SIMD_ALIGNMENT) std::int32_t aj[GMX_SIMD_REAL_WIDTH];
1986     alignas(GMX_SIMD_ALIGNMENT) std::int32_t ak[GMX_SIMD_REAL_WIDTH];
1987     alignas(GMX_SIMD_ALIGNMENT) std::int32_t al[GMX_SIMD_REAL_WIDTH];
1988     alignas(GMX_SIMD_ALIGNMENT) real         buf[3 * GMX_SIMD_REAL_WIDTH];
1989     real *                                   cp, *phi0, *mult;
1990     SimdReal                                 deg2rad_S(DEG2RAD);
1991     SimdReal                                 p_S, q_S;
1992     SimdReal                                 phi0_S, phi_S;
1993     SimdReal                                 mx_S, my_S, mz_S;
1994     SimdReal                                 nx_S, ny_S, nz_S;
1995     SimdReal                                 nrkj_m2_S, nrkj_n2_S;
1996     SimdReal                                 cp_S, mdphi_S, mult_S;
1997     SimdReal                                 sin_S, cos_S;
1998     SimdReal                                 mddphi_S;
1999     SimdReal                                 sf_i_S, msf_l_S;
2000     alignas(GMX_SIMD_ALIGNMENT) real         pbc_simd[9 * GMX_SIMD_REAL_WIDTH];
2001
2002     /* Extract aligned pointer for parameters and variables */
2003     cp   = buf + 0 * GMX_SIMD_REAL_WIDTH;
2004     phi0 = buf + 1 * GMX_SIMD_REAL_WIDTH;
2005     mult = buf + 2 * GMX_SIMD_REAL_WIDTH;
2006
2007     set_pbc_simd(pbc, pbc_simd);
2008
2009     /* nbonds is the number of dihedrals times nfa1, here we step GMX_SIMD_REAL_WIDTH dihs */
2010     for (i = 0; (i < nbonds); i += GMX_SIMD_REAL_WIDTH * nfa1)
2011     {
2012         /* Collect atoms quadruplets for GMX_SIMD_REAL_WIDTH dihedrals.
2013          * iu indexes into forceatoms, we should not let iu go beyond nbonds.
2014          */
2015         iu = i;
2016         for (s = 0; s < GMX_SIMD_REAL_WIDTH; s++)
2017         {
2018             type  = forceatoms[iu];
2019             ai[s] = forceatoms[iu + 1];
2020             aj[s] = forceatoms[iu + 2];
2021             ak[s] = forceatoms[iu + 3];
2022             al[s] = forceatoms[iu + 4];
2023
2024             /* At the end fill the arrays with the last atoms and 0 params */
2025             if (i + s * nfa1 < nbonds)
2026             {
2027                 cp[s]   = forceparams[type].pdihs.cpA;
2028                 phi0[s] = forceparams[type].pdihs.phiA;
2029                 mult[s] = forceparams[type].pdihs.mult;
2030
2031                 if (iu + nfa1 < nbonds)
2032                 {
2033                     iu += nfa1;
2034                 }
2035             }
2036             else
2037             {
2038                 cp[s]   = 0;
2039                 phi0[s] = 0;
2040                 mult[s] = 0;
2041             }
2042         }
2043
2044         /* Calculate GMX_SIMD_REAL_WIDTH dihedral angles at once */
2045         dih_angle_simd(x, ai, aj, ak, al, pbc_simd, &phi_S, &mx_S, &my_S, &mz_S, &nx_S, &ny_S,
2046                        &nz_S, &nrkj_m2_S, &nrkj_n2_S, &p_S, &q_S);
2047
2048         cp_S   = load<SimdReal>(cp);
2049         phi0_S = load<SimdReal>(phi0) * deg2rad_S;
2050         mult_S = load<SimdReal>(mult);
2051
2052         mdphi_S = fms(mult_S, phi_S, phi0_S);
2053
2054         /* Calculate GMX_SIMD_REAL_WIDTH sines at once */
2055         sincos(mdphi_S, &sin_S, &cos_S);
2056         mddphi_S = cp_S * mult_S * sin_S;
2057         sf_i_S   = mddphi_S * nrkj_m2_S;
2058         msf_l_S  = mddphi_S * nrkj_n2_S;
2059
2060         /* After this m?_S will contain f[i] */
2061         mx_S = sf_i_S * mx_S;
2062         my_S = sf_i_S * my_S;
2063         mz_S = sf_i_S * mz_S;
2064
2065         /* After this m?_S will contain -f[l] */
2066         nx_S = msf_l_S * nx_S;
2067         ny_S = msf_l_S * ny_S;
2068         nz_S = msf_l_S * nz_S;
2069
2070         do_dih_fup_noshiftf_simd(ai, aj, ak, al, p_S, q_S, mx_S, my_S, mz_S, nx_S, ny_S, nz_S, f);
2071     }
2072
2073     return 0;
2074 }
2075
2076 /* This is mostly a copy of the SIMD flavor of pdihs above, but with using
2077  * the RB potential instead of a harmonic potential.
2078  * This function can replace rbdihs() when no energy and virial are needed.
2079  */
2080 template<BondedKernelFlavor flavor>
2081 std::enable_if_t<flavor == BondedKernelFlavor::ForcesSimdWhenAvailable, real>
2082 rbdihs(int             nbonds,
2083        const t_iatom   forceatoms[],
2084        const t_iparams forceparams[],
2085        const rvec      x[],
2086        rvec4           f[],
2087        rvec gmx_unused fshift[],
2088        const t_pbc*    pbc,
2089        real gmx_unused lambda,
2090        real gmx_unused* dvdlambda,
2091        const t_mdatoms gmx_unused* md,
2092        t_fcdata gmx_unused* fcd,
2093        int gmx_unused* global_atom_index)
2094 {
2095     const int                                nfa1 = 5;
2096     int                                      i, iu, s, j;
2097     int                                      type;
2098     alignas(GMX_SIMD_ALIGNMENT) std::int32_t ai[GMX_SIMD_REAL_WIDTH];
2099     alignas(GMX_SIMD_ALIGNMENT) std::int32_t aj[GMX_SIMD_REAL_WIDTH];
2100     alignas(GMX_SIMD_ALIGNMENT) std::int32_t ak[GMX_SIMD_REAL_WIDTH];
2101     alignas(GMX_SIMD_ALIGNMENT) std::int32_t al[GMX_SIMD_REAL_WIDTH];
2102     alignas(GMX_SIMD_ALIGNMENT) real         parm[NR_RBDIHS * GMX_SIMD_REAL_WIDTH];
2103
2104     SimdReal                         p_S, q_S;
2105     SimdReal                         phi_S;
2106     SimdReal                         ddphi_S, cosfac_S;
2107     SimdReal                         mx_S, my_S, mz_S;
2108     SimdReal                         nx_S, ny_S, nz_S;
2109     SimdReal                         nrkj_m2_S, nrkj_n2_S;
2110     SimdReal                         parm_S, c_S;
2111     SimdReal                         sin_S, cos_S;
2112     SimdReal                         sf_i_S, msf_l_S;
2113     alignas(GMX_SIMD_ALIGNMENT) real pbc_simd[9 * GMX_SIMD_REAL_WIDTH];
2114
2115     SimdReal pi_S(M_PI);
2116     SimdReal one_S(1.0);
2117
2118     set_pbc_simd(pbc, pbc_simd);
2119
2120     /* nbonds is the number of dihedrals times nfa1, here we step GMX_SIMD_REAL_WIDTH dihs */
2121     for (i = 0; (i < nbonds); i += GMX_SIMD_REAL_WIDTH * nfa1)
2122     {
2123         /* Collect atoms quadruplets for GMX_SIMD_REAL_WIDTH dihedrals.
2124          * iu indexes into forceatoms, we should not let iu go beyond nbonds.
2125          */
2126         iu = i;
2127         for (s = 0; s < GMX_SIMD_REAL_WIDTH; s++)
2128         {
2129             type  = forceatoms[iu];
2130             ai[s] = forceatoms[iu + 1];
2131             aj[s] = forceatoms[iu + 2];
2132             ak[s] = forceatoms[iu + 3];
2133             al[s] = forceatoms[iu + 4];
2134
2135             /* At the end fill the arrays with the last atoms and 0 params */
2136             if (i + s * nfa1 < nbonds)
2137             {
2138                 /* We don't need the first parameter, since that's a constant
2139                  * which only affects the energies, not the forces.
2140                  */
2141                 for (j = 1; j < NR_RBDIHS; j++)
2142                 {
2143                     parm[j * GMX_SIMD_REAL_WIDTH + s] = forceparams[type].rbdihs.rbcA[j];
2144                 }
2145
2146                 if (iu + nfa1 < nbonds)
2147                 {
2148                     iu += nfa1;
2149                 }
2150             }
2151             else
2152             {
2153                 for (j = 1; j < NR_RBDIHS; j++)
2154                 {
2155                     parm[j * GMX_SIMD_REAL_WIDTH + s] = 0;
2156                 }
2157             }
2158         }
2159
2160         /* Calculate GMX_SIMD_REAL_WIDTH dihedral angles at once */
2161         dih_angle_simd(x, ai, aj, ak, al, pbc_simd, &phi_S, &mx_S, &my_S, &mz_S, &nx_S, &ny_S,
2162                        &nz_S, &nrkj_m2_S, &nrkj_n2_S, &p_S, &q_S);
2163
2164         /* Change to polymer convention */
2165         phi_S = phi_S - pi_S;
2166
2167         sincos(phi_S, &sin_S, &cos_S);
2168
2169         ddphi_S  = setZero();
2170         c_S      = one_S;
2171         cosfac_S = one_S;
2172         for (j = 1; j < NR_RBDIHS; j++)
2173         {
2174             parm_S   = load<SimdReal>(parm + j * GMX_SIMD_REAL_WIDTH);
2175             ddphi_S  = fma(c_S * parm_S, cosfac_S, ddphi_S);
2176             cosfac_S = cosfac_S * cos_S;
2177             c_S      = c_S + one_S;
2178         }
2179
2180         /* Note that here we do not use the minus sign which is present
2181          * in the normal RB code. This is corrected for through (m)sf below.
2182          */
2183         ddphi_S = ddphi_S * sin_S;
2184
2185         sf_i_S  = ddphi_S * nrkj_m2_S;
2186         msf_l_S = ddphi_S * nrkj_n2_S;
2187
2188         /* After this m?_S will contain f[i] */
2189         mx_S = sf_i_S * mx_S;
2190         my_S = sf_i_S * my_S;
2191         mz_S = sf_i_S * mz_S;
2192
2193         /* After this m?_S will contain -f[l] */
2194         nx_S = msf_l_S * nx_S;
2195         ny_S = msf_l_S * ny_S;
2196         nz_S = msf_l_S * nz_S;
2197
2198         do_dih_fup_noshiftf_simd(ai, aj, ak, al, p_S, q_S, mx_S, my_S, mz_S, nx_S, ny_S, nz_S, f);
2199     }
2200
2201     return 0;
2202 }
2203
2204 #endif // GMX_SIMD_HAVE_REAL
2205
2206
2207 template<BondedKernelFlavor flavor>
2208 real idihs(int             nbonds,
2209            const t_iatom   forceatoms[],
2210            const t_iparams forceparams[],
2211            const rvec      x[],
2212            rvec4           f[],
2213            rvec            fshift[],
2214            const t_pbc*    pbc,
2215            real            lambda,
2216            real*           dvdlambda,
2217            const t_mdatoms gmx_unused* md,
2218            t_fcdata gmx_unused* fcd,
2219            int gmx_unused* global_atom_index)
2220 {
2221     int  i, type, ai, aj, ak, al;
2222     int  t1, t2, t3;
2223     real phi, phi0, dphi0, ddphi, vtot;
2224     rvec r_ij, r_kj, r_kl, m, n;
2225     real L1, kk, dp, dp2, kA, kB, pA, pB, dvdl_term;
2226
2227     L1        = 1.0 - lambda;
2228     dvdl_term = 0;
2229     vtot      = 0.0;
2230     for (i = 0; (i < nbonds);)
2231     {
2232         type = forceatoms[i++];
2233         ai   = forceatoms[i++];
2234         aj   = forceatoms[i++];
2235         ak   = forceatoms[i++];
2236         al   = forceatoms[i++];
2237
2238         phi = dih_angle(x[ai], x[aj], x[ak], x[al], pbc, r_ij, r_kj, r_kl, m, n, &t1, &t2, &t3); /*  84         */
2239
2240         /* phi can jump if phi0 is close to Pi/-Pi, which will cause huge
2241          * force changes if we just apply a normal harmonic.
2242          * Instead, we first calculate phi-phi0 and take it modulo (-Pi,Pi).
2243          * This means we will never have the periodicity problem, unless
2244          * the dihedral is Pi away from phiO, which is very unlikely due to
2245          * the potential.
2246          */
2247         kA = forceparams[type].harmonic.krA;
2248         kB = forceparams[type].harmonic.krB;
2249         pA = forceparams[type].harmonic.rA;
2250         pB = forceparams[type].harmonic.rB;
2251
2252         kk    = L1 * kA + lambda * kB;
2253         phi0  = (L1 * pA + lambda * pB) * DEG2RAD;
2254         dphi0 = (pB - pA) * DEG2RAD;
2255
2256         dp = phi - phi0;
2257
2258         make_dp_periodic(&dp);
2259
2260         dp2 = dp * dp;
2261
2262         vtot += 0.5 * kk * dp2;
2263         ddphi = -kk * dp;
2264
2265         dvdl_term += 0.5 * (kB - kA) * dp2 - kk * dphi0 * dp;
2266
2267         do_dih_fup<flavor>(ai, aj, ak, al, -ddphi, r_ij, r_kj, r_kl, m, n, f, fshift, pbc, x, t1,
2268                            t2, t3); /* 112              */
2269         /* 218 TOTAL    */
2270     }
2271
2272     *dvdlambda += dvdl_term;
2273     return vtot;
2274 }
2275
2276 /*! \brief Computes angle restraints of two different types */
2277 template<BondedKernelFlavor flavor>
2278 real low_angres(int             nbonds,
2279                 const t_iatom   forceatoms[],
2280                 const t_iparams forceparams[],
2281                 const rvec      x[],
2282                 rvec4           f[],
2283                 rvec            fshift[],
2284                 const t_pbc*    pbc,
2285                 real            lambda,
2286                 real*           dvdlambda,
2287                 gmx_bool        bZAxis)
2288 {
2289     int  i, m, type, ai, aj, ak, al;
2290     int  t1, t2;
2291     real phi, cos_phi, cos_phi2, vid, vtot, dVdphi;
2292     rvec r_ij, r_kl, f_i, f_k = { 0, 0, 0 };
2293     real st, sth, nrij2, nrkl2, c, cij, ckl;
2294
2295     t2 = 0; /* avoid warning with gcc-3.3. It is never used uninitialized */
2296
2297     vtot = 0.0;
2298     ak = al = 0; /* to avoid warnings */
2299     for (i = 0; i < nbonds;)
2300     {
2301         type = forceatoms[i++];
2302         ai   = forceatoms[i++];
2303         aj   = forceatoms[i++];
2304         t1   = pbc_rvec_sub(pbc, x[aj], x[ai], r_ij); /*  3             */
2305         if (!bZAxis)
2306         {
2307             ak = forceatoms[i++];
2308             al = forceatoms[i++];
2309             t2 = pbc_rvec_sub(pbc, x[al], x[ak], r_kl); /*  3           */
2310         }
2311         else
2312         {
2313             r_kl[XX] = 0;
2314             r_kl[YY] = 0;
2315             r_kl[ZZ] = 1;
2316         }
2317
2318         cos_phi = cos_angle(r_ij, r_kl); /* 25          */
2319         phi     = std::acos(cos_phi);    /* 10           */
2320
2321         *dvdlambda += dopdihs_min(forceparams[type].pdihs.cpA, forceparams[type].pdihs.cpB,
2322                                   forceparams[type].pdihs.phiA, forceparams[type].pdihs.phiB,
2323                                   forceparams[type].pdihs.mult, phi, lambda, &vid, &dVdphi); /*  40 */
2324
2325         vtot += vid;
2326
2327         cos_phi2 = gmx::square(cos_phi); /*   1         */
2328         if (cos_phi2 < 1)
2329         {
2330             st    = -dVdphi * gmx::invsqrt(1 - cos_phi2); /*  12                */
2331             sth   = st * cos_phi;                         /*   1                */
2332             nrij2 = iprod(r_ij, r_ij);                    /*   5                */
2333             nrkl2 = iprod(r_kl, r_kl);                    /*   5          */
2334
2335             c   = st * gmx::invsqrt(nrij2 * nrkl2); /*  11              */
2336             cij = sth / nrij2;                      /*  10              */
2337             ckl = sth / nrkl2;                      /*  10              */
2338
2339             for (m = 0; m < DIM; m++) /*  18+18       */
2340             {
2341                 f_i[m] = (c * r_kl[m] - cij * r_ij[m]);
2342                 f[ai][m] += f_i[m];
2343                 f[aj][m] -= f_i[m];
2344                 if (!bZAxis)
2345                 {
2346                     f_k[m] = (c * r_ij[m] - ckl * r_kl[m]);
2347                     f[ak][m] += f_k[m];
2348                     f[al][m] -= f_k[m];
2349                 }
2350             }
2351
2352             if (computeVirial(flavor))
2353             {
2354                 rvec_inc(fshift[t1], f_i);
2355                 rvec_dec(fshift[CENTRAL], f_i);
2356                 if (!bZAxis)
2357                 {
2358                     rvec_inc(fshift[t2], f_k);
2359                     rvec_dec(fshift[CENTRAL], f_k);
2360                 }
2361             }
2362         }
2363     }
2364
2365     return vtot; /*  184 / 157 (bZAxis)  total  */
2366 }
2367
2368 template<BondedKernelFlavor flavor>
2369 real angres(int             nbonds,
2370             const t_iatom   forceatoms[],
2371             const t_iparams forceparams[],
2372             const rvec      x[],
2373             rvec4           f[],
2374             rvec            fshift[],
2375             const t_pbc*    pbc,
2376             real            lambda,
2377             real*           dvdlambda,
2378             const t_mdatoms gmx_unused* md,
2379             t_fcdata gmx_unused* fcd,
2380             int gmx_unused* global_atom_index)
2381 {
2382     return low_angres<flavor>(nbonds, forceatoms, forceparams, x, f, fshift, pbc, lambda, dvdlambda, FALSE);
2383 }
2384
2385 template<BondedKernelFlavor flavor>
2386 real angresz(int             nbonds,
2387              const t_iatom   forceatoms[],
2388              const t_iparams forceparams[],
2389              const rvec      x[],
2390              rvec4           f[],
2391              rvec            fshift[],
2392              const t_pbc*    pbc,
2393              real            lambda,
2394              real*           dvdlambda,
2395              const t_mdatoms gmx_unused* md,
2396              t_fcdata gmx_unused* fcd,
2397              int gmx_unused* global_atom_index)
2398 {
2399     return low_angres<flavor>(nbonds, forceatoms, forceparams, x, f, fshift, pbc, lambda, dvdlambda, TRUE);
2400 }
2401
2402 template<BondedKernelFlavor flavor>
2403 real dihres(int             nbonds,
2404             const t_iatom   forceatoms[],
2405             const t_iparams forceparams[],
2406             const rvec      x[],
2407             rvec4           f[],
2408             rvec            fshift[],
2409             const t_pbc*    pbc,
2410             real            lambda,
2411             real*           dvdlambda,
2412             const t_mdatoms gmx_unused* md,
2413             t_fcdata gmx_unused* fcd,
2414             int gmx_unused* global_atom_index)
2415 {
2416     real vtot = 0;
2417     int  ai, aj, ak, al, i, type, t1, t2, t3;
2418     real phi0A, phi0B, dphiA, dphiB, kfacA, kfacB, phi0, dphi, kfac;
2419     real phi, ddphi, ddp, ddp2, dp, d2r, L1;
2420     rvec r_ij, r_kj, r_kl, m, n;
2421
2422     L1 = 1.0 - lambda;
2423
2424     d2r = DEG2RAD;
2425
2426     for (i = 0; (i < nbonds);)
2427     {
2428         type = forceatoms[i++];
2429         ai   = forceatoms[i++];
2430         aj   = forceatoms[i++];
2431         ak   = forceatoms[i++];
2432         al   = forceatoms[i++];
2433
2434         phi0A = forceparams[type].dihres.phiA * d2r;
2435         dphiA = forceparams[type].dihres.dphiA * d2r;
2436         kfacA = forceparams[type].dihres.kfacA;
2437
2438         phi0B = forceparams[type].dihres.phiB * d2r;
2439         dphiB = forceparams[type].dihres.dphiB * d2r;
2440         kfacB = forceparams[type].dihres.kfacB;
2441
2442         phi0 = L1 * phi0A + lambda * phi0B;
2443         dphi = L1 * dphiA + lambda * dphiB;
2444         kfac = L1 * kfacA + lambda * kfacB;
2445
2446         phi = dih_angle(x[ai], x[aj], x[ak], x[al], pbc, r_ij, r_kj, r_kl, m, n, &t1, &t2, &t3);
2447         /* 84 flops */
2448
2449         /* phi can jump if phi0 is close to Pi/-Pi, which will cause huge
2450          * force changes if we just apply a normal harmonic.
2451          * Instead, we first calculate phi-phi0 and take it modulo (-Pi,Pi).
2452          * This means we will never have the periodicity problem, unless
2453          * the dihedral is Pi away from phiO, which is very unlikely due to
2454          * the potential.
2455          */
2456         dp = phi - phi0;
2457         make_dp_periodic(&dp);
2458
2459         if (dp > dphi)
2460         {
2461             ddp = dp - dphi;
2462         }
2463         else if (dp < -dphi)
2464         {
2465             ddp = dp + dphi;
2466         }
2467         else
2468         {
2469             ddp = 0;
2470         }
2471
2472         if (ddp != 0.0)
2473         {
2474             ddp2 = ddp * ddp;
2475             vtot += 0.5 * kfac * ddp2;
2476             ddphi = kfac * ddp;
2477
2478             *dvdlambda += 0.5 * (kfacB - kfacA) * ddp2;
2479             /* lambda dependence from changing restraint distances */
2480             if (ddp > 0)
2481             {
2482                 *dvdlambda -= kfac * ddp * ((dphiB - dphiA) + (phi0B - phi0A));
2483             }
2484             else if (ddp < 0)
2485             {
2486                 *dvdlambda += kfac * ddp * ((dphiB - dphiA) - (phi0B - phi0A));
2487             }
2488             do_dih_fup<flavor>(ai, aj, ak, al, ddphi, r_ij, r_kj, r_kl, m, n, f, fshift, pbc, x, t1,
2489                                t2, t3); /* 112          */
2490         }
2491     }
2492     return vtot;
2493 }
2494
2495
2496 real unimplemented(int gmx_unused nbonds,
2497                    const t_iatom gmx_unused forceatoms[],
2498                    const t_iparams gmx_unused forceparams[],
2499                    const rvec gmx_unused x[],
2500                    rvec4 gmx_unused f[],
2501                    rvec gmx_unused fshift[],
2502                    const t_pbc gmx_unused* pbc,
2503                    real gmx_unused lambda,
2504                    real gmx_unused* dvdlambda,
2505                    const t_mdatoms gmx_unused* md,
2506                    t_fcdata gmx_unused* fcd,
2507                    int gmx_unused* global_atom_index)
2508 {
2509     gmx_impl("*** you are using a not implemented function");
2510 }
2511
2512 template<BondedKernelFlavor flavor>
2513 real restrangles(int             nbonds,
2514                  const t_iatom   forceatoms[],
2515                  const t_iparams forceparams[],
2516                  const rvec      x[],
2517                  rvec4           f[],
2518                  rvec            fshift[],
2519                  const t_pbc*    pbc,
2520                  real gmx_unused lambda,
2521                  real gmx_unused* dvdlambda,
2522                  const t_mdatoms gmx_unused* md,
2523                  t_fcdata gmx_unused* fcd,
2524                  int gmx_unused* global_atom_index)
2525 {
2526     int    i, d, ai, aj, ak, type, m;
2527     int    t1, t2;
2528     real   v, vtot;
2529     rvec   f_i, f_j, f_k;
2530     double prefactor, ratio_ante, ratio_post;
2531     rvec   delta_ante, delta_post, vec_temp;
2532
2533     vtot = 0.0;
2534     for (i = 0; (i < nbonds);)
2535     {
2536         type = forceatoms[i++];
2537         ai   = forceatoms[i++];
2538         aj   = forceatoms[i++];
2539         ak   = forceatoms[i++];
2540
2541         t1 = pbc_rvec_sub(pbc, x[ai], x[aj], vec_temp);
2542         pbc_rvec_sub(pbc, x[aj], x[ai], delta_ante);
2543         t2 = pbc_rvec_sub(pbc, x[ak], x[aj], delta_post);
2544
2545
2546         /* This function computes factors needed for restricted angle potential.
2547          * The restricted angle potential is used in coarse-grained simulations to avoid singularities
2548          * when three particles align and the dihedral angle and dihedral potential
2549          * cannot be calculated. This potential is calculated using the formula:
2550            real restrangles(int nbonds,
2551             const t_iatom forceatoms[],const t_iparams forceparams[],
2552             const rvec x[],rvec4 f[],rvec fshift[],
2553             const t_pbc *pbc,
2554             real gmx_unused lambda,real gmx_unused *dvdlambda,
2555             const t_mdatoms gmx_unused *md,t_fcdata gmx_unused *fcd,
2556             int gmx_unused *global_atom_index)
2557            {
2558            int  i, d, ai, aj, ak, type, m;
2559            int t1, t2;
2560            rvec r_ij,r_kj;
2561            real v, vtot;
2562            rvec f_i, f_j, f_k;
2563            real prefactor, ratio_ante, ratio_post;
2564            rvec delta_ante, delta_post, vec_temp;
2565
2566            vtot = 0.0;
2567            for(i=0; (i<nbonds); )
2568            {
2569            type = forceatoms[i++];
2570            ai   = forceatoms[i++];
2571            aj   = forceatoms[i++];
2572            ak   = forceatoms[i++];
2573
2574          * \f[V_{\rm ReB}(\theta_i) = \frac{1}{2} k_{\theta} \frac{(\cos\theta_i - \cos\theta_0)^2}
2575          * {\sin^2\theta_i}\f] ({eq:ReB} and ref \cite{MonicaGoga2013} from the manual).
2576          * For more explanations see comments file "restcbt.h". */
2577
2578         compute_factors_restangles(type, forceparams, delta_ante, delta_post, &prefactor,
2579                                    &ratio_ante, &ratio_post, &v);
2580
2581         /*   Forces are computed per component */
2582         for (d = 0; d < DIM; d++)
2583         {
2584             f_i[d] = prefactor * (ratio_ante * delta_ante[d] - delta_post[d]);
2585             f_j[d] = prefactor
2586                      * ((ratio_post + 1.0) * delta_post[d] - (ratio_ante + 1.0) * delta_ante[d]);
2587             f_k[d] = prefactor * (delta_ante[d] - ratio_post * delta_post[d]);
2588         }
2589
2590         /*   Computation of potential energy   */
2591
2592         vtot += v;
2593
2594         /*   Update forces */
2595
2596         for (m = 0; (m < DIM); m++)
2597         {
2598             f[ai][m] += f_i[m];
2599             f[aj][m] += f_j[m];
2600             f[ak][m] += f_k[m];
2601         }
2602
2603         if (computeVirial(flavor))
2604         {
2605             rvec_inc(fshift[t1], f_i);
2606             rvec_inc(fshift[CENTRAL], f_j);
2607             rvec_inc(fshift[t2], f_k);
2608         }
2609     }
2610     return vtot;
2611 }
2612
2613
2614 template<BondedKernelFlavor flavor>
2615 real restrdihs(int             nbonds,
2616                const t_iatom   forceatoms[],
2617                const t_iparams forceparams[],
2618                const rvec      x[],
2619                rvec4           f[],
2620                rvec            fshift[],
2621                const t_pbc*    pbc,
2622                real gmx_unused lambda,
2623                real gmx_unused* dvlambda,
2624                const t_mdatoms gmx_unused* md,
2625                t_fcdata gmx_unused* fcd,
2626                int gmx_unused* global_atom_index)
2627 {
2628     int  i, d, type, ai, aj, ak, al;
2629     rvec f_i, f_j, f_k, f_l;
2630     rvec dx_jl;
2631     int  t1, t2, t3;
2632     real v, vtot;
2633     rvec delta_ante, delta_crnt, delta_post, vec_temp;
2634     real factor_phi_ai_ante, factor_phi_ai_crnt, factor_phi_ai_post;
2635     real factor_phi_aj_ante, factor_phi_aj_crnt, factor_phi_aj_post;
2636     real factor_phi_ak_ante, factor_phi_ak_crnt, factor_phi_ak_post;
2637     real factor_phi_al_ante, factor_phi_al_crnt, factor_phi_al_post;
2638     real prefactor_phi;
2639
2640
2641     vtot = 0.0;
2642     for (i = 0; (i < nbonds);)
2643     {
2644         type = forceatoms[i++];
2645         ai   = forceatoms[i++];
2646         aj   = forceatoms[i++];
2647         ak   = forceatoms[i++];
2648         al   = forceatoms[i++];
2649
2650         t1 = pbc_rvec_sub(pbc, x[ai], x[aj], vec_temp);
2651         pbc_rvec_sub(pbc, x[aj], x[ai], delta_ante);
2652         t2 = pbc_rvec_sub(pbc, x[ak], x[aj], delta_crnt);
2653         pbc_rvec_sub(pbc, x[ak], x[al], vec_temp);
2654         pbc_rvec_sub(pbc, x[al], x[ak], delta_post);
2655
2656         /* This function computes factors needed for restricted angle potential.
2657          * The restricted angle potential is used in coarse-grained simulations to avoid singularities
2658          * when three particles align and the dihedral angle and dihedral potential cannot be calculated.
2659          * This potential is calculated using the formula:
2660          * \f[V_{\rm ReB}(\theta_i) = \frac{1}{2} k_{\theta}
2661          * \frac{(\cos\theta_i - \cos\theta_0)^2}{\sin^2\theta_i}\f]
2662          * ({eq:ReB} and ref \cite{MonicaGoga2013} from the manual).
2663          * For more explanations see comments file "restcbt.h" */
2664
2665         compute_factors_restrdihs(
2666                 type, forceparams, delta_ante, delta_crnt, delta_post, &factor_phi_ai_ante,
2667                 &factor_phi_ai_crnt, &factor_phi_ai_post, &factor_phi_aj_ante, &factor_phi_aj_crnt,
2668                 &factor_phi_aj_post, &factor_phi_ak_ante, &factor_phi_ak_crnt, &factor_phi_ak_post,
2669                 &factor_phi_al_ante, &factor_phi_al_crnt, &factor_phi_al_post, &prefactor_phi, &v);
2670
2671
2672         /*      Computation of forces per component */
2673         for (d = 0; d < DIM; d++)
2674         {
2675             f_i[d] = prefactor_phi
2676                      * (factor_phi_ai_ante * delta_ante[d] + factor_phi_ai_crnt * delta_crnt[d]
2677                         + factor_phi_ai_post * delta_post[d]);
2678             f_j[d] = prefactor_phi
2679                      * (factor_phi_aj_ante * delta_ante[d] + factor_phi_aj_crnt * delta_crnt[d]
2680                         + factor_phi_aj_post * delta_post[d]);
2681             f_k[d] = prefactor_phi
2682                      * (factor_phi_ak_ante * delta_ante[d] + factor_phi_ak_crnt * delta_crnt[d]
2683                         + factor_phi_ak_post * delta_post[d]);
2684             f_l[d] = prefactor_phi
2685                      * (factor_phi_al_ante * delta_ante[d] + factor_phi_al_crnt * delta_crnt[d]
2686                         + factor_phi_al_post * delta_post[d]);
2687         }
2688         /*      Computation of the energy */
2689
2690         vtot += v;
2691
2692
2693         /*    Updating the forces */
2694
2695         rvec_inc(f[ai], f_i);
2696         rvec_inc(f[aj], f_j);
2697         rvec_inc(f[ak], f_k);
2698         rvec_inc(f[al], f_l);
2699
2700         if (computeVirial(flavor))
2701         {
2702             if (pbc)
2703             {
2704                 t3 = pbc_rvec_sub(pbc, x[al], x[aj], dx_jl);
2705             }
2706             else
2707             {
2708                 t3 = CENTRAL;
2709             }
2710
2711             rvec_inc(fshift[t1], f_i);
2712             rvec_inc(fshift[CENTRAL], f_j);
2713             rvec_inc(fshift[t2], f_k);
2714             rvec_inc(fshift[t3], f_l);
2715         }
2716     }
2717
2718     return vtot;
2719 }
2720
2721
2722 template<BondedKernelFlavor flavor>
2723 real cbtdihs(int             nbonds,
2724              const t_iatom   forceatoms[],
2725              const t_iparams forceparams[],
2726              const rvec      x[],
2727              rvec4           f[],
2728              rvec            fshift[],
2729              const t_pbc*    pbc,
2730              real gmx_unused lambda,
2731              real gmx_unused* dvdlambda,
2732              const t_mdatoms gmx_unused* md,
2733              t_fcdata gmx_unused* fcd,
2734              int gmx_unused* global_atom_index)
2735 {
2736     int  type, ai, aj, ak, al, i, d;
2737     int  t1, t2, t3;
2738     real v, vtot;
2739     rvec vec_temp;
2740     rvec f_i, f_j, f_k, f_l;
2741     rvec dx_jl;
2742     rvec delta_ante, delta_crnt, delta_post;
2743     rvec f_phi_ai, f_phi_aj, f_phi_ak, f_phi_al;
2744     rvec f_theta_ante_ai, f_theta_ante_aj, f_theta_ante_ak;
2745     rvec f_theta_post_aj, f_theta_post_ak, f_theta_post_al;
2746
2747
2748     vtot = 0.0;
2749     for (i = 0; (i < nbonds);)
2750     {
2751         type = forceatoms[i++];
2752         ai   = forceatoms[i++];
2753         aj   = forceatoms[i++];
2754         ak   = forceatoms[i++];
2755         al   = forceatoms[i++];
2756
2757
2758         t1 = pbc_rvec_sub(pbc, x[ai], x[aj], vec_temp);
2759         pbc_rvec_sub(pbc, x[aj], x[ai], delta_ante);
2760         t2 = pbc_rvec_sub(pbc, x[ak], x[aj], vec_temp);
2761         pbc_rvec_sub(pbc, x[ak], x[aj], delta_crnt);
2762         pbc_rvec_sub(pbc, x[ak], x[al], vec_temp);
2763         pbc_rvec_sub(pbc, x[al], x[ak], delta_post);
2764
2765         /* \brief Compute factors for CBT potential
2766          * The combined bending-torsion potential goes to zero in a very smooth manner, eliminating the numerical
2767          * instabilities, when three coarse-grained particles align and the dihedral angle and standard
2768          * dihedral potentials cannot be calculated. The CBT potential is calculated using the formula:
2769          * \f[V_{\rm CBT}(\theta_{i-1}, \theta_i, \phi_i) = k_{\phi} \sin^3\theta_{i-1} \sin^3\theta_{i}
2770          * \sum_{n=0}^4 { a_n \cos^n\phi_i}\f] ({eq:CBT} and ref \cite{MonicaGoga2013} from the manual).
2771          * It contains in its expression not only the dihedral angle \f$\phi\f$
2772          * but also \f[\theta_{i-1}\f] (theta_ante bellow) and \f[\theta_{i}\f] (theta_post bellow)
2773          * --- the adjacent bending angles.
2774          * For more explanations see comments file "restcbt.h". */
2775
2776         compute_factors_cbtdihs(type, forceparams, delta_ante, delta_crnt, delta_post, f_phi_ai,
2777                                 f_phi_aj, f_phi_ak, f_phi_al, f_theta_ante_ai, f_theta_ante_aj,
2778                                 f_theta_ante_ak, f_theta_post_aj, f_theta_post_ak, f_theta_post_al, &v);
2779
2780
2781         /*      Acumulate the resuts per beads */
2782         for (d = 0; d < DIM; d++)
2783         {
2784             f_i[d] = f_phi_ai[d] + f_theta_ante_ai[d];
2785             f_j[d] = f_phi_aj[d] + f_theta_ante_aj[d] + f_theta_post_aj[d];
2786             f_k[d] = f_phi_ak[d] + f_theta_ante_ak[d] + f_theta_post_ak[d];
2787             f_l[d] = f_phi_al[d] + f_theta_post_al[d];
2788         }
2789
2790         /*      Compute the potential energy */
2791
2792         vtot += v;
2793
2794
2795         /*  Updating the forces */
2796         rvec_inc(f[ai], f_i);
2797         rvec_inc(f[aj], f_j);
2798         rvec_inc(f[ak], f_k);
2799         rvec_inc(f[al], f_l);
2800
2801
2802         if (computeVirial(flavor))
2803         {
2804             if (pbc)
2805             {
2806                 t3 = pbc_rvec_sub(pbc, x[al], x[aj], dx_jl);
2807             }
2808             else
2809             {
2810                 t3 = CENTRAL;
2811             }
2812
2813             rvec_inc(fshift[t1], f_i);
2814             rvec_inc(fshift[CENTRAL], f_j);
2815             rvec_inc(fshift[t2], f_k);
2816             rvec_inc(fshift[t3], f_l);
2817         }
2818     }
2819
2820     return vtot;
2821 }
2822
2823 template<BondedKernelFlavor flavor>
2824 std::enable_if_t<flavor != BondedKernelFlavor::ForcesSimdWhenAvailable || !GMX_SIMD_HAVE_REAL, real>
2825 rbdihs(int             nbonds,
2826        const t_iatom   forceatoms[],
2827        const t_iparams forceparams[],
2828        const rvec      x[],
2829        rvec4           f[],
2830        rvec            fshift[],
2831        const t_pbc*    pbc,
2832        real            lambda,
2833        real*           dvdlambda,
2834        const t_mdatoms gmx_unused* md,
2835        t_fcdata gmx_unused* fcd,
2836        int gmx_unused* global_atom_index)
2837 {
2838     const real c0 = 0.0, c1 = 1.0, c2 = 2.0, c3 = 3.0, c4 = 4.0, c5 = 5.0;
2839     int        type, ai, aj, ak, al, i, j;
2840     int        t1, t2, t3;
2841     rvec       r_ij, r_kj, r_kl, m, n;
2842     real       parmA[NR_RBDIHS];
2843     real       parmB[NR_RBDIHS];
2844     real       parm[NR_RBDIHS];
2845     real       cos_phi, phi, rbp, rbpBA;
2846     real       v, ddphi, sin_phi;
2847     real       cosfac, vtot;
2848     real       L1        = 1.0 - lambda;
2849     real       dvdl_term = 0;
2850
2851     vtot = 0.0;
2852     for (i = 0; (i < nbonds);)
2853     {
2854         type = forceatoms[i++];
2855         ai   = forceatoms[i++];
2856         aj   = forceatoms[i++];
2857         ak   = forceatoms[i++];
2858         al   = forceatoms[i++];
2859
2860         phi = dih_angle(x[ai], x[aj], x[ak], x[al], pbc, r_ij, r_kj, r_kl, m, n, &t1, &t2, &t3); /*  84         */
2861
2862         /* Change to polymer convention */
2863         if (phi < c0)
2864         {
2865             phi += M_PI;
2866         }
2867         else
2868         {
2869             phi -= M_PI; /*   1         */
2870         }
2871         cos_phi = std::cos(phi);
2872         /* Beware of accuracy loss, cannot use 1-sqrt(cos^2) ! */
2873         sin_phi = std::sin(phi);
2874
2875         for (j = 0; (j < NR_RBDIHS); j++)
2876         {
2877             parmA[j] = forceparams[type].rbdihs.rbcA[j];
2878             parmB[j] = forceparams[type].rbdihs.rbcB[j];
2879             parm[j]  = L1 * parmA[j] + lambda * parmB[j];
2880         }
2881         /* Calculate cosine powers */
2882         /* Calculate the energy */
2883         /* Calculate the derivative */
2884
2885         v = parm[0];
2886         dvdl_term += (parmB[0] - parmA[0]);
2887         ddphi  = c0;
2888         cosfac = c1;
2889
2890         rbp   = parm[1];
2891         rbpBA = parmB[1] - parmA[1];
2892         ddphi += rbp * cosfac;
2893         cosfac *= cos_phi;
2894         v += cosfac * rbp;
2895         dvdl_term += cosfac * rbpBA;
2896         rbp   = parm[2];
2897         rbpBA = parmB[2] - parmA[2];
2898         ddphi += c2 * rbp * cosfac;
2899         cosfac *= cos_phi;
2900         v += cosfac * rbp;
2901         dvdl_term += cosfac * rbpBA;
2902         rbp   = parm[3];
2903         rbpBA = parmB[3] - parmA[3];
2904         ddphi += c3 * rbp * cosfac;
2905         cosfac *= cos_phi;
2906         v += cosfac * rbp;
2907         dvdl_term += cosfac * rbpBA;
2908         rbp   = parm[4];
2909         rbpBA = parmB[4] - parmA[4];
2910         ddphi += c4 * rbp * cosfac;
2911         cosfac *= cos_phi;
2912         v += cosfac * rbp;
2913         dvdl_term += cosfac * rbpBA;
2914         rbp   = parm[5];
2915         rbpBA = parmB[5] - parmA[5];
2916         ddphi += c5 * rbp * cosfac;
2917         cosfac *= cos_phi;
2918         v += cosfac * rbp;
2919         dvdl_term += cosfac * rbpBA;
2920
2921         ddphi = -ddphi * sin_phi; /*  11                */
2922
2923         do_dih_fup<flavor>(ai, aj, ak, al, ddphi, r_ij, r_kj, r_kl, m, n, f, fshift, pbc, x, t1, t2,
2924                            t3); /* 112          */
2925         vtot += v;
2926     }
2927     *dvdlambda += dvdl_term;
2928
2929     return vtot;
2930 }
2931
2932 //! \endcond
2933
2934 /*! \brief Mysterious undocumented function */
2935 int cmap_setup_grid_index(int ip, int grid_spacing, int* ipm1, int* ipp1, int* ipp2)
2936 {
2937     int im1, ip1, ip2;
2938
2939     if (ip < 0)
2940     {
2941         ip = ip + grid_spacing - 1;
2942     }
2943     else if (ip > grid_spacing)
2944     {
2945         ip = ip - grid_spacing - 1;
2946     }
2947
2948     im1 = ip - 1;
2949     ip1 = ip + 1;
2950     ip2 = ip + 2;
2951
2952     if (ip == 0)
2953     {
2954         im1 = grid_spacing - 1;
2955     }
2956     else if (ip == grid_spacing - 2)
2957     {
2958         ip2 = 0;
2959     }
2960     else if (ip == grid_spacing - 1)
2961     {
2962         ip1 = 0;
2963         ip2 = 1;
2964     }
2965
2966     *ipm1 = im1;
2967     *ipp1 = ip1;
2968     *ipp2 = ip2;
2969
2970     return ip;
2971 }
2972
2973 } // namespace
2974
2975 real cmap_dihs(int                 nbonds,
2976                const t_iatom       forceatoms[],
2977                const t_iparams     forceparams[],
2978                const gmx_cmap_t*   cmap_grid,
2979                const rvec          x[],
2980                rvec4               f[],
2981                rvec                fshift[],
2982                const struct t_pbc* pbc,
2983                real gmx_unused lambda,
2984                real gmx_unused* dvdlambda,
2985                const t_mdatoms gmx_unused* md,
2986                t_fcdata gmx_unused* fcd,
2987                int gmx_unused* global_atom_index)
2988 {
2989     int i, n;
2990     int ai, aj, ak, al, am;
2991     int a1i, a1j, a1k, a1l, a2i, a2j, a2k, a2l;
2992     int type;
2993     int t11, t21, t31, t12, t22, t32;
2994     int iphi1, ip1m1, ip1p1, ip1p2;
2995     int iphi2, ip2m1, ip2p1, ip2p2;
2996     int l1, l2, l3;
2997     int pos1, pos2, pos3, pos4;
2998
2999     real ty[4], ty1[4], ty2[4], ty12[4], tx[16];
3000     real phi1, cos_phi1, sin_phi1, xphi1;
3001     real phi2, cos_phi2, sin_phi2, xphi2;
3002     real dx, tt, tu, e, df1, df2, vtot;
3003     real ra21, rb21, rg21, rg1, rgr1, ra2r1, rb2r1, rabr1;
3004     real ra22, rb22, rg22, rg2, rgr2, ra2r2, rb2r2, rabr2;
3005     real fg1, hg1, fga1, hgb1, gaa1, gbb1;
3006     real fg2, hg2, fga2, hgb2, gaa2, gbb2;
3007     real fac;
3008
3009     rvec r1_ij, r1_kj, r1_kl, m1, n1;
3010     rvec r2_ij, r2_kj, r2_kl, m2, n2;
3011     rvec f1_i, f1_j, f1_k, f1_l;
3012     rvec f2_i, f2_j, f2_k, f2_l;
3013     rvec a1, b1, a2, b2;
3014     rvec f1, g1, h1, f2, g2, h2;
3015     rvec dtf1, dtg1, dth1, dtf2, dtg2, dth2;
3016
3017     int loop_index[4][4] = { { 0, 4, 8, 12 }, { 1, 5, 9, 13 }, { 2, 6, 10, 14 }, { 3, 7, 11, 15 } };
3018
3019     /* Total CMAP energy */
3020     vtot = 0;
3021
3022     for (n = 0; n < nbonds;)
3023     {
3024         /* Five atoms are involved in the two torsions */
3025         type = forceatoms[n++];
3026         ai   = forceatoms[n++];
3027         aj   = forceatoms[n++];
3028         ak   = forceatoms[n++];
3029         al   = forceatoms[n++];
3030         am   = forceatoms[n++];
3031
3032         /* Which CMAP type is this */
3033         const int   cmapA = forceparams[type].cmap.cmapA;
3034         const real* cmapd = cmap_grid->cmapdata[cmapA].cmap.data();
3035
3036         /* First torsion */
3037         a1i = ai;
3038         a1j = aj;
3039         a1k = ak;
3040         a1l = al;
3041
3042         phi1 = dih_angle(x[a1i], x[a1j], x[a1k], x[a1l], pbc, r1_ij, r1_kj, r1_kl, m1, n1, &t11,
3043                          &t21, &t31); /* 84 */
3044
3045         cos_phi1 = std::cos(phi1);
3046
3047         a1[0] = r1_ij[1] * r1_kj[2] - r1_ij[2] * r1_kj[1];
3048         a1[1] = r1_ij[2] * r1_kj[0] - r1_ij[0] * r1_kj[2];
3049         a1[2] = r1_ij[0] * r1_kj[1] - r1_ij[1] * r1_kj[0]; /* 9 */
3050
3051         b1[0] = r1_kl[1] * r1_kj[2] - r1_kl[2] * r1_kj[1];
3052         b1[1] = r1_kl[2] * r1_kj[0] - r1_kl[0] * r1_kj[2];
3053         b1[2] = r1_kl[0] * r1_kj[1] - r1_kl[1] * r1_kj[0]; /* 9 */
3054
3055         pbc_rvec_sub(pbc, x[a1l], x[a1k], h1);
3056
3057         ra21 = iprod(a1, a1);       /* 5 */
3058         rb21 = iprod(b1, b1);       /* 5 */
3059         rg21 = iprod(r1_kj, r1_kj); /* 5 */
3060         rg1  = sqrt(rg21);
3061
3062         rgr1  = 1.0 / rg1;
3063         ra2r1 = 1.0 / ra21;
3064         rb2r1 = 1.0 / rb21;
3065         rabr1 = sqrt(ra2r1 * rb2r1);
3066
3067         sin_phi1 = rg1 * rabr1 * iprod(a1, h1) * (-1);
3068
3069         if (cos_phi1 < -0.5 || cos_phi1 > 0.5)
3070         {
3071             phi1 = std::asin(sin_phi1);
3072
3073             if (cos_phi1 < 0)
3074             {
3075                 if (phi1 > 0)
3076                 {
3077                     phi1 = M_PI - phi1;
3078                 }
3079                 else
3080                 {
3081                     phi1 = -M_PI - phi1;
3082                 }
3083             }
3084         }
3085         else
3086         {
3087             phi1 = std::acos(cos_phi1);
3088
3089             if (sin_phi1 < 0)
3090             {
3091                 phi1 = -phi1;
3092             }
3093         }
3094
3095         xphi1 = phi1 + M_PI; /* 1 */
3096
3097         /* Second torsion */
3098         a2i = aj;
3099         a2j = ak;
3100         a2k = al;
3101         a2l = am;
3102
3103         phi2 = dih_angle(x[a2i], x[a2j], x[a2k], x[a2l], pbc, r2_ij, r2_kj, r2_kl, m2, n2, &t12,
3104                          &t22, &t32); /* 84 */
3105
3106         cos_phi2 = std::cos(phi2);
3107
3108         a2[0] = r2_ij[1] * r2_kj[2] - r2_ij[2] * r2_kj[1];
3109         a2[1] = r2_ij[2] * r2_kj[0] - r2_ij[0] * r2_kj[2];
3110         a2[2] = r2_ij[0] * r2_kj[1] - r2_ij[1] * r2_kj[0]; /* 9 */
3111
3112         b2[0] = r2_kl[1] * r2_kj[2] - r2_kl[2] * r2_kj[1];
3113         b2[1] = r2_kl[2] * r2_kj[0] - r2_kl[0] * r2_kj[2];
3114         b2[2] = r2_kl[0] * r2_kj[1] - r2_kl[1] * r2_kj[0]; /* 9 */
3115
3116         pbc_rvec_sub(pbc, x[a2l], x[a2k], h2);
3117
3118         ra22 = iprod(a2, a2);       /* 5 */
3119         rb22 = iprod(b2, b2);       /* 5 */
3120         rg22 = iprod(r2_kj, r2_kj); /* 5 */
3121         rg2  = sqrt(rg22);
3122
3123         rgr2  = 1.0 / rg2;
3124         ra2r2 = 1.0 / ra22;
3125         rb2r2 = 1.0 / rb22;
3126         rabr2 = sqrt(ra2r2 * rb2r2);
3127
3128         sin_phi2 = rg2 * rabr2 * iprod(a2, h2) * (-1);
3129
3130         if (cos_phi2 < -0.5 || cos_phi2 > 0.5)
3131         {
3132             phi2 = std::asin(sin_phi2);
3133
3134             if (cos_phi2 < 0)
3135             {
3136                 if (phi2 > 0)
3137                 {
3138                     phi2 = M_PI - phi2;
3139                 }
3140                 else
3141                 {
3142                     phi2 = -M_PI - phi2;
3143                 }
3144             }
3145         }
3146         else
3147         {
3148             phi2 = std::acos(cos_phi2);
3149
3150             if (sin_phi2 < 0)
3151             {
3152                 phi2 = -phi2;
3153             }
3154         }
3155
3156         xphi2 = phi2 + M_PI; /* 1 */
3157
3158         /* Range mangling */
3159         if (xphi1 < 0)
3160         {
3161             xphi1 = xphi1 + 2 * M_PI;
3162         }
3163         else if (xphi1 >= 2 * M_PI)
3164         {
3165             xphi1 = xphi1 - 2 * M_PI;
3166         }
3167
3168         if (xphi2 < 0)
3169         {
3170             xphi2 = xphi2 + 2 * M_PI;
3171         }
3172         else if (xphi2 >= 2 * M_PI)
3173         {
3174             xphi2 = xphi2 - 2 * M_PI;
3175         }
3176
3177         /* Number of grid points */
3178         dx = 2 * M_PI / cmap_grid->grid_spacing;
3179
3180         /* Where on the grid are we */
3181         iphi1 = static_cast<int>(xphi1 / dx);
3182         iphi2 = static_cast<int>(xphi2 / dx);
3183
3184         iphi1 = cmap_setup_grid_index(iphi1, cmap_grid->grid_spacing, &ip1m1, &ip1p1, &ip1p2);
3185         iphi2 = cmap_setup_grid_index(iphi2, cmap_grid->grid_spacing, &ip2m1, &ip2p1, &ip2p2);
3186
3187         pos1 = iphi1 * cmap_grid->grid_spacing + iphi2;
3188         pos2 = ip1p1 * cmap_grid->grid_spacing + iphi2;
3189         pos3 = ip1p1 * cmap_grid->grid_spacing + ip2p1;
3190         pos4 = iphi1 * cmap_grid->grid_spacing + ip2p1;
3191
3192         ty[0] = cmapd[pos1 * 4];
3193         ty[1] = cmapd[pos2 * 4];
3194         ty[2] = cmapd[pos3 * 4];
3195         ty[3] = cmapd[pos4 * 4];
3196
3197         ty1[0] = cmapd[pos1 * 4 + 1];
3198         ty1[1] = cmapd[pos2 * 4 + 1];
3199         ty1[2] = cmapd[pos3 * 4 + 1];
3200         ty1[3] = cmapd[pos4 * 4 + 1];
3201
3202         ty2[0] = cmapd[pos1 * 4 + 2];
3203         ty2[1] = cmapd[pos2 * 4 + 2];
3204         ty2[2] = cmapd[pos3 * 4 + 2];
3205         ty2[3] = cmapd[pos4 * 4 + 2];
3206
3207         ty12[0] = cmapd[pos1 * 4 + 3];
3208         ty12[1] = cmapd[pos2 * 4 + 3];
3209         ty12[2] = cmapd[pos3 * 4 + 3];
3210         ty12[3] = cmapd[pos4 * 4 + 3];
3211
3212         /* Switch to degrees */
3213         dx    = 360.0 / cmap_grid->grid_spacing;
3214         xphi1 = xphi1 * RAD2DEG;
3215         xphi2 = xphi2 * RAD2DEG;
3216
3217         for (i = 0; i < 4; i++) /* 16 */
3218         {
3219             tx[i]      = ty[i];
3220             tx[i + 4]  = ty1[i] * dx;
3221             tx[i + 8]  = ty2[i] * dx;
3222             tx[i + 12] = ty12[i] * dx * dx;
3223         }
3224
3225         real tc[16] = { 0 };
3226         for (int idx = 0; idx < 16; idx++) /* 1056 */
3227         {
3228             for (int k = 0; k < 16; k++)
3229             {
3230                 tc[idx] += cmap_coeff_matrix[k * 16 + idx] * tx[k];
3231             }
3232         }
3233
3234         tt = (xphi1 - iphi1 * dx) / dx;
3235         tu = (xphi2 - iphi2 * dx) / dx;
3236
3237         e   = 0;
3238         df1 = 0;
3239         df2 = 0;
3240
3241         for (i = 3; i >= 0; i--)
3242         {
3243             l1 = loop_index[i][3];
3244             l2 = loop_index[i][2];
3245             l3 = loop_index[i][1];
3246
3247             e = tt * e + ((tc[i * 4 + 3] * tu + tc[i * 4 + 2]) * tu + tc[i * 4 + 1]) * tu + tc[i * 4];
3248             df1 = tu * df1 + (3.0 * tc[l1] * tt + 2.0 * tc[l2]) * tt + tc[l3];
3249             df2 = tt * df2 + (3.0 * tc[i * 4 + 3] * tu + 2.0 * tc[i * 4 + 2]) * tu + tc[i * 4 + 1];
3250         }
3251
3252         fac = RAD2DEG / dx;
3253         df1 = df1 * fac;
3254         df2 = df2 * fac;
3255
3256         /* CMAP energy */
3257         vtot += e;
3258
3259         /* Do forces - first torsion */
3260         fg1  = iprod(r1_ij, r1_kj);
3261         hg1  = iprod(r1_kl, r1_kj);
3262         fga1 = fg1 * ra2r1 * rgr1;
3263         hgb1 = hg1 * rb2r1 * rgr1;
3264         gaa1 = -ra2r1 * rg1;
3265         gbb1 = rb2r1 * rg1;
3266
3267         for (i = 0; i < DIM; i++)
3268         {
3269             dtf1[i] = gaa1 * a1[i];
3270             dtg1[i] = fga1 * a1[i] - hgb1 * b1[i];
3271             dth1[i] = gbb1 * b1[i];
3272
3273             f1[i] = df1 * dtf1[i];
3274             g1[i] = df1 * dtg1[i];
3275             h1[i] = df1 * dth1[i];
3276
3277             f1_i[i] = f1[i];
3278             f1_j[i] = -f1[i] - g1[i];
3279             f1_k[i] = h1[i] + g1[i];
3280             f1_l[i] = -h1[i];
3281
3282             f[a1i][i] = f[a1i][i] + f1_i[i];
3283             f[a1j][i] = f[a1j][i] + f1_j[i]; /* - f1[i] - g1[i] */
3284             f[a1k][i] = f[a1k][i] + f1_k[i]; /* h1[i] + g1[i] */
3285             f[a1l][i] = f[a1l][i] + f1_l[i]; /* h1[i] */
3286         }
3287
3288         /* Do forces - second torsion */
3289         fg2  = iprod(r2_ij, r2_kj);
3290         hg2  = iprod(r2_kl, r2_kj);
3291         fga2 = fg2 * ra2r2 * rgr2;
3292         hgb2 = hg2 * rb2r2 * rgr2;
3293         gaa2 = -ra2r2 * rg2;
3294         gbb2 = rb2r2 * rg2;
3295
3296         for (i = 0; i < DIM; i++)
3297         {
3298             dtf2[i] = gaa2 * a2[i];
3299             dtg2[i] = fga2 * a2[i] - hgb2 * b2[i];
3300             dth2[i] = gbb2 * b2[i];
3301
3302             f2[i] = df2 * dtf2[i];
3303             g2[i] = df2 * dtg2[i];
3304             h2[i] = df2 * dth2[i];
3305
3306             f2_i[i] = f2[i];
3307             f2_j[i] = -f2[i] - g2[i];
3308             f2_k[i] = h2[i] + g2[i];
3309             f2_l[i] = -h2[i];
3310
3311             f[a2i][i] = f[a2i][i] + f2_i[i]; /* f2[i] */
3312             f[a2j][i] = f[a2j][i] + f2_j[i]; /* - f2[i] - g2[i] */
3313             f[a2k][i] = f[a2k][i] + f2_k[i]; /* h2[i] + g2[i] */
3314             f[a2l][i] = f[a2l][i] + f2_l[i]; /* - h2[i] */
3315         }
3316
3317         /* Shift forces */
3318         if (fshift != nullptr)
3319         {
3320             if (pbc)
3321             {
3322                 t31 = pbc_rvec_sub(pbc, x[a1l], x[a1j], h1);
3323                 t32 = pbc_rvec_sub(pbc, x[a2l], x[a2j], h2);
3324             }
3325             else
3326             {
3327                 t31 = CENTRAL;
3328                 t32 = CENTRAL;
3329             }
3330
3331             rvec_inc(fshift[t11], f1_i);
3332             rvec_inc(fshift[CENTRAL], f1_j);
3333             rvec_inc(fshift[t21], f1_k);
3334             rvec_inc(fshift[t31], f1_l);
3335
3336             rvec_inc(fshift[t12], f2_i);
3337             rvec_inc(fshift[CENTRAL], f2_j);
3338             rvec_inc(fshift[t22], f2_k);
3339             rvec_inc(fshift[t32], f2_l);
3340         }
3341     }
3342     return vtot;
3343 }
3344
3345 namespace
3346 {
3347
3348 //! \cond
3349 /***********************************************************
3350  *
3351  *   G R O M O S  9 6   F U N C T I O N S
3352  *
3353  ***********************************************************/
3354 real g96harmonic(real kA, real kB, real xA, real xB, real x, real lambda, real* V, real* F)
3355 {
3356     const real half = 0.5;
3357     real       L1, kk, x0, dx, dx2;
3358     real       v, f, dvdlambda;
3359
3360     L1 = 1.0 - lambda;
3361     kk = L1 * kA + lambda * kB;
3362     x0 = L1 * xA + lambda * xB;
3363
3364     dx  = x - x0;
3365     dx2 = dx * dx;
3366
3367     f         = -kk * dx;
3368     v         = half * kk * dx2;
3369     dvdlambda = half * (kB - kA) * dx2 + (xA - xB) * kk * dx;
3370
3371     *F = f;
3372     *V = v;
3373
3374     return dvdlambda;
3375
3376     /* That was 21 flops */
3377 }
3378
3379 template<BondedKernelFlavor flavor>
3380 real g96bonds(int             nbonds,
3381               const t_iatom   forceatoms[],
3382               const t_iparams forceparams[],
3383               const rvec      x[],
3384               rvec4           f[],
3385               rvec            fshift[],
3386               const t_pbc*    pbc,
3387               real            lambda,
3388               real*           dvdlambda,
3389               const t_mdatoms gmx_unused* md,
3390               t_fcdata gmx_unused* fcd,
3391               int gmx_unused* global_atom_index)
3392 {
3393     int  i, ki, ai, aj, type;
3394     real dr2, fbond, vbond, vtot;
3395     rvec dx;
3396
3397     vtot = 0.0;
3398     for (i = 0; (i < nbonds);)
3399     {
3400         type = forceatoms[i++];
3401         ai   = forceatoms[i++];
3402         aj   = forceatoms[i++];
3403
3404         ki  = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /*   3      */
3405         dr2 = iprod(dx, dx);                       /*   5               */
3406
3407         *dvdlambda += g96harmonic(forceparams[type].harmonic.krA, forceparams[type].harmonic.krB,
3408                                   forceparams[type].harmonic.rA, forceparams[type].harmonic.rB, dr2,
3409                                   lambda, &vbond, &fbond);
3410
3411         vtot += 0.5 * vbond; /* 1*/
3412
3413         spreadBondForces<flavor>(fbond, dx, ai, aj, f, ki, fshift); /* 15 */
3414     }                                                               /* 44 TOTAL */
3415     return vtot;
3416 }
3417
3418 real g96bond_angle(const rvec xi, const rvec xj, const rvec xk, const t_pbc* pbc, rvec r_ij, rvec r_kj, int* t1, int* t2)
3419 /* Return value is the angle between the bonds i-j and j-k */
3420 {
3421     real costh;
3422
3423     *t1 = pbc_rvec_sub(pbc, xi, xj, r_ij); /*  3                */
3424     *t2 = pbc_rvec_sub(pbc, xk, xj, r_kj); /*  3                */
3425
3426     costh = cos_angle(r_ij, r_kj); /* 25                */
3427     /* 41 TOTAL */
3428     return costh;
3429 }
3430
3431 template<BondedKernelFlavor flavor>
3432 real g96angles(int             nbonds,
3433                const t_iatom   forceatoms[],
3434                const t_iparams forceparams[],
3435                const rvec      x[],
3436                rvec4           f[],
3437                rvec            fshift[],
3438                const t_pbc*    pbc,
3439                real            lambda,
3440                real*           dvdlambda,
3441                const t_mdatoms gmx_unused* md,
3442                t_fcdata gmx_unused* fcd,
3443                int gmx_unused* global_atom_index)
3444 {
3445     int  i, ai, aj, ak, type, m, t1, t2;
3446     rvec r_ij, r_kj;
3447     real cos_theta, dVdt, va, vtot;
3448     real rij_1, rij_2, rkj_1, rkj_2, rijrkj_1;
3449     rvec f_i, f_j, f_k;
3450
3451     vtot = 0.0;
3452     for (i = 0; (i < nbonds);)
3453     {
3454         type = forceatoms[i++];
3455         ai   = forceatoms[i++];
3456         aj   = forceatoms[i++];
3457         ak   = forceatoms[i++];
3458
3459         cos_theta = g96bond_angle(x[ai], x[aj], x[ak], pbc, r_ij, r_kj, &t1, &t2);
3460
3461         *dvdlambda += g96harmonic(forceparams[type].harmonic.krA, forceparams[type].harmonic.krB,
3462                                   forceparams[type].harmonic.rA, forceparams[type].harmonic.rB,
3463                                   cos_theta, lambda, &va, &dVdt);
3464         vtot += va;
3465
3466         rij_1    = gmx::invsqrt(iprod(r_ij, r_ij));
3467         rkj_1    = gmx::invsqrt(iprod(r_kj, r_kj));
3468         rij_2    = rij_1 * rij_1;
3469         rkj_2    = rkj_1 * rkj_1;
3470         rijrkj_1 = rij_1 * rkj_1; /* 23 */
3471
3472         for (m = 0; (m < DIM); m++) /*  42      */
3473         {
3474             f_i[m] = dVdt * (r_kj[m] * rijrkj_1 - r_ij[m] * rij_2 * cos_theta);
3475             f_k[m] = dVdt * (r_ij[m] * rijrkj_1 - r_kj[m] * rkj_2 * cos_theta);
3476             f_j[m] = -f_i[m] - f_k[m];
3477             f[ai][m] += f_i[m];
3478             f[aj][m] += f_j[m];
3479             f[ak][m] += f_k[m];
3480         }
3481
3482         if (computeVirial(flavor))
3483         {
3484             rvec_inc(fshift[t1], f_i);
3485             rvec_inc(fshift[CENTRAL], f_j);
3486             rvec_inc(fshift[t2], f_k); /* 9 */
3487         }
3488         /* 163 TOTAL    */
3489     }
3490     return vtot;
3491 }
3492
3493 template<BondedKernelFlavor flavor>
3494 real cross_bond_bond(int             nbonds,
3495                      const t_iatom   forceatoms[],
3496                      const t_iparams forceparams[],
3497                      const rvec      x[],
3498                      rvec4           f[],
3499                      rvec            fshift[],
3500                      const t_pbc*    pbc,
3501                      real gmx_unused lambda,
3502                      real gmx_unused* dvdlambda,
3503                      const t_mdatoms gmx_unused* md,
3504                      t_fcdata gmx_unused* fcd,
3505                      int gmx_unused* global_atom_index)
3506 {
3507     /* Potential from Lawrence and Skimmer, Chem. Phys. Lett. 372 (2003)
3508      * pp. 842-847
3509      */
3510     int  i, ai, aj, ak, type, m, t1, t2;
3511     rvec r_ij, r_kj;
3512     real vtot, vrr, s1, s2, r1, r2, r1e, r2e, krr;
3513     rvec f_i, f_j, f_k;
3514
3515     vtot = 0.0;
3516     for (i = 0; (i < nbonds);)
3517     {
3518         type = forceatoms[i++];
3519         ai   = forceatoms[i++];
3520         aj   = forceatoms[i++];
3521         ak   = forceatoms[i++];
3522         r1e  = forceparams[type].cross_bb.r1e;
3523         r2e  = forceparams[type].cross_bb.r2e;
3524         krr  = forceparams[type].cross_bb.krr;
3525
3526         /* Compute distance vectors ... */
3527         t1 = pbc_rvec_sub(pbc, x[ai], x[aj], r_ij);
3528         t2 = pbc_rvec_sub(pbc, x[ak], x[aj], r_kj);
3529
3530         /* ... and their lengths */
3531         r1 = norm(r_ij);
3532         r2 = norm(r_kj);
3533
3534         /* Deviations from ideality */
3535         s1 = r1 - r1e;
3536         s2 = r2 - r2e;
3537
3538         /* Energy (can be negative!) */
3539         vrr = krr * s1 * s2;
3540         vtot += vrr;
3541
3542         /* Forces */
3543         svmul(-krr * s2 / r1, r_ij, f_i);
3544         svmul(-krr * s1 / r2, r_kj, f_k);
3545
3546         for (m = 0; (m < DIM); m++) /*  12      */
3547         {
3548             f_j[m] = -f_i[m] - f_k[m];
3549             f[ai][m] += f_i[m];
3550             f[aj][m] += f_j[m];
3551             f[ak][m] += f_k[m];
3552         }
3553
3554         if (computeVirial(flavor))
3555         {
3556             rvec_inc(fshift[t1], f_i);
3557             rvec_inc(fshift[CENTRAL], f_j);
3558             rvec_inc(fshift[t2], f_k); /* 9 */
3559         }
3560         /* 163 TOTAL    */
3561     }
3562     return vtot;
3563 }
3564
3565 template<BondedKernelFlavor flavor>
3566 real cross_bond_angle(int             nbonds,
3567                       const t_iatom   forceatoms[],
3568                       const t_iparams forceparams[],
3569                       const rvec      x[],
3570                       rvec4           f[],
3571                       rvec            fshift[],
3572                       const t_pbc*    pbc,
3573                       real gmx_unused lambda,
3574                       real gmx_unused* dvdlambda,
3575                       const t_mdatoms gmx_unused* md,
3576                       t_fcdata gmx_unused* fcd,
3577                       int gmx_unused* global_atom_index)
3578 {
3579     /* Potential from Lawrence and Skimmer, Chem. Phys. Lett. 372 (2003)
3580      * pp. 842-847
3581      */
3582     int  i, ai, aj, ak, type, m, t1, t2;
3583     rvec r_ij, r_kj, r_ik;
3584     real vtot, vrt, s1, s2, s3, r1, r2, r3, r1e, r2e, r3e, krt, k1, k2, k3;
3585     rvec f_i, f_j, f_k;
3586
3587     vtot = 0.0;
3588     for (i = 0; (i < nbonds);)
3589     {
3590         type = forceatoms[i++];
3591         ai   = forceatoms[i++];
3592         aj   = forceatoms[i++];
3593         ak   = forceatoms[i++];
3594         r1e  = forceparams[type].cross_ba.r1e;
3595         r2e  = forceparams[type].cross_ba.r2e;
3596         r3e  = forceparams[type].cross_ba.r3e;
3597         krt  = forceparams[type].cross_ba.krt;
3598
3599         /* Compute distance vectors ... */
3600         t1 = pbc_rvec_sub(pbc, x[ai], x[aj], r_ij);
3601         t2 = pbc_rvec_sub(pbc, x[ak], x[aj], r_kj);
3602         pbc_rvec_sub(pbc, x[ai], x[ak], r_ik);
3603
3604         /* ... and their lengths */
3605         r1 = norm(r_ij);
3606         r2 = norm(r_kj);
3607         r3 = norm(r_ik);
3608
3609         /* Deviations from ideality */
3610         s1 = r1 - r1e;
3611         s2 = r2 - r2e;
3612         s3 = r3 - r3e;
3613
3614         /* Energy (can be negative!) */
3615         vrt = krt * s3 * (s1 + s2);
3616         vtot += vrt;
3617
3618         /* Forces */
3619         k1 = -krt * (s3 / r1);
3620         k2 = -krt * (s3 / r2);
3621         k3 = -krt * (s1 + s2) / r3;
3622         for (m = 0; (m < DIM); m++)
3623         {
3624             f_i[m] = k1 * r_ij[m] + k3 * r_ik[m];
3625             f_k[m] = k2 * r_kj[m] - k3 * r_ik[m];
3626             f_j[m] = -f_i[m] - f_k[m];
3627         }
3628
3629         for (m = 0; (m < DIM); m++) /*  12      */
3630         {
3631             f[ai][m] += f_i[m];
3632             f[aj][m] += f_j[m];
3633             f[ak][m] += f_k[m];
3634         }
3635
3636         if (computeVirial(flavor))
3637         {
3638             rvec_inc(fshift[t1], f_i);
3639             rvec_inc(fshift[CENTRAL], f_j);
3640             rvec_inc(fshift[t2], f_k); /* 9 */
3641         }
3642         /* 163 TOTAL    */
3643     }
3644     return vtot;
3645 }
3646
3647 /*! \brief Computes the potential and force for a tabulated potential */
3648 real bonded_tab(const char*          type,
3649                 int                  table_nr,
3650                 const bondedtable_t* table,
3651                 real                 kA,
3652                 real                 kB,
3653                 real                 r,
3654                 real                 lambda,
3655                 real*                V,
3656                 real*                F)
3657 {
3658     real k, tabscale, rt, eps, eps2, Yt, Ft, Geps, Heps2, Fp, VV, FF;
3659     int  n0, nnn;
3660     real dvdlambda;
3661
3662     k = (1.0 - lambda) * kA + lambda * kB;
3663
3664     tabscale          = table->scale;
3665     const real* VFtab = table->data.data();
3666
3667     rt = r * tabscale;
3668     n0 = static_cast<int>(rt);
3669     if (n0 >= table->n)
3670     {
3671         gmx_fatal(FARGS,
3672                   "A tabulated %s interaction table number %d is out of the table range: r %f, "
3673                   "between table indices %d and %d, table length %d",
3674                   type, table_nr, r, n0, n0 + 1, table->n);
3675     }
3676     eps   = rt - n0;
3677     eps2  = eps * eps;
3678     nnn   = 4 * n0;
3679     Yt    = VFtab[nnn];
3680     Ft    = VFtab[nnn + 1];
3681     Geps  = VFtab[nnn + 2] * eps;
3682     Heps2 = VFtab[nnn + 3] * eps2;
3683     Fp    = Ft + Geps + Heps2;
3684     VV    = Yt + Fp * eps;
3685     FF    = Fp + Geps + 2.0 * Heps2;
3686
3687     *F        = -k * FF * tabscale;
3688     *V        = k * VV;
3689     dvdlambda = (kB - kA) * VV;
3690
3691     return dvdlambda;
3692
3693     /* That was 22 flops */
3694 }
3695
3696 template<BondedKernelFlavor flavor>
3697 real tab_bonds(int             nbonds,
3698                const t_iatom   forceatoms[],
3699                const t_iparams forceparams[],
3700                const rvec      x[],
3701                rvec4           f[],
3702                rvec            fshift[],
3703                const t_pbc*    pbc,
3704                real            lambda,
3705                real*           dvdlambda,
3706                const t_mdatoms gmx_unused* md,
3707                t_fcdata*                   fcd,
3708                int gmx_unused* global_atom_index)
3709 {
3710     int  i, ki, ai, aj, type, table;
3711     real dr, dr2, fbond, vbond, vtot;
3712     rvec dx;
3713
3714     vtot = 0.0;
3715     for (i = 0; (i < nbonds);)
3716     {
3717         type = forceatoms[i++];
3718         ai   = forceatoms[i++];
3719         aj   = forceatoms[i++];
3720
3721         ki  = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /*   3      */
3722         dr2 = iprod(dx, dx);                       /*   5               */
3723         dr  = dr2 * gmx::invsqrt(dr2);             /*  10               */
3724
3725         table = forceparams[type].tab.table;
3726
3727         *dvdlambda += bonded_tab("bond", table, &fcd->bondtab[table], forceparams[type].tab.kA,
3728                                  forceparams[type].tab.kB, dr, lambda, &vbond, &fbond); /*  22 */
3729
3730         if (dr2 == 0.0)
3731         {
3732             continue;
3733         }
3734
3735
3736         vtot += vbond;              /* 1*/
3737         fbond *= gmx::invsqrt(dr2); /*   6              */
3738
3739         spreadBondForces<flavor>(fbond, dx, ai, aj, f, ki, fshift); /* 15 */
3740     }                                                               /* 62 TOTAL */
3741     return vtot;
3742 }
3743
3744 template<BondedKernelFlavor flavor>
3745 real tab_angles(int             nbonds,
3746                 const t_iatom   forceatoms[],
3747                 const t_iparams forceparams[],
3748                 const rvec      x[],
3749                 rvec4           f[],
3750                 rvec            fshift[],
3751                 const t_pbc*    pbc,
3752                 real            lambda,
3753                 real*           dvdlambda,
3754                 const t_mdatoms gmx_unused* md,
3755                 t_fcdata*                   fcd,
3756                 int gmx_unused* global_atom_index)
3757 {
3758     int  i, ai, aj, ak, t1, t2, type, table;
3759     rvec r_ij, r_kj;
3760     real cos_theta, cos_theta2, theta, dVdt, va, vtot;
3761
3762     vtot = 0.0;
3763     for (i = 0; (i < nbonds);)
3764     {
3765         type = forceatoms[i++];
3766         ai   = forceatoms[i++];
3767         aj   = forceatoms[i++];
3768         ak   = forceatoms[i++];
3769
3770         theta = bond_angle(x[ai], x[aj], x[ak], pbc, r_ij, r_kj, &cos_theta, &t1, &t2); /*  41 */
3771
3772         table = forceparams[type].tab.table;
3773
3774         *dvdlambda += bonded_tab("angle", table, &fcd->angletab[table], forceparams[type].tab.kA,
3775                                  forceparams[type].tab.kB, theta, lambda, &va, &dVdt); /*  22  */
3776         vtot += va;
3777
3778         cos_theta2 = gmx::square(cos_theta); /*   1             */
3779         if (cos_theta2 < 1)
3780         {
3781             int  m;
3782             real st, sth;
3783             real cik, cii, ckk;
3784             real nrkj2, nrij2;
3785             rvec f_i, f_j, f_k;
3786
3787             st    = dVdt * gmx::invsqrt(1 - cos_theta2); /*  12         */
3788             sth   = st * cos_theta;                      /*   1         */
3789             nrkj2 = iprod(r_kj, r_kj);                   /*   5         */
3790             nrij2 = iprod(r_ij, r_ij);
3791
3792             cik = st * gmx::invsqrt(nrkj2 * nrij2); /*  12              */
3793             cii = sth / nrij2;                      /*  10              */
3794             ckk = sth / nrkj2;                      /*  10              */
3795
3796             for (m = 0; (m < DIM); m++) /*  39          */
3797             {
3798                 f_i[m] = -(cik * r_kj[m] - cii * r_ij[m]);
3799                 f_k[m] = -(cik * r_ij[m] - ckk * r_kj[m]);
3800                 f_j[m] = -f_i[m] - f_k[m];
3801                 f[ai][m] += f_i[m];
3802                 f[aj][m] += f_j[m];
3803                 f[ak][m] += f_k[m];
3804             }
3805
3806             if (computeVirial(flavor))
3807             {
3808                 rvec_inc(fshift[t1], f_i);
3809                 rvec_inc(fshift[CENTRAL], f_j);
3810                 rvec_inc(fshift[t2], f_k);
3811             }
3812         } /* 169 TOTAL  */
3813     }
3814     return vtot;
3815 }
3816
3817 template<BondedKernelFlavor flavor>
3818 real tab_dihs(int             nbonds,
3819               const t_iatom   forceatoms[],
3820               const t_iparams forceparams[],
3821               const rvec      x[],
3822               rvec4           f[],
3823               rvec            fshift[],
3824               const t_pbc*    pbc,
3825               real            lambda,
3826               real*           dvdlambda,
3827               const t_mdatoms gmx_unused* md,
3828               t_fcdata*                   fcd,
3829               int gmx_unused* global_atom_index)
3830 {
3831     int  i, type, ai, aj, ak, al, table;
3832     int  t1, t2, t3;
3833     rvec r_ij, r_kj, r_kl, m, n;
3834     real phi, ddphi, vpd, vtot;
3835
3836     vtot = 0.0;
3837     for (i = 0; (i < nbonds);)
3838     {
3839         type = forceatoms[i++];
3840         ai   = forceatoms[i++];
3841         aj   = forceatoms[i++];
3842         ak   = forceatoms[i++];
3843         al   = forceatoms[i++];
3844
3845         phi = dih_angle(x[ai], x[aj], x[ak], x[al], pbc, r_ij, r_kj, r_kl, m, n, &t1, &t2, &t3); /*  84  */
3846
3847         table = forceparams[type].tab.table;
3848
3849         /* Hopefully phi+M_PI never results in values < 0 */
3850         *dvdlambda += bonded_tab("dihedral", table, &fcd->dihtab[table], forceparams[type].tab.kA,
3851                                  forceparams[type].tab.kB, phi + M_PI, lambda, &vpd, &ddphi);
3852
3853         vtot += vpd;
3854         do_dih_fup<flavor>(ai, aj, ak, al, -ddphi, r_ij, r_kj, r_kl, m, n, f, fshift, pbc, x, t1,
3855                            t2, t3); /* 112      */
3856
3857     } /* 227 TOTAL  */
3858
3859     return vtot;
3860 }
3861
3862 struct BondedInteractions
3863 {
3864     BondedFunction function;
3865     int            nrnbIndex;
3866 };
3867
3868 // Bug in old clang versions prevents constexpr. constexpr is needed for MSVC.
3869 #if defined(__clang__) && __clang_major__ < 6
3870 #    define CONSTEXPR_EXCL_OLD_CLANG const
3871 #else
3872 #    define CONSTEXPR_EXCL_OLD_CLANG constexpr
3873 #endif
3874
3875 /*! \brief Lookup table of bonded interaction functions
3876  *
3877  * This must have as many entries as interaction_function in ifunc.cpp */
3878 template<BondedKernelFlavor flavor>
3879 CONSTEXPR_EXCL_OLD_CLANG std::array<BondedInteractions, F_NRE> c_bondedInteractionFunctions = {
3880     BondedInteractions{ bonds<flavor>, eNR_BONDS },                       // F_BONDS
3881     BondedInteractions{ g96bonds<flavor>, eNR_BONDS },                    // F_G96BONDS
3882     BondedInteractions{ morse_bonds<flavor>, eNR_MORSE },                 // F_MORSE
3883     BondedInteractions{ cubic_bonds<flavor>, eNR_CUBICBONDS },            // F_CUBICBONDS
3884     BondedInteractions{ unimplemented, -1 },                              // F_CONNBONDS
3885     BondedInteractions{ bonds<flavor>, eNR_BONDS },                       // F_HARMONIC
3886     BondedInteractions{ FENE_bonds<flavor>, eNR_FENEBONDS },              // F_FENEBONDS
3887     BondedInteractions{ tab_bonds<flavor>, eNR_TABBONDS },                // F_TABBONDS
3888     BondedInteractions{ tab_bonds<flavor>, eNR_TABBONDS },                // F_TABBONDSNC
3889     BondedInteractions{ restraint_bonds<flavor>, eNR_RESTRBONDS },        // F_RESTRBONDS
3890     BondedInteractions{ angles<flavor>, eNR_ANGLES },                     // F_ANGLES
3891     BondedInteractions{ g96angles<flavor>, eNR_ANGLES },                  // F_G96ANGLES
3892     BondedInteractions{ restrangles<flavor>, eNR_ANGLES },                // F_RESTRANGLES
3893     BondedInteractions{ linear_angles<flavor>, eNR_ANGLES },              // F_LINEAR_ANGLES
3894     BondedInteractions{ cross_bond_bond<flavor>, eNR_CROSS_BOND_BOND },   // F_CROSS_BOND_BONDS
3895     BondedInteractions{ cross_bond_angle<flavor>, eNR_CROSS_BOND_ANGLE }, // F_CROSS_BOND_ANGLES
3896     BondedInteractions{ urey_bradley<flavor>, eNR_UREY_BRADLEY },         // F_UREY_BRADLEY
3897     BondedInteractions{ quartic_angles<flavor>, eNR_QANGLES },            // F_QUARTIC_ANGLES
3898     BondedInteractions{ tab_angles<flavor>, eNR_TABANGLES },              // F_TABANGLES
3899     BondedInteractions{ pdihs<flavor>, eNR_PROPER },                      // F_PDIHS
3900     BondedInteractions{ rbdihs<flavor>, eNR_RB },                         // F_RBDIHS
3901     BondedInteractions{ restrdihs<flavor>, eNR_PROPER },                  // F_RESTRDIHS
3902     BondedInteractions{ cbtdihs<flavor>, eNR_RB },                        // F_CBTDIHS
3903     BondedInteractions{ rbdihs<flavor>, eNR_FOURDIH },                    // F_FOURDIHS
3904     BondedInteractions{ idihs<flavor>, eNR_IMPROPER },                    // F_IDIHS
3905     BondedInteractions{ pdihs<flavor>, eNR_IMPROPER },                    // F_PIDIHS
3906     BondedInteractions{ tab_dihs<flavor>, eNR_TABDIHS },                  // F_TABDIHS
3907     BondedInteractions{ unimplemented, eNR_CMAP },                        // F_CMAP
3908     BondedInteractions{ unimplemented, -1 },                              // F_GB12_NOLONGERUSED
3909     BondedInteractions{ unimplemented, -1 },                              // F_GB13_NOLONGERUSED
3910     BondedInteractions{ unimplemented, -1 },                              // F_GB14_NOLONGERUSED
3911     BondedInteractions{ unimplemented, -1 },                              // F_GBPOL_NOLONGERUSED
3912     BondedInteractions{ unimplemented, -1 },                       // F_NPSOLVATION_NOLONGERUSED
3913     BondedInteractions{ unimplemented, eNR_NB14 },                 // F_LJ14
3914     BondedInteractions{ unimplemented, -1 },                       // F_COUL14
3915     BondedInteractions{ unimplemented, eNR_NB14 },                 // F_LJC14_Q
3916     BondedInteractions{ unimplemented, eNR_NB14 },                 // F_LJC_PAIRS_NB
3917     BondedInteractions{ unimplemented, -1 },                       // F_LJ
3918     BondedInteractions{ unimplemented, -1 },                       // F_BHAM
3919     BondedInteractions{ unimplemented, -1 },                       // F_LJ_LR_NOLONGERUSED
3920     BondedInteractions{ unimplemented, -1 },                       // F_BHAM_LR_NOLONGERUSED
3921     BondedInteractions{ unimplemented, -1 },                       // F_DISPCORR
3922     BondedInteractions{ unimplemented, -1 },                       // F_COUL_SR
3923     BondedInteractions{ unimplemented, -1 },                       // F_COUL_LR_NOLONGERUSED
3924     BondedInteractions{ unimplemented, -1 },                       // F_RF_EXCL
3925     BondedInteractions{ unimplemented, -1 },                       // F_COUL_RECIP
3926     BondedInteractions{ unimplemented, -1 },                       // F_LJ_RECIP
3927     BondedInteractions{ unimplemented, -1 },                       // F_DPD
3928     BondedInteractions{ polarize<flavor>, eNR_POLARIZE },          // F_POLARIZATION
3929     BondedInteractions{ water_pol<flavor>, eNR_WPOL },             // F_WATER_POL
3930     BondedInteractions{ thole_pol<flavor>, eNR_THOLE },            // F_THOLE_POL
3931     BondedInteractions{ anharm_polarize<flavor>, eNR_ANHARM_POL }, // F_ANHARM_POL
3932     BondedInteractions{ unimplemented, -1 },                       // F_POSRES
3933     BondedInteractions{ unimplemented, -1 },                       // F_FBPOSRES
3934     BondedInteractions{ ta_disres, eNR_DISRES },                   // F_DISRES
3935     BondedInteractions{ unimplemented, -1 },                       // F_DISRESVIOL
3936     BondedInteractions{ orires, eNR_ORIRES },                      // F_ORIRES
3937     BondedInteractions{ unimplemented, -1 },                       // F_ORIRESDEV
3938     BondedInteractions{ angres<flavor>, eNR_ANGRES },              // F_ANGRES
3939     BondedInteractions{ angresz<flavor>, eNR_ANGRESZ },            // F_ANGRESZ
3940     BondedInteractions{ dihres<flavor>, eNR_DIHRES },              // F_DIHRES
3941     BondedInteractions{ unimplemented, -1 },                       // F_DIHRESVIOL
3942     BondedInteractions{ unimplemented, -1 },                       // F_CONSTR
3943     BondedInteractions{ unimplemented, -1 },                       // F_CONSTRNC
3944     BondedInteractions{ unimplemented, -1 },                       // F_SETTLE
3945     BondedInteractions{ unimplemented, -1 },                       // F_VSITE2
3946     BondedInteractions{ unimplemented, -1 },                       // F_VSITE3
3947     BondedInteractions{ unimplemented, -1 },                       // F_VSITE3FD
3948     BondedInteractions{ unimplemented, -1 },                       // F_VSITE3FAD
3949     BondedInteractions{ unimplemented, -1 },                       // F_VSITE3OUT
3950     BondedInteractions{ unimplemented, -1 },                       // F_VSITE4FD
3951     BondedInteractions{ unimplemented, -1 },                       // F_VSITE4FDN
3952     BondedInteractions{ unimplemented, -1 },                       // F_VSITEN
3953     BondedInteractions{ unimplemented, -1 },                       // F_COM_PULL
3954     BondedInteractions{ unimplemented, -1 },                       // F_DENSITYFITTING
3955     BondedInteractions{ unimplemented, -1 },                       // F_EQM
3956     BondedInteractions{ unimplemented, -1 },                       // F_EPOT
3957     BondedInteractions{ unimplemented, -1 },                       // F_EKIN
3958     BondedInteractions{ unimplemented, -1 },                       // F_ETOT
3959     BondedInteractions{ unimplemented, -1 },                       // F_ECONSERVED
3960     BondedInteractions{ unimplemented, -1 },                       // F_TEMP
3961     BondedInteractions{ unimplemented, -1 },                       // F_VTEMP_NOLONGERUSED
3962     BondedInteractions{ unimplemented, -1 },                       // F_PDISPCORR
3963     BondedInteractions{ unimplemented, -1 },                       // F_PRES
3964     BondedInteractions{ unimplemented, -1 },                       // F_DVDL_CONSTR
3965     BondedInteractions{ unimplemented, -1 },                       // F_DVDL
3966     BondedInteractions{ unimplemented, -1 },                       // F_DKDL
3967     BondedInteractions{ unimplemented, -1 },                       // F_DVDL_COUL
3968     BondedInteractions{ unimplemented, -1 },                       // F_DVDL_VDW
3969     BondedInteractions{ unimplemented, -1 },                       // F_DVDL_BONDED
3970     BondedInteractions{ unimplemented, -1 },                       // F_DVDL_RESTRAINT
3971     BondedInteractions{ unimplemented, -1 },                       // F_DVDL_TEMPERATURE
3972 };
3973
3974 /*! \brief List of instantiated BondedInteractions list */
3975 CONSTEXPR_EXCL_OLD_CLANG
3976 gmx::EnumerationArray<BondedKernelFlavor, std::array<BondedInteractions, F_NRE>> c_bondedInteractionFunctionsPerFlavor = {
3977     c_bondedInteractionFunctions<BondedKernelFlavor::ForcesSimdWhenAvailable>,
3978     c_bondedInteractionFunctions<BondedKernelFlavor::ForcesNoSimd>,
3979     c_bondedInteractionFunctions<BondedKernelFlavor::ForcesAndVirialAndEnergy>,
3980     c_bondedInteractionFunctions<BondedKernelFlavor::ForcesAndEnergy>
3981 };
3982
3983 //! \endcond
3984
3985 } // namespace
3986
3987 real calculateSimpleBond(const int           ftype,
3988                          const int           numForceatoms,
3989                          const t_iatom       forceatoms[],
3990                          const t_iparams     forceparams[],
3991                          const rvec          x[],
3992                          rvec4               f[],
3993                          rvec                fshift[],
3994                          const struct t_pbc* pbc,
3995                          const real          lambda,
3996                          real*               dvdlambda,
3997                          const t_mdatoms*    md,
3998                          t_fcdata*           fcd,
3999                          int gmx_unused*          global_atom_index,
4000                          const BondedKernelFlavor bondedKernelFlavor)
4001 {
4002     const BondedInteractions& bonded = c_bondedInteractionFunctionsPerFlavor[bondedKernelFlavor][ftype];
4003
4004     real v = bonded.function(numForceatoms, forceatoms, forceparams, x, f, fshift, pbc, lambda,
4005                              dvdlambda, md, fcd, global_atom_index);
4006
4007     return v;
4008 }
4009
4010 int nrnbIndex(int ftype)
4011 {
4012     return c_bondedInteractionFunctions<BondedKernelFlavor::ForcesAndVirialAndEnergy>[ftype].nrnbIndex;
4013 }