diff --git a/include/openmc/capi.h b/include/openmc/capi.h
index 05512b379..5441c1026 100644
--- a/include/openmc/capi.h
+++ b/include/openmc/capi.h
@@ -95,7 +95,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[]);
+  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/container_util.h b/include/openmc/container_util.h
new file mode 100644
index 000000000..1e65a18df
--- /dev/null
+++ b/include/openmc/container_util.h
@@ -0,0 +1,17 @@
+#ifndef OPENMC_CONTAINER_UTIL_H
+#define OPENMC_CONTAINER_UTIL_H
+
+#include <algorithm> // for find
+#include <iterator> // for begin, end
+
+namespace openmc {
+
+template<class C, class T>
+inline bool contains(const C& v, const T& x)
+{
+  return std::end(v) != std::find(std::begin(v), std::end(v), x);
+}
+
+}
+
+#endif // OPENMC_CONTAINER_UTIL_H
diff --git a/include/openmc/settings.h b/include/openmc/settings.h
index 2a3111fe6..32eac43fb 100644
--- a/include/openmc/settings.h
+++ b/include/openmc/settings.h
@@ -7,6 +7,7 @@
 #include <array>
 #include <cstdint>
 #include <string>
+#include <unordered_set>
 #include <vector>
 
 #include "pugixml.hpp"
@@ -76,6 +77,8 @@ extern "C" int res_scat_method;          //!< resonance upscattering method
 extern "C" double res_scat_energy_min;   //!< Min energy in [eV] for res. upscattering
 extern "C" double res_scat_energy_max;   //!< Max energy in [eV] for res. upscattering
 extern "C" int run_mode;                 //!< Run mode (eigenvalue, fixed src, etc.)
+extern std::unordered_set<int> sourcepoint_batch; //!< Batches when source should be written
+extern std::unordered_set<int> statepoint_batch; //!< Batches when state should be written
 extern "C" int temperature_method;       //!< method for choosing temperatures
 extern "C" double temperature_tolerance; //!< Tolerance in [K] on choosing temperatures
 extern "C" double temperature_default;   //!< Default T in [K]
diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h
index 71b2e7e74..02e434aba 100644
--- a/include/openmc/simulation.h
+++ b/include/openmc/simulation.h
@@ -66,6 +66,13 @@ extern "C" void initialize_generation();
 
 extern "C" void initialize_history(Particle* p, int64_t index_source);
 
+//! Finalize a batch
+//!
+//! Handles synchronization and accumulation of tallies, calculation of Shannon
+//! entropy, getting single-batch estimate of keff, and turning on tallies when
+//! appropriate
+extern "C" void finalize_batch();
+
 //! Finalize a fission generation
 extern "C" void finalize_generation();
 
@@ -74,7 +81,6 @@ extern "C" int overall_generation();
 
 #ifdef OPENMC_MPI
 extern "C" void broadcast_results();
-extern "C" void broadcast_triggers();
 #endif
 
 } // namespace openmc
diff --git a/include/openmc/state_point.h b/include/openmc/state_point.h
index 8f98fccd8..6b019e7b8 100644
--- a/include/openmc/state_point.h
+++ b/include/openmc/state_point.h
@@ -9,6 +9,7 @@
 
 namespace openmc {
 
+void write_source_point(const char* filename);
 extern "C" void write_source_bank(hid_t group_id, Bank* source_bank);
 extern "C" void read_source_bank(hid_t group_id, Bank* source_bank);
 extern "C" void write_tally_results_nr(hid_t file_id);
diff --git a/include/openmc/timer.h b/include/openmc/timer.h
index 3d632019b..699e0fcbe 100644
--- a/include/openmc/timer.h
+++ b/include/openmc/timer.h
@@ -45,6 +45,7 @@ extern Timer time_bank_sendrecv;
 extern Timer time_finalize;
 extern Timer time_inactive;
 extern Timer time_initialize;
+extern Timer time_tallies;
 extern Timer time_total;
 
 } // namespace openmc
diff --git a/src/api.F90 b/src/api.F90
index 869660ae8..fa8b81475 100644
--- a/src/api.F90
+++ b/src/api.F90
@@ -132,8 +132,6 @@ contains
     n_batch_interval = 1
     n_lost_particles = 0
     n_particles = -1
