Fix some segfaults and add function to update single atoms

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

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));
 }

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 buildClusters(Atom*);
-extern void binGhostClusters(Parameter*, Atom*);
+extern void binClusters(Atom*);
+extern void updateSingleAtoms(Atom*);
 #endif

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);
+    binClusters(atom);
-    buildNeighbor(param, atom, neighbor);
+    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(&param, &atom, &neighbor, n + 1);
 #endif
+    /*
     if(param.force_field == FF_EAM) {
         timer[FORCE] = computeForceEam(&eam, &param, &atom, &neighbor, &stats);
     } else {
         timer[FORCE] = computeForceLJ(&param, &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(&param, &atom);
+        //initialIntegrate(&param, &atom);
 
         if((n + 1) % param.every) {
             updatePbc(&atom, &param);
@@ -279,13 +285,15 @@ int main(int argc, char** argv)
         traceAddresses(&param, &atom, &neighbor, n + 1);
 #endif
 
+        /*
         if(param.force_field == FF_EAM) {
             timer[FORCE] += computeForceEam(&eam, &param, &atom, &neighbor, &stats);
         } else {
             timer[FORCE] += computeForceLJ(&param, &atom, &neighbor, &stats);
         }
+        */
 
-        finalIntegrate(&param, &atom);
+        //finalIntegrate(&param, &atom);
 
         if(!((n + 1) % param.nstat) && (n+1) < param.ntimes) {
             computeThermo(n + 1, &param, &atom);

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);
-    if(param->input_file != NULL) {
+    nbiny = (int)((yhi - ylo) / binsizey);
-        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) {
+void buildNeighbor(Atom *atom, Neighbor *neighbor) {
-    int nall = atom->Nlocal + atom->Nghost;
+    printf("buildNeighbor start\n");
+    int nall = atom->Nclusters_local + atom->Nclusters_ghost;
 
     /* extend atom arrays if necessary */
     if(nall > nmax) {
@@ -239,15 +239,17 @@ 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];
 
+                if(c > 0) {
                     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);
+                    } while(m + 1 < c && (ibb_zmin - jbb_zmax) * (ibb_zmin - jbb_zmax) > cutneighsq);
 
-                while(m < cluster_bincount[jbin] && jbb_zmax - ibb_zmax < cutneighsq) {
+                    while(m < c && jbb_zmax - ibb_zmax < cutneighsq) {
                         if((jbb_zmin - ibb_zmax) * (jbb_zmin - ibb_zmax) > cutneighsq) {
                             break;
                         }
@@ -255,19 +257,19 @@ void buildNeighbor(Parameter *param, Atom *atom, Neighbor *neighbor) {
                         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
-                            }
                                 neighptr[n++] = cj;
                             }
                         }
 
                         m++;
+                        if(m < c) {
                             cj = loc_bin[m];
                             jbb_zmin = atom->clusters[cj].bbminz;
                             jbb_zmax = atom->clusters[cj].bbmaxz;
                         }
                     }
+                }
+            }
 
             neighbor->numneigh[ci] = n;
 
@@ -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,10 +384,12 @@ 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) {
+void sortAtomsByZCoord(Atom *atom) {
+    for(int bin = 0; bin < mbins; bin++) {
         int c = bincount[bin];
         int *bin_ptr = &bins[bin * atoms_per_bin];
 
@@ -406,27 +412,28 @@ void sortBinAtomsByZCoord(Parameter *param, Atom *atom, int bin) {
             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 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(nclusters < clusters_per_bin) {
         for(int cl = 0; cl < nclusters; cl++) {
             const int ci = atom->Nclusters_local;
+
+            if(ci >= atom->Nclusters_max) {
+                growClusters(atom);
+            }
+
             MD_FLOAT *cptr = cluster_ptr(ci);
             MD_FLOAT bbminx = INFINITY, bbmaxx = -INFINITY;
             MD_FLOAT bbminy = INFINITY, bbmaxy = -INFINITY;
@@ -463,7 +470,6 @@ void buildClusters(Parameter *param, Atom *atom) {
                 ac++;
             }
 
-                    cluster_bins[bin * clusters_per_bin + cl] = ci;
             atom->clusters[ci].bin = bin;
             atom->clusters[ci].bbminx = bbminx;
             atom->clusters[ci].bbmaxx = bbmaxx;
@@ -473,25 +479,37 @@ void buildClusters(Parameter *param, Atom *atom) {
             atom->clusters[ci].bbmaxz = bbmaxz;
             atom->Nclusters_local++;
         }
+    }
 
-                cluster_bincount[bin] = nclusters;
+    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;
             }
         }
 
-        if(resize) {
+        printf("ghost\n");
-            free(cluster_bins);
+        for(int ci = 0; ci < atom->Nclusters_ghost && !resize; ci++) {
-            clusters_per_bin *= 2;
-            cluster_bins = (int*) malloc(mbins * clusters_per_bin * sizeof(int));
-        }
-    }
-}
-
-void binGhostClusters(Parameter *param, Atom *atom) {
-    int nlocal = atom->Nclusters_local;
-
-    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;
@@ -515,7 +533,40 @@ void binGhostClusters(Parameter *param, Atom *atom) {
 
             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 updateSingleAtoms(Atom *atom) {
+    int Natom = 0;
+
+    for(int ci = 0; ci < atom->Nclusters_local; ci++) {
+        MD_FLOAT *cptr = cluster_ptr(ci);
+
+        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++;
+        }
+    }
+
+    if(Natom != atom->Nlocal) {
+        fprintf(stderr, "updateSingleAtoms(): Number of atoms changed!\n");
     }
 }

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);
 }