diff --git a/src/cell.cpp b/src/cell.cpp
index 9e640ce4e..d09cde948 100644
--- a/src/cell.cpp
+++ b/src/cell.cpp
@@ -9,7 +9,7 @@
 #include <set>
 #include <string>
 
-
+#include <fmt/format.h>
 #include <gsl/gsl>
 #include <fmt/format.h>
 
@@ -158,10 +158,8 @@ generate_rpn(int32_t cell_id, std::vector<int32_t> infix)
         // If we run out of operators without finding a left parenthesis, it
         // means there are mismatched parentheses.
         if (it == stack.rend()) {
-          std::stringstream err_msg;
-          err_msg << "Mismatched parentheses in region specification for cell "
-                  << cell_id;
-          fatal_error(err_msg);
+          fatal_error(fmt::format(
+            "Mismatched parentheses in region specification for cell {}", cell_id));
         }
         rpn.push_back(stack.back());
         stack.pop_back();
@@ -177,10 +175,8 @@ generate_rpn(int32_t cell_id, std::vector<int32_t> infix)
 
     // If the operator is a parenthesis it is mismatched.
     if (op >= OP_RIGHT_PAREN) {
-      std::stringstream err_msg;
-      err_msg << "Mismatched parentheses in region specification for cell "
-              << cell_id;
-      fatal_error(err_msg);
+      fatal_error(fmt::format(
+        "Mismatched parentheses in region specification for cell {}", cell_id));
     }
 
     rpn.push_back(stack.back());
@@ -198,9 +194,7 @@ void
 Universe::to_hdf5(hid_t universes_group) const
 {
   // Create a group for this universe.
-  std::stringstream group_name;
-  group_name << "universe " << id_;
-  auto group = create_group(universes_group, group_name);
+  auto group = create_group(universes_group, fmt::format("universe {}", id_));
 
   // Write the contained cells.
   if (cells_.size() > 0) {
@@ -301,22 +295,19 @@ CSGCell::CSGCell(pugi::xml_node cell_node)
   bool fill_present = check_for_node(cell_node, "fill");
   bool material_present = check_for_node(cell_node, "material");
   if (!(fill_present || material_present)) {
-    std::stringstream err_msg;
-    err_msg << "Neither material nor fill was specified for cell " << id_;
-    fatal_error(err_msg);
+    fatal_error(fmt::format(
+      "Neither material nor fill was specified for cell {}", id_));
   }
   if (fill_present && material_present) {
-    std::stringstream err_msg;
-    err_msg << "Cell " << id_ << " has both a material and a fill specified; "
-            << "only one can be specified per cell";
-    fatal_error(err_msg);
+    fatal_error(fmt::format("Cell {} has both a material and a fill specified; "
+      "only one can be specified per cell", id_));
   }
 
   if (fill_present) {
     fill_ = std::stoi(get_node_value(cell_node, "fill"));
     if (fill_ == universe_) {
-      fatal_error("Cell " + std::to_string(id_) +
-        " is filled with the same universe that it is contained in.");
+      fatal_error(fmt::format("Cell {} is filled with the same universe that"
+        "it is contained in.", id_));
     }
   } else {
     fill_ = C_NONE;
@@ -338,9 +329,8 @@ CSGCell::CSGCell(pugi::xml_node cell_node)
         }
       }
     } else {
-      std::stringstream err_msg;
-      err_msg << "An empty material element was specified for cell " << id_;
-      fatal_error(err_msg);
+      fatal_error(fmt::format("An empty material element was specified for cell {}",
+        id_));
     }
   }
 
@@ -351,20 +341,16 @@ CSGCell::CSGCell(pugi::xml_node cell_node)
 
     // Make sure this is a material-filled cell.
     if (material_.size() == 0) {
-      std::stringstream err_msg;
-      err_msg << "Cell " << id_ << " was specified with a temperature but "
-           "no material. Temperature specification is only valid for cells "
-           "filled with a material.";
-      fatal_error(err_msg);
+      fatal_error(fmt::format(
+        "Cell {} was specified with a temperature but no material. Temperature"
+        "specification is only valid for cells filled with a material.", id_));
     }
 
     // Make sure all temperatures are non-negative.
     for (auto T : sqrtkT_) {
       if (T < 0) {
-        std::stringstream err_msg;
-        err_msg << "Cell " << id_
-                << " was specified with a negative temperature";
-        fatal_error(err_msg);
+        fatal_error(fmt::format(
+          "Cell {} was specified with a negative temperature", id_));
       }
     }
 
@@ -426,17 +412,14 @@ CSGCell::CSGCell(pugi::xml_node cell_node)
   // Read the translation vector.
   if (check_for_node(cell_node, "translation")) {
     if (fill_ == C_NONE) {
-      std::stringstream err_msg;
-      err_msg << "Cannot apply a translation to cell " << id_
-              << " because it is not filled with another universe";
-      fatal_error(err_msg);
+      fatal_error(fmt::format("Cannot apply a translation to cell {}"
+        " because it is not filled with another universe", id_));
     }
 
     auto xyz {get_node_array<double>(cell_node, "translation")};
     if (xyz.size() != 3) {
-      std::stringstream err_msg;
-      err_msg << "Non-3D translation vector applied to cell " << id_;
-      fatal_error(err_msg);
+      fatal_error(fmt::format(
+        "Non-3D translation vector applied to cell {}", id_));
     }
     translation_ = xyz;
   }
@@ -444,17 +427,14 @@ CSGCell::CSGCell(pugi::xml_node cell_node)
   // Read the rotation transform.
   if (check_for_node(cell_node, "rotation")) {
     if (fill_ == C_NONE) {
-      std::stringstream err_msg;
-      err_msg << "Cannot apply a rotation to cell " << id_
-              << " because it is not filled with another universe";
-      fatal_error(err_msg);
+      fatal_error(fmt::format("Cannot apply a rotation to cell {}"
+        " because it is not filled with another universe", id_));
     }
 
     auto rot {get_node_array<double>(cell_node, "rotation")};
     if (rot.size() != 3 && rot.size() != 9) {
-      std::stringstream err_msg;
-      err_msg << "Non-3D rotation vector applied to cell " << id_;
-      fatal_error(err_msg);
+      fatal_error(fmt::format(
+        "Non-3D rotation vector applied to cell {}", id_));
     }
 
     // Compute and store the rotation matrix.
@@ -534,9 +514,7 @@ void
 CSGCell::to_hdf5(hid_t cell_group) const
 {
   // Create a group for this cell.
-  std::stringstream group_name;
-  group_name << "cell " << id_;
-  auto group = create_group(cell_group, group_name);
+  auto group = create_group(cell_group, fmt::format("cell {}", id_));
 
   if (!name_.empty()) {
     write_string(group, "name", name_, false);
@@ -1013,9 +991,7 @@ void read_cells(pugi::xml_node node)
     if (search == model::cell_map.end()) {
       model::cell_map[id] = i;
     } else {
-      std::stringstream err_msg;
-      err_msg << "Two or more cells use the same unique ID: " << id;
-      fatal_error(err_msg);
+      fatal_error(fmt::format("Two or more cells use the same unique ID: {}", id));
     }
   }
 
diff --git a/src/dagmc.cpp b/src/dagmc.cpp
index e14a6062a..3f16799a5 100644
--- a/src/dagmc.cpp
+++ b/src/dagmc.cpp
@@ -12,11 +12,10 @@
 #include "openmc/surface.h"
 
 #ifdef DAGMC
-
 #include "uwuw.hpp"
 #include "dagmcmetadata.hpp"
-
 #endif
+#include <fmt/format.h>
 
 #include <string>
 #include <sstream>
@@ -110,11 +109,10 @@ void legacy_assign_material(const std::string& mat_string, DAGCell* c)
         c->material_.push_back(m->id_);
       // report error if more than one material is found
       } else {
-        std::stringstream err_msg;
-        err_msg << "More than one material found with name " << mat_string
-                << ". Please ensure materials have unique names if using this"
-                << " property to assign materials.";
-        fatal_error(err_msg);
+        fatal_error(fmt::format(
+          "More than one material found with name {}. Please ensure materials "
+          "have unique names if using this property to assign materials.",
+          mat_string));
       }
     }
   }
@@ -125,10 +123,8 @@ void legacy_assign_material(const std::string& mat_string, DAGCell* c)
       auto id = std::stoi(mat_string);
       c->material_.emplace_back(id);
     } catch (const std::invalid_argument&) {
-      std::stringstream err_msg;
-      err_msg << "No material " << mat_string
-              << " found for volume (cell) " << c->id_;
-      fatal_error(err_msg);
+      fatal_error(fmt::format(
+        "No material {} found for volume (cell) {}", mat_string, c->id_));
     }
   }
 
@@ -153,7 +149,6 @@ void load_dagmc_geometry()
     model::DAG = new moab::DagMC();
   }
 
-
   std::string filename = settings::path_input + DAGMC_FILENAME;
   // --- Materials ---
 
@@ -253,9 +248,7 @@ void load_dagmc_geometry()
       rval = model::DAG->prop_value(vol_handle, "mat", mat_value);
       MB_CHK_ERR_CONT(rval);
     } else {
-      std::stringstream err_msg;
-      err_msg << "Volume " << c->id_ << " has no material assignment.";
-      fatal_error(err_msg.str());
+      fatal_error(fmt::format("Volume {} has no material assignment.", c->id_));
     }
 
     std::string cmp_str = mat_value;
@@ -277,10 +270,8 @@ void load_dagmc_geometry()
           int mat_number = uwuw.material_library[uwuw_mat].metadata["mat_number"].asInt();
           c->material_.push_back(mat_number);
         } else {
-          std::stringstream err_msg;
-          err_msg << "Material with value " << mat_value << " not found ";
-          err_msg << "in the UWUW material library";
-          fatal_error(err_msg);
+          fatal_error(fmt::format("Material with value {} not found in the "
+            "UWUW material library", mat_value));
         }
       } else {
         legacy_assign_material(mat_value, c);
@@ -348,10 +339,8 @@ void load_dagmc_geometry()
       } else if (bc_value == "periodic") {
         fatal_error("Periodic boundary condition not supported in DAGMC.");
       } else {
-        std::stringstream err_msg;
-        err_msg << "Unknown boundary condition \"" << bc_value
-                << "\" specified on surface " << s->id_;
-        fatal_error(err_msg);
+        fatal_error(fmt::format("Unknown boundary condition \"{}\" specified "
+          "on surface {}", bc_value, s->id_));
       }
     } else {
       // if no condition is found, set to transmit
diff --git a/src/distribution_spatial.cpp b/src/distribution_spatial.cpp
index 12432baf7..125d8afa6 100644
--- a/src/distribution_spatial.cpp
+++ b/src/distribution_spatial.cpp
@@ -96,9 +96,7 @@ CylindricalIndependent::CylindricalIndependent(pugi::xml_node node)
     if (origin.size() == 3) {
       origin_ = origin;
     } else {
-      std::stringstream err_msg;
-      err_msg << "Origin for cylindrical source distribution must be length 3";
-      fatal_error(err_msg);
+      fatal_error("Origin for cylindrical source distribution must be length 3");
     }
   } else {
     // If no coordinates were specified, default to (0, 0, 0)
@@ -162,9 +160,7 @@ SphericalIndependent::SphericalIndependent(pugi::xml_node node)
     if (origin.size() == 3) {
       origin_ = origin;
     } else {
-      std::stringstream err_msg;
-      err_msg << "Origin for spherical source distribution must be length 3";
-      fatal_error(err_msg);
+      fatal_error("Origin for spherical source distribution must be length 3");
     }
   } else {
     // If no coordinates were specified, default to (0, 0, 0)
diff --git a/src/geometry.cpp b/src/geometry.cpp
index 8bbc4d90c..d2012d901 100644
--- a/src/geometry.cpp
+++ b/src/geometry.cpp
@@ -1,7 +1,9 @@
 #include "openmc/geometry.h"
 
 #include <array>
-#include <sstream>
+
+#include <fmt/format.h>
+#include <fmt/ostream.h>
 
 #include "openmc/cell.h"
 #include "openmc/constants.h"
@@ -47,11 +49,9 @@ bool check_cell_overlap(Particle* p, bool error)
       if (c.contains(p->coord_[j].r, p->coord_[j].u, p->surface_)) {
         if (index_cell != p->coord_[j].cell) {
           if (error) {
-            std::stringstream err_msg;
-            err_msg << "Overlapping cells detected: " << c.id_ << ", "
-                    << model::cells[p->coord_[j].cell]->id_ << " on universe "
-                    << univ.id_;
-            fatal_error(err_msg);
+            fatal_error(fmt::format(
+              "Overlapping cells detected: {}, {} on universe {}",
+              c.id_, model::cells[p->coord_[j].cell]->id_, univ.id_));
           }
           return true;
         }
@@ -120,8 +120,7 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list)
 
   // Announce the cell that the particle is entering.
   if (found && (settings::verbosity >= 10 || p->trace_)) {
-    std::stringstream msg;
-    msg << "    Entering cell " << model::cells[i_cell]->id_;
+    auto msg = fmt::format("    Entering cell {}", model::cells[i_cell]->id_);
     write_message(msg, 1);
   }
 
@@ -229,11 +228,8 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list)
         if (lat.outer_ != NO_OUTER_UNIVERSE) {
           coord.universe = lat.outer_;
         } else {
-          std::stringstream err_msg;
-          err_msg << "Particle " << p->id_ << " is outside lattice "
-                  << lat.id_ << " but the lattice has no defined outer "
-                  "universe.";
-          warning(err_msg);
+          warning(fmt::format("Particle {} is outside lattice {} but the "
+            "lattice has no defined outer universe.", p->id_, lat.id_));
           return false;
         }
       }
@@ -298,11 +294,9 @@ cross_lattice(Particle* p, const BoundaryInfo& boundary)
   auto& lat {*model::lattices[coord.lattice]};
 
   if (settings::verbosity >= 10 || p->trace_) {
-    std::stringstream msg;
-    msg << "    Crossing lattice " << lat.id_ << ". Current position ("
-        << coord.lattice_x << "," << coord.lattice_y << ","
-        << coord.lattice_z << "). r=" << p->r();
-    write_message(msg, 1);
+    write_message(fmt::format(
+      "    Crossing lattice {}. Current position ({},{},{}). r={}",
+      lat.id_, coord.lattice_x, coord.lattice_y, coord.lattice_z, p->r()), 1);
   }
 
   // Set the lattice indices.
@@ -326,10 +320,8 @@ cross_lattice(Particle* p, const BoundaryInfo& boundary)
     p->n_coord_ = 1;
     bool found = find_cell(p, 0);
     if (!found && p->alive_) {
-      std::stringstream err_msg;
-      err_msg << "Could not locate particle " << p->id_
-              << " after crossing a lattice boundary";
-      p->mark_as_lost(err_msg);
+      p->mark_as_lost(fmt::format("Could not locate particle {} after "
+        "crossing a lattice boundary", p->id_));
     }
 
   } else {
@@ -343,10 +335,8 @@ cross_lattice(Particle* p, const BoundaryInfo& boundary)
       p->n_coord_ = 1;
       bool found = find_cell(p, 0);
       if (!found && p->alive_) {
-        std::stringstream err_msg;
-        err_msg << "Could not locate particle " << p->id_
-                << " after crossing a lattice boundary";
-        p->mark_as_lost(err_msg);
+        p->mark_as_lost(fmt::format("Could not locate particle {} after "
+          "crossing a lattice boundary", p->id_));
       }
     }
   }
@@ -400,10 +390,8 @@ BoundaryInfo distance_to_boundary(Particle* p)
       level_lat_trans = lattice_distance.second;
 
       if (d_lat < 0) {
-        std::stringstream err_msg;
-        err_msg << "Particle " << p->id_
-                << " had a negative distance to a lattice boundary";
-        p->mark_as_lost(err_msg);
+        p->mark_as_lost(fmt::format(
+          "Particle {} had a negative distance to a lattice boundary", p->id_));
       }
     }
 
