src/gromacs/gmxpreprocess/grompp.cpp

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