diff --git a/singularity-opac/constants/constants.hpp b/singularity-opac/constants/constants.hpp
index 580037a..b99fe95 100644
--- a/singularity-opac/constants/constants.hpp
+++ b/singularity-opac/constants/constants.hpp
@@ -1,5 +1,5 @@
 // ======================================================================
-// © 2021. Triad National Security, LLC. All rights reserved.  This
+// © 2021-2024. Triad National Security, LLC. All rights reserved.  This
 // program was produced under U.S. Government contract
 // 89233218CNA000001 for Los Alamos National Laboratory (LANL), which
 // is operated by Triad National Security, LLC for the U.S.
@@ -194,6 +194,45 @@ struct PhysicalConstants {
   static constexpr Real gA = axial_vector_coupling_constant;
 };
 
+struct RuntimePhysicalConstants {
+  template <typename T>
+  RuntimePhysicalConstants(T pc)
+      : na(pc.na), alpha(pc.alpha), h(pc.h), hbar(pc.hbar), kb(pc.kb),
+        r_gas(pc.r_gas), qe(pc.qe), c(pc.c), g_newt(pc.g_newt),
+        g_accel(pc.g_accel), me(pc.me), mp(pc.mp), mn(pc.mn), amu(pc.amu),
+        sb(pc.sb), ar(pc.ar), faraday(pc.faraday), mu0(pc.mu0), eps0(pc.eps0),
+        eV(pc.eV), Fc(pc.Fc), nu_sigma0(pc.nu_sigma0), gA(pc.gA) {}
+
+  const Real na;
+  const Real alpha;
+  const Real h;
+  const Real hbar;
+  const Real kb;
+  const Real r_gas;
+  const Real qe;
+  const Real c;
+  const Real g_newt;
+  const Real g_accel;
+  const Real me;
+  const Real mp;
+  const Real mn;
+  const Real amu;
+  const Real sb;
+  const Real ar;
+  const Real faraday;
+  const Real mu0;
+  const Real eps0;
+  const Real eV;
+  const Real Fc;
+  const Real nu_sigma0;
+  const Real gA;
+};
+
+template <typename T>
+RuntimePhysicalConstants GetRuntimePhysicalConstants(T phys_constants) {
+  return RuntimePhysicalConstants(phys_constants);
+}
+
 using PhysicalConstantsUnity =
     PhysicalConstants<BaseUnity, UnitConversionDefault>;
 using PhysicalConstantsSI =
diff --git a/singularity-opac/neutrinos/brt_neutrinos.hpp b/singularity-opac/neutrinos/brt_neutrinos.hpp
index 5721021..9138084 100644
--- a/singularity-opac/neutrinos/brt_neutrinos.hpp
+++ b/singularity-opac/neutrinos/brt_neutrinos.hpp
@@ -1,5 +1,5 @@
 // ======================================================================
-// © 2021. Triad National Security, LLC. All rights reserved.  This
+// © 2021-2024. Triad National Security, LLC. All rights reserved.  This
 // program was produced under U.S. Government contract
 // 89233218CNA000001 for Los Alamos National Laboratory (LANL), which
 // is operated by Triad National Security, LLC for the U.S.
