From: Berk Hess Date: Mon, 1 Dec 2014 21:40:40 +0000 (+0100) Subject: Make domdec a proper module X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=a56a53084c3929172f20471cc7e2c7cf88ad29ba;p=alexxy%2Fgromacs.git Make domdec a proper module Now has * Doxygen group definition * initial Doxygen comments * correct #include guards * no export "C" for header used only in C++ code * all relevant current files are declared as such in the module * new domdec_constraints.h and domdec_internal.h * split off domdec_specatomcomm.cpp and domdec_vsite.c from domdec_constraints.cpp Change-Id: I7ed79c1654850df31e242aa1e907afd124e3ae1b --- diff --git a/src/gromacs/domdec/domdec.cpp b/src/gromacs/domdec/domdec.cpp index f70b5c902d..fa932d308f 100644 --- a/src/gromacs/domdec/domdec.cpp +++ b/src/gromacs/domdec/domdec.cpp @@ -107,6 +107,10 @@ #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) @@ -2093,7 +2097,7 @@ void write_dd_pdb(const char *fn, gmx_int64_t step, const char *title, 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; @@ -2130,11 +2134,11 @@ real dd_cutoff_mbody(gmx_domdec_t *dd) 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); } @@ -2339,22 +2343,6 @@ void get_pme_nnodes(const gmx_domdec_t *dd, } } -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) { diff --git a/src/gromacs/domdec/domdec.h b/src/gromacs/domdec/domdec.h index d8b9e6373f..645dbcc951 100644 --- a/src/gromacs/domdec/domdec.h +++ b/src/gromacs/domdec/domdec.h @@ -32,9 +32,31 @@ * 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 + * + */ -#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 + * \inlibraryapi + * \ingroup module_domdec + */ + +#ifndef GMX_DOMDEC_DOMDEC_H +#define GMX_DOMDEC_DOMDEC_H #include @@ -60,57 +82,55 @@ 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, @@ -124,92 +144,111 @@ init_domain_decomposition(FILE *fplog, 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, @@ -229,62 +268,52 @@ void dd_partition_system(FILE *fplog, 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, @@ -293,56 +322,64 @@ void dd_make_local_top(gmx_domdec_t *dd, gmx_domdec_zones_t *zones, 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); @@ -351,4 +388,4 @@ void set_ddbox_cr(t_commrec *cr, ivec *dd_nc, } #endif -#endif /* _domdec_h */ +#endif diff --git a/src/gromacs/domdec/domdec_box.cpp b/src/gromacs/domdec/domdec_box.cpp index 8e6e5829fb..8dc0e518f8 100644 --- a/src/gromacs/domdec/domdec_box.cpp +++ b/src/gromacs/domdec/domdec_box.cpp @@ -33,6 +33,15 @@ * 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 + * \ingroup module_domdec + */ + #include "gmxpre.h" #include "gromacs/domdec/domdec.h" @@ -45,6 +54,7 @@ #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) { @@ -116,6 +126,7 @@ static void calc_cgcm_av_stddev(t_block *cgs, int n, rvec *x, rvec av, rvec stdd } } +/*! \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; @@ -225,6 +236,7 @@ static void set_tric_dir(ivec *dd_nc, gmx_ddbox_t *ddbox, matrix box) } } +/*! \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, diff --git a/src/gromacs/domdec/domdec_constraints.cpp b/src/gromacs/domdec/domdec_constraints.cpp index dacbc47be7..84141fc5c7 100644 --- a/src/gromacs/domdec/domdec_constraints.cpp +++ b/src/gromacs/domdec/domdec_constraints.cpp @@ -33,14 +33,24 @@ * 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 + * \ingroup module_domdec + */ + #include "gmxpre.h" +#include "domdec_constraints.h" + #include #include #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" @@ -54,373 +64,29 @@ #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) { @@ -430,14 +96,6 @@ void dd_move_x_constraints(gmx_domdec_t *dd, matrix box, } } -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) @@ -476,272 +134,7 @@ void dd_clear_local_vsite_indices(gmx_domdec_t *dd) } } -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, @@ -839,6 +232,7 @@ static void walk_out(int con, int con_offset, int a, int offset, int nrec, } } +/*! \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, @@ -941,6 +335,7 @@ static void atoms_to_settles(gmx_domdec_t *dd, 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, @@ -1263,99 +658,6 @@ int dd_make_local_constraints(gmx_domdec_t *dd, int at_start, 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) { @@ -1405,19 +707,3 @@ void init_domdec_constraints(gmx_domdec_t *dd, 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); -} diff --git a/src/gromacs/domdec/domdec_constraints.h b/src/gromacs/domdec/domdec_constraints.h new file mode 100644 index 0000000000..d59e2ceb93 --- /dev/null +++ b/src/gromacs/domdec/domdec_constraints.h @@ -0,0 +1,66 @@ +/* + * 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 + * \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 diff --git a/src/gromacs/domdec/domdec_internal.h b/src/gromacs/domdec/domdec_internal.h new file mode 100644 index 0000000000..fb560877c8 --- /dev/null +++ b/src/gromacs/domdec/domdec_internal.h @@ -0,0 +1,55 @@ +/* + * 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 + * \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 diff --git a/src/gromacs/domdec/domdec_network.cpp b/src/gromacs/domdec/domdec_network.cpp index 85c24e6094..3710f598b6 100644 --- a/src/gromacs/domdec/domdec_network.cpp +++ b/src/gromacs/domdec/domdec_network.cpp @@ -33,6 +33,15 @@ * 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 + * \ingroup module_domdec + */ + #include "gmxpre.h" #include "domdec_network.h" @@ -45,6 +54,7 @@ #include "gromacs/utility/gmxmpi.h" +/*! \brief Returns the MPI rank of the domain decomposition master rank */ #define DDMASTERRANK(dd) (dd->masterrank) diff --git a/src/gromacs/domdec/domdec_network.h b/src/gromacs/domdec/domdec_network.h index 5700800764..efadf4327f 100644 --- a/src/gromacs/domdec/domdec_network.h +++ b/src/gromacs/domdec/domdec_network.h @@ -33,22 +33,33 @@ * 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 + * \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, @@ -56,9 +67,10 @@ 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, @@ -66,9 +78,10 @@ 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, @@ -77,9 +90,10 @@ 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, @@ -96,33 +110,40 @@ 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 */ diff --git a/src/gromacs/domdec/domdec_setup.cpp b/src/gromacs/domdec/domdec_setup.cpp index 087b471ac3..8d0bf79b7b 100644 --- a/src/gromacs/domdec/domdec_setup.cpp +++ b/src/gromacs/domdec/domdec_setup.cpp @@ -33,6 +33,15 @@ * 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 + * \ingroup module_domdec + */ + #include "gmxpre.h" #include @@ -49,9 +58,16 @@ #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; @@ -84,6 +100,7 @@ static int factorize(int n, int **fac, int **mfac) 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; @@ -96,6 +113,7 @@ static gmx_bool largest_divisor(int n) return ldiv; } +/*! \brief Compute largest common divisor of \p n1 and \b n2 */ static int lcd(int n1, int n2) { int d, i; @@ -112,50 +130,53 @@ static int lcd(int n1, int n2) 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(std::pow(nnodes-npme, 1.0/3.0) + 0.5); + npp_root3 = static_cast(std::pow(ntot - npme, 1.0/3.0) + 0.5); npme_root2 = static_cast(std::sqrt(static_cast(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; @@ -167,59 +188,59 @@ static int guess_npme(FILE *fplog, gmx_mtop_t *mtop, t_inputrec *ir, matrix box, 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; } @@ -230,17 +251,18 @@ static int guess_npme(FILE *fplog, gmx_mtop_t *mtop, t_inputrec *ir, matrix box, 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; @@ -284,15 +306,23 @@ real comm_box_frac(ivec dd_nc, real cutoff, gmx_ddbox_t *ddbox) 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; @@ -300,9 +330,11 @@ static float comm_pme_cost_vol(int npme, int a, int b, int c) 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, @@ -462,12 +494,6 @@ static float comm_cost_est(real limit, real cutoff, } } - /* 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); @@ -498,6 +524,7 @@ static float comm_cost_est(real limit, real cutoff, 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, @@ -560,6 +587,7 @@ static void assign_factors(gmx_domdec_t *dd, } } +/*! \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, diff --git a/src/gromacs/domdec/domdec_specatomcomm.cpp b/src/gromacs/domdec/domdec_specatomcomm.cpp new file mode 100644 index 0000000000..abc52ba315 --- /dev/null +++ b/src/gromacs/domdec/domdec_specatomcomm.cpp @@ -0,0 +1,629 @@ +/* + * 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 + * \ingroup module_domdec + */ + +#include "gmxpre.h" + +#include "domdec_specatomcomm.h" + +#include + +#include + +#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; +} diff --git a/src/gromacs/domdec/domdec_specatomcomm.h b/src/gromacs/domdec/domdec_specatomcomm.h new file mode 100644 index 0000000000..d60a98f5b5 --- /dev/null +++ b/src/gromacs/domdec/domdec_specatomcomm.h @@ -0,0 +1,131 @@ +/* + * 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 + * \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 diff --git a/src/gromacs/domdec/domdec_topology.cpp b/src/gromacs/domdec/domdec_topology.cpp index 32f19f1eb0..39e40f66fe 100644 --- a/src/gromacs/domdec/domdec_topology.cpp +++ b/src/gromacs/domdec/domdec_topology.cpp @@ -33,6 +33,16 @@ * 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 + * \ingroup module_domdec + */ + #include "gmxpre.h" #include @@ -58,50 +68,62 @@ #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); @@ -116,7 +138,7 @@ static int nral_rt(int 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) { @@ -127,6 +149,7 @@ static gmx_bool dd_check_ftype(int ftype, gmx_bool bBCheck, (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); @@ -139,10 +162,11 @@ static void print_error_header(FILE *fplog, char *moltypename, int nprint) 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) @@ -276,6 +300,7 @@ static void print_missing_interactions_mb(FILE *fplog, t_commrec *cr, 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) { @@ -399,18 +424,19 @@ void dd_print_missing_interactions(FILE *fplog, t_commrec *cr, int local_count, } 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) @@ -438,7 +464,9 @@ static void global_atomnr_to_moltype_ind(gmx_reverse_top_t *rt, int i_gl, *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; @@ -468,6 +496,7 @@ static int count_excls(t_block *cgs, t_blocka *excls, int *n_intercg_excl) 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, @@ -570,12 +599,13 @@ static int low_make_reverse_ilist(t_ilist *il_mt, t_atom *atom, 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; @@ -609,12 +639,14 @@ static int make_reverse_ilist(gmx_moltype_t *molt, 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, @@ -688,20 +720,15 @@ static gmx_reverse_top_t *make_reverse_top(gmx_mtop_t *mtop, gmx_bool bFE, } 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; } @@ -790,7 +817,9 @@ void dd_make_reverse_top(FILE *fplog, } } -/* 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 @@ -844,6 +873,7 @@ add_ifunc_for_vsites(t_iatom *tiatoms, gmx_ga2la_t ga2la, } } +/*! \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; @@ -862,6 +892,7 @@ static gmx_inline void add_ifunc(int nral, t_iatom *tiatoms, t_ilist *il) 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) @@ -907,6 +938,7 @@ static void add_posres(int mol, int a_mol, const gmx_molblock_t *molb, 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) @@ -944,6 +976,7 @@ static void add_fbposres(int mol, int a_mol, const gmx_molblock_t *molb, 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, @@ -1035,6 +1068,7 @@ static void add_vsite(gmx_ga2la_t ga2la, int *index, int *rtil, } } +/*! \brief Update the local atom to local charge group index */ static void make_la2lc(gmx_domdec_t *dd) { int *cgindex, *la2lc, cg, a; @@ -1058,6 +1092,7 @@ static void make_la2lc(gmx_domdec_t *dd) } } +/*! \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; @@ -1074,16 +1109,16 @@ static real dd_dist2(t_pbc *pbc_null, rvec *cg_cm, const int *la2lc, int i, int 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) { @@ -1095,42 +1130,42 @@ static void combine_blocka(t_blocka *dest, const t_blocka *src, int nsrc) 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) @@ -1139,6 +1174,9 @@ static void combine_idef(t_idef *dest, const t_idef *src, int nsrc, 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) @@ -1154,9 +1192,12 @@ static void combine_idef(t_idef *dest, const t_idef *src, int nsrc, } } - 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]; @@ -1166,7 +1207,7 @@ static void combine_idef(t_idef *dest, const t_idef *src, int nsrc, 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]; } } @@ -1178,6 +1219,8 @@ static void combine_idef(t_idef *dest, const t_idef *src, int nsrc, /* 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) { @@ -1185,25 +1228,26 @@ static void combine_idef(t_idef *dest, const t_idef *src, int nsrc, 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. */ @@ -1430,6 +1474,7 @@ static int make_bondeds_zone(gmx_domdec_t *dd, 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) { @@ -1444,6 +1489,7 @@ static void set_no_exclusions_zone(gmx_domdec_t *dd, gmx_domdec_zones_t *zones, } } +/*! \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, @@ -1582,6 +1628,7 @@ static int make_exclusions_zone(gmx_domdec_t *dd, gmx_domdec_zones_t *zones, 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) @@ -1591,6 +1638,7 @@ static void check_alloc_index(t_blocka *ba, int nindex_max) } } +/*! \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) { @@ -1603,10 +1651,11 @@ static void check_exclusions_alloc(gmx_domdec_t *dd, gmx_domdec_zones_t *zones, 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) { @@ -1632,6 +1681,7 @@ static void finish_local_exclusions(gmx_domdec_t *dd, gmx_domdec_zones_t *zones, } } +/*! \brief Clear a t_idef struct */ static void clear_idef(t_idef *idef) { int ftype; @@ -1643,6 +1693,7 @@ static void clear_idef(t_idef *idef) } } +/*! \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, @@ -1720,7 +1771,7 @@ static int make_local_bondeds_excls(gmx_domdec_t *dd, } else { - idef_t = &rt->idef_thread[thread]; + idef_t = &rt->th_work[thread].idef; clear_idef(idef_t); } @@ -1733,8 +1784,8 @@ static int make_local_bondeds_excls(gmx_domdec_t *dd, } 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 @@ -1743,7 +1794,7 @@ static int make_local_bondeds_excls(gmx_domdec_t *dd, vsite_pbc_nalloc = NULL; } - rt->nbonded_thread[thread] = + rt->th_work[thread].nbonded = make_bondeds_zone(dd, zones, mtop->molblock, bRCheckMB, rcheck, bRCheck2B, rc2, @@ -1761,12 +1812,12 @@ static int make_local_bondeds_excls(gmx_domdec_t *dd, } 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, @@ -1778,25 +1829,24 @@ static int make_local_bondeds_excls(gmx_domdec_t *dd, 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; } } } @@ -1984,6 +2034,7 @@ void dd_init_local_state(gmx_domdec_t *dd, 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; @@ -2013,6 +2064,7 @@ static void check_link(t_blocka *link, int cg_gl, int cg_gl_j) } } +/*! \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; @@ -2039,7 +2091,7 @@ t_blocka *make_charge_group_links(gmx_mtop_t *mtop, gmx_domdec_t *dd, 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; @@ -2167,6 +2219,7 @@ t_blocka *make_charge_group_links(gmx_mtop_t *mtop, gmx_domdec_t *dd, 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, @@ -2244,6 +2297,7 @@ static void bonded_cg_distance_mol(gmx_moltype_t *molt, int *at2cg, *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, @@ -2286,6 +2340,7 @@ static void get_cgcm_mol(gmx_moltype_t *molt, gmx_ffparams_t *ffparams, 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; diff --git a/src/gromacs/domdec/domdec_vsite.cpp b/src/gromacs/domdec/domdec_vsite.cpp new file mode 100644 index 0000000000..e9d063f382 --- /dev/null +++ b/src/gromacs/domdec/domdec_vsite.cpp @@ -0,0 +1,193 @@ +/* + * 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 + * \ingroup module_domdec + */ + +#include "gmxpre.h" + +#include "domdec_vsite.h" + +#include + +#include + +#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); +} diff --git a/src/gromacs/domdec/domdec_vsite.h b/src/gromacs/domdec/domdec_vsite.h new file mode 100644 index 0000000000..7c42b391f3 --- /dev/null +++ b/src/gromacs/domdec/domdec_vsite.h @@ -0,0 +1,60 @@ +/* + * 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 + * \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