2 * This file is part of the GROMACS molecular simulation package.
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.
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/vec.h"
51 #include "nblib/basicdefinitions.h"
52 #include "bondtypes.h"
57 /*! \brief kernel to calculate the scalar part for simple harmonic bond forces
60 * \param k spring constant
61 * \param x0 equilibrium distance
62 * \param x input bond length
64 * \return tuple<force, potential energy>
67 inline std::tuple<T, T> harmonicScalarForce(T k, T x0, T x)
73 T epot = 0.5 * k * dx2;
75 return std::make_tuple(force, epot);
77 /* That was 6 flops */
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
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
90 * \return tuple<force, potential energy, lambda-interpolated energy>
93 inline std::tuple<T, T, T> harmonicScalarForce(T kA, T kB, T xA, T xB, T x, T lambda)
95 // code unchanged relative to Gromacs
98 T kk = L1 * kA + lambda * kB;
99 T x0 = L1 * xA + lambda * xB;
105 T epot = 0.5 * kk * dx2;
106 T dvdlambda = 0.5 * (kB - kA) * dx2 + (xA - xB) * kk * dx;
108 return std::make_tuple(force, epot, dvdlambda);
110 /* That was 19 flops */
113 //! abstraction layer for different 2-center bonds
115 inline auto bondKernel(T dr, const HarmonicBondType& bond)
117 return harmonicScalarForce(bond.forceConstant(), bond.equilConstant(), dr);
121 /*! \brief kernel to calculate the scalar part for the forth power pontential bond forces
124 * \param k spring constant
125 * \param x0 squared equilibrium distance
126 * \param x squared input bond length
128 * \return tuple<force, potential energy>
131 inline std::tuple<T, T> g96ScalarForce(T k, T x0, T x)
136 T force = -k * dx * x;
137 T epot = 0.25 * k * dx2;
139 return std::make_tuple(force, epot);
141 /* That was 7 flops */
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
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
154 * \return tuple<force, potential energy, lambda-interpolated energy>
157 inline std::tuple<T, T, T> g96ScalarForce(T kA, T kB, T xA, T xB, T x, T lambda)
160 T kk = L1 * kA + lambda * kB;
161 T x0 = L1 * xA + lambda * xB;
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;
171 return std::make_tuple(force, epot, dvdlambda);
173 /* That was 21 flops */
176 //! Abstraction layer for different 2-center bonds. Forth power case
178 inline auto bondKernel(T dr, const G96BondType& bond)
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);
185 /*! \brief kernel to calculate the scalar part for the morse pontential bond forces
188 * \param k force constant
189 * \param beta beta exponent
190 * \param x0 equilibrium distance
191 * \param x input bond length
193 * \return tuple<force, potential energy>
196 inline std::tuple<T, T> morseScalarForce(T k, T beta, T x0, T x)
198 T exponent = std::exp(-beta * (x - x0)); /* 12 */
199 T omexp = 1.0 - exponent; /* 1 */
200 T kexp = k * omexp; /* 1 */
202 T epot = kexp * omexp; /* 1 */
203 T force = -2.0 * beta * exponent * omexp; /* 4 */
205 return std::make_tuple(force, epot);
207 /* That was 23 flops */
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
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
222 * \return tuple<force, potential energy, lambda-interpolated energy>
225 inline std::tuple<T, T, T> morseScalarForce(T kA, T kB, T betaA, T betaB, T xA, T xB, T x, T lambda)
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 */
232 T exponent = std::exp(-beta * (x - x0)); /* 12 */
233 T omexp = 1.0 - exponent; /* 1 */
234 T kexp = k * omexp; /* 1 */
236 T epot = kexp * omexp; /* 1 */
237 T force = -2.0 * beta * exponent * omexp; /* 4 */
239 T dvdlambda = (kB - kA) * omexp * omexp
240 - (2.0 - 2.0 * omexp) * omexp * k
241 * ((xB - xA) * beta - (betaB - betaA) * (x - x0)); /* 15 */
243 return std::make_tuple(force, epot, dvdlambda);
245 /* That was 44 flops */
248 //! Abstraction layer for different 2-center bonds. Morse case
250 inline auto bondKernel(T dr, const MorseBondType& bond)
252 return morseScalarForce(bond.forceConstant(), bond.exponent(), bond.equilDistance(), dr);
256 /*! \brief kernel to calculate the scalar part for the FENE pontential bond forces
259 * \param k spring constant
260 * \param x0 equilibrium distance
261 * \param x input bond length
263 * \return tuple<force, potential energy>
266 inline std::tuple<T, T> FENEScalarForce(T k, T x0, T x)
271 T omx2_ox02 = 1.0 - (x2 / x02);
273 T epot = -0.5 * k * x02 * std::log(omx2_ox02);
274 T force = -k * x / omx2_ox02;
276 return std::make_tuple(force, epot);
278 /* That was 24 flops */
281 // TODO: Implement the free energy version of FENE (finitely extensible nonlinear elastic) bond types
283 //! Abstraction layer for different 2-center bonds. FENE case
285 inline auto bondKernel(T dr, const FENEBondType& bond)
287 return FENEScalarForce(bond.forceConstant(), bond.equilConstant(), dr);
291 /*! \brief kernel to calculate the scalar part for cubic potential bond forces
294 * \param kc cubic spring constant
295 * \param kq quadratic spring constant
296 * \param x0 equilibrium distance
297 * \param x input bond length
299 * \return tuple<force, potential energy>
302 inline std::tuple<T, T> cubicScalarForce(T kc, T kq, T x0, T x)
307 T kdist2 = kdist * dx;
309 T epot = kdist2 + (kc * kdist2 * dx);
310 T force = -((2.0 * kdist) + (3.0 * kdist2 * kc));
312 return std::make_tuple(force, epot);
314 /* That was 16 flops */
317 // TODO: Implement the free energy version of Cubic bond types
320 inline auto bondKernel(T dr, const CubicBondType& bond)
322 return cubicScalarForce(bond.cubicForceConstant(), bond.quadraticForceConstant(), bond.equilDistance(), dr);
326 /*! \brief kernel to calculate the scalar part for half attractive potential bond forces
329 * \param k spring constant
330 * \param x0 equilibrium distance
331 * \param x input bond length
333 * \return tuple<force, potential energy>
336 inline std::tuple<T, T> halfAttractiveScalarForce(T k, T x0, T x)
343 T epot = -0.5 * k * dx4;
344 T force = -2.0 * k * dx3;
346 return std::make_tuple(force, epot);
348 /* That was 10 flops */
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
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
361 * \return tuple<force, potential energy, lambda-interpolated energy>
364 inline std::tuple<T, T, T> halfAttractiveScalarForce(T kA, T kB, T xA, T xB, T x, T lambda)
367 T kk = L1 * kA + lambda * kB;
368 T x0 = L1 * xA + lambda * xB;
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);
379 return std::make_tuple(force, epot, dvdlambda);
381 /* That was 29 flops */
385 inline auto bondKernel(T dr, const HalfAttractiveQuarticBondType& bond)
387 return halfAttractiveScalarForce(bond.forceConstant(), bond.equilConstant(), dr);
391 //! Three-center interaction type kernels
393 // linear angles go here
394 // quartic angles go here
396 //! Three-center interaction type dispatch
399 inline auto threeCenterKernel(T dr, const HarmonicAngle& angle)
401 return harmonicScalarForce(angle.forceConstant(), angle.equilConstant(), dr);
405 //! \brief Computes and returns the proper dihedral force
407 inline auto fourCenterKernel(T phi, const ProperDihedral& properDihedral)
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);
416 //! \brief Ensure that a geometric quantity lies in (-pi, pi)
417 static inline void makeAnglePeriodic(real& angle)
423 else if (angle < -M_PI)
429 /*! \brief calculate the cosine of the angle between a and b
431 * \tparam T float or double
432 * \param a a 3D vector
433 * \param b another 3D vector
434 * \return the cosine of the angle between a and b
436 * ax*bx + ay*by + az*bz
437 * cos-vec (a,b) = ---------------------
441 inline T basicVectorCosAngle(gmx::BasicVector<T> a, gmx::BasicVector<T> b)
443 gmx::BasicVector<double> a_double(a[0], a[1], a[2]);
444 gmx::BasicVector<double> b_double(b[0], b[1], b[2]);
446 double numerator = dot(a_double, b_double);
447 double denominatorSq = dot(a_double, a_double) * dot(b_double, b_double);
449 T cosval = (denominatorSq > 0) ? static_cast<T>(numerator * gmx::invsqrt(denominatorSq)) : 1;
450 cosval = std::min(cosval, T(1.0));
453 return std::max(cosval, T(-1.0));
456 //! \brief Computes and returns a dihedral phi angle
457 static inline real dihedralPhi(rvec dxIJ, rvec dxKJ, rvec dxKL, rvec m, rvec n)
459 cprod(dxIJ, dxKJ, m);
460 cprod(dxKJ, dxKL, n);
461 real phi = gmx_angle(m, n);
462 real ipr = iprod(dxIJ, n);
463 real sign = (ipr < 0.0) ? -1.0 : 1.0;
468 //! \brief Computes and returns the improper dihedral force
470 inline auto fourCenterKernel(T phi, const ImproperDihedral& improperDihedral)
472 T deltaPhi = phi - improperDihedral.equilConstant();
473 /* deltaPhi cannot be outside (-pi,pi) */
474 makeAnglePeriodic(deltaPhi);
475 const T force = -improperDihedral.forceConstant() * deltaPhi;
476 const T ePot = 0.5 * improperDihedral.forceConstant() * deltaPhi * deltaPhi;
477 return std::make_tuple(force, ePot);
480 //! \brief Computes and returns the Ryckaert-Belleman dihedral force
482 inline auto fourCenterKernel(T phi, const RyckaertBellemanDihedral& ryckaertBellemanDihedral)
484 /* Change to polymer convention */
485 const T localPhi = (phi < 0) ? (phi += M_PI) : (phi -= M_PI);
486 T cos_phi = std::cos(localPhi);
487 T ePot = ryckaertBellemanDihedral[0];
491 for (int i = 1; i < int(ryckaertBellemanDihedral.size()); i++)
493 force += ryckaertBellemanDihedral[i] * cosineFactor * i;
494 cosineFactor *= cos_phi;
495 ePot += cosineFactor * ryckaertBellemanDihedral[i];
497 /* Beware of accuracy loss, cannot use 1-sqrt(cos^2) ! */
498 force = -force * std::sin(localPhi);
499 return std::make_tuple(force, ePot);
502 //! Two-center category common
504 /*! \brief Spreads and accumulates the bonded forces to the two atoms and adds the virial contribution when needed
506 * \p shiftIndex is used as the periodic shift.
509 inline void spreadTwoCenterForces(const T bondForce,
514 for (int m = 0; m < dimSize; m++) /* 15 */
516 const T fij = bondForce * dx[m];
517 (*force_i)[m] += fij;
518 (*force_j)[m] -= fij;
522 //! Three-center category common
524 /*! \brief spread force to 3 centers based on scalar force and angle
536 inline void spreadThreeCenterForces(T cos_theta,
538 const gmx::RVec& r_ij,
539 const gmx::RVec& r_kj,
544 T cos_theta2 = cos_theta * cos_theta;
547 T st = force / std::sqrt(1 - cos_theta2); /* 12 */
548 T sth = st * cos_theta; /* 1 */
549 T nrij2 = dot(r_ij, r_ij); /* 5 */
550 T nrkj2 = dot(r_kj, r_kj); /* 5 */
552 T nrij_1 = 1.0 / std::sqrt(nrij2); /* 10 */
553 T nrkj_1 = 1.0 / std::sqrt(nrkj2); /* 10 */
555 T cik = st * nrij_1 * nrkj_1; /* 2 */
556 T cii = sth * nrij_1 * nrij_1; /* 2 */
557 T ckk = sth * nrkj_1 * nrkj_1; /* 2 */
559 gmx::RVec f_i{0, 0, 0};
560 gmx::RVec f_j{0, 0, 0};
561 gmx::RVec f_k{0, 0, 0};
562 for (int m = 0; m < dimSize; m++)
564 f_i[m] = -(cik * r_kj[m] - cii * r_ij[m]);
565 f_k[m] = -(cik * r_ij[m] - ckk * r_kj[m]);
566 f_j[m] = -f_i[m] - f_k[m];
567 (*force_i)[m] += f_i[m];
568 (*force_j)[m] += f_j[m];
569 (*force_k)[m] += f_k[m];
574 //! Four-center category common
576 inline void spreadFourCenterForces(T force, rvec dxIJ, rvec dxJK, rvec dxKL, rvec m, rvec n,
582 rvec f_i, f_j, f_k, f_l;
583 rvec uvec, vvec, svec;
584 T iprm = iprod(m, m); /* 5 */
585 T iprn = iprod(n, n); /* 5 */
586 T nrkj2 = iprod(dxJK, dxJK); /* 5 */
587 T toler = nrkj2 * GMX_REAL_EPS;
588 if ((iprm > toler) && (iprn > toler))
590 T nrkj_1 = gmx::invsqrt(nrkj2); /* 10 */
591 T nrkj_2 = nrkj_1 * nrkj_1; /* 1 */
592 T nrkj = nrkj2 * nrkj_1; /* 1 */
593 T a = -force * nrkj / iprm; /* 11 */
594 svmul(a, m, f_i); /* 3 */
595 T b = force * nrkj / iprn; /* 11 */
596 svmul(b, n, f_l); /* 3 */
597 T p = iprod(dxIJ, dxJK); /* 5 */
599 T q = iprod(dxKL, dxJK); /* 5 */
601 svmul(p, f_i, uvec); /* 3 */
602 svmul(q, f_l, vvec); /* 3 */
603 rvec_sub(uvec, vvec, svec); /* 3 */
604 rvec_sub(f_i, svec, f_j); /* 3 */
605 rvec_add(f_l, svec, f_k); /* 3 */
606 rvec_inc(force_i->as_vec(), f_i); /* 3 */
607 rvec_dec(force_j->as_vec(), f_j); /* 3 */
608 rvec_dec(force_k->as_vec(), f_k); /* 3 */
609 rvec_inc(force_l->as_vec(), f_l); /* 3 */
614 #endif // NBLIB_LISTEDFORCES_KERNELS_HPP