diff --git a/lammps/force_lj.c b/lammps/force_lj.c
index 82fe953..ec1297b 100644
--- a/lammps/force_lj.c
+++ b/lammps/force_lj.c
@@ -136,6 +136,8 @@ double computeForceLJHalfNeigh(Parameter *param, Atom *atom, Neighbor *neighbor,
         const int type_i = atom->type[i];
 #endif
 
+        // Pragma required to vectorize the inner loop
+        #pragma simd reduction(+: fix,fiy,fiz)
         for(int k = 0; k < numneighs; k++) {
             int j = neighs[k];
             MD_FLOAT delx = xtmp - atom_x(j);