2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2014, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
40 #include "gromacs/math/utilities.h"
43 #include "gromacs/utility/smalloc.h"
53 #include "nonbonded.h"
55 /* This function computes factors needed for restricted angle potential.
56 * For explanations on formula used see file "restcbt.h" */
58 void compute_factors_restangles(int type
, const t_iparams forceparams
[],
59 rvec delta_ante
, rvec delta_post
,
60 real
*prefactor
, real
*ratio_ante
, real
*ratio_post
, real
*v
)
62 real theta_equil
, k_bending
;
63 real cosine_theta_equil
;
64 real c_ante
, c_cros
, c_post
;
66 real delta_cosine
, cosine_theta
;
68 real term_theta_theta_equil
;
70 k_bending
= forceparams
[type
].harmonic
.krA
;
71 theta_equil
= forceparams
[type
].harmonic
.rA
*DEG2RAD
;
72 theta_equil
= M_PI
- theta_equil
;
73 cosine_theta_equil
= cos(theta_equil
);
75 c_ante
= iprod(delta_ante
, delta_ante
);
76 c_cros
= iprod(delta_ante
, delta_post
);
77 c_post
= iprod(delta_post
, delta_post
);
79 norm
= gmx_invsqrt(c_ante
* c_post
);
80 cosine_theta
= c_cros
* norm
;
81 sine_theta_sq
= 1 - cosine_theta
* cosine_theta
;
83 *ratio_ante
= c_cros
/ c_ante
;
84 *ratio_post
= c_cros
/ c_post
;
86 delta_cosine
= cosine_theta
- cosine_theta_equil
;
87 term_theta_theta_equil
= 1 - cosine_theta
* cosine_theta_equil
;
88 *prefactor
= -(k_bending
) * delta_cosine
* norm
* term_theta_theta_equil
/ (sine_theta_sq
* sine_theta_sq
);
90 *v
= k_bending
* 0.5 * delta_cosine
* delta_cosine
/ sine_theta_sq
;
95 /* Compute factors for restricted dihedral potential
96 * For explanations on formula used see file "restcbt.h" */
97 void compute_factors_restrdihs(int type
, const t_iparams forceparams
[],
98 rvec delta_ante
, rvec delta_crnt
, rvec delta_post
,
99 real
*factor_phi_ai_ante
, real
*factor_phi_ai_crnt
, real
*factor_phi_ai_post
,
100 real
*factor_phi_aj_ante
, real
*factor_phi_aj_crnt
, real
*factor_phi_aj_post
,
101 real
*factor_phi_ak_ante
, real
*factor_phi_ak_crnt
, real
*factor_phi_ak_post
,
102 real
*factor_phi_al_ante
, real
*factor_phi_al_crnt
, real
*factor_phi_al_post
,
103 real
*prefactor_phi
, real
*v
)
106 real phi0
, sine_phi0
, cosine_phi0
;
108 real c_self_ante
, c_self_crnt
, c_self_post
;
109 real c_cros_ante
, c_cros_acrs
, c_cros_post
;
110 real c_prod
, d_post
, d_ante
;
111 real sine_phi_sq
, cosine_phi
;
112 real delta_cosine
, term_phi_phi0
;
113 real ratio_phi_ante
, ratio_phi_post
;
114 real cos_phi
, norm_phi
;
116 /* Read parameters phi0 and k_torsion */
117 phi0
= forceparams
[type
].pdihs
.phiA
* DEG2RAD
;
118 cosine_phi0
= cos(phi0
);
119 sine_phi0
= sin(phi0
);
120 k_torsion
= forceparams
[type
].pdihs
.cpA
;
122 /* Computation of the cosine of the dihedral angle. The scalar ("dot") product method
123 * is used. c_*_* cummulate the scalar products of the differences of particles
124 * positions while c_prod, d_ante and d_post are differences of products of scalar
125 * terms that are parts of the derivatives of forces */
126 c_self_ante
= iprod(delta_ante
, delta_ante
);
127 c_self_crnt
= iprod(delta_crnt
, delta_crnt
);
128 c_self_post
= iprod(delta_post
, delta_post
);
129 c_cros_ante
= iprod(delta_ante
, delta_crnt
);
130 c_cros_acrs
= iprod(delta_ante
, delta_post
);
131 c_cros_post
= iprod(delta_crnt
, delta_post
);
132 c_prod
= c_cros_ante
* c_cros_post
- c_self_crnt
* c_cros_acrs
;
133 d_ante
= c_self_ante
* c_self_crnt
- c_cros_ante
* c_cros_ante
;
134 d_post
= c_self_post
* c_self_crnt
- c_cros_post
* c_cros_post
;
136 /* When three consecutive beads align, we obtain values close to zero.
137 * Here we avoid small values to prevent round-off errors. */
138 if (d_ante
< GMX_REAL_EPS
)
140 d_ante
= GMX_REAL_EPS
;
142 if (d_post
< GMX_REAL_EPS
)
144 d_post
= GMX_REAL_EPS
;
147 /* Computes the square of the sinus of phi in sine_phi_sq */
148 norm_phi
= gmx_invsqrt(d_ante
* d_post
);
149 cosine_phi
= c_prod
* norm_phi
;
150 sine_phi_sq
= 1.0 - cosine_phi
* cosine_phi
;
152 /* It is possible that cosine_phi is slightly bigger than 1.0 due to round-off errors. */
153 if (sine_phi_sq
< 0.0)
158 /* Computation of the differences of cosines (delta_cosine) and a term (term_phi_phi0)
159 * that is part of the common prefactor_phi */
161 delta_cosine
= cosine_phi
- cosine_phi0
;
162 term_phi_phi0
= 1 - cosine_phi
* cosine_phi0
;
165 /* Computation of ratios */
166 ratio_phi_ante
= c_prod
/ d_ante
;
167 ratio_phi_post
= c_prod
/ d_post
;
169 /* Computation of the prefactor - common term for all forces */
170 *prefactor_phi
= -(k_torsion
) * delta_cosine
* norm_phi
* term_phi_phi0
/ (sine_phi_sq
* sine_phi_sq
);
172 /* Computation of force factors. Factors factor_phi_* are coming from the
173 * derivatives of the torsion angle (phi) with respect to the beads ai, aj, al, ak,
174 * (four) coordinates and they are multiplied in the force computations with the
175 * differences of the particles positions stored in parameters delta_ante,
176 * delta_crnt, delta_post. For formulas see file "restcbt.h" */
178 *factor_phi_ai_ante
= ratio_phi_ante
* c_self_crnt
;
179 *factor_phi_ai_crnt
= -c_cros_post
- ratio_phi_ante
* c_cros_ante
;
180 *factor_phi_ai_post
= c_self_crnt
;
181 *factor_phi_aj_ante
= -c_cros_post
- ratio_phi_ante
* (c_self_crnt
+ c_cros_ante
);
182 *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
;
183 *factor_phi_aj_post
= -(c_cros_ante
+ c_self_crnt
) - ratio_phi_post
* c_cros_post
;
184 *factor_phi_ak_ante
= c_cros_post
+ c_self_crnt
+ ratio_phi_ante
* c_cros_ante
;
185 *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
);
186 *factor_phi_ak_post
= c_cros_ante
+ ratio_phi_post
* (c_self_crnt
+ c_cros_post
);
187 *factor_phi_al_ante
= -c_self_crnt
;
188 *factor_phi_al_crnt
= c_cros_ante
+ ratio_phi_post
* c_cros_post
;
189 *factor_phi_al_post
= -ratio_phi_post
* c_self_crnt
;
191 /* Contribution to energy - see formula in file "restcbt.h"*/
192 *v
= k_torsion
* 0.5 * delta_cosine
* delta_cosine
/ sine_phi_sq
;
197 /* Compute factors for CBT potential
198 * For explanations on formula used see file "restcbt.h" */
200 void compute_factors_cbtdihs(int type
, const t_iparams forceparams
[],
201 rvec delta_ante
, rvec delta_crnt
, rvec delta_post
,
202 rvec f_phi_ai
, rvec f_phi_aj
, rvec f_phi_ak
, rvec f_phi_al
,
203 rvec f_theta_ante_ai
, rvec f_theta_ante_aj
, rvec f_theta_ante_ak
,
204 rvec f_theta_post_aj
, rvec f_theta_post_ak
, rvec f_theta_post_al
,
208 real torsion_coef
[NR_CBTDIHS
];
209 real c_self_ante
, c_self_crnt
, c_self_post
;
210 real c_cros_ante
, c_cros_acrs
, c_cros_post
;
211 real c_prod
, d_ante
, d_post
;
212 real norm_phi
, norm_theta_ante
, norm_theta_post
;
213 real cosine_phi
, cosine_theta_ante
, cosine_theta_post
;
214 real sine_theta_ante_sq
, sine_theta_post_sq
;
215 real sine_theta_ante
, sine_theta_post
;
217 real ratio_phi_ante
, ratio_phi_post
;
219 real factor_phi_ai_ante
, factor_phi_ai_crnt
, factor_phi_ai_post
;
220 real factor_phi_aj_ante
, factor_phi_aj_crnt
, factor_phi_aj_post
;
221 real factor_phi_ak_ante
, factor_phi_ak_crnt
, factor_phi_ak_post
;
222 real factor_phi_al_ante
, factor_phi_al_crnt
, factor_phi_al_post
;
223 real prefactor_theta_ante
, ratio_theta_ante_ante
, ratio_theta_ante_crnt
;
224 real prefactor_theta_post
, ratio_theta_post_crnt
, ratio_theta_post_post
;
226 /* The formula for combined bending-torsion potential (see file "restcbt.h") contains
227 * in its expression not only the dihedral angle \f[\phi\f] but also \f[\theta_{i-1}\f]
228 * (theta_ante bellow) and \f[\theta_{i}\f] (theta_post bellow)--- the adjacent bending
229 * angles. The forces for the particles ai, aj, ak, al have components coming from the
230 * derivatives of the potential with respect to all three angles.
231 * This function is organised in 4 parts
232 * PART 1 - Computes force factors common to all the derivatives for the four particles
233 * PART 2 - Computes the force components due to the derivatives of dihedral angle Phi
234 * PART 3 - Computes the force components due to the derivatives of bending angle Theta_Ante
235 * PART 4 - Computes the force components due to the derivatives of bending angle Theta_Post
236 * Bellow we will respct thuis structure */
239 /* PART 1 - COMPUTES FORCE FACTORS COMMON TO ALL DERIVATIVES FOR THE FOUR PARTICLES */
242 for (j
= 0; (j
< NR_CBTDIHS
); j
++)
244 torsion_coef
[j
] = forceparams
[type
].cbtdihs
.cbtcA
[j
];
247 /* Computation of the cosine of the dihedral angle. The scalar ("dot") product method
248 * is used. c_*_* cummulate the scalar products of the differences of particles
249 * positions while c_prod, d_ante and d_post are differences of products of scalar
250 * terms that are parts of the derivatives of forces */
252 c_self_ante
= iprod(delta_ante
, delta_ante
);
253 c_self_crnt
= iprod(delta_crnt
, delta_crnt
);
254 c_self_post
= iprod(delta_post
, delta_post
);
255 c_cros_ante
= iprod(delta_ante
, delta_crnt
);
256 c_cros_acrs
= iprod(delta_ante
, delta_post
);
257 c_cros_post
= iprod(delta_crnt
, delta_post
);
258 c_prod
= c_cros_ante
* c_cros_post
- c_self_crnt
* c_cros_acrs
;
259 d_ante
= c_self_ante
* c_self_crnt
- c_cros_ante
* c_cros_ante
;
260 d_post
= c_self_post
* c_self_crnt
- c_cros_post
* c_cros_post
;
262 /* When three consecutive beads align, we obtain values close to zero.
263 Here we avoid small values to prevent round-off errors. */
264 if (d_ante
< GMX_REAL_EPS
)
266 d_ante
= GMX_REAL_EPS
;
268 if (d_post
< GMX_REAL_EPS
)
270 d_post
= GMX_REAL_EPS
;
273 /* Computations of cosines */
274 norm_phi
= gmx_invsqrt(d_ante
* d_post
);
275 norm_theta_ante
= gmx_invsqrt(c_self_ante
* c_self_crnt
);
276 norm_theta_post
= gmx_invsqrt(c_self_crnt
* c_self_post
);
277 cosine_phi
= c_prod
* norm_phi
;
278 cosine_theta_ante
= c_cros_ante
* norm_theta_ante
;
279 cosine_theta_post
= c_cros_post
* norm_theta_post
;
280 sine_theta_ante_sq
= 1 - cosine_theta_ante
* cosine_theta_ante
;
281 sine_theta_post_sq
= 1 - cosine_theta_post
* cosine_theta_post
;
283 /* It is possible that cosine_theta is slightly bigger than 1.0 due to round-off errors. */
284 if (sine_theta_ante_sq
< 0.0)
286 sine_theta_ante_sq
= 0.0;
288 if (sine_theta_post_sq
< 0.0)
290 sine_theta_post_sq
= 0.0;
293 sine_theta_ante
= sqrt(sine_theta_ante_sq
);
294 sine_theta_post
= sqrt(sine_theta_post_sq
);
296 /* PART 2 - COMPUTES FORCE COMPONENTS DUE TO DERIVATIVES TO DIHEDRAL ANGLE PHI */
298 /* Computation of ratios */
299 ratio_phi_ante
= c_prod
/ d_ante
;
300 ratio_phi_post
= c_prod
/ d_post
;
302 /* Computation of the prefactor */
303 /* Computing 2nd power */
306 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
) *
307 sine_theta_ante_sq
* sine_theta_ante
* sine_theta_post_sq
* sine_theta_post
;
309 /* Computation of factors (important for gaining speed). Factors factor_phi_* are coming from the
310 * derivatives of the torsion angle (phi) with respect to the beads ai, aj, al, ak,
311 * (four) coordinates and they are multiplied in the force computations with the
312 * differences of the particles positions stored in parameters delta_ante,
313 * delta_crnt, delta_post. For formulas see file "restcbt.h" */
315 factor_phi_ai_ante
= ratio_phi_ante
* c_self_crnt
;
316 factor_phi_ai_crnt
= -c_cros_post
- ratio_phi_ante
* c_cros_ante
;
317 factor_phi_ai_post
= c_self_crnt
;
318 factor_phi_aj_ante
= -c_cros_post
- ratio_phi_ante
* (c_self_crnt
+ c_cros_ante
);
319 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
;
320 factor_phi_aj_post
= -(c_cros_ante
+ c_self_crnt
) - ratio_phi_post
* c_cros_post
;
321 factor_phi_ak_ante
= c_cros_post
+ c_self_crnt
+ ratio_phi_ante
* c_cros_ante
;
322 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
);
323 factor_phi_ak_post
= c_cros_ante
+ ratio_phi_post
* (c_self_crnt
+ c_cros_post
);
324 factor_phi_al_ante
= -c_self_crnt
;
325 factor_phi_al_crnt
= c_cros_ante
+ ratio_phi_post
* c_cros_post
;
326 factor_phi_al_post
= -ratio_phi_post
* c_self_crnt
;
328 /* Computation of forces due to the derivatives of dihedral angle phi*/
329 for (d
= 0; d
< DIM
; d
++)
331 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
]);
332 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
]);
333 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
]);
334 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
]);
337 /* PART 3 - COMPUTES THE FORCE COMPONENTS DUE TO THE DERIVATIVES OF BENDING ANGLE THETHA_ANTHE */
338 /* Computation of ratios */
339 ratio_theta_ante_ante
= c_cros_ante
/ c_self_ante
;
340 ratio_theta_ante_crnt
= c_cros_ante
/ c_self_crnt
;
342 /* Computation of the prefactor */
343 /* Computing 2nd power */
345 /* Computing 3rd power */
348 prefactor_theta_ante
= -torsion_coef
[0] * norm_theta_ante
* ( torsion_coef
[1] + torsion_coef
[2] * cosine_phi
+ torsion_coef
[3] * (r1
* r1
) +
349 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
;
352 /* Computation of forces due to the derivatives of bending angle theta_ante */
353 for (d
= 0; d
< DIM
; d
++)
355 f_theta_ante_ai
[d
] = prefactor_theta_ante
* (ratio_theta_ante_ante
* delta_ante
[d
] - delta_crnt
[d
]);
356 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
]);
357 f_theta_ante_ak
[d
] = prefactor_theta_ante
* (delta_ante
[d
] - ratio_theta_ante_crnt
* delta_crnt
[d
]);
360 /* PART 4 - COMPUTES THE FORCE COMPONENTS DUE TO THE DERIVATIVES OF THE BENDING ANGLE THETA_POST */
362 /* Computation of ratios */
363 ratio_theta_post_crnt
= c_cros_post
/ c_self_crnt
;
364 ratio_theta_post_post
= c_cros_post
/ c_self_post
;
366 /* Computation of the prefactor */
367 /* Computing 2nd power */
369 /* Computing 3rd power */
372 prefactor_theta_post
= -torsion_coef
[0] * norm_theta_post
* (torsion_coef
[1] + torsion_coef
[2] * cosine_phi
+ torsion_coef
[3] * (r1
* r1
) +
373 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
;
376 /* Computation of forces due to the derivatives of bending angle Theta_Post */
377 for (d
= 0; d
< DIM
; d
++)
379 f_theta_post_aj
[d
] = prefactor_theta_post
* (ratio_theta_post_crnt
* delta_crnt
[d
] - delta_post
[d
]);
380 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
]);
381 f_theta_post_al
[d
] = prefactor_theta_post
* (delta_crnt
[d
] - ratio_theta_post_post
* delta_post
[d
]);
386 /* Contribution to energy - for formula see file "restcbt.h" */
387 *v
= torsion_coef
[0] * (torsion_coef
[1] + torsion_coef
[2] * cosine_phi
+ torsion_coef
[3] * (r1
* r1
) +
388 torsion_coef
[4] * (r2
* (r2
* r2
)) + torsion_coef
[5] * (r2
* (r2
* (r2
* r2
)))) * sine_theta_ante_sq
*
389 sine_theta_ante
* sine_theta_post_sq
* sine_theta_post
;