src/gromacs/mdlib/constr.cpp

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