diff --git a/include/openmc/bank.h b/include/openmc/bank.h
index 95386514d7b..fd8fbd73ee5 100644
--- a/include/openmc/bank.h
+++ b/include/openmc/bank.h
@@ -22,6 +22,14 @@ extern SharedArray<SourceSite> surf_source_bank;
 
 extern SharedArray<SourceSite> fission_bank;
 
+extern vector<vector<int>> ifp_source_delayed_group_bank;
+
+extern vector<vector<double>> ifp_source_lifetime_bank;
+
+extern vector<vector<int>> ifp_fission_delayed_group_bank;
+
+extern vector<vector<double>> ifp_fission_lifetime_bank;
+
 extern vector<int64_t> progeny_per_particle;
 
 } // namespace simulation
diff --git a/include/openmc/constants.h b/include/openmc/constants.h
index 605ae1839d8..97d0218d1a0 100644
--- a/include/openmc/constants.h
+++ b/include/openmc/constants.h
@@ -308,7 +308,10 @@ enum TallyScore {
   SCORE_FISS_Q_PROMPT = -14,      // prompt fission Q-value
   SCORE_FISS_Q_RECOV = -15,       // recoverable fission Q-value
   SCORE_DECAY_RATE = -16,         // delayed neutron precursor decay rate
-  SCORE_PULSE_HEIGHT = -17        // pulse-height
+  SCORE_PULSE_HEIGHT = -17,       // pulse-height
+  SCORE_IFP_TIME_NUM = -18,       // IFP lifetime numerator
+  SCORE_IFP_BETA_NUM = -19,       // IFP delayed fraction numerator
+  SCORE_IFP_DENOM = -20           // IFP common denominator
 };
 
 // Global tally parameters
diff --git a/include/openmc/eigenvalue.h b/include/openmc/eigenvalue.h
index b456fee21eb..7d818ac16f7 100644
--- a/include/openmc/eigenvalue.h
+++ b/include/openmc/eigenvalue.h
@@ -79,6 +79,16 @@ void write_eigenvalue_hdf5(hid_t group);
 //! \param[in] group HDF5 group
 void read_eigenvalue_hdf5(hid_t group);
 
+//==============================================================================
+// Type definitions
+//==============================================================================
+
+//! Deserialization info for IFP
+struct DeserializationInfo {
+  int64_t index_local;
+  int64_t n;
+};
+
 } // namespace openmc
 
 #endif // OPENMC_EIGENVALUE_H
diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h
index 164148cce10..fbe75f6d57c 100644
--- a/include/openmc/particle_data.h
+++ b/include/openmc/particle_data.h
@@ -424,6 +424,9 @@ class ParticleData : public GeometryState {
 
   int cell_born_ {-1};
 
+  // Iterated Fission Probability
+  double lifetime_ {0.0}; //!< neutron lifetime [s]
+
   int n_collision_ {0};
 
   bool write_track_ {false};
@@ -522,6 +525,10 @@ class ParticleData : public GeometryState {
   double& time_last() { return time_last_; }
   const double& time_last() const { return time_last_; }
 
+  // Particle lifetime
+  double& lifetime() { return lifetime_; }
+  const double& lifetime() const { return lifetime_; }
+
   // What event took place, described in greater detail below
   TallyEvent& event() { return event_; }
   const TallyEvent& event() const { return event_; }
diff --git a/include/openmc/settings.h b/include/openmc/settings.h
index a95f1ced9f1..92600d81c04 100644
--- a/include/openmc/settings.h
+++ b/include/openmc/settings.h
@@ -24,6 +24,14 @@ enum class SSWCellType {
   To,
 };
 
+// Type of IFP parameters
+enum class IFPParameter {
+  None,
+  Both,
+  BetaEffective,
+  GenerationTime,
+};
+
 //==============================================================================
 // Global variable declarations
 //==============================================================================
@@ -42,7 +50,8 @@ extern bool
   delayed_photon_scaling;   //!< Scale fission photon yield to include delayed
 extern "C" bool entropy_on; //!< calculate Shannon entropy?
 extern "C" bool
-  event_based; //!< use event-based mode (instead of history-based)
+  event_based;   //!< use event-based mode (instead of history-based)
+extern bool ifp; //!< Use IFP for kinetics parameters?
 extern bool legendre_to_tabular; //!< convert Legendre distributions to tabular?
 extern bool material_cell_offsets; //!< create material cells offsets?
 extern "C" bool output_summary;    //!< write summary.h5?
@@ -109,6 +118,10 @@ extern array<double, 4>
   energy_cutoff; //!< Energy cutoff in [eV] for each particle type
 extern array<double, 4>
   time_cutoff; //!< Time cutoff in [s] for each particle type
+extern int
+  ifp_n_generation; //!< Number of generation for Iterated Fission Probability
+extern IFPParameter
+  ifp_parameter; //!< Parameter to calculate for Iterated Fission Probability
 extern int
   legendre_to_tabular_points; //!< number of points to convert Legendres
 extern int max_order;         //!< Maximum Legendre order for multigroup data
diff --git a/openmc/settings.py b/openmc/settings.py
index ddaac040191..6bf04654ea2 100644
--- a/openmc/settings.py
+++ b/openmc/settings.py
@@ -79,6 +79,13 @@ class Settings:
         .. versionadded:: 0.12
     generations_per_batch : int
         Number of generations per batch
+    iterated_fission_probability: dict
+        Dictionary indicating the Iterated Fission Probability parameters.
+        Acceptable keys are:
+
+        :parameter: Kinetics parameter to calculate between 'beta_effective' and
+                    'generation_time', or 'both' (str)
+        :n_generation: Number of generation (int)
     max_lost_particles : int
         Maximum number of lost particles
 
@@ -335,6 +342,9 @@ def __init__(self, **kwargs):
 
         self._output = None
 
+        # Iterated Fission Probability
+        self._iterated_fission_probability = {}
+
         # Output options
         self._statepoint = {}
         self._sourcepoint = {}
@@ -760,6 +770,27 @@ def verbosity(self, verbosity: int):
         cv.check_less_than('verbosity', verbosity, 10, True)
         self._verbosity = verbosity
 
+    @property
+    def iterated_fission_probability(self) -> dict:
+        return self._iterated_fission_probability
+
+    @iterated_fission_probability.setter
+    def iterated_fission_probability(self, iterated_fission_probability: dict):
+        cv.check_type(
+            "Iterated Fission Probability options",
+            iterated_fission_probability,
+            Mapping,
+        )
+        for key, value in iterated_fission_probability.items():
+            cv.check_value("Iterated Fission Probability key", key, ["parameter", "n_generation"])
+            if key == "parameter":
+                cv.check_value("parameter", value, ["beta_effective", "generation_time", "both"])
+            elif key == "n_generation":
+                cv.check_type("number of generation", value, Integral)
+                cv.check_greater_than("number of generation", value, 0)
+
+        self._iterated_fission_probability = iterated_fission_probability
+
     @property
     def tabular_legendre(self) -> dict:
         return self._tabular_legendre
@@ -1347,6 +1378,16 @@ def _create_no_reduce_subelement(self, root):
             element = ET.SubElement(root, "no_reduce")
             element.text = str(self._no_reduce).lower()
 
+    def _create_iterated_fission_probability_subelements(self, root):
+        if self.iterated_fission_probability:
+            element = ET.SubElement(root, "iterated_fission_probability")
+            if 'parameter' in self._iterated_fission_probability:
+                subelement = ET.SubElement(element, "parameter")
+                subelement.text = str(self._iterated_fission_probability['parameter'])
+            if 'n_generation' in self._iterated_fission_probability:
+                subelement = ET.SubElement(element, "n_generation")
+                subelement.text = str(self._iterated_fission_probability['n_generation'])
+
     def _create_tabular_legendre_subelements(self, root):
         if self.tabular_legendre:
             element = ET.SubElement(root, "tabular_legendre")
@@ -1747,6 +1788,16 @@ def _verbosity_from_xml_element(self, root):
         if text is not None:
             self.verbosity = int(text)
 
+    def _iterated_fission_probability_from_xml_element(self, root):
+        elem = root.find('iterated_fission_probability')
+        if elem is not None:
+            text = get_text(elem, 'parameter')
+            if text is not None:
+                self.iterated_fission_probability['parameter'] = text
+            text = get_text(elem, 'n_generation')
+            if text is not None:
+                self.iterated_fission_probability['n_generation'] = int(text)
+
     def _tabular_legendre_from_xml_element(self, root):
         elem = root.find('tabular_legendre')
         if elem is not None:
@@ -1946,6 +1997,7 @@ def to_xml_element(self, mesh_memo=None):
         self._create_trigger_subelement(element)
         self._create_no_reduce_subelement(element)
         self._create_verbosity_subelement(element)
+        self._create_iterated_fission_probability_subelements(element)
         self._create_tabular_legendre_subelements(element)
         self._create_temperature_subelements(element)
         self._create_trace_subelement(element)
@@ -2052,6 +2104,7 @@ def from_xml_element(cls, elem, meshes=None):
         settings._trigger_from_xml_element(elem)
         settings._no_reduce_from_xml_element(elem)
         settings._verbosity_from_xml_element(elem)
+        settings._iterated_fission_probability_from_xml_element(elem)
         settings._tabular_legendre_from_xml_element(elem)
         settings._temperature_from_xml_element(elem)
         settings._trace_from_xml_element(elem)
diff --git a/src/bank.cpp b/src/bank.cpp
index 8d00d544090..5fff1e3458f 100644
--- a/src/bank.cpp
+++ b/src/bank.cpp
@@ -26,6 +26,14 @@ SharedArray<SourceSite> surf_source_bank;
 // function.
 SharedArray<SourceSite> fission_bank;
 
+vector<vector<int>> ifp_source_delayed_group_bank;
+
+vector<vector<double>> ifp_source_lifetime_bank;
+
+vector<vector<int>> ifp_fission_delayed_group_bank;
+
+vector<vector<double>> ifp_fission_lifetime_bank;
+
 // Each entry in this vector corresponds to the number of progeny produced
 // this generation for the particle located at that index. This vector is
 // used to efficiently sort the fission bank after each iteration.
@@ -82,6 +90,10 @@ void sort_fission_bank()
   // over provisioned, so we can use that as scratch space.
   SourceSite* sorted_bank;
   vector<SourceSite> sorted_bank_holder;
+  vector<int>* sorted_ifp_delayed_group_bank;
+  vector<vector<int>> sorted_ifp_delayed_group_bank_holder;
+  vector<double>* sorted_ifp_lifetime_bank;
+  vector<vector<double>> sorted_ifp_lifetime_bank_holder;
 
   // If there is not enough space, allocate a temporary vector and point to it
   if (simulation::fission_bank.size() >
@@ -92,6 +104,21 @@ void sort_fission_bank()
     sorted_bank = &simulation::fission_bank[simulation::fission_bank.size()];
   }
 
+  if (settings::ifp) {
+    if (settings::ifp_parameter == IFPParameter::BetaEffective ||
+        settings::ifp_parameter == IFPParameter::Both) {
+      sorted_ifp_delayed_group_bank_holder.resize(
+        simulation::fission_bank.size());
+      sorted_ifp_delayed_group_bank =
+        sorted_ifp_delayed_group_bank_holder.data();
+    }
+    if (settings::ifp_parameter == IFPParameter::GenerationTime ||
+        settings::ifp_parameter == IFPParameter::Both) {
+      sorted_ifp_lifetime_bank_holder.resize(simulation::fission_bank.size());
+      sorted_ifp_lifetime_bank = sorted_ifp_lifetime_bank_holder.data();
+    }
+  }
+
   // Use parent and progeny indices to sort fission bank
   for (int64_t i = 0; i < simulation::fission_bank.size(); i++) {
     const auto& site = simulation::fission_bank[i];
@@ -102,11 +129,38 @@ void sort_fission_bank()
                   "shared fission bank size.");
     }
     sorted_bank[idx] = site;
+    if (settings::ifp) {
+      if (settings::ifp_parameter == IFPParameter::BetaEffective ||
+          settings::ifp_parameter == IFPParameter::Both) {
+        const auto& ifp_delayed_groups =
+          simulation::ifp_fission_delayed_group_bank[i];
+        sorted_ifp_delayed_group_bank[idx] = ifp_delayed_groups;
+      }
+      if (settings::ifp_parameter == IFPParameter::GenerationTime ||
+          settings::ifp_parameter == IFPParameter::Both) {
+        const auto& ifp_lifetimes = simulation::ifp_fission_lifetime_bank[i];
+        sorted_ifp_lifetime_bank[idx] = ifp_lifetimes;
+      }
+    }
   }
 
   // Copy sorted bank into the fission bank
   std::copy(sorted_bank, sorted_bank + simulation::fission_bank.size(),
     simulation::fission_bank.data());
