From 6291709ae701159646c4b5967c4bc6e8034685c5 Mon Sep 17 00:00:00 2001 From: Rafael Ravedutti Date: Tue, 25 Jan 2022 00:43:10 +0100 Subject: [PATCH] Add first draft code with GROMACS approach Signed-off-by: Rafael Ravedutti --- gromacs/atom.c | 14 ++ gromacs/includes/atom.h | 24 +++ gromacs/includes/neighbor.h | 3 +- gromacs/neighbor.c | 416 ++++++++++++++++++++++-------------- gromacs/pbc.c | 127 ++++++----- 5 files changed, 368 insertions(+), 216 deletions(-) diff --git a/gromacs/atom.c b/gromacs/atom.c index 5d41cd2..3c63821 100644 --- a/gromacs/atom.c +++ b/gromacs/atom.c @@ -46,16 +46,22 @@ void initAtom(Atom *atom) atom->x = NULL; atom->y = NULL; atom->z = NULL; atom->vx = NULL; atom->vy = NULL; atom->vz = NULL; atom->fx = NULL; atom->fy = NULL; atom->fz = NULL; + atom->cl_x = NULL; atom->Natoms = 0; atom->Nlocal = 0; atom->Nghost = 0; atom->Nmax = 0; + atom->Nclusters = 0; + atom->Nclusters_local = 0; + atom->Nclusters_ghost = 0; + atom->Nclusters_max = 0; atom->type = NULL; atom->ntypes = 0; atom->epsilon = NULL; atom->sigma6 = NULL; atom->cutforcesq = NULL; atom->cutneighsq = NULL; + atom->clusters = NULL; } void createAtom(Atom *atom, Parameter *param) @@ -274,3 +280,11 @@ void growAtom(Atom *atom) atom->fz = (MD_FLOAT*) reallocate(atom->fz, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT)); atom->type = (int *) reallocate(atom->type, ALIGNMENT, atom->Nmax * sizeof(int), nold * sizeof(int)); } + +void growClusters(Atom *atom) +{ + int nold = atom->Nclusters_max; + atom->Nclusters_max += DELTA; + atom->clusters = (Cluster*) reallocate(atom->clusters, ALIGNMENT, atom->Nclusters_max * sizeof(Cluster), nold * sizeof(Cluster)); + atom->cl_x = (MD_FLOAT*) reallocate(atom->cl_x, ALIGNMENT, atom->Nclusters_max * CLUSTER_DIM_N * 3, nold * CLUSTER_DIM_N * 3); +} diff --git a/gromacs/includes/atom.h b/gromacs/includes/atom.h index 70cef10..1978dbb 100644 --- a/gromacs/includes/atom.h +++ b/gromacs/includes/atom.h @@ -25,11 +25,25 @@ #ifndef __ATOM_H_ #define __ATOM_H_ +#define CLUSTER_DIM_N 4 +#define CLUSTER_DIM_M 4 + +typedef struct { + int bin; + int natoms; + int type[CLUSTER_DIM_N]; + MD_FLOAT bbminx, bbmaxx; + MD_FLOAT bbminy, bbmaxy; + MD_FLOAT bbminz, bbmaxz; +} Cluster; + typedef struct { int Natoms, Nlocal, Nghost, Nmax; + int Nclusters, Nclusters_local, Nclusters_ghost, Nclusters_max; MD_FLOAT *x, *y, *z; MD_FLOAT *vx, *vy, *vz; MD_FLOAT *fx, *fy, *fz; + MD_FLOAT *cl_x; int *border_map; int *type; int ntypes; @@ -37,23 +51,33 @@ typedef struct { MD_FLOAT *sigma6; MD_FLOAT *cutforcesq; MD_FLOAT *cutneighsq; + Cluster *clusters; } Atom; extern void initAtom(Atom*); extern void createAtom(Atom*, Parameter*); extern int readAtom(Atom*, Parameter*); extern void growAtom(Atom*); +extern void growClusters(Atom*); + +#define cluster_ptr(ci) &(atom->cl_x[ci * CLUSTER_DIM_N * 3]) #ifdef AOS #define POS_DATA_LAYOUT "AoS" #define atom_x(i) atom->x[(i) * 3 + 0] #define atom_y(i) atom->x[(i) * 3 + 1] #define atom_z(i) atom->x[(i) * 3 + 2] +#define cluster_x(cptr, i) cptr[(i) * 3 + 0] +#define cluster_y(cptr, i) cptr[(i) * 3 + 1] +#define cluster_z(cptr, i) cptr[(i) * 3 + 2] #else #define POS_DATA_LAYOUT "SoA" #define atom_x(i) atom->x[i] #define atom_y(i) atom->y[i] #define atom_z(i) atom->z[i] +#define cluster_x(cptr, i) cptr[0 * CLUSTER_DIM_N + (i)] +#define cluster_y(cptr, i) cptr[1 * CLUSTER_DIM_N + (i)] +#define cluster_z(cptr, i) cptr[2 * CLUSTER_DIM_N + (i)] #endif #endif diff --git a/gromacs/includes/neighbor.h b/gromacs/includes/neighbor.h index 949d515..0509f36 100644 --- a/gromacs/includes/neighbor.h +++ b/gromacs/includes/neighbor.h @@ -36,6 +36,7 @@ typedef struct { extern void initNeighbor(Neighbor*, Parameter*); extern void setupNeighbor(Parameter*); extern void binatoms(Atom*); -extern void buildNeighbor(Atom*, Neighbor*); +extern void buildNeighbor(Parameter*, Atom*, Neighbor*); extern void sortAtom(Atom*); +extern void buildClusters(Parameter*, Atom*); #endif diff --git a/gromacs/neighbor.c b/gromacs/neighbor.c index 6b512a0..a45f284 100644 --- a/gromacs/neighbor.c +++ b/gromacs/neighbor.c @@ -27,28 +27,34 @@ #include #include #include +#include #define SMALL 1.0e-6 #define FACTOR 0.999 -static MD_FLOAT xprd, yprd, zprd; -static MD_FLOAT bininvx, bininvy, bininvz; -static int mbinxlo, mbinylo, mbinzlo; -static int nbinx, nbiny, nbinz; -static int mbinx, mbiny, mbinz; // n bins in x, y, z +static MD_FLOAT xprd, yprd; +static MD_FLOAT bininvx, bininvy; +static int mbinxlo, mbinylo; +static int nbinx, nbiny; +static int mbinx, mbiny; // n bins in x, y static int *bincount; static int *bins; +static int *cluster_bincount; +static int *cluster_bins; static int mbins; //total number of bins static int atoms_per_bin; // max atoms per bin +static int clusters_per_bin; // max clusters per bin static MD_FLOAT cutneigh; static MD_FLOAT cutneighsq; // neighbor cutoff squared static int nmax; static int nstencil; // # of bins in stencil static int* stencil; // stencil list of bin offsets -static MD_FLOAT binsizex, binsizey, binsizez; +static MD_FLOAT binsizex, binsizey; -static int coord2bin(MD_FLOAT, MD_FLOAT , MD_FLOAT); -static MD_FLOAT bindist(int, int, int); +static int coord2bin(MD_FLOAT, MD_FLOAT); +static MD_FLOAT bindist(int, int); + +extern int *PBCx, *PBCy, *PBCz; /* exported subroutines */ void initNeighbor(Neighbor *neighbor, Parameter *param) @@ -56,60 +62,52 @@ void initNeighbor(Neighbor *neighbor, Parameter *param) MD_FLOAT neighscale = 5.0 / 6.0; xprd = param->nx * param->lattice; yprd = param->ny * param->lattice; - zprd = param->nz * param->lattice; cutneigh = param->cutneigh; nbinx = neighscale * param->nx; nbiny = neighscale * param->ny; - nbinz = neighscale * param->nz; nmax = 0; atoms_per_bin = 8; + clusters_per_bin = (atoms_per_bin / CLUSTER_DIM_N) + 4; stencil = NULL; bins = NULL; bincount = NULL; + cluster_bins = NULL; + cluster_bincount = NULL; neighbor->maxneighs = 100; neighbor->numneigh = NULL; neighbor->neighbors = NULL; } -void setupNeighbor(Parameter* param) -{ +void setupNeighbor(Parameter* param) { MD_FLOAT coord; - int mbinxhi, mbinyhi, mbinzhi; + int mbinxhi, mbinyhi; int nextx, nexty, nextz; if(param->input_file != NULL) { xprd = param->xprd; yprd = param->yprd; - zprd = param->zprd; } // TODO: update lo and hi for standard case and use them here instead MD_FLOAT xlo = 0.0; MD_FLOAT xhi = xprd; MD_FLOAT ylo = 0.0; MD_FLOAT yhi = yprd; - MD_FLOAT zlo = 0.0; MD_FLOAT zhi = zprd; cutneighsq = cutneigh * cutneigh; if(param->input_file != NULL) { binsizex = cutneigh * 0.5; binsizey = cutneigh * 0.5; - binsizez = cutneigh * 0.5; nbinx = (int)((param->xhi - param->xlo) / binsizex); nbiny = (int)((param->yhi - param->ylo) / binsizey); - nbinz = (int)((param->zhi - param->zlo) / binsizez); if(nbinx == 0) { nbinx = 1; } if(nbiny == 0) { nbiny = 1; } - if(nbinz == 0) { nbinz = 1; } bininvx = nbinx / (param->xhi - param->xlo); bininvy = nbiny / (param->yhi - param->ylo); - bininvz = nbinz / (param->zhi - param->zlo); } else { binsizex = xprd / nbinx; binsizey = yprd / nbiny; - binsizez = zprd / nbinz; bininvx = 1.0 / binsizex; bininvy = 1.0 / binsizey; - bininvz = 1.0 / binsizez; } coord = xlo - cutneigh - SMALL * xprd; @@ -128,14 +126,6 @@ void setupNeighbor(Parameter* param) coord = yhi + cutneigh + SMALL * yprd; mbinyhi = (int) (coord * bininvy); - coord = zlo - cutneigh - SMALL * zprd; - mbinzlo = (int) (coord * bininvz); - if (coord < 0.0) { - mbinzlo = mbinzlo - 1; - } - coord = zhi + cutneigh + SMALL * zprd; - mbinzhi = (int) (coord * bininvz); - mbinxlo = mbinxlo - 1; mbinxhi = mbinxhi + 1; mbinx = mbinxhi - mbinxlo + 1; @@ -144,55 +134,82 @@ void setupNeighbor(Parameter* param) mbinyhi = mbinyhi + 1; mbiny = mbinyhi - mbinylo + 1; - mbinzlo = mbinzlo - 1; - mbinzhi = mbinzhi + 1; - mbinz = mbinzhi - mbinzlo + 1; - nextx = (int) (cutneigh * bininvx); if(nextx * binsizex < FACTOR * cutneigh) nextx++; nexty = (int) (cutneigh * bininvy); if(nexty * binsizey < FACTOR * cutneigh) nexty++; - nextz = (int) (cutneigh * bininvz); - if(nextz * binsizez < FACTOR * cutneigh) nextz++; - if (stencil) { free(stencil); } - stencil = (int*) malloc( - (2 * nextz + 1) * (2 * nexty + 1) * (2 * nextx + 1) * sizeof(int)); - + stencil = (int *) malloc((2 * nexty + 1) * (2 * nextx + 1) * sizeof(int)); nstencil = 0; - int kstart = -nextz; - for(int k = kstart; k <= nextz; k++) { - for(int j = -nexty; j <= nexty; j++) { - for(int i = -nextx; i <= nextx; i++) { - if(bindist(i, j, k) < cutneighsq) { - stencil[nstencil++] = - k * mbiny * mbinx + j * mbinx + i; - } + for(int j = -nexty; j <= nexty; j++) { + for(int i = -nextx; i <= nextx; i++) { + if(bindist(i, j) < cutneighsq) { + stencil[nstencil++] = j * mbinx + i; } } } - mbins = mbinx * mbiny * mbinz; + mbins = mbinx * mbiny; - if (bincount) { - free(bincount); - } + if (bincount) { free(bincount); } bincount = (int*) malloc(mbins * sizeof(int)); - if (bins) { - free(bins); - } + if (bins) { free(bins); } bins = (int*) malloc(mbins * atoms_per_bin * sizeof(int)); + + if (cluster_bincount) { free(cluster_bincount); } + cluster_bincount = (int*) malloc(mbins * sizeof(int)); + + if (cluster_bins) { free(cluster_bins); } + cluster_bins = (int*) malloc(mbins * clusters_per_bin * sizeof(int)); } -void buildNeighbor(Atom *atom, Neighbor *neighbor) -{ +MD_FLOAT getBoundingBoxDistanceSq(Atom *atom, int ci, int cj) { + MD_FLOAT dl = atom->clusters[ci].bbminx - atom->clusters[cj].bbmaxx; + MD_FLOAT dh = atom->clusters[cj].bbminx - atom->clusters[ci].bbmaxx; + MD_FLOAT dm = MAX(dl, dh); + MD_FLOAT dm0 = MAX(dm, 0.0); + MD_FLOAT d2 = dm0 * dm0; + + dl = atom->clusters[ci].bbminy - atom->clusters[cj].bbmaxy; + dh = atom->clusters[cj].bbminy - atom->clusters[ci].bbmaxy; + dm = MAX(dl, dh); + dm0 = MAX(dm, 0.0); + d2 += dm0 * dm0; + + dl = atom->clusters[ci].bbminz - atom->clusters[cj].bbmaxz; + dh = atom->clusters[cj].bbminz - atom->clusters[ci].bbmaxz; + dm = MAX(dl, dh); + dm0 = MAX(dm, 0.0); + d2 += dm0 * dm0; + return d2; +} + +int atomDistanceInRange(Atom *atom, int ci, int cj, MD_FLOAT rsq) { + MD_FLOAT *ciptr = cluster_ptr(ci); + MD_FLOAT *cjptr = cluster_ptr(cj); + + for(int cii = 0; cii < atom->clusters[ci].natoms; cii++) { + for(int cjj = 0; cjj < atom->clusters[cj].natoms; cjj++) { + MD_FLOAT delx = cluster_x(ciptr, cii) - cluster_x(cjptr, cjj); + MD_FLOAT dely = cluster_y(ciptr, cii) - cluster_y(cjptr, cjj); + MD_FLOAT delz = cluster_z(ciptr, cii) - cluster_z(cjptr, cjj); + if(delx * delx + dely * dely + delz * delz < rsq) { + return 1; + } + } + } + + return 0; +} + +void buildNeighbor(Parameter *param, Atom *atom, Neighbor *neighbor) { int nall = atom->Nlocal + atom->Nghost; /* extend atom arrays if necessary */ @@ -206,6 +223,9 @@ void buildNeighbor(Atom *atom, Neighbor *neighbor) /* bin local & ghost atoms */ binatoms(atom); + buildClusters(param, atom); + + const MD_FLOAT rBB = cutneighsq / 2.0; // TODO: change this int resize = 1; /* loop over each atom, storing neighbors */ @@ -213,45 +233,49 @@ void buildNeighbor(Atom *atom, Neighbor *neighbor) int new_maxneighs = neighbor->maxneighs; resize = 0; - for(int i = 0; i < atom->Nlocal; i++) { - int* neighptr = &(neighbor->neighbors[i * neighbor->maxneighs]); + for(int ci = 0; ci < atom->Nclusters_local; ci++) { + int* neighptr = &(neighbor->neighbors[ci * neighbor->maxneighs]); int n = 0; - MD_FLOAT xtmp = atom_x(i); - MD_FLOAT ytmp = atom_y(i); - MD_FLOAT ztmp = atom_z(i); - int ibin = coord2bin(xtmp, ytmp, ztmp); - #ifdef EXPLICIT_TYPES - int type_i = atom->type[i]; - #endif + int ibin = atom->clusters[ci].bin; + MD_FLOAT ibb_zmin = atom->clusters[ci].bbminz; + MD_FLOAT ibb_zmax = atom->clusters[ci].bbmaxz; + for(int k = 0; k < nstencil; k++) { int jbin = ibin + stencil[k]; - int* loc_bin = &bins[jbin * atoms_per_bin]; + int *loc_bin = &cluster_bins[jbin * clusters_per_bin]; + int cj, m = -1; + MD_FLOAT jbb_zmin, jbb_zmax; - for(int m = 0; m < bincount[jbin]; m++) { - int j = loc_bin[m]; + do { + m++; + cj = loc_bin[m]; + jbb_zmin = atom->clusters[cj].bbminz; + jbb_zmax = atom->clusters[cj].bbmaxz; + } while((ibb_zmin - jbb_zmax) * (ibb_zmin - jbb_zmax) > cutneighsq); - if ( j == i ){ - continue; + while(m < cluster_bincount[jbin] && jbb_zmax - ibb_zmax < cutneighsq) { + if((jbb_zmin - ibb_zmax) * (jbb_zmin - ibb_zmax) > cutneighsq) { + break; } - MD_FLOAT delx = xtmp - atom_x(j); - MD_FLOAT dely = ytmp - atom_y(j); - MD_FLOAT delz = ztmp - atom_z(j); - MD_FLOAT rsq = delx * delx + dely * dely + delz * delz; - - #ifdef EXPLICIT_TYPES - int type_j = atom->type[j]; - const MD_FLOAT cutoff = atom->cutneighsq[type_i * atom->ntypes + type_j]; - #else - const MD_FLOAT cutoff = cutneighsq; - #endif - if( rsq <= cutoff ) { - neighptr[n++] = j; + double d_bb_sq = getBoundingBoxDistanceSq(atom, ci, cj); + if(d_bb_sq < cutneighsq) { + if(d_bb_sq < rBB || atomDistanceInRange(atom, ci, cj, cutneighsq)) { + if(cj == ci) { + // Add to neighbor list with diagonal interaction mask + } + neighptr[n++] = cj; + } } + + m++; + cj = loc_bin[m]; + jbb_zmin = atom->clusters[cj].bbminz; + jbb_zmax = atom->clusters[cj].bbmaxz; } } - neighbor->numneigh[i] = n; + neighbor->numneigh[ci] = n; if(n >= neighbor->maxneighs) { resize = 1; @@ -272,8 +296,7 @@ void buildNeighbor(Atom *atom, Neighbor *neighbor) } /* internal subroutines */ -MD_FLOAT bindist(int i, int j, int k) -{ +MD_FLOAT bindist(int i, int j) { MD_FLOAT delx, dely, delz; if(i > 0) { @@ -292,20 +315,11 @@ MD_FLOAT bindist(int i, int j, int k) dely = (j + 1) * binsizey; } - if(k > 0) { - delz = (k - 1) * binsizez; - } else if(k == 0) { - delz = 0.0; - } else { - delz = (k + 1) * binsizez; - } - - return (delx * delx + dely * dely + delz * delz); + return (delx * delx + dely * dely); } -int coord2bin(MD_FLOAT xin, MD_FLOAT yin, MD_FLOAT zin) -{ - int ix, iy, iz; +int coord2bin(MD_FLOAT xin, MD_FLOAT yin) { + int ix, iy; if(xin >= xprd) { ix = (int)((xin - xprd) * bininvx) + nbinx - mbinxlo; @@ -323,19 +337,28 @@ int coord2bin(MD_FLOAT xin, MD_FLOAT yin, MD_FLOAT zin) iy = (int)(yin * bininvy) - mbinylo - 1; } - if(zin >= zprd) { - iz = (int)((zin - zprd) * bininvz) + nbinz - mbinzlo; - } else if(zin >= 0.0) { - iz = (int)(zin * bininvz) - mbinzlo; - } else { - iz = (int)(zin * bininvz) - mbinzlo - 1; - } - - return (iz * mbiny * mbinx + iy * mbinx + ix + 1); + return (iy * mbinx + ix + 1); } -void binatoms(Atom *atom) -{ +void coord2bin2D(MD_FLOAT xin, MD_FLOAT yin, int *ix, int *iy) { + if(xin >= xprd) { + *ix = (int)((xin - xprd) * bininvx) + nbinx - mbinxlo; + } else if(xin >= 0.0) { + *ix = (int)(xin * bininvx) - mbinxlo; + } else { + *ix = (int)(xin * bininvx) - mbinxlo - 1; + } + + if(yin >= yprd) { + *iy = (int)((yin - yprd) * bininvy) + nbiny - mbinylo; + } else if(yin >= 0.0) { + *iy = (int)(yin * bininvy) - mbinylo; + } else { + *iy = (int)(yin * bininvy) - mbinylo - 1; + } +} + +void binatoms(Atom *atom) { int nall = atom->Nlocal + atom->Nghost; int resize = 1; @@ -347,7 +370,7 @@ void binatoms(Atom *atom) } for(int i = 0; i < nall; i++) { - int ibin = coord2bin(atom_x(i), atom_y(i), atom_z(i)); + int ibin = coord2bin(atom_x(i), atom_y(i)); if(bincount[ibin] < atoms_per_bin) { int ac = bincount[ibin]++; @@ -365,56 +388,137 @@ void binatoms(Atom *atom) } } -void sortAtom(Atom* atom) { - binatoms(atom); - int Nmax = atom->Nmax; - int* binpos = bincount; +// TODO: Use pigeonhole sorting +void sortBinAtomsByZCoord(Parameter *param, Atom *atom, int bin) { + int c = bincount[bin]; + int *bin_ptr = &bins[bin * atoms_per_bin]; - for(int i=1; ix; double* old_y = atom->y; double* old_z = atom->z; - double* old_vx = atom->vx; double* old_vy = atom->vy; double* old_vz = atom->vz; + int resize = 1; + while(resize > 0) { + resize = 0; + for(int bin = 0; bin < mbins; bin++) { + int c = bincount[bin]; + int ac = 0; + const int nclusters = ((c + CLUSTER_DIM_N - 1) / CLUSTER_DIM_N); - for(int mybin = 0; mybin0?binpos[mybin-1]:0; - int count = binpos[mybin] - start; - for(int k=0; kNclusters_local; + MD_FLOAT *cptr = cluster_ptr(ci); + MD_FLOAT bbminx = INFINITY, bbmaxx = -INFINITY; + MD_FLOAT bbminy = INFINITY, bbmaxy = -INFINITY; + MD_FLOAT bbminz = INFINITY, bbmaxz = -INFINITY; + atom->clusters[ci].natoms = 0; + + for(int cii = 0; cii < CLUSTER_DIM_N; cii++) { + if(ac < c) { + int i = bins[bin * atoms_per_bin + ac]; + MD_FLOAT xtmp = atom_x(i); + MD_FLOAT ytmp = atom_y(i); + MD_FLOAT ztmp = atom_z(i); + + cluster_x(cptr, cii) = xtmp; + cluster_y(cptr, cii) = ytmp; + cluster_z(cptr, cii) = ztmp; + + // TODO: To create the bounding boxes faster, we can use SIMD operations + if(bbminx > xtmp) { bbminx = xtmp; } + if(bbmaxx < xtmp) { bbmaxx = xtmp; } + if(bbminy > ytmp) { bbminy = ytmp; } + if(bbmaxy < ytmp) { bbmaxy = ytmp; } + if(bbminz > ytmp) { bbminz = ytmp; } + if(bbmaxz < ytmp) { bbmaxz = ytmp; } + + atom->clusters[ci].type[cii] = atom->type[i]; + atom->clusters[ci].natoms++; + } else { + cluster_x(cptr, cii) = INFINITY; + cluster_y(cptr, cii) = INFINITY; + cluster_z(cptr, cii) = INFINITY; + } + + ac++; + } + + cluster_bins[bin * clusters_per_bin + cl] = ci; + atom->clusters[ci].bin = bin; + atom->clusters[ci].bbminx = bbminx; + atom->clusters[ci].bbmaxx = bbmaxx; + atom->clusters[ci].bbminy = bbminy; + atom->clusters[ci].bbmaxy = bbmaxy; + atom->clusters[ci].bbminz = bbminz; + atom->clusters[ci].bbmaxz = bbmaxz; + atom->Nclusters_local++; + } + + cluster_bincount[bin] = nclusters; + } else { + resize = 1; + } + } + + if(resize) { + free(cluster_bins); + clusters_per_bin *= 2; + cluster_bins = (int*) malloc(mbins * clusters_per_bin * sizeof(int)); } } - - free(atom->x); - atom->x = new_x; -#ifndef AOS - free(atom->y); - free(atom->z); - atom->y = new_y; atom->z = new_z; -#endif - free(atom->vx); free(atom->vy); free(atom->vz); - atom->vx = new_vx; atom->vy = new_vy; atom->vz = new_vz; +} + +void binGhostClusters(Parameter *param, Atom *atom) { + int nlocal = atom->Nclusters_local; + + for(int ci = 0; ci < atom->Nclusters_ghost; ci++) { + MD_FLOAT *cptr = cluster_ptr(nlocal + ci); + MD_FLOAT xtmp, ytmp; + int ix = -1, iy = -1; + + xtmp = cluster_x(cptr, 0); + ytmp = cluster_y(cptr, 0); + coord2bin2D(xtmp, ytmp, &ix, &iy); + + for(int cii = 1; cii < atom->clusters[ci].natoms; cii++) { + int nix, niy; + xtmp = cluster_x(cptr, cii); + ytmp = cluster_y(cptr, cii); + coord2bin2D(xtmp, ytmp, &nix, &niy); + // Always put the cluster on the bin of its innermost atom so + // the cluster should be closer to local clusters + if(PBCx[ci] > 0 && ix > nix) { ix = nix; } + if(PBCx[ci] < 0 && ix < nix) { ix = nix; } + if(PBCy[ci] > 0 && iy > niy) { iy = niy; } + if(PBCy[ci] < 0 && iy < niy) { iy = niy; } + } + + int bin = iy * mbinx + ix + 1; + int c = cluster_bincount[bin]; + cluster_bins[bin * clusters_per_bin + c] = nlocal + ci; + cluster_bincount[bin]++; + } } diff --git a/gromacs/pbc.c b/gromacs/pbc.c index 3610701..24db091 100644 --- a/gromacs/pbc.c +++ b/gromacs/pbc.c @@ -35,8 +35,7 @@ static int *PBCx, *PBCy, *PBCz; static void growPbc(Atom*); /* exported subroutines */ -void initPbc(Atom* atom) -{ +void initPbc(Atom* atom) { NmaxGhost = 0; atom->border_map = NULL; PBCx = NULL; PBCy = NULL; PBCz = NULL; @@ -44,36 +43,37 @@ void initPbc(Atom* atom) /* update coordinates of ghost atoms */ /* uses mapping created in setupPbc */ -void updatePbc(Atom *atom, Parameter *param) -{ +void updatePbc(Atom *atom, Parameter *param) { int *border_map = atom->border_map; - int nlocal = atom->Nlocal; + int nlocal = atom->Nclusters_local; MD_FLOAT xprd = param->xprd; MD_FLOAT yprd = param->yprd; MD_FLOAT zprd = param->zprd; - for(int i = 0; i < atom->Nghost; i++) { - atom_x(nlocal + i) = atom_x(border_map[i]) + PBCx[i] * xprd; - atom_y(nlocal + i) = atom_y(border_map[i]) + PBCy[i] * yprd; - atom_z(nlocal + i) = atom_z(border_map[i]) + PBCz[i] * zprd; + for(int ci = 0; ci < atom->Nclusters_ghost; ci++) { + MD_FLOAT *cptr = cluster_ptr(nlocal + ci); + MD_FLOAT *bmap_cptr = cluster_ptr(border_map[ci]); + + for(int cii = 0; cii < atom->clusters[ci].natoms; cii++) { + cluster_x(cptr, cii) = cluster_x(bmap_cptr, cii) + PBCx[ci] * xprd; + cluster_y(cptr, cii) = cluster_y(bmap_cptr, cii) + PBCy[ci] * yprd; + cluster_z(cptr, cii) = cluster_z(bmap_cptr, cii) + PBCz[ci] * zprd; + } } } /* relocate atoms that have left domain according * to periodic boundary conditions */ -void updateAtomsPbc(Atom *atom, Parameter *param) -{ +void updateAtomsPbc(Atom *atom, Parameter *param) { MD_FLOAT xprd = param->xprd; MD_FLOAT yprd = param->yprd; MD_FLOAT zprd = param->zprd; for(int i = 0; i < atom->Nlocal; i++) { - if(atom_x(i) < 0.0) { atom_x(i) += xprd; } else if(atom_x(i) >= xprd) { atom_x(i) -= xprd; - } if(atom_y(i) < 0.0) { atom_y(i) += yprd; @@ -93,16 +93,21 @@ void updateAtomsPbc(Atom *atom, Parameter *param) * defining ghost atoms around domain * only creates mapping and coordinate corrections * that are then enforced in updatePbc */ -#define ADDGHOST(dx,dy,dz) \ - Nghost++; \ - border_map[Nghost] = i; \ - PBCx[Nghost] = dx; \ - PBCy[Nghost] = dy; \ - PBCz[Nghost] = dz; \ - atom->type[atom->Nlocal + Nghost] = atom->type[i] +#define ADDGHOST(dx,dy,dz) \ + Nghost++; \ + border_map[Nghost] = ci; \ + PBCx[Nghost] = dx; \ + PBCy[Nghost] = dy; \ + PBCz[Nghost] = dz; \ + copy_cluster_types(atom, atom->Nclusters_local + Nghost, ci) -void setupPbc(Atom *atom, Parameter *param) -{ +void copy_cluster_types(Atom *atom, int dest, int src) { + for(int cii = 0; cii < atom->clusters[src].natoms; cii++) { + atom->clusters[dest].type[cii] = atom->clusters[src].type[cii]; + } +} + +void setupPbc(Atom *atom, Parameter *param) { int *border_map = atom->border_map; MD_FLOAT xprd = param->xprd; MD_FLOAT yprd = param->yprd; @@ -110,58 +115,62 @@ void setupPbc(Atom *atom, Parameter *param) MD_FLOAT Cutneigh = param->cutneigh; int Nghost = -1; - for(int i = 0; i < atom->Nlocal; i++) { - - if (atom->Nlocal + Nghost + 7 >= atom->Nmax) { - growAtom(atom); + for(int ci = 0; ci < atom->Nclusters_local; ci++) { + if (atom->Nclusters_local + atom->Nclusters_ghost + 7 >= atom->Nclusters_max) { + growClusters(atom); } - if (Nghost + 7 >= NmaxGhost) { + + if (atom->Nclusters_ghost + 7 >= NmaxGhost) { growPbc(atom); border_map = atom->border_map; } - MD_FLOAT x = atom_x(i); - MD_FLOAT y = atom_y(i); - MD_FLOAT z = atom_z(i); + MD_FLOAT bbminx = atom->clusters[ci].bbminx; + MD_FLOAT bbmaxx = atom->clusters[ci].bbmaxx; + MD_FLOAT bbminy = atom->clusters[ci].bbminy; + MD_FLOAT bbmaxy = atom->clusters[ci].bbmaxy; + MD_FLOAT bbminz = atom->clusters[ci].bbminz; + MD_FLOAT bbmaxz = atom->clusters[ci].bbmaxz; /* Setup ghost atoms */ /* 6 planes */ - if (x < Cutneigh) { ADDGHOST(+1,0,0); } - if (x >= (xprd-Cutneigh)) { ADDGHOST(-1,0,0); } - if (y < Cutneigh) { ADDGHOST(0,+1,0); } - if (y >= (yprd-Cutneigh)) { ADDGHOST(0,-1,0); } - if (z < Cutneigh) { ADDGHOST(0,0,+1); } - if (z >= (zprd-Cutneigh)) { ADDGHOST(0,0,-1); } + if (bbminx < Cutneigh) { ADDGHOST(+1,0,0); } + if (bbmaxx >= (xprd-Cutneigh)) { ADDGHOST(-1,0,0); } + if (bbminy < Cutneigh) { ADDGHOST(0,+1,0); } + if (bbmaxy >= (yprd-Cutneigh)) { ADDGHOST(0,-1,0); } + if (bbminz < Cutneigh) { ADDGHOST(0,0,+1); } + if (bbmaxz >= (zprd-Cutneigh)) { ADDGHOST(0,0,-1); } /* 8 corners */ - if (x < Cutneigh && y < Cutneigh && z < Cutneigh) { ADDGHOST(+1,+1,+1); } - if (x < Cutneigh && y >= (yprd-Cutneigh) && z < Cutneigh) { ADDGHOST(+1,-1,+1); } - if (x < Cutneigh && y >= Cutneigh && z >= (zprd-Cutneigh)) { ADDGHOST(+1,+1,-1); } - if (x < Cutneigh && y >= (yprd-Cutneigh) && z >= (zprd-Cutneigh)) { ADDGHOST(+1,-1,-1); } - if (x >= (xprd-Cutneigh) && y < Cutneigh && z < Cutneigh) { ADDGHOST(-1,+1,+1); } - if (x >= (xprd-Cutneigh) && y >= (yprd-Cutneigh) && z < Cutneigh) { ADDGHOST(-1,-1,+1); } - if (x >= (xprd-Cutneigh) && y < Cutneigh && z >= (zprd-Cutneigh)) { ADDGHOST(-1,+1,-1); } - if (x >= (xprd-Cutneigh) && y >= (yprd-Cutneigh) && z >= (zprd-Cutneigh)) { ADDGHOST(-1,-1,-1); } + if (bbminx < Cutneigh && bbminy < Cutneigh && bbminz < Cutneigh) { ADDGHOST(+1,+1,+1); } + if (bbminx < Cutneigh && bbmaxy >= (yprd-Cutneigh) && bbminz < Cutneigh) { ADDGHOST(+1,-1,+1); } + if (bbminx < Cutneigh && bbminy < Cutneigh && bbmaxz >= (zprd-Cutneigh)) { ADDGHOST(+1,+1,-1); } + if (bbminx < Cutneigh && bbmaxy >= (yprd-Cutneigh) && bbmaxz >= (zprd-Cutneigh)) { ADDGHOST(+1,-1,-1); } + if (bbmaxx >= (xprd-Cutneigh) && bbminy < Cutneigh && bbminz < Cutneigh) { ADDGHOST(-1,+1,+1); } + if (bbmaxx >= (xprd-Cutneigh) && bbmaxy >= (yprd-Cutneigh) && bbminz < Cutneigh) { ADDGHOST(-1,-1,+1); } + if (bbmaxx >= (xprd-Cutneigh) && bbminy < Cutneigh && bbmaxz >= (zprd-Cutneigh)) { ADDGHOST(-1,+1,-1); } + if (bbmaxx >= (xprd-Cutneigh) && bbmaxy >= (yprd-Cutneigh) && bbmaxz >= (zprd-Cutneigh)) { ADDGHOST(-1,-1,-1); } /* 12 edges */ - if (x < Cutneigh && z < Cutneigh) { ADDGHOST(+1,0,+1); } - if (x < Cutneigh && z >= (zprd-Cutneigh)) { ADDGHOST(+1,0,-1); } - if (x >= (xprd-Cutneigh) && z < Cutneigh) { ADDGHOST(-1,0,+1); } - if (x >= (xprd-Cutneigh) && z >= (zprd-Cutneigh)) { ADDGHOST(-1,0,-1); } - if (y < Cutneigh && z < Cutneigh) { ADDGHOST(0,+1,+1); } - if (y < Cutneigh && z >= (zprd-Cutneigh)) { ADDGHOST(0,+1,-1); } - if (y >= (yprd-Cutneigh) && z < Cutneigh) { ADDGHOST(0,-1,+1); } - if (y >= (yprd-Cutneigh) && z >= (zprd-Cutneigh)) { ADDGHOST(0,-1,-1); } - if (y < Cutneigh && x < Cutneigh) { ADDGHOST(+1,+1,0); } - if (y < Cutneigh && x >= (xprd-Cutneigh)) { ADDGHOST(-1,+1,0); } - if (y >= (yprd-Cutneigh) && x < Cutneigh) { ADDGHOST(+1,-1,0); } - if (y >= (yprd-Cutneigh) && x >= (xprd-Cutneigh)) { ADDGHOST(-1,-1,0); } + if (bbminx < Cutneigh && bbminz < Cutneigh) { ADDGHOST(+1,0,+1); } + if (bbminx < Cutneigh && bbmaxz >= (zprd-Cutneigh)) { ADDGHOST(+1,0,-1); } + if (bbmaxx >= (xprd-Cutneigh) && bbminz < Cutneigh) { ADDGHOST(-1,0,+1); } + if (bbmaxx >= (xprd-Cutneigh) && bbmaxz >= (zprd-Cutneigh)) { ADDGHOST(-1,0,-1); } + if (bbminy < Cutneigh && bbminz < Cutneigh) { ADDGHOST(0,+1,+1); } + if (bbminy < Cutneigh && bbmaxz >= (zprd-Cutneigh)) { ADDGHOST(0,+1,-1); } + if (bbmaxy >= (yprd-Cutneigh) && bbminz < Cutneigh) { ADDGHOST(0,-1,+1); } + if (bbmaxy >= (yprd-Cutneigh) && bbmaxz >= (zprd-Cutneigh)) { ADDGHOST(0,-1,-1); } + if (bbminy < Cutneigh && bbminx < Cutneigh) { ADDGHOST(+1,+1,0); } + if (bbminy < Cutneigh && bbmaxx >= (xprd-Cutneigh)) { ADDGHOST(-1,+1,0); } + if (bbmaxy >= (yprd-Cutneigh) && bbminx < Cutneigh) { ADDGHOST(+1,-1,0); } + if (bbmaxy >= (yprd-Cutneigh) && bbmaxx >= (xprd-Cutneigh)) { ADDGHOST(-1,-1,0); } } + // increase by one to make it the ghost atom count - atom->Nghost = Nghost + 1; + atom->Nclusters_ghost = Nghost + 1; + atom->Nclusters = atom->Nclusters_local + Nghost + 1; } /* internal subroutines */ -void growPbc(Atom* atom) -{ +void growPbc(Atom* atom) { int nold = NmaxGhost; NmaxGhost += DELTA;