src/gromacs/mdlib/constr.cpp

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