diff --git a/gromacs/force_lj.c b/gromacs/force_lj.c
index 46cdd94..7ff61cb 100644
--- a/gromacs/force_lj.c
+++ b/gromacs/force_lj.c
@@ -67,7 +67,7 @@ double computeForceLJ(Parameter *param, Atom *atom, Neighbor *neighbor, Stats *s
                 MD_FLOAT fiy = 0;
                 MD_FLOAT fiz = 0;
 
-                for(int cjj = 0; cjj < CLUSTER_DIM_N; cjj++) {
+                for(int cjj = 0; cjj < atom->clusters[cj].natoms; cjj++) {
                     if(ci != cj || cii != cjj) {
                         MD_FLOAT delx = xtmp - cluster_x(cjptr, cjj);
                         MD_FLOAT dely = ytmp - cluster_y(cjptr, cjj);
@@ -80,6 +80,11 @@ 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);
+                            }
+                            */
                         }
                     }
                 }
@@ -87,6 +92,10 @@ 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);
+                }*/
             }
         }
 
diff --git a/gromacs/includes/atom.h b/gromacs/includes/atom.h
index 17065cf..4b04fda 100644
--- a/gromacs/includes/atom.h
+++ b/gromacs/includes/atom.h
@@ -25,8 +25,8 @@
 #ifndef __ATOM_H_
 #define __ATOM_H_
 
-#define CLUSTER_DIM_N       4
 #define CLUSTER_DIM_M       4
+#define CLUSTER_DIM_N       4
 
 typedef struct {
     int bin;
diff --git a/gromacs/includes/vtk.h b/gromacs/includes/vtk.h
index 0f03a74..66b3b32 100644
--- a/gromacs/includes/vtk.h
+++ b/gromacs/includes/vtk.h
@@ -24,5 +24,6 @@
 
 #ifndef __VTK_H_
 #define __VTK_H_
-extern int write_atoms_to_vtk_file(const char* filename, Atom* atom, int timestep);
+extern int write_local_atoms_to_vtk_file(const char* filename, Atom* atom, int timestep);
+extern int write_ghost_atoms_to_vtk_file(const char* filename, Atom* atom, int timestep);
 #endif
diff --git a/gromacs/main.c b/gromacs/main.c
index b4c1bd9..3556c6a 100644
--- a/gromacs/main.c
+++ b/gromacs/main.c
@@ -260,7 +260,8 @@ int main(int argc, char** argv) {
     timer[TOTAL] = getTimeStamp();
 
     if(param.vtk_file != NULL) {
-        write_atoms_to_vtk_file(param.vtk_file, &atom, 0);
+        write_local_atoms_to_vtk_file(param.vtk_file, &atom, 0);
+        write_ghost_atoms_to_vtk_file(param.vtk_file, &atom, 0);
     }
 
     for(int n = 0; n < param.ntimes; n++) {
@@ -289,7 +290,8 @@ int main(int argc, char** argv) {
         }
 
         if(param.vtk_file != NULL) {
-            write_atoms_to_vtk_file(param.vtk_file, &atom, n + 1);
+            write_local_atoms_to_vtk_file(param.vtk_file, &atom, n + 1);
+            write_ghost_atoms_to_vtk_file(param.vtk_file, &atom, n + 1);
         }
     }
 
diff --git a/gromacs/neighbor.c b/gromacs/neighbor.c
index 1e65b33..ad7455a 100644
--- a/gromacs/neighbor.c
+++ b/gromacs/neighbor.c
@@ -289,7 +289,7 @@ void buildNeighbor(Atom *atom, Neighbor *neighbor) {
 
     /*
     printf("\ncutneighsq = %f, rbb_sq = %f\n", cutneighsq, rbb_sq);
-    for(int ci = 4; ci < 6; ci++) {
+    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);
@@ -377,7 +377,6 @@ void coord2bin2D(MD_FLOAT xin, MD_FLOAT yin, int *ix, int *iy) {
 
 void binAtoms(Atom *atom) {
     printf("binAtoms start\n");
-    int nall = atom->Nlocal + atom->Nghost;
     int resize = 1;
 
     while(resize > 0) {
@@ -387,7 +386,7 @@ void binAtoms(Atom *atom) {
             bincount[i] = 0;
         }
 
-        for(int i = 0; i < nall; i++) {
+        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) {
@@ -507,8 +506,15 @@ void buildClusters(Atom *atom) {
         }
     }
 
+    printf("buildClusters end\n");
+}
+
+void binClusters(Atom *atom) {
+    printf("binClusters start\n");
+
     /*
-    for(int ci = 4; ci < 9; ci++) {
+    fprintf(stdout, "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);
@@ -518,11 +524,6 @@ void buildClusters(Atom *atom) {
     }
     */
 
-    printf("buildClusters end\n");
-}
-
-void binClusters(Atom *atom) {
-    printf("binClusters start\n");
     const int nlocal = atom->Nclusters_local;
     int resize = 1;
 
diff --git a/gromacs/pbc.c b/gromacs/pbc.c
index 54c0bf8..3d93c3b 100644
--- a/gromacs/pbc.c
+++ b/gromacs/pbc.c
@@ -122,19 +122,18 @@ void updateAtomsPbc(Atom *atom, Parameter *param) {
  * defining ghost atoms around domain
  * only creates mapping and coordinate corrections
  * that are then enforced in updatePbc */
-#define ADDGHOST(dx,dy,dz)                                          \
-    Nghost++;                                                       \
-    border_map[Nghost] = ci;                                        \
-    atom->PBCx[Nghost] = dx;                                        \
-    atom->PBCy[Nghost] = dy;                                        \
-    atom->PBCz[Nghost] = dz;                                        \
-    copy_cluster_types(atom, atom->Nclusters_local + Nghost, ci)
-
-void copy_cluster_types(Atom *atom, int dest, int src) {
-    for(int cii = 0; cii < atom->clusters[src].natoms; cii++) {
-        atom->clusters[dest].type[cii] = atom->clusters[src].type[cii];
+#define ADDGHOST(dx,dy,dz);                                                     \
+    Nghost++;                                                                   \
+    const int g_atom_idx = atom->Nclusters_local + Nghost;                      \
+    border_map[Nghost] = ci;                                                    \
+    atom->PBCx[Nghost] = dx;                                                    \
+    atom->PBCy[Nghost] = dy;                                                    \
+    atom->PBCz[Nghost] = dz;                                                    \
+    atom->clusters[g_atom_idx].natoms = atom->clusters[ci].natoms;              \
+    Nghost_atoms += atom->clusters[g_atom_idx].natoms;                          \
+    for(int cii = 0; cii < atom->clusters[ci].natoms; cii++) {                  \
+        atom->clusters[g_atom_idx].type[cii] = atom->clusters[ci].type[cii];    \
     }
-}
 
 /* internal subroutines */
 void growPbc(Atom* atom) {
@@ -154,6 +153,7 @@ void setupPbc(Atom *atom, Parameter *param) {
     MD_FLOAT zprd = param->zprd;
     MD_FLOAT Cutneigh = param->cutneigh;
     int Nghost = -1;
+    int Nghost_atoms = 0;
 
     for(int ci = 0; ci < atom->Nclusters_local; ci++) {
         if (atom->Nclusters_local + Nghost + 7 >= atom->Nclusters_max) {
@@ -205,6 +205,7 @@ void setupPbc(Atom *atom, Parameter *param) {
     }
 
     // increase by one to make it the ghost atom count
+    atom->Nghost = Nghost_atoms;
     atom->Nclusters_ghost = Nghost + 1;
     atom->Nclusters = atom->Nclusters_local + Nghost + 1;
 
diff --git a/gromacs/vtk.c b/gromacs/vtk.c
index f339ccd..4b788ea 100644
--- a/gromacs/vtk.c
+++ b/gromacs/vtk.c
@@ -3,9 +3,9 @@
 
 #include <atom.h>
 
-int write_atoms_to_vtk_file(const char* filename, Atom* atom, int timestep) {
+int write_local_atoms_to_vtk_file(const char* filename, Atom* atom, int timestep) {
     char timestep_filename[128];
-    snprintf(timestep_filename, sizeof timestep_filename, "%s_%d.vtk", filename, timestep);
+    snprintf(timestep_filename, sizeof timestep_filename, "%s_local_%d.vtk", filename, timestep);
     FILE* fp = fopen(timestep_filename, "wb");
 
     if(fp == NULL) {
@@ -45,3 +45,46 @@ int write_atoms_to_vtk_file(const char* filename, Atom* atom, int timestep) {
     fclose(fp);
     return 0;
 }
+
+int write_ghost_atoms_to_vtk_file(const char* filename, Atom* atom, int timestep) {
+    char timestep_filename[128];
+    snprintf(timestep_filename, sizeof timestep_filename, "%s_ghost_%d.vtk", filename, timestep);
+    FILE* fp = fopen(timestep_filename, "wb");
+
+    if(fp == NULL) {
+        fprintf(stderr, "Could not open VTK file for writing!\n");
+        return -1;
+    }
+
+    fprintf(fp, "# vtk DataFile Version 2.0\n");
+    fprintf(fp, "Particle data\n");
+    fprintf(fp, "ASCII\n");
+    fprintf(fp, "DATASET UNSTRUCTURED_GRID\n");
+    fprintf(fp, "POINTS %d double\n", atom->Nghost);
+    for(int ci = atom->Nclusters_local; ci < atom->Nclusters_local + atom->Nclusters_ghost; ++ci) {
+        MD_FLOAT *cptr = cluster_pos_ptr(ci);
+        for(int cii = 0; cii < atom->clusters[ci].natoms; ++cii) {
+            fprintf(fp, "%.4f %.4f %.4f\n", cluster_x(cptr, cii), cluster_y(cptr, cii), cluster_z(cptr, cii));
+        }
+    }
+    fprintf(fp, "\n\n");
+    fprintf(fp, "CELLS %d %d\n", atom->Nghost, atom->Nghost * 2);
+    for(int i = 0; i < atom->Nghost; ++i) {
+        fprintf(fp, "1 %d\n", i);
+    }
+    fprintf(fp, "\n\n");
+    fprintf(fp, "CELL_TYPES %d\n", atom->Nghost);
+    for(int i = 0; i < atom->Nghost; ++i) {
+        fprintf(fp, "1\n");
+    }
+    fprintf(fp, "\n\n");
+    fprintf(fp, "POINT_DATA %d\n", atom->Nghost);
+    fprintf(fp, "SCALARS mass double\n");
+    fprintf(fp, "LOOKUP_TABLE default\n");
+    for(int i = 0; i < atom->Nghost; i++) {
+        fprintf(fp, "1.0\n");
+    }
+    fprintf(fp, "\n\n");
+    fclose(fp);
+    return 0;
+}