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