diff --git a/CMakeLists.txt b/CMakeLists.txt
index 9c3e6a9..5e1e684 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -416,6 +416,7 @@ list(APPEND libopenmc_SOURCES
   src/particle_type.cpp
   src/photon.cpp
   src/physics.cpp
+  src/sensitivity_ce.cpp
   src/physics_common.cpp
   src/physics_mg.cpp
   src/plot.cpp
diff --git a/include/openmc/sensitivity_ce.h b/include/openmc/sensitivity_ce.h
new file mode 100644
index 0000000..187c138
--- /dev/null
+++ b/include/openmc/sensitivity_ce.h
@@ -0,0 +1,47 @@
+//! \file sensitivity_ce.h
+//! Continuous-energy k-eigenvalue sensitivity coefficients (local addition to OpenMC 0.16.0).
+//!
+//! Collision-history estimator with a fission-source importance F*(cell) (CLUTCH-type adjoint
+//! weighting). For one generation, the importance-weighted fission production is
+//!   Y = sum over collisions c of P_c,  P_c = (fission sites banked at c) * F*(cell of c),
+//! and to first order dk/k = dY/Y for a change of the data with the fission source held fixed.
+//! The derivative of each P_c with respect to ln(data) is obtained by differential-operator
+//! sampling along the history that leads to it:
+//!   - every track segment of length s in a material: -N_j sigma_{x,j}(E) s  for each nuclide j,
+//!     reaction x (loss: elastic, inelastic, n,2n, fission, capture, n,p, n,alpha);
+//!   - every collision: +1 for the nuclide and scattering reaction that was sampled (gain);
+//!   - the production at c itself: +1 for fission and nu-bar of the nuclide sampled at c;
+//!   - a history that starts from a fission site: +1 for chi of the parent nuclide in the birth
+//!     group (unconstrained; the normalisation constraint is applied in post-processing).
+//! S_x = sum_c P_c D_x(c) / sum_c P_c, with D_x(c) the running sum of these terms up to c.
+//!
+//! Enabled when the environment variable OPENMC_SENS names a configuration file:
+//!   groups G   <G+1 ascending boundaries in eV>
+//!   nuclides N <N OpenMC nuclide names>
+//!   fstar M    <M pairs: cell id, importance>      (optional; default importance 1)
+//! Results (active batches only) go to sensitivity.txt in the output directory.
+
+#ifndef OPENMC_SENSITIVITY_CE_H
+#define OPENMC_SENSITIVITY_CE_H
+
+namespace openmc {
+
+class Particle;
+struct SourceSite;
+
+namespace sens {
+
+extern bool on;
+
+void init();
+void begin_history(Particle& p, const SourceSite& site);
+void track(Particle& p, double distance);
+void production(Particle& p, int i_nuclide, int n_sites, double E_in);
+void scatter_gain(Particle& p, int i_nuclide, int mt, double E_in);
+void end_batch(bool active);
+void write();
+
+} // namespace sens
+} // namespace openmc
+
+#endif // OPENMC_SENSITIVITY_CE_H
diff --git a/src/particle.cpp b/src/particle.cpp
index 30ccb90..5dda237 100644
--- a/src/particle.cpp
+++ b/src/particle.cpp
@@ -1,4 +1,5 @@
 #include "openmc/particle.h"
+#include "openmc/sensitivity_ce.h"
 
 #include <algorithm> // copy, min
 #include <cmath>     // log, abs
@@ -295,6 +296,10 @@ void Particle::event_advance()
 
   // Advance particle in space and time
   this->move_distance(distance);
+
+  // CE sensitivity: loss terms along the track
+  if (sens::on && type().is_neutron() && material() != MATERIAL_VOID)
+    sens::track(*this, distance);
   double dt = distance / speed;
   this->time() += dt;
   this->lifetime() += dt;
diff --git a/src/physics.cpp b/src/physics.cpp
index 4bf459b..ba63f32 100644
--- a/src/physics.cpp
+++ b/src/physics.cpp
@@ -26,6 +26,7 @@
 #include "openmc/string_utils.h"
 #include "openmc/tallies/tally.h"
 #include "openmc/thermal.h"
+#include "openmc/sensitivity_ce.h"
 #include "openmc/weight_windows.h"
 
 #include <fmt/core.h>
@@ -104,6 +105,9 @@ void collision(Particle& p)
 
 void sample_neutron_reaction(Particle& p)
 {
+  // Incoming energy, for the CE sensitivity terms (sensitivity_ce)
+  double E_in_sens = p.E();
+
   // Sample a nuclide within the material
   int i_nuclide = sample_nuclide(p);
 
@@ -159,6 +163,8 @@ void sample_neutron_reaction(Particle& p)
     ncrystal_mat.scatter(p);
   } else {
     scatter(p, i_nuclide);
+    if (sens::on)
+      sens::scatter_gain(p, i_nuclide, p.event_mt(), E_in_sens);
   }
 
   // Advance URR seed stream 'N' times after energy changes
@@ -216,6 +222,7 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx)
     site.time = p.time();
     site.wgt = 1. / weight;
     site.surf_id = 0;
