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