2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
7 * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
40 * \brief This file defines low-level functions necessary for
41 * computing energies and forces for listed interactions.
43 * \author Mark Abraham <mark.j.abraham@gmail.com>
45 * \ingroup module_listed_forces
58 #include "gromacs/gmxlib/nrnb.h"
59 #include "gromacs/listed_forces/disre.h"
60 #include "gromacs/listed_forces/orires.h"
61 #include "gromacs/math/functions.h"
62 #include "gromacs/math/units.h"
63 #include "gromacs/math/utilities.h"
64 #include "gromacs/math/vec.h"
65 #include "gromacs/mdtypes/fcdata.h"
66 #include "gromacs/mdtypes/mdatom.h"
67 #include "gromacs/pbcutil/ishift.h"
68 #include "gromacs/pbcutil/pbc.h"
69 #include "gromacs/pbcutil/pbc_simd.h"
70 #include "gromacs/simd/simd.h"
71 #include "gromacs/simd/simd_math.h"
72 #include "gromacs/simd/vector_operations.h"
73 #include "gromacs/utility/basedefinitions.h"
74 #include "gromacs/utility/enumerationhelpers.h"
75 #include "gromacs/utility/fatalerror.h"
76 #include "gromacs/utility/real.h"
77 #include "gromacs/utility/smalloc.h"
79 #include "listed_internal.h"
82 using namespace gmx; // TODO: Remove when this file is moved into gmx namespace
84 const EnumerationArray<BondedKernelFlavor, std::string> c_bondedKernelFlavorStrings = {
85 "forces, using SIMD when available", "forces, not using SIMD",
86 "forces, virial, and energy (ie. not using SIMD)", "forces and energy (ie. not using SIMD)"
91 //! Type of CPU function to compute a bonded interaction.
92 using BondedFunction = real (*)(int nbonds,
93 const t_iatom iatoms[],
94 const t_iparams iparams[],
105 /*! \brief Mysterious CMAP coefficient matrix */
106 const int cmap_coeff_matrix[] = {
107 1, 0, -3, 2, 0, 0, 0, 0, -3, 0, 9, -6, 2, 0, -6, 4, 0, 0, 0, 0, 0, 0, 0, 0,
108 3, 0, -9, 6, -2, 0, 6, -4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 9, -6, 0, 0, -6, 4,
109 0, 0, 3, -2, 0, 0, 0, 0, 0, 0, -9, 6, 0, 0, 6, -4, 0, 0, 0, 0, 1, 0, -3, 2,
110 -2, 0, 6, -4, 1, 0, -3, 2, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 3, -2, 1, 0, -3, 2,
111 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 2, 0, 0, 3, -2, 0, 0, 0, 0, 0, 0, 3, -2,
112 0, 0, -6, 4, 0, 0, 3, -2, 0, 1, -2, 1, 0, 0, 0, 0, 0, -3, 6, -3, 0, 2, -4, 2,
113 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, -6, 3, 0, -2, 4, -2, 0, 0, 0, 0, 0, 0, 0, 0,
114 0, 0, -3, 3, 0, 0, 2, -2, 0, 0, -1, 1, 0, 0, 0, 0, 0, 0, 3, -3, 0, 0, -2, 2,
115 0, 0, 0, 0, 0, 1, -2, 1, 0, -2, 4, -2, 0, 1, -2, 1, 0, 0, 0, 0, 0, 0, 0, 0,
116 0, -1, 2, -1, 0, 1, -2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, -1, 1,
117 0, 0, 0, 0, 0, 0, -1, 1, 0, 0, 2, -2, 0, 0, -1, 1
121 /*! \brief Compute dx = xi - xj, modulo PBC if non-NULL
123 * \todo This kind of code appears in many places. Consolidate it */
124 int pbc_rvec_sub(const t_pbc* pbc, const rvec xi, const rvec xj, rvec dx)
128 return pbc_dx_aiuc(pbc, xi, xj, dx);
132 rvec_sub(xi, xj, dx);
140 real bond_angle(const rvec xi,
149 /* Return value is the angle between the bonds i-j and j-k */
154 *t1 = pbc_rvec_sub(pbc, xi, xj, r_ij); /* 3 */
155 *t2 = pbc_rvec_sub(pbc, xk, xj, r_kj); /* 3 */
157 *costh = cos_angle(r_ij, r_kj); /* 25 */
158 th = std::acos(*costh); /* 10 */
163 real dih_angle(const rvec xi,
177 *t1 = pbc_rvec_sub(pbc, xi, xj, r_ij); /* 3 */
178 *t2 = pbc_rvec_sub(pbc, xk, xj, r_kj); /* 3 */
179 *t3 = pbc_rvec_sub(pbc, xk, xl, r_kl); /* 3 */
181 cprod(r_ij, r_kj, m); /* 9 */
182 cprod(r_kj, r_kl, n); /* 9 */
183 real phi = gmx_angle(m, n); /* 49 (assuming 25 for atan2) */
184 real ipr = iprod(r_ij, n); /* 5 */
185 real sign = (ipr < 0.0) ? -1.0 : 1.0;
186 phi = sign * phi; /* 1 */
192 void make_dp_periodic(real* dp) /* 1 flop? */
194 /* dp cannot be outside (-pi,pi) */
199 else if (*dp < -M_PI)
208 /*! \brief Spreads and accumulates the bonded forces to the two atoms and adds the virial contribution when needed
210 * \p shiftIndex is used as the periodic shift.
212 template<BondedKernelFlavor flavor>
213 inline void spreadBondForces(const real bondForce,
221 for (int m = 0; m < DIM; m++) /* 15 */
223 const real fij = bondForce * dx[m];
226 if (computeVirial(flavor))
228 fshift[shiftIndex][m] += fij;
229 fshift[CENTRAL][m] -= fij;
234 /*! \brief Morse potential bond
236 * By Frank Everdij. Three parameters needed:
238 * b0 = equilibrium distance in nm
239 * be = beta in nm^-1 (actually, it's nu_e*Sqrt(2*pi*pi*mu/D_e))
240 * cb = well depth in kJ/mol
242 * Note: the potential is referenced to be +cb at infinite separation
243 * and zero at the equilibrium distance!
245 template<BondedKernelFlavor flavor>
246 real morse_bonds(int nbonds,
247 const t_iatom forceatoms[],
248 const t_iparams forceparams[],
255 const t_mdatoms gmx_unused* md,
256 t_fcdata gmx_unused* fcd,
257 int gmx_unused* global_atom_index)
259 const real one = 1.0;
260 const real two = 2.0;
261 real dr, dr2, temp, omtemp, cbomtemp, fbond, vbond, vtot;
262 real b0, be, cb, b0A, beA, cbA, b0B, beB, cbB, L1;
264 int i, ki, type, ai, aj;
267 for (i = 0; (i < nbonds);)
269 type = forceatoms[i++];
270 ai = forceatoms[i++];
271 aj = forceatoms[i++];
273 b0A = forceparams[type].morse.b0A;
274 beA = forceparams[type].morse.betaA;
275 cbA = forceparams[type].morse.cbA;
277 b0B = forceparams[type].morse.b0B;
278 beB = forceparams[type].morse.betaB;
279 cbB = forceparams[type].morse.cbB;
281 L1 = one - lambda; /* 1 */
282 b0 = L1 * b0A + lambda * b0B; /* 3 */
283 be = L1 * beA + lambda * beB; /* 3 */
284 cb = L1 * cbA + lambda * cbB; /* 3 */
286 ki = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /* 3 */
287 dr2 = iprod(dx, dx); /* 5 */
288 dr = dr2 * gmx::invsqrt(dr2); /* 10 */
289 temp = std::exp(-be * (dr - b0)); /* 12 */
293 /* bonds are constrainted. This may _not_ include bond constraints if they are lambda dependent */
294 *dvdlambda += cbB - cbA;
298 omtemp = one - temp; /* 1 */
299 cbomtemp = cb * omtemp; /* 1 */
300 vbond = cbomtemp * omtemp; /* 1 */
301 fbond = -two * be * temp * cbomtemp * gmx::invsqrt(dr2); /* 9 */
302 vtot += vbond; /* 1 */
304 *dvdlambda += (cbB - cbA) * omtemp * omtemp
305 - (2 - 2 * omtemp) * omtemp * cb * ((b0B - b0A) * be - (beB - beA) * (dr - b0)); /* 15 */
307 spreadBondForces<flavor>(fbond, dx, ai, aj, f, ki, fshift); /* 15 */
313 template<BondedKernelFlavor flavor>
314 real cubic_bonds(int nbonds,
315 const t_iatom forceatoms[],
316 const t_iparams forceparams[],
321 real gmx_unused lambda,
322 real gmx_unused* dvdlambda,
323 const t_mdatoms gmx_unused* md,
324 t_fcdata gmx_unused* fcd,
325 int gmx_unused* global_atom_index)
327 const real three = 3.0;
328 const real two = 2.0;
330 real dr, dr2, dist, kdist, kdist2, fbond, vbond, vtot;
332 int i, ki, type, ai, aj;
335 for (i = 0; (i < nbonds);)
337 type = forceatoms[i++];
338 ai = forceatoms[i++];
339 aj = forceatoms[i++];
341 b0 = forceparams[type].cubic.b0;
342 kb = forceparams[type].cubic.kb;
343 kcub = forceparams[type].cubic.kcub;
345 ki = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /* 3 */
346 dr2 = iprod(dx, dx); /* 5 */
353 dr = dr2 * gmx::invsqrt(dr2); /* 10 */
356 kdist2 = kdist * dist;
358 vbond = kdist2 + kcub * kdist2 * dist;
359 fbond = -(two * kdist + three * kdist2 * kcub) / dr;
361 vtot += vbond; /* 21 */
363 spreadBondForces<flavor>(fbond, dx, ai, aj, f, ki, fshift); /* 15 */
368 template<BondedKernelFlavor flavor>
369 real FENE_bonds(int nbonds,
370 const t_iatom forceatoms[],
371 const t_iparams forceparams[],
376 real gmx_unused lambda,
377 real gmx_unused* dvdlambda,
378 const t_mdatoms gmx_unused* md,
379 t_fcdata gmx_unused* fcd,
380 int* global_atom_index)
382 const real half = 0.5;
383 const real one = 1.0;
385 real dr2, bm2, omdr2obm2, fbond, vbond, vtot;
387 int i, ki, type, ai, aj;
390 for (i = 0; (i < nbonds);)
392 type = forceatoms[i++];
393 ai = forceatoms[i++];
394 aj = forceatoms[i++];
396 bm = forceparams[type].fene.bm;
397 kb = forceparams[type].fene.kb;
399 ki = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /* 3 */
400 dr2 = iprod(dx, dx); /* 5 */
411 gmx_fatal(FARGS, "r^2 (%f) >= bm^2 (%f) in FENE bond between atoms %d and %d", dr2, bm2,
412 glatnr(global_atom_index, ai), glatnr(global_atom_index, aj));
415 omdr2obm2 = one - dr2 / bm2;
417 vbond = -half * kb * bm2 * std::log(omdr2obm2);
418 fbond = -kb / omdr2obm2;
420 vtot += vbond; /* 35 */
422 spreadBondForces<flavor>(fbond, dx, ai, aj, f, ki, fshift); /* 15 */
427 real harmonic(real kA, real kB, real xA, real xB, real x, real lambda, real* V, real* F)
429 const real half = 0.5;
430 real L1, kk, x0, dx, dx2;
431 real v, f, dvdlambda;
434 kk = L1 * kA + lambda * kB;
435 x0 = L1 * xA + lambda * xB;
442 dvdlambda = half * (kB - kA) * dx2 + (xA - xB) * kk * dx;
449 /* That was 19 flops */
452 template<BondedKernelFlavor flavor>
453 std::enable_if_t<flavor != BondedKernelFlavor::ForcesSimdWhenAvailable || !GMX_SIMD_HAVE_REAL, real>
455 const t_iatom forceatoms[],
456 const t_iparams forceparams[],
463 const t_mdatoms gmx_unused* md,
464 t_fcdata gmx_unused* fcd,
465 int gmx_unused* global_atom_index)
467 int i, ki, ai, aj, type;
468 real dr, dr2, fbond, vbond, vtot;
472 for (i = 0; (i < nbonds);)
474 type = forceatoms[i++];
475 ai = forceatoms[i++];
476 aj = forceatoms[i++];
478 ki = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /* 3 */
479 dr2 = iprod(dx, dx); /* 5 */
480 dr = std::sqrt(dr2); /* 10 */
482 *dvdlambda += harmonic(forceparams[type].harmonic.krA, forceparams[type].harmonic.krB,
483 forceparams[type].harmonic.rA, forceparams[type].harmonic.rB, dr,
484 lambda, &vbond, &fbond); /* 19 */
492 vtot += vbond; /* 1*/
493 fbond *= gmx::invsqrt(dr2); /* 6 */
495 spreadBondForces<flavor>(fbond, dx, ai, aj, f, ki, fshift); /* 15 */
500 #if GMX_SIMD_HAVE_REAL
502 /*! \brief Computes forces (only) for harmonic bonds using SIMD intrinsics
504 * As plain-C bonds(), but using SIMD to calculate many bonds at once.
505 * This routines does not calculate energies and shift forces.
507 template<BondedKernelFlavor flavor>
508 std::enable_if_t<flavor == BondedKernelFlavor::ForcesSimdWhenAvailable, real>
510 const t_iatom forceatoms[],
511 const t_iparams forceparams[],
514 rvec gmx_unused fshift[],
516 real gmx_unused lambda,
517 real gmx_unused* dvdlambda,
518 const t_mdatoms gmx_unused* md,
519 t_fcdata gmx_unused* fcd,
520 int gmx_unused* global_atom_index)
522 constexpr int nfa1 = 3;
523 alignas(GMX_SIMD_ALIGNMENT) std::int32_t ai[GMX_SIMD_REAL_WIDTH];
524 alignas(GMX_SIMD_ALIGNMENT) std::int32_t aj[GMX_SIMD_REAL_WIDTH];
525 alignas(GMX_SIMD_ALIGNMENT) real coeff[2 * GMX_SIMD_REAL_WIDTH];
527 alignas(GMX_SIMD_ALIGNMENT) real pbc_simd[9 * GMX_SIMD_REAL_WIDTH];
529 set_pbc_simd(pbc, pbc_simd);
531 const SimdReal real_eps = SimdReal(GMX_REAL_EPS);
533 /* nbonds is the number of bonds times nfa1, here we step GMX_SIMD_REAL_WIDTH angles */
534 for (int i = 0; i < nbonds; i += GMX_SIMD_REAL_WIDTH * nfa1)
536 /* Collect atoms for GMX_SIMD_REAL_WIDTH angles.
537 * iu indexes into forceatoms, we should not let iu go beyond nbonds.
540 for (int s = 0; s < GMX_SIMD_REAL_WIDTH; s++)
542 const int type = forceatoms[iu];
543 ai[s] = forceatoms[iu + 1];
544 aj[s] = forceatoms[iu + 2];
546 /* At the end fill the arrays with the last atoms and 0 params */
547 if (i + s * nfa1 < nbonds)
549 coeff[s] = forceparams[type].harmonic.krA;
550 coeff[GMX_SIMD_REAL_WIDTH + s] = forceparams[type].harmonic.rA;
552 if (iu + nfa1 < nbonds)
560 coeff[GMX_SIMD_REAL_WIDTH + s] = 0;
564 /* Store the non PBC corrected distances packed and aligned */
567 gatherLoadUTranspose<3>(reinterpret_cast<const real*>(x), ai, &xi, &yi, &zi);
568 gatherLoadUTranspose<3>(reinterpret_cast<const real*>(x), aj, &xj, &yj, &zj);
569 SimdReal rij_x = xi - xj;
570 SimdReal rij_y = yi - yj;
571 SimdReal rij_z = zi - zj;
573 pbc_correct_dx_simd(&rij_x, &rij_y, &rij_z, pbc_simd);
575 const SimdReal dist2 = rij_x * rij_x + rij_y * rij_y + rij_z * rij_z;
576 // Here we avoid sqrt(0), the force will be zero because rij=0
577 const SimdReal invDist = invsqrt(max(dist2, real_eps));
579 const SimdReal k = load<SimdReal>(coeff);
580 const SimdReal r0 = load<SimdReal>(coeff + GMX_SIMD_REAL_WIDTH);
582 // Compute the force divided by the distance
583 const SimdReal forceOverR = -k * (dist2 * invDist - r0) * invDist;
585 const SimdReal f_x = forceOverR * rij_x;
586 const SimdReal f_y = forceOverR * rij_y;
587 const SimdReal f_z = forceOverR * rij_z;
589 transposeScatterIncrU<4>(reinterpret_cast<real*>(f), ai, f_x, f_y, f_z);
590 transposeScatterDecrU<4>(reinterpret_cast<real*>(f), aj, f_x, f_y, f_z);
596 #endif // GMX_SIMD_HAVE_REAL
598 template<BondedKernelFlavor flavor>
599 real restraint_bonds(int nbonds,
600 const t_iatom forceatoms[],
601 const t_iparams forceparams[],
608 const t_mdatoms gmx_unused* md,
609 t_fcdata gmx_unused* fcd,
610 int gmx_unused* global_atom_index)
612 int i, ki, ai, aj, type;
613 real dr, dr2, fbond, vbond, vtot;
615 real low, dlow, up1, dup1, up2, dup2, k, dk;
622 for (i = 0; (i < nbonds);)
624 type = forceatoms[i++];
625 ai = forceatoms[i++];
626 aj = forceatoms[i++];
628 ki = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /* 3 */
629 dr2 = iprod(dx, dx); /* 5 */
630 dr = dr2 * gmx::invsqrt(dr2); /* 10 */
632 low = L1 * forceparams[type].restraint.lowA + lambda * forceparams[type].restraint.lowB;
633 dlow = -forceparams[type].restraint.lowA + forceparams[type].restraint.lowB;
634 up1 = L1 * forceparams[type].restraint.up1A + lambda * forceparams[type].restraint.up1B;
635 dup1 = -forceparams[type].restraint.up1A + forceparams[type].restraint.up1B;
636 up2 = L1 * forceparams[type].restraint.up2A + lambda * forceparams[type].restraint.up2B;
637 dup2 = -forceparams[type].restraint.up2A + forceparams[type].restraint.up2B;
638 k = L1 * forceparams[type].restraint.kA + lambda * forceparams[type].restraint.kB;
639 dk = -forceparams[type].restraint.kA + forceparams[type].restraint.kB;
646 vbond = 0.5 * k * drh2;
648 *dvdlambda += 0.5 * dk * drh2 - k * dlow * drh;
659 vbond = 0.5 * k * drh2;
661 *dvdlambda += 0.5 * dk * drh2 - k * dup1 * drh;
666 vbond = k * (up2 - up1) * (0.5 * (up2 - up1) + drh);
667 fbond = -k * (up2 - up1);
668 *dvdlambda += dk * (up2 - up1) * (0.5 * (up2 - up1) + drh)
669 + k * (dup2 - dup1) * (up2 - up1 + drh) - k * (up2 - up1) * dup2;
677 vtot += vbond; /* 1*/
678 fbond *= gmx::invsqrt(dr2); /* 6 */
680 spreadBondForces<flavor>(fbond, dx, ai, aj, f, ki, fshift); /* 15 */
686 template<BondedKernelFlavor flavor>
687 real polarize(int nbonds,
688 const t_iatom forceatoms[],
689 const t_iparams forceparams[],
697 t_fcdata gmx_unused* fcd,
698 int gmx_unused* global_atom_index)
700 int i, ki, ai, aj, type;
701 real dr, dr2, fbond, vbond, vtot, ksh;
705 for (i = 0; (i < nbonds);)
707 type = forceatoms[i++];
708 ai = forceatoms[i++];
709 aj = forceatoms[i++];
710 ksh = gmx::square(md->chargeA[aj]) * ONE_4PI_EPS0 / forceparams[type].polarize.alpha;
712 ki = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /* 3 */
713 dr2 = iprod(dx, dx); /* 5 */
714 dr = std::sqrt(dr2); /* 10 */
716 *dvdlambda += harmonic(ksh, ksh, 0, 0, dr, lambda, &vbond, &fbond); /* 19 */
723 vtot += vbond; /* 1*/
724 fbond *= gmx::invsqrt(dr2); /* 6 */
726 spreadBondForces<flavor>(fbond, dx, ai, aj, f, ki, fshift); /* 15 */
731 template<BondedKernelFlavor flavor>
732 real anharm_polarize(int nbonds,
733 const t_iatom forceatoms[],
734 const t_iparams forceparams[],
742 t_fcdata gmx_unused* fcd,
743 int gmx_unused* global_atom_index)
745 int i, ki, ai, aj, type;
746 real dr, dr2, fbond, vbond, vtot, ksh, khyp, drcut, ddr, ddr3;
750 for (i = 0; (i < nbonds);)
752 type = forceatoms[i++];
753 ai = forceatoms[i++];
754 aj = forceatoms[i++];
755 ksh = gmx::square(md->chargeA[aj]) * ONE_4PI_EPS0 / forceparams[type].anharm_polarize.alpha; /* 7*/
756 khyp = forceparams[type].anharm_polarize.khyp;
757 drcut = forceparams[type].anharm_polarize.drcut;
759 ki = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /* 3 */
760 dr2 = iprod(dx, dx); /* 5 */
761 dr = dr2 * gmx::invsqrt(dr2); /* 10 */
763 *dvdlambda += harmonic(ksh, ksh, 0, 0, dr, lambda, &vbond, &fbond); /* 19 */
773 ddr3 = ddr * ddr * ddr;
774 vbond += khyp * ddr * ddr3;
775 fbond -= 4 * khyp * ddr3;
777 fbond *= gmx::invsqrt(dr2); /* 6 */
778 vtot += vbond; /* 1*/
780 spreadBondForces<flavor>(fbond, dx, ai, aj, f, ki, fshift); /* 15 */
785 template<BondedKernelFlavor flavor>
786 real water_pol(int nbonds,
787 const t_iatom forceatoms[],
788 const t_iparams forceparams[],
791 rvec gmx_unused fshift[],
792 const t_pbc gmx_unused* pbc,
793 real gmx_unused lambda,
794 real gmx_unused* dvdlambda,
795 const t_mdatoms gmx_unused* md,
796 t_fcdata gmx_unused* fcd,
797 int gmx_unused* global_atom_index)
799 /* This routine implements anisotropic polarizibility for water, through
800 * a shell connected to a dummy with spring constant that differ in the
801 * three spatial dimensions in the molecular frame.
803 int i, m, aO, aH1, aH2, aD, aS, type, type0, ki;
804 rvec dOH1, dOH2, dHH, dOD, dDS, nW, kk, dx, kdx, proj;
805 real vtot, fij, r_HH, r_OD, r_nW, tx, ty, tz, qS;
810 type0 = forceatoms[0];
812 qS = md->chargeA[aS];
813 kk[XX] = gmx::square(qS) * ONE_4PI_EPS0 / forceparams[type0].wpol.al_x;
814 kk[YY] = gmx::square(qS) * ONE_4PI_EPS0 / forceparams[type0].wpol.al_y;
815 kk[ZZ] = gmx::square(qS) * ONE_4PI_EPS0 / forceparams[type0].wpol.al_z;
816 r_HH = 1.0 / forceparams[type0].wpol.rHH;
817 for (i = 0; (i < nbonds); i += 6)
819 type = forceatoms[i];
822 gmx_fatal(FARGS, "Sorry, type = %d, type0 = %d, file = %s, line = %d", type, type0,
825 aO = forceatoms[i + 1];
826 aH1 = forceatoms[i + 2];
827 aH2 = forceatoms[i + 3];
828 aD = forceatoms[i + 4];
829 aS = forceatoms[i + 5];
831 /* Compute vectors describing the water frame */
832 pbc_rvec_sub(pbc, x[aH1], x[aO], dOH1);
833 pbc_rvec_sub(pbc, x[aH2], x[aO], dOH2);
834 pbc_rvec_sub(pbc, x[aH2], x[aH1], dHH);
835 pbc_rvec_sub(pbc, x[aD], x[aO], dOD);
836 ki = pbc_rvec_sub(pbc, x[aS], x[aD], dDS);
837 cprod(dOH1, dOH2, nW);
839 /* Compute inverse length of normal vector
840 * (this one could be precomputed, but I'm too lazy now)
842 r_nW = gmx::invsqrt(iprod(nW, nW));
843 /* This is for precision, but does not make a big difference,
846 r_OD = gmx::invsqrt(iprod(dOD, dOD));
848 /* Normalize the vectors in the water frame */
850 svmul(r_HH, dHH, dHH);
851 svmul(r_OD, dOD, dOD);
853 /* Compute displacement of shell along components of the vector */
854 dx[ZZ] = iprod(dDS, dOD);
855 /* Compute projection on the XY plane: dDS - dx[ZZ]*dOD */
856 for (m = 0; (m < DIM); m++)
858 proj[m] = dDS[m] - dx[ZZ] * dOD[m];
861 /*dx[XX] = iprod(dDS,nW);
862 dx[YY] = iprod(dDS,dHH);*/
863 dx[XX] = iprod(proj, nW);
864 for (m = 0; (m < DIM); m++)
866 proj[m] -= dx[XX] * nW[m];
868 dx[YY] = iprod(proj, dHH);
869 /* Now compute the forces and energy */
870 kdx[XX] = kk[XX] * dx[XX];
871 kdx[YY] = kk[YY] * dx[YY];
872 kdx[ZZ] = kk[ZZ] * dx[ZZ];
873 vtot += iprod(dx, kdx);
875 for (m = 0; (m < DIM); m++)
877 /* This is a tensor operation but written out for speed */
878 tx = nW[m] * kdx[XX];
879 ty = dHH[m] * kdx[YY];
880 tz = dOD[m] * kdx[ZZ];
884 if (computeVirial(flavor))
886 fshift[ki][m] += fij;
887 fshift[CENTRAL][m] -= fij;
895 template<BondedKernelFlavor flavor>
897 do_1_thole(const rvec xi, const rvec xj, rvec fi, rvec fj, const t_pbc* pbc, real qq, rvec fshift[], real afac)
900 real r12sq, r12_1, r12bar, v0, v1, fscal, ebar, fff;
903 t = pbc_rvec_sub(pbc, xi, xj, r12); /* 3 */
905 r12sq = iprod(r12, r12); /* 5 */
906 r12_1 = gmx::invsqrt(r12sq); /* 5 */
907 r12bar = afac / r12_1; /* 5 */
908 v0 = qq * ONE_4PI_EPS0 * r12_1; /* 2 */
909 ebar = std::exp(-r12bar); /* 5 */
910 v1 = (1 - (1 + 0.5 * r12bar) * ebar); /* 4 */
911 fscal = ((v0 * r12_1) * v1 - v0 * 0.5 * afac * ebar * (r12bar + 1)) * r12_1; /* 9 */
913 for (m = 0; (m < DIM); m++)
915 fff = fscal * r12[m];
918 if (computeVirial(flavor))
921 fshift[CENTRAL][m] -= fff;
925 return v0 * v1; /* 1 */
929 template<BondedKernelFlavor flavor>
930 real thole_pol(int nbonds,
931 const t_iatom forceatoms[],
932 const t_iparams forceparams[],
937 real gmx_unused lambda,
938 real gmx_unused* dvdlambda,
940 t_fcdata gmx_unused* fcd,
941 int gmx_unused* global_atom_index)
943 /* Interaction between two pairs of particles with opposite charge */
944 int i, type, a1, da1, a2, da2;
945 real q1, q2, qq, a, al1, al2, afac;
948 for (i = 0; (i < nbonds);)
950 type = forceatoms[i++];
951 a1 = forceatoms[i++];
952 da1 = forceatoms[i++];
953 a2 = forceatoms[i++];
954 da2 = forceatoms[i++];
955 q1 = md->chargeA[da1];
956 q2 = md->chargeA[da2];
957 a = forceparams[type].thole.a;
958 al1 = forceparams[type].thole.alpha1;
959 al2 = forceparams[type].thole.alpha2;
961 afac = a * gmx::invsixthroot(al1 * al2);
962 V += do_1_thole<flavor>(x[a1], x[a2], f[a1], f[a2], pbc, qq, fshift, afac);
963 V += do_1_thole<flavor>(x[da1], x[a2], f[da1], f[a2], pbc, -qq, fshift, afac);
964 V += do_1_thole<flavor>(x[a1], x[da2], f[a1], f[da2], pbc, -qq, fshift, afac);
965 V += do_1_thole<flavor>(x[da1], x[da2], f[da1], f[da2], pbc, qq, fshift, afac);
971 // Avoid gcc 386 -O3 code generation bug in this function (see Issue
972 // #3205 for more information)
973 #if defined(__GNUC__) && defined(__i386__) && defined(__OPTIMIZE__)
974 # pragma GCC push_options
975 # pragma GCC optimize("O1")
976 # define avoid_gcc_i386_o3_code_generation_bug
979 template<BondedKernelFlavor flavor>
980 std::enable_if_t<flavor != BondedKernelFlavor::ForcesSimdWhenAvailable || !GMX_SIMD_HAVE_REAL, real>
982 const t_iatom forceatoms[],
983 const t_iparams forceparams[],
990 const t_mdatoms gmx_unused* md,
991 t_fcdata gmx_unused* fcd,
992 int gmx_unused* global_atom_index)
994 int i, ai, aj, ak, t1, t2, type;
996 real cos_theta, cos_theta2, theta, dVdt, va, vtot;
999 for (i = 0; i < nbonds;)
1001 type = forceatoms[i++];
1002 ai = forceatoms[i++];
1003 aj = forceatoms[i++];
1004 ak = forceatoms[i++];
1006 theta = bond_angle(x[ai], x[aj], x[ak], pbc, r_ij, r_kj, &cos_theta, &t1, &t2); /* 41 */
1008 *dvdlambda += harmonic(forceparams[type].harmonic.krA, forceparams[type].harmonic.krB,
1009 forceparams[type].harmonic.rA * DEG2RAD,
1010 forceparams[type].harmonic.rB * DEG2RAD, theta, lambda, &va, &dVdt); /* 21 */
1013 cos_theta2 = gmx::square(cos_theta);
1020 real nrkj_1, nrij_1;
1023 st = dVdt * gmx::invsqrt(1 - cos_theta2); /* 12 */
1024 sth = st * cos_theta; /* 1 */
1025 nrij2 = iprod(r_ij, r_ij); /* 5 */
1026 nrkj2 = iprod(r_kj, r_kj); /* 5 */
1028 nrij_1 = gmx::invsqrt(nrij2); /* 10 */
1029 nrkj_1 = gmx::invsqrt(nrkj2); /* 10 */
1031 cik = st * nrij_1 * nrkj_1; /* 2 */
1032 cii = sth * nrij_1 * nrij_1; /* 2 */
1033 ckk = sth * nrkj_1 * nrkj_1; /* 2 */
1035 for (m = 0; m < DIM; m++)
1037 f_i[m] = -(cik * r_kj[m] - cii * r_ij[m]);
1038 f_k[m] = -(cik * r_ij[m] - ckk * r_kj[m]);
1039 f_j[m] = -f_i[m] - f_k[m];
1044 if (computeVirial(flavor))
1046 rvec_inc(fshift[t1], f_i);
1047 rvec_inc(fshift[CENTRAL], f_j);
1048 rvec_inc(fshift[t2], f_k);
1056 #ifdef avoid_gcc_i386_o3_code_generation_bug
1057 # pragma GCC pop_options
1058 # undef avoid_gcc_i386_o3_code_generation_bug
1061 #if GMX_SIMD_HAVE_REAL
1063 /* As angles, but using SIMD to calculate many angles at once.
1064 * This routines does not calculate energies and shift forces.
1066 template<BondedKernelFlavor flavor>
1067 std::enable_if_t<flavor == BondedKernelFlavor::ForcesSimdWhenAvailable, real>
1069 const t_iatom forceatoms[],
1070 const t_iparams forceparams[],
1073 rvec gmx_unused fshift[],
1075 real gmx_unused lambda,
1076 real gmx_unused* dvdlambda,
1077 const t_mdatoms gmx_unused* md,
1078 t_fcdata gmx_unused* fcd,
1079 int gmx_unused* global_atom_index)
1084 alignas(GMX_SIMD_ALIGNMENT) std::int32_t ai[GMX_SIMD_REAL_WIDTH];
1085 alignas(GMX_SIMD_ALIGNMENT) std::int32_t aj[GMX_SIMD_REAL_WIDTH];
1086 alignas(GMX_SIMD_ALIGNMENT) std::int32_t ak[GMX_SIMD_REAL_WIDTH];
1087 alignas(GMX_SIMD_ALIGNMENT) real coeff[2 * GMX_SIMD_REAL_WIDTH];
1088 SimdReal deg2rad_S(DEG2RAD);
1089 SimdReal xi_S, yi_S, zi_S;
1090 SimdReal xj_S, yj_S, zj_S;
1091 SimdReal xk_S, yk_S, zk_S;
1092 SimdReal k_S, theta0_S;
1093 SimdReal rijx_S, rijy_S, rijz_S;
1094 SimdReal rkjx_S, rkjy_S, rkjz_S;
1095 SimdReal one_S(1.0);
1096 SimdReal one_min_eps_S(1.0_real - GMX_REAL_EPS); // Largest number < 1
1099 SimdReal nrij2_S, nrij_1_S;
1100 SimdReal nrkj2_S, nrkj_1_S;
1101 SimdReal cos_S, invsin_S;
1104 SimdReal st_S, sth_S;
1105 SimdReal cik_S, cii_S, ckk_S;
1106 SimdReal f_ix_S, f_iy_S, f_iz_S;
1107 SimdReal f_kx_S, f_ky_S, f_kz_S;
1108 alignas(GMX_SIMD_ALIGNMENT) real pbc_simd[9 * GMX_SIMD_REAL_WIDTH];
1110 set_pbc_simd(pbc, pbc_simd);
1112 /* nbonds is the number of angles times nfa1, here we step GMX_SIMD_REAL_WIDTH angles */
1113 for (i = 0; (i < nbonds); i += GMX_SIMD_REAL_WIDTH * nfa1)
1115 /* Collect atoms for GMX_SIMD_REAL_WIDTH angles.
1116 * iu indexes into forceatoms, we should not let iu go beyond nbonds.
1119 for (s = 0; s < GMX_SIMD_REAL_WIDTH; s++)
1121 type = forceatoms[iu];
1122 ai[s] = forceatoms[iu + 1];
1123 aj[s] = forceatoms[iu + 2];
1124 ak[s] = forceatoms[iu + 3];
1126 /* At the end fill the arrays with the last atoms and 0 params */
1127 if (i + s * nfa1 < nbonds)
1129 coeff[s] = forceparams[type].harmonic.krA;
1130 coeff[GMX_SIMD_REAL_WIDTH + s] = forceparams[type].harmonic.rA;
1132 if (iu + nfa1 < nbonds)
1140 coeff[GMX_SIMD_REAL_WIDTH + s] = 0;
1144 /* Store the non PBC corrected distances packed and aligned */
1145 gatherLoadUTranspose<3>(reinterpret_cast<const real*>(x), ai, &xi_S, &yi_S, &zi_S);
1146 gatherLoadUTranspose<3>(reinterpret_cast<const real*>(x), aj, &xj_S, &yj_S, &zj_S);
1147 gatherLoadUTranspose<3>(reinterpret_cast<const real*>(x), ak, &xk_S, &yk_S, &zk_S);
1148 rijx_S = xi_S - xj_S;
1149 rijy_S = yi_S - yj_S;
1150 rijz_S = zi_S - zj_S;
1151 rkjx_S = xk_S - xj_S;
1152 rkjy_S = yk_S - yj_S;
1153 rkjz_S = zk_S - zj_S;
1155 k_S = load<SimdReal>(coeff);
1156 theta0_S = load<SimdReal>(coeff + GMX_SIMD_REAL_WIDTH) * deg2rad_S;
1158 pbc_correct_dx_simd(&rijx_S, &rijy_S, &rijz_S, pbc_simd);
1159 pbc_correct_dx_simd(&rkjx_S, &rkjy_S, &rkjz_S, pbc_simd);
1161 rij_rkj_S = iprod(rijx_S, rijy_S, rijz_S, rkjx_S, rkjy_S, rkjz_S);
1163 nrij2_S = norm2(rijx_S, rijy_S, rijz_S);
1164 nrkj2_S = norm2(rkjx_S, rkjy_S, rkjz_S);
1166 nrij_1_S = invsqrt(nrij2_S);
1167 nrkj_1_S = invsqrt(nrkj2_S);
1169 cos_S = rij_rkj_S * nrij_1_S * nrkj_1_S;
1171 /* We compute cos^2 using a division instead of squaring cos_S,
1172 * as we would then get additional error contributions from
1173 * the two invsqrt operations, which get amplified by
1174 * the 1/sqrt(1-cos^2) for cos values close to 1.
1176 * Note that the non-SIMD version of angles() avoids this issue
1177 * by computing the cosine using double precision.
1179 cos2_S = rij_rkj_S * rij_rkj_S / (nrij2_S * nrkj2_S);
1181 /* To allow for 0 and 180 degrees, we need to avoid issues when
1182 * the cosine is close to -1 and 1. For acos() the argument needs
1183 * to be in that range. For cos^2 we take the min of cos^2 and 1 - 1bit,
1184 * so we can safely compute 1/sin as 1/sqrt(1 - cos^2).
1186 cos_S = max(cos_S, -one_S);
1187 cos_S = min(cos_S, one_S);
1188 cos2_S = min(cos2_S, one_min_eps_S);
1190 theta_S = acos(cos_S);
1192 invsin_S = invsqrt(one_S - cos2_S);
1194 st_S = k_S * (theta0_S - theta_S) * invsin_S;
1195 sth_S = st_S * cos_S;
1197 cik_S = st_S * nrij_1_S * nrkj_1_S;
1198 cii_S = sth_S * nrij_1_S * nrij_1_S;
1199 ckk_S = sth_S * nrkj_1_S * nrkj_1_S;
1201 f_ix_S = cii_S * rijx_S;
1202 f_ix_S = fnma(cik_S, rkjx_S, f_ix_S);
1203 f_iy_S = cii_S * rijy_S;
1204 f_iy_S = fnma(cik_S, rkjy_S, f_iy_S);
1205 f_iz_S = cii_S * rijz_S;
1206 f_iz_S = fnma(cik_S, rkjz_S, f_iz_S);
1207 f_kx_S = ckk_S * rkjx_S;
1208 f_kx_S = fnma(cik_S, rijx_S, f_kx_S);
1209 f_ky_S = ckk_S * rkjy_S;
1210 f_ky_S = fnma(cik_S, rijy_S, f_ky_S);
1211 f_kz_S = ckk_S * rkjz_S;
1212 f_kz_S = fnma(cik_S, rijz_S, f_kz_S);
1214 transposeScatterIncrU<4>(reinterpret_cast<real*>(f), ai, f_ix_S, f_iy_S, f_iz_S);
1215 transposeScatterDecrU<4>(reinterpret_cast<real*>(f), aj, f_ix_S + f_kx_S, f_iy_S + f_ky_S,
1217 transposeScatterIncrU<4>(reinterpret_cast<real*>(f), ak, f_kx_S, f_ky_S, f_kz_S);
1223 #endif // GMX_SIMD_HAVE_REAL
1225 template<BondedKernelFlavor flavor>
1226 real linear_angles(int nbonds,
1227 const t_iatom forceatoms[],
1228 const t_iparams forceparams[],
1235 const t_mdatoms gmx_unused* md,
1236 t_fcdata gmx_unused* fcd,
1237 int gmx_unused* global_atom_index)
1239 int i, m, ai, aj, ak, t1, t2, type;
1241 real L1, kA, kB, aA, aB, dr, dr2, va, vtot, a, b, klin;
1242 rvec r_ij, r_kj, r_ik, dx;
1246 for (i = 0; (i < nbonds);)
1248 type = forceatoms[i++];
1249 ai = forceatoms[i++];
1250 aj = forceatoms[i++];
1251 ak = forceatoms[i++];
1253 kA = forceparams[type].linangle.klinA;
1254 kB = forceparams[type].linangle.klinB;
1255 klin = L1 * kA + lambda * kB;
1257 aA = forceparams[type].linangle.aA;
1258 aB = forceparams[type].linangle.aB;
1259 a = L1 * aA + lambda * aB;
1262 t1 = pbc_rvec_sub(pbc, x[ai], x[aj], r_ij);
1263 t2 = pbc_rvec_sub(pbc, x[ak], x[aj], r_kj);
1264 rvec_sub(r_ij, r_kj, r_ik);
1267 for (m = 0; (m < DIM); m++)
1269 dr = -a * r_ij[m] - b * r_kj[m];
1272 f_i[m] = a * klin * dr;
1273 f_k[m] = b * klin * dr;
1274 f_j[m] = -(f_i[m] + f_k[m]);
1279 va = 0.5 * klin * dr2;
1280 *dvdlambda += 0.5 * (kB - kA) * dr2 + klin * (aB - aA) * iprod(dx, r_ik);
1284 if (computeVirial(flavor))
1286 rvec_inc(fshift[t1], f_i);
1287 rvec_inc(fshift[CENTRAL], f_j);
1288 rvec_inc(fshift[t2], f_k);
1294 template<BondedKernelFlavor flavor>
1295 std::enable_if_t<flavor != BondedKernelFlavor::ForcesSimdWhenAvailable || !GMX_SIMD_HAVE_REAL, real>
1296 urey_bradley(int nbonds,
1297 const t_iatom forceatoms[],
1298 const t_iparams forceparams[],
1305 const t_mdatoms gmx_unused* md,
1306 t_fcdata gmx_unused* fcd,
1307 int gmx_unused* global_atom_index)
1309 int i, m, ai, aj, ak, t1, t2, type, ki;
1310 rvec r_ij, r_kj, r_ik;
1311 real cos_theta, cos_theta2, theta;
1312 real dVdt, va, vtot, dr, dr2, vbond, fbond, fik;
1313 real kthA, th0A, kUBA, r13A, kthB, th0B, kUBB, r13B;
1316 for (i = 0; (i < nbonds);)
1318 type = forceatoms[i++];
1319 ai = forceatoms[i++];
1320 aj = forceatoms[i++];
1321 ak = forceatoms[i++];
1322 th0A = forceparams[type].u_b.thetaA * DEG2RAD;
1323 kthA = forceparams[type].u_b.kthetaA;
1324 r13A = forceparams[type].u_b.r13A;
1325 kUBA = forceparams[type].u_b.kUBA;
1326 th0B = forceparams[type].u_b.thetaB * DEG2RAD;
1327 kthB = forceparams[type].u_b.kthetaB;
1328 r13B = forceparams[type].u_b.r13B;
1329 kUBB = forceparams[type].u_b.kUBB;
1331 theta = bond_angle(x[ai], x[aj], x[ak], pbc, r_ij, r_kj, &cos_theta, &t1, &t2); /* 41 */
1333 *dvdlambda += harmonic(kthA, kthB, th0A, th0B, theta, lambda, &va, &dVdt); /* 21 */
1336 ki = pbc_rvec_sub(pbc, x[ai], x[ak], r_ik); /* 3 */
1337 dr2 = iprod(r_ik, r_ik); /* 5 */
1338 dr = dr2 * gmx::invsqrt(dr2); /* 10 */
1340 *dvdlambda += harmonic(kUBA, kUBB, r13A, r13B, dr, lambda, &vbond, &fbond); /* 19 */
1342 cos_theta2 = gmx::square(cos_theta); /* 1 */
1350 st = dVdt * gmx::invsqrt(1 - cos_theta2); /* 12 */
1351 sth = st * cos_theta; /* 1 */
1352 nrkj2 = iprod(r_kj, r_kj); /* 5 */
1353 nrij2 = iprod(r_ij, r_ij);
1355 cik = st * gmx::invsqrt(nrkj2 * nrij2); /* 12 */
1356 cii = sth / nrij2; /* 10 */
1357 ckk = sth / nrkj2; /* 10 */
1359 for (m = 0; (m < DIM); m++) /* 39 */
1361 f_i[m] = -(cik * r_kj[m] - cii * r_ij[m]);
1362 f_k[m] = -(cik * r_ij[m] - ckk * r_kj[m]);
1363 f_j[m] = -f_i[m] - f_k[m];
1368 if (computeVirial(flavor))
1370 rvec_inc(fshift[t1], f_i);
1371 rvec_inc(fshift[CENTRAL], f_j);
1372 rvec_inc(fshift[t2], f_k);
1375 /* Time for the bond calculations */
1381 vtot += vbond; /* 1*/
1382 fbond *= gmx::invsqrt(dr2); /* 6 */
1384 for (m = 0; (m < DIM); m++) /* 15 */
1386 fik = fbond * r_ik[m];
1389 if (computeVirial(flavor))
1391 fshift[ki][m] += fik;
1392 fshift[CENTRAL][m] -= fik;
1399 #if GMX_SIMD_HAVE_REAL
1401 /* As urey_bradley, but using SIMD to calculate many potentials at once.
1402 * This routines does not calculate energies and shift forces.
1404 template<BondedKernelFlavor flavor>
1405 std::enable_if_t<flavor == BondedKernelFlavor::ForcesSimdWhenAvailable, real>
1406 urey_bradley(int nbonds,
1407 const t_iatom forceatoms[],
1408 const t_iparams forceparams[],
1411 rvec gmx_unused fshift[],
1413 real gmx_unused lambda,
1414 real gmx_unused* dvdlambda,
1415 const t_mdatoms gmx_unused* md,
1416 t_fcdata gmx_unused* fcd,
1417 int gmx_unused* global_atom_index)
1419 constexpr int nfa1 = 4;
1420 alignas(GMX_SIMD_ALIGNMENT) std::int32_t ai[GMX_SIMD_REAL_WIDTH];
1421 alignas(GMX_SIMD_ALIGNMENT) std::int32_t aj[GMX_SIMD_REAL_WIDTH];
1422 alignas(GMX_SIMD_ALIGNMENT) std::int32_t ak[GMX_SIMD_REAL_WIDTH];
1423 alignas(GMX_SIMD_ALIGNMENT) real coeff[4 * GMX_SIMD_REAL_WIDTH];
1424 alignas(GMX_SIMD_ALIGNMENT) real pbc_simd[9 * GMX_SIMD_REAL_WIDTH];
1426 set_pbc_simd(pbc, pbc_simd);
1428 /* nbonds is the number of angles times nfa1, here we step GMX_SIMD_REAL_WIDTH angles */
1429 for (int i = 0; i < nbonds; i += GMX_SIMD_REAL_WIDTH * nfa1)
1431 /* Collect atoms for GMX_SIMD_REAL_WIDTH angles.
1432 * iu indexes into forceatoms, we should not let iu go beyond nbonds.
1435 for (int s = 0; s < GMX_SIMD_REAL_WIDTH; s++)
1437 const int type = forceatoms[iu];
1438 ai[s] = forceatoms[iu + 1];
1439 aj[s] = forceatoms[iu + 2];
1440 ak[s] = forceatoms[iu + 3];
1442 /* At the end fill the arrays with the last atoms and 0 params */
1443 if (i + s * nfa1 < nbonds)
1445 coeff[s] = forceparams[type].u_b.kthetaA;
1446 coeff[GMX_SIMD_REAL_WIDTH + s] = forceparams[type].u_b.thetaA;
1447 coeff[GMX_SIMD_REAL_WIDTH * 2 + s] = forceparams[type].u_b.kUBA;
1448 coeff[GMX_SIMD_REAL_WIDTH * 3 + s] = forceparams[type].u_b.r13A;
1450 if (iu + nfa1 < nbonds)
1458 coeff[GMX_SIMD_REAL_WIDTH + s] = 0;
1459 coeff[GMX_SIMD_REAL_WIDTH * 2 + s] = 0;
1460 coeff[GMX_SIMD_REAL_WIDTH * 3 + s] = 0;
1464 SimdReal xi_S, yi_S, zi_S;
1465 SimdReal xj_S, yj_S, zj_S;
1466 SimdReal xk_S, yk_S, zk_S;
1468 /* Store the non PBC corrected distances packed and aligned */
1469 gatherLoadUTranspose<3>(reinterpret_cast<const real*>(x), ai, &xi_S, &yi_S, &zi_S);
1470 gatherLoadUTranspose<3>(reinterpret_cast<const real*>(x), aj, &xj_S, &yj_S, &zj_S);
1471 gatherLoadUTranspose<3>(reinterpret_cast<const real*>(x), ak, &xk_S, &yk_S, &zk_S);
1472 SimdReal rijx_S = xi_S - xj_S;
1473 SimdReal rijy_S = yi_S - yj_S;
1474 SimdReal rijz_S = zi_S - zj_S;
1475 SimdReal rkjx_S = xk_S - xj_S;
1476 SimdReal rkjy_S = yk_S - yj_S;
1477 SimdReal rkjz_S = zk_S - zj_S;
1478 SimdReal rikx_S = xi_S - xk_S;
1479 SimdReal riky_S = yi_S - yk_S;
1480 SimdReal rikz_S = zi_S - zk_S;
1482 const SimdReal ktheta_S = load<SimdReal>(coeff);
1483 const SimdReal theta0_S = load<SimdReal>(coeff + GMX_SIMD_REAL_WIDTH) * DEG2RAD;
1484 const SimdReal kUB_S = load<SimdReal>(coeff + 2 * GMX_SIMD_REAL_WIDTH);
1485 const SimdReal r13_S = load<SimdReal>(coeff + 3 * GMX_SIMD_REAL_WIDTH);
1487 pbc_correct_dx_simd(&rijx_S, &rijy_S, &rijz_S, pbc_simd);
1488 pbc_correct_dx_simd(&rkjx_S, &rkjy_S, &rkjz_S, pbc_simd);
1489 pbc_correct_dx_simd(&rikx_S, &riky_S, &rikz_S, pbc_simd);
1491 const SimdReal rij_rkj_S = iprod(rijx_S, rijy_S, rijz_S, rkjx_S, rkjy_S, rkjz_S);
1493 const SimdReal dr2_S = iprod(rikx_S, riky_S, rikz_S, rikx_S, riky_S, rikz_S);
1495 const SimdReal nrij2_S = norm2(rijx_S, rijy_S, rijz_S);
1496 const SimdReal nrkj2_S = norm2(rkjx_S, rkjy_S, rkjz_S);
1498 const SimdReal nrij_1_S = invsqrt(nrij2_S);
1499 const SimdReal nrkj_1_S = invsqrt(nrkj2_S);
1500 const SimdReal invdr2_S = invsqrt(dr2_S);
1501 const SimdReal dr_S = dr2_S * invdr2_S;
1503 constexpr real min_one_plus_eps = -1.0 + 2.0 * GMX_REAL_EPS; // Smallest number > -1
1505 /* To allow for 180 degrees, we take the max of cos and -1 + 1bit,
1506 * so we can safely get the 1/sin from 1/sqrt(1 - cos^2).
1507 * This also ensures that rounding errors would cause the argument
1508 * of simdAcos to be < -1.
1509 * Note that we do not take precautions for cos(0)=1, so the bonds
1510 * in an angle should not align at an angle of 0 degrees.
1512 const SimdReal cos_S = max(rij_rkj_S * nrij_1_S * nrkj_1_S, min_one_plus_eps);
1514 const SimdReal theta_S = acos(cos_S);
1515 const SimdReal invsin_S = invsqrt(1.0 - cos_S * cos_S);
1516 const SimdReal st_S = ktheta_S * (theta0_S - theta_S) * invsin_S;
1517 const SimdReal sth_S = st_S * cos_S;
1519 const SimdReal cik_S = st_S * nrij_1_S * nrkj_1_S;
1520 const SimdReal cii_S = sth_S * nrij_1_S * nrij_1_S;
1521 const SimdReal ckk_S = sth_S * nrkj_1_S * nrkj_1_S;
1523 const SimdReal sUB_S = kUB_S * (r13_S - dr_S) * invdr2_S;
1525 const SimdReal f_ikx_S = sUB_S * rikx_S;
1526 const SimdReal f_iky_S = sUB_S * riky_S;
1527 const SimdReal f_ikz_S = sUB_S * rikz_S;
1529 const SimdReal f_ix_S = fnma(cik_S, rkjx_S, cii_S * rijx_S) + f_ikx_S;
1530 const SimdReal f_iy_S = fnma(cik_S, rkjy_S, cii_S * rijy_S) + f_iky_S;
1531 const SimdReal f_iz_S = fnma(cik_S, rkjz_S, cii_S * rijz_S) + f_ikz_S;
1532 const SimdReal f_kx_S = fnma(cik_S, rijx_S, ckk_S * rkjx_S) - f_ikx_S;
1533 const SimdReal f_ky_S = fnma(cik_S, rijy_S, ckk_S * rkjy_S) - f_iky_S;
1534 const SimdReal f_kz_S = fnma(cik_S, rijz_S, ckk_S * rkjz_S) - f_ikz_S;
1536 transposeScatterIncrU<4>(reinterpret_cast<real*>(f), ai, f_ix_S, f_iy_S, f_iz_S);
1537 transposeScatterDecrU<4>(reinterpret_cast<real*>(f), aj, f_ix_S + f_kx_S, f_iy_S + f_ky_S,
1539 transposeScatterIncrU<4>(reinterpret_cast<real*>(f), ak, f_kx_S, f_ky_S, f_kz_S);
1545 #endif // GMX_SIMD_HAVE_REAL
1547 template<BondedKernelFlavor flavor>
1548 real quartic_angles(int nbonds,
1549 const t_iatom forceatoms[],
1550 const t_iparams forceparams[],
1555 real gmx_unused lambda,
1556 real gmx_unused* dvdlambda,
1557 const t_mdatoms gmx_unused* md,
1558 t_fcdata gmx_unused* fcd,
1559 int gmx_unused* global_atom_index)
1561 int i, j, ai, aj, ak, t1, t2, type;
1563 real cos_theta, cos_theta2, theta, dt, dVdt, va, dtp, c, vtot;
1566 for (i = 0; (i < nbonds);)
1568 type = forceatoms[i++];
1569 ai = forceatoms[i++];
1570 aj = forceatoms[i++];
1571 ak = forceatoms[i++];
1573 theta = bond_angle(x[ai], x[aj], x[ak], pbc, r_ij, r_kj, &cos_theta, &t1, &t2); /* 41 */
1575 dt = theta - forceparams[type].qangle.theta * DEG2RAD; /* 2 */
1578 va = forceparams[type].qangle.c[0];
1580 for (j = 1; j <= 4; j++)
1582 c = forceparams[type].qangle.c[j];
1583 dVdt -= j * c * dtp;
1591 cos_theta2 = gmx::square(cos_theta); /* 1 */
1600 st = dVdt * gmx::invsqrt(1 - cos_theta2); /* 12 */
1601 sth = st * cos_theta; /* 1 */
1602 nrkj2 = iprod(r_kj, r_kj); /* 5 */
1603 nrij2 = iprod(r_ij, r_ij);
1605 cik = st * gmx::invsqrt(nrkj2 * nrij2); /* 12 */
1606 cii = sth / nrij2; /* 10 */
1607 ckk = sth / nrkj2; /* 10 */
1609 for (m = 0; (m < DIM); m++) /* 39 */
1611 f_i[m] = -(cik * r_kj[m] - cii * r_ij[m]);
1612 f_k[m] = -(cik * r_ij[m] - ckk * r_kj[m]);
1613 f_j[m] = -f_i[m] - f_k[m];
1619 if (computeVirial(flavor))
1621 rvec_inc(fshift[t1], f_i);
1622 rvec_inc(fshift[CENTRAL], f_j);
1623 rvec_inc(fshift[t2], f_k);
1631 #if GMX_SIMD_HAVE_REAL
1633 /* As dih_angle above, but calculates 4 dihedral angles at once using SIMD,
1634 * also calculates the pre-factor required for the dihedral force update.
1635 * Note that bv and buf should be register aligned.
1637 inline void dih_angle_simd(const rvec* x,
1642 const real* pbc_simd,
1650 SimdReal* nrkj_m2_S,
1651 SimdReal* nrkj_n2_S,
1655 SimdReal xi_S, yi_S, zi_S;
1656 SimdReal xj_S, yj_S, zj_S;
1657 SimdReal xk_S, yk_S, zk_S;
1658 SimdReal xl_S, yl_S, zl_S;
1659 SimdReal rijx_S, rijy_S, rijz_S;
1660 SimdReal rkjx_S, rkjy_S, rkjz_S;
1661 SimdReal rklx_S, rkly_S, rklz_S;
1662 SimdReal cx_S, cy_S, cz_S;
1666 SimdReal iprm_S, iprn_S;
1667 SimdReal nrkj2_S, nrkj_1_S, nrkj_2_S, nrkj_S;
1669 SimdReal nrkj2_min_S;
1670 SimdReal real_eps_S;
1672 /* Used to avoid division by zero.
1673 * We take into acount that we multiply the result by real_eps_S.
1675 nrkj2_min_S = SimdReal(GMX_REAL_MIN / (2 * GMX_REAL_EPS));
1677 /* The value of the last significant bit (GMX_REAL_EPS is half of that) */
1678 real_eps_S = SimdReal(2 * GMX_REAL_EPS);
1680 /* Store the non PBC corrected distances packed and aligned */
1681 gatherLoadUTranspose<3>(reinterpret_cast<const real*>(x), ai, &xi_S, &yi_S, &zi_S);
1682 gatherLoadUTranspose<3>(reinterpret_cast<const real*>(x), aj, &xj_S, &yj_S, &zj_S);
1683 gatherLoadUTranspose<3>(reinterpret_cast<const real*>(x), ak, &xk_S, &yk_S, &zk_S);
1684 gatherLoadUTranspose<3>(reinterpret_cast<const real*>(x), al, &xl_S, &yl_S, &zl_S);
1685 rijx_S = xi_S - xj_S;
1686 rijy_S = yi_S - yj_S;
1687 rijz_S = zi_S - zj_S;
1688 rkjx_S = xk_S - xj_S;
1689 rkjy_S = yk_S - yj_S;
1690 rkjz_S = zk_S - zj_S;
1691 rklx_S = xk_S - xl_S;
1692 rkly_S = yk_S - yl_S;
1693 rklz_S = zk_S - zl_S;
1695 pbc_correct_dx_simd(&rijx_S, &rijy_S, &rijz_S, pbc_simd);
1696 pbc_correct_dx_simd(&rkjx_S, &rkjy_S, &rkjz_S, pbc_simd);
1697 pbc_correct_dx_simd(&rklx_S, &rkly_S, &rklz_S, pbc_simd);
1699 cprod(rijx_S, rijy_S, rijz_S, rkjx_S, rkjy_S, rkjz_S, mx_S, my_S, mz_S);
1701 cprod(rkjx_S, rkjy_S, rkjz_S, rklx_S, rkly_S, rklz_S, nx_S, ny_S, nz_S);
1703 cprod(*mx_S, *my_S, *mz_S, *nx_S, *ny_S, *nz_S, &cx_S, &cy_S, &cz_S);
1705 cn_S = sqrt(norm2(cx_S, cy_S, cz_S));
1707 s_S = iprod(*mx_S, *my_S, *mz_S, *nx_S, *ny_S, *nz_S);
1709 /* Determine the dihedral angle, the sign might need correction */
1710 *phi_S = atan2(cn_S, s_S);
1712 ipr_S = iprod(rijx_S, rijy_S, rijz_S, *nx_S, *ny_S, *nz_S);
1714 iprm_S = norm2(*mx_S, *my_S, *mz_S);
1715 iprn_S = norm2(*nx_S, *ny_S, *nz_S);
1717 nrkj2_S = norm2(rkjx_S, rkjy_S, rkjz_S);
1719 /* Avoid division by zero. When zero, the result is multiplied by 0
1720 * anyhow, so the 3 max below do not affect the final result.
1722 nrkj2_S = max(nrkj2_S, nrkj2_min_S);
1723 nrkj_1_S = invsqrt(nrkj2_S);
1724 nrkj_2_S = nrkj_1_S * nrkj_1_S;
1725 nrkj_S = nrkj2_S * nrkj_1_S;
1727 toler_S = nrkj2_S * real_eps_S;
1729 /* Here the plain-C code uses a conditional, but we can't do that in SIMD.
1730 * So we take a max with the tolerance instead. Since we multiply with
1731 * m or n later, the max does not affect the results.
1733 iprm_S = max(iprm_S, toler_S);
1734 iprn_S = max(iprn_S, toler_S);
1735 *nrkj_m2_S = nrkj_S * inv(iprm_S);
1736 *nrkj_n2_S = nrkj_S * inv(iprn_S);
1738 /* Set sign of phi_S with the sign of ipr_S; phi_S is currently positive */
1739 *phi_S = copysign(*phi_S, ipr_S);
1740 *p_S = iprod(rijx_S, rijy_S, rijz_S, rkjx_S, rkjy_S, rkjz_S);
1741 *p_S = *p_S * nrkj_2_S;
1743 *q_S = iprod(rklx_S, rkly_S, rklz_S, rkjx_S, rkjy_S, rkjz_S);
1744 *q_S = *q_S * nrkj_2_S;
1747 #endif // GMX_SIMD_HAVE_REAL
1751 template<BondedKernelFlavor flavor>
1752 void do_dih_fup(int i,
1771 rvec f_i, f_j, f_k, f_l;
1772 rvec uvec, vvec, svec, dx_jl;
1773 real iprm, iprn, nrkj, nrkj2, nrkj_1, nrkj_2;
1774 real a, b, p, q, toler;
1776 iprm = iprod(m, m); /* 5 */
1777 iprn = iprod(n, n); /* 5 */
1778 nrkj2 = iprod(r_kj, r_kj); /* 5 */
1779 toler = nrkj2 * GMX_REAL_EPS;
1780 if ((iprm > toler) && (iprn > toler))
1782 nrkj_1 = gmx::invsqrt(nrkj2); /* 10 */
1783 nrkj_2 = nrkj_1 * nrkj_1; /* 1 */
1784 nrkj = nrkj2 * nrkj_1; /* 1 */
1785 a = -ddphi * nrkj / iprm; /* 11 */
1786 svmul(a, m, f_i); /* 3 */
1787 b = ddphi * nrkj / iprn; /* 11 */
1788 svmul(b, n, f_l); /* 3 */
1789 p = iprod(r_ij, r_kj); /* 5 */
1790 p *= nrkj_2; /* 1 */
1791 q = iprod(r_kl, r_kj); /* 5 */
1792 q *= nrkj_2; /* 1 */
1793 svmul(p, f_i, uvec); /* 3 */
1794 svmul(q, f_l, vvec); /* 3 */
1795 rvec_sub(uvec, vvec, svec); /* 3 */
1796 rvec_sub(f_i, svec, f_j); /* 3 */
1797 rvec_add(f_l, svec, f_k); /* 3 */
1798 rvec_inc(f[i], f_i); /* 3 */
1799 rvec_dec(f[j], f_j); /* 3 */
1800 rvec_dec(f[k], f_k); /* 3 */
1801 rvec_inc(f[l], f_l); /* 3 */
1803 if (computeVirial(flavor))
1807 t3 = pbc_rvec_sub(pbc, x[l], x[j], dx_jl);
1814 rvec_inc(fshift[t1], f_i);
1815 rvec_dec(fshift[CENTRAL], f_j);
1816 rvec_dec(fshift[t2], f_k);
1817 rvec_inc(fshift[t3], f_l);
1826 #if GMX_SIMD_HAVE_REAL
1827 /* As do_dih_fup_noshiftf above, but with SIMD and pre-calculated pre-factors */
1828 inline void gmx_simdcall do_dih_fup_noshiftf_simd(const int* ai,
1842 SimdReal sx = p * f_i_x + q * mf_l_x;
1843 SimdReal sy = p * f_i_y + q * mf_l_y;
1844 SimdReal sz = p * f_i_z + q * mf_l_z;
1845 SimdReal f_j_x = f_i_x - sx;
1846 SimdReal f_j_y = f_i_y - sy;
1847 SimdReal f_j_z = f_i_z - sz;
1848 SimdReal f_k_x = mf_l_x - sx;
1849 SimdReal f_k_y = mf_l_y - sy;
1850 SimdReal f_k_z = mf_l_z - sz;
1851 transposeScatterIncrU<4>(reinterpret_cast<real*>(f), ai, f_i_x, f_i_y, f_i_z);
1852 transposeScatterDecrU<4>(reinterpret_cast<real*>(f), aj, f_j_x, f_j_y, f_j_z);
1853 transposeScatterIncrU<4>(reinterpret_cast<real*>(f), ak, f_k_x, f_k_y, f_k_z);
1854 transposeScatterDecrU<4>(reinterpret_cast<real*>(f), al, mf_l_x, mf_l_y, mf_l_z);
1856 #endif // GMX_SIMD_HAVE_REAL
1858 /*! \brief Computes and returns the proper dihedral force
1860 * With the appropriate kernel flavor, also computes and accumulates
1861 * the energy and dV/dlambda.
1863 template<BondedKernelFlavor flavor>
1864 real dopdihs(real cpA, real cpB, real phiA, real phiB, int mult, real phi, real lambda, real* V, real* dvdlambda)
1866 const real L1 = 1.0 - lambda;
1867 const real ph0 = (L1 * phiA + lambda * phiB) * DEG2RAD;
1868 const real dph0 = (phiB - phiA) * DEG2RAD;
1869 const real cp = L1 * cpA + lambda * cpB;
1871 const real mdphi = mult * phi - ph0;
1872 const real sdphi = std::sin(mdphi);
1873 const real ddphi = -cp * mult * sdphi;
1874 if (computeEnergy(flavor))
1876 const real v1 = 1 + std::cos(mdphi);
1879 *dvdlambda += (cpB - cpA) * v1 + cp * dph0 * sdphi;
1883 /* That was 40 flops */
1886 /*! \brief Similar to \p dopdihs(), except for a minus sign and a different treatment of mult/phi0 */
1887 real dopdihs_min(real cpA, real cpB, real phiA, real phiB, int mult, real phi, real lambda, real* V, real* F)
1889 real v, dvdlambda, mdphi, v1, sdphi, ddphi;
1890 real L1 = 1.0 - lambda;
1891 real ph0 = (L1 * phiA + lambda * phiB) * DEG2RAD;
1892 real dph0 = (phiB - phiA) * DEG2RAD;
1893 real cp = L1 * cpA + lambda * cpB;
1895 mdphi = mult * (phi - ph0);
1896 sdphi = std::sin(mdphi);
1897 ddphi = cp * mult * sdphi;
1898 v1 = 1.0 - std::cos(mdphi);
1901 dvdlambda = (cpB - cpA) * v1 + cp * dph0 * sdphi;
1908 /* That was 40 flops */
1911 template<BondedKernelFlavor flavor>
1912 std::enable_if_t<flavor != BondedKernelFlavor::ForcesSimdWhenAvailable || !GMX_SIMD_HAVE_REAL, real>
1914 const t_iatom forceatoms[],
1915 const t_iparams forceparams[],
1922 const t_mdatoms gmx_unused* md,
1923 t_fcdata gmx_unused* fcd,
1924 int gmx_unused* global_atom_index)
1927 rvec r_ij, r_kj, r_kl, m, n;
1931 for (int i = 0; i < nbonds;)
1933 const int ai = forceatoms[i + 1];
1934 const int aj = forceatoms[i + 2];
1935 const int ak = forceatoms[i + 3];
1936 const int al = forceatoms[i + 4];
1938 const real phi = dih_angle(x[ai], x[aj], x[ak], x[al], pbc, r_ij, r_kj, r_kl, m, n, &t1,
1941 /* Loop over dihedrals working on the same atoms,
1942 * so we avoid recalculating angles and distributing forces.
1947 const int type = forceatoms[i];
1948 ddphi_tot += dopdihs<flavor>(forceparams[type].pdihs.cpA, forceparams[type].pdihs.cpB,
1949 forceparams[type].pdihs.phiA, forceparams[type].pdihs.phiB,
1950 forceparams[type].pdihs.mult, phi, lambda, &vtot, dvdlambda);
1953 } while (i < nbonds && forceatoms[i + 1] == ai && forceatoms[i + 2] == aj
1954 && forceatoms[i + 3] == ak && forceatoms[i + 4] == al);
1956 do_dih_fup<flavor>(ai, aj, ak, al, ddphi_tot, r_ij, r_kj, r_kl, m, n, f, fshift, pbc, x, t1,
1963 #if GMX_SIMD_HAVE_REAL
1965 /* As pdihs above, but using SIMD to calculate multiple dihedrals at once */
1966 template<BondedKernelFlavor flavor>
1967 std::enable_if_t<flavor == BondedKernelFlavor::ForcesSimdWhenAvailable, real>
1969 const t_iatom forceatoms[],
1970 const t_iparams forceparams[],
1973 rvec gmx_unused fshift[],
1975 real gmx_unused lambda,
1976 real gmx_unused* dvdlambda,
1977 const t_mdatoms gmx_unused* md,
1978 t_fcdata gmx_unused* fcd,
1979 int gmx_unused* global_atom_index)
1984 alignas(GMX_SIMD_ALIGNMENT) std::int32_t ai[GMX_SIMD_REAL_WIDTH];
1985 alignas(GMX_SIMD_ALIGNMENT) std::int32_t aj[GMX_SIMD_REAL_WIDTH];
1986 alignas(GMX_SIMD_ALIGNMENT) std::int32_t ak[GMX_SIMD_REAL_WIDTH];
1987 alignas(GMX_SIMD_ALIGNMENT) std::int32_t al[GMX_SIMD_REAL_WIDTH];
1988 alignas(GMX_SIMD_ALIGNMENT) real buf[3 * GMX_SIMD_REAL_WIDTH];
1989 real * cp, *phi0, *mult;
1990 SimdReal deg2rad_S(DEG2RAD);
1992 SimdReal phi0_S, phi_S;
1993 SimdReal mx_S, my_S, mz_S;
1994 SimdReal nx_S, ny_S, nz_S;
1995 SimdReal nrkj_m2_S, nrkj_n2_S;
1996 SimdReal cp_S, mdphi_S, mult_S;
1997 SimdReal sin_S, cos_S;
1999 SimdReal sf_i_S, msf_l_S;
2000 alignas(GMX_SIMD_ALIGNMENT) real pbc_simd[9 * GMX_SIMD_REAL_WIDTH];
2002 /* Extract aligned pointer for parameters and variables */
2003 cp = buf + 0 * GMX_SIMD_REAL_WIDTH;
2004 phi0 = buf + 1 * GMX_SIMD_REAL_WIDTH;
2005 mult = buf + 2 * GMX_SIMD_REAL_WIDTH;
2007 set_pbc_simd(pbc, pbc_simd);
2009 /* nbonds is the number of dihedrals times nfa1, here we step GMX_SIMD_REAL_WIDTH dihs */
2010 for (i = 0; (i < nbonds); i += GMX_SIMD_REAL_WIDTH * nfa1)
2012 /* Collect atoms quadruplets for GMX_SIMD_REAL_WIDTH dihedrals.
2013 * iu indexes into forceatoms, we should not let iu go beyond nbonds.
2016 for (s = 0; s < GMX_SIMD_REAL_WIDTH; s++)
2018 type = forceatoms[iu];
2019 ai[s] = forceatoms[iu + 1];
2020 aj[s] = forceatoms[iu + 2];
2021 ak[s] = forceatoms[iu + 3];
2022 al[s] = forceatoms[iu + 4];
2024 /* At the end fill the arrays with the last atoms and 0 params */
2025 if (i + s * nfa1 < nbonds)
2027 cp[s] = forceparams[type].pdihs.cpA;
2028 phi0[s] = forceparams[type].pdihs.phiA;
2029 mult[s] = forceparams[type].pdihs.mult;
2031 if (iu + nfa1 < nbonds)
2044 /* Calculate GMX_SIMD_REAL_WIDTH dihedral angles at once */
2045 dih_angle_simd(x, ai, aj, ak, al, pbc_simd, &phi_S, &mx_S, &my_S, &mz_S, &nx_S, &ny_S,
2046 &nz_S, &nrkj_m2_S, &nrkj_n2_S, &p_S, &q_S);
2048 cp_S = load<SimdReal>(cp);
2049 phi0_S = load<SimdReal>(phi0) * deg2rad_S;
2050 mult_S = load<SimdReal>(mult);
2052 mdphi_S = fms(mult_S, phi_S, phi0_S);
2054 /* Calculate GMX_SIMD_REAL_WIDTH sines at once */
2055 sincos(mdphi_S, &sin_S, &cos_S);
2056 mddphi_S = cp_S * mult_S * sin_S;
2057 sf_i_S = mddphi_S * nrkj_m2_S;
2058 msf_l_S = mddphi_S * nrkj_n2_S;
2060 /* After this m?_S will contain f[i] */
2061 mx_S = sf_i_S * mx_S;
2062 my_S = sf_i_S * my_S;
2063 mz_S = sf_i_S * mz_S;
2065 /* After this m?_S will contain -f[l] */
2066 nx_S = msf_l_S * nx_S;
2067 ny_S = msf_l_S * ny_S;
2068 nz_S = msf_l_S * nz_S;
2070 do_dih_fup_noshiftf_simd(ai, aj, ak, al, p_S, q_S, mx_S, my_S, mz_S, nx_S, ny_S, nz_S, f);
2076 /* This is mostly a copy of the SIMD flavor of pdihs above, but with using
2077 * the RB potential instead of a harmonic potential.
2078 * This function can replace rbdihs() when no energy and virial are needed.
2080 template<BondedKernelFlavor flavor>
2081 std::enable_if_t<flavor == BondedKernelFlavor::ForcesSimdWhenAvailable, real>
2083 const t_iatom forceatoms[],
2084 const t_iparams forceparams[],
2087 rvec gmx_unused fshift[],
2089 real gmx_unused lambda,
2090 real gmx_unused* dvdlambda,
2091 const t_mdatoms gmx_unused* md,
2092 t_fcdata gmx_unused* fcd,
2093 int gmx_unused* global_atom_index)
2098 alignas(GMX_SIMD_ALIGNMENT) std::int32_t ai[GMX_SIMD_REAL_WIDTH];
2099 alignas(GMX_SIMD_ALIGNMENT) std::int32_t aj[GMX_SIMD_REAL_WIDTH];
2100 alignas(GMX_SIMD_ALIGNMENT) std::int32_t ak[GMX_SIMD_REAL_WIDTH];
2101 alignas(GMX_SIMD_ALIGNMENT) std::int32_t al[GMX_SIMD_REAL_WIDTH];
2102 alignas(GMX_SIMD_ALIGNMENT) real parm[NR_RBDIHS * GMX_SIMD_REAL_WIDTH];
2106 SimdReal ddphi_S, cosfac_S;
2107 SimdReal mx_S, my_S, mz_S;
2108 SimdReal nx_S, ny_S, nz_S;
2109 SimdReal nrkj_m2_S, nrkj_n2_S;
2110 SimdReal parm_S, c_S;
2111 SimdReal sin_S, cos_S;
2112 SimdReal sf_i_S, msf_l_S;
2113 alignas(GMX_SIMD_ALIGNMENT) real pbc_simd[9 * GMX_SIMD_REAL_WIDTH];
2115 SimdReal pi_S(M_PI);
2116 SimdReal one_S(1.0);
2118 set_pbc_simd(pbc, pbc_simd);
2120 /* nbonds is the number of dihedrals times nfa1, here we step GMX_SIMD_REAL_WIDTH dihs */
2121 for (i = 0; (i < nbonds); i += GMX_SIMD_REAL_WIDTH * nfa1)
2123 /* Collect atoms quadruplets for GMX_SIMD_REAL_WIDTH dihedrals.
2124 * iu indexes into forceatoms, we should not let iu go beyond nbonds.
2127 for (s = 0; s < GMX_SIMD_REAL_WIDTH; s++)
2129 type = forceatoms[iu];
2130 ai[s] = forceatoms[iu + 1];
2131 aj[s] = forceatoms[iu + 2];
2132 ak[s] = forceatoms[iu + 3];
2133 al[s] = forceatoms[iu + 4];
2135 /* At the end fill the arrays with the last atoms and 0 params */
2136 if (i + s * nfa1 < nbonds)
2138 /* We don't need the first parameter, since that's a constant
2139 * which only affects the energies, not the forces.
2141 for (j = 1; j < NR_RBDIHS; j++)
2143 parm[j * GMX_SIMD_REAL_WIDTH + s] = forceparams[type].rbdihs.rbcA[j];
2146 if (iu + nfa1 < nbonds)
2153 for (j = 1; j < NR_RBDIHS; j++)
2155 parm[j * GMX_SIMD_REAL_WIDTH + s] = 0;
2160 /* Calculate GMX_SIMD_REAL_WIDTH dihedral angles at once */
2161 dih_angle_simd(x, ai, aj, ak, al, pbc_simd, &phi_S, &mx_S, &my_S, &mz_S, &nx_S, &ny_S,
2162 &nz_S, &nrkj_m2_S, &nrkj_n2_S, &p_S, &q_S);
2164 /* Change to polymer convention */
2165 phi_S = phi_S - pi_S;
2167 sincos(phi_S, &sin_S, &cos_S);
2169 ddphi_S = setZero();
2172 for (j = 1; j < NR_RBDIHS; j++)
2174 parm_S = load<SimdReal>(parm + j * GMX_SIMD_REAL_WIDTH);
2175 ddphi_S = fma(c_S * parm_S, cosfac_S, ddphi_S);
2176 cosfac_S = cosfac_S * cos_S;
2180 /* Note that here we do not use the minus sign which is present
2181 * in the normal RB code. This is corrected for through (m)sf below.
2183 ddphi_S = ddphi_S * sin_S;
2185 sf_i_S = ddphi_S * nrkj_m2_S;
2186 msf_l_S = ddphi_S * nrkj_n2_S;
2188 /* After this m?_S will contain f[i] */
2189 mx_S = sf_i_S * mx_S;
2190 my_S = sf_i_S * my_S;
2191 mz_S = sf_i_S * mz_S;
2193 /* After this m?_S will contain -f[l] */
2194 nx_S = msf_l_S * nx_S;
2195 ny_S = msf_l_S * ny_S;
2196 nz_S = msf_l_S * nz_S;
2198 do_dih_fup_noshiftf_simd(ai, aj, ak, al, p_S, q_S, mx_S, my_S, mz_S, nx_S, ny_S, nz_S, f);
2204 #endif // GMX_SIMD_HAVE_REAL
2207 template<BondedKernelFlavor flavor>
2208 real idihs(int nbonds,
2209 const t_iatom forceatoms[],
2210 const t_iparams forceparams[],
2217 const t_mdatoms gmx_unused* md,
2218 t_fcdata gmx_unused* fcd,
2219 int gmx_unused* global_atom_index)
2221 int i, type, ai, aj, ak, al;
2223 real phi, phi0, dphi0, ddphi, vtot;
2224 rvec r_ij, r_kj, r_kl, m, n;
2225 real L1, kk, dp, dp2, kA, kB, pA, pB, dvdl_term;
2230 for (i = 0; (i < nbonds);)
2232 type = forceatoms[i++];
2233 ai = forceatoms[i++];
2234 aj = forceatoms[i++];
2235 ak = forceatoms[i++];
2236 al = forceatoms[i++];
2238 phi = dih_angle(x[ai], x[aj], x[ak], x[al], pbc, r_ij, r_kj, r_kl, m, n, &t1, &t2, &t3); /* 84 */
2240 /* phi can jump if phi0 is close to Pi/-Pi, which will cause huge
2241 * force changes if we just apply a normal harmonic.
2242 * Instead, we first calculate phi-phi0 and take it modulo (-Pi,Pi).
2243 * This means we will never have the periodicity problem, unless
2244 * the dihedral is Pi away from phiO, which is very unlikely due to
2247 kA = forceparams[type].harmonic.krA;
2248 kB = forceparams[type].harmonic.krB;
2249 pA = forceparams[type].harmonic.rA;
2250 pB = forceparams[type].harmonic.rB;
2252 kk = L1 * kA + lambda * kB;
2253 phi0 = (L1 * pA + lambda * pB) * DEG2RAD;
2254 dphi0 = (pB - pA) * DEG2RAD;
2258 make_dp_periodic(&dp);
2262 vtot += 0.5 * kk * dp2;
2265 dvdl_term += 0.5 * (kB - kA) * dp2 - kk * dphi0 * dp;
2267 do_dih_fup<flavor>(ai, aj, ak, al, -ddphi, r_ij, r_kj, r_kl, m, n, f, fshift, pbc, x, t1,
2272 *dvdlambda += dvdl_term;
2276 /*! \brief Computes angle restraints of two different types */
2277 template<BondedKernelFlavor flavor>
2278 real low_angres(int nbonds,
2279 const t_iatom forceatoms[],
2280 const t_iparams forceparams[],
2289 int i, m, type, ai, aj, ak, al;
2291 real phi, cos_phi, cos_phi2, vid, vtot, dVdphi;
2292 rvec r_ij, r_kl, f_i, f_k = { 0, 0, 0 };
2293 real st, sth, nrij2, nrkl2, c, cij, ckl;
2295 t2 = 0; /* avoid warning with gcc-3.3. It is never used uninitialized */
2298 ak = al = 0; /* to avoid warnings */
2299 for (i = 0; i < nbonds;)
2301 type = forceatoms[i++];
2302 ai = forceatoms[i++];
2303 aj = forceatoms[i++];
2304 t1 = pbc_rvec_sub(pbc, x[aj], x[ai], r_ij); /* 3 */
2307 ak = forceatoms[i++];
2308 al = forceatoms[i++];
2309 t2 = pbc_rvec_sub(pbc, x[al], x[ak], r_kl); /* 3 */
2318 cos_phi = cos_angle(r_ij, r_kl); /* 25 */
2319 phi = std::acos(cos_phi); /* 10 */
2321 *dvdlambda += dopdihs_min(forceparams[type].pdihs.cpA, forceparams[type].pdihs.cpB,
2322 forceparams[type].pdihs.phiA, forceparams[type].pdihs.phiB,
2323 forceparams[type].pdihs.mult, phi, lambda, &vid, &dVdphi); /* 40 */
2327 cos_phi2 = gmx::square(cos_phi); /* 1 */
2330 st = -dVdphi * gmx::invsqrt(1 - cos_phi2); /* 12 */
2331 sth = st * cos_phi; /* 1 */
2332 nrij2 = iprod(r_ij, r_ij); /* 5 */
2333 nrkl2 = iprod(r_kl, r_kl); /* 5 */
2335 c = st * gmx::invsqrt(nrij2 * nrkl2); /* 11 */
2336 cij = sth / nrij2; /* 10 */
2337 ckl = sth / nrkl2; /* 10 */
2339 for (m = 0; m < DIM; m++) /* 18+18 */
2341 f_i[m] = (c * r_kl[m] - cij * r_ij[m]);
2346 f_k[m] = (c * r_ij[m] - ckl * r_kl[m]);
2352 if (computeVirial(flavor))
2354 rvec_inc(fshift[t1], f_i);
2355 rvec_dec(fshift[CENTRAL], f_i);
2358 rvec_inc(fshift[t2], f_k);
2359 rvec_dec(fshift[CENTRAL], f_k);
2365 return vtot; /* 184 / 157 (bZAxis) total */
2368 template<BondedKernelFlavor flavor>
2369 real angres(int nbonds,
2370 const t_iatom forceatoms[],
2371 const t_iparams forceparams[],
2378 const t_mdatoms gmx_unused* md,
2379 t_fcdata gmx_unused* fcd,
2380 int gmx_unused* global_atom_index)
2382 return low_angres<flavor>(nbonds, forceatoms, forceparams, x, f, fshift, pbc, lambda, dvdlambda, FALSE);
2385 template<BondedKernelFlavor flavor>
2386 real angresz(int nbonds,
2387 const t_iatom forceatoms[],
2388 const t_iparams forceparams[],
2395 const t_mdatoms gmx_unused* md,
2396 t_fcdata gmx_unused* fcd,
2397 int gmx_unused* global_atom_index)
2399 return low_angres<flavor>(nbonds, forceatoms, forceparams, x, f, fshift, pbc, lambda, dvdlambda, TRUE);
2402 template<BondedKernelFlavor flavor>
2403 real dihres(int nbonds,
2404 const t_iatom forceatoms[],
2405 const t_iparams forceparams[],
2412 const t_mdatoms gmx_unused* md,
2413 t_fcdata gmx_unused* fcd,
2414 int gmx_unused* global_atom_index)
2417 int ai, aj, ak, al, i, type, t1, t2, t3;
2418 real phi0A, phi0B, dphiA, dphiB, kfacA, kfacB, phi0, dphi, kfac;
2419 real phi, ddphi, ddp, ddp2, dp, d2r, L1;
2420 rvec r_ij, r_kj, r_kl, m, n;
2426 for (i = 0; (i < nbonds);)
2428 type = forceatoms[i++];
2429 ai = forceatoms[i++];
2430 aj = forceatoms[i++];
2431 ak = forceatoms[i++];
2432 al = forceatoms[i++];
2434 phi0A = forceparams[type].dihres.phiA * d2r;
2435 dphiA = forceparams[type].dihres.dphiA * d2r;
2436 kfacA = forceparams[type].dihres.kfacA;
2438 phi0B = forceparams[type].dihres.phiB * d2r;
2439 dphiB = forceparams[type].dihres.dphiB * d2r;
2440 kfacB = forceparams[type].dihres.kfacB;
2442 phi0 = L1 * phi0A + lambda * phi0B;
2443 dphi = L1 * dphiA + lambda * dphiB;
2444 kfac = L1 * kfacA + lambda * kfacB;
2446 phi = dih_angle(x[ai], x[aj], x[ak], x[al], pbc, r_ij, r_kj, r_kl, m, n, &t1, &t2, &t3);
2449 /* phi can jump if phi0 is close to Pi/-Pi, which will cause huge
2450 * force changes if we just apply a normal harmonic.
2451 * Instead, we first calculate phi-phi0 and take it modulo (-Pi,Pi).
2452 * This means we will never have the periodicity problem, unless
2453 * the dihedral is Pi away from phiO, which is very unlikely due to
2457 make_dp_periodic(&dp);
2463 else if (dp < -dphi)
2475 vtot += 0.5 * kfac * ddp2;
2478 *dvdlambda += 0.5 * (kfacB - kfacA) * ddp2;
2479 /* lambda dependence from changing restraint distances */
2482 *dvdlambda -= kfac * ddp * ((dphiB - dphiA) + (phi0B - phi0A));
2486 *dvdlambda += kfac * ddp * ((dphiB - dphiA) - (phi0B - phi0A));
2488 do_dih_fup<flavor>(ai, aj, ak, al, ddphi, r_ij, r_kj, r_kl, m, n, f, fshift, pbc, x, t1,
2496 real unimplemented(int gmx_unused nbonds,
2497 const t_iatom gmx_unused forceatoms[],
2498 const t_iparams gmx_unused forceparams[],
2499 const rvec gmx_unused x[],
2500 rvec4 gmx_unused f[],
2501 rvec gmx_unused fshift[],
2502 const t_pbc gmx_unused* pbc,
2503 real gmx_unused lambda,
2504 real gmx_unused* dvdlambda,
2505 const t_mdatoms gmx_unused* md,
2506 t_fcdata gmx_unused* fcd,
2507 int gmx_unused* global_atom_index)
2509 gmx_impl("*** you are using a not implemented function");
2512 template<BondedKernelFlavor flavor>
2513 real restrangles(int nbonds,
2514 const t_iatom forceatoms[],
2515 const t_iparams forceparams[],
2520 real gmx_unused lambda,
2521 real gmx_unused* dvdlambda,
2522 const t_mdatoms gmx_unused* md,
2523 t_fcdata gmx_unused* fcd,
2524 int gmx_unused* global_atom_index)
2526 int i, d, ai, aj, ak, type, m;
2530 double prefactor, ratio_ante, ratio_post;
2531 rvec delta_ante, delta_post, vec_temp;
2534 for (i = 0; (i < nbonds);)
2536 type = forceatoms[i++];
2537 ai = forceatoms[i++];
2538 aj = forceatoms[i++];
2539 ak = forceatoms[i++];
2541 t1 = pbc_rvec_sub(pbc, x[ai], x[aj], vec_temp);
2542 pbc_rvec_sub(pbc, x[aj], x[ai], delta_ante);
2543 t2 = pbc_rvec_sub(pbc, x[ak], x[aj], delta_post);
2546 /* This function computes factors needed for restricted angle potential.
2547 * The restricted angle potential is used in coarse-grained simulations to avoid singularities
2548 * when three particles align and the dihedral angle and dihedral potential
2549 * cannot be calculated. This potential is calculated using the formula:
2550 real restrangles(int nbonds,
2551 const t_iatom forceatoms[],const t_iparams forceparams[],
2552 const rvec x[],rvec4 f[],rvec fshift[],
2554 real gmx_unused lambda,real gmx_unused *dvdlambda,
2555 const t_mdatoms gmx_unused *md,t_fcdata gmx_unused *fcd,
2556 int gmx_unused *global_atom_index)
2558 int i, d, ai, aj, ak, type, m;
2563 real prefactor, ratio_ante, ratio_post;
2564 rvec delta_ante, delta_post, vec_temp;
2567 for(i=0; (i<nbonds); )
2569 type = forceatoms[i++];
2570 ai = forceatoms[i++];
2571 aj = forceatoms[i++];
2572 ak = forceatoms[i++];
2574 * \f[V_{\rm ReB}(\theta_i) = \frac{1}{2} k_{\theta} \frac{(\cos\theta_i - \cos\theta_0)^2}
2575 * {\sin^2\theta_i}\f] ({eq:ReB} and ref \cite{MonicaGoga2013} from the manual).
2576 * For more explanations see comments file "restcbt.h". */
2578 compute_factors_restangles(type, forceparams, delta_ante, delta_post, &prefactor,
2579 &ratio_ante, &ratio_post, &v);
2581 /* Forces are computed per component */
2582 for (d = 0; d < DIM; d++)
2584 f_i[d] = prefactor * (ratio_ante * delta_ante[d] - delta_post[d]);
2586 * ((ratio_post + 1.0) * delta_post[d] - (ratio_ante + 1.0) * delta_ante[d]);
2587 f_k[d] = prefactor * (delta_ante[d] - ratio_post * delta_post[d]);
2590 /* Computation of potential energy */
2596 for (m = 0; (m < DIM); m++)
2603 if (computeVirial(flavor))
2605 rvec_inc(fshift[t1], f_i);
2606 rvec_inc(fshift[CENTRAL], f_j);
2607 rvec_inc(fshift[t2], f_k);
2614 template<BondedKernelFlavor flavor>
2615 real restrdihs(int nbonds,
2616 const t_iatom forceatoms[],
2617 const t_iparams forceparams[],
2622 real gmx_unused lambda,
2623 real gmx_unused* dvlambda,
2624 const t_mdatoms gmx_unused* md,
2625 t_fcdata gmx_unused* fcd,
2626 int gmx_unused* global_atom_index)
2628 int i, d, type, ai, aj, ak, al;
2629 rvec f_i, f_j, f_k, f_l;
2633 rvec delta_ante, delta_crnt, delta_post, vec_temp;
2634 real factor_phi_ai_ante, factor_phi_ai_crnt, factor_phi_ai_post;
2635 real factor_phi_aj_ante, factor_phi_aj_crnt, factor_phi_aj_post;
2636 real factor_phi_ak_ante, factor_phi_ak_crnt, factor_phi_ak_post;
2637 real factor_phi_al_ante, factor_phi_al_crnt, factor_phi_al_post;
2642 for (i = 0; (i < nbonds);)
2644 type = forceatoms[i++];
2645 ai = forceatoms[i++];
2646 aj = forceatoms[i++];
2647 ak = forceatoms[i++];
2648 al = forceatoms[i++];
2650 t1 = pbc_rvec_sub(pbc, x[ai], x[aj], vec_temp);
2651 pbc_rvec_sub(pbc, x[aj], x[ai], delta_ante);
2652 t2 = pbc_rvec_sub(pbc, x[ak], x[aj], delta_crnt);
2653 pbc_rvec_sub(pbc, x[ak], x[al], vec_temp);
2654 pbc_rvec_sub(pbc, x[al], x[ak], delta_post);
2656 /* This function computes factors needed for restricted angle potential.
2657 * The restricted angle potential is used in coarse-grained simulations to avoid singularities
2658 * when three particles align and the dihedral angle and dihedral potential cannot be calculated.
2659 * This potential is calculated using the formula:
2660 * \f[V_{\rm ReB}(\theta_i) = \frac{1}{2} k_{\theta}
2661 * \frac{(\cos\theta_i - \cos\theta_0)^2}{\sin^2\theta_i}\f]
2662 * ({eq:ReB} and ref \cite{MonicaGoga2013} from the manual).
2663 * For more explanations see comments file "restcbt.h" */
2665 compute_factors_restrdihs(
2666 type, forceparams, delta_ante, delta_crnt, delta_post, &factor_phi_ai_ante,
2667 &factor_phi_ai_crnt, &factor_phi_ai_post, &factor_phi_aj_ante, &factor_phi_aj_crnt,
2668 &factor_phi_aj_post, &factor_phi_ak_ante, &factor_phi_ak_crnt, &factor_phi_ak_post,
2669 &factor_phi_al_ante, &factor_phi_al_crnt, &factor_phi_al_post, &prefactor_phi, &v);
2672 /* Computation of forces per component */
2673 for (d = 0; d < DIM; d++)
2675 f_i[d] = prefactor_phi
2676 * (factor_phi_ai_ante * delta_ante[d] + factor_phi_ai_crnt * delta_crnt[d]
2677 + factor_phi_ai_post * delta_post[d]);
2678 f_j[d] = prefactor_phi
2679 * (factor_phi_aj_ante * delta_ante[d] + factor_phi_aj_crnt * delta_crnt[d]
2680 + factor_phi_aj_post * delta_post[d]);
2681 f_k[d] = prefactor_phi
2682 * (factor_phi_ak_ante * delta_ante[d] + factor_phi_ak_crnt * delta_crnt[d]
2683 + factor_phi_ak_post * delta_post[d]);
2684 f_l[d] = prefactor_phi
2685 * (factor_phi_al_ante * delta_ante[d] + factor_phi_al_crnt * delta_crnt[d]
2686 + factor_phi_al_post * delta_post[d]);
2688 /* Computation of the energy */
2693 /* Updating the forces */
2695 rvec_inc(f[ai], f_i);
2696 rvec_inc(f[aj], f_j);
2697 rvec_inc(f[ak], f_k);
2698 rvec_inc(f[al], f_l);
2700 if (computeVirial(flavor))
2704 t3 = pbc_rvec_sub(pbc, x[al], x[aj], dx_jl);
2711 rvec_inc(fshift[t1], f_i);
2712 rvec_inc(fshift[CENTRAL], f_j);
2713 rvec_inc(fshift[t2], f_k);
2714 rvec_inc(fshift[t3], f_l);
2722 template<BondedKernelFlavor flavor>
2723 real cbtdihs(int nbonds,
2724 const t_iatom forceatoms[],
2725 const t_iparams forceparams[],
2730 real gmx_unused lambda,
2731 real gmx_unused* dvdlambda,
2732 const t_mdatoms gmx_unused* md,
2733 t_fcdata gmx_unused* fcd,
2734 int gmx_unused* global_atom_index)
2736 int type, ai, aj, ak, al, i, d;
2740 rvec f_i, f_j, f_k, f_l;
2742 rvec delta_ante, delta_crnt, delta_post;
2743 rvec f_phi_ai, f_phi_aj, f_phi_ak, f_phi_al;
2744 rvec f_theta_ante_ai, f_theta_ante_aj, f_theta_ante_ak;
2745 rvec f_theta_post_aj, f_theta_post_ak, f_theta_post_al;
2749 for (i = 0; (i < nbonds);)
2751 type = forceatoms[i++];
2752 ai = forceatoms[i++];
2753 aj = forceatoms[i++];
2754 ak = forceatoms[i++];
2755 al = forceatoms[i++];
2758 t1 = pbc_rvec_sub(pbc, x[ai], x[aj], vec_temp);
2759 pbc_rvec_sub(pbc, x[aj], x[ai], delta_ante);
2760 t2 = pbc_rvec_sub(pbc, x[ak], x[aj], vec_temp);
2761 pbc_rvec_sub(pbc, x[ak], x[aj], delta_crnt);
2762 pbc_rvec_sub(pbc, x[ak], x[al], vec_temp);
2763 pbc_rvec_sub(pbc, x[al], x[ak], delta_post);
2765 /* \brief Compute factors for CBT potential
2766 * The combined bending-torsion potential goes to zero in a very smooth manner, eliminating the numerical
2767 * instabilities, when three coarse-grained particles align and the dihedral angle and standard
2768 * dihedral potentials cannot be calculated. The CBT potential is calculated using the formula:
2769 * \f[V_{\rm CBT}(\theta_{i-1}, \theta_i, \phi_i) = k_{\phi} \sin^3\theta_{i-1} \sin^3\theta_{i}
2770 * \sum_{n=0}^4 { a_n \cos^n\phi_i}\f] ({eq:CBT} and ref \cite{MonicaGoga2013} from the manual).
2771 * It contains in its expression not only the dihedral angle \f$\phi\f$
2772 * but also \f[\theta_{i-1}\f] (theta_ante bellow) and \f[\theta_{i}\f] (theta_post bellow)
2773 * --- the adjacent bending angles.
2774 * For more explanations see comments file "restcbt.h". */
2776 compute_factors_cbtdihs(type, forceparams, delta_ante, delta_crnt, delta_post, f_phi_ai,
2777 f_phi_aj, f_phi_ak, f_phi_al, f_theta_ante_ai, f_theta_ante_aj,
2778 f_theta_ante_ak, f_theta_post_aj, f_theta_post_ak, f_theta_post_al, &v);
2781 /* Acumulate the resuts per beads */
2782 for (d = 0; d < DIM; d++)
2784 f_i[d] = f_phi_ai[d] + f_theta_ante_ai[d];
2785 f_j[d] = f_phi_aj[d] + f_theta_ante_aj[d] + f_theta_post_aj[d];
2786 f_k[d] = f_phi_ak[d] + f_theta_ante_ak[d] + f_theta_post_ak[d];
2787 f_l[d] = f_phi_al[d] + f_theta_post_al[d];
2790 /* Compute the potential energy */
2795 /* Updating the forces */
2796 rvec_inc(f[ai], f_i);
2797 rvec_inc(f[aj], f_j);
2798 rvec_inc(f[ak], f_k);
2799 rvec_inc(f[al], f_l);
2802 if (computeVirial(flavor))
2806 t3 = pbc_rvec_sub(pbc, x[al], x[aj], dx_jl);
2813 rvec_inc(fshift[t1], f_i);
2814 rvec_inc(fshift[CENTRAL], f_j);
2815 rvec_inc(fshift[t2], f_k);
2816 rvec_inc(fshift[t3], f_l);
2823 template<BondedKernelFlavor flavor>
2824 std::enable_if_t<flavor != BondedKernelFlavor::ForcesSimdWhenAvailable || !GMX_SIMD_HAVE_REAL, real>
2826 const t_iatom forceatoms[],
2827 const t_iparams forceparams[],
2834 const t_mdatoms gmx_unused* md,
2835 t_fcdata gmx_unused* fcd,
2836 int gmx_unused* global_atom_index)
2838 const real c0 = 0.0, c1 = 1.0, c2 = 2.0, c3 = 3.0, c4 = 4.0, c5 = 5.0;
2839 int type, ai, aj, ak, al, i, j;
2841 rvec r_ij, r_kj, r_kl, m, n;
2842 real parmA[NR_RBDIHS];
2843 real parmB[NR_RBDIHS];
2844 real parm[NR_RBDIHS];
2845 real cos_phi, phi, rbp, rbpBA;
2846 real v, ddphi, sin_phi;
2848 real L1 = 1.0 - lambda;
2852 for (i = 0; (i < nbonds);)
2854 type = forceatoms[i++];
2855 ai = forceatoms[i++];
2856 aj = forceatoms[i++];
2857 ak = forceatoms[i++];
2858 al = forceatoms[i++];
2860 phi = dih_angle(x[ai], x[aj], x[ak], x[al], pbc, r_ij, r_kj, r_kl, m, n, &t1, &t2, &t3); /* 84 */
2862 /* Change to polymer convention */
2869 phi -= M_PI; /* 1 */
2871 cos_phi = std::cos(phi);
2872 /* Beware of accuracy loss, cannot use 1-sqrt(cos^2) ! */
2873 sin_phi = std::sin(phi);
2875 for (j = 0; (j < NR_RBDIHS); j++)
2877 parmA[j] = forceparams[type].rbdihs.rbcA[j];
2878 parmB[j] = forceparams[type].rbdihs.rbcB[j];
2879 parm[j] = L1 * parmA[j] + lambda * parmB[j];
2881 /* Calculate cosine powers */
2882 /* Calculate the energy */
2883 /* Calculate the derivative */
2886 dvdl_term += (parmB[0] - parmA[0]);
2891 rbpBA = parmB[1] - parmA[1];
2892 ddphi += rbp * cosfac;
2895 dvdl_term += cosfac * rbpBA;
2897 rbpBA = parmB[2] - parmA[2];
2898 ddphi += c2 * rbp * cosfac;
2901 dvdl_term += cosfac * rbpBA;
2903 rbpBA = parmB[3] - parmA[3];
2904 ddphi += c3 * rbp * cosfac;
2907 dvdl_term += cosfac * rbpBA;
2909 rbpBA = parmB[4] - parmA[4];
2910 ddphi += c4 * rbp * cosfac;
2913 dvdl_term += cosfac * rbpBA;
2915 rbpBA = parmB[5] - parmA[5];
2916 ddphi += c5 * rbp * cosfac;
2919 dvdl_term += cosfac * rbpBA;
2921 ddphi = -ddphi * sin_phi; /* 11 */
2923 do_dih_fup<flavor>(ai, aj, ak, al, ddphi, r_ij, r_kj, r_kl, m, n, f, fshift, pbc, x, t1, t2,
2927 *dvdlambda += dvdl_term;
2934 /*! \brief Mysterious undocumented function */
2935 int cmap_setup_grid_index(int ip, int grid_spacing, int* ipm1, int* ipp1, int* ipp2)
2941 ip = ip + grid_spacing - 1;
2943 else if (ip > grid_spacing)
2945 ip = ip - grid_spacing - 1;
2954 im1 = grid_spacing - 1;
2956 else if (ip == grid_spacing - 2)
2960 else if (ip == grid_spacing - 1)
2975 real cmap_dihs(int nbonds,
2976 const t_iatom forceatoms[],
2977 const t_iparams forceparams[],
2978 const gmx_cmap_t* cmap_grid,
2982 const struct t_pbc* pbc,
2983 real gmx_unused lambda,
2984 real gmx_unused* dvdlambda,
2985 const t_mdatoms gmx_unused* md,
2986 t_fcdata gmx_unused* fcd,
2987 int gmx_unused* global_atom_index)
2990 int ai, aj, ak, al, am;
2991 int a1i, a1j, a1k, a1l, a2i, a2j, a2k, a2l;
2993 int t11, t21, t31, t12, t22, t32;
2994 int iphi1, ip1m1, ip1p1, ip1p2;
2995 int iphi2, ip2m1, ip2p1, ip2p2;
2997 int pos1, pos2, pos3, pos4;
2999 real ty[4], ty1[4], ty2[4], ty12[4], tx[16];
3000 real phi1, cos_phi1, sin_phi1, xphi1;
3001 real phi2, cos_phi2, sin_phi2, xphi2;
3002 real dx, tt, tu, e, df1, df2, vtot;
3003 real ra21, rb21, rg21, rg1, rgr1, ra2r1, rb2r1, rabr1;
3004 real ra22, rb22, rg22, rg2, rgr2, ra2r2, rb2r2, rabr2;
3005 real fg1, hg1, fga1, hgb1, gaa1, gbb1;
3006 real fg2, hg2, fga2, hgb2, gaa2, gbb2;
3009 rvec r1_ij, r1_kj, r1_kl, m1, n1;
3010 rvec r2_ij, r2_kj, r2_kl, m2, n2;
3011 rvec f1_i, f1_j, f1_k, f1_l;
3012 rvec f2_i, f2_j, f2_k, f2_l;
3013 rvec a1, b1, a2, b2;
3014 rvec f1, g1, h1, f2, g2, h2;
3015 rvec dtf1, dtg1, dth1, dtf2, dtg2, dth2;
3017 int loop_index[4][4] = { { 0, 4, 8, 12 }, { 1, 5, 9, 13 }, { 2, 6, 10, 14 }, { 3, 7, 11, 15 } };
3019 /* Total CMAP energy */
3022 for (n = 0; n < nbonds;)
3024 /* Five atoms are involved in the two torsions */
3025 type = forceatoms[n++];
3026 ai = forceatoms[n++];
3027 aj = forceatoms[n++];
3028 ak = forceatoms[n++];
3029 al = forceatoms[n++];
3030 am = forceatoms[n++];
3032 /* Which CMAP type is this */
3033 const int cmapA = forceparams[type].cmap.cmapA;
3034 const real* cmapd = cmap_grid->cmapdata[cmapA].cmap.data();
3042 phi1 = dih_angle(x[a1i], x[a1j], x[a1k], x[a1l], pbc, r1_ij, r1_kj, r1_kl, m1, n1, &t11,
3043 &t21, &t31); /* 84 */
3045 cos_phi1 = std::cos(phi1);
3047 a1[0] = r1_ij[1] * r1_kj[2] - r1_ij[2] * r1_kj[1];
3048 a1[1] = r1_ij[2] * r1_kj[0] - r1_ij[0] * r1_kj[2];
3049 a1[2] = r1_ij[0] * r1_kj[1] - r1_ij[1] * r1_kj[0]; /* 9 */
3051 b1[0] = r1_kl[1] * r1_kj[2] - r1_kl[2] * r1_kj[1];
3052 b1[1] = r1_kl[2] * r1_kj[0] - r1_kl[0] * r1_kj[2];
3053 b1[2] = r1_kl[0] * r1_kj[1] - r1_kl[1] * r1_kj[0]; /* 9 */
3055 pbc_rvec_sub(pbc, x[a1l], x[a1k], h1);
3057 ra21 = iprod(a1, a1); /* 5 */
3058 rb21 = iprod(b1, b1); /* 5 */
3059 rg21 = iprod(r1_kj, r1_kj); /* 5 */
3065 rabr1 = sqrt(ra2r1 * rb2r1);
3067 sin_phi1 = rg1 * rabr1 * iprod(a1, h1) * (-1);
3069 if (cos_phi1 < -0.5 || cos_phi1 > 0.5)
3071 phi1 = std::asin(sin_phi1);
3081 phi1 = -M_PI - phi1;
3087 phi1 = std::acos(cos_phi1);
3095 xphi1 = phi1 + M_PI; /* 1 */
3097 /* Second torsion */
3103 phi2 = dih_angle(x[a2i], x[a2j], x[a2k], x[a2l], pbc, r2_ij, r2_kj, r2_kl, m2, n2, &t12,
3104 &t22, &t32); /* 84 */
3106 cos_phi2 = std::cos(phi2);
3108 a2[0] = r2_ij[1] * r2_kj[2] - r2_ij[2] * r2_kj[1];
3109 a2[1] = r2_ij[2] * r2_kj[0] - r2_ij[0] * r2_kj[2];
3110 a2[2] = r2_ij[0] * r2_kj[1] - r2_ij[1] * r2_kj[0]; /* 9 */
3112 b2[0] = r2_kl[1] * r2_kj[2] - r2_kl[2] * r2_kj[1];
3113 b2[1] = r2_kl[2] * r2_kj[0] - r2_kl[0] * r2_kj[2];
3114 b2[2] = r2_kl[0] * r2_kj[1] - r2_kl[1] * r2_kj[0]; /* 9 */
3116 pbc_rvec_sub(pbc, x[a2l], x[a2k], h2);
3118 ra22 = iprod(a2, a2); /* 5 */
3119 rb22 = iprod(b2, b2); /* 5 */
3120 rg22 = iprod(r2_kj, r2_kj); /* 5 */
3126 rabr2 = sqrt(ra2r2 * rb2r2);
3128 sin_phi2 = rg2 * rabr2 * iprod(a2, h2) * (-1);
3130 if (cos_phi2 < -0.5 || cos_phi2 > 0.5)
3132 phi2 = std::asin(sin_phi2);
3142 phi2 = -M_PI - phi2;
3148 phi2 = std::acos(cos_phi2);
3156 xphi2 = phi2 + M_PI; /* 1 */
3158 /* Range mangling */
3161 xphi1 = xphi1 + 2 * M_PI;
3163 else if (xphi1 >= 2 * M_PI)
3165 xphi1 = xphi1 - 2 * M_PI;
3170 xphi2 = xphi2 + 2 * M_PI;
3172 else if (xphi2 >= 2 * M_PI)
3174 xphi2 = xphi2 - 2 * M_PI;
3177 /* Number of grid points */
3178 dx = 2 * M_PI / cmap_grid->grid_spacing;
3180 /* Where on the grid are we */
3181 iphi1 = static_cast<int>(xphi1 / dx);
3182 iphi2 = static_cast<int>(xphi2 / dx);
3184 iphi1 = cmap_setup_grid_index(iphi1, cmap_grid->grid_spacing, &ip1m1, &ip1p1, &ip1p2);
3185 iphi2 = cmap_setup_grid_index(iphi2, cmap_grid->grid_spacing, &ip2m1, &ip2p1, &ip2p2);
3187 pos1 = iphi1 * cmap_grid->grid_spacing + iphi2;
3188 pos2 = ip1p1 * cmap_grid->grid_spacing + iphi2;
3189 pos3 = ip1p1 * cmap_grid->grid_spacing + ip2p1;
3190 pos4 = iphi1 * cmap_grid->grid_spacing + ip2p1;
3192 ty[0] = cmapd[pos1 * 4];
3193 ty[1] = cmapd[pos2 * 4];
3194 ty[2] = cmapd[pos3 * 4];
3195 ty[3] = cmapd[pos4 * 4];
3197 ty1[0] = cmapd[pos1 * 4 + 1];
3198 ty1[1] = cmapd[pos2 * 4 + 1];
3199 ty1[2] = cmapd[pos3 * 4 + 1];
3200 ty1[3] = cmapd[pos4 * 4 + 1];
3202 ty2[0] = cmapd[pos1 * 4 + 2];
3203 ty2[1] = cmapd[pos2 * 4 + 2];
3204 ty2[2] = cmapd[pos3 * 4 + 2];
3205 ty2[3] = cmapd[pos4 * 4 + 2];
3207 ty12[0] = cmapd[pos1 * 4 + 3];
3208 ty12[1] = cmapd[pos2 * 4 + 3];
3209 ty12[2] = cmapd[pos3 * 4 + 3];
3210 ty12[3] = cmapd[pos4 * 4 + 3];
3212 /* Switch to degrees */
3213 dx = 360.0 / cmap_grid->grid_spacing;
3214 xphi1 = xphi1 * RAD2DEG;
3215 xphi2 = xphi2 * RAD2DEG;
3217 for (i = 0; i < 4; i++) /* 16 */
3220 tx[i + 4] = ty1[i] * dx;
3221 tx[i + 8] = ty2[i] * dx;
3222 tx[i + 12] = ty12[i] * dx * dx;
3225 real tc[16] = { 0 };
3226 for (int idx = 0; idx < 16; idx++) /* 1056 */
3228 for (int k = 0; k < 16; k++)
3230 tc[idx] += cmap_coeff_matrix[k * 16 + idx] * tx[k];
3234 tt = (xphi1 - iphi1 * dx) / dx;
3235 tu = (xphi2 - iphi2 * dx) / dx;
3241 for (i = 3; i >= 0; i--)
3243 l1 = loop_index[i][3];
3244 l2 = loop_index[i][2];
3245 l3 = loop_index[i][1];
3247 e = tt * e + ((tc[i * 4 + 3] * tu + tc[i * 4 + 2]) * tu + tc[i * 4 + 1]) * tu + tc[i * 4];
3248 df1 = tu * df1 + (3.0 * tc[l1] * tt + 2.0 * tc[l2]) * tt + tc[l3];
3249 df2 = tt * df2 + (3.0 * tc[i * 4 + 3] * tu + 2.0 * tc[i * 4 + 2]) * tu + tc[i * 4 + 1];
3259 /* Do forces - first torsion */
3260 fg1 = iprod(r1_ij, r1_kj);
3261 hg1 = iprod(r1_kl, r1_kj);
3262 fga1 = fg1 * ra2r1 * rgr1;
3263 hgb1 = hg1 * rb2r1 * rgr1;
3264 gaa1 = -ra2r1 * rg1;
3267 for (i = 0; i < DIM; i++)
3269 dtf1[i] = gaa1 * a1[i];
3270 dtg1[i] = fga1 * a1[i] - hgb1 * b1[i];
3271 dth1[i] = gbb1 * b1[i];
3273 f1[i] = df1 * dtf1[i];
3274 g1[i] = df1 * dtg1[i];
3275 h1[i] = df1 * dth1[i];
3278 f1_j[i] = -f1[i] - g1[i];
3279 f1_k[i] = h1[i] + g1[i];
3282 f[a1i][i] = f[a1i][i] + f1_i[i];
3283 f[a1j][i] = f[a1j][i] + f1_j[i]; /* - f1[i] - g1[i] */
3284 f[a1k][i] = f[a1k][i] + f1_k[i]; /* h1[i] + g1[i] */
3285 f[a1l][i] = f[a1l][i] + f1_l[i]; /* h1[i] */
3288 /* Do forces - second torsion */
3289 fg2 = iprod(r2_ij, r2_kj);
3290 hg2 = iprod(r2_kl, r2_kj);
3291 fga2 = fg2 * ra2r2 * rgr2;
3292 hgb2 = hg2 * rb2r2 * rgr2;
3293 gaa2 = -ra2r2 * rg2;
3296 for (i = 0; i < DIM; i++)
3298 dtf2[i] = gaa2 * a2[i];
3299 dtg2[i] = fga2 * a2[i] - hgb2 * b2[i];
3300 dth2[i] = gbb2 * b2[i];
3302 f2[i] = df2 * dtf2[i];
3303 g2[i] = df2 * dtg2[i];
3304 h2[i] = df2 * dth2[i];
3307 f2_j[i] = -f2[i] - g2[i];
3308 f2_k[i] = h2[i] + g2[i];
3311 f[a2i][i] = f[a2i][i] + f2_i[i]; /* f2[i] */
3312 f[a2j][i] = f[a2j][i] + f2_j[i]; /* - f2[i] - g2[i] */
3313 f[a2k][i] = f[a2k][i] + f2_k[i]; /* h2[i] + g2[i] */
3314 f[a2l][i] = f[a2l][i] + f2_l[i]; /* - h2[i] */
3318 if (fshift != nullptr)
3322 t31 = pbc_rvec_sub(pbc, x[a1l], x[a1j], h1);
3323 t32 = pbc_rvec_sub(pbc, x[a2l], x[a2j], h2);
3331 rvec_inc(fshift[t11], f1_i);
3332 rvec_inc(fshift[CENTRAL], f1_j);
3333 rvec_inc(fshift[t21], f1_k);
3334 rvec_inc(fshift[t31], f1_l);
3336 rvec_inc(fshift[t12], f2_i);
3337 rvec_inc(fshift[CENTRAL], f2_j);
3338 rvec_inc(fshift[t22], f2_k);
3339 rvec_inc(fshift[t32], f2_l);
3349 /***********************************************************
3351 * G R O M O S 9 6 F U N C T I O N S
3353 ***********************************************************/
3354 real g96harmonic(real kA, real kB, real xA, real xB, real x, real lambda, real* V, real* F)
3356 const real half = 0.5;
3357 real L1, kk, x0, dx, dx2;
3358 real v, f, dvdlambda;
3361 kk = L1 * kA + lambda * kB;
3362 x0 = L1 * xA + lambda * xB;
3368 v = half * kk * dx2;
3369 dvdlambda = half * (kB - kA) * dx2 + (xA - xB) * kk * dx;
3376 /* That was 21 flops */
3379 template<BondedKernelFlavor flavor>
3380 real g96bonds(int nbonds,
3381 const t_iatom forceatoms[],
3382 const t_iparams forceparams[],
3389 const t_mdatoms gmx_unused* md,
3390 t_fcdata gmx_unused* fcd,
3391 int gmx_unused* global_atom_index)
3393 int i, ki, ai, aj, type;
3394 real dr2, fbond, vbond, vtot;
3398 for (i = 0; (i < nbonds);)
3400 type = forceatoms[i++];
3401 ai = forceatoms[i++];
3402 aj = forceatoms[i++];
3404 ki = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /* 3 */
3405 dr2 = iprod(dx, dx); /* 5 */
3407 *dvdlambda += g96harmonic(forceparams[type].harmonic.krA, forceparams[type].harmonic.krB,
3408 forceparams[type].harmonic.rA, forceparams[type].harmonic.rB, dr2,
3409 lambda, &vbond, &fbond);
3411 vtot += 0.5 * vbond; /* 1*/
3413 spreadBondForces<flavor>(fbond, dx, ai, aj, f, ki, fshift); /* 15 */
3418 real g96bond_angle(const rvec xi, const rvec xj, const rvec xk, const t_pbc* pbc, rvec r_ij, rvec r_kj, int* t1, int* t2)
3419 /* Return value is the angle between the bonds i-j and j-k */
3423 *t1 = pbc_rvec_sub(pbc, xi, xj, r_ij); /* 3 */
3424 *t2 = pbc_rvec_sub(pbc, xk, xj, r_kj); /* 3 */
3426 costh = cos_angle(r_ij, r_kj); /* 25 */
3431 template<BondedKernelFlavor flavor>
3432 real g96angles(int nbonds,
3433 const t_iatom forceatoms[],
3434 const t_iparams forceparams[],
3441 const t_mdatoms gmx_unused* md,
3442 t_fcdata gmx_unused* fcd,
3443 int gmx_unused* global_atom_index)
3445 int i, ai, aj, ak, type, m, t1, t2;
3447 real cos_theta, dVdt, va, vtot;
3448 real rij_1, rij_2, rkj_1, rkj_2, rijrkj_1;
3452 for (i = 0; (i < nbonds);)
3454 type = forceatoms[i++];
3455 ai = forceatoms[i++];
3456 aj = forceatoms[i++];
3457 ak = forceatoms[i++];
3459 cos_theta = g96bond_angle(x[ai], x[aj], x[ak], pbc, r_ij, r_kj, &t1, &t2);
3461 *dvdlambda += g96harmonic(forceparams[type].harmonic.krA, forceparams[type].harmonic.krB,
3462 forceparams[type].harmonic.rA, forceparams[type].harmonic.rB,
3463 cos_theta, lambda, &va, &dVdt);
3466 rij_1 = gmx::invsqrt(iprod(r_ij, r_ij));
3467 rkj_1 = gmx::invsqrt(iprod(r_kj, r_kj));
3468 rij_2 = rij_1 * rij_1;
3469 rkj_2 = rkj_1 * rkj_1;
3470 rijrkj_1 = rij_1 * rkj_1; /* 23 */
3472 for (m = 0; (m < DIM); m++) /* 42 */
3474 f_i[m] = dVdt * (r_kj[m] * rijrkj_1 - r_ij[m] * rij_2 * cos_theta);
3475 f_k[m] = dVdt * (r_ij[m] * rijrkj_1 - r_kj[m] * rkj_2 * cos_theta);
3476 f_j[m] = -f_i[m] - f_k[m];
3482 if (computeVirial(flavor))
3484 rvec_inc(fshift[t1], f_i);
3485 rvec_inc(fshift[CENTRAL], f_j);
3486 rvec_inc(fshift[t2], f_k); /* 9 */
3493 template<BondedKernelFlavor flavor>
3494 real cross_bond_bond(int nbonds,
3495 const t_iatom forceatoms[],
3496 const t_iparams forceparams[],
3501 real gmx_unused lambda,
3502 real gmx_unused* dvdlambda,
3503 const t_mdatoms gmx_unused* md,
3504 t_fcdata gmx_unused* fcd,
3505 int gmx_unused* global_atom_index)
3507 /* Potential from Lawrence and Skimmer, Chem. Phys. Lett. 372 (2003)
3510 int i, ai, aj, ak, type, m, t1, t2;
3512 real vtot, vrr, s1, s2, r1, r2, r1e, r2e, krr;
3516 for (i = 0; (i < nbonds);)
3518 type = forceatoms[i++];
3519 ai = forceatoms[i++];
3520 aj = forceatoms[i++];
3521 ak = forceatoms[i++];
3522 r1e = forceparams[type].cross_bb.r1e;
3523 r2e = forceparams[type].cross_bb.r2e;
3524 krr = forceparams[type].cross_bb.krr;
3526 /* Compute distance vectors ... */
3527 t1 = pbc_rvec_sub(pbc, x[ai], x[aj], r_ij);
3528 t2 = pbc_rvec_sub(pbc, x[ak], x[aj], r_kj);
3530 /* ... and their lengths */
3534 /* Deviations from ideality */
3538 /* Energy (can be negative!) */
3539 vrr = krr * s1 * s2;
3543 svmul(-krr * s2 / r1, r_ij, f_i);
3544 svmul(-krr * s1 / r2, r_kj, f_k);
3546 for (m = 0; (m < DIM); m++) /* 12 */
3548 f_j[m] = -f_i[m] - f_k[m];
3554 if (computeVirial(flavor))
3556 rvec_inc(fshift[t1], f_i);
3557 rvec_inc(fshift[CENTRAL], f_j);
3558 rvec_inc(fshift[t2], f_k); /* 9 */
3565 template<BondedKernelFlavor flavor>
3566 real cross_bond_angle(int nbonds,
3567 const t_iatom forceatoms[],
3568 const t_iparams forceparams[],
3573 real gmx_unused lambda,
3574 real gmx_unused* dvdlambda,
3575 const t_mdatoms gmx_unused* md,
3576 t_fcdata gmx_unused* fcd,
3577 int gmx_unused* global_atom_index)
3579 /* Potential from Lawrence and Skimmer, Chem. Phys. Lett. 372 (2003)
3582 int i, ai, aj, ak, type, m, t1, t2;
3583 rvec r_ij, r_kj, r_ik;
3584 real vtot, vrt, s1, s2, s3, r1, r2, r3, r1e, r2e, r3e, krt, k1, k2, k3;
3588 for (i = 0; (i < nbonds);)
3590 type = forceatoms[i++];
3591 ai = forceatoms[i++];
3592 aj = forceatoms[i++];
3593 ak = forceatoms[i++];
3594 r1e = forceparams[type].cross_ba.r1e;
3595 r2e = forceparams[type].cross_ba.r2e;
3596 r3e = forceparams[type].cross_ba.r3e;
3597 krt = forceparams[type].cross_ba.krt;
3599 /* Compute distance vectors ... */
3600 t1 = pbc_rvec_sub(pbc, x[ai], x[aj], r_ij);
3601 t2 = pbc_rvec_sub(pbc, x[ak], x[aj], r_kj);
3602 pbc_rvec_sub(pbc, x[ai], x[ak], r_ik);
3604 /* ... and their lengths */
3609 /* Deviations from ideality */
3614 /* Energy (can be negative!) */
3615 vrt = krt * s3 * (s1 + s2);
3619 k1 = -krt * (s3 / r1);
3620 k2 = -krt * (s3 / r2);
3621 k3 = -krt * (s1 + s2) / r3;
3622 for (m = 0; (m < DIM); m++)
3624 f_i[m] = k1 * r_ij[m] + k3 * r_ik[m];
3625 f_k[m] = k2 * r_kj[m] - k3 * r_ik[m];
3626 f_j[m] = -f_i[m] - f_k[m];
3629 for (m = 0; (m < DIM); m++) /* 12 */
3636 if (computeVirial(flavor))
3638 rvec_inc(fshift[t1], f_i);
3639 rvec_inc(fshift[CENTRAL], f_j);
3640 rvec_inc(fshift[t2], f_k); /* 9 */
3647 /*! \brief Computes the potential and force for a tabulated potential */
3648 real bonded_tab(const char* type,
3650 const bondedtable_t* table,
3658 real k, tabscale, rt, eps, eps2, Yt, Ft, Geps, Heps2, Fp, VV, FF;
3662 k = (1.0 - lambda) * kA + lambda * kB;
3664 tabscale = table->scale;
3665 const real* VFtab = table->data.data();
3668 n0 = static_cast<int>(rt);
3672 "A tabulated %s interaction table number %d is out of the table range: r %f, "
3673 "between table indices %d and %d, table length %d",
3674 type, table_nr, r, n0, n0 + 1, table->n);
3680 Ft = VFtab[nnn + 1];
3681 Geps = VFtab[nnn + 2] * eps;
3682 Heps2 = VFtab[nnn + 3] * eps2;
3683 Fp = Ft + Geps + Heps2;
3685 FF = Fp + Geps + 2.0 * Heps2;
3687 *F = -k * FF * tabscale;
3689 dvdlambda = (kB - kA) * VV;
3693 /* That was 22 flops */
3696 template<BondedKernelFlavor flavor>
3697 real tab_bonds(int nbonds,
3698 const t_iatom forceatoms[],
3699 const t_iparams forceparams[],
3706 const t_mdatoms gmx_unused* md,
3708 int gmx_unused* global_atom_index)
3710 int i, ki, ai, aj, type, table;
3711 real dr, dr2, fbond, vbond, vtot;
3715 for (i = 0; (i < nbonds);)
3717 type = forceatoms[i++];
3718 ai = forceatoms[i++];
3719 aj = forceatoms[i++];
3721 ki = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /* 3 */
3722 dr2 = iprod(dx, dx); /* 5 */
3723 dr = dr2 * gmx::invsqrt(dr2); /* 10 */
3725 table = forceparams[type].tab.table;
3727 *dvdlambda += bonded_tab("bond", table, &fcd->bondtab[table], forceparams[type].tab.kA,
3728 forceparams[type].tab.kB, dr, lambda, &vbond, &fbond); /* 22 */
3736 vtot += vbond; /* 1*/
3737 fbond *= gmx::invsqrt(dr2); /* 6 */
3739 spreadBondForces<flavor>(fbond, dx, ai, aj, f, ki, fshift); /* 15 */
3744 template<BondedKernelFlavor flavor>
3745 real tab_angles(int nbonds,
3746 const t_iatom forceatoms[],
3747 const t_iparams forceparams[],
3754 const t_mdatoms gmx_unused* md,
3756 int gmx_unused* global_atom_index)
3758 int i, ai, aj, ak, t1, t2, type, table;
3760 real cos_theta, cos_theta2, theta, dVdt, va, vtot;
3763 for (i = 0; (i < nbonds);)
3765 type = forceatoms[i++];
3766 ai = forceatoms[i++];
3767 aj = forceatoms[i++];
3768 ak = forceatoms[i++];
3770 theta = bond_angle(x[ai], x[aj], x[ak], pbc, r_ij, r_kj, &cos_theta, &t1, &t2); /* 41 */
3772 table = forceparams[type].tab.table;
3774 *dvdlambda += bonded_tab("angle", table, &fcd->angletab[table], forceparams[type].tab.kA,
3775 forceparams[type].tab.kB, theta, lambda, &va, &dVdt); /* 22 */
3778 cos_theta2 = gmx::square(cos_theta); /* 1 */
3787 st = dVdt * gmx::invsqrt(1 - cos_theta2); /* 12 */
3788 sth = st * cos_theta; /* 1 */
3789 nrkj2 = iprod(r_kj, r_kj); /* 5 */
3790 nrij2 = iprod(r_ij, r_ij);
3792 cik = st * gmx::invsqrt(nrkj2 * nrij2); /* 12 */
3793 cii = sth / nrij2; /* 10 */
3794 ckk = sth / nrkj2; /* 10 */
3796 for (m = 0; (m < DIM); m++) /* 39 */
3798 f_i[m] = -(cik * r_kj[m] - cii * r_ij[m]);
3799 f_k[m] = -(cik * r_ij[m] - ckk * r_kj[m]);
3800 f_j[m] = -f_i[m] - f_k[m];
3806 if (computeVirial(flavor))
3808 rvec_inc(fshift[t1], f_i);
3809 rvec_inc(fshift[CENTRAL], f_j);
3810 rvec_inc(fshift[t2], f_k);
3817 template<BondedKernelFlavor flavor>
3818 real tab_dihs(int nbonds,
3819 const t_iatom forceatoms[],
3820 const t_iparams forceparams[],
3827 const t_mdatoms gmx_unused* md,
3829 int gmx_unused* global_atom_index)
3831 int i, type, ai, aj, ak, al, table;
3833 rvec r_ij, r_kj, r_kl, m, n;
3834 real phi, ddphi, vpd, vtot;
3837 for (i = 0; (i < nbonds);)
3839 type = forceatoms[i++];
3840 ai = forceatoms[i++];
3841 aj = forceatoms[i++];
3842 ak = forceatoms[i++];
3843 al = forceatoms[i++];
3845 phi = dih_angle(x[ai], x[aj], x[ak], x[al], pbc, r_ij, r_kj, r_kl, m, n, &t1, &t2, &t3); /* 84 */
3847 table = forceparams[type].tab.table;
3849 /* Hopefully phi+M_PI never results in values < 0 */
3850 *dvdlambda += bonded_tab("dihedral", table, &fcd->dihtab[table], forceparams[type].tab.kA,
3851 forceparams[type].tab.kB, phi + M_PI, lambda, &vpd, &ddphi);
3854 do_dih_fup<flavor>(ai, aj, ak, al, -ddphi, r_ij, r_kj, r_kl, m, n, f, fshift, pbc, x, t1,
3862 struct BondedInteractions
3864 BondedFunction function;
3868 // Bug in old clang versions prevents constexpr. constexpr is needed for MSVC.
3869 #if defined(__clang__) && __clang_major__ < 6
3870 # define CONSTEXPR_EXCL_OLD_CLANG const
3872 # define CONSTEXPR_EXCL_OLD_CLANG constexpr
3875 /*! \brief Lookup table of bonded interaction functions
3877 * This must have as many entries as interaction_function in ifunc.cpp */
3878 template<BondedKernelFlavor flavor>
3879 CONSTEXPR_EXCL_OLD_CLANG std::array<BondedInteractions, F_NRE> c_bondedInteractionFunctions = {
3880 BondedInteractions{ bonds<flavor>, eNR_BONDS }, // F_BONDS
3881 BondedInteractions{ g96bonds<flavor>, eNR_BONDS }, // F_G96BONDS
3882 BondedInteractions{ morse_bonds<flavor>, eNR_MORSE }, // F_MORSE
3883 BondedInteractions{ cubic_bonds<flavor>, eNR_CUBICBONDS }, // F_CUBICBONDS
3884 BondedInteractions{ unimplemented, -1 }, // F_CONNBONDS
3885 BondedInteractions{ bonds<flavor>, eNR_BONDS }, // F_HARMONIC
3886 BondedInteractions{ FENE_bonds<flavor>, eNR_FENEBONDS }, // F_FENEBONDS
3887 BondedInteractions{ tab_bonds<flavor>, eNR_TABBONDS }, // F_TABBONDS
3888 BondedInteractions{ tab_bonds<flavor>, eNR_TABBONDS }, // F_TABBONDSNC
3889 BondedInteractions{ restraint_bonds<flavor>, eNR_RESTRBONDS }, // F_RESTRBONDS
3890 BondedInteractions{ angles<flavor>, eNR_ANGLES }, // F_ANGLES
3891 BondedInteractions{ g96angles<flavor>, eNR_ANGLES }, // F_G96ANGLES
3892 BondedInteractions{ restrangles<flavor>, eNR_ANGLES }, // F_RESTRANGLES
3893 BondedInteractions{ linear_angles<flavor>, eNR_ANGLES }, // F_LINEAR_ANGLES
3894 BondedInteractions{ cross_bond_bond<flavor>, eNR_CROSS_BOND_BOND }, // F_CROSS_BOND_BONDS
3895 BondedInteractions{ cross_bond_angle<flavor>, eNR_CROSS_BOND_ANGLE }, // F_CROSS_BOND_ANGLES
3896 BondedInteractions{ urey_bradley<flavor>, eNR_UREY_BRADLEY }, // F_UREY_BRADLEY
3897 BondedInteractions{ quartic_angles<flavor>, eNR_QANGLES }, // F_QUARTIC_ANGLES
3898 BondedInteractions{ tab_angles<flavor>, eNR_TABANGLES }, // F_TABANGLES
3899 BondedInteractions{ pdihs<flavor>, eNR_PROPER }, // F_PDIHS
3900 BondedInteractions{ rbdihs<flavor>, eNR_RB }, // F_RBDIHS
3901 BondedInteractions{ restrdihs<flavor>, eNR_PROPER }, // F_RESTRDIHS
3902 BondedInteractions{ cbtdihs<flavor>, eNR_RB }, // F_CBTDIHS
3903 BondedInteractions{ rbdihs<flavor>, eNR_FOURDIH }, // F_FOURDIHS
3904 BondedInteractions{ idihs<flavor>, eNR_IMPROPER }, // F_IDIHS
3905 BondedInteractions{ pdihs<flavor>, eNR_IMPROPER }, // F_PIDIHS
3906 BondedInteractions{ tab_dihs<flavor>, eNR_TABDIHS }, // F_TABDIHS
3907 BondedInteractions{ unimplemented, eNR_CMAP }, // F_CMAP
3908 BondedInteractions{ unimplemented, -1 }, // F_GB12_NOLONGERUSED
3909 BondedInteractions{ unimplemented, -1 }, // F_GB13_NOLONGERUSED
3910 BondedInteractions{ unimplemented, -1 }, // F_GB14_NOLONGERUSED
3911 BondedInteractions{ unimplemented, -1 }, // F_GBPOL_NOLONGERUSED
3912 BondedInteractions{ unimplemented, -1 }, // F_NPSOLVATION_NOLONGERUSED
3913 BondedInteractions{ unimplemented, eNR_NB14 }, // F_LJ14
3914 BondedInteractions{ unimplemented, -1 }, // F_COUL14
3915 BondedInteractions{ unimplemented, eNR_NB14 }, // F_LJC14_Q
3916 BondedInteractions{ unimplemented, eNR_NB14 }, // F_LJC_PAIRS_NB
3917 BondedInteractions{ unimplemented, -1 }, // F_LJ
3918 BondedInteractions{ unimplemented, -1 }, // F_BHAM
3919 BondedInteractions{ unimplemented, -1 }, // F_LJ_LR_NOLONGERUSED
3920 BondedInteractions{ unimplemented, -1 }, // F_BHAM_LR_NOLONGERUSED
3921 BondedInteractions{ unimplemented, -1 }, // F_DISPCORR
3922 BondedInteractions{ unimplemented, -1 }, // F_COUL_SR
3923 BondedInteractions{ unimplemented, -1 }, // F_COUL_LR_NOLONGERUSED
3924 BondedInteractions{ unimplemented, -1 }, // F_RF_EXCL
3925 BondedInteractions{ unimplemented, -1 }, // F_COUL_RECIP
3926 BondedInteractions{ unimplemented, -1 }, // F_LJ_RECIP
3927 BondedInteractions{ unimplemented, -1 }, // F_DPD
3928 BondedInteractions{ polarize<flavor>, eNR_POLARIZE }, // F_POLARIZATION
3929 BondedInteractions{ water_pol<flavor>, eNR_WPOL }, // F_WATER_POL
3930 BondedInteractions{ thole_pol<flavor>, eNR_THOLE }, // F_THOLE_POL
3931 BondedInteractions{ anharm_polarize<flavor>, eNR_ANHARM_POL }, // F_ANHARM_POL
3932 BondedInteractions{ unimplemented, -1 }, // F_POSRES
3933 BondedInteractions{ unimplemented, -1 }, // F_FBPOSRES
3934 BondedInteractions{ ta_disres, eNR_DISRES }, // F_DISRES
3935 BondedInteractions{ unimplemented, -1 }, // F_DISRESVIOL
3936 BondedInteractions{ orires, eNR_ORIRES }, // F_ORIRES
3937 BondedInteractions{ unimplemented, -1 }, // F_ORIRESDEV
3938 BondedInteractions{ angres<flavor>, eNR_ANGRES }, // F_ANGRES
3939 BondedInteractions{ angresz<flavor>, eNR_ANGRESZ }, // F_ANGRESZ
3940 BondedInteractions{ dihres<flavor>, eNR_DIHRES }, // F_DIHRES
3941 BondedInteractions{ unimplemented, -1 }, // F_DIHRESVIOL
3942 BondedInteractions{ unimplemented, -1 }, // F_CONSTR
3943 BondedInteractions{ unimplemented, -1 }, // F_CONSTRNC
3944 BondedInteractions{ unimplemented, -1 }, // F_SETTLE
3945 BondedInteractions{ unimplemented, -1 }, // F_VSITE2
3946 BondedInteractions{ unimplemented, -1 }, // F_VSITE3
3947 BondedInteractions{ unimplemented, -1 }, // F_VSITE3FD
3948 BondedInteractions{ unimplemented, -1 }, // F_VSITE3FAD
3949 BondedInteractions{ unimplemented, -1 }, // F_VSITE3OUT
3950 BondedInteractions{ unimplemented, -1 }, // F_VSITE4FD
3951 BondedInteractions{ unimplemented, -1 }, // F_VSITE4FDN
3952 BondedInteractions{ unimplemented, -1 }, // F_VSITEN
3953 BondedInteractions{ unimplemented, -1 }, // F_COM_PULL
3954 BondedInteractions{ unimplemented, -1 }, // F_DENSITYFITTING
3955 BondedInteractions{ unimplemented, -1 }, // F_EQM
3956 BondedInteractions{ unimplemented, -1 }, // F_EPOT
3957 BondedInteractions{ unimplemented, -1 }, // F_EKIN
3958 BondedInteractions{ unimplemented, -1 }, // F_ETOT
3959 BondedInteractions{ unimplemented, -1 }, // F_ECONSERVED
3960 BondedInteractions{ unimplemented, -1 }, // F_TEMP
3961 BondedInteractions{ unimplemented, -1 }, // F_VTEMP_NOLONGERUSED
3962 BondedInteractions{ unimplemented, -1 }, // F_PDISPCORR
3963 BondedInteractions{ unimplemented, -1 }, // F_PRES
3964 BondedInteractions{ unimplemented, -1 }, // F_DVDL_CONSTR
3965 BondedInteractions{ unimplemented, -1 }, // F_DVDL
3966 BondedInteractions{ unimplemented, -1 }, // F_DKDL
3967 BondedInteractions{ unimplemented, -1 }, // F_DVDL_COUL
3968 BondedInteractions{ unimplemented, -1 }, // F_DVDL_VDW
3969 BondedInteractions{ unimplemented, -1 }, // F_DVDL_BONDED
3970 BondedInteractions{ unimplemented, -1 }, // F_DVDL_RESTRAINT
3971 BondedInteractions{ unimplemented, -1 }, // F_DVDL_TEMPERATURE
3974 /*! \brief List of instantiated BondedInteractions list */
3975 CONSTEXPR_EXCL_OLD_CLANG
3976 gmx::EnumerationArray<BondedKernelFlavor, std::array<BondedInteractions, F_NRE>> c_bondedInteractionFunctionsPerFlavor = {
3977 c_bondedInteractionFunctions<BondedKernelFlavor::ForcesSimdWhenAvailable>,
3978 c_bondedInteractionFunctions<BondedKernelFlavor::ForcesNoSimd>,
3979 c_bondedInteractionFunctions<BondedKernelFlavor::ForcesAndVirialAndEnergy>,
3980 c_bondedInteractionFunctions<BondedKernelFlavor::ForcesAndEnergy>
3987 real calculateSimpleBond(const int ftype,
3988 const int numForceatoms,
3989 const t_iatom forceatoms[],
3990 const t_iparams forceparams[],
3994 const struct t_pbc* pbc,
3997 const t_mdatoms* md,
3999 int gmx_unused* global_atom_index,
4000 const BondedKernelFlavor bondedKernelFlavor)
4002 const BondedInteractions& bonded = c_bondedInteractionFunctionsPerFlavor[bondedKernelFlavor][ftype];
4004 real v = bonded.function(numForceatoms, forceatoms, forceparams, x, f, fshift, pbc, lambda,
4005 dvdlambda, md, fcd, global_atom_index);
4010 int nrnbIndex(int ftype)
4012 return c_bondedInteractionFunctions<BondedKernelFlavor::ForcesAndVirialAndEnergy>[ftype].nrnbIndex;