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 *           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         int gmx_unused numThreads = gmx_omp_nthreads_get(emntUpdate);
 690
 691 #pragma omp parallel for num_threads(numThreads) schedule(static)
 692         for (int i = 0; i < md.homenr; i++)
 693         {
 694             int freezeGroup = md.cFREEZE[i];
 695
 696             for (int d = 0; d < DIM; d++)
 697             {
 698                 if (ir.opts.nFreeze[freezeGroup][d])
 699                 {
 700                     v[i][d] = 0;
 701                 }
 702             }
 703         }
 704     }
 705
 706     return bOK;
 707 }
 708
 709 ArrayRef<real> Constraints::rmsdData() const
 710 {
 711     if (impl_->lincsd)
 712     {
 713         return lincs_rmsdData(impl_->lincsd);
 714     }
 715     else
 716     {
 717         return EmptyArrayRef();
 718     }
 719 }
 720
 721 real Constraints::rmsd() const
 722 {
 723     if (impl_->lincsd)
 724     {
 725         return lincs_rmsd(impl_->lincsd);
 726     }
 727     else
 728     {
 729         return 0;
 730     }
 731 }
 732
 733 FlexibleConstraintTreatment flexibleConstraintTreatment(bool haveDynamicsIntegrator)
 734 {
 735     if (haveDynamicsIntegrator)
 736     {
 737         return FlexibleConstraintTreatment::Include;
 738     }
 739     else
 740     {
 741         return FlexibleConstraintTreatment::Exclude;
 742     }
 743 }
 744
 745 t_blocka make_at2con(int                          numAtoms,
 746                      const t_ilist               *ilists,
 747                      const t_iparams             *iparams,
 748                      FlexibleConstraintTreatment  flexibleConstraintTreatment)
 749 {
 750     GMX_ASSERT(flexibleConstraintTreatment == FlexibleConstraintTreatment::Include || iparams != nullptr, "With flexible constraint detection we need valid iparams");
 751
 752     std::vector<int> count(numAtoms);
 753
 754     for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
 755     {
 756         const t_ilist &ilist  = ilists[ftype];
 757         const int      stride = 1 + NRAL(ftype);
 758         for (int i = 0; i < ilist.nr; i += stride)
 759         {
 760             if (flexibleConstraintTreatment == FlexibleConstraintTreatment::Include ||
 761                 !isConstraintFlexible(iparams, ilist.iatoms[i]))
 762             {
 763                 for (int j = 1; j < 3; j++)
 764                 {
 765                     int a = ilist.iatoms[i + j];
 766                     count[a]++;
 767                 }
 768             }
 769         }
 770     }
 771
 772     t_blocka at2con;
 773     at2con.nr           = numAtoms;
 774     at2con.nalloc_index = at2con.nr + 1;
 775     snew(at2con.index, at2con.nalloc_index);
 776     at2con.index[0] = 0;
 777     for (int a = 0; a < numAtoms; a++)
 778     {
 779         at2con.index[a + 1] = at2con.index[a] + count[a];
 780         count[a]            = 0;
 781     }
 782     at2con.nra      = at2con.index[at2con.nr];
 783     at2con.nalloc_a = at2con.nra;
 784     snew(at2con.a, at2con.nalloc_a);
 785
 786     /* The F_CONSTRNC constraints have constraint numbers
 787      * that continue after the last F_CONSTR constraint.
 788      */
 789     int numConstraints = 0;
 790     for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
 791     {
 792         const t_ilist &ilist  = ilists[ftype];
 793         const int      stride = 1 + NRAL(ftype);
 794         for (int i = 0; i < ilist.nr; i += stride)
 795         {
 796             if (flexibleConstraintTreatment == FlexibleConstraintTreatment::Include ||
 797                 !isConstraintFlexible(iparams, ilist.iatoms[i]))
 798             {
 799                 for (int j = 1; j < 3; j++)
 800                 {
 801                     int a = ilist.iatoms[i + j];
 802                     at2con.a[at2con.index[a] + count[a]++] = numConstraints;
 803                 }
 804             }
 805             numConstraints++;
 806         }
 807     }
 808
 809     return at2con;
 810 }
 811
 812 int countFlexibleConstraints(const t_ilist   *ilist,
 813                              const t_iparams *iparams)
 814 {
 815     int nflexcon = 0;
 816     for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
 817     {
 818         const int numIatomsPerConstraint = 3;
 819         int       ncon                   = ilist[ftype].nr /  numIatomsPerConstraint;
 820         t_iatom  *ia                     = ilist[ftype].iatoms;
 821         for (int con = 0; con < ncon; con++)
 822         {
 823             if (iparams[ia[0]].constr.dA == 0 &&
 824                 iparams[ia[0]].constr.dB == 0)
 825             {
 826                 nflexcon++;
 827             }
 828             ia += numIatomsPerConstraint;
 829         }
 830     }
 831
 832     return nflexcon;
 833 }
 834
 835 //! Returns the index of the settle to which each atom belongs.
 836 static int *make_at2settle(int natoms, const t_ilist *ilist)
 837 {
 838     int *at2s;
 839     int  a, stride, s;
 840
 841     snew(at2s, natoms);
 842     /* Set all to no settle */
 843     for (a = 0; a < natoms; a++)
 844     {
 845         at2s[a] = -1;
 846     }
 847
 848     stride = 1 + NRAL(F_SETTLE);
 849
 850     for (s = 0; s < ilist->nr; s += stride)
 851     {
 852         at2s[ilist->iatoms[s+1]] = s/stride;
 853         at2s[ilist->iatoms[s+2]] = s/stride;
 854         at2s[ilist->iatoms[s+3]] = s/stride;
 855     }
 856
 857     return at2s;
 858 }
 859
 860 void
 861 Constraints::Impl::setConstraints(const gmx_localtop_t &top,
 862                                   const t_mdatoms      &md)
 863 {
 864     idef = &top.idef;
 865
 866     if (ncon_tot > 0)
 867     {
 868         /* With DD we might also need to call LINCS on a domain no constraints for
 869          * communicating coordinates to other nodes that do have constraints.
 870          */
 871         if (ir.eConstrAlg == econtLINCS)
 872         {
 873             set_lincs(top.idef, md, EI_DYNAMICS(ir.eI), cr, lincsd);
 874         }
 875         if (ir.eConstrAlg == econtSHAKE)
 876         {
 877             if (cr->dd)
 878             {
 879                 // We are using the local topology, so there are only
 880                 // F_CONSTR constraints.
 881                 make_shake_sblock_dd(shaked, &idef->il[F_CONSTR], &top.cgs, cr->dd);
 882             }
 883             else
 884             {
 885                 make_shake_sblock_serial(shaked, idef, md);
 886             }
 887         }
 888     }
 889
 890     if (settled)
 891     {
 892         settle_set_constraints(settled,
 893                                &idef->il[F_SETTLE], md);
 894     }
 895
 896     /* Make a selection of the local atoms for essential dynamics */
 897     if (ed && cr->dd)
 898     {
 899         dd_make_local_ed_indices(cr->dd, ed);
 900     }
 901 }
 902
 903 void
 904 Constraints::setConstraints(const gmx_localtop_t &top,
 905                             const t_mdatoms      &md)
 906 {
 907     impl_->setConstraints(top, md);
 908 }
 909
 910 /*! \brief Makes a per-moleculetype container of mappings from atom
 911  * indices to constraint indices.
 912  *
 913  * Note that flexible constraints are only enabled with a dynamical integrator. */
 914 static std::vector<t_blocka>
 915 makeAtomToConstraintMappings(const gmx_mtop_t            &mtop,
 916                              FlexibleConstraintTreatment  flexibleConstraintTreatment)
 917 {
 918     std::vector<t_blocka> mapping;
 919     mapping.reserve(mtop.moltype.size());
 920     for (const gmx_moltype_t &moltype : mtop.moltype)
 921     {
 922         mapping.push_back(make_at2con(moltype.atoms.nr,
 923                                       moltype.ilist,
 924                                       mtop.ffparams.iparams,
 925                                       flexibleConstraintTreatment));
 926     }
 927     return mapping;
 928 }
 929
 930 Constraints::Constraints(const gmx_mtop_t     &mtop,
 931                          const t_inputrec     &ir,
 932                          FILE                 *log,
 933                          const t_mdatoms      &md,
 934                          const t_commrec      *cr,
 935                          const gmx_multisim_t &ms,
 936                          t_nrnb               *nrnb,
 937                          gmx_wallcycle        *wcycle,
 938                          bool                  pbcHandlingRequired,
 939                          int                   numConstraints,
 940                          int                   numSettles)
 941     : impl_(new Impl(mtop,
 942                      ir,
 943                      log,
 944                      md,
 945                      cr,
 946                      ms,
 947                      nrnb,
 948                      wcycle,
 949                      pbcHandlingRequired,
 950                      numConstraints,
 951                      numSettles))
 952 {
 953 }
 954
 955 Constraints::Impl::Impl(const gmx_mtop_t     &mtop_p,
 956                         const t_inputrec     &ir_p,
 957                         FILE                 *log_p,
 958                         const t_mdatoms      &md_p,
 959                         const t_commrec      *cr_p,
 960                         const gmx_multisim_t &ms_p,
 961                         t_nrnb               *nrnb_p,
 962                         gmx_wallcycle        *wcycle_p,
 963                         bool                  pbcHandlingRequired,
 964                         int                   numConstraints,
 965                         int                   numSettles)
 966     : ncon_tot(numConstraints),
 967       mtop(mtop_p),
 968       md(md_p),
 969       pbcHandlingRequired_(pbcHandlingRequired),
 970       log(log_p),
 971       cr(cr_p),
 972       ms(ms_p),
 973       ir(ir_p),
 974       nrnb(nrnb_p),
 975       wcycle(wcycle_p)
 976 {
 977     if (numConstraints + numSettles > 0 && ir.epc == epcMTTK)
 978     {
 979         gmx_fatal(FARGS, "Constraints are not implemented with MTTK pressure control.");
 980     }
 981
 982     nflexcon = 0;
 983     if (numConstraints > 0)
 984     {
 985         at2con_mt = makeAtomToConstraintMappings(mtop,
 986                                                  flexibleConstraintTreatment(EI_DYNAMICS(ir.eI)));
 987
 988         for (const gmx_molblock_t &molblock : mtop.molblock)
 989         {
 990             int count = countFlexibleConstraints(mtop.moltype[molblock.type].ilist,
 991                                                  mtop.ffparams.iparams);
 992             nflexcon += molblock.nmol*count;
 993         }
 994
 995         if (nflexcon > 0)
 996         {
 997             if (log)
 998             {
 999                 fprintf(log, "There are %d flexible constraints\n",
1000                         nflexcon);
1001                 if (ir.fc_stepsize == 0)
1002                 {
1003                     fprintf(log, "\n"
1004                             "WARNING: step size for flexible constraining = 0\n"
1005                             "         All flexible constraints will be rigid.\n"
1006                             "         Will try to keep all flexible constraints at their original length,\n"
1007                             "         but the lengths may exhibit some drift.\n\n");
1008                     nflexcon = 0;
1009                 }
1010             }
1011             if (nflexcon > 0)
1012             {
1013                 please_cite(log, "Hess2002");
1014             }
1015         }
1016
1017         if (ir.eConstrAlg == econtLINCS)
1018         {
1019             lincsd = init_lincs(log, mtop,
1020                                 nflexcon, at2con_mt,
1021                                 DOMAINDECOMP(cr) && cr->dd->bInterCGcons,
1022                                 ir.nLincsIter, ir.nProjOrder);
1023         }
1024
1025         if (ir.eConstrAlg == econtSHAKE)
1026         {
1027             if (DOMAINDECOMP(cr) && cr->dd->bInterCGcons)
1028             {
1029                 gmx_fatal(FARGS, "SHAKE is not supported with domain decomposition and constraint that cross charge group boundaries, use LINCS");
1030             }
1031             if (nflexcon)
1032             {
1033                 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");
1034             }
1035             please_cite(log, "Ryckaert77a");
1036             if (ir.bShakeSOR)
1037             {
1038                 please_cite(log, "Barth95a");
1039             }
1040
1041             shaked = shake_init();
1042         }
1043     }
1044
1045     if (numSettles > 0)
1046     {
1047         please_cite(log, "Miyamoto92a");
1048
1049         bInterCGsettles = inter_charge_group_settles(mtop);
1050
1051         settled         = settle_init(mtop);
1052
1053         /* Make an atom to settle index for use in domain decomposition */
1054         n_at2settle_mt = mtop.moltype.size();
1055         snew(at2settle_mt, n_at2settle_mt);
1056         for (size_t mt = 0; mt < mtop.moltype.size(); mt++)
1057         {
1058             at2settle_mt[mt] =
1059                 make_at2settle(mtop.moltype[mt].atoms.nr,
1060                                &mtop.moltype[mt].ilist[F_SETTLE]);
1061         }
1062
1063         /* Allocate thread-local work arrays */
1064         int nthreads = gmx_omp_nthreads_get(emntSETTLE);
1065         if (nthreads > 1 && vir_r_m_dr_th == nullptr)
1066         {
1067             snew(vir_r_m_dr_th, nthreads);
1068             snew(bSettleErrorHasOccurred, nthreads);
1069         }
1070     }
1071
1072     maxwarn = 999;
1073     char *env       = getenv("GMX_MAXCONSTRWARN");
1074     if (env)
1075     {
1076         maxwarn = 0;
1077         sscanf(env, "%8d", &maxwarn);
1078         if (maxwarn < 0)
1079         {
1080             maxwarn = INT_MAX;
1081         }
1082         if (log)
1083         {
1084             fprintf(log,
1085                     "Setting the maximum number of constraint warnings to %d\n",
1086                     maxwarn);
1087         }
1088         if (MASTER(cr))
1089         {
1090             fprintf(stderr,
1091                     "Setting the maximum number of constraint warnings to %d\n",
1092                     maxwarn);
1093         }
1094     }
1095     warncount_lincs  = 0;
1096     warncount_settle = 0;
1097 }
1098
1099 Constraints::Impl::~Impl()
1100 {
1101     done_lincs(lincsd);
1102 }
1103
1104 void Constraints::saveEdsamPointer(gmx_edsam * ed)
1105 {
1106     impl_->ed = ed;
1107 }
1108
1109 const ArrayRef<const t_blocka>
1110 Constraints::atom2constraints_moltype() const
1111 {
1112     return impl_->at2con_mt;
1113 }
1114
1115 int *const* Constraints::atom2settle_moltype() const
1116 {
1117     return impl_->at2settle_mt;
1118 }
1119
1120
1121 bool inter_charge_group_constraints(const gmx_mtop_t &mtop)
1122 {
1123     const gmx_moltype_t *molt;
1124     const t_block       *cgs;
1125     const t_ilist       *il;
1126     int                 *at2cg, cg, a, ftype, i;
1127     bool                 bInterCG;
1128
1129     bInterCG = FALSE;
1130     for (size_t mb = 0; mb < mtop.molblock.size() && !bInterCG; mb++)
1131     {
1132         molt = &mtop.moltype[mtop.molblock[mb].type];
1133
1134         if (molt->ilist[F_CONSTR].nr   > 0 ||
1135             molt->ilist[F_CONSTRNC].nr > 0 ||
1136             molt->ilist[F_SETTLE].nr > 0)
1137         {
1138             cgs  = &molt->cgs;
1139             snew(at2cg, molt->atoms.nr);
1140             for (cg = 0; cg < cgs->nr; cg++)
1141             {
1142                 for (a = cgs->index[cg]; a < cgs->index[cg+1]; a++)
1143                 {
1144                     at2cg[a] = cg;
1145                 }
1146             }
1147
1148             for (ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
1149             {
1150                 il = &molt->ilist[ftype];
1151                 for (i = 0; i < il->nr && !bInterCG; i += 1+NRAL(ftype))
1152                 {
1153                     if (at2cg[il->iatoms[i+1]] != at2cg[il->iatoms[i+2]])
1154                     {
1155                         bInterCG = TRUE;
1156                     }
1157                 }
1158             }
1159
1160             sfree(at2cg);
1161         }
1162     }
1163
1164     return bInterCG;
1165 }
1166
1167 bool inter_charge_group_settles(const gmx_mtop_t &mtop)
1168 {
1169     const gmx_moltype_t *molt;
1170     const t_block       *cgs;
1171     const t_ilist       *il;
1172     int                 *at2cg, cg, a, ftype, i;
1173     bool                 bInterCG;
1174
1175     bInterCG = FALSE;
1176     for (size_t mb = 0; mb < mtop.molblock.size() && !bInterCG; mb++)
1177     {
1178         molt = &mtop.moltype[mtop.molblock[mb].type];
1179
1180         if (molt->ilist[F_SETTLE].nr > 0)
1181         {
1182             cgs  = &molt->cgs;
1183             snew(at2cg, molt->atoms.nr);
1184             for (cg = 0; cg < cgs->nr; cg++)
1185             {
1186                 for (a = cgs->index[cg]; a < cgs->index[cg+1]; a++)
1187                 {
1188                     at2cg[a] = cg;
1189                 }
1190             }
1191
1192             for (ftype = F_SETTLE; ftype <= F_SETTLE; ftype++)
1193             {
1194                 il = &molt->ilist[ftype];
1195                 for (i = 0; i < il->nr && !bInterCG; i += 1+NRAL(F_SETTLE))
1196                 {
1197                     if (at2cg[il->iatoms[i+1]] != at2cg[il->iatoms[i+2]] ||
1198                         at2cg[il->iatoms[i+1]] != at2cg[il->iatoms[i+3]])
1199                     {
1200                         bInterCG = TRUE;
1201                     }
1202                 }
1203             }
1204
1205             sfree(at2cg);
1206         }
1207     }
1208
1209     return bInterCG;
1210 }
1211
1212 }  // namespace gmx