moebiusband/MD-Bench
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Add first draft code with GROMACS approach

Browse Source 
Signed-off-by: Rafael Ravedutti <rafaelravedutti@gmail.com>
This commit is contained in:
Rafael Ravedutti 2022-01-25 00:43:10 +01:00
parent 72730bc27b
commit 6291709ae7
5 changed files with 368 additions and 216 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

14
gromacs/atom.c
View File

@ -46,16 +46,22 @@ void initAtom(Atom *atom)
atom->x = NULL; atom->y = NULL; atom->z = NULL; atom->x = NULL; atom->y = NULL; atom->z = NULL;
atom->vx = NULL; atom->vy = NULL; atom->vz = NULL; atom->vx = NULL; atom->vy = NULL; atom->vz = NULL;
atom->fx = NULL; atom->fy = NULL; atom->fz = NULL; atom->fx = NULL; atom->fy = NULL; atom->fz = NULL;
atom->cl_x = NULL;
atom->Natoms = 0; atom->Natoms = 0;
atom->Nlocal = 0; atom->Nlocal = 0;
atom->Nghost = 0; atom->Nghost = 0;
atom->Nmax = 0; atom->Nmax = 0;
atom->Nclusters = 0;
atom->Nclusters_local = 0;
atom->Nclusters_ghost = 0;
atom->Nclusters_max = 0;
atom->type = NULL; atom->type = NULL;
atom->ntypes = 0; atom->ntypes = 0;
atom->epsilon = NULL; atom->epsilon = NULL;
atom->sigma6 = NULL; atom->sigma6 = NULL;
atom->cutforcesq = NULL; atom->cutforcesq = NULL;
atom->cutneighsq = NULL; atom->cutneighsq = NULL;
atom->clusters = NULL;
} }
void createAtom(Atom *atom, Parameter *param) 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->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)); 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);
}

24
gromacs/includes/atom.h
View File

@ -25,11 +25,25 @@
#ifndef __ATOM_H_ #ifndef __ATOM_H_
#define __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 { typedef struct {
int Natoms, Nlocal, Nghost, Nmax; int Natoms, Nlocal, Nghost, Nmax;
int Nclusters, Nclusters_local, Nclusters_ghost, Nclusters_max;
MD_FLOAT *x, *y, *z; MD_FLOAT *x, *y, *z;
MD_FLOAT *vx, *vy, *vz; MD_FLOAT *vx, *vy, *vz;
MD_FLOAT *fx, *fy, *fz; MD_FLOAT *fx, *fy, *fz;
MD_FLOAT *cl_x;
int *border_map; int *border_map;
int *type; int *type;
int ntypes; int ntypes;
@ -37,23 +51,33 @@ typedef struct {
MD_FLOAT *sigma6; MD_FLOAT *sigma6;
MD_FLOAT *cutforcesq; MD_FLOAT *cutforcesq;
MD_FLOAT *cutneighsq; MD_FLOAT *cutneighsq;
Cluster *clusters;
} Atom; } Atom;
extern void initAtom(Atom*); extern void initAtom(Atom*);
extern void createAtom(Atom*, Parameter*); extern void createAtom(Atom*, Parameter*);
extern int readAtom(Atom*, Parameter*); extern int readAtom(Atom*, Parameter*);
extern void growAtom(Atom*); extern void growAtom(Atom*);
extern void growClusters(Atom*);
#define cluster_ptr(ci) &(atom->cl_x[ci * CLUSTER_DIM_N * 3])
#ifdef AOS #ifdef AOS
#define POS_DATA_LAYOUT "AoS" #define POS_DATA_LAYOUT "AoS"
#define atom_x(i) atom->x[(i) * 3 + 0] #define atom_x(i) atom->x[(i) * 3 + 0]
#define atom_y(i) atom->x[(i) * 3 + 1] #define atom_y(i) atom->x[(i) * 3 + 1]
#define atom_z(i) atom->x[(i) * 3 + 2] #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 #else
#define POS_DATA_LAYOUT "SoA" #define POS_DATA_LAYOUT "SoA"
#define atom_x(i) atom->x[i] #define atom_x(i) atom->x[i]
#define atom_y(i) atom->y[i] #define atom_y(i) atom->y[i]
#define atom_z(i) atom->z[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
#endif #endif

3
gromacs/includes/neighbor.h
View File

