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, 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", warncount,
 247             (eConstrAlg == econtLINCS) ? "adjust the lincs warning threshold in your mdp file\nor " : "\n");
 248 }
 249
 250 //! Writes out coordinates.
 251 static void write_constr_pdb(const char*          fn,
 252                              const char*          title,
 253                              const gmx_mtop_t&    mtop,
 254                              int                  start,
 255                              int                  homenr,
 256                              const t_commrec*     cr,
 257                              ArrayRef<const RVec> x,
 258                              const matrix         box)
 259 {
 260     char          fname[STRLEN];
 261     FILE*         out;
 262     int           dd_ac0 = 0, dd_ac1 = 0, i, ii, resnr;
 263     gmx_domdec_t* dd;
 264     const char *  anm, *resnm;
 265
 266     dd = nullptr;
 267     if (DOMAINDECOMP(cr))
 268     {
 269         dd = cr->dd;
 270         dd_get_constraint_range(dd, &dd_ac0, &dd_ac1);
 271         start  = 0;
 272         homenr = dd_ac1;
 273     }
 274
 275     if (PAR(cr))
 276     {
 277         sprintf(fname, "%s_n%d.pdb", fn, cr->sim_nodeid);
 278     }
 279     else
 280     {
 281         sprintf(fname, "%s.pdb", fn);
 282     }
 283
 284     out = gmx_fio_fopen(fname, "w");
 285
 286     fprintf(out, "TITLE     %s\n", title);
 287     gmx_write_pdb_box(out, PbcType::Unset, box);
 288     int molb = 0;
 289     for (i = start; i < start + homenr; i++)
 290     {
 291         if (dd != nullptr)
 292         {
 293             if (i >= dd_numHomeAtoms(*dd) && i < dd_ac0)
 294             {
 295                 continue;
 296             }
 297             ii = dd->globalAtomIndices[i];
 298         }
 299         else
 300         {
 301             ii = i;
 302         }
 303         mtopGetAtomAndResidueName(mtop, ii, &molb, &anm, &resnr, &resnm, nullptr);
 304         gmx_fprintf_pdb_atomline(out, epdbATOM, ii + 1, anm, ' ', resnm, ' ', resnr, ' ',
 305                                  10 * x[i][XX], 10 * x[i][YY], 10 * x[i][ZZ], 1.0, 0.0, "");
 306     }
 307     fprintf(out, "TER\n");
 308
 309     gmx_fio_fclose(out);
 310 }
 311
 312 //! Writes out domain contents to help diagnose crashes.
 313 static void dump_confs(FILE*                log,
 314                        int64_t              step,
 315                        const gmx_mtop_t&    mtop,
 316                        int                  start,
 317                        int                  homenr,
 318                        const t_commrec*     cr,
 319                        ArrayRef<const RVec> x,
 320                        ArrayRef<const RVec> xprime,
 321                        const matrix         box)
 322 {
 323     char buf[STRLEN], buf2[22];
 324
 325     char* env = getenv("GMX_SUPPRESS_DUMP");
 326     if (env)
 327     {
 328         return;
 329     }
 330
 331     sprintf(buf, "step%sb", gmx_step_str(step, buf2));
 332     write_constr_pdb(buf, "initial coordinates", mtop, start, homenr, cr, x, box);
 333     sprintf(buf, "step%sc", gmx_step_str(step, buf2));
 334     write_constr_pdb(buf, "coordinates after constraining", mtop, start, homenr, cr, xprime, box);
 335     if (log)
 336     {
 337         fprintf(log, "Wrote pdb files with previous and current coordinates\n");
 338     }
 339     fprintf(stderr, "Wrote pdb files with previous and current coordinates\n");
 340 }
 341
 342 bool Constraints::apply(bool                      bLog,
 343                         bool                      bEner,
 344                         int64_t                   step,
 345                         int                       delta_step,
 346                         real                      step_scaling,
 347                         ArrayRefWithPadding<RVec> x,
 348                         ArrayRefWithPadding<RVec> xprime,
 349                         ArrayRef<RVec>            min_proj,
 350                         const matrix              box,
 351                         real                      lambda,
 352                         real*                     dvdlambda,
 353                         ArrayRefWithPadding<RVec> v,
 354                         bool                      computeVirial,
 355                         tensor                    constraintsVirial,
 356                         ConstraintVariable        econq)
 357 {
 358     return impl_->apply(bLog, bEner, step, delta_step, step_scaling, std::move(x),
 359                         std::move(xprime), min_proj, box, lambda, dvdlambda, std::move(v),
 360                         computeVirial, constraintsVirial, econq);
 361 }
 362
 363 bool Constraints::Impl::apply(bool                      bLog,
 364                               bool                      bEner,
 365                               int64_t                   step,
 366                               int                       delta_step,
 367                               real                      step_scaling,
 368                               ArrayRefWithPadding<RVec> x,
 369                               ArrayRefWithPadding<RVec> xprime,
 370                               ArrayRef<RVec>            min_proj,
 371                               const matrix              box,
 372                               real                      lambda,
 373                               real*                     dvdlambda,
 374                               ArrayRefWithPadding<RVec> v,
 375                               bool                      computeVirial,
 376                               tensor                    constraintsVirial,
 377                               ConstraintVariable        econq)
 378 {
 379     bool  bOK, bDump;
 380     int   start;
 381     real  scaled_delta_t;
 382     real  invdt, vir_fac = 0, t;
 383     int   nsettle;
 384     t_pbc pbc, *pbc_null;
 385     char  buf[22];
 386     int   nth;
 387
 388     wallcycle_start(wcycle, ewcCONSTR);
 389
 390     if (econq == ConstraintVariable::ForceDispl && !EI_ENERGY_MINIMIZATION(ir.eI))
 391     {
 392         gmx_incons(
 393                 "constrain called for forces displacements while not doing energy minimization, "
 394                 "can not do this while the LINCS and SETTLE constraint connection matrices are "
 395                 "mass weighted");
 396     }
 397
 398     bOK   = TRUE;
 399     bDump = FALSE;
 400
 401     start = 0;
 402
 403     scaled_delta_t = step_scaling * ir.delta_t;
 404
 405     /* Prepare time step for use in constraint implementations, and
 406        avoid generating inf when ir.delta_t = 0. */
 407     if (ir.delta_t == 0)
 408     {
 409         invdt = 0.0;
 410     }
 411     else
 412     {
 413         invdt = 1.0 / scaled_delta_t;
 414     }
 415
 416     if (ir.efep != efepNO && EI_DYNAMICS(ir.eI))
 417     {
 418         /* Set the constraint lengths for the step at which this configuration
 419          * is meant to be. The invmasses should not be changed.
 420          */
 421         lambda += delta_step * ir.fepvals->delta_lambda;
 422     }
 423
 424     if (computeVirial)
 425     {
 426         clear_mat(constraintsVirial);
 427     }
 428     const InteractionList& settle = idef->il[F_SETTLE];
 429     nsettle                       = settle.size() / (1 + NRAL(F_SETTLE));
 430
 431     if (nsettle > 0)
 432     {
 433         nth = gmx_omp_nthreads_get(emntSETTLE);
 434     }
 435     else
 436     {
 437         nth = 1;
 438     }
 439
 440     /* We do not need full pbc when constraints do not cross update groups
 441      * i.e. when dd->constraint_comm==NULL.
 442      * Note that PBC for constraints is different from PBC for bondeds.
 443      * For constraints there is both forward and backward communication.
 444      */
 445     if (ir.pbcType != PbcType::No && (cr->dd || pbcHandlingRequired_)
 446         && !(cr->dd && cr->dd->constraint_comm == nullptr))
 447     {
 448         /* With pbc=screw the screw has been changed to a shift
 449          * by the constraint coordinate communication routine,
 450          * so that here we can use normal pbc.
 451          */
 452         pbc_null = set_pbc_dd(&pbc, ir.pbcType, DOMAINDECOMP(cr) ? cr->dd->numCells : nullptr, FALSE, box);
 453     }
 454     else
 455     {
 456         pbc_null = nullptr;
 457     }
 458
 459     /* Communicate the coordinates required for the non-local constraints
 460      * for LINCS and/or SETTLE.
 461      */
 462     if (cr->dd)
 463     {
 464         dd_move_x_constraints(cr->dd, box, x.unpaddedArrayRef(), xprime.unpaddedArrayRef(),
 465                               econq == ConstraintVariable::Positions);
 466
 467         if (!v.empty())
 468         {
 469             /* We need to initialize the non-local components of v.
 470              * We never actually use these values, but we do increment them,
 471              * so we should avoid uninitialized variables and overflows.
 472              */
 473             clear_constraint_quantity_nonlocal(cr->dd, v.unpaddedArrayRef());
 474         }
 475     }
 476
 477     if (lincsd != nullptr)
 478     {
 479         bOK = constrain_lincs(bLog || bEner, ir, step, lincsd, inverseMasses_, cr, ms, x, xprime,
 480                               min_proj, box, pbc_null, hasMassPerturbedAtoms_, lambda, dvdlambda,
 481                               invdt, v.unpaddedArrayRef(), computeVirial, constraintsVirial, econq,
 482                               nrnb, maxwarn, &warncount_lincs);
 483         if (!bOK && maxwarn < INT_MAX)
 484         {
 485             if (log != nullptr)
 486             {
 487                 fprintf(log, "Constraint error in algorithm %s at step %s\n",
 488                         econstr_names[econtLINCS], gmx_step_str(step, buf));
 489             }
 490             bDump = TRUE;
 491         }
 492     }
 493
 494     if (shaked != nullptr)
 495     {
 496         bOK = constrain_shake(log, shaked.get(), inverseMasses_, *idef, ir, x.unpaddedArrayRef(),
 497                               xprime.unpaddedArrayRef(), min_proj, pbc_null, nrnb, lambda,
 498                               dvdlambda, invdt, v.unpaddedArrayRef(), computeVirial,
 499                               constraintsVirial, maxwarn < INT_MAX, econq);
 500
 501         if (!bOK && maxwarn < INT_MAX)
 502         {
 503             if (log != nullptr)
 504             {
 505                 fprintf(log, "Constraint error in algorithm %s at step %s\n",
 506                         econstr_names[econtSHAKE], gmx_step_str(step, buf));
 507             }
 508             bDump = TRUE;
 509         }
 510     }
 511
 512     if (nsettle > 0)
 513     {
 514         bool bSettleErrorHasOccurred0 = false;
 515
 516         switch (econq)
 517         {
 518             case ConstraintVariable::Positions:
 519 #pragma omp parallel for num_threads(nth) schedule(static)
 520                 for (int th = 0; th < nth; th++)
 521                 {
 522                     try
 523                     {
 524                         if (th > 0)
 525                         {
 526                             clear_mat(threadConstraintsVirial[th]);
 527                         }
 528
 529                         csettle(*settled, nth, th, pbc_null, x, xprime, invdt, v, computeVirial,
 530                                 th == 0 ? constraintsVirial : threadConstraintsVirial[th],
 531                                 th == 0 ? &bSettleErrorHasOccurred0 : &bSettleErrorHasOccurred[th]);
 532                     }
 533                     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
 534                 }
 535                 inc_nrnb(nrnb, eNR_SETTLE, nsettle);
 536                 if (!v.empty())
 537                 {
 538                     inc_nrnb(nrnb, eNR_CONSTR_V, nsettle * 3);
 539                 }
 540                 if (computeVirial)
 541                 {
 542                     inc_nrnb(nrnb, eNR_CONSTR_VIR, nsettle * 3);
 543                 }
 544                 break;
 545             case ConstraintVariable::Velocities:
 546             case ConstraintVariable::Derivative:
 547             case ConstraintVariable::Force:
 548             case ConstraintVariable::ForceDispl:
 549 #pragma omp parallel for num_threads(nth) schedule(static)
 550                 for (int th = 0; th < nth; th++)
 551                 {
 552                     try
 553                     {
 554                         int calcvir_atom_end;
 555
 556                         if (!computeVirial)
 557                         {
 558                             calcvir_atom_end = 0;
 559                         }
 560                         else
 561                         {
 562                             calcvir_atom_end = numHomeAtoms_;
 563                         }
 564
 565                         if (th > 0)
 566                         {
 567                             clear_mat(threadConstraintsVirial[th]);
 568                         }
 569
 570                         int start_th = (nsettle * th) / nth;
 571                         int end_th   = (nsettle * (th + 1)) / nth;
 572
 573                         if (start_th >= 0 && end_th - start_th > 0)
 574                         {
 575                             settle_proj(*settled, econq, end_th - start_th,
 576                                         settle.iatoms.data() + start_th * (1 + NRAL(F_SETTLE)),
 577                                         pbc_null, x.unpaddedArrayRef(), xprime.unpaddedArrayRef(),
 578                                         min_proj, calcvir_atom_end,
 579                                         th == 0 ? constraintsVirial : threadConstraintsVirial[th]);
 580                         }
 581                     }
 582                     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
 583                 }
 584                 /* This is an overestimate */
 585                 inc_nrnb(nrnb, eNR_SETTLE, nsettle);
 586                 break;
 587             case ConstraintVariable::Deriv_FlexCon:
 588                 /* Nothing to do, since the are no flexible constraints in settles */
 589                 break;
 590             default: gmx_incons("Unknown constraint quantity for settle");
 591         }
 592
 593         if (computeVirial)
 594         {
 595             /* Reduce the virial contributions over the threads */
 596             for (int th = 1; th < nth; th++)
 597             {
 598                 m_add(constraintsVirial, threadConstraintsVirial[th], constraintsVirial);
 599             }
 600         }
 601
 602         if (econq == ConstraintVariable::Positions)
 603         {
 604             for (int th = 1; th < nth; th++)
 605             {
 606                 bSettleErrorHasOccurred0 = bSettleErrorHasOccurred0 || bSettleErrorHasOccurred[th];
 607             }
 608
 609             if (bSettleErrorHasOccurred0)
 610             {
 611                 char buf[STRLEN];
 612                 sprintf(buf,
 613                         "\nstep "
 614                         "%" PRId64
 615                         ": 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 (computeVirial)
 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: vir_fac = 0.5 / (ir.delta_t * ir.delta_t); break;
 648             case ConstraintVariable::Velocities: vir_fac = 0.5 / ir.delta_t; break;
 649             case ConstraintVariable::Force:
 650             case ConstraintVariable::ForceDispl: vir_fac = 0.5; break;
 651             default: gmx_incons("Unsupported constraint quantity for virial");
 652         }
 653
 654         if (EI_VV(ir.eI))
 655         {
 656             vir_fac *= 2; /* only constraining over half the distance here */
 657         }
 658         for (int i = 0; i < DIM; i++)
 659         {
 660             for (int j = 0; j < DIM; j++)
 661             {
 662                 constraintsVirial[i][j] *= vir_fac;
 663             }
 664         }
 665     }
 666
 667     if (bDump)
 668     {
 669         dump_confs(log, step, mtop, start, numHomeAtoms_, cr, x.unpaddedArrayRef(),
 670                    xprime.unpaddedArrayRef(), box);
 671     }
 672
 673     if (econq == ConstraintVariable::Positions)
 674     {
 675         if (ir.bPull && pull_have_constraint(*pull_work))
 676         {
 677             if (EI_DYNAMICS(ir.eI))
 678             {
 679                 t = ir.init_t + (step + delta_step) * ir.delta_t;
 680             }
 681             else
 682             {
 683                 t = ir.init_t;
 684             }
 685             set_pbc(&pbc, ir.pbcType, box);
 686             pull_constraint(pull_work, masses_, &pbc, cr, ir.delta_t, t,
 687                             as_rvec_array(x.unpaddedArrayRef().data()),
 688                             as_rvec_array(xprime.unpaddedArrayRef().data()),
 689                             as_rvec_array(v.unpaddedArrayRef().data()), constraintsVirial);
 690         }
 691         if (ed && delta_step > 0)
 692         {
 693             /* apply the essential dynamics constraints here */
 694             do_edsam(&ir, step, cr, as_rvec_array(xprime.unpaddedArrayRef().data()),
 695                      as_rvec_array(v.unpaddedArrayRef().data()), box, ed);
 696         }
 697     }
 698     wallcycle_stop(wcycle, ewcCONSTR);
 699
 700     if (!v.empty() && cFREEZE_)
 701     {
 702         /* Set the velocities of frozen dimensions to zero */
 703         ArrayRef<RVec> vRef = v.unpaddedArrayRef();
 704
 705         int gmx_unused numThreads = gmx_omp_nthreads_get(emntUpdate);
 706
 707 #pragma omp parallel for num_threads(numThreads) schedule(static)
 708         for (int i = 0; i < numHomeAtoms_; i++)
 709         {
 710             int freezeGroup = cFREEZE_[i];
 711
 712             for (int d = 0; d < DIM; d++)
 713             {
 714                 if (ir.opts.nFreeze[freezeGroup][d])
 715                 {
 716                     vRef[i][d] = 0;
 717                 }
 718             }
 719         }
 720     }
 721
 722     return bOK;
 723 } // namespace gmx
 724
 725 ArrayRef<real> Constraints::rmsdData() const
 726 {
 727     if (impl_->lincsd)
 728     {
 729         return lincs_rmsdData(impl_->lincsd);
 730     }
 731     else
 732     {
 733         return {};
 734     }
 735 }
 736
 737 real Constraints::rmsd() const
 738 {
 739     if (impl_->lincsd)
 740     {
 741         return lincs_rmsd(impl_->lincsd);
 742     }
 743     else
 744     {
 745         return 0;
 746     }
 747 }
 748
 749 int Constraints::numConstraintsTotal()
 750 {
 751     return impl_->ncon_tot;
 752 }
 753
 754 FlexibleConstraintTreatment flexibleConstraintTreatment(bool haveDynamicsIntegrator)
 755 {
 756     if (haveDynamicsIntegrator)
 757     {
 758         return FlexibleConstraintTreatment::Include;
 759     }
 760     else
 761     {
 762         return FlexibleConstraintTreatment::Exclude;
 763     }
 764 }
 765
 766 /*! \brief Returns a block struct to go from atoms to constraints
 767  *
 768  * The block struct will contain constraint indices with lower indices
 769  * directly matching the order in F_CONSTR and higher indices matching
 770  * the order in F_CONSTRNC offset by the number of constraints in F_CONSTR.
 771  *
 772  * \param[in]  numAtoms  The number of atoms to construct the list for
 773  * \param[in]  ilists    The interaction lists, size F_NRE
 774  * \param[in]  iparams   Interaction parameters, can be null when
 775  *                       \p flexibleConstraintTreatment==Include
 776  * \param[in]  flexibleConstraintTreatment  The flexible constraint treatment,
 777  *                                          see enum above
 778  *
 779  * \returns a block struct with all constraints for each atom
 780  */
 781 static ListOfLists<int> makeAtomsToConstraintsList(int                             numAtoms,
 782                                                    ArrayRef<const InteractionList> ilists,
 783                                                    ArrayRef<const t_iparams>       iparams,
 784                                                    FlexibleConstraintTreatment flexibleConstraintTreatment)
 785 {
 786     GMX_ASSERT(flexibleConstraintTreatment == FlexibleConstraintTreatment::Include || !iparams.empty(),
 787                "With flexible constraint detection we need valid iparams");
 788
 789     std::vector<int> count(numAtoms);
 790
 791     for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
 792     {
 793         const InteractionList& ilist  = ilists[ftype];
 794         const int              stride = 1 + NRAL(ftype);
 795         for (int i = 0; i < ilist.size(); i += stride)
 796         {
 797             if (flexibleConstraintTreatment == FlexibleConstraintTreatment::Include
 798                 || !isConstraintFlexible(iparams, ilist.iatoms[i]))
 799             {
 800                 for (int j = 1; j < 3; j++)
 801                 {
 802                     int a = ilist.iatoms[i + j];
 803                     count[a]++;
 804                 }
 805             }
 806         }
 807     }
 808
 809     std::vector<int> listRanges(numAtoms + 1);
 810     for (int a = 0; a < numAtoms; a++)
 811     {
 812         listRanges[a + 1] = listRanges[a] + count[a];
 813         count[a]          = 0;
 814     }
 815     std::vector<int> elements(listRanges[numAtoms]);
 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 InteractionList& 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                     const int a                          = ilist.iatoms[i + j];
 833                     elements[listRanges[a] + count[a]++] = numConstraints;
 834                 }
 835             }
 836             numConstraints++;
 837         }
 838     }
 839
 840     return ListOfLists<int>(std::move(listRanges), std::move(elements));
 841 }
 842
 843 ListOfLists<int> make_at2con(int                             numAtoms,
 844                              ArrayRef<const InteractionList> ilist,
 845                              ArrayRef<const t_iparams>       iparams,
 846                              FlexibleConstraintTreatment     flexibleConstraintTreatment)
 847 {
 848     return makeAtomsToConstraintsList(numAtoms, ilist, iparams, flexibleConstraintTreatment);
 849 }
 850
 851 ListOfLists<int> make_at2con(const gmx_moltype_t&           moltype,
 852                              gmx::ArrayRef<const t_iparams> iparams,
 853                              FlexibleConstraintTreatment    flexibleConstraintTreatment)
 854 {
 855     return makeAtomsToConstraintsList(moltype.atoms.nr, makeConstArrayRef(moltype.ilist), iparams,
 856                                       flexibleConstraintTreatment);
 857 }
 858
 859 //! Return the number of flexible constraints in the \c ilist and \c iparams.
 860 int countFlexibleConstraints(ArrayRef<const InteractionList> ilist, ArrayRef<const t_iparams> iparams)
 861 {
 862     int nflexcon = 0;
 863     for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
 864     {
 865         const int numIatomsPerConstraint = 3;
 866         for (int i = 0; i < ilist[ftype].size(); i += numIatomsPerConstraint)
 867         {
 868             const int type = ilist[ftype].iatoms[i];
 869             if (iparams[type].constr.dA == 0 && iparams[type].constr.dB == 0)
 870             {
 871                 nflexcon++;
 872             }
 873         }
 874     }
 875
 876     return nflexcon;
 877 }
 878
 879 //! Returns the index of the settle to which each atom belongs.
 880 static std::vector<int> make_at2settle(int natoms, const InteractionList& ilist)
 881 {
 882     /* Set all to no settle */
 883     std::vector<int> at2s(natoms, -1);
 884
 885     const int stride = 1 + NRAL(F_SETTLE);
 886
 887     for (int s = 0; s < ilist.size(); s += stride)
 888     {
 889         at2s[ilist.iatoms[s + 1]] = s / stride;
 890         at2s[ilist.iatoms[s + 2]] = s / stride;
 891         at2s[ilist.iatoms[s + 3]] = s / stride;
 892     }
 893
 894     return at2s;
 895 }
 896
 897 void Constraints::Impl::setConstraints(gmx_localtop_t* top,
 898                                        int             numAtoms,
 899                                        int             numHomeAtoms,
 900                                        real*           masses,
 901                                        real*           inverseMasses,
 902                                        const bool      hasMassPerturbedAtoms,
 903                                        const real      lambda,
 904                                        unsigned short* cFREEZE)
 905 {
 906     numAtoms_              = numAtoms;
 907     numHomeAtoms_          = numHomeAtoms;
 908     masses_                = masses;
 909     inverseMasses_         = inverseMasses;
 910     hasMassPerturbedAtoms_ = hasMassPerturbedAtoms;
 911     lambda_                = lambda;
 912     cFREEZE_               = cFREEZE;
 913
 914     idef = &top->idef;
 915
 916     if (ncon_tot > 0)
 917     {
 918         /* With DD we might also need to call LINCS on a domain no constraints for
 919          * communicating coordinates to other nodes that do have constraints.
 920          */
 921         if (ir.eConstrAlg == econtLINCS)
 922         {
 923             set_lincs(*idef, numAtoms_, inverseMasses_, lambda_, EI_DYNAMICS(ir.eI), cr, lincsd);
 924         }
 925         if (ir.eConstrAlg == econtSHAKE)
 926         {
 927             if (cr->dd)
 928             {
 929                 // We are using the local topology, so there are only
 930                 // F_CONSTR constraints.
 931                 GMX_RELEASE_ASSERT(idef->il[F_CONSTRNC].empty(),
 932                                    "Here we should not have no-connect constraints");
 933                 make_shake_sblock_dd(shaked.get(), idef->il[F_CONSTR]);
 934             }
 935             else
 936             {
 937                 make_shake_sblock_serial(shaked.get(), &top->idef, numAtoms_);
 938             }
 939         }
 940     }
 941
 942     if (settled)
 943     {
 944         settled->setConstraints(idef->il[F_SETTLE], numHomeAtoms_, masses_, inverseMasses_);
 945     }
 946
 947     /* Make a selection of the local atoms for essential dynamics */
 948     if (ed && cr->dd)
 949     {
 950         dd_make_local_ed_indices(cr->dd, ed);
 951     }
 952 }
 953
 954 void Constraints::setConstraints(gmx_localtop_t* top,
 955                                  const int       numAtoms,
 956                                  const int       numHomeAtoms,
 957                                  real*           masses,
 958                                  real*           inverseMasses,
 959                                  const bool      hasMassPerturbedAtoms,
 960                                  const real      lambda,
 961                                  unsigned short* cFREEZE)
 962 {
 963     impl_->setConstraints(top, numAtoms, numHomeAtoms, masses, inverseMasses, hasMassPerturbedAtoms,
 964                           lambda, cFREEZE);
 965 }
 966
 967 /*! \brief Makes a per-moleculetype container of mappings from atom
 968  * indices to constraint indices.
 969  *
 970  * Note that flexible constraints are only enabled with a dynamical integrator. */
 971 static std::vector<ListOfLists<int>> makeAtomToConstraintMappings(const gmx_mtop_t& mtop,
 972                                                                   FlexibleConstraintTreatment flexibleConstraintTreatment)
 973 {
 974     std::vector<ListOfLists<int>> mapping;
 975     mapping.reserve(mtop.moltype.size());
 976     for (const gmx_moltype_t& moltype : mtop.moltype)
 977     {
 978         mapping.push_back(make_at2con(moltype, mtop.ffparams.iparams, flexibleConstraintTreatment));
 979     }
 980     return mapping;
 981 }
 982
 983 Constraints::Constraints(const gmx_mtop_t&     mtop,
 984                          const t_inputrec&     ir,
 985                          pull_t*               pull_work,
 986                          FILE*                 log,
 987                          const t_commrec*      cr,
 988                          const gmx_multisim_t* ms,
 989                          t_nrnb*               nrnb,
 990                          gmx_wallcycle*        wcycle,
 991                          bool                  pbcHandlingRequired,
 992                          int                   numConstraints,
 993                          int                   numSettles) :
 994     impl_(new Impl(mtop, ir, pull_work, log, cr, ms, nrnb, wcycle, pbcHandlingRequired, numConstraints, numSettles))
 995 {
 996 }
 997
 998 Constraints::Impl::Impl(const gmx_mtop_t&     mtop_p,
 999                         const t_inputrec&     ir_p,
1000                         pull_t*               pull_work,
1001                         FILE*                 log_p,
1002                         const t_commrec*      cr_p,
1003                         const gmx_multisim_t* ms_p,
1004                         t_nrnb*               nrnb_p,
1005                         gmx_wallcycle*        wcycle_p,
1006                         bool                  pbcHandlingRequired,
1007                         int                   numConstraints,
1008                         int                   numSettles) :
1009     ncon_tot(numConstraints),
1010     mtop(mtop_p),
1011     pbcHandlingRequired_(pbcHandlingRequired),
1012     log(log_p),
1013     cr(cr_p),
1014     ms(ms_p),
1015     pull_work(pull_work),
1016     ir(ir_p),
1017     nrnb(nrnb_p),
1018     wcycle(wcycle_p)
1019 {
1020     if (numConstraints + numSettles > 0 && ir.epc == epcMTTK)
1021     {
1022         gmx_fatal(FARGS, "Constraints are not implemented with MTTK pressure control.");
1023     }
1024
1025     nflexcon = 0;
1026     if (numConstraints > 0)
1027     {
1028         at2con_mt = makeAtomToConstraintMappings(mtop, flexibleConstraintTreatment(EI_DYNAMICS(ir.eI)));
1029
1030         for (const gmx_molblock_t& molblock : mtop.molblock)
1031         {
1032             int count = countFlexibleConstraints(mtop.moltype[molblock.type].ilist, mtop.ffparams.iparams);
1033             nflexcon += molblock.nmol * count;
1034         }
1035
1036         if (nflexcon > 0)
1037         {
1038             if (log)
1039             {
1040                 fprintf(log, "There are %d flexible constraints\n", nflexcon);
1041                 if (ir.fc_stepsize == 0)
1042                 {
1043                     fprintf(log,
1044                             "\n"
1045                             "WARNING: step size for flexible constraining = 0\n"
1046                             "         All flexible constraints will be rigid.\n"
1047                             "         Will try to keep all flexible constraints at their original "
1048                             "length,\n"
1049                             "         but the lengths may exhibit some drift.\n\n");
1050                     nflexcon = 0;
1051                 }
1052             }
1053             if (nflexcon > 0)
1054             {
1055                 please_cite(log, "Hess2002");
1056             }
1057         }
1058
1059         if (ir.eConstrAlg == econtLINCS)
1060         {
1061             lincsd = init_lincs(log, mtop, nflexcon, at2con_mt,
1062                                 DOMAINDECOMP(cr) && ddHaveSplitConstraints(*cr->dd), ir.nLincsIter,
1063                                 ir.nProjOrder);
1064         }
1065
1066         if (ir.eConstrAlg == econtSHAKE)
1067         {
1068             if (DOMAINDECOMP(cr) && ddHaveSplitConstraints(*cr->dd))
1069             {
1070                 gmx_fatal(FARGS,
1071                           "SHAKE is not supported with domain decomposition and constraint that "
1072                           "cross domain boundaries, use LINCS");
1073             }
1074             if (nflexcon)
1075             {
1076                 gmx_fatal(FARGS,
1077                           "For this system also velocities and/or forces need to be constrained, "
1078                           "this can not be done with SHAKE, you should select LINCS");
1079             }
1080             please_cite(log, "Ryckaert77a");
1081             if (ir.bShakeSOR)
1082             {
1083                 please_cite(log, "Barth95a");
1084             }
1085
1086             shaked = std::make_unique<shakedata>();
1087         }
1088     }
1089
1090     if (numSettles > 0)
1091     {
1092         please_cite(log, "Miyamoto92a");
1093
1094         settled = std::make_unique<SettleData>(mtop);
1095
1096         /* Make an atom to settle index for use in domain decomposition */
1097         for (size_t mt = 0; mt < mtop.moltype.size(); mt++)
1098         {
1099             at2settle_mt.emplace_back(
1100                     make_at2settle(mtop.moltype[mt].atoms.nr, mtop.moltype[mt].ilist[F_SETTLE]));
1101         }
1102
1103         /* Allocate thread-local work arrays */
1104         int nthreads = gmx_omp_nthreads_get(emntSETTLE);
1105         if (nthreads > 1 && threadConstraintsVirial == nullptr)
1106         {
1107             snew(threadConstraintsVirial, nthreads);
1108             snew(bSettleErrorHasOccurred, nthreads);
1109         }
1110     }
1111
1112     maxwarn   = 999;
1113     char* env = getenv("GMX_MAXCONSTRWARN");
1114     if (env)
1115     {
1116         maxwarn = 0;
1117         sscanf(env, "%8d", &maxwarn);
1118         if (maxwarn < 0)
1119         {
1120             maxwarn = INT_MAX;
1121         }
1122         if (log)
1123         {
1124             fprintf(log, "Setting the maximum number of constraint warnings to %d\n", maxwarn);
1125         }
1126         if (MASTER(cr))
1127         {
1128             fprintf(stderr, "Setting the maximum number of constraint warnings to %d\n", maxwarn);
1129         }
1130     }
1131     warncount_lincs  = 0;
1132     warncount_settle = 0;
1133 }
1134
1135 Constraints::Impl::~Impl()
1136 {
1137     if (bSettleErrorHasOccurred != nullptr)
1138     {
1139         sfree(bSettleErrorHasOccurred);
1140     }
1141     if (threadConstraintsVirial != nullptr)
1142     {
1143         sfree(threadConstraintsVirial);
1144     }
1145     done_lincs(lincsd);
1146 }
1147
1148 void Constraints::saveEdsamPointer(gmx_edsam* ed)
1149 {
1150     impl_->ed = ed;
1151 }
1152
1153 ArrayRef<const ListOfLists<int>> Constraints::atom2constraints_moltype() const
1154 {
1155     return impl_->at2con_mt;
1156 }
1157
1158 ArrayRef<const std::vector<int>> Constraints::atom2settle_moltype() const
1159 {
1160     return impl_->at2settle_mt;
1161 }
1162
1163 void do_constrain_first(FILE*                     fplog,
1164                         gmx::Constraints*         constr,
1165                         const t_inputrec*         ir,
1166                         int                       numAtoms,
1167                         int                       numHomeAtoms,
1168                         ArrayRefWithPadding<RVec> x,
1169                         ArrayRefWithPadding<RVec> v,
1170                         const matrix              box,
1171                         real                      lambda)
1172 {
1173     int     i, m, start, end;
1174     int64_t step;
1175     real    dt = ir->delta_t;
1176     real    dvdl_dum;
1177
1178     PaddedVector<RVec> savex(numAtoms);
1179
1180     start = 0;
1181     end   = numHomeAtoms;
1182
1183     if (debug)
1184     {
1185         fprintf(debug, "vcm: start=%d, homenr=%d, end=%d\n", start, numHomeAtoms, end);
1186     }
1187     /* Do a first constrain to reset particles... */
1188     step = ir->init_step;
1189     if (fplog)
1190     {
1191         char buf[STEPSTRSIZE];
1192         fprintf(fplog, "\nConstraining the starting coordinates (step %s)\n", gmx_step_str(step, buf));
1193     }
1194     dvdl_dum = 0;
1195
1196     bool needsLogging  = true;
1197     bool computeEnergy = false;
1198     bool computeVirial = false;
1199     /* constrain the current position */
1200     constr->apply(needsLogging, computeEnergy, step, 0, 1.0, x, x, {}, box, lambda, &dvdl_dum, {},
1201                   computeVirial, nullptr, gmx::ConstraintVariable::Positions);
1202     if (EI_VV(ir->eI))
1203     {
1204         /* constrain the inital velocity, and save it */
1205         /* also may be useful if we need the ekin from the halfstep for velocity verlet */
1206         constr->apply(needsLogging, computeEnergy, step, 0, 1.0, x, v, v.unpaddedArrayRef(), box, lambda,
1207                       &dvdl_dum, {}, computeVirial, nullptr, gmx::ConstraintVariable::Velocities);
1208     }
1209     /* constrain the inital velocities at t-dt/2 */
1210     if (EI_STATE_VELOCITY(ir->eI) && ir->eI != eiVV)
1211     {
1212         auto subX = x.paddedArrayRef().subArray(start, end);
1213         auto subV = v.paddedArrayRef().subArray(start, end);
1214         for (i = start; (i < end); i++)
1215         {
1216             for (m = 0; (m < DIM); m++)
1217             {
1218                 /* Reverse the velocity */
1219                 subV[i][m] = -subV[i][m];
1220                 /* Store the position at t-dt in buf */
1221                 savex[i][m] = subX[i][m] + dt * subV[i][m];
1222             }
1223         }
1224         /* Shake the positions at t=-dt with the positions at t=0
1225          * as reference coordinates.
1226          */
1227         if (fplog)
1228         {
1229             char buf[STEPSTRSIZE];
1230             fprintf(fplog, "\nConstraining the coordinates at t0-dt (step %s)\n", gmx_step_str(step, buf));
1231         }
1232         dvdl_dum = 0;
1233         constr->apply(needsLogging, computeEnergy, step, -1, 1.0, x, savex.arrayRefWithPadding(),
1234                       {}, box, lambda, &dvdl_dum, v, computeVirial, nullptr,
1235                       gmx::ConstraintVariable::Positions);
1236
1237         for (i = start; i < end; i++)
1238         {
1239             for (m = 0; m < DIM; m++)
1240             {
1241                 /* Re-reverse the velocities */
1242                 subV[i][m] = -subV[i][m];
1243             }
1244         }
1245     }
1246 }
1247
1248 void constrain_velocities(gmx::Constraints* constr,
1249                           bool              do_log,
1250                           bool              do_ene,
1251                           int64_t           step,
1252                           t_state*          state,
1253                           real*             dvdlambda,
1254                           gmx_bool          computeVirial,
1255                           tensor            constraintsVirial)
1256 {
1257     if (constr != nullptr)
1258     {
1259         constr->apply(do_log, do_ene, step, 1, 1.0, state->x.arrayRefWithPadding(),
1260                       state->v.arrayRefWithPadding(), state->v.arrayRefWithPadding().unpaddedArrayRef(),
1261                       state->box, state->lambda[efptBONDED], dvdlambda, ArrayRefWithPadding<RVec>(),
1262                       computeVirial, constraintsVirial, ConstraintVariable::Velocities);
1263     }
1264 }
1265
1266 void constrain_coordinates(gmx::Constraints*         constr,
1267                            bool                      do_log,
1268                            bool                      do_ene,
1269                            int64_t                   step,
1270                            t_state*                  state,
1271                            ArrayRefWithPadding<RVec> xp,
1272                            real*                     dvdlambda,
1273                            gmx_bool                  computeVirial,
1274                            tensor                    constraintsVirial)
1275 {
1276     if (constr != nullptr)
1277     {
1278         constr->apply(do_log, do_ene, step, 1, 1.0, state->x.arrayRefWithPadding(), std::move(xp),
1279                       ArrayRef<RVec>(), state->box, state->lambda[efptBONDED], dvdlambda,
1280                       state->v.arrayRefWithPadding(), computeVirial, constraintsVirial,
1281                       ConstraintVariable::Positions);
1282     }
1283 }
1284
1285 } // namespace gmx