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