src/gromacs/gmxpreprocess/grompp.cpp

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