*/
#include "gmxpre.h"
-#include "gromacs/utility/enumerationhelpers.h"
#include "readir.h"
#include <cctype>
#include "gromacs/gmxpreprocess/toputil.h"
#include "gromacs/math/functions.h"
#include "gromacs/math/units.h"
+#include "gromacs/math/utilities.h"
#include "gromacs/math/vec.h"
#include "gromacs/mdlib/calc_verletbuf.h"
+#include "gromacs/mdlib/vcm.h"
#include "gromacs/mdrun/mdmodules.h"
+#include "gromacs/mdtypes/awh_params.h"
#include "gromacs/mdtypes/inputrec.h"
#include "gromacs/mdtypes/md_enums.h"
#include "gromacs/mdtypes/multipletimestepping.h"
#include "gromacs/utility/keyvaluetreebuilder.h"
#include "gromacs/utility/keyvaluetreemdpwriter.h"
#include "gromacs/utility/keyvaluetreetransform.h"
-#include "gromacs/utility/mdmodulenotification.h"
+#include "gromacs/utility/mdmodulesnotifiers.h"
#include "gromacs/utility/smalloc.h"
#include "gromacs/utility/strconvert.h"
#include "gromacs/utility/stringcompare.h"
#include "gromacs/utility/stringutil.h"
#include "gromacs/utility/textwriter.h"
-#define MAXPTR 254
#define NOGID 255
+using gmx::BasicVector;
+
/* Resource parameters
* Do not change any of these until you read the instruction
* in readinp.h. Some cpp's do not take spaces after the backslash
struct gmx_inputrec_strings
{
- char tcgrps[STRLEN], tau_t[STRLEN], ref_t[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 tcgrps[STRLEN], tau_t[STRLEN], ref_t[STRLEN], accelerationGroups[STRLEN],
+ acceleration[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];
gmx::EnumerationArray<FreeEnergyPerturbationCouplingType, std::string> fep_lambda;
char lambda_weights[STRLEN];
std::vector<std::string> pullGroupNames;
char anneal[STRLEN], anneal_npoints[STRLEN], anneal_time[STRLEN], anneal_temp[STRLEN];
};
+// NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
static gmx_inputrec_strings* inputrecStrings = nullptr;
void init_inputrec_strings()
}
-enum
+//! How to treat coverage of the whole system for a set of atom groupsx
+enum class GroupCoverage
{
- egrptpALL, /* All particles have to be a member of a group. */
- egrptpALL_GENREST, /* A rest group with name is generated for particles *
- * that are not part of any group. */
- egrptpPART, /* As egrptpALL_GENREST, but no name is generated *
- * for the rest group. */
- egrptpONE /* Merge all selected groups into one group, *
- * make a rest group for the remaining particles. */
+ All, //!< All particles have to be a member of a group
+ AllGenerateRest, //<! A rest group with name is generated for particles not part of any group
+ Partial, //<! As \p AllGenerateRest, but no name for the rest group is generated
+ OneGroup //<! Merge all selected groups into one group, make a rest group for the remaining particles
};
+// NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
static const char* constraints[eshNR + 1] = { "none", "h-bonds", "all-bonds",
"h-angles", "all-angles", nullptr };
+// NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
static const char* couple_lam[ecouplamNR + 1] = { "vdw-q", "vdw", "q", "none", nullptr };
static void getSimTemps(int ntemps, t_simtemp* simtemp, gmx::ArrayRef<double> temperature_lambdas)
}
}
-void check_ir(const char* mdparin,
- const gmx::MdModulesNotifier& mdModulesNotifier,
- t_inputrec* ir,
- t_gromppopts* opts,
- warninp_t wi)
+void check_ir(const char* mdparin,
+ const gmx::MDModulesNotifiers& mdModulesNotifiers,
+ t_inputrec* ir,
+ t_gromppopts* opts,
+ warninp_t wi)
/* Check internal consistency.
* NOTE: index groups are not set here yet, don't check things
* like temperature coupling group options here, but in triple_check
check_nst("nstcalcenergy", ir->nstcalcenergy, "nstenergy", &ir->nstenergy, wi);
}
- // Inquire all MdModules, if their parameters match with the energy
+ // Inquire all MDModules, if their parameters match with the energy
// calculation frequency
gmx::EnergyCalculationFrequencyErrors energyCalculationFrequencyErrors(ir->nstcalcenergy);
- mdModulesNotifier.preProcessingNotifications_.notify(&energyCalculationFrequencyErrors);
+ mdModulesNotifiers.preProcessingNotifier_.notify(&energyCalculationFrequencyErrors);
// Emit all errors from the energy calculation frequency checks
for (const std::string& energyFrequencyErrorMessage :
CHECK((fep->all_lambda[i][j] < 0) || (fep->all_lambda[i][j] > 1));
}
}
+
+ if (fep->softcoreFunction == SoftcoreType::Gapsys)
+ {
+ if (fep->scGapsysScaleLinpointQ < 0.0)
+ {
+ sprintf(warn_buf,
+ "sc_scale_linpoint_Q_gapsys is equal %g but must be >= 0",
+ fep->scGapsysScaleLinpointQ);
+ warning_note(wi, warn_buf);
+ }
+
+ if ((fep->scGapsysScaleLinpointLJ < 0.0) || (fep->scGapsysScaleLinpointLJ >= 1.0))
+ {
+ sprintf(warn_buf,
+ "sc_scale_linpoint_LJ_gapsys is equal %g but must be in [0,1) when used "
+ "with "
+ "sc_function=gapsys.",
+ fep->scGapsysScaleLinpointLJ);
+ warning_note(wi, warn_buf);
+ }
+ }
}
if ((ir->bSimTemp) || (ir->efep == FreeEnergyPerturbationType::Expanded))
ir->nstcomm = abs(ir->nstcomm);
}
- if (ir->nstcalcenergy > 0 && ir->nstcomm < ir->nstcalcenergy)
+ if (ir->nstcalcenergy > 0 && ir->nstcomm < ir->nstcalcenergy
+ && ir->comm_mode != ComRemovalAlgorithm::LinearAccelerationCorrection)
{
warning_note(wi,
- "nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting "
- "nstcomm to nstcalcenergy");
- ir->nstcomm = ir->nstcalcenergy;
+ "nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, consider "
+ "setting nstcomm equal to nstcalcenergy for less overhead");
}
if (ir->comm_mode == ComRemovalAlgorithm::Angular)
static void do_simtemp_params(t_inputrec* ir)
{
-
- snew(ir->simtempvals->temperatures, ir->fepvals->n_lambda);
+ ir->simtempvals->temperatures.resize(ir->fepvals->n_lambda);
getSimTemps(ir->fepvals->n_lambda,
ir->simtempvals.get(),
ir->fepvals->all_lambda[FreeEnergyPerturbationCouplingType::Temperature]);
}
}
+static void convertRvecs(warninp_t wi, gmx::ArrayRef<const std::string> inputs, const char* name, rvec* outputs)
+{
+ int i = 0, d = 0;
+ for (const auto& input : inputs)
+ {
+ try
+ {
+ outputs[i][d] = gmx::fromString<real>(input);
+ }
+ catch (gmx::GromacsException&)
+ {
+ auto message = gmx::formatString(
+ "Invalid value for mdp option %s. %s should only consist of real numbers "
+ "separated by spaces.",
+ name,
+ name);
+ warning_error(wi, message);
+ }
+ ++d;
+ if (d == DIM)
+ {
+ d = 0;
+ ++i;
+ }
+ }
+}
+
static void do_wall_params(t_inputrec* ir, char* wall_atomtype, char* wall_density, t_gromppopts* opts, warninp_t wi)
{
opts->wall_atomtype[0] = nullptr;
ir->bDoAwh = (getEnum<Boolean>(&inp, "awh", wi) != Boolean::No);
if (ir->bDoAwh)
{
- ir->awhParams = gmx::readAwhParams(&inp, wi);
+ ir->awhParams = std::make_unique<gmx::AwhParams>(&inp, wi);
}
/* Enforced rotation */
setStringEntry(&inp, "temperature-lambdas", "");
fep->lambda_neighbors = get_eint(&inp, "calc-lambda-neighbors", 1, wi);
setStringEntry(&inp, "init-lambda-weights", inputrecStrings->lambda_weights, nullptr);
- fep->edHdLPrintEnergy = getEnum<FreeEnergyPrintEnergy>(&inp, "dhdl-print-energy", wi);
- fep->sc_alpha = get_ereal(&inp, "sc-alpha", 0.0, wi);
- fep->sc_power = get_eint(&inp, "sc-power", 1, wi);
- fep->sc_r_power = get_ereal(&inp, "sc-r-power", 6.0, wi);
- fep->sc_sigma = get_ereal(&inp, "sc-sigma", 0.3, wi);
- fep->bScCoul = (getEnum<Boolean>(&inp, "sc-coul", wi) != Boolean::No);
- fep->dh_hist_size = get_eint(&inp, "dh_hist_size", 0, wi);
- fep->dh_hist_spacing = get_ereal(&inp, "dh_hist_spacing", 0.1, wi);
- fep->separate_dhdl_file = getEnum<SeparateDhdlFile>(&inp, "separate-dhdl-file", wi);
- fep->dhdl_derivatives = getEnum<DhDlDerivativeCalculation>(&inp, "dhdl-derivatives", wi);
- fep->dh_hist_size = get_eint(&inp, "dh_hist_size", 0, wi);
- fep->dh_hist_spacing = get_ereal(&inp, "dh_hist_spacing", 0.1, wi);
+ fep->edHdLPrintEnergy = getEnum<FreeEnergyPrintEnergy>(&inp, "dhdl-print-energy", wi);
+ fep->softcoreFunction = getEnum<SoftcoreType>(&inp, "sc-function", wi);
+ fep->sc_alpha = get_ereal(&inp, "sc-alpha", 0.0, wi);
+ fep->sc_power = get_eint(&inp, "sc-power", 1, wi);
+ fep->sc_r_power = get_ereal(&inp, "sc-r-power", 6.0, wi);
+ fep->sc_sigma = get_ereal(&inp, "sc-sigma", 0.3, wi);
+ fep->bScCoul = (getEnum<Boolean>(&inp, "sc-coul", wi) != Boolean::No);
+ fep->scGapsysScaleLinpointLJ = get_ereal(&inp, "sc-gapsys-scale-linpoint-lj", 0.85, wi);
+ fep->scGapsysScaleLinpointQ = get_ereal(&inp, "sc-gapsys-scale-linpoint-q", 0.3, wi);
+ fep->scGapsysSigmaLJ = get_ereal(&inp, "sc-gapsys-sigma-lj", 0.3, wi);
+ fep->dh_hist_size = get_eint(&inp, "dh_hist_size", 0, wi);
+ fep->dh_hist_spacing = get_ereal(&inp, "dh_hist_spacing", 0.1, wi);
+ fep->separate_dhdl_file = getEnum<SeparateDhdlFile>(&inp, "separate-dhdl-file", wi);
+ fep->dhdl_derivatives = getEnum<DhDlDerivativeCalculation>(&inp, "dhdl-derivatives", wi);
+ fep->dh_hist_size = get_eint(&inp, "dh_hist_size", 0, wi);
+ fep->dh_hist_spacing = get_ereal(&inp, "dh_hist_spacing", 0.1, wi);
/* Non-equilibrium MD stuff */
printStringNewline(&inp, "Non-equilibrium MD stuff");
+ setStringEntry(&inp, "acc-grps", inputrecStrings->accelerationGroups, nullptr);
+ setStringEntry(&inp, "accelerate", inputrecStrings->acceleration, nullptr);
setStringEntry(&inp, "freezegrps", inputrecStrings->freeze, nullptr);
setStringEntry(&inp, "freezedim", inputrecStrings->frdim, nullptr);
ir->cos_accel = get_ereal(&inp, "cos-acceleration", 0, wi);
ir->userreal4 = get_ereal(&inp, "userreal4", 0, wi);
#undef CTYPE
+ if (mdparout)
{
gmx::TextOutputFile stream(mdparout);
write_inpfile(&stream, mdparout, &inp, FALSE, writeMdpHeader, wi);
if (ir->bDoAwh)
{
- gmx::checkAwhParams(ir->awhParams, ir, wi);
+ gmx::checkAwhParams(*ir->awhParams, *ir, wi);
}
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
- string in a collection), so should be refactored and located more
- centrally. */
-int search_string(const char* s, int ng, char* gn[])
+int search_string(const char* s, int ng, char* const gn[])
{
int i;
}
}
-static void do_numbering(int natoms,
- SimulationGroups* groups,
- gmx::ArrayRef<std::string> groupsFromMdpFile,
- t_blocka* block,
- char* gnames[],
- SimulationAtomGroupType gtype,
- int restnm,
- int grptp,
- bool bVerbose,
- warninp_t wi)
+/*! Creates the groups of atom indices for group type \p gtype
+ *
+ * \param[in] natoms The total number of atoms in the system
+ * \param[in,out] groups Index \p gtype in this list of list of groups will be set
+ * \param[in] groupsFromMdpFile The list of group names set for \p gtype in the mdp file
+ * \param[in] block The list of atom indices for all available index groups
+ * \param[in] gnames The list of names for all available index groups
+ * \param[in] gtype The group type to creates groups for
+ * \param[in] restnm The index of rest group name in \p gnames
+ * \param[in] coverage How to treat coverage of all atoms in the system
+ * \param[in] bVerbose Whether to print when we make a rest group
+ * \param[in,out] wi List of warnings
+ */
+static void do_numbering(const int natoms,
+ SimulationGroups* groups,
+ gmx::ArrayRef<const std::string> groupsFromMdpFile,
+ const t_blocka* block,
+ char* const gnames[],
+ const SimulationAtomGroupType gtype,
+ const int restnm,
+ const GroupCoverage coverage,
+ const bool bVerbose,
+ warninp_t wi)
{
unsigned short* cbuf;
AtomGroupIndices* grps = &(groups->groups[gtype]);
{
/* Lookup the group name in the block structure */
const int gid = search_string(groupsFromMdpFile[i].c_str(), block->nr, gnames);
- if ((grptp != egrptpONE) || (i == 0))
+ if ((coverage != GroupCoverage::OneGroup) || (i == 0))
{
grps->emplace_back(gid);
}
else
{
/* Store the group number in buffer */
- if (grptp == egrptpONE)
+ if (coverage == GroupCoverage::OneGroup)
{
cbuf[aj] = 0;
}
/* Now check whether we have done all atoms */
if (ntot != natoms)
{
- if (grptp == egrptpALL)
+ if (coverage == GroupCoverage::All)
{
gmx_fatal(FARGS, "%d atoms are not part of any of the %s groups", natoms - ntot, title);
}
- else if (grptp == egrptpPART)
+ else if (coverage == GroupCoverage::Partial)
{
sprintf(warn_buf, "%d atoms are not part of any of the %s groups", natoms - ntot, title);
warning_note(wi, warn_buf);
cbuf[j] = grps->size();
}
}
- if (grptp != egrptpPART)
+ if (coverage != GroupCoverage::Partial)
{
if (bVerbose)
{
const t_atom& local = atomP.atom();
int i = atomP.globalAtomNumber();
nrdf2[i] = 0;
- if (local.ptype == eptAtom || local.ptype == eptNucleus)
+ if (local.ptype == ParticleType::Atom || local.ptype == ParticleType::Nucleus)
{
int g = getGroupType(groups, SimulationAtomGroupType::Freeze, i);
for (int d = 0; d < DIM; d++)
*/
int ai = as + ia[i + 1];
int aj = as + ia[i + 2];
- if (((atom[ia[i + 1]].ptype == eptNucleus) || (atom[ia[i + 1]].ptype == eptAtom))
- && ((atom[ia[i + 2]].ptype == eptNucleus) || (atom[ia[i + 2]].ptype == eptAtom)))
+ if (((atom[ia[i + 1]].ptype == ParticleType::Nucleus)
+ || (atom[ia[i + 1]].ptype == ParticleType::Atom))
+ && ((atom[ia[i + 2]].ptype == ParticleType::Nucleus)
+ || (atom[ia[i + 2]].ptype == ParticleType::Atom)))
{
if (nrdf2[ai] > 0)
{
* But since this is much larger than STRLEN, such a line can not be parsed.
* The real maximum is the number of names that fit in a string: STRLEN/2.
*/
-#define EGP_MAX (STRLEN / 2)
int j, k, nr;
bool bSet;
j = 0;
while ((j < nr)
&& gmx_strcasecmp(
- names[2 * i].c_str(),
- *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][j]])))
+ names[2 * i].c_str(),
+ *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][j]])))
{
j++;
}
k = 0;
while ((k < nr)
&& gmx_strcasecmp(
- names[2 * i + 1].c_str(),
- *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][k]])))
+ names[2 * i + 1].c_str(),
+ *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][k]])))
{
k++;
}
}
}
-void do_index(const char* mdparin,
- const char* ndx,
- gmx_mtop_t* mtop,
- bool bVerbose,
- const gmx::MdModulesNotifier& notifier,
- t_inputrec* ir,
- warninp_t wi)
+void do_index(const char* mdparin,
+ const char* ndx,
+ gmx_mtop_t* mtop,
+ bool bVerbose,
+ const gmx::MDModulesNotifiers& mdModulesNotifiers,
+ t_inputrec* ir,
+ warninp_t wi)
{
t_blocka* defaultIndexGroups;
int natoms;
snew(defaultIndexGroups, 1);
snew(defaultIndexGroups->index, 1);
snew(gnames, 1);
- atoms_all = gmx_mtop_global_atoms(mtop);
+ atoms_all = gmx_mtop_global_atoms(*mtop);
analyse(&atoms_all, defaultIndexGroups, &gnames, FALSE, TRUE);
done_atom(&atoms_all);
}
gnames,
SimulationAtomGroupType::TemperatureCoupling,
restnm,
- useReferenceTemperature ? egrptpALL : egrptpALL_GENREST,
+ useReferenceTemperature ? GroupCoverage::All : GroupCoverage::AllGenerateRest,
bVerbose,
wi);
nr = groups->groups[SimulationAtomGroupType::TemperatureCoupling].size();
gmx::IndexGroupsAndNames defaultIndexGroupsAndNames(
*defaultIndexGroups, gmx::arrayRefFromArray(gnames, defaultIndexGroups->nr));
- notifier.preProcessingNotifications_.notify(defaultIndexGroupsAndNames);
+ mdModulesNotifiers.preProcessingNotifier_.notify(defaultIndexGroupsAndNames);
+
+ auto accelerations = gmx::splitString(inputrecStrings->acceleration);
+ auto accelerationGroupNames = gmx::splitString(inputrecStrings->accelerationGroups);
+ if (accelerationGroupNames.size() * DIM != accelerations.size())
+ {
+ gmx_fatal(FARGS,
+ "Invalid Acceleration input: %zu groups and %zu acc. values",
+ accelerationGroupNames.size(),
+ accelerations.size());
+ }
+ do_numbering(natoms,
+ groups,
+ accelerationGroupNames,
+ defaultIndexGroups,
+ gnames,
+ SimulationAtomGroupType::Acceleration,
+ restnm,
+ GroupCoverage::AllGenerateRest,
+ bVerbose,
+ wi);
+ nr = groups->groups[SimulationAtomGroupType::Acceleration].size();
+ snew(ir->opts.acceleration, nr);
+ ir->opts.ngacc = nr;
+
+ convertRvecs(wi, accelerations, "accelerations", ir->opts.acceleration);
auto freezeDims = gmx::splitString(inputrecStrings->frdim);
auto freezeGroupNames = gmx::splitString(inputrecStrings->freeze);
gnames,
SimulationAtomGroupType::Freeze,
restnm,
- egrptpALL_GENREST,
+ GroupCoverage::AllGenerateRest,
bVerbose,
wi);
nr = groups->groups[SimulationAtomGroupType::Freeze].size();
gnames,
SimulationAtomGroupType::EnergyOutput,
restnm,
- egrptpALL_GENREST,
+ GroupCoverage::AllGenerateRest,
bVerbose,
wi);
add_wall_energrps(groups, ir->nwall, symtab);
gnames,
SimulationAtomGroupType::MassCenterVelocityRemoval,
restnm,
- vcmGroupNames.empty() ? egrptpALL_GENREST : egrptpPART,
+ vcmGroupNames.empty() ? GroupCoverage::AllGenerateRest : GroupCoverage::Partial,
bVerbose,
wi);
gnames,
SimulationAtomGroupType::User1,
restnm,
- egrptpALL_GENREST,
+ GroupCoverage::AllGenerateRest,
bVerbose,
wi);
auto user2GroupNames = gmx::splitString(inputrecStrings->user2);
gnames,
SimulationAtomGroupType::User2,
restnm,
- egrptpALL_GENREST,
+ GroupCoverage::AllGenerateRest,
bVerbose,
wi);
auto compressedXGroupNames = gmx::splitString(inputrecStrings->x_compressed_groups);
gnames,
SimulationAtomGroupType::CompressedPositionOutput,
restnm,
- egrptpONE,
+ GroupCoverage::OneGroup,
bVerbose,
wi);
auto orirefFitGroupNames = gmx::splitString(inputrecStrings->orirefitgrp);
gnames,
SimulationAtomGroupType::OrientationRestraintsFit,
restnm,
- egrptpALL_GENREST,
+ GroupCoverage::AllGenerateRest,
bVerbose,
wi);
gnames,
SimulationAtomGroupType::QuantumMechanics,
restnm,
- egrptpALL_GENREST,
+ GroupCoverage::AllGenerateRest,
bVerbose,
wi);
ir->opts.ngQM = qmGroupNames.size();
}
-static void check_disre(const gmx_mtop_t* mtop)
+static void check_disre(const gmx_mtop_t& mtop)
{
if (gmx_mtop_ftype_count(mtop, F_DISRES) > 0)
{
- const gmx_ffparams_t& ffparams = mtop->ffparams;
+ const gmx_ffparams_t& ffparams = mtop.ffparams;
int ndouble = 0;
int old_label = -1;
for (int i = 0; i < ffparams.numTypes(); i++)
}
}
-static bool absolute_reference(const t_inputrec* ir, const gmx_mtop_t* sys, const bool posres_only, ivec AbsRef)
+//! Returns whether dimensions have an absolute reference due to walls, pbc or freezing
+static BasicVector<bool> haveAbsoluteReference(const t_inputrec& ir)
{
- int d, g, i;
- gmx_mtop_ilistloop_t iloop;
- int nmol;
- const t_iparams* pr;
+ BasicVector<bool> absRef = { false, false, false };
- clear_ivec(AbsRef);
-
- if (!posres_only)
+ /* Check the degrees of freedom of the COM (not taking COMM removal into account) */
+ for (int d = 0; d < DIM; d++)
{
- /* Check the COM */
- for (d = 0; d < DIM; d++)
- {
- AbsRef[d] = (d < ndof_com(ir) ? 0 : 1);
- }
- /* Check for freeze groups */
- for (g = 0; g < ir->opts.ngfrz; g++)
+ absRef[d] = (d >= ndof_com(&ir));
+ }
+ /* Check for freeze groups */
+ for (int g = 0; g < ir.opts.ngfrz; g++)
+ {
+ for (int d = 0; d < DIM; d++)
{
- for (d = 0; d < DIM; d++)
+ if (ir.opts.nFreeze[g][d] != 0)
{
- if (ir->opts.nFreeze[g][d] != 0)
- {
- AbsRef[d] = 1;
- }
+ absRef[d] = true;
}
}
}
- /* Check for position restraints */
- iloop = gmx_mtop_ilistloop_init(sys);
- while (const InteractionLists* ilist = gmx_mtop_ilistloop_next(iloop, &nmol))
+ return absRef;
+}
+
+//! Returns whether position restraints are used for dimensions
+static BasicVector<bool> havePositionRestraints(const gmx_mtop_t& sys)
+{
+ BasicVector<bool> havePosres = { false, false, false };
+
+ for (const auto ilists : IListRange(sys))
{
- if (nmol > 0 && (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
+ const auto& posResList = ilists.list()[F_POSRES];
+ const auto& fbPosResList = ilists.list()[F_FBPOSRES];
+ if (ilists.nmol() > 0 && (!havePosres[XX] || !havePosres[YY] || !havePosres[ZZ]))
{
- for (i = 0; i < (*ilist)[F_POSRES].size(); i += 2)
+ for (int i = 0; i < posResList.size(); i += 2)
{
- pr = &sys->ffparams.iparams[(*ilist)[F_POSRES].iatoms[i]];
- for (d = 0; d < DIM; d++)
+ const t_iparams& pr = sys.ffparams.iparams[posResList.iatoms[i]];
+ for (int d = 0; d < DIM; d++)
{
- if (pr->posres.fcA[d] != 0)
+ if (pr.posres.fcA[d] != 0)
{
- AbsRef[d] = 1;
+ havePosres[d] = true;
}
}
}
- for (i = 0; i < (*ilist)[F_FBPOSRES].size(); i += 2)
+ for (int i = 0; i < fbPosResList.size(); i += 2)
{
/* Check for flat-bottom posres */
- pr = &sys->ffparams.iparams[(*ilist)[F_FBPOSRES].iatoms[i]];
- if (pr->fbposres.k != 0)
+ const t_iparams& pr = sys.ffparams.iparams[fbPosResList.iatoms[i]];
+ if (pr.fbposres.k != 0)
{
- switch (pr->fbposres.geom)
+ switch (pr.fbposres.geom)
{
- case efbposresSPHERE: AbsRef[XX] = AbsRef[YY] = AbsRef[ZZ] = 1; break;
- case efbposresCYLINDERX: AbsRef[YY] = AbsRef[ZZ] = 1; break;
- case efbposresCYLINDERY: AbsRef[XX] = AbsRef[ZZ] = 1; break;
+ case efbposresSPHERE: havePosres = { true, true, true }; break;
+ case efbposresCYLINDERX: havePosres[YY] = havePosres[ZZ] = true; break;
+ case efbposresCYLINDERY: havePosres[XX] = havePosres[ZZ] = true; break;
case efbposresCYLINDER:
/* efbposres is a synonym for efbposresCYLINDERZ for backwards compatibility */
- case efbposresCYLINDERZ: AbsRef[XX] = AbsRef[YY] = 1; break;
+ case efbposresCYLINDERZ: havePosres[XX] = havePosres[YY] = true; break;
case efbposresX: /* d=XX */
case efbposresY: /* d=YY */
case efbposresZ: /* d=ZZ */
- d = pr->fbposres.geom - efbposresX;
- AbsRef[d] = 1;
+ havePosres[pr.fbposres.geom - efbposresX] = true;
break;
default:
gmx_fatal(FARGS,
- " Invalid geometry for flat-bottom position restraint.\n"
+ "Invalid geometry for flat-bottom position restraint.\n"
"Expected nr between 1 and %d. Found %d\n",
efbposresNR - 1,
- pr->fbposres.geom);
+ pr.fbposres.geom);
}
}
}
}
}
- return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
+ return havePosres;
}
-static void check_combination_rule_differences(const gmx_mtop_t* mtop,
+static void check_combination_rule_differences(const gmx_mtop_t& mtop,
int state,
bool* bC6ParametersWorkWithGeometricRules,
bool* bC6ParametersWorkWithLBRules,
ptr = getenv("GMX_LJCOMB_TOL");
if (ptr != nullptr)
{
- double dbl;
+ double dbl;
double gmx_unused canary;
if (sscanf(ptr, "%lf%lf", &dbl, &canary) != 1)
*bC6ParametersWorkWithLBRules = TRUE;
*bC6ParametersWorkWithGeometricRules = TRUE;
bCanDoLBRules = TRUE;
- ntypes = mtop->ffparams.atnr;
+ ntypes = mtop.ffparams.atnr;
snew(typecount, ntypes);
gmx_mtop_count_atomtypes(mtop, state, typecount);
*bLBRulesPossible = TRUE;
for (tpi = 0; tpi < ntypes; ++tpi)
{
- c6i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c6;
- c12i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c12;
+ c6i = mtop.ffparams.iparams[(ntypes + 1) * tpi].lj.c6;
+ c12i = mtop.ffparams.iparams[(ntypes + 1) * tpi].lj.c12;
for (tpj = tpi; tpj < ntypes; ++tpj)
{
- c6j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c6;
- c12j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c12;
- c6 = mtop->ffparams.iparams[ntypes * tpi + tpj].lj.c6;
+ c6j = mtop.ffparams.iparams[(ntypes + 1) * tpj].lj.c6;
+ c12j = mtop.ffparams.iparams[(ntypes + 1) * tpj].lj.c12;
+ c6 = mtop.ffparams.iparams[ntypes * tpi + tpj].lj.c6;
c6_geometric = std::sqrt(c6i * c6j);
if (!gmx_numzero(c6_geometric))
{
sfree(typecount);
}
-static void check_combination_rules(const t_inputrec* ir, const gmx_mtop_t* mtop, warninp_t wi)
+static void check_combination_rules(const t_inputrec* ir, const gmx_mtop_t& mtop, warninp_t wi)
{
bool bLBRulesPossible, bC6ParametersWorkWithGeometricRules, bC6ParametersWorkWithLBRules;
}
}
+static bool allTrue(const BasicVector<bool>& boolVector)
+{
+ return boolVector[0] && boolVector[1] && boolVector[2];
+}
+
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
char err_buf[STRLEN];
int i, m, c, nmol;
bool bCharge;
+ real * mgrp, mt;
gmx_mtop_atomloop_block_t aloopb;
- ivec AbsRef;
char warn_buf[STRLEN];
set_warning_line(wi, mdparin, -1);
- if (absolute_reference(ir, sys, false, AbsRef))
+ if (ir->comm_mode != ComRemovalAlgorithm::No && allTrue(havePositionRestraints(*sys)))
{
warning_note(wi,
"Removing center of mass motion in the presence of position restraints might "
* of errors. The factor 0.5 is because energy distributes
* equally over Ekin and Epot.
*/
- max_T_error = 0.5 * tau * ir->verletbuf_tol / (nrdf_at * BOLTZ * T);
+ max_T_error = 0.5 * tau * ir->verletbuf_tol / (nrdf_at * gmx::c_boltz * T);
if (max_T_error > T_error_warn)
{
sprintf(warn_buf,
if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != IntegrationAlgorithm::BD
&& ir->comm_mode == ComRemovalAlgorithm::No
- && !(absolute_reference(ir, sys, FALSE, AbsRef) || ir->nsteps <= 10) && !ETC_ANDERSEN(ir->etc))
+ && !(allTrue(haveAbsoluteReference(*ir)) || allTrue(havePositionRestraints(*sys)) || ir->nsteps <= 10)
+ && !ETC_ANDERSEN(ir->etc))
{
warning(wi,
"You are not using center of mass motion removal (mdp option comm-mode), numerical "
/* Check for pressure coupling with absolute position restraints */
if (ir->epc != PressureCoupling::No && ir->refcoord_scaling == RefCoordScaling::No)
{
- absolute_reference(ir, sys, TRUE, AbsRef);
+ const BasicVector<bool> havePosres = havePositionRestraints(*sys);
{
for (m = 0; m < DIM; m++)
{
- if (AbsRef[m] && norm2(ir->compress[m]) > 0)
+ if (havePosres[m] && norm2(ir->compress[m]) > 0)
{
warning(wi,
"You are using pressure coupling with absolute position restraints, "
}
bCharge = FALSE;
- aloopb = gmx_mtop_atomloop_block_init(sys);
+ aloopb = gmx_mtop_atomloop_block_init(*sys);
const t_atom* atom;
while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
{
/* Check if combination rules used in LJ-PME are the same as in the force field */
if (EVDW_PME(ir->vdwtype))
{
- check_combination_rules(ir, sys, wi);
+ check_combination_rules(ir, *sys, wi);
}
/* Generalized reaction field */
"constant by hand.");
}
+ ir->useConstantAcceleration = false;
+ for (int i = 0; (i < gmx::ssize(sys->groups.groups[SimulationAtomGroupType::Acceleration])); i++)
+ {
+ if (norm2(ir->opts.acceleration[i]) != 0)
+ {
+ ir->useConstantAcceleration = true;
+ }
+ }
+ if (ir->useConstantAcceleration)
+ {
+ gmx::RVec acceleration = { 0.0_real, 0.0_real, 0.0_real };
+ snew(mgrp, sys->groups.groups[SimulationAtomGroupType::Acceleration].size());
+ for (const AtomProxy atomP : AtomRange(*sys))
+ {
+ const t_atom& local = atomP.atom();
+ int i = atomP.globalAtomNumber();
+ mgrp[getGroupType(sys->groups, SimulationAtomGroupType::Acceleration, i)] += local.m;
+ }
+ mt = 0.0;
+ for (i = 0; (i < gmx::ssize(sys->groups.groups[SimulationAtomGroupType::Acceleration])); i++)
+ {
+ for (m = 0; (m < DIM); m++)
+ {
+ acceleration[m] += ir->opts.acceleration[i][m] * mgrp[i];
+ }
+ mt += mgrp[i];
+ }
+ for (m = 0; (m < DIM); m++)
+ {
+ if (fabs(acceleration[m]) > 1e-6)
+ {
+ const char* dim[DIM] = { "X", "Y", "Z" };
+ fprintf(stderr,
+ "Net Acceleration in %s direction, will %s be corrected\n",
+ dim[m],
+ ir->nstcomm != 0 ? "" : "not");
+ if (ir->nstcomm != 0 && m < ndof_com(ir))
+ {
+ acceleration[m] /= mt;
+ for (i = 0;
+ (i < gmx::ssize(sys->groups.groups[SimulationAtomGroupType::Acceleration]));
+ i++)
+ {
+ ir->opts.acceleration[i][m] -= acceleration[m];
+ }
+ }
+ }
+ }
+ sfree(mgrp);
+ }
+
if (ir->efep != FreeEnergyPerturbationType::No && ir->fepvals->sc_alpha != 0
&& !gmx_within_tol(sys->ffparams.reppow, 12.0, 10 * GMX_DOUBLE_EPS))
{
bWarned = FALSE;
for (i = 0; i < ir->pull->ncoord && !bWarned; i++)
{
- if (ir->pull->coord[i].group[0] == 0 || ir->pull->coord[i].group[1] == 0)
+ if (ir->pull->coord[i].eGeom != PullGroupGeometry::Transformation
+ && (ir->pull->coord[i].group[0] == 0 || ir->pull->coord[i].group[1] == 0))
{
- absolute_reference(ir, sys, FALSE, AbsRef);
+ const auto absRef = haveAbsoluteReference(*ir);
+ const auto havePosres = havePositionRestraints(*sys);
for (m = 0; m < DIM; m++)
{
- if (ir->pull->coord[i].dim[m] && !AbsRef[m])
+ if (ir->pull->coord[i].dim[m] && !(absRef[m] || havePosres[m]))
{
warning(wi,
"You are using an absolute reference for pulling, but the rest of "
}
}
- check_disre(sys);
+ check_disre(*sys);
}
void double_check(t_inputrec* ir, matrix box, bool bHasNormalConstraints, bool bHasAnyConstraints, warninp_t wi)
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff