src/gromacs/listed_forces/restcbt.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2014,2015,2017,2019,2021, 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 /*! \internal \file
  36  *
  37  * \brief
  38  * This file contains function definitions necessary
  39  * for computations of forces due to restricted angle, restricted dihedral and
  40  * combined bending-torsion potentials.
  41  *
  42  * \author Nicolae Goga
  43  *
  44  * \ingroup module_listed_forces
  45  */
  46 #include "gmxpre.h"
  47
  48 #include "restcbt.h"
  49
  50 #include <cmath>
  51
  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"
  57
  58 /* This function computes factors needed for restricted angle potential.
  59  * For explanations on formula used see file "restcbt.h" */
  60
  61 void compute_factors_restangles(int             type,
  62                                 const t_iparams forceparams[],
  63                                 rvec            delta_ante,
  64                                 rvec            delta_post,
  65                                 double*         prefactor,
  66                                 double*         ratio_ante,
  67                                 double*         ratio_post,
  68                                 real*           v)
  69 {
  70     // These variables are double to make the code
  71     // reproducible.
  72     double theta_equil, k_bending;
  73     double cosine_theta_equil;
  74     double c_ante, c_cros, c_post;
  75     double norm;
  76     double delta_cosine, cosine_theta;
  77     double sine_theta_sq;
  78     double term_theta_theta_equil;
  79
  80     k_bending          = forceparams[type].harmonic.krA;
  81     theta_equil        = forceparams[type].harmonic.rA * gmx::c_deg2Rad;
  82     theta_equil        = M_PI - theta_equil;
  83     cosine_theta_equil = cos(theta_equil);
  84
  85     c_ante = iprod(delta_ante, delta_ante);
  86     c_cros = iprod(delta_ante, delta_post);
  87     c_post = iprod(delta_post, delta_post);
  88
  89     norm          = gmx::invsqrt(c_ante * c_post);
  90     cosine_theta  = c_cros * norm;
  91     sine_theta_sq = 1 - cosine_theta * cosine_theta;
  92
  93     *ratio_ante = c_cros / c_ante;
  94     *ratio_post = c_cros / c_post;
  95
  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);
  99
 100     *v = k_bending * 0.5 * delta_cosine * delta_cosine / sine_theta_sq;
 101 }
 102
 103
 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[],
 108                                rvec            delta_ante,
 109                                rvec            delta_crnt,
 110                                rvec            delta_post,
 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,
 123                                real*           prefactor_phi,
 124                                real*           v)
 125 {
 126
 127     real phi0, cosine_phi0;
 128     real k_torsion;
 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;
 135     real norm_phi;
 136
 137     /* Read parameters phi0 and k_torsion */
 138     phi0        = forceparams[type].pdihs.phiA * gmx::c_deg2Rad;
 139     cosine_phi0 = cos(phi0);
 140     k_torsion   = forceparams[type].pdihs.cpA;
 141
 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;
 155
 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)
 159     {
 160         d_ante = GMX_REAL_EPS;
 161     }
 162     if (d_post < GMX_REAL_EPS)
 163     {
 164         d_post = GMX_REAL_EPS;
 165     }
 166
 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;
 171
 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)
 174     {
 175         sine_phi_sq = 0.0;
 176     }
 177
 178     /* Computation of the differences of cosines (delta_cosine) and a term (term_phi_phi0)
 179      * that is part of the common prefactor_phi */
 180
 181     delta_cosine  = cosine_phi - cosine_phi0;
 182     term_phi_phi0 = 1 - cosine_phi * cosine_phi0;
 183
 184
 185     /*      Computation of ratios */
 186     ratio_phi_ante = c_prod / d_ante;
 187     ratio_phi_post = c_prod / d_post;
 188
 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);
 191
 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" */
 197
 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;
 212
 213     /* Contribution to energy  - see formula in file "restcbt.h"*/
 214     *v = k_torsion * 0.5 * delta_cosine * delta_cosine / sine_phi_sq;
 215 }
 216
 217
 218 /* Compute factors for CBT potential
 219  * For explanations on formula used see file "restcbt.h" */
 220
 221 void compute_factors_cbtdihs(int             type,
 222                              const t_iparams forceparams[],
 223                              rvec            delta_ante,
 224                              rvec            delta_crnt,
 225                              rvec            delta_post,
 226                              rvec            f_phi_ai,
 227                              rvec            f_phi_aj,
 228                              rvec            f_phi_ak,
 229                              rvec            f_phi_al,
 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,
 236                              real*           v)
 237 {
 238     int  j, d;
 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;
 247     real prefactor_phi;
 248     real ratio_phi_ante, ratio_phi_post;
 249     real r1, r2;
 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;
 256
 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 */
 268
 269
 270     /* PART 1 - COMPUTES FORCE FACTORS COMMON TO ALL DERIVATIVES FOR THE FOUR PARTICLES */
 271
 272
 273     for (j = 0; (j < NR_CBTDIHS); j++)
 274     {
 275         torsion_coef[j] = forceparams[type].cbtdihs.cbtcA[j];
 276     }
 277
 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 */
 282
 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;
 292
 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)
 296     {
 297         d_ante = GMX_REAL_EPS;
 298     }
 299     if (d_post < GMX_REAL_EPS)
 300     {
 301         d_post = GMX_REAL_EPS;
 302     }
 303
 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;
 313
 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)
 316     {
 317         sine_theta_ante_sq = 0.0;
 318     }
 319     if (sine_theta_post_sq < 0.0)
 320     {
 321         sine_theta_post_sq = 0.0;
 322     }
 323
 324     sine_theta_ante = std::sqrt(sine_theta_ante_sq);
 325     sine_theta_post = std::sqrt(sine_theta_post_sq);
 326
 327     /* PART 2 - COMPUTES FORCE COMPONENTS DUE TO DERIVATIVES TO DIHEDRAL ANGLE PHI */
 328
 329     /*      Computation of ratios */
 330     ratio_phi_ante = c_prod / d_ante;
 331     ratio_phi_post = c_prod / d_post;
 332
 333     /*       Computation of the prefactor */
 334     /*      Computing 2nd power */
 335     r1 = cosine_phi;
 336
 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;
 341
 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" */
 347
 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;
 362
 363     /* Computation of forces due to the derivatives of dihedral angle phi*/
 364     for (d = 0; d < DIM; d++)
 365     {
 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]);
 378     }
 379
 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;
 384
 385     /*      Computation of the prefactor */
 386     /*      Computing 2nd power */
 387     r1 = cosine_phi;
 388     /*      Computing 3rd power */
 389     r2 = cosine_phi;
 390
 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;
 396
 397
 398     /*      Computation of forces due to the derivatives of bending angle theta_ante */
 399     for (d = 0; d < DIM; d++)
 400     {
 401         f_theta_ante_ai[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]);
 406         f_theta_ante_ak[d] =
 407                 prefactor_theta_ante * (delta_ante[d] - ratio_theta_ante_crnt * delta_crnt[d]);
 408     }
 409
 410     /* PART 4 - COMPUTES THE FORCE COMPONENTS DUE TO THE DERIVATIVES OF THE BENDING ANGLE THETA_POST */
 411
 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;
 415
 416     /*     Computation of the prefactor */
 417     /*      Computing 2nd power */
 418     r1 = cosine_phi;
 419     /*      Computing 3rd power */
 420     r2 = cosine_phi;
 421
 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;
 427
 428
 429     /*      Computation of forces due to the derivatives of bending angle Theta_Post */
 430     for (d = 0; d < DIM; d++)
 431     {
 432         f_theta_post_aj[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]);
 437         f_theta_post_al[d] =
 438                 prefactor_theta_post * (delta_crnt[d] - ratio_theta_post_post * delta_post[d]);
 439     }
 440     r1 = cosine_phi;
 441     r2 = cosine_phi;
 442
 443     /* Contribution to energy - for formula see file "restcbt.h" */
 444     *v = torsion_coef[0]
 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;
 448 }