src/gromacs/gmxpreprocess/grompp.cpp

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