#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/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
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)
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)
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");
ir->userreal4 = get_ereal(&inp, "userreal4", 0, wi);
#undef CTYPE
+ if (mdparout)
{
gmx::TextOutputFile stream(mdparout);
write_inpfile(&stream, mdparout, &inp, FALSE, writeMdpHeader, wi);
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)
{
* 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++;
}
gnames,
SimulationAtomGroupType::TemperatureCoupling,
restnm,
- useReferenceTemperature ? egrptpALL : egrptpALL_GENREST,
+ useReferenceTemperature ? GroupCoverage::All : GroupCoverage::AllGenerateRest,
bVerbose,
wi);
nr = groups->groups[SimulationAtomGroupType::TemperatureCoupling].size();
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 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;
- 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 */
- for (const auto ilist : IListRange(sys))
+ 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 (ilist.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 (int i = 0; i < ilist.list()[F_POSRES].size(); i += 2)
+ for (int i = 0; i < posResList.size(); i += 2)
{
- pr = &sys.ffparams.iparams[ilist.list()[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.list()[F_FBPOSRES].size(); i += 2)
+ for (int i = 0; i < fbPosResList.size(); i += 2)
{
/* Check for flat-bottom posres */
- pr = &sys.ffparams.iparams[ilist.list()[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,
ptr = getenv("GMX_LJCOMB_TOL");
if (ptr != nullptr)
{
- double dbl;
+ double dbl;
double gmx_unused canary;
if (sscanf(ptr, "%lf%lf", &dbl, &canary) != 1)
}
}
+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
int i, m, c, nmol;
bool bCharge;
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 "
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, "
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 "
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff