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            gpp_atomtype *atype, 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, nrmols, 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                        atype, &nrmols, 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, atype);
 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(gpp_atomtype *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] = get_atomtype_type(opts->wall_atomtype[i], at);
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     gpp_atomtype                        *atype;
1696     int                                  nvsite, comb;
1697     real                                 fudgeQQ;
1698     double                               reppow;
1699     const char                          *mdparin;
1700     int                                  ntype;
1701     bool                                 bNeedVel, bGenVel;
1702     gmx_bool                             have_atomnumber;
1703     gmx_output_env_t                    *oenv;
1704     gmx_bool                             bVerbose = FALSE;
1705     warninp                             *wi;
1706     char                                 warn_buf[STRLEN];
1707
1708     t_filenm                             fnm[] = {
1709         { efMDP, nullptr,  nullptr,        ffREAD  },
1710         { efMDP, "-po", "mdout",     ffWRITE },
1711         { efSTX, "-c",  nullptr,        ffREAD  },
1712         { efSTX, "-r",  "restraint", ffOPTRD },
1713         { efSTX, "-rb", "restraint", ffOPTRD },
1714         { efNDX, nullptr,  nullptr,        ffOPTRD },
1715         { efTOP, nullptr,  nullptr,        ffREAD  },
1716         { efTOP, "-pp", "processed", ffOPTWR },
1717         { efTPR, "-o",  nullptr,        ffWRITE },
1718         { efTRN, "-t",  nullptr,        ffOPTRD },
1719         { efEDR, "-e",  nullptr,        ffOPTRD },
1720         /* This group is needed by the VMD viewer as the start configuration for IMD sessions: */
1721         { efGRO, "-imd", "imdgroup", ffOPTWR },
1722         { efTRN, "-ref", "rotref",   ffOPTRW }
1723     };
1724 #define NFILE asize(fnm)
1725
1726     /* Command line options */
1727     gmx_bool            bRenum   = TRUE;
1728     gmx_bool            bRmVSBds = TRUE, bZero = FALSE;
1729     int                 i, maxwarn = 0;
1730     real                fr_time = -1;
1731     t_pargs             pa[]    = {
1732         { "-v",       FALSE, etBOOL, {&bVerbose},
1733           "Be loud and noisy" },
1734         { "-time",    FALSE, etREAL, {&fr_time},
1735           "Take frame at or first after this time." },
1736         { "-rmvsbds", FALSE, etBOOL, {&bRmVSBds},
1737           "Remove constant bonded interactions with virtual sites" },
1738         { "-maxwarn", FALSE, etINT,  {&maxwarn},
1739           "Number of allowed warnings during input processing. Not for normal use and may generate unstable systems" },
1740         { "-zero",    FALSE, etBOOL, {&bZero},
1741           "Set parameters for bonded interactions without defaults to zero instead of generating an error" },
1742         { "-renum",   FALSE, etBOOL, {&bRenum},
1743           "Renumber atomtypes and minimize number of atomtypes" }
1744     };
1745
1746     /* Parse the command line */
1747     if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa,
1748                            asize(desc), desc, 0, nullptr, &oenv))
1749     {
1750         return 0;
1751     }
1752
1753     /* Initiate some variables */
1754     gmx::MDModules mdModules;
1755     t_inputrec     irInstance;
1756     t_inputrec    *ir = &irInstance;
1757     snew(opts, 1);
1758     snew(opts->include, STRLEN);
1759     snew(opts->define, STRLEN);
1760
1761     wi = init_warning(TRUE, maxwarn);
1762
1763     /* PARAMETER file processing */
1764     mdparin = opt2fn("-f", NFILE, fnm);
1765     set_warning_line(wi, mdparin, -1);
1766     try
1767     {
1768         get_ir(mdparin, opt2fn("-po", NFILE, fnm), &mdModules, ir, opts, WriteMdpHeader::yes, wi);
1769     }
1770     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
1771
1772     if (bVerbose)
1773     {
1774         fprintf(stderr, "checking input for internal consistency...\n");
1775     }
1776     check_ir(mdparin, ir, opts, wi);
1777
1778     if (ir->ld_seed == -1)
1779     {
1780         ir->ld_seed = static_cast<int>(gmx::makeRandomSeed());
1781         fprintf(stderr, "Setting the LD random seed to %" PRId64 "\n", ir->ld_seed);
1782     }
1783
1784     if (ir->expandedvals->lmc_seed == -1)
1785     {
1786         ir->expandedvals->lmc_seed = static_cast<int>(gmx::makeRandomSeed());
1787         fprintf(stderr, "Setting the lambda MC random seed to %d\n", ir->expandedvals->lmc_seed);
1788     }
1789
1790     bNeedVel = EI_STATE_VELOCITY(ir->eI);
1791     bGenVel  = (bNeedVel && opts->bGenVel);
1792     if (bGenVel && ir->bContinuation)
1793     {
1794         sprintf(warn_buf,
1795                 "Generating velocities is inconsistent with attempting "
1796                 "to continue a previous run. Choose only one of "
1797                 "gen-vel = yes and continuation = yes.");
1798         warning_error(wi, warn_buf);
1799     }
1800
1801     std::array<InteractionTypeParameters, F_NRE> plist;
1802     init_plist(plist);
1803     gmx_mtop_t sys;
1804     atype = init_atomtype();
1805     if (debug)
1806     {
1807         pr_symtab(debug, 0, "Just opened", &sys.symtab);
1808     }
1809
1810     const char *fn = ftp2fn(efTOP, NFILE, fnm);
1811     if (!gmx_fexist(fn))
1812     {
1813         gmx_fatal(FARGS, "%s does not exist", fn);
1814     }
1815
1816     t_state state;
1817     new_status(fn, opt2fn_null("-pp", NFILE, fnm), opt2fn("-c", NFILE, fnm),
1818                opts, ir, bZero, bGenVel, bVerbose, &state,
1819                atype, &sys, &mi, &intermolecular_interactions,
1820                plist, &comb, &reppow, &fudgeQQ,
1821                opts->bMorse,
1822                wi);
1823
1824     if (debug)
1825     {
1826         pr_symtab(debug, 0, "After new_status", &sys.symtab);
1827     }
1828
1829     nvsite = 0;
1830     /* set parameters for virtual site construction (not for vsiten) */
1831     for (size_t mt = 0; mt < sys.moltype.size(); mt++)
1832     {
1833         nvsite +=
1834             set_vsites(bVerbose, &sys.moltype[mt].atoms, atype, mi[mt].plist);
1835     }
1836     /* now throw away all obsolete bonds, angles and dihedrals: */
1837     /* note: constraints are ALWAYS removed */
1838     if (nvsite)
1839     {
1840         for (size_t mt = 0; mt < sys.moltype.size(); mt++)
1841         {
1842             clean_vsite_bondeds(mi[mt].plist, sys.moltype[mt].atoms.nr, bRmVSBds);
1843         }
1844     }
1845
1846     if (ir->cutoff_scheme == ecutsVERLET)
1847     {
1848         fprintf(stderr, "Removing all charge groups because cutoff-scheme=%s\n",
1849                 ecutscheme_names[ir->cutoff_scheme]);
1850
1851         /* Remove all charge groups */
1852         gmx_mtop_remove_chargegroups(&sys);
1853     }
1854
1855     if ((count_constraints(&sys, mi, wi) != 0) && (ir->eConstrAlg == econtSHAKE))
1856     {
1857         if (ir->eI == eiCG || ir->eI == eiLBFGS)
1858         {
1859             sprintf(warn_buf, "Can not do %s with %s, use %s",
1860                     EI(ir->eI), econstr_names[econtSHAKE], econstr_names[econtLINCS]);
1861             warning_error(wi, warn_buf);
1862         }
1863         if (ir->bPeriodicMols)
1864         {
1865             sprintf(warn_buf, "Can not do periodic molecules with %s, use %s",
1866                     econstr_names[econtSHAKE], econstr_names[econtLINCS]);
1867             warning_error(wi, warn_buf);
1868         }
1869     }
1870
1871     if (EI_SD (ir->eI) &&  ir->etc != etcNO)
1872     {
1873         warning_note(wi, "Temperature coupling is ignored with SD integrators.");
1874     }
1875
1876     /* If we are doing QM/MM, check that we got the atom numbers */
1877     have_atomnumber = TRUE;
1878     for (i = 0; i < get_atomtype_ntypes(atype); i++)
1879     {
1880         have_atomnumber = have_atomnumber && (get_atomtype_atomnumber(i, atype) >= 0);
1881     }
1882     if (!have_atomnumber && ir->bQMMM)
1883     {
1884         warning_error(wi,
1885                       "\n"
1886                       "It appears as if you are trying to run a QM/MM calculation, but the force\n"
1887                       "field you are using does not contain atom numbers fields. This is an\n"
1888                       "optional field (introduced in GROMACS 3.3) for general runs, but mandatory\n"
1889                       "for QM/MM. The good news is that it is easy to add - put the atom number as\n"
1890                       "an integer just before the mass column in ffXXXnb.itp.\n"
1891                       "NB: United atoms have the same atom numbers as normal ones.\n\n");
1892     }
1893
1894     /* Check for errors in the input now, since they might cause problems
1895      * during processing further down.
1896      */
1897     check_warning_error(wi, FARGS);
1898
1899     if (nint_ftype(&sys, mi, F_POSRES) > 0 ||
1900         nint_ftype(&sys, mi, F_FBPOSRES) > 0)
1901     {
1902         if (ir->epc == epcPARRINELLORAHMAN || ir->epc == epcMTTK)
1903         {
1904             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.",
1905                     EPCOUPLTYPE(ir->epc), EPCOUPLTYPE(epcBERENDSEN));
1906             warning_note(wi, warn_buf);
1907         }
1908
1909         const char *fn = opt2fn("-r", NFILE, fnm);
1910         const char *fnB;
1911
1912         if (!gmx_fexist(fn))
1913         {
1914             gmx_fatal(FARGS,
1915                       "Cannot find position restraint file %s (option -r).\n"
1916                       "From GROMACS-2018, you need to specify the position restraint "
1917                       "coordinate files explicitly to avoid mistakes, although you can "
1918                       "still use the same file as you specify for the -c option.", fn);
1919         }
1920
1921         if (opt2bSet("-rb", NFILE, fnm))
1922         {
1923             fnB = opt2fn("-rb", NFILE, fnm);
1924             if (!gmx_fexist(fnB))
1925             {
1926                 gmx_fatal(FARGS,
1927                           "Cannot find B-state position restraint file %s (option -rb).\n"
1928                           "From GROMACS-2018, you need to specify the position restraint "
1929                           "coordinate files explicitly to avoid mistakes, although you can "
1930                           "still use the same file as you specify for the -c option.", fn);
1931             }
1932         }
1933         else
1934         {
1935             fnB = fn;
1936         }
1937
1938         if (bVerbose)
1939         {
1940             fprintf(stderr, "Reading position restraint coords from %s", fn);
1941             if (strcmp(fn, fnB) == 0)
1942             {
1943                 fprintf(stderr, "\n");
1944             }
1945             else
1946             {
1947                 fprintf(stderr, " and %s\n", fnB);
1948             }
1949         }
1950         gen_posres(&sys, mi, fn, fnB,
1951                    ir->refcoord_scaling, ir->ePBC,
1952                    ir->posres_com, ir->posres_comB,
1953                    wi);
1954     }
1955
1956     /* If we are using CMAP, setup the pre-interpolation grid */
1957     if (plist[F_CMAP].ncmap() > 0)
1958     {
1959         init_cmap_grid(&sys.ffparams.cmap_grid, plist[F_CMAP].cmapAngles, plist[F_CMAP].cmakeGridSpacing);
1960         setup_cmap(plist[F_CMAP].cmakeGridSpacing, plist[F_CMAP].cmapAngles, plist[F_CMAP].cmap, &sys.ffparams.cmap_grid);
1961     }
1962
1963     set_wall_atomtype(atype, opts, ir, wi);
1964     if (bRenum)
1965     {
1966         renum_atype(plist, &sys, ir->wall_atomtype, atype, bVerbose);
1967         get_atomtype_ntypes(atype);
1968     }
1969
1970     if (ir->implicit_solvent)
1971     {
1972         gmx_fatal(FARGS, "Implicit solvation is no longer supported");
1973     }
1974
1975     /* PELA: Copy the atomtype data to the topology atomtype list */
1976     copy_atomtype_atomtypes(atype, &(sys.atomtypes));
1977
1978     if (debug)
1979     {
1980         pr_symtab(debug, 0, "After renum_atype", &sys.symtab);
1981     }
1982
1983     if (bVerbose)
1984     {
1985         fprintf(stderr, "converting bonded parameters...\n");
1986     }
1987
1988     ntype = get_atomtype_ntypes(atype);
1989     convertInteractionTypeParameters(ntype, plist, mi, intermolecular_interactions.get(),
1990                                      comb, reppow, fudgeQQ, &sys);
1991
1992     if (debug)
1993     {
1994         pr_symtab(debug, 0, "After converInteractionTypeParameters", &sys.symtab);
1995     }
1996
1997     /* set ptype to VSite for virtual sites */
1998     for (gmx_moltype_t &moltype : sys.moltype)
1999     {
2000         set_vsites_ptype(FALSE, &moltype);
2001     }
2002     if (debug)
2003     {
2004         pr_symtab(debug, 0, "After virtual sites", &sys.symtab);
2005     }
2006     /* Check velocity for virtual sites and shells */
2007     if (bGenVel)
2008     {
2009         check_vel(&sys, state.v.rvec_array());
2010     }
2011
2012     /* check for shells and inpurecs */
2013     check_shells_inputrec(&sys, ir, wi);
2014
2015     /* check masses */
2016     check_mol(&sys, wi);
2017
2018     checkForUnboundAtoms(&sys, bVerbose, wi);
2019
2020     for (const gmx_moltype_t &moltype : sys.moltype)
2021     {
2022         check_cg_sizes(ftp2fn(efTOP, NFILE, fnm), &moltype.cgs, wi);
2023     }
2024
2025     if (EI_DYNAMICS(ir->eI) && ir->eI != eiBD)
2026     {
2027         check_bonds_timestep(&sys, ir->delta_t, wi);
2028     }
2029
2030     checkDecoupledModeAccuracy(&sys, ir, wi);
2031
2032     if (EI_ENERGY_MINIMIZATION(ir->eI) && 0 == ir->nsteps)
2033     {
2034         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.");
2035     }
2036
2037     check_warning_error(wi, FARGS);
2038
2039     if (bVerbose)
2040     {
2041         fprintf(stderr, "initialising group options...\n");
2042     }
2043     do_index(mdparin, ftp2fn_null(efNDX, NFILE, fnm),
2044              &sys, bVerbose, ir,
2045              wi);
2046
2047     if (ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0)
2048     {
2049         if (EI_DYNAMICS(ir->eI) && inputrec2nboundeddim(ir) == 3)
2050         {
2051             real buffer_temp;
2052
2053             if (EI_MD(ir->eI) && ir->etc == etcNO)
2054             {
2055                 if (bGenVel)
2056                 {
2057                     buffer_temp = opts->tempi;
2058                 }
2059                 else
2060                 {
2061                     buffer_temp = calc_temp(&sys, ir, state.v.rvec_array());
2062                 }
2063                 if (buffer_temp > 0)
2064                 {
2065                     sprintf(warn_buf, "NVE simulation: will use the initial temperature of %.3f K for determining the Verlet buffer size", buffer_temp);
2066                     warning_note(wi, warn_buf);
2067                 }
2068                 else
2069                 {
2070                     sprintf(warn_buf, "NVE simulation with an initial temperature of zero: will use a Verlet buffer of %d%%. Check your energy drift!",
2071                             gmx::roundToInt(verlet_buffer_ratio_NVE_T0*100));
2072                     warning_note(wi, warn_buf);
2073                 }
2074             }
2075             else
2076             {
2077                 buffer_temp = get_max_reference_temp(ir, wi);
2078             }
2079
2080             if (EI_MD(ir->eI) && ir->etc == etcNO && buffer_temp == 0)
2081             {
2082                 /* NVE with initial T=0: we add a fixed ratio to rlist.
2083                  * Since we don't actually use verletbuf_tol, we set it to -1
2084                  * so it can't be misused later.
2085                  */
2086                 ir->rlist         *= 1.0 + verlet_buffer_ratio_NVE_T0;
2087                 ir->verletbuf_tol  = -1;
2088             }
2089             else
2090             {
2091                 /* We warn for NVE simulations with a drift tolerance that
2092                  * might result in a 1(.1)% drift over the total run-time.
2093                  * Note that we can't warn when nsteps=0, since we don't
2094                  * know how many steps the user intends to run.
2095                  */
2096                 if (EI_MD(ir->eI) && ir->etc == etcNO && ir->nstlist > 1 &&
2097                     ir->nsteps > 0)
2098                 {
2099                     const real driftTolerance = 0.01;
2100                     /* We use 2 DOF per atom = 2kT pot+kin energy,
2101                      * to be on the safe side with constraints.
2102                      */
2103                     const real totalEnergyDriftPerAtomPerPicosecond = 2*BOLTZ*buffer_temp/(ir->nsteps*ir->delta_t);
2104
2105                     if (ir->verletbuf_tol > 1.1*driftTolerance*totalEnergyDriftPerAtomPerPicosecond)
2106                     {
2107                         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.",
2108                                 ir->verletbuf_tol, ir->nsteps*ir->delta_t,
2109                                 gmx::roundToInt(ir->verletbuf_tol/totalEnergyDriftPerAtomPerPicosecond*100),
2110                                 gmx::roundToInt(100*driftTolerance),
2111                                 driftTolerance*totalEnergyDriftPerAtomPerPicosecond);
2112                         warning_note(wi, warn_buf);
2113                     }
2114                 }
2115
2116                 set_verlet_buffer(&sys, ir, buffer_temp, state.box, wi);
2117             }
2118         }
2119     }
2120
2121     /* Init the temperature coupling state */
2122     init_gtc_state(&state, ir->opts.ngtc, 0, ir->opts.nhchainlength); /* need to add nnhpres here? */
2123
2124     if (bVerbose)
2125     {
2126         fprintf(stderr, "Checking consistency between energy and charge groups...\n");
2127     }
2128     check_eg_vs_cg(&sys);
2129
2130     if (debug)
2131     {
2132         pr_symtab(debug, 0, "After index", &sys.symtab);
2133     }
2134
2135     triple_check(mdparin, ir, &sys, wi);
2136     close_symtab(&sys.symtab);
2137     if (debug)
2138     {
2139         pr_symtab(debug, 0, "After close", &sys.symtab);
2140     }
2141
2142     /* make exclusions between QM atoms and remove charges if needed */
2143     if (ir->bQMMM)
2144     {
2145         if (ir->QMMMscheme == eQMMMschemenormal && ir->ns_type == ensSIMPLE)
2146         {
2147             gmx_fatal(FARGS, "electrostatic embedding only works with grid neighboursearching, use ns-type=grid instead\n");
2148         }
2149         else
2150         {
2151             generate_qmexcl(&sys, ir, wi, GmxQmmmMode::GMX_QMMM_ORIGINAL);
2152         }
2153         if (ir->QMMMscheme != eQMMMschemeoniom)
2154         {
2155             std::vector<int> qmmmAtoms = qmmmAtomIndices(*ir, sys);
2156             removeQmmmAtomCharges(&sys, qmmmAtoms);
2157         }
2158     }
2159
2160     if (ir->eI == eiMimic)
2161     {
2162         generate_qmexcl(&sys, ir, wi, GmxQmmmMode::GMX_QMMM_MIMIC);
2163     }
2164
2165     if (ftp2bSet(efTRN, NFILE, fnm))
2166     {
2167         if (bVerbose)
2168         {
2169             fprintf(stderr, "getting data from old trajectory ...\n");
2170         }
2171         cont_status(ftp2fn(efTRN, NFILE, fnm), ftp2fn_null(efEDR, NFILE, fnm),
2172                     bNeedVel, bGenVel, fr_time, ir, &state, &sys, oenv);
2173     }
2174
2175     if (ir->ePBC == epbcXY && ir->nwall != 2)
2176     {
2177         clear_rvec(state.box[ZZ]);
2178     }
2179
2180     if (ir->cutoff_scheme != ecutsVERLET && ir->rlist > 0)
2181     {
2182         set_warning_line(wi, mdparin, -1);
2183         check_chargegroup_radii(&sys, ir, state.x.rvec_array(), wi);
2184     }
2185
2186     if (EEL_FULL(ir->coulombtype) || EVDW_PME(ir->vdwtype))
2187     {
2188         /* Calculate the optimal grid dimensions */
2189         matrix          scaledBox;
2190         EwaldBoxZScaler boxScaler(*ir);
2191         boxScaler.scaleBox(state.box, scaledBox);
2192
2193         if (ir->nkx > 0 && ir->nky > 0 && ir->nkz > 0)
2194         {
2195             /* Mark fourier_spacing as not used */
2196             ir->fourier_spacing = 0;
2197         }
2198         else if (ir->nkx != 0 && ir->nky != 0 && ir->nkz != 0)
2199         {
2200             set_warning_line(wi, mdparin, -1);
2201             warning_error(wi, "Some of the Fourier grid sizes are set, but all of them need to be set.");
2202         }
2203         const int minGridSize = minimalPmeGridSize(ir->pme_order);
2204         calcFftGrid(stdout, scaledBox, ir->fourier_spacing, minGridSize,
2205                     &(ir->nkx), &(ir->nky), &(ir->nkz));
2206         if (ir->nkx < minGridSize ||
2207             ir->nky < minGridSize ||
2208             ir->nkz < minGridSize)
2209         {
2210             warning_error(wi, "The PME grid size should be >= 2*(pme-order - 1); either manually increase the grid size or decrease pme-order");
2211         }
2212     }
2213
2214     /* MRS: eventually figure out better logic for initializing the fep
2215        values that makes declaring the lambda and declaring the state not
2216        potentially conflict if not handled correctly. */
2217     if (ir->efep != efepNO)
2218     {
2219         state.fep_state = ir->fepvals->init_fep_state;
2220         for (i = 0; i < efptNR; i++)
2221         {
2222             /* init_lambda trumps state definitions*/
2223             if (ir->fepvals->init_lambda >= 0)
2224             {
2225                 state.lambda[i] = ir->fepvals->init_lambda;
2226             }
2227             else
2228             {
2229                 if (ir->fepvals->all_lambda[i] == nullptr)
2230                 {
2231                     gmx_fatal(FARGS, "Values of lambda not set for a free energy calculation!");
2232                 }
2233                 else
2234                 {
2235                     state.lambda[i] = ir->fepvals->all_lambda[i][state.fep_state];
2236                 }
2237             }
2238         }
2239     }
2240
2241     pull_t *pull = nullptr;
2242
2243     if (ir->bPull)
2244     {
2245         pull = set_pull_init(ir, &sys, state.x.rvec_array(), state.box, state.lambda[efptMASS], wi);
2246     }
2247
2248     /* Modules that supply external potential for pull coordinates
2249      * should register those potentials here. finish_pull() will check
2250      * that providers have been registerd for all external potentials.
2251      */
2252
2253     if (ir->bDoAwh)
2254     {
2255         setStateDependentAwhParams(ir->awhParams, ir->pull, pull,
2256                                    state.box, ir->ePBC, &ir->opts, wi);
2257     }
2258
2259     if (ir->bPull)
2260     {
2261         finish_pull(pull);
2262     }
2263
2264     if (ir->bRot)
2265     {
2266         set_reference_positions(ir->rot, state.x.rvec_array(), state.box,
2267                                 opt2fn("-ref", NFILE, fnm), opt2bSet("-ref", NFILE, fnm),
2268                                 wi);
2269     }
2270
2271     /*  reset_multinr(sys); */
2272
2273     if (EEL_PME(ir->coulombtype))
2274     {
2275         float ratio = pme_load_estimate(&sys, ir, state.box);
2276         fprintf(stderr, "Estimate for the relative computational load of the PME mesh part: %.2f\n", ratio);
2277         /* With free energy we might need to do PME both for the A and B state
2278          * charges. This will double the cost, but the optimal performance will
2279          * then probably be at a slightly larger cut-off and grid spacing.
2280          */
2281         if ((ir->efep == efepNO && ratio > 1.0/2.0) ||
2282             (ir->efep != efepNO && ratio > 2.0/3.0))
2283         {
2284             warning_note(wi,
2285                          "The optimal PME mesh load for parallel simulations is below 0.5\n"
2286                          "and for highly parallel simulations between 0.25 and 0.33,\n"
2287                          "for higher performance, increase the cut-off and the PME grid spacing.\n");
2288             if (ir->efep != efepNO)
2289             {
2290                 warning_note(wi,
2291                              "For free energy simulations, the optimal load limit increases from 0.5 to 0.667\n");
2292             }
2293         }
2294     }
2295
2296     {
2297         char   warn_buf[STRLEN];
2298         double cio = compute_io(ir, sys.natoms, &sys.groups, F_NRE, 1);
2299         sprintf(warn_buf, "This run will generate roughly %.0f Mb of data", cio);
2300         if (cio > 2000)
2301         {
2302             set_warning_line(wi, mdparin, -1);
2303             warning_note(wi, warn_buf);
2304         }
2305         else
2306         {
2307             printf("%s\n", warn_buf);
2308         }
2309     }
2310
2311     if (bVerbose)
2312     {
2313         fprintf(stderr, "writing run input file...\n");
2314     }
2315
2316     done_warning(wi, FARGS);
2317     write_tpx_state(ftp2fn(efTPR, NFILE, fnm), ir, &state, &sys);
2318
2319     /* Output IMD group, if bIMD is TRUE */
2320     write_IMDgroup_to_file(ir->bIMD, ir, &state, &sys, NFILE, fnm);
2321
2322     sfree(opts->define);
2323     sfree(opts->include);
2324     sfree(opts);
2325     done_plist(plist);
2326     for (auto &mol : mi)
2327     {
2328         // Some of the contents of molinfo have been stolen, so
2329         // fullCleanUp can't be called.
2330         mol.partialCleanUp();
2331     }
2332     done_atomtype(atype);
2333     done_inputrec_strings();
2334     output_env_done(oenv);
2335
2336     return 0;
2337 }