Make domdec a proper module
authorBerk Hess <hess@kth.se>
Mon, 1 Dec 2014 21:40:40 +0000 (22:40 +0100)
committerGerrit Code Review <gerrit@gerrit.gromacs.org>
Fri, 5 Dec 2014 17:38:16 +0000 (18:38 +0100)
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

14 files changed:
src/gromacs/domdec/domdec.cpp
src/gromacs/domdec/domdec.h
src/gromacs/domdec/domdec_box.cpp
src/gromacs/domdec/domdec_constraints.cpp
src/gromacs/domdec/domdec_constraints.h [new file with mode: 0644]
src/gromacs/domdec/domdec_internal.h [new file with mode: 0644]
src/gromacs/domdec/domdec_network.cpp
src/gromacs/domdec/domdec_network.h
src/gromacs/domdec/domdec_setup.cpp
src/gromacs/domdec/domdec_specatomcomm.cpp [new file with mode: 0644]
src/gromacs/domdec/domdec_specatomcomm.h [new file with mode: 0644]
src/gromacs/domdec/domdec_topology.cpp
src/gromacs/domdec/domdec_vsite.cpp [new file with mode: 0644]
src/gromacs/domdec/domdec_vsite.h [new file with mode: 0644]

index f70b5c902d483c1398212b4ff11c54ce97efa2bd..fa932d308f8559693bd6a773d79ef5492562edc6 100644 (file)
 #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)
 {
index d8b9e6373f496f61380dd1793f73d448e94142c3..645dbcc951e42a61247ea01596f85829281f5f37 100644 (file)
  * To help us fund GROMACS development, we humbly ask that you cite
  * the research papers on the package. Check out http://www.gromacs.org.
  */
+/*! \defgroup module_domdec Spatial domain decomposition (for parallelization over MPI)
+ * \ingroup group_mdrun
+ *
+ * \brief Manages the decomposition of the simulation volume over MPI
+ * ranks to try to distribute work evenly with minimal communication
+ * overheads.
+ *
+ * \todo Get domdec stuff out of legacyheaders/types/commrec.h
+ *
+ * \author Berk Hess <hess@kth.se>
+ *
+ */
 
-#ifndef _domdec_h
-#define _domdec_h
+/*! \libinternal \file
+ *
+ * \brief This file declares functions for mdrun to call to manage the
+ * details of its domain decomposition.
+ *
+ * \author Berk Hess <hess@kth.se>
+ * \inlibraryapi
+ * \ingroup module_domdec
+ */
+
+#ifndef GMX_DOMDEC_DOMDEC_H
+#define GMX_DOMDEC_DOMDEC_H
 
 #include <stdio.h>
 
 extern "C" {
 #endif
 
-int ddglatnr(gmx_domdec_t *dd, int i);
-/* Returns the global topology atom number belonging to local atom index i.
- * This function is intended for writing ascii output
+/*! \brief Returns the global topology atom number belonging to local atom index i.
+ *
+ * This function is intended for writing ASCII output
  * and returns atom numbers starting at 1.
  * When dd=NULL returns i+1.
  */
+int ddglatnr(gmx_domdec_t *dd, int i);
 
+/*! \brief Return a block struct for the charge groups of the whole system */
 t_block *dd_charge_groups_global(gmx_domdec_t *dd);
-/* Return a block struct for the charge groups of the whole system */
-
-gmx_bool dd_filled_nsgrid_home(gmx_domdec_t *dd);
-/* Is the ns grid already filled with the home particles? */
 
-void dd_store_state(gmx_domdec_t *dd, t_state *state);
-/* Store the global cg indices of the home cgs in state,
- * so it can be reset, even after a new DD partitioning.
+/*! \brief Store the global cg indices of the home cgs in state,
+ *
+ * This means it can be reset, even after a new DD partitioning.
  */
+void dd_store_state(gmx_domdec_t *dd, t_state *state);
 
+/*! \brief Returns a pointer to the gmx_domdec_zones_t struct */
 gmx_domdec_zones_t *domdec_zones(gmx_domdec_t *dd);
 
+/*! \brief Sets the j-charge-group range for i-charge-group \p icg */
 void dd_get_ns_ranges(gmx_domdec_t *dd, int icg,
                       int *jcg0, int *jcg1, ivec shift0, ivec shift1);
 
+/*! \brief Returns the atom range in the local state for atoms involved in virtual sites */
 int dd_natoms_vsite(gmx_domdec_t *dd);
 
+/*! \brief Sets the atom range for atom in the local state for atoms received in constraints communication */
 void dd_get_constraint_range(gmx_domdec_t *dd,
                              int *at_start, int *at_end);
 
-real dd_cutoff_mbody(gmx_domdec_t *dd);
-
-real dd_cutoff_twobody(gmx_domdec_t *dd);
-
+/*! \brief Get the number of PME nodes along x and y, can be called with dd=NULL */
 void get_pme_nnodes(const gmx_domdec_t *dd,
                     int *npmenodes_x, int *npmenodes_y);
-/* Get the number of PME nodes along x and y, can be called with dd=NULL */
-
-gmx_bool gmx_pmeonlynode(t_commrec *cr, int nodeid);
-/* Return if nodeid in cr->mpi_comm_mysim is a PME-only node */
 
+/*! \brief Returns the set of DD nodes that communicate with pme node cr->nodeid */
 void get_pme_ddnodes(t_commrec *cr, int pmenodeid,
                      int *nmy_ddnodes, int **my_ddnodes, int *node_peer);
-/* Returns the set of DD nodes that communicate with pme node cr->nodeid */
 
+/*! \brief Returns the maximum shift for coordinate communication in PME, dim x */
 int dd_pme_maxshift_x(gmx_domdec_t *dd);
-/* Returns the maximum shift for coordinate communication in PME, dim x */
 
+/*! \brief Returns the maximum shift for coordinate communication in PME, dim y */
 int dd_pme_maxshift_y(gmx_domdec_t *dd);
-/* Returns the maximum shift for coordinate communication in PME, dim y */
 
+/*! \brief Generates the MPI communicators for domain decomposition */
 void make_dd_communicators(FILE *fplog, t_commrec *cr, int dd_node_order);
 
+/*! \brief Initialized the domain decomposition, chooses the DD grid and PME ranks */
 gmx_domdec_t *
 init_domain_decomposition(FILE *fplog,
                           t_commrec *cr,
@@ -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
index 8e6e5829fb22771175e3940768bfecb87faac4c9..8dc0e518f80a75decf8f02b0333e6da8ea6735c3 100644 (file)
  * the research papers on the package. Check out http://www.gromacs.org.
  */
 
+/*! \internal \file
+ *
+ * \brief This file defines functions used by the domdec module
+ * for (bounding) box and pbc information generation.
+ *
+ * \author Berk Hess <hess@kth.se>
+ * \ingroup module_domdec
+ */
+
 #include "gmxpre.h"
 
 #include "gromacs/domdec/domdec.h"
@@ -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,
index dacbc47be7f91fefd134b9c632ed7122cbcfdaee..84141fc5c72180891db135421a09905e203589ad 100644 (file)
  * the research papers on the package. Check out http://www.gromacs.org.
  */
 
+/*! \internal \file
+ *
+ * \brief This file implements functions for domdec to use
+ * while managing inter-atomic constraints.
+ *
+ * \author Berk Hess <hess@kth.se>
+ * \ingroup module_domdec
+ */
+
 #include "gmxpre.h"
 
+#include "domdec_constraints.h"
+
 #include <assert.h>
 
 #include <algorithm>
 
 #include "gromacs/domdec/domdec.h"
-#include "gromacs/domdec/domdec_network.h"
 #include "gromacs/legacyheaders/constr.h"
 #include "gromacs/legacyheaders/gmx_ga2la.h"
 #include "gromacs/legacyheaders/gmx_hash.h"
 #include "gromacs/utility/gmxassert.h"
 #include "gromacs/utility/smalloc.h"
 
-typedef struct {
-    int  nsend;
-    int *a;
-    int  a_nalloc;
-    int  nrecv;
-} gmx_specatsend_t;
-
-typedef struct {
-    int *ind;
-    int  nalloc;
-    int  n;
-} ind_req_t;
-
-typedef struct gmx_domdec_specat_comm {
-    /* The number of indices to receive during the setup */
-    int              nreq[DIM][2][2];
-    /* The atoms to send */
-    gmx_specatsend_t spas[DIM][2];
-    gmx_bool        *bSendAtom;
-    int              bSendAtom_nalloc;
-    /* Send buffers */
-    int             *ibuf;
-    int              ibuf_nalloc;
-    rvec            *vbuf;
-    int              vbuf_nalloc;
-    rvec            *vbuf2;
-    int              vbuf2_nalloc;
-    /* The range in the local buffer(s) for received atoms */
-    int              at_start;
-    int              at_end;
-
-    /* The atom indices we need from the surrounding cells.
-     * We can gather the indices over nthread threads.
-     */
-    int        nthread;
-    ind_req_t *ireq;
-} gmx_domdec_specat_comm_t;
+#include "domdec_specatomcomm.h"
 
+/*! \brief Struct used during constraint setup with domain decomposition */
 typedef struct gmx_domdec_constraints {
-    int       *molb_con_offset;
-    int       *molb_ncon_mol;
-    /* The fully local and connected constraints */
-    int        ncon;
+    //! @cond Doxygen_Suppress
+    int       *molb_con_offset; /**< Offset in the constraint array for each molblock */
+    int       *molb_ncon_mol;   /**< The number of constraints per molecule for each molblock */
+
+    int        ncon;            /**< The fully local and conneced constraints */
     /* The global constraint number, only required for clearing gc_req */
-    int       *con_gl;
-    int       *con_nlocat;
-    int        con_nalloc;
-    /* Boolean that tells if a global constraint index has been requested */
-    char      *gc_req;
-    /* Global to local communicated constraint atom only index */
-    gmx_hash_t ga2la;
+    int       *con_gl;          /**< Global constraint indices for local constraints */
+    int       *con_nlocat;      /**< Number of local atoms (2/1/0) for each constraint */
+    int        con_nalloc;      /**< Allocation size for \p con_gl and \p con_nlocat */
+
+    char      *gc_req;          /**< Boolean that tells if a global constraint index has been requested; note: size global #constraints */
+    gmx_hash_t ga2la;           /**< Global to local communicated constraint atom only index */
 
     /* Multi-threading stuff */
-    int      nthread;
-    t_ilist *ils;
+    int      nthread;           /**< Number of threads used for DD constraint setup */
+    t_ilist *ils;               /**< Constraint ilist working arrays, size \p nthread */
+    //! @endcond
 } gmx_domdec_constraints_t;
 
-
-static void dd_move_f_specat(gmx_domdec_t *dd, gmx_domdec_specat_comm_t *spac,
-                             rvec *f, rvec *fshift)
-{
-    gmx_specatsend_t *spas;
-    rvec             *vbuf;
-    int               n, n0, n1, d, dim, dir, i;
-    ivec              vis;
-    int               is;
-    gmx_bool          bPBC, bScrew;
-
-    n = spac->at_end;
-    for (d = dd->ndim-1; d >= 0; d--)
-    {
-        dim = dd->dim[d];
-        if (dd->nc[dim] > 2)
-        {
-            /* Pulse the grid forward and backward */
-            spas = spac->spas[d];
-            n0   = spas[0].nrecv;
-            n1   = spas[1].nrecv;
-            n   -= n1 + n0;
-            vbuf = spac->vbuf;
-            /* Send and receive the coordinates */
-            dd_sendrecv2_rvec(dd, d,
-                              f+n+n1, n0, vbuf, spas[0].nsend,
-                              f+n, n1, vbuf+spas[0].nsend, spas[1].nsend);
-            for (dir = 0; dir < 2; dir++)
-            {
-                bPBC   = ((dir == 0 && dd->ci[dim] == 0) ||
-                          (dir == 1 && dd->ci[dim] == dd->nc[dim]-1));
-                bScrew = (bPBC && dd->bScrewPBC && dim == XX);
-
-                spas = &spac->spas[d][dir];
-                /* Sum the buffer into the required forces */
-                if (!bPBC || (!bScrew && fshift == NULL))
-                {
-                    for (i = 0; i < spas->nsend; i++)
-                    {
-                        rvec_inc(f[spas->a[i]], *vbuf);
-                        vbuf++;
-                    }
-                }
-                else
-                {
-                    clear_ivec(vis);
-                    vis[dim] = (dir == 0 ? 1 : -1);
-                    is       = IVEC2IS(vis);
-                    if (!bScrew)
-                    {
-                        /* Sum and add to shift forces */
-                        for (i = 0; i < spas->nsend; i++)
-                        {
-                            rvec_inc(f[spas->a[i]], *vbuf);
-                            rvec_inc(fshift[is], *vbuf);
-                            vbuf++;
-                        }
-                    }
-                    else
-                    {
-                        /* Rotate the forces */
-                        for (i = 0; i < spas->nsend; i++)
-                        {
-                            f[spas->a[i]][XX] += (*vbuf)[XX];
-                            f[spas->a[i]][YY] -= (*vbuf)[YY];
-                            f[spas->a[i]][ZZ] -= (*vbuf)[ZZ];
-                            if (fshift)
-                            {
-                                rvec_inc(fshift[is], *vbuf);
-                            }
-                            vbuf++;
-                        }
-                    }
-                }
-            }
-        }
-        else
-        {
-            /* Two cells, so we only need to communicate one way */
-            spas = &spac->spas[d][0];
-            n   -= spas->nrecv;
-            /* Send and receive the coordinates */
-            dd_sendrecv_rvec(dd, d, dddirForward,
-                             f+n, spas->nrecv, spac->vbuf, spas->nsend);
-            /* Sum the buffer into the required forces */
-            if (dd->bScrewPBC && dim == XX &&
-                (dd->ci[dim] == 0 ||
-                 dd->ci[dim] == dd->nc[dim]-1))
-            {
-                for (i = 0; i < spas->nsend; i++)
-                {
-                    /* Rotate the force */
-                    f[spas->a[i]][XX] += spac->vbuf[i][XX];
-                    f[spas->a[i]][YY] -= spac->vbuf[i][YY];
-                    f[spas->a[i]][ZZ] -= spac->vbuf[i][ZZ];
-                }
-            }
-            else
-            {
-                for (i = 0; i < spas->nsend; i++)
-                {
-                    rvec_inc(f[spas->a[i]], spac->vbuf[i]);
-                }
-            }
-        }
-    }
-}
-
-void dd_move_f_vsites(gmx_domdec_t *dd, rvec *f, rvec *fshift)
-{
-    if (dd->vsite_comm)
-    {
-        dd_move_f_specat(dd, dd->vsite_comm, f, fshift);
-    }
-}
-
-void dd_clear_f_vsites(gmx_domdec_t *dd, rvec *f)
-{
-    int i;
-
-    if (dd->vsite_comm)
-    {
-        for (i = dd->vsite_comm->at_start; i < dd->vsite_comm->at_end; i++)
-        {
-            clear_rvec(f[i]);
-        }
-    }
-}
-
-static void dd_move_x_specat(gmx_domdec_t *dd, gmx_domdec_specat_comm_t *spac,
-                             matrix box,
-                             rvec *x0,
-                             rvec *x1, gmx_bool bX1IsCoord)
-{
-    gmx_specatsend_t *spas;
-    rvec             *x, *vbuf, *rbuf;
-    int               nvec, v, n, nn, ns0, ns1, nr0, nr1, nr, d, dim, dir, i;
-    gmx_bool          bPBC, bScrew = FALSE;
-    rvec              shift = {0, 0, 0};
-
-    nvec = 1;
-    if (x1 != NULL)
-    {
-        nvec++;
-    }
-
-    n = spac->at_start;
-    for (d = 0; d < dd->ndim; d++)
-    {
-        dim = dd->dim[d];
-        if (dd->nc[dim] > 2)
-        {
-            /* Pulse the grid forward and backward */
-            vbuf = spac->vbuf;
-            for (dir = 0; dir < 2; dir++)
-            {
-                if (dir == 0 && dd->ci[dim] == 0)
-                {
-                    bPBC   = TRUE;
-                    bScrew = (dd->bScrewPBC && dim == XX);
-                    copy_rvec(box[dim], shift);
-                }
-                else if (dir == 1 && dd->ci[dim] == dd->nc[dim]-1)
-                {
-                    bPBC   = TRUE;
-                    bScrew = (dd->bScrewPBC && dim == XX);
-                    for (i = 0; i < DIM; i++)
-                    {
-                        shift[i] = -box[dim][i];
-                    }
-                }
-                else
-                {
-                    bPBC   = FALSE;
-                    bScrew = FALSE;
-                }
-                spas = &spac->spas[d][dir];
-                for (v = 0; v < nvec; v++)
-                {
-                    x = (v == 0 ? x0 : x1);
-                    /* Copy the required coordinates to the send buffer */
-                    if (!bPBC || (v == 1 && !bX1IsCoord))
-                    {
-                        /* Only copy */
-                        for (i = 0; i < spas->nsend; i++)
-                        {
-                            copy_rvec(x[spas->a[i]], *vbuf);
-                            vbuf++;
-                        }
-                    }
-                    else if (!bScrew)
-                    {
-                        /* Shift coordinates */
-                        for (i = 0; i < spas->nsend; i++)
-                        {
-                            rvec_add(x[spas->a[i]], shift, *vbuf);
-                            vbuf++;
-                        }
-                    }
-                    else
-                    {
-                        /* Shift and rotate coordinates */
-                        for (i = 0; i < spas->nsend; i++)
-                        {
-                            (*vbuf)[XX] =               x[spas->a[i]][XX] + shift[XX];
-                            (*vbuf)[YY] = box[YY][YY] - x[spas->a[i]][YY] + shift[YY];
-                            (*vbuf)[ZZ] = box[ZZ][ZZ] - x[spas->a[i]][ZZ] + shift[ZZ];
-                            vbuf++;
-                        }
-                    }
-                }
-            }
-            /* Send and receive the coordinates */
-            spas = spac->spas[d];
-            ns0  = spas[0].nsend;
-            nr0  = spas[0].nrecv;
-            ns1  = spas[1].nsend;
-            nr1  = spas[1].nrecv;
-            if (nvec == 1)
-            {
-                dd_sendrecv2_rvec(dd, d,
-                                  spac->vbuf+ns0, ns1, x0+n, nr1,
-                                  spac->vbuf, ns0, x0+n+nr1, nr0);
-            }
-            else
-            {
-                /* Communicate both vectors in one buffer */
-                rbuf = spac->vbuf2;
-                dd_sendrecv2_rvec(dd, d,
-                                  spac->vbuf+2*ns0, 2*ns1, rbuf, 2*nr1,
-                                  spac->vbuf, 2*ns0, rbuf+2*nr1, 2*nr0);
-                /* Split the buffer into the two vectors */
-                nn = n;
-                for (dir = 1; dir >= 0; dir--)
-                {
-                    nr = spas[dir].nrecv;
-                    for (v = 0; v < 2; v++)
-                    {
-                        x = (v == 0 ? x0 : x1);
-                        for (i = 0; i < nr; i++)
-                        {
-                            copy_rvec(*rbuf, x[nn+i]);
-                            rbuf++;
-                        }
-                    }
-                    nn += nr;
-                }
-            }
-            n += nr0 + nr1;
-        }
-        else
-        {
-            spas = &spac->spas[d][0];
-            /* Copy the required coordinates to the send buffer */
-            vbuf = spac->vbuf;
-            for (v = 0; v < nvec; v++)
-            {
-                x = (v == 0 ? x0 : x1);
-                if (dd->bScrewPBC && dim == XX &&
-                    (dd->ci[XX] == 0 || dd->ci[XX] == dd->nc[XX]-1))
-                {
-                    /* Here we only perform the rotation, the rest of the pbc
-                     * is handled in the constraint or viste routines.
-                     */
-                    for (i = 0; i < spas->nsend; i++)
-                    {
-                        (*vbuf)[XX] =               x[spas->a[i]][XX];
-                        (*vbuf)[YY] = box[YY][YY] - x[spas->a[i]][YY];
-                        (*vbuf)[ZZ] = box[ZZ][ZZ] - x[spas->a[i]][ZZ];
-                        vbuf++;
-                    }
-                }
-                else
-                {
-                    for (i = 0; i < spas->nsend; i++)
-                    {
-                        copy_rvec(x[spas->a[i]], *vbuf);
-                        vbuf++;
-                    }
-                }
-            }
-            /* Send and receive the coordinates */
-            if (nvec == 1)
-            {
-                dd_sendrecv_rvec(dd, d, dddirBackward,
-                                 spac->vbuf, spas->nsend, x0+n, spas->nrecv);
-            }
-            else
-            {
-                /* Communicate both vectors in one buffer */
-                rbuf = spac->vbuf2;
-                dd_sendrecv_rvec(dd, d, dddirBackward,
-                                 spac->vbuf, 2*spas->nsend, rbuf, 2*spas->nrecv);
-                /* Split the buffer into the two vectors */
-                nr = spas[0].nrecv;
-                for (v = 0; v < 2; v++)
-                {
-                    x = (v == 0 ? x0 : x1);
-                    for (i = 0; i < nr; i++)
-                    {
-                        copy_rvec(*rbuf, x[n+i]);
-                        rbuf++;
-                    }
-                }
-            }
-            n += spas->nrecv;
-        }
-    }
-}
-
 void dd_move_x_constraints(gmx_domdec_t *dd, matrix box,
                            rvec *x0, rvec *x1, gmx_bool bX1IsCoord)
 {
@@ -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 (file)
index 0000000..d59e2ce
--- /dev/null
@@ -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 <hess@kth.se>
+ * \ingroup module_domdec
+ */
+
+#ifndef GMX_DOMDEC_DOMDEC_CONSTRAINTS_H
+#define GMX_DOMDEC_DOMDEC_CONSTRAINTS_H
+
+#include "gromacs/legacyheaders/types/commrec_fwd.h"
+#include "gromacs/legacyheaders/types/constr.h"
+#include "gromacs/topology/idef.h"
+#include "gromacs/topology/topology.h"
+
+/*! \brief Clears the local indices for the constraint communication setup */
+void dd_clear_local_constraint_indices(gmx_domdec_t *dd);
+
+/*! \brief Sets up communication and atom indices for all local+connected constraints */
+int dd_make_local_constraints(gmx_domdec_t *dd, int at_start,
+                              const gmx_mtop_t *mtop,
+                              const int *cginfo,
+                              gmx_constr_t constr, int nrec,
+                              t_ilist *il_local);
+
+/*! \brief Initializes the data structures for constraint communication */
+void init_domdec_constraints(gmx_domdec_t *dd,
+                             gmx_mtop_t   *mtop);
+
+#endif
diff --git a/src/gromacs/domdec/domdec_internal.h b/src/gromacs/domdec/domdec_internal.h
new file mode 100644 (file)
index 0000000..fb56087
--- /dev/null
@@ -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 <hess@kth.se>
+ * \ingroup module_domdec
+ */
+#ifndef GMX_DOMDEC_DOMDEC_INTERNAL_H
+#define GMX_DOMDEC_DOMDEC_INTERNAL_H
+
+#include "gromacs/legacyheaders/types/commrec.h"
+
+/*! \brief Returns the DD cut-off distance for multi-body interactions */
+real dd_cutoff_multibody(const gmx_domdec_t *dd);
+
+/*! \brief Returns the DD cut-off distance for two-body interactions */
+real dd_cutoff_twobody(const gmx_domdec_t *dd);
+
+#endif
index 85c24e60946cec994ff13926593eebd25667cbb7..3710f598b645b794cb01717a77c4245bf761c480 100644 (file)
  * the research papers on the package. Check out http://www.gromacs.org.
  */
 
+/*! \internal \file
+ *
+ * \brief This file defines functions for (mostly) the domdec module
+ * to use MPI functionality
+ *
+ * \author Berk Hess <hess@kth.se>
+ * \ingroup module_domdec
+ */
+
 #include "gmxpre.h"
 
 #include "domdec_network.h"
@@ -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)
 
 
index 570080076496123ba6da625bca917f6a9bcb6d97..efadf4327f7586108293867fd52244485b606562 100644 (file)
  * the research papers on the package. Check out http://www.gromacs.org.
  */
 
-#ifndef _domdec_network_h
-#define _domdec_network_h
+/*! \libinternal \file
+ *
+ * \brief This file declares functions for (mostly) the domdec module
+ * to use MPI functionality
+ *
+ * \todo Wrap the raw dd_bcast in md.cpp into a higher-level function
+ * in the domdec module, then this file can be module-internal.
+ *
+ * \author Berk Hess <hess@kth.se>
+ * \inlibraryapi
+ * \ingroup module_domdec
+ */
+#ifndef GMX_DOMDEC_DOMDEC_NETWORK_H
+#define GMX_DOMDEC_DOMDEC_NETWORK_H
 
 #include "gromacs/legacyheaders/typedefs.h"
 
-#ifdef __cplusplus
-extern "C" {
-#endif
-
+/* \brief */
 enum {
     dddirForward, dddirBackward
 };
 
-/* Move integers in the comm. region one cell along the domain decomposition
- * in the dimension indexed by ddimind
- * forward (direction=dddirFoward) or backward (direction=dddirBackward).
+/*! \brief Move integers in the communication region one cell along
+ * the domain decomposition
+ *
+ * Moves in the dimension indexed by ddimind, either forward
+ * (direction=dddirFoward) or backward (direction=dddirBackward).
  */
 void
 dd_sendrecv_int(const gmx_domdec_t *dd,
@@ -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 */
index 087b471ac3e7357f161e97b28dbb7c8bb8635831..8d0bf79b7b2e6e1d753728f442876b0b3f42c8b4 100644 (file)
  * the research papers on the package. Check out http://www.gromacs.org.
  */
 
+/*! \internal \file
+ *
+ * \brief This file defines functions used by the domdec module
+ * in its initial setup phase.
+ *
+ * \author Berk Hess <hess@kth.se>
+ * \ingroup module_domdec
+ */
+
 #include "gmxpre.h"
 
 #include <assert.h>
 #include "gromacs/math/vec.h"
 #include "gromacs/utility/smalloc.h"
 
-/* Margin for setting up the DD grid */
+/*! \brief Margin for setting up the DD grid */
 #define DD_GRID_MARGIN_PRES_SCALE 1.05
 
+/*! \brief Factorize \p n.
+ *
+ * \param[in]    n     Value to factorize
+ * \param[out]   fac   Pointer to array of factors (to be allocated in this function)
+ * \param[out]   mfac  Pointer to array of the number of times each factor repeats in the factorization (to be allocated in this function)
+ * \return             The number of unique factors
+ */
 static int factorize(int n, int **fac, int **mfac)
 {
     int d, ndiv;
@@ -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<int>(std::pow(nnodes-npme, 1.0/3.0) + 0.5);
+    npp_root3  = static_cast<int>(std::pow(ntot - npme, 1.0/3.0) + 0.5);
     npme_root2 = static_cast<int>(std::sqrt(static_cast<double>(npme)) + 0.5);
 
     /* The check below gives a reasonable division:
-     * factor 5 allowed at 5 or more PP nodes,
-     * factor 7 allowed at 49 or more PP nodes.
+     * factor 5 allowed at 5 or more PP ranks,
+     * factor 7 allowed at 49 or more PP ranks.
      */
     if (ldiv > 3 + npp_root3)
     {
         return FALSE;
     }
 
-    /* Check if the number of PP and PME nodes have a reasonable sized
+    /* Check if the number of PP and PME ranks have a reasonable sized
      * denominator in common, such that we can use 2D PME decomposition
      * when required (which requires nx_pp == nx_pme).
      * The factor of 2 allows for a maximum ratio of 2^2=4
      * between nx_pme and ny_pme.
      */
-    if (lcd(nnodes-npme, npme)*2 < npme_root2)
+    if (lcd(ntot - npme, npme)*2 < npme_root2)
     {
         return FALSE;
     }
 
     /* Does this division gives a reasonable PME load? */
-    return fits_pme_ratio(nnodes, npme, ratio);
+    return fits_pme_ratio(ntot, npme, ratio);
 }
 
+/*! \brief Make a guess for the number of PME ranks to use. */
 static int guess_npme(FILE *fplog, gmx_mtop_t *mtop, t_inputrec *ir, matrix box,
-                      int nnodes)
+                      int nrank_tot)
 {
     float      ratio;
     int        npme;
@@ -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 (file)
index 0000000..abc52ba
--- /dev/null
@@ -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 <hess@kth.se>
+ * \ingroup module_domdec
+ */
+
+#include "gmxpre.h"
+
+#include "domdec_specatomcomm.h"
+
+#include <assert.h>
+
+#include <algorithm>
+
+#include "gromacs/domdec/domdec.h"
+#include "gromacs/domdec/domdec_network.h"
+#include "gromacs/legacyheaders/gmx_ga2la.h"
+#include "gromacs/legacyheaders/gmx_hash.h"
+#include "gromacs/legacyheaders/gmx_omp_nthreads.h"
+#include "gromacs/legacyheaders/macros.h"
+#include "gromacs/legacyheaders/types/commrec.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/pbcutil/ishift.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/gmxassert.h"
+#include "gromacs/utility/smalloc.h"
+
+void dd_move_f_specat(gmx_domdec_t *dd, gmx_domdec_specat_comm_t *spac,
+                      rvec *f, rvec *fshift)
+{
+    gmx_specatsend_t *spas;
+    rvec             *vbuf;
+    int               n, n0, n1, d, dim, dir, i;
+    ivec              vis;
+    int               is;
+    gmx_bool          bPBC, bScrew;
+
+    n = spac->at_end;
+    for (d = dd->ndim-1; d >= 0; d--)
+    {
+        dim = dd->dim[d];
+        if (dd->nc[dim] > 2)
+        {
+            /* Pulse the grid forward and backward */
+            spas = spac->spas[d];
+            n0   = spas[0].nrecv;
+            n1   = spas[1].nrecv;
+            n   -= n1 + n0;
+            vbuf = spac->vbuf;
+            /* Send and receive the coordinates */
+            dd_sendrecv2_rvec(dd, d,
+                              f+n+n1, n0, vbuf, spas[0].nsend,
+                              f+n, n1, vbuf+spas[0].nsend, spas[1].nsend);
+            for (dir = 0; dir < 2; dir++)
+            {
+                bPBC   = ((dir == 0 && dd->ci[dim] == 0) ||
+                          (dir == 1 && dd->ci[dim] == dd->nc[dim]-1));
+                bScrew = (bPBC && dd->bScrewPBC && dim == XX);
+
+                spas = &spac->spas[d][dir];
+                /* Sum the buffer into the required forces */
+                if (!bPBC || (!bScrew && fshift == NULL))
+                {
+                    for (i = 0; i < spas->nsend; i++)
+                    {
+                        rvec_inc(f[spas->a[i]], *vbuf);
+                        vbuf++;
+                    }
+                }
+                else
+                {
+                    clear_ivec(vis);
+                    vis[dim] = (dir == 0 ? 1 : -1);
+                    is       = IVEC2IS(vis);
+                    if (!bScrew)
+                    {
+                        /* Sum and add to shift forces */
+                        for (i = 0; i < spas->nsend; i++)
+                        {
+                            rvec_inc(f[spas->a[i]], *vbuf);
+                            rvec_inc(fshift[is], *vbuf);
+                            vbuf++;
+                        }
+                    }
+                    else
+                    {
+                        /* Rotate the forces */
+                        for (i = 0; i < spas->nsend; i++)
+                        {
+                            f[spas->a[i]][XX] += (*vbuf)[XX];
+                            f[spas->a[i]][YY] -= (*vbuf)[YY];
+                            f[spas->a[i]][ZZ] -= (*vbuf)[ZZ];
+                            if (fshift)
+                            {
+                                rvec_inc(fshift[is], *vbuf);
+                            }
+                            vbuf++;
+                        }
+                    }
+                }
+            }
+        }
+        else
+        {
+            /* Two cells, so we only need to communicate one way */
+            spas = &spac->spas[d][0];
+            n   -= spas->nrecv;
+            /* Send and receive the coordinates */
+            dd_sendrecv_rvec(dd, d, dddirForward,
+                             f+n, spas->nrecv, spac->vbuf, spas->nsend);
+            /* Sum the buffer into the required forces */
+            if (dd->bScrewPBC && dim == XX &&
+                (dd->ci[dim] == 0 ||
+                 dd->ci[dim] == dd->nc[dim]-1))
+            {
+                for (i = 0; i < spas->nsend; i++)
+                {
+                    /* Rotate the force */
+                    f[spas->a[i]][XX] += spac->vbuf[i][XX];
+                    f[spas->a[i]][YY] -= spac->vbuf[i][YY];
+                    f[spas->a[i]][ZZ] -= spac->vbuf[i][ZZ];
+                }
+            }
+            else
+            {
+                for (i = 0; i < spas->nsend; i++)
+                {
+                    rvec_inc(f[spas->a[i]], spac->vbuf[i]);
+                }
+            }
+        }
+    }
+}
+
+void dd_move_x_specat(gmx_domdec_t *dd, gmx_domdec_specat_comm_t *spac,
+                      matrix box,
+                      rvec *x0,
+                      rvec *x1, gmx_bool bX1IsCoord)
+{
+    gmx_specatsend_t *spas;
+    rvec             *x, *vbuf, *rbuf;
+    int               nvec, v, n, nn, ns0, ns1, nr0, nr1, nr, d, dim, dir, i;
+    gmx_bool          bPBC, bScrew = FALSE;
+    rvec              shift = {0, 0, 0};
+
+    nvec = 1;
+    if (x1 != NULL)
+    {
+        nvec++;
+    }
+
+    n = spac->at_start;
+    for (d = 0; d < dd->ndim; d++)
+    {
+        dim = dd->dim[d];
+        if (dd->nc[dim] > 2)
+        {
+            /* Pulse the grid forward and backward */
+            vbuf = spac->vbuf;
+            for (dir = 0; dir < 2; dir++)
+            {
+                if (dir == 0 && dd->ci[dim] == 0)
+                {
+                    bPBC   = TRUE;
+                    bScrew = (dd->bScrewPBC && dim == XX);
+                    copy_rvec(box[dim], shift);
+                }
+                else if (dir == 1 && dd->ci[dim] == dd->nc[dim]-1)
+                {
+                    bPBC   = TRUE;
+                    bScrew = (dd->bScrewPBC && dim == XX);
+                    for (i = 0; i < DIM; i++)
+                    {
+                        shift[i] = -box[dim][i];
+                    }
+                }
+                else
+                {
+                    bPBC   = FALSE;
+                    bScrew = FALSE;
+                }
+                spas = &spac->spas[d][dir];
+                for (v = 0; v < nvec; v++)
+                {
+                    x = (v == 0 ? x0 : x1);
+                    /* Copy the required coordinates to the send buffer */
+                    if (!bPBC || (v == 1 && !bX1IsCoord))
+                    {
+                        /* Only copy */
+                        for (i = 0; i < spas->nsend; i++)
+                        {
+                            copy_rvec(x[spas->a[i]], *vbuf);
+                            vbuf++;
+                        }
+                    }
+                    else if (!bScrew)
+                    {
+                        /* Shift coordinates */
+                        for (i = 0; i < spas->nsend; i++)
+                        {
+                            rvec_add(x[spas->a[i]], shift, *vbuf);
+                            vbuf++;
+                        }
+                    }
+                    else
+                    {
+                        /* Shift and rotate coordinates */
+                        for (i = 0; i < spas->nsend; i++)
+                        {
+                            (*vbuf)[XX] =               x[spas->a[i]][XX] + shift[XX];
+                            (*vbuf)[YY] = box[YY][YY] - x[spas->a[i]][YY] + shift[YY];
+                            (*vbuf)[ZZ] = box[ZZ][ZZ] - x[spas->a[i]][ZZ] + shift[ZZ];
+                            vbuf++;
+                        }
+                    }
+                }
+            }
+            /* Send and receive the coordinates */
+            spas = spac->spas[d];
+            ns0  = spas[0].nsend;
+            nr0  = spas[0].nrecv;
+            ns1  = spas[1].nsend;
+            nr1  = spas[1].nrecv;
+            if (nvec == 1)
+            {
+                dd_sendrecv2_rvec(dd, d,
+                                  spac->vbuf+ns0, ns1, x0+n, nr1,
+                                  spac->vbuf, ns0, x0+n+nr1, nr0);
+            }
+            else
+            {
+                /* Communicate both vectors in one buffer */
+                rbuf = spac->vbuf2;
+                dd_sendrecv2_rvec(dd, d,
+                                  spac->vbuf+2*ns0, 2*ns1, rbuf, 2*nr1,
+                                  spac->vbuf, 2*ns0, rbuf+2*nr1, 2*nr0);
+                /* Split the buffer into the two vectors */
+                nn = n;
+                for (dir = 1; dir >= 0; dir--)
+                {
+                    nr = spas[dir].nrecv;
+                    for (v = 0; v < 2; v++)
+                    {
+                        x = (v == 0 ? x0 : x1);
+                        for (i = 0; i < nr; i++)
+                        {
+                            copy_rvec(*rbuf, x[nn+i]);
+                            rbuf++;
+                        }
+                    }
+                    nn += nr;
+                }
+            }
+            n += nr0 + nr1;
+        }
+        else
+        {
+            spas = &spac->spas[d][0];
+            /* Copy the required coordinates to the send buffer */
+            vbuf = spac->vbuf;
+            for (v = 0; v < nvec; v++)
+            {
+                x = (v == 0 ? x0 : x1);
+                if (dd->bScrewPBC && dim == XX &&
+                    (dd->ci[XX] == 0 || dd->ci[XX] == dd->nc[XX]-1))
+                {
+                    /* Here we only perform the rotation, the rest of the pbc
+                     * is handled in the constraint or viste routines.
+                     */
+                    for (i = 0; i < spas->nsend; i++)
+                    {
+                        (*vbuf)[XX] =               x[spas->a[i]][XX];
+                        (*vbuf)[YY] = box[YY][YY] - x[spas->a[i]][YY];
+                        (*vbuf)[ZZ] = box[ZZ][ZZ] - x[spas->a[i]][ZZ];
+                        vbuf++;
+                    }
+                }
+                else
+                {
+                    for (i = 0; i < spas->nsend; i++)
+                    {
+                        copy_rvec(x[spas->a[i]], *vbuf);
+                        vbuf++;
+                    }
+                }
+            }
+            /* Send and receive the coordinates */
+            if (nvec == 1)
+            {
+                dd_sendrecv_rvec(dd, d, dddirBackward,
+                                 spac->vbuf, spas->nsend, x0+n, spas->nrecv);
+            }
+            else
+            {
+                /* Communicate both vectors in one buffer */
+                rbuf = spac->vbuf2;
+                dd_sendrecv_rvec(dd, d, dddirBackward,
+                                 spac->vbuf, 2*spas->nsend, rbuf, 2*spas->nrecv);
+                /* Split the buffer into the two vectors */
+                nr = spas[0].nrecv;
+                for (v = 0; v < 2; v++)
+                {
+                    x = (v == 0 ? x0 : x1);
+                    for (i = 0; i < nr; i++)
+                    {
+                        copy_rvec(*rbuf, x[n+i]);
+                        rbuf++;
+                    }
+                }
+            }
+            n += spas->nrecv;
+        }
+    }
+}
+
+int setup_specat_communication(gmx_domdec_t             *dd,
+                               ind_req_t                *ireq,
+                               gmx_domdec_specat_comm_t *spac,
+                               gmx_hash_t                ga2la_specat,
+                               int                       at_start,
+                               int                       vbuf_fac,
+                               const char               *specat_type,
+                               const char               *add_err)
+{
+    int               nsend[2], nlast, nsend_zero[2] = {0, 0}, *nsend_ptr;
+    int               d, dim, ndir, dir, nr, ns, i, nrecv_local, n0, start, indr, ind, buf[2];
+    int               nat_tot_specat, nat_tot_prev, nalloc_old;
+    gmx_bool          bPBC;
+    gmx_specatsend_t *spas;
+
+    if (debug)
+    {
+        fprintf(debug, "Begin setup_specat_communication for %s\n", specat_type);
+    }
+
+    /* nsend[0]: the number of atoms requested by this node only,
+     *           we communicate this for more efficients checks
+     * nsend[1]: the total number of requested atoms
+     */
+    nsend[0] = ireq->n;
+    nsend[1] = nsend[0];
+    nlast    = nsend[1];
+    for (d = dd->ndim-1; d >= 0; d--)
+    {
+        /* Pulse the grid forward and backward */
+        dim  = dd->dim[d];
+        bPBC = (dim < dd->npbcdim);
+        if (dd->nc[dim] == 2)
+        {
+            /* Only 2 cells, so we only need to communicate once */
+            ndir = 1;
+        }
+        else
+        {
+            ndir = 2;
+        }
+        for (dir = 0; dir < ndir; dir++)
+        {
+            if (!bPBC &&
+                dd->nc[dim] > 2 &&
+                ((dir == 0 && dd->ci[dim] == dd->nc[dim] - 1) ||
+                 (dir == 1 && dd->ci[dim] == 0)))
+            {
+                /* No pbc: the fist/last cell should not request atoms */
+                nsend_ptr = nsend_zero;
+            }
+            else
+            {
+                nsend_ptr = nsend;
+            }
+            /* Communicate the number of indices */
+            dd_sendrecv_int(dd, d, dir == 0 ? dddirForward : dddirBackward,
+                            nsend_ptr, 2, spac->nreq[d][dir], 2);
+            nr = spac->nreq[d][dir][1];
+            if (nlast+nr > ireq->nalloc)
+            {
+                ireq->nalloc = over_alloc_dd(nlast+nr);
+                srenew(ireq->ind, ireq->nalloc);
+            }
+            /* Communicate the indices */
+            dd_sendrecv_int(dd, d, dir == 0 ? dddirForward : dddirBackward,
+                            ireq->ind, nsend_ptr[1], ireq->ind+nlast, nr);
+            nlast += nr;
+        }
+        nsend[1] = nlast;
+    }
+    if (debug)
+    {
+        fprintf(debug, "Communicated the counts\n");
+    }
+
+    /* Search for the requested atoms and communicate the indices we have */
+    nat_tot_specat = at_start;
+    nrecv_local    = 0;
+    for (d = 0; d < dd->ndim; d++)
+    {
+        /* Pulse the grid forward and backward */
+        if (dd->dim[d] >= dd->npbcdim || dd->nc[dd->dim[d]] > 2)
+        {
+            ndir = 2;
+        }
+        else
+        {
+            ndir = 1;
+        }
+        nat_tot_prev = nat_tot_specat;
+        for (dir = ndir-1; dir >= 0; dir--)
+        {
+            if (nat_tot_specat > spac->bSendAtom_nalloc)
+            {
+                nalloc_old             = spac->bSendAtom_nalloc;
+                spac->bSendAtom_nalloc = over_alloc_dd(nat_tot_specat);
+                srenew(spac->bSendAtom, spac->bSendAtom_nalloc);
+                for (i = nalloc_old; i < spac->bSendAtom_nalloc; i++)
+                {
+                    spac->bSendAtom[i] = FALSE;
+                }
+            }
+            spas = &spac->spas[d][dir];
+            n0   = spac->nreq[d][dir][0];
+            nr   = spac->nreq[d][dir][1];
+            if (debug)
+            {
+                fprintf(debug, "dim=%d, dir=%d, searching for %d atoms\n",
+                        d, dir, nr);
+            }
+            start       = nlast - nr;
+            spas->nsend = 0;
+            nsend[0]    = 0;
+            for (i = 0; i < nr; i++)
+            {
+                indr = ireq->ind[start+i];
+                ind  = -1;
+                /* Check if this is a home atom and if so ind will be set */
+                if (!ga2la_get_home(dd->ga2la, indr, &ind))
+                {
+                    /* Search in the communicated atoms */
+                    ind = gmx_hash_get_minone(ga2la_specat, indr);
+                }
+                if (ind >= 0)
+                {
+                    if (i < n0 || !spac->bSendAtom[ind])
+                    {
+                        if (spas->nsend+1 > spas->a_nalloc)
+                        {
+                            spas->a_nalloc = over_alloc_large(spas->nsend+1);
+                            srenew(spas->a, spas->a_nalloc);
+                        }
+                        /* Store the local index so we know which coordinates
+                         * to send out later.
+                         */
+                        spas->a[spas->nsend] = ind;
+                        spac->bSendAtom[ind] = TRUE;
+                        if (spas->nsend+1 > spac->ibuf_nalloc)
+                        {
+                            spac->ibuf_nalloc = over_alloc_large(spas->nsend+1);
+                            srenew(spac->ibuf, spac->ibuf_nalloc);
+                        }
+                        /* Store the global index so we can send it now */
+                        spac->ibuf[spas->nsend] = indr;
+                        if (i < n0)
+                        {
+                            nsend[0]++;
+                        }
+                        spas->nsend++;
+                    }
+                }
+            }
+            nlast = start;
+            /* Clear the local flags */
+            for (i = 0; i < spas->nsend; i++)
+            {
+                spac->bSendAtom[spas->a[i]] = FALSE;
+            }
+            /* Send and receive the number of indices to communicate */
+            nsend[1] = spas->nsend;
+            dd_sendrecv_int(dd, d, dir == 0 ? dddirBackward : dddirForward,
+                            nsend, 2, buf, 2);
+            if (debug)
+            {
+                fprintf(debug, "Send to rank %d, %d (%d) indices, "
+                        "receive from rank %d, %d (%d) indices\n",
+                        dd->neighbor[d][1-dir], nsend[1], nsend[0],
+                        dd->neighbor[d][dir], buf[1], buf[0]);
+                if (gmx_debug_at)
+                {
+                    for (i = 0; i < spas->nsend; i++)
+                    {
+                        fprintf(debug, " %d", spac->ibuf[i]+1);
+                    }
+                    fprintf(debug, "\n");
+                }
+            }
+            nrecv_local += buf[0];
+            spas->nrecv  = buf[1];
+            if (nat_tot_specat + spas->nrecv > dd->gatindex_nalloc)
+            {
+                dd->gatindex_nalloc =
+                    over_alloc_dd(nat_tot_specat + spas->nrecv);
+                srenew(dd->gatindex, dd->gatindex_nalloc);
+            }
+            /* Send and receive the indices */
+            dd_sendrecv_int(dd, d, dir == 0 ? dddirBackward : dddirForward,
+                            spac->ibuf, spas->nsend,
+                            dd->gatindex+nat_tot_specat, spas->nrecv);
+            nat_tot_specat += spas->nrecv;
+        }
+
+        /* Allocate the x/f communication buffers */
+        ns = spac->spas[d][0].nsend;
+        nr = spac->spas[d][0].nrecv;
+        if (ndir == 2)
+        {
+            ns += spac->spas[d][1].nsend;
+            nr += spac->spas[d][1].nrecv;
+        }
+        if (vbuf_fac*ns > spac->vbuf_nalloc)
+        {
+            spac->vbuf_nalloc = over_alloc_dd(vbuf_fac*ns);
+            srenew(spac->vbuf, spac->vbuf_nalloc);
+        }
+        if (vbuf_fac == 2 && vbuf_fac*nr > spac->vbuf2_nalloc)
+        {
+            spac->vbuf2_nalloc = over_alloc_dd(vbuf_fac*nr);
+            srenew(spac->vbuf2, spac->vbuf2_nalloc);
+        }
+
+        /* Make a global to local index for the communication atoms */
+        for (i = nat_tot_prev; i < nat_tot_specat; i++)
+        {
+            gmx_hash_change_or_set(ga2la_specat, dd->gatindex[i], i);
+        }
+    }
+
+    /* Check that in the end we got the number of atoms we asked for */
+    if (nrecv_local != ireq->n)
+    {
+        if (debug)
+        {
+            fprintf(debug, "Requested %d, received %d (tot recv %d)\n",
+                    ireq->n, nrecv_local, nat_tot_specat-at_start);
+            if (gmx_debug_at)
+            {
+                for (i = 0; i < ireq->n; i++)
+                {
+                    ind = gmx_hash_get_minone(ga2la_specat, ireq->ind[i]);
+                    fprintf(debug, " %s%d",
+                            (ind >= 0) ? "" : "!",
+                            ireq->ind[i]+1);
+                }
+                fprintf(debug, "\n");
+            }
+        }
+        fprintf(stderr, "\nDD cell %d %d %d: Neighboring cells do not have atoms:",
+                dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
+        for (i = 0; i < ireq->n; i++)
+        {
+            if (gmx_hash_get_minone(ga2la_specat, ireq->ind[i]) < 0)
+            {
+                fprintf(stderr, " %d", ireq->ind[i]+1);
+            }
+        }
+        fprintf(stderr, "\n");
+        gmx_fatal(FARGS, "DD cell %d %d %d could only obtain %d of the %d atoms that are connected via %ss from the neighboring cells. This probably means your %s lengths are too long compared to the domain decomposition cell size. Decrease the number of domain decomposition grid cells%s%s.",
+                  dd->ci[XX], dd->ci[YY], dd->ci[ZZ],
+                  nrecv_local, ireq->n, specat_type,
+                  specat_type, add_err,
+                  dd->bGridJump ? " or use the -rcon option of mdrun" : "");
+    }
+
+    spac->at_start = at_start;
+    spac->at_end   = nat_tot_specat;
+
+    if (debug)
+    {
+        fprintf(debug, "Done setup_specat_communication\n");
+    }
+
+    return nat_tot_specat;
+}
+
+gmx_domdec_specat_comm_t *specat_comm_init(int nthread)
+{
+    gmx_domdec_specat_comm_t *spac;
+
+    snew(spac, 1);
+    spac->nthread = nthread;
+    snew(spac->ireq, spac->nthread);
+
+    return spac;
+}
diff --git a/src/gromacs/domdec/domdec_specatomcomm.h b/src/gromacs/domdec/domdec_specatomcomm.h
new file mode 100644 (file)
index 0000000..d60a98f
--- /dev/null
@@ -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 <hess@kth.se>
+ * \ingroup module_domdec
+ */
+
+#ifndef GMX_DOMDEC_DOMDEC_SPECATOMCOMM_H
+#define GMX_DOMDEC_DOMDEC_SPECATOMCOMM_H
+
+#include "gromacs/legacyheaders/types/commrec.h"
+#include "gromacs/math/vectypes.h"
+#include "gromacs/utility/basedefinitions.h"
+
+typedef struct {
+    int  nsend;
+    int *a;
+    int  a_nalloc;
+    int  nrecv;
+} gmx_specatsend_t;
+
+typedef struct {
+    int *ind;
+    int  nalloc;
+    int  n;
+} ind_req_t;
+
+/*! \internal \brief Struct with setup and buffers for special atom communication */
+typedef struct gmx_domdec_specat_comm {
+    /* The number of indices to receive during the setup */
+    int              nreq[DIM][2][2];  /**< The nr. of atoms requested, per DIM, direction and direct/indirect */
+    /* The atoms to send */
+    gmx_specatsend_t spas[DIM][2];     /**< The communication setup per DIM, direction */
+    gmx_bool        *bSendAtom;        /**< Work buffer that tells if spec.atoms should be sent */
+    int              bSendAtom_nalloc; /**< Allocation size of \p bSendAtom */
+    /* Send buffers */
+    int             *ibuf;             /**< Integer send buffer */
+    int              ibuf_nalloc;      /**< Allocation size of \p ibuf */
+    rvec            *vbuf;             /**< rvec send buffer */
+    int              vbuf_nalloc;      /**< Allocation size of \p vbuf */
+    rvec            *vbuf2;            /**< rvec send buffer */
+    int              vbuf2_nalloc;     /**< Allocation size of \p vbuf2 */
+    /* The range in the local buffer(s) for received atoms */
+    int              at_start;         /**< Start index of received atoms */
+    int              at_end;           /**< End index of received atoms */
+
+    /* The atom indices we need from the surrounding cells.
+     * We can gather the indices over nthread threads.
+     */
+    int        nthread;                /**< Number of threads used for spec.atom communication */
+    ind_req_t *ireq;                   /**< Index request buffer per thread, allocation size \p nthread */
+} gmx_domdec_specat_comm_t;
+
+/*! \brief Communicates the force for special atoms, the shift forces are reduced with \p fshift != NULL */
+void dd_move_f_specat(gmx_domdec_t *dd, gmx_domdec_specat_comm_t *spac,
+                      rvec *f, rvec *fshift);
+
+/*! \brief Communicates the coordinates for special atoms
+ *
+ * \param[in]     dd         Domain decomposition struct
+ * \param[in]     spac       Special atom communication struct
+ * \param[in]     box        Box, used for pbc
+ * \param[in,out] x0         Vector to communicate
+ * \param[in,out] x1         Vector to communicate, when != NULL
+ * \param[in]     bX1IsCoord Tells is \p x1 is a coordinate vector (needs pbc)
+ */
+void dd_move_x_specat(gmx_domdec_t *dd, gmx_domdec_specat_comm_t *spac,
+                      matrix box,
+                      rvec *x0,
+                      rvec *x1, gmx_bool bX1IsCoord);
+
+/*! \brief Sets up the communication for special atoms
+ *
+ * \param[in]     dd         Domain decomposition struct
+ * \param[in]     ireq List of requested atom indices
+ * \param[in,out] spac   Special atom communication struct
+ * \param[out]    ga2la_specat Global to local special atom index
+ * \param[in]     at_start     Index in local state where to start storing communicated atoms
+ * \param[in]     vbuf_fac     Buffer factor, 1 or 2 for communicating 1 or 2 vectors
+ * \param[in]     specat_type  Name of the special atom, used for error message
+ * \param[in]     add_err      Text to add at the end of error message when atoms can't be found
+ */
+int setup_specat_communication(gmx_domdec_t             *dd,
+                               ind_req_t                *ireq,
+                               gmx_domdec_specat_comm_t *spac,
+                               gmx_hash_t                ga2la_specat,
+                               int                       at_start,
+                               int                       vbuf_fac,
+                               const char               *specat_type,
+                               const char               *add_err);
+
+/*! \brief Initialize a special communication struct */
+gmx_domdec_specat_comm_t *specat_comm_init(int nthread);
+
+#endif
index 32f19f1eb0f6d35c82d7f8aa0bcf480280034694..39e40f66fe873f1dc6fecaac4fa1dac477a4dfc4 100644 (file)
  * the research papers on the package. Check out http://www.gromacs.org.
  */
 
+/*! \internal \file
+ *
+ * \brief This file defines functions used by the domdec module
+ * while managing the construction, use and error checking for
+ * topologies local to a DD rank.
+ *
+ * \author Berk Hess <hess@kth.se>
+ * \ingroup module_domdec
+ */
+
 #include "gmxpre.h"
 
 #include <string.h>
 #include "gromacs/utility/gmxassert.h"
 #include "gromacs/utility/smalloc.h"
 
-/* for dd_init_local_state */
+#include "domdec_constraints.h"
+#include "domdec_internal.h"
+#include "domdec_vsite.h"
+
+/*! \brief The number of integer item in the local state, used for broadcasting of the state */
 #define NITEM_DD_INIT_LOCAL_STATE 5
 
 typedef struct {
     int  *index;  /* Index for each atom into il                  */
     int  *il;     /* ftype|type|a0|...|an|ftype|...               */
-} gmx_reverse_ilist_t;
+} reverse_ilist_t;
 
 typedef struct {
     int  a_start;
     int  a_end;
     int  natoms_mol;
     int  type;
-} gmx_molblock_ind_t;
+} molblock_ind_t;
 
+/*! \brief Struct for thread local work data for local topology generation */
+typedef struct {
+    t_idef     idef;             /**< Partial local topology */
+    int      **vsite_pbc;        /**< vsite PBC structure */
+    int       *vsite_pbc_nalloc; /**< Allocation sizes for vsite_pbc */
+    int        nbonded;          /**< The number of bondeds in this struct */
+    t_blocka   excl;             /**< List of exclusions */
+    int        excl_count;       /**< The total exclusion count for \p excl */
+} thread_work_t;
+
+/*! \brief Struct for the reverse topology: links bonded interactions to atomsx */
 typedef struct gmx_reverse_top {
-    gmx_bool             bExclRequired; /* Do we require all exclusions to be assigned? */
-    gmx_bool             bConstr;       /* Are there constraints in this revserse top?  */
-    gmx_bool             bSettle;       /* Are there settles in this revserse top?  */
-    gmx_bool             bBCheck;       /* All bonded interactions have to be assigned? */
-    gmx_bool             bMultiCGmols;  /* Are the multi charge-group molecules?        */
-    gmx_reverse_ilist_t *ril_mt;        /* Reverse ilist for all moltypes      */
-    int                  ril_mt_tot_size;
-    int                  ilsort;        /* The sorting state of bondeds for free energy */
-    gmx_molblock_ind_t  *mbi;
-    int                  nmolblock;
+    //! @cond Doxygen_Suppress
+    gmx_bool         bExclRequired;   /**< Do we require all exclusions to be assigned? */
+    gmx_bool         bConstr;         /**< Are there constraints in this revserse top?  */
+    gmx_bool         bSettle;         /**< Are there settles in this revserse top?  */
+    gmx_bool         bBCheck;         /**< All bonded interactions have to be assigned? */
+    gmx_bool         bMultiCGmols;    /**< Are the multi charge-group molecules?        */
+    reverse_ilist_t *ril_mt;          /**< Reverse ilist for all moltypes      */
+    int              ril_mt_tot_size; /**< The size of ril_mt[?].index summed over all entries */
+    int              ilsort;          /**< The sorting state of bondeds for free energy */
+    molblock_ind_t  *mbi;             /**< molblock to global atom index for quick lookup of molblocks on atom index */
+    int              nmolblock;       /**< The number of molblocks, size of \p mbi */
 
     /* Work data structures for multi-threading */
-    int         nthread;
-    t_idef     *idef_thread;
-    int      ***vsite_pbc;
-    int       **vsite_pbc_nalloc;
-    int        *nbonded_thread;
-    t_blocka   *excl_thread;
-    int        *excl_count_thread;
+    int            nthread;           /**< The number of threads to be used */
+    thread_work_t *th_work;           /**< Thread work array for local topology generation */
 
     /* Pointers only used for an error message */
-    gmx_mtop_t     *err_top_global;
-    gmx_localtop_t *err_top_local;
+    gmx_mtop_t     *err_top_global; /**< Pointer to the global top, only used for error reporting */
+    gmx_localtop_t *err_top_local;  /**< Pointer to the local top, only used for error reporting */
+    //! @endcond
 } gmx_reverse_top_t;
 
+/*! \brief Returns the number of atom entries for il in gmx_reverse_top_t */
 static int nral_rt(int ftype)
 {
-    /* Returns the number of atom entries for il in gmx_reverse_top_t */
     int nral;
 
     nral = NRAL(ftype);
@@ -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 (file)
index 0000000..e9d063f
--- /dev/null
@@ -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 <hess@kth.se>
+ * \ingroup module_domdec
+ */
+
+#include "gmxpre.h"
+
+#include "domdec_vsite.h"
+
+#include <assert.h>
+
+#include <algorithm>
+
+#include "gromacs/domdec/domdec.h"
+#include "gromacs/legacyheaders/gmx_ga2la.h"
+#include "gromacs/legacyheaders/gmx_hash.h"
+#include "gromacs/legacyheaders/gmx_omp_nthreads.h"
+#include "gromacs/legacyheaders/macros.h"
+#include "gromacs/legacyheaders/types/commrec.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/pbcutil/ishift.h"
+#include "gromacs/topology/mtop_util.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/gmxassert.h"
+#include "gromacs/utility/smalloc.h"
+
+#include "domdec_specatomcomm.h"
+
+void dd_move_f_vsites(gmx_domdec_t *dd, rvec *f, rvec *fshift)
+{
+    if (dd->vsite_comm)
+    {
+        dd_move_f_specat(dd, dd->vsite_comm, f, fshift);
+    }
+}
+
+void dd_clear_f_vsites(gmx_domdec_t *dd, rvec *f)
+{
+    int i;
+
+    if (dd->vsite_comm)
+    {
+        for (i = dd->vsite_comm->at_start; i < dd->vsite_comm->at_end; i++)
+        {
+            clear_rvec(f[i]);
+        }
+    }
+}
+
+void dd_move_x_vsites(gmx_domdec_t *dd, matrix box, rvec *x)
+{
+    if (dd->vsite_comm)
+    {
+        dd_move_x_specat(dd, dd->vsite_comm, box, x, NULL, FALSE);
+    }
+}
+
+int dd_make_local_vsites(gmx_domdec_t *dd, int at_start, t_ilist *lil)
+{
+    gmx_domdec_specat_comm_t *spac;
+    ind_req_t                *ireq;
+    gmx_hash_t                ga2la_specat;
+    int  ftype, nral, i, j, a;
+    t_ilist                  *lilf;
+    t_iatom                  *iatoms;
+    int  at_end;
+
+    spac         = dd->vsite_comm;
+    ireq         = &spac->ireq[0];
+    ga2la_specat = dd->ga2la_vsite;
+
+    ireq->n = 0;
+    /* Loop over all the home vsites */
+    for (ftype = 0; ftype < F_NRE; ftype++)
+    {
+        if (interaction_function[ftype].flags & IF_VSITE)
+        {
+            nral = NRAL(ftype);
+            lilf = &lil[ftype];
+            for (i = 0; i < lilf->nr; i += 1+nral)
+            {
+                iatoms = lilf->iatoms + i;
+                /* Check if we have the other atoms */
+                for (j = 1; j < 1+nral; j++)
+                {
+                    if (iatoms[j] < 0)
+                    {
+                        /* This is not a home atom,
+                         * we need to ask our neighbors.
+                         */
+                        a = -iatoms[j] - 1;
+                        /* Check to not ask for the same atom more than once */
+                        if (gmx_hash_get_minone(dd->ga2la_vsite, a) == -1)
+                        {
+                            /* Add this non-home atom to the list */
+                            if (ireq->n+1 > ireq->nalloc)
+                            {
+                                ireq->nalloc = over_alloc_large(ireq->n+1);
+                                srenew(ireq->ind, ireq->nalloc);
+                            }
+                            ireq->ind[ireq->n++] = a;
+                            /* Temporarily mark with -2,
+                             * we get the index later.
+                             */
+                            gmx_hash_set(ga2la_specat, a, -2);
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    at_end = setup_specat_communication(dd, ireq, dd->vsite_comm, ga2la_specat,
+                                        at_start, 1, "vsite", "");
+
+    /* Fill in the missing indices */
+    for (ftype = 0; ftype < F_NRE; ftype++)
+    {
+        if (interaction_function[ftype].flags & IF_VSITE)
+        {
+            nral = NRAL(ftype);
+            lilf = &lil[ftype];
+            for (i = 0; i < lilf->nr; i += 1+nral)
+            {
+                iatoms = lilf->iatoms + i;
+                for (j = 1; j < 1+nral; j++)
+                {
+                    if (iatoms[j] < 0)
+                    {
+                        iatoms[j] = gmx_hash_get_minone(ga2la_specat, -iatoms[j]-1);
+                    }
+                }
+            }
+        }
+    }
+
+    return at_end;
+}
+
+void init_domdec_vsites(gmx_domdec_t *dd, int n_intercg_vsite)
+{
+    if (debug)
+    {
+        fprintf(debug, "Begin init_domdec_vsites\n");
+    }
+
+    /* Use a hash table for the global to local index.
+     * The number of keys is a rough estimate, it will be optimized later.
+     */
+    dd->ga2la_vsite = gmx_hash_init(std::min(n_intercg_vsite/20,
+                                             n_intercg_vsite/(2*dd->nnodes)));
+
+    dd->vsite_comm = specat_comm_init(1);
+}
diff --git a/src/gromacs/domdec/domdec_vsite.h b/src/gromacs/domdec/domdec_vsite.h
new file mode 100644 (file)
index 0000000..7c42b39
--- /dev/null
@@ -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 <hess@kth.se>
+ * \ingroup module_domdec
+ */
+
+#ifndef GMX_DOMDEC_DOMDEC_VSITE_H
+#define GMX_DOMDEC_DOMDEC_VSITE_H
+
+#include "gromacs/legacyheaders/types/commrec_fwd.h"
+#include "gromacs/topology/idef.h"
+#include "gromacs/topology/topology.h"
+
+/*! \brief Clears the local indices for the virtual site communication setup */
+void dd_clear_local_vsite_indices(gmx_domdec_t *dd);
+
+/*! \brief Sets up communication and atom indices for all local vsites */
+int dd_make_local_vsites(gmx_domdec_t *dd, int at_start, t_ilist *lil);
+
+/*! \brief Initializes the data structures for virtual site communication */
+void init_domdec_vsites(gmx_domdec_t *dd, int n_intercg_vsite);
+
+#endif