diff --git a/gromacs/atom.c b/gromacs/atom.c
index 1567f74..9e0e830 100644
--- a/gromacs/atom.c
+++ b/gromacs/atom.c
@@ -45,8 +45,9 @@ void initAtom(Atom *atom)
 {
     atom->x  = NULL; atom->y  = NULL; atom->z  = NULL;
     atom->vx = NULL; atom->vy = NULL; atom->vz = NULL;
-    atom->fx = NULL; atom->fy = NULL; atom->fz = NULL;
     atom->cl_x = NULL;
+    atom->cl_v = NULL;
+    atom->cl_f = NULL;
     atom->Natoms = 0;
     atom->Nlocal = 0;
     atom->Nghost = 0;
@@ -275,9 +276,6 @@ void growAtom(Atom *atom)
     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));
-    atom->fx = (MD_FLOAT*) reallocate(atom->fx, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT));
-    atom->fy = (MD_FLOAT*) reallocate(atom->fy, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT));
-    atom->fz = (MD_FLOAT*) reallocate(atom->fz, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT));
     atom->type = (int *) reallocate(atom->type, ALIGNMENT, atom->Nmax * sizeof(int), nold * sizeof(int));
 }
 
@@ -287,4 +285,6 @@ void growClusters(Atom *atom)
     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 * sizeof(MD_FLOAT), nold * CLUSTER_DIM_N * 3 * sizeof(MD_FLOAT));
+    atom->cl_f = (MD_FLOAT*) reallocate(atom->cl_f, ALIGNMENT, atom->Nclusters_max * CLUSTER_DIM_N * 3 * sizeof(MD_FLOAT), nold * CLUSTER_DIM_N * 3 * sizeof(MD_FLOAT));
+    atom->cl_v = (MD_FLOAT*) reallocate(atom->cl_v, ALIGNMENT, atom->Nclusters_max * CLUSTER_DIM_N * 3 * sizeof(MD_FLOAT), nold * CLUSTER_DIM_N * 3 * sizeof(MD_FLOAT));
 }
