src/gromacs/gmxpreprocess/grompp.cpp

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