diff --git a/Makefile b/Makefile
index 9d3cc90..1738a3d 100644
--- a/Makefile
+++ b/Makefile
@@ -51,6 +51,10 @@ ifeq ($(strip $(USE_REFERENCE_VERSION)),true)
     DEFINES += -DUSE_REFERENCE_VERSION
 endif
 
+ifeq ($(strip $(HALF_NEIGHBOR_LISTS_CHECK_CJ)),true)
+    DEFINES += -DHALF_NEIGHBOR_LISTS_CHECK_CJ
+endif
+
 ifeq ($(strip $(DEBUG)),true)
     DEFINES += -DDEBUG
 endif
diff --git a/config.mk b/config.mk
index 9afe7a0..1e8318b 100644
--- a/config.mk
+++ b/config.mk
@@ -29,6 +29,8 @@ COMPUTE_STATS ?= true
 USE_REFERENCE_VERSION ?= false
 # Enable XTC output
 XTC_OUTPUT ?= false
+# Check if cj is local when decreasing reaction force
+HALF_NEIGHBOR_LISTS_CHECK_CJ ?= false
 
 #Feature options
 OPTIONS =  -DALIGNMENT=64
diff --git a/data/argon_1000/mdbench_params.conf b/data/argon_1000/mdbench_params.conf
index 474b886..836e678 100644
--- a/data/argon_1000/mdbench_params.conf
+++ b/data/argon_1000/mdbench_params.conf
@@ -9,4 +9,3 @@ v_out_freq 5
 cutforce 0.9
 reneigh_every 100
 nstat 125000
-half_neigh 0
diff --git a/gromacs/force_lj.c b/gromacs/force_lj.c
index 58aad50..6754506 100644
--- a/gromacs/force_lj.c
+++ b/gromacs/force_lj.c
@@ -263,7 +263,13 @@ double computeForceLJ_2xnn_half(Parameter *param, Atom *atom, Neighbor *neighbor
             fiy2 = simd_add(fiy2, ty2);
             fiz2 = simd_add(fiz2, tz2);
 
+            #ifdef HALF_NEIGHBOR_LISTS_CHECK_CJ
+            if(cj < CJ1_FROM_CI(atom->Nlocal)) {
+                simd_h_decr3(cj_f, tx0 + tx2, ty0 + ty2, tz0 + tz2);
+            }
+            #else
             simd_h_decr3(cj_f, tx0 + tx2, ty0 + ty2, tz0 + tz2);
+            #endif
         }
 
         simd_h_dual_incr_reduced_sum(&ci_f[CL_X_OFFSET], fix0, fix2);
@@ -576,9 +582,17 @@ double computeForceLJ_4xn_half(Parameter *param, Atom *atom, Neighbor *neighbor,
             fiy3 = simd_add(fiy3, ty3);
             fiz3 = simd_add(fiz3, tz3);
 
+            #ifdef HALF_NEIGHBOR_LISTS_CHECK_CJ
+            if(cj < CJ1_FROM_CI(atom->Nlocal)) {
+                simd_store(&cj_f[CL_X_OFFSET], simd_load(&cj_f[CL_X_OFFSET]) - (tx0 + tx1 + tx2 + tx3));
+                simd_store(&cj_f[CL_Y_OFFSET], simd_load(&cj_f[CL_Y_OFFSET]) - (ty0 + ty1 + ty2 + ty3));
+                simd_store(&cj_f[CL_Z_OFFSET], simd_load(&cj_f[CL_Z_OFFSET]) - (tz0 + tz1 + tz2 + tz3));
+            }
+            #else
             simd_store(&cj_f[CL_X_OFFSET], simd_load(&cj_f[CL_X_OFFSET]) - (tx0 + tx1 + tx2 + tx3));
             simd_store(&cj_f[CL_Y_OFFSET], simd_load(&cj_f[CL_Y_OFFSET]) - (ty0 + ty1 + ty2 + ty3));
             simd_store(&cj_f[CL_Z_OFFSET], simd_load(&cj_f[CL_Z_OFFSET]) - (tz0 + tz1 + tz2 + tz3));
+            #endif
         }
 
         simd_incr_reduced_sum(&ci_f[CL_X_OFFSET], fix0, fix1, fix2, fix3);