diff --git a/src/force.c b/src/force.c
new file mode 100644
index 0000000..5dd56bf
--- /dev/null
+++ b/src/force.c
@@ -0,0 +1,89 @@
+/*
+ * =======================================================================================
+ *
+ *   Author:   Jan Eitzinger (je), jan.eitzinger@fau.de
+ *   Copyright (c) 2021 RRZE, University Erlangen-Nuremberg
+ *
+ *   This file is part of MD-Bench.
+ *
+ *   MD-Bench is free software: you can redistribute it and/or modify it
+ *   under the terms of the GNU Lesser General Public License as published
+ *   by the Free Software Foundation, either version 3 of the License, or
+ *   (at your option) any later version.
+ *
+ *   MD-Bench is distributed in the hope that it will be useful, but WITHOUT ANY
+ *   WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
+ *   PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more
+ *   details.
+ *
+ *   You should have received a copy of the GNU Lesser General Public License along
+ *   with MD-Bench.  If not, see <https://www.gnu.org/licenses/>.
+ * =======================================================================================
+ */
+#include <likwid-marker.h>
+
+#include <timing.h>
+#include <neighbor.h>
+#include <parameter.h>
+#include <atom.h>
+
+double computeForce(
+        Parameter *param,
+        Atom *atom,
+        Neighbor *neighbor)
+{
+    int Nlocal = atom->Nlocal;
+    int* neighs;
+    MD_FLOAT cutforcesq = param->cutforce * param->cutforce;
+    MD_FLOAT sigma6 = param->sigma6;
+    MD_FLOAT epsilon = param->epsilon;
+    MD_FLOAT* fx = atom->fx; MD_FLOAT* fy = atom->fy; MD_FLOAT* fz = atom->fz;
+    MD_FLOAT S, E;
+
+    S = getTimeStamp();
+    for(int i = 0; i < Nlocal; i++) {
+        fx[i] = 0.0;
+        fy[i] = 0.0;
+        fz[i] = 0.0;
+    }
+
+    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;
+
+        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;
+
+            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;
+            }
+        }
+
+        fx[i] += fix;
+        fy[i] += fiy;
+        fz[i] += fiz;
+    }
+    LIKWID_MARKER_STOP("force");
+    E = getTimeStamp();
+
+    return E-S;
+}
diff --git a/src/main.c b/src/main.c
index c064815..15d71fd 100644
--- a/src/main.c
+++ b/src/main.c
@@ -47,6 +47,7 @@ typedef enum {
     NUMTIMER
 } timertype;
 
+extern double computeForce( Parameter*, Atom*, Neighbor*);
 
 void init(Parameter *param)
 {
@@ -141,65 +142,6 @@ void finalIntegrate(Parameter *param, Atom *atom)
     }
 }
 
-double computeForce(Parameter *param, Atom *atom, Neighbor *neighbor)
-{
-    int Nlocal = atom->Nlocal;
-    int* neighs;
-    MD_FLOAT cutforcesq = param->cutforce * param->cutforce;
-    MD_FLOAT sigma6 = param->sigma6;
-    MD_FLOAT epsilon = param->epsilon;
-    MD_FLOAT* fx = atom->fx; MD_FLOAT* fy = atom->fy; MD_FLOAT* fz = atom->fz;
-    MD_FLOAT S, E;
-
-    S = getTimeStamp();
-    for(int i = 0; i < Nlocal; i++) {
-        fx[i] = 0.0;
-        fy[i] = 0.0;
-        fz[i] = 0.0;
-    }
-
-    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;
-
-        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;
-
-            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;
-            }
-        }
-
-        fx[i] += fix;
-        fy[i] += fiy;
-        fz[i] += fiz;
-    }
-    LIKWID_MARKER_STOP("force");
-    E = getTimeStamp();
-
-    return E-S;
-}
-
-
 void printAtomState(Atom *atom)
 {
     printf("Atom counts: Natoms=%d Nlocal=%d Nghost=%d Nmax=%d\n",