-    n_source_points = 0
-    n_state_points = 0
     n_tallies = 0
     output_summary = .true.
     output_tallies = .true.
diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90
index 418c2f2e1..681d08c5c 100644
--- a/src/cmfd_execute.F90
+++ b/src/cmfd_execute.F90
@@ -29,7 +29,7 @@ contains
 ! EXECUTE_CMFD runs the CMFD calculation
 !==============================================================================
 
-  subroutine execute_cmfd()
+  subroutine execute_cmfd() bind(C)
 
     use cmfd_data,              only: set_up_cmfd
     use cmfd_solver,            only: cmfd_solver_execute
diff --git a/src/cmfd_header.F90 b/src/cmfd_header.F90
index 036e193b0..f4b14b6f2 100644
--- a/src/cmfd_header.F90
+++ b/src/cmfd_header.F90
@@ -143,7 +143,7 @@ module cmfd_header
   logical, public :: cmfd_run_adjoint = .false.
 
   ! CMFD run logicals
-  logical, public :: cmfd_on = .false.
+  logical(C_BOOL), public, bind(C) :: cmfd_on = .false.
 
   ! CMFD display info
   character(len=25), public :: cmfd_display = 'balance'
diff --git a/src/input_xml.F90 b/src/input_xml.F90
index 2a8ded733..cb51672da 100644
--- a/src/input_xml.F90
+++ b/src/input_xml.F90
@@ -210,9 +210,7 @@ contains
 
     integer :: i
     integer :: n
-    integer, allocatable :: temp_int_array(:)
     type(XMLNode) :: root
-    type(XMLNode) :: node_sp
     type(XMLNode) :: node_res_scat
     type(XMLNode) :: node_vol
     type(XMLNode), allocatable :: node_vol_list(:)
@@ -220,89 +218,6 @@ contains
     ! Get proper XMLNode type given pointer
     root % ptr = root_ptr
 
