From a119fcdfdd4bc41848094bba45eeffb819650c51 Mon Sep 17 00:00:00 2001 From: Rafael Ravedutti Date: Thu, 27 Jan 2022 03:07:31 +0100 Subject: [PATCH] Fix some segfaults and add function to update single atoms Signed-off-by: Rafael Ravedutti --- gromacs/atom.c | 2 +- gromacs/includes/neighbor.h | 9 +- gromacs/main.c | 20 ++- gromacs/neighbor.c | 341 +++++++++++++++++++++--------------- gromacs/pbc.c | 13 +- 5 files changed, 225 insertions(+), 160 deletions(-) diff --git a/gromacs/atom.c b/gromacs/atom.c index 3c63821..1567f74 100644 --- a/gromacs/atom.c +++ b/gromacs/atom.c @@ -286,5 +286,5 @@ 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); + atom->cl_x = (MD_FLOAT*) reallocate(atom->cl_x, ALIGNMENT, atom->Nclusters_max * CLUSTER_DIM_N * 3 * sizeof(MD_FLOAT), nold * CLUSTER_DIM_N * 3 * sizeof(MD_FLOAT)); } diff --git a/gromacs/includes/neighbor.h b/gromacs/includes/neighbor.h index 04b7e3b..1b83efb 100644 --- a/gromacs/includes/neighbor.h +++ b/gromacs/includes/neighbor.h @@ -34,10 +34,11 @@ typedef struct { } Neighbor; extern void initNeighbor(Neighbor*, Parameter*); -extern void setupNeighbor(Parameter*); +extern void setupNeighbor(Parameter*, Atom*); extern void binatoms(Atom*); -extern void buildNeighbor(Parameter*, Atom*, Neighbor*); +extern void buildNeighbor(Atom*, Neighbor*); extern void sortAtom(Atom*); -extern void buildClusters(Parameter*, Atom*); -extern void binGhostClusters(Parameter*, Atom*); +extern void buildClusters(Atom*); +extern void binClusters(Atom*); +extern void updateSingleAtoms(Atom*); #endif diff --git a/gromacs/main.c b/gromacs/main.c index 15f606b..ca18756 100644 --- a/gromacs/main.c +++ b/gromacs/main.c @@ -96,12 +96,13 @@ double setup( } else { readAtom(atom, param); } - setupNeighbor(param); + setupNeighbor(param, atom); setupThermo(param, atom->Natoms); if(param->input_file == NULL) { adjustThermo(param, atom); } + buildClusters(atom); setupPbc(atom, param); - buildClusters(param, atom); - buildNeighbor(param, atom, neighbor); + binClusters(atom); + buildNeighbor(atom, neighbor); E = getTimeStamp(); return E-S; @@ -116,9 +117,12 @@ double reneighbour( S = getTimeStamp(); LIKWID_MARKER_START("reneighbour"); + updateSingleAtoms(atom); updateAtomsPbc(atom, param); + buildClusters(atom); setupPbc(atom, param); - buildNeighbor(param, atom, neighbor); + binClusters(atom); + buildNeighbor(atom, neighbor); LIKWID_MARKER_STOP("reneighbour"); E = getTimeStamp(); @@ -253,11 +257,13 @@ int main(int argc, char** argv) #if defined(MEM_TRACER) || defined(INDEX_TRACER) traceAddresses(¶m, &atom, &neighbor, n + 1); #endif + /* if(param.force_field == FF_EAM) { timer[FORCE] = computeForceEam(&eam, ¶m, &atom, &neighbor, &stats); } else { timer[FORCE] = computeForceLJ(¶m, &atom, &neighbor, &stats); } + */ timer[NEIGH] = 0.0; timer[TOTAL] = getTimeStamp(); @@ -267,7 +273,7 @@ int main(int argc, char** argv) } for(int n = 0; n < param.ntimes; n++) { - initialIntegrate(¶m, &atom); + //initialIntegrate(¶m, &atom); if((n + 1) % param.every) { updatePbc(&atom, ¶m); @@ -279,13 +285,15 @@ int main(int argc, char** argv) traceAddresses(¶m, &atom, &neighbor, n + 1); #endif + /* if(param.force_field == FF_EAM) { timer[FORCE] += computeForceEam(&eam, ¶m, &atom, &neighbor, &stats); } else { timer[FORCE] += computeForceLJ(¶m, &atom, &neighbor, &stats); } + */ - finalIntegrate(¶m, &atom); + //finalIntegrate(¶m, &atom); if(!((n + 1) % param.nstat) && (n+1) < param.ntimes) { computeThermo(n + 1, ¶m, &atom); diff --git a/gromacs/neighbor.c b/gromacs/neighbor.c index 2f2ebc9..11a438c 100644 --- a/gromacs/neighbor.c +++ b/gromacs/neighbor.c @@ -32,7 +32,7 @@ #define SMALL 1.0e-6 #define FACTOR 0.999 -static MD_FLOAT xprd, yprd; +static MD_FLOAT xprd, yprd, zprd; static MD_FLOAT bininvx, bininvy; static int mbinxlo, mbinylo; static int nbinx, nbiny; @@ -55,14 +55,12 @@ static int coord2bin(MD_FLOAT, MD_FLOAT); static MD_FLOAT bindist(int, int); /* exported subroutines */ -void initNeighbor(Neighbor *neighbor, Parameter *param) -{ +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; nmax = 0; atoms_per_bin = 8; clusters_per_bin = (atoms_per_bin / CLUSTER_DIM_N) + 4; @@ -76,7 +74,7 @@ void initNeighbor(Neighbor *neighbor, Parameter *param) neighbor->neighbors = NULL; } -void setupNeighbor(Parameter* param) { +void setupNeighbor(Parameter *param, Atom *atom) { MD_FLOAT coord; int mbinxhi, mbinyhi; int nextx, nexty, nextz; @@ -84,30 +82,31 @@ void setupNeighbor(Parameter* param) { 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; + MD_FLOAT atom_density = ((MD_FLOAT)(atom->Nlocal)) / ((xhi - xlo) * (yhi - ylo) * (zhi - zlo)); + MD_FLOAT atoms_in_cell = MAX(CLUSTER_DIM_M, CLUSTER_DIM_N); + binsizex = cbrt(atoms_in_cell / atom_density); + binsizey = cbrt(atoms_in_cell / atom_density); cutneighsq = cutneigh * cutneigh; + nbinx = (int)((xhi - xlo) / binsizex); + nbiny = (int)((yhi - ylo) / binsizey); + if(nbinx == 0) { nbinx = 1; } + if(nbiny == 0) { nbiny = 1; } + bininvx = 1.0 / binsizex; + bininvy = 1.0 / binsizey; - if(param->input_file != NULL) { - binsizex = cutneigh * 0.5; - binsizey = cutneigh * 0.5; - nbinx = (int)((param->xhi - param->xlo) / binsizex); - nbiny = (int)((param->yhi - param->ylo) / binsizey); - if(nbinx == 0) { nbinx = 1; } - if(nbiny == 0) { nbiny = 1; } - bininvx = nbinx / (param->xhi - param->xlo); - bininvy = nbiny / (param->yhi - param->ylo); - } else { - binsizex = xprd / nbinx; - binsizey = yprd / nbiny; - bininvx = 1.0 / binsizex; - bininvy = 1.0 / binsizey; - } - + printf("hi = %f, %f\n", xhi, yhi); + printf("atom_density = %f\n", atom_density); + printf("atoms_in_cell = %f\n", atoms_in_cell); + printf("binsize = %f, %f\n", binsizex, binsizey); + printf("nbin = %d, %d\n", nbinx, nbiny); coord = xlo - cutneigh - SMALL * xprd; mbinxlo = (int) (coord * bininvx); if (coord < 0.0) { @@ -207,8 +206,9 @@ int atomDistanceInRange(Atom *atom, int ci, int cj, MD_FLOAT rsq) { return 0; } -void buildNeighbor(Parameter *param, Atom *atom, Neighbor *neighbor) { - int nall = atom->Nlocal + atom->Nghost; +void buildNeighbor(Atom *atom, Neighbor *neighbor) { + printf("buildNeighbor start\n"); + int nall = atom->Nclusters_local + atom->Nclusters_ghost; /* extend atom arrays if necessary */ if(nall > nmax) { @@ -239,33 +239,35 @@ void buildNeighbor(Parameter *param, Atom *atom, Neighbor *neighbor) { int *loc_bin = &cluster_bins[jbin * clusters_per_bin]; int cj, m = -1; MD_FLOAT jbb_zmin, jbb_zmax; + const int c = cluster_bincount[jbin]; - 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(c > 0) { + 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 < cluster_bincount[jbin] && jbb_zmax - ibb_zmax < cutneighsq) { - if((jbb_zmin - ibb_zmax) * (jbb_zmin - ibb_zmax) > cutneighsq) { - break; - } + while(m < c && jbb_zmax - 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) { - if(d_bb_sq < rBB || atomDistanceInRange(atom, ci, cj, cutneighsq)) { - if(cj == ci) { - // Add to neighbor list with diagonal interaction mask + double d_bb_sq = getBoundingBoxDistanceSq(atom, ci, cj); + if(d_bb_sq < cutneighsq) { + if(d_bb_sq < rBB || atomDistanceInRange(atom, ci, cj, cutneighsq)) { + neighptr[n++] = cj; } - neighptr[n++] = cj; + } + + m++; + if(m < c) { + cj = loc_bin[m]; + jbb_zmin = atom->clusters[cj].bbminz; + jbb_zmax = atom->clusters[cj].bbmaxz; } } - - m++; - cj = loc_bin[m]; - jbb_zmin = atom->clusters[cj].bbminz; - jbb_zmax = atom->clusters[cj].bbmaxz; } } @@ -287,6 +289,7 @@ void buildNeighbor(Parameter *param, Atom *atom, Neighbor *neighbor) { neighbor->neighbors = (int*) malloc(atom->Nmax * neighbor->maxneighs * sizeof(int)); } } + printf("buildNeighbor end\n"); } /* internal subroutines */ @@ -352,7 +355,8 @@ void coord2bin2D(MD_FLOAT xin, MD_FLOAT yin, int *ix, int *iy) { } } -void binatoms(Atom *atom) { +void binAtoms(Atom *atom) { + printf("binatoms start\n"); int nall = atom->Nlocal + atom->Nghost; int resize = 1; @@ -380,142 +384,189 @@ void binatoms(Atom *atom) { bins = (int*) malloc(mbins * atoms_per_bin * sizeof(int)); } } + printf("binatoms end\n"); } // TODO: Use pigeonhole sorting -void sortBinAtomsByZCoord(Parameter *param, Atom *atom, int bin) { - int c = bincount[bin]; - int *bin_ptr = &bins[bin * atoms_per_bin]; +void sortAtomsByZCoord(Atom *atom) { + for(int bin = 0; bin < mbins; bin++) { + int c = bincount[bin]; + int *bin_ptr = &bins[bin * atoms_per_bin]; - for(int ac_i = 0; ac_i < c; ac_i++) { - int i = bin_ptr[ac_i]; - int min_ac = ac_i; - int min_idx = i; - MD_FLOAT min_z = atom_z(i); + for(int ac_i = 0; ac_i < c; ac_i++) { + int i = bin_ptr[ac_i]; + int min_ac = ac_i; + int min_idx = i; + MD_FLOAT min_z = atom_z(i); - for(int ac_j = ac_i + 1; ac_j < c; ac_j++) { - int j = bin_ptr[ac_j]; - MD_FLOAT zj = atom_z(j); - if(zj < min_z) { - min_ac = zj; - min_idx = j; - min_z = zj; + for(int ac_j = ac_i + 1; ac_j < c; ac_j++) { + int j = bin_ptr[ac_j]; + MD_FLOAT zj = atom_z(j); + if(zj < min_z) { + min_ac = zj; + min_idx = j; + min_z = zj; + } } - } - bin_ptr[ac_i] = min_idx; - bin_ptr[min_ac] = i; + bin_ptr[ac_i] = min_idx; + bin_ptr[min_ac] = i; + } } } -void buildClusters(Parameter *param, Atom *atom) { +void buildClusters(Atom *atom) { + printf("buildClusters start\n"); + atom->Nclusters_local = 0; + /* bin local atoms */ - binatoms(atom); + binAtoms(atom); + sortAtomsByZCoord(atom); for(int bin = 0; bin < mbins; bin++) { - sortBinAtomsByZCoord(param, atom, bin); - } + int c = bincount[bin]; + int ac = 0; + const int nclusters = ((c + CLUSTER_DIM_N - 1) / CLUSTER_DIM_N); + for(int cl = 0; cl < nclusters; cl++) { + const int ci = atom->Nclusters_local; - 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(ci >= atom->Nclusters_max) { + growClusters(atom); + } - 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; + 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); + 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; + 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; } + // 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++; + 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; } - cluster_bincount[bin] = nclusters; + ac++; + } + + 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++; + } + } + + printf("buildClusters end\n"); +} + +void binClusters(Atom *atom) { + printf("binClusters start\n"); + const int nlocal = atom->Nclusters_local; + int resize = 1; + + while(resize > 0) { + resize = 0; + + for(int bin = 0; bin < mbins; bin++) { + cluster_bincount[bin] = 0; + } + + printf("local\n"); + for(int ci = 0; ci < nlocal && !resize; ci++) { + int bin = atom->clusters[ci].bin; + int c = cluster_bincount[bin]; + if(c < clusters_per_bin) { + cluster_bins[bin * clusters_per_bin + c] = ci; + cluster_bincount[bin]++; } else { resize = 1; } } + printf("ghost\n"); + for(int ci = 0; ci < atom->Nclusters_ghost && !resize; 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(atom->PBCx[ci] > 0 && ix > nix) { ix = nix; } + if(atom->PBCx[ci] < 0 && ix < nix) { ix = nix; } + if(atom->PBCy[ci] > 0 && iy > niy) { iy = niy; } + if(atom->PBCy[ci] < 0 && iy < niy) { iy = niy; } + } + + int bin = iy * mbinx + ix + 1; + int c = cluster_bincount[bin]; + if(c < clusters_per_bin) { + cluster_bins[bin * clusters_per_bin + c] = nlocal + ci; + cluster_bincount[bin]++; + } else { + resize = 1; + } + } + + printf("end\n"); if(resize) { free(cluster_bins); clusters_per_bin *= 2; cluster_bins = (int*) malloc(mbins * clusters_per_bin * sizeof(int)); } } + + printf("binClusters end\n"); } -void binGhostClusters(Parameter *param, Atom *atom) { - int nlocal = atom->Nclusters_local; +void updateSingleAtoms(Atom *atom) { + int Natom = 0; - 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; + for(int ci = 0; ci < atom->Nclusters_local; ci++) { + MD_FLOAT *cptr = cluster_ptr(ci); - 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(atom->PBCx[ci] > 0 && ix > nix) { ix = nix; } - if(atom->PBCx[ci] < 0 && ix < nix) { ix = nix; } - if(atom->PBCy[ci] > 0 && iy > niy) { iy = niy; } - if(atom->PBCy[ci] < 0 && iy < niy) { iy = niy; } + for(int cii = 0; cii < atom->clusters[ci].natoms; cii++) { + atom_x(Natom) = cluster_x(cptr, cii); + atom_y(Natom) = cluster_y(cptr, cii); + atom_z(Natom) = cluster_z(cptr, cii); + Natom++; } + } - int bin = iy * mbinx + ix + 1; - int c = cluster_bincount[bin]; - cluster_bins[bin * clusters_per_bin + c] = nlocal + ci; - cluster_bincount[bin]++; + if(Natom != atom->Nlocal) { + fprintf(stderr, "updateSingleAtoms(): Number of atoms changed!\n"); } } diff --git a/gromacs/pbc.c b/gromacs/pbc.c index 6b62708..b2a3d2b 100644 --- a/gromacs/pbc.c +++ b/gromacs/pbc.c @@ -44,6 +44,7 @@ void initPbc(Atom* atom) { /* update coordinates of ghost atoms */ /* uses mapping created in setupPbc */ void updatePbc(Atom *atom, Parameter *param) { + printf("updatePbc start\n"); int *border_map = atom->border_map; int nlocal = atom->Nclusters_local; MD_FLOAT xprd = param->xprd; @@ -60,11 +61,13 @@ void updatePbc(Atom *atom, Parameter *param) { cluster_z(cptr, cii) = cluster_z(bmap_cptr, cii) + atom->PBCz[ci] * zprd; } } + printf("updatePbc end\n"); } /* relocate atoms that have left domain according * to periodic boundary conditions */ void updateAtomsPbc(Atom *atom, Parameter *param) { + printf("updateAtomsPbc start\n"); MD_FLOAT xprd = param->xprd; MD_FLOAT yprd = param->yprd; MD_FLOAT zprd = param->zprd; @@ -88,6 +91,7 @@ void updateAtomsPbc(Atom *atom, Parameter *param) { atom_z(i) -= zprd; } } + printf("updateAtomsPbc end\n"); } /* setup periodic boundary conditions by @@ -120,6 +124,7 @@ void growPbc(Atom* atom) { } void setupPbc(Atom *atom, Parameter *param) { + printf("setupPbc start\n"); int *border_map = atom->border_map; MD_FLOAT xprd = param->xprd; MD_FLOAT yprd = param->yprd; @@ -128,11 +133,11 @@ void setupPbc(Atom *atom, Parameter *param) { int Nghost = -1; for(int ci = 0; ci < atom->Nclusters_local; ci++) { - if (atom->Nclusters_local + atom->Nclusters_ghost + 7 >= atom->Nclusters_max) { + if (atom->Nclusters_local + Nghost + 7 >= atom->Nclusters_max) { growClusters(atom); } - if (atom->Nclusters_ghost + 7 >= NmaxGhost) { + if (Nghost + 7 >= NmaxGhost) { growPbc(atom); border_map = atom->border_map; } @@ -143,6 +148,7 @@ void setupPbc(Atom *atom, Parameter *param) { 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 (bbminx < Cutneigh) { ADDGHOST(+1,0,0); } @@ -179,7 +185,6 @@ void setupPbc(Atom *atom, Parameter *param) { atom->Nclusters_ghost = Nghost + 1; atom->Nclusters = atom->Nclusters_local + Nghost + 1; - // Update and bin created ghost clusters + // Update created ghost clusters positions updatePbc(atom, param); - binGhostClusters(param, atom); }