diff --git a/include/openmc/capi.h b/include/openmc/capi.h
index 6c3338854c..e674dcc2ac 100644
--- a/include/openmc/capi.h
+++ b/include/openmc/capi.h
@@ -100,7 +100,7 @@ extern "C" {
   int openmc_sphharm_filter_get_cosine(int32_t index, char cosine[]);
   int openmc_sphharm_filter_set_order(int32_t index, int order);
   int openmc_sphharm_filter_set_cosine(int32_t index, const char cosine[]);
-  int openmc_statepoint_write(const char filename[], bool* write_source);
+  int openmc_statepoint_write(const char* filename, bool* write_source);
   int openmc_tally_allocate(int32_t index, const char* type);
   int openmc_tally_get_active(int32_t index, bool* active);
   int openmc_tally_get_estimator(int32_t index, int32_t* estimator);
diff --git a/include/openmc/constants.h b/include/openmc/constants.h
index 375ebc2767..956e38ca29 100644
--- a/include/openmc/constants.h
+++ b/include/openmc/constants.h
@@ -28,6 +28,7 @@ constexpr std::array<int, 3> VERSION {VERSION_MAJOR, VERSION_MINOR, VERSION_RELE
 constexpr int HDF5_VERSION[] {2, 0};
 
 // Version numbers for binary files
+constexpr std::array<int, 2> VERSION_STATEPOINT {17, 0};
 constexpr std::array<int, 2> VERSION_PARTICLE_RESTART {2, 0};
 constexpr std::array<int, 2> VERSION_TRACK {2, 0};
 constexpr std::array<int, 2> VERSION_SUMMARY {6, 0};
diff --git a/include/openmc/eigenvalue.h b/include/openmc/eigenvalue.h
index f7e93ec2b7..f89bdc8bf7 100644
--- a/include/openmc/eigenvalue.h
+++ b/include/openmc/eigenvalue.h
@@ -9,6 +9,7 @@
 #include <vector>
 
 #include "xtensor/xtensor.hpp"
+#include <hdf5.h>
 
 #include "openmc/particle.h"
 
@@ -78,6 +79,10 @@ void ufs_count_sites();
 //! Get UFS weight corresponding to particle's location
 extern "C" double ufs_get_weight(const Particle* p);
 
+//! Write data related to k-eigenvalue to statepoint
+//! \param[in] group HDF5 group
+void write_eigenvalue_hdf5(hid_t group);
+
 } // namespace openmc
 
 #endif // OPENMC_EIGENVALUE_H
diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h
index 7fcd58a06d..1ca4d57b4a 100644
--- a/include/openmc/mesh.h
+++ b/include/openmc/mesh.h
@@ -128,7 +128,7 @@ extern "C" void read_meshes(pugi::xml_node* root);
 
 //! Write mesh data to an HDF5 group
 //! \param[in] group HDF5 group
-extern "C" void meshes_to_hdf5(hid_t group);
+void meshes_to_hdf5(hid_t group);
 
 } // namespace openmc
 
diff --git a/openmc/capi/math.py b/openmc/capi/math.py
index e5a4b12c61..1baa9cf640 100644
--- a/openmc/capi/math.py
+++ b/openmc/capi/math.py
@@ -6,8 +6,8 @@ from numpy.ctypeslib import ndpointer
 from . import _dll
 
 
-_dll.t_percentile_c.restype = c_double
-_dll.t_percentile_c.argtypes = [c_double, c_int]
+_dll.t_percentile.restype = c_double
+_dll.t_percentile.argtypes = [c_double, c_int]
 
 _dll.calc_pn_c.restype = None
 _dll.calc_pn_c.argtypes = [c_int, c_double, ndpointer(c_double)]
