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