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