From e7737e915140827530dae991f23340e4fe7f4cff Mon Sep 17 00:00:00 2001 From: Rafael Ravedutti Date: Fri, 18 Mar 2022 01:28:11 +0100 Subject: [PATCH] Refactor half neighbor lists code Signed-off-by: Rafael Ravedutti --- lammps/force_lj.c | 16 +++++----- lammps/includes/neighbor.h | 4 +-- lammps/main.c | 4 +++ lammps/neighbor.c | 63 ++++++++++---------------------------- lammps/parameter.c | 2 +- 5 files changed, 33 insertions(+), 56 deletions(-) diff --git a/lammps/force_lj.c b/lammps/force_lj.c index e6b6ec0..68ae157 100644 --- a/lammps/force_lj.c +++ b/lammps/force_lj.c @@ -28,8 +28,7 @@ #include #include -double computeForceLJFullNeigh(Parameter *param, Atom *atom, Neighbor *neighbor, Stats *stats) -{ +double computeForceLJFullNeigh(Parameter *param, Atom *atom, Neighbor *neighbor, Stats *stats) { int Nlocal = atom->Nlocal; int* neighs; MD_FLOAT* fx = atom->fx; @@ -108,8 +107,7 @@ double computeForceLJFullNeigh(Parameter *param, Atom *atom, Neighbor *neighbor, return E-S; } -double computeForceLJHalfNeigh(Parameter *param, Atom *atom, Neighbor *neighbor, Stats *stats) -{ +double computeForceLJHalfNeigh(Parameter *param, Atom *atom, Neighbor *neighbor, Stats *stats) { int Nlocal = atom->Nlocal; int* neighs; MD_FLOAT* fx = atom->fx; @@ -167,9 +165,13 @@ double computeForceLJHalfNeigh(Parameter *param, Atom *atom, Neighbor *neighbor, fiy += dely * force; fiz += delz * force; - fx[j] -= delx * force; - fy[j] -= dely * force; - fz[j] -= delz * force; + // We do not need to update forces for ghost atoms + if(j < Nlocal) { + // TODO: Experiment with AoS layout for forces + fx[j] -= delx * force; + fy[j] -= dely * force; + fz[j] -= delz * force; + } } } diff --git a/lammps/includes/neighbor.h b/lammps/includes/neighbor.h index 86df386..193b848 100644 --- a/lammps/includes/neighbor.h +++ b/lammps/includes/neighbor.h @@ -28,10 +28,10 @@ typedef struct { int every; int ncalls; - int* neighbors; int maxneighs; + int half_neigh; + int* neighbors; int* numneigh; - int halfneigh; } Neighbor; extern void initNeighbor(Neighbor*, Parameter*); diff --git a/lammps/main.c b/lammps/main.c index 8bda40f..9347369 100644 --- a/lammps/main.c +++ b/lammps/main.c @@ -181,6 +181,10 @@ int main(int argc, char** argv) { param.nz = atoi(argv[++i]); continue; } + if((strcmp(argv[i], "-half") == 0)) { + param.half_neigh = atoi(argv[++i]); + continue; + } if((strcmp(argv[i], "-r") == 0) || (strcmp(argv[i], "--radius") == 0)) { param.cutforce = atof(argv[++i]); continue; diff --git a/lammps/neighbor.c b/lammps/neighbor.c index a7c684e..8900a04 100644 --- a/lammps/neighbor.c +++ b/lammps/neighbor.c @@ -46,14 +46,11 @@ 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 int halfneigh; //use half neighbour list - static int coord2bin(MD_FLOAT, MD_FLOAT , MD_FLOAT); static MD_FLOAT bindist(int, 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; @@ -70,11 +67,10 @@ void initNeighbor(Neighbor *neighbor, Parameter *param) neighbor->maxneighs = 100; neighbor->numneigh = NULL; neighbor->neighbors = NULL; - neighbor->halfneigh = param->half_neigh; + neighbor->half_neigh = param->half_neigh; } -void setupNeighbor(Parameter* param) -{ +void setupNeighbor(Parameter* param) { MD_FLOAT coord; int mbinxhi, mbinyhi, mbinzhi; int nextx, nexty, nextz; @@ -116,25 +112,19 @@ void setupNeighbor(Parameter* param) coord = xlo - cutneigh - SMALL * xprd; mbinxlo = (int) (coord * bininvx); - if (coord < 0.0) { - mbinxlo = mbinxlo - 1; - } + if (coord < 0.0) { mbinxlo = mbinxlo - 1; } coord = xhi + cutneigh + SMALL * xprd; mbinxhi = (int) (coord * bininvx); coord = ylo - cutneigh - SMALL * yprd; mbinylo = (int) (coord * bininvy); - if (coord < 0.0) { - mbinylo = mbinylo - 1; - } + if (coord < 0.0) { mbinylo = mbinylo - 1; } 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; - } + if (coord < 0.0) { mbinzlo = mbinzlo - 1; } coord = zhi + cutneigh + SMALL * zprd; mbinzhi = (int) (coord * bininvz); @@ -152,20 +142,13 @@ void setupNeighbor(Parameter* param) 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)); - + if (stencil) { free(stencil); } + stencil = (int*) malloc((2 * nextz + 1) * (2 * nexty + 1) * (2 * nextx + 1) * sizeof(int)); nstencil = 0; int kstart = -nextz; @@ -173,28 +156,20 @@ void setupNeighbor(Parameter* param) 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; + stencil[nstencil++] = k * mbiny * mbinx + j * mbinx + i; } } } } mbins = mbinx * mbiny * mbinz; - - 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)); } -void buildNeighbor(Atom *atom, Neighbor *neighbor) -{ +void buildNeighbor(Atom *atom, Neighbor *neighbor) { int nall = atom->Nlocal + atom->Nghost; /* extend atom arrays if necessary */ @@ -231,8 +206,7 @@ void buildNeighbor(Atom *atom, Neighbor *neighbor) for(int m = 0; m < bincount[jbin]; m++) { int j = loc_bin[m]; - - if ( (j == i) || (neighbor->halfneigh && (j < i))) { + if((j == i) || (neighbor->half_neigh && (j < i))) { continue; } @@ -247,7 +221,7 @@ void buildNeighbor(Atom *atom, Neighbor *neighbor) #else const MD_FLOAT cutoff = cutneighsq; #endif - if( rsq <= cutoff ) { + if(rsq <= cutoff) { neighptr[n++] = j; } } @@ -274,8 +248,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, int k) { MD_FLOAT delx, dely, delz; if(i > 0) { @@ -305,8 +278,7 @@ MD_FLOAT bindist(int i, int j, int k) return (delx * delx + dely * dely + delz * delz); } -int coord2bin(MD_FLOAT xin, MD_FLOAT yin, MD_FLOAT zin) -{ +int coord2bin(MD_FLOAT xin, MD_FLOAT yin, MD_FLOAT zin) { int ix, iy, iz; if(xin >= xprd) { @@ -336,8 +308,7 @@ int coord2bin(MD_FLOAT xin, MD_FLOAT yin, MD_FLOAT zin) return (iz * mbiny * mbinx + iy * mbinx + ix + 1); } -void binatoms(Atom *atom) -{ +void binatoms(Atom *atom) { int nall = atom->Nlocal + atom->Nghost; int resize = 1; diff --git a/lammps/parameter.c b/lammps/parameter.c index 5bb9632..90fc3cf 100644 --- a/lammps/parameter.c +++ b/lammps/parameter.c @@ -52,7 +52,7 @@ void initParameter(Parameter *param) { param->reneigh_every = 20; param->x_out_every = 20; param->v_out_every = 5; - param->half_neigh = 1; + param->half_neigh = 0; param->proc_freq = 2.4; }