src/gromacs/gmxpreprocess/grompp.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 #include "gmxpre.h"
  38
  39 #include "grompp.h"
  40
  41 #include <cerrno>
  42 #include <climits>
  43 #include <cmath>
  44 #include <cstring>
  45
  46 #include <algorithm>
  47 #include <memory>
  48 #include <vector>
  49
  50 #include <sys/types.h>
  51
  52 #include "gromacs/awh/read_params.h"
  53 #include "gromacs/commandline/pargs.h"
  54 #include "gromacs/ewald/ewald_utils.h"
  55 #include "gromacs/ewald/pme.h"
  56 #include "gromacs/fft/calcgrid.h"
  57 #include "gromacs/fileio/confio.h"
  58 #include "gromacs/fileio/enxio.h"
  59 #include "gromacs/fileio/tpxio.h"
  60 #include "gromacs/fileio/trxio.h"
  61 #include "gromacs/fileio/warninp.h"
  62 #include "gromacs/gmxpreprocess/add_par.h"
  63 #include "gromacs/gmxpreprocess/convparm.h"
  64 #include "gromacs/gmxpreprocess/gen_maxwell_velocities.h"
  65 #include "gromacs/gmxpreprocess/gpp_atomtype.h"
  66 #include "gromacs/gmxpreprocess/grompp_impl.h"
  67 #include "gromacs/gmxpreprocess/notset.h"
  68 #include "gromacs/gmxpreprocess/readir.h"
  69 #include "gromacs/gmxpreprocess/tomorse.h"
  70 #include "gromacs/gmxpreprocess/topio.h"
  71 #include "gromacs/gmxpreprocess/toputil.h"
  72 #include "gromacs/gmxpreprocess/vsite_parm.h"
  73 #include "gromacs/imd/imd.h"
  74 #include "gromacs/math/functions.h"
  75 #include "gromacs/math/invertmatrix.h"
  76 #include "gromacs/math/units.h"
  77 #include "gromacs/math/vec.h"
  78 #include "gromacs/mdlib/calc_verletbuf.h"
  79 #include "gromacs/mdlib/compute_io.h"
  80 #include "gromacs/mdlib/constr.h"
  81 #include "gromacs/mdlib/perf_est.h"
  82 #include "gromacs/mdlib/qmmm.h"
  83 #include "gromacs/mdlib/vsite.h"
  84 #include "gromacs/mdrun/mdmodules.h"
  85 #include "gromacs/mdtypes/inputrec.h"
  86 #include "gromacs/mdtypes/md_enums.h"
  87 #include "gromacs/mdtypes/nblist.h"
  88 #include "gromacs/mdtypes/state.h"
  89 #include "gromacs/pbcutil/boxutilities.h"
  90 #include "gromacs/pbcutil/pbc.h"
  91 #include "gromacs/pulling/pull.h"
  92 #include "gromacs/random/seed.h"
  93 #include "gromacs/topology/ifunc.h"
  94 #include "gromacs/topology/mtop_util.h"
  95 #include "gromacs/topology/symtab.h"
  96 #include "gromacs/topology/topology.h"
  97 #include "gromacs/trajectory/trajectoryframe.h"
  98 #include "gromacs/utility/arraysize.h"
  99 #include "gromacs/utility/cstringutil.h"
 100 #include "gromacs/utility/exceptions.h"
 101 #include "gromacs/utility/fatalerror.h"
 102 #include "gromacs/utility/futil.h"
 103 #include "gromacs/utility/gmxassert.h"
 104 #include "gromacs/utility/keyvaluetreebuilder.h"
 105 #include "gromacs/utility/mdmodulenotification.h"
 106 #include "gromacs/utility/smalloc.h"
 107 #include "gromacs/utility/snprintf.h"
 108
 109 /* TODO The implementation details should move to their own source file. */
 110 InteractionOfType::InteractionOfType(gmx::ArrayRef<const int>  atoms,
 111                                      gmx::ArrayRef<const real> params,
 112                                      const std::string&        name) :
 113     atoms_(atoms.begin(), atoms.end()),
 114     interactionTypeName_(name)
 115 {
 116     GMX_RELEASE_ASSERT(
 117             params.size() <= forceParam_.size(),
 118             gmx::formatString("Cannot have more parameters than the maximum number possible (%d)", MAXFORCEPARAM)
 119                     .c_str());
 120     auto forceParamIt = forceParam_.begin();
 121     for (const auto param : params)
 122     {
 123         *forceParamIt++ = param;
 124     }
 125     while (forceParamIt != forceParam_.end())
 126     {
 127         *forceParamIt++ = NOTSET;
 128     }
 129 }
 130
 131 const int& InteractionOfType::ai() const
 132 {
 133     GMX_RELEASE_ASSERT(!atoms_.empty(), "Need to have at least one atom set");
 134     return atoms_[0];
 135 }
 136
 137 const int& InteractionOfType::aj() const
 138 {
 139     GMX_RELEASE_ASSERT(atoms_.size() > 1, "Need to have at least two atoms set");
 140     return atoms_[1];
 141 }
 142
 143 const int& InteractionOfType::ak() const
 144 {
 145     GMX_RELEASE_ASSERT(atoms_.size() > 2, "Need to have at least three atoms set");
 146     return atoms_[2];
 147 }
 148
 149 const int& InteractionOfType::al() const
 150 {
 151     GMX_RELEASE_ASSERT(atoms_.size() > 3, "Need to have at least four atoms set");
 152     return atoms_[3];
 153 }
 154
 155 const int& InteractionOfType::am() const
 156 {
 157     GMX_RELEASE_ASSERT(atoms_.size() > 4, "Need to have at least five atoms set");
 158     return atoms_[4];
 159 }
 160
 161 const real& InteractionOfType::c0() const
 162 {
 163     return forceParam_[0];
 164 }
 165
 166 const real& InteractionOfType::c1() const
 167 {
 168     return forceParam_[1];
 169 }
 170
 171 const real& InteractionOfType::c2() const
 172 {
 173     return forceParam_[2];
 174 }
 175
 176 const std::string& InteractionOfType::interactionTypeName() const
 177 {
 178     return interactionTypeName_;
 179 }
 180
 181 void InteractionOfType::sortBondAtomIds()
 182 {
 183     if (aj() < ai())
 184     {
 185         std::swap(atoms_[0], atoms_[1]);
 186     }
 187 }
 188
 189 void InteractionOfType::sortAngleAtomIds()
 190 {
 191     if (ak() < ai())
 192     {
 193         std::swap(atoms_[0], atoms_[2]);
 194     }
 195 }
 196
 197 void InteractionOfType::sortDihedralAtomIds()
 198 {
 199     if (al() < ai())
 200     {
 201         std::swap(atoms_[0], atoms_[3]);
 202         std::swap(atoms_[1], atoms_[2]);
 203     }
 204 }
 205
 206 void InteractionOfType::sortAtomIds()
 207 {
 208     if (isBond())
 209     {
 210         sortBondAtomIds();
 211     }
 212     else if (isAngle())
 213     {
 214         sortAngleAtomIds();
 215     }
 216     else if (isDihedral())
 217     {
 218         sortDihedralAtomIds();
 219     }
 220     else
 221     {
 222         GMX_THROW(gmx::InternalError(
 223                 "Can not sort parameters other than bonds, angles or dihedrals"));
 224     }
 225 };
 226
 227 void InteractionOfType::setForceParameter(int pos, real value)
 228 {
 229     GMX_RELEASE_ASSERT(pos < MAXFORCEPARAM,
 230                        "Can't set parameter beyond the maximum number of parameters");
 231     forceParam_[pos] = value;
 232 }
 233
 234 void MoleculeInformation::initMolInfo()
 235 {
 236     init_block(&mols);
 237     init_blocka(&excls);
 238     init_t_atoms(&atoms, 0, FALSE);
 239 }
 240
 241 void MoleculeInformation::partialCleanUp()
 242 {
 243     done_block(&mols);
 244 }
 245
 246 void MoleculeInformation::fullCleanUp()
 247 {
 248     done_atom(&atoms);
 249     done_block(&mols);
 250 }
 251
 252 static int rm_interactions(int ifunc, gmx::ArrayRef<MoleculeInformation> mols)
 253 {
 254     int n = 0;
 255     /* For all the molecule types */
 256     for (auto& mol : mols)
 257     {
 258         n += mol.interactions[ifunc].size();
 259         mol.interactions[ifunc].interactionTypes.clear();
 260     }
 261     return n;
 262 }
 263
 264 static int check_atom_names(const char* fn1, const char* fn2, gmx_mtop_t* mtop, const t_atoms* at)
 265 {
 266     int      m, i, j, nmismatch;
 267     t_atoms* tat;
 268 #define MAXMISMATCH 20
 269
 270     if (mtop->natoms != at->nr)
 271     {
 272         gmx_incons("comparing atom names");
 273     }
 274
 275     nmismatch = 0;
 276     i         = 0;
 277     for (const gmx_molblock_t& molb : mtop->molblock)
 278     {
 279         tat = &mtop->moltype[molb.type].atoms;
 280         for (m = 0; m < molb.nmol; m++)
 281         {
 282             for (j = 0; j < tat->nr; j++)
 283             {
 284                 if (strcmp(*(tat->atomname[j]), *(at->atomname[i])) != 0)
 285                 {
 286                     if (nmismatch < MAXMISMATCH)
 287                     {
 288                         fprintf(stderr,
 289                                 "Warning: atom name %d in %s and %s does not match (%s - %s)\n",
 290                                 i + 1, fn1, fn2, *(tat->atomname[j]), *(at->atomname[i]));
 291                     }
 292                     else if (nmismatch == MAXMISMATCH)
 293                     {
 294                         fprintf(stderr, "(more than %d non-matching atom names)\n", MAXMISMATCH);
 295                     }
 296                     nmismatch++;
 297                 }
 298                 i++;
 299             }
 300         }
 301     }
 302
 303     return nmismatch;
 304 }
 305
 306 static void check_bonds_timestep(const gmx_mtop_t* mtop, double dt, warninp* wi)
 307 {
 308     /* This check is not intended to ensure accurate integration,
 309      * rather it is to signal mistakes in the mdp settings.
 310      * A common mistake is to forget to turn on constraints
 311      * for MD after energy minimization with flexible bonds.
 312      * This check can also detect too large time steps for flexible water
 313      * models, but such errors will often be masked by the constraints
 314      * mdp options, which turns flexible water into water with bond constraints,
 315      * but without an angle constraint. Unfortunately such incorrect use
 316      * of water models can not easily be detected without checking
 317      * for specific model names.
 318      *
 319      * The stability limit of leap-frog or velocity verlet is 4.44 steps
 320      * per oscillational period.
 321      * But accurate bonds distributions are lost far before that limit.
 322      * To allow relatively common schemes (although not common with Gromacs)
 323      * of dt=1 fs without constraints and dt=2 fs with only H-bond constraints
 324      * we set the note limit to 10.
 325      */
 326     int  min_steps_warn = 5;
 327     int  min_steps_note = 10;
 328     int  ftype;
 329     int  i, a1, a2, w_a1, w_a2, j;
 330     real twopi2, limit2, fc, re, m1, m2, period2, w_period2;
 331     bool bFound, bWater, bWarn;
 332     char warn_buf[STRLEN];
 333
 334     /* Get the interaction parameters */
 335     gmx::ArrayRef<const t_iparams> ip = mtop->ffparams.iparams;
 336
 337     twopi2 = gmx::square(2 * M_PI);
 338
 339     limit2 = gmx::square(min_steps_note * dt);
 340
 341     w_a1 = w_a2 = -1;
 342     w_period2   = -1.0;
 343
 344     const gmx_moltype_t* w_moltype = nullptr;
 345     for (const gmx_moltype_t& moltype : mtop->moltype)
 346     {
 347         const t_atom*           atom  = moltype.atoms.atom;
 348         const InteractionLists& ilist = moltype.ilist;
 349         const InteractionList&  ilc   = ilist[F_CONSTR];
 350         const InteractionList&  ils   = ilist[F_SETTLE];
 351         for (ftype = 0; ftype < F_NRE; ftype++)
 352         {
 353             if (!(ftype == F_BONDS || ftype == F_G96BONDS || ftype == F_HARMONIC))
 354             {
 355                 continue;
 356             }
 357
 358             const InteractionList& ilb = ilist[ftype];
 359             for (i = 0; i < ilb.size(); i += 3)
 360             {
 361                 fc = ip[ilb.iatoms[i]].harmonic.krA;
 362                 re = ip[ilb.iatoms[i]].harmonic.rA;
 363                 if (ftype == F_G96BONDS)
 364                 {
 365                     /* Convert squared sqaure fc to harmonic fc */
 366                     fc = 2 * fc * re;
 367                 }
 368                 a1 = ilb.iatoms[i + 1];
 369                 a2 = ilb.iatoms[i + 2];
 370                 m1 = atom[a1].m;
 371                 m2 = atom[a2].m;
 372                 if (fc > 0 && m1 > 0 && m2 > 0)
 373                 {
 374                     period2 = twopi2 * m1 * m2 / ((m1 + m2) * fc);
 375                 }
 376                 else
 377                 {
 378                     period2 = GMX_FLOAT_MAX;
 379                 }
 380                 if (debug)
 381                 {
 382                     fprintf(debug, "fc %g m1 %g m2 %g period %g\n", fc, m1, m2, std::sqrt(period2));
 383                 }
 384                 if (period2 < limit2)
 385                 {
 386                     bFound = false;
 387                     for (j = 0; j < ilc.size(); j += 3)
 388                     {
 389                         if ((ilc.iatoms[j + 1] == a1 && ilc.iatoms[j + 2] == a2)
 390                             || (ilc.iatoms[j + 1] == a2 && ilc.iatoms[j + 2] == a1))
 391                         {
 392                             bFound = true;
 393                         }
 394                     }
 395                     for (j = 0; j < ils.size(); j += 4)
 396                     {
 397                         if ((a1 == ils.iatoms[j + 1] || a1 == ils.iatoms[j + 2] || a1 == ils.iatoms[j + 3])
 398                             && (a2 == ils.iatoms[j + 1] || a2 == ils.iatoms[j + 2]
 399                                 || a2 == ils.iatoms[j + 3]))
 400                         {
 401                             bFound = true;
 402                         }
 403                     }
 404                     if (!bFound && (w_moltype == nullptr || period2 < w_period2))
 405                     {
 406                         w_moltype = &moltype;
 407                         w_a1      = a1;
 408                         w_a2      = a2;
 409                         w_period2 = period2;
 410                     }
 411                 }
 412             }
 413         }
 414     }
 415
 416     if (w_moltype != nullptr)
 417     {
 418         bWarn = (w_period2 < gmx::square(min_steps_warn * dt));
 419         /* A check that would recognize most water models */
 420         bWater = ((*w_moltype->atoms.atomname[0])[0] == 'O' && w_moltype->atoms.nr <= 5);
 421         sprintf(warn_buf,
 422                 "The bond in molecule-type %s between atoms %d %s and %d %s has an estimated "
 423                 "oscillational period of %.1e ps, which is less than %d times the time step of "
 424                 "%.1e ps.\n"
 425                 "%s",
 426                 *w_moltype->name, w_a1 + 1, *w_moltype->atoms.atomname[w_a1], w_a2 + 1,
 427                 *w_moltype->atoms.atomname[w_a2], std::sqrt(w_period2),
 428                 bWarn ? min_steps_warn : min_steps_note, dt,
 429                 bWater ? "Maybe you asked for fexible water."
 430                        : "Maybe you forgot to change the constraints mdp option.");
 431         if (bWarn)
 432         {
 433             warning(wi, warn_buf);
 434         }
 435         else
 436         {
 437             warning_note(wi, warn_buf);
 438         }
 439     }
 440 }
 441
 442 static void check_vel(gmx_mtop_t* mtop, rvec v[])
 443 {
 444     for (const AtomProxy atomP : AtomRange(*mtop))
 445     {
 446         const t_atom& local = atomP.atom();
 447         int           i     = atomP.globalAtomNumber();
 448         if (local.ptype == eptShell || local.ptype == eptBond || local.ptype == eptVSite)
 449         {
 450             clear_rvec(v[i]);
 451         }
 452     }
 453 }
 454
 455 static void check_shells_inputrec(gmx_mtop_t* mtop, t_inputrec* ir, warninp* wi)
 456 {
 457     int  nshells = 0;
 458     char warn_buf[STRLEN];
 459
 460     for (const AtomProxy atomP : AtomRange(*mtop))
 461     {
 462         const t_atom& local = atomP.atom();
 463         if (local.ptype == eptShell || local.ptype == eptBond)
 464         {
 465             nshells++;
 466         }
 467     }
 468     if ((nshells > 0) && (ir->nstcalcenergy != 1))
 469     {
 470         set_warning_line(wi, "unknown", -1);
 471         snprintf(warn_buf, STRLEN, "There are %d shells, changing nstcalcenergy from %d to 1",
 472                  nshells, ir->nstcalcenergy);
 473         ir->nstcalcenergy = 1;
 474         warning(wi, warn_buf);
 475     }
 476 }
 477
 478 /* TODO Decide whether this function can be consolidated with
 479  * gmx_mtop_ftype_count */
 480 static int nint_ftype(gmx_mtop_t* mtop, gmx::ArrayRef<const MoleculeInformation> mi, int ftype)
 481 {
 482     int nint = 0;
 483     for (const gmx_molblock_t& molb : mtop->molblock)
 484     {
 485         nint += molb.nmol * mi[molb.type].interactions[ftype].size();
 486     }
 487
 488     return nint;
 489 }
 490
 491 /* This routine reorders the molecule type array
 492  * in the order of use in the molblocks,
 493  * unused molecule types are deleted.
 494  */
 495 static void renumber_moltypes(gmx_mtop_t* sys, std::vector<MoleculeInformation>* molinfo)
 496 {
 497
 498     std::vector<int> order;
 499     for (gmx_molblock_t& molblock : sys->molblock)
 500     {
 501         const auto found = std::find_if(order.begin(), order.end(), [&molblock](const auto& entry) {
 502             return molblock.type == entry;
 503         });
 504         if (found == order.end())
 505         {
 506             /* This type did not occur yet, add it */
 507             order.push_back(molblock.type);
 508             molblock.type = order.size() - 1;
 509         }
 510         else
 511         {
 512             molblock.type = std::distance(order.begin(), found);
 513         }
 514     }
 515
 516     /* We still need to reorder the molinfo structs */
 517     std::vector<MoleculeInformation> minew(order.size());
 518     int                              index = 0;
 519     for (auto& mi : *molinfo)
 520     {
 521         const auto found = std::find(order.begin(), order.end(), index);
 522         if (found != order.end())
 523         {
 524             int pos    = std::distance(order.begin(), found);
 525             minew[pos] = mi;
 526         }
 527         else
 528         {
 529             // Need to manually clean up memory ....
 530             mi.fullCleanUp();
 531         }
 532         index++;
 533     }
 534
 535     *molinfo = minew;
 536 }
 537
 538 static void molinfo2mtop(gmx::ArrayRef<const MoleculeInformation> mi, gmx_mtop_t* mtop)
 539 {
 540     mtop->moltype.resize(mi.size());
 541     int pos = 0;
 542     for (const auto& mol : mi)
 543     {
 544         gmx_moltype_t& molt = mtop->moltype[pos];
 545         molt.name           = mol.name;
 546         molt.atoms          = mol.atoms;
 547         /* ilists are copied later */
 548         molt.excls = mol.excls;
 549         pos++;
 550     }
 551 }
 552
 553 static void new_status(const char*                           topfile,
 554                        const char*                           topppfile,
 555                        const char*                           confin,
 556                        t_gromppopts*                         opts,
 557                        t_inputrec*                           ir,
 558                        gmx_bool                              bZero,
 559                        bool                                  bGenVel,
 560                        bool                                  bVerbose,
 561                        t_state*                              state,
 562                        PreprocessingAtomTypes*               atypes,
 563                        gmx_mtop_t*                           sys,
 564                        std::vector<MoleculeInformation>*     mi,
 565                        std::unique_ptr<MoleculeInformation>* intermolecular_interactions,
 566                        gmx::ArrayRef<InteractionsOfType>     interactions,
 567                        int*                                  comb,
 568                        double*                               reppow,
 569                        real*                                 fudgeQQ,
 570                        gmx_bool                              bMorse,
 571                        warninp*                              wi)
 572 {
 573     std::vector<gmx_molblock_t> molblock;
 574     int                         i, nmismatch;
 575     bool                        ffParametrizedWithHBondConstraints;
 576     char                        buf[STRLEN];
 577     char                        warn_buf[STRLEN];
 578
 579     /* TOPOLOGY processing */
 580     sys->name = do_top(bVerbose, topfile, topppfile, opts, bZero, &(sys->symtab), interactions,
 581                        comb, reppow, fudgeQQ, atypes, mi, intermolecular_interactions, ir,
 582                        &molblock, &ffParametrizedWithHBondConstraints, wi);
 583
 584     sys->molblock.clear();
 585
 586     sys->natoms = 0;
 587     for (const gmx_molblock_t& molb : molblock)
 588     {
 589         if (!sys->molblock.empty() && molb.type == sys->molblock.back().type)
 590         {
 591             /* Merge consecutive blocks with the same molecule type */
 592             sys->molblock.back().nmol += molb.nmol;
 593         }
 594         else if (molb.nmol > 0)
 595         {
 596             /* Add a new molblock to the topology */
 597             sys->molblock.push_back(molb);
 598         }
 599         sys->natoms += molb.nmol * (*mi)[sys->molblock.back().type].atoms.nr;
 600     }
 601     if (sys->molblock.empty())
 602     {
 603         gmx_fatal(FARGS, "No molecules were defined in the system");
 604     }
 605
 606     renumber_moltypes(sys, mi);
 607
 608     if (bMorse)
 609     {
 610         convert_harmonics(*mi, atypes);
 611     }
 612
 613     if (ir->eDisre == edrNone)
 614     {
 615         i = rm_interactions(F_DISRES, *mi);
 616         if (i > 0)
 617         {
 618             set_warning_line(wi, "unknown", -1);
 619             sprintf(warn_buf, "disre = no, removed %d distance restraints", i);
 620             warning_note(wi, warn_buf);
 621         }
 622     }
 623     if (!opts->bOrire)
 624     {
 625         i = rm_interactions(F_ORIRES, *mi);
 626         if (i > 0)
 627         {
 628             set_warning_line(wi, "unknown", -1);
 629             sprintf(warn_buf, "orire = no, removed %d orientation restraints", i);
 630             warning_note(wi, warn_buf);
 631         }
 632     }
 633
 634     /* Copy structures from msys to sys */
 635     molinfo2mtop(*mi, sys);
 636
 637     gmx_mtop_finalize(sys);
 638
 639     /* Discourage using the, previously common, setup of applying constraints
 640      * to all bonds with force fields that have been parametrized with H-bond
 641      * constraints only. Do not print note with large timesteps or vsites.
 642      */
 643     if (opts->nshake == eshALLBONDS && ffParametrizedWithHBondConstraints && ir->delta_t < 0.0026
 644         && gmx_mtop_ftype_count(sys, F_VSITE3FD) == 0)
 645     {
 646         set_warning_line(wi, "unknown", -1);
 647         warning_note(wi,
 648                      "You are using constraints on all bonds, whereas the forcefield "
 649                      "has been parametrized only with constraints involving hydrogen atoms. "
 650                      "We suggest using constraints = h-bonds instead, this will also improve "
 651                      "performance.");
 652     }
 653
 654     /* COORDINATE file processing */
 655     if (bVerbose)
 656     {
 657         fprintf(stderr, "processing coordinates...\n");
 658     }
 659
 660     t_topology* conftop;
 661     rvec*       x = nullptr;
 662     rvec*       v = nullptr;
 663     snew(conftop, 1);
 664     read_tps_conf(confin, conftop, nullptr, &x, &v, state->box, FALSE);
 665     state->natoms = conftop->atoms.nr;
 666     if (state->natoms != sys->natoms)
 667     {
 668         gmx_fatal(FARGS,
 669                   "number of coordinates in coordinate file (%s, %d)\n"
 670                   "             does not match topology (%s, %d)",
 671                   confin, state->natoms, topfile, sys->natoms);
 672     }
 673     /* It would be nice to get rid of the copies below, but we don't know
 674      * a priori if the number of atoms in confin matches what we expect.
 675      */
 676     state->flags |= (1 << estX);
 677     if (v != nullptr)
 678     {
 679         state->flags |= (1 << estV);
 680     }
 681     state_change_natoms(state, state->natoms);
 682     std::copy(x, x + state->natoms, state->x.data());
 683     sfree(x);
 684     if (v != nullptr)
 685     {
 686         std::copy(v, v + state->natoms, state->v.data());
 687         sfree(v);
 688     }
 689     /* This call fixes the box shape for runs with pressure scaling */
 690     set_box_rel(ir, state);
 691
 692     nmismatch = check_atom_names(topfile, confin, sys, &conftop->atoms);
 693     done_top(conftop);
 694     sfree(conftop);
 695
 696     if (nmismatch)
 697     {
 698         sprintf(buf,
 699                 "%d non-matching atom name%s\n"
 700                 "atom names from %s will be used\n"
 701                 "atom names from %s will be ignored\n",
 702                 nmismatch, (nmismatch == 1) ? "" : "s", topfile, confin);
 703         warning(wi, buf);
 704     }
 705
 706     /* Do more checks, mostly related to constraints */
 707     if (bVerbose)
 708     {
 709         fprintf(stderr, "double-checking input for internal consistency...\n");
 710     }
 711     {
 712         bool bHasNormalConstraints =
 713                 0 < (nint_ftype(sys, *mi, F_CONSTR) + nint_ftype(sys, *mi, F_CONSTRNC));
 714         bool bHasAnyConstraints = bHasNormalConstraints || 0 < nint_ftype(sys, *mi, F_SETTLE);
 715         double_check(ir, state->box, bHasNormalConstraints, bHasAnyConstraints, wi);
 716     }
 717
 718     if (bGenVel)
 719     {
 720         real* mass;
 721
 722         snew(mass, state->natoms);
 723         for (const AtomProxy atomP : AtomRange(*sys))
 724         {
 725             const t_atom& local = atomP.atom();
 726             int           i     = atomP.globalAtomNumber();
 727             mass[i]             = local.m;
 728         }
 729
 730         if (opts->seed == -1)
 731         {
 732             opts->seed = static_cast<int>(gmx::makeRandomSeed());
 733             fprintf(stderr, "Setting gen_seed to %d\n", opts->seed);
 734         }
 735         state->flags |= (1 << estV);
 736         maxwell_speed(opts->tempi, opts->seed, sys, state->v.rvec_array());
 737
 738         stop_cm(stdout, state->natoms, mass, state->x.rvec_array(), state->v.rvec_array());
 739         sfree(mass);
 740     }
 741 }
 742
 743 static void copy_state(const char* slog, t_trxframe* fr, bool bReadVel, t_state* state, double* use_time)
 744 {
 745     if (fr->not_ok & FRAME_NOT_OK)
 746     {
 747         gmx_fatal(FARGS, "Can not start from an incomplete frame");
 748     }
 749     if (!fr->bX)
 750     {
 751         gmx_fatal(FARGS, "Did not find a frame with coordinates in file %s", slog);
 752     }
 753
 754     std::copy(fr->x, fr->x + state->natoms, state->x.data());
 755     if (bReadVel)
 756     {
 757         if (!fr->bV)
 758         {
 759             gmx_incons("Trajecory frame unexpectedly does not contain velocities");
 760         }
 761         std::copy(fr->v, fr->v + state->natoms, state->v.data());
 762     }
 763     if (fr->bBox)
 764     {
 765         copy_mat(fr->box, state->box);
 766     }
 767
 768     *use_time = fr->time;
 769 }
 770
 771 static void cont_status(const char*             slog,
 772                         const char*             ener,
 773                         bool                    bNeedVel,
 774                         bool                    bGenVel,
 775                         real                    fr_time,
 776                         t_inputrec*             ir,
 777                         t_state*                state,
 778                         gmx_mtop_t*             sys,
 779                         const gmx_output_env_t* oenv)
 780 /* If fr_time == -1 read the last frame available which is complete */
 781 {
 782     bool         bReadVel;
 783     t_trxframe   fr;
 784     t_trxstatus* fp;
 785     double       use_time;
 786
 787     bReadVel = (bNeedVel && !bGenVel);
 788
 789     fprintf(stderr, "Reading Coordinates%s and Box size from old trajectory\n",
 790             bReadVel ? ", Velocities" : "");
 791     if (fr_time == -1)
 792     {
 793         fprintf(stderr, "Will read whole trajectory\n");
 794     }
 795     else
 796     {
 797         fprintf(stderr, "Will read till time %g\n", fr_time);
 798     }
 799     if (!bReadVel)
 800     {
 801         if (bGenVel)
 802         {
 803             fprintf(stderr,
 804                     "Velocities generated: "
 805                     "ignoring velocities in input trajectory\n");
 806         }
 807         read_first_frame(oenv, &fp, slog, &fr, TRX_NEED_X);
 808     }
 809     else
 810     {
 811         read_first_frame(oenv, &fp, slog, &fr, TRX_NEED_X | TRX_NEED_V);
 812
 813         if (!fr.bV)
 814         {
 815             fprintf(stderr,
 816                     "\n"
 817                     "WARNING: Did not find a frame with velocities in file %s,\n"
 818                     "         all velocities will be set to zero!\n\n",
 819                     slog);
 820             for (auto& vi : makeArrayRef(state->v))
 821             {
 822                 vi = { 0, 0, 0 };
 823             }
 824             close_trx(fp);
 825             /* Search for a frame without velocities */
 826             bReadVel = false;
 827             read_first_frame(oenv, &fp, slog, &fr, TRX_NEED_X);
 828         }
 829     }
 830
 831     state->natoms = fr.natoms;
 832
 833     if (sys->natoms != state->natoms)
 834     {
 835         gmx_fatal(FARGS,
 836                   "Number of atoms in Topology "
 837                   "is not the same as in Trajectory");
 838     }
 839     copy_state(slog, &fr, bReadVel, state, &use_time);
 840
 841     /* Find the appropriate frame */
 842     while ((fr_time == -1 || fr.time < fr_time) && read_next_frame(oenv, fp, &fr))
 843     {
 844         copy_state(slog, &fr, bReadVel, state, &use_time);
 845     }
 846
 847     close_trx(fp);
 848
 849     /* Set the relative box lengths for preserving the box shape.
 850      * Note that this call can lead to differences in the last bit
 851      * with respect to using gmx convert-tpr to create a [REF].tpx[ref] file.
 852      */
 853     set_box_rel(ir, state);
 854
 855     fprintf(stderr, "Using frame at t = %g ps\n", use_time);
 856     fprintf(stderr, "Starting time for run is %g ps\n", ir->init_t);
 857
 858     if ((ir->epc != epcNO || ir->etc == etcNOSEHOOVER) && ener)
 859     {
 860         get_enx_state(ener, use_time, sys->groups, ir, state);
 861         preserve_box_shape(ir, state->box_rel, state->boxv);
 862     }
 863 }
 864
 865 static void read_posres(gmx_mtop_t*                              mtop,
 866                         gmx::ArrayRef<const MoleculeInformation> molinfo,
 867                         gmx_bool                                 bTopB,
 868                         const char*                              fn,
 869                         int                                      rc_scaling,
 870                         int                                      ePBC,
 871                         rvec                                     com,
 872                         warninp*                                 wi)
 873 {
 874     gmx_bool*   hadAtom;
 875     rvec *      x, *v;
 876     dvec        sum;
 877     double      totmass;
 878     t_topology* top;
 879     matrix      box, invbox;
 880     int         natoms, npbcdim = 0;
 881     char        warn_buf[STRLEN];
 882     int         a, nat_molb;
 883     t_atom*     atom;
 884
 885     snew(top, 1);
 886     read_tps_conf(fn, top, nullptr, &x, &v, box, FALSE);
 887     natoms = top->atoms.nr;
 888     done_top(top);
 889     sfree(top);
 890     if (natoms != mtop->natoms)
 891     {
 892         sprintf(warn_buf,
 893                 "The number of atoms in %s (%d) does not match the number of atoms in the topology "
 894                 "(%d). Will assume that the first %d atoms in the topology and %s match.",
 895                 fn, natoms, mtop->natoms, std::min(mtop->natoms, natoms), fn);
 896         warning(wi, warn_buf);
 897     }
 898
 899     npbcdim = ePBC2npbcdim(ePBC);
 900     GMX_RELEASE_ASSERT(npbcdim <= DIM, "Invalid npbcdim");
 901     clear_rvec(com);
 902     if (rc_scaling != erscNO)
 903     {
 904         copy_mat(box, invbox);
 905         for (int j = npbcdim; j < DIM; j++)
 906         {
 907             clear_rvec(invbox[j]);
 908             invbox[j][j] = 1;
 909         }
 910         gmx::invertBoxMatrix(invbox, invbox);
 911     }
 912
 913     /* Copy the reference coordinates to mtop */
 914     clear_dvec(sum);
 915     totmass = 0;
 916     a       = 0;
 917     snew(hadAtom, natoms);
 918     for (gmx_molblock_t& molb : mtop->molblock)
 919     {
 920         nat_molb                       = molb.nmol * mtop->moltype[molb.type].atoms.nr;
 921         const InteractionsOfType* pr   = &(molinfo[molb.type].interactions[F_POSRES]);
 922         const InteractionsOfType* prfb = &(molinfo[molb.type].interactions[F_FBPOSRES]);
 923         if (pr->size() > 0 || prfb->size() > 0)
 924         {
 925             atom = mtop->moltype[molb.type].atoms.atom;
 926             for (const auto& restraint : pr->interactionTypes)
 927             {
 928                 int ai = restraint.ai();
 929                 if (ai >= natoms)
 930                 {
 931                     gmx_fatal(FARGS,
 932                               "Position restraint atom index (%d) in moltype '%s' is larger than "
 933                               "number of atoms in %s (%d).\n",
 934                               ai + 1, *molinfo[molb.type].name, fn, natoms);
 935                 }
 936                 hadAtom[ai] = TRUE;
 937                 if (rc_scaling == erscCOM)
 938                 {
 939                     /* Determine the center of mass of the posres reference coordinates */
 940                     for (int j = 0; j < npbcdim; j++)
 941                     {
 942                         sum[j] += atom[ai].m * x[a + ai][j];
 943                     }
 944                     totmass += atom[ai].m;
 945                 }
 946             }
 947             /* Same for flat-bottomed posres, but do not count an atom twice for COM */
 948             for (const auto& restraint : prfb->interactionTypes)
 949             {
 950                 int ai = restraint.ai();
 951                 if (ai >= natoms)
 952                 {
 953                     gmx_fatal(FARGS,
 954                               "Position restraint atom index (%d) in moltype '%s' is larger than "
 955                               "number of atoms in %s (%d).\n",
 956                               ai + 1, *molinfo[molb.type].name, fn, natoms);
 957                 }
 958                 if (rc_scaling == erscCOM && !hadAtom[ai])
 959                 {
 960                     /* Determine the center of mass of the posres reference coordinates */
 961                     for (int j = 0; j < npbcdim; j++)
 962                     {
 963                         sum[j] += atom[ai].m * x[a + ai][j];
 964                     }
 965                     totmass += atom[ai].m;
 966                 }
 967             }
 968             if (!bTopB)
 969             {
 970                 molb.posres_xA.resize(nat_molb);
 971                 for (int i = 0; i < nat_molb; i++)
 972                 {
 973                     copy_rvec(x[a + i], molb.posres_xA[i]);
 974                 }
 975             }
 976             else
 977             {
 978                 molb.posres_xB.resize(nat_molb);
 979                 for (int i = 0; i < nat_molb; i++)
 980                 {
 981                     copy_rvec(x[a + i], molb.posres_xB[i]);
 982                 }
 983             }
 984         }
 985         a += nat_molb;
 986     }
 987     if (rc_scaling == erscCOM)
 988     {
 989         if (totmass == 0)
 990         {
 991             gmx_fatal(FARGS, "The total mass of the position restraint atoms is 0");
 992         }
 993         for (int j = 0; j < npbcdim; j++)
 994         {
 995             com[j] = sum[j] / totmass;
 996         }
 997         fprintf(stderr,
 998                 "The center of mass of the position restraint coord's is %6.3f %6.3f %6.3f\n",
 999                 com[XX], com[YY], com[ZZ]);
1000     }
1001
1002     if (rc_scaling != erscNO)
1003     {
1004         GMX_ASSERT(npbcdim <= DIM, "Only DIM dimensions can have PBC");
1005
1006         for (gmx_molblock_t& molb : mtop->molblock)
1007         {
1008             nat_molb = molb.nmol * mtop->moltype[molb.type].atoms.nr;
1009             if (!molb.posres_xA.empty() || !molb.posres_xB.empty())
1010             {
1011                 std::vector<gmx::RVec>& xp = (!bTopB ? molb.posres_xA : molb.posres_xB);
1012                 for (int i = 0; i < nat_molb; i++)
1013                 {
1014                     for (int j = 0; j < npbcdim; j++)
1015                     {
1016                         if (rc_scaling == erscALL)
1017                         {
1018                             /* Convert from Cartesian to crystal coordinates */
1019                             xp[i][j] *= invbox[j][j];
1020                             for (int k = j + 1; k < npbcdim; k++)
1021                             {
1022                                 xp[i][j] += invbox[k][j] * xp[i][k];
1023                             }
1024                         }
1025                         else if (rc_scaling == erscCOM)
1026                         {
1027                             /* Subtract the center of mass */
1028                             xp[i][j] -= com[j];
1029                         }
1030                     }
1031                 }
1032             }
1033         }
1034
1035         if (rc_scaling == erscCOM)
1036         {
1037             /* Convert the COM from Cartesian to crystal coordinates */
1038             for (int j = 0; j < npbcdim; j++)
1039             {
1040                 com[j] *= invbox[j][j];
1041                 for (int k = j + 1; k < npbcdim; k++)
1042                 {
1043                     com[j] += invbox[k][j] * com[k];
1044                 }
1045             }
1046         }
1047     }
1048
1049     sfree(x);
1050     sfree(v);
1051     sfree(hadAtom);
1052 }
1053
1054 static void gen_posres(gmx_mtop_t*                              mtop,
1055                        gmx::ArrayRef<const MoleculeInformation> mi,
1056                        const char*                              fnA,
1057                        const char*                              fnB,
1058                        int                                      rc_scaling,
1059                        int                                      ePBC,
1060                        rvec                                     com,
1061                        rvec                                     comB,
1062                        warninp*                                 wi)
1063 {
1064     read_posres(mtop, mi, FALSE, fnA, rc_scaling, ePBC, com, wi);
1065     /* It is safer to simply read the b-state posres rather than trying
1066      * to be smart and copy the positions.
1067      */
1068     read_posres(mtop, mi, TRUE, fnB, rc_scaling, ePBC, comB, wi);
1069 }
1070
1071 static void set_wall_atomtype(PreprocessingAtomTypes* at, t_gromppopts* opts, t_inputrec* ir, warninp* wi)
1072 {
1073     int  i;
1074     char warn_buf[STRLEN];
1075
1076     if (ir->nwall > 0)
1077     {
1078         fprintf(stderr, "Searching the wall atom type(s)\n");
1079     }
1080     for (i = 0; i < ir->nwall; i++)
1081     {
1082         ir->wall_atomtype[i] = at->atomTypeFromName(opts->wall_atomtype[i]);
1083         if (ir->wall_atomtype[i] == NOTSET)
1084         {
1085             sprintf(warn_buf, "Specified wall atom type %s is not defined", opts->wall_atomtype[i]);
1086             warning_error(wi, warn_buf);
1087         }
1088     }
1089 }
1090
1091 static int nrdf_internal(const t_atoms* atoms)
1092 {
1093     int i, nmass, nrdf;
1094
1095     nmass = 0;
1096     for (i = 0; i < atoms->nr; i++)
1097     {
1098         /* Vsite ptype might not be set here yet, so also check the mass */
1099         if ((atoms->atom[i].ptype == eptAtom || atoms->atom[i].ptype == eptNucleus) && atoms->atom[i].m > 0)
1100         {
1101             nmass++;
1102         }
1103     }
1104     switch (nmass)
1105     {
1106         case 0: nrdf = 0; break;
1107         case 1: nrdf = 0; break;
1108         case 2: nrdf = 1; break;
1109         default: nrdf = nmass * 3 - 6; break;
1110     }
1111
1112     return nrdf;
1113 }
1114
1115 static void spline1d(double dx, const double* y, int n, double* u, double* y2)
1116 {
1117     int    i;
1118     double p, q;
1119
1120     y2[0] = 0.0;
1121     u[0]  = 0.0;
1122
1123     for (i = 1; i < n - 1; i++)
1124     {
1125         p     = 0.5 * y2[i - 1] + 2.0;
1126         y2[i] = -0.5 / p;
1127         q     = (y[i + 1] - 2.0 * y[i] + y[i - 1]) / dx;
1128         u[i]  = (3.0 * q / dx - 0.5 * u[i - 1]) / p;
1129     }
1130
1131     y2[n - 1] = 0.0;
1132
1133     for (i = n - 2; i >= 0; i--)
1134     {
1135         y2[i] = y2[i] * y2[i + 1] + u[i];
1136     }
1137 }
1138
1139
1140 static void
1141 interpolate1d(double xmin, double dx, const double* ya, const double* y2a, double x, double* y, double* y1)
1142 {
1143     int    ix;
1144     double a, b;
1145
1146     ix = static_cast<int>((x - xmin) / dx);
1147
1148     a = (xmin + (ix + 1) * dx - x) / dx;
1149     b = (x - xmin - ix * dx) / dx;
1150
1151     *y = a * ya[ix] + b * ya[ix + 1]
1152          + ((a * a * a - a) * y2a[ix] + (b * b * b - b) * y2a[ix + 1]) * (dx * dx) / 6.0;
1153     *y1 = (ya[ix + 1] - ya[ix]) / dx - (3.0 * a * a - 1.0) / 6.0 * dx * y2a[ix]
1154           + (3.0 * b * b - 1.0) / 6.0 * dx * y2a[ix + 1];
1155 }
1156
1157
1158 static void setup_cmap(int grid_spacing, int nc, gmx::ArrayRef<const real> grid, gmx_cmap_t* cmap_grid)
1159 {
1160     int    i, j, k, ii, jj, kk, idx;
1161     int    offset;
1162     double dx, xmin, v, v1, v2, v12;
1163     double phi, psi;
1164
1165     std::vector<double> tmp_u(2 * grid_spacing, 0.0);
1166     std::vector<double> tmp_u2(2 * grid_spacing, 0.0);
1167     std::vector<double> tmp_yy(2 * grid_spacing, 0.0);
1168     std::vector<double> tmp_y1(2 * grid_spacing, 0.0);
1169     std::vector<double> tmp_t2(2 * grid_spacing * 2 * grid_spacing, 0.0);
1170     std::vector<double> tmp_grid(2 * grid_spacing * 2 * grid_spacing, 0.0);
1171
1172     dx   = 360.0 / grid_spacing;
1173     xmin = -180.0 - dx * grid_spacing / 2;
1174
1175     for (kk = 0; kk < nc; kk++)
1176     {
1177         /* Compute an offset depending on which cmap we are using
1178          * Offset will be the map number multiplied with the
1179          * grid_spacing * grid_spacing * 2
1180          */
1181         offset = kk * grid_spacing * grid_spacing * 2;
1182
1183         for (i = 0; i < 2 * grid_spacing; i++)
1184         {
1185             ii = (i + grid_spacing - grid_spacing / 2) % grid_spacing;
1186
1187             for (j = 0; j < 2 * grid_spacing; j++)
1188             {
1189                 jj = (j + grid_spacing - grid_spacing / 2) % grid_spacing;
1190                 tmp_grid[i * grid_spacing * 2 + j] = grid[offset + ii * grid_spacing + jj];
1191             }
1192         }
1193
1194         for (i = 0; i < 2 * grid_spacing; i++)
1195         {
1196             spline1d(dx, &(tmp_grid[2 * grid_spacing * i]), 2 * grid_spacing, tmp_u.data(),
1197                      &(tmp_t2[2 * grid_spacing * i]));
1198         }
1199
1200         for (i = grid_spacing / 2; i < grid_spacing + grid_spacing / 2; i++)
1201         {
1202             ii  = i - grid_spacing / 2;
1203             phi = ii * dx - 180.0;
1204
1205             for (j = grid_spacing / 2; j < grid_spacing + grid_spacing / 2; j++)
1206             {
1207                 jj  = j - grid_spacing / 2;
1208                 psi = jj * dx - 180.0;
1209
1210                 for (k = 0; k < 2 * grid_spacing; k++)
1211                 {
1212                     interpolate1d(xmin, dx, &(tmp_grid[2 * grid_spacing * k]),
1213                                   &(tmp_t2[2 * grid_spacing * k]), psi, &tmp_yy[k], &tmp_y1[k]);
1214                 }
1215
1216                 spline1d(dx, tmp_yy.data(), 2 * grid_spacing, tmp_u.data(), tmp_u2.data());
1217                 interpolate1d(xmin, dx, tmp_yy.data(), tmp_u2.data(), phi, &v, &v1);
1218                 spline1d(dx, tmp_y1.data(), 2 * grid_spacing, tmp_u.data(), tmp_u2.data());
1219                 interpolate1d(xmin, dx, tmp_y1.data(), tmp_u2.data(), phi, &v2, &v12);
1220
1221                 idx                                       = ii * grid_spacing + jj;
1222                 cmap_grid->cmapdata[kk].cmap[idx * 4]     = grid[offset + ii * grid_spacing + jj];
1223                 cmap_grid->cmapdata[kk].cmap[idx * 4 + 1] = v1;
1224                 cmap_grid->cmapdata[kk].cmap[idx * 4 + 2] = v2;
1225                 cmap_grid->cmapdata[kk].cmap[idx * 4 + 3] = v12;
1226             }
1227         }
1228     }
1229 }
1230
1231 static void init_cmap_grid(gmx_cmap_t* cmap_grid, int ngrid, int grid_spacing)
1232 {
1233     int i, nelem;
1234
1235     cmap_grid->grid_spacing = grid_spacing;
1236     nelem                   = cmap_grid->grid_spacing * cmap_grid->grid_spacing;
1237
1238     cmap_grid->cmapdata.resize(ngrid);
1239
1240     for (i = 0; i < ngrid; i++)
1241     {
1242         cmap_grid->cmapdata[i].cmap.resize(4 * nelem);
1243     }
1244 }
1245
1246
1247 static int count_constraints(const gmx_mtop_t* mtop, gmx::ArrayRef<const MoleculeInformation> mi, warninp* wi)
1248 {
1249     int  count, count_mol;
1250     char buf[STRLEN];
1251
1252     count = 0;
1253     for (const gmx_molblock_t& molb : mtop->molblock)
1254     {
1255         count_mol                                            = 0;
1256         gmx::ArrayRef<const InteractionsOfType> interactions = mi[molb.type].interactions;
1257
1258         for (int i = 0; i < F_NRE; i++)
1259         {
1260             if (i == F_SETTLE)
1261             {
1262                 count_mol += 3 * interactions[i].size();
1263             }
1264             else if (interaction_function[i].flags & IF_CONSTRAINT)
1265             {
1266                 count_mol += interactions[i].size();
1267             }
1268         }
1269
1270         if (count_mol > nrdf_internal(&mi[molb.type].atoms))
1271         {
1272             sprintf(buf,
1273                     "Molecule type '%s' has %d constraints.\n"
1274                     "For stability and efficiency there should not be more constraints than "
1275                     "internal number of degrees of freedom: %d.\n",
1276                     *mi[molb.type].name, count_mol, nrdf_internal(&mi[molb.type].atoms));
1277             warning(wi, buf);
1278         }
1279         count += molb.nmol * count_mol;
1280     }
1281
1282     return count;
1283 }
1284
1285 static real calc_temp(const gmx_mtop_t* mtop, const t_inputrec* ir, rvec* v)
1286 {
1287     double sum_mv2 = 0;
1288     for (const AtomProxy atomP : AtomRange(*mtop))
1289     {
1290         const t_atom& local = atomP.atom();
1291         int           i     = atomP.globalAtomNumber();
1292         sum_mv2 += local.m * norm2(v[i]);
1293     }
1294
1295     double nrdf = 0;
1296     for (int g = 0; g < ir->opts.ngtc; g++)
1297     {
1298         nrdf += ir->opts.nrdf[g];
1299     }
1300
1301     return sum_mv2 / (nrdf * BOLTZ);
1302 }
1303
1304 static real get_max_reference_temp(const t_inputrec* ir, warninp* wi)
1305 {
1306     real ref_t;
1307     int  i;
1308     bool bNoCoupl;
1309
1310     ref_t    = 0;
1311     bNoCoupl = false;
1312     for (i = 0; i < ir->opts.ngtc; i++)
1313     {
1314         if (ir->opts.tau_t[i] < 0)
1315         {
1316             bNoCoupl = true;
1317         }
1318         else
1319         {
1320             ref_t = std::max(ref_t, ir->opts.ref_t[i]);
1321         }
1322     }
1323
1324     if (bNoCoupl)
1325     {
1326         char buf[STRLEN];
1327
1328         sprintf(buf,
1329                 "Some temperature coupling groups do not use temperature coupling. We will assume "
1330                 "their temperature is not more than %.3f K. If their temperature is higher, the "
1331                 "energy error and the Verlet buffer might be underestimated.",
1332                 ref_t);
1333         warning(wi, buf);
1334     }
1335
1336     return ref_t;
1337 }
1338
1339 /* Checks if there are unbound atoms in moleculetype molt.
1340  * Prints a note for each unbound atoms and a warning if any is present.
1341  */
1342 static void checkForUnboundAtoms(const gmx_moltype_t* molt, gmx_bool bVerbose, warninp* wi)
1343 {
1344     const t_atoms* atoms = &molt->atoms;
1345
1346     if (atoms->nr == 1)
1347     {
1348         /* Only one atom, there can't be unbound atoms */
1349         return;
1350     }
1351
1352     std::vector<int> count(atoms->nr, 0);
1353
1354     for (int ftype = 0; ftype < F_NRE; ftype++)
1355     {
1356         if (((interaction_function[ftype].flags & IF_BOND) && ftype != F_CONNBONDS)
1357             || (interaction_function[ftype].flags & IF_CONSTRAINT) || ftype == F_SETTLE)
1358         {
1359             const InteractionList& il   = molt->ilist[ftype];
1360             const int              nral = NRAL(ftype);
1361
1362             for (int i = 0; i < il.size(); i += 1 + nral)
1363             {
1364                 for (int j = 0; j < nral; j++)
1365                 {
1366                     count[il.iatoms[i + 1 + j]]++;
1367                 }
1368             }
1369         }
1370     }
1371
1372     int numDanglingAtoms = 0;
1373     for (int a = 0; a < atoms->nr; a++)
1374     {
1375         if (atoms->atom[a].ptype != eptVSite && count[a] == 0)
1376         {
1377             if (bVerbose)
1378             {
1379                 fprintf(stderr,
1380                         "\nAtom %d '%s' in moleculetype '%s' is not bound by a potential or "
1381                         "constraint to any other atom in the same moleculetype.\n",
1382                         a + 1, *atoms->atomname[a], *molt->name);
1383             }
1384             numDanglingAtoms++;
1385         }
1386     }
1387
1388     if (numDanglingAtoms > 0)
1389     {
1390         char buf[STRLEN];
1391         sprintf(buf,
1392                 "In moleculetype '%s' %d atoms are not bound by a potential or constraint to any "
1393                 "other atom in the same moleculetype. Although technically this might not cause "
1394                 "issues in a simulation, this often means that the user forgot to add a "
1395                 "bond/potential/constraint or put multiple molecules in the same moleculetype "
1396                 "definition by mistake. Run with -v to get information for each atom.",
1397                 *molt->name, numDanglingAtoms);
1398         warning_note(wi, buf);
1399     }
1400 }
1401
1402 /* Checks all moleculetypes for unbound atoms */
1403 static void checkForUnboundAtoms(const gmx_mtop_t* mtop, gmx_bool bVerbose, warninp* wi)
1404 {
1405     for (const gmx_moltype_t& molt : mtop->moltype)
1406     {
1407         checkForUnboundAtoms(&molt, bVerbose, wi);
1408     }
1409 }
1410
1411 /*! \brief Checks if there are decoupled modes in moleculetype \p molt.
1412  *
1413  * The specific decoupled modes this routine check for are angle modes
1414  * where the two bonds are constrained and the atoms a both ends are only
1415  * involved in a single constraint; the mass of the two atoms needs to
1416  * differ by more than \p massFactorThreshold.
1417  */
1418 static bool haveDecoupledModeInMol(const gmx_moltype_t&           molt,
1419                                    gmx::ArrayRef<const t_iparams> iparams,
1420                                    real                           massFactorThreshold)
1421 {
1422     if (molt.ilist[F_CONSTR].size() == 0 && molt.ilist[F_CONSTRNC].size() == 0)
1423     {
1424         return false;
1425     }
1426
1427     const t_atom* atom = molt.atoms.atom;
1428
1429     t_blocka atomToConstraints =
1430             gmx::make_at2con(molt, iparams, gmx::FlexibleConstraintTreatment::Exclude);
1431
1432     bool haveDecoupledMode = false;
1433     for (int ftype = 0; ftype < F_NRE; ftype++)
1434     {
1435         if (interaction_function[ftype].flags & IF_ATYPE)
1436         {
1437             const int              nral = NRAL(ftype);
1438             const InteractionList& il   = molt.ilist[ftype];
1439             for (int i = 0; i < il.size(); i += 1 + nral)
1440             {
1441                 /* Here we check for the mass difference between the atoms
1442                  * at both ends of the angle, that the atoms at the ends have
1443                  * 1 contraint and the atom in the middle at least 3; we check
1444                  * that the 3 atoms are linked by constraints below.
1445                  * We check for at least three constraints for the middle atom,
1446                  * since with only the two bonds in the angle, we have 3-atom
1447                  * molecule, which has much more thermal exhange in this single
1448                  * angle mode than molecules with more atoms.
1449                  * Note that this check also catches molecules where more atoms
1450                  * are connected to one or more atoms in the angle, but by
1451                  * bond potentials instead of angles. But such cases will not
1452                  * occur in "normal" molecules and it doesn't hurt running
1453                  * those with higher accuracy settings as well.
1454                  */
1455                 int a0 = il.iatoms[1 + i];
1456                 int a1 = il.iatoms[1 + i + 1];
1457                 int a2 = il.iatoms[1 + i + 2];
1458                 if ((atom[a0].m > atom[a2].m * massFactorThreshold || atom[a2].m > atom[a0].m * massFactorThreshold)
1459                     && atomToConstraints.index[a0 + 1] - atomToConstraints.index[a0] == 1
1460                     && atomToConstraints.index[a2 + 1] - atomToConstraints.index[a2] == 1
1461                     && atomToConstraints.index[a1 + 1] - atomToConstraints.index[a1] >= 3)
1462                 {
1463                     int constraint0 = atomToConstraints.a[atomToConstraints.index[a0]];
1464                     int constraint2 = atomToConstraints.a[atomToConstraints.index[a2]];
1465
1466                     bool foundAtom0 = false;
1467                     bool foundAtom2 = false;
1468                     for (int conIndex = atomToConstraints.index[a1];
1469                          conIndex < atomToConstraints.index[a1 + 1]; conIndex++)
1470                     {
1471                         if (atomToConstraints.a[conIndex] == constraint0)
1472                         {
1473                             foundAtom0 = true;
1474                         }
1475                         if (atomToConstraints.a[conIndex] == constraint2)
1476                         {
1477                             foundAtom2 = true;
1478                         }
1479                     }
1480                     if (foundAtom0 && foundAtom2)
1481                     {
1482                         haveDecoupledMode = true;
1483                     }
1484                 }
1485             }
1486         }
1487     }
1488
1489     done_blocka(&atomToConstraints);
1490
1491     return haveDecoupledMode;
1492 }
1493
1494 /*! \brief Checks if the Verlet buffer and constraint accuracy is sufficient for decoupled dynamic modes.
1495  *
1496  * When decoupled modes are present and the accuray in insufficient,
1497  * this routine issues a warning if the accuracy is insufficient.
1498  */
1499 static void checkDecoupledModeAccuracy(const gmx_mtop_t* mtop, const t_inputrec* ir, warninp* wi)
1500 {
1501     /* We only have issues with decoupled modes with normal MD.
1502      * With stochastic dynamics equipartitioning is enforced strongly.
1503      */
1504     if (!EI_MD(ir->eI))
1505     {
1506         return;
1507     }
1508
1509     /* When atoms of very different mass are involved in an angle potential
1510      * and both bonds in the angle are constrained, the dynamic modes in such
1511      * angles have very different periods and significant energy exchange
1512      * takes several nanoseconds. Thus even a small amount of error in
1513      * different algorithms can lead to violation of equipartitioning.
1514      * The parameters below are mainly based on an all-atom chloroform model
1515      * with all bonds constrained. Equipartitioning is off by more than 1%
1516      * (need to run 5-10 ns) when the difference in mass between H and Cl
1517      * is more than a factor 13 and the accuracy is less than the thresholds
1518      * given below. This has been verified on some other molecules.
1519      *
1520      * Note that the buffer and shake parameters have unit length and
1521      * energy/time, respectively, so they will "only" work correctly
1522      * for atomistic force fields using MD units.
1523      */
1524     const real massFactorThreshold      = 13.0;
1525     const real bufferToleranceThreshold = 1e-4;
1526     const int  lincsIterationThreshold  = 2;
1527     const int  lincsOrderThreshold      = 4;
1528     const real shakeToleranceThreshold  = 0.005 * ir->delta_t;
1529
1530     bool lincsWithSufficientTolerance = (ir->eConstrAlg == econtLINCS && ir->nLincsIter >= lincsIterationThreshold
1531                                          && ir->nProjOrder >= lincsOrderThreshold);
1532     bool shakeWithSufficientTolerance =
1533             (ir->eConstrAlg == econtSHAKE && ir->shake_tol <= 1.1 * shakeToleranceThreshold);
1534     if (ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol <= 1.1 * bufferToleranceThreshold
1535         && (lincsWithSufficientTolerance || shakeWithSufficientTolerance))
1536     {
1537         return;
1538     }
1539
1540     bool haveDecoupledMode = false;
1541     for (const gmx_moltype_t& molt : mtop->moltype)
1542     {
1543         if (haveDecoupledModeInMol(molt, mtop->ffparams.iparams, massFactorThreshold))
1544         {
1545             haveDecoupledMode = true;
1546         }
1547     }
1548
1549     if (haveDecoupledMode)
1550     {
1551         std::string message = gmx::formatString(
1552                 "There are atoms at both ends of an angle, connected by constraints "
1553                 "and with masses that differ by more than a factor of %g. This means "
1554                 "that there are likely dynamic modes that are only very weakly coupled.",
1555                 massFactorThreshold);
1556         if (ir->cutoff_scheme == ecutsVERLET)
1557         {
1558             message += gmx::formatString(
1559                     " To ensure good equipartitioning, you need to either not use "
1560                     "constraints on all bonds (but, if possible, only on bonds involving "
1561                     "hydrogens) or use integrator = %s or decrease one or more tolerances: "
1562                     "verlet-buffer-tolerance <= %g, LINCS iterations >= %d, LINCS order "
1563                     ">= %d or SHAKE tolerance <= %g",
1564                     ei_names[eiSD1], bufferToleranceThreshold, lincsIterationThreshold,
1565                     lincsOrderThreshold, shakeToleranceThreshold);
1566         }
1567         else
1568         {
1569             message += gmx::formatString(
1570                     " To ensure good equipartitioning, we suggest to switch to the %s "
1571                     "cutoff-scheme, since that allows for better control over the Verlet "
1572                     "buffer size and thus over the energy drift.",
1573                     ecutscheme_names[ecutsVERLET]);
1574         }
1575         warning(wi, message);
1576     }
1577 }
1578
1579 static void set_verlet_buffer(const gmx_mtop_t* mtop, t_inputrec* ir, real buffer_temp, matrix box, warninp* wi)
1580 {
1581     char warn_buf[STRLEN];
1582
1583     printf("Determining Verlet buffer for a tolerance of %g kJ/mol/ps at %g K\n", ir->verletbuf_tol,
1584            buffer_temp);
1585
1586     /* Calculate the buffer size for simple atom vs atoms list */
1587     VerletbufListSetup listSetup1x1;
1588     listSetup1x1.cluster_size_i = 1;
1589     listSetup1x1.cluster_size_j = 1;
1590     const real rlist_1x1 = calcVerletBufferSize(*mtop, det(box), *ir, ir->nstlist, ir->nstlist - 1,
1591                                                 buffer_temp, listSetup1x1);
1592
1593     /* Set the pair-list buffer size in ir */
1594     VerletbufListSetup listSetup4x4 = verletbufGetSafeListSetup(ListSetupType::CpuNoSimd);
1595     ir->rlist = calcVerletBufferSize(*mtop, det(box), *ir, ir->nstlist, ir->nstlist - 1,
1596                                      buffer_temp, listSetup4x4);
1597
1598     const int n_nonlin_vsite = countNonlinearVsites(*mtop);
1599     if (n_nonlin_vsite > 0)
1600     {
1601         sprintf(warn_buf,
1602                 "There are %d non-linear virtual site constructions. Their contribution to the "
1603                 "energy error is approximated. In most cases this does not affect the error "
1604                 "significantly.",
1605                 n_nonlin_vsite);
1606         warning_note(wi, warn_buf);
1607     }
1608
1609     printf("Calculated rlist for %dx%d atom pair-list as %.3f nm, buffer size %.3f nm\n", 1, 1,
1610            rlist_1x1, rlist_1x1 - std::max(ir->rvdw, ir->rcoulomb));
1611
1612     printf("Set rlist, assuming %dx%d atom pair-list, to %.3f nm, buffer size %.3f nm\n",
1613            listSetup4x4.cluster_size_i, listSetup4x4.cluster_size_j, ir->rlist,
1614            ir->rlist - std::max(ir->rvdw, ir->rcoulomb));
1615
1616     printf("Note that mdrun will redetermine rlist based on the actual pair-list setup\n");
1617
1618     if (gmx::square(ir->rlist) >= max_cutoff2(ir->ePBC, box))
1619     {
1620         gmx_fatal(FARGS,
1621                   "The pair-list cut-off (%g nm) is longer than half the shortest box vector or "
1622                   "longer than the smallest box diagonal element (%g nm). Increase the box size or "
1623                   "decrease nstlist or increase verlet-buffer-tolerance.",
1624                   ir->rlist, std::sqrt(max_cutoff2(ir->ePBC, box)));
1625     }
1626 }
1627
1628 int gmx_grompp(int argc, char* argv[])
1629 {
1630     const char* desc[] = {
1631         "[THISMODULE] (the gromacs preprocessor)",
1632         "reads a molecular topology file, checks the validity of the",
1633         "file, expands the topology from a molecular description to an atomic",
1634         "description. The topology file contains information about",
1635         "molecule types and the number of molecules, the preprocessor",
1636         "copies each molecule as needed. ",
1637         "There is no limitation on the number of molecule types. ",
1638         "Bonds and bond-angles can be converted into constraints, separately",
1639         "for hydrogens and heavy atoms.",
1640         "Then a coordinate file is read and velocities can be generated",
1641         "from a Maxwellian distribution if requested.",
1642         "[THISMODULE] also reads parameters for [gmx-mdrun] ",
1643         "(eg. number of MD steps, time step, cut-off), and others such as",
1644         "NEMD parameters, which are corrected so that the net acceleration",
1645         "is zero.",
1646         "Eventually a binary file is produced that can serve as the sole input",
1647         "file for the MD program.[PAR]",
1648
1649         "[THISMODULE] uses the atom names from the topology file. The atom names",
1650         "in the coordinate file (option [TT]-c[tt]) are only read to generate",
1651         "warnings when they do not match the atom names in the topology.",
1652         "Note that the atom names are irrelevant for the simulation as",
1653         "only the atom types are used for generating interaction parameters.[PAR]",
1654
1655         "[THISMODULE] uses a built-in preprocessor to resolve includes, macros, ",
1656         "etc. The preprocessor supports the following keywords::",
1657         "",
1658         "    #ifdef VARIABLE",
1659         "    #ifndef VARIABLE",
1660         "    #else",
1661         "    #endif",
1662         "    #define VARIABLE",
1663         "    #undef VARIABLE",
1664         "    #include \"filename\"",
1665         "    #include <filename>",
1666         "",
1667         "The functioning of these statements in your topology may be modulated by",
1668         "using the following two flags in your [REF].mdp[ref] file::",
1669         "",
1670         "    define = -DVARIABLE1 -DVARIABLE2",
1671         "    include = -I/home/john/doe",
1672         "",
1673         "For further information a C-programming textbook may help you out.",
1674         "Specifying the [TT]-pp[tt] flag will get the pre-processed",
1675         "topology file written out so that you can verify its contents.[PAR]",
1676
1677         "When using position restraints, a file with restraint coordinates",
1678         "must be supplied with [TT]-r[tt] (can be the same file as supplied",
1679         "for [TT]-c[tt]). For free energy calculations, separate reference",
1680         "coordinates for the B topology can be supplied with [TT]-rb[tt],",
1681         "otherwise they will be equal to those of the A topology.[PAR]",
1682
1683         "Starting coordinates can be read from trajectory with [TT]-t[tt].",
1684         "The last frame with coordinates and velocities will be read,",
1685         "unless the [TT]-time[tt] option is used. Only if this information",
1686         "is absent will the coordinates in the [TT]-c[tt] file be used.",
1687         "Note that these velocities will not be used when [TT]gen_vel = yes[tt]",
1688         "in your [REF].mdp[ref] file. An energy file can be supplied with",
1689         "[TT]-e[tt] to read Nose-Hoover and/or Parrinello-Rahman coupling",
1690         "variables.[PAR]",
1691
1692         "[THISMODULE] can be used to restart simulations (preserving",
1693         "continuity) by supplying just a checkpoint file with [TT]-t[tt].",
1694         "However, for simply changing the number of run steps to extend",
1695         "a run, using [gmx-convert-tpr] is more convenient than [THISMODULE].",
1696         "You then supply the old checkpoint file directly to [gmx-mdrun]",
1697         "with [TT]-cpi[tt]. If you wish to change the ensemble or things",
1698         "like output frequency, then supplying the checkpoint file to",
1699         "[THISMODULE] with [TT]-t[tt] along with a new [REF].mdp[ref] file",
1700         "with [TT]-f[tt] is the recommended procedure. Actually preserving",
1701         "the ensemble (if possible) still requires passing the checkpoint",
1702         "file to [gmx-mdrun] [TT]-cpi[tt].[PAR]",
1703
1704         "By default, all bonded interactions which have constant energy due to",
1705         "virtual site constructions will be removed. If this constant energy is",
1706         "not zero, this will result in a shift in the total energy. All bonded",
1707         "interactions can be kept by turning off [TT]-rmvsbds[tt]. Additionally,",
1708         "all constraints for distances which will be constant anyway because",
1709         "of virtual site constructions will be removed. If any constraints remain",
1710         "which involve virtual sites, a fatal error will result.[PAR]",
1711
1712         "To verify your run input file, please take note of all warnings",
1713         "on the screen, and correct where necessary. Do also look at the contents",
1714         "of the [TT]mdout.mdp[tt] file; this contains comment lines, as well as",
1715         "the input that [THISMODULE] has read. If in doubt, you can start [THISMODULE]",
1716         "with the [TT]-debug[tt] option which will give you more information",
1717         "in a file called [TT]grompp.log[tt] (along with real debug info). You",
1718         "can see the contents of the run input file with the [gmx-dump]",
1719         "program. [gmx-check] can be used to compare the contents of two",
1720         "run input files.[PAR]",
1721
1722         "The [TT]-maxwarn[tt] option can be used to override warnings printed",
1723         "by [THISMODULE] that otherwise halt output. In some cases, warnings are",
1724         "harmless, but usually they are not. The user is advised to carefully",
1725         "interpret the output messages before attempting to bypass them with",
1726         "this option."
1727     };
1728     t_gromppopts*                        opts;
1729     std::vector<MoleculeInformation>     mi;
1730     std::unique_ptr<MoleculeInformation> intermolecular_interactions;
1731     int                                  nvsite, comb;
1732     real                                 fudgeQQ;
1733     double                               reppow;
1734     const char*                          mdparin;
1735     bool                                 bNeedVel, bGenVel;
1736     gmx_bool                             have_atomnumber;
1737     gmx_output_env_t*                    oenv;
1738     gmx_bool                             bVerbose = FALSE;
1739     warninp*                             wi;
1740     char                                 warn_buf[STRLEN];
1741
1742     t_filenm fnm[] = { { efMDP, nullptr, nullptr, ffREAD },
1743                        { efMDP, "-po", "mdout", ffWRITE },
1744                        { efSTX, "-c", nullptr, ffREAD },
1745                        { efSTX, "-r", "restraint", ffOPTRD },
1746                        { efSTX, "-rb", "restraint", ffOPTRD },
1747                        { efNDX, nullptr, nullptr, ffOPTRD },
1748                        { efTOP, nullptr, nullptr, ffREAD },
1749                        { efTOP, "-pp", "processed", ffOPTWR },
1750                        { efTPR, "-o", nullptr, ffWRITE },
1751                        { efTRN, "-t", nullptr, ffOPTRD },
1752                        { efEDR, "-e", nullptr, ffOPTRD },
1753                        /* This group is needed by the VMD viewer as the start configuration for IMD sessions: */
1754                        { efGRO, "-imd", "imdgroup", ffOPTWR },
1755                        { efTRN, "-ref", "rotref", ffOPTRW } };
1756 #define NFILE asize(fnm)
1757
1758     /* Command line options */
1759     gmx_bool bRenum   = TRUE;
1760     gmx_bool bRmVSBds = TRUE, bZero = FALSE;
1761     int      i, maxwarn             = 0;
1762     real     fr_time = -1;
1763     t_pargs  pa[]    = {
1764         { "-v", FALSE, etBOOL, { &bVerbose }, "Be loud and noisy" },
1765         { "-time", FALSE, etREAL, { &fr_time }, "Take frame at or first after this time." },
1766         { "-rmvsbds",
1767           FALSE,
1768           etBOOL,
1769           { &bRmVSBds },
1770           "Remove constant bonded interactions with virtual sites" },
1771         { "-maxwarn",
1772           FALSE,
1773           etINT,
1774           { &maxwarn },
1775           "Number of allowed warnings during input processing. Not for normal use and may "
1776           "generate unstable systems" },
1777         { "-zero",
1778           FALSE,
1779           etBOOL,
1780           { &bZero },
1781           "Set parameters for bonded interactions without defaults to zero instead of "
1782           "generating an error" },
1783         { "-renum",
1784           FALSE,
1785           etBOOL,
1786           { &bRenum },
1787           "Renumber atomtypes and minimize number of atomtypes" }
1788     };
1789
1790     /* Parse the command line */
1791     if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
1792     {
1793         return 0;
1794     }
1795
1796     /* Initiate some variables */
1797     gmx::MDModules mdModules;
1798     t_inputrec     irInstance;
1799     t_inputrec*    ir = &irInstance;
1800     snew(opts, 1);
1801     snew(opts->include, STRLEN);
1802     snew(opts->define, STRLEN);
1803
1804     wi = init_warning(TRUE, maxwarn);
1805
1806     /* PARAMETER file processing */
1807     mdparin = opt2fn("-f", NFILE, fnm);
1808     set_warning_line(wi, mdparin, -1);
1809     try
1810     {
1811         get_ir(mdparin, opt2fn("-po", NFILE, fnm), &mdModules, ir, opts, WriteMdpHeader::yes, wi);
1812     }
1813     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
1814
1815     if (bVerbose)
1816     {
1817         fprintf(stderr, "checking input for internal consistency...\n");
1818     }
1819     check_ir(mdparin, mdModules.notifier(), ir, opts, wi);
1820
1821     if (ir->ld_seed == -1)
1822     {
1823         ir->ld_seed = static_cast<int>(gmx::makeRandomSeed());
1824         fprintf(stderr, "Setting the LD random seed to %" PRId64 "\n", ir->ld_seed);
1825     }
1826
1827     if (ir->expandedvals->lmc_seed == -1)
1828     {
1829         ir->expandedvals->lmc_seed = static_cast<int>(gmx::makeRandomSeed());
1830         fprintf(stderr, "Setting the lambda MC random seed to %d\n", ir->expandedvals->lmc_seed);
1831     }
1832
1833     bNeedVel = EI_STATE_VELOCITY(ir->eI);
1834     bGenVel  = (bNeedVel && opts->bGenVel);
1835     if (bGenVel && ir->bContinuation)
1836     {
1837         sprintf(warn_buf,
1838                 "Generating velocities is inconsistent with attempting "
1839                 "to continue a previous run. Choose only one of "
1840                 "gen-vel = yes and continuation = yes.");
1841         warning_error(wi, warn_buf);
1842     }
1843
1844     std::array<InteractionsOfType, F_NRE> interactions;
1845     gmx_mtop_t                            sys;
1846     PreprocessingAtomTypes                atypes;
1847     if (debug)
1848     {
1849         pr_symtab(debug, 0, "Just opened", &sys.symtab);
1850     }
1851
1852     const char* fn = ftp2fn(efTOP, NFILE, fnm);
1853     if (!gmx_fexist(fn))
1854     {
1855         gmx_fatal(FARGS, "%s does not exist", fn);
1856     }
1857
1858     t_state state;
1859     new_status(fn, opt2fn_null("-pp", NFILE, fnm), opt2fn("-c", NFILE, fnm), opts, ir, bZero,
1860                bGenVel, bVerbose, &state, &atypes, &sys, &mi, &intermolecular_interactions,
1861                interactions, &comb, &reppow, &fudgeQQ, opts->bMorse, wi);
1862
1863     if (debug)
1864     {
1865         pr_symtab(debug, 0, "After new_status", &sys.symtab);
1866     }
1867
1868     nvsite = 0;
1869     /* set parameters for virtual site construction (not for vsiten) */
1870     for (size_t mt = 0; mt < sys.moltype.size(); mt++)
1871     {
1872         nvsite += set_vsites(bVerbose, &sys.moltype[mt].atoms, &atypes, mi[mt].interactions);
1873     }
1874     /* now throw away all obsolete bonds, angles and dihedrals: */
1875     /* note: constraints are ALWAYS removed */
1876     if (nvsite)
1877     {
1878         for (size_t mt = 0; mt < sys.moltype.size(); mt++)
1879         {
1880             clean_vsite_bondeds(mi[mt].interactions, sys.moltype[mt].atoms.nr, bRmVSBds);
1881         }
1882     }
1883
1884     if ((count_constraints(&sys, mi, wi) != 0) && (ir->eConstrAlg == econtSHAKE))
1885     {
1886         if (ir->eI == eiCG || ir->eI == eiLBFGS)
1887         {
1888             sprintf(warn_buf, "Can not do %s with %s, use %s", EI(ir->eI),
1889                     econstr_names[econtSHAKE], econstr_names[econtLINCS]);
1890             warning_error(wi, warn_buf);
1891         }
1892         if (ir->bPeriodicMols)
1893         {
1894             sprintf(warn_buf, "Can not do periodic molecules with %s, use %s",
1895                     econstr_names[econtSHAKE], econstr_names[econtLINCS]);
1896             warning_error(wi, warn_buf);
1897         }
1898     }
1899
1900     if (EI_SD(ir->eI) && ir->etc != etcNO)
1901     {
1902         warning_note(wi, "Temperature coupling is ignored with SD integrators.");
1903     }
1904
1905     /* If we are doing QM/MM, check that we got the atom numbers */
1906     have_atomnumber = TRUE;
1907     for (gmx::index i = 0; i < gmx::ssize(atypes); i++)
1908     {
1909         have_atomnumber = have_atomnumber && (atypes.atomNumberFromAtomType(i) >= 0);
1910     }
1911     if (!have_atomnumber && ir->bQMMM)
1912     {
1913         warning_error(
1914                 wi,
1915                 "\n"
1916                 "It appears as if you are trying to run a QM/MM calculation, but the force\n"
1917                 "field you are using does not contain atom numbers fields. This is an\n"
1918                 "optional field (introduced in GROMACS 3.3) for general runs, but mandatory\n"
1919                 "for QM/MM. The good news is that it is easy to add - put the atom number as\n"
1920                 "an integer just before the mass column in ffXXXnb.itp.\n"
1921                 "NB: United atoms have the same atom numbers as normal ones.\n\n");
1922     }
1923
1924     /* Check for errors in the input now, since they might cause problems
1925      * during processing further down.
1926      */
1927     check_warning_error(wi, FARGS);
1928
1929     if (nint_ftype(&sys, mi, F_POSRES) > 0 || nint_ftype(&sys, mi, F_FBPOSRES) > 0)
1930     {
1931         if (ir->epc == epcPARRINELLORAHMAN || ir->epc == epcMTTK)
1932         {
1933             sprintf(warn_buf,
1934                     "You are combining position restraints with %s pressure coupling, which can "
1935                     "lead to instabilities. If you really want to combine position restraints with "
1936                     "pressure coupling, we suggest to use %s pressure coupling instead.",
1937                     EPCOUPLTYPE(ir->epc), EPCOUPLTYPE(epcBERENDSEN));
1938             warning_note(wi, warn_buf);
1939         }
1940
1941         const char* fn = opt2fn("-r", NFILE, fnm);
1942         const char* fnB;
1943
1944         if (!gmx_fexist(fn))
1945         {
1946             gmx_fatal(FARGS,
1947                       "Cannot find position restraint file %s (option -r).\n"
1948                       "From GROMACS-2018, you need to specify the position restraint "
1949                       "coordinate files explicitly to avoid mistakes, although you can "
1950                       "still use the same file as you specify for the -c option.",
1951                       fn);
1952         }
1953
1954         if (opt2bSet("-rb", NFILE, fnm))
1955         {
1956             fnB = opt2fn("-rb", NFILE, fnm);
1957             if (!gmx_fexist(fnB))
1958             {
1959                 gmx_fatal(FARGS,
1960                           "Cannot find B-state position restraint file %s (option -rb).\n"
1961                           "From GROMACS-2018, you need to specify the position restraint "
1962                           "coordinate files explicitly to avoid mistakes, although you can "
1963                           "still use the same file as you specify for the -c option.",
1964                           fn);
1965             }
1966         }
1967         else
1968         {
1969             fnB = fn;
1970         }
1971
1972         if (bVerbose)
1973         {
1974             fprintf(stderr, "Reading position restraint coords from %s", fn);
1975             if (strcmp(fn, fnB) == 0)
1976             {
1977                 fprintf(stderr, "\n");
1978             }
1979             else
1980             {
1981                 fprintf(stderr, " and %s\n", fnB);
1982             }
1983         }
1984         gen_posres(&sys, mi, fn, fnB, ir->refcoord_scaling, ir->ePBC, ir->posres_com, ir->posres_comB, wi);
1985     }
1986
1987     /* If we are using CMAP, setup the pre-interpolation grid */
1988     if (interactions[F_CMAP].ncmap() > 0)
1989     {
1990         init_cmap_grid(&sys.ffparams.cmap_grid, interactions[F_CMAP].cmapAngles,
1991                        interactions[F_CMAP].cmakeGridSpacing);
1992         setup_cmap(interactions[F_CMAP].cmakeGridSpacing, interactions[F_CMAP].cmapAngles,
1993                    interactions[F_CMAP].cmap, &sys.ffparams.cmap_grid);
1994     }
1995
1996     set_wall_atomtype(&atypes, opts, ir, wi);
1997     if (bRenum)
1998     {
1999         atypes.renumberTypes(interactions, &sys, ir->wall_atomtype, bVerbose);
2000     }
2001
2002     if (ir->implicit_solvent)
2003     {
2004         gmx_fatal(FARGS, "Implicit solvation is no longer supported");
2005     }
2006
2007     /* PELA: Copy the atomtype data to the topology atomtype list */
2008     atypes.copyTot_atomtypes(&(sys.atomtypes));
2009
2010     if (debug)
2011     {
2012         pr_symtab(debug, 0, "After atype.renumberTypes", &sys.symtab);
2013     }
2014
2015     if (bVerbose)
2016     {
2017         fprintf(stderr, "converting bonded parameters...\n");
2018     }
2019
2020     const int ntype = atypes.size();
2021     convertInteractionsOfType(ntype, interactions, mi, intermolecular_interactions.get(), comb,
2022                               reppow, fudgeQQ, &sys);
2023
2024     if (debug)
2025     {
2026         pr_symtab(debug, 0, "After converInteractionsOfType", &sys.symtab);
2027     }
2028
2029     /* set ptype to VSite for virtual sites */
2030     for (gmx_moltype_t& moltype : sys.moltype)
2031     {
2032         set_vsites_ptype(FALSE, &moltype);
2033     }
2034     if (debug)
2035     {
2036         pr_symtab(debug, 0, "After virtual sites", &sys.symtab);
2037     }
2038     /* Check velocity for virtual sites and shells */
2039     if (bGenVel)
2040     {
2041         check_vel(&sys, state.v.rvec_array());
2042     }
2043
2044     /* check for shells and inpurecs */
2045     check_shells_inputrec(&sys, ir, wi);
2046
2047     /* check masses */
2048     check_mol(&sys, wi);
2049
2050     checkForUnboundAtoms(&sys, bVerbose, wi);
2051
2052     if (EI_DYNAMICS(ir->eI) && ir->eI != eiBD)
2053     {
2054         check_bonds_timestep(&sys, ir->delta_t, wi);
2055     }
2056
2057     checkDecoupledModeAccuracy(&sys, ir, wi);
2058
2059     if (EI_ENERGY_MINIMIZATION(ir->eI) && 0 == ir->nsteps)
2060     {
2061         warning_note(
2062                 wi,
2063                 "Zero-step energy minimization will alter the coordinates before calculating the "
2064                 "energy. If you just want the energy of a single point, try zero-step MD (with "
2065                 "unconstrained_start = yes). To do multiple single-point energy evaluations of "
2066                 "different configurations of the same topology, use mdrun -rerun.");
2067     }
2068
2069     check_warning_error(wi, FARGS);
2070
2071     if (bVerbose)
2072     {
2073         fprintf(stderr, "initialising group options...\n");
2074     }
2075     do_index(mdparin, ftp2fn_null(efNDX, NFILE, fnm), &sys, bVerbose, mdModules.notifier(), ir, wi);
2076
2077     if (ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0)
2078     {
2079         if (EI_DYNAMICS(ir->eI) && inputrec2nboundeddim(ir) == 3)
2080         {
2081             real buffer_temp;
2082
2083             if (EI_MD(ir->eI) && ir->etc == etcNO)
2084             {
2085                 if (bGenVel)
2086                 {
2087                     buffer_temp = opts->tempi;
2088                 }
2089                 else
2090                 {
2091                     buffer_temp = calc_temp(&sys, ir, state.v.rvec_array());
2092                 }
2093                 if (buffer_temp > 0)
2094                 {
2095                     sprintf(warn_buf,
2096                             "NVE simulation: will use the initial temperature of %.3f K for "
2097                             "determining the Verlet buffer size",
2098                             buffer_temp);
2099                     warning_note(wi, warn_buf);
2100                 }
2101                 else
2102                 {
2103                     sprintf(warn_buf,
2104                             "NVE simulation with an initial temperature of zero: will use a Verlet "
2105                             "buffer of %d%%. Check your energy drift!",
2106                             gmx::roundToInt(verlet_buffer_ratio_NVE_T0 * 100));
2107                     warning_note(wi, warn_buf);
2108                 }
2109             }
2110             else
2111             {
2112                 buffer_temp = get_max_reference_temp(ir, wi);
2113             }
2114
2115             if (EI_MD(ir->eI) && ir->etc == etcNO && buffer_temp == 0)
2116             {
2117                 /* NVE with initial T=0: we add a fixed ratio to rlist.
2118                  * Since we don't actually use verletbuf_tol, we set it to -1
2119                  * so it can't be misused later.
2120                  */
2121                 ir->rlist *= 1.0 + verlet_buffer_ratio_NVE_T0;
2122                 ir->verletbuf_tol = -1;
2123             }
2124             else
2125             {
2126                 /* We warn for NVE simulations with a drift tolerance that
2127                  * might result in a 1(.1)% drift over the total run-time.
2128                  * Note that we can't warn when nsteps=0, since we don't
2129                  * know how many steps the user intends to run.
2130                  */
2131                 if (EI_MD(ir->eI) && ir->etc == etcNO && ir->nstlist > 1 && ir->nsteps > 0)
2132                 {
2133                     const real driftTolerance = 0.01;
2134                     /* We use 2 DOF per atom = 2kT pot+kin energy,
2135                      * to be on the safe side with constraints.
2136                      */
2137                     const real totalEnergyDriftPerAtomPerPicosecond =
2138                             2 * BOLTZ * buffer_temp / (ir->nsteps * ir->delta_t);
2139
2140                     if (ir->verletbuf_tol > 1.1 * driftTolerance * totalEnergyDriftPerAtomPerPicosecond)
2141                     {
2142                         sprintf(warn_buf,
2143                                 "You are using a Verlet buffer tolerance of %g kJ/mol/ps for an "
2144                                 "NVE simulation of length %g ps, which can give a final drift of "
2145                                 "%d%%. For conserving energy to %d%% when using constraints, you "
2146                                 "might need to set verlet-buffer-tolerance to %.1e.",
2147                                 ir->verletbuf_tol, ir->nsteps * ir->delta_t,
2148                                 gmx::roundToInt(ir->verletbuf_tol / totalEnergyDriftPerAtomPerPicosecond * 100),
2149                                 gmx::roundToInt(100 * driftTolerance),
2150                                 driftTolerance * totalEnergyDriftPerAtomPerPicosecond);
2151                         warning_note(wi, warn_buf);
2152                     }
2153                 }
2154
2155                 set_verlet_buffer(&sys, ir, buffer_temp, state.box, wi);
2156             }
2157         }
2158     }
2159
2160     /* Init the temperature coupling state */
2161     init_gtc_state(&state, ir->opts.ngtc, 0, ir->opts.nhchainlength); /* need to add nnhpres here? */
2162
2163     if (debug)
2164     {
2165         pr_symtab(debug, 0, "After index", &sys.symtab);
2166     }
2167
2168     triple_check(mdparin, ir, &sys, wi);
2169     close_symtab(&sys.symtab);
2170     if (debug)
2171     {
2172         pr_symtab(debug, 0, "After close", &sys.symtab);
2173     }
2174
2175     /* make exclusions between QM atoms and remove charges if needed */
2176     if (ir->bQMMM)
2177     {
2178         generate_qmexcl(&sys, ir, wi, GmxQmmmMode::GMX_QMMM_ORIGINAL);
2179         if (ir->QMMMscheme != eQMMMschemeoniom)
2180         {
2181             std::vector<int> qmmmAtoms = qmmmAtomIndices(*ir, sys);
2182             removeQmmmAtomCharges(&sys, qmmmAtoms);
2183         }
2184     }
2185
2186     if (ir->eI == eiMimic)
2187     {
2188         generate_qmexcl(&sys, ir, wi, GmxQmmmMode::GMX_QMMM_MIMIC);
2189     }
2190
2191     if (ftp2bSet(efTRN, NFILE, fnm))
2192     {
2193         if (bVerbose)
2194         {
2195             fprintf(stderr, "getting data from old trajectory ...\n");
2196         }
2197         cont_status(ftp2fn(efTRN, NFILE, fnm), ftp2fn_null(efEDR, NFILE, fnm), bNeedVel, bGenVel,
2198                     fr_time, ir, &state, &sys, oenv);
2199     }
2200
2201     if (ir->ePBC == epbcXY && ir->nwall != 2)
2202     {
2203         clear_rvec(state.box[ZZ]);
2204     }
2205
2206     if (EEL_FULL(ir->coulombtype) || EVDW_PME(ir->vdwtype))
2207     {
2208         /* Calculate the optimal grid dimensions */
2209         matrix          scaledBox;
2210         EwaldBoxZScaler boxScaler(*ir);
2211         boxScaler.scaleBox(state.box, scaledBox);
2212
2213         if (ir->nkx > 0 && ir->nky > 0 && ir->nkz > 0)
2214         {
2215             /* Mark fourier_spacing as not used */
2216             ir->fourier_spacing = 0;
2217         }
2218         else if (ir->nkx != 0 && ir->nky != 0 && ir->nkz != 0)
2219         {
2220             set_warning_line(wi, mdparin, -1);
2221             warning_error(
2222                     wi, "Some of the Fourier grid sizes are set, but all of them need to be set.");
2223         }
2224         const int minGridSize = minimalPmeGridSize(ir->pme_order);
2225         calcFftGrid(stdout, scaledBox, ir->fourier_spacing, minGridSize, &(ir->nkx), &(ir->nky),
2226                     &(ir->nkz));
2227         if (ir->nkx < minGridSize || ir->nky < minGridSize || ir->nkz < minGridSize)
2228         {
2229             warning_error(wi,
2230                           "The PME grid size should be >= 2*(pme-order - 1); either manually "
2231                           "increase the grid size or decrease pme-order");
2232         }
2233     }
2234
2235     /* MRS: eventually figure out better logic for initializing the fep
2236        values that makes declaring the lambda and declaring the state not
2237        potentially conflict if not handled correctly. */
2238     if (ir->efep != efepNO)
2239     {
2240         state.fep_state = ir->fepvals->init_fep_state;
2241         for (i = 0; i < efptNR; i++)
2242         {
2243             /* init_lambda trumps state definitions*/
2244             if (ir->fepvals->init_lambda >= 0)
2245             {
2246                 state.lambda[i] = ir->fepvals->init_lambda;
2247             }
2248             else
2249             {
2250                 if (ir->fepvals->all_lambda[i] == nullptr)
2251                 {
2252                     gmx_fatal(FARGS, "Values of lambda not set for a free energy calculation!");
2253                 }
2254                 else
2255                 {
2256                     state.lambda[i] = ir->fepvals->all_lambda[i][state.fep_state];
2257                 }
2258             }
2259         }
2260     }
2261
2262     pull_t* pull = nullptr;
2263
2264     if (ir->bPull)
2265     {
2266         pull = set_pull_init(ir, &sys, state.x.rvec_array(), state.box, state.lambda[efptMASS], wi);
2267     }
2268
2269     /* Modules that supply external potential for pull coordinates
2270      * should register those potentials here. finish_pull() will check
2271      * that providers have been registerd for all external potentials.
2272      */
2273
2274     if (ir->bDoAwh)
2275     {
2276         tensor compressibility = { { 0 } };
2277         if (ir->epc != epcNO)
2278         {
2279             copy_mat(ir->compress, compressibility);
2280         }
2281         setStateDependentAwhParams(ir->awhParams, ir->pull, pull, state.box, ir->ePBC,
2282                                    compressibility, &ir->opts, wi);
2283     }
2284
2285     if (ir->bPull)
2286     {
2287         finish_pull(pull);
2288     }
2289
2290     if (ir->bRot)
2291     {
2292         set_reference_positions(ir->rot, state.x.rvec_array(), state.box,
2293                                 opt2fn("-ref", NFILE, fnm), opt2bSet("-ref", NFILE, fnm), wi);
2294     }
2295
2296     /*  reset_multinr(sys); */
2297
2298     if (EEL_PME(ir->coulombtype))
2299     {
2300         float ratio = pme_load_estimate(sys, *ir, state.box);
2301         fprintf(stderr, "Estimate for the relative computational load of the PME mesh part: %.2f\n", ratio);
2302         /* With free energy we might need to do PME both for the A and B state
2303          * charges. This will double the cost, but the optimal performance will
2304          * then probably be at a slightly larger cut-off and grid spacing.
2305          */
2306         if ((ir->efep == efepNO && ratio > 1.0 / 2.0) || (ir->efep != efepNO && ratio > 2.0 / 3.0))
2307         {
2308             warning_note(
2309                     wi,
2310                     "The optimal PME mesh load for parallel simulations is below 0.5\n"
2311                     "and for highly parallel simulations between 0.25 and 0.33,\n"
2312                     "for higher performance, increase the cut-off and the PME grid spacing.\n");
2313             if (ir->efep != efepNO)
2314             {
2315                 warning_note(wi,
2316                              "For free energy simulations, the optimal load limit increases from "
2317                              "0.5 to 0.667\n");
2318             }
2319         }
2320     }
2321
2322     {
2323         char   warn_buf[STRLEN];
2324         double cio = compute_io(ir, sys.natoms, sys.groups, F_NRE, 1);
2325         sprintf(warn_buf, "This run will generate roughly %.0f Mb of data", cio);
2326         if (cio > 2000)
2327         {
2328             set_warning_line(wi, mdparin, -1);
2329             warning_note(wi, warn_buf);
2330         }
2331         else
2332         {
2333             printf("%s\n", warn_buf);
2334         }
2335     }
2336
2337     // Add the md modules internal parameters that are not mdp options
2338     // e.g., atom indices
2339
2340     {
2341         gmx::KeyValueTreeBuilder internalParameterBuilder;
2342         mdModules.notifier().notifier_.notify(internalParameterBuilder.rootObject());
2343         ir->internalParameters =
2344                 std::make_unique<gmx::KeyValueTreeObject>(internalParameterBuilder.build());
2345     }
2346
2347     if (bVerbose)
2348     {
2349         fprintf(stderr, "writing run input file...\n");
2350     }
2351
2352     done_warning(wi, FARGS);
2353     write_tpx_state(ftp2fn(efTPR, NFILE, fnm), ir, &state, &sys);
2354
2355     /* Output IMD group, if bIMD is TRUE */
2356     gmx::write_IMDgroup_to_file(ir->bIMD, ir, &state, &sys, NFILE, fnm);
2357
2358     sfree(opts->define);
2359     sfree(opts->include);
2360     sfree(opts);
2361     for (auto& mol : mi)
2362     {
2363         // Some of the contents of molinfo have been stolen, so
2364         // fullCleanUp can't be called.
2365         mol.partialCleanUp();
2366     }
2367     done_inputrec_strings();
2368     output_env_done(oenv);
2369
2370     return 0;
2371 }