diff --git a/include/openmc/hdf5_interface.h b/include/openmc/hdf5_interface.h
index ee3419fa2..ddb2319ae 100644
--- a/include/openmc/hdf5_interface.h
+++ b/include/openmc/hdf5_interface.h
@@ -26,13 +26,18 @@ namespace openmc {
 //==============================================================================
 
 void read_attr(hid_t obj_id, const char* name, hid_t mem_type_id,
-               void* buffer);
+  void* buffer);
+
 void write_attr(hid_t obj_id, int ndim, const hsize_t* dims, const char* name,
-                hid_t mem_type_id, const void* buffer);
-void read_dataset(hid_t obj_id, const char* name, hid_t mem_type_id,
-                  void* buffer, bool indep);
-void write_dataset(hid_t group_id, int ndim, const hsize_t* dims, const char* name,
-                   hid_t mem_type_id, const void* buffer, bool indep);
+  hid_t mem_type_id, const void* buffer);
+
+void read_dataset_lowlevel(hid_t obj_id, const char* name, hid_t mem_type_id,
+  hid_t mem_space_id, bool indep, void* buffer);
+
+void write_dataset_lowlevel(hid_t group_id, int ndim, const hsize_t* dims,
+  const char* name, hid_t mem_type_id, hid_t mem_space_id, bool indep,
+  const void* buffer);
+
 bool using_mpio_device(hid_t obj_id);
 
 //==============================================================================
@@ -86,34 +91,32 @@ extern "C" {
   void read_attr_string(hid_t obj_id, const char* name, size_t slen,
                         char* buffer);
   void read_complex(hid_t obj_id, const char* name,
-                             std::complex<double>* buffer, bool indep);
-  void read_double(hid_t obj_id, const char* name, double* buffer,
-                            bool indep);
-  void read_int(hid_t obj_id, const char* name, int* buffer,
-                         bool indep);
+                    std::complex<double>* buffer, bool indep);
+  void read_double(hid_t obj_id, const char* name, double* buffer, bool indep);
+  void read_int(hid_t obj_id, const char* name, int* buffer, bool indep);
   void read_llong(hid_t obj_id, const char* name, long long* buffer,
-                           bool indep);
-  void read_string(hid_t obj_id, const char* name, size_t slen,
-                            char* buffer, bool indep);
+                  bool indep);
+  void read_string(hid_t obj_id, const char* name, size_t slen, char* buffer,
+                   bool indep);
 
 
-  void read_tally_results(hid_t group_id, hsize_t n_filter,
-                                   hsize_t n_score, double* results);
+  void read_tally_results(hid_t group_id, hsize_t n_filter, hsize_t n_score,
+                          double* results);
   void write_attr_double(hid_t obj_id, int ndim, const hsize_t* dims,
-                                  const char* name, const double* buffer);
+                         const char* name, const double* buffer);
   void write_attr_int(hid_t obj_id, int ndim, const hsize_t* dims,
-                               const char* name, const int* buffer);
+                      const char* name, const int* buffer);
   void write_attr_string(hid_t obj_id, const char* name, const char* buffer);
   void write_double(hid_t group_id, int ndim, const hsize_t* dims,
-                             const char* name, const double* buffer, bool indep);
+                    const char* name, const double* buffer, bool indep);
   void write_int(hid_t group_id, int ndim, const hsize_t* dims,
-                          const char* name, const int* buffer, bool indep);
+                 const char* name, const int* buffer, bool indep);
   void write_llong(hid_t group_id, int ndim, const hsize_t* dims,
-                            const char* name, const long long* buffer, bool indep);
+                   const char* name, const long long* buffer, bool indep);
   void write_string(hid_t group_id, int ndim, const hsize_t* dims, size_t slen,
-                             const char* name, char const* buffer, bool indep);
+                    const char* name, char const* buffer, bool indep);
   void write_tally_results(hid_t group_id, hsize_t n_filter, hsize_t n_score,
-                                    const double* results);
+                           const double* results);
 } // extern "C"
 
 //==============================================================================
