diff --git a/CHANGELOG.md b/CHANGELOG.md
index 16b4f4c7e52..86c08afc012 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -3,6 +3,7 @@
 ## Current develop
 
 ### Added (new features/APIs/variables/...)
+- [[PR618]](https://github.com/lanl/singularity-eos/pull/618) Add PTDerivativesFromPreferred for computing derivatives of a mixture in a cell
 
 ### Fixed (Repair bugs, etc)
 
diff --git a/doc/sphinx/src/using-eos.rst b/doc/sphinx/src/using-eos.rst
index faae9c7fbe7..b196e526f9a 100644
--- a/doc/sphinx/src/using-eos.rst
+++ b/doc/sphinx/src/using-eos.rst
@@ -1334,6 +1334,31 @@ specific internal energy in :math:`erg/g`.
 
 The function
 
+.. code-block:: cpp
+
+  template <typename Lambda_t = Real *>
+  PORTABLE_INLINE_FUNCTION void
+  PTDerivativesFromPreferred(const Real rho, const Real sie, const Real P, const Real T,
+                             Lambda_t &&lambda, Real &dedP_T, Real &drdP_T, Real &dedT_P,
+                             Real &drdT_P) const;
+
+computes the partial derivatives of density and specific internal
+energy with respect to pressure and temperature, with either pressure
+or temperature fixed. Each EOS model expects consistent density, 
+energy, pressure and temperature values to be provided so that it
+can perform this calculation performantly. The intended
+use of this method is to compute the cell-averaged thermodynamic
+derivatives in a mixed cell in pressure-temperature equilibrium.
+
+.. warning::
+
+  Not all equations of state provide implementations of
+  ``PTDerivativesFromPreferred``, when an implementation is not
+  available, ``singularity-eos`` uses a finite differences
+  approximation.
+
+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 46840465337..ce94ce24e55 100644
--- a/singularity-eos/eos/eos_base.hpp
+++ b/singularity-eos/eos/eos_base.hpp
@@ -206,6 +206,33 @@ struct MeanAtomicProperties {
   }
 };
 
