diff --git a/include/openmc/cell.h b/include/openmc/cell.h
index 6343af605..f436592cd 100644
--- a/include/openmc/cell.h
+++ b/include/openmc/cell.h
@@ -3,6 +3,7 @@
 
 #include <cstdint>
 #include <limits>
+#include <memory> // for unique_ptr
 #include <string>
 #include <unordered_map>
 #include <vector>
@@ -44,15 +45,11 @@ class Cell;
 class Universe;
 
 namespace model {
+  extern std::vector<std::unique_ptr<Cell>> cells;
+  extern std::unordered_map<int32_t, int32_t> cell_map;
 
-extern "C" int32_t n_cells;
-
-extern std::vector<Cell*> cells;
-extern std::unordered_map<int32_t, int32_t> cell_map;
-
-extern std::vector<Universe*> universes;
-extern std::unordered_map<int32_t, int32_t> universe_map;
-
+  extern std::vector<std::unique_ptr<Universe>> universes;
+  extern std::unordered_map<int32_t, int32_t> universe_map;
 } // namespace model
 
 //==============================================================================
diff --git a/include/openmc/lattice.h b/include/openmc/lattice.h
index 68728de5f..ba7cbdf5e 100644
--- a/include/openmc/lattice.h
+++ b/include/openmc/lattice.h
@@ -3,6 +3,7 @@
 
 #include <array>
 #include <cstdint>
+#include <memory> // for unique_ptr
 #include <string>
 #include <unordered_map>
 #include <vector>
@@ -33,10 +34,8 @@ enum class LatticeType {
 class Lattice;
 
 namespace model {
-
-extern std::vector<Lattice*> lattices;
-extern std::unordered_map<int32_t, int32_t> lattice_map;
-
+  extern std::vector<std::unique_ptr<Lattice>> lattices;
+  extern std::unordered_map<int32_t, int32_t> lattice_map;
 } // namespace model
 
 //==============================================================================
diff --git a/include/openmc/material.h b/include/openmc/material.h
index f43db5928..8a09dd9a2 100644
--- a/include/openmc/material.h
+++ b/include/openmc/material.h
@@ -23,7 +23,7 @@ class Material;
 
 namespace model {
 
-extern std::vector<Material*> materials;
+extern std::vector<std::unique_ptr<Material>> materials;
 extern std::unordered_map<int32_t, int32_t> material_map;
 
 } // namespace model
diff --git a/include/openmc/surface.h b/include/openmc/surface.h
index 69944ef92..e1acefa23 100644
--- a/include/openmc/surface.h
+++ b/include/openmc/surface.h
@@ -1,9 +1,10 @@
 #ifndef OPENMC_SURFACE_H
 #define OPENMC_SURFACE_H
 
-#include <map>
+#include <memory>  // for unique_ptr
 #include <limits>  // For numeric_limits
 #include <string>
+#include <unordered_map>
 #include <vector>
 
 #include "hdf5.h"
@@ -35,10 +36,8 @@ extern "C" const int BC_PERIODIC;
 class Surface;
 
 namespace model {
-
-extern std::vector<Surface*> surfaces;
-extern std::map<int, int> surface_map;
-
+  extern std::vector<std::unique_ptr<Surface>> surfaces;
+  extern std::unordered_map<int, int> surface_map;
 } // namespace model
 
 //==============================================================================
diff --git a/src/cell.cpp b/src/cell.cpp
index 2608d8f6f..d4712b9d1 100644
--- a/src/cell.cpp
+++ b/src/cell.cpp
@@ -22,13 +22,11 @@ namespace openmc {
 //==============================================================================
 
 namespace model {
+  std::vector<std::unique_ptr<Cell>> cells;
+  std::unordered_map<int32_t, int32_t> cell_map;
 
-std::vector<Cell*> cells;
-std::unordered_map<int32_t, int32_t> cell_map;
-
-std::vector<Universe*> universes;
-std::unordered_map<int32_t, int32_t> universe_map;
-
+  std::vector<std::unique_ptr<Universe>> universes;
+  std::unordered_map<int32_t, int32_t> universe_map;
 } // namespace model
 
 //==============================================================================
@@ -653,7 +651,7 @@ void read_cells(pugi::xml_node node)
   // Loop over XML cell elements and populate the array.
   model::cells.reserve(n_cells);
   for (pugi::xml_node cell_node : node.children("cell")) {
-    model::cells.push_back(new CSGCell(cell_node));
+    model::cells.push_back(std::make_unique<CSGCell>(cell_node));
   }
 
   // Fill the cell map.
@@ -674,7 +672,7 @@ void read_cells(pugi::xml_node node)
     int32_t uid = model::cells[i]->universe_;
     auto it = model::universe_map.find(uid);
     if (it == model::universe_map.end()) {
-      model::universes.push_back(new Universe());
+      model::universes.push_back(std::make_unique<Universe>());
       model::universes.back()->id_ = uid;
       model::universes.back()->cells_.push_back(i);
       model::universe_map[uid] = model::universes.size() - 1;
@@ -823,7 +821,7 @@ openmc_extend_cells(int32_t n, int32_t* index_start, int32_t* index_end)
   if (index_start) *index_start = model::cells.size();
   if (index_end) *index_end = model::cells.size() + n - 1;
   for (int32_t i = 0; i < n; i++) {
-    model::cells.push_back(new CSGCell());
+    model::cells.push_back(std::make_unique<CSGCell>());
   }
   return 0;
 }
diff --git a/src/dagmc.cpp b/src/dagmc.cpp
index f417c12f5..a9d116dfa 100644
--- a/src/dagmc.cpp
+++ b/src/dagmc.cpp
@@ -149,8 +149,8 @@ void load_dagmc_geometry()
     // Populate the Universe vector and dict
     auto it = model::universe_map.find(dagmc_univ_id);
     if (it == model::universe_map.end()) {
-      model::universes.push_back(new Universe());
-      model::universes.back()-> id_ = dagmc_univ_id;
+      model::universes.push_back(std::make_unique<Universe>());
+      model::universes.back()->id_ = dagmc_univ_id;
       model::universes.back()->cells_.push_back(i);
       model::universe_map[dagmc_univ_id] = model::universes.size() - 1;
     } else {
@@ -255,7 +255,6 @@ void load_dagmc_geometry()
 
   // initialize surface objects
   int n_surfaces = model::DAG->num_entities(2);
-  model::surfaces.resize(n_surfaces);
 
   for (int i = 0; i < n_surfaces; i++) {
     moab::EntityHandle surf_handle = model::DAG->entity_by_index(2, i+1);
@@ -303,7 +302,7 @@ void load_dagmc_geometry()
     }
 
     // add to global array and map
-    model::surfaces[i] = s;
+    model::surfaces.push_back(s);
     model::surface_map[s->id_] = s->id_;
   }
 
diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp
index c118c04ee..3e64aae2a 100644
--- a/src/geometry_aux.cpp
+++ b/src/geometry_aux.cpp
@@ -66,7 +66,7 @@ void
 adjust_indices()
 {
   // Adjust material/fill idices.
-  for (Cell* c : model::cells) {
+  for (auto& c : model::cells) {
     if (c->fill_ != C_NONE) {
       int32_t id = c->fill_;
       auto search_univ = model::universe_map.find(id);
@@ -102,7 +102,7 @@ adjust_indices()
   }
 
   // Change cell.universe values from IDs to indices.
-  for (Cell* c : model::cells) {
+  for (auto& c : model::cells) {
     auto search = model::universe_map.find(c->universe_);
     if (search != model::universe_map.end()) {
       c->universe_ = search->second;
@@ -115,7 +115,7 @@ adjust_indices()
   }
 
   // Change all lattice universe values from IDs to indices.
-  for (Lattice* l : model::lattices) {
+  for (auto& l : model::lattices) {
     l->adjust_indices();
   }
 }
@@ -125,7 +125,7 @@ adjust_indices()
 void
 assign_temperatures()
 {
-  for (Cell* c : model::cells) {
+  for (auto& c : model::cells) {
     // Ignore non-material cells and cells with defined temperature.
     if (c->material_.size() == 0) continue;
     if (c->sqrtkT_.size() > 0) continue;
@@ -224,12 +224,12 @@ find_root_universe()
 {
   // Find all the universes listed as a cell fill.
   std::unordered_set<int32_t> fill_univ_ids;
-  for (Cell* c : model::cells) {
+  for (const auto& c : model::cells) {
     fill_univ_ids.insert(c->fill_);
   }
 
   // Find all the universes contained in a lattice.
-  for (Lattice* lat : model::lattices) {
+  for (const auto& lat : model::lattices) {
     for (auto it = lat->begin(); it != lat->end(); ++it) {
       fill_univ_ids.insert(*it);
     }
@@ -308,7 +308,7 @@ prepare_distribcell()
   // unique distribcell array index.
   int distribcell_index = 0;
   std::vector<int32_t> target_univ_ids;
-  for (Universe* u : model::universes) {
+  for (const auto& u : model::universes) {
     for (auto cell_indx : u->cells_) {
       if (distribcells.find(cell_indx) != distribcells.end()) {
         model::cells[cell_indx]->distribcell_index_ = distribcell_index;
@@ -320,19 +320,19 @@ prepare_distribcell()
 
   // Allocate the cell and lattice offset tables.
   int n_maps = target_univ_ids.size();
-  for (Cell* c : model::cells) {
+  for (auto& c : model::cells) {
     if (c->type_ != FILL_MATERIAL) {
       c->offset_.resize(n_maps, C_NONE);
     }
   }
-  for (Lattice* lat : model::lattices) {
+  for (auto& lat : model::lattices) {
     lat->allocate_offset_table(n_maps);
   }
 
   // Fill the cell and lattice offset tables.
   for (int map = 0; map < target_univ_ids.size(); map++) {
     auto target_univ_id = target_univ_ids[map];
-    for (Universe* univ : model::universes) {
+    for (const auto& univ : model::universes) {
       int32_t offset {0};  // TODO: is this a bug?  It matches F90 implementation.
       for (int32_t cell_indx : univ->cells_) {
         Cell& c = *model::cells[cell_indx];
@@ -515,15 +515,12 @@ maximum_levels(int32_t univ)
 void
 free_memory_geometry()
 {
-  for (Cell* c : model::cells) {delete c;}
   model::cells.clear();
   model::cell_map.clear();
 
-  for (Universe* u : model::universes) {delete u;}
   model::universes.clear();
   model::universe_map.clear();
 
-  for (Lattice* lat : model::lattices) {delete lat;}
   model::lattices.clear();
   model::lattice_map.clear();
 
diff --git a/src/lattice.cpp b/src/lattice.cpp
index d8185de62..780300b28 100644
--- a/src/lattice.cpp
+++ b/src/lattice.cpp
@@ -19,10 +19,8 @@ namespace openmc {
 //==============================================================================
 
 namespace model {
-
-std::vector<Lattice*> lattices;
-std::unordered_map<int32_t, int32_t> lattice_map;
-
+  std::vector<std::unique_ptr<Lattice>> lattices;
+  std::unordered_map<int32_t, int32_t> lattice_map;
 }
 
 //==============================================================================
@@ -867,10 +865,10 @@ HexLattice::to_hdf5_inner(hid_t lat_group) const
 void read_lattices(pugi::xml_node node)
 {
   for (pugi::xml_node lat_node : node.children("lattice")) {
-    model::lattices.push_back(new RectLattice(lat_node));
+    model::lattices.push_back(std::make_unique<RectLattice>(lat_node));
   }
   for (pugi::xml_node lat_node : node.children("hex_lattice")) {
-    model::lattices.push_back(new HexLattice(lat_node));
+    model::lattices.push_back(std::make_unique<HexLattice>(lat_node));
   }
 
   // Fill the lattice map.
diff --git a/src/material.cpp b/src/material.cpp
index 8236a4ef0..2eed3cc71 100644
--- a/src/material.cpp
+++ b/src/material.cpp
@@ -37,7 +37,7 @@ namespace openmc {
 
 namespace model {
 
-std::vector<Material*> materials;
+std::vector<std::unique_ptr<Material>> materials;
 std::unordered_map<int32_t, int32_t> material_map;
 
 } // namespace model
@@ -895,7 +895,7 @@ void read_materials_xml()
   // Loop over XML material elements and populate the array.
   pugi::xml_node root = doc.document_element();
   for (pugi::xml_node material_node : root.children("material")) {
-    model::materials.push_back(new Material(material_node));
+    model::materials.push_back(std::make_unique<Material>(material_node));
   }
   model::materials.shrink_to_fit();
 
@@ -915,7 +915,6 @@ void read_materials_xml()
 
 void free_memory_material()
 {
-  for (Material *mat : model::materials) {delete mat;}
   model::materials.clear();
   model::material_map.clear();
 }
@@ -942,7 +941,7 @@ openmc_material_add_nuclide(int32_t index, const char* name, double density)
 {
   int err = 0;
   if (index >= 0 && index < model::materials.size()) {
-    Material* m = model::materials[index];
+    auto& m = model::materials[index];
 
     // Check if nuclide is already in material
     for (int i = 0; i < m->nuclide_.size(); ++i) {
@@ -1032,7 +1031,7 @@ extern "C" int
 openmc_material_get_volume(int32_t index, double* volume)
 {
   if (index >= 0 && index < model::materials.size()) {
-    Material* m = model::materials[index];
+    auto& m = model::materials[index];
     if (m->volume_ >= 0.0) {
       *volume = m->volume_;
       return 0;
@@ -1112,7 +1111,7 @@ extern "C" int
 openmc_material_set_volume(int32_t index, double volume)
 {
   if (index >= 0 && index < model::materials.size()) {
-    Material* m = model::materials[index];
+    auto& m {model::materials[index]};
     if (volume >= 0.0) {
       m->volume_ = volume;
       return 0;
@@ -1132,7 +1131,7 @@ openmc_extend_materials(int32_t n, int32_t* index_start, int32_t* index_end)
   if (index_start) *index_start = model::materials.size();
   if (index_end) *index_end = model::materials.size() + n - 1;
   for (int32_t i = 0; i < n; i++) {
-    model::materials.push_back(new Material());
+    model::materials.push_back(std::make_unique<Material>());
   }
   return 0;
 }
diff --git a/src/mesh.cpp b/src/mesh.cpp
index ab263cd77..068d02918 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -728,7 +728,7 @@ openmc_extend_meshes(int32_t n, int32_t* index_start, int32_t* index_end)
 {
   if (index_start) *index_start = model::meshes.size();
   for (int i = 0; i < n; ++i) {
-    model::meshes.emplace_back(new RegularMesh{});
+    model::meshes.push_back(std::make_unique<RegularMesh>());
   }
   if (index_end) *index_end = model::meshes.size() - 1;
 
@@ -866,7 +866,7 @@ void read_meshes(pugi::xml_node root)
 {
   for (auto node : root.children("mesh")) {
     // Read mesh and add to vector
-    model::meshes.emplace_back(new RegularMesh{node});
+    model::meshes.push_back(std::make_unique<RegularMesh>(node));
 
     // Map ID to position in vector
     model::mesh_map[model::meshes.back()->id_] = model::meshes.size() - 1;
diff --git a/src/mgxs_interface.cpp b/src/mgxs_interface.cpp
index d2913a17a..fb56b702f 100644
--- a/src/mgxs_interface.cpp
+++ b/src/mgxs_interface.cpp
@@ -136,7 +136,7 @@ void create_macro_xs()
   for (int i = 0; i < model::materials.size(); ++i) {
     if (kTs[i].size() > 0) {
       // Convert atom_densities to a vector
-      Material* mat = model::materials[i];
+      auto& mat {model::materials[i]};
       std::vector<double> atom_densities(mat->atom_density_.begin(),
         mat->atom_density_.end());
 
diff --git a/src/particle.cpp b/src/particle.cpp
index ccce1502f..dfb9fca61 100644
--- a/src/particle.cpp
+++ b/src/particle.cpp
@@ -402,7 +402,7 @@ Particle::cross_surface()
 {
   int i_surface = std::abs(surface);
   // TODO: off-by-one
-  const auto& surf {model::surfaces[i_surface - 1]};
+  const auto& surf {model::surfaces[i_surface - 1].get()};
   if (settings::verbosity >= 10 || simulation::trace) {
     write_message("    Crossing surface " + std::to_string(surf->id_));
   }
@@ -532,7 +532,7 @@ Particle::cross_surface()
     // Get a pointer to the partner periodic surface
     auto surf_p = dynamic_cast<PeriodicSurface*>(surf);
     auto other = dynamic_cast<PeriodicSurface*>(
-      model::surfaces[surf_p->i_periodic_]);
+      model::surfaces[surf_p->i_periodic_].get());
 
     // Adjust the particle's location and direction.
     Position r {coord[0].xyz};
diff --git a/src/plot.cpp b/src/plot.cpp
index ece9c76e7..e5f5092ea 100644
--- a/src/plot.cpp
+++ b/src/plot.cpp
@@ -666,7 +666,7 @@ void position_rgb(Particle p, Plot pl, RGBColor& rgb, int& id)
   } else {
     if (PlotColorBy::mats == pl.color_by_) {
       // Assign color based on material
-      Cell* c = model::cells[p.coord[j].cell];
+      const auto& c = model::cells[p.coord[j].cell];
       if (c->type_ == FILL_UNIVERSE) {
         // If we stopped on a middle universe level, treat as if not found
         rgb = pl.not_found_;
diff --git a/src/reaction_product.cpp b/src/reaction_product.cpp
index 7e3175343..585cff912 100644
--- a/src/reaction_product.cpp
+++ b/src/reaction_product.cpp
@@ -69,13 +69,13 @@ ReactionProduct::ReactionProduct(hid_t group)
     // Determine distribution type and read data
     read_attribute(dgroup, "type", temp);
     if (temp == "uncorrelated") {
-      distribution_.emplace_back(new UncorrelatedAngleEnergy{dgroup});
+      distribution_.push_back(std::make_unique<UncorrelatedAngleEnergy>(dgroup));
     } else if (temp == "correlated") {
-      distribution_.emplace_back(new CorrelatedAngleEnergy{dgroup});
+      distribution_.push_back(std::make_unique<CorrelatedAngleEnergy>(dgroup));
     } else if (temp == "nbody") {
-      distribution_.emplace_back(new NBodyPhaseSpace{dgroup});
+      distribution_.push_back(std::make_unique<NBodyPhaseSpace>(dgroup));
     } else if (temp == "kalbach-mann") {
-      distribution_.emplace_back(new KalbachMann{dgroup});
+      distribution_.push_back(std::make_unique<KalbachMann>(dgroup));
     }
 
     close_group(dgroup);
diff --git a/src/settings.cpp b/src/settings.cpp
index 92e85ab17..dfa5ebb7e 100644
--- a/src/settings.cpp
+++ b/src/settings.cpp
@@ -509,7 +509,7 @@ void read_settings_xml()
 
     // Read entropy mesh from <entropy>
     auto node_entropy = root.child("entropy");
-    model::meshes.emplace_back(new RegularMesh{node_entropy});
+    model::meshes.push_back(std::make_unique<RegularMesh>(node_entropy));
 
     // Set entropy mesh index
     index_entropy_mesh = model::meshes.size() - 1;
@@ -555,7 +555,7 @@ void read_settings_xml()
 
     // Read entropy mesh from <entropy>
     auto node_ufs = root.child("uniform_fs");
-    model::meshes.emplace_back(new RegularMesh{node_ufs});
+    model::meshes.push_back(std::make_unique<RegularMesh>(node_ufs));
 
     // Set entropy mesh index
     index_ufs_mesh = model::meshes.size() - 1;
diff --git a/src/source.cpp b/src/source.cpp
index 1ebdc235b..720bd07ea 100644
--- a/src/source.cpp
+++ b/src/source.cpp
@@ -172,7 +172,7 @@ Bank SourceDistribution::sample() const
       if (space_box) {
         if (space_box->only_fissionable()) {
           // Determine material
-          auto c = model::cells[cell_index];
+          const auto& c = model::cells[cell_index];
           auto mat_index = c->material_.size() == 1
             ? c->material_[0] : c->material_[instance];
 
diff --git a/src/summary.cpp b/src/summary.cpp
index 1eeab5d04..500a5c15c 100644
--- a/src/summary.cpp
+++ b/src/summary.cpp
@@ -102,19 +102,19 @@ void write_geometry(hid_t file)
   write_attribute(geom_group, "n_lattices", model::lattices.size());
 
   auto cells_group = create_group(geom_group, "cells");
-  for (Cell* c : model::cells) c->to_hdf5(cells_group);
+  for (const auto& c : model::cells) c->to_hdf5(cells_group);
   close_group(cells_group);
 
   auto surfaces_group = create_group(geom_group, "surfaces");
-  for (Surface* surf : model::surfaces) surf->to_hdf5(surfaces_group);
+  for (const auto& surf : model::surfaces) surf->to_hdf5(surfaces_group);
   close_group(surfaces_group);
 
   auto universes_group = create_group(geom_group, "universes");
-  for (Universe* u : model::universes) u->to_hdf5(universes_group);
+  for (const auto& u : model::universes) u->to_hdf5(universes_group);
   close_group(universes_group);
 
   auto lattices_group = create_group(geom_group, "lattices");
-  for (Lattice* lat : model::lattices) lat->to_hdf5(lattices_group);
+  for (const auto& lat : model::lattices) lat->to_hdf5(lattices_group);
   close_group(lattices_group);
 
   close_group(geom_group);
diff --git a/src/surface.cpp b/src/surface.cpp
index d56adbfd9..bf0315698 100644
--- a/src/surface.cpp
+++ b/src/surface.cpp
@@ -27,10 +27,8 @@ extern "C" const int BC_PERIODIC {3};
 //==============================================================================
 
 namespace model {
-
-std::vector<Surface*> surfaces;
-std::map<int, int> surface_map;
-
+  std::vector<std::unique_ptr<Surface>> surfaces;
+  std::unordered_map<int, int> surface_map;
 } // namespace model
 
 //==============================================================================
@@ -1080,40 +1078,40 @@ void read_surfaces(pugi::xml_node node)
       std::string surf_type = get_node_value(surf_node, "type", true, true);
 
       if (surf_type == "x-plane") {
-        model::surfaces.push_back(new SurfaceXPlane(surf_node));
+        model::surfaces.push_back(std::make_unique<SurfaceXPlane>(surf_node));
 
       } else if (surf_type == "y-plane") {
-        model::surfaces.push_back(new SurfaceYPlane(surf_node));
+        model::surfaces.push_back(std::make_unique<SurfaceYPlane>(surf_node));
 
       } else if (surf_type == "z-plane") {
-        model::surfaces.push_back(new SurfaceZPlane(surf_node));
+        model::surfaces.push_back(std::make_unique<SurfaceZPlane>(surf_node));
 
       } else if (surf_type == "plane") {
-        model::surfaces.push_back(new SurfacePlane(surf_node));
+        model::surfaces.push_back(std::make_unique<SurfacePlane>(surf_node));
 
       } else if (surf_type == "x-cylinder") {
-        model::surfaces.push_back(new SurfaceXCylinder(surf_node));
+        model::surfaces.push_back(std::make_unique<SurfaceXCylinder>(surf_node));
 
       } else if (surf_type == "y-cylinder") {
-        model::surfaces.push_back(new SurfaceYCylinder(surf_node));
+        model::surfaces.push_back(std::make_unique<SurfaceYCylinder>(surf_node));
 
       } else if (surf_type == "z-cylinder") {
-        model::surfaces.push_back(new SurfaceZCylinder(surf_node));
+        model::surfaces.push_back(std::make_unique<SurfaceZCylinder>(surf_node));
 
       } else if (surf_type == "sphere") {
-        model::surfaces.push_back(new SurfaceSphere(surf_node));
+        model::surfaces.push_back(std::make_unique<SurfaceSphere>(surf_node));
 
       } else if (surf_type == "x-cone") {
-        model::surfaces.push_back(new SurfaceXCone(surf_node));
+        model::surfaces.push_back(std::make_unique<SurfaceXCone>(surf_node));
 
       } else if (surf_type == "y-cone") {
-        model::surfaces.push_back(new SurfaceYCone(surf_node));
+        model::surfaces.push_back(std::make_unique<SurfaceYCone>(surf_node));
 
       } else if (surf_type == "z-cone") {
-        model::surfaces.push_back(new SurfaceZCone(surf_node));
+        model::surfaces.push_back(std::make_unique<SurfaceZCone>(surf_node));
 
       } else if (surf_type == "quadric") {
-        model::surfaces.push_back(new SurfaceQuadric(surf_node));
+        model::surfaces.push_back(std::make_unique<SurfaceQuadric>(surf_node));
 
       } else {
         std::stringstream err_msg;
@@ -1143,8 +1141,8 @@ void read_surfaces(pugi::xml_node node)
   for (int i_surf = 0; i_surf < model::surfaces.size(); i_surf++) {
     if (model::surfaces[i_surf]->bc_ == BC_PERIODIC) {
       // Downcast to the PeriodicSurface type.
-      Surface* surf_base = model::surfaces[i_surf];
-      PeriodicSurface* surf = dynamic_cast<PeriodicSurface*>(surf_base);
+      Surface* surf_base = model::surfaces[i_surf].get();
+      auto surf = dynamic_cast<PeriodicSurface*>(surf_base);
 
       // Make sure this surface inherits from PeriodicSurface.
       if (!surf) {
@@ -1188,8 +1186,8 @@ void read_surfaces(pugi::xml_node node)
   for (int i_surf = 0; i_surf < model::surfaces.size(); i_surf++) {
     if (model::surfaces[i_surf]->bc_ == BC_PERIODIC) {
       // Downcast to the PeriodicSurface type.
-      Surface* surf_base = model::surfaces[i_surf];
-      PeriodicSurface* surf = dynamic_cast<PeriodicSurface*>(surf_base);
+      Surface* surf_base = model::surfaces[i_surf].get();
+      auto surf = dynamic_cast<PeriodicSurface*>(surf_base);
 
       // Also try downcasting to the SurfacePlane type (which must be handled
       // differently).
@@ -1260,7 +1258,6 @@ void read_surfaces(pugi::xml_node node)
 
 void free_memory_surfaces()
 {
-  for (Surface* surf : model::surfaces) {delete surf;}
   model::surfaces.clear();
   model::surface_map.clear();
 }