-    ! Check if the user has specified to write state points
-    if (check_for_node(root, "state_point")) then
-
-      ! Get pointer to state_point node
-      node_sp = root % child("state_point")
-
-      ! Determine number of batches at which to store state points
-      if (check_for_node(node_sp, "batches")) then
-        n_state_points = node_word_count(node_sp, "batches")
-      else
-        n_state_points = 0
-      end if
-
-      if (n_state_points > 0) then
-        ! User gave specific batches to write state points
-        allocate(temp_int_array(n_state_points))
-        call get_node_array(node_sp, "batches", temp_int_array)
-        do i = 1, n_state_points
-          call statepoint_batch % add(temp_int_array(i))
-        end do
-        deallocate(temp_int_array)
-      else
-        ! If neither were specified, write state point at last batch
-        n_state_points = 1
-        call statepoint_batch % add(n_batches)
-      end if
-    else
-      ! If no <state_point> tag was present, by default write state point at
-      ! last batch only
-      n_state_points = 1
-      call statepoint_batch % add(n_batches)
-    end if
-
-    ! Check if the user has specified to write source points
-    if (check_for_node(root, "source_point")) then
-
-      ! Get pointer to source_point node
-      node_sp = root % child("source_point")
-
-      ! Determine number of batches at which to store source points
-      if (check_for_node(node_sp, "batches")) then
-        n_source_points = node_word_count(node_sp, "batches")
-      else
-        n_source_points = 0
-      end if
-
-      if (n_source_points > 0) then
-        ! User gave specific batches to write source points
-        allocate(temp_int_array(n_source_points))
-        call get_node_array(node_sp, "batches", temp_int_array)
-        do i = 1, n_source_points
-          call sourcepoint_batch % add(temp_int_array(i))
-        end do
-        deallocate(temp_int_array)
-      else
-        ! If neither were specified, write source points with state points
-        n_source_points = n_state_points
-        do i = 1, n_state_points
-          call sourcepoint_batch % add(statepoint_batch % get_item(i))
-        end do
-      end if
-    else
-      ! If no <source_point> tag was present, by default we keep source bank in
-      ! statepoint file and write it out at statepoints intervals
-      n_source_points = n_state_points
-      do i = 1, n_state_points
-        call sourcepoint_batch % add(statepoint_batch % get_item(i))
-      end do
-    end if
-
-    ! If source is not seperate and is to be written out in the statepoint file,
-    ! make sure that the sourcepoint batch numbers are contained in the
-    ! statepoint list
-    if (.not. source_separate) then
-      do i = 1, n_source_points
-        if (.not. statepoint_batch % contains(sourcepoint_batch % &
-             get_item(i))) then
-          call fatal_error('Sourcepoint batches are not a subset&
-               & of statepoint batches.')
-        end if
-      end do
-    end if
-
     ! Resonance scattering parameters
     if (check_for_node(root, "resonance_scattering")) then
       node_res_scat = root % child("resonance_scattering")
diff --git a/src/output.F90 b/src/output.F90
index 3b3bd3c4d..8b4e2507b 100644
--- a/src/output.F90
+++ b/src/output.F90
@@ -360,7 +360,7 @@ contains
 ! multiplication factor as well as the average value if we're in active batches
 !===============================================================================
 
-  subroutine print_batch_keff()
+  subroutine print_batch_keff() bind(C)
 
     integer :: i  ! overall generation
     integer :: n  ! number of active generations
diff --git a/src/settings.F90 b/src/settings.F90
index d0353342b..a64d60b70 100644
--- a/src/settings.F90
+++ b/src/settings.F90
@@ -109,14 +109,6 @@ module settings
   ! Whether create fission neutrons or not. Only applied for MODE_FIXEDSOURCE
   logical(C_BOOL), bind(C) :: create_fission_neutrons
 
-  ! Information about state points to be written
-  integer :: n_state_points = 0
-  type(SetInt) :: statepoint_batch
-
-  ! Information about source points to be written
-  integer :: n_source_points = 0
-  type(SetInt) :: sourcepoint_batch
-
   character(MAX_FILE_LEN) :: path_input               ! Path to input file
   character(MAX_FILE_LEN) :: path_cross_sections = '' ! Path to cross_sections.xml
   character(MAX_FILE_LEN) :: path_multipole           ! Path to wmp library
@@ -149,10 +141,14 @@ contains
 !===============================================================================
 
   subroutine free_memory_settings()
+    interface
+      subroutine free_memory_settings_c() bind(C)
+      end subroutine
+    end interface
+
     if (allocated(res_scat_nuclides)) deallocate(res_scat_nuclides)
 
-    call statepoint_batch % clear()
-    call sourcepoint_batch % clear()
+    call free_memory_settings_c()
   end subroutine free_memory_settings
 
 end module settings
diff --git a/src/settings.cpp b/src/settings.cpp
index 3442c9a5a..ee1f21dd5 100644
--- a/src/settings.cpp
+++ b/src/settings.cpp
@@ -9,6 +9,7 @@
 
 #include "openmc/capi.h"
 #include "openmc/constants.h"
+#include "openmc/container_util.h"
 #include "openmc/distribution.h"
 #include "openmc/distribution_multi.h"
 #include "openmc/distribution_spatial.h"
@@ -87,6 +88,8 @@ int res_scat_method {RES_SCAT_ARES};
 double res_scat_energy_min {0.01};
 double res_scat_energy_max {1000.0};
 int run_mode {-1};
+std::unordered_set<int> sourcepoint_batch;
+std::unordered_set<int> statepoint_batch;
 int temperature_method {TEMPERATURE_NEAREST};
 double temperature_tolerance {10.0};
 double temperature_default {293.6};
@@ -507,7 +510,7 @@ void read_settings_xml()
     // Reshape into track_identifiers
     int n_tracks = temp.size() / 3;
     for (int i = 0; i < n_tracks; ++i) {
-      settings::track_identifiers.push_back({temp[3*i], temp[3*i + 1],
+      track_identifiers.push_back({temp[3*i], temp[3*i + 1],
         temp[3*i + 2]});
     }
   }
@@ -548,7 +551,7 @@ void read_settings_xml()
       // If the user did not specify how many mesh cells are to be used in
       // each direction, we automatically determine an appropriate number of
       // cells
-      int n = std::ceil(std::pow(settings::n_particles / 20.0, 1.0/3.0));
+      int n = std::ceil(std::pow(n_particles / 20.0, 1.0/3.0));
       m.shape_ = {n, n, n};
       m.n_dimension_ = 3;
 
@@ -557,7 +560,7 @@ void read_settings_xml()
     }
 
     // Turn on Shannon entropy calculation
-    settings::entropy_on = true;
+    entropy_on = true;
   }
 
   // Uniform fission source weighting mesh