@@ -462,10 +450,7 @@ openmc_find_cell(const double* xyz, int32_t* index, int32_t* instance)
   p.u() = {0.0, 0.0, 1.0};
 
   if (!find_cell(&p, false)) {
-    std::stringstream msg;
-    msg << "Could not find cell at position (" << p.r().x << ", " << p.r().y
-      << ", " << p.r().z << ").";
-    set_errmsg(msg);
+    set_errmsg(fmt::format("Could not find cell at position {}.", p.r()));
     return OPENMC_E_GEOMETRY;
   }
 
diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp
index 6ee289949..0df1198b1 100644
--- a/src/geometry_aux.cpp
+++ b/src/geometry_aux.cpp
@@ -4,7 +4,8 @@
 #include <sstream>
 #include <unordered_set>
 
-#include "pugixml.hpp"
+#include <fmt/format.h>
+#include <pugixml.hpp>
 
 #include "openmc/cell.h"
 #include "openmc/constants.h"
@@ -94,10 +95,8 @@ adjust_indices()
         c->type_ = Fill::LATTICE;
         c->fill_ = search_lat->second;
       } else {
-        std::stringstream err_msg;
-        err_msg << "Specified fill " << id << " on cell " << c->id_
-                << " is neither a universe nor a lattice.";
-        fatal_error(err_msg);
+        fatal_error(fmt::format("Specified fill {} on cell {} is neither a "
+          "universe nor a lattice.", id, c->id_));
       }
     } else {
       c->type_ = Fill::MATERIAL;
@@ -105,10 +104,9 @@ adjust_indices()
         if (mat_id != MATERIAL_VOID) {
           auto search = model::material_map.find(mat_id);
           if (search == model::material_map.end()) {
-            std::stringstream err_msg;
-            err_msg << "Could not find material " << mat_id
-                    << " specified on cell " << c->id_;
-            fatal_error(err_msg);
+            fatal_error(fmt::format(
+              "Could not find material {} specified on cell {}",
+              mat_id, c->id_));
           }
           // Change from ID to index
           mat_id = search->second;
@@ -123,10 +121,8 @@ adjust_indices()
     if (search != model::universe_map.end()) {
       c->universe_ = search->second;
     } else {
-      std::stringstream err_msg;
-      err_msg << "Could not find universe " << c->universe_
-              << " specified on cell " << c->id_;
-      fatal_error(err_msg);
+      fatal_error(fmt::format("Could not find universe {} specified on cell {}",
+        c->universe_, c->id_));
     }
   }
 
@@ -345,23 +341,21 @@ prepare_distribcell()
 
     if (c.material_.size() > 1) {
       if (c.material_.size() != c.n_instances_) {
-        std::stringstream err_msg;
-        err_msg <<  "Cell " << c.id_ <<  " was specified with "
-              << c.material_.size() << " materials but has " << c.n_instances_
-              << " distributed instances. The number of materials must equal "
-              "one or the number of instances.";
-        fatal_error(err_msg);
+        fatal_error(fmt::format(
+          "Cell {} was specified with {} materials but has {} distributed "
+          "instances. The number of materials must equal one or the number "
+          "of instances.", c.id_, c.material_.size(), c.n_instances_
+        ));
       }
     }
 
     if (c.sqrtkT_.size() > 1) {
       if (c.sqrtkT_.size() != c.n_instances_) {
-        std::stringstream err_msg;
-        err_msg <<  "Cell " << c.id_ <<  " was specified with "
-          << c.sqrtkT_.size() << " temperatures but has " << c.n_instances_
-          << " distributed instances. The number of temperatures must equal "
-          "one or the number of instances.";
-        fatal_error(err_msg);
+        fatal_error(fmt::format(
+          "Cell {} was specified with {} temperatures but has {} distributed "
+          "instances. The number of temperatures must equal one or the number "
+          "of instances.", c.id_, c.sqrtkT_.size(), c.n_instances_
+        ));
       }
     }
   }
diff --git a/src/hdf5_interface.cpp b/src/hdf5_interface.cpp
index bd73e05fa..f207f53b7 100644
--- a/src/hdf5_interface.cpp
+++ b/src/hdf5_interface.cpp
@@ -2,12 +2,12 @@
 
 #include <array>
 #include <cstring>
-#include <sstream>
 #include <stdexcept>
 #include <string>
 
 #include "xtensor/xtensor.hpp"
 #include "xtensor/xarray.hpp"
+#include <fmt/format.h>
 
 #include "hdf5.h"
 #include "hdf5_hl.h"
@@ -105,9 +105,7 @@ create_group(hid_t parent_id, char const *name)
 {
   hid_t out = H5Gcreate(parent_id, name, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
   if (out < 0) {
-    std::stringstream err_msg;
-    err_msg << "Failed to create HDF5 group \"" << name << "\"";
-    fatal_error(err_msg);
+    fatal_error(fmt::format("Failed to create HDF5 group \"{}\"", name));
   }
   return out;
 }
@@ -161,17 +159,13 @@ ensure_exists(hid_t obj_id, const char* name, bool attribute)
 {
   if (attribute) {
     if (!attribute_exists(obj_id, name)) {
-      std::stringstream err_msg;
-      err_msg << "Attribute \"" << name << "\" does not exist in object "
-              << object_name(obj_id);
-      fatal_error(err_msg);
+      fatal_error(fmt::format("Attribute \"{}\" does not exist in object {}",
+        name, object_name(obj_id)));
     }
   } else {
     if (!object_exists(obj_id, name)) {
-      std::stringstream err_msg;
-      err_msg << "Object \"" << name << "\" does not exist in object "
-              << object_name(obj_id);
-      fatal_error(err_msg);
+      fatal_error(fmt::format("Object \"{}\" does not exist in object {}",
+        name, object_name(obj_id)));
     }
   }
 }
@@ -194,9 +188,7 @@ file_open(const char* filename, char mode, bool parallel)
       flags = (mode == 'x' ? H5F_ACC_EXCL : H5F_ACC_TRUNC);
       break;
     default:
-      std::stringstream err_msg;
-      err_msg <<  "Invalid file mode: " << mode;
-      fatal_error(err_msg);
+      fatal_error(fmt::format("Invalid file mode: ", mode));
   }
 
   hid_t plist = H5P_DEFAULT;
@@ -216,9 +208,8 @@ file_open(const char* filename, char mode, bool parallel)
     file_id = H5Fopen(filename, flags, plist);
   }
   if (file_id < 0) {
-    std::stringstream msg;
-    msg << "Failed to open HDF5 file with mode '" << mode << "': " << filename;
-    fatal_error(msg);
+    fatal_error(fmt::format(
+      "Failed to open HDF5 file with mode '{}': {}", mode, filename));
   }
 
 #ifdef PHDF5
@@ -394,9 +385,7 @@ object_exists(hid_t object_id, const char* name)
 {
   htri_t out = H5LTpath_valid(object_id, name, true);
   if (out < 0) {
-    std::stringstream err_msg;
-    err_msg << "Failed to check if object \"" << name << "\" exists.";
-    fatal_error(err_msg);
+    fatal_error(fmt::format("Failed to check if object \"{}\" exists.", name));
   }
   return (out > 0);
 }
diff --git a/src/initialize.cpp b/src/initialize.cpp
index 22479ae29..a5b516f68 100644
--- a/src/initialize.cpp
+++ b/src/initialize.cpp
@@ -3,13 +3,13 @@
 #include <cstddef>
 #include <cstdlib> // for getenv
 #include <cstring>
-#include <sstream>
 #include <string>
 #include <vector>
 
 #ifdef _OPENMP
 #include <omp.h>
 #endif
+#include <fmt/format.h>
 
 #include "openmc/capi.h"
 #include "openmc/constants.h"
