diff --git a/src/USER-MEAMC/meam.h b/src/USER-MEAMC/meam.h
index 944b3fdf28..c562cb52d9 100644
--- a/src/USER-MEAMC/meam.h
+++ b/src/USER-MEAMC/meam.h
@@ -248,5 +248,13 @@ static inline void setall3d(TYPE (&arr)[maxi][maxj][maxk], const TYPE v) {
         arr[i][j][k] = v;
 }
 
+// Helper functions
+
+static inline double fdiv_zero(const double n, const double d) {
+  if (iszero(d))
+    return 0.0;
+  return n / d;
+}
+
 };
 #endif
diff --git a/src/USER-MEAMC/meam_dens_final.cpp b/src/USER-MEAMC/meam_dens_final.cpp
index e25c2e7299..de188e497d 100644
--- a/src/USER-MEAMC/meam_dens_final.cpp
+++ b/src/USER-MEAMC/meam_dens_final.cpp
@@ -33,9 +33,9 @@ MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_
 
       if (rho0[i] > 0.0) {
         if (this->ialloy == 1) {
-          t_ave[i][0] = (tsq_ave[i][0] != 0.0 ) ? t_ave[i][0] / tsq_ave[i][0] : 0.0;
-          t_ave[i][1] = (tsq_ave[i][1] != 0.0 ) ? t_ave[i][1] / tsq_ave[i][1] : 0.0;
-          t_ave[i][0] = (tsq_ave[i][2] != 0.0 ) ? t_ave[i][2] / tsq_ave[i][2] : 0.0;
+          t_ave[i][0] = fdiv_zero(t_ave[i][0], tsq_ave[i][0]);
+          t_ave[i][1] = fdiv_zero(t_ave[i][1], tsq_ave[i][1]);
+          t_ave[i][2] = fdiv_zero(t_ave[i][2], tsq_ave[i][2]);
         } else if (this->ialloy == 2) {
           t_ave[i][0] = this->t1_meam[elti];
           t_ave[i][1] = this->t2_meam[elti];
diff --git a/src/USER-MEAMC/meam_force.cpp b/src/USER-MEAMC/meam_force.cpp
index 41eee1e6c5..85314dd8a2 100644
--- a/src/USER-MEAMC/meam_force.cpp
+++ b/src/USER-MEAMC/meam_force.cpp
@@ -245,31 +245,19 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int
 
         if (this->ialloy == 1) {
 
-          a1i = 0.0;
-          a1j = 0.0;
-          a2i = 0.0;
-          a2j = 0.0;
-          a3i = 0.0;
-          a3j = 0.0;
-          if (!iszero(tsq_ave[i][0]))
-            a1i = drhoa0j * sij / tsq_ave[i][0];
-          if (!iszero(tsq_ave[j][0]))
-            a1j = drhoa0i * sij / tsq_ave[j][0];
-          if (!iszero(tsq_ave[i][1]))
-            a2i = drhoa0j * sij / tsq_ave[i][1];
-          if (!iszero(tsq_ave[j][1]))
-            a2j = drhoa0i * sij / tsq_ave[j][1];
-          if (!iszero(tsq_ave[i][2]))
-            a3i = drhoa0j * sij / tsq_ave[i][2];
-          if (!iszero(tsq_ave[j][2]))
-            a3j = drhoa0i * sij / tsq_ave[j][2];
+          a1i = fdiv_zero(drhoa0j * sij, tsq_ave[i][0]);
+          a1j = fdiv_zero(drhoa0i * sij, tsq_ave[j][0]);
+          a2i = fdiv_zero(drhoa0j * sij, tsq_ave[i][1]);
+          a2j = fdiv_zero(drhoa0i * sij, tsq_ave[j][1]);
+          a3i = fdiv_zero(drhoa0j * sij, tsq_ave[i][2]);
+          a3j = fdiv_zero(drhoa0i * sij, tsq_ave[j][2]);
 
-          dt1dr1 = a1i * (t1mj - t1i * t1mj*t1mj);
-          dt1dr2 = a1j * (t1mi - t1j * t1mi*t1mi);
-          dt2dr1 = a2i * (t2mj - t2i * t2mj*t2mj);
-          dt2dr2 = a2j * (t2mi - t2j * t2mi*t2mi);
-          dt3dr1 = a3i * (t3mj - t3i * t3mj*t3mj);
-          dt3dr2 = a3j * (t3mi - t3j * t3mi*t3mi);
+          dt1dr1 = a1i * (t1mj - t1i * MathSpecial::square(t1mj));
+          dt1dr2 = a1j * (t1mi - t1j * MathSpecial::square(t1mi));
+          dt2dr1 = a2i * (t2mj - t2i * MathSpecial::square(t2mj));
+          dt2dr2 = a2j * (t2mi - t2j * MathSpecial::square(t2mi));
+          dt3dr1 = a3i * (t3mj - t3i * MathSpecial::square(t3mj));
+          dt3dr2 = a3j * (t3mi - t3j * MathSpecial::square(t3mi));
 
         } else if (this->ialloy == 2) {
 
@@ -331,25 +319,12 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int
           drho3ds2 = a3 * rhoa3i * arg1j3 - a3a * rhoa3i * arg3j3;
 
           if (this->ialloy == 1) {
-
-            a1i = 0.0;
-            a1j = 0.0;
-            a2i = 0.0;
-            a2j = 0.0;
-            a3i = 0.0;
-            a3j = 0.0;
-            if (!iszero(tsq_ave[i][0]))
-              a1i = rhoa0j / tsq_ave[i][0];
-            if (!iszero(tsq_ave[j][0]))
-              a1j = rhoa0i / tsq_ave[j][0];
-            if (!iszero(tsq_ave[i][1]))
-              a2i = rhoa0j / tsq_ave[i][1];
-            if (!iszero(tsq_ave[j][1]))
-              a2j = rhoa0i / tsq_ave[j][1];
-            if (!iszero(tsq_ave[i][2]))
-              a3i = rhoa0j / tsq_ave[i][2];
-            if (!iszero(tsq_ave[j][2]))
-              a3j = rhoa0i / tsq_ave[j][2];
+            a1i = fdiv_zero(rhoa0j, tsq_ave[i][0]);
+            a1j = fdiv_zero(rhoa0i, tsq_ave[j][0]);
+            a2i = fdiv_zero(rhoa0j, tsq_ave[i][1]);
+            a2j = fdiv_zero(rhoa0i, tsq_ave[j][1]);
+            a3i = fdiv_zero(rhoa0j, tsq_ave[i][2]);
+            a3j = fdiv_zero(rhoa0i, tsq_ave[j][2]);
 
             dt1ds1 = a1i * (t1mj - t1i * MathSpecial::square(t1mj));
             dt1ds2 = a1j * (t1mi - t1j * MathSpecial::square(t1mi));