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