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