+    site.parent_nuclide = i_nuclide;
 
     // Sample delayed group and angle/energy for fission reaction
     sample_fission_neutron(i_nuclide, rx, &site, p);
@@ -284,6 +291,10 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx)
   // bank was not found to be full then these values are already equivalent.
   nu = n_sites_stored;
 
+  // CE sensitivity: importance-weighted production at this collision
+  if (sens::on && use_fission_bank)
+    sens::production(p, i_nuclide, nu, p.E());
+
   // Store the total weight banked for analog fission tallies
   p.n_bank() = nu;
   p.wgt_bank() = nu / weight;
diff --git a/src/sensitivity_ce.cpp b/src/sensitivity_ce.cpp
new file mode 100644
index 0000000..91ec995
--- /dev/null
+++ b/src/sensitivity_ce.cpp
@@ -0,0 +1,315 @@
+//! \file sensitivity_ce.cpp
+//! Continuous-energy k-eigenvalue sensitivity coefficients; see sensitivity_ce.h.
+
+#include "openmc/sensitivity_ce.h"
+
+#include <algorithm>
+#include <cmath>
+#include <cstdlib>
+#include <fstream>
+#include <memory>
+#include <mutex>
+#include <sstream>
+#include <string>
+#include <unordered_map>
+
+#include "openmc/cell.h"
+#include "openmc/error.h"
+#include "openmc/material.h"
+#include "openmc/nuclide.h"
+#include "openmc/particle.h"
+#include "openmc/settings.h"
+#include "openmc/simulation.h"
+
+namespace openmc {
+namespace sens {
+
+bool on {false};
+
+namespace {
+
+enum Rx { EL, INEL, N2N, FIS, CAP, NP, NA, NU, CHI, NRX };
+const int RX_MT[NRX] = {2, 4, 16, 18, 102, 103, 107, 452, 1018};
+
+bool initialized {false};
+int G {0};
+vector<double> bounds;                // ascending, G+1
+vector<std::string> names;            // sensitivity nuclides
+vector<int> slot_of;                  // global nuclide index -> slot (or -1)
+vector<double> fstar;                 // cell index -> fission-source importance
+size_t nbins {0};
+
+struct Buf {
+  vector<double> R;                   // running derivative of the history
+  vector<char> flag;
+  vector<int> touched;
+  vector<double> num;                 // sum_c P_c D(c), this batch
+  vector<double> births;              // fission-born source neutrons per (slot, group)
+  double den {0.0};                   // sum_c P_c, this batch
+  void resize()
+  {
+    R.assign(nbins, 0.0);
+    flag.assign(nbins, 0);
+    num.assign(nbins, 0.0);
+    births.assign(names.size() * G, 0.0);
+  }
+  void add(size_t i, double v)
+  {
+    if (!flag[i]) {
+      flag[i] = 1;
+      touched.push_back(static_cast<int>(i));
+    }
+    R[i] += v;
+  }
+  void reset()
+  {
+    for (int i : touched) {
+      R[i] = 0.0;
+      flag[i] = 0;
+    }
+    touched.clear();
+  }
+};
+
+std::mutex reg_mutex;
+vector<std::unique_ptr<Buf>> buffers;
+thread_local Buf* tbuf {nullptr};
+
+// accumulated over active batches
+vector<double> tot_num, sum_s, sum_s2, tot_births;
+double tot_den {0.0};
+int n_real {0};
+
+Buf& buf()
+{
+  if (!tbuf) {
+    std::lock_guard<std::mutex> lock(reg_mutex);
+    buffers.push_back(std::make_unique<Buf>());
+    tbuf = buffers.back().get();
+    tbuf->resize();
+  }
+  return *tbuf;
+}
+
+inline size_t bin(int slot, int rx, int g)
+{
+  return (static_cast<size_t>(slot) * NRX + rx) * G + g;
+}
+
+inline int group(double E)
+{
+  if (E < bounds.front() || E >= bounds.back())
+    return -1;
+  auto it = std::upper_bound(bounds.begin(), bounds.end(), E);
+  return static_cast<int>(it - bounds.begin()) - 1;
+}
+
+} // namespace
+
+void init()
+{
+  if (initialized)
+    return;
+  initialized = true;
+  const char* path = std::getenv("OPENMC_SENS");
+  if (!path || settings::run_mode != RunMode::EIGENVALUE)
+    return;
+  std::ifstream in(path);
+  if (!in)
+    fatal_error(fmt::format("OPENMC_SENS: cannot open {}", path));
+  std::string key;
+  std::unordered_map<int, double> fs;
+  while (in >> key) {
+    if (key == "groups") {
+      in >> G;
+      bounds.resize(G + 1);
+      for (auto& b : bounds)
+        in >> b;
+    } else if (key == "nuclides") {
+      int n;
+      in >> n;
+      names.resize(n);
+      for (auto& s : names)
+        in >> s;
+    } else if (key == "fstar") {
+      int n;
+      in >> n;
+      for (int i = 0; i < n; ++i) {
+        int id;
+        double v;
+        in >> id >> v;
+        fs[id] = v;
+      }
+    } else {
+      fatal_error(fmt::format("OPENMC_SENS: unknown keyword {}", key));
+    }
+  }
+  if (G <= 0 || names.empty())
+    fatal_error("OPENMC_SENS: groups and nuclides are required");
+  slot_of.assign(data::nuclides.size(), -1);
+  for (int s = 0; s < static_cast<int>(names.size()); ++s) {
+    auto it = data::nuclide_map.find(names[s]);
+    if (it == data::nuclide_map.end())
+      fatal_error(fmt::format("OPENMC_SENS: nuclide {} is not in the problem", names[s]));
+    slot_of[it->second] = s;
+  }
+  fstar.assign(model::cells.size(), 1.0);
+  for (auto& [id, v] : fs) {
+    auto it = model::cell_map.find(id);
+    if (it == model::cell_map.end())
+      fatal_error(fmt::format("OPENMC_SENS: cell {} not found", id));
+    fstar[it->second] = v;
+  }
+  nbins = names.size() * NRX * G;
+  tot_num.assign(nbins, 0.0);
+  sum_s.assign(nbins, 0.0);
+  sum_s2.assign(nbins, 0.0);
+  tot_births.assign(names.size() * G, 0.0);
+  simulation::need_depletion_rx = true; // (n,gamma), (n,p), (n,a), (n,2n) per nuclide
+  on = true;
+  write_message(fmt::format("Sensitivity (CE, collision history): {} nuclides, {} groups, "
+                            "{} cells with importance", names.size(), G, fs.size()),
+    5);
+}
+
+void begin_history(Particle& p, const SourceSite& site)
+{
+  Buf& b = buf();
+  b.reset();
+  if (site.parent_nuclide < 0)
+    return;
+  int s = slot_of[site.parent_nuclide];
+  int g = group(site.E);
+  if (s < 0 || g < 0)
+    return;
+  b.add(bin(s, CHI, g), 1.0);
+  b.births[static_cast<size_t>(s) * G + g] += 1.0;
+}
+
+void track(Particle& p, double distance)
+{
+  int g = group(p.E());
+  if (g < 0)
+    return;
+  Buf& b = buf();
+  const auto& mat = *model::materials[p.material()];
+  for (size_t i = 0; i < mat.nuclide_.size(); ++i) {
+    int j = mat.nuclide_[i];
+    int s = slot_of[j];
+    if (s < 0)
+      continue;
+    double ns = mat.atom_density(i, p.density_mult()) * distance;
+    // OpenMC evaluates the elastic cross section lazily
+    if (p.neutron_xs(j).elastic == CACHE_INVALID)
+      data::nuclides[j]->calculate_elastic_xs(p);
+    const auto& m = p.neutron_xs(j);
+    double np = m.reaction[1], na = m.reaction[2], n2n = m.reaction[3];
+    double cap = m.absorption - m.fission - np - na;
+    double inel = m.total - m.elastic - m.absorption - n2n;
+    double xs[7] = {m.elastic, inel, n2n, m.fission, cap, np, na};
+    for (int r = 0; r < 7; ++r) {
+      if (xs[r] != 0.0)
+        b.add(bin(s, r, g), -ns * xs[r]);
+    }
+  }
+}
+
+void production(Particle& p, int i_nuclide, int n_sites, double E_in)
+{
+  if (n_sites <= 0)
+    return;
+  Buf& b = buf();
+  double P = n_sites * fstar[p.lowest_coord().cell()];
+  b.den += P;
+  for (int i : b.touched)
+    b.num[i] += P * b.R[i];
+  int s = slot_of[i_nuclide];
+  int g = group(E_in);
+  if (s >= 0 && g >= 0) {
+    b.num[bin(s, FIS, g)] += P;
+    b.num[bin(s, NU, g)] += P;
+  }
+}
+
+void scatter_gain(Particle& p, int i_nuclide, int mt, double E_in)
+{
+  int s = slot_of[i_nuclide];
+  int g = group(E_in);
+  if (s < 0 || g < 0)
+    return;
+  int r = (mt == 2) ? EL : (mt == 16 ? N2N : INEL);
+  buf().add(bin(s, r, g), 1.0);
+}
+
+void end_batch(bool active)
+{
+  vector<double> num(nbins, 0.0), births(names.size() * G, 0.0);
+  double den = 0.0;
+  for (auto& b : buffers) {
+    for (size_t i = 0; i < nbins; ++i) {
+      num[i] += b->num[i];
+      b->num[i] = 0.0;
+    }
+    for (size_t i = 0; i < births.size(); ++i) {
+      births[i] += b->births[i];
+      b->births[i] = 0.0;
+    }
+    den += b->den;
+    b->den = 0.0;
+  }
+  if (!active || den <= 0.0)
+    return;
+  for (size_t i = 0; i < nbins; ++i) {
+    double s = num[i] / den;
+    tot_num[i] += num[i];
+    sum_s[i] += s;
+    sum_s2[i] += s * s;
+  }
+  for (size_t i = 0; i < births.size(); ++i)
+    tot_births[i] += births[i];
+  tot_den += den;
+  ++n_real;
+}
+
+void write()
+{
+  if (!on)
+    return;
+  std::ofstream out(settings::path_output + "sensitivity.txt");
+  out.precision(9);
+  out << "# continuous-energy k sensitivity (collision history, fission-source importance)\n";
+  out << "realizations " << n_real << "\n";
+  out << "production " << tot_den << "\n";
+  out << "groups " << G << "\n";
+  for (double e : bounds)
+    out << e << " ";
+  out << "\nnuclides " << names.size() << "\n";
+  for (auto& s : names)
+    out << s << " ";
+  out << "\n";
+  for (size_t s = 0; s < names.size(); ++s) {
+    out << "births " << names[s];
+    for (int g = 0; g < G; ++g)
+      out << " " << tot_births[s * G + g];
+    out << "\n";
+    for (int r = 0; r < NRX; ++r) {
+      out << "S " << names[s] << " " << RX_MT[r];
+      for (int g = 0; g < G; ++g)
+        out << " " << (tot_den > 0 ? tot_num[bin(s, r, g)] / tot_den : 0.0);
+      out << "\nU " << names[s] << " " << RX_MT[r];
+      for (int g = 0; g < G; ++g) {
+        size_t i = bin(s, r, g);
+        double var = 0.0;
+        if (n_real > 1) {
+          double m = sum_s[i] / n_real;
+          var = std::max(0.0, (sum_s2[i] / n_real - m * m) / (n_real - 1));
+        }
+        out << " " << std::sqrt(var);
+      }
+      out << "\n";
+    }
+  }
+}
+
+} // namespace sens
+} // namespace openmc
diff --git a/src/simulation.cpp b/src/simulation.cpp
index 03f40a7..80cb730 100644
--- a/src/simulation.cpp
+++ b/src/simulation.cpp
@@ -1,4 +1,5 @@
 #include "openmc/simulation.h"
