src/gromacs/gmxlib/restcbt.c

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