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