diff --git a/lammps/force_lj.c b/lammps/force_lj.c
index 73c0269..ecc1719 100644
--- a/lammps/force_lj.c
+++ b/lammps/force_lj.c
@@ -28,7 +28,8 @@
 #include <atom.h>
 #include <stats.h>
 
-double computeForceLJ(Parameter *param, Atom *atom, Neighbor *neighbor, Stats *stats) {
+double computeForceLJFullNeigh(Parameter *param, Atom *atom, Neighbor *neighbor, Stats *stats)
+{
     int Nlocal = atom->Nlocal;
     int* neighs;
     MD_FLOAT* fx = atom->fx;
@@ -64,7 +65,6 @@ double computeForceLJ(Parameter *param, Atom *atom, Neighbor *neighbor, Stats *s
         const int type_i = atom->type[i];
 #endif
 
-
         for(int k = 0; k < numneighs; k++) {
             int j = neighs[k];
             MD_FLOAT delx = xtmp - atom_x(j);
@@ -102,3 +102,81 @@ double computeForceLJ(Parameter *param, Atom *atom, Neighbor *neighbor, Stats *s
     double E = getTimeStamp();
     return E-S;
 }
+
+double computeForceLJHalfNeigh(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;
+    }
+
+    double S = getTimeStamp();
+    LIKWID_MARKER_START("forceLJ-halfneigh");
+
+    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 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;
+                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[j] -= delx * force;
+                fy[j] -= dely * force;
+                fz[j] -= delz * force;
+            }
+        }
+
+        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("forceLJ-halfneigh");
+    double E = getTimeStamp();
+    return E-S;
+}
diff --git a/lammps/includes/parameter.h b/lammps/includes/parameter.h
index bff46dd..31aa062 100644
--- a/lammps/includes/parameter.h
+++ b/lammps/includes/parameter.h
@@ -42,6 +42,7 @@ typedef struct {
     int ntimes;
     int nstat;
     int every;
+    int halfneigh;
     MD_FLOAT dt;
     MD_FLOAT dtforce;
     MD_FLOAT cutforce;
diff --git a/lammps/main.c b/lammps/main.c
index ebb2dbd..ffdb863 100644
--- a/lammps/main.c
+++ b/lammps/main.c
@@ -45,7 +45,8 @@
 
 #define HLINE "----------------------------------------------------------------------------\n"
 
-extern double computeForceLJ(Parameter*, Atom*, Neighbor*, Stats*);
+extern double computeForceLJFullNeigh(Parameter*, Atom*, Neighbor*, Stats*);
+extern double computeForceLJHalfNeigh(Parameter*, Atom*, Neighbor*, Stats*);
 extern double computeForceEam(Eam*, Parameter*, Atom*, Neighbor*, Stats*);
 
 void init(Parameter *param)
@@ -70,6 +71,7 @@ void init(Parameter *param)
     param->dtforce = 0.5 * param->dt;
     param->every = 20;
     param->proc_freq = 2.4;
+    param->halfneigh = 0;
 }
 
 double setup(
@@ -258,7 +260,11 @@ int main(int argc, char** argv)
     if(param.force_field == FF_EAM) {
         timer[FORCE] = computeForceEam(&eam, &param, &atom, &neighbor, &stats);
     } else {
-        timer[FORCE] = computeForceLJ(&param, &atom, &neighbor, &stats);
+        if( param.halfneigh ) {
+            timer[FORCE] = computeForceLJHalfNeigh(&param, &atom, &neighbor, &stats);
+        } else {
+            timer[FORCE] = computeForceLJFullNeigh(&param, &atom, &neighbor, &stats);
+        }
     }
 
     timer[NEIGH] = 0.0;
@@ -284,7 +290,11 @@ int main(int argc, char** argv)
         if(param.force_field == FF_EAM) {
             timer[FORCE] += computeForceEam(&eam, &param, &atom, &neighbor, &stats);
         } else {
-            timer[FORCE] += computeForceLJ(&param, &atom, &neighbor, &stats);
+            if( param.halfneigh ) {
+                timer[FORCE] = computeForceLJHalfNeigh(&param, &atom, &neighbor, &stats);
+            } else {
+                timer[FORCE] = computeForceLJFullNeigh(&param, &atom, &neighbor, &stats);
+            }
         }
 
         finalIntegrate(&param, &atom);
diff --git a/lammps/neighbor.c b/lammps/neighbor.c
index 6b512a0..549ddbc 100644
--- a/lammps/neighbor.c
+++ b/lammps/neighbor.c
@@ -46,6 +46,7 @@ static int nmax;
 static int nstencil;      // # of bins in stencil
 static int* stencil;      // stencil list of bin offsets
 static MD_FLOAT binsizex, binsizey, binsizez;
+static int halfneigh; //use half neighbour list
 
 static int coord2bin(MD_FLOAT, MD_FLOAT , MD_FLOAT);
 static MD_FLOAT bindist(int, int, int);
@@ -69,6 +70,7 @@ void initNeighbor(Neighbor *neighbor, Parameter *param)
     neighbor->maxneighs = 100;
     neighbor->numneigh = NULL;
     neighbor->neighbors = NULL;
+    neighbor->halfneigh = param->halfneigh;
 }
 
 void setupNeighbor(Parameter* param)