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 InteractionType::InteractionType(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 &InteractionType::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 &InteractionType::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 &InteractionType::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 &InteractionType::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 &InteractionType::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 &InteractionType::c0() const
 158 {
 159     return forceParam_[0];
 160 }
 161
 162 const real &InteractionType::c1() const
 163 {
 164     return forceParam_[1];
 165 }
 166
 167 const real &InteractionType::c2() const
 168 {
 169     return forceParam_[2];
 170 }
 171
 172 const std::string &InteractionType::interactionTypeName() const
 173 {
 174     return interactionTypeName_;
 175 }
 176
 177 void InteractionType::sortBondAtomIds()
 178 {
 179     if (aj() < ai())
 180     {
 181         std::swap(atoms_[0], atoms_[1]);
 182     }
 183 }
 184
 185 void InteractionType::sortAngleAtomIds()
 186 {
 187     if (ak() < ai())
 188     {
 189         std::swap(atoms_[0], atoms_[2]);
 190     }
 191 }
 192
 193 void InteractionType::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 InteractionType::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 InteractionType::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.plist[ifunc].size();
 255         mol.plist[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].plist[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<InteractionTypeParameters> plist,
 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                        plist, 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     /* If using the group scheme, make sure charge groups are made whole to avoid errors
 772      * in calculating charge group size later on
 773      */
 774     if (ir->cutoff_scheme == ecutsGROUP && ir->ePBC != epbcNONE)
 775     {
 776         // Need temporary rvec for coordinates
 777         do_pbc_first_mtop(nullptr, ir->ePBC, state->box, sys, state->x.rvec_array());
 778     }
 779
 780     /* Do more checks, mostly related to constraints */
 781     if (bVerbose)
 782     {
 783         fprintf(stderr, "double-checking input for internal consistency...\n");
 784     }
 785     {
 786         bool bHasNormalConstraints = 0 < (nint_ftype(sys, *mi, F_CONSTR) +
 787                                           nint_ftype(sys, *mi, F_CONSTRNC));
 788         bool bHasAnyConstraints = bHasNormalConstraints || 0 < nint_ftype(sys, *mi, F_SETTLE);
 789         double_check(ir, state->box,
 790                      bHasNormalConstraints,
 791                      bHasAnyConstraints,
 792                      wi);
 793     }
 794
 795     if (bGenVel)
 796     {
 797         real                   *mass;
 798
 799         snew(mass, state->natoms);
 800         for (const AtomProxy atomP : AtomRange(*sys))
 801         {
 802             const t_atom &local = atomP.atom();
 803             int           i     = atomP.globalAtomNumber();
 804             mass[i] = local.m;
 805         }
 806
 807         if (opts->seed == -1)
 808         {
 809             opts->seed = static_cast<int>(gmx::makeRandomSeed());
 810             fprintf(stderr, "Setting gen_seed to %d\n", opts->seed);
 811         }
 812         state->flags |= (1 << estV);
 813         maxwell_speed(opts->tempi, opts->seed, sys, state->v.rvec_array());
 814
 815         stop_cm(stdout, state->natoms, mass, state->x.rvec_array(), state->v.rvec_array());
 816         sfree(mass);
 817     }
 818 }
 819
 820 static void copy_state(const char *slog, t_trxframe *fr,
 821                        bool bReadVel, t_state *state,
 822                        double *use_time)
 823 {
 824     if (fr->not_ok & FRAME_NOT_OK)
 825     {
 826         gmx_fatal(FARGS, "Can not start from an incomplete frame");
 827     }
 828     if (!fr->bX)
 829     {
 830         gmx_fatal(FARGS, "Did not find a frame with coordinates in file %s",
 831                   slog);
 832     }
 833
 834     std::copy(fr->x, fr->x+state->natoms, state->x.data());
 835     if (bReadVel)
 836     {
 837         if (!fr->bV)
 838         {
 839             gmx_incons("Trajecory frame unexpectedly does not contain velocities");
 840         }
 841         std::copy(fr->v, fr->v+state->natoms, state->v.data());
 842     }
 843     if (fr->bBox)
 844     {
 845         copy_mat(fr->box, state->box);
 846     }
 847
 848     *use_time = fr->time;
 849 }
 850
 851 static void cont_status(const char *slog, const char *ener,
 852                         bool bNeedVel, bool bGenVel, real fr_time,
 853                         t_inputrec *ir, t_state *state,
 854                         gmx_mtop_t *sys,
 855                         const gmx_output_env_t *oenv)
 856 /* If fr_time == -1 read the last frame available which is complete */
 857 {
 858     bool         bReadVel;
 859     t_trxframe   fr;
 860     t_trxstatus *fp;
 861     double       use_time;
 862
 863     bReadVel = (bNeedVel && !bGenVel);
 864
 865     fprintf(stderr,
 866             "Reading Coordinates%s and Box size from old trajectory\n",
 867             bReadVel ? ", Velocities" : "");
 868     if (fr_time == -1)
 869     {
 870         fprintf(stderr, "Will read whole trajectory\n");
 871     }
 872     else
 873     {
 874         fprintf(stderr, "Will read till time %g\n", fr_time);
 875     }
 876     if (!bReadVel)
 877     {
 878         if (bGenVel)
 879         {
 880             fprintf(stderr, "Velocities generated: "
 881                     "ignoring velocities in input trajectory\n");
 882         }
 883         read_first_frame(oenv, &fp, slog, &fr, TRX_NEED_X);
 884     }
 885     else
 886     {
 887         read_first_frame(oenv, &fp, slog, &fr, TRX_NEED_X | TRX_NEED_V);
 888
 889         if (!fr.bV)
 890         {
 891             fprintf(stderr,
 892                     "\n"
 893                     "WARNING: Did not find a frame with velocities in file %s,\n"
 894                     "         all velocities will be set to zero!\n\n", slog);
 895             for (auto &vi : makeArrayRef(state->v))
 896             {
 897                 vi = {0, 0, 0};
 898             }
 899             close_trx(fp);
 900             /* Search for a frame without velocities */
 901             bReadVel = false;
 902             read_first_frame(oenv, &fp, slog, &fr, TRX_NEED_X);
 903         }
 904     }
 905
 906     state->natoms = fr.natoms;
 907
 908     if (sys->natoms != state->natoms)
 909     {
 910         gmx_fatal(FARGS, "Number of atoms in Topology "
 911                   "is not the same as in Trajectory");
 912     }
 913     copy_state(slog, &fr, bReadVel, state, &use_time);
 914
 915     /* Find the appropriate frame */
 916     while ((fr_time == -1 || fr.time < fr_time) &&
 917            read_next_frame(oenv, fp, &fr))
 918     {
 919         copy_state(slog, &fr, bReadVel, state, &use_time);
 920     }
 921
 922     close_trx(fp);
 923
 924     /* Set the relative box lengths for preserving the box shape.
 925      * Note that this call can lead to differences in the last bit
 926      * with respect to using gmx convert-tpr to create a [REF].tpx[ref] file.
 927      */
 928     set_box_rel(ir, state);
 929
 930     fprintf(stderr, "Using frame at t = %g ps\n", use_time);
 931     fprintf(stderr, "Starting time for run is %g ps\n", ir->init_t);
 932
 933     if ((ir->epc != epcNO  || ir->etc == etcNOSEHOOVER) && ener)
 934     {
 935         get_enx_state(ener, use_time, sys->groups, ir, state);
 936         preserve_box_shape(ir, state->box_rel, state->boxv);
 937     }
 938 }
 939
 940 static void read_posres(gmx_mtop_t *mtop,
 941                         gmx::ArrayRef<const MoleculeInformation> molinfo,
 942                         gmx_bool bTopB,
 943                         const char *fn,
 944                         int rc_scaling, int ePBC,
 945                         rvec com,
 946                         warninp *wi)
 947 {
 948     gmx_bool           *hadAtom;
 949     rvec               *x, *v;
 950     dvec                sum;
 951     double              totmass;
 952     t_topology         *top;
 953     matrix              box, invbox;
 954     int                 natoms, npbcdim = 0;
 955     char                warn_buf[STRLEN];
 956     int                 a, nat_molb;
 957     t_atom             *atom;
 958
 959     snew(top, 1);
 960     read_tps_conf(fn, top, nullptr, &x, &v, box, FALSE);
 961     natoms = top->atoms.nr;
 962     done_top(top);
 963     sfree(top);
 964     if (natoms != mtop->natoms)
 965     {
 966         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);
 967         warning(wi, warn_buf);
 968     }
 969
 970     npbcdim = ePBC2npbcdim(ePBC);
 971     clear_rvec(com);
 972     if (rc_scaling != erscNO)
 973     {
 974         copy_mat(box, invbox);
 975         for (int j = npbcdim; j < DIM; j++)
 976         {
 977             clear_rvec(invbox[j]);
 978             invbox[j][j] = 1;
 979         }
 980         gmx::invertBoxMatrix(invbox, invbox);
 981     }
 982
 983     /* Copy the reference coordinates to mtop */
 984     clear_dvec(sum);
 985     totmass = 0;
 986     a       = 0;
 987     snew(hadAtom, natoms);
 988     for (gmx_molblock_t &molb : mtop->molblock)
 989     {
 990         nat_molb = molb.nmol*mtop->moltype[molb.type].atoms.nr;
 991         const InteractionTypeParameters *pr   = &(molinfo[molb.type].plist[F_POSRES]);
 992         const InteractionTypeParameters *prfb = &(molinfo[molb.type].plist[F_FBPOSRES]);
 993         if (pr->size() > 0 || prfb->size() > 0)
 994         {
 995             atom = mtop->moltype[molb.type].atoms.atom;
 996             for (const auto &restraint : pr->interactionTypes)
 997             {
 998                 int ai = restraint.ai();
 999                 if (ai >= natoms)
1000                 {
1001                     gmx_fatal(FARGS, "Position restraint atom index (%d) in moltype '%s' is larger than number of atoms in %s (%d).\n",
1002                               ai+1, *molinfo[molb.type].name, fn, natoms);
1003                 }
1004                 hadAtom[ai] = TRUE;
1005                 if (rc_scaling == erscCOM)
1006                 {
1007                     /* Determine the center of mass of the posres reference coordinates */
1008                     for (int j = 0; j < npbcdim; j++)
1009                     {
1010                         sum[j] += atom[ai].m*x[a+ai][j];
1011                     }
1012                     totmass  += atom[ai].m;
1013                 }
1014             }
1015             /* Same for flat-bottomed posres, but do not count an atom twice for COM */
1016             for (const auto &restraint : prfb->interactionTypes)
1017             {
1018                 int ai = restraint.ai();
1019                 if (ai >= natoms)
1020                 {
1021                     gmx_fatal(FARGS, "Position restraint atom index (%d) in moltype '%s' is larger than number of atoms in %s (%d).\n",
1022                               ai+1, *molinfo[molb.type].name, fn, natoms);
1023                 }
1024                 if (rc_scaling == erscCOM && !hadAtom[ai])
1025                 {
1026                     /* Determine the center of mass of the posres reference coordinates */
1027                     for (int j = 0; j < npbcdim; j++)
1028                     {
1029                         sum[j] += atom[ai].m*x[a+ai][j];
1030                     }
1031                     totmass  += atom[ai].m;
1032                 }
1033             }
1034             if (!bTopB)
1035             {
1036                 molb.posres_xA.resize(nat_molb);
1037                 for (int i = 0; i < nat_molb; i++)
1038                 {
1039                     copy_rvec(x[a+i], molb.posres_xA[i]);
1040                 }
1041             }
1042             else
1043             {
1044                 molb.posres_xB.resize(nat_molb);
1045                 for (int i = 0; i < nat_molb; i++)
1046                 {
1047                     copy_rvec(x[a+i], molb.posres_xB[i]);
1048                 }
1049             }
1050         }
1051         a += nat_molb;
1052     }
1053     if (rc_scaling == erscCOM)
1054     {
1055         if (totmass == 0)
1056         {
1057             gmx_fatal(FARGS, "The total mass of the position restraint atoms is 0");
1058         }
1059         for (int j = 0; j < npbcdim; j++)
1060         {
1061             com[j] = sum[j]/totmass;
1062         }
1063         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]);
1064     }
1065
1066     if (rc_scaling != erscNO)
1067     {
1068         GMX_ASSERT(npbcdim <= DIM, "Only DIM dimensions can have PBC");
1069
1070         for (gmx_molblock_t &molb : mtop->molblock)
1071         {
1072             nat_molb = molb.nmol*mtop->moltype[molb.type].atoms.nr;
1073             if (!molb.posres_xA.empty() || !molb.posres_xB.empty())
1074             {
1075                 std::vector<gmx::RVec> &xp = (!bTopB ? molb.posres_xA : molb.posres_xB);
1076                 for (int i = 0; i < nat_molb; i++)
1077                 {
1078                     for (int j = 0; j < npbcdim; j++)
1079                     {
1080                         if (rc_scaling == erscALL)
1081                         {
1082                             /* Convert from Cartesian to crystal coordinates */
1083                             xp[i][j] *= invbox[j][j];
1084                             for (int k = j+1; k < npbcdim; k++)
1085                             {
1086                                 xp[i][j] += invbox[k][j]*xp[i][k];
1087                             }
1088                         }
1089                         else if (rc_scaling == erscCOM)
1090                         {
1091                             /* Subtract the center of mass */
1092                             xp[i][j] -= com[j];
1093                         }
1094                     }
1095                 }
1096             }
1097         }
1098
1099         if (rc_scaling == erscCOM)
1100         {
1101             /* Convert the COM from Cartesian to crystal coordinates */
1102             for (int j = 0; j < npbcdim; j++)
1103             {
1104                 com[j] *= invbox[j][j];
1105                 for (int k = j+1; k < npbcdim; k++)
1106                 {
1107                     com[j] += invbox[k][j]*com[k];
1108                 }
1109             }
1110         }
1111     }
1112
1113     sfree(x);
1114     sfree(v);
1115     sfree(hadAtom);
1116 }
1117
1118 static void gen_posres(gmx_mtop_t *mtop,
1119                        gmx::ArrayRef<const MoleculeInformation> mi,
1120                        const char *fnA, const char *fnB,
1121                        int rc_scaling, int ePBC,
1122                        rvec com, rvec comB,
1123                        warninp *wi)
1124 {
1125     read_posres(mtop, mi, FALSE, fnA, rc_scaling, ePBC, com, wi);
1126     /* It is safer to simply read the b-state posres rather than trying
1127      * to be smart and copy the positions.
1128      */
1129     read_posres(mtop, mi, TRUE, fnB, rc_scaling, ePBC, comB, wi);
1130 }
1131
1132 static void set_wall_atomtype(PreprocessingAtomTypes *at, t_gromppopts *opts,
1133                               t_inputrec *ir, warninp *wi)
1134 {
1135     int  i;
1136     char warn_buf[STRLEN];
1137
1138     if (ir->nwall > 0)
1139     {
1140         fprintf(stderr, "Searching the wall atom type(s)\n");
1141     }
1142     for (i = 0; i < ir->nwall; i++)
1143     {
1144         ir->wall_atomtype[i] = at->atomTypeFromName(opts->wall_atomtype[i]);
1145         if (ir->wall_atomtype[i] == NOTSET)
1146         {
1147             sprintf(warn_buf, "Specified wall atom type %s is not defined", opts->wall_atomtype[i]);
1148             warning_error(wi, warn_buf);
1149         }
1150     }
1151 }
1152
1153 static int nrdf_internal(const t_atoms *atoms)
1154 {
1155     int i, nmass, nrdf;
1156
1157     nmass = 0;
1158     for (i = 0; i < atoms->nr; i++)
1159     {
1160         /* Vsite ptype might not be set here yet, so also check the mass */
1161         if ((atoms->atom[i].ptype == eptAtom ||
1162              atoms->atom[i].ptype == eptNucleus)
1163             && atoms->atom[i].m > 0)
1164         {
1165             nmass++;
1166         }
1167     }
1168     switch (nmass)
1169     {
1170         case 0:  nrdf = 0; break;
1171         case 1:  nrdf = 0; break;
1172         case 2:  nrdf = 1; break;
1173         default: nrdf = nmass*3 - 6; break;
1174     }
1175
1176     return nrdf;
1177 }
1178
1179 static void
1180 spline1d( double              dx,
1181           const double *      y,
1182           int                 n,
1183           double       *      u,
1184           double       *      y2 )
1185 {
1186     int    i;
1187     double p, q;
1188
1189     y2[0] = 0.0;
1190     u[0]  = 0.0;
1191
1192     for (i = 1; i < n-1; i++)
1193     {
1194         p     = 0.5*y2[i-1]+2.0;
1195         y2[i] = -0.5/p;
1196         q     = (y[i+1]-2.0*y[i]+y[i-1])/dx;
1197         u[i]  = (3.0*q/dx-0.5*u[i-1])/p;
1198     }
1199
1200     y2[n-1] = 0.0;
1201
1202     for (i = n-2; i >= 0; i--)
1203     {
1204         y2[i] = y2[i]*y2[i+1]+u[i];
1205     }
1206 }
1207
1208
1209 static void
1210 interpolate1d( double           xmin,
1211                double           dx,
1212                const double *   ya,
1213                const double *   y2a,
1214                double           x,
1215                double       *   y,
1216                double       *   y1)
1217 {
1218     int    ix;
1219     double a, b;
1220
1221     ix = static_cast<int>((x-xmin)/dx);
1222
1223     a = (xmin+(ix+1)*dx-x)/dx;
1224     b = (x-xmin-ix*dx)/dx;
1225
1226     *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;
1227     *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];
1228 }
1229
1230
1231 static void
1232 setup_cmap (int                       grid_spacing,
1233             int                       nc,
1234             gmx::ArrayRef<const real> grid,
1235             gmx_cmap_t          *     cmap_grid)
1236 {
1237     int                 i, j, k, ii, jj, kk, idx;
1238     int                 offset;
1239     double              dx, xmin, v, v1, v2, v12;
1240     double              phi, psi;
1241
1242     std::vector<double> tmp_u(2*grid_spacing, 0.0);
1243     std::vector<double> tmp_u2(2*grid_spacing, 0.0);
1244     std::vector<double> tmp_yy(2*grid_spacing, 0.0);
1245     std::vector<double> tmp_y1(2*grid_spacing, 0.0);
1246     std::vector<double> tmp_t2(2*grid_spacing*2*grid_spacing, 0.0);
1247     std::vector<double> tmp_grid(2*grid_spacing*2*grid_spacing, 0.0);
1248
1249     dx   = 360.0/grid_spacing;
1250     xmin = -180.0-dx*grid_spacing/2;
1251
1252     for (kk = 0; kk < nc; kk++)
1253     {
1254         /* Compute an offset depending on which cmap we are using
1255          * Offset will be the map number multiplied with the
1256          * grid_spacing * grid_spacing * 2
1257          */
1258         offset = kk * grid_spacing * grid_spacing * 2;
1259
1260         for (i = 0; i < 2*grid_spacing; i++)
1261         {
1262             ii = (i+grid_spacing-grid_spacing/2)%grid_spacing;
1263
1264             for (j = 0; j < 2*grid_spacing; j++)
1265             {
1266                 jj = (j+grid_spacing-grid_spacing/2)%grid_spacing;
1267                 tmp_grid[i*grid_spacing*2+j] = grid[offset+ii*grid_spacing+jj];
1268             }
1269         }
1270
1271         for (i = 0; i < 2*grid_spacing; i++)
1272         {
1273             spline1d(dx, &(tmp_grid[2*grid_spacing*i]), 2*grid_spacing, tmp_u.data(), &(tmp_t2[2*grid_spacing*i]));
1274         }
1275
1276         for (i = grid_spacing/2; i < grid_spacing+grid_spacing/2; i++)
1277         {
1278             ii  = i-grid_spacing/2;
1279             phi = ii*dx-180.0;
1280
1281             for (j = grid_spacing/2; j < grid_spacing+grid_spacing/2; j++)
1282             {
1283                 jj  = j-grid_spacing/2;
1284                 psi = jj*dx-180.0;
1285
1286                 for (k = 0; k < 2*grid_spacing; k++)
1287                 {
1288                     interpolate1d(xmin, dx, &(tmp_grid[2*grid_spacing*k]),
1289                                   &(tmp_t2[2*grid_spacing*k]), psi, &tmp_yy[k], &tmp_y1[k]);
1290                 }
1291
1292                 spline1d(dx, tmp_yy.data(), 2*grid_spacing, tmp_u.data(), tmp_u2.data());
1293                 interpolate1d(xmin, dx, tmp_yy.data(), tmp_u2.data(), phi, &v, &v1);
1294                 spline1d(dx, tmp_y1.data(), 2*grid_spacing, tmp_u.data(), tmp_u2.data());
1295                 interpolate1d(xmin, dx, tmp_y1.data(), tmp_u2.data(), phi, &v2, &v12);
1296
1297                 idx = ii*grid_spacing+jj;
1298                 cmap_grid->cmapdata[kk].cmap[idx*4]   = grid[offset+ii*grid_spacing+jj];
1299                 cmap_grid->cmapdata[kk].cmap[idx*4+1] = v1;
1300                 cmap_grid->cmapdata[kk].cmap[idx*4+2] = v2;
1301                 cmap_grid->cmapdata[kk].cmap[idx*4+3] = v12;
1302             }
1303         }
1304     }
1305 }
1306
1307 static void init_cmap_grid(gmx_cmap_t *cmap_grid, int ngrid, int grid_spacing)
1308 {
1309     int i, nelem;
1310
1311     cmap_grid->grid_spacing = grid_spacing;
1312     nelem                   = cmap_grid->grid_spacing*cmap_grid->grid_spacing;
1313
1314     cmap_grid->cmapdata.resize(ngrid);
1315
1316     for (i = 0; i < ngrid; i++)
1317     {
1318         cmap_grid->cmapdata[i].cmap.resize(4*nelem);
1319     }
1320 }
1321
1322
1323 static int count_constraints(const gmx_mtop_t                        *mtop,
1324                              gmx::ArrayRef<const MoleculeInformation> mi,
1325                              warninp                                 *wi)
1326 {
1327     int             count, count_mol;
1328     char            buf[STRLEN];
1329
1330     count = 0;
1331     for (const gmx_molblock_t &molb : mtop->molblock)
1332     {
1333         count_mol = 0;
1334         gmx::ArrayRef<const InteractionTypeParameters> plist = mi[molb.type].plist;
1335
1336         for (int i = 0; i < F_NRE; i++)
1337         {
1338             if (i == F_SETTLE)
1339             {
1340                 count_mol += 3*plist[i].size();
1341             }
1342             else if (interaction_function[i].flags & IF_CONSTRAINT)
1343             {
1344                 count_mol += plist[i].size();
1345             }
1346         }
1347
1348         if (count_mol > nrdf_internal(&mi[molb.type].atoms))
1349         {
1350             sprintf(buf,
1351                     "Molecule type '%s' has %d constraints.\n"
1352                     "For stability and efficiency there should not be more constraints than internal number of degrees of freedom: %d.\n",
1353                     *mi[molb.type].name, count_mol,
1354                     nrdf_internal(&mi[molb.type].atoms));
1355             warning(wi, buf);
1356         }
1357         count += molb.nmol*count_mol;
1358     }
1359
1360     return count;
1361 }
1362
1363 static real calc_temp(const gmx_mtop_t *mtop,
1364                       const t_inputrec *ir,
1365                       rvec             *v)
1366 {
1367     double                     sum_mv2 = 0;
1368     for (const AtomProxy atomP : AtomRange(*mtop))
1369     {
1370         const t_atom &local = atomP.atom();
1371         int           i     = atomP.globalAtomNumber();
1372         sum_mv2 += local.m*norm2(v[i]);
1373     }
1374
1375     double nrdf = 0;
1376     for (int g = 0; g < ir->opts.ngtc; g++)
1377     {
1378         nrdf += ir->opts.nrdf[g];
1379     }
1380
1381     return sum_mv2/(nrdf*BOLTZ);
1382 }
1383
1384 static real get_max_reference_temp(const t_inputrec *ir,
1385                                    warninp          *wi)
1386 {
1387     real         ref_t;
1388     int          i;
1389     bool         bNoCoupl;
1390
1391     ref_t    = 0;
1392     bNoCoupl = false;
1393     for (i = 0; i < ir->opts.ngtc; i++)
1394     {
1395         if (ir->opts.tau_t[i] < 0)
1396         {
1397             bNoCoupl = true;
1398         }
1399         else
1400         {
1401             ref_t = std::max(ref_t, ir->opts.ref_t[i]);
1402         }
1403     }
1404
1405     if (bNoCoupl)
1406     {
1407         char buf[STRLEN];
1408
1409         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.",
1410                 ref_t);
1411         warning(wi, buf);
1412     }
1413
1414     return ref_t;
1415 }
1416
1417 /* Checks if there are unbound atoms in moleculetype molt.
1418  * Prints a note for each unbound atoms and a warning if any is present.
1419  */
1420 static void checkForUnboundAtoms(const gmx_moltype_t     *molt,
1421                                  gmx_bool                 bVerbose,
1422                                  warninp                 *wi)
1423 {
1424     const t_atoms *atoms = &molt->atoms;
1425
1426     if (atoms->nr == 1)
1427     {
1428         /* Only one atom, there can't be unbound atoms */
1429         return;
1430     }
1431
1432     std::vector<int> count(atoms->nr, 0);
1433
1434     for (int ftype = 0; ftype < F_NRE; ftype++)
1435     {
1436         if (((interaction_function[ftype].flags & IF_BOND) && ftype != F_CONNBONDS) ||
1437             (interaction_function[ftype].flags & IF_CONSTRAINT) ||
1438             ftype == F_SETTLE)
1439         {
1440             const InteractionList &il   = molt->ilist[ftype];
1441             const int              nral = NRAL(ftype);
1442
1443             for (int i = 0; i < il.size(); i += 1 + nral)
1444             {
1445                 for (int j = 0; j < nral; j++)
1446                 {
1447                     count[il.iatoms[i + 1 + j]]++;
1448                 }
1449             }
1450         }
1451     }
1452
1453     int numDanglingAtoms = 0;
1454     for (int a = 0; a < atoms->nr; a++)
1455     {
1456         if (atoms->atom[a].ptype != eptVSite &&
1457             count[a] == 0)
1458         {
1459             if (bVerbose)
1460             {
1461                 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",
1462                         a + 1, *atoms->atomname[a], *molt->name);
1463             }
1464             numDanglingAtoms++;
1465         }
1466     }
1467
1468     if (numDanglingAtoms > 0)
1469     {
1470         char buf[STRLEN];
1471         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.",
1472                 *molt->name, numDanglingAtoms);
1473         warning_note(wi, buf);
1474     }
1475 }
1476
1477 /* Checks all moleculetypes for unbound atoms */
1478 static void checkForUnboundAtoms(const gmx_mtop_t     *mtop,
1479                                  gmx_bool              bVerbose,
1480                                  warninp              *wi)
1481 {
1482     for (const gmx_moltype_t &molt : mtop->moltype)
1483     {
1484         checkForUnboundAtoms(&molt, bVerbose, wi);
1485     }
1486 }
1487
1488 /*! \brief Checks if there are decoupled modes in moleculetype \p molt.
1489  *
1490  * The specific decoupled modes this routine check for are angle modes
1491  * where the two bonds are constrained and the atoms a both ends are only
1492  * involved in a single constraint; the mass of the two atoms needs to
1493  * differ by more than \p massFactorThreshold.
1494  */
1495 static bool haveDecoupledModeInMol(const gmx_moltype_t            &molt,
1496                                    gmx::ArrayRef<const t_iparams>  iparams,
1497                                    real                            massFactorThreshold)
1498 {
1499     if (molt.ilist[F_CONSTR].size() == 0 &&
1500         molt.ilist[F_CONSTRNC].size() == 0)
1501     {
1502         return false;
1503     }
1504
1505     const t_atom * atom = molt.atoms.atom;
1506
1507     t_blocka       atomToConstraints =
1508         gmx::make_at2con(molt, iparams,
1509                          gmx::FlexibleConstraintTreatment::Exclude);
1510
1511     bool           haveDecoupledMode = false;
1512     for (int ftype = 0; ftype < F_NRE; ftype++)
1513     {
1514         if (interaction_function[ftype].flags & IF_ATYPE)
1515         {
1516             const int              nral = NRAL(ftype);
1517             const InteractionList &il   = molt.ilist[ftype];
1518             for (int i = 0; i < il.size(); i += 1 + nral)
1519             {
1520                 /* Here we check for the mass difference between the atoms
1521                  * at both ends of the angle, that the atoms at the ends have
1522                  * 1 contraint and the atom in the middle at least 3; we check
1523                  * that the 3 atoms are linked by constraints below.
1524                  * We check for at least three constraints for the middle atom,
1525                  * since with only the two bonds in the angle, we have 3-atom
1526                  * molecule, which has much more thermal exhange in this single
1527                  * angle mode than molecules with more atoms.
1528                  * Note that this check also catches molecules where more atoms
1529                  * are connected to one or more atoms in the angle, but by
1530                  * bond potentials instead of angles. But such cases will not
1531                  * occur in "normal" molecules and it doesn't hurt running
1532                  * those with higher accuracy settings as well.
1533                  */
1534                 int a0 = il.iatoms[1 + i];
1535                 int a1 = il.iatoms[1 + i + 1];
1536                 int a2 = il.iatoms[1 + i + 2];
1537                 if ((atom[a0].m > atom[a2].m*massFactorThreshold ||
1538                      atom[a2].m > atom[a0].m*massFactorThreshold) &&
1539                     atomToConstraints.index[a0 + 1] - atomToConstraints.index[a0] == 1 &&
1540                     atomToConstraints.index[a2 + 1] - atomToConstraints.index[a2] == 1 &&
1541                     atomToConstraints.index[a1 + 1] - atomToConstraints.index[a1] >= 3)
1542                 {
1543                     int      constraint0 = atomToConstraints.a[atomToConstraints.index[a0]];
1544                     int      constraint2 = atomToConstraints.a[atomToConstraints.index[a2]];
1545
1546                     bool     foundAtom0  = false;
1547                     bool     foundAtom2  = false;
1548                     for (int conIndex = atomToConstraints.index[a1]; conIndex < atomToConstraints.index[a1 + 1]; conIndex++)
1549                     {
1550                         if (atomToConstraints.a[conIndex] == constraint0)
1551                         {
1552                             foundAtom0 = true;
1553                         }
1554                         if (atomToConstraints.a[conIndex] == constraint2)
1555                         {
1556                             foundAtom2 = true;
1557                         }
1558                     }
1559                     if (foundAtom0 && foundAtom2)
1560                     {
1561                         haveDecoupledMode = true;
1562                     }
1563                 }
1564             }
1565         }
1566     }
1567
1568     done_blocka(&atomToConstraints);
1569
1570     return haveDecoupledMode;
1571 }
1572
1573 /*! \brief Checks if the Verlet buffer and constraint accuracy is sufficient for decoupled dynamic modes.
1574  *
1575  * When decoupled modes are present and the accuray in insufficient,
1576  * this routine issues a warning if the accuracy is insufficient.
1577  */
1578 static void checkDecoupledModeAccuracy(const gmx_mtop_t *mtop,
1579                                        const t_inputrec *ir,
1580                                        warninp          *wi)
1581 {
1582     /* We only have issues with decoupled modes with normal MD.
1583      * With stochastic dynamics equipartitioning is enforced strongly.
1584      */
1585     if (!EI_MD(ir->eI))
1586     {
1587         return;
1588     }
1589
1590     /* When atoms of very different mass are involved in an angle potential
1591      * and both bonds in the angle are constrained, the dynamic modes in such
1592      * angles have very different periods and significant energy exchange
1593      * takes several nanoseconds. Thus even a small amount of error in
1594      * different algorithms can lead to violation of equipartitioning.
1595      * The parameters below are mainly based on an all-atom chloroform model
1596      * with all bonds constrained. Equipartitioning is off by more than 1%
1597      * (need to run 5-10 ns) when the difference in mass between H and Cl
1598      * is more than a factor 13 and the accuracy is less than the thresholds
1599      * given below. This has been verified on some other molecules.
1600      *
1601      * Note that the buffer and shake parameters have unit length and
1602      * energy/time, respectively, so they will "only" work correctly
1603      * for atomistic force fields using MD units.
1604      */
1605     const real     massFactorThreshold      = 13.0;
1606     const real     bufferToleranceThreshold = 1e-4;
1607     const int      lincsIterationThreshold  = 2;
1608     const int      lincsOrderThreshold      = 4;
1609     const real     shakeToleranceThreshold  = 0.005*ir->delta_t;
1610
1611     bool           lincsWithSufficientTolerance = (ir->eConstrAlg == econtLINCS && ir->nLincsIter >= lincsIterationThreshold && ir->nProjOrder >= lincsOrderThreshold);
1612     bool           shakeWithSufficientTolerance = (ir->eConstrAlg == econtSHAKE && ir->shake_tol <= 1.1*shakeToleranceThreshold);
1613     if (ir->cutoff_scheme == ecutsVERLET &&
1614         ir->verletbuf_tol <= 1.1*bufferToleranceThreshold &&
1615         (lincsWithSufficientTolerance || shakeWithSufficientTolerance))
1616     {
1617         return;
1618     }
1619
1620     bool haveDecoupledMode = false;
1621     for (const gmx_moltype_t &molt : mtop->moltype)
1622     {
1623         if (haveDecoupledModeInMol(molt, mtop->ffparams.iparams,
1624                                    massFactorThreshold))
1625         {
1626             haveDecoupledMode = true;
1627         }
1628     }
1629
1630     if (haveDecoupledMode)
1631     {
1632         std::string message = gmx::formatString(
1633                     "There are atoms at both ends of an angle, connected by constraints "
1634                     "and with masses that differ by more than a factor of %g. This means "
1635                     "that there are likely dynamic modes that are only very weakly coupled.",
1636                     massFactorThreshold);
1637         if (ir->cutoff_scheme == ecutsVERLET)
1638         {
1639             message += gmx::formatString(
1640                         " To ensure good equipartitioning, you need to either not use "
1641                         "constraints on all bonds (but, if possible, only on bonds involving "
1642                         "hydrogens) or use integrator = %s or decrease one or more tolerances: "
1643                         "verlet-buffer-tolerance <= %g, LINCS iterations >= %d, LINCS order "
1644                         ">= %d or SHAKE tolerance <= %g",
1645                         ei_names[eiSD1],
1646                         bufferToleranceThreshold,
1647                         lincsIterationThreshold, lincsOrderThreshold,
1648                         shakeToleranceThreshold);
1649         }
1650         else
1651         {
1652             message += gmx::formatString(
1653                         " To ensure good equipartitioning, we suggest to switch to the %s "
1654                         "cutoff-scheme, since that allows for better control over the Verlet "
1655                         "buffer size and thus over the energy drift.",
1656                         ecutscheme_names[ecutsVERLET]);
1657         }
1658         warning(wi, message);
1659     }
1660 }
1661
1662 static void set_verlet_buffer(const gmx_mtop_t *mtop,
1663                               t_inputrec       *ir,
1664                               real              buffer_temp,
1665                               matrix            box,
1666                               warninp          *wi)
1667 {
1668     char                   warn_buf[STRLEN];
1669
1670     printf("Determining Verlet buffer for a tolerance of %g kJ/mol/ps at %g K\n", ir->verletbuf_tol, buffer_temp);
1671
1672     /* Calculate the buffer size for simple atom vs atoms list */
1673     VerletbufListSetup listSetup1x1;
1674     listSetup1x1.cluster_size_i = 1;
1675     listSetup1x1.cluster_size_j = 1;
1676     const real rlist_1x1 =
1677         calcVerletBufferSize(*mtop, det(box), *ir, ir->nstlist, ir->nstlist - 1,
1678                              buffer_temp, listSetup1x1);
1679
1680     /* Set the pair-list buffer size in ir */
1681     VerletbufListSetup listSetup4x4 =
1682         verletbufGetSafeListSetup(ListSetupType::CpuNoSimd);
1683     ir->rlist =
1684         calcVerletBufferSize(*mtop, det(box), *ir, ir->nstlist, ir->nstlist - 1,
1685                              buffer_temp, listSetup4x4);
1686
1687     const int n_nonlin_vsite = countNonlinearVsites(*mtop);
1688     if (n_nonlin_vsite > 0)
1689     {
1690         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);
1691         warning_note(wi, warn_buf);
1692     }
1693
1694     printf("Calculated rlist for %dx%d atom pair-list as %.3f nm, buffer size %.3f nm\n",
1695            1, 1, rlist_1x1, rlist_1x1-std::max(ir->rvdw, ir->rcoulomb));
1696
1697     printf("Set rlist, assuming %dx%d atom pair-list, to %.3f nm, buffer size %.3f nm\n",
1698            listSetup4x4.cluster_size_i, listSetup4x4.cluster_size_j,
1699            ir->rlist, ir->rlist-std::max(ir->rvdw, ir->rcoulomb));
1700
1701     printf("Note that mdrun will redetermine rlist based on the actual pair-list setup\n");
1702
1703     if (gmx::square(ir->rlist) >= max_cutoff2(ir->ePBC, box))
1704     {
1705         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)));
1706     }
1707 }
1708
1709 int gmx_grompp(int argc, char *argv[])
1710 {
1711     const char                          *desc[] = {
1712         "[THISMODULE] (the gromacs preprocessor)",
1713         "reads a molecular topology file, checks the validity of the",
1714         "file, expands the topology from a molecular description to an atomic",
1715         "description. The topology file contains information about",
1716         "molecule types and the number of molecules, the preprocessor",
1717         "copies each molecule as needed. ",
1718         "There is no limitation on the number of molecule types. ",
1719         "Bonds and bond-angles can be converted into constraints, separately",
1720         "for hydrogens and heavy atoms.",
1721         "Then a coordinate file is read and velocities can be generated",
1722         "from a Maxwellian distribution if requested.",
1723         "[THISMODULE] also reads parameters for [gmx-mdrun] ",
1724         "(eg. number of MD steps, time step, cut-off), and others such as",
1725         "NEMD parameters, which are corrected so that the net acceleration",
1726         "is zero.",
1727         "Eventually a binary file is produced that can serve as the sole input",
1728         "file for the MD program.[PAR]",
1729
1730         "[THISMODULE] uses the atom names from the topology file. The atom names",
1731         "in the coordinate file (option [TT]-c[tt]) are only read to generate",
1732         "warnings when they do not match the atom names in the topology.",
1733         "Note that the atom names are irrelevant for the simulation as",
1734         "only the atom types are used for generating interaction parameters.[PAR]",
1735
1736         "[THISMODULE] uses a built-in preprocessor to resolve includes, macros, ",
1737         "etc. The preprocessor supports the following keywords::",
1738         "",
1739         "    #ifdef VARIABLE",
1740         "    #ifndef VARIABLE",
1741         "    #else",
1742         "    #endif",
1743         "    #define VARIABLE",
1744         "    #undef VARIABLE",
1745         "    #include \"filename\"",
1746         "    #include <filename>",
1747         "",
1748         "The functioning of these statements in your topology may be modulated by",
1749         "using the following two flags in your [REF].mdp[ref] file::",
1750         "",
1751         "    define = -DVARIABLE1 -DVARIABLE2",
1752         "    include = -I/home/john/doe",
1753         "",
1754         "For further information a C-programming textbook may help you out.",
1755         "Specifying the [TT]-pp[tt] flag will get the pre-processed",
1756         "topology file written out so that you can verify its contents.[PAR]",
1757
1758         "When using position restraints, a file with restraint coordinates",
1759         "must be supplied with [TT]-r[tt] (can be the same file as supplied",
1760         "for [TT]-c[tt]). For free energy calculations, separate reference",
1761         "coordinates for the B topology can be supplied with [TT]-rb[tt],",
1762         "otherwise they will be equal to those of the A topology.[PAR]",
1763
1764         "Starting coordinates can be read from trajectory with [TT]-t[tt].",
1765         "The last frame with coordinates and velocities will be read,",
1766         "unless the [TT]-time[tt] option is used. Only if this information",
1767         "is absent will the coordinates in the [TT]-c[tt] file be used.",
1768         "Note that these velocities will not be used when [TT]gen_vel = yes[tt]",
1769         "in your [REF].mdp[ref] file. An energy file can be supplied with",
1770         "[TT]-e[tt] to read Nose-Hoover and/or Parrinello-Rahman coupling",
1771         "variables.[PAR]",
1772
1773         "[THISMODULE] can be used to restart simulations (preserving",
1774         "continuity) by supplying just a checkpoint file with [TT]-t[tt].",
1775         "However, for simply changing the number of run steps to extend",
1776         "a run, using [gmx-convert-tpr] is more convenient than [THISMODULE].",
1777         "You then supply the old checkpoint file directly to [gmx-mdrun]",
1778         "with [TT]-cpi[tt]. If you wish to change the ensemble or things",
1779         "like output frequency, then supplying the checkpoint file to",
1780         "[THISMODULE] with [TT]-t[tt] along with a new [REF].mdp[ref] file",
1781         "with [TT]-f[tt] is the recommended procedure. Actually preserving",
1782         "the ensemble (if possible) still requires passing the checkpoint",
1783         "file to [gmx-mdrun] [TT]-cpi[tt].[PAR]",
1784
1785         "By default, all bonded interactions which have constant energy due to",
1786         "virtual site constructions will be removed. If this constant energy is",
1787         "not zero, this will result in a shift in the total energy. All bonded",
1788         "interactions can be kept by turning off [TT]-rmvsbds[tt]. Additionally,",
1789         "all constraints for distances which will be constant anyway because",
1790         "of virtual site constructions will be removed. If any constraints remain",
1791         "which involve virtual sites, a fatal error will result.[PAR]",
1792
1793         "To verify your run input file, please take note of all warnings",
1794         "on the screen, and correct where necessary. Do also look at the contents",
1795         "of the [TT]mdout.mdp[tt] file; this contains comment lines, as well as",
1796         "the input that [THISMODULE] has read. If in doubt, you can start [THISMODULE]",
1797         "with the [TT]-debug[tt] option which will give you more information",
1798         "in a file called [TT]grompp.log[tt] (along with real debug info). You",
1799         "can see the contents of the run input file with the [gmx-dump]",
1800         "program. [gmx-check] can be used to compare the contents of two",
1801         "run input files.[PAR]",
1802
1803         "The [TT]-maxwarn[tt] option can be used to override warnings printed",
1804         "by [THISMODULE] that otherwise halt output. In some cases, warnings are",
1805         "harmless, but usually they are not. The user is advised to carefully",
1806         "interpret the output messages before attempting to bypass them with",
1807         "this option."
1808     };
1809     t_gromppopts                        *opts;
1810     std::vector<MoleculeInformation>     mi;
1811     std::unique_ptr<MoleculeInformation> intermolecular_interactions;
1812     int                                  nvsite, comb;
1813     real                                 fudgeQQ;
1814     double                               reppow;
1815     const char                          *mdparin;
1816     bool                                 bNeedVel, bGenVel;
1817     gmx_bool                             have_atomnumber;
1818     gmx_output_env_t                    *oenv;
1819     gmx_bool                             bVerbose = FALSE;
1820     warninp                             *wi;
1821     char                                 warn_buf[STRLEN];
1822
1823     t_filenm                             fnm[] = {
1824         { efMDP, nullptr,  nullptr,        ffREAD  },
1825         { efMDP, "-po", "mdout",     ffWRITE },
1826         { efSTX, "-c",  nullptr,        ffREAD  },
1827         { efSTX, "-r",  "restraint", ffOPTRD },
1828         { efSTX, "-rb", "restraint", ffOPTRD },
1829         { efNDX, nullptr,  nullptr,        ffOPTRD },
1830         { efTOP, nullptr,  nullptr,        ffREAD  },
1831         { efTOP, "-pp", "processed", ffOPTWR },
1832         { efTPR, "-o",  nullptr,        ffWRITE },
1833         { efTRN, "-t",  nullptr,        ffOPTRD },
1834         { efEDR, "-e",  nullptr,        ffOPTRD },
1835         /* This group is needed by the VMD viewer as the start configuration for IMD sessions: */
1836         { efGRO, "-imd", "imdgroup", ffOPTWR },
1837         { efTRN, "-ref", "rotref",   ffOPTRW }
1838     };
1839 #define NFILE asize(fnm)
1840
1841     /* Command line options */
1842     gmx_bool            bRenum   = TRUE;
1843     gmx_bool            bRmVSBds = TRUE, bZero = FALSE;
1844     int                 i, maxwarn = 0;
1845     real                fr_time = -1;
1846     t_pargs             pa[]    = {
1847         { "-v",       FALSE, etBOOL, {&bVerbose},
1848           "Be loud and noisy" },
1849         { "-time",    FALSE, etREAL, {&fr_time},
1850           "Take frame at or first after this time." },
1851         { "-rmvsbds", FALSE, etBOOL, {&bRmVSBds},
1852           "Remove constant bonded interactions with virtual sites" },
1853         { "-maxwarn", FALSE, etINT,  {&maxwarn},
1854           "Number of allowed warnings during input processing. Not for normal use and may generate unstable systems" },
1855         { "-zero",    FALSE, etBOOL, {&bZero},
1856           "Set parameters for bonded interactions without defaults to zero instead of generating an error" },
1857         { "-renum",   FALSE, etBOOL, {&bRenum},
1858           "Renumber atomtypes and minimize number of atomtypes" }
1859     };
1860
1861     /* Parse the command line */
1862     if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa,
1863                            asize(desc), desc, 0, nullptr, &oenv))
1864     {
1865         return 0;
1866     }
1867
1868     /* Initiate some variables */
1869     gmx::MDModules mdModules;
1870     t_inputrec     irInstance;
1871     t_inputrec    *ir = &irInstance;
1872     snew(opts, 1);
1873     snew(opts->include, STRLEN);
1874     snew(opts->define, STRLEN);
1875
1876     wi = init_warning(TRUE, maxwarn);
1877
1878     /* PARAMETER file processing */
1879     mdparin = opt2fn("-f", NFILE, fnm);
1880     set_warning_line(wi, mdparin, -1);
1881     try
1882     {
1883         get_ir(mdparin, opt2fn("-po", NFILE, fnm), &mdModules, ir, opts, WriteMdpHeader::yes, wi);
1884     }
1885     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
1886
1887     if (bVerbose)
1888     {
1889         fprintf(stderr, "checking input for internal consistency...\n");
1890     }
1891     check_ir(mdparin, ir, opts, wi);
1892
1893     if (ir->ld_seed == -1)
1894     {
1895         ir->ld_seed = static_cast<int>(gmx::makeRandomSeed());
1896         fprintf(stderr, "Setting the LD random seed to %" PRId64 "\n", ir->ld_seed);
1897     }
1898
1899     if (ir->expandedvals->lmc_seed == -1)
1900     {
1901         ir->expandedvals->lmc_seed = static_cast<int>(gmx::makeRandomSeed());
1902         fprintf(stderr, "Setting the lambda MC random seed to %d\n", ir->expandedvals->lmc_seed);
1903     }
1904
1905     bNeedVel = EI_STATE_VELOCITY(ir->eI);
1906     bGenVel  = (bNeedVel && opts->bGenVel);
1907     if (bGenVel && ir->bContinuation)
1908     {
1909         sprintf(warn_buf,
1910                 "Generating velocities is inconsistent with attempting "
1911                 "to continue a previous run. Choose only one of "
1912                 "gen-vel = yes and continuation = yes.");
1913         warning_error(wi, warn_buf);
1914     }
1915
1916     std::array<InteractionTypeParameters, F_NRE> plist;
1917     gmx_mtop_t             sys;
1918     PreprocessingAtomTypes atypes;
1919     if (debug)
1920     {
1921         pr_symtab(debug, 0, "Just opened", &sys.symtab);
1922     }
1923
1924     const char *fn = ftp2fn(efTOP, NFILE, fnm);
1925     if (!gmx_fexist(fn))
1926     {
1927         gmx_fatal(FARGS, "%s does not exist", fn);
1928     }
1929
1930     t_state state;
1931     new_status(fn, opt2fn_null("-pp", NFILE, fnm), opt2fn("-c", NFILE, fnm),
1932                opts, ir, bZero, bGenVel, bVerbose, &state,
1933                &atypes, &sys, &mi, &intermolecular_interactions,
1934                plist, &comb, &reppow, &fudgeQQ,
1935                opts->bMorse,
1936                wi);
1937
1938     if (debug)
1939     {
1940         pr_symtab(debug, 0, "After new_status", &sys.symtab);
1941     }
1942
1943     nvsite = 0;
1944     /* set parameters for virtual site construction (not for vsiten) */
1945     for (size_t mt = 0; mt < sys.moltype.size(); mt++)
1946     {
1947         nvsite +=
1948             set_vsites(bVerbose, &sys.moltype[mt].atoms, &atypes, mi[mt].plist);
1949     }
1950     /* now throw away all obsolete bonds, angles and dihedrals: */
1951     /* note: constraints are ALWAYS removed */
1952     if (nvsite)
1953     {
1954         for (size_t mt = 0; mt < sys.moltype.size(); mt++)
1955         {
1956             clean_vsite_bondeds(mi[mt].plist, sys.moltype[mt].atoms.nr, bRmVSBds);
1957         }
1958     }
1959
1960     if (ir->cutoff_scheme == ecutsVERLET)
1961     {
1962         fprintf(stderr, "Removing all charge groups because cutoff-scheme=%s\n",
1963                 ecutscheme_names[ir->cutoff_scheme]);
1964
1965         /* Remove all charge groups */
1966         gmx_mtop_remove_chargegroups(&sys);
1967     }
1968
1969     if ((count_constraints(&sys, mi, wi) != 0) && (ir->eConstrAlg == econtSHAKE))
1970     {
1971         if (ir->eI == eiCG || ir->eI == eiLBFGS)
1972         {
1973             sprintf(warn_buf, "Can not do %s with %s, use %s",
1974                     EI(ir->eI), econstr_names[econtSHAKE], econstr_names[econtLINCS]);
1975             warning_error(wi, warn_buf);
1976         }
1977         if (ir->bPeriodicMols)
1978         {
1979             sprintf(warn_buf, "Can not do periodic molecules with %s, use %s",
1980                     econstr_names[econtSHAKE], econstr_names[econtLINCS]);
1981             warning_error(wi, warn_buf);
1982         }
1983     }
1984
1985     if (EI_SD (ir->eI) &&  ir->etc != etcNO)
1986     {
1987         warning_note(wi, "Temperature coupling is ignored with SD integrators.");
1988     }
1989
1990     /* If we are doing QM/MM, check that we got the atom numbers */
1991     have_atomnumber = TRUE;
1992     for (int i = 0; i < gmx::ssize(atypes); i++)
1993     {
1994         have_atomnumber = have_atomnumber && (atypes.atomNumberFromAtomType(i) >= 0);
1995     }
1996     if (!have_atomnumber && ir->bQMMM)
1997     {
1998         warning_error(wi,
1999                       "\n"
2000                       "It appears as if you are trying to run a QM/MM calculation, but the force\n"
2001                       "field you are using does not contain atom numbers fields. This is an\n"
2002                       "optional field (introduced in GROMACS 3.3) for general runs, but mandatory\n"
2003                       "for QM/MM. The good news is that it is easy to add - put the atom number as\n"
2004                       "an integer just before the mass column in ffXXXnb.itp.\n"
2005                       "NB: United atoms have the same atom numbers as normal ones.\n\n");
2006     }
2007
2008     /* Check for errors in the input now, since they might cause problems
2009      * during processing further down.
2010      */
2011     check_warning_error(wi, FARGS);
2012
2013     if (nint_ftype(&sys, mi, F_POSRES) > 0 ||
2014         nint_ftype(&sys, mi, F_FBPOSRES) > 0)
2015     {
2016         if (ir->epc == epcPARRINELLORAHMAN || ir->epc == epcMTTK)
2017         {
2018             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.",
2019                     EPCOUPLTYPE(ir->epc), EPCOUPLTYPE(epcBERENDSEN));
2020             warning_note(wi, warn_buf);
2021         }
2022
2023         const char *fn = opt2fn("-r", NFILE, fnm);
2024         const char *fnB;
2025
2026         if (!gmx_fexist(fn))
2027         {
2028             gmx_fatal(FARGS,
2029                       "Cannot find position restraint file %s (option -r).\n"
2030                       "From GROMACS-2018, you need to specify the position restraint "
2031                       "coordinate files explicitly to avoid mistakes, although you can "
2032                       "still use the same file as you specify for the -c option.", fn);
2033         }
2034
2035         if (opt2bSet("-rb", NFILE, fnm))
2036         {
2037             fnB = opt2fn("-rb", NFILE, fnm);
2038             if (!gmx_fexist(fnB))
2039             {
2040                 gmx_fatal(FARGS,
2041                           "Cannot find B-state position restraint file %s (option -rb).\n"
2042                           "From GROMACS-2018, you need to specify the position restraint "
2043                           "coordinate files explicitly to avoid mistakes, although you can "
2044                           "still use the same file as you specify for the -c option.", fn);
2045             }
2046         }
2047         else
2048         {
2049             fnB = fn;
2050         }
2051
2052         if (bVerbose)
2053         {
2054             fprintf(stderr, "Reading position restraint coords from %s", fn);
2055             if (strcmp(fn, fnB) == 0)
2056             {
2057                 fprintf(stderr, "\n");
2058             }
2059             else
2060             {
2061                 fprintf(stderr, " and %s\n", fnB);
2062             }
2063         }
2064         gen_posres(&sys, mi, fn, fnB,
2065                    ir->refcoord_scaling, ir->ePBC,
2066                    ir->posres_com, ir->posres_comB,
2067                    wi);
2068     }
2069
2070     /* If we are using CMAP, setup the pre-interpolation grid */
2071     if (plist[F_CMAP].ncmap() > 0)
2072     {
2073         init_cmap_grid(&sys.ffparams.cmap_grid, plist[F_CMAP].cmapAngles, plist[F_CMAP].cmakeGridSpacing);
2074         setup_cmap(plist[F_CMAP].cmakeGridSpacing, plist[F_CMAP].cmapAngles, plist[F_CMAP].cmap, &sys.ffparams.cmap_grid);
2075     }
2076
2077     set_wall_atomtype(&atypes, opts, ir, wi);
2078     if (bRenum)
2079     {
2080         atypes.renumberTypes(plist, &sys, ir->wall_atomtype, bVerbose);
2081     }
2082
2083     if (ir->implicit_solvent)
2084     {
2085         gmx_fatal(FARGS, "Implicit solvation is no longer supported");
2086     }
2087
2088     /* PELA: Copy the atomtype data to the topology atomtype list */
2089     atypes.copyTot_atomtypes(&(sys.atomtypes));
2090
2091     if (debug)
2092     {
2093         pr_symtab(debug, 0, "After atype.renumberTypes", &sys.symtab);
2094     }
2095
2096     if (bVerbose)
2097     {
2098         fprintf(stderr, "converting bonded parameters...\n");
2099     }
2100
2101     const int ntype = atypes.size();
2102     convertInteractionTypeParameters(ntype, plist, mi, intermolecular_interactions.get(),
2103                                      comb, reppow, fudgeQQ, &sys);
2104
2105     if (debug)
2106     {
2107         pr_symtab(debug, 0, "After converInteractionTypeParameters", &sys.symtab);
2108     }
2109
2110     /* set ptype to VSite for virtual sites */
2111     for (gmx_moltype_t &moltype : sys.moltype)
2112     {
2113         set_vsites_ptype(FALSE, &moltype);
2114     }
2115     if (debug)
2116     {
2117         pr_symtab(debug, 0, "After virtual sites", &sys.symtab);
2118     }
2119     /* Check velocity for virtual sites and shells */
2120     if (bGenVel)
2121     {
2122         check_vel(&sys, state.v.rvec_array());
2123     }
2124
2125     /* check for shells and inpurecs */
2126     check_shells_inputrec(&sys, ir, wi);
2127
2128     /* check masses */
2129     check_mol(&sys, wi);
2130
2131     checkForUnboundAtoms(&sys, bVerbose, wi);
2132
2133     for (const gmx_moltype_t &moltype : sys.moltype)
2134     {
2135         check_cg_sizes(ftp2fn(efTOP, NFILE, fnm), &moltype.cgs, wi);
2136     }
2137
2138     if (EI_DYNAMICS(ir->eI) && ir->eI != eiBD)
2139     {
2140         check_bonds_timestep(&sys, ir->delta_t, wi);
2141     }
2142
2143     checkDecoupledModeAccuracy(&sys, ir, wi);
2144
2145     if (EI_ENERGY_MINIMIZATION(ir->eI) && 0 == ir->nsteps)
2146     {
2147         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.");
2148     }
2149
2150     check_warning_error(wi, FARGS);
2151
2152     if (bVerbose)
2153     {
2154         fprintf(stderr, "initialising group options...\n");
2155     }
2156     do_index(mdparin, ftp2fn_null(efNDX, NFILE, fnm),
2157              &sys, bVerbose, ir,
2158              wi);
2159
2160     if (ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0)
2161     {
2162         if (EI_DYNAMICS(ir->eI) && inputrec2nboundeddim(ir) == 3)
2163         {
2164             real buffer_temp;
2165
2166             if (EI_MD(ir->eI) && ir->etc == etcNO)
2167             {
2168                 if (bGenVel)
2169                 {
2170                     buffer_temp = opts->tempi;
2171                 }
2172                 else
2173                 {
2174                     buffer_temp = calc_temp(&sys, ir, state.v.rvec_array());
2175                 }
2176                 if (buffer_temp > 0)
2177                 {
2178                     sprintf(warn_buf, "NVE simulation: will use the initial temperature of %.3f K for determining the Verlet buffer size", buffer_temp);
2179                     warning_note(wi, warn_buf);
2180                 }
2181                 else
2182                 {
2183                     sprintf(warn_buf, "NVE simulation with an initial temperature of zero: will use a Verlet buffer of %d%%. Check your energy drift!",
2184                             gmx::roundToInt(verlet_buffer_ratio_NVE_T0*100));
2185                     warning_note(wi, warn_buf);
2186                 }
2187             }
2188             else
2189             {
2190                 buffer_temp = get_max_reference_temp(ir, wi);
2191             }
2192
2193             if (EI_MD(ir->eI) && ir->etc == etcNO && buffer_temp == 0)
2194             {
2195                 /* NVE with initial T=0: we add a fixed ratio to rlist.
2196                  * Since we don't actually use verletbuf_tol, we set it to -1
2197                  * so it can't be misused later.
2198                  */
2199                 ir->rlist         *= 1.0 + verlet_buffer_ratio_NVE_T0;
2200                 ir->verletbuf_tol  = -1;
2201             }
2202             else
2203             {
2204                 /* We warn for NVE simulations with a drift tolerance that
2205                  * might result in a 1(.1)% drift over the total run-time.
2206                  * Note that we can't warn when nsteps=0, since we don't
2207                  * know how many steps the user intends to run.
2208                  */
2209                 if (EI_MD(ir->eI) && ir->etc == etcNO && ir->nstlist > 1 &&
2210                     ir->nsteps > 0)
2211                 {
2212                     const real driftTolerance = 0.01;
2213                     /* We use 2 DOF per atom = 2kT pot+kin energy,
2214                      * to be on the safe side with constraints.
2215                      */
2216                     const real totalEnergyDriftPerAtomPerPicosecond = 2*BOLTZ*buffer_temp/(ir->nsteps*ir->delta_t);
2217
2218                     if (ir->verletbuf_tol > 1.1*driftTolerance*totalEnergyDriftPerAtomPerPicosecond)
2219                     {
2220                         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.",
2221                                 ir->verletbuf_tol, ir->nsteps*ir->delta_t,
2222                                 gmx::roundToInt(ir->verletbuf_tol/totalEnergyDriftPerAtomPerPicosecond*100),
2223                                 gmx::roundToInt(100*driftTolerance),
2224                                 driftTolerance*totalEnergyDriftPerAtomPerPicosecond);
2225                         warning_note(wi, warn_buf);
2226                     }
2227                 }
2228
2229                 set_verlet_buffer(&sys, ir, buffer_temp, state.box, wi);
2230             }
2231         }
2232     }
2233
2234     /* Init the temperature coupling state */
2235     init_gtc_state(&state, ir->opts.ngtc, 0, ir->opts.nhchainlength); /* need to add nnhpres here? */
2236
2237     if (bVerbose)
2238     {
2239         fprintf(stderr, "Checking consistency between energy and charge groups...\n");
2240     }
2241     check_eg_vs_cg(&sys);
2242
2243     if (debug)
2244     {
2245         pr_symtab(debug, 0, "After index", &sys.symtab);
2246     }
2247
2248     triple_check(mdparin, ir, &sys, wi);
2249     close_symtab(&sys.symtab);
2250     if (debug)
2251     {
2252         pr_symtab(debug, 0, "After close", &sys.symtab);
2253     }
2254
2255     /* make exclusions between QM atoms and remove charges if needed */
2256     if (ir->bQMMM)
2257     {
2258         if (ir->QMMMscheme == eQMMMschemenormal && ir->ns_type == ensSIMPLE)
2259         {
2260             gmx_fatal(FARGS, "electrostatic embedding only works with grid neighboursearching, use ns-type=grid instead\n");
2261         }
2262         else
2263         {
2264             generate_qmexcl(&sys, ir, wi, GmxQmmmMode::GMX_QMMM_ORIGINAL);
2265         }
2266         if (ir->QMMMscheme != eQMMMschemeoniom)
2267         {
2268             std::vector<int> qmmmAtoms = qmmmAtomIndices(*ir, sys);
2269             removeQmmmAtomCharges(&sys, qmmmAtoms);
2270         }
2271     }
2272
2273     if (ir->eI == eiMimic)
2274     {
2275         generate_qmexcl(&sys, ir, wi, GmxQmmmMode::GMX_QMMM_MIMIC);
2276     }
2277
2278     if (ftp2bSet(efTRN, NFILE, fnm))
2279     {
2280         if (bVerbose)
2281         {
2282             fprintf(stderr, "getting data from old trajectory ...\n");
2283         }
2284         cont_status(ftp2fn(efTRN, NFILE, fnm), ftp2fn_null(efEDR, NFILE, fnm),
2285                     bNeedVel, bGenVel, fr_time, ir, &state, &sys, oenv);
2286     }
2287
2288     if (ir->ePBC == epbcXY && ir->nwall != 2)
2289     {
2290         clear_rvec(state.box[ZZ]);
2291     }
2292
2293     if (ir->cutoff_scheme != ecutsVERLET && ir->rlist > 0)
2294     {
2295         set_warning_line(wi, mdparin, -1);
2296         check_chargegroup_radii(&sys, ir, state.x.rvec_array(), wi);
2297     }
2298
2299     if (EEL_FULL(ir->coulombtype) || EVDW_PME(ir->vdwtype))
2300     {
2301         /* Calculate the optimal grid dimensions */
2302         matrix          scaledBox;
2303         EwaldBoxZScaler boxScaler(*ir);
2304         boxScaler.scaleBox(state.box, scaledBox);
2305
2306         if (ir->nkx > 0 && ir->nky > 0 && ir->nkz > 0)
2307         {
2308             /* Mark fourier_spacing as not used */
2309             ir->fourier_spacing = 0;
2310         }
2311         else if (ir->nkx != 0 && ir->nky != 0 && ir->nkz != 0)
2312         {
2313             set_warning_line(wi, mdparin, -1);
2314             warning_error(wi, "Some of the Fourier grid sizes are set, but all of them need to be set.");
2315         }
2316         const int minGridSize = minimalPmeGridSize(ir->pme_order);
2317         calcFftGrid(stdout, scaledBox, ir->fourier_spacing, minGridSize,
2318                     &(ir->nkx), &(ir->nky), &(ir->nkz));
2319         if (ir->nkx < minGridSize ||
2320             ir->nky < minGridSize ||
2321             ir->nkz < minGridSize)
2322         {
2323             warning_error(wi, "The PME grid size should be >= 2*(pme-order - 1); either manually increase the grid size or decrease pme-order");
2324         }
2325     }
2326
2327     /* MRS: eventually figure out better logic for initializing the fep
2328        values that makes declaring the lambda and declaring the state not
2329        potentially conflict if not handled correctly. */
2330     if (ir->efep != efepNO)
2331     {
2332         state.fep_state = ir->fepvals->init_fep_state;
2333         for (i = 0; i < efptNR; i++)
2334         {
2335             /* init_lambda trumps state definitions*/
2336             if (ir->fepvals->init_lambda >= 0)
2337             {
2338                 state.lambda[i] = ir->fepvals->init_lambda;
2339             }
2340             else
2341             {
2342                 if (ir->fepvals->all_lambda[i] == nullptr)
2343                 {
2344                     gmx_fatal(FARGS, "Values of lambda not set for a free energy calculation!");
2345                 }
2346                 else
2347                 {
2348                     state.lambda[i] = ir->fepvals->all_lambda[i][state.fep_state];
2349                 }
2350             }
2351         }
2352     }
2353
2354     pull_t *pull = nullptr;
2355
2356     if (ir->bPull)
2357     {
2358         pull = set_pull_init(ir, &sys, state.x.rvec_array(), state.box, state.lambda[efptMASS], wi);
2359     }
2360
2361     /* Modules that supply external potential for pull coordinates
2362      * should register those potentials here. finish_pull() will check
2363      * that providers have been registerd for all external potentials.
2364      */
2365
2366     if (ir->bDoAwh)
2367     {
2368         setStateDependentAwhParams(ir->awhParams, ir->pull, pull,
2369                                    state.box, ir->ePBC, &ir->opts, wi);
2370     }
2371
2372     if (ir->bPull)
2373     {
2374         finish_pull(pull);
2375     }
2376
2377     if (ir->bRot)
2378     {
2379         set_reference_positions(ir->rot, state.x.rvec_array(), state.box,
2380                                 opt2fn("-ref", NFILE, fnm), opt2bSet("-ref", NFILE, fnm),
2381                                 wi);
2382     }
2383
2384     /*  reset_multinr(sys); */
2385
2386     if (EEL_PME(ir->coulombtype))
2387     {
2388         float ratio = pme_load_estimate(&sys, ir, state.box);
2389         fprintf(stderr, "Estimate for the relative computational load of the PME mesh part: %.2f\n", ratio);
2390         /* With free energy we might need to do PME both for the A and B state
2391          * charges. This will double the cost, but the optimal performance will
2392          * then probably be at a slightly larger cut-off and grid spacing.
2393          */
2394         if ((ir->efep == efepNO && ratio > 1.0/2.0) ||
2395             (ir->efep != efepNO && ratio > 2.0/3.0))
2396         {
2397             warning_note(wi,
2398                          "The optimal PME mesh load for parallel simulations is below 0.5\n"
2399                          "and for highly parallel simulations between 0.25 and 0.33,\n"
2400                          "for higher performance, increase the cut-off and the PME grid spacing.\n");
2401             if (ir->efep != efepNO)
2402             {
2403                 warning_note(wi,
2404                              "For free energy simulations, the optimal load limit increases from 0.5 to 0.667\n");
2405             }
2406         }
2407     }
2408
2409     {
2410         char   warn_buf[STRLEN];
2411         double cio = compute_io(ir, sys.natoms, sys.groups, F_NRE, 1);
2412         sprintf(warn_buf, "This run will generate roughly %.0f Mb of data", cio);
2413         if (cio > 2000)
2414         {
2415             set_warning_line(wi, mdparin, -1);
2416             warning_note(wi, warn_buf);
2417         }
2418         else
2419         {
2420             printf("%s\n", warn_buf);
2421         }
2422     }
2423
2424     if (bVerbose)
2425     {
2426         fprintf(stderr, "writing run input file...\n");
2427     }
2428
2429     done_warning(wi, FARGS);
2430     write_tpx_state(ftp2fn(efTPR, NFILE, fnm), ir, &state, &sys);
2431
2432     /* Output IMD group, if bIMD is TRUE */
2433     gmx::write_IMDgroup_to_file(ir->bIMD, ir, &state, &sys, NFILE, fnm);
2434
2435     sfree(opts->define);
2436     sfree(opts->include);
2437     sfree(opts);
2438     for (auto &mol : mi)
2439     {
2440         // Some of the contents of molinfo have been stolen, so
2441         // fullCleanUp can't be called.
2442         mol.partialCleanUp();
2443     }
2444     done_inputrec_strings();
2445     output_env_done(oenv);
2446
2447     return 0;
2448 }