2 * This file is part of the GROMACS molecular simulation package.
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,2021, 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.
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.
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.
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.
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.
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.
52 #include <sys/types.h>
54 #include "gromacs/applied_forces/awh/read_params.h"
55 #include "gromacs/commandline/pargs.h"
56 #include "gromacs/ewald/ewald_utils.h"
57 #include "gromacs/ewald/pme.h"
58 #include "gromacs/fft/calcgrid.h"
59 #include "gromacs/fileio/confio.h"
60 #include "gromacs/fileio/enxio.h"
61 #include "gromacs/fileio/tpxio.h"
62 #include "gromacs/fileio/trxio.h"
63 #include "gromacs/fileio/warninp.h"
64 #include "gromacs/gmxpreprocess/add_par.h"
65 #include "gromacs/gmxpreprocess/convparm.h"
66 #include "gromacs/gmxpreprocess/gen_maxwell_velocities.h"
67 #include "gromacs/gmxpreprocess/gpp_atomtype.h"
68 #include "gromacs/gmxpreprocess/grompp_impl.h"
69 #include "gromacs/gmxpreprocess/notset.h"
70 #include "gromacs/gmxpreprocess/readir.h"
71 #include "gromacs/gmxpreprocess/tomorse.h"
72 #include "gromacs/gmxpreprocess/topio.h"
73 #include "gromacs/gmxpreprocess/toputil.h"
74 #include "gromacs/gmxpreprocess/vsite_parm.h"
75 #include "gromacs/imd/imd.h"
76 #include "gromacs/math/functions.h"
77 #include "gromacs/math/invertmatrix.h"
78 #include "gromacs/math/units.h"
79 #include "gromacs/math/vec.h"
80 #include "gromacs/mdlib/calc_verletbuf.h"
81 #include "gromacs/mdlib/compute_io.h"
82 #include "gromacs/mdlib/constr.h"
83 #include "gromacs/mdlib/md_support.h"
84 #include "gromacs/mdlib/perf_est.h"
85 #include "gromacs/mdlib/vsite.h"
86 #include "gromacs/mdrun/mdmodules.h"
87 #include "gromacs/mdtypes/inputrec.h"
88 #include "gromacs/mdtypes/md_enums.h"
89 #include "gromacs/mdtypes/nblist.h"
90 #include "gromacs/mdtypes/state.h"
91 #include "gromacs/pbcutil/boxutilities.h"
92 #include "gromacs/pbcutil/pbc.h"
93 #include "gromacs/pulling/pull.h"
94 #include "gromacs/random/seed.h"
95 #include "gromacs/topology/ifunc.h"
96 #include "gromacs/topology/mtop_util.h"
97 #include "gromacs/topology/symtab.h"
98 #include "gromacs/topology/topology.h"
99 #include "gromacs/trajectory/trajectoryframe.h"
100 #include "gromacs/utility/arraysize.h"
101 #include "gromacs/utility/cstringutil.h"
102 #include "gromacs/utility/exceptions.h"
103 #include "gromacs/utility/fatalerror.h"
104 #include "gromacs/utility/filestream.h"
105 #include "gromacs/utility/futil.h"
106 #include "gromacs/utility/gmxassert.h"
107 #include "gromacs/utility/keyvaluetreebuilder.h"
108 #include "gromacs/utility/listoflists.h"
109 #include "gromacs/utility/logger.h"
110 #include "gromacs/utility/loggerbuilder.h"
111 #include "gromacs/utility/mdmodulesnotifiers.h"
112 #include "gromacs/utility/smalloc.h"
113 #include "gromacs/utility/snprintf.h"
115 /* TODO The implementation details should move to their own source file. */
116 InteractionOfType::InteractionOfType(gmx::ArrayRef<const int> atoms,
117 gmx::ArrayRef<const real> params,
118 const std::string& name) :
119 atoms_(atoms.begin(), atoms.end()),
120 interactionTypeName_(name)
123 params.size() <= forceParam_.size(),
124 gmx::formatString("Cannot have more parameters than the maximum number possible (%d)", MAXFORCEPARAM)
126 std::array<real, MAXFORCEPARAM>::iterator forceParamIt = forceParam_.begin();
127 for (const auto param : params)
129 *forceParamIt++ = param;
131 while (forceParamIt != forceParam_.end())
133 *forceParamIt++ = NOTSET;
137 const int& InteractionOfType::ai() const
139 GMX_RELEASE_ASSERT(!atoms_.empty(), "Need to have at least one atom set");
143 const int& InteractionOfType::aj() const
145 GMX_RELEASE_ASSERT(atoms_.size() > 1, "Need to have at least two atoms set");
149 const int& InteractionOfType::ak() const
151 GMX_RELEASE_ASSERT(atoms_.size() > 2, "Need to have at least three atoms set");
155 const int& InteractionOfType::al() const
157 GMX_RELEASE_ASSERT(atoms_.size() > 3, "Need to have at least four atoms set");
161 const int& InteractionOfType::am() const
163 GMX_RELEASE_ASSERT(atoms_.size() > 4, "Need to have at least five atoms set");
167 const real& InteractionOfType::c0() const
169 return forceParam_[0];
172 const real& InteractionOfType::c1() const
174 return forceParam_[1];
177 const real& InteractionOfType::c2() const
179 return forceParam_[2];
182 const std::string& InteractionOfType::interactionTypeName() const
184 return interactionTypeName_;
187 void InteractionOfType::sortBondAtomIds()
191 std::swap(atoms_[0], atoms_[1]);
195 void InteractionOfType::sortAngleAtomIds()
199 std::swap(atoms_[0], atoms_[2]);
203 void InteractionOfType::sortDihedralAtomIds()
207 std::swap(atoms_[0], atoms_[3]);
208 std::swap(atoms_[1], atoms_[2]);
212 void InteractionOfType::sortAtomIds()
222 else if (isDihedral())
224 sortDihedralAtomIds();
228 GMX_THROW(gmx::InternalError(
229 "Can not sort parameters other than bonds, angles or dihedrals"));
233 void InteractionOfType::setForceParameter(int pos, real value)
235 GMX_RELEASE_ASSERT(pos < MAXFORCEPARAM,
236 "Can't set parameter beyond the maximum number of parameters");
237 forceParam_[pos] = value;
240 void MoleculeInformation::initMolInfo()
244 init_t_atoms(&atoms, 0, FALSE);
247 void MoleculeInformation::partialCleanUp()
252 void MoleculeInformation::fullCleanUp()
258 static int rm_interactions(int ifunc, gmx::ArrayRef<MoleculeInformation> mols)
261 /* For all the molecule types */
262 for (auto& mol : mols)
264 n += mol.interactions[ifunc].size();
265 mol.interactions[ifunc].interactionTypes.clear();
270 static int check_atom_names(const char* fn1,
274 const gmx::MDLogger& logger)
276 int m, i, j, nmismatch;
279 constexpr int c_maxNumberOfMismatches = 20;
281 if (mtop->natoms != at->nr)
283 gmx_incons("comparing atom names");
288 for (const gmx_molblock_t& molb : mtop->molblock)
290 tat = &mtop->moltype[molb.type].atoms;
291 for (m = 0; m < molb.nmol; m++)
293 for (j = 0; j < tat->nr; j++)
295 if (strcmp(*(tat->atomname[j]), *(at->atomname[i])) != 0)
297 if (nmismatch < c_maxNumberOfMismatches)
299 GMX_LOG(logger.warning)
301 .appendTextFormatted(
302 "atom name %d in %s and %s does not match (%s - %s)",
309 else if (nmismatch == c_maxNumberOfMismatches)
311 GMX_LOG(logger.warning)
313 .appendTextFormatted("(more than %d non-matching atom names)",
314 c_maxNumberOfMismatches);
326 static void check_bonds_timestep(const gmx_mtop_t* mtop, double dt, warninp* wi)
328 /* This check is not intended to ensure accurate integration,
329 * rather it is to signal mistakes in the mdp settings.
330 * A common mistake is to forget to turn on constraints
331 * for MD after energy minimization with flexible bonds.
332 * This check can also detect too large time steps for flexible water
333 * models, but such errors will often be masked by the constraints
334 * mdp options, which turns flexible water into water with bond constraints,
335 * but without an angle constraint. Unfortunately such incorrect use
336 * of water models can not easily be detected without checking
337 * for specific model names.
339 * The stability limit of leap-frog or velocity verlet is 4.44 steps
340 * per oscillational period.
341 * But accurate bonds distributions are lost far before that limit.
342 * To allow relatively common schemes (although not common with Gromacs)
343 * of dt=1 fs without constraints and dt=2 fs with only H-bond constraints
344 * we set the note limit to 10.
346 int min_steps_warn = 5;
347 int min_steps_note = 10;
349 int i, a1, a2, w_a1, w_a2, j;
350 real twopi2, limit2, fc, re, m1, m2, period2, w_period2;
351 bool bFound, bWater, bWarn;
353 /* Get the interaction parameters */
354 gmx::ArrayRef<const t_iparams> ip = mtop->ffparams.iparams;
356 twopi2 = gmx::square(2 * M_PI);
358 limit2 = gmx::square(min_steps_note * dt);
363 const gmx_moltype_t* w_moltype = nullptr;
364 for (const gmx_moltype_t& moltype : mtop->moltype)
366 const t_atom* atom = moltype.atoms.atom;
367 const InteractionLists& ilist = moltype.ilist;
368 const InteractionList& ilc = ilist[F_CONSTR];
369 const InteractionList& ils = ilist[F_SETTLE];
370 for (ftype = 0; ftype < F_NRE; ftype++)
372 if (!(ftype == F_BONDS || ftype == F_G96BONDS || ftype == F_HARMONIC))
377 const InteractionList& ilb = ilist[ftype];
378 for (i = 0; i < ilb.size(); i += 3)
380 fc = ip[ilb.iatoms[i]].harmonic.krA;
381 re = ip[ilb.iatoms[i]].harmonic.rA;
382 if (ftype == F_G96BONDS)
384 /* Convert squared sqaure fc to harmonic fc */
387 a1 = ilb.iatoms[i + 1];
388 a2 = ilb.iatoms[i + 2];
391 if (fc > 0 && m1 > 0 && m2 > 0)
393 period2 = twopi2 * m1 * m2 / ((m1 + m2) * fc);
397 period2 = GMX_FLOAT_MAX;
401 fprintf(debug, "fc %g m1 %g m2 %g period %g\n", fc, m1, m2, std::sqrt(period2));
403 if (period2 < limit2)
406 for (j = 0; j < ilc.size(); j += 3)
408 if ((ilc.iatoms[j + 1] == a1 && ilc.iatoms[j + 2] == a2)
409 || (ilc.iatoms[j + 1] == a2 && ilc.iatoms[j + 2] == a1))
414 for (j = 0; j < ils.size(); j += 4)
416 if ((a1 == ils.iatoms[j + 1] || a1 == ils.iatoms[j + 2] || a1 == ils.iatoms[j + 3])
417 && (a2 == ils.iatoms[j + 1] || a2 == ils.iatoms[j + 2]
418 || a2 == ils.iatoms[j + 3]))
423 if (!bFound && (w_moltype == nullptr || period2 < w_period2))
425 w_moltype = &moltype;
435 if (w_moltype != nullptr)
437 bWarn = (w_period2 < gmx::square(min_steps_warn * dt));
438 /* A check that would recognize most water models */
439 bWater = ((*w_moltype->atoms.atomname[0])[0] == 'O' && w_moltype->atoms.nr <= 5);
440 std::string warningMessage = gmx::formatString(
441 "The bond in molecule-type %s between atoms %d %s and %d %s has an estimated "
442 "oscillational period of %.1e ps, which is less than %d times the time step of "
447 *w_moltype->atoms.atomname[w_a1],
449 *w_moltype->atoms.atomname[w_a2],
450 std::sqrt(w_period2),
451 bWarn ? min_steps_warn : min_steps_note,
453 bWater ? "Maybe you asked for fexible water."
454 : "Maybe you forgot to change the constraints mdp option.");
457 warning(wi, warningMessage.c_str());
461 warning_note(wi, warningMessage.c_str());
466 static void check_vel(gmx_mtop_t* mtop, rvec v[])
468 for (const AtomProxy atomP : AtomRange(*mtop))
470 const t_atom& local = atomP.atom();
471 int i = atomP.globalAtomNumber();
472 if (local.ptype == ParticleType::Shell || local.ptype == ParticleType::Bond
473 || local.ptype == ParticleType::VSite)
480 static void check_shells_inputrec(gmx_mtop_t* mtop, t_inputrec* ir, warninp* wi)
484 for (const AtomProxy atomP : AtomRange(*mtop))
486 const t_atom& local = atomP.atom();
487 if (local.ptype == ParticleType::Shell || local.ptype == ParticleType::Bond)
492 if ((nshells > 0) && (ir->nstcalcenergy != 1))
494 set_warning_line(wi, "unknown", -1);
495 std::string warningMessage = gmx::formatString(
496 "There are %d shells, changing nstcalcenergy from %d to 1", nshells, ir->nstcalcenergy);
497 ir->nstcalcenergy = 1;
498 warning(wi, warningMessage.c_str());
502 /* TODO Decide whether this function can be consolidated with
503 * gmx_mtop_ftype_count */
504 static int nint_ftype(gmx_mtop_t* mtop, gmx::ArrayRef<const MoleculeInformation> mi, int ftype)
507 for (const gmx_molblock_t& molb : mtop->molblock)
509 nint += molb.nmol * mi[molb.type].interactions[ftype].size();
515 /* This routine reorders the molecule type array
516 * in the order of use in the molblocks,
517 * unused molecule types are deleted.
519 static void renumber_moltypes(gmx_mtop_t* sys, std::vector<MoleculeInformation>* molinfo)
522 std::vector<int> order;
523 for (gmx_molblock_t& molblock : sys->molblock)
525 const auto found = std::find_if(order.begin(), order.end(), [&molblock](const auto& entry) {
526 return molblock.type == entry;
528 if (found == order.end())
530 /* This type did not occur yet, add it */
531 order.push_back(molblock.type);
532 molblock.type = order.size() - 1;
536 molblock.type = std::distance(order.begin(), found);
540 /* We still need to reorder the molinfo structs */
541 std::vector<MoleculeInformation> minew(order.size());
543 for (auto& mi : *molinfo)
545 const auto found = std::find(order.begin(), order.end(), index);
546 if (found != order.end())
548 int pos = std::distance(order.begin(), found);
553 // Need to manually clean up memory ....
562 static void molinfo2mtop(gmx::ArrayRef<const MoleculeInformation> mi, gmx_mtop_t* mtop)
564 mtop->moltype.resize(mi.size());
566 for (const auto& mol : mi)
568 gmx_moltype_t& molt = mtop->moltype[pos];
569 molt.name = mol.name;
570 molt.atoms = mol.atoms;
571 /* ilists are copied later */
572 molt.excls = mol.excls;
577 static void new_status(const char* topfile,
578 const char* topppfile,
586 PreprocessingAtomTypes* atypes,
588 std::vector<MoleculeInformation>* mi,
589 std::unique_ptr<MoleculeInformation>* intermolecular_interactions,
590 gmx::ArrayRef<InteractionsOfType> interactions,
591 CombinationRule* comb,
596 const gmx::MDLogger& logger)
598 std::vector<gmx_molblock_t> molblock;
600 bool ffParametrizedWithHBondConstraints;
602 /* TOPOLOGY processing */
603 sys->name = do_top(bVerbose,
615 intermolecular_interactions,
618 &ffParametrizedWithHBondConstraints,
622 sys->molblock.clear();
625 for (const gmx_molblock_t& molb : molblock)
627 if (!sys->molblock.empty() && molb.type == sys->molblock.back().type)
629 /* Merge consecutive blocks with the same molecule type */
630 sys->molblock.back().nmol += molb.nmol;
632 else if (molb.nmol > 0)
634 /* Add a new molblock to the topology */
635 sys->molblock.push_back(molb);
637 sys->natoms += molb.nmol * (*mi)[sys->molblock.back().type].atoms.nr;
639 if (sys->molblock.empty())
641 gmx_fatal(FARGS, "No molecules were defined in the system");
644 renumber_moltypes(sys, mi);
648 convert_harmonics(*mi, atypes);
651 if (ir->eDisre == DistanceRestraintRefinement::None)
653 i = rm_interactions(F_DISRES, *mi);
656 set_warning_line(wi, "unknown", -1);
657 std::string warningMessage =
658 gmx::formatString("disre = no, removed %d distance restraints", i);
659 warning_note(wi, warningMessage.c_str());
664 i = rm_interactions(F_ORIRES, *mi);
667 set_warning_line(wi, "unknown", -1);
668 std::string warningMessage =
669 gmx::formatString("orire = no, removed %d orientation restraints", i);
670 warning_note(wi, warningMessage.c_str());
674 /* Copy structures from msys to sys */
675 molinfo2mtop(*mi, sys);
679 /* Discourage using the, previously common, setup of applying constraints
680 * to all bonds with force fields that have been parametrized with H-bond
681 * constraints only. Do not print note with large timesteps or vsites.
683 if (opts->nshake == eshALLBONDS && ffParametrizedWithHBondConstraints && ir->delta_t < 0.0026
684 && gmx_mtop_ftype_count(*sys, F_VSITE3FD) == 0)
686 set_warning_line(wi, "unknown", -1);
688 "You are using constraints on all bonds, whereas the forcefield "
689 "has been parametrized only with constraints involving hydrogen atoms. "
690 "We suggest using constraints = h-bonds instead, this will also improve "
694 /* COORDINATE file processing */
697 GMX_LOG(logger.info).asParagraph().appendTextFormatted("processing coordinates...");
704 read_tps_conf(confin, conftop, nullptr, &x, &v, state->box, FALSE);
705 state->natoms = conftop->atoms.nr;
706 if (state->natoms != sys->natoms)
709 "number of coordinates in coordinate file (%s, %d)\n"
710 " does not match topology (%s, %d)",
716 /* It would be nice to get rid of the copies below, but we don't know
717 * a priori if the number of atoms in confin matches what we expect.
719 state->flags |= enumValueToBitMask(StateEntry::X);
722 state->flags |= enumValueToBitMask(StateEntry::V);
724 state_change_natoms(state, state->natoms);
725 std::copy(x, x + state->natoms, state->x.data());
729 std::copy(v, v + state->natoms, state->v.data());
732 /* This call fixes the box shape for runs with pressure scaling */
733 set_box_rel(ir, state);
735 nmismatch = check_atom_names(topfile, confin, sys, &conftop->atoms, logger);
741 std::string warningMessage = gmx::formatString(
742 "%d non-matching atom name%s\n"
743 "atom names from %s will be used\n"
744 "atom names from %s will be ignored\n",
746 (nmismatch == 1) ? "" : "s",
749 warning(wi, warningMessage.c_str());
752 /* Do more checks, mostly related to constraints */
757 .appendTextFormatted("double-checking input for internal consistency...");
760 bool bHasNormalConstraints =
761 0 < (nint_ftype(sys, *mi, F_CONSTR) + nint_ftype(sys, *mi, F_CONSTRNC));
762 bool bHasAnyConstraints = bHasNormalConstraints || 0 < nint_ftype(sys, *mi, F_SETTLE);
763 double_check(ir, state->box, bHasNormalConstraints, bHasAnyConstraints, wi);
770 snew(mass, state->natoms);
771 for (const AtomProxy atomP : AtomRange(*sys))
773 const t_atom& local = atomP.atom();
774 int i = atomP.globalAtomNumber();
778 if (opts->seed == -1)
780 opts->seed = static_cast<int>(gmx::makeRandomSeed());
781 GMX_LOG(logger.info).asParagraph().appendTextFormatted("Setting gen_seed to %d", opts->seed);
783 state->flags |= enumValueToBitMask(StateEntry::V);
784 maxwell_speed(opts->tempi, opts->seed, sys, state->v.rvec_array(), logger);
786 stop_cm(logger, state->natoms, mass, state->x.rvec_array(), state->v.rvec_array());
791 static void copy_state(const char* slog, t_trxframe* fr, bool bReadVel, t_state* state, double* use_time)
793 if (fr->not_ok & FRAME_NOT_OK)
795 gmx_fatal(FARGS, "Can not start from an incomplete frame");
799 gmx_fatal(FARGS, "Did not find a frame with coordinates in file %s", slog);
802 std::copy(fr->x, fr->x + state->natoms, state->x.data());
807 gmx_incons("Trajecory frame unexpectedly does not contain velocities");
809 std::copy(fr->v, fr->v + state->natoms, state->v.data());
813 copy_mat(fr->box, state->box);
816 *use_time = fr->time;
819 static void cont_status(const char* slog,
827 const gmx_output_env_t* oenv,
828 const gmx::MDLogger& logger)
829 /* If fr_time == -1 read the last frame available which is complete */
836 bReadVel = (bNeedVel && !bGenVel);
840 .appendTextFormatted("Reading Coordinates%s and Box size from old trajectory",
841 bReadVel ? ", Velocities" : "");
844 GMX_LOG(logger.info).asParagraph().appendTextFormatted("Will read whole trajectory");
848 GMX_LOG(logger.info).asParagraph().appendTextFormatted("Will read till time %g", fr_time);
856 .appendTextFormatted(
857 "Velocities generated: "
858 "ignoring velocities in input trajectory");
860 read_first_frame(oenv, &fp, slog, &fr, TRX_NEED_X);
864 read_first_frame(oenv, &fp, slog, &fr, TRX_NEED_X | TRX_NEED_V);
868 GMX_LOG(logger.warning)
870 .appendTextFormatted(
871 "WARNING: Did not find a frame with velocities in file %s,\n"
872 " all velocities will be set to zero!",
874 for (auto& vi : makeArrayRef(state->v))
879 /* Search for a frame without velocities */
881 read_first_frame(oenv, &fp, slog, &fr, TRX_NEED_X);
885 state->natoms = fr.natoms;
887 if (sys->natoms != state->natoms)
890 "Number of atoms in Topology "
891 "is not the same as in Trajectory");
893 copy_state(slog, &fr, bReadVel, state, &use_time);
895 /* Find the appropriate frame */
896 while ((fr_time == -1 || fr.time < fr_time) && read_next_frame(oenv, fp, &fr))
898 copy_state(slog, &fr, bReadVel, state, &use_time);
903 /* Set the relative box lengths for preserving the box shape.
904 * Note that this call can lead to differences in the last bit
905 * with respect to using gmx convert-tpr to create a [REF].tpx[ref] file.
907 set_box_rel(ir, state);
909 GMX_LOG(logger.info).asParagraph().appendTextFormatted("Using frame at t = %g ps", use_time);
910 GMX_LOG(logger.info).asParagraph().appendTextFormatted("Starting time for run is %g ps", ir->init_t);
912 if ((ir->epc != PressureCoupling::No || ir->etc == TemperatureCoupling::NoseHoover) && ener)
914 get_enx_state(ener, use_time, sys->groups, ir, state);
915 preserve_box_shape(ir, state->box_rel, state->boxv);
919 static void read_posres(gmx_mtop_t* mtop,
920 gmx::ArrayRef<const MoleculeInformation> molinfo,
923 RefCoordScaling rc_scaling,
927 const gmx::MDLogger& logger)
935 int natoms, npbcdim = 0;
940 read_tps_conf(fn, top, nullptr, &x, &v, box, FALSE);
941 natoms = top->atoms.nr;
944 if (natoms != mtop->natoms)
946 std::string warningMessage = gmx::formatString(
947 "The number of atoms in %s (%d) does not match the number of atoms in the topology "
948 "(%d). Will assume that the first %d atoms in the topology and %s match.",
952 std::min(mtop->natoms, natoms),
954 warning(wi, warningMessage.c_str());
957 npbcdim = numPbcDimensions(pbcType);
958 GMX_RELEASE_ASSERT(npbcdim <= DIM, "Invalid npbcdim");
960 if (rc_scaling != RefCoordScaling::No)
962 copy_mat(box, invbox);
963 for (int j = npbcdim; j < DIM; j++)
965 clear_rvec(invbox[j]);
968 gmx::invertBoxMatrix(invbox, invbox);
971 /* Copy the reference coordinates to mtop */
975 snew(hadAtom, natoms);
976 for (gmx_molblock_t& molb : mtop->molblock)
978 nat_molb = molb.nmol * mtop->moltype[molb.type].atoms.nr;
979 const InteractionsOfType* pr = &(molinfo[molb.type].interactions[F_POSRES]);
980 const InteractionsOfType* prfb = &(molinfo[molb.type].interactions[F_FBPOSRES]);
981 if (pr->size() > 0 || prfb->size() > 0)
983 atom = mtop->moltype[molb.type].atoms.atom;
984 for (const auto& restraint : pr->interactionTypes)
986 int ai = restraint.ai();
990 "Position restraint atom index (%d) in moltype '%s' is larger than "
991 "number of atoms in %s (%d).\n",
993 *molinfo[molb.type].name,
998 if (rc_scaling == RefCoordScaling::Com)
1000 /* Determine the center of mass of the posres reference coordinates */
1001 for (int j = 0; j < npbcdim; j++)
1003 sum[j] += atom[ai].m * x[a + ai][j];
1005 totmass += atom[ai].m;
1008 /* Same for flat-bottomed posres, but do not count an atom twice for COM */
1009 for (const auto& restraint : prfb->interactionTypes)
1011 int ai = restraint.ai();
1015 "Position restraint atom index (%d) in moltype '%s' is larger than "
1016 "number of atoms in %s (%d).\n",
1018 *molinfo[molb.type].name,
1022 if (rc_scaling == RefCoordScaling::Com && !hadAtom[ai])
1024 /* Determine the center of mass of the posres reference coordinates */
1025 for (int j = 0; j < npbcdim; j++)
1027 sum[j] += atom[ai].m * x[a + ai][j];
1029 totmass += atom[ai].m;
1034 molb.posres_xA.resize(nat_molb);
1035 for (int i = 0; i < nat_molb; i++)
1037 copy_rvec(x[a + i], molb.posres_xA[i]);
1042 molb.posres_xB.resize(nat_molb);
1043 for (int i = 0; i < nat_molb; i++)
1045 copy_rvec(x[a + i], molb.posres_xB[i]);
1051 if (rc_scaling == RefCoordScaling::Com)
1055 gmx_fatal(FARGS, "The total mass of the position restraint atoms is 0");
1057 for (int j = 0; j < npbcdim; j++)
1059 com[j] = sum[j] / totmass;
1061 GMX_LOG(logger.info)
1063 .appendTextFormatted(
1064 "The center of mass of the position restraint coord's is %6.3f %6.3f %6.3f",
1070 if (rc_scaling != RefCoordScaling::No)
1072 GMX_ASSERT(npbcdim <= DIM, "Only DIM dimensions can have PBC");
1074 for (gmx_molblock_t& molb : mtop->molblock)
1076 nat_molb = molb.nmol * mtop->moltype[molb.type].atoms.nr;
1077 if (!molb.posres_xA.empty() || !molb.posres_xB.empty())
1079 std::vector<gmx::RVec>& xp = (!bTopB ? molb.posres_xA : molb.posres_xB);
1080 for (int i = 0; i < nat_molb; i++)
1082 for (int j = 0; j < npbcdim; j++)
1084 if (rc_scaling == RefCoordScaling::All)
1086 /* Convert from Cartesian to crystal coordinates */
1087 xp[i][j] *= invbox[j][j];
1088 for (int k = j + 1; k < npbcdim; k++)
1090 xp[i][j] += invbox[k][j] * xp[i][k];
1093 else if (rc_scaling == RefCoordScaling::Com)
1095 /* Subtract the center of mass */
1103 if (rc_scaling == RefCoordScaling::Com)
1105 /* Convert the COM from Cartesian to crystal coordinates */
1106 for (int j = 0; j < npbcdim; j++)
1108 com[j] *= invbox[j][j];
1109 for (int k = j + 1; k < npbcdim; k++)
1111 com[j] += invbox[k][j] * com[k];
1122 static void gen_posres(gmx_mtop_t* mtop,
1123 gmx::ArrayRef<const MoleculeInformation> mi,
1126 RefCoordScaling rc_scaling,
1131 const gmx::MDLogger& logger)
1133 read_posres(mtop, mi, FALSE, fnA, rc_scaling, pbcType, com, wi, logger);
1134 /* It is safer to simply read the b-state posres rather than trying
1135 * to be smart and copy the positions.
1137 read_posres(mtop, mi, TRUE, fnB, rc_scaling, pbcType, comB, wi, logger);
1140 static void set_wall_atomtype(PreprocessingAtomTypes* at,
1144 const gmx::MDLogger& logger)
1150 GMX_LOG(logger.info).asParagraph().appendTextFormatted("Searching the wall atom type(s)");
1152 for (i = 0; i < ir->nwall; i++)
1154 auto atomType = at->atomTypeFromName(opts->wall_atomtype[i]);
1155 if (!atomType.has_value())
1157 std::string warningMessage = gmx::formatString(
1158 "Specified wall atom type %s is not defined", opts->wall_atomtype[i]);
1159 warning_error(wi, warningMessage.c_str());
1163 ir->wall_atomtype[i] = *atomType;
1168 static int nrdf_internal(const t_atoms* atoms)
1173 for (i = 0; i < atoms->nr; i++)
1175 /* Vsite ptype might not be set here yet, so also check the mass */
1176 if ((atoms->atom[i].ptype == ParticleType::Atom || atoms->atom[i].ptype == ParticleType::Nucleus)
1177 && atoms->atom[i].m > 0)
1184 case 0: // Fall through intended
1185 case 1: nrdf = 0; break;
1186 case 2: nrdf = 1; break;
1187 default: nrdf = nmass * 3 - 6; break;
1193 static void spline1d(double dx, const double* y, int n, double* u, double* y2)
1201 for (i = 1; i < n - 1; i++)
1203 p = 0.5 * y2[i - 1] + 2.0;
1205 q = (y[i + 1] - 2.0 * y[i] + y[i - 1]) / dx;
1206 u[i] = (3.0 * q / dx - 0.5 * u[i - 1]) / p;
1211 for (i = n - 2; i >= 0; i--)
1213 y2[i] = y2[i] * y2[i + 1] + u[i];
1219 interpolate1d(double xmin, double dx, const double* ya, const double* y2a, double x, double* y, double* y1)
1224 ix = static_cast<int>((x - xmin) / dx);
1226 a = (xmin + (ix + 1) * dx - x) / dx;
1227 b = (x - xmin - ix * dx) / dx;
1229 *y = a * ya[ix] + b * ya[ix + 1]
1230 + ((a * a * a - a) * y2a[ix] + (b * b * b - b) * y2a[ix + 1]) * (dx * dx) / 6.0;
1231 *y1 = (ya[ix + 1] - ya[ix]) / dx - (3.0 * a * a - 1.0) / 6.0 * dx * y2a[ix]
1232 + (3.0 * b * b - 1.0) / 6.0 * dx * y2a[ix + 1];
1236 static void setup_cmap(int grid_spacing, int nc, gmx::ArrayRef<const real> grid, gmx_cmap_t* cmap_grid)
1238 int i, j, k, ii, jj, kk, idx;
1240 double dx, xmin, v, v1, v2, v12;
1243 std::vector<double> tmp_u(2 * grid_spacing, 0.0);
1244 std::vector<double> tmp_u2(2 * grid_spacing, 0.0);
1245 std::vector<double> tmp_yy(2 * grid_spacing, 0.0);
1246 std::vector<double> tmp_y1(2 * grid_spacing, 0.0);
1247 std::vector<double> tmp_t2(2 * grid_spacing * 2 * grid_spacing, 0.0);
1248 std::vector<double> tmp_grid(2 * grid_spacing * 2 * grid_spacing, 0.0);
1250 dx = 360.0 / grid_spacing;
1251 xmin = -180.0 - dx * grid_spacing / 2;
1253 for (kk = 0; kk < nc; kk++)
1255 /* Compute an offset depending on which cmap we are using
1256 * Offset will be the map number multiplied with the
1257 * grid_spacing * grid_spacing * 2
1259 offset = kk * grid_spacing * grid_spacing * 2;
1261 for (i = 0; i < 2 * grid_spacing; i++)
1263 ii = (i + grid_spacing - grid_spacing / 2) % grid_spacing;
1265 for (j = 0; j < 2 * grid_spacing; j++)
1267 jj = (j + grid_spacing - grid_spacing / 2) % grid_spacing;
1268 tmp_grid[i * grid_spacing * 2 + j] = grid[offset + ii * grid_spacing + jj];
1272 for (i = 0; i < 2 * grid_spacing; i++)
1275 &(tmp_grid[2 * grid_spacing * i]),
1278 &(tmp_t2[2 * grid_spacing * i]));
1281 for (i = grid_spacing / 2; i < grid_spacing + grid_spacing / 2; i++)
1283 ii = i - grid_spacing / 2;
1284 phi = ii * dx - 180.0;
1286 for (j = grid_spacing / 2; j < grid_spacing + grid_spacing / 2; j++)
1288 jj = j - grid_spacing / 2;
1289 psi = jj * dx - 180.0;
1291 for (k = 0; k < 2 * grid_spacing; k++)
1295 &(tmp_grid[2 * grid_spacing * k]),
1296 &(tmp_t2[2 * grid_spacing * k]),
1302 spline1d(dx, tmp_yy.data(), 2 * grid_spacing, tmp_u.data(), tmp_u2.data());
1303 interpolate1d(xmin, dx, tmp_yy.data(), tmp_u2.data(), phi, &v, &v1);
1304 spline1d(dx, tmp_y1.data(), 2 * grid_spacing, tmp_u.data(), tmp_u2.data());
1305 interpolate1d(xmin, dx, tmp_y1.data(), tmp_u2.data(), phi, &v2, &v12);
1307 idx = ii * grid_spacing + jj;
1308 cmap_grid->cmapdata[kk].cmap[idx * 4] = grid[offset + ii * grid_spacing + jj];
1309 cmap_grid->cmapdata[kk].cmap[idx * 4 + 1] = v1;
1310 cmap_grid->cmapdata[kk].cmap[idx * 4 + 2] = v2;
1311 cmap_grid->cmapdata[kk].cmap[idx * 4 + 3] = v12;
1317 static void init_cmap_grid(gmx_cmap_t* cmap_grid, int ngrid, int grid_spacing)
1321 cmap_grid->grid_spacing = grid_spacing;
1322 nelem = cmap_grid->grid_spacing * cmap_grid->grid_spacing;
1324 cmap_grid->cmapdata.resize(ngrid);
1326 for (i = 0; i < ngrid; i++)
1328 cmap_grid->cmapdata[i].cmap.resize(4 * nelem);
1333 static int count_constraints(const gmx_mtop_t* mtop, gmx::ArrayRef<const MoleculeInformation> mi, warninp* wi)
1335 int count, count_mol;
1338 for (const gmx_molblock_t& molb : mtop->molblock)
1341 gmx::ArrayRef<const InteractionsOfType> interactions = mi[molb.type].interactions;
1343 for (int i = 0; i < F_NRE; i++)
1347 count_mol += 3 * interactions[i].size();
1349 else if (interaction_function[i].flags & IF_CONSTRAINT)
1351 count_mol += interactions[i].size();
1355 if (count_mol > nrdf_internal(&mi[molb.type].atoms))
1357 std::string warningMessage = gmx::formatString(
1358 "Molecule type '%s' has %d constraints.\n"
1359 "For stability and efficiency there should not be more constraints than "
1360 "internal number of degrees of freedom: %d.\n",
1361 *mi[molb.type].name,
1363 nrdf_internal(&mi[molb.type].atoms));
1364 warning(wi, warningMessage.c_str());
1366 count += molb.nmol * count_mol;
1372 static real calc_temp(const gmx_mtop_t* mtop, const t_inputrec* ir, rvec* v)
1375 for (const AtomProxy atomP : AtomRange(*mtop))
1377 const t_atom& local = atomP.atom();
1378 int i = atomP.globalAtomNumber();
1379 sum_mv2 += local.m * norm2(v[i]);
1383 for (int g = 0; g < ir->opts.ngtc; g++)
1385 nrdf += ir->opts.nrdf[g];
1388 return sum_mv2 / (nrdf * gmx::c_boltz);
1391 static real get_max_reference_temp(const t_inputrec* ir, warninp* wi)
1399 for (i = 0; i < ir->opts.ngtc; i++)
1401 if (ir->opts.tau_t[i] < 0)
1407 ref_t = std::max(ref_t, ir->opts.ref_t[i]);
1413 std::string warningMessage = gmx::formatString(
1414 "Some temperature coupling groups do not use temperature coupling. We will assume "
1415 "their temperature is not more than %.3f K. If their temperature is higher, the "
1416 "energy error and the Verlet buffer might be underestimated.",
1418 warning(wi, warningMessage.c_str());
1424 /* Checks if there are unbound atoms in moleculetype molt.
1425 * Prints a note for each unbound atoms and a warning if any is present.
1427 static void checkForUnboundAtoms(const gmx_moltype_t* molt, gmx_bool bVerbose, warninp* wi, const gmx::MDLogger& logger)
1429 const t_atoms* atoms = &molt->atoms;
1433 /* Only one atom, there can't be unbound atoms */
1437 std::vector<int> count(atoms->nr, 0);
1439 for (int ftype = 0; ftype < F_NRE; ftype++)
1441 if (((interaction_function[ftype].flags & IF_BOND) && NRAL(ftype) == 2 && ftype != F_CONNBONDS)
1442 || (interaction_function[ftype].flags & IF_CONSTRAINT) || ftype == F_SETTLE)
1444 const InteractionList& il = molt->ilist[ftype];
1445 const int nral = NRAL(ftype);
1447 for (int i = 0; i < il.size(); i += 1 + nral)
1449 for (int j = 0; j < nral; j++)
1451 count[il.iatoms[i + 1 + j]]++;
1457 int numDanglingAtoms = 0;
1458 for (int a = 0; a < atoms->nr; a++)
1460 if (atoms->atom[a].ptype != ParticleType::VSite && count[a] == 0)
1464 GMX_LOG(logger.warning)
1466 .appendTextFormatted(
1467 "Atom %d '%s' in moleculetype '%s' is not bound by a potential or "
1468 "constraint to any other atom in the same moleculetype.",
1470 *atoms->atomname[a],
1477 if (numDanglingAtoms > 0)
1479 std::string warningMessage = gmx::formatString(
1480 "In moleculetype '%s' %d atoms are not bound by a potential or constraint to any "
1481 "other atom in the same moleculetype. Although technically this might not cause "
1482 "issues in a simulation, this often means that the user forgot to add a "
1483 "bond/potential/constraint or put multiple molecules in the same moleculetype "
1484 "definition by mistake. Run with -v to get information for each atom.",
1487 warning_note(wi, warningMessage.c_str());
1491 /* Checks all moleculetypes for unbound atoms */
1492 static void checkForUnboundAtoms(const gmx_mtop_t* mtop, gmx_bool bVerbose, warninp* wi, const gmx::MDLogger& logger)
1494 for (const gmx_moltype_t& molt : mtop->moltype)
1496 checkForUnboundAtoms(&molt, bVerbose, wi, logger);
1500 /*! \brief Checks if there are decoupled modes in moleculetype \p molt.
1502 * The specific decoupled modes this routine check for are angle modes
1503 * where the two bonds are constrained and the atoms a both ends are only
1504 * involved in a single constraint; the mass of the two atoms needs to
1505 * differ by more than \p massFactorThreshold.
1507 static bool haveDecoupledModeInMol(const gmx_moltype_t& molt,
1508 gmx::ArrayRef<const t_iparams> iparams,
1509 real massFactorThreshold)
1511 if (molt.ilist[F_CONSTR].empty() && molt.ilist[F_CONSTRNC].empty())
1516 const t_atom* atom = molt.atoms.atom;
1518 const auto atomToConstraints =
1519 gmx::make_at2con(molt, iparams, gmx::FlexibleConstraintTreatment::Exclude);
1521 bool haveDecoupledMode = false;
1522 for (int ftype = 0; ftype < F_NRE; ftype++)
1524 if (interaction_function[ftype].flags & IF_ATYPE)
1526 const int nral = NRAL(ftype);
1527 const InteractionList& il = molt.ilist[ftype];
1528 for (int i = 0; i < il.size(); i += 1 + nral)
1530 /* Here we check for the mass difference between the atoms
1531 * at both ends of the angle, that the atoms at the ends have
1532 * 1 contraint and the atom in the middle at least 3; we check
1533 * that the 3 atoms are linked by constraints below.
1534 * We check for at least three constraints for the middle atom,
1535 * since with only the two bonds in the angle, we have 3-atom
1536 * molecule, which has much more thermal exhange in this single
1537 * angle mode than molecules with more atoms.
1538 * Note that this check also catches molecules where more atoms
1539 * are connected to one or more atoms in the angle, but by
1540 * bond potentials instead of angles. But such cases will not
1541 * occur in "normal" molecules and it doesn't hurt running
1542 * those with higher accuracy settings as well.
1544 int a0 = il.iatoms[1 + i];
1545 int a1 = il.iatoms[1 + i + 1];
1546 int a2 = il.iatoms[1 + i + 2];
1547 if ((atom[a0].m > atom[a2].m * massFactorThreshold || atom[a2].m > atom[a0].m * massFactorThreshold)
1548 && atomToConstraints[a0].ssize() == 1 && atomToConstraints[a2].ssize() == 1
1549 && atomToConstraints[a1].ssize() >= 3)
1551 int constraint0 = atomToConstraints[a0][0];
1552 int constraint2 = atomToConstraints[a2][0];
1554 bool foundAtom0 = false;
1555 bool foundAtom2 = false;
1556 for (const int constraint : atomToConstraints[a1])
1558 if (constraint == constraint0)
1562 if (constraint == constraint2)
1567 if (foundAtom0 && foundAtom2)
1569 haveDecoupledMode = true;
1576 return haveDecoupledMode;
1579 /*! \brief Checks if the Verlet buffer and constraint accuracy is sufficient for decoupled dynamic modes.
1581 * When decoupled modes are present and the accuray in insufficient,
1582 * this routine issues a warning if the accuracy is insufficient.
1584 static void checkDecoupledModeAccuracy(const gmx_mtop_t* mtop, const t_inputrec* ir, warninp* wi)
1586 /* We only have issues with decoupled modes with normal MD.
1587 * With stochastic dynamics equipartitioning is enforced strongly.
1594 /* When atoms of very different mass are involved in an angle potential
1595 * and both bonds in the angle are constrained, the dynamic modes in such
1596 * angles have very different periods and significant energy exchange
1597 * takes several nanoseconds. Thus even a small amount of error in
1598 * different algorithms can lead to violation of equipartitioning.
1599 * The parameters below are mainly based on an all-atom chloroform model
1600 * with all bonds constrained. Equipartitioning is off by more than 1%
1601 * (need to run 5-10 ns) when the difference in mass between H and Cl
1602 * is more than a factor 13 and the accuracy is less than the thresholds
1603 * given below. This has been verified on some other molecules.
1605 * Note that the buffer and shake parameters have unit length and
1606 * energy/time, respectively, so they will "only" work correctly
1607 * for atomistic force fields using MD units.
1609 const real massFactorThreshold = 13.0;
1610 const real bufferToleranceThreshold = 1e-4;
1611 const int lincsIterationThreshold = 2;
1612 const int lincsOrderThreshold = 4;
1613 const real shakeToleranceThreshold = 0.005 * ir->delta_t;
1615 bool lincsWithSufficientTolerance =
1616 (ir->eConstrAlg == ConstraintAlgorithm::Lincs
1617 && ir->nLincsIter >= lincsIterationThreshold && ir->nProjOrder >= lincsOrderThreshold);
1618 bool shakeWithSufficientTolerance = (ir->eConstrAlg == ConstraintAlgorithm::Shake
1619 && ir->shake_tol <= 1.1 * shakeToleranceThreshold);
1620 if (ir->cutoff_scheme == CutoffScheme::Verlet && ir->verletbuf_tol <= 1.1 * bufferToleranceThreshold
1621 && (lincsWithSufficientTolerance || shakeWithSufficientTolerance))
1626 bool haveDecoupledMode = false;
1627 for (const gmx_moltype_t& molt : mtop->moltype)
1629 if (haveDecoupledModeInMol(molt, mtop->ffparams.iparams, massFactorThreshold))
1631 haveDecoupledMode = true;
1635 if (haveDecoupledMode)
1637 std::string message = gmx::formatString(
1638 "There are atoms at both ends of an angle, connected by constraints "
1639 "and with masses that differ by more than a factor of %g. This means "
1640 "that there are likely dynamic modes that are only very weakly coupled.",
1641 massFactorThreshold);
1642 if (ir->cutoff_scheme == CutoffScheme::Verlet)
1644 message += gmx::formatString(
1645 " To ensure good equipartitioning, you need to either not use "
1646 "constraints on all bonds (but, if possible, only on bonds involving "
1647 "hydrogens) or use integrator = %s or decrease one or more tolerances: "
1648 "verlet-buffer-tolerance <= %g, LINCS iterations >= %d, LINCS order "
1649 ">= %d or SHAKE tolerance <= %g",
1650 enumValueToString(IntegrationAlgorithm::SD1),
1651 bufferToleranceThreshold,
1652 lincsIterationThreshold,
1653 lincsOrderThreshold,
1654 shakeToleranceThreshold);
1658 message += gmx::formatString(
1659 " To ensure good equipartitioning, we suggest to switch to the %s "
1660 "cutoff-scheme, since that allows for better control over the Verlet "
1661 "buffer size and thus over the energy drift.",
1662 enumValueToString(CutoffScheme::Verlet));
1664 warning(wi, message);
1668 static void set_verlet_buffer(const gmx_mtop_t* mtop,
1673 const gmx::MDLogger& logger)
1675 GMX_LOG(logger.info)
1677 .appendTextFormatted(
1678 "Determining Verlet buffer for a tolerance of %g kJ/mol/ps at %g K",
1682 /* Calculate the buffer size for simple atom vs atoms list */
1683 VerletbufListSetup listSetup1x1;
1684 listSetup1x1.cluster_size_i = 1;
1685 listSetup1x1.cluster_size_j = 1;
1686 const real rlist_1x1 = calcVerletBufferSize(
1687 *mtop, det(box), *ir, ir->nstlist, ir->nstlist - 1, buffer_temp, listSetup1x1);
1689 /* Set the pair-list buffer size in ir */
1690 VerletbufListSetup listSetup4x4 = verletbufGetSafeListSetup(ListSetupType::CpuNoSimd);
1691 ir->rlist = calcVerletBufferSize(
1692 *mtop, det(box), *ir, ir->nstlist, ir->nstlist - 1, buffer_temp, listSetup4x4);
1694 const int n_nonlin_vsite = gmx::countNonlinearVsites(*mtop);
1695 if (n_nonlin_vsite > 0)
1697 std::string warningMessage = gmx::formatString(
1698 "There are %d non-linear virtual site constructions. Their contribution to the "
1699 "energy error is approximated. In most cases this does not affect the error "
1702 warning_note(wi, warningMessage);
1705 GMX_LOG(logger.info)
1707 .appendTextFormatted(
1708 "Calculated rlist for %dx%d atom pair-list as %.3f nm, buffer size %.3f nm",
1712 rlist_1x1 - std::max(ir->rvdw, ir->rcoulomb));
1714 GMX_LOG(logger.info)
1716 .appendTextFormatted(
1717 "Set rlist, assuming %dx%d atom pair-list, to %.3f nm, buffer size %.3f nm",
1718 listSetup4x4.cluster_size_i,
1719 listSetup4x4.cluster_size_j,
1721 ir->rlist - std::max(ir->rvdw, ir->rcoulomb));
1723 GMX_LOG(logger.info)
1725 .appendTextFormatted(
1726 "Note that mdrun will redetermine rlist based on the actual pair-list setup");
1728 if (gmx::square(ir->rlist) >= max_cutoff2(ir->pbcType, box))
1731 "The pair-list cut-off (%g nm) is longer than half the shortest box vector or "
1732 "longer than the smallest box diagonal element (%g nm). Increase the box size or "
1733 "decrease nstlist or increase verlet-buffer-tolerance.",
1735 std::sqrt(max_cutoff2(ir->pbcType, box)));
1739 int gmx_grompp(int argc, char* argv[])
1741 const char* desc[] = {
1742 "[THISMODULE] (the gromacs preprocessor)",
1743 "reads a molecular topology file, checks the validity of the",
1744 "file, expands the topology from a molecular description to an atomic",
1745 "description. The topology file contains information about",
1746 "molecule types and the number of molecules, the preprocessor",
1747 "copies each molecule as needed. ",
1748 "There is no limitation on the number of molecule types. ",
1749 "Bonds and bond-angles can be converted into constraints, separately",
1750 "for hydrogens and heavy atoms.",
1751 "Then a coordinate file is read and velocities can be generated",
1752 "from a Maxwellian distribution if requested.",
1753 "[THISMODULE] also reads parameters for [gmx-mdrun] ",
1754 "(eg. number of MD steps, time step, cut-off).",
1755 "Eventually a binary file is produced that can serve as the sole input",
1756 "file for the MD program.[PAR]",
1758 "[THISMODULE] uses the atom names from the topology file. The atom names",
1759 "in the coordinate file (option [TT]-c[tt]) are only read to generate",
1760 "warnings when they do not match the atom names in the topology.",
1761 "Note that the atom names are irrelevant for the simulation as",
1762 "only the atom types are used for generating interaction parameters.[PAR]",
1764 "[THISMODULE] uses a built-in preprocessor to resolve includes, macros, ",
1765 "etc. The preprocessor supports the following keywords::",
1768 " #ifndef VARIABLE",
1771 " #define VARIABLE",
1773 " #include \"filename\"",
1774 " #include <filename>",
1776 "The functioning of these statements in your topology may be modulated by",
1777 "using the following two flags in your [REF].mdp[ref] file::",
1779 " define = -DVARIABLE1 -DVARIABLE2",
1780 " include = -I/home/john/doe",
1782 "For further information a C-programming textbook may help you out.",
1783 "Specifying the [TT]-pp[tt] flag will get the pre-processed",
1784 "topology file written out so that you can verify its contents.[PAR]",
1786 "When using position restraints, a file with restraint coordinates",
1787 "must be supplied with [TT]-r[tt] (can be the same file as supplied",
1788 "for [TT]-c[tt]). For free energy calculations, separate reference",
1789 "coordinates for the B topology can be supplied with [TT]-rb[tt],",
1790 "otherwise they will be equal to those of the A topology.[PAR]",
1792 "Starting coordinates can be read from trajectory with [TT]-t[tt].",
1793 "The last frame with coordinates and velocities will be read,",
1794 "unless the [TT]-time[tt] option is used. Only if this information",
1795 "is absent will the coordinates in the [TT]-c[tt] file be used.",
1796 "Note that these velocities will not be used when [TT]gen_vel = yes[tt]",
1797 "in your [REF].mdp[ref] file. An energy file can be supplied with",
1798 "[TT]-e[tt] to read Nose-Hoover and/or Parrinello-Rahman coupling",
1801 "[THISMODULE] can be used to restart simulations (preserving",
1802 "continuity) by supplying just a checkpoint file with [TT]-t[tt].",
1803 "However, for simply changing the number of run steps to extend",
1804 "a run, using [gmx-convert-tpr] is more convenient than [THISMODULE].",
1805 "You then supply the old checkpoint file directly to [gmx-mdrun]",
1806 "with [TT]-cpi[tt]. If you wish to change the ensemble or things",
1807 "like output frequency, then supplying the checkpoint file to",
1808 "[THISMODULE] with [TT]-t[tt] along with a new [REF].mdp[ref] file",
1809 "with [TT]-f[tt] is the recommended procedure. Actually preserving",
1810 "the ensemble (if possible) still requires passing the checkpoint",
1811 "file to [gmx-mdrun] [TT]-cpi[tt].[PAR]",
1813 "By default, all bonded interactions which have constant energy due to",
1814 "virtual site constructions will be removed. If this constant energy is",
1815 "not zero, this will result in a shift in the total energy. All bonded",
1816 "interactions can be kept by turning off [TT]-rmvsbds[tt]. Additionally,",
1817 "all constraints for distances which will be constant anyway because",
1818 "of virtual site constructions will be removed. If any constraints remain",
1819 "which involve virtual sites, a fatal error will result.[PAR]",
1821 "To verify your run input file, please take note of all warnings",
1822 "on the screen, and correct where necessary. Do also look at the contents",
1823 "of the [TT]mdout.mdp[tt] file; this contains comment lines, as well as",
1824 "the input that [THISMODULE] has read. If in doubt, you can start [THISMODULE]",
1825 "with the [TT]-debug[tt] option which will give you more information",
1826 "in a file called [TT]grompp.log[tt] (along with real debug info). You",
1827 "can see the contents of the run input file with the [gmx-dump]",
1828 "program. [gmx-check] can be used to compare the contents of two",
1829 "run input files.[PAR]",
1831 "The [TT]-maxwarn[tt] option can be used to override warnings printed",
1832 "by [THISMODULE] that otherwise halt output. In some cases, warnings are",
1833 "harmless, but usually they are not. The user is advised to carefully",
1834 "interpret the output messages before attempting to bypass them with",
1837 std::vector<MoleculeInformation> mi;
1838 std::unique_ptr<MoleculeInformation> intermolecular_interactions;
1840 CombinationRule comb;
1843 const char* mdparin;
1844 bool bNeedVel, bGenVel;
1845 gmx_output_env_t* oenv;
1846 gmx_bool bVerbose = FALSE;
1849 t_filenm fnm[] = { { efMDP, nullptr, nullptr, ffREAD },
1850 { efMDP, "-po", "mdout", ffWRITE },
1851 { efSTX, "-c", nullptr, ffREAD },
1852 { efSTX, "-r", "restraint", ffOPTRD },
1853 { efSTX, "-rb", "restraint", ffOPTRD },
1854 { efNDX, nullptr, nullptr, ffOPTRD },
1855 { efTOP, nullptr, nullptr, ffREAD },
1856 { efTOP, "-pp", "processed", ffOPTWR },
1857 { efTPR, "-o", nullptr, ffWRITE },
1858 { efTRN, "-t", nullptr, ffOPTRD },
1859 { efEDR, "-e", nullptr, ffOPTRD },
1860 /* This group is needed by the VMD viewer as the start configuration for IMD sessions: */
1861 { efGRO, "-imd", "imdgroup", ffOPTWR },
1862 { efTRN, "-ref", "rotref", ffOPTRW | ffALLOW_MISSING } };
1863 #define NFILE asize(fnm)
1865 /* Command line options */
1866 gmx_bool bRenum = TRUE;
1867 gmx_bool bRmVSBds = TRUE, bZero = FALSE;
1871 { "-v", FALSE, etBOOL, { &bVerbose }, "Be loud and noisy" },
1872 { "-time", FALSE, etREAL, { &fr_time }, "Take frame at or first after this time." },
1877 "Remove constant bonded interactions with virtual sites" },
1882 "Number of allowed warnings during input processing. Not for normal use and may "
1883 "generate unstable systems" },
1888 "Set parameters for bonded interactions without defaults to zero instead of "
1889 "generating an error" },
1894 "Renumber atomtypes and minimize number of atomtypes" }
1897 /* Parse the command line */
1898 if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
1903 /* Initiate some variables */
1904 gmx::MDModules mdModules;
1905 t_inputrec irInstance;
1906 t_inputrec* ir = &irInstance;
1907 t_gromppopts optsInstance;
1908 t_gromppopts* opts = &optsInstance;
1909 snew(opts->include, STRLEN);
1910 snew(opts->define, STRLEN);
1912 gmx::LoggerBuilder builder;
1913 builder.addTargetStream(gmx::MDLogger::LogLevel::Info, &gmx::TextOutputFile::standardOutput());
1914 builder.addTargetStream(gmx::MDLogger::LogLevel::Warning, &gmx::TextOutputFile::standardError());
1915 gmx::LoggerOwner logOwner(builder.build());
1916 const gmx::MDLogger logger(logOwner.logger());
1919 wi = init_warning(TRUE, maxwarn);
1921 /* PARAMETER file processing */
1922 mdparin = opt2fn("-f", NFILE, fnm);
1923 set_warning_line(wi, mdparin, -1);
1926 get_ir(mdparin, opt2fn("-po", NFILE, fnm), &mdModules, ir, opts, WriteMdpHeader::yes, wi);
1928 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
1930 // Now that the MDModules have their options assigned from get_ir, subscribe
1931 // to eventual notifications during pre-processing their data
1932 mdModules.subscribeToPreProcessingNotifications();
1937 GMX_LOG(logger.info)
1939 .appendTextFormatted("checking input for internal consistency...");
1941 check_ir(mdparin, mdModules.notifiers(), ir, opts, wi);
1943 if (ir->ld_seed == -1)
1945 ir->ld_seed = static_cast<int>(gmx::makeRandomSeed());
1946 GMX_LOG(logger.info)
1948 .appendTextFormatted("Setting the LD random seed to %" PRId64 "", ir->ld_seed);
1951 if (ir->expandedvals->lmc_seed == -1)
1953 ir->expandedvals->lmc_seed = static_cast<int>(gmx::makeRandomSeed());
1954 GMX_LOG(logger.info)
1956 .appendTextFormatted("Setting the lambda MC random seed to %d", ir->expandedvals->lmc_seed);
1959 bNeedVel = EI_STATE_VELOCITY(ir->eI);
1960 bGenVel = (bNeedVel && opts->bGenVel);
1961 if (bGenVel && ir->bContinuation)
1963 std::string warningMessage = gmx::formatString(
1964 "Generating velocities is inconsistent with attempting "
1965 "to continue a previous run. Choose only one of "
1966 "gen-vel = yes and continuation = yes.");
1967 warning_error(wi, warningMessage);
1970 std::array<InteractionsOfType, F_NRE> interactions;
1972 PreprocessingAtomTypes atypes;
1975 pr_symtab(debug, 0, "Just opened", &sys.symtab);
1978 const char* fn = ftp2fn(efTOP, NFILE, fnm);
1979 if (!gmx_fexist(fn))
1981 gmx_fatal(FARGS, "%s does not exist", fn);
1986 opt2fn_null("-pp", NFILE, fnm),
1987 opt2fn("-c", NFILE, fnm),
1997 &intermolecular_interactions,
2008 pr_symtab(debug, 0, "After new_status", &sys.symtab);
2012 /* set parameters for virtual site construction (not for vsiten) */
2013 for (size_t mt = 0; mt < sys.moltype.size(); mt++)
2015 nvsite += set_vsites(bVerbose, &sys.moltype[mt].atoms, &atypes, mi[mt].interactions, logger);
2017 /* now throw away all obsolete bonds, angles and dihedrals: */
2018 /* note: constraints are ALWAYS removed */
2021 for (size_t mt = 0; mt < sys.moltype.size(); mt++)
2023 clean_vsite_bondeds(mi[mt].interactions, sys.moltype[mt].atoms.nr, bRmVSBds, logger);
2027 if ((count_constraints(&sys, mi, wi) != 0) && (ir->eConstrAlg == ConstraintAlgorithm::Shake))
2029 if (ir->eI == IntegrationAlgorithm::CG || ir->eI == IntegrationAlgorithm::LBFGS)
2031 std::string warningMessage =
2032 gmx::formatString("Can not do %s with %s, use %s",
2033 enumValueToString(ir->eI),
2034 enumValueToString(ConstraintAlgorithm::Shake),
2035 enumValueToString(ConstraintAlgorithm::Lincs));
2036 warning_error(wi, warningMessage);
2038 if (ir->bPeriodicMols)
2040 std::string warningMessage =
2041 gmx::formatString("Can not do periodic molecules with %s, use %s",
2042 enumValueToString(ConstraintAlgorithm::Shake),
2043 enumValueToString(ConstraintAlgorithm::Lincs));
2044 warning_error(wi, warningMessage);
2048 if (EI_SD(ir->eI) && ir->etc != TemperatureCoupling::No)
2050 warning_note(wi, "Temperature coupling is ignored with SD integrators.");
2053 /* Check for errors in the input now, since they might cause problems
2054 * during processing further down.
2056 check_warning_error(wi, FARGS);
2058 if (nint_ftype(&sys, mi, F_POSRES) > 0 || nint_ftype(&sys, mi, F_FBPOSRES) > 0)
2060 if (ir->epc == PressureCoupling::ParrinelloRahman || ir->epc == PressureCoupling::Mttk)
2062 std::string warningMessage = gmx::formatString(
2063 "You are combining position restraints with %s pressure coupling, which can "
2064 "lead to instabilities. If you really want to combine position restraints with "
2065 "pressure coupling, we suggest to use %s pressure coupling instead.",
2066 enumValueToString(ir->epc),
2067 enumValueToString(PressureCoupling::Berendsen));
2068 warning_note(wi, warningMessage);
2071 const char* fn = opt2fn("-r", NFILE, fnm);
2074 if (!gmx_fexist(fn))
2077 "Cannot find position restraint file %s (option -r).\n"
2078 "From GROMACS-2018, you need to specify the position restraint "
2079 "coordinate files explicitly to avoid mistakes, although you can "
2080 "still use the same file as you specify for the -c option.",
2084 if (opt2bSet("-rb", NFILE, fnm))
2086 fnB = opt2fn("-rb", NFILE, fnm);
2087 if (!gmx_fexist(fnB))
2090 "Cannot find B-state position restraint file %s (option -rb).\n"
2091 "From GROMACS-2018, you need to specify the position restraint "
2092 "coordinate files explicitly to avoid mistakes, although you can "
2093 "still use the same file as you specify for the -c option.",
2104 std::string message = gmx::formatString("Reading position restraint coords from %s", fn);
2105 if (strcmp(fn, fnB) != 0)
2107 message += gmx::formatString(" and %s", fnB);
2109 GMX_LOG(logger.info).asParagraph().appendText(message);
2111 gen_posres(&sys, mi, fn, fnB, ir->refcoord_scaling, ir->pbcType, ir->posres_com, ir->posres_comB, wi, logger);
2114 /* If we are using CMAP, setup the pre-interpolation grid */
2115 if (interactions[F_CMAP].ncmap() > 0)
2117 init_cmap_grid(&sys.ffparams.cmap_grid,
2118 interactions[F_CMAP].cmapAngles,
2119 interactions[F_CMAP].cmakeGridSpacing);
2120 setup_cmap(interactions[F_CMAP].cmakeGridSpacing,
2121 interactions[F_CMAP].cmapAngles,
2122 interactions[F_CMAP].cmap,
2123 &sys.ffparams.cmap_grid);
2126 set_wall_atomtype(&atypes, opts, ir, wi, logger);
2129 atypes.renumberTypes(interactions, &sys, ir->wall_atomtype, bVerbose);
2132 if (ir->implicit_solvent)
2134 gmx_fatal(FARGS, "Implicit solvation is no longer supported");
2137 /* PELA: Copy the atomtype data to the topology atomtype list */
2138 atypes.copyTot_atomtypes(&(sys.atomtypes));
2142 pr_symtab(debug, 0, "After atype.renumberTypes", &sys.symtab);
2147 GMX_LOG(logger.info).asParagraph().appendTextFormatted("converting bonded parameters...");
2150 const int ntype = atypes.size();
2151 convertInteractionsOfType(
2152 ntype, interactions, mi, intermolecular_interactions.get(), comb, reppow, fudgeQQ, &sys);
2156 pr_symtab(debug, 0, "After converInteractionsOfType", &sys.symtab);
2159 /* set ptype to VSite for virtual sites */
2160 for (gmx_moltype_t& moltype : sys.moltype)
2162 set_vsites_ptype(FALSE, &moltype, logger);
2166 pr_symtab(debug, 0, "After virtual sites", &sys.symtab);
2168 /* Check velocity for virtual sites and shells */
2171 check_vel(&sys, state.v.rvec_array());
2174 /* check for shells and inpurecs */
2175 check_shells_inputrec(&sys, ir, wi);
2178 check_mol(&sys, wi);
2180 checkForUnboundAtoms(&sys, bVerbose, wi, logger);
2182 if (EI_DYNAMICS(ir->eI) && ir->eI != IntegrationAlgorithm::BD)
2184 check_bonds_timestep(&sys, ir->delta_t, wi);
2187 checkDecoupledModeAccuracy(&sys, ir, wi);
2189 if (EI_ENERGY_MINIMIZATION(ir->eI) && 0 == ir->nsteps)
2193 "Zero-step energy minimization will alter the coordinates before calculating the "
2194 "energy. If you just want the energy of a single point, try zero-step MD (with "
2195 "unconstrained_start = yes). To do multiple single-point energy evaluations of "
2196 "different configurations of the same topology, use mdrun -rerun.");
2199 check_warning_error(wi, FARGS);
2203 GMX_LOG(logger.info).asParagraph().appendTextFormatted("initialising group options...");
2205 do_index(mdparin, ftp2fn_null(efNDX, NFILE, fnm), &sys, bVerbose, mdModules.notifiers(), ir, wi);
2207 if (ir->cutoff_scheme == CutoffScheme::Verlet && ir->verletbuf_tol > 0)
2209 if (EI_DYNAMICS(ir->eI) && inputrec2nboundeddim(ir) == 3)
2213 if (EI_MD(ir->eI) && ir->etc == TemperatureCoupling::No)
2217 buffer_temp = opts->tempi;
2221 buffer_temp = calc_temp(&sys, ir, state.v.rvec_array());
2223 if (buffer_temp > 0)
2225 std::string warningMessage = gmx::formatString(
2226 "NVE simulation: will use the initial temperature of %.3f K for "
2227 "determining the Verlet buffer size",
2229 warning_note(wi, warningMessage);
2233 std::string warningMessage = gmx::formatString(
2234 "NVE simulation with an initial temperature of zero: will use a Verlet "
2235 "buffer of %d%%. Check your energy drift!",
2236 gmx::roundToInt(verlet_buffer_ratio_NVE_T0 * 100));
2237 warning_note(wi, warningMessage);
2242 buffer_temp = get_max_reference_temp(ir, wi);
2245 if (EI_MD(ir->eI) && ir->etc == TemperatureCoupling::No && buffer_temp == 0)
2247 /* NVE with initial T=0: we add a fixed ratio to rlist.
2248 * Since we don't actually use verletbuf_tol, we set it to -1
2249 * so it can't be misused later.
2251 ir->rlist *= 1.0 + verlet_buffer_ratio_NVE_T0;
2252 ir->verletbuf_tol = -1;
2256 /* We warn for NVE simulations with a drift tolerance that
2257 * might result in a 1(.1)% drift over the total run-time.
2258 * Note that we can't warn when nsteps=0, since we don't
2259 * know how many steps the user intends to run.
2261 if (EI_MD(ir->eI) && ir->etc == TemperatureCoupling::No && ir->nstlist > 1 && ir->nsteps > 0)
2263 const real driftTolerance = 0.01;
2264 /* We use 2 DOF per atom = 2kT pot+kin energy,
2265 * to be on the safe side with constraints.
2267 const real totalEnergyDriftPerAtomPerPicosecond =
2268 2 * gmx::c_boltz * buffer_temp / (ir->nsteps * ir->delta_t);
2270 if (ir->verletbuf_tol > 1.1 * driftTolerance * totalEnergyDriftPerAtomPerPicosecond)
2272 std::string warningMessage = gmx::formatString(
2273 "You are using a Verlet buffer tolerance of %g kJ/mol/ps for an "
2274 "NVE simulation of length %g ps, which can give a final drift of "
2275 "%d%%. For conserving energy to %d%% when using constraints, you "
2276 "might need to set verlet-buffer-tolerance to %.1e.",
2278 ir->nsteps * ir->delta_t,
2279 gmx::roundToInt(ir->verletbuf_tol / totalEnergyDriftPerAtomPerPicosecond * 100),
2280 gmx::roundToInt(100 * driftTolerance),
2281 driftTolerance * totalEnergyDriftPerAtomPerPicosecond);
2282 warning_note(wi, warningMessage);
2286 set_verlet_buffer(&sys, ir, buffer_temp, state.box, wi, logger);
2291 /* Init the temperature coupling state */
2292 init_gtc_state(&state, ir->opts.ngtc, 0, ir->opts.nhchainlength); /* need to add nnhpres here? */
2296 pr_symtab(debug, 0, "After index", &sys.symtab);
2299 triple_check(mdparin, ir, &sys, wi);
2300 close_symtab(&sys.symtab);
2303 pr_symtab(debug, 0, "After close", &sys.symtab);
2306 if (ir->eI == IntegrationAlgorithm::Mimic)
2308 generate_qmexcl(&sys, ir, logger);
2311 if (ftp2bSet(efTRN, NFILE, fnm))
2315 GMX_LOG(logger.info)
2317 .appendTextFormatted("getting data from old trajectory ...");
2319 cont_status(ftp2fn(efTRN, NFILE, fnm),
2320 ftp2fn_null(efEDR, NFILE, fnm),
2331 if (ir->pbcType == PbcType::XY && ir->nwall != 2)
2333 clear_rvec(state.box[ZZ]);
2336 if (EEL_FULL(ir->coulombtype) || EVDW_PME(ir->vdwtype))
2338 /* Calculate the optimal grid dimensions */
2340 EwaldBoxZScaler boxScaler(*ir);
2341 boxScaler.scaleBox(state.box, scaledBox);
2343 if (ir->nkx > 0 && ir->nky > 0 && ir->nkz > 0)
2345 /* Mark fourier_spacing as not used */
2346 ir->fourier_spacing = 0;
2348 else if (ir->nkx != 0 && ir->nky != 0 && ir->nkz != 0)
2350 set_warning_line(wi, mdparin, -1);
2352 wi, "Some of the Fourier grid sizes are set, but all of them need to be set.");
2354 const int minGridSize = minimalPmeGridSize(ir->pme_order);
2355 calcFftGrid(stdout, scaledBox, ir->fourier_spacing, minGridSize, &(ir->nkx), &(ir->nky), &(ir->nkz));
2356 if (ir->nkx < minGridSize || ir->nky < minGridSize || ir->nkz < minGridSize)
2359 "The PME grid size should be >= 2*(pme-order - 1); either manually "
2360 "increase the grid size or decrease pme-order");
2364 /* MRS: eventually figure out better logic for initializing the fep
2365 values that makes declaring the lambda and declaring the state not
2366 potentially conflict if not handled correctly. */
2367 if (ir->efep != FreeEnergyPerturbationType::No)
2369 state.fep_state = ir->fepvals->init_fep_state;
2370 for (i = 0; i < static_cast<int>(FreeEnergyPerturbationCouplingType::Count); i++)
2372 /* init_lambda trumps state definitions*/
2373 if (ir->fepvals->init_lambda >= 0)
2375 state.lambda[i] = ir->fepvals->init_lambda;
2379 if (ir->fepvals->all_lambda[i].empty())
2381 gmx_fatal(FARGS, "Values of lambda not set for a free energy calculation!");
2385 state.lambda[i] = ir->fepvals->all_lambda[i][state.fep_state];
2391 pull_t* pull = nullptr;
2395 pull = set_pull_init(
2396 ir, sys, state.x, state.box, state.lambda[FreeEnergyPerturbationCouplingType::Mass], wi);
2399 /* Modules that supply external potential for pull coordinates
2400 * should register those potentials here. finish_pull() will check
2401 * that providers have been registerd for all external potentials.
2406 tensor compressibility = { { 0 } };
2407 if (ir->epc != PressureCoupling::No)
2409 copy_mat(ir->compress, compressibility);
2411 setStateDependentAwhParams(
2412 ir->awhParams.get(),
2419 ir->efep != FreeEnergyPerturbationType::No
2420 ? ir->fepvals->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Fep)]
2421 [ir->fepvals->init_fep_state]
2434 set_reference_positions(ir->rot,
2435 state.x.rvec_array(),
2437 opt2fn("-ref", NFILE, fnm),
2438 opt2bSet("-ref", NFILE, fnm),
2442 /* reset_multinr(sys); */
2444 if (EEL_PME(ir->coulombtype))
2446 float ratio = pme_load_estimate(sys, *ir, state.box);
2447 GMX_LOG(logger.info)
2449 .appendTextFormatted(
2450 "Estimate for the relative computational load of the PME mesh part: %.2f", ratio);
2451 /* With free energy we might need to do PME both for the A and B state
2452 * charges. This will double the cost, but the optimal performance will
2453 * then probably be at a slightly larger cut-off and grid spacing.
2455 if ((ir->efep == FreeEnergyPerturbationType::No && ratio > 1.0 / 2.0)
2456 || (ir->efep != FreeEnergyPerturbationType::No && ratio > 2.0 / 3.0))
2460 "The optimal PME mesh load for parallel simulations is below 0.5\n"
2461 "and for highly parallel simulations between 0.25 and 0.33,\n"
2462 "for higher performance, increase the cut-off and the PME grid spacing.\n");
2463 if (ir->efep != FreeEnergyPerturbationType::No)
2466 "For free energy simulations, the optimal load limit increases from "
2473 double cio = compute_io(ir, sys.natoms, sys.groups, F_NRE, 1);
2474 std::string warningMessage =
2475 gmx::formatString("This run will generate roughly %.0f Mb of data", cio);
2478 set_warning_line(wi, mdparin, -1);
2479 warning_note(wi, warningMessage);
2483 GMX_LOG(logger.info).asParagraph().appendText(warningMessage);
2487 // Add the md modules internal parameters that are not mdp options
2488 // e.g., atom indices
2491 gmx::KeyValueTreeBuilder internalParameterBuilder;
2492 mdModules.notifiers().preProcessingNotifier_.notify(internalParameterBuilder.rootObject());
2493 ir->internalParameters =
2494 std::make_unique<gmx::KeyValueTreeObject>(internalParameterBuilder.build());
2497 if (ir->comm_mode != ComRemovalAlgorithm::No)
2499 const int nstglobalcomm = computeGlobalCommunicationPeriod(ir);
2500 if (ir->nstcomm % nstglobalcomm != 0)
2505 "COM removal frequency is set to (%d).\n"
2506 "Other settings require a global communication frequency of %d.\n"
2507 "Note that this will require additional global communication steps,\n"
2508 "which will reduce performance when using multiple ranks.\n"
2509 "Consider setting nstcomm to a multiple of %d.",
2518 GMX_LOG(logger.info).asParagraph().appendTextFormatted("writing run input file...");
2521 done_warning(wi, FARGS);
2522 write_tpx_state(ftp2fn(efTPR, NFILE, fnm), ir, &state, sys);
2524 /* Output IMD group, if bIMD is TRUE */
2525 gmx::write_IMDgroup_to_file(ir->bIMD, ir, &state, sys, NFILE, fnm);
2527 sfree(opts->define);
2528 sfree(opts->wall_atomtype[0]);
2529 sfree(opts->wall_atomtype[1]);
2530 sfree(opts->include);
2531 sfree(opts->couple_moltype);
2533 for (auto& mol : mi)
2535 // Some of the contents of molinfo have been stolen, so
2536 // fullCleanUp can't be called.
2537 mol.partialCleanUp();
2539 done_inputrec_strings();
2540 output_env_done(oenv);