@@ -155,9 +155,8 @@ parse_command_line(int argc, char* argv[])
           settings::path_particle_restart = argv[i];
           settings::particle_restart_run = true;
         } else {
-          std::stringstream msg;
-          msg << "Unrecognized file after restart flag: " << filetype << ".";
-          strcpy(openmc_err_msg, msg.str().c_str());
+          auto msg = fmt::format("Unrecognized file after restart flag: {}.", filetype);
+          strcpy(openmc_err_msg, msg.c_str());
           return OPENMC_E_INVALID_ARGUMENT;
         }
 
@@ -224,7 +223,7 @@ parse_command_line(int argc, char* argv[])
         settings::write_all_tracks = true;
 
       } else {
-        std::cerr << "Unknown option: " << argv[i] << '\n';
+        fmt::print(stderr, "Unknown option: {}\n", argv[i]);
         print_usage();
         return OPENMC_E_UNASSIGNED;
       }
diff --git a/src/lattice.cpp b/src/lattice.cpp
index 092ff3d11..a8dd23aa8 100644
--- a/src/lattice.cpp
+++ b/src/lattice.cpp
@@ -1,9 +1,11 @@
 #include "openmc/lattice.h"
 
 #include <cmath>
-#include <sstream>
+#include <string>
 #include <vector>
 
+#include <fmt/format.h>
+
 #include "openmc/cell.h"
 #include "openmc/error.h"
 #include "openmc/geometry.h"
@@ -71,10 +73,8 @@ Lattice::adjust_indices()
     if (search != model::universe_map.end()) {
       *it = search->second;
     } else {
-      std::stringstream err_msg;
-      err_msg << "Invalid universe number " << uid << " specified on "
-           "lattice " << id_;
-      fatal_error(err_msg);
+      fatal_error(fmt::format(
+        "Invalid universe number {} specified on lattice {}", uid, id_));
     }
   }
 
@@ -84,10 +84,8 @@ Lattice::adjust_indices()
     if (search != model::universe_map.end()) {
       outer_ = search->second;
     } else {
-      std::stringstream err_msg;
-      err_msg << "Invalid universe number " << outer_ << " specified on "
-           "lattice " << id_;
-      fatal_error(err_msg);
+      fatal_error(fmt::format(
+        "Invalid universe number {} specified on lattice {}", outer_, id_));
     }
   }
 }
@@ -184,12 +182,9 @@ RectLattice::RectLattice(pugi::xml_node lat_node)
   std::string univ_str {get_node_value(lat_node, "universes")};
   std::vector<std::string> univ_words {split(univ_str)};
   if (univ_words.size() != nx*ny*nz) {
-    std::stringstream err_msg;
-    err_msg << "Expected " << nx*ny*nz
-            << " universes for a rectangular lattice of size "
-            << nx << "x" << ny << "x" << nz << " but " << univ_words.size()
-            << " were specified.";
-    fatal_error(err_msg);
+    fatal_error(fmt::format(
+      "Expected {} universes for a rectangular lattice of size {}x{]x{} but {} "
+      "were specified.", nx*ny*nz,  nx, ny, nz, univ_words.size()));
   }
 
   // Parse the universes.
@@ -487,12 +482,10 @@ HexLattice::HexLattice(pugi::xml_node lat_node)
   std::string univ_str {get_node_value(lat_node, "universes")};
   std::vector<std::string> univ_words {split(univ_str)};
   if (univ_words.size() != n_univ) {
-    std::stringstream err_msg;
-    err_msg << "Expected " << n_univ
-            << " universes for a hexagonal lattice with " << n_rings_
-            << " rings and " << n_axial_ << " axial levels" << " but "
-            << univ_words.size() << " were specified.";
-    fatal_error(err_msg);
+    fatal_error(fmt::format(
+      "Expected {} universes for a hexagonal lattice with {} rings and {} "
+      "axial levels but {} were specified.", n_univ, n_rings_, n_axial_,
+      univ_words.size()));
   }
 
   // Parse the universes.
@@ -1069,9 +1062,8 @@ void read_lattices(pugi::xml_node node)
     if (in_map == model::lattice_map.end()) {
       model::lattice_map[id] = i_lat;
     } else {
-      std::stringstream err_msg;
-      err_msg << "Two or more lattices use the same unique ID: " << id;
-      fatal_error(err_msg);
+      fatal_error(fmt::format(
+        "Two or more lattices use the same unique ID: {}", id));
     }
   }
 }
diff --git a/src/mgxs.cpp b/src/mgxs.cpp
index 55be267bc..094361dcb 100644
--- a/src/mgxs.cpp
+++ b/src/mgxs.cpp
@@ -9,6 +9,7 @@
 #include <omp.h>
 #endif
 
+#include <fmt/format.h>
 #include "xtensor/xmath.hpp"
 #include "xtensor/xsort.hpp"
 #include "xtensor/xadapt.hpp"
@@ -121,10 +122,9 @@ Mgxs::metadata_from_hdf5(hid_t xs_id, const std::vector<double>& temperature,
             temps_to_read.push_back(std::round(temp_actual));
           }
         } else {
-          std::stringstream msg;
-          msg << "MGXS library does not contain cross sections for "
-            << in_name << " at or near " << std::round(T) << " K.";
-          fatal_error(msg);
+          fatal_error(fmt::format(
+            "MGXS library does not contain cross sections for {} at or near {} K.",
+            in_name, std::round(T)));
         }
       }
       break;
@@ -350,10 +350,9 @@ Mgxs::Mgxs(const std::string& in_name, const std::vector<double>& mat_kTs,
           auto temp_actual = micros[m]->kTs[micro_t[m]];
 
           if (std::abs(temp_actual - temp_desired) >= K_BOLTZMANN * settings::temperature_tolerance) {
-            std::stringstream msg;
-            msg << "MGXS Library does not contain cross section for " << name
-              << " at or near " << std::round(temp_desired / K_BOLTZMANN) << "K.";
-            fatal_error(msg);
+            fatal_error(fmt::format(
+              "MGXS Library does not contain cross section for {} at or near {} K.",
+              name, std::round(temp_desired / K_BOLTZMANN)));
           }
         }
         break;
diff --git a/src/output.cpp b/src/output.cpp
index f32d4c259..ecfdb8a95 100644
--- a/src/output.cpp
+++ b/src/output.cpp
@@ -12,6 +12,7 @@
 #include <utility> // for pair
 
 #include <fmt/format.h>
+#include <fmt/ostream.h>
 #ifdef _OPENMP
 #include <omp.h>
 #endif
@@ -184,10 +185,8 @@ extern "C" void print_particle(Particle* p)
         p->coord_[i].lattice_y, p->coord_[i].lattice_z);
     }
 
-    fmt::print("    r = ");
-    std::cout << p->coord_[i].r << '\n';
-    fmt::print("    u = ");
-    std::cout << p->coord_[i].u << '\n';
+    fmt::print("    r = {}\n", p->coord_[i].r);
+    fmt::print("    u = {}\n", p->coord_[i].u);
   }
 
   // Display miscellaneous info.
@@ -586,7 +585,6 @@ write_tallies()
   // Open the tallies.out file.
   std::ofstream tallies_out;
   tallies_out.open("tallies.out", std::ios::out | std::ios::trunc);
