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     t_gromppopts*                        opts;
1766     std::vector<MoleculeInformation>     mi;
1767     std::unique_ptr<MoleculeInformation> intermolecular_interactions;
1768     int                                  nvsite, comb;
1769     real                                 fudgeQQ;
1770     double                               reppow;
1771     const char*                          mdparin;
1772     bool                                 bNeedVel, bGenVel;
1773     gmx_output_env_t*                    oenv;
1774     gmx_bool                             bVerbose = FALSE;
1775     warninp*                             wi;
1776
1777     t_filenm fnm[] = { { efMDP, nullptr, nullptr, ffREAD },
1778                        { efMDP, "-po", "mdout", ffWRITE },
1779                        { efSTX, "-c", nullptr, ffREAD },
1780                        { efSTX, "-r", "restraint", ffOPTRD },
1781                        { efSTX, "-rb", "restraint", ffOPTRD },
1782                        { efNDX, nullptr, nullptr, ffOPTRD },
1783                        { efTOP, nullptr, nullptr, ffREAD },
1784                        { efTOP, "-pp", "processed", ffOPTWR },
1785                        { efTPR, "-o", nullptr, ffWRITE },
1786                        { efTRN, "-t", nullptr, ffOPTRD },
1787                        { efEDR, "-e", nullptr, ffOPTRD },
1788                        /* This group is needed by the VMD viewer as the start configuration for IMD sessions: */
1789                        { efGRO, "-imd", "imdgroup", ffOPTWR },
1790                        { efTRN, "-ref", "rotref", ffOPTRW | ffALLOW_MISSING } };
1791 #define NFILE asize(fnm)
1792
1793     /* Command line options */
1794     gmx_bool bRenum   = TRUE;
1795     gmx_bool bRmVSBds = TRUE, bZero = FALSE;
1796     int      i, maxwarn             = 0;
1797     real     fr_time = -1;
1798     t_pargs  pa[]    = {
1799         { "-v", FALSE, etBOOL, { &bVerbose }, "Be loud and noisy" },
1800         { "-time", FALSE, etREAL, { &fr_time }, "Take frame at or first after this time." },
1801         { "-rmvsbds",
1802           FALSE,
1803           etBOOL,
1804           { &bRmVSBds },
1805           "Remove constant bonded interactions with virtual sites" },
1806         { "-maxwarn",
1807           FALSE,
1808           etINT,
1809           { &maxwarn },
1810           "Number of allowed warnings during input processing. Not for normal use and may "
1811           "generate unstable systems" },
1812         { "-zero",
1813           FALSE,
1814           etBOOL,
1815           { &bZero },
1816           "Set parameters for bonded interactions without defaults to zero instead of "
1817           "generating an error" },
1818         { "-renum",
1819           FALSE,
1820           etBOOL,
1821           { &bRenum },
1822           "Renumber atomtypes and minimize number of atomtypes" }
1823     };
1824
1825     /* Parse the command line */
1826     if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
1827     {
1828         return 0;
1829     }
1830
1831     /* Initiate some variables */
1832     gmx::MDModules mdModules;
1833     t_inputrec     irInstance;
1834     t_inputrec*    ir = &irInstance;
1835     snew(opts, 1);
1836     snew(opts->include, STRLEN);
1837     snew(opts->define, STRLEN);
1838     snew(opts->mtsLevel2Forces, STRLEN);
1839
1840     gmx::LoggerBuilder builder;
1841     builder.addTargetStream(gmx::MDLogger::LogLevel::Info, &gmx::TextOutputFile::standardOutput());
1842     builder.addTargetStream(gmx::MDLogger::LogLevel::Warning, &gmx::TextOutputFile::standardError());
1843     gmx::LoggerOwner    logOwner(builder.build());
1844     const gmx::MDLogger logger(logOwner.logger());
1845
1846
1847     wi = init_warning(TRUE, maxwarn);
1848
1849     /* PARAMETER file processing */
1850     mdparin = opt2fn("-f", NFILE, fnm);
1851     set_warning_line(wi, mdparin, -1);
1852     try
1853     {
1854         get_ir(mdparin, opt2fn("-po", NFILE, fnm), &mdModules, ir, opts, WriteMdpHeader::yes, wi);
1855     }
1856     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
1857
1858     // Now that the MdModules have their options assigned from get_ir, subscribe
1859     // to eventual notifications during pre-processing their data
1860     mdModules.subscribeToPreProcessingNotifications();
1861
1862
1863     if (bVerbose)
1864     {
1865         GMX_LOG(logger.info)
1866                 .asParagraph()
1867                 .appendTextFormatted("checking input for internal consistency...");
1868     }
1869     check_ir(mdparin, mdModules.notifier(), ir, opts, wi);
1870
1871     if (ir->ld_seed == -1)
1872     {
1873         ir->ld_seed = static_cast<int>(gmx::makeRandomSeed());
1874         GMX_LOG(logger.info)
1875                 .asParagraph()
1876                 .appendTextFormatted("Setting the LD random seed to %" PRId64 "", ir->ld_seed);
1877     }
1878
1879     if (ir->expandedvals->lmc_seed == -1)
1880     {
1881         ir->expandedvals->lmc_seed = static_cast<int>(gmx::makeRandomSeed());
1882         GMX_LOG(logger.info)
1883                 .asParagraph()
1884                 .appendTextFormatted("Setting the lambda MC random seed to %d", ir->expandedvals->lmc_seed);
1885     }
1886
1887     bNeedVel = EI_STATE_VELOCITY(ir->eI);
1888     bGenVel  = (bNeedVel && opts->bGenVel);
1889     if (bGenVel && ir->bContinuation)
1890     {
1891         std::string warningMessage = gmx::formatString(
1892                 "Generating velocities is inconsistent with attempting "
1893                 "to continue a previous run. Choose only one of "
1894                 "gen-vel = yes and continuation = yes.");
1895         warning_error(wi, warningMessage);
1896     }
1897
1898     std::array<InteractionsOfType, F_NRE> interactions;
1899     gmx_mtop_t                            sys;
1900     PreprocessingAtomTypes                atypes;
1901     if (debug)
1902     {
1903         pr_symtab(debug, 0, "Just opened", &sys.symtab);
1904     }
1905
1906     const char* fn = ftp2fn(efTOP, NFILE, fnm);
1907     if (!gmx_fexist(fn))
1908     {
1909         gmx_fatal(FARGS, "%s does not exist", fn);
1910     }
1911
1912     t_state state;
1913     new_status(fn, opt2fn_null("-pp", NFILE, fnm), opt2fn("-c", NFILE, fnm), opts, ir, bZero,
1914                bGenVel, bVerbose, &state, &atypes, &sys, &mi, &intermolecular_interactions,
1915                interactions, &comb, &reppow, &fudgeQQ, opts->bMorse, wi, logger);
1916
1917     if (debug)
1918     {
1919         pr_symtab(debug, 0, "After new_status", &sys.symtab);
1920     }
1921
1922     nvsite = 0;
1923     /* set parameters for virtual site construction (not for vsiten) */
1924     for (size_t mt = 0; mt < sys.moltype.size(); mt++)
1925     {
1926         nvsite += set_vsites(bVerbose, &sys.moltype[mt].atoms, &atypes, mi[mt].interactions, logger);
1927     }
1928     /* now throw away all obsolete bonds, angles and dihedrals: */
1929     /* note: constraints are ALWAYS removed */
1930     if (nvsite)
1931     {
1932         for (size_t mt = 0; mt < sys.moltype.size(); mt++)
1933         {
1934             clean_vsite_bondeds(mi[mt].interactions, sys.moltype[mt].atoms.nr, bRmVSBds, logger);
1935         }
1936     }
1937
1938     if ((count_constraints(&sys, mi, wi) != 0) && (ir->eConstrAlg == econtSHAKE))
1939     {
1940         if (ir->eI == eiCG || ir->eI == eiLBFGS)
1941         {
1942             std::string warningMessage =
1943                     gmx::formatString("Can not do %s with %s, use %s", EI(ir->eI),
1944                                       econstr_names[econtSHAKE], econstr_names[econtLINCS]);
1945             warning_error(wi, warningMessage);
1946         }
1947         if (ir->bPeriodicMols)
1948         {
1949             std::string warningMessage =
1950                     gmx::formatString("Can not do periodic molecules with %s, use %s",
1951                                       econstr_names[econtSHAKE], econstr_names[econtLINCS]);
1952             warning_error(wi, warningMessage);
1953         }
1954     }
1955
1956     if (EI_SD(ir->eI) && ir->etc != etcNO)
1957     {
1958         warning_note(wi, "Temperature coupling is ignored with SD integrators.");
1959     }
1960
1961     /* Check for errors in the input now, since they might cause problems
1962      * during processing further down.
1963      */
1964     check_warning_error(wi, FARGS);
1965
1966     if (nint_ftype(&sys, mi, F_POSRES) > 0 || nint_ftype(&sys, mi, F_FBPOSRES) > 0)
1967     {
1968         if (ir->epc == epcPARRINELLORAHMAN || ir->epc == epcMTTK)
1969         {
1970             std::string warningMessage = gmx::formatString(
1971                     "You are combining position restraints with %s pressure coupling, which can "
1972                     "lead to instabilities. If you really want to combine position restraints with "
1973                     "pressure coupling, we suggest to use %s pressure coupling instead.",
1974                     EPCOUPLTYPE(ir->epc), EPCOUPLTYPE(epcBERENDSEN));
1975             warning_note(wi, warningMessage);
1976         }
1977
1978         const char* fn = opt2fn("-r", NFILE, fnm);
1979         const char* fnB;
1980
1981         if (!gmx_fexist(fn))
1982         {
1983             gmx_fatal(FARGS,
1984                       "Cannot find position restraint file %s (option -r).\n"
1985                       "From GROMACS-2018, you need to specify the position restraint "
1986                       "coordinate files explicitly to avoid mistakes, although you can "
1987                       "still use the same file as you specify for the -c option.",
1988                       fn);
1989         }
1990
1991         if (opt2bSet("-rb", NFILE, fnm))
1992         {
1993             fnB = opt2fn("-rb", NFILE, fnm);
1994             if (!gmx_fexist(fnB))
1995             {
1996                 gmx_fatal(FARGS,
1997                           "Cannot find B-state position restraint file %s (option -rb).\n"
1998                           "From GROMACS-2018, you need to specify the position restraint "
1999                           "coordinate files explicitly to avoid mistakes, although you can "
2000                           "still use the same file as you specify for the -c option.",
2001                           fn);
2002             }
2003         }
2004         else
2005         {
2006             fnB = fn;
2007         }
2008
2009         if (bVerbose)
2010         {
2011             std::string message = gmx::formatString("Reading position restraint coords from %s", fn);
2012             if (strcmp(fn, fnB) != 0)
2013             {
2014                 message += gmx::formatString(" and %s", fnB);
2015             }
2016             GMX_LOG(logger.info).asParagraph().appendText(message);
2017         }
2018         gen_posres(&sys, mi, fn, fnB, ir->refcoord_scaling, ir->pbcType, ir->posres_com,
2019                    ir->posres_comB, wi, logger);
2020     }
2021
2022     /* If we are using CMAP, setup the pre-interpolation grid */
2023     if (interactions[F_CMAP].ncmap() > 0)
2024     {
2025         init_cmap_grid(&sys.ffparams.cmap_grid, interactions[F_CMAP].cmapAngles,
2026                        interactions[F_CMAP].cmakeGridSpacing);
2027         setup_cmap(interactions[F_CMAP].cmakeGridSpacing, interactions[F_CMAP].cmapAngles,
2028                    interactions[F_CMAP].cmap, &sys.ffparams.cmap_grid);
2029     }
2030
2031     set_wall_atomtype(&atypes, opts, ir, wi, logger);
2032     if (bRenum)
2033     {
2034         atypes.renumberTypes(interactions, &sys, ir->wall_atomtype, bVerbose);
2035     }
2036
2037     if (ir->implicit_solvent)
2038     {
2039         gmx_fatal(FARGS, "Implicit solvation is no longer supported");
2040     }
2041
2042     /* PELA: Copy the atomtype data to the topology atomtype list */
2043     atypes.copyTot_atomtypes(&(sys.atomtypes));
2044
2045     if (debug)
2046     {
2047         pr_symtab(debug, 0, "After atype.renumberTypes", &sys.symtab);
2048     }
2049
2050     if (bVerbose)
2051     {
2052         GMX_LOG(logger.info).asParagraph().appendTextFormatted("converting bonded parameters...");
2053     }
2054
2055     const int ntype = atypes.size();
2056     convertInteractionsOfType(ntype, interactions, mi, intermolecular_interactions.get(), comb,
2057                               reppow, fudgeQQ, &sys);
2058
2059     if (debug)
2060     {
2061         pr_symtab(debug, 0, "After converInteractionsOfType", &sys.symtab);
2062     }
2063
2064     /* set ptype to VSite for virtual sites */
2065     for (gmx_moltype_t& moltype : sys.moltype)
2066     {
2067         set_vsites_ptype(FALSE, &moltype, logger);
2068     }
2069     if (debug)
2070     {
2071         pr_symtab(debug, 0, "After virtual sites", &sys.symtab);
2072     }
2073     /* Check velocity for virtual sites and shells */
2074     if (bGenVel)
2075     {
2076         check_vel(&sys, state.v.rvec_array());
2077     }
2078
2079     /* check for shells and inpurecs */
2080     check_shells_inputrec(&sys, ir, wi);
2081
2082     /* check masses */
2083     check_mol(&sys, wi);
2084
2085     checkForUnboundAtoms(&sys, bVerbose, wi, logger);
2086
2087     if (EI_DYNAMICS(ir->eI) && ir->eI != eiBD)
2088     {
2089         check_bonds_timestep(&sys, ir->delta_t, wi);
2090     }
2091
2092     checkDecoupledModeAccuracy(&sys, ir, wi);
2093
2094     if (EI_ENERGY_MINIMIZATION(ir->eI) && 0 == ir->nsteps)
2095     {
2096         warning_note(
2097                 wi,
2098                 "Zero-step energy minimization will alter the coordinates before calculating the "
2099                 "energy. If you just want the energy of a single point, try zero-step MD (with "
2100                 "unconstrained_start = yes). To do multiple single-point energy evaluations of "
2101                 "different configurations of the same topology, use mdrun -rerun.");
2102     }
2103
2104     check_warning_error(wi, FARGS);
2105
2106     if (bVerbose)
2107     {
2108         GMX_LOG(logger.info).asParagraph().appendTextFormatted("initialising group options...");
2109     }
2110     do_index(mdparin, ftp2fn_null(efNDX, NFILE, fnm), &sys, bVerbose, mdModules.notifier(), ir, wi);
2111
2112     if (ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0)
2113     {
2114         if (EI_DYNAMICS(ir->eI) && inputrec2nboundeddim(ir) == 3)
2115         {
2116             real buffer_temp;
2117
2118             if (EI_MD(ir->eI) && ir->etc == etcNO)
2119             {
2120                 if (bGenVel)
2121                 {
2122                     buffer_temp = opts->tempi;
2123                 }
2124                 else
2125                 {
2126                     buffer_temp = calc_temp(&sys, ir, state.v.rvec_array());
2127                 }
2128                 if (buffer_temp > 0)
2129                 {
2130                     std::string warningMessage = gmx::formatString(
2131                             "NVE simulation: will use the initial temperature of %.3f K for "
2132                             "determining the Verlet buffer size",
2133                             buffer_temp);
2134                     warning_note(wi, warningMessage);
2135                 }
2136                 else
2137                 {
2138                     std::string warningMessage = gmx::formatString(
2139                             "NVE simulation with an initial temperature of zero: will use a Verlet "
2140                             "buffer of %d%%. Check your energy drift!",
2141                             gmx::roundToInt(verlet_buffer_ratio_NVE_T0 * 100));
2142                     warning_note(wi, warningMessage);
2143                 }
2144             }
2145             else
2146             {
2147                 buffer_temp = get_max_reference_temp(ir, wi);
2148             }
2149
2150             if (EI_MD(ir->eI) && ir->etc == etcNO && buffer_temp == 0)
2151             {
2152                 /* NVE with initial T=0: we add a fixed ratio to rlist.
2153                  * Since we don't actually use verletbuf_tol, we set it to -1
2154                  * so it can't be misused later.
2155                  */
2156                 ir->rlist *= 1.0 + verlet_buffer_ratio_NVE_T0;
2157                 ir->verletbuf_tol = -1;
2158             }
2159             else
2160             {
2161                 /* We warn for NVE simulations with a drift tolerance that
2162                  * might result in a 1(.1)% drift over the total run-time.
2163                  * Note that we can't warn when nsteps=0, since we don't
2164                  * know how many steps the user intends to run.
2165                  */
2166                 if (EI_MD(ir->eI) && ir->etc == etcNO && ir->nstlist > 1 && ir->nsteps > 0)
2167                 {
2168                     const real driftTolerance = 0.01;
2169                     /* We use 2 DOF per atom = 2kT pot+kin energy,
2170                      * to be on the safe side with constraints.
2171                      */
2172                     const real totalEnergyDriftPerAtomPerPicosecond =
2173                             2 * BOLTZ * buffer_temp / (ir->nsteps * ir->delta_t);
2174
2175                     if (ir->verletbuf_tol > 1.1 * driftTolerance * totalEnergyDriftPerAtomPerPicosecond)
2176                     {
2177                         std::string warningMessage = gmx::formatString(
2178                                 "You are using a Verlet buffer tolerance of %g kJ/mol/ps for an "
2179                                 "NVE simulation of length %g ps, which can give a final drift of "
2180                                 "%d%%. For conserving energy to %d%% when using constraints, you "
2181                                 "might need to set verlet-buffer-tolerance to %.1e.",
2182                                 ir->verletbuf_tol, ir->nsteps * ir->delta_t,
2183                                 gmx::roundToInt(ir->verletbuf_tol / totalEnergyDriftPerAtomPerPicosecond * 100),
2184                                 gmx::roundToInt(100 * driftTolerance),
2185                                 driftTolerance * totalEnergyDriftPerAtomPerPicosecond);
2186                         warning_note(wi, warningMessage);
2187                     }
2188                 }
2189
2190                 set_verlet_buffer(&sys, ir, buffer_temp, state.box, wi, logger);
2191             }
2192         }
2193     }
2194
2195     /* Init the temperature coupling state */
2196     init_gtc_state(&state, ir->opts.ngtc, 0, ir->opts.nhchainlength); /* need to add nnhpres here? */
2197
2198     if (debug)
2199     {
2200         pr_symtab(debug, 0, "After index", &sys.symtab);
2201     }
2202
2203     triple_check(mdparin, ir, &sys, wi);
2204     close_symtab(&sys.symtab);
2205     if (debug)
2206     {
2207         pr_symtab(debug, 0, "After close", &sys.symtab);
2208     }
2209
2210     if (ir->eI == eiMimic)
2211     {
2212         generate_qmexcl(&sys, ir, logger);
2213     }
2214
2215     if (ftp2bSet(efTRN, NFILE, fnm))
2216     {
2217         if (bVerbose)
2218         {
2219             GMX_LOG(logger.info)
2220                     .asParagraph()
2221                     .appendTextFormatted("getting data from old trajectory ...");
2222         }
2223         cont_status(ftp2fn(efTRN, NFILE, fnm), ftp2fn_null(efEDR, NFILE, fnm), bNeedVel, bGenVel,
2224                     fr_time, ir, &state, &sys, oenv, logger);
2225     }
2226
2227     if (ir->pbcType == PbcType::XY && ir->nwall != 2)
2228     {
2229         clear_rvec(state.box[ZZ]);
2230     }
2231
2232     if (EEL_FULL(ir->coulombtype) || EVDW_PME(ir->vdwtype))
2233     {
2234         /* Calculate the optimal grid dimensions */
2235         matrix          scaledBox;
2236         EwaldBoxZScaler boxScaler(*ir);
2237         boxScaler.scaleBox(state.box, scaledBox);
2238
2239         if (ir->nkx > 0 && ir->nky > 0 && ir->nkz > 0)
2240         {
2241             /* Mark fourier_spacing as not used */
2242             ir->fourier_spacing = 0;
2243         }
2244         else if (ir->nkx != 0 && ir->nky != 0 && ir->nkz != 0)
2245         {
2246             set_warning_line(wi, mdparin, -1);
2247             warning_error(
2248                     wi, "Some of the Fourier grid sizes are set, but all of them need to be set.");
2249         }
2250         const int minGridSize = minimalPmeGridSize(ir->pme_order);
2251         calcFftGrid(stdout, scaledBox, ir->fourier_spacing, minGridSize, &(ir->nkx), &(ir->nky),
2252                     &(ir->nkz));
2253         if (ir->nkx < minGridSize || ir->nky < minGridSize || ir->nkz < minGridSize)
2254         {
2255             warning_error(wi,
2256                           "The PME grid size should be >= 2*(pme-order - 1); either manually "
2257                           "increase the grid size or decrease pme-order");
2258         }
2259     }
2260
2261     /* MRS: eventually figure out better logic for initializing the fep
2262        values that makes declaring the lambda and declaring the state not
2263        potentially conflict if not handled correctly. */
2264     if (ir->efep != efepNO)
2265     {
2266         state.fep_state = ir->fepvals->init_fep_state;
2267         for (i = 0; i < efptNR; i++)
2268         {
2269             /* init_lambda trumps state definitions*/
2270             if (ir->fepvals->init_lambda >= 0)
2271             {
2272                 state.lambda[i] = ir->fepvals->init_lambda;
2273             }
2274             else
2275             {
2276                 if (ir->fepvals->all_lambda[i] == nullptr)
2277                 {
2278                     gmx_fatal(FARGS, "Values of lambda not set for a free energy calculation!");
2279                 }
2280                 else
2281                 {
2282                     state.lambda[i] = ir->fepvals->all_lambda[i][state.fep_state];
2283                 }
2284             }
2285         }
2286     }
2287
2288     pull_t* pull = nullptr;
2289
2290     if (ir->bPull)
2291     {
2292         pull = set_pull_init(ir, &sys, state.x.rvec_array(), state.box, state.lambda[efptMASS], wi);
2293     }
2294
2295     /* Modules that supply external potential for pull coordinates
2296      * should register those potentials here. finish_pull() will check
2297      * that providers have been registerd for all external potentials.
2298      */
2299
2300     if (ir->bDoAwh)
2301     {
2302         tensor compressibility = { { 0 } };
2303         if (ir->epc != epcNO)
2304         {
2305             copy_mat(ir->compress, compressibility);
2306         }
2307         setStateDependentAwhParams(
2308                 ir->awhParams, ir->pull, pull, state.box, ir->pbcType, compressibility, &ir->opts,
2309                 ir->efep != efepNO ? ir->fepvals->all_lambda[efptFEP][ir->fepvals->init_fep_state] : 0,
2310                 sys, wi);
2311     }
2312
2313     if (ir->bPull)
2314     {
2315         finish_pull(pull);
2316     }
2317
2318     if (ir->bRot)
2319     {
2320         set_reference_positions(ir->rot, state.x.rvec_array(), state.box,
2321                                 opt2fn("-ref", NFILE, fnm), opt2bSet("-ref", NFILE, fnm), wi);
2322     }
2323
2324     /*  reset_multinr(sys); */
2325
2326     if (EEL_PME(ir->coulombtype))
2327     {
2328         float ratio = pme_load_estimate(sys, *ir, state.box);
2329         GMX_LOG(logger.info)
2330                 .asParagraph()
2331                 .appendTextFormatted(
2332                         "Estimate for the relative computational load of the PME mesh part: %.2f", ratio);
2333         /* With free energy we might need to do PME both for the A and B state
2334          * charges. This will double the cost, but the optimal performance will
2335          * then probably be at a slightly larger cut-off and grid spacing.
2336          */
2337         if ((ir->efep == efepNO && ratio > 1.0 / 2.0) || (ir->efep != efepNO && ratio > 2.0 / 3.0))
2338         {
2339             warning_note(
2340                     wi,
2341                     "The optimal PME mesh load for parallel simulations is below 0.5\n"
2342                     "and for highly parallel simulations between 0.25 and 0.33,\n"
2343                     "for higher performance, increase the cut-off and the PME grid spacing.\n");
2344             if (ir->efep != efepNO)
2345             {
2346                 warning_note(wi,
2347                              "For free energy simulations, the optimal load limit increases from "
2348                              "0.5 to 0.667\n");
2349             }
2350         }
2351     }
2352
2353     {
2354         double      cio = compute_io(ir, sys.natoms, sys.groups, F_NRE, 1);
2355         std::string warningMessage =
2356                 gmx::formatString("This run will generate roughly %.0f Mb of data", cio);
2357         if (cio > 2000)
2358         {
2359             set_warning_line(wi, mdparin, -1);
2360             warning_note(wi, warningMessage);
2361         }
2362         else
2363         {
2364             GMX_LOG(logger.info).asParagraph().appendText(warningMessage);
2365         }
2366     }
2367
2368     // Add the md modules internal parameters that are not mdp options
2369     // e.g., atom indices
2370
2371     {
2372         gmx::KeyValueTreeBuilder internalParameterBuilder;
2373         mdModules.notifier().preProcessingNotifications_.notify(internalParameterBuilder.rootObject());
2374         ir->internalParameters =
2375                 std::make_unique<gmx::KeyValueTreeObject>(internalParameterBuilder.build());
2376     }
2377
2378     if (bVerbose)
2379     {
2380         GMX_LOG(logger.info).asParagraph().appendTextFormatted("writing run input file...");
2381     }
2382
2383     done_warning(wi, FARGS);
2384     write_tpx_state(ftp2fn(efTPR, NFILE, fnm), ir, &state, &sys);
2385
2386     /* Output IMD group, if bIMD is TRUE */
2387     gmx::write_IMDgroup_to_file(ir->bIMD, ir, &state, &sys, NFILE, fnm);
2388
2389     sfree(opts->define);
2390     sfree(opts->wall_atomtype[0]);
2391     sfree(opts->wall_atomtype[1]);
2392     sfree(opts->include);
2393     sfree(opts->mtsLevel2Forces);
2394     sfree(opts);
2395     for (auto& mol : mi)
2396     {
2397         // Some of the contents of molinfo have been stolen, so
2398         // fullCleanUp can't be called.
2399         mol.partialCleanUp();
2400     }
2401     done_inputrec_strings();
2402     output_env_done(oenv);
2403
2404     return 0;
2405 }