src/gromacs/mdlib/constr.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,2016,2017,2018,2019, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 /*! \internal \file
  38  * \brief Defines the high-level constraint code.
  39  *
  40  * \author Berk Hess <hess@kth.se>
  41  * \author Mark Abraham <mark.j.abraham@gmail.com>
  42  * \ingroup module_mdlib
  43  */
  44 #include "gmxpre.h"
  45
  46 #include "constr.h"
  47
  48 #include <cassert>
  49 #include <cmath>
  50 #include <cstdlib>
  51
  52 #include <algorithm>
  53
  54 #include "gromacs/domdec/domdec.h"
  55 #include "gromacs/domdec/domdec_struct.h"
  56 #include "gromacs/essentialdynamics/edsam.h"
  57 #include "gromacs/fileio/confio.h"
  58 #include "gromacs/fileio/gmxfio.h"
  59 #include "gromacs/fileio/pdbio.h"
  60 #include "gromacs/gmxlib/nrnb.h"
  61 #include "gromacs/math/utilities.h"
  62 #include "gromacs/math/vec.h"
  63 #include "gromacs/mdlib/gmx_omp_nthreads.h"
  64 #include "gromacs/mdlib/lincs.h"
  65 #include "gromacs/mdlib/settle.h"
  66 #include "gromacs/mdlib/shake.h"
  67 #include "gromacs/mdtypes/commrec.h"
  68 #include "gromacs/mdtypes/inputrec.h"
  69 #include "gromacs/mdtypes/md_enums.h"
  70 #include "gromacs/mdtypes/mdatom.h"
  71 #include "gromacs/mdtypes/state.h"
  72 #include "gromacs/pbcutil/pbc.h"
  73 #include "gromacs/pulling/pull.h"
  74 #include "gromacs/timing/wallcycle.h"
  75 #include "gromacs/topology/block.h"
  76 #include "gromacs/topology/ifunc.h"
  77 #include "gromacs/topology/mtop_lookup.h"
  78 #include "gromacs/topology/mtop_util.h"
  79 #include "gromacs/utility/arrayref.h"
  80 #include "gromacs/utility/exceptions.h"
  81 #include "gromacs/utility/fatalerror.h"
  82 #include "gromacs/utility/gmxassert.h"
  83 #include "gromacs/utility/pleasecite.h"
  84 #include "gromacs/utility/smalloc.h"
  85 #include "gromacs/utility/txtdump.h"
  86
  87 namespace gmx
  88 {
  89
  90 /* \brief Impl class for Constraints
  91  *
  92  * \todo Members like md, idef are valid only for the lifetime of a
  93  * domain, which would be good to make clearer in the structure of the
  94  * code. It should not be possible to call apply() if setConstraints()
  95  * has not been called. For example, this could be achieved if
  96  * setConstraints returned a valid object with such a method.  That
  97  * still requires that we manage the lifetime of that object
  98  * correctly, however. */
  99 class Constraints::Impl
 100 {
 101     public:
 102         Impl(const gmx_mtop_t     &mtop_p,
 103              const t_inputrec     &ir_p,
 104              pull_t               *pull_work,
 105              FILE                 *log_p,
 106              const t_mdatoms      &md_p,
 107              const t_commrec      *cr_p,
 108              const gmx_multisim_t *ms,
 109              t_nrnb               *nrnb,
 110              gmx_wallcycle        *wcycle_p,
 111              bool                  pbcHandlingRequired,
 112              int                   numConstraints,
 113              int                   numSettles);
 114         ~Impl();
 115         void setConstraints(const gmx_localtop_t &top,
 116                             const t_mdatoms      &md);
 117         bool apply(bool                  bLog,
 118                    bool                  bEner,
 119                    int64_t               step,
 120                    int                   delta_step,
 121                    real                  step_scaling,
 122                    rvec                 *x,
 123                    rvec                 *xprime,
 124                    rvec                 *min_proj,
 125                    matrix                box,
 126                    real                  lambda,
 127                    real                 *dvdlambda,
 128                    rvec                 *v,
 129                    tensor               *vir,
 130                    ConstraintVariable    econq);
 131         //! The total number of constraints.
 132         int                   ncon_tot = 0;
 133         //! The number of flexible constraints.
 134         int                   nflexcon = 0;
 135         //! A list of atoms to constraints for each moleculetype.
 136         std::vector<t_blocka> at2con_mt;
 137         //! A list of atoms to settles for each moleculetype
 138         std::vector < std::vector < int>> at2settle_mt;
 139         //! Whether any SETTLES cross charge-group boundaries.
 140         bool                  bInterCGsettles = false;
 141         //! LINCS data.
 142         Lincs                *lincsd = nullptr; // TODO this should become a unique_ptr
 143         //! SHAKE data.
 144         shakedata            *shaked = nullptr;
 145         //! SETTLE data.
 146         settledata           *settled = nullptr;
 147         //! The maximum number of warnings.
 148         int                   maxwarn = 0;
 149         //! The number of warnings for LINCS.
 150         int                   warncount_lincs = 0;
 151         //! The number of warnings for SETTLE.
 152         int                   warncount_settle = 0;
 153         //! The essential dynamics data.
 154         gmx_edsam *           ed = nullptr;
 155
 156         //! Thread-local virial contribution.
 157         tensor            *vir_r_m_dr_th = {nullptr};
 158         //! Did a settle error occur?
 159         bool              *bSettleErrorHasOccurred = nullptr;
 160
 161         //! Pointer to the global topology - only used for printing warnings.
 162         const gmx_mtop_t  &mtop;
 163         //! Parameters for the interactions in this domain.
 164         const t_idef      *idef = nullptr;
 165         //! Data about atoms in this domain.
 166         const t_mdatoms   &md;
 167         //! Whether we need to do pbc for handling bonds.
 168         bool               pbcHandlingRequired_ = false;
 169
 170         //! Logging support.
 171         FILE                 *log = nullptr;
 172         //! Communication support.
 173         const t_commrec      *cr = nullptr;
 174         //! Multi-sim support.
 175         const gmx_multisim_t *ms = nullptr;
 176         //! Pulling code object, if any.
 177         pull_t               *pull_work = nullptr;
 178         /*!\brief Input options.
 179          *
 180          * \todo Replace with IMdpOptions */
 181         const t_inputrec &ir;
 182         //! Flop counting support.
 183         t_nrnb           *nrnb = nullptr;
 184         //! Tracks wallcycle usage.
 185         gmx_wallcycle    *wcycle;
 186 };
 187
 188 Constraints::~Constraints() = default;
 189
 190 int Constraints::numFlexibleConstraints() const
 191 {
 192     return impl_->nflexcon;
 193 }
 194
 195 bool Constraints::havePerturbedConstraints() const
 196 {
 197     const gmx_ffparams_t &ffparams = impl_->mtop.ffparams;
 198
 199     for (size_t i = 0; i < ffparams.functype.size(); i++)
 200     {
 201         if ((ffparams.functype[i] == F_CONSTR ||
 202              ffparams.functype[i] == F_CONSTRNC) &&
 203             ffparams.iparams[i].constr.dA != ffparams.iparams[i].constr.dB)
 204         {
 205             return true;
 206         }
 207     }
 208
 209     return false;
 210 }
 211
 212 //! Clears constraint quantities for atoms in nonlocal region.
 213 static void clear_constraint_quantity_nonlocal(gmx_domdec_t *dd, rvec *q)
 214 {
 215     int nonlocal_at_start, nonlocal_at_end, at;
 216
 217     dd_get_constraint_range(dd, &nonlocal_at_start, &nonlocal_at_end);
 218
 219     for (at = nonlocal_at_start; at < nonlocal_at_end; at++)
 220     {
 221         clear_rvec(q[at]);
 222     }
 223 }
 224
 225 void too_many_constraint_warnings(int eConstrAlg, int warncount)
 226 {
 227     gmx_fatal(FARGS,
 228               "Too many %s warnings (%d)\n"
 229               "If you know what you are doing you can %s"
 230               "set the environment variable GMX_MAXCONSTRWARN to -1,\n"
 231               "but normally it is better to fix the problem",
 232               (eConstrAlg == econtLINCS) ? "LINCS" : "SETTLE", warncount,
 233               (eConstrAlg == econtLINCS) ?
 234               "adjust the lincs warning threshold in your mdp file\nor " : "\n");
 235 }
 236
 237 //! Writes out coordinates.
 238 static void write_constr_pdb(const char *fn, const char *title,
 239                              const gmx_mtop_t &mtop,
 240                              int start, int homenr, const t_commrec *cr,
 241                              const rvec x[], matrix box)
 242 {
 243     char          fname[STRLEN];
 244     FILE         *out;
 245     int           dd_ac0 = 0, dd_ac1 = 0, i, ii, resnr;
 246     gmx_domdec_t *dd;
 247     const char   *anm, *resnm;
 248
 249     dd = nullptr;
 250     if (DOMAINDECOMP(cr))
 251     {
 252         dd = cr->dd;
 253         dd_get_constraint_range(dd, &dd_ac0, &dd_ac1);
 254         start  = 0;
 255         homenr = dd_ac1;
 256     }
 257
 258     if (PAR(cr))
 259     {
 260         sprintf(fname, "%s_n%d.pdb", fn, cr->sim_nodeid);
 261     }
 262     else
 263     {
 264         sprintf(fname, "%s.pdb", fn);
 265     }
 266
 267     out = gmx_fio_fopen(fname, "w");
 268
 269     fprintf(out, "TITLE     %s\n", title);
 270     gmx_write_pdb_box(out, -1, box);
 271     int molb = 0;
 272     for (i = start; i < start+homenr; i++)
 273     {
 274         if (dd != nullptr)
 275         {
 276             if (i >= dd_numHomeAtoms(*dd) && i < dd_ac0)
 277             {
 278                 continue;
 279             }
 280             ii = dd->globalAtomIndices[i];
 281         }
 282         else
 283         {
 284             ii = i;
 285         }
 286         mtopGetAtomAndResidueName(mtop, ii, &molb, &anm, &resnr, &resnm, nullptr);
 287         gmx_fprintf_pdb_atomline(out, epdbATOM, ii+1, anm, ' ', resnm, ' ', resnr, ' ',
 288                                  10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ], 1.0, 0.0, "");
 289     }
 290     fprintf(out, "TER\n");
 291
 292     gmx_fio_fclose(out);
 293 }
 294
 295 //! Writes out domain contents to help diagnose crashes.
 296 static void dump_confs(FILE *log, int64_t step, const gmx_mtop_t &mtop,
 297                        int start, int homenr, const t_commrec *cr,
 298                        const rvec x[], rvec xprime[], matrix box)
 299 {
 300     char  buf[STRLEN], buf2[22];
 301
 302     char *env = getenv("GMX_SUPPRESS_DUMP");
 303     if (env)
 304     {
 305         return;
 306     }
 307
 308     sprintf(buf, "step%sb", gmx_step_str(step, buf2));
 309     write_constr_pdb(buf, "initial coordinates",
 310                      mtop, start, homenr, cr, x, box);
 311     sprintf(buf, "step%sc", gmx_step_str(step, buf2));
 312     write_constr_pdb(buf, "coordinates after constraining",
 313                      mtop, start, homenr, cr, xprime, box);
 314     if (log)
 315     {
 316         fprintf(log, "Wrote pdb files with previous and current coordinates\n");
 317     }
 318     fprintf(stderr, "Wrote pdb files with previous and current coordinates\n");
 319 }
 320
 321 bool
 322 Constraints::apply(bool                  bLog,
 323                    bool                  bEner,
 324                    int64_t               step,
 325                    int                   delta_step,
 326                    real                  step_scaling,
 327                    rvec                 *x,
 328                    rvec                 *xprime,
 329                    rvec                 *min_proj,
 330                    matrix                box,
 331                    real                  lambda,
 332                    real                 *dvdlambda,
 333                    rvec                 *v,
 334                    tensor               *vir,
 335                    ConstraintVariable    econq)
 336 {
 337     return impl_->apply(bLog,
 338                         bEner,
 339                         step,
 340                         delta_step,
 341                         step_scaling,
 342                         x,
 343                         xprime,
 344                         min_proj,
 345                         box,
 346                         lambda,
 347                         dvdlambda,
 348                         v,
 349                         vir,
 350                         econq);
 351 }
 352
 353 bool
 354 Constraints::Impl::apply(bool                  bLog,
 355                          bool                  bEner,
 356                          int64_t               step,
 357                          int                   delta_step,
 358                          real                  step_scaling,
 359                          rvec                 *x,
 360                          rvec                 *xprime,
 361                          rvec                 *min_proj,
 362                          matrix                box,
 363                          real                  lambda,
 364                          real                 *dvdlambda,
 365                          rvec                 *v,
 366                          tensor               *vir,
 367                          ConstraintVariable    econq)
 368 {
 369     bool        bOK, bDump;
 370     int         start, homenr;
 371     tensor      vir_r_m_dr;
 372     real        scaled_delta_t;
 373     real        invdt, vir_fac = 0, t;
 374     int         nsettle;
 375     t_pbc       pbc, *pbc_null;
 376     char        buf[22];
 377     int         nth;
 378
 379     wallcycle_start(wcycle, ewcCONSTR);
 380
 381     if (econq == ConstraintVariable::ForceDispl && !EI_ENERGY_MINIMIZATION(ir.eI))
 382     {
 383         gmx_incons("constrain called for forces displacements while not doing energy minimization, can not do this while the LINCS and SETTLE constraint connection matrices are mass weighted");
 384     }
 385
 386     bOK   = TRUE;
 387     bDump = FALSE;
 388
 389     start  = 0;
 390     homenr = md.homenr;
 391
 392     scaled_delta_t = step_scaling * ir.delta_t;
 393
 394     /* Prepare time step for use in constraint implementations, and
 395        avoid generating inf when ir.delta_t = 0. */
 396     if (ir.delta_t == 0)
 397     {
 398         invdt = 0.0;
 399     }
 400     else
 401     {
 402         invdt = 1.0/scaled_delta_t;
 403     }
 404
 405     if (ir.efep != efepNO && EI_DYNAMICS(ir.eI))
 406     {
 407         /* Set the constraint lengths for the step at which this configuration
 408          * is meant to be. The invmasses should not be changed.
 409          */
 410         lambda += delta_step*ir.fepvals->delta_lambda;
 411     }
 412
 413     if (vir != nullptr)
 414     {
 415         clear_mat(vir_r_m_dr);
 416     }
 417     const t_ilist *settle = &idef->il[F_SETTLE];
 418     nsettle = settle->nr/(1+NRAL(F_SETTLE));
 419
 420     if (nsettle > 0)
 421     {
 422         nth = gmx_omp_nthreads_get(emntSETTLE);
 423     }
 424     else
 425     {
 426         nth = 1;
 427     }
 428
 429     /* We do not need full pbc when constraints do not cross charge groups,
 430      * i.e. when dd->constraint_comm==NULL.
 431      * Note that PBC for constraints is different from PBC for bondeds.
 432      * For constraints there is both forward and backward communication.
 433      */
 434     if (ir.ePBC != epbcNONE &&
 435         (cr->dd || pbcHandlingRequired_) && !(cr->dd && cr->dd->constraint_comm == nullptr))
 436     {
 437         /* With pbc=screw the screw has been changed to a shift
 438          * by the constraint coordinate communication routine,
 439          * so that here we can use normal pbc.
 440          */
 441         pbc_null = set_pbc_dd(&pbc, ir.ePBC,
 442                               DOMAINDECOMP(cr) ? cr->dd->nc : nullptr,
 443                               FALSE, box);
 444     }
 445     else
 446     {
 447         pbc_null = nullptr;
 448     }
 449
 450     /* Communicate the coordinates required for the non-local constraints
 451      * for LINCS and/or SETTLE.
 452      */
 453     if (cr->dd)
 454     {
 455         dd_move_x_constraints(cr->dd, box, x, xprime, econq == ConstraintVariable::Positions);
 456
 457         if (v != nullptr)
 458         {
 459             /* We need to initialize the non-local components of v.
 460              * We never actually use these values, but we do increment them,
 461              * so we should avoid uninitialized variables and overflows.
 462              */
 463             clear_constraint_quantity_nonlocal(cr->dd, v);
 464         }
 465     }
 466
 467     if (lincsd != nullptr)
 468     {
 469         bOK = constrain_lincs(bLog || bEner, ir, step, lincsd, md, cr, ms,
 470                               x, xprime, min_proj,
 471                               box, pbc_null, lambda, dvdlambda,
 472                               invdt, v, vir != nullptr, vir_r_m_dr,
 473                               econq, nrnb,
 474                               maxwarn, &warncount_lincs);
 475         if (!bOK && maxwarn < INT_MAX)
 476         {
 477             if (log != nullptr)
 478             {
 479                 fprintf(log, "Constraint error in algorithm %s at step %s\n",
 480                         econstr_names[econtLINCS], gmx_step_str(step, buf));
 481             }
 482             bDump = TRUE;
 483         }
 484     }
 485
 486     if (shaked != nullptr)
 487     {
 488         bOK = constrain_shake(log, shaked,
 489                               md.invmass,
 490                               *idef, ir, x, xprime, min_proj, nrnb,
 491                               lambda, dvdlambda,
 492                               invdt, v, vir != nullptr, vir_r_m_dr,
 493                               maxwarn < INT_MAX, econq);
 494
 495         if (!bOK && maxwarn < INT_MAX)
 496         {
 497             if (log != nullptr)
 498             {
 499                 fprintf(log, "Constraint error in algorithm %s at step %s\n",
 500                         econstr_names[econtSHAKE], gmx_step_str(step, buf));
 501             }
 502             bDump = TRUE;
 503         }
 504     }
 505
 506     if (nsettle > 0)
 507     {
 508         bool bSettleErrorHasOccurred0 = false;
 509
 510         switch (econq)
 511         {
 512             case ConstraintVariable::Positions:
 513 #pragma omp parallel for num_threads(nth) schedule(static)
 514                 for (int th = 0; th < nth; th++)
 515                 {
 516                     try
 517                     {
 518                         if (th > 0)
 519                         {
 520                             clear_mat(vir_r_m_dr_th[th]);
 521                         }
 522
 523                         csettle(settled,
 524                                 nth, th,
 525                                 pbc_null,
 526                                 x[0], xprime[0],
 527                                 invdt, v ? v[0] : nullptr,
 528                                 vir != nullptr,
 529                                 th == 0 ? vir_r_m_dr : vir_r_m_dr_th[th],
 530                                 th == 0 ? &bSettleErrorHasOccurred0 : &bSettleErrorHasOccurred[th]);
 531                     }
 532                     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
 533                 }
 534                 inc_nrnb(nrnb, eNR_SETTLE, nsettle);
 535                 if (v != nullptr)
 536                 {
 537                     inc_nrnb(nrnb, eNR_CONSTR_V, nsettle*3);
 538                 }
 539                 if (vir != nullptr)
 540                 {
 541                     inc_nrnb(nrnb, eNR_CONSTR_VIR, nsettle*3);
 542                 }
 543                 break;
 544             case ConstraintVariable::Velocities:
 545             case ConstraintVariable::Derivative:
 546             case ConstraintVariable::Force:
 547             case ConstraintVariable::ForceDispl:
 548 #pragma omp parallel for num_threads(nth) schedule(static)
 549                 for (int th = 0; th < nth; th++)
 550                 {
 551                     try
 552                     {
 553                         int calcvir_atom_end;
 554
 555                         if (vir == nullptr)
 556                         {
 557                             calcvir_atom_end = 0;
 558                         }
 559                         else
 560                         {
 561                             calcvir_atom_end = md.homenr;
 562                         }
 563
 564                         if (th > 0)
 565                         {
 566                             clear_mat(vir_r_m_dr_th[th]);
 567                         }
 568
 569                         int start_th = (nsettle* th   )/nth;
 570                         int end_th   = (nsettle*(th+1))/nth;
 571
 572                         if (start_th >= 0 && end_th - start_th > 0)
 573                         {
 574                             settle_proj(settled, econq,
 575                                         end_th-start_th,
 576                                         settle->iatoms+start_th*(1+NRAL(F_SETTLE)),
 577                                         pbc_null,
 578                                         x,
 579                                         xprime, min_proj, calcvir_atom_end,
 580                                         th == 0 ? vir_r_m_dr : vir_r_m_dr_th[th]);
 581                         }
 582                     }
 583                     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
 584                 }
 585                 /* This is an overestimate */
 586                 inc_nrnb(nrnb, eNR_SETTLE, nsettle);
 587                 break;
 588             case ConstraintVariable::Deriv_FlexCon:
 589                 /* Nothing to do, since the are no flexible constraints in settles */
 590                 break;
 591             default:
 592                 gmx_incons("Unknown constraint quantity for settle");
 593         }
 594
 595         if (vir != nullptr)
 596         {
 597             /* Reduce the virial contributions over the threads */
 598             for (int th = 1; th < nth; th++)
 599             {
 600                 m_add(vir_r_m_dr, vir_r_m_dr_th[th], vir_r_m_dr);
 601             }
 602         }
 603
 604         if (econq == ConstraintVariable::Positions)
 605         {
 606             for (int th = 1; th < nth; th++)
 607             {
 608                 bSettleErrorHasOccurred0 = bSettleErrorHasOccurred0 || bSettleErrorHasOccurred[th];
 609             }
 610
 611             if (bSettleErrorHasOccurred0)
 612             {
 613                 char buf[STRLEN];
 614                 sprintf(buf,
 615                         "\nstep " "%" PRId64 ": One or more water molecules can not be settled.\n"
 616                         "Check for bad contacts and/or reduce the timestep if appropriate.\n",
 617                         step);
 618                 if (log)
 619                 {
 620                     fprintf(log, "%s", buf);
 621                 }
 622                 fprintf(stderr, "%s", buf);
 623                 warncount_settle++;
 624                 if (warncount_settle > maxwarn)
 625                 {
 626                     too_many_constraint_warnings(-1, warncount_settle);
 627                 }
 628                 bDump = TRUE;
 629
 630                 bOK   = FALSE;
 631             }
 632         }
 633     }
 634
 635     if (vir != nullptr)
 636     {
 637         /* The normal uses of constrain() pass step_scaling = 1.0.
 638          * The call to constrain() for SD1 that passes step_scaling =
 639          * 0.5 also passes vir = NULL, so cannot reach this
 640          * assertion. This assertion should remain until someone knows
 641          * that this path works for their intended purpose, and then
 642          * they can use scaled_delta_t instead of ir.delta_t
 643          * below. */
 644         assert(gmx_within_tol(step_scaling, 1.0, GMX_REAL_EPS));
 645         switch (econq)
 646         {
 647             case ConstraintVariable::Positions:
 648                 vir_fac = 0.5/(ir.delta_t*ir.delta_t);
 649                 break;
 650             case ConstraintVariable::Velocities:
 651                 vir_fac = 0.5/ir.delta_t;
 652                 break;
 653             case ConstraintVariable::Force:
 654             case ConstraintVariable::ForceDispl:
 655                 vir_fac = 0.5;
 656                 break;
 657             default:
 658                 gmx_incons("Unsupported constraint quantity for virial");
 659         }
 660
 661         if (EI_VV(ir.eI))
 662         {
 663             vir_fac *= 2;  /* only constraining over half the distance here */
 664         }
 665         for (int i = 0; i < DIM; i++)
 666         {
 667             for (int j = 0; j < DIM; j++)
 668             {
 669                 (*vir)[i][j] = vir_fac*vir_r_m_dr[i][j];
 670             }
 671         }
 672     }
 673
 674     if (bDump)
 675     {
 676         dump_confs(log, step, mtop, start, homenr, cr, x, xprime, box);
 677     }
 678
 679     if (econq == ConstraintVariable::Positions)
 680     {
 681         if (ir.bPull && pull_have_constraint(pull_work))
 682         {
 683             if (EI_DYNAMICS(ir.eI))
 684             {
 685                 t = ir.init_t + (step + delta_step)*ir.delta_t;
 686             }
 687             else
 688             {
 689                 t = ir.init_t;
 690             }
 691             set_pbc(&pbc, ir.ePBC, box);
 692             pull_constraint(pull_work, &md, &pbc, cr, ir.delta_t, t, x, xprime, v, *vir);
 693         }
 694         if (ed && delta_step > 0)
 695         {
 696             /* apply the essential dynamics constraints here */
 697             do_edsam(&ir, step, cr, xprime, v, box, ed);
 698         }
 699     }
 700     wallcycle_stop(wcycle, ewcCONSTR);
 701
 702     if (v != nullptr && md.cFREEZE)
 703     {
 704         /* Set the velocities of frozen dimensions to zero */
 705
 706         int gmx_unused numThreads = gmx_omp_nthreads_get(emntUpdate);
 707
 708 #pragma omp parallel for num_threads(numThreads) schedule(static)
 709         for (int i = 0; i < md.homenr; i++)
 710         {
 711             int freezeGroup = md.cFREEZE[i];
 712
 713             for (int d = 0; d < DIM; d++)
 714             {
 715                 if (ir.opts.nFreeze[freezeGroup][d])
 716                 {
 717                     v[i][d] = 0;
 718                 }
 719             }
 720         }
 721     }
 722
 723     return bOK;
 724 }
 725
 726 ArrayRef<real> Constraints::rmsdData() const
 727 {
 728     if (impl_->lincsd)
 729     {
 730         return lincs_rmsdData(impl_->lincsd);
 731     }
 732     else
 733     {
 734         return {};
 735     }
 736 }
 737
 738 real Constraints::rmsd() const
 739 {
 740     if (impl_->lincsd)
 741     {
 742         return lincs_rmsd(impl_->lincsd);
 743     }
 744     else
 745     {
 746         return 0;
 747     }
 748 }
 749
 750 FlexibleConstraintTreatment flexibleConstraintTreatment(bool haveDynamicsIntegrator)
 751 {
 752     if (haveDynamicsIntegrator)
 753     {
 754         return FlexibleConstraintTreatment::Include;
 755     }
 756     else
 757     {
 758         return FlexibleConstraintTreatment::Exclude;
 759     }
 760 }
 761
 762 /*! \brief Returns a block struct to go from atoms to constraints
 763  *
 764  * The block struct will contain constraint indices with lower indices
 765  * directly matching the order in F_CONSTR and higher indices matching
 766  * the order in F_CONSTRNC offset by the number of constraints in F_CONSTR.
 767  *
 768  * \param[in]  numAtoms  The number of atoms to construct the list for
 769  * \param[in]  ilists    The interaction lists, size F_NRE
 770  * \param[in]  iparams   Interaction parameters, can be null when flexibleConstraintTreatment=Include
 771  * \param[in]  flexibleConstraintTreatment  The flexible constraint treatment, see enum above
 772  * \returns a block struct with all constraints for each atom
 773  */
 774 template <typename T>
 775 static t_blocka
 776 makeAtomsToConstraintsList(int                          numAtoms,
 777                            const T                     *ilists,
 778                            const t_iparams             *iparams,
 779                            FlexibleConstraintTreatment  flexibleConstraintTreatment)
 780 {
 781     GMX_ASSERT(flexibleConstraintTreatment == FlexibleConstraintTreatment::Include || iparams != nullptr, "With flexible constraint detection we need valid iparams");
 782
 783     std::vector<int> count(numAtoms);
 784
 785     for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
 786     {
 787         const T   &ilist  = ilists[ftype];
 788         const int  stride = 1 + NRAL(ftype);
 789         for (int i = 0; i < ilist.size(); i += stride)
 790         {
 791             if (flexibleConstraintTreatment == FlexibleConstraintTreatment::Include ||
 792                 !isConstraintFlexible(iparams, ilist.iatoms[i]))
 793             {
 794                 for (int j = 1; j < 3; j++)
 795                 {
 796                     int a = ilist.iatoms[i + j];
 797                     count[a]++;
 798                 }
 799             }
 800         }
 801     }
 802
 803     t_blocka at2con;
 804     at2con.nr           = numAtoms;
 805     at2con.nalloc_index = at2con.nr + 1;
 806     snew(at2con.index, at2con.nalloc_index);
 807     at2con.index[0] = 0;
 808     for (int a = 0; a < numAtoms; a++)
 809     {
 810         at2con.index[a + 1] = at2con.index[a] + count[a];
 811         count[a]            = 0;
 812     }
 813     at2con.nra      = at2con.index[at2con.nr];
 814     at2con.nalloc_a = at2con.nra;
 815     snew(at2con.a, at2con.nalloc_a);
 816
 817     /* The F_CONSTRNC constraints have constraint numbers
 818      * that continue after the last F_CONSTR constraint.
 819      */
 820     int numConstraints = 0;
 821     for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
 822     {
 823         const T   &ilist  = ilists[ftype];
 824         const int  stride = 1 + NRAL(ftype);
 825         for (int i = 0; i < ilist.size(); i += stride)
 826         {
 827             if (flexibleConstraintTreatment == FlexibleConstraintTreatment::Include ||
 828                 !isConstraintFlexible(iparams, ilist.iatoms[i]))
 829             {
 830                 for (int j = 1; j < 3; j++)
 831                 {
 832                     int a = ilist.iatoms[i + j];
 833                     at2con.a[at2con.index[a] + count[a]++] = numConstraints;
 834                 }
 835             }
 836             numConstraints++;
 837         }
 838     }
 839
 840     return at2con;
 841 }
 842
 843 t_blocka make_at2con(int                          numAtoms,
 844                      const t_ilist               *ilist,
 845                      const t_iparams             *iparams,
 846                      FlexibleConstraintTreatment  flexibleConstraintTreatment)
 847 {
 848     return makeAtomsToConstraintsList(numAtoms, ilist, iparams, flexibleConstraintTreatment);
 849 }
 850
 851 t_blocka make_at2con(const gmx_moltype_t            &moltype,
 852                      gmx::ArrayRef<const t_iparams>  iparams,
 853                      FlexibleConstraintTreatment     flexibleConstraintTreatment)
 854 {
 855     return makeAtomsToConstraintsList(moltype.atoms.nr, moltype.ilist.data(), iparams.data(), flexibleConstraintTreatment);
 856 }
 857
 858 //! Return the number of flexible constraints in the \c ilist and \c iparams.
 859 template <typename T>
 860 static int
 861 countFlexibleConstraintsTemplate(const T         *ilist,
 862                                  const t_iparams *iparams)
 863 {
 864     int nflexcon = 0;
 865     for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
 866     {
 867         const int numIatomsPerConstraint = 3;
 868         for (int i = 0; i < ilist[ftype].size(); i += numIatomsPerConstraint)
 869         {
 870             const int type = ilist[ftype].iatoms[i];
 871             if (iparams[type].constr.dA == 0 &&
 872                 iparams[type].constr.dB == 0)
 873             {
 874                 nflexcon++;
 875             }
 876         }
 877     }
 878
 879     return nflexcon;
 880 }
 881
 882 int countFlexibleConstraints(const t_ilist   *ilist,
 883                              const t_iparams *iparams)
 884 {
 885     return countFlexibleConstraintsTemplate(ilist, iparams);
 886 }
 887
 888 //! Returns the index of the settle to which each atom belongs.
 889 static std::vector<int> make_at2settle(int                    natoms,
 890                                        const InteractionList &ilist)
 891 {
 892     /* Set all to no settle */
 893     std::vector<int> at2s(natoms, -1);
 894
 895     const int        stride = 1 + NRAL(F_SETTLE);
 896
 897     for (int s = 0; s < ilist.size(); s += stride)
 898     {
 899         at2s[ilist.iatoms[s + 1]] = s/stride;
 900         at2s[ilist.iatoms[s + 2]] = s/stride;
 901         at2s[ilist.iatoms[s + 3]] = s/stride;
 902     }
 903
 904     return at2s;
 905 }
 906
 907 void
 908 Constraints::Impl::setConstraints(const gmx_localtop_t &top,
 909                                   const t_mdatoms      &md)
 910 {
 911     idef = &top.idef;
 912
 913     if (ncon_tot > 0)
 914     {
 915         /* With DD we might also need to call LINCS on a domain no constraints for
 916          * communicating coordinates to other nodes that do have constraints.
 917          */
 918         if (ir.eConstrAlg == econtLINCS)
 919         {
 920             set_lincs(top.idef, md, EI_DYNAMICS(ir.eI), cr, lincsd);
 921         }
 922         if (ir.eConstrAlg == econtSHAKE)
 923         {
 924             if (cr->dd)
 925             {
 926                 // We are using the local topology, so there are only
 927                 // F_CONSTR constraints.
 928                 make_shake_sblock_dd(shaked, &idef->il[F_CONSTR], cr->dd);
 929             }
 930             else
 931             {
 932                 make_shake_sblock_serial(shaked, idef, md);
 933             }
 934         }
 935     }
 936
 937     if (settled)
 938     {
 939         settle_set_constraints(settled,
 940                                &idef->il[F_SETTLE], md);
 941     }
 942
 943     /* Make a selection of the local atoms for essential dynamics */
 944     if (ed && cr->dd)
 945     {
 946         dd_make_local_ed_indices(cr->dd, ed);
 947     }
 948 }
 949
 950 void
 951 Constraints::setConstraints(const gmx_localtop_t &top,
 952                             const t_mdatoms      &md)
 953 {
 954     impl_->setConstraints(top, md);
 955 }
 956
 957 /*! \brief Makes a per-moleculetype container of mappings from atom
 958  * indices to constraint indices.
 959  *
 960  * Note that flexible constraints are only enabled with a dynamical integrator. */
 961 static std::vector<t_blocka>
 962 makeAtomToConstraintMappings(const gmx_mtop_t            &mtop,
 963                              FlexibleConstraintTreatment  flexibleConstraintTreatment)
 964 {
 965     std::vector<t_blocka> mapping;
 966     mapping.reserve(mtop.moltype.size());
 967     for (const gmx_moltype_t &moltype : mtop.moltype)
 968     {
 969         mapping.push_back(make_at2con(moltype,
 970                                       mtop.ffparams.iparams,
 971                                       flexibleConstraintTreatment));
 972     }
 973     return mapping;
 974 }
 975
 976 Constraints::Constraints(const gmx_mtop_t     &mtop,
 977                          const t_inputrec     &ir,
 978                          pull_t               *pull_work,
 979                          FILE                 *log,
 980                          const t_mdatoms      &md,
 981                          const t_commrec      *cr,
 982                          const gmx_multisim_t *ms,
 983                          t_nrnb               *nrnb,
 984                          gmx_wallcycle        *wcycle,
 985                          bool                  pbcHandlingRequired,
 986                          int                   numConstraints,
 987                          int                   numSettles)
 988     : impl_(new Impl(mtop,
 989                      ir,
 990                      pull_work,
 991                      log,
 992                      md,
 993                      cr,
 994                      ms,
 995                      nrnb,
 996                      wcycle,
 997                      pbcHandlingRequired,
 998                      numConstraints,
 999                      numSettles))
