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