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,2018,2019, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
50 #include <sys/types.h>
52 #include "gromacs/awh/read_params.h"
53 #include "gromacs/commandline/pargs.h"
54 #include "gromacs/ewald/ewald_utils.h"
55 #include "gromacs/ewald/pme.h"
56 #include "gromacs/fft/calcgrid.h"
57 #include "gromacs/fileio/confio.h"
58 #include "gromacs/fileio/enxio.h"
59 #include "gromacs/fileio/tpxio.h"
60 #include "gromacs/fileio/trxio.h"
61 #include "gromacs/fileio/warninp.h"
62 #include "gromacs/gmxpreprocess/add_par.h"
63 #include "gromacs/gmxpreprocess/convparm.h"
64 #include "gromacs/gmxpreprocess/gen_maxwell_velocities.h"
65 #include "gromacs/gmxpreprocess/gpp_atomtype.h"
66 #include "gromacs/gmxpreprocess/grompp_impl.h"
67 #include "gromacs/gmxpreprocess/notset.h"
68 #include "gromacs/gmxpreprocess/readir.h"
69 #include "gromacs/gmxpreprocess/tomorse.h"
70 #include "gromacs/gmxpreprocess/topio.h"
71 #include "gromacs/gmxpreprocess/toputil.h"
72 #include "gromacs/gmxpreprocess/vsite_parm.h"
73 #include "gromacs/imd/imd.h"
74 #include "gromacs/math/functions.h"
75 #include "gromacs/math/invertmatrix.h"
76 #include "gromacs/math/units.h"
77 #include "gromacs/math/vec.h"
78 #include "gromacs/mdlib/calc_verletbuf.h"
79 #include "gromacs/mdlib/compute_io.h"
80 #include "gromacs/mdlib/constr.h"
81 #include "gromacs/mdlib/perf_est.h"
82 #include "gromacs/mdlib/qmmm.h"
83 #include "gromacs/mdlib/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/futil.h"
103 #include "gromacs/utility/gmxassert.h"
104 #include "gromacs/utility/smalloc.h"
105 #include "gromacs/utility/snprintf.h"
107 /* TODO The implementation details should move to their own source file. */
108 InteractionOfType::InteractionOfType(gmx::ArrayRef<const int> atoms,
109 gmx::ArrayRef<const real> params,
110 const std::string &name)
111 : atoms_(atoms.begin(), atoms.end()),
112 interactionTypeName_(name)
114 GMX_RELEASE_ASSERT(params.size() <= forceParam_.size(),
115 gmx::formatString("Cannot have more parameters than the maximum number possible (%d)", MAXFORCEPARAM).c_str());
116 auto forceParamIt = forceParam_.begin();
117 for (const auto param : params)
119 *forceParamIt++ = param;
121 while (forceParamIt != forceParam_.end())
123 *forceParamIt++ = NOTSET;
127 const int &InteractionOfType::ai() const
129 GMX_RELEASE_ASSERT(!atoms_.empty(), "Need to have at least one atom set");
133 const int &InteractionOfType::aj() const
135 GMX_RELEASE_ASSERT(atoms_.size() > 1, "Need to have at least two atoms set");
139 const int &InteractionOfType::ak() const
141 GMX_RELEASE_ASSERT(atoms_.size() > 2, "Need to have at least three atoms set");
145 const int &InteractionOfType::al() const
147 GMX_RELEASE_ASSERT(atoms_.size() > 3, "Need to have at least four atoms set");
151 const int &InteractionOfType::am() const
153 GMX_RELEASE_ASSERT(atoms_.size() > 4, "Need to have at least five atoms set");
157 const real &InteractionOfType::c0() const
159 return forceParam_[0];
162 const real &InteractionOfType::c1() const
164 return forceParam_[1];
167 const real &InteractionOfType::c2() const
169 return forceParam_[2];
172 const std::string &InteractionOfType::interactionTypeName() const
174 return interactionTypeName_;
177 void InteractionOfType::sortBondAtomIds()
181 std::swap(atoms_[0], atoms_[1]);
185 void InteractionOfType::sortAngleAtomIds()
189 std::swap(atoms_[0], atoms_[2]);
193 void InteractionOfType::sortDihedralAtomIds()
197 std::swap(atoms_[0], atoms_[3]);
198 std::swap(atoms_[1], atoms_[2]);
202 void InteractionOfType::sortAtomIds()
212 else if (isDihedral())
214 sortDihedralAtomIds();
218 GMX_THROW(gmx::InternalError("Can not sort parameters other than bonds, angles or dihedrals"));
222 void InteractionOfType::setForceParameter(int pos, real value)
224 GMX_RELEASE_ASSERT(pos < MAXFORCEPARAM, "Can't set parameter beyond the maximum number of parameters");
225 forceParam_[pos] = value;
228 void MoleculeInformation::initMolInfo()
233 init_t_atoms(&atoms, 0, FALSE);
236 void MoleculeInformation::partialCleanUp()
241 void MoleculeInformation::fullCleanUp()
248 static int rm_interactions(int ifunc, gmx::ArrayRef<MoleculeInformation> mols)
251 /* For all the molecule types */
252 for (auto &mol : mols)
254 n += mol.interactions[ifunc].size();
255 mol.interactions[ifunc].interactionTypes.clear();
260 static int check_atom_names(const char *fn1, const char *fn2,
261 gmx_mtop_t *mtop, const t_atoms *at)
263 int m, i, j, nmismatch;
265 #define MAXMISMATCH 20
267 if (mtop->natoms != at->nr)
269 gmx_incons("comparing atom names");
274 for (const gmx_molblock_t &molb : mtop->molblock)
276 tat = &mtop->moltype[molb.type].atoms;
277 for (m = 0; m < molb.nmol; m++)
279 for (j = 0; j < tat->nr; j++)
281 if (strcmp( *(tat->atomname[j]), *(at->atomname[i]) ) != 0)
283 if (nmismatch < MAXMISMATCH)
286 "Warning: atom name %d in %s and %s does not match (%s - %s)\n",
287 i+1, fn1, fn2, *(tat->atomname[j]), *(at->atomname[i]));
289 else if (nmismatch == MAXMISMATCH)
291 fprintf(stderr, "(more than %d non-matching atom names)\n", MAXMISMATCH);
303 static void check_eg_vs_cg(gmx_mtop_t *mtop)
305 int astart, m, cg, j, firstj;
306 unsigned char firsteg, eg;
309 /* Go through all the charge groups and make sure all their
310 * atoms are in the same energy group.
314 for (const gmx_molblock_t &molb : mtop->molblock)
316 molt = &mtop->moltype[molb.type];
317 for (m = 0; m < molb.nmol; m++)
319 for (cg = 0; cg < molt->cgs.nr; cg++)
321 /* Get the energy group of the first atom in this charge group */
322 firstj = astart + molt->cgs.index[cg];
323 firsteg = getGroupType(mtop->groups, SimulationAtomGroupType::EnergyOutput, firstj);
324 for (j = molt->cgs.index[cg]+1; j < molt->cgs.index[cg+1]; j++)
326 eg = getGroupType(mtop->groups, SimulationAtomGroupType::EnergyOutput, astart+j);
329 gmx_fatal(FARGS, "atoms %d and %d in charge group %d of molecule type '%s' are in different energy groups",
330 firstj+1, astart+j+1, cg+1, *molt->name);
334 astart += molt->atoms.nr;
339 static void check_cg_sizes(const char *topfn, const t_block *cgs, warninp *wi)
342 char warn_buf[STRLEN];
345 for (cg = 0; cg < cgs->nr; cg++)
347 maxsize = std::max(maxsize, cgs->index[cg+1]-cgs->index[cg]);
350 if (maxsize > MAX_CHARGEGROUP_SIZE)
352 gmx_fatal(FARGS, "The largest charge group contains %d atoms. The maximum is %d.", maxsize, MAX_CHARGEGROUP_SIZE);
354 else if (maxsize > 10)
356 set_warning_line(wi, topfn, -1);
358 "The largest charge group contains %d atoms.\n"
359 "Since atoms only see each other when the centers of geometry of the charge groups they belong to are within the cut-off distance, too large charge groups can lead to serious cut-off artifacts.\n"
360 "For efficiency and accuracy, charge group should consist of a few atoms.\n"
361 "For all-atom force fields use: CH3, CH2, CH, NH2, NH, OH, CO2, CO, etc.",
363 warning_note(wi, warn_buf);
367 static void check_bonds_timestep(const gmx_mtop_t *mtop, double dt, warninp *wi)
369 /* This check is not intended to ensure accurate integration,
370 * rather it is to signal mistakes in the mdp settings.
371 * A common mistake is to forget to turn on constraints
372 * for MD after energy minimization with flexible bonds.
373 * This check can also detect too large time steps for flexible water
374 * models, but such errors will often be masked by the constraints
375 * mdp options, which turns flexible water into water with bond constraints,
376 * but without an angle constraint. Unfortunately such incorrect use
377 * of water models can not easily be detected without checking
378 * for specific model names.
380 * The stability limit of leap-frog or velocity verlet is 4.44 steps
381 * per oscillational period.
382 * But accurate bonds distributions are lost far before that limit.
383 * To allow relatively common schemes (although not common with Gromacs)
384 * of dt=1 fs without constraints and dt=2 fs with only H-bond constraints
385 * we set the note limit to 10.
387 int min_steps_warn = 5;
388 int min_steps_note = 10;
390 int i, a1, a2, w_a1, w_a2, j;
391 real twopi2, limit2, fc, re, m1, m2, period2, w_period2;
392 bool bFound, bWater, bWarn;
393 char warn_buf[STRLEN];
395 /* Get the interaction parameters */
396 gmx::ArrayRef<const t_iparams> ip = mtop->ffparams.iparams;
398 twopi2 = gmx::square(2*M_PI);
400 limit2 = gmx::square(min_steps_note*dt);
405 const gmx_moltype_t *w_moltype = nullptr;
406 for (const gmx_moltype_t &moltype : mtop->moltype)
408 const t_atom *atom = moltype.atoms.atom;
409 const InteractionLists &ilist = moltype.ilist;
410 const InteractionList &ilc = ilist[F_CONSTR];
411 const InteractionList &ils = ilist[F_SETTLE];
412 for (ftype = 0; ftype < F_NRE; ftype++)
414 if (!(ftype == F_BONDS || ftype == F_G96BONDS || ftype == F_HARMONIC))
419 const InteractionList &ilb = ilist[ftype];
420 for (i = 0; i < ilb.size(); i += 3)
422 fc = ip[ilb.iatoms[i]].harmonic.krA;
423 re = ip[ilb.iatoms[i]].harmonic.rA;
424 if (ftype == F_G96BONDS)
426 /* Convert squared sqaure fc to harmonic fc */
429 a1 = ilb.iatoms[i+1];
430 a2 = ilb.iatoms[i+2];
433 if (fc > 0 && m1 > 0 && m2 > 0)
435 period2 = twopi2*m1*m2/((m1 + m2)*fc);
439 period2 = GMX_FLOAT_MAX;
443 fprintf(debug, "fc %g m1 %g m2 %g period %g\n",
444 fc, m1, m2, std::sqrt(period2));
446 if (period2 < limit2)
449 for (j = 0; j < ilc.size(); j += 3)
451 if ((ilc.iatoms[j+1] == a1 && ilc.iatoms[j+2] == a2) ||
452 (ilc.iatoms[j+1] == a2 && ilc.iatoms[j+2] == a1))
457 for (j = 0; j < ils.size(); j += 4)
459 if ((a1 == ils.iatoms[j+1] || a1 == ils.iatoms[j+2] || a1 == ils.iatoms[j+3]) &&
460 (a2 == ils.iatoms[j+1] || a2 == ils.iatoms[j+2] || a2 == ils.iatoms[j+3]))
466 (w_moltype == nullptr || period2 < w_period2))
468 w_moltype = &moltype;
478 if (w_moltype != nullptr)
480 bWarn = (w_period2 < gmx::square(min_steps_warn*dt));
481 /* A check that would recognize most water models */
482 bWater = ((*w_moltype->atoms.atomname[0])[0] == 'O' &&
483 w_moltype->atoms.nr <= 5);
484 sprintf(warn_buf, "The bond in molecule-type %s between atoms %d %s and %d %s has an estimated oscillational period of %.1e ps, which is less than %d times the time step of %.1e ps.\n"
487 w_a1+1, *w_moltype->atoms.atomname[w_a1],
488 w_a2+1, *w_moltype->atoms.atomname[w_a2],
489 std::sqrt(w_period2), bWarn ? min_steps_warn : min_steps_note, dt,
491 "Maybe you asked for fexible water." :
492 "Maybe you forgot to change the constraints mdp option.");
495 warning(wi, warn_buf);
499 warning_note(wi, warn_buf);
504 static void check_vel(gmx_mtop_t *mtop, rvec v[])
506 for (const AtomProxy atomP : AtomRange(*mtop))
508 const t_atom &local = atomP.atom();
509 int i = atomP.globalAtomNumber();
510 if (local.ptype == eptShell ||
511 local.ptype == eptBond ||
512 local.ptype == eptVSite)
519 static void check_shells_inputrec(gmx_mtop_t *mtop,
524 char warn_buf[STRLEN];
526 for (const AtomProxy atomP : AtomRange(*mtop))
528 const t_atom &local = atomP.atom();
529 if (local.ptype == eptShell ||
530 local.ptype == eptBond)
535 if ((nshells > 0) && (ir->nstcalcenergy != 1))
537 set_warning_line(wi, "unknown", -1);
538 snprintf(warn_buf, STRLEN,
539 "There are %d shells, changing nstcalcenergy from %d to 1",
540 nshells, ir->nstcalcenergy);
541 ir->nstcalcenergy = 1;
542 warning(wi, warn_buf);
546 /* TODO Decide whether this function can be consolidated with
547 * gmx_mtop_ftype_count */
548 static int nint_ftype(gmx_mtop_t *mtop, gmx::ArrayRef<const MoleculeInformation> mi, int ftype)
551 for (const gmx_molblock_t &molb : mtop->molblock)
553 nint += molb.nmol*mi[molb.type].interactions[ftype].size();
559 /* This routine reorders the molecule type array
560 * in the order of use in the molblocks,
561 * unused molecule types are deleted.
563 static void renumber_moltypes(gmx_mtop_t *sys,
564 std::vector<MoleculeInformation> *molinfo)
567 std::vector<int> order;
568 for (gmx_molblock_t &molblock : sys->molblock)
570 const auto found = std::find_if(order.begin(), order.end(),
571 [&molblock](const auto &entry)
572 { return molblock.type == entry; });
573 if (found == order.end())
575 /* This type did not occur yet, add it */
576 order.push_back(molblock.type);
577 molblock.type = order.size() - 1;
581 molblock.type = std::distance(order.begin(), found);
585 /* We still need to reorder the molinfo structs */
586 std::vector<MoleculeInformation> minew(order.size());
588 for (auto &mi : *molinfo)
590 const auto found = std::find(order.begin(), order.end(), index);
591 if (found != order.end())
593 int pos = std::distance(order.begin(), found);
598 // Need to manually clean up memory ....
607 static void molinfo2mtop(gmx::ArrayRef<const MoleculeInformation> mi, gmx_mtop_t *mtop)
609 mtop->moltype.resize(mi.size());
611 for (const auto &mol : mi)
613 gmx_moltype_t &molt = mtop->moltype[pos];
614 molt.name = mol.name;
615 molt.atoms = mol.atoms;
616 /* ilists are copied later */
618 molt.excls = mol.excls;
624 new_status(const char *topfile, const char *topppfile, const char *confin,
625 t_gromppopts *opts, t_inputrec *ir, gmx_bool bZero,
626 bool bGenVel, bool bVerbose, t_state *state,
627 PreprocessingAtomTypes *atypes, gmx_mtop_t *sys,
628 std::vector<MoleculeInformation> *mi,
629 std::unique_ptr<MoleculeInformation> *intermolecular_interactions,
630 gmx::ArrayRef<InteractionsOfType> interactions,
631 int *comb, double *reppow, real *fudgeQQ,
635 std::vector<gmx_molblock_t> molblock;
637 bool ffParametrizedWithHBondConstraints;
639 char warn_buf[STRLEN];
641 /* TOPOLOGY processing */
642 sys->name = do_top(bVerbose, topfile, topppfile, opts, bZero, &(sys->symtab),
643 interactions, comb, reppow, fudgeQQ,
644 atypes, mi, intermolecular_interactions,
647 &ffParametrizedWithHBondConstraints,
650 sys->molblock.clear();
653 for (const gmx_molblock_t &molb : molblock)
655 if (!sys->molblock.empty() &&
656 molb.type == sys->molblock.back().type)
658 /* Merge consecutive blocks with the same molecule type */
659 sys->molblock.back().nmol += molb.nmol;
661 else if (molb.nmol > 0)
663 /* Add a new molblock to the topology */
664 sys->molblock.push_back(molb);
666 sys->natoms += molb.nmol*(*mi)[sys->molblock.back().type].atoms.nr;
668 if (sys->molblock.empty())
670 gmx_fatal(FARGS, "No molecules were defined in the system");
673 renumber_moltypes(sys, mi);
677 convert_harmonics(*mi, atypes);
680 if (ir->eDisre == edrNone)
682 i = rm_interactions(F_DISRES, *mi);
685 set_warning_line(wi, "unknown", -1);
686 sprintf(warn_buf, "disre = no, removed %d distance restraints", i);
687 warning_note(wi, warn_buf);
692 i = rm_interactions(F_ORIRES, *mi);
695 set_warning_line(wi, "unknown", -1);
696 sprintf(warn_buf, "orire = no, removed %d orientation restraints", i);
697 warning_note(wi, warn_buf);
701 /* Copy structures from msys to sys */
702 molinfo2mtop(*mi, sys);
704 gmx_mtop_finalize(sys);
706 /* Discourage using the, previously common, setup of applying constraints
707 * to all bonds with force fields that have been parametrized with H-bond
708 * constraints only. Do not print note with large timesteps or vsites.
710 if (opts->nshake == eshALLBONDS &&
711 ffParametrizedWithHBondConstraints &&
712 ir->delta_t < 0.0026 &&
713 gmx_mtop_ftype_count(sys, F_VSITE3FD) == 0)
715 set_warning_line(wi, "unknown", -1);
716 warning_note(wi, "You are using constraints on all bonds, whereas the forcefield "
717 "has been parametrized only with constraints involving hydrogen atoms. "
718 "We suggest using constraints = h-bonds instead, this will also improve performance.");
721 /* COORDINATE file processing */
724 fprintf(stderr, "processing coordinates...\n");
731 read_tps_conf(confin, conftop, nullptr, &x, &v, state->box, FALSE);
732 state->natoms = conftop->atoms.nr;
733 if (state->natoms != sys->natoms)
735 gmx_fatal(FARGS, "number of coordinates in coordinate file (%s, %d)\n"
736 " does not match topology (%s, %d)",
737 confin, state->natoms, topfile, sys->natoms);
739 /* It would be nice to get rid of the copies below, but we don't know
740 * a priori if the number of atoms in confin matches what we expect.
742 state->flags |= (1 << estX);
745 state->flags |= (1 << estV);
747 state_change_natoms(state, state->natoms);
748 std::copy(x, x+state->natoms, state->x.data());
752 std::copy(v, v+state->natoms, state->v.data());
755 /* This call fixes the box shape for runs with pressure scaling */
756 set_box_rel(ir, state);
758 nmismatch = check_atom_names(topfile, confin, sys, &conftop->atoms);
764 sprintf(buf, "%d non-matching atom name%s\n"
765 "atom names from %s will be used\n"
766 "atom names from %s will be ignored\n",
767 nmismatch, (nmismatch == 1) ? "" : "s", topfile, confin);
771 /* Do more checks, mostly related to constraints */
774 fprintf(stderr, "double-checking input for internal consistency...\n");
777 bool bHasNormalConstraints = 0 < (nint_ftype(sys, *mi, F_CONSTR) +
778 nint_ftype(sys, *mi, F_CONSTRNC));
779 bool bHasAnyConstraints = bHasNormalConstraints || 0 < nint_ftype(sys, *mi, F_SETTLE);
780 double_check(ir, state->box,
781 bHasNormalConstraints,
790 snew(mass, state->natoms);
791 for (const AtomProxy atomP : AtomRange(*sys))
793 const t_atom &local = atomP.atom();
794 int i = atomP.globalAtomNumber();
798 if (opts->seed == -1)
800 opts->seed = static_cast<int>(gmx::makeRandomSeed());
801 fprintf(stderr, "Setting gen_seed to %d\n", opts->seed);
803 state->flags |= (1 << estV);
804 maxwell_speed(opts->tempi, opts->seed, sys, state->v.rvec_array());
806 stop_cm(stdout, state->natoms, mass, state->x.rvec_array(), state->v.rvec_array());
811 static void copy_state(const char *slog, t_trxframe *fr,
812 bool bReadVel, t_state *state,
815 if (fr->not_ok & FRAME_NOT_OK)
817 gmx_fatal(FARGS, "Can not start from an incomplete frame");
821 gmx_fatal(FARGS, "Did not find a frame with coordinates in file %s",
825 std::copy(fr->x, fr->x+state->natoms, state->x.data());
830 gmx_incons("Trajecory frame unexpectedly does not contain velocities");
832 std::copy(fr->v, fr->v+state->natoms, state->v.data());
836 copy_mat(fr->box, state->box);
839 *use_time = fr->time;
842 static void cont_status(const char *slog, const char *ener,
843 bool bNeedVel, bool bGenVel, real fr_time,
844 t_inputrec *ir, t_state *state,
846 const gmx_output_env_t *oenv)
847 /* If fr_time == -1 read the last frame available which is complete */
854 bReadVel = (bNeedVel && !bGenVel);
857 "Reading Coordinates%s and Box size from old trajectory\n",
858 bReadVel ? ", Velocities" : "");
861 fprintf(stderr, "Will read whole trajectory\n");
865 fprintf(stderr, "Will read till time %g\n", fr_time);
871 fprintf(stderr, "Velocities generated: "
872 "ignoring velocities in input trajectory\n");
874 read_first_frame(oenv, &fp, slog, &fr, TRX_NEED_X);
878 read_first_frame(oenv, &fp, slog, &fr, TRX_NEED_X | TRX_NEED_V);
884 "WARNING: Did not find a frame with velocities in file %s,\n"
885 " all velocities will be set to zero!\n\n", slog);
886 for (auto &vi : makeArrayRef(state->v))
891 /* Search for a frame without velocities */
893 read_first_frame(oenv, &fp, slog, &fr, TRX_NEED_X);
897 state->natoms = fr.natoms;
899 if (sys->natoms != state->natoms)
901 gmx_fatal(FARGS, "Number of atoms in Topology "
902 "is not the same as in Trajectory");
904 copy_state(slog, &fr, bReadVel, state, &use_time);
906 /* Find the appropriate frame */
907 while ((fr_time == -1 || fr.time < fr_time) &&
908 read_next_frame(oenv, fp, &fr))
910 copy_state(slog, &fr, bReadVel, state, &use_time);
915 /* Set the relative box lengths for preserving the box shape.
916 * Note that this call can lead to differences in the last bit
917 * with respect to using gmx convert-tpr to create a [REF].tpx[ref] file.
919 set_box_rel(ir, state);
921 fprintf(stderr, "Using frame at t = %g ps\n", use_time);
922 fprintf(stderr, "Starting time for run is %g ps\n", ir->init_t);
924 if ((ir->epc != epcNO || ir->etc == etcNOSEHOOVER) && ener)
926 get_enx_state(ener, use_time, sys->groups, ir, state);
927 preserve_box_shape(ir, state->box_rel, state->boxv);
931 static void read_posres(gmx_mtop_t *mtop,
932 gmx::ArrayRef<const MoleculeInformation> molinfo,
935 int rc_scaling, int ePBC,
945 int natoms, npbcdim = 0;
946 char warn_buf[STRLEN];
951 read_tps_conf(fn, top, nullptr, &x, &v, box, FALSE);
952 natoms = top->atoms.nr;
955 if (natoms != mtop->natoms)
957 sprintf(warn_buf, "The number of atoms in %s (%d) does not match the number of atoms in the topology (%d). Will assume that the first %d atoms in the topology and %s match.", fn, natoms, mtop->natoms, std::min(mtop->natoms, natoms), fn);
958 warning(wi, warn_buf);
961 npbcdim = ePBC2npbcdim(ePBC);
963 if (rc_scaling != erscNO)
965 copy_mat(box, invbox);
966 for (int j = npbcdim; j < DIM; j++)
968 clear_rvec(invbox[j]);
971 gmx::invertBoxMatrix(invbox, invbox);
974 /* Copy the reference coordinates to mtop */
978 snew(hadAtom, natoms);
979 for (gmx_molblock_t &molb : mtop->molblock)
981 nat_molb = molb.nmol*mtop->moltype[molb.type].atoms.nr;
982 const InteractionsOfType *pr = &(molinfo[molb.type].interactions[F_POSRES]);
983 const InteractionsOfType *prfb = &(molinfo[molb.type].interactions[F_FBPOSRES]);
984 if (pr->size() > 0 || prfb->size() > 0)
986 atom = mtop->moltype[molb.type].atoms.atom;
987 for (const auto &restraint : pr->interactionTypes)
989 int ai = restraint.ai();
992 gmx_fatal(FARGS, "Position restraint atom index (%d) in moltype '%s' is larger than number of atoms in %s (%d).\n",
993 ai+1, *molinfo[molb.type].name, fn, natoms);
996 if (rc_scaling == erscCOM)
998 /* Determine the center of mass of the posres reference coordinates */
999 for (int j = 0; j < npbcdim; j++)
1001 sum[j] += atom[ai].m*x[a+ai][j];
1003 totmass += atom[ai].m;
1006 /* Same for flat-bottomed posres, but do not count an atom twice for COM */
1007 for (const auto &restraint : prfb->interactionTypes)
1009 int ai = restraint.ai();
1012 gmx_fatal(FARGS, "Position restraint atom index (%d) in moltype '%s' is larger than number of atoms in %s (%d).\n",
1013 ai+1, *molinfo[molb.type].name, fn, natoms);
1015 if (rc_scaling == erscCOM && !hadAtom[ai])
1017 /* Determine the center of mass of the posres reference coordinates */
1018 for (int j = 0; j < npbcdim; j++)
1020 sum[j] += atom[ai].m*x[a+ai][j];
1022 totmass += atom[ai].m;
1027 molb.posres_xA.resize(nat_molb);
1028 for (int i = 0; i < nat_molb; i++)
1030 copy_rvec(x[a+i], molb.posres_xA[i]);
1035 molb.posres_xB.resize(nat_molb);
1036 for (int i = 0; i < nat_molb; i++)
1038 copy_rvec(x[a+i], molb.posres_xB[i]);
1044 if (rc_scaling == erscCOM)
1048 gmx_fatal(FARGS, "The total mass of the position restraint atoms is 0");
1050 for (int j = 0; j < npbcdim; j++)
1052 com[j] = sum[j]/totmass;
1054 fprintf(stderr, "The center of mass of the position restraint coord's is %6.3f %6.3f %6.3f\n", com[XX], com[YY], com[ZZ]);
1057 if (rc_scaling != erscNO)
1059 GMX_ASSERT(npbcdim <= DIM, "Only DIM dimensions can have PBC");
1061 for (gmx_molblock_t &molb : mtop->molblock)
1063 nat_molb = molb.nmol*mtop->moltype[molb.type].atoms.nr;
1064 if (!molb.posres_xA.empty() || !molb.posres_xB.empty())
1066 std::vector<gmx::RVec> &xp = (!bTopB ? molb.posres_xA : molb.posres_xB);
1067 for (int i = 0; i < nat_molb; i++)
1069 for (int j = 0; j < npbcdim; j++)
1071 if (rc_scaling == erscALL)
1073 /* Convert from Cartesian to crystal coordinates */
1074 xp[i][j] *= invbox[j][j];
1075 for (int k = j+1; k < npbcdim; k++)
1077 xp[i][j] += invbox[k][j]*xp[i][k];
1080 else if (rc_scaling == erscCOM)
1082 /* Subtract the center of mass */
1090 if (rc_scaling == erscCOM)
1092 /* Convert the COM from Cartesian to crystal coordinates */
1093 for (int j = 0; j < npbcdim; j++)
1095 com[j] *= invbox[j][j];
1096 for (int k = j+1; k < npbcdim; k++)
1098 com[j] += invbox[k][j]*com[k];
1109 static void gen_posres(gmx_mtop_t *mtop,
1110 gmx::ArrayRef<const MoleculeInformation> mi,
1111 const char *fnA, const char *fnB,
1112 int rc_scaling, int ePBC,
1113 rvec com, rvec comB,
1116 read_posres(mtop, mi, FALSE, fnA, rc_scaling, ePBC, com, wi);
1117 /* It is safer to simply read the b-state posres rather than trying
1118 * to be smart and copy the positions.
1120 read_posres(mtop, mi, TRUE, fnB, rc_scaling, ePBC, comB, wi);
1123 static void set_wall_atomtype(PreprocessingAtomTypes *at, t_gromppopts *opts,
1124 t_inputrec *ir, warninp *wi)
1127 char warn_buf[STRLEN];
1131 fprintf(stderr, "Searching the wall atom type(s)\n");
1133 for (i = 0; i < ir->nwall; i++)
1135 ir->wall_atomtype[i] = at->atomTypeFromName(opts->wall_atomtype[i]);
1136 if (ir->wall_atomtype[i] == NOTSET)
1138 sprintf(warn_buf, "Specified wall atom type %s is not defined", opts->wall_atomtype[i]);
1139 warning_error(wi, warn_buf);
1144 static int nrdf_internal(const t_atoms *atoms)
1149 for (i = 0; i < atoms->nr; i++)
1151 /* Vsite ptype might not be set here yet, so also check the mass */
1152 if ((atoms->atom[i].ptype == eptAtom ||
1153 atoms->atom[i].ptype == eptNucleus)
1154 && atoms->atom[i].m > 0)
1161 case 0: nrdf = 0; break;
1162 case 1: nrdf = 0; break;
1163 case 2: nrdf = 1; break;
1164 default: nrdf = nmass*3 - 6; break;
1171 spline1d( double dx,
1183 for (i = 1; i < n-1; i++)
1185 p = 0.5*y2[i-1]+2.0;
1187 q = (y[i+1]-2.0*y[i]+y[i-1])/dx;
1188 u[i] = (3.0*q/dx-0.5*u[i-1])/p;
1193 for (i = n-2; i >= 0; i--)
1195 y2[i] = y2[i]*y2[i+1]+u[i];
1201 interpolate1d( double xmin,
1212 ix = static_cast<int>((x-xmin)/dx);
1214 a = (xmin+(ix+1)*dx-x)/dx;
1215 b = (x-xmin-ix*dx)/dx;
1217 *y = a*ya[ix]+b*ya[ix+1]+((a*a*a-a)*y2a[ix]+(b*b*b-b)*y2a[ix+1])*(dx*dx)/6.0;
1218 *y1 = (ya[ix+1]-ya[ix])/dx-(3.0*a*a-1.0)/6.0*dx*y2a[ix]+(3.0*b*b-1.0)/6.0*dx*y2a[ix+1];
1223 setup_cmap (int grid_spacing,
1225 gmx::ArrayRef<const real> grid,
1226 gmx_cmap_t * cmap_grid)
1228 int i, j, k, ii, jj, kk, idx;
1230 double dx, xmin, v, v1, v2, v12;
1233 std::vector<double> tmp_u(2*grid_spacing, 0.0);
1234 std::vector<double> tmp_u2(2*grid_spacing, 0.0);
1235 std::vector<double> tmp_yy(2*grid_spacing, 0.0);
1236 std::vector<double> tmp_y1(2*grid_spacing, 0.0);
1237 std::vector<double> tmp_t2(2*grid_spacing*2*grid_spacing, 0.0);
1238 std::vector<double> tmp_grid(2*grid_spacing*2*grid_spacing, 0.0);
1240 dx = 360.0/grid_spacing;
1241 xmin = -180.0-dx*grid_spacing/2;
1243 for (kk = 0; kk < nc; kk++)
1245 /* Compute an offset depending on which cmap we are using
1246 * Offset will be the map number multiplied with the
1247 * grid_spacing * grid_spacing * 2
1249 offset = kk * grid_spacing * grid_spacing * 2;
1251 for (i = 0; i < 2*grid_spacing; i++)
1253 ii = (i+grid_spacing-grid_spacing/2)%grid_spacing;
1255 for (j = 0; j < 2*grid_spacing; j++)
1257 jj = (j+grid_spacing-grid_spacing/2)%grid_spacing;
1258 tmp_grid[i*grid_spacing*2+j] = grid[offset+ii*grid_spacing+jj];
1262 for (i = 0; i < 2*grid_spacing; i++)
1264 spline1d(dx, &(tmp_grid[2*grid_spacing*i]), 2*grid_spacing, tmp_u.data(), &(tmp_t2[2*grid_spacing*i]));
1267 for (i = grid_spacing/2; i < grid_spacing+grid_spacing/2; i++)
1269 ii = i-grid_spacing/2;
1272 for (j = grid_spacing/2; j < grid_spacing+grid_spacing/2; j++)
1274 jj = j-grid_spacing/2;
1277 for (k = 0; k < 2*grid_spacing; k++)
1279 interpolate1d(xmin, dx, &(tmp_grid[2*grid_spacing*k]),
1280 &(tmp_t2[2*grid_spacing*k]), psi, &tmp_yy[k], &tmp_y1[k]);
1283 spline1d(dx, tmp_yy.data(), 2*grid_spacing, tmp_u.data(), tmp_u2.data());
1284 interpolate1d(xmin, dx, tmp_yy.data(), tmp_u2.data(), phi, &v, &v1);
1285 spline1d(dx, tmp_y1.data(), 2*grid_spacing, tmp_u.data(), tmp_u2.data());
1286 interpolate1d(xmin, dx, tmp_y1.data(), tmp_u2.data(), phi, &v2, &v12);
1288 idx = ii*grid_spacing+jj;
1289 cmap_grid->cmapdata[kk].cmap[idx*4] = grid[offset+ii*grid_spacing+jj];
1290 cmap_grid->cmapdata[kk].cmap[idx*4+1] = v1;
1291 cmap_grid->cmapdata[kk].cmap[idx*4+2] = v2;
1292 cmap_grid->cmapdata[kk].cmap[idx*4+3] = v12;
1298 static void init_cmap_grid(gmx_cmap_t *cmap_grid, int ngrid, int grid_spacing)
1302 cmap_grid->grid_spacing = grid_spacing;
1303 nelem = cmap_grid->grid_spacing*cmap_grid->grid_spacing;
1305 cmap_grid->cmapdata.resize(ngrid);
1307 for (i = 0; i < ngrid; i++)
1309 cmap_grid->cmapdata[i].cmap.resize(4*nelem);
1314 static int count_constraints(const gmx_mtop_t *mtop,
1315 gmx::ArrayRef<const MoleculeInformation> mi,
1318 int count, count_mol;
1322 for (const gmx_molblock_t &molb : mtop->molblock)
1325 gmx::ArrayRef<const InteractionsOfType> interactions = mi[molb.type].interactions;
1327 for (int i = 0; i < F_NRE; i++)
1331 count_mol += 3*interactions[i].size();
1333 else if (interaction_function[i].flags & IF_CONSTRAINT)
1335 count_mol += interactions[i].size();
1339 if (count_mol > nrdf_internal(&mi[molb.type].atoms))
1342 "Molecule type '%s' has %d constraints.\n"
1343 "For stability and efficiency there should not be more constraints than internal number of degrees of freedom: %d.\n",
1344 *mi[molb.type].name, count_mol,
1345 nrdf_internal(&mi[molb.type].atoms));
1348 count += molb.nmol*count_mol;
1354 static real calc_temp(const gmx_mtop_t *mtop,
1355 const t_inputrec *ir,
1359 for (const AtomProxy atomP : AtomRange(*mtop))
1361 const t_atom &local = atomP.atom();
1362 int i = atomP.globalAtomNumber();
1363 sum_mv2 += local.m*norm2(v[i]);
1367 for (int g = 0; g < ir->opts.ngtc; g++)
1369 nrdf += ir->opts.nrdf[g];
1372 return sum_mv2/(nrdf*BOLTZ);
1375 static real get_max_reference_temp(const t_inputrec *ir,
1384 for (i = 0; i < ir->opts.ngtc; i++)
1386 if (ir->opts.tau_t[i] < 0)
1392 ref_t = std::max(ref_t, ir->opts.ref_t[i]);
1400 sprintf(buf, "Some temperature coupling groups do not use temperature coupling. We will assume their temperature is not more than %.3f K. If their temperature is higher, the energy error and the Verlet buffer might be underestimated.",
1408 /* Checks if there are unbound atoms in moleculetype molt.
1409 * Prints a note for each unbound atoms and a warning if any is present.
1411 static void checkForUnboundAtoms(const gmx_moltype_t *molt,
1415 const t_atoms *atoms = &molt->atoms;
1419 /* Only one atom, there can't be unbound atoms */
1423 std::vector<int> count(atoms->nr, 0);
1425 for (int ftype = 0; ftype < F_NRE; ftype++)
1427 if (((interaction_function[ftype].flags & IF_BOND) && ftype != F_CONNBONDS) ||
1428 (interaction_function[ftype].flags & IF_CONSTRAINT) ||
1431 const InteractionList &il = molt->ilist[ftype];
1432 const int nral = NRAL(ftype);
1434 for (int i = 0; i < il.size(); i += 1 + nral)
1436 for (int j = 0; j < nral; j++)
1438 count[il.iatoms[i + 1 + j]]++;
1444 int numDanglingAtoms = 0;
1445 for (int a = 0; a < atoms->nr; a++)
1447 if (atoms->atom[a].ptype != eptVSite &&
1452 fprintf(stderr, "\nAtom %d '%s' in moleculetype '%s' is not bound by a potential or constraint to any other atom in the same moleculetype.\n",
1453 a + 1, *atoms->atomname[a], *molt->name);
1459 if (numDanglingAtoms > 0)
1462 sprintf(buf, "In moleculetype '%s' %d atoms are not bound by a potential or constraint to any other atom in the same moleculetype. Although technically this might not cause issues in a simulation, this often means that the user forgot to add a bond/potential/constraint or put multiple molecules in the same moleculetype definition by mistake. Run with -v to get information for each atom.",
1463 *molt->name, numDanglingAtoms);
1464 warning_note(wi, buf);
1468 /* Checks all moleculetypes for unbound atoms */
1469 static void checkForUnboundAtoms(const gmx_mtop_t *mtop,
1473 for (const gmx_moltype_t &molt : mtop->moltype)
1475 checkForUnboundAtoms(&molt, bVerbose, wi);
1479 /*! \brief Checks if there are decoupled modes in moleculetype \p molt.
1481 * The specific decoupled modes this routine check for are angle modes
1482 * where the two bonds are constrained and the atoms a both ends are only
1483 * involved in a single constraint; the mass of the two atoms needs to
1484 * differ by more than \p massFactorThreshold.
1486 static bool haveDecoupledModeInMol(const gmx_moltype_t &molt,
1487 gmx::ArrayRef<const t_iparams> iparams,
1488 real massFactorThreshold)
1490 if (molt.ilist[F_CONSTR].size() == 0 &&
1491 molt.ilist[F_CONSTRNC].size() == 0)
1496 const t_atom * atom = molt.atoms.atom;
1498 t_blocka atomToConstraints =
1499 gmx::make_at2con(molt, iparams,
1500 gmx::FlexibleConstraintTreatment::Exclude);
1502 bool haveDecoupledMode = false;
1503 for (int ftype = 0; ftype < F_NRE; ftype++)
1505 if (interaction_function[ftype].flags & IF_ATYPE)
1507 const int nral = NRAL(ftype);
1508 const InteractionList &il = molt.ilist[ftype];
1509 for (int i = 0; i < il.size(); i += 1 + nral)
1511 /* Here we check for the mass difference between the atoms
1512 * at both ends of the angle, that the atoms at the ends have
1513 * 1 contraint and the atom in the middle at least 3; we check
1514 * that the 3 atoms are linked by constraints below.
1515 * We check for at least three constraints for the middle atom,
1516 * since with only the two bonds in the angle, we have 3-atom
1517 * molecule, which has much more thermal exhange in this single
1518 * angle mode than molecules with more atoms.
1519 * Note that this check also catches molecules where more atoms
1520 * are connected to one or more atoms in the angle, but by
1521 * bond potentials instead of angles. But such cases will not
1522 * occur in "normal" molecules and it doesn't hurt running
1523 * those with higher accuracy settings as well.
1525 int a0 = il.iatoms[1 + i];
1526 int a1 = il.iatoms[1 + i + 1];
1527 int a2 = il.iatoms[1 + i + 2];
1528 if ((atom[a0].m > atom[a2].m*massFactorThreshold ||
1529 atom[a2].m > atom[a0].m*massFactorThreshold) &&
1530 atomToConstraints.index[a0 + 1] - atomToConstraints.index[a0] == 1 &&
1531 atomToConstraints.index[a2 + 1] - atomToConstraints.index[a2] == 1 &&
1532 atomToConstraints.index[a1 + 1] - atomToConstraints.index[a1] >= 3)
1534 int constraint0 = atomToConstraints.a[atomToConstraints.index[a0]];
1535 int constraint2 = atomToConstraints.a[atomToConstraints.index[a2]];
1537 bool foundAtom0 = false;
1538 bool foundAtom2 = false;
1539 for (int conIndex = atomToConstraints.index[a1]; conIndex < atomToConstraints.index[a1 + 1]; conIndex++)
1541 if (atomToConstraints.a[conIndex] == constraint0)
1545 if (atomToConstraints.a[conIndex] == constraint2)
1550 if (foundAtom0 && foundAtom2)
1552 haveDecoupledMode = true;
1559 done_blocka(&atomToConstraints);
1561 return haveDecoupledMode;
1564 /*! \brief Checks if the Verlet buffer and constraint accuracy is sufficient for decoupled dynamic modes.
1566 * When decoupled modes are present and the accuray in insufficient,
1567 * this routine issues a warning if the accuracy is insufficient.
1569 static void checkDecoupledModeAccuracy(const gmx_mtop_t *mtop,
1570 const t_inputrec *ir,
1573 /* We only have issues with decoupled modes with normal MD.
1574 * With stochastic dynamics equipartitioning is enforced strongly.
1581 /* When atoms of very different mass are involved in an angle potential
1582 * and both bonds in the angle are constrained, the dynamic modes in such
1583 * angles have very different periods and significant energy exchange
1584 * takes several nanoseconds. Thus even a small amount of error in
1585 * different algorithms can lead to violation of equipartitioning.
1586 * The parameters below are mainly based on an all-atom chloroform model
1587 * with all bonds constrained. Equipartitioning is off by more than 1%
1588 * (need to run 5-10 ns) when the difference in mass between H and Cl
1589 * is more than a factor 13 and the accuracy is less than the thresholds
1590 * given below. This has been verified on some other molecules.
1592 * Note that the buffer and shake parameters have unit length and
1593 * energy/time, respectively, so they will "only" work correctly
1594 * for atomistic force fields using MD units.
1596 const real massFactorThreshold = 13.0;
1597 const real bufferToleranceThreshold = 1e-4;
1598 const int lincsIterationThreshold = 2;
1599 const int lincsOrderThreshold = 4;
1600 const real shakeToleranceThreshold = 0.005*ir->delta_t;
1602 bool lincsWithSufficientTolerance = (ir->eConstrAlg == econtLINCS && ir->nLincsIter >= lincsIterationThreshold && ir->nProjOrder >= lincsOrderThreshold);
1603 bool shakeWithSufficientTolerance = (ir->eConstrAlg == econtSHAKE && ir->shake_tol <= 1.1*shakeToleranceThreshold);
1604 if (ir->cutoff_scheme == ecutsVERLET &&
1605 ir->verletbuf_tol <= 1.1*bufferToleranceThreshold &&
1606 (lincsWithSufficientTolerance || shakeWithSufficientTolerance))
1611 bool haveDecoupledMode = false;
1612 for (const gmx_moltype_t &molt : mtop->moltype)
1614 if (haveDecoupledModeInMol(molt, mtop->ffparams.iparams,
1615 massFactorThreshold))
1617 haveDecoupledMode = true;
1621 if (haveDecoupledMode)
1623 std::string message = gmx::formatString(
1624 "There are atoms at both ends of an angle, connected by constraints "
1625 "and with masses that differ by more than a factor of %g. This means "
1626 "that there are likely dynamic modes that are only very weakly coupled.",
1627 massFactorThreshold);
1628 if (ir->cutoff_scheme == ecutsVERLET)
1630 message += gmx::formatString(
1631 " To ensure good equipartitioning, you need to either not use "
1632 "constraints on all bonds (but, if possible, only on bonds involving "
1633 "hydrogens) or use integrator = %s or decrease one or more tolerances: "
1634 "verlet-buffer-tolerance <= %g, LINCS iterations >= %d, LINCS order "
1635 ">= %d or SHAKE tolerance <= %g",
1637 bufferToleranceThreshold,
1638 lincsIterationThreshold, lincsOrderThreshold,
1639 shakeToleranceThreshold);
1643 message += gmx::formatString(
1644 " To ensure good equipartitioning, we suggest to switch to the %s "
1645 "cutoff-scheme, since that allows for better control over the Verlet "
1646 "buffer size and thus over the energy drift.",
1647 ecutscheme_names[ecutsVERLET]);
1649 warning(wi, message);
1653 static void set_verlet_buffer(const gmx_mtop_t *mtop,
1659 char warn_buf[STRLEN];
1661 printf("Determining Verlet buffer for a tolerance of %g kJ/mol/ps at %g K\n", ir->verletbuf_tol, buffer_temp);
1663 /* Calculate the buffer size for simple atom vs atoms list */
1664 VerletbufListSetup listSetup1x1;
1665 listSetup1x1.cluster_size_i = 1;
1666 listSetup1x1.cluster_size_j = 1;
1667 const real rlist_1x1 =
1668 calcVerletBufferSize(*mtop, det(box), *ir, ir->nstlist, ir->nstlist - 1,
1669 buffer_temp, listSetup1x1);
1671 /* Set the pair-list buffer size in ir */
1672 VerletbufListSetup listSetup4x4 =
1673 verletbufGetSafeListSetup(ListSetupType::CpuNoSimd);
1675 calcVerletBufferSize(*mtop, det(box), *ir, ir->nstlist, ir->nstlist - 1,
1676 buffer_temp, listSetup4x4);
1678 const int n_nonlin_vsite = countNonlinearVsites(*mtop);
1679 if (n_nonlin_vsite > 0)
1681 sprintf(warn_buf, "There are %d non-linear virtual site constructions. Their contribution to the energy error is approximated. In most cases this does not affect the error significantly.", n_nonlin_vsite);
1682 warning_note(wi, warn_buf);
1685 printf("Calculated rlist for %dx%d atom pair-list as %.3f nm, buffer size %.3f nm\n",
1686 1, 1, rlist_1x1, rlist_1x1-std::max(ir->rvdw, ir->rcoulomb));
1688 printf("Set rlist, assuming %dx%d atom pair-list, to %.3f nm, buffer size %.3f nm\n",
1689 listSetup4x4.cluster_size_i, listSetup4x4.cluster_size_j,
1690 ir->rlist, ir->rlist-std::max(ir->rvdw, ir->rcoulomb));
1692 printf("Note that mdrun will redetermine rlist based on the actual pair-list setup\n");
1694 if (gmx::square(ir->rlist) >= max_cutoff2(ir->ePBC, box))
1696 gmx_fatal(FARGS, "The pair-list cut-off (%g nm) is longer than half the shortest box vector or longer than the smallest box diagonal element (%g nm). Increase the box size or decrease nstlist or increase verlet-buffer-tolerance.", ir->rlist, std::sqrt(max_cutoff2(ir->ePBC, box)));
1700 int gmx_grompp(int argc, char *argv[])
1702 const char *desc[] = {
1703 "[THISMODULE] (the gromacs preprocessor)",
1704 "reads a molecular topology file, checks the validity of the",
1705 "file, expands the topology from a molecular description to an atomic",
1706 "description. The topology file contains information about",
1707 "molecule types and the number of molecules, the preprocessor",
1708 "copies each molecule as needed. ",
1709 "There is no limitation on the number of molecule types. ",
1710 "Bonds and bond-angles can be converted into constraints, separately",
1711 "for hydrogens and heavy atoms.",
1712 "Then a coordinate file is read and velocities can be generated",
1713 "from a Maxwellian distribution if requested.",
1714 "[THISMODULE] also reads parameters for [gmx-mdrun] ",
1715 "(eg. number of MD steps, time step, cut-off), and others such as",
1716 "NEMD parameters, which are corrected so that the net acceleration",
1718 "Eventually a binary file is produced that can serve as the sole input",
1719 "file for the MD program.[PAR]",
1721 "[THISMODULE] uses the atom names from the topology file. The atom names",
1722 "in the coordinate file (option [TT]-c[tt]) are only read to generate",
1723 "warnings when they do not match the atom names in the topology.",
1724 "Note that the atom names are irrelevant for the simulation as",
1725 "only the atom types are used for generating interaction parameters.[PAR]",
1727 "[THISMODULE] uses a built-in preprocessor to resolve includes, macros, ",
1728 "etc. The preprocessor supports the following keywords::",
1731 " #ifndef VARIABLE",
1734 " #define VARIABLE",
1736 " #include \"filename\"",
1737 " #include <filename>",
1739 "The functioning of these statements in your topology may be modulated by",
1740 "using the following two flags in your [REF].mdp[ref] file::",
1742 " define = -DVARIABLE1 -DVARIABLE2",
1743 " include = -I/home/john/doe",
1745 "For further information a C-programming textbook may help you out.",
1746 "Specifying the [TT]-pp[tt] flag will get the pre-processed",
1747 "topology file written out so that you can verify its contents.[PAR]",
1749 "When using position restraints, a file with restraint coordinates",
1750 "must be supplied with [TT]-r[tt] (can be the same file as supplied",
1751 "for [TT]-c[tt]). For free energy calculations, separate reference",
1752 "coordinates for the B topology can be supplied with [TT]-rb[tt],",
1753 "otherwise they will be equal to those of the A topology.[PAR]",
1755 "Starting coordinates can be read from trajectory with [TT]-t[tt].",
1756 "The last frame with coordinates and velocities will be read,",
1757 "unless the [TT]-time[tt] option is used. Only if this information",
1758 "is absent will the coordinates in the [TT]-c[tt] file be used.",
1759 "Note that these velocities will not be used when [TT]gen_vel = yes[tt]",
1760 "in your [REF].mdp[ref] file. An energy file can be supplied with",
1761 "[TT]-e[tt] to read Nose-Hoover and/or Parrinello-Rahman coupling",
1764 "[THISMODULE] can be used to restart simulations (preserving",
1765 "continuity) by supplying just a checkpoint file with [TT]-t[tt].",
1766 "However, for simply changing the number of run steps to extend",
1767 "a run, using [gmx-convert-tpr] is more convenient than [THISMODULE].",
1768 "You then supply the old checkpoint file directly to [gmx-mdrun]",
1769 "with [TT]-cpi[tt]. If you wish to change the ensemble or things",
1770 "like output frequency, then supplying the checkpoint file to",
1771 "[THISMODULE] with [TT]-t[tt] along with a new [REF].mdp[ref] file",
1772 "with [TT]-f[tt] is the recommended procedure. Actually preserving",
1773 "the ensemble (if possible) still requires passing the checkpoint",
1774 "file to [gmx-mdrun] [TT]-cpi[tt].[PAR]",
1776 "By default, all bonded interactions which have constant energy due to",
1777 "virtual site constructions will be removed. If this constant energy is",
1778 "not zero, this will result in a shift in the total energy. All bonded",
1779 "interactions can be kept by turning off [TT]-rmvsbds[tt]. Additionally,",
1780 "all constraints for distances which will be constant anyway because",
1781 "of virtual site constructions will be removed. If any constraints remain",
1782 "which involve virtual sites, a fatal error will result.[PAR]",
1784 "To verify your run input file, please take note of all warnings",
1785 "on the screen, and correct where necessary. Do also look at the contents",
1786 "of the [TT]mdout.mdp[tt] file; this contains comment lines, as well as",
1787 "the input that [THISMODULE] has read. If in doubt, you can start [THISMODULE]",
1788 "with the [TT]-debug[tt] option which will give you more information",
1789 "in a file called [TT]grompp.log[tt] (along with real debug info). You",
1790 "can see the contents of the run input file with the [gmx-dump]",
1791 "program. [gmx-check] can be used to compare the contents of two",
1792 "run input files.[PAR]",
1794 "The [TT]-maxwarn[tt] option can be used to override warnings printed",
1795 "by [THISMODULE] that otherwise halt output. In some cases, warnings are",
1796 "harmless, but usually they are not. The user is advised to carefully",
1797 "interpret the output messages before attempting to bypass them with",
1801 std::vector<MoleculeInformation> mi;
1802 std::unique_ptr<MoleculeInformation> intermolecular_interactions;
1806 const char *mdparin;
1807 bool bNeedVel, bGenVel;
1808 gmx_bool have_atomnumber;
1809 gmx_output_env_t *oenv;
1810 gmx_bool bVerbose = FALSE;
1812 char warn_buf[STRLEN];
1815 { efMDP, nullptr, nullptr, ffREAD },
1816 { efMDP, "-po", "mdout", ffWRITE },
1817 { efSTX, "-c", nullptr, ffREAD },
1818 { efSTX, "-r", "restraint", ffOPTRD },
1819 { efSTX, "-rb", "restraint", ffOPTRD },
1820 { efNDX, nullptr, nullptr, ffOPTRD },
1821 { efTOP, nullptr, nullptr, ffREAD },
1822 { efTOP, "-pp", "processed", ffOPTWR },
1823 { efTPR, "-o", nullptr, ffWRITE },
1824 { efTRN, "-t", nullptr, ffOPTRD },
1825 { efEDR, "-e", nullptr, ffOPTRD },
1826 /* This group is needed by the VMD viewer as the start configuration for IMD sessions: */
1827 { efGRO, "-imd", "imdgroup", ffOPTWR },
1828 { efTRN, "-ref", "rotref", ffOPTRW }
1830 #define NFILE asize(fnm)
1832 /* Command line options */
1833 gmx_bool bRenum = TRUE;
1834 gmx_bool bRmVSBds = TRUE, bZero = FALSE;
1838 { "-v", FALSE, etBOOL, {&bVerbose},
1839 "Be loud and noisy" },
1840 { "-time", FALSE, etREAL, {&fr_time},
1841 "Take frame at or first after this time." },
1842 { "-rmvsbds", FALSE, etBOOL, {&bRmVSBds},
1843 "Remove constant bonded interactions with virtual sites" },
1844 { "-maxwarn", FALSE, etINT, {&maxwarn},
1845 "Number of allowed warnings during input processing. Not for normal use and may generate unstable systems" },
1846 { "-zero", FALSE, etBOOL, {&bZero},
1847 "Set parameters for bonded interactions without defaults to zero instead of generating an error" },
1848 { "-renum", FALSE, etBOOL, {&bRenum},
1849 "Renumber atomtypes and minimize number of atomtypes" }
1852 /* Parse the command line */
1853 if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa,
1854 asize(desc), desc, 0, nullptr, &oenv))
1859 /* Initiate some variables */
1860 gmx::MDModules mdModules;
1861 t_inputrec irInstance;
1862 t_inputrec *ir = &irInstance;
1864 snew(opts->include, STRLEN);
1865 snew(opts->define, STRLEN);
1867 wi = init_warning(TRUE, maxwarn);
1869 /* PARAMETER file processing */
1870 mdparin = opt2fn("-f", NFILE, fnm);
1871 set_warning_line(wi, mdparin, -1);
1874 get_ir(mdparin, opt2fn("-po", NFILE, fnm), &mdModules, ir, opts, WriteMdpHeader::yes, wi);
1876 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
1880 fprintf(stderr, "checking input for internal consistency...\n");
1882 check_ir(mdparin, ir, opts, wi);
1884 if (ir->ld_seed == -1)
1886 ir->ld_seed = static_cast<int>(gmx::makeRandomSeed());
1887 fprintf(stderr, "Setting the LD random seed to %" PRId64 "\n", ir->ld_seed);
1890 if (ir->expandedvals->lmc_seed == -1)
1892 ir->expandedvals->lmc_seed = static_cast<int>(gmx::makeRandomSeed());
1893 fprintf(stderr, "Setting the lambda MC random seed to %d\n", ir->expandedvals->lmc_seed);
1896 bNeedVel = EI_STATE_VELOCITY(ir->eI);
1897 bGenVel = (bNeedVel && opts->bGenVel);
1898 if (bGenVel && ir->bContinuation)
1901 "Generating velocities is inconsistent with attempting "
1902 "to continue a previous run. Choose only one of "
1903 "gen-vel = yes and continuation = yes.");
1904 warning_error(wi, warn_buf);
1907 std::array<InteractionsOfType, F_NRE> interactions;
1909 PreprocessingAtomTypes atypes;
1912 pr_symtab(debug, 0, "Just opened", &sys.symtab);
1915 const char *fn = ftp2fn(efTOP, NFILE, fnm);
1916 if (!gmx_fexist(fn))
1918 gmx_fatal(FARGS, "%s does not exist", fn);
1922 new_status(fn, opt2fn_null("-pp", NFILE, fnm), opt2fn("-c", NFILE, fnm),
1923 opts, ir, bZero, bGenVel, bVerbose, &state,
1924 &atypes, &sys, &mi, &intermolecular_interactions,
1925 interactions, &comb, &reppow, &fudgeQQ,
1931 pr_symtab(debug, 0, "After new_status", &sys.symtab);
1935 /* set parameters for virtual site construction (not for vsiten) */
1936 for (size_t mt = 0; mt < sys.moltype.size(); mt++)
1939 set_vsites(bVerbose, &sys.moltype[mt].atoms, &atypes, mi[mt].interactions);
1941 /* now throw away all obsolete bonds, angles and dihedrals: */
1942 /* note: constraints are ALWAYS removed */
1945 for (size_t mt = 0; mt < sys.moltype.size(); mt++)
1947 clean_vsite_bondeds(mi[mt].interactions, sys.moltype[mt].atoms.nr, bRmVSBds);
1951 if (ir->cutoff_scheme == ecutsVERLET)
1953 fprintf(stderr, "Removing all charge groups because cutoff-scheme=%s\n",
1954 ecutscheme_names[ir->cutoff_scheme]);
1956 /* Remove all charge groups */
1957 gmx_mtop_remove_chargegroups(&sys);
1960 if ((count_constraints(&sys, mi, wi) != 0) && (ir->eConstrAlg == econtSHAKE))
1962 if (ir->eI == eiCG || ir->eI == eiLBFGS)
1964 sprintf(warn_buf, "Can not do %s with %s, use %s",
1965 EI(ir->eI), econstr_names[econtSHAKE], econstr_names[econtLINCS]);
1966 warning_error(wi, warn_buf);
1968 if (ir->bPeriodicMols)
1970 sprintf(warn_buf, "Can not do periodic molecules with %s, use %s",
1971 econstr_names[econtSHAKE], econstr_names[econtLINCS]);
1972 warning_error(wi, warn_buf);
1976 if (EI_SD (ir->eI) && ir->etc != etcNO)
1978 warning_note(wi, "Temperature coupling is ignored with SD integrators.");
1981 /* If we are doing QM/MM, check that we got the atom numbers */
1982 have_atomnumber = TRUE;
1983 for (int i = 0; i < gmx::ssize(atypes); i++)
1985 have_atomnumber = have_atomnumber && (atypes.atomNumberFromAtomType(i) >= 0);
1987 if (!have_atomnumber && ir->bQMMM)
1991 "It appears as if you are trying to run a QM/MM calculation, but the force\n"
1992 "field you are using does not contain atom numbers fields. This is an\n"
1993 "optional field (introduced in GROMACS 3.3) for general runs, but mandatory\n"
1994 "for QM/MM. The good news is that it is easy to add - put the atom number as\n"
1995 "an integer just before the mass column in ffXXXnb.itp.\n"
1996 "NB: United atoms have the same atom numbers as normal ones.\n\n");
1999 /* Check for errors in the input now, since they might cause problems
2000 * during processing further down.
2002 check_warning_error(wi, FARGS);
2004 if (nint_ftype(&sys, mi, F_POSRES) > 0 ||
2005 nint_ftype(&sys, mi, F_FBPOSRES) > 0)
2007 if (ir->epc == epcPARRINELLORAHMAN || ir->epc == epcMTTK)
2009 sprintf(warn_buf, "You are combining position restraints with %s pressure coupling, which can lead to instabilities. If you really want to combine position restraints with pressure coupling, we suggest to use %s pressure coupling instead.",
2010 EPCOUPLTYPE(ir->epc), EPCOUPLTYPE(epcBERENDSEN));
2011 warning_note(wi, warn_buf);
2014 const char *fn = opt2fn("-r", NFILE, fnm);
2017 if (!gmx_fexist(fn))
2020 "Cannot find position restraint file %s (option -r).\n"
2021 "From GROMACS-2018, you need to specify the position restraint "
2022 "coordinate files explicitly to avoid mistakes, although you can "
2023 "still use the same file as you specify for the -c option.", fn);
2026 if (opt2bSet("-rb", NFILE, fnm))
2028 fnB = opt2fn("-rb", NFILE, fnm);
2029 if (!gmx_fexist(fnB))
2032 "Cannot find B-state position restraint file %s (option -rb).\n"
2033 "From GROMACS-2018, you need to specify the position restraint "
2034 "coordinate files explicitly to avoid mistakes, although you can "
2035 "still use the same file as you specify for the -c option.", fn);
2045 fprintf(stderr, "Reading position restraint coords from %s", fn);
2046 if (strcmp(fn, fnB) == 0)
2048 fprintf(stderr, "\n");
2052 fprintf(stderr, " and %s\n", fnB);
2055 gen_posres(&sys, mi, fn, fnB,
2056 ir->refcoord_scaling, ir->ePBC,
2057 ir->posres_com, ir->posres_comB,
2061 /* If we are using CMAP, setup the pre-interpolation grid */
2062 if (interactions[F_CMAP].ncmap() > 0)
2064 init_cmap_grid(&sys.ffparams.cmap_grid, interactions[F_CMAP].cmapAngles, interactions[F_CMAP].cmakeGridSpacing);
2065 setup_cmap(interactions[F_CMAP].cmakeGridSpacing, interactions[F_CMAP].cmapAngles, interactions[F_CMAP].cmap, &sys.ffparams.cmap_grid);
2068 set_wall_atomtype(&atypes, opts, ir, wi);
2071 atypes.renumberTypes(interactions, &sys, ir->wall_atomtype, bVerbose);
2074 if (ir->implicit_solvent)
2076 gmx_fatal(FARGS, "Implicit solvation is no longer supported");
2079 /* PELA: Copy the atomtype data to the topology atomtype list */
2080 atypes.copyTot_atomtypes(&(sys.atomtypes));
2084 pr_symtab(debug, 0, "After atype.renumberTypes", &sys.symtab);
2089 fprintf(stderr, "converting bonded parameters...\n");
2092 const int ntype = atypes.size();
2093 convertInteractionsOfType(ntype, interactions, mi, intermolecular_interactions.get(),
2094 comb, reppow, fudgeQQ, &sys);
2098 pr_symtab(debug, 0, "After converInteractionsOfType", &sys.symtab);
2101 /* set ptype to VSite for virtual sites */
2102 for (gmx_moltype_t &moltype : sys.moltype)
2104 set_vsites_ptype(FALSE, &moltype);
2108 pr_symtab(debug, 0, "After virtual sites", &sys.symtab);
2110 /* Check velocity for virtual sites and shells */
2113 check_vel(&sys, state.v.rvec_array());
2116 /* check for shells and inpurecs */
2117 check_shells_inputrec(&sys, ir, wi);
2120 check_mol(&sys, wi);
2122 checkForUnboundAtoms(&sys, bVerbose, wi);
2124 for (const gmx_moltype_t &moltype : sys.moltype)
2126 check_cg_sizes(ftp2fn(efTOP, NFILE, fnm), &moltype.cgs, wi);
2129 if (EI_DYNAMICS(ir->eI) && ir->eI != eiBD)
2131 check_bonds_timestep(&sys, ir->delta_t, wi);
2134 checkDecoupledModeAccuracy(&sys, ir, wi);
2136 if (EI_ENERGY_MINIMIZATION(ir->eI) && 0 == ir->nsteps)
2138 warning_note(wi, "Zero-step energy minimization will alter the coordinates before calculating the energy. If you just want the energy of a single point, try zero-step MD (with unconstrained_start = yes). To do multiple single-point energy evaluations of different configurations of the same topology, use mdrun -rerun.");
2141 check_warning_error(wi, FARGS);
2145 fprintf(stderr, "initialising group options...\n");
2147 do_index(mdparin, ftp2fn_null(efNDX, NFILE, fnm),
2151 if (ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0)
2153 if (EI_DYNAMICS(ir->eI) && inputrec2nboundeddim(ir) == 3)
2157 if (EI_MD(ir->eI) && ir->etc == etcNO)
2161 buffer_temp = opts->tempi;
2165 buffer_temp = calc_temp(&sys, ir, state.v.rvec_array());
2167 if (buffer_temp > 0)
2169 sprintf(warn_buf, "NVE simulation: will use the initial temperature of %.3f K for determining the Verlet buffer size", buffer_temp);
2170 warning_note(wi, warn_buf);
2174 sprintf(warn_buf, "NVE simulation with an initial temperature of zero: will use a Verlet buffer of %d%%. Check your energy drift!",
2175 gmx::roundToInt(verlet_buffer_ratio_NVE_T0*100));
2176 warning_note(wi, warn_buf);
2181 buffer_temp = get_max_reference_temp(ir, wi);
2184 if (EI_MD(ir->eI) && ir->etc == etcNO && buffer_temp == 0)
2186 /* NVE with initial T=0: we add a fixed ratio to rlist.
2187 * Since we don't actually use verletbuf_tol, we set it to -1
2188 * so it can't be misused later.
2190 ir->rlist *= 1.0 + verlet_buffer_ratio_NVE_T0;
2191 ir->verletbuf_tol = -1;
2195 /* We warn for NVE simulations with a drift tolerance that
2196 * might result in a 1(.1)% drift over the total run-time.
2197 * Note that we can't warn when nsteps=0, since we don't
2198 * know how many steps the user intends to run.
2200 if (EI_MD(ir->eI) && ir->etc == etcNO && ir->nstlist > 1 &&
2203 const real driftTolerance = 0.01;
2204 /* We use 2 DOF per atom = 2kT pot+kin energy,
2205 * to be on the safe side with constraints.
2207 const real totalEnergyDriftPerAtomPerPicosecond = 2*BOLTZ*buffer_temp/(ir->nsteps*ir->delta_t);
2209 if (ir->verletbuf_tol > 1.1*driftTolerance*totalEnergyDriftPerAtomPerPicosecond)
2211 sprintf(warn_buf, "You are using a Verlet buffer tolerance of %g kJ/mol/ps for an NVE simulation of length %g ps, which can give a final drift of %d%%. For conserving energy to %d%% when using constraints, you might need to set verlet-buffer-tolerance to %.1e.",
2212 ir->verletbuf_tol, ir->nsteps*ir->delta_t,
2213 gmx::roundToInt(ir->verletbuf_tol/totalEnergyDriftPerAtomPerPicosecond*100),
2214 gmx::roundToInt(100*driftTolerance),
2215 driftTolerance*totalEnergyDriftPerAtomPerPicosecond);
2216 warning_note(wi, warn_buf);
2220 set_verlet_buffer(&sys, ir, buffer_temp, state.box, wi);
2225 /* Init the temperature coupling state */
2226 init_gtc_state(&state, ir->opts.ngtc, 0, ir->opts.nhchainlength); /* need to add nnhpres here? */
2230 fprintf(stderr, "Checking consistency between energy and charge groups...\n");
2232 check_eg_vs_cg(&sys);
2236 pr_symtab(debug, 0, "After index", &sys.symtab);
2239 triple_check(mdparin, ir, &sys, wi);
2240 close_symtab(&sys.symtab);
2243 pr_symtab(debug, 0, "After close", &sys.symtab);
2246 /* make exclusions between QM atoms and remove charges if needed */
2249 if (ir->QMMMscheme == eQMMMschemenormal && ir->ns_type == ensSIMPLE)
2251 gmx_fatal(FARGS, "electrostatic embedding only works with grid neighboursearching, use ns-type=grid instead\n");
2255 generate_qmexcl(&sys, ir, wi, GmxQmmmMode::GMX_QMMM_ORIGINAL);
2257 if (ir->QMMMscheme != eQMMMschemeoniom)
2259 std::vector<int> qmmmAtoms = qmmmAtomIndices(*ir, sys);
2260 removeQmmmAtomCharges(&sys, qmmmAtoms);
2264 if (ir->eI == eiMimic)
2266 generate_qmexcl(&sys, ir, wi, GmxQmmmMode::GMX_QMMM_MIMIC);
2269 if (ftp2bSet(efTRN, NFILE, fnm))
2273 fprintf(stderr, "getting data from old trajectory ...\n");
2275 cont_status(ftp2fn(efTRN, NFILE, fnm), ftp2fn_null(efEDR, NFILE, fnm),
2276 bNeedVel, bGenVel, fr_time, ir, &state, &sys, oenv);
2279 if (ir->ePBC == epbcXY && ir->nwall != 2)
2281 clear_rvec(state.box[ZZ]);
2284 if (ir->cutoff_scheme != ecutsVERLET && ir->rlist > 0)
2286 set_warning_line(wi, mdparin, -1);
2287 check_chargegroup_radii(&sys, ir, state.x.rvec_array(), wi);
2290 if (EEL_FULL(ir->coulombtype) || EVDW_PME(ir->vdwtype))
2292 /* Calculate the optimal grid dimensions */
2294 EwaldBoxZScaler boxScaler(*ir);
2295 boxScaler.scaleBox(state.box, scaledBox);
2297 if (ir->nkx > 0 && ir->nky > 0 && ir->nkz > 0)
2299 /* Mark fourier_spacing as not used */
2300 ir->fourier_spacing = 0;
2302 else if (ir->nkx != 0 && ir->nky != 0 && ir->nkz != 0)
2304 set_warning_line(wi, mdparin, -1);
2305 warning_error(wi, "Some of the Fourier grid sizes are set, but all of them need to be set.");
2307 const int minGridSize = minimalPmeGridSize(ir->pme_order);
2308 calcFftGrid(stdout, scaledBox, ir->fourier_spacing, minGridSize,
2309 &(ir->nkx), &(ir->nky), &(ir->nkz));
2310 if (ir->nkx < minGridSize ||
2311 ir->nky < minGridSize ||
2312 ir->nkz < minGridSize)
2314 warning_error(wi, "The PME grid size should be >= 2*(pme-order - 1); either manually increase the grid size or decrease pme-order");
2318 /* MRS: eventually figure out better logic for initializing the fep
2319 values that makes declaring the lambda and declaring the state not
2320 potentially conflict if not handled correctly. */
2321 if (ir->efep != efepNO)
2323 state.fep_state = ir->fepvals->init_fep_state;
2324 for (i = 0; i < efptNR; i++)
2326 /* init_lambda trumps state definitions*/
2327 if (ir->fepvals->init_lambda >= 0)
2329 state.lambda[i] = ir->fepvals->init_lambda;
2333 if (ir->fepvals->all_lambda[i] == nullptr)
2335 gmx_fatal(FARGS, "Values of lambda not set for a free energy calculation!");
2339 state.lambda[i] = ir->fepvals->all_lambda[i][state.fep_state];
2345 pull_t *pull = nullptr;
2349 pull = set_pull_init(ir, &sys, state.x.rvec_array(), state.box, state.lambda[efptMASS], wi);
2352 /* Modules that supply external potential for pull coordinates
2353 * should register those potentials here. finish_pull() will check
2354 * that providers have been registerd for all external potentials.
2359 setStateDependentAwhParams(ir->awhParams, ir->pull, pull,
2360 state.box, ir->ePBC, &ir->opts, wi);
2370 set_reference_positions(ir->rot, state.x.rvec_array(), state.box,
2371 opt2fn("-ref", NFILE, fnm), opt2bSet("-ref", NFILE, fnm),
2375 /* reset_multinr(sys); */
2377 if (EEL_PME(ir->coulombtype))
2379 float ratio = pme_load_estimate(&sys, ir, state.box);
2380 fprintf(stderr, "Estimate for the relative computational load of the PME mesh part: %.2f\n", ratio);
2381 /* With free energy we might need to do PME both for the A and B state
2382 * charges. This will double the cost, but the optimal performance will
2383 * then probably be at a slightly larger cut-off and grid spacing.
2385 if ((ir->efep == efepNO && ratio > 1.0/2.0) ||
2386 (ir->efep != efepNO && ratio > 2.0/3.0))
2389 "The optimal PME mesh load for parallel simulations is below 0.5\n"
2390 "and for highly parallel simulations between 0.25 and 0.33,\n"
2391 "for higher performance, increase the cut-off and the PME grid spacing.\n");
2392 if (ir->efep != efepNO)
2395 "For free energy simulations, the optimal load limit increases from 0.5 to 0.667\n");
2401 char warn_buf[STRLEN];
2402 double cio = compute_io(ir, sys.natoms, sys.groups, F_NRE, 1);
2403 sprintf(warn_buf, "This run will generate roughly %.0f Mb of data", cio);
2406 set_warning_line(wi, mdparin, -1);
2407 warning_note(wi, warn_buf);
2411 printf("%s\n", warn_buf);
2417 fprintf(stderr, "writing run input file...\n");
2420 done_warning(wi, FARGS);
2421 write_tpx_state(ftp2fn(efTPR, NFILE, fnm), ir, &state, &sys);
2423 /* Output IMD group, if bIMD is TRUE */
2424 gmx::write_IMDgroup_to_file(ir->bIMD, ir, &state, &sys, NFILE, fnm);
2426 sfree(opts->define);
2427 sfree(opts->include);
2429 for (auto &mol : mi)
2431 // Some of the contents of molinfo have been stolen, so
2432 // fullCleanUp can't be called.
2433 mol.partialCleanUp();
2435 done_inputrec_strings();
2436 output_env_done(oenv);