@@ -233,7 +236,8 @@ template<typename T> inline
 std::enable_if_t<std::is_scalar<std::decay_t<T>>::value>
 read_dataset(hid_t obj_id, const char* name, T& buffer, bool indep=false)
 {
-  read_dataset(obj_id, name, H5TypeMap<T>::type_id, &buffer, indep);
+  read_dataset_lowlevel(obj_id, name, H5TypeMap<T>::type_id, H5S_ALL, indep,
+                        &buffer);
 }
 
 // overload for std::string
@@ -251,9 +255,11 @@ read_dataset(hid_t obj_id, const char* name, std::string& str, bool indep=false)
 
 // array version
 template<typename T, std::size_t N> inline void
-read_dataset(hid_t dset, const char* name, std::array<T, N>& buffer, bool indep=false)
+read_dataset(hid_t dset, const char* name, std::array<T, N>& buffer,
+             bool indep=false)
 {
-  read_dataset(dset, name, H5TypeMap<T>::type_id, buffer.data(), indep);
+  read_dataset_lowlevel(dset, name, H5TypeMap<T>::type_id, H5S_ALL, indep,
+                        buffer.data());
 }
 
 // vector version
@@ -267,11 +273,13 @@ void read_dataset(hid_t dset, std::vector<T>& vec, bool indep=false)
   vec.resize(shape[0]);
 
   // Read data into vector
-  read_dataset(dset, nullptr, H5TypeMap<T>::type_id, vec.data(), indep);
+  read_dataset_lowlevel(dset, nullptr, H5TypeMap<T>::type_id, H5S_ALL, indep,
+                        vec.data());
 }
 
 template <typename T>
-void read_dataset(hid_t obj_id, const char* name, std::vector<T>& vec, bool indep=false)
+void read_dataset(hid_t obj_id, const char* name, std::vector<T>& vec,
+                  bool indep=false)
 {
   hid_t dset = open_dataset(obj_id, name);
   read_dataset(dset, vec, indep);
@@ -291,14 +299,17 @@ void read_dataset(hid_t dset, xt::xarray<T>& arr, bool indep=false)
   arr.resize(shape);
 
   // Read data from attribute
-  read_dataset(dset, nullptr, H5TypeMap<T>::type_id, arr.data(), indep);
+  read_dataset_lowlevel(dset, nullptr, H5TypeMap<T>::type_id, H5S_ALL, indep,
+                        arr.data());
 }
 
 template<>
-void read_dataset(hid_t dset, xt::xarray<std::complex<double>>& arr, bool indep);
+void read_dataset(hid_t dset, xt::xarray<std::complex<double>>& arr,
+                  bool indep);
 
 template <typename T>
