api/nblib/listed_forces/kernels.hpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2020,2021, 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/vec.h"
  51 #include "nblib/basicdefinitions.h"
  52 #include "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 * x;
 137     T epot = 0.25 * k * dx2;
 138
 139     return std::make_tuple(force, epot);
 140
 141     /* That was 7 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 * x;
 167     T epot = 0.25 * 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. Forth power case
 177 template <class T>
 178 inline auto bondKernel(T dr, const G96BondType& bond)
 179 {
 180     // NOTE: Not assuming GROMACS' convention of storing squared bond.equilConstant() for this type
 181     return g96ScalarForce(bond.forceConstant(), bond.equilConstant() * bond.equilConstant(), dr * dr);
 182 }
 183
 184
 185 /*! \brief kernel to calculate the scalar part for the morse pontential bond forces
 186  *         for lambda = 0
 187  *
 188  * \param k force constant
 189  * \param beta beta exponent
 190  * \param x0 equilibrium distance
 191  * \param x  input bond length
 192  *
 193  * \return tuple<force, potential energy>
 194  */
 195 template <class T>
 196 inline std::tuple<T, T> morseScalarForce(T k, T beta, T x0, T x)
 197 {
 198     T exponent = std::exp(-beta * (x - x0));      /* 12 */
 199     T omexp = 1.0 - exponent;                     /*  1 */
 200     T kexp = k * omexp;                           /*  1 */
 201
 202     T epot = kexp * omexp;                        /*  1 */
 203     T force = -2.0 * beta * exponent * omexp;     /*  4 */
 204
 205     return std::make_tuple(force, epot);
 206
 207     /* That was 23 flops */
 208 }
 209
 210 /*! \brief kernel to calculate the scalar part for morse potential bond forces
 211  *         for non-zero lambda to interpolate between A and B states
 212  *
 213  * \param kA force constant state A
 214  * \param kB force constant state B
 215  * \param betaA beta exponent state A
 216  * \param betaB beta exponent state B
 217  * \param xA equilibrium distance state A
 218  * \param xB equilibrium distance state B
 219  * \param x  input bond length
 220  * \param lambda interpolation factor between A and B state
 221  *
 222  * \return tuple<force, potential energy, lambda-interpolated energy>
 223  */
 224 template <class T>
 225 inline std::tuple<T, T, T> morseScalarForce(T kA, T kB, T betaA, T betaB, T xA, T xB, T x, T lambda)
 226 {
 227     T L1 = 1.0 - lambda;                          /* 1 */
 228     T x0 = L1 * xA + lambda * xB;                 /* 3 */
 229     T beta = L1 * betaA + lambda * betaB;         /* 3 */
 230     T k = L1 * kA + lambda * kB;                  /* 3 */
 231
 232     T exponent = std::exp(-beta * (x - x0));      /* 12 */
 233     T omexp = 1.0 - exponent;                     /*  1 */
 234     T kexp = k * omexp;                           /*  1 */
 235
 236     T epot = kexp * omexp;                        /*  1 */
 237     T force = -2.0 * beta * exponent * omexp;     /*  4 */
 238
 239     T dvdlambda = (kB - kA) * omexp * omexp
 240                     - (2.0 - 2.0 * omexp) * omexp * k
 241                     * ((xB - xA) * beta - (betaB - betaA) * (x - x0)); /* 15 */
 242
 243     return std::make_tuple(force, epot, dvdlambda);
 244
 245     /* That was 44 flops */
 246 }
 247
 248 //! Abstraction layer for different 2-center bonds. Morse case
 249 template <class T>
 250 inline auto bondKernel(T dr, const MorseBondType& bond)
 251 {
 252     return morseScalarForce(bond.forceConstant(), bond.exponent(), bond.equilDistance(), dr);
 253 }
 254
 255
 256 /*! \brief kernel to calculate the scalar part for the FENE pontential bond forces
 257  *         for lambda = 0
 258  *
 259  * \param k spring constant
 260  * \param x0 equilibrium distance
 261  * \param x  input bond length
 262  *
 263  * \return tuple<force, potential energy>
 264  */
 265 template <class T>
 266 inline std::tuple<T, T> FENEScalarForce(T k, T x0, T x)
 267 {
 268     T x02 = x0 * x0;
 269     T x2 = x * x;
 270
 271     T omx2_ox02 = 1.0 - (x2 / x02);
 272
 273     T epot = -0.5 * k * x02 * std::log(omx2_ox02);
 274     T force = -k * x / omx2_ox02;
 275
 276     return std::make_tuple(force, epot);
 277
 278     /* That was 24 flops */
 279 }
 280
 281 // TODO: Implement the free energy version of FENE (finitely extensible nonlinear elastic) bond types
 282
 283 //! Abstraction layer for different 2-center bonds. FENE case
 284 template <class T>
 285 inline auto bondKernel(T dr, const FENEBondType& bond)
 286 {
 287     return FENEScalarForce(bond.forceConstant(), bond.equilConstant(), dr);
 288 }
 289
 290
 291 /*! \brief kernel to calculate the scalar part for cubic potential bond forces
 292  *         for lambda = 0
 293  *
 294  * \param kc cubic spring constant
 295  * \param kq quadratic 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> cubicScalarForce(T kc, T kq, T x0, T x)
 303 {
 304     T dx = x - x0;
 305
 306     T kdist  = kq * dx;
 307     T kdist2 = kdist * dx;
 308
 309     T epot = kdist2 + (kc * kdist2 * dx);
 310     T force = -((2.0 * kdist) + (3.0 * kdist2 * kc));
 311
 312     return std::make_tuple(force, epot);
 313
 314     /* That was 16 flops */
 315 }
 316
 317 // TODO: Implement the free energy version of Cubic bond types
 318
 319 template <class T>
 320 inline auto bondKernel(T dr, const CubicBondType& bond)
 321 {
 322     return cubicScalarForce(bond.cubicForceConstant(), bond.quadraticForceConstant(), bond.equilDistance(), dr);
 323 }
 324
 325
 326 /*! \brief kernel to calculate the scalar part for half attractive potential bond forces
 327  *         for lambda = 0
 328  *
 329  * \param k spring constant
 330  * \param x0 equilibrium distance
 331  * \param x  input bond length
 332  *
 333  * \return tuple<force, potential energy>
 334  */
 335 template <class T>
 336 inline std::tuple<T, T> halfAttractiveScalarForce(T k, T x0, T x)
 337 {
 338     T dx = x - x0;
 339     T dx2 = dx * dx;
 340     T dx3 = dx2 * dx;
 341     T dx4 = dx2 * dx2;
 342
 343     T epot = -0.5 * k * dx4;
 344     T force = -2.0 * k * dx3;
 345
 346     return std::make_tuple(force, epot);
 347
 348     /* That was 10 flops */
 349 }
 350
 351 /*! \brief kernel to calculate the scalar part for half attractive potential bond forces
 352  *         for non-zero lambda to interpolate between A and B states
 353  *
 354  * \param kA spring constant state A
 355  * \param kB spring constant state B
 356  * \param xA equilibrium distance state A
 357  * \param xB equilibrium distance state B
 358  * \param x  input bond length
 359  * \param lambda interpolation factor between A and B state
 360  *
 361  * \return tuple<force, potential energy, lambda-interpolated energy>
 362  */
 363 template <class T>
 364 inline std::tuple<T, T, T> halfAttractiveScalarForce(T kA, T kB, T xA, T xB, T x, T lambda)
 365 {
 366     T L1 = 1.0 - lambda;
 367     T kk = L1 * kA + lambda * kB;
 368     T x0 = L1 * xA + lambda * xB;
 369
 370     T dx  = x - x0;
 371     T dx2 = dx * dx;
 372     T dx3 = dx2 * dx;
 373     T dx4 = dx2 * dx2;
 374
 375     T epot = -0.5 * kk * dx4;
 376     T force = -2.0 * kk * dx3;
 377     T dvdlambda = 0.5 * (kB - kA) * dx4 + (2.0 * (xA - xB) * kk * dx3);
 378
 379     return std::make_tuple(force, epot, dvdlambda);
 380
 381     /* That was 29 flops */
 382 }
 383
 384 template <class T>
 385 inline auto bondKernel(T dr, const HalfAttractiveQuarticBondType& bond)
 386 {
 387     return halfAttractiveScalarForce(bond.forceConstant(), bond.equilConstant(), dr);
 388 }
 389
 390
 391 //! Three-center interaction type kernels
 392
 393 // linear angles go here
 394 // quartic angles go here
 395
 396 //! Three-center interaction type dispatch
 397
 398 template <class T>
 399 inline auto threeCenterKernel(T dr, const HarmonicAngle& angle)
 400 {
 401     return harmonicScalarForce(angle.forceConstant(), angle.equilConstant(), dr);
 402 }
 403
 404
 405 //! \brief Computes and returns the proper dihedral force
 406 template <class T>
 407 inline auto fourCenterKernel(T phi, const ProperDihedral& properDihedral)
 408 {
 409     const T deltaPhi = properDihedral.multiplicity() * phi - properDihedral.equilDistance();
 410     const T force = -properDihedral.forceConstant() * properDihedral.multiplicity() * std::sin(deltaPhi);
 411     const T ePot = properDihedral.forceConstant() * ( 1 + std::cos(deltaPhi) );
 412     return std::make_tuple(force, ePot);
 413 }
 414
 415
 416 //! \brief Ensure that a geometric quantity lies in (-pi, pi)
 417 static inline void makeAnglePeriodic(real& angle)
 418 {
 419     if (angle >= M_PI)
 420     {
 421         angle -= 2 * M_PI;
 422     }
 423     else if (angle < -M_PI)
 424     {
 425         angle += 2 * M_PI;
 426     }
 427 }
 428
 429 //! \brief Computes and returns a dihedral phi angle
 430 static inline real dihedralPhi(rvec dxIJ, rvec dxKJ, rvec dxKL, rvec m, rvec n)
 431 {
 432     cprod(dxIJ, dxKJ, m);
 433     cprod(dxKJ, dxKL, n);
 434     real phi  = gmx_angle(m, n);
 435     real ipr  = iprod(dxIJ, n);
 436     real sign = (ipr < 0.0) ? -1.0 : 1.0;
 437     phi       = sign * phi;
 438     return phi;
 439 }
 440
 441 //! \brief Computes and returns the improper dihedral force
 442 template <class T>
 443 inline auto fourCenterKernel(T phi, const ImproperDihedral& improperDihedral)
 444 {
 445     T deltaPhi = phi - improperDihedral.equilConstant();
 446     /* deltaPhi cannot be outside (-pi,pi) */
 447     makeAnglePeriodic(deltaPhi);
 448     const T force = -improperDihedral.forceConstant()  * deltaPhi;
 449     const T ePot = 0.5 * improperDihedral.forceConstant() * deltaPhi * deltaPhi;
 450     return std::make_tuple(force, ePot);
 451 }
 452
 453 //! \brief Computes and returns the Ryckaert-Belleman dihedral force
 454 template <class T>
 455 inline auto fourCenterKernel(T phi, const RyckaertBellemanDihedral& ryckaertBellemanDihedral)
 456 {
 457     /* Change to polymer convention */
 458     const T localPhi = (phi < 0) ? (phi += M_PI) : (phi -= M_PI);
 459     T cos_phi = std::cos(localPhi);
 460     T ePot = ryckaertBellemanDihedral[0];
 461     T force = 0;
 462     T cosineFactor = 1;
 463
 464     for (int i = 1; i < int(ryckaertBellemanDihedral.size()); i++)
 465     {
 466         force += ryckaertBellemanDihedral[i] * cosineFactor * i;
 467         cosineFactor *= cos_phi;
 468         ePot += cosineFactor * ryckaertBellemanDihedral[i];
 469     }
 470     /* Beware of accuracy loss, cannot use 1-sqrt(cos^2) ! */
 471     force = -force * std::sin(localPhi);
 472     return std::make_tuple(force, ePot);
 473 }
 474
 475 //! Two-center category common
 476
 477 /*! \brief Spreads and accumulates the bonded forces to the two atoms and adds the virial contribution when needed
 478  *
 479  * \p shiftIndex is used as the periodic shift.
 480  */
 481 template <class T>
 482 inline void spreadTwoCenterForces(const T bondForce,
 483                                   const gmx::RVec& dx,
 484                                   gmx::RVec* force_i,
 485                                   gmx::RVec* force_j)
 486 {
 487     for (int m = 0; m < dimSize; m++) /*  15          */
 488     {
 489         const T fij = bondForce * dx[m];
 490         (*force_i)[m] += fij;
 491         (*force_j)[m] -= fij;
 492     }
 493 }
 494
 495 //! Three-center category common
 496
 497 /*! \brief spread force to 3 centers based on scalar force and angle
 498  *
 499  * @tparam T
 500  * @param cos_theta
 501  * @param force
 502  * @param r_ij
 503  * @param r_kj
 504  * @param force_i
 505  * @param force_j
 506  * @param force_k
 507  */
 508 template <class T>
 509 inline void spreadThreeCenterForces(T cos_theta,
 510                              T force,
 511                              const gmx::RVec& r_ij,
 512                              const gmx::RVec& r_kj,
 513                              gmx::RVec* force_i,
 514                              gmx::RVec* force_j,
 515                              gmx::RVec* force_k)
 516 {
 517     T cos_theta2 = cos_theta * cos_theta;
 518     if (cos_theta2 < 1)
 519     {
 520         T st    = force / std::sqrt(1 - cos_theta2); /*  12             */
 521         T sth   = st * cos_theta;                      /*   1           */
 522         T nrij2 = dot(r_ij, r_ij);                   /*   5             */
 523         T nrkj2 = dot(r_kj, r_kj);                   /*   5             */
 524
 525         T nrij_1 = 1.0 / std::sqrt(nrij2); /*  10               */
 526         T nrkj_1 = 1.0 / std::sqrt(nrkj2); /*  10               */
 527
 528         T cik = st * nrij_1 * nrkj_1;  /*   2           */
 529         T cii = sth * nrij_1 * nrij_1; /*   2           */
 530         T ckk = sth * nrkj_1 * nrkj_1; /*   2           */
 531
 532         gmx::RVec f_i{0, 0, 0};
 533         gmx::RVec f_j{0, 0, 0};
 534         gmx::RVec f_k{0, 0, 0};
 535         for (int m = 0; m < dimSize; m++)
 536         { /*  39                */
 537             f_i[m] = -(cik * r_kj[m] - cii * r_ij[m]);
 538             f_k[m] = -(cik * r_ij[m] - ckk * r_kj[m]);
 539             f_j[m] = -f_i[m] - f_k[m];
 540             (*force_i)[m] += f_i[m];
 541             (*force_j)[m] += f_j[m];
 542             (*force_k)[m] += f_k[m];
 543         }
 544     } /* 161 TOTAL      */
 545 }
 546
 547 //! Four-center category common
 548 template <class T>
 549 inline void spreadFourCenterForces(T force, rvec dxIJ, rvec dxJK, rvec dxKL, rvec m, rvec n,
 550                             gmx::RVec* force_i,
 551                             gmx::RVec* force_j,
 552                             gmx::RVec* force_k,
 553                             gmx::RVec* force_l)
 554 {
 555     rvec f_i, f_j, f_k, f_l;
 556     rvec uvec, vvec, svec;
 557     T iprm  = iprod(m, m);       /*  5    */
 558     T iprn  = iprod(n, n);       /*  5  */
 559     T nrkj2 = iprod(dxJK, dxJK); /*  5  */
 560     T toler = nrkj2 * GMX_REAL_EPS;
 561     if ((iprm > toler) && (iprn > toler))
 562     {
 563         T nrkj_1 = gmx::invsqrt(nrkj2);  /* 10  */
 564         T nrkj_2 = nrkj_1 * nrkj_1;      /*  1  */
 565         T nrkj   = nrkj2 * nrkj_1;       /*  1  */
 566         T a      = -force * nrkj / iprm; /* 11  */
 567         svmul(a, m, f_i);              /*  3    */
 568         T b = force * nrkj / iprn;       /* 11  */
 569         svmul(b, n, f_l);              /*  3  */
 570         T p = iprod(dxIJ, dxJK);         /*  5  */
 571         p *= nrkj_2;                   /*  1    */
 572         T q = iprod(dxKL, dxJK);         /*  5  */
 573         q *= nrkj_2;                   /*  1    */
 574         svmul(p, f_i, uvec);           /*  3    */
 575         svmul(q, f_l, vvec);           /*  3    */
 576         rvec_sub(uvec, vvec, svec);    /*  3    */
 577         rvec_sub(f_i, svec, f_j);      /*  3    */
 578         rvec_add(f_l, svec, f_k);      /*  3    */
 579         rvec_inc(force_i->as_vec(), f_i);           /*  3       */
 580         rvec_dec(force_j->as_vec(), f_j);           /*  3       */
 581         rvec_dec(force_k->as_vec(), f_k);           /*  3       */
 582         rvec_inc(force_l->as_vec(), f_l);           /*  3       */
 583     }
 584 }
 585
 586 } // namespace nblib
 587 #endif // NBLIB_LISTEDFORCES_KERNELS_HPP