src/gromacs/gmxpreprocess/grompp.cpp

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