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