-  tallies_out << std::setprecision(6);
 
   // Loop over each tally.
   for (auto i_tally = 0; i_tally < model::tallies.size(); ++i_tally) {
@@ -595,10 +593,10 @@ write_tallies()
     // Write header block.
     std::string tally_header("TALLY " + std::to_string(tally.id_));
     if (!tally.name_.empty()) tally_header += ": " + tally.name_;
-    tallies_out << header(tally_header) << "\n\n";
+    fmt::print(tallies_out, "{}\n\n", header(tally_header));
 
     if (!tally.writable_) {
-      tallies_out << " Internal\n\n";
+      fmt::print(tallies_out, " Internal\n\n");
       continue;
     }
 
@@ -614,21 +612,20 @@ write_tallies()
       const auto& deriv {model::tally_derivs[tally.deriv_]};
       switch (deriv.variable) {
       case DerivativeVariable::DENSITY:
-        tallies_out << " Density derivative Material "
-          << std::to_string(deriv.diff_material) << "\n";
+        fmt::print(tallies_out, " Density derivative Material {}\n",
+          deriv.diff_material);
         break;
       case DerivativeVariable::NUCLIDE_DENSITY:
-        tallies_out << " Nuclide density derivative Material "
-          << std::to_string(deriv.diff_material) << "  Nuclide "
-          << data::nuclides[deriv.diff_nuclide]->name_ << "\n";
+        fmt::print(tallies_out, " Nuclide density derivative Material {} Nuclide {}\n",
+          deriv.diff_material, data::nuclides[deriv.diff_nuclide]->name_);
         break;
       case DerivativeVariable::TEMPERATURE:
-        tallies_out << " Temperature derivative Material "
-          << std::to_string(deriv.diff_material) << "\n";
+        fmt::print(tallies_out, " Temperature derivative Material {}\n",
+          deriv.diff_material);
         break;
       default:
-        fatal_error("Differential tally dependent variable for tally "
-          + std::to_string(tally.id_) + " not defined in output.cpp");
+        fatal_error(fmt::format("Differential tally dependent variable for "
+          "tally {} not defined in output.cpp", tally.id_));
       }
     }
 
@@ -651,8 +648,8 @@ write_tallies()
           auto i_filt = tally.filters(i);
           const auto& filt {*model::tally_filters[i_filt]};
           auto& match {filter_matches[i_filt]};
-          tallies_out << std::string(indent+1, ' ')
-            << filt.text_label(match.i_bin_) << "\n";
+          fmt::print(tallies_out, "{0:{1}}{2}\n", "", indent + 1,
+            filt.text_label(match.i_bin_));
         }
         indent += 2;
       }
@@ -662,14 +659,14 @@ write_tallies()
       for (auto i_nuclide : tally.nuclides_) {
         // Write label for this nuclide bin.
         if (i_nuclide == -1) {
-          tallies_out << std::string(indent+1, ' ') << "Total Material\n";
+          fmt::print(tallies_out, "{0:{1}}Total Material\n", "", indent + 1);
         } else {
           if (settings::run_CE) {
-            tallies_out << std::string(indent+1, ' ')
-              << data::nuclides[i_nuclide]->name_ << "\n";
+            fmt::print(tallies_out, "{0:{1}}{2}\n", "", indent + 1,
+              data::nuclides[i_nuclide]->name_);
           } else {
-            tallies_out << std::string(indent+1, ' ')
-              << data::mg.nuclides_[i_nuclide].name << "\n";
+            fmt::print(tallies_out, "{0:{1}}{2}\n", "", indent + 1,
+              data::mg.nuclides_[i_nuclide].name);
           }
         }
 
@@ -681,9 +678,8 @@ write_tallies()
           double mean, stdev;
           std::tie(mean, stdev) = mean_stdev(
             &tally.results_(filter_index, score_index, 0), tally.n_realizations_);
-          tallies_out << std::string(indent+1, ' ')  << std::left
-            << std::setw(36) << score_name << " " << mean << " +/- "
-            << t_value * stdev << "\n";
+          fmt::print(tallies_out, "{0:{1}}{2:<36} {3:.6} +/- {4:.6}\n",
+            "", indent + 1, score_name, mean, t_value * stdev);
           score_index += 1;
         }
         indent -= 2;
diff --git a/src/particle.cpp b/src/particle.cpp
index bd9a61139..276fc197d 100644
--- a/src/particle.cpp
+++ b/src/particle.cpp
@@ -2,7 +2,8 @@
 
 #include <algorithm> // copy, min
 #include <cmath>     // log, abs, copysign
-#include <sstream>
+
+#include <fmt/format.h>
 
 #include "openmc/bank.h"
 #include "openmc/capi.h"
@@ -643,14 +644,13 @@ Particle::write_restart() const
   if (settings::run_mode == RunMode::PARTICLE) return;
 
   // Set up file name
-  std::stringstream filename;
-  filename << settings::path_output << "particle_" << simulation::current_batch
-    << '_' << id_ << ".h5";
+  auto filename = fmt::format("{}particle_{}_{}.h5", settings::path_output,
+    simulation::current_batch, id_);
 
   #pragma omp critical (WriteParticleRestart)
   {
     // Create file
-    hid_t file_id = file_open(filename.str(), 'w');
+    hid_t file_id = file_open(filename, 'w');
 
     // Write filetype and version info
     write_attribute(file_id, "filetype", "particle restart");
diff --git a/src/physics.cpp b/src/physics.cpp
index 6c46bf429..a0f929125 100644
--- a/src/physics.cpp
+++ b/src/physics.cpp
@@ -21,9 +21,10 @@
 #include "openmc/thermal.h"
 #include "openmc/tallies/tally.h"
 
+#include <fmt/format.h>
+
 #include <algorithm> // for max, min, max_element
 #include <cmath> // for sqrt, exp, log, abs, copysign
-#include <sstream>
 
 namespace openmc {
 
@@ -61,17 +62,19 @@ void collision(Particle* p)
 
   // Display information about collision
   if (settings::verbosity >= 10 || p->trace_) {
-    std::stringstream msg;
+    std::string msg;
     if (p->event_ == TallyEvent::KILL) {
-      msg << "    Killed. Energy = " << p->E_ << " eV.";
+      msg = fmt::format("    Killed. Energy = {} eV.", p->E_);
     } else if (p->type_ == Particle::Type::neutron) {
-      msg << "    " << reaction_name(p->event_mt_) << " with " <<
-        data::nuclides[p->event_nuclide_]->name_ << ". Energy = " << p->E_ << " eV.";
+      msg = fmt::format("    {} with {}. Energy = {} eV.",
+        reaction_name(p->event_mt_), data::nuclides[p->event_nuclide_]->name_,
+        p->E_);
     } else if (p->type_ == Particle::Type::photon) {
-      msg << "    " << reaction_name(p->event_mt_) << " with " <<
-        to_element(data::nuclides[p->event_nuclide_]->name_) << ". Energy = " << p->E_ << " eV.";
+      msg = fmt::format("    {} with {}. Energy = {} eV.",
+        reaction_name(p->event_mt_),
+        to_element(data::nuclides[p->event_nuclide_]->name_), p->E_);
     } else {
-      msg << "    Disappeared. Energy = " << p->E_ << " eV.";
+      msg = fmt::format("    Disappeared. Energy = {} eV.", p->E_);
     }
     write_message(msg, 1);
   }
@@ -189,7 +192,7 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx)
 
     // Sample delayed group and angle/energy for fission reaction
     sample_fission_neutron(i_nuclide, rx, p->E_, &site, p->current_seed());
-    
+
     // Store fission site in bank
     if (use_fission_bank) {
       int64_t idx = simulation::fission_bank.thread_safe_append(site);
@@ -210,7 +213,7 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx)
     if (p->delayed_group_ > 0) {
       nu_d[p->delayed_group_-1]++;
     }
-    
+
     // Write fission particles to nuBank
     if (use_fission_bank) {
       p->nu_bank_.emplace_back();
@@ -220,7 +223,7 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx)
       nu_bank_entry->delayed_group    = site.delayed_group;
     }
   }
