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

Add first working version with 4x8 config (ref kernel)

Browse Source 
Signed-off-by: Rafael Ravedutti <rafaelravedutti@gmail.com>
This commit is contained in:
Rafael Ravedutti 2022-03-10 22:33:41 +01:00
parent 22d0f0b958
commit d79c3c2a1d
7 changed files with 191 additions and 187 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
gromacs/atom.c
View File

@ -429,10 +429,10 @@ void growAtom(Atom *atom) {
void growClusters(Atom *atom) { void growClusters(Atom *atom) {
int nold = atom->Nclusters_max; int nold = atom->Nclusters_max;
int jfac = MIN(1, CLUSTER_M / CLUSTER_N); // If M>N, we need to allocate more j-clusters int jterm = MAX(1, CLUSTER_M / CLUSTER_N); // If M>N, we need to allocate more j-clusters
atom->Nclusters_max += DELTA; atom->Nclusters_max += DELTA;
atom->iclusters = (Cluster*) reallocate(atom->iclusters, ALIGNMENT, atom->Nclusters_max * sizeof(Cluster), nold * sizeof(Cluster)); atom->iclusters = (Cluster*) reallocate(atom->iclusters, ALIGNMENT, atom->Nclusters_max * sizeof(Cluster), nold * sizeof(Cluster));
atom->jclusters = (Cluster*) reallocate(atom->jclusters, ALIGNMENT, atom->Nclusters_max * jfac * sizeof(Cluster), nold * jfac * sizeof(Cluster)); atom->jclusters = (Cluster*) reallocate(atom->jclusters, ALIGNMENT, atom->Nclusters_max * jterm * sizeof(Cluster), nold * jterm * sizeof(Cluster));
atom->icluster_bin = (int*) reallocate(atom->icluster_bin, ALIGNMENT, atom->Nclusters_max * sizeof(int), nold * sizeof(int)); atom->icluster_bin = (int*) reallocate(atom->icluster_bin, ALIGNMENT, atom->Nclusters_max * sizeof(int), nold * sizeof(int));
atom->cl_x = (MD_FLOAT*) reallocate(atom->cl_x, ALIGNMENT, atom->Nclusters_max * CLUSTER_M * 3 * sizeof(MD_FLOAT), nold * CLUSTER_M * 3 * sizeof(MD_FLOAT)); atom->cl_x = (MD_FLOAT*) reallocate(atom->cl_x, ALIGNMENT, atom->Nclusters_max * CLUSTER_M * 3 * sizeof(MD_FLOAT), nold * CLUSTER_M * 3 * sizeof(MD_FLOAT));
atom->cl_f = (MD_FLOAT*) reallocate(atom->cl_f, ALIGNMENT, atom->Nclusters_max * CLUSTER_M * 3 * sizeof(MD_FLOAT), nold * CLUSTER_M * 3 * sizeof(MD_FLOAT)); atom->cl_f = (MD_FLOAT*) reallocate(atom->cl_f, ALIGNMENT, atom->Nclusters_max * CLUSTER_M * 3 * sizeof(MD_FLOAT), nold * CLUSTER_M * 3 * sizeof(MD_FLOAT));

10
gromacs/force_lj.c
View File

@ -55,6 +55,8 @@ double computeForceLJ_ref(Parameter *param, Atom *atom, Neighbor *neighbor, Stat
#pragma omp parallel for #pragma omp parallel for
for(int ci = 0; ci < atom->Nclusters_local; ci++) { for(int ci = 0; ci < atom->Nclusters_local; ci++) {
int ci_cj0 = CJ0_FROM_CI(ci);
int ci_cj1 = CJ1_FROM_CI(ci);
int ci_vec_base = CI_VECTOR_BASE_INDEX(ci); int ci_vec_base = CI_VECTOR_BASE_INDEX(ci);
MD_FLOAT *ci_x = &atom->cl_x[ci_vec_base]; MD_FLOAT *ci_x = &atom->cl_x[ci_vec_base];
MD_FLOAT *ci_f = &atom->cl_f[ci_vec_base]; MD_FLOAT *ci_f = &atom->cl_f[ci_vec_base];
@ -76,7 +78,13 @@ double computeForceLJ_ref(Parameter *param, Atom *atom, Neighbor *neighbor, Stat
MD_FLOAT fiz = 0; MD_FLOAT fiz = 0;
for(int cjj = 0; cjj < CLUSTER_N; cjj++) { for(int cjj = 0; cjj < CLUSTER_N; cjj++) {
if(ci != cj || cii != cjj) { #if CLUSTER_M == CLUSTER_N
if(ci_cj0 != cj || cii != cjj) {
#elif CLUSTER_M < CLUSTER_N
if(ci_cj0 != cj || cii + CLUSTER_M * (ci & 0x1) != cjj) {
#else
if((ci_cj0 != cj || cii != cjj) && (ci_cj1 != cj || cii != cjj + CLUSTER_N)) {
#endif
MD_FLOAT delx = xtmp - cj_x[CL_X_OFFSET + cjj]; MD_FLOAT delx = xtmp - cj_x[CL_X_OFFSET + cjj];
MD_FLOAT dely = ytmp - cj_x[CL_Y_OFFSET + cjj]; MD_FLOAT dely = ytmp - cj_x[CL_Y_OFFSET + cjj];
MD_FLOAT delz = ztmp - cj_x[CL_Z_OFFSET + cjj]; MD_FLOAT delz = ztmp - cj_x[CL_Z_OFFSET + cjj];

1
gromacs/includes/atom.h
View File

@ -117,6 +117,7 @@ typedef struct {
int *cl_type; int *cl_type;
Cluster *iclusters, *jclusters; Cluster *iclusters, *jclusters;
int *icluster_bin; int *icluster_bin;
int dummy_cj;
} Atom; } Atom;
extern void initAtom(Atom*); extern void initAtom(Atom*);

1
gromacs/includes/neighbor.h
View File

@ -40,6 +40,7 @@ extern void buildNeighbor(Atom*, Neighbor*);
extern void pruneNeighbor(Parameter*, Atom*, Neighbor*); extern void pruneNeighbor(Parameter*, Atom*, Neighbor*);
extern void sortAtom(Atom*); extern void sortAtom(Atom*);
extern void buildClusters(Atom*); extern void buildClusters(Atom*);
extern void defineJClusters(Atom*);
extern void binClusters(Atom*); extern void binClusters(Atom*);
extern void updateSingleAtoms(Atom*); extern void updateSingleAtoms(Atom*);
#endif #endif

3
gromacs/main.c
View File

@ -70,6 +70,7 @@ double setup(Parameter *param, Eam *eam, Atom *atom, Neighbor *neighbor, Stats *
setupThermo(param, atom->Natoms); setupThermo(param, atom->Natoms);
if(param->input_file == NULL) { adjustThermo(param, atom); } if(param->input_file == NULL) { adjustThermo(param, atom); }
buildClusters(atom); buildClusters(atom);
defineJClusters(atom);
setupPbc(atom, param); setupPbc(atom, param);
binClusters(atom); binClusters(atom);
buildNeighbor(atom, neighbor); buildNeighbor(atom, neighbor);
@ -84,6 +85,7 @@ double reneighbour(Parameter *param, Atom *atom, Neighbor *neighbor) {
updateSingleAtoms(atom); updateSingleAtoms(atom);
updateAtomsPbc(atom, param); updateAtomsPbc(atom, param);
buildClusters(atom); buildClusters(atom);
defineJClusters(atom);
setupPbc(atom, param); setupPbc(atom, param);
binClusters(atom); binClusters(atom);
buildNeighbor(atom, neighbor); buildNeighbor(atom, neighbor);
@ -313,6 +315,7 @@ int main(int argc, char** argv) {
} }
timer[TOTAL] = getTimeStamp() - timer[TOTAL]; timer[TOTAL] = getTimeStamp() - timer[TOTAL];
updateSingleAtoms(&atom);
computeThermo(-1, &param, &atom); computeThermo(-1, &param, &atom);
if(param.xtc_file != NULL) { if(param.xtc_file != NULL) {

184
gromacs/neighbor.c
View File

@ -201,11 +201,10 @@ int atomDistanceInRange(Atom *atom, int ci, int cj, MD_FLOAT rsq) {
void buildNeighbor(Atom *atom, Neighbor *neighbor) { void buildNeighbor(Atom *atom, Neighbor *neighbor) {
DEBUG_MESSAGE("buildNeighbor start\n"); DEBUG_MESSAGE("buildNeighbor start\n");
int nall = atom->Nclusters_local + atom->Nclusters_ghost;
/* extend atom arrays if necessary */ /* extend atom arrays if necessary */
if(nall > nmax) { if(atom->Nclusters_local > nmax) {
nmax = nall; nmax = atom->Nclusters_local;
if(neighbor->numneigh) free(neighbor->numneigh); if(neighbor->numneigh) free(neighbor->numneigh);
if(neighbor->neighbors) free(neighbor->neighbors); if(neighbor->neighbors) free(neighbor->neighbors);
neighbor->numneigh = (int*) malloc(nmax * sizeof(int)); neighbor->numneigh = (int*) malloc(nmax * sizeof(int));
@ -307,7 +306,7 @@ void buildNeighbor(Atom *atom, Neighbor *neighbor) {
// Fill neighbor list with dummy values to fit vector width // Fill neighbor list with dummy values to fit vector width
if(CLUSTER_N < VECTOR_WIDTH) { if(CLUSTER_N < VECTOR_WIDTH) {
while(n % (VECTOR_WIDTH / CLUSTER_N)) { while(n % (VECTOR_WIDTH / CLUSTER_N)) {
neighptr[n++] = CJ0_FROM_CI(nall - 1); // Last cluster is always a dummy cluster neighptr[n++] = atom->dummy_cj; // Last cluster is always a dummy cluster
} }
} }
@ -332,38 +331,40 @@ void buildNeighbor(Atom *atom, Neighbor *neighbor) {
/* /*
DEBUG_MESSAGE("\ncutneighsq = %f, rbb_sq = %f\n", cutneighsq, rbb_sq); DEBUG_MESSAGE("\ncutneighsq = %f, rbb_sq = %f\n", cutneighsq, rbb_sq);
for(int ci = 0; ci < 6; ci++) { for(int ci = 0; ci < 6; ci++) {
MD_FLOAT *ciptr = cluster_pos_ptr(ci); int ci_vec_base = CI_VECTOR_BASE_INDEX(ci);
MD_FLOAT *ci_x = &atom->cl_x[ci_vec_base];
int* neighptr = &(neighbor->neighbors[ci * neighbor->maxneighs]); int* neighptr = &(neighbor->neighbors[ci * neighbor->maxneighs]);
DEBUG_MESSAGE("Cluster %d, bbx = {%f, %f}, bby = {%f, %f}, bbz = {%f, %f}\n", DEBUG_MESSAGE("Cluster %d, bbx = {%f, %f}, bby = {%f, %f}, bbz = {%f, %f}\n",
ci, ci,
atom->clusters[ci].bbminx, atom->iclusters[ci].bbminx,
atom->clusters[ci].bbmaxx, atom->iclusters[ci].bbmaxx,
atom->clusters[ci].bbminy, atom->iclusters[ci].bbminy,
atom->clusters[ci].bbmaxy, atom->iclusters[ci].bbmaxy,
atom->clusters[ci].bbminz, atom->iclusters[ci].bbminz,
atom->clusters[ci].bbmaxz); atom->iclusters[ci].bbmaxz);
for(int cii = 0; cii < CLUSTER_M; cii++) { for(int cii = 0; cii < CLUSTER_M; cii++) {
DEBUG_MESSAGE("%f, %f, %f\n", cluster_x(ciptr, cii), cluster_y(ciptr, cii), cluster_z(ciptr, cii)); DEBUG_MESSAGE("%f, %f, %f\n", ci_x[CL_X_OFFSET + cii], ci_x[CL_Y_OFFSET + cii], ci_x[CL_Z_OFFSET + cii]);
} }
DEBUG_MESSAGE("Neighbors:\n"); DEBUG_MESSAGE("Neighbors:\n");
for(int k = 0; k < neighbor->numneigh[ci]; k++) { for(int k = 0; k < neighbor->numneigh[ci]; k++) {
const int cj = neighptr[k]; int cj = neighptr[k];
MD_FLOAT *cjptr = cluster_pos_ptr(cj); int cj_vec_base = CJ_VECTOR_BASE_INDEX(cj);
MD_FLOAT *cj_x = &atom->cl_x[cj_vec_base];
DEBUG_MESSAGE(" Cluster %d, bbx = {%f, %f}, bby = {%f, %f}, bbz = {%f, %f}\n", DEBUG_MESSAGE(" Cluster %d, bbx = {%f, %f}, bby = {%f, %f}, bbz = {%f, %f}\n",
cj, cj,
atom->clusters[cj].bbminx, atom->jclusters[cj].bbminx,
atom->clusters[cj].bbmaxx, atom->jclusters[cj].bbmaxx,
atom->clusters[cj].bbminy, atom->jclusters[cj].bbminy,
atom->clusters[cj].bbmaxy, atom->jclusters[cj].bbmaxy,
atom->clusters[cj].bbminz, atom->jclusters[cj].bbminz,
atom->clusters[cj].bbmaxz); atom->jclusters[cj].bbmaxz);
for(int cjj = 0; cjj < CLUSTER_M; cjj++) { for(int cjj = 0; cjj < CLUSTER_N; cjj++) {
DEBUG_MESSAGE(" %f, %f, %f\n", cluster_x(cjptr, cjj), cluster_y(cjptr, cjj), cluster_z(cjptr, cjj)); DEBUG_MESSAGE(" %f, %f, %f\n", cj_x[CL_X_OFFSET + cjj], cj_x[CL_Y_OFFSET + cjj], cj_x[CL_Z_OFFSET + cjj]);
} }
} }
} }
@ -374,7 +375,6 @@ void buildNeighbor(Atom *atom, Neighbor *neighbor) {
void pruneNeighbor(Parameter *param, Atom *atom, Neighbor *neighbor) { void pruneNeighbor(Parameter *param, Atom *atom, Neighbor *neighbor) {
DEBUG_MESSAGE("pruneNeighbor start\n"); DEBUG_MESSAGE("pruneNeighbor start\n");
int nall = atom->Nclusters_local + atom->Nclusters_ghost;
//MD_FLOAT cutsq = param->cutforce * param->cutforce; //MD_FLOAT cutsq = param->cutforce * param->cutforce;
MD_FLOAT cutsq = cutneighsq; MD_FLOAT cutsq = cutneighsq;
@ -385,7 +385,7 @@ void pruneNeighbor(Parameter *param, Atom *atom, Neighbor *neighbor) {
// Remove dummy clusters if necessary // Remove dummy clusters if necessary
if(CLUSTER_N < VECTOR_WIDTH) { if(CLUSTER_N < VECTOR_WIDTH) {
while(neighs[numneighs - 1] == CJ0_FROM_CI(nall - 1)) { while(neighs[numneighs - 1] == atom->dummy_cj) {
numneighs--; numneighs--;
} }
} }
@ -403,7 +403,7 @@ void pruneNeighbor(Parameter *param, Atom *atom, Neighbor *neighbor) {
// Readd dummy clusters if necessary // Readd dummy clusters if necessary
if(CLUSTER_N < VECTOR_WIDTH) { if(CLUSTER_N < VECTOR_WIDTH) {
while(numneighs % (VECTOR_WIDTH / CLUSTER_N)) { while(numneighs % (VECTOR_WIDTH / CLUSTER_N)) {
neighs[numneighs++] = CJ0_FROM_CI(nall - 1); // Last cluster is always a dummy cluster neighs[numneighs++] = atom->dummy_cj; // Last cluster is always a dummy cluster
} }
} }
@ -617,7 +617,7 @@ void buildClusters(Atom *atom) {
void defineJClusters(Atom *atom) { void defineJClusters(Atom *atom) {
DEBUG_MESSAGE("defineJClusters start\n"); DEBUG_MESSAGE("defineJClusters start\n");
for(int ci = 0; ci < atom->Nclusters_local + atom->Nclusters_ghost; ci++) { for(int ci = 0; ci < atom->Nclusters_local; ci++) {
int cj0 = CJ0_FROM_CI(ci); int cj0 = CJ0_FROM_CI(ci);
if(CLUSTER_M == CLUSTER_N) { if(CLUSTER_M == CLUSTER_N) {
@ -687,13 +687,14 @@ void defineJClusters(Atom *atom) {
} else { } else {
if(ci % 2 == 0) { if(ci % 2 == 0) {
atom->jclusters[cj0].bbminx = MIN(atom->iclusters[ci].bbminx, atom->iclusters[ci + 1].bbminx); const int ci1 = ci + 1;
atom->jclusters[cj0].bbmaxx = MAX(atom->iclusters[ci].bbmaxx, atom->iclusters[ci + 1].bbmaxx); atom->jclusters[cj0].bbminx = MIN(atom->iclusters[ci].bbminx, atom->iclusters[ci1].bbminx);
atom->jclusters[cj0].bbminy = MIN(atom->iclusters[ci].bbminy, atom->iclusters[ci + 1].bbminy); atom->jclusters[cj0].bbmaxx = MAX(atom->iclusters[ci].bbmaxx, atom->iclusters[ci1].bbmaxx);
atom->jclusters[cj0].bbmaxy = MAX(atom->iclusters[ci].bbmaxy, atom->iclusters[ci + 1].bbmaxy); atom->jclusters[cj0].bbminy = MIN(atom->iclusters[ci].bbminy, atom->iclusters[ci1].bbminy);
atom->jclusters[cj0].bbminz = MIN(atom->iclusters[ci].bbminz, atom->iclusters[ci + 1].bbminz); atom->jclusters[cj0].bbmaxy = MAX(atom->iclusters[ci].bbmaxy, atom->iclusters[ci1].bbmaxy);
atom->jclusters[cj0].bbmaxz = MAX(atom->iclusters[ci].bbmaxz, atom->iclusters[ci + 1].bbmaxz); atom->jclusters[cj0].bbminz = MIN(atom->iclusters[ci].bbminz, atom->iclusters[ci1].bbminz);
atom->jclusters[cj0].natoms = atom->iclusters[ci].natoms + atom->iclusters[ci + 1].natoms; atom->jclusters[cj0].bbmaxz = MAX(atom->iclusters[ci].bbmaxz, atom->iclusters[ci1].bbmaxz);
atom->jclusters[cj0].natoms = atom->iclusters[ci].natoms + atom->iclusters[ci1].natoms;
} }
} }
} }
@ -726,12 +727,8 @@ void binClusters(Atom *atom) {
*/ */
const int nlocal = atom->Nclusters_local; const int nlocal = atom->Nclusters_local;
const int ncji = MAX(1, CLUSTER_M / CLUSTER_N); const int jfac = MAX(1, CLUSTER_N / CLUSTER_M);
if(ncji > 2) { const int ncj = atom->Nclusters_local / jfac;
fprintf(stderr, "binClusters(): ncji cannot be higher than 2!");
}
defineJClusters(atom);
int resize = 1; int resize = 1;
while(resize > 0) { while(resize > 0) {
@ -764,74 +761,65 @@ void binClusters(Atom *atom) {
} }
for(int cg = 0; cg < atom->Nclusters_ghost && !resize; cg++) { for(int cg = 0; cg < atom->Nclusters_ghost && !resize; cg++) {
// Assure we add this j-cluster only once in the bin const int cj = ncj + cg;
if(!(CLUSTER_M < CLUSTER_N && cg % 2)) { int ix = -1, iy = -1;
const int ci = nlocal + cg; MD_FLOAT xtmp, ytmp;
int ix = -1, iy = -1;
MD_FLOAT xtmp, ytmp;
for(int cji = 0; cji < ncji && !resize; cji++) { if(atom->jclusters[cj].natoms > 0) {
int cj = (cji == 0) ? CJ0_FROM_CI(ci) : CJ1_FROM_CI(ci); int cj_vec_base = CJ_VECTOR_BASE_INDEX(cj);
const int n = atom->jclusters[cj].natoms; MD_FLOAT *cj_x = &atom->cl_x[cj_vec_base];
MD_FLOAT cj_minz = atom->jclusters[cj].bbminz;
if(n > 0) { xtmp = cj_x[CL_X_OFFSET + 0];
int cj_vec_base = CJ_VECTOR_BASE_INDEX(cj); ytmp = cj_x[CL_Y_OFFSET + 0];
MD_FLOAT *cj_x = &atom->cl_x[cj_vec_base]; coord2bin2D(xtmp, ytmp, &ix, &iy);
MD_FLOAT cj_minz = atom->jclusters[cj].bbminz; ix = MAX(MIN(ix, mbinx - 1), 0);
iy = MAX(MIN(iy, mbiny - 1), 0);
for(int cjj = 1; cjj < atom->jclusters[cj].natoms; cjj++) {
int nix, niy;
xtmp = cj_x[CL_X_OFFSET + cjj];
ytmp = cj_x[CL_Y_OFFSET + cjj];
coord2bin2D(xtmp, ytmp, &nix, &niy);
nix = MAX(MIN(nix, mbinx - 1), 0);
niy = MAX(MIN(niy, mbiny - 1), 0);
xtmp = cj_x[CL_X_OFFSET + 0]; // Always put the cluster on the bin of its innermost atom so
ytmp = cj_x[CL_Y_OFFSET + 0]; // the cluster should be closer to local clusters
coord2bin2D(xtmp, ytmp, &ix, &iy); if(atom->PBCx[cg] > 0 && ix > nix) { ix = nix; }
ix = MAX(MIN(ix, mbinx - 1), 0); if(atom->PBCx[cg] < 0 && ix < nix) { ix = nix; }
iy = MAX(MIN(iy, mbiny - 1), 0); if(atom->PBCy[cg] > 0 && iy > niy) { iy = niy; }
for(int cjj = 1; cjj < n; cjj++) { if(atom->PBCy[cg] < 0 && iy < niy) { iy = niy; }
int nix, niy; }
xtmp = cj_x[CL_X_OFFSET + cjj];
ytmp = cj_x[CL_Y_OFFSET + cjj];
coord2bin2D(xtmp, ytmp, &nix, &niy);
nix = MAX(MIN(nix, mbinx - 1), 0);
niy = MAX(MIN(niy, mbiny - 1), 0);
// Always put the cluster on the bin of its innermost atom so int bin = iy * mbinx + ix + 1;
// the cluster should be closer to local clusters int c = bin_nclusters[bin];
if(atom->PBCx[cg] > 0 && ix > nix) { ix = nix; } if(c < clusters_per_bin) {
if(atom->PBCx[cg] < 0 && ix < nix) { ix = nix; } // Insert the current ghost cluster in the bin keeping clusters
if(atom->PBCy[cg] > 0 && iy > niy) { iy = niy; } // sorted by z coordinate
if(atom->PBCy[cg] < 0 && iy < niy) { iy = niy; } int inserted = 0;
} for(int i = 0; i < c; i++) {
int last_cl = bin_clusters[bin * clusters_per_bin + i];
if(atom->jclusters[last_cl].bbminz > cj_minz) {
bin_clusters[bin * clusters_per_bin + i] = cj;
int bin = iy * mbinx + ix + 1; for(int j = i + 1; j <= c; j++) {
int c = bin_nclusters[bin]; int tmp = bin_clusters[bin * clusters_per_bin + j];
if(c < clusters_per_bin) { bin_clusters[bin * clusters_per_bin + j] = last_cl;
// Insert the current ghost cluster in the bin keeping clusters last_cl = tmp;
// sorted by z coordinate
int inserted = 0;
for(int i = 0; i < c; i++) {
int last_cl = bin_clusters[bin * clusters_per_bin + i];
if(atom->jclusters[last_cl].bbminz > cj_minz) {
bin_clusters[bin * clusters_per_bin + i] = cj;
for(int j = i + 1; j <= c; j++) {
int tmp = bin_clusters[bin * clusters_per_bin + j];
bin_clusters[bin * clusters_per_bin + j] = last_cl;
last_cl = tmp;
}
inserted = 1;
break;
}
} }
if(!inserted) { inserted = 1;
bin_clusters[bin * clusters_per_bin + c] = cj; break;
}
atom->icluster_bin[ci] = bin;
bin_nclusters[bin]++;
} else {
resize = 1;
} }
} }
if(!inserted) {
bin_clusters[bin * clusters_per_bin + c] = cj;
}
bin_nclusters[bin]++;
} else {
resize = 1;
} }
} }
} }

175
gromacs/pbc.c
View File

@ -47,30 +47,30 @@ void initPbc(Atom* atom) {
/* uses mapping created in setupPbc */ /* uses mapping created in setupPbc */
void updatePbc(Atom *atom, Parameter *param, int firstUpdate) { void updatePbc(Atom *atom, Parameter *param, int firstUpdate) {
DEBUG_MESSAGE("updatePbc start\n"); DEBUG_MESSAGE("updatePbc start\n");
int *border_map = atom->border_map; int jfac = MAX(1, CLUSTER_N / CLUSTER_M);
int nlocal = atom->Nclusters_local; int ncj = atom->Nclusters_local / jfac;
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 cg = 0; cg < atom->Nclusters_ghost; cg++) { for(int cg = 0; cg < atom->Nclusters_ghost; cg++) {
const int ci = nlocal + cg; const int cj = ncj + cg;
int ci_vec_base = CI_VECTOR_BASE_INDEX(ci); int cj_vec_base = CJ_VECTOR_BASE_INDEX(cj);
int bmap_vec_base = CI_VECTOR_BASE_INDEX(border_map[cg]); int bmap_vec_base = CJ_VECTOR_BASE_INDEX(atom->border_map[cg]);
MD_FLOAT *ci_x = &atom->cl_x[ci_vec_base]; MD_FLOAT *cj_x = &atom->cl_x[cj_vec_base];
MD_FLOAT *bmap_x = &atom->cl_x[bmap_vec_base]; MD_FLOAT *bmap_x = &atom->cl_x[bmap_vec_base];
MD_FLOAT bbminx = INFINITY, bbmaxx = -INFINITY; MD_FLOAT bbminx = INFINITY, bbmaxx = -INFINITY;
MD_FLOAT bbminy = INFINITY, bbmaxy = -INFINITY; MD_FLOAT bbminy = INFINITY, bbmaxy = -INFINITY;
MD_FLOAT bbminz = INFINITY, bbmaxz = -INFINITY; MD_FLOAT bbminz = INFINITY, bbmaxz = -INFINITY;
for(int cii = 0; cii < atom->iclusters[ci].natoms; cii++) { for(int cjj = 0; cjj < atom->jclusters[cj].natoms; cjj++) {
MD_FLOAT xtmp = bmap_x[CL_X_OFFSET + cii] + atom->PBCx[cg] * xprd; MD_FLOAT xtmp = bmap_x[CL_X_OFFSET + cjj] + atom->PBCx[cg] * xprd;
MD_FLOAT ytmp = bmap_x[CL_Y_OFFSET + cii] + atom->PBCy[cg] * yprd; MD_FLOAT ytmp = bmap_x[CL_Y_OFFSET + cjj] + atom->PBCy[cg] * yprd;
MD_FLOAT ztmp = bmap_x[CL_Z_OFFSET + cii] + atom->PBCz[cg] * zprd; MD_FLOAT ztmp = bmap_x[CL_Z_OFFSET + cjj] + atom->PBCz[cg] * zprd;
ci_x[CL_X_OFFSET + cii] = xtmp; cj_x[CL_X_OFFSET + cjj] = xtmp;
ci_x[CL_Y_OFFSET + cii] = ytmp; cj_x[CL_Y_OFFSET + cjj] = ytmp;
ci_x[CL_Z_OFFSET + cii] = ztmp; cj_x[CL_Z_OFFSET + cjj] = ztmp;
if(firstUpdate) { if(firstUpdate) {
// TODO: To create the bounding boxes faster, we can use SIMD operations // TODO: To create the bounding boxes faster, we can use SIMD operations
@ -84,18 +84,18 @@ void updatePbc(Atom *atom, Parameter *param, int firstUpdate) {
} }
if(firstUpdate) { if(firstUpdate) {
for(int cii = atom->iclusters[ci].natoms; cii < CLUSTER_M; cii++) { for(int cjj = atom->jclusters[cj].natoms; cjj < CLUSTER_N; cjj++) {
ci_x[CL_X_OFFSET + cii] = INFINITY; cj_x[CL_X_OFFSET + cjj] = INFINITY;
ci_x[CL_Y_OFFSET + cii] = INFINITY; cj_x[CL_Y_OFFSET + cjj] = INFINITY;
ci_x[CL_Z_OFFSET + cii] = INFINITY; cj_x[CL_Z_OFFSET + cjj] = INFINITY;
} }
atom->iclusters[ci].bbminx = bbminx; atom->jclusters[cj].bbminx = bbminx;
atom->iclusters[ci].bbmaxx = bbmaxx; atom->jclusters[cj].bbmaxx = bbmaxx;
atom->iclusters[ci].bbminy = bbminy; atom->jclusters[cj].bbminy = bbminy;
atom->iclusters[ci].bbmaxy = bbmaxy; atom->jclusters[cj].bbmaxy = bbmaxy;
atom->iclusters[ci].bbminz = bbminz; atom->jclusters[cj].bbminz = bbminz;
atom->iclusters[ci].bbmaxz = bbmaxz; atom->jclusters[cj].bbmaxz = bbmaxz;
} }
} }
@ -136,18 +136,18 @@ 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++; \
const int cg = atom->Nclusters_local + Nghost; \ const int cg = ncj + Nghost; \
const int natoms_ci = atom->iclusters[ci].natoms; \ const int cj_natoms = atom->jclusters[cj].natoms; \
border_map[Nghost] = ci; \ atom->border_map[Nghost] = cj; \
atom->PBCx[Nghost] = dx; \ atom->PBCx[Nghost] = dx; \
atom->PBCy[Nghost] = dy; \ atom->PBCy[Nghost] = dy; \
atom->PBCz[Nghost] = dz; \ atom->PBCz[Nghost] = dz; \
atom->iclusters[cg].natoms = natoms_ci; \ atom->jclusters[cg].natoms = cj_natoms; \
Nghost_atoms += natoms_ci; \ Nghost_atoms += cj_natoms; \
int ci_sca_base = CI_SCALAR_BASE_INDEX(ci); \ int cj_sca_base = CJ_SCALAR_BASE_INDEX(cj); \
int cg_sca_base = CI_SCALAR_BASE_INDEX(cg); \ int cg_sca_base = CJ_SCALAR_BASE_INDEX(cg); \
for(int cii = 0; cii < natoms_ci; cii++) { \ for(int cjj = 0; cjj < cj_natoms; cjj++) { \
atom->cl_type[cg_sca_base + cii] = atom->cl_type[ci_sca_base + cii]; \ atom->cl_type[cg_sca_base + cjj] = atom->cl_type[cj_sca_base + cjj]; \
} }
/* internal subroutines */ /* internal subroutines */
@ -163,77 +163,80 @@ void growPbc(Atom* atom) {
void setupPbc(Atom *atom, Parameter *param) { void setupPbc(Atom *atom, Parameter *param) {
DEBUG_MESSAGE("setupPbc start\n"); DEBUG_MESSAGE("setupPbc start\n");
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;
MD_FLOAT zprd = param->zprd; MD_FLOAT zprd = param->zprd;
MD_FLOAT Cutneigh = param->cutneigh; MD_FLOAT Cutneigh = param->cutneigh;
int jfac = MAX(1, CLUSTER_N / CLUSTER_M);
int ncj = atom->Nclusters_local / jfac;
int Nghost = -1; int Nghost = -1;
int Nghost_atoms = 0; int Nghost_atoms = 0;
for(int ci = 0; ci < atom->Nclusters_local; ci++) { for(int cj = 0; cj < ncj; cj++) {
if (atom->Nclusters_local + Nghost + 7 >= atom->Nclusters_max) { if(atom->jclusters[cj].natoms > 0) {
growClusters(atom); if(atom->Nclusters_local + (Nghost + 7) * jfac >= atom->Nclusters_max) {
growClusters(atom);
}
if((Nghost + 7) * jfac >= NmaxGhost) {
growPbc(atom);
}
MD_FLOAT bbminx = atom->jclusters[cj].bbminx;
MD_FLOAT bbmaxx = atom->jclusters[cj].bbmaxx;
MD_FLOAT bbminy = atom->jclusters[cj].bbminy;
MD_FLOAT bbmaxy = atom->jclusters[cj].bbmaxy;
MD_FLOAT bbminz = atom->jclusters[cj].bbminz;
MD_FLOAT bbmaxz = atom->jclusters[cj].bbmaxz;
/* Setup ghost atoms */
/* 6 planes */
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 (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 (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); }
} }
if (Nghost + 7 >= NmaxGhost) {
growPbc(atom);
border_map = atom->border_map;
}
MD_FLOAT bbminx = atom->iclusters[ci].bbminx;
MD_FLOAT bbmaxx = atom->iclusters[ci].bbmaxx;
MD_FLOAT bbminy = atom->iclusters[ci].bbminy;
MD_FLOAT bbmaxy = atom->iclusters[ci].bbmaxy;
MD_FLOAT bbminz = atom->iclusters[ci].bbminz;
MD_FLOAT bbmaxz = atom->iclusters[ci].bbmaxz;
/* Setup ghost atoms */
/* 6 planes */
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 (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 (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); }
} }
if(atom->Nclusters_local + Nghost + 2 >= atom->Nclusters_max) { if(ncj + (Nghost + 1) * jfac >= atom->Nclusters_max) {
growClusters(atom); growClusters(atom);
} }
// Add dummy cluster at the end // Add dummy cluster at the end
int ci_vec_base = CI_VECTOR_BASE_INDEX(atom->Nclusters_local + Nghost + 1); int cj_vec_base = CJ_VECTOR_BASE_INDEX(ncj + Nghost + 1);
MD_FLOAT *ci_x = &atom->cl_x[ci_vec_base]; MD_FLOAT *cj_x = &atom->cl_x[cj_vec_base];
for(int cii = 0; cii < MAX(CLUSTER_M, CLUSTER_N); cii++) { for(int cjj = 0; cjj < CLUSTER_N; cjj++) {
ci_x[CL_X_OFFSET + cii] = INFINITY; cj_x[CL_X_OFFSET + cjj] = INFINITY;
ci_x[CL_Y_OFFSET + cii] = INFINITY; cj_x[CL_Y_OFFSET + cjj] = INFINITY;
ci_x[CL_Z_OFFSET + cii] = INFINITY; cj_x[CL_Z_OFFSET + cjj] = INFINITY;
} }
// increase by one to make it the ghost atom count // increase by one to make it the ghost atom count
atom->dummy_cj = ncj + Nghost + 1;
atom->Nghost = Nghost_atoms; atom->Nghost = Nghost_atoms;
atom->Nclusters_ghost = Nghost + 1; atom->Nclusters_ghost = Nghost + 1;
atom->Nclusters = atom->Nclusters_local + Nghost + 1; atom->Nclusters = atom->Nclusters_local + Nghost + 1;