diff --git a/Makefile b/Makefile
index 2b480c0..da24a29 100644
--- a/Makefile
+++ b/Makefile
@@ -33,9 +33,6 @@ ifneq ($(NEIGHBORS_LOOP_RUNS),)
     DEFINES += -DNEIGHBORS_LOOP_RUNS=$(NEIGHBORS_LOOP_RUNS)
 endif
 
-ifeq ($(strip $(PRINT_STATS)),true)
-    DEFINES += -DPRINT_STATS
-endif
 ifeq ($(strip $(EXPLICIT_TYPES)),true)
     DEFINES += -DEXPLICIT_TYPES
 endif
diff --git a/config.mk b/config.mk
index 312c19d..0765626 100644
--- a/config.mk
+++ b/config.mk
@@ -9,8 +9,6 @@ DATA_LAYOUT ?= AOS
 # Assembly syntax to generate (ATT/INTEL)
 ASM_SYNTAX ?= ATT
 
-# Output detailed statistics
-PRINT_STATS ?= true
 # Number of times to run the atoms loop on stubbed variant
 ATOMS_LOOP_RUNS ?= 1
 # Number of times to run the neighbors loop on stubbed variant
diff --git a/src/atom.c b/src/atom.c
index 4e28d7a..9964a14 100644
--- a/src/atom.c
+++ b/src/atom.c
@@ -41,14 +41,12 @@ void initAtom(Atom *atom)
     atom->Nlocal = 0;
     atom->Nghost = 0;
     atom->Nmax   = 0;
-    #ifdef EXPLICIT_TYPES
     atom->type = NULL;
     atom->ntypes = 0;
     atom->epsilon = NULL;
     atom->sigma6 = NULL;
     atom->cutforcesq = NULL;
     atom->cutneighsq = NULL;
-    #endif
 }
 
 void createAtom(Atom *atom, Parameter *param)
@@ -58,8 +56,6 @@ void createAtom(Atom *atom, Parameter *param)
     MD_FLOAT zlo = 0.0; MD_FLOAT zhi = param->zprd;
     atom->Natoms = 4 * param->nx * param->ny * param->nz;
     atom->Nlocal = 0;
-
-    #ifdef EXPLICIT_TYPES
     atom->ntypes = param->ntypes;
     atom->epsilon = allocate(ALIGNMENT, atom->ntypes * atom->ntypes * sizeof(MD_FLOAT));
     atom->sigma6 = allocate(ALIGNMENT, atom->ntypes * atom->ntypes * sizeof(MD_FLOAT));
@@ -71,7 +67,6 @@ void createAtom(Atom *atom, Parameter *param)
         atom->cutneighsq[i] = param->cutneigh * param->cutneigh;
         atom->cutforcesq[i] = param->cutforce * param->cutforce;
     }
-    #endif
 
     MD_FLOAT alat = pow((4.0 / param->rho), (1.0 / 3.0));
     int ilo = (int) (xlo / (0.5 * alat) - 1);
@@ -142,9 +137,7 @@ void createAtom(Atom *atom, Parameter *param)
                 atom->vx[atom->Nlocal] = vxtmp;
                 atom->vy[atom->Nlocal] = vytmp;
                 atom->vz[atom->Nlocal] = vztmp;
-                #ifdef EXPLICIT_TYPES
                 atom->type[atom->Nlocal] = rand() % atom->ntypes;
-                #endif
                 atom->Nlocal++;
             }
         }
@@ -177,9 +170,7 @@ void growAtom(Atom *atom)
     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));
-    #ifdef EXPLICIT_TYPES
     atom->type = (int *) reallocate(atom->type, ALIGNMENT, atom->Nmax * sizeof(int), nold * sizeof(int));
-    #endif
 }
 
 
