diff --git a/openmc/capi.py b/openmc/capi.py
index bb9241d7d3..f4c9a7f3ce 100644
--- a/openmc/capi.py
+++ b/openmc/capi.py
@@ -43,6 +43,8 @@ class OpenMCLibrary(object):
         self._dll.openmc_find.restype = None
         self._dll.openmc_init.argtypes = [POINTER(c_int)]
         self._dll.openmc_init.restype = None
+        self._dll.openmc_get_keff.argtypes = [POINTER(c_double*2)]
+        self._dll.openmc_get_keff.restype = c_int
         self._dll.openmc_load_nuclide.argtypes = [c_char_p]
         self._dll.openmc_load_nuclide.restype = c_int
         self._dll.openmc_material_add_nuclide.argtypes = [
@@ -141,7 +143,21 @@ class OpenMCLibrary(object):
         else:
             return self._dll.openmc_init(None)
 
+    def keff(self):
+        """Return the calculated k-eigenvalue and its standard deviation.
+
+        Returns
+        -------
+        tuple
+            Mean k-eigenvalue and standard deviation of the mean
+
+        """
+        k = (c_double*2)()
+        err = self._dll.openmc_get_keff(k)
+        return tuple(k)
+
     def load_nuclide(self, name):
+
         """Load cross section data for a nuclide.
 
         Parameters
diff --git a/src/api.F90 b/src/api.F90
index b00193fcff..45f511f9cc 100644
--- a/src/api.F90
+++ b/src/api.F90
@@ -5,7 +5,7 @@ module openmc_api
   use hdf5, only: HID_T
 
   use constants,       only: K_BOLTZMANN
-  use eigenvalue,      only: k_sum
+  use eigenvalue,      only: k_sum, openmc_get_keff
   use finalize,        only: openmc_finalize
   use geometry,        only: find_cell
   use global
@@ -23,6 +23,7 @@ module openmc_api
   public :: openmc_cell_set_temperature
   public :: openmc_finalize
   public :: openmc_find
+  public :: openmc_get_keff
   public :: openmc_init
   public :: openmc_load_nuclide
   public :: openmc_material_add_nuclide
diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90
index f370ab38eb..9f4da7c82f 100644
--- a/src/eigenvalue.F90
+++ b/src/eigenvalue.F90
@@ -1,5 +1,6 @@
 module eigenvalue
 
+  use, intrinsic :: ISO_C_BINDING
 
   use algorithm,   only: binary_search
   use constants,   only: ZERO
@@ -439,8 +440,8 @@ contains
   end subroutine calculate_average_keff
 
 !===============================================================================
-! CALCULATE_COMBINED_KEFF calculates a minimum variance estimate of k-effective
-! based on a linear combination of the collision, absorption, and tracklength
+! OPENMC_GET_KEFF calculates a minimum variance estimate of k-effective based on
+! a linear combination of the collision, absorption, and tracklength
 ! estimates. The theory behind this can be found in M. Halperin, "Almost
 ! linearly-optimum combination of unbiased estimates," J. Am. Stat. Assoc., 56,
 ! 36-43 (1961), doi:10.1080/01621459.1961.10482088. The implementation here
@@ -448,7 +449,9 @@ contains
 ! of keff confidence intervals in MCNP," Nucl. Technol., 111, 169-182 (1995).
 !===============================================================================
 
-  subroutine calculate_combined_keff()
+  function openmc_get_keff(k_combined) result(err) bind(C)
+    real(C_DOUBLE), intent(out) :: k_combined(2)
+    integer(C_INT) :: err
 
     integer :: l        ! loop index
     integer :: i, j, k  ! indices referring to collision, absorption, or track
@@ -459,9 +462,14 @@ contains
     real(8) :: g        ! sum of weighting factors
     real(8) :: S(3)     ! sums used for variance calculation
 
+    k_combined = ZERO
+
     ! Make sure we have at least four realizations. Notice that at the end,
     ! there is a N-3 term in a denominator.
-    if (n_realizations <= 3) return
+    if (n_realizations <= 3) then
+      err = -1
+      return
+    end if
 
     ! Initialize variables
     n = real(n_realizations, 8)
@@ -523,7 +531,6 @@ contains
       ! Initialize variables
       g = ZERO
       S = ZERO
-      k_combined = ZERO
 
       do l = 1, 3
         ! Permutations of estimates
@@ -590,8 +597,9 @@ contains
       k_combined(2) = sqrt(k_combined(2))
 
     end if
+    err = 0
 
-  end subroutine calculate_combined_keff
+  end function openmc_get_keff
 
 !===============================================================================
 ! COUNT_SOURCE_FOR_UFS determines the source fraction in each UFS mesh cell and
diff --git a/src/global.F90 b/src/global.F90
index 2690ed1518..41f8654694 100644
--- a/src/global.F90
+++ b/src/global.F90
@@ -256,7 +256,6 @@ module global
   real(8) :: k_col_abs = ZERO ! sum over batches of k_collision * k_absorption
   real(8) :: k_col_tra = ZERO ! sum over batches of k_collision * k_tracklength
   real(8) :: k_abs_tra = ZERO ! sum over batches of k_absorption * k_tracklength
-  real(8) :: k_combined(2)    ! combined best estimate of k-effective
 
   ! Shannon entropy
   logical :: entropy_on = .false.
diff --git a/src/output.F90 b/src/output.F90
index 828af717d8..9835225e06 100644
--- a/src/output.F90
+++ b/src/output.F90
@@ -1,8 +1,10 @@
 module output
 
+  use, intrinsic :: ISO_C_BINDING
   use, intrinsic :: ISO_FORTRAN_ENV
 
   use constants
+  use eigenvalue,      only: openmc_get_keff
   use endf,            only: reaction_name
   use error,           only: fatal_error, warning
   use geometry_header, only: Cell, Universe, Lattice, RectLattice, &
@@ -608,6 +610,8 @@ contains
     real(8) :: t_n1    ! t-value with N-1 degrees of freedom
     real(8) :: t_n3    ! t-value with N-3 degrees of freedom
     real(8) :: x(2)    ! mean and standard deviation
+    real(C_DOUBLE) :: k_combined(2)
+    integer(C_INT) :: err
 
     ! display header block for results
     call header("Results", 4)
@@ -634,8 +638,11 @@ contains
           write(ou,102) "k-effective (Track-length)", x(1), t_n1 * x(2)
           x(:) = mean_stdev(r(:, K_ABSORPTION), n)
           write(ou,102) "k-effective (Absorption)", x(1), t_n1 * x(2)
-          if (n > 3) write(ou,102) "Combined k-effective", k_combined(1), &
-               t_n3 * k_combined(2)
+          if (n > 3) then
+            err = openmc_get_keff(k_combined)
+            write(ou,102) "Combined k-effective", k_combined(1), &
+                 t_n3 * k_combined(2)
+          end if
         end if
         x(:) = mean_stdev(global_tallies(:, LEAKAGE), n)
         write(ou,102) "Leakage Fraction", x(1), t_n1 * x(2)
diff --git a/src/simulation.F90 b/src/simulation.F90
index 87a17452d9..505d5b528f 100644
--- a/src/simulation.F90
+++ b/src/simulation.F90
@@ -5,9 +5,8 @@ module simulation
   use cmfd_execute,    only: cmfd_init_batch, execute_cmfd
   use constants,       only: ZERO
   use eigenvalue,      only: count_source_for_ufs, calculate_average_keff, &
-                             calculate_combined_keff, calculate_generation_keff, &
-                             shannon_entropy, synchronize_bank, keff_generation, &
-                             k_sum
+                             calculate_generation_keff, shannon_entropy, &
+                             synchronize_bank, keff_generation, k_sum
 #ifdef _OPENMP
   use eigenvalue,      only: join_bank_from_threads
 #endif
@@ -300,9 +299,6 @@ contains
     if (run_mode == MODE_EIGENVALUE) then
       ! Perform CMFD calculation if on
       if (cmfd_on) call execute_cmfd()
-
-      ! Calculate combined estimate of k-effective
-      if (master) call calculate_combined_keff()
     end if
 
     ! Check_triggers
@@ -327,13 +323,6 @@ contains
       call write_source_point()
     end if
 
-    if (master .and. current_batch == n_max_batches .and. &
-         run_mode == MODE_EIGENVALUE) then
-      ! Make sure combined estimate of k-effective is calculated at the last
-      ! batch in case no state point is written
-      call calculate_combined_keff()
-    end if
-
   end subroutine finalize_batch
 
 !===============================================================================
diff --git a/src/state_point.F90 b/src/state_point.F90
index e4c9a43ae7..f48b5638b4 100644
--- a/src/state_point.F90
+++ b/src/state_point.F90
@@ -11,11 +11,12 @@ module state_point
 ! intervals, using the <state_point ... /> tag.
 !===============================================================================
 
-  use, intrinsic :: ISO_C_BINDING, only: c_loc, c_ptr
+  use, intrinsic :: ISO_C_BINDING
 
   use hdf5
 
   use constants
+  use eigenvalue,         only: openmc_get_keff
   use endf,               only: reaction_name
   use error,              only: fatal_error, warning
   use global
@@ -46,6 +47,8 @@ contains
     integer(HID_T) :: cmfd_group, tallies_group, tally_group, meshes_group, &
                       mesh_group, filters_group, filter_group, derivs_group, &
                       deriv_group, runtime_group
+    integer(C_INT) :: err
+    real(C_DOUBLE) :: k_combined(2)
     character(MAX_WORD_LEN), allocatable :: str_array(:)
     character(MAX_FILE_LEN)    :: filename
     type(TallyObject), pointer    :: tally
@@ -117,6 +120,7 @@ contains
         call write_dataset(file_id, "k_col_abs", k_col_abs)
         call write_dataset(file_id, "k_col_tra", k_col_tra)
         call write_dataset(file_id, "k_abs_tra", k_abs_tra)
+        err = openmc_get_keff(k_combined)
         call write_dataset(file_id, "k_combined", k_combined)
 
         ! Write out CMFD info
@@ -647,7 +651,6 @@ contains
     integer(HID_T) :: cmfd_group
     integer(HID_T) :: tallies_group
     integer(HID_T) :: tally_group
-    real(8) :: real_array(3)
     logical :: source_present
     character(MAX_WORD_LEN) :: word
 
@@ -727,7 +730,6 @@ contains
       call read_dataset(k_col_abs, file_id, "k_col_abs")
       call read_dataset(k_col_tra, file_id, "k_col_tra")
       call read_dataset(k_abs_tra, file_id, "k_abs_tra")
-      call read_dataset(real_array(1:2), file_id, "k_combined")
 
       ! Take maximum of statepoint n_inactive and input n_inactive
       n_inactive = max(n_inactive, int_array(1))
diff --git a/src/trigger.F90 b/src/trigger.F90
index 9eb39e9cc2..dc0a10b153 100644
--- a/src/trigger.F90
+++ b/src/trigger.F90
@@ -1,10 +1,13 @@
 module trigger
 
+  use, intrinsic :: ISO_C_BINDING
+
 #ifdef MPI
   use message_passing
 #endif
 
   use constants
+  use eigenvalue,     only: openmc_get_keff
   use global
   use string,         only: to_str
   use output,         only: warning, write_message
@@ -101,10 +104,12 @@ contains
     integer :: score_index    ! scoring bin index
     integer :: n_order        ! loop index for moment orders
     integer :: nm_order       ! loop index for Ynm moment orders
+    integer(C_INT) :: err
     real(8) :: uncertainty    ! trigger uncertainty
     real(8) :: std_dev = ZERO ! trigger standard deviation
     real(8) :: rel_err = ZERO ! trigger relative error
     real(8) :: ratio          ! ratio of the uncertainty/trigger threshold
+    real(C_DOUBLE) :: k_combined(2)
     type(TallyObject), pointer     :: t               ! tally pointer
     type(TriggerObject), pointer   :: trigger         ! tally trigger
 
@@ -119,6 +124,7 @@ contains
       ! Check eigenvalue trigger
       if (run_mode == MODE_EIGENVALUE) then
         if (keff_trigger % trigger_type /= 0) then
+          err = openmc_get_keff(k_combined)
           select case (keff_trigger % trigger_type)
           case(VARIANCE)
             uncertainty = k_combined(2) ** 2