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