src/gromacs/gmxpreprocess/grompp.cpp

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