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