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