-  
+
   // If shared fission bank was full, and no fissions could be added,
   // set the particle fission flag to false.
   if (nu == skipped) {
diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp
index 1ac7c1979..526a21439 100644
--- a/src/physics_mg.cpp
+++ b/src/physics_mg.cpp
@@ -1,8 +1,8 @@
 #include "openmc/physics_mg.h"
 
 #include <stdexcept>
-#include <sstream>
 
+#include <fmt/format.h>
 #include "xtensor/xarray.hpp"
 
 #include "openmc/bank.h"
@@ -31,10 +31,8 @@ collision_mg(Particle* p)
   sample_reaction(p);
 
   // Display information about collision
-  if ((settings::verbosity >= 10) || (p->trace_)) {
-    std::stringstream msg;
-    msg << "    Energy Group = " << p->g_;
-    write_message(msg, 1);
+  if ((settings::verbosity >= 10) || p->trace_) {
+    write_message(fmt::format("    Energy Group = {}", p->g_), 1);
   }
 }
 
@@ -113,13 +111,13 @@ create_fission_sites(Particle* p)
   // Initialize the counter of delayed neutrons encountered for each delayed
   // group.
   double nu_d[MAX_DELAYED_GROUPS] = {0.};
-  
+
   // Clear out particle's nu fission bank
   p->nu_bank_.clear();
 
   p->fission_ = true;
   int skipped = 0;
-  
+
   // Determine whether to place fission sites into the shared fission bank
   // or the secondary particle bank.
   bool use_fission_bank = (settings::run_mode == RunMode::EIGENVALUE);
@@ -176,7 +174,7 @@ create_fission_sites(Particle* p)
     if (p->delayed_group_ > 0) {
       nu_d[dg]++;
     }
-    
+
     // Write fission particles to nuBank
     if (use_fission_bank) {
       p->nu_bank_.emplace_back();
@@ -186,7 +184,7 @@ create_fission_sites(Particle* p)
       nu_bank_entry->delayed_group    = site.delayed_group;
     }
   }
-  
+
   // If shared fission bank was full, and no fissions could be added,
   // set the particle fission flag to false.
   if (nu == skipped) {
diff --git a/src/plot.cpp b/src/plot.cpp
index 5aa6a61b8..0ff785712 100644
--- a/src/plot.cpp
+++ b/src/plot.cpp
@@ -4,6 +4,8 @@
 #include <fstream>
 #include <sstream>
 
+#include <fmt/format.h>
+#include <fmt/ostream.h>
 #include "xtensor/xview.hpp"
 
 #include "openmc/constants.h"
@@ -92,10 +94,8 @@ extern "C"
 int openmc_plot_geometry()
 {
   for (auto pl : model::plots) {
-    std::stringstream ss;
-    ss << "Processing plot " << pl.id_ << ": "
-       << pl.path_plot_ << "...";
-    write_message(ss.str(), 5);
+    write_message(fmt::format("Processing plot {}: {}...",
+      pl.id_, pl.path_plot_), 5);
 
     if (PlotType::slice == pl.type_) {
       // create 2D image
@@ -188,9 +188,7 @@ Plot::set_id(pugi::xml_node plot_node)
 
   // Check to make sure 'id' hasn't been used
   if (model::plot_map.find(id_) != model::plot_map.end()) {
-    std::stringstream err_msg;
-    err_msg << "Two or more plots use the same unique ID: " << id_;
-    fatal_error(err_msg.str());
+    fatal_error(fmt::format("Two or more plots use the same unique ID: {}", id_));
   }
 }
 
@@ -210,11 +208,9 @@ Plot::set_type(pugi::xml_node plot_node)
     else if (type_str == "voxel") {
       type_ = PlotType::voxel;
     } else {
-    // if we're here, something is wrong
-    std::stringstream err_msg;
-    err_msg << "Unsupported plot type '" << type_str
-            << "' in plot " << id_;
-    fatal_error(err_msg.str());
+      // if we're here, something is wrong
+      fatal_error(fmt::format("Unsupported plot type '{}' in plot {}",
+        type_str, id_));
     }
   }
 }
@@ -223,24 +219,24 @@ void
 Plot::set_output_path(pugi::xml_node plot_node)
 {
   // Set output file path
-  std::stringstream filename;
+  std::string filename;
 
   if (check_for_node(plot_node, "filename")) {
-    filename << get_node_value(plot_node, "filename");
+    filename = get_node_value(plot_node, "filename");
   } else {
-    filename << "plot_" << id_;
+    filename = fmt::format("plot_{}", id_);
   }
   // add appropriate file extension to name
   switch(type_) {
   case PlotType::slice:
-    filename << ".ppm";
+    filename.append(".ppm");
     break;
   case PlotType::voxel:
-    filename << ".h5";
+    filename.append(".h5");
     break;
   }
 
-  path_plot_ = filename.str();
+  path_plot_ = filename;
 
   // Copy plot pixel size
   std::vector<int> pxls = get_node_array<int>(plot_node, "pixels");
@@ -249,10 +245,7 @@ Plot::set_output_path(pugi::xml_node plot_node)
       pixels_[0] = pxls[0];
       pixels_[1] = pxls[1];
     } else {
-      std::stringstream err_msg;
-      err_msg << "<pixels> must be length 2 in slice plot "
-              << id_;
-      fatal_error(err_msg.str());
+      fatal_error(fmt::format("<pixels> must be length 2 in slice plot {}", id_));
     }
   } else if (PlotType::voxel == type_) {
     if (pxls.size() == 3) {
@@ -260,10 +253,7 @@ Plot::set_output_path(pugi::xml_node plot_node)
       pixels_[1] = pxls[1];
       pixels_[2] = pxls[2];
     } else {
-      std::stringstream err_msg;
-      err_msg << "<pixels> must be length 3 in voxel plot "
-              << id_;
-      fatal_error(err_msg.str());
+      fatal_error(fmt::format("<pixels> must be length 3 in voxel plot {}", id_));
     }
   }
 }
@@ -276,19 +266,13 @@ Plot::set_bg_color(pugi::xml_node plot_node)
     std::vector<int> bg_rgb = get_node_array<int>(plot_node, "background");
     if (PlotType::voxel == type_) {
       if (mpi::master) {
-        std::stringstream err_msg;
-        err_msg << "Background color ignored in voxel plot "
-                << id_;
-        warning(err_msg.str());
+        warning(fmt::format("Background color ignored in voxel plot {}", id_));
       }
     }
     if (bg_rgb.size() == 3) {
       not_found_ = bg_rgb;
     } else {
-      std::stringstream err_msg;
-      err_msg << "Bad background RGB in plot "
-              << id_;
-      fatal_error(err_msg);
+      fatal_error(fmt::format("Bad background RGB in plot {}", id_));
     }
   }
 }
@@ -309,10 +293,8 @@ Plot::set_basis(pugi::xml_node plot_node)
     } else if ("yz" == pl_basis) {
       basis_ = PlotBasis::yz;
     } else {
-      std::stringstream err_msg;
-      err_msg << "Unsupported plot basis '" << pl_basis
-              << "' in plot " << id_;
-      fatal_error(err_msg);
+      fatal_error(fmt::format("Unsupported plot basis '{}' in plot {}",
+        pl_basis, id_));
     }
   }
 }
@@ -325,10 +307,7 @@ Plot::set_origin(pugi::xml_node plot_node)
   if (pl_origin.size() == 3) {
     origin_ = pl_origin;
   } else {
-    std::stringstream err_msg;
-    err_msg << "Origin must be length 3 in plot "
-            << id_;
-    fatal_error(err_msg);
+    fatal_error(fmt::format("Origin must be length 3 in plot {}", id_));
   }
 }
 
@@ -342,20 +321,14 @@ Plot::set_width(pugi::xml_node plot_node)
       width_.x = pl_width[0];
       width_.y = pl_width[1];
     } else {
-      std::stringstream err_msg;
-      err_msg << "<width> must be length 2 in slice plot "
-              << id_;
-      fatal_error(err_msg);
+      fatal_error(fmt::format("<width> must be length 2 in slice plot {}", id_));
     }
   } else if (PlotType::voxel == type_) {
     if (pl_width.size() == 3) {
       pl_width = get_node_array<double>(plot_node, "width");
       width_ = pl_width;
     } else {
-      std::stringstream err_msg;
-      err_msg << "<width> must be length 3 in voxel plot "
-              << id_;
-      fatal_error(err_msg);
+      fatal_error(fmt::format("<width> must be length 3 in voxel plot {}", id_));
     }
   }
 }
