api/nblib/listed_forces/kernels.hpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2020, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 /*! \inpublicapi \file
  36  * \brief
  37  * Implements kernels for nblib supported bondtypes
  38  *
  39  * \author Victor Holanda <victor.holanda@cscs.ch>
  40  * \author Joe Jordan <ejjordan@kth.se>
  41  * \author Prashanth Kanduri <kanduri@cscs.ch>
  42  * \author Sebastian Keller <keller@cscs.ch>
  43  * \author Artem Zhmurov <zhmurov@gmail.com>
  44  */
  45 #ifndef NBLIB_LISTEDFORCES_KERNELS_HPP
  46 #define NBLIB_LISTEDFORCES_KERNELS_HPP
  47
  48 #include <tuple>
  49
  50 #include "gromacs/math/functions.h"
  51 #include "gromacs/math/vectypes.h"
  52 #include "nblib/listed_forces/bondtypes.h"
  53
  54 #define NBLIB_ALWAYS_INLINE __attribute((always_inline))
  55
  56 namespace nblib
  57 {
  58
  59 /*! \brief kernel to calculate the scalar part for simple harmonic bond forces
  60  *         for lambda = 0
  61  *
  62  * \param k spring constant
  63  * \param x0 equilibrium distance
  64  * \param x  input bond length
  65  *
  66  * \return tuple<force, potential energy>
  67  */
  68 template <class T>
  69 inline std::tuple<T, T> harmonicScalarForce(T k, T x0, T x)
  70 {
  71     T dx  = x - x0;
  72     T dx2 = dx * dx;
  73
  74     T force = -k * dx;
  75     T epot = 0.5 * k * dx2;
  76
  77     return std::make_tuple(force, epot);
  78
  79     /* That was 6 flops */
  80 }
  81
  82 /*! \brief kernel to calculate the scalar part for simple harmonic bond forces
  83  *         for non-zero lambda to interpolate between A and B states
  84  *
  85  * \param kA spring constant state A
  86  * \param kB spring constant state B
  87  * \param xA equilibrium distance state A
  88  * \param xB equilibrium distance state B
  89  * \param x  input bond length
  90  * \param lambda interpolation factor between A and B state
  91  *
  92  * \return tuple<force, potential energy, lambda-interpolated energy>
  93  */
  94 template <class T>
  95 inline std::tuple<T, T, T> harmonicScalarForce(T kA, T kB, T xA, T xB, T x, T lambda)
  96 {
  97     // code unchanged relative to Gromacs
  98
  99     T L1 = 1.0 - lambda;
 100     T kk = L1 * kA + lambda * kB;
 101     T x0 = L1 * xA + lambda * xB;
 102
 103     T dx  = x - x0;
 104     T dx2 = dx * dx;
 105
 106     T force     = -kk * dx;
 107     T epot      = 0.5 * kk * dx2;
 108     T dvdlambda = 0.5 * (kB - kA) * dx2 + (xA - xB) * kk * dx;
 109
 110     return std::make_tuple(force, epot, dvdlambda);
 111
 112     /* That was 19 flops */
 113 }
 114
 115 //! abstraction layer for different 2-center bonds
 116 template <class T>
 117 inline auto bondKernel(T dr, const HarmonicBondType& bond)
 118 {
 119     return harmonicScalarForce(bond.forceConstant(), bond.equilConstant(), dr);
 120 }
 121
 122
 123 /*! \brief kernel to calculate the scalar part for the forth power pontential bond forces
 124  *         for lambda = 0
 125  *
 126  * \param k spring constant
 127  * \param x0 squared equilibrium distance
 128  * \param x  squared input bond length
 129  *
 130  * \return tuple<force, potential energy>
 131  */
 132 template <class T>
 133 inline std::tuple<T, T> g96ScalarForce(T k, T x0, T x)
 134 {
 135     T dx  = x - x0;
 136     T dx2 = dx * dx;
 137
 138     T force = -k * dx;
 139     T epot = 0.5 * k * dx2;
 140
 141     return std::make_tuple(force, epot);
 142
 143     /* That was 6 flops */
 144 }
 145
 146 /*! \brief kernel to calculate the scalar part for forth power pontential bond forces
 147  *         for non-zero lambda to interpolate between A and B states
 148  *
 149  * \param kA spring constant state A
 150  * \param kB spring constant state B
 151  * \param xA squared equilibrium distance state A
 152  * \param xB squared equilibrium distance state B
 153  * \param x  squared input bond length
 154  * \param lambda interpolation factor between A and B state
 155  *
 156  * \return tuple<force, potential energy, lambda-interpolated energy>
 157  */
 158 template <class T>
 159 inline std::tuple<T, T, T> g96ScalarForce(T kA, T kB, T xA, T xB, T x, T lambda)
 160 {
 161     T L1 = 1.0 - lambda;
 162     T kk = L1 * kA + lambda * kB;
 163     T x0 = L1 * xA + lambda * xB;
 164
 165     T dx  = x - x0;
 166     T dx2 = dx * dx;
 167
 168     T force = -kk * dx;
 169     T epot = 0.5 * kk * dx2;
 170     // TODO: Check if this is 1/2 or 1/4
 171     T dvdlambda = 0.5 * (kB - kA) * dx2 + (xA - xB) * kk * dx;
 172
 173     return std::make_tuple(force, epot, dvdlambda);
 174
 175     /* That was 21 flops */
 176 }
 177
 178 //! Abstraction layer for different 2-center bonds. Fourth power case
 179 template <class T>
 180 inline auto bondKernel(T dr, const G96BondType& bond)
 181 {
 182     auto [force, ePot] = g96ScalarForce(bond.forceConstant(), bond.equilConstant(), dr*dr);
 183     force *= dr;
 184     ePot  *= 0.5;
 185     return std::make_tuple(force, ePot);
 186 }
 187
 188
 189 /*! \brief kernel to calculate the scalar part for the morse pontential bond forces
 190  *         for lambda = 0
 191  *
 192  * \param k force constant
 193  * \param beta beta exponent
 194  * \param x0 equilibrium distance
 195  * \param x  input bond length
 196  *
 197  * \return tuple<force, potential energy>
 198  */
 199 template <class T>
 200 inline std::tuple<T, T> morseScalarForce(T k, T beta, T x0, T x)
 201 {
 202     T exponent = std::exp(-beta * (x - x0));      /* 12 */
 203     T omexp = 1.0 - exponent;                     /*  1 */
 204     T kexp = k * omexp;                           /*  1 */
 205
 206     T epot = kexp * omexp;                        /*  1 */
 207     T force = -2.0 * beta * exponent * kexp;      /*  4 */
 208
 209     return std::make_tuple(force, epot);
 210
 211     /* That was 20 flops */
 212 }
 213
 214 /*! \brief kernel to calculate the scalar part for morse potential bond forces
 215  *         for non-zero lambda to interpolate between A and B states
 216  *
 217  * \param kA force constant state A
 218  * \param kB force constant state B
 219  * \param betaA beta exponent state A
 220  * \param betaB beta exponent state B
 221  * \param xA equilibrium distance state A
 222  * \param xB equilibrium distance state B
 223  * \param x  input bond length
 224  * \param lambda interpolation factor between A and B state
 225  *
 226  * \return tuple<force, potential energy, lambda-interpolated energy>
 227  */
 228 template <class T>
 229 inline std::tuple<T, T, T> morseScalarForce(T kA, T kB, T betaA, T betaB, T xA, T xB, T x, T lambda)
 230 {
 231     T L1 = 1.0 - lambda;                          /* 1 */
 232     T x0 = L1 * xA + lambda * xB;                 /* 3 */
 233     T beta = L1 * betaA + lambda * betaB;         /* 3 */
 234     T k = L1 * kA + lambda * kB;                  /* 3 */
 235
 236     T exponent = std::exp(-beta * (x - x0));      /* 12 */
 237     T omexp = 1.0 - exponent;                     /*  1 */
 238     T kexp = k * omexp;                           /*  1 */
 239
 240     T epot = kexp * omexp;                        /*  1 */
 241     T force = -2.0 * beta * exponent * kexp;      /*  4 */
 242
 243     T dvdlambda = (kB - kA) * omexp * omexp
 244                     - (2.0 - 2.0 * omexp) * omexp * k
 245                     * ((xB - xA) * beta - (betaB - betaA) * (x - x0)); /* 15 */
 246
 247     return std::make_tuple(force, epot, dvdlambda);
 248
 249     /* That was 44 flops */
 250 }
 251
 252 //! Abstraction layer for different 2-center bonds. Morse case
 253 template <class T>
 254 inline auto bondKernel(T dr, const MorseBondType& bond)
 255 {
 256     return morseScalarForce(bond.forceConstant(), bond.exponent(), bond.equilDistance(), dr);
 257 }
 258
 259 /*! \brief kernel to calculate the scalar part for the 1-4 LJ non-bonded forces
 260  *
 261  * \param c6  C6 parameter of LJ potential
 262  * \param c12 C12 parameter of LJ potential
 263  * \param r   distance between the atoms
 264  *
 265  * \return tuple<force, potential energy>
 266  */
 267 template <class T>
 268 inline std::tuple<T, T> pairLJScalarForce(C6 c6, C12 c12, T r)
 269 {
 270     T rinv  = 1./r;                           /* 1 */
 271     T rinv2 = rinv * rinv;                    /* 1 */
 272     T rinv6 = rinv2 * rinv2 * rinv2;          /* 2 */
 273
 274     T epot  = rinv6*(c12*rinv6 - c6);         /* 3 */
 275
 276     T c6_     = 6.*c6;                        /* 1 */
 277     T c12_    = 12.*c12;                      /* 1 */
 278
 279     T force = rinv6*(c12_*rinv6 - c6_)*rinv;  /* 4 */
 280
 281     return std::make_tuple(force, epot);
 282
 283     /* That was 13 flops */
 284 }
 285
 286 //! Abstraction layer for different 2-center bonds. 1-4 LJ pair interactions case
 287 template <class T>
 288 inline auto bondKernel(T dr, const PairLJType& bond)
 289 {
 290     return pairLJScalarForce(bond.c6(), bond.c12(), dr);
 291 }
 292
 293
 294 /*! \brief kernel to calculate the scalar part for the FENE pontential bond forces
 295  *         for lambda = 0
 296  *
 297  * \param k spring constant
 298  * \param x0 equilibrium distance
 299  * \param x  input bond length
 300  *
 301  * \return tuple<force, potential energy>
 302  */
 303 template <class T>
 304 inline std::tuple<T, T> FENEScalarForce(T k, T x0, T x)
 305 {
 306     T x02 = x0 * x0;
 307     T x2 = x * x;
 308
 309     T omx2_ox02 = 1.0 - (x2 / x02);
 310
 311     T epot = -0.5 * k * x02 * std::log(omx2_ox02);
 312     T force = -k / omx2_ox02;
 313
 314     return std::make_tuple(force, epot);
 315
 316     /* That was 24 flops */
 317 }
 318
 319 // TODO: Implement the free energy version of FENE (finitely extensible nonlinear elastic) bond types
 320
 321 //! Abstraction layer for different 2-center bonds. FENE case
 322 template <class T>
 323 inline auto bondKernel(T dr, const FENEBondType& bond)
 324 {
 325     auto [force, ePot] = FENEScalarForce(bond.forceConstant(), bond.equilConstant(), dr);
 326     force *= dr;
 327     return std::make_tuple(force, ePot);
 328 }
 329
 330
 331 /*! \brief kernel to calculate the scalar part for cubic potential bond forces
 332  *         for lambda = 0
 333  *
 334  * \param kc cubic spring constant
 335  * \param kq quadratic spring constant
 336  * \param x0 equilibrium distance
 337  * \param x  input bond length
 338  *
 339  * \return tuple<force, potential energy>
 340  */
 341 template <class T>
 342 inline std::tuple<T, T> cubicScalarForce(T kc, T kq, T x0, T x)
 343 {
 344     T dx = x - x0;
 345     //T dx2 = dx * dx;
 346
 347     T kdist  = kq * dx;
 348     T kdist2 = kdist * dx;
 349
 350     T epot = kdist2 + (kc * kdist2 * dx);
 351     T force = -((2.0 * kdist) + (3.0 * kdist2 * kc));
 352
 353     return std::make_tuple(force, epot);
 354
 355     /* That was 16 flops */
 356 }
 357
 358 // TODO: Implement the free energy version of Cubic bond types
 359
 360 template <class T>
 361 inline auto bondKernel(T dr, const CubicBondType& bond)
 362 {
 363     return cubicScalarForce(bond.cubicForceConstant(), bond.quadraticForceConstant(), bond.equilDistance(), dr);
 364 }
 365
 366
 367 /*! \brief kernel to calculate the scalar part for half attractive potential bond forces
 368  *         for lambda = 0
 369  *
 370  * \param k spring constant
 371  * \param x0 equilibrium distance
 372  * \param x  input bond length
 373  *
 374  * \return tuple<force, potential energy>
 375  */
 376 template <class T>
 377 inline std::tuple<T, T> halfAttractiveScalarForce(T k, T x0, T x)
 378 {
 379     T dx = x - x0;
 380     T dx2 = dx * dx;
 381     T dx3 = dx2 * dx;
 382     T dx4 = dx2 * dx2;
 383
 384     T epot = -0.5 * k * dx4;
 385     T force = -2.0 * k * dx3;
 386
 387     return std::make_tuple(force, epot);
 388
 389     /* That was 10 flops */
 390 }
 391
 392 /*! \brief kernel to calculate the scalar part for half attractive potential bond forces
 393  *         for non-zero lambda to interpolate between A and B states
 394  *
 395  * \param kA spring constant state A
 396  * \param kB spring constant state B
 397  * \param xA equilibrium distance state A
 398  * \param xB equilibrium distance state B
 399  * \param x  input bond length
 400  * \param lambda interpolation factor between A and B state
 401  *
 402  * \return tuple<force, potential energy, lambda-interpolated energy>
 403  */
 404 template <class T>
 405 inline std::tuple<T, T, T> halfAttractiveScalarForce(T kA, T kB, T xA, T xB, T x, T lambda)
 406 {
 407     T L1 = 1.0 - lambda;
 408     T kk = L1 * kA + lambda * kB;
 409     T x0 = L1 * xA + lambda * xB;
 410
 411     T dx  = x - x0;
 412     T dx2 = dx * dx;
 413     T dx3 = dx2 * dx;
 414     T dx4 = dx2 * dx2;
 415
 416     T epot = -0.5 * kk * dx4;
 417     T force = -2.0 * kk * dx3;
 418     T dvdlambda = 0.5 * (kB - kA) * dx4 + (2.0 * (xA - xB) * kk * dx3);
 419
 420     return std::make_tuple(force, epot, dvdlambda);
 421
 422     /* That was 29 flops */
 423 }
 424
 425 template <class T>
 426 inline auto bondKernel(T dr, const HalfAttractiveQuarticBondType& bond)
 427 {
 428     return halfAttractiveScalarForce(bond.forceConstant(), bond.equilConstant(), dr);
 429 }
 430
 431
 432 //! Three-center interaction type kernels
 433
 434 /*! \brief kernel to calculate the scalar part for linear angle forces
 435  *         for lambda = 0
 436  *
 437  * \param k  force constant
 438  * \param a  contribution of rij vector
 439  * \param dr weighted distance metric
 440  *
 441  * \return tuple<force, potential energy>
 442  */
 443 template <class T>
 444 inline std::tuple<T, T, T> linearAnglesScalar(T k, T a, T dr)
 445 {
 446     T b = T(1.0) - a;
 447
 448     T kdr  = k * dr;
 449     T epot = 0.5 * kdr * dr;
 450
 451     T ci = a * k;
 452     T ck = b * k;
 453
 454     return std::make_tuple(ci, ck, epot);
 455
 456     /* That was 5 flops */
 457 }
 458
 459 template <class T>
 460 inline auto threeCenterKernel(T dr, const LinearAngle& angle)
 461 {
 462     return linearAnglesScalar(angle.forceConstant(), angle.equilConstant(), dr);
 463 }
 464
 465 //! Harmonic Angle
 466 template <class T>
 467 inline auto threeCenterKernel(T dr, const HarmonicAngle& angle)
 468 {
 469     return harmonicScalarForce(angle.forceConstant(), angle.equilConstant(), dr);
 470 }
 471
 472 //! Cosine based (GROMOS-96) Angle
 473 template <class T>
 474 inline auto threeCenterKernel(T dr, const G96Angle& angle)
 475 {
 476     auto costheta = std::cos(dr);
 477     auto feTuple = g96ScalarForce(angle.forceConstant(), angle.equilConstant(), costheta);
 478
 479     // The above kernel call effectively computes the derivative of the potential with respect to
 480     // cos(theta). However, we need the derivative with respect to theta. We use this extra -sin(theta)
 481     // factor to account for this before the forces are spread between the particles.
 482
 483     std::get<0>(feTuple) *= -std::sqrt(1 - costheta*costheta);
 484     return feTuple;
 485 }
 486
 487 /*! \brief kernel to calculate the scalar part for cross bond-bond forces
 488  *         for lambda = 0
 489  *
 490  * \param k force constant
 491  * \param r0ij equilibrium distance between particles i & j
 492  * \param r0kj equilibrium distance between particles k & j
 493  * \param rij  input bond length between particles i & j
 494  * \param rkj  input bond length between particles k & j
 495  *
 496  * \return tuple<force scalar i, force scalar k, potential energy>
 497  */
 498
 499 template <class T>
 500 inline std::tuple<T, T, T> crossBondBondScalarForce(T k, T r0ij, T r0kj, T rij, T rkj)
 501 {
 502     T si = rij - r0ij;
 503     T sk = rkj - r0kj;
 504
 505     T epot = k * si * sk;
 506
 507     T ci = -k * sk / rij;
 508     T ck = -k * si / rkj;
 509
 510     return std::make_tuple(ci, ck, epot);
 511
 512     /* That was 8 flops */
 513 }
 514
 515 //! Cross bond-bond interaction
 516 template <class T>
 517 inline auto threeCenterKernel(T drij, T drkj, const CrossBondBond& crossBondBond)
 518 {
 519     return crossBondBondScalarForce(crossBondBond.forceConstant(), crossBondBond.equilDistanceIJ(), crossBondBond.equilDistanceKJ(), drij, drkj);
 520 }
 521
 522 /*! \brief kernel to calculate the scalar part for cross bond-angle forces
 523  *         for lambda = 0
 524  *
 525  * \param k force constant
 526  * \param r0ij equilibrium distance between particles i & j
 527  * \param r0kj equilibrium distance between particles k & j
 528  * \param r0ik equilibrium distance between particles i & k
 529  * \param rij  input bond length between particles i & j
 530  * \param rkj  input bond length between particles k & j
 531  * \param rik  input bond length between particles i & k
 532  *
 533  * \return tuple<force, potential energy>
 534  */
 535
 536 template <class T>
 537 inline std::tuple<T, T, T, T> crossBondAngleScalarForce(T k, T r0ij, T r0kj, T r0ik, T rij, T rkj, T rik)
 538 {
 539     T sij = rij - r0ij;
 540     T skj = rkj - r0kj;
 541     T sik = rik - r0ik;
 542
 543     T epot = k * sik * (sij + skj);
 544
 545     T ci = -k * sik / rij;
 546     T cj = -k * sik / rkj;
 547     T ck = -k * (sij + skj) / rik;
 548
 549     return std::make_tuple(ci, cj, ck, epot);
 550
 551     /* That was 13 flops */
 552 }
 553
 554 //! Cross bond-bond interaction
 555 template <class T>
 556 inline auto crossBondAngleKernel(T drij, T drkj, T drik, const CrossBondAngle& crossBondAngle)
 557 {
 558     return crossBondAngleScalarForce(crossBondAngle.forceConstant(), crossBondAngle.equilDistanceIJ(), crossBondAngle.equilDistanceKJ(), crossBondAngle.equilDistanceIK(), drij, drkj, drik);
 559 }
 560
 561 //! Quartic Angle
 562 template <class T>
 563 inline auto threeCenterKernel(T dr, const QuarticAngle& angle)
 564 {
 565     T dt = dr - angle.equilConstant();       /*  1          */
 566
 567     T force  = 0;
 568     T energy = angle.forceConstant(0);
 569     T dtp  = 1.0;
 570     for (auto j = 1; j <= 4; j++)
 571     { /* 24     */
 572         T c = angle.forceConstant(j);
 573         force -= j * c * dtp;               /*  3          */
 574         dtp *= dt;                          /*  1          */
 575         energy += c * dtp;                  /*  2          */
 576     }
 577
 578     /* TOTAL 25 */
 579     return std::make_tuple(force, energy);
 580 }
 581
 582 //! \brief Restricted Angle potential. Returns scalar force and energy
 583 template <class T>
 584 inline auto threeCenterKernel(T theta, const RestrictedAngle& angle)
 585 {
 586     T costheta = std::cos(theta);
 587     auto [force, ePot] = harmonicScalarForce(angle.forceConstant(), angle.equilConstant(), costheta);
 588
 589     // The above kernel call effectively computes the derivative of the potential with respect to
 590     // cos(theta). However, we need the derivative with respect to theta.
 591     // This introduces the extra (cos(theta)*cos(eqAngle) - 1)/(sin(theta)^3 factor for the force
 592     // The call also computes the potential energy without the sin(theta)^-2 factor
 593
 594     T sintheta2 = (1 - costheta*costheta);
 595     T sintheta  = std::sqrt(sintheta2);
 596     force *= (costheta*angle.equilConstant() - 1)/(sintheta2*sintheta);
 597     ePot /= sintheta2;
 598     return std::make_tuple(force, ePot);
 599 }
 600
 601 //! \brief Computes and returns the proper dihedral force
 602 template <class T>
 603 inline auto fourCenterKernel(T phi, const ProperDihedral& properDihedral)
 604 {
 605     T deltaPhi = properDihedral.multiplicity() * phi - properDihedral.equilDistance();
 606     T force = -properDihedral.forceConstant() * properDihedral.multiplicity() * std::sin(deltaPhi);
 607     T ePot = properDihedral.forceConstant() * ( 1 + std::cos(deltaPhi) );
 608     return std::make_tuple(force, ePot);
 609 }
 610
 611
 612 //! \brief Ensure that a geometric quantity lies in (-pi, pi)
 613 template<class T>
 614 inline void makeAnglePeriodic(T& angle)
 615 {
 616     if (angle >= M_PI)
 617     {
 618         angle -= 2 * M_PI;
 619     }
 620     else if (angle < -M_PI)
 621     {
 622         angle += 2 * M_PI;
 623     }
 624 }
 625
 626 /*! \brief calculate the cosine of the angle between aInput and bInput
 627  *
 628  * \tparam T       float or double
 629  * \param aInput   aInput 3D vector
 630  * \param bInput   another 3D vector
 631  * \return         the cosine of the angle between aInput and bInput
 632  *
 633  *                       ax*bx + ay*by + az*bz
 634  * cos(aInput,bInput) = -----------------------, where aInput = (ax, ay, az)
 635  *                      ||aInput|| * ||bInput||
 636  */
 637 template<class T>
 638 inline T basicVectorCosAngle(gmx::BasicVector<T> aInput, gmx::BasicVector<T> bInput)
 639 {
 640     gmx::BasicVector<double> a_double(aInput[0], aInput[1], aInput[2]);
 641     gmx::BasicVector<double> b_double(bInput[0], bInput[1], bInput[2]);
 642
 643     double numerator     = dot(a_double, b_double);
 644     double denominatorSq = dot(a_double, a_double) * dot(b_double, b_double);
 645
 646     T cosval = (denominatorSq > 0) ? static_cast<T>(numerator * gmx::invsqrt(denominatorSq)) : 1;
 647     cosval   = std::min(cosval, T(1.0));
 648
 649     /* 25 TOTAL */
 650     return std::max(cosval, T(-1.0));
 651 }
 652
 653 /*! \brief compute the angle between vectors a and b
 654  *
 655  * \tparam T    scalar type (float, double, or similar)
 656  * \param a     a 3D vector
 657  * \param b     another 3D vector
 658  * \return      the angle between a and b
 659  *
 660  * This routine calculates the angle between a & b without any loss of accuracy close to 0/PI.
 661  *
 662  * Note: This function is not (yet) implemented for the C++ replacement of the
 663  * deprecated rvec and dvec.
 664  */
 665 template<class T>
 666 inline T basicVectorAngle(gmx::BasicVector<T> a, gmx::BasicVector<T> b)
 667 {
 668     gmx::BasicVector<T> w = cross(a, b);
 669
 670     T wlen = norm(w);
 671     T s    = dot(a, b);
 672
 673     return std::atan2(wlen, s);
 674 }
 675
 676 /*! \brief Computes the dihedral phi angle and two cross products
 677  *
 678  * \tparam T        scalar type (float or double, or similar)
 679  * \param[in] dxIJ
 680  * \param[in] dxKJ
 681  * \param[in] dxKL
 682  * \param[out] m    output for \p dxIJ x \p dxKJ
 683  * \param[out] n    output for \p dxKJ x \p dxKL
 684  * \return          the angle between m and n
 685  */
 686 template<class T>
 687 inline T dihedralPhi(gmx::BasicVector<T> dxIJ,
 688                      gmx::BasicVector<T> dxKJ,
 689                      gmx::BasicVector<T> dxKL,
 690                      gmx::BasicVector<T>* m,
 691                      gmx::BasicVector<T>* n)
 692 {
 693     *m = cross(dxIJ, dxKJ);
 694     *n = cross(dxKJ, dxKL);
 695     T phi  = basicVectorAngle(*m, *n);
 696     T ipr  = dot(dxIJ, *n);
 697     T sign = (ipr < 0.0) ? -1.0 : 1.0;
 698     phi    = sign * phi;
 699     return phi;
 700 }
 701
 702 //! \brief Computes and returns the improper dihedral force
 703 template <class T>
 704 inline auto fourCenterKernel(T phi, const ImproperDihedral& improperDihedral)
 705 {
 706     T deltaPhi = phi - improperDihedral.equilConstant();
 707     /* deltaPhi cannot be outside (-pi,pi) */
 708     makeAnglePeriodic(deltaPhi);
 709     const T force = -improperDihedral.forceConstant()  * deltaPhi;
 710     const T ePot = 0.5 * improperDihedral.forceConstant() * deltaPhi * deltaPhi;
 711     return std::make_tuple(force, ePot);
 712 }
 713
 714 //! \brief Computes and returns the Ryckaert-Belleman dihedral force
 715 template <class T>
 716 inline auto fourCenterKernel(T phi, const RyckaertBellemanDihedral& ryckaertBellemanDihedral)
 717 {
 718     /* Change to polymer convention */
 719     const T localPhi = (phi < 0) ? (phi += M_PI) : (phi -= M_PI);
 720     T cos_phi = std::cos(localPhi);
 721     T ePot = ryckaertBellemanDihedral[0];
 722     T force = 0;
 723     T cosineFactor = 1;
 724
 725     for (int i = 1; i < int(ryckaertBellemanDihedral.size()); i++)
 726     {
 727         force += ryckaertBellemanDihedral[i] * cosineFactor * i;
 728         cosineFactor *= cos_phi;
 729         ePot += cosineFactor * ryckaertBellemanDihedral[i];
 730     }
 731     /* Beware of accuracy loss, cannot use 1-sqrt(cos^2) ! */
 732     force = -force * std::sin(localPhi);
 733     return std::make_tuple(force, ePot);
 734 }
 735
 736 //! Two-center category common
 737
 738 //! \brief add shift forces, if requested (static compiler decision)
 739 template<class T, class ShiftForce>
 740 inline void
 741 addShiftForce(const gmx::BasicVector<T>& interactionForce, ShiftForce* shiftForce)
 742 {
 743     *shiftForce += interactionForce;
 744 }
 745
 746 //! \brief this will be called if shift forces are not computed (and removed by the compiler)
 747 template<class T>
 748 inline void addShiftForce([[maybe_unused]] const gmx::BasicVector<T>& fij,
 749                           [[maybe_unused]] std::nullptr_t*)
 750 {
 751 }
 752
 753 /*! \brief Spreads and accumulates the bonded forces to the two atoms and adds the virial contribution when needed
 754  *
 755  * \p shiftIndex is used as the periodic shift.
 756  */
 757 template <class T, class ShiftForce>
 758 inline void spreadTwoCenterForces(const T bondForce,
 759                                   const gmx::BasicVector<T>& dx,
 760                                   gmx::BasicVector<T>*       force_i,
 761                                   gmx::BasicVector<T>*       force_j,
 762                                   ShiftForce*                shf_ij,
 763                                   ShiftForce*                shf_c)
 764 {
 765     gmx::BasicVector<T> fij = bondForce * dx;
 766     *force_i += fij;
 767     *force_j -= fij;
 768
 769     addShiftForce(fij, shf_ij);
 770     addShiftForce(T(-1.0)*fij, shf_c);
 771     /* 15 Total */
 772 }
 773
 774 //! Three-center category common
 775
 776 /*! \brief spread force to 3 centers based on scalar force and angle
 777  *
 778  * \tparam T
 779  * \param cos_theta
 780  * \param force
 781  * \param r_ij
 782  * \param r_kj
 783  * \param force_i
 784  * \param force_j
 785  * \param force_k
 786  */
 787 template <class T, class ShiftForce>
 788 inline void spreadThreeCenterForces(T                          cos_theta,
 789                                     T                          force,
 790                                     const gmx::BasicVector<T>& r_ij,
 791                                     const gmx::BasicVector<T>& r_kj,
 792                                     gmx::BasicVector<T>*       force_i,
 793                                     gmx::BasicVector<T>*       force_j,
 794                                     gmx::BasicVector<T>*       force_k,
 795                                     ShiftForce*                shf_ij,
 796                                     ShiftForce*                shf_kj,
 797                                     ShiftForce*                shf_c)
 798 {
 799     T cos_theta2 = cos_theta * cos_theta;
 800     if (cos_theta2 < 1)                              /*   1             */
 801     {
 802         T st    = force / std::sqrt(1 - cos_theta2); /*  12             */
 803         T sth   = st * cos_theta;                    /*   1             */
 804         T nrij2 = dot(r_ij, r_ij);                   /*   5             */
 805         T nrkj2 = dot(r_kj, r_kj);                   /*   5             */
 806
 807         T cik = st / std::sqrt(nrij2 * nrkj2);       /*  11             */
 808         T cii = sth / nrij2;                         /*   1             */
 809         T ckk = sth / nrkj2;                         /*   1             */
 810
 811         /*  33 */
 812         gmx::BasicVector<T> f_i = cii * r_ij - cik * r_kj;
 813         gmx::BasicVector<T> f_k = ckk * r_kj - cik * r_ij;
 814         gmx::BasicVector<T> f_j = T(-1.0) * (f_i + f_k);
 815         *force_i += f_i;
 816         *force_j += f_j;
 817         *force_k += f_k;
 818
 819         addShiftForce(f_i, shf_ij);
 820         addShiftForce(f_j, shf_c);
 821         addShiftForce(f_k, shf_kj);
 822
 823     } /* 70 TOTAL       */
 824 }
 825
 826 //! Four-center category common
 827 template <class T, class ShiftForce>
 828 inline void spreadFourCenterForces(T force,
 829                                    const gmx::BasicVector<T>& dxIJ,
 830                                    const gmx::BasicVector<T>& dxJK,
 831                                    const gmx::BasicVector<T>& dxKL,
 832                                    const gmx::BasicVector<T>& m,
 833                                    const gmx::BasicVector<T>& n,
 834                                    gmx::BasicVector<T>* force_i,
 835                                    gmx::BasicVector<T>* force_j,
 836                                    gmx::BasicVector<T>* force_k,
 837                                    gmx::BasicVector<T>* force_l,
 838                                    ShiftForce* shf_ij,
 839                                    ShiftForce* shf_kj,
 840                                    ShiftForce* shf_lj,
 841                                    ShiftForce* shf_c)
 842 {
 843     T norm2_m = dot(m, m);                             /* 5 */
 844     T norm2_n = dot(n, n);                             /* 5 */
 845     T norm2_jk = dot(dxJK, dxJK);                      /* 5 */
 846     T toler = norm2_jk * GMX_REAL_EPS;
 847     if ((norm2_m > toler) && (norm2_n > toler))
 848     {
 849         T rcp_norm2_jk = 1.0f / norm2_jk;              /* 1 */
 850         T norm_jk = std::sqrt(norm2_jk);               /* 10 */
 851
 852         T a = -force * norm_jk / norm2_m;              /* 11 */
 853         gmx::BasicVector<T> f_i = a * m;               /* 3 */
 854
 855         T b = force * norm_jk / norm2_n;               /* 11 */
 856         gmx::BasicVector<T> f_l = b * n;               /* 3 */
 857
 858         T p = rcp_norm2_jk * dot(dxIJ, dxJK);          /* 6 */
 859         T q = rcp_norm2_jk * dot(dxKL, dxJK);          /* 6 */
 860         gmx::BasicVector<T> svec = p * f_i - q * f_l;  /* 9 */
 861
 862         gmx::BasicVector<T> f_j = svec - f_i;          /* 3 */
 863         gmx::BasicVector<T> f_k = T(-1.0)*svec - f_l;  /* 6 */
 864
 865         *force_i += f_i;                               /* 3 */
 866         *force_j += f_j;                               /* 3 */
 867         *force_k += f_k;                               /* 3 */
 868         *force_l += f_l;                               /* 3 */
 869
 870         addShiftForce(f_i, shf_ij);
 871         addShiftForce(f_j, shf_c);
 872         addShiftForce(f_k, shf_kj);
 873         addShiftForce(f_l, shf_lj);
 874     }
 875 }
 876
 877 } // namespace nblib
 878
 879 #undef NBLIB_ALWAYS_INLINE
 880
 881 #endif // NBLIB_LISTEDFORCES_KERNELS_HPP