src/gromacs/gmxpreprocess/grompp.cpp

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