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