+  if (settings::ifp) {
+    if (settings::ifp_parameter == IFPParameter::BetaEffective ||
+        settings::ifp_parameter == IFPParameter::Both) {
+      std::copy(sorted_ifp_delayed_group_bank,
+        sorted_ifp_delayed_group_bank + simulation::fission_bank.size(),
+        simulation::ifp_fission_delayed_group_bank.data());
+    }
+    if (settings::ifp_parameter == IFPParameter::GenerationTime ||
+        settings::ifp_parameter == IFPParameter::Both) {
+      std::copy(sorted_ifp_lifetime_bank,
+        sorted_ifp_lifetime_bank + simulation::fission_bank.size(),
+        simulation::ifp_fission_lifetime_bank.data());
+    }
+  }
 }
 
 //==============================================================================
diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp
index 8669d76f947..69a05cf4a0a 100644
--- a/src/eigenvalue.cpp
+++ b/src/eigenvalue.cpp
@@ -153,11 +153,41 @@ void synchronize_bank()
   // Allocate temporary source bank -- we don't really know how many fission
   // sites were created, so overallocate by a factor of 3
   int64_t index_temp = 0;
+
   vector<SourceSite> temp_sites(3 * simulation::work_per_rank);
+  vector<vector<int>> temp_delayed_groups;
+  vector<vector<double>> temp_lifetimes;
+  if (settings::ifp) {
+    if (settings::ifp_parameter == IFPParameter::BetaEffective ||
+        settings::ifp_parameter == IFPParameter::Both) {
+      temp_delayed_groups.resize(3 * simulation::work_per_rank);
+    }
+    if (settings::ifp_parameter == IFPParameter::GenerationTime ||
+        settings::ifp_parameter == IFPParameter::Both) {
+      temp_lifetimes.resize(3 * simulation::work_per_rank);
+    }
+  }
 
   for (int64_t i = 0; i < simulation::fission_bank.size(); i++) {
     const auto& site = simulation::fission_bank[i];
 
+    // Declare pointer to constant IFP data that will be initialized if
+    // ifp is requested by the user.
+    const vector<int>* delayed_groups_ptr;
+    const vector<double>* lifetimes_ptr;
+
+    // Initialize IFP data pointer
+    if (settings::ifp) {
+      if (settings::ifp_parameter == IFPParameter::BetaEffective ||
+          settings::ifp_parameter == IFPParameter::Both) {
+        delayed_groups_ptr = &simulation::ifp_fission_delayed_group_bank[i];
+      }
+      if (settings::ifp_parameter == IFPParameter::GenerationTime ||
+          settings::ifp_parameter == IFPParameter::Both) {
+        lifetimes_ptr = &simulation::ifp_fission_lifetime_bank[i];
+      }
+    }
+
     // If there are less than n_particles particles banked, automatically add
     // int(n_particles/total) sites to temp_sites. For example, if you need
     // 1000 and 300 were banked, this would add 3 source sites per banked site
@@ -165,6 +195,16 @@ void synchronize_bank()
     if (total < settings::n_particles) {
       for (int64_t j = 1; j <= settings::n_particles / total; ++j) {
         temp_sites[index_temp] = site;
+        if (settings::ifp) {
+          if (settings::ifp_parameter == IFPParameter::BetaEffective ||
+              settings::ifp_parameter == IFPParameter::Both) {
+            temp_delayed_groups[index_temp] = *delayed_groups_ptr;
+          }
+          if (settings::ifp_parameter == IFPParameter::GenerationTime ||
+              settings::ifp_parameter == IFPParameter::Both) {
+            temp_lifetimes[index_temp] = *lifetimes_ptr;
+          }
+        }
         ++index_temp;
       }
     }
@@ -172,6 +212,16 @@ void synchronize_bank()
     // Randomly sample sites needed
     if (prn(&seed) < p_sample) {
       temp_sites[index_temp] = site;
+      if (settings::ifp) {
+        if (settings::ifp_parameter == IFPParameter::BetaEffective ||
+            settings::ifp_parameter == IFPParameter::Both) {
+          temp_delayed_groups[index_temp] = *delayed_groups_ptr;
+        }
+        if (settings::ifp_parameter == IFPParameter::GenerationTime ||
+            settings::ifp_parameter == IFPParameter::Both) {
+          temp_lifetimes[index_temp] = *lifetimes_ptr;
+        }
+      }
       ++index_temp;
     }
   }
@@ -188,6 +238,8 @@ void synchronize_bank()
   MPI_Exscan(&index_temp, &start, 1, MPI_INT64_T, MPI_SUM, mpi::intracomm);
   finish = start + index_temp;
 
+  // TODO: protect for MPI_Exscan at rank 0
+
   // Allocate space for bank_position if this hasn't been done yet
   int64_t bank_position[mpi::n_procs];
   MPI_Allgather(
@@ -211,9 +263,23 @@ void synchronize_bank()
       // If we have too few sites, repeat sites from the very end of the
       // fission bank
       sites_needed = settings::n_particles - finish;
+      // TODO: sites_needed > simulation::fission_bank.size() or other test to
+      // make sure we don't need info from other proc
       for (int i = 0; i < sites_needed; ++i) {
         int i_bank = simulation::fission_bank.size() - sites_needed + i;
         temp_sites[index_temp] = simulation::fission_bank[i_bank];
+        if (settings::ifp) {
+          if (settings::ifp_parameter == IFPParameter::BetaEffective ||
+              settings::ifp_parameter == IFPParameter::Both) {
+            temp_delayed_groups[index_temp] =
+              simulation::ifp_fission_delayed_group_bank[i_bank];
+          }
+          if (settings::ifp_parameter == IFPParameter::GenerationTime ||
+              settings::ifp_parameter == IFPParameter::Both) {
+            temp_lifetimes[index_temp] =
+              simulation::ifp_fission_lifetime_bank[i_bank];
+          }
+        }
         ++index_temp;
       }
     }
