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