+#include "openmc/sensitivity_ce.h"
 
 #include "openmc/bank.h"
 #include "openmc/capi.h"
@@ -182,6 +183,9 @@ int openmc_simulation_finalize()
   simulation::time_active.stop();
   simulation::time_finalize.start();
 
+  // CE sensitivity results
+  sens::write();
+
   // Clear material nuclide mapping
   for (auto& mat : model::materials) {
     mat->mat_nuclide_index_.clear();
@@ -477,6 +481,9 @@ void allocate_banks()
 
 void initialize_batch()
 {
+  // CE sensitivity (enabled by OPENMC_SENS; reads its configuration once)
+  sens::init();
+
   // Increment current batch
   ++simulation::current_batch;
   if (settings::run_mode == RunMode::FIXED_SOURCE) {
@@ -520,6 +527,10 @@ void initialize_batch()
 
 void finalize_batch()
 {
+  // CE sensitivity: collect this batch
+  if (sens::on)
+    sens::end_batch(simulation::current_batch > settings::n_inactive);
+
   // Reduce tallies onto master process and accumulate
   simulation::time_tallies.start();
   accumulate_tallies();
@@ -694,6 +705,8 @@ void sample_source_particle(Particle& p, int64_t index_source)
   // Sample a particle from the source bank
   if (settings::run_mode == RunMode::EIGENVALUE) {
     p.from_source(&simulation::source_bank[index_source - 1]);
+    if (sens::on)
+      sens::begin_history(p, simulation::source_bank[index_source - 1]);
   } else if (settings::run_mode == RunMode::FIXED_SOURCE) {
     // initialize random number seed
     int64_t id = compute_transport_seed(compute_particle_id(index_source));