*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013-2020,2021, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2016,2017, The GROMACS development team.
+ * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
#include <cstdlib>
#include <algorithm>
+#include <memory>
#include <string>
-#include "gromacs/awh/read_params.h"
+#include "gromacs/applied_forces/awh/read_params.h"
#include "gromacs/fileio/readinp.h"
#include "gromacs/fileio/warninp.h"
#include "gromacs/gmxlib/network.h"
#include "gromacs/mdrun/mdmodules.h"
#include "gromacs/mdtypes/inputrec.h"
#include "gromacs/mdtypes/md_enums.h"
+#include "gromacs/mdtypes/multipletimestepping.h"
#include "gromacs/mdtypes/pull_params.h"
#include "gromacs/options/options.h"
#include "gromacs/options/treesupport.h"
* message.
*/
-typedef struct t_inputrec_strings
+struct gmx_inputrec_strings
{
char tcgrps[STRLEN], tau_t[STRLEN], ref_t[STRLEN], acc[STRLEN], accgrps[STRLEN], freeze[STRLEN],
frdim[STRLEN], energy[STRLEN], user1[STRLEN], user2[STRLEN], vcm[STRLEN],
x_compressed_groups[STRLEN], couple_moltype[STRLEN], orirefitgrp[STRLEN],
egptable[STRLEN], egpexcl[STRLEN], wall_atomtype[STRLEN], wall_density[STRLEN],
deform[STRLEN], QMMM[STRLEN], imd_grp[STRLEN];
- char fep_lambda[efptNR][STRLEN];
- char lambda_weights[STRLEN];
- char** pull_grp;
- char** rot_grp;
- char anneal[STRLEN], anneal_npoints[STRLEN], anneal_time[STRLEN], anneal_temp[STRLEN];
- char QMmethod[STRLEN], QMbasis[STRLEN], QMcharge[STRLEN], QMmult[STRLEN], bSH[STRLEN],
- CASorbitals[STRLEN], CASelectrons[STRLEN], SAon[STRLEN], SAoff[STRLEN], SAsteps[STRLEN];
-
-} gmx_inputrec_strings;
+ char fep_lambda[efptNR][STRLEN];
+ char lambda_weights[STRLEN];
+ std::vector<std::string> pullGroupNames;
+ std::vector<std::string> rotateGroupNames;
+ char anneal[STRLEN], anneal_npoints[STRLEN], anneal_time[STRLEN], anneal_temp[STRLEN];
+};
-static gmx_inputrec_strings* is = nullptr;
+static gmx_inputrec_strings* inputrecStrings = nullptr;
void init_inputrec_strings()
{
- if (is)
+ if (inputrecStrings)
{
gmx_incons(
"Attempted to call init_inputrec_strings before calling done_inputrec_strings. "
"Only one inputrec (i.e. .mdp file) can be parsed at a time.");
}
- snew(is, 1);
+ inputrecStrings = new gmx_inputrec_strings();
}
void done_inputrec_strings()
{
- sfree(is);
- is = nullptr;
+ delete inputrecStrings;
+ inputrecStrings = nullptr;
}
}
}
+static void checkMtsRequirement(const t_inputrec& ir, const char* param, const int nstValue, warninp_t wi)
+{
+ GMX_RELEASE_ASSERT(ir.mtsLevels.size() >= 2, "Need at least two levels for MTS");
+ const int mtsFactor = ir.mtsLevels.back().stepFactor;
+ if (nstValue % mtsFactor != 0)
+ {
+ auto message = gmx::formatString(
+ "With MTS, %s = %d should be a multiple of mts-factor = %d", param, nstValue, mtsFactor);
+ warning_error(wi, message.c_str());
+ }
+}
+
+static void setupMtsLevels(gmx::ArrayRef<gmx::MtsLevel> mtsLevels,
+ const t_inputrec& ir,
+ const t_gromppopts& opts,
+ warninp_t wi)
+{
+ /* MD-VV has no MTS support yet.
+ * SD1 needs different scaling coefficients for the different MTS forces
+ * and the different forces are currently not available in ForceBuffers.
+ */
+ if (ir.eI != eiMD)
+ {
+ auto message = gmx::formatString(
+ "Multiple time stepping is only supported with integrator %s", ei_names[eiMD]);
+ warning_error(wi, message.c_str());
+ }
+ if (opts.numMtsLevels != 2)
+ {
+ warning_error(wi, "Only mts-levels = 2 is supported");
+ }
+ else
+ {
+ const std::vector<std::string> inputForceGroups = gmx::splitString(opts.mtsLevel2Forces);
+ auto& forceGroups = mtsLevels[1].forceGroups;
+ for (const auto& inputForceGroup : inputForceGroups)
+ {
+ bool found = false;
+ int nameIndex = 0;
+ for (const auto& forceGroupName : gmx::mtsForceGroupNames)
+ {
+ if (gmx::equalCaseInsensitive(inputForceGroup, forceGroupName))
+ {
+ forceGroups.set(nameIndex);
+ found = true;
+ }
+ nameIndex++;
+ }
+ if (!found)
+ {
+ auto message =
+ gmx::formatString("Unknown MTS force group '%s'", inputForceGroup.c_str());
+ warning_error(wi, message.c_str());
+ }
+ }
+
+ if (mtsLevels[1].stepFactor <= 1)
+ {
+ gmx_fatal(FARGS, "mts-factor should be larger than 1");
+ }
+
+ // Make the level 0 use the complement of the force groups of group 1
+ mtsLevels[0].forceGroups = ~mtsLevels[1].forceGroups;
+ mtsLevels[0].stepFactor = 1;
+
+ if ((EEL_FULL(ir.coulombtype) || EVDW_PME(ir.vdwtype))
+ && !mtsLevels[1].forceGroups[static_cast<int>(gmx::MtsForceGroups::LongrangeNonbonded)])
+ {
+ warning_error(wi,
+ "With long-range electrostatics and/or LJ treatment, the long-range part "
+ "has to be part of the mts-level2-forces");
+ }
+
+ if (ir.nstcalcenergy > 0)
+ {
+ checkMtsRequirement(ir, "nstcalcenergy", ir.nstcalcenergy, wi);
+ }
+ checkMtsRequirement(ir, "nstenergy", ir.nstenergy, wi);
+ checkMtsRequirement(ir, "nstlog", ir.nstlog, wi);
+ if (ir.efep != efepNO)
+ {
+ checkMtsRequirement(ir, "nstdhdl", ir.fepvals->nstdhdl, wi);
+ }
+
+ if (ir.bPull)
+ {
+ const int pullMtsLevel = gmx::forceGroupMtsLevel(ir.mtsLevels, gmx::MtsForceGroups::Pull);
+ if (ir.pull->nstxout % ir.mtsLevels[pullMtsLevel].stepFactor != 0)
+ {
+ warning_error(wi, "pull-nstxout should be a multiple of mts-factor");
+ }
+ if (ir.pull->nstfout % ir.mtsLevels[pullMtsLevel].stepFactor != 0)
+ {
+ warning_error(wi, "pull-nstfout should be a multiple of mts-factor");
+ }
+ }
+ }
+}
+
void check_ir(const char* mdparin,
const gmx::MdModulesNotifier& mdModulesNotifier,
t_inputrec* ir,
{
ir->nstpcouple = ir_optimal_nstpcouple(ir);
}
+ if (ir->useMts && ir->nstpcouple % ir->mtsLevels.back().stepFactor != 0)
+ {
+ warning_error(wi,
+ "With multiple time stepping, nstpcouple should be a mutiple of "
+ "mts-factor");
+ }
}
if (ir->nstcalcenergy > 0)
// Inquire all MdModules, if their parameters match with the energy
// calculation frequency
gmx::EnergyCalculationFrequencyErrors energyCalculationFrequencyErrors(ir->nstcalcenergy);
- mdModulesNotifier.notifier_.notify(&energyCalculationFrequencyErrors);
+ mdModulesNotifier.preProcessingNotifications_.notify(&energyCalculationFrequencyErrors);
// Emit all errors from the energy calculation frequency checks
for (const std::string& energyFrequencyErrorMessage :
/* TPI STUFF */
if (EI_TPI(ir->eI))
{
- sprintf(err_buf, "TPI only works with pbc = %s", epbc_names[epbcXYZ]);
- CHECK(ir->ePBC != epbcXYZ);
+ sprintf(err_buf, "TPI only works with pbc = %s", c_pbcTypeNames[PbcType::Xyz].c_str());
+ CHECK(ir->pbcType != PbcType::Xyz);
sprintf(err_buf, "with TPI nstlist should be larger than zero");
CHECK(ir->nstlist <= 0);
sprintf(err_buf, "TPI does not work with full electrostatics other than PME");
sprintf(err_buf, "The soft-core power is %d and can only be 1 or 2", fep->sc_power);
CHECK(fep->sc_alpha != 0 && fep->sc_power != 1 && fep->sc_power != 2);
- sprintf(err_buf, "The soft-core sc-r-power is %d and can only be 6 or 48",
+ sprintf(err_buf,
+ "The soft-core sc-r-power is %d and can only be 6. (sc-r-power 48 is no longer "
+ "supported.)",
static_cast<int>(fep->sc_r_power));
- CHECK(fep->sc_alpha != 0 && fep->sc_r_power != 6.0 && fep->sc_r_power != 48.0);
+ CHECK(fep->sc_alpha != 0 && fep->sc_r_power != 6.0);
sprintf(err_buf,
"Can't use positive delta-lambda (%g) if initial state/lambda does not start at "
}
/* PBC/WALLS */
- sprintf(err_buf, "walls only work with pbc=%s", epbc_names[epbcXY]);
- CHECK(ir->nwall && ir->ePBC != epbcXY);
+ sprintf(err_buf, "walls only work with pbc=%s", c_pbcTypeNames[PbcType::XY].c_str());
+ CHECK(ir->nwall && ir->pbcType != PbcType::XY);
/* VACUUM STUFF */
- if (ir->ePBC != epbcXYZ && ir->nwall != 2)
+ if (ir->pbcType != PbcType::Xyz && ir->nwall != 2)
{
- if (ir->ePBC == epbcNONE)
+ if (ir->pbcType == PbcType::No)
{
if (ir->epc != epcNO)
{
}
else
{
- sprintf(err_buf, "Can not have pressure coupling with pbc=%s", epbc_names[ir->ePBC]);
+ sprintf(err_buf, "Can not have pressure coupling with pbc=%s",
+ c_pbcTypeNames[ir->pbcType].c_str());
CHECK(ir->epc != epcNO);
}
- sprintf(err_buf, "Can not have Ewald with pbc=%s", epbc_names[ir->ePBC]);
+ sprintf(err_buf, "Can not have Ewald with pbc=%s", c_pbcTypeNames[ir->pbcType].c_str());
CHECK(EEL_FULL(ir->coulombtype));
- sprintf(err_buf, "Can not have dispersion correction with pbc=%s", epbc_names[ir->ePBC]);
+ sprintf(err_buf, "Can not have dispersion correction with pbc=%s",
+ c_pbcTypeNames[ir->pbcType].c_str());
CHECK(ir->eDispCorr != edispcNO);
}
"with coulombtype = %s or coulombtype = %s\n"
"without periodic boundary conditions (pbc = %s) and\n"
"rcoulomb and rvdw set to zero",
- eel_names[eelCUT], eel_names[eelUSER], epbc_names[epbcNONE]);
+ eel_names[eelCUT], eel_names[eelUSER], c_pbcTypeNames[PbcType::No].c_str());
CHECK(((ir->coulombtype != eelCUT) && (ir->coulombtype != eelUSER))
- || (ir->ePBC != epbcNONE) || (ir->rcoulomb != 0.0) || (ir->rvdw != 0.0));
+ || (ir->pbcType != PbcType::No) || (ir->rcoulomb != 0.0) || (ir->rvdw != 0.0));
if (ir->nstlist > 0)
{
"Can not remove the rotation around the center of mass with periodic "
"molecules");
CHECK(ir->bPeriodicMols);
- if (ir->ePBC != epbcNONE)
+ if (ir->pbcType != PbcType::No)
{
warning(wi,
"Removing the rotation around the center of mass in a periodic system, "
}
}
- if (EI_STATE_VELOCITY(ir->eI) && !EI_SD(ir->eI) && ir->ePBC == epbcNONE && ir->comm_mode != ecmANGULAR)
+ if (EI_STATE_VELOCITY(ir->eI) && !EI_SD(ir->eI) && ir->pbcType == PbcType::No && ir->comm_mode != ecmANGULAR)
{
sprintf(warn_buf,
"Tumbling and flying ice-cubes: We are not removing rotation around center of mass "
{
if (ir->ewald_geometry == eewg3D)
{
- sprintf(warn_buf, "With pbc=%s you should use ewald-geometry=%s", epbc_names[ir->ePBC],
- eewg_names[eewg3DC]);
+ sprintf(warn_buf, "With pbc=%s you should use ewald-geometry=%s",
+ c_pbcTypeNames[ir->pbcType].c_str(), eewg_names[eewg3DC]);
warning(wi, warn_buf);
}
/* This check avoids extra pbc coding for exclusion corrections */
sprintf(err_buf, "wall-ewald-zfac should be >= 2");
CHECK(ir->wall_ewald_zfac < 2);
}
- if ((ir->ewald_geometry == eewg3DC) && (ir->ePBC != epbcXY) && EEL_FULL(ir->coulombtype))
+ if ((ir->ewald_geometry == eewg3DC) && (ir->pbcType != PbcType::XY) && EEL_FULL(ir->coulombtype))
{
sprintf(warn_buf, "With %s and ewald_geometry = %s you should use pbc = %s",
- eel_names[ir->coulombtype], eewg_names[eewg3DC], epbc_names[epbcXY]);
+ eel_names[ir->coulombtype], eewg_names[eewg3DC], c_pbcTypeNames[PbcType::XY].c_str());
warning(wi, warn_buf);
}
if ((ir->epsilon_surface != 0) && EEL_FULL(ir->coulombtype))
if (ir->bQMMM)
{
- warning_error(wi, "QMMM is currently not supported");
- if (!EI_DYNAMICS(ir->eI))
- {
- char buf[STRLEN];
- sprintf(buf, "QMMM is only supported with dynamics, not with integrator %s", ei_names[ir->eI]);
- warning_error(wi, buf);
- }
+ warning_error(wi, "The QMMM integration you are trying to use is no longer supported");
}
if (ir->bAdress)
{
gmx_fatal(FARGS, "AdResS simulations are no longer supported");
}
+
+ // cosine acceleration is only supported in leap-frog
+ if (ir->cos_accel != 0.0 && ir->eI != eiMD)
+ {
+ warning_error(wi, "cos-acceleration is only supported by integrator = md");
+ }
}
/* interpret a number of doubles from a string and put them in an array,
}
- /* make it easier if sc_r_power = 48 by increasing it to the 4th power, to be in the right scale. */
- if (fep->sc_r_power == 48)
- {
- if (fep->sc_alpha > 0.1)
- {
- gmx_fatal(FARGS,
- "sc_alpha (%f) for sc_r_power = 48 should usually be between 0.001 and 0.004",
- fep->sc_alpha);
- }
- }
-
/* now read in the weights */
parse_n_real(weights, &nweights, &(expand->init_lambda_weights), wi);
if (nweights == 0)
GetSimTemps(ir->fepvals->n_lambda, ir->simtempvals, ir->fepvals->all_lambda[efptTEMPERATURE]);
}
-static void convertYesNos(warninp_t /*wi*/,
- gmx::ArrayRef<const std::string> inputs,
- const char* /*name*/,
- gmx_bool* outputs)
-{
- int i = 0;
- for (const auto& input : inputs)
- {
- outputs[i] = gmx::equalCaseInsensitive(input, "Y", 1);
- ++i;
- }
-}
-
template<typename T>
void convertInts(warninp_t wi, gmx::ArrayRef<const std::string> inputs, const char* name, T* outputs)
{
gmx_fatal(FARGS, "Expected %d elements for wall_atomtype, found %zu", ir->nwall,
wallAtomTypes.size());
}
+ GMX_RELEASE_ASSERT(ir->nwall < 3, "Invalid number of walls");
for (int i = 0; i < ir->nwall; i++)
{
opts->wall_atomtype[i] = gmx_strdup(wallAtomTypes[i].c_str());
printStringNoNewline(
&inp, "Part index is updated automatically on checkpointing (keeps files separate)");
ir->simulation_part = get_eint(&inp, "simulation-part", 1, wi);
+ printStringNoNewline(&inp, "Multiple time-stepping");
+ ir->useMts = (get_eeenum(&inp, "mts", yesno_names, wi) != 0);
+ if (ir->useMts)
+ {
+ opts->numMtsLevels = get_eint(&inp, "mts-levels", 2, wi);
+ ir->mtsLevels.resize(2);
+ gmx::MtsLevel& mtsLevel = ir->mtsLevels[1];
+ opts->mtsLevel2Forces = setStringEntry(&inp, "mts-level2-forces", "longrange-nonbonded");
+ mtsLevel.stepFactor = get_eint(&inp, "mts-level2-factor", 2, wi);
+
+ // We clear after reading without dynamics to not force the user to remove MTS mdp options
+ if (!EI_DYNAMICS(ir->eI))
+ {
+ ir->useMts = false;
+ ir->mtsLevels.clear();
+ }
+ }
printStringNoNewline(&inp, "mode for center of mass motion removal");
ir->comm_mode = get_eeenum(&inp, "comm-mode", ecm_names, wi);
printStringNoNewline(&inp, "number of steps for center of mass motion removal");
ir->nstcomm = get_eint(&inp, "nstcomm", 100, wi);
printStringNoNewline(&inp, "group(s) for center of mass motion removal");
- setStringEntry(&inp, "comm-grps", is->vcm, nullptr);
+ setStringEntry(&inp, "comm-grps", inputrecStrings->vcm, nullptr);
printStringNewline(&inp, "LANGEVIN DYNAMICS OPTIONS");
printStringNoNewline(&inp, "Friction coefficient (amu/ps) and random seed");
printStringNoNewline(&inp, "This selects the subset of atoms for the compressed");
printStringNoNewline(&inp, "trajectory file. You can select multiple groups. By");
printStringNoNewline(&inp, "default, all atoms will be written.");
- setStringEntry(&inp, "compressed-x-grps", is->x_compressed_groups, nullptr);
+ setStringEntry(&inp, "compressed-x-grps", inputrecStrings->x_compressed_groups, nullptr);
printStringNoNewline(&inp, "Selection of energy groups");
- setStringEntry(&inp, "energygrps", is->energy, nullptr);
+ setStringEntry(&inp, "energygrps", inputrecStrings->energy, nullptr);
/* Neighbor searching */
printStringNewline(&inp, "NEIGHBORSEARCHING PARAMETERS");
printStringNoNewline(&inp, "nblist update frequency");
ir->nstlist = get_eint(&inp, "nstlist", 10, wi);
printStringNoNewline(&inp, "Periodic boundary conditions: xyz, no, xy");
- ir->ePBC = get_eeenum(&inp, "pbc", epbc_names, wi);
+ // TODO This conversion should be removed when proper std:string handling will be added to get_eeenum(...), etc.
+ std::vector<const char*> pbcTypesNamesChar;
+ for (const auto& pbcTypeName : c_pbcTypeNames)
+ {
+ pbcTypesNamesChar.push_back(pbcTypeName.c_str());
+ }
+ ir->pbcType = static_cast<PbcType>(get_eeenum(&inp, "pbc", pbcTypesNamesChar.data(), wi));
ir->bPeriodicMols = get_eeenum(&inp, "periodic-molecules", yesno_names, wi) != 0;
printStringNoNewline(&inp,
"Allowed energy error due to the Verlet buffer in kJ/mol/ps per atom,");
printStringNoNewline(&inp, "Extension of the potential lookup tables beyond the cut-off");
ir->tabext = get_ereal(&inp, "table-extension", 1.0, wi);
printStringNoNewline(&inp, "Separate tables between energy group pairs");
- setStringEntry(&inp, "energygrp-table", is->egptable, nullptr);
+ setStringEntry(&inp, "energygrp-table", inputrecStrings->egptable, nullptr);
printStringNoNewline(&inp, "Spacing for the PME/PPPM FFT grid");
ir->fourier_spacing = get_ereal(&inp, "fourierspacing", 0.12, wi);
printStringNoNewline(&inp, "FFT grid size, when a value is 0 fourierspacing will be used");
ir->opts.nhchainlength = get_eint(&inp, "nh-chain-length", 10, wi);
ir->bPrintNHChains = (get_eeenum(&inp, "print-nose-hoover-chain-variables", yesno_names, wi) != 0);
printStringNoNewline(&inp, "Groups to couple separately");
- setStringEntry(&inp, "tc-grps", is->tcgrps, nullptr);
+ setStringEntry(&inp, "tc-grps", inputrecStrings->tcgrps, nullptr);
printStringNoNewline(&inp, "Time constant (ps) and reference temperature (K)");
- setStringEntry(&inp, "tau-t", is->tau_t, nullptr);
- setStringEntry(&inp, "ref-t", is->ref_t, nullptr);
+ setStringEntry(&inp, "tau-t", inputrecStrings->tau_t, nullptr);
+ setStringEntry(&inp, "ref-t", inputrecStrings->ref_t, nullptr);
printStringNoNewline(&inp, "pressure coupling");
ir->epc = get_eeenum(&inp, "pcoupl", epcoupl_names, wi);
ir->epct = get_eeenum(&inp, "pcoupltype", epcoupltype_names, wi);
/* QMMM */
printStringNewline(&inp, "OPTIONS FOR QMMM calculations");
ir->bQMMM = (get_eeenum(&inp, "QMMM", yesno_names, wi) != 0);
- printStringNoNewline(&inp, "Groups treated Quantum Mechanically");
- setStringEntry(&inp, "QMMM-grps", is->QMMM, nullptr);
- printStringNoNewline(&inp, "QM method");
- setStringEntry(&inp, "QMmethod", is->QMmethod, nullptr);
- printStringNoNewline(&inp, "QMMM scheme");
- ir->QMMMscheme = get_eeenum(&inp, "QMMMscheme", eQMMMscheme_names, wi);
- printStringNoNewline(&inp, "QM basisset");
- setStringEntry(&inp, "QMbasis", is->QMbasis, nullptr);
- printStringNoNewline(&inp, "QM charge");
- setStringEntry(&inp, "QMcharge", is->QMcharge, nullptr);
- printStringNoNewline(&inp, "QM multiplicity");
- setStringEntry(&inp, "QMmult", is->QMmult, nullptr);
- printStringNoNewline(&inp, "Surface Hopping");
- setStringEntry(&inp, "SH", is->bSH, nullptr);
- printStringNoNewline(&inp, "CAS space options");
- setStringEntry(&inp, "CASorbitals", is->CASorbitals, nullptr);
- setStringEntry(&inp, "CASelectrons", is->CASelectrons, nullptr);
- setStringEntry(&inp, "SAon", is->SAon, nullptr);
- setStringEntry(&inp, "SAoff", is->SAoff, nullptr);
- setStringEntry(&inp, "SAsteps", is->SAsteps, nullptr);
- printStringNoNewline(&inp, "Scale factor for MM charges");
- ir->scalefactor = get_ereal(&inp, "MMChargeScaleFactor", 1.0, wi);
+ printStringNoNewline(&inp, "Groups treated with MiMiC");
+ setStringEntry(&inp, "QMMM-grps", inputrecStrings->QMMM, nullptr);
/* Simulated annealing */
printStringNewline(&inp, "SIMULATED ANNEALING");
printStringNoNewline(&inp, "Type of annealing for each temperature group (no/single/periodic)");
- setStringEntry(&inp, "annealing", is->anneal, nullptr);
+ setStringEntry(&inp, "annealing", inputrecStrings->anneal, nullptr);
printStringNoNewline(&inp,
"Number of time points to use for specifying annealing in each group");
- setStringEntry(&inp, "annealing-npoints", is->anneal_npoints, nullptr);
+ setStringEntry(&inp, "annealing-npoints", inputrecStrings->anneal_npoints, nullptr);
printStringNoNewline(&inp, "List of times at the annealing points for each group");
- setStringEntry(&inp, "annealing-time", is->anneal_time, nullptr);
+ setStringEntry(&inp, "annealing-time", inputrecStrings->anneal_time, nullptr);
printStringNoNewline(&inp, "Temp. at each annealing point, for each group.");
- setStringEntry(&inp, "annealing-temp", is->anneal_temp, nullptr);
+ setStringEntry(&inp, "annealing-temp", inputrecStrings->anneal_temp, nullptr);
/* Startup run */
printStringNewline(&inp, "GENERATE VELOCITIES FOR STARTUP RUN");
printStringNewline(&inp, "ENERGY GROUP EXCLUSIONS");
printStringNoNewline(
&inp, "Pairs of energy groups for which all non-bonded interactions are excluded");
- setStringEntry(&inp, "energygrp-excl", is->egpexcl, nullptr);
+ setStringEntry(&inp, "energygrp-excl", inputrecStrings->egpexcl, nullptr);
/* Walls */
printStringNewline(&inp, "WALLS");
ir->nwall = get_eint(&inp, "nwall", 0, wi);
ir->wall_type = get_eeenum(&inp, "wall-type", ewt_names, wi);
ir->wall_r_linpot = get_ereal(&inp, "wall-r-linpot", -1, wi);
- setStringEntry(&inp, "wall-atomtype", is->wall_atomtype, nullptr);
- setStringEntry(&inp, "wall-density", is->wall_density, nullptr);
+ setStringEntry(&inp, "wall-atomtype", inputrecStrings->wall_atomtype, nullptr);
+ setStringEntry(&inp, "wall-density", inputrecStrings->wall_density, nullptr);
ir->wall_ewald_zfac = get_ereal(&inp, "wall-ewald-zfac", 3, wi);
/* COM pulling */
ir->bPull = (get_eeenum(&inp, "pull", yesno_names, wi) != 0);
if (ir->bPull)
{
- snew(ir->pull, 1);
- is->pull_grp = read_pullparams(&inp, ir->pull, wi);
+ ir->pull = std::make_unique<pull_params_t>();
+ inputrecStrings->pullGroupNames = read_pullparams(&inp, ir->pull.get(), wi);
+
+ if (ir->useMts)
+ {
+ for (int c = 0; c < ir->pull->ncoord; c++)
+ {
+ if (ir->pull->coord[c].eType == epullCONSTRAINT)
+ {
+ warning_error(wi,
+ "Constraint COM pulling is not supported in combination with "
+ "multiple time stepping");
+ break;
+ }
+ }
+ }
}
/* AWH biasing
- NOTE: needs COM pulling input */
+ NOTE: needs COM pulling or free energy input */
printStringNewline(&inp, "AWH biasing");
ir->bDoAwh = (get_eeenum(&inp, "awh", yesno_names, wi) != 0);
if (ir->bDoAwh)
{
- if (ir->bPull)
- {
- ir->awhParams = gmx::readAndCheckAwhParams(&inp, ir, wi);
- }
- else
- {
- gmx_fatal(FARGS, "AWH biasing is only compatible with COM pulling turned on");
- }
+ ir->awhParams = gmx::readAwhParams(&inp, wi);
}
/* Enforced rotation */
if (ir->bRot)
{
snew(ir->rot, 1);
- is->rot_grp = read_rotparams(&inp, ir->rot, wi);
+ inputrecStrings->rotateGroupNames = read_rotparams(&inp, ir->rot, wi);
}
/* Interactive MD */
ir->bIMD = FALSE;
printStringNewline(&inp, "Group to display and/or manipulate in interactive MD session");
- setStringEntry(&inp, "IMD-group", is->imd_grp, nullptr);
- if (is->imd_grp[0] != '\0')
+ setStringEntry(&inp, "IMD-group", inputrecStrings->imd_grp, nullptr);
+ if (inputrecStrings->imd_grp[0] != '\0')
{
snew(ir->imd, 1);
ir->bIMD = TRUE;
printStringNoNewline(&inp, "Orientation restraints force constant and tau for time averaging");
ir->orires_fc = get_ereal(&inp, "orire-fc", 0.0, wi);
ir->orires_tau = get_ereal(&inp, "orire-tau", 0.0, wi);
- setStringEntry(&inp, "orire-fitgrp", is->orirefitgrp, nullptr);
+ setStringEntry(&inp, "orire-fitgrp", inputrecStrings->orirefitgrp, nullptr);
printStringNoNewline(&inp, "Output frequency for trace(SD) and S to energy file");
ir->nstorireout = get_eint(&inp, "nstorireout", 100, wi);
/* free energy variables */
printStringNewline(&inp, "Free energy variables");
ir->efep = get_eeenum(&inp, "free-energy", efep_names, wi);
- setStringEntry(&inp, "couple-moltype", is->couple_moltype, nullptr);
+ setStringEntry(&inp, "couple-moltype", inputrecStrings->couple_moltype, nullptr);
opts->couple_lam0 = get_eeenum(&inp, "couple-lambda0", couple_lam, wi);
opts->couple_lam1 = get_eeenum(&inp, "couple-lambda1", couple_lam, wi);
opts->bCoupleIntra = (get_eeenum(&inp, "couple-intramol", yesno_names, wi) != 0);
fep->init_fep_state = get_eint(&inp, "init-lambda-state", -1, wi);
fep->delta_lambda = get_ereal(&inp, "delta-lambda", 0.0, wi);
fep->nstdhdl = get_eint(&inp, "nstdhdl", 50, wi);
- setStringEntry(&inp, "fep-lambdas", is->fep_lambda[efptFEP], nullptr);
- setStringEntry(&inp, "mass-lambdas", is->fep_lambda[efptMASS], nullptr);
- setStringEntry(&inp, "coul-lambdas", is->fep_lambda[efptCOUL], nullptr);
- setStringEntry(&inp, "vdw-lambdas", is->fep_lambda[efptVDW], nullptr);
- setStringEntry(&inp, "bonded-lambdas", is->fep_lambda[efptBONDED], nullptr);
- setStringEntry(&inp, "restraint-lambdas", is->fep_lambda[efptRESTRAINT], nullptr);
- setStringEntry(&inp, "temperature-lambdas", is->fep_lambda[efptTEMPERATURE], nullptr);
+ setStringEntry(&inp, "fep-lambdas", inputrecStrings->fep_lambda[efptFEP], nullptr);
+ setStringEntry(&inp, "mass-lambdas", inputrecStrings->fep_lambda[efptMASS], nullptr);
+ setStringEntry(&inp, "coul-lambdas", inputrecStrings->fep_lambda[efptCOUL], nullptr);
+ setStringEntry(&inp, "vdw-lambdas", inputrecStrings->fep_lambda[efptVDW], nullptr);
+ setStringEntry(&inp, "bonded-lambdas", inputrecStrings->fep_lambda[efptBONDED], nullptr);
+ setStringEntry(&inp, "restraint-lambdas", inputrecStrings->fep_lambda[efptRESTRAINT], nullptr);
+ setStringEntry(&inp, "temperature-lambdas", inputrecStrings->fep_lambda[efptTEMPERATURE], nullptr);
fep->lambda_neighbors = get_eint(&inp, "calc-lambda-neighbors", 1, wi);
- setStringEntry(&inp, "init-lambda-weights", is->lambda_weights, nullptr);
+ setStringEntry(&inp, "init-lambda-weights", inputrecStrings->lambda_weights, nullptr);
fep->edHdLPrintEnergy = get_eeenum(&inp, "dhdl-print-energy", edHdLPrintEnergy_names, wi);
fep->sc_alpha = get_ereal(&inp, "sc-alpha", 0.0, wi);
fep->sc_power = get_eint(&inp, "sc-power", 1, wi);
/* Non-equilibrium MD stuff */
printStringNewline(&inp, "Non-equilibrium MD stuff");
- setStringEntry(&inp, "acc-grps", is->accgrps, nullptr);
- setStringEntry(&inp, "accelerate", is->acc, nullptr);
- setStringEntry(&inp, "freezegrps", is->freeze, nullptr);
- setStringEntry(&inp, "freezedim", is->frdim, nullptr);
+ setStringEntry(&inp, "acc-grps", inputrecStrings->accgrps, nullptr);
+ setStringEntry(&inp, "accelerate", inputrecStrings->acc, nullptr);
+ setStringEntry(&inp, "freezegrps", inputrecStrings->freeze, nullptr);
+ setStringEntry(&inp, "freezedim", inputrecStrings->frdim, nullptr);
ir->cos_accel = get_ereal(&inp, "cos-acceleration", 0, wi);
- setStringEntry(&inp, "deform", is->deform, nullptr);
+ setStringEntry(&inp, "deform", inputrecStrings->deform, nullptr);
/* simulated tempering variables */
printStringNewline(&inp, "simulated tempering variables");
/* User defined thingies */
printStringNewline(&inp, "User defined thingies");
- setStringEntry(&inp, "user1-grps", is->user1, nullptr);
- setStringEntry(&inp, "user2-grps", is->user2, nullptr);
+ setStringEntry(&inp, "user1-grps", inputrecStrings->user1, nullptr);
+ setStringEntry(&inp, "user2-grps", inputrecStrings->user2, nullptr);
ir->userint1 = get_eint(&inp, "userint1", 0, wi);
ir->userint2 = get_eint(&inp, "userint2", 0, wi);
ir->userint3 = get_eint(&inp, "userint3", 0, wi);
}
opts->couple_moltype = nullptr;
- if (strlen(is->couple_moltype) > 0)
+ if (strlen(inputrecStrings->couple_moltype) > 0)
{
if (ir->efep != efepNO)
{
- opts->couple_moltype = gmx_strdup(is->couple_moltype);
+ opts->couple_moltype = gmx_strdup(inputrecStrings->couple_moltype);
if (opts->couple_lam0 == opts->couple_lam1)
{
warning(wi, "The lambda=0 and lambda=1 states for coupling are identical");
/* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
if (ir->efep != efepNO)
{
- if (fep->delta_lambda > 0)
+ if (fep->delta_lambda != 0)
{
ir->efep = efepSLOWGROWTH;
}
{
ir->bExpanded = TRUE;
}
- do_fep_params(ir, is->fep_lambda, is->lambda_weights, wi);
+ do_fep_params(ir, inputrecStrings->fep_lambda, inputrecStrings->lambda_weights, wi);
if (ir->bSimTemp) /* done after fep params */
{
do_simtemp_params(ir);
/* WALL PARAMETERS */
- do_wall_params(ir, is->wall_atomtype, is->wall_density, opts, wi);
+ do_wall_params(ir, inputrecStrings->wall_atomtype, inputrecStrings->wall_density, opts, wi);
/* ORIENTATION RESTRAINT PARAMETERS */
- if (opts->bOrire && gmx::splitString(is->orirefitgrp).size() != 1)
+ if (opts->bOrire && gmx::splitString(inputrecStrings->orirefitgrp).size() != 1)
{
warning_error(wi, "ERROR: Need one orientation restraint fit group\n");
}
}
double gmx_unused canary;
- int ndeform = sscanf(is->deform, "%lf %lf %lf %lf %lf %lf %lf", &(dumdub[0][0]), &(dumdub[0][1]),
- &(dumdub[0][2]), &(dumdub[0][3]), &(dumdub[0][4]), &(dumdub[0][5]), &canary);
+ int ndeform = sscanf(inputrecStrings->deform, "%lf %lf %lf %lf %lf %lf %lf", &(dumdub[0][0]),
+ &(dumdub[0][1]), &(dumdub[0][2]), &(dumdub[0][3]), &(dumdub[0][4]),
+ &(dumdub[0][5]), &canary);
- if (strlen(is->deform) > 0 && ndeform != 6)
+ if (strlen(inputrecStrings->deform) > 0 && ndeform != 6)
{
- warning_error(
- wi, gmx::formatString(
- "Cannot parse exactly 6 box deformation velocities from string '%s'", is->deform)
- .c_str());
+ warning_error(wi,
+ gmx::formatString(
+ "Cannot parse exactly 6 box deformation velocities from string '%s'",
+ inputrecStrings->deform)
+ .c_str());
}
for (i = 0; i < 3; i++)
{
}
}
- sfree(dumstr[0]);
- sfree(dumstr[1]);
-}
-
-static int search_QMstring(const char* s, int ng, const char* gn[])
-{
- /* same as normal search_string, but this one searches QM strings */
- int i;
+ /* Set up MTS levels, this needs to happen before checking AWH parameters */
+ if (ir->useMts)
+ {
+ setupMtsLevels(ir->mtsLevels, *ir, *opts, wi);
+ }
- for (i = 0; (i < ng); i++)
+ if (ir->bDoAwh)
{
- if (gmx_strcasecmp(s, gn[i]) == 0)
- {
- return i;
- }
+ gmx::checkAwhParams(ir->awhParams, ir, wi);
}
- gmx_fatal(FARGS, "this QM method or basisset (%s) is not implemented\n!", s);
-} /* search_QMstring */
+ sfree(dumstr[0]);
+ sfree(dumstr[1]);
+}
/* We would like gn to be const as well, but C doesn't allow this */
/* TODO this is utility functionality (search for the index of a
* belong to different TC or VCM groups it is anyhow difficult
* to determine the optimal nrdf assignment.
*/
- pull = ir->pull;
+ pull = ir->pull.get();
for (int i = 0; i < pull->ncoord; i++)
{
pgrp = &pull->group[pull->coord[i].group[j]];
- if (pgrp->nat > 0)
+ if (!pgrp->ind.empty())
{
/* Subtract 1/2 dof from each group */
int ai = pgrp->ind[0];
set_warning_line(wi, mdparin, -1);
- auto temperatureCouplingTauValues = gmx::splitString(is->tau_t);
- auto temperatureCouplingReferenceValues = gmx::splitString(is->ref_t);
- auto temperatureCouplingGroupNames = gmx::splitString(is->tcgrps);
+ auto temperatureCouplingTauValues = gmx::splitString(inputrecStrings->tau_t);
+ auto temperatureCouplingReferenceValues = gmx::splitString(inputrecStrings->ref_t);
+ auto temperatureCouplingGroupNames = gmx::splitString(inputrecStrings->tcgrps);
if (temperatureCouplingTauValues.size() != temperatureCouplingGroupNames.size()
|| temperatureCouplingReferenceValues.size() != temperatureCouplingGroupNames.size())
{
}
/* Simulated annealing for each group. There are nr groups */
- auto simulatedAnnealingGroupNames = gmx::splitString(is->anneal);
+ auto simulatedAnnealingGroupNames = gmx::splitString(inputrecStrings->anneal);
if (simulatedAnnealingGroupNames.size() == 1
&& gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[0], "N", 1))
{
if (bAnneal)
{
/* Read the other fields too */
- auto simulatedAnnealingPoints = gmx::splitString(is->anneal_npoints);
+ auto simulatedAnnealingPoints = gmx::splitString(inputrecStrings->anneal_npoints);
if (simulatedAnnealingPoints.size() != simulatedAnnealingGroupNames.size())
{
gmx_fatal(FARGS, "Found %zu annealing-npoints values for %zu groups\n",
numSimulatedAnnealingFields += ir->opts.anneal_npoints[i];
}
- auto simulatedAnnealingTimes = gmx::splitString(is->anneal_time);
+ auto simulatedAnnealingTimes = gmx::splitString(inputrecStrings->anneal_time);
if (simulatedAnnealingTimes.size() != numSimulatedAnnealingFields)
{
gmx_fatal(FARGS, "Found %zu annealing-time values, wanted %zu\n",
simulatedAnnealingTimes.size(), numSimulatedAnnealingFields);
}
- auto simulatedAnnealingTemperatures = gmx::splitString(is->anneal_temp);
+ auto simulatedAnnealingTemperatures = gmx::splitString(inputrecStrings->anneal_temp);
if (simulatedAnnealingTemperatures.size() != numSimulatedAnnealingFields)
{
gmx_fatal(FARGS, "Found %zu annealing-temp values, wanted %zu\n",
{
for (int i = 1; i < ir->pull->ngroup; i++)
{
- const int gid = search_string(is->pull_grp[i], defaultIndexGroups->nr, gnames);
+ const int gid = search_string(inputrecStrings->pullGroupNames[i].c_str(),
+ defaultIndexGroups->nr, gnames);
GMX_ASSERT(defaultIndexGroups, "Must have initialized default index groups");
atomGroupRangeValidation(natoms, gid, *defaultIndexGroups);
}
- make_pull_groups(ir->pull, is->pull_grp, defaultIndexGroups, gnames);
+ process_pull_groups(ir->pull->group, inputrecStrings->pullGroupNames, defaultIndexGroups, gnames);
- make_pull_coords(ir->pull);
+ checkPullCoords(ir->pull->group, ir->pull->coord);
}
if (ir->bRot)
{
- make_rotation_groups(ir->rot, is->rot_grp, defaultIndexGroups, gnames);
+ make_rotation_groups(ir->rot, inputrecStrings->rotateGroupNames, defaultIndexGroups, gnames);
}
if (ir->eSwapCoords != eswapNO)
/* Make indices for IMD session */
if (ir->bIMD)
{
- make_IMD_group(ir->imd, is->imd_grp, defaultIndexGroups, gnames);
+ make_IMD_group(ir->imd, inputrecStrings->imd_grp, defaultIndexGroups, gnames);
}
gmx::IndexGroupsAndNames defaultIndexGroupsAndNames(
*defaultIndexGroups, gmx::arrayRefFromArray(gnames, defaultIndexGroups->nr));
- notifier.notifier_.notify(defaultIndexGroupsAndNames);
+ notifier.preProcessingNotifications_.notify(defaultIndexGroupsAndNames);
- auto accelerations = gmx::splitString(is->acc);
- auto accelerationGroupNames = gmx::splitString(is->accgrps);
+ auto accelerations = gmx::splitString(inputrecStrings->acc);
+ auto accelerationGroupNames = gmx::splitString(inputrecStrings->accgrps);
if (accelerationGroupNames.size() * DIM != accelerations.size())
{
gmx_fatal(FARGS, "Invalid Acceleration input: %zu groups and %zu acc. values",
convertRvecs(wi, accelerations, "anneal-time", ir->opts.acc);
- auto freezeDims = gmx::splitString(is->frdim);
- auto freezeGroupNames = gmx::splitString(is->freeze);
+ auto freezeDims = gmx::splitString(inputrecStrings->frdim);
+ auto freezeGroupNames = gmx::splitString(inputrecStrings->freeze);
if (freezeDims.size() != DIM * freezeGroupNames.size())
{
gmx_fatal(FARGS, "Invalid Freezing input: %zu groups and %zu freeze values",
}
}
- auto energyGroupNames = gmx::splitString(is->energy);
+ auto energyGroupNames = gmx::splitString(inputrecStrings->energy);
do_numbering(natoms, groups, energyGroupNames, defaultIndexGroups, gnames,
SimulationAtomGroupType::EnergyOutput, restnm, egrptpALL_GENREST, bVerbose, wi);
add_wall_energrps(groups, ir->nwall, symtab);
ir->opts.ngener = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
- auto vcmGroupNames = gmx::splitString(is->vcm);
+ auto vcmGroupNames = gmx::splitString(inputrecStrings->vcm);
do_numbering(natoms, groups, vcmGroupNames, defaultIndexGroups, gnames,
SimulationAtomGroupType::MassCenterVelocityRemoval, restnm,
vcmGroupNames.empty() ? egrptpALL_GENREST : egrptpPART, bVerbose, wi);
/* Now we have filled the freeze struct, so we can calculate NRDF */
calc_nrdf(mtop, ir, gnames);
- auto user1GroupNames = gmx::splitString(is->user1);
+ auto user1GroupNames = gmx::splitString(inputrecStrings->user1);
do_numbering(natoms, groups, user1GroupNames, defaultIndexGroups, gnames,
SimulationAtomGroupType::User1, restnm, egrptpALL_GENREST, bVerbose, wi);
- auto user2GroupNames = gmx::splitString(is->user2);
+ auto user2GroupNames = gmx::splitString(inputrecStrings->user2);
do_numbering(natoms, groups, user2GroupNames, defaultIndexGroups, gnames,
SimulationAtomGroupType::User2, restnm, egrptpALL_GENREST, bVerbose, wi);
- auto compressedXGroupNames = gmx::splitString(is->x_compressed_groups);
+ auto compressedXGroupNames = gmx::splitString(inputrecStrings->x_compressed_groups);
do_numbering(natoms, groups, compressedXGroupNames, defaultIndexGroups, gnames,
SimulationAtomGroupType::CompressedPositionOutput, restnm, egrptpONE, bVerbose, wi);
- auto orirefFitGroupNames = gmx::splitString(is->orirefitgrp);
+ auto orirefFitGroupNames = gmx::splitString(inputrecStrings->orirefitgrp);
do_numbering(natoms, groups, orirefFitGroupNames, defaultIndexGroups, gnames,
SimulationAtomGroupType::OrientationRestraintsFit, restnm, egrptpALL_GENREST,
bVerbose, wi);
- /* QMMM input processing */
- auto qmGroupNames = gmx::splitString(is->QMMM);
- auto qmMethods = gmx::splitString(is->QMmethod);
- auto qmBasisSets = gmx::splitString(is->QMbasis);
- if (ir->eI != eiMimic)
- {
- if (qmMethods.size() != qmGroupNames.size() || qmBasisSets.size() != qmGroupNames.size())
- {
- gmx_fatal(FARGS,
- "Invalid QMMM input: %zu groups %zu basissets"
- " and %zu methods\n",
- qmGroupNames.size(), qmBasisSets.size(), qmMethods.size());
- }
- /* group rest, if any, is always MM! */
- do_numbering(natoms, groups, qmGroupNames, defaultIndexGroups, gnames,
- SimulationAtomGroupType::QuantumMechanics, restnm, egrptpALL_GENREST, bVerbose, wi);
- nr = qmGroupNames.size(); /*atoms->grps[egcQMMM].nr;*/
- ir->opts.ngQM = qmGroupNames.size();
- snew(ir->opts.QMmethod, nr);
- snew(ir->opts.QMbasis, nr);
- for (i = 0; i < nr; i++)
- {
- /* input consists of strings: RHF CASSCF PM3 .. These need to be
- * converted to the corresponding enum in names.c
- */
- ir->opts.QMmethod[i] = search_QMstring(qmMethods[i].c_str(), eQMmethodNR, eQMmethod_names);
- ir->opts.QMbasis[i] = search_QMstring(qmBasisSets[i].c_str(), eQMbasisNR, eQMbasis_names);
- }
- auto qmMultiplicities = gmx::splitString(is->QMmult);
- auto qmCharges = gmx::splitString(is->QMcharge);
- auto qmbSH = gmx::splitString(is->bSH);
- snew(ir->opts.QMmult, nr);
- snew(ir->opts.QMcharge, nr);
- snew(ir->opts.bSH, nr);
- convertInts(wi, qmMultiplicities, "QMmult", ir->opts.QMmult);
- convertInts(wi, qmCharges, "QMcharge", ir->opts.QMcharge);
- convertYesNos(wi, qmbSH, "bSH", ir->opts.bSH);
-
- auto CASelectrons = gmx::splitString(is->CASelectrons);
- auto CASorbitals = gmx::splitString(is->CASorbitals);
- snew(ir->opts.CASelectrons, nr);
- snew(ir->opts.CASorbitals, nr);
- convertInts(wi, CASelectrons, "CASelectrons", ir->opts.CASelectrons);
- convertInts(wi, CASorbitals, "CASOrbitals", ir->opts.CASorbitals);
-
- auto SAon = gmx::splitString(is->SAon);
- auto SAoff = gmx::splitString(is->SAoff);
- auto SAsteps = gmx::splitString(is->SAsteps);
- snew(ir->opts.SAon, nr);
- snew(ir->opts.SAoff, nr);
- snew(ir->opts.SAsteps, nr);
- convertInts(wi, SAon, "SAon", ir->opts.SAon);
- convertInts(wi, SAoff, "SAoff", ir->opts.SAoff);
- convertInts(wi, SAsteps, "SAsteps", ir->opts.SAsteps);
- }
- else
+ /* MiMiC QMMM input processing */
+ auto qmGroupNames = gmx::splitString(inputrecStrings->QMMM);
+ if (qmGroupNames.size() > 1)
{
- /* MiMiC */
- if (qmGroupNames.size() > 1)
- {
- gmx_fatal(FARGS, "Currently, having more than one QM group in MiMiC is not supported");
- }
- /* group rest, if any, is always MM! */
- do_numbering(natoms, groups, qmGroupNames, defaultIndexGroups, gnames,
- SimulationAtomGroupType::QuantumMechanics, restnm, egrptpALL_GENREST, bVerbose, wi);
-
- ir->opts.ngQM = qmGroupNames.size();
+ gmx_fatal(FARGS, "Currently, having more than one QM group in MiMiC is not supported");
}
+ /* group rest, if any, is always MM! */
+ do_numbering(natoms, groups, qmGroupNames, defaultIndexGroups, gnames,
+ SimulationAtomGroupType::QuantumMechanics, restnm, egrptpALL_GENREST, bVerbose, wi);
+ ir->opts.ngQM = qmGroupNames.size();
- /* end of QMMM input */
+ /* end of MiMiC QMMM input */
if (bVerbose)
{
nr = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
snew(ir->opts.egp_flags, nr * nr);
- bExcl = do_egp_flag(ir, groups, "energygrp-excl", is->egpexcl, EGP_EXCL);
+ bExcl = do_egp_flag(ir, groups, "energygrp-excl", inputrecStrings->egpexcl, EGP_EXCL);
if (bExcl && ir->cutoff_scheme == ecutsVERLET)
{
warning_error(wi, "Energy group exclusions are currently not supported");
warning(wi, "Can not exclude the lattice Coulomb energy between energy groups");
}
- bTable = do_egp_flag(ir, groups, "energygrp-table", is->egptable, EGP_TABLE);
+ bTable = do_egp_flag(ir, groups, "energygrp-table", inputrecStrings->egptable, EGP_TABLE);
if (bTable && !(ir->vdwtype == evdwUSER) && !(ir->coulombtype == eelUSER)
&& !(ir->coulombtype == eelPMEUSER) && !(ir->coulombtype == eelPMEUSERSWITCH))
{
void triple_check(const char* mdparin, t_inputrec* ir, gmx_mtop_t* sys, warninp_t wi)
{
+ // Not meeting MTS requirements should have resulted in a fatal error, so we can assert here
+ gmx::assertMtsRequirements(*ir);
+
char err_buf[STRLEN];
int i, m, c, nmol;
bool bCharge, bAcc;
char warn_buf[STRLEN];
const char* ptr;
- ptr = check_box(ir->ePBC, box);
+ ptr = check_box(ir->pbcType, box);
if (ptr)
{
warning_error(wi, ptr);
ir->LincsWarnAngle = 90.0;
}
- if (ir->ePBC != epbcNONE)
+ if (ir->pbcType != PbcType::No)
{
if (ir->nstlist == 0)
{
"With nstlist=0 atoms are only put into the box at step 0, therefore drifting "
"atoms might cause the simulation to crash.");
}
- if (gmx::square(ir->rlist) >= max_cutoff2(ir->ePBC, box))
+ if (gmx::square(ir->rlist) >= max_cutoff2(ir->pbcType, box))
{
sprintf(warn_buf,
"ERROR: The cut-off length is longer than half the shortest box vector or "
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2011-2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2011-2019,2020,2021, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
#include "gromacs/domdec/localatomsetmanager.h"
#include "gromacs/domdec/partition.h"
#include "gromacs/ewald/ewald_utils.h"
-#include "gromacs/ewald/pme.h"
#include "gromacs/ewald/pme_gpu_program.h"
+#include "gromacs/ewald/pme_only.h"
#include "gromacs/ewald/pme_pp_comm_gpu.h"
#include "gromacs/fileio/checkpoint.h"
#include "gromacs/fileio/gmxfio.h"
#include "gromacs/fileio/tpxio.h"
#include "gromacs/gmxlib/network.h"
#include "gromacs/gmxlib/nrnb.h"
-#include "gromacs/gpu_utils/gpu_utils.h"
+#include "gromacs/gpu_utils/device_stream_manager.h"
#include "gromacs/hardware/cpuinfo.h"
#include "gromacs/hardware/detecthardware.h"
+#include "gromacs/hardware/device_management.h"
+#include "gromacs/hardware/hardwaretopology.h"
#include "gromacs/hardware/printhardware.h"
#include "gromacs/imd/imd.h"
#include "gromacs/listed_forces/disre.h"
#include "gromacs/listed_forces/gpubonded.h"
+#include "gromacs/listed_forces/listed_forces.h"
#include "gromacs/listed_forces/orires.h"
#include "gromacs/math/functions.h"
#include "gromacs/math/utilities.h"
#include "gromacs/mdlib/force.h"
#include "gromacs/mdlib/forcerec.h"
#include "gromacs/mdlib/gmx_omp_nthreads.h"
+#include "gromacs/mdlib/gpuforcereduction.h"
#include "gromacs/mdlib/makeconstraints.h"
#include "gromacs/mdlib/md_support.h"
#include "gromacs/mdlib/mdatoms.h"
-#include "gromacs/mdlib/membed.h"
-#include "gromacs/mdlib/qmmm.h"
#include "gromacs/mdlib/sighandler.h"
#include "gromacs/mdlib/stophandler.h"
+#include "gromacs/mdlib/tgroup.h"
#include "gromacs/mdlib/updategroups.h"
+#include "gromacs/mdlib/vsite.h"
#include "gromacs/mdrun/mdmodules.h"
#include "gromacs/mdrun/simulationcontext.h"
+#include "gromacs/mdrun/simulationinput.h"
+#include "gromacs/mdrun/simulationinputhandle.h"
#include "gromacs/mdrunutility/handlerestart.h"
#include "gromacs/mdrunutility/logging.h"
#include "gromacs/mdrunutility/multisim.h"
#include "gromacs/mdrunutility/printtime.h"
#include "gromacs/mdrunutility/threadaffinity.h"
+#include "gromacs/mdtypes/checkpointdata.h"
#include "gromacs/mdtypes/commrec.h"
#include "gromacs/mdtypes/enerdata.h"
#include "gromacs/mdtypes/fcdata.h"
+#include "gromacs/mdtypes/forcerec.h"
#include "gromacs/mdtypes/group.h"
#include "gromacs/mdtypes/inputrec.h"
+#include "gromacs/mdtypes/interaction_const.h"
#include "gromacs/mdtypes/md_enums.h"
+#include "gromacs/mdtypes/mdatom.h"
#include "gromacs/mdtypes/mdrunoptions.h"
#include "gromacs/mdtypes/observableshistory.h"
#include "gromacs/mdtypes/simulation_workload.h"
#include "gromacs/mdtypes/state.h"
#include "gromacs/mdtypes/state_propagator_data_gpu.h"
+#include "gromacs/modularsimulator/modularsimulator.h"
#include "gromacs/nbnxm/gpu_data_mgmt.h"
#include "gromacs/nbnxm/nbnxm.h"
#include "gromacs/nbnxm/pairlist_tuning.h"
#include "gromacs/utility/stringutil.h"
#include "isimulator.h"
+#include "membedholder.h"
#include "replicaexchange.h"
#include "simulatorbuilder.h"
// getenv results are ignored when clearly they are used.
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wunused-result"
- devFlags.enableGpuBufferOps = (getenv("GMX_USE_GPU_BUFFER_OPS") != nullptr)
- && (GMX_GPU == GMX_GPU_CUDA) && useGpuForNonbonded;
+
+ devFlags.enableGpuBufferOps =
+ GMX_GPU_CUDA && useGpuForNonbonded && (getenv("GMX_USE_GPU_BUFFER_OPS") != nullptr);
+ devFlags.enableGpuHaloExchange = GMX_GPU_CUDA && GMX_THREAD_MPI && getenv("GMX_GPU_DD_COMMS") != nullptr;
devFlags.forceGpuUpdateDefault = (getenv("GMX_FORCE_UPDATE_DEFAULT_GPU") != nullptr) || GMX_FAHCORE;
- devFlags.enableGpuHaloExchange =
- (getenv("GMX_GPU_DD_COMMS") != nullptr && GMX_THREAD_MPI && (GMX_GPU == GMX_GPU_CUDA));
devFlags.enableGpuPmePPComm =
- (getenv("GMX_GPU_PME_PP_COMMS") != nullptr && GMX_THREAD_MPI && (GMX_GPU == GMX_GPU_CUDA));
+ GMX_GPU_CUDA && GMX_THREAD_MPI && getenv("GMX_GPU_PME_PP_COMMS") != nullptr;
+
#pragma GCC diagnostic pop
if (devFlags.enableGpuBufferOps)
GMX_LOG(mdlog.warning)
.asParagraph()
.appendTextFormatted(
- "This run uses the 'GPU halo exchange' feature, enabled by the "
+ "This run has requested the 'GPU halo exchange' feature, enabled by "
+ "the "
"GMX_GPU_DD_COMMS environment variable.");
}
else
{
if (pmeRunMode == PmeRunMode::GPU)
{
+ if (!devFlags.enableGpuBufferOps)
+ {
+ GMX_LOG(mdlog.warning)
+ .asParagraph()
+ .appendTextFormatted(
+ "Enabling GPU buffer operations required by GMX_GPU_PME_PP_COMMS "
+ "(equivalent with GMX_USE_GPU_BUFFER_OPS=1).");
+ devFlags.enableGpuBufferOps = true;
+ }
GMX_LOG(mdlog.warning)
.asParagraph()
.appendTextFormatted(
// Copy members of master runner.
// \todo Replace with builder when Simulation context and/or runner phases are better defined.
- // Ref https://redmine.gromacs.org/issues/2587 and https://redmine.gromacs.org/issues/2375
+ // Ref https://gitlab.com/gromacs/gromacs/-/issues/2587 and https://gitlab.com/gromacs/gromacs/-/issues/2375
newRunner.hw_opt = hw_opt;
newRunner.filenames = filenames;
+ newRunner.hwinfo_ = hwinfo_;
newRunner.oenv = oenv;
newRunner.mdrunOptions = mdrunOptions;
newRunner.domdecOptions = domdecOptions;
newRunner.pforce = pforce;
// Give the spawned thread the newly created valid communicator
// for the simulation.
- newRunner.communicator = MPI_COMM_WORLD;
- newRunner.ms = ms;
- newRunner.startingBehavior = startingBehavior;
- newRunner.stopHandlerBuilder_ = std::make_unique<StopHandlerBuilder>(*stopHandlerBuilder_);
+ newRunner.libraryWorldCommunicator = MPI_COMM_WORLD;
+ newRunner.simulationCommunicator = MPI_COMM_WORLD;
+ newRunner.ms = ms;
+ newRunner.startingBehavior = startingBehavior;
+ newRunner.stopHandlerBuilder_ = std::make_unique<StopHandlerBuilder>(*stopHandlerBuilder_);
+ newRunner.inputHolder_ = inputHolder_;
threadMpiMdrunnerAccessBarrier();
// Give the master thread the newly created valid communicator for
// the simulation.
- communicator = MPI_COMM_WORLD;
+ libraryWorldCommunicator = MPI_COMM_WORLD;
+ simulationCommunicator = MPI_COMM_WORLD;
threadMpiMdrunnerAccessBarrier();
#else
GMX_UNUSED_VALUE(numThreadsToLaunch);
bool makeGpuPairList,
const gmx::CpuInfo& cpuinfo)
{
+ // We checked the cut-offs in grompp, but double-check here.
+ // We have PME+LJcutoff kernels for rcoulomb>rvdw.
+ if (EEL_PME_EWALD(ir->coulombtype) && ir->vdwtype == eelCUT)
+ {
+ GMX_RELEASE_ASSERT(ir->rcoulomb >= ir->rvdw,
+ "With Verlet lists and PME we should have rcoulomb>=rvdw");
+ }
+ else
+ {
+ GMX_RELEASE_ASSERT(ir->rcoulomb == ir->rvdw,
+ "With Verlet lists and no PME rcoulomb and rvdw should be identical");
+ }
/* For NVE simulations, we will retain the initial list buffer */
if (EI_DYNAMICS(ir->eI) && ir->verletbuf_tol > 0 && !(EI_MD(ir->eI) && ir->etc == etcNO))
{
{
matrix box;
t_forcerec* fr = nullptr;
- t_fcdata* fcd = nullptr;
real ewaldcoeff_q = 0;
real ewaldcoeff_lj = 0;
int nChargePerturbed = -1, nTypePerturbed = 0;
gmx_wallcycle_t wcycle;
gmx_walltime_accounting_t walltime_accounting = nullptr;
- gmx_membed_t* membed = nullptr;
- gmx_hw_info_t* hwinfo = nullptr;
+ MembedHolder membedHolder(filenames.size(), filenames.data());
/* CAUTION: threads may be started later on in this function, so
cr doesn't reflect the final parallel state right now */
/* TODO: inputrec should tell us whether we use an algorithm, not a file option */
const bool doEssentialDynamics = opt2bSet("-ei", filenames.size(), filenames.data());
- const bool doMembed = opt2bSet("-membed", filenames.size(), filenames.data());
const bool doRerun = mdrunOptions.rerun;
// Handle task-assignment related user options.
{
fplog = gmx_fio_getfp(logFileHandle);
}
- const bool isSimulationMasterRank = findIsSimulationMasterRank(ms, communicator);
+ const bool isSimulationMasterRank = findIsSimulationMasterRank(ms, simulationCommunicator);
gmx::LoggerOwner logOwner(buildLogger(fplog, isSimulationMasterRank));
gmx::MDLogger mdlog(logOwner.logger());
- // TODO The thread-MPI master rank makes a working
- // PhysicalNodeCommunicator here, but it gets rebuilt by all ranks
- // after the threads have been launched. This works because no use
- // is made of that communicator until after the execution paths
- // have rejoined. But it is likely that we can improve the way
- // this is expressed, e.g. by expressly running detection only the
- // master rank for thread-MPI, rather than relying on the mutex
- // and reference count.
- PhysicalNodeCommunicator physicalNodeComm(communicator, gmx_physicalnode_id_hash());
- hwinfo = gmx_detect_hardware(mdlog, physicalNodeComm);
-
- gmx_print_detected_hardware(fplog, isSimulationMasterRank && isMasterSim(ms), mdlog, hwinfo);
+ gmx_print_detected_hardware(fplog, isSimulationMasterRank && isMasterSim(ms), mdlog, hwinfo_);
- std::vector<int> gpuIdsToUse = makeGpuIdsToUse(hwinfo->gpu_info, hw_opt.gpuIdsAvailable);
+ std::vector<int> gpuIdsToUse = makeGpuIdsToUse(hwinfo_->deviceInfoList, hw_opt.gpuIdsAvailable);
+ const int numDevicesToUse = gmx::ssize(gpuIdsToUse);
// Print citation requests after all software/hardware printing
pleaseCiteGromacs(fplog);
- // TODO Replace this by unique_ptr once t_inputrec is C++
- t_inputrec inputrecInstance;
- t_inputrec* inputrec = nullptr;
- std::unique_ptr<t_state> globalState;
+ // Note: legacy program logic relies on checking whether these pointers are assigned.
+ // Objects may or may not be allocated later.
+ std::unique_ptr<t_inputrec> inputrec;
+ std::unique_ptr<t_state> globalState;
auto partialDeserializedTpr = std::make_unique<PartialDeserializedTprFile>();
if (isSimulationMasterRank)
{
+ // Allocate objects to be initialized by later function calls.
/* Only the master rank has the global state */
globalState = std::make_unique<t_state>();
+ inputrec = std::make_unique<t_inputrec>();
/* Read (nearly) all data required for the simulation
* and keep the partly serialized tpr contents to send to other ranks later
*/
- *partialDeserializedTpr = read_tpx_state(ftp2fn(efTPR, filenames.size(), filenames.data()),
- &inputrecInstance, globalState.get(), &mtop);
- inputrec = &inputrecInstance;
+ applyGlobalSimulationState(*inputHolder_.get(), partialDeserializedTpr.get(),
+ globalState.get(), inputrec.get(), &mtop);
}
/* Check and update the hardware options for internal consistency */
checkAndUpdateHardwareOptions(mdlog, &hw_opt, isSimulationMasterRank, domdecOptions.numPmeRanks,
- inputrec);
+ inputrec.get());
if (GMX_THREAD_MPI && isSimulationMasterRank)
{
// the number of GPUs to choose the number of ranks.
auto canUseGpuForNonbonded = buildSupportsNonbondedOnGpu(nullptr);
useGpuForNonbonded = decideWhetherToUseGpusForNonbondedWithThreadMpi(
- nonbondedTarget, gpuIdsToUse, userGpuTaskAssignment, emulateGpuNonbonded,
+ nonbondedTarget, numDevicesToUse, userGpuTaskAssignment, emulateGpuNonbonded,
canUseGpuForNonbonded,
gpuAccelerationOfNonbondedIsUseful(mdlog, *inputrec, GMX_THREAD_MPI),
hw_opt.nthreads_tmpi);
useGpuForPme = decideWhetherToUseGpusForPmeWithThreadMpi(
- useGpuForNonbonded, pmeTarget, gpuIdsToUse, userGpuTaskAssignment, *hwinfo,
- *inputrec, mtop, hw_opt.nthreads_tmpi, domdecOptions.numPmeRanks);
+ useGpuForNonbonded, pmeTarget, numDevicesToUse, userGpuTaskAssignment, *hwinfo_,
+ *inputrec, hw_opt.nthreads_tmpi, domdecOptions.numPmeRanks);
}
GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
* TODO Over-writing the user-supplied value here does
* prevent any possible subsequent checks from working
* correctly. */
- hw_opt.nthreads_tmpi = get_nthreads_mpi(hwinfo, &hw_opt, gpuIdsToUse, useGpuForNonbonded,
- useGpuForPme, inputrec, &mtop, mdlog, doMembed);
+ hw_opt.nthreads_tmpi =
+ get_nthreads_mpi(hwinfo_, &hw_opt, numDevicesToUse, useGpuForNonbonded, useGpuForPme,
+ inputrec.get(), &mtop, mdlog, membedHolder.doMembed());
// Now start the threads for thread MPI.
spawnThreads(hw_opt.nthreads_tmpi);
// The spawned threads enter mdrunner() and execution of
// master and spawned threads joins at the end of this block.
- physicalNodeComm = PhysicalNodeCommunicator(communicator, gmx_physicalnode_id_hash());
}
- GMX_RELEASE_ASSERT(communicator == MPI_COMM_WORLD, "Must have valid world communicator");
- CommrecHandle crHandle = init_commrec(communicator, ms);
+ GMX_RELEASE_ASSERT(ms || simulationCommunicator != MPI_COMM_NULL,
+ "Must have valid communicator unless running a multi-simulation");
+ CommrecHandle crHandle = init_commrec(simulationCommunicator);
t_commrec* cr = crHandle.get();
GMX_RELEASE_ASSERT(cr != nullptr, "Must have valid commrec");
+ PhysicalNodeCommunicator physicalNodeComm(libraryWorldCommunicator, gmx_physicalnode_id_hash());
+
+ // If we detected the topology on this system, double-check that it makes sense
+ if (hwinfo_->hardwareTopology->isThisSystem())
+ {
+ hardwareTopologyDoubleCheckDetection(mdlog, *hwinfo_->hardwareTopology);
+ }
+
if (PAR(cr))
{
/* now broadcast everything to the non-master nodes/threads: */
if (!isSimulationMasterRank)
{
- inputrec = &inputrecInstance;
+ // Until now, only the master rank has a non-null pointer.
+ // On non-master ranks, allocate the object that will receive data in the following call.
+ inputrec = std::make_unique<t_inputrec>();
}
- init_parallel(cr, inputrec, &mtop, partialDeserializedTpr.get());
+ init_parallel(cr->mpiDefaultCommunicator, MASTER(cr), inputrec.get(), &mtop,
+ partialDeserializedTpr.get());
}
GMX_RELEASE_ASSERT(inputrec != nullptr, "All ranks should have a valid inputrec now");
partialDeserializedTpr.reset(nullptr);
// the inputrec read by the master rank. The ranks can now all run
// the task-deciding functions and will agree on the result
// without needing to communicate.
- //
- // TODO Should we do the communication in debug mode to support
- // having an assertion?
const bool useDomainDecomposition = (PAR(cr) && !(EI_TPI(inputrec->eI) || inputrec->eI == eiNM));
// Note that these variables describe only their own node.
bool useGpuForPme = false;
bool useGpuForBonded = false;
bool useGpuForUpdate = false;
- bool gpusWereDetected = hwinfo->ngpu_compatible_tot > 0;
+ bool gpusWereDetected = hwinfo_->ngpu_compatible_tot > 0;
try
{
// It's possible that there are different numbers of GPUs on
nonbondedTarget, userGpuTaskAssignment, emulateGpuNonbonded, canUseGpuForNonbonded,
gpuAccelerationOfNonbondedIsUseful(mdlog, *inputrec, !GMX_THREAD_MPI), gpusWereDetected);
useGpuForPme = decideWhetherToUseGpusForPme(
- useGpuForNonbonded, pmeTarget, userGpuTaskAssignment, *hwinfo, *inputrec, mtop,
- cr->nnodes, domdecOptions.numPmeRanks, gpusWereDetected);
- auto canUseGpuForBonded = buildSupportsGpuBondeds(nullptr)
- && inputSupportsGpuBondeds(*inputrec, mtop, nullptr);
- useGpuForBonded = decideWhetherToUseGpusForBonded(
- useGpuForNonbonded, useGpuForPme, bondedTarget, canUseGpuForBonded,
- EVDW_PME(inputrec->vdwtype), EEL_PME_EWALD(inputrec->coulombtype),
- domdecOptions.numPmeRanks, gpusWereDetected);
+ useGpuForNonbonded, pmeTarget, userGpuTaskAssignment, *hwinfo_, *inputrec,
+ cr->sizeOfDefaultCommunicator, domdecOptions.numPmeRanks, gpusWereDetected);
+ useGpuForBonded = decideWhetherToUseGpusForBonded(useGpuForNonbonded, useGpuForPme,
+ bondedTarget, *inputrec, mtop,
+ domdecOptions.numPmeRanks, gpusWereDetected);
}
GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
const DevelopmentFeatureFlags devFlags =
manageDevelopmentFeatures(mdlog, useGpuForNonbonded, pmeRunMode);
- const bool inputIsCompatibleWithModularSimulator = ModularSimulator::isInputCompatible(
- false, inputrec, doRerun, mtop, ms, replExParams, nullptr, doEssentialDynamics, doMembed);
- const bool useModularSimulator = inputIsCompatibleWithModularSimulator
- && !(getenv("GMX_DISABLE_MODULAR_SIMULATOR") != nullptr);
+ const bool useModularSimulator =
+ checkUseModularSimulator(false, inputrec.get(), doRerun, mtop, ms, replExParams,
+ nullptr, doEssentialDynamics, membedHolder.doMembed());
// Build restraints.
// TODO: hide restraint implementation details from Mdrunner.
// TODO: Error handling
mdModules_->assignOptionsToModules(*inputrec->params, nullptr);
- const auto& mdModulesNotifier = mdModules_->notifier().notifier_;
+ // now that the MdModules know their options, they know which callbacks to sign up to
+ mdModules_->subscribeToSimulationSetupNotifications();
+ const auto& mdModulesNotifier = mdModules_->notifier().simulationSetupNotifications_;
if (inputrec->internalParameters != nullptr)
{
if (fplog != nullptr)
{
- pr_inputrec(fplog, 0, "Input Parameters", inputrec, FALSE);
+ pr_inputrec(fplog, 0, "Input Parameters", inputrec.get(), FALSE);
fprintf(fplog, "\n");
}
}
/* now make sure the state is initialized and propagated */
- set_state_entries(globalState.get(), inputrec, useModularSimulator);
+ set_state_entries(globalState.get(), inputrec.get(), useModularSimulator);
}
/* NM and TPI parallelize over force/energy calculations, not atoms,
{
globalState = std::make_unique<t_state>();
}
- broadcastStateWithoutDynamics(cr, globalState.get());
+ broadcastStateWithoutDynamics(cr->mpiDefaultCommunicator, DOMAINDECOMP(cr), PAR(cr),
+ globalState.get());
}
/* A parallel command line option consistency check that we can
{
if (domdecOptions.numPmeRanks > 0)
{
- gmx_fatal_collective(FARGS, cr->mpi_comm_mysim, MASTER(cr),
+ gmx_fatal_collective(FARGS, cr->mpiDefaultCommunicator, MASTER(cr),
"PME-only ranks are requested, but the system does not use PME "
"for electrostatics or LJ");
}
* So the PME-only nodes (if present) will also initialize
* the distance restraints.
*/
- snew(fcd, 1);
/* This needs to be called before read_checkpoint to extend the state */
- init_disres(fplog, &mtop, inputrec, cr, ms, fcd, globalState.get(), replExParams.exchangeInterval > 0);
+ t_disresdata* disresdata;
+ snew(disresdata, 1);
+ init_disres(fplog, &mtop, inputrec.get(), DisResRunMode::MDRun,
+ MASTER(cr) ? DDRole::Master : DDRole::Agent,
+ PAR(cr) ? NumRanks::Multiple : NumRanks::Single, cr->mpi_comm_mysim, ms, disresdata,
+ globalState.get(), replExParams.exchangeInterval > 0);
- init_orires(fplog, &mtop, inputrec, cr, ms, globalState.get(), &(fcd->orires));
+ t_oriresdata* oriresdata;
+ snew(oriresdata, 1);
+ init_orires(fplog, &mtop, inputrec.get(), cr, ms, globalState.get(), oriresdata);
- auto deform = prepareBoxDeformation(globalState->box, cr, *inputrec);
+ auto deform = prepareBoxDeformation(
+ globalState != nullptr ? globalState->box : box, MASTER(cr) ? DDRole::Master : DDRole::Agent,
+ PAR(cr) ? NumRanks::Multiple : NumRanks::Single, cr->mpi_comm_mygroup, *inputrec);
#if GMX_FAHCORE
/* We have to remember the generation's first step before reading checkpoint.
ObservablesHistory observablesHistory = {};
+ auto modularSimulatorCheckpointData = std::make_unique<ReadCheckpointDataHolder>();
if (startingBehavior != StartingBehavior::NewSimulation)
{
/* Check if checkpoint file exists before doing continuation.
inputrec->nsteps = -1;
}
- load_checkpoint(opt2fn_master("-cpi", filenames.size(), filenames.data(), cr),
- logFileHandle, cr, domdecOptions.numCells, inputrec, globalState.get(),
- &observablesHistory, mdrunOptions.reproducible, mdModules_->notifier());
+ // Finish applying initial simulation state information from external sources on all ranks.
+ // Reconcile checkpoint file data with Mdrunner state established up to this point.
+ applyLocalState(*inputHolder_.get(), logFileHandle, cr, domdecOptions.numCells,
+ inputrec.get(), globalState.get(), &observablesHistory,
+ mdrunOptions.reproducible, mdModules_->notifier(),
+ modularSimulatorCheckpointData.get(), useModularSimulator);
+ // TODO: (#3652) Synchronize filesystem state, SimulationInput contents, and program
+ // invariants
+ // on all code paths.
+ // Write checkpoint or provide hook to update SimulationInput.
+ // If there was a checkpoint file, SimulationInput contains more information
+ // than if there wasn't. At this point, we have synchronized the in-memory
+ // state with the filesystem state only for restarted simulations. We should
+ // be calling applyLocalState unconditionally and expect that the completeness
+ // of SimulationInput is not dependent on its creation method.
if (startingBehavior == StartingBehavior::RestartWithAppending && logFileHandle)
{
"file field.");
}
/* override nsteps with value set on the commandline */
- override_nsteps_cmdline(mdlog, mdrunOptions.numStepsCommandline, inputrec);
+ override_nsteps_cmdline(mdlog, mdrunOptions.numStepsCommandline, inputrec.get());
- if (SIMMASTER(cr))
+ if (isSimulationMasterRank)
{
copy_mat(globalState->box, box);
}
if (PAR(cr))
{
- gmx_bcast(sizeof(box), box, cr);
+ gmx_bcast(sizeof(box), box, cr->mpiDefaultCommunicator);
}
if (inputrec->cutoff_scheme != ecutsVERLET)
"Verlet scheme, or use an earlier version of GROMACS if necessary.");
}
/* Update rlist and nstlist. */
- prepare_verlet_scheme(fplog, cr, inputrec, nstlist_cmdline, &mtop, box,
+ /* Note: prepare_verlet_scheme is calling increaseNstlist(...), which (while attempting to
+ * increase rlist) tries to check if the newly chosen value fits with the DD scheme. As this is
+ * run before any DD scheme is set up, this check is never executed. See #3334 for more details.
+ */
+ prepare_verlet_scheme(fplog, cr, inputrec.get(), nstlist_cmdline, &mtop, box,
useGpuForNonbonded || (emulateGpuNonbonded == EmulateGpuNonbonded::Yes),
- *hwinfo->cpuInfo);
+ *hwinfo_->cpuInfo);
- const bool prefer1DAnd1PulseDD = (devFlags.enableGpuHaloExchange && useGpuForNonbonded);
// This builder is necessary while we have multi-part construction
// of DD. Before DD is constructed, we use the existence of
// the builder object to indicate that further construction of DD
if (useDomainDecomposition)
{
ddBuilder = std::make_unique<DomainDecompositionBuilder>(
- mdlog, cr, domdecOptions, mdrunOptions, prefer1DAnd1PulseDD, mtop, *inputrec, box,
+ mdlog, cr, domdecOptions, mdrunOptions, mtop, *inputrec, box,
positionsFromStatePointer(globalState.get()));
}
else
{
/* PME, if used, is done on all nodes with 1D decomposition */
- cr->npmenodes = 0;
- cr->duty = (DUTY_PP | DUTY_PME);
+ cr->nnodes = cr->sizeOfDefaultCommunicator;
+ cr->sim_nodeid = cr->rankInDefaultCommunicator;
+ cr->nodeid = cr->rankInDefaultCommunicator;
+ cr->npmenodes = 0;
+ cr->duty = (DUTY_PP | DUTY_PME);
- if (inputrec->ePBC == epbcSCREW)
+ if (inputrec->pbcType == PbcType::Screw)
{
gmx_fatal(FARGS, "pbc=screw is only implemented with domain decomposition");
}
}
- // Produce the task assignment for this rank.
- GpuTaskAssignmentsBuilder gpuTaskAssignmentsBuilder;
- GpuTaskAssignments gpuTaskAssignments = gpuTaskAssignmentsBuilder.build(
- gpuIdsToUse, userGpuTaskAssignment, *hwinfo, communicator, physicalNodeComm,
+ // Produce the task assignment for this rank - done after DD is constructed
+ GpuTaskAssignments gpuTaskAssignments = GpuTaskAssignmentsBuilder::build(
+ gpuIdsToUse, userGpuTaskAssignment, *hwinfo_, simulationCommunicator, physicalNodeComm,
nonbondedTarget, pmeTarget, bondedTarget, updateTarget, useGpuForNonbonded,
useGpuForPme, thisRankHasDuty(cr, DUTY_PP),
// TODO cr->duty & DUTY_PME should imply that a PME
EEL_PME(inputrec->coulombtype) && thisRankHasDuty(cr, DUTY_PME));
// Get the device handles for the modules, nullptr when no task is assigned.
- gmx_device_info_t* nonbondedDeviceInfo = gpuTaskAssignments.initNonbondedDevice(cr);
- gmx_device_info_t* pmeDeviceInfo = gpuTaskAssignments.initPmeDevice();
+ int deviceId = -1;
+ DeviceInformation* deviceInfo = gpuTaskAssignments.initDevice(&deviceId);
- // TODO Initialize GPU streams here.
+ // timing enabling - TODO put this in gpu_utils (even though generally this is just option handling?)
+ bool useTiming = true;
+
+ if (GMX_GPU_CUDA)
+ {
+ /* WARNING: CUDA timings are incorrect with multiple streams.
+ * This is the main reason why they are disabled by default.
+ */
+ // TODO: Consider turning on by default when we can detect nr of streams.
+ useTiming = (getenv("GMX_ENABLE_GPU_TIMING") != nullptr);
+ }
+ else if (GMX_GPU_OPENCL)
+ {
+ useTiming = (getenv("GMX_DISABLE_GPU_TIMING") == nullptr);
+ }
// TODO Currently this is always built, yet DD partition code
// checks if it is built before using it. Probably it should
const bool printHostName = (cr->nnodes > 1);
gpuTaskAssignments.reportGpuUsage(mdlog, printHostName, useGpuForBonded, pmeRunMode, useGpuForUpdate);
+ const bool disableNonbondedCalculation = (getenv("GMX_NO_NONBONDED") != nullptr);
+ if (disableNonbondedCalculation)
+ {
+ /* turn off non-bonded calculations */
+ GMX_LOG(mdlog.warning)
+ .asParagraph()
+ .appendText(
+ "Found environment variable GMX_NO_NONBONDED.\n"
+ "Disabling nonbonded calculations.");
+ }
+
+ MdrunScheduleWorkload runScheduleWork;
+
+ bool useGpuDirectHalo = decideWhetherToUseGpuForHalo(
+ devFlags, havePPDomainDecomposition(cr), useGpuForNonbonded, useModularSimulator,
+ doRerun, EI_ENERGY_MINIMIZATION(inputrec->eI));
+
+ // Also populates the simulation constant workload description.
+ runScheduleWork.simulationWork = createSimulationWorkload(
+ *inputrec, disableNonbondedCalculation, devFlags, useGpuForNonbonded, pmeRunMode,
+ useGpuForBonded, useGpuForUpdate, useGpuDirectHalo);
+
+ std::unique_ptr<DeviceStreamManager> deviceStreamManager = nullptr;
+
+ if (deviceInfo != nullptr)
+ {
+ if (DOMAINDECOMP(cr) && thisRankHasDuty(cr, DUTY_PP))
+ {
+ dd_setup_dlb_resource_sharing(cr, deviceId);
+ }
+ deviceStreamManager = std::make_unique<DeviceStreamManager>(
+ *deviceInfo, havePPDomainDecomposition(cr), runScheduleWork.simulationWork, useTiming);
+ }
+
// If the user chose a task assignment, give them some hints
// where appropriate.
if (!userGpuTaskAssignment.empty())
{
- gpuTaskAssignments.logPerformanceHints(mdlog, ssize(gpuIdsToUse));
+ gpuTaskAssignments.logPerformanceHints(mdlog, numDevicesToUse);
}
if (PAR(cr))
.appendTextFormatted(
"This is simulation %d out of %d running as a composite GROMACS\n"
"multi-simulation job. Setup for this simulation:\n",
- ms->sim, ms->nsim);
+ ms->simulationIndex_, ms->numSimulations_);
}
GMX_LOG(mdlog.warning)
.appendTextFormatted("Using %d MPI %s\n", cr->nnodes,
// that existing affinity setting was from OpenMP or something
// else, so we run this code both before and after we initialize
// the OpenMP support.
- gmx_check_thread_affinity_set(mdlog, &hw_opt, hwinfo->nthreads_hw_avail, FALSE);
+ gmx_check_thread_affinity_set(mdlog, &hw_opt, hwinfo_->nthreads_hw_avail, FALSE);
/* Check and update the number of OpenMP threads requested */
- checkAndUpdateRequestedNumOpenmpThreads(&hw_opt, *hwinfo, cr, ms, physicalNodeComm.size_,
+ checkAndUpdateRequestedNumOpenmpThreads(&hw_opt, *hwinfo_, cr, ms, physicalNodeComm.size_,
pmeRunMode, mtop, *inputrec);
- gmx_omp_nthreads_init(mdlog, cr, hwinfo->nthreads_hw_avail, physicalNodeComm.size_,
+ gmx_omp_nthreads_init(mdlog, cr, hwinfo_->nthreads_hw_avail, physicalNodeComm.size_,
hw_opt.nthreads_omp, hw_opt.nthreads_omp_pme, !thisRankHasDuty(cr, DUTY_PP));
// Enable FP exception detection, but not in
}
/* Now that we know the setup is consistent, check for efficiency */
- check_resource_division_efficiency(hwinfo, gpuTaskAssignments.thisRankHasAnyGpuTask(),
+ check_resource_division_efficiency(hwinfo_, gpuTaskAssignments.thisRankHasAnyGpuTask(),
mdrunOptions.ntompOptionIsSet, cr, mdlog);
/* getting number of PP/PME threads on this MPI / tMPI rank.
*/
const int numThreadsOnThisRank = thisRankHasDuty(cr, DUTY_PP) ? gmx_omp_nthreads_get(emntNonbonded)
: gmx_omp_nthreads_get(emntPME);
- checkHardwareOversubscription(numThreadsOnThisRank, cr->nodeid, *hwinfo->hardwareTopology,
+ checkHardwareOversubscription(numThreadsOnThisRank, cr->nodeid, *hwinfo_->hardwareTopology,
physicalNodeComm, mdlog);
// Enable Peer access between GPUs where available
// Only for DD, only master PP rank needs to perform setup, and only if thread MPI plus
// any of the GPU communication features are active.
if (DOMAINDECOMP(cr) && MASTER(cr) && thisRankHasDuty(cr, DUTY_PP) && GMX_THREAD_MPI
- && (devFlags.enableGpuHaloExchange || devFlags.enableGpuPmePPComm))
+ && (runScheduleWork.simulationWork.useGpuHaloExchange
+ || runScheduleWork.simulationWork.useGpuPmePpCommunication))
{
setupGpuDevicePeerAccess(gpuIdsToUse, mdlog);
}
* - which indicates that probably the OpenMP library has changed it
* since we first checked).
*/
- gmx_check_thread_affinity_set(mdlog, &hw_opt, hwinfo->nthreads_hw_avail, TRUE);
+ gmx_check_thread_affinity_set(mdlog, &hw_opt, hwinfo_->nthreads_hw_avail, TRUE);
int numThreadsOnThisNode, intraNodeThreadOffset;
analyzeThreadsOnThisNode(physicalNodeComm, numThreadsOnThisRank, &numThreadsOnThisNode,
&intraNodeThreadOffset);
/* Set the CPU affinity */
- gmx_set_thread_affinity(mdlog, cr, &hw_opt, *hwinfo->hardwareTopology, numThreadsOnThisRank,
+ gmx_set_thread_affinity(mdlog, cr, &hw_opt, *hwinfo_->hardwareTopology, numThreadsOnThisRank,
numThreadsOnThisNode, intraNodeThreadOffset, nullptr);
}
/* Master synchronizes its value of reset_counters with all nodes
* including PME only nodes */
int64_t reset_counters = wcycle_get_reset_counters(wcycle);
- gmx_bcast_sim(sizeof(reset_counters), &reset_counters, cr);
+ gmx_bcast(sizeof(reset_counters), &reset_counters, cr->mpi_comm_mysim);
wcycle_set_reset_counters(wcycle, reset_counters);
}
// Membrane embedding must be initialized before we call init_forcerec()
- if (doMembed)
- {
- if (MASTER(cr))
- {
- fprintf(stderr, "Initializing membed");
- }
- /* Note that membed cannot work in parallel because mtop is
- * changed here. Fix this if we ever want to make it run with
- * multiple ranks. */
- membed = init_membed(fplog, filenames.size(), filenames.data(), &mtop, inputrec,
- globalState.get(), cr, &mdrunOptions.checkpointOptions.period);
- }
+ membedHolder.initializeMembed(fplog, filenames.size(), filenames.data(), &mtop, inputrec.get(),
+ globalState.get(), cr, &mdrunOptions.checkpointOptions.period);
- const bool thisRankHasPmeGpuTask = gpuTaskAssignments.thisRankHasPmeGpuTask();
- std::unique_ptr<MDAtoms> mdAtoms;
- std::unique_ptr<gmx_vsite_t> vsite;
+ const bool thisRankHasPmeGpuTask = gpuTaskAssignments.thisRankHasPmeGpuTask();
+ std::unique_ptr<MDAtoms> mdAtoms;
+ std::unique_ptr<VirtualSitesHandler> vsite;
+ std::unique_ptr<GpuBonded> gpuBonded;
t_nrnb nrnb;
if (thisRankHasDuty(cr, DUTY_PP))
{
mdModulesNotifier.notify(*cr);
mdModulesNotifier.notify(&atomSets);
- mdModulesNotifier.notify(PeriodicBoundaryConditionType{ inputrec->ePBC });
+ mdModulesNotifier.notify(inputrec->pbcType);
mdModulesNotifier.notify(SimulationTimeStep{ inputrec->delta_t });
/* Initiate forcerecord */
fr = new t_forcerec;
fr->forceProviders = mdModules_->initForceProviders();
- init_forcerec(fplog, mdlog, fr, fcd, inputrec, &mtop, cr, box,
+ init_forcerec(fplog, mdlog, fr, inputrec.get(), &mtop, cr, box,
opt2fn("-table", filenames.size(), filenames.data()),
opt2fn("-tablep", filenames.size(), filenames.data()),
- opt2fns("-tableb", filenames.size(), filenames.data()), *hwinfo,
- nonbondedDeviceInfo, useGpuForBonded,
- pmeRunMode == PmeRunMode::GPU && !thisRankHasDuty(cr, DUTY_PME), pforce, wcycle);
-
- // TODO Move this to happen during domain decomposition setup,
- // once stream and event handling works well with that.
- // TODO remove need to pass local stream into GPU halo exchange - Redmine #3093
- if (havePPDomainDecomposition(cr) && prefer1DAnd1PulseDD && is1DAnd1PulseDD(*cr->dd))
+ opt2fns("-tableb", filenames.size(), filenames.data()), pforce);
+ // Dirty hack, for fixing disres and orires should be made mdmodules
+ fr->fcdata->disres = disresdata;
+ fr->fcdata->orires = oriresdata;
+
+ // Save a handle to device stream manager to use elsewhere in the code
+ // TODO: Forcerec is not a correct place to store it.
+ fr->deviceStreamManager = deviceStreamManager.get();
+
+ if (runScheduleWork.simulationWork.useGpuPmePpCommunication && !thisRankHasDuty(cr, DUTY_PME))
{
- GMX_RELEASE_ASSERT(devFlags.enableGpuBufferOps,
- "Must use GMX_USE_GPU_BUFFER_OPS=1 to use GMX_GPU_DD_COMMS=1");
- void* streamLocal =
- Nbnxm::gpu_get_command_stream(fr->nbv->gpu_nbv, InteractionLocality::Local);
- void* streamNonLocal =
- Nbnxm::gpu_get_command_stream(fr->nbv->gpu_nbv, InteractionLocality::NonLocal);
- GMX_LOG(mdlog.warning)
- .asParagraph()
- .appendTextFormatted(
- "NOTE: This run uses the 'GPU halo exchange' feature, enabled by the "
- "GMX_GPU_DD_COMMS environment variable.");
- cr->dd->gpuHaloExchange = std::make_unique<GpuHaloExchange>(
- cr->dd, cr->mpi_comm_mysim, streamLocal, streamNonLocal);
+ GMX_RELEASE_ASSERT(
+ deviceStreamManager != nullptr,
+ "GPU device stream manager should be valid in order to use PME-PP direct "
+ "communications.");
+ GMX_RELEASE_ASSERT(
+ deviceStreamManager->streamIsValid(DeviceStreamType::PmePpTransfer),
+ "GPU PP-PME stream should be valid in order to use GPU PME-PP direct "
+ "communications.");
+ fr->pmePpCommGpu = std::make_unique<gmx::PmePpCommGpu>(
+ cr->mpi_comm_mysim, cr->dd->pme_nodeid, deviceStreamManager->context(),
+ deviceStreamManager->stream(DeviceStreamType::PmePpTransfer));
+ }
+
+ fr->nbv = Nbnxm::init_nb_verlet(mdlog, inputrec.get(), fr, cr, *hwinfo_,
+ runScheduleWork.simulationWork.useGpuNonbonded,
+ deviceStreamManager.get(), &mtop, box, wcycle);
+ // TODO: Move the logic below to a GPU bonded builder
+ if (runScheduleWork.simulationWork.useGpuBonded)
+ {
+ GMX_RELEASE_ASSERT(deviceStreamManager != nullptr,
+ "GPU device stream manager should be valid in order to use GPU "
+ "version of bonded forces.");
+ gpuBonded = std::make_unique<GpuBonded>(
+ mtop.ffparams, fr->ic->epsfac * fr->fudgeQQ, deviceStreamManager->context(),
+ deviceStreamManager->bondedStream(havePPDomainDecomposition(cr)), wcycle);
+ fr->gpuBonded = gpuBonded.get();
}
/* Initialize the mdAtoms structure.
}
/* Initialize the virtual site communication */
- vsite = initVsite(mtop, cr);
+ vsite = makeVirtualSitesHandler(mtop, cr, fr->pbcType);
calc_shifts(box, fr->shift_vec);
/* With periodic molecules the charge groups should be whole at start up
* and the virtual sites should not be far from their proper positions.
*/
- if (!inputrec->bContinuation && MASTER(cr) && !(inputrec->ePBC != epbcNONE && inputrec->bPeriodicMols))
+ if (!inputrec->bContinuation && MASTER(cr)
+ && !(inputrec->pbcType != PbcType::No && inputrec->bPeriodicMols))
{
/* Make molecules whole at start of run */
- if (fr->ePBC != epbcNONE)
+ if (fr->pbcType != PbcType::No)
{
- do_pbc_first_mtop(fplog, inputrec->ePBC, box, &mtop, globalState->x.rvec_array());
+ do_pbc_first_mtop(fplog, inputrec->pbcType, box, &mtop, globalState->x.rvec_array());
}
if (vsite)
{
* for the initial distribution in the domain decomposition
* and for the initial shell prediction.
*/
- constructVsitesGlobal(mtop, globalState->x);
+ constructVirtualSitesGlobal(mtop, globalState->x);
}
}
PmeGpuProgramStorage pmeGpuProgram;
if (thisRankHasPmeGpuTask)
{
- pmeGpuProgram = buildPmeGpuProgram(pmeDeviceInfo);
+ GMX_RELEASE_ASSERT(
+ (deviceStreamManager != nullptr),
+ "GPU device stream manager should be initialized in order to use GPU for PME.");
+ GMX_RELEASE_ASSERT((deviceInfo != nullptr),
+ "GPU device should be initialized in order to use GPU for PME.");
+ pmeGpuProgram = buildPmeGpuProgram(deviceStreamManager->context());
}
/* Initiate PME if necessary,
if (cr->npmenodes > 0)
{
/* The PME only nodes need to know nChargePerturbed(FEP on Q) and nTypePerturbed(FEP on LJ)*/
- gmx_bcast_sim(sizeof(nChargePerturbed), &nChargePerturbed, cr);
- gmx_bcast_sim(sizeof(nTypePerturbed), &nTypePerturbed, cr);
+ gmx_bcast(sizeof(nChargePerturbed), &nChargePerturbed, cr->mpi_comm_mysim);
+ gmx_bcast(sizeof(nTypePerturbed), &nTypePerturbed, cr->mpi_comm_mysim);
}
if (thisRankHasDuty(cr, DUTY_PME))
{
try
{
- pmedata = gmx_pme_init(cr, getNumPmeDomains(cr->dd), inputrec, nChargePerturbed != 0,
- nTypePerturbed != 0, mdrunOptions.reproducible, ewaldcoeff_q,
- ewaldcoeff_lj, gmx_omp_nthreads_get(emntPME), pmeRunMode,
- nullptr, pmeDeviceInfo, pmeGpuProgram.get(), mdlog);
+ // TODO: This should be in the builder.
+ GMX_RELEASE_ASSERT(!runScheduleWork.simulationWork.useGpuPme
+ || (deviceStreamManager != nullptr),
+ "Device stream manager should be valid in order to use GPU "
+ "version of PME.");
+ GMX_RELEASE_ASSERT(
+ !runScheduleWork.simulationWork.useGpuPme
+ || deviceStreamManager->streamIsValid(DeviceStreamType::Pme),
+ "GPU PME stream should be valid in order to use GPU version of PME.");
+
+ const DeviceContext* deviceContext = runScheduleWork.simulationWork.useGpuPme
+ ? &deviceStreamManager->context()
+ : nullptr;
+ const DeviceStream* pmeStream =
+ runScheduleWork.simulationWork.useGpuPme
+ ? &deviceStreamManager->stream(DeviceStreamType::Pme)
+ : nullptr;
+
+ pmedata = gmx_pme_init(cr, getNumPmeDomains(cr->dd), inputrec.get(),
+ nChargePerturbed != 0, nTypePerturbed != 0,
+ mdrunOptions.reproducible, ewaldcoeff_q, ewaldcoeff_lj,
+ gmx_omp_nthreads_get(emntPME), pmeRunMode, nullptr,
+ deviceContext, pmeStream, pmeGpuProgram.get(), mdlog);
}
GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
}
if (inputrec->bPull)
{
/* Initialize pull code */
- pull_work = init_pull(fplog, inputrec->pull, inputrec, &mtop, cr, &atomSets,
+ pull_work = init_pull(fplog, inputrec->pull.get(), inputrec.get(), &mtop, cr, &atomSets,
inputrec->fepvals->init_lambda);
if (inputrec->pull->bXOutAverage || inputrec->pull->bFOutAverage)
{
if (inputrec->bRot)
{
/* Initialize enforced rotation code */
- enforcedRotation =
- init_rot(fplog, inputrec, filenames.size(), filenames.data(), cr, &atomSets,
- globalState.get(), &mtop, oenv, mdrunOptions, startingBehavior);
+ enforcedRotation = init_rot(fplog, inputrec.get(), filenames.size(), filenames.data(),
+ cr, &atomSets, globalState.get(), &mtop, oenv, mdrunOptions,
+ startingBehavior);
}
t_swap* swap = nullptr;
if (inputrec->eSwapCoords != eswapNO)
{
/* Initialize ion swapping code */
- swap = init_swapcoords(fplog, inputrec,
+ swap = init_swapcoords(fplog, inputrec.get(),
opt2fn_master("-swap", filenames.size(), filenames.data(), cr),
&mtop, globalState.get(), &observablesHistory, cr, &atomSets,
oenv, mdrunOptions, startingBehavior);
}
/* Let makeConstraints know whether we have essential dynamics constraints. */
- auto constr = makeConstraints(mtop, *inputrec, pull_work, doEssentialDynamics, fplog,
- *mdAtoms->mdatoms(), cr, ms, &nrnb, wcycle, fr->bMolPBC);
+ auto constr = makeConstraints(mtop, *inputrec, pull_work, doEssentialDynamics, fplog, cr,
+ ms, &nrnb, wcycle, fr->bMolPBC);
/* Energy terms and groups */
gmx_enerdata_t enerd(mtop.groups.groups[SimulationAtomGroupType::EnergyOutput].size(),
inputrec->fepvals->n_lambda);
+ // cos acceleration is only supported by md, but older tpr
+ // files might still combine it with other integrators
+ GMX_RELEASE_ASSERT(inputrec->cos_accel == 0.0 || inputrec->eI == eiMD,
+ "cos_acceleration is only supported by integrator=md");
+
/* Kinetic energy data */
gmx_ekindata_t ekind;
init_ekindata(fplog, &mtop, &(inputrec->opts), &ekind, inputrec->cos_accel);
/* Set up interactive MD (IMD) */
auto imdSession =
- makeImdSession(inputrec, cr, wcycle, &enerd, ms, &mtop, mdlog,
+ makeImdSession(inputrec.get(), cr, wcycle, &enerd, ms, &mtop, mdlog,
MASTER(cr) ? globalState->x.rvec_array() : nullptr, filenames.size(),
filenames.data(), oenv, mdrunOptions.imdOptions, startingBehavior);
/* This call is not included in init_domain_decomposition mainly
* because fr->cginfo_mb is set later.
*/
- dd_init_bondeds(fplog, cr->dd, &mtop, vsite.get(), inputrec,
+ dd_init_bondeds(fplog, cr->dd, mtop, vsite.get(), inputrec.get(),
domdecOptions.checkBondedInteractions, fr->cginfo_mb);
}
- // TODO This is not the right place to manage the lifetime of
- // this data structure, but currently it's the easiest way to
- // make it work.
- MdrunScheduleWorkload runScheduleWork;
- // Also populates the simulation constant workload description.
- runScheduleWork.simulationWork = createSimulationWorkload(
- useGpuForNonbonded, pmeRunMode, useGpuForBonded, useGpuForUpdate,
- devFlags.enableGpuBufferOps, devFlags.enableGpuHaloExchange,
- devFlags.enableGpuPmePPComm, haveEwaldSurfaceContribution(*inputrec));
+ if (runScheduleWork.simulationWork.useGpuBufferOps)
+ {
+ fr->gpuForceReduction[gmx::AtomLocality::Local] = std::make_unique<gmx::GpuForceReduction>(
+ deviceStreamManager->context(),
+ deviceStreamManager->stream(gmx::DeviceStreamType::NonBondedLocal), wcycle);
+ fr->gpuForceReduction[gmx::AtomLocality::NonLocal] = std::make_unique<gmx::GpuForceReduction>(
+ deviceStreamManager->context(),
+ deviceStreamManager->stream(gmx::DeviceStreamType::NonBondedNonLocal), wcycle);
+ }
std::unique_ptr<gmx::StatePropagatorDataGpu> stateGpu;
if (gpusWereDetected
- && ((useGpuForPme && thisRankHasDuty(cr, DUTY_PME))
+ && ((runScheduleWork.simulationWork.useGpuPme && thisRankHasDuty(cr, DUTY_PME))
|| runScheduleWork.simulationWork.useGpuBufferOps))
{
- const void* pmeStream = pme_gpu_get_device_stream(fr->pmedata);
- const void* localStream =
- fr->nbv->gpu_nbv != nullptr
- ? Nbnxm::gpu_get_command_stream(fr->nbv->gpu_nbv, InteractionLocality::Local)
- : nullptr;
- const void* nonLocalStream =
- fr->nbv->gpu_nbv != nullptr
- ? Nbnxm::gpu_get_command_stream(fr->nbv->gpu_nbv, InteractionLocality::NonLocal)
- : nullptr;
- const void* deviceContext = pme_gpu_get_device_context(fr->pmedata);
- const int paddingSize = pme_gpu_get_padding_size(fr->pmedata);
GpuApiCallBehavior transferKind = (inputrec->eI == eiMD && !doRerun && !useModularSimulator)
? GpuApiCallBehavior::Async
: GpuApiCallBehavior::Sync;
-
+ GMX_RELEASE_ASSERT(deviceStreamManager != nullptr,
+ "GPU device stream manager should be initialized to use GPU.");
stateGpu = std::make_unique<gmx::StatePropagatorDataGpu>(
- pmeStream, localStream, nonLocalStream, deviceContext, transferKind, paddingSize, wcycle);
+ *deviceStreamManager, transferKind, pme_gpu_get_block_size(fr->pmedata), wcycle);
fr->stateGpu = stateGpu.get();
}
GMX_ASSERT(stopHandlerBuilder_, "Runner must provide StopHandlerBuilder to simulator.");
SimulatorBuilder simulatorBuilder;
+ simulatorBuilder.add(SimulatorStateData(globalState.get(), &observablesHistory, &enerd, &ekind));
+ simulatorBuilder.add(std::move(membedHolder));
+ simulatorBuilder.add(std::move(stopHandlerBuilder_));
+ simulatorBuilder.add(SimulatorConfig(mdrunOptions, startingBehavior, &runScheduleWork));
+
+
+ simulatorBuilder.add(SimulatorEnv(fplog, cr, ms, mdlog, oenv));
+ simulatorBuilder.add(Profiling(&nrnb, walltime_accounting, wcycle));
+ simulatorBuilder.add(ConstraintsParam(
+ constr.get(), enforcedRotation ? enforcedRotation->getLegacyEnfrot() : nullptr,
+ vsite.get()));
+ // TODO: Separate `fr` to a separate add, and make the `build` handle the coupling sensibly.
+ simulatorBuilder.add(LegacyInput(static_cast<int>(filenames.size()), filenames.data(),
+ inputrec.get(), fr));
+ simulatorBuilder.add(ReplicaExchangeParameters(replExParams));
+ simulatorBuilder.add(InteractiveMD(imdSession.get()));
+ simulatorBuilder.add(SimulatorModules(mdModules_->outputProvider(), mdModules_->notifier()));
+ simulatorBuilder.add(CenterOfMassPulling(pull_work));
+ // Todo move to an MDModule
+ simulatorBuilder.add(IonSwapping(swap));
+ simulatorBuilder.add(TopologyData(&mtop, mdAtoms.get()));
+ simulatorBuilder.add(BoxDeformationHandle(deform.get()));
+ simulatorBuilder.add(std::move(modularSimulatorCheckpointData));
+
// build and run simulator object based on user-input
- auto simulator = simulatorBuilder.build(
- inputIsCompatibleWithModularSimulator, fplog, cr, ms, mdlog,
- static_cast<int>(filenames.size()), filenames.data(), oenv, mdrunOptions,
- startingBehavior, vsite.get(), constr.get(),
- enforcedRotation ? enforcedRotation->getLegacyEnfrot() : nullptr, deform.get(),
- mdModules_->outputProvider(), mdModules_->notifier(), inputrec, imdSession.get(),
- pull_work, swap, &mtop, fcd, globalState.get(), &observablesHistory, mdAtoms.get(),
- &nrnb, wcycle, fr, &enerd, &ekind, &runScheduleWork, replExParams, membed,
- walltime_accounting, std::move(stopHandlerBuilder_), doRerun);
+ auto simulator = simulatorBuilder.build(useModularSimulator);
simulator->run();
if (fr->pmePpCommGpu)
GMX_RELEASE_ASSERT(pmedata, "pmedata was NULL while cr->duty was not DUTY_PP");
/* do PME only */
walltime_accounting = walltime_accounting_init(gmx_omp_nthreads_get(emntPME));
- gmx_pmeonly(pmedata, cr, &nrnb, wcycle, walltime_accounting, inputrec, pmeRunMode);
+ gmx_pmeonly(pmedata, cr, &nrnb, wcycle, walltime_accounting, inputrec.get(), pmeRunMode,
+ deviceStreamManager.get());
}
wallcycle_stop(wcycle, ewcRUN);
/* Finish up, write some stuff
* if rerunMD, don't write last frame again
*/
- finish_run(fplog, mdlog, cr, inputrec, &nrnb, wcycle, walltime_accounting,
+ finish_run(fplog, mdlog, cr, inputrec.get(), &nrnb, wcycle, walltime_accounting,
fr ? fr->nbv.get() : nullptr, pmedata, EI_DYNAMICS(inputrec->eI) && !isMultiSim(ms));
// clean up cycle counter
wallcycle_destroy(wcycle);
+ deviceStreamManager.reset(nullptr);
// Free PME data
if (pmedata)
{
}
// FIXME: this is only here to manually unpin mdAtoms->chargeA_ and state->x,
- // before we destroy the GPU context(s) in free_gpu_resources().
+ // before we destroy the GPU context(s)
// Pinned buffers are associated with contexts in CUDA.
// As soon as we destroy GPU contexts after mdrunner() exits, these lines should go.
mdAtoms.reset(nullptr);
globalState.reset(nullptr);
mdModules_.reset(nullptr); // destruct force providers here as they might also use the GPU
+ gpuBonded.reset(nullptr);
+ /* Free pinned buffers in *fr */
+ delete fr;
+ fr = nullptr;
+ // TODO convert to C++ so we can get rid of these frees
+ sfree(disresdata);
+ sfree(oriresdata);
- /* Free GPU memory and set a physical node tMPI barrier (which should eventually go away) */
- free_gpu_resources(fr, physicalNodeComm, hwinfo->gpu_info);
- free_gpu(nonbondedDeviceInfo);
- free_gpu(pmeDeviceInfo);
- done_forcerec(fr, mtop.molblock.size());
- sfree(fcd);
+ if (!hwinfo_->deviceInfoList.empty())
+ {
+ /* stop the GPU profiler (only CUDA) */
+ stopGpuProfiler();
+ }
- if (doMembed)
+ /* With tMPI we need to wait for all ranks to finish deallocation before
+ * destroying the CUDA context as some tMPI ranks may be sharing
+ * GPU and context.
+ *
+ * This is not a concern in OpenCL where we use one context per rank.
+ *
+ * Note: it is safe to not call the barrier on the ranks which do not use GPU,
+ * but it is easier and more futureproof to call it on the whole node.
+ *
+ * Note that this function needs to be called even if GPUs are not used
+ * in this run because the PME ranks have no knowledge of whether GPUs
+ * are used or not, but all ranks need to enter the barrier below.
+ * \todo Remove this physical node barrier after making sure
+ * that it's not needed anymore (with a shared GPU run).
+ */
+ if (GMX_THREAD_MPI)
{
- free_membed(membed);
+ physicalNodeComm.barrier();
}
+ releaseDevice(deviceInfo);
/* Does what it says */
print_date_and_time(fplog, cr->nodeid, "Finished mdrun", gmx_gettime());
/* we need to join all threads. The sub-threads join when they
exit this function, but the master thread needs to be told to
wait for that. */
- if (PAR(cr) && MASTER(cr))
+ if (MASTER(cr))
{
tMPI_Finalize();
}
#endif
return rc;
-}
+} // namespace gmx
Mdrunner::~Mdrunner()
{
Mdrunner::Mdrunner(Mdrunner&&) noexcept = default;
-//NOLINTNEXTLINE(performance-noexcept-move-constructor) working around GCC bug 58265
+//NOLINTNEXTLINE(performance-noexcept-move-constructor) working around GCC bug 58265 in CentOS 7
Mdrunner& Mdrunner::operator=(Mdrunner&& /*handle*/) noexcept(BUGFREE_NOEXCEPT_STRING) = default;
class Mdrunner::BuilderImplementation
real forceWarningThreshold,
StartingBehavior startingBehavior);
+ void addHardwareDetectionResult(const gmx_hw_info_t* hwinfo);
+
void addDomdec(const DomdecOptions& options);
+ void addInput(SimulationInputHandle inputHolder);
+
void addVerletList(int nstlist);
void addReplicaExchange(const ReplicaExchangeParameters& params);
//! Command-line override for the duration of a neighbor list with the Verlet scheme.
int nstlist_ = 0;
+ //! World communicator, used for hardware detection and task assignment
+ MPI_Comm libraryWorldCommunicator_ = MPI_COMM_NULL;
+
//! Multisim communicator handle.
gmx_multisim_t* multiSimulation_;
//! mdrun communicator
- MPI_Comm communicator_ = MPI_COMM_NULL;
+ MPI_Comm simulationCommunicator_ = MPI_COMM_NULL;
//! Print a warning if any force is larger than this (in kJ/mol nm).
real forceWarningThreshold_ = -1;
//! The modules that comprise the functionality of mdrun.
std::unique_ptr<MDModules> mdModules_;
+ //! Detected hardware.
+ const gmx_hw_info_t* hwinfo_ = nullptr;
+
//! \brief Parallelism information.
gmx_hw_opt_t hardwareOptions_;
* \brief Builder for simulation stop signal handler.
*/
std::unique_ptr<StopHandlerBuilder> stopHandlerBuilder_ = nullptr;
+
+ /*!
+ * \brief Sources for initial simulation state.
+ *
+ * See issue #3652 for near-term refinements to the SimulationInput interface.
+ *
+ * See issue #3379 for broader discussion on API aspects of simulation inputs and outputs.
+ */
+ SimulationInputHandle inputHolder_;
};
Mdrunner::BuilderImplementation::BuilderImplementation(std::unique_ptr<MDModules> mdModules,
compat::not_null<SimulationContext*> context) :
mdModules_(std::move(mdModules))
{
- communicator_ = context->communicator_;
- multiSimulation_ = context->multiSimulation_.get();
+ libraryWorldCommunicator_ = context->libraryWorldCommunicator_;
+ simulationCommunicator_ = context->simulationCommunicator_;
+ multiSimulation_ = context->multiSimulation_.get();
}
Mdrunner::BuilderImplementation::~BuilderImplementation() = default;
newRunner.filenames = filenames_;
- newRunner.communicator = communicator_;
+ newRunner.libraryWorldCommunicator = libraryWorldCommunicator_;
+
+ newRunner.simulationCommunicator = simulationCommunicator_;
// nullptr is a valid value for the multisim handle
newRunner.ms = multiSimulation_;
+ if (hwinfo_)
+ {
+ newRunner.hwinfo_ = hwinfo_;
+ }
+ else
+ {
+ GMX_THROW(gmx::APIError(
+ "MdrunnerBuilder::addHardwareDetectionResult() is required before build()"));
+ }
+
+ if (inputHolder_)
+ {
+ newRunner.inputHolder_ = std::move(inputHolder_);
+ }
+ else
+ {
+ GMX_THROW(gmx::APIError("MdrunnerBuilder::addInput() is required before build()."));
+ }
+
// \todo Clarify ownership and lifetime management for gmx_output_env_t
// \todo Update sanity checking when output environment has clearly specified invariants.
// Initialization and default values for oenv are not well specified in the current version.
return newRunner;
}
+void Mdrunner::BuilderImplementation::addHardwareDetectionResult(const gmx_hw_info_t* hwinfo)
+{
+ hwinfo_ = hwinfo;
+}
+
void Mdrunner::BuilderImplementation::addNonBonded(const char* nbpu_opt)
{
nbpu_opt_ = nbpu_opt;
stopHandlerBuilder_ = std::move(builder);
}
+void Mdrunner::BuilderImplementation::addInput(SimulationInputHandle inputHolder)
+{
+ inputHolder_ = std::move(inputHolder);
+}
+
MdrunnerBuilder::MdrunnerBuilder(std::unique_ptr<MDModules> mdModules,
compat::not_null<SimulationContext*> context) :
impl_{ std::make_unique<Mdrunner::BuilderImplementation>(std::move(mdModules), context) }
MdrunnerBuilder::~MdrunnerBuilder() = default;
+MdrunnerBuilder& MdrunnerBuilder::addHardwareDetectionResult(const gmx_hw_info_t* hwinfo)
+{
+ impl_->addHardwareDetectionResult(hwinfo);
+ return *this;
+}
+
MdrunnerBuilder& MdrunnerBuilder::addSimulationMethod(const MdrunOptions& options,
real forceWarningThreshold,
const StartingBehavior startingBehavior)
return *this;
}
+MdrunnerBuilder& MdrunnerBuilder::addInput(SimulationInputHandle input)
+{
+ impl_->addInput(std::move(input));
+ return *this;
+}
+
MdrunnerBuilder::MdrunnerBuilder(MdrunnerBuilder&&) noexcept = default;
MdrunnerBuilder& MdrunnerBuilder::operator=(MdrunnerBuilder&&) noexcept = default;
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff