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