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