diff --git a/src/eam_utils.c b/src/eam_utils.c
new file mode 100644
index 0000000..8ffbad8
--- /dev/null
+++ b/src/eam_utils.c
@@ -0,0 +1,284 @@
+/*
+ * =======================================================================================
+ *
+ *   Author:   Jan Eitzinger (je), jan.eitzinger@fau.de
+ *   Copyright (c) 2020 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 <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include <allocate.h>
+#include <atom.h>
+#include <eam.h>
+#include <parameter.h>
+#include <util.h>
+
+#ifndef MAXLINE
+#define MAXLINE 4096
+#endif
+
+void initEam(Eam* eam, const char* input_file, int ntypes) {
+    eam->nmax = 0;
+    eam->fp = NULL;
+    eam->ntypes = ntypes;
+    eam->cutforcesq = (MD_FLOAT *) allocate(ALIGNMENT, ntypes * ntypes * sizeof(MD_FLOAT));
+    coeff(eam, input_file);
+    init_style(eam);
+}
+
+void coeff(Eam* eam, const char* arg) {
+    read_file(&eam->file, arg);
+    int n = strlen(arg) + 1;
+    int ntypes = eam->ntypes;
+    double cutmax = eam->file.cut;
+    for(int i=0; i<ntypes*ntypes; i++)
+        eam->cutforcesq[i] = cutmax * cutmax;
+}
+
+void init_style(Eam* eam) {
+    // convert read-in file(s) to arrays and spline them
+    file2array(eam);
+    array2spline(eam);
+}
+
+void read_file(Funcfl* file, const char* filename) {
+    FILE* fptr;
+    char line[MAXLINE];
+
+    fptr = fopen(filename, "r");
+    if(fptr == NULL) {
+        printf("Can't open EAM Potential file: %s\n", filename);
+        exit(0);
+    }
+
+    int tmp;
+    fgets(line, MAXLINE, fptr);
+    fgets(line, MAXLINE, fptr);
+    sscanf(line, "%d %lg", &tmp, &(file->mass));
+    fgets(line, MAXLINE, fptr);
+    sscanf(line, "%d %lg %d %lg %lg", &file->nrho, &file->drho, &file->nr, &file->dr, &file->cut);
+
+    //printf("Read: %lf %i %lf %i %lf %lf\n",file->mass,file->nrho,file->drho,file->nr,file->dr,file->cut);
+    file->frho = (MD_FLOAT *) allocate(ALIGNMENT, (file->nrho + 1) * sizeof(MD_FLOAT));
+    file->rhor = (MD_FLOAT *) allocate(ALIGNMENT, (file->nr + 1) * sizeof(MD_FLOAT));
+    file->zr = (MD_FLOAT *) allocate(ALIGNMENT, (file->nr + 1) * sizeof(MD_FLOAT));
+    grab(fptr, file->nrho, file->frho);
+    grab(fptr, file->nr, file->zr);
+    grab(fptr, file->nr, file->rhor);
+    for(int i = file->nrho; i > 0; i--) file->frho[i] = file->frho[i - 1];
+    for(int i = file->nr; i > 0; i--) file->rhor[i] = file->rhor[i - 1];
+    for(int i = file->nr; i > 0; i--) file->zr[i] = file->zr[i - 1];
+    fclose(fptr);
+}
+
+void file2array(Eam* eam) {
+    int i, j, k, m, n;
+    double sixth = 1.0 / 6.0;
+
+    // determine max function params from all active funcfl files
+    // active means some element is pointing at it via map
+    int active;
+    double rmax, rhomax;
+    eam->dr = eam->drho = rmax = rhomax = 0.0;
+    active = 0;
+    Funcfl* file = &eam->file;
+    eam->dr = MAX(eam->dr, file->dr);
+    eam->drho = MAX(eam->drho, file->drho);
+    rmax = MAX(rmax, (file->nr - 1) * file->dr);
+    rhomax = MAX(rhomax, (file->nrho - 1) * file->drho);
+
+    // set nr,nrho from cutoff and spacings
+    // 0.5 is for round-off in divide
+    eam->nr = (int)(rmax / eam->dr + 0.5);
+    eam->nrho = (int)(rhomax / eam->drho + 0.5);
+
+    // ------------------------------------------------------------------
+    // setup frho arrays
+    // ------------------------------------------------------------------
+
+    // allocate frho arrays
+    // nfrho = # of funcfl files + 1 for zero array
+    eam->frho = (MD_FLOAT *) allocate(ALIGNMENT, (eam->nrho + 1) * sizeof(MD_FLOAT));
+
+    // interpolate each file's frho to a single grid and cutoff
+    double r, p, cof1, cof2, cof3, cof4;
+    n = 0;
+
+    for(m = 1; m <= eam->nrho; m++) {
+        r = (m - 1) * eam->drho;
+        p = r / file->drho + 1.0;
+        k = (int)(p);
+        k = MIN(k, file->nrho - 2);
+        k = MAX(k, 2);
+        p -= k;
+        p = MIN(p, 2.0);
+        cof1 = -sixth * p * (p - 1.0) * (p - 2.0);
+        cof2 = 0.5 * (p * p - 1.0) * (p - 2.0);
+        cof3 = -0.5 * p * (p + 1.0) * (p - 2.0);
+        cof4 = sixth * p * (p * p - 1.0);
+        eam->frho[m] = cof1 * file->frho[k - 1] + cof2 * file->frho[k] +
+        cof3 * file->frho[k + 1] + cof4 * file->frho[k + 2];
+    }
+
+
+    // ------------------------------------------------------------------
+    // setup rhor arrays
+    // ------------------------------------------------------------------
+
+    // allocate rhor arrays
+    // nrhor = # of funcfl files
+    eam->rhor = (MD_FLOAT *) allocate(ALIGNMENT, (eam->nr + 1) * sizeof(MD_FLOAT));
+
+    // interpolate each file's rhor to a single grid and cutoff
+    for(m = 1; m <= eam->nr; m++) {
+        r = (m - 1) * eam->dr;
+        p = r / file->dr + 1.0;
+        k = (int)(p);
+        k = MIN(k, file->nr - 2);
+        k = MAX(k, 2);
+        p -= k;
+        p = MIN(p, 2.0);
+        cof1 = -sixth * p * (p - 1.0) * (p - 2.0);
+        cof2 = 0.5 * (p * p - 1.0) * (p - 2.0);
+        cof3 = -0.5 * p * (p + 1.0) * (p - 2.0);
+        cof4 = sixth * p * (p * p - 1.0);
+        eam->rhor[m] = cof1 * file->rhor[k - 1] + cof2 * file->rhor[k] +
+        cof3 * file->rhor[k + 1] + cof4 * file->rhor[k + 2];
+        //if(m==119)printf("BuildRho: %e %e %e %e %e %e\n",rhor[m],cof1,cof2,cof3,cof4,file->rhor[k]);
+    }
+
+    // type2rhor[i][j] = which rhor array (0 to nrhor-1) each type pair maps to
+    // for funcfl files, I,J mapping only depends on I
+    // OK if map = -1 (non-EAM atom in pair hybrid) b/c type2rhor not used
+
+    // ------------------------------------------------------------------
+    // setup z2r arrays
+    // ------------------------------------------------------------------
+
+    // allocate z2r arrays
+    // nz2r = N*(N+1)/2 where N = # of funcfl files
+    eam->z2r = (MD_FLOAT *) allocate(ALIGNMENT, (eam->nr + 1) * sizeof(MD_FLOAT));
+
+    // create a z2r array for each file against other files, only for I >= J
+    // interpolate zri and zrj to a single grid and cutoff
+    double zri, zrj;
+    Funcfl* ifile = &eam->file;
+    Funcfl* jfile = &eam->file;
+
+    for(m = 1; m <= eam->nr; m++) {
+        r = (m - 1) * eam->dr;
+        p = r / ifile->dr + 1.0;
+        k = (int)(p);
+        k = MIN(k, ifile->nr - 2);
+        k = MAX(k, 2);
+        p -= k;
+        p = MIN(p, 2.0);
+        cof1 = -sixth * p * (p - 1.0) * (p - 2.0);
+        cof2 = 0.5 * (p * p - 1.0) * (p - 2.0);
+        cof3 = -0.5 * p * (p + 1.0) * (p - 2.0);
+        cof4 = sixth * p * (p * p - 1.0);
+        zri = cof1 * ifile->zr[k - 1] + cof2 * ifile->zr[k] +
+        cof3 * ifile->zr[k + 1] + cof4 * ifile->zr[k + 2];
+
+        p = r / jfile->dr + 1.0;
+        k = (int)(p);
+        k = MIN(k, jfile->nr - 2);
+        k = MAX(k, 2);
+        p -= k;
+        p = MIN(p, 2.0);
+        cof1 = -sixth * p * (p - 1.0) * (p - 2.0);
+        cof2 = 0.5 * (p * p - 1.0) * (p - 2.0);
+        cof3 = -0.5 * p * (p + 1.0) * (p - 2.0);
+        cof4 = sixth * p * (p * p - 1.0);
+        zrj = cof1 * jfile->zr[k - 1] + cof2 * jfile->zr[k] +
+        cof3 * jfile->zr[k + 1] + cof4 * jfile->zr[k + 2];
+
+        eam->z2r[m] = 27.2 * 0.529 * zri * zrj;
+    }
+}
+
+void array2spline(Eam* eam) {
+    eam->rdr = 1.0 / eam->dr;
+    eam->rdrho = 1.0 / eam->drho;
+    eam->nrho_tot = (eam->nrho + 1) * 7 + 64;
+    eam->nr_tot = (eam->nr + 1) * 7 + 64;
+    eam->nrho_tot -= eam->nrho_tot%64;
+    eam->nr_tot -= eam->nr_tot%64;
+
+    int ntypes = eam->ntypes;
+    eam->frho_spline = (MD_FLOAT *) allocate(ALIGNMENT, ntypes * ntypes * eam->nrho_tot * sizeof(MD_FLOAT));
+    eam->rhor_spline = (MD_FLOAT *) allocate(ALIGNMENT, ntypes * ntypes * eam->nr_tot * sizeof(MD_FLOAT));
+    eam->z2r_spline = (MD_FLOAT *) allocate(ALIGNMENT, ntypes * ntypes * eam->nr_tot * sizeof(MD_FLOAT));
+    interpolate(eam->nrho, eam->drho, eam->frho, eam->frho_spline);
+    interpolate(eam->nr, eam->dr, eam->rhor, eam->rhor_spline);
+    interpolate(eam->nr, eam->dr, eam->z2r, eam->z2r_spline);
+
+    // replicate data for multiple types;
+    for(int tt = 0; tt < ntypes * ntypes; tt++) {
+        for(int k = 0; k < eam->nrho_tot; k++)
+            eam->frho_spline[tt*eam->nrho_tot + k] = eam->frho_spline[k];
+        for(int k = 0; k < eam->nr_tot; k++)
+            eam->rhor_spline[tt*eam->nr_tot + k] = eam->rhor_spline[k];
+        for(int k = 0; k < eam->nr_tot; k++)
+            eam->z2r_spline[tt*eam->nr_tot + k] = eam->z2r_spline[k];
+    }
+}
+
+void interpolate(int n, MD_FLOAT delta, MD_FLOAT* f, MD_FLOAT* spline) {
+    for(int m = 1; m <= n; m++) spline[m * 7 + 6] = f[m];
+
+    spline[1 * 7 + 5] = spline[2 * 7 + 6] - spline[1 * 7 + 6];
+    spline[2 * 7 + 5] = 0.5 * (spline[3 * 7 + 6] - spline[1 * 7 + 6]);
+    spline[(n - 1) * 7 + 5] = 0.5 * (spline[n * 7 + 6] - spline[(n - 2) * 7 + 6]);
+    spline[n * 7 + 5] = spline[n * 7 + 6] - spline[(n - 1) * 7 + 6];
+
+    for(int m = 3; m <= n - 2; m++)
+        spline[m * 7 + 5] = ((spline[(m - 2) * 7 + 6] - spline[(m + 2) * 7 + 6]) +
+                            8.0 * (spline[(m + 1) * 7 + 6] - spline[(m - 1) * 7 + 6])) / 12.0;
+
+    for(int m = 1; m <= n - 1; m++) {
+        spline[m * 7 + 4] = 3.0 * (spline[(m + 1) * 7 + 6] - spline[m * 7 + 6]) -
+                            2.0 * spline[m * 7 + 5] - spline[(m + 1) * 7 + 5];
+        spline[m * 7 + 3] = spline[m * 7 + 5] + spline[(m + 1) * 7 + 5] -
+                            2.0 * (spline[(m + 1) * 7 + 6] - spline[m * 7 + 6]);
+    }
+
+    spline[n * 7 + 4] = 0.0;
+    spline[n * 7 + 3] = 0.0;
+
+    for(int m = 1; m <= n; m++) {
+        spline[m * 7 + 2] = spline[m * 7 + 5] / delta;
+        spline[m * 7 + 1] = 2.0 * spline[m * 7 + 4] / delta;
+        spline[m * 7 + 0] = 3.0 * spline[m * 7 + 3] / delta;
+    }
+}
+
+void grab(FILE* fptr, int n, MD_FLOAT* list) {
+    char* ptr;
+    char line[MAXLINE];
+    int i = 0;
+
+    while(i < n) {
+        fgets(line, MAXLINE, fptr);
+        ptr = strtok(line, " \t\n\r\f");
+        list[i++] = atof(ptr);
+        while(ptr = strtok(NULL, " \t\n\r\f")) list[i++] = atof(ptr);
+    }
+}
diff --git a/src/force_eam.c b/src/force_eam.c
new file mode 100644
index 0000000..ce32d91
--- /dev/null
+++ b/src/force_eam.c
@@ -0,0 +1,177 @@
+/*
+ * =======================================================================================
+ *
+ *   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 <math.h>
+
+#include <allocate.h>
+#include <timing.h>
+#include <neighbor.h>
+#include <parameter.h>
+#include <atom.h>
+#include <stats.h>
+#include <eam.h>
+#include <util.h>
+
+double computeForceEam(Eam* eam, Atom *atom, Neighbor *neighbor, Stats *stats, int first_exec, int timestep) {
+    if(eam->nmax < atom->Nmax) {
+        eam->nmax = atom->Nmax;
+        if(eam->fp != NULL) { free(eam->fp); }
+        eam->fp = (MD_FLOAT *) allocate(ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT));
+    }
+
+    int Nlocal = atom->Nlocal;
+    int* neighs;
+    MD_FLOAT* fx = atom->fx; MD_FLOAT* fy = atom->fy; MD_FLOAT* fz = atom->fz; MD_FLOAT* fp = eam->fp;
+    MD_FLOAT* rhor_spline = eam->rhor_spline; MD_FLOAT* frho_spline = eam->frho_spline; MD_FLOAT* z2r_spline = eam->z2r_spline;
+    int rdr = eam->rdr; int nr = eam->nr; int nr_tot = eam->nr_tot; int 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];
+        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 rhoi = 0;
+        const int type_i = atom->type[i];
+
+        #pragma ivdep
+        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;
+            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];
+
+            if(rsq < cutforcesq) {
+                MD_FLOAT p = sqrt(rsq) * rdr + 1.0;
+                int m = (int)(p);
+                m = m < nr - 1 ? m : nr - 1;
+                p -= m;
+                p = p < 1.0 ? p : 1.0;
+
+                rhoi += ((rhor_spline[type_ij * nr_tot + m * 7 + 3] * p +
+                          rhor_spline[type_ij * nr_tot + m * 7 + 4]) * p +
+                          rhor_spline[type_ij * nr_tot + m * 7 + 5]) * p +
+                          rhor_spline[type_ij * nr_tot + m * 7 + 6];
+            }
+        }
+
+        const int type_ii = type_i * type_i;
+        MD_FLOAT p = 1.0 * rhoi * rdrho + 1.0;
+        int m = (int)(p);
+        m = MAX(1, MIN(m, nrho - 1));
+        p -= m;
+        p = MIN(p, 1.0);
+        fp[i] = (frho_spline[type_ii * nrho_tot + m * 7 + 0] * p +
+                 frho_spline[type_ii * nrho_tot + m * 7 + 1]) * p +
+                 frho_spline[type_ii * nrho_tot + m * 7 + 2];
+    }
+
+    LIKWID_MARKER_STOP("force_eam_fp");
+    LIKWID_MARKER_START("force_eam");
+
+    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;
+        const int type_i = atom->type[i];
+
+        #pragma ivdep
+        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;
+            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];
+
+            if(rsq < cutforcesq) {
+                MD_FLOAT r = sqrt(rsq);
+                MD_FLOAT p = r * rdr + 1.0;
+                int m = (int)(p);
+                m = m < nr - 1 ? m : nr - 1;
+                p -= m;
+                p = p < 1.0 ? p : 1.0;
+
+
+                // rhoip = derivative of (density at atom j due to atom i)
+                // rhojp = derivative of (density at atom i due to atom j)
+                // phi = pair potential energy
+                // phip = phi'
+                // z2 = phi * r
+                // z2p = (phi * r)' = (phi' r) + phi
+                // psip needs both fp[i] and fp[j] terms since r_ij appears in two
+                //   terms of embed eng: Fi(sum rho_ij) and Fj(sum rho_ji)
+                //   hence embed' = Fi(sum rho_ij) rhojp + Fj(sum rho_ji) rhoip
+
+                MD_FLOAT rhoip = (rhor_spline[type_ij * nr_tot + m * 7 + 0] * p +
+                                  rhor_spline[type_ij * nr_tot + m * 7 + 1]) * p +
+                                  rhor_spline[type_ij * nr_tot + m * 7 + 2];
+
+                MD_FLOAT z2p = (z2r_spline[type_ij * nr_tot + m * 7 + 0] * p +
+                                z2r_spline[type_ij * nr_tot + m * 7 + 1]) * p +
+                                z2r_spline[type_ij * nr_tot + m * 7 + 2];
+
+                MD_FLOAT z2 = ((z2r_spline[type_ij * nr_tot + m * 7 + 3] * p +
+                                z2r_spline[type_ij * nr_tot + m * 7 + 4]) * p +
+                                z2r_spline[type_ij * nr_tot + m * 7 + 5]) * p +
+                                z2r_spline[type_ij * nr_tot + m * 7 + 6];
+
+                MD_FLOAT recip = 1.0 / r;
+                MD_FLOAT phi = z2 * recip;
+                MD_FLOAT phip = z2p * recip - phi * recip;
+                MD_FLOAT psip = fp[i] * rhoip + fp[j] * rhoip + phip;
+                MD_FLOAT fpair = -psip * recip;
+
+                fix += delx * fpair;
+                fiy += dely * fpair;
+                fiz += delz * fpair;
+                //fpair *= 0.5;
+            }
+        }
+
+        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);
+    }
+
+    LIKWID_MARKER_STOP("force_eam");
+    double E = getTimeStamp();
+    return E-S;
+}
diff --git a/src/includes/atom.h b/src/includes/atom.h
index ed1a0f9..98d2117 100644
--- a/src/includes/atom.h
+++ b/src/includes/atom.h
@@ -30,14 +30,12 @@ typedef struct {
     MD_FLOAT *x, *y, *z;
     MD_FLOAT *vx, *vy, *vz;
     MD_FLOAT *fx, *fy, *fz;
-    #ifdef EXPLICIT_TYPES
     int *type;
     int ntypes;
     MD_FLOAT *epsilon;
     MD_FLOAT *sigma6;
     MD_FLOAT *cutforcesq;
     MD_FLOAT *cutneighsq;
-    #endif
 } Atom;
 
 extern void initAtom(Atom*);
diff --git a/src/includes/eam.h b/src/includes/eam.h
new file mode 100644
index 0000000..280f419
--- /dev/null
+++ b/src/includes/eam.h
@@ -0,0 +1,57 @@
+/*
+ * =======================================================================================
+ *
+ *   Author:   Jan Eitzinger (je), jan.eitzinger@fau.de
+ *   Copyright (c) 2020 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 <stdio.h>
+
+#include <atom.h>
+#include <parameter.h>
+
+#ifndef __EAM_H_
+#define __EAM_H_
+typedef struct {
+    int nrho, nr;
+    MD_FLOAT drho, dr, cut, mass;
+    MD_FLOAT *frho, *rhor, *zr;
+} Funcfl;
+
+typedef struct {
+    MD_FLOAT* fp;
+    int nmax;
+    int ntypes;
+    int nrho, nr;
+    int nrho_tot, nr_tot;
+    MD_FLOAT dr, rdr, drho, rdrho;
+    MD_FLOAT *frho, *rhor, *z2r;
+    MD_FLOAT *rhor_spline, *frho_spline, *z2r_spline;
+    MD_FLOAT *cutforcesq;
+    Funcfl file;
+} Eam;
+
+void initEam(Eam* eam, const char* input_file, int ntypes);
+void coeff(Eam* eam, const char* arg);
+void init_style(Eam* eam);
+void read_file(Funcfl* file, const char* filename);
+void file2array(Eam* eam);
+void array2spline(Eam* eam);
+void interpolate(int n, MD_FLOAT delta, MD_FLOAT* f, MD_FLOAT* spline);
+void grab(FILE* fptr, int n, MD_FLOAT* list);
+#endif
diff --git a/src/includes/parameter.h b/src/includes/parameter.h
index e95e1e8..ac6c737 100644
--- a/src/includes/parameter.h
+++ b/src/includes/parameter.h
@@ -23,6 +23,9 @@
 #ifndef __PARAMETER_H_
 #define __PARAMETER_H_
 
+#define FF_LJ   0
+#define FF_EAM  1
+
 #if PRECISION == 1
 #define MD_FLOAT float
 #else
@@ -30,6 +33,8 @@
 #endif
 
 typedef struct {
+    int force_field;
+    char* input_file;
     MD_FLOAT epsilon;
     MD_FLOAT sigma6;
     MD_FLOAT temp;
diff --git a/src/main-stub.c b/src/main-stub.c
index 72521f2..07fbbe8 100644
--- a/src/main-stub.c
+++ b/src/main-stub.c
@@ -128,7 +128,6 @@ int main(int argc, const char *argv[]) {
     initAtom(atom);
     initStats(&stats);
 
-    #ifdef EXPLICIT_TYPES
     atom->ntypes = param.ntypes;
     atom->epsilon = allocate(ALIGNMENT, atom->ntypes * atom->ntypes * sizeof(MD_FLOAT));
     atom->sigma6 = allocate(ALIGNMENT, atom->ntypes * atom->ntypes * sizeof(MD_FLOAT));
@@ -140,7 +139,6 @@ int main(int argc, const char *argv[]) {
         atom->cutneighsq[i] = param.cutneigh * param.cutneigh;
         atom->cutforcesq[i] = param.cutforce * param.cutforce;
     }
-    #endif
 
     DEBUG("Creating atoms...\n");
     for(int i = 0; i < param.nx; ++i) {
@@ -176,9 +174,7 @@ int main(int argc, const char *argv[]) {
                     for(int jj = 0; jj < fac_y; ++jj) {
                         for(int kk = 0; kk < fac_z; ++kk) {
                             if(added_atoms < atoms_per_unit_cell) {
-                                #ifdef EXPLICIT_TYPES
                                 atom->type[atom->Nlocal] = rand() % atom->ntypes;
-                                #endif
                                 ADD_ATOM(ii * offset_x, jj * offset_y, kk * offset_z, vx, vy, vz);
                                 added_atoms++;
                             }
diff --git a/src/main.c b/src/main.c
index 703ca9f..43391d5 100644
--- a/src/main.c
+++ b/src/main.c
@@ -39,14 +39,18 @@
 #include <thermo.h>
 #include <pbc.h>
 #include <timers.h>
+#include <eam.h>
 
 #define HLINE "----------------------------------------------------------------------------\n"
 
 extern double computeForce(Parameter*, Atom*, Neighbor*);
 extern double computeForceTracing(Parameter*, Atom*, Neighbor*, Stats*, int, int);
+extern double computeForceEam(Eam* eam, Atom *atom, Neighbor *neighbor, Stats *stats, int first_exec, int timestep);
 
 void init(Parameter *param)
 {
+    param->input_file = NULL;
+    param->force_field = FF_LJ;
     param->epsilon = 1.0;
     param->sigma6 = 1.0;
     param->rho = 0.8442;
@@ -68,6 +72,7 @@ void init(Parameter *param)
 
 double setup(
         Parameter *param,
+        Eam *eam,
         Atom *atom,
         Neighbor *neighbor,
         Stats *stats)
@@ -79,6 +84,7 @@ double setup(
     param->zprd = param->nz * param->lattice;
 
     S = getTimeStamp();
+    if(param->force_field == FF_EAM) { initEam(eam, param->input_file, param->ntypes); }
     initAtom(atom);
     initNeighbor(neighbor, param);
     initPbc();
@@ -147,16 +153,31 @@ void printAtomState(Atom *atom)
     printf("Atom counts: Natoms=%d Nlocal=%d Nghost=%d Nmax=%d\n",
             atom->Natoms, atom->Nlocal, atom->Nghost, atom->Nmax);
 
-/*     int nall = atom->Nlocal + atom->Nghost; */
+    /*     int nall = atom->Nlocal + atom->Nghost; */
 