+/* A version of PT Derivatives that uses only finite
+   differences. Factored out of the base class so it can be used to
+   cross-validate per-class implementations.
+*/
+template <typename EOS_t, typename Lambda_t>
+PORTABLE_INLINE_FUNCTION void
+PTDerivativesByFiniteDifferences(const EOS_t &eos, const Real rho, const Real sie,
+                                 const Real P, const Real T, const Real derivative_eps,
+                                 Lambda_t &&lambda, Real &dedP_T, Real &drdP_T,
+                                 Real &dedT_P, Real &drdT_P) {
+  Real r_pert, e_pert;
+
+  const Real absP = std::max(std::abs(P), robust::EPS() / derivative_eps);
+  const Real dP = -robust::sgn(P) * absP * derivative_eps; // might be negative
+
+  const Real absT = std::max(std::abs(T), robust::EPS() / derivative_eps);
+  const Real dT = absT * derivative_eps; // can't be negative
+
+  eos.DensityEnergyFromPressureTemperature(P + dP, T, lambda, r_pert, e_pert);
+  drdP_T = robust::ratio(r_pert - rho, dP);
+  dedP_T = robust::ratio(e_pert - sie, dP);
+
+  eos.DensityEnergyFromPressureTemperature(P, T + dT, lambda, r_pert, e_pert);
+  drdT_P = robust::ratio(r_pert - rho, dT);
+  dedT_P = robust::ratio(e_pert - sie, dT);
+}
+
 /*
 This is a CRTP that allows for static inheritance so that default behavior for
 various member functions can be defined.
@@ -288,6 +315,23 @@ class EosBase {
     return sie + (P / rho) - T * S;
   }
 
+  // The base class computes these differences by finite
+  // differences. When available, this should absolutely be overloaded
+  // by an individual EOS, but this provides a default implementation
+  // at the very least.
+  //
+  // TODO(JMM): Should we have a vectorized version of this function?
+  template <typename Lambda_t = Real *>
+  PORTABLE_INLINE_FUNCTION void
+  PTDerivativesFromPreferred(const Real rho, const Real sie, const Real P, const Real T,
+                             Lambda_t &&lambda, Real &dedP_T, Real &drdP_T, Real &dedT_P,
+                             Real &drdT_P) const {
+    constexpr Real derivative_eps = 3.0e-6;
+    const CRTP &copy = *(static_cast<CRTP const *>(this));
+    PTDerivativesByFiniteDifferences(copy, rho, sie, P, T, derivative_eps, lambda, dedP_T,
+                                     drdP_T, dedT_P, drdT_P);
+  }
+
   // Vector member functions
   template <typename RealIndexer, typename ConstRealIndexer, typename LambdaIndexer>
   inline void
diff --git a/singularity-eos/eos/eos_eospac.hpp b/singularity-eos/eos/eos_eospac.hpp
index 82f2f09e880..a2292743c70 100644
--- a/singularity-eos/eos/eos_eospac.hpp
+++ b/singularity-eos/eos/eos_eospac.hpp
@@ -1,5 +1,5 @@
 //------------------------------------------------------------------------------
-// © 2021-2025. Triad National Security, LLC. All rights reserved.  This
+// © 2021-2026. 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.  Department of Energy/National
@@ -231,6 +231,11 @@ class EOSPAC : public EosBase<EOSPAC> {
   PORTABLE_INLINE_FUNCTION Real GruneisenParamFromDensityInternalEnergy(
       const Real rho, const Real sie,
       Indexer_t &&lambda = static_cast<Real *>(nullptr)) const;
+  template <typename Lambda_t>
+  PORTABLE_INLINE_FUNCTION void
+  PTDerivativesFromPreferred(const Real rho, const Real sie, const Real P, const Real T,
+                             Lambda_t &&lambda, Real &dedP_T, Real &drdP_T, Real &dedT_P,
+                             Real &drdT_P) const;
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION void
   FillEos(Real &rho, Real &temp, Real &energy, Real &press, Real &cv, Real &bmod,
@@ -1710,6 +1715,52 @@ PORTABLE_INLINE_FUNCTION void EOSPAC::DensityEnergyFromPressureTemperature(
 #endif // ON DEVICE
 }
 
+template <typename Lambda_t>
+PORTABLE_INLINE_FUNCTION void
+EOSPAC::PTDerivativesFromPreferred(const Real rho, const Real sie, const Real press,
+                                   const Real temp, Lambda_t &&lambda, Real &dedP_T,
+                                   Real &drdP_T, Real &dedT_P, Real &drdT_P) const {
+#if SINGULARITY_ON_DEVICE
+  PORTABLE_ALWAYS_ABORT("EOSPAC calls not supported on device\n");
+#else
+  PORTABLE_REQUIRE(split_ == TableSplit::Total,
+                   "Density of pressure and temperature only supported for total "
+                   "tables at this time");
+  using namespace EospacWrapper;
+  EOS_REAL P[1] = {pressureToSesame(press)};
+  EOS_REAL T[1] = {temperatureToSesame(temp)};
+  EOS_REAL R[1] = {rho};
+  EOS_REAL E[1] = {sieToSesame(sie)};
+  EOS_REAL z[1], dx[1], dy[1];
+  EOS_INTEGER nxypairs = 1;
+  EOS_INTEGER table;
+
+  table = EofRT_table_;
+  eosSafeInterpolate(&table, nxypairs, R, T, z, dx, dy, "EofRT", Verbosity::Quiet);
+  const Real dedr_T = sieFromSesame(dx[0]);
+  const Real dedT_r = sieFromSesame(temperatureToSesame(dy[0]));
+
+  table = PofRT_table_;
+  eosSafeInterpolate(&table, nxypairs, R, T, z, dx, dy, "PofRT", Verbosity::Quiet);
+  const Real dPdr_T = pressureFromSesame(dx[0]);
+  const Real dPdT_r = pressureFromSesame(temperatureToSesame(dy[0]));
+
+  dedP_T = robust::ratio(dedr_T, dPdr_T);
+  dedT_P = dedT_r - robust::ratio(dedr_T * dPdT_r, dPdr_T);
+  drdP_T = robust::ratio(1., dPdr_T);
+  drdT_P = -robust::ratio(dPdT_r, dPdr_T);
+  /*
+  // Alternative approach
+  table = RofPT_table_;
+  eosSafeInterpolate(&table, nxypairs, P, T, z, dx, dy, "RofPT", Verbosity::Quiet);
+  drdP_T = pressureToSesame(dx[0]);
+  drdT_P = temperatureToSesame(dy[0]);
+  dedP_T = dedr_T * drdP_T;
+  dedT_P = dedT_r + dedr_T * drdT_P;
+  */
+#endif // ON DEVICE
+}
+
 template <typename Indexer_t>
 PORTABLE_INLINE_FUNCTION void
 EOSPAC::ValuesAtReferenceState(Real &rho, Real &temp, Real &sie, Real &press, Real &cv,
diff --git a/singularity-eos/eos/eos_ideal.hpp b/singularity-eos/eos/eos_ideal.hpp
index 146a2b66a05..0885ae4658e 100644
--- a/singularity-eos/eos/eos_ideal.hpp
+++ b/singularity-eos/eos/eos_ideal.hpp
@@ -148,6 +148,25 @@ class IdealGas : public EosBase<IdealGas> {
       Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
     return _gm1;
   }
+
+  template <typename Lambda_t = Real *>
+  PORTABLE_INLINE_FUNCTION void
+  PTDerivativesFromPreferred(const Real rho, const Real sie, const Real P, const Real T,
+                             Lambda_t &&lambda, Real &dedP_T, Real &drdP_T, Real &dedT_P,
+                             Real &drdT_P) const {
+    // P = (gm1) rho cv T = (gm1) rho sie
+
+    dedP_T = 0;
+    dedT_P = _Cv;
+
+    // => rho(P, T) = P / (gm1 cv T) = P / (gm1 sie)
+    drdP_T = robust::ratio(1.0, _gm1 * sie);
+
+    // rho(P, T) = (P / (gm1 cv)) T^{-1}
+    // => drdT_P = -1 (P / (gm1 cv)) T^{-2}
+    drdT_P = -robust::ratio(P, _gm1 * _Cv * T * T);
+  }
+
   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_spiner_rho_sie.hpp b/singularity-eos/eos/eos_spiner_rho_sie.hpp
index f0c2986fa4d..617ca8ea732 100644
--- a/singularity-eos/eos/eos_spiner_rho_sie.hpp
+++ b/singularity-eos/eos/eos_spiner_rho_sie.hpp
@@ -215,6 +215,14 @@ class SpinerEOSDependsRhoSieTransformable
   PORTABLE_INLINE_FUNCTION void
   DensityEnergyFromPressureTemperature(const Real press, const Real temp,
                                        Indexer_t &&lambda, Real &rho, Real &sie) const;
+  /*
+  // TODO(JMM): For now using FD. Fix this.
+  template <typename Lambda_t = Real *>
+  PORTABLE_INLINE_FUNCTION void
+  PTDerivativesFromPreferred(const Real rho, const Real sie, const Real P, const Real T,
+                             Lambda_t &&lambda, Real &dedP_T, Real &drdP_T, Real &dedT_P,
+                             Real &drdT_P) const;
+  */
   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_spiner_rho_temp.hpp b/singularity-eos/eos/eos_spiner_rho_temp.hpp
index 1ffb45ab8a4..0f8b49e75a0 100644
--- a/singularity-eos/eos/eos_spiner_rho_temp.hpp
+++ b/singularity-eos/eos/eos_spiner_rho_temp.hpp
@@ -196,6 +196,14 @@ class SpinerEOSDependsRhoT : public EosBase<SpinerEOSDependsRhoT> {
   PORTABLE_INLINE_FUNCTION void
   DensityEnergyFromPressureTemperature(const Real press, const Real temp,
                                        Indexer_t &&lambda, Real &rho, Real &sie) const;
+  /*
+  // TODO(JMM): For now using FD. Fix this.
+  template <typename Lambda_t = Real *>
+  PORTABLE_INLINE_FUNCTION void
+  PTDerivativesFromPreferred(const Real rho, const Real sie, const Real P, const Real T,
+                             Lambda_t &&lambda, Real &dedP_T, Real &drdP_T, Real &dedT_P,
+                             Real &drdT_P) const;
+  */
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION void
   FillEos(Real &rho, Real &temp, Real &energy, Real &press, Real &cv, Real &bmod,
@@ -1005,6 +1013,42 @@ PORTABLE_INLINE_FUNCTION void SpinerEOSDependsRhoT::DensityEnergyFromPressureTem
   sie = InternalEnergyFromDensityTemperature(rho, temp, lambda);
 }
 
+/*
+// TODO(JMM): I would like to use this, but shockingly we don't have
+// dPdT stored. So need to add that to sesame2spiner first. Possibly
+// should be a different MR.
+template <typename Lambda_t>
+PORTABLE_INLINE_FUNCTION void SpinerEOSDependsRhoT::PTDerivativesFromPreferred(
+    const Real rho, const Real sie, const Real P, const Real T, Lambda_t &&lambda,
+    Real &dedP_T, Real &drdP_T, Real &dedT_P, Real &drdT_P) const {
+  const Real lT = lT_(T);
+  const Real lRho = lRho_(rho);
+  const TableStatus whereAmI = getLocDependsRhoT_(lRho, lT);
+
+  if (whereAmI == TableStatus::OffBottom) {
+    dedP_T = robust::ratio(1., dPdECold_.interpToReal(lRho));
+    drdP_T = robust::ratio(1., dPdRhoCold_.interpToReal(lRho));
+    dedT_P = dEdTCold_.interpToReal(lRho);
+    drdT_P = 0;                                 // TODO(JMM) is this right?
+  } else if (whereAmI == TableStatus::OffTop) { // idela gas
+    const Real gm1 = gm1Max_.interpToReal(lRho);
+    const Real Cv = dEdTMax_.interpToReal(lRho);
+    dedP_T = 0;
+    dedT_P = Cv;
+    drdP_T = robust::ratio(1.0, gm1 * sie);
+    drdT_P = -robust::ratio(P, gm1 * Cv * T * T);
+  } else {
+    dedP_T = robust::ratio(1., dPdE_.interpToReal(lRho, lT));
+    drdP_T = robust::ratio(1., dPdRho_.interpToReal(lRho, lT));
+    dedT_P = dEdT_.interpToReal(lRho, lT);
+    // dP = (dP/dr)_T dr + (dP/dT)_r dT
+    // dP = 0
+    // => 0 = (dP/dr)_T dr + (dP/dT)_r dT
+    drdT_P = -robust::ratio(dPdT_.interpToReal(lRho, lT), dPdRho_.interpToReal(lRho, lT));
+  }
+}
+*/
+
 template <typename Indexer_t>
 PORTABLE_INLINE_FUNCTION void
 SpinerEOSDependsRhoT::FillEos(Real &rho, Real &temp, Real &energy, Real &press, Real &cv,
diff --git a/singularity-eos/eos/eos_variant.hpp b/singularity-eos/eos/eos_variant.hpp
index 691b6ad91ca..463835ecc37 100644
--- a/singularity-eos/eos/eos_variant.hpp
+++ b/singularity-eos/eos/eos_variant.hpp
@@ -1,5 +1,5 @@
 //------------------------------------------------------------------------------
-// © 2021-2025. Triad National Security, LLC. All rights reserved.  This
+// © 2021-2026. 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.  Department of Energy/National
@@ -280,6 +280,20 @@ class Variant {
         eos_);
   }
 
+  // TODO(JMM): Do we need a vectorized version of this call?
+  template <typename Lambda_t = Real *>
+  PORTABLE_INLINE_FUNCTION void
+  PTDerivativesFromPreferred(const Real rho, const Real sie, const Real P, const Real T,
+                             Lambda_t &&lambda, Real &dedP_T, Real &drdP_T, Real &dedT_P,
+                             Real &drdT_P) const {
+    return PortsOfCall::visit(
+        [&](const auto &eos) {
+          return eos.PTDerivativesFromPreferred(rho, sie, P, T, lambda, dedP_T, drdP_T,
+                                                dedT_P, drdT_P);
+        },
+        eos_);
+  }
+
   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/modifiers/eos_unitsystem.hpp b/singularity-eos/eos/modifiers/eos_unitsystem.hpp
index 5154ea43120..a1b1bf4ffc6 100644
--- a/singularity-eos/eos/modifiers/eos_unitsystem.hpp
+++ b/singularity-eos/eos/modifiers/eos_unitsystem.hpp
@@ -214,6 +214,21 @@ class UnitSystem : public EosBase<UnitSystem<T>> {
     sie *= inv_sie_unit_;
   }
 
+  template <typename Lambda_t = Real *>
+  PORTABLE_INLINE_FUNCTION void
+  PTDerivativesFromPreferred(const Real rho, const Real sie, const Real P,
+                             const Real temp, Lambda_t &&lambda, Real &dedP_T,
+                             Real &drdP_T, Real &dedT_P, Real &drdT_P) const {
+    t_.PTDerivativesFromPreferred(rho * rho_unit_, sie * sie_unit_, P * press_unit_,
+                                  temp * temp_unit_, lambda, dedP_T, drdP_T, dedT_P,
+                                  drdT_P);
+    // scale outputs
+    dedP_T *= robust::ratio(press_unit_, sie_unit_);
+    drdP_T *= robust::ratio(press_unit_, rho_unit_);
+    dedT_P *= robust::ratio(temp_unit_, sie_unit_);
+    drdT_P *= robust::ratio(temp_unit_, rho_unit_);
+  }
+
   template <typename Indexer_t = Real *>
   PORTABLE_FUNCTION void FillEos(Real &rho, Real &temp, Real &energy, Real &press,
                                  Real &cv, Real &bmod, const unsigned long output,
diff --git a/singularity-eos/eos/modifiers/floored_energy.hpp b/singularity-eos/eos/modifiers/floored_energy.hpp
index 51debed081b..125c205f302 100644
--- a/singularity-eos/eos/modifiers/floored_energy.hpp
+++ b/singularity-eos/eos/modifiers/floored_energy.hpp
@@ -126,6 +126,17 @@ class FlooredEnergy : public EosBase<FlooredEnergy<T>> {
                                     Indexer_t &&lambda = nullptr) const {
     t_.InternalEnergyFromDensityPressure(rho, P, sie, lambda);
   }
+
+  template <typename Lambda_t = Real *>
+  PORTABLE_INLINE_FUNCTION void
+  PTDerivativesFromPreferred(const Real rho, const Real sie, const Real P,
+                             const Real temp, Lambda_t &&lambda, Real &dedP_T,
+                             Real &drdP_T, Real &dedT_P, Real &drdT_P) const {
+    const Real min_sie = t_.MinInternalEnergyFromDensity(rho);
+    t_.PTDerivativesFromPreferred(rho, std::max(sie, min_sie), P, temp, lambda, dedP_T,
+                                  drdP_T, dedT_P, drdT_P);
+  }
+
   template <typename Indexer_t = Real *>
   PORTABLE_FUNCTION void FillEos(Real &rho, Real &temp, Real &energy, Real &press,
                                  Real &cv, Real &bmod, const unsigned long output,
diff --git a/singularity-eos/eos/modifiers/ramps_eos.hpp b/singularity-eos/eos/modifiers/ramps_eos.hpp
index 4c9cdb866df..75a9bee2950 100644
--- a/singularity-eos/eos/modifiers/ramps_eos.hpp
+++ b/singularity-eos/eos/modifiers/ramps_eos.hpp
@@ -213,6 +213,9 @@ class BilinearRampEOS : public EosBase<BilinearRampEOS<T>> {
     // TODO(JMM): This doesn't seem right. dpdrho relevant.
     return t_.GruneisenParamFromDensityTemperature(rho, temperature, lambda);
   }
+
+  /* JMM: Use the base class PTDerivativesFromPreferred for now */
+
   template <typename Indexer_t = Real *>
   PORTABLE_FUNCTION void FillEos(Real &rho, Real &temp, Real &energy, Real &press,
                                  Real &cv, Real &bmod, const unsigned long output,
diff --git a/singularity-eos/eos/modifiers/relativistic_eos.hpp b/singularity-eos/eos/modifiers/relativistic_eos.hpp
index e74ee889260..00771536148 100644
--- a/singularity-eos/eos/modifiers/relativistic_eos.hpp
+++ b/singularity-eos/eos/modifiers/relativistic_eos.hpp
@@ -145,6 +145,14 @@ class RelativisticEOS : public EosBase<RelativisticEOS<T>> {
                                     Indexer_t &&lambda = nullptr) const {
     t_.InternalEnergyFromDensityPressure(rho, P, sie, lambda);
   }
+  template <typename Lambda_t = Real *>
+  PORTABLE_INLINE_FUNCTION void
+  PTDerivativesFromPreferred(const Real rho, const Real sie, const Real P,
+                             const Real temp, Lambda_t &&lambda, Real &dedP_T,
+                             Real &drdP_T, Real &dedT_P, Real &drdT_P) const {
+    t_.PTDerivativesFromPreferred(rho, sie, P, temp, lambda, dedP_T, drdP_T, dedT_P,
+                                  drdT_P);
+  }
   template <typename Indexer_t = Real *>
   PORTABLE_FUNCTION void FillEos(Real &rho, Real &temp, Real &energy, Real &press,
                                  Real &cv, Real &bmod, const unsigned long output,
diff --git a/singularity-eos/eos/modifiers/scaled_eos.hpp b/singularity-eos/eos/modifiers/scaled_eos.hpp
index 67a95ffc6c6..3e388854f55 100644
--- a/singularity-eos/eos/modifiers/scaled_eos.hpp
+++ b/singularity-eos/eos/modifiers/scaled_eos.hpp
@@ -141,6 +141,20 @@ class ScaledEOS : public EosBase<ScaledEOS<T>> {
                                     Indexer_t &&lambda = nullptr) const {
     t_.InternalEnergyFromDensityPressure(scale_ * rho, P, sie, lambda);
   }
+
+  template <typename Lambda_t = Real *>
+  PORTABLE_INLINE_FUNCTION void
+  PTDerivativesFromPreferred(const Real rho, const Real sie, const Real P,
+                             const Real temp, Lambda_t &&lambda, Real &dedP_T,
+                             Real &drdP_T, Real &dedT_P, Real &drdT_P) const {
+    t_.PTDerivativesFromPreferred(scale_ * rho, inv_scale_ * sie, P, temp, lambda, dedP_T,
+                                  drdP_T, dedT_P, drdT_P);
+    dedP_T *= scale_;
+    drdP_T *= inv_scale_;
+    dedT_P *= scale_;
+    drdT_P *= inv_scale_;
+  }
+
   template <typename Indexer_t = Real *>
   PORTABLE_FUNCTION void FillEos(Real &rho, Real &temp, Real &energy, Real &press,
                                  Real &cv, Real &bmod, const unsigned long output,
diff --git a/singularity-eos/eos/modifiers/shifted_eos.hpp b/singularity-eos/eos/modifiers/shifted_eos.hpp
index 3b636f4ab3e..1afb4c1a2e3 100644
--- a/singularity-eos/eos/modifiers/shifted_eos.hpp
+++ b/singularity-eos/eos/modifiers/shifted_eos.hpp
@@ -134,6 +134,16 @@ class ShiftedEOS : public EosBase<ShiftedEOS<T>> {
     t_.InternalEnergyFromDensityPressure(rho, P, sie, lambda);
     sie += shift_;
   }
+
+  template <typename Lambda_t = Real *>
+  PORTABLE_INLINE_FUNCTION void
+  PTDerivativesFromPreferred(const Real rho, const Real sie, const Real P,
+                             const Real temp, Lambda_t &&lambda, Real &dedP_T,
+                             Real &drdP_T, Real &dedT_P, Real &drdT_P) const {
+    t_.PTDerivativesFromPreferred(rho, sie - shift_, P, temp, lambda, dedP_T, drdP_T,
+                                  dedT_P, drdT_P);
+  }
+
   template <typename Indexer_t = Real *>
   PORTABLE_FUNCTION void FillEos(Real &rho, Real &temp, Real &energy, Real &press,
                                  Real &cv, Real &bmod, const unsigned long output,
diff --git a/singularity-eos/eos/modifiers/zsplit_eos.hpp b/singularity-eos/eos/modifiers/zsplit_eos.hpp
index a2d98b972a9..8602b2742ad 100644
--- a/singularity-eos/eos/modifiers/zsplit_eos.hpp
+++ b/singularity-eos/eos/modifiers/zsplit_eos.hpp
@@ -180,6 +180,21 @@ class ZSplit : public EosBase<ZSplit<ztype, T>> {
     sie *= scale;
   }
 
+  template <typename Lambda_t = Real *>
+  PORTABLE_INLINE_FUNCTION void
+  PTDerivativesFromPreferred(const Real rho, const Real sie, const Real P,
+                             const Real temp, Lambda_t &&lambda, Real &dedP_T,
+                             Real &drdP_T, Real &dedT_P, Real &drdT_P) const {
+    const Real scale = GetScale_(lambda);
+    const Real iscale = GetInvScale_(lambda);
+    t_.PTDerivativesFromPreferred(rho, sie * iscale, P * iscale, temp, lambda, dedP_T,
+                                  drdP_T, dedT_P, drdT_P);
+    // dedP scale cancels
+    drdP_T *= iscale;
+    dedT_P *= scale;
+    // drdT_P does not need scaling
+  }
+
   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/test/test_eos_ideal.cpp b/test/test_eos_ideal.cpp
index dd2f5ab8d0d..2faf1f6b33a 100644
--- a/test/test_eos_ideal.cpp
+++ b/test/test_eos_ideal.cpp
@@ -1,5 +1,5 @@
 //------------------------------------------------------------------------------
-// © 2021-2025. Triad National Security, LLC. All rights reserved.  This
+// © 2021-2026. 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.  Department of Energy/National
@@ -198,6 +198,53 @@ SCENARIO("Ideal gas density energy from prssure temperature",
           nwrong);
       THEN("There are no errors") { REQUIRE(nwrong == 0); }
     }
+    WHEN("We check PT derivatives from preferred") {
+      constexpr int N = 100;
+      constexpr Real derivative_eps = 3e-6;
+      int nwrong = 0;
+      portableReduce(
+          "Check PT derivatives from preferred", 1, N,
+          PORTABLE_LAMBDA(const int i, int &nw) {
+            Real rho = i;
+            Real T = 100 * i;
+            Real P = device_eos.PressureFromDensityTemperature(rho, T);
+            Real sie = device_eos.InternalEnergyFromDensityTemperature(rho, T);
+            Real dedP_T, drdP_T, dedT_P, drdT_P;
+            Real dedP_T_fd, drdP_T_fd, dedT_P_fd, drdT_P_fd;
+            device_eos.PTDerivativesFromPreferred(rho, sie, P, T,
+                                                  static_cast<Real *>(nullptr), dedP_T,
+                                                  drdP_T, dedT_P, drdT_P);
+            singularity::eos_base::PTDerivativesByFiniteDifferences(
+                device_eos, rho, sie, P, T, derivative_eps, static_cast<Real *>(nullptr),
+                dedP_T_fd, drdP_T_fd, dedT_P_fd, drdT_P_fd);
+            if (!isClose(dedP_T, dedP_T_fd)) {
+              printf("rho, sie, P, T, dedP_T expected, true, diff:\n\t"
+                     "%.14e %.14e %.14e %.14e %.14e %.14e %.14e\n",
+                     rho, sie, P, T, dedP_T, dedP_T_fd, dedP_T - dedP_T_fd);
+              nw += 1;
+            }
+            if (!isClose(drdP_T, drdP_T_fd)) {
+              printf("rho, sie, P, T, drdP_T expected, true, diff:\n\t"
+                     "%.14e %.14e %.14e %.14e %.14e %.14e %.14e\n",
+                     rho, sie, P, T, drdP_T, drdP_T_fd, drdP_T - drdP_T_fd);
+              nw += 1;
+            }
+            if (!isClose(dedT_P, dedT_P_fd)) {
+              printf("rho, sie, P, T, dedT_P expected, true, diff:\n\t"
+                     "%.14e %.14e %.14e %.14e %.14e %.14e %.14e\n",
+                     rho, sie, P, T, dedT_P, dedT_P_fd, dedT_P - dedT_P_fd);
+              nw += 1;
+            }
+            if (!isClose(drdT_P, drdT_P_fd)) {
+              printf("rho, sie, P, T, drdT_P expected, true, diff:\n\t"
+                     "%.14e %.14e %.14e %.14e %.14e %.14e %.14e\n",
+                     rho, sie, P, T, drdT_P, drdT_P_fd, drdT_P - drdT_P_fd);
+              nw += 1;
+            }
+          },
+          nwrong);
+      THEN("There are no errors") { REQUIRE(nwrong == 0); }
+    }
   }
 }
 
diff --git a/test/test_eos_modifiers.cpp b/test/test_eos_modifiers.cpp
index 7e7ea9ba055..d37252633ec 100644
--- a/test/test_eos_modifiers.cpp
+++ b/test/test_eos_modifiers.cpp
@@ -156,6 +156,63 @@ SCENARIO("EOS Builder and Modifiers", "[EOSBuilder][Modifiers][IdealGas]") {
           eos_ramped.PrintParams();
         }
       }
+
+      WHEN("We compute PT derivatives from preferred for the modified EOS") {
+        // Choose a representative state
+        constexpr Real test_rho = 1.5;
+        constexpr Real test_T = 2.0;
+
+        // Compute primitive quantities
+        const Real P = eos.PressureFromDensityTemperature(test_rho, test_T);
+        const Real sie_val = eos.InternalEnergyFromDensityTemperature(test_rho, test_T);
+
+        // Derivatives from the EOS preferred interface
+        Real dedP_T, drdP_T, dedT_P, drdT_P;
+        eos.PTDerivativesFromPreferred(test_rho, sie_val, P, test_T,
+                                       static_cast<Real *>(nullptr), dedP_T, drdP_T,
+                                       dedT_P, drdT_P);
+
+        // Finite‑difference reference
+        constexpr Real derivative_eps = 3e-6;
+        Real dedP_T_fd, drdP_T_fd, dedT_P_fd, drdT_P_fd;
+        singularity::eos_base::PTDerivativesByFiniteDifferences(
+            eos, test_rho, sie_val, P, test_T, derivative_eps,
+            static_cast<Real *>(nullptr), dedP_T_fd, drdP_T_fd, dedT_P_fd, drdT_P_fd);
+
+        // Compare each component
+        THEN("dedP_T matches finite‑difference reference") {
+          bool close = isClose(dedP_T, dedP_T_fd);
+          if (!close) {
+            printf("dedP_T mismatch: eos = %.14e, fd = %.14e, diff = %.14e\n", dedP_T,
+                   dedP_T_fd, dedP_T - dedP_T_fd);
+          }
+          REQUIRE(close);
+        }
+        THEN("drdP_T matches finite‑difference reference") {
+          bool close = isClose(drdP_T, drdP_T_fd);
+          if (!close) {
+            printf("drdP_T mismatch: eos = %.14e, fd = %.14e, diff = %.14e\n", drdP_T,
+                   drdP_T_fd, drdP_T - drdP_T_fd);
+          }
+          REQUIRE(close);
+        }
+        THEN("dedT_P matches finite‑difference reference") {
+          bool close = isClose(dedT_P, dedT_P_fd);
+          if (!close) {
+            printf("dedT_P mismatch: eos = %.14e, fd = %.14e, diff = %.14e\n", dedT_P,
+                   dedT_P_fd, dedT_P - dedT_P_fd);
+          }
+          REQUIRE(close);
+        }
+        THEN("drdT_P matches finite‑difference reference") {
+          bool close = isClose(drdT_P, drdT_P_fd);
+          if (!close) {
+            printf("drdT_P mismatch: eos = %.14e, fd = %.14e, diff = %.14e\n", drdT_P,
+                   drdT_P_fd, drdT_P - drdT_P_fd);
+          }
+          REQUIRE(close);
+        }
+      }
     }
 #ifdef SINGULARITY_BUILD_CLOSURE
     WHEN("We construct a non-modifying modifier") {
@@ -204,56 +261,53 @@ SCENARIO("EOS Builder and Modifiers", "[EOSBuilder][Modifiers][IdealGas]") {
       // bmod (rho_t2)
       const Real bmodrt2 = rho_t2 * b / r0;
       THEN("P_eos(alpha_0*rmid, T0) = P_ramp(rmid,T0)") {
-        INFO("P_eos(alpha_0*rmid, T0): "
-             << P_armid
-             << " P_ramp(rmid, T0): " << igra.PressureFromDensityTemperature(rmid, T0));
+        printf("P_eos(alpha_0*rmid, T0): %.14e  P_ramp(rmid, T0): %.14e\n", P_armid,
+               igra.PressureFromDensityTemperature(rmid, T0));
         REQUIRE(isClose(P_armid, igra.PressureFromDensityTemperature(rmid, T0), 1.e-12));
       }
       THEN("We obtain correct ramp behavior in P(rho) for rho <r0, [r0,rmid], [rmid,r1] "
            "and >r1") {
         // also check pressures on ramp
-        INFO("reference P((r0+rmid)/2, T0): "
-             << Prhot1 << " test P((r0+rmid)/2, T0): "
-             << igra.PressureFromDensityTemperature(rho_t1, T0));
+        printf("reference P((r0+rmid)/2, T0): %.14e  test P((r0+rmid)/2, T0): %.14e\n",
+               Prhot1, igra.PressureFromDensityTemperature(rho_t1, T0));
         REQUIRE(isClose(Prhot1, igra.PressureFromDensityTemperature(rho_t1, T0), 1.e-12));
-        INFO("reference P((rmid+r1)/2, T0): "
-             << Prhot2 << " test P((rmid+r1)/2, T0): "
-             << igra.PressureFromDensityTemperature(rho_t2, T0));
+        printf("reference P((rmid+r1)/2, T0): %.14e  test P((rmid+r1)/2, T0): %.14e\n",
+               Prhot2, igra.PressureFromDensityTemperature(rho_t2, T0));
         REQUIRE(isClose(Prhot2, igra.PressureFromDensityTemperature(rho_t2, T0), 1.e-12));
         // check pressure below and beyond ramp matches unmodified ideal gas
-        INFO("reference P(0.8*r0, T0): "
-             << 0.8 * r0 * gm1 * Cv * T0 << " test P(0.8*r0, T0): "
-             << igra.PressureFromDensityTemperature(0.8 * r0, T0));
+        printf("reference P(0.8*r0, T0): %.14e  test P(0.8*r0, T0): %.14e\n",
+               0.8 * r0 * gm1 * Cv * T0,
+               igra.PressureFromDensityTemperature(0.8 * r0, T0));
         REQUIRE(isClose(0.8 * r0 * gm1 * Cv * T0,
                         igra.PressureFromDensityTemperature(0.8 * r0, T0), 1.e-12));
-        INFO("reference P(1.2*r1, T0): "
-             << 1.2 * r1 * gm1 * Cv * T0 << " test P(1.2*r1, T0): "
-             << igra.PressureFromDensityTemperature(1.2 * r1, T0));
-        REQUIRE(isClose(1.2 * r1 * gm1 * Cv * T0,
-                        igra.PressureFromDensityTemperature(1.2 * r1, T0), 1.e-12));
+        printf("reference P(1.2*r1, T0): %.14e  test P(1.2*r1, T0): %.14e\n",
+               1.2 * r0 * gm1 * Cv * T0,
+               igra.PressureFromDensityTemperature(1.2 * r0, T0));
+        REQUIRE(isClose(1.2 * r0 * gm1 * Cv * T0,
+                        igra.PressureFromDensityTemperature(1.2 * r0, T0), 1.e-12));
       }
       THEN("We obtain correct ramp behavior in bmod(rho) for rho <r0, [r0,rmid], "
            "[rmid,r1] and >r0") {
         // check bulk moduli on both pieces of ramp
-        INFO("reference bmod((r0+rmid)/2, T0): "
-             << bmodrt1 << " test bmod((r0+rmid)/2, T0): "
-             << igra.BulkModulusFromDensityTemperature(rho_t1, T0));
+        printf(
+            "reference bmod((r0+rmid)/2, T0): %.14e  test bmod((r0+rmid)/2, T0): %.14e\n",
+            bmodrt1, igra.BulkModulusFromDensityTemperature(rho_t1, T0));
         REQUIRE(
             isClose(bmodrt1, igra.BulkModulusFromDensityTemperature(rho_t1, T0), 1.e-12));
-        INFO("reference bmod((rmid+r1)/2, T0): "
-             << bmodrt2 << " test bmod((rmid+r1)/2, T0): "
-             << igra.BulkModulusFromDensityTemperature(rho_t2, T0));
+        printf(
+            "reference bmod((rmid+r1)/2, T0): %.14e  test bmod((rmid+r1)/2, T0): %.14e\n",
+            bmodrt2, igra.BulkModulusFromDensityTemperature(rho_t2, T0));
         REQUIRE(
             isClose(bmodrt2, igra.BulkModulusFromDensityTemperature(rho_t2, T0), 1.e-12));
         // check bulk modulus below and beyond ramp matches unmodified ideal gas
-        INFO("reference bmod(0.8*r0, T0): "
-             << 0.8 * r0 * gm1 * (gm1 + 1.0) * Cv * T0 << " test bmod(0.8*r0, T0): "
-             << igra.BulkModulusFromDensityTemperature(0.8 * r0, T0));
+        printf("reference bmod(0.8*r0, T0): %.14e  test bmod(0.8*r0, T0): %.14e\n",
+               0.8 * r0 * gm1 * (gm1 + 1.0) * Cv * T0,
+               igra.BulkModulusFromDensityTemperature(0.8 * r0, T0));
         REQUIRE(isClose(0.8 * r0 * gm1 * (gm1 + 1.0) * Cv * T0,
                         igra.BulkModulusFromDensityTemperature(0.8 * r0, T0), 1.e-12));
-        INFO("reference bmod(1.2*r1, T0): "
-             << 1.2 * r1 * gm1 * (gm1 + 1.0) * Cv * T0 << " test bmod(1.2*r1, T0): "
-             << igra.BulkModulusFromDensityTemperature(1.2 * r1, T0));
+        printf("reference bmod(1.2*r1, T0): %.14e  test bmod(1.2*r1, T0): %.14e\n",
+               1.2 * r1 * gm1 * (gm1 + 1.0) * Cv * T0,
+               igra.BulkModulusFromDensityTemperature(1.2 * r1, T0));
         REQUIRE(isClose(1.2 * r1 * gm1 * (gm1 + 1.0) * Cv * T0,
                         igra.BulkModulusFromDensityTemperature(1.2 * r1, T0), 1.e-12));
       }
diff --git a/test/test_eos_tabulated.cpp b/test/test_eos_tabulated.cpp
index 0ad2c0d32ab..e61fb6b8185 100644
--- a/test/test_eos_tabulated.cpp
+++ b/test/test_eos_tabulated.cpp
@@ -97,6 +97,45 @@ SCENARIO("SpinerEOS depends on Rho and T", "[SpinerEOS][DependsRhoT][EOSPAC]") {
     auto steelEOS_host = steelEOS_host_polymorphic.get<SpinerEOSDependsRhoT>();
 
     EOS eospac = EOSPAC(steelID);
+    THEN("The EOSPAC model can provide derivatives from preferred") {
+      Real rho = 50;
+      Real T = 500;
+      Real P = eospac.PressureFromDensityTemperature(rho, T);
+      Real sie = eospac.InternalEnergyFromDensityTemperature(rho, T);
+
+      constexpr Real derivative_eps = 1e-8;
+      Real dedP_T, drdP_T, dedT_P, drdT_P;
+      Real dedP_T_fd, drdP_T_fd, dedT_P_fd, drdT_P_fd;
+      eospac.PTDerivativesFromPreferred(rho, sie, P, T, static_cast<Real *>(nullptr),
+                                        dedP_T, drdP_T, dedT_P, drdT_P);
+      singularity::eos_base::PTDerivativesByFiniteDifferences(
+          eospac, rho, sie, P, T, derivative_eps, static_cast<Real *>(nullptr), dedP_T_fd,
+          drdP_T_fd, dedT_P_fd, drdT_P_fd);
+      bool is_close_dedp = isClose(dedP_T, dedP_T_fd);
+      if (!is_close_dedp) {
+        printf("dedP_T not close! %.14e %.14e %.14e\n", dedP_T, dedP_T_fd,
+               dedP_T - dedP_T_fd);
+      }
+      REQUIRE(is_close_dedp);
+      bool is_close_drdP = isClose(drdP_T, drdP_T_fd);
+      if (!is_close_drdP) {
+        printf("drdP_T not close! %.14e %.14e %.14e\n", drdP_T, drdP_T_fd,
+               drdP_T - drdP_T_fd);
+      }
+      REQUIRE(is_close_drdP);
+      bool is_close_dedT = isClose(dedT_P, dedT_P_fd);
+      if (!is_close_dedT) {
+        printf("dedT_P not close! %.14e %.14e %.14e\n", dedT_P, dedT_P_fd,
+               dedT_P - dedT_P_fd);
+      }
+      REQUIRE(is_close_dedT);
+      bool is_close_drdT = isClose(drdT_P, drdT_P_fd);
+      if (!is_close_drdT) {
+        printf("drdT_P not close! %.14e %.14e %.14e\n", drdT_P, drdT_P_fd,
+               drdT_P - drdT_P_fd);
+      }
+      REQUIRE(is_close_drdT);
+    }
 
     THEN("The correct metadata is read in") {
       REQUIRE(steelEOS_host.matid() == steelID);