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