@@ -229,15 +295,47 @@ void synchronize_bank()
   // ==========================================================================
   // SEND BANK SITES TO NEIGHBORS
 
+  // If IFP, broadcast the number of generation from the size of the first
+  // element
+  int ifp_n_generation;
+  if (settings::ifp) {
+    // TODO:
+    if (mpi::rank == 0) {
+      if (settings::ifp_parameter == IFPParameter::BetaEffective ||
+          settings::ifp_parameter == IFPParameter::Both) {
+        ifp_n_generation = static_cast<int>(temp_delayed_groups[0].size());
+      } else {
+        ifp_n_generation = static_cast<int>(temp_lifetimes[0].size());
+      }
+    }
+    MPI_Bcast(&ifp_n_generation, 1, MPI_INT, 0, mpi::intracomm);
+  }
+
   int64_t index_local = 0;
   vector<MPI_Request> requests;
 
+  vector<int> send_delayed_groups;
+  vector<double> send_lifetimes;
+
   if (start < settings::n_particles) {
     // Determine the index of the processor which has the first part of the
     // source_bank for the local processor
     int neighbor = upper_bound_index(
       simulation::work_index.begin(), simulation::work_index.end(), start);
 
+    // Resize IFP send buffers
+    if (settings::ifp && mpi::n_procs > 1) {
+      if (settings::ifp_parameter == IFPParameter::BetaEffective ||
+          settings::ifp_parameter == IFPParameter::Both) {
+        send_delayed_groups.resize(
+          ifp_n_generation * 3 * simulation::work_per_rank);
+      }
+      if (settings::ifp_parameter == IFPParameter::GenerationTime ||
+          settings::ifp_parameter == IFPParameter::Both) {
+        send_lifetimes.resize(ifp_n_generation * 3 * simulation::work_per_rank);
+      }
+    }
+
     while (start < finish) {
       // Determine the number of sites to send
       int64_t n =
@@ -250,6 +348,41 @@ void synchronize_bank()
         MPI_Isend(&temp_sites[index_local], static_cast<int>(n),
           mpi::source_site, neighbor, mpi::rank, mpi::intracomm,
           &requests.back());
+
+        if (settings::ifp) {
+
+          for (int i = index_local; i < index_local + n; i++) {
+            if (settings::ifp_parameter == IFPParameter::BetaEffective ||
+                settings::ifp_parameter == IFPParameter::Both) {
+              std::copy(temp_delayed_groups[i].begin(),
+                temp_delayed_groups[i].end(),
+                send_delayed_groups.begin() + i * ifp_n_generation);
+            }
+            if (settings::ifp_parameter == IFPParameter::GenerationTime ||
+                settings::ifp_parameter == IFPParameter::Both) {
+              std::copy(temp_lifetimes[i].begin(), temp_lifetimes[i].end(),
+                send_lifetimes.begin() + i * ifp_n_generation);
+            }
+          }
+
+          // Send delayed groups
+          if (settings::ifp_parameter == IFPParameter::BetaEffective ||
+              settings::ifp_parameter == IFPParameter::Both) {
+            requests.emplace_back();
+            MPI_Isend(&send_delayed_groups[ifp_n_generation * index_local],
+              ifp_n_generation * static_cast<int>(n), MPI_INT, neighbor,
+              mpi::rank, mpi::intracomm, &requests.back());
+          }
+
+          // Send lifetimes
+          if (settings::ifp_parameter == IFPParameter::GenerationTime ||
+              settings::ifp_parameter == IFPParameter::Both) {
+            requests.emplace_back();
+            MPI_Isend(&send_lifetimes[ifp_n_generation * index_local],
+              ifp_n_generation * static_cast<int>(n), MPI_DOUBLE, neighbor,
+              mpi::rank, mpi::intracomm, &requests.back());
+          }
+        }
       }
 
       // Increment all indices
@@ -271,6 +404,10 @@ void synchronize_bank()
   start = simulation::work_index[mpi::rank];
   index_local = 0;
 
+  vector<int> recv_delayed_groups;
+  vector<double> recv_lifetimes;
+  vector<DeserializationInfo> deserialization_info;
+
   // Determine what process has the source sites that will need to be stored at
   // the beginning of this processor's source bank.
 
@@ -282,6 +419,18 @@ void synchronize_bank()
       upper_bound_index(bank_position, bank_position + mpi::n_procs, start);
   }
 
+  // Resize IFP receive buffers
+  if (settings::ifp && mpi::n_procs > 1) {
+    if (settings::ifp_parameter == IFPParameter::BetaEffective ||
+        settings::ifp_parameter == IFPParameter::Both) {
+      recv_delayed_groups.resize(ifp_n_generation * simulation::work_per_rank);
+    }
+    if (settings::ifp_parameter == IFPParameter::GenerationTime ||
+        settings::ifp_parameter == IFPParameter::Both) {
+      recv_lifetimes.resize(ifp_n_generation * simulation::work_per_rank);
+    }
+  }
+
   while (start < simulation::work_index[mpi::rank + 1]) {
     // Determine how many sites need to be received
     int64_t n;
@@ -301,13 +450,53 @@ void synchronize_bank()
       MPI_Irecv(&simulation::source_bank[index_local], static_cast<int>(n),
         mpi::source_site, neighbor, neighbor, mpi::intracomm, &requests.back());
 
+      if (settings::ifp) {
+
+        // Receive delayed groups
+        if (settings::ifp_parameter == IFPParameter::BetaEffective ||
+            settings::ifp_parameter == IFPParameter::Both) {
+          requests.emplace_back();
+          MPI_Irecv(&recv_delayed_groups[ifp_n_generation * index_local],
+            ifp_n_generation * static_cast<int>(n), MPI_INT, neighbor, neighbor,
+            mpi::intracomm, &requests.back());
+        }
+
+        // Receive lifetimes
+        if (settings::ifp_parameter == IFPParameter::GenerationTime ||
+            settings::ifp_parameter == IFPParameter::Both) {
+          requests.emplace_back();
+          MPI_Irecv(&recv_lifetimes[ifp_n_generation * index_local],
+            ifp_n_generation * static_cast<int>(n), MPI_DOUBLE, neighbor,
+            neighbor, mpi::intracomm, &requests.back());
+        }
+
+        // Deserialization info to reconstruct data later
+        DeserializationInfo info = {index_local, n};
+        deserialization_info.push_back(info);
+      }
+
     } else {
-      // If the source sites are on this procesor, we can simply copy them
+      // If the source sites are on this processor, we can simply copy them
       // from the temp_sites bank
 
       index_temp = start - bank_position[mpi::rank];
       std::copy(&temp_sites[index_temp], &temp_sites[index_temp + n],
         &simulation::source_bank[index_local]);
+
+      if (settings::ifp) {
+        if (settings::ifp_parameter == IFPParameter::BetaEffective ||
+            settings::ifp_parameter == IFPParameter::Both) {
+          std::copy(&temp_delayed_groups[index_temp],
+            &temp_delayed_groups[index_temp + n],
+            &simulation::ifp_source_delayed_group_bank[index_local]);
+        }
+        if (settings::ifp_parameter == IFPParameter::GenerationTime ||
+            settings::ifp_parameter == IFPParameter::Both) {
+          std::copy(&temp_lifetimes[index_temp],
+            &temp_lifetimes[index_temp + n],
+            &simulation::ifp_source_lifetime_bank[index_local]);
+        }
+      }
     }
 
     // Increment all indices
@@ -323,10 +512,65 @@ void synchronize_bank()
   int n_request = requests.size();
   MPI_Waitall(n_request, requests.data(), MPI_STATUSES_IGNORE);
 
+  if (settings::ifp) {
+
+    // Deserialize
+    int64_t n;
+    for (auto info : deserialization_info) {
+      index_local = info.index_local;
+      n = info.n;
+
+      // Store deserialized data in banks
+      for (int i = index_local; i < index_local + n; i++) {
+        if (settings::ifp_parameter == IFPParameter::BetaEffective ||
+            settings::ifp_parameter == IFPParameter::Both) {
+          vector<int> delayed_groups_received(
+            recv_delayed_groups.begin() + ifp_n_generation * i,
+            recv_delayed_groups.begin() + ifp_n_generation * (i + 1));
+          simulation::ifp_source_delayed_group_bank[i] =
+            delayed_groups_received;
+        }
+        if (settings::ifp_parameter == IFPParameter::GenerationTime ||
+            settings::ifp_parameter == IFPParameter::Both) {
+          vector<double> lifetimes_received(
+            recv_lifetimes.begin() + ifp_n_generation * i,
+            recv_lifetimes.begin() + ifp_n_generation * (i + 1));
+          simulation::ifp_source_lifetime_bank[i] = lifetimes_received;
+        }
+      }
+    }
+
+    // Clear IFP buffers
+    send_delayed_groups.clear();
+    send_lifetimes.clear();
+    recv_delayed_groups.clear();
+    recv_lifetimes.clear();
+
+    // Clear deserialization info
+    deserialization_info.clear();
+  }
+
 #else
   std::copy(temp_sites.data(), temp_sites.data() + settings::n_particles,
     simulation::source_bank.begin());
+  if (settings::ifp) {
+    if (settings::ifp_parameter == IFPParameter::BetaEffective ||
+        settings::ifp_parameter == IFPParameter::Both) {
+      std::copy(temp_delayed_groups.data(),
+        temp_delayed_groups.data() + settings::n_particles,
+        simulation::ifp_source_delayed_group_bank.begin());
+    }
+    if (settings::ifp_parameter == IFPParameter::GenerationTime ||
+        settings::ifp_parameter == IFPParameter::Both) {
+      std::copy(temp_lifetimes.data(),
+        temp_lifetimes.data() + settings::n_particles,
+        simulation::ifp_source_lifetime_bank.begin());
+    }
+  }
 #endif
+  temp_sites.clear();
+  temp_delayed_groups.clear();
+  temp_lifetimes.clear();
 
   simulation::time_bank_sendrecv.stop();
   simulation::time_bank.stop();
diff --git a/src/finalize.cpp b/src/finalize.cpp
index 2bde0e3387c..a574cc5b840 100644
--- a/src/finalize.cpp
+++ b/src/finalize.cpp
@@ -163,8 +163,9 @@ int openmc_finalize()
 
   // Free all MPI types
 #ifdef OPENMC_MPI
-  if (mpi::source_site != MPI_DATATYPE_NULL)
+  if (mpi::source_site != MPI_DATATYPE_NULL) {
     MPI_Type_free(&mpi::source_site);
+  }
 #endif
 
   return 0;
diff --git a/src/output.cpp b/src/output.cpp
index 5fdbea1304e..ae5d3658aeb 100644
--- a/src/output.cpp
+++ b/src/output.cpp
@@ -601,6 +601,9 @@ const std::unordered_map<int, const char*> score_names = {
   {SCORE_FISS_Q_RECOV, "Recoverable fission power"},
   {SCORE_CURRENT, "Current"},
   {SCORE_PULSE_HEIGHT, "pulse-height"},
+  {SCORE_IFP_TIME_NUM, "IFP lifetime numerator"},
+  {SCORE_IFP_BETA_NUM, "IFP delayed fraction numerator"},
+  {SCORE_IFP_DENOM, "IFP common denominator"},
 };
 
 //! Create an ASCII output file showing all tally results.
diff --git a/src/particle.cpp b/src/particle.cpp
index 64c50c9438f..ad9a10c1ca3 100644
--- a/src/particle.cpp
+++ b/src/particle.cpp
@@ -114,6 +114,7 @@ void Particle::from_source(const SourceSite* src)
   n_collision() = 0;
   fission() = false;
   zero_flux_derivs();
+  lifetime() = 0.0;
 
   // Copy attributes from source bank site
   type() = src->particle;
@@ -233,6 +234,7 @@ void Particle::event_advance()
     coord(j).r += distance * coord(j).u;
   }
   this->time() += distance / this->speed();
+  this->lifetime() += distance / this->speed();
 
   // Kill particle if its time exceeds the cutoff
   bool hit_time_boundary = false;
