diff --git a/Makefile b/Makefile
index 5370032..90e3ce1 100644
--- a/Makefile
+++ b/Makefile
@@ -1,5 +1,3 @@
-TAG = CLANG
-
 #CONFIGURE BUILD SYSTEM
 TARGET	   = MDBench-$(TAG)
 BUILD_DIR  = ./$(TAG)
@@ -8,14 +6,23 @@ MAKE_DIR   = ./
 Q         ?= @
 
 #DO NOT EDIT BELOW
+include $(MAKE_DIR)/config.mk
 include $(MAKE_DIR)/include_$(TAG).mk
 INCLUDES  += -I./src/includes
 
+ifeq ($(strip $(DATA_LAYOUT)),AOS)
+DEFINES +=  -DAOS
+endif
+ifeq ($(strip $(DATA_TYPE)),SP)
+DEFINES +=  -DPRECISION=1
+else
+DEFINES +=  -DPRECISION=2
+endif
+
 VPATH     = $(SRC_DIR)
 ASM       = $(patsubst $(SRC_DIR)/%.c, $(BUILD_DIR)/%.s,$(wildcard $(SRC_DIR)/*.c))
 OBJ       = $(patsubst $(SRC_DIR)/%.c, $(BUILD_DIR)/%.o,$(wildcard $(SRC_DIR)/*.c))
-CPPFLAGS := $(CPPFLAGS) $(DEFINES) $(INCLUDES)
-
+CPPFLAGS := $(CPPFLAGS) $(DEFINES) $(OPTIONS) $(INCLUDES)
 
 ${TARGET}: $(BUILD_DIR) $(OBJ)
 	@echo "===>  LINKING  $(TARGET)"
diff --git a/README.md b/README.md
index ba47bcb..94e2875 100644
--- a/README.md
+++ b/README.md
@@ -4,9 +4,10 @@ A simple, sequential  C implementation of the [Mantevo miniMD](https://github.co
 
 ## Build
 
-1. Open the `Makefile` and edit the `TAG` value according to the tool chain used. Currently supported is GCC, CLANG (LLVM), and ICC (Intel).
-2. Open and adapt the compiler flags in `<include_<TOOLCHAIN>.mk`, e.g. in `include_ICC.mk` for the Intel tool chain.
-3. Build the binary calling `make`.
+1. Open `config.mk` and edit the `TAG` value according to the tool chain used. Currently supported is GCC, CLANG (LLVM), and ICC (Intel).
+2. Change `DATA_LAYOUT` and `DATA_TYPE` if desired in config.mk.
+3. Open and adapt the compiler flags in `<include_<TOOLCHAIN>.mk`, e.g. in `include_ICC.mk` for the Intel tool chain.
+4. Build the binary calling `make`.
 
 You can clean intermediate build results with `make clean`, and all build results with `make distclean`.
 You have to call `make clean` before `make` if you changed the build settings.
diff --git a/config.mk b/config.mk
new file mode 100644
index 0000000..9cf51d0
--- /dev/null
+++ b/config.mk
@@ -0,0 +1,7 @@
+# Supported: GCC, CLANG, ICC
+TAG ?= GCC
+DATA_TYPE ?= SP#SP or DP
+DATA_LAYOUT ?= SoA#AOS or SOA
+
+#Feature options
+OPTIONS +=  -DALIGNMENT=64
diff --git a/include_CLANG.mk b/include_CLANG.mk
index 68e0076..59fa579 100644
--- a/include_CLANG.mk
+++ b/include_CLANG.mk
@@ -1,5 +1,4 @@
 CC  = cc
-CXX = cc
 LINKER = $(CC)
 
 ANSI_CFLAGS  = -ansi
@@ -7,11 +6,9 @@ ANSI_CFLAGS += -std=c99
 ANSI_CFLAGS += -pedantic
 ANSI_CFLAGS += -Wextra
 
-CFLAGS   = -Ofast $(ANSI_CFLAGS) -Xpreprocessor -fopenmp #-g
+CFLAGS   = -Ofast $(ANSI_CFLAGS) -g #-Xpreprocessor -fopenmp -g
 ASFLAGS  = -masm=intel
-CXXFLAGS = $(CFLAGS)
-FCFLAGS  =
 LFLAGS   =
-DEFINES  = -D_GNU_SOURCE -DALIGNMENT=64 -DPRECISION=2
+DEFINES  = -D_GNU_SOURCE
 INCLUDES =
-LIBS     = -lomp
+LIBS     = -lm #-lomp
diff --git a/include_GCC.mk b/include_GCC.mk
index 83b5c5b..d29cca9 100644
--- a/include_GCC.mk
+++ b/include_GCC.mk
@@ -1,7 +1,5 @@
 CC  = gcc
-CXX = g++
-FC  = gfortran
-LINKER = $(CXX)
+LINKER = $(CC)
 
 ANSI_CFLAGS  = -ansi
 ANSI_CFLAGS += -std=c99
@@ -9,13 +7,9 @@ ANSI_CFLAGS += -pedantic
 ANSI_CFLAGS += -Wextra
 
 # CFLAGS   = -O0 -g  -std=c99 -fargument-noalias
-CFLAGS   = -O3 -march=znver1  -ffast-math -funroll-loops -fopenmp
-CXXFLAGS = $(CFLAGS)
+CFLAGS   = -O3 -march=znver1  -ffast-math -funroll-loops # -fopenmp
 ASFLAGS  =  -masm=intel
-FCFLAGS  =
 LFLAGS   =
 DEFINES  = -D_GNU_SOURCE
 INCLUDES =
-LIBS     =
-
-
+LIBS     = -lm
diff --git a/include_ICC.mk b/include_ICC.mk
index 25207af..ddb61ac 100644
--- a/include_ICC.mk
+++ b/include_ICC.mk
@@ -1,5 +1,4 @@
 CC  = icc
-CXX = icpc
 LINKER = $(CC)
 
 OPENMP  = #-qopenmp
@@ -11,10 +10,8 @@ PROFILE  = #-profile-functions -g  -pg
 #OPTS     = -fast -no-vec $(PROFILE)
 OPTS     = -fast -xHost $(PROFILE)
 CFLAGS   = $(PROFILE) -restrict $(OPENMP) $(OPTS)
-CXXFLAGS = $(CFLAGS)
 ASFLAGS  = -masm=intel
-FCFLAGS  =
 LFLAGS   = $(PROFILE) $(OPTS) $(OPENMP)
-DEFINES  = -D_GNU_SOURCE  -DALIGNMENT=64 # -DLIKWID_PERFMON   -DPRECISION=1
+DEFINES  = -D_GNU_SOURCE # -DALIGNMENT=64  -DLIKWID_PERFMON   -DPRECISION=1
 INCLUDES = #$(LIKWID_INC)
-LIBS     = #$(LIKWID_LIB) -llikwid
+LIBS     = -lm #$(LIKWID_LIB) -llikwid
diff --git a/src/atom.c b/src/atom.c
index 313ffdc..687de55 100644
--- a/src/atom.c
+++ b/src/atom.c
@@ -1,8 +1,10 @@
 /*
  * =======================================================================================
  *
- *   Author:   Jan Eitzinger (je), jan.eitzinger@fau.de
- *   Copyright (c) 2020 RRZE, University Erlangen-Nuremberg
+ *   Authors:   Jan Eitzinger (je), jan.eitzinger@fau.de
+ *              Rafael Ravedutti (rr), rafaelravedutti@gmail.com
+ *
+ *   Copyright (c) 2021 RRZE, University Erlangen-Nuremberg
  *
  *   This file is part of MD-Bench.
  *
@@ -111,9 +113,9 @@ void createAtom(Atom *atom, Parameter *param)
                     growAtom(atom);
                 }
 
-                atom->x[atom->Nlocal] = xtmp;
-                atom->y[atom->Nlocal] = ytmp;
-                atom->z[atom->Nlocal] = ztmp;
+                atom_x(atom->Nlocal) = xtmp;
+                atom_y(atom->Nlocal) = ytmp;
+                atom_z(atom->Nlocal) = ztmp;
                 atom->vx[atom->Nlocal] = vxtmp;
                 atom->vy[atom->Nlocal] = vytmp;
                 atom->vz[atom->Nlocal] = vztmp;
@@ -136,9 +138,13 @@ void growAtom(Atom *atom)
     int nold = atom->Nmax;
     atom->Nmax += DELTA;
 
+    #ifdef AOS
+    atom->x  = (MD_FLOAT*) reallocate(atom->x,  ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT) * 3, nold * sizeof(MD_FLOAT) * 3);
+    #else
     atom->x  = (MD_FLOAT*) reallocate(atom->x,  ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT));
     atom->y  = (MD_FLOAT*) reallocate(atom->y,  ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT));
     atom->z  = (MD_FLOAT*) reallocate(atom->z,  ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT));
+    #endif
     atom->vx = (MD_FLOAT*) reallocate(atom->vx, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT));
     atom->vy = (MD_FLOAT*) reallocate(atom->vy, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT));
     atom->vz = (MD_FLOAT*) reallocate(atom->vz, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT));
diff --git a/src/includes/atom.h b/src/includes/atom.h
index 97e0429..e090936 100644
--- a/src/includes/atom.h
+++ b/src/includes/atom.h
@@ -2,7 +2,7 @@
  * =======================================================================================
  *
  *   Author:   Jan Eitzinger (je), jan.eitzinger@fau.de
- *   Copyright (c) 2020 RRZE, University Erlangen-Nuremberg
+ *   Copyright (c) 2021 RRZE, University Erlangen-Nuremberg
  *
  *   This file is part of MD-Bench.
  *
@@ -35,4 +35,17 @@ typedef struct {
 extern void initAtom(Atom*);
 extern void createAtom(Atom*, Parameter*);
 extern void growAtom(Atom*);
+
+#ifdef AOS
+#define POS_DATA_LAYOUT     "AoS"
+#define atom_x(i)           atom->x[(i) * 3 + 0]
+#define atom_y(i)           atom->x[(i) * 3 + 1]
+#define atom_z(i)           atom->x[(i) * 3 + 2]
+#else
+#define POS_DATA_LAYOUT     "SoA"
+#define atom_x(i)           atom->x[i]
+#define atom_y(i)           atom->y[i]
+#define atom_z(i)           atom->z[i]
+#endif
+
 #endif
diff --git a/src/includes/parameter.h b/src/includes/parameter.h
index 87da5f8..4f6f914 100644
--- a/src/includes/parameter.h
+++ b/src/includes/parameter.h
@@ -23,9 +23,6 @@
 #ifndef __PARAMETER_H_
 #define __PARAMETER_H_
 
-#ifndef PRECISION
-#define PRECISION 2
-#endif
 #if PRECISION == 1
 #define MD_FLOAT float
 #else
diff --git a/src/main.c b/src/main.c
index 7fb68d1..c064815 100644
--- a/src/main.c
+++ b/src/main.c
@@ -116,7 +116,6 @@ double reneighbour(
 
 void initialIntegrate(Parameter *param, Atom *atom)
 {
-    MD_FLOAT* x = atom->x; MD_FLOAT* y = atom->y; MD_FLOAT* z = atom->z;
     MD_FLOAT* fx = atom->fx; MD_FLOAT* fy = atom->fy; MD_FLOAT* fz = atom->fz;
     MD_FLOAT* vx = atom->vx; MD_FLOAT* vy = atom->vy; MD_FLOAT* vz = atom->vz;
 
@@ -124,9 +123,9 @@ void initialIntegrate(Parameter *param, Atom *atom)
         vx[i] += param->dtforce * fx[i];
         vy[i] += param->dtforce * fy[i];
         vz[i] += param->dtforce * fz[i];
-        x[i]  += param->dt * vx[i];
-        y[i]  += param->dt * vy[i];
-        z[i]  += param->dt * vz[i];
+        atom_x(i) = atom_x(i) + param->dt * vx[i];
+        atom_y(i) = atom_y(i) + param->dt * vy[i];
+        atom_z(i) = atom_z(i) + param->dt * vz[i];
     }
 }
 
@@ -149,7 +148,6 @@ double computeForce(Parameter *param, Atom *atom, Neighbor *neighbor)
     MD_FLOAT cutforcesq = param->cutforce * param->cutforce;
     MD_FLOAT sigma6 = param->sigma6;
     MD_FLOAT epsilon = param->epsilon;
-    MD_FLOAT* x = atom->x; MD_FLOAT* y = atom->y; MD_FLOAT* z = atom->z;
     MD_FLOAT* fx = atom->fx; MD_FLOAT* fy = atom->fy; MD_FLOAT* fz = atom->fz;
     MD_FLOAT S, E;
 
@@ -166,9 +164,9 @@ double computeForce(Parameter *param, Atom *atom, Neighbor *neighbor)
     for(int i = 0; i < Nlocal; i++) {
         neighs = &neighbor->neighbors[i * neighbor->maxneighs];
         int numneighs = neighbor->numneigh[i];
-        MD_FLOAT xtmp = x[i];
-        MD_FLOAT ytmp = y[i];
-        MD_FLOAT ztmp = z[i];
+        MD_FLOAT xtmp = atom_x(i);
+        MD_FLOAT ytmp = atom_y(i);
+        MD_FLOAT ztmp = atom_z(i);
 
         MD_FLOAT fix = 0;
         MD_FLOAT fiy = 0;
@@ -176,9 +174,9 @@ double computeForce(Parameter *param, Atom *atom, Neighbor *neighbor)
 
         for(int k = 0; k < numneighs; k++) {
             int j = neighs[k];
-            MD_FLOAT delx = xtmp - x[j];
-            MD_FLOAT dely = ytmp - y[j];
-            MD_FLOAT delz = ztmp - z[j];
+            MD_FLOAT delx = xtmp - atom_x(j);
+            MD_FLOAT dely = ytmp - atom_y(j);
+            MD_FLOAT delz = ztmp - atom_z(j);
             MD_FLOAT rsq = delx * delx + dely * dely + delz * delz;
 
             if(rsq < cutforcesq) {
@@ -293,6 +291,7 @@ int main (int argc, char** argv)
     computeThermo(-1, &param, &atom);
 
     printf(HLINE);
+    printf("Data layout for positions: %s\n", POS_DATA_LAYOUT);
 #if PRECISION == 1
     printf("Using single precision floating point.\n");
 #else
diff --git a/src/neighbor.c b/src/neighbor.c
index 7460efb..05591ed 100644
--- a/src/neighbor.c
+++ b/src/neighbor.c
@@ -2,7 +2,7 @@
  * =======================================================================================
  *
  *   Author:   Jan Eitzinger (je), jan.eitzinger@fau.de
- *   Copyright (c) 2020 RRZE, University Erlangen-Nuremberg
+ *   Copyright (c) 2021 RRZE, University Erlangen-Nuremberg
  *
  *   This file is part of MD-Bench.
  *
@@ -184,9 +184,6 @@ void buildNeighbor(Atom *atom, Neighbor *neighbor)
     /* bin local & ghost atoms */
     binatoms(atom);
     int resize = 1;
-    MD_FLOAT* x = atom->x;
-    MD_FLOAT* y = atom->y;
-    MD_FLOAT* z = atom->z;
 
     /* loop over each atom, storing neighbors */
     while(resize) {
@@ -196,9 +193,9 @@ void buildNeighbor(Atom *atom, Neighbor *neighbor)
         for(int i = 0; i < atom->Nlocal; i++) {
             int* neighptr = &(neighbor->neighbors[i * neighbor->maxneighs]);
             int n = 0;
-            MD_FLOAT xtmp = x[i];
-            MD_FLOAT ytmp = y[i];
-            MD_FLOAT ztmp = z[i];
+            MD_FLOAT xtmp = atom_x(i);
+            MD_FLOAT ytmp = atom_y(i);
+            MD_FLOAT ztmp = atom_z(i);
             int ibin = coord2bin(xtmp, ytmp, ztmp);
 
             for(int k = 0; k < nstencil; k++) {
@@ -212,9 +209,9 @@ void buildNeighbor(Atom *atom, Neighbor *neighbor)
                         continue;
                     }
 
-                    MD_FLOAT delx = xtmp - x[j];
-                    MD_FLOAT dely = ytmp - y[j];
-                    MD_FLOAT delz = ztmp - z[j];
+                    MD_FLOAT delx = xtmp - atom_x(j);
+                    MD_FLOAT dely = ytmp - atom_y(j);
+                    MD_FLOAT delz = ztmp - atom_z(j);
                     MD_FLOAT rsq = delx * delx + dely * dely + delz * delz;
 
                     if( rsq <= cutneighsq ) {
@@ -309,9 +306,6 @@ int coord2bin(MD_FLOAT xin, MD_FLOAT yin, MD_FLOAT zin)
 void binatoms(Atom *atom)
 {
     int nall = atom->Nlocal + atom->Nghost;
-    MD_FLOAT* x = atom->x;
-    MD_FLOAT* y = atom->y;
-    MD_FLOAT* z = atom->z;
     int resize = 1;
 
     while(resize > 0) {
@@ -322,7 +316,7 @@ void binatoms(Atom *atom)
         }
 
         for(int i = 0; i < nall; i++) {
-            int ibin = coord2bin(x[i], y[i], z[i]);
+            int ibin = coord2bin(atom_x(i), atom_y(i), atom_z(i));
 
             if(bincount[ibin] < atoms_per_bin) {
                 int ac = bincount[ibin]++;
diff --git a/src/pbc.c b/src/pbc.c
index 54ef7f8..6428669 100644
--- a/src/pbc.c
+++ b/src/pbc.c
@@ -48,17 +48,14 @@ void initPbc()
 void updatePbc(Atom *atom, Parameter *param)
 {
     int nlocal = atom->Nlocal;
-    MD_FLOAT* x = atom->x;
-    MD_FLOAT* y = atom->y;
-    MD_FLOAT* z = atom->z;
     MD_FLOAT xprd = param->xprd;
     MD_FLOAT yprd = param->yprd;
     MD_FLOAT zprd = param->zprd;
 
     for(int i = 0; i < atom->Nghost; i++) {
-        x[nlocal + i] = x[BorderMap[i]] + PBCx[i] * xprd;
-        y[nlocal + i] = y[BorderMap[i]] + PBCy[i] * yprd;
-        z[nlocal + i] = z[BorderMap[i]] + PBCz[i] * zprd;
+        atom_x(nlocal + i) = atom_x(BorderMap[i]) + PBCx[i] * xprd;
+        atom_y(nlocal + i) = atom_y(BorderMap[i]) + PBCy[i] * yprd;
+        atom_z(nlocal + i) = atom_z(BorderMap[i]) + PBCz[i] * zprd;
     }
 }
 
@@ -66,31 +63,28 @@ void updatePbc(Atom *atom, Parameter *param)
  * to periodic boundary conditions */
 void updateAtomsPbc(Atom *atom, Parameter *param)
 {
-    MD_FLOAT* x = atom->x;
-    MD_FLOAT* y = atom->y;
-    MD_FLOAT* z = atom->z;
     MD_FLOAT xprd = param->xprd;
     MD_FLOAT yprd = param->yprd;
     MD_FLOAT zprd = param->zprd;
 
     for(int i = 0; i < atom->Nlocal; i++) {
 
-        if(x[i] < 0.0) {
-            x[i] += xprd;
-        } else if(x[i] >= xprd) {
-            x[i] -= xprd;
+        if(atom_x(i) < 0.0) {
+            atom_x(i) += xprd;
+        } else if(atom_x(i) >= xprd) {
+            atom_x(i) -= xprd;
         }
 
-        if(y[i] < 0.0) {
-            y[i] += yprd;
-        } else if(y[i] >= yprd) {
-            y[i] -= yprd;
+        if(atom_y(i) < 0.0) {
+            atom_y(i) += yprd;
+        } else if(atom_y(i) >= yprd) {
+            atom_y(i) -= yprd;
         }
 
-        if(z[i] < 0.0) {
-            z[i] += zprd;
-        } else if(z[i] >= zprd) {
-            z[i] -= zprd;
+        if(atom_z(i) < 0.0) {
+            atom_z(i) += zprd;
+        } else if(atom_z(i) >= zprd) {
+            atom_z(i) -= zprd;
         }
     }
 }
@@ -99,10 +93,14 @@ 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++; BorderMap[Nghost] = i; PBCx[Nghost] = dx; PBCy[Nghost] = dy; PBCz[Nghost] = dz;
+#define ADDGHOST(dx,dy,dz) \
+    Nghost++;              \
+    BorderMap[Nghost] = i; \
+    PBCx[Nghost] = dx;     \
+    PBCy[Nghost] = dy;     \
+    PBCz[Nghost] = dz;
 void setupPbc(Atom *atom, Parameter *param)
 {
-    MD_FLOAT* x = atom->x; MD_FLOAT* y = atom->y; MD_FLOAT* z = atom->z;
     MD_FLOAT xprd = param->xprd;
     MD_FLOAT yprd = param->yprd;
     MD_FLOAT zprd = param->zprd;
@@ -113,42 +111,45 @@ void setupPbc(Atom *atom, Parameter *param)
 
         if (atom->Nlocal + Nghost + 7 >= atom->Nmax) {
             growAtom(atom);
-            x = atom->x;  y = atom->y;  z = atom->z;
         }
         if (Nghost + 7 >= NmaxGhost) {
             growPbc();
         }
 
+        MD_FLOAT x = atom_x(i);
+        MD_FLOAT y = atom_y(i);
+        MD_FLOAT z = atom_z(i);
+
         /* Setup ghost atoms */
         /* 6 planes */
-        if (x[i] < Cutneigh)         { ADDGHOST(+1,0,0); }
-        if (x[i] >= (xprd-Cutneigh)) { ADDGHOST(-1,0,0); }
-        if (y[i] < Cutneigh)         { ADDGHOST(0,+1,0); }
-        if (y[i] >= (yprd-Cutneigh)) { ADDGHOST(0,-1,0); }
-        if (z[i] < Cutneigh)         { ADDGHOST(0,0,+1); }
-        if (z[i] >= (zprd-Cutneigh)) { ADDGHOST(0,0,-1); }
+        if (x < Cutneigh)         { ADDGHOST(+1,0,0); }
+        if (x >= (xprd-Cutneigh)) { ADDGHOST(-1,0,0); }
+        if (y < Cutneigh)         { ADDGHOST(0,+1,0); }
+        if (y >= (yprd-Cutneigh)) { ADDGHOST(0,-1,0); }
+        if (z < Cutneigh)         { ADDGHOST(0,0,+1); }
+        if (z >= (zprd-Cutneigh)) { ADDGHOST(0,0,-1); }
         /* 8 corners */
-        if (x[i] < Cutneigh         && y[i] < Cutneigh         && z[i] < Cutneigh)         { ADDGHOST(+1,+1,+1); }
-        if (x[i] < Cutneigh         && y[i] >= (yprd-Cutneigh) && z[i] < Cutneigh)         { ADDGHOST(+1,-1,+1); }
-        if (x[i] < Cutneigh         && y[i] >= Cutneigh        && z[i] >= (zprd-Cutneigh)) { ADDGHOST(+1,+1,-1); }
-        if (x[i] < Cutneigh         && y[i] >= (yprd-Cutneigh) && z[i] >= (zprd-Cutneigh)) { ADDGHOST(+1,-1,-1); }
-        if (x[i] >= (xprd-Cutneigh) && y[i] < Cutneigh         && z[i] < Cutneigh)         { ADDGHOST(-1,+1,+1); }
-        if (x[i] >= (xprd-Cutneigh) && y[i] >= (yprd-Cutneigh) && z[i] < Cutneigh)         { ADDGHOST(-1,-1,+1); }
-        if (x[i] >= (xprd-Cutneigh) && y[i] < Cutneigh         && z[i] >= (zprd-Cutneigh)) { ADDGHOST(-1,+1,-1); }
-        if (x[i] >= (xprd-Cutneigh) && y[i] >= (yprd-Cutneigh) && z[i] >= (zprd-Cutneigh)) { ADDGHOST(-1,-1,-1); }
+        if (x < Cutneigh         && y < Cutneigh         && z < Cutneigh)         { ADDGHOST(+1,+1,+1); }
+        if (x < Cutneigh         && y >= (yprd-Cutneigh) && z < Cutneigh)         { ADDGHOST(+1,-1,+1); }
+        if (x < Cutneigh         && y >= Cutneigh        && z >= (zprd-Cutneigh)) { ADDGHOST(+1,+1,-1); }
+        if (x < Cutneigh         && y >= (yprd-Cutneigh) && z >= (zprd-Cutneigh)) { ADDGHOST(+1,-1,-1); }
+        if (x >= (xprd-Cutneigh) && y < Cutneigh         && z < Cutneigh)         { ADDGHOST(-1,+1,+1); }
+        if (x >= (xprd-Cutneigh) && y >= (yprd-Cutneigh) && z < Cutneigh)         { ADDGHOST(-1,-1,+1); }
+        if (x >= (xprd-Cutneigh) && y < Cutneigh         && z >= (zprd-Cutneigh)) { ADDGHOST(-1,+1,-1); }
+        if (x >= (xprd-Cutneigh) && y >= (yprd-Cutneigh) && z >= (zprd-Cutneigh)) { ADDGHOST(-1,-1,-1); }
         /* 12 edges */
-        if (x[i] < Cutneigh         && z[i] < Cutneigh)         { ADDGHOST(+1,0,+1); }
-        if (x[i] < Cutneigh         && z[i] >= (zprd-Cutneigh)) { ADDGHOST(+1,0,-1); }
-        if (x[i] >= (xprd-Cutneigh) && z[i] < Cutneigh)         { ADDGHOST(-1,0,+1); }
-        if (x[i] >= (xprd-Cutneigh) && z[i] >= (zprd-Cutneigh)) { ADDGHOST(-1,0,-1); }
-        if (y[i] < Cutneigh         && z[i] < Cutneigh)         { ADDGHOST(0,+1,+1); }
-        if (y[i] < Cutneigh         && z[i] >= (zprd-Cutneigh)) { ADDGHOST(0,+1,-1); }
-        if (y[i] >= (yprd-Cutneigh) && z[i] < Cutneigh)         { ADDGHOST(0,-1,+1); }
-        if (y[i] >= (yprd-Cutneigh) && z[i] >= (zprd-Cutneigh)) { ADDGHOST(0,-1,-1); }
-        if (y[i] < Cutneigh         && x[i] < Cutneigh)         { ADDGHOST(+1,+1,0); }
-        if (y[i] < Cutneigh         && x[i] >= (xprd-Cutneigh)) { ADDGHOST(-1,+1,0); }
-        if (y[i] >= (yprd-Cutneigh) && x[i] < Cutneigh)         { ADDGHOST(+1,-1,0); }
-        if (y[i] >= (yprd-Cutneigh) && x[i] >= (xprd-Cutneigh)) { ADDGHOST(-1,-1,0); }
+        if (x < Cutneigh         && z < Cutneigh)         { ADDGHOST(+1,0,+1); }
+        if (x < Cutneigh         && z >= (zprd-Cutneigh)) { ADDGHOST(+1,0,-1); }
+        if (x >= (xprd-Cutneigh) && z < Cutneigh)         { ADDGHOST(-1,0,+1); }
+        if (x >= (xprd-Cutneigh) && z >= (zprd-Cutneigh)) { ADDGHOST(-1,0,-1); }
+        if (y < Cutneigh         && z < Cutneigh)         { ADDGHOST(0,+1,+1); }
+        if (y < Cutneigh         && z >= (zprd-Cutneigh)) { ADDGHOST(0,+1,-1); }
+        if (y >= (yprd-Cutneigh) && z < Cutneigh)         { ADDGHOST(0,-1,+1); }
+        if (y >= (yprd-Cutneigh) && z >= (zprd-Cutneigh)) { ADDGHOST(0,-1,-1); }
+        if (y < Cutneigh         && x < Cutneigh)         { ADDGHOST(+1,+1,0); }
+        if (y < Cutneigh         && x >= (xprd-Cutneigh)) { ADDGHOST(-1,+1,0); }
+        if (y >= (yprd-Cutneigh) && x < Cutneigh)         { ADDGHOST(+1,-1,0); }
+        if (y >= (yprd-Cutneigh) && x >= (xprd-Cutneigh)) { ADDGHOST(-1,-1,0); }
     }
     // increase by one to make it the ghost atom count
     atom->Nghost = Nghost + 1;