src/gromacs/mdlib/constr.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 /*! \internal \file
  39  * \brief Defines the high-level constraint code.
  40  *
  41  * \author Berk Hess <hess@kth.se>
  42  * \author Mark Abraham <mark.j.abraham@gmail.com>
  43  * \ingroup module_mdlib
  44  */
  45 #include "gmxpre.h"
  46
  47 #include "constr.h"
  48
  49 #include <cassert>
  50 #include <cmath>
  51 #include <cstdlib>
  52
  53 #include <algorithm>
  54
  55 #include "gromacs/domdec/domdec.h"
  56 #include "gromacs/domdec/domdec_struct.h"
  57 #include "gromacs/essentialdynamics/edsam.h"
  58 #include "gromacs/fileio/confio.h"
  59 #include "gromacs/fileio/gmxfio.h"
  60 #include "gromacs/fileio/pdbio.h"
  61 #include "gromacs/gmxlib/nrnb.h"
  62 #include "gromacs/math/utilities.h"
  63 #include "gromacs/math/vec.h"
  64 #include "gromacs/mdlib/gmx_omp_nthreads.h"
  65 #include "gromacs/mdlib/lincs.h"
  66 #include "gromacs/mdlib/settle.h"
  67 #include "gromacs/mdlib/shake.h"
  68 #include "gromacs/mdtypes/commrec.h"
  69 #include "gromacs/mdtypes/inputrec.h"
  70 #include "gromacs/mdtypes/md_enums.h"
  71 #include "gromacs/mdtypes/mdatom.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_mdatoms&      md_p,
 107          const t_commrec*      cr_p,
 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, const t_mdatoms& md);
 116     bool apply(bool                      bLog,
 117                bool                      bEner,
 118                int64_t                   step,
 119                int                       delta_step,
 120                real                      step_scaling,
 121                ArrayRefWithPadding<RVec> x,
 122                ArrayRefWithPadding<RVec> xprime,
 123                ArrayRef<RVec>            min_proj,
 124                const matrix              box,
 125                real                      lambda,
 126                real*                     dvdlambda,
 127                ArrayRefWithPadding<RVec> v,
 128                tensor*                   vir,
 129                ConstraintVariable        econq);
 130     //! The total number of constraints.
 131     int ncon_tot = 0;
 132     //! The number of flexible constraints.
 133     int nflexcon = 0;
 134     //! A list of atoms to constraints for each moleculetype.
 135     std::vector<ListOfLists<int>> at2con_mt;
 136     //! A list of atoms to settles for each moleculetype
 137     std::vector<std::vector<int>> at2settle_mt;
 138     //! LINCS data.
 139     Lincs* lincsd = nullptr; // TODO this should become a unique_ptr
 140     //! SHAKE data.
 141     shakedata* shaked = nullptr;
 142     //! SETTLE data.
 143     settledata* settled = nullptr;
 144     //! The maximum number of warnings.
 145     int maxwarn = 0;
 146     //! The number of warnings for LINCS.
 147     int warncount_lincs = 0;
 148     //! The number of warnings for SETTLE.
 149     int warncount_settle = 0;
 150     //! The essential dynamics data.
 151     gmx_edsam* ed = nullptr;
 152
 153     //! Thread-local virial contribution.
 154     tensor* vir_r_m_dr_th = { nullptr };
 155     //! Did a settle error occur?
 156     bool* bSettleErrorHasOccurred = nullptr;
 157
 158     //! Pointer to the global topology - only used for printing warnings.
 159     const gmx_mtop_t& mtop;
 160     //! Parameters for the interactions in this domain.
 161     const InteractionDefinitions* idef = nullptr;
 162     //! Data about atoms in this domain.
 163     const t_mdatoms& md;
 164     //! Whether we need to do pbc for handling bonds.
 165     bool pbcHandlingRequired_ = false;
 166
 167     //! Logging support.
 168     FILE* log = nullptr;
 169     //! Communication support.
 170     const t_commrec* cr = nullptr;
 171     //! Multi-sim support.
 172     const gmx_multisim_t* ms = nullptr;
 173     //! Pulling code object, if any.
 174     pull_t* pull_work = nullptr;
 175     /*!\brief Input options.
 176      *
 177      * \todo Replace with IMdpOptions */
 178     const t_inputrec& ir;
 179     //! Flop counting support.
 180     t_nrnb* nrnb = nullptr;
 181     //! Tracks wallcycle usage.
 182     gmx_wallcycle* wcycle;
 183 };
 184
 185 Constraints::~Constraints() = default;
 186
 187 int Constraints::numFlexibleConstraints() const
 188 {
 189     return impl_->nflexcon;
 190 }
 191
 192 bool Constraints::havePerturbedConstraints() const
 193 {
 194     const gmx_ffparams_t& ffparams = impl_->mtop.ffparams;
 195
 196     for (size_t i = 0; i < ffparams.functype.size(); i++)
 197     {
 198         if ((ffparams.functype[i] == F_CONSTR || ffparams.functype[i] == F_CONSTRNC)
 199             && ffparams.iparams[i].constr.dA != ffparams.iparams[i].constr.dB)
 200         {
 201             return true;
 202         }
 203     }
 204
 205     return false;
 206 }
 207
 208 //! Clears constraint quantities for atoms in nonlocal region.
 209 static void clear_constraint_quantity_nonlocal(gmx_domdec_t* dd, ArrayRef<RVec> q)
 210 {
 211     int nonlocal_at_start, nonlocal_at_end, at;
 212
 213     dd_get_constraint_range(dd, &nonlocal_at_start, &nonlocal_at_end);
 214
 215     for (at = nonlocal_at_start; at < nonlocal_at_end; at++)
 216     {
 217         clear_rvec(q[at]);
 218     }
 219 }
 220
 221 void too_many_constraint_warnings(int eConstrAlg, int warncount)
 222 {
 223     gmx_fatal(
 224             FARGS,
 225             "Too many %s warnings (%d)\n"
 226             "If you know what you are doing you can %s"
 227             "set the environment variable GMX_MAXCONSTRWARN to -1,\n"
 228             "but normally it is better to fix the problem",
 229             (eConstrAlg == econtLINCS) ? "LINCS" : "SETTLE", warncount,
 230             (eConstrAlg == econtLINCS) ? "adjust the lincs warning threshold in your mdp file\nor " : "\n");
 231 }
 232
 233 //! Writes out coordinates.
 234 static void write_constr_pdb(const char*          fn,
 235                              const char*          title,
 236                              const gmx_mtop_t&    mtop,
 237                              int                  start,
 238                              int                  homenr,
 239                              const t_commrec*     cr,
 240                              ArrayRef<const RVec> x,
 241                              const matrix         box)
 242 {
 243     char          fname[STRLEN];
 244     FILE*         out;
 245     int           dd_ac0 = 0, dd_ac1 = 0, i, ii, resnr;
 246     gmx_domdec_t* dd;
 247     const char *  anm, *resnm;
 248
 249     dd = nullptr;
 250     if (DOMAINDECOMP(cr))
 251     {
 252         dd = cr->dd;
 253         dd_get_constraint_range(dd, &dd_ac0, &dd_ac1);
 254         start  = 0;
 255         homenr = dd_ac1;
 256     }
 257
 258     if (PAR(cr))
 259     {
 260         sprintf(fname, "%s_n%d.pdb", fn, cr->sim_nodeid);
 261     }
 262     else
 263     {
 264         sprintf(fname, "%s.pdb", fn);
 265     }
 266
 267     out = gmx_fio_fopen(fname, "w");
 268
 269     fprintf(out, "TITLE     %s\n", title);
 270     gmx_write_pdb_box(out, PbcType::Unset, box);
 271     int molb = 0;
 272     for (i = start; i < start + homenr; i++)
 273     {
 274         if (dd != nullptr)
 275         {
 276             if (i >= dd_numHomeAtoms(*dd) && i < dd_ac0)
 277             {
 278                 continue;
 279             }
 280             ii = dd->globalAtomIndices[i];
 281         }
 282         else
 283         {
 284             ii = i;
 285         }
 286         mtopGetAtomAndResidueName(mtop, ii, &molb, &anm, &resnr, &resnm, nullptr);
 287         gmx_fprintf_pdb_atomline(out, epdbATOM, ii + 1, anm, ' ', resnm, ' ', resnr, ' ',
 288                                  10 * x[i][XX], 10 * x[i][YY], 10 * x[i][ZZ], 1.0, 0.0, "");
 289     }
 290     fprintf(out, "TER\n");
 291
 292     gmx_fio_fclose(out);
 293 }
 294
 295 //! Writes out domain contents to help diagnose crashes.
 296 static void dump_confs(FILE*                log,
 297                        int64_t              step,
 298                        const gmx_mtop_t&    mtop,
 299                        int                  start,
 300                        int                  homenr,
 301                        const t_commrec*     cr,
 302                        ArrayRef<const RVec> x,
 303                        ArrayRef<const RVec> xprime,
 304                        const matrix         box)
 305 {
 306     char buf[STRLEN], buf2[22];
 307
 308     char* env = getenv("GMX_SUPPRESS_DUMP");
 309     if (env)
 310     {
 311         return;
 312     }
 313
 314     sprintf(buf, "step%sb", gmx_step_str(step, buf2));
 315     write_constr_pdb(buf, "initial coordinates", mtop, start, homenr, cr, x, box);
 316     sprintf(buf, "step%sc", gmx_step_str(step, buf2));
 317     write_constr_pdb(buf, "coordinates after constraining", mtop, start, homenr, cr, xprime, box);
 318     if (log)
 319     {
 320         fprintf(log, "Wrote pdb files with previous and current coordinates\n");
 321     }
 322     fprintf(stderr, "Wrote pdb files with previous and current coordinates\n");
 323 }
 324
 325 bool Constraints::apply(bool                      bLog,
 326                         bool                      bEner,
 327                         int64_t                   step,
 328                         int                       delta_step,
 329                         real                      step_scaling,
 330                         ArrayRefWithPadding<RVec> x,
 331                         ArrayRefWithPadding<RVec> xprime,
 332                         ArrayRef<RVec>            min_proj,
 333                         const matrix              box,
 334                         real                      lambda,
 335                         real*                     dvdlambda,
 336                         ArrayRefWithPadding<RVec> v,
 337                         tensor*                   vir,
 338                         ConstraintVariable        econq)
 339 {
 340     return impl_->apply(bLog, bEner, step, delta_step, step_scaling, std::move(x), std::move(xprime),
 341                         min_proj, box, lambda, dvdlambda, std::move(v), vir, econq);
 342 }
 343
 344 bool Constraints::Impl::apply(bool                      bLog,
 345                               bool                      bEner,
 346                               int64_t                   step,
 347                               int                       delta_step,
 348                               real                      step_scaling,
 349                               ArrayRefWithPadding<RVec> x,
 350                               ArrayRefWithPadding<RVec> xprime,
 351                               ArrayRef<RVec>            min_proj,
 352                               const matrix              box,
 353                               real                      lambda,
 354                               real*                     dvdlambda,
 355                               ArrayRefWithPadding<RVec> v,
 356                               tensor*                   vir,
 357                               ConstraintVariable        econq)
 358 {
 359     bool   bOK, bDump;
 360     int    start, homenr;
 361     tensor vir_r_m_dr;
 362     real   scaled_delta_t;
 363     real   invdt, vir_fac = 0, t;
 364     int    nsettle;
 365     t_pbc  pbc, *pbc_null;
 366     char   buf[22];
 367     int    nth;
 368
 369     wallcycle_start(wcycle, ewcCONSTR);
 370
 371     if (econq == ConstraintVariable::ForceDispl && !EI_ENERGY_MINIMIZATION(ir.eI))
 372     {
 373         gmx_incons(
 374                 "constrain called for forces displacements while not doing energy minimization, "
 375                 "can not do this while the LINCS and SETTLE constraint connection matrices are "
 376                 "mass weighted");
 377     }
 378
 379     bOK   = TRUE;
 380     bDump = FALSE;
 381
 382     start  = 0;
 383     homenr = md.homenr;
 384
 385     scaled_delta_t = step_scaling * ir.delta_t;
 386
 387     /* Prepare time step for use in constraint implementations, and
 388        avoid generating inf when ir.delta_t = 0. */
 389     if (ir.delta_t == 0)
 390     {
 391         invdt = 0.0;
 392     }
 393     else
 394     {
 395         invdt = 1.0 / scaled_delta_t;
 396     }
 397
 398     if (ir.efep != efepNO && EI_DYNAMICS(ir.eI))
 399     {
 400         /* Set the constraint lengths for the step at which this configuration
 401          * is meant to be. The invmasses should not be changed.
 402          */
 403         lambda += delta_step * ir.fepvals->delta_lambda;
 404     }
 405
 406     if (vir != nullptr)
 407     {
 408         clear_mat(vir_r_m_dr);
 409     }
 410     const InteractionList& settle = idef->il[F_SETTLE];
 411     nsettle                       = settle.size() / (1 + NRAL(F_SETTLE));
 412
 413     if (nsettle > 0)
 414     {
 415         nth = gmx_omp_nthreads_get(emntSETTLE);
 416     }
 417     else
 418     {
 419         nth = 1;
 420     }
 421
 422     /* We do not need full pbc when constraints do not cross update groups
 423      * i.e. when dd->constraint_comm==NULL.
 424      * Note that PBC for constraints is different from PBC for bondeds.
 425      * For constraints there is both forward and backward communication.
 426      */
 427     if (ir.pbcType != PbcType::No && (cr->dd || pbcHandlingRequired_)
 428         && !(cr->dd && cr->dd->constraint_comm == nullptr))
 429     {
 430         /* With pbc=screw the screw has been changed to a shift
 431          * by the constraint coordinate communication routine,
 432          * so that here we can use normal pbc.
 433          */
 434         pbc_null = set_pbc_dd(&pbc, ir.pbcType, DOMAINDECOMP(cr) ? cr->dd->numCells : nullptr, FALSE, box);
 435     }
 436     else
 437     {
 438         pbc_null = nullptr;
 439     }
 440
 441     /* Communicate the coordinates required for the non-local constraints
 442      * for LINCS and/or SETTLE.
 443      */
 444     if (cr->dd)
 445     {
 446         dd_move_x_constraints(cr->dd, box, x.unpaddedArrayRef(), xprime.unpaddedArrayRef(),
 447                               econq == ConstraintVariable::Positions);
 448
 449         if (!v.empty())
 450         {
 451             /* We need to initialize the non-local components of v.
 452              * We never actually use these values, but we do increment them,
 453              * so we should avoid uninitialized variables and overflows.
 454              */
 455             clear_constraint_quantity_nonlocal(cr->dd, v.unpaddedArrayRef());
 456         }
 457     }
 458
 459     if (lincsd != nullptr)
 460     {
 461         bOK = constrain_lincs(bLog || bEner, ir, step, lincsd, md, cr, ms, x, xprime, min_proj, box,
 462                               pbc_null, lambda, dvdlambda, invdt, v.unpaddedArrayRef(),
 463                               vir != nullptr, vir_r_m_dr, econq, nrnb, maxwarn, &warncount_lincs);
 464         if (!bOK && maxwarn < INT_MAX)
 465         {
 466             if (log != nullptr)
 467             {
 468                 fprintf(log, "Constraint error in algorithm %s at step %s\n",
 469                         econstr_names[econtLINCS], gmx_step_str(step, buf));
 470             }
 471             bDump = TRUE;
 472         }
 473     }
 474
 475     if (shaked != nullptr)
 476     {
 477         bOK = constrain_shake(
 478                 log, shaked, md.invmass, *idef, ir, as_rvec_array(x.unpaddedArrayRef().data()),
 479                 as_rvec_array(xprime.unpaddedArrayRef().data()), as_rvec_array(min_proj.data()),
 480                 nrnb, lambda, dvdlambda, invdt, as_rvec_array(v.unpaddedArrayRef().data()),
 481                 vir != nullptr, vir_r_m_dr, maxwarn < INT_MAX, econq);
 482
 483         if (!bOK && maxwarn < INT_MAX)
 484         {
 485             if (log != nullptr)
 486             {
 487                 fprintf(log, "Constraint error in algorithm %s at step %s\n",
 488                         econstr_names[econtSHAKE], gmx_step_str(step, buf));
 489             }
 490             bDump = TRUE;
 491         }
 492     }
 493
 494     if (nsettle > 0)
 495     {
 496         bool bSettleErrorHasOccurred0 = false;
 497
 498         switch (econq)
 499         {
 500             case ConstraintVariable::Positions:
 501 #pragma omp parallel for num_threads(nth) schedule(static)
 502                 for (int th = 0; th < nth; th++)
 503                 {
 504                     try
 505                     {
 506                         if (th > 0)
 507                         {
 508                             clear_mat(vir_r_m_dr_th[th]);
 509                         }
 510
 511                         csettle(settled, nth, th, pbc_null, x, xprime, invdt, v, vir != nullptr,
 512                                 th == 0 ? vir_r_m_dr : vir_r_m_dr_th[th],
 513                                 th == 0 ? &bSettleErrorHasOccurred0 : &bSettleErrorHasOccurred[th]);
 514                     }
 515                     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
 516                 }
 517                 inc_nrnb(nrnb, eNR_SETTLE, nsettle);
 518                 if (!v.empty())
 519                 {
 520                     inc_nrnb(nrnb, eNR_CONSTR_V, nsettle * 3);
 521                 }
 522                 if (vir != nullptr)
 523                 {
 524                     inc_nrnb(nrnb, eNR_CONSTR_VIR, nsettle * 3);
 525                 }
 526                 break;
 527             case ConstraintVariable::Velocities:
 528             case ConstraintVariable::Derivative:
 529             case ConstraintVariable::Force:
 530             case ConstraintVariable::ForceDispl:
 531 #pragma omp parallel for num_threads(nth) schedule(static)
 532                 for (int th = 0; th < nth; th++)
 533                 {
 534                     try
 535                     {
 536                         int calcvir_atom_end;
 537
 538                         if (vir == nullptr)
 539                         {
 540                             calcvir_atom_end = 0;
 541                         }
 542                         else
 543                         {
 544                             calcvir_atom_end = md.homenr;
 545                         }
 546
 547                         if (th > 0)
 548                         {
 549                             clear_mat(vir_r_m_dr_th[th]);
 550                         }
 551
 552                         int start_th = (nsettle * th) / nth;
 553                         int end_th   = (nsettle * (th + 1)) / nth;
 554
 555                         if (start_th >= 0 && end_th - start_th > 0)
 556                         {
 557                             settle_proj(settled, econq, end_th - start_th,
 558                                         settle.iatoms.data() + start_th * (1 + NRAL(F_SETTLE)), pbc_null,
 559                                         x.unpaddedArrayRef(), xprime.unpaddedArrayRef(), min_proj,
 560                                         calcvir_atom_end, th == 0 ? vir_r_m_dr : vir_r_m_dr_th[th]);
 561                         }
 562                     }
 563                     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
 564                 }
 565                 /* This is an overestimate */
 566                 inc_nrnb(nrnb, eNR_SETTLE, nsettle);
 567                 break;
 568             case ConstraintVariable::Deriv_FlexCon:
 569                 /* Nothing to do, since the are no flexible constraints in settles */
 570                 break;
 571             default: gmx_incons("Unknown constraint quantity for settle");
 572         }
 573
 574         if (vir != nullptr)
 575         {
 576             /* Reduce the virial contributions over the threads */
 577             for (int th = 1; th < nth; th++)
 578             {
 579                 m_add(vir_r_m_dr, vir_r_m_dr_th[th], vir_r_m_dr);
 580             }
 581         }
 582
 583         if (econq == ConstraintVariable::Positions)
 584         {
 585             for (int th = 1; th < nth; th++)
 586             {
 587                 bSettleErrorHasOccurred0 = bSettleErrorHasOccurred0 || bSettleErrorHasOccurred[th];
 588             }
 589
 590             if (bSettleErrorHasOccurred0)
 591             {
 592                 char buf[STRLEN];
 593                 sprintf(buf,
 594                         "\nstep "
 595                         "%" PRId64
 596                         ": One or more water molecules can not be settled.\n"
 597                         "Check for bad contacts and/or reduce the timestep if appropriate.\n",
 598                         step);
 599                 if (log)
 600                 {
 601                     fprintf(log, "%s", buf);
 602                 }
 603                 fprintf(stderr, "%s", buf);
 604                 warncount_settle++;
 605                 if (warncount_settle > maxwarn)
 606                 {
 607                     too_many_constraint_warnings(-1, warncount_settle);
 608                 }
 609                 bDump = TRUE;
 610
 611                 bOK = FALSE;
 612             }
 613         }
 614     }
 615
 616     if (vir != nullptr)
 617     {
 618         /* The normal uses of constrain() pass step_scaling = 1.0.
 619          * The call to constrain() for SD1 that passes step_scaling =
 620          * 0.5 also passes vir = NULL, so cannot reach this
 621          * assertion. This assertion should remain until someone knows
 622          * that this path works for their intended purpose, and then
 623          * they can use scaled_delta_t instead of ir.delta_t
 624          * below. */
 625         assert(gmx_within_tol(step_scaling, 1.0, GMX_REAL_EPS));
 626         switch (econq)
 627         {
 628             case ConstraintVariable::Positions: vir_fac = 0.5 / (ir.delta_t * ir.delta_t); break;
 629             case ConstraintVariable::Velocities: vir_fac = 0.5 / ir.delta_t; break;
 630             case ConstraintVariable::Force:
 631             case ConstraintVariable::ForceDispl: vir_fac = 0.5; break;
 632             default: gmx_incons("Unsupported constraint quantity for virial");
 633         }
 634
 635         if (EI_VV(ir.eI))
 636         {
 637             vir_fac *= 2; /* only constraining over half the distance here */
 638         }
 639         for (int i = 0; i < DIM; i++)
 640         {
 641             for (int j = 0; j < DIM; j++)
 642             {
 643                 (*vir)[i][j] = vir_fac * vir_r_m_dr[i][j];
 644             }
 645         }
 646     }
 647
 648     if (bDump)
 649     {
 650         dump_confs(log, step, mtop, start, homenr, cr, x.unpaddedArrayRef(),
 651                    xprime.unpaddedArrayRef(), box);
 652     }
 653
 654     if (econq == ConstraintVariable::Positions)
 655     {
 656         if (ir.bPull && pull_have_constraint(pull_work))
 657         {
 658             if (EI_DYNAMICS(ir.eI))
 659             {
 660                 t = ir.init_t + (step + delta_step) * ir.delta_t;
 661             }
 662             else
 663             {
 664                 t = ir.init_t;
 665             }
 666             set_pbc(&pbc, ir.pbcType, box);
 667             pull_constraint(pull_work, &md, &pbc, cr, ir.delta_t, t,
 668                             as_rvec_array(x.unpaddedArrayRef().data()),
 669                             as_rvec_array(xprime.unpaddedArrayRef().data()),
 670                             as_rvec_array(v.unpaddedArrayRef().data()), *vir);
 671         }
 672         if (ed && delta_step > 0)
 673         {
 674             /* apply the essential dynamics constraints here */
 675             do_edsam(&ir, step, cr, as_rvec_array(xprime.unpaddedArrayRef().data()),
 676                      as_rvec_array(v.unpaddedArrayRef().data()), box, ed);
 677         }
 678     }
 679     wallcycle_stop(wcycle, ewcCONSTR);
 680
 681     if (!v.empty() && md.cFREEZE)
 682     {
 683         /* Set the velocities of frozen dimensions to zero */
 684         ArrayRef<RVec> vRef = v.unpaddedArrayRef();
 685
 686         int gmx_unused numThreads = gmx_omp_nthreads_get(emntUpdate);
 687
 688 #pragma omp parallel for num_threads(numThreads) schedule(static)
 689         for (int i = 0; i < md.homenr; i++)
 690         {
 691             int freezeGroup = md.cFREEZE[i];
 692
 693             for (int d = 0; d < DIM; d++)
 694             {
 695                 if (ir.opts.nFreeze[freezeGroup][d])
 696                 {
 697                     vRef[i][d] = 0;
 698                 }
 699             }
 700         }
 701     }
 702
 703     return bOK;
 704 }
 705
 706 ArrayRef<real> Constraints::rmsdData() const
 707 {
 708     if (impl_->lincsd)
 709     {
 710         return lincs_rmsdData(impl_->lincsd);
 711     }
 712     else
 713     {
 714         return {};
 715     }
 716 }
 717
 718 real Constraints::rmsd() const
 719 {
 720     if (impl_->lincsd)
 721     {
 722         return lincs_rmsd(impl_->lincsd);
 723     }
 724     else
 725     {
 726         return 0;
 727     }
 728 }
 729
 730 int Constraints::numConstraintsTotal()
 731 {
 732     return impl_->ncon_tot;
 733 }
 734
 735 FlexibleConstraintTreatment flexibleConstraintTreatment(bool haveDynamicsIntegrator)
 736 {
 737     if (haveDynamicsIntegrator)
 738     {
 739         return FlexibleConstraintTreatment::Include;
 740     }
 741     else
 742     {
 743         return FlexibleConstraintTreatment::Exclude;
 744     }
 745 }
 746
 747 /*! \brief Returns a block struct to go from atoms to constraints
 748  *
 749  * The block struct will contain constraint indices with lower indices
 750  * directly matching the order in F_CONSTR and higher indices matching
 751  * the order in F_CONSTRNC offset by the number of constraints in F_CONSTR.
 752  *
 753  * \param[in]  numAtoms  The number of atoms to construct the list for
 754  * \param[in]  ilists    The interaction lists, size F_NRE
 755  * \param[in]  iparams   Interaction parameters, can be null when
 756  *                       \p flexibleConstraintTreatment==Include
 757  * \param[in]  flexibleConstraintTreatment  The flexible constraint treatment,
 758  *                                          see enum above
 759  *
 760  * \returns a block struct with all constraints for each atom
 761  */
 762 static ListOfLists<int> makeAtomsToConstraintsList(int                             numAtoms,
 763                                                    ArrayRef<const InteractionList> ilists,
 764                                                    ArrayRef<const t_iparams>       iparams,
 765                                                    FlexibleConstraintTreatment flexibleConstraintTreatment)
 766 {
 767     GMX_ASSERT(flexibleConstraintTreatment == FlexibleConstraintTreatment::Include || !iparams.empty(),
 768                "With flexible constraint detection we need valid iparams");
 769
 770     std::vector<int> count(numAtoms);
 771
 772     for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
 773     {
 774         const InteractionList& ilist  = ilists[ftype];
 775         const int              stride = 1 + NRAL(ftype);
 776         for (int i = 0; i < ilist.size(); i += stride)
 777         {
 778             if (flexibleConstraintTreatment == FlexibleConstraintTreatment::Include
 779                 || !isConstraintFlexible(iparams, ilist.iatoms[i]))
 780             {
 781                 for (int j = 1; j < 3; j++)
 782                 {
 783                     int a = ilist.iatoms[i + j];
 784                     count[a]++;
 785                 }
 786             }
 787         }
 788     }
 789
 790     std::vector<int> listRanges(numAtoms + 1);
 791     for (int a = 0; a < numAtoms; a++)
 792     {
 793         listRanges[a + 1] = listRanges[a] + count[a];
 794         count[a]          = 0;
 795     }
 796     std::vector<int> elements(listRanges[numAtoms]);
 797
 798     /* The F_CONSTRNC constraints have constraint numbers
 799      * that continue after the last F_CONSTR constraint.
 800      */
 801     int numConstraints = 0;
 802     for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
 803     {
 804         const InteractionList& ilist  = ilists[ftype];
 805         const int              stride = 1 + NRAL(ftype);
 806         for (int i = 0; i < ilist.size(); i += stride)
 807         {
 808             if (flexibleConstraintTreatment == FlexibleConstraintTreatment::Include
 809                 || !isConstraintFlexible(iparams, ilist.iatoms[i]))
 810             {
 811                 for (int j = 1; j < 3; j++)
 812                 {
 813                     const int a                          = ilist.iatoms[i + j];
 814                     elements[listRanges[a] + count[a]++] = numConstraints;
 815                 }
 816             }
 817             numConstraints++;
 818         }
 819     }
 820
 821     return ListOfLists<int>(std::move(listRanges), std::move(elements));
 822 }
 823
 824 ListOfLists<int> make_at2con(int                             numAtoms,
 825                              ArrayRef<const InteractionList> ilist,
 826                              ArrayRef<const t_iparams>       iparams,
 827                              FlexibleConstraintTreatment     flexibleConstraintTreatment)
 828 {
 829     return makeAtomsToConstraintsList(numAtoms, ilist, iparams, flexibleConstraintTreatment);
 830 }
 831
 832 ListOfLists<int> make_at2con(const gmx_moltype_t&           moltype,
 833                              gmx::ArrayRef<const t_iparams> iparams,
 834                              FlexibleConstraintTreatment    flexibleConstraintTreatment)
 835 {
 836     return makeAtomsToConstraintsList(moltype.atoms.nr, makeConstArrayRef(moltype.ilist), iparams,
 837                                       flexibleConstraintTreatment);
 838 }
 839
 840 //! Return the number of flexible constraints in the \c ilist and \c iparams.
 841 int countFlexibleConstraints(ArrayRef<const InteractionList> ilist, ArrayRef<const t_iparams> iparams)
 842 {
 843     int nflexcon = 0;
 844     for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
 845     {
 846         const int numIatomsPerConstraint = 3;
 847         for (int i = 0; i < ilist[ftype].size(); i += numIatomsPerConstraint)
 848         {
 849             const int type = ilist[ftype].iatoms[i];
 850             if (iparams[type].constr.dA == 0 && iparams[type].constr.dB == 0)
 851             {
 852                 nflexcon++;
 853             }
 854         }
 855     }
 856
 857     return nflexcon;
 858 }
 859
 860 //! Returns the index of the settle to which each atom belongs.
 861 static std::vector<int> make_at2settle(int natoms, const InteractionList& ilist)
 862 {
 863     /* Set all to no settle */
 864     std::vector<int> at2s(natoms, -1);
 865
 866     const int stride = 1 + NRAL(F_SETTLE);
 867
 868     for (int s = 0; s < ilist.size(); s += stride)
 869     {
 870         at2s[ilist.iatoms[s + 1]] = s / stride;
 871         at2s[ilist.iatoms[s + 2]] = s / stride;
 872         at2s[ilist.iatoms[s + 3]] = s / stride;
 873     }
 874
 875     return at2s;
 876 }
 877
 878 void Constraints::Impl::setConstraints(gmx_localtop_t* top, const t_mdatoms& md)
 879 {
 880     idef = &top->idef;
 881
 882     if (ncon_tot > 0)
 883     {
 884         /* With DD we might also need to call LINCS on a domain no constraints for
 885          * communicating coordinates to other nodes that do have constraints.
 886          */
 887         if (ir.eConstrAlg == econtLINCS)
 888         {
 889             set_lincs(*idef, md, EI_DYNAMICS(ir.eI), cr, lincsd);
 890         }
 891         if (ir.eConstrAlg == econtSHAKE)
 892         {
 893             if (cr->dd)
 894             {
 895                 // We are using the local topology, so there are only
 896                 // F_CONSTR constraints.
 897                 make_shake_sblock_dd(shaked, idef->il[F_CONSTR], cr->dd);
 898             }
 899             else
 900             {
 901                 make_shake_sblock_serial(shaked, &top->idef, md);
 902             }
 903         }
 904     }
 905
 906     if (settled)
 907     {
 908         settle_set_constraints(settled, idef->il[F_SETTLE], md);
 909     }
 910
 911     /* Make a selection of the local atoms for essential dynamics */
 912     if (ed && cr->dd)
 913     {
 914         dd_make_local_ed_indices(cr->dd, ed);
 915     }
 916 }
 917
 918 void Constraints::setConstraints(gmx_localtop_t* top, const t_mdatoms& md)
 919 {
 920     impl_->setConstraints(top, md);
 921 }
 922
 923 /*! \brief Makes a per-moleculetype container of mappings from atom
 924  * indices to constraint indices.
 925  *
 926  * Note that flexible constraints are only enabled with a dynamical integrator. */
 927 static std::vector<ListOfLists<int>> makeAtomToConstraintMappings(const gmx_mtop_t& mtop,
 928                                                                   FlexibleConstraintTreatment flexibleConstraintTreatment)
 929 {
 930     std::vector<ListOfLists<int>> mapping;
 931     mapping.reserve(mtop.moltype.size());
 932     for (const gmx_moltype_t& moltype : mtop.moltype)
 933     {
 934         mapping.push_back(make_at2con(moltype, mtop.ffparams.iparams, flexibleConstraintTreatment));
 935     }
 936     return mapping;
 937 }
 938
 939 Constraints::Constraints(const gmx_mtop_t&     mtop,
 940                          const t_inputrec&     ir,
 941                          pull_t*               pull_work,
 942                          FILE*                 log,
 943                          const t_mdatoms&      md,
 944                          const t_commrec*      cr,
 945                          const gmx_multisim_t* ms,
 946                          t_nrnb*               nrnb,
 947                          gmx_wallcycle*        wcycle,
 948                          bool                  pbcHandlingRequired,
 949                          int                   numConstraints,
 950                          int                   numSettles) :
 951     impl_(new Impl(mtop, ir, pull_work, log, md, cr, ms, nrnb, wcycle, pbcHandlingRequired, numConstraints, numSettles))
 952 {
 953 }
 954
 955 Constraints::Impl::Impl(const gmx_mtop_t&     mtop_p,
 956                         const t_inputrec&     ir_p,
 957                         pull_t*               pull_work,
 958                         FILE*                 log_p,
 959                         const t_mdatoms&      md_p,
 960                         const t_commrec*      cr_p,
 961                         const gmx_multisim_t* ms_p,
 962                         t_nrnb*               nrnb_p,
 963                         gmx_wallcycle*        wcycle_p,
 964                         bool                  pbcHandlingRequired,
 965                         int                   numConstraints,
 966                         int                   numSettles) :
 967     ncon_tot(numConstraints),
 968     mtop(mtop_p),
 969     md(md_p),
 970     pbcHandlingRequired_(pbcHandlingRequired),
 971     log(log_p),
 972     cr(cr_p),
 973     ms(ms_p),
 974     pull_work(pull_work),
 975     ir(ir_p),
 976     nrnb(nrnb_p),
 977     wcycle(wcycle_p)
 978 {
 979     if (numConstraints + numSettles > 0 && ir.epc == epcMTTK)
 980     {
 981         gmx_fatal(FARGS, "Constraints are not implemented with MTTK pressure control.");
 982     }
 983
 984     nflexcon = 0;
 985     if (numConstraints > 0)
 986     {
 987         at2con_mt = makeAtomToConstraintMappings(mtop, flexibleConstraintTreatment(EI_DYNAMICS(ir.eI)));
 988
 989         for (const gmx_molblock_t& molblock : mtop.molblock)
 990         {
 991             int count = countFlexibleConstraints(mtop.moltype[molblock.type].ilist, mtop.ffparams.iparams);
 992             nflexcon += molblock.nmol * count;
 993         }
 994
 995         if (nflexcon > 0)
 996         {
 997             if (log)
 998             {
 999                 fprintf(log, "There are %d flexible constraints\n", nflexcon);
1000                 if (ir.fc_stepsize == 0)
1001                 {
1002                     fprintf(log,
1003                             "\n"
1004                             "WARNING: step size for flexible constraining = 0\n"
1005                             "         All flexible constraints will be rigid.\n"
1006                             "         Will try to keep all flexible constraints at their original "
1007                             "length,\n"
1008                             "         but the lengths may exhibit some drift.\n\n");
1009                     nflexcon = 0;
1010                 }
1011             }
1012             if (nflexcon > 0)
1013             {
1014                 please_cite(log, "Hess2002");
1015             }
1016         }
1017
1018         if (ir.eConstrAlg == econtLINCS)
1019         {
1020             lincsd = init_lincs(log, mtop, nflexcon, at2con_mt,
1021                                 DOMAINDECOMP(cr) && ddHaveSplitConstraints(*cr->dd), ir.nLincsIter,
1022                                 ir.nProjOrder);
1023         }
1024
1025         if (ir.eConstrAlg == econtSHAKE)
1026         {
1027             if (DOMAINDECOMP(cr) && ddHaveSplitConstraints(*cr->dd))
1028             {
1029                 gmx_fatal(FARGS,
1030                           "SHAKE is not supported with domain decomposition and constraint that "
1031                           "cross domain boundaries, use LINCS");
1032             }
1033             if (nflexcon)
1034             {
1035                 gmx_fatal(FARGS,
1036                           "For this system also velocities and/or forces need to be constrained, "
1037                           "this can not be done with SHAKE, you should select LINCS");
1038             }
1039             please_cite(log, "Ryckaert77a");
1040             if (ir.bShakeSOR)
1041             {
1042                 please_cite(log, "Barth95a");
1043             }
1044
1045             shaked = shake_init();
1046         }
1047     }
1048
1049     if (numSettles > 0)
1050     {
1051         please_cite(log, "Miyamoto92a");
1052
1053         settled = settle_init(mtop);
1054
1055         /* Make an atom to settle index for use in domain decomposition */
1056         for (size_t mt = 0; mt < mtop.moltype.size(); mt++)
1057         {
1058             at2settle_mt.emplace_back(
1059                     make_at2settle(mtop.moltype[mt].atoms.nr, mtop.moltype[mt].ilist[F_SETTLE]));
1060         }
1061
1062         /* Allocate thread-local work arrays */
1063         int nthreads = gmx_omp_nthreads_get(emntSETTLE);
1064         if (nthreads > 1 && vir_r_m_dr_th == nullptr)
1065         {
1066             snew(vir_r_m_dr_th, nthreads);
1067             snew(bSettleErrorHasOccurred, nthreads);
1068         }
1069     }
1070
1071     maxwarn   = 999;
1072     char* env = getenv("GMX_MAXCONSTRWARN");
1073     if (env)
1074     {
1075         maxwarn = 0;
1076         sscanf(env, "%8d", &maxwarn);
1077         if (maxwarn < 0)
1078         {
1079             maxwarn = INT_MAX;
1080         }
1081         if (log)
1082         {
1083             fprintf(log, "Setting the maximum number of constraint warnings to %d\n", maxwarn);
1084         }
1085         if (MASTER(cr))
1086         {
1087             fprintf(stderr, "Setting the maximum number of constraint warnings to %d\n", maxwarn);
1088         }
1089     }
1090     warncount_lincs  = 0;
1091     warncount_settle = 0;
1092 }
1093
1094 Constraints::Impl::~Impl()
1095 {
1096     if (bSettleErrorHasOccurred != nullptr)
1097     {
1098         sfree(bSettleErrorHasOccurred);
1099     }
1100     if (vir_r_m_dr_th != nullptr)
1101     {
1102         sfree(vir_r_m_dr_th);
1103     }
1104     if (settled != nullptr)
1105     {
1106         settle_free(settled);
1107     }
1108     done_lincs(lincsd);
1109 }
1110
1111 void Constraints::saveEdsamPointer(gmx_edsam* ed)
1112 {
1113     impl_->ed = ed;
1114 }
1115
1116 ArrayRef<const ListOfLists<int>> Constraints::atom2constraints_moltype() const
1117 {
1118     return impl_->at2con_mt;
1119 }
1120
1121 ArrayRef<const std::vector<int>> Constraints::atom2settle_moltype() const
1122 {
1123     return impl_->at2settle_mt;
1124 }
1125
1126 void do_constrain_first(FILE*                     fplog,
1127                         gmx::Constraints*         constr,
1128                         const t_inputrec*         ir,
1129                         const t_mdatoms*          md,
1130                         int                       natoms,
1131                         ArrayRefWithPadding<RVec> x,
1132                         ArrayRefWithPadding<RVec> v,
1133                         const matrix              box,
1134                         real                      lambda)
1135 {
1136     int     i, m, start, end;
1137     int64_t step;
1138     real    dt = ir->delta_t;
1139     real    dvdl_dum;
1140
1141     PaddedVector<RVec> savex(natoms);
1142
1143     start = 0;
1144     end   = md->homenr;
1145
1146     if (debug)
1147     {
1148         fprintf(debug, "vcm: start=%d, homenr=%d, end=%d\n", start, md->homenr, end);
1149     }
1150     /* Do a first constrain to reset particles... */
1151     step = ir->init_step;
1152     if (fplog)
1153     {
1154         char buf[STEPSTRSIZE];
1155         fprintf(fplog, "\nConstraining the starting coordinates (step %s)\n", gmx_step_str(step, buf));
1156     }
1157     dvdl_dum = 0;
1158
1159     /* constrain the current position */
1160     constr->apply(TRUE, FALSE, step, 0, 1.0, x, x, {}, box, lambda, &dvdl_dum, {}, nullptr,
1161                   gmx::ConstraintVariable::Positions);
1162     if (EI_VV(ir->eI))
1163     {
1164         /* constrain the inital velocity, and save it */
1165         /* also may be useful if we need the ekin from the halfstep for velocity verlet */
1166         constr->apply(TRUE, FALSE, step, 0, 1.0, x, v, v.unpaddedArrayRef(), box, lambda, &dvdl_dum,
1167                       {}, nullptr, gmx::ConstraintVariable::Velocities);
1168     }
1169     /* constrain the inital velocities at t-dt/2 */
1170     if (EI_STATE_VELOCITY(ir->eI) && ir->eI != eiVV)
1171     {
1172         auto subX = x.paddedArrayRef().subArray(start, end);
1173         auto subV = v.paddedArrayRef().subArray(start, end);
1174         for (i = start; (i < end); i++)
1175         {
1176             for (m = 0; (m < DIM); m++)
1177             {
1178                 /* Reverse the velocity */
1179                 subV[i][m] = -subV[i][m];
1180                 /* Store the position at t-dt in buf */
1181                 savex[i][m] = subX[i][m] + dt * subV[i][m];
1182             }
1183         }
1184         /* Shake the positions at t=-dt with the positions at t=0
1185          * as reference coordinates.
1186          */
1187         if (fplog)
1188         {
1189             char buf[STEPSTRSIZE];
1190             fprintf(fplog, "\nConstraining the coordinates at t0-dt (step %s)\n", gmx_step_str(step, buf));
1191         }
1192         dvdl_dum = 0;
1193         constr->apply(TRUE, FALSE, step, -1, 1.0, x, savex.arrayRefWithPadding(), {}, box, lambda,
1194                       &dvdl_dum, v, nullptr, gmx::ConstraintVariable::Positions);
1195
1196         for (i = start; i < end; i++)
1197         {
1198             for (m = 0; m < DIM; m++)
1199             {
1200                 /* Re-reverse the velocities */
1201                 subV[i][m] = -subV[i][m];
1202             }
1203         }
1204     }
1205 }
1206
1207 } // namespace gmx