2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
35 /*! \inpublicapi \file
37 * Implements kernels for nblib supported bondtypes
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>
45 #ifndef NBLIB_LISTEDFORCES_KERNELS_HPP
46 #define NBLIB_LISTEDFORCES_KERNELS_HPP
50 #include "gromacs/math/functions.h"
51 #include "gromacs/math/vectypes.h"
52 #include "nblib/listed_forces/bondtypes.h"
54 #define NBLIB_ALWAYS_INLINE __attribute((always_inline))
59 /*! \brief kernel to calculate the scalar part for simple harmonic bond forces
62 * \param k spring constant
63 * \param x0 equilibrium distance
64 * \param x input bond length
66 * \return tuple<force, potential energy>
69 inline std::tuple<T, T> harmonicScalarForce(T k, T x0, T x)
75 T epot = 0.5 * k * dx2;
77 return std::make_tuple(force, epot);
79 /* That was 6 flops */
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
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
92 * \return tuple<force, potential energy, lambda-interpolated energy>
95 inline std::tuple<T, T, T> harmonicScalarForce(T kA, T kB, T xA, T xB, T x, T lambda)
97 // code unchanged relative to Gromacs
100 T kk = L1 * kA + lambda * kB;
101 T x0 = L1 * xA + lambda * xB;
107 T epot = 0.5 * kk * dx2;
108 T dvdlambda = 0.5 * (kB - kA) * dx2 + (xA - xB) * kk * dx;
110 return std::make_tuple(force, epot, dvdlambda);
112 /* That was 19 flops */
115 //! abstraction layer for different 2-center bonds
117 inline auto bondKernel(T dr, const HarmonicBondType& bond)
119 return harmonicScalarForce(bond.forceConstant(), bond.equilConstant(), dr);
123 /*! \brief kernel to calculate the scalar part for the forth power pontential bond forces
126 * \param k spring constant
127 * \param x0 squared equilibrium distance
128 * \param x squared input bond length
130 * \return tuple<force, potential energy>
133 inline std::tuple<T, T> g96ScalarForce(T k, T x0, T x)
139 T epot = 0.5 * k * dx2;
141 return std::make_tuple(force, epot);
143 /* That was 6 flops */
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
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
156 * \return tuple<force, potential energy, lambda-interpolated energy>
159 inline std::tuple<T, T, T> g96ScalarForce(T kA, T kB, T xA, T xB, T x, T lambda)
162 T kk = L1 * kA + lambda * kB;
163 T x0 = L1 * xA + lambda * xB;
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;
173 return std::make_tuple(force, epot, dvdlambda);
175 /* That was 21 flops */
178 //! Abstraction layer for different 2-center bonds. Fourth power case
180 inline auto bondKernel(T dr, const G96BondType& bond)
182 auto [force, ePot] = g96ScalarForce(bond.forceConstant(), bond.equilConstant(), dr*dr);
185 return std::make_tuple(force, ePot);
189 /*! \brief kernel to calculate the scalar part for the morse pontential bond forces
192 * \param k force constant
193 * \param beta beta exponent
194 * \param x0 equilibrium distance
195 * \param x input bond length
197 * \return tuple<force, potential energy>
200 inline std::tuple<T, T> morseScalarForce(T k, T beta, T x0, T x)
202 T exponent = std::exp(-beta * (x - x0)); /* 12 */
203 T omexp = 1.0 - exponent; /* 1 */
204 T kexp = k * omexp; /* 1 */
206 T epot = kexp * omexp; /* 1 */
207 T force = -2.0 * beta * exponent * kexp; /* 4 */
209 return std::make_tuple(force, epot);
211 /* That was 20 flops */
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
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
226 * \return tuple<force, potential energy, lambda-interpolated energy>
229 inline std::tuple<T, T, T> morseScalarForce(T kA, T kB, T betaA, T betaB, T xA, T xB, T x, T lambda)
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 */
236 T exponent = std::exp(-beta * (x - x0)); /* 12 */
237 T omexp = 1.0 - exponent; /* 1 */
238 T kexp = k * omexp; /* 1 */
240 T epot = kexp * omexp; /* 1 */
241 T force = -2.0 * beta * exponent * kexp; /* 4 */
243 T dvdlambda = (kB - kA) * omexp * omexp
244 - (2.0 - 2.0 * omexp) * omexp * k
245 * ((xB - xA) * beta - (betaB - betaA) * (x - x0)); /* 15 */
247 return std::make_tuple(force, epot, dvdlambda);
249 /* That was 44 flops */
252 //! Abstraction layer for different 2-center bonds. Morse case
254 inline auto bondKernel(T dr, const MorseBondType& bond)
256 return morseScalarForce(bond.forceConstant(), bond.exponent(), bond.equilDistance(), dr);
260 /*! \brief kernel to calculate the scalar part for the FENE pontential bond forces
263 * \param k spring constant
264 * \param x0 equilibrium distance
265 * \param x input bond length
267 * \return tuple<force, potential energy>
270 inline std::tuple<T, T> FENEScalarForce(T k, T x0, T x)
275 T omx2_ox02 = 1.0 - (x2 / x02);
277 T epot = -0.5 * k * x02 * std::log(omx2_ox02);
278 T force = -k / omx2_ox02;
280 return std::make_tuple(force, epot);
282 /* That was 24 flops */
285 // TODO: Implement the free energy version of FENE (finitely extensible nonlinear elastic) bond types
287 //! Abstraction layer for different 2-center bonds. FENE case
289 inline auto bondKernel(T dr, const FENEBondType& bond)
291 auto [force, ePot] = FENEScalarForce(bond.forceConstant(), bond.equilConstant(), dr);
293 return std::make_tuple(force, ePot);
297 /*! \brief kernel to calculate the scalar part for cubic potential bond forces
300 * \param kc cubic spring constant
301 * \param kq quadratic spring constant
302 * \param x0 equilibrium distance
303 * \param x input bond length
305 * \return tuple<force, potential energy>
308 inline std::tuple<T, T> cubicScalarForce(T kc, T kq, T x0, T x)
314 T kdist2 = kdist * dx;
316 T epot = kdist2 + (kc * kdist2 * dx);
317 T force = -((2.0 * kdist) + (3.0 * kdist2 * kc));
319 return std::make_tuple(force, epot);
321 /* That was 16 flops */
324 // TODO: Implement the free energy version of Cubic bond types
327 inline auto bondKernel(T dr, const CubicBondType& bond)
329 return cubicScalarForce(bond.cubicForceConstant(), bond.quadraticForceConstant(), bond.equilDistance(), dr);
333 /*! \brief kernel to calculate the scalar part for half attractive potential bond forces
336 * \param k spring constant
337 * \param x0 equilibrium distance
338 * \param x input bond length
340 * \return tuple<force, potential energy>
343 inline std::tuple<T, T> halfAttractiveScalarForce(T k, T x0, T x)
350 T epot = -0.5 * k * dx4;
351 T force = -2.0 * k * dx3;
353 return std::make_tuple(force, epot);
355 /* That was 10 flops */
358 /*! \brief kernel to calculate the scalar part for half attractive potential bond forces
359 * for non-zero lambda to interpolate between A and B states
361 * \param kA spring constant state A
362 * \param kB spring constant state B
363 * \param xA equilibrium distance state A
364 * \param xB equilibrium distance state B
365 * \param x input bond length
366 * \param lambda interpolation factor between A and B state
368 * \return tuple<force, potential energy, lambda-interpolated energy>
371 inline std::tuple<T, T, T> halfAttractiveScalarForce(T kA, T kB, T xA, T xB, T x, T lambda)
374 T kk = L1 * kA + lambda * kB;
375 T x0 = L1 * xA + lambda * xB;
382 T epot = -0.5 * kk * dx4;
383 T force = -2.0 * kk * dx3;
384 T dvdlambda = 0.5 * (kB - kA) * dx4 + (2.0 * (xA - xB) * kk * dx3);
386 return std::make_tuple(force, epot, dvdlambda);
388 /* That was 29 flops */
392 inline auto bondKernel(T dr, const HalfAttractiveQuarticBondType& bond)
394 return halfAttractiveScalarForce(bond.forceConstant(), bond.equilConstant(), dr);
398 //! Three-center interaction type kernels
403 inline auto threeCenterKernel(T dr, const HarmonicAngle& angle)
405 return harmonicScalarForce(angle.forceConstant(), angle.equilConstant(), dr);
409 //! \brief Computes and returns the proper dihedral force
411 inline auto fourCenterKernel(T phi, const ProperDihedral& properDihedral)
413 T deltaPhi = properDihedral.multiplicity() * phi - properDihedral.equilDistance();
414 T force = -properDihedral.forceConstant() * properDihedral.multiplicity() * std::sin(deltaPhi);
415 T ePot = properDihedral.forceConstant() * ( 1 + std::cos(deltaPhi) );
416 return std::make_tuple(force, ePot);
420 //! \brief Ensure that a geometric quantity lies in (-pi, pi)
422 inline void makeAnglePeriodic(T& angle)
428 else if (angle < -M_PI)
434 /*! \brief calculate the cosine of the angle between aInput and bInput
436 * \tparam T float or double
437 * \param aInput aInput 3D vector
438 * \param bInput another 3D vector
439 * \return the cosine of the angle between aInput and bInput
441 * ax*bx + ay*by + az*bz
442 * cos(aInput,bInput) = -----------------------, where aInput = (ax, ay, az)
443 * ||aInput|| * ||bInput||
446 inline T basicVectorCosAngle(gmx::BasicVector<T> aInput, gmx::BasicVector<T> bInput)
448 gmx::BasicVector<double> a_double(aInput[0], aInput[1], aInput[2]);
449 gmx::BasicVector<double> b_double(bInput[0], bInput[1], bInput[2]);
451 double numerator = dot(a_double, b_double);
452 double denominatorSq = dot(a_double, a_double) * dot(b_double, b_double);
454 T cosval = (denominatorSq > 0) ? static_cast<T>(numerator * gmx::invsqrt(denominatorSq)) : 1;
455 cosval = std::min(cosval, T(1.0));
458 return std::max(cosval, T(-1.0));
461 /*! \brief compute the angle between vectors a and b
463 * \tparam T scalar type (float, double, or similar)
464 * \param a a 3D vector
465 * \param b another 3D vector
466 * \return the angle between a and b
468 * This routine calculates the angle between a & b without any loss of accuracy close to 0/PI.
470 * Note: This function is not (yet) implemented for the C++ replacement of the
471 * deprecated rvec and dvec.
474 inline T basicVectorAngle(gmx::BasicVector<T> a, gmx::BasicVector<T> b)
476 gmx::BasicVector<T> w = cross(a, b);
481 return std::atan2(wlen, s);
484 /*! \brief Computes the dihedral phi angle and two cross products
486 * \tparam T scalar type (float or double, or similar)
490 * \param[out] m output for \p dxIJ x \p dxKJ
491 * \param[out] n output for \p dxKJ x \p dxKL
492 * \return the angle between m and n
495 inline T dihedralPhi(gmx::BasicVector<T> dxIJ,
496 gmx::BasicVector<T> dxKJ,
497 gmx::BasicVector<T> dxKL,
498 gmx::BasicVector<T>* m,
499 gmx::BasicVector<T>* n)
501 *m = cross(dxIJ, dxKJ);
502 *n = cross(dxKJ, dxKL);
503 T phi = basicVectorAngle(*m, *n);
504 T ipr = dot(dxIJ, *n);
505 T sign = (ipr < 0.0) ? -1.0 : 1.0;
510 //! \brief Computes and returns the improper dihedral force
512 inline auto fourCenterKernel(T phi, const ImproperDihedral& improperDihedral)
514 T deltaPhi = phi - improperDihedral.equilConstant();
515 /* deltaPhi cannot be outside (-pi,pi) */
516 makeAnglePeriodic(deltaPhi);
517 const T force = -improperDihedral.forceConstant() * deltaPhi;
518 const T ePot = 0.5 * improperDihedral.forceConstant() * deltaPhi * deltaPhi;
519 return std::make_tuple(force, ePot);
522 //! \brief Computes and returns the Ryckaert-Belleman dihedral force
524 inline auto fourCenterKernel(T phi, const RyckaertBellemanDihedral& ryckaertBellemanDihedral)
526 /* Change to polymer convention */
527 const T localPhi = (phi < 0) ? (phi += M_PI) : (phi -= M_PI);
528 T cos_phi = std::cos(localPhi);
529 T ePot = ryckaertBellemanDihedral[0];
533 for (int i = 1; i < int(ryckaertBellemanDihedral.size()); i++)
535 force += ryckaertBellemanDihedral[i] * cosineFactor * i;
536 cosineFactor *= cos_phi;
537 ePot += cosineFactor * ryckaertBellemanDihedral[i];
539 /* Beware of accuracy loss, cannot use 1-sqrt(cos^2) ! */
540 force = -force * std::sin(localPhi);
541 return std::make_tuple(force, ePot);
544 //! Two-center category common
546 //! \brief add shift forces, if requested (static compiler decision)
547 template<class T, class ShiftForce>
549 addShiftForce(const gmx::BasicVector<T>& interactionForce, ShiftForce* shiftForce)
551 *shiftForce += interactionForce;
554 //! \brief this will be called if shift forces are not computed (and removed by the compiler)
556 inline void addShiftForce([[maybe_unused]] const gmx::BasicVector<T>& fij,
557 [[maybe_unused]] std::nullptr_t*)
561 /*! \brief Spreads and accumulates the bonded forces to the two atoms and adds the virial contribution when needed
563 * \p shiftIndex is used as the periodic shift.
565 template <class T, class ShiftForce>
566 inline void spreadTwoCenterForces(const T bondForce,
567 const gmx::BasicVector<T>& dx,
568 gmx::BasicVector<T>* force_i,
569 gmx::BasicVector<T>* force_j,
573 gmx::BasicVector<T> fij = bondForce * dx;
577 addShiftForce(fij, shf_ij);
578 addShiftForce(T(-1.0)*fij, shf_c);
582 //! Three-center category common
584 /*! \brief spread force to 3 centers based on scalar force and angle
595 template <class T, class ShiftForce>
596 inline void spreadThreeCenterForces(T cos_theta,
598 const gmx::BasicVector<T>& r_ij,
599 const gmx::BasicVector<T>& r_kj,
600 gmx::BasicVector<T>* force_i,
601 gmx::BasicVector<T>* force_j,
602 gmx::BasicVector<T>* force_k,
607 T cos_theta2 = cos_theta * cos_theta;
608 if (cos_theta2 < 1) /* 1 */
610 T st = force / std::sqrt(1 - cos_theta2); /* 12 */
611 T sth = st * cos_theta; /* 1 */
612 T nrij2 = dot(r_ij, r_ij); /* 5 */
613 T nrkj2 = dot(r_kj, r_kj); /* 5 */
615 T cik = st / std::sqrt(nrij2 * nrkj2); /* 11 */
616 T cii = sth / nrij2; /* 1 */
617 T ckk = sth / nrkj2; /* 1 */
620 gmx::BasicVector<T> f_i = cii * r_ij - cik * r_kj;
621 gmx::BasicVector<T> f_k = ckk * r_kj - cik * r_ij;
622 gmx::BasicVector<T> f_j = T(-1.0) * (f_i + f_k);
627 addShiftForce(f_i, shf_ij);
628 addShiftForce(f_j, shf_c);
629 addShiftForce(f_k, shf_kj);
634 //! Four-center category common
635 template <class T, class ShiftForce>
636 inline void spreadFourCenterForces(T force,
637 const gmx::BasicVector<T>& dxIJ,
638 const gmx::BasicVector<T>& dxJK,
639 const gmx::BasicVector<T>& dxKL,
640 const gmx::BasicVector<T>& m,
641 const gmx::BasicVector<T>& n,
642 gmx::BasicVector<T>* force_i,
643 gmx::BasicVector<T>* force_j,
644 gmx::BasicVector<T>* force_k,
645 gmx::BasicVector<T>* force_l,
651 T norm2_m = dot(m, m); /* 5 */
652 T norm2_n = dot(n, n); /* 5 */
653 T norm2_jk = dot(dxJK, dxJK); /* 5 */
654 T toler = norm2_jk * GMX_REAL_EPS;
655 if ((norm2_m > toler) && (norm2_n > toler))
657 T rcp_norm2_jk = 1.0f / norm2_jk; /* 1 */
658 T norm_jk = std::sqrt(norm2_jk); /* 10 */
660 T a = -force * norm_jk / norm2_m; /* 11 */
661 gmx::BasicVector<T> f_i = a * m; /* 3 */
663 T b = force * norm_jk / norm2_n; /* 11 */
664 gmx::BasicVector<T> f_l = b * n; /* 3 */
666 T p = rcp_norm2_jk * dot(dxIJ, dxJK); /* 6 */
667 T q = rcp_norm2_jk * dot(dxKL, dxJK); /* 6 */
668 gmx::BasicVector<T> svec = p * f_i - q * f_l; /* 9 */
670 gmx::BasicVector<T> f_j = svec - f_i; /* 3 */
671 gmx::BasicVector<T> f_k = T(-1.0)*svec - f_l; /* 6 */
673 *force_i += f_i; /* 3 */
674 *force_j += f_j; /* 3 */
675 *force_k += f_k; /* 3 */
676 *force_l += f_l; /* 3 */
678 addShiftForce(f_i, shf_ij);
679 addShiftForce(f_j, shf_c);
680 addShiftForce(f_k, shf_kj);
681 addShiftForce(f_l, shf_lj);
687 #undef NBLIB_ALWAYS_INLINE
689 #endif // NBLIB_LISTEDFORCES_KERNELS_HPP