used as a reference atom, but the user can also select any other atom if
it would be closer to center of the group.
+When there are large pull groups, such as a
+lipid bilayer, ``pull-pbc-ref-prev-step-com`` can be used to avoid potential
+large movements of the center of mass in case that atoms in the pull group
+move so much that the reference atom is too far from the intended center of mass.
+With this option enabled the center of mass from the previous step is used,
+instead of the position of the reference atom, to determine the reference position.
+The position of the reference atom is still used for the first step. For large pull
+groups it is important to select a reference atom that is close to the intended
+center of mass, i.e. do not use ``pull-group?-pbcatom = 0``.
+
For a layered system, for instance a lipid bilayer, it may be of
interest to calculate the PMF of a lipid as function of its distance
from the whole bilayer. The whole bilayer can be taken as reference
to the GPU the same way as offload is done to other GPU architectures.
This can have performance benefits, in particular on modern desktop and mobile
Intel CPUs this offload can give up to 20% higher simulation performance.
+
+Allow using COM of previous step as PBC reference
+"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+Added an option (``pull-pbc-ref-from-prev-step-com``), when pulling, to use
+the COM of the group of the previous step, to calculate PBC jumps instead of a
+reference atom, which can sometimes move a lot during the simulation.
+With this option the PBC reference atom is only used at initialization.
+This can be of use when using large pull groups or groups with potentially
+large relative movement of atoms.
frequency for writing out the force of all the pulled group
(0 is never)
+.. mdp:: pull-pbc-ref-prev-step-com
+
+ .. mdp-value:: no
+
+ Use the reference atom (:mdp:`pull-group1-pbcatom`) for the
+ treatment of periodic boundary conditions.
+
+ .. mdp-value:: yes
+
+ Use the COM of the previous step as reference for the treatment
+ of periodic boundary conditions. The reference is initialized
+ using the reference atom (:mdp:`pull-group1-pbcatom`), which should
+ be located centrally in the group. Using the COM from the
+ previous step can be useful if one or more pull groups are large.
.. mdp:: pull-ngroups
vector. For determining the COM, all atoms in the group are put at
their periodic image which is closest to
:mdp:`pull-group1-pbcatom`. A value of 0 means that the middle
- atom (number wise) is used. This parameter is not used with
+ atom (number wise) is used, which is only safe for small groups.
+ :ref:`gmx grompp` checks that the maximum distance from the reference
+ atom (specifically chosen, or not) to the other atoms in the group
+ is not too large. This parameter is not used with
:mdp:`pull-coord1-geometry` cylinder. A value of -1 turns on cosine
weighting, which is useful for a group of molecules in a periodic
system, *e.g.* a water slab (see Engin et al. J. Chem. Phys. B
enum cptv {
cptv_Unknown = 17, /**< Version before numbering scheme */
cptv_RemoveBuildMachineInformation, /**< remove functionality that makes mdrun builds non-reproducible */
+ cptv_ComPrevStepAsPullGroupReference, /**< Allow using COM of previous step as pull group PBC reference */
cptv_Count /**< the total number of cptv versions */
};
"awh_forceCorrelationGrid"
};
+enum {
+ epullsPREVSTEPCOM,
+ epullsNR
+};
+
+static const char *epull_prev_step_com_names[epullsNR] =
+{
+ "Pull groups prev step COM"
+};
+
+
//! Higher level vector element type, only used for formatting checkpoint dumps
enum class CptElementType
{
kineticEnergy, //!< Kinetic energy, needed for T/P-coupling state
energyHistory, //!< Energy observable statistics
freeEnergyHistory, //!< Free-energy state and observable statistics
- accWeightHistogram //!< Accelerated weight histogram method state
+ accWeightHistogram, //!< Accelerated weight histogram method state
+ pullState //!< COM of previous step.
};
//! \brief Return the name of a checkpoint entry based on part and part entry
case StatePart::energyHistory: return eenh_names[ecpt];
case StatePart::freeEnergyHistory: return edfh_names[ecpt];
case StatePart::accWeightHistogram: return eawhh_names[ecpt];
+ case StatePart::pullState: return epull_prev_step_com_names[ecpt];
}
return nullptr;
case estDISRE_RM3TAV: ret = do_cpte_n_reals(xd, part, i, sflags, &state->hist.ndisrepairs, &state->hist.disre_rm3tav, list); break;
case estORIRE_INITF: ret = do_cpte_real (xd, part, i, sflags, &state->hist.orire_initf, list); break;
case estORIRE_DTAV: ret = do_cpte_n_reals(xd, part, i, sflags, &state->hist.norire_Dtav, &state->hist.orire_Dtav, list); break;
+ case estPREVSTEPCOM: ret = doVector<double>(xd, part, i, sflags, &state->com_prev_step, list); break;
default:
gmx_fatal(FARGS, "Unknown state entry %d\n"
"You are reading a checkpoint file written by different code, which is not supported", i);
}
}
}
-
return ret;
}
tpxv_GenericParamsForElectricField, /**< Introduced KeyValueTree and moved electric field parameters */
tpxv_AcceleratedWeightHistogram, /**< sampling with accelerated weight histogram method (AWH) */
tpxv_RemoveImplicitSolvation, /**< removed support for implicit solvation */
+ tpxv_PullPrevStepCOMAsReference, /**< Enabled using the COM of the pull group of the last frame as reference for PBC */
tpxv_Count /**< the total number of tpxv versions */
};
}
gmx_fio_do_int(fio, pull->nstxout);
gmx_fio_do_int(fio, pull->nstfout);
+ if (file_version >= tpxv_PullPrevStepCOMAsReference)
+ {
+ gmx_fio_do_gmx_bool(fio, pull->bSetPbcRefToPrevStepCOM);
+ }
+ else
+ {
+ pull->bSetPbcRefToPrevStepCOM = FALSE;
+ }
if (bRead)
{
snew(pull->group, pull->ngroup);
/* read pull parameters */
printStringNoNewline(inp, "Cylinder radius for dynamic reaction force groups (nm)");
- pull->cylinder_r = get_ereal(inp, "pull-cylinder-r", 1.5, wi);
- pull->constr_tol = get_ereal(inp, "pull-constr-tol", 1E-6, wi);
- pull->bPrintCOM = (get_eeenum(inp, "pull-print-com", yesno_names, wi) != 0);
- pull->bPrintRefValue = (get_eeenum(inp, "pull-print-ref-value", yesno_names, wi) != 0);
- pull->bPrintComp = (get_eeenum(inp, "pull-print-components", yesno_names, wi) != 0);
- pull->nstxout = get_eint(inp, "pull-nstxout", 50, wi);
- pull->nstfout = get_eint(inp, "pull-nstfout", 50, wi);
+ pull->cylinder_r = get_ereal(inp, "pull-cylinder-r", 1.5, wi);
+ pull->constr_tol = get_ereal(inp, "pull-constr-tol", 1E-6, wi);
+ pull->bPrintCOM = (get_eeenum(inp, "pull-print-com", yesno_names, wi) != 0);
+ pull->bPrintRefValue = (get_eeenum(inp, "pull-print-ref-value", yesno_names, wi) != 0);
+ pull->bPrintComp = (get_eeenum(inp, "pull-print-components", yesno_names, wi) != 0);
+ pull->nstxout = get_eint(inp, "pull-nstxout", 50, wi);
+ pull->nstfout = get_eint(inp, "pull-nstfout", 50, wi);
+ pull->bSetPbcRefToPrevStepCOM = (get_eeenum(inp, "pull-pbc-ref-prev-step-com", yesno_names, wi) != 0);
printStringNoNewline(inp, "Number of pull groups");
pull->ngroup = get_eint(inp, "pull-ngroups", 1, wi);
printStringNoNewline(inp, "Number of pull coordinates");
t_pull_group *pgrp;
/* Absolute reference group (might not be used) is special */
- pgrp = &pull->group[0];
- pgrp->nat = 0;
- pgrp->pbcatom = -1;
+ pgrp = &pull->group[0];
+ pgrp->nat = 0;
+ pgrp->pbcatom = -1;
+ pgrp->pbcatom_input = -1;
for (g = 1; g < pull->ngroup; g++)
{
- pgrp = &pull->group[g];
+ pgrp = &pull->group[g];
+ pgrp->pbcatom_input = pgrp->pbcatom;
if (strcmp(pgnames[g], "") == 0)
{
t_start = ir->init_t + ir->init_step*ir->delta_t;
+ if (pull->bSetPbcRefToPrevStepCOM)
+ {
+ initPullComFromPrevStep(nullptr, pull_work, md, &pbc, x);
+ }
pull_calc_coms(nullptr, pull_work, md, &pbc, t_start, x, nullptr);
- int groupThatFailsPbc = pullCheckPbcWithinGroups(*pull_work, x, pbc, c_pullGroupPbcMargin);
- if (groupThatFailsPbc >= 0)
+ for (int g = 0; g < pull->ngroup; g++)
{
- char buf[STRLEN];
- sprintf(buf,
- "Pull group %d has atoms at a distance larger than %g times half the box size from the PBC atom (%d). If atoms are or will more beyond half the box size from the PBC atom, the COM will be ill defined.",
- groupThatFailsPbc,
- c_pullGroupPbcMargin,
- pull->group[groupThatFailsPbc].pbcatom);
- set_warning_line(wi, nullptr, -1);
- warning(wi, buf);
+ bool groupObeysPbc =
+ pullCheckPbcWithinGroup(*pull_work,
+ gmx::arrayRefFromArray(reinterpret_cast<gmx::RVec *>(x),
+ mtop->natoms),
+ pbc, g, c_pullGroupSmallGroupThreshold);
+ if (!groupObeysPbc)
+ {
+ char buf[STRLEN];
+ if (pull->group[g].pbcatom_input == 0)
+ {
+ sprintf(buf, "When the maximum distance from a pull group reference atom to other atoms in the "
+ "group is larger than %g times half the box size a centrally placed "
+ "atom should be chosen as pbcatom. Pull group %d is larger than that and does not have "
+ "a specific atom selected as reference atom.", c_pullGroupSmallGroupThreshold, g);
+ warning_error(wi, buf);
+ }
+ else if (!pull->bSetPbcRefToPrevStepCOM)
+ {
+ sprintf(buf, "The maximum distance from the chosen PBC atom (%d) of pull group %d to other "
+ "atoms in the group is larger than %g times half the box size. "
+ "Set the pull-pbc-ref-prev-step-com option to yes.", pull->group[g].pbcatom + 1,
+ g, c_pullGroupSmallGroupThreshold);
+ warning_error(wi, buf);
+ }
+ }
+ if (groupObeysPbc)
+ {
+ groupObeysPbc =
+ pullCheckPbcWithinGroup(*pull_work,
+ gmx::arrayRefFromArray(reinterpret_cast<gmx::RVec *>(x),
+ mtop->natoms),
+ pbc, g, c_pullGroupPbcMargin);
+ if (!groupObeysPbc)
+ {
+ char buf[STRLEN];
+ sprintf(buf,
+ "Pull group %d has atoms at a distance larger than %g times half the box size from the PBC atom (%d). "
+ "If atoms are or will more beyond half the box size from the PBC atom, the COM will be ill defined.",
+ g, c_pullGroupPbcMargin, pull->group[g].pbcatom + 1);
+ set_warning_line(wi, nullptr, -1);
+ warning(wi, buf);
+ }
+ }
}
fprintf(stderr, "Pull group natoms pbc atom distance at start reference at t=0\n");
#include "gromacs/mdtypes/md_enums.h"
#include "gromacs/mdtypes/state.h"
#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/pulling/pull.h"
#include "gromacs/timing/wallcycle.h"
#include "gromacs/topology/mtop_util.h"
#include "gromacs/trajectory/trajectoryframe.h"
snew(state->dfhist, 1);
init_df_history(state->dfhist, ir->fepvals->n_lambda);
}
+
+ if (ir->pull && ir->pull->bSetPbcRefToPrevStepCOM)
+ {
+ state->flags |= (1<<estPREVSTEPCOM);
+ }
}
*/
observablesHistory->energyHistory = {};
}
+ if (ir->pull && ir->pull->bSetPbcRefToPrevStepCOM)
+ {
+ /* Copy the pull group COM of the previous step from the checkpoint state to the pull state */
+ setPrevStepPullComFromState(ir->pull_work, state);
+ }
+ }
+ else if (ir->pull && ir->pull->bSetPbcRefToPrevStepCOM)
+ {
+ allocStatePrevStepPullCom(state, ir->pull_work);
+ t_pbc pbc;
+ set_pbc(&pbc, ir->ePBC, state->box);
+ initPullComFromPrevStep(cr, ir->pull_work, mdatoms, &pbc, as_rvec_array(state->x.data()));
+ updatePrevStepCom(ir->pull_work);
+ setStatePrevStepPullCom(ir->pull_work, state);
}
if (observablesHistory->energyHistory == nullptr)
{
state, graph,
nrnb, wcycle, upd, constr);
+ if (MASTER(cr) && ir->bPull && ir->pull->bSetPbcRefToPrevStepCOM)
+ {
+ updatePrevStepCom(ir->pull_work);
+ setStatePrevStepPullCom(ir->pull_work, state);
+ }
+
if (ir->eI == eiVVAK)
{
/* erase F_EKIN and F_TEMP here? */
PS("pull-print-components", EBOOL(pull->bPrintComp));
PI("pull-nstxout", pull->nstxout);
PI("pull-nstfout", pull->nstfout);
+ PS("pull-pbc-ref-prev-step-com", EBOOL(pull->bSetPbcRefToPrevStepCOM));
PI("pull-ngroups", pull->ngroup);
for (g = 0; g < pull->ngroup; g++)
{
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2015,2016, by the GROMACS development team, led by
+ * Copyright (c) 2015,2016,2018, 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.
/*! \brief Struct that defines a pull group */
typedef struct {
- int nat; /**< Number of atoms in the pull group */
- int *ind; /**< The global atoms numbers */
- int nweight; /**< The number of weights (0 or nat) */
- real *weight; /**< Weights (use all 1 when weight==NULL) */
- int pbcatom; /**< The reference atom for pbc (global number) */
+ int nat; /**< Number of atoms in the pull group */
+ int *ind; /**< The global atoms numbers */
+ int nweight; /**< The number of weights (0 or nat) */
+ real *weight; /**< Weights (use all 1 when weight==NULL) */
+ int pbcatom; /**< The reference atom for pbc (global number) */
+ int pbcatom_input; /**< The reference atom for pbc (global number) as specified in the input parameters */
} t_pull_group;
/*! Maximum number of pull groups that can be used in a pull coordinate */
/*! \brief Struct containing all pull parameters */
typedef struct pull_params_t {
- int ngroup; /**< Number of pull groups */
- int ncoord; /**< Number of pull coordinates */
- real cylinder_r; /**< Radius of cylinder for dynamic COM (nm) */
- real constr_tol; /**< Absolute tolerance for constraints in (nm) */
- gmx_bool bPrintCOM; /**< Print coordinates of COM for each coord */
- gmx_bool bPrintRefValue; /**< Print the reference value for each coord */
- gmx_bool bPrintComp; /**< Print cartesian components for each coord with geometry=distance */
- int nstxout; /**< Output interval for pull x */
- int nstfout; /**< Output interval for pull f */
+ int ngroup; /**< Number of pull groups */
+ int ncoord; /**< Number of pull coordinates */
+ real cylinder_r; /**< Radius of cylinder for dynamic COM (nm) */
+ real constr_tol; /**< Absolute tolerance for constraints in (nm) */
+ gmx_bool bPrintCOM; /**< Print coordinates of COM for each coord */
+ gmx_bool bPrintRefValue; /**< Print the reference value for each coord */
+ gmx_bool bPrintComp; /**< Print cartesian components for each coord with geometry=distance */
+ gmx_bool bSetPbcRefToPrevStepCOM; /**< Use the COM of each group from the previous step as reference */
+ int nstxout; /**< Output interval for pull x */
+ int nstfout; /**< Output interval for pull f */
- t_pull_group *group; /**< groups to pull/restrain/etc/ */
- t_pull_coord *coord; /**< the pull coordinates */
+ t_pull_group *group; /**< groups to pull/restrain/etc/ */
+ t_pull_coord *coord; /**< the pull coordinates */
} pull_params_t;
/*! \endcond */
#include "gromacs/mdtypes/df_history.h"
#include "gromacs/mdtypes/inputrec.h"
#include "gromacs/mdtypes/md_enums.h"
+#include "gromacs/mdtypes/pull-params.h"
#include "gromacs/mdtypes/swaphistory.h"
#include "gromacs/pbcutil/boxutilities.h"
#include "gromacs/pbcutil/pbc.h"
estORIRE_INITF, estORIRE_DTAV,
estSVIR_PREV, estNH_VXI, estVETA, estVOL0, estNHPRES_XI, estNHPRES_VXI, estFVIR_PREV,
estFEPSTATE, estMC_RNG_NOTSUPPORTED, estMC_RNGI_NOTSUPPORTED,
- estBAROS_INT,
+ estBAROS_INT, estPREVSTEPCOM,
estNR
};
int ddp_count; //!< The DD partitioning count for this state
int ddp_count_cg_gl; //!< The DD partitioning count for index_gl
std::vector<int> cg_gl; //!< The global cg number of the local cgs
+
+ std::vector<double> com_prev_step; //!< The COM of the previous step of each pull group
};
#ifndef DOXYGEN
#include "pull_internal.h"
-static int groupPbcFromParams(const t_pull_group ¶ms)
+static int groupPbcFromParams(const t_pull_group ¶ms, bool setPbcRefToPrevStepCOM)
{
if (params.nat <= 1)
{
}
else if (params.pbcatom >= 0)
{
- return epgrppbcREFAT;
+ if (setPbcRefToPrevStepCOM)
+ {
+ return epgrppbcPREVSTEPCOM;
+ }
+ else
+ {
+ return epgrppbcREFAT;
+ }
}
else
{
* This will be fixed when the pointers are replacted by std::vector.
*/
pull_group_work_t::pull_group_work_t(const t_pull_group ¶ms,
- gmx::LocalAtomSet atomSet) :
+ gmx::LocalAtomSet atomSet,
+ bool bSetPbcRefToPrevStepCOM) :
params(params),
- epgrppbc(groupPbcFromParams(params)),
+ epgrppbc(groupPbcFromParams(params, bSetPbcRefToPrevStepCOM)),
needToCalcCom(false),
atomSet(atomSet),
mwscale(0),
/* We should participate if we have pull or pbc atoms */
if (!bMustParticipate &&
- (group.atomSet.numAtomsLocal() > 0 ||
- (group.epgrppbc == epgrppbcREFAT &&
- group.pbcAtomSet->numAtomsLocal() > 0)))
+ (group.atomSet.numAtomsLocal() > 0 || (group.epgrppbc == epgrppbcREFAT &&
+ group.pbcAtomSet->numAtomsLocal() > 0)))
{
bMustParticipate = TRUE;
}
for (int i = 0; i < pull_params->ngroup; ++i)
{
- pull->group.emplace_back(pull_params->group[i], atomSets->add({pull_params->group[i].ind, pull_params->group[i].ind+pull_params->group[i].nat}));
+ pull->group.emplace_back(pull_params->group[i], atomSets->add({pull_params->group[i].ind, pull_params->group[i].ind+pull_params->group[i].nat}),
+ pull_params->bSetPbcRefToPrevStepCOM);
}
if (cr != nullptr && DOMAINDECOMP(cr))
pull->bAngle = FALSE;
GMX_RELEASE_ASSERT(pull->group[0].params.nat == 0, "pull group 0 is an absolute reference group and should not contain atoms");
+ pull->group[0].x_prev_step[XX] = NAN;
pull->numCoordinatesWithExternalPotential = 0;
init_pull_group_index(fplog, cr, g, pgrp,
bConstraint, pulldim_con,
mtop, ir, lambda);
+
+ pgrp->x_prev_step[XX] = NAN;
}
else
{
}
}
const auto &referenceGroup = pull->group[coord.params.group[0]];
- pull->dyna.emplace_back(referenceGroup.params, referenceGroup.atomSet);
+ pull->dyna.emplace_back(referenceGroup.params, referenceGroup.atomSet, pull->params.bSetPbcRefToPrevStepCOM);
}
}
struct t_inputrec;
struct t_mdatoms;
struct t_pbc;
+class t_state;
namespace gmx
{
/*! \brief Margin for checking pull group PBC distances compared to half the box size */
static constexpr real c_pullGroupPbcMargin = 0.9;
+/*! \brief Threshold (as a factor of half the box size) for accepting pull groups without explicitly set refatom */
+static constexpr real c_pullGroupSmallGroupThreshold = 0.5;
/*! \brief Checks whether all groups that use a reference atom are within PBC restrictions
*
real max_pull_distance2(const pull_coord_work_t *pcrd,
const t_pbc *pbc);
+/*! \brief Copies the COM from the previous step of all pull groups to the checkpoint state container
+ *
+ * \param[in] pull The COM pull force calculation data structure
+ * \param[in] state The global state container
+ */
+void setStatePrevStepPullCom(const struct pull_t *pull, t_state *state);
+
+/*! \brief Copies the pull group COM of the previous step from the checkpoint state to the pull state
+ *
+ * \param[in] pull The COM pull force calculation data structure
+ * \param[in] state The global state container
+ */
+void setPrevStepPullComFromState(struct pull_t *pull, const t_state *state);
+
+/*! \brief Sets the previous step COM to the current COM
+ *
+ * \param[in] pull The COM pull force calculation data structure
+ */
+void updatePrevStepCom(struct pull_t *pull);
+
+/*! \brief Resizes the vector, in the state container, containing the COMs from the previous step
+ *
+ * \param[in] state The global state container
+ * \param[in] pull The COM pull force calculation data structure
+ */
+void allocStatePrevStepPullCom(t_state *state, pull_t *pull);
+
+/*! \brief Initializes the COM of the previous step (set to initial COM)
+ *
+ * \param[in] cr Struct for communication info.
+ * \param[in] pull The pull data structure.
+ * \param[in] md All atoms.
+ * \param[in] pbc Information struct about periodicity.
+ * \param[in] x The local positions.
+ */
+void initPullComFromPrevStep(const t_commrec *cr,
+ pull_t *pull,
+ const t_mdatoms *md,
+ t_pbc *pbc,
+ const rvec x[]);
+
#endif
#endif
enum {
- epgrppbcNONE, epgrppbcREFAT, epgrppbcCOS
+ epgrppbcNONE, epgrppbcREFAT, epgrppbcCOS, epgrppbcPREVSTEPCOM
};
/*! \internal
{
/*! \brief Constructor
*
- * \param[in] params The group parameters set by the user
- * \param[in] atomSet The global to local atom set manager
+ * \param[in] params The group parameters set by the user
+ * \param[in] atomSet The global to local atom set manager
+ * \param[in] setPbcRefToPrevStepCOM Does this pull group use the COM from the previous step as reference position?
*/
pull_group_work_t(const t_pull_group ¶ms,
- gmx::LocalAtomSet atomSet);
+ gmx::LocalAtomSet atomSet,
+ bool setPbcRefToPrevStepCOM);
/* Data only modified at initialization */
const t_pull_group params; /**< The pull group parameters */
When no pbc refence atom is used, this pointer shall be null. */
/* Data, potentially, changed at every pull call */
- real mwscale; /**< mass*weight scaling factor 1/sum w m */
- real wscale; /**< scaling factor for the weights: sum w m/sum w w m */
- real invtm; /**< inverse total mass of the group: 1/wscale sum w m */
- std::vector < gmx::BasicVector < double>> mdw; /**< mass*gradient(weight) for atoms */
- std::vector<double> dv; /**< distance to the other group(s) along vec */
- dvec x; /**< COM before update */
- dvec xp; /**< COM after update before constraining */
+ real mwscale; /**< mass*weight scaling factor 1/sum w m */
+ real wscale; /**< scaling factor for the weights: sum w m/sum w w m */
+ real invtm; /**< inverse total mass of the group: 1/wscale sum w m */
+ std::vector < gmx::BasicVector < double>> mdw; /**< mass*gradient(weight) for atoms */
+ std::vector<double> dv; /**< distance to the other group(s) along vec */
+ dvec x; /**< COM before update */
+ dvec xp; /**< COM after update before constraining */
+ dvec x_prev_step; /**< center of mass of the previous step */
};
/* Struct describing the instantaneous spatial layout of a pull coordinate */
#include "gromacs/math/utilities.h"
#include "gromacs/math/vec.h"
#include "gromacs/mdtypes/commrec.h"
+#include "gromacs/mdtypes/inputrec.h"
#include "gromacs/mdtypes/md_enums.h"
#include "gromacs/mdtypes/mdatom.h"
+#include "gromacs/mdtypes/state.h"
#include "gromacs/pbcutil/pbc.h"
#include "gromacs/pulling/pull.h"
#include "gromacs/utility/fatalerror.h"
for (size_t g = 0; g < pull->group.size(); g++)
{
const pull_group_work_t &group = pull->group[g];
- if (group.needToCalcCom && group.epgrppbc == epgrppbcREFAT)
+ if (group.needToCalcCom && (group.epgrppbc == epgrppbcREFAT || group.epgrppbc == epgrppbcPREVSTEPCOM))
{
setPbcAtomCoords(pull->group[g], x, x_pbc[g]);
numPbcAtoms++;
{
if (pgrp->epgrppbc != epgrppbcCOS)
{
- rvec x_pbc = { 0, 0, 0 };
+ rvec x_pbc = { 0, 0, 0 };
- if (pgrp->epgrppbc == epgrppbcREFAT)
+ switch (pgrp->epgrppbc)
{
- /* Set the pbc atom */
- copy_rvec(comm->pbcAtomBuffer[g], x_pbc);
+ case epgrppbcREFAT:
+ /* Set the pbc atom */
+ copy_rvec(comm->pbcAtomBuffer[g], x_pbc);
+ break;
+ case epgrppbcPREVSTEPCOM:
+ /* Set the pbc reference to the COM of the group of the last step */
+ copy_dvec_to_rvec(pgrp->x_prev_step, comm->pbcAtomBuffer[g]);
+ copy_dvec_to_rvec(pgrp->x_prev_step, x_pbc);
}
/* The final sums should end up in comSums[0] */
{
pgrp->xp[m] = dvecBuffer[1][m]*pgrp->mwscale;
}
- if (pgrp->epgrppbc == epgrppbcREFAT)
+ if (pgrp->epgrppbc == epgrppbcREFAT || pgrp->epgrppbc == epgrppbcPREVSTEPCOM)
{
pgrp->x[m] += comm->pbcAtomBuffer[g][m];
if (xp)
for (size_t g = 0; g < pull.group.size(); g++)
{
const pull_group_work_t &group = pull.group[g];
- if (group.epgrppbc == epgrppbcREFAT &&
+ if ((group.epgrppbc == epgrppbcREFAT || group.epgrppbc == epgrppbcPREVSTEPCOM) &&
!pullGroupObeysPbcRestrictions(group, dimUsed[g], x, pbc, pull.comm.pbcAtomBuffer[g], pbcMargin))
{
return g;
/* Check PBC if the group uses a PBC reference atom treatment. */
const pull_group_work_t &group = pull.group[groupNr];
- if (group.epgrppbc != epgrppbcREFAT)
+ if (group.epgrppbc != epgrppbcREFAT && group.epgrppbc != epgrppbcPREVSTEPCOM)
{
return true;
}
return (pullGroupObeysPbcRestrictions(group, dimUsed, as_rvec_array(x.data()), pbc, pull.comm.pbcAtomBuffer[groupNr], pbcMargin));
}
+
+void setStatePrevStepPullCom(const struct pull_t *pull, t_state *state)
+{
+ for (size_t i = 0; i < state->com_prev_step.size()/DIM; i++)
+ {
+ for (int j = 0; j < DIM; j++)
+ {
+ state->com_prev_step[i*DIM+j] = pull->group[i].x_prev_step[j];
+ }
+ }
+}
+
+void setPrevStepPullComFromState(struct pull_t *pull, const t_state *state)
+{
+ for (size_t i = 0; i < state->com_prev_step.size()/DIM; i++)
+ {
+ for (int j = 0; j < DIM; j++)
+ {
+ pull->group[i].x_prev_step[j] = state->com_prev_step[i*DIM+j];
+ }
+ }
+}
+
+void updatePrevStepCom(struct pull_t *pull)
+{
+ for (size_t g = 0; g < pull->group.size(); g++)
+ {
+ if (pull->group[g].needToCalcCom)
+ {
+ for (int j = 0; j < DIM; j++)
+ {
+ pull->group[g].x_prev_step[j] = pull->group[g].x[j];
+ }
+ }
+ }
+}
+
+void allocStatePrevStepPullCom(t_state *state, pull_t *pull)
+{
+ if (!pull)
+ {
+ state->com_prev_step.clear();
+ return;
+ }
+ size_t ngroup = pull->group.size();
+ if (state->com_prev_step.size()/DIM != ngroup)
+ {
+ state->com_prev_step.resize(ngroup * DIM, NAN);
+ }
+}
+
+void initPullComFromPrevStep(const t_commrec *cr,
+ pull_t *pull,
+ const t_mdatoms *md,
+ t_pbc *pbc,
+ const rvec x[])
+{
+ pull_comm_t *comm = &pull->comm;
+ size_t ngroup = pull->group.size();
+
+ comm->pbcAtomBuffer.resize(ngroup);
+ comm->comBuffer.resize(ngroup*DIM);
+
+ for (size_t g = 0; g < ngroup; g++)
+ {
+ pull_group_work_t *pgrp;
+
+ pgrp = &pull->group[g];
+
+ if (pgrp->needToCalcCom && pgrp->epgrppbc == epgrppbcPREVSTEPCOM)
+ {
+ GMX_ASSERT(pgrp->params.nat > 1, "Groups with no atoms, or only one atom, should not "
+ "use the COM from the previous step as reference.");
+
+ rvec x_pbc = { 0, 0, 0 };
+ pull_set_pbcatoms(cr, pull, x, comm->pbcAtomBuffer);
+ copy_rvec(comm->pbcAtomBuffer[g], x_pbc);
+
+ if (debug)
+ {
+ fprintf(debug, "Initialising prev step COM of pull group %zu. x_pbc =", g);
+ for (int m = 0; m < DIM; m++)
+ {
+ fprintf(debug, " %f", x_pbc[m]);
+ }
+ fprintf(debug, "\n");
+ }
+
+ /* The following is to a large extent similar to pull_calc_coms() */
+
+ /* The final sums should end up in sum_com[0] */
+ ComSums &comSumsTotal = pull->comSums[0];
+
+ if (pgrp->atomSet.numAtomsLocal() <= c_pullMaxNumLocalAtomsSingleThreaded)
+ {
+ sum_com_part(pgrp, 0, pgrp->atomSet.numAtomsLocal(),
+ x, nullptr, md->massT,
+ pbc, x_pbc,
+ &comSumsTotal);
+ }
+ else
+ {
+#pragma omp parallel for num_threads(pull->nthreads) schedule(static)
+ for (int t = 0; t < pull->nthreads; t++)
+ {
+ int ind_start = (pgrp->atomSet.numAtomsLocal()*(t + 0))/pull->nthreads;
+ int ind_end = (pgrp->atomSet.numAtomsLocal()*(t + 1))/pull->nthreads;
+ sum_com_part(pgrp, ind_start, ind_end,
+ x, nullptr, md->massT,
+ pbc, x_pbc,
+ &pull->comSums[t]);
+ }
+
+ /* Reduce the thread contributions to sum_com[0] */
+ for (int t = 1; t < pull->nthreads; t++)
+ {
+ comSumsTotal.sum_wm += pull->comSums[t].sum_wm;
+ comSumsTotal.sum_wwm += pull->comSums[t].sum_wwm;
+ dvec_inc(comSumsTotal.sum_wmx, pull->comSums[t].sum_wmx);
+ dvec_inc(comSumsTotal.sum_wmxp, pull->comSums[t].sum_wmxp);
+ }
+ }
+
+ if (pgrp->localWeights.empty())
+ {
+ comSumsTotal.sum_wwm = comSumsTotal.sum_wm;
+ }
+
+ /* Copy local sums to a buffer for global summing */
+ copy_dvec(comSumsTotal.sum_wmx, comm->comBuffer[g*3]);
+ copy_dvec(comSumsTotal.sum_wmxp, comm->comBuffer[g*3 + 1]);
+ comm->comBuffer[g*3 + 2][0] = comSumsTotal.sum_wm;
+ comm->comBuffer[g*3 + 2][1] = comSumsTotal.sum_wwm;
+ comm->comBuffer[g*3 + 2][2] = 0;
+ }
+ }
+
+ pullAllReduce(cr, comm, ngroup*3*DIM, static_cast<double *>(comm->comBuffer[0]));
+
+ for (size_t g = 0; g < ngroup; g++)
+ {
+ pull_group_work_t *pgrp;
+
+ pgrp = &pull->group[g];
+ if (pgrp->needToCalcCom)
+ {
+ if (pgrp->epgrppbc == epgrppbcPREVSTEPCOM)
+ {
+ double wmass, wwmass;
+
+ /* Determine the inverse mass */
+ wmass = comm->comBuffer[g*3+2][0];
+ wwmass = comm->comBuffer[g*3+2][1];
+ pgrp->mwscale = 1.0/wmass;
+ /* invtm==0 signals a frozen group, so then we should keep it zero */
+ if (pgrp->invtm != 0)
+ {
+ pgrp->wscale = wmass/wwmass;
+ pgrp->invtm = wwmass/(wmass*wmass);
+ }
+ /* Divide by the total mass */
+ for (int m = 0; m < DIM; m++)
+ {
+ pgrp->x[m] = comm->comBuffer[g*3 ][m]*pgrp->mwscale;
+ if (pgrp->epgrppbc == epgrppbcREFAT || pgrp->epgrppbc == epgrppbcPREVSTEPCOM)
+ {
+ pgrp->x[m] += comm->pbcAtomBuffer[g][m];
+ }
+ }
+ if (debug)
+ {
+ fprintf(debug, "Pull group %zu wmass %f invtm %f\n",
+ g, 1.0/pgrp->mwscale, pgrp->invtm);
+ fprintf(debug, "Initialising prev step COM of pull group %zu to", g);
+ for (int m = 0; m < DIM; m++)
+ {
+ fprintf(debug, " %f", pgrp->x[m]);
+ }
+ fprintf(debug, "\n");
+ }
+ copy_dvec(pgrp->x, pgrp->x_prev_step);
+ }
+ }
+ }
+}
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff