* the research papers on the package. Check out http://www.gromacs.org.
*/
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
+#include "gmxpre.h"
-#include <stdio.h>
-#include <time.h>
+#include "domdec.h"
+
+#include "config.h"
+
+#include <assert.h>
+#include <limits.h>
#include <math.h>
-#include <string.h>
+#include <stdio.h>
#include <stdlib.h>
-#include <assert.h>
+#include <string.h>
-#include "typedefs.h"
-#include "gromacs/utility/smalloc.h"
-#include "gmx_fatal.h"
-#include "gmx_fatal_collective.h"
-#include "vec.h"
-#include "domdec.h"
-#include "domdec_network.h"
-#include "nrnb.h"
-#include "pbc.h"
-#include "chargegroup.h"
-#include "constr.h"
-#include "mdatoms.h"
-#include "names.h"
-#include "force.h"
-#include "pme.h"
-#include "mdrun.h"
-#include "nsgrid.h"
-#include "shellfc.h"
-#include "mtop_util.h"
-#include "gmx_ga2la.h"
-#include "macros.h"
-#include "nbnxn_search.h"
-#include "bondf.h"
-#include "gmx_omp_nthreads.h"
-#include "gpu_utils.h"
-
-#include "gromacs/fileio/futil.h"
+#include <algorithm>
+
+#include "gromacs/domdec/domdec_network.h"
+#include "gromacs/ewald/pme.h"
#include "gromacs/fileio/gmxfio.h"
#include "gromacs/fileio/pdbio.h"
+#include "gromacs/imd/imd.h"
+#include "gromacs/legacyheaders/chargegroup.h"
+#include "gromacs/legacyheaders/constr.h"
+#include "gromacs/legacyheaders/force.h"
+#include "gromacs/legacyheaders/genborn.h"
+#include "gromacs/legacyheaders/gmx_ga2la.h"
+#include "gromacs/legacyheaders/gmx_omp_nthreads.h"
+#include "gromacs/legacyheaders/gpu_utils.h"
+#include "gromacs/legacyheaders/mdatoms.h"
+#include "gromacs/legacyheaders/mdrun.h"
+#include "gromacs/legacyheaders/names.h"
+#include "gromacs/legacyheaders/network.h"
+#include "gromacs/legacyheaders/nrnb.h"
+#include "gromacs/legacyheaders/nsgrid.h"
+#include "gromacs/legacyheaders/shellfc.h"
+#include "gromacs/legacyheaders/typedefs.h"
+#include "gromacs/legacyheaders/vsite.h"
+#include "gromacs/legacyheaders/types/commrec.h"
+#include "gromacs/legacyheaders/types/constr.h"
+#include "gromacs/legacyheaders/types/enums.h"
+#include "gromacs/legacyheaders/types/forcerec.h"
+#include "gromacs/legacyheaders/types/hw_info.h"
+#include "gromacs/legacyheaders/types/ifunc.h"
+#include "gromacs/legacyheaders/types/inputrec.h"
+#include "gromacs/legacyheaders/types/mdatom.h"
+#include "gromacs/legacyheaders/types/nrnb.h"
+#include "gromacs/legacyheaders/types/ns.h"
+#include "gromacs/legacyheaders/types/nsgrid.h"
+#include "gromacs/legacyheaders/types/shellfc.h"
+#include "gromacs/legacyheaders/types/simple.h"
+#include "gromacs/legacyheaders/types/state.h"
+#include "gromacs/listed-forces/manage-threading.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/math/vectypes.h"
+#include "gromacs/mdlib/nb_verlet.h"
+#include "gromacs/mdlib/nbnxn_search.h"
+#include "gromacs/pbcutil/ishift.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/pulling/pull.h"
+#include "gromacs/pulling/pull_rotation.h"
+#include "gromacs/swap/swapcoords.h"
#include "gromacs/timing/wallcycle.h"
+#include "gromacs/topology/block.h"
+#include "gromacs/topology/idef.h"
+#include "gromacs/topology/mtop_util.h"
+#include "gromacs/topology/topology.h"
+#include "gromacs/utility/basedefinitions.h"
+#include "gromacs/utility/basenetwork.h"
+#include "gromacs/utility/cstringutil.h"
+#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/gmxmpi.h"
-#include "gromacs/swap/swapcoords.h"
#include "gromacs/utility/qsort_threadsafe.h"
-#include "gromacs/pulling/pull.h"
-#include "gromacs/pulling/pull_rotation.h"
-#include "gromacs/imd/imd.h"
+#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)
rvec *buf, *sbuf;
ivec vis;
int is;
- gmx_bool bPBC, bScrew;
+ gmx_bool bShiftForcesNeedPbc, bScrew;
comm = dd->comm;
buf = comm->vbuf.v;
- n = 0;
nzone = comm->zones.n/2;
nat_tot = dd->nat_tot;
for (d = dd->ndim-1; d >= 0; d--)
{
- bPBC = (dd->ci[dd->dim[d]] == 0);
- bScrew = (bPBC && dd->bScrewPBC && dd->dim[d] == XX);
+ /* Only forces in domains near the PBC boundaries need to
+ consider PBC in the treatment of fshift */
+ bShiftForcesNeedPbc = (dd->ci[dd->dim[d]] == 0);
+ bScrew = (bShiftForcesNeedPbc && dd->bScrewPBC && dd->dim[d] == XX);
if (fshift == NULL && !bScrew)
{
- bPBC = FALSE;
+ bShiftForcesNeedPbc = FALSE;
}
/* Determine which shift vector we need */
clear_ivec(vis);
index = ind->index;
/* Add the received forces */
n = 0;
- if (!bPBC)
+ if (!bShiftForcesNeedPbc)
{
for (i = 0; i < ind->nsend[nzone]; i++)
{
}
else if (!bScrew)
{
+ /* fshift should always be defined if this function is
+ * called when bShiftForcesNeedPbc is true */
+ assert(NULL != fshift);
for (i = 0; i < ind->nsend[nzone]; i++)
{
at0 = cgindex[index[i]];
buf = &comm->vbuf.v[0][0];
- n = 0;
nzone = comm->zones.n/2;
nat_tot = dd->nat_tot;
for (d = dd->ndim-1; d >= 0; d--)
static void dd_move_cellx(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
rvec cell_ns_x0, rvec cell_ns_x1)
{
- int d, d1, dim, dim1, pos, buf_size, i, j, k, p, npulse, npulse_min;
+ int d, d1, dim, pos, buf_size, i, j, p, npulse, npulse_min;
gmx_ddzone_t *zp;
gmx_ddzone_t buf_s[DDZONECOMM_MAXZONE];
gmx_ddzone_t buf_r[DDZONECOMM_MAXZONE];
if (bPBC)
{
/* Take the minimum to avoid double communication */
- npulse_min = min(npulse, dd->nc[dim]-1-npulse);
+ npulse_min = std::min(npulse, dd->nc[dim]-1-npulse);
}
else
{
{
for (d1 = d; d1 < dd->ndim-1; d1++)
{
- extr_s[d1][0] = max(extr_s[d1][0], extr_r[d1][0]);
- extr_s[d1][1] = min(extr_s[d1][1], extr_r[d1][1]);
- extr_s[d1][2] = min(extr_s[d1][2], extr_r[d1][2]);
+ extr_s[d1][0] = std::max(extr_s[d1][0], extr_r[d1][0]);
+ extr_s[d1][1] = std::min(extr_s[d1][1], extr_r[d1][1]);
+ extr_s[d1][2] = std::min(extr_s[d1][2], extr_r[d1][2]);
}
}
}
{
if (bUse)
{
- buf_e[i].min0 = min(buf_e[i].min0, buf_r[i].min0);
- buf_e[i].max1 = max(buf_e[i].max1, buf_r[i].max1);
- buf_e[i].min1 = min(buf_e[i].min1, buf_r[i].min1);
+ buf_e[i].min0 = std::min(buf_e[i].min0, buf_r[i].min0);
+ buf_e[i].max1 = std::max(buf_e[i].max1, buf_r[i].max1);
+ buf_e[i].min1 = std::min(buf_e[i].min1, buf_r[i].min1);
}
if (dd->ndim == 3 && d == 0 && i == buf_size - 1)
}
if (bUse && dh[d1] >= 0)
{
- buf_e[i].mch0 = max(buf_e[i].mch0, buf_r[i].mch0-dh[d1]);
- buf_e[i].mch1 = max(buf_e[i].mch1, buf_r[i].mch1-dh[d1]);
+ buf_e[i].mch0 = std::max(buf_e[i].mch0, buf_r[i].mch0-dh[d1]);
+ buf_e[i].mch1 = std::max(buf_e[i].mch1, buf_r[i].mch1-dh[d1]);
}
}
/* Copy the received buffer to the send buffer,
for (d1 = d; d1 < dd->ndim-1; d1++)
{
- extr_s[d1][1] = min(extr_s[d1][1], buf_e[pos].min0);
- extr_s[d1][0] = max(extr_s[d1][0], buf_e[pos].max1);
- extr_s[d1][2] = min(extr_s[d1][2], buf_e[pos].min1);
+ extr_s[d1][1] = std::min(extr_s[d1][1], buf_e[pos].min0);
+ extr_s[d1][0] = std::max(extr_s[d1][0], buf_e[pos].max1);
+ extr_s[d1][2] = std::min(extr_s[d1][2], buf_e[pos].min1);
pos++;
}
{
print_ddzone(debug, 1, i, 0, &comm->zone_d1[i]);
}
- cell_ns_x0[dim] = min(cell_ns_x0[dim], comm->zone_d1[i].min0);
- cell_ns_x1[dim] = max(cell_ns_x1[dim], comm->zone_d1[i].max1);
+ cell_ns_x0[dim] = std::min(cell_ns_x0[dim], comm->zone_d1[i].min0);
+ cell_ns_x1[dim] = std::max(cell_ns_x1[dim], comm->zone_d1[i].max1);
}
}
if (dd->ndim >= 3)
{
print_ddzone(debug, 2, i, j, &comm->zone_d2[i][j]);
}
- cell_ns_x0[dim] = min(cell_ns_x0[dim], comm->zone_d2[i][j].min0);
- cell_ns_x1[dim] = max(cell_ns_x1[dim], comm->zone_d2[i][j].max1);
+ cell_ns_x0[dim] = std::min(cell_ns_x0[dim], comm->zone_d2[i][j].min0);
+ cell_ns_x1[dim] = std::max(cell_ns_x1[dim], comm->zone_d2[i][j].max1);
}
}
}
void dd_collect_vec(gmx_domdec_t *dd,
t_state *state_local, rvec *lv, rvec *v)
{
- gmx_domdec_master_t *ma;
- int n, i, c, a, nalloc = 0;
- rvec *buf = NULL;
-
dd_collect_cg(dd, state_local);
if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
{
gmx_domdec_master_t *ma;
int *scounts = NULL, *disps = NULL;
- int n, i, c, a, nalloc = 0;
+ int n, i, c, a;
rvec *buf = NULL;
if (DDMASTER(dd))
snew(grid_r, 2*dd->nnodes);
}
- dd_gather(dd, 2*sizeof(rvec), grid_s[0], DDMASTER(dd) ? grid_r[0] : NULL);
+ dd_gather(dd, 2*sizeof(rvec), grid_s, DDMASTER(dd) ? grid_r : NULL);
if (DDMASTER(dd))
{
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;
r = comm->cellsize_min[dd->dim[0]];
for (di = 1; di < dd->ndim; di++)
{
- r = min(r, comm->cellsize_min[dd->dim[di]]);
+ r = std::min(r, comm->cellsize_min[dd->dim[di]]);
}
if (comm->bBondComm)
{
- r = max(r, comm->cutoff_mbody);
+ r = std::max(r, comm->cutoff_mbody);
}
else
{
- r = min(r, comm->cutoff);
+ r = std::min(r, comm->cutoff);
}
}
}
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 max(dd->comm->cutoff, r_mb);
+ return std::max(dd->comm->cutoff, r_mb);
}
static int gmx_ddcoord2pmeindex(t_commrec *cr, int x, int y, int z)
{
gmx_domdec_t *dd;
- ivec coords, coords_pme, nc;
+ ivec coords;
int slab;
dd = cr->dd;
{
gmx_domdec_t *dd;
gmx_domdec_comm_t *comm;
- ivec coord, coord_pme;
int i;
int pmenode = -1;
if (comm->bCartesianPP_PME)
{
#ifdef GMX_MPI
+ ivec coord, coord_pme;
MPI_Cart_coords(cr->mpi_comm_mysim, sim_nodeid, DIM, coord);
if (coord[comm->cartpmedim] < dd->nc[comm->cartpmedim])
{
}
}
-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)
{
static gmx_bool receive_vir_ener(t_commrec *cr)
{
gmx_domdec_comm_t *comm;
- int pmenode, coords[DIM], rank;
+ int pmenode;
gmx_bool bReceive;
bReceive = TRUE;
{
pmenode = dd_simnode2pmenode(cr, cr->sim_nodeid);
#ifdef GMX_MPI
+ ivec coords;
MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, coords);
coords[comm->cartpmedim]++;
if (coords[comm->cartpmedim] < cr->dd->nc[comm->cartpmedim])
{
+ int rank;
MPI_Cart_rank(cr->mpi_comm_mysim, coords, &rank);
if (dd_simnode2pmenode(cr, rank) == pmenode)
{
{
int nzone, zone, zone1, cg0, cg1, cg1_p1, cg, cg_gl, a, a_gl;
int *zone2cg, *zone_ncg1, *index_gl, *gatindex;
- gmx_ga2la_t *ga2la;
- char *bLocalCG;
gmx_bool bCGs;
- bLocalCG = dd->comm->bLocalCG;
-
if (dd->nat_tot > dd->gatindex_nalloc)
{
dd->gatindex_nalloc = over_alloc_dd(dd->nat_tot);
static int check_bLocalCG(gmx_domdec_t *dd, int ncg_sys, const char *bLocalCG,
const char *where)
{
- int ncg, i, ngl, nerr;
+ int i, ngl, nerr;
nerr = 0;
if (bLocalCG == NULL)
/* The cut-off might have changed, e.g. by PME load balacning,
* from the value used to set comm->cellsize_min, so check it.
*/
- cellsize_min = max(cellsize_min, comm->cutoff/comm->cd[dim_ind].np_dlb);
+ cellsize_min = std::max(cellsize_min, comm->cutoff/comm->cd[dim_ind].np_dlb);
if (comm->bPMELoadBalDLBLimits)
{
/* Check for the cut-off limit set by the PME load balancing */
- cellsize_min = max(cellsize_min, comm->PMELoadBal_max_cutoff/comm->cd[dim_ind].np_dlb);
+ cellsize_min = std::max(cellsize_min, comm->PMELoadBal_max_cutoff/comm->cd[dim_ind].np_dlb);
}
}
{
if (comm->bPMELoadBalDLBLimits)
{
- cutoff = max(cutoff, comm->PMELoadBal_max_cutoff);
+ cutoff = std::max(cutoff, comm->PMELoadBal_max_cutoff);
}
- grid_jump_limit = max(grid_jump_limit,
- cutoff/comm->cd[dim_ind].np);
+ grid_jump_limit = std::max(grid_jump_limit,
+ cutoff/comm->cd[dim_ind].np);
}
return grid_jump_limit;
{
slab = pmeindex % ddpme->nslab;
}
- ddpme->pp_min[slab] = min(ddpme->pp_min[slab], xyz[dimind]);
- ddpme->pp_max[slab] = max(ddpme->pp_max[slab], xyz[dimind]);
+ ddpme->pp_min[slab] = std::min(ddpme->pp_min[slab], xyz[dimind]);
+ ddpme->pp_max[slab] = std::max(ddpme->pp_max[slab], xyz[dimind]);
}
}
{
npulse[d]++;
}
- cellsize_min[d] = min(cellsize_min[d], cellsize);
+ cellsize_min[d] = std::min(cellsize_min[d], cellsize);
}
if (setmode == setcellsizeslbLOCAL)
{
gmx_bool bUniform, gmx_int64_t step)
{
gmx_domdec_comm_t *comm;
- int ncd, d1, i, j, pos;
+ int ncd, d1, i, pos;
real *cell_size;
real load_aver, load_i, imbalance, change, change_max, sc;
real cellsize_limit_f, dist_min_f, dist_min_f_hard, space;
imbalance = (load_i - load_aver)/(load_aver > 0 ? load_aver : 1);
/* Determine the change of the cell size using underrelaxation */
change = -relax*imbalance;
- change_max = max(change_max, max(change, -change));
+ change_max = std::max(change_max, std::max(change, -change));
}
/* Limit the amount of scaling.
* We need to use the same rescaling for all cells in one row,
gmx_ddbox_t *ddbox)
{
gmx_domdec_comm_t *comm;
- int d1, dim1, pos;
+ int d1, pos;
comm = dd->comm;
{
gmx_domdec_master_t *ma;
int **tmp_ind = NULL, *tmp_nalloc = NULL;
- int i, icg, j, k, k0, k1, d, npbcdim;
+ int i, icg, j, k, k0, k1, d;
matrix tcm;
- rvec box_size, cg_cm;
+ rvec cg_cm;
ivec ind;
real nrcg, inv_ncg, pos_d;
atom_id *cgindex;
- gmx_bool bUnbounded, bScrew;
+ gmx_bool bScrew;
ma = dd->ma;
}
dd_scatterv(dd,
- DDMASTER(dd) ? ma->ibuf : NULL,
- DDMASTER(dd) ? ma->ibuf+dd->nnodes : NULL,
- DDMASTER(dd) ? ma->cg : NULL,
+ bMaster ? ma->ibuf : NULL,
+ bMaster ? ma->ibuf+dd->nnodes : NULL,
+ bMaster ? ma->cg : NULL,
dd->ncg_home*sizeof(int), dd->index_gl);
/* Determine the home charge group sizes */
int *move)
{
int npbcdim;
- int c, i, cg, k, k0, k1, d, dim, dim2, dir, d2, d3, d4, cell_d;
- int mc, cdd, nrcg, ncg_recv, nat_recv, nvs, nvr, nvec, vec;
+ int cg, k, k0, k1, d, dim, d2;
+ int mc, nrcg;
int flag;
gmx_bool bScrew;
ivec dev;
int *move;
int npbcdim;
int ncg[DIM*2], nat[DIM*2];
- int c, i, cg, k, k0, k1, d, dim, dim2, dir, d2, d3, d4, cell_d;
- int mc, cdd, nrcg, ncg_recv, nat_recv, nvs, nvr, nvec, vec;
+ int c, i, cg, k, d, dim, dim2, dir, d2, d3;
+ int mc, cdd, nrcg, ncg_recv, nvs, nvr, nvec, vec;
int sbuf[2], rbuf[2];
int home_pos_cg, home_pos_at, buf_pos;
int flag;
gmx_bool bV = FALSE, bSDX = FALSE, bCGP = FALSE;
- gmx_bool bScrew;
- ivec dev;
- real inv_ncg, pos_d;
+ real pos_d;
matrix tcm;
- rvec *cg_cm = NULL, cell_x0, cell_x1, limitd, limit0, limit1, cm_new;
+ rvec *cg_cm = NULL, cell_x0, cell_x1, limitd, limit0, limit1;
atom_id *cgindex;
cginfo_mb_t *cginfo_mb;
gmx_domdec_comm_t *comm;
{
dim = dd->dim[d];
ncg_recv = 0;
- nat_recv = 0;
nvr = 0;
for (dir = 0; dir < (dd->nc[dim] == 2 ? 1 : 2); dir++)
{
comm->cgcm_state[cdd], nvs,
comm->vbuf.v+nvr, i);
ncg_recv += rbuf[0];
- nat_recv += rbuf[1];
nvr += i;
}
gmx_domdec_comm_t *comm;
gmx_domdec_load_t *load;
gmx_domdec_root_t *root = NULL;
- int d, dim, cid, i, pos;
+ int d, dim, i, pos;
float cell_frac = 0, sbuf[DD_NLOAD_MAX];
gmx_bool bSepPME;
for (i = 0; i < dd->nc[dim]; i++)
{
load->sum += load->load[pos++];
- load->max = max(load->max, load->load[pos]);
+ load->max = std::max(load->max, load->load[pos]);
pos++;
if (dd->bGridJump)
{
/* This direction could not be load balanced properly,
* therefore we need to use the maximum iso the average load.
*/
- load->sum_m = max(load->sum_m, load->load[pos]);
+ load->sum_m = std::max(load->sum_m, load->load[pos]);
}
else
{
load->sum_m += load->load[pos];
}
pos++;
- load->cvol_min = min(load->cvol_min, load->load[pos]);
+ load->cvol_min = std::min(load->cvol_min, load->load[pos]);
pos++;
if (d < dd->ndim-1)
{
}
if (bSepPME)
{
- load->mdf = max(load->mdf, load->load[pos]);
+ load->mdf = std::max(load->mdf, load->load[pos]);
pos++;
- load->pme = max(load->pme, load->load[pos]);
+ load->pme = std::max(load->pme, load->load[pos]);
pos++;
}
}
void setup_dd_grid(FILE *fplog, gmx_domdec_t *dd)
{
- gmx_bool bZYX;
int d, dim, i, j, m;
ivec tmp, s;
int nzone, nzonep;
static void make_pp_communicator(FILE *fplog, t_commrec *cr, int gmx_unused reorder)
{
gmx_domdec_t *dd;
+ dd = cr->dd;
+
+#ifdef GMX_MPI
gmx_domdec_comm_t *comm;
- int i, rank, *buf;
+ int rank, *buf;
ivec periods;
-#ifdef GMX_MPI
MPI_Comm comm_cart;
-#endif
- dd = cr->dd;
comm = dd->comm;
-#ifdef GMX_MPI
if (comm->bCartesianPP)
{
/* Set up cartesian communication for the particle-particle part */
dd->nc[XX], dd->nc[YY], dd->nc[ZZ]);
}
- for (i = 0; i < DIM; i++)
+ for (int i = 0; i < DIM; i++)
{
periods[i] = TRUE;
}
/* Determine the master coordinates and rank.
* The DD master should be the same node as the master of this sim.
*/
- for (i = 0; i < dd->nnodes; i++)
+ for (int i = 0; i < dd->nnodes; i++)
{
if (comm->ddindex2simnodeid[i] == 0)
{
}
}
-static void receive_ddindex2simnodeid(t_commrec *cr)
+static void receive_ddindex2simnodeid(t_commrec gmx_unused *cr)
{
+#ifdef GMX_MPI
gmx_domdec_t *dd;
-
gmx_domdec_comm_t *comm;
- int *buf;
dd = cr->dd;
comm = dd->comm;
-#ifdef GMX_MPI
if (!comm->bCartesianPP_PME && comm->bCartesianPP)
{
+ int *buf;
snew(comm->ddindex2simnodeid, dd->nnodes);
snew(buf, dd->nnodes);
if (cr->duty & DUTY_PP)
{
buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
}
-#ifdef GMX_MPI
/* Communicate the ddindex to simulation nodeid index */
MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
cr->mpi_comm_mysim);
-#endif
sfree(buf);
}
#endif
{
gmx_domdec_t *dd;
gmx_domdec_comm_t *comm;
- int i, rank;
+ int i;
gmx_bool bDiv[DIM];
- ivec periods;
#ifdef GMX_MPI
MPI_Comm comm_cart;
#endif
#ifdef GMX_MPI
if (comm->bCartesianPP_PME)
{
+ int rank;
+ ivec periods;
+
if (fplog)
{
fprintf(fplog, "Will use a Cartesian communicator for PP <-> PME: %d x %d x %d\n", comm->ntot[XX], comm->ntot[YY], comm->ntot[ZZ]);
}
MPI_Cart_create(cr->mpi_comm_mysim, DIM, comm->ntot, periods, reorder,
&comm_cart);
-
MPI_Comm_rank(comm_cart, &rank);
if (MASTERNODE(cr) && rank != 0)
{
for (i = 0; i < nc; i++)
{
dbl = 0;
- sscanf(size_string, "%lf%n", &dbl, &n);
+ sscanf(size_string, "%20lf%n", &dbl, &n);
if (dbl == 0)
{
gmx_fatal(FARGS, "Incorrect or not enough DD cell size entries for direction %s: '%s'", dir, size_string);
val = getenv(env_var);
if (val)
{
- if (sscanf(val, "%d", &nst) <= 0)
+ if (sscanf(val, "%20d", &nst) <= 0)
{
nst = 1;
}
{
d = dd->dim[di];
/* Check using the initial average cell size */
- r = min(r, ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
+ r = std::min(r, ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
}
return r;
const char *dlb_opt, gmx_bool bRecordLoad,
unsigned long Flags, t_inputrec *ir)
{
- gmx_domdec_t *dd;
int eDLB = -1;
char buf[STRLEN];
gmx_domdec_t *dd;
gmx_domdec_comm_t *comm;
int recload;
- int d, i, j;
real r_2b, r_mb, r_bonded = -1, r_bonded_limit = -1, limit, acs;
gmx_bool bC;
char buf[STRLEN];
+ const real tenPercentMargin = 1.1;
if (fplog)
{
* within nstlist steps. Since boundaries are allowed to displace by half
* a cell size, DD cells should be at least the size of the list buffer.
*/
- comm->cellsize_limit = max(comm->cellsize_limit,
- ir->rlistlong - max(ir->rvdw, ir->rcoulomb));
+ comm->cellsize_limit = std::max(comm->cellsize_limit,
+ ir->rlistlong - std::max(ir->rvdw, ir->rcoulomb));
if (comm->bInterCGBondeds)
{
}
else
{
- comm->cutoff = max(comm->cutoff, comm->cutoff_mbody);
+ comm->cutoff = std::max(comm->cutoff, comm->cutoff_mbody);
}
r_bonded_limit = comm->cutoff_mbody;
}
*/
if (Flags & MD_DDBONDCOMM)
{
- if (max(r_2b, r_mb) > comm->cutoff)
+ if (std::max(r_2b, r_mb) > comm->cutoff)
{
- r_bonded = max(r_2b, r_mb);
- r_bonded_limit = 1.1*r_bonded;
+ r_bonded = std::max(r_2b, r_mb);
+ r_bonded_limit = tenPercentMargin*r_bonded;
comm->bBondComm = TRUE;
}
else
{
r_bonded = r_mb;
- r_bonded_limit = min(1.1*r_bonded, comm->cutoff);
+ r_bonded_limit = std::min(tenPercentMargin*r_bonded, comm->cutoff);
}
/* We determine cutoff_mbody later */
}
/* No special bonded communication,
* simply increase the DD cut-off.
*/
- r_bonded_limit = 1.1*max(r_2b, r_mb);
+ r_bonded_limit = tenPercentMargin*std::max(r_2b, r_mb);
comm->cutoff_mbody = r_bonded_limit;
- comm->cutoff = max(comm->cutoff, comm->cutoff_mbody);
+ comm->cutoff = std::max(comm->cutoff, comm->cutoff_mbody);
}
}
if (fplog)
"Minimum cell size due to bonded interactions: %.3f nm\n",
r_bonded_limit);
}
- comm->cellsize_limit = max(comm->cellsize_limit, r_bonded_limit);
+ comm->cellsize_limit = std::max(comm->cellsize_limit, r_bonded_limit);
}
if (dd->bInterCGcons && rconstr <= 0)
"User supplied maximum distance required for P-LINCS: %.3f nm\n",
rconstr);
}
- comm->cellsize_limit = max(comm->cellsize_limit, rconstr);
+ comm->cellsize_limit = std::max(comm->cellsize_limit, rconstr);
comm->cgs_gl = gmx_mtop_global_cgs(mtop);
if (comm->eDLB != edlbNO)
{
/* Check if this does not limit the scaling */
- comm->cutoff_mbody = min(comm->cutoff_mbody, dlb_scale*acs);
+ comm->cutoff_mbody = std::min(comm->cutoff_mbody, dlb_scale*acs);
}
if (!comm->bBondComm)
{
/* Without bBondComm do not go beyond the n.b. cut-off */
- comm->cutoff_mbody = min(comm->cutoff_mbody, comm->cutoff);
+ comm->cutoff_mbody = std::min(comm->cutoff_mbody, comm->cutoff);
if (comm->cellsize_limit >= comm->cutoff)
{
/* We don't loose a lot of efficieny
}
}
/* Check if we did not end up below our original limit */
- comm->cutoff_mbody = max(comm->cutoff_mbody, r_bonded_limit);
+ comm->cutoff_mbody = std::max(comm->cutoff_mbody, r_bonded_limit);
if (comm->cutoff_mbody > comm->cellsize_limit)
{
cellsize_min = comm->cellsize_min[dd->dim[0]];
for (d = 1; d < dd->ndim; d++)
{
- cellsize_min = min(cellsize_min, comm->cellsize_min[dd->dim[d]]);
+ cellsize_min = std::min(cellsize_min, comm->cellsize_min[dd->dim[d]]);
}
if (cellsize_min < comm->cellsize_limit*1.05)
t_inputrec *ir, gmx_bool bBCheck, cginfo_mb_t *cginfo_mb)
{
gmx_domdec_comm_t *comm;
- gmx_bool bBondComm;
- int d;
dd_make_reverse_top(fplog, dd, mtop, vsite, ir, bBCheck);
limit = dd->comm->cellsize_min[XX];
for (d = 1; d < DIM; d++)
{
- limit = min(limit, dd->comm->cellsize_min[d]);
+ limit = std::min(limit, dd->comm->cellsize_min[d]);
}
}
{
fprintf(fplog, "%40s %-7s %6.3f nm\n",
"two-body bonded interactions", "(-rdd)",
- max(comm->cutoff, comm->cutoff_mbody));
+ std::max(comm->cutoff, comm->cutoff_mbody));
fprintf(fplog, "%40s %-7s %6.3f nm\n",
"multi-body bonded interactions", "(-rdd)",
- (comm->bBondComm || dd->bGridJump) ? comm->cutoff_mbody : min(comm->cutoff, limit));
+ (comm->bBondComm || dd->bGridJump) ? comm->cutoff_mbody : std::min(comm->cutoff, limit));
}
if (dd->vsite_comm)
{
/* Determine the maximum number of comm. pulses in one dimension */
- comm->cellsize_limit = max(comm->cellsize_limit, comm->cutoff_mbody);
+ comm->cellsize_limit = std::max(comm->cellsize_limit, comm->cutoff_mbody);
/* Determine the maximum required number of grid pulses */
if (comm->cellsize_limit >= comm->cutoff)
else
{
/* There is no cell size limit */
- npulse = max(dd->nc[XX]-1, max(dd->nc[YY]-1, dd->nc[ZZ]-1));
+ npulse = std::max(dd->nc[XX]-1, std::max(dd->nc[YY]-1, dd->nc[ZZ]-1));
}
if (!bNoCutOff && npulse > 1)
dim = dd->dim[d];
npulse_d = (int)(1 + dd->nc[dim]*comm->cutoff
/(ddbox->box_size[dim]*ddbox->skew_fac[dim]*dlb_scale));
- npulse_d_max = max(npulse_d_max, npulse_d);
+ npulse_d_max = std::max(npulse_d_max, npulse_d);
}
- npulse = min(npulse, npulse_d_max);
+ npulse = std::min(npulse, npulse_d_max);
}
/* This env var can override npulse */
comm->bVacDLBNoLimit = (ir->ePBC == epbcNONE);
for (d = 0; d < dd->ndim; d++)
{
- comm->cd[d].np_dlb = min(npulse, dd->nc[dd->dim[d]]-1);
+ comm->cd[d].np_dlb = std::min(npulse, dd->nc[dd->dim[d]]-1);
comm->cd[d].np_nalloc = comm->cd[d].np_dlb;
snew(comm->cd[d].ind, comm->cd[d].np_nalloc);
- comm->maxpulse = max(comm->maxpulse, comm->cd[d].np_dlb);
+ comm->maxpulse = std::max(comm->maxpulse, comm->cd[d].np_dlb);
if (comm->cd[d].np_dlb < dd->nc[dd->dim[d]]-1)
{
comm->bVacDLBNoLimit = FALSE;
/* cellsize_limit is set for LINCS in init_domain_decomposition */
if (!comm->bVacDLBNoLimit)
{
- comm->cellsize_limit = max(comm->cellsize_limit,
- comm->cutoff/comm->maxpulse);
+ comm->cellsize_limit = std::max(comm->cellsize_limit,
+ comm->cutoff/comm->maxpulse);
}
- comm->cellsize_limit = max(comm->cellsize_limit, comm->cutoff_mbody);
+ comm->cellsize_limit = std::max(comm->cellsize_limit, comm->cutoff_mbody);
/* Set the minimum cell size for each DD dimension */
for (d = 0; d < dd->ndim; d++)
{
}
if (comm->cutoff_mbody <= 0)
{
- comm->cutoff_mbody = min(comm->cutoff, comm->cellsize_limit);
+ comm->cutoff_mbody = std::min(comm->cutoff, comm->cellsize_limit);
}
if (comm->bDynLoadBal)
{
}
natoms_tot = comm->cgs_gl.index[comm->cgs_gl.nr];
- dd->ga2la = ga2la_init(natoms_tot, vol_frac*natoms_tot);
+ dd->ga2la = ga2la_init(natoms_tot, static_cast<int>(vol_frac*natoms_tot));
}
static gmx_bool test_dd_cutoff(t_commrec *cr,
c->c[1][1] = comm->cell_x0[dim1];
if (dd->bGridJump)
{
- c->c[1][1] = max(comm->cell_x0[dim1], comm->zone_d1[1].mch0);
+ c->c[1][1] = std::max(comm->cell_x0[dim1], comm->zone_d1[1].mch0);
if (bDistMB)
{
/* For the multi-body distance we need the maximum */
- c->bc[1] = max(comm->cell_x0[dim1], comm->zone_d1[1].p1_0);
+ c->bc[1] = std::max(comm->cell_x0[dim1], comm->zone_d1[1].p1_0);
}
}
/* Set the upper-right corner for rounding */
if (j >= 4)
{
c->c[2][j-4] =
- max(c->c[2][j-4],
- comm->zone_d2[zones->shift[i][dim0]][zones->shift[i][dim1]].mch0);
+ std::max(c->c[2][j-4],
+ comm->zone_d2[zones->shift[i][dim0]][zones->shift[i][dim1]].mch0);
}
}
}
{
for (j = 0; j < 2; j++)
{
- c->bc[2] = max(c->bc[2], comm->zone_d2[i][j].p1_0);
+ c->bc[2] = std::max(c->bc[2], comm->zone_d2[i][j].p1_0);
}
}
}
c->cr1[3] = comm->cell_x1[dim1];
if (dd->bGridJump)
{
- c->cr1[0] = max(comm->cell_x1[dim1], comm->zone_d1[1].mch1);
+ c->cr1[0] = std::max(comm->cell_x1[dim1], comm->zone_d1[1].mch1);
if (bDistMB)
{
/* For the multi-body distance we need the maximum */
- c->bcr1 = max(comm->cell_x1[dim1], comm->zone_d1[1].p1_1);
+ c->bcr1 = std::max(comm->cell_x1[dim1], comm->zone_d1[1].p1_1);
}
}
}
{
int dim_ind, dim, dim0, dim1, dim2, dimd, p, nat_tot;
int nzone, nzone_send, zone, zonei, cg0, cg1;
- int c, i, j, cg, cg_gl, nrcg;
+ int c, i, cg, cg_gl, nrcg;
int *zone_cg_range, pos_cg, *index_gl, *cgindex, *recv_i;
gmx_domdec_comm_t *comm;
gmx_domdec_zones_t *zones;
gmx_domdec_ind_t *ind;
cginfo_mb_t *cginfo_mb;
gmx_bool bBondComm, bDist2B, bDistMB, bDistBonded;
- real r_mb, r_comm2, r_scomm2, r_bcomm2, r_0, r_1, r2inc, inv_ncg;
+ real r_comm2, r_bcomm2;
dd_corners_t corners;
ivec tric_dist;
rvec *cg_cm, *normal, *v_d, *v_0 = NULL, *v_1 = NULL, *recv_vr;
for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
{
- dim = dd->dim[dim_ind];
-
/* Check if we need to use triclinic distances */
tric_dist[dim_ind] = 0;
for (i = 0; i <= dim_ind; i++)
/* The rvec buffer is also required for atom buffers
* of size nrecv in dd_move_x and dd_move_f.
*/
- i = max(cd->ind[0].nrecv[nzone+1], ind->nrecv[nzone+1]);
+ i = std::max(cd->ind[0].nrecv[nzone+1], ind->nrecv[nzone+1]);
vec_rvec_check_alloc(&comm->vbuf2, i);
}
}
gmx_domdec_comm_t *comm;
gmx_domdec_zones_t *zones;
gmx_bool bDistMB;
- int z, zi, zj0, zj1, d, dim;
+ int z, zi, d, dim;
real rcs, rcmbs;
int i, j;
- real size_j, add_tric;
real vol;
comm = dd->comm;
* With multiple pulses this will lead
* to a larger zone then strictly necessary.
*/
- zones->size[z].x1[dim] = max(zones->size[z].x1[dim],
- zones->size[zi].x1[dim]+rcmbs);
+ zones->size[z].x1[dim] = std::max(zones->size[z].x1[dim],
+ zones->size[zi].x1[dim]+rcmbs);
}
}
}
{
if (zones->shift[z][dim] > 0)
{
- zones->size[z].x1[dim] = max(zones->size[z].x1[dim],
- zones->size[zi].x1[dim]+rcs);
+ zones->size[z].x1[dim] = std::max(zones->size[z].x1[dim],
+ zones->size[zi].x1[dim]+rcs);
}
}
}
{
corner[ZZ] = zones->size[z].x1[ZZ];
}
- if (dd->ndim == 1 && box[ZZ][YY] != 0)
- {
- /* With 1D domain decomposition the cg's are not in
- * the triclinic box, but triclinic x-y and rectangular y-z.
- * Shift y back, so it will later end up at 0.
- */
- corner[YY] -= corner[ZZ]*box[ZZ][YY]/box[ZZ][ZZ];
- }
/* Apply the triclinic couplings */
assert(ddbox->npbcdim <= DIM);
for (i = YY; i < ddbox->npbcdim; i++)
{
for (j = XX; j < i; j++)
{
- corner[j] += corner[i]*box[i][j]/box[i][i];
+ /* With 1D domain decomposition the cg's are not in
+ * a triclinic box, but triclinic x-y and rectangular y/x-z.
+ * So we should ignore the coupling for the non
+ * domain-decomposed dimension of the pair x and y.
+ */
+ if (!(dd->ndim == 1 && ((dd->dim[0] == XX && j == YY) ||
+ (dd->dim[0] == YY && j == XX))))
+ {
+ corner[j] += corner[i]*box[i][j]/box[i][i];
+ }
}
}
if (c == 0)
{
for (i = 0; i < DIM; i++)
{
- corner_min[i] = min(corner_min[i], corner[i]);
- corner_max[i] = max(corner_max[i], corner[i]);
+ corner_min[i] = std::min(corner_min[i], corner[i]);
+ corner_max[i] = std::max(corner_max[i], corner[i]);
}
}
}
{
gmx_domdec_sort_t *sort;
gmx_cgsort_t *cgsort, *sort_i;
- int ncg_new, nsort2, nsort_new, i, *a, moved, *ibuf;
- int sort_last, sort_skip;
+ int ncg_new, nsort2, nsort_new, i, *a, moved;
sort = dd->comm->sort;
int ncg_home_old)
{
gmx_domdec_sort_t *sort;
- gmx_cgsort_t *cgsort, *sort_i;
+ gmx_cgsort_t *cgsort;
int *cgindex;
int ncg_new, i, *ibuf, cgsize;
rvec *vbuf;
t_block *cgs_gl;
gmx_int64_t step_pcoupl;
rvec cell_ns_x0, cell_ns_x1;
- int i, j, n, ncgindex_set, ncg_home_old = -1, ncg_moved, nat_f_novirsum;
+ int i, n, ncgindex_set, ncg_home_old = -1, ncg_moved, nat_f_novirsum;
gmx_bool bBoxChanged, bNStGlobalComm, bDoDLB, bCheckDLB, bTurnOnDLB, bLogLoad;
gmx_bool bRedist, bSortCG, bResortAll;
ivec ncells_old = {0, 0, 0}, ncells_new = {0, 0, 0}, np;
* so the extra communication cost is negligible.
*/
const int nddp_chk_dlb = 100;
-
bCheckDLB = (comm->n_load_collect == 0 ||
comm->n_load_have % nddp_chk_dlb == nddp_chk_dlb - 1);
}
* the research papers on the package. Check out http://www.gromacs.org.
*/
/* This file is completely threadsafe - keep it that way! */
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
+#include "gmxpre.h"
#include <math.h>
-#include "main.h"
-#include "constr.h"
-#include "copyrite.h"
-#include "physics.h"
-#include "vec.h"
-#include "pbc.h"
-#include "gromacs/utility/smalloc.h"
-#include "mdrun.h"
-#include "nrnb.h"
-#include "domdec.h"
-#include "mtop_util.h"
-#include "gmx_omp_nthreads.h"
+#include <stdlib.h>
+#include "gromacs/domdec/domdec.h"
#include "gromacs/fileio/gmxfio.h"
-#include "gmx_fatal.h"
+#include "gromacs/legacyheaders/constr.h"
+#include "gromacs/legacyheaders/copyrite.h"
+#include "gromacs/legacyheaders/gmx_omp_nthreads.h"
+#include "gromacs/legacyheaders/mdrun.h"
+#include "gromacs/legacyheaders/nrnb.h"
+#include "gromacs/legacyheaders/types/commrec.h"
+#include "gromacs/math/units.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/mdlib/bitmask.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/topology/block.h"
+#include "gromacs/topology/mtop_util.h"
+#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/gmxomp.h"
+#include "gromacs/utility/smalloc.h"
typedef struct {
int b0; /* first constraint for this thread */
int *ind_r; /* constraint index for updating atom data */
int ind_nalloc; /* allocation size of ind and ind_r */
tensor vir_r_m_dr; /* temporary variable for virial calculation */
+ real dhdlambda; /* temporary variable for lambda derivative */
} lincs_thread_t;
typedef struct gmx_lincsdata {
real *bllen; /* the reference bond length */
int nth; /* The number of threads doing LINCS */
lincs_thread_t *th; /* LINCS thread division */
- unsigned *atf; /* atom flags for thread parallelization */
+ gmx_bitmask_t *atf; /* atom flags for thread parallelization */
int atf_nalloc; /* allocation size of atf */
/* arrays for temporary storage in the LINCS algorithm */
rvec *tmpv;
static void do_lincsp(rvec *x, rvec *f, rvec *fp, t_pbc *pbc,
struct gmx_lincsdata *lincsd, int th,
real *invmass,
- int econq, real *dvdlambda,
+ int econq, gmx_bool bCalcDHDL,
gmx_bool bCalcVir, tensor rmdf)
{
int b0, b1, b, i, j, k, n;
- real tmp0, tmp1, tmp2, im1, im2, mvb, rlen, len, wfac, lam;
+ real tmp0, tmp1, tmp2, mvb;
rvec dx;
int *bla, *blnr, *blbnb;
rvec *r;
lincs_update_atoms(lincsd, th, 1.0, sol, r,
(econq != econqForce) ? invmass : NULL, fp);
- if (dvdlambda != NULL)
+ if (bCalcDHDL)
{
- #pragma omp barrier
+ real dhdlambda;
+
+ dhdlambda = 0;
for (b = b0; b < b1; b++)
{
- *dvdlambda -= sol[b]*lincsd->ddist[b];
+ dhdlambda -= sol[b]*lincsd->ddist[b];
}
- /* 10 ncons flops */
+
+ lincsd->th[th].dhdlambda = dhdlambda;
}
if (bCalcVir)
struct gmx_lincsdata *lincsd, int th,
real *invmass,
t_commrec *cr,
- gmx_bool bCalcLambda,
+ gmx_bool bCalcDHDL,
real wangle, int *warn,
real invdt, rvec *v,
gmx_bool bCalcVir, tensor vir_r_m_dr)
{
int b0, b1, b, i, j, k, n, iter;
- real tmp0, tmp1, tmp2, im1, im2, mvb, rlen, len, len2, dlen2, wfac;
+ real tmp0, tmp1, tmp2, mvb, rlen, len, len2, dlen2, wfac;
rvec dx;
int *bla, *blnr, *blbnb;
rvec *r;
dlen2 = 2*len2 - norm2(dx);
if (dlen2 < wfac*len2 && (nlocat == NULL || nlocat[b]))
{
+ /* not race free - see detailed comment in caller */
*warn = b;
}
if (dlen2 > 0)
/* 16 ncons flops */
}
- if (nlocat != NULL && bCalcLambda)
+ if (nlocat != NULL && (bCalcDHDL || bCalcVir))
{
- /* In lincs_update_atoms thread might cross-read mlambda */
+ /* In lincs_update_atoms threads might cross-read mlambda */
#pragma omp barrier
/* Only account for local atoms */
}
}
+ if (bCalcDHDL)
+ {
+ real dhdl;
+
+ dhdl = 0;
+ for (b = b0; b < b1; b++)
+ {
+ /* Note that this this is dhdl*dt^2, the dt^2 factor is corrected
+ * later after the contributions are reduced over the threads.
+ */
+ dhdl -= lincsd->mlambda[b]*lincsd->ddist[b];
+ }
+ lincsd->th[th].dhdlambda = dhdl;
+ }
+
if (bCalcVir)
{
/* Constraint virial */
li->triangle[li->ntriangle] = i;
li->tri_bits[li->ntriangle] = 0;
li->ntriangle++;
- if (li->blnr[i+1] - li->blnr[i] > sizeof(li->tri_bits[0])*8 - 1)
+ if (li->blnr[i+1] - li->blnr[i] > static_cast<int>(sizeof(li->tri_bits[0])*8 - 1))
{
gmx_fatal(FARGS, "A constraint is connected to %d constraints, this is more than the %d allowed for constraints participating in triangles",
li->blnr[i+1] - li->blnr[i],
{
lincs_thread_t *li_m;
int th;
- unsigned *atf;
+ gmx_bitmask_t *atf;
int a;
if (natoms > li->atf_nalloc)
/* Clear the atom flags */
for (a = 0; a < natoms; a++)
{
- atf[a] = 0;
+ bitmask_clear(&atf[a]);
+ }
+
+ if (li->nth > BITMASK_SIZE)
+ {
+ gmx_fatal(FARGS, "More than %d threads is not supported for LINCS.", BITMASK_SIZE);
}
for (th = 0; th < li->nth; th++)
li_th->b0 = (li->nc* th )/li->nth;
li_th->b1 = (li->nc*(th+1))/li->nth;
- if (th < sizeof(*atf)*8)
+ /* For each atom set a flag for constraints from each */
+ for (b = li_th->b0; b < li_th->b1; b++)
{
- /* For each atom set a flag for constraints from each */
- for (b = li_th->b0; b < li_th->b1; b++)
- {
- atf[li->bla[b*2] ] |= (1U<<th);
- atf[li->bla[b*2+1]] |= (1U<<th);
- }
+ bitmask_set_bit(&atf[li->bla[b*2]], th);
+ bitmask_set_bit(&atf[li->bla[b*2+1]], th);
}
}
for (th = 0; th < li->nth; th++)
{
lincs_thread_t *li_th;
- unsigned mask;
+ gmx_bitmask_t mask;
int b;
li_th = &li->th[th];
srenew(li_th->ind_r, li_th->ind_nalloc);
}
- if (th < sizeof(*atf)*8)
- {
- mask = (1U<<th) - 1U;
+ bitmask_init_low_bits(&mask, th);
- li_th->nind = 0;
- li_th->nind_r = 0;
- for (b = li_th->b0; b < li_th->b1; b++)
+ li_th->nind = 0;
+ li_th->nind_r = 0;
+ for (b = li_th->b0; b < li_th->b1; b++)
+ {
+ /* We let the constraint with the lowest thread index
+ * operate on atoms with constraints from multiple threads.
+ */
+ if (bitmask_is_disjoint(atf[li->bla[b*2]], mask) &&
+ bitmask_is_disjoint(atf[li->bla[b*2+1]], mask))
{
- /* We let the constraint with the lowest thread index
- * operate on atoms with constraints from multiple threads.
- */
- if (((atf[li->bla[b*2]] & mask) == 0) &&
- ((atf[li->bla[b*2+1]] & mask) == 0))
- {
- /* Add the constraint to the local atom update index */
- li_th->ind[li_th->nind++] = b;
- }
- else
- {
- /* Add the constraint to the rest block */
- li_th->ind_r[li_th->nind_r++] = b;
- }
+ /* Add the constraint to the local atom update index */
+ li_th->ind[li_th->nind++] = b;
}
- }
- else
- {
- /* We are out of bits, assign all constraints to rest */
- for (b = li_th->b0; b < li_th->b1; b++)
+ else
{
+ /* Add the constraint to the rest block */
li_th->ind_r[li_th->nind_r++] = b;
}
}
int i, k, ncc_alloc, ni, con, nconnect, concon;
int type, a1, a2;
real lenA = 0, lenB;
- gmx_bool bLocal;
li->nc = 0;
li->ncc = 0;
li->blnr[con] = nconnect;
for (i = 0; i < ni; i++)
{
- bLocal = TRUE;
type = iatom[3*i];
a1 = iatom[3*i+1];
a2 = iatom[3*i+2];
t_nrnb *nrnb,
int maxwarn, int *warncount)
{
+ gmx_bool bCalcDHDL;
char buf[STRLEN], buf2[22], buf3[STRLEN];
- int i, warn, p_imax, error;
+ int i, warn, p_imax;
real ncons_loc, p_ssd, p_max = 0;
rvec dx;
gmx_bool bOK;
bOK = TRUE;
+ /* This boolean should be set by a flag passed to this routine.
+ * We can also easily check if any constraint length is changed,
+ * if not dH/dlambda=0 and we can also set the boolean to FALSE.
+ */
+ bCalcDHDL = (ir->efep != efepNO && dvdlambda != NULL);
+
if (lincsd->nc == 0 && cr->dd == NULL)
{
if (bLog || bEner)
if (econq == econqCoord)
{
+ /* We can't use bCalcDHDL here, since NULL can be passed for dvdlambda
+ * also with efep!=fepNO.
+ */
if (ir->efep != efepNO)
{
if (md->nMassPerturbed && lincsd->matlam != md->lambda)
do_lincs(x, xprime, box, pbc, lincsd, th,
md->invmass, cr,
- bCalcVir || (ir->efep != efepNO),
+ bCalcDHDL,
ir->LincsWarnAngle, &warn,
invdt, v, bCalcVir,
th == 0 ? vir_r_m_dr : lincsd->th[th].vir_r_m_dr);
}
- if (ir->efep != efepNO)
- {
- real dt_2, dvdl = 0;
-
- /* TODO This should probably use invdt, so that sd integrator scaling works properly */
- dt_2 = 1.0/(ir->delta_t*ir->delta_t);
- for (i = 0; (i < lincsd->nc); i++)
- {
- dvdl -= lincsd->mlambda[i]*dt_2*lincsd->ddist[i];
- }
- *dvdlambda += dvdl;
- }
-
if (bLog && fplog && lincsd->nc > 0)
{
fprintf(fplog, " Rel. Constraint Deviation: RMS MAX between atoms\n");
int th = gmx_omp_get_thread_num();
do_lincsp(x, xprime, min_proj, pbc, lincsd, th,
- md->invmass, econq, ir->efep != efepNO ? dvdlambda : NULL,
+ md->invmass, econq, bCalcDHDL,
bCalcVir, th == 0 ? vir_r_m_dr : lincsd->th[th].vir_r_m_dr);
}
}
+ if (bCalcDHDL)
+ {
+ /* Reduce the dH/dlambda contributions over the threads */
+ real dhdlambda;
+ int th;
+
+ dhdlambda = 0;
+ for (th = 0; th < lincsd->nth; th++)
+ {
+ dhdlambda += lincsd->th[th].dhdlambda;
+ }
+ if (econqCoord)
+ {
+ /* dhdlambda contains dH/dlambda*dt^2, correct for this */
++ /* TODO This should probably use invdt, so that sd integrator scaling works properly */
+ dhdlambda /= ir->delta_t*ir->delta_t;
+ }
+ *dvdlambda += dhdlambda;
+ }
+
if (bCalcVir && lincsd->nth > 1)
{
for (i = 1; i < lincsd->nth; i++)
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff