?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Precision fix for rescbt code.
[gromacs.git] / src / gromacs / listed_forces / bonded.cpp
blob7c4f176ee8fe713082e3e90670c01293c3011420
1 /*
2 * This file is part of the GROMACS molecular simulation package.
3 *
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.
10 *
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.
15 *
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.
20 *
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.
25 *
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.
33 *
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.
36 */
37 /*! \internal \file
38 *
39 * \brief This file defines low-level functions necessary for
40 * computing energies and forces for listed interactions.
41 *
42 * \author Mark Abraham <mark.j.abraham@gmail.com>
43 *
44 * \ingroup module_listed_forces
45 */
46 #include "gmxpre.h"
47
48 #include "bonded.h"
49
50 #include "config.h"
51
52 #include <cassert>
53 #include <cmath>
54
55 #include <algorithm>
56
57 #include "gromacs/listed_forces/pairs.h"
58 #include "gromacs/math/functions.h"
59 #include "gromacs/math/units.h"
60 #include "gromacs/math/utilities.h"
61 #include "gromacs/math/vec.h"
62 #include "gromacs/mdtypes/fcdata.h"
63 #include "gromacs/pbcutil/ishift.h"
64 #include "gromacs/pbcutil/mshift.h"
65 #include "gromacs/pbcutil/pbc.h"
66 #include "gromacs/pbcutil/pbc_simd.h"
67 #include "gromacs/simd/simd.h"
68 #include "gromacs/simd/simd_math.h"
69 #include "gromacs/simd/vector_operations.h"
70 #include "gromacs/utility/basedefinitions.h"
71 #include "gromacs/utility/fatalerror.h"
72 #include "gromacs/utility/real.h"
73 #include "gromacs/utility/smalloc.h"
74
75 #include "listed_internal.h"
76 #include "restcbt.h"
77
78 using namespace gmx; // TODO: Remove when this file is moved into gmx namespace
79
80 /*! \brief Mysterious CMAP coefficient matrix */
81 const int cmap_coeff_matrix[] = {
82 1, 0, -3, 2, 0, 0, 0, 0, -3, 0, 9, -6, 2, 0, -6, 4,
83 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, -9, 6, -2, 0, 6, -4,
84 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 9, -6, 0, 0, -6, 4,
85 0, 0, 3, -2, 0, 0, 0, 0, 0, 0, -9, 6, 0, 0, 6, -4,
86 0, 0, 0, 0, 1, 0, -3, 2, -2, 0, 6, -4, 1, 0, -3, 2,
87 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 3, -2, 1, 0, -3, 2,
88 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 2, 0, 0, 3, -2,
89 0, 0, 0, 0, 0, 0, 3, -2, 0, 0, -6, 4, 0, 0, 3, -2,
90 0, 1, -2, 1, 0, 0, 0, 0, 0, -3, 6, -3, 0, 2, -4, 2,
91 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, -6, 3, 0, -2, 4, -2,
92 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 3, 0, 0, 2, -2,
93 0, 0, -1, 1, 0, 0, 0, 0, 0, 0, 3, -3, 0, 0, -2, 2,
94 0, 0, 0, 0, 0, 1, -2, 1, 0, -2, 4, -2, 0, 1, -2, 1,
95 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 2, -1, 0, 1, -2, 1,
96 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, -1, 1,
97 0, 0, 0, 0, 0, 0, -1, 1, 0, 0, 2, -2, 0, 0, -1, 1
98 };
99
100
101 /*! \brief Compute dx = xi - xj, modulo PBC if non-NULL
102 *
103 * \todo This kind of code appears in many places. Consolidate it */
104 static int pbc_rvec_sub(const t_pbc *pbc, const rvec xi, const rvec xj, rvec dx)
105 {
106 if (pbc)
107 {
108 return pbc_dx_aiuc(pbc, xi, xj, dx);
109 }
110 else
111 {
112 rvec_sub(xi, xj, dx);
113 return CENTRAL;
114 }
115 }
116
117 /*! \brief Morse potential bond
118 *
119 * By Frank Everdij. Three parameters needed:
120 *
121 * b0 = equilibrium distance in nm
122 * be = beta in nm^-1 (actually, it's nu_e*Sqrt(2*pi*pi*mu/D_e))
123 * cb = well depth in kJ/mol
124 *
125 * Note: the potential is referenced to be +cb at infinite separation
126 * and zero at the equilibrium distance!
127 */
128 real morse_bonds(int nbonds,
129 const t_iatom forceatoms[], const t_iparams forceparams[],
130 const rvec x[], rvec4 f[], rvec fshift[],
131 const t_pbc *pbc, const t_graph *g,
132 real lambda, real *dvdlambda,
133 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
134 int gmx_unused *global_atom_index)
135 {
136 const real one = 1.0;
137 const real two = 2.0;
138 real dr, dr2, temp, omtemp, cbomtemp, fbond, vbond, fij, vtot;
139 real b0, be, cb, b0A, beA, cbA, b0B, beB, cbB, L1;
140 rvec dx;
141 int i, m, ki, type, ai, aj;
142 ivec dt;
143
144 vtot = 0.0;
145 for (i = 0; (i < nbonds); )
146 {
147 type = forceatoms[i++];
148 ai = forceatoms[i++];
149 aj = forceatoms[i++];
150
151 b0A = forceparams[type].morse.b0A;
152 beA = forceparams[type].morse.betaA;
153 cbA = forceparams[type].morse.cbA;
154
155 b0B = forceparams[type].morse.b0B;
156 beB = forceparams[type].morse.betaB;
157 cbB = forceparams[type].morse.cbB;
158
159 L1 = one-lambda; /* 1 */
160 b0 = L1*b0A + lambda*b0B; /* 3 */
161 be = L1*beA + lambda*beB; /* 3 */
162 cb = L1*cbA + lambda*cbB; /* 3 */
163
164 ki = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /* 3 */
165 dr2 = iprod(dx, dx); /* 5 */
166 dr = dr2*gmx::invsqrt(dr2); /* 10 */
167 temp = std::exp(-be*(dr-b0)); /* 12 */
168
169 if (temp == one)
170 {
171 /* bonds are constrainted. This may _not_ include bond constraints if they are lambda dependent */
172 *dvdlambda += cbB-cbA;
173 continue;
174 }
175
176 omtemp = one-temp; /* 1 */
177 cbomtemp = cb*omtemp; /* 1 */
178 vbond = cbomtemp*omtemp; /* 1 */
179 fbond = -two*be*temp*cbomtemp*gmx::invsqrt(dr2); /* 9 */
180 vtot += vbond; /* 1 */
181
182 *dvdlambda += (cbB - cbA) * omtemp * omtemp - (2-2*omtemp)*omtemp * cb * ((b0B-b0A)*be - (beB-beA)*(dr-b0)); /* 15 */
183
184 if (g)
185 {
186 ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, aj), dt);
187 ki = IVEC2IS(dt);
188 }
189
190 for (m = 0; (m < DIM); m++) /* 15 */
191 {
192 fij = fbond*dx[m];
193 f[ai][m] += fij;
194 f[aj][m] -= fij;
195 fshift[ki][m] += fij;
196 fshift[CENTRAL][m] -= fij;
197 }
198 } /* 83 TOTAL */
199 return vtot;
200 }
201
202 //! \cond
203 real cubic_bonds(int nbonds,
204 const t_iatom forceatoms[], const t_iparams forceparams[],
205 const rvec x[], rvec4 f[], rvec fshift[],
206 const t_pbc *pbc, const t_graph *g,
207 real gmx_unused lambda, real gmx_unused *dvdlambda,
208 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
209 int gmx_unused *global_atom_index)
210 {
211 const real three = 3.0;
212 const real two = 2.0;
213 real kb, b0, kcub;
214 real dr, dr2, dist, kdist, kdist2, fbond, vbond, fij, vtot;
215 rvec dx;
216 int i, m, ki, type, ai, aj;
217 ivec dt;
218
219 vtot = 0.0;
220 for (i = 0; (i < nbonds); )
221 {
222 type = forceatoms[i++];
223 ai = forceatoms[i++];
224 aj = forceatoms[i++];
225
226 b0 = forceparams[type].cubic.b0;
227 kb = forceparams[type].cubic.kb;
228 kcub = forceparams[type].cubic.kcub;
229
230 ki = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /* 3 */
231 dr2 = iprod(dx, dx); /* 5 */
232
233 if (dr2 == 0.0)
234 {
235 continue;
236 }
237
238 dr = dr2*gmx::invsqrt(dr2); /* 10 */
239 dist = dr-b0;
240 kdist = kb*dist;
241 kdist2 = kdist*dist;
242
243 vbond = kdist2 + kcub*kdist2*dist;
244 fbond = -(two*kdist + three*kdist2*kcub)/dr;
245
246 vtot += vbond; /* 21 */
247
248 if (g)
249 {
250 ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, aj), dt);
251 ki = IVEC2IS(dt);
252 }
253 for (m = 0; (m < DIM); m++) /* 15 */
254 {
255 fij = fbond*dx[m];
256 f[ai][m] += fij;
257 f[aj][m] -= fij;
258 fshift[ki][m] += fij;
259 fshift[CENTRAL][m] -= fij;
260 }
261 } /* 54 TOTAL */
262 return vtot;
263 }
264
265 real FENE_bonds(int nbonds,
266 const t_iatom forceatoms[], const t_iparams forceparams[],
267 const rvec x[], rvec4 f[], rvec fshift[],
268 const t_pbc *pbc, const t_graph *g,
269 real gmx_unused lambda, real gmx_unused *dvdlambda,
270 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
271 int *global_atom_index)
272 {
273 const real half = 0.5;
274 const real one = 1.0;
275 real bm, kb;
276 real dr2, bm2, omdr2obm2, fbond, vbond, fij, vtot;
277 rvec dx;
278 int i, m, ki, type, ai, aj;
279 ivec dt;
280
281 vtot = 0.0;
282 for (i = 0; (i < nbonds); )
283 {
284 type = forceatoms[i++];
285 ai = forceatoms[i++];
286 aj = forceatoms[i++];
287
288 bm = forceparams[type].fene.bm;
289 kb = forceparams[type].fene.kb;
290
291 ki = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /* 3 */
292 dr2 = iprod(dx, dx); /* 5 */
293
294 if (dr2 == 0.0)
295 {
296 continue;
297 }
298
299 bm2 = bm*bm;
300
301 if (dr2 >= bm2)
302 {
303 gmx_fatal(FARGS,
304 "r^2 (%f) >= bm^2 (%f) in FENE bond between atoms %d and %d",
305 dr2, bm2,
306 glatnr(global_atom_index, ai),
307 glatnr(global_atom_index, aj));
308 }
309
310 omdr2obm2 = one - dr2/bm2;
311
312 vbond = -half*kb*bm2*std::log(omdr2obm2);
313 fbond = -kb/omdr2obm2;
314
315 vtot += vbond; /* 35 */
316
317 if (g)
318 {
319 ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, aj), dt);
320 ki = IVEC2IS(dt);
321 }
322 for (m = 0; (m < DIM); m++) /* 15 */
323 {
324 fij = fbond*dx[m];
325 f[ai][m] += fij;
326 f[aj][m] -= fij;
327 fshift[ki][m] += fij;
328 fshift[CENTRAL][m] -= fij;
329 }
330 } /* 58 TOTAL */
331 return vtot;
332 }
333
334 static real harmonic(real kA, real kB, real xA, real xB, real x, real lambda,
335 real *V, real *F)
336 {
337 const real half = 0.5;
338 real L1, kk, x0, dx, dx2;
339 real v, f, dvdlambda;
340
341 L1 = 1.0-lambda;
342 kk = L1*kA+lambda*kB;
343 x0 = L1*xA+lambda*xB;
344
345 dx = x-x0;
346 dx2 = dx*dx;
347
348 f = -kk*dx;
349 v = half*kk*dx2;
350 dvdlambda = half*(kB-kA)*dx2 + (xA-xB)*kk*dx;
351
352 *F = f;
353 *V = v;
354
355 return dvdlambda;
356
357 /* That was 19 flops */
358 }
359
360
361 real bonds(int nbonds,
362 const t_iatom forceatoms[], const t_iparams forceparams[],
363 const rvec x[], rvec4 f[], rvec fshift[],
364 const t_pbc *pbc, const t_graph *g,
365 real lambda, real *dvdlambda,
366 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
367 int gmx_unused *global_atom_index)
368 {
369 int i, m, ki, ai, aj, type;
370 real dr, dr2, fbond, vbond, fij, vtot;
371 rvec dx;
372 ivec dt;
373
374 vtot = 0.0;
375 for (i = 0; (i < nbonds); )
376 {
377 type = forceatoms[i++];
378 ai = forceatoms[i++];
379 aj = forceatoms[i++];
380
381 ki = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /* 3 */
382 dr2 = iprod(dx, dx); /* 5 */
383 dr = std::sqrt(dr2); /* 10 */
384
385 *dvdlambda += harmonic(forceparams[type].harmonic.krA,
386 forceparams[type].harmonic.krB,
387 forceparams[type].harmonic.rA,
388 forceparams[type].harmonic.rB,
389 dr, lambda, &vbond, &fbond); /* 19 */
390
391 if (dr2 == 0.0)
392 {
393 continue;
394 }
395
396
397 vtot += vbond; /* 1*/
398 fbond *= gmx::invsqrt(dr2); /* 6 */
399 if (g)
400 {
401 ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, aj), dt);
402 ki = IVEC2IS(dt);
403 }
404 for (m = 0; (m < DIM); m++) /* 15 */
405 {
406 fij = fbond*dx[m];
407 f[ai][m] += fij;
408 f[aj][m] -= fij;
409 fshift[ki][m] += fij;
410 fshift[CENTRAL][m] -= fij;
411 }
412 } /* 59 TOTAL */
413 return vtot;
414 }
415
416 real restraint_bonds(int nbonds,
417 const t_iatom forceatoms[], const t_iparams forceparams[],
418 const rvec x[], rvec4 f[], rvec fshift[],
419 const t_pbc *pbc, const t_graph *g,
420 real lambda, real *dvdlambda,
421 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
422 int gmx_unused *global_atom_index)
423 {
424 int i, m, ki, ai, aj, type;
425 real dr, dr2, fbond, vbond, fij, vtot;
426 real L1;
427 real low, dlow, up1, dup1, up2, dup2, k, dk;
428 real drh, drh2;
429 rvec dx;
430 ivec dt;
431
432 L1 = 1.0 - lambda;
433
434 vtot = 0.0;
435 for (i = 0; (i < nbonds); )
436 {
437 type = forceatoms[i++];
438 ai = forceatoms[i++];
439 aj = forceatoms[i++];
440
441 ki = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /* 3 */
442 dr2 = iprod(dx, dx); /* 5 */
443 dr = dr2*gmx::invsqrt(dr2); /* 10 */
444
445 low = L1*forceparams[type].restraint.lowA + lambda*forceparams[type].restraint.lowB;
446 dlow = -forceparams[type].restraint.lowA + forceparams[type].restraint.lowB;
447 up1 = L1*forceparams[type].restraint.up1A + lambda*forceparams[type].restraint.up1B;
448 dup1 = -forceparams[type].restraint.up1A + forceparams[type].restraint.up1B;
449 up2 = L1*forceparams[type].restraint.up2A + lambda*forceparams[type].restraint.up2B;
450 dup2 = -forceparams[type].restraint.up2A + forceparams[type].restraint.up2B;
451 k = L1*forceparams[type].restraint.kA + lambda*forceparams[type].restraint.kB;
452 dk = -forceparams[type].restraint.kA + forceparams[type].restraint.kB;
453 /* 24 */
454
455 if (dr < low)
456 {
457 drh = dr - low;
458 drh2 = drh*drh;
459 vbond = 0.5*k*drh2;
460 fbond = -k*drh;
461 *dvdlambda += 0.5*dk*drh2 - k*dlow*drh;
462 } /* 11 */
463 else if (dr <= up1)
464 {
465 vbond = 0;
466 fbond = 0;
467 }
468 else if (dr <= up2)
469 {
470 drh = dr - up1;
471 drh2 = drh*drh;
472 vbond = 0.5*k*drh2;
473 fbond = -k*drh;
474 *dvdlambda += 0.5*dk*drh2 - k*dup1*drh;
475 } /* 11 */
476 else
477 {
478 drh = dr - up2;
479 vbond = k*(up2 - up1)*(0.5*(up2 - up1) + drh);
480 fbond = -k*(up2 - up1);
481 *dvdlambda += dk*(up2 - up1)*(0.5*(up2 - up1) + drh)
482 + k*(dup2 - dup1)*(up2 - up1 + drh)
483 - k*(up2 - up1)*dup2;
484 }
485
486 if (dr2 == 0.0)
487 {
488 continue;
489 }
490
491 vtot += vbond; /* 1*/
492 fbond *= gmx::invsqrt(dr2); /* 6 */
493 if (g)
494 {
495 ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, aj), dt);
496 ki = IVEC2IS(dt);
497 }
498 for (m = 0; (m < DIM); m++) /* 15 */
499 {
500 fij = fbond*dx[m];
501 f[ai][m] += fij;
502 f[aj][m] -= fij;
503 fshift[ki][m] += fij;
504 fshift[CENTRAL][m] -= fij;
505 }
506 } /* 59 TOTAL */
507
508 return vtot;
509 }
510
511 real polarize(int nbonds,
512 const t_iatom forceatoms[], const t_iparams forceparams[],
513 const rvec x[], rvec4 f[], rvec fshift[],
514 const t_pbc *pbc, const t_graph *g,
515 real lambda, real *dvdlambda,
516 const t_mdatoms *md, t_fcdata gmx_unused *fcd,
517 int gmx_unused *global_atom_index)
518 {
519 int i, m, ki, ai, aj, type;
520 real dr, dr2, fbond, vbond, fij, vtot, ksh;
521 rvec dx;
522 ivec dt;
523
524 vtot = 0.0;
525 for (i = 0; (i < nbonds); )
526 {
527 type = forceatoms[i++];
528 ai = forceatoms[i++];
529 aj = forceatoms[i++];
530 ksh = gmx::square(md->chargeA[aj])*ONE_4PI_EPS0/forceparams[type].polarize.alpha;
531
532 ki = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /* 3 */
533 dr2 = iprod(dx, dx); /* 5 */
534 dr = std::sqrt(dr2); /* 10 */
535
536 *dvdlambda += harmonic(ksh, ksh, 0, 0, dr, lambda, &vbond, &fbond); /* 19 */
537
538 if (dr2 == 0.0)
539 {
540 continue;
541 }
542
543 vtot += vbond; /* 1*/
544 fbond *= gmx::invsqrt(dr2); /* 6 */
545
546 if (g)
547 {
548 ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, aj), dt);
549 ki = IVEC2IS(dt);
550 }
551 for (m = 0; (m < DIM); m++) /* 15 */
552 {
553 fij = fbond*dx[m];
554 f[ai][m] += fij;
555 f[aj][m] -= fij;
556 fshift[ki][m] += fij;
557 fshift[CENTRAL][m] -= fij;
558 }
559 } /* 59 TOTAL */
560 return vtot;
561 }
562
563 real anharm_polarize(int nbonds,
564 const t_iatom forceatoms[], const t_iparams forceparams[],
565 const rvec x[], rvec4 f[], rvec fshift[],
566 const t_pbc *pbc, const t_graph *g,
567 real lambda, real *dvdlambda,
568 const t_mdatoms *md, t_fcdata gmx_unused *fcd,
569 int gmx_unused *global_atom_index)
570 {
571 int i, m, ki, ai, aj, type;
572 real dr, dr2, fbond, vbond, fij, vtot, ksh, khyp, drcut, ddr, ddr3;
573 rvec dx;
574 ivec dt;
575
576 vtot = 0.0;
577 for (i = 0; (i < nbonds); )
578 {
579 type = forceatoms[i++];
580 ai = forceatoms[i++];
581 aj = forceatoms[i++];
582 ksh = gmx::square(md->chargeA[aj])*ONE_4PI_EPS0/forceparams[type].anharm_polarize.alpha; /* 7*/
583 khyp = forceparams[type].anharm_polarize.khyp;
584 drcut = forceparams[type].anharm_polarize.drcut;
585
586 ki = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /* 3 */
587 dr2 = iprod(dx, dx); /* 5 */
588 dr = dr2*gmx::invsqrt(dr2); /* 10 */
589
590 *dvdlambda += harmonic(ksh, ksh, 0, 0, dr, lambda, &vbond, &fbond); /* 19 */
591
592 if (dr2 == 0.0)
593 {
594 continue;
595 }
596
597 if (dr > drcut)
598 {
599 ddr = dr-drcut;
600 ddr3 = ddr*ddr*ddr;
601 vbond += khyp*ddr*ddr3;
602 fbond -= 4*khyp*ddr3;
603 }
604 fbond *= gmx::invsqrt(dr2); /* 6 */
605 vtot += vbond; /* 1*/
606
607 if (g)
608 {
609 ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, aj), dt);
610 ki = IVEC2IS(dt);
611 }
612 for (m = 0; (m < DIM); m++) /* 15 */
613 {
614 fij = fbond*dx[m];
615 f[ai][m] += fij;
616 f[aj][m] -= fij;
617 fshift[ki][m] += fij;
618 fshift[CENTRAL][m] -= fij;
619 }
620 } /* 72 TOTAL */
621 return vtot;
622 }
623
624 real water_pol(int nbonds,
625 const t_iatom forceatoms[], const t_iparams forceparams[],
626 const rvec x[], rvec4 f[], rvec gmx_unused fshift[],
627 const t_pbc gmx_unused *pbc, const t_graph gmx_unused *g,
628 real gmx_unused lambda, real gmx_unused *dvdlambda,
629 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
630 int gmx_unused *global_atom_index)
631 {
632 /* This routine implements anisotropic polarizibility for water, through
633 * a shell connected to a dummy with spring constant that differ in the
634 * three spatial dimensions in the molecular frame.
635 */
636 int i, m, aO, aH1, aH2, aD, aS, type, type0, ki;
637 ivec dt;
638 rvec dOH1, dOH2, dHH, dOD, dDS, nW, kk, dx, kdx, proj;
639 real vtot, fij, r_HH, r_OD, r_nW, tx, ty, tz, qS;
640
641 vtot = 0.0;
642 if (nbonds > 0)
643 {
644 type0 = forceatoms[0];
645 aS = forceatoms[5];
646 qS = md->chargeA[aS];
647 kk[XX] = gmx::square(qS)*ONE_4PI_EPS0/forceparams[type0].wpol.al_x;
648 kk[YY] = gmx::square(qS)*ONE_4PI_EPS0/forceparams[type0].wpol.al_y;
649 kk[ZZ] = gmx::square(qS)*ONE_4PI_EPS0/forceparams[type0].wpol.al_z;
650 r_HH = 1.0/forceparams[type0].wpol.rHH;
651 for (i = 0; (i < nbonds); i += 6)
652 {
653 type = forceatoms[i];
654 if (type != type0)
655 {
656 gmx_fatal(FARGS, "Sorry, type = %d, type0 = %d, file = %s, line = %d",
657 type, type0, __FILE__, __LINE__);
658 }
659 aO = forceatoms[i+1];
660 aH1 = forceatoms[i+2];
661 aH2 = forceatoms[i+3];
662 aD = forceatoms[i+4];
663 aS = forceatoms[i+5];
664
665 /* Compute vectors describing the water frame */
666 pbc_rvec_sub(pbc, x[aH1], x[aO], dOH1);
667 pbc_rvec_sub(pbc, x[aH2], x[aO], dOH2);
668 pbc_rvec_sub(pbc, x[aH2], x[aH1], dHH);
669 pbc_rvec_sub(pbc, x[aD], x[aO], dOD);
670 ki = pbc_rvec_sub(pbc, x[aS], x[aD], dDS);
671 cprod(dOH1, dOH2, nW);
672
673 /* Compute inverse length of normal vector
674 * (this one could be precomputed, but I'm too lazy now)
675 */
676 r_nW = gmx::invsqrt(iprod(nW, nW));
677 /* This is for precision, but does not make a big difference,
678 * it can go later.
679 */
680 r_OD = gmx::invsqrt(iprod(dOD, dOD));
681
682 /* Normalize the vectors in the water frame */
683 svmul(r_nW, nW, nW);
684 svmul(r_HH, dHH, dHH);
685 svmul(r_OD, dOD, dOD);
686
687 /* Compute displacement of shell along components of the vector */
688 dx[ZZ] = iprod(dDS, dOD);
689 /* Compute projection on the XY plane: dDS - dx[ZZ]*dOD */
690 for (m = 0; (m < DIM); m++)
691 {
692 proj[m] = dDS[m]-dx[ZZ]*dOD[m];
693 }
694
695 /*dx[XX] = iprod(dDS,nW);
696 dx[YY] = iprod(dDS,dHH);*/
697 dx[XX] = iprod(proj, nW);
698 for (m = 0; (m < DIM); m++)
699 {
700 proj[m] -= dx[XX]*nW[m];
701 }
702 dx[YY] = iprod(proj, dHH);
703 /* Now compute the forces and energy */
704 kdx[XX] = kk[XX]*dx[XX];
705 kdx[YY] = kk[YY]*dx[YY];
706 kdx[ZZ] = kk[ZZ]*dx[ZZ];
707 vtot += iprod(dx, kdx);
708
709 if (g)
710 {
711 ivec_sub(SHIFT_IVEC(g, aS), SHIFT_IVEC(g, aD), dt);
712 ki = IVEC2IS(dt);
713 }
714
715 for (m = 0; (m < DIM); m++)
716 {
717 /* This is a tensor operation but written out for speed */
718 tx = nW[m]*kdx[XX];
719 ty = dHH[m]*kdx[YY];
720 tz = dOD[m]*kdx[ZZ];
721 fij = -tx-ty-tz;
722 f[aS][m] += fij;
723 f[aD][m] -= fij;
724 fshift[ki][m] += fij;
725 fshift[CENTRAL][m] -= fij;
726 }
727 }
728 }
729 return 0.5*vtot;
730 }
731
732 static real do_1_thole(const rvec xi, const rvec xj, rvec fi, rvec fj,
733 const t_pbc *pbc, real qq,
734 rvec fshift[], real afac)
735 {
736 rvec r12;
737 real r12sq, r12_1, r12bar, v0, v1, fscal, ebar, fff;
738 int m, t;
739
740 t = pbc_rvec_sub(pbc, xi, xj, r12); /* 3 */
741
742 r12sq = iprod(r12, r12); /* 5 */
743 r12_1 = gmx::invsqrt(r12sq); /* 5 */
744 r12bar = afac/r12_1; /* 5 */
745 v0 = qq*ONE_4PI_EPS0*r12_1; /* 2 */
746 ebar = std::exp(-r12bar); /* 5 */
747 v1 = (1-(1+0.5*r12bar)*ebar); /* 4 */
748 fscal = ((v0*r12_1)*v1 - v0*0.5*afac*ebar*(r12bar+1))*r12_1; /* 9 */
749
750 for (m = 0; (m < DIM); m++)
751 {
752 fff = fscal*r12[m];
753 fi[m] += fff;
754 fj[m] -= fff;
755 fshift[t][m] += fff;
756 fshift[CENTRAL][m] -= fff;
757 } /* 15 */
758
759 return v0*v1; /* 1 */
760 /* 54 */
761 }
762
763 real thole_pol(int nbonds,
764 const t_iatom forceatoms[], const t_iparams forceparams[],
765 const rvec x[], rvec4 f[], rvec fshift[],
766 const t_pbc *pbc, const t_graph gmx_unused *g,
767 real gmx_unused lambda, real gmx_unused *dvdlambda,
768 const t_mdatoms *md, t_fcdata gmx_unused *fcd,
769 int gmx_unused *global_atom_index)
770 {
771 /* Interaction between two pairs of particles with opposite charge */
772 int i, type, a1, da1, a2, da2;
773 real q1, q2, qq, a, al1, al2, afac;
774 real V = 0;
775
776 for (i = 0; (i < nbonds); )
777 {
778 type = forceatoms[i++];
779 a1 = forceatoms[i++];
780 da1 = forceatoms[i++];
781 a2 = forceatoms[i++];
782 da2 = forceatoms[i++];
783 q1 = md->chargeA[da1];
784 q2 = md->chargeA[da2];
785 a = forceparams[type].thole.a;
786 al1 = forceparams[type].thole.alpha1;
787 al2 = forceparams[type].thole.alpha2;
788 qq = q1*q2;
789 afac = a*gmx::invsixthroot(al1*al2);
790 V += do_1_thole(x[a1], x[a2], f[a1], f[a2], pbc, qq, fshift, afac);
791 V += do_1_thole(x[da1], x[a2], f[da1], f[a2], pbc, -qq, fshift, afac);
792 V += do_1_thole(x[a1], x[da2], f[a1], f[da2], pbc, -qq, fshift, afac);
793 V += do_1_thole(x[da1], x[da2], f[da1], f[da2], pbc, qq, fshift, afac);
794 }
795 /* 290 flops */
796 return V;
797 }
798
799 real bond_angle(const rvec xi, const rvec xj, const rvec xk, const t_pbc *pbc,
800 rvec r_ij, rvec r_kj, real *costh,
801 int *t1, int *t2)
802 /* Return value is the angle between the bonds i-j and j-k */
803 {
804 /* 41 FLOPS */
805 real th;
806
807 *t1 = pbc_rvec_sub(pbc, xi, xj, r_ij); /* 3 */
808 *t2 = pbc_rvec_sub(pbc, xk, xj, r_kj); /* 3 */
809
810 *costh = cos_angle(r_ij, r_kj); /* 25 */
811 th = std::acos(*costh); /* 10 */
812 /* 41 TOTAL */
813 return th;
814 }
815
816 real angles(int nbonds,
817 const t_iatom forceatoms[], const t_iparams forceparams[],
818 const rvec x[], rvec4 f[], rvec fshift[],
819 const t_pbc *pbc, const t_graph *g,
820 real lambda, real *dvdlambda,
821 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
822 int gmx_unused *global_atom_index)
823 {
824 int i, ai, aj, ak, t1, t2, type;
825 rvec r_ij, r_kj;
826 real cos_theta, cos_theta2, theta, dVdt, va, vtot;
827 ivec jt, dt_ij, dt_kj;
828
829 vtot = 0.0;
830 for (i = 0; i < nbonds; )
831 {
832 type = forceatoms[i++];
833 ai = forceatoms[i++];
834 aj = forceatoms[i++];
835 ak = forceatoms[i++];
836
837 theta = bond_angle(x[ai], x[aj], x[ak], pbc,
838 r_ij, r_kj, &cos_theta, &t1, &t2); /* 41 */
839
840 *dvdlambda += harmonic(forceparams[type].harmonic.krA,
841 forceparams[type].harmonic.krB,
842 forceparams[type].harmonic.rA*DEG2RAD,
843 forceparams[type].harmonic.rB*DEG2RAD,
844 theta, lambda, &va, &dVdt); /* 21 */
845 vtot += va;
846
847 cos_theta2 = gmx::square(cos_theta);
848 if (cos_theta2 < 1)
849 {
850 int m;
851 real st, sth;
852 real cik, cii, ckk;
853 real nrkj2, nrij2;
854 real nrkj_1, nrij_1;
855 rvec f_i, f_j, f_k;
856
857 st = dVdt*gmx::invsqrt(1 - cos_theta2); /* 12 */
858 sth = st*cos_theta; /* 1 */
859 nrij2 = iprod(r_ij, r_ij); /* 5 */
860 nrkj2 = iprod(r_kj, r_kj); /* 5 */
861
862 nrij_1 = gmx::invsqrt(nrij2); /* 10 */
863 nrkj_1 = gmx::invsqrt(nrkj2); /* 10 */
864
865 cik = st*nrij_1*nrkj_1; /* 2 */
866 cii = sth*nrij_1*nrij_1; /* 2 */
867 ckk = sth*nrkj_1*nrkj_1; /* 2 */
868
869 for (m = 0; m < DIM; m++)
870 { /* 39 */
871 f_i[m] = -(cik*r_kj[m] - cii*r_ij[m]);
872 f_k[m] = -(cik*r_ij[m] - ckk*r_kj[m]);
873 f_j[m] = -f_i[m] - f_k[m];
874 f[ai][m] += f_i[m];
875 f[aj][m] += f_j[m];
876 f[ak][m] += f_k[m];
877 }
878 if (g != nullptr)
879 {
880 copy_ivec(SHIFT_IVEC(g, aj), jt);
881
882 ivec_sub(SHIFT_IVEC(g, ai), jt, dt_ij);
883 ivec_sub(SHIFT_IVEC(g, ak), jt, dt_kj);
884 t1 = IVEC2IS(dt_ij);
885 t2 = IVEC2IS(dt_kj);
886 }
887 rvec_inc(fshift[t1], f_i);
888 rvec_inc(fshift[CENTRAL], f_j);
889 rvec_inc(fshift[t2], f_k);
890 } /* 161 TOTAL */
891 }
892
893 return vtot;
894 }
895
896 #if GMX_SIMD_HAVE_REAL
897
898 /* As angles, but using SIMD to calculate many angles at once.
899 * This routines does not calculate energies and shift forces.
900 */
901 void
902 angles_noener_simd(int nbonds,
903 const t_iatom forceatoms[], const t_iparams forceparams[],
904 const rvec x[], rvec4 f[],
905 const t_pbc *pbc, const t_graph gmx_unused *g,
906 real gmx_unused lambda,
907 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
908 int gmx_unused *global_atom_index)
909 {
910 const int nfa1 = 4;
911 int i, iu, s;
912 int type;
913 alignas(GMX_SIMD_ALIGNMENT) std::int32_t ai[GMX_SIMD_REAL_WIDTH];
914 alignas(GMX_SIMD_ALIGNMENT) std::int32_t aj[GMX_SIMD_REAL_WIDTH];
915 alignas(GMX_SIMD_ALIGNMENT) std::int32_t ak[GMX_SIMD_REAL_WIDTH];
916 alignas(GMX_SIMD_ALIGNMENT) real coeff[2*GMX_SIMD_REAL_WIDTH];
917 SimdReal deg2rad_S(DEG2RAD);
918 SimdReal xi_S, yi_S, zi_S;
919 SimdReal xj_S, yj_S, zj_S;
920 SimdReal xk_S, yk_S, zk_S;
921 SimdReal k_S, theta0_S;
922 SimdReal rijx_S, rijy_S, rijz_S;
923 SimdReal rkjx_S, rkjy_S, rkjz_S;
924 SimdReal one_S(1.0);
925 SimdReal min_one_plus_eps_S(-1.0 + 2.0*GMX_REAL_EPS); // Smallest number > -1
926
927 SimdReal rij_rkj_S;
928 SimdReal nrij2_S, nrij_1_S;
929 SimdReal nrkj2_S, nrkj_1_S;
930 SimdReal cos_S, invsin_S;
931 SimdReal theta_S;
932 SimdReal st_S, sth_S;
933 SimdReal cik_S, cii_S, ckk_S;
934 SimdReal f_ix_S, f_iy_S, f_iz_S;
935 SimdReal f_kx_S, f_ky_S, f_kz_S;
936 alignas(GMX_SIMD_ALIGNMENT) real pbc_simd[9*GMX_SIMD_REAL_WIDTH];
937
938 set_pbc_simd(pbc, pbc_simd);
939
940 /* nbonds is the number of angles times nfa1, here we step GMX_SIMD_REAL_WIDTH angles */
941 for (i = 0; (i < nbonds); i += GMX_SIMD_REAL_WIDTH*nfa1)
942 {
943 /* Collect atoms for GMX_SIMD_REAL_WIDTH angles.
944 * iu indexes into forceatoms, we should not let iu go beyond nbonds.
945 */
946 iu = i;
947 for (s = 0; s < GMX_SIMD_REAL_WIDTH; s++)
948 {
949 type = forceatoms[iu];
950 ai[s] = forceatoms[iu+1];
951 aj[s] = forceatoms[iu+2];
952 ak[s] = forceatoms[iu+3];
953
954 /* At the end fill the arrays with the last atoms and 0 params */
955 if (i + s*nfa1 < nbonds)
956 {
957 coeff[s] = forceparams[type].harmonic.krA;
958 coeff[GMX_SIMD_REAL_WIDTH+s] = forceparams[type].harmonic.rA;
959
960 if (iu + nfa1 < nbonds)
961 {
962 iu += nfa1;
963 }
964 }
965 else
966 {
967 coeff[s] = 0;
968 coeff[GMX_SIMD_REAL_WIDTH+s] = 0;
969 }
970 }
971
972 /* Store the non PBC corrected distances packed and aligned */
973 gatherLoadUTranspose<3>(reinterpret_cast<const real *>(x), ai, &xi_S, &yi_S, &zi_S);
974 gatherLoadUTranspose<3>(reinterpret_cast<const real *>(x), aj, &xj_S, &yj_S, &zj_S);
975 gatherLoadUTranspose<3>(reinterpret_cast<const real *>(x), ak, &xk_S, &yk_S, &zk_S);
976 rijx_S = xi_S - xj_S;
977 rijy_S = yi_S - yj_S;
978 rijz_S = zi_S - zj_S;
979 rkjx_S = xk_S - xj_S;
980 rkjy_S = yk_S - yj_S;
981 rkjz_S = zk_S - zj_S;
982
983 k_S = load<SimdReal>(coeff);
984 theta0_S = load<SimdReal>(coeff+GMX_SIMD_REAL_WIDTH) * deg2rad_S;
985
986 pbc_correct_dx_simd(&rijx_S, &rijy_S, &rijz_S, pbc_simd);
987 pbc_correct_dx_simd(&rkjx_S, &rkjy_S, &rkjz_S, pbc_simd);
988
989 rij_rkj_S = iprod(rijx_S, rijy_S, rijz_S,
990 rkjx_S, rkjy_S, rkjz_S);
991
992 nrij2_S = norm2(rijx_S, rijy_S, rijz_S);
993 nrkj2_S = norm2(rkjx_S, rkjy_S, rkjz_S);
994
995 nrij_1_S = invsqrt(nrij2_S);
996 nrkj_1_S = invsqrt(nrkj2_S);
997
998 cos_S = rij_rkj_S * nrij_1_S * nrkj_1_S;
999
1000 /* To allow for 180 degrees, we take the max of cos and -1 + 1bit,
1001 * so we can safely get the 1/sin from 1/sqrt(1 - cos^2).
1002 * This also ensures that rounding errors would cause the argument
1003 * of simdAcos to be < -1.
1004 * Note that we do not take precautions for cos(0)=1, so the outer
1005 * atoms in an angle should not be on top of each other.
1006 */
1007 cos_S = max(cos_S, min_one_plus_eps_S);
1008
1009 theta_S = acos(cos_S);
1010
1011 invsin_S = invsqrt( one_S - cos_S * cos_S );
1012
1013 st_S = k_S * (theta0_S - theta_S) * invsin_S;
1014 sth_S = st_S * cos_S;
1015
1016 cik_S = st_S * nrij_1_S * nrkj_1_S;
1017 cii_S = sth_S * nrij_1_S * nrij_1_S;
1018 ckk_S = sth_S * nrkj_1_S * nrkj_1_S;
1019
1020 f_ix_S = cii_S * rijx_S;
1021 f_ix_S = fnma(cik_S, rkjx_S, f_ix_S);
1022 f_iy_S = cii_S * rijy_S;
1023 f_iy_S = fnma(cik_S, rkjy_S, f_iy_S);
1024 f_iz_S = cii_S * rijz_S;
1025 f_iz_S = fnma(cik_S, rkjz_S, f_iz_S);
1026 f_kx_S = ckk_S * rkjx_S;
1027 f_kx_S = fnma(cik_S, rijx_S, f_kx_S);
1028 f_ky_S = ckk_S * rkjy_S;
1029 f_ky_S = fnma(cik_S, rijy_S, f_ky_S);
1030 f_kz_S = ckk_S * rkjz_S;
1031 f_kz_S = fnma(cik_S, rijz_S, f_kz_S);
1032
1033 transposeScatterIncrU<4>(reinterpret_cast<real *>(f), ai, f_ix_S, f_iy_S, f_iz_S);
1034 transposeScatterDecrU<4>(reinterpret_cast<real *>(f), aj, f_ix_S + f_kx_S, f_iy_S + f_ky_S, f_iz_S + f_kz_S);
1035 transposeScatterIncrU<4>(reinterpret_cast<real *>(f), ak, f_kx_S, f_ky_S, f_kz_S);
1036 }
1037 }
1038
1039 #endif // GMX_SIMD_HAVE_REAL
1040
1041 real linear_angles(int nbonds,
1042 const t_iatom forceatoms[], const t_iparams forceparams[],
1043 const rvec x[], rvec4 f[], rvec fshift[],
1044 const t_pbc *pbc, const t_graph *g,
1045 real lambda, real *dvdlambda,
1046 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
1047 int gmx_unused *global_atom_index)
1048 {
1049 int i, m, ai, aj, ak, t1, t2, type;
1050 rvec f_i, f_j, f_k;
1051 real L1, kA, kB, aA, aB, dr, dr2, va, vtot, a, b, klin;
1052 ivec jt, dt_ij, dt_kj;
1053 rvec r_ij, r_kj, r_ik, dx;
1054
1055 L1 = 1-lambda;
1056 vtot = 0.0;
1057 for (i = 0; (i < nbonds); )
1058 {
1059 type = forceatoms[i++];
1060 ai = forceatoms[i++];
1061 aj = forceatoms[i++];
1062 ak = forceatoms[i++];
1063
1064 kA = forceparams[type].linangle.klinA;
1065 kB = forceparams[type].linangle.klinB;
1066 klin = L1*kA + lambda*kB;
1067
1068 aA = forceparams[type].linangle.aA;
1069 aB = forceparams[type].linangle.aB;
1070 a = L1*aA+lambda*aB;
1071 b = 1-a;
1072
1073 t1 = pbc_rvec_sub(pbc, x[ai], x[aj], r_ij);
1074 t2 = pbc_rvec_sub(pbc, x[ak], x[aj], r_kj);
1075 rvec_sub(r_ij, r_kj, r_ik);
1076
1077 dr2 = 0;
1078 for (m = 0; (m < DIM); m++)
1079 {
1080 dr = -a * r_ij[m] - b * r_kj[m];
1081 dr2 += dr*dr;
1082 dx[m] = dr;
1083 f_i[m] = a*klin*dr;
1084 f_k[m] = b*klin*dr;
1085 f_j[m] = -(f_i[m]+f_k[m]);
1086 f[ai][m] += f_i[m];
1087 f[aj][m] += f_j[m];
1088 f[ak][m] += f_k[m];
1089 }
1090 va = 0.5*klin*dr2;
1091 *dvdlambda += 0.5*(kB-kA)*dr2 + klin*(aB-aA)*iprod(dx, r_ik);
1092
1093 vtot += va;
1094
1095 if (g)
1096 {
1097 copy_ivec(SHIFT_IVEC(g, aj), jt);
1098
1099 ivec_sub(SHIFT_IVEC(g, ai), jt, dt_ij);
1100 ivec_sub(SHIFT_IVEC(g, ak), jt, dt_kj);
1101 t1 = IVEC2IS(dt_ij);
1102 t2 = IVEC2IS(dt_kj);
1103 }
1104 rvec_inc(fshift[t1], f_i);
1105 rvec_inc(fshift[CENTRAL], f_j);
1106 rvec_inc(fshift[t2], f_k);
1107 } /* 57 TOTAL */
1108 return vtot;
1109 }
1110
1111 real urey_bradley(int nbonds,
1112 const t_iatom forceatoms[], const t_iparams forceparams[],
1113 const rvec x[], rvec4 f[], rvec fshift[],
1114 const t_pbc *pbc, const t_graph *g,
1115 real lambda, real *dvdlambda,
1116 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
1117 int gmx_unused *global_atom_index)
1118 {
1119 int i, m, ai, aj, ak, t1, t2, type, ki;
1120 rvec r_ij, r_kj, r_ik;
1121 real cos_theta, cos_theta2, theta;
1122 real dVdt, va, vtot, dr, dr2, vbond, fbond, fik;
1123 real kthA, th0A, kUBA, r13A, kthB, th0B, kUBB, r13B;
1124 ivec jt, dt_ij, dt_kj, dt_ik;
1125
1126 vtot = 0.0;
1127 for (i = 0; (i < nbonds); )
1128 {
1129 type = forceatoms[i++];
1130 ai = forceatoms[i++];
1131 aj = forceatoms[i++];
1132 ak = forceatoms[i++];
1133 th0A = forceparams[type].u_b.thetaA*DEG2RAD;
1134 kthA = forceparams[type].u_b.kthetaA;
1135 r13A = forceparams[type].u_b.r13A;
1136 kUBA = forceparams[type].u_b.kUBA;
1137 th0B = forceparams[type].u_b.thetaB*DEG2RAD;
1138 kthB = forceparams[type].u_b.kthetaB;
1139 r13B = forceparams[type].u_b.r13B;
1140 kUBB = forceparams[type].u_b.kUBB;
1141
1142 theta = bond_angle(x[ai], x[aj], x[ak], pbc,
1143 r_ij, r_kj, &cos_theta, &t1, &t2); /* 41 */
1144
1145 *dvdlambda += harmonic(kthA, kthB, th0A, th0B, theta, lambda, &va, &dVdt); /* 21 */
1146 vtot += va;
1147
1148 ki = pbc_rvec_sub(pbc, x[ai], x[ak], r_ik); /* 3 */
1149 dr2 = iprod(r_ik, r_ik); /* 5 */
1150 dr = dr2*gmx::invsqrt(dr2); /* 10 */
1151
1152 *dvdlambda += harmonic(kUBA, kUBB, r13A, r13B, dr, lambda, &vbond, &fbond); /* 19 */
1153
1154 cos_theta2 = gmx::square(cos_theta); /* 1 */
1155 if (cos_theta2 < 1)
1156 {
1157 real st, sth;
1158 real cik, cii, ckk;
1159 real nrkj2, nrij2;
1160 rvec f_i, f_j, f_k;
1161
1162 st = dVdt*gmx::invsqrt(1 - cos_theta2); /* 12 */
1163 sth = st*cos_theta; /* 1 */
1164 nrkj2 = iprod(r_kj, r_kj); /* 5 */
1165 nrij2 = iprod(r_ij, r_ij);
1166
1167 cik = st*gmx::invsqrt(nrkj2*nrij2); /* 12 */
1168 cii = sth/nrij2; /* 10 */
1169 ckk = sth/nrkj2; /* 10 */
1170
1171 for (m = 0; (m < DIM); m++) /* 39 */
1172 {
1173 f_i[m] = -(cik*r_kj[m]-cii*r_ij[m]);
1174 f_k[m] = -(cik*r_ij[m]-ckk*r_kj[m]);
1175 f_j[m] = -f_i[m]-f_k[m];
1176 f[ai][m] += f_i[m];
1177 f[aj][m] += f_j[m];
1178 f[ak][m] += f_k[m];
1179 }
1180 if (g)
1181 {
1182 copy_ivec(SHIFT_IVEC(g, aj), jt);
1183
1184 ivec_sub(SHIFT_IVEC(g, ai), jt, dt_ij);
1185 ivec_sub(SHIFT_IVEC(g, ak), jt, dt_kj);
1186 t1 = IVEC2IS(dt_ij);
1187 t2 = IVEC2IS(dt_kj);
1188 }
1189 rvec_inc(fshift[t1], f_i);
1190 rvec_inc(fshift[CENTRAL], f_j);
1191 rvec_inc(fshift[t2], f_k);
1192 } /* 161 TOTAL */
1193 /* Time for the bond calculations */
1194 if (dr2 == 0.0)
1195 {
1196 continue;
1197 }
1198
1199 vtot += vbond; /* 1*/
1200 fbond *= gmx::invsqrt(dr2); /* 6 */
1201
1202 if (g)
1203 {
1204 ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, ak), dt_ik);
1205 ki = IVEC2IS(dt_ik);
1206 }
1207 for (m = 0; (m < DIM); m++) /* 15 */
1208 {
1209 fik = fbond*r_ik[m];
1210 f[ai][m] += fik;
1211 f[ak][m] -= fik;
1212 fshift[ki][m] += fik;
1213 fshift[CENTRAL][m] -= fik;
1214 }
1215 }
1216 return vtot;
1217 }
1218
1219 #if GMX_SIMD_HAVE_REAL
1220
1221 /* As urey_bradley, but using SIMD to calculate many potentials at once.
1222 * This routines does not calculate energies and shift forces.
1223 */
1224 void urey_bradley_noener_simd(int nbonds,
1225 const t_iatom forceatoms[], const t_iparams forceparams[],
1226 const rvec x[], rvec4 f[],
1227 const t_pbc *pbc, const t_graph gmx_unused *g,
1228 real gmx_unused lambda,
1229 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
1230 int gmx_unused *global_atom_index)
1231 {
1232 constexpr int nfa1 = 4;
1233 alignas(GMX_SIMD_ALIGNMENT) std::int32_t ai[GMX_SIMD_REAL_WIDTH];
1234 alignas(GMX_SIMD_ALIGNMENT) std::int32_t aj[GMX_SIMD_REAL_WIDTH];
1235 alignas(GMX_SIMD_ALIGNMENT) std::int32_t ak[GMX_SIMD_REAL_WIDTH];
1236 alignas(GMX_SIMD_ALIGNMENT) real coeff[4*GMX_SIMD_REAL_WIDTH];
1237 alignas(GMX_SIMD_ALIGNMENT) real pbc_simd[9*GMX_SIMD_REAL_WIDTH];
1238
1239 set_pbc_simd(pbc, pbc_simd);
1240
1241 /* nbonds is the number of angles times nfa1, here we step GMX_SIMD_REAL_WIDTH angles */
1242 for (int i = 0; i < nbonds; i += GMX_SIMD_REAL_WIDTH*nfa1)
1243 {
1244 /* Collect atoms for GMX_SIMD_REAL_WIDTH angles.
1245 * iu indexes into forceatoms, we should not let iu go beyond nbonds.
1246 */
1247 int iu = i;
1248 for (int s = 0; s < GMX_SIMD_REAL_WIDTH; s++)
1249 {
1250 const int type = forceatoms[iu];
1251 ai[s] = forceatoms[iu+1];
1252 aj[s] = forceatoms[iu+2];
1253 ak[s] = forceatoms[iu+3];
1254
1255 /* At the end fill the arrays with the last atoms and 0 params */
1256 if (i + s*nfa1 < nbonds)
1257 {
1258 coeff[s] = forceparams[type].u_b.kthetaA;
1259 coeff[GMX_SIMD_REAL_WIDTH+s] = forceparams[type].u_b.thetaA;
1260 coeff[GMX_SIMD_REAL_WIDTH*2+s] = forceparams[type].u_b.kUBA;
1261 coeff[GMX_SIMD_REAL_WIDTH*3+s] = forceparams[type].u_b.r13A;
1262
1263 if (iu + nfa1 < nbonds)
1264 {
1265 iu += nfa1;
1266 }
1267 }
1268 else
1269 {
1270 coeff[s] = 0;
1271 coeff[GMX_SIMD_REAL_WIDTH+s] = 0;
1272 coeff[GMX_SIMD_REAL_WIDTH*2+s] = 0;
1273 coeff[GMX_SIMD_REAL_WIDTH*3+s] = 0;
1274 }
1275 }
1276
1277 SimdReal xi_S, yi_S, zi_S;
1278 SimdReal xj_S, yj_S, zj_S;
1279 SimdReal xk_S, yk_S, zk_S;
1280
1281 /* Store the non PBC corrected distances packed and aligned */
1282 gatherLoadUTranspose<3>(reinterpret_cast<const real *>(x), ai, &xi_S, &yi_S, &zi_S);
1283 gatherLoadUTranspose<3>(reinterpret_cast<const real *>(x), aj, &xj_S, &yj_S, &zj_S);
1284 gatherLoadUTranspose<3>(reinterpret_cast<const real *>(x), ak, &xk_S, &yk_S, &zk_S);
1285 SimdReal rijx_S = xi_S - xj_S;
1286 SimdReal rijy_S = yi_S - yj_S;
1287 SimdReal rijz_S = zi_S - zj_S;
1288 SimdReal rkjx_S = xk_S - xj_S;
1289 SimdReal rkjy_S = yk_S - yj_S;
1290 SimdReal rkjz_S = zk_S - zj_S;
1291 SimdReal rikx_S = xi_S - xk_S;
1292 SimdReal riky_S = yi_S - yk_S;
1293 SimdReal rikz_S = zi_S - zk_S;
1294
1295 const SimdReal ktheta_S = load<SimdReal>(coeff);
1296 const SimdReal theta0_S = load<SimdReal>(coeff+GMX_SIMD_REAL_WIDTH) * DEG2RAD;
1297 const SimdReal kUB_S = load<SimdReal>(coeff+2*GMX_SIMD_REAL_WIDTH);
1298 const SimdReal r13_S = load<SimdReal>(coeff+3*GMX_SIMD_REAL_WIDTH);
1299
1300 pbc_correct_dx_simd(&rijx_S, &rijy_S, &rijz_S, pbc_simd);
1301 pbc_correct_dx_simd(&rkjx_S, &rkjy_S, &rkjz_S, pbc_simd);
1302 pbc_correct_dx_simd(&rikx_S, &riky_S, &rikz_S, pbc_simd);
1303
1304 const SimdReal rij_rkj_S = iprod(rijx_S, rijy_S, rijz_S,
1305 rkjx_S, rkjy_S, rkjz_S);
1306
1307 const SimdReal dr2_S = iprod(rikx_S, riky_S, rikz_S,
1308 rikx_S, riky_S, rikz_S);
1309
1310 const SimdReal nrij2_S = norm2(rijx_S, rijy_S, rijz_S);
1311 const SimdReal nrkj2_S = norm2(rkjx_S, rkjy_S, rkjz_S);
1312
1313 const SimdReal nrij_1_S = invsqrt(nrij2_S);
1314 const SimdReal nrkj_1_S = invsqrt(nrkj2_S);
1315 const SimdReal invdr2_S = invsqrt(dr2_S);
1316 const SimdReal dr_S = dr2_S*invdr2_S;
1317
1318 constexpr real min_one_plus_eps = -1.0 + 2.0*GMX_REAL_EPS; // Smallest number > -1
1319
1320 /* To allow for 180 degrees, we take the max of cos and -1 + 1bit,
1321 * so we can safely get the 1/sin from 1/sqrt(1 - cos^2).
1322 * This also ensures that rounding errors would cause the argument
1323 * of simdAcos to be < -1.
1324 * Note that we do not take precautions for cos(0)=1, so the bonds
1325 * in an angle should not align at an angle of 0 degrees.
1326 */
1327 const SimdReal cos_S = max(rij_rkj_S * nrij_1_S * nrkj_1_S, min_one_plus_eps);
1328
1329 const SimdReal theta_S = acos(cos_S);
1330 const SimdReal invsin_S = invsqrt( 1.0 - cos_S * cos_S );
1331 const SimdReal st_S = ktheta_S * (theta0_S - theta_S) * invsin_S;
1332 const SimdReal sth_S = st_S * cos_S;
1333
1334 const SimdReal cik_S = st_S * nrij_1_S * nrkj_1_S;
1335 const SimdReal cii_S = sth_S * nrij_1_S * nrij_1_S;
1336 const SimdReal ckk_S = sth_S * nrkj_1_S * nrkj_1_S;
1337
1338 const SimdReal sUB_S = kUB_S * (r13_S - dr_S) * invdr2_S;
1339
1340 const SimdReal f_ikx_S = sUB_S * rikx_S;
1341 const SimdReal f_iky_S = sUB_S * riky_S;
1342 const SimdReal f_ikz_S = sUB_S * rikz_S;
1343
1344 const SimdReal f_ix_S = fnma(cik_S, rkjx_S, cii_S * rijx_S) + f_ikx_S;
1345 const SimdReal f_iy_S = fnma(cik_S, rkjy_S, cii_S * rijy_S) + f_iky_S;
1346 const SimdReal f_iz_S = fnma(cik_S, rkjz_S, cii_S * rijz_S) + f_ikz_S;
1347 const SimdReal f_kx_S = fnma(cik_S, rijx_S, ckk_S * rkjx_S) - f_ikx_S;
1348 const SimdReal f_ky_S = fnma(cik_S, rijy_S, ckk_S * rkjy_S) - f_iky_S;
1349 const SimdReal f_kz_S = fnma(cik_S, rijz_S, ckk_S * rkjz_S) - f_ikz_S;
1350
1351 transposeScatterIncrU<4>(reinterpret_cast<real *>(f), ai, f_ix_S, f_iy_S, f_iz_S);
1352 transposeScatterDecrU<4>(reinterpret_cast<real *>(f), aj, f_ix_S + f_kx_S, f_iy_S + f_ky_S, f_iz_S + f_kz_S);
1353 transposeScatterIncrU<4>(reinterpret_cast<real *>(f), ak, f_kx_S, f_ky_S, f_kz_S);
1354 }
1355 }
1356
1357 #endif // GMX_SIMD_HAVE_REAL
1358
1359 real quartic_angles(int nbonds,
1360 const t_iatom forceatoms[], const t_iparams forceparams[],
1361 const rvec x[], rvec4 f[], rvec fshift[],
1362 const t_pbc *pbc, const t_graph *g,
1363 real gmx_unused lambda, real gmx_unused *dvdlambda,
1364 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
1365 int gmx_unused *global_atom_index)
1366 {
1367 int i, j, ai, aj, ak, t1, t2, type;
1368 rvec r_ij, r_kj;
1369 real cos_theta, cos_theta2, theta, dt, dVdt, va, dtp, c, vtot;
1370 ivec jt, dt_ij, dt_kj;
1371
1372 vtot = 0.0;
1373 for (i = 0; (i < nbonds); )
1374 {
1375 type = forceatoms[i++];
1376 ai = forceatoms[i++];
1377 aj = forceatoms[i++];
1378 ak = forceatoms[i++];
1379
1380 theta = bond_angle(x[ai], x[aj], x[ak], pbc,
1381 r_ij, r_kj, &cos_theta, &t1, &t2); /* 41 */
1382
1383 dt = theta - forceparams[type].qangle.theta*DEG2RAD; /* 2 */
1384
1385 dVdt = 0;
1386 va = forceparams[type].qangle.c[0];
1387 dtp = 1.0;
1388 for (j = 1; j <= 4; j++)
1389 {
1390 c = forceparams[type].qangle.c[j];
1391 dVdt -= j*c*dtp;
1392 dtp *= dt;
1393 va += c*dtp;
1394 }
1395 /* 20 */
1396
1397 vtot += va;
1398
1399 cos_theta2 = gmx::square(cos_theta); /* 1 */
1400 if (cos_theta2 < 1)
1401 {
1402 int m;
1403 real st, sth;
1404 real cik, cii, ckk;
1405 real nrkj2, nrij2;
1406 rvec f_i, f_j, f_k;
1407
1408 st = dVdt*gmx::invsqrt(1 - cos_theta2); /* 12 */
1409 sth = st*cos_theta; /* 1 */
1410 nrkj2 = iprod(r_kj, r_kj); /* 5 */
1411 nrij2 = iprod(r_ij, r_ij);
1412
1413 cik = st*gmx::invsqrt(nrkj2*nrij2); /* 12 */
1414 cii = sth/nrij2; /* 10 */
1415 ckk = sth/nrkj2; /* 10 */
1416
1417 for (m = 0; (m < DIM); m++) /* 39 */
1418 {
1419 f_i[m] = -(cik*r_kj[m]-cii*r_ij[m]);
1420 f_k[m] = -(cik*r_ij[m]-ckk*r_kj[m]);
1421 f_j[m] = -f_i[m]-f_k[m];
1422 f[ai][m] += f_i[m];
1423 f[aj][m] += f_j[m];
1424 f[ak][m] += f_k[m];
1425 }
1426 if (g)
1427 {
1428 copy_ivec(SHIFT_IVEC(g, aj), jt);
1429
1430 ivec_sub(SHIFT_IVEC(g, ai), jt, dt_ij);
1431 ivec_sub(SHIFT_IVEC(g, ak), jt, dt_kj);
1432 t1 = IVEC2IS(dt_ij);
1433 t2 = IVEC2IS(dt_kj);
1434 }
1435 rvec_inc(fshift[t1], f_i);
1436 rvec_inc(fshift[CENTRAL], f_j);
1437 rvec_inc(fshift[t2], f_k);
1438 } /* 153 TOTAL */
1439 }
1440 return vtot;
1441 }
1442
1443 real dih_angle(const rvec xi, const rvec xj, const rvec xk, const rvec xl,
1444 const t_pbc *pbc,
1445 rvec r_ij, rvec r_kj, rvec r_kl, rvec m, rvec n,
1446 int *t1, int *t2, int *t3)
1447 {
1448 *t1 = pbc_rvec_sub(pbc, xi, xj, r_ij); /* 3 */
1449 *t2 = pbc_rvec_sub(pbc, xk, xj, r_kj); /* 3 */
1450 *t3 = pbc_rvec_sub(pbc, xk, xl, r_kl); /* 3 */
1451
1452 cprod(r_ij, r_kj, m); /* 9 */
1453 cprod(r_kj, r_kl, n); /* 9 */
1454 real phi = gmx_angle(m, n); /* 49 (assuming 25 for atan2) */
1455 real ipr = iprod(r_ij, n); /* 5 */
1456 real sign = (ipr < 0.0) ? -1.0 : 1.0;
1457 phi = sign*phi; /* 1 */
1458 /* 82 TOTAL */
1459 return phi;
1460 }
1461
1462
1463 #if GMX_SIMD_HAVE_REAL
1464
1465 /* As dih_angle above, but calculates 4 dihedral angles at once using SIMD,
1466 * also calculates the pre-factor required for the dihedral force update.
1467 * Note that bv and buf should be register aligned.
1468 */
1469 static inline void
1470 dih_angle_simd(const rvec *x,
1471 const int *ai, const int *aj, const int *ak, const int *al,
1472 const real *pbc_simd,
1473 SimdReal *phi_S,
1474 SimdReal *mx_S, SimdReal *my_S, SimdReal *mz_S,
1475 SimdReal *nx_S, SimdReal *ny_S, SimdReal *nz_S,
1476 SimdReal *nrkj_m2_S,
1477 SimdReal *nrkj_n2_S,
1478 SimdReal *p_S,
1479 SimdReal *q_S)
1480 {
1481 SimdReal xi_S, yi_S, zi_S;
1482 SimdReal xj_S, yj_S, zj_S;
1483 SimdReal xk_S, yk_S, zk_S;
1484 SimdReal xl_S, yl_S, zl_S;
1485 SimdReal rijx_S, rijy_S, rijz_S;
1486 SimdReal rkjx_S, rkjy_S, rkjz_S;
1487 SimdReal rklx_S, rkly_S, rklz_S;
1488 SimdReal cx_S, cy_S, cz_S;
1489 SimdReal cn_S;
1490 SimdReal s_S;
1491 SimdReal ipr_S;
1492 SimdReal iprm_S, iprn_S;
1493 SimdReal nrkj2_S, nrkj_1_S, nrkj_2_S, nrkj_S;
1494 SimdReal toler_S;
1495 SimdReal nrkj2_min_S;
1496 SimdReal real_eps_S;
1497
1498 /* Used to avoid division by zero.
1499 * We take into acount that we multiply the result by real_eps_S.
1500 */
1501 nrkj2_min_S = SimdReal(GMX_REAL_MIN/(2*GMX_REAL_EPS));
1502
1503 /* The value of the last significant bit (GMX_REAL_EPS is half of that) */
1504 real_eps_S = SimdReal(2*GMX_REAL_EPS);
1505
1506 /* Store the non PBC corrected distances packed and aligned */
1507 gatherLoadUTranspose<3>(reinterpret_cast<const real *>(x), ai, &xi_S, &yi_S, &zi_S);
1508 gatherLoadUTranspose<3>(reinterpret_cast<const real *>(x), aj, &xj_S, &yj_S, &zj_S);
1509 gatherLoadUTranspose<3>(reinterpret_cast<const real *>(x), ak, &xk_S, &yk_S, &zk_S);
1510 gatherLoadUTranspose<3>(reinterpret_cast<const real *>(x), al, &xl_S, &yl_S, &zl_S);
1511 rijx_S = xi_S - xj_S;
1512 rijy_S = yi_S - yj_S;
1513 rijz_S = zi_S - zj_S;
1514 rkjx_S = xk_S - xj_S;
1515 rkjy_S = yk_S - yj_S;
1516 rkjz_S = zk_S - zj_S;
1517 rklx_S = xk_S - xl_S;
1518 rkly_S = yk_S - yl_S;
1519 rklz_S = zk_S - zl_S;
1520
1521 pbc_correct_dx_simd(&rijx_S, &rijy_S, &rijz_S, pbc_simd);
1522 pbc_correct_dx_simd(&rkjx_S, &rkjy_S, &rkjz_S, pbc_simd);
1523 pbc_correct_dx_simd(&rklx_S, &rkly_S, &rklz_S, pbc_simd);
1524
1525 cprod(rijx_S, rijy_S, rijz_S,
1526 rkjx_S, rkjy_S, rkjz_S,
1527 mx_S, my_S, mz_S);
1528
1529 cprod(rkjx_S, rkjy_S, rkjz_S,
1530 rklx_S, rkly_S, rklz_S,
1531 nx_S, ny_S, nz_S);
1532
1533 cprod(*mx_S, *my_S, *mz_S,
1534 *nx_S, *ny_S, *nz_S,
1535 &cx_S, &cy_S, &cz_S);
1536
1537 cn_S = sqrt(norm2(cx_S, cy_S, cz_S));
1538
1539 s_S = iprod(*mx_S, *my_S, *mz_S, *nx_S, *ny_S, *nz_S);
1540
1541 /* Determine the dihedral angle, the sign might need correction */
1542 *phi_S = atan2(cn_S, s_S);
1543
1544 ipr_S = iprod(rijx_S, rijy_S, rijz_S,
1545 *nx_S, *ny_S, *nz_S);
1546
1547 iprm_S = norm2(*mx_S, *my_S, *mz_S);
1548 iprn_S = norm2(*nx_S, *ny_S, *nz_S);
1549
1550 nrkj2_S = norm2(rkjx_S, rkjy_S, rkjz_S);
1551
1552 /* Avoid division by zero. When zero, the result is multiplied by 0
1553 * anyhow, so the 3 max below do not affect the final result.
1554 */
1555 nrkj2_S = max(nrkj2_S, nrkj2_min_S);
1556 nrkj_1_S = invsqrt(nrkj2_S);
1557 nrkj_2_S = nrkj_1_S * nrkj_1_S;
1558 nrkj_S = nrkj2_S * nrkj_1_S;
1559
1560 toler_S = nrkj2_S * real_eps_S;
1561
1562 /* Here the plain-C code uses a conditional, but we can't do that in SIMD.
1563 * So we take a max with the tolerance instead. Since we multiply with
1564 * m or n later, the max does not affect the results.
1565 */
1566 iprm_S = max(iprm_S, toler_S);
1567 iprn_S = max(iprn_S, toler_S);
1568 *nrkj_m2_S = nrkj_S * inv(iprm_S);
1569 *nrkj_n2_S = nrkj_S * inv(iprn_S);
1570
1571 /* Set sign of phi_S with the sign of ipr_S; phi_S is currently positive */
1572 *phi_S = copysign(*phi_S, ipr_S);
1573 *p_S = iprod(rijx_S, rijy_S, rijz_S, rkjx_S, rkjy_S, rkjz_S);
1574 *p_S = *p_S * nrkj_2_S;
1575
1576 *q_S = iprod(rklx_S, rkly_S, rklz_S, rkjx_S, rkjy_S, rkjz_S);
1577 *q_S = *q_S * nrkj_2_S;
1578 }
1579
1580 #endif // GMX_SIMD_HAVE_REAL
1581
1582 void do_dih_fup(int i, int j, int k, int l, real ddphi,
1583 rvec r_ij, rvec r_kj, rvec r_kl,
1584 rvec m, rvec n, rvec4 f[], rvec fshift[],
1585 const t_pbc *pbc, const t_graph *g,
1586 const rvec x[], int t1, int t2, int t3)
1587 {
1588 /* 143 FLOPS */
1589 rvec f_i, f_j, f_k, f_l;
1590 rvec uvec, vvec, svec, dx_jl;
1591 real iprm, iprn, nrkj, nrkj2, nrkj_1, nrkj_2;
1592 real a, b, p, q, toler;
1593 ivec jt, dt_ij, dt_kj, dt_lj;
1594
1595 iprm = iprod(m, m); /* 5 */
1596 iprn = iprod(n, n); /* 5 */
1597 nrkj2 = iprod(r_kj, r_kj); /* 5 */
1598 toler = nrkj2*GMX_REAL_EPS;
1599 if ((iprm > toler) && (iprn > toler))
1600 {
1601 nrkj_1 = gmx::invsqrt(nrkj2); /* 10 */
1602 nrkj_2 = nrkj_1*nrkj_1; /* 1 */
1603 nrkj = nrkj2*nrkj_1; /* 1 */
1604 a = -ddphi*nrkj/iprm; /* 11 */
1605 svmul(a, m, f_i); /* 3 */
1606 b = ddphi*nrkj/iprn; /* 11 */
1607 svmul(b, n, f_l); /* 3 */
1608 p = iprod(r_ij, r_kj); /* 5 */
1609 p *= nrkj_2; /* 1 */
1610 q = iprod(r_kl, r_kj); /* 5 */
1611 q *= nrkj_2; /* 1 */
1612 svmul(p, f_i, uvec); /* 3 */
1613 svmul(q, f_l, vvec); /* 3 */
1614 rvec_sub(uvec, vvec, svec); /* 3 */
1615 rvec_sub(f_i, svec, f_j); /* 3 */
1616 rvec_add(f_l, svec, f_k); /* 3 */
1617 rvec_inc(f[i], f_i); /* 3 */
1618 rvec_dec(f[j], f_j); /* 3 */
1619 rvec_dec(f[k], f_k); /* 3 */
1620 rvec_inc(f[l], f_l); /* 3 */
1621
1622 if (g)
1623 {
1624 copy_ivec(SHIFT_IVEC(g, j), jt);
1625 ivec_sub(SHIFT_IVEC(g, i), jt, dt_ij);
1626 ivec_sub(SHIFT_IVEC(g, k), jt, dt_kj);
1627 ivec_sub(SHIFT_IVEC(g, l), jt, dt_lj);
1628 t1 = IVEC2IS(dt_ij);
1629 t2 = IVEC2IS(dt_kj);
1630 t3 = IVEC2IS(dt_lj);
1631 }
1632 else if (pbc)
1633 {
1634 t3 = pbc_rvec_sub(pbc, x[l], x[j], dx_jl);
1635 }
1636 else
1637 {
1638 t3 = CENTRAL;
1639 }
1640
1641 rvec_inc(fshift[t1], f_i);
1642 rvec_dec(fshift[CENTRAL], f_j);
1643 rvec_dec(fshift[t2], f_k);
1644 rvec_inc(fshift[t3], f_l);
1645 }
1646 /* 112 TOTAL */
1647 }
1648
1649 /* As do_dih_fup above, but without shift forces */
1650 static void
1651 do_dih_fup_noshiftf(int i, int j, int k, int l, real ddphi,
1652 rvec r_ij, rvec r_kj, rvec r_kl,
1653 rvec m, rvec n, rvec4 f[])
1654 {
1655 rvec f_i, f_j, f_k, f_l;
1656 rvec uvec, vvec, svec;
1657 real iprm, iprn, nrkj, nrkj2, nrkj_1, nrkj_2;
1658 real a, b, p, q, toler;
1659
1660 iprm = iprod(m, m); /* 5 */
1661 iprn = iprod(n, n); /* 5 */
1662 nrkj2 = iprod(r_kj, r_kj); /* 5 */
1663 toler = nrkj2*GMX_REAL_EPS;
1664 if ((iprm > toler) && (iprn > toler))
1665 {
1666 nrkj_1 = gmx::invsqrt(nrkj2); /* 10 */
1667 nrkj_2 = nrkj_1*nrkj_1; /* 1 */
1668 nrkj = nrkj2*nrkj_1; /* 1 */
1669 a = -ddphi*nrkj/iprm; /* 11 */
1670 svmul(a, m, f_i); /* 3 */
1671 b = ddphi*nrkj/iprn; /* 11 */
1672 svmul(b, n, f_l); /* 3 */
1673 p = iprod(r_ij, r_kj); /* 5 */
1674 p *= nrkj_2; /* 1 */
1675 q = iprod(r_kl, r_kj); /* 5 */
1676 q *= nrkj_2; /* 1 */
1677 svmul(p, f_i, uvec); /* 3 */
1678 svmul(q, f_l, vvec); /* 3 */
1679 rvec_sub(uvec, vvec, svec); /* 3 */
1680 rvec_sub(f_i, svec, f_j); /* 3 */
1681 rvec_add(f_l, svec, f_k); /* 3 */
1682 rvec_inc(f[i], f_i); /* 3 */
1683 rvec_dec(f[j], f_j); /* 3 */
1684 rvec_dec(f[k], f_k); /* 3 */
1685 rvec_inc(f[l], f_l); /* 3 */
1686 }
1687 }
1688
1689 #if GMX_SIMD_HAVE_REAL
1690 /* As do_dih_fup_noshiftf above, but with SIMD and pre-calculated pre-factors */
1691 static inline void gmx_simdcall
1692 do_dih_fup_noshiftf_simd(const int *ai, const int *aj, const int *ak, const int *al,
1693 SimdReal p, SimdReal q,
1694 SimdReal f_i_x, SimdReal f_i_y, SimdReal f_i_z,
1695 SimdReal mf_l_x, SimdReal mf_l_y, SimdReal mf_l_z,
1696 rvec4 f[])
1697 {
1698 SimdReal sx = p * f_i_x + q * mf_l_x;
1699 SimdReal sy = p * f_i_y + q * mf_l_y;
1700 SimdReal sz = p * f_i_z + q * mf_l_z;
1701 SimdReal f_j_x = f_i_x - sx;
1702 SimdReal f_j_y = f_i_y - sy;
1703 SimdReal f_j_z = f_i_z - sz;
1704 SimdReal f_k_x = mf_l_x - sx;
1705 SimdReal f_k_y = mf_l_y - sy;
1706 SimdReal f_k_z = mf_l_z - sz;
1707 transposeScatterIncrU<4>(reinterpret_cast<real *>(f), ai, f_i_x, f_i_y, f_i_z);
1708 transposeScatterDecrU<4>(reinterpret_cast<real *>(f), aj, f_j_x, f_j_y, f_j_z);
1709 transposeScatterIncrU<4>(reinterpret_cast<real *>(f), ak, f_k_x, f_k_y, f_k_z);
1710 transposeScatterDecrU<4>(reinterpret_cast<real *>(f), al, mf_l_x, mf_l_y, mf_l_z);
1711 }
1712 #endif // GMX_SIMD_HAVE_REAL
1713
1714 static real dopdihs(real cpA, real cpB, real phiA, real phiB, int mult,
1715 real phi, real lambda, real *V, real *F)
1716 {
1717 real v, dvdlambda, mdphi, v1, sdphi, ddphi;
1718 real L1 = 1.0 - lambda;
1719 real ph0 = (L1*phiA + lambda*phiB)*DEG2RAD;
1720 real dph0 = (phiB - phiA)*DEG2RAD;
1721 real cp = L1*cpA + lambda*cpB;
1722
1723 mdphi = mult*phi - ph0;
1724 sdphi = std::sin(mdphi);
1725 ddphi = -cp*mult*sdphi;
1726 v1 = 1.0 + std::cos(mdphi);
1727 v = cp*v1;
1728
1729 dvdlambda = (cpB - cpA)*v1 + cp*dph0*sdphi;
1730
1731 *V = v;
1732 *F = ddphi;
1733
1734 return dvdlambda;
1735
1736 /* That was 40 flops */
1737 }
1738
1739 static void
1740 dopdihs_noener(real cpA, real cpB, real phiA, real phiB, int mult,
1741 real phi, real lambda, real *F)
1742 {
1743 real mdphi, sdphi, ddphi;
1744 real L1 = 1.0 - lambda;
1745 real ph0 = (L1*phiA + lambda*phiB)*DEG2RAD;
1746 real cp = L1*cpA + lambda*cpB;
1747
1748 mdphi = mult*phi - ph0;
1749 sdphi = std::sin(mdphi);
1750 ddphi = -cp*mult*sdphi;
1751
1752 *F = ddphi;
1753
1754 /* That was 20 flops */
1755 }
1756
1757 static real dopdihs_min(real cpA, real cpB, real phiA, real phiB, int mult,
1758 real phi, real lambda, real *V, real *F)
1759 /* similar to dopdihs, except for a minus sign *
1760 * and a different treatment of mult/phi0 */
1761 {
1762 real v, dvdlambda, mdphi, v1, sdphi, ddphi;
1763 real L1 = 1.0 - lambda;
1764 real ph0 = (L1*phiA + lambda*phiB)*DEG2RAD;
1765 real dph0 = (phiB - phiA)*DEG2RAD;
1766 real cp = L1*cpA + lambda*cpB;
1767
1768 mdphi = mult*(phi-ph0);
1769 sdphi = std::sin(mdphi);
1770 ddphi = cp*mult*sdphi;
1771 v1 = 1.0-std::cos(mdphi);
1772 v = cp*v1;
1773
1774 dvdlambda = (cpB-cpA)*v1 + cp*dph0*sdphi;
1775
1776 *V = v;
1777 *F = ddphi;
1778
1779 return dvdlambda;
1780
1781 /* That was 40 flops */
1782 }
1783
1784 real pdihs(int nbonds,
1785 const t_iatom forceatoms[], const t_iparams forceparams[],
1786 const rvec x[], rvec4 f[], rvec fshift[],
1787 const t_pbc *pbc, const t_graph *g,
1788 real lambda, real *dvdlambda,
1789 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
1790 int gmx_unused *global_atom_index)
1791 {
1792 int i, type, ai, aj, ak, al;
1793 int t1, t2, t3;
1794 rvec r_ij, r_kj, r_kl, m, n;
1795 real phi, ddphi, vpd, vtot;
1796
1797 vtot = 0.0;
1798
1799 for (i = 0; (i < nbonds); )
1800 {
1801 type = forceatoms[i++];
1802 ai = forceatoms[i++];
1803 aj = forceatoms[i++];
1804 ak = forceatoms[i++];
1805 al = forceatoms[i++];
1806
1807 phi = dih_angle(x[ai], x[aj], x[ak], x[al], pbc, r_ij, r_kj, r_kl, m, n,
1808 &t1, &t2, &t3); /* 84 */
1809 *dvdlambda += dopdihs(forceparams[type].pdihs.cpA,
1810 forceparams[type].pdihs.cpB,
1811 forceparams[type].pdihs.phiA,
1812 forceparams[type].pdihs.phiB,
1813 forceparams[type].pdihs.mult,
1814 phi, lambda, &vpd, &ddphi);
1815
1816 vtot += vpd;
1817 do_dih_fup(ai, aj, ak, al, ddphi, r_ij, r_kj, r_kl, m, n,
1818 f, fshift, pbc, g, x, t1, t2, t3); /* 112 */
1819
1820 } /* 223 TOTAL */
1821
1822 return vtot;
1823 }
1824
1825 void make_dp_periodic(real *dp) /* 1 flop? */
1826 {
1827 /* dp cannot be outside (-pi,pi) */
1828 if (*dp >= M_PI)
1829 {
1830 *dp -= 2*M_PI;
1831 }
1832 else if (*dp < -M_PI)
1833 {
1834 *dp += 2*M_PI;
1835 }
1836 }
1837
1838 /* As pdihs above, but without calculating energies and shift forces */
1839 void
1840 pdihs_noener(int nbonds,
1841 const t_iatom forceatoms[], const t_iparams forceparams[],
1842 const rvec x[], rvec4 f[],
1843 const t_pbc gmx_unused *pbc, const t_graph gmx_unused *g,
1844 real lambda,
1845 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
1846 int gmx_unused *global_atom_index)
1847 {
1848 int i, type, ai, aj, ak, al;
1849 int t1, t2, t3;
1850 rvec r_ij, r_kj, r_kl, m, n;
1851 real phi, ddphi_tot, ddphi;
1852
1853 for (i = 0; (i < nbonds); )
1854 {
1855 ai = forceatoms[i+1];
1856 aj = forceatoms[i+2];
1857 ak = forceatoms[i+3];
1858 al = forceatoms[i+4];
1859
1860 phi = dih_angle(x[ai], x[aj], x[ak], x[al], pbc, r_ij, r_kj, r_kl, m, n,
1861 &t1, &t2, &t3);
1862
1863 ddphi_tot = 0;
1864
1865 /* Loop over dihedrals working on the same atoms,
1866 * so we avoid recalculating angles and force distributions.
1867 */
1868 do
1869 {
1870 type = forceatoms[i];
1871 dopdihs_noener(forceparams[type].pdihs.cpA,
1872 forceparams[type].pdihs.cpB,
1873 forceparams[type].pdihs.phiA,
1874 forceparams[type].pdihs.phiB,
1875 forceparams[type].pdihs.mult,
1876 phi, lambda, &ddphi);
1877 ddphi_tot += ddphi;
1878
1879 i += 5;
1880 }
1881 while (i < nbonds &&
1882 forceatoms[i+1] == ai &&
1883 forceatoms[i+2] == aj &&
1884 forceatoms[i+3] == ak &&
1885 forceatoms[i+4] == al);
1886
1887 do_dih_fup_noshiftf(ai, aj, ak, al, ddphi_tot, r_ij, r_kj, r_kl, m, n, f);
1888 }
1889 }
1890
1891
1892 #if GMX_SIMD_HAVE_REAL
1893
1894 /* As pdihs_noner above, but using SIMD to calculate many dihedrals at once */
1895 void
1896 pdihs_noener_simd(int nbonds,
1897 const t_iatom forceatoms[], const t_iparams forceparams[],
1898 const rvec x[], rvec4 f[],
1899 const t_pbc *pbc, const t_graph gmx_unused *g,
1900 real gmx_unused lambda,
1901 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
1902 int gmx_unused *global_atom_index)
1903 {
1904 const int nfa1 = 5;
1905 int i, iu, s;
1906 int type;
1907 alignas(GMX_SIMD_ALIGNMENT) std::int32_t ai[GMX_SIMD_REAL_WIDTH];
1908 alignas(GMX_SIMD_ALIGNMENT) std::int32_t aj[GMX_SIMD_REAL_WIDTH];
1909 alignas(GMX_SIMD_ALIGNMENT) std::int32_t ak[GMX_SIMD_REAL_WIDTH];
1910 alignas(GMX_SIMD_ALIGNMENT) std::int32_t al[GMX_SIMD_REAL_WIDTH];
1911 alignas(GMX_SIMD_ALIGNMENT) real buf[3*GMX_SIMD_REAL_WIDTH];
1912 real *cp, *phi0, *mult;
1913 SimdReal deg2rad_S(DEG2RAD);
1914 SimdReal p_S, q_S;
1915 SimdReal phi0_S, phi_S;
1916 SimdReal mx_S, my_S, mz_S;
1917 SimdReal nx_S, ny_S, nz_S;
1918 SimdReal nrkj_m2_S, nrkj_n2_S;
1919 SimdReal cp_S, mdphi_S, mult_S;
1920 SimdReal sin_S, cos_S;
1921 SimdReal mddphi_S;
1922 SimdReal sf_i_S, msf_l_S;
1923 alignas(GMX_SIMD_ALIGNMENT) real pbc_simd[9*GMX_SIMD_REAL_WIDTH];
1924
1925 /* Extract aligned pointer for parameters and variables */
1926 cp = buf + 0*GMX_SIMD_REAL_WIDTH;
1927 phi0 = buf + 1*GMX_SIMD_REAL_WIDTH;
1928 mult = buf + 2*GMX_SIMD_REAL_WIDTH;
1929
1930 set_pbc_simd(pbc, pbc_simd);
1931
1932 /* nbonds is the number of dihedrals times nfa1, here we step GMX_SIMD_REAL_WIDTH dihs */
1933 for (i = 0; (i < nbonds); i += GMX_SIMD_REAL_WIDTH*nfa1)
1934 {
1935 /* Collect atoms quadruplets for GMX_SIMD_REAL_WIDTH dihedrals.
1936 * iu indexes into forceatoms, we should not let iu go beyond nbonds.
1937 */
1938 iu = i;
1939 for (s = 0; s < GMX_SIMD_REAL_WIDTH; s++)
1940 {
1941 type = forceatoms[iu];
1942 ai[s] = forceatoms[iu+1];
1943 aj[s] = forceatoms[iu+2];
1944 ak[s] = forceatoms[iu+3];
1945 al[s] = forceatoms[iu+4];
1946
1947 /* At the end fill the arrays with the last atoms and 0 params */
1948 if (i + s*nfa1 < nbonds)
1949 {
1950 cp[s] = forceparams[type].pdihs.cpA;
1951 phi0[s] = forceparams[type].pdihs.phiA;
1952 mult[s] = forceparams[type].pdihs.mult;
1953
1954 if (iu + nfa1 < nbonds)
1955 {
1956 iu += nfa1;
1957 }
1958 }
1959 else
1960 {
1961 cp[s] = 0;
1962 phi0[s] = 0;
1963 mult[s] = 0;
1964 }
1965 }
1966
1967 /* Caclulate GMX_SIMD_REAL_WIDTH dihedral angles at once */
1968 dih_angle_simd(x, ai, aj, ak, al, pbc_simd,
1969 &phi_S,
1970 &mx_S, &my_S, &mz_S,
1971 &nx_S, &ny_S, &nz_S,
1972 &nrkj_m2_S,
1973 &nrkj_n2_S,
1974 &p_S, &q_S);
1975
1976 cp_S = load<SimdReal>(cp);
1977 phi0_S = load<SimdReal>(phi0) * deg2rad_S;
1978 mult_S = load<SimdReal>(mult);
1979
1980 mdphi_S = fms(mult_S, phi_S, phi0_S);
1981
1982 /* Calculate GMX_SIMD_REAL_WIDTH sines at once */
1983 sincos(mdphi_S, &sin_S, &cos_S);
1984 mddphi_S = cp_S * mult_S * sin_S;
1985 sf_i_S = mddphi_S * nrkj_m2_S;
1986 msf_l_S = mddphi_S * nrkj_n2_S;
1987
1988 /* After this m?_S will contain f[i] */
1989 mx_S = sf_i_S * mx_S;
1990 my_S = sf_i_S * my_S;
1991 mz_S = sf_i_S * mz_S;
1992
1993 /* After this m?_S will contain -f[l] */
1994 nx_S = msf_l_S * nx_S;
1995 ny_S = msf_l_S * ny_S;
1996 nz_S = msf_l_S * nz_S;
1997
1998 do_dih_fup_noshiftf_simd(ai, aj, ak, al,
1999 p_S, q_S,
2000 mx_S, my_S, mz_S,
2001 nx_S, ny_S, nz_S,
2002 f);
2003 }
2004 }
2005
2006 /* This is mostly a copy of pdihs_noener_simd above, but with using
2007 * the RB potential instead of a harmonic potential.
2008 * This function can replace rbdihs() when no energy and virial are needed.
2009 */
2010 void
2011 rbdihs_noener_simd(int nbonds,
2012 const t_iatom forceatoms[], const t_iparams forceparams[],
2013 const rvec x[], rvec4 f[],
2014 const t_pbc *pbc, const t_graph gmx_unused *g,
2015 real gmx_unused lambda,
2016 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
2017 int gmx_unused *global_atom_index)
2018 {
2019 const int nfa1 = 5;
2020 int i, iu, s, j;
2021 int type;
2022 alignas(GMX_SIMD_ALIGNMENT) std::int32_t ai[GMX_SIMD_REAL_WIDTH];
2023 alignas(GMX_SIMD_ALIGNMENT) std::int32_t aj[GMX_SIMD_REAL_WIDTH];
2024 alignas(GMX_SIMD_ALIGNMENT) std::int32_t ak[GMX_SIMD_REAL_WIDTH];
2025 alignas(GMX_SIMD_ALIGNMENT) std::int32_t al[GMX_SIMD_REAL_WIDTH];
2026 alignas(GMX_SIMD_ALIGNMENT) real parm[NR_RBDIHS*GMX_SIMD_REAL_WIDTH];
2027
2028 SimdReal p_S, q_S;
2029 SimdReal phi_S;
2030 SimdReal ddphi_S, cosfac_S;
2031 SimdReal mx_S, my_S, mz_S;
2032 SimdReal nx_S, ny_S, nz_S;
2033 SimdReal nrkj_m2_S, nrkj_n2_S;
2034 SimdReal parm_S, c_S;
2035 SimdReal sin_S, cos_S;
2036 SimdReal sf_i_S, msf_l_S;
2037 alignas(GMX_SIMD_ALIGNMENT) real pbc_simd[9*GMX_SIMD_REAL_WIDTH];
2038
2039 SimdReal pi_S(M_PI);
2040 SimdReal one_S(1.0);
2041
2042 set_pbc_simd(pbc, pbc_simd);
2043
2044 /* nbonds is the number of dihedrals times nfa1, here we step GMX_SIMD_REAL_WIDTH dihs */
2045 for (i = 0; (i < nbonds); i += GMX_SIMD_REAL_WIDTH*nfa1)
2046 {
2047 /* Collect atoms quadruplets for GMX_SIMD_REAL_WIDTH dihedrals.
2048 * iu indexes into forceatoms, we should not let iu go beyond nbonds.
2049 */
2050 iu = i;
2051 for (s = 0; s < GMX_SIMD_REAL_WIDTH; s++)
2052 {
2053 type = forceatoms[iu];
2054 ai[s] = forceatoms[iu+1];
2055 aj[s] = forceatoms[iu+2];
2056 ak[s] = forceatoms[iu+3];
2057 al[s] = forceatoms[iu+4];
2058
2059 /* At the end fill the arrays with the last atoms and 0 params */
2060 if (i + s*nfa1 < nbonds)
2061 {
2062 /* We don't need the first parameter, since that's a constant
2063 * which only affects the energies, not the forces.
2064 */
2065 for (j = 1; j < NR_RBDIHS; j++)
2066 {
2067 parm[j*GMX_SIMD_REAL_WIDTH + s] =
2068 forceparams[type].rbdihs.rbcA[j];
2069 }
2070
2071 if (iu + nfa1 < nbonds)
2072 {
2073 iu += nfa1;
2074 }
2075 }
2076 else
2077 {
2078 for (j = 1; j < NR_RBDIHS; j++)
2079 {
2080 parm[j*GMX_SIMD_REAL_WIDTH + s] = 0;
2081 }
2082 }
2083 }
2084
2085 /* Caclulate GMX_SIMD_REAL_WIDTH dihedral angles at once */
2086 dih_angle_simd(x, ai, aj, ak, al, pbc_simd,
2087 &phi_S,
2088 &mx_S, &my_S, &mz_S,
2089 &nx_S, &ny_S, &nz_S,
2090 &nrkj_m2_S,
2091 &nrkj_n2_S,
2092 &p_S, &q_S);
2093
2094 /* Change to polymer convention */
2095 phi_S = phi_S - pi_S;
2096
2097 sincos(phi_S, &sin_S, &cos_S);
2098
2099 ddphi_S = setZero();
2100 c_S = one_S;
2101 cosfac_S = one_S;
2102 for (j = 1; j < NR_RBDIHS; j++)
2103 {
2104 parm_S = load<SimdReal>(parm + j*GMX_SIMD_REAL_WIDTH);
2105 ddphi_S = fma(c_S * parm_S, cosfac_S, ddphi_S);
2106 cosfac_S = cosfac_S * cos_S;
2107 c_S = c_S + one_S;
2108 }
2109
2110 /* Note that here we do not use the minus sign which is present
2111 * in the normal RB code. This is corrected for through (m)sf below.
2112 */
2113 ddphi_S = ddphi_S * sin_S;
2114
2115 sf_i_S = ddphi_S * nrkj_m2_S;
2116 msf_l_S = ddphi_S * nrkj_n2_S;
2117
2118 /* After this m?_S will contain f[i] */
2119 mx_S = sf_i_S * mx_S;
2120 my_S = sf_i_S * my_S;
2121 mz_S = sf_i_S * mz_S;
2122
2123 /* After this m?_S will contain -f[l] */
2124 nx_S = msf_l_S * nx_S;
2125 ny_S = msf_l_S * ny_S;
2126 nz_S = msf_l_S * nz_S;
2127
2128 do_dih_fup_noshiftf_simd(ai, aj, ak, al,
2129 p_S, q_S,
2130 mx_S, my_S, mz_S,
2131 nx_S, ny_S, nz_S,
2132 f);
2133 }
2134 }
2135
2136 #endif // GMX_SIMD_HAVE_REAL
2137
2138
2139 real idihs(int nbonds,
2140 const t_iatom forceatoms[], const t_iparams forceparams[],
2141 const rvec x[], rvec4 f[], rvec fshift[],
2142 const t_pbc *pbc, const t_graph *g,
2143 real lambda, real *dvdlambda,
2144 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
2145 int gmx_unused *global_atom_index)
2146 {
2147 int i, type, ai, aj, ak, al;
2148 int t1, t2, t3;
2149 real phi, phi0, dphi0, ddphi, vtot;
2150 rvec r_ij, r_kj, r_kl, m, n;
2151 real L1, kk, dp, dp2, kA, kB, pA, pB, dvdl_term;
2152
2153 L1 = 1.0-lambda;
2154 dvdl_term = 0;
2155 vtot = 0.0;
2156 for (i = 0; (i < nbonds); )
2157 {
2158 type = forceatoms[i++];
2159 ai = forceatoms[i++];
2160 aj = forceatoms[i++];
2161 ak = forceatoms[i++];
2162 al = forceatoms[i++];
2163
2164 phi = dih_angle(x[ai], x[aj], x[ak], x[al], pbc, r_ij, r_kj, r_kl, m, n,
2165 &t1, &t2, &t3); /* 84 */
2166
2167 /* phi can jump if phi0 is close to Pi/-Pi, which will cause huge
2168 * force changes if we just apply a normal harmonic.
2169 * Instead, we first calculate phi-phi0 and take it modulo (-Pi,Pi).
2170 * This means we will never have the periodicity problem, unless
2171 * the dihedral is Pi away from phiO, which is very unlikely due to
2172 * the potential.
2173 */
2174 kA = forceparams[type].harmonic.krA;
2175 kB = forceparams[type].harmonic.krB;
2176 pA = forceparams[type].harmonic.rA;
2177 pB = forceparams[type].harmonic.rB;
2178
2179 kk = L1*kA + lambda*kB;
2180 phi0 = (L1*pA + lambda*pB)*DEG2RAD;
2181 dphi0 = (pB - pA)*DEG2RAD;
2182
2183 dp = phi-phi0;
2184
2185 make_dp_periodic(&dp);
2186
2187 dp2 = dp*dp;
2188
2189 vtot += 0.5*kk*dp2;
2190 ddphi = -kk*dp;
2191
2192 dvdl_term += 0.5*(kB - kA)*dp2 - kk*dphi0*dp;
2193
2194 do_dih_fup(ai, aj, ak, al, -ddphi, r_ij, r_kj, r_kl, m, n,
2195 f, fshift, pbc, g, x, t1, t2, t3); /* 112 */
2196 /* 218 TOTAL */
2197 }
2198
2199 *dvdlambda += dvdl_term;
2200 return vtot;
2201 }
2202
2203 static real low_angres(int nbonds,
2204 const t_iatom forceatoms[], const t_iparams forceparams[],
2205 const rvec x[], rvec4 f[], rvec fshift[],
2206 const t_pbc *pbc, const t_graph *g,
2207 real lambda, real *dvdlambda,
2208 gmx_bool bZAxis)
2209 {
2210 int i, m, type, ai, aj, ak, al;
2211 int t1, t2;
2212 real phi, cos_phi, cos_phi2, vid, vtot, dVdphi;
2213 rvec r_ij, r_kl, f_i, f_k = {0, 0, 0};
2214 real st, sth, nrij2, nrkl2, c, cij, ckl;
2215
2216 ivec dt;
2217 t2 = 0; /* avoid warning with gcc-3.3. It is never used uninitialized */
2218
2219 vtot = 0.0;
2220 ak = al = 0; /* to avoid warnings */
2221 for (i = 0; i < nbonds; )
2222 {
2223 type = forceatoms[i++];
2224 ai = forceatoms[i++];
2225 aj = forceatoms[i++];
2226 t1 = pbc_rvec_sub(pbc, x[aj], x[ai], r_ij); /* 3 */
2227 if (!bZAxis)
2228 {
2229 ak = forceatoms[i++];
2230 al = forceatoms[i++];
2231 t2 = pbc_rvec_sub(pbc, x[al], x[ak], r_kl); /* 3 */
2232 }
2233 else
2234 {
2235 r_kl[XX] = 0;
2236 r_kl[YY] = 0;
2237 r_kl[ZZ] = 1;
2238 }
2239
2240 cos_phi = cos_angle(r_ij, r_kl); /* 25 */
2241 phi = std::acos(cos_phi); /* 10 */
2242
2243 *dvdlambda += dopdihs_min(forceparams[type].pdihs.cpA,
2244 forceparams[type].pdihs.cpB,
2245 forceparams[type].pdihs.phiA,
2246 forceparams[type].pdihs.phiB,
2247 forceparams[type].pdihs.mult,
2248 phi, lambda, &vid, &dVdphi); /* 40 */
2249
2250 vtot += vid;
2251
2252 cos_phi2 = gmx::square(cos_phi); /* 1 */
2253 if (cos_phi2 < 1)
2254 {
2255 st = -dVdphi*gmx::invsqrt(1 - cos_phi2); /* 12 */
2256 sth = st*cos_phi; /* 1 */
2257 nrij2 = iprod(r_ij, r_ij); /* 5 */
2258 nrkl2 = iprod(r_kl, r_kl); /* 5 */
2259
2260 c = st*gmx::invsqrt(nrij2*nrkl2); /* 11 */
2261 cij = sth/nrij2; /* 10 */
2262 ckl = sth/nrkl2; /* 10 */
2263
2264 for (m = 0; m < DIM; m++) /* 18+18 */
2265 {
2266 f_i[m] = (c*r_kl[m]-cij*r_ij[m]);
2267 f[ai][m] += f_i[m];
2268 f[aj][m] -= f_i[m];
2269 if (!bZAxis)
2270 {
2271 f_k[m] = (c*r_ij[m]-ckl*r_kl[m]);
2272 f[ak][m] += f_k[m];
2273 f[al][m] -= f_k[m];
2274 }
2275 }
2276
2277 if (g)
2278 {
2279 ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, aj), dt);
2280 t1 = IVEC2IS(dt);
2281 }
2282 rvec_inc(fshift[t1], f_i);
2283 rvec_dec(fshift[CENTRAL], f_i);
2284 if (!bZAxis)
2285 {
2286 if (g)
2287 {
2288 ivec_sub(SHIFT_IVEC(g, ak), SHIFT_IVEC(g, al), dt);
2289 t2 = IVEC2IS(dt);
2290 }
2291 rvec_inc(fshift[t2], f_k);
2292 rvec_dec(fshift[CENTRAL], f_k);
2293 }
2294 }
2295 }
2296
2297 return vtot; /* 184 / 157 (bZAxis) total */
2298 }
2299
2300 real angres(int nbonds,
2301 const t_iatom forceatoms[], const t_iparams forceparams[],
2302 const rvec x[], rvec4 f[], rvec fshift[],
2303 const t_pbc *pbc, const t_graph *g,
2304 real lambda, real *dvdlambda,
2305 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
2306 int gmx_unused *global_atom_index)
2307 {
2308 return low_angres(nbonds, forceatoms, forceparams, x, f, fshift, pbc, g,
2309 lambda, dvdlambda, FALSE);
2310 }
2311
2312 real angresz(int nbonds,
2313 const t_iatom forceatoms[], const t_iparams forceparams[],
2314 const rvec x[], rvec4 f[], rvec fshift[],
2315 const t_pbc *pbc, const t_graph *g,
2316 real lambda, real *dvdlambda,
2317 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
2318 int gmx_unused *global_atom_index)
2319 {
2320 return low_angres(nbonds, forceatoms, forceparams, x, f, fshift, pbc, g,
2321 lambda, dvdlambda, TRUE);
2322 }
2323
2324 real dihres(int nbonds,
2325 const t_iatom forceatoms[], const t_iparams forceparams[],
2326 const rvec x[], rvec4 f[], rvec fshift[],
2327 const t_pbc *pbc, const t_graph *g,
2328 real lambda, real *dvdlambda,
2329 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
2330 int gmx_unused *global_atom_index)
2331 {
2332 real vtot = 0;
2333 int ai, aj, ak, al, i, type, t1, t2, t3;
2334 real phi0A, phi0B, dphiA, dphiB, kfacA, kfacB, phi0, dphi, kfac;
2335 real phi, ddphi, ddp, ddp2, dp, d2r, L1;
2336 rvec r_ij, r_kj, r_kl, m, n;
2337
2338 L1 = 1.0-lambda;
2339
2340 d2r = DEG2RAD;
2341
2342 for (i = 0; (i < nbonds); )
2343 {
2344 type = forceatoms[i++];
2345 ai = forceatoms[i++];
2346 aj = forceatoms[i++];
2347 ak = forceatoms[i++];
2348 al = forceatoms[i++];
2349
2350 phi0A = forceparams[type].dihres.phiA*d2r;
2351 dphiA = forceparams[type].dihres.dphiA*d2r;
2352 kfacA = forceparams[type].dihres.kfacA;
2353
2354 phi0B = forceparams[type].dihres.phiB*d2r;
2355 dphiB = forceparams[type].dihres.dphiB*d2r;
2356 kfacB = forceparams[type].dihres.kfacB;
2357
2358 phi0 = L1*phi0A + lambda*phi0B;
2359 dphi = L1*dphiA + lambda*dphiB;
2360 kfac = L1*kfacA + lambda*kfacB;
2361
2362 phi = dih_angle(x[ai], x[aj], x[ak], x[al], pbc, r_ij, r_kj, r_kl, m, n,
2363 &t1, &t2, &t3);
2364 /* 84 flops */
2365
2366 /* phi can jump if phi0 is close to Pi/-Pi, which will cause huge
2367 * force changes if we just apply a normal harmonic.
2368 * Instead, we first calculate phi-phi0 and take it modulo (-Pi,Pi).
2369 * This means we will never have the periodicity problem, unless
2370 * the dihedral is Pi away from phiO, which is very unlikely due to
2371 * the potential.
2372 */
2373 dp = phi-phi0;
2374 make_dp_periodic(&dp);
2375
2376 if (dp > dphi)
2377 {
2378 ddp = dp-dphi;
2379 }
2380 else if (dp < -dphi)
2381 {
2382 ddp = dp+dphi;
2383 }
2384 else
2385 {
2386 ddp = 0;
2387 }
2388
2389 if (ddp != 0.0)
2390 {
2391 ddp2 = ddp*ddp;
2392 vtot += 0.5*kfac*ddp2;
2393 ddphi = kfac*ddp;
2394
2395 *dvdlambda += 0.5*(kfacB - kfacA)*ddp2;
2396 /* lambda dependence from changing restraint distances */
2397 if (ddp > 0)
2398 {
2399 *dvdlambda -= kfac*ddp*((dphiB - dphiA)+(phi0B - phi0A));
2400 }
2401 else if (ddp < 0)
2402 {
2403 *dvdlambda += kfac*ddp*((dphiB - dphiA)-(phi0B - phi0A));
2404 }
2405 do_dih_fup(ai, aj, ak, al, ddphi, r_ij, r_kj, r_kl, m, n,
2406 f, fshift, pbc, g, x, t1, t2, t3); /* 112 */
2407 }
2408 }
2409 return vtot;
2410 }
2411
2412
2413 real unimplemented(int gmx_unused nbonds,
2414 const t_iatom gmx_unused forceatoms[], const t_iparams gmx_unused forceparams[],
2415 const rvec gmx_unused x[], rvec4 gmx_unused f[], rvec gmx_unused fshift[],
2416 const t_pbc gmx_unused *pbc, const t_graph gmx_unused *g,
2417 real gmx_unused lambda, real gmx_unused *dvdlambda,
2418 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
2419 int gmx_unused *global_atom_index)
2420 {
2421 gmx_impl("*** you are using a not implemented function");
2422 }
2423
2424 real restrangles(int nbonds,
2425 const t_iatom forceatoms[], const t_iparams forceparams[],
2426 const rvec x[], rvec4 f[], rvec fshift[],
2427 const t_pbc *pbc, const t_graph *g,
2428 real gmx_unused lambda, real gmx_unused *dvdlambda,
2429 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
2430 int gmx_unused *global_atom_index)
2431 {
2432 int i, d, ai, aj, ak, type, m;
2433 int t1, t2;
2434 real v, vtot;
2435 ivec jt, dt_ij, dt_kj;
2436 rvec f_i, f_j, f_k;
2437 double prefactor, ratio_ante, ratio_post;
2438 rvec delta_ante, delta_post, vec_temp;
2439
2440 vtot = 0.0;
2441 for (i = 0; (i < nbonds); )
2442 {
2443 type = forceatoms[i++];
2444 ai = forceatoms[i++];
2445 aj = forceatoms[i++];
2446 ak = forceatoms[i++];
2447
2448 t1 = pbc_rvec_sub(pbc, x[ai], x[aj], vec_temp);
2449 pbc_rvec_sub(pbc, x[aj], x[ai], delta_ante);
2450 t2 = pbc_rvec_sub(pbc, x[ak], x[aj], delta_post);
2451
2452
2453 /* This function computes factors needed for restricted angle potential.
2454 * The restricted angle potential is used in coarse-grained simulations to avoid singularities
2455 * when three particles align and the dihedral angle and dihedral potential
2456 * cannot be calculated. This potential is calculated using the formula:
2457 real restrangles(int nbonds,
2458 const t_iatom forceatoms[],const t_iparams forceparams[],
2459 const rvec x[],rvec4 f[],rvec fshift[],
2460 const t_pbc *pbc,const t_graph *g,
2461 real gmx_unused lambda,real gmx_unused *dvdlambda,
2462 const t_mdatoms gmx_unused *md,t_fcdata gmx_unused *fcd,
2463 int gmx_unused *global_atom_index)
2464 {
2465 int i, d, ai, aj, ak, type, m;
2466 int t1, t2;
2467 rvec r_ij,r_kj;
2468 real v, vtot;
2469 ivec jt,dt_ij,dt_kj;
2470 rvec f_i, f_j, f_k;
2471 real prefactor, ratio_ante, ratio_post;
2472 rvec delta_ante, delta_post, vec_temp;
2473
2474 vtot = 0.0;
2475 for(i=0; (i<nbonds); )
2476 {
2477 type = forceatoms[i++];
2478 ai = forceatoms[i++];
2479 aj = forceatoms[i++];
2480 ak = forceatoms[i++];
2481
2482 * \f[V_{\rm ReB}(\theta_i) = \frac{1}{2} k_{\theta} \frac{(\cos\theta_i - \cos\theta_0)^2}
2483 * {\sin^2\theta_i}\f] ({eq:ReB} and ref \cite{MonicaGoga2013} from the manual).
2484 * For more explanations see comments file "restcbt.h". */
2485
2486 compute_factors_restangles(type, forceparams, delta_ante, delta_post,
2487 &prefactor, &ratio_ante, &ratio_post, &v);
2488
2489 /* Forces are computed per component */
2490 for (d = 0; d < DIM; d++)
2491 {
2492 f_i[d] = prefactor * (ratio_ante * delta_ante[d] - delta_post[d]);
2493 f_j[d] = prefactor * ((ratio_post + 1.0) * delta_post[d] - (ratio_ante + 1.0) * delta_ante[d]);
2494 f_k[d] = prefactor * (delta_ante[d] - ratio_post * delta_post[d]);
2495 }
2496
2497 /* Computation of potential energy */
2498
2499 vtot += v;
2500
2501 /* Update forces */
2502
2503 for (m = 0; (m < DIM); m++)
2504 {
2505 f[ai][m] += f_i[m];
2506 f[aj][m] += f_j[m];
2507 f[ak][m] += f_k[m];
2508 }
2509
2510 if (g)
2511 {
2512 copy_ivec(SHIFT_IVEC(g, aj), jt);
2513 ivec_sub(SHIFT_IVEC(g, ai), jt, dt_ij);
2514 ivec_sub(SHIFT_IVEC(g, ak), jt, dt_kj);
2515 t1 = IVEC2IS(dt_ij);
2516 t2 = IVEC2IS(dt_kj);
2517 }
2518
2519 rvec_inc(fshift[t1], f_i);
2520 rvec_inc(fshift[CENTRAL], f_j);
2521 rvec_inc(fshift[t2], f_k);
2522 }
2523 return vtot;
2524 }
2525
2526
2527 real restrdihs(int nbonds,
2528 const t_iatom forceatoms[], const t_iparams forceparams[],
2529 const rvec x[], rvec4 f[], rvec fshift[],
2530 const t_pbc *pbc, const t_graph *g,
2531 real gmx_unused lambda, real gmx_unused *dvlambda,
2532 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
2533 int gmx_unused *global_atom_index)
2534 {
2535 int i, d, type, ai, aj, ak, al;
2536 rvec f_i, f_j, f_k, f_l;
2537 rvec dx_jl;
2538 ivec jt, dt_ij, dt_kj, dt_lj;
2539 int t1, t2, t3;
2540 real v, vtot;
2541 rvec delta_ante, delta_crnt, delta_post, vec_temp;
2542 real factor_phi_ai_ante, factor_phi_ai_crnt, factor_phi_ai_post;
2543 real factor_phi_aj_ante, factor_phi_aj_crnt, factor_phi_aj_post;
2544 real factor_phi_ak_ante, factor_phi_ak_crnt, factor_phi_ak_post;
2545 real factor_phi_al_ante, factor_phi_al_crnt, factor_phi_al_post;
2546 real prefactor_phi;
2547
2548
2549 vtot = 0.0;
2550 for (i = 0; (i < nbonds); )
2551 {
2552 type = forceatoms[i++];
2553 ai = forceatoms[i++];
2554 aj = forceatoms[i++];
2555 ak = forceatoms[i++];
2556 al = forceatoms[i++];
2557
2558 t1 = pbc_rvec_sub(pbc, x[ai], x[aj], vec_temp);
2559 pbc_rvec_sub(pbc, x[aj], x[ai], delta_ante);
2560 t2 = pbc_rvec_sub(pbc, x[ak], x[aj], delta_crnt);
2561 pbc_rvec_sub(pbc, x[ak], x[al], vec_temp);
2562 pbc_rvec_sub(pbc, x[al], x[ak], delta_post);
2563
2564 /* This function computes factors needed for restricted angle potential.
2565 * The restricted angle potential is used in coarse-grained simulations to avoid singularities
2566 * when three particles align and the dihedral angle and dihedral potential cannot be calculated.
2567 * This potential is calculated using the formula:
2568 * \f[V_{\rm ReB}(\theta_i) = \frac{1}{2} k_{\theta}
2569 * \frac{(\cos\theta_i - \cos\theta_0)^2}{\sin^2\theta_i}\f]
2570 * ({eq:ReB} and ref \cite{MonicaGoga2013} from the manual).
2571 * For more explanations see comments file "restcbt.h" */
2572
2573 compute_factors_restrdihs(type, forceparams,
2574 delta_ante, delta_crnt, delta_post,
2575 &factor_phi_ai_ante, &factor_phi_ai_crnt, &factor_phi_ai_post,
2576 &factor_phi_aj_ante, &factor_phi_aj_crnt, &factor_phi_aj_post,
2577 &factor_phi_ak_ante, &factor_phi_ak_crnt, &factor_phi_ak_post,
2578 &factor_phi_al_ante, &factor_phi_al_crnt, &factor_phi_al_post,
2579 &prefactor_phi, &v);
2580
2581
2582 /* Computation of forces per component */
2583 for (d = 0; d < DIM; d++)
2584 {
2585 f_i[d] = prefactor_phi * (factor_phi_ai_ante * delta_ante[d] + factor_phi_ai_crnt * delta_crnt[d] + factor_phi_ai_post * delta_post[d]);
2586 f_j[d] = prefactor_phi * (factor_phi_aj_ante * delta_ante[d] + factor_phi_aj_crnt * delta_crnt[d] + factor_phi_aj_post * delta_post[d]);
2587 f_k[d] = prefactor_phi * (factor_phi_ak_ante * delta_ante[d] + factor_phi_ak_crnt * delta_crnt[d] + factor_phi_ak_post * delta_post[d]);
2588 f_l[d] = prefactor_phi * (factor_phi_al_ante * delta_ante[d] + factor_phi_al_crnt * delta_crnt[d] + factor_phi_al_post * delta_post[d]);
2589 }
2590 /* Computation of the energy */
2591
2592 vtot += v;
2593
2594
2595
2596 /* Updating the forces */
2597
2598 rvec_inc(f[ai], f_i);
2599 rvec_inc(f[aj], f_j);
2600 rvec_inc(f[ak], f_k);
2601 rvec_inc(f[al], f_l);
2602
2603
2604 /* Updating the fshift forces for the pressure coupling */
2605 if (g)
2606 {
2607 copy_ivec(SHIFT_IVEC(g, aj), jt);
2608 ivec_sub(SHIFT_IVEC(g, ai), jt, dt_ij);
2609 ivec_sub(SHIFT_IVEC(g, ak), jt, dt_kj);
2610 ivec_sub(SHIFT_IVEC(g, al), jt, dt_lj);
2611 t1 = IVEC2IS(dt_ij);
2612 t2 = IVEC2IS(dt_kj);
2613 t3 = IVEC2IS(dt_lj);
2614 }
2615 else if (pbc)
2616 {
2617 t3 = pbc_rvec_sub(pbc, x[al], x[aj], dx_jl);
2618 }
2619 else
2620 {
2621 t3 = CENTRAL;
2622 }
2623
2624 rvec_inc(fshift[t1], f_i);
2625 rvec_inc(fshift[CENTRAL], f_j);
2626 rvec_inc(fshift[t2], f_k);
2627 rvec_inc(fshift[t3], f_l);
2628
2629 }
2630
2631 return vtot;
2632 }
2633
2634
2635 real cbtdihs(int nbonds,
2636 const t_iatom forceatoms[], const t_iparams forceparams[],
2637 const rvec x[], rvec4 f[], rvec fshift[],
2638 const t_pbc *pbc, const t_graph *g,
2639 real gmx_unused lambda, real gmx_unused *dvdlambda,
2640 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
2641 int gmx_unused *global_atom_index)
2642 {
2643 int type, ai, aj, ak, al, i, d;
2644 int t1, t2, t3;
2645 real v, vtot;
2646 rvec vec_temp;
2647 rvec f_i, f_j, f_k, f_l;
2648 ivec jt, dt_ij, dt_kj, dt_lj;
2649 rvec dx_jl;
2650 rvec delta_ante, delta_crnt, delta_post;
2651 rvec f_phi_ai, f_phi_aj, f_phi_ak, f_phi_al;
2652 rvec f_theta_ante_ai, f_theta_ante_aj, f_theta_ante_ak;
2653 rvec f_theta_post_aj, f_theta_post_ak, f_theta_post_al;
2654
2655
2656
2657
2658 vtot = 0.0;
2659 for (i = 0; (i < nbonds); )
2660 {
2661 type = forceatoms[i++];
2662 ai = forceatoms[i++];
2663 aj = forceatoms[i++];
2664 ak = forceatoms[i++];
2665 al = forceatoms[i++];
2666
2667
2668 t1 = pbc_rvec_sub(pbc, x[ai], x[aj], vec_temp);
2669 pbc_rvec_sub(pbc, x[aj], x[ai], delta_ante);
2670 t2 = pbc_rvec_sub(pbc, x[ak], x[aj], vec_temp);
2671 pbc_rvec_sub(pbc, x[ak], x[aj], delta_crnt);
2672 pbc_rvec_sub(pbc, x[ak], x[al], vec_temp);
2673 pbc_rvec_sub(pbc, x[al], x[ak], delta_post);
2674
2675 /* \brief Compute factors for CBT potential
2676 * The combined bending-torsion potential goes to zero in a very smooth manner, eliminating the numerical
2677 * instabilities, when three coarse-grained particles align and the dihedral angle and standard
2678 * dihedral potentials cannot be calculated. The CBT potential is calculated using the formula:
2679 * \f[V_{\rm CBT}(\theta_{i-1}, \theta_i, \phi_i) = k_{\phi} \sin^3\theta_{i-1} \sin^3\theta_{i}
2680 * \sum_{n=0}^4 { a_n \cos^n\phi_i}\f] ({eq:CBT} and ref \cite{MonicaGoga2013} from the manual).
2681 * It contains in its expression not only the dihedral angle \f$\phi\f$
2682 * but also \f[\theta_{i-1}\f] (theta_ante bellow) and \f[\theta_{i}\f] (theta_post bellow)
2683 * --- the adjacent bending angles.
2684 * For more explanations see comments file "restcbt.h". */
2685
2686 compute_factors_cbtdihs(type, forceparams, delta_ante, delta_crnt, delta_post,
2687 f_phi_ai, f_phi_aj, f_phi_ak, f_phi_al,
2688 f_theta_ante_ai, f_theta_ante_aj, f_theta_ante_ak,
2689 f_theta_post_aj, f_theta_post_ak, f_theta_post_al,
2690 &v);
2691
2692
2693 /* Acumulate the resuts per beads */
2694 for (d = 0; d < DIM; d++)
2695 {
2696 f_i[d] = f_phi_ai[d] + f_theta_ante_ai[d];
2697 f_j[d] = f_phi_aj[d] + f_theta_ante_aj[d] + f_theta_post_aj[d];
2698 f_k[d] = f_phi_ak[d] + f_theta_ante_ak[d] + f_theta_post_ak[d];
2699 f_l[d] = f_phi_al[d] + f_theta_post_al[d];
2700 }
2701
2702 /* Compute the potential energy */
2703
2704 vtot += v;
2705
2706
2707 /* Updating the forces */
2708 rvec_inc(f[ai], f_i);
2709 rvec_inc(f[aj], f_j);
2710 rvec_inc(f[ak], f_k);
2711 rvec_inc(f[al], f_l);
2712
2713
2714 /* Updating the fshift forces for the pressure coupling */
2715 if (g)
2716 {
2717 copy_ivec(SHIFT_IVEC(g, aj), jt);
2718 ivec_sub(SHIFT_IVEC(g, ai), jt, dt_ij);
2719 ivec_sub(SHIFT_IVEC(g, ak), jt, dt_kj);
2720 ivec_sub(SHIFT_IVEC(g, al), jt, dt_lj);
2721 t1 = IVEC2IS(dt_ij);
2722 t2 = IVEC2IS(dt_kj);
2723 t3 = IVEC2IS(dt_lj);
2724 }
2725 else if (pbc)
2726 {
2727 t3 = pbc_rvec_sub(pbc, x[al], x[aj], dx_jl);
2728 }
2729 else
2730 {
2731 t3 = CENTRAL;
2732 }
2733
2734 rvec_inc(fshift[t1], f_i);
2735 rvec_inc(fshift[CENTRAL], f_j);
2736 rvec_inc(fshift[t2], f_k);
2737 rvec_inc(fshift[t3], f_l);
2738 }
2739
2740 return vtot;
2741 }
2742
2743 real rbdihs(int nbonds,
2744 const t_iatom forceatoms[], const t_iparams forceparams[],
2745 const rvec x[], rvec4 f[], rvec fshift[],
2746 const t_pbc *pbc, const t_graph *g,
2747 real lambda, real *dvdlambda,
2748 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
2749 int gmx_unused *global_atom_index)
2750 {
2751 const real c0 = 0.0, c1 = 1.0, c2 = 2.0, c3 = 3.0, c4 = 4.0, c5 = 5.0;
2752 int type, ai, aj, ak, al, i, j;
2753 int t1, t2, t3;
2754 rvec r_ij, r_kj, r_kl, m, n;
2755 real parmA[NR_RBDIHS];
2756 real parmB[NR_RBDIHS];
2757 real parm[NR_RBDIHS];
2758 real cos_phi, phi, rbp, rbpBA;
2759 real v, ddphi, sin_phi;
2760 real cosfac, vtot;
2761 real L1 = 1.0-lambda;
2762 real dvdl_term = 0;
2763
2764 vtot = 0.0;
2765 for (i = 0; (i < nbonds); )
2766 {
2767 type = forceatoms[i++];
2768 ai = forceatoms[i++];
2769 aj = forceatoms[i++];
2770 ak = forceatoms[i++];
2771 al = forceatoms[i++];
2772
2773 phi = dih_angle(x[ai], x[aj], x[ak], x[al], pbc, r_ij, r_kj, r_kl, m, n,
2774 &t1, &t2, &t3); /* 84 */
2775
2776 /* Change to polymer convention */
2777 if (phi < c0)
2778 {
2779 phi += M_PI;
2780 }
2781 else
2782 {
2783 phi -= M_PI; /* 1 */
2784
2785 }
2786 cos_phi = std::cos(phi);
2787 /* Beware of accuracy loss, cannot use 1-sqrt(cos^2) ! */
2788 sin_phi = std::sin(phi);
2789
2790 for (j = 0; (j < NR_RBDIHS); j++)
2791 {
2792 parmA[j] = forceparams[type].rbdihs.rbcA[j];
2793 parmB[j] = forceparams[type].rbdihs.rbcB[j];
2794 parm[j] = L1*parmA[j]+lambda*parmB[j];
2795 }
2796 /* Calculate cosine powers */
2797 /* Calculate the energy */
2798 /* Calculate the derivative */
2799
2800 v = parm[0];
2801 dvdl_term += (parmB[0]-parmA[0]);
2802 ddphi = c0;
2803 cosfac = c1;
2804
2805 rbp = parm[1];
2806 rbpBA = parmB[1]-parmA[1];
2807 ddphi += rbp*cosfac;
2808 cosfac *= cos_phi;
2809 v += cosfac*rbp;
2810 dvdl_term += cosfac*rbpBA;
2811 rbp = parm[2];
2812 rbpBA = parmB[2]-parmA[2];
2813 ddphi += c2*rbp*cosfac;
2814 cosfac *= cos_phi;
2815 v += cosfac*rbp;
2816 dvdl_term += cosfac*rbpBA;
2817 rbp = parm[3];
2818 rbpBA = parmB[3]-parmA[3];
2819 ddphi += c3*rbp*cosfac;
2820 cosfac *= cos_phi;
2821 v += cosfac*rbp;
2822 dvdl_term += cosfac*rbpBA;
2823 rbp = parm[4];
2824 rbpBA = parmB[4]-parmA[4];
2825 ddphi += c4*rbp*cosfac;
2826 cosfac *= cos_phi;
2827 v += cosfac*rbp;
2828 dvdl_term += cosfac*rbpBA;
2829 rbp = parm[5];
2830 rbpBA = parmB[5]-parmA[5];
2831 ddphi += c5*rbp*cosfac;
2832 cosfac *= cos_phi;
2833 v += cosfac*rbp;
2834 dvdl_term += cosfac*rbpBA;
2835
2836 ddphi = -ddphi*sin_phi; /* 11 */
2837
2838 do_dih_fup(ai, aj, ak, al, ddphi, r_ij, r_kj, r_kl, m, n,
2839 f, fshift, pbc, g, x, t1, t2, t3); /* 112 */
2840 vtot += v;
2841 }
2842 *dvdlambda += dvdl_term;
2843
2844 return vtot;
2845 }
2846
2847 //! \endcond
2848
2849 /*! \brief Mysterious undocumented function */
2850 static int
2851 cmap_setup_grid_index(int ip, int grid_spacing, int *ipm1, int *ipp1, int *ipp2)
2852 {
2853 int im1, ip1, ip2;
2854
2855 if (ip < 0)
2856 {
2857 ip = ip + grid_spacing - 1;
2858 }
2859 else if (ip > grid_spacing)
2860 {
2861 ip = ip - grid_spacing - 1;
2862 }
2863
2864 im1 = ip - 1;
2865 ip1 = ip + 1;
2866 ip2 = ip + 2;
2867
2868 if (ip == 0)
2869 {
2870 im1 = grid_spacing - 1;
2871 }
2872 else if (ip == grid_spacing-2)
2873 {
2874 ip2 = 0;
2875 }
2876 else if (ip == grid_spacing-1)
2877 {
2878 ip1 = 0;
2879 ip2 = 1;
2880 }
2881
2882 *ipm1 = im1;
2883 *ipp1 = ip1;
2884 *ipp2 = ip2;
2885
2886 return ip;
2887
2888 }
2889
2890 real
2891 cmap_dihs(int nbonds,
2892 const t_iatom forceatoms[], const t_iparams forceparams[],
2893 const gmx_cmap_t *cmap_grid,
2894 const rvec x[], rvec4 f[], rvec fshift[],
2895 const struct t_pbc *pbc, const struct t_graph *g,
2896 real gmx_unused lambda, real gmx_unused *dvdlambda,
2897 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
2898 int gmx_unused *global_atom_index)
2899 {
2900 int i, n;
2901 int ai, aj, ak, al, am;
2902 int a1i, a1j, a1k, a1l, a2i, a2j, a2k, a2l;
2903 int type;
2904 int t11, t21, t31, t12, t22, t32;
2905 int iphi1, ip1m1, ip1p1, ip1p2;
2906 int iphi2, ip2m1, ip2p1, ip2p2;
2907 int l1, l2, l3;
2908 int pos1, pos2, pos3, pos4;
2909
2910 real ty[4], ty1[4], ty2[4], ty12[4], tx[16];
2911 real phi1, cos_phi1, sin_phi1, xphi1;
2912 real phi2, cos_phi2, sin_phi2, xphi2;
2913 real dx, tt, tu, e, df1, df2, vtot;
2914 real ra21, rb21, rg21, rg1, rgr1, ra2r1, rb2r1, rabr1;
2915 real ra22, rb22, rg22, rg2, rgr2, ra2r2, rb2r2, rabr2;
2916 real fg1, hg1, fga1, hgb1, gaa1, gbb1;
2917 real fg2, hg2, fga2, hgb2, gaa2, gbb2;
2918 real fac;
2919
2920 rvec r1_ij, r1_kj, r1_kl, m1, n1;
2921 rvec r2_ij, r2_kj, r2_kl, m2, n2;
2922 rvec f1_i, f1_j, f1_k, f1_l;
2923 rvec f2_i, f2_j, f2_k, f2_l;
2924 rvec a1, b1, a2, b2;
2925 rvec f1, g1, h1, f2, g2, h2;
2926 rvec dtf1, dtg1, dth1, dtf2, dtg2, dth2;
2927 ivec jt1, dt1_ij, dt1_kj, dt1_lj;
2928 ivec jt2, dt2_ij, dt2_kj, dt2_lj;
2929
2930 int loop_index[4][4] = {
2931 {0, 4, 8, 12},
2932 {1, 5, 9, 13},
2933 {2, 6, 10, 14},
2934 {3, 7, 11, 15}
2935 };
2936
2937 /* Total CMAP energy */
2938 vtot = 0;
2939
2940 for (n = 0; n < nbonds; )
2941 {
2942 /* Five atoms are involved in the two torsions */
2943 type = forceatoms[n++];
2944 ai = forceatoms[n++];
2945 aj = forceatoms[n++];
2946 ak = forceatoms[n++];
2947 al = forceatoms[n++];
2948 am = forceatoms[n++];
2949
2950 /* Which CMAP type is this */
2951 const int cmapA = forceparams[type].cmap.cmapA;
2952 const real *cmapd = cmap_grid->cmapdata[cmapA].cmap.data();
2953
2954 /* First torsion */
2955 a1i = ai;
2956 a1j = aj;
2957 a1k = ak;
2958 a1l = al;
2959
2960 phi1 = dih_angle(x[a1i], x[a1j], x[a1k], x[a1l], pbc, r1_ij, r1_kj, r1_kl, m1, n1,
2961 &t11, &t21, &t31); /* 84 */
2962
2963 cos_phi1 = std::cos(phi1);
2964
2965 a1[0] = r1_ij[1]*r1_kj[2]-r1_ij[2]*r1_kj[1];
2966 a1[1] = r1_ij[2]*r1_kj[0]-r1_ij[0]*r1_kj[2];
2967 a1[2] = r1_ij[0]*r1_kj[1]-r1_ij[1]*r1_kj[0]; /* 9 */
2968
2969 b1[0] = r1_kl[1]*r1_kj[2]-r1_kl[2]*r1_kj[1];
2970 b1[1] = r1_kl[2]*r1_kj[0]-r1_kl[0]*r1_kj[2];
2971 b1[2] = r1_kl[0]*r1_kj[1]-r1_kl[1]*r1_kj[0]; /* 9 */
2972
2973 pbc_rvec_sub(pbc, x[a1l], x[a1k], h1);
2974
2975 ra21 = iprod(a1, a1); /* 5 */
2976 rb21 = iprod(b1, b1); /* 5 */
2977 rg21 = iprod(r1_kj, r1_kj); /* 5 */
2978 rg1 = sqrt(rg21);
2979
2980 rgr1 = 1.0/rg1;
2981 ra2r1 = 1.0/ra21;
2982 rb2r1 = 1.0/rb21;
2983 rabr1 = sqrt(ra2r1*rb2r1);
2984
2985 sin_phi1 = rg1 * rabr1 * iprod(a1, h1) * (-1);
2986
2987 if (cos_phi1 < -0.5 || cos_phi1 > 0.5)
2988 {
2989 phi1 = std::asin(sin_phi1);
2990
2991 if (cos_phi1 < 0)
2992 {
2993 if (phi1 > 0)
2994 {
2995 phi1 = M_PI - phi1;
2996 }
2997 else
2998 {
2999 phi1 = -M_PI - phi1;
3000 }
3001 }
3002 }
3003 else
3004 {
3005 phi1 = std::acos(cos_phi1);
3006
3007 if (sin_phi1 < 0)
3008 {
3009 phi1 = -phi1;
3010 }
3011 }
3012
3013 xphi1 = phi1 + M_PI; /* 1 */
3014
3015 /* Second torsion */
3016 a2i = aj;
3017 a2j = ak;
3018 a2k = al;
3019 a2l = am;
3020
3021 phi2 = dih_angle(x[a2i], x[a2j], x[a2k], x[a2l], pbc, r2_ij, r2_kj, r2_kl, m2, n2,
3022 &t12, &t22, &t32); /* 84 */
3023
3024 cos_phi2 = std::cos(phi2);
3025
3026 a2[0] = r2_ij[1]*r2_kj[2]-r2_ij[2]*r2_kj[1];
3027 a2[1] = r2_ij[2]*r2_kj[0]-r2_ij[0]*r2_kj[2];
3028 a2[2] = r2_ij[0]*r2_kj[1]-r2_ij[1]*r2_kj[0]; /* 9 */
3029
3030 b2[0] = r2_kl[1]*r2_kj[2]-r2_kl[2]*r2_kj[1];
3031 b2[1] = r2_kl[2]*r2_kj[0]-r2_kl[0]*r2_kj[2];
3032 b2[2] = r2_kl[0]*r2_kj[1]-r2_kl[1]*r2_kj[0]; /* 9 */
3033
3034 pbc_rvec_sub(pbc, x[a2l], x[a2k], h2);
3035
3036 ra22 = iprod(a2, a2); /* 5 */
3037 rb22 = iprod(b2, b2); /* 5 */
3038 rg22 = iprod(r2_kj, r2_kj); /* 5 */
3039 rg2 = sqrt(rg22);
3040
3041 rgr2 = 1.0/rg2;
3042 ra2r2 = 1.0/ra22;
3043 rb2r2 = 1.0/rb22;
3044 rabr2 = sqrt(ra2r2*rb2r2);
3045
3046 sin_phi2 = rg2 * rabr2 * iprod(a2, h2) * (-1);
3047
3048 if (cos_phi2 < -0.5 || cos_phi2 > 0.5)
3049 {
3050 phi2 = std::asin(sin_phi2);
3051
3052 if (cos_phi2 < 0)
3053 {
3054 if (phi2 > 0)
3055 {
3056 phi2 = M_PI - phi2;
3057 }
3058 else
3059 {
3060 phi2 = -M_PI - phi2;
3061 }
3062 }
3063 }
3064 else
3065 {
3066 phi2 = std::acos(cos_phi2);
3067
3068 if (sin_phi2 < 0)
3069 {
3070 phi2 = -phi2;
3071 }
3072 }
3073
3074 xphi2 = phi2 + M_PI; /* 1 */
3075
3076 /* Range mangling */
3077 if (xphi1 < 0)
3078 {
3079 xphi1 = xphi1 + 2*M_PI;
3080 }
3081 else if (xphi1 >= 2*M_PI)
3082 {
3083 xphi1 = xphi1 - 2*M_PI;
3084 }
3085
3086 if (xphi2 < 0)
3087 {
3088 xphi2 = xphi2 + 2*M_PI;
3089 }
3090 else if (xphi2 >= 2*M_PI)
3091 {
3092 xphi2 = xphi2 - 2*M_PI;
3093 }
3094
3095 /* Number of grid points */
3096 dx = 2*M_PI / cmap_grid->grid_spacing;
3097
3098 /* Where on the grid are we */
3099 iphi1 = static_cast<int>(xphi1/dx);
3100 iphi2 = static_cast<int>(xphi2/dx);
3101
3102 iphi1 = cmap_setup_grid_index(iphi1, cmap_grid->grid_spacing, &ip1m1, &ip1p1, &ip1p2);
3103 iphi2 = cmap_setup_grid_index(iphi2, cmap_grid->grid_spacing, &ip2m1, &ip2p1, &ip2p2);
3104
3105 pos1 = iphi1*cmap_grid->grid_spacing+iphi2;
3106 pos2 = ip1p1*cmap_grid->grid_spacing+iphi2;
3107 pos3 = ip1p1*cmap_grid->grid_spacing+ip2p1;
3108 pos4 = iphi1*cmap_grid->grid_spacing+ip2p1;
3109
3110 ty[0] = cmapd[pos1*4];
3111 ty[1] = cmapd[pos2*4];
3112 ty[2] = cmapd[pos3*4];
3113 ty[3] = cmapd[pos4*4];
3114
3115 ty1[0] = cmapd[pos1*4+1];
3116 ty1[1] = cmapd[pos2*4+1];
3117 ty1[2] = cmapd[pos3*4+1];
3118 ty1[3] = cmapd[pos4*4+1];
3119
3120 ty2[0] = cmapd[pos1*4+2];
3121 ty2[1] = cmapd[pos2*4+2];
3122 ty2[2] = cmapd[pos3*4+2];
3123 ty2[3] = cmapd[pos4*4+2];
3124
3125 ty12[0] = cmapd[pos1*4+3];
3126 ty12[1] = cmapd[pos2*4+3];
3127 ty12[2] = cmapd[pos3*4+3];
3128 ty12[3] = cmapd[pos4*4+3];
3129
3130 /* Switch to degrees */
3131 dx = 360.0 / cmap_grid->grid_spacing;
3132 xphi1 = xphi1 * RAD2DEG;
3133 xphi2 = xphi2 * RAD2DEG;
3134
3135 for (i = 0; i < 4; i++) /* 16 */
3136 {
3137 tx[i] = ty[i];
3138 tx[i+4] = ty1[i]*dx;
3139 tx[i+8] = ty2[i]*dx;
3140 tx[i+12] = ty12[i]*dx*dx;
3141 }
3142
3143 real tc[16] = {0};
3144 for (int idx = 0; idx < 16; idx++) /* 1056 */
3145 {
3146 for (int k = 0; k < 16; k++)
3147 {
3148 tc[idx] += cmap_coeff_matrix[k*16+idx]*tx[k];
3149 }
3150 }
3151
3152 tt = (xphi1-iphi1*dx)/dx;
3153 tu = (xphi2-iphi2*dx)/dx;
3154
3155 e = 0;
3156 df1 = 0;
3157 df2 = 0;
3158
3159 for (i = 3; i >= 0; i--)
3160 {
3161 l1 = loop_index[i][3];
3162 l2 = loop_index[i][2];
3163 l3 = loop_index[i][1];
3164
3165 e = tt * e + ((tc[i*4+3]*tu+tc[i*4+2])*tu + tc[i*4+1])*tu+tc[i*4];
3166 df1 = tu * df1 + (3.0*tc[l1]*tt+2.0*tc[l2])*tt+tc[l3];
3167 df2 = tt * df2 + (3.0*tc[i*4+3]*tu+2.0*tc[i*4+2])*tu+tc[i*4+1];
3168 }
3169
3170 fac = RAD2DEG/dx;
3171 df1 = df1 * fac;
3172 df2 = df2 * fac;
3173
3174 /* CMAP energy */
3175 vtot += e;
3176
3177 /* Do forces - first torsion */
3178 fg1 = iprod(r1_ij, r1_kj);
3179 hg1 = iprod(r1_kl, r1_kj);
3180 fga1 = fg1*ra2r1*rgr1;
3181 hgb1 = hg1*rb2r1*rgr1;
3182 gaa1 = -ra2r1*rg1;
3183 gbb1 = rb2r1*rg1;
3184
3185 for (i = 0; i < DIM; i++)
3186 {
3187 dtf1[i] = gaa1 * a1[i];
3188 dtg1[i] = fga1 * a1[i] - hgb1 * b1[i];
3189 dth1[i] = gbb1 * b1[i];
3190
3191 f1[i] = df1 * dtf1[i];
3192 g1[i] = df1 * dtg1[i];
3193 h1[i] = df1 * dth1[i];
3194
3195 f1_i[i] = f1[i];
3196 f1_j[i] = -f1[i] - g1[i];
3197 f1_k[i] = h1[i] + g1[i];
3198 f1_l[i] = -h1[i];
3199
3200 f[a1i][i] = f[a1i][i] + f1_i[i];
3201 f[a1j][i] = f[a1j][i] + f1_j[i]; /* - f1[i] - g1[i] */
3202 f[a1k][i] = f[a1k][i] + f1_k[i]; /* h1[i] + g1[i] */
3203 f[a1l][i] = f[a1l][i] + f1_l[i]; /* h1[i] */
3204 }
3205
3206 /* Do forces - second torsion */
3207 fg2 = iprod(r2_ij, r2_kj);
3208 hg2 = iprod(r2_kl, r2_kj);
3209 fga2 = fg2*ra2r2*rgr2;
3210 hgb2 = hg2*rb2r2*rgr2;
3211 gaa2 = -ra2r2*rg2;
3212 gbb2 = rb2r2*rg2;
3213
3214 for (i = 0; i < DIM; i++)
3215 {
3216 dtf2[i] = gaa2 * a2[i];
3217 dtg2[i] = fga2 * a2[i] - hgb2 * b2[i];
3218 dth2[i] = gbb2 * b2[i];
3219
3220 f2[i] = df2 * dtf2[i];
3221 g2[i] = df2 * dtg2[i];
3222 h2[i] = df2 * dth2[i];
3223
3224 f2_i[i] = f2[i];
3225 f2_j[i] = -f2[i] - g2[i];
3226 f2_k[i] = h2[i] + g2[i];
3227 f2_l[i] = -h2[i];
3228
3229 f[a2i][i] = f[a2i][i] + f2_i[i]; /* f2[i] */
3230 f[a2j][i] = f[a2j][i] + f2_j[i]; /* - f2[i] - g2[i] */
3231 f[a2k][i] = f[a2k][i] + f2_k[i]; /* h2[i] + g2[i] */
3232 f[a2l][i] = f[a2l][i] + f2_l[i]; /* - h2[i] */
3233 }
3234
3235 /* Shift forces */
3236 if (g)
3237 {
3238 copy_ivec(SHIFT_IVEC(g, a1j), jt1);
3239 ivec_sub(SHIFT_IVEC(g, a1i), jt1, dt1_ij);
3240 ivec_sub(SHIFT_IVEC(g, a1k), jt1, dt1_kj);
3241 ivec_sub(SHIFT_IVEC(g, a1l), jt1, dt1_lj);
3242 t11 = IVEC2IS(dt1_ij);
3243 t21 = IVEC2IS(dt1_kj);
3244 t31 = IVEC2IS(dt1_lj);
3245
3246 copy_ivec(SHIFT_IVEC(g, a2j), jt2);
3247 ivec_sub(SHIFT_IVEC(g, a2i), jt2, dt2_ij);
3248 ivec_sub(SHIFT_IVEC(g, a2k), jt2, dt2_kj);
3249 ivec_sub(SHIFT_IVEC(g, a2l), jt2, dt2_lj);
3250 t12 = IVEC2IS(dt2_ij);
3251 t22 = IVEC2IS(dt2_kj);
3252 t32 = IVEC2IS(dt2_lj);
3253 }
3254 else if (pbc)
3255 {
3256 t31 = pbc_rvec_sub(pbc, x[a1l], x[a1j], h1);
3257 t32 = pbc_rvec_sub(pbc, x[a2l], x[a2j], h2);
3258 }
3259 else
3260 {
3261 t31 = CENTRAL;
3262 t32 = CENTRAL;
3263 }
3264
3265 rvec_inc(fshift[t11], f1_i);
3266 rvec_inc(fshift[CENTRAL], f1_j);
3267 rvec_inc(fshift[t21], f1_k);
3268 rvec_inc(fshift[t31], f1_l);
3269
3270 rvec_inc(fshift[t21], f2_i);
3271 rvec_inc(fshift[CENTRAL], f2_j);
3272 rvec_inc(fshift[t22], f2_k);
3273 rvec_inc(fshift[t32], f2_l);
3274 }
3275 return vtot;
3276 }
3277
3278
3279 //! \cond
3280 /***********************************************************
3281 *
3282 * G R O M O S 9 6 F U N C T I O N S
3283 *
3284 ***********************************************************/
3285 static real g96harmonic(real kA, real kB, real xA, real xB, real x, real lambda,
3286 real *V, real *F)
3287 {
3288 const real half = 0.5;
3289 real L1, kk, x0, dx, dx2;
3290 real v, f, dvdlambda;
3291
3292 L1 = 1.0-lambda;
3293 kk = L1*kA+lambda*kB;
3294 x0 = L1*xA+lambda*xB;
3295
3296 dx = x-x0;
3297 dx2 = dx*dx;
3298
3299 f = -kk*dx;
3300 v = half*kk*dx2;
3301 dvdlambda = half*(kB-kA)*dx2 + (xA-xB)*kk*dx;
3302
3303 *F = f;
3304 *V = v;
3305
3306 return dvdlambda;
3307
3308 /* That was 21 flops */
3309 }
3310
3311 real g96bonds(int nbonds,
3312 const t_iatom forceatoms[], const t_iparams forceparams[],
3313 const rvec x[], rvec4 f[], rvec fshift[],
3314 const t_pbc *pbc, const t_graph *g,
3315 real lambda, real *dvdlambda,
3316 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
3317 int gmx_unused *global_atom_index)
3318 {
3319 int i, m, ki, ai, aj, type;
3320 real dr2, fbond, vbond, fij, vtot;
3321 rvec dx;
3322 ivec dt;
3323
3324 vtot = 0.0;
3325 for (i = 0; (i < nbonds); )
3326 {
3327 type = forceatoms[i++];
3328 ai = forceatoms[i++];
3329 aj = forceatoms[i++];
3330
3331 ki = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /* 3 */
3332 dr2 = iprod(dx, dx); /* 5 */
3333
3334 *dvdlambda += g96harmonic(forceparams[type].harmonic.krA,
3335 forceparams[type].harmonic.krB,
3336 forceparams[type].harmonic.rA,
3337 forceparams[type].harmonic.rB,
3338 dr2, lambda, &vbond, &fbond);
3339
3340 vtot += 0.5*vbond; /* 1*/
3341
3342 if (g)
3343 {
3344 ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, aj), dt);
3345 ki = IVEC2IS(dt);
3346 }
3347 for (m = 0; (m < DIM); m++) /* 15 */
3348 {
3349 fij = fbond*dx[m];
3350 f[ai][m] += fij;
3351 f[aj][m] -= fij;
3352 fshift[ki][m] += fij;
3353 fshift[CENTRAL][m] -= fij;
3354 }
3355 } /* 44 TOTAL */
3356 return vtot;
3357 }
3358
3359 static real g96bond_angle(const rvec xi, const rvec xj, const rvec xk, const t_pbc *pbc,
3360 rvec r_ij, rvec r_kj,
3361 int *t1, int *t2)
3362 /* Return value is the angle between the bonds i-j and j-k */
3363 {
3364 real costh;
3365
3366 *t1 = pbc_rvec_sub(pbc, xi, xj, r_ij); /* 3 */
3367 *t2 = pbc_rvec_sub(pbc, xk, xj, r_kj); /* 3 */
3368
3369 costh = cos_angle(r_ij, r_kj); /* 25 */
3370 /* 41 TOTAL */
3371 return costh;
3372 }
3373
3374 real g96angles(int nbonds,
3375 const t_iatom forceatoms[], const t_iparams forceparams[],
3376 const rvec x[], rvec4 f[], rvec fshift[],
3377 const t_pbc *pbc, const t_graph *g,
3378 real lambda, real *dvdlambda,
3379 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
3380 int gmx_unused *global_atom_index)
3381 {
3382 int i, ai, aj, ak, type, m, t1, t2;
3383 rvec r_ij, r_kj;
3384 real cos_theta, dVdt, va, vtot;
3385 real rij_1, rij_2, rkj_1, rkj_2, rijrkj_1;
3386 rvec f_i, f_j, f_k;
3387 ivec jt, dt_ij, dt_kj;
3388
3389 vtot = 0.0;
3390 for (i = 0; (i < nbonds); )
3391 {
3392 type = forceatoms[i++];
3393 ai = forceatoms[i++];
3394 aj = forceatoms[i++];
3395 ak = forceatoms[i++];
3396
3397 cos_theta = g96bond_angle(x[ai], x[aj], x[ak], pbc, r_ij, r_kj, &t1, &t2);
3398
3399 *dvdlambda += g96harmonic(forceparams[type].harmonic.krA,
3400 forceparams[type].harmonic.krB,
3401 forceparams[type].harmonic.rA,
3402 forceparams[type].harmonic.rB,
3403 cos_theta, lambda, &va, &dVdt);
3404 vtot += va;
3405
3406 rij_1 = gmx::invsqrt(iprod(r_ij, r_ij));
3407 rkj_1 = gmx::invsqrt(iprod(r_kj, r_kj));
3408 rij_2 = rij_1*rij_1;
3409 rkj_2 = rkj_1*rkj_1;
3410 rijrkj_1 = rij_1*rkj_1; /* 23 */
3411
3412 for (m = 0; (m < DIM); m++) /* 42 */
3413 {
3414 f_i[m] = dVdt*(r_kj[m]*rijrkj_1 - r_ij[m]*rij_2*cos_theta);
3415 f_k[m] = dVdt*(r_ij[m]*rijrkj_1 - r_kj[m]*rkj_2*cos_theta);
3416 f_j[m] = -f_i[m]-f_k[m];
3417 f[ai][m] += f_i[m];
3418 f[aj][m] += f_j[m];
3419 f[ak][m] += f_k[m];
3420 }
3421
3422 if (g)
3423 {
3424 copy_ivec(SHIFT_IVEC(g, aj), jt);
3425
3426 ivec_sub(SHIFT_IVEC(g, ai), jt, dt_ij);
3427 ivec_sub(SHIFT_IVEC(g, ak), jt, dt_kj);
3428 t1 = IVEC2IS(dt_ij);
3429 t2 = IVEC2IS(dt_kj);
3430 }
3431 rvec_inc(fshift[t1], f_i);
3432 rvec_inc(fshift[CENTRAL], f_j);
3433 rvec_inc(fshift[t2], f_k); /* 9 */
3434 /* 163 TOTAL */
3435 }
3436 return vtot;
3437 }
3438
3439 real cross_bond_bond(int nbonds,
3440 const t_iatom forceatoms[], const t_iparams forceparams[],
3441 const rvec x[], rvec4 f[], rvec fshift[],
3442 const t_pbc *pbc, const t_graph *g,
3443 real gmx_unused lambda, real gmx_unused *dvdlambda,
3444 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
3445 int gmx_unused *global_atom_index)
3446 {
3447 /* Potential from Lawrence and Skimmer, Chem. Phys. Lett. 372 (2003)
3448 * pp. 842-847
3449 */
3450 int i, ai, aj, ak, type, m, t1, t2;
3451 rvec r_ij, r_kj;
3452 real vtot, vrr, s1, s2, r1, r2, r1e, r2e, krr;
3453 rvec f_i, f_j, f_k;
3454 ivec jt, dt_ij, dt_kj;
3455
3456 vtot = 0.0;
3457 for (i = 0; (i < nbonds); )
3458 {
3459 type = forceatoms[i++];
3460 ai = forceatoms[i++];
3461 aj = forceatoms[i++];
3462 ak = forceatoms[i++];
3463 r1e = forceparams[type].cross_bb.r1e;
3464 r2e = forceparams[type].cross_bb.r2e;
3465 krr = forceparams[type].cross_bb.krr;
3466
3467 /* Compute distance vectors ... */
3468 t1 = pbc_rvec_sub(pbc, x[ai], x[aj], r_ij);
3469 t2 = pbc_rvec_sub(pbc, x[ak], x[aj], r_kj);
3470
3471 /* ... and their lengths */
3472 r1 = norm(r_ij);
3473 r2 = norm(r_kj);
3474
3475 /* Deviations from ideality */
3476 s1 = r1-r1e;
3477 s2 = r2-r2e;
3478
3479 /* Energy (can be negative!) */
3480 vrr = krr*s1*s2;
3481 vtot += vrr;
3482
3483 /* Forces */
3484 svmul(-krr*s2/r1, r_ij, f_i);
3485 svmul(-krr*s1/r2, r_kj, f_k);
3486
3487 for (m = 0; (m < DIM); m++) /* 12 */
3488 {
3489 f_j[m] = -f_i[m] - f_k[m];
3490 f[ai][m] += f_i[m];
3491 f[aj][m] += f_j[m];
3492 f[ak][m] += f_k[m];
3493 }
3494
3495 /* Virial stuff */
3496 if (g)
3497 {
3498 copy_ivec(SHIFT_IVEC(g, aj), jt);
3499
3500 ivec_sub(SHIFT_IVEC(g, ai), jt, dt_ij);
3501 ivec_sub(SHIFT_IVEC(g, ak), jt, dt_kj);
3502 t1 = IVEC2IS(dt_ij);
3503 t2 = IVEC2IS(dt_kj);
3504 }
3505 rvec_inc(fshift[t1], f_i);
3506 rvec_inc(fshift[CENTRAL], f_j);
3507 rvec_inc(fshift[t2], f_k); /* 9 */
3508 /* 163 TOTAL */
3509 }
3510 return vtot;
3511 }
3512
3513 real cross_bond_angle(int nbonds,
3514 const t_iatom forceatoms[], const t_iparams forceparams[],
3515 const rvec x[], rvec4 f[], rvec fshift[],
3516 const t_pbc *pbc, const t_graph *g,
3517 real gmx_unused lambda, real gmx_unused *dvdlambda,
3518 const t_mdatoms gmx_unused *md, t_fcdata gmx_unused *fcd,
3519 int gmx_unused *global_atom_index)
3520 {
3521 /* Potential from Lawrence and Skimmer, Chem. Phys. Lett. 372 (2003)
3522 * pp. 842-847
3523 */
3524 int i, ai, aj, ak, type, m, t1, t2;
3525 rvec r_ij, r_kj, r_ik;
3526 real vtot, vrt, s1, s2, s3, r1, r2, r3, r1e, r2e, r3e, krt, k1, k2, k3;
3527 rvec f_i, f_j, f_k;
3528 ivec jt, dt_ij, dt_kj;
3529
3530 vtot = 0.0;
3531 for (i = 0; (i < nbonds); )
3532 {
3533 type = forceatoms[i++];
3534 ai = forceatoms[i++];
3535 aj = forceatoms[i++];
3536 ak = forceatoms[i++];
3537 r1e = forceparams[type].cross_ba.r1e;
3538 r2e = forceparams[type].cross_ba.r2e;
3539 r3e = forceparams[type].cross_ba.r3e;
3540 krt = forceparams[type].cross_ba.krt;
3541
3542 /* Compute distance vectors ... */
3543 t1 = pbc_rvec_sub(pbc, x[ai], x[aj], r_ij);
3544 t2 = pbc_rvec_sub(pbc, x[ak], x[aj], r_kj);
3545 pbc_rvec_sub(pbc, x[ai], x[ak], r_ik);
3546
3547 /* ... and their lengths */
3548 r1 = norm(r_ij);
3549 r2 = norm(r_kj);
3550 r3 = norm(r_ik);
3551
3552 /* Deviations from ideality */
3553 s1 = r1-r1e;
3554 s2 = r2-r2e;
3555 s3 = r3-r3e;
3556
3557 /* Energy (can be negative!) */
3558 vrt = krt*s3*(s1+s2);
3559 vtot += vrt;
3560
3561 /* Forces */
3562 k1 = -krt*(s3/r1);
3563 k2 = -krt*(s3/r2);
3564 k3 = -krt*(s1+s2)/r3;
3565 for (m = 0; (m < DIM); m++)
3566 {
3567 f_i[m] = k1*r_ij[m] + k3*r_ik[m];
3568 f_k[m] = k2*r_kj[m] - k3*r_ik[m];
3569 f_j[m] = -f_i[m] - f_k[m];
3570 }
3571
3572 for (m = 0; (m < DIM); m++) /* 12 */
3573 {
3574 f[ai][m] += f_i[m];
3575 f[aj][m] += f_j[m];
3576 f[ak][m] += f_k[m];
3577 }
3578
3579 /* Virial stuff */
3580 if (g)
3581 {
3582 copy_ivec(SHIFT_IVEC(g, aj), jt);
3583
3584 ivec_sub(SHIFT_IVEC(g, ai), jt, dt_ij);
3585 ivec_sub(SHIFT_IVEC(g, ak), jt, dt_kj);
3586 t1 = IVEC2IS(dt_ij);
3587 t2 = IVEC2IS(dt_kj);
3588 }
3589 rvec_inc(fshift[t1], f_i);
3590 rvec_inc(fshift[CENTRAL], f_j);
3591 rvec_inc(fshift[t2], f_k); /* 9 */
3592 /* 163 TOTAL */
3593 }
3594 return vtot;
3595 }
3596
3597 static real bonded_tab(const char *type, int table_nr,
3598 const bondedtable_t *table, real kA, real kB, real r,
3599 real lambda, real *V, real *F)
3600 {
3601 real k, tabscale, *VFtab, rt, eps, eps2, Yt, Ft, Geps, Heps2, Fp, VV, FF;
3602 int n0, nnn;
3603 real dvdlambda;
3604
3605 k = (1.0 - lambda)*kA + lambda*kB;
3606
3607 tabscale = table->scale;
3608 VFtab = table->data;
3609
3610 rt = r*tabscale;
3611 n0 = static_cast<int>(rt);
3612 if (n0 >= table->n)
3613 {
3614 gmx_fatal(FARGS, "A tabulated %s interaction table number %d is out of the table range: r %f, between table indices %d and %d, table length %d",
3615 type, table_nr, r, n0, n0+1, table->n);
3616 }
3617 eps = rt - n0;
3618 eps2 = eps*eps;
3619 nnn = 4*n0;
3620 Yt = VFtab[nnn];
3621 Ft = VFtab[nnn+1];
3622 Geps = VFtab[nnn+2]*eps;
3623 Heps2 = VFtab[nnn+3]*eps2;
3624 Fp = Ft + Geps + Heps2;
3625 VV = Yt + Fp*eps;
3626 FF = Fp + Geps + 2.0*Heps2;
3627
3628 *F = -k*FF*tabscale;
3629 *V = k*VV;
3630 dvdlambda = (kB - kA)*VV;
3631
3632 return dvdlambda;
3633
3634 /* That was 22 flops */
3635 }
3636
3637 real tab_bonds(int nbonds,
3638 const t_iatom forceatoms[], const t_iparams forceparams[],
3639 const rvec x[], rvec4 f[], rvec fshift[],
3640 const t_pbc *pbc, const t_graph *g,
3641 real lambda, real *dvdlambda,
3642 const t_mdatoms gmx_unused *md, t_fcdata *fcd,
3643 int gmx_unused *global_atom_index)
3644 {
3645 int i, m, ki, ai, aj, type, table;
3646 real dr, dr2, fbond, vbond, fij, vtot;
3647 rvec dx;
3648 ivec dt;
3649
3650 vtot = 0.0;
3651 for (i = 0; (i < nbonds); )
3652 {
3653 type = forceatoms[i++];
3654 ai = forceatoms[i++];
3655 aj = forceatoms[i++];
3656
3657 ki = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /* 3 */
3658 dr2 = iprod(dx, dx); /* 5 */
3659 dr = dr2*gmx::invsqrt(dr2); /* 10 */
3660
3661 table = forceparams[type].tab.table;
3662
3663 *dvdlambda += bonded_tab("bond", table,
3664 &fcd->bondtab[table],
3665 forceparams[type].tab.kA,
3666 forceparams[type].tab.kB,
3667 dr, lambda, &vbond, &fbond); /* 22 */
3668
3669 if (dr2 == 0.0)
3670 {
3671 continue;
3672 }
3673
3674
3675 vtot += vbond; /* 1*/
3676 fbond *= gmx::invsqrt(dr2); /* 6 */
3677 if (g)
3678 {
3679 ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, aj), dt);
3680 ki = IVEC2IS(dt);
3681 }
3682 for (m = 0; (m < DIM); m++) /* 15 */
3683 {
3684 fij = fbond*dx[m];
3685 f[ai][m] += fij;
3686 f[aj][m] -= fij;
3687 fshift[ki][m] += fij;
3688 fshift[CENTRAL][m] -= fij;
3689 }
3690 } /* 62 TOTAL */
3691 return vtot;
3692 }
3693
3694 real tab_angles(int nbonds,
3695 const t_iatom forceatoms[], const t_iparams forceparams[],
3696 const rvec x[], rvec4 f[], rvec fshift[],
3697 const t_pbc *pbc, const t_graph *g,
3698 real lambda, real *dvdlambda,
3699 const t_mdatoms gmx_unused *md, t_fcdata *fcd,
3700 int gmx_unused *global_atom_index)
3701 {
3702 int i, ai, aj, ak, t1, t2, type, table;
3703 rvec r_ij, r_kj;
3704 real cos_theta, cos_theta2, theta, dVdt, va, vtot;
3705 ivec jt, dt_ij, dt_kj;
3706
3707 vtot = 0.0;
3708 for (i = 0; (i < nbonds); )
3709 {
3710 type = forceatoms[i++];
3711 ai = forceatoms[i++];
3712 aj = forceatoms[i++];
3713 ak = forceatoms[i++];
3714
3715 theta = bond_angle(x[ai], x[aj], x[ak], pbc,
3716 r_ij, r_kj, &cos_theta, &t1, &t2); /* 41 */
3717
3718 table = forceparams[type].tab.table;
3719
3720 *dvdlambda += bonded_tab("angle", table,
3721 &fcd->angletab[table],
3722 forceparams[type].tab.kA,
3723 forceparams[type].tab.kB,
3724 theta, lambda, &va, &dVdt); /* 22 */
3725 vtot += va;
3726
3727 cos_theta2 = gmx::square(cos_theta); /* 1 */
3728 if (cos_theta2 < 1)
3729 {
3730 int m;
3731 real st, sth;
3732 real cik, cii, ckk;
3733 real nrkj2, nrij2;
3734 rvec f_i, f_j, f_k;
3735
3736 st = dVdt*gmx::invsqrt(1 - cos_theta2); /* 12 */
3737 sth = st*cos_theta; /* 1 */
3738 nrkj2 = iprod(r_kj, r_kj); /* 5 */
3739 nrij2 = iprod(r_ij, r_ij);
3740
3741 cik = st*gmx::invsqrt(nrkj2*nrij2); /* 12 */
3742 cii = sth/nrij2; /* 10 */
3743 ckk = sth/nrkj2; /* 10 */
3744
3745 for (m = 0; (m < DIM); m++) /* 39 */
3746 {
3747 f_i[m] = -(cik*r_kj[m]-cii*r_ij[m]);
3748 f_k[m] = -(cik*r_ij[m]-ckk*r_kj[m]);
3749 f_j[m] = -f_i[m]-f_k[m];
3750 f[ai][m] += f_i[m];
3751 f[aj][m] += f_j[m];
3752 f[ak][m] += f_k[m];
3753 }
3754 if (g)
3755 {
3756 copy_ivec(SHIFT_IVEC(g, aj), jt);
3757
3758 ivec_sub(SHIFT_IVEC(g, ai), jt, dt_ij);
3759 ivec_sub(SHIFT_IVEC(g, ak), jt, dt_kj);
3760 t1 = IVEC2IS(dt_ij);
3761 t2 = IVEC2IS(dt_kj);
3762 }
3763 rvec_inc(fshift[t1], f_i);
3764 rvec_inc(fshift[CENTRAL], f_j);
3765 rvec_inc(fshift[t2], f_k);
3766 } /* 169 TOTAL */
3767 }
3768 return vtot;
3769 }
3770
3771 real tab_dihs(int nbonds,
3772 const t_iatom forceatoms[], const t_iparams forceparams[],
3773 const rvec x[], rvec4 f[], rvec fshift[],
3774 const t_pbc *pbc, const t_graph *g,
3775 real lambda, real *dvdlambda,
3776 const t_mdatoms gmx_unused *md, t_fcdata *fcd,
3777 int gmx_unused *global_atom_index)
3778 {
3779 int i, type, ai, aj, ak, al, table;
3780 int t1, t2, t3;
3781 rvec r_ij, r_kj, r_kl, m, n;
3782 real phi, ddphi, vpd, vtot;
3783
3784 vtot = 0.0;
3785 for (i = 0; (i < nbonds); )
3786 {
3787 type = forceatoms[i++];
3788 ai = forceatoms[i++];
3789 aj = forceatoms[i++];
3790 ak = forceatoms[i++];
3791 al = forceatoms[i++];
3792
3793 phi = dih_angle(x[ai], x[aj], x[ak], x[al], pbc, r_ij, r_kj, r_kl, m, n,
3794 &t1, &t2, &t3); /* 84 */
3795
3796 table = forceparams[type].tab.table;
3797
3798 /* Hopefully phi+M_PI never results in values < 0 */
3799 *dvdlambda += bonded_tab("dihedral", table,
3800 &fcd->dihtab[table],
3801 forceparams[type].tab.kA,
3802 forceparams[type].tab.kB,
3803 phi+M_PI, lambda, &vpd, &ddphi);
3804
3805 vtot += vpd;
3806 do_dih_fup(ai, aj, ak, al, -ddphi, r_ij, r_kj, r_kl, m, n,
3807 f, fshift, pbc, g, x, t1, t2, t3); /* 112 */
3808
3809 } /* 227 TOTAL */
3810
3811 return vtot;
3812 }
3813
3814 //! \endcond
The ultimate molecular dynamics simulation package