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