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.
37 /*! \libinternal \file
41 * This file contains function declarations necessary
42 for computations of forces due to restricted angle, restricted dihedral and
43 combined bending-torsion potentials.
46 * \author Nicolae Goga
51 #ifndef GMX_BONDED_RESTCBT_H
52 #define GMX_BONDED_RESTCBT_H
54 #include "gromacs/legacyheaders/types/simple.h"
55 #include "gromacs/topology/idef.h"
56 #include "gromacs/math/vec.h"
63 /*! \brief This function computes factors needed for restricted angle potentials.
65 * The restricted angle potential is used in coarse-grained simulations to avoid singularities
66 * when three particles align and the dihedral angle and dihedral potential cannot be calculated.
67 * This potential is calculated using the formula:
68 * \f[V_{\rm ReB}(\theta_i) = \frac{1}{2} k_{\theta} \frac{(\cos\theta_i - \cos\theta_0)^2}{\sin^2\theta_i}\f]
69 * (see section "Restricted Bending Potential" from the manual).
70 * The derivative of the restricted angle potential is calculated as:
71 * \f[\frac{\partial V_{\rm ReB}(\theta_i)} {\partial \vec r_{k}} = \frac{dV_{\rm ReB}(\theta_i)}{dcos\theta_i} \frac{\partial cos\theta_i}{\partial \vec r_{k}}\f]
72 * where all the derivatives of the bending angle with respect to Cartesian coordinates are calculated as in Allen & Tildesley (pp. 330-332)
74 * \param[in] type type of force parameters
75 * \param[in] forceparams array of parameters for force computations
76 * \param[in] delta_ante distance between the first two particles
77 * \param[in] delta_post distance between the last two particles
78 * \param[out] prefactor common term that comes in front of each force
79 * \param[out] ratio_ante ratio of scalar products of delta_ante with delta_post
80 and delta_ante with delta_ante
81 * \param[out] ratio_post ratio of scalar products of delta_ante with delta_post
82 and delta_post with delta_ante
83 * \param[out] v contribution to energy (see formula above)
87 void compute_factors_restangles(int type, const t_iparams forceparams[],
88 rvec delta_ante, rvec delta_post,
89 real *prefactor, real *ratio_ante, real *ratio_post, real *v);
92 /*! \brief Compute factors for restricted dihedral potentials.
94 * The restricted dihedral potential is the equivalent of the restricted bending potential
95 * for the dihedral angle. It imposes the dihedral angle to have only one equilibrium value.
96 * This potential is calculated using the formula:
97 * \f[V_{\rm ReT}(\phi_i) = \frac{1}{2} k_{\phi} \frac{(\cos\phi_i - \cos\phi_0)^2}{\sin^2\phi_i}\f]
98 * (see section "Proper dihedrals: Restricted torsion potential" from the manual).
99 * The derivative of the restricted dihedral potential is calculated as:
100 * \f[\frac{\partial V_{\rm ReT}(\phi_i)} {\partial \vec r_{k}} = \frac{dV_{\rm ReT}(\phi_i)}{dcos\phi_i} \frac{\partial cos\phi_i}{\partial \vec r_{k}}\f]
101 * where all the derivatives of the dihedral angle with respect to Cartesian coordinates
102 * are calculated as in Allen & Tildesley (pp. 330-332). Factors factor_phi_* are coming from the
103 * derivatives of the torsion angle (phi) with respect to the beads ai, aj, ak, al, (four) coordinates
104 * and they are multiplied in the force computations with the particle distance
105 * stored in parameters delta_ante, delta_crnt, delta_post.
107 * \param[in] type type of force parameters
108 * \param[in] forceparams array of parameters for force computations
109 * \param[in] delta_ante distance between the first two particles
110 * \param[in] delta_crnt distance between the middle pair of particles
111 * \param[in] delta_post distance between the last two particles
112 * \param[out] factor_phi_ai_ante force factor for particle ai multiplied with delta_ante
113 * \param[out] factor_phi_ai_crnt force factor for particle ai multiplied with delta_crnt
114 * \param[out] factor_phi_ai_post force factor for particle ai multiplied with delta_post
115 * \param[out] factor_phi_aj_ante force factor for particle aj multiplied with delta_ante
116 * \param[out] factor_phi_aj_crnt force factor for particle aj multiplied with delta_crnt
117 * \param[out] factor_phi_aj_post force factor for particle aj multiplied with delta_post
118 * \param[out] factor_phi_ak_ante force factor for particle ak multiplied with delta_ante
119 * \param[out] factor_phi_ak_crnt force factor for particle ak multiplied with delta_crnt
120 * \param[out] factor_phi_ak_post force factor for particle ak multiplied with delta_post
121 * \param[out] factor_phi_al_ante force factor for particle al multiplied with delta_ante
122 * \param[out] factor_phi_al_crnt force factor for particle al multiplied with delta_crnt
123 * \param[out] factor_phi_al_post force factor for particle al multiplied with delta_post
124 * \param[out] prefactor_phi multiplication constant of the torsion force
125 * \param[out] v contribution to energy (see formula above)
128 void compute_factors_restrdihs(int type, const t_iparams forceparams[],
129 rvec delta_ante, rvec delta_crnt, rvec delta_post,
130 real *factor_phi_ai_ante, real *factor_phi_ai_crnt, real *factor_phi_ai_post,
131 real *factor_phi_aj_ante, real *factor_phi_aj_crnt, real *factor_phi_aj_post,
132 real *factor_phi_ak_ante, real *factor_phi_ak_crnt, real *factor_phi_ak_post,
133 real *factor_phi_al_ante, real *factor_phi_al_crnt, real *factor_phi_al_post,
134 real *prefactor_phi, real *v);
136 /*! \brief Compute factors for combined bending-torsion (CBT) potentials.
138 * The combined bending-torsion potential goes to zero in a very smooth manner, eliminating the numerical
139 * instabilities, when three coarse-grained particles align and the dihedral angle and standard
140 * dihedral potentials cannot be calculated. The CBT potential is calculated using the formula:
141 * \f[V_{\rm CBT}(\theta_{i-1}, \theta_i, \phi_i) = k_{\phi} \sin^3\theta_{i-1} \sin^3\theta_{i}
142 * \sum_{n=0}^4 { a_n \cos^n\phi_i}\f] (see section "Proper dihedrals: Combined bending-torsion potential" from the manual).
143 * It contains in its expression not only the dihedral angle \f$\phi\f$
144 * but also \f$\theta_{i-1}\f$ (denoted as theta_ante below) and \f$\theta_{i}\f$ (denoted as theta_post below)
145 * --- the adjacent bending angles. The derivative of the CBT potential is calculated as:
146 * \f[\frac{\partial V_{\rm CBT}(\theta_{i-1},\theta_i,\phi_i)} {\partial \vec r_{l}} = \frac{\partial V_
147 * {\rm CBT}}{\partial \theta_{i-1}} \frac{\partial \theta_{i-1}}{\partial \vec r_{l}} +
148 * \frac{\partial V_{\rm CBT}}{\partial \phi_{i }} \frac{\partial \phi_{i }}{\partial \vec r_{l}}\f]
149 * where all the derivatives of the angles with respect to Cartesian coordinates are calculated as
150 * in Allen & Tildesley (pp. 330-332). Factors f_phi_* come from the derivatives of the torsion angle
151 * with respect to the beads ai, aj, ak, al (four) coordinates; f_theta_ante_* come from the derivatives of
152 * the bending angle theta_ante (theta_{i-1} in formula above) with respect to the beads ai, aj, ak (three
153 * particles) coordinates and f_theta_post_* come from the derivatives of the bending angle theta_post
154 * (theta_{i} in formula above) with respect to the beads aj, ak, al (three particles) coordinates.
156 * \param[in] type type of force parameters
157 * \param[in] forceparams array of parameters for force computations
158 * \param[in] delta_ante distance between the first two particles
159 * \param[in] delta_crnt distance between the middle pair of particles
160 * \param[in] delta_post distance between the last two particles
161 * \param[out] f_phi_ai force for particle ai due to derivative in phi angle
162 * \param[out] f_phi_aj force for particle aj due to derivative in phi angle
163 * \param[out] f_phi_ak force for particle ak due to derivative in phi angle
164 * \param[out] f_phi_al force for particle al due to derivative in phi angle
165 * \param[out] f_theta_ante_ai force for particle ai due to derivative in theta_ante angle
166 * \param[out] f_theta_ante_aj force for particle aj due to derivative in theta_ante angle
167 * \param[out] f_theta_ante_ak force for particle ak due to derivative in theta_ante angle
168 * \param[out] f_theta_post_aj force for particle aj due to derivative in theta_post angle
169 * \param[out] f_theta_post_ak force for particle ak due to derivative in theta_post angle
170 * \param[out] f_theta_post_al force for particle al due to derivative in theta_psot angle
171 * \param[out] v contribution to energy (see formula above)
174 void compute_factors_cbtdihs(int type, const t_iparams forceparams[],
175 rvec delta_ante, rvec delta_crnt, rvec delta_post,
176 rvec f_phi_ai, rvec f_phi_aj, rvec f_phi_ak, rvec f_phi_al,
177 rvec f_theta_ante_ai, rvec f_theta_ante_aj, rvec f_theta_ante_ak,
178 rvec f_theta_post_aj, rvec f_theta_post_ak, rvec f_theta_post_al,