1 | |
2 | |
3 | |
4 | |
5 | |
6 | |
7 | |
8 | |
9 | |
10 | |
11 | |
12 | |
13 | |
14 | |
15 | |
16 | |
17 | |
18 | |
19 | |
20 | |
21 | |
22 | |
23 | |
24 | |
25 | |
26 | |
27 | |
28 | |
29 | |
30 | |
31 | |
32 | |
33 | |
34 | |
35 | |
36 | |
37 | |
38 | #include <math.h> |
39 | #include <assert.h> |
40 | #include "physics.h" |
41 | #include "gromacs/math/vec.h" |
42 | #include "gromacs/math/utilities.h" |
43 | #include "txtdump.h" |
44 | #include "bondf.h" |
45 | #include "gromacs/utility/smalloc.h" |
46 | #include "pbc.h" |
47 | #include "ns.h" |
48 | #include "macros.h" |
49 | #include "names.h" |
50 | #include "mshift.h" |
51 | #include "disre.h" |
52 | #include "orires.h" |
53 | #include "force.h" |
54 | #include "nonbonded.h" |
55 | |
56 | |
57 | |
58 | |
59 | void compute_factors_restangles(int type, const t_iparams forceparams[], |
60 | rvec delta_ante, rvec delta_post, |
61 | real *prefactor, real *ratio_ante, real *ratio_post, real *v) |
62 | { |
63 | real theta_equil, k_bending; |
64 | real cosine_theta_equil; |
65 | real c_ante, c_cros, c_post; |
66 | real norm; |
67 | real delta_cosine, cosine_theta; |
68 | real sine_theta_sq; |
69 | real term_theta_theta_equil; |
70 | |
71 | k_bending = forceparams[type].harmonic.krA; |
72 | theta_equil = forceparams[type].harmonic.rA*DEG2RAD(3.14159265358979323846/180.0); |
73 | theta_equil = M_PI3.14159265358979323846 - theta_equil; |
74 | cosine_theta_equil = cos(theta_equil); |
75 | |
76 | c_ante = iprod(delta_ante, delta_ante); |
77 | c_cros = iprod(delta_ante, delta_post); |
78 | c_post = iprod(delta_post, delta_post); |
79 | |
80 | norm = gmx_invsqrt(c_ante * c_post)gmx_software_invsqrt(c_ante * c_post); |
81 | cosine_theta = c_cros * norm; |
82 | sine_theta_sq = 1 - cosine_theta * cosine_theta; |
83 | |
84 | *ratio_ante = c_cros / c_ante; |
85 | *ratio_post = c_cros / c_post; |
86 | |
87 | delta_cosine = cosine_theta - cosine_theta_equil; |
88 | term_theta_theta_equil = 1 - cosine_theta * cosine_theta_equil; |
89 | *prefactor = -(k_bending) * delta_cosine * norm * term_theta_theta_equil / (sine_theta_sq * sine_theta_sq); |
90 | |
91 | *v = k_bending * 0.5 * delta_cosine * delta_cosine / sine_theta_sq; |
92 | |
93 | } |
94 | |
95 | |
96 | |
97 | |
98 | void compute_factors_restrdihs(int type, const t_iparams forceparams[], |
99 | rvec delta_ante, rvec delta_crnt, rvec delta_post, |
100 | real *factor_phi_ai_ante, real *factor_phi_ai_crnt, real *factor_phi_ai_post, |
101 | real *factor_phi_aj_ante, real *factor_phi_aj_crnt, real *factor_phi_aj_post, |
102 | real *factor_phi_ak_ante, real *factor_phi_ak_crnt, real *factor_phi_ak_post, |
103 | real *factor_phi_al_ante, real *factor_phi_al_crnt, real *factor_phi_al_post, |
104 | real *prefactor_phi, real *v) |
105 | { |
106 | |
107 | real phi0, sine_phi0, cosine_phi0; |
108 | real k_torsion; |
109 | real c_self_ante, c_self_crnt, c_self_post; |
110 | real c_cros_ante, c_cros_acrs, c_cros_post; |
111 | real c_prod, d_post, d_ante; |
112 | real sine_phi_sq, cosine_phi; |
113 | real delta_cosine, term_phi_phi0; |
114 | real ratio_phi_ante, ratio_phi_post; |
115 | real cos_phi, norm_phi; |
116 | |
117 | |
118 | phi0 = forceparams[type].pdihs.phiA * DEG2RAD(3.14159265358979323846/180.0); |
119 | cosine_phi0 = cos(phi0); |
120 | sine_phi0 = sin(phi0); |
| Value stored to 'sine_phi0' is never read |
121 | k_torsion = forceparams[type].pdihs.cpA; |
122 | |
123 | |
124 | |
125 | |
126 | |
127 | c_self_ante = iprod(delta_ante, delta_ante); |
128 | c_self_crnt = iprod(delta_crnt, delta_crnt); |
129 | c_self_post = iprod(delta_post, delta_post); |
130 | c_cros_ante = iprod(delta_ante, delta_crnt); |
131 | c_cros_acrs = iprod(delta_ante, delta_post); |
132 | c_cros_post = iprod(delta_crnt, delta_post); |
133 | c_prod = c_cros_ante * c_cros_post - c_self_crnt * c_cros_acrs; |
134 | d_ante = c_self_ante * c_self_crnt - c_cros_ante * c_cros_ante; |
135 | d_post = c_self_post * c_self_crnt - c_cros_post * c_cros_post; |
136 | |
137 | |
138 | |
139 | if (d_ante < GMX_REAL_EPS5.96046448E-08) |
140 | { |
141 | d_ante = GMX_REAL_EPS5.96046448E-08; |
142 | } |
143 | if (d_post < GMX_REAL_EPS5.96046448E-08) |
144 | { |
145 | d_post = GMX_REAL_EPS5.96046448E-08; |
146 | } |
147 | |
148 | |
149 | norm_phi = gmx_invsqrt(d_ante * d_post)gmx_software_invsqrt(d_ante * d_post); |
150 | cosine_phi = c_prod * norm_phi; |
151 | sine_phi_sq = 1.0 - cosine_phi * cosine_phi; |
152 | |
153 | |
154 | if (sine_phi_sq < 0.0) |
155 | { |
156 | sine_phi_sq = 0.0; |
157 | } |
158 | |
159 | |
160 | |
161 | |
162 | delta_cosine = cosine_phi - cosine_phi0; |
163 | term_phi_phi0 = 1 - cosine_phi * cosine_phi0; |
164 | |
165 | |
166 | |
167 | ratio_phi_ante = c_prod / d_ante; |
168 | ratio_phi_post = c_prod / d_post; |
169 | |
170 | |
171 | *prefactor_phi = -(k_torsion) * delta_cosine * norm_phi * term_phi_phi0 / (sine_phi_sq * sine_phi_sq); |
172 | |
173 | |
174 | |
175 | |
176 | |
177 | |
178 | |
179 | *factor_phi_ai_ante = ratio_phi_ante * c_self_crnt; |
180 | *factor_phi_ai_crnt = -c_cros_post - ratio_phi_ante * c_cros_ante; |
181 | *factor_phi_ai_post = c_self_crnt; |
182 | *factor_phi_aj_ante = -c_cros_post - ratio_phi_ante * (c_self_crnt + c_cros_ante); |
183 | *factor_phi_aj_crnt = c_cros_post + c_cros_acrs * 2.0 + ratio_phi_ante * (c_self_ante + c_cros_ante) + ratio_phi_post * c_self_post; |
184 | *factor_phi_aj_post = -(c_cros_ante + c_self_crnt) - ratio_phi_post * c_cros_post; |
185 | *factor_phi_ak_ante = c_cros_post + c_self_crnt + ratio_phi_ante * c_cros_ante; |
186 | *factor_phi_ak_crnt = -(c_cros_ante + c_cros_acrs * 2.0)- ratio_phi_ante * c_self_ante - ratio_phi_post * (c_self_post + c_cros_post); |
187 | *factor_phi_ak_post = c_cros_ante + ratio_phi_post * (c_self_crnt + c_cros_post); |
188 | *factor_phi_al_ante = -c_self_crnt; |
189 | *factor_phi_al_crnt = c_cros_ante + ratio_phi_post * c_cros_post; |
190 | *factor_phi_al_post = -ratio_phi_post * c_self_crnt; |
191 | |
192 | |
193 | *v = k_torsion * 0.5 * delta_cosine * delta_cosine / sine_phi_sq; |
194 | } |
195 | |
196 | |
197 | |
198 | |
199 | |
200 | |
201 | void compute_factors_cbtdihs(int type, const t_iparams forceparams[], |
202 | rvec delta_ante, rvec delta_crnt, rvec delta_post, |
203 | rvec f_phi_ai, rvec f_phi_aj, rvec f_phi_ak, rvec f_phi_al, |
204 | rvec f_theta_ante_ai, rvec f_theta_ante_aj, rvec f_theta_ante_ak, |
205 | rvec f_theta_post_aj, rvec f_theta_post_ak, rvec f_theta_post_al, |
206 | real * v) |
207 | { |
208 | int j, d; |
209 | real torsion_coef[NR_CBTDIHS6]; |
210 | real c_self_ante, c_self_crnt, c_self_post; |
211 | real c_cros_ante, c_cros_acrs, c_cros_post; |
212 | real c_prod, d_ante, d_post; |
213 | real norm_phi, norm_theta_ante, norm_theta_post; |
214 | real cosine_phi, cosine_theta_ante, cosine_theta_post; |
215 | real sine_theta_ante_sq, sine_theta_post_sq; |
216 | real sine_theta_ante, sine_theta_post; |
217 | real prefactor_phi; |
218 | real ratio_phi_ante, ratio_phi_post; |
219 | real r1, r2; |
220 | real factor_phi_ai_ante, factor_phi_ai_crnt, factor_phi_ai_post; |
221 | real factor_phi_aj_ante, factor_phi_aj_crnt, factor_phi_aj_post; |
222 | real factor_phi_ak_ante, factor_phi_ak_crnt, factor_phi_ak_post; |
223 | real factor_phi_al_ante, factor_phi_al_crnt, factor_phi_al_post; |
224 | real prefactor_theta_ante, ratio_theta_ante_ante, ratio_theta_ante_crnt; |
225 | real prefactor_theta_post, ratio_theta_post_crnt, ratio_theta_post_post; |
226 | |
227 | |
228 | |
229 | |
230 | |
231 | |
232 | |
233 | |
234 | |
235 | |
236 | |
237 | |
238 | |
239 | |
240 | |
241 | |
242 | |
243 | for (j = 0; (j < NR_CBTDIHS6); j++) |
244 | { |
245 | torsion_coef[j] = forceparams[type].cbtdihs.cbtcA[j]; |
246 | } |
247 | |
248 | |
249 | |
250 | |
251 | |
252 | |
253 | c_self_ante = iprod(delta_ante, delta_ante); |
254 | c_self_crnt = iprod(delta_crnt, delta_crnt); |
255 | c_self_post = iprod(delta_post, delta_post); |
256 | c_cros_ante = iprod(delta_ante, delta_crnt); |
257 | c_cros_acrs = iprod(delta_ante, delta_post); |
258 | c_cros_post = iprod(delta_crnt, delta_post); |
259 | c_prod = c_cros_ante * c_cros_post - c_self_crnt * c_cros_acrs; |
260 | d_ante = c_self_ante * c_self_crnt - c_cros_ante * c_cros_ante; |
261 | d_post = c_self_post * c_self_crnt - c_cros_post * c_cros_post; |
262 | |
263 | |
264 | |
265 | if (d_ante < GMX_REAL_EPS5.96046448E-08) |
266 | { |
267 | d_ante = GMX_REAL_EPS5.96046448E-08; |
268 | } |
269 | if (d_post < GMX_REAL_EPS5.96046448E-08) |
270 | { |
271 | d_post = GMX_REAL_EPS5.96046448E-08; |
272 | } |
273 | |
274 | |
275 | norm_phi = gmx_invsqrt(d_ante * d_post)gmx_software_invsqrt(d_ante * d_post); |
276 | norm_theta_ante = gmx_invsqrt(c_self_ante * c_self_crnt)gmx_software_invsqrt(c_self_ante * c_self_crnt); |
277 | norm_theta_post = gmx_invsqrt(c_self_crnt * c_self_post)gmx_software_invsqrt(c_self_crnt * c_self_post); |
278 | cosine_phi = c_prod * norm_phi; |
279 | cosine_theta_ante = c_cros_ante * norm_theta_ante; |
280 | cosine_theta_post = c_cros_post * norm_theta_post; |
281 | sine_theta_ante_sq = 1 - cosine_theta_ante * cosine_theta_ante; |
282 | sine_theta_post_sq = 1 - cosine_theta_post * cosine_theta_post; |
283 | |
284 | |
285 | if (sine_theta_ante_sq < 0.0) |
286 | { |
287 | sine_theta_ante_sq = 0.0; |
288 | } |
289 | if (sine_theta_post_sq < 0.0) |
290 | { |
291 | sine_theta_post_sq = 0.0; |
292 | } |
293 | |
294 | sine_theta_ante = sqrt(sine_theta_ante_sq); |
295 | sine_theta_post = sqrt(sine_theta_post_sq); |
296 | |
297 | |
298 | |
299 | |
300 | ratio_phi_ante = c_prod / d_ante; |
301 | ratio_phi_post = c_prod / d_post; |
302 | |
303 | |
304 | |
305 | r1 = cosine_phi; |
306 | |
307 | prefactor_phi = -torsion_coef[0] * norm_phi * (torsion_coef[2] + torsion_coef[3] * 2.0 * cosine_phi + torsion_coef[4] * 3.0 * (r1 * r1) + 4*torsion_coef[5]*r1*r1*r1) * |
308 | sine_theta_ante_sq * sine_theta_ante * sine_theta_post_sq * sine_theta_post; |
309 | |
310 | |
311 | |
312 | |
313 | |
314 | |
315 | |
316 | factor_phi_ai_ante = ratio_phi_ante * c_self_crnt; |
317 | factor_phi_ai_crnt = -c_cros_post - ratio_phi_ante * c_cros_ante; |
318 | factor_phi_ai_post = c_self_crnt; |
319 | factor_phi_aj_ante = -c_cros_post - ratio_phi_ante * (c_self_crnt + c_cros_ante); |
320 | factor_phi_aj_crnt = c_cros_post + c_cros_acrs * 2.0 + ratio_phi_ante * (c_self_ante + c_cros_ante) + ratio_phi_post * c_self_post; |
321 | factor_phi_aj_post = -(c_cros_ante + c_self_crnt) - ratio_phi_post * c_cros_post; |
322 | factor_phi_ak_ante = c_cros_post + c_self_crnt + ratio_phi_ante * c_cros_ante; |
323 | factor_phi_ak_crnt = -(c_cros_ante + c_cros_acrs * 2.0) - ratio_phi_ante * c_self_ante - ratio_phi_post * (c_self_post + c_cros_post); |
324 | factor_phi_ak_post = c_cros_ante + ratio_phi_post * (c_self_crnt + c_cros_post); |
325 | factor_phi_al_ante = -c_self_crnt; |
326 | factor_phi_al_crnt = c_cros_ante + ratio_phi_post * c_cros_post; |
327 | factor_phi_al_post = -ratio_phi_post * c_self_crnt; |
328 | |
329 | |
330 | for (d = 0; d < DIM3; d++) |
331 | { |
332 | f_phi_ai[d] = prefactor_phi * (factor_phi_ai_ante * delta_ante[d] + factor_phi_ai_crnt * delta_crnt[d] + factor_phi_ai_post * delta_post[d]); |
333 | f_phi_aj[d] = prefactor_phi * (factor_phi_aj_ante * delta_ante[d] + factor_phi_aj_crnt * delta_crnt[d] + factor_phi_aj_post * delta_post[d]); |
334 | f_phi_ak[d] = prefactor_phi * (factor_phi_ak_ante * delta_ante[d] + factor_phi_ak_crnt * delta_crnt[d] + factor_phi_ak_post * delta_post[d]); |
335 | f_phi_al[d] = prefactor_phi * (factor_phi_al_ante * delta_ante[d] + factor_phi_al_crnt * delta_crnt[d] + factor_phi_al_post * delta_post[d]); |
336 | } |
337 | |
338 | |
339 | |
340 | ratio_theta_ante_ante = c_cros_ante / c_self_ante; |
341 | ratio_theta_ante_crnt = c_cros_ante / c_self_crnt; |
342 | |
343 | |
344 | |
345 | r1 = cosine_phi; |
346 | |
347 | r2 = cosine_phi; |
348 | |
349 | prefactor_theta_ante = -torsion_coef[0] * norm_theta_ante * ( torsion_coef[1] + torsion_coef[2] * cosine_phi + torsion_coef[3] * (r1 * r1) + |
350 | torsion_coef[4] * (r2 * (r2 * r2))+ torsion_coef[5] * (r2 * (r2 * (r2 * r2)))) * (-3.0) * cosine_theta_ante * sine_theta_ante * sine_theta_post_sq * sine_theta_post; |
351 | |
352 | |
353 | |
354 | for (d = 0; d < DIM3; d++) |
355 | { |
356 | f_theta_ante_ai[d] = prefactor_theta_ante * (ratio_theta_ante_ante * delta_ante[d] - delta_crnt[d]); |
357 | f_theta_ante_aj[d] = prefactor_theta_ante * ((ratio_theta_ante_crnt + 1.0) * delta_crnt[d] - (ratio_theta_ante_ante + 1.0) * delta_ante[d]); |
358 | f_theta_ante_ak[d] = prefactor_theta_ante * (delta_ante[d] - ratio_theta_ante_crnt * delta_crnt[d]); |
359 | } |
360 | |
361 | |
362 | |
363 | |
364 | ratio_theta_post_crnt = c_cros_post / c_self_crnt; |
365 | ratio_theta_post_post = c_cros_post / c_self_post; |
366 | |
367 | |
368 | |
369 | r1 = cosine_phi; |
370 | |
371 | r2 = cosine_phi; |
372 | |
373 | prefactor_theta_post = -torsion_coef[0] * norm_theta_post * (torsion_coef[1] + torsion_coef[2] * cosine_phi + torsion_coef[3] * (r1 * r1) + |
374 | torsion_coef[4] * (r2 * (r2 * r2)) + torsion_coef[5] * (r2 * (r2 * (r2 * r2)))) * sine_theta_ante_sq * sine_theta_ante * (-3.0) * cosine_theta_post * sine_theta_post; |
375 | |
376 | |
377 | |
378 | for (d = 0; d < DIM3; d++) |
379 | { |
380 | f_theta_post_aj[d] = prefactor_theta_post * (ratio_theta_post_crnt * delta_crnt[d] - delta_post[d]); |
381 | f_theta_post_ak[d] = prefactor_theta_post * ((ratio_theta_post_post + 1.0) * delta_post[d] - (ratio_theta_post_crnt + 1.0) * delta_crnt[d]); |
382 | f_theta_post_al[d] = prefactor_theta_post * (delta_crnt[d] - ratio_theta_post_post * delta_post[d]); |
383 | } |
384 | r1 = cosine_phi; |
385 | r2 = cosine_phi; |
386 | |
387 | |
388 | *v = torsion_coef[0] * (torsion_coef[1] + torsion_coef[2] * cosine_phi + torsion_coef[3] * (r1 * r1) + |
389 | torsion_coef[4] * (r2 * (r2 * r2)) + torsion_coef[5] * (r2 * (r2 * (r2 * r2)))) * sine_theta_ante_sq * |
390 | sine_theta_ante * sine_theta_post_sq * sine_theta_post; |
391 | |
392 | |
393 | } |