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.
   7  * Copyright (c) 2019,2020, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 /*! \internal \file
  39  * \brief Defines SHAKE code.
  40  *
  41  * \author David van der Spoel <david.vanderspoel@icm.uu.se>
  42  * \author Berk Hess <hess@kth.se>
  43  * \author Mark Abraham <mark.j.abraham@gmail.com>
  44  * \ingroup module_mdlib
  45  */
  46 #include "gmxpre.h"
  47
  48 #include "shake.h"
  49
  50 #include <cmath>
  51
  52 #include <algorithm>
  53
  54 #include "gromacs/domdec/domdec_struct.h"
  55 #include "gromacs/gmxlib/nrnb.h"
  56 #include "gromacs/math/functions.h"
  57 #include "gromacs/math/vec.h"
  58 #include "gromacs/mdlib/constr.h"
  59 #include "gromacs/mdlib/splitter.h"
  60 #include "gromacs/mdtypes/inputrec.h"
  61 #include "gromacs/mdtypes/md_enums.h"
  62 #include "gromacs/mdtypes/mdatom.h"
  63 #include "gromacs/topology/invblock.h"
  64 #include "gromacs/utility/fatalerror.h"
  65 #include "gromacs/utility/smalloc.h"
  66
  67 namespace gmx
  68 {
  69
  70 struct shakedata
  71 {
  72     rvec* rij;
  73     real* half_of_reduced_mass;
  74     real* distance_squared_tolerance;
  75     real* constraint_distance_squared;
  76     int   nalloc;
  77     /* SOR stuff */
  78     real delta;
  79     real omega;
  80     real gamma;
  81     int  nblocks;       /* The number of SHAKE blocks         */
  82     int* sblock;        /* The SHAKE blocks                   */
  83     int  sblock_nalloc; /* The allocation size of sblock      */
  84     /*! \brief Scaled Lagrange multiplier for each constraint.
  85      *
  86      * Value is -2 * eta from p. 336 of the paper, divided by the
  87      * constraint distance. */
  88     real* scaled_lagrange_multiplier;
  89     int   lagr_nalloc; /* The allocation size of scaled_lagrange_multiplier */
  90 };
  91
  92 shakedata* shake_init()
  93 {
  94     shakedata* d;
  95
  96     snew(d, 1);
  97
  98     d->nalloc                      = 0;
  99     d->rij                         = nullptr;
 100     d->half_of_reduced_mass        = nullptr;
 101     d->distance_squared_tolerance  = nullptr;
 102     d->constraint_distance_squared = nullptr;
 103
 104     /* SOR initialization */
 105     d->delta = 0.1;
 106     d->omega = 1.0;
 107     d->gamma = 1000000;
 108
 109     return d;
 110 }
 111
 112 void done_shake(shakedata* d)
 113 {
 114     sfree(d->rij);
 115     sfree(d->half_of_reduced_mass);
 116     sfree(d->distance_squared_tolerance);
 117     sfree(d->constraint_distance_squared);
 118     sfree(d->sblock);
 119     sfree(d->scaled_lagrange_multiplier);
 120     sfree(d);
 121 }
 122
 123 typedef struct
 124 {
 125     int iatom[3];
 126     int blocknr;
 127 } t_sortblock;
 128
 129 //! Compares sort blocks.
 130 static int pcomp(const void* p1, const void* p2)
 131 {
 132     int                db;
 133     int                min1, min2, max1, max2;
 134     const t_sortblock* a1 = reinterpret_cast<const t_sortblock*>(p1);
 135     const t_sortblock* a2 = reinterpret_cast<const t_sortblock*>(p2);
 136
 137     db = a1->blocknr - a2->blocknr;
 138
 139     if (db != 0)
 140     {
 141         return db;
 142     }
 143
 144     min1 = std::min(a1->iatom[1], a1->iatom[2]);
 145     max1 = std::max(a1->iatom[1], a1->iatom[2]);
 146     min2 = std::min(a2->iatom[1], a2->iatom[2]);
 147     max2 = std::max(a2->iatom[1], a2->iatom[2]);
 148
 149     if (min1 == min2)
 150     {
 151         return max1 - max2;
 152     }
 153     else
 154     {
 155         return min1 - min2;
 156     }
 157 }
 158
 159 //! Prints sortblocks
 160 static void pr_sortblock(FILE* fp, const char* title, int nsb, t_sortblock sb[])
 161 {
 162     int i;
 163
 164     fprintf(fp, "%s\n", title);
 165     for (i = 0; (i < nsb); i++)
 166     {
 167         fprintf(fp, "i: %5d, iatom: (%5d %5d %5d), blocknr: %5d\n", i, sb[i].iatom[0],
 168                 sb[i].iatom[1], sb[i].iatom[2], sb[i].blocknr);
 169     }
 170 }
 171
 172 //! Reallocates a vector.
 173 static void resizeLagrangianData(shakedata* shaked, int ncons)
 174 {
 175     if (ncons > shaked->lagr_nalloc)
 176     {
 177         shaked->lagr_nalloc = over_alloc_dd(ncons);
 178         srenew(shaked->scaled_lagrange_multiplier, shaked->lagr_nalloc);
 179     }
 180 }
 181
 182 void make_shake_sblock_serial(shakedata* shaked, InteractionDefinitions* idef, const t_mdatoms& md)
 183 {
 184     int          i, j, m, ncons;
 185     int          bstart, bnr;
 186     t_blocka     sblocks;
 187     t_sortblock* sb;
 188     int*         inv_sblock;
 189
 190     /* Since we are processing the local topology,
 191      * the F_CONSTRNC ilist has been concatenated to the F_CONSTR ilist.
 192      */
 193     ncons = idef->il[F_CONSTR].size() / 3;
 194
 195     init_blocka(&sblocks);
 196     sfree(sblocks.index); // To solve memory leak
 197     gen_sblocks(nullptr, 0, md.homenr, *idef, &sblocks, FALSE);
 198
 199     /*
 200        bstart=(idef->nodeid > 0) ? blocks->multinr[idef->nodeid-1] : 0;
 201        nblocks=blocks->multinr[idef->nodeid] - bstart;
 202      */
 203     bstart          = 0;
 204     shaked->nblocks = sblocks.nr;
 205     if (debug)
 206     {
 207         fprintf(debug, "ncons: %d, bstart: %d, nblocks: %d\n", ncons, bstart, shaked->nblocks);
 208     }
 209
 210     /* Calculate block number for each atom */
 211     inv_sblock = make_invblocka(&sblocks, md.nr);
 212
 213     done_blocka(&sblocks);
 214
 215     /* Store the block number in temp array and
 216      * sort the constraints in order of the sblock number
 217      * and the atom numbers, really sorting a segment of the array!
 218      */
 219     int* iatom = idef->il[F_CONSTR].iatoms.data();
 220     snew(sb, ncons);
 221     for (i = 0; (i < ncons); i++, iatom += 3)
 222     {
 223         for (m = 0; (m < 3); m++)
 224         {
 225             sb[i].iatom[m] = iatom[m];
 226         }
 227         sb[i].blocknr = inv_sblock[iatom[1]];
 228     }
 229
 230     /* Now sort the blocks */
 231     if (debug)
 232     {
 233         pr_sortblock(debug, "Before sorting", ncons, sb);
 234         fprintf(debug, "Going to sort constraints\n");
 235     }
 236
 237     std::qsort(sb, ncons, sizeof(*sb), pcomp);
 238
 239     if (debug)
 240     {
 241         pr_sortblock(debug, "After sorting", ncons, sb);
 242     }
 243
 244     iatom = idef->il[F_CONSTR].iatoms.data();
 245     for (i = 0; (i < ncons); i++, iatom += 3)
 246     {
 247         for (m = 0; (m < 3); m++)
 248         {
 249             iatom[m] = sb[i].iatom[m];
 250         }
 251     }
 252
 253     j = 0;
 254     snew(shaked->sblock, shaked->nblocks + 1);
 255     bnr = -2;
 256     for (i = 0; (i < ncons); i++)
 257     {
 258         if (sb[i].blocknr != bnr)
 259         {
 260             bnr                 = sb[i].blocknr;
 261             shaked->sblock[j++] = 3 * i;
 262         }
 263     }
 264     /* Last block... */
 265     shaked->sblock[j++] = 3 * ncons;
 266
 267     if (j != (shaked->nblocks + 1))
 268     {
 269         fprintf(stderr, "bstart: %d\n", bstart);
 270         fprintf(stderr, "j: %d, nblocks: %d, ncons: %d\n", j, shaked->nblocks, ncons);
 271         for (i = 0; (i < ncons); i++)
 272         {
 273             fprintf(stderr, "i: %5d  sb[i].blocknr: %5d\n", i, sb[i].blocknr);
 274         }
 275         for (j = 0; (j <= shaked->nblocks); j++)
 276         {
 277             fprintf(stderr, "sblock[%3d]=%5d\n", j, shaked->sblock[j]);
 278         }
 279         gmx_fatal(FARGS,
 280                   "DEATH HORROR: "
 281                   "sblocks does not match idef->il[F_CONSTR]");
 282     }
 283     sfree(sb);
 284     sfree(inv_sblock);
 285     resizeLagrangianData(shaked, ncons);
 286 }
 287
 288 // TODO: Check if this code is useful. It might never be called.
 289 void make_shake_sblock_dd(shakedata* shaked, const InteractionList& ilcon, const gmx_domdec_t* dd)
 290 {
 291     int ncons, c, cg;
 292
 293     if (dd->ncg_home + 1 > shaked->sblock_nalloc)
 294     {
 295         shaked->sblock_nalloc = over_alloc_dd(dd->ncg_home + 1);
 296         srenew(shaked->sblock, shaked->sblock_nalloc);
 297     }
 298
 299     ncons            = ilcon.size() / 3;
 300     const int* iatom = ilcon.iatoms.data();
 301     shaked->nblocks  = 0;
 302     cg               = 0;
 303     for (c = 0; c < ncons; c++)
 304     {
 305         if (c == 0 || iatom[1] >= cg + 1)
 306         {
 307             shaked->sblock[shaked->nblocks++] = 3 * c;
 308             while (iatom[1] >= cg + 1)
 309             {
 310                 cg++;
 311             }
 312         }
 313         iatom += 3;
 314     }
 315     shaked->sblock[shaked->nblocks] = 3 * ncons;
 316     resizeLagrangianData(shaked, ncons);
 317 }
 318
 319 /*! \brief Inner kernel for SHAKE constraints
 320  *
 321  * Original implementation from R.C. van Schaik and W.F. van Gunsteren
 322  * (ETH Zuerich, June 1992), adapted for GROMACS by David van der
 323  * Spoel November 1992.
 324  *
 325  * The algorithm here is based section five of Ryckaert, Ciccotti and
 326  * Berendsen, J Comp Phys, 23, 327, 1977.
 327  *
 328  * \param[in]    iatom                         Mini-topology of triples of constraint type (unused in this
 329  *                                             function) and indices of the two atoms involved
 330  * \param[in]    ncon                          Number of constraints
 331  * \param[out]   nnit                          Number of iterations performed
 332  * \param[in]    maxnit                        Maximum number of iterations permitted
 333  * \param[in]    constraint_distance_squared   The objective value for each constraint
 334  * \param[inout] positions                     The initial (and final) values of the positions of all atoms
 335  * \param[in]    initial_displacements         The initial displacements of each constraint
 336  * \param[in]    half_of_reduced_mass          Half of the reduced mass for each constraint
 337  * \param[in]    omega                         SHAKE over-relaxation factor (set non-1.0 by
 338  *                                             using shake-sor=yes in the .mdp, but there is no documentation anywhere)
 339  * \param[in]    invmass                       Inverse mass of each atom
 340  * \param[in]    distance_squared_tolerance    Multiplicative tolerance on the difference in the
 341  *                                             square of the constrained distance (see code)
 342  * \param[out]   scaled_lagrange_multiplier    Scaled Lagrange multiplier for each constraint (-2 * eta from p. 336
 343  *                                             of the paper, divided by the constraint distance)
 344  * \param[out]   nerror                        Zero upon success, returns one more than the index of the
 345  *                                             problematic constraint if the input was malformed
 346  *
 347  * \todo Make SHAKE use better data structures, in particular for iatom. */
 348 void cshake(const int  iatom[],
 349             int        ncon,
 350             int*       nnit,
 351             int        maxnit,
 352             const real constraint_distance_squared[],
 353             real       positions[],
 354             const real initial_displacements[],
 355             const real half_of_reduced_mass[],
 356             real       omega,
 357             const real invmass[],
 358             const real distance_squared_tolerance[],
 359             real       scaled_lagrange_multiplier[],
 360             int*       nerror)
 361 {
 362     /* default should be increased! MRS 8/4/2009 */
 363     const real mytol = 1e-10;
 364
 365     int  ll, i, j, i3, j3, l3;
 366     int  ix, iy, iz, jx, jy, jz;
 367     real r_dot_r_prime;
 368     real constraint_distance_squared_ll;
 369     real r_prime_squared;
 370     real scaled_lagrange_multiplier_ll;
 371     real r_prime_x, r_prime_y, r_prime_z, diff, im, jm;
 372     real xh, yh, zh, rijx, rijy, rijz;
 373     int  nit, error, nconv;
 374     real iconvf;
 375
 376     // TODO nconv is used solely as a boolean, so we should write the
 377     // code like that
 378     error = 0;
 379     nconv = 1;
 380     for (nit = 0; (nit < maxnit) && (nconv != 0) && (error == 0); nit++)
 381     {
 382         nconv = 0;
 383         for (ll = 0; (ll < ncon) && (error == 0); ll++)
 384         {
 385             l3   = 3 * ll;
 386             rijx = initial_displacements[l3 + XX];
 387             rijy = initial_displacements[l3 + YY];
 388             rijz = initial_displacements[l3 + ZZ];
 389             i    = iatom[l3 + 1];
 390             j    = iatom[l3 + 2];
 391             i3   = 3 * i;
 392             j3   = 3 * j;
 393             ix   = i3 + XX;
 394             iy   = i3 + YY;
 395             iz   = i3 + ZZ;
 396             jx   = j3 + XX;
 397             jy   = j3 + YY;
 398             jz   = j3 + ZZ;
 399
 400             /* Compute r prime between atoms i and j, which is the
 401                displacement *before* this update stage */
 402             r_prime_x       = positions[ix] - positions[jx];
 403             r_prime_y       = positions[iy] - positions[jy];
 404             r_prime_z       = positions[iz] - positions[jz];
 405             r_prime_squared = (r_prime_x * r_prime_x + r_prime_y * r_prime_y + r_prime_z * r_prime_z);
 406             constraint_distance_squared_ll = constraint_distance_squared[ll];
 407             diff                           = constraint_distance_squared_ll - r_prime_squared;
 408
 409             /* iconvf is less than 1 when the error is smaller than a bound */
 410             iconvf = fabs(diff) * distance_squared_tolerance[ll];
 411
 412             if (iconvf > 1.0)
 413             {
 414                 nconv         = static_cast<int>(iconvf);
 415                 r_dot_r_prime = (rijx * r_prime_x + rijy * r_prime_y + rijz * r_prime_z);
 416
 417                 if (r_dot_r_prime < constraint_distance_squared_ll * mytol)
 418                 {
 419                     error = ll + 1;
 420                 }
 421                 else
 422                 {
 423                     /* The next line solves equation 5.6 (neglecting
 424                        the term in g^2), for g */
 425                     scaled_lagrange_multiplier_ll = omega * diff * half_of_reduced_mass[ll] / r_dot_r_prime;
 426                     scaled_lagrange_multiplier[ll] += scaled_lagrange_multiplier_ll;
 427                     xh = rijx * scaled_lagrange_multiplier_ll;
 428                     yh = rijy * scaled_lagrange_multiplier_ll;
 429                     zh = rijz * scaled_lagrange_multiplier_ll;
 430                     im = invmass[i];
 431                     jm = invmass[j];
 432                     positions[ix] += xh * im;
 433                     positions[iy] += yh * im;
 434                     positions[iz] += zh * im;
 435                     positions[jx] -= xh * jm;
 436                     positions[jy] -= yh * jm;
 437                     positions[jz] -= zh * jm;
 438                 }
 439             }
 440         }
 441     }
 442     *nnit   = nit;
 443     *nerror = error;
 444 }
 445
 446 //! Implements RATTLE (ie. SHAKE for velocity verlet integrators)
 447 static void crattle(const int  iatom[],
 448                     int        ncon,
 449                     int*       nnit,
 450                     int        maxnit,
 451                     const real constraint_distance_squared[],
 452                     real       vp[],
 453                     const real rij[],
 454                     const real m2[],
 455                     real       omega,
 456                     const real invmass[],
 457                     const real distance_squared_tolerance[],
 458                     real       scaled_lagrange_multiplier[],
 459                     int*       nerror,
 460                     real       invdt)
 461 {
 462     /*
 463      *     r.c. van schaik and w.f. van gunsteren
 464      *     eth zuerich
 465      *     june 1992
 466      *     Adapted for use with Gromacs by David van der Spoel november 92 and later.
 467      *     rattle added by M.R. Shirts, April 2004, from code written by Jay Ponder in TINKER
 468      *     second part of rattle algorithm
 469      */
 470
 471     int  ll, i, j, i3, j3, l3;
 472     int  ix, iy, iz, jx, jy, jz;
 473     real constraint_distance_squared_ll;
 474     real vpijd, vx, vy, vz, acor, fac, im, jm;
 475     real xh, yh, zh, rijx, rijy, rijz;
 476     int  nit, error, nconv;
 477     real iconvf;
 478
 479     // TODO nconv is used solely as a boolean, so we should write the
 480     // code like that
 481     error = 0;
 482     nconv = 1;
 483     for (nit = 0; (nit < maxnit) && (nconv != 0) && (error == 0); nit++)
 484     {
 485         nconv = 0;
 486         for (ll = 0; (ll < ncon) && (error == 0); ll++)
 487         {
 488             l3   = 3 * ll;
 489             rijx = rij[l3 + XX];
 490             rijy = rij[l3 + YY];
 491             rijz = rij[l3 + ZZ];
 492             i    = iatom[l3 + 1];
 493             j    = iatom[l3 + 2];
 494             i3   = 3 * i;
 495             j3   = 3 * j;
 496             ix   = i3 + XX;
 497             iy   = i3 + YY;
 498             iz   = i3 + ZZ;
 499             jx   = j3 + XX;
 500             jy   = j3 + YY;
 501             jz   = j3 + ZZ;
 502             vx   = vp[ix] - vp[jx];
 503             vy   = vp[iy] - vp[jy];
 504             vz   = vp[iz] - vp[jz];
 505
 506             vpijd                          = vx * rijx + vy * rijy + vz * rijz;
 507             constraint_distance_squared_ll = constraint_distance_squared[ll];
 508
 509             /* iconv is zero when the error is smaller than a bound */
 510             iconvf = fabs(vpijd) * (distance_squared_tolerance[ll] / invdt);
 511
 512             if (iconvf > 1)
 513             {
 514                 nconv = static_cast<int>(iconvf);
 515                 fac   = omega * 2.0 * m2[ll] / constraint_distance_squared_ll;
 516                 acor  = -fac * vpijd;
 517                 scaled_lagrange_multiplier[ll] += acor;
 518                 xh = rijx * acor;
 519                 yh = rijy * acor;
 520                 zh = rijz * acor;
 521
 522                 im = invmass[i];
 523                 jm = invmass[j];
 524
 525                 vp[ix] += xh * im;
 526                 vp[iy] += yh * im;
 527                 vp[iz] += zh * im;
 528                 vp[jx] -= xh * jm;
 529                 vp[jy] -= yh * jm;
 530                 vp[jz] -= zh * jm;
 531             }
 532         }
 533     }
 534     *nnit   = nit;
 535     *nerror = error;
 536 }
 537
 538 //! Applies SHAKE
 539 static int vec_shakef(FILE*                     fplog,
 540                       shakedata*                shaked,
 541                       const real                invmass[],
 542                       int                       ncon,
 543                       ArrayRef<const t_iparams> ip,
 544                       const int*                iatom,
 545                       real                      tol,
 546                       const rvec                x[],
 547                       rvec                      prime[],
 548                       real                      omega,
 549                       bool                      bFEP,
 550                       real                      lambda,
 551                       real                      scaled_lagrange_multiplier[],
 552                       real                      invdt,
 553                       rvec*                     v,
 554                       bool                      bCalcVir,
 555                       tensor                    vir_r_m_dr,
 556                       ConstraintVariable        econq)
 557 {
 558     rvec* rij;
 559     real *half_of_reduced_mass, *distance_squared_tolerance, *constraint_distance_squared;
 560     int   maxnit = 1000;
 561     int   nit    = 0, ll, i, j, d, d2, type;
 562     real  L1;
 563     real  mm    = 0., tmp;
 564     int   error = 0;
 565     real  constraint_distance;
 566
 567     if (ncon > shaked->nalloc)
 568     {
 569         shaked->nalloc = over_alloc_dd(ncon);
 570         srenew(shaked->rij, shaked->nalloc);
 571         srenew(shaked->half_of_reduced_mass, shaked->nalloc);
 572         srenew(shaked->distance_squared_tolerance, shaked->nalloc);
 573         srenew(shaked->constraint_distance_squared, shaked->nalloc);
 574     }
 575     rij                         = shaked->rij;
 576     half_of_reduced_mass        = shaked->half_of_reduced_mass;
 577     distance_squared_tolerance  = shaked->distance_squared_tolerance;
 578     constraint_distance_squared = shaked->constraint_distance_squared;
 579
 580     L1            = 1.0 - lambda;
 581     const int* ia = iatom;
 582     for (ll = 0; (ll < ncon); ll++, ia += 3)
 583     {
 584         type = ia[0];
 585         i    = ia[1];
 586         j    = ia[2];
 587
 588         mm                       = 2.0 * (invmass[i] + invmass[j]);
 589         rij[ll][XX]              = x[i][XX] - x[j][XX];
 590         rij[ll][YY]              = x[i][YY] - x[j][YY];
 591         rij[ll][ZZ]              = x[i][ZZ] - x[j][ZZ];
 592         half_of_reduced_mass[ll] = 1.0 / mm;
 593         if (bFEP)
 594         {
 595             constraint_distance = L1 * ip[type].constr.dA + lambda * ip[type].constr.dB;
 596         }
 597         else
 598         {
 599             constraint_distance = ip[type].constr.dA;
 600         }
 601         constraint_distance_squared[ll] = gmx::square(constraint_distance);
 602         distance_squared_tolerance[ll]  = 0.5 / (constraint_distance_squared[ll] * tol);
 603     }
 604
 605     switch (econq)
 606     {
 607         case ConstraintVariable::Positions:
 608             cshake(iatom, ncon, &nit, maxnit, constraint_distance_squared, prime[0], rij[0],
 609                    half_of_reduced_mass, omega, invmass, distance_squared_tolerance,
 610                    scaled_lagrange_multiplier, &error);
 611             break;
 612         case ConstraintVariable::Velocities:
 613             crattle(iatom, ncon, &nit, maxnit, constraint_distance_squared, prime[0], rij[0],
 614                     half_of_reduced_mass, omega, invmass, distance_squared_tolerance,
 615                     scaled_lagrange_multiplier, &error, invdt);
 616             break;
 617         default: gmx_incons("Unknown constraint quantity for SHAKE");
 618     }
 619
 620     if (nit >= maxnit)
 621     {
 622         if (fplog)
 623         {
 624             fprintf(fplog, "Shake did not converge in %d steps\n", maxnit);
 625         }
 626         fprintf(stderr, "Shake did not converge in %d steps\n", maxnit);
 627         nit = 0;
 628     }
 629     else if (error != 0)
 630     {
 631         if (fplog)
 632         {
 633             fprintf(fplog,
 634                     "Inner product between old and new vector <= 0.0!\n"
 635                     "constraint #%d atoms %d and %d\n",
 636                     error - 1, iatom[3 * (error - 1) + 1] + 1, iatom[3 * (error - 1) + 2] + 1);
 637         }
 638         fprintf(stderr,
 639                 "Inner product between old and new vector <= 0.0!\n"
 640                 "constraint #%d atoms %d and %d\n",
 641                 error - 1, iatom[3 * (error - 1) + 1] + 1, iatom[3 * (error - 1) + 2] + 1);
 642         nit = 0;
 643     }
 644
 645     /* Constraint virial and correct the Lagrange multipliers for the length */
 646
 647     ia = iatom;
 648
 649     for (ll = 0; (ll < ncon); ll++, ia += 3)
 650     {
 651         type = ia[0];
 652         i    = ia[1];
 653         j    = ia[2];
 654
 655         if ((econq == ConstraintVariable::Positions) && v != nullptr)
 656         {
 657             /* Correct the velocities */
 658             mm = scaled_lagrange_multiplier[ll] * invmass[i] * invdt;
 659             for (d = 0; d < DIM; d++)
 660             {
 661                 v[ia[1]][d] += mm * rij[ll][d];
 662             }
 663             mm = scaled_lagrange_multiplier[ll] * invmass[j] * invdt;
 664             for (d = 0; d < DIM; d++)
 665             {
 666                 v[ia[2]][d] -= mm * rij[ll][d];
 667             }
 668             /* 16 flops */
 669         }
 670
 671         /* constraint virial */
 672         if (bCalcVir)
 673         {
 674             mm = scaled_lagrange_multiplier[ll];
 675             for (d = 0; d < DIM; d++)
 676             {
 677                 tmp = mm * rij[ll][d];
 678                 for (d2 = 0; d2 < DIM; d2++)
 679                 {
 680                     vir_r_m_dr[d][d2] -= tmp * rij[ll][d2];
 681                 }
 682             }
 683             /* 21 flops */
 684         }
 685
 686         /* cshake and crattle produce Lagrange multipliers scaled by
 687            the reciprocal of the constraint length, so fix that */
 688         if (bFEP)
 689         {
 690             constraint_distance = L1 * ip[type].constr.dA + lambda * ip[type].constr.dB;
 691         }
 692         else
 693         {
 694             constraint_distance = ip[type].constr.dA;
 695         }
 696         scaled_lagrange_multiplier[ll] *= constraint_distance;
 697     }
 698
 699     return nit;
 700 }
 701
 702 //! Check that constraints are satisfied.
 703 static void check_cons(FILE*                     log,
 704                        int                       nc,
 705                        const rvec                x[],
 706                        rvec                      prime[],
 707                        rvec                      v[],
 708                        ArrayRef<const t_iparams> ip,
 709                        const int*                iatom,
 710                        const real                invmass[],
 711                        ConstraintVariable        econq)
 712 {
 713     int  ai, aj;
 714     int  i;
 715     real d, dp;
 716     rvec dx, dv;
 717
 718     GMX_ASSERT(v, "Input has to be non-null");
 719     fprintf(log, "    i     mi      j     mj      before       after   should be\n");
 720     const int* ia = iatom;
 721     for (i = 0; (i < nc); i++, ia += 3)
 722     {
 723         ai = ia[1];
 724         aj = ia[2];
 725         rvec_sub(x[ai], x[aj], dx);
 726         d = norm(dx);
 727
 728         switch (econq)
 729         {
 730             case ConstraintVariable::Positions:
 731                 rvec_sub(prime[ai], prime[aj], dx);
 732                 dp = norm(dx);
 733                 fprintf(log, "%5d  %5.2f  %5d  %5.2f  %10.5f  %10.5f  %10.5f\n", ai + 1,
 734                         1.0 / invmass[ai], aj + 1, 1.0 / invmass[aj], d, dp, ip[ia[0]].constr.dA);
 735                 break;
 736             case ConstraintVariable::Velocities:
 737                 rvec_sub(v[ai], v[aj], dv);
 738                 d = iprod(dx, dv);
 739                 rvec_sub(prime[ai], prime[aj], dv);
 740                 dp = iprod(dx, dv);
 741                 fprintf(log, "%5d  %5.2f  %5d  %5.2f  %10.5f  %10.5f  %10.5f\n", ai + 1,
 742                         1.0 / invmass[ai], aj + 1, 1.0 / invmass[aj], d, dp, 0.);
 743                 break;
 744             default: gmx_incons("Unknown constraint quantity for SHAKE");
 745         }
 746     }
 747 }
 748
 749 //! Applies SHAKE.
 750 static bool bshakef(FILE*                         log,
 751                     shakedata*                    shaked,
 752                     const real                    invmass[],
 753                     const InteractionDefinitions& idef,
 754                     const t_inputrec&             ir,
 755                     const rvec                    x_s[],
 756                     rvec                          prime[],
 757                     t_nrnb*                       nrnb,
 758                     real                          lambda,
 759                     real*                         dvdlambda,
 760                     real                          invdt,
 761                     rvec*                         v,
 762                     bool                          bCalcVir,
 763                     tensor                        vir_r_m_dr,
 764                     bool                          bDumpOnError,
 765                     ConstraintVariable            econq)
 766 {
 767     real *lam, dt_2, dvdl;
 768     int   i, n0, ncon, blen, type, ll;
 769     int   tnit = 0, trij = 0;
 770
 771     ncon = idef.il[F_CONSTR].size() / 3;
 772
 773     for (ll = 0; ll < ncon; ll++)
 774     {
 775         shaked->scaled_lagrange_multiplier[ll] = 0;
 776     }
 777
 778     // TODO Rewrite this block so that it is obvious that i, iatoms
 779     // and lam are all iteration variables. Is this easier if the
 780     // sblock data structure is organized differently?
 781     const int* iatoms = &(idef.il[F_CONSTR].iatoms[shaked->sblock[0]]);
 782     lam               = shaked->scaled_lagrange_multiplier;
 783     for (i = 0; (i < shaked->nblocks);)
 784     {
 785         blen = (shaked->sblock[i + 1] - shaked->sblock[i]);
 786         blen /= 3;
 787         n0 = vec_shakef(log, shaked, invmass, blen, idef.iparams, iatoms, ir.shake_tol, x_s, prime,
 788                         shaked->omega, ir.efep != efepNO, lambda, lam, invdt, v, bCalcVir,
 789                         vir_r_m_dr, econq);
 790
 791         if (n0 == 0)
 792         {
 793             if (bDumpOnError && log)
 794             {
 795                 {
 796                     check_cons(log, blen, x_s, prime, v, idef.iparams, iatoms, invmass, econq);
 797                 }
 798             }
 799             return FALSE;
 800         }
 801         tnit += n0 * blen;
 802         trij += blen;
 803         iatoms += 3 * blen; /* Increment pointer! */
 804         lam += blen;
 805         i++;
 806     }
 807     /* only for position part? */
 808     if (econq == ConstraintVariable::Positions)
 809     {
 810         if (ir.efep != efepNO)
 811         {
 812             ArrayRef<const t_iparams> iparams = idef.iparams;
 813
 814             /* TODO This should probably use invdt, so that sd integrator scaling works properly */
 815             dt_2 = 1 / gmx::square(ir.delta_t);
 816             dvdl = 0;
 817             for (ll = 0; ll < ncon; ll++)
 818             {
 819                 type = idef.il[F_CONSTR].iatoms[3 * ll];
 820
 821                 /* Per equations in the manual, dv/dl = -2 \sum_ll lagrangian_ll * r_ll * (d_B - d_A) */
 822                 /* The vector scaled_lagrange_multiplier[ll] contains the value -2 r_ll eta_ll
 823                    (eta_ll is the estimate of the Langrangian, definition on page 336 of Ryckaert et
 824                    al 1977), so the pre-factors are already present. */
 825                 const real bondA = iparams[type].constr.dA;
 826                 const real bondB = iparams[type].constr.dB;
 827                 dvdl += shaked->scaled_lagrange_multiplier[ll] * dt_2 * (bondB - bondA);
 828             }
 829             *dvdlambda += dvdl;
 830         }
 831     }
 832     if (ir.bShakeSOR)
 833     {
 834         if (tnit > shaked->gamma)
 835         {
 836             shaked->delta *= -0.5;
 837         }
 838         shaked->omega += shaked->delta;
 839         shaked->gamma = tnit;
 840     }
 841     inc_nrnb(nrnb, eNR_SHAKE, tnit);
 842     inc_nrnb(nrnb, eNR_SHAKE_RIJ, trij);
 843     if (v)
 844     {
 845         inc_nrnb(nrnb, eNR_CONSTR_V, trij * 2);
 846     }
 847     if (bCalcVir)
 848     {
 849         inc_nrnb(nrnb, eNR_CONSTR_VIR, trij);
 850     }
 851
 852     return TRUE;
 853 }
 854
 855 bool constrain_shake(FILE*                         log,
 856                      shakedata*                    shaked,
 857                      const real                    invmass[],
 858                      const InteractionDefinitions& idef,
 859                      const t_inputrec&             ir,
 860                      const rvec                    x_s[],
 861                      rvec                          xprime[],
 862                      rvec                          vprime[],
 863                      t_nrnb*                       nrnb,
 864                      real                          lambda,
 865                      real*                         dvdlambda,
 866                      real                          invdt,
 867                      rvec*                         v,
 868                      bool                          bCalcVir,
 869                      tensor                        vir_r_m_dr,
 870                      bool                          bDumpOnError,
 871                      ConstraintVariable            econq)
 872 {
 873     if (shaked->nblocks == 0)
 874     {
 875         return true;
 876     }
 877     bool bOK;
 878     switch (econq)
 879     {
 880         case (ConstraintVariable::Positions):
 881             bOK = bshakef(log, shaked, invmass, idef, ir, x_s, xprime, nrnb, lambda, dvdlambda,
 882                           invdt, v, bCalcVir, vir_r_m_dr, bDumpOnError, econq);
 883             break;
 884         case (ConstraintVariable::Velocities):
 885             bOK = bshakef(log, shaked, invmass, idef, ir, x_s, vprime, nrnb, lambda, dvdlambda,
 886                           invdt, nullptr, bCalcVir, vir_r_m_dr, bDumpOnError, econq);
 887             break;
 888         default:
 889             gmx_fatal(FARGS,
 890                       "Internal error, SHAKE called for constraining something else than "
 891                       "coordinates");
 892     }
 893     return bOK;
 894 }
 895
 896 } // namespace gmx