diff --git a/CHANGELOG.md b/CHANGELOG.md
index e2e855636a..6f04083d0b 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -5,7 +5,8 @@
 ### Added (new features/APIs/variables/...)
 
 ### Fixed (Repair bugs, etc)
-- [[OR437]](https://github.com/lanl/singularity-eos/pull/437) Fix segfault on HIP, clean up warnings, add strict sanitizer test
+- [[PR439]](https://github.com/lanl/singularity-eos/pull/439) Add mean atomic mass and number to EOS API
+- [[PR437]](https://github.com/lanl/singularity-eos/pull/437) Fix segfault on HIP, clean up warnings, add strict sanitizer test
 
 ### Changed (changing behavior/API/variables/...)
 
diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst
index 1087807b2b..5f3b914856 100644
--- a/doc/sphinx/src/models.rst
+++ b/doc/sphinx/src/models.rst
@@ -338,6 +338,38 @@ Note: sometimes temperatures are measured in eV for which the conversion is
 
 Sesame units are equivalent to the mm-mg-µs unit system.
 
+The ``MeanAtomicProperties`` struct
+------------------------------------
+
+Several analytic equations of state optionally accept mean atomic mass
+and number as physics properties. These are the average number of
+nucleons and protons in a constituent nucleus respectively. They are
+not necessarily integers, as a given material may be made up of
+multiple kinds of atom. For example, dry air contains both nitrogen
+and oxygen.
+
+The mean atomic mass and number are frequently carried in the
+container struct
+
+.. code-block:: cpp
+
+  struct MeanAtomicProperties {
+    Real Abar, Zbar;
+
+    // default is hydrogen
+    static constexpr Real DEFAULT_ABAR = 1.0;
+    static constexpr Real DEFAULT_ZBAR = 1.0;
+
+    PORTABLE_INLINE_FUNCTION
+    MeanAtomicProperties(Real Abar_, Real Zbar_) : Abar(Abar_), Zbar(Zbar_) {}
+    PORTABLE_INLINE_FUNCTION
+    MeanAtomicProperties() : Abar(DEFAULT_ABAR), Zbar(DEFAULT_ZBAR) {}
+  };
+
+which owns the atomic mass ``Abar`` and atomic number ``Zbar``. You
+may set these by constructing the struct or by setting the fields in a
+pre-constructed struct. The defaults are for hydrogen.
+
 Implemented EOS models
 ----------------------
 
@@ -429,6 +461,14 @@ these values are not set, they will be the same as those returned by the
 conditions are given, the return values of the :code:`ValuesAtReferenceState()`
 function will not be the same.
 
+Both constructors also optionally accept `MeanAtomicProperties` for
+the atomic mass and number as a final optional parameter, e.g.,
+
+.. code-block:: cpp
+
+   IdealGas(Real gm1, Real Cv, MeanAtomicProperties(Abar, Zbar));
+   IdealGas(Real gm1, Real Cv, Real EntropyT0, Real EntropyRho0, MeanAtomicProperties(Abar, Zbar));
+
 Stiffened Gas
 `````````````
 
@@ -486,6 +526,15 @@ these values are not set, they will be the same as those returned by the
 conditions are given, the return values of the :code:`ValuesAtReferenceState()`
 function will not be the same.
 
+Both constructors also optionally accept `MeanAtomicProperties` for
+the atomic mass and number as a final optional parameter, e.g.,
+
+.. code-block:: cpp
+
+   StiffGas(Real gm1, Real Cv, Real Pinf, Real q, MeanAtomicProperties(Abar, Zbar));
+   StiffGas(Real gm1, Real Cv, Real Pinf, Real q, Real qp, Real T0, Real P0,
+            EntropyRho0, MeanAtomicProperties(Abar, Zbar));
+
 Noble-Abel
 ``````````
 
@@ -544,6 +593,17 @@ these values are not set, they will be the same as those returned by the
 conditions are given, the return values of the :code:`ValuesAtReferenceState()`
 function will not be the same.
 
+Both constructors also optionally accept `MeanAtomicProperties` for
+the atomic mass and number as a final optional parameter, e.g.,
+
+.. code-block:: cpp
+
+   NobleAbel(Real gm1, Real Cv, Real b, Real q,
+             MeanAtomicProperties(Abar, Zbar));
+   NobleAbel(Real gm1, Real Cv, Real b, Real q, Real qp, Real T0, Real P0,
+             MeanAtomicProperties(Abar, Zbar));
+
+
 Carnahan-Starling
 `````````````````
 
@@ -630,6 +690,16 @@ these values are not set, they will be the same as those returned by the
 conditions are given, the return values of the :code:`ValuesAtReferenceState()`
 function will not be the same.
 
+Both constructors also optionally accept `MeanAtomicProperties` for
+the atomic mass and number as a final optional parameter, e.g.,
+
+.. code-block:: cpp
+
+    CarnahanStarling(Real gm1, Real Cv, Real b, Real q,
+                     MeanAtomicProperties(Abar, Zbar))
+    CarnahanStarling(Real gm1, Real Cv, Real b, Real q, Real qp, Real T0, Real P0,
+                     MeanAtomicProperties(Abar, Zbar))
+
 Gruneisen EOS
 `````````````
 
@@ -767,6 +837,9 @@ There is an overload of the ``Gruneisen`` class which computes
   Gruneisen(const Real C0, const Real s1, const Real s2, const Real s3, const Real G0,
             const Real b, const Real rho0, const Real T0, const Real P0, const Real Cv)
 
+Both constructors also optionally accept `MeanAtomicProperties` for
+the atomic mass and number as a final optional parameter.
+
 Extendended Vinet EOS
 `````````````````````
 
@@ -849,6 +922,9 @@ is :math:`S_0`, and ``expconsts`` is a pointer to the constant array of length
 :math:`d_2` to :math:`d_{40}`. Expansion coefficients not used should be set to
 0.0.
 
+This constructor also optionally accepts `MeanAtomicProperties` for
+the atomic mass and number as a final optional parameter.
+
 Mie-Gruneisen linear :math:`U_s`- :math:`u_p` EOS
 `````````````````````````````````````````````````
 
@@ -995,6 +1071,9 @@ where
 ``G0`` is :math:`\Gamma(\rho_0)`, ``Cv0`` is :math:`C_V`,
 ``E0`` is :math:`E_0`, and ``S0`` is :math:`S_0`. 
 
+This constructor also optionally accepts `MeanAtomicProperties` for
+the atomic mass and number as a final optional parameter.
+
 Mie-Gruneisen power expansion EOS
 `````````````````````````````````
 As we noted above, the assumption of a linear :math:`U_s`- :math:`u_p` relation is simply not valid at large compressions. At 
@@ -1075,6 +1154,8 @@ where
 :math:`K_0` to :math:`K_{40}`. Expansion coefficients not used should be set to
 :math:`0.0`.
 
+This constructor also optionally accepts `MeanAtomicProperties` for
+the atomic mass and number as a final optional parameter.
 
 
 JWL EOS
@@ -1127,6 +1208,9 @@ where ``A`` is :math:`A`, ``B`` is :math:`B`, ``R1`` is :math:`R_1`,
 ``R2`` is :math:`R_2`, ``w`` is :math:`w`, ``rho0`` is :math:`\rho_0`,
 and ``Cv`` is :math:`C_V`.
 
+This constructor also optionally accepts `MeanAtomicProperties` for
+the atomic mass and number as a final optional parameter.
+
 Davis EOS
 ```````````
 
@@ -1274,6 +1358,9 @@ where ``rho0`` is :math:`\rho_0`, ``e0`` is :math:`e_0`, ``P0`` is
 :math:`Z`, ``alpha`` is :math:`\alpha`, and ``Cv0`` is the specific
 heat capacity at the reference state.
 
+This constructor also optionally accepts `MeanAtomicProperties` for
+the atomic mass and number as a final optional parameter.
+
 Davis Products EOS
 '''''''''''''''''''
 
@@ -1381,6 +1468,9 @@ where ``a`` is :math:`a`, ``b`` is :math:`b`, ``k`` is :math:`k`,
 ``n`` is :math:`n`, ``vc`` is :math:`V_\mathrm{C}`, ``pc`` is
 :math:`P_\mathrm{C}`, ``Cv`` is :math:`C_{V,0}`.
 
+This constructor also optionally accepts `MeanAtomicProperties` for
+the atomic mass and number as a final optional parameter.
+
 Spiner EOS
 ````````````
 
@@ -1443,6 +1533,12 @@ of the material in the file, and ``reproducability_mode`` is a boolean
 which slightly changes how initial guesses for root finds are
 computed. The constructor for ``SpinerEOSDependsRhoSie`` is identical.
 
+.. note::
+
+  Mean atomic mass and number are loaded from input tables. The
+  ``SpinerEOS`` model does **not** support the
+  ``MeanAtomicProperties`` struct.
+
 ``sp5`` files and ``sesame2spiner``
 `````````````````````````````````````
 
@@ -1630,6 +1726,8 @@ values in expansion and compression.
 
 and similar expressions for :math:`b_2^*`.
 
+The SAP polynomial EOS also optionally accepts a
+``MeanAtomicProperties`` struct.
 
 
 Stellar Collapse EOS
@@ -1711,6 +1809,17 @@ where ``filename`` is the file containing the tabulated model,
 original `Stellar Collapse`_ format, and ``filter_bmod`` specifies
 whether or not to apply the above-described median filter.
 
+.. note::
+
+  The ``StellarCollapse`` EOS assumes nuclear statistical equilibrium
+  and as such mean atomic mass and number are state variables. As such
+  class does not accept the ``MeanAtomicProperties`` struct. The
+  ``MeanAtomicMassFromDensityTemperature`` and
+  ``MeanAtomicNumberFromDensityTemperature`` functions return the
+  relevant quantities for some thermodynamic state. The
+  ``MeanAtomicMass()`` and ``MeanAtomicNumber()`` functions raise an
+  error.
+
 ``StellarCollapse`` also provides 
 
 .. cpp:function:: void Save(const std::string &filename)
@@ -1796,6 +1905,14 @@ of the Helmholtz free energy (hence the name Helmholtz EOS). The free
 energy is pre-computed via integrals over the Fermi sphere and
 tabulated in a file provided from `Frank Timmes's website`_.
 
+.. note::
+
+  Since mean atomic mass and number are required inputs, the
+  ``MeanAtomicMassFromDensityTemperature`` and
+  ``MeanAtomicNumberFromDensityAndTemperature`` functions simply
+  return the input values. The ``MeanAtomicMass()`` and
+  ``MeanAtomicNumber`` functions produce an error.
+
 The table is a simple small ascii file. To ensure thermodyanic
 consistency, the table is interpolated using either biquintic or
 bicubic Hermite polynomials, which are sufficiently high order that
@@ -1897,6 +2014,13 @@ and :math:`splitCowan` uses the cold curve plus Cowan-nuclear model
 for ions and the final option ``linear_interp`` uses linear instead of 
 bilinear interpolation. 
 
+.. note::
+
+  Mean atomic mass and number are loaded from input tables. The
+  ``EOSPAC`` model does **not** support the ``MeanAtomicProperties``
+  struct.
+
+
 Note for performance reasons this EOS uses a slightly different vector API.
 See :ref:`EOSPAC Vector Functions <eospac_vector>` for more details.
 
diff --git a/doc/sphinx/src/using-eos.rst b/doc/sphinx/src/using-eos.rst
index 8810914adb..e9a2ecf848 100644
--- a/doc/sphinx/src/using-eos.rst
+++ b/doc/sphinx/src/using-eos.rst
@@ -1037,6 +1037,56 @@ specific internal energy in :math:`erg/g`.
 
 The function
 
+.. code-block:: cpp
+
+   template<typename Indexer_t Real*>
+   void MeanAtomicMassFromDensityTemperature(const Real rho, const Real T,
+                                             Indexer_t &&lambda = nullpter) const;
+
+returns the mean atomic mass (i.e., the number of nucleons) of a
+material given density in :math:`g/cm^3` and temperature in
+Kelvin. The reason this is allowed to vary with density and
+temperature is that some equations of state, such as the Stellar
+Collapse and Helmholtz equations of state encapsulate reactive flows
+where the average nucleus may depend on thermodynamic variables. For
+most materials, however, this is not the case and a convenience
+function that drops the dependence is available:
+
+.. code-block:: cpp
+
+   Real MeanAtomicMass() const;
+
+The function
+
+.. code-block:: cpp
+
+   template<typename Indexer_t Real*>
+   void MeanAtomicNumberFromDensityTemperature(const Real rho, const Real T,
+                                               Indexer_t &&lambda = nullpter) const;
+
+returns the mean atomic number (i.e., the number of protons in the
+nucleus) of a material given density in :math:`g/cm^3` and temperature
+in Kelvin. The reason this is allowed to vary with density and
+temperature is that some equations of state, such as the Stellar
+Collapse and Helmholtz equations of state encapsulate reactive flows
+where the average nucleus may depend on thermodynamic variables. For
+most materials, however, this is not the case and a convenience
+function that drops the dependence is available:
+
+.. code-block:: cpp
+
+   Real MeanAtomicNumber() const;
+
+
+
+.. warning::
+
+  For materials where the mean atomic mass and number **do** vary with
+  density and temperature, the convenience call without this
+  dependence will produce an error.
+
+The function
+
 .. code-block:: cpp
 
    template <typename Indexer_t = Real*>
diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp
index 83d0bf03a9..8122925f2a 100644
--- a/singularity-eos/eos/eos_base.hpp
+++ b/singularity-eos/eos/eos_base.hpp
@@ -84,6 +84,16 @@ char *StrCat(char *destination, const char *source) {
   using EosBase<EOSDERIVED>::GibbsFreeEnergyFromDensityTemperature;                      \
   using EosBase<EOSDERIVED>::GibbsFreeEnergyFromDensityInternalEnergy;
 
+// This macro adds these methods to a derived class. Due to scope,
+// these can't be implemented in the base class, unless we make
+// _AZbar public. Not all EOS's may want these default functions
+// TODO(JMM): Should we go the alternate route and make _AZbar public?
+#define SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(_AZbar)                                     \
+  PORTABLE_INLINE_FUNCTION                                                               \
+  Real MeanAtomicMass() const { return _AZbar.Abar; }                                    \
+  PORTABLE_INLINE_FUNCTION                                                               \
+  Real MeanAtomicNumber() const { return _AZbar.Zbar; }
+
 // This macro adds several methods that most modifiers will
 // want. Not ALL modifiers will want these methods as written here,
 // so use this macro with care.
@@ -101,7 +111,11 @@ char *StrCat(char *destination, const char *source) {
   }                                                                                      \
   constexpr bool AllDynamicMemoryIsShareable() const {                                   \
     return t_.AllDynamicMemoryIsShareable();                                             \
-  }
+  }                                                                                      \
+  PORTABLE_INLINE_FUNCTION                                                               \
+  Real MeanAtomicMass() const { return t_.MeanAtomicMass(); }                            \
+  PORTABLE_INLINE_FUNCTION                                                               \
+  Real MeanAtomicNumber() const { return t_.MeanAtomicNumber(); }
 
 class Factor {
   Real value_ = 1.0;
@@ -132,6 +146,34 @@ struct Transform {
   Factor x, y, f;
 };
 
+/*
+  This is a utility struct used to bundle mean atomic
+  mass/number. Used in the default implementations of MeanAtomicMass
+  and MeanAtomicNumber provided by the base class.
+ */
+struct MeanAtomicProperties {
+  Real Abar, Zbar;
+
+  // default is hydrogen
+  static constexpr Real DEFAULT_ABAR = 1.0;
+  static constexpr Real DEFAULT_ZBAR = 1.0;
+
+  PORTABLE_INLINE_FUNCTION
+  MeanAtomicProperties(Real Abar_, Real Zbar_) : Abar(Abar_), Zbar(Zbar_) {}
+  PORTABLE_INLINE_FUNCTION
+  MeanAtomicProperties() : Abar(DEFAULT_ABAR), Zbar(DEFAULT_ZBAR) {}
+  PORTABLE_INLINE_FUNCTION
+  void CheckParams() const {
+    PORTABLE_ALWAYS_REQUIRE(Abar > 0, "Positive mean atomic mass");
+    PORTABLE_ALWAYS_REQUIRE(Zbar > 0, "Positive mean atomic number");
+  }
+  PORTABLE_INLINE_FUNCTION
+  void PrintParams() const {
+    printf("      Abar  = %g\n", Abar);
+    printf("      Zbar  = %g\n", Zbar);
+  }
+};
+
 /*
 This is a CRTP that allows for static inheritance so that default behavior for
 various member functions can be defined.
@@ -739,6 +781,29 @@ class EosBase {
   PORTABLE_INLINE_FUNCTION
   Real RhoPmin(const Real temp) const { return 0.0; }
 
+  // JMM: EOS's which encapsulate a mix or reactions may wish to vary
+  // this.  For example, Helmholtz and StellarCollapse. This isn't the
+  // default, so by default the base class provides a specialization.
+  // for models where density and temperature are required, the EOS
+  // developer is in charge of either throwing an error or choosing
+  // reasonable defaults.
+  // TODO(JMM): Should we provide vector implementations if we depend
+  // on rho, T, etc?
+  template <typename Indexer_t = Real *>
+  PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature(
+      const Real rho, const Real T,
+      Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
+    CRTP copy = *(static_cast<CRTP const *>(this));
+    return copy.MeanAtomicMass();
+  }
+  template <typename Indexer_t = Real *>
+  PORTABLE_INLINE_FUNCTION Real MeanAtomicNumberFromDensityTemperature(
+      const Real rho, const Real T,
+      Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
+    CRTP copy = *(static_cast<CRTP const *>(this));
+    return copy.MeanAtomicNumber();
+  }
+
   // Default entropy behavior is to cause an error
   PORTABLE_FORCEINLINE_FUNCTION
   void EntropyIsNotEnabled(const char *eosname) const {
diff --git a/singularity-eos/eos/eos_carnahan_starling.hpp b/singularity-eos/eos/eos_carnahan_starling.hpp
index e77e844332..ba07356441 100644
--- a/singularity-eos/eos/eos_carnahan_starling.hpp
+++ b/singularity-eos/eos/eos_carnahan_starling.hpp
@@ -38,7 +38,9 @@ using namespace eos_base;
 class CarnahanStarling : public EosBase<CarnahanStarling> {
  public:
   CarnahanStarling() = default;
-  PORTABLE_INLINE_FUNCTION CarnahanStarling(Real gm1, Real Cv, Real bb, Real qq)
+  PORTABLE_INLINE_FUNCTION
+  CarnahanStarling(Real gm1, Real Cv, Real bb, Real qq,
+                   const MeanAtomicProperties &AZbar = MeanAtomicProperties())
       : _Cv(Cv), _gm1(gm1), _bb(bb), _qq(qq), _T0(ROOM_TEMPERATURE),
         _P0(ATMOSPHERIC_PRESSURE), _qp(0.0),
         _rho0(DensityFromPressureTemperature(_P0, _T0)), _vol0(robust::ratio(1.0, _rho0)),
@@ -46,18 +48,19 @@ class CarnahanStarling : public EosBase<CarnahanStarling> {
         _bmod0(_rho0 * Cv * _T0 *
                (PartialRhoZedFromDensity(_rho0) +
                 ZedFromDensity(_rho0) * ZedFromDensity(_rho0) * gm1)),
-        _dpde0(_rho0 * ZedFromDensity(_rho0) * gm1) {
+        _dpde0(_rho0 * ZedFromDensity(_rho0) * gm1), _AZbar(AZbar) {
     CheckParams();
   }
-  PORTABLE_INLINE_FUNCTION CarnahanStarling(Real gm1, Real Cv, Real bb, Real qq, Real qp,
-                                            Real T0, Real P0)
+  PORTABLE_INLINE_FUNCTION
+  CarnahanStarling(Real gm1, Real Cv, Real bb, Real qq, Real qp, Real T0, Real P0,
+                   const MeanAtomicProperties &AZbar = MeanAtomicProperties())
       : _Cv(Cv), _gm1(gm1), _bb(bb), _qq(qq), _T0(T0), _P0(P0), _qp(qp),
         _rho0(DensityFromPressureTemperature(P0, T0)), _vol0(robust::ratio(1.0, _rho0)),
         _sie0(Cv * T0 + qq),
         _bmod0(_rho0 * Cv * T0 *
                (PartialRhoZedFromDensity(_rho0) +
                 ZedFromDensity(_rho0) * ZedFromDensity(_rho0) * gm1)),
-        _dpde0(_rho0 * ZedFromDensity(_rho0) * gm1) {
+        _dpde0(_rho0 * ZedFromDensity(_rho0) * gm1), _AZbar(AZbar) {
     CheckParams();
   }
   CarnahanStarling GetOnDevice() { return *this; }
@@ -71,6 +74,7 @@ class CarnahanStarling : public EosBase<CarnahanStarling> {
     PORTABLE_ALWAYS_REQUIRE(_Cv >= 0, "Heat capacity must be positive");
     PORTABLE_ALWAYS_REQUIRE(_gm1 >= 0, "Gruneisen parameter must be positive");
     PORTABLE_ALWAYS_REQUIRE(_bb >= 0, "Covolume must be positive");
+    _AZbar.CheckParams();
   }
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION Real ZedFromDensity(
@@ -221,6 +225,8 @@ class CarnahanStarling : public EosBase<CarnahanStarling> {
     bmod = _bmod0;
     dpde = _dpde0;
   }
+  // Default accessors for mean atomic mass/number
+  SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(_AZbar)
   // Generic functions provided by the base class. These contain e.g. the vector
   // overloads that use the scalar versions declared here
   SG_ADD_BASE_CLASS_USINGS(CarnahanStarling)
@@ -235,6 +241,7 @@ class CarnahanStarling : public EosBase<CarnahanStarling> {
     printf("Carnahan-Starling Parameters:\nGamma = %g\nCv    = %g\nb     = %g\nq     = "
            "%g\n",
            _gm1 + 1.0, _Cv, _bb, _qq);
+    _AZbar.PrintParams();
   }
   template <typename Indexer_t>
   PORTABLE_INLINE_FUNCTION void
@@ -253,6 +260,7 @@ class CarnahanStarling : public EosBase<CarnahanStarling> {
   Real _T0, _P0, _qp;
   // reference values
   Real _rho0, _vol0, _sie0, _bmod0, _dpde0;
+  MeanAtomicProperties _AZbar;
   // static constexpr const Real _T0 = ROOM_TEMPERATURE;
   // static constexpr const Real _P0 = ATMOSPHERIC_PRESSURE;
   static constexpr const unsigned long _preferred_input =
diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp
index 2bfd96ad76..ff7091c969 100644
--- a/singularity-eos/eos/eos_davis.hpp
+++ b/singularity-eos/eos/eos_davis.hpp
@@ -35,9 +35,10 @@ class DavisReactants : public EosBase<DavisReactants> {
   PORTABLE_INLINE_FUNCTION
   DavisReactants(const Real rho0, const Real e0, const Real P0, const Real T0,
                  const Real A, const Real B, const Real C, const Real G0, const Real Z,
-                 const Real alpha, const Real Cv0)
+                 const Real alpha, const Real Cv0,
+                 const MeanAtomicProperties &AZbar = MeanAtomicProperties())
       : _rho0(rho0), _e0(e0), _P0(P0), _T0(T0), _A(A), _B(B), _C(C), _G0(G0), _Z(Z),
-        _alpha(alpha), _Cv0(Cv0) {
+        _alpha(alpha), _Cv0(Cv0), _AZbar(AZbar) {
     CheckParams();
   }
   DavisReactants GetOnDevice() { return *this; }
@@ -45,6 +46,7 @@ class DavisReactants : public EosBase<DavisReactants> {
   void CheckParams() const {
     PORTABLE_REQUIRE(_rho0 >= 0, "Density must be positive");
     PORTABLE_REQUIRE(_T0 >= 0, "Temperature must be positive");
+    _AZbar.CheckParams();
   }
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION Real TemperatureFromDensityInternalEnergy(
@@ -160,6 +162,8 @@ class DavisReactants : public EosBase<DavisReactants> {
   PORTABLE_INLINE_FUNCTION void
   DensityEnergyFromPressureTemperature(const Real press, const Real temp,
                                        Indexer_t &&lambda, Real &rho, Real &sie) const;
+  // Default accessors for mean atomic mass/number
+  SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(_AZbar)
   // Generic functions provided by the base class. These contain e.g. the vector
   // overloads that use the scalar versions declared here
   SG_ADD_BASE_CLASS_USINGS(DavisReactants)
@@ -175,6 +179,7 @@ class DavisReactants : public EosBase<DavisReactants> {
     printf("%srho0:%e e0:%e P0:%e\nT0:%e A:%e B:%e\nC:%e G0:%e Z:%e\nalpha:%e "
            "Cv0:%e\n",
            s1, _rho0, _e0, _P0, _T0, _A, _B, _C, _G0, _Z, _alpha, _Cv0);
+    _AZbar.PrintParams();
   }
   void inline Finalize() {}
   static std::string EosType() { return std::string("DavisReactants"); }
@@ -183,6 +188,7 @@ class DavisReactants : public EosBase<DavisReactants> {
  private:
   static constexpr Real onethird = 1.0 / 3.0;
   Real _rho0, _e0, _P0, _T0, _A, _B, _C, _G0, _Z, _alpha, _Cv0;
+  MeanAtomicProperties _AZbar;
   // static constexpr const char _eos_type[] = "DavisReactants";
   static constexpr unsigned long _preferred_input =
       thermalqs::density | thermalqs::specific_internal_energy;
@@ -198,10 +204,12 @@ class DavisProducts : public EosBase<DavisProducts> {
   DavisProducts() = default;
   PORTABLE_INLINE_FUNCTION
   DavisProducts(const Real a, const Real b, const Real k, const Real n, const Real vc,
-                const Real pc, const Real Cv)
-      : _a(a), _b(b), _k(k), _n(n), _vc(vc), _pc(pc), _Cv(Cv) {}
+                const Real pc, const Real Cv,
+                const MeanAtomicProperties &AZbar = MeanAtomicProperties())
+      : _a(a), _b(b), _k(k), _n(n), _vc(vc), _pc(pc), _Cv(Cv), _AZbar(AZbar) {}
   PORTABLE_INLINE_FUNCTION
   void CheckParams() const {
+    _AZbar.CheckParams();
     // TODO(JMM): Stub.
   }
   DavisProducts GetOnDevice() { return *this; }
@@ -299,6 +307,8 @@ class DavisProducts : public EosBase<DavisProducts> {
   ValuesAtReferenceState(Real &rho, Real &temp, Real &sie, Real &press, Real &cv,
                          Real &bmod, Real &dpde, Real &dvdt,
                          Indexer_t &&lambda = static_cast<Real *>(nullptr)) const;
+  // Default accessors for mean atomic mass/number
+  SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(_AZbar)
   // Generic functions provided by the base class. These contain e.g. the vector
   // overloads that use the scalar versions declared here
   SG_ADD_BASE_CLASS_USINGS(DavisProducts)
@@ -313,6 +323,7 @@ class DavisProducts : public EosBase<DavisProducts> {
     static constexpr char s1[]{"DavisProducts Params: "};
     printf("%sa:%e b:%e k:%e\nn:%e vc:%e pc:%e\nCv:%e \n", s1, _a, _b, _k, _n, _vc, _pc,
            _Cv);
+    _AZbar.PrintParams();
   }
   inline void Finalize() {}
   static std::string EosType() { return std::string("DavisProducts"); }
@@ -321,6 +332,7 @@ class DavisProducts : public EosBase<DavisProducts> {
  private:
   static constexpr Real onethird = 1.0 / 3.0;
   Real _a, _b, _k, _n, _vc, _pc, _Cv;
+  MeanAtomicProperties _AZbar;
   // static constexpr const char _eos_type[] = "DavisProducts";
   static constexpr const unsigned long _preferred_input =
       thermalqs::density | thermalqs::specific_internal_energy;
diff --git a/singularity-eos/eos/eos_eospac.hpp b/singularity-eos/eos/eos_eospac.hpp
index 1f08a65873..763cb53695 100644
--- a/singularity-eos/eos/eos_eospac.hpp
+++ b/singularity-eos/eos/eos_eospac.hpp
@@ -141,6 +141,7 @@ class EOSPAC : public EosBase<EOSPAC> {
     // TODO(JMM): More validation checks?
     PORTABLE_ALWAYS_REQUIRE(rho_min_ >= 0, "Non-negative minimum density");
     PORTABLE_ALWAYS_REQUIRE(temp_min_ >= 0, "Non-negative minimum temperature");
+    AZbar_.CheckParams();
   }
   inline EOSPAC GetOnDevice() { return *this; }
 
@@ -457,6 +458,8 @@ class EOSPAC : public EosBase<EOSPAC> {
   using EosBase<EOSPAC>::FillEos;
   using EosBase<EOSPAC>::EntropyIsNotEnabled;
 
+  SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(AZbar_)
+
   // EOSPAC vector implementations
   template <typename LambdaIndexer>
   inline void
@@ -1097,6 +1100,7 @@ class EOSPAC : public EosBase<EOSPAC> {
   static std::string EosPyType() { return EosType(); }
   PORTABLE_INLINE_FUNCTION void PrintParams() const {
     printf("EOSPAC parameters:\nmatid = %i\n", matid_);
+    AZbar_.PrintParams();
   }
   PORTABLE_FORCEINLINE_FUNCTION Real MinimumDensity() const { return rho_min_; }
   PORTABLE_FORCEINLINE_FUNCTION Real MinimumTemperature() const { return temp_min_; }
@@ -1127,6 +1131,7 @@ class EOSPAC : public EosBase<EOSPAC> {
   // TODO(JMM): Is the fact that EOS_INTEGER isn't a size_t a
   // problem? Could it ever realistically overflow?
   EOS_INTEGER shared_size_, packed_size_;
+  MeanAtomicProperties AZbar_;
 
   static inline std::map<std::string, unsigned int> &scratch_nbuffers() {
     static std::map<std::string, unsigned int> nbuffers = {
@@ -1194,6 +1199,9 @@ inline EOSPAC::EOSPAC(const int matid, bool invert_at_setup, Real insert_data,
   rho_ref_ = m.normalDensity;
   rho_min_ = m.rhoMin;
   temp_min_ = m.TMin;
+  // use std::max to hydrogen, in case of bad table
+  AZbar_.Abar = std::max(1.0, m.meanAtomicMass);
+  AZbar_.Zbar = std::max(1.0, m.meanAtomicNumber);
 
   EOS_REAL R[1] = {rho_ref_};
   EOS_REAL T[1] = {temperatureToSesame(temp_ref_)};
diff --git a/singularity-eos/eos/eos_gruneisen.hpp b/singularity-eos/eos/eos_gruneisen.hpp
index 93570a37e7..a2892e9a2a 100644
--- a/singularity-eos/eos/eos_gruneisen.hpp
+++ b/singularity-eos/eos/eos_gruneisen.hpp
@@ -39,9 +39,10 @@ class Gruneisen : public EosBase<Gruneisen> {
   PORTABLE_INLINE_FUNCTION
   Gruneisen(const Real C0, const Real s1, const Real s2, const Real s3, const Real G0,
             const Real b, const Real rho0, const Real T0, const Real P0, const Real Cv,
-            const Real rho_max)
+            const Real rho_max,
+            const MeanAtomicProperties &AZbar = MeanAtomicProperties())
       : _C0(C0), _s1(s1), _s2(s2), _s3(s3), _G0(G0), _b(b), _rho0(rho0), _T0(T0), _P0(P0),
-        _Cv(Cv), _rho_max(rho_max) {
+        _Cv(Cv), _rho_max(rho_max), _AZbar(AZbar) {
     // Warn user when provided rho_max is greater than the computed rho_max
 #ifndef NDEBUG
     const Real computed_rho_max = ComputeRhoMax(s1, s2, s3, rho0);
@@ -58,9 +59,11 @@ class Gruneisen : public EosBase<Gruneisen> {
   // Constructor when rho_max isn't specified automatically determines _rho_max
   PORTABLE_INLINE_FUNCTION
   Gruneisen(const Real C0, const Real s1, const Real s2, const Real s3, const Real G0,
-            const Real b, const Real rho0, const Real T0, const Real P0, const Real Cv)
+            const Real b, const Real rho0, const Real T0, const Real P0, const Real Cv,
+            const MeanAtomicProperties &AZbar = MeanAtomicProperties())
       : _C0(C0), _s1(s1), _s2(s2), _s3(s3), _G0(G0), _b(b), _rho0(rho0), _T0(T0), _P0(P0),
-        _Cv(Cv), _rho_max(RHOMAX_SAFETY * ComputeRhoMax(s1, s2, s3, rho0)) {}
+        _Cv(Cv), _rho_max(RHOMAX_SAFETY * ComputeRhoMax(s1, s2, s3, rho0)),
+        _AZbar(AZbar) {}
   static PORTABLE_INLINE_FUNCTION Real ComputeRhoMax(const Real s1, const Real s2,
                                                      const Real s3, const Real rho0);
   PORTABLE_INLINE_FUNCTION
@@ -71,6 +74,7 @@ class Gruneisen : public EosBase<Gruneisen> {
     PORTABLE_ALWAYS_REQUIRE(_Cv >= 0, "Non-negative heat capacity required");
     PORTABLE_ALWAYS_REQUIRE(_rho_max > _rho0,
                             "Maximum density must be greater than reference");
+    _AZbar.CheckParams();
   }
   PORTABLE_INLINE_FUNCTION Real
   MaxStableDensityAtTemperature(const Real temperature) const;
@@ -150,6 +154,7 @@ class Gruneisen : public EosBase<Gruneisen> {
                          Indexer_t &&lambda = static_cast<Real *>(nullptr)) const;
   // Generic functions provided by the base class. These contain e.g. the vector
   // overloads that use the scalar versions declared here
+  SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(_AZbar)
   SG_ADD_BASE_CLASS_USINGS(Gruneisen)
   PORTABLE_INLINE_FUNCTION
   int nlambda() const noexcept { return 0; }
@@ -163,6 +168,7 @@ class Gruneisen : public EosBase<Gruneisen> {
     printf("%s C0:%e s1:%e s2:%e s3:%e\n  G0:%e b:%e rho0:%e T0:%e\n  P0:%eCv:%e "
            "rho_max:%e\n",
            s1, _C0, _s1, _s2, _s3, _G0, _b, _rho0, _T0, _P0, _Cv, _rho_max);
+    _AZbar.PrintParams();
   }
   template <typename Indexer_t>
   PORTABLE_INLINE_FUNCTION void
@@ -174,6 +180,7 @@ class Gruneisen : public EosBase<Gruneisen> {
 
  private:
   Real _C0, _s1, _s2, _s3, _G0, _b, _rho0, _T0, _P0, _Cv, _rho_max;
+  MeanAtomicProperties _AZbar;
   // static constexpr const char _eos_type[] = {"Gruneisen"};
   PORTABLE_INLINE_FUNCTION
   Real Gamma(const Real rho_in) const {
diff --git a/singularity-eos/eos/eos_helmholtz.hpp b/singularity-eos/eos/eos_helmholtz.hpp
index e5fbb0be78..b6ade5ac19 100644
--- a/singularity-eos/eos/eos_helmholtz.hpp
+++ b/singularity-eos/eos/eos_helmholtz.hpp
@@ -640,6 +640,31 @@ class Helmholtz : public EosBase<Helmholtz> {
     return gamma3 - 1.0;
   }
 
+  PORTABLE_INLINE_FUNCTION
+  Real MeanAtomicMass() const {
+    PORTABLE_THROW_OR_ABORT("For Helmholtz EOS, mean atomic mass is an input!\n");
+    return 1.0;
+  }
+  PORTABLE_INLINE_FUNCTION
+  Real MeanAtomicNumber() const {
+    PORTABLE_THROW_OR_ABORT("For Helmholtz EOS, mean atomic number is an input!\n");
+    return 1.0;
+  }
+  template <typename Indexer_t = Real *>
+  PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature(
+      const Real rho, const Real T,
+      Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
+    using namespace HelmUtils;
+    return lambda[Lambda::Abar];
+  }
+  template <typename Indexer_t = Real *>
+  PORTABLE_INLINE_FUNCTION Real MeanAtomicNumberFromDensityTemperature(
+      const Real rho, const Real T,
+      Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
+    using namespace HelmUtils;
+    return lambda[Lambda::Zbar];
+  }
+
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION void
   FillEos(Real &rho, Real &temp, Real &energy, Real &press, Real &cv, Real &bmod,
diff --git a/singularity-eos/eos/eos_ideal.hpp b/singularity-eos/eos/eos_ideal.hpp
index d395928404..02a0514674 100644
--- a/singularity-eos/eos/eos_ideal.hpp
+++ b/singularity-eos/eos/eos_ideal.hpp
@@ -38,16 +38,20 @@ using namespace eos_base;
 class IdealGas : public EosBase<IdealGas> {
  public:
   IdealGas() = default;
-  PORTABLE_INLINE_FUNCTION IdealGas(Real gm1, Real Cv)
+  PORTABLE_INLINE_FUNCTION
+  IdealGas(Real gm1, Real Cv, const MeanAtomicProperties &AZbar = MeanAtomicProperties())
       : _Cv(Cv), _gm1(gm1), _rho0(_P0 / (_gm1 * _Cv * _T0)), _sie0(_Cv * _T0),
         _bmod0((_gm1 + 1) * _gm1 * _rho0 * _Cv * _T0), _dpde0(_gm1 * _rho0),
-        _dvdt0(1. / (_rho0 * _T0)), _EntropyT0(_T0), _EntropyRho0(_rho0) {
+        _dvdt0(1. / (_rho0 * _T0)), _EntropyT0(_T0), _EntropyRho0(_rho0), _AZbar(AZbar) {
     CheckParams();
   }
-  PORTABLE_INLINE_FUNCTION IdealGas(Real gm1, Real Cv, Real EntropyT0, Real EntropyRho0)
+  PORTABLE_INLINE_FUNCTION
+  IdealGas(Real gm1, Real Cv, Real EntropyT0, Real EntropyRho0,
+           const MeanAtomicProperties &AZbar = MeanAtomicProperties())
       : _Cv(Cv), _gm1(gm1), _rho0(_P0 / (_gm1 * _Cv * _T0)), _sie0(_Cv * _T0),
         _bmod0((_gm1 + 1) * _gm1 * _rho0 * _Cv * _T0), _dpde0(_gm1 * _rho0),
-        _dvdt0(1. / (_rho0 * _T0)), _EntropyT0(EntropyT0), _EntropyRho0(EntropyRho0) {
+        _dvdt0(1. / (_rho0 * _T0)), _EntropyT0(EntropyT0), _EntropyRho0(EntropyRho0),
+        _AZbar(AZbar) {
     CheckParams();
   }
 
@@ -67,6 +71,7 @@ class IdealGas : public EosBase<IdealGas> {
                             "Entropy reference temperature must be positive");
     PORTABLE_ALWAYS_REQUIRE(_EntropyRho0 >= 0,
                             "Entropy reference density must be positive");
+    _AZbar.CheckParams();
     return;
   }
   template <typename Indexer_t = Real *>
@@ -148,6 +153,7 @@ class IdealGas : public EosBase<IdealGas> {
   FillEos(Real &rho, Real &temp, Real &energy, Real &press, Real &cv, Real &bmod,
           const unsigned long output,
           Indexer_t &&lambda = static_cast<Real *>(nullptr)) const;
+
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION void
   ValuesAtReferenceState(Real &rho, Real &temp, Real &sie, Real &press, Real &cv,
@@ -165,6 +171,7 @@ class IdealGas : public EosBase<IdealGas> {
   }
   // Generic functions provided by the base class. These contain e.g. the vector
   // overloads that use the scalar versions declared here
+  SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(_AZbar)
   SG_ADD_BASE_CLASS_USINGS(IdealGas)
   PORTABLE_INLINE_FUNCTION
   int nlambda() const noexcept { return 0; }
@@ -175,6 +182,7 @@ class IdealGas : public EosBase<IdealGas> {
   static inline unsigned long max_scratch_size(unsigned int nelements) { return 0; }
   PORTABLE_INLINE_FUNCTION void PrintParams() const {
     printf("Ideal Gas Parameters:\nGamma = %g\nCv    = %g\n", _gm1 + 1.0, _Cv);
+    _AZbar.PrintParams();
   }
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION void
@@ -198,6 +206,8 @@ class IdealGas : public EosBase<IdealGas> {
       thermalqs::density | thermalqs::specific_internal_energy;
   // optional entropy reference state variables
   Real _EntropyT0, _EntropyRho0;
+  // optional mean atomic mass and number
+  MeanAtomicProperties _AZbar;
 };
 
 template <typename Indexer_t>
diff --git a/singularity-eos/eos/eos_jwl.hpp b/singularity-eos/eos/eos_jwl.hpp
index 1931541630..31cf5477b6 100644
--- a/singularity-eos/eos/eos_jwl.hpp
+++ b/singularity-eos/eos/eos_jwl.hpp
@@ -40,17 +40,20 @@ class JWL : public EosBase<JWL> {
  public:
   JWL() = default;
   PORTABLE_INLINE_FUNCTION JWL(const Real A, const Real B, const Real R1, const Real R2,
-                               const Real w, const Real rho0, const Real Cv)
+                               const Real w, const Real rho0, const Real Cv,
+                               const MeanAtomicProperties &AZbar = MeanAtomicProperties())
       : _A(A), _B(B), _R1(R1), _R2(R2), _w(w), _rho0(rho0), _Cv(Cv),
-        _c1(robust::ratio(A, rho0 * R1)), _c2(robust::ratio(B, rho0 * R2)) {
-    assert(R1 > 0.0);
-    assert(R2 > 0.0);
+        _c1(robust::ratio(A, rho0 * R1)), _c2(robust::ratio(B, rho0 * R2)),
+        _AZbar(AZbar) {
+    PORTABLE_REQUIRE(R1 > 0.0, "Ratio of rho0 to A > 0");
+    PORTABLE_REQUIRE(R2 > 0.0, "Ratio of rho0 to B > 0");
     CheckParams();
   }
   PORTABLE_INLINE_FUNCTION void CheckParams() const {
     PORTABLE_ALWAYS_REQUIRE(_w > 0.0, "w > 0");
     PORTABLE_ALWAYS_REQUIRE(_rho0 > 0.0, "Positive reference density");
     PORTABLE_ALWAYS_REQUIRE(_Cv > 0.0, "Positive specific heat");
+    _AZbar.CheckParams();
   }
   JWL GetOnDevice() { return *this; }
   template <typename Indexer_t = Real *>
@@ -111,6 +114,7 @@ class JWL : public EosBase<JWL> {
           Indexer_t &&lambda = static_cast<Real *>(nullptr)) const;
   // Generic functions provided by the base class. These contain e.g. the vector
   // overloads that use the scalar versions declared here
+  SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(_AZbar)
   SG_ADD_BASE_CLASS_USINGS(JWL)
   PORTABLE_INLINE_FUNCTION
   int nlambda() const noexcept { return 0; }
@@ -128,6 +132,7 @@ class JWL : public EosBase<JWL> {
     static constexpr char s1[]{"JWL Params: "};
     printf("%sA:%e B:%e R1: %e\nR2:%e w:%e rho0:%e\nCv:%e\n", s1, _A, _B, _R1, _R2, _w,
            _rho0, _Cv);
+    _AZbar.PrintParams();
   }
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION void
@@ -139,6 +144,7 @@ class JWL : public EosBase<JWL> {
 
  private:
   Real _A, _B, _R1, _R2, _w, _rho0, _Cv, _c1, _c2;
+  MeanAtomicProperties _AZbar;
   PORTABLE_INLINE_FUNCTION Real ReferenceEnergy(const Real rho) const;
   PORTABLE_INLINE_FUNCTION Real ReferencePressure(const Real rho) const;
   // static constexpr const char _eos_type[] = "JWL";
diff --git a/singularity-eos/eos/eos_mgusup.hpp b/singularity-eos/eos/eos_mgusup.hpp
index 409e596fe1..c315c1619a 100644
--- a/singularity-eos/eos/eos_mgusup.hpp
+++ b/singularity-eos/eos/eos_mgusup.hpp
@@ -38,8 +38,10 @@ class MGUsup : public EosBase<MGUsup> {
   MGUsup() = default;
   // Constructor
   MGUsup(const Real rho0, const Real T0, const Real Cs, const Real s, const Real G0,
-         const Real Cv0, const Real E0, const Real S0)
-      : _rho0(rho0), _T0(T0), _Cs(Cs), _s(s), _G0(G0), _Cv0(Cv0), _E0(E0), _S0(S0) {
+         const Real Cv0, const Real E0, const Real S0,
+         const MeanAtomicProperties &AZbar = MeanAtomicProperties())
+      : _rho0(rho0), _T0(T0), _Cs(Cs), _s(s), _G0(G0), _Cv0(Cv0), _E0(E0), _S0(S0),
+        _AZbar(AZbar) {
     CheckParams();
   }
   PORTABLE_INLINE_FUNCTION void CheckParams() const;
@@ -130,6 +132,7 @@ class MGUsup : public EosBase<MGUsup> {
                          Indexer_t &&lambda = static_cast<Real *>(nullptr)) const;
   // Generic functions provided by the base class. These contain e.g. the vector
   // overloads that use the scalar versions declared here
+  SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(_AZbar)
   SG_ADD_BASE_CLASS_USINGS(MGUsup)
   PORTABLE_INLINE_FUNCTION
   int nlambda() const noexcept { return 0; }
@@ -142,6 +145,7 @@ class MGUsup : public EosBase<MGUsup> {
     static constexpr char st[]{"MGUsup Params: "};
     printf("%s rho0:%e T0:%e Cs:%e s:%e\n  G0:%e Cv0:%e E0:%e S0:%e\n", st, _rho0, _T0,
            _Cs, _s, _G0, _Cv0, _E0, _S0);
+    _AZbar.PrintParams();
     printf("\n\n");
   }
   // Density/Energy from P/T not unique, if used will give error
@@ -157,6 +161,7 @@ class MGUsup : public EosBase<MGUsup> {
   static constexpr const unsigned long _preferred_input =
       thermalqs::density | thermalqs::specific_internal_energy;
   Real _rho0, _T0, _Cs, _s, _G0, _Cv0, _E0, _S0;
+  MeanAtomicProperties _AZbar;
 };
 
 PORTABLE_INLINE_FUNCTION void MGUsup::CheckParams() const {
@@ -178,6 +183,7 @@ PORTABLE_INLINE_FUNCTION void MGUsup::CheckParams() const {
   if (_Cv0 < 0.0) {
     PORTABLE_ALWAYS_THROW_OR_ABORT("Required MGUsup model parameter Cv0 < 0");
   }
+  _AZbar.CheckParams();
 }
 
 PORTABLE_INLINE_FUNCTION Real MGUsup::HugPressureFromDensity(Real rho) const {
diff --git a/singularity-eos/eos/eos_noble_abel.hpp b/singularity-eos/eos/eos_noble_abel.hpp
index a39816996b..e937d153b9 100644
--- a/singularity-eos/eos/eos_noble_abel.hpp
+++ b/singularity-eos/eos/eos_noble_abel.hpp
@@ -36,24 +36,27 @@ using namespace eos_base;
 class NobleAbel : public EosBase<NobleAbel> {
  public:
   NobleAbel() = default;
-  PORTABLE_INLINE_FUNCTION NobleAbel(Real gm1, Real Cv, Real bb, Real qq)
+  PORTABLE_INLINE_FUNCTION
+  NobleAbel(Real gm1, Real Cv, Real bb, Real qq,
+            const MeanAtomicProperties &AZbar = MeanAtomicProperties())
       : _Cv(Cv), _gm1(gm1), _bb(bb), _qq(qq), _T0(ROOM_TEMPERATURE),
         _P0(ATMOSPHERIC_PRESSURE), _qp(0.0),
         _rho0(robust::ratio(_P0, gm1 * Cv * _T0 + bb * _P0)),
         _vol0(robust::ratio(gm1 * Cv * _T0 + bb * _P0, _P0)), _sie0(Cv * _T0 + qq),
         _bmod0(robust::ratio(_rho0 * Cv * _T0 * gm1 * (gm1 + 1.0),
                              (1.0 - bb * _rho0) * (1.0 - bb * _rho0))),
-        _dpde0(robust::ratio(_rho0 * gm1, 1.0 - bb * _rho0)) {
+        _dpde0(robust::ratio(_rho0 * gm1, 1.0 - bb * _rho0)), _AZbar(AZbar) {
     CheckParams();
   }
-  PORTABLE_INLINE_FUNCTION NobleAbel(Real gm1, Real Cv, Real bb, Real qq, Real qp,
-                                     Real T0, Real P0)
+  PORTABLE_INLINE_FUNCTION
+  NobleAbel(Real gm1, Real Cv, Real bb, Real qq, Real qp, Real T0, Real P0,
+            const MeanAtomicProperties &AZbar = MeanAtomicProperties())
       : _Cv(Cv), _gm1(gm1), _bb(bb), _qq(qq), _T0(T0), _P0(P0), _qp(qp),
         _rho0(robust::ratio(P0, gm1 * Cv * T0 + bb * P0)),
         _vol0(robust::ratio(gm1 * Cv * T0 + bb * P0, P0)), _sie0(Cv * T0 + qq),
         _bmod0(robust::ratio(_rho0 * Cv * T0 * gm1 * (gm1 + 1.0),
                              (1.0 - bb * _rho0) * (1.0 - bb * _rho0))),
-        _dpde0(robust::ratio(_rho0 * gm1, 1.0 - bb * _rho0)) {
+        _dpde0(robust::ratio(_rho0 * gm1, 1.0 - bb * _rho0)), _AZbar(AZbar) {
     CheckParams();
   }
   NobleAbel GetOnDevice() { return *this; }
@@ -67,6 +70,7 @@ class NobleAbel : public EosBase<NobleAbel> {
     PORTABLE_ALWAYS_REQUIRE(_Cv > 0, "Heat capacity must be positive");
     PORTABLE_ALWAYS_REQUIRE(_gm1 >= 0, "Gruneisen parameter must be positive");
     PORTABLE_ALWAYS_REQUIRE(_bb >= 0, "Covolume must be positive");
+    _AZbar.CheckParams();
   }
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION Real InternalEnergyFromDensityTemperature(
@@ -175,6 +179,7 @@ class NobleAbel : public EosBase<NobleAbel> {
   }
   // Generic functions provided by the base class. These contain e.g. the vector
   // overloads that use the scalar versions declared here
+  SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(_AZbar)
   SG_ADD_BASE_CLASS_USINGS(NobleAbel)
   PORTABLE_INLINE_FUNCTION
   int nlambda() const noexcept { return 0; }
@@ -187,6 +192,7 @@ class NobleAbel : public EosBase<NobleAbel> {
     printf("Noble-Abel Parameters:\nGamma = %g\nCv    = %g\nb     = %g\nq     = "
            "%g\n",
            _gm1 + 1.0, _Cv, _bb, _qq);
+    _AZbar.PrintParams();
   }
   template <typename Indexer_t>
   PORTABLE_INLINE_FUNCTION void
@@ -206,6 +212,7 @@ class NobleAbel : public EosBase<NobleAbel> {
   Real _T0, _P0, _qp;
   // reference values
   Real _rho0, _vol0, _sie0, _bmod0, _dpde0;
+  MeanAtomicProperties _AZbar;
   // static constexpr const Real _T0 = ROOM_TEMPERATURE;
   // static constexpr const Real _P0 = ATMOSPHERIC_PRESSURE;
   static constexpr const unsigned long _preferred_input =
diff --git a/singularity-eos/eos/eos_powermg.hpp b/singularity-eos/eos/eos_powermg.hpp
index ea8904285a..7e400b89a7 100644
--- a/singularity-eos/eos/eos_powermg.hpp
+++ b/singularity-eos/eos/eos_powermg.hpp
@@ -38,8 +38,10 @@ class PowerMG : public EosBase<PowerMG> {
   PowerMG() = default;
   // Constructor
   PowerMG(const Real rho0, const Real T0, const Real G0, const Real Cv0, const Real E0,
-          const Real S0, const Real Pmin, const Real *expconsts)
-      : _rho0(rho0), _T0(T0), _G0(G0), _Cv0(Cv0), _E0(E0), _S0(S0), _Pmin(Pmin) {
+          const Real S0, const Real Pmin, const Real *expconsts,
+          const MeanAtomicProperties &AZbar = MeanAtomicProperties())
+      : _rho0(rho0), _T0(T0), _G0(G0), _Cv0(Cv0), _E0(E0), _S0(S0), _Pmin(Pmin),
+        _AZbar(AZbar) {
     _InitializePowerMG(expconsts);
     CheckParams();
     if (_Pmin >= 0.0) {
@@ -136,6 +138,7 @@ class PowerMG : public EosBase<PowerMG> {
                          Indexer_t &&lambda = static_cast<Real *>(nullptr)) const;
   // Generic functions provided by the base class. These contain e.g. the vector
   // overloads that use the scalar versions declared here
+  SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(_AZbar)
   SG_ADD_BASE_CLASS_USINGS(PowerMG)
   PORTABLE_INLINE_FUNCTION
   int nlambda() const noexcept { return 0; }
@@ -167,6 +170,7 @@ class PowerMG : public EosBase<PowerMG> {
   static constexpr const unsigned long _preferred_input =
       thermalqs::density | thermalqs::specific_internal_energy;
   Real _rho0, _T0, _G0, _Cv0, _E0, _S0, _Pmin;
+  MeanAtomicProperties _AZbar;
   static constexpr const int PressureCoeffsK0toK40Size = 41;
   int _M;
   Real _K0toK40[PressureCoeffsK0toK40Size];
@@ -194,6 +198,7 @@ PORTABLE_INLINE_FUNCTION void PowerMG::CheckParams() const {
   if (_K0toK40[0] < 0.0) {
     PORTABLE_ALWAYS_THROW_OR_ABORT("Required PowerMG model parameter K0 < 0");
   }
+  _AZbar.CheckParams();
 }
 
 inline void PowerMG::_InitializePowerMG(const Real *K0toK40input) {
diff --git a/singularity-eos/eos/eos_sap_polynomial.hpp b/singularity-eos/eos/eos_sap_polynomial.hpp
index c26e77d429..7d361190c9 100644
--- a/singularity-eos/eos/eos_sap_polynomial.hpp
+++ b/singularity-eos/eos/eos_sap_polynomial.hpp
@@ -39,15 +39,17 @@ class SAP_Polynomial : public EosBase<SAP_Polynomial> {
   PORTABLE_INLINE_FUNCTION
   SAP_Polynomial(const Real rho0, const Real a0, const Real a1, const Real a2c,
                  const Real a2e, const Real a3, const Real b0, const Real b1,
-                 const Real b2c, const Real b2e, const Real b3)
+                 const Real b2c, const Real b2e, const Real b3,
+                 const MeanAtomicProperties &AZbar = MeanAtomicProperties())
       : _a0(a0), _a1(a1), _a2c(a2c), _a2e(a2e), _a3(a3), _b0(b0), _b1(b1), _b2c(b2c),
-        _b2e(b2e), _b3(b3), _rho0(rho0) {
+        _b2e(b2e), _b3(b3), _rho0(rho0), _AZbar(AZbar) {
     CheckParams();
   }
 
   SAP_Polynomial GetOnDevice() { return *this; }
   PORTABLE_INLINE_FUNCTION void CheckParams() const {
     PORTABLE_ALWAYS_REQUIRE(_rho0 >= 0, "Reference density must be non-negative");
+    _AZbar.CheckParams();
   }
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION Real TemperatureFromDensityInternalEnergy(
@@ -183,6 +185,7 @@ class SAP_Polynomial : public EosBase<SAP_Polynomial> {
   }
   // Generic functions provided by the base class. These contain e.g. the vector
   // overloads that use the scalar versions declared here
+  SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(_AZbar)
   SG_ADD_BASE_CLASS_USINGS(SAP_Polynomial)
   PORTABLE_INLINE_FUNCTION
   int nlambda() const noexcept { return 0; }
@@ -204,6 +207,7 @@ class SAP_Polynomial : public EosBase<SAP_Polynomial> {
     printf("      b2c = %g\n", _b2c);
     printf("      b2e = %g\n", _b2e);
     printf("      b3  = %g\n", _b3);
+    _AZbar.PrintParams();
   }
   template <typename Indexer_t>
   PORTABLE_INLINE_FUNCTION void
@@ -221,7 +225,7 @@ class SAP_Polynomial : public EosBase<SAP_Polynomial> {
   Real _a0, _a1, _a2c, _a2e, _a3, _b0, _b1, _b2c, _b2e, _b3;
   // reference values
   Real _rho0;
-
+  MeanAtomicProperties _AZbar;
   static constexpr const unsigned long _preferred_input =
       thermalqs::density | thermalqs::specific_internal_energy;
 };
diff --git a/singularity-eos/eos/eos_spiner.hpp b/singularity-eos/eos/eos_spiner.hpp
index 6da8e47c4a..7a9bb2e285 100644
--- a/singularity-eos/eos/eos_spiner.hpp
+++ b/singularity-eos/eos/eos_spiner.hpp
@@ -80,6 +80,7 @@ class SpinerEOSDependsRhoT : public EosBase<SpinerEOSDependsRhoT> {
   struct Lambda {
     enum Index { lRho = 0, lT = 1 };
   };
+  SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(AZbar_)
   SG_ADD_BASE_CLASS_USINGS(SpinerEOSDependsRhoT);
   inline SpinerEOSDependsRhoT(const std::string &filename, int matid,
                               bool reproduciblity_mode = false);
@@ -307,6 +308,7 @@ class SpinerEOSDependsRhoT : public EosBase<SpinerEOSDependsRhoT> {
   Real rhoNormal_, TNormal_, sieNormal_, PNormal_;
   Real CvNormal_, bModNormal_, dPdENormal_, dVdTNormal_;
   Real lRhoOffset_, lTOffset_; // offsets must be non-negative
+  MeanAtomicProperties AZbar_;
   int matid_;
   bool reproducible_;
   // whereAmI_ and status_ used only for reporting. They are not thread-safe.
@@ -349,6 +351,7 @@ class SpinerEOSDependsRhoSie : public EosBase<SpinerEOSDependsRhoSie> {
   };
   using STricks = table_utils::SpinerTricks<SpinerEOSDependsRhoSie>;
 
+  SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(AZbar_)
   SG_ADD_BASE_CLASS_USINGS(SpinerEOSDependsRhoSie);
   PORTABLE_INLINE_FUNCTION SpinerEOSDependsRhoSie()
       : memoryStatus_(DataStatus::Deallocated) {}
@@ -538,6 +541,7 @@ class SpinerEOSDependsRhoSie : public EosBase<SpinerEOSDependsRhoSie> {
       thermalqs::density | thermalqs::temperature;
   // static constexpr const char _eos_type[] = "SpinerEOSDependsRhoSie";
   int matid_;
+  MeanAtomicProperties AZbar_;
   bool reproducible_;
   mutable RootFinding1D::Status status_;
   static constexpr const int _n_lambda = 1;
@@ -716,6 +720,11 @@ inline herr_t SpinerEOSDependsRhoT::loadDataboxes_(const std::string &matid_str,
   status += H5LTget_attribute_double(file, matid_str.c_str(),
                                      SP5::Material::normalDensity, &rhoNormal_);
   rhoNormal_ = std::abs(rhoNormal_);
+  // Mean atomic mass and mean atomic number
+  status += H5LTget_attribute_double(file, matid_str.c_str(),
+                                     SP5::Material::meanAtomicMass, &(AZbar_.Abar));
+  status += H5LTget_attribute_double(file, matid_str.c_str(),
+                                     SP5::Material::meanAtomicNumber, &(AZbar_.Zbar));
 
   // tables
   status += P_.loadHDF(lTGroup, SP5::Fields::P);
@@ -1552,6 +1561,11 @@ herr_t SpinerEOSDependsRhoSie::loadDataboxes_(const std::string &matid_str, hid_
   status += H5LTget_attribute_double(file, matid_str.c_str(),
                                      SP5::Material::normalDensity, &rhoNormal_);
   rhoNormal_ = std::abs(rhoNormal_);
+  // Mean atomic mass and mean atomic number
+  status += H5LTget_attribute_double(file, matid_str.c_str(),
+                                     SP5::Material::meanAtomicMass, &(AZbar_.Abar));
+  status += H5LTget_attribute_double(file, matid_str.c_str(),
+                                     SP5::Material::meanAtomicNumber, &(AZbar_.Zbar));
 
   // sometimes independent variables
   status += sie_.loadHDF(lTGroup, SP5::Fields::sie);
diff --git a/singularity-eos/eos/eos_stellar_collapse.hpp b/singularity-eos/eos/eos_stellar_collapse.hpp
index 9901fffa8c..4ff7788dda 100644
--- a/singularity-eos/eos/eos_stellar_collapse.hpp
+++ b/singularity-eos/eos/eos_stellar_collapse.hpp
@@ -162,6 +162,34 @@ class StellarCollapse : public EosBase<StellarCollapse> {
                                        Indexer_t &&lambda, Real &rho, Real &sie) const;
 
   // Properties of an NSE EOS
+  PORTABLE_INLINE_FUNCTION
+  Real MeanAtomicMass() const {
+    PORTABLE_THROW_OR_ABORT(
+        "For Stellar Collapse EOS, mean atomic mass is a state variable!\n");
+    return 1.0;
+  }
+  PORTABLE_INLINE_FUNCTION
+  Real MeanAtomicNumber() const {
+    PORTABLE_THROW_OR_ABORT(
+        "For Stellar Collapse EOS, mean atomic number is a state variable!\n");
+    return 1.0;
+  }
+  template <typename Indexer_t = Real *>
+  PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature(
+      const Real rho, const Real T,
+      Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
+    Real lRho, lT, Ye;
+    getLogsFromRhoT_(rho, T, lambda, lRho, lT, Ye);
+    return Abar_.interpToReal(Ye, lT, lRho);
+  }
+  template <typename Indexer_t = Real *>
+  PORTABLE_INLINE_FUNCTION Real MeanAtomicNumberFromDensityTemperature(
+      const Real rho, const Real T,
+      Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
+    Real lRho, lT, Ye;
+    getLogsFromRhoT_(rho, T, lambda, lRho, lT, Ye);
+    return Zbar_.interpToReal(Ye, lT, lRho);
+  }
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION void MassFractionsFromDensityTemperature(
       const Real rho, const Real temperature, Real &Xa, Real &Xh, Real &Xn, Real &Xp,
diff --git a/singularity-eos/eos/eos_stiff.hpp b/singularity-eos/eos/eos_stiff.hpp
index 4a64aeffc5..4c33ca77bc 100644
--- a/singularity-eos/eos/eos_stiff.hpp
+++ b/singularity-eos/eos/eos_stiff.hpp
@@ -36,24 +36,27 @@ using namespace eos_base;
 class StiffGas : public EosBase<StiffGas> {
  public:
   StiffGas() = default;
-  PORTABLE_INLINE_FUNCTION StiffGas(Real gm1, Real Cv, Real Pinf, Real qq)
+  PORTABLE_INLINE_FUNCTION
+  StiffGas(Real gm1, Real Cv, Real Pinf, Real qq,
+           const MeanAtomicProperties &AZbar = MeanAtomicProperties())
       : _Cv(Cv), _gm1(gm1), _Pinf(Pinf), _qq(qq), _T0(ROOM_TEMPERATURE),
         _P0(ATMOSPHERIC_PRESSURE), _qp(0.0),
         _rho0(robust::ratio(_P0 + Pinf, gm1 * Cv * _T0)),
         _vol0(robust::ratio(gm1 * Cv * _T0, _P0 + Pinf)),
         _sie0(robust::ratio(_P0 + (gm1 + 1.0) * Pinf, _P0 + Pinf) * Cv * _T0 + qq),
         _bmod0(gm1 * (gm1 + 1.0) * robust::ratio(_P0 + Pinf, gm1 * Cv * _T0) * Cv * _T0),
-        _dpde0(_rho0 * _gm1) {
+        _dpde0(_rho0 * _gm1), _AZbar(AZbar) {
     CheckParams();
   }
-  PORTABLE_INLINE_FUNCTION StiffGas(Real gm1, Real Cv, Real Pinf, Real qq, Real qp,
-                                    Real T0, Real P0)
+  PORTABLE_INLINE_FUNCTION
+  StiffGas(Real gm1, Real Cv, Real Pinf, Real qq, Real qp, Real T0, Real P0,
+           const MeanAtomicProperties &AZbar = MeanAtomicProperties())
       : _Cv(Cv), _gm1(gm1), _Pinf(Pinf), _qq(qq), _T0(T0), _P0(P0), _qp(qp),
         _rho0(robust::ratio(P0 + Pinf, gm1 * Cv * T0)),
         _vol0(robust::ratio(gm1 * Cv * _T0, _P0 + Pinf)),
         _sie0(robust::ratio(P0 + (gm1 + 1.0) * Pinf, P0 + Pinf) * Cv * T0 + qq),
         _bmod0(gm1 * (gm1 + 1.0) * robust::ratio(P0 + Pinf, gm1 * Cv * T0) * Cv * T0),
-        _dpde0(_rho0 * _gm1) {
+        _dpde0(_rho0 * _gm1), _AZbar(AZbar) {
     CheckParams();
   }
   StiffGas GetOnDevice() { return *this; }
@@ -66,6 +69,7 @@ class StiffGas : public EosBase<StiffGas> {
   PORTABLE_INLINE_FUNCTION void CheckParams() const {
     PORTABLE_ALWAYS_REQUIRE(_Cv >= 0, "Heat capacity must be positive");
     PORTABLE_ALWAYS_REQUIRE(_gm1 >= 0, "Gruneisen parameter must be positive");
+    _AZbar.CheckParams();
   }
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION Real InternalEnergyFromDensityTemperature(
@@ -168,6 +172,7 @@ class StiffGas : public EosBase<StiffGas> {
   // Generic functions provided by the base class. These contain e.g. the vector
   // overloads that use the scalar versions declared here
   SG_ADD_BASE_CLASS_USINGS(StiffGas)
+  SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(_AZbar)
   PORTABLE_INLINE_FUNCTION
   int nlambda() const noexcept { return 0; }
   static constexpr unsigned long PreferredInput() { return _preferred_input; }
@@ -179,6 +184,7 @@ class StiffGas : public EosBase<StiffGas> {
     printf("Stiff Gas Parameters:\nGamma = %g\nCv    = %g\nPinf  = %g\nq     = "
            "%g\n",
            _gm1 + 1.0, _Cv, _Pinf, _qq);
+    _AZbar.PrintParams();
   }
   template <typename Indexer_t>
   PORTABLE_INLINE_FUNCTION void
@@ -199,6 +205,7 @@ class StiffGas : public EosBase<StiffGas> {
   Real _T0, _P0, _qp;
   // reference values
   Real _rho0, _vol0, _sie0, _bmod0, _dpde0;
+  MeanAtomicProperties _AZbar;
   // static constexpr const Real _T0 = ROOM_TEMPERATURE;
   // static constexpr const Real _P0 = ATMOSPHERIC_PRESSURE;
   static constexpr const unsigned long _preferred_input =
diff --git a/singularity-eos/eos/eos_variant.hpp b/singularity-eos/eos/eos_variant.hpp
index 36711df58f..31dcb7867a 100644
--- a/singularity-eos/eos/eos_variant.hpp
+++ b/singularity-eos/eos/eos_variant.hpp
@@ -333,6 +333,36 @@ class Variant {
     return mpark::visit([](const auto &eos) { return eos.MinimumTemperature(); }, eos_);
   }
 
+  // Atomic mass/atomic number functions
+  PORTABLE_INLINE_FUNCTION
+  Real MeanAtomicMass() const {
+    return mpark::visit([](const auto &eos) { return eos.MeanAtomicMass(); }, eos_);
+  }
+  PORTABLE_INLINE_FUNCTION
+  Real MeanAtomicNumber() const {
+    return mpark::visit([](const auto &eos) { return eos.MeanAtomicNumber(); }, eos_);
+  }
+  template <typename Indexer_t = Real *>
+  PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature(
+      const Real rho, const Real T,
+      Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
+    return mpark::visit(
+        [&rho, &T, &lambda](const auto &eos) {
+          return eos.MeanAtomicMassFromDensityTemperature(rho, T, lambda);
+        },
+        eos_);
+  }
+  template <typename Indexer_t = Real *>
+  PORTABLE_INLINE_FUNCTION Real MeanAtomicNumberFromDensityTemperature(
+      const Real rho, const Real T,
+      Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
+    return mpark::visit(
+        [&rho, &T, &lambda](const auto &eos) {
+          return eos.MeanAtomicNumberFromDensityTemperature(rho, T, lambda);
+        },
+        eos_);
+  }
+
   /*
   Vector versions of the member functions run on the host but the scalar
   lookups will run on the device
diff --git a/singularity-eos/eos/eos_vinet.hpp b/singularity-eos/eos/eos_vinet.hpp
index 570d1bb319..8ea913e404 100644
--- a/singularity-eos/eos/eos_vinet.hpp
+++ b/singularity-eos/eos/eos_vinet.hpp
@@ -38,10 +38,13 @@ using namespace robust;
 class Vinet : public EosBase<Vinet> {
  public:
   Vinet() = default;
-  // Constructor
+  // TODO(JMM): Should Vinet and MGPower take AZbar before
+  // the series coefficients? Would require an overload...
   Vinet(const Real rho0, const Real T0, const Real B0, const Real BP0, const Real A0,
-        const Real Cv0, const Real E0, const Real S0, const Real *expconsts)
-      : _rho0(rho0), _T0(T0), _B0(B0), _BP0(BP0), _A0(A0), _Cv0(Cv0), _E0(E0), _S0(S0) {
+        const Real Cv0, const Real E0, const Real S0, const Real *expconsts,
+        const MeanAtomicProperties &AZbar = MeanAtomicProperties())
+      : _rho0(rho0), _T0(T0), _B0(B0), _BP0(BP0), _A0(A0), _Cv0(Cv0), _E0(E0), _S0(S0),
+        _AZbar(AZbar) {
     CheckParams();
     InitializeVinet(expconsts);
   }
@@ -131,6 +134,7 @@ class Vinet : public EosBase<Vinet> {
   // Generic functions provided by the base class. These contain e.g. the vector
   // overloads that use the scalar versions declared here
   SG_ADD_BASE_CLASS_USINGS(Vinet)
+  SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(_AZbar)
   PORTABLE_INLINE_FUNCTION
   int nlambda() const noexcept { return 0; }
   static constexpr unsigned long PreferredInput() { return _preferred_input; }
@@ -161,6 +165,7 @@ class Vinet : public EosBase<Vinet> {
   static constexpr const unsigned long _preferred_input =
       thermalqs::density | thermalqs::temperature;
   Real _rho0, _T0, _B0, _BP0, _A0, _Cv0, _E0, _S0;
+  MeanAtomicProperties _AZbar;
   static constexpr const int PressureCoeffsd2tod40Size = 39;
   static constexpr const int VinetInternalParametersSize = PressureCoeffsd2tod40Size + 4;
   Real _VIP[VinetInternalParametersSize], _d2tod40[PressureCoeffsd2tod40Size];
@@ -189,6 +194,7 @@ PORTABLE_INLINE_FUNCTION void Vinet::CheckParams() const {
   if (_Cv0 < 0.0) {
     PORTABLE_ALWAYS_THROW_OR_ABORT("Required Vinet model parameter Cv0 < 0");
   }
+  _AZbar.CheckParams();
 }
 
 inline void Vinet::InitializeVinet(const Real *d2tod40input) {
diff --git a/singularity-eos/eos/modifiers/eos_unitsystem.hpp b/singularity-eos/eos/modifiers/eos_unitsystem.hpp
index 53d9435aab..957aa5eed6 100644
--- a/singularity-eos/eos/modifiers/eos_unitsystem.hpp
+++ b/singularity-eos/eos/modifiers/eos_unitsystem.hpp
@@ -249,6 +249,21 @@ class UnitSystem : public EosBase<UnitSystem<T>> {
     return inv_temp_unit_ * t_.MinimumTemperature();
   }
 
+  template <typename Indexer_t = Real *>
+  PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature(
+      const Real rho, const Real temperature,
+      Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
+    return t_.MeanAtomicMassFromDensityTemperature(rho * rho_unit_,
+                                                   temperature * temp_unit_, lambda);
+  }
+  template <typename Indexer_t = Real *>
+  PORTABLE_INLINE_FUNCTION Real MeanAtomicNumberFromDensityTemperature(
+      const Real rho, const Real temperature,
+      Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
+    return t_.MeanAtomicNumberFromDensityTemperature(rho * rho_unit_,
+                                                     temperature * temp_unit_, lambda);
+  }
+
   // vector implementations
   template <typename LambdaIndexer>
   inline void TemperatureFromDensityInternalEnergy(
diff --git a/singularity-eos/eos/modifiers/ramps_eos.hpp b/singularity-eos/eos/modifiers/ramps_eos.hpp
index a0ce08db2e..f316c57dec 100644
--- a/singularity-eos/eos/modifiers/ramps_eos.hpp
+++ b/singularity-eos/eos/modifiers/ramps_eos.hpp
@@ -210,6 +210,7 @@ class BilinearRampEOS : public EosBase<BilinearRampEOS<T>> {
   template <typename Indexer_t = Real *>
   PORTABLE_FUNCTION Real GruneisenParamFromDensityTemperature(
       const Real rho, const Real temperature, Indexer_t &&lambda = nullptr) const {
+    // TODO(JMM): This doesn't seem right. dpdrho relevant.
     return t_.GruneisenParamFromDensityTemperature(rho, temperature, lambda);
   }
   template <typename Indexer_t = Real *>
@@ -242,6 +243,19 @@ class BilinearRampEOS : public EosBase<BilinearRampEOS<T>> {
     return;
   }
 
+  template <typename Indexer_t = Real *>
+  PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature(
+      const Real rho, const Real temperature,
+      Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
+    return t_.MeanAtomicMassFromDensityTemperature(rho, temperature, lambda);
+  }
+  template <typename Indexer_t = Real *>
+  PORTABLE_INLINE_FUNCTION Real MeanAtomicNumberFromDensityTemperature(
+      const Real rho, const Real temperature,
+      Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
+    return t_.MeanAtomicNumberFromDensityTemperature(rho, temperature, lambda);
+  }
+
   // vector implementations
   template <typename LambdaIndexer>
   inline void TemperatureFromDensityInternalEnergy(
diff --git a/singularity-eos/eos/modifiers/relativistic_eos.hpp b/singularity-eos/eos/modifiers/relativistic_eos.hpp
index 0fd296d78d..8729c30a83 100644
--- a/singularity-eos/eos/modifiers/relativistic_eos.hpp
+++ b/singularity-eos/eos/modifiers/relativistic_eos.hpp
@@ -182,6 +182,19 @@ class RelativisticEOS : public EosBase<RelativisticEOS<T>> {
     t_.ValuesAtReferenceState(rho, temp, sie, press, cv, bmod, dpde, dvdt, lambda);
   }
 
+  template <typename Indexer_t = Real *>
+  PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature(
+      const Real rho, const Real temperature,
+      Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
+    return t_.MeanAtomicMassFromDensityTemperature(rho, temperature, lambda);
+  }
+  template <typename Indexer_t = Real *>
+  PORTABLE_INLINE_FUNCTION Real MeanAtomicNumberFromDensityTemperature(
+      const Real rho, const Real temperature,
+      Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
+    return t_.MeanAtomicNumberFromDensityTemperature(rho, temperature, lambda);
+  }
+
   SG_ADD_MODIFIER_METHODS(T, t_);
 
  private:
diff --git a/singularity-eos/eos/modifiers/scaled_eos.hpp b/singularity-eos/eos/modifiers/scaled_eos.hpp
index bf64e90a04..9576c306c4 100644
--- a/singularity-eos/eos/modifiers/scaled_eos.hpp
+++ b/singularity-eos/eos/modifiers/scaled_eos.hpp
@@ -342,6 +342,19 @@ class ScaledEOS : public EosBase<ScaledEOS<T>> {
     return t_.MinimumTemperature();
   }
 
+  template <typename Indexer_t = Real *>
+  PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature(
+      const Real rho, const Real temperature,
+      Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
+    return t_.MeanAtomicMassFromDensityTemperature(scale_ * rho, temperature, lambda);
+  }
+  template <typename Indexer_t = Real *>
+  PORTABLE_INLINE_FUNCTION Real MeanAtomicNumberFromDensityTemperature(
+      const Real rho, const Real temperature,
+      Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
+    return t_.MeanAtomicNumberFromDensityTemperature(scale_ * rho, temperature, lambda);
+  }
+
   SG_ADD_MODIFIER_METHODS(T, t_);
 
  private:
diff --git a/singularity-eos/eos/modifiers/shifted_eos.hpp b/singularity-eos/eos/modifiers/shifted_eos.hpp
index 79eaad6e4e..c943bec517 100644
--- a/singularity-eos/eos/modifiers/shifted_eos.hpp
+++ b/singularity-eos/eos/modifiers/shifted_eos.hpp
@@ -352,6 +352,19 @@ class ShiftedEOS : public EosBase<ShiftedEOS<T>> {
     return t_.MinimumTemperature();
   }
 
+  template <typename Indexer_t = Real *>
+  PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature(
+      const Real rho, const Real temperature,
+      Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
+    return t_.MeanAtomicMassFromDensityTemperature(rho, temperature, lambda);
+  }
+  template <typename Indexer_t = Real *>
+  PORTABLE_INLINE_FUNCTION Real MeanAtomicNumberFromDensityTemperature(
+      const Real rho, const Real temperature,
+      Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
+    return t_.MeanAtomicNumberFromDensityTemperature(rho, temperature, lambda);
+  }
+
   SG_ADD_MODIFIER_METHODS(T, t_);
 
  private:
diff --git a/test/test_eos_ideal.cpp b/test/test_eos_ideal.cpp
index 11b8b6447a..dafc287b73 100644
--- a/test/test_eos_ideal.cpp
+++ b/test/test_eos_ideal.cpp
@@ -30,6 +30,7 @@
 #include <test/eos_unit_test_helpers.hpp>
 
 using singularity::IdealGas;
+using singularity::MeanAtomicProperties;
 using EOS = singularity::Variant<IdealGas>;
 
 SCENARIO("Ideal gas entropy", "[IdealGas][Entropy][GibbsFreeEnergy]") {
@@ -80,6 +81,44 @@ SCENARIO("Ideal gas entropy", "[IdealGas][Entropy][GibbsFreeEnergy]") {
   }
 }
 
+SCENARIO("Ideal gas mean atomic properties",
+         "[IdealGas][MeanAtomicMass][MeanAtomicNumber]") {
+  constexpr Real Cv = 5.0;
+  constexpr Real gm1 = 0.4;
+  constexpr Real Abar = 4.0; // Helium
+  constexpr Real Zbar = 2.0;
+  const MeanAtomicProperties azbar(Abar, Zbar);
+  GIVEN("An ideal gas initialized with mean atomic poroperties") {
+    EOS host_eos = IdealGas(gm1, Cv, azbar);
+    WHEN("We evaluate it on host") {
+      Real Ab_eval = host_eos.MeanAtomicMass();
+      Real Zb_eval = host_eos.MeanAtomicNumber();
+      THEN("We get the right answer") {
+        REQUIRE(isClose(Ab_eval, Abar, 1e-12));
+        REQUIRE(isClose(Zb_eval, Zbar, 1e-12));
+      }
+    }
+    WHEN("We evaluate it on device, using a loop") {
+      constexpr int N = 100;
+      auto device_eos = host_eos.GetOnDevice();
+      int nwrong = 0;
+      portableReduce(
+          "Check mean atomic number", 0, N,
+          PORTABLE_LAMBDA(const int i, int &nw) {
+            double rho = i;
+            double T = 100.0 * i;
+            double Ab_eval = device_eos.MeanAtomicMassFromDensityTemperature(rho, T);
+            double Zb_eval = device_eos.MeanAtomicNumberFromDensityTemperature(rho, T);
+            nw += !(isClose(Ab_eval, Abar, 1e-12)) + !(isClose(Zb_eval, Zbar, 1e-12));
+          },
+          nwrong);
+      REQUIRE(nwrong == 0);
+      device_eos.Finalize();
+    }
+    host_eos.Finalize();
+  }
+}
+
 // A non-standard pattern where we do a reduction
 class CheckPofRE {
  public:
diff --git a/test/test_eos_tabulated.cpp b/test/test_eos_tabulated.cpp
index 9d0dcb72d4..6a31ed9451 100644
--- a/test/test_eos_tabulated.cpp
+++ b/test/test_eos_tabulated.cpp
@@ -86,6 +86,11 @@ SCENARIO("SpinerEOS depends on Rho and T", "[SpinerEOS][DependsRhoT][EOSPAC]") {
         REQUIRE(isClose(T, T_pac));
       }
 
+      AND_THEN("We get the correct mean atomic mass and number") {
+        REQUIRE(steelEOS_host.MeanAtomicMass() == eospac.MeanAtomicMass());
+        REQUIRE(steelEOS_host.MeanAtomicNumber() == eospac.MeanAtomicNumber());
+      }
+
       // TODO: this needs to be a much more rigorous test
       AND_THEN("Quantities can be read from density and temperature") {
         const Real sie_pac = eospac.InternalEnergyFromDensityTemperature(1e0, 1e6);
@@ -233,6 +238,11 @@ SCENARIO("SpinerEOS depends on rho and sie", "[SpinerEOS][DependsRhoSie]") {
     THEN("The correct metadata is read in") {
       REQUIRE(steelEOS_host.matid() == steelID);
 
+      AND_THEN("We get the correct mean atomic mass and number") {
+        REQUIRE(isClose(steelEOS_host.MeanAtomicMass(), 55.37));
+        REQUIRE(isClose(steelEOS_host.MeanAtomicNumber(), 25.80));
+      }
+
       int nw_ie2{0}, nw_te2{0};
 #ifdef PORTABILITY_STRATEGY_KOKKOS
       Kokkos::fence();