diff --git a/docs/source/methods/random_ray.rst b/docs/source/methods/random_ray.rst
index 3d98747e4af..c827c351dd0 100644
--- a/docs/source/methods/random_ray.rst
+++ b/docs/source/methods/random_ray.rst
@@ -542,7 +542,7 @@ in that cell for the iteration from Equation :eq:`phi_naive` to:
 .. math::
     :label: phi_missed_one
 
-    \phi_{i,g,n}^{missed} = \frac{Q_{i,g,n} }{\Sigma_{t,i,g}} 
+    \phi_{i,g,n}^{missed} = \frac{Q_{i,g,n} }{\Sigma_{t,i,g}}
 
 as the streaming operator has gone to zero. While this is obviously innacurate
 as it ignores transport, for most problems where the region is only occasionally
@@ -1060,6 +1060,49 @@ random ray and Monte Carlo, however.
   develop the scattering source by way of inactive batches before beginning
   active batches.
 
+------------------------
+Adjoint Flux Solver Mode
+------------------------
+
+The random ray solver in OpenMC can also be used to solve for the adjoint flux,
+:math:`\psi^{\dagger}`. In combination with the regular (forward) flux solution,
+the adjoint flux is useful for perturbation methods as well as for computing
+weight windows for subsequent Monte Carlo simulations. The adjoint flux can be
+thought of as the "backwards" flux, representing the flux where a particle is
+born at an absoprtion point (and typical absorption energy), and then undergoes
+transport with a transposed scattering matrix. That is, instead of sampling a
+particle and seeing where it might go as in a standard forward solve, we will
+sample an absorption location and see where the particle that was absorbed there
+might have come from. Notably, for typical neutron absorption at low energy
+levels, this means that adjoint flux particles are typically sampled at a low
+energy and then upscatter (via a transposed scattering matrix) over their
+lifetimes.
+
+In OpenMC, the random ray adjoint solver is implemented simply by transposing
+the scattering matrix, swapping :math:`\nu\Sigma_f` and :math:`\chi`, and then
+running a normal transport solve. When no external fixed source is present, no
+additional changes are needed in the transport process. However, if an external
+fixed forward source is present in the simulation problem, then an additional
+step is taken to compute the accompanying fixed adjoint source. In OpenMC, the
+adjoint flux does *not* represent a response function for a particular detector
+region. Rather, the adjoint flux is the global response, making it appropriate
+for use with weight window generation schemes for global variance reduction.
+Thus, if using a fixed source, the external source for the adjoint mode is
+simply computed as being :math:`1 / \phi`, where :math:`\phi` is the forward
+scalar flux that results from a normal forward solve (which OpenMC will run
+first automatically when in adjoint mode). The adjoint external source will be
+computed for each source region in the simulation mesh, independent of any
+tallies. The adjoint external source is always flat, even when a linear
+scattering and fission source shape is used. When in adjoint mode, all reported
+results (e.g., tallies, eigenvalues, etc.) are derived from the adjoint flux,
+even when the physical meaning is not necessarily obvious. These values are
+still reported, though we emphasize that the primary use case for adjoint mode
+is for producing adjoint flux tallies to support subsequent perturbation studies
+and weight window generation.
+
+Note that the adjoint :math:`k_{eff}` is statistically the same as the forward
+:math:`k_{eff}`, despite the flux distributions taking different shapes.
+
 ---------------------------
 Fundamental Sources of Bias
 ---------------------------
diff --git a/docs/source/usersguide/random_ray.rst b/docs/source/usersguide/random_ray.rst
index 5ca0ab6bed9..2b9cf67240b 100644
--- a/docs/source/usersguide/random_ray.rst
+++ b/docs/source/usersguide/random_ray.rst
@@ -558,7 +558,7 @@ following methods are currently available in OpenMC:
      - Cons
    * - ``simulation_averaged``
      - Accumulates total active ray lengths in each FSR over all iterations,
-       improving the estimate of the volume in each cell each iteration. 
+       improving the estimate of the volume in each cell each iteration.
      - * Virtually unbiased after several iterations
        * Asymptotically approaches the true analytical volume
        * Typically most efficient in terms of speed vs. accuracy
@@ -593,6 +593,33 @@ estimator, the following code would be used:
 
     settings.random_ray['volume_estimator'] = 'naive'
 
+-----------------
+Adjoint Flux Mode
+-----------------
+
+The adjoint flux random ray solver mode can be enabled as:
+entire
+::
+
+    settings.random_ray['adjoint'] = True
+
+When enabled, OpenMC will first run a forward transport simulation followed by
+an adjoint transport simulation. The purpose of the forward solve is to compute
+the adjoint external source when an external source is present in the
+simulation. Simulation settings (e.g., number of rays, batches, etc.) will be
+identical for both simulations. At the conclusion of the run, all results (e.g.,
+tallies, plots, etc.) will be derived from the adjoint flux rather than the
+forward flux but are not labeled any differently. The initial forward flux
+solution will not be stored or available in the final statepoint file. Those
+wishing to do analysis requiring both the forward and adjoint solutions will
+need to run two separate simulations and load both statepoint files.
+
+.. note::
+    When adjoint mode is selected, OpenMC will always perform a full forward
+    solve and then run a full adjoint solve immediately afterwards. Statepoint
+    and tally results will be derived from the adjoint flux, but will not be
+    labeled any differently.
+
 ---------------------------------------
 Putting it All Together: Example Inputs
 ---------------------------------------
