2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2014,2015,2017,2019, 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.
38 * This file contains function definitions necessary
39 * for computations of forces due to restricted angle, restricted dihedral and
40 * combined bending-torsion potentials.
42 * \author Nicolae Goga
44 * \ingroup module_listed_forces
52 #include "gromacs/math/functions.h"
53 #include "gromacs/math/units.h"
54 #include "gromacs/math/utilities.h"
55 #include "gromacs/math/vec.h"
56 #include "gromacs/topology/idef.h"
58 /* This function computes factors needed for restricted angle potential.
59 * For explanations on formula used see file "restcbt.h" */
61 void compute_factors_restangles(int type,
62 const t_iparams forceparams[],
70 // These variables are double to make the code
72 double theta_equil, k_bending;
73 double cosine_theta_equil;
74 double c_ante, c_cros, c_post;
76 double delta_cosine, cosine_theta;
78 double term_theta_theta_equil;
80 k_bending = forceparams[type].harmonic.krA;
81 theta_equil = forceparams[type].harmonic.rA * DEG2RAD;
82 theta_equil = M_PI - theta_equil;
83 cosine_theta_equil = cos(theta_equil);
85 c_ante = iprod(delta_ante, delta_ante);
86 c_cros = iprod(delta_ante, delta_post);
87 c_post = iprod(delta_post, delta_post);
89 norm = gmx::invsqrt(c_ante * c_post);
90 cosine_theta = c_cros * norm;
91 sine_theta_sq = 1 - cosine_theta * cosine_theta;
93 *ratio_ante = c_cros / c_ante;
94 *ratio_post = c_cros / c_post;
96 delta_cosine = cosine_theta - cosine_theta_equil;
97 term_theta_theta_equil = 1 - cosine_theta * cosine_theta_equil;
98 *prefactor = -(k_bending)*delta_cosine * norm * term_theta_theta_equil / (sine_theta_sq * sine_theta_sq);
100 *v = k_bending * 0.5 * delta_cosine * delta_cosine / sine_theta_sq;
104 /* Compute factors for restricted dihedral potential
105 * For explanations on formula used see file "restcbt.h" */
106 void compute_factors_restrdihs(int type,
107 const t_iparams forceparams[],
111 real* factor_phi_ai_ante,
112 real* factor_phi_ai_crnt,
113 real* factor_phi_ai_post,
114 real* factor_phi_aj_ante,
115 real* factor_phi_aj_crnt,
116 real* factor_phi_aj_post,
117 real* factor_phi_ak_ante,
118 real* factor_phi_ak_crnt,
119 real* factor_phi_ak_post,
120 real* factor_phi_al_ante,
121 real* factor_phi_al_crnt,
122 real* factor_phi_al_post,
127 real phi0, cosine_phi0;
129 real c_self_ante, c_self_crnt, c_self_post;
130 real c_cros_ante, c_cros_acrs, c_cros_post;
131 real c_prod, d_post, d_ante;
132 real sine_phi_sq, cosine_phi;
133 real delta_cosine, term_phi_phi0;
134 real ratio_phi_ante, ratio_phi_post;
137 /* Read parameters phi0 and k_torsion */
138 phi0 = forceparams[type].pdihs.phiA * DEG2RAD;
139 cosine_phi0 = cos(phi0);
140 k_torsion = forceparams[type].pdihs.cpA;
142 /* Computation of the cosine of the dihedral angle. The scalar ("dot") product method
143 * is used. c_*_* cummulate the scalar products of the differences of particles
144 * positions while c_prod, d_ante and d_post are differences of products of scalar
145 * terms that are parts of the derivatives of forces */
146 c_self_ante = iprod(delta_ante, delta_ante);
147 c_self_crnt = iprod(delta_crnt, delta_crnt);
148 c_self_post = iprod(delta_post, delta_post);
149 c_cros_ante = iprod(delta_ante, delta_crnt);
150 c_cros_acrs = iprod(delta_ante, delta_post);
151 c_cros_post = iprod(delta_crnt, delta_post);
152 c_prod = c_cros_ante * c_cros_post - c_self_crnt * c_cros_acrs;
153 d_ante = c_self_ante * c_self_crnt - c_cros_ante * c_cros_ante;
154 d_post = c_self_post * c_self_crnt - c_cros_post * c_cros_post;
156 /* When three consecutive beads align, we obtain values close to zero.
157 * Here we avoid small values to prevent round-off errors. */
158 if (d_ante < GMX_REAL_EPS)
160 d_ante = GMX_REAL_EPS;
162 if (d_post < GMX_REAL_EPS)
164 d_post = GMX_REAL_EPS;
167 /* Computes the square of the sinus of phi in sine_phi_sq */
168 norm_phi = gmx::invsqrt(d_ante * d_post);
169 cosine_phi = c_prod * norm_phi;
170 sine_phi_sq = 1.0 - cosine_phi * cosine_phi;
172 /* It is possible that cosine_phi is slightly bigger than 1.0 due to round-off errors. */
173 if (sine_phi_sq < 0.0)
178 /* Computation of the differences of cosines (delta_cosine) and a term (term_phi_phi0)
179 * that is part of the common prefactor_phi */
181 delta_cosine = cosine_phi - cosine_phi0;
182 term_phi_phi0 = 1 - cosine_phi * cosine_phi0;
185 /* Computation of ratios */
186 ratio_phi_ante = c_prod / d_ante;
187 ratio_phi_post = c_prod / d_post;
189 /* Computation of the prefactor - common term for all forces */
190 *prefactor_phi = -(k_torsion)*delta_cosine * norm_phi * term_phi_phi0 / (sine_phi_sq * sine_phi_sq);
192 /* Computation of force factors. Factors factor_phi_* are coming from the
193 * derivatives of the torsion angle (phi) with respect to the beads ai, aj, al, ak,
194 * (four) coordinates and they are multiplied in the force computations with the
195 * differences of the particles positions stored in parameters delta_ante,
196 * delta_crnt, delta_post. For formulas see file "restcbt.h" */
198 *factor_phi_ai_ante = ratio_phi_ante * c_self_crnt;
199 *factor_phi_ai_crnt = -c_cros_post - ratio_phi_ante * c_cros_ante;
200 *factor_phi_ai_post = c_self_crnt;
201 *factor_phi_aj_ante = -c_cros_post - ratio_phi_ante * (c_self_crnt + c_cros_ante);
202 *factor_phi_aj_crnt = c_cros_post + c_cros_acrs * 2.0
203 + ratio_phi_ante * (c_self_ante + c_cros_ante) + ratio_phi_post * c_self_post;
204 *factor_phi_aj_post = -(c_cros_ante + c_self_crnt) - ratio_phi_post * c_cros_post;
205 *factor_phi_ak_ante = c_cros_post + c_self_crnt + ratio_phi_ante * c_cros_ante;
206 *factor_phi_ak_crnt = -(c_cros_ante + c_cros_acrs * 2.0) - ratio_phi_ante * c_self_ante
207 - ratio_phi_post * (c_self_post + c_cros_post);
208 *factor_phi_ak_post = c_cros_ante + ratio_phi_post * (c_self_crnt + c_cros_post);
209 *factor_phi_al_ante = -c_self_crnt;
210 *factor_phi_al_crnt = c_cros_ante + ratio_phi_post * c_cros_post;
211 *factor_phi_al_post = -ratio_phi_post * c_self_crnt;
213 /* Contribution to energy - see formula in file "restcbt.h"*/
214 *v = k_torsion * 0.5 * delta_cosine * delta_cosine / sine_phi_sq;
218 /* Compute factors for CBT potential
219 * For explanations on formula used see file "restcbt.h" */
221 void compute_factors_cbtdihs(int type,
222 const t_iparams forceparams[],
230 rvec f_theta_ante_ai,
231 rvec f_theta_ante_aj,
232 rvec f_theta_ante_ak,
233 rvec f_theta_post_aj,
234 rvec f_theta_post_ak,
235 rvec f_theta_post_al,
239 real torsion_coef[NR_CBTDIHS];
240 real c_self_ante, c_self_crnt, c_self_post;
241 real c_cros_ante, c_cros_acrs, c_cros_post;
242 real c_prod, d_ante, d_post;
243 real norm_phi, norm_theta_ante, norm_theta_post;
244 real cosine_phi, cosine_theta_ante, cosine_theta_post;
245 real sine_theta_ante_sq, sine_theta_post_sq;
246 real sine_theta_ante, sine_theta_post;
248 real ratio_phi_ante, ratio_phi_post;
250 real factor_phi_ai_ante, factor_phi_ai_crnt, factor_phi_ai_post;
251 real factor_phi_aj_ante, factor_phi_aj_crnt, factor_phi_aj_post;
252 real factor_phi_ak_ante, factor_phi_ak_crnt, factor_phi_ak_post;
253 real factor_phi_al_ante, factor_phi_al_crnt, factor_phi_al_post;
254 real prefactor_theta_ante, ratio_theta_ante_ante, ratio_theta_ante_crnt;
255 real prefactor_theta_post, ratio_theta_post_crnt, ratio_theta_post_post;
257 /* The formula for combined bending-torsion potential (see file "restcbt.h") contains
258 * in its expression not only the dihedral angle \f[\phi\f] but also \f[\theta_{i-1}\f]
259 * (theta_ante bellow) and \f[\theta_{i}\f] (theta_post bellow)--- the adjacent bending
260 * angles. The forces for the particles ai, aj, ak, al have components coming from the
261 * derivatives of the potential with respect to all three angles.
262 * This function is organised in 4 parts
263 * PART 1 - Computes force factors common to all the derivatives for the four particles
264 * PART 2 - Computes the force components due to the derivatives of dihedral angle Phi
265 * PART 3 - Computes the force components due to the derivatives of bending angle Theta_Ante
266 * PART 4 - Computes the force components due to the derivatives of bending angle Theta_Post
267 * Bellow we will respct thuis structure */
270 /* PART 1 - COMPUTES FORCE FACTORS COMMON TO ALL DERIVATIVES FOR THE FOUR PARTICLES */
273 for (j = 0; (j < NR_CBTDIHS); j++)
275 torsion_coef[j] = forceparams[type].cbtdihs.cbtcA[j];
278 /* Computation of the cosine of the dihedral angle. The scalar ("dot") product method
279 * is used. c_*_* cummulate the scalar products of the differences of particles
280 * positions while c_prod, d_ante and d_post are differences of products of scalar
281 * terms that are parts of the derivatives of forces */
283 c_self_ante = iprod(delta_ante, delta_ante);
284 c_self_crnt = iprod(delta_crnt, delta_crnt);
285 c_self_post = iprod(delta_post, delta_post);
286 c_cros_ante = iprod(delta_ante, delta_crnt);
287 c_cros_acrs = iprod(delta_ante, delta_post);
288 c_cros_post = iprod(delta_crnt, delta_post);
289 c_prod = c_cros_ante * c_cros_post - c_self_crnt * c_cros_acrs;
290 d_ante = c_self_ante * c_self_crnt - c_cros_ante * c_cros_ante;
291 d_post = c_self_post * c_self_crnt - c_cros_post * c_cros_post;
293 /* When three consecutive beads align, we obtain values close to zero.
294 Here we avoid small values to prevent round-off errors. */
295 if (d_ante < GMX_REAL_EPS)
297 d_ante = GMX_REAL_EPS;
299 if (d_post < GMX_REAL_EPS)
301 d_post = GMX_REAL_EPS;
304 /* Computations of cosines */
305 norm_phi = gmx::invsqrt(d_ante * d_post);
306 norm_theta_ante = gmx::invsqrt(c_self_ante * c_self_crnt);
307 norm_theta_post = gmx::invsqrt(c_self_crnt * c_self_post);
308 cosine_phi = c_prod * norm_phi;
309 cosine_theta_ante = c_cros_ante * norm_theta_ante;
310 cosine_theta_post = c_cros_post * norm_theta_post;
311 sine_theta_ante_sq = 1 - cosine_theta_ante * cosine_theta_ante;
312 sine_theta_post_sq = 1 - cosine_theta_post * cosine_theta_post;
314 /* It is possible that cosine_theta is slightly bigger than 1.0 due to round-off errors. */
315 if (sine_theta_ante_sq < 0.0)
317 sine_theta_ante_sq = 0.0;
319 if (sine_theta_post_sq < 0.0)
321 sine_theta_post_sq = 0.0;
324 sine_theta_ante = std::sqrt(sine_theta_ante_sq);
325 sine_theta_post = std::sqrt(sine_theta_post_sq);
327 /* PART 2 - COMPUTES FORCE COMPONENTS DUE TO DERIVATIVES TO DIHEDRAL ANGLE PHI */
329 /* Computation of ratios */
330 ratio_phi_ante = c_prod / d_ante;
331 ratio_phi_post = c_prod / d_post;
333 /* Computation of the prefactor */
334 /* Computing 2nd power */
337 prefactor_phi = -torsion_coef[0] * norm_phi
338 * (torsion_coef[2] + torsion_coef[3] * 2.0 * cosine_phi
339 + torsion_coef[4] * 3.0 * (r1 * r1) + 4 * torsion_coef[5] * r1 * r1 * r1)
340 * sine_theta_ante_sq * sine_theta_ante * sine_theta_post_sq * sine_theta_post;
342 /* Computation of factors (important for gaining speed). Factors factor_phi_* are coming from the
343 * derivatives of the torsion angle (phi) with respect to the beads ai, aj, al, ak,
344 * (four) coordinates and they are multiplied in the force computations with the
345 * differences of the particles positions stored in parameters delta_ante,
346 * delta_crnt, delta_post. For formulas see file "restcbt.h" */
348 factor_phi_ai_ante = ratio_phi_ante * c_self_crnt;
349 factor_phi_ai_crnt = -c_cros_post - ratio_phi_ante * c_cros_ante;
350 factor_phi_ai_post = c_self_crnt;
351 factor_phi_aj_ante = -c_cros_post - ratio_phi_ante * (c_self_crnt + c_cros_ante);
352 factor_phi_aj_crnt = c_cros_post + c_cros_acrs * 2.0
353 + ratio_phi_ante * (c_self_ante + c_cros_ante) + ratio_phi_post * c_self_post;
354 factor_phi_aj_post = -(c_cros_ante + c_self_crnt) - ratio_phi_post * c_cros_post;
355 factor_phi_ak_ante = c_cros_post + c_self_crnt + ratio_phi_ante * c_cros_ante;
356 factor_phi_ak_crnt = -(c_cros_ante + c_cros_acrs * 2.0) - ratio_phi_ante * c_self_ante
357 - ratio_phi_post * (c_self_post + c_cros_post);
358 factor_phi_ak_post = c_cros_ante + ratio_phi_post * (c_self_crnt + c_cros_post);
359 factor_phi_al_ante = -c_self_crnt;
360 factor_phi_al_crnt = c_cros_ante + ratio_phi_post * c_cros_post;
361 factor_phi_al_post = -ratio_phi_post * c_self_crnt;
363 /* Computation of forces due to the derivatives of dihedral angle phi*/
364 for (d = 0; d < DIM; d++)
366 f_phi_ai[d] = prefactor_phi
367 * (factor_phi_ai_ante * delta_ante[d] + factor_phi_ai_crnt * delta_crnt[d]
368 + factor_phi_ai_post * delta_post[d]);
369 f_phi_aj[d] = prefactor_phi
370 * (factor_phi_aj_ante * delta_ante[d] + factor_phi_aj_crnt * delta_crnt[d]
371 + factor_phi_aj_post * delta_post[d]);
372 f_phi_ak[d] = prefactor_phi
373 * (factor_phi_ak_ante * delta_ante[d] + factor_phi_ak_crnt * delta_crnt[d]
374 + factor_phi_ak_post * delta_post[d]);
375 f_phi_al[d] = prefactor_phi
376 * (factor_phi_al_ante * delta_ante[d] + factor_phi_al_crnt * delta_crnt[d]
377 + factor_phi_al_post * delta_post[d]);
380 /* PART 3 - COMPUTES THE FORCE COMPONENTS DUE TO THE DERIVATIVES OF BENDING ANGLE THETHA_ANTHE */
381 /* Computation of ratios */
382 ratio_theta_ante_ante = c_cros_ante / c_self_ante;
383 ratio_theta_ante_crnt = c_cros_ante / c_self_crnt;
385 /* Computation of the prefactor */
386 /* Computing 2nd power */
388 /* Computing 3rd power */
391 prefactor_theta_ante =
392 -torsion_coef[0] * norm_theta_ante
393 * (torsion_coef[1] + torsion_coef[2] * cosine_phi + torsion_coef[3] * (r1 * r1)
394 + torsion_coef[4] * (r2 * (r2 * r2)) + torsion_coef[5] * (r2 * (r2 * (r2 * r2))))
395 * (-3.0) * cosine_theta_ante * sine_theta_ante * sine_theta_post_sq * sine_theta_post;
398 /* Computation of forces due to the derivatives of bending angle theta_ante */
399 for (d = 0; d < DIM; d++)
402 prefactor_theta_ante * (ratio_theta_ante_ante * delta_ante[d] - delta_crnt[d]);
403 f_theta_ante_aj[d] = prefactor_theta_ante
404 * ((ratio_theta_ante_crnt + 1.0) * delta_crnt[d]
405 - (ratio_theta_ante_ante + 1.0) * delta_ante[d]);
407 prefactor_theta_ante * (delta_ante[d] - ratio_theta_ante_crnt * delta_crnt[d]);
410 /* PART 4 - COMPUTES THE FORCE COMPONENTS DUE TO THE DERIVATIVES OF THE BENDING ANGLE THETA_POST */
412 /* Computation of ratios */
413 ratio_theta_post_crnt = c_cros_post / c_self_crnt;
414 ratio_theta_post_post = c_cros_post / c_self_post;
416 /* Computation of the prefactor */
417 /* Computing 2nd power */
419 /* Computing 3rd power */
422 prefactor_theta_post =
423 -torsion_coef[0] * norm_theta_post
424 * (torsion_coef[1] + torsion_coef[2] * cosine_phi + torsion_coef[3] * (r1 * r1)
425 + torsion_coef[4] * (r2 * (r2 * r2)) + torsion_coef[5] * (r2 * (r2 * (r2 * r2))))
426 * sine_theta_ante_sq * sine_theta_ante * (-3.0) * cosine_theta_post * sine_theta_post;
429 /* Computation of forces due to the derivatives of bending angle Theta_Post */
430 for (d = 0; d < DIM; d++)
433 prefactor_theta_post * (ratio_theta_post_crnt * delta_crnt[d] - delta_post[d]);
434 f_theta_post_ak[d] = prefactor_theta_post
435 * ((ratio_theta_post_post + 1.0) * delta_post[d]
436 - (ratio_theta_post_crnt + 1.0) * delta_crnt[d]);
438 prefactor_theta_post * (delta_crnt[d] - ratio_theta_post_post * delta_post[d]);
443 /* Contribution to energy - for formula see file "restcbt.h" */
445 * (torsion_coef[1] + torsion_coef[2] * cosine_phi + torsion_coef[3] * (r1 * r1)
446 + torsion_coef[4] * (r2 * (r2 * r2)) + torsion_coef[5] * (r2 * (r2 * (r2 * r2))))
447 * sine_theta_ante_sq * sine_theta_ante * sine_theta_post_sq * sine_theta_post;