diff --git a/src/create_atoms.cpp b/src/create_atoms.cpp
index b1229f6a68..2f128470ab 100644
--- a/src/create_atoms.cpp
+++ b/src/create_atoms.cpp
@@ -271,6 +271,14 @@ void CreateAtoms::add_many()
   //     and is shifted away from true boundary
   //     which is where atoms are likely to be generated
 
+  double epsilon[3];
+  if (triclinic) epsilon[0] = epsilon[1] = epsilon[2] = EPSILON;
+  else {
+    epsilon[0] = domain->prd[0] * EPSILON;
+    epsilon[1] = domain->prd[1] * EPSILON;
+    epsilon[2] = domain->prd[2] * EPSILON;
+  }
+
   double sublo[3],subhi[3];
 
   if (triclinic == 0) {
@@ -284,16 +292,16 @@ void CreateAtoms::add_many()
   }
 
   if (domain->xperiodic) {
-    if (comm->myloc[0] == 0) sublo[0] -= EPSILON;
-    if (comm->myloc[0] == comm->procgrid[0]-1) subhi[0] -= 2.0*EPSILON;
+    if (comm->myloc[0] == 0) sublo[0] -= epsilon[0];
+    if (comm->myloc[0] == comm->procgrid[0]-1) subhi[0] -= 2.0*epsilon[0];
   }
   if (domain->yperiodic) {
-    if (comm->myloc[1] == 0) sublo[1] -= EPSILON;
-    if (comm->myloc[1] == comm->procgrid[1]-1) subhi[1] -= 2.0*EPSILON;
+    if (comm->myloc[1] == 0) sublo[1] -= epsilon[1];
+    if (comm->myloc[1] == comm->procgrid[1]-1) subhi[1] -= 2.0*epsilon[1];
   }
   if (domain->zperiodic) {
-    if (comm->myloc[2] == 0) sublo[2] -= EPSILON;
-    if (comm->myloc[2] == comm->procgrid[2]-1) subhi[2] -= 2.0*EPSILON;
+    if (comm->myloc[2] == 0) sublo[2] -= epsilon[2];
+    if (comm->myloc[2] == comm->procgrid[2]-1) subhi[2] -= 2.0*epsilon[2];
   }
 
   // iterate on 3d periodic lattice of unit cells using loop bounds