@ -36,6 +36,7 @@ typedef struct {
extern void initNeighbor(Neighbor*, Parameter*); extern void initNeighbor(Neighbor*, Parameter*);
extern void setupNeighbor(Parameter*); extern void setupNeighbor(Parameter*);
extern void binatoms(Atom*); extern void binatoms(Atom*);
extern void buildNeighbor(Atom*, Neighbor*); extern void buildNeighbor(Parameter*, Atom*, Neighbor*);
extern void sortAtom(Atom*); extern void sortAtom(Atom*);
extern void buildClusters(Parameter*, Atom*);
#endif #endif

408
gromacs/neighbor.c
View File

@ -27,28 +27,34 @@
#include <neighbor.h> #include <neighbor.h>
#include <parameter.h> #include <parameter.h>
#include <atom.h> #include <atom.h>
#include <util.h>
#define SMALL 1.0e-6 #define SMALL 1.0e-6
#define FACTOR 0.999 #define FACTOR 0.999
static MD_FLOAT xprd, yprd, zprd; static MD_FLOAT xprd, yprd;
static MD_FLOAT bininvx, bininvy, bininvz; static MD_FLOAT bininvx, bininvy;
static int mbinxlo, mbinylo, mbinzlo; static int mbinxlo, mbinylo;
static int nbinx, nbiny, nbinz; static int nbinx, nbiny;
static int mbinx, mbiny, mbinz; // n bins in x, y, z static int mbinx, mbiny; // n bins in x, y
static int *bincount; static int *bincount;
static int *bins; static int *bins;
static int *cluster_bincount;
static int *cluster_bins;
static int mbins; //total number of bins static int mbins; //total number of bins
static int atoms_per_bin; // max atoms per bin 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 cutneigh;
static MD_FLOAT cutneighsq; // neighbor cutoff squared static MD_FLOAT cutneighsq; // neighbor cutoff squared
static int nmax; static int nmax;
static int nstencil; // # of bins in stencil static int nstencil; // # of bins in stencil
static int* stencil; // stencil list of bin offsets 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 int coord2bin(MD_FLOAT, MD_FLOAT);
static MD_FLOAT bindist(int, int, int); static MD_FLOAT bindist(int, int);
extern int *PBCx, *PBCy, *PBCz;
/* exported subroutines */ /* exported subroutines */
void initNeighbor(Neighbor *neighbor, Parameter *param) void initNeighbor(Neighbor *neighbor, Parameter *param)
@ -56,60 +62,52 @@ void initNeighbor(Neighbor *neighbor, Parameter *param)
MD_FLOAT neighscale = 5.0 / 6.0; MD_FLOAT neighscale = 5.0 / 6.0;
xprd = param->nx * param->lattice; xprd = param->nx * param->lattice;
yprd = param->ny * param->lattice; yprd = param->ny * param->lattice;
zprd = param->nz * param->lattice;
cutneigh = param->cutneigh; cutneigh = param->cutneigh;
nbinx = neighscale * param->nx; nbinx = neighscale * param->nx;
nbiny = neighscale * param->ny; nbiny = neighscale * param->ny;
nbinz = neighscale * param->nz;
nmax = 0; nmax = 0;
atoms_per_bin = 8; atoms_per_bin = 8;
clusters_per_bin = (atoms_per_bin / CLUSTER_DIM_N) + 4;
stencil = NULL; stencil = NULL;
bins = NULL; bins = NULL;
bincount = NULL; bincount = NULL;
cluster_bins = NULL;
cluster_bincount = NULL;
neighbor->maxneighs = 100; neighbor->maxneighs = 100;
neighbor->numneigh = NULL; neighbor->numneigh = NULL;
neighbor->neighbors = NULL; neighbor->neighbors = NULL;
} }
void setupNeighbor(Parameter* param) void setupNeighbor(Parameter* param) {
{
MD_FLOAT coord; MD_FLOAT coord;
int mbinxhi, mbinyhi, mbinzhi; int mbinxhi, mbinyhi;
int nextx, nexty, nextz; int nextx, nexty, nextz;
if(param->input_file != NULL) { if(param->input_file != NULL) {
xprd = param->xprd; xprd = param->xprd;
yprd = param->yprd; yprd = param->yprd;
zprd = param->zprd;
} }
// TODO: update lo and hi for standard case and use them here instead // TODO: update lo and hi for standard case and use them here instead
MD_FLOAT xlo = 0.0; MD_FLOAT xhi = xprd; MD_FLOAT xlo = 0.0; MD_FLOAT xhi = xprd;
MD_FLOAT ylo = 0.0; MD_FLOAT yhi = yprd; MD_FLOAT ylo = 0.0; MD_FLOAT yhi = yprd;
MD_FLOAT zlo = 0.0; MD_FLOAT zhi = zprd;
cutneighsq = cutneigh * cutneigh; cutneighsq = cutneigh * cutneigh;
if(param->input_file != NULL) { if(param->input_file != NULL) {
binsizex = cutneigh * 0.5; binsizex = cutneigh * 0.5;
binsizey = cutneigh * 0.5; binsizey = cutneigh * 0.5;
binsizez = cutneigh * 0.5;
nbinx = (int)((param->xhi - param->xlo) / binsizex); nbinx = (int)((param->xhi - param->xlo) / binsizex);
nbiny = (int)((param->yhi - param->ylo) / binsizey); nbiny = (int)((param->yhi - param->ylo) / binsizey);
nbinz = (int)((param->zhi - param->zlo) / binsizez);
if(nbinx == 0) { nbinx = 1; } if(nbinx == 0) { nbinx = 1; }
if(nbiny == 0) { nbiny = 1; } if(nbiny == 0) { nbiny = 1; }
if(nbinz == 0) { nbinz = 1; }
bininvx = nbinx / (param->xhi - param->xlo); bininvx = nbinx / (param->xhi - param->xlo);
bininvy = nbiny / (param->yhi - param->ylo); bininvy = nbiny / (param->yhi - param->ylo);
bininvz = nbinz / (param->zhi - param->zlo);
} else { } else {
binsizex = xprd / nbinx; binsizex = xprd / nbinx;
binsizey = yprd / nbiny; binsizey = yprd / nbiny;
binsizez = zprd / nbinz;
bininvx = 1.0 / binsizex; bininvx = 1.0 / binsizex;
bininvy = 1.0 / binsizey; bininvy = 1.0 / binsizey;
bininvz = 1.0 / binsizez;
} }
coord = xlo - cutneigh - SMALL * xprd; coord = xlo - cutneigh - SMALL * xprd;
@ -128,14 +126,6 @@ void setupNeighbor(Parameter* param)
coord = yhi + cutneigh + SMALL * yprd; coord = yhi + cutneigh + SMALL * yprd;
mbinyhi = (int) (coord * bininvy); 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; mbinxlo = mbinxlo - 1;
mbinxhi = mbinxhi + 1; mbinxhi = mbinxhi + 1;
mbinx = mbinxhi - mbinxlo + 1; mbinx = mbinxhi - mbinxlo + 1;
@ -144,55 +134,82 @@ void setupNeighbor(Parameter* param)
mbinyhi = mbinyhi + 1; mbinyhi = mbinyhi + 1;
mbiny = mbinyhi - mbinylo + 1; mbiny = mbinyhi - mbinylo + 1;
mbinzlo = mbinzlo - 1;
mbinzhi = mbinzhi + 1;
mbinz = mbinzhi - mbinzlo + 1;
nextx = (int) (cutneigh * bininvx); nextx = (int) (cutneigh * bininvx);
if(nextx * binsizex < FACTOR * cutneigh) nextx++; if(nextx * binsizex < FACTOR * cutneigh) nextx++;
nexty = (int) (cutneigh * bininvy); nexty = (int) (cutneigh * bininvy);
if(nexty * binsizey < FACTOR * cutneigh) nexty++; if(nexty * binsizey < FACTOR * cutneigh) nexty++;
nextz = (int) (cutneigh * bininvz);
if(nextz * binsizez < FACTOR * cutneigh) nextz++;
if (stencil) { if (stencil) {
free(stencil); free(stencil);
} }
stencil = (int*) malloc( stencil = (int *) malloc((2 * nexty + 1) * (2 * nextx + 1) * sizeof(int));
(2 * nextz + 1) * (2 * nexty + 1) * (2 * nextx + 1) * sizeof(int));
nstencil = 0; nstencil = 0;
int kstart = -nextz;
for(int k = kstart; k <= nextz; k++) {
for(int j = -nexty; j <= nexty; j++) { for(int j = -nexty; j <= nexty; j++) {
for(int i = -nextx; i <= nextx; i++) { for(int i = -nextx; i <= nextx; i++) {
if(bindist(i, j, k) < cutneighsq) { if(bindist(i, j) < cutneighsq) {
stencil[nstencil++] = stencil[nstencil++] = j * mbinx + i;
k * mbiny * mbinx + j * mbinx + i;
}
} }
} }
} }
mbins = mbinx * mbiny * mbinz; mbins = mbinx * mbiny;
if (bincount) { if (bincount) { free(bincount); }
free(bincount);
}
bincount = (int*) malloc(mbins * sizeof(int)); bincount = (int*) malloc(mbins * sizeof(int));
if (bins) { if (bins) { free(bins); }
free(bins);
}
bins = (int*) malloc(mbins * atoms_per_bin * sizeof(int)); 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; int nall = atom->Nlocal + atom->Nghost;
/* extend atom arrays if necessary */ /* extend atom arrays if necessary */
@ -206,6 +223,9 @@ void buildNeighbor(Atom *atom, Neighbor *neighbor)
/* bin local & ghost atoms */ /* bin local & ghost atoms */
binatoms(atom); binatoms(atom);
buildClusters(param, atom);
const MD_FLOAT rBB = cutneighsq / 2.0; // TODO: change this
int resize = 1; int resize = 1;
/* loop over each atom, storing neighbors */ /* loop over each atom, storing neighbors */
@ -213,45 +233,49 @@ void buildNeighbor(Atom *atom, Neighbor *neighbor)
int new_maxneighs = neighbor->maxneighs; int new_maxneighs = neighbor->maxneighs;
resize = 0; resize = 0;
for(int i = 0; i < atom->Nlocal; i++) { for(int ci = 0; ci < atom->Nclusters_local; ci++) {
int* neighptr = &(neighbor->neighbors[i * neighbor->maxneighs]); int* neighptr = &(neighbor->neighbors[ci * neighbor->maxneighs]);
int n = 0; int n = 0;
MD_FLOAT xtmp = atom_x(i); int ibin = atom->clusters[ci].bin;
MD_FLOAT ytmp = atom_y(i); MD_FLOAT ibb_zmin = atom->clusters[ci].bbminz;
MD_FLOAT ztmp = atom_z(i); MD_FLOAT ibb_zmax = atom->clusters[ci].bbmaxz;
int ibin = coord2bin(xtmp, ytmp, ztmp);
#ifdef EXPLICIT_TYPES
int type_i = atom->type[i];
#endif
for(int k = 0; k < nstencil; k++) { for(int k = 0; k < nstencil; k++) {
int jbin = ibin + stencil[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++) { do {
int j = loc_bin[m]; 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 ){ while(m < cluster_bincount[jbin] && jbb_zmax - ibb_zmax < cutneighsq) {
continue; if((jbb_zmin - ibb_zmax) * (jbb_zmin - ibb_zmax) > cutneighsq) {
break;
} }
MD_FLOAT delx = xtmp - atom_x(j); double d_bb_sq = getBoundingBoxDistanceSq(atom, ci, cj);
MD_FLOAT dely = ytmp - atom_y(j); if(d_bb_sq < cutneighsq) {
MD_FLOAT delz = ztmp - atom_z(j); if(d_bb_sq < rBB || atomDistanceInRange(atom, ci, cj, cutneighsq)) {
MD_FLOAT rsq = delx * delx + dely * dely + delz * delz; if(cj == ci) {
// Add to neighbor list with diagonal interaction mask
#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;
} }
neighptr[n++] = cj;
} }
} }
neighbor->numneigh[i] = n; m++;
cj = loc_bin[m];
jbb_zmin = atom->clusters[cj].bbminz;
jbb_zmax = atom->clusters[cj].bbmaxz;
}
}
neighbor->numneigh[ci] = n;
if(n >= neighbor->maxneighs) { if(n >= neighbor->maxneighs) {
resize = 1; resize = 1;
@ -272,8 +296,7 @@ void buildNeighbor(Atom *atom, Neighbor *neighbor)
} }
/* internal subroutines */ /* internal subroutines */
MD_FLOAT bindist(int i, int j, int k) MD_FLOAT bindist(int i, int j) {
{
MD_FLOAT delx, dely, delz; MD_FLOAT delx, dely, delz;
if(i > 0) { if(i > 0) {
@ -292,20 +315,11 @@ MD_FLOAT bindist(int i, int j, int k)
dely = (j + 1) * binsizey; dely = (j + 1) * binsizey;
} }
if(k > 0) { return (delx * delx + dely * dely);
delz = (k - 1) * binsizez;
} else if(k == 0) {
delz = 0.0;
} else {
delz = (k + 1) * binsizez;
} }
return (delx * delx + dely * dely + delz * delz); int coord2bin(MD_FLOAT xin, MD_FLOAT yin) {
} int ix, iy;
int coord2bin(MD_FLOAT xin, MD_FLOAT yin, MD_FLOAT zin)
{
int ix, iy, iz;
if(xin >= xprd) { if(xin >= xprd) {
ix = (int)((xin - xprd) * bininvx) + nbinx - mbinxlo; 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; iy = (int)(yin * bininvy) - mbinylo - 1;
} }
if(zin >= zprd) { return (iy * mbinx + ix + 1);
iz = (int)((zin - zprd) * bininvz) + nbinz - mbinzlo; }
} else if(zin >= 0.0) {
iz = (int)(zin * bininvz) - mbinzlo; 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 { } else {
iz = (int)(zin * bininvz) - mbinzlo - 1; *ix = (int)(xin * bininvx) - mbinxlo - 1;
} }
return (iz * mbiny * mbinx + iy * mbinx + ix + 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) void binatoms(Atom *atom) {
{
int nall = atom->Nlocal + atom->Nghost; int nall = atom->Nlocal + atom->Nghost;
int resize = 1; int resize = 1;
@ -347,7 +370,7 @@ void binatoms(Atom *atom)
} }
for(int i = 0; i < nall; i++) { 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) { if(bincount[ibin] < atoms_per_bin) {
int ac = bincount[ibin]++; int ac = bincount[ibin]++;
@ -365,56 +388,137 @@ void binatoms(Atom *atom)
} }
} }
void sortAtom(Atom* atom) { // TODO: Use pigeonhole sorting
binatoms(atom); void sortBinAtomsByZCoord(Parameter *param, Atom *atom, int bin) {
int Nmax = atom->Nmax; int c = bincount[bin];
int* binpos = bincount; int *bin_ptr = &bins[bin * atoms_per_bin];
for(int i=1; i<mbins; i++) { for(int ac_i = 0; ac_i < c; ac_i++) {
binpos[i] += binpos[i-1]; int i = bin_ptr[ac_i];
} int min_ac = ac_i;
int min_idx = i;
MD_FLOAT min_z = atom_z(i);
#ifdef AOS for(int ac_j = ac_i + 1; ac_j < c; ac_j++) {
double* new_x = (double*) malloc(Nmax * sizeof(MD_FLOAT) * 3); int j = bin_ptr[ac_j];
#else MD_FLOAT zj = atom_z(j);
double* new_x = (double*) malloc(Nmax * sizeof(MD_FLOAT)); if(zj < min_z) {
double* new_y = (double*) malloc(Nmax * sizeof(MD_FLOAT)); min_ac = zj;
double* new_z = (double*) malloc(Nmax * sizeof(MD_FLOAT)); min_idx = j;
#endif min_z = zj;
double* new_vx = (double*) malloc(Nmax * sizeof(MD_FLOAT));
double* new_vy = (double*) malloc(Nmax * sizeof(MD_FLOAT));
double* new_vz = (double*) malloc(Nmax * sizeof(MD_FLOAT));
double* old_x = atom->x; 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;
for(int mybin = 0; mybin<mbins; mybin++) {
int start = mybin>0?binpos[mybin-1]:0;
int count = binpos[mybin] - start;
for(int k=0; k<count; k++) {
int new_i = start + k;
int old_i = bins[mybin * atoms_per_bin + k];
#ifdef AOS
new_x[new_i * 3 + 0] = old_x[old_i * 3 + 0];
new_x[new_i * 3 + 1] = old_x[old_i * 3 + 1];
new_x[new_i * 3 + 2] = old_x[old_i * 3 + 2];
#else
new_x[new_i] = old_x[old_i];
new_y[new_i] = old_y[old_i];
new_z[new_i] = old_z[old_i];
#endif
new_vx[new_i] = old_vx[old_i];
new_vy[new_i] = old_vy[old_i];
new_vz[new_i] = old_vz[old_i];
} }
} }
free(atom->x); bin_ptr[ac_i] = min_idx;
atom->x = new_x; bin_ptr[min_ac] = i;
#ifndef AOS }
free(atom->y); }
free(atom->z);
atom->y = new_y; atom->z = new_z; void buildClusters(Parameter *param, Atom *atom) {
#endif for(int bin = 0; bin < mbins; bin++) {
free(atom->vx); free(atom->vy); free(atom->vz); sortBinAtomsByZCoord(param, atom, bin);
atom->vx = new_vx; atom->vy = new_vy; atom->vz = new_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);
if(nclusters < clusters_per_bin) {
for(int cl = 0; cl < nclusters; cl++) {
const int ci = atom->Nclusters_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));
}
}
}
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]++;
}
} }

