src/gromacs/gmxpreprocess/grompp.cpp

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