@@ -240,6 +242,7 @@ void Particle::event_advance()
   if (time() > time_cutoff) {
     double dt = time() - time_cutoff;
     time() = time_cutoff;
+    lifetime() = time_cutoff;
 
     double push_back_distance = speed() * dt;
     this->move_distance(-push_back_distance);
diff --git a/src/physics.cpp b/src/physics.cpp
index 69e74f2eafd..5909d8bda56 100644
--- a/src/physics.cpp
+++ b/src/physics.cpp
@@ -209,6 +209,59 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx)
     // Sample delayed group and angle/energy for fission reaction
     sample_fission_neutron(i_nuclide, rx, &site, p);
 
+    // Iterated Fission Probability (IFP) method
+    // Needs to be done after the delayed group is found
+
+    vector<int> updated_ifp_delayed_groups;
+    vector<double> updated_ifp_lifetimes;
+
+    if (settings::ifp) {
+
+      if (settings::ifp_parameter == IFPParameter::BetaEffective ||
+          settings::ifp_parameter == IFPParameter::Both) {
+
+        const auto& ifp_delayed_groups =
+          simulation::ifp_source_delayed_group_bank[p.current_work() - 1];
+        size_t idx = ifp_delayed_groups.size();
+
+        if (idx < settings::ifp_n_generation) {
+          updated_ifp_delayed_groups.resize(idx + 1);
+          for (size_t i = 0; i < idx; i++) {
+            updated_ifp_delayed_groups[i] = ifp_delayed_groups[i];
+          }
+          updated_ifp_delayed_groups[idx] = site.delayed_group;
+        } else if (idx == settings::ifp_n_generation) {
+          updated_ifp_delayed_groups.resize(idx);
+          for (size_t i = 0; i < idx - 1; i++) {
+            updated_ifp_delayed_groups[i] = ifp_delayed_groups[i + 1];
+          }
+          updated_ifp_delayed_groups[idx - 1] = site.delayed_group;
+        }
+      }
+
+      if (settings::ifp_parameter == IFPParameter::GenerationTime ||
+          settings::ifp_parameter == IFPParameter::Both) {
+
+        const auto& ifp_lifetimes =
+          simulation::ifp_source_lifetime_bank[p.current_work() - 1];
+        size_t idx = ifp_lifetimes.size();
+
+        if (idx < settings::ifp_n_generation) {
+          updated_ifp_lifetimes.resize(idx + 1);
+          for (size_t i = 0; i < idx; i++) {
+            updated_ifp_lifetimes[i] = ifp_lifetimes[i];
+          }
+          updated_ifp_lifetimes[idx] = p.lifetime();
+        } else if (idx == settings::ifp_n_generation) {
+          updated_ifp_lifetimes.resize(idx);
+          for (size_t i = 0; i < idx - 1; i++) {
+            updated_ifp_lifetimes[i] = ifp_lifetimes[i + 1];
+          }
+          updated_ifp_lifetimes[idx - 1] = p.lifetime();
+        }
+      }
+    }
+
     // Store fission site in bank
     if (use_fission_bank) {
       int64_t idx = simulation::fission_bank.thread_safe_append(site);
@@ -226,6 +279,21 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx)
         // Break out of loop as no more sites can be added to fission bank
         break;
       }
+      // If IFP, add the IFP information in the IFP banks using the same index
+      // as the one used to append the fission site to the fission bank.
+      // Multithreading protection is guaranteed by the index returned by the
+      // previous thread_safe_append call.
+      if (settings::ifp) {
+        if (settings::ifp_parameter == IFPParameter::BetaEffective ||
+            settings::ifp_parameter == IFPParameter::Both) {
+          simulation::ifp_fission_delayed_group_bank[idx] =
+            updated_ifp_delayed_groups;
+        }
+        if (settings::ifp_parameter == IFPParameter::GenerationTime ||
+            settings::ifp_parameter == IFPParameter::Both) {
+          simulation::ifp_fission_lifetime_bank[idx] = updated_ifp_lifetimes;
+        }
+      }
     } else {
       p.secondary_bank().push_back(site);
     }
diff --git a/src/reaction.cpp b/src/reaction.cpp
index 0643cb08bdd..2ac170dd5b7 100644
--- a/src/reaction.cpp
+++ b/src/reaction.cpp
@@ -172,6 +172,9 @@ std::unordered_map<int, std::string> REACTION_NAME_MAP {
   {SCORE_FISS_Q_PROMPT, "fission-q-prompt"},
   {SCORE_FISS_Q_RECOV, "fission-q-recoverable"},
   {SCORE_PULSE_HEIGHT, "pulse-height"},
+  {SCORE_IFP_TIME_NUM, "ifp-time-numerator"},
+  {SCORE_IFP_BETA_NUM, "ifp-beta-numerator"},
+  {SCORE_IFP_DENOM, "ifp-denominator"},
   // Normal ENDF-based reactions
   {TOTAL_XS, "(n,total)"},
   {ELASTIC, "(n,elastic)"},
diff --git a/src/settings.cpp b/src/settings.cpp
index ba451e219c5..90c11459df5 100644
--- a/src/settings.cpp
+++ b/src/settings.cpp
@@ -51,6 +51,7 @@ bool create_fission_neutrons {true};
 bool delayed_photon_scaling {true};
 bool entropy_on {false};
 bool event_based {false};
+bool ifp {false};
 bool legendre_to_tabular {true};
 bool material_cell_offsets {true};
 bool output_summary {true};
@@ -102,6 +103,8 @@ int max_particle_events {1000000};
 ElectronTreatment electron_treatment {ElectronTreatment::TTB};
 array<double, 4> energy_cutoff {0.0, 1000.0, 0.0, 0.0};
 array<double, 4> time_cutoff {INFTY, INFTY, INFTY, INFTY};
+int ifp_n_generation {0};
+IFPParameter ifp_parameter {IFPParameter::None};
 int legendre_to_tabular_points {C_NONE};
 int max_order {0};
 int n_log_bins {8000};
@@ -976,6 +979,55 @@ void read_settings_xml(pugi::xml_node root)
     temperature_range[1] = range.at(1);
   }
 