@@ -58,7 +58,7 @@ def t_percentile(p, df):
 
     """
 
-    return _dll.t_percentile_c(p, df)
+    return _dll.t_percentile(p, df)
 
 
 def calc_pn(n, x):
@@ -261,7 +261,7 @@ def normal_variate(mean_value, std_dev):
     Parameters
     ----------
     mean_value : float
-        Mean of the Normal distribution 
+        Mean of the Normal distribution
     std_dev : float
         Standard deviation of the normal distribution
 
diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp
index 2b0f731cb2..446052fbb9 100644
--- a/src/eigenvalue.cpp
+++ b/src/eigenvalue.cpp
@@ -622,7 +622,7 @@ double ufs_get_weight(const Particle* p)
   }
 }
 
-extern "C" void write_eigenvalue_hdf5(hid_t group)
+void write_eigenvalue_hdf5(hid_t group)
 {
   write_dataset(group, "n_inactive", settings::n_inactive);
   write_dataset(group, "generations_per_batch", settings::gen_per_batch);
diff --git a/src/state_point.F90 b/src/state_point.F90
index e7d1611270..04e354142e 100644
--- a/src/state_point.F90
+++ b/src/state_point.F90
@@ -51,294 +51,172 @@ contains
 ! OPENMC_STATEPOINT_WRITE writes an HDF5 statepoint file to disk
 !===============================================================================
 
-  function openmc_statepoint_write(filename, write_source) result(err) bind(C)
-    type(C_PTR),     value                :: filename
-    logical(C_BOOL), intent(in), optional :: write_source
-    integer(C_INT) :: err
+  subroutine statepoint_write_f(file_id) bind(C)
+    integer(HID_T), value :: file_id
 
-    logical :: write_source_
     integer :: i, j
     integer :: i_xs
     integer, allocatable :: id_array(:)
-    integer(HID_T) :: file_id
     integer(HID_T) :: tallies_group, tally_group, &
                       filters_group, filter_group, derivs_group, &
-                      deriv_group, runtime_group
+                      deriv_group
     character(MAX_WORD_LEN), allocatable :: str_array(:)
-    character(C_CHAR), pointer :: string(:)
-    character(len=:, kind=C_CHAR), allocatable :: filename_
     character(MAX_WORD_LEN, kind=C_CHAR) :: temp_name
-    logical :: parallel
 
-    interface
-      subroutine meshes_to_hdf5(group) bind(C)
-        import HID_T
-        integer(HID_T), value :: group
-      end subroutine
-      subroutine write_eigenvalue_hdf5(group) bind(C)
-        import HID_T
-        integer(HID_T), value :: group
-      end subroutine
-      subroutine write_tally_results_nr(file_id) bind(C)
-        import HID_T
-        integer(HID_T), value :: file_id
-      end subroutine
-    end interface
+    ! Open tallies group
+    tallies_group = open_group(file_id, "tallies")
 
-    err = 0
-
-    ! Set the filename
-    if (c_associated(filename)) then
-      call c_f_pointer(filename, string, [MAX_FILE_LEN])
-      filename_ = to_f_string(string)
-    else
-      ! Set filename for state point
-      filename_ = trim(path_output) // 'statepoint.' // &
-           & zero_padded(current_batch, count_digits(n_max_batches))
-      filename_ = trim(filename_) // '.h5'
+    ! Write information for derivatives.
+    if (size(tally_derivs) > 0) then
+      derivs_group = create_group(tallies_group, "derivatives")
+      do i = 1, size(tally_derivs)
+        associate(deriv => tally_derivs(i))
+          deriv_group = create_group(derivs_group, "derivative " &
+               // trim(to_str(deriv % id)))
+          select case (deriv % variable)
+          case (DIFF_DENSITY)
+            call write_dataset(deriv_group, "independent variable", "density")
+            call write_dataset(deriv_group, "material", deriv % diff_material)
+          case (DIFF_NUCLIDE_DENSITY)
+            call write_dataset(deriv_group, "independent variable", &
+                 "nuclide_density")
+            call write_dataset(deriv_group, "material", deriv % diff_material)
+            call write_dataset(deriv_group, "nuclide", &
+                 nuclides(deriv % diff_nuclide) % name)
+          case (DIFF_TEMPERATURE)
+            call write_dataset(deriv_group, "independent variable", &
+                 "temperature")
+            call write_dataset(deriv_group, "material", deriv % diff_material)
+          case default
+            call fatal_error("Independent variable for derivative " &
+                 // trim(to_str(deriv % id)) // " not defined in &
+                 &state_point.F90.")
+          end select
+          call close_group(deriv_group)
+        end associate
+      end do
+      call close_group(derivs_group)
     end if
 
-    ! Determine whether or not to write the source bank
-    if (present(write_source)) then
-      write_source_ = write_source
-    else
-      write_source_ = .true.
+    ! Write number of filters
+    filters_group = create_group(tallies_group, "filters")
+    call write_attribute(filters_group, "n_filters", n_filters)
+
+    if (n_filters > 0) then
+      ! Write IDs of filters
+      allocate(id_array(n_filters))
+      do i = 1, n_filters
+        id_array(i) = filters(i) % obj % id
+      end do
+      call write_attribute(filters_group, "ids", id_array)
+      deallocate(id_array)
+
+      ! Write filter information
+      FILTER_LOOP: do i = 1, n_filters
+        filter_group = create_group(filters_group, "filter " // &
+             trim(to_str(filters(i) % obj % id)))
+        call filters(i) % obj % to_statepoint(filter_group)
+        call close_group(filter_group)
+      end do FILTER_LOOP
     end if
 
-    ! Write message
-    call write_message("Creating state point " // trim(filename_) // "...", 5)
-
-    if (master) then
-      ! Create statepoint file
-      file_id = file_open(filename_, 'w')
-
-      ! Write file type
-      call write_attribute(file_id, "filetype", "statepoint")
-
-      ! Write revision number for state point file
-      call write_attribute(file_id, "version", VERSION_STATEPOINT)
-
-      ! Write OpenMC version
-      call write_attribute(file_id, "openmc_version", VERSION)
-#ifdef GIT_SHA1
-      call write_attribute(file_id, "git_sha1", GIT_SHA1)
-#endif
-
-      ! Write current date and time
-      call write_attribute(file_id, "date_and_time", time_stamp())
-
-      ! Write path to input
-      call write_attribute(file_id, "path", path_input)
-
-      ! Write out random number seed
-      call write_dataset(file_id, "seed", openmc_get_seed())
-
-      ! Write run information
-      if (run_CE) then
-        call write_dataset(file_id, "energy_mode", "continuous-energy")
-      else
-        call write_dataset(file_id, "energy_mode", "multi-group")
-      end if
-      select case(run_mode)
-      case (MODE_FIXEDSOURCE)
-        call write_dataset(file_id, "run_mode", "fixed source")
-      case (MODE_EIGENVALUE)
-        call write_dataset(file_id, "run_mode", "eigenvalue")
-      end select
-      if (photon_transport) then
-        call write_attribute(file_id, "photon_transport", 1)
-      else
-        call write_attribute(file_id, "photon_transport", 0)
-      end if
-      call write_dataset(file_id, "n_particles", n_particles)
-      call write_dataset(file_id, "n_batches", n_batches)
-
-      ! Write out current batch number
-      call write_dataset(file_id, "current_batch", current_batch)
-
-      ! Indicate whether source bank is stored in statepoint
-      if (write_source_) then
-        call write_attribute(file_id, "source_present", 1)
-      else
-        call write_attribute(file_id, "source_present", 0)
-      end if
-
-      ! Write out information for eigenvalue run
-      if (run_mode == MODE_EIGENVALUE) then
-        call write_eigenvalue_hdf5(file_id)
-      end if
-
-      tallies_group = create_group(file_id, "tallies")
-
-      ! Write meshes
-      call meshes_to_hdf5(tallies_group)
-
-      ! Write information for derivatives.
-      if (size(tally_derivs) > 0) then
-        derivs_group = create_group(tallies_group, "derivatives")
-        do i = 1, size(tally_derivs)
-          associate(deriv => tally_derivs(i))
-            deriv_group = create_group(derivs_group, "derivative " &
-                 // trim(to_str(deriv % id)))
-            select case (deriv % variable)
-            case (DIFF_DENSITY)
-              call write_dataset(deriv_group, "independent variable", "density")
-              call write_dataset(deriv_group, "material", deriv % diff_material)
-            case (DIFF_NUCLIDE_DENSITY)
-              call write_dataset(deriv_group, "independent variable", &
-                   "nuclide_density")
-              call write_dataset(deriv_group, "material", deriv % diff_material)
-              call write_dataset(deriv_group, "nuclide", &
-                   nuclides(deriv % diff_nuclide) % name)
-            case (DIFF_TEMPERATURE)
-              call write_dataset(deriv_group, "independent variable", &
-                   "temperature")
-              call write_dataset(deriv_group, "material", deriv % diff_material)
-            case default
-              call fatal_error("Independent variable for derivative " &
-                   // trim(to_str(deriv % id)) // " not defined in &
-                   &state_point.F90.")
-            end select
-            call close_group(deriv_group)
-          end associate
-        end do
-        call close_group(derivs_group)
-      end if
-
-      ! Write number of filters
-      filters_group = create_group(tallies_group, "filters")
-      call write_attribute(filters_group, "n_filters", n_filters)
-
-      if (n_filters > 0) then
-        ! Write IDs of filters
-        allocate(id_array(n_filters))
-        do i = 1, n_filters
-          id_array(i) = filters(i) % obj % id
-        end do
-        call write_attribute(filters_group, "ids", id_array)
-        deallocate(id_array)
-
-        ! Write filter information
-        FILTER_LOOP: do i = 1, n_filters
-          filter_group = create_group(filters_group, "filter " // &
-               trim(to_str(filters(i) % obj % id)))
-          call filters(i) % obj % to_statepoint(filter_group)
-          call close_group(filter_group)
-        end do FILTER_LOOP
-      end if
-
-      call close_group(filters_group)
+    call close_group(filters_group)
 
       ! Write number of tallies
-      call write_attribute(tallies_group, "n_tallies", n_tallies)
+    call write_attribute(tallies_group, "n_tallies", n_tallies)
 
-      if (n_tallies > 0) then
-        ! Write array of tally IDs
-        allocate(id_array(n_tallies))
-        do i = 1, n_tallies
-          id_array(i) = tallies(i) % obj % id
-        end do
-        call write_attribute(tallies_group, "ids", id_array)
-        deallocate(id_array)
+    if (n_tallies > 0) then
+      ! Write array of tally IDs
+      allocate(id_array(n_tallies))
+      do i = 1, n_tallies
+        id_array(i) = tallies(i) % obj % id
+      end do
+      call write_attribute(tallies_group, "ids", id_array)
+      deallocate(id_array)
 
-        ! Write all tally information except results
-        TALLY_METADATA: do i = 1, n_tallies
+      ! Write all tally information except results
+      TALLY_METADATA: do i = 1, n_tallies
 
-          ! Get pointer to tally
-          associate (tally => tallies(i) % obj)
-          tally_group = create_group(tallies_group, "tally " // &
-               trim(to_str(tally % id)))
+        ! Get pointer to tally
+        associate (tally => tallies(i) % obj)
+        tally_group = create_group(tallies_group, "tally " // &
+             trim(to_str(tally % id)))
 
-          ! Write the name for this tally
-          call write_dataset(tally_group, "name", tally % name)
+        ! Write the name for this tally
+        call write_dataset(tally_group, "name", tally % name)
 
-          select case(tally % estimator)
-          case (ESTIMATOR_ANALOG)
-            call write_dataset(tally_group, "estimator", "analog")
-          case (ESTIMATOR_TRACKLENGTH)
-            call write_dataset(tally_group, "estimator", "tracklength")
-          case (ESTIMATOR_COLLISION)
-            call write_dataset(tally_group, "estimator", "collision")
-          end select
-          call write_dataset(tally_group, "n_realizations", &
-               tally % n_realizations)
+        select case(tally % estimator)
+        case (ESTIMATOR_ANALOG)
+          call write_dataset(tally_group, "estimator", "analog")
+        case (ESTIMATOR_TRACKLENGTH)
+          call write_dataset(tally_group, "estimator", "tracklength")
+        case (ESTIMATOR_COLLISION)
+          call write_dataset(tally_group, "estimator", "collision")
+        end select
+        call write_dataset(tally_group, "n_realizations", &
+             tally % n_realizations)
 
-          call write_dataset(tally_group, "n_filters", size(tally % filter))
-          if (size(tally % filter) > 0) then
-            ! Write IDs of filters
-            allocate(id_array(size(tally % filter)))
-            do j = 1, size(tally % filter)
-              id_array(j) = filters(tally % filter(j)) % obj % id
-            end do
-            call write_dataset(tally_group, "filters", id_array)
-            deallocate(id_array)
-          end if
+        call write_dataset(tally_group, "n_filters", size(tally % filter))
+        if (size(tally % filter) > 0) then
+          ! Write IDs of filters
+          allocate(id_array(size(tally % filter)))
+          do j = 1, size(tally % filter)
+            id_array(j) = filters(tally % filter(j)) % obj % id
+          end do
+          call write_dataset(tally_group, "filters", id_array)
+          deallocate(id_array)
+        end if
 
-          ! Set up nuclide bin array and then write
-          allocate(str_array(tally % n_nuclide_bins))
-          NUCLIDE_LOOP: do j = 1, tally % n_nuclide_bins
-            if (tally % nuclide_bins(j) > 0) then
-              if (run_CE) then
-                i_xs = index(nuclides(tally % nuclide_bins(j)) % name, '.')
-                if (i_xs > 0) then
-                  str_array(j) = nuclides(tally % nuclide_bins(j)) % name(1 : i_xs-1)
-                else
-                  str_array(j) = nuclides(tally % nuclide_bins(j)) % name
-                end if
+        ! Set up nuclide bin array and then write
+        allocate(str_array(tally % n_nuclide_bins))
+        NUCLIDE_LOOP: do j = 1, tally % n_nuclide_bins
+          if (tally % nuclide_bins(j) > 0) then
+            if (run_CE) then
+              i_xs = index(nuclides(tally % nuclide_bins(j)) % name, '.')
+              if (i_xs > 0) then
+                str_array(j) = nuclides(tally % nuclide_bins(j)) % name(1 : i_xs-1)
               else
-                call get_name_c(tally % nuclide_bins(j), len(temp_name), &
-                                temp_name)
-                i_xs = index(temp_name, '.')
-                if (i_xs > 0) then
-                  str_array(j) = trim(temp_name(1 : i_xs-1))
-                else
-                  str_array(j) = trim(temp_name)
-                end if
+                str_array(j) = nuclides(tally % nuclide_bins(j)) % name
               end if
             else
-              str_array(j) = 'total'
+              call get_name_c(tally % nuclide_bins(j), len(temp_name), &
+                              temp_name)
+              i_xs = index(temp_name, '.')
+              if (i_xs > 0) then
+                str_array(j) = trim(temp_name(1 : i_xs-1))
+              else
+                str_array(j) = trim(temp_name)
+              end if
             end if
-          end do NUCLIDE_LOOP
-          call write_dataset(tally_group, "nuclides", str_array)
-          deallocate(str_array)
-
-          ! Write derivative information.
-          if (tally % deriv /= NONE) then
-            call write_dataset(tally_group, "derivative", &
-                 tally_derivs(tally % deriv) % id)
+          else
+            str_array(j) = 'total'
           end if
+        end do NUCLIDE_LOOP
+        call write_dataset(tally_group, "nuclides", str_array)
+        deallocate(str_array)
 
-          ! Write scores.
-          call write_dataset(tally_group, "n_score_bins", tally % n_score_bins)
-          allocate(str_array(size(tally % score_bins)))
-          do j = 1, size(tally % score_bins)
-            str_array(j) = reaction_name(tally % score_bins(j))
-          end do
-          call write_dataset(tally_group, "score_bins", str_array)
+        ! Write derivative information.
+        if (tally % deriv /= NONE) then
+          call write_dataset(tally_group, "derivative", &
+               tally_derivs(tally % deriv) % id)
+        end if
 
-          deallocate(str_array)
+        ! Write scores.
+        call write_dataset(tally_group, "n_score_bins", tally % n_score_bins)
+        allocate(str_array(size(tally % score_bins)))
+        do j = 1, size(tally % score_bins)
+          str_array(j) = reaction_name(tally % score_bins(j))
+        end do
+        call write_dataset(tally_group, "score_bins", str_array)
 
-          call close_group(tally_group)
-          end associate
-        end do TALLY_METADATA
+        deallocate(str_array)
 
-      end if
-
-      call close_group(tallies_group)
+        call close_group(tally_group)
+        end associate
+      end do TALLY_METADATA
     end if
 
-    ! Check for the no-tally-reduction method
-    if (.not. reduce_tallies) then
-      ! If using the no-tally-reduction method, we need to collect tally
-      ! results before writing them to the state point file.
-
-      call write_tally_results_nr(file_id)
-
-    elseif (master) then
-
-      ! Write number of global realizations
-      call write_dataset(file_id, "n_realizations", n_realizations)
-
+    if (reduce_tallies) then
       ! Write global tallies
       call write_dataset(file_id, "global_tallies", global_tallies)
 
@@ -347,8 +225,6 @@ contains
         ! Indicate that tallies are on
         call write_attribute(file_id, "tallies_present", 1)
 
-        tallies_group = open_group(file_id, "tallies")
-
         ! Write all tally results
         TALLY_RESULTS: do i = 1, n_tallies
           associate (tally => tallies(i) % obj)
@@ -359,61 +235,14 @@ contains
             call close_group(tally_group)
           end associate
         end do TALLY_RESULTS
-
-        call close_group(tallies_group)
       else
         ! Indicate tallies are off
         call write_attribute(file_id, "tallies_present", 0)
       end if
-
-
-      ! Write out the runtime metrics.
-      runtime_group = create_group(file_id, "runtime")
-      call write_dataset(runtime_group, "total initialization", &
-           time_initialize_elapsed())
-      call write_dataset(runtime_group, "reading cross sections", &
-           time_read_xs_elapsed())
-      call write_dataset(runtime_group, "simulation", &
-           time_inactive_elapsed() + time_active_elapsed())
-      call write_dataset(runtime_group, "transport", &
-           time_transport_elapsed())
-      if (run_mode == MODE_EIGENVALUE) then
-        call write_dataset(runtime_group, "inactive batches", &
-             time_inactive_elapsed())
-      end if
-      call write_dataset(runtime_group, "active batches", &
-           time_active_elapsed())
-      if (run_mode == MODE_EIGENVALUE) then
-        call write_dataset(runtime_group, "synchronizing fission bank", &
-             time_bank_elapsed())
-        call write_dataset(runtime_group, "sampling source sites", &
-             time_bank_sample_elapsed())
-        call write_dataset(runtime_group, "SEND-RECV source sites", &
-             time_bank_sendrecv_elapsed())
-      end if
-      call write_dataset(runtime_group, "accumulating tallies", &
-           time_tallies_elapsed())
-      call write_dataset(runtime_group, "total", time_total_elapsed())
-      call close_group(runtime_group)
-
-      call file_close(file_id)
     end if
 
-#ifdef PHDF5
-    parallel = .true.
-#else
-    parallel = .false.
-#endif
-
-    ! Write the source bank if desired
-    if (write_source_) then
-      if (master .or. parallel) then
-        file_id = file_open(filename_, 'a', parallel=.true.)
-      end if
-      call write_source_bank(file_id)
-      if (master .or. parallel) call file_close(file_id)
-    end if
-  end function openmc_statepoint_write
+    call close_group(tallies_group)
+  end subroutine
 
 !===============================================================================
 ! LOAD_STATE_POINT
diff --git a/src/state_point.cpp b/src/state_point.cpp
index a808e3ef59..3668b7b7c4 100644
--- a/src/state_point.cpp
+++ b/src/state_point.cpp
@@ -14,13 +14,152 @@
 #include "openmc/eigenvalue.h"
 #include "openmc/error.h"
 #include "openmc/hdf5_interface.h"
+#include "openmc/mesh.h"
 #include "openmc/message_passing.h"
+#include "openmc/output.h"
 #include "openmc/settings.h"
 #include "openmc/simulation.h"
+#include "openmc/tallies/filter.h"
 #include "openmc/tallies/tally.h"
+#include "openmc/timer.h"
 
 namespace openmc {
 
+extern "C" void statepoint_write_f(hid_t file_id);
+
+extern "C" int
+openmc_statepoint_write(const char* filename, bool* write_source)
+{
+  // Set the filename
+  std::string filename_;
+  if (filename) {
+    filename_ = filename;
+  } else {
+    // Determine width for zero padding
+    int w = std::to_string(settings::n_max_batches).size();
+
+    // Set filename for state point
+    std::stringstream ss;
+    ss << settings::path_output << "statepoint." << std::setfill('0')
+      << std::setw(w) << simulation::current_batch << ".h5";
+    filename_ = ss.str();
+  }
+
+  // Determine whether or not to write the source bank
+  bool write_source_ = write_source ? *write_source : true;
+
+  // Write message
+  write_message("Creating state point " + filename_ + "...", 5);
+
+  hid_t file_id;
+  if (mpi::master) {
+    // Create statepoint file
+    file_id = file_open(filename_, 'w');
+
+    // Write file type
+    write_attribute(file_id, "filetype", "statepoint");
+
+    // Write revision number for state point file
+    write_attribute(file_id, "version", VERSION_STATEPOINT);
+
+    // Write OpenMC version
+    write_attribute(file_id, "openmc_version", VERSION);
+#ifdef GIT_SHA1
+    write_attribute(file_id, "git_sha1", GIT_SHA1);
+#endif
+
+    // Write current date and time
+    write_attribute(file_id, "date_and_time", time_stamp());
+
+    // Write path to input
+    write_attribute(file_id, "path", settings::path_input);
+
+    // Write out random number seed
+    write_dataset(file_id, "seed", openmc_get_seed());
+
+    // Write run information
+    write_dataset(file_id, "energy_mode", settings::run_CE ?
+      "continuous-energy" : "multi-group");
+    switch (settings::run_mode) {
+    case RUN_MODE_FIXEDSOURCE:
+      write_dataset(file_id, "run_mode", "fixed source");
+      break;
+    case RUN_MODE_EIGENVALUE:
+      write_dataset(file_id, "run_mode", "eigenvalue");
+      break;
+    }
+    write_attribute(file_id, "photon_transport", settings::photon_transport);
+    write_dataset(file_id, "n_particles", settings::n_particles);
+    write_dataset(file_id, "n_batches", settings::n_batches);
+
+    // Write out current batch number
+    write_dataset(file_id, "current_batch", simulation::current_batch);
+
+    // Indicate whether source bank is stored in statepoint
+    write_attribute(file_id, "source_present", write_source_);
+
+    // Write out information for eigenvalue run
+    if (settings::run_mode == RUN_MODE_EIGENVALUE) write_eigenvalue_hdf5(file_id);
+
+    // Write tally-related data
+    hid_t tallies_group = create_group(file_id, "tallies");
+    meshes_to_hdf5(tallies_group);
+    statepoint_write_f(file_id);
+    close_group(tallies_group);
+  }
+
+  // Check for the no-tally-reduction method
+  if (!settings::reduce_tallies) {
+    // If using the no-tally-reduction method, we need to collect tally
+    // results before writing them to the state point file.
+
+    write_tally_results_nr(file_id);
+
+  } else if (mpi::master) {
+
+    // Write number of global realizations
+    write_dataset(file_id, "n_realizations", n_realizations);
+
+    // Write out the runtime metrics.
+    using namespace simulation;
+    hid_t runtime_group = create_group(file_id, "runtime");
+    write_dataset(runtime_group, "total initialization",  time_initialize.elapsed());
+    write_dataset(runtime_group, "reading cross sections", time_read_xs.elapsed());
+    write_dataset(runtime_group, "simulation", time_inactive.elapsed()
+      + time_active.elapsed());
+    write_dataset(runtime_group, "transport", time_transport.elapsed());
+    if (settings::run_mode == RUN_MODE_EIGENVALUE) {
+      write_dataset(runtime_group, "inactive batches", time_inactive.elapsed());
+    }
+    write_dataset(runtime_group, "active batches", time_active.elapsed());
+    if (settings::run_mode == RUN_MODE_EIGENVALUE) {
+      write_dataset(runtime_group, "synchronizing fission bank", time_bank.elapsed());
+      write_dataset(runtime_group, "sampling source sites", time_bank_sample.elapsed());
+      write_dataset(runtime_group, "SEND-RECV source sites", time_bank_sendrecv.elapsed());
+    }
+    write_dataset(runtime_group, "accumulating tallies", time_tallies.elapsed());
+    write_dataset(runtime_group, "total", time_total.elapsed());
+    close_group(runtime_group);
+
+    file_close(file_id);
+  }
+
+#ifdef PHDF5
+  bool parallel = true;
+#else
+  bool parallel = false;
+#endif
+
+  // Write the source bank if desired
+  if (write_source_) {
+    if (mpi::master || parallel) file_id = file_open(filename_, 'a', true);
+    write_source_bank(file_id);
+    if (mpi::master || parallel) file_close(file_id);
+  }
+
+  return 0;
+}
+
 hid_t h5banktype() {
   // Create type for array of 3 reals
   hsize_t dims[] {3};