src/gromacs/mdlib/shakef.cpp

   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) 2013,2014,2015, 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 #include "gmxpre.h"
  38
  39 #include <math.h>
  40
  41 #include "gromacs/gmxlib/nrnb.h"
  42 #include "gromacs/math/functions.h"
  43 #include "gromacs/math/vec.h"
  44 #include "gromacs/mdlib/constr.h"
  45 #include "gromacs/mdtypes/inputrec.h"
  46 #include "gromacs/mdtypes/md_enums.h"
  47 #include "gromacs/utility/smalloc.h"
  48
  49 typedef struct gmx_shakedata
  50 {
  51     rvec *rij;
  52     real *half_of_reduced_mass;
  53     real *distance_squared_tolerance;
  54     real *constraint_distance_squared;
  55     int   nalloc;
  56     /* SOR stuff */
  57     real  delta;
  58     real  omega;
  59     real  gamma;
  60 } t_gmx_shakedata;
  61
  62 gmx_shakedata_t shake_init()
  63 {
  64     gmx_shakedata_t d;
  65
  66     snew(d, 1);
  67
  68     d->nalloc                      = 0;
  69     d->rij                         = NULL;
  70     d->half_of_reduced_mass        = NULL;
  71     d->distance_squared_tolerance  = NULL;
  72     d->constraint_distance_squared = NULL;
  73
  74     /* SOR initialization */
  75     d->delta = 0.1;
  76     d->omega = 1.0;
  77     d->gamma = 1000000;
  78
  79     return d;
  80 }
  81
  82 /*! \brief Inner kernel for SHAKE constraints
  83  *
  84  * Original implementation from R.C. van Schaik and W.F. van Gunsteren
  85  * (ETH Zuerich, June 1992), adapted for GROMACS by David van der
  86  * Spoel November 1992.
  87  *
  88  * The algorithm here is based section five of Ryckaert, Ciccotti and
  89  * Berendsen, J Comp Phys, 23, 327, 1977.
  90  *
  91  * \param[in]    iatom                         Mini-topology of triples of constraint type (unused in this
  92  *                                             function) and indices of the two atoms involved
  93  * \param[in]    ncon                          Number of constraints
  94  * \param[out]   nnit                          Number of iterations performed
  95  * \param[in]    maxnit                        Maximum number of iterations permitted
  96  * \param[in]    constraint_distance_squared   The objective value for each constraint
  97  * \param[inout] positions                     The initial (and final) values of the positions of all atoms
  98  * \param[in]    initial_displacements         The initial displacements of each constraint
  99  * \param[in]    half_of_reduced_mass          Half of the reduced mass for each constraint
 100  * \param[in]    omega                         SHAKE over-relaxation factor (set non-1.0 by
 101  *                                             using shake-sor=yes in the .mdp, but there is no documentation anywhere)
 102  * \param[in]    invmass                       Inverse mass of each atom
 103  * \param[in]    distance_squared_tolerance    Multiplicative tolerance on the difference in the
 104  *                                             square of the constrained distance (see code)
 105  * \param[out]   scaled_lagrange_multiplier    Scaled Lagrange multiplier for each constraint (-2 * eta from p. 336
 106  *                                             of the paper, divided by the constraint distance)
 107  * \param[out]   nerror                        Zero upon success, returns one more than the index of the
 108  *                                             problematic constraint if the input was malformed
 109  *
 110  * \todo Make SHAKE use better data structures, in particular for iatom. */
 111 void cshake(const int iatom[], int ncon, int *nnit, int maxnit,
 112             const real constraint_distance_squared[], real positions[],
 113             const real initial_displacements[], const real half_of_reduced_mass[], real omega,
 114             const real invmass[], const real distance_squared_tolerance[],
 115             real scaled_lagrange_multiplier[], int *nerror)
 116 {
 117     /* default should be increased! MRS 8/4/2009 */
 118     const real mytol = 1e-10;
 119
 120     int        ll, i, j, i3, j3, l3;
 121     int        ix, iy, iz, jx, jy, jz;
 122     real       r_dot_r_prime;
 123     real       constraint_distance_squared_ll;
 124     real       r_prime_squared;
 125     real       scaled_lagrange_multiplier_ll;
 126     real       r_prime_x, r_prime_y, r_prime_z, diff, im, jm;
 127     real       xh, yh, zh, rijx, rijy, rijz;
 128     int        nit, error, nconv;
 129     real       iconvf;
 130
 131     // TODO nconv is used solely as a boolean, so we should write the
 132     // code like that
 133     error = 0;
 134     nconv = 1;
 135     for (nit = 0; (nit < maxnit) && (nconv != 0) && (error == 0); nit++)
 136     {
 137         nconv = 0;
 138         for (ll = 0; (ll < ncon) && (error == 0); ll++)
 139         {
 140             l3    = 3*ll;
 141             rijx  = initial_displacements[l3+XX];
 142             rijy  = initial_displacements[l3+YY];
 143             rijz  = initial_displacements[l3+ZZ];
 144             i     = iatom[l3+1];
 145             j     = iatom[l3+2];
 146             i3    = 3*i;
 147             j3    = 3*j;
 148             ix    = i3+XX;
 149             iy    = i3+YY;
 150             iz    = i3+ZZ;
 151             jx    = j3+XX;
 152             jy    = j3+YY;
 153             jz    = j3+ZZ;
 154
 155             /* Compute r prime between atoms i and j, which is the
 156                displacement *before* this update stage */
 157             r_prime_x       = positions[ix]-positions[jx];
 158             r_prime_y       = positions[iy]-positions[jy];
 159             r_prime_z       = positions[iz]-positions[jz];
 160             r_prime_squared = (r_prime_x * r_prime_x +
 161                                r_prime_y * r_prime_y +
 162                                r_prime_z * r_prime_z);
 163             constraint_distance_squared_ll = constraint_distance_squared[ll];
 164             diff    = constraint_distance_squared_ll - r_prime_squared;
 165
 166             /* iconvf is less than 1 when the error is smaller than a bound */
 167             iconvf = fabs(diff) * distance_squared_tolerance[ll];
 168
 169             if (iconvf > 1.0)
 170             {
 171                 nconv         = static_cast<int>(iconvf);
 172                 r_dot_r_prime = (rijx * r_prime_x +
 173                                  rijy * r_prime_y +
 174                                  rijz * r_prime_z);
 175
 176                 if (r_dot_r_prime < constraint_distance_squared_ll * mytol)
 177                 {
 178                     error = ll+1;
 179                 }
 180                 else
 181                 {
 182                     /* The next line solves equation 5.6 (neglecting
 183                        the term in g^2), for g */
 184                     scaled_lagrange_multiplier_ll   = omega*diff*half_of_reduced_mass[ll]/r_dot_r_prime;
 185                     scaled_lagrange_multiplier[ll] += scaled_lagrange_multiplier_ll;
 186                     xh                              = rijx * scaled_lagrange_multiplier_ll;
 187                     yh                              = rijy * scaled_lagrange_multiplier_ll;
 188                     zh                              = rijz * scaled_lagrange_multiplier_ll;
 189                     im                              = invmass[i];
 190                     jm                              = invmass[j];
 191                     positions[ix]                  += xh*im;
 192                     positions[iy]                  += yh*im;
 193                     positions[iz]                  += zh*im;
 194                     positions[jx]                  -= xh*jm;
 195                     positions[jy]                  -= yh*jm;
 196                     positions[jz]                  -= zh*jm;
 197                 }
 198             }
 199         }
 200     }
 201     *nnit   = nit;
 202     *nerror = error;
 203 }
 204
 205 int vec_shakef(FILE *fplog, gmx_shakedata_t shaked,
 206                real invmass[], int ncon,
 207                t_iparams ip[], t_iatom *iatom,
 208                real tol, rvec x[], rvec prime[], real omega,
 209                gmx_bool bFEP, real lambda, real scaled_lagrange_multiplier[],
 210                real invdt, rvec *v,
 211                gmx_bool bCalcVir, tensor vir_r_m_dr, int econq)
 212 {
 213     rvec    *rij;
 214     real    *half_of_reduced_mass, *distance_squared_tolerance, *constraint_distance_squared;
 215     int      maxnit = 1000;
 216     int      nit    = 0, ll, i, j, d, d2, type;
 217     t_iatom *ia;
 218     real     L1;
 219     real     mm    = 0., tmp;
 220     int      error = 0;
 221     real     constraint_distance;
 222
 223     if (ncon > shaked->nalloc)
 224     {
 225         shaked->nalloc = over_alloc_dd(ncon);
 226         srenew(shaked->rij, shaked->nalloc);
 227         srenew(shaked->half_of_reduced_mass, shaked->nalloc);
 228         srenew(shaked->distance_squared_tolerance, shaked->nalloc);
 229         srenew(shaked->constraint_distance_squared, shaked->nalloc);
 230     }
 231     rij                          = shaked->rij;
 232     half_of_reduced_mass         = shaked->half_of_reduced_mass;
 233     distance_squared_tolerance   = shaked->distance_squared_tolerance;
 234     constraint_distance_squared  = shaked->constraint_distance_squared;
 235
 236     L1   = 1.0-lambda;
 237     ia   = iatom;
 238     for (ll = 0; (ll < ncon); ll++, ia += 3)
 239     {
 240         type  = ia[0];
 241         i     = ia[1];
 242         j     = ia[2];
 243
 244         mm                       = 2.0*(invmass[i]+invmass[j]);
 245         rij[ll][XX]              = x[i][XX]-x[j][XX];
 246         rij[ll][YY]              = x[i][YY]-x[j][YY];
 247         rij[ll][ZZ]              = x[i][ZZ]-x[j][ZZ];
 248         half_of_reduced_mass[ll] = 1.0/mm;
 249         if (bFEP)
 250         {
 251             constraint_distance = L1*ip[type].constr.dA + lambda*ip[type].constr.dB;
 252         }
 253         else
 254         {
 255             constraint_distance = ip[type].constr.dA;
 256         }
 257         constraint_distance_squared[ll]  = gmx::square(constraint_distance);
 258         distance_squared_tolerance[ll]   = 0.5/(constraint_distance_squared[ll]*tol);
 259     }
 260
 261     switch (econq)
 262     {
 263         case econqCoord:
 264             cshake(iatom, ncon, &nit, maxnit, constraint_distance_squared, prime[0], rij[0], half_of_reduced_mass, omega, invmass, distance_squared_tolerance, scaled_lagrange_multiplier, &error);
 265             break;
 266         case econqVeloc:
 267             crattle(iatom, ncon, &nit, maxnit, constraint_distance_squared, prime[0], rij[0], half_of_reduced_mass, omega, invmass, distance_squared_tolerance, scaled_lagrange_multiplier, &error, invdt);
 268             break;
 269     }
 270
 271     if (nit >= maxnit)
 272     {
 273         if (fplog)
 274         {
 275             fprintf(fplog, "Shake did not converge in %d steps\n", maxnit);
 276         }
 277         fprintf(stderr, "Shake did not converge in %d steps\n", maxnit);
 278         nit = 0;
 279     }
 280     else if (error != 0)
 281     {
 282         if (fplog)
 283         {
 284             fprintf(fplog, "Inner product between old and new vector <= 0.0!\n"
 285                     "constraint #%d atoms %d and %d\n",
 286                     error-1, iatom[3*(error-1)+1]+1, iatom[3*(error-1)+2]+1);
 287         }
 288         fprintf(stderr, "Inner product between old and new vector <= 0.0!\n"
 289                 "constraint #%d atoms %d and %d\n",
 290                 error-1, iatom[3*(error-1)+1]+1, iatom[3*(error-1)+2]+1);
 291         nit = 0;
 292     }
 293
 294     /* Constraint virial and correct the Lagrange multipliers for the length */
 295
 296     ia = iatom;
 297
 298     for (ll = 0; (ll < ncon); ll++, ia += 3)
 299     {
 300         type  = ia[0];
 301         i     = ia[1];
 302         j     = ia[2];
 303
 304         if ((econq == econqCoord) && v != NULL)
 305         {
 306             /* Correct the velocities */
 307             mm = scaled_lagrange_multiplier[ll]*invmass[i]*invdt;
 308             for (d = 0; d < DIM; d++)
 309             {
 310                 v[ia[1]][d] += mm*rij[ll][d];
 311             }
 312             mm = scaled_lagrange_multiplier[ll]*invmass[j]*invdt;
 313             for (d = 0; d < DIM; d++)
 314             {
 315                 v[ia[2]][d] -= mm*rij[ll][d];
 316             }
 317             /* 16 flops */
 318         }
 319
 320         /* constraint virial */
 321         if (bCalcVir)
 322         {
 323             mm = scaled_lagrange_multiplier[ll];
 324             for (d = 0; d < DIM; d++)
 325             {
 326                 tmp = mm*rij[ll][d];
 327                 for (d2 = 0; d2 < DIM; d2++)
 328                 {
 329                     vir_r_m_dr[d][d2] -= tmp*rij[ll][d2];
 330                 }
 331             }
 332             /* 21 flops */
 333         }
 334
 335         /* cshake and crattle produce Lagrange multipliers scaled by
 336            the reciprocal of the constraint length, so fix that */
 337         if (bFEP)
 338         {
 339             constraint_distance = L1*ip[type].constr.dA + lambda*ip[type].constr.dB;
 340         }
 341         else
 342         {
 343             constraint_distance = ip[type].constr.dA;
 344         }
 345         scaled_lagrange_multiplier[ll] *= constraint_distance;
 346     }
 347
 348     return nit;
 349 }
 350
 351 static void check_cons(FILE *log, int nc, rvec x[], rvec prime[], rvec v[],
 352                        t_iparams ip[], t_iatom *iatom,
 353                        real invmass[], int econq)
 354 {
 355     t_iatom *ia;
 356     int      ai, aj;
 357     int      i;
 358     real     d, dp;
 359     rvec     dx, dv;
 360
 361     fprintf(log,
 362             "    i     mi      j     mj      before       after   should be\n");
 363     ia = iatom;
 364     for (i = 0; (i < nc); i++, ia += 3)
 365     {
 366         ai = ia[1];
 367         aj = ia[2];
 368         rvec_sub(x[ai], x[aj], dx);
 369         d = norm(dx);
 370
 371         switch (econq)
 372         {
 373             case econqCoord:
 374                 rvec_sub(prime[ai], prime[aj], dx);
 375                 dp = norm(dx);
 376                 fprintf(log, "%5d  %5.2f  %5d  %5.2f  %10.5f  %10.5f  %10.5f\n",
 377                         ai+1, 1.0/invmass[ai],
 378                         aj+1, 1.0/invmass[aj], d, dp, ip[ia[0]].constr.dA);
 379                 break;
 380             case econqVeloc:
 381                 rvec_sub(v[ai], v[aj], dv);
 382                 d = iprod(dx, dv);
 383                 rvec_sub(prime[ai], prime[aj], dv);
 384                 dp = iprod(dx, dv);
 385                 fprintf(log, "%5d  %5.2f  %5d  %5.2f  %10.5f  %10.5f  %10.5f\n",
 386                         ai+1, 1.0/invmass[ai],
 387                         aj+1, 1.0/invmass[aj], d, dp, 0.);
 388                 break;
 389         }
 390     }
 391 }
 392
 393 gmx_bool bshakef(FILE *log, gmx_shakedata_t shaked,
 394                  real invmass[], int nblocks, int sblock[],
 395                  t_idef *idef, t_inputrec *ir, rvec x_s[], rvec prime[],
 396                  t_nrnb *nrnb, real *scaled_lagrange_multiplier, real lambda, real *dvdlambda,
 397                  real invdt, rvec *v, gmx_bool bCalcVir, tensor vir_r_m_dr,
 398                  gmx_bool bDumpOnError, int econq)
 399 {
 400     t_iatom *iatoms;
 401     real     dt_2, dvdl;
 402     int      i, n0, ncon, blen, type, ll;
 403     int      tnit = 0, trij = 0;
 404
 405 #ifdef DEBUG
 406     fprintf(log, "nblocks=%d, sblock[0]=%d\n", nblocks, sblock[0]);
 407 #endif
 408
 409     ncon = idef->il[F_CONSTR].nr/3;
 410
 411     for (ll = 0; ll < ncon; ll++)
 412     {
 413         scaled_lagrange_multiplier[ll] = 0;
 414     }
 415
 416     iatoms = &(idef->il[F_CONSTR].iatoms[sblock[0]]);
 417     for (i = 0; (i < nblocks); )
 418     {
 419         blen  = (sblock[i+1]-sblock[i]);
 420         blen /= 3;
 421         n0    = vec_shakef(log, shaked, invmass, blen, idef->iparams,
 422                            iatoms, ir->shake_tol, x_s, prime, shaked->omega,
 423                            ir->efep != efepNO, lambda, scaled_lagrange_multiplier, invdt, v, bCalcVir, vir_r_m_dr,
 424                            econq);
 425
 426 #ifdef DEBUGSHAKE
 427         check_cons(log, blen, x_s, prime, v, idef->iparams, iatoms, invmass, econq);
 428 #endif
 429
 430         if (n0 == 0)
 431         {
 432             if (bDumpOnError && log)
 433             {
 434                 {
 435                     check_cons(log, blen, x_s, prime, v, idef->iparams, iatoms, invmass, econq);
 436                 }
 437             }
 438             return FALSE;
 439         }
 440         tnit                       += n0*blen;
 441         trij                       += blen;
 442         iatoms                     += 3*blen; /* Increment pointer! */
 443         scaled_lagrange_multiplier += blen;
 444         i++;
 445     }
 446     /* only for position part? */
 447     if (econq == econqCoord)
 448     {
 449         if (ir->efep != efepNO)
 450         {
 451             real bondA, bondB;
 452             /* TODO This should probably use invdt, so that sd integrator scaling works properly */
 453             dt_2 = 1/gmx::square(ir->delta_t);
 454             dvdl = 0;
 455             for (ll = 0; ll < ncon; ll++)
 456             {
 457                 type  = idef->il[F_CONSTR].iatoms[3*ll];
 458
 459                 /* Per equations in the manual, dv/dl = -2 \sum_ll lagrangian_ll * r_ll * (d_B - d_A) */
 460                 /* The vector scaled_lagrange_multiplier[ll] contains the value -2 r_ll eta_ll (eta_ll is the
 461                    estimate of the Langrangian, definition on page 336 of Ryckaert et al 1977),
 462                    so the pre-factors are already present. */
 463                 bondA = idef->iparams[type].constr.dA;
 464                 bondB = idef->iparams[type].constr.dB;
 465                 dvdl += scaled_lagrange_multiplier[ll] * dt_2 * (bondB - bondA);
 466             }
 467             *dvdlambda += dvdl;
 468         }
 469     }
 470 #ifdef DEBUG
 471     fprintf(log, "tnit: %5d  omega: %10.5f\n", tnit, omega);
 472 #endif
 473     if (ir->bShakeSOR)
 474     {
 475         if (tnit > shaked->gamma)
 476         {
 477             shaked->delta *= -0.5;
 478         }
 479         shaked->omega += shaked->delta;
 480         shaked->gamma  = tnit;
 481     }
 482     inc_nrnb(nrnb, eNR_SHAKE, tnit);
 483     inc_nrnb(nrnb, eNR_SHAKE_RIJ, trij);
 484     if (v)
 485     {
 486         inc_nrnb(nrnb, eNR_CONSTR_V, trij*2);
 487     }
 488     if (bCalcVir)
 489     {
 490         inc_nrnb(nrnb, eNR_CONSTR_VIR, trij);
 491     }
 492
 493     return TRUE;
 494 }
 495
 496 void crattle(int iatom[], int ncon, int *nnit, int maxnit,
 497              real constraint_distance_squared[], real vp[], real rij[], real m2[], real omega,
 498              real invmass[], real distance_squared_tolerance[], real scaled_lagrange_multiplier[],
 499              int *nerror, real invdt)
 500 {
 501     /*
 502      *     r.c. van schaik and w.f. van gunsteren
 503      *     eth zuerich
 504      *     june 1992
 505      *     Adapted for use with Gromacs by David van der Spoel november 92 and later.
 506      *     rattle added by M.R. Shirts, April 2004, from code written by Jay Ponder in TINKER
 507      *     second part of rattle algorithm
 508      */
 509
 510     int          ll, i, j, i3, j3, l3;
 511     int          ix, iy, iz, jx, jy, jz;
 512     real         constraint_distance_squared_ll;
 513     real         vpijd, vx, vy, vz, acor, fac, im, jm;
 514     real         xh, yh, zh, rijx, rijy, rijz;
 515     int          nit, error, nconv;
 516     real         iconvf;
 517
 518     // TODO nconv is used solely as a boolean, so we should write the
 519     // code like that
 520     error = 0;
 521     nconv = 1;
 522     for (nit = 0; (nit < maxnit) && (nconv != 0) && (error == 0); nit++)
 523     {
 524         nconv = 0;
 525         for (ll = 0; (ll < ncon) && (error == 0); ll++)
 526         {
 527             l3      = 3*ll;
 528             rijx    = rij[l3+XX];
 529             rijy    = rij[l3+YY];
 530             rijz    = rij[l3+ZZ];
 531             i       = iatom[l3+1];
 532             j       = iatom[l3+2];
 533             i3      = 3*i;
 534             j3      = 3*j;
 535             ix      = i3+XX;
 536             iy      = i3+YY;
 537             iz      = i3+ZZ;
 538             jx      = j3+XX;
 539             jy      = j3+YY;
 540             jz      = j3+ZZ;
 541             vx      = vp[ix]-vp[jx];
 542             vy      = vp[iy]-vp[jy];
 543             vz      = vp[iz]-vp[jz];
 544
 545             vpijd   = vx*rijx+vy*rijy+vz*rijz;
 546             constraint_distance_squared_ll = constraint_distance_squared[ll];
 547
 548             /* iconv is zero when the error is smaller than a bound */
 549             iconvf   = fabs(vpijd)*(distance_squared_tolerance[ll]/invdt);
 550
 551             if (iconvf > 1)
 552             {
 553                 nconv     = static_cast<int>(iconvf);
 554                 fac       = omega*2.0*m2[ll]/constraint_distance_squared_ll;
 555                 acor      = -fac*vpijd;
 556                 scaled_lagrange_multiplier[ll] += acor;
 557                 xh        = rijx*acor;
 558                 yh        = rijy*acor;
 559                 zh        = rijz*acor;
 560
 561                 im        = invmass[i];
 562                 jm        = invmass[j];
 563
 564                 vp[ix] += xh*im;
 565                 vp[iy] += yh*im;
 566                 vp[iz] += zh*im;
 567                 vp[jx] -= xh*jm;
 568                 vp[jy] -= yh*jm;
 569                 vp[jz] -= zh*jm;
 570             }
 571         }
 572     }
 573     *nnit   = nit;
 574     *nerror = error;
 575 }