-void read_dataset(hid_t obj_id, const char* name, xt::xarray<T>& arr, bool indep=false)
+void read_dataset(hid_t obj_id, const char* name, xt::xarray<T>& arr,
+                  bool indep=false)
 {
   // Open dataset and read array
   hid_t dset = open_dataset(obj_id, name);
@@ -308,7 +319,8 @@ void read_dataset(hid_t obj_id, const char* name, xt::xarray<T>& arr, bool indep
 
 
 template <typename T, std::size_t N>
-void read_dataset(hid_t obj_id, const char* name, xt::xtensor<T, N>& arr, bool indep=false)
+void read_dataset(hid_t obj_id, const char* name, xt::xtensor<T, N>& arr,
+                  bool indep=false)
 {
   // Open dataset and read array
   hid_t dset = open_dataset(obj_id, name);
@@ -346,7 +358,7 @@ read_dataset(hid_t obj_id, const char* name, Position& r, bool indep=false)
 
 template <typename T, std::size_t N>
 inline void read_dataset_as_shape(hid_t obj_id, const char* name,
-                           xt::xtensor<T, N>& arr, bool indep=false)
+                                  xt::xtensor<T, N>& arr, bool indep=false)
 {
   hid_t dset = open_dataset(obj_id, name);
 
@@ -357,7 +369,8 @@ inline void read_dataset_as_shape(hid_t obj_id, const char* name,
   std::vector<T> buffer(size);
 
   // Read data from attribute
-  read_dataset(dset, nullptr, H5TypeMap<T>::type_id, buffer.data(), indep);
+  read_dataset_lowlevel(dset, nullptr, H5TypeMap<T>::type_id, H5S_ALL, indep,
+                        buffer.data());
 
   // Adapt into xarray
   arr = xt::adapt(buffer, arr.shape());
@@ -367,8 +380,8 @@ inline void read_dataset_as_shape(hid_t obj_id, const char* name,
 
 
 template <typename T, std::size_t N>
-inline void read_nd_vector(hid_t obj_id, const char* name, xt::xtensor<T, N>& result,
-     bool must_have=false)
+inline void read_nd_vector(hid_t obj_id, const char* name,
+                           xt::xtensor<T, N>& result, bool must_have=false)
 {
   if (object_exists(obj_id, name)) {
     read_dataset_as_shape(obj_id, name, result, true);
@@ -431,7 +444,8 @@ template<typename T> inline
 std::enable_if_t<std::is_scalar<std::decay_t<T>>::value>
 write_dataset(hid_t obj_id, const char* name, T buffer)
 {
-  write_dataset(obj_id, 0, nullptr, name, H5TypeMap<T>::type_id, &buffer, false);
+  write_dataset_lowlevel(obj_id, 0, nullptr, name, H5TypeMap<T>::type_id,
+                         H5S_ALL, false, &buffer);
 }
 
 inline void
@@ -444,11 +458,13 @@ template<typename T, std::size_t N> inline void
 write_dataset(hid_t obj_id, const char* name, const std::array<T, N>& buffer)
 {
   hsize_t dims[] {N};
-  write_dataset(obj_id, 1, dims, name, H5TypeMap<T>::type_id, buffer.data(), false);
+  write_dataset_lowlevel(obj_id, 1, dims, name, H5TypeMap<T>::type_id,
+                         H5S_ALL, false, buffer.data());
 }
 
 inline void
-write_dataset(hid_t obj_id, const char* name, const std::vector<std::string>& buffer)
+write_dataset(hid_t obj_id, const char* name,
+              const std::vector<std::string>& buffer)
 {
   auto n {buffer.size()};
   hsize_t dims[] {n};
@@ -477,7 +493,8 @@ template<typename T> inline void
 write_dataset(hid_t obj_id, const char* name, const std::vector<T>& buffer)
 {
   hsize_t dims[] {buffer.size()};
-  write_dataset(obj_id, 1, dims, name, H5TypeMap<T>::type_id, buffer.data(), false);
+  write_dataset_lowlevel(obj_id, 1, dims, name, H5TypeMap<T>::type_id,
+                         H5S_ALL, false, buffer.data());
 }
 
 // Template for xarray, xtensor, etc.
@@ -487,8 +504,8 @@ write_dataset(hid_t obj_id, const char* name, const xt::xcontainer<D>& arr)
   using T = typename D::value_type;
   auto s = arr.shape();
   std::vector<hsize_t> dims {s.cbegin(), s.cend()};
-  write_dataset(obj_id, dims.size(), dims.data(), name, H5TypeMap<T>::type_id,
-                arr.data(), false);
+  write_dataset_lowlevel(obj_id, dims.size(), dims.data(), name,
+                         H5TypeMap<T>::type_id, H5S_ALL, false, arr.data());
 }
 
 inline void
diff --git a/src/hdf5_interface.cpp b/src/hdf5_interface.cpp
index b1a90cb8f..bd73e05fa 100644
--- a/src/hdf5_interface.cpp
+++ b/src/hdf5_interface.cpp
@@ -49,7 +49,8 @@ get_shape(hid_t obj_id, hsize_t* dims)
   } else if (type == H5I_ATTR) {
     dspace = H5Aget_space(obj_id);
   } else {
-    throw std::runtime_error{"Expected dataset or attribute in call to get_shape."};
+    throw std::runtime_error{
+      "Expected dataset or attribute in call to get_shape."};
   }
   H5Sget_simple_extent_dims(dspace, dims, nullptr);
   H5Sclose(dspace);
@@ -74,7 +75,8 @@ std::vector<hsize_t> object_shape(hid_t obj_id)
   } else if (type == H5I_ATTR) {
     dspace = H5Aget_space(obj_id);
   } else {
-    throw std::runtime_error{"Expected dataset or attribute in call to object_shape."};
+    throw std::runtime_error{
+      "Expected dataset or attribute in call to object_shape."};
   }
   int n = H5Sget_simple_extent_ndims(dspace);
 
@@ -471,8 +473,8 @@ read_attr_string(hid_t obj_id, const char* name, size_t slen, char* buffer)
 
 
 void
-read_dataset(hid_t obj_id, const char* name, hid_t mem_type_id,
-             void* buffer, bool indep)
+read_dataset_lowlevel(hid_t obj_id, const char* name, hid_t mem_type_id,
+                      hid_t mem_space_id, bool indep, void* buffer)
 {
   hid_t dset = obj_id;
   if (name) dset = open_dataset(obj_id, name);
@@ -487,11 +489,11 @@ read_dataset(hid_t obj_id, const char* name, hid_t mem_type_id,
     H5Pset_dxpl_mpio(plist, data_xfer_mode);
 
     // Read data
-    H5Dread(dset, mem_type_id, H5S_ALL, H5S_ALL, plist, buffer);
+    H5Dread(dset, mem_type_id, mem_space_id, H5S_ALL, plist, buffer);
     H5Pclose(plist);
 #endif
   } else {
-    H5Dread(dset, mem_type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, buffer);
+    H5Dread(dset, mem_type_id, mem_space_id, H5S_ALL, H5P_DEFAULT, buffer);
   }
 
   if (name) H5Dclose(dset);
@@ -520,26 +522,28 @@ void read_dataset(hid_t dset, xt::xarray<std::complex<double>>& arr, bool indep)
 void
 read_double(hid_t obj_id, const char* name, double* buffer, bool indep)
 {
-  read_dataset(obj_id, name, H5T_NATIVE_DOUBLE, buffer, indep);
+  read_dataset_lowlevel(obj_id, name, H5T_NATIVE_DOUBLE, H5S_ALL, indep,
+                        buffer);
 }
 
 
 void
 read_int(hid_t obj_id, const char* name, int* buffer, bool indep)
 {
-  read_dataset(obj_id, name, H5T_NATIVE_INT, buffer, indep);
+  read_dataset_lowlevel(obj_id, name, H5T_NATIVE_INT, H5S_ALL, indep, buffer);
 }
 
 
 void
 read_llong(hid_t obj_id, const char* name, long long* buffer, bool indep)
 {
-  read_dataset(obj_id, name, H5T_NATIVE_LLONG, buffer, indep);
+  read_dataset_lowlevel(obj_id, name, H5T_NATIVE_LLONG, H5S_ALL, indep, buffer);
 }
 
 
 void
-read_string(hid_t obj_id, const char* name, size_t slen, char* buffer, bool indep)
+read_string(hid_t obj_id, const char* name, size_t slen, char* buffer,
+            bool indep)
 {
   // Create datatype for a string
   hid_t datatype = H5Tcopy(H5T_C_S1);
@@ -548,7 +552,7 @@ read_string(hid_t obj_id, const char* name, size_t slen, char* buffer, bool inde
   H5Tset_strpad(datatype, H5T_STR_NULLPAD);
 
   // Read data into buffer
-  read_dataset(obj_id, name, datatype, buffer, indep);
+  read_dataset_lowlevel(obj_id, name, datatype, H5S_ALL, indep, buffer);
 
   // Free resources
   H5Tclose(datatype);
@@ -556,7 +560,8 @@ read_string(hid_t obj_id, const char* name, size_t slen, char* buffer, bool inde
 
 
 void
-read_complex(hid_t obj_id, const char* name, std::complex<double>* buffer, bool indep)
+read_complex(hid_t obj_id, const char* name, std::complex<double>* buffer,
+             bool indep)
 {
   // Create compound datatype for complex numbers
   struct complex_t {
@@ -569,7 +574,7 @@ read_complex(hid_t obj_id, const char* name, std::complex<double>* buffer, bool
   H5Tinsert(complex_id, "i", HOFFSET(complex_t, im), H5T_NATIVE_DOUBLE);
 
   // Read data
-  read_dataset(obj_id, name, complex_id, buffer, indep);
+  read_dataset_lowlevel(obj_id, name, complex_id, H5S_ALL, indep, buffer);
 
   // Free resources
   H5Tclose(complex_id);
@@ -577,21 +582,22 @@ read_complex(hid_t obj_id, const char* name, std::complex<double>* buffer, bool
 
 
 void
-read_tally_results(hid_t group_id, hsize_t n_filter, hsize_t n_score, double* results)
+read_tally_results(hid_t group_id, hsize_t n_filter, hsize_t n_score,
+                   double* results)
 {
   // Create dataspace for hyperslab in memory
-  hsize_t dims[] {n_filter, n_score, 3};
-  hsize_t start[] {0, 0, 1};
-  hsize_t count[] {n_filter, n_score, 2};
-  hid_t memspace = H5Screate_simple(3, dims, nullptr);
+  constexpr int ndim = 3;
+  hsize_t dims[ndim] {n_filter, n_score, 3};
+  hsize_t start[ndim] {0, 0, 1};
+  hsize_t count[ndim] {n_filter, n_score, 2};
+  hid_t memspace = H5Screate_simple(ndim, dims, nullptr);
   H5Sselect_hyperslab(memspace, H5S_SELECT_SET, start, nullptr, count, nullptr);
 
-  // Create and write dataset
-  hid_t dset = H5Dopen(group_id, "results", H5P_DEFAULT);
-  H5Dread(dset, H5T_NATIVE_DOUBLE, memspace, H5S_ALL, H5P_DEFAULT, results);
+  // Read the dataset
+  read_dataset_lowlevel(group_id, "results", H5T_NATIVE_DOUBLE, memspace,
+                        false, results);
 
   // Free resources
-  H5Dclose(dset);
   H5Sclose(memspace);
 }
 
@@ -654,8 +660,9 @@ write_attr_string(hid_t obj_id, const char* name, const char* buffer)
 
 
 void
-write_dataset(hid_t group_id, int ndim, const hsize_t* dims, const char* name,
-              hid_t mem_type_id, const void* buffer, bool indep)
+write_dataset_lowlevel(hid_t group_id, int ndim, const hsize_t* dims,
+  const char* name, hid_t mem_type_id, hid_t mem_space_id, bool indep,
+  const void* buffer)
 {
   // If array is given, create a simple dataspace. Otherwise, create a scalar
   // datascape.
@@ -679,11 +686,11 @@ write_dataset(hid_t group_id, int ndim, const hsize_t* dims, const char* name,
     H5Pset_dxpl_mpio(plist, data_xfer_mode);
 
     // Write data
-    H5Dwrite(dset, mem_type_id, H5S_ALL, H5S_ALL, plist, buffer);
+    H5Dwrite(dset, mem_type_id, mem_space_id, H5S_ALL, plist, buffer);
     H5Pclose(plist);
 #endif
   } else {
-    H5Dwrite(dset, mem_type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, buffer);
+    H5Dwrite(dset, mem_type_id, mem_space_id, H5S_ALL, H5P_DEFAULT, buffer);
   }
 
   // Free resources
@@ -696,7 +703,8 @@ void
 write_double(hid_t group_id, int ndim, const hsize_t* dims, const char* name,
              const double* buffer, bool indep)
 {
-  write_dataset(group_id, ndim, dims, name, H5T_NATIVE_DOUBLE, buffer, indep);
+  write_dataset_lowlevel(group_id, ndim, dims, name, H5T_NATIVE_DOUBLE, H5S_ALL,
+                         indep, buffer);
 }
 
 
@@ -704,7 +712,8 @@ void
 write_int(hid_t group_id, int ndim, const hsize_t* dims, const char* name,
           const int* buffer, bool indep)
 {
-  write_dataset(group_id, ndim, dims, name, H5T_NATIVE_INT, buffer, indep);
+  write_dataset_lowlevel(group_id, ndim, dims, name, H5T_NATIVE_INT, H5S_ALL,
+                         indep, buffer);
 }
 
 
@@ -712,7 +721,8 @@ void
 write_llong(hid_t group_id, int ndim, const hsize_t* dims, const char* name,
             const long long* buffer, bool indep)
 {
-  write_dataset(group_id, ndim, dims, name, H5T_NATIVE_LLONG, buffer, indep);
+  write_dataset_lowlevel(group_id, ndim, dims, name, H5T_NATIVE_LLONG, H5S_ALL,
+                         indep, buffer);
 }
 
 
@@ -725,7 +735,8 @@ write_string(hid_t group_id, int ndim, const hsize_t* dims, size_t slen,
     hid_t datatype = H5Tcopy(H5T_C_S1);
     H5Tset_size(datatype, slen);
 
-    write_dataset(group_id, ndim, dims, name, datatype, buffer, indep);
+    write_dataset_lowlevel(group_id, ndim, dims, name, datatype, H5S_ALL, indep,
+                           buffer);
 
     // Free resources
     H5Tclose(datatype);
@@ -734,34 +745,34 @@ write_string(hid_t group_id, int ndim, const hsize_t* dims, size_t slen,
 
 
 void
-write_string(hid_t group_id, const char* name, const std::string& buffer, bool indep)
+write_string(hid_t group_id, const char* name, const std::string& buffer,
+             bool indep)
 {
-  write_string(group_id, 0, nullptr, buffer.length(), name, buffer.c_str(), indep);
+  write_string(group_id, 0, nullptr, buffer.length(), name, buffer.c_str(),
+               indep);
 }
 
 
 void
-write_tally_results(hid_t group_id, hsize_t n_filter, hsize_t n_score, const double* results)
+write_tally_results(hid_t group_id, hsize_t n_filter, hsize_t n_score,
+                    const double* results)
 {
   // Set dimensions of sum/sum_sq hyperslab to store
-  hsize_t count[] {n_filter, n_score, 2};
-  hid_t dspace = H5Screate_simple(3, count, nullptr);
+  constexpr int ndim = 3;
+  hsize_t count[ndim] {n_filter, n_score, 2};
 
   // Set dimensions of results array
-  hsize_t dims[] {n_filter, n_score, 3};
-  hsize_t start[] {0, 0, 1};
-  hid_t memspace = H5Screate_simple(3, dims, nullptr);
+  hsize_t dims[ndim] {n_filter, n_score, 3};
+  hsize_t start[ndim] {0, 0, 1};
+  hid_t memspace = H5Screate_simple(ndim, dims, nullptr);
   H5Sselect_hyperslab(memspace, H5S_SELECT_SET, start, nullptr, count, nullptr);
 
   // Create and write dataset
-  hid_t dset = H5Dcreate(group_id, "results", H5T_NATIVE_DOUBLE, dspace,
-                         H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
-  H5Dwrite(dset, H5T_NATIVE_DOUBLE, memspace, H5S_ALL, H5P_DEFAULT, results);
+  write_dataset_lowlevel(group_id, ndim, count, "results", H5T_NATIVE_DOUBLE,
+                         memspace, false, results);
 
   // Free resources
-  H5Dclose(dset);
   H5Sclose(memspace);
-  H5Sclose(dspace);
 }
 
 
diff --git a/src/state_point.cpp b/src/state_point.cpp
index 39f7a198c..60fd3826d 100644
--- a/src/state_point.cpp
+++ b/src/state_point.cpp
@@ -420,8 +420,8 @@ void load_state_point()
   if (mpi::master) {
 #endif
     // Read global tally data
-    read_dataset(file_id, "global_tallies", H5T_NATIVE_DOUBLE,
-      simulation::global_tallies.data(), false);
+    read_dataset_lowlevel(file_id, "global_tallies", H5T_NATIVE_DOUBLE,
+      H5S_ALL, false, simulation::global_tallies.data());
 
     // Check if tally results are present
     bool present;