if (ir->bPull)
{
/* Update the local pull groups */
- dd_make_local_pull_groups(cr, ir->pull_work, mdatoms);
+ dd_make_local_pull_groups(cr, ir->pull_work);
}
if (ir->bRot)
#include <stdlib.h>
#include <string.h>
+#include "gromacs/domdec/localatomsetmanager.h"
#include "gromacs/fileio/readinp.h"
#include "gromacs/fileio/warninp.h"
#include "gromacs/gmxpreprocess/readir.h"
double t_start;
pull = ir->pull;
- pull_work = init_pull(nullptr, pull, ir, mtop, nullptr, lambda);
- auto mdAtoms = gmx::makeMDAtoms(nullptr, *mtop, *ir, false);
- auto md = mdAtoms->mdatoms();
+ gmx::LocalAtomSetManager atomSets;
+ pull_work = init_pull(nullptr, pull, ir, mtop, nullptr, &atomSets, lambda);
+ auto mdAtoms = gmx::makeMDAtoms(nullptr, *mtop, *ir, false);
+ auto md = mdAtoms->mdatoms();
atoms2md(mtop, ir, -1, nullptr, mtop->natoms, mdAtoms.get());
if (ir->efep)
{
/* Initialize pull code */
inputrec->pull_work =
init_pull(fplog, inputrec->pull, inputrec,
- &mtop, cr, inputrec->fepvals->init_lambda);
+ &mtop, cr, &atomSets, inputrec->fepvals->init_lambda);
if (EI_DYNAMICS(inputrec->eI) && MASTER(cr))
{
init_pull_output_files(inputrec->pull_work,
#include "gromacs/commandline/filenm.h"
#include "gromacs/domdec/domdec_struct.h"
#include "gromacs/domdec/ga2la.h"
+#include "gromacs/domdec/localatomset.h"
+#include "gromacs/domdec/localatomsetmanager.h"
#include "gromacs/fileio/gmxfio.h"
#include "gromacs/fileio/xvgr.h"
#include "gromacs/gmxlib/network.h"
#include "gromacs/mdtypes/md_enums.h"
#include "gromacs/mdtypes/mdatom.h"
#include "gromacs/pbcutil/pbc.h"
-#include "gromacs/pulling/pull_internal.h"
#include "gromacs/topology/mtop_lookup.h"
#include "gromacs/topology/topology.h"
#include "gromacs/utility/cstringutil.h"
#include "gromacs/utility/real.h"
#include "gromacs/utility/smalloc.h"
+#include "pull_internal.h"
+
+static int groupPbcFromParams(const t_pull_group ¶ms)
+{
+ if (params.nat <= 1)
+ {
+ /* no PBC required */
+ return epgrppbcNONE;
+ }
+ else if (params.pbcatom >= 0)
+ {
+ return epgrppbcREFAT;
+ }
+ else
+ {
+ return epgrppbcCOS;
+ }
+}
+
+/* NOTE: The params initialization currently copies pointers.
+ * So the lifetime of the source, currently always inputrec,
+ * should not end before that of this object.
+ * 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) :
+ params(params),
+ epgrppbc(groupPbcFromParams(params)),
+ needToCalcCom(false),
+ atomSet(atomSet),
+ mwscale(0),
+ wscale(0),
+ invtm(0)
+{
+ clear_dvec(x);
+ clear_dvec(xp);
+};
+
bool pull_coordinate_is_angletype(const t_pull_coord *pcrd)
{
return (pcrd->eGeom == epullgANGLE ||
{
double inv_wm = pgrp->mwscale;
+ auto localAtomIndices = pgrp->atomSet.localIndex();
for (int i = ind_start; i < ind_end; i++)
{
- int ii = pgrp->ind_loc[i];
+ int ii = localAtomIndices[i];
double wmass = md->massT[ii];
- if (pgrp->weight_loc)
+ if (!pgrp->localWeights.empty())
{
- wmass *= pgrp->weight_loc[i];
+ wmass *= pgrp->localWeights[i];
}
for (int d = 0; d < DIM; d++)
const dvec f_pull, int sign, rvec *f,
int nthreads)
{
- if (pgrp->params.nat == 1 && pgrp->nat_loc == 1)
+ auto localAtomIndices = pgrp->atomSet.localIndex();
+
+ if (pgrp->params.nat == 1 && pgrp->atomSet.numAtomsLocal() == 1)
{
/* Only one atom and our rank has this atom: we can skip
* the mass weighting, which means that this code also works
*/
for (int d = 0; d < DIM; d++)
{
- f[pgrp->ind_loc[0]][d] += sign*f_pull[d];
+ f[localAtomIndices[0]][d] += sign*f_pull[d];
}
}
else
{
- if (pgrp->nat_loc <= c_pullMaxNumLocalAtomsSingleThreaded)
+ if (localAtomIndices.size() <= c_pullMaxNumLocalAtomsSingleThreaded)
{
- apply_forces_grp_part(pgrp, 0, pgrp->nat_loc, md, f_pull, sign, f);
+ apply_forces_grp_part(pgrp, 0, localAtomIndices.size(), md, f_pull, sign, f);
}
else
{
#pragma omp parallel for num_threads(nthreads) schedule(static)
for (int th = 0; th < nthreads; th++)
{
- int ind_start = (pgrp->nat_loc*(th + 0))/nthreads;
- int ind_end = (pgrp->nat_loc*(th + 1))/nthreads;
+ int ind_start = (localAtomIndices.size()*(th + 0))/nthreads;
+ int ind_end = (localAtomIndices.size()*(th + 1))/nthreads;
apply_forces_grp_part(pgrp, ind_start, ind_end,
md, f_pull, sign, f);
}
{
double inv_wm = pgrp->mwscale;
+ auto localAtomIndices = pgrp->atomSet.localIndex();
+
/* The cylinder group is always a slab in the system, thus large.
* Therefore we always thread-parallelize this group.
*/
+ int numAtomsLocal = localAtomIndices.size();
#pragma omp parallel for num_threads(nthreads) schedule(static)
- for (int i = 0; i < pgrp->nat_loc; i++)
+ for (int i = 0; i < numAtomsLocal; i++)
{
- int ii = pgrp->ind_loc[i];
+ double weight = pgrp->localWeights[i];
+ if (weight == 0)
+ {
+ continue;
+ }
+ int ii = localAtomIndices[i];
double mass = md->massT[ii];
- double weight = pgrp->weight_loc[i];
/* The stored axial distance from the cylinder center (dv) needs
* to be corrected for an offset (dv_corr), which was unknown when
* we calculated dv.
double inpr;
double lambda, rm, invdt = 0;
gmx_bool bConverged_all, bConverged = FALSE;
- int niter = 0, g, ii, j, m, max_iter = 100;
+ int niter = 0, ii, j, m, max_iter = 100;
double a;
dvec tmp, tmp3;
snew(r_ij, pull->coord.size());
snew(dr_tot, pull->coord.size());
- snew(rnew, pull->ngroup);
+ snew(rnew, pull->group.size());
/* copy the current unconstrained positions for use in iterations. We
iterate until rinew[i] and rjnew[j] obey the constraints. Then
rinew - pull.x_unc[i] is the correction dr to group i */
- for (g = 0; g < pull->ngroup; g++)
+ for (size_t g = 0; g < pull->group.size(); g++)
{
copy_dvec(pull->group[g].xp, rnew[g]);
}
{
if (debug)
{
- fprintf(debug, "NOT CONVERGED YET: Group %d:"
+ fprintf(debug, "Pull constraint not converged: "
+ "groups %d %d,"
"d_ref = %f, current d = %f\n",
- g, coord.value_ref, dnorm(unc_ij));
+ coord.params.group[0], coord.params.group[1],
+ coord.value_ref, dnorm(unc_ij));
}
bConverged_all = FALSE;
}
/* update atoms in the groups */
- for (g = 0; g < pull->ngroup; g++)
+ for (size_t g = 0; g < pull->group.size(); g++)
{
const pull_group_work_t *pgrp;
dvec dr;
}
/* update the atom positions */
+ auto localAtomIndices = pgrp->atomSet.localIndex();
copy_dvec(dr, tmp);
- for (j = 0; j < pgrp->nat_loc; j++)
+ for (gmx::index j = 0; j < localAtomIndices.size(); j++)
{
- ii = pgrp->ind_loc[j];
- if (pgrp->weight_loc)
+ ii = localAtomIndices[j];
+ if (!pgrp->localWeights.empty())
{
- dsvmul(pgrp->wscale*pgrp->weight_loc[j], dr, tmp);
+ dsvmul(pgrp->wscale*pgrp->localWeights[j], dr, tmp);
}
for (m = 0; m < DIM; m++)
{
}
}
-static void make_local_pull_group(gmx_ga2la_t *ga2la,
- pull_group_work_t *pg, int start, int end)
-{
- int i, ii;
-
- pg->nat_loc = 0;
- for (i = 0; i < pg->params.nat; i++)
- {
- ii = pg->params.ind[i];
- if (ga2la)
- {
- if (!ga2la_get_home(ga2la, ii, &ii))
- {
- ii = -1;
- }
- }
- if (ii >= start && ii < end)
- {
- /* This is a home atom, add it to the local pull group */
- if (pg->nat_loc >= pg->nalloc_loc)
- {
- pg->nalloc_loc = over_alloc_dd(pg->nat_loc+1);
- srenew(pg->ind_loc, pg->nalloc_loc);
- if (pg->epgrppbc == epgrppbcCOS || pg->params.weight != nullptr)
- {
- srenew(pg->weight_loc, pg->nalloc_loc);
- }
- }
- pg->ind_loc[pg->nat_loc] = ii;
- if (pg->params.weight != nullptr)
- {
- pg->weight_loc[pg->nat_loc] = pg->params.weight[i];
- }
- pg->nat_loc++;
- }
- }
-}
-
-void dd_make_local_pull_groups(const t_commrec *cr, struct pull_t *pull, t_mdatoms *md)
+void dd_make_local_pull_groups(const t_commrec *cr, struct pull_t *pull)
{
gmx_domdec_t *dd;
pull_comm_t *comm;
- gmx_ga2la_t *ga2la;
gmx_bool bMustParticipate;
- int g;
dd = cr->dd;
comm = &pull->comm;
- if (dd)
- {
- ga2la = dd->ga2la;
- }
- else
- {
- ga2la = nullptr;
- }
-
/* We always make the master node participate, such that it can do i/o,
* add the virial and to simplify MC type extensions people might have.
*/
bMustParticipate = (comm->bParticipateAll || comm->isMasterRank);
- for (g = 0; g < pull->ngroup; g++)
+ for (pull_group_work_t &group : pull->group)
{
- int a;
-
- make_local_pull_group(ga2la, &pull->group[g],
- 0, md->homenr);
+ if (group.epgrppbc == epgrppbcCOS || !group.globalWeights.empty())
+ {
+ group.localWeights.resize(group.atomSet.numAtomsLocal());
+ for (size_t i = 0; i < group.atomSet.numAtomsLocal(); ++i)
+ {
+ group.localWeights[i] = group.globalWeights[group.atomSet.collectiveIndex()[i]];
+ }
+ }
GMX_ASSERT(bMustParticipate || dd != nullptr, "Either all ranks (including this rank) participate, or we use DD and need to have access to dd here");
/* We should participate if we have pull or pbc atoms */
if (!bMustParticipate &&
- (pull->group[g].nat_loc > 0 ||
- (pull->group[g].params.pbcatom >= 0 &&
- ga2la_get_home(dd->ga2la, pull->group[g].params.pbcatom, &a))))
+ (group.atomSet.numAtomsLocal() > 0 ||
+ (group.params.pbcatom >= 0 &&
+ ga2la_is_home(dd->ga2la, group.params.pbcatom))))
{
bMustParticipate = TRUE;
}
const gmx_mtop_t *mtop,
const t_inputrec *ir, real lambda)
{
- if (EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD)
- {
- /* There are no masses in the integrator.
- * But we still want to have the correct mass-weighted COMs.
- * So we store the real masses in the weights.
- * We do not set nweight, so these weights do not end up in the tpx file.
- */
- if (pg->params.nweight == 0)
- {
- snew(pg->params.weight, pg->params.nat);
- }
- }
+ /* With EM and BD there are no masses in the integrator.
+ * But we still want to have the correct mass-weighted COMs.
+ * So we store the real masses in the weights.
+ */
+ const bool setWeights = (pg->params.nweight > 0 ||
+ EI_ENERGY_MINIMIZATION(ir->eI) ||
+ ir->eI == eiBD);
- if (cr && PAR(cr))
- {
- pg->nat_loc = 0;
- pg->nalloc_loc = 0;
- pg->ind_loc = nullptr;
- pg->weight_loc = nullptr;
- }
- else
- {
- pg->nat_loc = pg->params.nat;
- pg->ind_loc = pg->params.ind;
- if (pg->epgrppbc == epgrppbcCOS)
- {
- snew(pg->weight_loc, pg->params.nat);
- }
- else
- {
- pg->weight_loc = pg->params.weight;
- }
- }
+ /* In parallel, store we need to extract localWeights from weights at DD time */
+ std::vector<real> &weights = ((cr && PAR(cr)) ? pg->globalWeights : pg->localWeights);
- const gmx_groups_t *groups = &mtop->groups;
+ const gmx_groups_t *groups = &mtop->groups;
/* Count frozen dimensions and (weighted) mass */
int nfrozen = 0;
/* Move the mass to the weight */
w *= m;
m = 1;
- pg->params.weight[i] = w;
}
else if (ir->eI == eiBD)
{
}
w *= m/mbd;
m = mbd;
- pg->params.weight[i] = w;
+ }
+ if (setWeights)
+ {
+ weights.push_back(w);
}
tmass += m;
wmass += m*w;
struct pull_t *
init_pull(FILE *fplog, const pull_params_t *pull_params, const t_inputrec *ir,
- const gmx_mtop_t *mtop, const t_commrec *cr,
+ const gmx_mtop_t *mtop, const t_commrec *cr, gmx::LocalAtomSetManager *atomSets,
real lambda)
{
struct pull_t *pull;
pull->params.group = nullptr;
pull->params.coord = nullptr;
- pull->ngroup = pull_params->ngroup;
- snew(pull->group, pull->ngroup);
+ 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->bPotential = FALSE;
pull->bConstraint = FALSE;
pull->bCylinder = FALSE;
pull->bAngle = FALSE;
- for (int g = 0; g < pull->ngroup; g++)
- {
- pull_group_work_t *pgrp = &pull->group[g];
-
- /* Copy the pull group parameters */
- pgrp->params = pull_params->group[g];
-
- if (g == 0)
- {
- GMX_RELEASE_ASSERT(pgrp->params.nat == 0, "pull group 0 is an absolute reference group and should not contain atoms");
- clear_dvec(pgrp->x);
- clear_dvec(pgrp->xp);
- pgrp->bCalcCOM = FALSE;
- }
-
- /* Avoid pointer copies by allocating and copying arrays */
- snew(pgrp->params.ind, pgrp->params.nat);
- for (int i = 0; i < pgrp->params.nat; i++)
- {
- pgrp->params.ind[i] = pull_params->group[g].ind[i];
- }
- if (pgrp->params.nweight > 0)
- {
- snew(pgrp->params.weight, pgrp->params.nweight);
- for (int i = 0; i < pgrp->params.nweight; i++)
- {
- pgrp->params.weight[i] = pull_params->group[g].weight[i];
- }
- }
- }
+ GMX_RELEASE_ASSERT(pull->group[0].params.nat == 0, "pull group 0 is an absolute reference group and should not contain atoms");
pull->numCoordinatesWithExternalPotential = 0;
/* We only need to calculate the plain COM of a group
* when it is not only used as a cylinder group.
+ * Also the absolute reference group 0 needs no COM computation.
*/
- int groupRangeStart = (pcrd->params.eGeom == epullgCYL ? 1 : 0);
- int groupRangeEnd = pcrd->params.ngroup + 1;
-
- for (int i = groupRangeStart; i < groupRangeEnd; i++)
+ for (int i = 0; i < pcrd->params.ngroup; i++)
{
int groupIndex = pcrd->params.group[i];
- if (groupIndex > 0)
+ if (groupIndex > 0 &&
+ !(pcrd->params.eGeom == epullgCYL && i == 0))
{
- pull->group[groupIndex].bCalcCOM = TRUE;
+ pull->group[groupIndex].needToCalcCom = true;
}
}
fprintf(fplog, "Will apply constraint COM pulling\n");
}
// Don't include the reference group 0 in output, so we report ngroup-1
- GMX_RELEASE_ASSERT(pull->ngroup - 1 > 0, "The reference absolute position pull group should always be present");
+ int numRealGroups = pull->group.size() - 1;
+ GMX_RELEASE_ASSERT(numRealGroups > 0, "The reference absolute position pull group should always be present");
fprintf(fplog, "with %zu pull coordinate%s and %d group%s\n",
pull->coord.size(), pull->coord.size() == 1 ? "" : "s",
- (pull->ngroup - 1), (pull->ngroup - 1) == 1 ? "" : "s");
+ numRealGroups, numRealGroups == 1 ? "" : "s");
if (bAbs)
{
fprintf(fplog, "with an absolute reference\n");
}
bCos = FALSE;
// Don't include the reference group 0 in loop
- for (int g = 1; g < pull->ngroup; g++)
+ for (size_t g = 1; g < pull->group.size(); g++)
{
if (pull->group[g].params.nat > 1 &&
pull->group[g].params.pbcatom < 0)
{
/* We are using cosine weighting */
- fprintf(fplog, "Cosine weighting is used for group %d\n", g);
+ fprintf(fplog, "Cosine weighting is used for group %zu\n", g);
bCos = TRUE;
}
}
pull->bRefAt = FALSE;
pull->cosdim = -1;
- for (int g = 0; g < pull->ngroup; g++)
+ for (size_t g = 0; g < pull->group.size(); g++)
{
pull_group_work_t *pgrp;
pgrp = &pull->group[g];
- pgrp->epgrppbc = epgrppbcNONE;
if (pgrp->params.nat > 0)
{
/* There is an issue when a group is used in multiple coordinates
gmx_bool bGroupUsed = FALSE;
for (int gi = 0; gi < coord.params.ngroup; gi++)
{
- if (coord.params.group[gi] == g)
+ if (coord.params.group[gi] == static_cast<int>(g))
{
bGroupUsed = TRUE;
}
/* Determine if we need to take PBC into account for calculating
* the COM's of the pull groups.
*/
- GMX_RELEASE_ASSERT(pull->npbcdim <= DIM, "npbcdim must be <= the number of dimensions");
- for (int m = 0; m < pull->npbcdim; m++)
+ switch (pgrp->epgrppbc)
{
- GMX_RELEASE_ASSERT(m <= DIM, "npbcdim must be <= the number of dimensions");
- if (pulldim[m] == 1 && pgrp->params.nat > 1)
- {
- if (pgrp->params.pbcatom >= 0)
+ case epgrppbcREFAT:
+ pull->bRefAt = TRUE;
+ break;
+ case epgrppbcCOS:
+ if (pgrp->params.weight != nullptr)
{
- pgrp->epgrppbc = epgrppbcREFAT;
- pull->bRefAt = TRUE;
+ gmx_fatal(FARGS, "Pull groups can not have relative weights and cosine weighting at same time");
}
- else
+ for (int m = 0; m < DIM; m++)
{
- if (pgrp->params.weight != nullptr)
+ if (m < pull->npbcdim && pulldim[m] == 1)
{
- gmx_fatal(FARGS, "Pull groups can not have relative weights and cosine weighting at same time");
- }
- pgrp->epgrppbc = epgrppbcCOS;
- if (pull->cosdim >= 0 && pull->cosdim != m)
- {
- gmx_fatal(FARGS, "Can only use cosine weighting with pulling in one dimension (use mdp option pull-coord?-dim)");
+ if (pull->cosdim >= 0 && pull->cosdim != m)
+ {
+ gmx_fatal(FARGS, "Can only use cosine weighting with pulling in one dimension (use mdp option pull-coord?-dim)");
+ }
+ pull->cosdim = m;
}
- pull->cosdim = m;
}
- }
+ break;
+ case epgrppbcNONE:
+ break;
}
+
/* Set the indices */
init_pull_group_index(fplog, cr, g, pgrp,
bConstraint, pulldim_con,
/* If we use cylinder coordinates, do some initialising for them */
if (pull->bCylinder)
{
- snew(pull->dyna, pull->coord.size());
-
for (const pull_coord_work_t &coord : pull->coord)
{
if (coord.params.eGeom == epullgCYL)
gmx_fatal(FARGS, "A cylinder pull group is not supported when using absolute reference!\n");
}
}
+ const auto &referenceGroup = pull->group[coord.params.group[0]];
+ pull->dyna.emplace_back(referenceGroup.params, referenceGroup.atomSet);
}
}
- else
- {
- pull->dyna = nullptr;
- }
/* The gmx_omp_nthreads module might not be initialized here, so max(1,) */
pull->nthreads = std::max(1, gmx_omp_nthreads_get(emntDefault));
}
}
-static void destroy_pull_group(pull_group_work_t *pgrp)
-{
- if (pgrp->ind_loc != pgrp->params.ind)
- {
- sfree(pgrp->ind_loc);
- }
- if (pgrp->weight_loc != pgrp->params.weight)
- {
- sfree(pgrp->weight_loc);
- }
- sfree(pgrp->mdw);
- sfree(pgrp->dv);
-
- sfree(pgrp->params.ind);
- sfree(pgrp->params.weight);
-}
-
static void destroy_pull(struct pull_t *pull)
{
- for (int i = 0; i < pull->ngroup; i++)
- {
- destroy_pull_group(&pull->group[i]);
- }
- for (size_t c = 0; c < pull->coord.size(); c++)
- {
- if (pull->coord[c].params.eGeom == epullgCYL)
- {
- destroy_pull_group(&(pull->dyna[c]));
- }
- }
- sfree(pull->group);
- sfree(pull->dyna);
-
#if GMX_MPI
if (pull->comm.mpi_comm_com != MPI_COMM_NULL)
{
#include "gromacs/math/vectypes.h"
#include "gromacs/mdtypes/pull-params.h"
-#include "gromacs/pulling/pull_internal.h"
#include "gromacs/utility/basedefinitions.h"
#include "gromacs/utility/real.h"
struct ContinuationOptions;
struct gmx_mtop_t;
struct gmx_output_env_t;
+struct pull_coord_work_t;
struct pull_params_t;
struct t_commrec;
struct t_filenm;
namespace gmx
{
class ForceWithVirial;
+class LocalAtomSetManager;
}
/*! \brief Returns if the pull coordinate is an angle
*
* \param cr Structure for communication info.
* \param pull The pull group.
- * \param md All atoms.
*/
-void dd_make_local_pull_groups(const t_commrec *cr,
- struct pull_t *pull, t_mdatoms *md);
+void dd_make_local_pull_groups(const t_commrec *cr, struct pull_t *pull);
/*! \brief Allocate, initialize and return a pull work struct.
* \param ir The inputrec.
* \param mtop The topology of the whole system.
* \param cr Struct for communication info.
+ * \param atomSets The manager that handles the pull atom sets
* \param lambda FEP lambda.
*/
struct pull_t *init_pull(FILE *fplog,
const t_inputrec *ir,
const gmx_mtop_t *mtop,
const t_commrec *cr,
+ gmx::LocalAtomSetManager *atomSets,
real lambda);
/*! \brief Set up and open the pull output files, when requested.
* the research papers on the package. Check out http://www.gromacs.org.
*/
-/*! \libinternal \file
+/*! \internal \file
*
*
* \brief
use in the pull code.
*
* \author Berk Hess
- *
- * \inlibraryapi
*/
#ifndef GMX_PULLING_PULL_INTERNAL_H
#include <vector>
+#include "gromacs/domdec/localatomset.h"
#include "gromacs/mdtypes/pull-params.h"
#include "gromacs/utility/gmxmpi.h"
-/*! \cond INTERNAL */
-
/*! \brief Determines up to what local atom count a pull group gets processed single-threaded.
*
* We set this limit to 1 with debug to catch bugs.
epgrppbcNONE, epgrppbcREFAT, epgrppbcCOS
};
-typedef struct
+/*! \internal
+ * \brief Pull group data used during pulling
+ */
+struct pull_group_work_t
{
- t_pull_group params;
-
- gmx_bool bCalcCOM; /* Calculate COM? Not if only used as cylinder group */
- int epgrppbc; /* The type of pbc for this pull group, see enum above */
-
- int nat_loc; /* Number of local pull atoms */
- int nalloc_loc; /* Allocation size for ind_loc and weight_loc */
- int *ind_loc; /* Local pull indices */
- real *weight_loc; /* Weights for the local indices */
-
- 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 */
- dvec *mdw; /* mass*gradient(weight) for atoms */
- double *dv; /* distance to the other group along vec */
- dvec x; /* center of mass before update */
- dvec xp; /* center of mass after update before constraining */
-}
-pull_group_work_t;
+ /*! \brief Constructor
+ *
+ * \param[in] params The group parameters set by the user
+ * \param[in] atomSet The global to local atom set manager
+ */
+ pull_group_work_t(const t_pull_group ¶ms,
+ gmx::LocalAtomSet atomSet);
+
+ /* Data only modified at initialization */
+ const t_pull_group params; /**< The pull group parameters */
+ const int epgrppbc; /**< The type of pbc for this pull group, see enum above */
+ bool needToCalcCom; /**< Do we need to calculate the COM? (Not for group 0 or if only used as cylinder group) */
+ std::vector<real> globalWeights; /**< Weights per atom set by the user and/or mass/friction coefficients, if empty all weights are equal */
+
+ /* Data modified only at init or at domain decomposition */
+ gmx::LocalAtomSet atomSet; /**< Global to local atom set mapper */
+ std::vector<real> localWeights; /**< Weights for the local atoms */
+
+ /* 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 */
+};
/* Struct describing the instantaneous spatial layout of a pull coordinate */
struct PullCoordSpatialData
gmx_bool bCylinder; /* Is group 0 a cylinder group? */
/* Parameters + dynamic data for groups */
- int ngroup; /* Number of pull groups */
- pull_group_work_t *group; /* The pull group param and work data */
- pull_group_work_t *dyna; /* Dynamic groups for geom=cylinder */
+ std::vector<pull_group_work_t> group; /* The pull group param and work data */
+ std::vector<pull_group_work_t> dyna; /* Dynamic groups for geom=cylinder */
/* Parameters + dynamic data for coordinates */
std::vector<pull_coord_work_t> coord; /* The pull group param and work data */
int numExternalPotentialsStillToBeAppliedThisStep;
};
-/*! \endcond */
-
#endif
#include "gromacs/mdtypes/mdatom.h"
#include "gromacs/pbcutil/pbc.h"
#include "gromacs/pulling/pull.h"
-#include "gromacs/pulling/pull_internal.h"
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/futil.h"
#include "gromacs/utility/gmxassert.h"
#include "gromacs/utility/real.h"
#include "gromacs/utility/smalloc.h"
+#include "pull_internal.h"
+
#if GMX_MPI
// Helper function to deduce MPI datatype from the type of data
const rvec *x,
rvec *x_pbc)
{
- int g, n;
-
- n = 0;
- for (g = 0; g < pull->ngroup; g++)
+ int n = 0;
+ for (size_t g = 0; g < pull->group.size(); g++)
{
- if (!pull->group[g].bCalcCOM || pull->group[g].params.pbcatom == -1)
+ if (!pull->group[g].needToCalcCom ||
+ pull->group[g].params.pbcatom == -1)
{
clear_rvec(x_pbc[g]);
}
* This can be very expensive at high parallelization, so we only
* do this after each DD repartitioning.
*/
- pullAllReduce(cr, &pull->comm, pull->ngroup*DIM, x_pbc[0]);
+ pullAllReduce(cr, &pull->comm, pull->group.size()*DIM, x_pbc[0]);
}
}
{
/* The size and stride per coord for the reduction buffer */
const int stride = 9;
- int i, ii, m, start, end;
- rvec g_x, dx, dir;
double inv_cyl_r2;
pull_comm_t *comm;
- gmx_ga2la_t *ga2la = nullptr;
comm = &pull->comm;
snew(comm->dbuf_cyl, pull->coord.size()*stride);
}
- if (cr && DOMAINDECOMP(cr))
- {
- ga2la = cr->dd->ga2la;
- }
-
- start = 0;
- end = md->homenr;
-
inv_cyl_r2 = 1.0/gmx::square(pull->params.cylinder_r);
/* loop over all groups to make a reference group for each*/
if (pcrd->params.eGeom == epullgCYL)
{
- pull_group_work_t *pref, *pgrp, *pdyna;
-
/* pref will be the same group for all pull coordinates */
- pref = &pull->group[pcrd->params.group[0]];
- pgrp = &pull->group[pcrd->params.group[1]];
- pdyna = &pull->dyna[c];
- copy_dvec_to_rvec(pcrd->spatialData.vec, dir);
- pdyna->nat_loc = 0;
-
- /* We calculate distances with respect to the reference location
- * of this cylinder group (g_x), which we already have now since
- * we reduced the other group COM over the ranks. This resolves
+ const pull_group_work_t &pref = pull->group[pcrd->params.group[0]];
+ const pull_group_work_t &pgrp = pull->group[pcrd->params.group[1]];
+ pull_group_work_t &pdyna = pull->dyna[c];
+ rvec direction;
+ copy_dvec_to_rvec(pcrd->spatialData.vec, direction);
+
+ /* Since we have not calculated the COM of the cylinder group yet,
+ * we calculate distances with respect to location of the pull
+ * group minus the reference position along the vector.
+ * here we already have the COM of the pull group. This resolves
* any PBC issues and we don't need to use a PBC-atom here.
*/
if (pcrd->params.rate != 0)
/* With rate=0, value_ref is set initially */
pcrd->value_ref = pcrd->params.init + pcrd->params.rate*t;
}
- for (m = 0; m < DIM; m++)
+ rvec reference;
+ for (int m = 0; m < DIM; m++)
{
- g_x[m] = pgrp->x[m] - pcrd->spatialData.vec[m]*pcrd->value_ref;
+ reference[m] = pgrp.x[m] - pcrd->spatialData.vec[m]*pcrd->value_ref;
}
+ auto localAtomIndices = pref.atomSet.localIndex();
+
+ /* This actually only needs to be done at init or DD time,
+ * but resizing with the same size does not cause much overhead.
+ */
+ pdyna.localWeights.resize(localAtomIndices.size());
+ pdyna.mdw.resize(localAtomIndices.size());
+ pdyna.dv.resize(localAtomIndices.size());
+
/* loop over all atoms in the main ref group */
- for (i = 0; i < pref->params.nat; i++)
+ for (gmx::index indexInSet = 0; indexInSet < localAtomIndices.size(); indexInSet++)
{
- ii = pref->params.ind[i];
- if (ga2la)
+ int atomIndex = localAtomIndices[indexInSet];
+ rvec dx;
+ pbc_dx_aiuc(pbc, x[atomIndex], reference, dx);
+ double axialLocation = iprod(direction, dx);
+ dvec radialLocation;
+ double dr2 = 0;
+ for (int m = 0; m < DIM; m++)
{
- if (!ga2la_get_home(ga2la, pref->params.ind[i], &ii))
- {
- ii = -1;
- }
+ /* Determine the radial components */
+ radialLocation[m] = dx[m] - axialLocation*direction[m];
+ dr2 += gmx::square(radialLocation[m]);
}
- if (ii >= start && ii < end)
- {
- double dr2, dr2_rel, inp;
- dvec dr;
+ double dr2_rel = dr2*inv_cyl_r2;
- pbc_dx_aiuc(pbc, x[ii], g_x, dx);
- inp = iprod(dir, dx);
- dr2 = 0;
- for (m = 0; m < DIM; m++)
+ if (dr2_rel < 1)
+ {
+ /* add atom to sum of COM and to weight array */
+
+ double mass = md->massT[atomIndex];
+ /* The radial weight function is 1-2x^2+x^4,
+ * where x=r/cylinder_r. Since this function depends
+ * on the radial component, we also get radial forces
+ * on both groups.
+ */
+ double weight = 1 + (-2 + dr2_rel)*dr2_rel;
+ double dweight_r = (-4 + 4*dr2_rel)*inv_cyl_r2;
+ pdyna.localWeights[indexInSet] = weight;
+ sum_a += mass*weight*axialLocation;
+ wmass += mass*weight;
+ wwmass += mass*weight*weight;
+ dvec mdw;
+ dsvmul(mass*dweight_r, radialLocation, mdw);
+ copy_dvec(mdw, pdyna.mdw[indexInSet]);
+ /* Currently we only have the axial component of the
+ * offset from the cylinder COM up to an unkown offset.
+ * We add this offset after the reduction needed
+ * for determining the COM of the cylinder group.
+ */
+ pdyna.dv[indexInSet] = axialLocation;
+ for (int m = 0; m < DIM; m++)
{
- /* Determine the radial components */
- dr[m] = dx[m] - inp*dir[m];
- dr2 += dr[m]*dr[m];
- }
- dr2_rel = dr2*inv_cyl_r2;
-
- if (dr2_rel < 1)
- {
- double mass, weight, dweight_r;
- dvec mdw;
-
- /* add to index, to sum of COM, to weight array */
- if (pdyna->nat_loc >= pdyna->nalloc_loc)
- {
- pdyna->nalloc_loc = over_alloc_large(pdyna->nat_loc+1);
- srenew(pdyna->ind_loc, pdyna->nalloc_loc);
- srenew(pdyna->weight_loc, pdyna->nalloc_loc);
- srenew(pdyna->mdw, pdyna->nalloc_loc);
- srenew(pdyna->dv, pdyna->nalloc_loc);
- }
- pdyna->ind_loc[pdyna->nat_loc] = ii;
-
- mass = md->massT[ii];
- /* The radial weight function is 1-2x^2+x^4,
- * where x=r/cylinder_r. Since this function depends
- * on the radial component, we also get radial forces
- * on both groups.
- */
- weight = 1 + (-2 + dr2_rel)*dr2_rel;
- dweight_r = (-4 + 4*dr2_rel)*inv_cyl_r2;
- pdyna->weight_loc[pdyna->nat_loc] = weight;
- sum_a += mass*weight*inp;
- wmass += mass*weight;
- wwmass += mass*weight*weight;
- dsvmul(mass*dweight_r, dr, mdw);
- copy_dvec(mdw, pdyna->mdw[pdyna->nat_loc]);
- /* Currently we only have the axial component of the
- * distance (inp) up to an unkown offset. We add this
- * offset after the reduction needs to determine the
- * COM of the cylinder group.
- */
- pdyna->dv[pdyna->nat_loc] = inp;
- for (m = 0; m < DIM; m++)
- {
- radf_fac0[m] += mdw[m];
- radf_fac1[m] += mdw[m]*inp;
- }
- pdyna->nat_loc++;
+ radf_fac0[m] += mdw[m];
+ radf_fac1[m] += mdw[m]*axialLocation;
}
}
+ else
+ {
+ pdyna.localWeights[indexInSet] = 0;
+ }
}
}
comm->dbuf_cyl[c*stride+0] = wmass;
if (pcrd->params.eGeom == epullgCYL)
{
- pull_group_work_t *pdyna = &pull->dyna[c];
- pull_group_work_t *pgrp = &pull->group[pcrd->params.group[1]];
- PullCoordSpatialData &spatialData = pcrd->spatialData;
+ pull_group_work_t *pdyna = &pull->dyna[c];
+ pull_group_work_t *pgrp = &pull->group[pcrd->params.group[1]];
+ PullCoordSpatialData &spatialData = pcrd->spatialData;
- double wmass = comm->dbuf_cyl[c*stride+0];
- double wwmass = comm->dbuf_cyl[c*stride+1];
+ double wmass = comm->dbuf_cyl[c*stride+0];
+ double wwmass = comm->dbuf_cyl[c*stride+1];
pdyna->mwscale = 1.0/wmass;
/* Cylinder pulling can't be used with constraints, but we set
* wscale and invtm anyhow, in case someone would like to use them.
* to the atoms in the cylinder group.
*/
spatialData.cyl_dev = 0;
- for (m = 0; m < DIM; m++)
+ for (int m = 0; m < DIM; m++)
{
- g_x[m] = pgrp->x[m] - spatialData.vec[m]*pcrd->value_ref;
+ double reference = pgrp->x[m] - spatialData.vec[m]*pcrd->value_ref;
double dist = -spatialData.vec[m]*comm->dbuf_cyl[c*stride+2]*pdyna->mwscale;
- pdyna->x[m] = g_x[m] - dist;
+ pdyna->x[m] = reference - dist;
spatialData.cyl_dev += dist;
}
/* Now we know the exact COM of the cylinder reference group,
* multiplied with the axial pull force will give the radial
* force on the pulled (non-cylinder) group.
*/
- for (m = 0; m < DIM; m++)
+ for (int m = 0; m < DIM; m++)
{
spatialData.ffrad[m] = (comm->dbuf_cyl[c*stride+6+m] +
comm->dbuf_cyl[c*stride+3+m]*spatialData.cyl_dev)/wmass;
dvec sum_wmx = { 0, 0, 0 };
dvec sum_wmxp = { 0, 0, 0 };
+ auto localAtomIndices = pgrp->atomSet.localIndex();
for (int i = ind_start; i < ind_end; i++)
{
- int ii = pgrp->ind_loc[i];
+ int ii = localAtomIndices[i];
real wm;
- if (pgrp->weight_loc == nullptr)
+ if (pgrp->localWeights.empty())
{
wm = mass[ii];
sum_wm += wm;
{
real w;
- w = pgrp->weight_loc[i];
+ w = pgrp->localWeights[i];
wm = w*mass[ii];
sum_wm += wm;
sum_wwm += wm*w;
double sum_cmp = 0;
double sum_smp = 0;
+ auto localAtomIndices = pgrp->atomSet.localIndex();
+
for (int i = ind_start; i < ind_end; i++)
{
- int ii = pgrp->ind_loc[i];
+ int ii = localAtomIndices[i];
real m = mass[ii];
/* Determine cos and sin sums */
real cw = std::cos(x[ii][cosdim]*twopi_box);
double t,
const rvec x[], rvec *xp)
{
- int g;
real twopi_box = 0;
pull_comm_t *comm;
if (comm->rbuf == nullptr)
{
- snew(comm->rbuf, pull->ngroup);
+ snew(comm->rbuf, pull->group.size());
}
if (comm->dbuf == nullptr)
{
- snew(comm->dbuf, 3*pull->ngroup);
+ snew(comm->dbuf, pull->group.size()*DIM);
}
if (pull->bRefAt && pull->bSetPBCatoms)
twopi_box = 2.0*M_PI/pbc->box[pull->cosdim][pull->cosdim];
}
- for (g = 0; g < pull->ngroup; g++)
+ for (size_t g = 0; g < pull->group.size(); g++)
{
pull_group_work_t *pgrp;
pgrp = &pull->group[g];
- if (pgrp->bCalcCOM)
+ if (pgrp->needToCalcCom)
{
if (pgrp->epgrppbc != epgrppbcCOS)
{
* in that case a check group mass != 0 has been done before.
*/
if (pgrp->params.nat == 1 &&
- pgrp->nat_loc == 1 &&
- md->massT[pgrp->ind_loc[0]] == 0)
+ pgrp->atomSet.numAtomsLocal() == 1 &&
+ md->massT[pgrp->atomSet.localIndex()[0]] == 0)
{
GMX_ASSERT(xp == nullptr, "We should not have groups with zero mass with constraints, i.e. xp!=NULL");
/* Copy the single atom coordinate */
for (int d = 0; d < DIM; d++)
{
- sum_com->sum_wmx[d] = x[pgrp->ind_loc[0]][d];
+ sum_com->sum_wmx[d] = x[pgrp->atomSet.localIndex()[0]][d];
}
/* Set all mass factors to 1 to get the correct COM */
sum_com->sum_wm = 1;
sum_com->sum_wwm = 1;
}
- else if (pgrp->nat_loc <= c_pullMaxNumLocalAtomsSingleThreaded)
+ else if (pgrp->atomSet.numAtomsLocal() <= c_pullMaxNumLocalAtomsSingleThreaded)
{
- sum_com_part(pgrp, 0, pgrp->nat_loc,
+ sum_com_part(pgrp, 0, pgrp->atomSet.numAtomsLocal(),
x, xp, md->massT,
pbc, x_pbc,
sum_com);
#pragma omp parallel for num_threads(pull->nthreads) schedule(static)
for (int t = 0; t < pull->nthreads; t++)
{
- int ind_start = (pgrp->nat_loc*(t + 0))/pull->nthreads;
- int ind_end = (pgrp->nat_loc*(t + 1))/pull->nthreads;
+ 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, xp, md->massT,
pbc, x_pbc,
}
}
- if (pgrp->weight_loc == nullptr)
+ if (pgrp->localWeights.empty())
{
sum_com->sum_wwm = sum_com->sum_wm;
}
#pragma omp parallel for num_threads(pull->nthreads) schedule(static)
for (int t = 0; t < pull->nthreads; t++)
{
- int ind_start = (pgrp->nat_loc*(t + 0))/pull->nthreads;
- int ind_end = (pgrp->nat_loc*(t + 1))/pull->nthreads;
+ int ind_start = (pgrp->atomSet.numAtomsLocal()*(t + 0))/pull->nthreads;
+ int ind_end = (pgrp->atomSet.numAtomsLocal()*(t + 1))/pull->nthreads;
sum_com_part_cosweight(pgrp, ind_start, ind_end,
pull->cosdim, twopi_box,
x, xp, md->massT,
}
}
- pullAllReduce(cr, comm, pull->ngroup*3*DIM, comm->dbuf[0]);
+ pullAllReduce(cr, comm, pull->group.size()*3*DIM, comm->dbuf[0]);
- for (g = 0; g < pull->ngroup; g++)
+ for (size_t g = 0; g < pull->group.size(); g++)
{
pull_group_work_t *pgrp;
pgrp = &pull->group[g];
- if (pgrp->bCalcCOM)
+ if (pgrp->needToCalcCom)
{
GMX_ASSERT(pgrp->params.nat > 0, "Normal pull groups should have atoms, only group 0, which should have bCalcCom=FALSE has nat=0");
{
/* Cosine weighting geometry */
double csw, snw, wmass, wwmass;
- int i, ii;
/* Determine the optimal location of the cosine weight */
csw = comm->dbuf[g*3][0];
/* Set the weights for the local atoms */
csw *= pgrp->invtm;
snw *= pgrp->invtm;
- for (i = 0; i < pgrp->nat_loc; i++)
+ for (size_t i = 0; i < pgrp->atomSet.numAtomsLocal(); i++)
{
- ii = pgrp->ind_loc[i];
- pgrp->weight_loc[i] = csw*std::cos(twopi_box*x[ii][pull->cosdim]) +
+ int ii = pgrp->atomSet.localIndex()[i];
+ pgrp->localWeights[i] = csw*std::cos(twopi_box*x[ii][pull->cosdim]) +
snw*std::sin(twopi_box*x[ii][pull->cosdim]);
}
if (xp)
}
if (debug)
{
- fprintf(debug, "Pull group %d wmass %f invtm %f\n",
+ fprintf(debug, "Pull group %zu wmass %f invtm %f\n",
g, 1.0/pgrp->mwscale, pgrp->invtm);
}
}
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff