diff --git a/src/OPENMP/pair_sw_omp.cpp b/src/OPENMP/pair_sw_omp.cpp
index f6d615b2a1..32df823e31 100644
--- a/src/OPENMP/pair_sw_omp.cpp
+++ b/src/OPENMP/pair_sw_omp.cpp
@@ -156,50 +156,51 @@ void PairSWOMP::eval(int iifrom, int iito, ThrData * const thr)
       if (EVFLAG) ev_tally_thr(this,i,j,nlocal,/* newton_pair */ 1,
                                evdwl,0.0,fpair,delx,dely,delz,thr);
     }
-
-    jnumm1 = numshort - 1;
-
-    for (jj = 0; jj < jnumm1; jj++) {
-      j = neighshort_thr[jj];
-      jtype = map[type[j]];
-      ijparam = elem3param[itype][jtype][jtype];
-      delr1[0] = x[j].x - xtmp;
-      delr1[1] = x[j].y - ytmp;
-      delr1[2] = x[j].z - ztmp;
-      rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2];
-
-      double fjxtmp,fjytmp,fjztmp;
-      fjxtmp = fjytmp = fjztmp = 0.0;
-
-      for (kk = jj+1; kk < numshort; kk++) {
-        k = neighshort_thr[kk];
-        ktype = map[type[k]];
-        ikparam = elem3param[itype][ktype][ktype];
-        ijkparam = elem3param[itype][jtype][ktype];
-
-        delr2[0] = x[k].x - xtmp;
-        delr2[1] = x[k].y - ytmp;
-        delr2[2] = x[k].z - ztmp;
-        rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2];
-
-        threebody(&params[ijparam],&params[ikparam],&params[ijkparam],
-                  rsq1,rsq2,delr1,delr2,fj,fk,EFLAG,evdwl);
-
-        fxtmp -= fj[0] + fk[0];
-        fytmp -= fj[1] + fk[1];
-        fztmp -= fj[2] + fk[2];
-        fjxtmp += fj[0];
-        fjytmp += fj[1];
-        fjztmp += fj[2];
-        f[k].x += fk[0];
-        f[k].y += fk[1];
-        f[k].z += fk[2];
-
-        if (EVFLAG) ev_tally3_thr(this,i,j,k,evdwl,0.0,fj,fk,delr1,delr2,thr);
+    if (threebody_on) {
+      jnumm1 = numshort - 1;
+    
+      for (jj = 0; jj < jnumm1; jj++) {
+        j = neighshort_thr[jj];
+        jtype = map[type[j]];
+        ijparam = elem3param[itype][jtype][jtype];
+        delr1[0] = x[j].x - xtmp;
+        delr1[1] = x[j].y - ytmp;
+        delr1[2] = x[j].z - ztmp;
+        rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2];
+    
+        double fjxtmp,fjytmp,fjztmp;
+        fjxtmp = fjytmp = fjztmp = 0.0;
+    
+        for (kk = jj+1; kk < numshort; kk++) {
+          k = neighshort_thr[kk];
+          ktype = map[type[k]];
+          ikparam = elem3param[itype][ktype][ktype];
+          ijkparam = elem3param[itype][jtype][ktype];
+    
+          delr2[0] = x[k].x - xtmp;
+          delr2[1] = x[k].y - ytmp;
+          delr2[2] = x[k].z - ztmp;
+          rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2];
+    
+          threebody(&params[ijparam],&params[ikparam],&params[ijkparam],
+                    rsq1,rsq2,delr1,delr2,fj,fk,EFLAG,evdwl);
+    
+          fxtmp -= fj[0] + fk[0];
+          fytmp -= fj[1] + fk[1];
+          fztmp -= fj[2] + fk[2];
+          fjxtmp += fj[0];
+          fjytmp += fj[1];
+          fjztmp += fj[2];
+          f[k].x += fk[0];
+          f[k].y += fk[1];
+          f[k].z += fk[2];
+    
+          if (EVFLAG) ev_tally3_thr(this,i,j,k,evdwl,0.0,fj,fk,delr1,delr2,thr);
+        }
+        f[j].x += fjxtmp;
+        f[j].y += fjytmp;
+        f[j].z += fjztmp;
       }
-      f[j].x += fjxtmp;
-      f[j].y += fjytmp;
-      f[j].z += fjztmp;
     }
     f[i].x += fxtmp;
     f[i].y += fytmp;