1000 {
1001 }
1002
1003 Constraints::Impl::Impl(const gmx_mtop_t     &mtop_p,
1004                         const t_inputrec     &ir_p,
1005                         pull_t               *pull_work,
1006                         FILE                 *log_p,
1007                         const t_mdatoms      &md_p,
1008                         const t_commrec      *cr_p,
1009                         const gmx_multisim_t *ms_p,
1010                         t_nrnb               *nrnb_p,
1011                         gmx_wallcycle        *wcycle_p,
1012                         bool                  pbcHandlingRequired,
1013                         int                   numConstraints,
1014                         int                   numSettles)
1015     : ncon_tot(numConstraints),
1016       mtop(mtop_p),
1017       md(md_p),
1018       pbcHandlingRequired_(pbcHandlingRequired),
1019       log(log_p),
1020       cr(cr_p),
1021       ms(ms_p),
1022       pull_work(pull_work),
1023       ir(ir_p),
1024       nrnb(nrnb_p),
1025       wcycle(wcycle_p)
1026 {
1027     if (numConstraints + numSettles > 0 && ir.epc == epcMTTK)
1028     {
1029         gmx_fatal(FARGS, "Constraints are not implemented with MTTK pressure control.");
1030     }
1031
1032     nflexcon = 0;
1033     if (numConstraints > 0)
1034     {
1035         at2con_mt = makeAtomToConstraintMappings(mtop,
1036                                                  flexibleConstraintTreatment(EI_DYNAMICS(ir.eI)));
1037
1038         for (const gmx_molblock_t &molblock : mtop.molblock)
1039         {
1040             int count =
1041                 countFlexibleConstraintsTemplate(mtop.moltype[molblock.type].ilist.data(),
1042                                                  mtop.ffparams.iparams.data());
1043             nflexcon += molblock.nmol*count;
1044         }
1045
1046         if (nflexcon > 0)
1047         {
1048             if (log)
1049             {
1050                 fprintf(log, "There are %d flexible constraints\n",
1051                         nflexcon);
1052                 if (ir.fc_stepsize == 0)
1053                 {
1054                     fprintf(log, "\n"
1055                             "WARNING: step size for flexible constraining = 0\n"
1056                             "         All flexible constraints will be rigid.\n"
1057                             "         Will try to keep all flexible constraints at their original length,\n"
1058                             "         but the lengths may exhibit some drift.\n\n");
1059                     nflexcon = 0;
1060                 }
1061             }
1062             if (nflexcon > 0)
1063             {
1064                 please_cite(log, "Hess2002");
1065             }
1066         }
1067
1068         if (ir.eConstrAlg == econtLINCS)
1069         {
1070             lincsd = init_lincs(log, mtop,
1071                                 nflexcon, at2con_mt,
1072                                 DOMAINDECOMP(cr) && cr->dd->splitConstraints,
1073                                 ir.nLincsIter, ir.nProjOrder);
1074         }
1075
1076         if (ir.eConstrAlg == econtSHAKE)
1077         {
1078             if (DOMAINDECOMP(cr) && cr->dd->splitConstraints)
1079             {
1080                 gmx_fatal(FARGS, "SHAKE is not supported with domain decomposition and constraint that cross charge group boundaries, use LINCS");
1081             }
1082             if (nflexcon)
1083             {
1084                 gmx_fatal(FARGS, "For this system also velocities and/or forces need to be constrained, this can not be done with SHAKE, you should select LINCS");
1085             }
1086             please_cite(log, "Ryckaert77a");
1087             if (ir.bShakeSOR)
1088             {
1089                 please_cite(log, "Barth95a");
1090             }
1091
1092             shaked = shake_init();
1093         }
1094     }
1095
1096     if (numSettles > 0)
1097     {
1098         please_cite(log, "Miyamoto92a");
1099
1100         bInterCGsettles = inter_charge_group_settles(mtop);
1101
1102         settled         = settle_init(mtop);
1103
1104         /* Make an atom to settle index for use in domain decomposition */
1105         for (size_t mt = 0; mt < mtop.moltype.size(); mt++)
1106         {
1107             at2settle_mt.emplace_back(make_at2settle(mtop.moltype[mt].atoms.nr,
1108                                                      mtop.moltype[mt].ilist[F_SETTLE]));
1109         }
1110
1111         /* Allocate thread-local work arrays */
1112         int nthreads = gmx_omp_nthreads_get(emntSETTLE);
1113         if (nthreads > 1 && vir_r_m_dr_th == nullptr)
1114         {
1115             snew(vir_r_m_dr_th, nthreads);
1116             snew(bSettleErrorHasOccurred, nthreads);
1117         }
1118     }
1119
1120     maxwarn = 999;
1121     char *env       = getenv("GMX_MAXCONSTRWARN");
1122     if (env)
1123     {
1124         maxwarn = 0;
1125         sscanf(env, "%8d", &maxwarn);
1126         if (maxwarn < 0)
1127         {
1128             maxwarn = INT_MAX;
1129         }
1130         if (log)
1131         {
1132             fprintf(log,
1133                     "Setting the maximum number of constraint warnings to %d\n",
1134                     maxwarn);
1135         }
1136         if (MASTER(cr))
1137         {
1138             fprintf(stderr,
1139                     "Setting the maximum number of constraint warnings to %d\n",
1140                     maxwarn);
1141         }
1142     }
1143     warncount_lincs  = 0;
1144     warncount_settle = 0;
1145 }
1146
1147 Constraints::Impl::~Impl()
1148 {
1149     for (auto blocka : at2con_mt)
1150     {
1151         done_blocka(&blocka);
1152     }
1153     if (bSettleErrorHasOccurred != nullptr)
1154     {
1155         sfree(bSettleErrorHasOccurred);
1156     }
1157     if (vir_r_m_dr_th != nullptr)
1158     {
1159         sfree(vir_r_m_dr_th);
1160     }
1161     if (settled != nullptr)
1162     {
1163         settle_free(settled);
1164     }
1165     done_lincs(lincsd);
1166 }
1167
1168 void Constraints::saveEdsamPointer(gmx_edsam * ed)
1169 {
1170     impl_->ed = ed;
1171 }
1172
1173 ArrayRef<const t_blocka>
1174 Constraints::atom2constraints_moltype() const
1175 {
1176     return impl_->at2con_mt;
1177 }
1178
1179 ArrayRef < const std::vector < int>> Constraints::atom2settle_moltype() const
1180 {
1181     return impl_->at2settle_mt;
1182 }
1183
1184
1185 bool inter_charge_group_constraints(const gmx_mtop_t &mtop)
1186 {
1187     const gmx_moltype_t *molt;
1188     const t_block       *cgs;
1189     int                 *at2cg, cg, a, ftype, i;
1190     bool                 bInterCG;
1191
1192     bInterCG = FALSE;
1193     for (size_t mb = 0; mb < mtop.molblock.size() && !bInterCG; mb++)
1194     {
1195         molt = &mtop.moltype[mtop.molblock[mb].type];
1196
1197         if (molt->ilist[F_CONSTR].size()   > 0 ||
1198             molt->ilist[F_CONSTRNC].size() > 0 ||
1199             molt->ilist[F_SETTLE].size()   > 0)
1200         {
1201             cgs  = &molt->cgs;
1202             snew(at2cg, molt->atoms.nr);
1203             for (cg = 0; cg < cgs->nr; cg++)
1204             {
1205                 for (a = cgs->index[cg]; a < cgs->index[cg+1]; a++)
1206                 {
1207                     at2cg[a] = cg;
1208                 }
1209             }
1210
1211             for (ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
1212             {
1213                 const InteractionList &il = molt->ilist[ftype];
1214                 for (i = 0; i < il.size() && !bInterCG; i += 1+NRAL(ftype))
1215                 {
1216                     if (at2cg[il.iatoms[i + 1]] != at2cg[il.iatoms[i + 2]])
1217                     {
1218                         bInterCG = TRUE;
1219                     }
1220                 }
1221             }
1222
1223             sfree(at2cg);
1224         }
1225     }
1226
1227     return bInterCG;
1228 }
1229
1230 bool inter_charge_group_settles(const gmx_mtop_t &mtop)
1231 {
1232     const gmx_moltype_t *molt;
1233     const t_block       *cgs;
1234     int                 *at2cg, cg, a, ftype, i;
1235     bool                 bInterCG;
1236
1237     bInterCG = FALSE;
1238     for (size_t mb = 0; mb < mtop.molblock.size() && !bInterCG; mb++)
1239     {
1240         molt = &mtop.moltype[mtop.molblock[mb].type];
1241
1242         if (molt->ilist[F_SETTLE].size() > 0)
1243         {
1244             cgs  = &molt->cgs;
1245             snew(at2cg, molt->atoms.nr);
1246             for (cg = 0; cg < cgs->nr; cg++)
1247             {
1248                 for (a = cgs->index[cg]; a < cgs->index[cg+1]; a++)
1249                 {
1250                     at2cg[a] = cg;
1251                 }
1252             }
1253
1254             for (ftype = F_SETTLE; ftype <= F_SETTLE; ftype++)
1255             {
1256                 const InteractionList &il = molt->ilist[ftype];
1257                 for (i = 0; i < il.size() && !bInterCG; i += 1+NRAL(F_SETTLE))
1258                 {
1259                     if (at2cg[il.iatoms[i + 1]] != at2cg[il.iatoms[i + 2]] ||
1260                         at2cg[il.iatoms[i + 1]] != at2cg[il.iatoms[i + 3]])
1261                     {
1262                         bInterCG = TRUE;
1263                     }
1264                 }
1265             }
1266
1267             sfree(at2cg);
1268         }
1269     }
1270
1271     return bInterCG;
1272 }
1273
1274 void do_constrain_first(FILE *fplog, gmx::Constraints *constr,
1275                         const t_inputrec *ir, const t_mdatoms *md,
1276                         t_state *state)
1277 {
1278     int             i, m, start, end;
1279     int64_t         step;
1280     real            dt = ir->delta_t;
1281     real            dvdl_dum;
1282     rvec           *savex;
1283
1284     /* We need to allocate one element extra, since we might use
1285      * (unaligned) 4-wide SIMD loads to access rvec entries.
1286      */
1287     snew(savex, state->natoms + 1);
1288
1289     start = 0;
1290     end   = md->homenr;
1291
1292     if (debug)
1293     {
1294         fprintf(debug, "vcm: start=%d, homenr=%d, end=%d\n",
1295                 start, md->homenr, end);
1296     }
1297     /* Do a first constrain to reset particles... */
1298     step = ir->init_step;
1299     if (fplog)
1300     {
1301         char buf[STEPSTRSIZE];
1302         fprintf(fplog, "\nConstraining the starting coordinates (step %s)\n",
1303                 gmx_step_str(step, buf));
1304     }
1305     dvdl_dum = 0;
1306
1307     /* constrain the current position */
1308     constr->apply(TRUE, FALSE,
1309                   step, 0, 1.0,
1310                   state->x.rvec_array(), state->x.rvec_array(), nullptr,
1311                   state->box,
1312                   state->lambda[efptBONDED], &dvdl_dum,
1313                   nullptr, nullptr, gmx::ConstraintVariable::Positions);
1314     if (EI_VV(ir->eI))
1315     {
1316         /* constrain the inital velocity, and save it */
1317         /* also may be useful if we need the ekin from the halfstep for velocity verlet */
1318         constr->apply(TRUE, FALSE,
1319                       step, 0, 1.0,
1320                       state->x.rvec_array(), state->v.rvec_array(), state->v.rvec_array(),
1321                       state->box,
1322                       state->lambda[efptBONDED], &dvdl_dum,
1323                       nullptr, nullptr, gmx::ConstraintVariable::Velocities);
1324     }
1325     /* constrain the inital velocities at t-dt/2 */
1326     if (EI_STATE_VELOCITY(ir->eI) && ir->eI != eiVV)
1327     {
1328         auto x = makeArrayRef(state->x).subArray(start, end);
1329         auto v = makeArrayRef(state->v).subArray(start, end);
1330         for (i = start; (i < end); i++)
1331         {
1332             for (m = 0; (m < DIM); m++)
1333             {
1334                 /* Reverse the velocity */
1335                 v[i][m] = -v[i][m];
1336                 /* Store the position at t-dt in buf */
1337                 savex[i][m] = x[i][m] + dt*v[i][m];
1338             }
1339         }
1340         /* Shake the positions at t=-dt with the positions at t=0
1341          * as reference coordinates.
1342          */
1343         if (fplog)
1344         {
1345             char buf[STEPSTRSIZE];
1346             fprintf(fplog, "\nConstraining the coordinates at t0-dt (step %s)\n",
1347                     gmx_step_str(step, buf));
1348         }
1349         dvdl_dum = 0;
1350         constr->apply(TRUE, FALSE,
1351                       step, -1, 1.0,
1352                       state->x.rvec_array(), savex, nullptr,
1353                       state->box,
1354                       state->lambda[efptBONDED], &dvdl_dum,
1355                       state->v.rvec_array(), nullptr, gmx::ConstraintVariable::Positions);
1356
1357         for (i = start; i < end; i++)
1358         {
1359             for (m = 0; m < DIM; m++)
1360             {
1361                 /* Re-reverse the velocities */
1362                 v[i][m] = -v[i][m];
1363             }
1364         }
1365     }
1366     sfree(savex);
1367 }
1368
1369 }  // namespace gmx