@@ -367,9 +340,7 @@ Plot::set_universe(pugi::xml_node plot_node)
   if (check_for_node(plot_node, "level")) {
     level_ = std::stoi(get_node_value(plot_node, "level"));
     if (level_ < 0) {
-      std::stringstream err_msg;
-      err_msg << "Bad universe level in plot " << id_;
-      fatal_error(err_msg);
+      fatal_error(fmt::format("Bad universe level in plot {}", id_));
     }
   } else {
     level_ = PLOT_LEVEL_LOWEST;
@@ -391,10 +362,8 @@ Plot::set_default_colors(pugi::xml_node plot_node)
     color_by_ = PlotColorBy::mats;
     colors_.resize(model::materials.size());
   } else {
-    std::stringstream err_msg;
-    err_msg << "Unsupported plot color type '" << pl_color_by
-            << "' in plot " << id_;
-    fatal_error(err_msg);
+    fatal_error(fmt::format("Unsupported plot color type '{}' in plot {}",
+      pl_color_by, id_));
   }
 
   for (auto& c : colors_) {
@@ -411,10 +380,7 @@ Plot::set_user_colors(pugi::xml_node plot_node)
 {
   if (!plot_node.select_nodes("color").empty() && PlotType::voxel == type_) {
     if (mpi::master) {
-      std::stringstream err_msg;
-      err_msg << "Color specifications ignored in voxel plot "
-              << id_;
-      warning(err_msg);
+      warning(fmt::format("Color specifications ignored in voxel plot {}", id_));
     }
   }
 
@@ -422,19 +388,15 @@ Plot::set_user_colors(pugi::xml_node plot_node)
     // Make sure 3 values are specified for RGB
     std::vector<int> user_rgb = get_node_array<int>(cn, "rgb");
     if (user_rgb.size() != 3) {
-      std::stringstream err_msg;
-      err_msg << "Bad RGB in plot " << id_;
-      fatal_error(err_msg);
+      fatal_error(fmt::format("Bad RGB in plot {}", id_));
     }
     // Ensure that there is an id for this color specification
     int col_id;
     if (check_for_node(cn, "id")) {
       col_id = std::stoi(get_node_value(cn, "id"));
     } else {
-      std::stringstream err_msg;
-      err_msg << "Must specify id for color specification in plot "
-              << id_;
-      fatal_error(err_msg);
+      fatal_error(fmt::format(
+        "Must specify id for color specification in plot {}", id_));
     }
     // Add RGB
     if (PlotColorBy::cells == color_by_) {
@@ -442,20 +404,16 @@ Plot::set_user_colors(pugi::xml_node plot_node)
         col_id = model::cell_map[col_id];
         colors_[col_id] = user_rgb;
       } else {
-        std::stringstream err_msg;
-        err_msg << "Could not find cell " << col_id
-                << " specified in plot " << id_;
-        fatal_error(err_msg);
+        fatal_error(fmt::format("Could not find cell {} specified in plot {}",
+          col_id, id_));
       }
     } else if (PlotColorBy::mats == color_by_) {
       if (model::material_map.find(col_id) != model::material_map.end()) {
         col_id = model::material_map[col_id];
         colors_[col_id] = user_rgb;
       } else {
-        std::stringstream err_msg;
-        err_msg << "Could not find material " << col_id
-                << " specified in plot " << id_;
-        fatal_error(err_msg);
+        fatal_error(fmt::format(
+          "Could not find material {} specified in plot {}", col_id, id_));
       }
     }
   } // color node loop
@@ -469,9 +427,7 @@ Plot::set_meshlines(pugi::xml_node plot_node)
 
   if (!mesh_line_nodes.empty()) {
     if (PlotType::voxel == type_) {
-      std::stringstream msg;
-      msg << "Meshlines ignored in voxel plot " << id_;
-      warning(msg);
+      warning(fmt::format("Meshlines ignored in voxel plot {}", id_));
     }
 
     if (mesh_line_nodes.size() == 1) {
@@ -483,9 +439,8 @@ Plot::set_meshlines(pugi::xml_node plot_node)
       if (check_for_node(meshlines_node, "meshtype")) {
         meshtype = get_node_value(meshlines_node, "meshtype");
       } else {
-        std::stringstream err_msg;
-        err_msg << "Must specify a meshtype for meshlines specification in plot " << id_;
-        fatal_error(err_msg);
+        fatal_error(fmt::format(
+          "Must specify a meshtype for meshlines specification in plot {}", id_));
       }
 
       // Ensure that there is a linewidth for this meshlines specification
@@ -494,9 +449,8 @@ Plot::set_meshlines(pugi::xml_node plot_node)
         meshline_width = get_node_value(meshlines_node, "linewidth");
         meshlines_width_ = std::stoi(meshline_width);
       } else {
-        std::stringstream err_msg;
-        err_msg << "Must specify a linewidth for meshlines specification in plot " << id_;
-        fatal_error(err_msg);
+        fatal_error(fmt::format(
+          "Must specify a linewidth for meshlines specification in plot {}", id_));
       }
 
       // Check for color
@@ -504,9 +458,7 @@ Plot::set_meshlines(pugi::xml_node plot_node)
         // Check and make sure 3 values are specified for RGB
         std::vector<int> ml_rgb = get_node_array<int>(meshlines_node, "color");
         if (ml_rgb.size() != 3) {
-          std::stringstream err_msg;
-          err_msg << "Bad RGB for meshlines color in plot " << id_;
-          fatal_error(err_msg);
+          fatal_error(fmt::format("Bad RGB for meshlines color in plot {}", id_));
         }
         meshlines_color_ = ml_rgb;
       }
@@ -514,9 +466,7 @@ Plot::set_meshlines(pugi::xml_node plot_node)
       // Set mesh based on type
       if ("ufs" == meshtype) {
         if (!simulation::ufs_mesh) {
-          std::stringstream err_msg;
-          err_msg << "No UFS mesh for meshlines on plot " << id_;
-          fatal_error(err_msg);
+          fatal_error(fmt::format("No UFS mesh for meshlines on plot {}", id_));
         } else {
           for (int i = 0; i < model::meshes.size(); ++i) {
             if (const auto* m
@@ -531,9 +481,7 @@ Plot::set_meshlines(pugi::xml_node plot_node)
         }
       } else if ("entropy" == meshtype) {
         if (!simulation::entropy_mesh) {
-          std::stringstream err_msg;
-          err_msg << "No entropy mesh for meshlines on plot " << id_;
-          fatal_error(err_msg);
+          fatal_error(fmt::format("No entropy mesh for meshlines on plot {}", id_));
         } else {
           for (int i = 0; i < model::meshes.size(); ++i) {
             if (const auto* m
@@ -553,29 +501,22 @@ Plot::set_meshlines(pugi::xml_node plot_node)
           tally_mesh_id = std::stoi(get_node_value(meshlines_node, "id"));
         } else {
           std::stringstream err_msg;
-          err_msg << "Must specify a mesh id for meshlines tally "
-                  << "mesh specification in plot " << id_;
-          fatal_error(err_msg);
+          fatal_error(fmt::format("Must specify a mesh id for meshlines tally "
+            "mesh specification in plot {}", id_));
         }
         // find the tally index
         int idx;
         int err = openmc_get_mesh_index(tally_mesh_id, &idx);
         if (err != 0) {
-          std::stringstream err_msg;
-          err_msg << "Could not find mesh " << tally_mesh_id
-                  << " specified in meshlines for plot " << id_;
-          fatal_error(err_msg);
+          fatal_error(fmt::format("Could not find mesh {} specified in "
+            "meshlines for plot {}", tally_mesh_id, id_));
         }
         index_meshlines_mesh_ = idx;
       } else {
-        std::stringstream err_msg;
-        err_msg << "Invalid type for meshlines on plot " << id_ ;
-        fatal_error(err_msg);
+        fatal_error(fmt::format("Invalid type for meshlines on plot {}", id_ ));
       }
     } else {
-      std::stringstream err_msg;
-      err_msg << "Mutliple meshlines specified in plot " << id_;
-      fatal_error(err_msg);
+      fatal_error(fmt::format("Mutliple meshlines specified in plot {}", id_));
     }
   }
 }
@@ -589,9 +530,7 @@ Plot::set_mask(pugi::xml_node plot_node)
   if (!mask_nodes.empty()) {
     if (PlotType::voxel == type_) {
       if (mpi::master) {
-        std::stringstream wrn_msg;
-        wrn_msg << "Mask ignored in voxel plot " << id_;
-        warning(wrn_msg);
+        warning(fmt::format("Mask ignored in voxel plot {}", id_));
       }
     }
 
@@ -602,9 +541,7 @@ Plot::set_mask(pugi::xml_node plot_node)
       // Determine how many components there are and allocate
       std::vector<int> iarray = get_node_array<int>(mask_node, "components");
       if (iarray.size() == 0) {
-        std::stringstream err_msg;
-        err_msg << "Missing <components> in mask of plot " << id_;
-        fatal_error(err_msg);
+        fatal_error(fmt::format("Missing <components> in mask of plot {}", id_));
       }
 
       // First we need to change the user-specified identifiers to indices
@@ -615,20 +552,16 @@ Plot::set_mask(pugi::xml_node plot_node)
             col_id  = model::cell_map[col_id];
           }
           else {
-            std::stringstream err_msg;
-            err_msg << "Could not find cell " << col_id
-                    << " specified in the mask in plot " << id_;
-            fatal_error(err_msg);
+            fatal_error(fmt::format("Could not find cell {} specified in the "
+              "mask in plot {}", col_id, id_));
           }
         } else if (PlotColorBy::mats == color_by_) {
           if (model::material_map.find(col_id) != model::material_map.end()) {
             col_id = model::material_map[col_id];
           }
           else {
-            std::stringstream err_msg;
-            err_msg << "Could not find material " << col_id
-                    << " specified in the mask in plot " << id_;
-            fatal_error(err_msg);
+            fatal_error(fmt::format("Could not find material {} specified in "
+              "the mask in plot {}", col_id, id_));
           }
         }
       }
