*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2016,2017,2018,2019, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2016,2017 by 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 "grompp.h"
+#include <array>
#include <cerrno>
#include <climits>
#include <cmath>
#include <sys/types.h>
-#include "gromacs/awh/read_params.h"
+#include "gromacs/applied_forces/awh/read_params.h"
#include "gromacs/commandline/pargs.h"
#include "gromacs/ewald/ewald_utils.h"
#include "gromacs/ewald/pme.h"
#include "gromacs/mdlib/calc_verletbuf.h"
#include "gromacs/mdlib/compute_io.h"
#include "gromacs/mdlib/constr.h"
+#include "gromacs/mdlib/md_support.h"
#include "gromacs/mdlib/perf_est.h"
-#include "gromacs/mdlib/qmmm.h"
#include "gromacs/mdlib/vsite.h"
#include "gromacs/mdrun/mdmodules.h"
#include "gromacs/mdtypes/inputrec.h"
#include "gromacs/utility/cstringutil.h"
#include "gromacs/utility/exceptions.h"
#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/filestream.h"
#include "gromacs/utility/futil.h"
#include "gromacs/utility/gmxassert.h"
#include "gromacs/utility/keyvaluetreebuilder.h"
-#include "gromacs/utility/mdmodulenotification.h"
+#include "gromacs/utility/listoflists.h"
+#include "gromacs/utility/logger.h"
+#include "gromacs/utility/loggerbuilder.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;
void MoleculeInformation::initMolInfo()
{
init_block(&mols);
- init_blocka(&excls);
+ excls.clear();
init_t_atoms(&atoms, 0, FALSE);
}
return n;
}
-static int check_atom_names(const char* fn1, const char* fn2, gmx_mtop_t* mtop, const t_atoms* at)
+static int check_atom_names(const char* fn1,
+ const char* fn2,
+ gmx_mtop_t* mtop,
+ const t_atoms* at,
+ const gmx::MDLogger& logger)
{
int m, i, j, nmismatch;
t_atoms* tat;
-#define MAXMISMATCH 20
+
+ constexpr int c_maxNumberOfMismatches = 20;
if (mtop->natoms != at->nr)
{
{
if (strcmp(*(tat->atomname[j]), *(at->atomname[i])) != 0)
{
- if (nmismatch < MAXMISMATCH)
+ if (nmismatch < c_maxNumberOfMismatches)
{
- fprintf(stderr,
- "Warning: atom name %d in %s and %s does not match (%s - %s)\n",
- i + 1, fn1, fn2, *(tat->atomname[j]), *(at->atomname[i]));
+ GMX_LOG(logger.warning)
+ .asParagraph()
+ .appendTextFormatted(
+ "atom name %d in %s and %s does not match (%s - %s)",
+ i + 1,
+ fn1,
+ fn2,
+ *(tat->atomname[j]),
+ *(at->atomname[i]));
}
- else if (nmismatch == MAXMISMATCH)
+ else if (nmismatch == c_maxNumberOfMismatches)
{
- fprintf(stderr, "(more than %d non-matching atom names)\n", MAXMISMATCH);
+ GMX_LOG(logger.warning)
+ .asParagraph()
+ .appendTextFormatted("(more than %d non-matching atom names)",
+ c_maxNumberOfMismatches);
}
nmismatch++;
}
int i, a1, a2, w_a1, w_a2, j;
real twopi2, limit2, fc, re, m1, m2, period2, w_period2;
bool bFound, bWater, bWarn;
- char warn_buf[STRLEN];
/* Get the interaction parameters */
gmx::ArrayRef<const t_iparams> ip = mtop->ffparams.iparams;
bWarn = (w_period2 < gmx::square(min_steps_warn * dt));
/* A check that would recognize most water models */
bWater = ((*w_moltype->atoms.atomname[0])[0] == 'O' && w_moltype->atoms.nr <= 5);
- sprintf(warn_buf,
+ std::string warningMessage = gmx::formatString(
"The bond in molecule-type %s between atoms %d %s and %d %s has an estimated "
"oscillational period of %.1e ps, which is less than %d times the time step of "
"%.1e ps.\n"
"%s",
- *w_moltype->name, w_a1 + 1, *w_moltype->atoms.atomname[w_a1], w_a2 + 1,
- *w_moltype->atoms.atomname[w_a2], std::sqrt(w_period2),
- bWarn ? min_steps_warn : min_steps_note, dt,
+ *w_moltype->name,
+ w_a1 + 1,
+ *w_moltype->atoms.atomname[w_a1],
+ w_a2 + 1,
+ *w_moltype->atoms.atomname[w_a2],
+ std::sqrt(w_period2),
+ bWarn ? min_steps_warn : min_steps_note,
+ dt,
bWater ? "Maybe you asked for fexible water."
: "Maybe you forgot to change the constraints mdp option.");
if (bWarn)
{
- warning(wi, warn_buf);
+ warning(wi, warningMessage.c_str());
}
else
{
- warning_note(wi, warn_buf);
+ warning_note(wi, warningMessage.c_str());
}
}
}
{
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]);
}
static void check_shells_inputrec(gmx_mtop_t* mtop, t_inputrec* ir, warninp* wi)
{
- int nshells = 0;
- char warn_buf[STRLEN];
+ int nshells = 0;
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++;
}
if ((nshells > 0) && (ir->nstcalcenergy != 1))
{
set_warning_line(wi, "unknown", -1);
- snprintf(warn_buf, STRLEN, "There are %d shells, changing nstcalcenergy from %d to 1",
- nshells, ir->nstcalcenergy);
+ std::string warningMessage = gmx::formatString(
+ "There are %d shells, changing nstcalcenergy from %d to 1", nshells, ir->nstcalcenergy);
ir->nstcalcenergy = 1;
- warning(wi, warn_buf);
+ warning(wi, warningMessage.c_str());
}
}
std::vector<MoleculeInformation>* mi,
std::unique_ptr<MoleculeInformation>* intermolecular_interactions,
gmx::ArrayRef<InteractionsOfType> interactions,
- int* comb,
+ CombinationRule* comb,
double* reppow,
real* fudgeQQ,
gmx_bool bMorse,
- warninp* wi)
+ warninp* wi,
+ const gmx::MDLogger& logger)
{
std::vector<gmx_molblock_t> molblock;
int i, nmismatch;
bool ffParametrizedWithHBondConstraints;
- char buf[STRLEN];
- char warn_buf[STRLEN];
/* TOPOLOGY processing */
- sys->name = do_top(bVerbose, topfile, topppfile, opts, bZero, &(sys->symtab), interactions,
- comb, reppow, fudgeQQ, atypes, mi, intermolecular_interactions, ir,
- &molblock, &ffParametrizedWithHBondConstraints, wi);
+ sys->name = do_top(bVerbose,
+ topfile,
+ topppfile,
+ opts,
+ bZero,
+ &(sys->symtab),
+ interactions,
+ comb,
+ reppow,
+ fudgeQQ,
+ atypes,
+ mi,
+ intermolecular_interactions,
+ ir,
+ &molblock,
+ &ffParametrizedWithHBondConstraints,
+ wi,
+ logger);
sys->molblock.clear();
convert_harmonics(*mi, atypes);
}
- if (ir->eDisre == edrNone)
+ if (ir->eDisre == DistanceRestraintRefinement::None)
{
i = rm_interactions(F_DISRES, *mi);
if (i > 0)
{
set_warning_line(wi, "unknown", -1);
- sprintf(warn_buf, "disre = no, removed %d distance restraints", i);
- warning_note(wi, warn_buf);
+ std::string warningMessage =
+ gmx::formatString("disre = no, removed %d distance restraints", i);
+ warning_note(wi, warningMessage.c_str());
}
}
if (!opts->bOrire)
if (i > 0)
{
set_warning_line(wi, "unknown", -1);
- sprintf(warn_buf, "orire = no, removed %d orientation restraints", i);
- warning_note(wi, warn_buf);
+ std::string warningMessage =
+ gmx::formatString("orire = no, removed %d orientation restraints", i);
+ warning_note(wi, warningMessage.c_str());
}
}
/* Copy structures from msys to sys */
molinfo2mtop(*mi, sys);
- gmx_mtop_finalize(sys);
+ sys->finalize();
/* Discourage using the, previously common, setup of applying constraints
* to all bonds with force fields that have been parametrized with H-bond
* 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,
/* COORDINATE file processing */
if (bVerbose)
{
- fprintf(stderr, "processing coordinates...\n");
+ GMX_LOG(logger.info).asParagraph().appendTextFormatted("processing coordinates...");
}
t_topology* conftop;
gmx_fatal(FARGS,
"number of coordinates in coordinate file (%s, %d)\n"
" does not match topology (%s, %d)",
- confin, state->natoms, topfile, sys->natoms);
+ confin,
+ state->natoms,
+ topfile,
+ sys->natoms);
}
/* 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());
/* This call fixes the box shape for runs with pressure scaling */
set_box_rel(ir, state);
- nmismatch = check_atom_names(topfile, confin, sys, &conftop->atoms);
+ nmismatch = check_atom_names(topfile, confin, sys, &conftop->atoms, logger);
done_top(conftop);
sfree(conftop);
if (nmismatch)
{
- sprintf(buf,
+ std::string warningMessage = gmx::formatString(
"%d non-matching atom name%s\n"
"atom names from %s will be used\n"
"atom names from %s will be ignored\n",
- nmismatch, (nmismatch == 1) ? "" : "s", topfile, confin);
- warning(wi, buf);
+ nmismatch,
+ (nmismatch == 1) ? "" : "s",
+ topfile,
+ confin);
+ warning(wi, warningMessage.c_str());
}
/* Do more checks, mostly related to constraints */
if (bVerbose)
{
- fprintf(stderr, "double-checking input for internal consistency...\n");
+ GMX_LOG(logger.info)
+ .asParagraph()
+ .appendTextFormatted("double-checking input for internal consistency...");
}
{
bool bHasNormalConstraints =
if (opts->seed == -1)
{
opts->seed = static_cast<int>(gmx::makeRandomSeed());
- fprintf(stderr, "Setting gen_seed to %d\n", opts->seed);
+ GMX_LOG(logger.info).asParagraph().appendTextFormatted("Setting gen_seed to %d", opts->seed);
}
- state->flags |= (1 << estV);
- maxwell_speed(opts->tempi, opts->seed, sys, state->v.rvec_array());
+ state->flags |= enumValueToBitMask(StateEntry::V);
+ maxwell_speed(opts->tempi, opts->seed, sys, state->v.rvec_array(), logger);
- stop_cm(stdout, state->natoms, mass, state->x.rvec_array(), state->v.rvec_array());
+ stop_cm(logger, state->natoms, mass, state->x.rvec_array(), state->v.rvec_array());
sfree(mass);
}
}
t_inputrec* ir,
t_state* state,
gmx_mtop_t* sys,
- const gmx_output_env_t* oenv)
+ const gmx_output_env_t* oenv,
+ const gmx::MDLogger& logger)
/* If fr_time == -1 read the last frame available which is complete */
{
bool bReadVel;
bReadVel = (bNeedVel && !bGenVel);
- fprintf(stderr, "Reading Coordinates%s and Box size from old trajectory\n",
- bReadVel ? ", Velocities" : "");
+ GMX_LOG(logger.info)
+ .asParagraph()
+ .appendTextFormatted("Reading Coordinates%s and Box size from old trajectory",
+ bReadVel ? ", Velocities" : "");
if (fr_time == -1)
{
- fprintf(stderr, "Will read whole trajectory\n");
+ GMX_LOG(logger.info).asParagraph().appendTextFormatted("Will read whole trajectory");
}
else
{
- fprintf(stderr, "Will read till time %g\n", fr_time);
+ GMX_LOG(logger.info).asParagraph().appendTextFormatted("Will read till time %g", fr_time);
}
if (!bReadVel)
{
if (bGenVel)
{
- fprintf(stderr,
- "Velocities generated: "
- "ignoring velocities in input trajectory\n");
+ GMX_LOG(logger.info)
+ .asParagraph()
+ .appendTextFormatted(
+ "Velocities generated: "
+ "ignoring velocities in input trajectory");
}
read_first_frame(oenv, &fp, slog, &fr, TRX_NEED_X);
}
if (!fr.bV)
{
- fprintf(stderr,
- "\n"
- "WARNING: Did not find a frame with velocities in file %s,\n"
- " all velocities will be set to zero!\n\n",
- slog);
+ GMX_LOG(logger.warning)
+ .asParagraph()
+ .appendTextFormatted(
+ "WARNING: Did not find a frame with velocities in file %s,\n"
+ " all velocities will be set to zero!",
+ slog);
for (auto& vi : makeArrayRef(state->v))
{
vi = { 0, 0, 0 };
*/
set_box_rel(ir, state);
- fprintf(stderr, "Using frame at t = %g ps\n", use_time);
- fprintf(stderr, "Starting time for run is %g ps\n", ir->init_t);
+ GMX_LOG(logger.info).asParagraph().appendTextFormatted("Using frame at t = %g ps", use_time);
+ GMX_LOG(logger.info).asParagraph().appendTextFormatted("Starting time for run is %g ps", ir->init_t);
- if ((ir->epc != epcNO || ir->etc == etcNOSEHOOVER) && ener)
+ if ((ir->epc != PressureCoupling::No || ir->etc == TemperatureCoupling::NoseHoover) && ener)
{
get_enx_state(ener, use_time, sys->groups, ir, state);
preserve_box_shape(ir, state->box_rel, state->boxv);
gmx::ArrayRef<const MoleculeInformation> molinfo,
gmx_bool bTopB,
const char* fn,
- int rc_scaling,
- int ePBC,
+ RefCoordScaling rc_scaling,
+ PbcType pbcType,
rvec com,
- warninp* wi)
+ warninp* wi,
+ const gmx::MDLogger& logger)
{
gmx_bool* hadAtom;
rvec * x, *v;
t_topology* top;
matrix box, invbox;
int natoms, npbcdim = 0;
- char warn_buf[STRLEN];
int a, nat_molb;
t_atom* atom;
sfree(top);
if (natoms != mtop->natoms)
{
- sprintf(warn_buf,
+ std::string warningMessage = gmx::formatString(
"The number of atoms in %s (%d) does not match the number of atoms in the topology "
"(%d). Will assume that the first %d atoms in the topology and %s match.",
- fn, natoms, mtop->natoms, std::min(mtop->natoms, natoms), fn);
- warning(wi, warn_buf);
+ fn,
+ natoms,
+ mtop->natoms,
+ std::min(mtop->natoms, natoms),
+ fn);
+ warning(wi, warningMessage.c_str());
}
- npbcdim = ePBC2npbcdim(ePBC);
+ npbcdim = numPbcDimensions(pbcType);
GMX_RELEASE_ASSERT(npbcdim <= DIM, "Invalid npbcdim");
clear_rvec(com);
- if (rc_scaling != erscNO)
+ if (rc_scaling != RefCoordScaling::No)
{
copy_mat(box, invbox);
for (int j = npbcdim; j < DIM; j++)
gmx_fatal(FARGS,
"Position restraint atom index (%d) in moltype '%s' is larger than "
"number of atoms in %s (%d).\n",
- ai + 1, *molinfo[molb.type].name, fn, natoms);
+ ai + 1,
+ *molinfo[molb.type].name,
+ fn,
+ natoms);
}
hadAtom[ai] = TRUE;
- if (rc_scaling == erscCOM)
+ if (rc_scaling == RefCoordScaling::Com)
{
/* Determine the center of mass of the posres reference coordinates */
for (int j = 0; j < npbcdim; j++)
gmx_fatal(FARGS,
"Position restraint atom index (%d) in moltype '%s' is larger than "
"number of atoms in %s (%d).\n",
- ai + 1, *molinfo[molb.type].name, fn, natoms);
+ ai + 1,
+ *molinfo[molb.type].name,
+ fn,
+ natoms);
}
- if (rc_scaling == erscCOM && !hadAtom[ai])
+ if (rc_scaling == RefCoordScaling::Com && !hadAtom[ai])
{
/* Determine the center of mass of the posres reference coordinates */
for (int j = 0; j < npbcdim; j++)
}
a += nat_molb;
}
- if (rc_scaling == erscCOM)
+ if (rc_scaling == RefCoordScaling::Com)
{
if (totmass == 0)
{
{
com[j] = sum[j] / totmass;
}
- fprintf(stderr,
- "The center of mass of the position restraint coord's is %6.3f %6.3f %6.3f\n",
- com[XX], com[YY], com[ZZ]);
+ GMX_LOG(logger.info)
+ .asParagraph()
+ .appendTextFormatted(
+ "The center of mass of the position restraint coord's is %6.3f %6.3f %6.3f",
+ com[XX],
+ com[YY],
+ com[ZZ]);
}
- if (rc_scaling != erscNO)
+ if (rc_scaling != RefCoordScaling::No)
{
GMX_ASSERT(npbcdim <= DIM, "Only DIM dimensions can have PBC");
{
for (int j = 0; j < npbcdim; j++)
{
- if (rc_scaling == erscALL)
+ if (rc_scaling == RefCoordScaling::All)
{
/* Convert from Cartesian to crystal coordinates */
xp[i][j] *= invbox[j][j];
xp[i][j] += invbox[k][j] * xp[i][k];
}
}
- else if (rc_scaling == erscCOM)
+ else if (rc_scaling == RefCoordScaling::Com)
{
/* Subtract the center of mass */
xp[i][j] -= com[j];
}
}
- if (rc_scaling == erscCOM)
+ if (rc_scaling == RefCoordScaling::Com)
{
/* Convert the COM from Cartesian to crystal coordinates */
for (int j = 0; j < npbcdim; j++)
gmx::ArrayRef<const MoleculeInformation> mi,
const char* fnA,
const char* fnB,
- int rc_scaling,
- int ePBC,
+ RefCoordScaling rc_scaling,
+ PbcType pbcType,
rvec com,
rvec comB,
- warninp* wi)
+ warninp* wi,
+ const gmx::MDLogger& logger)
{
- read_posres(mtop, mi, FALSE, fnA, rc_scaling, ePBC, com, wi);
+ read_posres(mtop, mi, FALSE, fnA, rc_scaling, pbcType, com, wi, logger);
/* It is safer to simply read the b-state posres rather than trying
* to be smart and copy the positions.
*/
- read_posres(mtop, mi, TRUE, fnB, rc_scaling, ePBC, comB, wi);
+ read_posres(mtop, mi, TRUE, fnB, rc_scaling, pbcType, comB, wi, logger);
}
-static void set_wall_atomtype(PreprocessingAtomTypes* at, t_gromppopts* opts, t_inputrec* ir, warninp* wi)
+static void set_wall_atomtype(PreprocessingAtomTypes* at,
+ t_gromppopts* opts,
+ t_inputrec* ir,
+ warninp* wi,
+ const gmx::MDLogger& logger)
{
- int i;
- char warn_buf[STRLEN];
+ int i;
if (ir->nwall > 0)
{
- fprintf(stderr, "Searching the wall atom type(s)\n");
+ GMX_LOG(logger.info).asParagraph().appendTextFormatted("Searching the wall atom type(s)");
}
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())
{
- sprintf(warn_buf, "Specified wall atom type %s is not defined", opts->wall_atomtype[i]);
- warning_error(wi, warn_buf);
+ 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++;
}
}
switch (nmass)
{
- case 0: nrdf = 0; break;
+ case 0: // Fall through intended
case 1: nrdf = 0; break;
case 2: nrdf = 1; break;
default: nrdf = nmass * 3 - 6; break;
for (i = 0; i < 2 * grid_spacing; i++)
{
- spline1d(dx, &(tmp_grid[2 * grid_spacing * i]), 2 * grid_spacing, tmp_u.data(),
+ spline1d(dx,
+ &(tmp_grid[2 * grid_spacing * i]),
+ 2 * grid_spacing,
+ tmp_u.data(),
&(tmp_t2[2 * grid_spacing * i]));
}
for (k = 0; k < 2 * grid_spacing; k++)
{
- interpolate1d(xmin, dx, &(tmp_grid[2 * grid_spacing * k]),
- &(tmp_t2[2 * grid_spacing * k]), psi, &tmp_yy[k], &tmp_y1[k]);
+ interpolate1d(xmin,
+ dx,
+ &(tmp_grid[2 * grid_spacing * k]),
+ &(tmp_t2[2 * grid_spacing * k]),
+ psi,
+ &tmp_yy[k],
+ &tmp_y1[k]);
}
spline1d(dx, tmp_yy.data(), 2 * grid_spacing, tmp_u.data(), tmp_u2.data());
static int count_constraints(const gmx_mtop_t* mtop, gmx::ArrayRef<const MoleculeInformation> mi, warninp* wi)
{
- int count, count_mol;
- char buf[STRLEN];
+ int count, count_mol;
count = 0;
for (const gmx_molblock_t& molb : mtop->molblock)
if (count_mol > nrdf_internal(&mi[molb.type].atoms))
{
- sprintf(buf,
+ std::string warningMessage = gmx::formatString(
"Molecule type '%s' has %d constraints.\n"
"For stability and efficiency there should not be more constraints than "
"internal number of degrees of freedom: %d.\n",
- *mi[molb.type].name, count_mol, nrdf_internal(&mi[molb.type].atoms));
- warning(wi, buf);
+ *mi[molb.type].name,
+ count_mol,
+ nrdf_internal(&mi[molb.type].atoms));
+ warning(wi, warningMessage.c_str());
}
count += molb.nmol * count_mol;
}
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)
if (bNoCoupl)
{
- char buf[STRLEN];
-
- sprintf(buf,
+ std::string warningMessage = gmx::formatString(
"Some temperature coupling groups do not use temperature coupling. We will assume "
"their temperature is not more than %.3f K. If their temperature is higher, the "
"energy error and the Verlet buffer might be underestimated.",
ref_t);
- warning(wi, buf);
+ warning(wi, warningMessage.c_str());
}
return ref_t;
/* Checks if there are unbound atoms in moleculetype molt.
* Prints a note for each unbound atoms and a warning if any is present.
*/
-static void checkForUnboundAtoms(const gmx_moltype_t* molt, gmx_bool bVerbose, warninp* wi)
+static void checkForUnboundAtoms(const gmx_moltype_t* molt, gmx_bool bVerbose, warninp* wi, const gmx::MDLogger& logger)
{
const t_atoms* atoms = &molt->atoms;
for (int ftype = 0; ftype < F_NRE; ftype++)
{
- if (((interaction_function[ftype].flags & IF_BOND) && ftype != F_CONNBONDS)
+ if (((interaction_function[ftype].flags & IF_BOND) && NRAL(ftype) == 2 && ftype != F_CONNBONDS)
|| (interaction_function[ftype].flags & IF_CONSTRAINT) || ftype == F_SETTLE)
{
const InteractionList& il = molt->ilist[ftype];
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)
{
- fprintf(stderr,
- "\nAtom %d '%s' in moleculetype '%s' is not bound by a potential or "
- "constraint to any other atom in the same moleculetype.\n",
- a + 1, *atoms->atomname[a], *molt->name);
+ GMX_LOG(logger.warning)
+ .asParagraph()
+ .appendTextFormatted(
+ "Atom %d '%s' in moleculetype '%s' is not bound by a potential or "
+ "constraint to any other atom in the same moleculetype.",
+ a + 1,
+ *atoms->atomname[a],
+ *molt->name);
}
numDanglingAtoms++;
}
if (numDanglingAtoms > 0)
{
- char buf[STRLEN];
- sprintf(buf,
+ std::string warningMessage = gmx::formatString(
"In moleculetype '%s' %d atoms are not bound by a potential or constraint to any "
"other atom in the same moleculetype. Although technically this might not cause "
"issues in a simulation, this often means that the user forgot to add a "
"bond/potential/constraint or put multiple molecules in the same moleculetype "
"definition by mistake. Run with -v to get information for each atom.",
- *molt->name, numDanglingAtoms);
- warning_note(wi, buf);
+ *molt->name,
+ numDanglingAtoms);
+ warning_note(wi, warningMessage.c_str());
}
}
/* Checks all moleculetypes for unbound atoms */
-static void checkForUnboundAtoms(const gmx_mtop_t* mtop, gmx_bool bVerbose, warninp* wi)
+static void checkForUnboundAtoms(const gmx_mtop_t* mtop, gmx_bool bVerbose, warninp* wi, const gmx::MDLogger& logger)
{
for (const gmx_moltype_t& molt : mtop->moltype)
{
- checkForUnboundAtoms(&molt, bVerbose, wi);
+ checkForUnboundAtoms(&molt, bVerbose, wi, logger);
}
}
gmx::ArrayRef<const t_iparams> iparams,
real massFactorThreshold)
{
- if (molt.ilist[F_CONSTR].size() == 0 && molt.ilist[F_CONSTRNC].size() == 0)
+ if (molt.ilist[F_CONSTR].empty() && molt.ilist[F_CONSTRNC].empty())
{
return false;
}
const t_atom* atom = molt.atoms.atom;
- t_blocka atomToConstraints =
+ const auto atomToConstraints =
gmx::make_at2con(molt, iparams, gmx::FlexibleConstraintTreatment::Exclude);
bool haveDecoupledMode = false;
int a1 = il.iatoms[1 + i + 1];
int a2 = il.iatoms[1 + i + 2];
if ((atom[a0].m > atom[a2].m * massFactorThreshold || atom[a2].m > atom[a0].m * massFactorThreshold)
- && atomToConstraints.index[a0 + 1] - atomToConstraints.index[a0] == 1
- && atomToConstraints.index[a2 + 1] - atomToConstraints.index[a2] == 1
- && atomToConstraints.index[a1 + 1] - atomToConstraints.index[a1] >= 3)
+ && atomToConstraints[a0].ssize() == 1 && atomToConstraints[a2].ssize() == 1
+ && atomToConstraints[a1].ssize() >= 3)
{
- int constraint0 = atomToConstraints.a[atomToConstraints.index[a0]];
- int constraint2 = atomToConstraints.a[atomToConstraints.index[a2]];
+ int constraint0 = atomToConstraints[a0][0];
+ int constraint2 = atomToConstraints[a2][0];
bool foundAtom0 = false;
bool foundAtom2 = false;
- for (int conIndex = atomToConstraints.index[a1];
- conIndex < atomToConstraints.index[a1 + 1]; conIndex++)
+ for (const int constraint : atomToConstraints[a1])
{
- if (atomToConstraints.a[conIndex] == constraint0)
+ if (constraint == constraint0)
{
foundAtom0 = true;
}
- if (atomToConstraints.a[conIndex] == constraint2)
+ if (constraint == constraint2)
{
foundAtom2 = true;
}
}
}
- done_blocka(&atomToConstraints);
-
return haveDecoupledMode;
}
/*! \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)
const int lincsOrderThreshold = 4;
const real shakeToleranceThreshold = 0.005 * ir->delta_t;
- bool lincsWithSufficientTolerance = (ir->eConstrAlg == econtLINCS && ir->nLincsIter >= lincsIterationThreshold
- && ir->nProjOrder >= lincsOrderThreshold);
- bool shakeWithSufficientTolerance =
- (ir->eConstrAlg == econtSHAKE && ir->shake_tol <= 1.1 * shakeToleranceThreshold);
- if (ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol <= 1.1 * bufferToleranceThreshold
+ bool lincsWithSufficientTolerance =
+ (ir->eConstrAlg == ConstraintAlgorithm::Lincs
+ && ir->nLincsIter >= lincsIterationThreshold && ir->nProjOrder >= lincsOrderThreshold);
+ bool shakeWithSufficientTolerance = (ir->eConstrAlg == ConstraintAlgorithm::Shake
+ && ir->shake_tol <= 1.1 * shakeToleranceThreshold);
+ if (ir->cutoff_scheme == CutoffScheme::Verlet && ir->verletbuf_tol <= 1.1 * bufferToleranceThreshold
&& (lincsWithSufficientTolerance || shakeWithSufficientTolerance))
{
return;
"and with masses that differ by more than a factor of %g. This means "
"that there are likely dynamic modes that are only very weakly coupled.",
massFactorThreshold);
- if (ir->cutoff_scheme == ecutsVERLET)
+ if (ir->cutoff_scheme == CutoffScheme::Verlet)
{
message += gmx::formatString(
" To ensure good equipartitioning, you need to either not use "
"hydrogens) or use integrator = %s or decrease one or more tolerances: "
"verlet-buffer-tolerance <= %g, LINCS iterations >= %d, LINCS order "
">= %d or SHAKE tolerance <= %g",
- ei_names[eiSD1], bufferToleranceThreshold, lincsIterationThreshold,
- lincsOrderThreshold, shakeToleranceThreshold);
+ enumValueToString(IntegrationAlgorithm::SD1),
+ bufferToleranceThreshold,
+ lincsIterationThreshold,
+ lincsOrderThreshold,
+ shakeToleranceThreshold);
}
else
{
" To ensure good equipartitioning, we suggest to switch to the %s "
"cutoff-scheme, since that allows for better control over the Verlet "
"buffer size and thus over the energy drift.",
- ecutscheme_names[ecutsVERLET]);
+ enumValueToString(CutoffScheme::Verlet));
}
warning(wi, message);
}
}
-static void set_verlet_buffer(const gmx_mtop_t* mtop, t_inputrec* ir, real buffer_temp, matrix box, warninp* wi)
+static void set_verlet_buffer(const gmx_mtop_t* mtop,
+ t_inputrec* ir,
+ real buffer_temp,
+ matrix box,
+ warninp* wi,
+ const gmx::MDLogger& logger)
{
- char warn_buf[STRLEN];
-
- printf("Determining Verlet buffer for a tolerance of %g kJ/mol/ps at %g K\n", ir->verletbuf_tol,
- buffer_temp);
+ GMX_LOG(logger.info)
+ .asParagraph()
+ .appendTextFormatted(
+ "Determining Verlet buffer for a tolerance of %g kJ/mol/ps at %g K",
+ ir->verletbuf_tol,
+ buffer_temp);
/* Calculate the buffer size for simple atom vs atoms list */
VerletbufListSetup listSetup1x1;
listSetup1x1.cluster_size_i = 1;
listSetup1x1.cluster_size_j = 1;
- const real rlist_1x1 = calcVerletBufferSize(*mtop, det(box), *ir, ir->nstlist, ir->nstlist - 1,
- buffer_temp, listSetup1x1);
+ const real rlist_1x1 = calcVerletBufferSize(
+ *mtop, det(box), *ir, ir->nstlist, ir->nstlist - 1, buffer_temp, listSetup1x1);
/* Set the pair-list buffer size in ir */
VerletbufListSetup listSetup4x4 = verletbufGetSafeListSetup(ListSetupType::CpuNoSimd);
- ir->rlist = calcVerletBufferSize(*mtop, det(box), *ir, ir->nstlist, ir->nstlist - 1,
- buffer_temp, listSetup4x4);
+ ir->rlist = calcVerletBufferSize(
+ *mtop, det(box), *ir, ir->nstlist, ir->nstlist - 1, buffer_temp, listSetup4x4);
- const int n_nonlin_vsite = countNonlinearVsites(*mtop);
+ const int n_nonlin_vsite = gmx::countNonlinearVsites(*mtop);
if (n_nonlin_vsite > 0)
{
- sprintf(warn_buf,
+ std::string warningMessage = gmx::formatString(
"There are %d non-linear virtual site constructions. Their contribution to the "
"energy error is approximated. In most cases this does not affect the error "
"significantly.",
n_nonlin_vsite);
- warning_note(wi, warn_buf);
- }
-
- printf("Calculated rlist for %dx%d atom pair-list as %.3f nm, buffer size %.3f nm\n", 1, 1,
- rlist_1x1, rlist_1x1 - std::max(ir->rvdw, ir->rcoulomb));
-
- printf("Set rlist, assuming %dx%d atom pair-list, to %.3f nm, buffer size %.3f nm\n",
- listSetup4x4.cluster_size_i, listSetup4x4.cluster_size_j, ir->rlist,
- ir->rlist - std::max(ir->rvdw, ir->rcoulomb));
-
- printf("Note that mdrun will redetermine rlist based on the actual pair-list setup\n");
-
- if (gmx::square(ir->rlist) >= max_cutoff2(ir->ePBC, box))
+ warning_note(wi, warningMessage);
+ }
+
+ GMX_LOG(logger.info)
+ .asParagraph()
+ .appendTextFormatted(
+ "Calculated rlist for %dx%d atom pair-list as %.3f nm, buffer size %.3f nm",
+ 1,
+ 1,
+ rlist_1x1,
+ rlist_1x1 - std::max(ir->rvdw, ir->rcoulomb));
+
+ GMX_LOG(logger.info)
+ .asParagraph()
+ .appendTextFormatted(
+ "Set rlist, assuming %dx%d atom pair-list, to %.3f nm, buffer size %.3f nm",
+ listSetup4x4.cluster_size_i,
+ listSetup4x4.cluster_size_j,
+ ir->rlist,
+ ir->rlist - std::max(ir->rvdw, ir->rcoulomb));
+
+ GMX_LOG(logger.info)
+ .asParagraph()
+ .appendTextFormatted(
+ "Note that mdrun will redetermine rlist based on the actual pair-list setup");
+
+ if (gmx::square(ir->rlist) >= max_cutoff2(ir->pbcType, box))
{
gmx_fatal(FARGS,
"The pair-list cut-off (%g nm) is longer than half the shortest box vector or "
"longer than the smallest box diagonal element (%g nm). Increase the box size or "
"decrease nstlist or increase verlet-buffer-tolerance.",
- ir->rlist, std::sqrt(max_cutoff2(ir->ePBC, box)));
+ ir->rlist,
+ std::sqrt(max_cutoff2(ir->pbcType, box)));
}
}
"Then a coordinate file is read and velocities can be generated",
"from a Maxwellian distribution if requested.",
"[THISMODULE] also reads parameters for [gmx-mdrun] ",
- "(eg. number of MD steps, time step, cut-off), and others such as",
- "NEMD parameters, which are corrected so that the net acceleration",
- "is zero.",
+ "(eg. number of MD steps, time step, cut-off).",
"Eventually a binary file is produced that can serve as the sole input",
"file for the MD program.[PAR]",
"interpret the output messages before attempting to bypass them with",
"this option."
};
- t_gromppopts* opts;
std::vector<MoleculeInformation> mi;
std::unique_ptr<MoleculeInformation> intermolecular_interactions;
- int nvsite, comb;
+ int nvsite;
+ CombinationRule comb;
real fudgeQQ;
double reppow;
const char* mdparin;
bool bNeedVel, bGenVel;
- gmx_bool have_atomnumber;
gmx_output_env_t* oenv;
gmx_bool bVerbose = FALSE;
warninp* wi;
- char warn_buf[STRLEN];
t_filenm fnm[] = { { efMDP, nullptr, nullptr, ffREAD },
{ efMDP, "-po", "mdout", ffWRITE },
{ 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 } };
+ { 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::MDModules mdModules;
t_inputrec irInstance;
t_inputrec* ir = &irInstance;
- snew(opts, 1);
+ t_gromppopts optsInstance;
+ t_gromppopts* opts = &optsInstance;
snew(opts->include, STRLEN);
snew(opts->define, STRLEN);
+ gmx::LoggerBuilder builder;
+ builder.addTargetStream(gmx::MDLogger::LogLevel::Info, &gmx::TextOutputFile::standardOutput());
+ builder.addTargetStream(gmx::MDLogger::LogLevel::Warning, &gmx::TextOutputFile::standardError());
+ gmx::LoggerOwner logOwner(builder.build());
+ const gmx::MDLogger logger(logOwner.logger());
+
+
wi = init_warning(TRUE, maxwarn);
/* PARAMETER file processing */
}
GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
+ // 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)
{
- fprintf(stderr, "checking input for internal consistency...\n");
+ GMX_LOG(logger.info)
+ .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)
{
ir->ld_seed = static_cast<int>(gmx::makeRandomSeed());
- fprintf(stderr, "Setting the LD random seed to %" PRId64 "\n", ir->ld_seed);
+ GMX_LOG(logger.info)
+ .asParagraph()
+ .appendTextFormatted("Setting the LD random seed to %" PRId64 "", ir->ld_seed);
}
if (ir->expandedvals->lmc_seed == -1)
{
ir->expandedvals->lmc_seed = static_cast<int>(gmx::makeRandomSeed());
- fprintf(stderr, "Setting the lambda MC random seed to %d\n", ir->expandedvals->lmc_seed);
+ GMX_LOG(logger.info)
+ .asParagraph()
+ .appendTextFormatted("Setting the lambda MC random seed to %d", ir->expandedvals->lmc_seed);
}
bNeedVel = EI_STATE_VELOCITY(ir->eI);
bGenVel = (bNeedVel && opts->bGenVel);
if (bGenVel && ir->bContinuation)
{
- sprintf(warn_buf,
+ std::string warningMessage = gmx::formatString(
"Generating velocities is inconsistent with attempting "
"to continue a previous run. Choose only one of "
"gen-vel = yes and continuation = yes.");
- warning_error(wi, warn_buf);
+ warning_error(wi, warningMessage);
}
std::array<InteractionsOfType, F_NRE> interactions;
}
t_state state;
- new_status(fn, opt2fn_null("-pp", NFILE, fnm), opt2fn("-c", NFILE, fnm), opts, ir, bZero,
- bGenVel, bVerbose, &state, &atypes, &sys, &mi, &intermolecular_interactions,
- interactions, &comb, &reppow, &fudgeQQ, opts->bMorse, wi);
+ new_status(fn,
+ opt2fn_null("-pp", NFILE, fnm),
+ opt2fn("-c", NFILE, fnm),
+ opts,
+ ir,
+ bZero,
+ bGenVel,
+ bVerbose,
+ &state,
+ &atypes,
+ &sys,
+ &mi,
+ &intermolecular_interactions,
+ interactions,
+ &comb,
+ &reppow,
+ &fudgeQQ,
+ opts->bMorse,
+ wi,
+ logger);
if (debug)
{
/* set parameters for virtual site construction (not for vsiten) */
for (size_t mt = 0; mt < sys.moltype.size(); mt++)
{
- nvsite += set_vsites(bVerbose, &sys.moltype[mt].atoms, &atypes, mi[mt].interactions);
+ nvsite += set_vsites(bVerbose, &sys.moltype[mt].atoms, &atypes, mi[mt].interactions, logger);
}
/* now throw away all obsolete bonds, angles and dihedrals: */
/* note: constraints are ALWAYS removed */
{
for (size_t mt = 0; mt < sys.moltype.size(); mt++)
{
- clean_vsite_bondeds(mi[mt].interactions, sys.moltype[mt].atoms.nr, bRmVSBds);
+ clean_vsite_bondeds(mi[mt].interactions, sys.moltype[mt].atoms.nr, bRmVSBds, logger);
}
}
- if ((count_constraints(&sys, mi, wi) != 0) && (ir->eConstrAlg == econtSHAKE))
+ if ((count_constraints(&sys, mi, wi) != 0) && (ir->eConstrAlg == ConstraintAlgorithm::Shake))
{
- if (ir->eI == eiCG || ir->eI == eiLBFGS)
+ if (ir->eI == IntegrationAlgorithm::CG || ir->eI == IntegrationAlgorithm::LBFGS)
{
- sprintf(warn_buf, "Can not do %s with %s, use %s", EI(ir->eI),
- econstr_names[econtSHAKE], econstr_names[econtLINCS]);
- warning_error(wi, warn_buf);
+ std::string warningMessage =
+ gmx::formatString("Can not do %s with %s, use %s",
+ enumValueToString(ir->eI),
+ enumValueToString(ConstraintAlgorithm::Shake),
+ enumValueToString(ConstraintAlgorithm::Lincs));
+ warning_error(wi, warningMessage);
}
if (ir->bPeriodicMols)
{
- sprintf(warn_buf, "Can not do periodic molecules with %s, use %s",
- econstr_names[econtSHAKE], econstr_names[econtLINCS]);
- warning_error(wi, warn_buf);
+ std::string warningMessage =
+ gmx::formatString("Can not do periodic molecules with %s, use %s",
+ enumValueToString(ConstraintAlgorithm::Shake),
+ enumValueToString(ConstraintAlgorithm::Lincs));
+ warning_error(wi, warningMessage);
}
}
- if (EI_SD(ir->eI) && ir->etc != etcNO)
+ if (EI_SD(ir->eI) && ir->etc != TemperatureCoupling::No)
{
warning_note(wi, "Temperature coupling is ignored with SD integrators.");
}
- /* If we are doing QM/MM, check that we got the atom numbers */
- have_atomnumber = TRUE;
- for (gmx::index i = 0; i < gmx::ssize(atypes); i++)
- {
- have_atomnumber = have_atomnumber && (atypes.atomNumberFromAtomType(i) >= 0);
- }
- if (!have_atomnumber && ir->bQMMM)
- {
- warning_error(
- wi,
- "\n"
- "It appears as if you are trying to run a QM/MM calculation, but the force\n"
- "field you are using does not contain atom numbers fields. This is an\n"
- "optional field (introduced in GROMACS 3.3) for general runs, but mandatory\n"
- "for QM/MM. The good news is that it is easy to add - put the atom number as\n"
- "an integer just before the mass column in ffXXXnb.itp.\n"
- "NB: United atoms have the same atom numbers as normal ones.\n\n");
- }
-
/* Check for errors in the input now, since they might cause problems
* during processing further down.
*/
if (nint_ftype(&sys, mi, F_POSRES) > 0 || nint_ftype(&sys, mi, F_FBPOSRES) > 0)
{
- if (ir->epc == epcPARRINELLORAHMAN || ir->epc == epcMTTK)
+ if (ir->epc == PressureCoupling::ParrinelloRahman || ir->epc == PressureCoupling::Mttk)
{
- sprintf(warn_buf,
+ std::string warningMessage = gmx::formatString(
"You are combining position restraints with %s pressure coupling, which can "
"lead to instabilities. If you really want to combine position restraints with "
"pressure coupling, we suggest to use %s pressure coupling instead.",
- EPCOUPLTYPE(ir->epc), EPCOUPLTYPE(epcBERENDSEN));
- warning_note(wi, warn_buf);
+ enumValueToString(ir->epc),
+ enumValueToString(PressureCoupling::Berendsen));
+ warning_note(wi, warningMessage);
}
const char* fn = opt2fn("-r", NFILE, fnm);
if (bVerbose)
{
- fprintf(stderr, "Reading position restraint coords from %s", fn);
- if (strcmp(fn, fnB) == 0)
+ std::string message = gmx::formatString("Reading position restraint coords from %s", fn);
+ if (strcmp(fn, fnB) != 0)
{
- fprintf(stderr, "\n");
- }
- else
- {
- fprintf(stderr, " and %s\n", fnB);
+ message += gmx::formatString(" and %s", fnB);
}
+ GMX_LOG(logger.info).asParagraph().appendText(message);
}
- gen_posres(&sys, mi, fn, fnB, ir->refcoord_scaling, ir->ePBC, ir->posres_com, ir->posres_comB, wi);
+ gen_posres(&sys, mi, fn, fnB, ir->refcoord_scaling, ir->pbcType, ir->posres_com, ir->posres_comB, wi, logger);
}
/* If we are using CMAP, setup the pre-interpolation grid */
if (interactions[F_CMAP].ncmap() > 0)
{
- init_cmap_grid(&sys.ffparams.cmap_grid, interactions[F_CMAP].cmapAngles,
+ init_cmap_grid(&sys.ffparams.cmap_grid,
+ interactions[F_CMAP].cmapAngles,
interactions[F_CMAP].cmakeGridSpacing);
- setup_cmap(interactions[F_CMAP].cmakeGridSpacing, interactions[F_CMAP].cmapAngles,
- interactions[F_CMAP].cmap, &sys.ffparams.cmap_grid);
+ setup_cmap(interactions[F_CMAP].cmakeGridSpacing,
+ interactions[F_CMAP].cmapAngles,
+ interactions[F_CMAP].cmap,
+ &sys.ffparams.cmap_grid);
}
- set_wall_atomtype(&atypes, opts, ir, wi);
+ set_wall_atomtype(&atypes, opts, ir, wi, logger);
if (bRenum)
{
atypes.renumberTypes(interactions, &sys, ir->wall_atomtype, bVerbose);
if (bVerbose)
{
- fprintf(stderr, "converting bonded parameters...\n");
+ GMX_LOG(logger.info).asParagraph().appendTextFormatted("converting bonded parameters...");
}
const int ntype = atypes.size();
- convertInteractionsOfType(ntype, interactions, mi, intermolecular_interactions.get(), comb,
- reppow, fudgeQQ, &sys);
+ convertInteractionsOfType(
+ ntype, interactions, mi, intermolecular_interactions.get(), comb, reppow, fudgeQQ, &sys);
if (debug)
{
/* set ptype to VSite for virtual sites */
for (gmx_moltype_t& moltype : sys.moltype)
{
- set_vsites_ptype(FALSE, &moltype);
+ set_vsites_ptype(FALSE, &moltype, logger);
}
if (debug)
{
/* check masses */
check_mol(&sys, wi);
- checkForUnboundAtoms(&sys, bVerbose, wi);
+ if (haveFepPerturbedMassesInSettles(sys))
+ {
+ warning_error(wi,
+ "SETTLE is not implemented for atoms whose mass is perturbed. "
+ "You might instead use normal constraints.");
+ }
- if (EI_DYNAMICS(ir->eI) && ir->eI != eiBD)
+ checkForUnboundAtoms(&sys, bVerbose, wi, logger);
+
+ if (EI_DYNAMICS(ir->eI) && ir->eI != IntegrationAlgorithm::BD)
{
check_bonds_timestep(&sys, ir->delta_t, wi);
}
if (bVerbose)
{
- fprintf(stderr, "initialising group options...\n");
+ 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 == ecutsVERLET && ir->verletbuf_tol > 0)
+ if (ir->cutoff_scheme == CutoffScheme::Verlet && ir->verletbuf_tol > 0)
{
if (EI_DYNAMICS(ir->eI) && inputrec2nboundeddim(ir) == 3)
{
real buffer_temp;
- if (EI_MD(ir->eI) && ir->etc == etcNO)
+ if (EI_MD(ir->eI) && ir->etc == TemperatureCoupling::No)
{
if (bGenVel)
{
}
if (buffer_temp > 0)
{
- sprintf(warn_buf,
+ std::string warningMessage = gmx::formatString(
"NVE simulation: will use the initial temperature of %.3f K for "
"determining the Verlet buffer size",
buffer_temp);
- warning_note(wi, warn_buf);
+ warning_note(wi, warningMessage);
}
else
{
- sprintf(warn_buf,
+ std::string warningMessage = gmx::formatString(
"NVE simulation with an initial temperature of zero: will use a Verlet "
"buffer of %d%%. Check your energy drift!",
gmx::roundToInt(verlet_buffer_ratio_NVE_T0 * 100));
- warning_note(wi, warn_buf);
+ warning_note(wi, warningMessage);
}
}
else
buffer_temp = get_max_reference_temp(ir, wi);
}
- if (EI_MD(ir->eI) && ir->etc == etcNO && buffer_temp == 0)
+ if (EI_MD(ir->eI) && ir->etc == TemperatureCoupling::No && buffer_temp == 0)
{
/* NVE with initial T=0: we add a fixed ratio to rlist.
* Since we don't actually use verletbuf_tol, we set it to -1
* Note that we can't warn when nsteps=0, since we don't
* know how many steps the user intends to run.
*/
- if (EI_MD(ir->eI) && ir->etc == etcNO && ir->nstlist > 1 && ir->nsteps > 0)
+ if (EI_MD(ir->eI) && ir->etc == TemperatureCoupling::No && ir->nstlist > 1 && ir->nsteps > 0)
{
const real driftTolerance = 0.01;
/* We use 2 DOF per atom = 2kT pot+kin energy,
* 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)
{
- sprintf(warn_buf,
+ std::string warningMessage = gmx::formatString(
"You are using a Verlet buffer tolerance of %g kJ/mol/ps for an "
"NVE simulation of length %g ps, which can give a final drift of "
"%d%%. For conserving energy to %d%% when using constraints, you "
"might need to set verlet-buffer-tolerance to %.1e.",
- ir->verletbuf_tol, ir->nsteps * ir->delta_t,
+ ir->verletbuf_tol,
+ ir->nsteps * ir->delta_t,
gmx::roundToInt(ir->verletbuf_tol / totalEnergyDriftPerAtomPerPicosecond * 100),
gmx::roundToInt(100 * driftTolerance),
driftTolerance * totalEnergyDriftPerAtomPerPicosecond);
- warning_note(wi, warn_buf);
+ warning_note(wi, warningMessage);
}
}
- set_verlet_buffer(&sys, ir, buffer_temp, state.box, wi);
+ set_verlet_buffer(&sys, ir, buffer_temp, state.box, wi, logger);
}
}
}
pr_symtab(debug, 0, "After close", &sys.symtab);
}
- /* make exclusions between QM atoms and remove charges if needed */
- if (ir->bQMMM)
+ if (ir->eI == IntegrationAlgorithm::Mimic)
{
- generate_qmexcl(&sys, ir, wi, GmxQmmmMode::GMX_QMMM_ORIGINAL);
- if (ir->QMMMscheme != eQMMMschemeoniom)
- {
- std::vector<int> qmmmAtoms = qmmmAtomIndices(*ir, sys);
- removeQmmmAtomCharges(&sys, qmmmAtoms);
- }
- }
-
- if (ir->eI == eiMimic)
- {
- generate_qmexcl(&sys, ir, wi, GmxQmmmMode::GMX_QMMM_MIMIC);
+ generate_qmexcl(&sys, ir, logger);
}
if (ftp2bSet(efTRN, NFILE, fnm))
{
if (bVerbose)
{
- fprintf(stderr, "getting data from old trajectory ...\n");
+ GMX_LOG(logger.info)
+ .asParagraph()
+ .appendTextFormatted("getting data from old trajectory ...");
}
- cont_status(ftp2fn(efTRN, NFILE, fnm), ftp2fn_null(efEDR, NFILE, fnm), bNeedVel, bGenVel,
- fr_time, ir, &state, &sys, oenv);
+ cont_status(ftp2fn(efTRN, NFILE, fnm),
+ ftp2fn_null(efEDR, NFILE, fnm),
+ bNeedVel,
+ bGenVel,
+ fr_time,
+ ir,
+ &state,
+ &sys,
+ oenv,
+ logger);
}
- if (ir->ePBC == epbcXY && ir->nwall != 2)
+ if (ir->pbcType == PbcType::XY && ir->nwall != 2)
{
clear_rvec(state.box[ZZ]);
}
{
/* 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)
wi, "Some of the Fourier grid sizes are set, but all of them need to be set.");
}
const int minGridSize = minimalPmeGridSize(ir->pme_order);
- calcFftGrid(stdout, scaledBox, ir->fourier_spacing, minGridSize, &(ir->nkx), &(ir->nky),
- &(ir->nkz));
+ calcFftGrid(stdout, scaledBox, ir->fourier_spacing, minGridSize, &(ir->nkx), &(ir->nky), &(ir->nkz));
if (ir->nkx < minGridSize || ir->nky < minGridSize || ir->nkz < minGridSize)
{
warning_error(wi,
/* MRS: eventually figure out better logic for initializing the fep
values that makes declaring the lambda and declaring the state not
potentially conflict if not handled correctly. */
- if (ir->efep != efepNO)
+ if (ir->efep != FreeEnergyPerturbationType::No)
{
state.fep_state = ir->fepvals->init_fep_state;
- for (i = 0; i < efptNR; i++)
+ for (i = 0; i < static_cast<int>(FreeEnergyPerturbationCouplingType::Count); i++)
{
/* init_lambda trumps state definitions*/
if (ir->fepvals->init_lambda >= 0)
}
else
{
- if (ir->fepvals->all_lambda[i] == nullptr)
+ if (ir->fepvals->all_lambda[i].empty())
{
gmx_fatal(FARGS, "Values of lambda not set for a free energy calculation!");
}
if (ir->bPull)
{
- pull = set_pull_init(ir, &sys, state.x.rvec_array(), state.box, state.lambda[efptMASS], wi);
+ pull = set_pull_init(
+ ir, sys, state.x, state.box, state.lambda[FreeEnergyPerturbationCouplingType::Mass], wi);
}
/* Modules that supply external potential for pull coordinates
if (ir->bDoAwh)
{
tensor compressibility = { { 0 } };
- if (ir->epc != epcNO)
+ if (ir->epc != PressureCoupling::No)
{
copy_mat(ir->compress, compressibility);
}
- setStateDependentAwhParams(ir->awhParams, ir->pull, pull, state.box, ir->ePBC,
- compressibility, &ir->opts, wi);
+ setStateDependentAwhParams(
+ 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,
+ sys,
+ wi);
}
if (ir->bPull)
if (ir->bRot)
{
- set_reference_positions(ir->rot, state.x.rvec_array(), state.box,
- opt2fn("-ref", NFILE, fnm), opt2bSet("-ref", NFILE, fnm), wi);
+ set_reference_positions(ir->rot,
+ state.x.rvec_array(),
+ state.box,
+ opt2fn("-ref", NFILE, fnm),
+ opt2bSet("-ref", NFILE, fnm),
+ wi);
}
/* reset_multinr(sys); */
if (EEL_PME(ir->coulombtype))
{
float ratio = pme_load_estimate(sys, *ir, state.box);
- fprintf(stderr, "Estimate for the relative computational load of the PME mesh part: %.2f\n", ratio);
+ GMX_LOG(logger.info)
+ .asParagraph()
+ .appendTextFormatted(
+ "Estimate for the relative computational load of the PME mesh part: %.2f", ratio);
/* With free energy we might need to do PME both for the A and B state
* charges. This will double the cost, but the optimal performance will
* then probably be at a slightly larger cut-off and grid spacing.
*/
- if ((ir->efep == efepNO && ratio > 1.0 / 2.0) || (ir->efep != efepNO && ratio > 2.0 / 3.0))
+ if ((ir->efep == FreeEnergyPerturbationType::No && ratio > 1.0 / 2.0)
+ || (ir->efep != FreeEnergyPerturbationType::No && ratio > 2.0 / 3.0))
{
warning_note(
wi,
"The optimal PME mesh load for parallel simulations is below 0.5\n"
"and for highly parallel simulations between 0.25 and 0.33,\n"
"for higher performance, increase the cut-off and the PME grid spacing.\n");
- if (ir->efep != efepNO)
+ if (ir->efep != FreeEnergyPerturbationType::No)
{
warning_note(wi,
"For free energy simulations, the optimal load limit increases from "
}
{
- char warn_buf[STRLEN];
- double cio = compute_io(ir, sys.natoms, sys.groups, F_NRE, 1);
- sprintf(warn_buf, "This run will generate roughly %.0f Mb of data", cio);
+ double cio = compute_io(ir, sys.natoms, sys.groups, F_NRE, 1);
+ std::string warningMessage =
+ gmx::formatString("This run will generate roughly %.0f Mb of data", cio);
if (cio > 2000)
{
set_warning_line(wi, mdparin, -1);
- warning_note(wi, warn_buf);
+ warning_note(wi, warningMessage);
}
else
{
- printf("%s\n", warn_buf);
+ GMX_LOG(logger.info).asParagraph().appendText(warningMessage);
}
}
+ // 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().notifier_.notify(internalParameterBuilder.rootObject());
+ mdModules.notifiers().preProcessingNotifier_.notify(internalParameterBuilder.rootObject());
ir->internalParameters =
std::make_unique<gmx::KeyValueTreeObject>(internalParameterBuilder.build());
}
+ if (ir->comm_mode != ComRemovalAlgorithm::No)
+ {
+ const int nstglobalcomm = computeGlobalCommunicationPeriod(ir);
+ if (ir->nstcomm % nstglobalcomm != 0)
+ {
+ warning_note(
+ wi,
+ gmx::formatString(
+ "COM removal frequency is set to (%d).\n"
+ "Other settings require a global communication frequency of %d.\n"
+ "Note that this will require additional global communication steps,\n"
+ "which will reduce performance when using multiple ranks.\n"
+ "Consider setting nstcomm to a multiple of %d.",
+ ir->nstcomm,
+ nstglobalcomm,
+ nstglobalcomm));
+ }
+ }
+
if (bVerbose)
{
- fprintf(stderr, "writing run input file...\n");
+ GMX_LOG(logger.info).asParagraph().appendTextFormatted("writing run input file...");
}
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);
+ sfree(opts->couple_moltype);
+
for (auto& mol : mi)
{
// Some of the contents of molinfo have been stolen, so
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff