#include "gromacs/utility/real.h"
#include "gromacs/utility/smalloc.h"
+#include "domdec_constraints.h"
+#include "domdec_internal.h"
+#include "domdec_vsite.h"
+
#define DDRANK(dd, rank) (rank)
#define DDMASTERRANK(dd) (dd->masterrank)
gmx_fio_fclose(out);
}
-real dd_cutoff_mbody(gmx_domdec_t *dd)
+real dd_cutoff_multibody(const gmx_domdec_t *dd)
{
gmx_domdec_comm_t *comm;
int di;
return r;
}
-real dd_cutoff_twobody(gmx_domdec_t *dd)
+real dd_cutoff_twobody(const gmx_domdec_t *dd)
{
real r_mb;
- r_mb = dd_cutoff_mbody(dd);
+ r_mb = dd_cutoff_multibody(dd);
return std::max(dd->comm->cutoff, r_mb);
}
}
}
-gmx_bool gmx_pmeonlynode(t_commrec *cr, int sim_nodeid)
-{
- gmx_bool bPMEOnlyNode;
-
- if (DOMAINDECOMP(cr))
- {
- bPMEOnlyNode = (dd_simnode2pmenode(cr, sim_nodeid) == -1);
- }
- else
- {
- bPMEOnlyNode = FALSE;
- }
-
- return bPMEOnlyNode;
-}
-
void get_pme_ddnodes(t_commrec *cr, int pmenodeid,
int *nmy_ddnodes, int **my_ddnodes, int *node_peer)
{
* To help us fund GROMACS development, we humbly ask that you cite
* the research papers on the package. Check out http://www.gromacs.org.
*/
+/*! \defgroup module_domdec Spatial domain decomposition (for parallelization over MPI)
+ * \ingroup group_mdrun
+ *
+ * \brief Manages the decomposition of the simulation volume over MPI
+ * ranks to try to distribute work evenly with minimal communication
+ * overheads.
+ *
+ * \todo Get domdec stuff out of legacyheaders/types/commrec.h
+ *
+ * \author Berk Hess <hess@kth.se>
+ *
+ */
-#ifndef _domdec_h
-#define _domdec_h
+/*! \libinternal \file
+ *
+ * \brief This file declares functions for mdrun to call to manage the
+ * details of its domain decomposition.
+ *
+ * \author Berk Hess <hess@kth.se>
+ * \inlibraryapi
+ * \ingroup module_domdec
+ */
+
+#ifndef GMX_DOMDEC_DOMDEC_H
+#define GMX_DOMDEC_DOMDEC_H
#include <stdio.h>
extern "C" {
#endif
-int ddglatnr(gmx_domdec_t *dd, int i);
-/* Returns the global topology atom number belonging to local atom index i.
- * This function is intended for writing ascii output
+/*! \brief Returns the global topology atom number belonging to local atom index i.
+ *
+ * This function is intended for writing ASCII output
* and returns atom numbers starting at 1.
* When dd=NULL returns i+1.
*/
+int ddglatnr(gmx_domdec_t *dd, int i);
+/*! \brief Return a block struct for the charge groups of the whole system */
t_block *dd_charge_groups_global(gmx_domdec_t *dd);
-/* Return a block struct for the charge groups of the whole system */
-
-gmx_bool dd_filled_nsgrid_home(gmx_domdec_t *dd);
-/* Is the ns grid already filled with the home particles? */
-void dd_store_state(gmx_domdec_t *dd, t_state *state);
-/* Store the global cg indices of the home cgs in state,
- * so it can be reset, even after a new DD partitioning.
+/*! \brief Store the global cg indices of the home cgs in state,
+ *
+ * This means it can be reset, even after a new DD partitioning.
*/
+void dd_store_state(gmx_domdec_t *dd, t_state *state);
+/*! \brief Returns a pointer to the gmx_domdec_zones_t struct */
gmx_domdec_zones_t *domdec_zones(gmx_domdec_t *dd);
+/*! \brief Sets the j-charge-group range for i-charge-group \p icg */
void dd_get_ns_ranges(gmx_domdec_t *dd, int icg,
int *jcg0, int *jcg1, ivec shift0, ivec shift1);
+/*! \brief Returns the atom range in the local state for atoms involved in virtual sites */
int dd_natoms_vsite(gmx_domdec_t *dd);
+/*! \brief Sets the atom range for atom in the local state for atoms received in constraints communication */
void dd_get_constraint_range(gmx_domdec_t *dd,
int *at_start, int *at_end);
-real dd_cutoff_mbody(gmx_domdec_t *dd);
-
-real dd_cutoff_twobody(gmx_domdec_t *dd);
-
+/*! \brief Get the number of PME nodes along x and y, can be called with dd=NULL */
void get_pme_nnodes(const gmx_domdec_t *dd,
int *npmenodes_x, int *npmenodes_y);
-/* Get the number of PME nodes along x and y, can be called with dd=NULL */
-
-gmx_bool gmx_pmeonlynode(t_commrec *cr, int nodeid);
-/* Return if nodeid in cr->mpi_comm_mysim is a PME-only node */
+/*! \brief Returns the set of DD nodes that communicate with pme node cr->nodeid */
void get_pme_ddnodes(t_commrec *cr, int pmenodeid,
int *nmy_ddnodes, int **my_ddnodes, int *node_peer);
-/* Returns the set of DD nodes that communicate with pme node cr->nodeid */
+/*! \brief Returns the maximum shift for coordinate communication in PME, dim x */
int dd_pme_maxshift_x(gmx_domdec_t *dd);
-/* Returns the maximum shift for coordinate communication in PME, dim x */
+/*! \brief Returns the maximum shift for coordinate communication in PME, dim y */
int dd_pme_maxshift_y(gmx_domdec_t *dd);
-/* Returns the maximum shift for coordinate communication in PME, dim y */
+/*! \brief Generates the MPI communicators for domain decomposition */
void make_dd_communicators(FILE *fplog, t_commrec *cr, int dd_node_order);
+/*! \brief Initialized the domain decomposition, chooses the DD grid and PME ranks */
gmx_domdec_t *
init_domain_decomposition(FILE *fplog,
t_commrec *cr,
gmx_ddbox_t *ddbox,
int *npme_x, int *npme_y);
+/*! \brief Initialize data structures for bonded interactions */
void dd_init_bondeds(FILE *fplog,
gmx_domdec_t *dd, gmx_mtop_t *mtop,
gmx_vsite_t *vsite,
t_inputrec *ir, gmx_bool bBCheck, cginfo_mb_t *cginfo_mb);
-/* Initialize data structures for bonded interactions */
+/*! \brief Returns if we need to do pbc for calculating bonded interactions */
gmx_bool dd_bonded_molpbc(gmx_domdec_t *dd, int ePBC);
-/* Returns if we need to do pbc for calculating bonded interactions */
+/*! \brief Set DD grid dimensions and limits
+ *
+ * Should be called after calling dd_init_bondeds.
+ */
void set_dd_parameters(FILE *fplog, gmx_domdec_t *dd, real dlb_scale,
t_inputrec *ir,
gmx_ddbox_t *ddbox);
-/* Set DD grid dimensions and limits,
- * should be called after calling dd_init_bondeds.
- */
-gmx_bool change_dd_cutoff(t_commrec *cr, t_state *state, t_inputrec *ir,
- real cutoff_req );
-/* Change the DD non-bonded communication cut-off.
+/*! \brief Change the DD non-bonded communication cut-off.
+ *
* This could fail when trying to increase the cut-off,
* then FALSE will be returned and the cut-off is not modified.
*/
+gmx_bool change_dd_cutoff(t_commrec *cr, t_state *state, t_inputrec *ir,
+ real cutoff_req );
-void change_dd_dlb_cutoff_limit(t_commrec *cr);
-/* Domain boundary changes due to the DD dynamic load balancing can limit
+/*! \brief Limit DLB to preserve the option of returning to the current cut-off.
+ *
+ * Domain boundary changes due to the DD dynamic load balancing can limit
* the cut-off distance that can be set in change_dd_cutoff. This function
* limits the DLB such that using the currently set cut-off should still be
* possible after subsequently setting a shorter cut-off with change_dd_cutoff.
*/
+void change_dd_dlb_cutoff_limit(t_commrec *cr);
+/*! \brief Return if the DLB lock is set */
gmx_bool dd_dlb_is_locked(const gmx_domdec_t *dd);
-/* Return if the DLB lock is set */
+/*! \brief Set a lock such that with DLB=auto DLB can (not) get turned on */
void dd_dlb_set_lock(gmx_domdec_t *dd, gmx_bool bValue);
-/* Set a lock such that with DLB=auto DLB can (not) get turned on */
-void dd_setup_dlb_resource_sharing(t_commrec *cr,
- const gmx_hw_info_t *hwinfo,
- const gmx_hw_opt_t *hw_opt);
-/* When domains (PP MPI ranks) share a GPU, the individual GPU wait times
+/*! \brief Set up communication for averaging GPU wait times over ranks
+ *
+ * When domains (PP MPI ranks) share a GPU, the individual GPU wait times
* are meaningless, as it depends on the order in which tasks on the same
* GPU finish. Therefore there wait times need to be averaged over the ranks
* sharing the same GPU. This function sets up the communication for that.
*/
+void dd_setup_dlb_resource_sharing(t_commrec *cr,
+ const gmx_hw_info_t *hwinfo,
+ const gmx_hw_opt_t *hw_opt);
+/*! \brief Sets up the DD communication setup */
void setup_dd_grid(FILE *fplog, gmx_domdec_t *dd);
+/*! \brief Collects local rvec arrays \p lv to \p v on the master rank */
void dd_collect_vec(gmx_domdec_t *dd,
t_state *state_local, rvec *lv, rvec *v);
+/*! \brief Collects the local state \p state_local to \p state on the master rank */
void dd_collect_state(gmx_domdec_t *dd,
t_state *state_local, t_state *state);
+/*! \brief Cycle counter indices used internally in the domain decomposition */
enum {
ddCyclStep, ddCyclPPduringPME, ddCyclF, ddCyclWaitGPU, ddCyclPME, ddCyclNr
};
+/*! \brief Add the wallcycle count to the DD counter */
void dd_cycles_add(gmx_domdec_t *dd, float cycles, int ddCycl);
-/* Add the wallcycle count to the DD counter */
+/*! \brief Start the force flop count */
void dd_force_flop_start(gmx_domdec_t *dd, t_nrnb *nrnb);
-/* Start the force flop count */
+/*! \brief Stop the force flop count */
void dd_force_flop_stop(gmx_domdec_t *dd, t_nrnb *nrnb);
-/* Stop the force flop count */
-float dd_pme_f_ratio(gmx_domdec_t *dd);
-/* Return the PME/PP force load ratio, or -1 if nothing was measured.
+/*! \brief Return the PME/PP force load ratio, or -1 if nothing was measured.
+ *
* Should only be called on the DD master node.
*/
+float dd_pme_f_ratio(gmx_domdec_t *dd);
+/*! \brief Communicate the coordinates to the neighboring cells and do pbc. */
void dd_move_x(gmx_domdec_t *dd, matrix box, rvec x[]);
-/* Communicate the coordinates to the neighboring cells and do pbc. */
-void dd_move_f(gmx_domdec_t *dd, rvec f[], rvec *fshift);
-/* Sum the forces over the neighboring cells.
+/*! \brief Sum the forces over the neighboring cells.
+ *
* When fshift!=NULL the shift forces are updated to obtain
* the correct virial from the single sum including f.
*/
+void dd_move_f(gmx_domdec_t *dd, rvec f[], rvec *fshift);
+/*! \brief Communicate a real for each atom to the neighboring cells. */
void dd_atom_spread_real(gmx_domdec_t *dd, real v[]);
-/* Communicate a real for each atom to the neighboring cells. */
+/*! \brief Sum the contributions to a real for each atom over the neighboring cells. */
void dd_atom_sum_real(gmx_domdec_t *dd, real v[]);
-/* Sum the contributions to a real for each atom over the neighboring cells. */
+/*! \brief Partition the system over the nodes.
+ *
+ * step is only used for printing error messages.
+ * If bMasterState==TRUE then state_global from the master node is used,
+ * else state_local is redistributed between the nodes.
+ * When f!=NULL, *f will be reallocated to the size of state_local.
+ */
void dd_partition_system(FILE *fplog,
gmx_int64_t step,
t_commrec *cr,
t_nrnb *nrnb,
gmx_wallcycle_t wcycle,
gmx_bool bVerbose);
-/* Partition the system over the nodes.
- * step is only used for printing error messages.
- * If bMasterState==TRUE then state_global from the master node is used,
- * else state_local is redistributed between the nodes.
- * When f!=NULL, *f will be reallocated to the size of state_local.
- */
+/*! \brief Reset all the statistics and counters for total run counting */
void reset_dd_statistics_counters(gmx_domdec_t *dd);
-/* Reset all the statistics and counters for total run counting */
+/*! \brief Print statistics for domain decomposition communication */
void print_dd_statistics(t_commrec *cr, t_inputrec *ir, FILE *fplog);
/* In domdec_con.c */
+/*! \brief Communicates the virtual site forces, reduces the shift forces when \p fshift != NULL */
void dd_move_f_vsites(gmx_domdec_t *dd, rvec *f, rvec *fshift);
+/*! \brief Clears the forces for virtual sites */
void dd_clear_f_vsites(gmx_domdec_t *dd, rvec *f);
+/*! \brief Move x0 and also x1 if x1!=NULL. bX1IsCoord tells if to do PBC on x1 */
void dd_move_x_constraints(gmx_domdec_t *dd, matrix box,
rvec *x0, rvec *x1, gmx_bool bX1IsCoord);
-/* Move x0 and also x1 if x1!=NULL. bX1IsCoord tells if to do PBC on x1 */
+/*! \brief Communicates the coordinates involved in virtual sites */
void dd_move_x_vsites(gmx_domdec_t *dd, matrix box, rvec *x);
+/*! \brief Returns the local atom count array for all constraints
+ *
+ * The local atom count for a constraint, possible values 2/1/0, is needed
+ * to avoid not/double-counting contributions linked to the Lagrange
+ * multiplier, such as the virial and free-energy derivatives.
+ */
int *dd_constraints_nlocalatoms(gmx_domdec_t *dd);
-void dd_clear_local_constraint_indices(gmx_domdec_t *dd);
-
-void dd_clear_local_vsite_indices(gmx_domdec_t *dd);
-
-int dd_make_local_vsites(gmx_domdec_t *dd, int at_start, t_ilist *lil);
-
-int dd_make_local_constraints(gmx_domdec_t *dd, int at_start,
- const gmx_mtop_t *mtop,
- const int *cginfo,
- gmx_constr_t constr, int nrec,
- t_ilist *il_local);
-
-void init_domdec_constraints(gmx_domdec_t *dd,
- gmx_mtop_t *mtop);
-
-void init_domdec_vsites(gmx_domdec_t *dd, int n_intercg_vsite);
-
-
/* In domdec_top.c */
+/*! \brief Print error output when interactions are missing */
void dd_print_missing_interactions(FILE *fplog, t_commrec *cr,
int local_count, gmx_mtop_t *top_global, t_state *state_local);
+/*! \brief Generate and store the reverse topology */
void dd_make_reverse_top(FILE *fplog,
gmx_domdec_t *dd, gmx_mtop_t *mtop,
gmx_vsite_t *vsite,
t_inputrec *ir, gmx_bool bBCheck);
+/*! \brief Store the local charge group index in \p lcgs */
void dd_make_local_cgs(gmx_domdec_t *dd, t_block *lcgs);
+/*! \brief Generate the local topology and virtual site data */
void dd_make_local_top(gmx_domdec_t *dd, gmx_domdec_zones_t *zones,
int npbcdim, matrix box,
rvec cellsize_min, ivec npulse,
gmx_vsite_t *vsite,
gmx_mtop_t *top, gmx_localtop_t *ltop);
+/*! \brief Sort ltop->ilist when we are doing free energy. */
void dd_sort_local_top(gmx_domdec_t *dd, t_mdatoms *mdatoms,
gmx_localtop_t *ltop);
-/* Sort ltop->ilist when we are doing free energy. */
+/*! \brief Construct local topology */
gmx_localtop_t *dd_init_local_top(gmx_mtop_t *top_global);
+/*! \brief Construct local state */
void dd_init_local_state(gmx_domdec_t *dd,
t_state *state_global, t_state *local_state);
+/*! \brief Generate a list of links between charge groups that are linked by bonded interactions */
t_blocka *make_charge_group_links(gmx_mtop_t *mtop, gmx_domdec_t *dd,
cginfo_mb_t *cginfo_mb);
+/*! \brief Calculate the maximum distance involved in 2-body and multi-body bonded interactions */
void dd_bonded_cg_distance(FILE *fplog, gmx_mtop_t *mtop,
t_inputrec *ir, rvec *x, matrix box,
gmx_bool bBCheck,
real *r_2b, real *r_mb);
+/*! \brief Dump a pdb file with the current DD home + communicated atoms.
+ *
+ * When natoms=-1, dump all known atoms.
+ */
void write_dd_pdb(const char *fn, gmx_int64_t step, const char *title,
gmx_mtop_t *mtop,
t_commrec *cr,
int natoms, rvec x[], matrix box);
-/* Dump a pdb file with the current DD home + communicated atoms.
- * When natoms=-1, dump all known atoms.
- */
/* In domdec_setup.c */
+/*! \brief Returns the volume fraction of the system that is communicated */
real comm_box_frac(ivec dd_nc, real cutoff, gmx_ddbox_t *ddbox);
-/* Returns the volume fraction of the system that is communicated */
+/*! \brief Determines the optimal DD cell setup dd->nc and possibly npmenodes
+ * for the system.
+ *
+ * On the master node returns the actual cellsize limit used.
+ */
real dd_choose_grid(FILE *fplog,
t_commrec *cr, gmx_domdec_t *dd, t_inputrec *ir,
gmx_mtop_t *mtop, matrix box, gmx_ddbox_t *ddbox,
gmx_bool bDynLoadBal, real dlb_scale,
real cellsize_limit, real cutoff_dd,
gmx_bool bInterCGBondeds);
-/* Determines the optimal DD cell setup dd->nc and possibly npmenodes
- * for the system.
- * On the master node returns the actual cellsize limit used.
- */
/* In domdec_box.c */
+/*! \brief Set the box and PBC data in \p ddbox */
void set_ddbox(gmx_domdec_t *dd, gmx_bool bMasterState, t_commrec *cr_sum,
t_inputrec *ir, matrix box,
gmx_bool bCalcUnboundedSize, t_block *cgs, rvec *x,
gmx_ddbox_t *ddbox);
+/*! \brief Set the box and PBC data in \p ddbox */
void set_ddbox_cr(t_commrec *cr, ivec *dd_nc,
t_inputrec *ir, matrix box, t_block *cgs, rvec *x,
gmx_ddbox_t *ddbox);
}
#endif
-#endif /* _domdec_h */
+#endif
* the research papers on the package. Check out http://www.gromacs.org.
*/
+/*! \internal \file
+ *
+ * \brief This file defines functions used by the domdec module
+ * for (bounding) box and pbc information generation.
+ *
+ * \author Berk Hess <hess@kth.se>
+ * \ingroup module_domdec
+ */
+
#include "gmxpre.h"
#include "gromacs/domdec/domdec.h"
#include "gromacs/pbcutil/pbc.h"
#include "gromacs/utility/fatalerror.h"
+/*! \brief Calculates the average and standard deviation in 3D of n charge groups */
static void calc_cgcm_av_stddev(t_block *cgs, int n, rvec *x, rvec av, rvec stddev,
t_commrec *cr_sum)
{
}
}
+/*! \brief Determines if dimensions require triclinic treatment and stores this info in ddbox */
static void set_tric_dir(ivec *dd_nc, gmx_ddbox_t *ddbox, matrix box)
{
int npbcdim, d, i, j;
}
}
+/*! \brief This function calculates and bounding box and pbc infor and populates ddbox */
static void low_set_ddbox(t_inputrec *ir, ivec *dd_nc, matrix box,
gmx_bool bCalcUnboundedSize, int ncg, t_block *cgs, rvec *x,
t_commrec *cr_sum,
* the research papers on the package. Check out http://www.gromacs.org.
*/
+/*! \internal \file
+ *
+ * \brief This file implements functions for domdec to use
+ * while managing inter-atomic constraints.
+ *
+ * \author Berk Hess <hess@kth.se>
+ * \ingroup module_domdec
+ */
+
#include "gmxpre.h"
+#include "domdec_constraints.h"
+
#include <assert.h>
#include <algorithm>
#include "gromacs/domdec/domdec.h"
-#include "gromacs/domdec/domdec_network.h"
#include "gromacs/legacyheaders/constr.h"
#include "gromacs/legacyheaders/gmx_ga2la.h"
#include "gromacs/legacyheaders/gmx_hash.h"
#include "gromacs/utility/gmxassert.h"
#include "gromacs/utility/smalloc.h"
-typedef struct {
- int nsend;
- int *a;
- int a_nalloc;
- int nrecv;
-} gmx_specatsend_t;
-
-typedef struct {
- int *ind;
- int nalloc;
- int n;
-} ind_req_t;
-
-typedef struct gmx_domdec_specat_comm {
- /* The number of indices to receive during the setup */
- int nreq[DIM][2][2];
- /* The atoms to send */
- gmx_specatsend_t spas[DIM][2];
- gmx_bool *bSendAtom;
- int bSendAtom_nalloc;
- /* Send buffers */
- int *ibuf;
- int ibuf_nalloc;
- rvec *vbuf;
- int vbuf_nalloc;
- rvec *vbuf2;
- int vbuf2_nalloc;
- /* The range in the local buffer(s) for received atoms */
- int at_start;
- int at_end;
-
- /* The atom indices we need from the surrounding cells.
- * We can gather the indices over nthread threads.
- */
- int nthread;
- ind_req_t *ireq;
-} gmx_domdec_specat_comm_t;
+#include "domdec_specatomcomm.h"
+/*! \brief Struct used during constraint setup with domain decomposition */
typedef struct gmx_domdec_constraints {
- int *molb_con_offset;
- int *molb_ncon_mol;
- /* The fully local and connected constraints */
- int ncon;
+ //! @cond Doxygen_Suppress
+ int *molb_con_offset; /**< Offset in the constraint array for each molblock */
+ int *molb_ncon_mol; /**< The number of constraints per molecule for each molblock */
+
+ int ncon; /**< The fully local and conneced constraints */
/* The global constraint number, only required for clearing gc_req */
- int *con_gl;
- int *con_nlocat;
- int con_nalloc;
- /* Boolean that tells if a global constraint index has been requested */
- char *gc_req;
- /* Global to local communicated constraint atom only index */
- gmx_hash_t ga2la;
+ int *con_gl; /**< Global constraint indices for local constraints */
+ int *con_nlocat; /**< Number of local atoms (2/1/0) for each constraint */
+ int con_nalloc; /**< Allocation size for \p con_gl and \p con_nlocat */
+
+ char *gc_req; /**< Boolean that tells if a global constraint index has been requested; note: size global #constraints */
+ gmx_hash_t ga2la; /**< Global to local communicated constraint atom only index */
/* Multi-threading stuff */
- int nthread;
- t_ilist *ils;
+ int nthread; /**< Number of threads used for DD constraint setup */
+ t_ilist *ils; /**< Constraint ilist working arrays, size \p nthread */
+ //! @endcond
} gmx_domdec_constraints_t;
-
-static void dd_move_f_specat(gmx_domdec_t *dd, gmx_domdec_specat_comm_t *spac,
- rvec *f, rvec *fshift)
-{
- gmx_specatsend_t *spas;
- rvec *vbuf;
- int n, n0, n1, d, dim, dir, i;
- ivec vis;
- int is;
- gmx_bool bPBC, bScrew;
-
- n = spac->at_end;
- for (d = dd->ndim-1; d >= 0; d--)
- {
- dim = dd->dim[d];
- if (dd->nc[dim] > 2)
- {
- /* Pulse the grid forward and backward */
- spas = spac->spas[d];
- n0 = spas[0].nrecv;
- n1 = spas[1].nrecv;
- n -= n1 + n0;
- vbuf = spac->vbuf;
- /* Send and receive the coordinates */
- dd_sendrecv2_rvec(dd, d,
- f+n+n1, n0, vbuf, spas[0].nsend,
- f+n, n1, vbuf+spas[0].nsend, spas[1].nsend);
- for (dir = 0; dir < 2; dir++)
- {
- bPBC = ((dir == 0 && dd->ci[dim] == 0) ||
- (dir == 1 && dd->ci[dim] == dd->nc[dim]-1));
- bScrew = (bPBC && dd->bScrewPBC && dim == XX);
-
- spas = &spac->spas[d][dir];
- /* Sum the buffer into the required forces */
- if (!bPBC || (!bScrew && fshift == NULL))
- {
- for (i = 0; i < spas->nsend; i++)
- {
- rvec_inc(f[spas->a[i]], *vbuf);
- vbuf++;
- }
- }
- else
- {
- clear_ivec(vis);
- vis[dim] = (dir == 0 ? 1 : -1);
- is = IVEC2IS(vis);
- if (!bScrew)
- {
- /* Sum and add to shift forces */
- for (i = 0; i < spas->nsend; i++)
- {
- rvec_inc(f[spas->a[i]], *vbuf);
- rvec_inc(fshift[is], *vbuf);
- vbuf++;
- }
- }
- else
- {
- /* Rotate the forces */
- for (i = 0; i < spas->nsend; i++)
- {
- f[spas->a[i]][XX] += (*vbuf)[XX];
- f[spas->a[i]][YY] -= (*vbuf)[YY];
- f[spas->a[i]][ZZ] -= (*vbuf)[ZZ];
- if (fshift)
- {
- rvec_inc(fshift[is], *vbuf);
- }
- vbuf++;
- }
- }
- }
- }
- }
- else
- {
- /* Two cells, so we only need to communicate one way */
- spas = &spac->spas[d][0];
- n -= spas->nrecv;
- /* Send and receive the coordinates */
- dd_sendrecv_rvec(dd, d, dddirForward,
- f+n, spas->nrecv, spac->vbuf, spas->nsend);
- /* Sum the buffer into the required forces */
- if (dd->bScrewPBC && dim == XX &&
- (dd->ci[dim] == 0 ||
- dd->ci[dim] == dd->nc[dim]-1))
- {
- for (i = 0; i < spas->nsend; i++)
- {
- /* Rotate the force */
- f[spas->a[i]][XX] += spac->vbuf[i][XX];
- f[spas->a[i]][YY] -= spac->vbuf[i][YY];
- f[spas->a[i]][ZZ] -= spac->vbuf[i][ZZ];
- }
- }
- else
- {
- for (i = 0; i < spas->nsend; i++)
- {
- rvec_inc(f[spas->a[i]], spac->vbuf[i]);
- }
- }
- }
- }
-}
-
-void dd_move_f_vsites(gmx_domdec_t *dd, rvec *f, rvec *fshift)
-{
- if (dd->vsite_comm)
- {
- dd_move_f_specat(dd, dd->vsite_comm, f, fshift);
- }
-}
-
-void dd_clear_f_vsites(gmx_domdec_t *dd, rvec *f)
-{
- int i;
-
- if (dd->vsite_comm)
- {
- for (i = dd->vsite_comm->at_start; i < dd->vsite_comm->at_end; i++)
- {
- clear_rvec(f[i]);
- }
- }
-}
-
-static void dd_move_x_specat(gmx_domdec_t *dd, gmx_domdec_specat_comm_t *spac,
- matrix box,
- rvec *x0,
- rvec *x1, gmx_bool bX1IsCoord)
-{
- gmx_specatsend_t *spas;
- rvec *x, *vbuf, *rbuf;
- int nvec, v, n, nn, ns0, ns1, nr0, nr1, nr, d, dim, dir, i;
- gmx_bool bPBC, bScrew = FALSE;
- rvec shift = {0, 0, 0};
-
- nvec = 1;
- if (x1 != NULL)
- {
- nvec++;
- }
-
- n = spac->at_start;
- for (d = 0; d < dd->ndim; d++)
- {
- dim = dd->dim[d];
- if (dd->nc[dim] > 2)
- {
- /* Pulse the grid forward and backward */
- vbuf = spac->vbuf;
- for (dir = 0; dir < 2; dir++)
- {
- if (dir == 0 && dd->ci[dim] == 0)
- {
- bPBC = TRUE;
- bScrew = (dd->bScrewPBC && dim == XX);
- copy_rvec(box[dim], shift);
- }
- else if (dir == 1 && dd->ci[dim] == dd->nc[dim]-1)
- {
- bPBC = TRUE;
- bScrew = (dd->bScrewPBC && dim == XX);
- for (i = 0; i < DIM; i++)
- {
- shift[i] = -box[dim][i];
- }
- }
- else
- {
- bPBC = FALSE;
- bScrew = FALSE;
- }
- spas = &spac->spas[d][dir];
- for (v = 0; v < nvec; v++)
- {
- x = (v == 0 ? x0 : x1);
- /* Copy the required coordinates to the send buffer */
- if (!bPBC || (v == 1 && !bX1IsCoord))
- {
- /* Only copy */
- for (i = 0; i < spas->nsend; i++)
- {
- copy_rvec(x[spas->a[i]], *vbuf);
- vbuf++;
- }
- }
- else if (!bScrew)
- {
- /* Shift coordinates */
- for (i = 0; i < spas->nsend; i++)
- {
- rvec_add(x[spas->a[i]], shift, *vbuf);
- vbuf++;
- }
- }
- else
- {
- /* Shift and rotate coordinates */
- for (i = 0; i < spas->nsend; i++)
- {
- (*vbuf)[XX] = x[spas->a[i]][XX] + shift[XX];
- (*vbuf)[YY] = box[YY][YY] - x[spas->a[i]][YY] + shift[YY];
- (*vbuf)[ZZ] = box[ZZ][ZZ] - x[spas->a[i]][ZZ] + shift[ZZ];
- vbuf++;
- }
- }
- }
- }
- /* Send and receive the coordinates */
- spas = spac->spas[d];
- ns0 = spas[0].nsend;
- nr0 = spas[0].nrecv;
- ns1 = spas[1].nsend;
- nr1 = spas[1].nrecv;
- if (nvec == 1)
- {
- dd_sendrecv2_rvec(dd, d,
- spac->vbuf+ns0, ns1, x0+n, nr1,
- spac->vbuf, ns0, x0+n+nr1, nr0);
- }
- else
- {
- /* Communicate both vectors in one buffer */
- rbuf = spac->vbuf2;
- dd_sendrecv2_rvec(dd, d,
- spac->vbuf+2*ns0, 2*ns1, rbuf, 2*nr1,
- spac->vbuf, 2*ns0, rbuf+2*nr1, 2*nr0);
- /* Split the buffer into the two vectors */
- nn = n;
- for (dir = 1; dir >= 0; dir--)
- {
- nr = spas[dir].nrecv;
- for (v = 0; v < 2; v++)
- {
- x = (v == 0 ? x0 : x1);
- for (i = 0; i < nr; i++)
- {
- copy_rvec(*rbuf, x[nn+i]);
- rbuf++;
- }
- }
- nn += nr;
- }
- }
- n += nr0 + nr1;
- }
- else
- {
- spas = &spac->spas[d][0];
- /* Copy the required coordinates to the send buffer */
- vbuf = spac->vbuf;
- for (v = 0; v < nvec; v++)
- {
- x = (v == 0 ? x0 : x1);
- if (dd->bScrewPBC && dim == XX &&
- (dd->ci[XX] == 0 || dd->ci[XX] == dd->nc[XX]-1))
- {
- /* Here we only perform the rotation, the rest of the pbc
- * is handled in the constraint or viste routines.
- */
- for (i = 0; i < spas->nsend; i++)
- {
- (*vbuf)[XX] = x[spas->a[i]][XX];
- (*vbuf)[YY] = box[YY][YY] - x[spas->a[i]][YY];
- (*vbuf)[ZZ] = box[ZZ][ZZ] - x[spas->a[i]][ZZ];
- vbuf++;
- }
- }
- else
- {
- for (i = 0; i < spas->nsend; i++)
- {
- copy_rvec(x[spas->a[i]], *vbuf);
- vbuf++;
- }
- }
- }
- /* Send and receive the coordinates */
- if (nvec == 1)
- {
- dd_sendrecv_rvec(dd, d, dddirBackward,
- spac->vbuf, spas->nsend, x0+n, spas->nrecv);
- }
- else
- {
- /* Communicate both vectors in one buffer */
- rbuf = spac->vbuf2;
- dd_sendrecv_rvec(dd, d, dddirBackward,
- spac->vbuf, 2*spas->nsend, rbuf, 2*spas->nrecv);
- /* Split the buffer into the two vectors */
- nr = spas[0].nrecv;
- for (v = 0; v < 2; v++)
- {
- x = (v == 0 ? x0 : x1);
- for (i = 0; i < nr; i++)
- {
- copy_rvec(*rbuf, x[n+i]);
- rbuf++;
- }
- }
- }
- n += spas->nrecv;
- }
- }
-}
-
void dd_move_x_constraints(gmx_domdec_t *dd, matrix box,
rvec *x0, rvec *x1, gmx_bool bX1IsCoord)
{
}
}
-void dd_move_x_vsites(gmx_domdec_t *dd, matrix box, rvec *x)
-{
- if (dd->vsite_comm)
- {
- dd_move_x_specat(dd, dd->vsite_comm, box, x, NULL, FALSE);
- }
-}
-
int *dd_constraints_nlocalatoms(gmx_domdec_t *dd)
{
if (dd->constraints)
}
}
-static int setup_specat_communication(gmx_domdec_t *dd,
- ind_req_t *ireq,
- gmx_domdec_specat_comm_t *spac,
- gmx_hash_t ga2la_specat,
- int at_start,
- int vbuf_fac,
- const char *specat_type,
- const char *add_err)
-{
- int nsend[2], nlast, nsend_zero[2] = {0, 0}, *nsend_ptr;
- int d, dim, ndir, dir, nr, ns, i, nrecv_local, n0, start, indr, ind, buf[2];
- int nat_tot_specat, nat_tot_prev, nalloc_old;
- gmx_bool bPBC;
- gmx_specatsend_t *spas;
-
- if (debug)
- {
- fprintf(debug, "Begin setup_specat_communication for %s\n", specat_type);
- }
-
- /* nsend[0]: the number of atoms requested by this node only,
- * we communicate this for more efficients checks
- * nsend[1]: the total number of requested atoms
- */
- nsend[0] = ireq->n;
- nsend[1] = nsend[0];
- nlast = nsend[1];
- for (d = dd->ndim-1; d >= 0; d--)
- {
- /* Pulse the grid forward and backward */
- dim = dd->dim[d];
- bPBC = (dim < dd->npbcdim);
- if (dd->nc[dim] == 2)
- {
- /* Only 2 cells, so we only need to communicate once */
- ndir = 1;
- }
- else
- {
- ndir = 2;
- }
- for (dir = 0; dir < ndir; dir++)
- {
- if (!bPBC &&
- dd->nc[dim] > 2 &&
- ((dir == 0 && dd->ci[dim] == dd->nc[dim] - 1) ||
- (dir == 1 && dd->ci[dim] == 0)))
- {
- /* No pbc: the fist/last cell should not request atoms */
- nsend_ptr = nsend_zero;
- }
- else
- {
- nsend_ptr = nsend;
- }
- /* Communicate the number of indices */
- dd_sendrecv_int(dd, d, dir == 0 ? dddirForward : dddirBackward,
- nsend_ptr, 2, spac->nreq[d][dir], 2);
- nr = spac->nreq[d][dir][1];
- if (nlast+nr > ireq->nalloc)
- {
- ireq->nalloc = over_alloc_dd(nlast+nr);
- srenew(ireq->ind, ireq->nalloc);
- }
- /* Communicate the indices */
- dd_sendrecv_int(dd, d, dir == 0 ? dddirForward : dddirBackward,
- ireq->ind, nsend_ptr[1], ireq->ind+nlast, nr);
- nlast += nr;
- }
- nsend[1] = nlast;
- }
- if (debug)
- {
- fprintf(debug, "Communicated the counts\n");
- }
-
- /* Search for the requested atoms and communicate the indices we have */
- nat_tot_specat = at_start;
- nrecv_local = 0;
- for (d = 0; d < dd->ndim; d++)
- {
- /* Pulse the grid forward and backward */
- if (dd->dim[d] >= dd->npbcdim || dd->nc[dd->dim[d]] > 2)
- {
- ndir = 2;
- }
- else
- {
- ndir = 1;
- }
- nat_tot_prev = nat_tot_specat;
- for (dir = ndir-1; dir >= 0; dir--)
- {
- if (nat_tot_specat > spac->bSendAtom_nalloc)
- {
- nalloc_old = spac->bSendAtom_nalloc;
- spac->bSendAtom_nalloc = over_alloc_dd(nat_tot_specat);
- srenew(spac->bSendAtom, spac->bSendAtom_nalloc);
- for (i = nalloc_old; i < spac->bSendAtom_nalloc; i++)
- {
- spac->bSendAtom[i] = FALSE;
- }
- }
- spas = &spac->spas[d][dir];
- n0 = spac->nreq[d][dir][0];
- nr = spac->nreq[d][dir][1];
- if (debug)
- {
- fprintf(debug, "dim=%d, dir=%d, searching for %d atoms\n",
- d, dir, nr);
- }
- start = nlast - nr;
- spas->nsend = 0;
- nsend[0] = 0;
- for (i = 0; i < nr; i++)
- {
- indr = ireq->ind[start+i];
- ind = -1;
- /* Check if this is a home atom and if so ind will be set */
- if (!ga2la_get_home(dd->ga2la, indr, &ind))
- {
- /* Search in the communicated atoms */
- ind = gmx_hash_get_minone(ga2la_specat, indr);
- }
- if (ind >= 0)
- {
- if (i < n0 || !spac->bSendAtom[ind])
- {
- if (spas->nsend+1 > spas->a_nalloc)
- {
- spas->a_nalloc = over_alloc_large(spas->nsend+1);
- srenew(spas->a, spas->a_nalloc);
- }
- /* Store the local index so we know which coordinates
- * to send out later.
- */
- spas->a[spas->nsend] = ind;
- spac->bSendAtom[ind] = TRUE;
- if (spas->nsend+1 > spac->ibuf_nalloc)
- {
- spac->ibuf_nalloc = over_alloc_large(spas->nsend+1);
- srenew(spac->ibuf, spac->ibuf_nalloc);
- }
- /* Store the global index so we can send it now */
- spac->ibuf[spas->nsend] = indr;
- if (i < n0)
- {
- nsend[0]++;
- }
- spas->nsend++;
- }
- }
- }
- nlast = start;
- /* Clear the local flags */
- for (i = 0; i < spas->nsend; i++)
- {
- spac->bSendAtom[spas->a[i]] = FALSE;
- }
- /* Send and receive the number of indices to communicate */
- nsend[1] = spas->nsend;
- dd_sendrecv_int(dd, d, dir == 0 ? dddirBackward : dddirForward,
- nsend, 2, buf, 2);
- if (debug)
- {
- fprintf(debug, "Send to rank %d, %d (%d) indices, "
- "receive from rank %d, %d (%d) indices\n",
- dd->neighbor[d][1-dir], nsend[1], nsend[0],
- dd->neighbor[d][dir], buf[1], buf[0]);
- if (gmx_debug_at)
- {
- for (i = 0; i < spas->nsend; i++)
- {
- fprintf(debug, " %d", spac->ibuf[i]+1);
- }
- fprintf(debug, "\n");
- }
- }
- nrecv_local += buf[0];
- spas->nrecv = buf[1];
- if (nat_tot_specat + spas->nrecv > dd->gatindex_nalloc)
- {
- dd->gatindex_nalloc =
- over_alloc_dd(nat_tot_specat + spas->nrecv);
- srenew(dd->gatindex, dd->gatindex_nalloc);
- }
- /* Send and receive the indices */
- dd_sendrecv_int(dd, d, dir == 0 ? dddirBackward : dddirForward,
- spac->ibuf, spas->nsend,
- dd->gatindex+nat_tot_specat, spas->nrecv);
- nat_tot_specat += spas->nrecv;
- }
-
- /* Allocate the x/f communication buffers */
- ns = spac->spas[d][0].nsend;
- nr = spac->spas[d][0].nrecv;
- if (ndir == 2)
- {
- ns += spac->spas[d][1].nsend;
- nr += spac->spas[d][1].nrecv;
- }
- if (vbuf_fac*ns > spac->vbuf_nalloc)
- {
- spac->vbuf_nalloc = over_alloc_dd(vbuf_fac*ns);
- srenew(spac->vbuf, spac->vbuf_nalloc);
- }
- if (vbuf_fac == 2 && vbuf_fac*nr > spac->vbuf2_nalloc)
- {
- spac->vbuf2_nalloc = over_alloc_dd(vbuf_fac*nr);
- srenew(spac->vbuf2, spac->vbuf2_nalloc);
- }
-
- /* Make a global to local index for the communication atoms */
- for (i = nat_tot_prev; i < nat_tot_specat; i++)
- {
- gmx_hash_change_or_set(ga2la_specat, dd->gatindex[i], i);
- }
- }
-
- /* Check that in the end we got the number of atoms we asked for */
- if (nrecv_local != ireq->n)
- {
- if (debug)
- {
- fprintf(debug, "Requested %d, received %d (tot recv %d)\n",
- ireq->n, nrecv_local, nat_tot_specat-at_start);
- if (gmx_debug_at)
- {
- for (i = 0; i < ireq->n; i++)
- {
- ind = gmx_hash_get_minone(ga2la_specat, ireq->ind[i]);
- fprintf(debug, " %s%d",
- (ind >= 0) ? "" : "!",
- ireq->ind[i]+1);
- }
- fprintf(debug, "\n");
- }
- }
- fprintf(stderr, "\nDD cell %d %d %d: Neighboring cells do not have atoms:",
- dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
- for (i = 0; i < ireq->n; i++)
- {
- if (gmx_hash_get_minone(ga2la_specat, ireq->ind[i]) < 0)
- {
- fprintf(stderr, " %d", ireq->ind[i]+1);
- }
- }
- fprintf(stderr, "\n");
- gmx_fatal(FARGS, "DD cell %d %d %d could only obtain %d of the %d atoms that are connected via %ss from the neighboring cells. This probably means your %s lengths are too long compared to the domain decomposition cell size. Decrease the number of domain decomposition grid cells%s%s.",
- dd->ci[XX], dd->ci[YY], dd->ci[ZZ],
- nrecv_local, ireq->n, specat_type,
- specat_type, add_err,
- dd->bGridJump ? " or use the -rcon option of mdrun" : "");
- }
-
- spac->at_start = at_start;
- spac->at_end = nat_tot_specat;
-
- if (debug)
- {
- fprintf(debug, "Done setup_specat_communication\n");
- }
-
- return nat_tot_specat;
-}
-
+/*! \brief Walks over the constraints out from the local atoms into the non-local atoms and adds them to a list */
static void walk_out(int con, int con_offset, int a, int offset, int nrec,
int ncon1, const t_iatom *ia1, const t_iatom *ia2,
const t_blocka *at2con,
}
}
+/*! \brief Looks up SETTLE constraints for a range of charge-groups */
static void atoms_to_settles(gmx_domdec_t *dd,
const gmx_mtop_t *mtop,
const int *cginfo,
gmx_mtop_atomlookup_destroy(alook);
}
+/*! \brief Looks up constraint for the local atoms */
static void atoms_to_constraints(gmx_domdec_t *dd,
const gmx_mtop_t *mtop,
const int *cginfo,
return at_end;
}
-int dd_make_local_vsites(gmx_domdec_t *dd, int at_start, t_ilist *lil)
-{
- gmx_domdec_specat_comm_t *spac;
- ind_req_t *ireq;
- gmx_hash_t ga2la_specat;
- int ftype, nral, i, j, a;
- t_ilist *lilf;
- t_iatom *iatoms;
- int at_end;
-
- spac = dd->vsite_comm;
- ireq = &spac->ireq[0];
- ga2la_specat = dd->ga2la_vsite;
-
- ireq->n = 0;
- /* Loop over all the home vsites */
- for (ftype = 0; ftype < F_NRE; ftype++)
- {
- if (interaction_function[ftype].flags & IF_VSITE)
- {
- nral = NRAL(ftype);
- lilf = &lil[ftype];
- for (i = 0; i < lilf->nr; i += 1+nral)
- {
- iatoms = lilf->iatoms + i;
- /* Check if we have the other atoms */
- for (j = 1; j < 1+nral; j++)
- {
- if (iatoms[j] < 0)
- {
- /* This is not a home atom,
- * we need to ask our neighbors.
- */
- a = -iatoms[j] - 1;
- /* Check to not ask for the same atom more than once */
- if (gmx_hash_get_minone(dd->ga2la_vsite, a) == -1)
- {
- /* Add this non-home atom to the list */
- if (ireq->n+1 > ireq->nalloc)
- {
- ireq->nalloc = over_alloc_large(ireq->n+1);
- srenew(ireq->ind, ireq->nalloc);
- }
- ireq->ind[ireq->n++] = a;
- /* Temporarily mark with -2,
- * we get the index later.
- */
- gmx_hash_set(ga2la_specat, a, -2);
- }
- }
- }
- }
- }
- }
-
- at_end = setup_specat_communication(dd, ireq, dd->vsite_comm, ga2la_specat,
- at_start, 1, "vsite", "");
-
- /* Fill in the missing indices */
- for (ftype = 0; ftype < F_NRE; ftype++)
- {
- if (interaction_function[ftype].flags & IF_VSITE)
- {
- nral = NRAL(ftype);
- lilf = &lil[ftype];
- for (i = 0; i < lilf->nr; i += 1+nral)
- {
- iatoms = lilf->iatoms + i;
- for (j = 1; j < 1+nral; j++)
- {
- if (iatoms[j] < 0)
- {
- iatoms[j] = gmx_hash_get_minone(ga2la_specat, -iatoms[j]-1);
- }
- }
- }
- }
- }
-
- return at_end;
-}
-
-static gmx_domdec_specat_comm_t *specat_comm_init(int nthread)
-{
- gmx_domdec_specat_comm_t *spac;
-
- snew(spac, 1);
- spac->nthread = nthread;
- snew(spac->ireq, spac->nthread);
-
- return spac;
-}
-
void init_domdec_constraints(gmx_domdec_t *dd,
gmx_mtop_t *mtop)
{
dd->constraint_comm = specat_comm_init(dc->nthread);
}
-
-void init_domdec_vsites(gmx_domdec_t *dd, int n_intercg_vsite)
-{
- if (debug)
- {
- fprintf(debug, "Begin init_domdec_vsites\n");
- }
-
- /* Use a hash table for the global to local index.
- * The number of keys is a rough estimate, it will be optimized later.
- */
- dd->ga2la_vsite = gmx_hash_init(std::min(n_intercg_vsite/20,
- n_intercg_vsite/(2*dd->nnodes)));
-
- dd->vsite_comm = specat_comm_init(1);
-}
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2005,2006,2007,2008,2009,2010,2012,2013,2014, 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.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ *
+ * \brief This file declares functions for domdec to use
+ * while managing inter-atomic constraints.
+ *
+ * \author Berk Hess <hess@kth.se>
+ * \ingroup module_domdec
+ */
+
+#ifndef GMX_DOMDEC_DOMDEC_CONSTRAINTS_H
+#define GMX_DOMDEC_DOMDEC_CONSTRAINTS_H
+
+#include "gromacs/legacyheaders/types/commrec_fwd.h"
+#include "gromacs/legacyheaders/types/constr.h"
+#include "gromacs/topology/idef.h"
+#include "gromacs/topology/topology.h"
+
+/*! \brief Clears the local indices for the constraint communication setup */
+void dd_clear_local_constraint_indices(gmx_domdec_t *dd);
+
+/*! \brief Sets up communication and atom indices for all local+connected constraints */
+int dd_make_local_constraints(gmx_domdec_t *dd, int at_start,
+ const gmx_mtop_t *mtop,
+ const int *cginfo,
+ gmx_constr_t constr, int nrec,
+ t_ilist *il_local);
+
+/*! \brief Initializes the data structures for constraint communication */
+void init_domdec_constraints(gmx_domdec_t *dd,
+ gmx_mtop_t *mtop);
+
+#endif
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 1991-2003 David van der Spoel, Erik Lindahl, University of Groningen.
+ * Copyright (c) 2013,2014, 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.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ *
+ * \brief Declares implementation functions and types for the domain
+ * decomposition module.
+ *
+ * \author Berk Hess <hess@kth.se>
+ * \ingroup module_domdec
+ */
+#ifndef GMX_DOMDEC_DOMDEC_INTERNAL_H
+#define GMX_DOMDEC_DOMDEC_INTERNAL_H
+
+#include "gromacs/legacyheaders/types/commrec.h"
+
+/*! \brief Returns the DD cut-off distance for multi-body interactions */
+real dd_cutoff_multibody(const gmx_domdec_t *dd);
+
+/*! \brief Returns the DD cut-off distance for two-body interactions */
+real dd_cutoff_twobody(const gmx_domdec_t *dd);
+
+#endif
* the research papers on the package. Check out http://www.gromacs.org.
*/
+/*! \internal \file
+ *
+ * \brief This file defines functions for (mostly) the domdec module
+ * to use MPI functionality
+ *
+ * \author Berk Hess <hess@kth.se>
+ * \ingroup module_domdec
+ */
+
#include "gmxpre.h"
#include "domdec_network.h"
#include "gromacs/utility/gmxmpi.h"
+/*! \brief Returns the MPI rank of the domain decomposition master rank */
#define DDMASTERRANK(dd) (dd->masterrank)
* the research papers on the package. Check out http://www.gromacs.org.
*/
-#ifndef _domdec_network_h
-#define _domdec_network_h
+/*! \libinternal \file
+ *
+ * \brief This file declares functions for (mostly) the domdec module
+ * to use MPI functionality
+ *
+ * \todo Wrap the raw dd_bcast in md.cpp into a higher-level function
+ * in the domdec module, then this file can be module-internal.
+ *
+ * \author Berk Hess <hess@kth.se>
+ * \inlibraryapi
+ * \ingroup module_domdec
+ */
+#ifndef GMX_DOMDEC_DOMDEC_NETWORK_H
+#define GMX_DOMDEC_DOMDEC_NETWORK_H
#include "gromacs/legacyheaders/typedefs.h"
-#ifdef __cplusplus
-extern "C" {
-#endif
-
+/* \brief */
enum {
dddirForward, dddirBackward
};
-/* Move integers in the comm. region one cell along the domain decomposition
- * in the dimension indexed by ddimind
- * forward (direction=dddirFoward) or backward (direction=dddirBackward).
+/*! \brief Move integers in the communication region one cell along
+ * the domain decomposition
+ *
+ * Moves in the dimension indexed by ddimind, either forward
+ * (direction=dddirFoward) or backward (direction=dddirBackward).
*/
void
dd_sendrecv_int(const gmx_domdec_t *dd,
int *buf_s, int n_s,
int *buf_r, int n_r);
-/* Move reals in the comm. region one cell along the domain decomposition
- * in the dimension indexed by ddimind
- * forward (direction=dddirFoward) or backward (direction=dddirBackward).
+/*! \brief Move reals in the comm. region one cell along the domain decomposition
+ *
+ * Moves in the dimension indexed by ddimind, either forward
+ * (direction=dddirFoward) or backward (direction=dddirBackward).
*/
void
dd_sendrecv_real(const gmx_domdec_t *dd,
real *buf_s, int n_s,
real *buf_r, int n_r);
-/* Move revc's in the comm. region one cell along the domain decomposition
- * in dimension indexed by ddimind
- * forward (direction=dddirFoward) or backward (direction=dddirBackward).
+/*! \brief Move revc's in the comm. region one cell along the domain decomposition
+ *
+ * Moves in dimension indexed by ddimind, either forward
+ * (direction=dddirFoward) or backward (direction=dddirBackward).
*/
void
dd_sendrecv_rvec(const gmx_domdec_t *dd,
rvec *buf_r, int n_r);
-/* Move revc's in the comm. region one cell along the domain decomposition
- * in dimension indexed by ddimind
- * simultaneously in the forward and backward directions.
+/*! \brief Move revc's in the comm. region one cell along the domain decomposition
+ *
+ * Moves in dimension indexed by ddimind, simultaneously in the forward
+ * and backward directions.
*/
void
dd_sendrecv2_rvec(const gmx_domdec_t *dd,
* The DD master node is the master for these operations.
*/
+/*! \brief Broadcasts \p nbytes from \p data on \p DDMASTERRANK to all PP ranks */
void
dd_bcast(gmx_domdec_t *dd, int nbytes, void *data);
-/* Copies src to dest on the master node and then broadcasts */
+/*! \brief Copies \p nbytes from \p src to \p dest on \p DDMASTERRANK
+ * and then broadcasts to \p dest on all PP ranks */
void
dd_bcastc(gmx_domdec_t *dd, int nbytes, void *src, void *dest);
+/*! \brief Scatters \p nbytes from \p src on \p DDMASTERRANK to all PP ranks, received in \p dest */
void
dd_scatter(gmx_domdec_t *dd, int nbytes, void *src, void *dest);
+/*! \brief Gathers \p nbytes from \p src on all PP ranks, received in \p dest on \p DDMASTERRANK */
void
dd_gather(gmx_domdec_t *dd, int nbytes, void *src, void *dest);
-/* If rcount==0, rbuf is allowed to be NULL */
+/*! \brief Scatters \p scounts bytes from \p src on \p DDMASTERRANK to all PP ranks, receiving \p rcount bytes in \p dest.
+ *
+ * See man MPI_Scatterv for details of how to construct scounts and disps.
+ * If rcount==0, rbuf is allowed to be NULL */
void
dd_scatterv(gmx_domdec_t *dd,
int *scounts, int *disps, void *sbuf,
int rcount, void *rbuf);
-/* If scount==0, sbuf is allowed to be NULL */
+/*! \brief Gathers \p rcount bytes from \p src on all PP ranks, received in \p scounts bytes in \p dest on \p DDMASTERRANK.
+ *
+ * See man MPI_Gatherv for details of how to construct scounts and disps.
+ *
+ * If scount==0, sbuf is allowed to be NULL */
void
dd_gatherv(gmx_domdec_t *dd,
int scount, void *sbuf,
int *rcounts, int *disps, void *rbuf);
-#ifdef __cplusplus
-}
#endif
-
-#endif /* _domdec_network_h */
* the research papers on the package. Check out http://www.gromacs.org.
*/
+/*! \internal \file
+ *
+ * \brief This file defines functions used by the domdec module
+ * in its initial setup phase.
+ *
+ * \author Berk Hess <hess@kth.se>
+ * \ingroup module_domdec
+ */
+
#include "gmxpre.h"
#include <assert.h>
#include "gromacs/math/vec.h"
#include "gromacs/utility/smalloc.h"
-/* Margin for setting up the DD grid */
+/*! \brief Margin for setting up the DD grid */
#define DD_GRID_MARGIN_PRES_SCALE 1.05
+/*! \brief Factorize \p n.
+ *
+ * \param[in] n Value to factorize
+ * \param[out] fac Pointer to array of factors (to be allocated in this function)
+ * \param[out] mfac Pointer to array of the number of times each factor repeats in the factorization (to be allocated in this function)
+ * \return The number of unique factors
+ */
static int factorize(int n, int **fac, int **mfac)
{
int d, ndiv;
return ndiv;
}
+/*! \brief Find largest divisor of \p n smaller than \p n*/
static gmx_bool largest_divisor(int n)
{
int ndiv, *div, *mdiv, ldiv;
return ldiv;
}
+/*! \brief Compute largest common divisor of \p n1 and \b n2 */
static int lcd(int n1, int n2)
{
int d, i;
return d;
}
-static gmx_bool fits_pme_ratio(int nnodes, int npme, float ratio)
+/*! \brief Returns TRUE when there are enough PME ranks for the ratio */
+static gmx_bool fits_pme_ratio(int nrank_tot, int nrank_pme, float ratio)
{
- return ((double)npme/(double)nnodes > 0.95*ratio);
+ return ((double)nrank_pme/(double)nrank_tot > 0.95*ratio);
}
-static gmx_bool fits_pp_pme_perf(int nnodes, int npme, float ratio)
+/*! \brief Returns TRUE when npme out of ntot ranks doing PME is expected to give reasonable performance */
+static gmx_bool fits_pp_pme_perf(int ntot, int npme, float ratio)
{
int ndiv, *div, *mdiv, ldiv;
int npp_root3, npme_root2;
- ndiv = factorize(nnodes-npme, &div, &mdiv);
+ ndiv = factorize(ntot - npme, &div, &mdiv);
ldiv = div[ndiv-1];
sfree(div);
sfree(mdiv);
- npp_root3 = static_cast<int>(std::pow(nnodes-npme, 1.0/3.0) + 0.5);
+ npp_root3 = static_cast<int>(std::pow(ntot - npme, 1.0/3.0) + 0.5);
npme_root2 = static_cast<int>(std::sqrt(static_cast<double>(npme)) + 0.5);
/* The check below gives a reasonable division:
- * factor 5 allowed at 5 or more PP nodes,
- * factor 7 allowed at 49 or more PP nodes.
+ * factor 5 allowed at 5 or more PP ranks,
+ * factor 7 allowed at 49 or more PP ranks.
*/
if (ldiv > 3 + npp_root3)
{
return FALSE;
}
- /* Check if the number of PP and PME nodes have a reasonable sized
+ /* Check if the number of PP and PME ranks have a reasonable sized
* denominator in common, such that we can use 2D PME decomposition
* when required (which requires nx_pp == nx_pme).
* The factor of 2 allows for a maximum ratio of 2^2=4
* between nx_pme and ny_pme.
*/
- if (lcd(nnodes-npme, npme)*2 < npme_root2)
+ if (lcd(ntot - npme, npme)*2 < npme_root2)
{
return FALSE;
}
/* Does this division gives a reasonable PME load? */
- return fits_pme_ratio(nnodes, npme, ratio);
+ return fits_pme_ratio(ntot, npme, ratio);
}
+/*! \brief Make a guess for the number of PME ranks to use. */
static int guess_npme(FILE *fplog, gmx_mtop_t *mtop, t_inputrec *ir, matrix box,
- int nnodes)
+ int nrank_tot)
{
float ratio;
int npme;
fprintf(fplog, "Guess for relative PME load: %.2f\n", ratio);
}
- /* We assume the optimal node ratio is close to the load ratio.
+ /* We assume the optimal rank ratio is close to the load ratio.
* The communication load is neglected,
* but (hopefully) this will balance out between PP and PME.
*/
- if (!fits_pme_ratio(nnodes, nnodes/2, ratio))
+ if (!fits_pme_ratio(nrank_tot, nrank_tot/2, ratio))
{
- /* We would need more than nnodes/2 PME only nodes,
+ /* We would need more than nrank_tot/2 PME only nodes,
* which is not possible. Since the PME load is very high,
- * we will not loose much performance when all nodes do PME.
+ * we will not loose much performance when all ranks do PME.
*/
return 0;
}
- /* First try to find npme as a factor of nnodes up to nnodes/3.
+ /* First try to find npme as a factor of nrank_tot up to nrank_tot/3.
* We start with a minimum PME node fraction of 1/16
* and avoid ratios which lead to large prime factors in nnodes-npme.
*/
- npme = (nnodes + 15)/16;
- while (npme <= nnodes/3)
+ npme = (nrank_tot + 15)/16;
+ while (npme <= nrank_tot/3)
{
- if (nnodes % npme == 0)
+ if (nrank_tot % npme == 0)
{
/* Note that fits_perf might change the PME grid,
* in the current implementation it does not.
*/
- if (fits_pp_pme_perf(nnodes, npme, ratio))
+ if (fits_pp_pme_perf(nrank_tot, npme, ratio))
{
break;
}
}
npme++;
}
- if (npme > nnodes/3)
+ if (npme > nrank_tot/3)
{
/* Try any possible number for npme */
npme = 1;
- while (npme <= nnodes/2)
+ while (npme <= nrank_tot/2)
{
/* Note that fits_perf may change the PME grid */
- if (fits_pp_pme_perf(nnodes, npme, ratio))
+ if (fits_pp_pme_perf(nrank_tot, npme, ratio))
{
break;
}
npme++;
}
}
- if (npme > nnodes/2)
+ if (npme > nrank_tot/2)
{
gmx_fatal(FARGS, "Could not find an appropriate number of separate PME ranks. i.e. >= %5f*#ranks (%d) and <= #ranks/2 (%d) and reasonable performance wise (grid_x=%d, grid_y=%d).\n"
"Use the -npme option of mdrun or change the number of ranks or the PME grid dimensions, see the manual for details.",
- ratio, (int)(0.95*ratio*nnodes+0.5), nnodes/2, ir->nkx, ir->nky);
+ ratio, (int)(0.95*ratio*nrank_tot + 0.5), nrank_tot/2, ir->nkx, ir->nky);
/* Keep the compiler happy */
npme = 0;
}
fprintf(fplog,
"Will use %d particle-particle and %d PME only ranks\n"
"This is a guess, check the performance at the end of the log file\n",
- nnodes-npme, npme);
+ nrank_tot - npme, npme);
}
fprintf(stderr, "\n"
"Will use %d particle-particle and %d PME only ranks\n"
"This is a guess, check the performance at the end of the log file\n",
- nnodes-npme, npme);
+ nrank_tot - npme, npme);
}
return npme;
}
+/*! \brief Return \p n divided by \p f rounded up to the next integer. */
static int div_up(int n, int f)
{
return (n + f - 1)/f;
return comm_vol;
}
+/*! \brief Return whether the DD inhomogeneous in the z direction */
static gmx_bool inhomogeneous_z(const t_inputrec *ir)
{
return ((EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD) &&
ir->ePBC == epbcXYZ && ir->ewald_geometry == eewg3DC);
}
-/* Avoid integer overflows */
+/*! \brief Estimate cost of PME FFT communication
+ *
+ * This only takes the communication into account and not imbalance
+ * in the calculation. But the imbalance in communication and calculation
+ * are similar and therefore these formulas also prefer load balance
+ * in the FFT and pme_solve calculation.
+ */
static float comm_pme_cost_vol(int npme, int a, int b, int c)
{
+ /* We use a float here, since an integer might overflow */
float comm_vol;
comm_vol = npme - 1;
comm_vol *= div_up(a, npme);
comm_vol *= div_up(b, npme);
comm_vol *= c;
+
return comm_vol;
}
+/*! \brief Estimate cost of communication for a possible domain decomposition. */
static float comm_cost_est(real limit, real cutoff,
matrix box, gmx_ddbox_t *ddbox,
int natoms, t_inputrec *ir,
}
}
- /* PME FFT communication volume.
- * This only takes the communication into account and not imbalance
- * in the calculation. But the imbalance in communication and calculation
- * are similar and therefore these formulas also prefer load balance
- * in the FFT and pme_solve calculation.
- */
comm_pme += comm_pme_cost_vol(npme[YY], ir->nky, ir->nkz, ir->nkx);
comm_pme += comm_pme_cost_vol(npme[XX], ir->nkx, ir->nky, ir->nkz);
return 3*natoms*(comm_vol + cost_pbcdx) + comm_pme;
}
+/*! \brief Assign penalty factors to possible domain decompositions, based on the estimated communication costs. */
static void assign_factors(gmx_domdec_t *dd,
real limit, real cutoff,
matrix box, gmx_ddbox_t *ddbox,
}
}
+/*! \brief Determine the optimal distribution of DD cells for the simulation system and number of MPI ranks */
static real optimize_ncells(FILE *fplog,
int nnodes_tot, int npme_only,
gmx_bool bDynLoadBal, real dlb_scale,
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2006,2007,2008,2009,2010,2012,2013,2014, 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.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+
+/*! \internal \file
+ *
+ * \brief This file implements functions for domdec to use
+ * while managing inter-atomic constraints.
+ *
+ * \author Berk Hess <hess@kth.se>
+ * \ingroup module_domdec
+ */
+
+#include "gmxpre.h"
+
+#include "domdec_specatomcomm.h"
+
+#include <assert.h>
+
+#include <algorithm>
+
+#include "gromacs/domdec/domdec.h"
+#include "gromacs/domdec/domdec_network.h"
+#include "gromacs/legacyheaders/gmx_ga2la.h"
+#include "gromacs/legacyheaders/gmx_hash.h"
+#include "gromacs/legacyheaders/gmx_omp_nthreads.h"
+#include "gromacs/legacyheaders/macros.h"
+#include "gromacs/legacyheaders/types/commrec.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/pbcutil/ishift.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/gmxassert.h"
+#include "gromacs/utility/smalloc.h"
+
+void dd_move_f_specat(gmx_domdec_t *dd, gmx_domdec_specat_comm_t *spac,
+ rvec *f, rvec *fshift)
+{
+ gmx_specatsend_t *spas;
+ rvec *vbuf;
+ int n, n0, n1, d, dim, dir, i;
+ ivec vis;
+ int is;
+ gmx_bool bPBC, bScrew;
+
+ n = spac->at_end;
+ for (d = dd->ndim-1; d >= 0; d--)
+ {
+ dim = dd->dim[d];
+ if (dd->nc[dim] > 2)
+ {
+ /* Pulse the grid forward and backward */
+ spas = spac->spas[d];
+ n0 = spas[0].nrecv;
+ n1 = spas[1].nrecv;
+ n -= n1 + n0;
+ vbuf = spac->vbuf;
+ /* Send and receive the coordinates */
+ dd_sendrecv2_rvec(dd, d,
+ f+n+n1, n0, vbuf, spas[0].nsend,
+ f+n, n1, vbuf+spas[0].nsend, spas[1].nsend);
+ for (dir = 0; dir < 2; dir++)
+ {
+ bPBC = ((dir == 0 && dd->ci[dim] == 0) ||
+ (dir == 1 && dd->ci[dim] == dd->nc[dim]-1));
+ bScrew = (bPBC && dd->bScrewPBC && dim == XX);
+
+ spas = &spac->spas[d][dir];
+ /* Sum the buffer into the required forces */
+ if (!bPBC || (!bScrew && fshift == NULL))
+ {
+ for (i = 0; i < spas->nsend; i++)
+ {
+ rvec_inc(f[spas->a[i]], *vbuf);
+ vbuf++;
+ }
+ }
+ else
+ {
+ clear_ivec(vis);
+ vis[dim] = (dir == 0 ? 1 : -1);
+ is = IVEC2IS(vis);
+ if (!bScrew)
+ {
+ /* Sum and add to shift forces */
+ for (i = 0; i < spas->nsend; i++)
+ {
+ rvec_inc(f[spas->a[i]], *vbuf);
+ rvec_inc(fshift[is], *vbuf);
+ vbuf++;
+ }
+ }
+ else
+ {
+ /* Rotate the forces */
+ for (i = 0; i < spas->nsend; i++)
+ {
+ f[spas->a[i]][XX] += (*vbuf)[XX];
+ f[spas->a[i]][YY] -= (*vbuf)[YY];
+ f[spas->a[i]][ZZ] -= (*vbuf)[ZZ];
+ if (fshift)
+ {
+ rvec_inc(fshift[is], *vbuf);
+ }
+ vbuf++;
+ }
+ }
+ }
+ }
+ }
+ else
+ {
+ /* Two cells, so we only need to communicate one way */
+ spas = &spac->spas[d][0];
+ n -= spas->nrecv;
+ /* Send and receive the coordinates */
+ dd_sendrecv_rvec(dd, d, dddirForward,
+ f+n, spas->nrecv, spac->vbuf, spas->nsend);
+ /* Sum the buffer into the required forces */
+ if (dd->bScrewPBC && dim == XX &&
+ (dd->ci[dim] == 0 ||
+ dd->ci[dim] == dd->nc[dim]-1))
+ {
+ for (i = 0; i < spas->nsend; i++)
+ {
+ /* Rotate the force */
+ f[spas->a[i]][XX] += spac->vbuf[i][XX];
+ f[spas->a[i]][YY] -= spac->vbuf[i][YY];
+ f[spas->a[i]][ZZ] -= spac->vbuf[i][ZZ];
+ }
+ }
+ else
+ {
+ for (i = 0; i < spas->nsend; i++)
+ {
+ rvec_inc(f[spas->a[i]], spac->vbuf[i]);
+ }
+ }
+ }
+ }
+}
+
+void dd_move_x_specat(gmx_domdec_t *dd, gmx_domdec_specat_comm_t *spac,
+ matrix box,
+ rvec *x0,
+ rvec *x1, gmx_bool bX1IsCoord)
+{
+ gmx_specatsend_t *spas;
+ rvec *x, *vbuf, *rbuf;
+ int nvec, v, n, nn, ns0, ns1, nr0, nr1, nr, d, dim, dir, i;
+ gmx_bool bPBC, bScrew = FALSE;
+ rvec shift = {0, 0, 0};
+
+ nvec = 1;
+ if (x1 != NULL)
+ {
+ nvec++;
+ }
+
+ n = spac->at_start;
+ for (d = 0; d < dd->ndim; d++)
+ {
+ dim = dd->dim[d];
+ if (dd->nc[dim] > 2)
+ {
+ /* Pulse the grid forward and backward */
+ vbuf = spac->vbuf;
+ for (dir = 0; dir < 2; dir++)
+ {
+ if (dir == 0 && dd->ci[dim] == 0)
+ {
+ bPBC = TRUE;
+ bScrew = (dd->bScrewPBC && dim == XX);
+ copy_rvec(box[dim], shift);
+ }
+ else if (dir == 1 && dd->ci[dim] == dd->nc[dim]-1)
+ {
+ bPBC = TRUE;
+ bScrew = (dd->bScrewPBC && dim == XX);
+ for (i = 0; i < DIM; i++)
+ {
+ shift[i] = -box[dim][i];
+ }
+ }
+ else
+ {
+ bPBC = FALSE;
+ bScrew = FALSE;
+ }
+ spas = &spac->spas[d][dir];
+ for (v = 0; v < nvec; v++)
+ {
+ x = (v == 0 ? x0 : x1);
+ /* Copy the required coordinates to the send buffer */
+ if (!bPBC || (v == 1 && !bX1IsCoord))
+ {
+ /* Only copy */
+ for (i = 0; i < spas->nsend; i++)
+ {
+ copy_rvec(x[spas->a[i]], *vbuf);
+ vbuf++;
+ }
+ }
+ else if (!bScrew)
+ {
+ /* Shift coordinates */
+ for (i = 0; i < spas->nsend; i++)
+ {
+ rvec_add(x[spas->a[i]], shift, *vbuf);
+ vbuf++;
+ }
+ }
+ else
+ {
+ /* Shift and rotate coordinates */
+ for (i = 0; i < spas->nsend; i++)
+ {
+ (*vbuf)[XX] = x[spas->a[i]][XX] + shift[XX];
+ (*vbuf)[YY] = box[YY][YY] - x[spas->a[i]][YY] + shift[YY];
+ (*vbuf)[ZZ] = box[ZZ][ZZ] - x[spas->a[i]][ZZ] + shift[ZZ];
+ vbuf++;
+ }
+ }
+ }
+ }
+ /* Send and receive the coordinates */
+ spas = spac->spas[d];
+ ns0 = spas[0].nsend;
+ nr0 = spas[0].nrecv;
+ ns1 = spas[1].nsend;
+ nr1 = spas[1].nrecv;
+ if (nvec == 1)
+ {
+ dd_sendrecv2_rvec(dd, d,
+ spac->vbuf+ns0, ns1, x0+n, nr1,
+ spac->vbuf, ns0, x0+n+nr1, nr0);
+ }
+ else
+ {
+ /* Communicate both vectors in one buffer */
+ rbuf = spac->vbuf2;
+ dd_sendrecv2_rvec(dd, d,
+ spac->vbuf+2*ns0, 2*ns1, rbuf, 2*nr1,
+ spac->vbuf, 2*ns0, rbuf+2*nr1, 2*nr0);
+ /* Split the buffer into the two vectors */
+ nn = n;
+ for (dir = 1; dir >= 0; dir--)
+ {
+ nr = spas[dir].nrecv;
+ for (v = 0; v < 2; v++)
+ {
+ x = (v == 0 ? x0 : x1);
+ for (i = 0; i < nr; i++)
+ {
+ copy_rvec(*rbuf, x[nn+i]);
+ rbuf++;
+ }
+ }
+ nn += nr;
+ }
+ }
+ n += nr0 + nr1;
+ }
+ else
+ {
+ spas = &spac->spas[d][0];
+ /* Copy the required coordinates to the send buffer */
+ vbuf = spac->vbuf;
+ for (v = 0; v < nvec; v++)
+ {
+ x = (v == 0 ? x0 : x1);
+ if (dd->bScrewPBC && dim == XX &&
+ (dd->ci[XX] == 0 || dd->ci[XX] == dd->nc[XX]-1))
+ {
+ /* Here we only perform the rotation, the rest of the pbc
+ * is handled in the constraint or viste routines.
+ */
+ for (i = 0; i < spas->nsend; i++)
+ {
+ (*vbuf)[XX] = x[spas->a[i]][XX];
+ (*vbuf)[YY] = box[YY][YY] - x[spas->a[i]][YY];
+ (*vbuf)[ZZ] = box[ZZ][ZZ] - x[spas->a[i]][ZZ];
+ vbuf++;
+ }
+ }
+ else
+ {
+ for (i = 0; i < spas->nsend; i++)
+ {
+ copy_rvec(x[spas->a[i]], *vbuf);
+ vbuf++;
+ }
+ }
+ }
+ /* Send and receive the coordinates */
+ if (nvec == 1)
+ {
+ dd_sendrecv_rvec(dd, d, dddirBackward,
+ spac->vbuf, spas->nsend, x0+n, spas->nrecv);
+ }
+ else
+ {
+ /* Communicate both vectors in one buffer */
+ rbuf = spac->vbuf2;
+ dd_sendrecv_rvec(dd, d, dddirBackward,
+ spac->vbuf, 2*spas->nsend, rbuf, 2*spas->nrecv);
+ /* Split the buffer into the two vectors */
+ nr = spas[0].nrecv;
+ for (v = 0; v < 2; v++)
+ {
+ x = (v == 0 ? x0 : x1);
+ for (i = 0; i < nr; i++)
+ {
+ copy_rvec(*rbuf, x[n+i]);
+ rbuf++;
+ }
+ }
+ }
+ n += spas->nrecv;
+ }
+ }
+}
+
+int setup_specat_communication(gmx_domdec_t *dd,
+ ind_req_t *ireq,
+ gmx_domdec_specat_comm_t *spac,
+ gmx_hash_t ga2la_specat,
+ int at_start,
+ int vbuf_fac,
+ const char *specat_type,
+ const char *add_err)
+{
+ int nsend[2], nlast, nsend_zero[2] = {0, 0}, *nsend_ptr;
+ int d, dim, ndir, dir, nr, ns, i, nrecv_local, n0, start, indr, ind, buf[2];
+ int nat_tot_specat, nat_tot_prev, nalloc_old;
+ gmx_bool bPBC;
+ gmx_specatsend_t *spas;
+
+ if (debug)
+ {
+ fprintf(debug, "Begin setup_specat_communication for %s\n", specat_type);
+ }
+
+ /* nsend[0]: the number of atoms requested by this node only,
+ * we communicate this for more efficients checks
+ * nsend[1]: the total number of requested atoms
+ */
+ nsend[0] = ireq->n;
+ nsend[1] = nsend[0];
+ nlast = nsend[1];
+ for (d = dd->ndim-1; d >= 0; d--)
+ {
+ /* Pulse the grid forward and backward */
+ dim = dd->dim[d];
+ bPBC = (dim < dd->npbcdim);
+ if (dd->nc[dim] == 2)
+ {
+ /* Only 2 cells, so we only need to communicate once */
+ ndir = 1;
+ }
+ else
+ {
+ ndir = 2;
+ }
+ for (dir = 0; dir < ndir; dir++)
+ {
+ if (!bPBC &&
+ dd->nc[dim] > 2 &&
+ ((dir == 0 && dd->ci[dim] == dd->nc[dim] - 1) ||
+ (dir == 1 && dd->ci[dim] == 0)))
+ {
+ /* No pbc: the fist/last cell should not request atoms */
+ nsend_ptr = nsend_zero;
+ }
+ else
+ {
+ nsend_ptr = nsend;
+ }
+ /* Communicate the number of indices */
+ dd_sendrecv_int(dd, d, dir == 0 ? dddirForward : dddirBackward,
+ nsend_ptr, 2, spac->nreq[d][dir], 2);
+ nr = spac->nreq[d][dir][1];
+ if (nlast+nr > ireq->nalloc)
+ {
+ ireq->nalloc = over_alloc_dd(nlast+nr);
+ srenew(ireq->ind, ireq->nalloc);
+ }
+ /* Communicate the indices */
+ dd_sendrecv_int(dd, d, dir == 0 ? dddirForward : dddirBackward,
+ ireq->ind, nsend_ptr[1], ireq->ind+nlast, nr);
+ nlast += nr;
+ }
+ nsend[1] = nlast;
+ }
+ if (debug)
+ {
+ fprintf(debug, "Communicated the counts\n");
+ }
+
+ /* Search for the requested atoms and communicate the indices we have */
+ nat_tot_specat = at_start;
+ nrecv_local = 0;
+ for (d = 0; d < dd->ndim; d++)
+ {
+ /* Pulse the grid forward and backward */
+ if (dd->dim[d] >= dd->npbcdim || dd->nc[dd->dim[d]] > 2)
+ {
+ ndir = 2;
+ }
+ else
+ {
+ ndir = 1;
+ }
+ nat_tot_prev = nat_tot_specat;
+ for (dir = ndir-1; dir >= 0; dir--)
+ {
+ if (nat_tot_specat > spac->bSendAtom_nalloc)
+ {
+ nalloc_old = spac->bSendAtom_nalloc;
+ spac->bSendAtom_nalloc = over_alloc_dd(nat_tot_specat);
+ srenew(spac->bSendAtom, spac->bSendAtom_nalloc);
+ for (i = nalloc_old; i < spac->bSendAtom_nalloc; i++)
+ {
+ spac->bSendAtom[i] = FALSE;
+ }
+ }
+ spas = &spac->spas[d][dir];
+ n0 = spac->nreq[d][dir][0];
+ nr = spac->nreq[d][dir][1];
+ if (debug)
+ {
+ fprintf(debug, "dim=%d, dir=%d, searching for %d atoms\n",
+ d, dir, nr);
+ }
+ start = nlast - nr;
+ spas->nsend = 0;
+ nsend[0] = 0;
+ for (i = 0; i < nr; i++)
+ {
+ indr = ireq->ind[start+i];
+ ind = -1;
+ /* Check if this is a home atom and if so ind will be set */
+ if (!ga2la_get_home(dd->ga2la, indr, &ind))
+ {
+ /* Search in the communicated atoms */
+ ind = gmx_hash_get_minone(ga2la_specat, indr);
+ }
+ if (ind >= 0)
+ {
+ if (i < n0 || !spac->bSendAtom[ind])
+ {
+ if (spas->nsend+1 > spas->a_nalloc)
+ {
+ spas->a_nalloc = over_alloc_large(spas->nsend+1);
+ srenew(spas->a, spas->a_nalloc);
+ }
+ /* Store the local index so we know which coordinates
+ * to send out later.
+ */
+ spas->a[spas->nsend] = ind;
+ spac->bSendAtom[ind] = TRUE;
+ if (spas->nsend+1 > spac->ibuf_nalloc)
+ {
+ spac->ibuf_nalloc = over_alloc_large(spas->nsend+1);
+ srenew(spac->ibuf, spac->ibuf_nalloc);
+ }
+ /* Store the global index so we can send it now */
+ spac->ibuf[spas->nsend] = indr;
+ if (i < n0)
+ {
+ nsend[0]++;
+ }
+ spas->nsend++;
+ }
+ }
+ }
+ nlast = start;
+ /* Clear the local flags */
+ for (i = 0; i < spas->nsend; i++)
+ {
+ spac->bSendAtom[spas->a[i]] = FALSE;
+ }
+ /* Send and receive the number of indices to communicate */
+ nsend[1] = spas->nsend;
+ dd_sendrecv_int(dd, d, dir == 0 ? dddirBackward : dddirForward,
+ nsend, 2, buf, 2);
+ if (debug)
+ {
+ fprintf(debug, "Send to rank %d, %d (%d) indices, "
+ "receive from rank %d, %d (%d) indices\n",
+ dd->neighbor[d][1-dir], nsend[1], nsend[0],
+ dd->neighbor[d][dir], buf[1], buf[0]);
+ if (gmx_debug_at)
+ {
+ for (i = 0; i < spas->nsend; i++)
+ {
+ fprintf(debug, " %d", spac->ibuf[i]+1);
+ }
+ fprintf(debug, "\n");
+ }
+ }
+ nrecv_local += buf[0];
+ spas->nrecv = buf[1];
+ if (nat_tot_specat + spas->nrecv > dd->gatindex_nalloc)
+ {
+ dd->gatindex_nalloc =
+ over_alloc_dd(nat_tot_specat + spas->nrecv);
+ srenew(dd->gatindex, dd->gatindex_nalloc);
+ }
+ /* Send and receive the indices */
+ dd_sendrecv_int(dd, d, dir == 0 ? dddirBackward : dddirForward,
+ spac->ibuf, spas->nsend,
+ dd->gatindex+nat_tot_specat, spas->nrecv);
+ nat_tot_specat += spas->nrecv;
+ }
+
+ /* Allocate the x/f communication buffers */
+ ns = spac->spas[d][0].nsend;
+ nr = spac->spas[d][0].nrecv;
+ if (ndir == 2)
+ {
+ ns += spac->spas[d][1].nsend;
+ nr += spac->spas[d][1].nrecv;
+ }
+ if (vbuf_fac*ns > spac->vbuf_nalloc)
+ {
+ spac->vbuf_nalloc = over_alloc_dd(vbuf_fac*ns);
+ srenew(spac->vbuf, spac->vbuf_nalloc);
+ }
+ if (vbuf_fac == 2 && vbuf_fac*nr > spac->vbuf2_nalloc)
+ {
+ spac->vbuf2_nalloc = over_alloc_dd(vbuf_fac*nr);
+ srenew(spac->vbuf2, spac->vbuf2_nalloc);
+ }
+
+ /* Make a global to local index for the communication atoms */
+ for (i = nat_tot_prev; i < nat_tot_specat; i++)
+ {
+ gmx_hash_change_or_set(ga2la_specat, dd->gatindex[i], i);
+ }
+ }
+
+ /* Check that in the end we got the number of atoms we asked for */
+ if (nrecv_local != ireq->n)
+ {
+ if (debug)
+ {
+ fprintf(debug, "Requested %d, received %d (tot recv %d)\n",
+ ireq->n, nrecv_local, nat_tot_specat-at_start);
+ if (gmx_debug_at)
+ {
+ for (i = 0; i < ireq->n; i++)
+ {
+ ind = gmx_hash_get_minone(ga2la_specat, ireq->ind[i]);
+ fprintf(debug, " %s%d",
+ (ind >= 0) ? "" : "!",
+ ireq->ind[i]+1);
+ }
+ fprintf(debug, "\n");
+ }
+ }
+ fprintf(stderr, "\nDD cell %d %d %d: Neighboring cells do not have atoms:",
+ dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
+ for (i = 0; i < ireq->n; i++)
+ {
+ if (gmx_hash_get_minone(ga2la_specat, ireq->ind[i]) < 0)
+ {
+ fprintf(stderr, " %d", ireq->ind[i]+1);
+ }
+ }
+ fprintf(stderr, "\n");
+ gmx_fatal(FARGS, "DD cell %d %d %d could only obtain %d of the %d atoms that are connected via %ss from the neighboring cells. This probably means your %s lengths are too long compared to the domain decomposition cell size. Decrease the number of domain decomposition grid cells%s%s.",
+ dd->ci[XX], dd->ci[YY], dd->ci[ZZ],
+ nrecv_local, ireq->n, specat_type,
+ specat_type, add_err,
+ dd->bGridJump ? " or use the -rcon option of mdrun" : "");
+ }
+
+ spac->at_start = at_start;
+ spac->at_end = nat_tot_specat;
+
+ if (debug)
+ {
+ fprintf(debug, "Done setup_specat_communication\n");
+ }
+
+ return nat_tot_specat;
+}
+
+gmx_domdec_specat_comm_t *specat_comm_init(int nthread)
+{
+ gmx_domdec_specat_comm_t *spac;
+
+ snew(spac, 1);
+ spac->nthread = nthread;
+ snew(spac->ireq, spac->nthread);
+
+ return spac;
+}
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2005,2006,2007,2008,2009,2010,2012,2013,2014, 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.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ *
+ * \brief This file declares functions for domdec to use
+ * while managing communication of atoms required for special purposes
+ *
+ * \author Berk Hess <hess@kth.se>
+ * \ingroup module_domdec
+ */
+
+#ifndef GMX_DOMDEC_DOMDEC_SPECATOMCOMM_H
+#define GMX_DOMDEC_DOMDEC_SPECATOMCOMM_H
+
+#include "gromacs/legacyheaders/types/commrec.h"
+#include "gromacs/math/vectypes.h"
+#include "gromacs/utility/basedefinitions.h"
+
+typedef struct {
+ int nsend;
+ int *a;
+ int a_nalloc;
+ int nrecv;
+} gmx_specatsend_t;
+
+typedef struct {
+ int *ind;
+ int nalloc;
+ int n;
+} ind_req_t;
+
+/*! \internal \brief Struct with setup and buffers for special atom communication */
+typedef struct gmx_domdec_specat_comm {
+ /* The number of indices to receive during the setup */
+ int nreq[DIM][2][2]; /**< The nr. of atoms requested, per DIM, direction and direct/indirect */
+ /* The atoms to send */
+ gmx_specatsend_t spas[DIM][2]; /**< The communication setup per DIM, direction */
+ gmx_bool *bSendAtom; /**< Work buffer that tells if spec.atoms should be sent */
+ int bSendAtom_nalloc; /**< Allocation size of \p bSendAtom */
+ /* Send buffers */
+ int *ibuf; /**< Integer send buffer */
+ int ibuf_nalloc; /**< Allocation size of \p ibuf */
+ rvec *vbuf; /**< rvec send buffer */
+ int vbuf_nalloc; /**< Allocation size of \p vbuf */
+ rvec *vbuf2; /**< rvec send buffer */
+ int vbuf2_nalloc; /**< Allocation size of \p vbuf2 */
+ /* The range in the local buffer(s) for received atoms */
+ int at_start; /**< Start index of received atoms */
+ int at_end; /**< End index of received atoms */
+
+ /* The atom indices we need from the surrounding cells.
+ * We can gather the indices over nthread threads.
+ */
+ int nthread; /**< Number of threads used for spec.atom communication */
+ ind_req_t *ireq; /**< Index request buffer per thread, allocation size \p nthread */
+} gmx_domdec_specat_comm_t;
+
+/*! \brief Communicates the force for special atoms, the shift forces are reduced with \p fshift != NULL */
+void dd_move_f_specat(gmx_domdec_t *dd, gmx_domdec_specat_comm_t *spac,
+ rvec *f, rvec *fshift);
+
+/*! \brief Communicates the coordinates for special atoms
+ *
+ * \param[in] dd Domain decomposition struct
+ * \param[in] spac Special atom communication struct
+ * \param[in] box Box, used for pbc
+ * \param[in,out] x0 Vector to communicate
+ * \param[in,out] x1 Vector to communicate, when != NULL
+ * \param[in] bX1IsCoord Tells is \p x1 is a coordinate vector (needs pbc)
+ */
+void dd_move_x_specat(gmx_domdec_t *dd, gmx_domdec_specat_comm_t *spac,
+ matrix box,
+ rvec *x0,
+ rvec *x1, gmx_bool bX1IsCoord);
+
+/*! \brief Sets up the communication for special atoms
+ *
+ * \param[in] dd Domain decomposition struct
+ * \param[in] ireq List of requested atom indices
+ * \param[in,out] spac Special atom communication struct
+ * \param[out] ga2la_specat Global to local special atom index
+ * \param[in] at_start Index in local state where to start storing communicated atoms
+ * \param[in] vbuf_fac Buffer factor, 1 or 2 for communicating 1 or 2 vectors
+ * \param[in] specat_type Name of the special atom, used for error message
+ * \param[in] add_err Text to add at the end of error message when atoms can't be found
+ */
+int setup_specat_communication(gmx_domdec_t *dd,
+ ind_req_t *ireq,
+ gmx_domdec_specat_comm_t *spac,
+ gmx_hash_t ga2la_specat,
+ int at_start,
+ int vbuf_fac,
+ const char *specat_type,
+ const char *add_err);
+
+/*! \brief Initialize a special communication struct */
+gmx_domdec_specat_comm_t *specat_comm_init(int nthread);
+
+#endif
* the research papers on the package. Check out http://www.gromacs.org.
*/
+/*! \internal \file
+ *
+ * \brief This file defines functions used by the domdec module
+ * while managing the construction, use and error checking for
+ * topologies local to a DD rank.
+ *
+ * \author Berk Hess <hess@kth.se>
+ * \ingroup module_domdec
+ */
+
#include "gmxpre.h"
#include <string.h>
#include "gromacs/utility/gmxassert.h"
#include "gromacs/utility/smalloc.h"
-/* for dd_init_local_state */
+#include "domdec_constraints.h"
+#include "domdec_internal.h"
+#include "domdec_vsite.h"
+
+/*! \brief The number of integer item in the local state, used for broadcasting of the state */
#define NITEM_DD_INIT_LOCAL_STATE 5
typedef struct {
int *index; /* Index for each atom into il */
int *il; /* ftype|type|a0|...|an|ftype|... */
-} gmx_reverse_ilist_t;
+} reverse_ilist_t;
typedef struct {
int a_start;
int a_end;
int natoms_mol;
int type;
-} gmx_molblock_ind_t;
+} molblock_ind_t;
+/*! \brief Struct for thread local work data for local topology generation */
+typedef struct {
+ t_idef idef; /**< Partial local topology */
+ int **vsite_pbc; /**< vsite PBC structure */
+ int *vsite_pbc_nalloc; /**< Allocation sizes for vsite_pbc */
+ int nbonded; /**< The number of bondeds in this struct */
+ t_blocka excl; /**< List of exclusions */
+ int excl_count; /**< The total exclusion count for \p excl */
+} thread_work_t;
+
+/*! \brief Struct for the reverse topology: links bonded interactions to atomsx */
typedef struct gmx_reverse_top {
- gmx_bool bExclRequired; /* Do we require all exclusions to be assigned? */
- gmx_bool bConstr; /* Are there constraints in this revserse top? */
- gmx_bool bSettle; /* Are there settles in this revserse top? */
- gmx_bool bBCheck; /* All bonded interactions have to be assigned? */
- gmx_bool bMultiCGmols; /* Are the multi charge-group molecules? */
- gmx_reverse_ilist_t *ril_mt; /* Reverse ilist for all moltypes */
- int ril_mt_tot_size;
- int ilsort; /* The sorting state of bondeds for free energy */
- gmx_molblock_ind_t *mbi;
- int nmolblock;
+ //! @cond Doxygen_Suppress
+ gmx_bool bExclRequired; /**< Do we require all exclusions to be assigned? */
+ gmx_bool bConstr; /**< Are there constraints in this revserse top? */
+ gmx_bool bSettle; /**< Are there settles in this revserse top? */
+ gmx_bool bBCheck; /**< All bonded interactions have to be assigned? */
+ gmx_bool bMultiCGmols; /**< Are the multi charge-group molecules? */
+ reverse_ilist_t *ril_mt; /**< Reverse ilist for all moltypes */
+ int ril_mt_tot_size; /**< The size of ril_mt[?].index summed over all entries */
+ int ilsort; /**< The sorting state of bondeds for free energy */
+ molblock_ind_t *mbi; /**< molblock to global atom index for quick lookup of molblocks on atom index */
+ int nmolblock; /**< The number of molblocks, size of \p mbi */
/* Work data structures for multi-threading */
- int nthread;
- t_idef *idef_thread;
- int ***vsite_pbc;
- int **vsite_pbc_nalloc;
- int *nbonded_thread;
- t_blocka *excl_thread;
- int *excl_count_thread;
+ int nthread; /**< The number of threads to be used */
+ thread_work_t *th_work; /**< Thread work array for local topology generation */
/* Pointers only used for an error message */
- gmx_mtop_t *err_top_global;
- gmx_localtop_t *err_top_local;
+ gmx_mtop_t *err_top_global; /**< Pointer to the global top, only used for error reporting */
+ gmx_localtop_t *err_top_local; /**< Pointer to the local top, only used for error reporting */
+ //! @endcond
} gmx_reverse_top_t;
+/*! \brief Returns the number of atom entries for il in gmx_reverse_top_t */
static int nral_rt(int ftype)
{
- /* Returns the number of atom entries for il in gmx_reverse_top_t */
int nral;
nral = NRAL(ftype);
return nral;
}
-/* This function tells which interactions need to be assigned exactly once */
+/*! \brief Return whether interactions of type \p ftype need to be assigned exactly once */
static gmx_bool dd_check_ftype(int ftype, gmx_bool bBCheck,
gmx_bool bConstr, gmx_bool bSettle)
{
(bSettle && ftype == F_SETTLE));
}
+/*! \brief Print a header on error messages */
static void print_error_header(FILE *fplog, char *moltypename, int nprint)
{
fprintf(fplog, "\nMolecule type '%s'\n", moltypename);
nprint);
}
+/*! \brief Help print error output when interactions are missing */
static void print_missing_interactions_mb(FILE *fplog, t_commrec *cr,
gmx_reverse_top_t *rt,
char *moltypename,
- gmx_reverse_ilist_t *ril,
+ reverse_ilist_t *ril,
int a_start, int a_end,
int nat_mol, int nmol,
t_idef *idef)
sfree(assigned);
}
+/*! \brief Help print error output when interactions are missing */
static void print_missing_interactions_atoms(FILE *fplog, t_commrec *cr,
gmx_mtop_t *mtop, t_idef *idef)
{
}
else
{
- gmx_fatal(FARGS, "%d of the %d bonded interactions could not be calculated because some atoms involved moved further apart than the multi-body cut-off distance (%g nm) or the two-body cut-off distance (%g nm), see option -rdd, for pairs and tabulated bonds also see option -ddcheck", -ndiff_tot, cr->dd->nbonded_global, dd_cutoff_mbody(cr->dd), dd_cutoff_twobody(cr->dd));
+ gmx_fatal(FARGS, "%d of the %d bonded interactions could not be calculated because some atoms involved moved further apart than the multi-body cut-off distance (%g nm) or the two-body cut-off distance (%g nm), see option -rdd, for pairs and tabulated bonds also see option -ddcheck", -ndiff_tot, cr->dd->nbonded_global, dd_cutoff_multibody(cr->dd), dd_cutoff_twobody(cr->dd));
}
}
}
+/*! \brief Return global topology molecule information for global atom index \p i_gl */
static void global_atomnr_to_moltype_ind(gmx_reverse_top_t *rt, int i_gl,
int *mb, int *mt, int *mol, int *i_mol)
{
- gmx_molblock_ind_t *mbi = rt->mbi;
- int start = 0;
- int end = rt->nmolblock; /* exclusive */
- int mid;
+ molblock_ind_t *mbi = rt->mbi;
+ int start = 0;
+ int end = rt->nmolblock; /* exclusive */
+ int mid;
/* binary search for molblock_ind */
while (TRUE)
*i_mol = (i_gl - mbi->a_start) - (*mol)*mbi->natoms_mol;
}
-static int count_excls(t_block *cgs, t_blocka *excls, int *n_intercg_excl)
+/*! \brief Count the exclusions for all atoms in \p cgs */
+static int count_excls(const t_block *cgs, const t_blocka *excls,
+ int *n_intercg_excl)
{
int n, cg, at0, at1, at, excl, atj;
return n;
}
+/*! \brief Run the reverse ilist generation and store it when \p bAssign = TRUE */
static int low_make_reverse_ilist(t_ilist *il_mt, t_atom *atom,
int **vsite_pbc,
int *count,
return nint;
}
+/*! \brief Make the reverse ilist: a list of bonded interactions linked to atoms */
static int make_reverse_ilist(gmx_moltype_t *molt,
int **vsite_pbc,
gmx_bool bConstr, gmx_bool bSettle,
gmx_bool bBCheck,
gmx_bool bLinkToAllAtoms,
- gmx_reverse_ilist_t *ril_mt)
+ reverse_ilist_t *ril_mt)
{
int nat_mt, *count, i, nint_mt;
return nint_mt;
}
-static void destroy_reverse_ilist(gmx_reverse_ilist_t *ril)
+/*! \brief Destroys a reverse_ilist_t struct */
+static void destroy_reverse_ilist(reverse_ilist_t *ril)
{
sfree(ril->index);
sfree(ril->il);
}
+/*! \brief Generate the reverse topology */
static gmx_reverse_top_t *make_reverse_top(gmx_mtop_t *mtop, gmx_bool bFE,
int ***vsite_pbc_molt,
gmx_bool bConstr, gmx_bool bSettle,
}
rt->nthread = gmx_omp_nthreads_get(emntDomdec);
- snew(rt->idef_thread, rt->nthread);
+ snew(rt->th_work, rt->nthread);
if (vsite_pbc_molt != NULL)
{
- snew(rt->vsite_pbc, rt->nthread);
- snew(rt->vsite_pbc_nalloc, rt->nthread);
for (thread = 0; thread < rt->nthread; thread++)
{
- snew(rt->vsite_pbc[thread], F_VSITEN-F_VSITE2+1);
- snew(rt->vsite_pbc_nalloc[thread], F_VSITEN-F_VSITE2+1);
+ snew(rt->th_work[thread].vsite_pbc, F_VSITEN-F_VSITE2+1);
+ snew(rt->th_work[thread].vsite_pbc_nalloc, F_VSITEN-F_VSITE2+1);
}
}
- snew(rt->nbonded_thread, rt->nthread);
- snew(rt->excl_thread, rt->nthread);
- snew(rt->excl_count_thread, rt->nthread);
return rt;
}
}
}
-/* This routine is very similar to add_ifunc, but vsites interactions
+/*! \brief Store a vsite interaction at the end of \p il
+ *
+ * This routine is very similar to add_ifunc, but vsites interactions
* have more work to do than other kinds of interactions, and the
* complex way nral (and thus vector contents) depends on ftype
* confuses static analysis tools unless we fuse the vsite
}
}
+/*! \brief Store a bonded interaction at the end of \p il */
static gmx_inline void add_ifunc(int nral, t_iatom *tiatoms, t_ilist *il)
{
t_iatom *liatoms;
il->nr += 1 + nral;
}
+/*! \brief Store a position restraint in idef and iatoms, complex because the parameters are different for each entry */
static void add_posres(int mol, int a_mol, const gmx_molblock_t *molb,
t_iatom *iatoms, const t_iparams *ip_in,
t_idef *idef)
iatoms[0] = n;
}
+/*! \brief Store a flat-bottomed position restraint in idef and iatoms, complex because the parameters are different for each entry */
static void add_fbposres(int mol, int a_mol, const gmx_molblock_t *molb,
t_iatom *iatoms, const t_iparams *ip_in,
t_idef *idef)
iatoms[0] = n;
}
+/*! \brief Store a virtual site interaction, complex because of PBC and recursion */
static void add_vsite(gmx_ga2la_t ga2la, int *index, int *rtil,
int ftype, int nral,
gmx_bool bHomeA, int a, int a_gl, int a_mol,
}
}
+/*! \brief Update the local atom to local charge group index */
static void make_la2lc(gmx_domdec_t *dd)
{
int *cgindex, *la2lc, cg, a;
}
}
+/*! \brief Returns the squared distance between charge groups \p i and \p j */
static real dd_dist2(t_pbc *pbc_null, rvec *cg_cm, const int *la2lc, int i, int j)
{
rvec dx;
return norm2(dx);
}
-/* Append the nsrc t_blocka block structures in src to *dest */
-static void combine_blocka(t_blocka *dest, const t_blocka *src, int nsrc)
+/*! \brief Append t_blocka block structures 1 to nsrc in src to *dest */
+static void combine_blocka(t_blocka *dest, const thread_work_t *src, int nsrc)
{
int ni, na, s, i;
- ni = src[nsrc-1].nr;
+ ni = src[nsrc-1].excl.nr;
na = 0;
for (s = 0; s < nsrc; s++)
{
- na += src[s].nra;
+ na += src[s].excl.nra;
}
if (ni + 1 > dest->nalloc_index)
{
dest->nalloc_a = over_alloc_large(dest->nra+na);
srenew(dest->a, dest->nalloc_a);
}
- for (s = 0; s < nsrc; s++)
+ for (s = 1; s < nsrc; s++)
{
- for (i = dest->nr+1; i < src[s].nr+1; i++)
+ for (i = dest->nr+1; i < src[s].excl.nr+1; i++)
{
- dest->index[i] = dest->nra + src[s].index[i];
+ dest->index[i] = dest->nra + src[s].excl.index[i];
}
- for (i = 0; i < src[s].nra; i++)
+ for (i = 0; i < src[s].excl.nra; i++)
{
- dest->a[dest->nra+i] = src[s].a[i];
+ dest->a[dest->nra+i] = src[s].excl.a[i];
}
- dest->nr = src[s].nr;
- dest->nra += src[s].nra;
+ dest->nr = src[s].excl.nr;
+ dest->nra += src[s].excl.nra;
}
}
-/* Append the nsrc t_idef structures in src to *dest,
+/*! \brief Append t_idef structures 1 to nsrc in src to *dest,
* virtual sites need special attention, as pbc info differs per vsite.
*/
-static void combine_idef(t_idef *dest, const t_idef *src, int nsrc,
- gmx_vsite_t *vsite, int ***vsite_pbc_t)
+static void combine_idef(t_idef *dest, const thread_work_t *src, int nsrc,
+ gmx_vsite_t *vsite)
{
- int ftype, n, s, i;
- t_ilist *ild;
- const t_ilist *ils;
- gmx_bool vpbc;
- int nral1 = 0, ftv = 0;
+ int ftype;
for (ftype = 0; ftype < F_NRE; ftype++)
{
+ int n, s;
+
n = 0;
- for (s = 0; s < nsrc; s++)
+ for (s = 1; s < nsrc; s++)
{
- n += src[s].il[ftype].nr;
+ n += src[s].idef.il[ftype].nr;
}
if (n > 0)
{
+ t_ilist *ild;
+
ild = &dest->il[ftype];
if (ild->nr + n > ild->nalloc)
srenew(ild->iatoms, ild->nalloc);
}
+ gmx_bool vpbc;
+ int nral1 = 0, ftv = 0;
+
vpbc = ((interaction_function[ftype].flags & IF_VSITE) &&
vsite->vsite_pbc_loc != NULL);
if (vpbc)
}
}
- for (s = 0; s < nsrc; s++)
+ for (s = 1; s < nsrc; s++)
{
- ils = &src[s].il[ftype];
+ const t_ilist *ils;
+ int i;
+
+ ils = &src[s].idef.il[ftype];
for (i = 0; i < ils->nr; i++)
{
ild->iatoms[ild->nr+i] = ils->iatoms[i];
for (i = 0; i < ils->nr; i += nral1)
{
vsite->vsite_pbc_loc[ftv][(ild->nr+i)/nral1] =
- vsite_pbc_t[s][ftv][i/nral1];
+ src[s].vsite_pbc[ftv][i/nral1];
}
}
/* Position restraints need an additional treatment */
if (dest->il[F_POSRES].nr > 0)
{
+ int n, s, i;
+
n = dest->il[F_POSRES].nr/2;
if (n > dest->iparams_posres_nalloc)
{
srenew(dest->iparams_posres, dest->iparams_posres_nalloc);
}
/* Set n to the number of original position restraints in dest */
- for (s = 0; s < nsrc; s++)
+ for (s = 1; s < nsrc; s++)
{
- n -= src[s].il[F_POSRES].nr/2;
+ n -= src[s].idef.il[F_POSRES].nr/2;
}
- for (s = 0; s < nsrc; s++)
+ for (s = 1; s < nsrc; s++)
{
- for (i = 0; i < src[s].il[F_POSRES].nr/2; i++)
+ for (i = 0; i < src[s].idef.il[F_POSRES].nr/2; i++)
{
/* Correct the index into iparams_posres */
dest->il[F_POSRES].iatoms[n*2] = n;
/* Copy the position restraint force parameters */
- dest->iparams_posres[n] = src[s].iparams_posres[i];
+ dest->iparams_posres[n] = src[s].idef.iparams_posres[i];
n++;
}
}
}
}
-/* This function looks up and assigns bonded interactions for zone iz.
+/*! \brief This function looks up and assigns bonded interactions for zone iz.
+ *
* With thread parallelizing each thread acts on a different atom range:
* at_start to at_end.
*/
return nbonded_local;
}
+/*! \brief Set the exclusion data for i-zone \p iz for the case of no exclusions */
static void set_no_exclusions_zone(gmx_domdec_t *dd, gmx_domdec_zones_t *zones,
int iz, t_blocka *lexcls)
{
}
}
+/*! \brief Set the exclusion data for i-zone \p iz */
static int make_exclusions_zone(gmx_domdec_t *dd, gmx_domdec_zones_t *zones,
const gmx_moltype_t *moltype,
gmx_bool bRCheck, real rc2,
return count;
}
+/*! \brief Ensure we have enough space in \p ba for \p nindex_max indices */
static void check_alloc_index(t_blocka *ba, int nindex_max)
{
if (nindex_max+1 > ba->nalloc_index)
}
}
+/*! \brief Ensure that we have enough space for exclusion storate in \p lexcls */
static void check_exclusions_alloc(gmx_domdec_t *dd, gmx_domdec_zones_t *zones,
t_blocka *lexcls)
{
for (thread = 1; thread < dd->reverse_top->nthread; thread++)
{
- check_alloc_index(&dd->reverse_top->excl_thread[thread], nr);
+ check_alloc_index(&dd->reverse_top->th_work[thread].excl, nr);
}
}
+/*! \brief Set the total count indexes for the local exclusions, needed by several functions */
static void finish_local_exclusions(gmx_domdec_t *dd, gmx_domdec_zones_t *zones,
t_blocka *lexcls)
{
}
}
+/*! \brief Clear a t_idef struct */
static void clear_idef(t_idef *idef)
{
int ftype;
}
}
+/*! \brief Generate and store all required local bonded interactions in \p idef and local exclusions in \p lexcls */
static int make_local_bondeds_excls(gmx_domdec_t *dd,
gmx_domdec_zones_t *zones,
const gmx_mtop_t *mtop,
}
else
{
- idef_t = &rt->idef_thread[thread];
+ idef_t = &rt->th_work[thread].idef;
clear_idef(idef_t);
}
}
else
{
- vsite_pbc = rt->vsite_pbc[thread];
- vsite_pbc_nalloc = rt->vsite_pbc_nalloc[thread];
+ vsite_pbc = rt->th_work[thread].vsite_pbc;
+ vsite_pbc_nalloc = rt->th_work[thread].vsite_pbc_nalloc;
}
}
else
vsite_pbc_nalloc = NULL;
}
- rt->nbonded_thread[thread] =
+ rt->th_work[thread].nbonded =
make_bondeds_zone(dd, zones,
mtop->molblock,
bRCheckMB, rcheck, bRCheck2B, rc2,
}
else
{
- excl_t = &rt->excl_thread[thread];
+ excl_t = &rt->th_work[thread].excl;
excl_t->nr = 0;
excl_t->nra = 0;
}
- rt->excl_count_thread[thread] =
+ rt->th_work[thread].excl_count =
make_exclusions_zone(dd, zones,
mtop->moltype, bRCheck2B, rc2,
la2lc, pbc_null, cg_cm, cginfo,
if (rt->nthread > 1)
{
- combine_idef(idef, rt->idef_thread+1, rt->nthread-1,
- vsite, rt->vsite_pbc+1);
+ combine_idef(idef, rt->th_work, rt->nthread, vsite);
}
for (thread = 0; thread < rt->nthread; thread++)
{
- nbonded_local += rt->nbonded_thread[thread];
+ nbonded_local += rt->th_work[thread].nbonded;
}
if (iz < nzone_excl)
{
if (rt->nthread > 1)
{
- combine_blocka(lexcls, rt->excl_thread+1, rt->nthread-1);
+ combine_blocka(lexcls, rt->th_work, rt->nthread);
}
for (thread = 0; thread < rt->nthread; thread++)
{
- *excl_count += rt->excl_count_thread[thread];
+ *excl_count += rt->th_work[thread].excl_count;
}
}
}
state_local->flags = buf[0];
}
+/*! \brief Check if a link is stored in \p link between charge groups \p cg_gl and \p cg_gl_j and if not so, store a link */
static void check_link(t_blocka *link, int cg_gl, int cg_gl_j)
{
int k;
}
}
+/*! \brief Allocate and return an array of the charge group index for all atoms */
static int *make_at2cg(t_block *cgs)
{
int *at2cg, cg, a;
t_block *cgs;
t_blocka *excls;
int *a2c;
- gmx_reverse_ilist_t ril;
+ reverse_ilist_t ril;
t_blocka *link;
cginfo_mb_t *cgi_mb;
return link;
}
+/*! \brief Set the distance, function type and atom indices for the longest distance between molecules of molecule type \p molt for two-body and multi-body bonded interactions */
static void bonded_cg_distance_mol(gmx_moltype_t *molt, int *at2cg,
gmx_bool bBCheck, gmx_bool bExcl, rvec *cg_cm,
real *r_2b, int *ft2b, int *a2_1, int *a2_2,
*r_mb = sqrt(r2_mb);
}
+//! Compute charge group centers of mass for molecule \p molt
static void get_cgcm_mol(gmx_moltype_t *molt, gmx_ffparams_t *ffparams,
int ePBC, t_graph *graph, matrix box,
gmx_vsite_t *vsite,
calc_cgcm(NULL, 0, molt->cgs.nr, &molt->cgs, xs, cg_cm);
}
+//! Returns whether \p molt has a virtual site
static int have_vsite_molt(gmx_moltype_t *molt)
{
int i;
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2006,2007,2008,2009,2010,2012,2013,2014, 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.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+
+/*! \internal \file
+ *
+ * \brief This file implements functions for domdec to use
+ * while managing inter-atomic constraints.
+ *
+ * \author Berk Hess <hess@kth.se>
+ * \ingroup module_domdec
+ */
+
+#include "gmxpre.h"
+
+#include "domdec_vsite.h"
+
+#include <assert.h>
+
+#include <algorithm>
+
+#include "gromacs/domdec/domdec.h"
+#include "gromacs/legacyheaders/gmx_ga2la.h"
+#include "gromacs/legacyheaders/gmx_hash.h"
+#include "gromacs/legacyheaders/gmx_omp_nthreads.h"
+#include "gromacs/legacyheaders/macros.h"
+#include "gromacs/legacyheaders/types/commrec.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/pbcutil/ishift.h"
+#include "gromacs/topology/mtop_util.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/gmxassert.h"
+#include "gromacs/utility/smalloc.h"
+
+#include "domdec_specatomcomm.h"
+
+void dd_move_f_vsites(gmx_domdec_t *dd, rvec *f, rvec *fshift)
+{
+ if (dd->vsite_comm)
+ {
+ dd_move_f_specat(dd, dd->vsite_comm, f, fshift);
+ }
+}
+
+void dd_clear_f_vsites(gmx_domdec_t *dd, rvec *f)
+{
+ int i;
+
+ if (dd->vsite_comm)
+ {
+ for (i = dd->vsite_comm->at_start; i < dd->vsite_comm->at_end; i++)
+ {
+ clear_rvec(f[i]);
+ }
+ }
+}
+
+void dd_move_x_vsites(gmx_domdec_t *dd, matrix box, rvec *x)
+{
+ if (dd->vsite_comm)
+ {
+ dd_move_x_specat(dd, dd->vsite_comm, box, x, NULL, FALSE);
+ }
+}
+
+int dd_make_local_vsites(gmx_domdec_t *dd, int at_start, t_ilist *lil)
+{
+ gmx_domdec_specat_comm_t *spac;
+ ind_req_t *ireq;
+ gmx_hash_t ga2la_specat;
+ int ftype, nral, i, j, a;
+ t_ilist *lilf;
+ t_iatom *iatoms;
+ int at_end;
+
+ spac = dd->vsite_comm;
+ ireq = &spac->ireq[0];
+ ga2la_specat = dd->ga2la_vsite;
+
+ ireq->n = 0;
+ /* Loop over all the home vsites */
+ for (ftype = 0; ftype < F_NRE; ftype++)
+ {
+ if (interaction_function[ftype].flags & IF_VSITE)
+ {
+ nral = NRAL(ftype);
+ lilf = &lil[ftype];
+ for (i = 0; i < lilf->nr; i += 1+nral)
+ {
+ iatoms = lilf->iatoms + i;
+ /* Check if we have the other atoms */
+ for (j = 1; j < 1+nral; j++)
+ {
+ if (iatoms[j] < 0)
+ {
+ /* This is not a home atom,
+ * we need to ask our neighbors.
+ */
+ a = -iatoms[j] - 1;
+ /* Check to not ask for the same atom more than once */
+ if (gmx_hash_get_minone(dd->ga2la_vsite, a) == -1)
+ {
+ /* Add this non-home atom to the list */
+ if (ireq->n+1 > ireq->nalloc)
+ {
+ ireq->nalloc = over_alloc_large(ireq->n+1);
+ srenew(ireq->ind, ireq->nalloc);
+ }
+ ireq->ind[ireq->n++] = a;
+ /* Temporarily mark with -2,
+ * we get the index later.
+ */
+ gmx_hash_set(ga2la_specat, a, -2);
+ }
+ }
+ }
+ }
+ }
+ }
+
+ at_end = setup_specat_communication(dd, ireq, dd->vsite_comm, ga2la_specat,
+ at_start, 1, "vsite", "");
+
+ /* Fill in the missing indices */
+ for (ftype = 0; ftype < F_NRE; ftype++)
+ {
+ if (interaction_function[ftype].flags & IF_VSITE)
+ {
+ nral = NRAL(ftype);
+ lilf = &lil[ftype];
+ for (i = 0; i < lilf->nr; i += 1+nral)
+ {
+ iatoms = lilf->iatoms + i;
+ for (j = 1; j < 1+nral; j++)
+ {
+ if (iatoms[j] < 0)
+ {
+ iatoms[j] = gmx_hash_get_minone(ga2la_specat, -iatoms[j]-1);
+ }
+ }
+ }
+ }
+ }
+
+ return at_end;
+}
+
+void init_domdec_vsites(gmx_domdec_t *dd, int n_intercg_vsite)
+{
+ if (debug)
+ {
+ fprintf(debug, "Begin init_domdec_vsites\n");
+ }
+
+ /* Use a hash table for the global to local index.
+ * The number of keys is a rough estimate, it will be optimized later.
+ */
+ dd->ga2la_vsite = gmx_hash_init(std::min(n_intercg_vsite/20,
+ n_intercg_vsite/(2*dd->nnodes)));
+
+ dd->vsite_comm = specat_comm_init(1);
+}
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2005,2006,2007,2008,2009,2010,2012,2013,2014, 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.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ *
+ * \brief This file declares functions for domdec to use
+ * while managing virtual sites.
+ *
+ * \author Berk Hess <hess@kth.se>
+ * \ingroup module_domdec
+ */
+
+#ifndef GMX_DOMDEC_DOMDEC_VSITE_H
+#define GMX_DOMDEC_DOMDEC_VSITE_H
+
+#include "gromacs/legacyheaders/types/commrec_fwd.h"
+#include "gromacs/topology/idef.h"
+#include "gromacs/topology/topology.h"
+
+/*! \brief Clears the local indices for the virtual site communication setup */
+void dd_clear_local_vsite_indices(gmx_domdec_t *dd);
+
+/*! \brief Sets up communication and atom indices for all local vsites */
+int dd_make_local_vsites(gmx_domdec_t *dd, int at_start, t_ilist *lil);
+
+/*! \brief Initializes the data structures for virtual site communication */
+void init_domdec_vsites(gmx_domdec_t *dd, int n_intercg_vsite);
+
+#endif