From: Erik Lindahl Date: Tue, 7 Jul 2015 20:02:24 +0000 (+0200) Subject: Move old neighborsearch code to C++ X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?p=alexxy%2Fgromacs.git;a=commitdiff_plain;h=9299a15137735835796c34afa383c626c1fe2cbd Move old neighborsearch code to C++ This will eventually be removed in favor of the new verlet code, but for now it will make life simpler to have everything moved to C++. Change-Id: I3242e5f345769a42ffb99d1953bc1bef39324ff5 --- diff --git a/src/gromacs/mdlib/ns.c b/src/gromacs/mdlib/ns.cpp similarity index 97% rename from src/gromacs/mdlib/ns.c rename to src/gromacs/mdlib/ns.cpp index d6e9314fb1..e3f37fc45a 100644 --- a/src/gromacs/mdlib/ns.c +++ b/src/gromacs/mdlib/ns.cpp @@ -42,6 +42,10 @@ #include #include +#include + +#include + #include "gromacs/domdec/domdec.h" #include "gromacs/legacyheaders/force.h" #include "gromacs/legacyheaders/macros.h" @@ -93,7 +97,7 @@ static gmx_bool NOTEXCL_(t_excl e[], atom_id i, atom_id j) static int round_up_to_simd_width(int length, int simd_width) { - int offset, newlength; + int offset; offset = (simd_width > 0) ? length % simd_width : 0; @@ -132,7 +136,7 @@ static void init_nblist(FILE *log, t_nblist *nl_sr, t_nblist *nl_lr, { t_nblist *nl; int homenr; - int i, nn; + int i; for (i = 0; (i < 2); i++) { @@ -199,7 +203,6 @@ void init_neighbor_list(FILE *log, t_forcerec *fr, int homenr) */ int maxsr, maxsr_wat, maxlr, maxlr_wat; int ielec, ivdw, ielecmod, ivdwmod, type; - int solvent; int igeometry_def, igeometry_w, igeometry_ww; int i; gmx_bool bElecAndVdwSwitchDiffers; @@ -217,11 +220,11 @@ void init_neighbor_list(FILE *log, t_forcerec *fr, int homenr) * all the nlist arrays many times in a row. * The numbers seem very accurate, but they are uncritical. */ - maxsr_wat = min(fr->nWatMol, (homenr+2)/3); + maxsr_wat = std::min(fr->nWatMol, (homenr+2)/3); if (fr->bTwinRange) { maxlr = 50; - maxlr_wat = min(maxsr_wat, maxlr); + maxlr_wat = std::min(maxsr_wat, maxlr); } else { @@ -363,9 +366,7 @@ static void reset_neighbor_lists(t_forcerec *fr, gmx_bool bResetSR, gmx_bool bRe static gmx_inline void new_i_nblist(t_nblist *nlist, atom_id i_atom, int shift, int gid) { - int i, k, nri, nshift; - - nri = nlist->nri; + int nri = nlist->nri; /* Check whether we have to increase the i counter */ if ((nri == -1) || @@ -594,14 +595,14 @@ put_in_list_at(gmx_bool bHaveVdW[], t_nblist * vdwc_ww = NULL; t_nblist * coul_ww = NULL; - int i, j, jcg, igid, gid, nbl_ind, ind_ij; + int i, j, jcg, igid, gid, nbl_ind; atom_id jj, jj0, jj1, i_atom; - int i0, nicg, len; + int i0, nicg; int *cginfo; int *type, *typeB; real *charge, *chargeB; - real qi, qiB, qq, rlj; + real qi, qiB; gmx_bool bFreeEnergy, bFree, bFreeJ, bNotEx, *bPert; gmx_bool bDoVdW_i, bDoCoul_i, bDoCoul_i_sol; int iwater, jwater; @@ -1094,7 +1095,7 @@ put_in_list_adress(gmx_bool bHaveVdW[], gmx_bool bLR, gmx_bool bDoVdW, gmx_bool bDoCoul, - int solvent_opt) + int gmx_unused solvent_opt) { /* The a[] index has been removed, * to put it back in i_atom should be a[i0] and jj should be a[jj]. @@ -1105,33 +1106,26 @@ put_in_list_adress(gmx_bool bHaveVdW[], t_nblist * vdwc_adress = NULL; t_nblist * vdw_adress = NULL; t_nblist * coul_adress = NULL; - t_nblist * vdwc_ww = NULL; - t_nblist * coul_ww = NULL; int i, j, jcg, igid, gid, nbl_ind, nbl_ind_adress; atom_id jj, jj0, jj1, i_atom; - int i0, nicg, len; + int i0, nicg; int *cginfo; - int *type, *typeB; - real *charge, *chargeB; + int *type; + real *charge; real *wf; - real qi, qiB, qq, rlj; - gmx_bool bFreeEnergy, bFree, bFreeJ, bNotEx, *bPert; - gmx_bool bDoVdW_i, bDoCoul_i, bDoCoul_i_sol; + real qi; + gmx_bool bNotEx; + gmx_bool bDoVdW_i, bDoCoul_i; gmx_bool b_hybrid; - gmx_bool j_all_atom; - int iwater, jwater; t_nblist *nlist, *nlist_adress; gmx_bool bEnergyGroupCG; /* Copy some pointers */ cginfo = fr->cginfo; charge = md->chargeA; - chargeB = md->chargeB; type = md->typeA; - typeB = md->typeB; - bPert = md->bPerturbed; wf = md->wf; /* Get atom range */ @@ -1141,8 +1135,6 @@ put_in_list_adress(gmx_bool bHaveVdW[], /* Get the i charge group info */ igid = GET_CGINFO_GID(cginfo[icg]); - iwater = (solvent_opt != esolNO) ? GET_CGINFO_SOLOPT(cginfo[icg]) : esolNO; - if (md->nPerturbed) { gmx_fatal(FARGS, "AdResS does not support free energy pertubation\n"); @@ -1585,18 +1577,18 @@ static real calc_image_tric(rvec xi, rvec xj, matrix box, /* Perform NINT operation, using trunc operation, therefore * we first add 2.5 then subtract 2 again */ - tz = dz*b_inv[ZZ] + h25; + tz = static_cast(dz*b_inv[ZZ] + h25); tz -= 2; dz -= tz*box[ZZ][ZZ]; dy -= tz*box[ZZ][YY]; dx -= tz*box[ZZ][XX]; - ty = dy*b_inv[YY] + h25; + ty = static_cast(dy*b_inv[YY] + h25); ty -= 2; dy -= ty*box[YY][YY]; dx -= ty*box[YY][XX]; - tx = dx*b_inv[XX]+h25; + tx = static_cast(dx*b_inv[XX]+h25); tx -= 2; dx -= tx*box[XX][XX]; @@ -1625,9 +1617,9 @@ static real calc_image_rect(rvec xi, rvec xj, rvec box_size, /* Perform NINT operation, using trunc operation, therefore * we first add 1.5 then subtract 1 again */ - tx = dx*b_inv[XX] + h15; - ty = dy*b_inv[YY] + h15; - tz = dz*b_inv[ZZ] + h15; + tx = static_cast(dx*b_inv[XX] + h15); + ty = static_cast(dy*b_inv[YY] + h15); + tz = static_cast(dz*b_inv[ZZ] + h15); tx--; ty--; tz--; @@ -1645,7 +1637,7 @@ static real calc_image_rect(rvec xi, rvec xj, rvec box_size, return r2; } -static void add_simple(t_ns_buf *nsbuf, int nrj, atom_id cg_j, +static void add_simple(t_ns_buf * nsbuf, int nrj, atom_id cg_j, gmx_bool bHaveVdW[], int ngid, t_mdatoms *md, int icg, int jgid, t_block *cgs, t_excl bexcl[], int shift, t_forcerec *fr, put_in_list_t *put_in_list) @@ -1673,7 +1665,6 @@ static void ns_inner_tric(rvec x[], int icg, int *i_egp_flags, int j, nrj, jgid; int *cginfo = fr->cginfo; atom_id cg_j, *cgindex; - t_ns_buf *nsbuf; cgindex = cgs->index; shift = CENTRAL; @@ -1706,7 +1697,6 @@ static void ns_inner_rect(rvec x[], int icg, int *i_egp_flags, int j, nrj, jgid; int *cginfo = fr->cginfo; atom_id cg_j, *cgindex; - t_ns_buf *nsbuf; cgindex = cgs->index; if (bBox) @@ -1760,7 +1750,7 @@ static int ns_simple_core(t_forcerec *fr, { int naaj, k; real rlist2; - int nsearch, icg, jcg, igid, i0, nri, nn; + int nsearch, icg, igid, nn; int *cginfo; t_ns_buf *nsbuf; /* atom_id *i_atoms; */ @@ -2044,8 +2034,8 @@ static void get_cutoff2(t_forcerec *fr, gmx_bool bDoLongRange, *rvdw2 = *rs2; *rcoul2 = *rs2; } - *rm2 = min(*rvdw2, *rcoul2); - *rl2 = max(*rvdw2, *rcoul2); + *rm2 = std::min(*rvdw2, *rcoul2); + *rl2 = std::max(*rvdw2, *rcoul2); } static void init_nsgrid_lists(t_forcerec *fr, int ngid, gmx_ns_t *ns) @@ -2107,15 +2097,15 @@ static int nsgrid_core(t_commrec *cr, t_forcerec *fr, #endif int dx0, dx1, dy0, dy1, dz0, dz1; int Nx, Ny, Nz, shift = -1, j, nrj, nns, nn = -1; - real gridx, gridy, gridz, grid_x, grid_y, grid_z; + real gridx, gridy, gridz, grid_x, grid_y; real *dcx2, *dcy2, *dcz2; int zgi, ygi, xgi; - int cg0, cg1, icg = -1, cgsnr, i0, igid, nri, naaj, max_jcg; + int cg0, cg1, icg = -1, cgsnr, i0, igid, naaj, max_jcg; int jcg0, jcg1, jjcg, cgj0, jgid; int *grida, *gridnra, *gridind; gmx_bool rvdw_lt_rcoul, rcoul_lt_rvdw; - rvec xi, *cgcm, grid_offset; - real r2, rs2, rvdw2, rcoul2, rm2, rl2, XI, YI, ZI, dcx, dcy, dcz, tmp1, tmp2; + rvec *cgcm, grid_offset; + real r2, rs2, rvdw2, rcoul2, rm2, rl2, XI, YI, ZI, tmp1, tmp2; int *i_egp_flags; gmx_bool bDomDec, bTriclinicX, bTriclinicY; ivec ncpddc; @@ -2167,7 +2157,6 @@ static int nsgrid_core(t_commrec *cr, t_forcerec *fr, gridz = grid->cell_size[ZZ]; grid_x = 1/gridx; grid_y = 1/gridy; - grid_z = 1/gridz; copy_rvec(grid->cell_offset, grid_offset); copy_ivec(grid->ncpddc, ncpddc); dcx2 = grid->dcx2; @@ -2213,7 +2202,7 @@ static int nsgrid_core(t_commrec *cr, t_forcerec *fr, else { if (d == XX && - box[XX][XX] - fabs(box[YY][XX]) - fabs(box[ZZ][XX]) < sqrt(rl2)) + box[XX][XX] - fabs(box[YY][XX]) - fabs(box[ZZ][XX]) < std::sqrt(rl2)) { shp[d] = 2; } @@ -2512,7 +2501,6 @@ static int nsgrid_core(t_commrec *cr, t_forcerec *fr, } } } - /* setexcl(nri,i_atoms,&top->atoms.excl,FALSE,bexcl); */ setexcl(cgs->index[icg], cgs->index[icg+1], &top->excls, FALSE, bexcl); } /* No need to perform any left-over force calculations anymore (as we used to do here) @@ -2548,7 +2536,6 @@ void init_ns(FILE *fplog, const t_commrec *cr, { int mt, icg, nr_in_cg, maxcg, i, j, jcg, ngid, ncg; t_block *cgs; - char *ptr; /* Compute largest charge groups size (# atoms) */ nr_in_cg = 1; @@ -2557,7 +2544,7 @@ void init_ns(FILE *fplog, const t_commrec *cr, cgs = &mtop->moltype[mt].cgs; for (icg = 0; (icg < cgs->nr); icg++) { - nr_in_cg = max(nr_in_cg, (int)(cgs->index[icg+1]-cgs->index[icg])); + nr_in_cg = std::max(nr_in_cg, (int)(cgs->index[icg+1]-cgs->index[icg])); } } @@ -2669,13 +2656,11 @@ int search_neighbours(FILE *log, t_forcerec *fr, { t_block *cgs = &(top->cgs); rvec box_size, grid_x0, grid_x1; - int i, j, m, ngid; + int m, ngid; real min_size, grid_dens; int nsearch; gmx_bool bGrid; - char *ptr; - gmx_bool *i_egp_flags; - int cg_start, cg_end, start, end; + int start, end; gmx_ns_t *ns; t_grid *grid; gmx_domdec_zones_t *dd_zones; @@ -2700,7 +2685,7 @@ int search_neighbours(FILE *log, t_forcerec *fr, } if (!bGrid) { - min_size = min(box_size[XX], min(box_size[YY], box_size[ZZ])); + min_size = std::min(box_size[XX], std::min(box_size[YY], box_size[ZZ])); if (2*fr->rlistlong >= min_size) { gmx_fatal(FARGS, "One of the box diagonal elements has become smaller than twice the cut-off length."); diff --git a/src/gromacs/mdlib/nsgrid.c b/src/gromacs/mdlib/nsgrid.cpp similarity index 97% rename from src/gromacs/mdlib/nsgrid.c rename to src/gromacs/mdlib/nsgrid.cpp index 75e689e8f4..b92d5a5c76 100644 --- a/src/gromacs/mdlib/nsgrid.c +++ b/src/gromacs/mdlib/nsgrid.cpp @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2013,2014, by the GROMACS development team, led by + * Copyright (c) 2013,2014,2015, 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. @@ -42,6 +42,10 @@ #include #include +#include + +#include + #include "gromacs/domdec/domdec.h" #include "gromacs/fileio/pdbio.h" #include "gromacs/legacyheaders/macros.h" @@ -90,7 +94,7 @@ static void calc_x_av_stddev(int n, rvec *x, rvec av, rvec stddev) for (d = 0; d < DIM; d++) { av[d] = s1[d]; - stddev[d] = sqrt(s2[d] - s1[d]*s1[d]); + stddev[d] = std::sqrt(s2[d] - s1[d]*s1[d]); } } @@ -202,7 +206,6 @@ static void set_grid_sizes(matrix box, rvec izones_x0, rvec izones_x1, real rlis { int i, j; gmx_bool bDD, bDDRect; - rvec av, stddev; rvec izones_size; real inv_r_ideal, size, add_tric, radd; @@ -218,7 +221,7 @@ static void set_grid_sizes(matrix box, rvec izones_x0, rvec izones_x1, real rlis } /* Use the ideal number of cg's per cell to set the ideal cell size */ - inv_r_ideal = pow(grid_density/grid->ncg_ideal, 1.0/3.0); + inv_r_ideal = std::pow((real)(grid_density/grid->ncg_ideal), (real)(1.0/3.0)); if (rlist > 0 && inv_r_ideal*rlist < 1) { inv_r_ideal = 1/rlist; @@ -273,8 +276,11 @@ static void set_grid_sizes(matrix box, rvec izones_x0, rvec izones_x1, real rlis /* Check if the cell boundary in this direction is * perpendicular to the Cartesian axis. + * Since grid->npbcdim isan integer that in principle can take + * any value, we help the compiler avoid warnings and potentially + * optimize by ensuring that j < DIM here. */ - for (j = i+1; j < grid->npbcdim; j++) + for (j = i+1; j < grid->npbcdim && j < DIM; j++) { if (box[j][i] != 0) { @@ -357,7 +363,6 @@ static void set_grid_sizes(matrix box, rvec izones_x0, rvec izones_x1, real rlis t_grid *init_grid(FILE *fplog, t_forcerec *fr) { - int d, m; char *ptr; t_grid *grid; @@ -475,7 +480,6 @@ void grid_first(FILE *fplog, t_grid *grid, real rlistlong, real grid_density) { int i, m; - ivec cx; set_grid_sizes(box, izones_x0, izones_x1, rlistlong, dd, ddbox, grid, grid_density); @@ -495,7 +499,7 @@ void grid_first(FILE *fplog, t_grid *grid, } } - m = max(grid->n[XX], max(grid->n[YY], grid->n[ZZ])); + m = std::max(grid->n[XX], std::max(grid->n[YY], grid->n[ZZ])); if (m > grid->dc_nalloc) { /* Allocate with double the initial size for box scaling */ @@ -642,7 +646,7 @@ void fill_grid(gmx_domdec_zones_t *dd_zones, int cg0, int cg1, rvec cg_cm[]) { int *cell_index; - int nrx, nry, nrz; + int nry, nrz; rvec n_box, offset; int zone, ccg0, ccg1, cg, d, not_used; ivec shift0, useall, b0, b1, ind; @@ -661,7 +665,6 @@ void fill_grid(gmx_domdec_zones_t *dd_zones, cell_index = grid->cell_index; /* Initiate cell borders */ - nrx = grid->n[XX]; nry = grid->n[YY]; nrz = grid->n[ZZ]; for (d = 0; d < DIM; d++) @@ -689,7 +692,7 @@ void fill_grid(gmx_domdec_zones_t *dd_zones, { for (d = 0; d < DIM; d++) { - ind[d] = (cg_cm[cg][d] - offset[d])*n_box[d]; + ind[d] = static_cast((cg_cm[cg][d] - offset[d])*n_box[d]); /* With pbc we should be done here. * Without pbc cg's outside the grid * should be assigned to the closest grid cell. @@ -764,7 +767,7 @@ void fill_grid(gmx_domdec_zones_t *dd_zones, bUse = TRUE; for (d = 0; d < DIM; d++) { - ind[d] = (cg_cm[cg][d] - offset[d])*n_box[d]; + ind[d] = static_cast((cg_cm[cg][d] - offset[d])*n_box[d]); /* Here we have to correct for rounding problems, * as this cg_cm to cell index operation is not necessarily * binary identical to the operation for the DD zone assignment