diff --git a/include/openmc/cell.h b/include/openmc/cell.h
index 3008948168..d2365c84c4 100644
--- a/include/openmc/cell.h
+++ b/include/openmc/cell.h
@@ -381,6 +381,9 @@ public:
   std::pair<double, int32_t> distance(Position r, Direction u,
     int32_t on_surface, GeometryState* p) const override;
 
+  std::pair<double, int32_t> distance_in_virtual_lattice(
+    Position r, Direction u, int32_t on_surface, GeometryState* p) const;
+
   bool contains(Position r, Direction u, int32_t on_surface) const override
   {
     return region_.contains(r, u, on_surface);
diff --git a/include/openmc/geometry.h b/include/openmc/geometry.h
index 107cc7d1f3..d41a98e301 100644
--- a/include/openmc/geometry.h
+++ b/include/openmc/geometry.h
@@ -64,6 +64,10 @@ bool exhaustive_find_cell(GeometryState& p, bool verbose = false);
 bool neighbor_list_find_cell(
   GeometryState& p, bool verbose = false); // Only usable on surface crossings
 
+bool find_cell_in_virtual_lattice(GeometryState& p,
+  bool verbose =
+    false); // Only usable on triso surface crossings in virtual lattice
+
 //==============================================================================
 //! Move a particle into a new lattice tile.
 //==============================================================================
diff --git a/include/openmc/universe.h b/include/openmc/universe.h
index 2afc4b97ac..1c9fcceabf 100644
--- a/include/openmc/universe.h
+++ b/include/openmc/universe.h
@@ -38,6 +38,8 @@ public:
 
   virtual bool find_cell(GeometryState& p) const;
 
+  virtual bool find_cell_in_virtual_lattice(GeometryState& p) const;
+
   BoundingBox bounding_box() const;
 
   /* By default, universes are CSG universes. The DAGMC
diff --git a/src/cell.cpp b/src/cell.cpp
index a5d68c8c56..f7d077fb51 100644
--- a/src/cell.cpp
+++ b/src/cell.cpp
@@ -505,96 +505,101 @@ vector<int32_t>::iterator CSGCell::find_left_parenthesis(
 }
 std::pair<double, int32_t> CSGCell::distance(
   Position r, Direction u, int32_t on_surface, GeometryState* p) const
+{
+  if (virtual_lattice_) {
+    return distance_in_virtual_lattice(r, u, on_surface, p);
+  } else {
+    return region_.distance(r, u, on_surface);
+  }
+}
+
+std::pair<double, int32_t> CSGCell::distance_in_virtual_lattice(
+  Position r, Direction u, int32_t on_surface, GeometryState* p) const
 {
   double min_dist {INFTY};
   int32_t i_surf {std::numeric_limits<int32_t>::max()};
   double min_dis_vl;
   int32_t i_surf_vl;
-  if (virtual_lattice_) {
-    double max_dis = p->collision_distance();
-    double tol_dis = 0;
-    vector<double> dis_to_bou(3), dis_to_bou_max(3);
-    double u_value = sqrt(pow(u.x, 2) + pow(u.y, 2) +
-                          pow(u.z, 2)); // don't know if u has been normalized
-    vector<double> norm_u = {u.x / u_value, u.y / u_value, u.z / u_value};
-    vector<int> lat_ind(3);
-    vector<double> temp_pos = {r.x, r.y, r.z};
-    int loop_time;
-    for (int i = 0; i < 3; i++) {
-      lat_ind[i] = floor((temp_pos[i] - vl_lower_left_[i]) / vl_pitch_[i]);
-      if (lat_ind[i] == vl_shape_[i] && norm_u[i] < 0) {
-        lat_ind[i] = vl_shape_[i] - 1;
-      }
-      if (lat_ind[i] == -1 && norm_u[i] > 0) {
-        lat_ind[i] = 0;
-      }
+
+  double max_dis = p->collision_distance();
+  double tol_dis = 0;
+  vector<double> dis_to_bou(3), dis_to_bou_max(3);
+  double u_value = sqrt(pow(u.x, 2) + pow(u.y, 2) +
+                        pow(u.z, 2)); // don't know if u has been normalized
+  vector<double> norm_u = {u.x / u_value, u.y / u_value, u.z / u_value};
+  vector<int> lat_ind(3);
+  vector<double> temp_pos = {r.x, r.y, r.z};
+  int loop_time;
+  for (int i = 0; i < 3; i++) {
+    lat_ind[i] = floor((temp_pos[i] - vl_lower_left_[i]) / vl_pitch_[i]);
+    if (lat_ind[i] == vl_shape_[i] && norm_u[i] < 0) {
+      lat_ind[i] = vl_shape_[i] - 1;
     }
-
-    dis_to_bou = {INFTY, INFTY, INFTY};
-    for (int i = 0; i < 3; i++) {
-      if (norm_u[i] > 0) {
-        dis_to_bou[i] = std::abs(
-          ((lat_ind[i] + 1) * vl_pitch_[i] + vl_lower_left_[i] - temp_pos[i]) /
-          norm_u[i]);
-        dis_to_bou_max[i] = vl_pitch_[i] / norm_u[i];
-      } else if (norm_u[i] < 0) {
-        dis_to_bou[i] = std::abs(
-          (lat_ind[i] * vl_pitch_[i] + vl_lower_left_[i] - temp_pos[i]) /
-          norm_u[i]);
-        dis_to_bou_max[i] = -vl_pitch_[i] / norm_u[i];
-      }
+    if (lat_ind[i] == -1 && norm_u[i] > 0) {
+      lat_ind[i] = 0;
     }
-
-    while (true) {
-      if (lat_ind[0] < 0 || lat_ind[0] >= vl_shape_[0] || lat_ind[1] < 0 ||
-          lat_ind[1] >= vl_shape_[1] || lat_ind[2] < 0 ||
-          lat_ind[2] >= vl_shape_[2])
-        break;
-
-      for (int token :
-        vl_triso_distribution_[lat_ind[0] + lat_ind[1] * vl_shape_[0] +
-                               lat_ind[2] * vl_shape_[0] * vl_shape_[1]]) {
-        bool coincident {std::abs(token) == std::abs(on_surface)};
-        double d {model::surfaces[abs(token) - 1]->distance(r, u, coincident)};
-        if (d < min_dist) {
-          if (min_dist - d >= FP_PRECISION * min_dist) {
-            min_dist = d;
-            i_surf = -token;
-          }
-        }
-      }
-
-      int mes_bou_crossed = 0;
-      if (dis_to_bou[1] < dis_to_bou[0]) {
-        mes_bou_crossed = 1;
-      }
-      if (dis_to_bou[2] < dis_to_bou[mes_bou_crossed]) {
-        mes_bou_crossed = 2;
-      }
-
-      tol_dis = dis_to_bou[mes_bou_crossed];
-
-      if (min_dist < tol_dis) {
-        break;
-      }
-
-      if (norm_u[mes_bou_crossed] > 0) {
-        lat_ind[mes_bou_crossed] += 1;
-      } else {
-        lat_ind[mes_bou_crossed] += -1;
-      }
-
-      dis_to_bou[mes_bou_crossed] += dis_to_bou_max[mes_bou_crossed];
-
-      if (tol_dis > max_dis) {
-        break;
-      }
-    }
-
-  } else {
-    return region_.distance(r, u, on_surface);
   }
 
+  dis_to_bou = {INFTY, INFTY, INFTY};
+  for (int i = 0; i < 3; i++) {
+    if (norm_u[i] > 0) {
+      dis_to_bou[i] = std::abs(
+        ((lat_ind[i] + 1) * vl_pitch_[i] + vl_lower_left_[i] - temp_pos[i]) /
+        norm_u[i]);
+      dis_to_bou_max[i] = vl_pitch_[i] / norm_u[i];
+    } else if (norm_u[i] < 0) {
+      dis_to_bou[i] =
+        std::abs((lat_ind[i] * vl_pitch_[i] + vl_lower_left_[i] - temp_pos[i]) /
+                 norm_u[i]);
+      dis_to_bou_max[i] = -vl_pitch_[i] / norm_u[i];
+    }
+  }
+
+  while (true) {
+    if (lat_ind[0] < 0 || lat_ind[0] >= vl_shape_[0] || lat_ind[1] < 0 ||
+        lat_ind[1] >= vl_shape_[1] || lat_ind[2] < 0 ||
+        lat_ind[2] >= vl_shape_[2])
+      break;
+
+    for (int token :
+      vl_triso_distribution_[lat_ind[0] + lat_ind[1] * vl_shape_[0] +
+                             lat_ind[2] * vl_shape_[0] * vl_shape_[1]]) {
+      bool coincident {std::abs(token) == std::abs(on_surface)};
+      double d {model::surfaces[abs(token) - 1]->distance(r, u, coincident)};
+      if (d < min_dist) {
+        if (min_dist - d >= FP_PRECISION * min_dist) {
+          min_dist = d;
+          i_surf = -token;
+        }
+      }
+    }
+
+    int mes_bou_crossed = 0;
+    if (dis_to_bou[1] < dis_to_bou[0]) {
+      mes_bou_crossed = 1;
+    }
+    if (dis_to_bou[2] < dis_to_bou[mes_bou_crossed]) {
+      mes_bou_crossed = 2;
+    }
+
+    tol_dis = dis_to_bou[mes_bou_crossed];
+
+    if (min_dist < tol_dis) {
+      break;
+    }
+
+    if (norm_u[mes_bou_crossed] > 0) {
+      lat_ind[mes_bou_crossed] += 1;
+    } else {
+      lat_ind[mes_bou_crossed] += -1;
+    }
+
+    dis_to_bou[mes_bou_crossed] += dis_to_bou_max[mes_bou_crossed];
+
+    if (tol_dis > max_dis) {
+      break;
+    }
+  }
   return {min_dist, i_surf};
 }
 
diff --git a/src/geometry.cpp b/src/geometry.cpp
index df48508916..f32b63ed7a 100644
--- a/src/geometry.cpp
+++ b/src/geometry.cpp
@@ -294,6 +294,98 @@ bool exhaustive_find_cell(GeometryState& p, bool verbose)
   return find_cell_inner(p, nullptr, verbose);
 }
 
+bool find_cell_in_virtual_lattice(GeometryState& p, bool verbose)
+{
+  int i_surface = std::abs(p.surface());
+  if (p.surface() > 0) {
+    for (int i = p.n_coord(); i < model::n_coord_levels; i++) {
+      p.coord(i).reset();
+    }
+    p.coord(p.n_coord() - 1).cell() =
+      model::cell_map[model::surfaces[i_surface - 1]->triso_base_index_];
+  } else if (p.surface() < 0) {
+    for (int i = p.n_coord(); i < model::n_coord_levels; i++) {
+      p.coord(i).reset();
+    }
+    if (model::surfaces[i_surface - 1]->triso_particle_index_ == -1) {
+      fatal_error(fmt::format("Particle cell of surface {} is not defined",
+        model::surfaces[i_surface - 1]->id_));
+    }
+    p.lowest_coord().cell() =
+      model::cell_map[model::surfaces[i_surface - 1]->triso_particle_index_];
+  }
+
+  // find material
+  bool found = true;
+  int i_cell = p.lowest_coord().cell();
+  for (;; ++p.n_coord()) {
+    if (i_cell == C_NONE) {
+      int i_universe = p.lowest_coord().universe();
+      const auto& univ {model::universes[i_universe]};
+
+      if (univ->filled_with_triso_base_ != -1) {
+        p.lowest_coord().cell() =
+          model::cell_map[univ->filled_with_triso_base_];
+        found = true;
+      } else {
+        found = univ->find_cell(p);
+      }
+      if (!found) {
+        return found;
+      }
+    }
+
+    i_cell = p.lowest_coord().cell();
+
+    Cell& c {*model::cells[i_cell]};
+    if (c.type_ == Fill::MATERIAL) {
+      // Found a material cell which means this is the lowest coord level.
+
+      p.cell_instance() = 0;
+      // Find the distribcell instance number.
+      if (c.distribcell_index_ >= 0) {
+        p.cell_instance() = cell_instance_at_level(p, p.n_coord() - 1);
+      }
+
+      // Set the material and temperature.
+      p.material_last() = p.material();
+      if (c.material_.size() > 1) {
+        p.material() = c.material_[p.cell_instance()];
+      } else {
+        p.material() = c.material_[0];
+      }
+      p.sqrtkT_last() = p.sqrtkT();
+      if (c.sqrtkT_.size() > 1) {
+        p.sqrtkT() = c.sqrtkT_[p.cell_instance()];
+      } else {
+        p.sqrtkT() = c.sqrtkT_[0];
+      }
+      return found;
+
+    } else if (c.type_ == Fill::UNIVERSE) {
+      //========================================================================
+      //! Found a lower universe, update this coord level then search the
+      //! next.
+
+      // Set the lower coordinate level universe.
+      auto& coor {p.coord(p.n_coord())};
+      coor.universe() = c.fill_;
+
+      // Set the position and direction.
+      coor.r() = p.r_local();
+      coor.u() = p.u_local();
+
+      // Apply translation.
+      coor.r() -= c.translation_;
+
+      // Apply rotation.
+      if (!c.rotation_.empty()) {
+        coor.rotate(c.rotation_);
+      }
+      i_cell = C_NONE;
+    }
+  }
+}
 //==============================================================================
 
 void cross_lattice(GeometryState& p, const BoundaryInfo& boundary, bool verbose)
diff --git a/src/particle.cpp b/src/particle.cpp
index dd813013a0..67cb144935 100644
--- a/src/particle.cpp
+++ b/src/particle.cpp
@@ -577,97 +577,13 @@ void Particle::cross_surface(const Surface& surf)
   int i_surface = std::abs(surface());
   bool verbose = settings::verbosity >= 10 || trace();
   if (surf.is_triso_surface_) {
-    if (surface() > 0) {
-      for (int i = n_coord(); i < model::n_coord_levels; i++) {
-        coord(i).reset();
-      }
-      coord(n_coord() - 1).cell() =
-        model::cell_map[model::surfaces[i_surface - 1]->triso_base_index_];
-    } else if (surface() < 0) {
-      for (int i = n_coord(); i < model::n_coord_levels; i++) {
-        coord(i).reset();
-      }
-      if (model::surfaces[i_surface - 1]->triso_particle_index_ == -1) {
-        fatal_error(fmt::format("Particle cell of surface {} is not defined",
-          model::surfaces[i_surface - 1]->id_));
-      }
-      coord(n_coord() - 1).cell() =
-        model::cell_map[model::surfaces[i_surface - 1]->triso_particle_index_];
-    }
-
-    // find material
-    bool found = true;
-    int i_cell = coord(n_coord() - 1).cell();
-    for (;; ++n_coord()) {
-      if (i_cell == C_NONE) {
-        int i_universe = coord(n_coord() - 1).universe();
-        const auto& univ {model::universes[i_universe]};
-
-        if (univ->filled_with_triso_base_ != -1) {
-          coord(n_coord() - 1).cell() =
-            model::cell_map[univ->filled_with_triso_base_];
-          found = true;
-        } else {
-          found = univ->find_cell(*this);
-        }
-        if (!found) {
-          break;
-        }
-      }
-
-      i_cell = coord(n_coord() - 1).cell();
-
-      Cell& c {*model::cells[i_cell]};
-      if (c.type_ == Fill::MATERIAL) {
-        // Found a material cell which means this is the lowest coord level.
-
-        cell_instance() = 0;
-        // Find the distribcell instance number.
-        if (c.distribcell_index_ >= 0) {
-          cell_instance() = cell_instance_at_level(*this, n_coord() - 1);
-        }
-
-        // Set the material and temperature.
-        material_last() = material();
-        if (c.material_.size() > 1) {
-          material() = c.material_[cell_instance()];
-        } else {
-          material() = c.material_[0];
-        }
-        sqrtkT_last() = sqrtkT();
-        if (c.sqrtkT_.size() > 1) {
-          sqrtkT() = c.sqrtkT_[cell_instance()];
-        } else {
-          sqrtkT() = c.sqrtkT_[0];
-        }
-        return;
-
-      } else if (c.type_ == Fill::UNIVERSE) {
-        //========================================================================
-        //! Found a lower universe, update this coord level then search the
-        //! next.
-
-        // Set the lower coordinate level universe.
-        auto& coor {coord(n_coord())};
-        coor.universe() = c.fill_;
-
-        // Set the position and direction.
-        coor.r() = r_local();
-        coor.u() = u_local();
-
-        // Apply translation.
-        coor.r() -= c.translation_;
-
-        // Apply rotation.
-        if (!c.rotation_.empty()) {
-          coor.rotate(c.rotation_);
-        }
-        i_cell = C_NONE;
-      }
+    if (find_cell_in_virtual_lattice(*this, verbose)) {
+      return;
     }
   } else {
-    if (neighbor_list_find_cell(*this, verbose))
+    if (neighbor_list_find_cell(*this, verbose)) {
       return;
+    }
   }
 
   // ==========================================================================
diff --git a/src/universe.cpp b/src/universe.cpp
index 5fa2892b27..054f57d060 100644
--- a/src/universe.cpp
+++ b/src/universe.cpp
@@ -40,59 +40,9 @@ void Universe::to_hdf5(hid_t universes_group) const
 bool Universe::find_cell(GeometryState& p) const
 {
   if (filled_with_triso_base_ != -1) {
-    Cell& c {*model::cells[model::cell_map[filled_with_triso_base_]]};
-    vector<int> lat_ind(3);
-    Position r {p.r_local()};
-    lat_ind[0] = std::max<int>(
-      std::min<int>(floor((r.x - c.vl_lower_left_[0]) / c.vl_pitch_[0]),
-        c.vl_shape_[0] - 1),
-      0);
-    lat_ind[1] = std::max<int>(
-      std::min<int>(floor((r.y - c.vl_lower_left_[1]) / c.vl_pitch_[1]),
-        c.vl_shape_[1] - 1),
-      0);
-    lat_ind[2] = std::max<int>(
-      std::min<int>(floor((r.z - c.vl_lower_left_[2]) / c.vl_pitch_[2]),
-        c.vl_shape_[2] - 1),
-      0);
-
-    int32_t i_univ = p.coord(p.n_coord() - 1).universe();
-    for (int token :
-      c.vl_triso_distribution_[lat_ind[0] + lat_ind[1] * c.vl_shape_[0] +
-                               lat_ind[2] * c.vl_shape_[0] * c.vl_shape_[1]]) {
-      vector<double> triso_center =
-        model::surfaces[abs(token) - 1]->get_center();
-      double triso_radius = model::surfaces[abs(token) - 1]->get_radius();
-      if (model::cells[model::cell_map[model::surfaces[abs(token) - 1]
-                                         ->triso_base_index_]]
-            ->universe_ != i_univ)
-        continue;
-      if (abs(token) == abs(p.surface())) {
-        if (p.surface() < 0) {
-          p.coord(p.n_coord() - 1).cell() =
-            model::cell_map[model::surfaces[abs(token) - 1]
-                              ->triso_particle_index_];
-          return true;
-        } else {
-          p.coord(p.n_coord() - 1).cell() =
-            model::cell_map[filled_with_triso_base_];
-          return true;
-        }
-      }
-      if (pow(r.x - triso_center[0], 2) + pow(r.y - triso_center[1], 2) +
-            pow(r.z - triso_center[2], 2) <
-          pow(triso_radius, 2)) {
-        p.coord(p.n_coord() - 1).cell() =
-          model::cell_map[model::surfaces[abs(token) - 1]
-                            ->triso_particle_index_];
-        return true;
-      }
-    }
-    if (model::cells[model::cell_map[filled_with_triso_base_]]->universe_ ==
-        i_univ) {
-      p.coord(p.n_coord() - 1).cell() =
-        model::cell_map[filled_with_triso_base_];
-      return true;
+    bool found = find_cell_in_virtual_lattice(p);
+    if (found) {
+      return found;
     }
   }
   const auto& cells {
@@ -114,6 +64,60 @@ bool Universe::find_cell(GeometryState& p) const
   }
   return false;
 }
+bool Universe::find_cell_in_virtual_lattice(GeometryState& p) const
+{
+  Cell& c {*model::cells[model::cell_map[filled_with_triso_base_]]};
+  vector<int> lat_ind(3);
+  Position r {p.r_local()};
+  lat_ind[0] = std::max<int>(
+    std::min<int>(
+      floor((r.x - c.vl_lower_left_[0]) / c.vl_pitch_[0]), c.vl_shape_[0] - 1),
+    0);
+  lat_ind[1] = std::max<int>(
+    std::min<int>(
+      floor((r.y - c.vl_lower_left_[1]) / c.vl_pitch_[1]), c.vl_shape_[1] - 1),
+    0);
+  lat_ind[2] = std::max<int>(
+    std::min<int>(
+      floor((r.z - c.vl_lower_left_[2]) / c.vl_pitch_[2]), c.vl_shape_[2] - 1),
+    0);
+
+  int32_t i_univ = p.lowest_coord().universe();
+  for (int token :
+    c.vl_triso_distribution_[lat_ind[0] + lat_ind[1] * c.vl_shape_[0] +
+                             lat_ind[2] * c.vl_shape_[0] * c.vl_shape_[1]]) {
+    vector<double> triso_center = model::surfaces[abs(token) - 1]->get_center();
+    double triso_radius = model::surfaces[abs(token) - 1]->get_radius();
+    if (model::cells
+          [model::cell_map[model::surfaces[abs(token) - 1]->triso_base_index_]]
+            ->universe_ != i_univ)
+      continue;
+    if (abs(token) == abs(p.surface())) {
+      if (p.surface() < 0) {
+        p.lowest_coord().cell() =
+          model::cell_map[model::surfaces[abs(token) - 1]
+                            ->triso_particle_index_];
+        return true;
+      } else {
+        p.lowest_coord().cell() = model::cell_map[filled_with_triso_base_];
+        return true;
+      }
+    }
+    if (pow(r.x - triso_center[0], 2) + pow(r.y - triso_center[1], 2) +
+          pow(r.z - triso_center[2], 2) <
+        pow(triso_radius, 2)) {
+      p.lowest_coord().cell() =
+        model::cell_map[model::surfaces[abs(token) - 1]->triso_particle_index_];
+      return true;
+    }
+  }
+  if (model::cells[model::cell_map[filled_with_triso_base_]]->universe_ ==
+      i_univ) {
+    p.lowest_coord().cell() = model::cell_map[filled_with_triso_base_];
+    return true;
+  }
+  return false;
+}
 
 BoundingBox Universe::bounding_box() const
 {