diff --git a/gromacs/force_eam.c b/gromacs/force_eam.c
index b686abf..1c1456b 100644
--- a/gromacs/force_eam.c
+++ b/gromacs/force_eam.c
@@ -33,6 +33,7 @@
 #include <util.h>
 
 double computeForceEam(Eam* eam, Parameter* param, Atom *atom, Neighbor *neighbor, Stats *stats) {
+    /*
     if(eam->nmax < atom->Nmax) {
         eam->nmax = atom->Nmax;
         if(eam->fp != NULL) { free(eam->fp); }
@@ -45,9 +46,11 @@ double computeForceEam(Eam* eam, Parameter* param, Atom *atom, Neighbor *neighbo
     MD_FLOAT* rhor_spline = eam->rhor_spline; MD_FLOAT* frho_spline = eam->frho_spline; MD_FLOAT* z2r_spline = eam->z2r_spline;
     MD_FLOAT rdr = eam->rdr; int nr = eam->nr; int nr_tot = eam->nr_tot; MD_FLOAT rdrho = eam->rdrho;
     int nrho = eam->nrho; int nrho_tot = eam->nrho_tot;
+    */
     double S = getTimeStamp();
 
     LIKWID_MARKER_START("force_eam_fp");
+    /*
     #pragma omp parallel for
     for(int i = 0; i < Nlocal; i++) {
         neighs = &neighbor->neighbors[i * neighbor->maxneighs];
@@ -207,6 +210,7 @@ double computeForceEam(Eam* eam, Parameter* param, Atom *atom, Neighbor *neighbo
         addStat(stats->total_force_iters, (numneighs + VECTOR_WIDTH - 1) / VECTOR_WIDTH);
     }
 
+    */
     LIKWID_MARKER_STOP("force_eam");
     double E = getTimeStamp();
     return E-S;
diff --git a/gromacs/force_lj.c b/gromacs/force_lj.c
index 73c0269..fd69675 100644
--- a/gromacs/force_lj.c
+++ b/gromacs/force_lj.c
@@ -31,69 +31,61 @@
 double computeForceLJ(Parameter *param, Atom *atom, Neighbor *neighbor, Stats *stats) {
     int Nlocal = atom->Nlocal;
     int* neighs;
-    MD_FLOAT* fx = atom->fx;
-    MD_FLOAT* fy = atom->fy;
-    MD_FLOAT* fz = atom->fz;
-#ifndef EXPLICIT_TYPES
     MD_FLOAT cutforcesq = param->cutforce * param->cutforce;
     MD_FLOAT sigma6 = param->sigma6;
     MD_FLOAT epsilon = param->epsilon;
-#endif
 
-    for(int i = 0; i < Nlocal; i++) {
-        fx[i] = 0.0;
-        fy[i] = 0.0;
-        fz[i] = 0.0;
+    for(int ci = 0; ci < atom->Nclusters_local; ci++) {
+        MD_FLOAT *fptr = cluster_force_ptr(ci);
+        for(int cii = 0; cii < atom->clusters[ci].natoms; cii++) {
+            cluster_x(fptr, cii) = 0.0;
+            cluster_y(fptr, cii) = 0.0;
+            cluster_z(fptr, cii) = 0.0;
+        }
     }
 
     double S = getTimeStamp();
     LIKWID_MARKER_START("force");
 
 #pragma omp parallel for
-    for(int i = 0; i < Nlocal; i++) {
-        neighs = &neighbor->neighbors[i * neighbor->maxneighs];
-        int numneighs = neighbor->numneigh[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;
-        MD_FLOAT fiz = 0;
-
-#ifdef EXPLICIT_TYPES
-        const int type_i = atom->type[i];
-#endif
-
+    for(int ci = 0; ci < atom->Nclusters_local; ci++) {
+        MD_FLOAT *ciptr = cluster_pos_ptr(ci);
+        MD_FLOAT *cifptr = cluster_force_ptr(ci);
+        neighs = &neighbor->neighbors[ci * neighbor->maxneighs];
+        int numneighs = neighbor->numneigh[ci];
 
         for(int k = 0; k < numneighs; k++) {
-            int j = neighs[k];
-            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;
+            int cj = neighs[k];
+            MD_FLOAT *cjptr = cluster_pos_ptr(cj);
+            for(int cii = 0; cii < CLUSTER_DIM_M; cii++) {
+                MD_FLOAT xtmp = cluster_x(ciptr, cii);
+                MD_FLOAT ytmp = cluster_y(ciptr, cii);
+                MD_FLOAT ztmp = cluster_z(ciptr, cii);
+                MD_FLOAT fix = 0;
+                MD_FLOAT fiy = 0;
+                MD_FLOAT fiz = 0;
 
-#ifdef EXPLICIT_TYPES
-            const int type_j = atom->type[j];
-            const int type_ij = type_i * atom->ntypes + type_j;
-            const MD_FLOAT cutforcesq = atom->cutforcesq[type_ij];
-            const MD_FLOAT sigma6 = atom->sigma6[type_ij];
-            const MD_FLOAT epsilon = atom->epsilon[type_ij];
-#endif
+                for(int cjj = 0; cjj < CLUSTER_DIM_N; cjj++) {
+                    MD_FLOAT delx = xtmp - cluster_x(cjptr, cjj);
+                    MD_FLOAT dely = ytmp - cluster_y(cjptr, cjj);
+                    MD_FLOAT delz = ztmp - cluster_z(cjptr, cjj);
+                    MD_FLOAT rsq = delx * delx + dely * dely + delz * delz;
+                    if(rsq < cutforcesq) {
+                        MD_FLOAT sr2 = 1.0 / rsq;
+                        MD_FLOAT sr6 = sr2 * sr2 * sr2 * sigma6;
+                        MD_FLOAT force = 48.0 * sr6 * (sr6 - 0.5) * sr2 * epsilon;
+                        fix += delx * force;
+                        fiy += dely * force;
+                        fiz += delz * force;
+                    }
+                }
 
-            if(rsq < cutforcesq) {
-                MD_FLOAT sr2 = 1.0 / rsq;
-                MD_FLOAT sr6 = sr2 * sr2 * sr2 * sigma6;
-                MD_FLOAT force = 48.0 * sr6 * (sr6 - 0.5) * sr2 * epsilon;
-                fix += delx * force;
-                fiy += dely * force;
-                fiz += delz * force;
+                cluster_x(cifptr, cii) += fix;
+                cluster_y(cifptr, cii) += fiy;
+                cluster_z(cifptr, cii) += fiz;
             }
         }
 
-        fx[i] += fix;
-        fy[i] += fiy;
-        fz[i] += fiz;
-
         addStat(stats->total_force_neighs, numneighs);
         addStat(stats->total_force_iters, (numneighs + VECTOR_WIDTH - 1) / VECTOR_WIDTH);
     }
diff --git a/gromacs/includes/atom.h b/gromacs/includes/atom.h
index 4b2a5bd..17065cf 100644
--- a/gromacs/includes/atom.h
+++ b/gromacs/includes/atom.h
@@ -42,8 +42,9 @@ typedef struct {
     int Nclusters, Nclusters_local, Nclusters_ghost, Nclusters_max;
     MD_FLOAT *x, *y, *z;
     MD_FLOAT *vx, *vy, *vz;
-    MD_FLOAT *fx, *fy, *fz;
     MD_FLOAT *cl_x;
+    MD_FLOAT *cl_v;
+    MD_FLOAT *cl_f;
     int *border_map;
     int *type;
     int ntypes;
@@ -61,7 +62,9 @@ extern int readAtom(Atom*, Parameter*);
 extern void growAtom(Atom*);
 extern void growClusters(Atom*);
 
-#define cluster_ptr(ci)     &(atom->cl_x[(ci) * CLUSTER_DIM_N * 3])
+#define cluster_pos_ptr(ci)       &(atom->cl_x[(ci) * CLUSTER_DIM_N * 3])
+#define cluster_velocity_ptr(ci)  &(atom->cl_v[(ci) * CLUSTER_DIM_N * 3])
+#define cluster_force_ptr(ci)     &(atom->cl_f[(ci) * CLUSTER_DIM_N * 3])
 
 #ifdef AOS
 #define POS_DATA_LAYOUT     "AoS"
diff --git a/gromacs/main.c b/gromacs/main.c
index ca18756..8a05002 100644
--- a/gromacs/main.c
+++ b/gromacs/main.c
@@ -129,30 +129,33 @@ double reneighbour(
     return E-S;
 }
 
-void initialIntegrate(Parameter *param, Atom *atom)
-{
-    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;
+void initialIntegrate(Parameter *param, Atom *atom) {
+    for(int ci = 0; ci < atom->Nclusters_local; ci++) {
+        MD_FLOAT *ciptr = cluster_pos_ptr(ci);
+        MD_FLOAT *civptr = cluster_velocity_ptr(ci);
+        MD_FLOAT *cifptr = cluster_force_ptr(ci);
 
-    for(int i = 0; i < atom->Nlocal; i++) {
-        vx[i] += param->dtforce * fx[i];
-        vy[i] += param->dtforce * fy[i];
-        vz[i] += param->dtforce * fz[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];
+        for(int cii = 0; cii < atom->clusters[ci].natoms; cii++) {
+            cluster_x(civptr, cii) += param->dtforce * cluster_x(cifptr, cii);
+            cluster_y(civptr, cii) += param->dtforce * cluster_y(cifptr, cii);
+            cluster_z(civptr, cii) += param->dtforce * cluster_z(cifptr, cii);
+            cluster_x(ciptr, cii) += param->dt * cluster_x(civptr, cii);
+            cluster_y(ciptr, cii) += param->dt * cluster_y(civptr, cii);
+            cluster_z(ciptr, cii) += param->dt * cluster_z(civptr, cii);
+        }
     }
 }
 
-void finalIntegrate(Parameter *param, Atom *atom)
-{
-    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;
+void finalIntegrate(Parameter *param, Atom *atom) {
+    for(int ci = 0; ci < atom->Nclusters_local; ci++) {
+        MD_FLOAT *civptr = cluster_velocity_ptr(ci);
+        MD_FLOAT *cifptr = cluster_force_ptr(ci);
 
-    for(int i = 0; i < atom->Nlocal; i++) {
-        vx[i] += param->dtforce * fx[i];
-        vy[i] += param->dtforce * fy[i];
-        vz[i] += param->dtforce * fz[i];
+        for(int cii = 0; cii < atom->clusters[ci].natoms; cii++) {
+            cluster_x(civptr, cii) += param->dtforce * cluster_x(cifptr, cii);
+            cluster_y(civptr, cii) += param->dtforce * cluster_y(cifptr, cii);
+            cluster_z(civptr, cii) += param->dtforce * cluster_z(cifptr, cii);
+        }
     }
 }
 
@@ -257,13 +260,11 @@ 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();
@@ -273,7 +274,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);
@@ -285,15 +286,13 @@ 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);
diff --git a/gromacs/neighbor.c b/gromacs/neighbor.c
index 8d84036..1ae533e 100644
--- a/gromacs/neighbor.c
+++ b/gromacs/neighbor.c
@@ -184,8 +184,8 @@ MD_FLOAT getBoundingBoxDistanceSq(Atom *atom, int ci, int cj) {
 }
 
 int atomDistanceInRange(Atom *atom, int ci, int cj, MD_FLOAT rsq) {
-    MD_FLOAT *ciptr = cluster_ptr(ci);
-    MD_FLOAT *cjptr = cluster_ptr(cj);
+    MD_FLOAT *ciptr = cluster_pos_ptr(ci);
+    MD_FLOAT *cjptr = cluster_pos_ptr(cj);
 
     for(int cii = 0; cii < atom->clusters[ci].natoms; cii++) {
         for(int cjj = 0; cjj < atom->clusters[cj].natoms; cjj++) {
@@ -422,7 +422,7 @@ void buildClusters(Atom *atom) {
                 growClusters(atom);
             }
 
-            MD_FLOAT *cptr = cluster_ptr(ci);
+            MD_FLOAT *cptr = cluster_pos_ptr(ci);
             MD_FLOAT bbminx = INFINITY, bbmaxx = -INFINITY;
             MD_FLOAT bbminy = INFINITY, bbmaxy = -INFINITY;
             MD_FLOAT bbminz = INFINITY, bbmaxz = -INFINITY;
@@ -493,7 +493,7 @@ void binClusters(Atom *atom) {
         }
 
         for(int ci = 0; ci < atom->Nclusters_ghost && !resize; ci++) {
-            MD_FLOAT *cptr = cluster_ptr(nlocal + ci);
+            MD_FLOAT *cptr = cluster_pos_ptr(nlocal + ci);
             MD_FLOAT xtmp, ytmp;
             int ix = -1, iy = -1;
 
@@ -539,7 +539,7 @@ void updateSingleAtoms(Atom *atom) {
     int Natom = 0;
 
     for(int ci = 0; ci < atom->Nclusters_local; ci++) {
-        MD_FLOAT *cptr = cluster_ptr(ci);
+        MD_FLOAT *cptr = cluster_pos_ptr(ci);
 
         for(int cii = 0; cii < atom->clusters[ci].natoms; cii++) {
             atom_x(Natom) = cluster_x(cptr, cii);
diff --git a/gromacs/pbc.c b/gromacs/pbc.c
index 1bef7c4..31d8031 100644
--- a/gromacs/pbc.c
+++ b/gromacs/pbc.c
@@ -51,8 +51,8 @@ void updatePbc(Atom *atom, Parameter *param) {
     MD_FLOAT zprd = param->zprd;
 
     for(int ci = 0; ci < atom->Nclusters_ghost; ci++) {
-        MD_FLOAT *cptr = cluster_ptr(nlocal + ci);
-        MD_FLOAT *bmap_cptr = cluster_ptr(border_map[ci]);
+        MD_FLOAT *cptr = cluster_pos_ptr(nlocal + ci);
+        MD_FLOAT *bmap_cptr = cluster_pos_ptr(border_map[ci]);
 
         for(int cii = 0; cii < atom->clusters[ci].natoms; cii++) {
             cluster_x(cptr, cii) = cluster_x(bmap_cptr, cii) + atom->PBCx[ci] * xprd;
diff --git a/gromacs/tracing.c b/gromacs/tracing.c
index d91a97a..0bf3bfc 100644
--- a/gromacs/tracing.c
+++ b/gromacs/tracing.c
@@ -30,7 +30,7 @@ void traceAddresses(Parameter *param, Atom *atom, Neighbor *neighbor, int timest
     INDEX_TRACER_INIT;
     int Nlocal = atom->Nlocal;
     int* neighs;
-    MD_FLOAT* fx = atom->fx; MD_FLOAT* fy = atom->fy; MD_FLOAT* fz = atom->fz;
+    //MD_FLOAT* fx = atom->fx; MD_FLOAT* fy = atom->fy; MD_FLOAT* fz = atom->fz;
 
     INDEX_TRACE_NATOMS(Nlocal, atom->Nghost, neighbor->maxneighs);
     for(int i = 0; i < Nlocal; i++) {
@@ -60,12 +60,14 @@ void traceAddresses(Parameter *param, Atom *atom, Neighbor *neighbor, int timest
             #endif
         }
 
+        /*
         MEM_TRACE(fx[i], 'R');
         MEM_TRACE(fx[i], 'W');
         MEM_TRACE(fy[i], 'R');
         MEM_TRACE(fy[i], 'W');
         MEM_TRACE(fz[i], 'R');
         MEM_TRACE(fz[i], 'W');
+        */
     }
 
     INDEX_TRACER_END;