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