-/*     for (int i=0; i<nall; i++) { */
-/*         printf("%d  %f %f %f\n", i, atom->x[i], atom->y[i], atom->z[i]); */
-/*     } */
+    /*     for (int i=0; i<nall; i++) { */
+    /*         printf("%d  %f %f %f\n", i, atom->x[i], atom->y[i], atom->z[i]); */
+    /*     } */
 }
 
-int main (int argc, char** argv)
+int str2ff(const char *string)
+{
+    if(strncmp(string, "lj", 2) == 0) return FF_LJ;
+    if(strncmp(string, "eam", 3) == 0) return FF_EAM;
+    return -1;
+}
+
+const char* ff2str(int ff)
+{
+    if(ff == FF_LJ) { return "lj"; }
+    if(ff == FF_EAM) { return "eam"; }
+    return "invalid";
+}
+
+int main(int argc, char** argv)
 {
     double timer[NUMTIMER];
+    Eam eam;
     Atom atom;
     Neighbor neighbor;
     Stats stats;
@@ -173,6 +194,19 @@ int main (int argc, char** argv)
 
     for(int i = 0; i < argc; i++)
     {
+        if((strcmp(argv[i], "-f") == 0))
+        {
+            if((param.force_field = str2ff(argv[++i])) < 0) {
+                fprintf(stderr, "Invalid force field!\n");
+                exit(-1);
+            }
+            continue;
+        }
+        if((strcmp(argv[i], "-i") == 0))
+        {
+            param.input_file = strdup(argv[++i]);
+            continue;
+        }
         if((strcmp(argv[i], "-n") == 0) || (strcmp(argv[i], "--nsteps") == 0))
         {
             param.ntimes = atoi(argv[++i]);
@@ -193,7 +227,7 @@ int main (int argc, char** argv)
             param.nz = atoi(argv[++i]);
             continue;
         }
-        if((strcmp(argv[i], "-f") == 0))
+        if((strcmp(argv[i], "--freq") == 0))
         {
             param.proc_freq = atof(argv[++i]);
             continue;
@@ -202,21 +236,27 @@ int main (int argc, char** argv)
         {
             printf("MD Bench: A minimalistic re-implementation of miniMD\n");
             printf(HLINE);
+            printf("-f <string>:          force field (lj or eam), default lj\n");
+            printf("-i <string>:          input file for EAM\n");
             printf("-n / --nsteps <int>:  set number of timesteps for simulation\n");
             printf("-nx/-ny/-nz <int>:    set linear dimension of systembox in x/y/z direction\n");
-            printf("-f <real>:            processor frequency (GHz)\n");
+            printf("--freq <real>:        processor frequency (GHz)\n");
             printf(HLINE);
             exit(EXIT_SUCCESS);
         }
     }
 
-    setup(&param, &atom, &neighbor, &stats);
+    setup(&param, &eam, &atom, &neighbor, &stats);
     computeThermo(0, &param, &atom);
-#if defined(MEM_TRACER) || defined(INDEX_TRACER) || defined(PRINT_STATS)
-    computeForceTracing(&param, &atom, &neighbor, &stats, 1, 0);
+    if(param.force_field == FF_EAM) {
+        computeForceEam(&eam, &atom, &neighbor, &stats, 1, 0);
+    } else {
+#if defined(MEM_TRACER) || defined(INDEX_TRACER) || defined(COMPUTE_STATS)
+        computeForceTracing(&param, &atom, &neighbor, &stats, 1, 0);
 #else
-    computeForce(&param, &atom, &neighbor);
+        computeForce(&param, &atom, &neighbor);
 #endif
+    }
 
     timer[FORCE] = 0.0;
     timer[NEIGH] = 0.0;
@@ -231,11 +271,15 @@ int main (int argc, char** argv)
             timer[NEIGH] += reneighbour(&param, &atom, &neighbor);
         }
 
-#if defined(MEM_TRACER) || defined(INDEX_TRACER) || defined(PRINT_STATS)
-        timer[FORCE] += computeForceTracing(&param, &atom, &neighbor, &stats, 0, n + 1);
+        if(param.force_field == FF_EAM) {
+            timer[FORCE] += computeForceEam(&eam, &atom, &neighbor, &stats, 0, n + 1);
+        } else {
+#if defined(MEM_TRACER) || defined(INDEX_TRACER) || defined(COMPUTE_STATS)
+            timer[FORCE] += computeForceTracing(&param, &atom, &neighbor, &stats, 0, n + 1);
 #else
-        timer[FORCE] += computeForce(&param, &atom, &neighbor);
+            timer[FORCE] += computeForce(&param, &atom, &neighbor);
 #endif
+        }
         finalIntegrate(&param, &atom);
 
         if(!((n + 1) % param.nstat) && (n+1) < param.ntimes) {
@@ -247,6 +291,7 @@ int main (int argc, char** argv)
     computeThermo(-1, &param, &atom);
 
     printf(HLINE);
+    printf("Force field: %s\n", ff2str(param.force_field));
     printf("Data layout for positions: %s\n", POS_DATA_LAYOUT);
 #if PRECISION == 1
     printf("Using single precision floating point.\n");
@@ -260,7 +305,7 @@ int main (int argc, char** argv)
     printf(HLINE);
     printf("Performance: %.2f million atom updates per second\n",
             1e-6 * (double) atom.Natoms * param.ntimes / timer[TOTAL]);
-#ifdef PRINT_STATS
+#ifdef COMPUTE_STATS
     displayStatistics(&atom, &param, &stats, timer);
 #endif
     LIKWID_MARKER_CLOSE;
diff --git a/src/stats.c b/src/stats.c
index b31b35b..71f7a00 100644
--- a/src/stats.c
+++ b/src/stats.c
@@ -14,11 +14,15 @@ void displayStatistics(Atom *atom, Parameter *param, Stats *stats, double *timer
 #ifdef COMPUTE_STATS
     double force_useful_volume = 1e-9 * ( (double)(atom->Nlocal * (param->ntimes + 1)) * (sizeof(MD_FLOAT) * 6 + sizeof(int)) +
                                           (double)(stats->total_force_neighs) * (sizeof(MD_FLOAT) * 3 + sizeof(int)) );
+    double avg_neigh = stats->total_force_neighs / (double)(atom->Nlocal * (param->ntimes + 1));
+    double avg_simd = stats->total_force_iters / (double)(atom->Nlocal * (param->ntimes + 1));
 #ifdef EXPLICIT_TYPES
     force_useful_volume += 1e-9 * (double)((atom.Nlocal * (param.ntimes + 1)) + stats.total_force_neighs) * sizeof(int);
 #endif
     printf("Statistics:\n");
     printf("\tVector width: %d, Processor frequency: %.4f GHz\n", VECTOR_WIDTH, param->proc_freq);
+    printf("\tAverage neighbors per atom: %.4f\n", avg_neigh);
+    printf("\tAverage SIMD iterations per atom: %.4f\n", avg_simd);
     printf("\tTotal number of computed pair interactions: %lld\n", stats->total_force_neighs);
     printf("\tTotal number of SIMD iterations: %lld\n", stats->total_force_iters);
     printf("\tUseful read data volume for force computation: %.2fGB\n", force_useful_volume);