@@ -646,9 +579,7 @@ Plot::set_mask(pugi::xml_node plot_node)
       }
 
     } else {
-      std::stringstream err_msg;
-      err_msg << "Mutliple masks specified in plot " << id_;
-      fatal_error(err_msg);
+      fatal_error(fmt::format("Mutliple masks specified in plot {}", id_));
     }
   }
 }
@@ -660,18 +591,14 @@ void Plot::set_overlap_color(pugi::xml_node plot_node) {
     // check for custom overlap color
     if (check_for_node(plot_node, "overlap_color")) {
       if (!color_overlaps_) {
-        std::stringstream wrn_msg;
-        wrn_msg << "Overlap color specified in plot " << id_
-                << " but overlaps won't be shown.";
-        warning(wrn_msg);
+        warning(fmt::format(
+          "Overlap color specified in plot {} but overlaps won't be shown.", id_));
       }
       std::vector<int> olap_clr = get_node_array<int>(plot_node, "overlap_color");
       if (olap_clr.size() == 3) {
         overlap_color_ = olap_clr;
       } else {
-        std::stringstream err_msg;
-        err_msg << "Bad overlap RGB in plot " << id_;
-        fatal_error(err_msg);
+        fatal_error(fmt::format("Bad overlap RGB in plot {}", id_));
       }
     }
   }
@@ -716,9 +643,9 @@ void output_ppm(Plot pl, const ImageData& data)
   of.open(fname);
 
   // Write header
-  of << "P6" << "\n";
+  of << "P6\n";
   of << pl.pixels_[0] << " " << pl.pixels_[1] << "\n";
-  of << "255" << "\n";
+  of << "255\n";
   of.close();
 
   of.open(fname, std::ios::binary | std::ios::app);
@@ -729,11 +656,7 @@ void output_ppm(Plot pl, const ImageData& data)
       of << rgb.red << rgb.green << rgb.blue;
     }
   }
-
-  // Close file
-  // THIS IS HERE TO MATCH FORTRAN VERSION, NOT TECHNICALLY NECESSARY
   of << "\n";
-  of.close();
 }
 
 //==============================================================================
diff --git a/src/reaction.cpp b/src/reaction.cpp
index f7f1e87ac..ff22540e1 100644
--- a/src/reaction.cpp
+++ b/src/reaction.cpp
@@ -3,6 +3,8 @@
 #include <string>
 #include <utility> // for move
 
+#include <fmt/format.h>
+
 #include "openmc/constants.h"
 #include "openmc/hdf5_interface.h"
 #include "openmc/endf.h"
@@ -35,8 +37,7 @@ Reaction::Reaction(hid_t group, const std::vector<int>& temperatures)
   // Read cross section and threshold_idx data
   for (auto t : temperatures) {
     // Get group corresponding to temperature
-    std::string temp_str {std::to_string(t) + "K"};
-    hid_t temp_group = open_group(group, temp_str.c_str());
+    hid_t temp_group = open_group(group, fmt::format("{}K", t).c_str());
     hid_t dset = open_dataset(temp_group, "xs");
 
     // Get threshold index
@@ -178,7 +179,7 @@ std::string reaction_name(int mt)
   } else if (mt == N_NPA) {
     return "(n,npa)";
   } else if (N_N1 <= mt && mt <= N_N40) {
-    return "(n,n" + std::to_string(mt-50) + ")";
+    return fmt::format("(n,n{})", mt - 50);
   } else if (mt == N_NC) {
     return "(n,nc)";
   } else if (mt == N_DISAPPEAR) {
@@ -244,33 +245,31 @@ std::string reaction_name(int mt)
   } else if (mt == PHOTOELECTRIC) {
     return "photoelectric";
   } else if (534 <= mt && mt <= 572) {
-    std::stringstream name;
-    name << "photoelectric, " << SUBSHELLS[mt - 534] << " subshell";
-    return name.str();
+    return fmt::format("photoelectric, {} subshell", SUBSHELLS[mt - 534]);
   } else if (600 <= mt && mt <= 648) {
-    return "(n,p" + std::to_string(mt-600) + ")";
+    return fmt::format("(n,p{})", mt - 600);
   } else if (mt == 649) {
     return "(n,pc)";
   } else if (650 <= mt && mt <= 698) {
-    return "(n,d" + std::to_string(mt-650) + ")";
+    return fmt::format("(n,d{})", mt - 650);
   } else if (mt == 699) {
     return "(n,dc)";
   } else if (700 <= mt && mt <= 748) {
-    return "(n,t" + std::to_string(mt-700) + ")";
+    return fmt::format("(n,t{})", mt - 700);
   } else if (mt == 749) {
     return "(n,tc)";
   } else if (750 <= mt && mt <= 798) {
-    return "(n,3He" + std::to_string(mt-750) + ")";
+    return fmt::format("(n,3He{})", mt - 750);
   } else if (mt == 799) {
     return "(n,3Hec)";
   } else if (800 <= mt && mt <= 848) {
-    return "(n,a" + std::to_string(mt-800) + ")";
+    return fmt::format("(n,a{})", mt - 800);
   } else if (mt == 849) {
     return "(n,ac)";
   } else if (mt == HEATING_LOCAL) {
     return "heating-local";
   } else {
-    return "MT=" + std::to_string(mt);
+    return fmt::format("MT={}", mt);
   }
 }
 
diff --git a/src/secondary_uncorrelated.cpp b/src/secondary_uncorrelated.cpp
index 1421fcf83..320761916 100644
--- a/src/secondary_uncorrelated.cpp
+++ b/src/secondary_uncorrelated.cpp
@@ -1,8 +1,9 @@
 #include "openmc/secondary_uncorrelated.h"
 
-#include <sstream> // for stringstream
 #include <string>  // for string
 
+#include <fmt/format.h>
+
 #include "openmc/error.h"
 #include "openmc/hdf5_interface.h"
 #include "openmc/random_lcg.h"
@@ -42,9 +43,7 @@ UncorrelatedAngleEnergy::UncorrelatedAngleEnergy(hid_t group)
     } else if (type == "watt") {
       energy_ = UPtrEDist{new WattEnergy{energy_group}};
     } else {
-      std::stringstream msg;
-      msg << "Energy distribution type '" << type << "' not implemented.";
-      warning(msg);
+      warning(fmt::format("Energy distribution type '{}' not implemented.", type));
     }
     close_group(energy_group);
   }