diff --git a/include/openmc/hdf5_interface.h b/include/openmc/hdf5_interface.h
index 72eaa3936..ddb2319ae 100644
--- a/include/openmc/hdf5_interface.h
+++ b/include/openmc/hdf5_interface.h
@@ -32,10 +32,11 @@ 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_lowlevel(hid_t obj_id, const char* name, hid_t mem_type_id,
-  void* buffer, bool indep);
+  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, const void* buffer, bool indep);
+  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);
 
@@ -235,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_lowlevel(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
@@ -256,8 +258,8 @@ 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_lowlevel(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
@@ -271,8 +273,8 @@ void read_dataset(hid_t dset, std::vector<T>& vec, bool indep=false)
   vec.resize(shape[0]);
 
   // Read data into vector
-  read_dataset_lowlevel(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>
@@ -297,8 +299,8 @@ void read_dataset(hid_t dset, xt::xarray<T>& arr, bool indep=false)
   arr.resize(shape);
 
   // Read data from attribute
-  read_dataset_lowlevel(dset, nullptr, H5TypeMap<T>::type_id, arr.data(),
-                        indep);
+  read_dataset_lowlevel(dset, nullptr, H5TypeMap<T>::type_id, H5S_ALL, indep,
+                        arr.data());
 }
 
 template<>
@@ -367,8 +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_lowlevel(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());
@@ -443,7 +445,7 @@ 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_lowlevel(obj_id, 0, nullptr, name, H5TypeMap<T>::type_id,
-                         &buffer, false);
+                         H5S_ALL, false, &buffer);
 }
 
 inline void
@@ -457,7 +459,7 @@ write_dataset(hid_t obj_id, const char* name, const std::array<T, N>& buffer)
 {
   hsize_t dims[] {N};
   write_dataset_lowlevel(obj_id, 1, dims, name, H5TypeMap<T>::type_id,
-                         buffer.data(), false);
+                         H5S_ALL, false, buffer.data());
 }
 
 inline void
@@ -492,7 +494,7 @@ write_dataset(hid_t obj_id, const char* name, const std::vector<T>& buffer)
 {
   hsize_t dims[] {buffer.size()};
   write_dataset_lowlevel(obj_id, 1, dims, name, H5TypeMap<T>::type_id,
-                         buffer.data(), false);
+                         H5S_ALL, false, buffer.data());
 }
 
 // Template for xarray, xtensor, etc.
@@ -503,7 +505,7 @@ write_dataset(hid_t obj_id, const char* name, const xt::xcontainer<D>& arr)
   auto s = arr.shape();
   std::vector<hsize_t> dims {s.cbegin(), s.cend()};
   write_dataset_lowlevel(obj_id, dims.size(), dims.data(), name,
-                         H5TypeMap<T>::type_id, arr.data(), false);
+                         H5TypeMap<T>::type_id, H5S_ALL, false, arr.data());
 }
 
 inline void
diff --git a/src/hdf5_interface.cpp b/src/hdf5_interface.cpp
index d4434225a..bd73e05fa 100644
--- a/src/hdf5_interface.cpp
+++ b/src/hdf5_interface.cpp
@@ -474,7 +474,7 @@ read_attr_string(hid_t obj_id, const char* name, size_t slen, char* buffer)
 
 void
 read_dataset_lowlevel(hid_t obj_id, const char* name, hid_t mem_type_id,
-                      void* buffer, bool indep)
+                      hid_t mem_space_id, bool indep, void* buffer)
 {
   hid_t dset = obj_id;
   if (name) dset = open_dataset(obj_id, name);
@@ -489,11 +489,11 @@ read_dataset_lowlevel(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);
@@ -522,21 +522,22 @@ 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_lowlevel(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_lowlevel(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_lowlevel(obj_id, name, H5T_NATIVE_LLONG, buffer, indep);
+  read_dataset_lowlevel(obj_id, name, H5T_NATIVE_LLONG, H5S_ALL, indep, buffer);
 }
 
 
@@ -551,7 +552,7 @@ read_string(hid_t obj_id, const char* name, size_t slen, char* buffer,
   H5Tset_strpad(datatype, H5T_STR_NULLPAD);
 
   // Read data into buffer
-  read_dataset_lowlevel(obj_id, name, datatype, buffer, indep);
+  read_dataset_lowlevel(obj_id, name, datatype, H5S_ALL, indep, buffer);
 
   // Free resources
   H5Tclose(datatype);
@@ -573,7 +574,7 @@ read_complex(hid_t obj_id, const char* name, std::complex<double>* buffer,
   H5Tinsert(complex_id, "i", HOFFSET(complex_t, im), H5T_NATIVE_DOUBLE);
 
   // Read data
-  read_dataset_lowlevel(obj_id, name, complex_id, buffer, indep);
+  read_dataset_lowlevel(obj_id, name, complex_id, H5S_ALL, indep, buffer);
 
   // Free resources
   H5Tclose(complex_id);
@@ -585,18 +586,18 @@ 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);
 }
 
@@ -660,7 +661,8 @@ write_attr_string(hid_t obj_id, const char* name, const char* buffer)
 
 void
 write_dataset_lowlevel(hid_t group_id, int ndim, const hsize_t* dims,
-  const char* name, hid_t mem_type_id, const void* buffer, bool indep)
+  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.
@@ -684,11 +686,11 @@ write_dataset_lowlevel(hid_t group_id, int ndim, const hsize_t* dims,
     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
@@ -701,8 +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_lowlevel(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);
 }
 
 
@@ -710,8 +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_lowlevel(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);
 }
 
 
@@ -719,8 +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_lowlevel(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);
 }
 
 
@@ -733,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_lowlevel(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);
@@ -755,24 +758,21 @@ 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 31454ec4d..fb82bb6d2 100644
--- a/src/state_point.cpp
+++ b/src/state_point.cpp
@@ -419,7 +419,7 @@ void load_state_point()
 #endif
     // Read global tally data
     read_dataset_lowlevel(file_id, "global_tallies", H5T_NATIVE_DOUBLE,
-      simulation::global_tallies.data(), false);
+      H5S_ALL, false, simulation::global_tallies.data());
 
     // Check if tally results are present
     bool present;