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,2019, 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/gmxlib/nrnb.h"
58 #include "gromacs/listed_forces/disre.h"
59 #include "gromacs/listed_forces/orires.h"
60 #include "gromacs/math/functions.h"
61 #include "gromacs/math/units.h"
62 #include "gromacs/math/utilities.h"
63 #include "gromacs/math/vec.h"
64 #include "gromacs/mdtypes/fcdata.h"
65 #include "gromacs/mdtypes/mdatom.h"
66 #include "gromacs/pbcutil/ishift.h"
67 #include "gromacs/pbcutil/mshift.h"
68 #include "gromacs/pbcutil/pbc.h"
69 #include "gromacs/pbcutil/pbc_simd.h"
70 #include "gromacs/simd/simd.h"
71 #include "gromacs/simd/simd_math.h"
72 #include "gromacs/simd/vector_operations.h"
73 #include "gromacs/utility/basedefinitions.h"
74 #include "gromacs/utility/enumerationhelpers.h"
75 #include "gromacs/utility/fatalerror.h"
76 #include "gromacs/utility/real.h"
77 #include "gromacs/utility/smalloc.h"
79 #include "listed_internal.h"
82 using namespace gmx
; // TODO: Remove when this file is moved into gmx namespace
87 //! Type of CPU function to compute a bonded interaction.
88 using BondedFunction
= real (*)(int nbonds
,
89 const t_iatom iatoms
[],
90 const t_iparams iparams
[],
102 /*! \brief Mysterious CMAP coefficient matrix */
103 const int cmap_coeff_matrix
[] = {
104 1, 0, -3, 2, 0, 0, 0, 0, -3, 0, 9, -6, 2, 0, -6, 4, 0, 0, 0, 0, 0, 0, 0, 0,
105 3, 0, -9, 6, -2, 0, 6, -4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 9, -6, 0, 0, -6, 4,
106 0, 0, 3, -2, 0, 0, 0, 0, 0, 0, -9, 6, 0, 0, 6, -4, 0, 0, 0, 0, 1, 0, -3, 2,
107 -2, 0, 6, -4, 1, 0, -3, 2, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 3, -2, 1, 0, -3, 2,
108 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 2, 0, 0, 3, -2, 0, 0, 0, 0, 0, 0, 3, -2,
109 0, 0, -6, 4, 0, 0, 3, -2, 0, 1, -2, 1, 0, 0, 0, 0, 0, -3, 6, -3, 0, 2, -4, 2,
110 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, -6, 3, 0, -2, 4, -2, 0, 0, 0, 0, 0, 0, 0, 0,
111 0, 0, -3, 3, 0, 0, 2, -2, 0, 0, -1, 1, 0, 0, 0, 0, 0, 0, 3, -3, 0, 0, -2, 2,
112 0, 0, 0, 0, 0, 1, -2, 1, 0, -2, 4, -2, 0, 1, -2, 1, 0, 0, 0, 0, 0, 0, 0, 0,
113 0, -1, 2, -1, 0, 1, -2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, -1, 1,
114 0, 0, 0, 0, 0, 0, -1, 1, 0, 0, 2, -2, 0, 0, -1, 1
118 /*! \brief Compute dx = xi - xj, modulo PBC if non-NULL
120 * \todo This kind of code appears in many places. Consolidate it */
121 int pbc_rvec_sub(const t_pbc
* pbc
, const rvec xi
, const rvec xj
, rvec dx
)
125 return pbc_dx_aiuc(pbc
, xi
, xj
, dx
);
129 rvec_sub(xi
, xj
, dx
);
137 real
bond_angle(const rvec xi
,
146 /* Return value is the angle between the bonds i-j and j-k */
151 *t1
= pbc_rvec_sub(pbc
, xi
, xj
, r_ij
); /* 3 */
152 *t2
= pbc_rvec_sub(pbc
, xk
, xj
, r_kj
); /* 3 */
154 *costh
= cos_angle(r_ij
, r_kj
); /* 25 */
155 th
= std::acos(*costh
); /* 10 */
160 real
dih_angle(const rvec xi
,
174 *t1
= pbc_rvec_sub(pbc
, xi
, xj
, r_ij
); /* 3 */
175 *t2
= pbc_rvec_sub(pbc
, xk
, xj
, r_kj
); /* 3 */
176 *t3
= pbc_rvec_sub(pbc
, xk
, xl
, r_kl
); /* 3 */
178 cprod(r_ij
, r_kj
, m
); /* 9 */
179 cprod(r_kj
, r_kl
, n
); /* 9 */
180 real phi
= gmx_angle(m
, n
); /* 49 (assuming 25 for atan2) */
181 real ipr
= iprod(r_ij
, n
); /* 5 */
182 real sign
= (ipr
< 0.0) ? -1.0 : 1.0;
183 phi
= sign
* phi
; /* 1 */
189 void make_dp_periodic(real
* dp
) /* 1 flop? */
191 /* dp cannot be outside (-pi,pi) */
196 else if (*dp
< -M_PI
)
205 /*! \brief Spreads and accumulates the bonded forces to the two atoms and adds the virial contribution when needed
207 * When \p g==nullptr, \p shiftIndex is used as the periodic shift.
208 * When \p g!=nullptr, the graph is used to compute the periodic shift.
210 template<BondedKernelFlavor flavor
>
211 inline void spreadBondForces(const real bondForce
,
220 if (computeVirial(flavor
) && g
)
223 ivec_sub(SHIFT_IVEC(g
, ai
), SHIFT_IVEC(g
, aj
), dt
);
224 shiftIndex
= IVEC2IS(dt
);
227 for (int m
= 0; m
< DIM
; m
++) /* 15 */
229 const real fij
= bondForce
* dx
[m
];
232 if (computeVirial(flavor
))
234 fshift
[shiftIndex
][m
] += fij
;
235 fshift
[CENTRAL
][m
] -= fij
;
240 /*! \brief Morse potential bond
242 * By Frank Everdij. Three parameters needed:
244 * b0 = equilibrium distance in nm
245 * be = beta in nm^-1 (actually, it's nu_e*Sqrt(2*pi*pi*mu/D_e))
246 * cb = well depth in kJ/mol
248 * Note: the potential is referenced to be +cb at infinite separation
249 * and zero at the equilibrium distance!
251 template<BondedKernelFlavor flavor
>
252 real
morse_bonds(int nbonds
,
253 const t_iatom forceatoms
[],
254 const t_iparams forceparams
[],
262 const t_mdatoms gmx_unused
* md
,
263 t_fcdata gmx_unused
* fcd
,
264 int gmx_unused
* global_atom_index
)
266 const real one
= 1.0;
267 const real two
= 2.0;
268 real dr
, dr2
, temp
, omtemp
, cbomtemp
, fbond
, vbond
, vtot
;
269 real b0
, be
, cb
, b0A
, beA
, cbA
, b0B
, beB
, cbB
, L1
;
271 int i
, ki
, type
, ai
, aj
;
274 for (i
= 0; (i
< nbonds
);)
276 type
= forceatoms
[i
++];
277 ai
= forceatoms
[i
++];
278 aj
= forceatoms
[i
++];
280 b0A
= forceparams
[type
].morse
.b0A
;
281 beA
= forceparams
[type
].morse
.betaA
;
282 cbA
= forceparams
[type
].morse
.cbA
;
284 b0B
= forceparams
[type
].morse
.b0B
;
285 beB
= forceparams
[type
].morse
.betaB
;
286 cbB
= forceparams
[type
].morse
.cbB
;
288 L1
= one
- lambda
; /* 1 */
289 b0
= L1
* b0A
+ lambda
* b0B
; /* 3 */
290 be
= L1
* beA
+ lambda
* beB
; /* 3 */
291 cb
= L1
* cbA
+ lambda
* cbB
; /* 3 */
293 ki
= pbc_rvec_sub(pbc
, x
[ai
], x
[aj
], dx
); /* 3 */
294 dr2
= iprod(dx
, dx
); /* 5 */
295 dr
= dr2
* gmx::invsqrt(dr2
); /* 10 */
296 temp
= std::exp(-be
* (dr
- b0
)); /* 12 */
300 /* bonds are constrainted. This may _not_ include bond constraints if they are lambda dependent */
301 *dvdlambda
+= cbB
- cbA
;
305 omtemp
= one
- temp
; /* 1 */
306 cbomtemp
= cb
* omtemp
; /* 1 */
307 vbond
= cbomtemp
* omtemp
; /* 1 */
308 fbond
= -two
* be
* temp
* cbomtemp
* gmx::invsqrt(dr2
); /* 9 */
309 vtot
+= vbond
; /* 1 */
311 *dvdlambda
+= (cbB
- cbA
) * omtemp
* omtemp
312 - (2 - 2 * omtemp
) * omtemp
* cb
* ((b0B
- b0A
) * be
- (beB
- beA
) * (dr
- b0
)); /* 15 */
314 spreadBondForces
<flavor
>(fbond
, dx
, ai
, aj
, f
, ki
, g
, fshift
); /* 15 */
320 template<BondedKernelFlavor flavor
>
321 real
cubic_bonds(int nbonds
,
322 const t_iatom forceatoms
[],
323 const t_iparams forceparams
[],
329 real gmx_unused lambda
,
330 real gmx_unused
* dvdlambda
,
331 const t_mdatoms gmx_unused
* md
,
332 t_fcdata gmx_unused
* fcd
,
333 int gmx_unused
* global_atom_index
)
335 const real three
= 3.0;
336 const real two
= 2.0;
338 real dr
, dr2
, dist
, kdist
, kdist2
, fbond
, vbond
, vtot
;
340 int i
, ki
, type
, ai
, aj
;
343 for (i
= 0; (i
< nbonds
);)
345 type
= forceatoms
[i
++];
346 ai
= forceatoms
[i
++];
347 aj
= forceatoms
[i
++];
349 b0
= forceparams
[type
].cubic
.b0
;
350 kb
= forceparams
[type
].cubic
.kb
;
351 kcub
= forceparams
[type
].cubic
.kcub
;
353 ki
= pbc_rvec_sub(pbc
, x
[ai
], x
[aj
], dx
); /* 3 */
354 dr2
= iprod(dx
, dx
); /* 5 */
361 dr
= dr2
* gmx::invsqrt(dr2
); /* 10 */
364 kdist2
= kdist
* dist
;
366 vbond
= kdist2
+ kcub
* kdist2
* dist
;
367 fbond
= -(two
* kdist
+ three
* kdist2
* kcub
) / dr
;
369 vtot
+= vbond
; /* 21 */
371 spreadBondForces
<flavor
>(fbond
, dx
, ai
, aj
, f
, ki
, g
, fshift
); /* 15 */
376 template<BondedKernelFlavor flavor
>
377 real
FENE_bonds(int nbonds
,
378 const t_iatom forceatoms
[],
379 const t_iparams forceparams
[],
385 real gmx_unused lambda
,
386 real gmx_unused
* dvdlambda
,
387 const t_mdatoms gmx_unused
* md
,
388 t_fcdata gmx_unused
* fcd
,
389 int* global_atom_index
)
391 const real half
= 0.5;
392 const real one
= 1.0;
394 real dr2
, bm2
, omdr2obm2
, fbond
, vbond
, vtot
;
396 int i
, ki
, type
, ai
, aj
;
399 for (i
= 0; (i
< nbonds
);)
401 type
= forceatoms
[i
++];
402 ai
= forceatoms
[i
++];
403 aj
= forceatoms
[i
++];
405 bm
= forceparams
[type
].fene
.bm
;
406 kb
= forceparams
[type
].fene
.kb
;
408 ki
= pbc_rvec_sub(pbc
, x
[ai
], x
[aj
], dx
); /* 3 */
409 dr2
= iprod(dx
, dx
); /* 5 */
420 gmx_fatal(FARGS
, "r^2 (%f) >= bm^2 (%f) in FENE bond between atoms %d and %d", dr2
, bm2
,
421 glatnr(global_atom_index
, ai
), glatnr(global_atom_index
, aj
));
424 omdr2obm2
= one
- dr2
/ bm2
;
426 vbond
= -half
* kb
* bm2
* std::log(omdr2obm2
);
427 fbond
= -kb
/ omdr2obm2
;
429 vtot
+= vbond
; /* 35 */
431 spreadBondForces
<flavor
>(fbond
, dx
, ai
, aj
, f
, ki
, g
, fshift
); /* 15 */
436 real
harmonic(real kA
, real kB
, real xA
, real xB
, real x
, real lambda
, real
* V
, real
* F
)
438 const real half
= 0.5;
439 real L1
, kk
, x0
, dx
, dx2
;
440 real v
, f
, dvdlambda
;
443 kk
= L1
* kA
+ lambda
* kB
;
444 x0
= L1
* xA
+ lambda
* xB
;
451 dvdlambda
= half
* (kB
- kA
) * dx2
+ (xA
- xB
) * kk
* dx
;
458 /* That was 19 flops */
462 template<BondedKernelFlavor flavor
>
463 real
bonds(int nbonds
,
464 const t_iatom forceatoms
[],
465 const t_iparams forceparams
[],
473 const t_mdatoms gmx_unused
* md
,
474 t_fcdata gmx_unused
* fcd
,
475 int gmx_unused
* global_atom_index
)
477 int i
, ki
, ai
, aj
, type
;
478 real dr
, dr2
, fbond
, vbond
, vtot
;
482 for (i
= 0; (i
< nbonds
);)
484 type
= forceatoms
[i
++];
485 ai
= forceatoms
[i
++];
486 aj
= forceatoms
[i
++];
488 ki
= pbc_rvec_sub(pbc
, x
[ai
], x
[aj
], dx
); /* 3 */
489 dr2
= iprod(dx
, dx
); /* 5 */
490 dr
= std::sqrt(dr2
); /* 10 */
492 *dvdlambda
+= harmonic(forceparams
[type
].harmonic
.krA
, forceparams
[type
].harmonic
.krB
,
493 forceparams
[type
].harmonic
.rA
, forceparams
[type
].harmonic
.rB
, dr
,
494 lambda
, &vbond
, &fbond
); /* 19 */
502 vtot
+= vbond
; /* 1*/
503 fbond
*= gmx::invsqrt(dr2
); /* 6 */
505 spreadBondForces
<flavor
>(fbond
, dx
, ai
, aj
, f
, ki
, g
, fshift
); /* 15 */
510 template<BondedKernelFlavor flavor
>
511 real
restraint_bonds(int nbonds
,
512 const t_iatom forceatoms
[],
513 const t_iparams forceparams
[],
521 const t_mdatoms gmx_unused
* md
,
522 t_fcdata gmx_unused
* fcd
,
523 int gmx_unused
* global_atom_index
)
525 int i
, ki
, ai
, aj
, type
;
526 real dr
, dr2
, fbond
, vbond
, vtot
;
528 real low
, dlow
, up1
, dup1
, up2
, dup2
, k
, dk
;
535 for (i
= 0; (i
< nbonds
);)
537 type
= forceatoms
[i
++];
538 ai
= forceatoms
[i
++];
539 aj
= forceatoms
[i
++];
541 ki
= pbc_rvec_sub(pbc
, x
[ai
], x
[aj
], dx
); /* 3 */
542 dr2
= iprod(dx
, dx
); /* 5 */
543 dr
= dr2
* gmx::invsqrt(dr2
); /* 10 */
545 low
= L1
* forceparams
[type
].restraint
.lowA
+ lambda
* forceparams
[type
].restraint
.lowB
;
546 dlow
= -forceparams
[type
].restraint
.lowA
+ forceparams
[type
].restraint
.lowB
;
547 up1
= L1
* forceparams
[type
].restraint
.up1A
+ lambda
* forceparams
[type
].restraint
.up1B
;
548 dup1
= -forceparams
[type
].restraint
.up1A
+ forceparams
[type
].restraint
.up1B
;
549 up2
= L1
* forceparams
[type
].restraint
.up2A
+ lambda
* forceparams
[type
].restraint
.up2B
;
550 dup2
= -forceparams
[type
].restraint
.up2A
+ forceparams
[type
].restraint
.up2B
;
551 k
= L1
* forceparams
[type
].restraint
.kA
+ lambda
* forceparams
[type
].restraint
.kB
;
552 dk
= -forceparams
[type
].restraint
.kA
+ forceparams
[type
].restraint
.kB
;
559 vbond
= 0.5 * k
* drh2
;
561 *dvdlambda
+= 0.5 * dk
* drh2
- k
* dlow
* drh
;
572 vbond
= 0.5 * k
* drh2
;
574 *dvdlambda
+= 0.5 * dk
* drh2
- k
* dup1
* drh
;
579 vbond
= k
* (up2
- up1
) * (0.5 * (up2
- up1
) + drh
);
580 fbond
= -k
* (up2
- up1
);
581 *dvdlambda
+= dk
* (up2
- up1
) * (0.5 * (up2
- up1
) + drh
)
582 + k
* (dup2
- dup1
) * (up2
- up1
+ drh
) - k
* (up2
- up1
) * dup2
;
590 vtot
+= vbond
; /* 1*/
591 fbond
*= gmx::invsqrt(dr2
); /* 6 */
593 spreadBondForces
<flavor
>(fbond
, dx
, ai
, aj
, f
, ki
, g
, fshift
); /* 15 */
599 template<BondedKernelFlavor flavor
>
600 real
polarize(int nbonds
,
601 const t_iatom forceatoms
[],
602 const t_iparams forceparams
[],
611 t_fcdata gmx_unused
* fcd
,
612 int gmx_unused
* global_atom_index
)
614 int i
, ki
, ai
, aj
, type
;
615 real dr
, dr2
, fbond
, vbond
, vtot
, ksh
;
619 for (i
= 0; (i
< nbonds
);)
621 type
= forceatoms
[i
++];
622 ai
= forceatoms
[i
++];
623 aj
= forceatoms
[i
++];
624 ksh
= gmx::square(md
->chargeA
[aj
]) * ONE_4PI_EPS0
/ forceparams
[type
].polarize
.alpha
;
626 ki
= pbc_rvec_sub(pbc
, x
[ai
], x
[aj
], dx
); /* 3 */
627 dr2
= iprod(dx
, dx
); /* 5 */
628 dr
= std::sqrt(dr2
); /* 10 */
630 *dvdlambda
+= harmonic(ksh
, ksh
, 0, 0, dr
, lambda
, &vbond
, &fbond
); /* 19 */
637 vtot
+= vbond
; /* 1*/
638 fbond
*= gmx::invsqrt(dr2
); /* 6 */
640 spreadBondForces
<flavor
>(fbond
, dx
, ai
, aj
, f
, ki
, g
, fshift
); /* 15 */
645 template<BondedKernelFlavor flavor
>
646 real
anharm_polarize(int nbonds
,
647 const t_iatom forceatoms
[],
648 const t_iparams forceparams
[],
657 t_fcdata gmx_unused
* fcd
,
658 int gmx_unused
* global_atom_index
)
660 int i
, ki
, ai
, aj
, type
;
661 real dr
, dr2
, fbond
, vbond
, vtot
, ksh
, khyp
, drcut
, ddr
, ddr3
;
665 for (i
= 0; (i
< nbonds
);)
667 type
= forceatoms
[i
++];
668 ai
= forceatoms
[i
++];
669 aj
= forceatoms
[i
++];
670 ksh
= gmx::square(md
->chargeA
[aj
]) * ONE_4PI_EPS0
/ forceparams
[type
].anharm_polarize
.alpha
; /* 7*/
671 khyp
= forceparams
[type
].anharm_polarize
.khyp
;
672 drcut
= forceparams
[type
].anharm_polarize
.drcut
;
674 ki
= pbc_rvec_sub(pbc
, x
[ai
], x
[aj
], dx
); /* 3 */
675 dr2
= iprod(dx
, dx
); /* 5 */
676 dr
= dr2
* gmx::invsqrt(dr2
); /* 10 */
678 *dvdlambda
+= harmonic(ksh
, ksh
, 0, 0, dr
, lambda
, &vbond
, &fbond
); /* 19 */
688 ddr3
= ddr
* ddr
* ddr
;
689 vbond
+= khyp
* ddr
* ddr3
;
690 fbond
-= 4 * khyp
* ddr3
;
692 fbond
*= gmx::invsqrt(dr2
); /* 6 */
693 vtot
+= vbond
; /* 1*/
695 spreadBondForces
<flavor
>(fbond
, dx
, ai
, aj
, f
, ki
, g
, fshift
); /* 15 */
700 template<BondedKernelFlavor flavor
>
701 real
water_pol(int nbonds
,
702 const t_iatom forceatoms
[],
703 const t_iparams forceparams
[],
706 rvec gmx_unused fshift
[],
707 const t_pbc gmx_unused
* pbc
,
708 const t_graph gmx_unused
* g
,
709 real gmx_unused lambda
,
710 real gmx_unused
* dvdlambda
,
711 const t_mdatoms gmx_unused
* md
,
712 t_fcdata gmx_unused
* fcd
,
713 int gmx_unused
* global_atom_index
)
715 /* This routine implements anisotropic polarizibility for water, through
716 * a shell connected to a dummy with spring constant that differ in the
717 * three spatial dimensions in the molecular frame.
719 int i
, m
, aO
, aH1
, aH2
, aD
, aS
, type
, type0
, ki
;
721 rvec dOH1
, dOH2
, dHH
, dOD
, dDS
, nW
, kk
, dx
, kdx
, proj
;
722 real vtot
, fij
, r_HH
, r_OD
, r_nW
, tx
, ty
, tz
, qS
;
727 type0
= forceatoms
[0];
729 qS
= md
->chargeA
[aS
];
730 kk
[XX
] = gmx::square(qS
) * ONE_4PI_EPS0
/ forceparams
[type0
].wpol
.al_x
;
731 kk
[YY
] = gmx::square(qS
) * ONE_4PI_EPS0
/ forceparams
[type0
].wpol
.al_y
;
732 kk
[ZZ
] = gmx::square(qS
) * ONE_4PI_EPS0
/ forceparams
[type0
].wpol
.al_z
;
733 r_HH
= 1.0 / forceparams
[type0
].wpol
.rHH
;
734 for (i
= 0; (i
< nbonds
); i
+= 6)
736 type
= forceatoms
[i
];
739 gmx_fatal(FARGS
, "Sorry, type = %d, type0 = %d, file = %s, line = %d", type
, type0
,
742 aO
= forceatoms
[i
+ 1];
743 aH1
= forceatoms
[i
+ 2];
744 aH2
= forceatoms
[i
+ 3];
745 aD
= forceatoms
[i
+ 4];
746 aS
= forceatoms
[i
+ 5];
748 /* Compute vectors describing the water frame */
749 pbc_rvec_sub(pbc
, x
[aH1
], x
[aO
], dOH1
);
750 pbc_rvec_sub(pbc
, x
[aH2
], x
[aO
], dOH2
);
751 pbc_rvec_sub(pbc
, x
[aH2
], x
[aH1
], dHH
);
752 pbc_rvec_sub(pbc
, x
[aD
], x
[aO
], dOD
);
753 ki
= pbc_rvec_sub(pbc
, x
[aS
], x
[aD
], dDS
);
754 cprod(dOH1
, dOH2
, nW
);
756 /* Compute inverse length of normal vector
757 * (this one could be precomputed, but I'm too lazy now)
759 r_nW
= gmx::invsqrt(iprod(nW
, nW
));
760 /* This is for precision, but does not make a big difference,
763 r_OD
= gmx::invsqrt(iprod(dOD
, dOD
));
765 /* Normalize the vectors in the water frame */
767 svmul(r_HH
, dHH
, dHH
);
768 svmul(r_OD
, dOD
, dOD
);
770 /* Compute displacement of shell along components of the vector */
771 dx
[ZZ
] = iprod(dDS
, dOD
);
772 /* Compute projection on the XY plane: dDS - dx[ZZ]*dOD */
773 for (m
= 0; (m
< DIM
); m
++)
775 proj
[m
] = dDS
[m
] - dx
[ZZ
] * dOD
[m
];
778 /*dx[XX] = iprod(dDS,nW);
779 dx[YY] = iprod(dDS,dHH);*/
780 dx
[XX
] = iprod(proj
, nW
);
781 for (m
= 0; (m
< DIM
); m
++)
783 proj
[m
] -= dx
[XX
] * nW
[m
];
785 dx
[YY
] = iprod(proj
, dHH
);
786 /* Now compute the forces and energy */
787 kdx
[XX
] = kk
[XX
] * dx
[XX
];
788 kdx
[YY
] = kk
[YY
] * dx
[YY
];
789 kdx
[ZZ
] = kk
[ZZ
] * dx
[ZZ
];
790 vtot
+= iprod(dx
, kdx
);
792 if (computeVirial(flavor
) && g
)
794 ivec_sub(SHIFT_IVEC(g
, aS
), SHIFT_IVEC(g
, aD
), dt
);
798 for (m
= 0; (m
< DIM
); m
++)
800 /* This is a tensor operation but written out for speed */
801 tx
= nW
[m
] * kdx
[XX
];
802 ty
= dHH
[m
] * kdx
[YY
];
803 tz
= dOD
[m
] * kdx
[ZZ
];
807 if (computeVirial(flavor
))
809 fshift
[ki
][m
] += fij
;
810 fshift
[CENTRAL
][m
] -= fij
;
818 template<BondedKernelFlavor flavor
>
820 do_1_thole(const rvec xi
, const rvec xj
, rvec fi
, rvec fj
, const t_pbc
* pbc
, real qq
, rvec fshift
[], real afac
)
823 real r12sq
, r12_1
, r12bar
, v0
, v1
, fscal
, ebar
, fff
;
826 t
= pbc_rvec_sub(pbc
, xi
, xj
, r12
); /* 3 */
828 r12sq
= iprod(r12
, r12
); /* 5 */
829 r12_1
= gmx::invsqrt(r12sq
); /* 5 */
830 r12bar
= afac
/ r12_1
; /* 5 */
831 v0
= qq
* ONE_4PI_EPS0
* r12_1
; /* 2 */
832 ebar
= std::exp(-r12bar
); /* 5 */
833 v1
= (1 - (1 + 0.5 * r12bar
) * ebar
); /* 4 */
834 fscal
= ((v0
* r12_1
) * v1
- v0
* 0.5 * afac
* ebar
* (r12bar
+ 1)) * r12_1
; /* 9 */
836 for (m
= 0; (m
< DIM
); m
++)
838 fff
= fscal
* r12
[m
];
841 if (computeVirial(flavor
))
844 fshift
[CENTRAL
][m
] -= fff
;
848 return v0
* v1
; /* 1 */
852 template<BondedKernelFlavor flavor
>
853 real
thole_pol(int nbonds
,
854 const t_iatom forceatoms
[],
855 const t_iparams forceparams
[],
860 const t_graph gmx_unused
* g
,
861 real gmx_unused lambda
,
862 real gmx_unused
* dvdlambda
,
864 t_fcdata gmx_unused
* fcd
,
865 int gmx_unused
* global_atom_index
)
867 /* Interaction between two pairs of particles with opposite charge */
868 int i
, type
, a1
, da1
, a2
, da2
;
869 real q1
, q2
, qq
, a
, al1
, al2
, afac
;
872 for (i
= 0; (i
< nbonds
);)
874 type
= forceatoms
[i
++];
875 a1
= forceatoms
[i
++];
876 da1
= forceatoms
[i
++];
877 a2
= forceatoms
[i
++];
878 da2
= forceatoms
[i
++];
879 q1
= md
->chargeA
[da1
];
880 q2
= md
->chargeA
[da2
];
881 a
= forceparams
[type
].thole
.a
;
882 al1
= forceparams
[type
].thole
.alpha1
;
883 al2
= forceparams
[type
].thole
.alpha2
;
885 afac
= a
* gmx::invsixthroot(al1
* al2
);
886 V
+= do_1_thole
<flavor
>(x
[a1
], x
[a2
], f
[a1
], f
[a2
], pbc
, qq
, fshift
, afac
);
887 V
+= do_1_thole
<flavor
>(x
[da1
], x
[a2
], f
[da1
], f
[a2
], pbc
, -qq
, fshift
, afac
);
888 V
+= do_1_thole
<flavor
>(x
[a1
], x
[da2
], f
[a1
], f
[da2
], pbc
, -qq
, fshift
, afac
);
889 V
+= do_1_thole
<flavor
>(x
[da1
], x
[da2
], f
[da1
], f
[da2
], pbc
, qq
, fshift
, afac
);
895 // Avoid gcc 386 -O3 code generation bug in this function (see Redmine
896 // #3205 for more information)
897 #if defined(__GNUC__) && defined(__i386__) && defined(__OPTIMIZE__)
898 # pragma GCC push_options
899 # pragma GCC optimize("O2")
900 # define avoid_gcc_i386_o3_code_generation_bug
903 template<BondedKernelFlavor flavor
>
904 std::enable_if_t
<flavor
!= BondedKernelFlavor::ForcesSimdWhenAvailable
|| !GMX_SIMD_HAVE_REAL
, real
>
906 const t_iatom forceatoms
[],
907 const t_iparams forceparams
[],
915 const t_mdatoms gmx_unused
* md
,
916 t_fcdata gmx_unused
* fcd
,
917 int gmx_unused
* global_atom_index
)
919 int i
, ai
, aj
, ak
, t1
, t2
, type
;
921 real cos_theta
, cos_theta2
, theta
, dVdt
, va
, vtot
;
922 ivec jt
, dt_ij
, dt_kj
;
925 for (i
= 0; i
< nbonds
;)
927 type
= forceatoms
[i
++];
928 ai
= forceatoms
[i
++];
929 aj
= forceatoms
[i
++];
930 ak
= forceatoms
[i
++];
932 theta
= bond_angle(x
[ai
], x
[aj
], x
[ak
], pbc
, r_ij
, r_kj
, &cos_theta
, &t1
, &t2
); /* 41 */
934 *dvdlambda
+= harmonic(forceparams
[type
].harmonic
.krA
, forceparams
[type
].harmonic
.krB
,
935 forceparams
[type
].harmonic
.rA
* DEG2RAD
,
936 forceparams
[type
].harmonic
.rB
* DEG2RAD
, theta
, lambda
, &va
, &dVdt
); /* 21 */
939 cos_theta2
= gmx::square(cos_theta
);
949 st
= dVdt
* gmx::invsqrt(1 - cos_theta2
); /* 12 */
950 sth
= st
* cos_theta
; /* 1 */
951 nrij2
= iprod(r_ij
, r_ij
); /* 5 */
952 nrkj2
= iprod(r_kj
, r_kj
); /* 5 */
954 nrij_1
= gmx::invsqrt(nrij2
); /* 10 */
955 nrkj_1
= gmx::invsqrt(nrkj2
); /* 10 */
957 cik
= st
* nrij_1
* nrkj_1
; /* 2 */
958 cii
= sth
* nrij_1
* nrij_1
; /* 2 */
959 ckk
= sth
* nrkj_1
* nrkj_1
; /* 2 */
961 for (m
= 0; m
< DIM
; m
++)
963 f_i
[m
] = -(cik
* r_kj
[m
] - cii
* r_ij
[m
]);
964 f_k
[m
] = -(cik
* r_ij
[m
] - ckk
* r_kj
[m
]);
965 f_j
[m
] = -f_i
[m
] - f_k
[m
];
970 if (computeVirial(flavor
))
974 copy_ivec(SHIFT_IVEC(g
, aj
), jt
);
976 ivec_sub(SHIFT_IVEC(g
, ai
), jt
, dt_ij
);
977 ivec_sub(SHIFT_IVEC(g
, ak
), jt
, dt_kj
);
981 rvec_inc(fshift
[t1
], f_i
);
982 rvec_inc(fshift
[CENTRAL
], f_j
);
983 rvec_inc(fshift
[t2
], f_k
);
991 #ifdef avoid_gcc_i386_o3_code_generation_bug
992 # pragma GCC pop_options
993 # undef avoid_gcc_i386_o3_code_generation_bug
996 #if GMX_SIMD_HAVE_REAL
998 /* As angles, but using SIMD to calculate many angles at once.
999 * This routines does not calculate energies and shift forces.
1001 template<BondedKernelFlavor flavor
>
1002 std::enable_if_t
<flavor
== BondedKernelFlavor::ForcesSimdWhenAvailable
, real
>
1004 const t_iatom forceatoms
[],
1005 const t_iparams forceparams
[],
1008 rvec gmx_unused fshift
[],
1010 const t_graph gmx_unused
* g
,
1011 real gmx_unused lambda
,
1012 real gmx_unused
* dvdlambda
,
1013 const t_mdatoms gmx_unused
* md
,
1014 t_fcdata gmx_unused
* fcd
,
1015 int gmx_unused
* global_atom_index
)
1020 alignas(GMX_SIMD_ALIGNMENT
) std::int32_t ai
[GMX_SIMD_REAL_WIDTH
];
1021 alignas(GMX_SIMD_ALIGNMENT
) std::int32_t aj
[GMX_SIMD_REAL_WIDTH
];
1022 alignas(GMX_SIMD_ALIGNMENT
) std::int32_t ak
[GMX_SIMD_REAL_WIDTH
];
1023 alignas(GMX_SIMD_ALIGNMENT
) real coeff
[2 * GMX_SIMD_REAL_WIDTH
];
1024 SimdReal
deg2rad_S(DEG2RAD
);
1025 SimdReal xi_S
, yi_S
, zi_S
;
1026 SimdReal xj_S
, yj_S
, zj_S
;
1027 SimdReal xk_S
, yk_S
, zk_S
;
1028 SimdReal k_S
, theta0_S
;
1029 SimdReal rijx_S
, rijy_S
, rijz_S
;
1030 SimdReal rkjx_S
, rkjy_S
, rkjz_S
;
1031 SimdReal
one_S(1.0);
1032 SimdReal
min_one_plus_eps_S(-1.0 + 2.0 * GMX_REAL_EPS
); // Smallest number > -1
1035 SimdReal nrij2_S
, nrij_1_S
;
1036 SimdReal nrkj2_S
, nrkj_1_S
;
1037 SimdReal cos_S
, invsin_S
;
1039 SimdReal st_S
, sth_S
;
1040 SimdReal cik_S
, cii_S
, ckk_S
;
1041 SimdReal f_ix_S
, f_iy_S
, f_iz_S
;
1042 SimdReal f_kx_S
, f_ky_S
, f_kz_S
;
1043 alignas(GMX_SIMD_ALIGNMENT
) real pbc_simd
[9 * GMX_SIMD_REAL_WIDTH
];
1045 set_pbc_simd(pbc
, pbc_simd
);
1047 /* nbonds is the number of angles times nfa1, here we step GMX_SIMD_REAL_WIDTH angles */
1048 for (i
= 0; (i
< nbonds
); i
+= GMX_SIMD_REAL_WIDTH
* nfa1
)
1050 /* Collect atoms for GMX_SIMD_REAL_WIDTH angles.
1051 * iu indexes into forceatoms, we should not let iu go beyond nbonds.
1054 for (s
= 0; s
< GMX_SIMD_REAL_WIDTH
; s
++)
1056 type
= forceatoms
[iu
];
1057 ai
[s
] = forceatoms
[iu
+ 1];
1058 aj
[s
] = forceatoms
[iu
+ 2];
1059 ak
[s
] = forceatoms
[iu
+ 3];
1061 /* At the end fill the arrays with the last atoms and 0 params */
1062 if (i
+ s
* nfa1
< nbonds
)
1064 coeff
[s
] = forceparams
[type
].harmonic
.krA
;
1065 coeff
[GMX_SIMD_REAL_WIDTH
+ s
] = forceparams
[type
].harmonic
.rA
;
1067 if (iu
+ nfa1
< nbonds
)
1075 coeff
[GMX_SIMD_REAL_WIDTH
+ s
] = 0;
1079 /* Store the non PBC corrected distances packed and aligned */
1080 gatherLoadUTranspose
<3>(reinterpret_cast<const real
*>(x
), ai
, &xi_S
, &yi_S
, &zi_S
);
1081 gatherLoadUTranspose
<3>(reinterpret_cast<const real
*>(x
), aj
, &xj_S
, &yj_S
, &zj_S
);
1082 gatherLoadUTranspose
<3>(reinterpret_cast<const real
*>(x
), ak
, &xk_S
, &yk_S
, &zk_S
);
1083 rijx_S
= xi_S
- xj_S
;
1084 rijy_S
= yi_S
- yj_S
;
1085 rijz_S
= zi_S
- zj_S
;
1086 rkjx_S
= xk_S
- xj_S
;
1087 rkjy_S
= yk_S
- yj_S
;
1088 rkjz_S
= zk_S
- zj_S
;
1090 k_S
= load
<SimdReal
>(coeff
);
1091 theta0_S
= load
<SimdReal
>(coeff
+ GMX_SIMD_REAL_WIDTH
) * deg2rad_S
;
1093 pbc_correct_dx_simd(&rijx_S
, &rijy_S
, &rijz_S
, pbc_simd
);
1094 pbc_correct_dx_simd(&rkjx_S
, &rkjy_S
, &rkjz_S
, pbc_simd
);
1096 rij_rkj_S
= iprod(rijx_S
, rijy_S
, rijz_S
, rkjx_S
, rkjy_S
, rkjz_S
);
1098 nrij2_S
= norm2(rijx_S
, rijy_S
, rijz_S
);
1099 nrkj2_S
= norm2(rkjx_S
, rkjy_S
, rkjz_S
);
1101 nrij_1_S
= invsqrt(nrij2_S
);
1102 nrkj_1_S
= invsqrt(nrkj2_S
);
1104 cos_S
= rij_rkj_S
* nrij_1_S
* nrkj_1_S
;
1106 /* To allow for 180 degrees, we take the max of cos and -1 + 1bit,
1107 * so we can safely get the 1/sin from 1/sqrt(1 - cos^2).
1108 * This also ensures that rounding errors would cause the argument
1109 * of simdAcos to be < -1.
1110 * Note that we do not take precautions for cos(0)=1, so the outer
1111 * atoms in an angle should not be on top of each other.
1113 cos_S
= max(cos_S
, min_one_plus_eps_S
);
1115 theta_S
= acos(cos_S
);
1117 invsin_S
= invsqrt(one_S
- cos_S
* cos_S
);
1119 st_S
= k_S
* (theta0_S
- theta_S
) * invsin_S
;
1120 sth_S
= st_S
* cos_S
;
1122 cik_S
= st_S
* nrij_1_S
* nrkj_1_S
;
1123 cii_S
= sth_S
* nrij_1_S
* nrij_1_S
;
1124 ckk_S
= sth_S
* nrkj_1_S
* nrkj_1_S
;
1126 f_ix_S
= cii_S
* rijx_S
;
1127 f_ix_S
= fnma(cik_S
, rkjx_S
, f_ix_S
);
1128 f_iy_S
= cii_S
* rijy_S
;
1129 f_iy_S
= fnma(cik_S
, rkjy_S
, f_iy_S
);
1130 f_iz_S
= cii_S
* rijz_S
;
1131 f_iz_S
= fnma(cik_S
, rkjz_S
, f_iz_S
);
1132 f_kx_S
= ckk_S
* rkjx_S
;
1133 f_kx_S
= fnma(cik_S
, rijx_S
, f_kx_S
);
1134 f_ky_S
= ckk_S
* rkjy_S
;
1135 f_ky_S
= fnma(cik_S
, rijy_S
, f_ky_S
);
1136 f_kz_S
= ckk_S
* rkjz_S
;
1137 f_kz_S
= fnma(cik_S
, rijz_S
, f_kz_S
);
1139 transposeScatterIncrU
<4>(reinterpret_cast<real
*>(f
), ai
, f_ix_S
, f_iy_S
, f_iz_S
);
1140 transposeScatterDecrU
<4>(reinterpret_cast<real
*>(f
), aj
, f_ix_S
+ f_kx_S
, f_iy_S
+ f_ky_S
,
1142 transposeScatterIncrU
<4>(reinterpret_cast<real
*>(f
), ak
, f_kx_S
, f_ky_S
, f_kz_S
);
1148 #endif // GMX_SIMD_HAVE_REAL
1150 template<BondedKernelFlavor flavor
>
1151 real
linear_angles(int nbonds
,
1152 const t_iatom forceatoms
[],
1153 const t_iparams forceparams
[],
1161 const t_mdatoms gmx_unused
* md
,
1162 t_fcdata gmx_unused
* fcd
,
1163 int gmx_unused
* global_atom_index
)
1165 int i
, m
, ai
, aj
, ak
, t1
, t2
, type
;
1167 real L1
, kA
, kB
, aA
, aB
, dr
, dr2
, va
, vtot
, a
, b
, klin
;
1168 ivec jt
, dt_ij
, dt_kj
;
1169 rvec r_ij
, r_kj
, r_ik
, dx
;
1173 for (i
= 0; (i
< nbonds
);)
1175 type
= forceatoms
[i
++];
1176 ai
= forceatoms
[i
++];
1177 aj
= forceatoms
[i
++];
1178 ak
= forceatoms
[i
++];
1180 kA
= forceparams
[type
].linangle
.klinA
;
1181 kB
= forceparams
[type
].linangle
.klinB
;
1182 klin
= L1
* kA
+ lambda
* kB
;
1184 aA
= forceparams
[type
].linangle
.aA
;
1185 aB
= forceparams
[type
].linangle
.aB
;
1186 a
= L1
* aA
+ lambda
* aB
;
1189 t1
= pbc_rvec_sub(pbc
, x
[ai
], x
[aj
], r_ij
);
1190 t2
= pbc_rvec_sub(pbc
, x
[ak
], x
[aj
], r_kj
);
1191 rvec_sub(r_ij
, r_kj
, r_ik
);
1194 for (m
= 0; (m
< DIM
); m
++)
1196 dr
= -a
* r_ij
[m
] - b
* r_kj
[m
];
1199 f_i
[m
] = a
* klin
* dr
;
1200 f_k
[m
] = b
* klin
* dr
;
1201 f_j
[m
] = -(f_i
[m
] + f_k
[m
]);
1206 va
= 0.5 * klin
* dr2
;
1207 *dvdlambda
+= 0.5 * (kB
- kA
) * dr2
+ klin
* (aB
- aA
) * iprod(dx
, r_ik
);
1211 if (computeVirial(flavor
))
1215 copy_ivec(SHIFT_IVEC(g
, aj
), jt
);
1217 ivec_sub(SHIFT_IVEC(g
, ai
), jt
, dt_ij
);
1218 ivec_sub(SHIFT_IVEC(g
, ak
), jt
, dt_kj
);
1219 t1
= IVEC2IS(dt_ij
);
1220 t2
= IVEC2IS(dt_kj
);
1222 rvec_inc(fshift
[t1
], f_i
);
1223 rvec_inc(fshift
[CENTRAL
], f_j
);
1224 rvec_inc(fshift
[t2
], f_k
);
1230 template<BondedKernelFlavor flavor
>
1231 std::enable_if_t
<flavor
!= BondedKernelFlavor::ForcesSimdWhenAvailable
|| !GMX_SIMD_HAVE_REAL
, real
>
1232 urey_bradley(int nbonds
,
1233 const t_iatom forceatoms
[],
1234 const t_iparams forceparams
[],
1242 const t_mdatoms gmx_unused
* md
,
1243 t_fcdata gmx_unused
* fcd
,
1244 int gmx_unused
* global_atom_index
)
1246 int i
, m
, ai
, aj
, ak
, t1
, t2
, type
, ki
;
1247 rvec r_ij
, r_kj
, r_ik
;
1248 real cos_theta
, cos_theta2
, theta
;
1249 real dVdt
, va
, vtot
, dr
, dr2
, vbond
, fbond
, fik
;
1250 real kthA
, th0A
, kUBA
, r13A
, kthB
, th0B
, kUBB
, r13B
;
1251 ivec jt
, dt_ij
, dt_kj
, dt_ik
;
1254 for (i
= 0; (i
< nbonds
);)
1256 type
= forceatoms
[i
++];
1257 ai
= forceatoms
[i
++];
1258 aj
= forceatoms
[i
++];
1259 ak
= forceatoms
[i
++];
1260 th0A
= forceparams
[type
].u_b
.thetaA
* DEG2RAD
;
1261 kthA
= forceparams
[type
].u_b
.kthetaA
;
1262 r13A
= forceparams
[type
].u_b
.r13A
;
1263 kUBA
= forceparams
[type
].u_b
.kUBA
;
1264 th0B
= forceparams
[type
].u_b
.thetaB
* DEG2RAD
;
1265 kthB
= forceparams
[type
].u_b
.kthetaB
;
1266 r13B
= forceparams
[type
].u_b
.r13B
;
1267 kUBB
= forceparams
[type
].u_b
.kUBB
;
1269 theta
= bond_angle(x
[ai
], x
[aj
], x
[ak
], pbc
, r_ij
, r_kj
, &cos_theta
, &t1
, &t2
); /* 41 */
1271 *dvdlambda
+= harmonic(kthA
, kthB
, th0A
, th0B
, theta
, lambda
, &va
, &dVdt
); /* 21 */
1274 ki
= pbc_rvec_sub(pbc
, x
[ai
], x
[ak
], r_ik
); /* 3 */
1275 dr2
= iprod(r_ik
, r_ik
); /* 5 */
1276 dr
= dr2
* gmx::invsqrt(dr2
); /* 10 */
1278 *dvdlambda
+= harmonic(kUBA
, kUBB
, r13A
, r13B
, dr
, lambda
, &vbond
, &fbond
); /* 19 */
1280 cos_theta2
= gmx::square(cos_theta
); /* 1 */
1288 st
= dVdt
* gmx::invsqrt(1 - cos_theta2
); /* 12 */
1289 sth
= st
* cos_theta
; /* 1 */
1290 nrkj2
= iprod(r_kj
, r_kj
); /* 5 */
1291 nrij2
= iprod(r_ij
, r_ij
);
1293 cik
= st
* gmx::invsqrt(nrkj2
* nrij2
); /* 12 */
1294 cii
= sth
/ nrij2
; /* 10 */
1295 ckk
= sth
/ nrkj2
; /* 10 */
1297 for (m
= 0; (m
< DIM
); m
++) /* 39 */
1299 f_i
[m
] = -(cik
* r_kj
[m
] - cii
* r_ij
[m
]);
1300 f_k
[m
] = -(cik
* r_ij
[m
] - ckk
* r_kj
[m
]);
1301 f_j
[m
] = -f_i
[m
] - f_k
[m
];
1306 if (computeVirial(flavor
))
1310 copy_ivec(SHIFT_IVEC(g
, aj
), jt
);
1312 ivec_sub(SHIFT_IVEC(g
, ai
), jt
, dt_ij
);
1313 ivec_sub(SHIFT_IVEC(g
, ak
), jt
, dt_kj
);
1314 t1
= IVEC2IS(dt_ij
);
1315 t2
= IVEC2IS(dt_kj
);
1317 rvec_inc(fshift
[t1
], f_i
);
1318 rvec_inc(fshift
[CENTRAL
], f_j
);
1319 rvec_inc(fshift
[t2
], f_k
);
1322 /* Time for the bond calculations */
1328 vtot
+= vbond
; /* 1*/
1329 fbond
*= gmx::invsqrt(dr2
); /* 6 */
1331 if (computeVirial(flavor
) && g
)
1333 ivec_sub(SHIFT_IVEC(g
, ai
), SHIFT_IVEC(g
, ak
), dt_ik
);
1334 ki
= IVEC2IS(dt_ik
);
1336 for (m
= 0; (m
< DIM
); m
++) /* 15 */
1338 fik
= fbond
* r_ik
[m
];
1341 if (computeVirial(flavor
))
1343 fshift
[ki
][m
] += fik
;
1344 fshift
[CENTRAL
][m
] -= fik
;
1351 #if GMX_SIMD_HAVE_REAL
1353 /* As urey_bradley, but using SIMD to calculate many potentials at once.
1354 * This routines does not calculate energies and shift forces.
1356 template<BondedKernelFlavor flavor
>
1357 std::enable_if_t
<flavor
== BondedKernelFlavor::ForcesSimdWhenAvailable
, real
>
1358 urey_bradley(int nbonds
,
1359 const t_iatom forceatoms
[],
1360 const t_iparams forceparams
[],
1363 rvec gmx_unused fshift
[],
1365 const t_graph gmx_unused
* g
,
1366 real gmx_unused lambda
,
1367 real gmx_unused
* dvdlambda
,
1368 const t_mdatoms gmx_unused
* md
,
1369 t_fcdata gmx_unused
* fcd
,
1370 int gmx_unused
* global_atom_index
)
1372 constexpr int nfa1
= 4;
1373 alignas(GMX_SIMD_ALIGNMENT
) std::int32_t ai
[GMX_SIMD_REAL_WIDTH
];
1374 alignas(GMX_SIMD_ALIGNMENT
) std::int32_t aj
[GMX_SIMD_REAL_WIDTH
];
1375 alignas(GMX_SIMD_ALIGNMENT
) std::int32_t ak
[GMX_SIMD_REAL_WIDTH
];
1376 alignas(GMX_SIMD_ALIGNMENT
) real coeff
[4 * GMX_SIMD_REAL_WIDTH
];
1377 alignas(GMX_SIMD_ALIGNMENT
) real pbc_simd
[9 * GMX_SIMD_REAL_WIDTH
];
1379 set_pbc_simd(pbc
, pbc_simd
);
1381 /* nbonds is the number of angles times nfa1, here we step GMX_SIMD_REAL_WIDTH angles */
1382 for (int i
= 0; i
< nbonds
; i
+= GMX_SIMD_REAL_WIDTH
* nfa1
)
1384 /* Collect atoms for GMX_SIMD_REAL_WIDTH angles.
1385 * iu indexes into forceatoms, we should not let iu go beyond nbonds.
1388 for (int s
= 0; s
< GMX_SIMD_REAL_WIDTH
; s
++)
1390 const int type
= forceatoms
[iu
];
1391 ai
[s
] = forceatoms
[iu
+ 1];
1392 aj
[s
] = forceatoms
[iu
+ 2];
1393 ak
[s
] = forceatoms
[iu
+ 3];
1395 /* At the end fill the arrays with the last atoms and 0 params */
1396 if (i
+ s
* nfa1
< nbonds
)
1398 coeff
[s
] = forceparams
[type
].u_b
.kthetaA
;
1399 coeff
[GMX_SIMD_REAL_WIDTH
+ s
] = forceparams
[type
].u_b
.thetaA
;
1400 coeff
[GMX_SIMD_REAL_WIDTH
* 2 + s
] = forceparams
[type
].u_b
.kUBA
;
1401 coeff
[GMX_SIMD_REAL_WIDTH
* 3 + s
] = forceparams
[type
].u_b
.r13A
;
1403 if (iu
+ nfa1
< nbonds
)
1411 coeff
[GMX_SIMD_REAL_WIDTH
+ s
] = 0;
1412 coeff
[GMX_SIMD_REAL_WIDTH
* 2 + s
] = 0;
1413 coeff
[GMX_SIMD_REAL_WIDTH
* 3 + s
] = 0;
1417 SimdReal xi_S
, yi_S
, zi_S
;
1418 SimdReal xj_S
, yj_S
, zj_S
;
1419 SimdReal xk_S
, yk_S
, zk_S
;
1421 /* Store the non PBC corrected distances packed and aligned */
1422 gatherLoadUTranspose
<3>(reinterpret_cast<const real
*>(x
), ai
, &xi_S
, &yi_S
, &zi_S
);
1423 gatherLoadUTranspose
<3>(reinterpret_cast<const real
*>(x
), aj
, &xj_S
, &yj_S
, &zj_S
);
1424 gatherLoadUTranspose
<3>(reinterpret_cast<const real
*>(x
), ak
, &xk_S
, &yk_S
, &zk_S
);
1425 SimdReal rijx_S
= xi_S
- xj_S
;
1426 SimdReal rijy_S
= yi_S
- yj_S
;
1427 SimdReal rijz_S
= zi_S
- zj_S
;
1428 SimdReal rkjx_S
= xk_S
- xj_S
;
1429 SimdReal rkjy_S
= yk_S
- yj_S
;
1430 SimdReal rkjz_S
= zk_S
- zj_S
;
1431 SimdReal rikx_S
= xi_S
- xk_S
;
1432 SimdReal riky_S
= yi_S
- yk_S
;
1433 SimdReal rikz_S
= zi_S
- zk_S
;
1435 const SimdReal ktheta_S
= load
<SimdReal
>(coeff
);
1436 const SimdReal theta0_S
= load
<SimdReal
>(coeff
+ GMX_SIMD_REAL_WIDTH
) * DEG2RAD
;
1437 const SimdReal kUB_S
= load
<SimdReal
>(coeff
+ 2 * GMX_SIMD_REAL_WIDTH
);
1438 const SimdReal r13_S
= load
<SimdReal
>(coeff
+ 3 * GMX_SIMD_REAL_WIDTH
);
1440 pbc_correct_dx_simd(&rijx_S
, &rijy_S
, &rijz_S
, pbc_simd
);
1441 pbc_correct_dx_simd(&rkjx_S
, &rkjy_S
, &rkjz_S
, pbc_simd
);
1442 pbc_correct_dx_simd(&rikx_S
, &riky_S
, &rikz_S
, pbc_simd
);
1444 const SimdReal rij_rkj_S
= iprod(rijx_S
, rijy_S
, rijz_S
, rkjx_S
, rkjy_S
, rkjz_S
);
1446 const SimdReal dr2_S
= iprod(rikx_S
, riky_S
, rikz_S
, rikx_S
, riky_S
, rikz_S
);
1448 const SimdReal nrij2_S
= norm2(rijx_S
, rijy_S
, rijz_S
);
1449 const SimdReal nrkj2_S
= norm2(rkjx_S
, rkjy_S
, rkjz_S
);
1451 const SimdReal nrij_1_S
= invsqrt(nrij2_S
);
1452 const SimdReal nrkj_1_S
= invsqrt(nrkj2_S
);
1453 const SimdReal invdr2_S
= invsqrt(dr2_S
);
1454 const SimdReal dr_S
= dr2_S
* invdr2_S
;
1456 constexpr real min_one_plus_eps
= -1.0 + 2.0 * GMX_REAL_EPS
; // Smallest number > -1
1458 /* To allow for 180 degrees, we take the max of cos and -1 + 1bit,
1459 * so we can safely get the 1/sin from 1/sqrt(1 - cos^2).
1460 * This also ensures that rounding errors would cause the argument
1461 * of simdAcos to be < -1.
1462 * Note that we do not take precautions for cos(0)=1, so the bonds
1463 * in an angle should not align at an angle of 0 degrees.
1465 const SimdReal cos_S
= max(rij_rkj_S
* nrij_1_S
* nrkj_1_S
, min_one_plus_eps
);
1467 const SimdReal theta_S
= acos(cos_S
);
1468 const SimdReal invsin_S
= invsqrt(1.0 - cos_S
* cos_S
);
1469 const SimdReal st_S
= ktheta_S
* (theta0_S
- theta_S
) * invsin_S
;
1470 const SimdReal sth_S
= st_S
* cos_S
;
1472 const SimdReal cik_S
= st_S
* nrij_1_S
* nrkj_1_S
;
1473 const SimdReal cii_S
= sth_S
* nrij_1_S
* nrij_1_S
;
1474 const SimdReal ckk_S
= sth_S
* nrkj_1_S
* nrkj_1_S
;
1476 const SimdReal sUB_S
= kUB_S
* (r13_S
- dr_S
) * invdr2_S
;
1478 const SimdReal f_ikx_S
= sUB_S
* rikx_S
;
1479 const SimdReal f_iky_S
= sUB_S
* riky_S
;
1480 const SimdReal f_ikz_S
= sUB_S
* rikz_S
;
1482 const SimdReal f_ix_S
= fnma(cik_S
, rkjx_S
, cii_S
* rijx_S
) + f_ikx_S
;
1483 const SimdReal f_iy_S
= fnma(cik_S
, rkjy_S
, cii_S
* rijy_S
) + f_iky_S
;
1484 const SimdReal f_iz_S
= fnma(cik_S
, rkjz_S
, cii_S
* rijz_S
) + f_ikz_S
;
1485 const SimdReal f_kx_S
= fnma(cik_S
, rijx_S
, ckk_S
* rkjx_S
) - f_ikx_S
;
1486 const SimdReal f_ky_S
= fnma(cik_S
, rijy_S
, ckk_S
* rkjy_S
) - f_iky_S
;
1487 const SimdReal f_kz_S
= fnma(cik_S
, rijz_S
, ckk_S
* rkjz_S
) - f_ikz_S
;
1489 transposeScatterIncrU
<4>(reinterpret_cast<real
*>(f
), ai
, f_ix_S
, f_iy_S
, f_iz_S
);
1490 transposeScatterDecrU
<4>(reinterpret_cast<real
*>(f
), aj
, f_ix_S
+ f_kx_S
, f_iy_S
+ f_ky_S
,
1492 transposeScatterIncrU
<4>(reinterpret_cast<real
*>(f
), ak
, f_kx_S
, f_ky_S
, f_kz_S
);
1498 #endif // GMX_SIMD_HAVE_REAL
1500 template<BondedKernelFlavor flavor
>
1501 real
quartic_angles(int nbonds
,
1502 const t_iatom forceatoms
[],
1503 const t_iparams forceparams
[],
1509 real gmx_unused lambda
,
1510 real gmx_unused
* dvdlambda
,
1511 const t_mdatoms gmx_unused
* md
,
1512 t_fcdata gmx_unused
* fcd
,
1513 int gmx_unused
* global_atom_index
)
1515 int i
, j
, ai
, aj
, ak
, t1
, t2
, type
;
1517 real cos_theta
, cos_theta2
, theta
, dt
, dVdt
, va
, dtp
, c
, vtot
;
1518 ivec jt
, dt_ij
, dt_kj
;
1521 for (i
= 0; (i
< nbonds
);)
1523 type
= forceatoms
[i
++];
1524 ai
= forceatoms
[i
++];
1525 aj
= forceatoms
[i
++];
1526 ak
= forceatoms
[i
++];
1528 theta
= bond_angle(x
[ai
], x
[aj
], x
[ak
], pbc
, r_ij
, r_kj
, &cos_theta
, &t1
, &t2
); /* 41 */
1530 dt
= theta
- forceparams
[type
].qangle
.theta
* DEG2RAD
; /* 2 */
1533 va
= forceparams
[type
].qangle
.c
[0];
1535 for (j
= 1; j
<= 4; j
++)
1537 c
= forceparams
[type
].qangle
.c
[j
];
1538 dVdt
-= j
* c
* dtp
;
1546 cos_theta2
= gmx::square(cos_theta
); /* 1 */
1555 st
= dVdt
* gmx::invsqrt(1 - cos_theta2
); /* 12 */
1556 sth
= st
* cos_theta
; /* 1 */
1557 nrkj2
= iprod(r_kj
, r_kj
); /* 5 */
1558 nrij2
= iprod(r_ij
, r_ij
);
1560 cik
= st
* gmx::invsqrt(nrkj2
* nrij2
); /* 12 */
1561 cii
= sth
/ nrij2
; /* 10 */
1562 ckk
= sth
/ nrkj2
; /* 10 */
1564 for (m
= 0; (m
< DIM
); m
++) /* 39 */
1566 f_i
[m
] = -(cik
* r_kj
[m
] - cii
* r_ij
[m
]);
1567 f_k
[m
] = -(cik
* r_ij
[m
] - ckk
* r_kj
[m
]);
1568 f_j
[m
] = -f_i
[m
] - f_k
[m
];
1574 if (computeVirial(flavor
))
1578 copy_ivec(SHIFT_IVEC(g
, aj
), jt
);
1580 ivec_sub(SHIFT_IVEC(g
, ai
), jt
, dt_ij
);
1581 ivec_sub(SHIFT_IVEC(g
, ak
), jt
, dt_kj
);
1582 t1
= IVEC2IS(dt_ij
);
1583 t2
= IVEC2IS(dt_kj
);
1585 rvec_inc(fshift
[t1
], f_i
);
1586 rvec_inc(fshift
[CENTRAL
], f_j
);
1587 rvec_inc(fshift
[t2
], f_k
);
1595 #if GMX_SIMD_HAVE_REAL
1597 /* As dih_angle above, but calculates 4 dihedral angles at once using SIMD,
1598 * also calculates the pre-factor required for the dihedral force update.
1599 * Note that bv and buf should be register aligned.
1601 inline void dih_angle_simd(const rvec
* x
,
1606 const real
* pbc_simd
,
1614 SimdReal
* nrkj_m2_S
,
1615 SimdReal
* nrkj_n2_S
,
1619 SimdReal xi_S
, yi_S
, zi_S
;
1620 SimdReal xj_S
, yj_S
, zj_S
;
1621 SimdReal xk_S
, yk_S
, zk_S
;
1622 SimdReal xl_S
, yl_S
, zl_S
;
1623 SimdReal rijx_S
, rijy_S
, rijz_S
;
1624 SimdReal rkjx_S
, rkjy_S
, rkjz_S
;
1625 SimdReal rklx_S
, rkly_S
, rklz_S
;
1626 SimdReal cx_S
, cy_S
, cz_S
;
1630 SimdReal iprm_S
, iprn_S
;
1631 SimdReal nrkj2_S
, nrkj_1_S
, nrkj_2_S
, nrkj_S
;
1633 SimdReal nrkj2_min_S
;
1634 SimdReal real_eps_S
;
1636 /* Used to avoid division by zero.
1637 * We take into acount that we multiply the result by real_eps_S.
1639 nrkj2_min_S
= SimdReal(GMX_REAL_MIN
/ (2 * GMX_REAL_EPS
));
1641 /* The value of the last significant bit (GMX_REAL_EPS is half of that) */
1642 real_eps_S
= SimdReal(2 * GMX_REAL_EPS
);
1644 /* Store the non PBC corrected distances packed and aligned */
1645 gatherLoadUTranspose
<3>(reinterpret_cast<const real
*>(x
), ai
, &xi_S
, &yi_S
, &zi_S
);
1646 gatherLoadUTranspose
<3>(reinterpret_cast<const real
*>(x
), aj
, &xj_S
, &yj_S
, &zj_S
);
1647 gatherLoadUTranspose
<3>(reinterpret_cast<const real
*>(x
), ak
, &xk_S
, &yk_S
, &zk_S
);
1648 gatherLoadUTranspose
<3>(reinterpret_cast<const real
*>(x
), al
, &xl_S
, &yl_S
, &zl_S
);
1649 rijx_S
= xi_S
- xj_S
;
1650 rijy_S
= yi_S
- yj_S
;
1651 rijz_S
= zi_S
- zj_S
;
1652 rkjx_S
= xk_S
- xj_S
;
1653 rkjy_S
= yk_S
- yj_S
;
1654 rkjz_S
= zk_S
- zj_S
;
1655 rklx_S
= xk_S
- xl_S
;
1656 rkly_S
= yk_S
- yl_S
;
1657 rklz_S
= zk_S
- zl_S
;
1659 pbc_correct_dx_simd(&rijx_S
, &rijy_S
, &rijz_S
, pbc_simd
);
1660 pbc_correct_dx_simd(&rkjx_S
, &rkjy_S
, &rkjz_S
, pbc_simd
);
1661 pbc_correct_dx_simd(&rklx_S
, &rkly_S
, &rklz_S
, pbc_simd
);
1663 cprod(rijx_S
, rijy_S
, rijz_S
, rkjx_S
, rkjy_S
, rkjz_S
, mx_S
, my_S
, mz_S
);
1665 cprod(rkjx_S
, rkjy_S
, rkjz_S
, rklx_S
, rkly_S
, rklz_S
, nx_S
, ny_S
, nz_S
);
1667 cprod(*mx_S
, *my_S
, *mz_S
, *nx_S
, *ny_S
, *nz_S
, &cx_S
, &cy_S
, &cz_S
);
1669 cn_S
= sqrt(norm2(cx_S
, cy_S
, cz_S
));
1671 s_S
= iprod(*mx_S
, *my_S
, *mz_S
, *nx_S
, *ny_S
, *nz_S
);
1673 /* Determine the dihedral angle, the sign might need correction */
1674 *phi_S
= atan2(cn_S
, s_S
);
1676 ipr_S
= iprod(rijx_S
, rijy_S
, rijz_S
, *nx_S
, *ny_S
, *nz_S
);
1678 iprm_S
= norm2(*mx_S
, *my_S
, *mz_S
);
1679 iprn_S
= norm2(*nx_S
, *ny_S
, *nz_S
);
1681 nrkj2_S
= norm2(rkjx_S
, rkjy_S
, rkjz_S
);
1683 /* Avoid division by zero. When zero, the result is multiplied by 0
1684 * anyhow, so the 3 max below do not affect the final result.
1686 nrkj2_S
= max(nrkj2_S
, nrkj2_min_S
);
1687 nrkj_1_S
= invsqrt(nrkj2_S
);
1688 nrkj_2_S
= nrkj_1_S
* nrkj_1_S
;
1689 nrkj_S
= nrkj2_S
* nrkj_1_S
;
1691 toler_S
= nrkj2_S
* real_eps_S
;
1693 /* Here the plain-C code uses a conditional, but we can't do that in SIMD.
1694 * So we take a max with the tolerance instead. Since we multiply with
1695 * m or n later, the max does not affect the results.
1697 iprm_S
= max(iprm_S
, toler_S
);
1698 iprn_S
= max(iprn_S
, toler_S
);
1699 *nrkj_m2_S
= nrkj_S
* inv(iprm_S
);
1700 *nrkj_n2_S
= nrkj_S
* inv(iprn_S
);
1702 /* Set sign of phi_S with the sign of ipr_S; phi_S is currently positive */
1703 *phi_S
= copysign(*phi_S
, ipr_S
);
1704 *p_S
= iprod(rijx_S
, rijy_S
, rijz_S
, rkjx_S
, rkjy_S
, rkjz_S
);
1705 *p_S
= *p_S
* nrkj_2_S
;
1707 *q_S
= iprod(rklx_S
, rkly_S
, rklz_S
, rkjx_S
, rkjy_S
, rkjz_S
);
1708 *q_S
= *q_S
* nrkj_2_S
;
1711 #endif // GMX_SIMD_HAVE_REAL
1715 template<BondedKernelFlavor flavor
>
1716 void do_dih_fup(int i
,
1736 rvec f_i
, f_j
, f_k
, f_l
;
1737 rvec uvec
, vvec
, svec
, dx_jl
;
1738 real iprm
, iprn
, nrkj
, nrkj2
, nrkj_1
, nrkj_2
;
1739 real a
, b
, p
, q
, toler
;
1740 ivec jt
, dt_ij
, dt_kj
, dt_lj
;
1742 iprm
= iprod(m
, m
); /* 5 */
1743 iprn
= iprod(n
, n
); /* 5 */
1744 nrkj2
= iprod(r_kj
, r_kj
); /* 5 */
1745 toler
= nrkj2
* GMX_REAL_EPS
;
1746 if ((iprm
> toler
) && (iprn
> toler
))
1748 nrkj_1
= gmx::invsqrt(nrkj2
); /* 10 */
1749 nrkj_2
= nrkj_1
* nrkj_1
; /* 1 */
1750 nrkj
= nrkj2
* nrkj_1
; /* 1 */
1751 a
= -ddphi
* nrkj
/ iprm
; /* 11 */
1752 svmul(a
, m
, f_i
); /* 3 */
1753 b
= ddphi
* nrkj
/ iprn
; /* 11 */
1754 svmul(b
, n
, f_l
); /* 3 */
1755 p
= iprod(r_ij
, r_kj
); /* 5 */
1756 p
*= nrkj_2
; /* 1 */
1757 q
= iprod(r_kl
, r_kj
); /* 5 */
1758 q
*= nrkj_2
; /* 1 */
1759 svmul(p
, f_i
, uvec
); /* 3 */
1760 svmul(q
, f_l
, vvec
); /* 3 */
1761 rvec_sub(uvec
, vvec
, svec
); /* 3 */
1762 rvec_sub(f_i
, svec
, f_j
); /* 3 */
1763 rvec_add(f_l
, svec
, f_k
); /* 3 */
1764 rvec_inc(f
[i
], f_i
); /* 3 */
1765 rvec_dec(f
[j
], f_j
); /* 3 */
1766 rvec_dec(f
[k
], f_k
); /* 3 */
1767 rvec_inc(f
[l
], f_l
); /* 3 */
1769 if (computeVirial(flavor
))
1773 copy_ivec(SHIFT_IVEC(g
, j
), jt
);
1774 ivec_sub(SHIFT_IVEC(g
, i
), jt
, dt_ij
);
1775 ivec_sub(SHIFT_IVEC(g
, k
), jt
, dt_kj
);
1776 ivec_sub(SHIFT_IVEC(g
, l
), jt
, dt_lj
);
1777 t1
= IVEC2IS(dt_ij
);
1778 t2
= IVEC2IS(dt_kj
);
1779 t3
= IVEC2IS(dt_lj
);
1783 t3
= pbc_rvec_sub(pbc
, x
[l
], x
[j
], dx_jl
);
1790 rvec_inc(fshift
[t1
], f_i
);
1791 rvec_dec(fshift
[CENTRAL
], f_j
);
1792 rvec_dec(fshift
[t2
], f_k
);
1793 rvec_inc(fshift
[t3
], f_l
);
1802 #if GMX_SIMD_HAVE_REAL
1803 /* As do_dih_fup_noshiftf above, but with SIMD and pre-calculated pre-factors */
1804 inline void gmx_simdcall
do_dih_fup_noshiftf_simd(const int* ai
,
1818 SimdReal sx
= p
* f_i_x
+ q
* mf_l_x
;
1819 SimdReal sy
= p
* f_i_y
+ q
* mf_l_y
;
1820 SimdReal sz
= p
* f_i_z
+ q
* mf_l_z
;
1821 SimdReal f_j_x
= f_i_x
- sx
;
1822 SimdReal f_j_y
= f_i_y
- sy
;
1823 SimdReal f_j_z
= f_i_z
- sz
;
1824 SimdReal f_k_x
= mf_l_x
- sx
;
1825 SimdReal f_k_y
= mf_l_y
- sy
;
1826 SimdReal f_k_z
= mf_l_z
- sz
;
1827 transposeScatterIncrU
<4>(reinterpret_cast<real
*>(f
), ai
, f_i_x
, f_i_y
, f_i_z
);
1828 transposeScatterDecrU
<4>(reinterpret_cast<real
*>(f
), aj
, f_j_x
, f_j_y
, f_j_z
);
1829 transposeScatterIncrU
<4>(reinterpret_cast<real
*>(f
), ak
, f_k_x
, f_k_y
, f_k_z
);
1830 transposeScatterDecrU
<4>(reinterpret_cast<real
*>(f
), al
, mf_l_x
, mf_l_y
, mf_l_z
);
1832 #endif // GMX_SIMD_HAVE_REAL
1834 /*! \brief Computes and returns the proper dihedral force
1836 * With the appropriate kernel flavor, also computes and accumulates
1837 * the energy and dV/dlambda.
1839 template<BondedKernelFlavor flavor
>
1840 real
dopdihs(real cpA
, real cpB
, real phiA
, real phiB
, int mult
, real phi
, real lambda
, real
* V
, real
* dvdlambda
)
1842 const real L1
= 1.0 - lambda
;
1843 const real ph0
= (L1
* phiA
+ lambda
* phiB
) * DEG2RAD
;
1844 const real dph0
= (phiB
- phiA
) * DEG2RAD
;
1845 const real cp
= L1
* cpA
+ lambda
* cpB
;
1847 const real mdphi
= mult
* phi
- ph0
;
1848 const real sdphi
= std::sin(mdphi
);
1849 const real ddphi
= -cp
* mult
* sdphi
;
1850 if (computeEnergy(flavor
))
1852 const real v1
= 1 + std::cos(mdphi
);
1855 *dvdlambda
+= (cpB
- cpA
) * v1
+ cp
* dph0
* sdphi
;
1859 /* That was 40 flops */
1862 /*! \brief Similar to \p dopdihs(), except for a minus sign and a different treatment of mult/phi0 */
1863 real
dopdihs_min(real cpA
, real cpB
, real phiA
, real phiB
, int mult
, real phi
, real lambda
, real
* V
, real
* F
)
1865 real v
, dvdlambda
, mdphi
, v1
, sdphi
, ddphi
;
1866 real L1
= 1.0 - lambda
;
1867 real ph0
= (L1
* phiA
+ lambda
* phiB
) * DEG2RAD
;
1868 real dph0
= (phiB
- phiA
) * DEG2RAD
;
1869 real cp
= L1
* cpA
+ lambda
* cpB
;
1871 mdphi
= mult
* (phi
- ph0
);
1872 sdphi
= std::sin(mdphi
);
1873 ddphi
= cp
* mult
* sdphi
;
1874 v1
= 1.0 - std::cos(mdphi
);
1877 dvdlambda
= (cpB
- cpA
) * v1
+ cp
* dph0
* sdphi
;
1884 /* That was 40 flops */
1887 template<BondedKernelFlavor flavor
>
1888 std::enable_if_t
<flavor
!= BondedKernelFlavor::ForcesSimdWhenAvailable
|| !GMX_SIMD_HAVE_REAL
, real
>
1890 const t_iatom forceatoms
[],
1891 const t_iparams forceparams
[],
1899 const t_mdatoms gmx_unused
* md
,
1900 t_fcdata gmx_unused
* fcd
,
1901 int gmx_unused
* global_atom_index
)
1904 rvec r_ij
, r_kj
, r_kl
, m
, n
;
1908 for (int i
= 0; i
< nbonds
;)
1910 const int ai
= forceatoms
[i
+ 1];
1911 const int aj
= forceatoms
[i
+ 2];
1912 const int ak
= forceatoms
[i
+ 3];
1913 const int al
= forceatoms
[i
+ 4];
1915 const real phi
= dih_angle(x
[ai
], x
[aj
], x
[ak
], x
[al
], pbc
, r_ij
, r_kj
, r_kl
, m
, n
, &t1
,
1918 /* Loop over dihedrals working on the same atoms,
1919 * so we avoid recalculating angles and distributing forces.
1924 const int type
= forceatoms
[i
];
1925 ddphi_tot
+= dopdihs
<flavor
>(forceparams
[type
].pdihs
.cpA
, forceparams
[type
].pdihs
.cpB
,
1926 forceparams
[type
].pdihs
.phiA
, forceparams
[type
].pdihs
.phiB
,
1927 forceparams
[type
].pdihs
.mult
, phi
, lambda
, &vtot
, dvdlambda
);
1930 } while (i
< nbonds
&& forceatoms
[i
+ 1] == ai
&& forceatoms
[i
+ 2] == aj
1931 && forceatoms
[i
+ 3] == ak
&& forceatoms
[i
+ 4] == al
);
1933 do_dih_fup
<flavor
>(ai
, aj
, ak
, al
, ddphi_tot
, r_ij
, r_kj
, r_kl
, m
, n
, f
, fshift
, pbc
, g
, x
,
1934 t1
, t2
, t3
); /* 112 */
1940 #if GMX_SIMD_HAVE_REAL
1942 /* As pdihs above, but using SIMD to calculate multiple dihedrals at once */
1943 template<BondedKernelFlavor flavor
>
1944 std::enable_if_t
<flavor
== BondedKernelFlavor::ForcesSimdWhenAvailable
, real
>
1946 const t_iatom forceatoms
[],
1947 const t_iparams forceparams
[],
1950 rvec gmx_unused fshift
[],
1952 const t_graph gmx_unused
* g
,
1953 real gmx_unused lambda
,
1954 real gmx_unused
* dvdlambda
,
1955 const t_mdatoms gmx_unused
* md
,
1956 t_fcdata gmx_unused
* fcd
,
1957 int gmx_unused
* global_atom_index
)
1962 alignas(GMX_SIMD_ALIGNMENT
) std::int32_t ai
[GMX_SIMD_REAL_WIDTH
];
1963 alignas(GMX_SIMD_ALIGNMENT
) std::int32_t aj
[GMX_SIMD_REAL_WIDTH
];
1964 alignas(GMX_SIMD_ALIGNMENT
) std::int32_t ak
[GMX_SIMD_REAL_WIDTH
];
1965 alignas(GMX_SIMD_ALIGNMENT
) std::int32_t al
[GMX_SIMD_REAL_WIDTH
];
1966 alignas(GMX_SIMD_ALIGNMENT
) real buf
[3 * GMX_SIMD_REAL_WIDTH
];
1967 real
* cp
, *phi0
, *mult
;
1968 SimdReal
deg2rad_S(DEG2RAD
);
1970 SimdReal phi0_S
, phi_S
;
1971 SimdReal mx_S
, my_S
, mz_S
;
1972 SimdReal nx_S
, ny_S
, nz_S
;
1973 SimdReal nrkj_m2_S
, nrkj_n2_S
;
1974 SimdReal cp_S
, mdphi_S
, mult_S
;
1975 SimdReal sin_S
, cos_S
;
1977 SimdReal sf_i_S
, msf_l_S
;
1978 alignas(GMX_SIMD_ALIGNMENT
) real pbc_simd
[9 * GMX_SIMD_REAL_WIDTH
];
1980 /* Extract aligned pointer for parameters and variables */
1981 cp
= buf
+ 0 * GMX_SIMD_REAL_WIDTH
;
1982 phi0
= buf
+ 1 * GMX_SIMD_REAL_WIDTH
;
1983 mult
= buf
+ 2 * GMX_SIMD_REAL_WIDTH
;
1985 set_pbc_simd(pbc
, pbc_simd
);
1987 /* nbonds is the number of dihedrals times nfa1, here we step GMX_SIMD_REAL_WIDTH dihs */
1988 for (i
= 0; (i
< nbonds
); i
+= GMX_SIMD_REAL_WIDTH
* nfa1
)
1990 /* Collect atoms quadruplets for GMX_SIMD_REAL_WIDTH dihedrals.
1991 * iu indexes into forceatoms, we should not let iu go beyond nbonds.
1994 for (s
= 0; s
< GMX_SIMD_REAL_WIDTH
; s
++)
1996 type
= forceatoms
[iu
];
1997 ai
[s
] = forceatoms
[iu
+ 1];
1998 aj
[s
] = forceatoms
[iu
+ 2];
1999 ak
[s
] = forceatoms
[iu
+ 3];
2000 al
[s
] = forceatoms
[iu
+ 4];
2002 /* At the end fill the arrays with the last atoms and 0 params */
2003 if (i
+ s
* nfa1
< nbonds
)
2005 cp
[s
] = forceparams
[type
].pdihs
.cpA
;
2006 phi0
[s
] = forceparams
[type
].pdihs
.phiA
;
2007 mult
[s
] = forceparams
[type
].pdihs
.mult
;
2009 if (iu
+ nfa1
< nbonds
)
2022 /* Caclulate GMX_SIMD_REAL_WIDTH dihedral angles at once */
2023 dih_angle_simd(x
, ai
, aj
, ak
, al
, pbc_simd
, &phi_S
, &mx_S
, &my_S
, &mz_S
, &nx_S
, &ny_S
,
2024 &nz_S
, &nrkj_m2_S
, &nrkj_n2_S
, &p_S
, &q_S
);
2026 cp_S
= load
<SimdReal
>(cp
);
2027 phi0_S
= load
<SimdReal
>(phi0
) * deg2rad_S
;
2028 mult_S
= load
<SimdReal
>(mult
);
2030 mdphi_S
= fms(mult_S
, phi_S
, phi0_S
);
2032 /* Calculate GMX_SIMD_REAL_WIDTH sines at once */
2033 sincos(mdphi_S
, &sin_S
, &cos_S
);
2034 mddphi_S
= cp_S
* mult_S
* sin_S
;
2035 sf_i_S
= mddphi_S
* nrkj_m2_S
;
2036 msf_l_S
= mddphi_S
* nrkj_n2_S
;
2038 /* After this m?_S will contain f[i] */
2039 mx_S
= sf_i_S
* mx_S
;
2040 my_S
= sf_i_S
* my_S
;
2041 mz_S
= sf_i_S
* mz_S
;
2043 /* After this m?_S will contain -f[l] */
2044 nx_S
= msf_l_S
* nx_S
;
2045 ny_S
= msf_l_S
* ny_S
;
2046 nz_S
= msf_l_S
* nz_S
;
2048 do_dih_fup_noshiftf_simd(ai
, aj
, ak
, al
, p_S
, q_S
, mx_S
, my_S
, mz_S
, nx_S
, ny_S
, nz_S
, f
);
2054 /* This is mostly a copy of the SIMD flavor of pdihs above, but with using
2055 * the RB potential instead of a harmonic potential.
2056 * This function can replace rbdihs() when no energy and virial are needed.
2058 template<BondedKernelFlavor flavor
>
2059 std::enable_if_t
<flavor
== BondedKernelFlavor::ForcesSimdWhenAvailable
, real
>
2061 const t_iatom forceatoms
[],
2062 const t_iparams forceparams
[],
2065 rvec gmx_unused fshift
[],
2067 const t_graph gmx_unused
* g
,
2068 real gmx_unused lambda
,
2069 real gmx_unused
* dvdlambda
,
2070 const t_mdatoms gmx_unused
* md
,
2071 t_fcdata gmx_unused
* fcd
,
2072 int gmx_unused
* global_atom_index
)
2077 alignas(GMX_SIMD_ALIGNMENT
) std::int32_t ai
[GMX_SIMD_REAL_WIDTH
];
2078 alignas(GMX_SIMD_ALIGNMENT
) std::int32_t aj
[GMX_SIMD_REAL_WIDTH
];
2079 alignas(GMX_SIMD_ALIGNMENT
) std::int32_t ak
[GMX_SIMD_REAL_WIDTH
];
2080 alignas(GMX_SIMD_ALIGNMENT
) std::int32_t al
[GMX_SIMD_REAL_WIDTH
];
2081 alignas(GMX_SIMD_ALIGNMENT
) real parm
[NR_RBDIHS
* GMX_SIMD_REAL_WIDTH
];
2085 SimdReal ddphi_S
, cosfac_S
;
2086 SimdReal mx_S
, my_S
, mz_S
;
2087 SimdReal nx_S
, ny_S
, nz_S
;
2088 SimdReal nrkj_m2_S
, nrkj_n2_S
;
2089 SimdReal parm_S
, c_S
;
2090 SimdReal sin_S
, cos_S
;
2091 SimdReal sf_i_S
, msf_l_S
;
2092 alignas(GMX_SIMD_ALIGNMENT
) real pbc_simd
[9 * GMX_SIMD_REAL_WIDTH
];
2094 SimdReal
pi_S(M_PI
);
2095 SimdReal
one_S(1.0);
2097 set_pbc_simd(pbc
, pbc_simd
);
2099 /* nbonds is the number of dihedrals times nfa1, here we step GMX_SIMD_REAL_WIDTH dihs */
2100 for (i
= 0; (i
< nbonds
); i
+= GMX_SIMD_REAL_WIDTH
* nfa1
)
2102 /* Collect atoms quadruplets for GMX_SIMD_REAL_WIDTH dihedrals.
2103 * iu indexes into forceatoms, we should not let iu go beyond nbonds.
2106 for (s
= 0; s
< GMX_SIMD_REAL_WIDTH
; s
++)
2108 type
= forceatoms
[iu
];
2109 ai
[s
] = forceatoms
[iu
+ 1];
2110 aj
[s
] = forceatoms
[iu
+ 2];
2111 ak
[s
] = forceatoms
[iu
+ 3];
2112 al
[s
] = forceatoms
[iu
+ 4];
2114 /* At the end fill the arrays with the last atoms and 0 params */
2115 if (i
+ s
* nfa1
< nbonds
)
2117 /* We don't need the first parameter, since that's a constant
2118 * which only affects the energies, not the forces.
2120 for (j
= 1; j
< NR_RBDIHS
; j
++)
2122 parm
[j
* GMX_SIMD_REAL_WIDTH
+ s
] = forceparams
[type
].rbdihs
.rbcA
[j
];
2125 if (iu
+ nfa1
< nbonds
)
2132 for (j
= 1; j
< NR_RBDIHS
; j
++)
2134 parm
[j
* GMX_SIMD_REAL_WIDTH
+ s
] = 0;
2139 /* Caclulate GMX_SIMD_REAL_WIDTH dihedral angles at once */
2140 dih_angle_simd(x
, ai
, aj
, ak
, al
, pbc_simd
, &phi_S
, &mx_S
, &my_S
, &mz_S
, &nx_S
, &ny_S
,
2141 &nz_S
, &nrkj_m2_S
, &nrkj_n2_S
, &p_S
, &q_S
);
2143 /* Change to polymer convention */
2144 phi_S
= phi_S
- pi_S
;
2146 sincos(phi_S
, &sin_S
, &cos_S
);
2148 ddphi_S
= setZero();
2151 for (j
= 1; j
< NR_RBDIHS
; j
++)
2153 parm_S
= load
<SimdReal
>(parm
+ j
* GMX_SIMD_REAL_WIDTH
);
2154 ddphi_S
= fma(c_S
* parm_S
, cosfac_S
, ddphi_S
);
2155 cosfac_S
= cosfac_S
* cos_S
;
2159 /* Note that here we do not use the minus sign which is present
2160 * in the normal RB code. This is corrected for through (m)sf below.
2162 ddphi_S
= ddphi_S
* sin_S
;
2164 sf_i_S
= ddphi_S
* nrkj_m2_S
;
2165 msf_l_S
= ddphi_S
* nrkj_n2_S
;
2167 /* After this m?_S will contain f[i] */
2168 mx_S
= sf_i_S
* mx_S
;
2169 my_S
= sf_i_S
* my_S
;
2170 mz_S
= sf_i_S
* mz_S
;
2172 /* After this m?_S will contain -f[l] */
2173 nx_S
= msf_l_S
* nx_S
;
2174 ny_S
= msf_l_S
* ny_S
;
2175 nz_S
= msf_l_S
* nz_S
;
2177 do_dih_fup_noshiftf_simd(ai
, aj
, ak
, al
, p_S
, q_S
, mx_S
, my_S
, mz_S
, nx_S
, ny_S
, nz_S
, f
);
2183 #endif // GMX_SIMD_HAVE_REAL
2186 template<BondedKernelFlavor flavor
>
2187 real
idihs(int nbonds
,
2188 const t_iatom forceatoms
[],
2189 const t_iparams forceparams
[],
2197 const t_mdatoms gmx_unused
* md
,
2198 t_fcdata gmx_unused
* fcd
,
2199 int gmx_unused
* global_atom_index
)
2201 int i
, type
, ai
, aj
, ak
, al
;
2203 real phi
, phi0
, dphi0
, ddphi
, vtot
;
2204 rvec r_ij
, r_kj
, r_kl
, m
, n
;
2205 real L1
, kk
, dp
, dp2
, kA
, kB
, pA
, pB
, dvdl_term
;
2210 for (i
= 0; (i
< nbonds
);)
2212 type
= forceatoms
[i
++];
2213 ai
= forceatoms
[i
++];
2214 aj
= forceatoms
[i
++];
2215 ak
= forceatoms
[i
++];
2216 al
= forceatoms
[i
++];
2218 phi
= dih_angle(x
[ai
], x
[aj
], x
[ak
], x
[al
], pbc
, r_ij
, r_kj
, r_kl
, m
, n
, &t1
, &t2
, &t3
); /* 84 */
2220 /* phi can jump if phi0 is close to Pi/-Pi, which will cause huge
2221 * force changes if we just apply a normal harmonic.
2222 * Instead, we first calculate phi-phi0 and take it modulo (-Pi,Pi).
2223 * This means we will never have the periodicity problem, unless
2224 * the dihedral is Pi away from phiO, which is very unlikely due to
2227 kA
= forceparams
[type
].harmonic
.krA
;
2228 kB
= forceparams
[type
].harmonic
.krB
;
2229 pA
= forceparams
[type
].harmonic
.rA
;
2230 pB
= forceparams
[type
].harmonic
.rB
;
2232 kk
= L1
* kA
+ lambda
* kB
;
2233 phi0
= (L1
* pA
+ lambda
* pB
) * DEG2RAD
;
2234 dphi0
= (pB
- pA
) * DEG2RAD
;
2238 make_dp_periodic(&dp
);
2242 vtot
+= 0.5 * kk
* dp2
;
2245 dvdl_term
+= 0.5 * (kB
- kA
) * dp2
- kk
* dphi0
* dp
;
2247 do_dih_fup
<flavor
>(ai
, aj
, ak
, al
, -ddphi
, r_ij
, r_kj
, r_kl
, m
, n
, f
, fshift
, pbc
, g
, x
, t1
,
2252 *dvdlambda
+= dvdl_term
;
2256 /*! \brief Computes angle restraints of two different types */
2257 template<BondedKernelFlavor flavor
>
2258 real
low_angres(int nbonds
,
2259 const t_iatom forceatoms
[],
2260 const t_iparams forceparams
[],
2270 int i
, m
, type
, ai
, aj
, ak
, al
;
2272 real phi
, cos_phi
, cos_phi2
, vid
, vtot
, dVdphi
;
2273 rvec r_ij
, r_kl
, f_i
, f_k
= { 0, 0, 0 };
2274 real st
, sth
, nrij2
, nrkl2
, c
, cij
, ckl
;
2277 t2
= 0; /* avoid warning with gcc-3.3. It is never used uninitialized */
2280 ak
= al
= 0; /* to avoid warnings */
2281 for (i
= 0; i
< nbonds
;)
2283 type
= forceatoms
[i
++];
2284 ai
= forceatoms
[i
++];
2285 aj
= forceatoms
[i
++];
2286 t1
= pbc_rvec_sub(pbc
, x
[aj
], x
[ai
], r_ij
); /* 3 */
2289 ak
= forceatoms
[i
++];
2290 al
= forceatoms
[i
++];
2291 t2
= pbc_rvec_sub(pbc
, x
[al
], x
[ak
], r_kl
); /* 3 */
2300 cos_phi
= cos_angle(r_ij
, r_kl
); /* 25 */
2301 phi
= std::acos(cos_phi
); /* 10 */
2303 *dvdlambda
+= dopdihs_min(forceparams
[type
].pdihs
.cpA
, forceparams
[type
].pdihs
.cpB
,
2304 forceparams
[type
].pdihs
.phiA
, forceparams
[type
].pdihs
.phiB
,
2305 forceparams
[type
].pdihs
.mult
, phi
, lambda
, &vid
, &dVdphi
); /* 40 */
2309 cos_phi2
= gmx::square(cos_phi
); /* 1 */
2312 st
= -dVdphi
* gmx::invsqrt(1 - cos_phi2
); /* 12 */
2313 sth
= st
* cos_phi
; /* 1 */
2314 nrij2
= iprod(r_ij
, r_ij
); /* 5 */
2315 nrkl2
= iprod(r_kl
, r_kl
); /* 5 */
2317 c
= st
* gmx::invsqrt(nrij2
* nrkl2
); /* 11 */
2318 cij
= sth
/ nrij2
; /* 10 */
2319 ckl
= sth
/ nrkl2
; /* 10 */
2321 for (m
= 0; m
< DIM
; m
++) /* 18+18 */
2323 f_i
[m
] = (c
* r_kl
[m
] - cij
* r_ij
[m
]);
2328 f_k
[m
] = (c
* r_ij
[m
] - ckl
* r_kl
[m
]);
2334 if (computeVirial(flavor
))
2338 ivec_sub(SHIFT_IVEC(g
, ai
), SHIFT_IVEC(g
, aj
), dt
);
2341 rvec_inc(fshift
[t1
], f_i
);
2342 rvec_dec(fshift
[CENTRAL
], f_i
);
2347 ivec_sub(SHIFT_IVEC(g
, ak
), SHIFT_IVEC(g
, al
), dt
);
2350 rvec_inc(fshift
[t2
], f_k
);
2351 rvec_dec(fshift
[CENTRAL
], f_k
);
2357 return vtot
; /* 184 / 157 (bZAxis) total */
2360 template<BondedKernelFlavor flavor
>
2361 real
angres(int nbonds
,
2362 const t_iatom forceatoms
[],
2363 const t_iparams forceparams
[],
2371 const t_mdatoms gmx_unused
* md
,
2372 t_fcdata gmx_unused
* fcd
,
2373 int gmx_unused
* global_atom_index
)
2375 return low_angres
<flavor
>(nbonds
, forceatoms
, forceparams
, x
, f
, fshift
, pbc
, g
, lambda
,
2379 template<BondedKernelFlavor flavor
>
2380 real
angresz(int nbonds
,
2381 const t_iatom forceatoms
[],
2382 const t_iparams forceparams
[],
2390 const t_mdatoms gmx_unused
* md
,
2391 t_fcdata gmx_unused
* fcd
,
2392 int gmx_unused
* global_atom_index
)
2394 return low_angres
<flavor
>(nbonds
, forceatoms
, forceparams
, x
, f
, fshift
, pbc
, g
, lambda
,
2398 template<BondedKernelFlavor flavor
>
2399 real
dihres(int nbonds
,
2400 const t_iatom forceatoms
[],
2401 const t_iparams forceparams
[],
2409 const t_mdatoms gmx_unused
* md
,
2410 t_fcdata gmx_unused
* fcd
,
2411 int gmx_unused
* global_atom_index
)
2414 int ai
, aj
, ak
, al
, i
, type
, t1
, t2
, t3
;
2415 real phi0A
, phi0B
, dphiA
, dphiB
, kfacA
, kfacB
, phi0
, dphi
, kfac
;
2416 real phi
, ddphi
, ddp
, ddp2
, dp
, d2r
, L1
;
2417 rvec r_ij
, r_kj
, r_kl
, m
, n
;
2423 for (i
= 0; (i
< nbonds
);)
2425 type
= forceatoms
[i
++];
2426 ai
= forceatoms
[i
++];
2427 aj
= forceatoms
[i
++];
2428 ak
= forceatoms
[i
++];
2429 al
= forceatoms
[i
++];
2431 phi0A
= forceparams
[type
].dihres
.phiA
* d2r
;
2432 dphiA
= forceparams
[type
].dihres
.dphiA
* d2r
;
2433 kfacA
= forceparams
[type
].dihres
.kfacA
;
2435 phi0B
= forceparams
[type
].dihres
.phiB
* d2r
;
2436 dphiB
= forceparams
[type
].dihres
.dphiB
* d2r
;
2437 kfacB
= forceparams
[type
].dihres
.kfacB
;
2439 phi0
= L1
* phi0A
+ lambda
* phi0B
;
2440 dphi
= L1
* dphiA
+ lambda
* dphiB
;
2441 kfac
= L1
* kfacA
+ lambda
* kfacB
;
2443 phi
= dih_angle(x
[ai
], x
[aj
], x
[ak
], x
[al
], pbc
, r_ij
, r_kj
, r_kl
, m
, n
, &t1
, &t2
, &t3
);
2446 /* phi can jump if phi0 is close to Pi/-Pi, which will cause huge
2447 * force changes if we just apply a normal harmonic.
2448 * Instead, we first calculate phi-phi0 and take it modulo (-Pi,Pi).
2449 * This means we will never have the periodicity problem, unless
2450 * the dihedral is Pi away from phiO, which is very unlikely due to
2454 make_dp_periodic(&dp
);
2460 else if (dp
< -dphi
)
2472 vtot
+= 0.5 * kfac
* ddp2
;
2475 *dvdlambda
+= 0.5 * (kfacB
- kfacA
) * ddp2
;
2476 /* lambda dependence from changing restraint distances */
2479 *dvdlambda
-= kfac
* ddp
* ((dphiB
- dphiA
) + (phi0B
- phi0A
));
2483 *dvdlambda
+= kfac
* ddp
* ((dphiB
- dphiA
) - (phi0B
- phi0A
));
2485 do_dih_fup
<flavor
>(ai
, aj
, ak
, al
, ddphi
, r_ij
, r_kj
, r_kl
, m
, n
, f
, fshift
, pbc
, g
, x
,
2486 t1
, t2
, t3
); /* 112 */
2493 real
unimplemented(int gmx_unused nbonds
,
2494 const t_iatom gmx_unused forceatoms
[],
2495 const t_iparams gmx_unused forceparams
[],
2496 const rvec gmx_unused x
[],
2497 rvec4 gmx_unused f
[],
2498 rvec gmx_unused fshift
[],
2499 const t_pbc gmx_unused
* pbc
,
2500 const t_graph gmx_unused
* g
,
2501 real gmx_unused lambda
,
2502 real gmx_unused
* dvdlambda
,
2503 const t_mdatoms gmx_unused
* md
,
2504 t_fcdata gmx_unused
* fcd
,
2505 int gmx_unused
* global_atom_index
)
2507 gmx_impl("*** you are using a not implemented function");
2510 template<BondedKernelFlavor flavor
>
2511 real
restrangles(int nbonds
,
2512 const t_iatom forceatoms
[],
2513 const t_iparams forceparams
[],
2519 real gmx_unused lambda
,
2520 real gmx_unused
* dvdlambda
,
2521 const t_mdatoms gmx_unused
* md
,
2522 t_fcdata gmx_unused
* fcd
,
2523 int gmx_unused
* global_atom_index
)
2525 int i
, d
, ai
, aj
, ak
, type
, m
;
2528 ivec jt
, dt_ij
, dt_kj
;
2530 double prefactor
, ratio_ante
, ratio_post
;
2531 rvec delta_ante
, delta_post
, vec_temp
;
2534 for (i
= 0; (i
< nbonds
);)
2536 type
= forceatoms
[i
++];
2537 ai
= forceatoms
[i
++];
2538 aj
= forceatoms
[i
++];
2539 ak
= forceatoms
[i
++];
2541 t1
= pbc_rvec_sub(pbc
, x
[ai
], x
[aj
], vec_temp
);
2542 pbc_rvec_sub(pbc
, x
[aj
], x
[ai
], delta_ante
);
2543 t2
= pbc_rvec_sub(pbc
, x
[ak
], x
[aj
], delta_post
);
2546 /* This function computes factors needed for restricted angle potential.
2547 * The restricted angle potential is used in coarse-grained simulations to avoid singularities
2548 * when three particles align and the dihedral angle and dihedral potential
2549 * cannot be calculated. This potential is calculated using the formula:
2550 real restrangles(int nbonds,
2551 const t_iatom forceatoms[],const t_iparams forceparams[],
2552 const rvec x[],rvec4 f[],rvec fshift[],
2553 const t_pbc *pbc,const t_graph *g,
2554 real gmx_unused lambda,real gmx_unused *dvdlambda,
2555 const t_mdatoms gmx_unused *md,t_fcdata gmx_unused *fcd,
2556 int gmx_unused *global_atom_index)
2558 int i, d, ai, aj, ak, type, m;
2562 ivec jt,dt_ij,dt_kj;
2564 real prefactor, ratio_ante, ratio_post;
2565 rvec delta_ante, delta_post, vec_temp;
2568 for(i=0; (i<nbonds); )
2570 type = forceatoms[i++];
2571 ai = forceatoms[i++];
2572 aj = forceatoms[i++];
2573 ak = forceatoms[i++];
2575 * \f[V_{\rm ReB}(\theta_i) = \frac{1}{2} k_{\theta} \frac{(\cos\theta_i - \cos\theta_0)^2}
2576 * {\sin^2\theta_i}\f] ({eq:ReB} and ref \cite{MonicaGoga2013} from the manual).
2577 * For more explanations see comments file "restcbt.h". */
2579 compute_factors_restangles(type
, forceparams
, delta_ante
, delta_post
, &prefactor
,
2580 &ratio_ante
, &ratio_post
, &v
);
2582 /* Forces are computed per component */
2583 for (d
= 0; d
< DIM
; d
++)
2585 f_i
[d
] = prefactor
* (ratio_ante
* delta_ante
[d
] - delta_post
[d
]);
2587 * ((ratio_post
+ 1.0) * delta_post
[d
] - (ratio_ante
+ 1.0) * delta_ante
[d
]);
2588 f_k
[d
] = prefactor
* (delta_ante
[d
] - ratio_post
* delta_post
[d
]);
2591 /* Computation of potential energy */
2597 for (m
= 0; (m
< DIM
); m
++)
2604 if (computeVirial(flavor
))
2608 copy_ivec(SHIFT_IVEC(g
, aj
), jt
);
2609 ivec_sub(SHIFT_IVEC(g
, ai
), jt
, dt_ij
);
2610 ivec_sub(SHIFT_IVEC(g
, ak
), jt
, dt_kj
);
2611 t1
= IVEC2IS(dt_ij
);
2612 t2
= IVEC2IS(dt_kj
);
2615 rvec_inc(fshift
[t1
], f_i
);
2616 rvec_inc(fshift
[CENTRAL
], f_j
);
2617 rvec_inc(fshift
[t2
], f_k
);
2624 template<BondedKernelFlavor flavor
>
2625 real
restrdihs(int nbonds
,
2626 const t_iatom forceatoms
[],
2627 const t_iparams forceparams
[],
2633 real gmx_unused lambda
,
2634 real gmx_unused
* dvlambda
,
2635 const t_mdatoms gmx_unused
* md
,
2636 t_fcdata gmx_unused
* fcd
,
2637 int gmx_unused
* global_atom_index
)
2639 int i
, d
, type
, ai
, aj
, ak
, al
;
2640 rvec f_i
, f_j
, f_k
, f_l
;
2642 ivec jt
, dt_ij
, dt_kj
, dt_lj
;
2645 rvec delta_ante
, delta_crnt
, delta_post
, vec_temp
;
2646 real factor_phi_ai_ante
, factor_phi_ai_crnt
, factor_phi_ai_post
;
2647 real factor_phi_aj_ante
, factor_phi_aj_crnt
, factor_phi_aj_post
;
2648 real factor_phi_ak_ante
, factor_phi_ak_crnt
, factor_phi_ak_post
;
2649 real factor_phi_al_ante
, factor_phi_al_crnt
, factor_phi_al_post
;
2654 for (i
= 0; (i
< nbonds
);)
2656 type
= forceatoms
[i
++];
2657 ai
= forceatoms
[i
++];
2658 aj
= forceatoms
[i
++];
2659 ak
= forceatoms
[i
++];
2660 al
= forceatoms
[i
++];
2662 t1
= pbc_rvec_sub(pbc
, x
[ai
], x
[aj
], vec_temp
);
2663 pbc_rvec_sub(pbc
, x
[aj
], x
[ai
], delta_ante
);
2664 t2
= pbc_rvec_sub(pbc
, x
[ak
], x
[aj
], delta_crnt
);
2665 pbc_rvec_sub(pbc
, x
[ak
], x
[al
], vec_temp
);
2666 pbc_rvec_sub(pbc
, x
[al
], x
[ak
], delta_post
);
2668 /* This function computes factors needed for restricted angle potential.
2669 * The restricted angle potential is used in coarse-grained simulations to avoid singularities
2670 * when three particles align and the dihedral angle and dihedral potential cannot be calculated.
2671 * This potential is calculated using the formula:
2672 * \f[V_{\rm ReB}(\theta_i) = \frac{1}{2} k_{\theta}
2673 * \frac{(\cos\theta_i - \cos\theta_0)^2}{\sin^2\theta_i}\f]
2674 * ({eq:ReB} and ref \cite{MonicaGoga2013} from the manual).
2675 * For more explanations see comments file "restcbt.h" */
2677 compute_factors_restrdihs(
2678 type
, forceparams
, delta_ante
, delta_crnt
, delta_post
, &factor_phi_ai_ante
,
2679 &factor_phi_ai_crnt
, &factor_phi_ai_post
, &factor_phi_aj_ante
, &factor_phi_aj_crnt
,
2680 &factor_phi_aj_post
, &factor_phi_ak_ante
, &factor_phi_ak_crnt
, &factor_phi_ak_post
,
2681 &factor_phi_al_ante
, &factor_phi_al_crnt
, &factor_phi_al_post
, &prefactor_phi
, &v
);
2684 /* Computation of forces per component */
2685 for (d
= 0; d
< DIM
; d
++)
2687 f_i
[d
] = prefactor_phi
2688 * (factor_phi_ai_ante
* delta_ante
[d
] + factor_phi_ai_crnt
* delta_crnt
[d
]
2689 + factor_phi_ai_post
* delta_post
[d
]);
2690 f_j
[d
] = prefactor_phi
2691 * (factor_phi_aj_ante
* delta_ante
[d
] + factor_phi_aj_crnt
* delta_crnt
[d
]
2692 + factor_phi_aj_post
* delta_post
[d
]);
2693 f_k
[d
] = prefactor_phi
2694 * (factor_phi_ak_ante
* delta_ante
[d
] + factor_phi_ak_crnt
* delta_crnt
[d
]
2695 + factor_phi_ak_post
* delta_post
[d
]);
2696 f_l
[d
] = prefactor_phi
2697 * (factor_phi_al_ante
* delta_ante
[d
] + factor_phi_al_crnt
* delta_crnt
[d
]
2698 + factor_phi_al_post
* delta_post
[d
]);
2700 /* Computation of the energy */
2705 /* Updating the forces */
2707 rvec_inc(f
[ai
], f_i
);
2708 rvec_inc(f
[aj
], f_j
);
2709 rvec_inc(f
[ak
], f_k
);
2710 rvec_inc(f
[al
], f_l
);
2712 if (computeVirial(flavor
))
2716 copy_ivec(SHIFT_IVEC(g
, aj
), jt
);
2717 ivec_sub(SHIFT_IVEC(g
, ai
), jt
, dt_ij
);
2718 ivec_sub(SHIFT_IVEC(g
, ak
), jt
, dt_kj
);
2719 ivec_sub(SHIFT_IVEC(g
, al
), jt
, dt_lj
);
2720 t1
= IVEC2IS(dt_ij
);
2721 t2
= IVEC2IS(dt_kj
);
2722 t3
= IVEC2IS(dt_lj
);
2726 t3
= pbc_rvec_sub(pbc
, x
[al
], x
[aj
], dx_jl
);
2733 rvec_inc(fshift
[t1
], f_i
);
2734 rvec_inc(fshift
[CENTRAL
], f_j
);
2735 rvec_inc(fshift
[t2
], f_k
);
2736 rvec_inc(fshift
[t3
], f_l
);
2744 template<BondedKernelFlavor flavor
>
2745 real
cbtdihs(int nbonds
,
2746 const t_iatom forceatoms
[],
2747 const t_iparams forceparams
[],
2753 real gmx_unused lambda
,
2754 real gmx_unused
* dvdlambda
,
2755 const t_mdatoms gmx_unused
* md
,
2756 t_fcdata gmx_unused
* fcd
,
2757 int gmx_unused
* global_atom_index
)
2759 int type
, ai
, aj
, ak
, al
, i
, d
;
2763 rvec f_i
, f_j
, f_k
, f_l
;
2764 ivec jt
, dt_ij
, dt_kj
, dt_lj
;
2766 rvec delta_ante
, delta_crnt
, delta_post
;
2767 rvec f_phi_ai
, f_phi_aj
, f_phi_ak
, f_phi_al
;
2768 rvec f_theta_ante_ai
, f_theta_ante_aj
, f_theta_ante_ak
;
2769 rvec f_theta_post_aj
, f_theta_post_ak
, f_theta_post_al
;
2773 for (i
= 0; (i
< nbonds
);)
2775 type
= forceatoms
[i
++];
2776 ai
= forceatoms
[i
++];
2777 aj
= forceatoms
[i
++];
2778 ak
= forceatoms
[i
++];
2779 al
= forceatoms
[i
++];
2782 t1
= pbc_rvec_sub(pbc
, x
[ai
], x
[aj
], vec_temp
);
2783 pbc_rvec_sub(pbc
, x
[aj
], x
[ai
], delta_ante
);
2784 t2
= pbc_rvec_sub(pbc
, x
[ak
], x
[aj
], vec_temp
);
2785 pbc_rvec_sub(pbc
, x
[ak
], x
[aj
], delta_crnt
);
2786 pbc_rvec_sub(pbc
, x
[ak
], x
[al
], vec_temp
);
2787 pbc_rvec_sub(pbc
, x
[al
], x
[ak
], delta_post
);
2789 /* \brief Compute factors for CBT potential
2790 * The combined bending-torsion potential goes to zero in a very smooth manner, eliminating the numerical
2791 * instabilities, when three coarse-grained particles align and the dihedral angle and standard
2792 * dihedral potentials cannot be calculated. The CBT potential is calculated using the formula:
2793 * \f[V_{\rm CBT}(\theta_{i-1}, \theta_i, \phi_i) = k_{\phi} \sin^3\theta_{i-1} \sin^3\theta_{i}
2794 * \sum_{n=0}^4 { a_n \cos^n\phi_i}\f] ({eq:CBT} and ref \cite{MonicaGoga2013} from the manual).
2795 * It contains in its expression not only the dihedral angle \f$\phi\f$
2796 * but also \f[\theta_{i-1}\f] (theta_ante bellow) and \f[\theta_{i}\f] (theta_post bellow)
2797 * --- the adjacent bending angles.
2798 * For more explanations see comments file "restcbt.h". */
2800 compute_factors_cbtdihs(type
, forceparams
, delta_ante
, delta_crnt
, delta_post
, f_phi_ai
,
2801 f_phi_aj
, f_phi_ak
, f_phi_al
, f_theta_ante_ai
, f_theta_ante_aj
,
2802 f_theta_ante_ak
, f_theta_post_aj
, f_theta_post_ak
, f_theta_post_al
, &v
);
2805 /* Acumulate the resuts per beads */
2806 for (d
= 0; d
< DIM
; d
++)
2808 f_i
[d
] = f_phi_ai
[d
] + f_theta_ante_ai
[d
];
2809 f_j
[d
] = f_phi_aj
[d
] + f_theta_ante_aj
[d
] + f_theta_post_aj
[d
];
2810 f_k
[d
] = f_phi_ak
[d
] + f_theta_ante_ak
[d
] + f_theta_post_ak
[d
];
2811 f_l
[d
] = f_phi_al
[d
] + f_theta_post_al
[d
];
2814 /* Compute the potential energy */
2819 /* Updating the forces */
2820 rvec_inc(f
[ai
], f_i
);
2821 rvec_inc(f
[aj
], f_j
);
2822 rvec_inc(f
[ak
], f_k
);
2823 rvec_inc(f
[al
], f_l
);
2826 if (computeVirial(flavor
))
2830 copy_ivec(SHIFT_IVEC(g
, aj
), jt
);
2831 ivec_sub(SHIFT_IVEC(g
, ai
), jt
, dt_ij
);
2832 ivec_sub(SHIFT_IVEC(g
, ak
), jt
, dt_kj
);
2833 ivec_sub(SHIFT_IVEC(g
, al
), jt
, dt_lj
);
2834 t1
= IVEC2IS(dt_ij
);
2835 t2
= IVEC2IS(dt_kj
);
2836 t3
= IVEC2IS(dt_lj
);
2840 t3
= pbc_rvec_sub(pbc
, x
[al
], x
[aj
], dx_jl
);
2847 rvec_inc(fshift
[t1
], f_i
);
2848 rvec_inc(fshift
[CENTRAL
], f_j
);
2849 rvec_inc(fshift
[t2
], f_k
);
2850 rvec_inc(fshift
[t3
], f_l
);
2857 template<BondedKernelFlavor flavor
>
2858 std::enable_if_t
<flavor
!= BondedKernelFlavor::ForcesSimdWhenAvailable
|| !GMX_SIMD_HAVE_REAL
, real
>
2860 const t_iatom forceatoms
[],
2861 const t_iparams forceparams
[],
2869 const t_mdatoms gmx_unused
* md
,
2870 t_fcdata gmx_unused
* fcd
,
2871 int gmx_unused
* global_atom_index
)
2873 const real c0
= 0.0, c1
= 1.0, c2
= 2.0, c3
= 3.0, c4
= 4.0, c5
= 5.0;
2874 int type
, ai
, aj
, ak
, al
, i
, j
;
2876 rvec r_ij
, r_kj
, r_kl
, m
, n
;
2877 real parmA
[NR_RBDIHS
];
2878 real parmB
[NR_RBDIHS
];
2879 real parm
[NR_RBDIHS
];
2880 real cos_phi
, phi
, rbp
, rbpBA
;
2881 real v
, ddphi
, sin_phi
;
2883 real L1
= 1.0 - lambda
;
2887 for (i
= 0; (i
< nbonds
);)
2889 type
= forceatoms
[i
++];
2890 ai
= forceatoms
[i
++];
2891 aj
= forceatoms
[i
++];
2892 ak
= forceatoms
[i
++];
2893 al
= forceatoms
[i
++];
2895 phi
= dih_angle(x
[ai
], x
[aj
], x
[ak
], x
[al
], pbc
, r_ij
, r_kj
, r_kl
, m
, n
, &t1
, &t2
, &t3
); /* 84 */
2897 /* Change to polymer convention */
2904 phi
-= M_PI
; /* 1 */
2906 cos_phi
= std::cos(phi
);
2907 /* Beware of accuracy loss, cannot use 1-sqrt(cos^2) ! */
2908 sin_phi
= std::sin(phi
);
2910 for (j
= 0; (j
< NR_RBDIHS
); j
++)
2912 parmA
[j
] = forceparams
[type
].rbdihs
.rbcA
[j
];
2913 parmB
[j
] = forceparams
[type
].rbdihs
.rbcB
[j
];
2914 parm
[j
] = L1
* parmA
[j
] + lambda
* parmB
[j
];
2916 /* Calculate cosine powers */
2917 /* Calculate the energy */
2918 /* Calculate the derivative */
2921 dvdl_term
+= (parmB
[0] - parmA
[0]);
2926 rbpBA
= parmB
[1] - parmA
[1];
2927 ddphi
+= rbp
* cosfac
;
2930 dvdl_term
+= cosfac
* rbpBA
;
2932 rbpBA
= parmB
[2] - parmA
[2];
2933 ddphi
+= c2
* rbp
* cosfac
;
2936 dvdl_term
+= cosfac
* rbpBA
;
2938 rbpBA
= parmB
[3] - parmA
[3];
2939 ddphi
+= c3
* rbp
* cosfac
;
2942 dvdl_term
+= cosfac
* rbpBA
;
2944 rbpBA
= parmB
[4] - parmA
[4];
2945 ddphi
+= c4
* rbp
* cosfac
;
2948 dvdl_term
+= cosfac
* rbpBA
;
2950 rbpBA
= parmB
[5] - parmA
[5];
2951 ddphi
+= c5
* rbp
* cosfac
;
2954 dvdl_term
+= cosfac
* rbpBA
;
2956 ddphi
= -ddphi
* sin_phi
; /* 11 */
2958 do_dih_fup
<flavor
>(ai
, aj
, ak
, al
, ddphi
, r_ij
, r_kj
, r_kl
, m
, n
, f
, fshift
, pbc
, g
, x
, t1
,
2962 *dvdlambda
+= dvdl_term
;
2969 /*! \brief Mysterious undocumented function */
2970 int cmap_setup_grid_index(int ip
, int grid_spacing
, int* ipm1
, int* ipp1
, int* ipp2
)
2976 ip
= ip
+ grid_spacing
- 1;
2978 else if (ip
> grid_spacing
)
2980 ip
= ip
- grid_spacing
- 1;
2989 im1
= grid_spacing
- 1;
2991 else if (ip
== grid_spacing
- 2)
2995 else if (ip
== grid_spacing
- 1)
3010 real
cmap_dihs(int nbonds
,
3011 const t_iatom forceatoms
[],
3012 const t_iparams forceparams
[],
3013 const gmx_cmap_t
* cmap_grid
,
3017 const struct t_pbc
* pbc
,
3018 const struct t_graph
* g
,
3019 real gmx_unused lambda
,
3020 real gmx_unused
* dvdlambda
,
3021 const t_mdatoms gmx_unused
* md
,
3022 t_fcdata gmx_unused
* fcd
,
3023 int gmx_unused
* global_atom_index
)
3026 int ai
, aj
, ak
, al
, am
;
3027 int a1i
, a1j
, a1k
, a1l
, a2i
, a2j
, a2k
, a2l
;
3029 int t11
, t21
, t31
, t12
, t22
, t32
;
3030 int iphi1
, ip1m1
, ip1p1
, ip1p2
;
3031 int iphi2
, ip2m1
, ip2p1
, ip2p2
;
3033 int pos1
, pos2
, pos3
, pos4
;
3035 real ty
[4], ty1
[4], ty2
[4], ty12
[4], tx
[16];
3036 real phi1
, cos_phi1
, sin_phi1
, xphi1
;
3037 real phi2
, cos_phi2
, sin_phi2
, xphi2
;
3038 real dx
, tt
, tu
, e
, df1
, df2
, vtot
;
3039 real ra21
, rb21
, rg21
, rg1
, rgr1
, ra2r1
, rb2r1
, rabr1
;
3040 real ra22
, rb22
, rg22
, rg2
, rgr2
, ra2r2
, rb2r2
, rabr2
;
3041 real fg1
, hg1
, fga1
, hgb1
, gaa1
, gbb1
;
3042 real fg2
, hg2
, fga2
, hgb2
, gaa2
, gbb2
;
3045 rvec r1_ij
, r1_kj
, r1_kl
, m1
, n1
;
3046 rvec r2_ij
, r2_kj
, r2_kl
, m2
, n2
;
3047 rvec f1_i
, f1_j
, f1_k
, f1_l
;
3048 rvec f2_i
, f2_j
, f2_k
, f2_l
;
3049 rvec a1
, b1
, a2
, b2
;
3050 rvec f1
, g1
, h1
, f2
, g2
, h2
;
3051 rvec dtf1
, dtg1
, dth1
, dtf2
, dtg2
, dth2
;
3052 ivec jt1
, dt1_ij
, dt1_kj
, dt1_lj
;
3053 ivec jt2
, dt2_ij
, dt2_kj
, dt2_lj
;
3055 int loop_index
[4][4] = { { 0, 4, 8, 12 }, { 1, 5, 9, 13 }, { 2, 6, 10, 14 }, { 3, 7, 11, 15 } };
3057 /* Total CMAP energy */
3060 for (n
= 0; n
< nbonds
;)
3062 /* Five atoms are involved in the two torsions */
3063 type
= forceatoms
[n
++];
3064 ai
= forceatoms
[n
++];
3065 aj
= forceatoms
[n
++];
3066 ak
= forceatoms
[n
++];
3067 al
= forceatoms
[n
++];
3068 am
= forceatoms
[n
++];
3070 /* Which CMAP type is this */
3071 const int cmapA
= forceparams
[type
].cmap
.cmapA
;
3072 const real
* cmapd
= cmap_grid
->cmapdata
[cmapA
].cmap
.data();
3080 phi1
= dih_angle(x
[a1i
], x
[a1j
], x
[a1k
], x
[a1l
], pbc
, r1_ij
, r1_kj
, r1_kl
, m1
, n1
, &t11
,
3081 &t21
, &t31
); /* 84 */
3083 cos_phi1
= std::cos(phi1
);
3085 a1
[0] = r1_ij
[1] * r1_kj
[2] - r1_ij
[2] * r1_kj
[1];
3086 a1
[1] = r1_ij
[2] * r1_kj
[0] - r1_ij
[0] * r1_kj
[2];
3087 a1
[2] = r1_ij
[0] * r1_kj
[1] - r1_ij
[1] * r1_kj
[0]; /* 9 */
3089 b1
[0] = r1_kl
[1] * r1_kj
[2] - r1_kl
[2] * r1_kj
[1];
3090 b1
[1] = r1_kl
[2] * r1_kj
[0] - r1_kl
[0] * r1_kj
[2];
3091 b1
[2] = r1_kl
[0] * r1_kj
[1] - r1_kl
[1] * r1_kj
[0]; /* 9 */
3093 pbc_rvec_sub(pbc
, x
[a1l
], x
[a1k
], h1
);
3095 ra21
= iprod(a1
, a1
); /* 5 */
3096 rb21
= iprod(b1
, b1
); /* 5 */
3097 rg21
= iprod(r1_kj
, r1_kj
); /* 5 */
3103 rabr1
= sqrt(ra2r1
* rb2r1
);
3105 sin_phi1
= rg1
* rabr1
* iprod(a1
, h1
) * (-1);
3107 if (cos_phi1
< -0.5 || cos_phi1
> 0.5)
3109 phi1
= std::asin(sin_phi1
);
3119 phi1
= -M_PI
- phi1
;
3125 phi1
= std::acos(cos_phi1
);
3133 xphi1
= phi1
+ M_PI
; /* 1 */
3135 /* Second torsion */
3141 phi2
= dih_angle(x
[a2i
], x
[a2j
], x
[a2k
], x
[a2l
], pbc
, r2_ij
, r2_kj
, r2_kl
, m2
, n2
, &t12
,
3142 &t22
, &t32
); /* 84 */
3144 cos_phi2
= std::cos(phi2
);
3146 a2
[0] = r2_ij
[1] * r2_kj
[2] - r2_ij
[2] * r2_kj
[1];
3147 a2
[1] = r2_ij
[2] * r2_kj
[0] - r2_ij
[0] * r2_kj
[2];
3148 a2
[2] = r2_ij
[0] * r2_kj
[1] - r2_ij
[1] * r2_kj
[0]; /* 9 */
3150 b2
[0] = r2_kl
[1] * r2_kj
[2] - r2_kl
[2] * r2_kj
[1];
3151 b2
[1] = r2_kl
[2] * r2_kj
[0] - r2_kl
[0] * r2_kj
[2];
3152 b2
[2] = r2_kl
[0] * r2_kj
[1] - r2_kl
[1] * r2_kj
[0]; /* 9 */
3154 pbc_rvec_sub(pbc
, x
[a2l
], x
[a2k
], h2
);
3156 ra22
= iprod(a2
, a2
); /* 5 */
3157 rb22
= iprod(b2
, b2
); /* 5 */
3158 rg22
= iprod(r2_kj
, r2_kj
); /* 5 */
3164 rabr2
= sqrt(ra2r2
* rb2r2
);
3166 sin_phi2
= rg2
* rabr2
* iprod(a2
, h2
) * (-1);
3168 if (cos_phi2
< -0.5 || cos_phi2
> 0.5)
3170 phi2
= std::asin(sin_phi2
);
3180 phi2
= -M_PI
- phi2
;
3186 phi2
= std::acos(cos_phi2
);
3194 xphi2
= phi2
+ M_PI
; /* 1 */
3196 /* Range mangling */
3199 xphi1
= xphi1
+ 2 * M_PI
;
3201 else if (xphi1
>= 2 * M_PI
)
3203 xphi1
= xphi1
- 2 * M_PI
;
3208 xphi2
= xphi2
+ 2 * M_PI
;
3210 else if (xphi2
>= 2 * M_PI
)
3212 xphi2
= xphi2
- 2 * M_PI
;
3215 /* Number of grid points */
3216 dx
= 2 * M_PI
/ cmap_grid
->grid_spacing
;
3218 /* Where on the grid are we */
3219 iphi1
= static_cast<int>(xphi1
/ dx
);
3220 iphi2
= static_cast<int>(xphi2
/ dx
);
3222 iphi1
= cmap_setup_grid_index(iphi1
, cmap_grid
->grid_spacing
, &ip1m1
, &ip1p1
, &ip1p2
);
3223 iphi2
= cmap_setup_grid_index(iphi2
, cmap_grid
->grid_spacing
, &ip2m1
, &ip2p1
, &ip2p2
);
3225 pos1
= iphi1
* cmap_grid
->grid_spacing
+ iphi2
;
3226 pos2
= ip1p1
* cmap_grid
->grid_spacing
+ iphi2
;
3227 pos3
= ip1p1
* cmap_grid
->grid_spacing
+ ip2p1
;
3228 pos4
= iphi1
* cmap_grid
->grid_spacing
+ ip2p1
;
3230 ty
[0] = cmapd
[pos1
* 4];
3231 ty
[1] = cmapd
[pos2
* 4];
3232 ty
[2] = cmapd
[pos3
* 4];
3233 ty
[3] = cmapd
[pos4
* 4];
3235 ty1
[0] = cmapd
[pos1
* 4 + 1];
3236 ty1
[1] = cmapd
[pos2
* 4 + 1];
3237 ty1
[2] = cmapd
[pos3
* 4 + 1];
3238 ty1
[3] = cmapd
[pos4
* 4 + 1];
3240 ty2
[0] = cmapd
[pos1
* 4 + 2];
3241 ty2
[1] = cmapd
[pos2
* 4 + 2];
3242 ty2
[2] = cmapd
[pos3
* 4 + 2];
3243 ty2
[3] = cmapd
[pos4
* 4 + 2];
3245 ty12
[0] = cmapd
[pos1
* 4 + 3];
3246 ty12
[1] = cmapd
[pos2
* 4 + 3];
3247 ty12
[2] = cmapd
[pos3
* 4 + 3];
3248 ty12
[3] = cmapd
[pos4
* 4 + 3];
3250 /* Switch to degrees */
3251 dx
= 360.0 / cmap_grid
->grid_spacing
;
3252 xphi1
= xphi1
* RAD2DEG
;
3253 xphi2
= xphi2
* RAD2DEG
;
3255 for (i
= 0; i
< 4; i
++) /* 16 */
3258 tx
[i
+ 4] = ty1
[i
] * dx
;
3259 tx
[i
+ 8] = ty2
[i
] * dx
;
3260 tx
[i
+ 12] = ty12
[i
] * dx
* dx
;
3263 real tc
[16] = { 0 };
3264 for (int idx
= 0; idx
< 16; idx
++) /* 1056 */
3266 for (int k
= 0; k
< 16; k
++)
3268 tc
[idx
] += cmap_coeff_matrix
[k
* 16 + idx
] * tx
[k
];
3272 tt
= (xphi1
- iphi1
* dx
) / dx
;
3273 tu
= (xphi2
- iphi2
* dx
) / dx
;
3279 for (i
= 3; i
>= 0; i
--)
3281 l1
= loop_index
[i
][3];
3282 l2
= loop_index
[i
][2];
3283 l3
= loop_index
[i
][1];
3285 e
= tt
* e
+ ((tc
[i
* 4 + 3] * tu
+ tc
[i
* 4 + 2]) * tu
+ tc
[i
* 4 + 1]) * tu
+ tc
[i
* 4];
3286 df1
= tu
* df1
+ (3.0 * tc
[l1
] * tt
+ 2.0 * tc
[l2
]) * tt
+ tc
[l3
];
3287 df2
= tt
* df2
+ (3.0 * tc
[i
* 4 + 3] * tu
+ 2.0 * tc
[i
* 4 + 2]) * tu
+ tc
[i
* 4 + 1];
3297 /* Do forces - first torsion */
3298 fg1
= iprod(r1_ij
, r1_kj
);
3299 hg1
= iprod(r1_kl
, r1_kj
);
3300 fga1
= fg1
* ra2r1
* rgr1
;
3301 hgb1
= hg1
* rb2r1
* rgr1
;
3302 gaa1
= -ra2r1
* rg1
;
3305 for (i
= 0; i
< DIM
; i
++)
3307 dtf1
[i
] = gaa1
* a1
[i
];
3308 dtg1
[i
] = fga1
* a1
[i
] - hgb1
* b1
[i
];
3309 dth1
[i
] = gbb1
* b1
[i
];
3311 f1
[i
] = df1
* dtf1
[i
];
3312 g1
[i
] = df1
* dtg1
[i
];
3313 h1
[i
] = df1
* dth1
[i
];
3316 f1_j
[i
] = -f1
[i
] - g1
[i
];
3317 f1_k
[i
] = h1
[i
] + g1
[i
];
3320 f
[a1i
][i
] = f
[a1i
][i
] + f1_i
[i
];
3321 f
[a1j
][i
] = f
[a1j
][i
] + f1_j
[i
]; /* - f1[i] - g1[i] */
3322 f
[a1k
][i
] = f
[a1k
][i
] + f1_k
[i
]; /* h1[i] + g1[i] */
3323 f
[a1l
][i
] = f
[a1l
][i
] + f1_l
[i
]; /* h1[i] */
3326 /* Do forces - second torsion */
3327 fg2
= iprod(r2_ij
, r2_kj
);
3328 hg2
= iprod(r2_kl
, r2_kj
);
3329 fga2
= fg2
* ra2r2
* rgr2
;
3330 hgb2
= hg2
* rb2r2
* rgr2
;
3331 gaa2
= -ra2r2
* rg2
;
3334 for (i
= 0; i
< DIM
; i
++)
3336 dtf2
[i
] = gaa2
* a2
[i
];
3337 dtg2
[i
] = fga2
* a2
[i
] - hgb2
* b2
[i
];
3338 dth2
[i
] = gbb2
* b2
[i
];
3340 f2
[i
] = df2
* dtf2
[i
];
3341 g2
[i
] = df2
* dtg2
[i
];
3342 h2
[i
] = df2
* dth2
[i
];
3345 f2_j
[i
] = -f2
[i
] - g2
[i
];
3346 f2_k
[i
] = h2
[i
] + g2
[i
];
3349 f
[a2i
][i
] = f
[a2i
][i
] + f2_i
[i
]; /* f2[i] */
3350 f
[a2j
][i
] = f
[a2j
][i
] + f2_j
[i
]; /* - f2[i] - g2[i] */
3351 f
[a2k
][i
] = f
[a2k
][i
] + f2_k
[i
]; /* h2[i] + g2[i] */
3352 f
[a2l
][i
] = f
[a2l
][i
] + f2_l
[i
]; /* - h2[i] */
3356 if (fshift
!= nullptr)
3360 copy_ivec(SHIFT_IVEC(g
, a1j
), jt1
);
3361 ivec_sub(SHIFT_IVEC(g
, a1i
), jt1
, dt1_ij
);
3362 ivec_sub(SHIFT_IVEC(g
, a1k
), jt1
, dt1_kj
);
3363 ivec_sub(SHIFT_IVEC(g
, a1l
), jt1
, dt1_lj
);
3364 t11
= IVEC2IS(dt1_ij
);
3365 t21
= IVEC2IS(dt1_kj
);
3366 t31
= IVEC2IS(dt1_lj
);
3368 copy_ivec(SHIFT_IVEC(g
, a2j
), jt2
);
3369 ivec_sub(SHIFT_IVEC(g
, a2i
), jt2
, dt2_ij
);
3370 ivec_sub(SHIFT_IVEC(g
, a2k
), jt2
, dt2_kj
);
3371 ivec_sub(SHIFT_IVEC(g
, a2l
), jt2
, dt2_lj
);
3372 t12
= IVEC2IS(dt2_ij
);
3373 t22
= IVEC2IS(dt2_kj
);
3374 t32
= IVEC2IS(dt2_lj
);
3378 t31
= pbc_rvec_sub(pbc
, x
[a1l
], x
[a1j
], h1
);
3379 t32
= pbc_rvec_sub(pbc
, x
[a2l
], x
[a2j
], h2
);
3387 rvec_inc(fshift
[t11
], f1_i
);
3388 rvec_inc(fshift
[CENTRAL
], f1_j
);
3389 rvec_inc(fshift
[t21
], f1_k
);
3390 rvec_inc(fshift
[t31
], f1_l
);
3392 rvec_inc(fshift
[t12
], f2_i
);
3393 rvec_inc(fshift
[CENTRAL
], f2_j
);
3394 rvec_inc(fshift
[t22
], f2_k
);
3395 rvec_inc(fshift
[t32
], f2_l
);
3405 /***********************************************************
3407 * G R O M O S 9 6 F U N C T I O N S
3409 ***********************************************************/
3410 real
g96harmonic(real kA
, real kB
, real xA
, real xB
, real x
, real lambda
, real
* V
, real
* F
)
3412 const real half
= 0.5;
3413 real L1
, kk
, x0
, dx
, dx2
;
3414 real v
, f
, dvdlambda
;
3417 kk
= L1
* kA
+ lambda
* kB
;
3418 x0
= L1
* xA
+ lambda
* xB
;
3424 v
= half
* kk
* dx2
;
3425 dvdlambda
= half
* (kB
- kA
) * dx2
+ (xA
- xB
) * kk
* dx
;
3432 /* That was 21 flops */
3435 template<BondedKernelFlavor flavor
>
3436 real
g96bonds(int nbonds
,
3437 const t_iatom forceatoms
[],
3438 const t_iparams forceparams
[],
3446 const t_mdatoms gmx_unused
* md
,
3447 t_fcdata gmx_unused
* fcd
,
3448 int gmx_unused
* global_atom_index
)
3450 int i
, ki
, ai
, aj
, type
;
3451 real dr2
, fbond
, vbond
, vtot
;
3455 for (i
= 0; (i
< nbonds
);)
3457 type
= forceatoms
[i
++];
3458 ai
= forceatoms
[i
++];
3459 aj
= forceatoms
[i
++];
3461 ki
= pbc_rvec_sub(pbc
, x
[ai
], x
[aj
], dx
); /* 3 */
3462 dr2
= iprod(dx
, dx
); /* 5 */
3464 *dvdlambda
+= g96harmonic(forceparams
[type
].harmonic
.krA
, forceparams
[type
].harmonic
.krB
,
3465 forceparams
[type
].harmonic
.rA
, forceparams
[type
].harmonic
.rB
, dr2
,
3466 lambda
, &vbond
, &fbond
);
3468 vtot
+= 0.5 * vbond
; /* 1*/
3470 spreadBondForces
<flavor
>(fbond
, dx
, ai
, aj
, f
, ki
, g
, fshift
); /* 15 */
3475 real
g96bond_angle(const rvec xi
, const rvec xj
, const rvec xk
, const t_pbc
* pbc
, rvec r_ij
, rvec r_kj
, int* t1
, int* t2
)
3476 /* Return value is the angle between the bonds i-j and j-k */
3480 *t1
= pbc_rvec_sub(pbc
, xi
, xj
, r_ij
); /* 3 */
3481 *t2
= pbc_rvec_sub(pbc
, xk
, xj
, r_kj
); /* 3 */
3483 costh
= cos_angle(r_ij
, r_kj
); /* 25 */
3488 template<BondedKernelFlavor flavor
>
3489 real
g96angles(int nbonds
,
3490 const t_iatom forceatoms
[],
3491 const t_iparams forceparams
[],
3499 const t_mdatoms gmx_unused
* md
,
3500 t_fcdata gmx_unused
* fcd
,
3501 int gmx_unused
* global_atom_index
)
3503 int i
, ai
, aj
, ak
, type
, m
, t1
, t2
;
3505 real cos_theta
, dVdt
, va
, vtot
;
3506 real rij_1
, rij_2
, rkj_1
, rkj_2
, rijrkj_1
;
3508 ivec jt
, dt_ij
, dt_kj
;
3511 for (i
= 0; (i
< nbonds
);)
3513 type
= forceatoms
[i
++];
3514 ai
= forceatoms
[i
++];
3515 aj
= forceatoms
[i
++];
3516 ak
= forceatoms
[i
++];
3518 cos_theta
= g96bond_angle(x
[ai
], x
[aj
], x
[ak
], pbc
, r_ij
, r_kj
, &t1
, &t2
);
3520 *dvdlambda
+= g96harmonic(forceparams
[type
].harmonic
.krA
, forceparams
[type
].harmonic
.krB
,
3521 forceparams
[type
].harmonic
.rA
, forceparams
[type
].harmonic
.rB
,
3522 cos_theta
, lambda
, &va
, &dVdt
);
3525 rij_1
= gmx::invsqrt(iprod(r_ij
, r_ij
));
3526 rkj_1
= gmx::invsqrt(iprod(r_kj
, r_kj
));
3527 rij_2
= rij_1
* rij_1
;
3528 rkj_2
= rkj_1
* rkj_1
;
3529 rijrkj_1
= rij_1
* rkj_1
; /* 23 */
3531 for (m
= 0; (m
< DIM
); m
++) /* 42 */
3533 f_i
[m
] = dVdt
* (r_kj
[m
] * rijrkj_1
- r_ij
[m
] * rij_2
* cos_theta
);
3534 f_k
[m
] = dVdt
* (r_ij
[m
] * rijrkj_1
- r_kj
[m
] * rkj_2
* cos_theta
);
3535 f_j
[m
] = -f_i
[m
] - f_k
[m
];
3541 if (computeVirial(flavor
))
3545 copy_ivec(SHIFT_IVEC(g
, aj
), jt
);
3547 ivec_sub(SHIFT_IVEC(g
, ai
), jt
, dt_ij
);
3548 ivec_sub(SHIFT_IVEC(g
, ak
), jt
, dt_kj
);
3549 t1
= IVEC2IS(dt_ij
);
3550 t2
= IVEC2IS(dt_kj
);
3552 rvec_inc(fshift
[t1
], f_i
);
3553 rvec_inc(fshift
[CENTRAL
], f_j
);
3554 rvec_inc(fshift
[t2
], f_k
); /* 9 */
3561 template<BondedKernelFlavor flavor
>
3562 real
cross_bond_bond(int nbonds
,
3563 const t_iatom forceatoms
[],
3564 const t_iparams forceparams
[],
3570 real gmx_unused lambda
,
3571 real gmx_unused
* dvdlambda
,
3572 const t_mdatoms gmx_unused
* md
,
3573 t_fcdata gmx_unused
* fcd
,
3574 int gmx_unused
* global_atom_index
)
3576 /* Potential from Lawrence and Skimmer, Chem. Phys. Lett. 372 (2003)
3579 int i
, ai
, aj
, ak
, type
, m
, t1
, t2
;
3581 real vtot
, vrr
, s1
, s2
, r1
, r2
, r1e
, r2e
, krr
;
3583 ivec jt
, dt_ij
, dt_kj
;
3586 for (i
= 0; (i
< nbonds
);)
3588 type
= forceatoms
[i
++];
3589 ai
= forceatoms
[i
++];
3590 aj
= forceatoms
[i
++];
3591 ak
= forceatoms
[i
++];
3592 r1e
= forceparams
[type
].cross_bb
.r1e
;
3593 r2e
= forceparams
[type
].cross_bb
.r2e
;
3594 krr
= forceparams
[type
].cross_bb
.krr
;
3596 /* Compute distance vectors ... */
3597 t1
= pbc_rvec_sub(pbc
, x
[ai
], x
[aj
], r_ij
);
3598 t2
= pbc_rvec_sub(pbc
, x
[ak
], x
[aj
], r_kj
);
3600 /* ... and their lengths */
3604 /* Deviations from ideality */
3608 /* Energy (can be negative!) */
3609 vrr
= krr
* s1
* s2
;
3613 svmul(-krr
* s2
/ r1
, r_ij
, f_i
);
3614 svmul(-krr
* s1
/ r2
, r_kj
, f_k
);
3616 for (m
= 0; (m
< DIM
); m
++) /* 12 */
3618 f_j
[m
] = -f_i
[m
] - f_k
[m
];
3624 if (computeVirial(flavor
))
3628 copy_ivec(SHIFT_IVEC(g
, aj
), jt
);
3630 ivec_sub(SHIFT_IVEC(g
, ai
), jt
, dt_ij
);
3631 ivec_sub(SHIFT_IVEC(g
, ak
), jt
, dt_kj
);
3632 t1
= IVEC2IS(dt_ij
);
3633 t2
= IVEC2IS(dt_kj
);
3635 rvec_inc(fshift
[t1
], f_i
);
3636 rvec_inc(fshift
[CENTRAL
], f_j
);
3637 rvec_inc(fshift
[t2
], f_k
); /* 9 */
3644 template<BondedKernelFlavor flavor
>
3645 real
cross_bond_angle(int nbonds
,
3646 const t_iatom forceatoms
[],
3647 const t_iparams forceparams
[],
3653 real gmx_unused lambda
,
3654 real gmx_unused
* dvdlambda
,
3655 const t_mdatoms gmx_unused
* md
,
3656 t_fcdata gmx_unused
* fcd
,
3657 int gmx_unused
* global_atom_index
)
3659 /* Potential from Lawrence and Skimmer, Chem. Phys. Lett. 372 (2003)
3662 int i
, ai
, aj
, ak
, type
, m
, t1
, t2
;
3663 rvec r_ij
, r_kj
, r_ik
;
3664 real vtot
, vrt
, s1
, s2
, s3
, r1
, r2
, r3
, r1e
, r2e
, r3e
, krt
, k1
, k2
, k3
;
3666 ivec jt
, dt_ij
, dt_kj
;
3669 for (i
= 0; (i
< nbonds
);)
3671 type
= forceatoms
[i
++];
3672 ai
= forceatoms
[i
++];
3673 aj
= forceatoms
[i
++];
3674 ak
= forceatoms
[i
++];
3675 r1e
= forceparams
[type
].cross_ba
.r1e
;
3676 r2e
= forceparams
[type
].cross_ba
.r2e
;
3677 r3e
= forceparams
[type
].cross_ba
.r3e
;
3678 krt
= forceparams
[type
].cross_ba
.krt
;
3680 /* Compute distance vectors ... */
3681 t1
= pbc_rvec_sub(pbc
, x
[ai
], x
[aj
], r_ij
);
3682 t2
= pbc_rvec_sub(pbc
, x
[ak
], x
[aj
], r_kj
);
3683 pbc_rvec_sub(pbc
, x
[ai
], x
[ak
], r_ik
);
3685 /* ... and their lengths */
3690 /* Deviations from ideality */
3695 /* Energy (can be negative!) */
3696 vrt
= krt
* s3
* (s1
+ s2
);
3700 k1
= -krt
* (s3
/ r1
);
3701 k2
= -krt
* (s3
/ r2
);
3702 k3
= -krt
* (s1
+ s2
) / r3
;
3703 for (m
= 0; (m
< DIM
); m
++)
3705 f_i
[m
] = k1
* r_ij
[m
] + k3
* r_ik
[m
];
3706 f_k
[m
] = k2
* r_kj
[m
] - k3
* r_ik
[m
];
3707 f_j
[m
] = -f_i
[m
] - f_k
[m
];
3710 for (m
= 0; (m
< DIM
); m
++) /* 12 */
3717 if (computeVirial(flavor
))
3721 copy_ivec(SHIFT_IVEC(g
, aj
), jt
);
3723 ivec_sub(SHIFT_IVEC(g
, ai
), jt
, dt_ij
);
3724 ivec_sub(SHIFT_IVEC(g
, ak
), jt
, dt_kj
);
3725 t1
= IVEC2IS(dt_ij
);
3726 t2
= IVEC2IS(dt_kj
);
3728 rvec_inc(fshift
[t1
], f_i
);
3729 rvec_inc(fshift
[CENTRAL
], f_j
);
3730 rvec_inc(fshift
[t2
], f_k
); /* 9 */
3737 /*! \brief Computes the potential and force for a tabulated potential */
3738 real
bonded_tab(const char* type
,
3740 const bondedtable_t
* table
,
3748 real k
, tabscale
, *VFtab
, rt
, eps
, eps2
, Yt
, Ft
, Geps
, Heps2
, Fp
, VV
, FF
;
3752 k
= (1.0 - lambda
) * kA
+ lambda
* kB
;
3754 tabscale
= table
->scale
;
3755 VFtab
= table
->data
;
3758 n0
= static_cast<int>(rt
);
3762 "A tabulated %s interaction table number %d is out of the table range: r %f, "
3763 "between table indices %d and %d, table length %d",
3764 type
, table_nr
, r
, n0
, n0
+ 1, table
->n
);
3770 Ft
= VFtab
[nnn
+ 1];
3771 Geps
= VFtab
[nnn
+ 2] * eps
;
3772 Heps2
= VFtab
[nnn
+ 3] * eps2
;
3773 Fp
= Ft
+ Geps
+ Heps2
;
3775 FF
= Fp
+ Geps
+ 2.0 * Heps2
;
3777 *F
= -k
* FF
* tabscale
;
3779 dvdlambda
= (kB
- kA
) * VV
;
3783 /* That was 22 flops */
3786 template<BondedKernelFlavor flavor
>
3787 real
tab_bonds(int nbonds
,
3788 const t_iatom forceatoms
[],
3789 const t_iparams forceparams
[],
3797 const t_mdatoms gmx_unused
* md
,
3799 int gmx_unused
* global_atom_index
)
3801 int i
, ki
, ai
, aj
, type
, table
;
3802 real dr
, dr2
, fbond
, vbond
, vtot
;
3806 for (i
= 0; (i
< nbonds
);)
3808 type
= forceatoms
[i
++];
3809 ai
= forceatoms
[i
++];
3810 aj
= forceatoms
[i
++];
3812 ki
= pbc_rvec_sub(pbc
, x
[ai
], x
[aj
], dx
); /* 3 */
3813 dr2
= iprod(dx
, dx
); /* 5 */
3814 dr
= dr2
* gmx::invsqrt(dr2
); /* 10 */
3816 table
= forceparams
[type
].tab
.table
;
3818 *dvdlambda
+= bonded_tab("bond", table
, &fcd
->bondtab
[table
], forceparams
[type
].tab
.kA
,
3819 forceparams
[type
].tab
.kB
, dr
, lambda
, &vbond
, &fbond
); /* 22 */
3827 vtot
+= vbond
; /* 1*/
3828 fbond
*= gmx::invsqrt(dr2
); /* 6 */
3830 spreadBondForces
<flavor
>(fbond
, dx
, ai
, aj
, f
, ki
, g
, fshift
); /* 15 */
3835 template<BondedKernelFlavor flavor
>
3836 real
tab_angles(int nbonds
,
3837 const t_iatom forceatoms
[],
3838 const t_iparams forceparams
[],
3846 const t_mdatoms gmx_unused
* md
,
3848 int gmx_unused
* global_atom_index
)
3850 int i
, ai
, aj
, ak
, t1
, t2
, type
, table
;
3852 real cos_theta
, cos_theta2
, theta
, dVdt
, va
, vtot
;
3853 ivec jt
, dt_ij
, dt_kj
;
3856 for (i
= 0; (i
< nbonds
);)
3858 type
= forceatoms
[i
++];
3859 ai
= forceatoms
[i
++];
3860 aj
= forceatoms
[i
++];
3861 ak
= forceatoms
[i
++];
3863 theta
= bond_angle(x
[ai
], x
[aj
], x
[ak
], pbc
, r_ij
, r_kj
, &cos_theta
, &t1
, &t2
); /* 41 */
3865 table
= forceparams
[type
].tab
.table
;
3867 *dvdlambda
+= bonded_tab("angle", table
, &fcd
->angletab
[table
], forceparams
[type
].tab
.kA
,
3868 forceparams
[type
].tab
.kB
, theta
, lambda
, &va
, &dVdt
); /* 22 */
3871 cos_theta2
= gmx::square(cos_theta
); /* 1 */
3880 st
= dVdt
* gmx::invsqrt(1 - cos_theta2
); /* 12 */
3881 sth
= st
* cos_theta
; /* 1 */
3882 nrkj2
= iprod(r_kj
, r_kj
); /* 5 */
3883 nrij2
= iprod(r_ij
, r_ij
);
3885 cik
= st
* gmx::invsqrt(nrkj2
* nrij2
); /* 12 */
3886 cii
= sth
/ nrij2
; /* 10 */
3887 ckk
= sth
/ nrkj2
; /* 10 */
3889 for (m
= 0; (m
< DIM
); m
++) /* 39 */
3891 f_i
[m
] = -(cik
* r_kj
[m
] - cii
* r_ij
[m
]);
3892 f_k
[m
] = -(cik
* r_ij
[m
] - ckk
* r_kj
[m
]);
3893 f_j
[m
] = -f_i
[m
] - f_k
[m
];
3899 if (computeVirial(flavor
))
3903 copy_ivec(SHIFT_IVEC(g
, aj
), jt
);
3905 ivec_sub(SHIFT_IVEC(g
, ai
), jt
, dt_ij
);
3906 ivec_sub(SHIFT_IVEC(g
, ak
), jt
, dt_kj
);
3907 t1
= IVEC2IS(dt_ij
);
3908 t2
= IVEC2IS(dt_kj
);
3910 rvec_inc(fshift
[t1
], f_i
);
3911 rvec_inc(fshift
[CENTRAL
], f_j
);
3912 rvec_inc(fshift
[t2
], f_k
);
3919 template<BondedKernelFlavor flavor
>
3920 real
tab_dihs(int nbonds
,
3921 const t_iatom forceatoms
[],
3922 const t_iparams forceparams
[],
3930 const t_mdatoms gmx_unused
* md
,
3932 int gmx_unused
* global_atom_index
)
3934 int i
, type
, ai
, aj
, ak
, al
, table
;
3936 rvec r_ij
, r_kj
, r_kl
, m
, n
;
3937 real phi
, ddphi
, vpd
, vtot
;
3940 for (i
= 0; (i
< nbonds
);)
3942 type
= forceatoms
[i
++];
3943 ai
= forceatoms
[i
++];
3944 aj
= forceatoms
[i
++];
3945 ak
= forceatoms
[i
++];
3946 al
= forceatoms
[i
++];
3948 phi
= dih_angle(x
[ai
], x
[aj
], x
[ak
], x
[al
], pbc
, r_ij
, r_kj
, r_kl
, m
, n
, &t1
, &t2
, &t3
); /* 84 */
3950 table
= forceparams
[type
].tab
.table
;
3952 /* Hopefully phi+M_PI never results in values < 0 */
3953 *dvdlambda
+= bonded_tab("dihedral", table
, &fcd
->dihtab
[table
], forceparams
[type
].tab
.kA
,
3954 forceparams
[type
].tab
.kB
, phi
+ M_PI
, lambda
, &vpd
, &ddphi
);
3957 do_dih_fup
<flavor
>(ai
, aj
, ak
, al
, -ddphi
, r_ij
, r_kj
, r_kl
, m
, n
, f
, fshift
, pbc
, g
, x
, t1
,
3965 struct BondedInteractions
3967 BondedFunction function
;
3971 /*! \brief Lookup table of bonded interaction functions
3973 * This must have as many entries as interaction_function in ifunc.cpp */
3974 template<BondedKernelFlavor flavor
>
3975 const std::array
<BondedInteractions
, F_NRE
> c_bondedInteractionFunctions
= {
3976 BondedInteractions
{ bonds
<flavor
>, eNR_BONDS
}, // F_BONDS
3977 BondedInteractions
{ g96bonds
<flavor
>, eNR_BONDS
}, // F_G96BONDS
3978 BondedInteractions
{ morse_bonds
<flavor
>, eNR_MORSE
}, // F_MORSE
3979 BondedInteractions
{ cubic_bonds
<flavor
>, eNR_CUBICBONDS
}, // F_CUBICBONDS
3980 BondedInteractions
{ unimplemented
, -1 }, // F_CONNBONDS
3981 BondedInteractions
{ bonds
<flavor
>, eNR_BONDS
}, // F_HARMONIC
3982 BondedInteractions
{ FENE_bonds
<flavor
>, eNR_FENEBONDS
}, // F_FENEBONDS
3983 BondedInteractions
{ tab_bonds
<flavor
>, eNR_TABBONDS
}, // F_TABBONDS
3984 BondedInteractions
{ tab_bonds
<flavor
>, eNR_TABBONDS
}, // F_TABBONDSNC
3985 BondedInteractions
{ restraint_bonds
<flavor
>, eNR_RESTRBONDS
}, // F_RESTRBONDS
3986 BondedInteractions
{ angles
<flavor
>, eNR_ANGLES
}, // F_ANGLES
3987 BondedInteractions
{ g96angles
<flavor
>, eNR_ANGLES
}, // F_G96ANGLES
3988 BondedInteractions
{ restrangles
<flavor
>, eNR_ANGLES
}, // F_RESTRANGLES
3989 BondedInteractions
{ linear_angles
<flavor
>, eNR_ANGLES
}, // F_LINEAR_ANGLES
3990 BondedInteractions
{ cross_bond_bond
<flavor
>, eNR_CROSS_BOND_BOND
}, // F_CROSS_BOND_BONDS
3991 BondedInteractions
{ cross_bond_angle
<flavor
>, eNR_CROSS_BOND_ANGLE
}, // F_CROSS_BOND_ANGLES
3992 BondedInteractions
{ urey_bradley
<flavor
>, eNR_UREY_BRADLEY
}, // F_UREY_BRADLEY
3993 BondedInteractions
{ quartic_angles
<flavor
>, eNR_QANGLES
}, // F_QUARTIC_ANGLES
3994 BondedInteractions
{ tab_angles
<flavor
>, eNR_TABANGLES
}, // F_TABANGLES
3995 BondedInteractions
{ pdihs
<flavor
>, eNR_PROPER
}, // F_PDIHS
3996 BondedInteractions
{ rbdihs
<flavor
>, eNR_RB
}, // F_RBDIHS
3997 BondedInteractions
{ restrdihs
<flavor
>, eNR_PROPER
}, // F_RESTRDIHS
3998 BondedInteractions
{ cbtdihs
<flavor
>, eNR_RB
}, // F_CBTDIHS
3999 BondedInteractions
{ rbdihs
<flavor
>, eNR_FOURDIH
}, // F_FOURDIHS
4000 BondedInteractions
{ idihs
<flavor
>, eNR_IMPROPER
}, // F_IDIHS
4001 BondedInteractions
{ pdihs
<flavor
>, eNR_IMPROPER
}, // F_PIDIHS
4002 BondedInteractions
{ tab_dihs
<flavor
>, eNR_TABDIHS
}, // F_TABDIHS
4003 BondedInteractions
{ unimplemented
, eNR_CMAP
}, // F_CMAP
4004 BondedInteractions
{ unimplemented
, -1 }, // F_GB12_NOLONGERUSED
4005 BondedInteractions
{ unimplemented
, -1 }, // F_GB13_NOLONGERUSED
4006 BondedInteractions
{ unimplemented
, -1 }, // F_GB14_NOLONGERUSED
4007 BondedInteractions
{ unimplemented
, -1 }, // F_GBPOL_NOLONGERUSED
4008 BondedInteractions
{ unimplemented
, -1 }, // F_NPSOLVATION_NOLONGERUSED
4009 BondedInteractions
{ unimplemented
, eNR_NB14
}, // F_LJ14
4010 BondedInteractions
{ unimplemented
, -1 }, // F_COUL14
4011 BondedInteractions
{ unimplemented
, eNR_NB14
}, // F_LJC14_Q
4012 BondedInteractions
{ unimplemented
, eNR_NB14
}, // F_LJC_PAIRS_NB
4013 BondedInteractions
{ unimplemented
, -1 }, // F_LJ
4014 BondedInteractions
{ unimplemented
, -1 }, // F_BHAM
4015 BondedInteractions
{ unimplemented
, -1 }, // F_LJ_LR_NOLONGERUSED
4016 BondedInteractions
{ unimplemented
, -1 }, // F_BHAM_LR_NOLONGERUSED
4017 BondedInteractions
{ unimplemented
, -1 }, // F_DISPCORR
4018 BondedInteractions
{ unimplemented
, -1 }, // F_COUL_SR
4019 BondedInteractions
{ unimplemented
, -1 }, // F_COUL_LR_NOLONGERUSED
4020 BondedInteractions
{ unimplemented
, -1 }, // F_RF_EXCL
4021 BondedInteractions
{ unimplemented
, -1 }, // F_COUL_RECIP
4022 BondedInteractions
{ unimplemented
, -1 }, // F_LJ_RECIP
4023 BondedInteractions
{ unimplemented
, -1 }, // F_DPD
4024 BondedInteractions
{ polarize
<flavor
>, eNR_POLARIZE
}, // F_POLARIZATION
4025 BondedInteractions
{ water_pol
<flavor
>, eNR_WPOL
}, // F_WATER_POL
4026 BondedInteractions
{ thole_pol
<flavor
>, eNR_THOLE
}, // F_THOLE_POL
4027 BondedInteractions
{ anharm_polarize
<flavor
>, eNR_ANHARM_POL
}, // F_ANHARM_POL
4028 BondedInteractions
{ unimplemented
, -1 }, // F_POSRES
4029 BondedInteractions
{ unimplemented
, -1 }, // F_FBPOSRES
4030 BondedInteractions
{ ta_disres
, eNR_DISRES
}, // F_DISRES
4031 BondedInteractions
{ unimplemented
, -1 }, // F_DISRESVIOL
4032 BondedInteractions
{ orires
, eNR_ORIRES
}, // F_ORIRES
4033 BondedInteractions
{ unimplemented
, -1 }, // F_ORIRESDEV
4034 BondedInteractions
{ angres
<flavor
>, eNR_ANGRES
}, // F_ANGRES
4035 BondedInteractions
{ angresz
<flavor
>, eNR_ANGRESZ
}, // F_ANGRESZ
4036 BondedInteractions
{ dihres
<flavor
>, eNR_DIHRES
}, // F_DIHRES
4037 BondedInteractions
{ unimplemented
, -1 }, // F_DIHRESVIOL
4038 BondedInteractions
{ unimplemented
, -1 }, // F_CONSTR
4039 BondedInteractions
{ unimplemented
, -1 }, // F_CONSTRNC
4040 BondedInteractions
{ unimplemented
, -1 }, // F_SETTLE
4041 BondedInteractions
{ unimplemented
, -1 }, // F_VSITE2
4042 BondedInteractions
{ unimplemented
, -1 }, // F_VSITE3
4043 BondedInteractions
{ unimplemented
, -1 }, // F_VSITE3FD
4044 BondedInteractions
{ unimplemented
, -1 }, // F_VSITE3FAD
4045 BondedInteractions
{ unimplemented
, -1 }, // F_VSITE3OUT
4046 BondedInteractions
{ unimplemented
, -1 }, // F_VSITE4FD
4047 BondedInteractions
{ unimplemented
, -1 }, // F_VSITE4FDN
4048 BondedInteractions
{ unimplemented
, -1 }, // F_VSITEN
4049 BondedInteractions
{ unimplemented
, -1 }, // F_COM_PULL
4050 BondedInteractions
{ unimplemented
, -1 }, // F_DENSITYFITTING
4051 BondedInteractions
{ unimplemented
, -1 }, // F_EQM
4052 BondedInteractions
{ unimplemented
, -1 }, // F_EPOT
4053 BondedInteractions
{ unimplemented
, -1 }, // F_EKIN
4054 BondedInteractions
{ unimplemented
, -1 }, // F_ETOT
4055 BondedInteractions
{ unimplemented
, -1 }, // F_ECONSERVED
4056 BondedInteractions
{ unimplemented
, -1 }, // F_TEMP
4057 BondedInteractions
{ unimplemented
, -1 }, // F_VTEMP_NOLONGERUSED
4058 BondedInteractions
{ unimplemented
, -1 }, // F_PDISPCORR
4059 BondedInteractions
{ unimplemented
, -1 }, // F_PRES
4060 BondedInteractions
{ unimplemented
, -1 }, // F_DVDL_CONSTR
4061 BondedInteractions
{ unimplemented
, -1 }, // F_DVDL
4062 BondedInteractions
{ unimplemented
, -1 }, // F_DKDL
4063 BondedInteractions
{ unimplemented
, -1 }, // F_DVDL_COUL
4064 BondedInteractions
{ unimplemented
, -1 }, // F_DVDL_VDW
4065 BondedInteractions
{ unimplemented
, -1 }, // F_DVDL_BONDED
4066 BondedInteractions
{ unimplemented
, -1 }, // F_DVDL_RESTRAINT
4067 BondedInteractions
{ unimplemented
, -1 }, // F_DVDL_TEMPERATURE
4070 /*! \brief List of instantiated BondedInteractions list */
4071 const gmx::EnumerationArray
<BondedKernelFlavor
, std::array
<BondedInteractions
, F_NRE
>> c_bondedInteractionFunctionsPerFlavor
= {
4072 c_bondedInteractionFunctions
<BondedKernelFlavor::ForcesSimdWhenAvailable
>,
4073 c_bondedInteractionFunctions
<BondedKernelFlavor::ForcesNoSimd
>,
4074 c_bondedInteractionFunctions
<BondedKernelFlavor::ForcesAndVirialAndEnergy
>,
4075 c_bondedInteractionFunctions
<BondedKernelFlavor::ForcesAndEnergy
>
4082 real
calculateSimpleBond(const int ftype
,
4083 const int numForceatoms
,
4084 const t_iatom forceatoms
[],
4085 const t_iparams forceparams
[],
4089 const struct t_pbc
* pbc
,
4090 const struct t_graph gmx_unused
* g
,
4093 const t_mdatoms
* md
,
4095 int gmx_unused
* global_atom_index
,
4096 const BondedKernelFlavor bondedKernelFlavor
)
4098 const BondedInteractions
& bonded
= c_bondedInteractionFunctionsPerFlavor
[bondedKernelFlavor
][ftype
];
4100 real v
= bonded
.function(numForceatoms
, forceatoms
, forceparams
, x
, f
, fshift
, pbc
, g
, lambda
,
4101 dvdlambda
, md
, fcd
, global_atom_index
);
4106 int nrnbIndex(int ftype
)
4108 return c_bondedInteractionFunctions
<BondedKernelFlavor::ForcesAndVirialAndEnergy
>[ftype
].nrnbIndex
;