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,2018, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
39 * \brief This file defines low-level functions necessary for
40 * computing energies and forces for listed interactions.
42 * \author Mark Abraham <mark.j.abraham@gmail.com>
44 * \ingroup module_listed-forces
57 #include "gromacs/listed-forces/pairs.h"
58 #include "gromacs/math/functions.h"
59 #include "gromacs/math/units.h"
60 #include "gromacs/math/utilities.h"
61 #include "gromacs/math/vec.h"
62 #include "gromacs/mdtypes/fcdata.h"
63 #include "gromacs/pbcutil/ishift.h"
64 #include "gromacs/pbcutil/mshift.h"
65 #include "gromacs/pbcutil/pbc.h"
66 #include "gromacs/pbcutil/pbc-simd.h"
67 #include "gromacs/simd/simd.h"
68 #include "gromacs/simd/simd_math.h"
69 #include "gromacs/simd/vector_operations.h"
70 #include "gromacs/utility/basedefinitions.h"
71 #include "gromacs/utility/fatalerror.h"
72 #include "gromacs/utility/real.h"
73 #include "gromacs/utility/smalloc.h"
75 #include "listed-internal.h"
78 using namespace gmx; // TODO: Remove when this file is moved into gmx namespace
80 /*! \brief Mysterious CMAP coefficient matrix */
81 const int cmap_coeff_matrix[] = {
82 1, 0, -3, 2, 0, 0, 0, 0, -3, 0, 9, -6, 2, 0, -6, 4,
83 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, -9, 6, -2, 0, 6, -4,
84 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 9, -6, 0, 0, -6, 4,
85 0, 0, 3, -2, 0, 0, 0, 0, 0, 0, -9, 6, 0, 0, 6, -4,
86 0, 0, 0, 0, 1, 0, -3, 2, -2, 0, 6, -4, 1, 0, -3, 2,
87 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 3, -2, 1, 0, -3, 2,
88 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 2, 0, 0, 3, -2,
89 0, 0, 0, 0, 0, 0, 3, -2, 0, 0, -6, 4, 0, 0, 3, -2,
90 0, 1, -2, 1, 0, 0, 0, 0, 0, -3, 6, -3, 0, 2, -4, 2,
91 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, -6, 3, 0, -2, 4, -2,
92 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 3, 0, 0, 2, -2,
93 0, 0, -1, 1, 0, 0, 0, 0, 0, 0, 3, -3, 0, 0, -2, 2,
94 0, 0, 0, 0, 0, 1, -2, 1, 0, -2, 4, -2, 0, 1, -2, 1,
95 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 2, -1, 0, 1, -2, 1,
96 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, -1, 1,
97 0, 0, 0, 0, 0, 0, -1, 1, 0, 0, 2, -2, 0, 0, -1, 1
101 /*! \brief Compute dx = xi - xj, modulo PBC if non-NULL
103 * \todo This kind of code appears in many places. Consolidate it */
104 static int pbc_rvec_sub(const t_pbc *pbc, const rvec xi, const rvec xj, rvec dx)
108 return pbc_dx_aiuc(pbc, xi, xj, dx);
112 rvec_sub(xi, xj, dx);
117 /*! \brief Morse potential bond
119 * By Frank Everdij. Three parameters needed:
121 * b0 = equilibrium distance in nm
122 * be = beta in nm^-1 (actually, it's nu_e*Sqrt(2*pi*pi*mu/D_e))
123 * cb = well depth in kJ/mol
125 * Note: the potential is referenced to be +cb at infinite separation
126 * and zero at the equilibrium distance!
128 real morse_bonds(int nbonds,
129 const t_iatom forceatoms[], const t_iparams forceparams[],
130 const rvec x[], rvec4 f[], rvec fshift[],
131 const t_pbc *pbc, const t_graph *g,
132 real lambda, real *dvdlambda,
133 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
134 int gmx_unused *global_atom_index)
136 const real one = 1.0;
137 const real two = 2.0;
138 real dr, dr2, temp, omtemp, cbomtemp, fbond, vbond, fij, vtot;
139 real b0, be, cb, b0A, beA, cbA, b0B, beB, cbB, L1;
141 int i, m, ki, type, ai, aj;
145 for (i = 0; (i < nbonds); )
147 type = forceatoms[i++];
148 ai = forceatoms[i++];
149 aj = forceatoms[i++];
151 b0A = forceparams[type].morse.b0A;
152 beA = forceparams[type].morse.betaA;
153 cbA = forceparams[type].morse.cbA;
155 b0B = forceparams[type].morse.b0B;
156 beB = forceparams[type].morse.betaB;
157 cbB = forceparams[type].morse.cbB;
159 L1 = one-lambda; /* 1 */
160 b0 = L1*b0A + lambda*b0B; /* 3 */
161 be = L1*beA + lambda*beB; /* 3 */
162 cb = L1*cbA + lambda*cbB; /* 3 */
164 ki = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /* 3 */
165 dr2 = iprod(dx, dx); /* 5 */
166 dr = dr2*gmx::invsqrt(dr2); /* 10 */
167 temp = std::exp(-be*(dr-b0)); /* 12 */
171 /* bonds are constrainted. This may _not_ include bond constraints if they are lambda dependent */
172 *dvdlambda += cbB-cbA;
176 omtemp = one-temp; /* 1 */
177 cbomtemp = cb*omtemp; /* 1 */
178 vbond = cbomtemp*omtemp; /* 1 */
179 fbond = -two*be*temp*cbomtemp*gmx::invsqrt(dr2); /* 9 */
180 vtot += vbond; /* 1 */
182 *dvdlambda += (cbB - cbA) * omtemp * omtemp - (2-2*omtemp)*omtemp * cb * ((b0B-b0A)*be - (beB-beA)*(dr-b0)); /* 15 */
186 ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, aj), dt);
190 for (m = 0; (m < DIM); m++) /* 15 */
195 fshift[ki][m] += fij;
196 fshift[CENTRAL][m] -= fij;
203 real cubic_bonds(int nbonds,
204 const t_iatom forceatoms[], const t_iparams forceparams[],
205 const rvec x[], rvec4 f[], rvec fshift[],
206 const t_pbc *pbc, const t_graph *g,
207 real gmx_unused lambda, real gmx_unused *dvdlambda,
208 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
209 int gmx_unused *global_atom_index)
211 const real three = 3.0;
212 const real two = 2.0;
214 real dr, dr2, dist, kdist, kdist2, fbond, vbond, fij, vtot;
216 int i, m, ki, type, ai, aj;
220 for (i = 0; (i < nbonds); )
222 type = forceatoms[i++];
223 ai = forceatoms[i++];
224 aj = forceatoms[i++];
226 b0 = forceparams[type].cubic.b0;
227 kb = forceparams[type].cubic.kb;
228 kcub = forceparams[type].cubic.kcub;
230 ki = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /* 3 */
231 dr2 = iprod(dx, dx); /* 5 */
238 dr = dr2*gmx::invsqrt(dr2); /* 10 */
243 vbond = kdist2 + kcub*kdist2*dist;
244 fbond = -(two*kdist + three*kdist2*kcub)/dr;
246 vtot += vbond; /* 21 */
250 ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, aj), dt);
253 for (m = 0; (m < DIM); m++) /* 15 */
258 fshift[ki][m] += fij;
259 fshift[CENTRAL][m] -= fij;
265 real FENE_bonds(int nbonds,
266 const t_iatom forceatoms[], const t_iparams forceparams[],
267 const rvec x[], rvec4 f[], rvec fshift[],
268 const t_pbc *pbc, const t_graph *g,
269 real gmx_unused lambda, real gmx_unused *dvdlambda,
270 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
271 int *global_atom_index)
273 const real half = 0.5;
274 const real one = 1.0;
276 real dr2, bm2, omdr2obm2, fbond, vbond, fij, vtot;
278 int i, m, ki, type, ai, aj;
282 for (i = 0; (i < nbonds); )
284 type = forceatoms[i++];
285 ai = forceatoms[i++];
286 aj = forceatoms[i++];
288 bm = forceparams[type].fene.bm;
289 kb = forceparams[type].fene.kb;
291 ki = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /* 3 */
292 dr2 = iprod(dx, dx); /* 5 */
304 "r^2 (%f) >= bm^2 (%f) in FENE bond between atoms %d and %d",
306 glatnr(global_atom_index, ai),
307 glatnr(global_atom_index, aj));
310 omdr2obm2 = one - dr2/bm2;
312 vbond = -half*kb*bm2*std::log(omdr2obm2);
313 fbond = -kb/omdr2obm2;
315 vtot += vbond; /* 35 */
319 ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, aj), dt);
322 for (m = 0; (m < DIM); m++) /* 15 */
327 fshift[ki][m] += fij;
328 fshift[CENTRAL][m] -= fij;
334 static real harmonic(real kA, real kB, real xA, real xB, real x, real lambda,
337 const real half = 0.5;
338 real L1, kk, x0, dx, dx2;
339 real v, f, dvdlambda;
342 kk = L1*kA+lambda*kB;
343 x0 = L1*xA+lambda*xB;
350 dvdlambda = half*(kB-kA)*dx2 + (xA-xB)*kk*dx;
357 /* That was 19 flops */
361 real bonds(int nbonds,
362 const t_iatom forceatoms[], const t_iparams forceparams[],
363 const rvec x[], rvec4 f[], rvec fshift[],
364 const t_pbc *pbc, const t_graph *g,
365 real lambda, real *dvdlambda,
366 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
367 int gmx_unused *global_atom_index)
369 int i, m, ki, ai, aj, type;
370 real dr, dr2, fbond, vbond, fij, vtot;
375 for (i = 0; (i < nbonds); )
377 type = forceatoms[i++];
378 ai = forceatoms[i++];
379 aj = forceatoms[i++];
381 ki = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /* 3 */
382 dr2 = iprod(dx, dx); /* 5 */
383 dr = std::sqrt(dr2); /* 10 */
385 *dvdlambda += harmonic(forceparams[type].harmonic.krA,
386 forceparams[type].harmonic.krB,
387 forceparams[type].harmonic.rA,
388 forceparams[type].harmonic.rB,
389 dr, lambda, &vbond, &fbond); /* 19 */
397 vtot += vbond; /* 1*/
398 fbond *= gmx::invsqrt(dr2); /* 6 */
401 ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, aj), dt);
404 for (m = 0; (m < DIM); m++) /* 15 */
409 fshift[ki][m] += fij;
410 fshift[CENTRAL][m] -= fij;
416 real restraint_bonds(int nbonds,
417 const t_iatom forceatoms[], const t_iparams forceparams[],
418 const rvec x[], rvec4 f[], rvec fshift[],
419 const t_pbc *pbc, const t_graph *g,
420 real lambda, real *dvdlambda,
421 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
422 int gmx_unused *global_atom_index)
424 int i, m, ki, ai, aj, type;
425 real dr, dr2, fbond, vbond, fij, vtot;
427 real low, dlow, up1, dup1, up2, dup2, k, dk;
435 for (i = 0; (i < nbonds); )
437 type = forceatoms[i++];
438 ai = forceatoms[i++];
439 aj = forceatoms[i++];
441 ki = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /* 3 */
442 dr2 = iprod(dx, dx); /* 5 */
443 dr = dr2*gmx::invsqrt(dr2); /* 10 */
445 low = L1*forceparams[type].restraint.lowA + lambda*forceparams[type].restraint.lowB;
446 dlow = -forceparams[type].restraint.lowA + forceparams[type].restraint.lowB;
447 up1 = L1*forceparams[type].restraint.up1A + lambda*forceparams[type].restraint.up1B;
448 dup1 = -forceparams[type].restraint.up1A + forceparams[type].restraint.up1B;
449 up2 = L1*forceparams[type].restraint.up2A + lambda*forceparams[type].restraint.up2B;
450 dup2 = -forceparams[type].restraint.up2A + forceparams[type].restraint.up2B;
451 k = L1*forceparams[type].restraint.kA + lambda*forceparams[type].restraint.kB;
452 dk = -forceparams[type].restraint.kA + forceparams[type].restraint.kB;
461 *dvdlambda += 0.5*dk*drh2 - k*dlow*drh;
474 *dvdlambda += 0.5*dk*drh2 - k*dup1*drh;
479 vbond = k*(up2 - up1)*(0.5*(up2 - up1) + drh);
480 fbond = -k*(up2 - up1);
481 *dvdlambda += dk*(up2 - up1)*(0.5*(up2 - up1) + drh)
482 + k*(dup2 - dup1)*(up2 - up1 + drh)
483 - k*(up2 - up1)*dup2;
491 vtot += vbond; /* 1*/
492 fbond *= gmx::invsqrt(dr2); /* 6 */
495 ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, aj), dt);
498 for (m = 0; (m < DIM); m++) /* 15 */
503 fshift[ki][m] += fij;
504 fshift[CENTRAL][m] -= fij;
511 real polarize(int nbonds,
512 const t_iatom forceatoms[], const t_iparams forceparams[],
513 const rvec x[], rvec4 f[], rvec fshift[],
514 const t_pbc *pbc, const t_graph *g,
515 real lambda, real *dvdlambda,
516 const t_mdatoms *md, t_fcdata gmx_unused *fcd,
517 int gmx_unused *global_atom_index)
519 int i, m, ki, ai, aj, type;
520 real dr, dr2, fbond, vbond, fij, vtot, ksh;
525 for (i = 0; (i < nbonds); )
527 type = forceatoms[i++];
528 ai = forceatoms[i++];
529 aj = forceatoms[i++];
530 ksh = gmx::square(md->chargeA[aj])*ONE_4PI_EPS0/forceparams[type].polarize.alpha;
532 ki = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /* 3 */
533 dr2 = iprod(dx, dx); /* 5 */
534 dr = std::sqrt(dr2); /* 10 */
536 *dvdlambda += harmonic(ksh, ksh, 0, 0, dr, lambda, &vbond, &fbond); /* 19 */
543 vtot += vbond; /* 1*/
544 fbond *= gmx::invsqrt(dr2); /* 6 */
548 ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, aj), dt);
551 for (m = 0; (m < DIM); m++) /* 15 */
556 fshift[ki][m] += fij;
557 fshift[CENTRAL][m] -= fij;
563 real anharm_polarize(int nbonds,
564 const t_iatom forceatoms[], const t_iparams forceparams[],
565 const rvec x[], rvec4 f[], rvec fshift[],
566 const t_pbc *pbc, const t_graph *g,
567 real lambda, real *dvdlambda,
568 const t_mdatoms *md, t_fcdata gmx_unused *fcd,
569 int gmx_unused *global_atom_index)
571 int i, m, ki, ai, aj, type;
572 real dr, dr2, fbond, vbond, fij, vtot, ksh, khyp, drcut, ddr, ddr3;
577 for (i = 0; (i < nbonds); )
579 type = forceatoms[i++];
580 ai = forceatoms[i++];
581 aj = forceatoms[i++];
582 ksh = gmx::square(md->chargeA[aj])*ONE_4PI_EPS0/forceparams[type].anharm_polarize.alpha; /* 7*/
583 khyp = forceparams[type].anharm_polarize.khyp;
584 drcut = forceparams[type].anharm_polarize.drcut;
586 ki = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /* 3 */
587 dr2 = iprod(dx, dx); /* 5 */
588 dr = dr2*gmx::invsqrt(dr2); /* 10 */
590 *dvdlambda += harmonic(ksh, ksh, 0, 0, dr, lambda, &vbond, &fbond); /* 19 */
601 vbond += khyp*ddr*ddr3;
602 fbond -= 4*khyp*ddr3;
604 fbond *= gmx::invsqrt(dr2); /* 6 */
605 vtot += vbond; /* 1*/
609 ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, aj), dt);
612 for (m = 0; (m < DIM); m++) /* 15 */
617 fshift[ki][m] += fij;
618 fshift[CENTRAL][m] -= fij;
624 real water_pol(int nbonds,
625 const t_iatom forceatoms[], const t_iparams forceparams[],
626 const rvec x[], rvec4 f[], rvec gmx_unused fshift[],
627 const t_pbc gmx_unused *pbc, const t_graph gmx_unused *g,
628 real gmx_unused lambda, real gmx_unused *dvdlambda,
629 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
630 int gmx_unused *global_atom_index)
632 /* This routine implements anisotropic polarizibility for water, through
633 * a shell connected to a dummy with spring constant that differ in the
634 * three spatial dimensions in the molecular frame.
636 int i, m, aO, aH1, aH2, aD, aS, type, type0, ki;
638 rvec dOH1, dOH2, dHH, dOD, dDS, nW, kk, dx, kdx, proj;
639 real vtot, fij, r_HH, r_OD, r_nW, tx, ty, tz, qS;
644 type0 = forceatoms[0];
646 qS = md->chargeA[aS];
647 kk[XX] = gmx::square(qS)*ONE_4PI_EPS0/forceparams[type0].wpol.al_x;
648 kk[YY] = gmx::square(qS)*ONE_4PI_EPS0/forceparams[type0].wpol.al_y;
649 kk[ZZ] = gmx::square(qS)*ONE_4PI_EPS0/forceparams[type0].wpol.al_z;
650 r_HH = 1.0/forceparams[type0].wpol.rHH;
651 for (i = 0; (i < nbonds); i += 6)
653 type = forceatoms[i];
656 gmx_fatal(FARGS, "Sorry, type = %d, type0 = %d, file = %s, line = %d",
657 type, type0, __FILE__, __LINE__);
659 aO = forceatoms[i+1];
660 aH1 = forceatoms[i+2];
661 aH2 = forceatoms[i+3];
662 aD = forceatoms[i+4];
663 aS = forceatoms[i+5];
665 /* Compute vectors describing the water frame */
666 pbc_rvec_sub(pbc, x[aH1], x[aO], dOH1);
667 pbc_rvec_sub(pbc, x[aH2], x[aO], dOH2);
668 pbc_rvec_sub(pbc, x[aH2], x[aH1], dHH);
669 pbc_rvec_sub(pbc, x[aD], x[aO], dOD);
670 ki = pbc_rvec_sub(pbc, x[aS], x[aD], dDS);
671 cprod(dOH1, dOH2, nW);
673 /* Compute inverse length of normal vector
674 * (this one could be precomputed, but I'm too lazy now)
676 r_nW = gmx::invsqrt(iprod(nW, nW));
677 /* This is for precision, but does not make a big difference,
680 r_OD = gmx::invsqrt(iprod(dOD, dOD));
682 /* Normalize the vectors in the water frame */
684 svmul(r_HH, dHH, dHH);
685 svmul(r_OD, dOD, dOD);
687 /* Compute displacement of shell along components of the vector */
688 dx[ZZ] = iprod(dDS, dOD);
689 /* Compute projection on the XY plane: dDS - dx[ZZ]*dOD */
690 for (m = 0; (m < DIM); m++)
692 proj[m] = dDS[m]-dx[ZZ]*dOD[m];
695 /*dx[XX] = iprod(dDS,nW);
696 dx[YY] = iprod(dDS,dHH);*/
697 dx[XX] = iprod(proj, nW);
698 for (m = 0; (m < DIM); m++)
700 proj[m] -= dx[XX]*nW[m];
702 dx[YY] = iprod(proj, dHH);
703 /* Now compute the forces and energy */
704 kdx[XX] = kk[XX]*dx[XX];
705 kdx[YY] = kk[YY]*dx[YY];
706 kdx[ZZ] = kk[ZZ]*dx[ZZ];
707 vtot += iprod(dx, kdx);
711 ivec_sub(SHIFT_IVEC(g, aS), SHIFT_IVEC(g, aD), dt);
715 for (m = 0; (m < DIM); m++)
717 /* This is a tensor operation but written out for speed */
724 fshift[ki][m] += fij;
725 fshift[CENTRAL][m] -= fij;
732 static real do_1_thole(const rvec xi, const rvec xj, rvec fi, rvec fj,
733 const t_pbc *pbc, real qq,
734 rvec fshift[], real afac)
737 real r12sq, r12_1, r12bar, v0, v1, fscal, ebar, fff;
740 t = pbc_rvec_sub(pbc, xi, xj, r12); /* 3 */
742 r12sq = iprod(r12, r12); /* 5 */
743 r12_1 = gmx::invsqrt(r12sq); /* 5 */
744 r12bar = afac/r12_1; /* 5 */
745 v0 = qq*ONE_4PI_EPS0*r12_1; /* 2 */
746 ebar = std::exp(-r12bar); /* 5 */
747 v1 = (1-(1+0.5*r12bar)*ebar); /* 4 */
748 fscal = ((v0*r12_1)*v1 - v0*0.5*afac*ebar*(r12bar+1))*r12_1; /* 9 */
750 for (m = 0; (m < DIM); m++)
756 fshift[CENTRAL][m] -= fff;
759 return v0*v1; /* 1 */
763 real thole_pol(int nbonds,
764 const t_iatom forceatoms[], const t_iparams forceparams[],
765 const rvec x[], rvec4 f[], rvec fshift[],
766 const t_pbc *pbc, const t_graph gmx_unused *g,
767 real gmx_unused lambda, real gmx_unused *dvdlambda,
768 const t_mdatoms *md, t_fcdata gmx_unused *fcd,
769 int gmx_unused *global_atom_index)
771 /* Interaction between two pairs of particles with opposite charge */
772 int i, type, a1, da1, a2, da2;
773 real q1, q2, qq, a, al1, al2, afac;
776 for (i = 0; (i < nbonds); )
778 type = forceatoms[i++];
779 a1 = forceatoms[i++];
780 da1 = forceatoms[i++];
781 a2 = forceatoms[i++];
782 da2 = forceatoms[i++];
783 q1 = md->chargeA[da1];
784 q2 = md->chargeA[da2];
785 a = forceparams[type].thole.a;
786 al1 = forceparams[type].thole.alpha1;
787 al2 = forceparams[type].thole.alpha2;
789 afac = a*gmx::invsixthroot(al1*al2);
790 V += do_1_thole(x[a1], x[a2], f[a1], f[a2], pbc, qq, fshift, afac);
791 V += do_1_thole(x[da1], x[a2], f[da1], f[a2], pbc, -qq, fshift, afac);
792 V += do_1_thole(x[a1], x[da2], f[a1], f[da2], pbc, -qq, fshift, afac);
793 V += do_1_thole(x[da1], x[da2], f[da1], f[da2], pbc, qq, fshift, afac);
799 real bond_angle(const rvec xi, const rvec xj, const rvec xk, const t_pbc *pbc,
800 rvec r_ij, rvec r_kj, real *costh,
802 /* Return value is the angle between the bonds i-j and j-k */
807 *t1 = pbc_rvec_sub(pbc, xi, xj, r_ij); /* 3 */
808 *t2 = pbc_rvec_sub(pbc, xk, xj, r_kj); /* 3 */
810 *costh = cos_angle(r_ij, r_kj); /* 25 */
811 th = std::acos(*costh); /* 10 */
816 real angles(int nbonds,
817 const t_iatom forceatoms[], const t_iparams forceparams[],
818 const rvec x[], rvec4 f[], rvec fshift[],
819 const t_pbc *pbc, const t_graph *g,
820 real lambda, real *dvdlambda,
821 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
822 int gmx_unused *global_atom_index)
824 int i, ai, aj, ak, t1, t2, type;
826 real cos_theta, cos_theta2, theta, dVdt, va, vtot;
827 ivec jt, dt_ij, dt_kj;
830 for (i = 0; i < nbonds; )
832 type = forceatoms[i++];
833 ai = forceatoms[i++];
834 aj = forceatoms[i++];
835 ak = forceatoms[i++];
837 theta = bond_angle(x[ai], x[aj], x[ak], pbc,
838 r_ij, r_kj, &cos_theta, &t1, &t2); /* 41 */
840 *dvdlambda += harmonic(forceparams[type].harmonic.krA,
841 forceparams[type].harmonic.krB,
842 forceparams[type].harmonic.rA*DEG2RAD,
843 forceparams[type].harmonic.rB*DEG2RAD,
844 theta, lambda, &va, &dVdt); /* 21 */
847 cos_theta2 = gmx::square(cos_theta);
857 st = dVdt*gmx::invsqrt(1 - cos_theta2); /* 12 */
858 sth = st*cos_theta; /* 1 */
859 nrij2 = iprod(r_ij, r_ij); /* 5 */
860 nrkj2 = iprod(r_kj, r_kj); /* 5 */
862 nrij_1 = gmx::invsqrt(nrij2); /* 10 */
863 nrkj_1 = gmx::invsqrt(nrkj2); /* 10 */
865 cik = st*nrij_1*nrkj_1; /* 2 */
866 cii = sth*nrij_1*nrij_1; /* 2 */
867 ckk = sth*nrkj_1*nrkj_1; /* 2 */
869 for (m = 0; m < DIM; m++)
871 f_i[m] = -(cik*r_kj[m] - cii*r_ij[m]);
872 f_k[m] = -(cik*r_ij[m] - ckk*r_kj[m]);
873 f_j[m] = -f_i[m] - f_k[m];
880 copy_ivec(SHIFT_IVEC(g, aj), jt);
882 ivec_sub(SHIFT_IVEC(g, ai), jt, dt_ij);
883 ivec_sub(SHIFT_IVEC(g, ak), jt, dt_kj);
887 rvec_inc(fshift[t1], f_i);
888 rvec_inc(fshift[CENTRAL], f_j);
889 rvec_inc(fshift[t2], f_k);
896 #if GMX_SIMD_HAVE_REAL
898 /* As angles, but using SIMD to calculate many angles at once.
899 * This routines does not calculate energies and shift forces.
902 angles_noener_simd(int nbonds,
903 const t_iatom forceatoms[], const t_iparams forceparams[],
904 const rvec x[], rvec4 f[],
905 const t_pbc *pbc, const t_graph gmx_unused *g,
906 real gmx_unused lambda,
907 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
908 int gmx_unused *global_atom_index)
913 alignas(GMX_SIMD_ALIGNMENT) std::int32_t ai[GMX_SIMD_REAL_WIDTH];
914 alignas(GMX_SIMD_ALIGNMENT) std::int32_t aj[GMX_SIMD_REAL_WIDTH];
915 alignas(GMX_SIMD_ALIGNMENT) std::int32_t ak[GMX_SIMD_REAL_WIDTH];
916 alignas(GMX_SIMD_ALIGNMENT) real coeff[2*GMX_SIMD_REAL_WIDTH];
917 SimdReal deg2rad_S(DEG2RAD);
918 SimdReal xi_S, yi_S, zi_S;
919 SimdReal xj_S, yj_S, zj_S;
920 SimdReal xk_S, yk_S, zk_S;
921 SimdReal k_S, theta0_S;
922 SimdReal rijx_S, rijy_S, rijz_S;
923 SimdReal rkjx_S, rkjy_S, rkjz_S;
925 SimdReal min_one_plus_eps_S(-1.0 + 2.0*GMX_REAL_EPS); // Smallest number > -1
928 SimdReal nrij2_S, nrij_1_S;
929 SimdReal nrkj2_S, nrkj_1_S;
930 SimdReal cos_S, invsin_S;
932 SimdReal st_S, sth_S;
933 SimdReal cik_S, cii_S, ckk_S;
934 SimdReal f_ix_S, f_iy_S, f_iz_S;
935 SimdReal f_kx_S, f_ky_S, f_kz_S;
936 alignas(GMX_SIMD_ALIGNMENT) real pbc_simd[9*GMX_SIMD_REAL_WIDTH];
938 set_pbc_simd(pbc, pbc_simd);
940 /* nbonds is the number of angles times nfa1, here we step GMX_SIMD_REAL_WIDTH angles */
941 for (i = 0; (i < nbonds); i += GMX_SIMD_REAL_WIDTH*nfa1)
943 /* Collect atoms for GMX_SIMD_REAL_WIDTH angles.
944 * iu indexes into forceatoms, we should not let iu go beyond nbonds.
947 for (s = 0; s < GMX_SIMD_REAL_WIDTH; s++)
949 type = forceatoms[iu];
950 ai[s] = forceatoms[iu+1];
951 aj[s] = forceatoms[iu+2];
952 ak[s] = forceatoms[iu+3];
954 /* At the end fill the arrays with the last atoms and 0 params */
955 if (i + s*nfa1 < nbonds)
957 coeff[s] = forceparams[type].harmonic.krA;
958 coeff[GMX_SIMD_REAL_WIDTH+s] = forceparams[type].harmonic.rA;
960 if (iu + nfa1 < nbonds)
968 coeff[GMX_SIMD_REAL_WIDTH+s] = 0;
972 /* Store the non PBC corrected distances packed and aligned */
973 gatherLoadUTranspose<3>(reinterpret_cast<const real *>(x), ai, &xi_S, &yi_S, &zi_S);
974 gatherLoadUTranspose<3>(reinterpret_cast<const real *>(x), aj, &xj_S, &yj_S, &zj_S);
975 gatherLoadUTranspose<3>(reinterpret_cast<const real *>(x), ak, &xk_S, &yk_S, &zk_S);
976 rijx_S = xi_S - xj_S;
977 rijy_S = yi_S - yj_S;
978 rijz_S = zi_S - zj_S;
979 rkjx_S = xk_S - xj_S;
980 rkjy_S = yk_S - yj_S;
981 rkjz_S = zk_S - zj_S;
983 k_S = load<SimdReal>(coeff);
984 theta0_S = load<SimdReal>(coeff+GMX_SIMD_REAL_WIDTH) * deg2rad_S;
986 pbc_correct_dx_simd(&rijx_S, &rijy_S, &rijz_S, pbc_simd);
987 pbc_correct_dx_simd(&rkjx_S, &rkjy_S, &rkjz_S, pbc_simd);
989 rij_rkj_S = iprod(rijx_S, rijy_S, rijz_S,
990 rkjx_S, rkjy_S, rkjz_S);
992 nrij2_S = norm2(rijx_S, rijy_S, rijz_S);
993 nrkj2_S = norm2(rkjx_S, rkjy_S, rkjz_S);
995 nrij_1_S = invsqrt(nrij2_S);
996 nrkj_1_S = invsqrt(nrkj2_S);
998 cos_S = rij_rkj_S * nrij_1_S * nrkj_1_S;
1000 /* To allow for 180 degrees, we take the max of cos and -1 + 1bit,
1001 * so we can safely get the 1/sin from 1/sqrt(1 - cos^2).
1002 * This also ensures that rounding errors would cause the argument
1003 * of simdAcos to be < -1.
1004 * Note that we do not take precautions for cos(0)=1, so the outer
1005 * atoms in an angle should not be on top of each other.
1007 cos_S = max(cos_S, min_one_plus_eps_S);
1009 theta_S = acos(cos_S);
1011 invsin_S = invsqrt( one_S - cos_S * cos_S );
1013 st_S = k_S * (theta0_S - theta_S) * invsin_S;
1014 sth_S = st_S * cos_S;
1016 cik_S = st_S * nrij_1_S * nrkj_1_S;
1017 cii_S = sth_S * nrij_1_S * nrij_1_S;
1018 ckk_S = sth_S * nrkj_1_S * nrkj_1_S;
1020 f_ix_S = cii_S * rijx_S;
1021 f_ix_S = fnma(cik_S, rkjx_S, f_ix_S);
1022 f_iy_S = cii_S * rijy_S;
1023 f_iy_S = fnma(cik_S, rkjy_S, f_iy_S);
1024 f_iz_S = cii_S * rijz_S;
1025 f_iz_S = fnma(cik_S, rkjz_S, f_iz_S);
1026 f_kx_S = ckk_S * rkjx_S;
1027 f_kx_S = fnma(cik_S, rijx_S, f_kx_S);
1028 f_ky_S = ckk_S * rkjy_S;
1029 f_ky_S = fnma(cik_S, rijy_S, f_ky_S);
1030 f_kz_S = ckk_S * rkjz_S;
1031 f_kz_S = fnma(cik_S, rijz_S, f_kz_S);
1033 transposeScatterIncrU<4>(reinterpret_cast<real *>(f), ai, f_ix_S, f_iy_S, f_iz_S);
1034 transposeScatterDecrU<4>(reinterpret_cast<real *>(f), aj, f_ix_S + f_kx_S, f_iy_S + f_ky_S, f_iz_S + f_kz_S);
1035 transposeScatterIncrU<4>(reinterpret_cast<real *>(f), ak, f_kx_S, f_ky_S, f_kz_S);
1039 #endif // GMX_SIMD_HAVE_REAL
1041 real linear_angles(int nbonds,
1042 const t_iatom forceatoms[], const t_iparams forceparams[],
1043 const rvec x[], rvec4 f[], rvec fshift[],
1044 const t_pbc *pbc, const t_graph *g,
1045 real lambda, real *dvdlambda,
1046 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
1047 int gmx_unused *global_atom_index)
1049 int i, m, ai, aj, ak, t1, t2, type;
1051 real L1, kA, kB, aA, aB, dr, dr2, va, vtot, a, b, klin;
1052 ivec jt, dt_ij, dt_kj;
1053 rvec r_ij, r_kj, r_ik, dx;
1057 for (i = 0; (i < nbonds); )
1059 type = forceatoms[i++];
1060 ai = forceatoms[i++];
1061 aj = forceatoms[i++];
1062 ak = forceatoms[i++];
1064 kA = forceparams[type].linangle.klinA;
1065 kB = forceparams[type].linangle.klinB;
1066 klin = L1*kA + lambda*kB;
1068 aA = forceparams[type].linangle.aA;
1069 aB = forceparams[type].linangle.aB;
1070 a = L1*aA+lambda*aB;
1073 t1 = pbc_rvec_sub(pbc, x[ai], x[aj], r_ij);
1074 t2 = pbc_rvec_sub(pbc, x[ak], x[aj], r_kj);
1075 rvec_sub(r_ij, r_kj, r_ik);
1078 for (m = 0; (m < DIM); m++)
1080 dr = -a * r_ij[m] - b * r_kj[m];
1085 f_j[m] = -(f_i[m]+f_k[m]);
1091 *dvdlambda += 0.5*(kB-kA)*dr2 + klin*(aB-aA)*iprod(dx, r_ik);
1097 copy_ivec(SHIFT_IVEC(g, aj), jt);
1099 ivec_sub(SHIFT_IVEC(g, ai), jt, dt_ij);
1100 ivec_sub(SHIFT_IVEC(g, ak), jt, dt_kj);
1101 t1 = IVEC2IS(dt_ij);
1102 t2 = IVEC2IS(dt_kj);
1104 rvec_inc(fshift[t1], f_i);
1105 rvec_inc(fshift[CENTRAL], f_j);
1106 rvec_inc(fshift[t2], f_k);
1111 real urey_bradley(int nbonds,
1112 const t_iatom forceatoms[], const t_iparams forceparams[],
1113 const rvec x[], rvec4 f[], rvec fshift[],
1114 const t_pbc *pbc, const t_graph *g,
1115 real lambda, real *dvdlambda,
1116 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
1117 int gmx_unused *global_atom_index)
1119 int i, m, ai, aj, ak, t1, t2, type, ki;
1120 rvec r_ij, r_kj, r_ik;
1121 real cos_theta, cos_theta2, theta;
1122 real dVdt, va, vtot, dr, dr2, vbond, fbond, fik;
1123 real kthA, th0A, kUBA, r13A, kthB, th0B, kUBB, r13B;
1124 ivec jt, dt_ij, dt_kj, dt_ik;
1127 for (i = 0; (i < nbonds); )
1129 type = forceatoms[i++];
1130 ai = forceatoms[i++];
1131 aj = forceatoms[i++];
1132 ak = forceatoms[i++];
1133 th0A = forceparams[type].u_b.thetaA*DEG2RAD;
1134 kthA = forceparams[type].u_b.kthetaA;
1135 r13A = forceparams[type].u_b.r13A;
1136 kUBA = forceparams[type].u_b.kUBA;
1137 th0B = forceparams[type].u_b.thetaB*DEG2RAD;
1138 kthB = forceparams[type].u_b.kthetaB;
1139 r13B = forceparams[type].u_b.r13B;
1140 kUBB = forceparams[type].u_b.kUBB;
1142 theta = bond_angle(x[ai], x[aj], x[ak], pbc,
1143 r_ij, r_kj, &cos_theta, &t1, &t2); /* 41 */
1145 *dvdlambda += harmonic(kthA, kthB, th0A, th0B, theta, lambda, &va, &dVdt); /* 21 */
1148 ki = pbc_rvec_sub(pbc, x[ai], x[ak], r_ik); /* 3 */
1149 dr2 = iprod(r_ik, r_ik); /* 5 */
1150 dr = dr2*gmx::invsqrt(dr2); /* 10 */
1152 *dvdlambda += harmonic(kUBA, kUBB, r13A, r13B, dr, lambda, &vbond, &fbond); /* 19 */
1154 cos_theta2 = gmx::square(cos_theta); /* 1 */
1162 st = dVdt*gmx::invsqrt(1 - cos_theta2); /* 12 */
1163 sth = st*cos_theta; /* 1 */
1164 nrkj2 = iprod(r_kj, r_kj); /* 5 */
1165 nrij2 = iprod(r_ij, r_ij);
1167 cik = st*gmx::invsqrt(nrkj2*nrij2); /* 12 */
1168 cii = sth/nrij2; /* 10 */
1169 ckk = sth/nrkj2; /* 10 */
1171 for (m = 0; (m < DIM); m++) /* 39 */
1173 f_i[m] = -(cik*r_kj[m]-cii*r_ij[m]);
1174 f_k[m] = -(cik*r_ij[m]-ckk*r_kj[m]);
1175 f_j[m] = -f_i[m]-f_k[m];
1182 copy_ivec(SHIFT_IVEC(g, aj), jt);
1184 ivec_sub(SHIFT_IVEC(g, ai), jt, dt_ij);
1185 ivec_sub(SHIFT_IVEC(g, ak), jt, dt_kj);
1186 t1 = IVEC2IS(dt_ij);
1187 t2 = IVEC2IS(dt_kj);
1189 rvec_inc(fshift[t1], f_i);
1190 rvec_inc(fshift[CENTRAL], f_j);
1191 rvec_inc(fshift[t2], f_k);
1193 /* Time for the bond calculations */
1199 vtot += vbond; /* 1*/
1200 fbond *= gmx::invsqrt(dr2); /* 6 */
1204 ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, ak), dt_ik);
1205 ki = IVEC2IS(dt_ik);
1207 for (m = 0; (m < DIM); m++) /* 15 */
1209 fik = fbond*r_ik[m];
1212 fshift[ki][m] += fik;
1213 fshift[CENTRAL][m] -= fik;
1219 #if GMX_SIMD_HAVE_REAL
1221 /* As urey_bradley, but using SIMD to calculate many potentials at once.
1222 * This routines does not calculate energies and shift forces.
1224 void urey_bradley_noener_simd(int nbonds,
1225 const t_iatom forceatoms[], const t_iparams forceparams[],
1226 const rvec x[], rvec4 f[],
1227 const t_pbc *pbc, const t_graph gmx_unused *g,
1228 real gmx_unused lambda,
1229 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
1230 int gmx_unused *global_atom_index)
1232 constexpr int nfa1 = 4;
1233 alignas(GMX_SIMD_ALIGNMENT) std::int32_t ai[GMX_SIMD_REAL_WIDTH];
1234 alignas(GMX_SIMD_ALIGNMENT) std::int32_t aj[GMX_SIMD_REAL_WIDTH];
1235 alignas(GMX_SIMD_ALIGNMENT) std::int32_t ak[GMX_SIMD_REAL_WIDTH];
1236 alignas(GMX_SIMD_ALIGNMENT) real coeff[4*GMX_SIMD_REAL_WIDTH];
1237 alignas(GMX_SIMD_ALIGNMENT) real pbc_simd[9*GMX_SIMD_REAL_WIDTH];
1239 set_pbc_simd(pbc, pbc_simd);
1241 /* nbonds is the number of angles times nfa1, here we step GMX_SIMD_REAL_WIDTH angles */
1242 for (int i = 0; i < nbonds; i += GMX_SIMD_REAL_WIDTH*nfa1)
1244 /* Collect atoms for GMX_SIMD_REAL_WIDTH angles.
1245 * iu indexes into forceatoms, we should not let iu go beyond nbonds.
1248 for (int s = 0; s < GMX_SIMD_REAL_WIDTH; s++)
1250 const int type = forceatoms[iu];
1251 ai[s] = forceatoms[iu+1];
1252 aj[s] = forceatoms[iu+2];
1253 ak[s] = forceatoms[iu+3];
1255 /* At the end fill the arrays with the last atoms and 0 params */
1256 if (i + s*nfa1 < nbonds)
1258 coeff[s] = forceparams[type].u_b.kthetaA;
1259 coeff[GMX_SIMD_REAL_WIDTH+s] = forceparams[type].u_b.thetaA;
1260 coeff[GMX_SIMD_REAL_WIDTH*2+s] = forceparams[type].u_b.kUBA;
1261 coeff[GMX_SIMD_REAL_WIDTH*3+s] = forceparams[type].u_b.r13A;
1263 if (iu + nfa1 < nbonds)
1271 coeff[GMX_SIMD_REAL_WIDTH+s] = 0;
1272 coeff[GMX_SIMD_REAL_WIDTH*2+s] = 0;
1273 coeff[GMX_SIMD_REAL_WIDTH*3+s] = 0;
1277 SimdReal xi_S, yi_S, zi_S;
1278 SimdReal xj_S, yj_S, zj_S;
1279 SimdReal xk_S, yk_S, zk_S;
1281 /* Store the non PBC corrected distances packed and aligned */
1282 gatherLoadUTranspose<3>(reinterpret_cast<const real *>(x), ai, &xi_S, &yi_S, &zi_S);
1283 gatherLoadUTranspose<3>(reinterpret_cast<const real *>(x), aj, &xj_S, &yj_S, &zj_S);
1284 gatherLoadUTranspose<3>(reinterpret_cast<const real *>(x), ak, &xk_S, &yk_S, &zk_S);
1285 SimdReal rijx_S = xi_S - xj_S;
1286 SimdReal rijy_S = yi_S - yj_S;
1287 SimdReal rijz_S = zi_S - zj_S;
1288 SimdReal rkjx_S = xk_S - xj_S;
1289 SimdReal rkjy_S = yk_S - yj_S;
1290 SimdReal rkjz_S = zk_S - zj_S;
1291 SimdReal rikx_S = xi_S - xk_S;
1292 SimdReal riky_S = yi_S - yk_S;
1293 SimdReal rikz_S = zi_S - zk_S;
1295 const SimdReal ktheta_S = load<SimdReal>(coeff);
1296 const SimdReal theta0_S = load<SimdReal>(coeff+GMX_SIMD_REAL_WIDTH) * DEG2RAD;
1297 const SimdReal kUB_S = load<SimdReal>(coeff+2*GMX_SIMD_REAL_WIDTH);
1298 const SimdReal r13_S = load<SimdReal>(coeff+3*GMX_SIMD_REAL_WIDTH);
1300 pbc_correct_dx_simd(&rijx_S, &rijy_S, &rijz_S, pbc_simd);
1301 pbc_correct_dx_simd(&rkjx_S, &rkjy_S, &rkjz_S, pbc_simd);
1302 pbc_correct_dx_simd(&rikx_S, &riky_S, &rikz_S, pbc_simd);
1304 const SimdReal rij_rkj_S = iprod(rijx_S, rijy_S, rijz_S,
1305 rkjx_S, rkjy_S, rkjz_S);
1307 const SimdReal dr2_S = iprod(rikx_S, riky_S, rikz_S,
1308 rikx_S, riky_S, rikz_S);
1310 const SimdReal nrij2_S = norm2(rijx_S, rijy_S, rijz_S);
1311 const SimdReal nrkj2_S = norm2(rkjx_S, rkjy_S, rkjz_S);
1313 const SimdReal nrij_1_S = invsqrt(nrij2_S);
1314 const SimdReal nrkj_1_S = invsqrt(nrkj2_S);
1315 const SimdReal invdr2_S = invsqrt(dr2_S);
1316 const SimdReal dr_S = dr2_S*invdr2_S;
1318 constexpr real min_one_plus_eps = -1.0 + 2.0*GMX_REAL_EPS; // Smallest number > -1
1320 /* To allow for 180 degrees, we take the max of cos and -1 + 1bit,
1321 * so we can safely get the 1/sin from 1/sqrt(1 - cos^2).
1322 * This also ensures that rounding errors would cause the argument
1323 * of simdAcos to be < -1.
1324 * Note that we do not take precautions for cos(0)=1, so the bonds
1325 * in an angle should not align at an angle of 0 degrees.
1327 const SimdReal cos_S = max(rij_rkj_S * nrij_1_S * nrkj_1_S, min_one_plus_eps);
1329 const SimdReal theta_S = acos(cos_S);
1330 const SimdReal invsin_S = invsqrt( 1.0 - cos_S * cos_S );
1331 const SimdReal st_S = ktheta_S * (theta0_S - theta_S) * invsin_S;
1332 const SimdReal sth_S = st_S * cos_S;
1334 const SimdReal cik_S = st_S * nrij_1_S * nrkj_1_S;
1335 const SimdReal cii_S = sth_S * nrij_1_S * nrij_1_S;
1336 const SimdReal ckk_S = sth_S * nrkj_1_S * nrkj_1_S;
1338 const SimdReal sUB_S = kUB_S * (r13_S - dr_S) * invdr2_S;
1340 const SimdReal f_ikx_S = sUB_S * rikx_S;
1341 const SimdReal f_iky_S = sUB_S * riky_S;
1342 const SimdReal f_ikz_S = sUB_S * rikz_S;
1344 const SimdReal f_ix_S = fnma(cik_S, rkjx_S, cii_S * rijx_S) + f_ikx_S;
1345 const SimdReal f_iy_S = fnma(cik_S, rkjy_S, cii_S * rijy_S) + f_iky_S;
1346 const SimdReal f_iz_S = fnma(cik_S, rkjz_S, cii_S * rijz_S) + f_ikz_S;
1347 const SimdReal f_kx_S = fnma(cik_S, rijx_S, ckk_S * rkjx_S) - f_ikx_S;
1348 const SimdReal f_ky_S = fnma(cik_S, rijy_S, ckk_S * rkjy_S) - f_iky_S;
1349 const SimdReal f_kz_S = fnma(cik_S, rijz_S, ckk_S * rkjz_S) - f_ikz_S;
1351 transposeScatterIncrU<4>(reinterpret_cast<real *>(f), ai, f_ix_S, f_iy_S, f_iz_S);
1352 transposeScatterDecrU<4>(reinterpret_cast<real *>(f), aj, f_ix_S + f_kx_S, f_iy_S + f_ky_S, f_iz_S + f_kz_S);
1353 transposeScatterIncrU<4>(reinterpret_cast<real *>(f), ak, f_kx_S, f_ky_S, f_kz_S);
1357 #endif // GMX_SIMD_HAVE_REAL
1359 real quartic_angles(int nbonds,
1360 const t_iatom forceatoms[], const t_iparams forceparams[],
1361 const rvec x[], rvec4 f[], rvec fshift[],
1362 const t_pbc *pbc, const t_graph *g,
1363 real gmx_unused lambda, real gmx_unused *dvdlambda,
1364 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
1365 int gmx_unused *global_atom_index)
1367 int i, j, ai, aj, ak, t1, t2, type;
1369 real cos_theta, cos_theta2, theta, dt, dVdt, va, dtp, c, vtot;
1370 ivec jt, dt_ij, dt_kj;
1373 for (i = 0; (i < nbonds); )
1375 type = forceatoms[i++];
1376 ai = forceatoms[i++];
1377 aj = forceatoms[i++];
1378 ak = forceatoms[i++];
1380 theta = bond_angle(x[ai], x[aj], x[ak], pbc,
1381 r_ij, r_kj, &cos_theta, &t1, &t2); /* 41 */
1383 dt = theta - forceparams[type].qangle.theta*DEG2RAD; /* 2 */
1386 va = forceparams[type].qangle.c[0];
1388 for (j = 1; j <= 4; j++)
1390 c = forceparams[type].qangle.c[j];
1399 cos_theta2 = gmx::square(cos_theta); /* 1 */
1408 st = dVdt*gmx::invsqrt(1 - cos_theta2); /* 12 */
1409 sth = st*cos_theta; /* 1 */
1410 nrkj2 = iprod(r_kj, r_kj); /* 5 */
1411 nrij2 = iprod(r_ij, r_ij);
1413 cik = st*gmx::invsqrt(nrkj2*nrij2); /* 12 */
1414 cii = sth/nrij2; /* 10 */
1415 ckk = sth/nrkj2; /* 10 */
1417 for (m = 0; (m < DIM); m++) /* 39 */
1419 f_i[m] = -(cik*r_kj[m]-cii*r_ij[m]);
1420 f_k[m] = -(cik*r_ij[m]-ckk*r_kj[m]);
1421 f_j[m] = -f_i[m]-f_k[m];
1428 copy_ivec(SHIFT_IVEC(g, aj), jt);
1430 ivec_sub(SHIFT_IVEC(g, ai), jt, dt_ij);
1431 ivec_sub(SHIFT_IVEC(g, ak), jt, dt_kj);
1432 t1 = IVEC2IS(dt_ij);
1433 t2 = IVEC2IS(dt_kj);
1435 rvec_inc(fshift[t1], f_i);
1436 rvec_inc(fshift[CENTRAL], f_j);
1437 rvec_inc(fshift[t2], f_k);
1443 real dih_angle(const rvec xi, const rvec xj, const rvec xk, const rvec xl,
1445 rvec r_ij, rvec r_kj, rvec r_kl, rvec m, rvec n,
1446 int *t1, int *t2, int *t3)
1448 *t1 = pbc_rvec_sub(pbc, xi, xj, r_ij); /* 3 */
1449 *t2 = pbc_rvec_sub(pbc, xk, xj, r_kj); /* 3 */
1450 *t3 = pbc_rvec_sub(pbc, xk, xl, r_kl); /* 3 */
1452 cprod(r_ij, r_kj, m); /* 9 */
1453 cprod(r_kj, r_kl, n); /* 9 */
1454 real phi = gmx_angle(m, n); /* 49 (assuming 25 for atan2) */
1455 real ipr = iprod(r_ij, n); /* 5 */
1456 real sign = (ipr < 0.0) ? -1.0 : 1.0;
1457 phi = sign*phi; /* 1 */
1463 #if GMX_SIMD_HAVE_REAL
1465 /* As dih_angle above, but calculates 4 dihedral angles at once using SIMD,
1466 * also calculates the pre-factor required for the dihedral force update.
1467 * Note that bv and buf should be register aligned.
1470 dih_angle_simd(const rvec *x,
1471 const int *ai, const int *aj, const int *ak, const int *al,
1472 const real *pbc_simd,
1474 SimdReal *mx_S, SimdReal *my_S, SimdReal *mz_S,
1475 SimdReal *nx_S, SimdReal *ny_S, SimdReal *nz_S,
1476 SimdReal *nrkj_m2_S,
1477 SimdReal *nrkj_n2_S,
1481 SimdReal xi_S, yi_S, zi_S;
1482 SimdReal xj_S, yj_S, zj_S;
1483 SimdReal xk_S, yk_S, zk_S;
1484 SimdReal xl_S, yl_S, zl_S;
1485 SimdReal rijx_S, rijy_S, rijz_S;
1486 SimdReal rkjx_S, rkjy_S, rkjz_S;
1487 SimdReal rklx_S, rkly_S, rklz_S;
1488 SimdReal cx_S, cy_S, cz_S;
1492 SimdReal iprm_S, iprn_S;
1493 SimdReal nrkj2_S, nrkj_1_S, nrkj_2_S, nrkj_S;
1495 SimdReal nrkj2_min_S;
1496 SimdReal real_eps_S;
1498 /* Used to avoid division by zero.
1499 * We take into acount that we multiply the result by real_eps_S.
1501 nrkj2_min_S = SimdReal(GMX_REAL_MIN/(2*GMX_REAL_EPS));
1503 /* The value of the last significant bit (GMX_REAL_EPS is half of that) */
1504 real_eps_S = SimdReal(2*GMX_REAL_EPS);
1506 /* Store the non PBC corrected distances packed and aligned */
1507 gatherLoadUTranspose<3>(reinterpret_cast<const real *>(x), ai, &xi_S, &yi_S, &zi_S);
1508 gatherLoadUTranspose<3>(reinterpret_cast<const real *>(x), aj, &xj_S, &yj_S, &zj_S);
1509 gatherLoadUTranspose<3>(reinterpret_cast<const real *>(x), ak, &xk_S, &yk_S, &zk_S);
1510 gatherLoadUTranspose<3>(reinterpret_cast<const real *>(x), al, &xl_S, &yl_S, &zl_S);
1511 rijx_S = xi_S - xj_S;
1512 rijy_S = yi_S - yj_S;
1513 rijz_S = zi_S - zj_S;
1514 rkjx_S = xk_S - xj_S;
1515 rkjy_S = yk_S - yj_S;
1516 rkjz_S = zk_S - zj_S;
1517 rklx_S = xk_S - xl_S;
1518 rkly_S = yk_S - yl_S;
1519 rklz_S = zk_S - zl_S;
1521 pbc_correct_dx_simd(&rijx_S, &rijy_S, &rijz_S, pbc_simd);
1522 pbc_correct_dx_simd(&rkjx_S, &rkjy_S, &rkjz_S, pbc_simd);
1523 pbc_correct_dx_simd(&rklx_S, &rkly_S, &rklz_S, pbc_simd);
1525 cprod(rijx_S, rijy_S, rijz_S,
1526 rkjx_S, rkjy_S, rkjz_S,
1529 cprod(rkjx_S, rkjy_S, rkjz_S,
1530 rklx_S, rkly_S, rklz_S,
1533 cprod(*mx_S, *my_S, *mz_S,
1534 *nx_S, *ny_S, *nz_S,
1535 &cx_S, &cy_S, &cz_S);
1537 cn_S = sqrt(norm2(cx_S, cy_S, cz_S));
1539 s_S = iprod(*mx_S, *my_S, *mz_S, *nx_S, *ny_S, *nz_S);
1541 /* Determine the dihedral angle, the sign might need correction */
1542 *phi_S = atan2(cn_S, s_S);
1544 ipr_S = iprod(rijx_S, rijy_S, rijz_S,
1545 *nx_S, *ny_S, *nz_S);
1547 iprm_S = norm2(*mx_S, *my_S, *mz_S);
1548 iprn_S = norm2(*nx_S, *ny_S, *nz_S);
1550 nrkj2_S = norm2(rkjx_S, rkjy_S, rkjz_S);
1552 /* Avoid division by zero. When zero, the result is multiplied by 0
1553 * anyhow, so the 3 max below do not affect the final result.
1555 nrkj2_S = max(nrkj2_S, nrkj2_min_S);
1556 nrkj_1_S = invsqrt(nrkj2_S);
1557 nrkj_2_S = nrkj_1_S * nrkj_1_S;
1558 nrkj_S = nrkj2_S * nrkj_1_S;
1560 toler_S = nrkj2_S * real_eps_S;
1562 /* Here the plain-C code uses a conditional, but we can't do that in SIMD.
1563 * So we take a max with the tolerance instead. Since we multiply with
1564 * m or n later, the max does not affect the results.
1566 iprm_S = max(iprm_S, toler_S);
1567 iprn_S = max(iprn_S, toler_S);
1568 *nrkj_m2_S = nrkj_S * inv(iprm_S);
1569 *nrkj_n2_S = nrkj_S * inv(iprn_S);
1571 /* Set sign of phi_S with the sign of ipr_S; phi_S is currently positive */
1572 *phi_S = copysign(*phi_S, ipr_S);
1573 *p_S = iprod(rijx_S, rijy_S, rijz_S, rkjx_S, rkjy_S, rkjz_S);
1574 *p_S = *p_S * nrkj_2_S;
1576 *q_S = iprod(rklx_S, rkly_S, rklz_S, rkjx_S, rkjy_S, rkjz_S);
1577 *q_S = *q_S * nrkj_2_S;
1580 #endif // GMX_SIMD_HAVE_REAL
1582 void do_dih_fup(int i, int j, int k, int l, real ddphi,
1583 rvec r_ij, rvec r_kj, rvec r_kl,
1584 rvec m, rvec n, rvec4 f[], rvec fshift[],
1585 const t_pbc *pbc, const t_graph *g,
1586 const rvec x[], int t1, int t2, int t3)
1589 rvec f_i, f_j, f_k, f_l;
1590 rvec uvec, vvec, svec, dx_jl;
1591 real iprm, iprn, nrkj, nrkj2, nrkj_1, nrkj_2;
1592 real a, b, p, q, toler;
1593 ivec jt, dt_ij, dt_kj, dt_lj;
1595 iprm = iprod(m, m); /* 5 */
1596 iprn = iprod(n, n); /* 5 */
1597 nrkj2 = iprod(r_kj, r_kj); /* 5 */
1598 toler = nrkj2*GMX_REAL_EPS;
1599 if ((iprm > toler) && (iprn > toler))
1601 nrkj_1 = gmx::invsqrt(nrkj2); /* 10 */
1602 nrkj_2 = nrkj_1*nrkj_1; /* 1 */
1603 nrkj = nrkj2*nrkj_1; /* 1 */
1604 a = -ddphi*nrkj/iprm; /* 11 */
1605 svmul(a, m, f_i); /* 3 */
1606 b = ddphi*nrkj/iprn; /* 11 */
1607 svmul(b, n, f_l); /* 3 */
1608 p = iprod(r_ij, r_kj); /* 5 */
1609 p *= nrkj_2; /* 1 */
1610 q = iprod(r_kl, r_kj); /* 5 */
1611 q *= nrkj_2; /* 1 */
1612 svmul(p, f_i, uvec); /* 3 */
1613 svmul(q, f_l, vvec); /* 3 */
1614 rvec_sub(uvec, vvec, svec); /* 3 */
1615 rvec_sub(f_i, svec, f_j); /* 3 */
1616 rvec_add(f_l, svec, f_k); /* 3 */
1617 rvec_inc(f[i], f_i); /* 3 */
1618 rvec_dec(f[j], f_j); /* 3 */
1619 rvec_dec(f[k], f_k); /* 3 */
1620 rvec_inc(f[l], f_l); /* 3 */
1624 copy_ivec(SHIFT_IVEC(g, j), jt);
1625 ivec_sub(SHIFT_IVEC(g, i), jt, dt_ij);
1626 ivec_sub(SHIFT_IVEC(g, k), jt, dt_kj);
1627 ivec_sub(SHIFT_IVEC(g, l), jt, dt_lj);
1628 t1 = IVEC2IS(dt_ij);
1629 t2 = IVEC2IS(dt_kj);
1630 t3 = IVEC2IS(dt_lj);
1634 t3 = pbc_rvec_sub(pbc, x[l], x[j], dx_jl);
1641 rvec_inc(fshift[t1], f_i);
1642 rvec_dec(fshift[CENTRAL], f_j);
1643 rvec_dec(fshift[t2], f_k);
1644 rvec_inc(fshift[t3], f_l);
1649 /* As do_dih_fup above, but without shift forces */
1651 do_dih_fup_noshiftf(int i, int j, int k, int l, real ddphi,
1652 rvec r_ij, rvec r_kj, rvec r_kl,
1653 rvec m, rvec n, rvec4 f[])
1655 rvec f_i, f_j, f_k, f_l;
1656 rvec uvec, vvec, svec;
1657 real iprm, iprn, nrkj, nrkj2, nrkj_1, nrkj_2;
1658 real a, b, p, q, toler;
1660 iprm = iprod(m, m); /* 5 */
1661 iprn = iprod(n, n); /* 5 */
1662 nrkj2 = iprod(r_kj, r_kj); /* 5 */
1663 toler = nrkj2*GMX_REAL_EPS;
1664 if ((iprm > toler) && (iprn > toler))
1666 nrkj_1 = gmx::invsqrt(nrkj2); /* 10 */
1667 nrkj_2 = nrkj_1*nrkj_1; /* 1 */
1668 nrkj = nrkj2*nrkj_1; /* 1 */
1669 a = -ddphi*nrkj/iprm; /* 11 */
1670 svmul(a, m, f_i); /* 3 */
1671 b = ddphi*nrkj/iprn; /* 11 */
1672 svmul(b, n, f_l); /* 3 */
1673 p = iprod(r_ij, r_kj); /* 5 */
1674 p *= nrkj_2; /* 1 */
1675 q = iprod(r_kl, r_kj); /* 5 */
1676 q *= nrkj_2; /* 1 */
1677 svmul(p, f_i, uvec); /* 3 */
1678 svmul(q, f_l, vvec); /* 3 */
1679 rvec_sub(uvec, vvec, svec); /* 3 */
1680 rvec_sub(f_i, svec, f_j); /* 3 */
1681 rvec_add(f_l, svec, f_k); /* 3 */
1682 rvec_inc(f[i], f_i); /* 3 */
1683 rvec_dec(f[j], f_j); /* 3 */
1684 rvec_dec(f[k], f_k); /* 3 */
1685 rvec_inc(f[l], f_l); /* 3 */
1689 #if GMX_SIMD_HAVE_REAL
1690 /* As do_dih_fup_noshiftf above, but with SIMD and pre-calculated pre-factors */
1691 static inline void gmx_simdcall
1692 do_dih_fup_noshiftf_simd(const int *ai, const int *aj, const int *ak, const int *al,
1693 SimdReal p, SimdReal q,
1694 SimdReal f_i_x, SimdReal f_i_y, SimdReal f_i_z,
1695 SimdReal mf_l_x, SimdReal mf_l_y, SimdReal mf_l_z,
1698 SimdReal sx = p * f_i_x + q * mf_l_x;
1699 SimdReal sy = p * f_i_y + q * mf_l_y;
1700 SimdReal sz = p * f_i_z + q * mf_l_z;
1701 SimdReal f_j_x = f_i_x - sx;
1702 SimdReal f_j_y = f_i_y - sy;
1703 SimdReal f_j_z = f_i_z - sz;
1704 SimdReal f_k_x = mf_l_x - sx;
1705 SimdReal f_k_y = mf_l_y - sy;
1706 SimdReal f_k_z = mf_l_z - sz;
1707 transposeScatterIncrU<4>(reinterpret_cast<real *>(f), ai, f_i_x, f_i_y, f_i_z);
1708 transposeScatterDecrU<4>(reinterpret_cast<real *>(f), aj, f_j_x, f_j_y, f_j_z);
1709 transposeScatterIncrU<4>(reinterpret_cast<real *>(f), ak, f_k_x, f_k_y, f_k_z);
1710 transposeScatterDecrU<4>(reinterpret_cast<real *>(f), al, mf_l_x, mf_l_y, mf_l_z);
1712 #endif // GMX_SIMD_HAVE_REAL
1714 static real dopdihs(real cpA, real cpB, real phiA, real phiB, int mult,
1715 real phi, real lambda, real *V, real *F)
1717 real v, dvdlambda, mdphi, v1, sdphi, ddphi;
1718 real L1 = 1.0 - lambda;
1719 real ph0 = (L1*phiA + lambda*phiB)*DEG2RAD;
1720 real dph0 = (phiB - phiA)*DEG2RAD;
1721 real cp = L1*cpA + lambda*cpB;
1723 mdphi = mult*phi - ph0;
1724 sdphi = std::sin(mdphi);
1725 ddphi = -cp*mult*sdphi;
1726 v1 = 1.0 + std::cos(mdphi);
1729 dvdlambda = (cpB - cpA)*v1 + cp*dph0*sdphi;
1736 /* That was 40 flops */
1740 dopdihs_noener(real cpA, real cpB, real phiA, real phiB, int mult,
1741 real phi, real lambda, real *F)
1743 real mdphi, sdphi, ddphi;
1744 real L1 = 1.0 - lambda;
1745 real ph0 = (L1*phiA + lambda*phiB)*DEG2RAD;
1746 real cp = L1*cpA + lambda*cpB;
1748 mdphi = mult*phi - ph0;
1749 sdphi = std::sin(mdphi);
1750 ddphi = -cp*mult*sdphi;
1754 /* That was 20 flops */
1757 static real dopdihs_min(real cpA, real cpB, real phiA, real phiB, int mult,
1758 real phi, real lambda, real *V, real *F)
1759 /* similar to dopdihs, except for a minus sign *
1760 * and a different treatment of mult/phi0 */
1762 real v, dvdlambda, mdphi, v1, sdphi, ddphi;
1763 real L1 = 1.0 - lambda;
1764 real ph0 = (L1*phiA + lambda*phiB)*DEG2RAD;
1765 real dph0 = (phiB - phiA)*DEG2RAD;
1766 real cp = L1*cpA + lambda*cpB;
1768 mdphi = mult*(phi-ph0);
1769 sdphi = std::sin(mdphi);
1770 ddphi = cp*mult*sdphi;
1771 v1 = 1.0-std::cos(mdphi);
1774 dvdlambda = (cpB-cpA)*v1 + cp*dph0*sdphi;
1781 /* That was 40 flops */
1784 real pdihs(int nbonds,
1785 const t_iatom forceatoms[], const t_iparams forceparams[],
1786 const rvec x[], rvec4 f[], rvec fshift[],
1787 const t_pbc *pbc, const t_graph *g,
1788 real lambda, real *dvdlambda,
1789 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
1790 int gmx_unused *global_atom_index)
1792 int i, type, ai, aj, ak, al;
1794 rvec r_ij, r_kj, r_kl, m, n;
1795 real phi, ddphi, vpd, vtot;
1799 for (i = 0; (i < nbonds); )
1801 type = forceatoms[i++];
1802 ai = forceatoms[i++];
1803 aj = forceatoms[i++];
1804 ak = forceatoms[i++];
1805 al = forceatoms[i++];
1807 phi = dih_angle(x[ai], x[aj], x[ak], x[al], pbc, r_ij, r_kj, r_kl, m, n,
1808 &t1, &t2, &t3); /* 84 */
1809 *dvdlambda += dopdihs(forceparams[type].pdihs.cpA,
1810 forceparams[type].pdihs.cpB,
1811 forceparams[type].pdihs.phiA,
1812 forceparams[type].pdihs.phiB,
1813 forceparams[type].pdihs.mult,
1814 phi, lambda, &vpd, &ddphi);
1817 do_dih_fup(ai, aj, ak, al, ddphi, r_ij, r_kj, r_kl, m, n,
1818 f, fshift, pbc, g, x, t1, t2, t3); /* 112 */
1825 void make_dp_periodic(real *dp) /* 1 flop? */
1827 /* dp cannot be outside (-pi,pi) */
1832 else if (*dp < -M_PI)
1838 /* As pdihs above, but without calculating energies and shift forces */
1840 pdihs_noener(int nbonds,
1841 const t_iatom forceatoms[], const t_iparams forceparams[],
1842 const rvec x[], rvec4 f[],
1843 const t_pbc gmx_unused *pbc, const t_graph gmx_unused *g,
1845 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
1846 int gmx_unused *global_atom_index)
1848 int i, type, ai, aj, ak, al;
1850 rvec r_ij, r_kj, r_kl, m, n;
1851 real phi, ddphi_tot, ddphi;
1853 for (i = 0; (i < nbonds); )
1855 ai = forceatoms[i+1];
1856 aj = forceatoms[i+2];
1857 ak = forceatoms[i+3];
1858 al = forceatoms[i+4];
1860 phi = dih_angle(x[ai], x[aj], x[ak], x[al], pbc, r_ij, r_kj, r_kl, m, n,
1865 /* Loop over dihedrals working on the same atoms,
1866 * so we avoid recalculating angles and force distributions.
1870 type = forceatoms[i];
1871 dopdihs_noener(forceparams[type].pdihs.cpA,
1872 forceparams[type].pdihs.cpB,
1873 forceparams[type].pdihs.phiA,
1874 forceparams[type].pdihs.phiB,
1875 forceparams[type].pdihs.mult,
1876 phi, lambda, &ddphi);
1881 while (i < nbonds &&
1882 forceatoms[i+1] == ai &&
1883 forceatoms[i+2] == aj &&
1884 forceatoms[i+3] == ak &&
1885 forceatoms[i+4] == al);
1887 do_dih_fup_noshiftf(ai, aj, ak, al, ddphi_tot, r_ij, r_kj, r_kl, m, n, f);
1892 #if GMX_SIMD_HAVE_REAL
1894 /* As pdihs_noner above, but using SIMD to calculate many dihedrals at once */
1896 pdihs_noener_simd(int nbonds,
1897 const t_iatom forceatoms[], const t_iparams forceparams[],
1898 const rvec x[], rvec4 f[],
1899 const t_pbc *pbc, const t_graph gmx_unused *g,
1900 real gmx_unused lambda,
1901 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
1902 int gmx_unused *global_atom_index)
1907 alignas(GMX_SIMD_ALIGNMENT) std::int32_t ai[GMX_SIMD_REAL_WIDTH];
1908 alignas(GMX_SIMD_ALIGNMENT) std::int32_t aj[GMX_SIMD_REAL_WIDTH];
1909 alignas(GMX_SIMD_ALIGNMENT) std::int32_t ak[GMX_SIMD_REAL_WIDTH];
1910 alignas(GMX_SIMD_ALIGNMENT) std::int32_t al[GMX_SIMD_REAL_WIDTH];
1911 alignas(GMX_SIMD_ALIGNMENT) real buf[3*GMX_SIMD_REAL_WIDTH];
1912 real *cp, *phi0, *mult;
1913 SimdReal deg2rad_S(DEG2RAD);
1915 SimdReal phi0_S, phi_S;
1916 SimdReal mx_S, my_S, mz_S;
1917 SimdReal nx_S, ny_S, nz_S;
1918 SimdReal nrkj_m2_S, nrkj_n2_S;
1919 SimdReal cp_S, mdphi_S, mult_S;
1920 SimdReal sin_S, cos_S;
1922 SimdReal sf_i_S, msf_l_S;
1923 alignas(GMX_SIMD_ALIGNMENT) real pbc_simd[9*GMX_SIMD_REAL_WIDTH];
1925 /* Extract aligned pointer for parameters and variables */
1926 cp = buf + 0*GMX_SIMD_REAL_WIDTH;
1927 phi0 = buf + 1*GMX_SIMD_REAL_WIDTH;
1928 mult = buf + 2*GMX_SIMD_REAL_WIDTH;
1930 set_pbc_simd(pbc, pbc_simd);
1932 /* nbonds is the number of dihedrals times nfa1, here we step GMX_SIMD_REAL_WIDTH dihs */
1933 for (i = 0; (i < nbonds); i += GMX_SIMD_REAL_WIDTH*nfa1)
1935 /* Collect atoms quadruplets for GMX_SIMD_REAL_WIDTH dihedrals.
1936 * iu indexes into forceatoms, we should not let iu go beyond nbonds.
1939 for (s = 0; s < GMX_SIMD_REAL_WIDTH; s++)
1941 type = forceatoms[iu];
1942 ai[s] = forceatoms[iu+1];
1943 aj[s] = forceatoms[iu+2];
1944 ak[s] = forceatoms[iu+3];
1945 al[s] = forceatoms[iu+4];
1947 /* At the end fill the arrays with the last atoms and 0 params */
1948 if (i + s*nfa1 < nbonds)
1950 cp[s] = forceparams[type].pdihs.cpA;
1951 phi0[s] = forceparams[type].pdihs.phiA;
1952 mult[s] = forceparams[type].pdihs.mult;
1954 if (iu + nfa1 < nbonds)
1967 /* Caclulate GMX_SIMD_REAL_WIDTH dihedral angles at once */
1968 dih_angle_simd(x, ai, aj, ak, al, pbc_simd,
1970 &mx_S, &my_S, &mz_S,
1971 &nx_S, &ny_S, &nz_S,
1976 cp_S = load<SimdReal>(cp);
1977 phi0_S = load<SimdReal>(phi0) * deg2rad_S;
1978 mult_S = load<SimdReal>(mult);
1980 mdphi_S = fms(mult_S, phi_S, phi0_S);
1982 /* Calculate GMX_SIMD_REAL_WIDTH sines at once */
1983 sincos(mdphi_S, &sin_S, &cos_S);
1984 mddphi_S = cp_S * mult_S * sin_S;
1985 sf_i_S = mddphi_S * nrkj_m2_S;
1986 msf_l_S = mddphi_S * nrkj_n2_S;
1988 /* After this m?_S will contain f[i] */
1989 mx_S = sf_i_S * mx_S;
1990 my_S = sf_i_S * my_S;
1991 mz_S = sf_i_S * mz_S;
1993 /* After this m?_S will contain -f[l] */
1994 nx_S = msf_l_S * nx_S;
1995 ny_S = msf_l_S * ny_S;
1996 nz_S = msf_l_S * nz_S;
1998 do_dih_fup_noshiftf_simd(ai, aj, ak, al,
2006 /* This is mostly a copy of pdihs_noener_simd above, but with using
2007 * the RB potential instead of a harmonic potential.
2008 * This function can replace rbdihs() when no energy and virial are needed.
2011 rbdihs_noener_simd(int nbonds,
2012 const t_iatom forceatoms[], const t_iparams forceparams[],
2013 const rvec x[], rvec4 f[],
2014 const t_pbc *pbc, const t_graph gmx_unused *g,
2015 real gmx_unused lambda,
2016 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
2017 int gmx_unused *global_atom_index)
2022 alignas(GMX_SIMD_ALIGNMENT) std::int32_t ai[GMX_SIMD_REAL_WIDTH];
2023 alignas(GMX_SIMD_ALIGNMENT) std::int32_t aj[GMX_SIMD_REAL_WIDTH];
2024 alignas(GMX_SIMD_ALIGNMENT) std::int32_t ak[GMX_SIMD_REAL_WIDTH];
2025 alignas(GMX_SIMD_ALIGNMENT) std::int32_t al[GMX_SIMD_REAL_WIDTH];
2026 alignas(GMX_SIMD_ALIGNMENT) real parm[NR_RBDIHS*GMX_SIMD_REAL_WIDTH];
2030 SimdReal ddphi_S, cosfac_S;
2031 SimdReal mx_S, my_S, mz_S;
2032 SimdReal nx_S, ny_S, nz_S;
2033 SimdReal nrkj_m2_S, nrkj_n2_S;
2034 SimdReal parm_S, c_S;
2035 SimdReal sin_S, cos_S;
2036 SimdReal sf_i_S, msf_l_S;
2037 alignas(GMX_SIMD_ALIGNMENT) real pbc_simd[9*GMX_SIMD_REAL_WIDTH];
2039 SimdReal pi_S(M_PI);
2040 SimdReal one_S(1.0);
2042 set_pbc_simd(pbc, pbc_simd);
2044 /* nbonds is the number of dihedrals times nfa1, here we step GMX_SIMD_REAL_WIDTH dihs */
2045 for (i = 0; (i < nbonds); i += GMX_SIMD_REAL_WIDTH*nfa1)
2047 /* Collect atoms quadruplets for GMX_SIMD_REAL_WIDTH dihedrals.
2048 * iu indexes into forceatoms, we should not let iu go beyond nbonds.
2051 for (s = 0; s < GMX_SIMD_REAL_WIDTH; s++)
2053 type = forceatoms[iu];
2054 ai[s] = forceatoms[iu+1];
2055 aj[s] = forceatoms[iu+2];
2056 ak[s] = forceatoms[iu+3];
2057 al[s] = forceatoms[iu+4];
2059 /* At the end fill the arrays with the last atoms and 0 params */
2060 if (i + s*nfa1 < nbonds)
2062 /* We don't need the first parameter, since that's a constant
2063 * which only affects the energies, not the forces.
2065 for (j = 1; j < NR_RBDIHS; j++)
2067 parm[j*GMX_SIMD_REAL_WIDTH + s] =
2068 forceparams[type].rbdihs.rbcA[j];
2071 if (iu + nfa1 < nbonds)
2078 for (j = 1; j < NR_RBDIHS; j++)
2080 parm[j*GMX_SIMD_REAL_WIDTH + s] = 0;
2085 /* Caclulate GMX_SIMD_REAL_WIDTH dihedral angles at once */
2086 dih_angle_simd(x, ai, aj, ak, al, pbc_simd,
2088 &mx_S, &my_S, &mz_S,
2089 &nx_S, &ny_S, &nz_S,
2094 /* Change to polymer convention */
2095 phi_S = phi_S - pi_S;
2097 sincos(phi_S, &sin_S, &cos_S);
2099 ddphi_S = setZero();
2102 for (j = 1; j < NR_RBDIHS; j++)
2104 parm_S = load<SimdReal>(parm + j*GMX_SIMD_REAL_WIDTH);
2105 ddphi_S = fma(c_S * parm_S, cosfac_S, ddphi_S);
2106 cosfac_S = cosfac_S * cos_S;
2110 /* Note that here we do not use the minus sign which is present
2111 * in the normal RB code. This is corrected for through (m)sf below.
2113 ddphi_S = ddphi_S * sin_S;
2115 sf_i_S = ddphi_S * nrkj_m2_S;
2116 msf_l_S = ddphi_S * nrkj_n2_S;
2118 /* After this m?_S will contain f[i] */
2119 mx_S = sf_i_S * mx_S;
2120 my_S = sf_i_S * my_S;
2121 mz_S = sf_i_S * mz_S;
2123 /* After this m?_S will contain -f[l] */
2124 nx_S = msf_l_S * nx_S;
2125 ny_S = msf_l_S * ny_S;
2126 nz_S = msf_l_S * nz_S;
2128 do_dih_fup_noshiftf_simd(ai, aj, ak, al,
2136 #endif // GMX_SIMD_HAVE_REAL
2139 real idihs(int nbonds,
2140 const t_iatom forceatoms[], const t_iparams forceparams[],
2141 const rvec x[], rvec4 f[], rvec fshift[],
2142 const t_pbc *pbc, const t_graph *g,
2143 real lambda, real *dvdlambda,
2144 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
2145 int gmx_unused *global_atom_index)
2147 int i, type, ai, aj, ak, al;
2149 real phi, phi0, dphi0, ddphi, vtot;
2150 rvec r_ij, r_kj, r_kl, m, n;
2151 real L1, kk, dp, dp2, kA, kB, pA, pB, dvdl_term;
2156 for (i = 0; (i < nbonds); )
2158 type = forceatoms[i++];
2159 ai = forceatoms[i++];
2160 aj = forceatoms[i++];
2161 ak = forceatoms[i++];
2162 al = forceatoms[i++];
2164 phi = dih_angle(x[ai], x[aj], x[ak], x[al], pbc, r_ij, r_kj, r_kl, m, n,
2165 &t1, &t2, &t3); /* 84 */
2167 /* phi can jump if phi0 is close to Pi/-Pi, which will cause huge
2168 * force changes if we just apply a normal harmonic.
2169 * Instead, we first calculate phi-phi0 and take it modulo (-Pi,Pi).
2170 * This means we will never have the periodicity problem, unless
2171 * the dihedral is Pi away from phiO, which is very unlikely due to
2174 kA = forceparams[type].harmonic.krA;
2175 kB = forceparams[type].harmonic.krB;
2176 pA = forceparams[type].harmonic.rA;
2177 pB = forceparams[type].harmonic.rB;
2179 kk = L1*kA + lambda*kB;
2180 phi0 = (L1*pA + lambda*pB)*DEG2RAD;
2181 dphi0 = (pB - pA)*DEG2RAD;
2185 make_dp_periodic(&dp);
2192 dvdl_term += 0.5*(kB - kA)*dp2 - kk*dphi0*dp;
2194 do_dih_fup(ai, aj, ak, al, -ddphi, r_ij, r_kj, r_kl, m, n,
2195 f, fshift, pbc, g, x, t1, t2, t3); /* 112 */
2199 *dvdlambda += dvdl_term;
2203 static real low_angres(int nbonds,
2204 const t_iatom forceatoms[], const t_iparams forceparams[],
2205 const rvec x[], rvec4 f[], rvec fshift[],
2206 const t_pbc *pbc, const t_graph *g,
2207 real lambda, real *dvdlambda,
2210 int i, m, type, ai, aj, ak, al;
2212 real phi, cos_phi, cos_phi2, vid, vtot, dVdphi;
2213 rvec r_ij, r_kl, f_i, f_k = {0, 0, 0};
2214 real st, sth, nrij2, nrkl2, c, cij, ckl;
2217 t2 = 0; /* avoid warning with gcc-3.3. It is never used uninitialized */
2220 ak = al = 0; /* to avoid warnings */
2221 for (i = 0; i < nbonds; )
2223 type = forceatoms[i++];
2224 ai = forceatoms[i++];
2225 aj = forceatoms[i++];
2226 t1 = pbc_rvec_sub(pbc, x[aj], x[ai], r_ij); /* 3 */
2229 ak = forceatoms[i++];
2230 al = forceatoms[i++];
2231 t2 = pbc_rvec_sub(pbc, x[al], x[ak], r_kl); /* 3 */
2240 cos_phi = cos_angle(r_ij, r_kl); /* 25 */
2241 phi = std::acos(cos_phi); /* 10 */
2243 *dvdlambda += dopdihs_min(forceparams[type].pdihs.cpA,
2244 forceparams[type].pdihs.cpB,
2245 forceparams[type].pdihs.phiA,
2246 forceparams[type].pdihs.phiB,
2247 forceparams[type].pdihs.mult,
2248 phi, lambda, &vid, &dVdphi); /* 40 */
2252 cos_phi2 = gmx::square(cos_phi); /* 1 */
2255 st = -dVdphi*gmx::invsqrt(1 - cos_phi2); /* 12 */
2256 sth = st*cos_phi; /* 1 */
2257 nrij2 = iprod(r_ij, r_ij); /* 5 */
2258 nrkl2 = iprod(r_kl, r_kl); /* 5 */
2260 c = st*gmx::invsqrt(nrij2*nrkl2); /* 11 */
2261 cij = sth/nrij2; /* 10 */
2262 ckl = sth/nrkl2; /* 10 */
2264 for (m = 0; m < DIM; m++) /* 18+18 */
2266 f_i[m] = (c*r_kl[m]-cij*r_ij[m]);
2271 f_k[m] = (c*r_ij[m]-ckl*r_kl[m]);
2279 ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, aj), dt);
2282 rvec_inc(fshift[t1], f_i);
2283 rvec_dec(fshift[CENTRAL], f_i);
2288 ivec_sub(SHIFT_IVEC(g, ak), SHIFT_IVEC(g, al), dt);
2291 rvec_inc(fshift[t2], f_k);
2292 rvec_dec(fshift[CENTRAL], f_k);
2297 return vtot; /* 184 / 157 (bZAxis) total */
2300 real angres(int nbonds,
2301 const t_iatom forceatoms[], const t_iparams forceparams[],
2302 const rvec x[], rvec4 f[], rvec fshift[],
2303 const t_pbc *pbc, const t_graph *g,
2304 real lambda, real *dvdlambda,
2305 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
2306 int gmx_unused *global_atom_index)
2308 return low_angres(nbonds, forceatoms, forceparams, x, f, fshift, pbc, g,
2309 lambda, dvdlambda, FALSE);
2312 real angresz(int nbonds,
2313 const t_iatom forceatoms[], const t_iparams forceparams[],
2314 const rvec x[], rvec4 f[], rvec fshift[],
2315 const t_pbc *pbc, const t_graph *g,
2316 real lambda, real *dvdlambda,
2317 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
2318 int gmx_unused *global_atom_index)
2320 return low_angres(nbonds, forceatoms, forceparams, x, f, fshift, pbc, g,
2321 lambda, dvdlambda, TRUE);
2324 real dihres(int nbonds,
2325 const t_iatom forceatoms[], const t_iparams forceparams[],
2326 const rvec x[], rvec4 f[], rvec fshift[],
2327 const t_pbc *pbc, const t_graph *g,
2328 real lambda, real *dvdlambda,
2329 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
2330 int gmx_unused *global_atom_index)
2333 int ai, aj, ak, al, i, type, t1, t2, t3;
2334 real phi0A, phi0B, dphiA, dphiB, kfacA, kfacB, phi0, dphi, kfac;
2335 real phi, ddphi, ddp, ddp2, dp, d2r, L1;
2336 rvec r_ij, r_kj, r_kl, m, n;
2342 for (i = 0; (i < nbonds); )
2344 type = forceatoms[i++];
2345 ai = forceatoms[i++];
2346 aj = forceatoms[i++];
2347 ak = forceatoms[i++];
2348 al = forceatoms[i++];
2350 phi0A = forceparams[type].dihres.phiA*d2r;
2351 dphiA = forceparams[type].dihres.dphiA*d2r;
2352 kfacA = forceparams[type].dihres.kfacA;
2354 phi0B = forceparams[type].dihres.phiB*d2r;
2355 dphiB = forceparams[type].dihres.dphiB*d2r;
2356 kfacB = forceparams[type].dihres.kfacB;
2358 phi0 = L1*phi0A + lambda*phi0B;
2359 dphi = L1*dphiA + lambda*dphiB;
2360 kfac = L1*kfacA + lambda*kfacB;
2362 phi = dih_angle(x[ai], x[aj], x[ak], x[al], pbc, r_ij, r_kj, r_kl, m, n,
2366 /* phi can jump if phi0 is close to Pi/-Pi, which will cause huge
2367 * force changes if we just apply a normal harmonic.
2368 * Instead, we first calculate phi-phi0 and take it modulo (-Pi,Pi).
2369 * This means we will never have the periodicity problem, unless
2370 * the dihedral is Pi away from phiO, which is very unlikely due to
2374 make_dp_periodic(&dp);
2380 else if (dp < -dphi)
2392 vtot += 0.5*kfac*ddp2;
2395 *dvdlambda += 0.5*(kfacB - kfacA)*ddp2;
2396 /* lambda dependence from changing restraint distances */
2399 *dvdlambda -= kfac*ddp*((dphiB - dphiA)+(phi0B - phi0A));
2403 *dvdlambda += kfac*ddp*((dphiB - dphiA)-(phi0B - phi0A));
2405 do_dih_fup(ai, aj, ak, al, ddphi, r_ij, r_kj, r_kl, m, n,
2406 f, fshift, pbc, g, x, t1, t2, t3); /* 112 */
2413 real unimplemented(int gmx_unused nbonds,
2414 const t_iatom gmx_unused forceatoms[], const t_iparams gmx_unused forceparams[],
2415 const rvec gmx_unused x[], rvec4 gmx_unused f[], rvec gmx_unused fshift[],
2416 const t_pbc gmx_unused *pbc, const t_graph gmx_unused *g,
2417 real gmx_unused lambda, real gmx_unused *dvdlambda,
2418 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
2419 int gmx_unused *global_atom_index)
2421 gmx_impl("*** you are using a not implemented function");
2424 real restrangles(int nbonds,
2425 const t_iatom forceatoms[], const t_iparams forceparams[],
2426 const rvec x[], rvec4 f[], rvec fshift[],
2427 const t_pbc *pbc, const t_graph *g,
2428 real gmx_unused lambda, real gmx_unused *dvdlambda,
2429 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
2430 int gmx_unused *global_atom_index)
2432 int i, d, ai, aj, ak, type, m;
2435 ivec jt, dt_ij, dt_kj;
2437 real prefactor, ratio_ante, ratio_post;
2438 rvec delta_ante, delta_post, vec_temp;
2441 for (i = 0; (i < nbonds); )
2443 type = forceatoms[i++];
2444 ai = forceatoms[i++];
2445 aj = forceatoms[i++];
2446 ak = forceatoms[i++];
2448 t1 = pbc_rvec_sub(pbc, x[ai], x[aj], vec_temp);
2449 pbc_rvec_sub(pbc, x[aj], x[ai], delta_ante);
2450 t2 = pbc_rvec_sub(pbc, x[ak], x[aj], delta_post);
2453 /* This function computes factors needed for restricted angle potential.
2454 * The restricted angle potential is used in coarse-grained simulations to avoid singularities
2455 * when three particles align and the dihedral angle and dihedral potential
2456 * cannot be calculated. This potential is calculated using the formula:
2457 real restrangles(int nbonds,
2458 const t_iatom forceatoms[],const t_iparams forceparams[],
2459 const rvec x[],rvec4 f[],rvec fshift[],
2460 const t_pbc *pbc,const t_graph *g,
2461 real gmx_unused lambda,real gmx_unused *dvdlambda,
2462 const t_mdatoms gmx_unused *md,t_fcdata gmx_unused *fcd,
2463 int gmx_unused *global_atom_index)
2465 int i, d, ai, aj, ak, type, m;
2469 ivec jt,dt_ij,dt_kj;
2471 real prefactor, ratio_ante, ratio_post;
2472 rvec delta_ante, delta_post, vec_temp;
2475 for(i=0; (i<nbonds); )
2477 type = forceatoms[i++];
2478 ai = forceatoms[i++];
2479 aj = forceatoms[i++];
2480 ak = forceatoms[i++];
2482 * \f[V_{\rm ReB}(\theta_i) = \frac{1}{2} k_{\theta} \frac{(\cos\theta_i - \cos\theta_0)^2}
2483 * {\sin^2\theta_i}\f] ({eq:ReB} and ref \cite{MonicaGoga2013} from the manual).
2484 * For more explanations see comments file "restcbt.h". */
2486 compute_factors_restangles(type, forceparams, delta_ante, delta_post,
2487 &prefactor, &ratio_ante, &ratio_post, &v);
2489 /* Forces are computed per component */
2490 for (d = 0; d < DIM; d++)
2492 f_i[d] = prefactor * (ratio_ante * delta_ante[d] - delta_post[d]);
2493 f_j[d] = prefactor * ((ratio_post + 1.0) * delta_post[d] - (ratio_ante + 1.0) * delta_ante[d]);
2494 f_k[d] = prefactor * (delta_ante[d] - ratio_post * delta_post[d]);
2497 /* Computation of potential energy */
2503 for (m = 0; (m < DIM); m++)
2512 copy_ivec(SHIFT_IVEC(g, aj), jt);
2513 ivec_sub(SHIFT_IVEC(g, ai), jt, dt_ij);
2514 ivec_sub(SHIFT_IVEC(g, ak), jt, dt_kj);
2515 t1 = IVEC2IS(dt_ij);
2516 t2 = IVEC2IS(dt_kj);
2519 rvec_inc(fshift[t1], f_i);
2520 rvec_inc(fshift[CENTRAL], f_j);
2521 rvec_inc(fshift[t2], f_k);
2527 real restrdihs(int nbonds,
2528 const t_iatom forceatoms[], const t_iparams forceparams[],
2529 const rvec x[], rvec4 f[], rvec fshift[],
2530 const t_pbc *pbc, const t_graph *g,
2531 real gmx_unused lambda, real gmx_unused *dvlambda,
2532 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
2533 int gmx_unused *global_atom_index)
2535 int i, d, type, ai, aj, ak, al;
2536 rvec f_i, f_j, f_k, f_l;
2538 ivec jt, dt_ij, dt_kj, dt_lj;
2541 rvec delta_ante, delta_crnt, delta_post, vec_temp;
2542 real factor_phi_ai_ante, factor_phi_ai_crnt, factor_phi_ai_post;
2543 real factor_phi_aj_ante, factor_phi_aj_crnt, factor_phi_aj_post;
2544 real factor_phi_ak_ante, factor_phi_ak_crnt, factor_phi_ak_post;
2545 real factor_phi_al_ante, factor_phi_al_crnt, factor_phi_al_post;
2550 for (i = 0; (i < nbonds); )
2552 type = forceatoms[i++];
2553 ai = forceatoms[i++];
2554 aj = forceatoms[i++];
2555 ak = forceatoms[i++];
2556 al = forceatoms[i++];
2558 t1 = pbc_rvec_sub(pbc, x[ai], x[aj], vec_temp);
2559 pbc_rvec_sub(pbc, x[aj], x[ai], delta_ante);
2560 t2 = pbc_rvec_sub(pbc, x[ak], x[aj], delta_crnt);
2561 pbc_rvec_sub(pbc, x[ak], x[al], vec_temp);
2562 pbc_rvec_sub(pbc, x[al], x[ak], delta_post);
2564 /* This function computes factors needed for restricted angle potential.
2565 * The restricted angle potential is used in coarse-grained simulations to avoid singularities
2566 * when three particles align and the dihedral angle and dihedral potential cannot be calculated.
2567 * This potential is calculated using the formula:
2568 * \f[V_{\rm ReB}(\theta_i) = \frac{1}{2} k_{\theta}
2569 * \frac{(\cos\theta_i - \cos\theta_0)^2}{\sin^2\theta_i}\f]
2570 * ({eq:ReB} and ref \cite{MonicaGoga2013} from the manual).
2571 * For more explanations see comments file "restcbt.h" */
2573 compute_factors_restrdihs(type, forceparams,
2574 delta_ante, delta_crnt, delta_post,
2575 &factor_phi_ai_ante, &factor_phi_ai_crnt, &factor_phi_ai_post,
2576 &factor_phi_aj_ante, &factor_phi_aj_crnt, &factor_phi_aj_post,
2577 &factor_phi_ak_ante, &factor_phi_ak_crnt, &factor_phi_ak_post,
2578 &factor_phi_al_ante, &factor_phi_al_crnt, &factor_phi_al_post,
2579 &prefactor_phi, &v);
2582 /* Computation of forces per component */
2583 for (d = 0; d < DIM; d++)
2585 f_i[d] = prefactor_phi * (factor_phi_ai_ante * delta_ante[d] + factor_phi_ai_crnt * delta_crnt[d] + factor_phi_ai_post * delta_post[d]);
2586 f_j[d] = prefactor_phi * (factor_phi_aj_ante * delta_ante[d] + factor_phi_aj_crnt * delta_crnt[d] + factor_phi_aj_post * delta_post[d]);
2587 f_k[d] = prefactor_phi * (factor_phi_ak_ante * delta_ante[d] + factor_phi_ak_crnt * delta_crnt[d] + factor_phi_ak_post * delta_post[d]);
2588 f_l[d] = prefactor_phi * (factor_phi_al_ante * delta_ante[d] + factor_phi_al_crnt * delta_crnt[d] + factor_phi_al_post * delta_post[d]);
2590 /* Computation of the energy */
2596 /* Updating the forces */
2598 rvec_inc(f[ai], f_i);
2599 rvec_inc(f[aj], f_j);
2600 rvec_inc(f[ak], f_k);
2601 rvec_inc(f[al], f_l);
2604 /* Updating the fshift forces for the pressure coupling */
2607 copy_ivec(SHIFT_IVEC(g, aj), jt);
2608 ivec_sub(SHIFT_IVEC(g, ai), jt, dt_ij);
2609 ivec_sub(SHIFT_IVEC(g, ak), jt, dt_kj);
2610 ivec_sub(SHIFT_IVEC(g, al), jt, dt_lj);
2611 t1 = IVEC2IS(dt_ij);
2612 t2 = IVEC2IS(dt_kj);
2613 t3 = IVEC2IS(dt_lj);
2617 t3 = pbc_rvec_sub(pbc, x[al], x[aj], dx_jl);
2624 rvec_inc(fshift[t1], f_i);
2625 rvec_inc(fshift[CENTRAL], f_j);
2626 rvec_inc(fshift[t2], f_k);
2627 rvec_inc(fshift[t3], f_l);
2635 real cbtdihs(int nbonds,
2636 const t_iatom forceatoms[], const t_iparams forceparams[],
2637 const rvec x[], rvec4 f[], rvec fshift[],
2638 const t_pbc *pbc, const t_graph *g,
2639 real gmx_unused lambda, real gmx_unused *dvdlambda,
2640 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
2641 int gmx_unused *global_atom_index)
2643 int type, ai, aj, ak, al, i, d;
2647 rvec f_i, f_j, f_k, f_l;
2648 ivec jt, dt_ij, dt_kj, dt_lj;
2650 rvec delta_ante, delta_crnt, delta_post;
2651 rvec f_phi_ai, f_phi_aj, f_phi_ak, f_phi_al;
2652 rvec f_theta_ante_ai, f_theta_ante_aj, f_theta_ante_ak;
2653 rvec f_theta_post_aj, f_theta_post_ak, f_theta_post_al;
2659 for (i = 0; (i < nbonds); )
2661 type = forceatoms[i++];
2662 ai = forceatoms[i++];
2663 aj = forceatoms[i++];
2664 ak = forceatoms[i++];
2665 al = forceatoms[i++];
2668 t1 = pbc_rvec_sub(pbc, x[ai], x[aj], vec_temp);
2669 pbc_rvec_sub(pbc, x[aj], x[ai], delta_ante);
2670 t2 = pbc_rvec_sub(pbc, x[ak], x[aj], vec_temp);
2671 pbc_rvec_sub(pbc, x[ak], x[aj], delta_crnt);
2672 pbc_rvec_sub(pbc, x[ak], x[al], vec_temp);
2673 pbc_rvec_sub(pbc, x[al], x[ak], delta_post);
2675 /* \brief Compute factors for CBT potential
2676 * The combined bending-torsion potential goes to zero in a very smooth manner, eliminating the numerical
2677 * instabilities, when three coarse-grained particles align and the dihedral angle and standard
2678 * dihedral potentials cannot be calculated. The CBT potential is calculated using the formula:
2679 * \f[V_{\rm CBT}(\theta_{i-1}, \theta_i, \phi_i) = k_{\phi} \sin^3\theta_{i-1} \sin^3\theta_{i}
2680 * \sum_{n=0}^4 { a_n \cos^n\phi_i}\f] ({eq:CBT} and ref \cite{MonicaGoga2013} from the manual).
2681 * It contains in its expression not only the dihedral angle \f$\phi\f$
2682 * but also \f[\theta_{i-1}\f] (theta_ante bellow) and \f[\theta_{i}\f] (theta_post bellow)
2683 * --- the adjacent bending angles.
2684 * For more explanations see comments file "restcbt.h". */
2686 compute_factors_cbtdihs(type, forceparams, delta_ante, delta_crnt, delta_post,
2687 f_phi_ai, f_phi_aj, f_phi_ak, f_phi_al,
2688 f_theta_ante_ai, f_theta_ante_aj, f_theta_ante_ak,
2689 f_theta_post_aj, f_theta_post_ak, f_theta_post_al,
2693 /* Acumulate the resuts per beads */
2694 for (d = 0; d < DIM; d++)
2696 f_i[d] = f_phi_ai[d] + f_theta_ante_ai[d];
2697 f_j[d] = f_phi_aj[d] + f_theta_ante_aj[d] + f_theta_post_aj[d];
2698 f_k[d] = f_phi_ak[d] + f_theta_ante_ak[d] + f_theta_post_ak[d];
2699 f_l[d] = f_phi_al[d] + f_theta_post_al[d];
2702 /* Compute the potential energy */
2707 /* Updating the forces */
2708 rvec_inc(f[ai], f_i);
2709 rvec_inc(f[aj], f_j);
2710 rvec_inc(f[ak], f_k);
2711 rvec_inc(f[al], f_l);
2714 /* Updating the fshift forces for the pressure coupling */
2717 copy_ivec(SHIFT_IVEC(g, aj), jt);
2718 ivec_sub(SHIFT_IVEC(g, ai), jt, dt_ij);
2719 ivec_sub(SHIFT_IVEC(g, ak), jt, dt_kj);
2720 ivec_sub(SHIFT_IVEC(g, al), jt, dt_lj);
2721 t1 = IVEC2IS(dt_ij);
2722 t2 = IVEC2IS(dt_kj);
2723 t3 = IVEC2IS(dt_lj);
2727 t3 = pbc_rvec_sub(pbc, x[al], x[aj], dx_jl);
2734 rvec_inc(fshift[t1], f_i);
2735 rvec_inc(fshift[CENTRAL], f_j);
2736 rvec_inc(fshift[t2], f_k);
2737 rvec_inc(fshift[t3], f_l);
2743 real rbdihs(int nbonds,
2744 const t_iatom forceatoms[], const t_iparams forceparams[],
2745 const rvec x[], rvec4 f[], rvec fshift[],
2746 const t_pbc *pbc, const t_graph *g,
2747 real lambda, real *dvdlambda,
2748 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
2749 int gmx_unused *global_atom_index)
2751 const real c0 = 0.0, c1 = 1.0, c2 = 2.0, c3 = 3.0, c4 = 4.0, c5 = 5.0;
2752 int type, ai, aj, ak, al, i, j;
2754 rvec r_ij, r_kj, r_kl, m, n;
2755 real parmA[NR_RBDIHS];
2756 real parmB[NR_RBDIHS];
2757 real parm[NR_RBDIHS];
2758 real cos_phi, phi, rbp, rbpBA;
2759 real v, ddphi, sin_phi;
2761 real L1 = 1.0-lambda;
2765 for (i = 0; (i < nbonds); )
2767 type = forceatoms[i++];
2768 ai = forceatoms[i++];
2769 aj = forceatoms[i++];
2770 ak = forceatoms[i++];
2771 al = forceatoms[i++];
2773 phi = dih_angle(x[ai], x[aj], x[ak], x[al], pbc, r_ij, r_kj, r_kl, m, n,
2774 &t1, &t2, &t3); /* 84 */
2776 /* Change to polymer convention */
2783 phi -= M_PI; /* 1 */
2786 cos_phi = std::cos(phi);
2787 /* Beware of accuracy loss, cannot use 1-sqrt(cos^2) ! */
2788 sin_phi = std::sin(phi);
2790 for (j = 0; (j < NR_RBDIHS); j++)
2792 parmA[j] = forceparams[type].rbdihs.rbcA[j];
2793 parmB[j] = forceparams[type].rbdihs.rbcB[j];
2794 parm[j] = L1*parmA[j]+lambda*parmB[j];
2796 /* Calculate cosine powers */
2797 /* Calculate the energy */
2798 /* Calculate the derivative */
2801 dvdl_term += (parmB[0]-parmA[0]);
2806 rbpBA = parmB[1]-parmA[1];
2807 ddphi += rbp*cosfac;
2810 dvdl_term += cosfac*rbpBA;
2812 rbpBA = parmB[2]-parmA[2];
2813 ddphi += c2*rbp*cosfac;
2816 dvdl_term += cosfac*rbpBA;
2818 rbpBA = parmB[3]-parmA[3];
2819 ddphi += c3*rbp*cosfac;
2822 dvdl_term += cosfac*rbpBA;
2824 rbpBA = parmB[4]-parmA[4];
2825 ddphi += c4*rbp*cosfac;
2828 dvdl_term += cosfac*rbpBA;
2830 rbpBA = parmB[5]-parmA[5];
2831 ddphi += c5*rbp*cosfac;
2834 dvdl_term += cosfac*rbpBA;
2836 ddphi = -ddphi*sin_phi; /* 11 */
2838 do_dih_fup(ai, aj, ak, al, ddphi, r_ij, r_kj, r_kl, m, n,
2839 f, fshift, pbc, g, x, t1, t2, t3); /* 112 */
2842 *dvdlambda += dvdl_term;
2849 /*! \brief Mysterious undocumented function */
2851 cmap_setup_grid_index(int ip, int grid_spacing, int *ipm1, int *ipp1, int *ipp2)
2857 ip = ip + grid_spacing - 1;
2859 else if (ip > grid_spacing)
2861 ip = ip - grid_spacing - 1;
2870 im1 = grid_spacing - 1;
2872 else if (ip == grid_spacing-2)
2876 else if (ip == grid_spacing-1)
2891 cmap_dihs(int nbonds,
2892 const t_iatom forceatoms[], const t_iparams forceparams[],
2893 const gmx_cmap_t *cmap_grid,
2894 const rvec x[], rvec4 f[], rvec fshift[],
2895 const struct t_pbc *pbc, const struct t_graph *g,
2896 real gmx_unused lambda, real gmx_unused *dvdlambda,
2897 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
2898 int gmx_unused *global_atom_index)
2901 int ai, aj, ak, al, am;
2902 int a1i, a1j, a1k, a1l, a2i, a2j, a2k, a2l;
2904 int t11, t21, t31, t12, t22, t32;
2905 int iphi1, ip1m1, ip1p1, ip1p2;
2906 int iphi2, ip2m1, ip2p1, ip2p2;
2908 int pos1, pos2, pos3, pos4;
2910 real ty[4], ty1[4], ty2[4], ty12[4], tx[16];
2911 real phi1, cos_phi1, sin_phi1, xphi1;
2912 real phi2, cos_phi2, sin_phi2, xphi2;
2913 real dx, tt, tu, e, df1, df2, vtot;
2914 real ra21, rb21, rg21, rg1, rgr1, ra2r1, rb2r1, rabr1;
2915 real ra22, rb22, rg22, rg2, rgr2, ra2r2, rb2r2, rabr2;
2916 real fg1, hg1, fga1, hgb1, gaa1, gbb1;
2917 real fg2, hg2, fga2, hgb2, gaa2, gbb2;
2920 rvec r1_ij, r1_kj, r1_kl, m1, n1;
2921 rvec r2_ij, r2_kj, r2_kl, m2, n2;
2922 rvec f1_i, f1_j, f1_k, f1_l;
2923 rvec f2_i, f2_j, f2_k, f2_l;
2924 rvec a1, b1, a2, b2;
2925 rvec f1, g1, h1, f2, g2, h2;
2926 rvec dtf1, dtg1, dth1, dtf2, dtg2, dth2;
2927 ivec jt1, dt1_ij, dt1_kj, dt1_lj;
2928 ivec jt2, dt2_ij, dt2_kj, dt2_lj;
2930 int loop_index[4][4] = {
2937 /* Total CMAP energy */
2940 for (n = 0; n < nbonds; )
2942 /* Five atoms are involved in the two torsions */
2943 type = forceatoms[n++];
2944 ai = forceatoms[n++];
2945 aj = forceatoms[n++];
2946 ak = forceatoms[n++];
2947 al = forceatoms[n++];
2948 am = forceatoms[n++];
2950 /* Which CMAP type is this */
2951 const int cmapA = forceparams[type].cmap.cmapA;
2952 const real *cmapd = cmap_grid->cmapdata[cmapA].cmap.data();
2960 phi1 = dih_angle(x[a1i], x[a1j], x[a1k], x[a1l], pbc, r1_ij, r1_kj, r1_kl, m1, n1,
2961 &t11, &t21, &t31); /* 84 */
2963 cos_phi1 = std::cos(phi1);
2965 a1[0] = r1_ij[1]*r1_kj[2]-r1_ij[2]*r1_kj[1];
2966 a1[1] = r1_ij[2]*r1_kj[0]-r1_ij[0]*r1_kj[2];
2967 a1[2] = r1_ij[0]*r1_kj[1]-r1_ij[1]*r1_kj[0]; /* 9 */
2969 b1[0] = r1_kl[1]*r1_kj[2]-r1_kl[2]*r1_kj[1];
2970 b1[1] = r1_kl[2]*r1_kj[0]-r1_kl[0]*r1_kj[2];
2971 b1[2] = r1_kl[0]*r1_kj[1]-r1_kl[1]*r1_kj[0]; /* 9 */
2973 pbc_rvec_sub(pbc, x[a1l], x[a1k], h1);
2975 ra21 = iprod(a1, a1); /* 5 */
2976 rb21 = iprod(b1, b1); /* 5 */
2977 rg21 = iprod(r1_kj, r1_kj); /* 5 */
2983 rabr1 = sqrt(ra2r1*rb2r1);
2985 sin_phi1 = rg1 * rabr1 * iprod(a1, h1) * (-1);
2987 if (cos_phi1 < -0.5 || cos_phi1 > 0.5)
2989 phi1 = std::asin(sin_phi1);
2999 phi1 = -M_PI - phi1;
3005 phi1 = std::acos(cos_phi1);
3013 xphi1 = phi1 + M_PI; /* 1 */
3015 /* Second torsion */
3021 phi2 = dih_angle(x[a2i], x[a2j], x[a2k], x[a2l], pbc, r2_ij, r2_kj, r2_kl, m2, n2,
3022 &t12, &t22, &t32); /* 84 */
3024 cos_phi2 = std::cos(phi2);
3026 a2[0] = r2_ij[1]*r2_kj[2]-r2_ij[2]*r2_kj[1];
3027 a2[1] = r2_ij[2]*r2_kj[0]-r2_ij[0]*r2_kj[2];
3028 a2[2] = r2_ij[0]*r2_kj[1]-r2_ij[1]*r2_kj[0]; /* 9 */
3030 b2[0] = r2_kl[1]*r2_kj[2]-r2_kl[2]*r2_kj[1];
3031 b2[1] = r2_kl[2]*r2_kj[0]-r2_kl[0]*r2_kj[2];
3032 b2[2] = r2_kl[0]*r2_kj[1]-r2_kl[1]*r2_kj[0]; /* 9 */
3034 pbc_rvec_sub(pbc, x[a2l], x[a2k], h2);
3036 ra22 = iprod(a2, a2); /* 5 */
3037 rb22 = iprod(b2, b2); /* 5 */
3038 rg22 = iprod(r2_kj, r2_kj); /* 5 */
3044 rabr2 = sqrt(ra2r2*rb2r2);
3046 sin_phi2 = rg2 * rabr2 * iprod(a2, h2) * (-1);
3048 if (cos_phi2 < -0.5 || cos_phi2 > 0.5)
3050 phi2 = std::asin(sin_phi2);
3060 phi2 = -M_PI - phi2;
3066 phi2 = std::acos(cos_phi2);
3074 xphi2 = phi2 + M_PI; /* 1 */
3076 /* Range mangling */
3079 xphi1 = xphi1 + 2*M_PI;
3081 else if (xphi1 >= 2*M_PI)
3083 xphi1 = xphi1 - 2*M_PI;
3088 xphi2 = xphi2 + 2*M_PI;
3090 else if (xphi2 >= 2*M_PI)
3092 xphi2 = xphi2 - 2*M_PI;
3095 /* Number of grid points */
3096 dx = 2*M_PI / cmap_grid->grid_spacing;
3098 /* Where on the grid are we */
3099 iphi1 = static_cast<int>(xphi1/dx);
3100 iphi2 = static_cast<int>(xphi2/dx);
3102 iphi1 = cmap_setup_grid_index(iphi1, cmap_grid->grid_spacing, &ip1m1, &ip1p1, &ip1p2);
3103 iphi2 = cmap_setup_grid_index(iphi2, cmap_grid->grid_spacing, &ip2m1, &ip2p1, &ip2p2);
3105 pos1 = iphi1*cmap_grid->grid_spacing+iphi2;
3106 pos2 = ip1p1*cmap_grid->grid_spacing+iphi2;
3107 pos3 = ip1p1*cmap_grid->grid_spacing+ip2p1;
3108 pos4 = iphi1*cmap_grid->grid_spacing+ip2p1;
3110 ty[0] = cmapd[pos1*4];
3111 ty[1] = cmapd[pos2*4];
3112 ty[2] = cmapd[pos3*4];
3113 ty[3] = cmapd[pos4*4];
3115 ty1[0] = cmapd[pos1*4+1];
3116 ty1[1] = cmapd[pos2*4+1];
3117 ty1[2] = cmapd[pos3*4+1];
3118 ty1[3] = cmapd[pos4*4+1];
3120 ty2[0] = cmapd[pos1*4+2];
3121 ty2[1] = cmapd[pos2*4+2];
3122 ty2[2] = cmapd[pos3*4+2];
3123 ty2[3] = cmapd[pos4*4+2];
3125 ty12[0] = cmapd[pos1*4+3];
3126 ty12[1] = cmapd[pos2*4+3];
3127 ty12[2] = cmapd[pos3*4+3];
3128 ty12[3] = cmapd[pos4*4+3];
3130 /* Switch to degrees */
3131 dx = 360.0 / cmap_grid->grid_spacing;
3132 xphi1 = xphi1 * RAD2DEG;
3133 xphi2 = xphi2 * RAD2DEG;
3135 for (i = 0; i < 4; i++) /* 16 */
3138 tx[i+4] = ty1[i]*dx;
3139 tx[i+8] = ty2[i]*dx;
3140 tx[i+12] = ty12[i]*dx*dx;
3144 for (int idx = 0; idx < 16; idx++) /* 1056 */
3146 for (int k = 0; k < 16; k++)
3148 tc[idx] += cmap_coeff_matrix[k*16+idx]*tx[k];
3152 tt = (xphi1-iphi1*dx)/dx;
3153 tu = (xphi2-iphi2*dx)/dx;
3159 for (i = 3; i >= 0; i--)
3161 l1 = loop_index[i][3];
3162 l2 = loop_index[i][2];
3163 l3 = loop_index[i][1];
3165 e = tt * e + ((tc[i*4+3]*tu+tc[i*4+2])*tu + tc[i*4+1])*tu+tc[i*4];
3166 df1 = tu * df1 + (3.0*tc[l1]*tt+2.0*tc[l2])*tt+tc[l3];
3167 df2 = tt * df2 + (3.0*tc[i*4+3]*tu+2.0*tc[i*4+2])*tu+tc[i*4+1];
3177 /* Do forces - first torsion */
3178 fg1 = iprod(r1_ij, r1_kj);
3179 hg1 = iprod(r1_kl, r1_kj);
3180 fga1 = fg1*ra2r1*rgr1;
3181 hgb1 = hg1*rb2r1*rgr1;
3185 for (i = 0; i < DIM; i++)
3187 dtf1[i] = gaa1 * a1[i];
3188 dtg1[i] = fga1 * a1[i] - hgb1 * b1[i];
3189 dth1[i] = gbb1 * b1[i];
3191 f1[i] = df1 * dtf1[i];
3192 g1[i] = df1 * dtg1[i];
3193 h1[i] = df1 * dth1[i];
3196 f1_j[i] = -f1[i] - g1[i];
3197 f1_k[i] = h1[i] + g1[i];
3200 f[a1i][i] = f[a1i][i] + f1_i[i];
3201 f[a1j][i] = f[a1j][i] + f1_j[i]; /* - f1[i] - g1[i] */
3202 f[a1k][i] = f[a1k][i] + f1_k[i]; /* h1[i] + g1[i] */
3203 f[a1l][i] = f[a1l][i] + f1_l[i]; /* h1[i] */
3206 /* Do forces - second torsion */
3207 fg2 = iprod(r2_ij, r2_kj);
3208 hg2 = iprod(r2_kl, r2_kj);
3209 fga2 = fg2*ra2r2*rgr2;
3210 hgb2 = hg2*rb2r2*rgr2;
3214 for (i = 0; i < DIM; i++)
3216 dtf2[i] = gaa2 * a2[i];
3217 dtg2[i] = fga2 * a2[i] - hgb2 * b2[i];
3218 dth2[i] = gbb2 * b2[i];
3220 f2[i] = df2 * dtf2[i];
3221 g2[i] = df2 * dtg2[i];
3222 h2[i] = df2 * dth2[i];
3225 f2_j[i] = -f2[i] - g2[i];
3226 f2_k[i] = h2[i] + g2[i];
3229 f[a2i][i] = f[a2i][i] + f2_i[i]; /* f2[i] */
3230 f[a2j][i] = f[a2j][i] + f2_j[i]; /* - f2[i] - g2[i] */
3231 f[a2k][i] = f[a2k][i] + f2_k[i]; /* h2[i] + g2[i] */
3232 f[a2l][i] = f[a2l][i] + f2_l[i]; /* - h2[i] */
3238 copy_ivec(SHIFT_IVEC(g, a1j), jt1);
3239 ivec_sub(SHIFT_IVEC(g, a1i), jt1, dt1_ij);
3240 ivec_sub(SHIFT_IVEC(g, a1k), jt1, dt1_kj);
3241 ivec_sub(SHIFT_IVEC(g, a1l), jt1, dt1_lj);
3242 t11 = IVEC2IS(dt1_ij);
3243 t21 = IVEC2IS(dt1_kj);
3244 t31 = IVEC2IS(dt1_lj);
3246 copy_ivec(SHIFT_IVEC(g, a2j), jt2);
3247 ivec_sub(SHIFT_IVEC(g, a2i), jt2, dt2_ij);
3248 ivec_sub(SHIFT_IVEC(g, a2k), jt2, dt2_kj);
3249 ivec_sub(SHIFT_IVEC(g, a2l), jt2, dt2_lj);
3250 t12 = IVEC2IS(dt2_ij);
3251 t22 = IVEC2IS(dt2_kj);
3252 t32 = IVEC2IS(dt2_lj);
3256 t31 = pbc_rvec_sub(pbc, x[a1l], x[a1j], h1);
3257 t32 = pbc_rvec_sub(pbc, x[a2l], x[a2j], h2);
3265 rvec_inc(fshift[t11], f1_i);
3266 rvec_inc(fshift[CENTRAL], f1_j);
3267 rvec_inc(fshift[t21], f1_k);
3268 rvec_inc(fshift[t31], f1_l);
3270 rvec_inc(fshift[t21], f2_i);
3271 rvec_inc(fshift[CENTRAL], f2_j);
3272 rvec_inc(fshift[t22], f2_k);
3273 rvec_inc(fshift[t32], f2_l);
3280 /***********************************************************
3282 * G R O M O S 9 6 F U N C T I O N S
3284 ***********************************************************/
3285 static real g96harmonic(real kA, real kB, real xA, real xB, real x, real lambda,
3288 const real half = 0.5;
3289 real L1, kk, x0, dx, dx2;
3290 real v, f, dvdlambda;
3293 kk = L1*kA+lambda*kB;
3294 x0 = L1*xA+lambda*xB;
3301 dvdlambda = half*(kB-kA)*dx2 + (xA-xB)*kk*dx;
3308 /* That was 21 flops */
3311 real g96bonds(int nbonds,
3312 const t_iatom forceatoms[], const t_iparams forceparams[],
3313 const rvec x[], rvec4 f[], rvec fshift[],
3314 const t_pbc *pbc, const t_graph *g,
3315 real lambda, real *dvdlambda,
3316 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
3317 int gmx_unused *global_atom_index)
3319 int i, m, ki, ai, aj, type;
3320 real dr2, fbond, vbond, fij, vtot;
3325 for (i = 0; (i < nbonds); )
3327 type = forceatoms[i++];
3328 ai = forceatoms[i++];
3329 aj = forceatoms[i++];
3331 ki = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /* 3 */
3332 dr2 = iprod(dx, dx); /* 5 */
3334 *dvdlambda += g96harmonic(forceparams[type].harmonic.krA,
3335 forceparams[type].harmonic.krB,
3336 forceparams[type].harmonic.rA,
3337 forceparams[type].harmonic.rB,
3338 dr2, lambda, &vbond, &fbond);
3340 vtot += 0.5*vbond; /* 1*/
3344 ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, aj), dt);
3347 for (m = 0; (m < DIM); m++) /* 15 */
3352 fshift[ki][m] += fij;
3353 fshift[CENTRAL][m] -= fij;
3359 static real g96bond_angle(const rvec xi, const rvec xj, const rvec xk, const t_pbc *pbc,
3360 rvec r_ij, rvec r_kj,
3362 /* Return value is the angle between the bonds i-j and j-k */
3366 *t1 = pbc_rvec_sub(pbc, xi, xj, r_ij); /* 3 */
3367 *t2 = pbc_rvec_sub(pbc, xk, xj, r_kj); /* 3 */
3369 costh = cos_angle(r_ij, r_kj); /* 25 */
3374 real g96angles(int nbonds,
3375 const t_iatom forceatoms[], const t_iparams forceparams[],
3376 const rvec x[], rvec4 f[], rvec fshift[],
3377 const t_pbc *pbc, const t_graph *g,
3378 real lambda, real *dvdlambda,
3379 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
3380 int gmx_unused *global_atom_index)
3382 int i, ai, aj, ak, type, m, t1, t2;
3384 real cos_theta, dVdt, va, vtot;
3385 real rij_1, rij_2, rkj_1, rkj_2, rijrkj_1;
3387 ivec jt, dt_ij, dt_kj;
3390 for (i = 0; (i < nbonds); )
3392 type = forceatoms[i++];
3393 ai = forceatoms[i++];
3394 aj = forceatoms[i++];
3395 ak = forceatoms[i++];
3397 cos_theta = g96bond_angle(x[ai], x[aj], x[ak], pbc, r_ij, r_kj, &t1, &t2);
3399 *dvdlambda += g96harmonic(forceparams[type].harmonic.krA,
3400 forceparams[type].harmonic.krB,
3401 forceparams[type].harmonic.rA,
3402 forceparams[type].harmonic.rB,
3403 cos_theta, lambda, &va, &dVdt);
3406 rij_1 = gmx::invsqrt(iprod(r_ij, r_ij));
3407 rkj_1 = gmx::invsqrt(iprod(r_kj, r_kj));
3408 rij_2 = rij_1*rij_1;
3409 rkj_2 = rkj_1*rkj_1;
3410 rijrkj_1 = rij_1*rkj_1; /* 23 */
3412 for (m = 0; (m < DIM); m++) /* 42 */
3414 f_i[m] = dVdt*(r_kj[m]*rijrkj_1 - r_ij[m]*rij_2*cos_theta);
3415 f_k[m] = dVdt*(r_ij[m]*rijrkj_1 - r_kj[m]*rkj_2*cos_theta);
3416 f_j[m] = -f_i[m]-f_k[m];
3424 copy_ivec(SHIFT_IVEC(g, aj), jt);
3426 ivec_sub(SHIFT_IVEC(g, ai), jt, dt_ij);
3427 ivec_sub(SHIFT_IVEC(g, ak), jt, dt_kj);
3428 t1 = IVEC2IS(dt_ij);
3429 t2 = IVEC2IS(dt_kj);
3431 rvec_inc(fshift[t1], f_i);
3432 rvec_inc(fshift[CENTRAL], f_j);
3433 rvec_inc(fshift[t2], f_k); /* 9 */
3439 real cross_bond_bond(int nbonds,
3440 const t_iatom forceatoms[], const t_iparams forceparams[],
3441 const rvec x[], rvec4 f[], rvec fshift[],
3442 const t_pbc *pbc, const t_graph *g,
3443 real gmx_unused lambda, real gmx_unused *dvdlambda,
3444 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
3445 int gmx_unused *global_atom_index)
3447 /* Potential from Lawrence and Skimmer, Chem. Phys. Lett. 372 (2003)
3450 int i, ai, aj, ak, type, m, t1, t2;
3452 real vtot, vrr, s1, s2, r1, r2, r1e, r2e, krr;
3454 ivec jt, dt_ij, dt_kj;
3457 for (i = 0; (i < nbonds); )
3459 type = forceatoms[i++];
3460 ai = forceatoms[i++];
3461 aj = forceatoms[i++];
3462 ak = forceatoms[i++];
3463 r1e = forceparams[type].cross_bb.r1e;
3464 r2e = forceparams[type].cross_bb.r2e;
3465 krr = forceparams[type].cross_bb.krr;
3467 /* Compute distance vectors ... */
3468 t1 = pbc_rvec_sub(pbc, x[ai], x[aj], r_ij);
3469 t2 = pbc_rvec_sub(pbc, x[ak], x[aj], r_kj);
3471 /* ... and their lengths */
3475 /* Deviations from ideality */
3479 /* Energy (can be negative!) */
3484 svmul(-krr*s2/r1, r_ij, f_i);
3485 svmul(-krr*s1/r2, r_kj, f_k);
3487 for (m = 0; (m < DIM); m++) /* 12 */
3489 f_j[m] = -f_i[m] - f_k[m];
3498 copy_ivec(SHIFT_IVEC(g, aj), jt);
3500 ivec_sub(SHIFT_IVEC(g, ai), jt, dt_ij);
3501 ivec_sub(SHIFT_IVEC(g, ak), jt, dt_kj);
3502 t1 = IVEC2IS(dt_ij);
3503 t2 = IVEC2IS(dt_kj);
3505 rvec_inc(fshift[t1], f_i);
3506 rvec_inc(fshift[CENTRAL], f_j);
3507 rvec_inc(fshift[t2], f_k); /* 9 */
3513 real cross_bond_angle(int nbonds,
3514 const t_iatom forceatoms[], const t_iparams forceparams[],
3515 const rvec x[], rvec4 f[], rvec fshift[],
3516 const t_pbc *pbc, const t_graph *g,
3517 real gmx_unused lambda, real gmx_unused *dvdlambda,
3518 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
3519 int gmx_unused *global_atom_index)
3521 /* Potential from Lawrence and Skimmer, Chem. Phys. Lett. 372 (2003)
3524 int i, ai, aj, ak, type, m, t1, t2;
3525 rvec r_ij, r_kj, r_ik;
3526 real vtot, vrt, s1, s2, s3, r1, r2, r3, r1e, r2e, r3e, krt, k1, k2, k3;
3528 ivec jt, dt_ij, dt_kj;
3531 for (i = 0; (i < nbonds); )
3533 type = forceatoms[i++];
3534 ai = forceatoms[i++];
3535 aj = forceatoms[i++];
3536 ak = forceatoms[i++];
3537 r1e = forceparams[type].cross_ba.r1e;
3538 r2e = forceparams[type].cross_ba.r2e;
3539 r3e = forceparams[type].cross_ba.r3e;
3540 krt = forceparams[type].cross_ba.krt;
3542 /* Compute distance vectors ... */
3543 t1 = pbc_rvec_sub(pbc, x[ai], x[aj], r_ij);
3544 t2 = pbc_rvec_sub(pbc, x[ak], x[aj], r_kj);
3545 pbc_rvec_sub(pbc, x[ai], x[ak], r_ik);
3547 /* ... and their lengths */
3552 /* Deviations from ideality */
3557 /* Energy (can be negative!) */
3558 vrt = krt*s3*(s1+s2);
3564 k3 = -krt*(s1+s2)/r3;
3565 for (m = 0; (m < DIM); m++)
3567 f_i[m] = k1*r_ij[m] + k3*r_ik[m];
3568 f_k[m] = k2*r_kj[m] - k3*r_ik[m];
3569 f_j[m] = -f_i[m] - f_k[m];
3572 for (m = 0; (m < DIM); m++) /* 12 */
3582 copy_ivec(SHIFT_IVEC(g, aj), jt);
3584 ivec_sub(SHIFT_IVEC(g, ai), jt, dt_ij);
3585 ivec_sub(SHIFT_IVEC(g, ak), jt, dt_kj);
3586 t1 = IVEC2IS(dt_ij);
3587 t2 = IVEC2IS(dt_kj);
3589 rvec_inc(fshift[t1], f_i);
3590 rvec_inc(fshift[CENTRAL], f_j);
3591 rvec_inc(fshift[t2], f_k); /* 9 */
3597 static real bonded_tab(const char *type, int table_nr,
3598 const bondedtable_t *table, real kA, real kB, real r,
3599 real lambda, real *V, real *F)
3601 real k, tabscale, *VFtab, rt, eps, eps2, Yt, Ft, Geps, Heps2, Fp, VV, FF;
3605 k = (1.0 - lambda)*kA + lambda*kB;
3607 tabscale = table->scale;
3608 VFtab = table->data;
3611 n0 = static_cast<int>(rt);
3614 gmx_fatal(FARGS, "A tabulated %s interaction table number %d is out of the table range: r %f, between table indices %d and %d, table length %d",
3615 type, table_nr, r, n0, n0+1, table->n);
3622 Geps = VFtab[nnn+2]*eps;
3623 Heps2 = VFtab[nnn+3]*eps2;
3624 Fp = Ft + Geps + Heps2;
3626 FF = Fp + Geps + 2.0*Heps2;
3628 *F = -k*FF*tabscale;
3630 dvdlambda = (kB - kA)*VV;
3634 /* That was 22 flops */
3637 real tab_bonds(int nbonds,
3638 const t_iatom forceatoms[], const t_iparams forceparams[],
3639 const rvec x[], rvec4 f[], rvec fshift[],
3640 const t_pbc *pbc, const t_graph *g,
3641 real lambda, real *dvdlambda,
3642 const t_mdatoms gmx_unused *md, t_fcdata *fcd,
3643 int gmx_unused *global_atom_index)
3645 int i, m, ki, ai, aj, type, table;
3646 real dr, dr2, fbond, vbond, fij, vtot;
3651 for (i = 0; (i < nbonds); )
3653 type = forceatoms[i++];
3654 ai = forceatoms[i++];
3655 aj = forceatoms[i++];
3657 ki = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /* 3 */
3658 dr2 = iprod(dx, dx); /* 5 */
3659 dr = dr2*gmx::invsqrt(dr2); /* 10 */
3661 table = forceparams[type].tab.table;
3663 *dvdlambda += bonded_tab("bond", table,
3664 &fcd->bondtab[table],
3665 forceparams[type].tab.kA,
3666 forceparams[type].tab.kB,
3667 dr, lambda, &vbond, &fbond); /* 22 */
3675 vtot += vbond; /* 1*/
3676 fbond *= gmx::invsqrt(dr2); /* 6 */
3679 ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, aj), dt);
3682 for (m = 0; (m < DIM); m++) /* 15 */
3687 fshift[ki][m] += fij;
3688 fshift[CENTRAL][m] -= fij;
3694 real tab_angles(int nbonds,
3695 const t_iatom forceatoms[], const t_iparams forceparams[],
3696 const rvec x[], rvec4 f[], rvec fshift[],
3697 const t_pbc *pbc, const t_graph *g,
3698 real lambda, real *dvdlambda,
3699 const t_mdatoms gmx_unused *md, t_fcdata *fcd,
3700 int gmx_unused *global_atom_index)
3702 int i, ai, aj, ak, t1, t2, type, table;
3704 real cos_theta, cos_theta2, theta, dVdt, va, vtot;
3705 ivec jt, dt_ij, dt_kj;
3708 for (i = 0; (i < nbonds); )
3710 type = forceatoms[i++];
3711 ai = forceatoms[i++];
3712 aj = forceatoms[i++];
3713 ak = forceatoms[i++];
3715 theta = bond_angle(x[ai], x[aj], x[ak], pbc,
3716 r_ij, r_kj, &cos_theta, &t1, &t2); /* 41 */
3718 table = forceparams[type].tab.table;
3720 *dvdlambda += bonded_tab("angle", table,
3721 &fcd->angletab[table],
3722 forceparams[type].tab.kA,
3723 forceparams[type].tab.kB,
3724 theta, lambda, &va, &dVdt); /* 22 */
3727 cos_theta2 = gmx::square(cos_theta); /* 1 */
3736 st = dVdt*gmx::invsqrt(1 - cos_theta2); /* 12 */
3737 sth = st*cos_theta; /* 1 */
3738 nrkj2 = iprod(r_kj, r_kj); /* 5 */
3739 nrij2 = iprod(r_ij, r_ij);
3741 cik = st*gmx::invsqrt(nrkj2*nrij2); /* 12 */
3742 cii = sth/nrij2; /* 10 */
3743 ckk = sth/nrkj2; /* 10 */
3745 for (m = 0; (m < DIM); m++) /* 39 */
3747 f_i[m] = -(cik*r_kj[m]-cii*r_ij[m]);
3748 f_k[m] = -(cik*r_ij[m]-ckk*r_kj[m]);
3749 f_j[m] = -f_i[m]-f_k[m];
3756 copy_ivec(SHIFT_IVEC(g, aj), jt);
3758 ivec_sub(SHIFT_IVEC(g, ai), jt, dt_ij);
3759 ivec_sub(SHIFT_IVEC(g, ak), jt, dt_kj);
3760 t1 = IVEC2IS(dt_ij);
3761 t2 = IVEC2IS(dt_kj);
3763 rvec_inc(fshift[t1], f_i);
3764 rvec_inc(fshift[CENTRAL], f_j);
3765 rvec_inc(fshift[t2], f_k);
3771 real tab_dihs(int nbonds,
3772 const t_iatom forceatoms[], const t_iparams forceparams[],
3773 const rvec x[], rvec4 f[], rvec fshift[],
3774 const t_pbc *pbc, const t_graph *g,
3775 real lambda, real *dvdlambda,
3776 const t_mdatoms gmx_unused *md, t_fcdata *fcd,
3777 int gmx_unused *global_atom_index)
3779 int i, type, ai, aj, ak, al, table;
3781 rvec r_ij, r_kj, r_kl, m, n;
3782 real phi, ddphi, vpd, vtot;
3785 for (i = 0; (i < nbonds); )
3787 type = forceatoms[i++];
3788 ai = forceatoms[i++];
3789 aj = forceatoms[i++];
3790 ak = forceatoms[i++];
3791 al = forceatoms[i++];
3793 phi = dih_angle(x[ai], x[aj], x[ak], x[al], pbc, r_ij, r_kj, r_kl, m, n,
3794 &t1, &t2, &t3); /* 84 */
3796 table = forceparams[type].tab.table;
3798 /* Hopefully phi+M_PI never results in values < 0 */
3799 *dvdlambda += bonded_tab("dihedral", table,
3800 &fcd->dihtab[table],
3801 forceparams[type].tab.kA,
3802 forceparams[type].tab.kB,
3803 phi+M_PI, lambda, &vpd, &ddphi);
3806 do_dih_fup(ai, aj, ak, al, -ddphi, r_ij, r_kj, r_kl, m, n,
3807 f, fshift, pbc, g, x, t1, t2, t3); /* 112 */