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