117
gromacs/pbc.c
View File

@ -35,8 +35,7 @@ static int *PBCx, *PBCy, *PBCz;
static void growPbc(Atom*); static void growPbc(Atom*);
/* exported subroutines */ /* exported subroutines */
void initPbc(Atom* atom) void initPbc(Atom* atom) {
{
NmaxGhost = 0; NmaxGhost = 0;
atom->border_map = NULL; atom->border_map = NULL;
PBCx = NULL; PBCy = NULL; PBCz = NULL; PBCx = NULL; PBCy = NULL; PBCz = NULL;
@ -44,36 +43,37 @@ void initPbc(Atom* atom)
/* update coordinates of ghost atoms */ /* update coordinates of ghost atoms */
/* uses mapping created in setupPbc */ /* uses mapping created in setupPbc */
void updatePbc(Atom *atom, Parameter *param) void updatePbc(Atom *atom, Parameter *param) {
{
int *border_map = atom->border_map; int *border_map = atom->border_map;
int nlocal = atom->Nlocal; int nlocal = atom->Nclusters_local;
MD_FLOAT xprd = param->xprd; MD_FLOAT xprd = param->xprd;
MD_FLOAT yprd = param->yprd; MD_FLOAT yprd = param->yprd;
MD_FLOAT zprd = param->zprd; MD_FLOAT zprd = param->zprd;
for(int i = 0; i < atom->Nghost; i++) { for(int ci = 0; ci < atom->Nclusters_ghost; ci++) {
atom_x(nlocal + i) = atom_x(border_map[i]) + PBCx[i] * xprd; MD_FLOAT *cptr = cluster_ptr(nlocal + ci);
atom_y(nlocal + i) = atom_y(border_map[i]) + PBCy[i] * yprd; MD_FLOAT *bmap_cptr = cluster_ptr(border_map[ci]);
atom_z(nlocal + i) = atom_z(border_map[i]) + PBCz[i] * zprd;
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 /* relocate atoms that have left domain according
* to periodic boundary conditions */ * to periodic boundary conditions */
void updateAtomsPbc(Atom *atom, Parameter *param) void updateAtomsPbc(Atom *atom, Parameter *param) {
{
MD_FLOAT xprd = param->xprd; MD_FLOAT xprd = param->xprd;
MD_FLOAT yprd = param->yprd; MD_FLOAT yprd = param->yprd;
MD_FLOAT zprd = param->zprd; MD_FLOAT zprd = param->zprd;
for(int i = 0; i < atom->Nlocal; i++) { for(int i = 0; i < atom->Nlocal; i++) {
if(atom_x(i) < 0.0) { if(atom_x(i) < 0.0) {
atom_x(i) += xprd; atom_x(i) += xprd;
} else if(atom_x(i) >= xprd) { } else if(atom_x(i) >= xprd) {
atom_x(i) -= xprd; atom_x(i) -= xprd;
}
if(atom_y(i) < 0.0) { if(atom_y(i) < 0.0) {
atom_y(i) += yprd; atom_y(i) += yprd;
@ -95,14 +95,19 @@ void updateAtomsPbc(Atom *atom, Parameter *param)
* that are then enforced in updatePbc */ * that are then enforced in updatePbc */
#define ADDGHOST(dx,dy,dz) \ #define ADDGHOST(dx,dy,dz) \
Nghost++; \ Nghost++; \
border_map[Nghost] = i; \ border_map[Nghost] = ci; \
PBCx[Nghost] = dx; \ PBCx[Nghost] = dx; \
PBCy[Nghost] = dy; \ PBCy[Nghost] = dy; \
PBCz[Nghost] = dz; \ PBCz[Nghost] = dz; \
atom->type[atom->Nlocal + Nghost] = atom->type[i] 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; int *border_map = atom->border_map;
MD_FLOAT xprd = param->xprd; MD_FLOAT xprd = param->xprd;
MD_FLOAT yprd = param->yprd; MD_FLOAT yprd = param->yprd;
@ -110,58 +115,62 @@ void setupPbc(Atom *atom, Parameter *param)
MD_FLOAT Cutneigh = param->cutneigh; MD_FLOAT Cutneigh = param->cutneigh;
int Nghost = -1; int Nghost = -1;
for(int i = 0; i < atom->Nlocal; i++) { for(int ci = 0; ci < atom->Nclusters_local; ci++) {
if (atom->Nclusters_local + atom->Nclusters_ghost + 7 >= atom->Nclusters_max) {
if (atom->Nlocal + Nghost + 7 >= atom->Nmax) { growClusters(atom);
growAtom(atom);
} }
if (Nghost + 7 >= NmaxGhost) {
if (atom->Nclusters_ghost + 7 >= NmaxGhost) {
growPbc(atom); growPbc(atom);
border_map = atom->border_map; border_map = atom->border_map;
} }
MD_FLOAT x = atom_x(i); MD_FLOAT bbminx = atom->clusters[ci].bbminx;
MD_FLOAT y = atom_y(i); MD_FLOAT bbmaxx = atom->clusters[ci].bbmaxx;
MD_FLOAT z = atom_z(i); 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 */ /* Setup ghost atoms */
/* 6 planes */ /* 6 planes */
if (x < Cutneigh) { ADDGHOST(+1,0,0); } if (bbminx < Cutneigh) { ADDGHOST(+1,0,0); }
if (x >= (xprd-Cutneigh)) { ADDGHOST(-1,0,0); } if (bbmaxx >= (xprd-Cutneigh)) { ADDGHOST(-1,0,0); }
if (y < Cutneigh) { ADDGHOST(0,+1,0); } if (bbminy < Cutneigh) { ADDGHOST(0,+1,0); }
if (y >= (yprd-Cutneigh)) { ADDGHOST(0,-1,0); } if (bbmaxy >= (yprd-Cutneigh)) { ADDGHOST(0,-1,0); }
if (z < Cutneigh) { ADDGHOST(0,0,+1); } if (bbminz < Cutneigh) { ADDGHOST(0,0,+1); }
if (z >= (zprd-Cutneigh)) { ADDGHOST(0,0,-1); } if (bbmaxz >= (zprd-Cutneigh)) { ADDGHOST(0,0,-1); }
/* 8 corners */ /* 8 corners */
if (x < Cutneigh && y < Cutneigh && z < Cutneigh) { ADDGHOST(+1,+1,+1); } if (bbminx < Cutneigh && bbminy < Cutneigh && bbminz < Cutneigh) { ADDGHOST(+1,+1,+1); }
if (x < Cutneigh && y >= (yprd-Cutneigh) && z < Cutneigh) { ADDGHOST(+1,-1,+1); } if (bbminx < Cutneigh && bbmaxy >= (yprd-Cutneigh) && bbminz < Cutneigh) { ADDGHOST(+1,-1,+1); }
if (x < Cutneigh && y >= Cutneigh && z >= (zprd-Cutneigh)) { ADDGHOST(+1,+1,-1); } if (bbminx < Cutneigh && bbminy < Cutneigh && bbmaxz >= (zprd-Cutneigh)) { ADDGHOST(+1,+1,-1); }
if (x < Cutneigh && y >= (yprd-Cutneigh) && z >= (zprd-Cutneigh)) { ADDGHOST(+1,-1,-1); } if (bbminx < Cutneigh && bbmaxy >= (yprd-Cutneigh) && bbmaxz >= (zprd-Cutneigh)) { ADDGHOST(+1,-1,-1); }
if (x >= (xprd-Cutneigh) && y < Cutneigh && z < Cutneigh) { ADDGHOST(-1,+1,+1); } if (bbmaxx >= (xprd-Cutneigh) && bbminy < Cutneigh && bbminz < Cutneigh) { ADDGHOST(-1,+1,+1); }
if (x >= (xprd-Cutneigh) && y >= (yprd-Cutneigh) && z < Cutneigh) { ADDGHOST(-1,-1,+1); } if (bbmaxx >= (xprd-Cutneigh) && bbmaxy >= (yprd-Cutneigh) && bbminz < Cutneigh) { ADDGHOST(-1,-1,+1); }
if (x >= (xprd-Cutneigh) && y < Cutneigh && z >= (zprd-Cutneigh)) { ADDGHOST(-1,+1,-1); } if (bbmaxx >= (xprd-Cutneigh) && bbminy < Cutneigh && bbmaxz >= (zprd-Cutneigh)) { ADDGHOST(-1,+1,-1); }
if (x >= (xprd-Cutneigh) && y >= (yprd-Cutneigh) && z >= (zprd-Cutneigh)) { ADDGHOST(-1,-1,-1); } if (bbmaxx >= (xprd-Cutneigh) && bbmaxy >= (yprd-Cutneigh) && bbmaxz >= (zprd-Cutneigh)) { ADDGHOST(-1,-1,-1); }
/* 12 edges */ /* 12 edges */
if (x < Cutneigh && z < Cutneigh) { ADDGHOST(+1,0,+1); } if (bbminx < Cutneigh && bbminz < Cutneigh) { ADDGHOST(+1,0,+1); }
if (x < Cutneigh && z >= (zprd-Cutneigh)) { ADDGHOST(+1,0,-1); } if (bbminx < Cutneigh && bbmaxz >= (zprd-Cutneigh)) { ADDGHOST(+1,0,-1); }
if (x >= (xprd-Cutneigh) && z < Cutneigh) { ADDGHOST(-1,0,+1); } if (bbmaxx >= (xprd-Cutneigh) && bbminz < Cutneigh) { ADDGHOST(-1,0,+1); }
if (x >= (xprd-Cutneigh) && z >= (zprd-Cutneigh)) { ADDGHOST(-1,0,-1); } if (bbmaxx >= (xprd-Cutneigh) && bbmaxz >= (zprd-Cutneigh)) { ADDGHOST(-1,0,-1); }
if (y < Cutneigh && z < Cutneigh) { ADDGHOST(0,+1,+1); } if (bbminy < Cutneigh && bbminz < Cutneigh) { ADDGHOST(0,+1,+1); }
if (y < Cutneigh && z >= (zprd-Cutneigh)) { ADDGHOST(0,+1,-1); } if (bbminy < Cutneigh && bbmaxz >= (zprd-Cutneigh)) { ADDGHOST(0,+1,-1); }
if (y >= (yprd-Cutneigh) && z < Cutneigh) { ADDGHOST(0,-1,+1); } if (bbmaxy >= (yprd-Cutneigh) && bbminz < Cutneigh) { ADDGHOST(0,-1,+1); }
if (y >= (yprd-Cutneigh) && z >= (zprd-Cutneigh)) { ADDGHOST(0,-1,-1); } if (bbmaxy >= (yprd-Cutneigh) && bbmaxz >= (zprd-Cutneigh)) { ADDGHOST(0,-1,-1); }
if (y < Cutneigh && x < Cutneigh) { ADDGHOST(+1,+1,0); } if (bbminy < Cutneigh && bbminx < Cutneigh) { ADDGHOST(+1,+1,0); }
if (y < Cutneigh && x >= (xprd-Cutneigh)) { ADDGHOST(-1,+1,0); } if (bbminy < Cutneigh && bbmaxx >= (xprd-Cutneigh)) { ADDGHOST(-1,+1,0); }
if (y >= (yprd-Cutneigh) && x < Cutneigh) { ADDGHOST(+1,-1,0); } if (bbmaxy >= (yprd-Cutneigh) && bbminx < Cutneigh) { ADDGHOST(+1,-1,0); }
if (y >= (yprd-Cutneigh) && x >= (xprd-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 // 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 */ /* internal subroutines */
void growPbc(Atom* atom) void growPbc(Atom* atom) {
{
int nold = NmaxGhost; int nold = NmaxGhost;
NmaxGhost += DELTA; NmaxGhost += DELTA;