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