diff --git a/include/openmc/mgxs.h b/include/openmc/mgxs.h
index 3fb0608104f..9b1602f299a 100644
--- a/include/openmc/mgxs.h
+++ b/include/openmc/mgxs.h
@@ -86,8 +86,10 @@ class Mgxs {
   bool fissionable; // Is this fissionable
   bool is_isotropic {
     true}; // used to skip search for angle indices if isotropic
+  bool exists_in_model {true}; // Is this present in model
 
   Mgxs() = default;
+  Mgxs(bool exists) : exists_in_model(exists) {}
 
   //! \brief Constructor that loads the Mgxs object from the HDF5 file
   //!
diff --git a/include/openmc/random_ray/flat_source_domain.h b/include/openmc/random_ray/flat_source_domain.h
index 5c50f7fb0a0..b5b5db7062e 100644
--- a/include/openmc/random_ray/flat_source_domain.h
+++ b/include/openmc/random_ray/flat_source_domain.h
@@ -111,13 +111,17 @@ class FlatSourceDomain {
   virtual void all_reduce_replicated_source_regions();
   void convert_external_sources();
   void count_external_source_regions();
+  void set_adjoint_sources(const vector<double>& forward_flux);
   virtual void flux_swap();
   virtual double evaluate_flux_at_point(Position r, int64_t sr, int g) const;
   double compute_fixed_source_normalization_factor() const;
+  void flatten_xs();
+  void transpose_scattering_matrix();
 
   //----------------------------------------------------------------------------
   // Static Data members
   static bool volume_normalized_flux_tallies_;
+  static bool adjoint_; // If the user wants outputs based on the adjoint flux
 
   //----------------------------------------------------------------------------
   // Static data members
@@ -150,6 +154,19 @@ class FlatSourceDomain {
   vector<float> source_;
   vector<float> external_source_;
   vector<bool> external_source_present_;
+  vector<double> scalar_flux_final_;
+
+  // 2D arrays stored in 1D representing values for all materials x energy
+  // groups
+  int n_materials_;
+  vector<double> sigma_t_;
+  vector<double> nu_sigma_f_;
+  vector<double> sigma_f_;
+  vector<double> chi_;
+
+  // 3D arrays stored in 1D representing values for all materials x energy
+  // groups x energy groups
+  vector<double> sigma_s_;
 
 protected:
   //----------------------------------------------------------------------------
@@ -190,10 +207,6 @@ class FlatSourceDomain {
   vector<int> material_;
   vector<double> volume_naive_;
 
-  // 2D arrays stored in 1D representing values for all source regions x energy
-  // groups
-  vector<float> scalar_flux_final_;
-
   // Volumes for each tally and bin/score combination. This intermediate data
   // structure is used when tallying quantities that must be normalized by
   // volume (i.e., flux). The vector is index by tally index, while the inner 2D
diff --git a/include/openmc/random_ray/random_ray_simulation.h b/include/openmc/random_ray/random_ray_simulation.h
index 55bac6905c6..01f12bffb17 100644
--- a/include/openmc/random_ray/random_ray_simulation.h
+++ b/include/openmc/random_ray/random_ray_simulation.h
@@ -20,6 +20,8 @@ class RandomRaySimulation {
   //----------------------------------------------------------------------------
   // Methods
   void compute_segment_correction_factors();
+  void prepare_fixed_sources();
+  void prepare_fixed_sources_adjoint(vector<double>& forward_flux);
   void simulate();
   void reduce_simulation_statistics();
   void output_simulation_results() const;
@@ -30,8 +32,13 @@ class RandomRaySimulation {
     int64_t n_external_source_regions) const;
 
   //----------------------------------------------------------------------------
-  // Data members
+  // Accessors
+  FlatSourceDomain* domain() const { return domain_.get(); }
+
 private:
+  //----------------------------------------------------------------------------
+  // Data members
+
   // Contains all flat source region data
   unique_ptr<FlatSourceDomain> domain_;
 
diff --git a/openmc/deplete/openmc_operator.py b/openmc/deplete/openmc_operator.py
index 5837cd2c187..049e0dd37ee 100644
--- a/openmc/deplete/openmc_operator.py
+++ b/openmc/deplete/openmc_operator.py
@@ -208,7 +208,7 @@ def _get_burnable_mats(self) -> tuple[list[str], dict[str, float], list[str]]:
                 if nuclide in self.nuclides_with_data or self._decay_nucs:
                     model_nuclides.add(nuclide)
                 else:
-                    msg = (f"Nuclilde {nuclide} in material {mat.id} is not "
+                    msg = (f"Nuclide {nuclide} in material {mat.id} is not "
                            "present in the depletion chain and has no cross "
                            "section data.")
                     warn(msg)
diff --git a/openmc/settings.py b/openmc/settings.py
index 0a78fb564f8..a350de72ecc 100644
--- a/openmc/settings.py
+++ b/openmc/settings.py
@@ -170,6 +170,9 @@ class Settings:
             cm/cm^3. When disabled, flux tallies will be reported in units
             of cm (i.e., total distance traveled by neutrons in the spatial
             tally region).
+        :adjoint:
+            Whether to run the random ray solver in adjoint mode (bool). The
+            default is 'False'.
 
         .. versionadded:: 0.15.0
     resonance_scattering : dict
@@ -1113,6 +1116,8 @@ def random_ray(self, random_ray: dict):
                                ('flat', 'linear', 'linear_xy'))
             elif key == 'volume_normalized_flux_tallies':
                 cv.check_type('volume normalized flux tallies', random_ray[key], bool)
+            elif key == 'adjoint':
+                cv.check_type('adjoint', random_ray[key], bool)
             else:
                 raise ValueError(f'Unable to set random ray to "{key}" which is '
                                  'unsupported by OpenMC')
@@ -1916,6 +1921,10 @@ def _random_ray_from_xml_element(self, root):
                     self.random_ray['volume_normalized_flux_tallies'] = (
                         child.text in ('true', '1')
                     )
+                elif child.tag == 'adjoint':
+                    self.random_ray['adjoint'] = (
+                        child.text in ('true', '1')
+                    )
 
     def to_xml_element(self, mesh_memo=None):
         """Create a 'settings' element to be written to an XML file.
diff --git a/src/mgxs_interface.cpp b/src/mgxs_interface.cpp
index 0febc612bf0..ed734d401ea 100644
--- a/src/mgxs_interface.cpp
+++ b/src/mgxs_interface.cpp
@@ -146,7 +146,7 @@ void MgxsInterface::create_macro_xs()
         num_energy_groups_, num_delayed_groups_);
     } else {
       // Preserve the ordering of materials by including a blank entry
-      macro_xs_.emplace_back();
+      macro_xs_.emplace_back(false);
     }
   }
 }
diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp
index 62768b55f0f..19913c03c16 100644
--- a/src/random_ray/flat_source_domain.cpp
+++ b/src/random_ray/flat_source_domain.cpp
@@ -27,6 +27,7 @@ namespace openmc {
 RandomRayVolumeEstimator FlatSourceDomain::volume_estimator_ {
   RandomRayVolumeEstimator::HYBRID};
 bool FlatSourceDomain::volume_normalized_flux_tallies_ {false};
+bool FlatSourceDomain::adjoint_ {false};
 
 FlatSourceDomain::FlatSourceDomain() : negroups_(data::mg.num_energy_groups_)
 {
@@ -134,31 +135,23 @@ void FlatSourceDomain::update_neutron_source(double k_eff)
 
   double inverse_k_eff = 1.0 / k_eff;
 
-  // Temperature and angle indices, if using multiple temperature
-  // data sets and/or anisotropic data sets.
-  // TODO: Currently assumes we are only using single temp/single angle data.
-  const int t = 0;
-  const int a = 0;
-
   // Add scattering source
 #pragma omp parallel for
   for (int sr = 0; sr < n_source_regions_; sr++) {
     int material = material_[sr];
 
-    for (int e_out = 0; e_out < negroups_; e_out++) {
-      double sigma_t = data::mg.macro_xs_[material].get_xs(
-        MgxsType::TOTAL, e_out, nullptr, nullptr, nullptr, t, a);
-      double scatter_source = 0.0f;
+    for (int g_out = 0; g_out < negroups_; g_out++) {
+      double sigma_t = sigma_t_[material * negroups_ + g_out];
+      double scatter_source = 0.0;
 
-      for (int e_in = 0; e_in < negroups_; e_in++) {
-        double scalar_flux = scalar_flux_old_[sr * negroups_ + e_in];
-
-        double sigma_s = data::mg.macro_xs_[material].get_xs(
-          MgxsType::NU_SCATTER, e_in, &e_out, nullptr, nullptr, t, a);
+      for (int g_in = 0; g_in < negroups_; g_in++) {
+        double scalar_flux = scalar_flux_old_[sr * negroups_ + g_in];
+        double sigma_s =
+          sigma_s_[material * negroups_ * negroups_ + g_out * negroups_ + g_in];
         scatter_source += sigma_s * scalar_flux;
       }
 
-      source_[sr * negroups_ + e_out] = scatter_source / sigma_t;
+      source_[sr * negroups_ + g_out] = scatter_source / sigma_t;
     }
   }
 
@@ -167,20 +160,17 @@ void FlatSourceDomain::update_neutron_source(double k_eff)
   for (int sr = 0; sr < n_source_regions_; sr++) {
     int material = material_[sr];
 
-    for (int e_out = 0; e_out < negroups_; e_out++) {
-      double sigma_t = data::mg.macro_xs_[material].get_xs(
-        MgxsType::TOTAL, e_out, nullptr, nullptr, nullptr, t, a);
-      double fission_source = 0.0f;
-
-      for (int e_in = 0; e_in < negroups_; e_in++) {
-        double scalar_flux = scalar_flux_old_[sr * negroups_ + e_in];
-        double nu_sigma_f = data::mg.macro_xs_[material].get_xs(
-          MgxsType::NU_FISSION, e_in, nullptr, nullptr, nullptr, t, a);
-        double chi = data::mg.macro_xs_[material].get_xs(
-          MgxsType::CHI_PROMPT, e_in, &e_out, nullptr, nullptr, t, a);
+    for (int g_out = 0; g_out < negroups_; g_out++) {
+      double sigma_t = sigma_t_[material * negroups_ + g_out];
+      double fission_source = 0.0;
+
+      for (int g_in = 0; g_in < negroups_; g_in++) {
+        double scalar_flux = scalar_flux_old_[sr * negroups_ + g_in];
+        double nu_sigma_f = nu_sigma_f_[material * negroups_ + g_in];
+        double chi = chi_[material * negroups_ + g_out];
         fission_source += nu_sigma_f * scalar_flux * chi;
       }
-      source_[sr * negroups_ + e_out] +=
+      source_[sr * negroups_ + g_out] +=
         fission_source * inverse_k_eff / sigma_t;
     }
   }
@@ -188,7 +178,7 @@ void FlatSourceDomain::update_neutron_source(double k_eff)
   // Add external source if in fixed source mode
   if (settings::run_mode == RunMode::FIXED_SOURCE) {
 #pragma omp parallel for
-    for (int se = 0; se < n_source_elements_; se++) {
+    for (int64_t se = 0; se < n_source_elements_; se++) {
       source_[se] += external_source_[se];
     }
   }
@@ -206,8 +196,8 @@ void FlatSourceDomain::normalize_scalar_flux_and_volumes(
 
 // Normalize scalar flux to total distance travelled by all rays this iteration
 #pragma omp parallel for
-  for (int64_t e = 0; e < scalar_flux_new_.size(); e++) {
-    scalar_flux_new_[e] *= normalization_factor;
+  for (int64_t se = 0; se < scalar_flux_new_.size(); se++) {
+    scalar_flux_new_[se] *= normalization_factor;
   }
 
 // Accumulate cell-wise ray length tallies collected this iteration, then
@@ -223,16 +213,7 @@ void FlatSourceDomain::normalize_scalar_flux_and_volumes(
 void FlatSourceDomain::set_flux_to_flux_plus_source(
   int64_t idx, double volume, int material, int g)
 {
-  // Temperature and angle indices, if using multiple temperature
-  // data sets and/or anisotropic data sets.
-  // TODO: Currently assumes we are only using single temp/single
-  // angle data.
-  const int t = 0;
-  const int a = 0;
-
-  double sigma_t = data::mg.macro_xs_[material].get_xs(
-    MgxsType::TOTAL, g, nullptr, nullptr, nullptr, t, a);
-
+  double sigma_t = sigma_t_[material * negroups_ + g];
   scalar_flux_new_[idx] /= (sigma_t * volume);
   scalar_flux_new_[idx] += source_[idx];
 }
@@ -337,13 +318,6 @@ double FlatSourceDomain::compute_k_eff(double k_eff_old) const
   double fission_rate_old = 0;
   double fission_rate_new = 0;
 
-  // Temperature and angle indices, if using multiple temperature
-  // data sets and/or anisotropic data sets.
-  // TODO: Currently assumes we are only using single temp/single
-  // angle data.
-  const int t = 0;
-  const int a = 0;
-
   // Vector for gathering fission source terms for Shannon entropy calculation
   vector<float> p(n_source_regions_, 0.0f);
 
@@ -363,8 +337,7 @@ double FlatSourceDomain::compute_k_eff(double k_eff_old) const
 
     for (int g = 0; g < negroups_; g++) {
       int64_t idx = (sr * negroups_) + g;
-      double nu_sigma_f = data::mg.macro_xs_[material].get_xs(
-        MgxsType::NU_FISSION, g, nullptr, nullptr, nullptr, t, a);
+      double nu_sigma_f = nu_sigma_f_[material * negroups_ + g];
       sr_fission_source_old += nu_sigma_f * scalar_flux_old_[idx];
       sr_fission_source_new += nu_sigma_f * scalar_flux_new_[idx];
     }
@@ -548,7 +521,7 @@ double FlatSourceDomain::compute_fixed_source_normalization_factor() const
 {
   // If we are not in fixed source mode, then there are no external sources
   // so no normalization is needed.
-  if (settings::run_mode != RunMode::FIXED_SOURCE) {
+  if (settings::run_mode != RunMode::FIXED_SOURCE || adjoint_) {
     return 1.0;
   }
 
@@ -559,17 +532,10 @@ double FlatSourceDomain::compute_fixed_source_normalization_factor() const
   for (int sr = 0; sr < n_source_regions_; sr++) {
     int material = material_[sr];
     double volume = volume_[sr] * simulation_volume_;
-    for (int e = 0; e < negroups_; e++) {
-      // Temperature and angle indices, if using multiple temperature
-      // data sets and/or anisotropic data sets.
-      // TODO: Currently assumes we are only using single temp/single
-      // angle data.
-      const int t = 0;
-      const int a = 0;
-      double sigma_t = data::mg.macro_xs_[material].get_xs(
-        MgxsType::TOTAL, e, nullptr, nullptr, nullptr, t, a);
+    for (int g = 0; g < negroups_; g++) {
+      double sigma_t = sigma_t_[material * negroups_ + g];
       simulation_external_source_strength +=
-        external_source_[sr * negroups_ + e] * sigma_t * volume;
+        external_source_[sr * negroups_ + g] * sigma_t * volume;
     }
   }
 
@@ -603,13 +569,6 @@ void FlatSourceDomain::random_ray_tally()
   // Reset our tally volumes to zero
   reset_tally_volumes();
 
-  // Temperature and angle indices, if using multiple temperature
-  // data sets and/or anisotropic data sets.
-  // TODO: Currently assumes we are only using single temp/single
-  // angle data.
-  const int t = 0;
-  const int a = 0;
-
   double source_normalization_factor =
     compute_fixed_source_normalization_factor();
 
@@ -644,21 +603,15 @@ void FlatSourceDomain::random_ray_tally()
           break;
 
         case SCORE_TOTAL:
-          score = flux * volume *
-                  data::mg.macro_xs_[material].get_xs(
-                    MgxsType::TOTAL, g, NULL, NULL, NULL, t, a);
+          score = flux * volume * sigma_t_[material * negroups_ + g];
           break;
 
         case SCORE_FISSION:
-          score = flux * volume *
-                  data::mg.macro_xs_[material].get_xs(
-                    MgxsType::FISSION, g, NULL, NULL, NULL, t, a);
+          score = flux * volume * sigma_f_[material * negroups_ + g];
           break;
 
         case SCORE_NU_FISSION:
-          score = flux * volume *
-                  data::mg.macro_xs_[material].get_xs(
-                    MgxsType::NU_FISSION, g, NULL, NULL, NULL, t, a);
+          score = flux * volume * nu_sigma_f_[material * negroups_ + g];
           break;
 
         case SCORE_EVENTS:
@@ -957,9 +910,8 @@ void FlatSourceDomain::output_to_vtk() const
       for (int g = 0; g < negroups_; g++) {
         int64_t source_element = fsr * negroups_ + g;
         float flux = evaluate_flux_at_point(voxel_positions[i], fsr, g);
-        float Sigma_f = data::mg.macro_xs_[mat].get_xs(
-          MgxsType::FISSION, g, nullptr, nullptr, nullptr, 0, 0);
-        total_fission += Sigma_f * flux;
+        double sigma_f = sigma_f_[mat * negroups_ + g];
+        total_fission += sigma_f * flux;
       }
       total_fission = convert_to_big_endian<float>(total_fission);
       std::fwrite(&total_fission, sizeof(float), 1, plot);
@@ -977,10 +929,10 @@ void FlatSourceDomain::apply_external_source_to_source_region(
   const auto& discrete_energies = discrete->x();
   const auto& discrete_probs = discrete->prob();
 
-  for (int e = 0; e < discrete_energies.size(); e++) {
-    int g = data::mg.get_group_index(discrete_energies[e]);
+  for (int i = 0; i < discrete_energies.size(); i++) {
+    int g = data::mg.get_group_index(discrete_energies[i]);
     external_source_[source_region * negroups_ + g] +=
-      discrete_probs[e] * strength_factor;
+      discrete_probs[i] * strength_factor;
   }
 }
 
@@ -1074,27 +1026,107 @@ void FlatSourceDomain::convert_external_sources()
     }
   } // End loop over external sources
 
+// Divide the fixed source term by sigma t (to save time when applying each
+// iteration)
+#pragma omp parallel for
+  for (int sr = 0; sr < n_source_regions_; sr++) {
+    int material = material_[sr];
+    for (int g = 0; g < negroups_; g++) {
+      double sigma_t = sigma_t_[material * negroups_ + g];
+      external_source_[sr * negroups_ + g] /= sigma_t;
+    }
+  }
+}
+
+void FlatSourceDomain::flux_swap()
+{
+  scalar_flux_old_.swap(scalar_flux_new_);
+}
+
+void FlatSourceDomain::flatten_xs()
+{
   // Temperature and angle indices, if using multiple temperature
   // data sets and/or anisotropic data sets.
   // TODO: Currently assumes we are only using single temp/single angle data.
   const int t = 0;
   const int a = 0;
 
-// Divide the fixed source term by sigma t (to save time when applying each
-// iteration)
+  n_materials_ = data::mg.macro_xs_.size();
+  for (auto& m : data::mg.macro_xs_) {
+    for (int g_out = 0; g_out < negroups_; g_out++) {
+      if (m.exists_in_model) {
+        double sigma_t =
+          m.get_xs(MgxsType::TOTAL, g_out, NULL, NULL, NULL, t, a);
+        sigma_t_.push_back(sigma_t);
+
+        double nu_Sigma_f =
+          m.get_xs(MgxsType::NU_FISSION, g_out, NULL, NULL, NULL, t, a);
+        nu_sigma_f_.push_back(nu_Sigma_f);
+
+        double sigma_f =
+          m.get_xs(MgxsType::FISSION, g_out, NULL, NULL, NULL, t, a);
+        sigma_f_.push_back(sigma_f);
+
+        double chi =
+          m.get_xs(MgxsType::CHI_PROMPT, g_out, &g_out, NULL, NULL, t, a);
+        chi_.push_back(chi);
+
+        for (int g_in = 0; g_in < negroups_; g_in++) {
+          double sigma_s =
+            m.get_xs(MgxsType::NU_SCATTER, g_in, &g_out, NULL, NULL, t, a);
+          sigma_s_.push_back(sigma_s);
+        }
+      } else {
+        sigma_t_.push_back(0);
+        nu_sigma_f_.push_back(0);
+        sigma_f_.push_back(0);
+        chi_.push_back(0);
+        for (int g_in = 0; g_in < negroups_; g_in++) {
+          sigma_s_.push_back(0);
+        }
+      }
+    }
+  }
+}
+
+void FlatSourceDomain::set_adjoint_sources(const vector<double>& forward_flux)
+{
+  // Set the external source to 1/forward_flux
+  // The forward flux is given in terms of total for the forward simulation
+  // so we must convert it to a "per batch" quantity
+#pragma omp parallel for
+  for (int64_t se = 0; se < n_source_elements_; se++) {
+    external_source_[se] = 1.0 / forward_flux[se];
+  }
+
+  // Divide the fixed source term by sigma t (to save time when applying each
+  // iteration)
 #pragma omp parallel for
   for (int sr = 0; sr < n_source_regions_; sr++) {
     int material = material_[sr];
-    for (int e = 0; e < negroups_; e++) {
-      double sigma_t = data::mg.macro_xs_[material].get_xs(
-        MgxsType::TOTAL, e, nullptr, nullptr, nullptr, t, a);
-      external_source_[sr * negroups_ + e] /= sigma_t;
+    for (int g = 0; g < negroups_; g++) {
+      double sigma_t = sigma_t_[material * negroups_ + g];
+      external_source_[sr * negroups_ + g] /= sigma_t;
     }
   }
 }
-void FlatSourceDomain::flux_swap()
+
+void FlatSourceDomain::transpose_scattering_matrix()
 {
-  scalar_flux_old_.swap(scalar_flux_new_);
+  // Transpose the inner two dimensions for each material
+  for (int m = 0; m < n_materials_; ++m) {
+    int material_offset = m * negroups_ * negroups_;
+    for (int i = 0; i < negroups_; ++i) {
+      for (int j = i + 1; j < negroups_; ++j) {
+        // Calculate indices of the elements to swap
+        int idx1 = material_offset + i * negroups_ + j;
+        int idx2 = material_offset + j * negroups_ + i;
+
+        // Swap the elements to transpose the matrix
+        std::swap(sigma_s_[idx1], sigma_s_[idx2]);
+      }
+    }
+  }
 }
 
 } // namespace openmc
diff --git a/src/random_ray/linear_source_domain.cpp b/src/random_ray/linear_source_domain.cpp
index 5c3fa91c182..f3f6fa0687d 100644
--- a/src/random_ray/linear_source_domain.cpp
+++ b/src/random_ray/linear_source_domain.cpp
@@ -56,40 +56,30 @@ void LinearSourceDomain::update_neutron_source(double k_eff)
 
   double inverse_k_eff = 1.0 / k_eff;
 
-  // Temperature and angle indices, if using multiple temperature
-  // data sets and/or anisotropic data sets.
-  // TODO: Currently assumes we are only using single temp/single
-  // angle data.
-  const int t = 0;
-  const int a = 0;
-
 #pragma omp parallel for
   for (int sr = 0; sr < n_source_regions_; sr++) {
 
     int material = material_[sr];
     MomentMatrix invM = mom_matrix_[sr].inverse();
 
-    for (int e_out = 0; e_out < negroups_; e_out++) {
-      double sigma_t = data::mg.macro_xs_[material].get_xs(
-        MgxsType::TOTAL, e_out, nullptr, nullptr, nullptr, t, a);
+    for (int g_out = 0; g_out < negroups_; g_out++) {
+      double sigma_t = sigma_t_[material * negroups_ + g_out];
 
       double scatter_flat = 0.0f;
       double fission_flat = 0.0f;
       MomentArray scatter_linear = {0.0, 0.0, 0.0};
       MomentArray fission_linear = {0.0, 0.0, 0.0};
 
-      for (int e_in = 0; e_in < negroups_; e_in++) {
+      for (int g_in = 0; g_in < negroups_; g_in++) {
         // Handles for the flat and linear components of the flux
-        double flux_flat = scalar_flux_old_[sr * negroups_ + e_in];
-        MomentArray flux_linear = flux_moments_old_[sr * negroups_ + e_in];
+        double flux_flat = scalar_flux_old_[sr * negroups_ + g_in];
+        MomentArray flux_linear = flux_moments_old_[sr * negroups_ + g_in];
 
         // Handles for cross sections
-        double sigma_s = data::mg.macro_xs_[material].get_xs(
-          MgxsType::NU_SCATTER, e_in, &e_out, nullptr, nullptr, t, a);
-        double nu_sigma_f = data::mg.macro_xs_[material].get_xs(
-          MgxsType::NU_FISSION, e_in, nullptr, nullptr, nullptr, t, a);
-        double chi = data::mg.macro_xs_[material].get_xs(
-          MgxsType::CHI_PROMPT, e_in, &e_out, nullptr, nullptr, t, a);
+        double sigma_s =
+          sigma_s_[material * negroups_ * negroups_ + g_out * negroups_ + g_in];
+        double nu_sigma_f = nu_sigma_f_[material * negroups_ + g_in];
+        double chi = chi_[material * negroups_ + g_out];
 
         // Compute source terms for flat and linear components of the flux
         scatter_flat += sigma_s * flux_flat;
@@ -99,7 +89,7 @@ void LinearSourceDomain::update_neutron_source(double k_eff)
       }
 
       // Compute the flat source term
-      source_[sr * negroups_ + e_out] =
+      source_[sr * negroups_ + g_out] =
         (scatter_flat + fission_flat * inverse_k_eff) / sigma_t;
 
       // Compute the linear source terms
@@ -107,7 +97,7 @@ void LinearSourceDomain::update_neutron_source(double k_eff)
       // are not well known, we will leave the source gradients as zero
       // so as to avoid causing any numerical instability.
       if (simulation::current_batch > 10) {
-        source_gradients_[sr * negroups_ + e_out] =
+        source_gradients_[sr * negroups_ + g_out] =
           invM * ((scatter_linear + fission_linear * inverse_k_eff) / sigma_t);
       }
     }
@@ -116,7 +106,7 @@ void LinearSourceDomain::update_neutron_source(double k_eff)
   if (settings::run_mode == RunMode::FIXED_SOURCE) {
 // Add external source to flat source term if in fixed source mode
 #pragma omp parallel for
-    for (int se = 0; se < n_source_elements_; se++) {
+    for (int64_t se = 0; se < n_source_elements_; se++) {
       source_[se] += external_source_[se];
     }
   }
@@ -133,9 +123,9 @@ void LinearSourceDomain::normalize_scalar_flux_and_volumes(
 
 // Normalize flux to total distance travelled by all rays this iteration
 #pragma omp parallel for
-  for (int64_t e = 0; e < scalar_flux_new_.size(); e++) {
-    scalar_flux_new_[e] *= normalization_factor;
-    flux_moments_new_[e] *= normalization_factor;
+  for (int64_t se = 0; se < scalar_flux_new_.size(); se++) {
+    scalar_flux_new_[se] *= normalization_factor;
+    flux_moments_new_[se] *= normalization_factor;
   }
 
 // Accumulate cell-wise ray length tallies collected this iteration, then
diff --git a/src/random_ray/random_ray.cpp b/src/random_ray/random_ray.cpp
index eb64cb7d26e..7a359f35664 100644
--- a/src/random_ray/random_ray.cpp
+++ b/src/random_ray/random_ray.cpp
@@ -316,17 +316,9 @@ void RandomRay::attenuate_flux_flat_source(double distance, bool is_active)
   int64_t source_element = source_region * negroups_;
   int material = this->material();
 
-  // Temperature and angle indices, if using multiple temperature
-  // data sets and/or anisotropic data sets.
-  // TODO: Currently assumes we are only using single temp/single
-  // angle data.
-  const int t = 0;
-  const int a = 0;
-
   // MOC incoming flux attenuation + source contribution/attenuation equation
   for (int g = 0; g < negroups_; g++) {
-    float sigma_t = data::mg.macro_xs_[material].get_xs(
-      MgxsType::TOTAL, g, NULL, NULL, NULL, t, a);
+    float sigma_t = domain_->sigma_t_[material * negroups_ + g];
     float tau = sigma_t * distance;
     float exponential = cjosey_exponential(tau); // exponential = 1 - exp(-tau)
     float new_delta_psi =
@@ -388,13 +380,6 @@ void RandomRay::attenuate_flux_linear_source(double distance, bool is_active)
   int64_t source_element = source_region * negroups_;
   int material = this->material();
 
-  // Temperature and angle indices, if using multiple temperature
-  // data sets and/or anisotropic data sets.
-  // TODO: Currently assumes we are only using single temp/single
-  // angle data.
-  const int t = 0;
-  const int a = 0;
-
   Position& centroid = domain->centroid_[source_region];
   Position midpoint = r() + u() * (distance / 2.0);
 
@@ -422,8 +407,7 @@ void RandomRay::attenuate_flux_linear_source(double distance, bool is_active)
   for (int g = 0; g < negroups_; g++) {
 
     // Compute tau, the optical thickness of the ray segment
-    float sigma_t = data::mg.macro_xs_[material].get_xs(
-      MgxsType::TOTAL, g, NULL, NULL, NULL, t, a);
+    float sigma_t = domain_->sigma_t_[material * negroups_ + g];
     float tau = sigma_t * distance;
 
     // If tau is very small, set it to zero to avoid numerical issues.
diff --git a/src/random_ray/random_ray_simulation.cpp b/src/random_ray/random_ray_simulation.cpp
index 57959e80179..a9180c68e7b 100644
--- a/src/random_ray/random_ray_simulation.cpp
+++ b/src/random_ray/random_ray_simulation.cpp
@@ -23,6 +23,22 @@ namespace openmc {
 
 void openmc_run_random_ray()
 {
+  //////////////////////////////////////////////////////////
+  // Run forward simulation
+  //////////////////////////////////////////////////////////
+
+  // Check if adjoint calculation is needed. If it is, we will run the forward
+  // calculation first and then the adjoint calculation later.
+  bool adjoint_needed = FlatSourceDomain::adjoint_;
+
+  // Configure the domain for forward simulation
+  FlatSourceDomain::adjoint_ = false;
+
+  // If we're going to do an adjoint simulation afterwards, report that this is
+  // the initial forward flux solve.
+  if (adjoint_needed && mpi::master)
+    header("FORWARD FLUX SOLVE", 3);
+
   // Initialize OpenMC general data structures
   openmc_simulation_init();
 
@@ -30,26 +46,93 @@ void openmc_run_random_ray()
   if (mpi::master)
     validate_random_ray_inputs();
 
-  // Initialize Random Ray Simulation Object
-  RandomRaySimulation sim;
+  // Declare forward flux so that it can be saved for later adjoint simulation
+  vector<double> forward_flux;
+
+  {
+    // Initialize Random Ray Simulation Object
+    RandomRaySimulation sim;
+
+    // Initialize fixed sources, if present
+    sim.prepare_fixed_sources();
+
+    // Begin main simulation timer
+    simulation::time_total.start();
+
+    // Execute random ray simulation
+    sim.simulate();
+
+    // End main simulation timer
+    simulation::time_total.stop();
+
+    // Normalize and save the final forward flux
+    forward_flux = sim.domain()->scalar_flux_final_;
+
+    double source_normalization_factor =
+      sim.domain()->compute_fixed_source_normalization_factor() /
+      (settings::n_batches - settings::n_inactive);
+
+#pragma omp parallel for
+    for (uint64_t i = 0; i < forward_flux.size(); i++) {
+      forward_flux[i] *= source_normalization_factor;
+    }
+
+    // Finalize OpenMC
+    openmc_simulation_finalize();
+
+    // Reduce variables across MPI ranks
+    sim.reduce_simulation_statistics();
+
+    // Output all simulation results
+    sim.output_simulation_results();
+  }
+
+  //////////////////////////////////////////////////////////
+  // Run adjoint simulation (if enabled)
+  //////////////////////////////////////////////////////////
+
+  if (adjoint_needed) {
+    reset_timers();
+
+    // Configure the domain for adjoint simulation
+    FlatSourceDomain::adjoint_ = true;
+
+    if (mpi::master)
+      header("ADJOINT FLUX SOLVE", 3);
+
+    // Initialize OpenMC general data structures
+    openmc_simulation_init();
 
-  // Begin main simulation timer
-  simulation::time_total.start();
+    // Initialize Random Ray Simulation Object
+    RandomRaySimulation adjoint_sim;
 
-  // Execute random ray simulation
-  sim.simulate();
+    // Initialize adjoint fixed sources, if present
+    adjoint_sim.prepare_fixed_sources_adjoint(forward_flux);
 
-  // End main simulation timer
-  openmc::simulation::time_total.stop();
+    // Transpose scattering matrix
+    adjoint_sim.domain()->transpose_scattering_matrix();
 
-  // Finalize OpenMC
-  openmc_simulation_finalize();
+    // Swap nu_sigma_f and chi
+    adjoint_sim.domain()->nu_sigma_f_.swap(adjoint_sim.domain()->chi_);
 
-  // Reduce variables across MPI ranks
-  sim.reduce_simulation_statistics();
+    // Begin main simulation timer
+    simulation::time_total.start();
 
-  // Output all simulation results
-  sim.output_simulation_results();
+    // Execute random ray simulation
+    adjoint_sim.simulate();
+
+    // End main simulation timer
+    simulation::time_total.stop();
+
+    // Finalize OpenMC
+    openmc_simulation_finalize();
+
+    // Reduce variables across MPI ranks
+    adjoint_sim.reduce_simulation_statistics();
+
+    // Output all simulation results
+    adjoint_sim.output_simulation_results();
+  }
 }
 
 // Enforces restrictions on inputs in random ray mode.  While there are
@@ -254,16 +337,31 @@ RandomRaySimulation::RandomRaySimulation()
   default:
     fatal_error("Unknown random ray source shape");
   }
+
+  // Convert OpenMC native MGXS into a more efficient format
+  // internal to the random ray solver
+  domain_->flatten_xs();
 }
 
-void RandomRaySimulation::simulate()
+void RandomRaySimulation::prepare_fixed_sources()
 {
   if (settings::run_mode == RunMode::FIXED_SOURCE) {
     // Transfer external source user inputs onto random ray source regions
     domain_->convert_external_sources();
     domain_->count_external_source_regions();
   }
+}
+
+void RandomRaySimulation::prepare_fixed_sources_adjoint(
+  vector<double>& forward_flux)
+{
+  if (settings::run_mode == RunMode::FIXED_SOURCE) {
+    domain_->set_adjoint_sources(forward_flux);
+  }
+}
 
+void RandomRaySimulation::simulate()
+{
   // Random ray power iteration loop
   while (simulation::current_batch < settings::n_batches) {
 
@@ -314,18 +412,20 @@ void RandomRaySimulation::simulate()
     }
 
     // Execute all tallying tasks, if this is an active batch
-    if (simulation::current_batch > settings::n_inactive && mpi::master) {
-
-      // Generate mapping between source regions and tallies
-      if (!domain_->mapped_all_tallies_) {
-        domain_->convert_source_regions_to_tallies();
-      }
-
-      // Use above mapping to contribute FSR flux data to appropriate tallies
-      domain_->random_ray_tally();
+    if (simulation::current_batch > settings::n_inactive) {
 
       // Add this iteration's scalar flux estimate to final accumulated estimate
       domain_->accumulate_iteration_flux();
+
+      if (mpi::master) {
+        // Generate mapping between source regions and tallies
+        if (!domain_->mapped_all_tallies_) {
+          domain_->convert_source_regions_to_tallies();
+        }
+
+        // Use above mapping to contribute FSR flux data to appropriate tallies
+        domain_->random_ray_tally();
+      }
     }
 
     // Set phi_old = phi_new
@@ -448,6 +548,9 @@ void RandomRaySimulation::print_results_random_ray(
     }
     fmt::print(" Volume Estimator Type             = {}\n", estimator);
 
+    std::string adjoint_true = (FlatSourceDomain::adjoint_) ? "ON" : "OFF";
+    fmt::print(" Adjoint Flux Mode                 = {}\n", adjoint_true);
+
     header("Timing Statistics", 4);
     show_time("Total time for initialization", time_initialize.elapsed());
     show_time("Reading cross sections", time_read_xs.elapsed(), 1);
diff --git a/src/settings.cpp b/src/settings.cpp
index 3a905bdf961..d52177ae88f 100644
--- a/src/settings.cpp
+++ b/src/settings.cpp
@@ -301,6 +301,10 @@ void get_run_parameters(pugi::xml_node node_base)
       FlatSourceDomain::volume_normalized_flux_tallies_ =
         get_node_value_bool(random_ray_node, "volume_normalized_flux_tallies");
     }
+    if (check_for_node(random_ray_node, "adjoint")) {
+      FlatSourceDomain::adjoint_ =
+        get_node_value_bool(random_ray_node, "adjoint");
+    }
   }
 }
 
diff --git a/tests/regression_tests/random_ray_adjoint_fixed_source/__init__.py b/tests/regression_tests/random_ray_adjoint_fixed_source/__init__.py
new file mode 100644
index 00000000000..e69de29bb2d
diff --git a/tests/regression_tests/random_ray_adjoint_fixed_source/inputs_true.dat b/tests/regression_tests/random_ray_adjoint_fixed_source/inputs_true.dat
new file mode 100644
index 00000000000..686987512a0
--- /dev/null
+++ b/tests/regression_tests/random_ray_adjoint_fixed_source/inputs_true.dat
@@ -0,0 +1,245 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <cross_sections>mgxs.h5</cross_sections>
+    <material id="1" name="source">
+      <density units="macro" value="1.0"/>
+      <macroscopic name="source"/>
+    </material>
+    <material id="2" name="void">
+      <density units="macro" value="1.0"/>
+      <macroscopic name="void"/>
+    </material>
+    <material id="3" name="absorber">
+      <density units="macro" value="1.0"/>
+      <macroscopic name="absorber"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" name="infinite source region" universe="1"/>
+    <cell id="2" material="2" name="infinite void region" universe="2"/>
+    <cell id="3" material="3" name="infinite absorber region" universe="3"/>
+    <cell fill="4" id="4" universe="5"/>
+    <cell fill="5" id="5" name="full domain" region="1 -2 3 -4 5 -6" universe="6"/>
+    <lattice id="4">
+      <pitch>2.5 2.5 2.5</pitch>
+      <dimension>12 12 12</dimension>
+      <lower_left>0.0 0.0 0.0</lower_left>
+      <universes>
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+1 1 2 2 2 2 2 2 2 2 3 3 
+1 1 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+1 1 2 2 2 2 2 2 2 2 3 3 
+1 1 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 </universes>
+    </lattice>
+    <surface boundary="reflective" coeffs="0.0" id="1" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="30.0" id="2" type="x-plane"/>
+    <surface boundary="reflective" coeffs="0.0" id="3" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="30.0" id="4" type="y-plane"/>
+    <surface boundary="reflective" coeffs="0.0" id="5" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="30.0" id="6" type="z-plane"/>
+  </geometry>
+  <settings>
+    <run_mode>fixed source</run_mode>
+    <particles>90</particles>
+    <batches>10</batches>
+    <inactive>5</inactive>
+    <source particle="neutron" strength="3.14" type="independent">
+      <energy type="discrete">
+        <parameters>100.0 1.0</parameters>
+      </energy>
+      <constraints>
+        <domain_type>universe</domain_type>
+        <domain_ids>1</domain_ids>
+      </constraints>
+    </source>
+    <energy_mode>multi-group</energy_mode>
+    <random_ray>
+      <distance_active>500.0</distance_active>
+      <distance_inactive>100.0</distance_inactive>
+      <source particle="neutron" strength="1.0" type="independent">
+        <space type="box">
+          <parameters>0.0 0.0 0.0 30.0 30.0 30.0</parameters>
+        </space>
+      </source>
+      <volume_normalized_flux_tallies>True</volume_normalized_flux_tallies>
+      <adjoint>True</adjoint>
+    </random_ray>
+  </settings>
+  <tallies>
+    <filter id="3" type="material">
+      <bins>1</bins>
+    </filter>
+    <filter id="2" type="material">
+      <bins>2</bins>
+    </filter>
+    <filter id="1" type="material">
+      <bins>3</bins>
+    </filter>
+    <tally id="3" name="Source Tally">
+      <filters>3</filters>
+      <scores>flux</scores>
+      <estimator>tracklength</estimator>
+    </tally>
+    <tally id="2" name="Void Tally">
+      <filters>2</filters>
+      <scores>flux</scores>
+      <estimator>tracklength</estimator>
+    </tally>
+    <tally id="1" name="Absorber Tally">
+      <filters>1</filters>
+      <scores>flux</scores>
+      <estimator>tracklength</estimator>
+    </tally>
+  </tallies>
+</model>
diff --git a/tests/regression_tests/random_ray_adjoint_fixed_source/results_true.dat b/tests/regression_tests/random_ray_adjoint_fixed_source/results_true.dat
new file mode 100644
index 00000000000..a216fa9fcc3
--- /dev/null
+++ b/tests/regression_tests/random_ray_adjoint_fixed_source/results_true.dat
@@ -0,0 +1,9 @@
+tally 1:
+-7.235364E+03
+3.367109E+09
+tally 2:
+4.818311E+05
+6.269371E+10
+tally 3:
+1.515641E+06
+4.598791E+11
diff --git a/tests/regression_tests/random_ray_adjoint_fixed_source/test.py b/tests/regression_tests/random_ray_adjoint_fixed_source/test.py
new file mode 100644
index 00000000000..0295e36e9ac
--- /dev/null
+++ b/tests/regression_tests/random_ray_adjoint_fixed_source/test.py
@@ -0,0 +1,20 @@
+import os
+
+from openmc.examples import random_ray_three_region_cube
+
+from tests.testing_harness import TolerantPyAPITestHarness
+
+
+class MGXSTestHarness(TolerantPyAPITestHarness):
+    def _cleanup(self):
+        super()._cleanup()
+        f = 'mgxs.h5'
+        if os.path.exists(f):
+            os.remove(f)
+
+
+def test_random_ray_adjoint_fixed_source():
+    model = random_ray_three_region_cube()
+    model.settings.random_ray['adjoint'] = True
+    harness = MGXSTestHarness('statepoint.10.h5', model)
+    harness.main()
diff --git a/tests/regression_tests/random_ray_adjoint_k_eff/__init__.py b/tests/regression_tests/random_ray_adjoint_k_eff/__init__.py
new file mode 100644
index 00000000000..e69de29bb2d
diff --git a/tests/regression_tests/random_ray_adjoint_k_eff/inputs_true.dat b/tests/regression_tests/random_ray_adjoint_k_eff/inputs_true.dat
new file mode 100644
index 00000000000..725702a4912
--- /dev/null
+++ b/tests/regression_tests/random_ray_adjoint_k_eff/inputs_true.dat
@@ -0,0 +1,110 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <cross_sections>mgxs.h5</cross_sections>
+    <material id="1" name="UO2 fuel">
+      <density units="macro" value="1.0"/>
+      <macroscopic name="UO2"/>
+    </material>
+    <material id="2" name="Water">
+      <density units="macro" value="1.0"/>
+      <macroscopic name="LWTR"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" name="fuel inner a" region="-2" universe="1"/>
+    <cell id="2" material="1" name="fuel inner b" region="2 -3" universe="1"/>
+    <cell id="3" material="1" name="fuel inner c" region="3 -1" universe="1"/>
+    <cell id="4" material="2" name="moderator inner a" region="1 -4" universe="1"/>
+    <cell id="5" material="2" name="moderator outer b" region="4 -5" universe="1"/>
+    <cell id="6" material="2" name="moderator outer c" region="5" universe="1"/>
+    <cell fill="1" id="7" name="azimuthal_cell_0" region="6 -7" universe="2"/>
+    <cell fill="1" id="8" name="azimuthal_cell_1" region="7 -8" universe="2"/>
+    <cell fill="1" id="9" name="azimuthal_cell_2" region="8 -9" universe="2"/>
+    <cell fill="1" id="10" name="azimuthal_cell_3" region="9 -10" universe="2"/>
+    <cell fill="1" id="11" name="azimuthal_cell_4" region="10 -11" universe="2"/>
+    <cell fill="1" id="12" name="azimuthal_cell_5" region="11 -12" universe="2"/>
+    <cell fill="1" id="13" name="azimuthal_cell_6" region="12 -13" universe="2"/>
+    <cell fill="1" id="14" name="azimuthal_cell_7" region="13 -6" universe="2"/>
+    <cell id="15" material="2" name="moderator infinite" universe="3"/>
+    <cell fill="4" id="16" universe="5"/>
+    <cell fill="6" id="17" name="assembly" region="14 -15 16 -17" universe="7"/>
+    <lattice id="4">
+      <pitch>0.126 0.126</pitch>
+      <dimension>10 10</dimension>
+      <lower_left>-0.63 -0.63</lower_left>
+      <universes>
+3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 </universes>
+    </lattice>
+    <lattice id="6">
+      <pitch>1.26 1.26</pitch>
+      <dimension>2 2</dimension>
+      <lower_left>-1.26 -1.26</lower_left>
+      <universes>
+2 2 
+2 5 </universes>
+    </lattice>
+    <surface coeffs="0.0 0.0 0.54" id="1" name="Fuel OR" type="z-cylinder"/>
+    <surface coeffs="0.0 0.0 0.33" id="2" name="inner ring a" type="z-cylinder"/>
+    <surface coeffs="0.0 0.0 0.45" id="3" name="inner ring b" type="z-cylinder"/>
+    <surface coeffs="0.0 0.0 0.6" id="4" name="outer ring a" type="z-cylinder"/>
+    <surface coeffs="0.0 0.0 0.69" id="5" name="outer ring b" type="z-cylinder"/>
+    <surface coeffs="-0.0 1.0 0 0" id="6" type="plane"/>
+    <surface coeffs="-0.7071067811865475 0.7071067811865476 0 0" id="7" type="plane"/>
+    <surface coeffs="-1.0 6.123233995736766e-17 0 0" id="8" type="plane"/>
+    <surface coeffs="-0.7071067811865476 -0.7071067811865475 0 0" id="9" type="plane"/>
+    <surface coeffs="-1.2246467991473532e-16 -1.0 0 0" id="10" type="plane"/>
+    <surface coeffs="0.7071067811865475 -0.7071067811865477 0 0" id="11" type="plane"/>
+    <surface coeffs="1.0 -1.8369701987210297e-16 0 0" id="12" type="plane"/>
+    <surface coeffs="0.7071067811865477 0.7071067811865474 0 0" id="13" type="plane"/>
+    <surface boundary="reflective" coeffs="-1.26" id="14" name="minimum x" type="x-plane"/>
+    <surface boundary="reflective" coeffs="1.26" id="15" name="maximum x" type="x-plane"/>
+    <surface boundary="reflective" coeffs="-1.26" id="16" name="minimum y" type="y-plane"/>
+    <surface boundary="reflective" coeffs="1.26" id="17" name="maximum y" type="y-plane"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>10</batches>
+    <inactive>5</inactive>
+    <energy_mode>multi-group</energy_mode>
+    <random_ray>
+      <distance_active>100.0</distance_active>
+      <distance_inactive>20.0</distance_inactive>
+      <source particle="neutron" strength="1.0" type="independent">
+        <space type="box">
+          <parameters>-1.26 -1.26 -1 1.26 1.26 1</parameters>
+        </space>
+      </source>
+      <volume_normalized_flux_tallies>True</volume_normalized_flux_tallies>
+      <adjoint>True</adjoint>
+    </random_ray>
+  </settings>
+  <tallies>
+    <mesh id="1">
+      <dimension>2 2</dimension>
+      <lower_left>-1.26 -1.26</lower_left>
+      <upper_right>1.26 1.26</upper_right>
+    </mesh>
+    <filter id="1" type="mesh">
+      <bins>1</bins>
+    </filter>
+    <filter id="2" type="energy">
+      <bins>1e-05 0.0635 10.0 100.0 1000.0 500000.0 1000000.0 20000000.0</bins>
+    </filter>
+    <tally id="1" name="Mesh tally">
+      <filters>1 2</filters>
+      <scores>flux fission nu-fission</scores>
+      <estimator>analog</estimator>
+    </tally>
+  </tallies>
+</model>
diff --git a/tests/regression_tests/random_ray_adjoint_k_eff/results_true.dat b/tests/regression_tests/random_ray_adjoint_k_eff/results_true.dat
new file mode 100644
index 00000000000..1690d46e966
--- /dev/null
+++ b/tests/regression_tests/random_ray_adjoint_k_eff/results_true.dat
@@ -0,0 +1,171 @@
+k-combined:
+1.006640E+00 1.812967E-03
+tally 1:
+6.684129E+00
+8.939821E+00
+2.685967E+00
+1.443592E+00
+0.000000E+00
+0.000000E+00
+6.358774E+00
+8.091444E+00
+9.687217E-01
+1.878029E-01
+0.000000E+00
+0.000000E+00
+5.963160E+00
+7.117108E+00
+1.932332E-01
+7.473914E-03
+0.000000E+00
+0.000000E+00
+5.137593E+00
+5.283310E+00
+1.714616E-01
+5.884834E-03
+1.086218E-06
+2.361752E-13
+4.857253E+00
+4.719856E+00
+5.689580E-02
+6.476286E-04
+2.989356E-03
+1.787808E-06
+4.830516E+00
+4.666801E+00
+7.203015E-03
+1.037676E-05
+3.620020E+00
+2.620927E+00
+5.161382E+00
+5.328124E+00
+6.786255E-02
+9.210763E-04
+5.531943E+00
+6.120553E+00
+5.414034E+00
+5.864661E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+5.632338E+00
+6.347626E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+5.682608E+00
+6.462382E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+5.310716E+00
+5.645180E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+4.945409E+00
+4.893171E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+4.842688E+00
+4.690352E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+5.117198E+00
+5.237280E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+6.938711E+00
+9.633345E+00
+2.835258E+00
+1.608212E+00
+0.000000E+00
+0.000000E+00
+6.549505E+00
+8.584036E+00
+1.015138E+00
+2.061993E-01
+0.000000E+00
+0.000000E+00
+6.050651E+00
+7.327711E+00
+1.992816E-01
+7.948424E-03
+0.000000E+00
+0.000000E+00
+5.113981E+00
+5.234801E+00
+1.732323E-01
+6.006619E-03
+1.097435E-06
+2.410627E-13
+4.837033E+00
+4.680541E+00
+5.760042E-02
+6.637112E-04
+3.026377E-03
+1.832205E-06
+4.827049E+00
+4.660105E+00
+7.319913E-03
+1.071647E-05
+3.678770E+00
+2.706730E+00
+5.175337E+00
+5.356957E+00
+6.923046E-02
+9.586177E-04
+5.643451E+00
+6.370016E+00
+6.693323E+00
+8.964322E+00
+2.753307E+00
+1.516683E+00
+0.000000E+00
+0.000000E+00
+6.358384E+00
+8.090233E+00
+9.912008E-01
+1.965868E-01
+0.000000E+00
+0.000000E+00
+5.957484E+00
+7.103246E+00
+1.974033E-01
+7.798286E-03
+0.000000E+00
+0.000000E+00
+5.130744E+00
+5.268844E+00
+1.749233E-01
+6.123348E-03
+1.108148E-06
+2.457474E-13
+4.857340E+00
+4.720019E+00
+5.816659E-02
+6.768049E-04
+3.056125E-03
+1.868351E-06
+4.830629E+00
+4.667018E+00
+7.366289E-03
+1.085264E-05
+3.702077E+00
+2.741125E+00
+5.164864E+00
+5.335279E+00
+6.947917E-02
+9.655086E-04
+5.663725E+00
+6.415806E+00
diff --git a/tests/regression_tests/random_ray_adjoint_k_eff/test.py b/tests/regression_tests/random_ray_adjoint_k_eff/test.py
new file mode 100644
index 00000000000..44cf1182ae6
--- /dev/null
+++ b/tests/regression_tests/random_ray_adjoint_k_eff/test.py
@@ -0,0 +1,20 @@
+import os
+
+from openmc.examples import random_ray_lattice
+
+from tests.testing_harness import TolerantPyAPITestHarness
+
+
+class MGXSTestHarness(TolerantPyAPITestHarness):
+    def _cleanup(self):
+        super()._cleanup()
+        f = 'mgxs.h5'
+        if os.path.exists(f):
+            os.remove(f)
+
+
+def test_random_ray_basic():
+    model = random_ray_lattice()
+    model.settings.random_ray['adjoint'] = True
+    harness = MGXSTestHarness('statepoint.10.h5', model)
+    harness.main()