@@ -590,18 +593,48 @@ void read_settings_xml()
 
   if (index_ufs_mesh >= 0) {
     // Turn on uniform fission source weighting
-    settings::ufs_on = true;
+    ufs_on = true;
   }
 
-  // TODO: Read <state_point>
+  // Check if the user has specified to write state points
+  if (check_for_node(root, "state_point")) {
 
+    // Get pointer to state_point node
+    auto node_sp = root.child("state_point");
+
+    // Determine number of batches at which to store state points
+    if (check_for_node(node_sp, "batches")) {
+      // User gave specific batches to write state points
+      auto temp = get_node_array<int>(node_sp, "batches");
+      for (const auto& b : temp) {
+        statepoint_batch.insert(b);
+      }
+    } else {
+      // If neither were specified, write state point at last batch
+      statepoint_batch.insert(n_batches);
+    }
+  } else {
+    // If no <state_point> tag was present, by default write state point at
+    // last batch only
+    statepoint_batch.insert(n_batches);
+  }
 
   // Check if the user has specified to write source points
   if (check_for_node(root, "source_point")) {
     // Get source_point node
     xml_node node_sp = root.child("source_point");
 
-    // TODO: Read source point batches
+    // Determine batches at which to store source points
+    if (check_for_node(node_sp, "batches")) {
+      // User gave specific batches to write source points
+      auto temp = get_node_array<int>(node_sp, "batches");
+      for (const auto& b : temp) {
+        sourcepoint_batch.insert(b);
+      }
+    } else {
+      // If neither were specified, write source points with state points
+      sourcepoint_batch = statepoint_batch;
+    }
 
     // Check if the user has specified to write binary source file
     if (check_for_node(node_sp, "separate")) {
@@ -618,10 +651,19 @@ void read_settings_xml()
     // If no <source_point> tag was present, by default we keep source bank in
     // statepoint file and write it out at statepoints intervals
     source_separate = false;
-    // TODO: add defaults
+    sourcepoint_batch = statepoint_batch;
   }
 
-  // TODO: Check source points are subset
+  // If source is not seperate and is to be written out in the statepoint file,
+  // make sure that the sourcepoint batch numbers are contained in the
+  // statepoint list
+  if (!source_separate) {
+    for (const auto& b : sourcepoint_batch) {
+      if (!contains(statepoint_batch, b)) {
+        fatal_error("Sourcepoint batches are not a subset of statepoint batches.");
+      }
+    }
+  }
 
   // Check if the user has specified to not reduce tallies at the end of every
   // batch
@@ -787,6 +829,11 @@ extern "C" {
   const char* openmc_path_particle_restart() {
     return settings::path_particle_restart.c_str();
   }
+
+  void free_memory_settings_c() {
+    settings::statepoint_batch.clear();
+    settings::sourcepoint_batch.clear();
+  }
 }
 
 } // namespace openmc
diff --git a/src/simulation.F90 b/src/simulation.F90
index 18fadc9a5..e8f56459c 100644
--- a/src/simulation.F90
+++ b/src/simulation.F90
@@ -26,7 +26,7 @@ module simulation
   use random_lcg,      only: set_particle_seed
   use settings
   use simulation_header
-  use state_point,     only: openmc_statepoint_write, write_source_point, load_state_point
+  use state_point,     only: openmc_statepoint_write, load_state_point
   use string,          only: to_str
   use tally,           only: accumulate_tallies, setup_active_tallies, &
                              init_tally_routines
@@ -55,6 +55,9 @@ module simulation
     subroutine finalize_generation() bind(C)
     end subroutine finalize_generation
 
+    subroutine finalize_batch() bind(C)
+    end subroutine
+
     subroutine initialize_history(p, index_source) bind(C)
       import Particle, C_INT64_T
       type(Particle), intent(inout) :: p
@@ -141,73 +144,6 @@ contains
 
   end function openmc_next_batch
 
-!===============================================================================
-! FINALIZE_BATCH handles synchronization and accumulation of tallies,
-! calculation of Shannon entropy, getting single-batch estimate of keff, and
-! turning on tallies when appropriate
-!===============================================================================
-
-  subroutine finalize_batch()
-
-    integer(C_INT) :: err
-    character(MAX_FILE_LEN) :: filename
-
-    interface
-      subroutine broadcast_triggers() bind(C)
-      end subroutine broadcast_triggers
-    end interface
-
-    ! Reduce tallies onto master process and accumulate
-    call time_tallies % start()
-    call accumulate_tallies()
-    call time_tallies % stop()
-
-    ! Reset global tally results
-    if (current_batch <= n_inactive) then
-      global_tallies(:,:) = ZERO
-      n_realizations = 0
-    end if
-
-    if (run_mode == MODE_EIGENVALUE) then
-      ! Perform CMFD calculation if on
-      if (cmfd_on) call execute_cmfd()
-      ! Write batch output
-      if (master .and. verbosity >= 7) call print_batch_keff()
-    end if
-
-    ! Check_triggers
-    if (master) call check_triggers()
-#ifdef OPENMC_MPI
-    call broadcast_triggers()
-#endif
-    if (satisfy_triggers .or. &
-         (trigger_on .and. current_batch == n_max_batches)) then
-      call statepoint_batch % add(current_batch)
-    end if
-
-    ! Write out state point if it's been specified for this batch
-    if (statepoint_batch % contains(current_batch)) then
-      if (sourcepoint_batch % contains(current_batch) .and. source_write &
-           .and. .not. source_separate) then
-        err = openmc_statepoint_write(write_source=.true._C_BOOL)
-      else
-        err = openmc_statepoint_write(write_source=.false._C_BOOL)
-      end if
-    end if
-
-    ! Write out a separate source point if it's been specified for this batch
-    if (sourcepoint_batch % contains(current_batch) .and. source_write &
-         .and. source_separate) call write_source_point()
-
-    ! Write a continously-overwritten source point if requested.
-    if (source_latest) then
-      filename = trim(path_output) // 'source' // '.h5'
-      call write_source_point(filename)
-    end if
-
-  end subroutine finalize_batch
-
-
 !===============================================================================
 ! INITIALIZE_SIMULATION
 !===============================================================================
diff --git a/src/simulation.cpp b/src/simulation.cpp
index c4f1f0387..2f154b8ce 100644
--- a/src/simulation.cpp
+++ b/src/simulation.cpp
@@ -1,6 +1,7 @@
 #include "openmc/simulation.h"
 
 #include "openmc/capi.h"
+#include "openmc/container_util.h"
 #include "openmc/eigenvalue.h"
 #include "openmc/error.h"
 #include "openmc/message_passing.h"
@@ -8,6 +9,7 @@
 #include "openmc/random_lcg.h"
 #include "openmc/settings.h"
 #include "openmc/source.h"
+#include "openmc/state_point.h"
 #include "openmc/timer.h"
 #include "openmc/tallies/tally.h"
 
@@ -15,13 +17,19 @@
 #include <string>
 
 // data/functions from Fortran side
+extern "C" bool cmfd_on;
+
+extern "C" void accumulate_tallies();
 extern "C" void allocate_banks();
+extern "C" void check_triggers();
 extern "C" void cmfd_init_batch();
 extern "C" void cmfd_tally_init();
 extern "C" void configure_tallies();
+extern "C" void execute_cmfd();
 extern "C" void init_tally_routines();
 extern "C" void join_bank_from_threads();
 extern "C" void load_state_point();
+extern "C" void print_batch_keff();
 extern "C" void print_columns();
 extern "C" void print_generation();
 extern "C" void print_results();
@@ -234,6 +242,62 @@ void initialize_batch()
   setup_active_tallies();
 }
 
+void finalize_batch()
+{
+  // Reduce tallies onto master process and accumulate
+  time_tallies.start();
+  accumulate_tallies();
+  time_tallies.stop();
+
+  // Reset global tally results
+  if (simulation::current_batch <= settings::n_inactive) {
+    auto gt = global_tallies();
+    std::fill(gt.begin(), gt.end(), 0.0);
+    n_realizations = 0;
+  }
+
+  if (settings::run_mode == RUN_MODE_EIGENVALUE) {
+    // Perform CMFD calculation if on
+    if (cmfd_on) execute_cmfd();
+    // Write batch output
+    if (mpi::master && settings::verbosity >= 7) print_batch_keff();
+  }
+
+  // Check_triggers
+  if (mpi::master) check_triggers();
+#ifdef OPENMC_MPI
+  MPI_Bcast(&simulation::satisfy_triggers, 1, MPI_C_BOOL, 0, mpi::intracomm);
+#endif
+  if (simulation::satisfy_triggers || (settings::trigger_on &&
+      simulation::current_batch == settings::n_max_batches)) {
+    settings::statepoint_batch.insert(simulation::current_batch);
+  }
+
+  // Write out state point if it's been specified for this batch
+  if (contains(settings::statepoint_batch, simulation::current_batch)) {
+    if (contains(settings::sourcepoint_batch, simulation::current_batch)
+        && settings::source_write && !settings::source_separate) {
+      bool b = true;
+      openmc_statepoint_write(nullptr, &b);
+    } else {
+      bool b = false;
+      openmc_statepoint_write(nullptr, &b);
+    }
+  }
+
+  // Write out a separate source point if it's been specified for this batch
+  if (contains(settings::sourcepoint_batch, simulation::current_batch)
+      && settings::source_write && settings::source_separate) {
+    write_source_point(nullptr);
+  }
+
+  // Write a continously-overwritten source point if requested.
+  if (settings::source_latest) {
+    auto filename = settings::path_output + "source.h5";
+    write_source_point(filename.c_str());
+  }
+}
+
 void initialize_generation()
 {
   if (settings::run_mode == RUN_MODE_EIGENVALUE) {
@@ -406,10 +470,6 @@ void broadcast_results() {
   simulation::k_abs_tra = temp[2];
 }
 
-void broadcast_triggers()
-{
-  MPI_Bcast(&simulation::satisfy_triggers, 1, MPI_C_BOOL, 0, mpi::intracomm);
-}
 #endif
 
 //==============================================================================
diff --git a/src/state_point.F90 b/src/state_point.F90
index 63a4d1dd4..d99bc1b84 100644
--- a/src/state_point.F90
+++ b/src/state_point.F90
@@ -443,44 +443,6 @@ contains
     end if
   end function openmc_statepoint_write
 
-!===============================================================================
-! WRITE_SOURCE_POINT
-!===============================================================================
-
-  subroutine write_source_point(filename)
-    character(MAX_FILE_LEN), intent(in), optional :: filename
-
-    logical :: parallel
-    integer(HID_T) :: file_id
-    character(MAX_FILE_LEN) :: filename_
-
-    ! When using parallel HDF5, the file is written to collectively by all
-    ! processes. With MPI-only, the file is opened and written by the master
-    ! (note that the call to write_source_bank is by all processes since slave
-    ! processes need to send source bank data to the master.
-#ifdef PHDF5
-    parallel = .true.
-#else
-    parallel = .false.
-#endif
-
-    if (present(filename)) then
-      filename_ = filename
-    else
-      filename_ = trim(path_output) // 'source.' // &
-           & zero_padded(current_batch, count_digits(n_max_batches))
-      filename_ = trim(filename_) // '.h5'
-    end if
-
-    if (master .or. parallel) then
-      file_id = file_open(filename_, 'w', parallel=.true.)
-      call write_attribute(file_id, "filetype", 'source')
-    end if
-    call write_source_bank(file_id, source_bank)
-    if (master .or. parallel) call file_close(file_id)
-
-  end subroutine write_source_point
-
 !===============================================================================
 ! LOAD_STATE_POINT
 !===============================================================================
diff --git a/src/state_point.cpp b/src/state_point.cpp
index ed679836b..620b5f71a 100644
--- a/src/state_point.cpp
+++ b/src/state_point.cpp
@@ -1,6 +1,7 @@
 #include "openmc/state_point.h"
 
 #include <algorithm>
+#include <iomanip> // for setfill, setw
 #include <string>
 #include <vector>
 
@@ -19,7 +20,6 @@
 
 namespace openmc {
 
-
 hid_t h5banktype() {
   // Create type for array of 3 reals
   hsize_t dims[] {3};
@@ -37,6 +37,46 @@ hid_t h5banktype() {
   return banktype;
 }
 
+void
+write_source_point(const char* filename)
+{
+  // When using parallel HDF5, the file is written to collectively by all
+  // processes. With MPI-only, the file is opened and written by the master
+  // (note that the call to write_source_bank is by all processes since slave
+  // processes need to send source bank data to the master.
+#ifdef PHDF5
+  bool parallel = true;
+#else
+  bool parallel = false;
+#endif
+
+  std::string filename_;
+  if (filename) {
+    filename_ = filename;
+  } else {
+    // Determine width for zero padding
+    int w = std::to_string(settings::n_max_batches).size();
+
+    std::stringstream s;
+    s << settings::path_output << "source." << std::setfill('0')
+      << std::setw(w) << simulation::current_batch << ".h5";
+    filename_ = s.str();
+  }
+
+  hid_t file_id;
+  if (mpi::master || parallel) {
+    file_id = file_open(filename_, 'w', true);
+    write_attribute(file_id, "filetype", "source");
+  }
+
+  // Get pointer to source bank and write to file
+  Bank* source_bank;
+  int64_t n;
+  openmc_source_bank(&source_bank, &n);
+  write_source_bank(file_id, source_bank);
+
+  if (mpi::master || parallel) file_close(file_id);
+}
 
 void
 write_source_bank(hid_t group_id, Bank* source_bank)
diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90
index 8dfcf5934..19d99a014 100644
--- a/src/tallies/tally.F90
+++ b/src/tallies/tally.F90
@@ -3769,7 +3769,7 @@ contains
 ! within the batch to a new random variable
 !===============================================================================
 
-  subroutine accumulate_tallies()
+  subroutine accumulate_tallies() bind(C)
 
     integer :: i
     real(C_DOUBLE) :: k_col ! Copy of batch collision estimate of keff
diff --git a/src/tallies/trigger.F90 b/src/tallies/trigger.F90
index 27aa5ec50..715b3c206 100644
--- a/src/tallies/trigger.F90
+++ b/src/tallies/trigger.F90
@@ -29,7 +29,7 @@ contains
 ! and predicts the number of remainining batches to convergence.
 !===============================================================================
 
-  subroutine check_triggers()
+  subroutine check_triggers() bind(C)
 
     implicit none
 
diff --git a/src/timer.cpp b/src/timer.cpp
index 154e23db1..df8f4fb92 100644
--- a/src/timer.cpp
+++ b/src/timer.cpp
@@ -45,6 +45,7 @@ Timer time_bank_sendrecv;
 Timer time_finalize;
 Timer time_inactive;
 Timer time_initialize;
+Timer time_tallies;
 Timer time_total;
 
 //==============================================================================
@@ -58,6 +59,7 @@ extern "C" double time_bank_sendrecv_elapsed() { return time_bank_sendrecv.elaps
 extern "C" double time_finalize_elapsed() { return time_finalize.elapsed(); }
 extern "C" double time_inactive_elapsed() { return time_inactive.elapsed(); }
 extern "C" double time_initialize_elapsed() { return time_initialize.elapsed(); }
+extern "C" double time_tallies_elapsed() { return time_tallies.elapsed(); }
 extern "C" double time_total_elapsed() { return time_total.elapsed(); }
 
 extern "C" void reset_timers()
diff --git a/src/timer_header.F90 b/src/timer_header.F90
index ccf872e13..00af1596f 100644
--- a/src/timer_header.F90
+++ b/src/timer_header.F90
@@ -35,6 +35,10 @@ module timer_header
       import C_DOUBLE
       real(C_DOUBLE) :: t
     end function
+    function time_tallies_elapsed() result(t) bind(C)
+      import C_DOUBLE
+      real(C_DOUBLE) :: t
+    end function
     function time_total_elapsed() result(t) bind(C)
       import C_DOUBLE
       real(C_DOUBLE) :: t