File: | gromacs/gmxlib/bondfree.c |
Location: | line 3818, column 9 |
Description: | Value stored to 't3' is never read |
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, 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 | #ifdef HAVE_CONFIG_H1 |
38 | #include <config.h> |
39 | #endif |
40 | |
41 | #include <math.h> |
42 | #include <assert.h> |
43 | #include "physics.h" |
44 | #include "gromacs/math/vec.h" |
45 | #include "gromacs/math/utilities.h" |
46 | #include "txtdump.h" |
47 | #include "bondf.h" |
48 | #include "gromacs/utility/smalloc.h" |
49 | #include "pbc.h" |
50 | #include "ns.h" |
51 | #include "macros.h" |
52 | #include "names.h" |
53 | #include "gromacs/utility/fatalerror.h" |
54 | #include "mshift.h" |
55 | #include "disre.h" |
56 | #include "orires.h" |
57 | #include "force.h" |
58 | #include "nonbonded.h" |
59 | #include "restcbt.h" |
60 | |
61 | #include "gromacs/simd/simd.h" |
62 | #include "gromacs/simd/simd_math.h" |
63 | #include "gromacs/simd/vector_operations.h" |
64 | |
65 | /* Find a better place for this? */ |
66 | const int cmap_coeff_matrix[] = { |
67 | 1, 0, -3, 2, 0, 0, 0, 0, -3, 0, 9, -6, 2, 0, -6, 4, |
68 | 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, -9, 6, -2, 0, 6, -4, |
69 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 9, -6, 0, 0, -6, 4, |
70 | 0, 0, 3, -2, 0, 0, 0, 0, 0, 0, -9, 6, 0, 0, 6, -4, |
71 | 0, 0, 0, 0, 1, 0, -3, 2, -2, 0, 6, -4, 1, 0, -3, 2, |
72 | 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 3, -2, 1, 0, -3, 2, |
73 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 2, 0, 0, 3, -2, |
74 | 0, 0, 0, 0, 0, 0, 3, -2, 0, 0, -6, 4, 0, 0, 3, -2, |
75 | 0, 1, -2, 1, 0, 0, 0, 0, 0, -3, 6, -3, 0, 2, -4, 2, |
76 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, -6, 3, 0, -2, 4, -2, |
77 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 3, 0, 0, 2, -2, |
78 | 0, 0, -1, 1, 0, 0, 0, 0, 0, 0, 3, -3, 0, 0, -2, 2, |
79 | 0, 0, 0, 0, 0, 1, -2, 1, 0, -2, 4, -2, 0, 1, -2, 1, |
80 | 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 2, -1, 0, 1, -2, 1, |
81 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, -1, 1, |
82 | 0, 0, 0, 0, 0, 0, -1, 1, 0, 0, 2, -2, 0, 0, -1, 1 |
83 | }; |
84 | |
85 | |
86 | |
87 | int glatnr(int *global_atom_index, int i) |
88 | { |
89 | int atnr; |
90 | |
91 | if (global_atom_index == NULL((void*)0)) |
92 | { |
93 | atnr = i + 1; |
94 | } |
95 | else |
96 | { |
97 | atnr = global_atom_index[i] + 1; |
98 | } |
99 | |
100 | return atnr; |
101 | } |
102 | |
103 | static int pbc_rvec_sub(const t_pbc *pbc, const rvec xi, const rvec xj, rvec dx) |
104 | { |
105 | if (pbc) |
106 | { |
107 | return pbc_dx_aiuc(pbc, xi, xj, dx); |
108 | } |
109 | else |
110 | { |
111 | rvec_sub(xi, xj, dx); |
112 | return CENTRAL(((2*1 +1)*(2*1 +1)*(2*2 +1))/2); |
113 | } |
114 | } |
115 | |
116 | #ifdef GMX_SIMD_HAVE_REAL |
117 | |
118 | /* SIMD PBC data structure, containing 1/boxdiag and the box vectors */ |
119 | typedef struct { |
120 | gmx_simd_real_t__m128 inv_bzz; |
121 | gmx_simd_real_t__m128 inv_byy; |
122 | gmx_simd_real_t__m128 inv_bxx; |
123 | gmx_simd_real_t__m128 bzx; |
124 | gmx_simd_real_t__m128 bzy; |
125 | gmx_simd_real_t__m128 bzz; |
126 | gmx_simd_real_t__m128 byx; |
127 | gmx_simd_real_t__m128 byy; |
128 | gmx_simd_real_t__m128 bxx; |
129 | } pbc_simd_t; |
130 | |
131 | /* Set the SIMD pbc data from a normal t_pbc struct */ |
132 | static void set_pbc_simd(const t_pbc *pbc, pbc_simd_t *pbc_simd) |
133 | { |
134 | rvec inv_bdiag; |
135 | int d; |
136 | |
137 | /* Setting inv_bdiag to 0 effectively turns off PBC */ |
138 | clear_rvec(inv_bdiag); |
139 | if (pbc != NULL((void*)0)) |
140 | { |
141 | for (d = 0; d < pbc->ndim_ePBC; d++) |
142 | { |
143 | inv_bdiag[d] = 1.0/pbc->box[d][d]; |
144 | } |
145 | } |
146 | |
147 | pbc_simd->inv_bzz = gmx_simd_set1_r_mm_set1_ps(inv_bdiag[ZZ2]); |
148 | pbc_simd->inv_byy = gmx_simd_set1_r_mm_set1_ps(inv_bdiag[YY1]); |
149 | pbc_simd->inv_bxx = gmx_simd_set1_r_mm_set1_ps(inv_bdiag[XX0]); |
150 | |
151 | if (pbc != NULL((void*)0)) |
152 | { |
153 | pbc_simd->bzx = gmx_simd_set1_r_mm_set1_ps(pbc->box[ZZ2][XX0]); |
154 | pbc_simd->bzy = gmx_simd_set1_r_mm_set1_ps(pbc->box[ZZ2][YY1]); |
155 | pbc_simd->bzz = gmx_simd_set1_r_mm_set1_ps(pbc->box[ZZ2][ZZ2]); |
156 | pbc_simd->byx = gmx_simd_set1_r_mm_set1_ps(pbc->box[YY1][XX0]); |
157 | pbc_simd->byy = gmx_simd_set1_r_mm_set1_ps(pbc->box[YY1][YY1]); |
158 | pbc_simd->bxx = gmx_simd_set1_r_mm_set1_ps(pbc->box[XX0][XX0]); |
159 | } |
160 | else |
161 | { |
162 | pbc_simd->bzx = gmx_simd_setzero_r_mm_setzero_ps(); |
163 | pbc_simd->bzy = gmx_simd_setzero_r_mm_setzero_ps(); |
164 | pbc_simd->bzz = gmx_simd_setzero_r_mm_setzero_ps(); |
165 | pbc_simd->byx = gmx_simd_setzero_r_mm_setzero_ps(); |
166 | pbc_simd->byy = gmx_simd_setzero_r_mm_setzero_ps(); |
167 | pbc_simd->bxx = gmx_simd_setzero_r_mm_setzero_ps(); |
168 | } |
169 | } |
170 | |
171 | /* Correct distance vector *dx,*dy,*dz for PBC using SIMD */ |
172 | static gmx_inlineinline void |
173 | pbc_dx_simd(gmx_simd_real_t__m128 *dx, gmx_simd_real_t__m128 *dy, gmx_simd_real_t__m128 *dz, |
174 | const pbc_simd_t *pbc) |
175 | { |
176 | gmx_simd_real_t__m128 sh; |
177 | |
178 | sh = gmx_simd_round_r(gmx_simd_mul_r(*dz, pbc->inv_bzz))__extension__ ({ __m128 __X = (_mm_mul_ps(*dz, pbc->inv_bzz )); (__m128) __builtin_ia32_roundps((__v4sf)__X, ((0x00 | 0x00 ))); }); |
179 | *dx = gmx_simd_fnmadd_r(sh, pbc->bzx, *dx)_mm_sub_ps(*dx, _mm_mul_ps(sh, pbc->bzx)); |
180 | *dy = gmx_simd_fnmadd_r(sh, pbc->bzy, *dy)_mm_sub_ps(*dy, _mm_mul_ps(sh, pbc->bzy)); |
181 | *dz = gmx_simd_fnmadd_r(sh, pbc->bzz, *dz)_mm_sub_ps(*dz, _mm_mul_ps(sh, pbc->bzz)); |
182 | |
183 | sh = gmx_simd_round_r(gmx_simd_mul_r(*dy, pbc->inv_byy))__extension__ ({ __m128 __X = (_mm_mul_ps(*dy, pbc->inv_byy )); (__m128) __builtin_ia32_roundps((__v4sf)__X, ((0x00 | 0x00 ))); }); |
184 | *dx = gmx_simd_fnmadd_r(sh, pbc->byx, *dx)_mm_sub_ps(*dx, _mm_mul_ps(sh, pbc->byx)); |
185 | *dy = gmx_simd_fnmadd_r(sh, pbc->byy, *dy)_mm_sub_ps(*dy, _mm_mul_ps(sh, pbc->byy)); |
186 | |
187 | sh = gmx_simd_round_r(gmx_simd_mul_r(*dx, pbc->inv_bxx))__extension__ ({ __m128 __X = (_mm_mul_ps(*dx, pbc->inv_bxx )); (__m128) __builtin_ia32_roundps((__v4sf)__X, ((0x00 | 0x00 ))); }); |
188 | *dx = gmx_simd_fnmadd_r(sh, pbc->bxx, *dx)_mm_sub_ps(*dx, _mm_mul_ps(sh, pbc->bxx)); |
189 | } |
190 | |
191 | #endif /* GMX_SIMD_HAVE_REAL */ |
192 | |
193 | /* |
194 | * Morse potential bond by Frank Everdij |
195 | * |
196 | * Three parameters needed: |
197 | * |
198 | * b0 = equilibrium distance in nm |
199 | * be = beta in nm^-1 (actually, it's nu_e*Sqrt(2*pi*pi*mu/D_e)) |
200 | * cb = well depth in kJ/mol |
201 | * |
202 | * Note: the potential is referenced to be +cb at infinite separation |
203 | * and zero at the equilibrium distance! |
204 | */ |
205 | |
206 | real morse_bonds(int nbonds, |
207 | const t_iatom forceatoms[], const t_iparams forceparams[], |
208 | const rvec x[], rvec f[], rvec fshift[], |
209 | const t_pbc *pbc, const t_graph *g, |
210 | real lambda, real *dvdlambda, |
211 | const t_mdatoms gmx_unused__attribute__ ((unused)) *md, t_fcdata gmx_unused__attribute__ ((unused)) *fcd, |
212 | int gmx_unused__attribute__ ((unused)) *global_atom_index) |
213 | { |
214 | const real one = 1.0; |
215 | const real two = 2.0; |
216 | real dr, dr2, temp, omtemp, cbomtemp, fbond, vbond, fij, vtot; |
217 | real b0, be, cb, b0A, beA, cbA, b0B, beB, cbB, L1; |
218 | rvec dx; |
219 | int i, m, ki, type, ai, aj; |
220 | ivec dt; |
221 | |
222 | vtot = 0.0; |
223 | for (i = 0; (i < nbonds); ) |
224 | { |
225 | type = forceatoms[i++]; |
226 | ai = forceatoms[i++]; |
227 | aj = forceatoms[i++]; |
228 | |
229 | b0A = forceparams[type].morse.b0A; |
230 | beA = forceparams[type].morse.betaA; |
231 | cbA = forceparams[type].morse.cbA; |
232 | |
233 | b0B = forceparams[type].morse.b0B; |
234 | beB = forceparams[type].morse.betaB; |
235 | cbB = forceparams[type].morse.cbB; |
236 | |
237 | L1 = one-lambda; /* 1 */ |
238 | b0 = L1*b0A + lambda*b0B; /* 3 */ |
239 | be = L1*beA + lambda*beB; /* 3 */ |
240 | cb = L1*cbA + lambda*cbB; /* 3 */ |
241 | |
242 | ki = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /* 3 */ |
243 | dr2 = iprod(dx, dx); /* 5 */ |
244 | dr = dr2*gmx_invsqrt(dr2)gmx_software_invsqrt(dr2); /* 10 */ |
245 | temp = exp(-be*(dr-b0)); /* 12 */ |
246 | |
247 | if (temp == one) |
248 | { |
249 | /* bonds are constrainted. This may _not_ include bond constraints if they are lambda dependent */ |
250 | *dvdlambda += cbB-cbA; |
251 | continue; |
252 | } |
253 | |
254 | omtemp = one-temp; /* 1 */ |
255 | cbomtemp = cb*omtemp; /* 1 */ |
256 | vbond = cbomtemp*omtemp; /* 1 */ |
257 | fbond = -two*be*temp*cbomtemp*gmx_invsqrt(dr2)gmx_software_invsqrt(dr2); /* 9 */ |
258 | vtot += vbond; /* 1 */ |
259 | |
260 | *dvdlambda += (cbB - cbA) * omtemp * omtemp - (2-2*omtemp)*omtemp * cb * ((b0B-b0A)*be - (beB-beA)*(dr-b0)); /* 15 */ |
261 | |
262 | if (g) |
263 | { |
264 | ivec_sub(SHIFT_IVEC(g, ai)((g)->ishift[ai]), SHIFT_IVEC(g, aj)((g)->ishift[aj]), dt); |
265 | ki = IVEC2IS(dt)(((2*2 +1)*((2*1 +1)*(((dt)[2])+1)+((dt)[1])+1)+((dt)[0])+2)); |
266 | } |
267 | |
268 | for (m = 0; (m < DIM3); m++) /* 15 */ |
269 | { |
270 | fij = fbond*dx[m]; |
271 | f[ai][m] += fij; |
272 | f[aj][m] -= fij; |
273 | fshift[ki][m] += fij; |
274 | fshift[CENTRAL(((2*1 +1)*(2*1 +1)*(2*2 +1))/2)][m] -= fij; |
275 | } |
276 | } /* 83 TOTAL */ |
277 | return vtot; |
278 | } |
279 | |
280 | real cubic_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 gmx_unused__attribute__ ((unused)) lambda, real gmx_unused__attribute__ ((unused)) *dvdlambda, |
285 | const t_mdatoms gmx_unused__attribute__ ((unused)) *md, t_fcdata gmx_unused__attribute__ ((unused)) *fcd, |
286 | int gmx_unused__attribute__ ((unused)) *global_atom_index) |
287 | { |
288 | const real three = 3.0; |
289 | const real two = 2.0; |
290 | real kb, b0, kcub; |
291 | real dr, dr2, dist, kdist, kdist2, fbond, vbond, fij, vtot; |
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 | b0 = forceparams[type].cubic.b0; |
304 | kb = forceparams[type].cubic.kb; |
305 | kcub = forceparams[type].cubic.kcub; |
306 | |
307 | ki = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /* 3 */ |
308 | dr2 = iprod(dx, dx); /* 5 */ |
309 | |
310 | if (dr2 == 0.0) |
311 | { |
312 | continue; |
313 | } |
314 | |
315 | dr = dr2*gmx_invsqrt(dr2)gmx_software_invsqrt(dr2); /* 10 */ |
316 | dist = dr-b0; |
317 | kdist = kb*dist; |
318 | kdist2 = kdist*dist; |
319 | |
320 | vbond = kdist2 + kcub*kdist2*dist; |
321 | fbond = -(two*kdist + three*kdist2*kcub)/dr; |
322 | |
323 | vtot += vbond; /* 21 */ |
324 | |
325 | if (g) |
326 | { |
327 | ivec_sub(SHIFT_IVEC(g, ai)((g)->ishift[ai]), SHIFT_IVEC(g, aj)((g)->ishift[aj]), dt); |
328 | ki = IVEC2IS(dt)(((2*2 +1)*((2*1 +1)*(((dt)[2])+1)+((dt)[1])+1)+((dt)[0])+2)); |
329 | } |
330 | for (m = 0; (m < DIM3); m++) /* 15 */ |
331 | { |
332 | fij = fbond*dx[m]; |
333 | f[ai][m] += fij; |
334 | f[aj][m] -= fij; |
335 | fshift[ki][m] += fij; |
336 | fshift[CENTRAL(((2*1 +1)*(2*1 +1)*(2*2 +1))/2)][m] -= fij; |
337 | } |
338 | } /* 54 TOTAL */ |
339 | return vtot; |
340 | } |
341 | |
342 | real FENE_bonds(int nbonds, |
343 | const t_iatom forceatoms[], const t_iparams forceparams[], |
344 | const rvec x[], rvec f[], rvec fshift[], |
345 | const t_pbc *pbc, const t_graph *g, |
346 | real gmx_unused__attribute__ ((unused)) lambda, real gmx_unused__attribute__ ((unused)) *dvdlambda, |
347 | const t_mdatoms gmx_unused__attribute__ ((unused)) *md, t_fcdata gmx_unused__attribute__ ((unused)) *fcd, |
348 | int *global_atom_index) |
349 | { |
350 | const real half = 0.5; |
351 | const real one = 1.0; |
352 | real bm, kb; |
353 | real dr, dr2, bm2, omdr2obm2, fbond, vbond, fij, vtot; |
354 | rvec dx; |
355 | int i, m, ki, type, ai, aj; |
356 | ivec dt; |
357 | |
358 | vtot = 0.0; |
359 | for (i = 0; (i < nbonds); ) |
360 | { |
361 | type = forceatoms[i++]; |
362 | ai = forceatoms[i++]; |
363 | aj = forceatoms[i++]; |
364 | |
365 | bm = forceparams[type].fene.bm; |
366 | kb = forceparams[type].fene.kb; |
367 | |
368 | ki = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /* 3 */ |
369 | dr2 = iprod(dx, dx); /* 5 */ |
370 | |
371 | if (dr2 == 0.0) |
372 | { |
373 | continue; |
374 | } |
375 | |
376 | bm2 = bm*bm; |
377 | |
378 | if (dr2 >= bm2) |
379 | { |
380 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxlib/bondfree.c" , 380, |
381 | "r^2 (%f) >= bm^2 (%f) in FENE bond between atoms %d and %d", |
382 | dr2, bm2, |
383 | glatnr(global_atom_index, ai), |
384 | glatnr(global_atom_index, aj)); |
385 | } |
386 | |
387 | omdr2obm2 = one - dr2/bm2; |
388 | |
389 | vbond = -half*kb*bm2*log(omdr2obm2); |
390 | fbond = -kb/omdr2obm2; |
391 | |
392 | vtot += vbond; /* 35 */ |
393 | |
394 | if (g) |
395 | { |
396 | ivec_sub(SHIFT_IVEC(g, ai)((g)->ishift[ai]), SHIFT_IVEC(g, aj)((g)->ishift[aj]), dt); |
397 | ki = IVEC2IS(dt)(((2*2 +1)*((2*1 +1)*(((dt)[2])+1)+((dt)[1])+1)+((dt)[0])+2)); |
398 | } |
399 | for (m = 0; (m < DIM3); m++) /* 15 */ |
400 | { |
401 | fij = fbond*dx[m]; |
402 | f[ai][m] += fij; |
403 | f[aj][m] -= fij; |
404 | fshift[ki][m] += fij; |
405 | fshift[CENTRAL(((2*1 +1)*(2*1 +1)*(2*2 +1))/2)][m] -= fij; |
406 | } |
407 | } /* 58 TOTAL */ |
408 | return vtot; |
409 | } |
410 | |
411 | real harmonic(real kA, real kB, real xA, real xB, real x, real lambda, |
412 | real *V, real *F) |
413 | { |
414 | const real half = 0.5; |
415 | real L1, kk, x0, dx, dx2; |
416 | real v, f, dvdlambda; |
417 | |
418 | L1 = 1.0-lambda; |
419 | kk = L1*kA+lambda*kB; |
420 | x0 = L1*xA+lambda*xB; |
421 | |
422 | dx = x-x0; |
423 | dx2 = dx*dx; |
424 | |
425 | f = -kk*dx; |
426 | v = half*kk*dx2; |
427 | dvdlambda = half*(kB-kA)*dx2 + (xA-xB)*kk*dx; |
428 | |
429 | *F = f; |
430 | *V = v; |
431 | |
432 | return dvdlambda; |
433 | |
434 | /* That was 19 flops */ |
435 | } |
436 | |
437 | |
438 | real bonds(int nbonds, |
439 | const t_iatom forceatoms[], const t_iparams forceparams[], |
440 | const rvec x[], rvec f[], rvec fshift[], |
441 | const t_pbc *pbc, const t_graph *g, |
442 | real lambda, real *dvdlambda, |
443 | const t_mdatoms gmx_unused__attribute__ ((unused)) *md, t_fcdata gmx_unused__attribute__ ((unused)) *fcd, |
444 | int gmx_unused__attribute__ ((unused)) *global_atom_index) |
445 | { |
446 | int i, m, ki, ai, aj, type; |
447 | real dr, dr2, fbond, vbond, fij, vtot; |
448 | rvec dx; |
449 | ivec dt; |
450 | |
451 | vtot = 0.0; |
452 | for (i = 0; (i < nbonds); ) |
453 | { |
454 | type = forceatoms[i++]; |
455 | ai = forceatoms[i++]; |
456 | aj = forceatoms[i++]; |
457 | |
458 | ki = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /* 3 */ |
459 | dr2 = iprod(dx, dx); /* 5 */ |
460 | dr = dr2*gmx_invsqrt(dr2)gmx_software_invsqrt(dr2); /* 10 */ |
461 | |
462 | *dvdlambda += harmonic(forceparams[type].harmonic.krA, |
463 | forceparams[type].harmonic.krB, |
464 | forceparams[type].harmonic.rA, |
465 | forceparams[type].harmonic.rB, |
466 | dr, lambda, &vbond, &fbond); /* 19 */ |
467 | |
468 | if (dr2 == 0.0) |
469 | { |
470 | continue; |
471 | } |
472 | |
473 | |
474 | vtot += vbond; /* 1*/ |
475 | fbond *= gmx_invsqrt(dr2)gmx_software_invsqrt(dr2); /* 6 */ |
476 | #ifdef DEBUG |
477 | if (debug) |
478 | { |
479 | fprintf(debug, "BONDS: dr = %10g vbond = %10g fbond = %10g\n", |
480 | dr, vbond, fbond); |
481 | } |
482 | #endif |
483 | if (g) |
484 | { |
485 | ivec_sub(SHIFT_IVEC(g, ai)((g)->ishift[ai]), SHIFT_IVEC(g, aj)((g)->ishift[aj]), dt); |
486 | ki = IVEC2IS(dt)(((2*2 +1)*((2*1 +1)*(((dt)[2])+1)+((dt)[1])+1)+((dt)[0])+2)); |
487 | } |
488 | for (m = 0; (m < DIM3); m++) /* 15 */ |
489 | { |
490 | fij = fbond*dx[m]; |
491 | f[ai][m] += fij; |
492 | f[aj][m] -= fij; |
493 | fshift[ki][m] += fij; |
494 | fshift[CENTRAL(((2*1 +1)*(2*1 +1)*(2*2 +1))/2)][m] -= fij; |
495 | } |
496 | } /* 59 TOTAL */ |
497 | return vtot; |
498 | } |
499 | |
500 | real restraint_bonds(int nbonds, |
501 | const t_iatom forceatoms[], const t_iparams forceparams[], |
502 | const rvec x[], rvec f[], rvec fshift[], |
503 | const t_pbc *pbc, const t_graph *g, |
504 | real lambda, real *dvdlambda, |
505 | const t_mdatoms gmx_unused__attribute__ ((unused)) *md, t_fcdata gmx_unused__attribute__ ((unused)) *fcd, |
506 | int gmx_unused__attribute__ ((unused)) *global_atom_index) |
507 | { |
508 | int i, m, ki, ai, aj, type; |
509 | real dr, dr2, fbond, vbond, fij, vtot; |
510 | real L1; |
511 | real low, dlow, up1, dup1, up2, dup2, k, dk; |
512 | real drh, drh2; |
513 | rvec dx; |
514 | ivec dt; |
515 | |
516 | L1 = 1.0 - lambda; |
517 | |
518 | vtot = 0.0; |
519 | for (i = 0; (i < nbonds); ) |
520 | { |
521 | type = forceatoms[i++]; |
522 | ai = forceatoms[i++]; |
523 | aj = forceatoms[i++]; |
524 | |
525 | ki = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /* 3 */ |
526 | dr2 = iprod(dx, dx); /* 5 */ |
527 | dr = dr2*gmx_invsqrt(dr2)gmx_software_invsqrt(dr2); /* 10 */ |
528 | |
529 | low = L1*forceparams[type].restraint.lowA + lambda*forceparams[type].restraint.lowB; |
530 | dlow = -forceparams[type].restraint.lowA + forceparams[type].restraint.lowB; |
531 | up1 = L1*forceparams[type].restraint.up1A + lambda*forceparams[type].restraint.up1B; |
532 | dup1 = -forceparams[type].restraint.up1A + forceparams[type].restraint.up1B; |
533 | up2 = L1*forceparams[type].restraint.up2A + lambda*forceparams[type].restraint.up2B; |
534 | dup2 = -forceparams[type].restraint.up2A + forceparams[type].restraint.up2B; |
535 | k = L1*forceparams[type].restraint.kA + lambda*forceparams[type].restraint.kB; |
536 | dk = -forceparams[type].restraint.kA + forceparams[type].restraint.kB; |
537 | /* 24 */ |
538 | |
539 | if (dr < low) |
540 | { |
541 | drh = dr - low; |
542 | drh2 = drh*drh; |
543 | vbond = 0.5*k*drh2; |
544 | fbond = -k*drh; |
545 | *dvdlambda += 0.5*dk*drh2 - k*dlow*drh; |
546 | } /* 11 */ |
547 | else if (dr <= up1) |
548 | { |
549 | vbond = 0; |
550 | fbond = 0; |
551 | } |
552 | else if (dr <= up2) |
553 | { |
554 | drh = dr - up1; |
555 | drh2 = drh*drh; |
556 | vbond = 0.5*k*drh2; |
557 | fbond = -k*drh; |
558 | *dvdlambda += 0.5*dk*drh2 - k*dup1*drh; |
559 | } /* 11 */ |
560 | else |
561 | { |
562 | drh = dr - up2; |
563 | vbond = k*(up2 - up1)*(0.5*(up2 - up1) + drh); |
564 | fbond = -k*(up2 - up1); |
565 | *dvdlambda += dk*(up2 - up1)*(0.5*(up2 - up1) + drh) |
566 | + k*(dup2 - dup1)*(up2 - up1 + drh) |
567 | - k*(up2 - up1)*dup2; |
568 | } |
569 | |
570 | if (dr2 == 0.0) |
571 | { |
572 | continue; |
573 | } |
574 | |
575 | vtot += vbond; /* 1*/ |
576 | fbond *= gmx_invsqrt(dr2)gmx_software_invsqrt(dr2); /* 6 */ |
577 | #ifdef DEBUG |
578 | if (debug) |
579 | { |
580 | fprintf(debug, "BONDS: dr = %10g vbond = %10g fbond = %10g\n", |
581 | dr, vbond, fbond); |
582 | } |
583 | #endif |
584 | if (g) |
585 | { |
586 | ivec_sub(SHIFT_IVEC(g, ai)((g)->ishift[ai]), SHIFT_IVEC(g, aj)((g)->ishift[aj]), dt); |
587 | ki = IVEC2IS(dt)(((2*2 +1)*((2*1 +1)*(((dt)[2])+1)+((dt)[1])+1)+((dt)[0])+2)); |
588 | } |
589 | for (m = 0; (m < DIM3); m++) /* 15 */ |
590 | { |
591 | fij = fbond*dx[m]; |
592 | f[ai][m] += fij; |
593 | f[aj][m] -= fij; |
594 | fshift[ki][m] += fij; |
595 | fshift[CENTRAL(((2*1 +1)*(2*1 +1)*(2*2 +1))/2)][m] -= fij; |
596 | } |
597 | } /* 59 TOTAL */ |
598 | |
599 | return vtot; |
600 | } |
601 | |
602 | real polarize(int nbonds, |
603 | const t_iatom forceatoms[], const t_iparams forceparams[], |
604 | const rvec x[], rvec f[], rvec fshift[], |
605 | const t_pbc *pbc, const t_graph *g, |
606 | real lambda, real *dvdlambda, |
607 | const t_mdatoms *md, t_fcdata gmx_unused__attribute__ ((unused)) *fcd, |
608 | int gmx_unused__attribute__ ((unused)) *global_atom_index) |
609 | { |
610 | int i, m, ki, ai, aj, type; |
611 | real dr, dr2, fbond, vbond, fij, vtot, ksh; |
612 | rvec dx; |
613 | ivec dt; |
614 | |
615 | vtot = 0.0; |
616 | for (i = 0; (i < nbonds); ) |
617 | { |
618 | type = forceatoms[i++]; |
619 | ai = forceatoms[i++]; |
620 | aj = forceatoms[i++]; |
621 | ksh = sqr(md->chargeA[aj])*ONE_4PI_EPS0((332.0636930*(4.184))*0.1)/forceparams[type].polarize.alpha; |
622 | if (debug) |
623 | { |
624 | fprintf(debug, "POL: local ai = %d aj = %d ksh = %.3f\n", ai, aj, ksh); |
625 | } |
626 | |
627 | ki = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /* 3 */ |
628 | dr2 = iprod(dx, dx); /* 5 */ |
629 | dr = dr2*gmx_invsqrt(dr2)gmx_software_invsqrt(dr2); /* 10 */ |
630 | |
631 | *dvdlambda += harmonic(ksh, ksh, 0, 0, dr, lambda, &vbond, &fbond); /* 19 */ |
632 | |
633 | if (dr2 == 0.0) |
634 | { |
635 | continue; |
636 | } |
637 | |
638 | vtot += vbond; /* 1*/ |
639 | fbond *= gmx_invsqrt(dr2)gmx_software_invsqrt(dr2); /* 6 */ |
640 | |
641 | if (g) |
642 | { |
643 | ivec_sub(SHIFT_IVEC(g, ai)((g)->ishift[ai]), SHIFT_IVEC(g, aj)((g)->ishift[aj]), dt); |
644 | ki = IVEC2IS(dt)(((2*2 +1)*((2*1 +1)*(((dt)[2])+1)+((dt)[1])+1)+((dt)[0])+2)); |
645 | } |
646 | for (m = 0; (m < DIM3); m++) /* 15 */ |
647 | { |
648 | fij = fbond*dx[m]; |
649 | f[ai][m] += fij; |
650 | f[aj][m] -= fij; |
651 | fshift[ki][m] += fij; |
652 | fshift[CENTRAL(((2*1 +1)*(2*1 +1)*(2*2 +1))/2)][m] -= fij; |
653 | } |
654 | } /* 59 TOTAL */ |
655 | return vtot; |
656 | } |
657 | |
658 | real anharm_polarize(int nbonds, |
659 | const t_iatom forceatoms[], const t_iparams forceparams[], |
660 | const rvec x[], rvec f[], rvec fshift[], |
661 | const t_pbc *pbc, const t_graph *g, |
662 | real lambda, real *dvdlambda, |
663 | const t_mdatoms *md, t_fcdata gmx_unused__attribute__ ((unused)) *fcd, |
664 | int gmx_unused__attribute__ ((unused)) *global_atom_index) |
665 | { |
666 | int i, m, ki, ai, aj, type; |
667 | real dr, dr2, fbond, vbond, fij, vtot, ksh, khyp, drcut, ddr, ddr3; |
668 | rvec dx; |
669 | ivec dt; |
670 | |
671 | vtot = 0.0; |
672 | for (i = 0; (i < nbonds); ) |
673 | { |
674 | type = forceatoms[i++]; |
675 | ai = forceatoms[i++]; |
676 | aj = forceatoms[i++]; |
677 | ksh = sqr(md->chargeA[aj])*ONE_4PI_EPS0((332.0636930*(4.184))*0.1)/forceparams[type].anharm_polarize.alpha; /* 7*/ |
678 | khyp = forceparams[type].anharm_polarize.khyp; |
679 | drcut = forceparams[type].anharm_polarize.drcut; |
680 | if (debug) |
681 | { |
682 | fprintf(debug, "POL: local ai = %d aj = %d ksh = %.3f\n", ai, aj, ksh); |
683 | } |
684 | |
685 | ki = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /* 3 */ |
686 | dr2 = iprod(dx, dx); /* 5 */ |
687 | dr = dr2*gmx_invsqrt(dr2)gmx_software_invsqrt(dr2); /* 10 */ |
688 | |
689 | *dvdlambda += harmonic(ksh, ksh, 0, 0, dr, lambda, &vbond, &fbond); /* 19 */ |
690 | |
691 | if (dr2 == 0.0) |
692 | { |
693 | continue; |
694 | } |
695 | |
696 | if (dr > drcut) |
697 | { |
698 | ddr = dr-drcut; |
699 | ddr3 = ddr*ddr*ddr; |
700 | vbond += khyp*ddr*ddr3; |
701 | fbond -= 4*khyp*ddr3; |
702 | } |
703 | fbond *= gmx_invsqrt(dr2)gmx_software_invsqrt(dr2); /* 6 */ |
704 | vtot += vbond; /* 1*/ |
705 | |
706 | if (g) |
707 | { |
708 | ivec_sub(SHIFT_IVEC(g, ai)((g)->ishift[ai]), SHIFT_IVEC(g, aj)((g)->ishift[aj]), dt); |
709 | ki = IVEC2IS(dt)(((2*2 +1)*((2*1 +1)*(((dt)[2])+1)+((dt)[1])+1)+((dt)[0])+2)); |
710 | } |
711 | for (m = 0; (m < DIM3); m++) /* 15 */ |
712 | { |
713 | fij = fbond*dx[m]; |
714 | f[ai][m] += fij; |
715 | f[aj][m] -= fij; |
716 | fshift[ki][m] += fij; |
717 | fshift[CENTRAL(((2*1 +1)*(2*1 +1)*(2*2 +1))/2)][m] -= fij; |
718 | } |
719 | } /* 72 TOTAL */ |
720 | return vtot; |
721 | } |
722 | |
723 | real water_pol(int nbonds, |
724 | const t_iatom forceatoms[], const t_iparams forceparams[], |
725 | const rvec x[], rvec f[], rvec gmx_unused__attribute__ ((unused)) fshift[], |
726 | const t_pbc gmx_unused__attribute__ ((unused)) *pbc, const t_graph gmx_unused__attribute__ ((unused)) *g, |
727 | real gmx_unused__attribute__ ((unused)) lambda, real gmx_unused__attribute__ ((unused)) *dvdlambda, |
728 | const t_mdatoms gmx_unused__attribute__ ((unused)) *md, t_fcdata gmx_unused__attribute__ ((unused)) *fcd, |
729 | int gmx_unused__attribute__ ((unused)) *global_atom_index) |
730 | { |
731 | /* This routine implements anisotropic polarizibility for water, through |
732 | * a shell connected to a dummy with spring constant that differ in the |
733 | * three spatial dimensions in the molecular frame. |
734 | */ |
735 | int i, m, aO, aH1, aH2, aD, aS, type, type0; |
736 | rvec dOH1, dOH2, dHH, dOD, dDS, nW, kk, dx, kdx, proj; |
737 | #ifdef DEBUG |
738 | rvec df; |
739 | #endif |
740 | real vtot, fij, r_HH, r_OD, r_nW, tx, ty, tz, qS; |
741 | |
742 | vtot = 0.0; |
743 | if (nbonds > 0) |
744 | { |
745 | type0 = forceatoms[0]; |
746 | aS = forceatoms[5]; |
747 | qS = md->chargeA[aS]; |
748 | kk[XX0] = sqr(qS)*ONE_4PI_EPS0((332.0636930*(4.184))*0.1)/forceparams[type0].wpol.al_x; |
749 | kk[YY1] = sqr(qS)*ONE_4PI_EPS0((332.0636930*(4.184))*0.1)/forceparams[type0].wpol.al_y; |
750 | kk[ZZ2] = sqr(qS)*ONE_4PI_EPS0((332.0636930*(4.184))*0.1)/forceparams[type0].wpol.al_z; |
751 | r_HH = 1.0/forceparams[type0].wpol.rHH; |
752 | r_OD = 1.0/forceparams[type0].wpol.rOD; |
753 | if (debug) |
754 | { |
755 | fprintf(debug, "WPOL: qS = %10.5f aS = %5d\n", qS, aS); |
756 | fprintf(debug, "WPOL: kk = %10.3f %10.3f %10.3f\n", |
757 | kk[XX0], kk[YY1], kk[ZZ2]); |
758 | fprintf(debug, "WPOL: rOH = %10.3f rHH = %10.3f rOD = %10.3f\n", |
759 | forceparams[type0].wpol.rOH, |
760 | forceparams[type0].wpol.rHH, |
761 | forceparams[type0].wpol.rOD); |
762 | } |
763 | for (i = 0; (i < nbonds); i += 6) |
764 | { |
765 | type = forceatoms[i]; |
766 | if (type != type0) |
767 | { |
768 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxlib/bondfree.c" , 768, "Sorry, type = %d, type0 = %d, file = %s, line = %d", |
769 | type, type0, __FILE__"/home/alexxy/Develop/gromacs/src/gromacs/gmxlib/bondfree.c", __LINE__769); |
770 | } |
771 | aO = forceatoms[i+1]; |
772 | aH1 = forceatoms[i+2]; |
773 | aH2 = forceatoms[i+3]; |
774 | aD = forceatoms[i+4]; |
775 | aS = forceatoms[i+5]; |
776 | |
777 | /* Compute vectors describing the water frame */ |
778 | rvec_sub(x[aH1], x[aO], dOH1); |
779 | rvec_sub(x[aH2], x[aO], dOH2); |
780 | rvec_sub(x[aH2], x[aH1], dHH); |
781 | rvec_sub(x[aD], x[aO], dOD); |
782 | rvec_sub(x[aS], x[aD], dDS); |
783 | cprod(dOH1, dOH2, nW); |
784 | |
785 | /* Compute inverse length of normal vector |
786 | * (this one could be precomputed, but I'm too lazy now) |
787 | */ |
788 | r_nW = gmx_invsqrt(iprod(nW, nW))gmx_software_invsqrt(iprod(nW, nW)); |
789 | /* This is for precision, but does not make a big difference, |
790 | * it can go later. |
791 | */ |
792 | r_OD = gmx_invsqrt(iprod(dOD, dOD))gmx_software_invsqrt(iprod(dOD, dOD)); |
793 | |
794 | /* Normalize the vectors in the water frame */ |
795 | svmul(r_nW, nW, nW); |
796 | svmul(r_HH, dHH, dHH); |
797 | svmul(r_OD, dOD, dOD); |
798 | |
799 | /* Compute displacement of shell along components of the vector */ |
800 | dx[ZZ2] = iprod(dDS, dOD); |
801 | /* Compute projection on the XY plane: dDS - dx[ZZ]*dOD */ |
802 | for (m = 0; (m < DIM3); m++) |
803 | { |
804 | proj[m] = dDS[m]-dx[ZZ2]*dOD[m]; |
805 | } |
806 | |
807 | /*dx[XX] = iprod(dDS,nW); |
808 | dx[YY] = iprod(dDS,dHH);*/ |
809 | dx[XX0] = iprod(proj, nW); |
810 | for (m = 0; (m < DIM3); m++) |
811 | { |
812 | proj[m] -= dx[XX0]*nW[m]; |
813 | } |
814 | dx[YY1] = iprod(proj, dHH); |
815 | /*#define DEBUG*/ |
816 | #ifdef DEBUG |
817 | if (debug) |
818 | { |
819 | fprintf(debug, "WPOL: dx2=%10g dy2=%10g dz2=%10g sum=%10g dDS^2=%10g\n", |
820 | sqr(dx[XX0]), sqr(dx[YY1]), sqr(dx[ZZ2]), iprod(dx, dx), iprod(dDS, dDS)); |
821 | fprintf(debug, "WPOL: dHH=(%10g,%10g,%10g)\n", dHH[XX0], dHH[YY1], dHH[ZZ2]); |
822 | fprintf(debug, "WPOL: dOD=(%10g,%10g,%10g), 1/r_OD = %10g\n", |
823 | dOD[XX0], dOD[YY1], dOD[ZZ2], 1/r_OD); |
824 | fprintf(debug, "WPOL: nW =(%10g,%10g,%10g), 1/r_nW = %10g\n", |
825 | nW[XX0], nW[YY1], nW[ZZ2], 1/r_nW); |
826 | fprintf(debug, "WPOL: dx =%10g, dy =%10g, dz =%10g\n", |
827 | dx[XX0], dx[YY1], dx[ZZ2]); |
828 | fprintf(debug, "WPOL: dDSx=%10g, dDSy=%10g, dDSz=%10g\n", |
829 | dDS[XX0], dDS[YY1], dDS[ZZ2]); |
830 | } |
831 | #endif |
832 | /* Now compute the forces and energy */ |
833 | kdx[XX0] = kk[XX0]*dx[XX0]; |
834 | kdx[YY1] = kk[YY1]*dx[YY1]; |
835 | kdx[ZZ2] = kk[ZZ2]*dx[ZZ2]; |
836 | vtot += iprod(dx, kdx); |
837 | for (m = 0; (m < DIM3); m++) |
838 | { |
839 | /* This is a tensor operation but written out for speed */ |
840 | tx = nW[m]*kdx[XX0]; |
841 | ty = dHH[m]*kdx[YY1]; |
842 | tz = dOD[m]*kdx[ZZ2]; |
843 | fij = -tx-ty-tz; |
844 | #ifdef DEBUG |
845 | df[m] = fij; |
846 | #endif |
847 | f[aS][m] += fij; |
848 | f[aD][m] -= fij; |
849 | } |
850 | #ifdef DEBUG |
851 | if (debug) |
852 | { |
853 | fprintf(debug, "WPOL: vwpol=%g\n", 0.5*iprod(dx, kdx)); |
854 | fprintf(debug, "WPOL: df = (%10g, %10g, %10g)\n", df[XX0], df[YY1], df[ZZ2]); |
855 | } |
856 | #endif |
857 | } |
858 | } |
859 | return 0.5*vtot; |
860 | } |
861 | |
862 | static real do_1_thole(const rvec xi, const rvec xj, rvec fi, rvec fj, |
863 | const t_pbc *pbc, real qq, |
864 | rvec fshift[], real afac) |
865 | { |
866 | rvec r12; |
867 | real r12sq, r12_1, r12n, r12bar, v0, v1, fscal, ebar, fff; |
868 | int m, t; |
869 | |
870 | t = pbc_rvec_sub(pbc, xi, xj, r12); /* 3 */ |
871 | |
872 | r12sq = iprod(r12, r12); /* 5 */ |
873 | r12_1 = gmx_invsqrt(r12sq)gmx_software_invsqrt(r12sq); /* 5 */ |
874 | r12bar = afac/r12_1; /* 5 */ |
875 | v0 = qq*ONE_4PI_EPS0((332.0636930*(4.184))*0.1)*r12_1; /* 2 */ |
876 | ebar = exp(-r12bar); /* 5 */ |
877 | v1 = (1-(1+0.5*r12bar)*ebar); /* 4 */ |
878 | fscal = ((v0*r12_1)*v1 - v0*0.5*afac*ebar*(r12bar+1))*r12_1; /* 9 */ |
879 | if (debug) |
880 | { |
881 | fprintf(debug, "THOLE: v0 = %.3f v1 = %.3f r12= % .3f r12bar = %.3f fscal = %.3f ebar = %.3f\n", v0, v1, 1/r12_1, r12bar, fscal, ebar); |
882 | } |
883 | |
884 | for (m = 0; (m < DIM3); m++) |
885 | { |
886 | fff = fscal*r12[m]; |
887 | fi[m] += fff; |
888 | fj[m] -= fff; |
889 | fshift[t][m] += fff; |
890 | fshift[CENTRAL(((2*1 +1)*(2*1 +1)*(2*2 +1))/2)][m] -= fff; |
891 | } /* 15 */ |
892 | |
893 | return v0*v1; /* 1 */ |
894 | /* 54 */ |
895 | } |
896 | |
897 | real thole_pol(int nbonds, |
898 | const t_iatom forceatoms[], const t_iparams forceparams[], |
899 | const rvec x[], rvec f[], rvec fshift[], |
900 | const t_pbc *pbc, const t_graph gmx_unused__attribute__ ((unused)) *g, |
901 | real gmx_unused__attribute__ ((unused)) lambda, real gmx_unused__attribute__ ((unused)) *dvdlambda, |
902 | const t_mdatoms *md, t_fcdata gmx_unused__attribute__ ((unused)) *fcd, |
903 | int gmx_unused__attribute__ ((unused)) *global_atom_index) |
904 | { |
905 | /* Interaction between two pairs of particles with opposite charge */ |
906 | int i, type, a1, da1, a2, da2; |
907 | real q1, q2, qq, a, al1, al2, afac; |
908 | real V = 0; |
909 | |
910 | for (i = 0; (i < nbonds); ) |
911 | { |
912 | type = forceatoms[i++]; |
913 | a1 = forceatoms[i++]; |
914 | da1 = forceatoms[i++]; |
915 | a2 = forceatoms[i++]; |
916 | da2 = forceatoms[i++]; |
917 | q1 = md->chargeA[da1]; |
918 | q2 = md->chargeA[da2]; |
919 | a = forceparams[type].thole.a; |
920 | al1 = forceparams[type].thole.alpha1; |
921 | al2 = forceparams[type].thole.alpha2; |
922 | qq = q1*q2; |
923 | afac = a*pow(al1*al2, -1.0/6.0); |
924 | V += do_1_thole(x[a1], x[a2], f[a1], f[a2], pbc, qq, fshift, afac); |
925 | V += do_1_thole(x[da1], x[a2], f[da1], f[a2], pbc, -qq, fshift, afac); |
926 | V += do_1_thole(x[a1], x[da2], f[a1], f[da2], pbc, -qq, fshift, afac); |
927 | V += do_1_thole(x[da1], x[da2], f[da1], f[da2], pbc, qq, fshift, afac); |
928 | } |
929 | /* 290 flops */ |
930 | return V; |
931 | } |
932 | |
933 | real bond_angle(const rvec xi, const rvec xj, const rvec xk, const t_pbc *pbc, |
934 | rvec r_ij, rvec r_kj, real *costh, |
935 | int *t1, int *t2) |
936 | /* Return value is the angle between the bonds i-j and j-k */ |
937 | { |
938 | /* 41 FLOPS */ |
939 | real th; |
940 | |
941 | *t1 = pbc_rvec_sub(pbc, xi, xj, r_ij); /* 3 */ |
942 | *t2 = pbc_rvec_sub(pbc, xk, xj, r_kj); /* 3 */ |
943 | |
944 | *costh = cos_angle(r_ij, r_kj); /* 25 */ |
945 | th = acos(*costh); /* 10 */ |
946 | /* 41 TOTAL */ |
947 | return th; |
948 | } |
949 | |
950 | real angles(int nbonds, |
951 | const t_iatom forceatoms[], const t_iparams forceparams[], |
952 | const rvec x[], rvec f[], rvec fshift[], |
953 | const t_pbc *pbc, const t_graph *g, |
954 | real lambda, real *dvdlambda, |
955 | const t_mdatoms gmx_unused__attribute__ ((unused)) *md, t_fcdata gmx_unused__attribute__ ((unused)) *fcd, |
956 | int gmx_unused__attribute__ ((unused)) *global_atom_index) |
957 | { |
958 | int i, ai, aj, ak, t1, t2, type; |
959 | rvec r_ij, r_kj; |
960 | real cos_theta, cos_theta2, theta, dVdt, va, vtot; |
961 | ivec jt, dt_ij, dt_kj; |
962 | |
963 | vtot = 0.0; |
964 | for (i = 0; i < nbonds; ) |
965 | { |
966 | type = forceatoms[i++]; |
967 | ai = forceatoms[i++]; |
968 | aj = forceatoms[i++]; |
969 | ak = forceatoms[i++]; |
970 | |
971 | theta = bond_angle(x[ai], x[aj], x[ak], pbc, |
972 | r_ij, r_kj, &cos_theta, &t1, &t2); /* 41 */ |
973 | |
974 | *dvdlambda += harmonic(forceparams[type].harmonic.krA, |
975 | forceparams[type].harmonic.krB, |
976 | forceparams[type].harmonic.rA*DEG2RAD(3.14159265358979323846/180.0), |
977 | forceparams[type].harmonic.rB*DEG2RAD(3.14159265358979323846/180.0), |
978 | theta, lambda, &va, &dVdt); /* 21 */ |
979 | vtot += va; |
980 | |
981 | cos_theta2 = sqr(cos_theta); |
982 | if (cos_theta2 < 1) |
983 | { |
984 | int m; |
985 | real st, sth; |
986 | real cik, cii, ckk; |
987 | real nrkj2, nrij2; |
988 | real nrkj_1, nrij_1; |
989 | rvec f_i, f_j, f_k; |
990 | |
991 | st = dVdt*gmx_invsqrt(1 - cos_theta2)gmx_software_invsqrt(1 - cos_theta2); /* 12 */ |
992 | sth = st*cos_theta; /* 1 */ |
993 | #ifdef DEBUG |
994 | if (debug) |
995 | { |
996 | fprintf(debug, "ANGLES: theta = %10g vth = %10g dV/dtheta = %10g\n", |
997 | theta*RAD2DEG(180.0/3.14159265358979323846), va, dVdt); |
998 | } |
999 | #endif |
1000 | nrij2 = iprod(r_ij, r_ij); /* 5 */ |
1001 | nrkj2 = iprod(r_kj, r_kj); /* 5 */ |
1002 | |
1003 | nrij_1 = gmx_invsqrt(nrij2)gmx_software_invsqrt(nrij2); /* 10 */ |
1004 | nrkj_1 = gmx_invsqrt(nrkj2)gmx_software_invsqrt(nrkj2); /* 10 */ |
1005 | |
1006 | cik = st*nrij_1*nrkj_1; /* 2 */ |
1007 | cii = sth*nrij_1*nrij_1; /* 2 */ |
1008 | ckk = sth*nrkj_1*nrkj_1; /* 2 */ |
1009 | |
1010 | for (m = 0; m < DIM3; m++) |
1011 | { /* 39 */ |
1012 | f_i[m] = -(cik*r_kj[m] - cii*r_ij[m]); |
1013 | f_k[m] = -(cik*r_ij[m] - ckk*r_kj[m]); |
1014 | f_j[m] = -f_i[m] - f_k[m]; |
1015 | f[ai][m] += f_i[m]; |
1016 | f[aj][m] += f_j[m]; |
1017 | f[ak][m] += f_k[m]; |
1018 | } |
1019 | if (g != NULL((void*)0)) |
1020 | { |
1021 | copy_ivec(SHIFT_IVEC(g, aj)((g)->ishift[aj]), jt); |
1022 | |
1023 | ivec_sub(SHIFT_IVEC(g, ai)((g)->ishift[ai]), jt, dt_ij); |
1024 | ivec_sub(SHIFT_IVEC(g, ak)((g)->ishift[ak]), jt, dt_kj); |
1025 | t1 = IVEC2IS(dt_ij)(((2*2 +1)*((2*1 +1)*(((dt_ij)[2])+1)+((dt_ij)[1])+1)+((dt_ij )[0])+2)); |
1026 | t2 = IVEC2IS(dt_kj)(((2*2 +1)*((2*1 +1)*(((dt_kj)[2])+1)+((dt_kj)[1])+1)+((dt_kj )[0])+2)); |
1027 | } |
1028 | rvec_inc(fshift[t1], f_i); |
1029 | rvec_inc(fshift[CENTRAL(((2*1 +1)*(2*1 +1)*(2*2 +1))/2)], f_j); |
1030 | rvec_inc(fshift[t2], f_k); |
1031 | } /* 161 TOTAL */ |
1032 | } |
1033 | |
1034 | return vtot; |
1035 | } |
1036 | |
1037 | #ifdef GMX_SIMD_HAVE_REAL |
1038 | |
1039 | /* As angles, but using SIMD to calculate many dihedrals at once. |
1040 | * This routines does not calculate energies and shift forces. |
1041 | */ |
1042 | static gmx_inlineinline void |
1043 | angles_noener_simd(int nbonds, |
1044 | const t_iatom forceatoms[], const t_iparams forceparams[], |
1045 | const rvec x[], rvec f[], |
1046 | const t_pbc *pbc, const t_graph gmx_unused__attribute__ ((unused)) *g, |
1047 | real gmx_unused__attribute__ ((unused)) lambda, |
1048 | const t_mdatoms gmx_unused__attribute__ ((unused)) *md, t_fcdata gmx_unused__attribute__ ((unused)) *fcd, |
1049 | int gmx_unused__attribute__ ((unused)) *global_atom_index) |
1050 | { |
1051 | const int nfa1 = 4; |
1052 | int i, iu, s, m; |
1053 | int type, ai[GMX_SIMD_REAL_WIDTH4], aj[GMX_SIMD_REAL_WIDTH4]; |
1054 | int ak[GMX_SIMD_REAL_WIDTH4]; |
1055 | real coeff_array[2*GMX_SIMD_REAL_WIDTH4+GMX_SIMD_REAL_WIDTH4], *coeff; |
1056 | real dr_array[2*DIM3*GMX_SIMD_REAL_WIDTH4+GMX_SIMD_REAL_WIDTH4], *dr; |
1057 | real f_buf_array[6*GMX_SIMD_REAL_WIDTH4+GMX_SIMD_REAL_WIDTH4], *f_buf; |
1058 | gmx_simd_real_t__m128 k_S, theta0_S; |
1059 | gmx_simd_real_t__m128 rijx_S, rijy_S, rijz_S; |
1060 | gmx_simd_real_t__m128 rkjx_S, rkjy_S, rkjz_S; |
1061 | gmx_simd_real_t__m128 one_S; |
1062 | gmx_simd_real_t__m128 min_one_plus_eps_S; |
1063 | gmx_simd_real_t__m128 rij_rkj_S; |
1064 | gmx_simd_real_t__m128 nrij2_S, nrij_1_S; |
1065 | gmx_simd_real_t__m128 nrkj2_S, nrkj_1_S; |
1066 | gmx_simd_real_t__m128 cos_S, invsin_S; |
1067 | gmx_simd_real_t__m128 theta_S; |
1068 | gmx_simd_real_t__m128 st_S, sth_S; |
1069 | gmx_simd_real_t__m128 cik_S, cii_S, ckk_S; |
1070 | gmx_simd_real_t__m128 f_ix_S, f_iy_S, f_iz_S; |
1071 | gmx_simd_real_t__m128 f_kx_S, f_ky_S, f_kz_S; |
1072 | pbc_simd_t pbc_simd; |
1073 | |
1074 | /* Ensure register memory alignment */ |
1075 | coeff = gmx_simd_align_rgmx_simd_align_f(coeff_array); |
1076 | dr = gmx_simd_align_rgmx_simd_align_f(dr_array); |
1077 | f_buf = gmx_simd_align_rgmx_simd_align_f(f_buf_array); |
1078 | |
1079 | set_pbc_simd(pbc, &pbc_simd); |
1080 | |
1081 | one_S = gmx_simd_set1_r_mm_set1_ps(1.0); |
1082 | |
1083 | /* The smallest number > -1 */ |
1084 | min_one_plus_eps_S = gmx_simd_set1_r_mm_set1_ps(-1.0 + 2*GMX_REAL_EPS5.96046448E-08); |
1085 | |
1086 | /* nbonds is the number of angles times nfa1, here we step GMX_SIMD_REAL_WIDTH angles */ |
1087 | for (i = 0; (i < nbonds); i += GMX_SIMD_REAL_WIDTH4*nfa1) |
1088 | { |
1089 | /* Collect atoms for GMX_SIMD_REAL_WIDTH angles. |
1090 | * iu indexes into forceatoms, we should not let iu go beyond nbonds. |
1091 | */ |
1092 | iu = i; |
1093 | for (s = 0; s < GMX_SIMD_REAL_WIDTH4; s++) |
1094 | { |
1095 | type = forceatoms[iu]; |
1096 | ai[s] = forceatoms[iu+1]; |
1097 | aj[s] = forceatoms[iu+2]; |
1098 | ak[s] = forceatoms[iu+3]; |
1099 | |
1100 | coeff[s] = forceparams[type].harmonic.krA; |
1101 | coeff[GMX_SIMD_REAL_WIDTH4+s] = forceparams[type].harmonic.rA*DEG2RAD(3.14159265358979323846/180.0); |
1102 | |
1103 | /* If you can't use pbc_dx_simd below for PBC, e.g. because |
1104 | * you can't round in SIMD, use pbc_rvec_sub here. |
1105 | */ |
1106 | /* Store the non PBC corrected distances packed and aligned */ |
1107 | for (m = 0; m < DIM3; m++) |
1108 | { |
1109 | dr[s + m *GMX_SIMD_REAL_WIDTH4] = x[ai[s]][m] - x[aj[s]][m]; |
1110 | dr[s + (DIM3+m)*GMX_SIMD_REAL_WIDTH4] = x[ak[s]][m] - x[aj[s]][m]; |
1111 | } |
1112 | |
1113 | /* At the end fill the arrays with identical entries */ |
1114 | if (iu + nfa1 < nbonds) |
1115 | { |
1116 | iu += nfa1; |
1117 | } |
1118 | } |
1119 | |
1120 | k_S = gmx_simd_load_r_mm_load_ps(coeff); |
1121 | theta0_S = gmx_simd_load_r_mm_load_ps(coeff+GMX_SIMD_REAL_WIDTH4); |
1122 | |
1123 | rijx_S = gmx_simd_load_r_mm_load_ps(dr + 0*GMX_SIMD_REAL_WIDTH4); |
1124 | rijy_S = gmx_simd_load_r_mm_load_ps(dr + 1*GMX_SIMD_REAL_WIDTH4); |
1125 | rijz_S = gmx_simd_load_r_mm_load_ps(dr + 2*GMX_SIMD_REAL_WIDTH4); |
1126 | rkjx_S = gmx_simd_load_r_mm_load_ps(dr + 3*GMX_SIMD_REAL_WIDTH4); |
1127 | rkjy_S = gmx_simd_load_r_mm_load_ps(dr + 4*GMX_SIMD_REAL_WIDTH4); |
1128 | rkjz_S = gmx_simd_load_r_mm_load_ps(dr + 5*GMX_SIMD_REAL_WIDTH4); |
1129 | |
1130 | pbc_dx_simd(&rijx_S, &rijy_S, &rijz_S, &pbc_simd); |
1131 | pbc_dx_simd(&rkjx_S, &rkjy_S, &rkjz_S, &pbc_simd); |
1132 | |
1133 | rij_rkj_S = gmx_simd_iprod_rgmx_simd_iprod_f(rijx_S, rijy_S, rijz_S, |
1134 | rkjx_S, rkjy_S, rkjz_S); |
1135 | |
1136 | nrij2_S = gmx_simd_norm2_rgmx_simd_norm2_f(rijx_S, rijy_S, rijz_S); |
1137 | nrkj2_S = gmx_simd_norm2_rgmx_simd_norm2_f(rkjx_S, rkjy_S, rkjz_S); |
1138 | |
1139 | nrij_1_S = gmx_simd_invsqrt_rgmx_simd_invsqrt_f(nrij2_S); |
1140 | nrkj_1_S = gmx_simd_invsqrt_rgmx_simd_invsqrt_f(nrkj2_S); |
1141 | |
1142 | cos_S = gmx_simd_mul_r_mm_mul_ps(rij_rkj_S, gmx_simd_mul_r_mm_mul_ps(nrij_1_S, nrkj_1_S)); |
1143 | |
1144 | /* To allow for 180 degrees, we take the max of cos and -1 + 1bit, |
1145 | * so we can safely get the 1/sin from 1/sqrt(1 - cos^2). |
1146 | * This also ensures that rounding errors would cause the argument |
1147 | * of gmx_simd_acos_r to be < -1. |
1148 | * Note that we do not take precautions for cos(0)=1, so the outer |
1149 | * atoms in an angle should not be on top of each other. |
1150 | */ |
1151 | cos_S = gmx_simd_max_r_mm_max_ps(cos_S, min_one_plus_eps_S); |
1152 | |
1153 | theta_S = gmx_simd_acos_rgmx_simd_acos_f(cos_S); |
1154 | |
1155 | invsin_S = gmx_simd_invsqrt_rgmx_simd_invsqrt_f(gmx_simd_sub_r_mm_sub_ps(one_S, gmx_simd_mul_r_mm_mul_ps(cos_S, cos_S))); |
1156 | |
1157 | st_S = gmx_simd_mul_r_mm_mul_ps(gmx_simd_mul_r_mm_mul_ps(k_S, gmx_simd_sub_r_mm_sub_ps(theta0_S, theta_S)), |
1158 | invsin_S); |
1159 | sth_S = gmx_simd_mul_r_mm_mul_ps(st_S, cos_S); |
1160 | |
1161 | cik_S = gmx_simd_mul_r_mm_mul_ps(st_S, gmx_simd_mul_r_mm_mul_ps(nrij_1_S, nrkj_1_S)); |
1162 | cii_S = gmx_simd_mul_r_mm_mul_ps(sth_S, gmx_simd_mul_r_mm_mul_ps(nrij_1_S, nrij_1_S)); |
1163 | ckk_S = gmx_simd_mul_r_mm_mul_ps(sth_S, gmx_simd_mul_r_mm_mul_ps(nrkj_1_S, nrkj_1_S)); |
1164 | |
1165 | f_ix_S = gmx_simd_mul_r_mm_mul_ps(cii_S, rijx_S); |
1166 | f_ix_S = gmx_simd_fnmadd_r(cik_S, rkjx_S, f_ix_S)_mm_sub_ps(f_ix_S, _mm_mul_ps(cik_S, rkjx_S)); |
1167 | f_iy_S = gmx_simd_mul_r_mm_mul_ps(cii_S, rijy_S); |
1168 | f_iy_S = gmx_simd_fnmadd_r(cik_S, rkjy_S, f_iy_S)_mm_sub_ps(f_iy_S, _mm_mul_ps(cik_S, rkjy_S)); |
1169 | f_iz_S = gmx_simd_mul_r_mm_mul_ps(cii_S, rijz_S); |
1170 | f_iz_S = gmx_simd_fnmadd_r(cik_S, rkjz_S, f_iz_S)_mm_sub_ps(f_iz_S, _mm_mul_ps(cik_S, rkjz_S)); |
1171 | f_kx_S = gmx_simd_mul_r_mm_mul_ps(ckk_S, rkjx_S); |
1172 | f_kx_S = gmx_simd_fnmadd_r(cik_S, rijx_S, f_kx_S)_mm_sub_ps(f_kx_S, _mm_mul_ps(cik_S, rijx_S)); |
1173 | f_ky_S = gmx_simd_mul_r_mm_mul_ps(ckk_S, rkjy_S); |
1174 | f_ky_S = gmx_simd_fnmadd_r(cik_S, rijy_S, f_ky_S)_mm_sub_ps(f_ky_S, _mm_mul_ps(cik_S, rijy_S)); |
1175 | f_kz_S = gmx_simd_mul_r_mm_mul_ps(ckk_S, rkjz_S); |
1176 | f_kz_S = gmx_simd_fnmadd_r(cik_S, rijz_S, f_kz_S)_mm_sub_ps(f_kz_S, _mm_mul_ps(cik_S, rijz_S)); |
1177 | |
1178 | gmx_simd_store_r_mm_store_ps(f_buf + 0*GMX_SIMD_REAL_WIDTH4, f_ix_S); |
1179 | gmx_simd_store_r_mm_store_ps(f_buf + 1*GMX_SIMD_REAL_WIDTH4, f_iy_S); |
1180 | gmx_simd_store_r_mm_store_ps(f_buf + 2*GMX_SIMD_REAL_WIDTH4, f_iz_S); |
1181 | gmx_simd_store_r_mm_store_ps(f_buf + 3*GMX_SIMD_REAL_WIDTH4, f_kx_S); |
1182 | gmx_simd_store_r_mm_store_ps(f_buf + 4*GMX_SIMD_REAL_WIDTH4, f_ky_S); |
1183 | gmx_simd_store_r_mm_store_ps(f_buf + 5*GMX_SIMD_REAL_WIDTH4, f_kz_S); |
1184 | |
1185 | iu = i; |
1186 | s = 0; |
1187 | do |
1188 | { |
1189 | for (m = 0; m < DIM3; m++) |
1190 | { |
1191 | f[ai[s]][m] += f_buf[s + m*GMX_SIMD_REAL_WIDTH4]; |
1192 | f[aj[s]][m] -= f_buf[s + m*GMX_SIMD_REAL_WIDTH4] + f_buf[s + (DIM3+m)*GMX_SIMD_REAL_WIDTH4]; |
1193 | f[ak[s]][m] += f_buf[s + (DIM3+m)*GMX_SIMD_REAL_WIDTH4]; |
1194 | } |
1195 | s++; |
1196 | iu += nfa1; |
1197 | } |
1198 | while (s < GMX_SIMD_REAL_WIDTH4 && iu < nbonds); |
1199 | } |
1200 | } |
1201 | |
1202 | #endif /* GMX_SIMD_HAVE_REAL */ |
1203 | |
1204 | real linear_angles(int nbonds, |
1205 | const t_iatom forceatoms[], const t_iparams forceparams[], |
1206 | const rvec x[], rvec f[], rvec fshift[], |
1207 | const t_pbc *pbc, const t_graph *g, |
1208 | real lambda, real *dvdlambda, |
1209 | const t_mdatoms gmx_unused__attribute__ ((unused)) *md, t_fcdata gmx_unused__attribute__ ((unused)) *fcd, |
1210 | int gmx_unused__attribute__ ((unused)) *global_atom_index) |
1211 | { |
1212 | int i, m, ai, aj, ak, t1, t2, type; |
1213 | rvec f_i, f_j, f_k; |
1214 | real L1, kA, kB, aA, aB, dr, dr2, va, vtot, a, b, klin; |
1215 | ivec jt, dt_ij, dt_kj; |
1216 | rvec r_ij, r_kj, r_ik, dx; |
1217 | |
1218 | L1 = 1-lambda; |
1219 | vtot = 0.0; |
1220 | for (i = 0; (i < nbonds); ) |
1221 | { |
1222 | type = forceatoms[i++]; |
1223 | ai = forceatoms[i++]; |
1224 | aj = forceatoms[i++]; |
1225 | ak = forceatoms[i++]; |
1226 | |
1227 | kA = forceparams[type].linangle.klinA; |
1228 | kB = forceparams[type].linangle.klinB; |
1229 | klin = L1*kA + lambda*kB; |
1230 | |
1231 | aA = forceparams[type].linangle.aA; |
1232 | aB = forceparams[type].linangle.aB; |
1233 | a = L1*aA+lambda*aB; |
1234 | b = 1-a; |
1235 | |
1236 | t1 = pbc_rvec_sub(pbc, x[ai], x[aj], r_ij); |
1237 | t2 = pbc_rvec_sub(pbc, x[ak], x[aj], r_kj); |
1238 | rvec_sub(r_ij, r_kj, r_ik); |
1239 | |
1240 | dr2 = 0; |
1241 | for (m = 0; (m < DIM3); m++) |
1242 | { |
1243 | dr = -a * r_ij[m] - b * r_kj[m]; |
1244 | dr2 += dr*dr; |
1245 | dx[m] = dr; |
1246 | f_i[m] = a*klin*dr; |
1247 | f_k[m] = b*klin*dr; |
1248 | f_j[m] = -(f_i[m]+f_k[m]); |
1249 | f[ai][m] += f_i[m]; |
1250 | f[aj][m] += f_j[m]; |
1251 | f[ak][m] += f_k[m]; |
1252 | } |
1253 | va = 0.5*klin*dr2; |
1254 | *dvdlambda += 0.5*(kB-kA)*dr2 + klin*(aB-aA)*iprod(dx, r_ik); |
1255 | |
1256 | vtot += va; |
1257 | |
1258 | if (g) |
1259 | { |
1260 | copy_ivec(SHIFT_IVEC(g, aj)((g)->ishift[aj]), jt); |
1261 | |
1262 | ivec_sub(SHIFT_IVEC(g, ai)((g)->ishift[ai]), jt, dt_ij); |
1263 | ivec_sub(SHIFT_IVEC(g, ak)((g)->ishift[ak]), jt, dt_kj); |
1264 | t1 = IVEC2IS(dt_ij)(((2*2 +1)*((2*1 +1)*(((dt_ij)[2])+1)+((dt_ij)[1])+1)+((dt_ij )[0])+2)); |
1265 | t2 = IVEC2IS(dt_kj)(((2*2 +1)*((2*1 +1)*(((dt_kj)[2])+1)+((dt_kj)[1])+1)+((dt_kj )[0])+2)); |
1266 | } |
1267 | rvec_inc(fshift[t1], f_i); |
1268 | rvec_inc(fshift[CENTRAL(((2*1 +1)*(2*1 +1)*(2*2 +1))/2)], f_j); |
1269 | rvec_inc(fshift[t2], f_k); |
1270 | } /* 57 TOTAL */ |
1271 | return vtot; |
1272 | } |
1273 | |
1274 | real urey_bradley(int nbonds, |
1275 | const t_iatom forceatoms[], const t_iparams forceparams[], |
1276 | const rvec x[], rvec f[], rvec fshift[], |
1277 | const t_pbc *pbc, const t_graph *g, |
1278 | real lambda, real *dvdlambda, |
1279 | const t_mdatoms gmx_unused__attribute__ ((unused)) *md, t_fcdata gmx_unused__attribute__ ((unused)) *fcd, |
1280 | int gmx_unused__attribute__ ((unused)) *global_atom_index) |
1281 | { |
1282 | int i, m, ai, aj, ak, t1, t2, type, ki; |
1283 | rvec r_ij, r_kj, r_ik; |
1284 | real cos_theta, cos_theta2, theta; |
1285 | real dVdt, va, vtot, dr, dr2, vbond, fbond, fik; |
1286 | real kthA, th0A, kUBA, r13A, kthB, th0B, kUBB, r13B; |
1287 | ivec jt, dt_ij, dt_kj, dt_ik; |
1288 | |
1289 | vtot = 0.0; |
1290 | for (i = 0; (i < nbonds); ) |
1291 | { |
1292 | type = forceatoms[i++]; |
1293 | ai = forceatoms[i++]; |
1294 | aj = forceatoms[i++]; |
1295 | ak = forceatoms[i++]; |
1296 | th0A = forceparams[type].u_b.thetaA*DEG2RAD(3.14159265358979323846/180.0); |
1297 | kthA = forceparams[type].u_b.kthetaA; |
1298 | r13A = forceparams[type].u_b.r13A; |
1299 | kUBA = forceparams[type].u_b.kUBA; |
1300 | th0B = forceparams[type].u_b.thetaB*DEG2RAD(3.14159265358979323846/180.0); |
1301 | kthB = forceparams[type].u_b.kthetaB; |
1302 | r13B = forceparams[type].u_b.r13B; |
1303 | kUBB = forceparams[type].u_b.kUBB; |
1304 | |
1305 | theta = bond_angle(x[ai], x[aj], x[ak], pbc, |
1306 | r_ij, r_kj, &cos_theta, &t1, &t2); /* 41 */ |
1307 | |
1308 | *dvdlambda += harmonic(kthA, kthB, th0A, th0B, theta, lambda, &va, &dVdt); /* 21 */ |
1309 | vtot += va; |
1310 | |
1311 | ki = pbc_rvec_sub(pbc, x[ai], x[ak], r_ik); /* 3 */ |
1312 | dr2 = iprod(r_ik, r_ik); /* 5 */ |
1313 | dr = dr2*gmx_invsqrt(dr2)gmx_software_invsqrt(dr2); /* 10 */ |
1314 | |
1315 | *dvdlambda += harmonic(kUBA, kUBB, r13A, r13B, dr, lambda, &vbond, &fbond); /* 19 */ |
1316 | |
1317 | cos_theta2 = sqr(cos_theta); /* 1 */ |
1318 | if (cos_theta2 < 1) |
1319 | { |
1320 | real st, sth; |
1321 | real cik, cii, ckk; |
1322 | real nrkj2, nrij2; |
1323 | rvec f_i, f_j, f_k; |
1324 | |
1325 | st = dVdt*gmx_invsqrt(1 - cos_theta2)gmx_software_invsqrt(1 - cos_theta2); /* 12 */ |
1326 | sth = st*cos_theta; /* 1 */ |
1327 | #ifdef DEBUG |
1328 | if (debug) |
1329 | { |
1330 | fprintf(debug, "ANGLES: theta = %10g vth = %10g dV/dtheta = %10g\n", |
1331 | theta*RAD2DEG(180.0/3.14159265358979323846), va, dVdt); |
1332 | } |
1333 | #endif |
1334 | nrkj2 = iprod(r_kj, r_kj); /* 5 */ |
1335 | nrij2 = iprod(r_ij, r_ij); |
1336 | |
1337 | cik = st*gmx_invsqrt(nrkj2*nrij2)gmx_software_invsqrt(nrkj2*nrij2); /* 12 */ |
1338 | cii = sth/nrij2; /* 10 */ |
1339 | ckk = sth/nrkj2; /* 10 */ |
1340 | |
1341 | for (m = 0; (m < DIM3); m++) /* 39 */ |
1342 | { |
1343 | f_i[m] = -(cik*r_kj[m]-cii*r_ij[m]); |
1344 | f_k[m] = -(cik*r_ij[m]-ckk*r_kj[m]); |
1345 | f_j[m] = -f_i[m]-f_k[m]; |
1346 | f[ai][m] += f_i[m]; |
1347 | f[aj][m] += f_j[m]; |
1348 | f[ak][m] += f_k[m]; |
1349 | } |
1350 | if (g) |
1351 | { |
1352 | copy_ivec(SHIFT_IVEC(g, aj)((g)->ishift[aj]), jt); |
1353 | |
1354 | ivec_sub(SHIFT_IVEC(g, ai)((g)->ishift[ai]), jt, dt_ij); |
1355 | ivec_sub(SHIFT_IVEC(g, ak)((g)->ishift[ak]), jt, dt_kj); |
1356 | t1 = IVEC2IS(dt_ij)(((2*2 +1)*((2*1 +1)*(((dt_ij)[2])+1)+((dt_ij)[1])+1)+((dt_ij )[0])+2)); |
1357 | t2 = IVEC2IS(dt_kj)(((2*2 +1)*((2*1 +1)*(((dt_kj)[2])+1)+((dt_kj)[1])+1)+((dt_kj )[0])+2)); |
1358 | } |
1359 | rvec_inc(fshift[t1], f_i); |
1360 | rvec_inc(fshift[CENTRAL(((2*1 +1)*(2*1 +1)*(2*2 +1))/2)], f_j); |
1361 | rvec_inc(fshift[t2], f_k); |
1362 | } /* 161 TOTAL */ |
1363 | /* Time for the bond calculations */ |
1364 | if (dr2 == 0.0) |
1365 | { |
1366 | continue; |
1367 | } |
1368 | |
1369 | vtot += vbond; /* 1*/ |
1370 | fbond *= gmx_invsqrt(dr2)gmx_software_invsqrt(dr2); /* 6 */ |
1371 | |
1372 | if (g) |
1373 | { |
1374 | ivec_sub(SHIFT_IVEC(g, ai)((g)->ishift[ai]), SHIFT_IVEC(g, ak)((g)->ishift[ak]), dt_ik); |
1375 | ki = IVEC2IS(dt_ik)(((2*2 +1)*((2*1 +1)*(((dt_ik)[2])+1)+((dt_ik)[1])+1)+((dt_ik )[0])+2)); |
1376 | } |
1377 | for (m = 0; (m < DIM3); m++) /* 15 */ |
1378 | { |
1379 | fik = fbond*r_ik[m]; |
1380 | f[ai][m] += fik; |
1381 | f[ak][m] -= fik; |
1382 | fshift[ki][m] += fik; |
1383 | fshift[CENTRAL(((2*1 +1)*(2*1 +1)*(2*2 +1))/2)][m] -= fik; |
1384 | } |
1385 | } |
1386 | return vtot; |
1387 | } |
1388 | |
1389 | real quartic_angles(int nbonds, |
1390 | const t_iatom forceatoms[], const t_iparams forceparams[], |
1391 | const rvec x[], rvec f[], rvec fshift[], |
1392 | const t_pbc *pbc, const t_graph *g, |
1393 | real gmx_unused__attribute__ ((unused)) lambda, real gmx_unused__attribute__ ((unused)) *dvdlambda, |
1394 | const t_mdatoms gmx_unused__attribute__ ((unused)) *md, t_fcdata gmx_unused__attribute__ ((unused)) *fcd, |
1395 | int gmx_unused__attribute__ ((unused)) *global_atom_index) |
1396 | { |
1397 | int i, j, ai, aj, ak, t1, t2, type; |
1398 | rvec r_ij, r_kj; |
1399 | real cos_theta, cos_theta2, theta, dt, dVdt, va, dtp, c, vtot; |
1400 | ivec jt, dt_ij, dt_kj; |
1401 | |
1402 | vtot = 0.0; |
1403 | for (i = 0; (i < nbonds); ) |
1404 | { |
1405 | type = forceatoms[i++]; |
1406 | ai = forceatoms[i++]; |
1407 | aj = forceatoms[i++]; |
1408 | ak = forceatoms[i++]; |
1409 | |
1410 | theta = bond_angle(x[ai], x[aj], x[ak], pbc, |
1411 | r_ij, r_kj, &cos_theta, &t1, &t2); /* 41 */ |
1412 | |
1413 | dt = theta - forceparams[type].qangle.theta*DEG2RAD(3.14159265358979323846/180.0); /* 2 */ |
1414 | |
1415 | dVdt = 0; |
1416 | va = forceparams[type].qangle.c[0]; |
1417 | dtp = 1.0; |
1418 | for (j = 1; j <= 4; j++) |
1419 | { |
1420 | c = forceparams[type].qangle.c[j]; |
1421 | dVdt -= j*c*dtp; |
1422 | dtp *= dt; |
1423 | va += c*dtp; |
1424 | } |
1425 | /* 20 */ |
1426 | |
1427 | vtot += va; |
1428 | |
1429 | cos_theta2 = sqr(cos_theta); /* 1 */ |
1430 | if (cos_theta2 < 1) |
1431 | { |
1432 | int m; |
1433 | real st, sth; |
1434 | real cik, cii, ckk; |
1435 | real nrkj2, nrij2; |
1436 | rvec f_i, f_j, f_k; |
1437 | |
1438 | st = dVdt*gmx_invsqrt(1 - cos_theta2)gmx_software_invsqrt(1 - cos_theta2); /* 12 */ |
1439 | sth = st*cos_theta; /* 1 */ |
1440 | #ifdef DEBUG |
1441 | if (debug) |
1442 | { |
1443 | fprintf(debug, "ANGLES: theta = %10g vth = %10g dV/dtheta = %10g\n", |
1444 | theta*RAD2DEG(180.0/3.14159265358979323846), va, dVdt); |
1445 | } |
1446 | #endif |
1447 | nrkj2 = iprod(r_kj, r_kj); /* 5 */ |
1448 | nrij2 = iprod(r_ij, r_ij); |
1449 | |
1450 | cik = st*gmx_invsqrt(nrkj2*nrij2)gmx_software_invsqrt(nrkj2*nrij2); /* 12 */ |
1451 | cii = sth/nrij2; /* 10 */ |
1452 | ckk = sth/nrkj2; /* 10 */ |
1453 | |
1454 | for (m = 0; (m < DIM3); m++) /* 39 */ |
1455 | { |
1456 | f_i[m] = -(cik*r_kj[m]-cii*r_ij[m]); |
1457 | f_k[m] = -(cik*r_ij[m]-ckk*r_kj[m]); |
1458 | f_j[m] = -f_i[m]-f_k[m]; |
1459 | f[ai][m] += f_i[m]; |
1460 | f[aj][m] += f_j[m]; |
1461 | f[ak][m] += f_k[m]; |
1462 | } |
1463 | if (g) |
1464 | { |
1465 | copy_ivec(SHIFT_IVEC(g, aj)((g)->ishift[aj]), jt); |
1466 | |
1467 | ivec_sub(SHIFT_IVEC(g, ai)((g)->ishift[ai]), jt, dt_ij); |
1468 | ivec_sub(SHIFT_IVEC(g, ak)((g)->ishift[ak]), jt, dt_kj); |
1469 | t1 = IVEC2IS(dt_ij)(((2*2 +1)*((2*1 +1)*(((dt_ij)[2])+1)+((dt_ij)[1])+1)+((dt_ij )[0])+2)); |
1470 | t2 = IVEC2IS(dt_kj)(((2*2 +1)*((2*1 +1)*(((dt_kj)[2])+1)+((dt_kj)[1])+1)+((dt_kj )[0])+2)); |
1471 | } |
1472 | rvec_inc(fshift[t1], f_i); |
1473 | rvec_inc(fshift[CENTRAL(((2*1 +1)*(2*1 +1)*(2*2 +1))/2)], f_j); |
1474 | rvec_inc(fshift[t2], f_k); |
1475 | } /* 153 TOTAL */ |
1476 | } |
1477 | return vtot; |
1478 | } |
1479 | |
1480 | real dih_angle(const rvec xi, const rvec xj, const rvec xk, const rvec xl, |
1481 | const t_pbc *pbc, |
1482 | rvec r_ij, rvec r_kj, rvec r_kl, rvec m, rvec n, |
1483 | real *sign, int *t1, int *t2, int *t3) |
1484 | { |
1485 | real ipr, phi; |
1486 | |
1487 | *t1 = pbc_rvec_sub(pbc, xi, xj, r_ij); /* 3 */ |
1488 | *t2 = pbc_rvec_sub(pbc, xk, xj, r_kj); /* 3 */ |
1489 | *t3 = pbc_rvec_sub(pbc, xk, xl, r_kl); /* 3 */ |
1490 | |
1491 | cprod(r_ij, r_kj, m); /* 9 */ |
1492 | cprod(r_kj, r_kl, n); /* 9 */ |
1493 | phi = gmx_angle(m, n); /* 49 (assuming 25 for atan2) */ |
1494 | ipr = iprod(r_ij, n); /* 5 */ |
1495 | (*sign) = (ipr < 0.0) ? -1.0 : 1.0; |
1496 | phi = (*sign)*phi; /* 1 */ |
1497 | /* 82 TOTAL */ |
1498 | return phi; |
1499 | } |
1500 | |
1501 | |
1502 | #ifdef GMX_SIMD_HAVE_REAL |
1503 | |
1504 | /* As dih_angle above, but calculates 4 dihedral angles at once using SIMD, |
1505 | * also calculates the pre-factor required for the dihedral force update. |
1506 | * Note that bv and buf should be register aligned. |
1507 | */ |
1508 | static gmx_inlineinline void |
1509 | dih_angle_simd(const rvec *x, |
1510 | const int *ai, const int *aj, const int *ak, const int *al, |
1511 | const pbc_simd_t *pbc, |
1512 | real *dr, |
1513 | gmx_simd_real_t__m128 *phi_S, |
1514 | gmx_simd_real_t__m128 *mx_S, gmx_simd_real_t__m128 *my_S, gmx_simd_real_t__m128 *mz_S, |
1515 | gmx_simd_real_t__m128 *nx_S, gmx_simd_real_t__m128 *ny_S, gmx_simd_real_t__m128 *nz_S, |
1516 | gmx_simd_real_t__m128 *nrkj_m2_S, |
1517 | gmx_simd_real_t__m128 *nrkj_n2_S, |
1518 | real *p, |
1519 | real *q) |
1520 | { |
1521 | int s, m; |
1522 | gmx_simd_real_t__m128 rijx_S, rijy_S, rijz_S; |
1523 | gmx_simd_real_t__m128 rkjx_S, rkjy_S, rkjz_S; |
1524 | gmx_simd_real_t__m128 rklx_S, rkly_S, rklz_S; |
1525 | gmx_simd_real_t__m128 cx_S, cy_S, cz_S; |
1526 | gmx_simd_real_t__m128 cn_S; |
1527 | gmx_simd_real_t__m128 s_S; |
1528 | gmx_simd_real_t__m128 ipr_S; |
1529 | gmx_simd_real_t__m128 iprm_S, iprn_S; |
1530 | gmx_simd_real_t__m128 nrkj2_S, nrkj_1_S, nrkj_2_S, nrkj_S; |
1531 | gmx_simd_real_t__m128 toler_S; |
1532 | gmx_simd_real_t__m128 p_S, q_S; |
1533 | gmx_simd_real_t__m128 nrkj2_min_S; |
1534 | gmx_simd_real_t__m128 real_eps_S; |
1535 | |
1536 | /* Used to avoid division by zero. |
1537 | * We take into acount that we multiply the result by real_eps_S. |
1538 | */ |
1539 | nrkj2_min_S = gmx_simd_set1_r_mm_set1_ps(GMX_REAL_MIN1.17549435E-38/(2*GMX_REAL_EPS5.96046448E-08)); |
1540 | |
1541 | /* The value of the last significant bit (GMX_REAL_EPS is half of that) */ |
1542 | real_eps_S = gmx_simd_set1_r_mm_set1_ps(2*GMX_REAL_EPS5.96046448E-08); |
1543 | |
1544 | for (s = 0; s < GMX_SIMD_REAL_WIDTH4; s++) |
1545 | { |
1546 | /* If you can't use pbc_dx_simd below for PBC, e.g. because |
1547 | * you can't round in SIMD, use pbc_rvec_sub here. |
1548 | */ |
1549 | for (m = 0; m < DIM3; m++) |
1550 | { |
1551 | dr[s + (0*DIM3 + m)*GMX_SIMD_REAL_WIDTH4] = x[ai[s]][m] - x[aj[s]][m]; |
1552 | dr[s + (1*DIM3 + m)*GMX_SIMD_REAL_WIDTH4] = x[ak[s]][m] - x[aj[s]][m]; |
1553 | dr[s + (2*DIM3 + m)*GMX_SIMD_REAL_WIDTH4] = x[ak[s]][m] - x[al[s]][m]; |
1554 | } |
1555 | } |
1556 | |
1557 | rijx_S = gmx_simd_load_r_mm_load_ps(dr + 0*GMX_SIMD_REAL_WIDTH4); |
1558 | rijy_S = gmx_simd_load_r_mm_load_ps(dr + 1*GMX_SIMD_REAL_WIDTH4); |
1559 | rijz_S = gmx_simd_load_r_mm_load_ps(dr + 2*GMX_SIMD_REAL_WIDTH4); |
1560 | rkjx_S = gmx_simd_load_r_mm_load_ps(dr + 3*GMX_SIMD_REAL_WIDTH4); |
1561 | rkjy_S = gmx_simd_load_r_mm_load_ps(dr + 4*GMX_SIMD_REAL_WIDTH4); |
1562 | rkjz_S = gmx_simd_load_r_mm_load_ps(dr + 5*GMX_SIMD_REAL_WIDTH4); |
1563 | rklx_S = gmx_simd_load_r_mm_load_ps(dr + 6*GMX_SIMD_REAL_WIDTH4); |
1564 | rkly_S = gmx_simd_load_r_mm_load_ps(dr + 7*GMX_SIMD_REAL_WIDTH4); |
1565 | rklz_S = gmx_simd_load_r_mm_load_ps(dr + 8*GMX_SIMD_REAL_WIDTH4); |
1566 | |
1567 | pbc_dx_simd(&rijx_S, &rijy_S, &rijz_S, pbc); |
1568 | pbc_dx_simd(&rkjx_S, &rkjy_S, &rkjz_S, pbc); |
1569 | pbc_dx_simd(&rklx_S, &rkly_S, &rklz_S, pbc); |
1570 | |
1571 | gmx_simd_cprod_rgmx_simd_cprod_f(rijx_S, rijy_S, rijz_S, |
1572 | rkjx_S, rkjy_S, rkjz_S, |
1573 | mx_S, my_S, mz_S); |
1574 | |
1575 | gmx_simd_cprod_rgmx_simd_cprod_f(rkjx_S, rkjy_S, rkjz_S, |
1576 | rklx_S, rkly_S, rklz_S, |
1577 | nx_S, ny_S, nz_S); |
1578 | |
1579 | gmx_simd_cprod_rgmx_simd_cprod_f(*mx_S, *my_S, *mz_S, |
1580 | *nx_S, *ny_S, *nz_S, |
1581 | &cx_S, &cy_S, &cz_S); |
1582 | |
1583 | cn_S = gmx_simd_sqrt_rgmx_simd_sqrt_f(gmx_simd_norm2_rgmx_simd_norm2_f(cx_S, cy_S, cz_S)); |
1584 | |
1585 | s_S = gmx_simd_iprod_rgmx_simd_iprod_f(*mx_S, *my_S, *mz_S, *nx_S, *ny_S, *nz_S); |
1586 | |
1587 | /* Determine the dihedral angle, the sign might need correction */ |
1588 | *phi_S = gmx_simd_atan2_rgmx_simd_atan2_f(cn_S, s_S); |
1589 | |
1590 | ipr_S = gmx_simd_iprod_rgmx_simd_iprod_f(rijx_S, rijy_S, rijz_S, |
1591 | *nx_S, *ny_S, *nz_S); |
1592 | |
1593 | iprm_S = gmx_simd_norm2_rgmx_simd_norm2_f(*mx_S, *my_S, *mz_S); |
1594 | iprn_S = gmx_simd_norm2_rgmx_simd_norm2_f(*nx_S, *ny_S, *nz_S); |
1595 | |
1596 | nrkj2_S = gmx_simd_norm2_rgmx_simd_norm2_f(rkjx_S, rkjy_S, rkjz_S); |
1597 | |
1598 | /* Avoid division by zero. When zero, the result is multiplied by 0 |
1599 | * anyhow, so the 3 max below do not affect the final result. |
1600 | */ |
1601 | nrkj2_S = gmx_simd_max_r_mm_max_ps(nrkj2_S, nrkj2_min_S); |
1602 | nrkj_1_S = gmx_simd_invsqrt_rgmx_simd_invsqrt_f(nrkj2_S); |
1603 | nrkj_2_S = gmx_simd_mul_r_mm_mul_ps(nrkj_1_S, nrkj_1_S); |
1604 | nrkj_S = gmx_simd_mul_r_mm_mul_ps(nrkj2_S, nrkj_1_S); |
1605 | |
1606 | toler_S = gmx_simd_mul_r_mm_mul_ps(nrkj2_S, real_eps_S); |
1607 | |
1608 | /* Here the plain-C code uses a conditional, but we can't do that in SIMD. |
1609 | * So we take a max with the tolerance instead. Since we multiply with |
1610 | * m or n later, the max does not affect the results. |
1611 | */ |
1612 | iprm_S = gmx_simd_max_r_mm_max_ps(iprm_S, toler_S); |
1613 | iprn_S = gmx_simd_max_r_mm_max_ps(iprn_S, toler_S); |
1614 | *nrkj_m2_S = gmx_simd_mul_r_mm_mul_ps(nrkj_S, gmx_simd_inv_rgmx_simd_inv_f(iprm_S)); |
1615 | *nrkj_n2_S = gmx_simd_mul_r_mm_mul_ps(nrkj_S, gmx_simd_inv_rgmx_simd_inv_f(iprn_S)); |
1616 | |
1617 | /* Set sign of phi_S with the sign of ipr_S; phi_S is currently positive */ |
1618 | *phi_S = gmx_simd_xor_sign_rgmx_simd_xor_sign_f(*phi_S, ipr_S); |
1619 | p_S = gmx_simd_iprod_rgmx_simd_iprod_f(rijx_S, rijy_S, rijz_S, |
1620 | rkjx_S, rkjy_S, rkjz_S); |
1621 | p_S = gmx_simd_mul_r_mm_mul_ps(p_S, nrkj_2_S); |
1622 | |
1623 | q_S = gmx_simd_iprod_rgmx_simd_iprod_f(rklx_S, rkly_S, rklz_S, |
1624 | rkjx_S, rkjy_S, rkjz_S); |
1625 | q_S = gmx_simd_mul_r_mm_mul_ps(q_S, nrkj_2_S); |
1626 | |
1627 | gmx_simd_store_r_mm_store_ps(p, p_S); |
1628 | gmx_simd_store_r_mm_store_ps(q, q_S); |
1629 | } |
1630 | |
1631 | #endif /* GMX_SIMD_HAVE_REAL */ |
1632 | |
1633 | |
1634 | void do_dih_fup(int i, int j, int k, int l, real ddphi, |
1635 | rvec r_ij, rvec r_kj, rvec r_kl, |
1636 | rvec m, rvec n, rvec f[], rvec fshift[], |
1637 | const t_pbc *pbc, const t_graph *g, |
1638 | const rvec x[], int t1, int t2, int t3) |
1639 | { |
1640 | /* 143 FLOPS */ |
1641 | rvec f_i, f_j, f_k, f_l; |
1642 | rvec uvec, vvec, svec, dx_jl; |
1643 | real iprm, iprn, nrkj, nrkj2, nrkj_1, nrkj_2; |
1644 | real a, b, p, q, toler; |
1645 | ivec jt, dt_ij, dt_kj, dt_lj; |
1646 | |
1647 | iprm = iprod(m, m); /* 5 */ |
1648 | iprn = iprod(n, n); /* 5 */ |
1649 | nrkj2 = iprod(r_kj, r_kj); /* 5 */ |
1650 | toler = nrkj2*GMX_REAL_EPS5.96046448E-08; |
1651 | if ((iprm > toler) && (iprn > toler)) |
1652 | { |
1653 | nrkj_1 = gmx_invsqrt(nrkj2)gmx_software_invsqrt(nrkj2); /* 10 */ |
1654 | nrkj_2 = nrkj_1*nrkj_1; /* 1 */ |
1655 | nrkj = nrkj2*nrkj_1; /* 1 */ |
1656 | a = -ddphi*nrkj/iprm; /* 11 */ |
1657 | svmul(a, m, f_i); /* 3 */ |
1658 | b = ddphi*nrkj/iprn; /* 11 */ |
1659 | svmul(b, n, f_l); /* 3 */ |
1660 | p = iprod(r_ij, r_kj); /* 5 */ |
1661 | p *= nrkj_2; /* 1 */ |
1662 | q = iprod(r_kl, r_kj); /* 5 */ |
1663 | q *= nrkj_2; /* 1 */ |
1664 | svmul(p, f_i, uvec); /* 3 */ |
1665 | svmul(q, f_l, vvec); /* 3 */ |
1666 | rvec_sub(uvec, vvec, svec); /* 3 */ |
1667 | rvec_sub(f_i, svec, f_j); /* 3 */ |
1668 | rvec_add(f_l, svec, f_k); /* 3 */ |
1669 | rvec_inc(f[i], f_i); /* 3 */ |
1670 | rvec_dec(f[j], f_j); /* 3 */ |
1671 | rvec_dec(f[k], f_k); /* 3 */ |
1672 | rvec_inc(f[l], f_l); /* 3 */ |
1673 | |
1674 | if (g) |
1675 | { |
1676 | copy_ivec(SHIFT_IVEC(g, j)((g)->ishift[j]), jt); |
1677 | ivec_sub(SHIFT_IVEC(g, i)((g)->ishift[i]), jt, dt_ij); |
1678 | ivec_sub(SHIFT_IVEC(g, k)((g)->ishift[k]), jt, dt_kj); |
1679 | ivec_sub(SHIFT_IVEC(g, l)((g)->ishift[l]), jt, dt_lj); |
1680 | t1 = IVEC2IS(dt_ij)(((2*2 +1)*((2*1 +1)*(((dt_ij)[2])+1)+((dt_ij)[1])+1)+((dt_ij )[0])+2)); |
1681 | t2 = IVEC2IS(dt_kj)(((2*2 +1)*((2*1 +1)*(((dt_kj)[2])+1)+((dt_kj)[1])+1)+((dt_kj )[0])+2)); |
1682 | t3 = IVEC2IS(dt_lj)(((2*2 +1)*((2*1 +1)*(((dt_lj)[2])+1)+((dt_lj)[1])+1)+((dt_lj )[0])+2)); |
1683 | } |
1684 | else if (pbc) |
1685 | { |
1686 | t3 = pbc_rvec_sub(pbc, x[l], x[j], dx_jl); |
1687 | } |
1688 | else |
1689 | { |
1690 | t3 = CENTRAL(((2*1 +1)*(2*1 +1)*(2*2 +1))/2); |
1691 | } |
1692 | |
1693 | rvec_inc(fshift[t1], f_i); |
1694 | rvec_dec(fshift[CENTRAL(((2*1 +1)*(2*1 +1)*(2*2 +1))/2)], f_j); |
1695 | rvec_dec(fshift[t2], f_k); |
1696 | rvec_inc(fshift[t3], f_l); |
1697 | } |
1698 | /* 112 TOTAL */ |
1699 | } |
1700 | |
1701 | /* As do_dih_fup above, but without shift forces */ |
1702 | static void |
1703 | do_dih_fup_noshiftf(int i, int j, int k, int l, real ddphi, |
1704 | rvec r_ij, rvec r_kj, rvec r_kl, |
1705 | rvec m, rvec n, rvec f[]) |
1706 | { |
1707 | rvec f_i, f_j, f_k, f_l; |
1708 | rvec uvec, vvec, svec, dx_jl; |
1709 | real iprm, iprn, nrkj, nrkj2, nrkj_1, nrkj_2; |
1710 | real a, b, p, q, toler; |
1711 | ivec jt, dt_ij, dt_kj, dt_lj; |
1712 | |
1713 | iprm = iprod(m, m); /* 5 */ |
1714 | iprn = iprod(n, n); /* 5 */ |
1715 | nrkj2 = iprod(r_kj, r_kj); /* 5 */ |
1716 | toler = nrkj2*GMX_REAL_EPS5.96046448E-08; |
1717 | if ((iprm > toler) && (iprn > toler)) |
1718 | { |
1719 | nrkj_1 = gmx_invsqrt(nrkj2)gmx_software_invsqrt(nrkj2); /* 10 */ |
1720 | nrkj_2 = nrkj_1*nrkj_1; /* 1 */ |
1721 | nrkj = nrkj2*nrkj_1; /* 1 */ |
1722 | a = -ddphi*nrkj/iprm; /* 11 */ |
1723 | svmul(a, m, f_i); /* 3 */ |
1724 | b = ddphi*nrkj/iprn; /* 11 */ |
1725 | svmul(b, n, f_l); /* 3 */ |
1726 | p = iprod(r_ij, r_kj); /* 5 */ |
1727 | p *= nrkj_2; /* 1 */ |
1728 | q = iprod(r_kl, r_kj); /* 5 */ |
1729 | q *= nrkj_2; /* 1 */ |
1730 | svmul(p, f_i, uvec); /* 3 */ |
1731 | svmul(q, f_l, vvec); /* 3 */ |
1732 | rvec_sub(uvec, vvec, svec); /* 3 */ |
1733 | rvec_sub(f_i, svec, f_j); /* 3 */ |
1734 | rvec_add(f_l, svec, f_k); /* 3 */ |
1735 | rvec_inc(f[i], f_i); /* 3 */ |
1736 | rvec_dec(f[j], f_j); /* 3 */ |
1737 | rvec_dec(f[k], f_k); /* 3 */ |
1738 | rvec_inc(f[l], f_l); /* 3 */ |
1739 | } |
1740 | } |
1741 | |
1742 | /* As do_dih_fup_noshiftf above, but with pre-calculated pre-factors */ |
1743 | static gmx_inlineinline void |
1744 | do_dih_fup_noshiftf_precalc(int i, int j, int k, int l, |
1745 | real p, real q, |
1746 | real f_i_x, real f_i_y, real f_i_z, |
1747 | real mf_l_x, real mf_l_y, real mf_l_z, |
1748 | rvec f[]) |
1749 | { |
1750 | rvec f_i, f_j, f_k, f_l; |
1751 | rvec uvec, vvec, svec; |
1752 | |
1753 | f_i[XX0] = f_i_x; |
1754 | f_i[YY1] = f_i_y; |
1755 | f_i[ZZ2] = f_i_z; |
1756 | f_l[XX0] = -mf_l_x; |
1757 | f_l[YY1] = -mf_l_y; |
1758 | f_l[ZZ2] = -mf_l_z; |
1759 | svmul(p, f_i, uvec); |
1760 | svmul(q, f_l, vvec); |
1761 | rvec_sub(uvec, vvec, svec); |
1762 | rvec_sub(f_i, svec, f_j); |
1763 | rvec_add(f_l, svec, f_k); |
1764 | rvec_inc(f[i], f_i); |
1765 | rvec_dec(f[j], f_j); |
1766 | rvec_dec(f[k], f_k); |
1767 | rvec_inc(f[l], f_l); |
1768 | } |
1769 | |
1770 | |
1771 | real dopdihs(real cpA, real cpB, real phiA, real phiB, int mult, |
1772 | real phi, real lambda, real *V, real *F) |
1773 | { |
1774 | real v, dvdlambda, mdphi, v1, sdphi, ddphi; |
1775 | real L1 = 1.0 - lambda; |
1776 | real ph0 = (L1*phiA + lambda*phiB)*DEG2RAD(3.14159265358979323846/180.0); |
1777 | real dph0 = (phiB - phiA)*DEG2RAD(3.14159265358979323846/180.0); |
1778 | real cp = L1*cpA + lambda*cpB; |
1779 | |
1780 | mdphi = mult*phi - ph0; |
1781 | sdphi = sin(mdphi); |
1782 | ddphi = -cp*mult*sdphi; |
1783 | v1 = 1.0 + cos(mdphi); |
1784 | v = cp*v1; |
1785 | |
1786 | dvdlambda = (cpB - cpA)*v1 + cp*dph0*sdphi; |
1787 | |
1788 | *V = v; |
1789 | *F = ddphi; |
1790 | |
1791 | return dvdlambda; |
1792 | |
1793 | /* That was 40 flops */ |
1794 | } |
1795 | |
1796 | static void |
1797 | dopdihs_noener(real cpA, real cpB, real phiA, real phiB, int mult, |
1798 | real phi, real lambda, real *F) |
1799 | { |
1800 | real mdphi, sdphi, ddphi; |
1801 | real L1 = 1.0 - lambda; |
1802 | real ph0 = (L1*phiA + lambda*phiB)*DEG2RAD(3.14159265358979323846/180.0); |
1803 | real cp = L1*cpA + lambda*cpB; |
1804 | |
1805 | mdphi = mult*phi - ph0; |
1806 | sdphi = sin(mdphi); |
1807 | ddphi = -cp*mult*sdphi; |
1808 | |
1809 | *F = ddphi; |
1810 | |
1811 | /* That was 20 flops */ |
1812 | } |
1813 | |
1814 | static void |
1815 | dopdihs_mdphi(real cpA, real cpB, real phiA, real phiB, int mult, |
1816 | real phi, real lambda, real *cp, real *mdphi) |
1817 | { |
1818 | real L1 = 1.0 - lambda; |
1819 | real ph0 = (L1*phiA + lambda*phiB)*DEG2RAD(3.14159265358979323846/180.0); |
1820 | |
1821 | *cp = L1*cpA + lambda*cpB; |
1822 | |
1823 | *mdphi = mult*phi - ph0; |
1824 | } |
1825 | |
1826 | static real dopdihs_min(real cpA, real cpB, real phiA, real phiB, int mult, |
1827 | real phi, real lambda, real *V, real *F) |
1828 | /* similar to dopdihs, except for a minus sign * |
1829 | * and a different treatment of mult/phi0 */ |
1830 | { |
1831 | real v, dvdlambda, mdphi, v1, sdphi, ddphi; |
1832 | real L1 = 1.0 - lambda; |
1833 | real ph0 = (L1*phiA + lambda*phiB)*DEG2RAD(3.14159265358979323846/180.0); |
1834 | real dph0 = (phiB - phiA)*DEG2RAD(3.14159265358979323846/180.0); |
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 | real pdihs(int nbonds, |
1854 | const t_iatom forceatoms[], const t_iparams forceparams[], |
1855 | const rvec x[], rvec f[], rvec fshift[], |
1856 | const t_pbc *pbc, const t_graph *g, |
1857 | real lambda, real *dvdlambda, |
1858 | const t_mdatoms gmx_unused__attribute__ ((unused)) *md, t_fcdata gmx_unused__attribute__ ((unused)) *fcd, |
1859 | int gmx_unused__attribute__ ((unused)) *global_atom_index) |
1860 | { |
1861 | int i, type, ai, aj, ak, al; |
1862 | int t1, t2, t3; |
1863 | rvec r_ij, r_kj, r_kl, m, n; |
1864 | real phi, sign, ddphi, vpd, vtot; |
1865 | |
1866 | vtot = 0.0; |
1867 | |
1868 | for (i = 0; (i < nbonds); ) |
1869 | { |
1870 | type = forceatoms[i++]; |
1871 | ai = forceatoms[i++]; |
1872 | aj = forceatoms[i++]; |
1873 | ak = forceatoms[i++]; |
1874 | al = forceatoms[i++]; |
1875 | |
1876 | phi = dih_angle(x[ai], x[aj], x[ak], x[al], pbc, r_ij, r_kj, r_kl, m, n, |
1877 | &sign, &t1, &t2, &t3); /* 84 */ |
1878 | *dvdlambda += dopdihs(forceparams[type].pdihs.cpA, |
1879 | forceparams[type].pdihs.cpB, |
1880 | forceparams[type].pdihs.phiA, |
1881 | forceparams[type].pdihs.phiB, |
1882 | forceparams[type].pdihs.mult, |
1883 | phi, lambda, &vpd, &ddphi); |
1884 | |
1885 | vtot += vpd; |
1886 | do_dih_fup(ai, aj, ak, al, ddphi, r_ij, r_kj, r_kl, m, n, |
1887 | f, fshift, pbc, g, x, t1, t2, t3); /* 112 */ |
1888 | |
1889 | #ifdef DEBUG |
1890 | fprintf(debug, "pdih: (%d,%d,%d,%d) phi=%g\n", |
1891 | ai, aj, ak, al, phi); |
1892 | #endif |
1893 | } /* 223 TOTAL */ |
1894 | |
1895 | return vtot; |
1896 | } |
1897 | |
1898 | void make_dp_periodic(real *dp) /* 1 flop? */ |
1899 | { |
1900 | /* dp cannot be outside (-pi,pi) */ |
1901 | if (*dp >= M_PI3.14159265358979323846) |
1902 | { |
1903 | *dp -= 2*M_PI3.14159265358979323846; |
1904 | } |
1905 | else if (*dp < -M_PI3.14159265358979323846) |
1906 | { |
1907 | *dp += 2*M_PI3.14159265358979323846; |
1908 | } |
1909 | return; |
1910 | } |
1911 | |
1912 | /* As pdihs above, but without calculating energies and shift forces */ |
1913 | static void |
1914 | pdihs_noener(int nbonds, |
1915 | const t_iatom forceatoms[], const t_iparams forceparams[], |
1916 | const rvec x[], rvec f[], |
1917 | const t_pbc gmx_unused__attribute__ ((unused)) *pbc, const t_graph gmx_unused__attribute__ ((unused)) *g, |
1918 | real lambda, |
1919 | const t_mdatoms gmx_unused__attribute__ ((unused)) *md, t_fcdata gmx_unused__attribute__ ((unused)) *fcd, |
1920 | int gmx_unused__attribute__ ((unused)) *global_atom_index) |
1921 | { |
1922 | int i, type, ai, aj, ak, al; |
1923 | int t1, t2, t3; |
1924 | rvec r_ij, r_kj, r_kl, m, n; |
1925 | real phi, sign, ddphi_tot, ddphi; |
1926 | |
1927 | for (i = 0; (i < nbonds); ) |
1928 | { |
1929 | ai = forceatoms[i+1]; |
1930 | aj = forceatoms[i+2]; |
1931 | ak = forceatoms[i+3]; |
1932 | al = forceatoms[i+4]; |
1933 | |
1934 | phi = dih_angle(x[ai], x[aj], x[ak], x[al], pbc, r_ij, r_kj, r_kl, m, n, |
1935 | &sign, &t1, &t2, &t3); |
1936 | |
1937 | ddphi_tot = 0; |
1938 | |
1939 | /* Loop over dihedrals working on the same atoms, |
1940 | * so we avoid recalculating angles and force distributions. |
1941 | */ |
1942 | do |
1943 | { |
1944 | type = forceatoms[i]; |
1945 | dopdihs_noener(forceparams[type].pdihs.cpA, |
1946 | forceparams[type].pdihs.cpB, |
1947 | forceparams[type].pdihs.phiA, |
1948 | forceparams[type].pdihs.phiB, |
1949 | forceparams[type].pdihs.mult, |
1950 | phi, lambda, &ddphi); |
1951 | ddphi_tot += ddphi; |
1952 | |
1953 | i += 5; |
1954 | } |
1955 | while (i < nbonds && |
1956 | forceatoms[i+1] == ai && |
1957 | forceatoms[i+2] == aj && |
1958 | forceatoms[i+3] == ak && |
1959 | forceatoms[i+4] == al); |
1960 | |
1961 | do_dih_fup_noshiftf(ai, aj, ak, al, ddphi_tot, r_ij, r_kj, r_kl, m, n, f); |
1962 | } |
1963 | } |
1964 | |
1965 | |
1966 | #ifdef GMX_SIMD_HAVE_REAL |
1967 | |
1968 | /* As pdihs_noner above, but using SIMD to calculate many dihedrals at once */ |
1969 | static void |
1970 | pdihs_noener_simd(int nbonds, |
1971 | const t_iatom forceatoms[], const t_iparams forceparams[], |
1972 | const rvec x[], rvec f[], |
1973 | const t_pbc *pbc, const t_graph gmx_unused__attribute__ ((unused)) *g, |
1974 | real gmx_unused__attribute__ ((unused)) lambda, |
1975 | const t_mdatoms gmx_unused__attribute__ ((unused)) *md, t_fcdata gmx_unused__attribute__ ((unused)) *fcd, |
1976 | int gmx_unused__attribute__ ((unused)) *global_atom_index) |
1977 | { |
1978 | const int nfa1 = 5; |
1979 | int i, iu, s; |
1980 | int type, ai[GMX_SIMD_REAL_WIDTH4], aj[GMX_SIMD_REAL_WIDTH4], ak[GMX_SIMD_REAL_WIDTH4], al[GMX_SIMD_REAL_WIDTH4]; |
1981 | int t1[GMX_SIMD_REAL_WIDTH4], t2[GMX_SIMD_REAL_WIDTH4], t3[GMX_SIMD_REAL_WIDTH4]; |
1982 | real ddphi; |
1983 | real dr_array[3*DIM3*GMX_SIMD_REAL_WIDTH4+GMX_SIMD_REAL_WIDTH4], *dr; |
1984 | real buf_array[7*GMX_SIMD_REAL_WIDTH4+GMX_SIMD_REAL_WIDTH4], *buf; |
1985 | real *cp, *phi0, *mult, *phi, *p, *q, *sf_i, *msf_l; |
1986 | gmx_simd_real_t__m128 phi0_S, phi_S; |
1987 | gmx_simd_real_t__m128 mx_S, my_S, mz_S; |
1988 | gmx_simd_real_t__m128 nx_S, ny_S, nz_S; |
1989 | gmx_simd_real_t__m128 nrkj_m2_S, nrkj_n2_S; |
1990 | gmx_simd_real_t__m128 cp_S, mdphi_S, mult_S; |
1991 | gmx_simd_real_t__m128 sin_S, cos_S; |
1992 | gmx_simd_real_t__m128 mddphi_S; |
1993 | gmx_simd_real_t__m128 sf_i_S, msf_l_S; |
1994 | pbc_simd_t pbc_simd; |
1995 | |
1996 | /* Ensure SIMD register alignment */ |
1997 | dr = gmx_simd_align_rgmx_simd_align_f(dr_array); |
1998 | buf = gmx_simd_align_rgmx_simd_align_f(buf_array); |
1999 | |
2000 | /* Extract aligned pointer for parameters and variables */ |
2001 | cp = buf + 0*GMX_SIMD_REAL_WIDTH4; |
2002 | phi0 = buf + 1*GMX_SIMD_REAL_WIDTH4; |
2003 | mult = buf + 2*GMX_SIMD_REAL_WIDTH4; |
2004 | p = buf + 3*GMX_SIMD_REAL_WIDTH4; |
2005 | q = buf + 4*GMX_SIMD_REAL_WIDTH4; |
2006 | sf_i = buf + 5*GMX_SIMD_REAL_WIDTH4; |
2007 | msf_l = buf + 6*GMX_SIMD_REAL_WIDTH4; |
2008 | |
2009 | set_pbc_simd(pbc, &pbc_simd); |
2010 | |
2011 | /* nbonds is the number of dihedrals times nfa1, here we step GMX_SIMD_REAL_WIDTH dihs */ |
2012 | for (i = 0; (i < nbonds); i += GMX_SIMD_REAL_WIDTH4*nfa1) |
2013 | { |
2014 | /* Collect atoms quadruplets for GMX_SIMD_REAL_WIDTH dihedrals. |
2015 | * iu indexes into forceatoms, we should not let iu go beyond nbonds. |
2016 | */ |
2017 | iu = i; |
2018 | for (s = 0; s < GMX_SIMD_REAL_WIDTH4; s++) |
2019 | { |
2020 | type = forceatoms[iu]; |
2021 | ai[s] = forceatoms[iu+1]; |
2022 | aj[s] = forceatoms[iu+2]; |
2023 | ak[s] = forceatoms[iu+3]; |
2024 | al[s] = forceatoms[iu+4]; |
2025 | |
2026 | cp[s] = forceparams[type].pdihs.cpA; |
2027 | phi0[s] = forceparams[type].pdihs.phiA*DEG2RAD(3.14159265358979323846/180.0); |
2028 | mult[s] = forceparams[type].pdihs.mult; |
2029 | |
2030 | /* At the end fill the arrays with identical entries */ |
2031 | if (iu + nfa1 < nbonds) |
2032 | { |
2033 | iu += nfa1; |
2034 | } |
2035 | } |
2036 | |
2037 | /* Caclulate GMX_SIMD_REAL_WIDTH dihedral angles at once */ |
2038 | dih_angle_simd(x, ai, aj, ak, al, &pbc_simd, |
2039 | dr, |
2040 | &phi_S, |
2041 | &mx_S, &my_S, &mz_S, |
2042 | &nx_S, &ny_S, &nz_S, |
2043 | &nrkj_m2_S, |
2044 | &nrkj_n2_S, |
2045 | p, q); |
2046 | |
2047 | cp_S = gmx_simd_load_r_mm_load_ps(cp); |
2048 | phi0_S = gmx_simd_load_r_mm_load_ps(phi0); |
2049 | mult_S = gmx_simd_load_r_mm_load_ps(mult); |
2050 | |
2051 | mdphi_S = gmx_simd_sub_r_mm_sub_ps(gmx_simd_mul_r_mm_mul_ps(mult_S, phi_S), phi0_S); |
2052 | |
2053 | /* Calculate GMX_SIMD_REAL_WIDTH sines at once */ |
2054 | gmx_simd_sincos_rgmx_simd_sincos_f(mdphi_S, &sin_S, &cos_S); |
2055 | mddphi_S = gmx_simd_mul_r_mm_mul_ps(gmx_simd_mul_r_mm_mul_ps(cp_S, mult_S), sin_S); |
2056 | sf_i_S = gmx_simd_mul_r_mm_mul_ps(mddphi_S, nrkj_m2_S); |
2057 | msf_l_S = gmx_simd_mul_r_mm_mul_ps(mddphi_S, nrkj_n2_S); |
2058 | |
2059 | /* After this m?_S will contain f[i] */ |
2060 | mx_S = gmx_simd_mul_r_mm_mul_ps(sf_i_S, mx_S); |
2061 | my_S = gmx_simd_mul_r_mm_mul_ps(sf_i_S, my_S); |
2062 | mz_S = gmx_simd_mul_r_mm_mul_ps(sf_i_S, mz_S); |
2063 | |
2064 | /* After this m?_S will contain -f[l] */ |
2065 | nx_S = gmx_simd_mul_r_mm_mul_ps(msf_l_S, nx_S); |
2066 | ny_S = gmx_simd_mul_r_mm_mul_ps(msf_l_S, ny_S); |
2067 | nz_S = gmx_simd_mul_r_mm_mul_ps(msf_l_S, nz_S); |
2068 | |
2069 | gmx_simd_store_r_mm_store_ps(dr + 0*GMX_SIMD_REAL_WIDTH4, mx_S); |
2070 | gmx_simd_store_r_mm_store_ps(dr + 1*GMX_SIMD_REAL_WIDTH4, my_S); |
2071 | gmx_simd_store_r_mm_store_ps(dr + 2*GMX_SIMD_REAL_WIDTH4, mz_S); |
2072 | gmx_simd_store_r_mm_store_ps(dr + 3*GMX_SIMD_REAL_WIDTH4, nx_S); |
2073 | gmx_simd_store_r_mm_store_ps(dr + 4*GMX_SIMD_REAL_WIDTH4, ny_S); |
2074 | gmx_simd_store_r_mm_store_ps(dr + 5*GMX_SIMD_REAL_WIDTH4, nz_S); |
2075 | |
2076 | iu = i; |
2077 | s = 0; |
2078 | do |
2079 | { |
2080 | do_dih_fup_noshiftf_precalc(ai[s], aj[s], ak[s], al[s], |
2081 | p[s], q[s], |
2082 | dr[ XX0 *GMX_SIMD_REAL_WIDTH4+s], |
2083 | dr[ YY1 *GMX_SIMD_REAL_WIDTH4+s], |
2084 | dr[ ZZ2 *GMX_SIMD_REAL_WIDTH4+s], |
2085 | dr[(DIM3+XX0)*GMX_SIMD_REAL_WIDTH4+s], |
2086 | dr[(DIM3+YY1)*GMX_SIMD_REAL_WIDTH4+s], |
2087 | dr[(DIM3+ZZ2)*GMX_SIMD_REAL_WIDTH4+s], |
2088 | f); |
2089 | s++; |
2090 | iu += nfa1; |
2091 | } |
2092 | while (s < GMX_SIMD_REAL_WIDTH4 && iu < nbonds); |
2093 | } |
2094 | } |
2095 | |
2096 | #endif /* GMX_SIMD_HAVE_REAL */ |
2097 | |
2098 | |
2099 | real idihs(int nbonds, |
2100 | const t_iatom forceatoms[], const t_iparams forceparams[], |
2101 | const rvec x[], rvec f[], rvec fshift[], |
2102 | const t_pbc *pbc, const t_graph *g, |
2103 | real lambda, real *dvdlambda, |
2104 | const t_mdatoms gmx_unused__attribute__ ((unused)) *md, t_fcdata gmx_unused__attribute__ ((unused)) *fcd, |
2105 | int gmx_unused__attribute__ ((unused)) *global_atom_index) |
2106 | { |
2107 | int i, type, ai, aj, ak, al; |
2108 | int t1, t2, t3; |
2109 | real phi, phi0, dphi0, ddphi, sign, vtot; |
2110 | rvec r_ij, r_kj, r_kl, m, n; |
2111 | real L1, kk, dp, dp2, kA, kB, pA, pB, dvdl_term; |
2112 | |
2113 | L1 = 1.0-lambda; |
2114 | dvdl_term = 0; |
2115 | vtot = 0.0; |
2116 | for (i = 0; (i < nbonds); ) |
2117 | { |
2118 | type = forceatoms[i++]; |
2119 | ai = forceatoms[i++]; |
2120 | aj = forceatoms[i++]; |
2121 | ak = forceatoms[i++]; |
2122 | al = forceatoms[i++]; |
2123 | |
2124 | phi = dih_angle(x[ai], x[aj], x[ak], x[al], pbc, r_ij, r_kj, r_kl, m, n, |
2125 | &sign, &t1, &t2, &t3); /* 84 */ |
2126 | |
2127 | /* phi can jump if phi0 is close to Pi/-Pi, which will cause huge |
2128 | * force changes if we just apply a normal harmonic. |
2129 | * Instead, we first calculate phi-phi0 and take it modulo (-Pi,Pi). |
2130 | * This means we will never have the periodicity problem, unless |
2131 | * the dihedral is Pi away from phiO, which is very unlikely due to |
2132 | * the potential. |
2133 | */ |
2134 | kA = forceparams[type].harmonic.krA; |
2135 | kB = forceparams[type].harmonic.krB; |
2136 | pA = forceparams[type].harmonic.rA; |
2137 | pB = forceparams[type].harmonic.rB; |
2138 | |
2139 | kk = L1*kA + lambda*kB; |
2140 | phi0 = (L1*pA + lambda*pB)*DEG2RAD(3.14159265358979323846/180.0); |
2141 | dphi0 = (pB - pA)*DEG2RAD(3.14159265358979323846/180.0); |
2142 | |
2143 | dp = phi-phi0; |
2144 | |
2145 | make_dp_periodic(&dp); |
2146 | |
2147 | dp2 = dp*dp; |
2148 | |
2149 | vtot += 0.5*kk*dp2; |
2150 | ddphi = -kk*dp; |
2151 | |
2152 | dvdl_term += 0.5*(kB - kA)*dp2 - kk*dphi0*dp; |
2153 | |
2154 | do_dih_fup(ai, aj, ak, al, (real)(-ddphi), r_ij, r_kj, r_kl, m, n, |
2155 | f, fshift, pbc, g, x, t1, t2, t3); /* 112 */ |
2156 | /* 218 TOTAL */ |
2157 | #ifdef DEBUG |
2158 | if (debug) |
2159 | { |
2160 | fprintf(debug, "idih: (%d,%d,%d,%d) phi=%g\n", |
2161 | ai, aj, ak, al, phi); |
2162 | } |
2163 | #endif |
2164 | } |
2165 | |
2166 | *dvdlambda += dvdl_term; |
2167 | return vtot; |
2168 | } |
2169 | |
2170 | |
2171 | /*! \brief returns dx, rdist, and dpdl for functions posres() and fbposres() |
2172 | */ |
2173 | static void posres_dx(const rvec x, const rvec pos0A, const rvec pos0B, |
2174 | const rvec comA_sc, const rvec comB_sc, |
2175 | real lambda, |
2176 | t_pbc *pbc, int refcoord_scaling, int npbcdim, |
2177 | rvec dx, rvec rdist, rvec dpdl) |
2178 | { |
2179 | int m, d; |
2180 | real posA, posB, L1, ref = 0.; |
2181 | rvec pos; |
2182 | |
2183 | L1 = 1.0-lambda; |
2184 | |
2185 | for (m = 0; m < DIM3; m++) |
2186 | { |
2187 | posA = pos0A[m]; |
2188 | posB = pos0B[m]; |
2189 | if (m < npbcdim) |
2190 | { |
2191 | switch (refcoord_scaling) |
2192 | { |
2193 | case erscNO: |
2194 | ref = 0; |
2195 | rdist[m] = L1*posA + lambda*posB; |
2196 | dpdl[m] = posB - posA; |
2197 | break; |
2198 | case erscALL: |
2199 | /* Box relative coordinates are stored for dimensions with pbc */ |
2200 | posA *= pbc->box[m][m]; |
2201 | posB *= pbc->box[m][m]; |
2202 | assert(npbcdim <= DIM)((void) (0)); |
2203 | for (d = m+1; d < npbcdim; d++) |
2204 | { |
2205 | posA += pos0A[d]*pbc->box[d][m]; |
2206 | posB += pos0B[d]*pbc->box[d][m]; |
2207 | } |
2208 | ref = L1*posA + lambda*posB; |
2209 | rdist[m] = 0; |
2210 | dpdl[m] = posB - posA; |
2211 | break; |
2212 | case erscCOM: |
2213 | ref = L1*comA_sc[m] + lambda*comB_sc[m]; |
2214 | rdist[m] = L1*posA + lambda*posB; |
2215 | dpdl[m] = comB_sc[m] - comA_sc[m] + posB - posA; |
2216 | break; |
2217 | default: |
2218 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxlib/bondfree.c" , 2218, "No such scaling method implemented"); |
2219 | } |
2220 | } |
2221 | else |
2222 | { |
2223 | ref = L1*posA + lambda*posB; |
2224 | rdist[m] = 0; |
2225 | dpdl[m] = posB - posA; |
2226 | } |
2227 | |
2228 | /* We do pbc_dx with ref+rdist, |
2229 | * since with only ref we can be up to half a box vector wrong. |
2230 | */ |
2231 | pos[m] = ref + rdist[m]; |
2232 | } |
2233 | |
2234 | if (pbc) |
2235 | { |
2236 | pbc_dx(pbc, x, pos, dx); |
2237 | } |
2238 | else |
2239 | { |
2240 | rvec_sub(x, pos, dx); |
2241 | } |
2242 | } |
2243 | |
2244 | /*! \brief Adds forces of flat-bottomed positions restraints to f[] |
2245 | * and fixes vir_diag. Returns the flat-bottomed potential. */ |
2246 | real fbposres(int nbonds, |
2247 | const t_iatom forceatoms[], const t_iparams forceparams[], |
2248 | const rvec x[], rvec f[], rvec vir_diag, |
2249 | t_pbc *pbc, |
2250 | int refcoord_scaling, int ePBC, rvec com) |
2251 | /* compute flat-bottomed positions restraints */ |
2252 | { |
2253 | int i, ai, m, d, type, npbcdim = 0, fbdim; |
2254 | const t_iparams *pr; |
2255 | real vtot, kk, v; |
2256 | real ref = 0, dr, dr2, rpot, rfb, rfb2, fact, invdr; |
2257 | rvec com_sc, rdist, pos, dx, dpdl, fm; |
2258 | gmx_bool bInvert; |
2259 | |
2260 | npbcdim = ePBC2npbcdim(ePBC); |
2261 | |
2262 | if (refcoord_scaling == erscCOM) |
2263 | { |
2264 | clear_rvec(com_sc); |
2265 | for (m = 0; m < npbcdim; m++) |
2266 | { |
2267 | assert(npbcdim <= DIM)((void) (0)); |
2268 | for (d = m; d < npbcdim; d++) |
2269 | { |
2270 | com_sc[m] += com[d]*pbc->box[d][m]; |
2271 | } |
2272 | } |
2273 | } |
2274 | |
2275 | vtot = 0.0; |
2276 | for (i = 0; (i < nbonds); ) |
2277 | { |
2278 | type = forceatoms[i++]; |
2279 | ai = forceatoms[i++]; |
2280 | pr = &forceparams[type]; |
2281 | |
2282 | /* same calculation as for normal posres, but with identical A and B states, and lambda==0 */ |
2283 | posres_dx(x[ai], forceparams[type].fbposres.pos0, forceparams[type].fbposres.pos0, |
2284 | com_sc, com_sc, 0.0, |
2285 | pbc, refcoord_scaling, npbcdim, |
2286 | dx, rdist, dpdl); |
2287 | |
2288 | clear_rvec(fm); |
2289 | v = 0.0; |
2290 | |
2291 | kk = pr->fbposres.k; |
2292 | rfb = pr->fbposres.r; |
2293 | rfb2 = sqr(rfb); |
2294 | |
2295 | /* with rfb<0, push particle out of the sphere/cylinder/layer */ |
2296 | bInvert = FALSE0; |
2297 | if (rfb < 0.) |
2298 | { |
2299 | bInvert = TRUE1; |
2300 | rfb = -rfb; |
2301 | } |
2302 | |
2303 | switch (pr->fbposres.geom) |
2304 | { |
2305 | case efbposresSPHERE: |
2306 | /* spherical flat-bottom posres */ |
2307 | dr2 = norm2(dx); |
2308 | if (dr2 > 0.0 && |
2309 | ( (dr2 > rfb2 && bInvert == FALSE0 ) || (dr2 < rfb2 && bInvert == TRUE1 ) ) |
2310 | ) |
2311 | { |
2312 | dr = sqrt(dr2); |
2313 | v = 0.5*kk*sqr(dr - rfb); |
2314 | fact = -kk*(dr-rfb)/dr; /* Force pointing to the center pos0 */ |
2315 | svmul(fact, dx, fm); |
2316 | } |
2317 | break; |
2318 | case efbposresCYLINDER: |
2319 | /* cylidrical flat-bottom posres in x-y plane. fm[ZZ] = 0. */ |
2320 | dr2 = sqr(dx[XX0])+sqr(dx[YY1]); |
2321 | if (dr2 > 0.0 && |
2322 | ( (dr2 > rfb2 && bInvert == FALSE0 ) || (dr2 < rfb2 && bInvert == TRUE1 ) ) |
2323 | ) |
2324 | { |
2325 | dr = sqrt(dr2); |
2326 | invdr = 1./dr; |
2327 | v = 0.5*kk*sqr(dr - rfb); |
2328 | fm[XX0] = -kk*(dr-rfb)*dx[XX0]*invdr; /* Force pointing to the center */ |
2329 | fm[YY1] = -kk*(dr-rfb)*dx[YY1]*invdr; |
2330 | } |
2331 | break; |
2332 | case efbposresX: /* fbdim=XX */ |
2333 | case efbposresY: /* fbdim=YY */ |
2334 | case efbposresZ: /* fbdim=ZZ */ |
2335 | /* 1D flat-bottom potential */ |
2336 | fbdim = pr->fbposres.geom - efbposresX; |
2337 | dr = dx[fbdim]; |
2338 | if ( ( dr > rfb && bInvert == FALSE0 ) || ( 0 < dr && dr < rfb && bInvert == TRUE1 ) ) |
2339 | { |
2340 | v = 0.5*kk*sqr(dr - rfb); |
2341 | fm[fbdim] = -kk*(dr - rfb); |
2342 | } |
2343 | else if ( (dr < (-rfb) && bInvert == FALSE0 ) || ( (-rfb) < dr && dr < 0 && bInvert == TRUE1 )) |
2344 | { |
2345 | v = 0.5*kk*sqr(dr + rfb); |
2346 | fm[fbdim] = -kk*(dr + rfb); |
2347 | } |
2348 | break; |
2349 | } |
2350 | |
2351 | vtot += v; |
2352 | |
2353 | for (m = 0; (m < DIM3); m++) |
2354 | { |
2355 | f[ai][m] += fm[m]; |
2356 | /* Here we correct for the pbc_dx which included rdist */ |
2357 | vir_diag[m] -= 0.5*(dx[m] + rdist[m])*fm[m]; |
2358 | } |
2359 | } |
2360 | |
2361 | return vtot; |
2362 | } |
2363 | |
2364 | |
2365 | real posres(int nbonds, |
2366 | const t_iatom forceatoms[], const t_iparams forceparams[], |
2367 | const rvec x[], rvec f[], rvec vir_diag, |
2368 | t_pbc *pbc, |
2369 | real lambda, real *dvdlambda, |
2370 | int refcoord_scaling, int ePBC, rvec comA, rvec comB) |
2371 | { |
2372 | int i, ai, m, d, type, ki, npbcdim = 0; |
2373 | const t_iparams *pr; |
2374 | real L1; |
2375 | real vtot, kk, fm; |
2376 | real posA, posB, ref = 0; |
2377 | rvec comA_sc, comB_sc, rdist, dpdl, pos, dx; |
2378 | gmx_bool bForceValid = TRUE1; |
2379 | |
2380 | if ((f == NULL((void*)0)) || (vir_diag == NULL((void*)0))) /* should both be null together! */ |
2381 | { |
2382 | bForceValid = FALSE0; |
2383 | } |
2384 | |
2385 | npbcdim = ePBC2npbcdim(ePBC); |
2386 | |
2387 | if (refcoord_scaling == erscCOM) |
2388 | { |
2389 | clear_rvec(comA_sc); |
2390 | clear_rvec(comB_sc); |
2391 | for (m = 0; m < npbcdim; m++) |
2392 | { |
2393 | assert(npbcdim <= DIM)((void) (0)); |
2394 | for (d = m; d < npbcdim; d++) |
2395 | { |
2396 | comA_sc[m] += comA[d]*pbc->box[d][m]; |
2397 | comB_sc[m] += comB[d]*pbc->box[d][m]; |
2398 | } |
2399 | } |
2400 | } |
2401 | |
2402 | L1 = 1.0 - lambda; |
2403 | |
2404 | vtot = 0.0; |
2405 | for (i = 0; (i < nbonds); ) |
2406 | { |
2407 | type = forceatoms[i++]; |
2408 | ai = forceatoms[i++]; |
2409 | pr = &forceparams[type]; |
2410 | |
2411 | /* return dx, rdist, and dpdl */ |
2412 | posres_dx(x[ai], forceparams[type].posres.pos0A, forceparams[type].posres.pos0B, |
2413 | comA_sc, comB_sc, lambda, |
2414 | pbc, refcoord_scaling, npbcdim, |
2415 | dx, rdist, dpdl); |
2416 | |
2417 | for (m = 0; (m < DIM3); m++) |
2418 | { |
2419 | kk = L1*pr->posres.fcA[m] + lambda*pr->posres.fcB[m]; |
2420 | fm = -kk*dx[m]; |
2421 | vtot += 0.5*kk*dx[m]*dx[m]; |
2422 | *dvdlambda += |
2423 | 0.5*(pr->posres.fcB[m] - pr->posres.fcA[m])*dx[m]*dx[m] |
2424 | -fm*dpdl[m]; |
2425 | |
2426 | /* Here we correct for the pbc_dx which included rdist */ |
2427 | if (bForceValid) |
2428 | { |
2429 | f[ai][m] += fm; |
2430 | vir_diag[m] -= 0.5*(dx[m] + rdist[m])*fm; |
2431 | } |
2432 | } |
2433 | } |
2434 | |
2435 | return vtot; |
2436 | } |
2437 | |
2438 | static real low_angres(int nbonds, |
2439 | const t_iatom forceatoms[], const t_iparams forceparams[], |
2440 | const rvec x[], rvec f[], rvec fshift[], |
2441 | const t_pbc *pbc, const t_graph *g, |
2442 | real lambda, real *dvdlambda, |
2443 | gmx_bool bZAxis) |
2444 | { |
2445 | int i, m, type, ai, aj, ak, al; |
2446 | int t1, t2; |
2447 | real phi, cos_phi, cos_phi2, vid, vtot, dVdphi; |
2448 | rvec r_ij, r_kl, f_i, f_k = {0, 0, 0}; |
2449 | real st, sth, nrij2, nrkl2, c, cij, ckl; |
2450 | |
2451 | ivec dt; |
2452 | t2 = 0; /* avoid warning with gcc-3.3. It is never used uninitialized */ |
2453 | |
2454 | vtot = 0.0; |
2455 | ak = al = 0; /* to avoid warnings */ |
2456 | for (i = 0; i < nbonds; ) |
2457 | { |
2458 | type = forceatoms[i++]; |
2459 | ai = forceatoms[i++]; |
2460 | aj = forceatoms[i++]; |
2461 | t1 = pbc_rvec_sub(pbc, x[aj], x[ai], r_ij); /* 3 */ |
2462 | if (!bZAxis) |
2463 | { |
2464 | ak = forceatoms[i++]; |
2465 | al = forceatoms[i++]; |
2466 | t2 = pbc_rvec_sub(pbc, x[al], x[ak], r_kl); /* 3 */ |
2467 | } |
2468 | else |
2469 | { |
2470 | r_kl[XX0] = 0; |
2471 | r_kl[YY1] = 0; |
2472 | r_kl[ZZ2] = 1; |
2473 | } |
2474 | |
2475 | cos_phi = cos_angle(r_ij, r_kl); /* 25 */ |
2476 | phi = acos(cos_phi); /* 10 */ |
2477 | |
2478 | *dvdlambda += dopdihs_min(forceparams[type].pdihs.cpA, |
2479 | forceparams[type].pdihs.cpB, |
2480 | forceparams[type].pdihs.phiA, |
2481 | forceparams[type].pdihs.phiB, |
2482 | forceparams[type].pdihs.mult, |
2483 | phi, lambda, &vid, &dVdphi); /* 40 */ |
2484 | |
2485 | vtot += vid; |
2486 | |
2487 | cos_phi2 = sqr(cos_phi); /* 1 */ |
2488 | if (cos_phi2 < 1) |
2489 | { |
2490 | st = -dVdphi*gmx_invsqrt(1 - cos_phi2)gmx_software_invsqrt(1 - cos_phi2); /* 12 */ |
2491 | sth = st*cos_phi; /* 1 */ |
2492 | nrij2 = iprod(r_ij, r_ij); /* 5 */ |
2493 | nrkl2 = iprod(r_kl, r_kl); /* 5 */ |
2494 | |
2495 | c = st*gmx_invsqrt(nrij2*nrkl2)gmx_software_invsqrt(nrij2*nrkl2); /* 11 */ |
2496 | cij = sth/nrij2; /* 10 */ |
2497 | ckl = sth/nrkl2; /* 10 */ |
2498 | |
2499 | for (m = 0; m < DIM3; m++) /* 18+18 */ |
2500 | { |
2501 | f_i[m] = (c*r_kl[m]-cij*r_ij[m]); |
2502 | f[ai][m] += f_i[m]; |
2503 | f[aj][m] -= f_i[m]; |
2504 | if (!bZAxis) |
2505 | { |
2506 | f_k[m] = (c*r_ij[m]-ckl*r_kl[m]); |
2507 | f[ak][m] += f_k[m]; |
2508 | f[al][m] -= f_k[m]; |
2509 | } |
2510 | } |
2511 | |
2512 | if (g) |
2513 | { |
2514 | ivec_sub(SHIFT_IVEC(g, ai)((g)->ishift[ai]), SHIFT_IVEC(g, aj)((g)->ishift[aj]), dt); |
2515 | t1 = IVEC2IS(dt)(((2*2 +1)*((2*1 +1)*(((dt)[2])+1)+((dt)[1])+1)+((dt)[0])+2)); |
2516 | } |
2517 | rvec_inc(fshift[t1], f_i); |
2518 | rvec_dec(fshift[CENTRAL(((2*1 +1)*(2*1 +1)*(2*2 +1))/2)], f_i); |
2519 | if (!bZAxis) |
2520 | { |
2521 | if (g) |
2522 | { |
2523 | ivec_sub(SHIFT_IVEC(g, ak)((g)->ishift[ak]), SHIFT_IVEC(g, al)((g)->ishift[al]), dt); |
2524 | t2 = IVEC2IS(dt)(((2*2 +1)*((2*1 +1)*(((dt)[2])+1)+((dt)[1])+1)+((dt)[0])+2)); |
2525 | } |
2526 | rvec_inc(fshift[t2], f_k); |
2527 | rvec_dec(fshift[CENTRAL(((2*1 +1)*(2*1 +1)*(2*2 +1))/2)], f_k); |
2528 | } |
2529 | } |
2530 | } |
2531 | |
2532 | return vtot; /* 184 / 157 (bZAxis) total */ |
2533 | } |
2534 | |
2535 | real angres(int nbonds, |
2536 | const t_iatom forceatoms[], const t_iparams forceparams[], |
2537 | const rvec x[], rvec f[], rvec fshift[], |
2538 | const t_pbc *pbc, const t_graph *g, |
2539 | real lambda, real *dvdlambda, |
2540 | const t_mdatoms gmx_unused__attribute__ ((unused)) *md, t_fcdata gmx_unused__attribute__ ((unused)) *fcd, |
2541 | int gmx_unused__attribute__ ((unused)) *global_atom_index) |
2542 | { |
2543 | return low_angres(nbonds, forceatoms, forceparams, x, f, fshift, pbc, g, |
2544 | lambda, dvdlambda, FALSE0); |
2545 | } |
2546 | |
2547 | real angresz(int nbonds, |
2548 | const t_iatom forceatoms[], const t_iparams forceparams[], |
2549 | const rvec x[], rvec f[], rvec fshift[], |
2550 | const t_pbc *pbc, const t_graph *g, |
2551 | real lambda, real *dvdlambda, |
2552 | const t_mdatoms gmx_unused__attribute__ ((unused)) *md, t_fcdata gmx_unused__attribute__ ((unused)) *fcd, |
2553 | int gmx_unused__attribute__ ((unused)) *global_atom_index) |
2554 | { |
2555 | return low_angres(nbonds, forceatoms, forceparams, x, f, fshift, pbc, g, |
2556 | lambda, dvdlambda, TRUE1); |
2557 | } |
2558 | |
2559 | real dihres(int nbonds, |
2560 | const t_iatom forceatoms[], const t_iparams forceparams[], |
2561 | const rvec x[], rvec f[], rvec fshift[], |
2562 | const t_pbc *pbc, const t_graph *g, |
2563 | real lambda, real *dvdlambda, |
2564 | const t_mdatoms gmx_unused__attribute__ ((unused)) *md, t_fcdata gmx_unused__attribute__ ((unused)) *fcd, |
2565 | int gmx_unused__attribute__ ((unused)) *global_atom_index) |
2566 | { |
2567 | real vtot = 0; |
2568 | int ai, aj, ak, al, i, k, type, t1, t2, t3; |
2569 | real phi0A, phi0B, dphiA, dphiB, kfacA, kfacB, phi0, dphi, kfac; |
2570 | real phi, ddphi, ddp, ddp2, dp, sign, d2r, fc, L1; |
2571 | rvec r_ij, r_kj, r_kl, m, n; |
2572 | |
2573 | L1 = 1.0-lambda; |
2574 | |
2575 | d2r = DEG2RAD(3.14159265358979323846/180.0); |
2576 | k = 0; |
2577 | |
2578 | for (i = 0; (i < nbonds); ) |
2579 | { |
2580 | type = forceatoms[i++]; |
2581 | ai = forceatoms[i++]; |
2582 | aj = forceatoms[i++]; |
2583 | ak = forceatoms[i++]; |
2584 | al = forceatoms[i++]; |
2585 | |
2586 | phi0A = forceparams[type].dihres.phiA*d2r; |
2587 | dphiA = forceparams[type].dihres.dphiA*d2r; |
2588 | kfacA = forceparams[type].dihres.kfacA; |
2589 | |
2590 | phi0B = forceparams[type].dihres.phiB*d2r; |
2591 | dphiB = forceparams[type].dihres.dphiB*d2r; |
2592 | kfacB = forceparams[type].dihres.kfacB; |
2593 | |
2594 | phi0 = L1*phi0A + lambda*phi0B; |
2595 | dphi = L1*dphiA + lambda*dphiB; |
2596 | kfac = L1*kfacA + lambda*kfacB; |
2597 | |
2598 | phi = dih_angle(x[ai], x[aj], x[ak], x[al], pbc, r_ij, r_kj, r_kl, m, n, |
2599 | &sign, &t1, &t2, &t3); |
2600 | /* 84 flops */ |
2601 | |
2602 | if (debug) |
2603 | { |
2604 | fprintf(debug, "dihres[%d]: %d %d %d %d : phi=%f, dphi=%f, kfac=%f\n", |
2605 | k++, ai, aj, ak, al, phi0, dphi, kfac); |
2606 | } |
2607 | /* phi can jump if phi0 is close to Pi/-Pi, which will cause huge |
2608 | * force changes if we just apply a normal harmonic. |
2609 | * Instead, we first calculate phi-phi0 and take it modulo (-Pi,Pi). |
2610 | * This means we will never have the periodicity problem, unless |
2611 | * the dihedral is Pi away from phiO, which is very unlikely due to |
2612 | * the potential. |
2613 | */ |
2614 | dp = phi-phi0; |
2615 | make_dp_periodic(&dp); |
2616 | |
2617 | if (dp > dphi) |
2618 | { |
2619 | ddp = dp-dphi; |
2620 | } |
2621 | else if (dp < -dphi) |
2622 | { |
2623 | ddp = dp+dphi; |
2624 | } |
2625 | else |
2626 | { |
2627 | ddp = 0; |
2628 | } |
2629 | |
2630 | if (ddp != 0.0) |
2631 | { |
2632 | ddp2 = ddp*ddp; |
2633 | vtot += 0.5*kfac*ddp2; |
2634 | ddphi = kfac*ddp; |
2635 | |
2636 | *dvdlambda += 0.5*(kfacB - kfacA)*ddp2; |
2637 | /* lambda dependence from changing restraint distances */ |
2638 | if (ddp > 0) |
2639 | { |
2640 | *dvdlambda -= kfac*ddp*((dphiB - dphiA)+(phi0B - phi0A)); |
2641 | } |
2642 | else if (ddp < 0) |
2643 | { |
2644 | *dvdlambda += kfac*ddp*((dphiB - dphiA)-(phi0B - phi0A)); |
2645 | } |
2646 | do_dih_fup(ai, aj, ak, al, ddphi, r_ij, r_kj, r_kl, m, n, |
2647 | f, fshift, pbc, g, x, t1, t2, t3); /* 112 */ |
2648 | } |
2649 | } |
2650 | return vtot; |
2651 | } |
2652 | |
2653 | |
2654 | real unimplemented(int gmx_unused__attribute__ ((unused)) nbonds, |
2655 | const t_iatom gmx_unused__attribute__ ((unused)) forceatoms[], const t_iparams gmx_unused__attribute__ ((unused)) forceparams[], |
2656 | const rvec gmx_unused__attribute__ ((unused)) x[], rvec gmx_unused__attribute__ ((unused)) f[], rvec gmx_unused__attribute__ ((unused)) fshift[], |
2657 | const t_pbc gmx_unused__attribute__ ((unused)) *pbc, const t_graph gmx_unused__attribute__ ((unused)) *g, |
2658 | real gmx_unused__attribute__ ((unused)) lambda, real gmx_unused__attribute__ ((unused)) *dvdlambda, |
2659 | const t_mdatoms gmx_unused__attribute__ ((unused)) *md, t_fcdata gmx_unused__attribute__ ((unused)) *fcd, |
2660 | int gmx_unused__attribute__ ((unused)) *global_atom_index) |
2661 | { |
2662 | gmx_impl("*** you are using a not implemented function")_gmx_error("impl", "*** you are using a not implemented function" , "/home/alexxy/Develop/gromacs/src/gromacs/gmxlib/bondfree.c" , 2662); |
2663 | |
2664 | return 0.0; /* To make the compiler happy */ |
2665 | } |
2666 | |
2667 | real restrangles(int nbonds, |
2668 | const t_iatom forceatoms[], const t_iparams forceparams[], |
2669 | const rvec x[], rvec f[], rvec fshift[], |
2670 | const t_pbc *pbc, const t_graph *g, |
2671 | real gmx_unused__attribute__ ((unused)) lambda, real gmx_unused__attribute__ ((unused)) *dvdlambda, |
2672 | const t_mdatoms gmx_unused__attribute__ ((unused)) *md, t_fcdata gmx_unused__attribute__ ((unused)) *fcd, |
2673 | int gmx_unused__attribute__ ((unused)) *global_atom_index) |
2674 | { |
2675 | int i, d, ai, aj, ak, type, m; |
2676 | int t1, t2; |
2677 | rvec r_ij, r_kj; |
2678 | real v, vtot; |
2679 | ivec jt, dt_ij, dt_kj; |
2680 | rvec f_i, f_j, f_k; |
2681 | real prefactor, ratio_ante, ratio_post; |
2682 | rvec delta_ante, delta_post, vec_temp; |
2683 | |
2684 | vtot = 0.0; |
2685 | for (i = 0; (i < nbonds); ) |
2686 | { |
2687 | type = forceatoms[i++]; |
2688 | ai = forceatoms[i++]; |
2689 | aj = forceatoms[i++]; |
2690 | ak = forceatoms[i++]; |
2691 | |
2692 | t1 = pbc_rvec_sub(pbc, x[ai], x[aj], vec_temp); |
2693 | pbc_rvec_sub(pbc, x[aj], x[ai], delta_ante); |
2694 | t2 = pbc_rvec_sub(pbc, x[ak], x[aj], delta_post); |
2695 | |
2696 | |
2697 | /* This function computes factors needed for restricted angle potential. |
2698 | * The restricted angle potential is used in coarse-grained simulations to avoid singularities |
2699 | * when three particles align and the dihedral angle and dihedral potential |
2700 | * cannot be calculated. This potential is calculated using the formula: |
2701 | real restrangles(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 *dvdlambda, |
2706 | const t_mdatoms gmx_unused *md,t_fcdata gmx_unused *fcd, |
2707 | int gmx_unused *global_atom_index) |
2708 | { |
2709 | int i, d, ai, aj, ak, type, m; |
2710 | int t1, t2; |
2711 | rvec r_ij,r_kj; |
2712 | real v, vtot; |
2713 | ivec jt,dt_ij,dt_kj; |
2714 | rvec f_i, f_j, f_k; |
2715 | real prefactor, ratio_ante, ratio_post; |
2716 | rvec delta_ante, delta_post, vec_temp; |
2717 | |
2718 | vtot = 0.0; |
2719 | for(i=0; (i<nbonds); ) |
2720 | { |
2721 | type = forceatoms[i++]; |
2722 | ai = forceatoms[i++]; |
2723 | aj = forceatoms[i++]; |
2724 | ak = forceatoms[i++]; |
2725 | |
2726 | * \f[V_{\rm ReB}(\theta_i) = \frac{1}{2} k_{\theta} \frac{(\cos\theta_i - \cos\theta_0)^2} |
2727 | * {\sin^2\theta_i}\f] ({eq:ReB} and ref \cite{MonicaGoga2013} from the manual). |
2728 | * For more explanations see comments file "restcbt.h". */ |
2729 | |
2730 | compute_factors_restangles(type, forceparams, delta_ante, delta_post, |
2731 | &prefactor, &ratio_ante, &ratio_post, &v); |
2732 | |
2733 | /* Forces are computed per component */ |
2734 | for (d = 0; d < DIM3; d++) |
2735 | { |
2736 | f_i[d] = prefactor * (ratio_ante * delta_ante[d] - delta_post[d]); |
2737 | f_j[d] = prefactor * ((ratio_post + 1.0) * delta_post[d] - (ratio_ante + 1.0) * delta_ante[d]); |
2738 | f_k[d] = prefactor * (delta_ante[d] - ratio_post * delta_post[d]); |
2739 | } |
2740 | |
2741 | /* Computation of potential energy */ |
2742 | |
2743 | vtot += v; |
2744 | |
2745 | /* Update forces */ |
2746 | |
2747 | for (m = 0; (m < DIM3); m++) |
2748 | { |
2749 | f[ai][m] += f_i[m]; |
2750 | f[aj][m] += f_j[m]; |
2751 | f[ak][m] += f_k[m]; |
2752 | } |
2753 | |
2754 | if (g) |
2755 | { |
2756 | copy_ivec(SHIFT_IVEC(g, aj)((g)->ishift[aj]), jt); |
2757 | ivec_sub(SHIFT_IVEC(g, ai)((g)->ishift[ai]), jt, dt_ij); |
2758 | ivec_sub(SHIFT_IVEC(g, ak)((g)->ishift[ak]), jt, dt_kj); |
2759 | t1 = IVEC2IS(dt_ij)(((2*2 +1)*((2*1 +1)*(((dt_ij)[2])+1)+((dt_ij)[1])+1)+((dt_ij )[0])+2)); |
2760 | t2 = IVEC2IS(dt_kj)(((2*2 +1)*((2*1 +1)*(((dt_kj)[2])+1)+((dt_kj)[1])+1)+((dt_kj )[0])+2)); |
2761 | } |
2762 | |
2763 | rvec_inc(fshift[t1], f_i); |
2764 | rvec_inc(fshift[CENTRAL(((2*1 +1)*(2*1 +1)*(2*2 +1))/2)], f_j); |
2765 | rvec_inc(fshift[t2], f_k); |
2766 | } |
2767 | return vtot; |
2768 | } |
2769 | |
2770 | |
2771 | real restrdihs(int nbonds, |
2772 | const t_iatom forceatoms[], const t_iparams forceparams[], |
2773 | const rvec x[], rvec f[], rvec fshift[], |
2774 | const t_pbc *pbc, const t_graph *g, |
2775 | real gmx_unused__attribute__ ((unused)) lambda, real gmx_unused__attribute__ ((unused)) *dvlambda, |
2776 | const t_mdatoms gmx_unused__attribute__ ((unused)) *md, t_fcdata gmx_unused__attribute__ ((unused)) *fcd, |
2777 | int gmx_unused__attribute__ ((unused)) *global_atom_index) |
2778 | { |
2779 | int i, d, type, ai, aj, ak, al; |
2780 | rvec f_i, f_j, f_k, f_l; |
2781 | rvec dx_jl; |
2782 | ivec jt, dt_ij, dt_kj, dt_lj; |
2783 | int t1, t2, t3; |
2784 | real v, vtot; |
2785 | rvec delta_ante, delta_crnt, delta_post, vec_temp; |
2786 | real factor_phi_ai_ante, factor_phi_ai_crnt, factor_phi_ai_post; |
2787 | real factor_phi_aj_ante, factor_phi_aj_crnt, factor_phi_aj_post; |
2788 | real factor_phi_ak_ante, factor_phi_ak_crnt, factor_phi_ak_post; |
2789 | real factor_phi_al_ante, factor_phi_al_crnt, factor_phi_al_post; |
2790 | real prefactor_phi; |
2791 | |
2792 | |
2793 | vtot = 0.0; |
2794 | for (i = 0; (i < nbonds); ) |
2795 | { |
2796 | type = forceatoms[i++]; |
2797 | ai = forceatoms[i++]; |
2798 | aj = forceatoms[i++]; |
2799 | ak = forceatoms[i++]; |
2800 | al = forceatoms[i++]; |
2801 | |
2802 | t1 = pbc_rvec_sub(pbc, x[ai], x[aj], vec_temp); |
2803 | pbc_rvec_sub(pbc, x[aj], x[ai], delta_ante); |
2804 | t2 = pbc_rvec_sub(pbc, x[ak], x[aj], delta_crnt); |
2805 | t3 = pbc_rvec_sub(pbc, x[ak], x[al], vec_temp); |
2806 | pbc_rvec_sub(pbc, x[al], x[ak], delta_post); |
2807 | |
2808 | /* This function computes factors needed for restricted angle potential. |
2809 | * The restricted angle potential is used in coarse-grained simulations to avoid singularities |
2810 | * when three particles align and the dihedral angle and dihedral potential cannot be calculated. |
2811 | * This potential is calculated using the formula: |
2812 | * \f[V_{\rm ReB}(\theta_i) = \frac{1}{2} k_{\theta} |
2813 | * \frac{(\cos\theta_i - \cos\theta_0)^2}{\sin^2\theta_i}\f] |
2814 | * ({eq:ReB} and ref \cite{MonicaGoga2013} from the manual). |
2815 | * For more explanations see comments file "restcbt.h" */ |
2816 | |
2817 | compute_factors_restrdihs(type, forceparams, |
2818 | delta_ante, delta_crnt, delta_post, |
2819 | &factor_phi_ai_ante, &factor_phi_ai_crnt, &factor_phi_ai_post, |
2820 | &factor_phi_aj_ante, &factor_phi_aj_crnt, &factor_phi_aj_post, |
2821 | &factor_phi_ak_ante, &factor_phi_ak_crnt, &factor_phi_ak_post, |
2822 | &factor_phi_al_ante, &factor_phi_al_crnt, &factor_phi_al_post, |
2823 | &prefactor_phi, &v); |
2824 | |
2825 | |
2826 | /* Computation of forces per component */ |
2827 | for (d = 0; d < DIM3; d++) |
2828 | { |
2829 | 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]); |
2830 | 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]); |
2831 | 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]); |
2832 | 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]); |
2833 | } |
2834 | /* Computation of the energy */ |
2835 | |
2836 | vtot += v; |
2837 | |
2838 | |
2839 | |
2840 | /* Updating the forces */ |
2841 | |
2842 | rvec_inc(f[ai], f_i); |
2843 | rvec_inc(f[aj], f_j); |
2844 | rvec_inc(f[ak], f_k); |
2845 | rvec_inc(f[al], f_l); |
2846 | |
2847 | |
2848 | /* Updating the fshift forces for the pressure coupling */ |
2849 | if (g) |
2850 | { |
2851 | copy_ivec(SHIFT_IVEC(g, aj)((g)->ishift[aj]), jt); |
2852 | ivec_sub(SHIFT_IVEC(g, ai)((g)->ishift[ai]), jt, dt_ij); |
2853 | ivec_sub(SHIFT_IVEC(g, ak)((g)->ishift[ak]), jt, dt_kj); |
2854 | ivec_sub(SHIFT_IVEC(g, al)((g)->ishift[al]), jt, dt_lj); |
2855 | t1 = IVEC2IS(dt_ij)(((2*2 +1)*((2*1 +1)*(((dt_ij)[2])+1)+((dt_ij)[1])+1)+((dt_ij )[0])+2)); |
2856 | t2 = IVEC2IS(dt_kj)(((2*2 +1)*((2*1 +1)*(((dt_kj)[2])+1)+((dt_kj)[1])+1)+((dt_kj )[0])+2)); |
2857 | t3 = IVEC2IS(dt_lj)(((2*2 +1)*((2*1 +1)*(((dt_lj)[2])+1)+((dt_lj)[1])+1)+((dt_lj )[0])+2)); |
2858 | } |
2859 | else if (pbc) |
2860 | { |
2861 | t3 = pbc_rvec_sub(pbc, x[al], x[aj], dx_jl); |
2862 | } |
2863 | else |
2864 | { |
2865 | t3 = CENTRAL(((2*1 +1)*(2*1 +1)*(2*2 +1))/2); |
2866 | } |
2867 | |
2868 | rvec_inc(fshift[t1], f_i); |
2869 | rvec_inc(fshift[CENTRAL(((2*1 +1)*(2*1 +1)*(2*2 +1))/2)], f_j); |
2870 | rvec_inc(fshift[t2], f_k); |
2871 | rvec_inc(fshift[t3], f_l); |
2872 | |
2873 | } |
2874 | |
2875 | return vtot; |
2876 | } |
2877 | |
2878 | |
2879 | real cbtdihs(int nbonds, |
2880 | const t_iatom forceatoms[], const t_iparams forceparams[], |
2881 | const rvec x[], rvec f[], rvec fshift[], |
2882 | const t_pbc *pbc, const t_graph *g, |
2883 | real gmx_unused__attribute__ ((unused)) lambda, real gmx_unused__attribute__ ((unused)) *dvdlambda, |
2884 | const t_mdatoms gmx_unused__attribute__ ((unused)) *md, t_fcdata gmx_unused__attribute__ ((unused)) *fcd, |
2885 | int gmx_unused__attribute__ ((unused)) *global_atom_index) |
2886 | { |
2887 | int type, ai, aj, ak, al, i, d; |
2888 | int t1, t2, t3; |
2889 | real v, vtot; |
2890 | rvec vec_temp; |
2891 | rvec f_i, f_j, f_k, f_l; |
2892 | ivec jt, dt_ij, dt_kj, dt_lj; |
2893 | rvec dx_jl; |
2894 | rvec delta_ante, delta_crnt, delta_post; |
2895 | rvec f_phi_ai, f_phi_aj, f_phi_ak, f_phi_al; |
2896 | rvec f_theta_ante_ai, f_theta_ante_aj, f_theta_ante_ak; |
2897 | rvec f_theta_post_aj, f_theta_post_ak, f_theta_post_al; |
2898 | |
2899 | |
2900 | |
2901 | |
2902 | vtot = 0.0; |
2903 | for (i = 0; (i < nbonds); ) |
2904 | { |
2905 | type = forceatoms[i++]; |
2906 | ai = forceatoms[i++]; |
2907 | aj = forceatoms[i++]; |
2908 | ak = forceatoms[i++]; |
2909 | al = forceatoms[i++]; |
2910 | |
2911 | |
2912 | t1 = pbc_rvec_sub(pbc, x[ai], x[aj], vec_temp); |
2913 | pbc_rvec_sub(pbc, x[aj], x[ai], delta_ante); |
2914 | t2 = pbc_rvec_sub(pbc, x[ak], x[aj], vec_temp); |
2915 | pbc_rvec_sub(pbc, x[ak], x[aj], delta_crnt); |
2916 | t3 = pbc_rvec_sub(pbc, x[ak], x[al], vec_temp); |
2917 | pbc_rvec_sub(pbc, x[al], x[ak], delta_post); |
2918 | |
2919 | /* \brief Compute factors for CBT potential |
2920 | * The combined bending-torsion potential goes to zero in a very smooth manner, eliminating the numerical |
2921 | * instabilities, when three coarse-grained particles align and the dihedral angle and standard |
2922 | * dihedral potentials cannot be calculated. The CBT potential is calculated using the formula: |
2923 | * \f[V_{\rm CBT}(\theta_{i-1}, \theta_i, \phi_i) = k_{\phi} \sin^3\theta_{i-1} \sin^3\theta_{i} |
2924 | * \sum_{n=0}^4 { a_n \cos^n\phi_i}\f] ({eq:CBT} and ref \cite{MonicaGoga2013} from the manual). |
2925 | * It contains in its expression not only the dihedral angle \f$\phi\f$ |
2926 | * but also \f[\theta_{i-1}\f] (theta_ante bellow) and \f[\theta_{i}\f] (theta_post bellow) |
2927 | * --- the adjacent bending angles. |
2928 | * For more explanations see comments file "restcbt.h". */ |
2929 | |
2930 | compute_factors_cbtdihs(type, forceparams, delta_ante, delta_crnt, delta_post, |
2931 | f_phi_ai, f_phi_aj, f_phi_ak, f_phi_al, |
2932 | f_theta_ante_ai, f_theta_ante_aj, f_theta_ante_ak, |
2933 | f_theta_post_aj, f_theta_post_ak, f_theta_post_al, |
2934 | &v); |
2935 | |
2936 | |
2937 | /* Acumulate the resuts per beads */ |
2938 | for (d = 0; d < DIM3; d++) |
2939 | { |
2940 | f_i[d] = f_phi_ai[d] + f_theta_ante_ai[d]; |
2941 | f_j[d] = f_phi_aj[d] + f_theta_ante_aj[d] + f_theta_post_aj[d]; |
2942 | f_k[d] = f_phi_ak[d] + f_theta_ante_ak[d] + f_theta_post_ak[d]; |
2943 | f_l[d] = f_phi_al[d] + f_theta_post_al[d]; |
2944 | } |
2945 | |
2946 | /* Compute the potential energy */ |
2947 | |
2948 | vtot += v; |
2949 | |
2950 | |
2951 | /* Updating the forces */ |
2952 | rvec_inc(f[ai], f_i); |
2953 | rvec_inc(f[aj], f_j); |
2954 | rvec_inc(f[ak], f_k); |
2955 | rvec_inc(f[al], f_l); |
2956 | |
2957 | |
2958 | /* Updating the fshift forces for the pressure coupling */ |
2959 | if (g) |
2960 | { |
2961 | copy_ivec(SHIFT_IVEC(g, aj)((g)->ishift[aj]), jt); |
2962 | ivec_sub(SHIFT_IVEC(g, ai)((g)->ishift[ai]), jt, dt_ij); |
2963 | ivec_sub(SHIFT_IVEC(g, ak)((g)->ishift[ak]), jt, dt_kj); |
2964 | ivec_sub(SHIFT_IVEC(g, al)((g)->ishift[al]), jt, dt_lj); |
2965 | t1 = IVEC2IS(dt_ij)(((2*2 +1)*((2*1 +1)*(((dt_ij)[2])+1)+((dt_ij)[1])+1)+((dt_ij )[0])+2)); |
2966 | t2 = IVEC2IS(dt_kj)(((2*2 +1)*((2*1 +1)*(((dt_kj)[2])+1)+((dt_kj)[1])+1)+((dt_kj )[0])+2)); |
2967 | t3 = IVEC2IS(dt_lj)(((2*2 +1)*((2*1 +1)*(((dt_lj)[2])+1)+((dt_lj)[1])+1)+((dt_lj )[0])+2)); |
2968 | } |
2969 | else if (pbc) |
2970 | { |
2971 | t3 = pbc_rvec_sub(pbc, x[al], x[aj], dx_jl); |
2972 | } |
2973 | else |
2974 | { |
2975 | t3 = CENTRAL(((2*1 +1)*(2*1 +1)*(2*2 +1))/2); |
2976 | } |
2977 | |
2978 | rvec_inc(fshift[t1], f_i); |
2979 | rvec_inc(fshift[CENTRAL(((2*1 +1)*(2*1 +1)*(2*2 +1))/2)], f_j); |
2980 | rvec_inc(fshift[t2], f_k); |
2981 | rvec_inc(fshift[t3], f_l); |
2982 | } |
2983 | |
2984 | return vtot; |
2985 | } |
2986 | |
2987 | real rbdihs(int nbonds, |
2988 | const t_iatom forceatoms[], const t_iparams forceparams[], |
2989 | const rvec x[], rvec f[], rvec fshift[], |
2990 | const t_pbc *pbc, const t_graph *g, |
2991 | real lambda, real *dvdlambda, |
2992 | const t_mdatoms gmx_unused__attribute__ ((unused)) *md, t_fcdata gmx_unused__attribute__ ((unused)) *fcd, |
2993 | int gmx_unused__attribute__ ((unused)) *global_atom_index) |
2994 | { |
2995 | const real c0 = 0.0, c1 = 1.0, c2 = 2.0, c3 = 3.0, c4 = 4.0, c5 = 5.0; |
2996 | int type, ai, aj, ak, al, i, j; |
2997 | int t1, t2, t3; |
2998 | rvec r_ij, r_kj, r_kl, m, n; |
2999 | real parmA[NR_RBDIHS6]; |
3000 | real parmB[NR_RBDIHS6]; |
3001 | real parm[NR_RBDIHS6]; |
3002 | real cos_phi, phi, rbp, rbpBA; |
3003 | real v, sign, ddphi, sin_phi; |
3004 | real cosfac, vtot; |
3005 | real L1 = 1.0-lambda; |
3006 | real dvdl_term = 0; |
3007 | |
3008 | vtot = 0.0; |
3009 | for (i = 0; (i < nbonds); ) |
3010 | { |
3011 | type = forceatoms[i++]; |
3012 | ai = forceatoms[i++]; |
3013 | aj = forceatoms[i++]; |
3014 | ak = forceatoms[i++]; |
3015 | al = forceatoms[i++]; |
3016 | |
3017 | phi = dih_angle(x[ai], x[aj], x[ak], x[al], pbc, r_ij, r_kj, r_kl, m, n, |
3018 | &sign, &t1, &t2, &t3); /* 84 */ |
3019 | |
3020 | /* Change to polymer convention */ |
3021 | if (phi < c0) |
3022 | { |
3023 | phi += M_PI3.14159265358979323846; |
3024 | } |
3025 | else |
3026 | { |
3027 | phi -= M_PI3.14159265358979323846; /* 1 */ |
3028 | |
3029 | } |
3030 | cos_phi = cos(phi); |
3031 | /* Beware of accuracy loss, cannot use 1-sqrt(cos^2) ! */ |
3032 | sin_phi = sin(phi); |
3033 | |
3034 | for (j = 0; (j < NR_RBDIHS6); j++) |
3035 | { |
3036 | parmA[j] = forceparams[type].rbdihs.rbcA[j]; |
3037 | parmB[j] = forceparams[type].rbdihs.rbcB[j]; |
3038 | parm[j] = L1*parmA[j]+lambda*parmB[j]; |
3039 | } |
3040 | /* Calculate cosine powers */ |
3041 | /* Calculate the energy */ |
3042 | /* Calculate the derivative */ |
3043 | |
3044 | v = parm[0]; |
3045 | dvdl_term += (parmB[0]-parmA[0]); |
3046 | ddphi = c0; |
3047 | cosfac = c1; |
3048 | |
3049 | rbp = parm[1]; |
3050 | rbpBA = parmB[1]-parmA[1]; |
3051 | ddphi += rbp*cosfac; |
3052 | cosfac *= cos_phi; |
3053 | v += cosfac*rbp; |
3054 | dvdl_term += cosfac*rbpBA; |
3055 | rbp = parm[2]; |
3056 | rbpBA = parmB[2]-parmA[2]; |
3057 | ddphi += c2*rbp*cosfac; |
3058 | cosfac *= cos_phi; |
3059 | v += cosfac*rbp; |
3060 | dvdl_term += cosfac*rbpBA; |
3061 | rbp = parm[3]; |
3062 | rbpBA = parmB[3]-parmA[3]; |
3063 | ddphi += c3*rbp*cosfac; |
3064 | cosfac *= cos_phi; |
3065 | v += cosfac*rbp; |
3066 | dvdl_term += cosfac*rbpBA; |
3067 | rbp = parm[4]; |
3068 | rbpBA = parmB[4]-parmA[4]; |
3069 | ddphi += c4*rbp*cosfac; |
3070 | cosfac *= cos_phi; |
3071 | v += cosfac*rbp; |
3072 | dvdl_term += cosfac*rbpBA; |
3073 | rbp = parm[5]; |
3074 | rbpBA = parmB[5]-parmA[5]; |
3075 | ddphi += c5*rbp*cosfac; |
3076 | cosfac *= cos_phi; |
3077 | v += cosfac*rbp; |
3078 | dvdl_term += cosfac*rbpBA; |
3079 | |
3080 | ddphi = -ddphi*sin_phi; /* 11 */ |
3081 | |
3082 | do_dih_fup(ai, aj, ak, al, ddphi, r_ij, r_kj, r_kl, m, n, |
3083 | f, fshift, pbc, g, x, t1, t2, t3); /* 112 */ |
3084 | vtot += v; |
3085 | } |
3086 | *dvdlambda += dvdl_term; |
3087 | |
3088 | return vtot; |
3089 | } |
3090 | |
3091 | int cmap_setup_grid_index(int ip, int grid_spacing, int *ipm1, int *ipp1, int *ipp2) |
3092 | { |
3093 | int im1, ip1, ip2; |
3094 | |
3095 | if (ip < 0) |
3096 | { |
3097 | ip = ip + grid_spacing - 1; |
3098 | } |
3099 | else if (ip > grid_spacing) |
3100 | { |
3101 | ip = ip - grid_spacing - 1; |
3102 | } |
3103 | |
3104 | im1 = ip - 1; |
3105 | ip1 = ip + 1; |
3106 | ip2 = ip + 2; |
3107 | |
3108 | if (ip == 0) |
3109 | { |
3110 | im1 = grid_spacing - 1; |
3111 | } |
3112 | else if (ip == grid_spacing-2) |
3113 | { |
3114 | ip2 = 0; |
3115 | } |
3116 | else if (ip == grid_spacing-1) |
3117 | { |
3118 | ip1 = 0; |
3119 | ip2 = 1; |
3120 | } |
3121 | |
3122 | *ipm1 = im1; |
3123 | *ipp1 = ip1; |
3124 | *ipp2 = ip2; |
3125 | |
3126 | return ip; |
3127 | |
3128 | } |
3129 | |
3130 | real cmap_dihs(int nbonds, |
3131 | const t_iatom forceatoms[], const t_iparams forceparams[], |
3132 | const gmx_cmap_t *cmap_grid, |
3133 | const rvec x[], rvec f[], rvec fshift[], |
3134 | const t_pbc *pbc, const t_graph *g, |
3135 | real gmx_unused__attribute__ ((unused)) lambda, real gmx_unused__attribute__ ((unused)) *dvdlambda, |
3136 | const t_mdatoms gmx_unused__attribute__ ((unused)) *md, t_fcdata gmx_unused__attribute__ ((unused)) *fcd, |
3137 | int gmx_unused__attribute__ ((unused)) *global_atom_index) |
3138 | { |
3139 | int i, j, k, n, idx; |
3140 | int ai, aj, ak, al, am; |
3141 | int a1i, a1j, a1k, a1l, a2i, a2j, a2k, a2l; |
3142 | int type, cmapA; |
3143 | int t11, t21, t31, t12, t22, t32; |
3144 | int iphi1, ip1m1, ip1p1, ip1p2; |
3145 | int iphi2, ip2m1, ip2p1, ip2p2; |
3146 | int l1, l2, l3, l4; |
3147 | int pos1, pos2, pos3, pos4, tmp; |
3148 | |
3149 | real ty[4], ty1[4], ty2[4], ty12[4], tc[16], tx[16]; |
3150 | real phi1, psi1, cos_phi1, sin_phi1, sign1, xphi1; |
3151 | real phi2, psi2, cos_phi2, sin_phi2, sign2, xphi2; |
3152 | real dx, xx, tt, tu, e, df1, df2, ddf1, ddf2, ddf12, vtot; |
3153 | real ra21, rb21, rg21, rg1, rgr1, ra2r1, rb2r1, rabr1; |
3154 | real ra22, rb22, rg22, rg2, rgr2, ra2r2, rb2r2, rabr2; |
3155 | real fg1, hg1, fga1, hgb1, gaa1, gbb1; |
3156 | real fg2, hg2, fga2, hgb2, gaa2, gbb2; |
3157 | real fac; |
3158 | |
3159 | rvec r1_ij, r1_kj, r1_kl, m1, n1; |
3160 | rvec r2_ij, r2_kj, r2_kl, m2, n2; |
3161 | rvec f1_i, f1_j, f1_k, f1_l; |
3162 | rvec f2_i, f2_j, f2_k, f2_l; |
3163 | rvec a1, b1, a2, b2; |
3164 | rvec f1, g1, h1, f2, g2, h2; |
3165 | rvec dtf1, dtg1, dth1, dtf2, dtg2, dth2; |
3166 | ivec jt1, dt1_ij, dt1_kj, dt1_lj; |
3167 | ivec jt2, dt2_ij, dt2_kj, dt2_lj; |
3168 | |
3169 | const real *cmapd; |
3170 | |
3171 | int loop_index[4][4] = { |
3172 | {0, 4, 8, 12}, |
3173 | {1, 5, 9, 13}, |
3174 | {2, 6, 10, 14}, |
3175 | {3, 7, 11, 15} |
3176 | }; |
3177 | |
3178 | /* Total CMAP energy */ |
3179 | vtot = 0; |
3180 | |
3181 | for (n = 0; n < nbonds; ) |
3182 | { |
3183 | /* Five atoms are involved in the two torsions */ |
3184 | type = forceatoms[n++]; |
3185 | ai = forceatoms[n++]; |
3186 | aj = forceatoms[n++]; |
3187 | ak = forceatoms[n++]; |
3188 | al = forceatoms[n++]; |
3189 | am = forceatoms[n++]; |
3190 | |
3191 | /* Which CMAP type is this */ |
3192 | cmapA = forceparams[type].cmap.cmapA; |
3193 | cmapd = cmap_grid->cmapdata[cmapA].cmap; |
3194 | |
3195 | /* First torsion */ |
3196 | a1i = ai; |
3197 | a1j = aj; |
3198 | a1k = ak; |
3199 | a1l = al; |
3200 | |
3201 | phi1 = dih_angle(x[a1i], x[a1j], x[a1k], x[a1l], pbc, r1_ij, r1_kj, r1_kl, m1, n1, |
3202 | &sign1, &t11, &t21, &t31); /* 84 */ |
3203 | |
3204 | cos_phi1 = cos(phi1); |
3205 | |
3206 | a1[0] = r1_ij[1]*r1_kj[2]-r1_ij[2]*r1_kj[1]; |
3207 | a1[1] = r1_ij[2]*r1_kj[0]-r1_ij[0]*r1_kj[2]; |
3208 | a1[2] = r1_ij[0]*r1_kj[1]-r1_ij[1]*r1_kj[0]; /* 9 */ |
3209 | |
3210 | b1[0] = r1_kl[1]*r1_kj[2]-r1_kl[2]*r1_kj[1]; |
3211 | b1[1] = r1_kl[2]*r1_kj[0]-r1_kl[0]*r1_kj[2]; |
3212 | b1[2] = r1_kl[0]*r1_kj[1]-r1_kl[1]*r1_kj[0]; /* 9 */ |
3213 | |
3214 | tmp = pbc_rvec_sub(pbc, x[a1l], x[a1k], h1); |
3215 | |
3216 | ra21 = iprod(a1, a1); /* 5 */ |
3217 | rb21 = iprod(b1, b1); /* 5 */ |
3218 | rg21 = iprod(r1_kj, r1_kj); /* 5 */ |
3219 | rg1 = sqrt(rg21); |
3220 | |
3221 | rgr1 = 1.0/rg1; |
3222 | ra2r1 = 1.0/ra21; |
3223 | rb2r1 = 1.0/rb21; |
3224 | rabr1 = sqrt(ra2r1*rb2r1); |
3225 | |
3226 | sin_phi1 = rg1 * rabr1 * iprod(a1, h1) * (-1); |
3227 | |
3228 | if (cos_phi1 < -0.5 || cos_phi1 > 0.5) |
3229 | { |
3230 | phi1 = asin(sin_phi1); |
3231 | |
3232 | if (cos_phi1 < 0) |
3233 | { |
3234 | if (phi1 > 0) |
3235 | { |
3236 | phi1 = M_PI3.14159265358979323846 - phi1; |
3237 | } |
3238 | else |
3239 | { |
3240 | phi1 = -M_PI3.14159265358979323846 - phi1; |
3241 | } |
3242 | } |
3243 | } |
3244 | else |
3245 | { |
3246 | phi1 = acos(cos_phi1); |
3247 | |
3248 | if (sin_phi1 < 0) |
3249 | { |
3250 | phi1 = -phi1; |
3251 | } |
3252 | } |
3253 | |
3254 | xphi1 = phi1 + M_PI3.14159265358979323846; /* 1 */ |
3255 | |
3256 | /* Second torsion */ |
3257 | a2i = aj; |
3258 | a2j = ak; |
3259 | a2k = al; |
3260 | a2l = am; |
3261 | |
3262 | phi2 = dih_angle(x[a2i], x[a2j], x[a2k], x[a2l], pbc, r2_ij, r2_kj, r2_kl, m2, n2, |
3263 | &sign2, &t12, &t22, &t32); /* 84 */ |
3264 | |
3265 | cos_phi2 = cos(phi2); |
3266 | |
3267 | a2[0] = r2_ij[1]*r2_kj[2]-r2_ij[2]*r2_kj[1]; |
3268 | a2[1] = r2_ij[2]*r2_kj[0]-r2_ij[0]*r2_kj[2]; |
3269 | a2[2] = r2_ij[0]*r2_kj[1]-r2_ij[1]*r2_kj[0]; /* 9 */ |
3270 | |
3271 | b2[0] = r2_kl[1]*r2_kj[2]-r2_kl[2]*r2_kj[1]; |
3272 | b2[1] = r2_kl[2]*r2_kj[0]-r2_kl[0]*r2_kj[2]; |
3273 | b2[2] = r2_kl[0]*r2_kj[1]-r2_kl[1]*r2_kj[0]; /* 9 */ |
3274 | |
3275 | tmp = pbc_rvec_sub(pbc, x[a2l], x[a2k], h2); |
3276 | |
3277 | ra22 = iprod(a2, a2); /* 5 */ |
3278 | rb22 = iprod(b2, b2); /* 5 */ |
3279 | rg22 = iprod(r2_kj, r2_kj); /* 5 */ |
3280 | rg2 = sqrt(rg22); |
3281 | |
3282 | rgr2 = 1.0/rg2; |
3283 | ra2r2 = 1.0/ra22; |
3284 | rb2r2 = 1.0/rb22; |
3285 | rabr2 = sqrt(ra2r2*rb2r2); |
3286 | |
3287 | sin_phi2 = rg2 * rabr2 * iprod(a2, h2) * (-1); |
3288 | |
3289 | if (cos_phi2 < -0.5 || cos_phi2 > 0.5) |
3290 | { |
3291 | phi2 = asin(sin_phi2); |
3292 | |
3293 | if (cos_phi2 < 0) |
3294 | { |
3295 | if (phi2 > 0) |
3296 | { |
3297 | phi2 = M_PI3.14159265358979323846 - phi2; |
3298 | } |
3299 | else |
3300 | { |
3301 | phi2 = -M_PI3.14159265358979323846 - phi2; |
3302 | } |
3303 | } |
3304 | } |
3305 | else |
3306 | { |
3307 | phi2 = acos(cos_phi2); |
3308 | |
3309 | if (sin_phi2 < 0) |
3310 | { |
3311 | phi2 = -phi2; |
3312 | } |
3313 | } |
3314 | |
3315 | xphi2 = phi2 + M_PI3.14159265358979323846; /* 1 */ |
3316 | |
3317 | /* Range mangling */ |
3318 | if (xphi1 < 0) |
3319 | { |
3320 | xphi1 = xphi1 + 2*M_PI3.14159265358979323846; |
3321 | } |
3322 | else if (xphi1 >= 2*M_PI3.14159265358979323846) |
3323 | { |
3324 | xphi1 = xphi1 - 2*M_PI3.14159265358979323846; |
3325 | } |
3326 | |
3327 | if (xphi2 < 0) |
3328 | { |
3329 | xphi2 = xphi2 + 2*M_PI3.14159265358979323846; |
3330 | } |
3331 | else if (xphi2 >= 2*M_PI3.14159265358979323846) |
3332 | { |
3333 | xphi2 = xphi2 - 2*M_PI3.14159265358979323846; |
3334 | } |
3335 | |
3336 | /* Number of grid points */ |
3337 | dx = 2*M_PI3.14159265358979323846 / cmap_grid->grid_spacing; |
3338 | |
3339 | /* Where on the grid are we */ |
3340 | iphi1 = (int)(xphi1/dx); |
3341 | iphi2 = (int)(xphi2/dx); |
3342 | |
3343 | iphi1 = cmap_setup_grid_index(iphi1, cmap_grid->grid_spacing, &ip1m1, &ip1p1, &ip1p2); |
3344 | iphi2 = cmap_setup_grid_index(iphi2, cmap_grid->grid_spacing, &ip2m1, &ip2p1, &ip2p2); |
3345 | |
3346 | pos1 = iphi1*cmap_grid->grid_spacing+iphi2; |
3347 | pos2 = ip1p1*cmap_grid->grid_spacing+iphi2; |
3348 | pos3 = ip1p1*cmap_grid->grid_spacing+ip2p1; |
3349 | pos4 = iphi1*cmap_grid->grid_spacing+ip2p1; |
3350 | |
3351 | ty[0] = cmapd[pos1*4]; |
3352 | ty[1] = cmapd[pos2*4]; |
3353 | ty[2] = cmapd[pos3*4]; |
3354 | ty[3] = cmapd[pos4*4]; |
3355 | |
3356 | ty1[0] = cmapd[pos1*4+1]; |
3357 | ty1[1] = cmapd[pos2*4+1]; |
3358 | ty1[2] = cmapd[pos3*4+1]; |
3359 | ty1[3] = cmapd[pos4*4+1]; |
3360 | |
3361 | ty2[0] = cmapd[pos1*4+2]; |
3362 | ty2[1] = cmapd[pos2*4+2]; |
3363 | ty2[2] = cmapd[pos3*4+2]; |
3364 | ty2[3] = cmapd[pos4*4+2]; |
3365 | |
3366 | ty12[0] = cmapd[pos1*4+3]; |
3367 | ty12[1] = cmapd[pos2*4+3]; |
3368 | ty12[2] = cmapd[pos3*4+3]; |
3369 | ty12[3] = cmapd[pos4*4+3]; |
3370 | |
3371 | /* Switch to degrees */ |
3372 | dx = 360.0 / cmap_grid->grid_spacing; |
3373 | xphi1 = xphi1 * RAD2DEG(180.0/3.14159265358979323846); |
3374 | xphi2 = xphi2 * RAD2DEG(180.0/3.14159265358979323846); |
3375 | |
3376 | for (i = 0; i < 4; i++) /* 16 */ |
3377 | { |
3378 | tx[i] = ty[i]; |
3379 | tx[i+4] = ty1[i]*dx; |
3380 | tx[i+8] = ty2[i]*dx; |
3381 | tx[i+12] = ty12[i]*dx*dx; |
3382 | } |
3383 | |
3384 | idx = 0; |
3385 | for (i = 0; i < 4; i++) /* 1056 */ |
3386 | { |
3387 | for (j = 0; j < 4; j++) |
3388 | { |
3389 | xx = 0; |
3390 | for (k = 0; k < 16; k++) |
3391 | { |
3392 | xx = xx + cmap_coeff_matrix[k*16+idx]*tx[k]; |
3393 | } |
3394 | |
3395 | idx++; |
3396 | tc[i*4+j] = xx; |
3397 | } |
3398 | } |
3399 | |
3400 | tt = (xphi1-iphi1*dx)/dx; |
3401 | tu = (xphi2-iphi2*dx)/dx; |
3402 | |
3403 | e = 0; |
3404 | df1 = 0; |
3405 | df2 = 0; |
3406 | ddf1 = 0; |
3407 | ddf2 = 0; |
3408 | ddf12 = 0; |
3409 | |
3410 | for (i = 3; i >= 0; i--) |
3411 | { |
3412 | l1 = loop_index[i][3]; |
3413 | l2 = loop_index[i][2]; |
3414 | l3 = loop_index[i][1]; |
3415 | |
3416 | e = tt * e + ((tc[i*4+3]*tu+tc[i*4+2])*tu + tc[i*4+1])*tu+tc[i*4]; |
3417 | df1 = tu * df1 + (3.0*tc[l1]*tt+2.0*tc[l2])*tt+tc[l3]; |
3418 | df2 = tt * df2 + (3.0*tc[i*4+3]*tu+2.0*tc[i*4+2])*tu+tc[i*4+1]; |
3419 | ddf1 = tu * ddf1 + 2.0*3.0*tc[l1]*tt+2.0*tc[l2]; |
3420 | ddf2 = tt * ddf2 + 2.0*3.0*tc[4*i+3]*tu+2.0*tc[4*i+2]; |
3421 | } |
3422 | |
3423 | ddf12 = tc[5] + 2.0*tc[9]*tt + 3.0*tc[13]*tt*tt + 2.0*tu*(tc[6]+2.0*tc[10]*tt+3.0*tc[14]*tt*tt) + |
3424 | 3.0*tu*tu*(tc[7]+2.0*tc[11]*tt+3.0*tc[15]*tt*tt); |
3425 | |
3426 | fac = RAD2DEG(180.0/3.14159265358979323846)/dx; |
3427 | df1 = df1 * fac; |
3428 | df2 = df2 * fac; |
3429 | ddf1 = ddf1 * fac * fac; |
3430 | ddf2 = ddf2 * fac * fac; |
3431 | ddf12 = ddf12 * fac * fac; |
3432 | |
3433 | /* CMAP energy */ |
3434 | vtot += e; |
3435 | |
3436 | /* Do forces - first torsion */ |
3437 | fg1 = iprod(r1_ij, r1_kj); |
3438 | hg1 = iprod(r1_kl, r1_kj); |
3439 | fga1 = fg1*ra2r1*rgr1; |
3440 | hgb1 = hg1*rb2r1*rgr1; |
3441 | gaa1 = -ra2r1*rg1; |
3442 | gbb1 = rb2r1*rg1; |
3443 | |
3444 | for (i = 0; i < DIM3; i++) |
3445 | { |
3446 | dtf1[i] = gaa1 * a1[i]; |
3447 | dtg1[i] = fga1 * a1[i] - hgb1 * b1[i]; |
3448 | dth1[i] = gbb1 * b1[i]; |
3449 | |
3450 | f1[i] = df1 * dtf1[i]; |
3451 | g1[i] = df1 * dtg1[i]; |
3452 | h1[i] = df1 * dth1[i]; |
3453 | |
3454 | f1_i[i] = f1[i]; |
3455 | f1_j[i] = -f1[i] - g1[i]; |
3456 | f1_k[i] = h1[i] + g1[i]; |
3457 | f1_l[i] = -h1[i]; |
3458 | |
3459 | f[a1i][i] = f[a1i][i] + f1_i[i]; |
3460 | f[a1j][i] = f[a1j][i] + f1_j[i]; /* - f1[i] - g1[i] */ |
3461 | f[a1k][i] = f[a1k][i] + f1_k[i]; /* h1[i] + g1[i] */ |
3462 | f[a1l][i] = f[a1l][i] + f1_l[i]; /* h1[i] */ |
3463 | } |
3464 | |
3465 | /* Do forces - second torsion */ |
3466 | fg2 = iprod(r2_ij, r2_kj); |
3467 | hg2 = iprod(r2_kl, r2_kj); |
3468 | fga2 = fg2*ra2r2*rgr2; |
3469 | hgb2 = hg2*rb2r2*rgr2; |
3470 | gaa2 = -ra2r2*rg2; |
3471 | gbb2 = rb2r2*rg2; |
3472 | |
3473 | for (i = 0; i < DIM3; i++) |
3474 | { |
3475 | dtf2[i] = gaa2 * a2[i]; |
3476 | dtg2[i] = fga2 * a2[i] - hgb2 * b2[i]; |
3477 | dth2[i] = gbb2 * b2[i]; |
3478 | |
3479 | f2[i] = df2 * dtf2[i]; |
3480 | g2[i] = df2 * dtg2[i]; |
3481 | h2[i] = df2 * dth2[i]; |
3482 | |
3483 | f2_i[i] = f2[i]; |
3484 | f2_j[i] = -f2[i] - g2[i]; |
3485 | f2_k[i] = h2[i] + g2[i]; |
3486 | f2_l[i] = -h2[i]; |
3487 | |
3488 | f[a2i][i] = f[a2i][i] + f2_i[i]; /* f2[i] */ |
3489 | f[a2j][i] = f[a2j][i] + f2_j[i]; /* - f2[i] - g2[i] */ |
3490 | f[a2k][i] = f[a2k][i] + f2_k[i]; /* h2[i] + g2[i] */ |
3491 | f[a2l][i] = f[a2l][i] + f2_l[i]; /* - h2[i] */ |
3492 | } |
3493 | |
3494 | /* Shift forces */ |
3495 | if (g) |
3496 | { |
3497 | copy_ivec(SHIFT_IVEC(g, a1j)((g)->ishift[a1j]), jt1); |
3498 | ivec_sub(SHIFT_IVEC(g, a1i)((g)->ishift[a1i]), jt1, dt1_ij); |
3499 | ivec_sub(SHIFT_IVEC(g, a1k)((g)->ishift[a1k]), jt1, dt1_kj); |
3500 | ivec_sub(SHIFT_IVEC(g, a1l)((g)->ishift[a1l]), jt1, dt1_lj); |
3501 | t11 = IVEC2IS(dt1_ij)(((2*2 +1)*((2*1 +1)*(((dt1_ij)[2])+1)+((dt1_ij)[1])+1)+((dt1_ij )[0])+2)); |
3502 | t21 = IVEC2IS(dt1_kj)(((2*2 +1)*((2*1 +1)*(((dt1_kj)[2])+1)+((dt1_kj)[1])+1)+((dt1_kj )[0])+2)); |
3503 | t31 = IVEC2IS(dt1_lj)(((2*2 +1)*((2*1 +1)*(((dt1_lj)[2])+1)+((dt1_lj)[1])+1)+((dt1_lj )[0])+2)); |
3504 | |
3505 | copy_ivec(SHIFT_IVEC(g, a2j)((g)->ishift[a2j]), jt2); |
3506 | ivec_sub(SHIFT_IVEC(g, a2i)((g)->ishift[a2i]), jt2, dt2_ij); |
3507 | ivec_sub(SHIFT_IVEC(g, a2k)((g)->ishift[a2k]), jt2, dt2_kj); |
3508 | ivec_sub(SHIFT_IVEC(g, a2l)((g)->ishift[a2l]), jt2, dt2_lj); |
3509 | t12 = IVEC2IS(dt2_ij)(((2*2 +1)*((2*1 +1)*(((dt2_ij)[2])+1)+((dt2_ij)[1])+1)+((dt2_ij )[0])+2)); |
3510 | t22 = IVEC2IS(dt2_kj)(((2*2 +1)*((2*1 +1)*(((dt2_kj)[2])+1)+((dt2_kj)[1])+1)+((dt2_kj )[0])+2)); |
3511 | t32 = IVEC2IS(dt2_lj)(((2*2 +1)*((2*1 +1)*(((dt2_lj)[2])+1)+((dt2_lj)[1])+1)+((dt2_lj )[0])+2)); |
3512 | } |
3513 | else if (pbc) |
3514 | { |
3515 | t31 = pbc_rvec_sub(pbc, x[a1l], x[a1j], h1); |
3516 | t32 = pbc_rvec_sub(pbc, x[a2l], x[a2j], h2); |
3517 | } |
3518 | else |
3519 | { |
3520 | t31 = CENTRAL(((2*1 +1)*(2*1 +1)*(2*2 +1))/2); |
3521 | t32 = CENTRAL(((2*1 +1)*(2*1 +1)*(2*2 +1))/2); |
3522 | } |
3523 | |
3524 | rvec_inc(fshift[t11], f1_i); |
3525 | rvec_inc(fshift[CENTRAL(((2*1 +1)*(2*1 +1)*(2*2 +1))/2)], f1_j); |
3526 | rvec_inc(fshift[t21], f1_k); |
3527 | rvec_inc(fshift[t31], f1_l); |
3528 | |
3529 | rvec_inc(fshift[t21], f2_i); |
3530 | rvec_inc(fshift[CENTRAL(((2*1 +1)*(2*1 +1)*(2*2 +1))/2)], f2_j); |
3531 | rvec_inc(fshift[t22], f2_k); |
3532 | rvec_inc(fshift[t32], f2_l); |
3533 | } |
3534 | return vtot; |
3535 | } |
3536 | |
3537 | |
3538 | |
3539 | /*********************************************************** |
3540 | * |
3541 | * G R O M O S 9 6 F U N C T I O N S |
3542 | * |
3543 | ***********************************************************/ |
3544 | real g96harmonic(real kA, real kB, real xA, real xB, real x, real lambda, |
3545 | real *V, real *F) |
3546 | { |
3547 | const real half = 0.5; |
3548 | real L1, kk, x0, dx, dx2; |
3549 | real v, f, dvdlambda; |
3550 | |
3551 | L1 = 1.0-lambda; |
3552 | kk = L1*kA+lambda*kB; |
3553 | x0 = L1*xA+lambda*xB; |
3554 | |
3555 | dx = x-x0; |
3556 | dx2 = dx*dx; |
3557 | |
3558 | f = -kk*dx; |
3559 | v = half*kk*dx2; |
3560 | dvdlambda = half*(kB-kA)*dx2 + (xA-xB)*kk*dx; |
3561 | |
3562 | *F = f; |
3563 | *V = v; |
3564 | |
3565 | return dvdlambda; |
3566 | |
3567 | /* That was 21 flops */ |
3568 | } |
3569 | |
3570 | real g96bonds(int nbonds, |
3571 | const t_iatom forceatoms[], const t_iparams forceparams[], |
3572 | const rvec x[], rvec f[], rvec fshift[], |
3573 | const t_pbc *pbc, const t_graph *g, |
3574 | real lambda, real *dvdlambda, |
3575 | const t_mdatoms gmx_unused__attribute__ ((unused)) *md, t_fcdata gmx_unused__attribute__ ((unused)) *fcd, |
3576 | int gmx_unused__attribute__ ((unused)) *global_atom_index) |
3577 | { |
3578 | int i, m, ki, ai, aj, type; |
3579 | real dr2, fbond, vbond, fij, vtot; |
3580 | rvec dx; |
3581 | ivec dt; |
3582 | |
3583 | vtot = 0.0; |
3584 | for (i = 0; (i < nbonds); ) |
3585 | { |
3586 | type = forceatoms[i++]; |
3587 | ai = forceatoms[i++]; |
3588 | aj = forceatoms[i++]; |
3589 | |
3590 | ki = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /* 3 */ |
3591 | dr2 = iprod(dx, dx); /* 5 */ |
3592 | |
3593 | *dvdlambda += g96harmonic(forceparams[type].harmonic.krA, |
3594 | forceparams[type].harmonic.krB, |
3595 | forceparams[type].harmonic.rA, |
3596 | forceparams[type].harmonic.rB, |
3597 | dr2, lambda, &vbond, &fbond); |
3598 | |
3599 | vtot += 0.5*vbond; /* 1*/ |
3600 | #ifdef DEBUG |
3601 | if (debug) |
3602 | { |
3603 | fprintf(debug, "G96-BONDS: dr = %10g vbond = %10g fbond = %10g\n", |
3604 | sqrt(dr2), vbond, fbond); |
3605 | } |
3606 | #endif |
3607 | |
3608 | if (g) |
3609 | { |
3610 | ivec_sub(SHIFT_IVEC(g, ai)((g)->ishift[ai]), SHIFT_IVEC(g, aj)((g)->ishift[aj]), dt); |
3611 | ki = IVEC2IS(dt)(((2*2 +1)*((2*1 +1)*(((dt)[2])+1)+((dt)[1])+1)+((dt)[0])+2)); |
3612 | } |
3613 | for (m = 0; (m < DIM3); m++) /* 15 */ |
3614 | { |
3615 | fij = fbond*dx[m]; |
3616 | f[ai][m] += fij; |
3617 | f[aj][m] -= fij; |
3618 | fshift[ki][m] += fij; |
3619 | fshift[CENTRAL(((2*1 +1)*(2*1 +1)*(2*2 +1))/2)][m] -= fij; |
3620 | } |
3621 | } /* 44 TOTAL */ |
3622 | return vtot; |
3623 | } |
3624 | |
3625 | real g96bond_angle(const rvec xi, const rvec xj, const rvec xk, const t_pbc *pbc, |
3626 | rvec r_ij, rvec r_kj, |
3627 | int *t1, int *t2) |
3628 | /* Return value is the angle between the bonds i-j and j-k */ |
3629 | { |
3630 | real costh; |
3631 | |
3632 | *t1 = pbc_rvec_sub(pbc, xi, xj, r_ij); /* 3 */ |
3633 | *t2 = pbc_rvec_sub(pbc, xk, xj, r_kj); /* 3 */ |
3634 | |
3635 | costh = cos_angle(r_ij, r_kj); /* 25 */ |
3636 | /* 41 TOTAL */ |
3637 | return costh; |
3638 | } |
3639 | |
3640 | real g96angles(int nbonds, |
3641 | const t_iatom forceatoms[], const t_iparams forceparams[], |
3642 | const rvec x[], rvec f[], rvec fshift[], |
3643 | const t_pbc *pbc, const t_graph *g, |
3644 | real lambda, real *dvdlambda, |
3645 | const t_mdatoms gmx_unused__attribute__ ((unused)) *md, t_fcdata gmx_unused__attribute__ ((unused)) *fcd, |
3646 | int gmx_unused__attribute__ ((unused)) *global_atom_index) |
3647 | { |
3648 | int i, ai, aj, ak, type, m, t1, t2; |
3649 | rvec r_ij, r_kj; |
3650 | real cos_theta, dVdt, va, vtot; |
3651 | real rij_1, rij_2, rkj_1, rkj_2, rijrkj_1; |
3652 | rvec f_i, f_j, f_k; |
3653 | ivec jt, dt_ij, dt_kj; |
3654 | |
3655 | vtot = 0.0; |
3656 | for (i = 0; (i < nbonds); ) |
3657 | { |
3658 | type = forceatoms[i++]; |
3659 | ai = forceatoms[i++]; |
3660 | aj = forceatoms[i++]; |
3661 | ak = forceatoms[i++]; |
3662 | |
3663 | cos_theta = g96bond_angle(x[ai], x[aj], x[ak], pbc, r_ij, r_kj, &t1, &t2); |
3664 | |
3665 | *dvdlambda += g96harmonic(forceparams[type].harmonic.krA, |
3666 | forceparams[type].harmonic.krB, |
3667 | forceparams[type].harmonic.rA, |
3668 | forceparams[type].harmonic.rB, |
3669 | cos_theta, lambda, &va, &dVdt); |
3670 | vtot += va; |
3671 | |
3672 | rij_1 = gmx_invsqrt(iprod(r_ij, r_ij))gmx_software_invsqrt(iprod(r_ij, r_ij)); |
3673 | rkj_1 = gmx_invsqrt(iprod(r_kj, r_kj))gmx_software_invsqrt(iprod(r_kj, r_kj)); |
3674 | rij_2 = rij_1*rij_1; |
3675 | rkj_2 = rkj_1*rkj_1; |
3676 | rijrkj_1 = rij_1*rkj_1; /* 23 */ |
3677 | |
3678 | #ifdef DEBUG |
3679 | if (debug) |
3680 | { |
3681 | fprintf(debug, "G96ANGLES: costheta = %10g vth = %10g dV/dct = %10g\n", |
3682 | cos_theta, va, dVdt); |
3683 | } |
3684 | #endif |
3685 | for (m = 0; (m < DIM3); m++) /* 42 */ |
3686 | { |
3687 | f_i[m] = dVdt*(r_kj[m]*rijrkj_1 - r_ij[m]*rij_2*cos_theta); |
3688 | f_k[m] = dVdt*(r_ij[m]*rijrkj_1 - r_kj[m]*rkj_2*cos_theta); |
3689 | f_j[m] = -f_i[m]-f_k[m]; |
3690 | f[ai][m] += f_i[m]; |
3691 | f[aj][m] += f_j[m]; |
3692 | f[ak][m] += f_k[m]; |
3693 | } |
3694 | |
3695 | if (g) |
3696 | { |
3697 | copy_ivec(SHIFT_IVEC(g, aj)((g)->ishift[aj]), jt); |
3698 | |
3699 | ivec_sub(SHIFT_IVEC(g, ai)((g)->ishift[ai]), jt, dt_ij); |
3700 | ivec_sub(SHIFT_IVEC(g, ak)((g)->ishift[ak]), jt, dt_kj); |
3701 | t1 = IVEC2IS(dt_ij)(((2*2 +1)*((2*1 +1)*(((dt_ij)[2])+1)+((dt_ij)[1])+1)+((dt_ij )[0])+2)); |
3702 | t2 = IVEC2IS(dt_kj)(((2*2 +1)*((2*1 +1)*(((dt_kj)[2])+1)+((dt_kj)[1])+1)+((dt_kj )[0])+2)); |
3703 | } |
3704 | rvec_inc(fshift[t1], f_i); |
3705 | rvec_inc(fshift[CENTRAL(((2*1 +1)*(2*1 +1)*(2*2 +1))/2)], f_j); |
3706 | rvec_inc(fshift[t2], f_k); /* 9 */ |
3707 | /* 163 TOTAL */ |
3708 | } |
3709 | return vtot; |
3710 | } |
3711 | |
3712 | real cross_bond_bond(int nbonds, |
3713 | const t_iatom forceatoms[], const t_iparams forceparams[], |
3714 | const rvec x[], rvec f[], rvec fshift[], |
3715 | const t_pbc *pbc, const t_graph *g, |
3716 | real gmx_unused__attribute__ ((unused)) lambda, real gmx_unused__attribute__ ((unused)) *dvdlambda, |
3717 | const t_mdatoms gmx_unused__attribute__ ((unused)) *md, t_fcdata gmx_unused__attribute__ ((unused)) *fcd, |
3718 | int gmx_unused__attribute__ ((unused)) *global_atom_index) |
3719 | { |
3720 | /* Potential from Lawrence and Skimmer, Chem. Phys. Lett. 372 (2003) |
3721 | * pp. 842-847 |
3722 | */ |
3723 | int i, ai, aj, ak, type, m, t1, t2; |
3724 | rvec r_ij, r_kj; |
3725 | real vtot, vrr, s1, s2, r1, r2, r1e, r2e, krr; |
3726 | rvec f_i, f_j, f_k; |
3727 | ivec jt, dt_ij, dt_kj; |
3728 | |
3729 | vtot = 0.0; |
3730 | for (i = 0; (i < nbonds); ) |
3731 | { |
3732 | type = forceatoms[i++]; |
3733 | ai = forceatoms[i++]; |
3734 | aj = forceatoms[i++]; |
3735 | ak = forceatoms[i++]; |
3736 | r1e = forceparams[type].cross_bb.r1e; |
3737 | r2e = forceparams[type].cross_bb.r2e; |
3738 | krr = forceparams[type].cross_bb.krr; |
3739 | |
3740 | /* Compute distance vectors ... */ |
3741 | t1 = pbc_rvec_sub(pbc, x[ai], x[aj], r_ij); |
3742 | t2 = pbc_rvec_sub(pbc, x[ak], x[aj], r_kj); |
3743 | |
3744 | /* ... and their lengths */ |
3745 | r1 = norm(r_ij); |
3746 | r2 = norm(r_kj); |
3747 | |
3748 | /* Deviations from ideality */ |
3749 | s1 = r1-r1e; |
3750 | s2 = r2-r2e; |
3751 | |
3752 | /* Energy (can be negative!) */ |
3753 | vrr = krr*s1*s2; |
3754 | vtot += vrr; |
3755 | |
3756 | /* Forces */ |
3757 | svmul(-krr*s2/r1, r_ij, f_i); |
3758 | svmul(-krr*s1/r2, r_kj, f_k); |
3759 | |
3760 | for (m = 0; (m < DIM3); m++) /* 12 */ |
3761 | { |
3762 | f_j[m] = -f_i[m] - f_k[m]; |
3763 | f[ai][m] += f_i[m]; |
3764 | f[aj][m] += f_j[m]; |
3765 | f[ak][m] += f_k[m]; |
3766 | } |
3767 | |
3768 | /* Virial stuff */ |
3769 | if (g) |
3770 | { |
3771 | copy_ivec(SHIFT_IVEC(g, aj)((g)->ishift[aj]), jt); |
3772 | |
3773 | ivec_sub(SHIFT_IVEC(g, ai)((g)->ishift[ai]), jt, dt_ij); |
3774 | ivec_sub(SHIFT_IVEC(g, ak)((g)->ishift[ak]), jt, dt_kj); |
3775 | t1 = IVEC2IS(dt_ij)(((2*2 +1)*((2*1 +1)*(((dt_ij)[2])+1)+((dt_ij)[1])+1)+((dt_ij )[0])+2)); |
3776 | t2 = IVEC2IS(dt_kj)(((2*2 +1)*((2*1 +1)*(((dt_kj)[2])+1)+((dt_kj)[1])+1)+((dt_kj )[0])+2)); |
3777 | } |
3778 | rvec_inc(fshift[t1], f_i); |
3779 | rvec_inc(fshift[CENTRAL(((2*1 +1)*(2*1 +1)*(2*2 +1))/2)], f_j); |
3780 | rvec_inc(fshift[t2], f_k); /* 9 */ |
3781 | /* 163 TOTAL */ |
3782 | } |
3783 | return vtot; |
3784 | } |
3785 | |
3786 | real cross_bond_angle(int nbonds, |
3787 | const t_iatom forceatoms[], const t_iparams forceparams[], |
3788 | const rvec x[], rvec f[], rvec fshift[], |
3789 | const t_pbc *pbc, const t_graph *g, |
3790 | real gmx_unused__attribute__ ((unused)) lambda, real gmx_unused__attribute__ ((unused)) *dvdlambda, |
3791 | const t_mdatoms gmx_unused__attribute__ ((unused)) *md, t_fcdata gmx_unused__attribute__ ((unused)) *fcd, |
3792 | int gmx_unused__attribute__ ((unused)) *global_atom_index) |
3793 | { |
3794 | /* Potential from Lawrence and Skimmer, Chem. Phys. Lett. 372 (2003) |
3795 | * pp. 842-847 |
3796 | */ |
3797 | int i, ai, aj, ak, type, m, t1, t2, t3; |
3798 | rvec r_ij, r_kj, r_ik; |
3799 | real vtot, vrt, s1, s2, s3, r1, r2, r3, r1e, r2e, r3e, krt, k1, k2, k3; |
3800 | rvec f_i, f_j, f_k; |
3801 | ivec jt, dt_ij, dt_kj; |
3802 | |
3803 | vtot = 0.0; |
3804 | for (i = 0; (i < nbonds); ) |
3805 | { |
3806 | type = forceatoms[i++]; |
3807 | ai = forceatoms[i++]; |
3808 | aj = forceatoms[i++]; |
3809 | ak = forceatoms[i++]; |
3810 | r1e = forceparams[type].cross_ba.r1e; |
3811 | r2e = forceparams[type].cross_ba.r2e; |
3812 | r3e = forceparams[type].cross_ba.r3e; |
3813 | krt = forceparams[type].cross_ba.krt; |
3814 | |
3815 | /* Compute distance vectors ... */ |
3816 | t1 = pbc_rvec_sub(pbc, x[ai], x[aj], r_ij); |
3817 | t2 = pbc_rvec_sub(pbc, x[ak], x[aj], r_kj); |
3818 | t3 = pbc_rvec_sub(pbc, x[ai], x[ak], r_ik); |
Value stored to 't3' is never read | |
3819 | |
3820 | /* ... and their lengths */ |
3821 | r1 = norm(r_ij); |
3822 | r2 = norm(r_kj); |
3823 | r3 = norm(r_ik); |
3824 | |
3825 | /* Deviations from ideality */ |
3826 | s1 = r1-r1e; |
3827 | s2 = r2-r2e; |
3828 | s3 = r3-r3e; |
3829 | |
3830 | /* Energy (can be negative!) */ |
3831 | vrt = krt*s3*(s1+s2); |
3832 | vtot += vrt; |
3833 | |
3834 | /* Forces */ |
3835 | k1 = -krt*(s3/r1); |
3836 | k2 = -krt*(s3/r2); |
3837 | k3 = -krt*(s1+s2)/r3; |
3838 | for (m = 0; (m < DIM3); m++) |
3839 | { |
3840 | f_i[m] = k1*r_ij[m] + k3*r_ik[m]; |
3841 | f_k[m] = k2*r_kj[m] - k3*r_ik[m]; |
3842 | f_j[m] = -f_i[m] - f_k[m]; |
3843 | } |
3844 | |
3845 | for (m = 0; (m < DIM3); m++) /* 12 */ |
3846 | { |
3847 | f[ai][m] += f_i[m]; |
3848 | f[aj][m] += f_j[m]; |
3849 | f[ak][m] += f_k[m]; |
3850 | } |
3851 | |
3852 | /* Virial stuff */ |
3853 | if (g) |
3854 | { |
3855 | copy_ivec(SHIFT_IVEC(g, aj)((g)->ishift[aj]), jt); |
3856 | |
3857 | ivec_sub(SHIFT_IVEC(g, ai)((g)->ishift[ai]), jt, dt_ij); |
3858 | ivec_sub(SHIFT_IVEC(g, ak)((g)->ishift[ak]), jt, dt_kj); |
3859 | t1 = IVEC2IS(dt_ij)(((2*2 +1)*((2*1 +1)*(((dt_ij)[2])+1)+((dt_ij)[1])+1)+((dt_ij )[0])+2)); |
3860 | t2 = IVEC2IS(dt_kj)(((2*2 +1)*((2*1 +1)*(((dt_kj)[2])+1)+((dt_kj)[1])+1)+((dt_kj )[0])+2)); |
3861 | } |
3862 | rvec_inc(fshift[t1], f_i); |
3863 | rvec_inc(fshift[CENTRAL(((2*1 +1)*(2*1 +1)*(2*2 +1))/2)], f_j); |
3864 | rvec_inc(fshift[t2], f_k); /* 9 */ |
3865 | /* 163 TOTAL */ |
3866 | } |
3867 | return vtot; |
3868 | } |
3869 | |
3870 | static real bonded_tab(const char *type, int table_nr, |
3871 | const bondedtable_t *table, real kA, real kB, real r, |
3872 | real lambda, real *V, real *F) |
3873 | { |
3874 | real k, tabscale, *VFtab, rt, eps, eps2, Yt, Ft, Geps, Heps2, Fp, VV, FF; |
3875 | int n0, nnn; |
3876 | real v, f, dvdlambda; |
3877 | |
3878 | k = (1.0 - lambda)*kA + lambda*kB; |
3879 | |
3880 | tabscale = table->scale; |
3881 | VFtab = table->data; |
3882 | |
3883 | rt = r*tabscale; |
3884 | n0 = rt; |
3885 | if (n0 >= table->n) |
3886 | { |
3887 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxlib/bondfree.c" , 3887, "A tabulated %s interaction table number %d is out of the table range: r %f, between table indices %d and %d, table length %d", |
3888 | type, table_nr, r, n0, n0+1, table->n); |
3889 | } |
3890 | eps = rt - n0; |
3891 | eps2 = eps*eps; |
3892 | nnn = 4*n0; |
3893 | Yt = VFtab[nnn]; |
3894 | Ft = VFtab[nnn+1]; |
3895 | Geps = VFtab[nnn+2]*eps; |
3896 | Heps2 = VFtab[nnn+3]*eps2; |
3897 | Fp = Ft + Geps + Heps2; |
3898 | VV = Yt + Fp*eps; |
3899 | FF = Fp + Geps + 2.0*Heps2; |
3900 | |
3901 | *F = -k*FF*tabscale; |
3902 | *V = k*VV; |
3903 | dvdlambda = (kB - kA)*VV; |
3904 | |
3905 | return dvdlambda; |
3906 | |
3907 | /* That was 22 flops */ |
3908 | } |
3909 | |
3910 | real tab_bonds(int nbonds, |
3911 | const t_iatom forceatoms[], const t_iparams forceparams[], |
3912 | const rvec x[], rvec f[], rvec fshift[], |
3913 | const t_pbc *pbc, const t_graph *g, |
3914 | real lambda, real *dvdlambda, |
3915 | const t_mdatoms gmx_unused__attribute__ ((unused)) *md, t_fcdata *fcd, |
3916 | int gmx_unused__attribute__ ((unused)) *global_atom_index) |
3917 | { |
3918 | int i, m, ki, ai, aj, type, table; |
3919 | real dr, dr2, fbond, vbond, fij, vtot; |
3920 | rvec dx; |
3921 | ivec dt; |
3922 | |
3923 | vtot = 0.0; |
3924 | for (i = 0; (i < nbonds); ) |
3925 | { |
3926 | type = forceatoms[i++]; |
3927 | ai = forceatoms[i++]; |
3928 | aj = forceatoms[i++]; |
3929 | |
3930 | ki = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /* 3 */ |
3931 | dr2 = iprod(dx, dx); /* 5 */ |
3932 | dr = dr2*gmx_invsqrt(dr2)gmx_software_invsqrt(dr2); /* 10 */ |
3933 | |
3934 | table = forceparams[type].tab.table; |
3935 | |
3936 | *dvdlambda += bonded_tab("bond", table, |
3937 | &fcd->bondtab[table], |
3938 | forceparams[type].tab.kA, |
3939 | forceparams[type].tab.kB, |
3940 | dr, lambda, &vbond, &fbond); /* 22 */ |
3941 | |
3942 | if (dr2 == 0.0) |
3943 | { |
3944 | continue; |
3945 | } |
3946 | |
3947 | |
3948 | vtot += vbond; /* 1*/ |
3949 | fbond *= gmx_invsqrt(dr2)gmx_software_invsqrt(dr2); /* 6 */ |
3950 | #ifdef DEBUG |
3951 | if (debug) |
3952 | { |
3953 | fprintf(debug, "TABBONDS: dr = %10g vbond = %10g fbond = %10g\n", |
3954 | dr, vbond, fbond); |
3955 | } |
3956 | #endif |
3957 | if (g) |
3958 | { |
3959 | ivec_sub(SHIFT_IVEC(g, ai)((g)->ishift[ai]), SHIFT_IVEC(g, aj)((g)->ishift[aj]), dt); |
3960 | ki = IVEC2IS(dt)(((2*2 +1)*((2*1 +1)*(((dt)[2])+1)+((dt)[1])+1)+((dt)[0])+2)); |
3961 | } |
3962 | for (m = 0; (m < DIM3); m++) /* 15 */ |
3963 | { |
3964 | fij = fbond*dx[m]; |
3965 | f[ai][m] += fij; |
3966 | f[aj][m] -= fij; |
3967 | fshift[ki][m] += fij; |
3968 | fshift[CENTRAL(((2*1 +1)*(2*1 +1)*(2*2 +1))/2)][m] -= fij; |
3969 | } |
3970 | } /* 62 TOTAL */ |
3971 | return vtot; |
3972 | } |
3973 | |
3974 | real tab_angles(int nbonds, |
3975 | const t_iatom forceatoms[], const t_iparams forceparams[], |
3976 | const rvec x[], rvec f[], rvec fshift[], |
3977 | const t_pbc *pbc, const t_graph *g, |
3978 | real lambda, real *dvdlambda, |
3979 | const t_mdatoms gmx_unused__attribute__ ((unused)) *md, t_fcdata *fcd, |
3980 | int gmx_unused__attribute__ ((unused)) *global_atom_index) |
3981 | { |
3982 | int i, ai, aj, ak, t1, t2, type, table; |
3983 | rvec r_ij, r_kj; |
3984 | real cos_theta, cos_theta2, theta, dVdt, va, vtot; |
3985 | ivec jt, dt_ij, dt_kj; |
3986 | |
3987 | vtot = 0.0; |
3988 | for (i = 0; (i < nbonds); ) |
3989 | { |
3990 | type = forceatoms[i++]; |
3991 | ai = forceatoms[i++]; |
3992 | aj = forceatoms[i++]; |
3993 | ak = forceatoms[i++]; |
3994 | |
3995 | theta = bond_angle(x[ai], x[aj], x[ak], pbc, |
3996 | r_ij, r_kj, &cos_theta, &t1, &t2); /* 41 */ |
3997 | |
3998 | table = forceparams[type].tab.table; |
3999 | |
4000 | *dvdlambda += bonded_tab("angle", table, |
4001 | &fcd->angletab[table], |
4002 | forceparams[type].tab.kA, |
4003 | forceparams[type].tab.kB, |
4004 | theta, lambda, &va, &dVdt); /* 22 */ |
4005 | vtot += va; |
4006 | |
4007 | cos_theta2 = sqr(cos_theta); /* 1 */ |
4008 | if (cos_theta2 < 1) |
4009 | { |
4010 | int m; |
4011 | real snt, st, sth; |
4012 | real cik, cii, ckk; |
4013 | real nrkj2, nrij2; |
4014 | rvec f_i, f_j, f_k; |
4015 | |
4016 | st = dVdt*gmx_invsqrt(1 - cos_theta2)gmx_software_invsqrt(1 - cos_theta2); /* 12 */ |
4017 | sth = st*cos_theta; /* 1 */ |
4018 | #ifdef DEBUG |
4019 | if (debug) |
4020 | { |
4021 | fprintf(debug, "ANGLES: theta = %10g vth = %10g dV/dtheta = %10g\n", |
4022 | theta*RAD2DEG(180.0/3.14159265358979323846), va, dVdt); |
4023 | } |
4024 | #endif |
4025 | nrkj2 = iprod(r_kj, r_kj); /* 5 */ |
4026 | nrij2 = iprod(r_ij, r_ij); |
4027 | |
4028 | cik = st*gmx_invsqrt(nrkj2*nrij2)gmx_software_invsqrt(nrkj2*nrij2); /* 12 */ |
4029 | cii = sth/nrij2; /* 10 */ |
4030 | ckk = sth/nrkj2; /* 10 */ |
4031 | |
4032 | for (m = 0; (m < DIM3); m++) /* 39 */ |
4033 | { |
4034 | f_i[m] = -(cik*r_kj[m]-cii*r_ij[m]); |
4035 | f_k[m] = -(cik*r_ij[m]-ckk*r_kj[m]); |
4036 | f_j[m] = -f_i[m]-f_k[m]; |
4037 | f[ai][m] += f_i[m]; |
4038 | f[aj][m] += f_j[m]; |
4039 | f[ak][m] += f_k[m]; |
4040 | } |
4041 | if (g) |
4042 | { |
4043 | copy_ivec(SHIFT_IVEC(g, aj)((g)->ishift[aj]), jt); |
4044 | |
4045 | ivec_sub(SHIFT_IVEC(g, ai)((g)->ishift[ai]), jt, dt_ij); |
4046 | ivec_sub(SHIFT_IVEC(g, ak)((g)->ishift[ak]), jt, dt_kj); |
4047 | t1 = IVEC2IS(dt_ij)(((2*2 +1)*((2*1 +1)*(((dt_ij)[2])+1)+((dt_ij)[1])+1)+((dt_ij )[0])+2)); |
4048 | t2 = IVEC2IS(dt_kj)(((2*2 +1)*((2*1 +1)*(((dt_kj)[2])+1)+((dt_kj)[1])+1)+((dt_kj )[0])+2)); |
4049 | } |
4050 | rvec_inc(fshift[t1], f_i); |
4051 | rvec_inc(fshift[CENTRAL(((2*1 +1)*(2*1 +1)*(2*2 +1))/2)], f_j); |
4052 | rvec_inc(fshift[t2], f_k); |
4053 | } /* 169 TOTAL */ |
4054 | } |
4055 | return vtot; |
4056 | } |
4057 | |
4058 | real tab_dihs(int nbonds, |
4059 | const t_iatom forceatoms[], const t_iparams forceparams[], |
4060 | const rvec x[], rvec f[], rvec fshift[], |
4061 | const t_pbc *pbc, const t_graph *g, |
4062 | real lambda, real *dvdlambda, |
4063 | const t_mdatoms gmx_unused__attribute__ ((unused)) *md, t_fcdata *fcd, |
4064 | int gmx_unused__attribute__ ((unused)) *global_atom_index) |
4065 | { |
4066 | int i, type, ai, aj, ak, al, table; |
4067 | int t1, t2, t3; |
4068 | rvec r_ij, r_kj, r_kl, m, n; |
4069 | real phi, sign, ddphi, vpd, vtot; |
4070 | |
4071 | vtot = 0.0; |
4072 | for (i = 0; (i < nbonds); ) |
4073 | { |
4074 | type = forceatoms[i++]; |
4075 | ai = forceatoms[i++]; |
4076 | aj = forceatoms[i++]; |
4077 | ak = forceatoms[i++]; |
4078 | al = forceatoms[i++]; |
4079 | |
4080 | phi = dih_angle(x[ai], x[aj], x[ak], x[al], pbc, r_ij, r_kj, r_kl, m, n, |
4081 | &sign, &t1, &t2, &t3); /* 84 */ |
4082 | |
4083 | table = forceparams[type].tab.table; |
4084 | |
4085 | /* Hopefully phi+M_PI never results in values < 0 */ |
4086 | *dvdlambda += bonded_tab("dihedral", table, |
4087 | &fcd->dihtab[table], |
4088 | forceparams[type].tab.kA, |
4089 | forceparams[type].tab.kB, |
4090 | phi+M_PI3.14159265358979323846, lambda, &vpd, &ddphi); |
4091 | |
4092 | vtot += vpd; |
4093 | do_dih_fup(ai, aj, ak, al, -ddphi, r_ij, r_kj, r_kl, m, n, |
4094 | f, fshift, pbc, g, x, t1, t2, t3); /* 112 */ |
4095 | |
4096 | #ifdef DEBUG |
4097 | fprintf(debug, "pdih: (%d,%d,%d,%d) phi=%g\n", |
4098 | ai, aj, ak, al, phi); |
4099 | #endif |
4100 | } /* 227 TOTAL */ |
4101 | |
4102 | return vtot; |
4103 | } |
4104 | |
4105 | /* Return if this is a potential calculated in bondfree.c, |
4106 | * i.e. an interaction that actually calculates a potential and |
4107 | * works on multiple atoms (not e.g. a connection or a position restraint). |
4108 | */ |
4109 | static gmx_inlineinline gmx_bool ftype_is_bonded_potential(int ftype) |
4110 | { |
4111 | return |
4112 | (interaction_function[ftype].flags & IF_BOND1) && |
4113 | !(ftype == F_CONNBONDS || ftype == F_POSRES || ftype == F_FBPOSRES) && |
4114 | (ftype < F_GB12 || ftype > F_GB14); |
4115 | } |
4116 | |
4117 | static void divide_bondeds_over_threads(t_idef *idef, int nthreads) |
4118 | { |
4119 | int ftype; |
4120 | int nat1; |
4121 | int t; |
4122 | int il_nr_thread; |
4123 | |
4124 | idef->nthreads = nthreads; |
4125 | |
4126 | if (F_NRE*(nthreads+1) > idef->il_thread_division_nalloc) |
4127 | { |
4128 | idef->il_thread_division_nalloc = F_NRE*(nthreads+1); |
4129 | snew(idef->il_thread_division, idef->il_thread_division_nalloc)(idef->il_thread_division) = save_calloc("idef->il_thread_division" , "/home/alexxy/Develop/gromacs/src/gromacs/gmxlib/bondfree.c" , 4129, (idef->il_thread_division_nalloc), sizeof(*(idef-> il_thread_division))); |
4130 | } |
4131 | |
4132 | for (ftype = 0; ftype < F_NRE; ftype++) |
4133 | { |
4134 | if (ftype_is_bonded_potential(ftype)) |
4135 | { |
4136 | nat1 = interaction_function[ftype].nratoms + 1; |
4137 | |
4138 | for (t = 0; t <= nthreads; t++) |
4139 | { |
4140 | /* Divide the interactions equally over the threads. |
4141 | * When the different types of bonded interactions |
4142 | * are distributed roughly equally over the threads, |
4143 | * this should lead to well localized output into |
4144 | * the force buffer on each thread. |
4145 | * If this is not the case, a more advanced scheme |
4146 | * (not implemented yet) will do better. |
4147 | */ |
4148 | il_nr_thread = (((idef->il[ftype].nr/nat1)*t)/nthreads)*nat1; |
4149 | |
4150 | /* Ensure that distance restraint pairs with the same label |
4151 | * end up on the same thread. |
4152 | * This is slighlty tricky code, since the next for iteration |
4153 | * may have an initial il_nr_thread lower than the final value |
4154 | * in the previous iteration, but this will anyhow be increased |
4155 | * to the approriate value again by this while loop. |
4156 | */ |
4157 | while (ftype == F_DISRES && |
4158 | il_nr_thread > 0 && |
4159 | il_nr_thread < idef->il[ftype].nr && |
4160 | idef->iparams[idef->il[ftype].iatoms[il_nr_thread]].disres.label == |
4161 | idef->iparams[idef->il[ftype].iatoms[il_nr_thread-nat1]].disres.label) |
4162 | { |
4163 | il_nr_thread += nat1; |
4164 | } |
4165 | |
4166 | idef->il_thread_division[ftype*(nthreads+1)+t] = il_nr_thread; |
4167 | } |
4168 | } |
4169 | } |
4170 | } |
4171 | |
4172 | static unsigned |
4173 | calc_bonded_reduction_mask(const t_idef *idef, |
4174 | int shift, |
4175 | int t, int nt) |
4176 | { |
4177 | unsigned mask; |
4178 | int ftype, nb, nat1, nb0, nb1, i, a; |
4179 | |
4180 | mask = 0; |
4181 | |
4182 | for (ftype = 0; ftype < F_NRE; ftype++) |
4183 | { |
4184 | if (ftype_is_bonded_potential(ftype)) |
4185 | { |
4186 | nb = idef->il[ftype].nr; |
4187 | if (nb > 0) |
4188 | { |
4189 | nat1 = interaction_function[ftype].nratoms + 1; |
4190 | |
4191 | /* Divide this interaction equally over the threads. |
4192 | * This is not stored: should match division in calc_bonds. |
4193 | */ |
4194 | nb0 = idef->il_thread_division[ftype*(nt+1)+t]; |
4195 | nb1 = idef->il_thread_division[ftype*(nt+1)+t+1]; |
4196 | |
4197 | for (i = nb0; i < nb1; i += nat1) |
4198 | { |
4199 | for (a = 1; a < nat1; a++) |
4200 | { |
4201 | mask |= (1U << (idef->il[ftype].iatoms[i+a]>>shift)); |
4202 | } |
4203 | } |
4204 | } |
4205 | } |
4206 | } |
4207 | |
4208 | return mask; |
4209 | } |
4210 | |
4211 | void setup_bonded_threading(t_forcerec *fr, t_idef *idef) |
4212 | { |
4213 | #define MAX_BLOCK_BITS32 32 |
4214 | int t; |
4215 | int ctot, c, b; |
4216 | |
4217 | assert(fr->nthreads >= 1)((void) (0)); |
4218 | |
4219 | /* Divide the bonded interaction over the threads */ |
4220 | divide_bondeds_over_threads(idef, fr->nthreads); |
4221 | |
4222 | if (fr->nthreads == 1) |
4223 | { |
4224 | fr->red_nblock = 0; |
4225 | |
4226 | return; |
4227 | } |
4228 | |
4229 | /* We divide the force array in a maximum of 32 blocks. |
4230 | * Minimum force block reduction size is 2^6=64. |
4231 | */ |
4232 | fr->red_ashift = 6; |
4233 | while (fr->natoms_force > (int)(MAX_BLOCK_BITS32*(1U<<fr->red_ashift))) |
4234 | { |
4235 | fr->red_ashift++; |
4236 | } |
4237 | if (debug) |
4238 | { |
4239 | fprintf(debug, "bonded force buffer block atom shift %d bits\n", |
4240 | fr->red_ashift); |
4241 | } |
4242 | |
4243 | /* Determine to which blocks each thread's bonded force calculation |
4244 | * contributes. Store this is a mask for each thread. |
4245 | */ |
4246 | #pragma omp parallel for num_threads(fr->nthreads) schedule(static) |
4247 | for (t = 1; t < fr->nthreads; t++) |
4248 | { |
4249 | fr->f_t[t].red_mask = |
4250 | calc_bonded_reduction_mask(idef, fr->red_ashift, t, fr->nthreads); |
4251 | } |
4252 | |
4253 | /* Determine the maximum number of blocks we need to reduce over */ |
4254 | fr->red_nblock = 0; |
4255 | ctot = 0; |
4256 | for (t = 0; t < fr->nthreads; t++) |
4257 | { |
4258 | c = 0; |
4259 | for (b = 0; b < MAX_BLOCK_BITS32; b++) |
4260 | { |
4261 | if (fr->f_t[t].red_mask & (1U<<b)) |
4262 | { |
4263 | fr->red_nblock = max(fr->red_nblock, b+1)(((fr->red_nblock) > (b+1)) ? (fr->red_nblock) : (b+ 1) ); |
4264 | c++; |
4265 | } |
4266 | } |
4267 | if (debug) |
4268 | { |
4269 | fprintf(debug, "thread %d flags %x count %d\n", |
4270 | t, fr->f_t[t].red_mask, c); |
4271 | } |
4272 | ctot += c; |
4273 | } |
4274 | if (debug) |
4275 | { |
4276 | fprintf(debug, "Number of blocks to reduce: %d of size %d\n", |
4277 | fr->red_nblock, 1<<fr->red_ashift); |
4278 | fprintf(debug, "Reduction density %.2f density/#thread %.2f\n", |
4279 | ctot*(1<<fr->red_ashift)/(double)fr->natoms_force, |
4280 | ctot*(1<<fr->red_ashift)/(double)(fr->natoms_force*fr->nthreads)); |
4281 | } |
4282 | } |
4283 | |
4284 | static void zero_thread_forces(f_thread_t *f_t, int n, |
4285 | int nblock, int blocksize) |
4286 | { |
4287 | int b, a0, a1, a, i, j; |
4288 | |
4289 | if (n > f_t->f_nalloc) |
4290 | { |
4291 | f_t->f_nalloc = over_alloc_large(n)(int)(1.19*(n) + 1000); |
4292 | srenew(f_t->f, f_t->f_nalloc)(f_t->f) = save_realloc("f_t->f", "/home/alexxy/Develop/gromacs/src/gromacs/gmxlib/bondfree.c" , 4292, (f_t->f), (f_t->f_nalloc), sizeof(*(f_t->f)) ); |
4293 | } |
4294 | |
4295 | if (f_t->red_mask != 0) |
4296 | { |
4297 | for (b = 0; b < nblock; b++) |
4298 | { |
4299 | if (f_t->red_mask && (1U<<b)) |
4300 | { |
4301 | a0 = b*blocksize; |
4302 | a1 = min((b+1)*blocksize, n)((((b+1)*blocksize) < (n)) ? ((b+1)*blocksize) : (n) ); |
4303 | for (a = a0; a < a1; a++) |
4304 | { |
4305 | clear_rvec(f_t->f[a]); |
4306 | } |
4307 | } |
4308 | } |
4309 | } |
4310 | for (i = 0; i < SHIFTS((2*1 +1)*(2*1 +1)*(2*2 +1)); i++) |
4311 | { |
4312 | clear_rvec(f_t->fshift[i]); |
4313 | } |
4314 | for (i = 0; i < F_NRE; i++) |
4315 | { |
4316 | f_t->ener[i] = 0; |
4317 | } |
4318 | for (i = 0; i < egNR; i++) |
4319 | { |
4320 | for (j = 0; j < f_t->grpp.nener; j++) |
4321 | { |
4322 | f_t->grpp.ener[i][j] = 0; |
4323 | } |
4324 | } |
4325 | for (i = 0; i < efptNR; i++) |
4326 | { |
4327 | f_t->dvdl[i] = 0; |
4328 | } |
4329 | } |
4330 | |
4331 | static void reduce_thread_force_buffer(int n, rvec *f, |
4332 | int nthreads, f_thread_t *f_t, |
4333 | int nblock, int block_size) |
4334 | { |
4335 | /* The max thread number is arbitrary, |
4336 | * we used a fixed number to avoid memory management. |
4337 | * Using more than 16 threads is probably never useful performance wise. |
4338 | */ |
4339 | #define MAX_BONDED_THREADS256 256 |
4340 | int b; |
4341 | |
4342 | if (nthreads > MAX_BONDED_THREADS256) |
4343 | { |
4344 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxlib/bondfree.c" , 4344, "Can not reduce bonded forces on more than %d threads", |
4345 | MAX_BONDED_THREADS256); |
4346 | } |
4347 | |
4348 | /* This reduction can run on any number of threads, |
4349 | * independently of nthreads. |
4350 | */ |
4351 | #pragma omp parallel for num_threads(nthreads) schedule(static) |
4352 | for (b = 0; b < nblock; b++) |
4353 | { |
4354 | rvec *fp[MAX_BONDED_THREADS256]; |
4355 | int nfb, ft, fb; |
4356 | int a0, a1, a; |
4357 | |
4358 | /* Determine which threads contribute to this block */ |
4359 | nfb = 0; |
4360 | for (ft = 1; ft < nthreads; ft++) |
4361 | { |
4362 | if (f_t[ft].red_mask & (1U<<b)) |
4363 | { |
4364 | fp[nfb++] = f_t[ft].f; |
4365 | } |
4366 | } |
4367 | if (nfb > 0) |
4368 | { |
4369 | /* Reduce force buffers for threads that contribute */ |
4370 | a0 = b *block_size; |
4371 | a1 = (b+1)*block_size; |
4372 | a1 = min(a1, n)(((a1) < (n)) ? (a1) : (n) ); |
4373 | for (a = a0; a < a1; a++) |
4374 | { |
4375 | for (fb = 0; fb < nfb; fb++) |
4376 | { |
4377 | rvec_inc(f[a], fp[fb][a]); |
4378 | } |
4379 | } |
4380 | } |
4381 | } |
4382 | } |
4383 | |
4384 | static void reduce_thread_forces(int n, rvec *f, rvec *fshift, |
4385 | real *ener, gmx_grppairener_t *grpp, real *dvdl, |
4386 | int nthreads, f_thread_t *f_t, |
4387 | int nblock, int block_size, |
4388 | gmx_bool bCalcEnerVir, |
4389 | gmx_bool bDHDL) |
4390 | { |
4391 | if (nblock > 0) |
4392 | { |
4393 | /* Reduce the bonded force buffer */ |
4394 | reduce_thread_force_buffer(n, f, nthreads, f_t, nblock, block_size); |
4395 | } |
4396 | |
4397 | /* When necessary, reduce energy and virial using one thread only */ |
4398 | if (bCalcEnerVir) |
4399 | { |
4400 | int t, i, j; |
4401 | |
4402 | for (i = 0; i < SHIFTS((2*1 +1)*(2*1 +1)*(2*2 +1)); i++) |
4403 | { |
4404 | for (t = 1; t < nthreads; t++) |
4405 | { |
4406 | rvec_inc(fshift[i], f_t[t].fshift[i]); |
4407 | } |
4408 | } |
4409 | for (i = 0; i < F_NRE; i++) |
4410 | { |
4411 | for (t = 1; t < nthreads; t++) |
4412 | { |
4413 | ener[i] += f_t[t].ener[i]; |
4414 | } |
4415 | } |
4416 | for (i = 0; i < egNR; i++) |
4417 | { |
4418 | for (j = 0; j < f_t[1].grpp.nener; j++) |
4419 | { |
4420 | for (t = 1; t < nthreads; t++) |
4421 | { |
4422 | |
4423 | grpp->ener[i][j] += f_t[t].grpp.ener[i][j]; |
4424 | } |
4425 | } |
4426 | } |
4427 | if (bDHDL) |
4428 | { |
4429 | for (i = 0; i < efptNR; i++) |
4430 | { |
4431 | |
4432 | for (t = 1; t < nthreads; t++) |
4433 | { |
4434 | dvdl[i] += f_t[t].dvdl[i]; |
4435 | } |
4436 | } |
4437 | } |
4438 | } |
4439 | } |
4440 | |
4441 | static real calc_one_bond(FILE *fplog, int thread, |
4442 | int ftype, const t_idef *idef, |
4443 | rvec x[], rvec f[], rvec fshift[], |
4444 | t_forcerec *fr, |
4445 | const t_pbc *pbc, const t_graph *g, |
4446 | gmx_grppairener_t *grpp, |
4447 | t_nrnb *nrnb, |
4448 | real *lambda, real *dvdl, |
4449 | const t_mdatoms *md, t_fcdata *fcd, |
4450 | gmx_bool bCalcEnerVir, |
4451 | int *global_atom_index, gmx_bool bPrintSepPot) |
4452 | { |
4453 | int nat1, nbonds, efptFTYPE; |
4454 | real v = 0; |
4455 | t_iatom *iatoms; |
4456 | int nb0, nbn; |
4457 | |
4458 | if (IS_RESTRAINT_TYPE(ftype)(((ftype == F_POSRES) || (ftype == F_FBPOSRES) || (ftype == F_DISRES ) || (ftype == F_RESTRBONDS) || (ftype == F_DISRESVIOL) || (ftype == F_ORIRES) || (ftype == F_ORIRESDEV) || (ftype == F_ANGRES ) || (ftype == F_ANGRESZ) || (ftype == F_DIHRES)))) |
4459 | { |
4460 | efptFTYPE = efptRESTRAINT; |
4461 | } |
4462 | else |
4463 | { |
4464 | efptFTYPE = efptBONDED; |
4465 | } |
4466 | |
4467 | nat1 = interaction_function[ftype].nratoms + 1; |
4468 | nbonds = idef->il[ftype].nr/nat1; |
4469 | iatoms = idef->il[ftype].iatoms; |
4470 | |
4471 | nb0 = idef->il_thread_division[ftype*(idef->nthreads+1)+thread]; |
4472 | nbn = idef->il_thread_division[ftype*(idef->nthreads+1)+thread+1] - nb0; |
4473 | |
4474 | if (!IS_LISTED_LJ_C(ftype)((ftype) >= F_LJ14 && (ftype) <= F_LJC_PAIRS_NB )) |
4475 | { |
4476 | if (ftype == F_CMAP) |
4477 | { |
4478 | v = cmap_dihs(nbn, iatoms+nb0, |
4479 | idef->iparams, &idef->cmap_grid, |
4480 | (const rvec*)x, f, fshift, |
4481 | pbc, g, lambda[efptFTYPE], &(dvdl[efptFTYPE]), |
4482 | md, fcd, global_atom_index); |
4483 | } |
4484 | #ifdef GMX_SIMD_HAVE_REAL |
4485 | else if (ftype == F_ANGLES && |
4486 | !bCalcEnerVir && fr->efep == efepNO) |
4487 | { |
4488 | /* No energies, shift forces, dvdl */ |
4489 | angles_noener_simd(nbn, idef->il[ftype].iatoms+nb0, |
4490 | idef->iparams, |
4491 | (const rvec*)x, f, |
4492 | pbc, g, lambda[efptFTYPE], md, fcd, |
4493 | global_atom_index); |
4494 | v = 0; |
4495 | } |
4496 | #endif |
4497 | else if (ftype == F_PDIHS && |
4498 | !bCalcEnerVir && fr->efep == efepNO) |
4499 | { |
4500 | /* No energies, shift forces, dvdl */ |
4501 | #ifdef GMX_SIMD_HAVE_REAL |
4502 | pdihs_noener_simd |
4503 | #else |
4504 | pdihs_noener |
4505 | #endif |
4506 | (nbn, idef->il[ftype].iatoms+nb0, |
4507 | idef->iparams, |
4508 | (const rvec*)x, f, |
4509 | pbc, g, lambda[efptFTYPE], md, fcd, |
4510 | global_atom_index); |
4511 | v = 0; |
4512 | } |
4513 | else |
4514 | { |
4515 | v = interaction_function[ftype].ifunc(nbn, iatoms+nb0, |
4516 | idef->iparams, |
4517 | (const rvec*)x, f, fshift, |
4518 | pbc, g, lambda[efptFTYPE], &(dvdl[efptFTYPE]), |
4519 | md, fcd, global_atom_index); |
4520 | } |
4521 | if (bPrintSepPot) |
4522 | { |
4523 | fprintf(fplog, " %-23s #%4d V %12.5e dVdl %12.5e\n", |
4524 | interaction_function[ftype].longname, |
4525 | nbonds, v, lambda[efptFTYPE]); |
4526 | } |
4527 | } |
4528 | else |
4529 | { |
4530 | v = do_nonbonded_listed(ftype, nbn, iatoms+nb0, idef->iparams, (const rvec*)x, f, fshift, |
4531 | pbc, g, lambda, dvdl, md, fr, grpp, global_atom_index); |
4532 | |
4533 | if (bPrintSepPot) |
4534 | { |
4535 | fprintf(fplog, " %-5s + %-15s #%4d dVdl %12.5e\n", |
4536 | interaction_function[ftype].longname, |
4537 | interaction_function[F_LJ14].longname, nbonds, dvdl[efptVDW]); |
4538 | fprintf(fplog, " %-5s + %-15s #%4d dVdl %12.5e\n", |
4539 | interaction_function[ftype].longname, |
4540 | interaction_function[F_COUL14].longname, nbonds, dvdl[efptCOUL]); |
4541 | } |
4542 | } |
4543 | |
4544 | if (thread == 0) |
4545 | { |
4546 | inc_nrnb(nrnb, interaction_function[ftype].nrnb_ind, nbonds)(nrnb)->n[interaction_function[ftype].nrnb_ind] += nbonds; |
4547 | } |
4548 | |
4549 | return v; |
4550 | } |
4551 | |
4552 | void calc_bonds(FILE *fplog, const gmx_multisim_t *ms, |
4553 | const t_idef *idef, |
4554 | rvec x[], history_t *hist, |
4555 | rvec f[], t_forcerec *fr, |
4556 | const t_pbc *pbc, const t_graph *g, |
4557 | gmx_enerdata_t *enerd, t_nrnb *nrnb, |
4558 | real *lambda, |
4559 | const t_mdatoms *md, |
4560 | t_fcdata *fcd, int *global_atom_index, |
4561 | t_atomtypes gmx_unused__attribute__ ((unused)) *atype, gmx_genborn_t gmx_unused__attribute__ ((unused)) *born, |
4562 | int force_flags, |
4563 | gmx_bool bPrintSepPot, gmx_int64_t step) |
4564 | { |
4565 | gmx_bool bCalcEnerVir; |
4566 | int i; |
4567 | real v, dvdl[efptNR], dvdl_dum[efptNR]; /* The dummy array is to have a place to store the dhdl at other values |
4568 | of lambda, which will be thrown away in the end*/ |
4569 | const t_pbc *pbc_null; |
4570 | char buf[22]; |
4571 | int thread; |
4572 | |
4573 | assert(fr->nthreads == idef->nthreads)((void) (0)); |
4574 | |
4575 | bCalcEnerVir = (force_flags & (GMX_FORCE_VIRIAL(1<<8) | GMX_FORCE_ENERGY(1<<9))); |
4576 | |
4577 | for (i = 0; i < efptNR; i++) |
4578 | { |
4579 | dvdl[i] = 0.0; |
4580 | } |
4581 | if (fr->bMolPBC) |
4582 | { |
4583 | pbc_null = pbc; |
4584 | } |
4585 | else |
4586 | { |
4587 | pbc_null = NULL((void*)0); |
4588 | } |
4589 | if (bPrintSepPot) |
4590 | { |
4591 | fprintf(fplog, "Step %s: bonded V and dVdl for this node\n", |
4592 | gmx_step_str(step, buf)); |
4593 | } |
4594 | |
4595 | #ifdef DEBUG |
4596 | if (g && debug) |
4597 | { |
4598 | p_graph(debug, "Bondage is fun", g); |
4599 | } |
4600 | #endif |
4601 | |
4602 | /* Do pre force calculation stuff which might require communication */ |
4603 | if (idef->il[F_ORIRES].nr) |
4604 | { |
4605 | enerd->term[F_ORIRESDEV] = |
4606 | calc_orires_dev(ms, idef->il[F_ORIRES].nr, |
4607 | idef->il[F_ORIRES].iatoms, |
4608 | idef->iparams, md, (const rvec*)x, |
4609 | pbc_null, fcd, hist); |
4610 | } |
4611 | if (idef->il[F_DISRES].nr) |
4612 | { |
4613 | calc_disres_R_6(idef->il[F_DISRES].nr, |
4614 | idef->il[F_DISRES].iatoms, |
4615 | idef->iparams, (const rvec*)x, pbc_null, |
4616 | fcd, hist); |
4617 | #ifdef GMX_MPI |
4618 | if (fcd->disres.nsystems > 1) |
4619 | { |
4620 | gmx_sum_simgmx_sumf_sim(2*fcd->disres.nres, fcd->disres.Rt_6, ms); |
4621 | } |
4622 | #endif |
4623 | } |
4624 | |
4625 | #pragma omp parallel for num_threads(fr->nthreads) schedule(static) |
4626 | for (thread = 0; thread < fr->nthreads; thread++) |
4627 | { |
4628 | int ftype; |
4629 | real *epot, v; |
4630 | /* thread stuff */ |
4631 | rvec *ft, *fshift; |
4632 | real *dvdlt; |
4633 | gmx_grppairener_t *grpp; |
4634 | |
4635 | if (thread == 0) |
4636 | { |
4637 | ft = f; |
4638 | fshift = fr->fshift; |
4639 | epot = enerd->term; |
4640 | grpp = &enerd->grpp; |
4641 | dvdlt = dvdl; |
4642 | } |
4643 | else |
4644 | { |
4645 | zero_thread_forces(&fr->f_t[thread], fr->natoms_force, |
4646 | fr->red_nblock, 1<<fr->red_ashift); |
4647 | |
4648 | ft = fr->f_t[thread].f; |
4649 | fshift = fr->f_t[thread].fshift; |
4650 | epot = fr->f_t[thread].ener; |
4651 | grpp = &fr->f_t[thread].grpp; |
4652 | dvdlt = fr->f_t[thread].dvdl; |
4653 | } |
4654 | /* Loop over all bonded force types to calculate the bonded forces */ |
4655 | for (ftype = 0; (ftype < F_NRE); ftype++) |
4656 | { |
4657 | if (idef->il[ftype].nr > 0 && ftype_is_bonded_potential(ftype)) |
4658 | { |
4659 | v = calc_one_bond(fplog, thread, ftype, idef, x, |
4660 | ft, fshift, fr, pbc_null, g, grpp, |
4661 | nrnb, lambda, dvdlt, |
4662 | md, fcd, bCalcEnerVir, |
4663 | global_atom_index, bPrintSepPot); |
4664 | epot[ftype] += v; |
4665 | } |
4666 | } |
4667 | } |
4668 | if (fr->nthreads > 1) |
4669 | { |
4670 | reduce_thread_forces(fr->natoms_force, f, fr->fshift, |
4671 | enerd->term, &enerd->grpp, dvdl, |
4672 | fr->nthreads, fr->f_t, |
4673 | fr->red_nblock, 1<<fr->red_ashift, |
4674 | bCalcEnerVir, |
4675 | force_flags & GMX_FORCE_DHDL(1<<10)); |
4676 | } |
4677 | if (force_flags & GMX_FORCE_DHDL(1<<10)) |
4678 | { |
4679 | for (i = 0; i < efptNR; i++) |
4680 | { |
4681 | enerd->dvdl_nonlin[i] += dvdl[i]; |
4682 | } |
4683 | } |
4684 | |
4685 | /* Copy the sum of violations for the distance restraints from fcd */ |
4686 | if (fcd) |
4687 | { |
4688 | enerd->term[F_DISRESVIOL] = fcd->disres.sumviol; |
4689 | |
4690 | } |
4691 | } |
4692 | |
4693 | void calc_bonds_lambda(FILE *fplog, |
4694 | const t_idef *idef, |
4695 | rvec x[], |
4696 | t_forcerec *fr, |
4697 | const t_pbc *pbc, const t_graph *g, |
4698 | gmx_grppairener_t *grpp, real *epot, t_nrnb *nrnb, |
4699 | real *lambda, |
4700 | const t_mdatoms *md, |
4701 | t_fcdata *fcd, |
4702 | int *global_atom_index) |
4703 | { |
4704 | int i, ftype, nr_nonperturbed, nr; |
4705 | real v; |
4706 | real dvdl_dum[efptNR]; |
4707 | rvec *f, *fshift; |
4708 | const t_pbc *pbc_null; |
4709 | t_idef idef_fe; |
4710 | |
4711 | if (fr->bMolPBC) |
4712 | { |
4713 | pbc_null = pbc; |
4714 | } |
4715 | else |
4716 | { |
4717 | pbc_null = NULL((void*)0); |
4718 | } |
4719 | |
4720 | /* Copy the whole idef, so we can modify the contents locally */ |
4721 | idef_fe = *idef; |
4722 | idef_fe.nthreads = 1; |
4723 | snew(idef_fe.il_thread_division, F_NRE*(idef_fe.nthreads+1))(idef_fe.il_thread_division) = save_calloc("idef_fe.il_thread_division" , "/home/alexxy/Develop/gromacs/src/gromacs/gmxlib/bondfree.c" , 4723, (F_NRE*(idef_fe.nthreads+1)), sizeof(*(idef_fe.il_thread_division ))); |
4724 | |
4725 | /* We already have the forces, so we use temp buffers here */ |
4726 | snew(f, fr->natoms_force)(f) = save_calloc("f", "/home/alexxy/Develop/gromacs/src/gromacs/gmxlib/bondfree.c" , 4726, (fr->natoms_force), sizeof(*(f))); |
4727 | snew(fshift, SHIFTS)(fshift) = save_calloc("fshift", "/home/alexxy/Develop/gromacs/src/gromacs/gmxlib/bondfree.c" , 4727, (((2*1 +1)*(2*1 +1)*(2*2 +1))), sizeof(*(fshift))); |
4728 | |
4729 | /* Loop over all bonded force types to calculate the bonded energies */ |
4730 | for (ftype = 0; (ftype < F_NRE); ftype++) |
4731 | { |
4732 | if (ftype_is_bonded_potential(ftype)) |
4733 | { |
4734 | /* Set the work range of thread 0 to the perturbed bondeds only */ |
4735 | nr_nonperturbed = idef->il[ftype].nr_nonperturbed; |
4736 | nr = idef->il[ftype].nr; |
4737 | idef_fe.il_thread_division[ftype*2+0] = nr_nonperturbed; |
4738 | idef_fe.il_thread_division[ftype*2+1] = nr; |
4739 | |
4740 | /* This is only to get the flop count correct */ |
4741 | idef_fe.il[ftype].nr = nr - nr_nonperturbed; |
4742 | |
4743 | if (nr - nr_nonperturbed > 0) |
4744 | { |
4745 | v = calc_one_bond(fplog, 0, ftype, &idef_fe, |
4746 | x, f, fshift, fr, pbc_null, g, |
4747 | grpp, nrnb, lambda, dvdl_dum, |
4748 | md, fcd, TRUE1, |
4749 | global_atom_index, FALSE0); |
4750 | epot[ftype] += v; |
4751 | } |
4752 | } |
4753 | } |
4754 | |
4755 | sfree(fshift)save_free("fshift", "/home/alexxy/Develop/gromacs/src/gromacs/gmxlib/bondfree.c" , 4755, (fshift)); |
4756 | sfree(f)save_free("f", "/home/alexxy/Develop/gromacs/src/gromacs/gmxlib/bondfree.c" , 4756, (f)); |
4757 | |
4758 | sfree(idef_fe.il_thread_division)save_free("idef_fe.il_thread_division", "/home/alexxy/Develop/gromacs/src/gromacs/gmxlib/bondfree.c" , 4758, (idef_fe.il_thread_division)); |
4759 | } |