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