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;
}
-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)
if (fep->softcoreFunction == SoftcoreType::Gapsys)
{
- if (fep->sc_alpha < 0.0)
+ if (fep->scGapsysScaleLinpointQ < 0.0)
{
sprintf(warn_buf,
- "sc_alpha is equal %g but must be >= 0 when used with sc_function=gapsys.",
- fep->sc_alpha);
+ "sc_scale_linpoint_Q_gapsys is equal %g but must be >= 0",
+ fep->scGapsysScaleLinpointQ);
warning_note(wi, warn_buf);
}
- if ((fep->sc_sigma < 0.0) || (fep->sc_sigma >= 1.0))
+ if ((fep->scGapsysScaleLinpointLJ < 0.0) || (fep->scGapsysScaleLinpointLJ >= 1.0))
{
sprintf(warn_buf,
- "sc_sigma is equal %g but must be in [0,1) when used with "
+ "sc_scale_linpoint_LJ_gapsys is equal %g but must be in [0,1) when used "
+ "with "
"sc_function=gapsys.",
- fep->sc_sigma);
+ fep->scGapsysScaleLinpointLJ);
warning_note(wi, warn_buf);
}
}
}
}
+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;
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->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->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);
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)
{
gnames,
SimulationAtomGroupType::TemperatureCoupling,
restnm,
- useReferenceTemperature ? egrptpALL : egrptpALL_GENREST,
+ useReferenceTemperature ? GroupCoverage::All : GroupCoverage::AllGenerateRest,
bVerbose,
wi);
nr = groups->groups[SimulationAtomGroupType::TemperatureCoupling].size();
*defaultIndexGroups, gmx::arrayRefFromArray(gnames, defaultIndexGroups->nr));
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);
if (freezeDims.size() != DIM * freezeGroupNames.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();
char err_buf[STRLEN];
int i, m, c, nmol;
bool bCharge;
+ real * mgrp, mt;
gmx_mtop_atomloop_block_t aloopb;
char warn_buf[STRLEN];
"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))
{
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff