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, const t_idef* 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     t_iatom*     iatom;
 189     int*         inv_sblock;
 190
 191     /* Since we are processing the local topology,
 192      * the F_CONSTRNC ilist has been concatenated to the F_CONSTR ilist.
 193      */
 194     ncons = idef->il[F_CONSTR].nr / 3;
 195
 196     init_blocka(&sblocks);
 197     sfree(sblocks.index); // To solve memory leak
 198     gen_sblocks(nullptr, 0, md.homenr, idef, &sblocks, FALSE);
 199
 200     /*
 201        bstart=(idef->nodeid > 0) ? blocks->multinr[idef->nodeid-1] : 0;
 202        nblocks=blocks->multinr[idef->nodeid] - bstart;
 203      */
 204     bstart          = 0;
 205     shaked->nblocks = sblocks.nr;
 206     if (debug)
 207     {
 208         fprintf(debug, "ncons: %d, bstart: %d, nblocks: %d\n", ncons, bstart, shaked->nblocks);
 209     }
 210
 211     /* Calculate block number for each atom */
 212     inv_sblock = make_invblocka(&sblocks, md.nr);
 213
 214     done_blocka(&sblocks);
 215
 216     /* Store the block number in temp array and
 217      * sort the constraints in order of the sblock number
 218      * and the atom numbers, really sorting a segment of the array!
 219      */
 220     iatom = idef->il[F_CONSTR].iatoms;
 221     snew(sb, ncons);
 222     for (i = 0; (i < ncons); i++, iatom += 3)
 223     {
 224         for (m = 0; (m < 3); m++)
 225         {
 226             sb[i].iatom[m] = iatom[m];
 227         }
 228         sb[i].blocknr = inv_sblock[iatom[1]];
 229     }
 230
 231     /* Now sort the blocks */
 232     if (debug)
 233     {
 234         pr_sortblock(debug, "Before sorting", ncons, sb);
 235         fprintf(debug, "Going to sort constraints\n");
 236     }
 237
 238     std::qsort(sb, ncons, sizeof(*sb), pcomp);
 239
 240     if (debug)
 241     {
 242         pr_sortblock(debug, "After sorting", ncons, sb);
 243     }
 244
 245     iatom = idef->il[F_CONSTR].iatoms;
 246     for (i = 0; (i < ncons); i++, iatom += 3)
 247     {
 248         for (m = 0; (m < 3); m++)
 249         {
 250             iatom[m] = sb[i].iatom[m];
 251         }
 252     }
 253
 254     j = 0;
 255     snew(shaked->sblock, shaked->nblocks + 1);
 256     bnr = -2;
 257     for (i = 0; (i < ncons); i++)
 258     {
 259         if (sb[i].blocknr != bnr)
 260         {
 261             bnr                 = sb[i].blocknr;
 262             shaked->sblock[j++] = 3 * i;
 263         }
 264     }
 265     /* Last block... */
 266     shaked->sblock[j++] = 3 * ncons;
 267
 268     if (j != (shaked->nblocks + 1))
 269     {
 270         fprintf(stderr, "bstart: %d\n", bstart);
 271         fprintf(stderr, "j: %d, nblocks: %d, ncons: %d\n", j, shaked->nblocks, ncons);
 272         for (i = 0; (i < ncons); i++)
 273         {
 274             fprintf(stderr, "i: %5d  sb[i].blocknr: %5d\n", i, sb[i].blocknr);
 275         }
 276         for (j = 0; (j <= shaked->nblocks); j++)
 277         {
 278             fprintf(stderr, "sblock[%3d]=%5d\n", j, shaked->sblock[j]);
 279         }
 280         gmx_fatal(FARGS,
 281                   "DEATH HORROR: "
 282                   "sblocks does not match idef->il[F_CONSTR]");
 283     }
 284     sfree(sb);
 285     sfree(inv_sblock);
 286     resizeLagrangianData(shaked, ncons);
 287 }
 288
 289 // TODO: Check if this code is useful. It might never be called.
 290 void make_shake_sblock_dd(shakedata* shaked, const t_ilist* ilcon, const gmx_domdec_t* dd)
 291 {
 292     int      ncons, c, cg;
 293     t_iatom* iatom;
 294
 295     if (dd->ncg_home + 1 > shaked->sblock_nalloc)
 296     {
 297         shaked->sblock_nalloc = over_alloc_dd(dd->ncg_home + 1);
 298         srenew(shaked->sblock, shaked->sblock_nalloc);
 299     }
 300
 301     ncons           = ilcon->nr / 3;
 302     iatom           = ilcon->iatoms;
 303     shaked->nblocks = 0;
 304     cg              = 0;
 305     for (c = 0; c < ncons; c++)
 306     {
 307         if (c == 0 || iatom[1] >= cg + 1)
 308         {
 309             shaked->sblock[shaked->nblocks++] = 3 * c;
 310             while (iatom[1] >= cg + 1)
 311             {
 312                 cg++;
 313             }
 314         }
 315         iatom += 3;
 316     }
 317     shaked->sblock[shaked->nblocks] = 3 * ncons;
 318     resizeLagrangianData(shaked, ncons);
 319 }
 320
 321 /*! \brief Inner kernel for SHAKE constraints
 322  *
 323  * Original implementation from R.C. van Schaik and W.F. van Gunsteren
 324  * (ETH Zuerich, June 1992), adapted for GROMACS by David van der
 325  * Spoel November 1992.
 326  *
 327  * The algorithm here is based section five of Ryckaert, Ciccotti and
 328  * Berendsen, J Comp Phys, 23, 327, 1977.
 329  *
 330  * \param[in]    iatom                         Mini-topology of triples of constraint type (unused in this
 331  *                                             function) and indices of the two atoms involved
 332  * \param[in]    ncon                          Number of constraints
 333  * \param[out]   nnit                          Number of iterations performed
 334  * \param[in]    maxnit                        Maximum number of iterations permitted
 335  * \param[in]    constraint_distance_squared   The objective value for each constraint
 336  * \param[inout] positions                     The initial (and final) values of the positions of all atoms
 337  * \param[in]    initial_displacements         The initial displacements of each constraint
 338  * \param[in]    half_of_reduced_mass          Half of the reduced mass for each constraint
 339  * \param[in]    omega                         SHAKE over-relaxation factor (set non-1.0 by
 340  *                                             using shake-sor=yes in the .mdp, but there is no documentation anywhere)
 341  * \param[in]    invmass                       Inverse mass of each atom
 342  * \param[in]    distance_squared_tolerance    Multiplicative tolerance on the difference in the
 343  *                                             square of the constrained distance (see code)
 344  * \param[out]   scaled_lagrange_multiplier    Scaled Lagrange multiplier for each constraint (-2 * eta from p. 336
 345  *                                             of the paper, divided by the constraint distance)
 346  * \param[out]   nerror                        Zero upon success, returns one more than the index of the
 347  *                                             problematic constraint if the input was malformed
 348  *
 349  * \todo Make SHAKE use better data structures, in particular for iatom. */
 350 void cshake(const int  iatom[],
 351             int        ncon,
 352             int*       nnit,
 353             int        maxnit,
 354             const real constraint_distance_squared[],
 355             real       positions[],
 356             const real initial_displacements[],
 357             const real half_of_reduced_mass[],
 358             real       omega,
 359             const real invmass[],
 360             const real distance_squared_tolerance[],
 361             real       scaled_lagrange_multiplier[],
 362             int*       nerror)
 363 {
 364     /* default should be increased! MRS 8/4/2009 */
 365     const real mytol = 1e-10;
 366
 367     int  ll, i, j, i3, j3, l3;
 368     int  ix, iy, iz, jx, jy, jz;
 369     real r_dot_r_prime;
 370     real constraint_distance_squared_ll;
 371     real r_prime_squared;
 372     real scaled_lagrange_multiplier_ll;
 373     real r_prime_x, r_prime_y, r_prime_z, diff, im, jm;
 374     real xh, yh, zh, rijx, rijy, rijz;
 375     int  nit, error, nconv;
 376     real iconvf;
 377
 378     // TODO nconv is used solely as a boolean, so we should write the
 379     // code like that
 380     error = 0;
 381     nconv = 1;
 382     for (nit = 0; (nit < maxnit) && (nconv != 0) && (error == 0); nit++)
 383     {
 384         nconv = 0;
 385         for (ll = 0; (ll < ncon) && (error == 0); ll++)
 386         {
 387             l3   = 3 * ll;
 388             rijx = initial_displacements[l3 + XX];
 389             rijy = initial_displacements[l3 + YY];
 390             rijz = initial_displacements[l3 + ZZ];
 391             i    = iatom[l3 + 1];
 392             j    = iatom[l3 + 2];
 393             i3   = 3 * i;
 394             j3   = 3 * j;
 395             ix   = i3 + XX;
 396             iy   = i3 + YY;
 397             iz   = i3 + ZZ;
 398             jx   = j3 + XX;
 399             jy   = j3 + YY;
 400             jz   = j3 + ZZ;
 401
 402             /* Compute r prime between atoms i and j, which is the
 403                displacement *before* this update stage */
 404             r_prime_x       = positions[ix] - positions[jx];
 405             r_prime_y       = positions[iy] - positions[jy];
 406             r_prime_z       = positions[iz] - positions[jz];
 407             r_prime_squared = (r_prime_x * r_prime_x + r_prime_y * r_prime_y + r_prime_z * r_prime_z);
 408             constraint_distance_squared_ll = constraint_distance_squared[ll];
 409             diff                           = constraint_distance_squared_ll - r_prime_squared;
 410
 411             /* iconvf is less than 1 when the error is smaller than a bound */
 412             iconvf = fabs(diff) * distance_squared_tolerance[ll];
 413
 414             if (iconvf > 1.0)
 415             {
 416                 nconv         = static_cast<int>(iconvf);
 417                 r_dot_r_prime = (rijx * r_prime_x + rijy * r_prime_y + rijz * r_prime_z);
 418
 419                 if (r_dot_r_prime < constraint_distance_squared_ll * mytol)
 420                 {
 421                     error = ll + 1;
 422                 }
 423                 else
 424                 {
 425                     /* The next line solves equation 5.6 (neglecting
 426                        the term in g^2), for g */
 427                     scaled_lagrange_multiplier_ll = omega * diff * half_of_reduced_mass[ll] / r_dot_r_prime;
 428                     scaled_lagrange_multiplier[ll] += scaled_lagrange_multiplier_ll;
 429                     xh = rijx * scaled_lagrange_multiplier_ll;
 430                     yh = rijy * scaled_lagrange_multiplier_ll;
 431                     zh = rijz * scaled_lagrange_multiplier_ll;
 432                     im = invmass[i];
 433                     jm = invmass[j];
 434                     positions[ix] += xh * im;
 435                     positions[iy] += yh * im;
 436                     positions[iz] += zh * im;
 437                     positions[jx] -= xh * jm;
 438                     positions[jy] -= yh * jm;
 439                     positions[jz] -= zh * jm;
 440                 }
 441             }
 442         }
 443     }
 444     *nnit   = nit;
 445     *nerror = error;
 446 }
 447
 448 //! Implements RATTLE (ie. SHAKE for velocity verlet integrators)
 449 static void crattle(const int  iatom[],
 450                     int        ncon,
 451                     int*       nnit,
 452                     int        maxnit,
 453                     const real constraint_distance_squared[],
 454                     real       vp[],
 455                     const real rij[],
 456                     const real m2[],
 457                     real       omega,
 458                     const real invmass[],
 459                     const real distance_squared_tolerance[],
 460                     real       scaled_lagrange_multiplier[],
 461                     int*       nerror,
 462                     real       invdt)
 463 {
 464     /*
 465      *     r.c. van schaik and w.f. van gunsteren
 466      *     eth zuerich
 467      *     june 1992
 468      *     Adapted for use with Gromacs by David van der Spoel november 92 and later.
 469      *     rattle added by M.R. Shirts, April 2004, from code written by Jay Ponder in TINKER
 470      *     second part of rattle algorithm
 471      */
 472
 473     int  ll, i, j, i3, j3, l3;
 474     int  ix, iy, iz, jx, jy, jz;
 475     real constraint_distance_squared_ll;
 476     real vpijd, vx, vy, vz, acor, fac, im, jm;
 477     real xh, yh, zh, rijx, rijy, rijz;
 478     int  nit, error, nconv;
 479     real iconvf;
 480
 481     // TODO nconv is used solely as a boolean, so we should write the
 482     // code like that
 483     error = 0;
 484     nconv = 1;
 485     for (nit = 0; (nit < maxnit) && (nconv != 0) && (error == 0); nit++)
 486     {
 487         nconv = 0;
 488         for (ll = 0; (ll < ncon) && (error == 0); ll++)
 489         {
 490             l3   = 3 * ll;
 491             rijx = rij[l3 + XX];
 492             rijy = rij[l3 + YY];
 493             rijz = rij[l3 + ZZ];
 494             i    = iatom[l3 + 1];
 495             j    = iatom[l3 + 2];
 496             i3   = 3 * i;
 497             j3   = 3 * j;
 498             ix   = i3 + XX;
 499             iy   = i3 + YY;
 500             iz   = i3 + ZZ;
 501             jx   = j3 + XX;
 502             jy   = j3 + YY;
 503             jz   = j3 + ZZ;
 504             vx   = vp[ix] - vp[jx];
 505             vy   = vp[iy] - vp[jy];
 506             vz   = vp[iz] - vp[jz];
 507
 508             vpijd                          = vx * rijx + vy * rijy + vz * rijz;
 509             constraint_distance_squared_ll = constraint_distance_squared[ll];
 510
 511             /* iconv is zero when the error is smaller than a bound */
 512             iconvf = fabs(vpijd) * (distance_squared_tolerance[ll] / invdt);
 513
 514             if (iconvf > 1)
 515             {
 516                 nconv = static_cast<int>(iconvf);
 517                 fac   = omega * 2.0 * m2[ll] / constraint_distance_squared_ll;
 518                 acor  = -fac * vpijd;
 519                 scaled_lagrange_multiplier[ll] += acor;
 520                 xh = rijx * acor;
 521                 yh = rijy * acor;
 522                 zh = rijz * acor;
 523
 524                 im = invmass[i];
 525                 jm = invmass[j];
 526
 527                 vp[ix] += xh * im;
 528                 vp[iy] += yh * im;
 529                 vp[iz] += zh * im;
 530                 vp[jx] -= xh * jm;
 531                 vp[jy] -= yh * jm;
 532                 vp[jz] -= zh * jm;
 533             }
 534         }
 535     }
 536     *nnit   = nit;
 537     *nerror = error;
 538 }
 539
 540 //! Applies SHAKE
 541 static int vec_shakef(FILE*              fplog,
 542                       shakedata*         shaked,
 543                       const real         invmass[],
 544                       int                ncon,
 545                       t_iparams          ip[],
 546                       t_iatom*           iatom,
 547                       real               tol,
 548                       const rvec         x[],
 549                       rvec               prime[],
 550                       real               omega,
 551                       bool               bFEP,
 552                       real               lambda,
 553                       real               scaled_lagrange_multiplier[],
 554                       real               invdt,
 555                       rvec*              v,
 556                       bool               bCalcVir,
 557                       tensor             vir_r_m_dr,
 558                       ConstraintVariable econq)
 559 {
 560     rvec*    rij;
 561     real *   half_of_reduced_mass, *distance_squared_tolerance, *constraint_distance_squared;
 562     int      maxnit = 1000;
 563     int      nit    = 0, ll, i, j, d, d2, type;
 564     t_iatom* ia;
 565     real     L1;
 566     real     mm    = 0., tmp;
 567     int      error = 0;
 568     real     constraint_distance;
 569
 570     if (ncon > shaked->nalloc)
 571     {
 572         shaked->nalloc = over_alloc_dd(ncon);
 573         srenew(shaked->rij, shaked->nalloc);
 574         srenew(shaked->half_of_reduced_mass, shaked->nalloc);
 575         srenew(shaked->distance_squared_tolerance, shaked->nalloc);
 576         srenew(shaked->constraint_distance_squared, shaked->nalloc);
 577     }
 578     rij                         = shaked->rij;
 579     half_of_reduced_mass        = shaked->half_of_reduced_mass;
 580     distance_squared_tolerance  = shaked->distance_squared_tolerance;
 581     constraint_distance_squared = shaked->constraint_distance_squared;
 582
 583     L1 = 1.0 - lambda;
 584     ia = iatom;
 585     for (ll = 0; (ll < ncon); ll++, ia += 3)
 586     {
 587         type = ia[0];
 588         i    = ia[1];
 589         j    = ia[2];
 590
 591         mm                       = 2.0 * (invmass[i] + invmass[j]);
 592         rij[ll][XX]              = x[i][XX] - x[j][XX];
 593         rij[ll][YY]              = x[i][YY] - x[j][YY];
 594         rij[ll][ZZ]              = x[i][ZZ] - x[j][ZZ];
 595         half_of_reduced_mass[ll] = 1.0 / mm;
 596         if (bFEP)
 597         {
 598             constraint_distance = L1 * ip[type].constr.dA + lambda * ip[type].constr.dB;
 599         }
 600         else
 601         {
 602             constraint_distance = ip[type].constr.dA;
 603         }
 604         constraint_distance_squared[ll] = gmx::square(constraint_distance);
 605         distance_squared_tolerance[ll]  = 0.5 / (constraint_distance_squared[ll] * tol);
 606     }
 607
 608     switch (econq)
 609     {
 610         case ConstraintVariable::Positions:
 611             cshake(iatom, ncon, &nit, maxnit, constraint_distance_squared, prime[0], rij[0],
 612                    half_of_reduced_mass, omega, invmass, distance_squared_tolerance,
 613                    scaled_lagrange_multiplier, &error);
 614             break;
 615         case ConstraintVariable::Velocities:
 616             crattle(iatom, ncon, &nit, maxnit, constraint_distance_squared, prime[0], rij[0],
 617                     half_of_reduced_mass, omega, invmass, distance_squared_tolerance,
 618                     scaled_lagrange_multiplier, &error, invdt);
 619             break;
 620         default: gmx_incons("Unknown constraint quantity for SHAKE");
 621     }
 622
 623     if (nit >= maxnit)
 624     {
 625         if (fplog)
 626         {
 627             fprintf(fplog, "Shake did not converge in %d steps\n", maxnit);
 628         }
 629         fprintf(stderr, "Shake did not converge in %d steps\n", maxnit);
 630         nit = 0;
 631     }
 632     else if (error != 0)
 633     {
 634         if (fplog)
 635         {
 636             fprintf(fplog,
 637                     "Inner product between old and new vector <= 0.0!\n"
 638                     "constraint #%d atoms %d and %d\n",
 639                     error - 1, iatom[3 * (error - 1) + 1] + 1, iatom[3 * (error - 1) + 2] + 1);
 640         }
 641         fprintf(stderr,
 642                 "Inner product between old and new vector <= 0.0!\n"
 643                 "constraint #%d atoms %d and %d\n",
 644                 error - 1, iatom[3 * (error - 1) + 1] + 1, iatom[3 * (error - 1) + 2] + 1);
 645         nit = 0;
 646     }
 647
 648     /* Constraint virial and correct the Lagrange multipliers for the length */
 649
 650     ia = iatom;
 651
 652     for (ll = 0; (ll < ncon); ll++, ia += 3)
 653     {
 654         type = ia[0];
 655         i    = ia[1];
 656         j    = ia[2];
 657
 658         if ((econq == ConstraintVariable::Positions) && v != nullptr)
 659         {
 660             /* Correct the velocities */
 661             mm = scaled_lagrange_multiplier[ll] * invmass[i] * invdt;
 662             for (d = 0; d < DIM; d++)
 663             {
 664                 v[ia[1]][d] += mm * rij[ll][d];
 665             }
 666             mm = scaled_lagrange_multiplier[ll] * invmass[j] * invdt;
 667             for (d = 0; d < DIM; d++)
 668             {
 669                 v[ia[2]][d] -= mm * rij[ll][d];
 670             }
 671             /* 16 flops */
 672         }
 673
 674         /* constraint virial */
 675         if (bCalcVir)
 676         {
 677             mm = scaled_lagrange_multiplier[ll];
 678             for (d = 0; d < DIM; d++)
 679             {
 680                 tmp = mm * rij[ll][d];
 681                 for (d2 = 0; d2 < DIM; d2++)
 682                 {
 683                     vir_r_m_dr[d][d2] -= tmp * rij[ll][d2];
 684                 }
 685             }
 686             /* 21 flops */
 687         }
 688
 689         /* cshake and crattle produce Lagrange multipliers scaled by
 690            the reciprocal of the constraint length, so fix that */
 691         if (bFEP)
 692         {
 693             constraint_distance = L1 * ip[type].constr.dA + lambda * ip[type].constr.dB;
 694         }
 695         else
 696         {
 697             constraint_distance = ip[type].constr.dA;
 698         }
 699         scaled_lagrange_multiplier[ll] *= constraint_distance;
 700     }
 701
 702     return nit;
 703 }
 704
 705 //! Check that constraints are satisfied.
 706 static void check_cons(FILE*              log,
 707                        int                nc,
 708                        const rvec         x[],
 709                        rvec               prime[],
 710                        rvec               v[],
 711                        t_iparams          ip[],
 712                        t_iatom*           iatom,
 713                        const real         invmass[],
 714                        ConstraintVariable econq)
 715 {
 716     t_iatom* ia;
 717     int      ai, aj;
 718     int      i;
 719     real     d, dp;
 720     rvec     dx, dv;
 721
 722     GMX_ASSERT(v, "Input has to be non-null");
 723     fprintf(log, "    i     mi      j     mj      before       after   should be\n");
 724     ia = iatom;
 725     for (i = 0; (i < nc); i++, ia += 3)
 726     {
 727         ai = ia[1];
 728         aj = ia[2];
 729         rvec_sub(x[ai], x[aj], dx);
 730         d = norm(dx);
 731
 732         switch (econq)
 733         {
 734             case ConstraintVariable::Positions:
 735                 rvec_sub(prime[ai], prime[aj], dx);
 736                 dp = norm(dx);
 737                 fprintf(log, "%5d  %5.2f  %5d  %5.2f  %10.5f  %10.5f  %10.5f\n", ai + 1,
 738                         1.0 / invmass[ai], aj + 1, 1.0 / invmass[aj], d, dp, ip[ia[0]].constr.dA);
 739                 break;
 740             case ConstraintVariable::Velocities:
 741                 rvec_sub(v[ai], v[aj], dv);
 742                 d = iprod(dx, dv);
 743                 rvec_sub(prime[ai], prime[aj], dv);
 744                 dp = iprod(dx, dv);
 745                 fprintf(log, "%5d  %5.2f  %5d  %5.2f  %10.5f  %10.5f  %10.5f\n", ai + 1,
 746                         1.0 / invmass[ai], aj + 1, 1.0 / invmass[aj], d, dp, 0.);
 747                 break;
 748             default: gmx_incons("Unknown constraint quantity for SHAKE");
 749         }
 750     }
 751 }
 752
 753 //! Applies SHAKE.
 754 static bool bshakef(FILE*              log,
 755                     shakedata*         shaked,
 756                     const real         invmass[],
 757                     const t_idef&      idef,
 758                     const t_inputrec&  ir,
 759                     const rvec         x_s[],
 760                     rvec               prime[],
 761                     t_nrnb*            nrnb,
 762                     real               lambda,
 763                     real*              dvdlambda,
 764                     real               invdt,
 765                     rvec*              v,
 766                     bool               bCalcVir,
 767                     tensor             vir_r_m_dr,
 768                     bool               bDumpOnError,
 769                     ConstraintVariable econq)
 770 {
 771     t_iatom* iatoms;
 772     real *   lam, dt_2, dvdl;
 773     int      i, n0, ncon, blen, type, ll;
 774     int      tnit = 0, trij = 0;
 775
 776     ncon = idef.il[F_CONSTR].nr / 3;
 777
 778     for (ll = 0; ll < ncon; ll++)
 779     {
 780         shaked->scaled_lagrange_multiplier[ll] = 0;
 781     }
 782
 783     // TODO Rewrite this block so that it is obvious that i, iatoms
 784     // and lam are all iteration variables. Is this easier if the
 785     // sblock data structure is organized differently?
 786     iatoms = &(idef.il[F_CONSTR].iatoms[shaked->sblock[0]]);
 787     lam    = shaked->scaled_lagrange_multiplier;
 788     for (i = 0; (i < shaked->nblocks);)
 789     {
 790         blen = (shaked->sblock[i + 1] - shaked->sblock[i]);
 791         blen /= 3;
 792         n0 = vec_shakef(log, shaked, invmass, blen, idef.iparams, iatoms, ir.shake_tol, x_s, prime,
 793                         shaked->omega, ir.efep != efepNO, lambda, lam, invdt, v, bCalcVir,
 794                         vir_r_m_dr, econq);
 795
 796         if (n0 == 0)
 797         {
 798             if (bDumpOnError && log)
 799             {
 800                 {
 801                     check_cons(log, blen, x_s, prime, v, idef.iparams, iatoms, invmass, econq);
 802                 }
 803             }
 804             return FALSE;
 805         }
 806         tnit += n0 * blen;
 807         trij += blen;
 808         iatoms += 3 * blen; /* Increment pointer! */
 809         lam += blen;
 810         i++;
 811     }
 812     /* only for position part? */
 813     if (econq == ConstraintVariable::Positions)
 814     {
 815         if (ir.efep != efepNO)
 816         {
 817             real bondA, bondB;
 818             /* TODO This should probably use invdt, so that sd integrator scaling works properly */
 819             dt_2 = 1 / gmx::square(ir.delta_t);
 820             dvdl = 0;
 821             for (ll = 0; ll < ncon; ll++)
 822             {
 823                 type = idef.il[F_CONSTR].iatoms[3 * ll];
 824
 825                 /* Per equations in the manual, dv/dl = -2 \sum_ll lagrangian_ll * r_ll * (d_B - d_A) */
 826                 /* The vector scaled_lagrange_multiplier[ll] contains the value -2 r_ll eta_ll
 827                    (eta_ll is the estimate of the Langrangian, definition on page 336 of Ryckaert et
 828                    al 1977), so the pre-factors are already present. */
 829                 bondA = idef.iparams[type].constr.dA;
 830                 bondB = idef.iparams[type].constr.dB;
 831                 dvdl += shaked->scaled_lagrange_multiplier[ll] * dt_2 * (bondB - bondA);
 832             }
 833             *dvdlambda += dvdl;
 834         }
 835     }
 836     if (ir.bShakeSOR)
 837     {
 838         if (tnit > shaked->gamma)
 839         {
 840             shaked->delta *= -0.5;
 841         }
 842         shaked->omega += shaked->delta;
 843         shaked->gamma = tnit;
 844     }
 845     inc_nrnb(nrnb, eNR_SHAKE, tnit);
 846     inc_nrnb(nrnb, eNR_SHAKE_RIJ, trij);
 847     if (v)
 848     {
 849         inc_nrnb(nrnb, eNR_CONSTR_V, trij * 2);
 850     }
 851     if (bCalcVir)
 852     {
 853         inc_nrnb(nrnb, eNR_CONSTR_VIR, trij);
 854     }
 855
 856     return TRUE;
 857 }
 858
 859 bool constrain_shake(FILE*              log,
 860                      shakedata*         shaked,
 861                      const real         invmass[],
 862                      const t_idef&      idef,
 863                      const t_inputrec&  ir,
 864                      const rvec         x_s[],
 865                      rvec               xprime[],
 866                      rvec               vprime[],
 867                      t_nrnb*            nrnb,
 868                      real               lambda,
 869                      real*              dvdlambda,
 870                      real               invdt,
 871                      rvec*              v,
 872                      bool               bCalcVir,
 873                      tensor             vir_r_m_dr,
 874                      bool               bDumpOnError,
 875                      ConstraintVariable econq)
 876 {
 877     if (shaked->nblocks == 0)
 878     {
 879         return true;
 880     }
 881     bool bOK;
 882     switch (econq)
 883     {
 884         case (ConstraintVariable::Positions):
 885             bOK = bshakef(log, shaked, invmass, idef, ir, x_s, xprime, nrnb, lambda, dvdlambda,
 886                           invdt, v, bCalcVir, vir_r_m_dr, bDumpOnError, econq);
 887             break;
 888         case (ConstraintVariable::Velocities):
 889             bOK = bshakef(log, shaked, invmass, idef, ir, x_s, vprime, nrnb, lambda, dvdlambda,
 890                           invdt, nullptr, bCalcVir, vir_r_m_dr, bDumpOnError, econq);
 891             break;
 892         default:
 893             gmx_fatal(FARGS,
 894                       "Internal error, SHAKE called for constraining something else than "
 895                       "coordinates");
 896     }
 897     return bOK;
 898 }
 899
 900 } // namespace gmx