diff --git a/src/atom.c b/src/atom.c
index 687de55..4e28d7a 100644
--- a/src/atom.c
+++ b/src/atom.c
@@ -41,6 +41,14 @@ 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)
@@ -50,6 +58,21 @@ 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));
+    atom->cutforcesq = allocate(ALIGNMENT, atom->ntypes * atom->ntypes * sizeof(MD_FLOAT));
+    atom->cutneighsq = allocate(ALIGNMENT, atom->ntypes * atom->ntypes * sizeof(MD_FLOAT));
+    for(int i = 0; i < atom->ntypes * atom->ntypes; i++) {
+        atom->epsilon[i] = param->epsilon;
+        atom->sigma6[i] = param->sigma6;
+        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);
     int ihi = (int) (xhi / (0.5 * alat) + 1);
@@ -119,6 +142,9 @@ 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++;
             }
         }
@@ -151,6 +177,9 @@ 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/force.c b/src/force.c
index bd471c2..d1dda58 100644
--- a/src/force.c
+++ b/src/force.c
@@ -27,14 +27,16 @@
 #include <parameter.h>
 #include <atom.h>
 
-double computeForce(Parameter *param, Atom *atom, Neighbor *neighbor) {
+double computeForce(Parameter *param, Atom *atom, Neighbor *neighbor, int ntimes) {
+    double S = getTimeStamp();
     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;
-    MD_FLOAT* fx = atom->fx; MD_FLOAT* fy = atom->fy; MD_FLOAT* fz = atom->fz;
-    double S, E;
+    #endif
 
     for(int i = 0; i < Nlocal; i++) {
         fx[i] = 0.0;
@@ -42,7 +44,7 @@ double computeForce(Parameter *param, Atom *atom, Neighbor *neighbor) {
         fz[i] = 0.0;
     }
 
-#pragma omp parallel for
+    #pragma omp parallel for
     for(int i = 0; i < Nlocal; i++) {
         neighs = &neighbor->neighbors[i * neighbor->maxneighs];
         int numneighs = neighbor->numneigh[i];
@@ -52,16 +54,24 @@ double computeForce(Parameter *param, Atom *atom, Neighbor *neighbor) {
         MD_FLOAT fix = 0;
         MD_FLOAT fiy = 0;
         MD_FLOAT fiz = 0;
+        #ifdef EXPLICIT_TYPES
+        const int type_i = atom->type[i];
+        #endif
 
-//	printf("%d: %d\n", i, numneighs);
-
-        for(int n = 0; n < INTERNAL_LOOP_NTIMES; n++) {
+        for(int n = 0; n < ntimes; n++) {
             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;
+                #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
 
                 if(rsq < cutforcesq) {
                     MD_FLOAT sr2 = 1.0 / rsq;
@@ -79,5 +89,6 @@ double computeForce(Parameter *param, Atom *atom, Neighbor *neighbor) {
         fz[i] += fiz;
     }
 
-    return 0.0;
+    double E = getTimeStamp();
+    return E-S;
 }
diff --git a/src/includes/atom.h b/src/includes/atom.h
index e090936..ed1a0f9 100644
--- a/src/includes/atom.h
+++ b/src/includes/atom.h
@@ -30,6 +30,14 @@ 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/parameter.h b/src/includes/parameter.h
index 9d51b12..0ddec95 100644
--- a/src/includes/parameter.h
+++ b/src/includes/parameter.h
@@ -40,6 +40,7 @@ typedef struct {
     MD_FLOAT temp;
     MD_FLOAT rho;
     MD_FLOAT mass;
+    int ntypes;
     int ntimes;
     int nstat;
     int every;
diff --git a/src/main-stub.c b/src/main-stub.c
index ad2e580..a0973f1 100644
--- a/src/main-stub.c
+++ b/src/main-stub.c
@@ -16,12 +16,13 @@
 #define LATTICE_DISTANCE    10.0
 #define NEIGH_DISTANCE      1.0
 
-extern double computeForce( Parameter*, Atom*, Neighbor*);
+extern double computeForce( Parameter*, Atom*, Neighbor*, int);
 
 void init(Parameter *param) {
     param->epsilon = 1.0;
     param->sigma6 = 1.0;
     param->rho = 0.8442;
+    param->ntypes = 4;
     param->ntimes = 200;
     param->nx = 4;
     param->ny = 4;
@@ -49,6 +50,7 @@ void init(Parameter *param) {
                                         atom->vx[atom->Nlocal] = vy;                        \
                                         atom->vy[atom->Nlocal] = vy;                        \
                                         atom->vz[atom->Nlocal] = vz;                        \
+                                        atom->type[atom->Nlocal] = rand() % atom->ntypes;   \
                                         atom->Nlocal++
 
 int main(int argc, const char *argv[]) {
@@ -123,6 +125,20 @@ int main(int argc, const char *argv[]) {
     DEBUG("Initializing atoms...\n");
     initAtom(atom);
 
+    #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));
+    atom->cutforcesq = allocate(ALIGNMENT, atom->ntypes * atom->ntypes * sizeof(MD_FLOAT));
+    atom->cutneighsq = allocate(ALIGNMENT, atom->ntypes * atom->ntypes * sizeof(MD_FLOAT));
+    for(int i = 0; i < atom->ntypes * atom->ntypes; i++) {
+        atom->epsilon[i] = param->epsilon;
+        atom->sigma6[i] = param->sigma6;
+        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) {
         for(int j = 0; j < param.ny; ++j) {
@@ -173,6 +189,7 @@ int main(int argc, const char *argv[]) {
 
     if(!csv) {
         printf("Number of timesteps: %d\n", param.ntimes);
+        printf("Number of times to compute the most internal loop: %d\n", INTERNAL_LOOP_NTIMES);
         printf("System size (unit cells): %dx%dx%d\n", param.nx, param.ny, param.nz);
         printf("Atoms per unit cell: %d\n", atoms_per_unit_cell);
         printf("Total number of atoms: %d\n", atom->Nlocal);
@@ -188,13 +205,13 @@ int main(int argc, const char *argv[]) {
     DEBUG("Building neighbor lists...\n");
     buildNeighbor(atom, &neighbor);
     DEBUG("Computing forces...\n");
-    computeForce(&param, atom, &neighbor);
+    computeForce(&param, atom, &neighbor, 1);
 
     double S, E;
     S = getTimeStamp();
     LIKWID_MARKER_START("force");
     for(int i = 0; i < param.ntimes; i++) {
-        computeForce(&param, atom, &neighbor);
+        computeForce(&param, atom, &neighbor, INTERNAL_LOOP_NTIMES);
     }
     LIKWID_MARKER_STOP("force");
     E = getTimeStamp();
diff --git a/src/main.c b/src/main.c
index 6b09c8e..061da97 100644
--- a/src/main.c
+++ b/src/main.c
@@ -47,13 +47,14 @@ typedef enum {
     NUMTIMER
 } timertype;
 
-extern double computeForce( Parameter*, Atom*, Neighbor*);
+extern double computeForce( Parameter*, Atom*, Neighbor*, int);
 
 void init(Parameter *param)
 {
     param->epsilon = 1.0;
     param->sigma6 = 1.0;
     param->rho = 0.8442;
+    param->ntypes = 4;
     param->ntimes = 200;
     param->dt = 0.005;
     param->nx = 32;
@@ -205,7 +206,7 @@ int main (int argc, char** argv)
 
     setup(&param, &atom, &neighbor);
     computeThermo(0, &param, &atom);
-    computeForce(&param, &atom, &neighbor);
+    computeForce(&param, &atom, &neighbor, 1);
 
     timer[FORCE] = 0.0;
     timer[NEIGH] = 0.0;
@@ -221,7 +222,7 @@ int main (int argc, char** argv)
             timer[NEIGH] += reneighbour(&param, &atom, &neighbor);
         }
 
-        timer[FORCE] += computeForce(&param, &atom, &neighbor);
+        timer[FORCE] += computeForce(&param, &atom, &neighbor, 1);
         finalIntegrate(&param, &atom);
 
         if(!((n + 1) % param.nstat) && (n+1) < param.ntimes) {
diff --git a/src/neighbor.c b/src/neighbor.c
index 72a6ee1..91c7ab8 100644
--- a/src/neighbor.c
+++ b/src/neighbor.c
@@ -196,7 +196,9 @@ void buildNeighbor(Atom *atom, Neighbor *neighbor)
             MD_FLOAT ytmp = atom_y(i);
             MD_FLOAT ztmp = atom_z(i);
             int ibin = coord2bin(xtmp, ytmp, ztmp);
-
+            #ifdef EXPLICIT_TYPES
+            int type_i = atom->type[i];
+            #endif
             for(int k = 0; k < nstencil; k++) {
                 int jbin = ibin + stencil[k];
                 int* loc_bin = &bins[jbin * atoms_per_bin];
@@ -213,7 +215,13 @@ void buildNeighbor(Atom *atom, Neighbor *neighbor)
                     MD_FLOAT delz = ztmp - atom_z(j);
                     MD_FLOAT rsq = delx * delx + dely * dely + delz * delz;
 
-                    if( rsq <= cutneighsq ) {
+                    #ifdef EXPLICIT_TYPES
+                    int type_j = atom->type[j];
+                    const MD_FLOAT cutoff = atom->cutneighsq[type_i * atom->ntypes + type_j];
+                    #else
+                    const MD_FLOAT cutoff = cutneighsq;
+                    #endif
+                    if( rsq <= cutoff ) {
                         neighptr[n++] = j;
                     }
                 }