Add prunning kernel

Signed-off-by: Rafael Ravedutti <rafaelravedutti@gmail.com>

gromacs/includes/neighbor.h

@@ -37,6 +37,7 @@ extern void initNeighbor(Neighbor*, Parameter*);
 extern void setupNeighbor(Parameter*, Atom*);
 extern void binatoms(Atom*);
 extern void buildNeighbor(Atom*, Neighbor*);
+extern void pruneNeighbor(Parameter*, Atom*, Neighbor*);
 extern void sortAtom(Atom*);
 extern void buildClusters(Atom*);
 extern void binClusters(Atom*);

gromacs/includes/parameter.h

@@ -42,6 +42,7 @@ typedef struct {
     int ntimes;
     int nstat;
     int every;
+    int prune_every;
     MD_FLOAT dt;
     MD_FLOAT dtforce;
     MD_FLOAT cutforce;

gromacs/main.c

@@ -75,6 +75,7 @@ void init(Parameter *param) {
     param->mass = 1.0;
     param->dtforce = 0.5 * param->dt;
     param->every = 20;
+    param->prune_every = 1000;
     param->proc_freq = 2.4;
 }
 
@@ -278,6 +279,10 @@ int main(int argc, char** argv) {
         initialIntegrate(&param, &atom);
 
         if((n + 1) % param.every) {
+            if(!((n + 1) % param.prune_every)) {
+                pruneNeighbor(&param, &atom, &neighbor);
+            }
+
             updatePbc(&atom, &param, 0);
         } else {
             timer[NEIGH] += reneighbour(&param, &atom, &neighbor);

gromacs/neighbor.c

@@ -381,6 +381,47 @@ void buildNeighbor(Atom *atom, Neighbor *neighbor) {
     DEBUG_MESSAGE("buildNeighbor end\n");
 }
 
+void pruneNeighbor(Parameter *param, Atom *atom, Neighbor *neighbor) {
+    DEBUG_MESSAGE("pruneNeighbor start\n");
+    int nall = atom->Nclusters_local + atom->Nclusters_ghost;
+    //MD_FLOAT cutsq = param->cutforce * param->cutforce;
+    MD_FLOAT cutsq = cutneighsq;
+
+    for(int ci = 0; ci < atom->Nclusters_local; ci++) {
+        int *neighs = &neighbor->neighbors[ci * neighbor->maxneighs];
+        int numneighs = neighbor->numneigh[ci];
+        int k = 0;
+
+        // Remove dummy clusters if necessary
+        if(CLUSTER_DIM_N > CLUSTER_DIM_M) {
+            while(neighs[numneighs - 1] == nall - 1) {
+                numneighs--;
+            }
+        }
+
+        while(k < numneighs) {
+            int cj = neighs[k];
+            if(atomDistanceInRange(atom, ci, cj, cutsq)) {
+                k++;
+            } else {
+                numneighs--;
+                neighs[k] = neighs[numneighs];
+            }
+        }
+
+        // Readd dummy clusters if necessary
+        if(CLUSTER_DIM_N > CLUSTER_DIM_M) {
+            while(numneighs % (CLUSTER_DIM_N / CLUSTER_DIM_M)) {
+                neighs[numneighs++] = nall - 1; // Last cluster is always a dummy cluster
+            }
+        }
+
+        neighbor->numneigh[ci] = numneighs;
+    }
+
+    DEBUG_MESSAGE("pruneNeighbor end\n");
+}
+
 /* internal subroutines */
 MD_FLOAT bindist(int i, int j) {
     MD_FLOAT delx, dely, delz;