+  // Check for Iterated Fission Probability (IFP)
+  if (check_for_node(root, "iterated_fission_probability")) {
+
+    // IFP only works with eigenvalue calculations
+    if (run_mode == RunMode::EIGENVALUE) {
+      ifp = true;
+      // Get iterated_fission_probability write node
+      xml_node node_ifp = root.child("iterated_fission_probability");
+
+      // Parameter to compute
+      if (check_for_node(node_ifp, "parameter")) {
+        std::string temp_str =
+          get_node_value(node_ifp, "parameter", true, true);
+        if (temp_str == "beta_effective") {
+          ifp_parameter = IFPParameter::BetaEffective;
+        } else if (temp_str == "generation_time") {
+          ifp_parameter = IFPParameter::GenerationTime;
+        } else if (temp_str == "both") {
+          ifp_parameter = IFPParameter::Both;
+        } else {
+          fatal_error("Unrecognized parameter for IFP: " + temp_str);
+        }
+      } else {
+        fatal_error("<parameter> should be declared for ifp.");
+      }
+
+      // Number of generation
+      if (check_for_node(node_ifp, "n_generation")) {
+        ifp_n_generation = std::stoi(get_node_value(node_ifp, "n_generation"));
+        if (ifp_n_generation <= 0) {
+          fatal_error("<n_generation> must be greater than 0.");
+        }
+        // Avoid tallying 0 if IFP logs are not complete when active cycles
+        // start
+        if (ifp_n_generation > n_inactive) {
+          fatal_error("<n_generation> must be lower than or equal to the "
+                      "number of inactive cycles.");
+        }
+      } else {
+        fatal_error("<n_generation> must be specified as a subelement of "
+                    "<iterated_fission_probability>.");
+      }
+    } else {
+      fatal_error(
+        "Iterated Fission Probability can only be used in an eigenvalue "
+        "calculation.");
+    }
+  }
+
   // Check for tabular_legendre options
   if (check_for_node(root, "tabular_legendre")) {
     // Get pointer to tabular_legendre node
diff --git a/src/simulation.cpp b/src/simulation.cpp
index be03e48553c..81c714993ff 100644
--- a/src/simulation.cpp
+++ b/src/simulation.cpp
@@ -333,6 +333,23 @@ void allocate_banks()
 
     // Allocate fission bank
     init_fission_bank(3 * simulation::work_per_rank);
+
+    // Allocate IFP bank
+    if (settings::ifp) {
+      if (settings::ifp_parameter == IFPParameter::BetaEffective ||
+          settings::ifp_parameter == IFPParameter::Both) {
+        simulation::ifp_source_delayed_group_bank.resize(
+          simulation::work_per_rank);
+        simulation::ifp_fission_delayed_group_bank.resize(
+          3 * simulation::work_per_rank);
+      }
+      if (settings::ifp_parameter == IFPParameter::GenerationTime ||
+          settings::ifp_parameter == IFPParameter::Both) {
+        simulation::ifp_source_lifetime_bank.resize(simulation::work_per_rank);
+        simulation::ifp_fission_lifetime_bank.resize(
+          3 * simulation::work_per_rank);
+      }
+    }
   }
 
   if (settings::surf_source_write) {
diff --git a/src/source.cpp b/src/source.cpp
index 15fe8433ba5..4a60a9df422 100644
--- a/src/source.cpp
+++ b/src/source.cpp
@@ -594,6 +594,20 @@ void initialize_source()
 
     // sample external source distribution
     simulation::source_bank[i] = sample_external_source(&seed);
+
+    // Initialize IFP data
+    if (settings::ifp) {
+      if (settings::ifp_parameter == IFPParameter::BetaEffective ||
+          settings::ifp_parameter == IFPParameter::Both) {
+        vector<int> ifp_delayed_groups;
+        simulation::ifp_source_delayed_group_bank[i] = ifp_delayed_groups;
+      }
+      if (settings::ifp_parameter == IFPParameter::GenerationTime ||
+          settings::ifp_parameter == IFPParameter::Both) {
+        vector<double> ifp_lifetimes;
+        simulation::ifp_source_lifetime_bank[i] = ifp_lifetimes;
+      }
+    }
   }
 
   // Write out initial source
diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp
index 674987b8f13..73dadfbc6f9 100644
--- a/src/tallies/tally.cpp
+++ b/src/tallies/tally.cpp
@@ -584,7 +584,12 @@ void Tally::set_scores(const vector<std::string>& scores)
           }
         }
       }
+      break;
 
+    case SCORE_IFP_TIME_NUM:
+    case SCORE_IFP_BETA_NUM:
+    case SCORE_IFP_DENOM:
+      estimator_ = TallyEstimator::COLLISION;
       break;
     }
 
diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp
index e798161ec2f..a96cafb66a8 100644
--- a/src/tallies/tally_scoring.cpp
+++ b/src/tallies/tally_scoring.cpp
@@ -890,6 +890,61 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
         score_fission_q(p, score_bin, tally, flux, i_nuclide, atom_density);
       break;
 
+    case SCORE_IFP_TIME_NUM:
+      if (settings::ifp) {
+        if ((p.type() == Type::neutron) && (p.fission())) {
+          if (settings::ifp_parameter == IFPParameter::GenerationTime ||
+              settings::ifp_parameter == IFPParameter::Both) {
+            const auto& lifetimes =
+              simulation::ifp_source_lifetime_bank[p.current_work() - 1];
+            int n_generation = static_cast<int>(lifetimes.size());
+            if (n_generation == settings::ifp_n_generation) {
+              score = lifetimes[0] * p.wgt_last();
+            }
+          }
+        }
+      }
+      break;
+
+    case SCORE_IFP_BETA_NUM:
+      if (settings::ifp) {
+        if ((p.type() == Type::neutron) && (p.fission())) {
+          if (settings::ifp_parameter == IFPParameter::BetaEffective ||
+              settings::ifp_parameter == IFPParameter::Both) {
+            const auto& delayed_groups =
+              simulation::ifp_source_delayed_group_bank[p.current_work() - 1];
+            int n_generation = static_cast<int>(delayed_groups.size());
+            if (n_generation == settings::ifp_n_generation) {
+              if (delayed_groups[0] > 0) {
+                score = p.wgt_last();
+              }
+            }
+          }
+        }
+      }
+      break;
+
+    case SCORE_IFP_DENOM:
+      if (settings::ifp) {
+        if ((p.type() == Type::neutron) && (p.fission())) {
+          int n_generation;
+          if (settings::ifp_parameter == IFPParameter::BetaEffective ||
+              settings::ifp_parameter == IFPParameter::Both) {
+            n_generation = static_cast<int>(
+              simulation::ifp_source_delayed_group_bank[p.current_work() - 1]
+                .size());
+          } else {
+            n_generation = static_cast<int>(
+              simulation::ifp_source_lifetime_bank[p.current_work() - 1]
+                .size());
+          }
+          if (n_generation == settings::ifp_n_generation) {
+            score = p.wgt_last();
+          }
+        }
+      }
+      break;
+
     case N_2N:
     case N_3N:
     case N_4N:
diff --git a/tests/regression_tests/ifp/__init__.py b/tests/regression_tests/ifp/__init__.py
new file mode 100644
index 00000000000..e69de29bb2d
diff --git a/tests/regression_tests/ifp/inputs_true.dat b/tests/regression_tests/ifp/inputs_true.dat
new file mode 100644
index 00000000000..84f5e91acce
--- /dev/null
+++ b/tests/regression_tests/ifp/inputs_true.dat
@@ -0,0 +1,34 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material depletable="true" id="1" name="core">
+      <density units="g/cm3" value="16.0"/>
+      <nuclide ao="1.0" name="U235"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="-1" universe="1"/>
+    <surface boundary="vacuum" coeffs="0.0 0.0 0.0 10.0" id="1" type="sphere"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>1000</particles>
+    <batches>20</batches>
+    <inactive>5</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <space type="fission">
+        <parameters>-10.0 -10.0 -10.0 10.0 10.0 10.0</parameters>
+      </space>
+    </source>
+    <photon_transport>false</photon_transport>
+    <iterated_fission_probability>
+      <parameter>both</parameter>
+      <n_generation>5</n_generation>
+    </iterated_fission_probability>
+  </settings>
+  <tallies>
+    <tally id="1" name="ifp-scores">
+      <scores>ifp-time-numerator ifp-beta-numerator ifp-denominator</scores>
+    </tally>
+  </tallies>
+</model>
diff --git a/tests/regression_tests/ifp/results_true.dat b/tests/regression_tests/ifp/results_true.dat
new file mode 100644
index 00000000000..a74e2bd78bc
--- /dev/null
+++ b/tests/regression_tests/ifp/results_true.dat
@@ -0,0 +1,9 @@
+k-combined:
+1.007452E+00 5.705278E-03
+tally 1:
+8.996235E-08
+5.461421E-16
+4.800000E-02
+4.680000E-04
+1.512000E+01
+1.526063E+01
diff --git a/tests/regression_tests/ifp/test.py b/tests/regression_tests/ifp/test.py
new file mode 100644
index 00000000000..7567f69b3a5
--- /dev/null
+++ b/tests/regression_tests/ifp/test.py
@@ -0,0 +1,56 @@
+"""Test the Iterated Fission Probability (IFP) method to compute adjoint-weighted
+kinetics parameters using dedicated tallies."""
+
+import openmc
+import pytest
+
+from tests.testing_harness import PyAPITestHarness
+
+
+@pytest.fixture()
+def ifp_model():
+    openmc.reset_auto_ids()
+    model = openmc.model.Model()
+
+    # Material
+    material = openmc.Material(name="core")
+    material.add_nuclide("U235", 1.0)
+    material.set_density('g/cm3', 16.0)
+
+    # Geometry
+    radius = 10.0
+    sphere = openmc.Sphere(r=radius, boundary_type="vacuum")
+    cell = openmc.Cell(region=-sphere, fill=material)
+    model.geometry = openmc.Geometry([cell])
+
+    # Settings
+    settings = openmc.Settings()
+    settings.run_mode = 'eigenvalue'
+    settings.particles = 1000
+    settings.batches = 20
+    settings.inactive = 5
+    settings.photon_transport = False
+    settings.iterated_fission_probability = {"parameter": "both", "n_generation": 5}
+
+    bounds = [
+        -radius, -radius, -radius,
+        radius, radius, radius
+    ]
+    distribution = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
+    settings.source = openmc.source.IndependentSource(space=distribution)
+
+    model.settings = settings
+
+    # Tally IFP scores
+    tally = openmc.Tally(name="ifp-scores")
+    tally.scores = ["ifp-time-numerator", "ifp-beta-numerator", "ifp-denominator"]
+
+    tallies = openmc.Tallies([tally])
+    model.tallies = tallies
+
+    return model
+
+
+def test_iterated_fission_probability(ifp_model):
+    harness = PyAPITestHarness("statepoint.20.h5", model=ifp_model)
+    harness.main()
diff --git a/tests/unit_tests/test_ifp.py b/tests/unit_tests/test_ifp.py
new file mode 100644
index 00000000000..34199704c8c
--- /dev/null
+++ b/tests/unit_tests/test_ifp.py
@@ -0,0 +1,43 @@
+"""Test the Iterated Fission Probability (IFP) method to compute adjoint-weighted
+kinetics parameters using dedicated tallies."""
+
+import pytest
+import openmc
+
+
+def test_xml_serialization(run_in_tmpdir):
+    """Check that a simple use case can be written and read in XML."""
+    parameter = {"n_generation": 5}
+    settings = openmc.Settings()
+    settings.iterated_fission_probability = parameter
+    settings.export_to_xml()
+
+    read_settings = openmc.Settings.from_xml()
+    assert read_settings.iterated_fission_probability == parameter
+
+
+@pytest.fixture(scope="module")
+def geometry():
+    openmc.reset_auto_ids()
+    material = openmc.Material()
+    sphere = openmc.Sphere(r=1.0, boundary_type="vacuum")
+    cell = openmc.Cell(region=-sphere, fill=material)
+    return openmc.Geometry([cell])
+
+
+@pytest.mark.parametrize(
+    "settings, error",
+    [
+        ({"iterated_fission_probability": {"n_generation": 0}}, ValueError),
+        ({"run_mode": "fixed source", "iterated_fission_probability": {"n_generation": 5}}, RuntimeError),
+        ({"n_inactive": 6, "iterated_fission_probability": {"n_generation": 5}}, RuntimeError)
+    ],
+)
+def test_exceptions(settings, error, run_in_tmpdir, geometry):
+    """Test settings configuration that should return an error."""
+    with pytest.raises(error):
+        settings = openmc.Settings(**settings)
+        settings.particles = 100
+        settings.batches = 10
+        model = openmc.Model(geometry=geometry, settings=settings)
+        model.run()