#include "grompp.h"
+#include <array>
#include <cerrno>
#include <climits>
#include <cmath>
#include "gromacs/utility/listoflists.h"
#include "gromacs/utility/logger.h"
#include "gromacs/utility/loggerbuilder.h"
-#include "gromacs/utility/mdmodulenotification.h"
+#include "gromacs/utility/mdmodulesnotifiers.h"
#include "gromacs/utility/smalloc.h"
#include "gromacs/utility/snprintf.h"
InteractionOfType::InteractionOfType(gmx::ArrayRef<const int> atoms,
gmx::ArrayRef<const real> params,
const std::string& name) :
- atoms_(atoms.begin(), atoms.end()),
- interactionTypeName_(name)
+ atoms_(atoms.begin(), atoms.end()), interactionTypeName_(name)
{
GMX_RELEASE_ASSERT(
params.size() <= forceParam_.size(),
gmx::formatString("Cannot have more parameters than the maximum number possible (%d)", MAXFORCEPARAM)
.c_str());
- auto forceParamIt = forceParam_.begin();
+ std::array<real, MAXFORCEPARAM>::iterator forceParamIt = forceParam_.begin();
for (const auto param : params)
{
*forceParamIt++ = param;
{
const t_atom& local = atomP.atom();
int i = atomP.globalAtomNumber();
- if (local.ptype == eptShell || local.ptype == eptBond || local.ptype == eptVSite)
+ if (local.ptype == ParticleType::Shell || local.ptype == ParticleType::Bond
+ || local.ptype == ParticleType::VSite)
{
clear_rvec(v[i]);
}
for (const AtomProxy atomP : AtomRange(*mtop))
{
const t_atom& local = atomP.atom();
- if (local.ptype == eptShell || local.ptype == eptBond)
+ if (local.ptype == ParticleType::Shell || local.ptype == ParticleType::Bond)
{
nshells++;
}
* constraints only. Do not print note with large timesteps or vsites.
*/
if (opts->nshake == eshALLBONDS && ffParametrizedWithHBondConstraints && ir->delta_t < 0.0026
- && gmx_mtop_ftype_count(sys, F_VSITE3FD) == 0)
+ && gmx_mtop_ftype_count(*sys, F_VSITE3FD) == 0)
{
set_warning_line(wi, "unknown", -1);
warning_note(wi,
/* It would be nice to get rid of the copies below, but we don't know
* a priori if the number of atoms in confin matches what we expect.
*/
- state->flags |= (1 << estX);
+ state->flags |= enumValueToBitMask(StateEntry::X);
if (v != nullptr)
{
- state->flags |= (1 << estV);
+ state->flags |= enumValueToBitMask(StateEntry::V);
}
state_change_natoms(state, state->natoms);
std::copy(x, x + state->natoms, state->x.data());
opts->seed = static_cast<int>(gmx::makeRandomSeed());
GMX_LOG(logger.info).asParagraph().appendTextFormatted("Setting gen_seed to %d", opts->seed);
}
- state->flags |= (1 << estV);
+ state->flags |= enumValueToBitMask(StateEntry::V);
maxwell_speed(opts->tempi, opts->seed, sys, state->v.rvec_array(), logger);
stop_cm(logger, state->natoms, mass, state->x.rvec_array(), state->v.rvec_array());
}
for (i = 0; i < ir->nwall; i++)
{
- ir->wall_atomtype[i] = at->atomTypeFromName(opts->wall_atomtype[i]);
- if (ir->wall_atomtype[i] == NOTSET)
+ auto atomType = at->atomTypeFromName(opts->wall_atomtype[i]);
+ if (!atomType.has_value())
{
std::string warningMessage = gmx::formatString(
"Specified wall atom type %s is not defined", opts->wall_atomtype[i]);
warning_error(wi, warningMessage.c_str());
}
+ else
+ {
+ ir->wall_atomtype[i] = *atomType;
+ }
}
}
for (i = 0; i < atoms->nr; i++)
{
/* Vsite ptype might not be set here yet, so also check the mass */
- if ((atoms->atom[i].ptype == eptAtom || atoms->atom[i].ptype == eptNucleus) && atoms->atom[i].m > 0)
+ if ((atoms->atom[i].ptype == ParticleType::Atom || atoms->atom[i].ptype == ParticleType::Nucleus)
+ && atoms->atom[i].m > 0)
{
nmass++;
}
nrdf += ir->opts.nrdf[g];
}
- return sum_mv2 / (nrdf * BOLTZ);
+ return sum_mv2 / (nrdf * gmx::c_boltz);
}
static real get_max_reference_temp(const t_inputrec* ir, warninp* wi)
int numDanglingAtoms = 0;
for (int a = 0; a < atoms->nr; a++)
{
- if (atoms->atom[a].ptype != eptVSite && count[a] == 0)
+ if (atoms->atom[a].ptype != ParticleType::VSite && count[a] == 0)
{
if (bVerbose)
{
/*! \brief Checks if the Verlet buffer and constraint accuracy is sufficient for decoupled dynamic modes.
*
- * When decoupled modes are present and the accuray in insufficient,
+ * When decoupled modes are present and the accuracy in insufficient,
* this routine issues a warning if the accuracy is insufficient.
*/
static void checkDecoupledModeAccuracy(const gmx_mtop_t* mtop, const t_inputrec* ir, warninp* wi)
{ efEDR, "-e", nullptr, ffOPTRD },
/* This group is needed by the VMD viewer as the start configuration for IMD sessions: */
{ efGRO, "-imd", "imdgroup", ffOPTWR },
- { efTRN, "-ref", "rotref", ffOPTRW | ffALLOW_MISSING } };
+ { efTRN, "-ref", "rotref", ffOPTRW | ffALLOW_MISSING },
+ /* This group is needed by the QMMM MDModule: */
+ { efQMI, "-qmi", nullptr, ffOPTRD } };
#define NFILE asize(fnm)
/* Command line options */
}
GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
- // Now that the MdModules have their options assigned from get_ir, subscribe
+ // Now that the MDModules have their options assigned from get_ir, subscribe
// to eventual notifications during pre-processing their data
mdModules.subscribeToPreProcessingNotifications();
+ // Notify MDModules of existing logger
+ mdModules.notifiers().preProcessingNotifier_.notify(logger);
+
+ // Notify MDModules of existing warninp
+ mdModules.notifiers().preProcessingNotifier_.notify(wi);
+
+ // Notify QMMM MDModule of external QM input file command-line option
+ {
+ gmx::QMInputFileName qmInputFileName = { ftp2bSet(efQMI, NFILE, fnm), ftp2fn(efQMI, NFILE, fnm) };
+ mdModules.notifiers().preProcessingNotifier_.notify(qmInputFileName);
+ }
if (bVerbose)
{
.asParagraph()
.appendTextFormatted("checking input for internal consistency...");
}
- check_ir(mdparin, mdModules.notifier(), ir, opts, wi);
+ check_ir(mdparin, mdModules.notifiers(), ir, opts, wi);
if (ir->ld_seed == -1)
{
/* check masses */
check_mol(&sys, wi);
+ if (haveFepPerturbedMassesInSettles(sys))
+ {
+ warning_error(wi,
+ "SETTLE is not implemented for atoms whose mass is perturbed. "
+ "You might instead use normal constraints.");
+ }
+
checkForUnboundAtoms(&sys, bVerbose, wi, logger);
if (EI_DYNAMICS(ir->eI) && ir->eI != IntegrationAlgorithm::BD)
{
GMX_LOG(logger.info).asParagraph().appendTextFormatted("initialising group options...");
}
- do_index(mdparin, ftp2fn_null(efNDX, NFILE, fnm), &sys, bVerbose, mdModules.notifier(), ir, wi);
+ do_index(mdparin, ftp2fn_null(efNDX, NFILE, fnm), &sys, bVerbose, mdModules.notifiers(), ir, wi);
+
+ // Notify topology to MdModules for pre-processing after all indexes were built
+ mdModules.notifiers().preProcessingNotifier_.notify(&sys);
if (ir->cutoff_scheme == CutoffScheme::Verlet && ir->verletbuf_tol > 0)
{
* to be on the safe side with constraints.
*/
const real totalEnergyDriftPerAtomPerPicosecond =
- 2 * BOLTZ * buffer_temp / (ir->nsteps * ir->delta_t);
+ 2 * gmx::c_boltz * buffer_temp / (ir->nsteps * ir->delta_t);
if (ir->verletbuf_tol > 1.1 * driftTolerance * totalEnergyDriftPerAtomPerPicosecond)
{
{
/* Calculate the optimal grid dimensions */
matrix scaledBox;
- EwaldBoxZScaler boxScaler(*ir);
+ EwaldBoxZScaler boxScaler(inputrecPbcXY2Walls(ir), ir->wall_ewald_zfac);
boxScaler.scaleBox(state.box, scaledBox);
if (ir->nkx > 0 && ir->nky > 0 && ir->nkz > 0)
if (ir->bPull)
{
- pull = set_pull_init(ir,
- &sys,
- state.x.rvec_array(),
- state.box,
- state.lambda[FreeEnergyPerturbationCouplingType::Mass],
- wi);
+ pull = set_pull_init(
+ ir, sys, state.x, state.box, state.lambda[FreeEnergyPerturbationCouplingType::Mass], wi);
}
/* Modules that supply external potential for pull coordinates
copy_mat(ir->compress, compressibility);
}
setStateDependentAwhParams(
- ir->awhParams,
+ ir->awhParams.get(),
*ir->pull,
pull,
state.box,
ir->pbcType,
compressibility,
&ir->opts,
- ir->efep != FreeEnergyPerturbationType::No
- ? ir->fepvals->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Fep)]
- [ir->fepvals->init_fep_state]
- : 0,
+ ir->efep != FreeEnergyPerturbationType::No ? ir->fepvals->all_lambda[static_cast<int>(
+ FreeEnergyPerturbationCouplingType::Fep)][ir->fepvals->init_fep_state]
+ : 0,
sys,
wi);
}
}
}
+ // Hand over box and coordiantes to MdModules before they evaluate their final parameters
+ {
+ gmx::CoordinatesAndBoxPreprocessed coordinatesAndBoxPreprocessed;
+ coordinatesAndBoxPreprocessed.coordinates_ = state.x.arrayRefWithPadding();
+ copy_mat(state.box, coordinatesAndBoxPreprocessed.box_);
+ coordinatesAndBoxPreprocessed.pbc_ = ir->pbcType;
+ mdModules.notifiers().preProcessingNotifier_.notify(coordinatesAndBoxPreprocessed);
+ }
+
// Add the md modules internal parameters that are not mdp options
// e.g., atom indices
{
gmx::KeyValueTreeBuilder internalParameterBuilder;
- mdModules.notifier().preProcessingNotifications_.notify(internalParameterBuilder.rootObject());
+ mdModules.notifiers().preProcessingNotifier_.notify(internalParameterBuilder.rootObject());
ir->internalParameters =
std::make_unique<gmx::KeyValueTreeObject>(internalParameterBuilder.build());
}
}
done_warning(wi, FARGS);
- write_tpx_state(ftp2fn(efTPR, NFILE, fnm), ir, &state, &sys);
+ write_tpx_state(ftp2fn(efTPR, NFILE, fnm), ir, &state, sys);
/* Output IMD group, if bIMD is TRUE */
- gmx::write_IMDgroup_to_file(ir->bIMD, ir, &state, &sys, NFILE, fnm);
+ gmx::write_IMDgroup_to_file(ir->bIMD, ir, &state, sys, NFILE, fnm);
sfree(opts->define);
sfree(opts->wall_atomtype[0]);
sfree(opts->wall_atomtype[1]);
sfree(opts->include);
+ sfree(opts->couple_moltype);
+
for (auto& mol : mi)
{
// Some of the contents of molinfo have been stolen, so