From 85e7954932c07f64b7a7066a8d7de07891a95435 Mon Sep 17 00:00:00 2001 From: Rafael Ravedutti Date: Tue, 1 Feb 2022 20:16:04 +0100 Subject: [PATCH] Check all clusters in cell when building neighbor lists because ghost clusters may not be sorted Signed-off-by: Rafael Ravedutti --- gromacs/force_lj.c | 15 ++----- gromacs/includes/util.h | 5 +++ gromacs/main.c | 12 ++++-- gromacs/neighbor.c | 91 +++++++++++++++++++++++++++-------------- 4 files changed, 77 insertions(+), 46 deletions(-) diff --git a/gromacs/force_lj.c b/gromacs/force_lj.c index 7ff61cb..f55dbf0 100644 --- a/gromacs/force_lj.c +++ b/gromacs/force_lj.c @@ -28,9 +28,11 @@ #include #include #include +#include + double computeForceLJ(Parameter *param, Atom *atom, Neighbor *neighbor, Stats *stats) { - fprintf(stdout, "computeForceLJ begin\n"); + DEBUG_MESSAGE("computeForceLJ begin\n"); int Nlocal = atom->Nlocal; int* neighs; MD_FLOAT cutforcesq = param->cutforce * param->cutforce; @@ -80,11 +82,6 @@ double computeForceLJ(Parameter *param, Atom *atom, Neighbor *neighbor, Stats *s fix += delx * force; fiy += dely * force; fiz += delz * force; - /* - if(force < -50.0 || force > 50.0) { - fprintf(stdout, "%d-%d/%d-%d: %f, %f, %f ---- %f, %f, %f, %f\n", ci, cii, cj, cjj, xtmp, ytmp, ztmp, fix, fiy, fiz, force); - } - */ } } } @@ -92,10 +89,6 @@ double computeForceLJ(Parameter *param, Atom *atom, Neighbor *neighbor, Stats *s cluster_x(cifptr, cii) += fix; cluster_y(cifptr, cii) += fiy; cluster_z(cifptr, cii) += fiz; - /* - if(fix < -100.0 || fix > 100.0 || fiy < -100.0 || fiy > 100.0 || fiz < -100.0 || fiz > 100.0) { - fprintf(stdout, "%d-%d: %f, %f, %f ---- %f, %f, %f\n", ci, cii, xtmp, ytmp, ztmp, fix, fiy, fiz); - }*/ } } @@ -105,6 +98,6 @@ double computeForceLJ(Parameter *param, Atom *atom, Neighbor *neighbor, Stats *s LIKWID_MARKER_STOP("force"); double E = getTimeStamp(); - fprintf(stdout, "computeForceLJ end\n"); + DEBUG_MESSAGE("computeForceLJ end\n"); return E-S; } diff --git a/gromacs/includes/util.h b/gromacs/includes/util.h index 20d8806..c777282 100644 --- a/gromacs/includes/util.h +++ b/gromacs/includes/util.h @@ -32,6 +32,11 @@ #ifndef ABS #define ABS(a) ((a) >= 0 ? (a) : -(a)) #endif +#ifdef DEBUG +#define DEBUG_MESSAGE(msg) printf +#else +#define DEBUG_MESSAGE(msg) +#endif #define FF_LJ 0 #define FF_EAM 1 diff --git a/gromacs/main.c b/gromacs/main.c index 3556c6a..27cb50e 100644 --- a/gromacs/main.c +++ b/gromacs/main.c @@ -118,7 +118,8 @@ double reneighbour(Parameter *param, Atom *atom, Neighbor *neighbor) { } void initialIntegrate(Parameter *param, Atom *atom) { - fprintf(stdout, "initialIntegrate start\n"); + DEBUG_MESSAGE("initialIntegrate start\n"); + for(int ci = 0; ci < atom->Nclusters_local; ci++) { MD_FLOAT *ciptr = cluster_pos_ptr(ci); MD_FLOAT *civptr = cluster_velocity_ptr(ci); @@ -133,11 +134,13 @@ void initialIntegrate(Parameter *param, Atom *atom) { cluster_z(ciptr, cii) += param->dt * cluster_z(civptr, cii); } } - fprintf(stdout, "initialIntegrate end\n"); + + DEBUG_MESSAGE("initialIntegrate end\n"); } void finalIntegrate(Parameter *param, Atom *atom) { - fprintf(stdout, "finalIntegrate start\n"); + DEBUG_MESSAGE("finalIntegrate start\n"); + for(int ci = 0; ci < atom->Nclusters_local; ci++) { MD_FLOAT *civptr = cluster_velocity_ptr(ci); MD_FLOAT *cifptr = cluster_force_ptr(ci); @@ -148,7 +151,8 @@ void finalIntegrate(Parameter *param, Atom *atom) { cluster_z(civptr, cii) += param->dtforce * cluster_z(cifptr, cii); } } - fprintf(stdout, "finalIntegrate end\n"); + + DEBUG_MESSAGE("finalIntegrate end\n"); } void printAtomState(Atom *atom) { diff --git a/gromacs/neighbor.c b/gromacs/neighbor.c index 5a3f88d..e5e2118 100644 --- a/gromacs/neighbor.c +++ b/gromacs/neighbor.c @@ -202,7 +202,7 @@ int atomDistanceInRange(Atom *atom, int ci, int cj, MD_FLOAT rsq) { } void buildNeighbor(Atom *atom, Neighbor *neighbor) { - printf("buildNeighbor start\n"); + DEBUG_MESSAGE("buildNeighbor start\n"); int nall = atom->Nclusters_local + atom->Nclusters_ghost; /* extend atom arrays if necessary */ @@ -240,17 +240,17 @@ void buildNeighbor(Atom *atom, Neighbor *neighbor) { const int c = cluster_bincount[jbin]; if(c > 0) { - do { + //do { m++; cj = loc_bin[m]; jbb_zmin = atom->clusters[cj].bbminz; jbb_zmax = atom->clusters[cj].bbmaxz; - } while(m + 1 < c && (ibb_zmin - jbb_zmax) * (ibb_zmin - jbb_zmax) > cutneighsq); + //} while(m + 1 < c && (ibb_zmin - jbb_zmax) * (ibb_zmin - jbb_zmax) > cutneighsq); while(m < c) { - if((jbb_zmin - ibb_zmax) * (jbb_zmin - ibb_zmax) > cutneighsq) { + /*if((jbb_zmin - ibb_zmax) * (jbb_zmin - ibb_zmax) > cutneighsq) { break; - } + }*/ double d_bb_sq = getBoundingBoxDistanceSq(atom, ci, cj); if(d_bb_sq < cutneighsq) { @@ -280,7 +280,7 @@ void buildNeighbor(Atom *atom, Neighbor *neighbor) { } if(resize) { - printf("RESIZE %d\n", neighbor->maxneighs); + fprintf(stdout, "RESIZE %d\n", neighbor->maxneighs); neighbor->maxneighs = new_maxneighs * 1.2; free(neighbor->neighbors); neighbor->neighbors = (int*) malloc(atom->Nmax * neighbor->maxneighs * sizeof(int)); @@ -288,28 +288,46 @@ void buildNeighbor(Atom *atom, Neighbor *neighbor) { } /* - printf("\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++) { MD_FLOAT *ciptr = cluster_pos_ptr(ci); int* neighptr = &(neighbor->neighbors[ci * neighbor->maxneighs]); - printf("Cluster %d, bbx = {%f, %f}, bby = {%f, %f}, bbz = {%f, %f}\n", ci, atom->clusters[ci].bbminx, atom->clusters[ci].bbmaxx, atom->clusters[ci].bbminy, atom->clusters[ci].bbmaxy, atom->clusters[ci].bbminz, atom->clusters[ci].bbmaxz); + + DEBUG_MESSAGE("Cluster %d, bbx = {%f, %f}, bby = {%f, %f}, bbz = {%f, %f}\n", + ci, + atom->clusters[ci].bbminx, + atom->clusters[ci].bbmaxx, + atom->clusters[ci].bbminy, + atom->clusters[ci].bbmaxy, + atom->clusters[ci].bbminz, + atom->clusters[ci].bbmaxz); + for(int cii = 0; cii < CLUSTER_DIM_N; cii++) { - fprintf(stdout, "%f, %f, %f\n", cluster_x(ciptr, cii), cluster_y(ciptr, cii), cluster_z(ciptr, cii)); + DEBUG_MESSAGE("%f, %f, %f\n", cluster_x(ciptr, cii), cluster_y(ciptr, cii), cluster_z(ciptr, cii)); } - printf("Neighbors:\n"); + DEBUG_MESSAGE("Neighbors:\n"); for(int k = 0; k < neighbor->numneigh[ci]; k++) { const int cj = neighptr[k]; MD_FLOAT *cjptr = cluster_pos_ptr(cj); - printf(" Cluster %d, bbx = {%f, %f}, bby = {%f, %f}, bbz = {%f, %f}\n", cj, atom->clusters[cj].bbminx, atom->clusters[cj].bbmaxx, atom->clusters[cj].bbminy, atom->clusters[cj].bbmaxy, atom->clusters[cj].bbminz, atom->clusters[cj].bbmaxz); + + DEBUG_MESSAGE(" Cluster %d, bbx = {%f, %f}, bby = {%f, %f}, bbz = {%f, %f}\n", + cj, + atom->clusters[cj].bbminx, + atom->clusters[cj].bbmaxx, + atom->clusters[cj].bbminy, + atom->clusters[cj].bbmaxy, + atom->clusters[cj].bbminz, + atom->clusters[cj].bbmaxz); + for(int cjj = 0; cjj < CLUSTER_DIM_N; cjj++) { - fprintf(stdout, " %f, %f, %f\n", cluster_x(cjptr, cjj), cluster_y(cjptr, cjj), cluster_z(cjptr, cjj)); + DEBUG_MESSAGE(" %f, %f, %f\n", cluster_x(cjptr, cjj), cluster_y(cjptr, cjj), cluster_z(cjptr, cjj)); } } } */ - printf("buildNeighbor end\n"); + DEBUG_MESSAGE("buildNeighbor end\n"); } /* internal subroutines */ @@ -376,7 +394,7 @@ void coord2bin2D(MD_FLOAT xin, MD_FLOAT yin, int *ix, int *iy) { } void binAtoms(Atom *atom) { - printf("binAtoms start\n"); + DEBUG_MESSAGE("binAtoms start\n"); int resize = 1; while(resize > 0) { @@ -388,7 +406,6 @@ void binAtoms(Atom *atom) { for(int i = 0; i < atom->Nlocal; i++) { int ibin = coord2bin(atom_x(i), atom_y(i)); - if(ibin < 0 || ibin > mbins) { fprintf(stderr, "%d: %f, %f\n", i, atom_x(i), atom_y(i)); } if(bincount[ibin] < atoms_per_bin) { int ac = bincount[ibin]++; bins[ibin * atoms_per_bin + ac] = i; @@ -403,12 +420,13 @@ void binAtoms(Atom *atom) { bins = (int*) malloc(mbins * atoms_per_bin * sizeof(int)); } } - printf("binAtoms end\n"); + + DEBUG_MESSAGE("binAtoms end\n"); } // TODO: Use pigeonhole sorting void sortAtomsByZCoord(Atom *atom) { - printf("sortAtomsByZCoord start\n"); + DEBUG_MESSAGE("sortAtomsByZCoord start\n"); for(int bin = 0; bin < mbins; bin++) { int c = bincount[bin]; int *bin_ptr = &bins[bin * atoms_per_bin]; @@ -433,11 +451,12 @@ void sortAtomsByZCoord(Atom *atom) { bin_ptr[min_ac] = i; } } - printf("sortAtomsByZCoord end\n"); + + DEBUG_MESSAGE("sortAtomsByZCoord end\n"); } void buildClusters(Atom *atom) { - printf("buildClusters start\n"); + DEBUG_MESSAGE("buildClusters start\n"); atom->Nclusters_local = 0; /* bin local atoms */ @@ -506,20 +525,29 @@ void buildClusters(Atom *atom) { } } - printf("buildClusters end\n"); + DEBUG_MESSAGE("buildClusters end\n"); } void binClusters(Atom *atom) { - printf("binClusters start\n"); + DEBUG_MESSAGE("binClusters start\n"); /* - fprintf(stdout, "Nghost = %d\n", atom->Nclusters_ghost); + DEBUG_MESSAGE("Nghost = %d\n", atom->Nclusters_ghost); for(int ci = atom->Nclusters_local; ci < atom->Nclusters_local + 4; ci++) { MD_FLOAT *cptr = cluster_pos_ptr(ci); - fprintf(stdout, "Cluster %d:\n", ci); - fprintf(stdout, "bin=%d, Natoms=%d, bbox={%f,%f},{%f,%f},{%f,%f}\n", atom->clusters[ci].bin, atom->clusters[ci].natoms, atom->clusters[ci].bbminx, atom->clusters[ci].bbmaxx, atom->clusters[ci].bbminy, atom->clusters[ci].bbmaxy, atom->clusters[ci].bbminz, atom->clusters[ci].bbmaxz); + DEBUG_MESSAGE("Cluster %d:\n", ci); + DEBUG_MESSAGE("bin=%d, Natoms=%d, bbox={%f,%f},{%f,%f},{%f,%f}\n", + atom->clusters[ci].bin, + atom->clusters[ci].natoms, + atom->clusters[ci].bbminx, + atom->clusters[ci].bbmaxx, + atom->clusters[ci].bbminy, + atom->clusters[ci].bbmaxy, + atom->clusters[ci].bbminz, + atom->clusters[ci].bbmaxz); + for(int cii = 0; cii < CLUSTER_DIM_N; cii++) { - fprintf(stdout, "%f, %f, %f\n", cluster_x(cptr, cii), cluster_y(cptr, cii), cluster_z(cptr, cii)); + DEBUG_MESSAGE("%f, %f, %f\n", cluster_x(cptr, cii), cluster_y(cptr, cii), cluster_z(cptr, cii)); } } */ @@ -590,16 +618,16 @@ void binClusters(Atom *atom) { } /* - fprintf(stdout, "cluster_bincount\n"); - for(int i = 0; i < mbins; i++) { fprintf(stdout, "%d, ", cluster_bincount[i]); } - fprintf(stdout, "\n"); + DEBUG_MESSAGE("cluster_bincount\n"); + for(int i = 0; i < mbins; i++) { DEBUG_MESSAGE("%d, ", cluster_bincount[i]); } + DEBUG_MESSAGE("\n"); */ - printf("binClusters stop\n"); + DEBUG_MESSAGE("binClusters stop\n"); } void updateSingleAtoms(Atom *atom) { - printf("updateSingleAtoms start\n"); + DEBUG_MESSAGE("updateSingleAtoms start\n"); int Natom = 0; for(int ci = 0; ci < atom->Nclusters_local; ci++) { @@ -620,5 +648,6 @@ void updateSingleAtoms(Atom *atom) { if(Natom != atom->Nlocal) { fprintf(stderr, "updateSingleAtoms(): Number of atoms changed!\n"); } - printf("updateSingleAtoms stop\n"); + + DEBUG_MESSAGE("updateSingleAtoms stop\n"); }