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