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