src/gromacs/listed_forces/restcbt.h

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