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