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