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