@@ -32,6 +32,8 @@ namespace neutrinos {
 template <typename ThermalDistribution, typename pc = PhysicalConstantsCGS>
 class BRTOpacity {
  public:
+  using PC = pc;
+
   BRTOpacity() = default;
   BRTOpacity(const ThermalDistribution &dist) : dist_(dist) {}
   BRTOpacity GetOnDevice() { return *this; }
diff --git a/singularity-opac/neutrinos/gray_opacity_neutrinos.hpp b/singularity-opac/neutrinos/gray_opacity_neutrinos.hpp
index fc9a57e..830615e 100644
--- a/singularity-opac/neutrinos/gray_opacity_neutrinos.hpp
+++ b/singularity-opac/neutrinos/gray_opacity_neutrinos.hpp
@@ -1,5 +1,5 @@
 // ======================================================================
-// © 2021. Triad National Security, LLC. All rights reserved.  This
+// © 2021-2024. Triad National Security, LLC. All rights reserved.  This
 // program was produced under U.S. Government contract
 // 89233218CNA000001 for Los Alamos National Laboratory (LANL), which
 // is operated by Triad National Security, LLC for the U.S.
@@ -30,6 +30,8 @@ namespace neutrinos {
 template <typename ThermalDistribution, typename pc = PhysicalConstantsCGS>
 class GrayOpacity {
  public:
+  using PC = pc;
+
   GrayOpacity() = default;
   GrayOpacity(const Real kappa) : kappa_(kappa) {}
   GrayOpacity(const ThermalDistribution &dist, const Real kappa)
diff --git a/singularity-opac/neutrinos/mean_neutrino_variant.hpp b/singularity-opac/neutrinos/mean_neutrino_variant.hpp
index fd50081..488085b 100644
--- a/singularity-opac/neutrinos/mean_neutrino_variant.hpp
+++ b/singularity-opac/neutrinos/mean_neutrino_variant.hpp
@@ -1,5 +1,5 @@
 // ======================================================================
-// © 2021. Triad National Security, LLC. All rights reserved.  This
+// © 2021-2024. Triad National Security, LLC. All rights reserved.  This
 // program was produced under U.S. Government contract
 // 89233218CNA000001 for Los Alamos National Laboratory (LANL), which
 // is operated by Triad National Security, LLC for the U.S.
@@ -69,6 +69,16 @@ class MeanVariant {
         opac_);
   }
 
+  PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
+  GetRuntimePhysicalConstants() const {
+    return mpark::visit(
+        [=](const auto &opac) {
+          using PC = typename std::decay_t<decltype(opac)>::PC;
+          return singularity::GetRuntimePhysicalConstants(PC());
+        },
+        opac_);
+  }
+
   PORTABLE_INLINE_FUNCTION Real PlanckMeanAbsorptionCoefficient(
       const Real rho, const Real temp, const Real Ye,
       const RadiationType type) const {
@@ -91,6 +101,12 @@ class MeanVariant {
   inline void Finalize() noexcept {
     return mpark::visit([](auto &opac) { return opac.Finalize(); }, opac_);
   }
+
+#ifdef SPINER_USE_HDF
+  void Save(const std::string &filename) const {
+    return mpark::visit([=](auto &opac) { return opac.Save(filename); }, opac_);
+  }
+#endif
 };
 
 } // namespace impl
diff --git a/singularity-opac/neutrinos/mean_opacity_neutrinos.hpp b/singularity-opac/neutrinos/mean_opacity_neutrinos.hpp
index 61a6396..52e9485 100644
--- a/singularity-opac/neutrinos/mean_opacity_neutrinos.hpp
+++ b/singularity-opac/neutrinos/mean_opacity_neutrinos.hpp
@@ -1,5 +1,5 @@
 // ======================================================================
-// © 2022. Triad National Security, LLC. All rights reserved.  This
+// © 2022-2024. Triad National Security, LLC. All rights reserved.  This
 // program was produced under U.S. Government contract
 // 89233218CNA000001 for Los Alamos National Laboratory (LANL), which
 // is operated by Triad National Security, LLC for the U.S.
@@ -37,7 +37,6 @@ namespace impl {
 // TODO(BRR) Note: It is assumed that lambda is constant for all densities,
 // temperatures, and Ye
 
-template <typename pc = PhysicalConstantsCGS>
 class MeanOpacity {
 
  public:
@@ -134,6 +133,8 @@ class MeanOpacity {
                         const Real lTMax, const int NT, const Real YeMin,
                         const Real YeMax, const int NYe, Real lNuMin,
                         Real lNuMax, const int NNu, Real *lambda = nullptr) {
+    using PC = typename Opacity::PC;
+
     lkappaPlanck_.resize(NRho, NT, NYe, NEUTRINO_NTYPES);
     // index 0 is the species and is not interpolatable
     lkappaPlanck_.setRange(1, YeMin, YeMax, NYe);
@@ -159,8 +160,8 @@ class MeanOpacity {
             // Choose default temperature-specific frequency grid if frequency
             // grid not specified
             if (AUTOFREQ) {
-              lNuMin = toLog_(1.e-3 * pc::kb * fromLog_(lTMin) / pc::h);
-              lNuMax = toLog_(1.e3 * pc::kb * fromLog_(lTMax) / pc::h);
+              lNuMin = toLog_(1.e-3 * PC::kb * fromLog_(lTMin) / PC::h);
+              lNuMax = toLog_(1.e3 * PC::kb * fromLog_(lTMax) / PC::h);
             }
             const Real dlnu = (lNuMax - lNuMin) / (NNu - 1);
             // Integrate over frequency
@@ -219,10 +220,9 @@ class MeanOpacity {
 
 } // namespace impl
 
-using MeanOpacityScaleFree = impl::MeanOpacity<PhysicalConstantsUnity>;
-using MeanOpacityCGS = impl::MeanOpacity<PhysicalConstantsCGS>;
-using MeanOpacity = impl::MeanVariant<MeanOpacityScaleFree, MeanOpacityCGS,
-                                      MeanNonCGSUnits<MeanOpacityCGS>>;
+using MeanOpacityBase = impl::MeanOpacity;
+using MeanOpacity =
+    impl::MeanVariant<MeanOpacityBase, MeanNonCGSUnits<MeanOpacityBase>>;
 
 } // namespace neutrinos
 } // namespace singularity
diff --git a/singularity-opac/neutrinos/neutrino_variant.hpp b/singularity-opac/neutrinos/neutrino_variant.hpp
index ae14525..5736077 100644
--- a/singularity-opac/neutrinos/neutrino_variant.hpp
+++ b/singularity-opac/neutrinos/neutrino_variant.hpp
@@ -1,5 +1,5 @@
 // ======================================================================
-// © 2021. Triad National Security, LLC. All rights reserved.  This
+// © 2021-2024. Triad National Security, LLC. All rights reserved.  This
 // program was produced under U.S. Government contract
 // 89233218CNA000001 for Los Alamos National Laboratory (LANL), which
 // is operated by Triad National Security, LLC for the U.S.
@@ -71,6 +71,16 @@ class Variant {
         opac_);
   }
 
+  PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
+  GetRuntimePhysicalConstants() const {
+    return mpark::visit(
+        [=](const auto &opac) {
+          using PC = typename std::decay_t<decltype(opac)>::PC;
+          return singularity::GetRuntimePhysicalConstants(PC());
+        },
+        opac_);
+  }
+
   // Directional absorption coefficient with units of 1/length
   // Signature should be at least
   // rho, temp, Ye, type, nu, lambda
diff --git a/singularity-opac/neutrinos/non_cgs_neutrinos.hpp b/singularity-opac/neutrinos/non_cgs_neutrinos.hpp
index 7a760c7..400cb2b 100644
--- a/singularity-opac/neutrinos/non_cgs_neutrinos.hpp
+++ b/singularity-opac/neutrinos/non_cgs_neutrinos.hpp
@@ -1,5 +1,5 @@
 // ======================================================================
-// © 2021. Triad National Security, LLC. All rights reserved.  This
+// © 2021-2024. Triad National Security, LLC. All rights reserved.  This
 // program was produced under U.S. Government contract
 // 89233218CNA000001 for Los Alamos National Laboratory (LANL), which
 // is operated by Triad National Security, LLC for the U.S.
@@ -29,6 +29,8 @@ namespace neutrinos {
 template <typename Opac>
 class NonCGSUnits {
  public:
+  using PC = typename Opac::PC;
+
   NonCGSUnits() = default;
   NonCGSUnits(Opac &&opac, const Real time_unit, const Real mass_unit,
               const Real length_unit, const Real temp_unit)
@@ -47,6 +49,7 @@ class NonCGSUnits {
     return NonCGSUnits<Opac>(opac_.GetOnDevice(), time_unit_, mass_unit_,
                              length_unit_, temp_unit_);
   }
+
   inline void Finalize() noexcept { opac_.Finalize(); }
 
   PORTABLE_INLINE_FUNCTION
@@ -66,10 +69,11 @@ class NonCGSUnits {
   }
 
   template <typename FrequencyIndexer, typename DataIndexer>
-  PORTABLE_INLINE_FUNCTION void AbsorptionCoefficient(
-      const Real rho, const Real temp, const Real Ye, RadiationType type,
-      FrequencyIndexer &nu_bins, DataIndexer &coeffs, const int nbins,
-      Real *lambda = nullptr) const {
+  PORTABLE_INLINE_FUNCTION void
+  AbsorptionCoefficient(const Real rho, const Real temp, const Real Ye,
+                        RadiationType type, FrequencyIndexer &nu_bins,
+                        DataIndexer &coeffs, const int nbins,
+                        Real *lambda = nullptr) const {
     for (int i = 0; i < nbins; ++i) {
       nu_bins[i] *= freq_unit_;
     }
@@ -262,6 +266,12 @@ class MeanNonCGSUnits {
   PORTABLE_INLINE_FUNCTION
   int nlambda() const noexcept { return mean_opac_.nlambda(); }
 
+#ifdef SPINER_USE_HDF
+  void Save(const std::string &filename) const {
+    return mean_opac_.Save(filename);
+  }
+#endif
+
   PORTABLE_INLINE_FUNCTION
   Real PlanckMeanAbsorptionCoefficient(const Real rho, const Real temp,
                                        const Real Ye,
diff --git a/singularity-opac/neutrinos/spiner_opac_neutrinos.hpp b/singularity-opac/neutrinos/spiner_opac_neutrinos.hpp
index 07c6a7b..906bb20 100644
--- a/singularity-opac/neutrinos/spiner_opac_neutrinos.hpp
+++ b/singularity-opac/neutrinos/spiner_opac_neutrinos.hpp
@@ -1,5 +1,5 @@
 // ======================================================================
-// © 2021. Triad National Security, LLC. All rights reserved.  This
+// © 2021-2024. Triad National Security, LLC. All rights reserved.  This
 // program was produced under U.S. Government contract
 // 89233218CNA000001 for Los Alamos National Laboratory (LANL), which
 // is operated by Triad National Security, LLC for the U.S.
@@ -61,7 +61,9 @@ enum class DataStatus { Deallocated, OnDevice, OnHost };
 template <typename ThermalDistribution, typename pc = PhysicalConstantsCGS>
 class SpinerOpacity {
  public:
+  using PC = pc;
   using DataBox = Spiner::DataBox<Real>;
+
   static constexpr Real EPS = 10.0 * std::numeric_limits<Real>::min();
   static constexpr Real Hz2MeV = pc::h / (1e6 * pc::eV);
   static constexpr Real MeV2Hz = 1 / Hz2MeV;
@@ -110,7 +112,8 @@ class SpinerOpacity {
             Real J = std::max(opac.Emissivity(rho, T * MeV2K, Ye, type), 0.0);
             Real lJ = toLog_(J);
             lJ_(iRho, iT, iYe, idx) = lJ;
-            Real JYe = std::max(opac.NumberEmissivity(rho, T * MeV2K, Ye, type), 0.0);
+            Real JYe =
+                std::max(opac.NumberEmissivity(rho, T * MeV2K, Ye, type), 0.0);
             lJYe_(iRho, iT, iYe, idx) = toLog_(JYe);
             for (int ie = 0; ie < Ne; ++ie) {
               Real lE = lalphanu_.range(0).x(ie);
@@ -119,8 +122,8 @@ class SpinerOpacity {
               Real alpha = std::max(
                   opac.AbsorptionCoefficient(rho, T, Ye, type, nu), 0.0);
               lalphanu_(iRho, iT, iYe, idx, ie) = toLog_(alpha);
-              Real j = std::max(opac.EmissivityPerNuOmega(rho, T * MeV2K, Ye, type, nu),
-                                0.0);
+              Real j = std::max(
+                  opac.EmissivityPerNuOmega(rho, T * MeV2K, Ye, type, nu), 0.0);
               ljnu_(iRho, iT, iYe, idx, ie) = toLog_(j);
             }
           }
@@ -131,8 +134,8 @@ class SpinerOpacity {
 
   // DataBox constructor. Note that this constructor *shallow* copies
   // the databoxes, so they must be managed externally.
-  SpinerOpacity(const DataBox &lalphanu, const DataBox ljnu,
-                const DataBox lJ, const DataBox lJYe)
+  SpinerOpacity(const DataBox &lalphanu, const DataBox ljnu, const DataBox lJ,
+                const DataBox lJYe)
       : memoryStatus_(impl::DataStatus::OnHost), lalphanu_(lalphanu),
         ljnu_(ljnu), lJ_(lJ), lJYe_(lJYe) {}
 
diff --git a/singularity-opac/neutrinos/tophat_emissivity_neutrinos.hpp b/singularity-opac/neutrinos/tophat_emissivity_neutrinos.hpp
index ee7b9b1..bfb024b 100644
--- a/singularity-opac/neutrinos/tophat_emissivity_neutrinos.hpp
+++ b/singularity-opac/neutrinos/tophat_emissivity_neutrinos.hpp
@@ -1,5 +1,5 @@
 // ======================================================================
-// © 2021. Triad National Security, LLC. All rights reserved.  This
+// © 2021-2024. Triad National Security, LLC. All rights reserved.  This
 // program was produced under U.S. Government contract
 // 89233218CNA000001 for Los Alamos National Laboratory (LANL), which
 // is operated by Triad National Security, LLC for the U.S.
@@ -32,6 +32,8 @@ namespace neutrinos {
 template <typename ThermalDistribution, typename pc = PhysicalConstantsCGS>
 class TophatEmissivity {
  public:
+  using PC = pc;
+
   TophatEmissivity(const Real C, const Real numin, const Real numax)
       : C_(C), numin_(numin), numax_(numax) {}
   TophatEmissivity(const ThermalDistribution &dist, const Real C,
@@ -46,6 +48,7 @@ class TophatEmissivity {
     printf("Tophat emissivity. C, numin, numax = %g, %g, %g\n", C_, numin_,
            numax_);
   }
+
   inline void Finalize() noexcept {}
 
   PORTABLE_INLINE_FUNCTION
diff --git a/singularity-opac/photons/epbremsstrahlung_opacity_photons.hpp b/singularity-opac/photons/epbremsstrahlung_opacity_photons.hpp
index e4e9ccc..725681d 100644
--- a/singularity-opac/photons/epbremsstrahlung_opacity_photons.hpp
+++ b/singularity-opac/photons/epbremsstrahlung_opacity_photons.hpp
@@ -1,5 +1,5 @@
 // ======================================================================
-// © 2022. Triad National Security, LLC. All rights reserved.  This
+// © 2022-2024. Triad National Security, LLC. All rights reserved.  This
 // program was produced under U.S. Government contract
 // 89233218CNA000001 for Los Alamos National Laboratory (LANL), which
 // is operated by Triad National Security, LLC for the U.S.
@@ -32,6 +32,8 @@ namespace photons {
 template <typename pc = PhysicalConstantsCGS>
 class EPBremsstrahlungOpacity {
  public:
+  using PC = pc;
+
   EPBremsstrahlungOpacity() = default;
   EPBremsstrahlungOpacity(const PlanckDistribution<pc> &dist) : dist_(dist) {}
 
@@ -42,6 +44,10 @@ class EPBremsstrahlungOpacity {
   void PrintParams() const noexcept {
     printf("Electron-proton bremsstrahlung opacity.\n");
   }
+
+  PORTABLE_INLINE_FUNCTION
+  pc GetPhysicalConstants() const { return pc(); }
+
   inline void Finalize() noexcept {}
 
   PORTABLE_INLINE_FUNCTION
diff --git a/singularity-opac/photons/gray_opacity_photons.hpp b/singularity-opac/photons/gray_opacity_photons.hpp
index 9bc0358..dcca68c 100644
--- a/singularity-opac/photons/gray_opacity_photons.hpp
+++ b/singularity-opac/photons/gray_opacity_photons.hpp
@@ -1,5 +1,5 @@
 // ======================================================================
-// © 2021. Triad National Security, LLC. All rights reserved.  This
+// © 2021-2024. Triad National Security, LLC. All rights reserved.  This
 // program was produced under U.S. Government contract
 // 89233218CNA000001 for Los Alamos National Laboratory (LANL), which
 // is operated by Triad National Security, LLC for the U.S.
@@ -30,6 +30,8 @@ namespace photons {
 template <typename pc = PhysicalConstantsCGS>
 class GrayOpacity {
  public:
+  using PC = pc;
+
   GrayOpacity() = default;
   GrayOpacity(const Real kappa) : kappa_(kappa) {}
   GrayOpacity(const PlanckDistribution<pc> &dist, const Real kappa)
diff --git a/singularity-opac/photons/mean_opacity_photons.hpp b/singularity-opac/photons/mean_opacity_photons.hpp
index 8a58859..d4ac76b 100644
--- a/singularity-opac/photons/mean_opacity_photons.hpp
+++ b/singularity-opac/photons/mean_opacity_photons.hpp
@@ -1,5 +1,5 @@
 // ======================================================================
-// © 2022. Triad National Security, LLC. All rights reserved.  This
+// © 2022-2024. Triad National Security, LLC. All rights reserved.  This
 // program was produced under U.S. Government contract
 // 89233218CNA000001 for Los Alamos National Laboratory (LANL), which
 // is operated by Triad National Security, LLC for the U.S.
@@ -37,11 +37,9 @@ namespace impl {
 // TODO(BRR) Note: It is assumed that lambda is constant for all densities and
 // temperatures
 
-template <typename pc = PhysicalConstantsCGS>
 class MeanOpacity {
 
  public:
-  using DataBox = Spiner::DataBox<Real>;
   MeanOpacity() = default;
   template <typename Opacity>
   MeanOpacity(const Opacity &opac, const Real lRhoMin, const Real lRhoMax,
@@ -125,6 +123,8 @@ class MeanOpacity {
                         const Real lRhoMax, const int NRho, const Real lTMin,
                         const Real lTMax, const int NT, Real lNuMin,
                         Real lNuMax, const int NNu, Real *lambda = nullptr) {
+    using PC = typename Opacity::PC;
+
     lkappaPlanck_.resize(NRho, NT);
     lkappaPlanck_.setRange(0, lTMin, lTMax, NT);
     lkappaPlanck_.setRange(1, lRhoMin, lRhoMax, NRho);
@@ -142,8 +142,8 @@ class MeanOpacity {
         Real kappaRosselandNum = 0.;
         Real kappaRosselandDenom = 0.;
         if (AUTOFREQ) {
-          lNuMin = toLog_(1.e-3 * pc::kb * fromLog_(lTMin) / pc::h);
-          lNuMax = toLog_(1.e3 * pc::kb * fromLog_(lTMax) / pc::h);
+          lNuMin = toLog_(1.e-3 * PC::kb * fromLog_(lTMin) / PC::h);
+          lNuMax = toLog_(1.e3 * PC::kb * fromLog_(lTMax) / PC::h);
         }
         const Real dlnu = (lNuMax - lNuMin) / (NNu - 1);
         // Integrate over frequency
@@ -187,8 +187,8 @@ class MeanOpacity {
   PORTABLE_INLINE_FUNCTION Real fromLog_(const Real lx) const {
     return std::pow(10., lx);
   }
-  DataBox lkappaPlanck_;
-  DataBox lkappaRosseland_;
+  Spiner::DataBox<Real> lkappaPlanck_;
+  Spiner::DataBox<Real> lkappaRosseland_;
   const char *filename_;
 };
 
@@ -196,10 +196,9 @@ class MeanOpacity {
 
 } // namespace impl
 
-using MeanOpacityScaleFree = impl::MeanOpacity<PhysicalConstantsUnity>;
-using MeanOpacityCGS = impl::MeanOpacity<PhysicalConstantsCGS>;
-using MeanOpacity = impl::MeanVariant<MeanOpacityScaleFree, MeanOpacityCGS,
-                                      MeanNonCGSUnits<MeanOpacityCGS>>;
+using MeanOpacityBase = impl::MeanOpacity;
+using MeanOpacity =
+    impl::MeanVariant<MeanOpacityBase, MeanNonCGSUnits<MeanOpacityBase>>;
 
 } // namespace photons
 } // namespace singularity
diff --git a/singularity-opac/photons/mean_photon_variant.hpp b/singularity-opac/photons/mean_photon_variant.hpp
index af292ea..a259db4 100644
--- a/singularity-opac/photons/mean_photon_variant.hpp
+++ b/singularity-opac/photons/mean_photon_variant.hpp
@@ -1,5 +1,5 @@
 // ======================================================================
-// © 2022. Triad National Security, LLC. All rights reserved.  This
+// © 2022-2024. Triad National Security, LLC. All rights reserved.  This
 // program was produced under U.S. Government contract
 // 89233218CNA000001 for Los Alamos National Laboratory (LANL), which
 // is operated by Triad National Security, LLC for the U.S.
@@ -69,6 +69,16 @@ class MeanVariant {
         opac_);
   }
 
+  PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
+  GetRuntimePhysicalConstants() const {
+    return mpark::visit(
+        [=](const auto &opac) {
+          using PC = typename std::decay_t<decltype(opac)>::PC;
+          return singularity::GetRuntimePhysicalConstants(PC());
+        },
+        opac_);
+  }
+
   PORTABLE_INLINE_FUNCTION Real
   PlanckMeanAbsorptionCoefficient(const Real rho, const Real temp) const {
     return mpark::visit(
@@ -89,6 +99,12 @@ class MeanVariant {
   inline void Finalize() noexcept {
     return mpark::visit([](auto &opac) { return opac.Finalize(); }, opac_);
   }
+
+#ifdef SPINER_USE_HDF
+  void Save(const std::string &filename) const {
+    return mpark::visit([=](auto &opac) { return opac.Save(filename); }, opac_);
+  }
+#endif
 };
 
 } // namespace impl
diff --git a/singularity-opac/photons/non_cgs_photons.hpp b/singularity-opac/photons/non_cgs_photons.hpp
index a63df14..e128cf9 100644
--- a/singularity-opac/photons/non_cgs_photons.hpp
+++ b/singularity-opac/photons/non_cgs_photons.hpp
@@ -29,6 +29,8 @@ namespace photons {
 template <typename Opac>
 class NonCGSUnits {
  public:
+  using PC = typename Opac::PC;
+
   NonCGSUnits() = default;
   NonCGSUnits(Opac &&opac, const Real time_unit, const Real mass_unit,
               const Real length_unit, const Real temp_unit)
@@ -65,9 +67,10 @@ class NonCGSUnits {
   }
 
   template <typename FrequencyIndexer, typename DataIndexer>
-  PORTABLE_INLINE_FUNCTION void AbsorptionCoefficient(
-      const Real rho, const Real temp, FrequencyIndexer &nu_bins,
-      DataIndexer &coeffs, const int nbins, Real *lambda = nullptr) const {
+  PORTABLE_INLINE_FUNCTION void
+  AbsorptionCoefficient(const Real rho, const Real temp,
+                        FrequencyIndexer &nu_bins, DataIndexer &coeffs,
+                        const int nbins, Real *lambda = nullptr) const {
     for (int i = 0; i < nbins; ++i) {
       nu_bins[i] *= freq_unit_;
     }
@@ -219,6 +222,11 @@ class NonCGSUnits {
     return NoH * mass_unit_ / rho_unit_;
   }
 
+  template <typename T>
+  PORTABLE_INLINE_FUNCTION T GetPhysicalConstants() const {
+    return opac_.GetPhysicalConstants();
+  }
+
  private:
   Opac opac_;
   Real time_unit_, mass_unit_, length_unit_, temp_unit_;
@@ -246,6 +254,12 @@ class MeanNonCGSUnits {
   PORTABLE_INLINE_FUNCTION
   int nlambda() const noexcept { return mean_opac_.nlambda(); }
 
+#ifdef SPINER_USE_HDF
+  void Save(const std::string &filename) const {
+    return mean_opac_.Save(filename);
+  }
+#endif
+
   PORTABLE_INLINE_FUNCTION
   Real PlanckMeanAbsorptionCoefficient(const Real rho, const Real temp) const {
     const Real alpha = mean_opac_.PlanckMeanAbsorptionCoefficient(
diff --git a/singularity-opac/photons/photon_variant.hpp b/singularity-opac/photons/photon_variant.hpp
index 440e5e8..3386b8d 100644
--- a/singularity-opac/photons/photon_variant.hpp
+++ b/singularity-opac/photons/photon_variant.hpp
@@ -1,5 +1,5 @@
 // ======================================================================
-// © 2021. Triad National Security, LLC. All rights reserved.  This
+// © 2021-2024. Triad National Security, LLC. All rights reserved.  This
 // program was produced under U.S. Government contract
 // 89233218CNA000001 for Los Alamos National Laboratory (LANL), which
 // is operated by Triad National Security, LLC for the U.S.
@@ -71,6 +71,16 @@ class Variant {
         opac_);
   }
 
+  PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
+  GetRuntimePhysicalConstants() const {
+    return mpark::visit(
+        [=](const auto &opac) {
+          using PC = typename std::decay_t<decltype(opac)>::PC;
+          return singularity::GetRuntimePhysicalConstants(PC());
+        },
+        opac_);
+  }
+
   // Directional absorption coefficient with units of 1/length
   // Signature should be at least
   // rho, temp, nu, lambda
diff --git a/singularity-opac/photons/powerlaw_opacity_photons.hpp b/singularity-opac/photons/powerlaw_opacity_photons.hpp
index 7efe5e6..33a4c3f 100644
--- a/singularity-opac/photons/powerlaw_opacity_photons.hpp
+++ b/singularity-opac/photons/powerlaw_opacity_photons.hpp
@@ -30,6 +30,8 @@ namespace photons {
 template <typename pc = PhysicalConstantsCGS>
 class PowerLawOpacity {
  public:
+  using PC = pc;
+
   PowerLawOpacity() = default;
   PowerLawOpacity(const Real kappa0, const Real rho_exp, const Real temp_exp)
       : kappa0_(kappa0), rho_exp_(rho_exp), temp_exp_(temp_exp) {}
diff --git a/test/test_gray_opacities.cpp b/test/test_gray_opacities.cpp
index 464dbd1..aba0158 100644
--- a/test/test_gray_opacities.cpp
+++ b/test/test_gray_opacities.cpp
@@ -1,5 +1,5 @@
 // ======================================================================
-// © 2021. Triad National Security, LLC. All rights reserved.  This
+// © 2021-2024. Triad National Security, LLC. All rights reserved.  This
 // program was produced under U.S. Government contract
 // 89233218CNA000001 for Los Alamos National Laboratory (LANL), which
 // is operated by Triad National Security, LLC for the U.S.
@@ -54,9 +54,25 @@ TEST_CASE("Gray neutrino opacities", "[GrayNeutrinos]") {
     constexpr RadiationType type = RadiationType::NU_ELECTRON;
     constexpr Real nu = 1.25 * MeV2Hz; // 1 MeV
 
-    neutrinos::Gray opac_host(1);
+    neutrinos::Opacity opac_host = neutrinos::Gray(1.);
     neutrinos::Opacity opac = opac_host.GetOnDevice();
 
+    // Check constants from opacity
+    THEN("Check constants from mean opacity for consistency") {
+      auto constants = opac_host.GetRuntimePhysicalConstants();
+      int n_wrong = 0;
+      if (FractionalDifference(pc::eV, constants.eV) > EPS_TEST) {
+        n_wrong += 1;
+      }
+      if (FractionalDifference(pc::kb, constants.kb) > EPS_TEST) {
+        n_wrong += 1;
+      }
+      if (FractionalDifference(pc::h, constants.h) > EPS_TEST) {
+        n_wrong += 1;
+      }
+      REQUIRE(n_wrong == 0);
+    }
+
     THEN("The emissivity per nu omega is consistent with the emissity per nu") {
       int n_wrong_h = 0;
 #ifdef PORTABILITY_STRATEGY_KOKKOS
@@ -203,6 +219,23 @@ TEST_CASE("Gray photon opacities", "[GrayPhotons]") {
 
     photons::Opacity opac_host = photons::Gray(1);
     photons::Opacity opac = opac_host.GetOnDevice();
+
+    // Check constants from mean opacity
+    THEN("Check constants from mean opacity for consistency") {
+      auto constants = opac_host.GetRuntimePhysicalConstants();
+      int n_wrong = 0;
+      if (FractionalDifference(pc::eV, constants.eV) > EPS_TEST) {
+        n_wrong += 1;
+      }
+      if (FractionalDifference(pc::kb, constants.kb) > EPS_TEST) {
+        n_wrong += 1;
+      }
+      if (FractionalDifference(pc::h, constants.h) > EPS_TEST) {
+        n_wrong += 1;
+      }
+      REQUIRE(n_wrong == 0);
+    }
+
     THEN("The emissivity per nu omega is consistent with the emissity per nu") {
       int n_wrong_h = 0;
 #ifdef PORTABILITY_STRATEGY_KOKKOS
@@ -283,8 +316,7 @@ TEST_CASE("Gray photon opacities", "[GrayPhotons]") {
 
       Real *nu_bins = (Real *)PORTABLE_MALLOC(nbins * sizeof(Real));
       Real *temp_bins = (Real *)PORTABLE_MALLOC(ntemps * sizeof(Real));
-      DataBox loglin_bins(Spiner::AllocationTarget::Device, ntemps,
-                                  nbins);
+      DataBox loglin_bins(Spiner::AllocationTarget::Device, ntemps, nbins);
 
       portableFor(
           "set nu bins", 0, nbins, PORTABLE_LAMBDA(const int &i) {
diff --git a/test/test_mean_opacities.cpp b/test/test_mean_opacities.cpp
index b6ad91f..932d3d1 100644
--- a/test/test_mean_opacities.cpp
+++ b/test/test_mean_opacities.cpp
@@ -1,5 +1,5 @@
 // ======================================================================
-// © 2022. Triad National Security, LLC. All rights reserved.  This
+// © 2022-2024. Triad National Security, LLC. All rights reserved.  This
 // program was produced under U.S. Government contract
 // 89233218CNA000001 for Los Alamos National Laboratory (LANL), which
 // is operated by Triad National Security, LLC for the U.S.
@@ -46,7 +46,7 @@ using atomic_view = Kokkos::MemoryTraits<Kokkos::Atomic>;
 
 template <typename T>
 PORTABLE_INLINE_FUNCTION T FractionalDifference(const T &a, const T &b) {
-  return 2 * std::abs(b - a) / (std::abs(a) + std::abs(b) + 1e-20);
+  return 2 * std::abs(b - a) / (std::abs(a) + std::abs(b) + 1e-100);
 }
 constexpr Real EPS_TEST = 1e-3;
 template <typename T>
@@ -87,9 +87,8 @@ TEST_CASE("Mean neutrino opacities", "[MeanNeutrinos]") {
     neutrinos::Gray opac_host(kappa);
     neutrinos::Opacity opac = opac_host.GetOnDevice();
 
-    neutrinos::MeanOpacityCGS mean_opac_host(
+    neutrinos::MeanOpacity mean_opac_host = neutrinos::MeanOpacityBase(
         opac_host, lRhoMin, lRhoMax, NRho, lTMin, lTMax, NT, YeMin, YeMax, NYe);
-    // neutrinos::MeanOpacity mean_opac = mean_opac_host.GetOnDevice();
     auto mean_opac = mean_opac_host.GetOnDevice();
 
     THEN("The emissivity per nu omega is consistent with the emissity per nu") {
@@ -123,7 +122,7 @@ TEST_CASE("Mean neutrino opacities", "[MeanNeutrinos]") {
 #ifdef SPINER_USE_HDF
     THEN("We can save to disk and reload") {
       mean_opac.Save(grayname);
-      neutrinos::MeanOpacityCGS mean_opac_host_load(grayname);
+      neutrinos::MeanOpacity mean_opac_host_load(grayname);
       AND_THEN("The reloaded table matches the gray opacities") {
 
         auto mean_opac_load = mean_opac_host_load.GetOnDevice();
@@ -253,7 +252,6 @@ TEST_CASE("Mean neutrino scattering opacities", "[MeanNeutrinosS]") {
 
     neutrinos::MeanSOpacityCGS mean_opac_host(
         opac_host, lRhoMin, lRhoMax, NRho, lTMin, lTMax, NT, YeMin, YeMax, NYe);
-    // neutrinos::MeanOpacity mean_opac = mean_opac_host.GetOnDevice();
     auto mean_opac = mean_opac_host.GetOnDevice();
 
     THEN("The emissivity per nu omega is consistent with the emissity per nu") {
@@ -414,8 +412,8 @@ TEST_CASE("Mean photon opacities", "[MeanPhotons]") {
     photons::Gray opac_host(kappa);
     photons::Opacity opac = opac_host.GetOnDevice();
 
-    photons::MeanOpacityCGS mean_opac_host(opac_host, lRhoMin, lRhoMax, NRho,
-                                           lTMin, lTMax, NT);
+    photons::MeanOpacity mean_opac_host = photons::MeanOpacityBase(
+        opac_host, lRhoMin, lRhoMax, NRho, lTMin, lTMax, NT);
     auto mean_opac = mean_opac_host.GetOnDevice();
 
     THEN("The emissivity per nu omega is consistent with the emissity per nu") {
@@ -449,7 +447,8 @@ TEST_CASE("Mean photon opacities", "[MeanPhotons]") {
 #ifdef SPINER_USE_HDF
     THEN("We can save to disk and reload") {
       mean_opac.Save(grayname);
-      photons::MeanOpacityCGS mean_opac_host_load(grayname);
+      photons::MeanOpacity mean_opac_host_load =
+          photons::MeanOpacityBase(grayname);
       AND_THEN("The reloaded table matches the gray opacities") {
 
         auto mean_opac_load = mean_opac_host_load.GetOnDevice();
@@ -611,7 +610,8 @@ TEST_CASE("Mean photon scattering opacities", "[MeanPhotonS]") {
         int n_wrong = 0;
         portableReduce(
             "rebuilt table vs gray", 0, NRho, 0, NT, 0, 0,
-            PORTABLE_LAMBDA(const int iRho, const int iT, const int igarbage, int &accumulate) {
+            PORTABLE_LAMBDA(const int iRho, const int iT, const int igarbage,
+                            int &accumulate) {
               const Real lRho =
                   lRhoMin + (lRhoMax - lRhoMin) / (NRho - 1) * iRho;
               const Real rho = std::pow(10, lRho);