diff --git a/Makefile b/Makefile
index 5370032..7727bfe 100644
--- a/Makefile
+++ b/Makefile
@@ -5,6 +5,7 @@ TARGET	   = MDBench-$(TAG)
 BUILD_DIR  = ./$(TAG)
 SRC_DIR    = ./src
 MAKE_DIR   = ./
+FLAGS      = #-DAOS
 Q         ?= @
 
 #DO NOT EDIT BELOW
@@ -14,7 +15,7 @@ INCLUDES  += -I./src/includes
 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) $(INCLUDES) ${FLAGS}
 
 
 ${TARGET}: $(BUILD_DIR) $(OBJ)
diff --git a/src/atom.c b/src/atom.c
index 313ffdc..f70bc35 100644
--- a/src/atom.c
+++ b/src/atom.c
@@ -111,9 +111,9 @@ void createAtom(Atom *atom, Parameter *param)
                     growAtom(atom);
                 }
 
-                atom->x[atom->Nlocal] = xtmp;
-                atom->y[atom->Nlocal] = ytmp;
-                atom->z[atom->Nlocal] = ztmp;
+                set_atom_x(atom, atom->Nlocal, xtmp);
+                set_atom_y(atom, atom->Nlocal, ytmp);
+                set_atom_z(atom, atom->Nlocal, ztmp);
                 atom->vx[atom->Nlocal] = vxtmp;
                 atom->vy[atom->Nlocal] = vytmp;
                 atom->vz[atom->Nlocal] = vztmp;
@@ -136,9 +136,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..016896f 100644
--- a/src/includes/atom.h
+++ b/src/includes/atom.h
@@ -35,4 +35,23 @@ typedef struct {
 extern void initAtom(Atom*);
 extern void createAtom(Atom*, Parameter*);
 extern void growAtom(Atom*);
+
+#ifdef AOS
+inline const char *posDataLayout() { return "AoS"; }
+inline void set_atom_x(Atom *atom, int i, MD_FLOAT x) { atom->x[i * 3 + 0] = x; }
+inline void set_atom_y(Atom *atom, int i, MD_FLOAT y) { atom->x[i * 3 + 1] = y; }
+inline void set_atom_z(Atom *atom, int i, MD_FLOAT z) { atom->x[i * 3 + 2] = z; }
+inline MD_FLOAT get_atom_x(Atom *atom, int i) { return atom->x[i * 3 + 0]; }
+inline MD_FLOAT get_atom_y(Atom *atom, int i) { return atom->x[i * 3 + 1]; }
+inline MD_FLOAT get_atom_z(Atom *atom, int i) { return atom->x[i * 3 + 2]; }
+#else
+inline const char *posDataLayout() { return "SoA"; }
+inline void set_atom_x(Atom *atom, int i, MD_FLOAT x) { atom->x[i] = x; }
+inline void set_atom_y(Atom *atom, int i, MD_FLOAT y) { atom->y[i] = y; }
+inline void set_atom_z(Atom *atom, int i, MD_FLOAT z) { atom->z[i] = z; }
+inline MD_FLOAT get_atom_x(Atom *atom, int i) { return atom->x[i]; }
+inline MD_FLOAT get_atom_y(Atom *atom, int i) { return atom->y[i]; }
+inline MD_FLOAT get_atom_z(Atom *atom, int i) { return atom->z[i]; }
+#endif
+
 #endif
diff --git a/src/main.c b/src/main.c
index 7fb68d1..8e507b7 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];
+        set_atom_x(atom, i, get_atom_x(atom, i) + param->dt * vx[i]);
+        set_atom_y(atom, i, get_atom_y(atom, i) + param->dt * vy[i]);
+        set_atom_z(atom, i, get_atom_z(atom, 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 = get_atom_x(atom, i);
+        MD_FLOAT ytmp = get_atom_y(atom, i);
+        MD_FLOAT ztmp = get_atom_z(atom, 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 - get_atom_x(atom, j);
+            MD_FLOAT dely = ytmp - get_atom_y(atom, j);
+            MD_FLOAT delz = ztmp - get_atom_z(atom, 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", posDataLayout());
 #if PRECISION == 1
     printf("Using single precision floating point.\n");
 #else
diff --git a/src/neighbor.c b/src/neighbor.c
index 7460efb..fc80909 100644
--- a/src/neighbor.c
+++ b/src/neighbor.c
@@ -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 = get_atom_x(atom, i);
+            MD_FLOAT ytmp = get_atom_y(atom, i);
+            MD_FLOAT ztmp = get_atom_z(atom, 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 - get_atom_x(atom, j);
+                    MD_FLOAT dely = ytmp - get_atom_y(atom, j);
+                    MD_FLOAT delz = ztmp - get_atom_z(atom, 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(get_atom_x(atom, i), get_atom_y(atom, i), get_atom_z(atom, i));
 
             if(bincount[ibin] < atoms_per_bin) {
                 int ac = bincount[ibin]++;
diff --git a/src/pbc.c b/src/pbc.c
index 54ef7f8..7631e16 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;
+        set_atom_x(atom, nlocal + i, get_atom_x(atom, BorderMap[i]) + PBCx[i] * xprd);
+        set_atom_y(atom, nlocal + i, get_atom_y(atom, BorderMap[i]) + PBCy[i] * yprd);
+        set_atom_z(atom, nlocal + i, get_atom_z(atom, BorderMap[i]) + PBCz[i] * zprd);
     }
 }
 
@@ -66,31 +63,31 @@ 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++) {
+        MD_FLOAT x = get_atom_x(atom, i);
+        MD_FLOAT y = get_atom_y(atom, i);
+        MD_FLOAT z = get_atom_z(atom, i);
 
-        if(x[i] < 0.0) {
-            x[i] += xprd;
-        } else if(x[i] >= xprd) {
-            x[i] -= xprd;
+        if(x < 0.0) {
+            set_atom_x(atom, i, x + xprd);
+        } else if(x >= xprd) {
+            set_atom_x(atom, i, x - xprd);
         }
 
-        if(y[i] < 0.0) {
-            y[i] += yprd;
-        } else if(y[i] >= yprd) {
-            y[i] -= yprd;
+        if(y < 0.0) {
+            set_atom_y(atom, i, y + yprd);
+        } else if(y >= yprd) {
+            set_atom_y(atom, i, y - yprd);
         }
 
-        if(z[i] < 0.0) {
-            z[i] += zprd;
-        } else if(z[i] >= zprd) {
-            z[i] -= zprd;
+        if(z < 0.0) {
+            set_atom_z(atom, i, z + zprd);
+        } else if(z >= zprd) {
+            set_atom_z(atom, i, z - zprd);
         }
     }
 }
@@ -102,7 +99,6 @@ void updateAtomsPbc(Atom *atom, Parameter *param)
 #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 +109,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 = get_atom_x(atom, i);
+        MD_FLOAT y = get_atom_y(atom, i);
+        MD_FLOAT z = get_atom_z(atom, 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;