diff --git a/src/LEPTON/fix_efield_lepton.cpp b/src/LEPTON/fix_efield_lepton.cpp
index 05912382a33..a055c2b15ae 100644
--- a/src/LEPTON/fix_efield_lepton.cpp
+++ b/src/LEPTON/fix_efield_lepton.cpp
@@ -179,37 +179,37 @@ void FixEfieldLepton::post_force(int vflag)
   auto dphi_x = parsed.differentiate("x").createCompiledExpression();
   auto dphi_y = parsed.differentiate("y").createCompiledExpression();
   auto dphi_z = parsed.differentiate("z").createCompiledExpression();
-  std::vector<Lepton::CompiledExpression*> dphis = {&dphi_x, &dphi_y, &dphi_z};
+  std::array<Lepton::CompiledExpression*, 3> dphis = {&dphi_x, &dphi_y, &dphi_z};
 
-  // check if reference to x, y, z exist
-  const std::array<std::string, 3> variableNames = {"x", "y", "z"};
-  std::array<bool, 3> phi_has_ref = {true, true, true};
+  // array of vectors of ptrs to Lepton variable references
+  std::array<std::vector<double *>, 3> var_ref_ptrs{};
+
+  // fill ptr-vectors with Lepton refs as needed
+  const char* DIM_NAMES[] = {"x", "y", "z"};
   if (atom->q_flag){
     phi = parsed.createCompiledExpression();
-    for (size_t i = 0; i < 3; i++) {
+    for (size_t d = 0; d < 3; d++) {
       try {
-        phi.getVariableReference(variableNames[i]);
-      }
-      catch (Lepton::Exception &) {
-        phi_has_ref[i] = false;
+        double *ptr = &(phi.getVariableReference(DIM_NAMES[d]));
+        var_ref_ptrs[d].push_back(ptr);
+      } catch (Lepton::Exception &) {
+        // do nothing
       }
     }
   }
-  std::vector<std::array<bool, 3>> dphis_has_ref;
+
   bool e_uniform = true;
-  for (auto &dphi : dphis) {
-    dphis_has_ref.push_back({false, false, false});
-    for (size_t i = 0; i < 3; i++) {
+  for (size_t j = 0; j < 3; j++)
+    for (size_t d = 0; d < 3; d++) {
       try {
-        (*dphi).getVariableReference(variableNames[i]);
-        dphis_has_ref.back()[i] = true;
+        double *ptr = &((*dphis[j]).getVariableReference(DIM_NAMES[d]));
+        var_ref_ptrs[d].push_back(ptr);
         e_uniform = false;
       }
       catch (Lepton::Exception &) {
         // do nothing
       }
     }
-  }
   if (!e_uniform && atom->mu_flag && h < 0) {
     error->all(FLERR, "Fix {} requires keyword `step' for dipoles in a non-uniform electric field", style);
   }
@@ -227,43 +227,45 @@ void FixEfieldLepton::post_force(int vflag)
 
   double ex, ey, ez;
   double fx, fy, fz;
-  double v[6], unwrap[3];
+  double v[6], unwrap[3], dstep[3];
   double xf, yf, zf, xb, yb, zb;
   double exf, eyf, ezf, exb, eyb, ezb;
   double mu_norm, h_mu;
 
-  if (atom->q_flag && atom->mu_flag) {
-    double *q = atom->q;
-    double **mu = atom->mu;
-    double **t = atom->torque;
-    for (int i = 0; i < nlocal; i++) {
-      if (mask[i] & groupbit) {
-        if (region && !region->match(x[i][0], x[i][1], x[i][2])) continue;
-        domain->unmap(x[i], image[i], unwrap);
-
-        // evaluate e-field, used by q and mu
-        for (size_t j = 0; j < 3; j++) {
-          if (dphis_has_ref[j][0]) (*dphis[j]).getVariableReference("x") = unwrap[0];
-          if (dphis_has_ref[j][1]) (*dphis[j]).getVariableReference("y") = unwrap[1];
-          if (dphis_has_ref[j][2]) (*dphis[j]).getVariableReference("z") = unwrap[2];
+  double *q = atom->q;
+  double **mu = atom->mu;
+  double **t = atom->torque;
+
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) {
+      if (region && !region->match(x[i][0], x[i][1], x[i][2])) continue;
+      fx = fy = fz = 0.0;
+      domain->unmap(x[i], image[i], unwrap);
+
+      // put unwrapped coords into Lepton variable refs
+      for (size_t d = 0; d < 3; d++) {
+        for (auto & var_ref_ptr : var_ref_ptrs[d]) {
+          *var_ref_ptr = unwrap[d];
         }
-        ex = -dphi_x.evaluate();
-        ey = -dphi_y.evaluate();
-        ez = -dphi_z.evaluate();
+      }
 
-        if (phi_has_ref[0]) phi.getVariableReference("x") = unwrap[0];
-        if (phi_has_ref[1]) phi.getVariableReference("y") = unwrap[1];
-        if (phi_has_ref[2]) phi.getVariableReference("z") = unwrap[2];
+      // evaluate e-field, used by q and mu
+      ex = -dphi_x.evaluate();
+      ey = -dphi_y.evaluate();
+      ez = -dphi_z.evaluate();
 
-        // charges
-        // force = q E
+      // charges
+      // force = q E
+      if (atom->q_flag) {
         fx = qe2f * q[i] * ex;
         fy = qe2f * q[i] * ey;
         fz = qe2f * q[i] * ez;
         // potential energy = q phi
         fsum[0] += qe2f * q[i] * phi.evaluate();
+      }
 
-        // dipoles
+      if (atom->mu_flag) {
+      // dipoles
         mu_norm = sqrt(mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1] + mu[i][2]*mu[i][2]);
         if (mu_norm > EPSILON) {
           // torque = mu cross E
@@ -273,181 +275,58 @@ void FixEfieldLepton::post_force(int vflag)
           // potential energy = - mu dot E
           fsum[0] -= mue2e * (mu[i][0] * ex + mu[i][1] * ey + mu[i][2] * ez);
 
-          // force = (mu dot D) E
+          // force = (mu dot D) E for non-uniform E
           // using central difference method
-          h_mu = h / mu_norm;
-
-          xf = unwrap[0] + h_mu * mu[i][0];
-          yf = unwrap[1] + h_mu * mu[i][1];
-          zf = unwrap[2] + h_mu * mu[i][2];
-          for (size_t j = 0; j < 3; j++) {
-            if (dphis_has_ref[j][0]) (*dphis[j]).getVariableReference("x") = xf;
-            if (dphis_has_ref[j][1]) (*dphis[j]).getVariableReference("y") = yf;
-            if (dphis_has_ref[j][2]) (*dphis[j]).getVariableReference("z") = zf;
-          }
-          exf = -dphi_x.evaluate();
-          eyf = -dphi_y.evaluate();
-          ezf = -dphi_z.evaluate();
-
-          xb = unwrap[0] - h_mu * mu[i][0];
-          yb = unwrap[1] - h_mu * mu[i][1];
-          zb = unwrap[2] - h_mu * mu[i][2];
-          for (size_t j = 0; j < 3; j++) {
-            if (dphis_has_ref[j][0]) (*dphis[j]).getVariableReference("x") = xb;
-            if (dphis_has_ref[j][1]) (*dphis[j]).getVariableReference("y") = yb;
-            if (dphis_has_ref[j][2]) (*dphis[j]).getVariableReference("z") = zb;
+          if (!e_uniform) {
+            h_mu = h / mu_norm;
+            dstep[0] = h_mu * mu[i][0];
+            dstep[1] = h_mu * mu[i][1];
+            dstep[2] = h_mu * mu[i][2];
+
+            // one step forwards, two steps back ;)
+            for (size_t d = 0; d < 3; d++) {
+              for (auto & var_ref_ptr : var_ref_ptrs[d]) {
+                *var_ref_ptr += dstep[d];
+              }
+            }
+
+            exf = -dphi_x.evaluate();
+            eyf = -dphi_y.evaluate();
+            ezf = -dphi_z.evaluate();
+
+            for (size_t d = 0; d < 3; d++) {
+              for (auto & var_ref_ptr : var_ref_ptrs[d]) {
+                *var_ref_ptr -= 2*dstep[d];
+              }
+            }
+
+            exb = -dphi_x.evaluate();
+            eyb = -dphi_y.evaluate();
+            ezb = -dphi_z.evaluate();
+
+            fx += qe2f * (exf - exb) / 2.0 / h_mu;
+            fy += qe2f * (eyf - eyb) / 2.0 / h_mu;
+            fz += qe2f * (ezf - ezb) / 2.0 / h_mu;
           }
-          exb = -dphi_x.evaluate();
-          eyb = -dphi_y.evaluate();
-          ezb = -dphi_z.evaluate();
-
-          fx += qe2f * (exf - exb) / 2.0 / h_mu;
-          fy += qe2f * (eyf - eyb) / 2.0 / h_mu;
-          fz += qe2f * (ezf - ezb) / 2.0 / h_mu;
-        }
-
-        f[i][0] += fx;
-        f[i][1] += fy;
-        f[i][2] += fz;
-
-        fsum[1] += fx;
-        fsum[2] += fy;
-        fsum[3] += fz;
-
-        if (evflag) {
-          v[0] = fx * unwrap[0];
-          v[1] = fy * unwrap[1];
-          v[2] = fz * unwrap[2];
-          v[3] = fx * unwrap[1];
-          v[4] = fx * unwrap[2];
-          v[5] = fy * unwrap[2];
-          v_tally(i, v);
-        }
-      }
-    }
-  } else if (atom->q_flag && !atom->mu_flag) {
-    double *q = atom->q;
-    for (int i = 0; i < nlocal; i++) {
-      if (mask[i] & groupbit) {
-        if (region && !region->match(x[i][0], x[i][1], x[i][2])) continue;
-        domain->unmap(x[i], image[i], unwrap);
-
-        for (size_t j = 0; j < 3; j++) {
-          if (dphis_has_ref[j][0]) (*dphis[j]).getVariableReference("x") = unwrap[0];
-          if (dphis_has_ref[j][1]) (*dphis[j]).getVariableReference("y") = unwrap[1];
-          if (dphis_has_ref[j][2]) (*dphis[j]).getVariableReference("z") = unwrap[2];
-        }
-        ex = -dphi_x.evaluate();
-        ey = -dphi_y.evaluate();
-        ez = -dphi_z.evaluate();
-
-        if (phi_has_ref[0]) phi.getVariableReference("x") = unwrap[0];
-        if (phi_has_ref[1]) phi.getVariableReference("y") = unwrap[1];
-        if (phi_has_ref[2]) phi.getVariableReference("z") = unwrap[2];
-
-        // force = q E
-        fx = qe2f * q[i] * ex;
-        fy = qe2f * q[i] * ey;
-        fz = qe2f * q[i] * ez;
-        // potential energy = q phi
-        fsum[0] += qe2f * q[i] * phi.evaluate();
-
-        f[i][0] += fx;
-        f[i][1] += fy;
-        f[i][2] += fz;
-
-        fsum[1] += fx;
-        fsum[2] += fy;
-        fsum[3] += fz;
-
-        if (evflag) {
-          v[0] = fx * unwrap[0];
-          v[1] = fy * unwrap[1];
-          v[2] = fz * unwrap[2];
-          v[3] = fx * unwrap[1];
-          v[4] = fx * unwrap[2];
-          v[5] = fy * unwrap[2];
-          v_tally(i, v);
         }
       }
-    }
-  } else if (!atom->q_flag && atom->mu_flag) {
-    double **mu = atom->mu;
-    double **t = atom->torque;
-    for (int i = 0; i < nlocal; i++) {
-      if (mask[i] & groupbit) {
-        if (region && !region->match(x[i][0], x[i][1], x[i][2])) continue;
-
-        mu_norm = sqrt(mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1] + mu[i][2]*mu[i][2]);
-        if (mu_norm > EPSILON) continue;
-
-        domain->unmap(x[i], image[i], unwrap);
 
-        for (size_t j = 0; j < 3; j++) {
-          if (dphis_has_ref[j][0]) (*dphis[j]).getVariableReference("x") = unwrap[0];
-          if (dphis_has_ref[j][1]) (*dphis[j]).getVariableReference("y") = unwrap[1];
-          if (dphis_has_ref[j][2]) (*dphis[j]).getVariableReference("z") = unwrap[2];
-        }
-        ex = -dphi_x.evaluate();
-        ey = -dphi_y.evaluate();
-        ez = -dphi_z.evaluate();
-
-        // torque = mu cross E
-        t[i][0] += mue2e * (ez * mu[i][1] - ey * mu[i][2]);
-        t[i][1] += mue2e * (ex * mu[i][2] - ez * mu[i][0]);
-        t[i][2] += mue2e * (ey * mu[i][0] - ex * mu[i][1]);
-        // potential energy = - mu dot E
-        fsum[0] -= mue2e * (mu[i][0] * ex + mu[i][1] * ey + mu[i][2] * ez);
-
-        // force = (mu dot D) E
-        // using central difference method
-        h_mu = h / sqrt(mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1] + mu[i][2]*mu[i][2]);
-
-        xf = unwrap[0] + h_mu * mu[i][0];
-        yf = unwrap[1] + h_mu * mu[i][1];
-        zf = unwrap[2] + h_mu * mu[i][2];
-        for (size_t j = 0; j < 3; j++) {
-          if (dphis_has_ref[j][0]) (*dphis[j]).getVariableReference("x") = xf;
-          if (dphis_has_ref[j][1]) (*dphis[j]).getVariableReference("y") = yf;
-          if (dphis_has_ref[j][2]) (*dphis[j]).getVariableReference("z") = zf;
-        }
-        exf = -dphi_x.evaluate();
-        eyf = -dphi_y.evaluate();
-        ezf = -dphi_z.evaluate();
-
-        xb = unwrap[0] - h_mu * mu[i][0];
-        yb = unwrap[1] - h_mu * mu[i][1];
-        zb = unwrap[2] - h_mu * mu[i][2];
-        for (size_t j = 0; j < 3; j++) {
-          if (dphis_has_ref[j][0]) (*dphis[j]).getVariableReference("x") = xb;
-          if (dphis_has_ref[j][1]) (*dphis[j]).getVariableReference("y") = yb;
-          if (dphis_has_ref[j][2]) (*dphis[j]).getVariableReference("z") = zb;
-        }
-        exb = -dphi_x.evaluate();
-        eyb = -dphi_y.evaluate();
-        ezb = -dphi_z.evaluate();
-
-        fx = qe2f * (exf - exb) / 2.0 / h_mu;
-        fy = qe2f * (eyf - eyb) / 2.0 / h_mu;
-        fz = qe2f * (ezf - ezb) / 2.0 / h_mu;
-
-        f[i][0] += fx;
-        f[i][1] += fy;
-        f[i][2] += fz;
-
-        fsum[1] += fx;
-        fsum[2] += fy;
-        fsum[3] += fz;
-
-        if (evflag) {
-          v[0] = fx * unwrap[0];
-          v[1] = fy * unwrap[1];
-          v[2] = fz * unwrap[2];
-          v[3] = fx * unwrap[1];
-          v[4] = fx * unwrap[2];
-          v[5] = fy * unwrap[2];
-          v_tally(i, v);
-        }
+      f[i][0] += fx;
+      f[i][1] += fy;
+      f[i][2] += fz;
+
+      fsum[1] += fx;
+      fsum[2] += fy;
+      fsum[3] += fz;
+
+      if (evflag) {
+        v[0] = fx * unwrap[0];
+        v[1] = fy * unwrap[1];
+        v[2] = fz * unwrap[2];
+        v[3] = fx * unwrap[1];
+        v[4] = fx * unwrap[2];
+        v[5] = fy * unwrap[2];
+        v_tally(i, v);
       }
     }
   }