diff --git a/CHANGELOG.md b/CHANGELOG.md
index a9f34c841bc..477c17796b4 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -3,6 +3,7 @@
 ## Current develop
 
 ### Added (new features/APIs/variables/...)
+- [[PR589]](https://github.com/lanl/singularity-eos/pull/589) InternalEnergyFromDensityPressure
 
 ### Fixed (Repair bugs, etc)
 
diff --git a/doc/sphinx/src/using-eos.rst b/doc/sphinx/src/using-eos.rst
index 326c4173367..ed42b204eb1 100644
--- a/doc/sphinx/src/using-eos.rst
+++ b/doc/sphinx/src/using-eos.rst
@@ -1181,6 +1181,25 @@ given density in :math:`g/cm^3` and temperature in Kelvin.
 returns the unitless Gruneisen parameter given density in
 :math:`g/cm^3` and specific internal energy in :math:`erg/g`.
 
+.. code-block:: cpp
+
+   template <typename Indexer_t = Real*>
+   Real InternalEnergyFromDensityPressure(const Real rho, const Real P, 
+                                          Indexer_t &&lambda = nullptr) const;
+
+returns the specific internal energy in :math:`erg/g` given a density
+in :math:`g/cm^3` and a pressure in Barye. This call is often a root
+find, and thus an alternative signature allows an initial guess to be
+passed in by reference:
+
+.. code-block:: cpp
+
+   template <typename Indexer_t = Real*>
+   void InternalEnergyFromDensityPressure(const Real rho, const Real P, Real &sie,
+                                          Indexer_t &&lambda = nullptr) const;
+
+In this case, the guess will be set to the new value by the function.
+
 .. code-block:: cpp
 
    template <typename Indexer_t = Real*>
diff --git a/python/module.hpp b/python/module.hpp
index 5520ddf25cd..278cff68277 100644
--- a/python/module.hpp
+++ b/python/module.hpp
@@ -169,6 +169,7 @@ EOS_VEC_FUNC_TMPL(BulkModulusFromDensityTemperature, rhos, temperatures, bmods)
 EOS_VEC_FUNC_TMPL(BulkModulusFromDensityInternalEnergy, rhos, sies, bmods)
 EOS_VEC_FUNC_TMPL(GruneisenParamFromDensityTemperature, rhos, temperatures, gm1s)
 EOS_VEC_FUNC_TMPL(GruneisenParamFromDensityInternalEnergy, rhos, sies, gm1s)
+EOS_VEC_FUNC_TMPL(InternalEnergyFromDensityPressure, rhos, Ps, sies)
 
 struct EOSState {
   Real density;
@@ -229,6 +230,7 @@ struct VectorFunctions {
     .def("BulkModulusFromDensityInternalEnergy", &BulkModulusFromDensityInternalEnergy<T>, py::arg("rhos"), py::arg("sies"), py::arg("bmods"), py::arg("lmbdas"))
     .def("GruneisenParamFromDensityTemperature", &GruneisenParamFromDensityTemperature<T>, py::arg("rhos"), py::arg("temperatures"), py::arg("gm1s"), py::arg("lmbdas"))
     .def("GruneisenParamFromDensityInternalEnergy", &GruneisenParamFromDensityInternalEnergy<T>, py::arg("rhos"), py::arg("sies"), py::arg("gm1s"), py::arg("lmbdas"))
+    .def("InternalEnergyFromDensityPressure", &InternalEnergyFromDensityPressure<T>, py::arg("rhos"), py::arg("Ps"), py::arg("sies"), py::arg("lmbdas"))
 
     .def("TemperatureFromDensityInternalEnergy", &TemperatureFromDensityInternalEnergyNoLambda<T>, py::arg("rhos"), py::arg("sies"), py::arg("temperatures"))
     .def("InternalEnergyFromDensityTemperature", &InternalEnergyFromDensityTemperatureNoLambda<T>, py::arg("rhos"), py::arg("temperatures"), py::arg("sies"))
@@ -240,7 +242,7 @@ struct VectorFunctions {
     .def("BulkModulusFromDensityInternalEnergy", &BulkModulusFromDensityInternalEnergyNoLambda<T>, py::arg("rhos"), py::arg("sies"), py::arg("bmods"))
     .def("GruneisenParamFromDensityTemperature", &GruneisenParamFromDensityTemperatureNoLambda<T>, py::arg("rhos"), py::arg("temperatures"), py::arg("gm1s"))
     .def("GruneisenParamFromDensityInternalEnergy", &GruneisenParamFromDensityInternalEnergyNoLambda<T>, py::arg("rhos"), py::arg("sies"), py::arg("gm1s"))
-
+    .def("InternalEnergyFromDensityPressure", &InternalEnergyFromDensityPressureNoLambda<T>, py::arg("rhos"), py::arg("Ps"), py::arg("sies"))
 
     .def("FillEos", [](const T & self, py::array_t<Real> rhos,
                        py::array_t<Real> temperatures, py::array_t<Real> sies,
@@ -300,6 +302,7 @@ struct VectorFunctions<T,true> {
     .def("BulkModulusFromDensityInternalEnergy", &BulkModulusFromDensityInternalEnergyWithScratch<T>, py::arg("rhos"), py::arg("sies"), py::arg("bmods"), py::arg("scratch"), py::arg("lmbdas"))
     .def("GruneisenParamFromDensityTemperature", &GruneisenParamFromDensityTemperatureWithScratch<T>, py::arg("rhos"), py::arg("temperatures"), py::arg("gm1s"), py::arg("scratch"), py::arg("lmbdas"))
     .def("GruneisenParamFromDensityInternalEnergy", &GruneisenParamFromDensityInternalEnergyWithScratch<T>, py::arg("rhos"), py::arg("sies"), py::arg("gm1s"), py::arg("scratch"), py::arg("lmbdas"))
+    .def("InternalEnergyFromDensityPressure", &InternalEnergyFromDensityPressureWithScratch<T>, py::arg("rhos"), py::arg("Ps"), py::arg("sies"), py::arg("scratch"), py::arg("lmbdas"))
 
     .def("TemperatureFromDensityInternalEnergy", &TemperatureFromDensityInternalEnergyNoLambdaWithScratch<T>, py::arg("rhos"), py::arg("sies"), py::arg("temperatures"), py::arg("scratch"))
     .def("InternalEnergyFromDensityTemperature", &InternalEnergyFromDensityTemperatureNoLambdaWithScratch<T>, py::arg("rhos"), py::arg("temperatures"), py::arg("sies"), py::arg("scratch"))
@@ -311,6 +314,7 @@ struct VectorFunctions<T,true> {
     .def("BulkModulusFromDensityInternalEnergy", &BulkModulusFromDensityInternalEnergyNoLambdaWithScratch<T>, py::arg("rhos"), py::arg("sies"), py::arg("bmods"), py::arg("scratch"))
     .def("GruneisenParamFromDensityTemperature", &GruneisenParamFromDensityTemperatureNoLambdaWithScratch<T>, py::arg("rhos"), py::arg("temperatures"), py::arg("gm1s"), py::arg("scratch"))
     .def("GruneisenParamFromDensityInternalEnergy", &GruneisenParamFromDensityInternalEnergyNoLambdaWithScratch<T>, py::arg("rhos"), py::arg("sies"), py::arg("gm1s"), py::arg("scratch"))
+    .def("InternalEnergyFromDensityPressure", &InternalEnergyFromDensityPressureNoLambdaWithScratch<T>, py::arg("rhos"), py::arg("Ps"), py::arg("sies"), py::arg("scratch"))
 
     .def("FillEos", [](const T & self, py::array_t<Real> rhos,
                        py::array_t<Real> temperatures, py::array_t<Real> sies,
@@ -452,6 +456,14 @@ py::class_<T> eos_class(py::module_ & m, std::string name) {
       self.DensityEnergyFromPressureTemperature(press, temp, np<Real>(), rho, sie);
       return std::pair<Real, Real>(rho, sie);
     }, py::arg("press"), py::arg("temp"))
+    .def("InternalEnergyFromDensityPressure", [](const T & self, const Real rho, const Real P, Real sie) {
+      self.InternalEnergyFromDensityPressure(rho, P, sie, np<Real>());
+      return sie;
+    }, py::arg("rho"), py::arg("P"), py::arg("sie"))
+    .def("InternalEnergyFromDensityPressure", [](const T & self, const Real rho, const Real P, Real sie, py::array_t<Real> lambda) {
+      self.InternalEnergyFromDensityPressure(rho, P, sie, lambda.mutable_data());
+      return sie;
+    }, py::arg("rho"), py::arg("P"), py::arg("sie"), py::arg("lmbda"))
     .def("Finalize", &T::Finalize)
     .def_property_readonly_static("EosType", [](py::object) { return T::EosType(); })
     .def_static("scratch_size", &T::scratch_size, py::arg("method_name"), py::arg("nelements"))
diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp
index 739dbc62764..2e8581454b4 100644
--- a/singularity-eos/eos/eos_base.hpp
+++ b/singularity-eos/eos/eos_base.hpp
@@ -82,6 +82,7 @@ char *StrCat(char *destination, const char *source) {
   using EosBase<__VA_ARGS__>::BulkModulusFromDensityInternalEnergy;                      \
   using EosBase<__VA_ARGS__>::GruneisenParamFromDensityTemperature;                      \
   using EosBase<__VA_ARGS__>::GruneisenParamFromDensityInternalEnergy;                   \
+  using EosBase<__VA_ARGS__>::InternalEnergyFromDensityPressure;                         \
   using EosBase<__VA_ARGS__>::FillEos;                                                   \
   using EosBase<__VA_ARGS__>::EntropyFromDensityTemperature;                             \
   using EosBase<__VA_ARGS__>::EntropyFromDensityInternalEnergy;                          \
@@ -928,6 +929,76 @@ class EosBase {
                                               rho, sie);
   }
 
+  // JMM: Another set of calls that are often overloaded for special cases
+  // TODO(JMM): Do we also want TemperatureFromDensityPressure? That's
+  // the more likely fundamental call, but the less likely "useful"
+  // call...
+  // We pass in sie by value to allow for initial guesses
+  template <typename Indexer_t = Real *>
+  PORTABLE_INLINE_FUNCTION void InternalEnergyFromDensityPressure(
+      const Real rho, const Real P, Real &sie,
+      Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
+    using RootFinding1D::regula_falsi;
+    using RootFinding1D::Status;
+    const CRTP &eos = *(static_cast<CRTP const *>(this));
+    auto f = [&](const Real sie) {
+      return eos.PressureFromDensityInternalEnergy(rho, sie, lambda);
+    };
+    const Real sie_min =
+        eos.InternalEnergyFromDensityTemperature(rho, eos.MinimumTemperature());
+    // temp not bounded. just pick something huge.
+    const Real sie_max = eos.InternalEnergyFromDensityTemperature(rho, 1e20);
+    Real sie_guess =
+        (((sie_min < sie) && (sie < sie_max)) ? 0.5 * (sie_min + sie_max) : sie);
+    auto status = regula_falsi(f, P, sie_guess, sie_min, sie_max, robust::EPS(),
+                               robust::EPS(), sie);
+    if (status == Status::FAIL) {
+      PORTABLE_WARN("InternalEnergyFromDensityPressure: failed to find root\n");
+    }
+  }
+  // This version foregoes the initial guess. You get what you get.
+  template <typename Indexer_t = Real *>
+  PORTABLE_INLINE_FUNCTION Real InternalEnergyFromDensityPressure(
+      const Real rho, const Real P,
+      Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
+    const CRTP &eos = *(static_cast<CRTP const *>(this));
+    const Real sie_min =
+        eos.InternalEnergyFromDensityTemperature(rho, eos.MinimumTemperature());
+    // temp not bounded. just pick something huge.
+    const Real sie_max = eos.InternalEnergyFromDensityTemperature(rho, 1e20);
+    Real sie = 0.5 * (sie_min + sie_max);
+    eos.InternalEnergyFromDensityPressure(rho, P, sie, lambda);
+    return sie;
+  }
+  // Classes like EOSPAC probably wish to overload/shadow this version
+  template <typename RealIndexer, typename ConstRealIndexer, typename LambdaIndexer>
+  inline void InternalEnergyFromDensityPressure(ConstRealIndexer &&rhos,
+                                                ConstRealIndexer &&Ps, RealIndexer &&sies,
+                                                const int num,
+                                                LambdaIndexer &&lambdas) const {
+    static auto const name = SG_MEMBER_FUNC_NAME();
+    static auto const cname = name.c_str();
+    CRTP copy = *(static_cast<CRTP const *>(this));
+    portableFor(
+        cname, 0, num, PORTABLE_LAMBDA(const int i) {
+          copy.InternalEnergyFromDensityPressure(rhos[i], Ps[i], sies[i], lambdas[i]);
+        });
+  }
+  template <typename RealIndexer, typename ConstRealIndexer, typename LambdaIndexer>
+  inline void InternalEnergyFromDensityPressure(ConstRealIndexer &&rhos,
+                                                ConstRealIndexer &&Ps, RealIndexer &&sies,
+                                                Real * /*scratch */, const int num,
+                                                LambdaIndexer &&lambdas) const {
+    return InternalEnergyFromDensityPressure(rhos, Ps, sies, num, lambdas);
+  }
+  template <typename LambdaIndexer>
+  inline void InternalEnergyFromDensityPressure(const Real *rhos, const Real *Ps,
+                                                Real *sies, Real * /*scratch */,
+                                                const int num, LambdaIndexer &&lambdas,
+                                                Transform && = Transform()) const {
+    return InternalEnergyFromDensityPressure(rhos, Ps, sies, num, lambdas);
+  }
+
   // Serialization
   /*
     The methodology here is there are *three* size methods all EOS's provide:
diff --git a/singularity-eos/eos/eos_eospac.hpp b/singularity-eos/eos/eos_eospac.hpp
index 032a5334b13..246144083ed 100644
--- a/singularity-eos/eos/eos_eospac.hpp
+++ b/singularity-eos/eos/eos_eospac.hpp
@@ -494,6 +494,8 @@ class EOSPAC : public EosBase<EOSPAC> {
 
   // TODO(JMM): Add performant Gibbs Free Energy
   using EosBase<EOSPAC>::FillEos;
+  // TODO(JMM): Add a performant internal energy call
+  using EosBase<EOSPAC>::InternalEnergyFromDensityPressure;
 
   SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(AZbar_)
 
diff --git a/singularity-eos/eos/eos_ideal.hpp b/singularity-eos/eos/eos_ideal.hpp
index 364bca4aec5..73c6bea53a3 100644
--- a/singularity-eos/eos/eos_ideal.hpp
+++ b/singularity-eos/eos/eos_ideal.hpp
@@ -186,6 +186,14 @@ class IdealGas : public EosBase<IdealGas> {
     sie = MYMAX(0.0, _Cv * temp);
     rho = MYMAX(0.0, press / (_gm1 * sie));
   }
+  template <typename Indexer_t = Real *>
+  PORTABLE_INLINE_FUNCTION void InternalEnergyFromDensityPressure(
+      const Real rho, const Real P, Real &sie,
+      Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
+    // P = gm1 rho sie
+    sie = P / (_gm1 * rho);
+  }
+
   inline void Finalize() {}
   static std::string EosType() { return std::string("IdealGas"); }
   static std::string EosPyType() { return EosType(); }
diff --git a/singularity-eos/eos/eos_variant.hpp b/singularity-eos/eos/eos_variant.hpp
index 232d531ec7a..d04b7c3674f 100644
--- a/singularity-eos/eos/eos_variant.hpp
+++ b/singularity-eos/eos/eos_variant.hpp
@@ -336,6 +336,17 @@ class Variant {
         eos_);
   }
 
+  template <typename Indexer_t = Real *>
+  PORTABLE_INLINE_FUNCTION void InternalEnergyFromDensityPressure(
+      const Real rho, const Real P, Real &sie,
+      Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
+    return PortsOfCall::visit(
+        [&rho, &P, &sie, &lambda](const auto &eos) {
+          return eos.InternalEnergyFromDensityPressure(rho, P, sie, lambda);
+        },
+        eos_);
+  }
+
   PORTABLE_INLINE_FUNCTION
   Real RhoPmin(const Real temp) const {
     return PortsOfCall::visit([&temp](const auto &eos) { return eos.RhoPmin(temp); },
@@ -1194,6 +1205,55 @@ class Variant {
         eos_);
   }
 
+  template <typename RealIndexer, typename ConstRealIndexer>
+  inline void InternalEnergyFromDensityPressure(ConstRealIndexer &&rhos,
+                                                ConstRealIndexer &&Ps, RealIndexer &&sies,
+                                                const int num) const {
+    NullIndexer lambdas{}; // Returns null pointer for every index
+    return InternalEnergyFromDensityPressure(
+        std::forward<ConstRealIndexer>(rhos), std::forward<ConstRealIndexer>(Ps),
+        std::forward<RealIndexer>(sies), num, lambdas);
+  }
+
+  template <typename RealIndexer, typename ConstRealIndexer, typename LambdaIndexer>
+  inline void InternalEnergyFromDensityPressure(ConstRealIndexer &&rhos,
+                                                ConstRealIndexer &&Ps, RealIndexer &&sies,
+                                                const int num,
+                                                LambdaIndexer &&lambdas) const {
+    return PortsOfCall::visit(
+        [&](const auto &eos) {
+          return eos.InternalEnergyFromDensityPressure(
+              std::forward<ConstRealIndexer>(rhos), std::forward<ConstRealIndexer>(Ps),
+              std::forward<RealIndexer>(sies), num, std::forward<LambdaIndexer>(lambdas));
+        },
+        eos_);
+  }
+
+  template <typename RealIndexer, typename ConstRealIndexer>
+  inline void InternalEnergyFromDensityPressure(ConstRealIndexer &&rhos,
+                                                ConstRealIndexer &&Ps, RealIndexer &&sies,
+                                                Real *scratch, const int num) const {
+    NullIndexer lambdas{}; // Returns null pointer for every index
+    return InternalEnergyFromDensityPressure(
+        std::forward<ConstRealIndexer>(rhos), std::forward<ConstRealIndexer>(Ps),
+        std::forward<RealIndexer>(sies), scratch, num, lambdas);
+  }
+
+  template <typename RealIndexer, typename ConstRealIndexer, typename LambdaIndexer>
+  inline void InternalEnergyFromDensityPressure(ConstRealIndexer &&rhos,
+                                                ConstRealIndexer &&Ps, RealIndexer &&sies,
+                                                Real *scratch, const int num,
+                                                LambdaIndexer &&lambdas) const {
+    return PortsOfCall::visit(
+        [&rhos, &Ps, &sies, &scratch, &num, &lambdas](const auto &eos) {
+          return eos.InternalEnergyFromDensityPressure(
+              std::forward<ConstRealIndexer>(rhos), std::forward<ConstRealIndexer>(Ps),
+              std::forward<RealIndexer>(sies), scratch, num,
+              std::forward<LambdaIndexer>(lambdas));
+        },
+        eos_);
+  }
+
   template <typename RealIndexer>
   inline void FillEos(RealIndexer &&rhos, RealIndexer &&temps, RealIndexer &&energies,
                       RealIndexer &&presses, RealIndexer &&cvs, RealIndexer &&bmods,
diff --git a/singularity-eos/eos/modifiers/eos_unitsystem.hpp b/singularity-eos/eos/modifiers/eos_unitsystem.hpp
index 98bf594ef84..8acbf19025b 100644
--- a/singularity-eos/eos/modifiers/eos_unitsystem.hpp
+++ b/singularity-eos/eos/modifiers/eos_unitsystem.hpp
@@ -193,6 +193,14 @@ class UnitSystem : public EosBase<UnitSystem<T>> {
                                                              temp * temp_unit_, lambda);
     return gm1;
   }
+  template <typename Indexer_t = Real *>
+  PORTABLE_FUNCTION void
+  InternalEnergyFromDensityPressure(const Real rho, const Real P, Real &sie,
+                                    Indexer_t &&lambda = nullptr) const {
+    sie *= sie_unit_;
+    t_.InternalEnergyFromDensityPressure(rho_unit_ * rho, press_unit_ * P, sie, lambda);
+    sie *= inv_sie_unit_;
+  }
 
   template <typename Indexer_t = Real *>
   PORTABLE_FUNCTION void
diff --git a/singularity-eos/eos/modifiers/floored_energy.hpp b/singularity-eos/eos/modifiers/floored_energy.hpp
index b3b926df57c..314cd33a432 100644
--- a/singularity-eos/eos/modifiers/floored_energy.hpp
+++ b/singularity-eos/eos/modifiers/floored_energy.hpp
@@ -121,6 +121,12 @@ class FlooredEnergy : public EosBase<FlooredEnergy<T>> {
     return t_.GruneisenParamFromDensityTemperature(rho, temperature, lambda);
   }
   template <typename Indexer_t = Real *>
+  PORTABLE_FUNCTION void
+  InternalEnergyFromDensityPressure(const Real rho, const Real P, Real &sie,
+                                    Indexer_t &&lambda = nullptr) const {
+    t_.InternalEnergyFromDensityPressure(rho, P, sie, lambda);
+  }
+  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,
                                  Indexer_t &&lambda = nullptr) const {
diff --git a/singularity-eos/eos/modifiers/relativistic_eos.hpp b/singularity-eos/eos/modifiers/relativistic_eos.hpp
index fafa2f95a62..8863a7295a9 100644
--- a/singularity-eos/eos/modifiers/relativistic_eos.hpp
+++ b/singularity-eos/eos/modifiers/relativistic_eos.hpp
@@ -140,6 +140,12 @@ class RelativisticEOS : public EosBase<RelativisticEOS<T>> {
     return t_.GruneisenParamFromDensityTemperature(rho, temperature, lambda);
   }
   template <typename Indexer_t = Real *>
+  PORTABLE_FUNCTION void
+  InternalEnergyFromDensityPressure(const Real rho, const Real P, Real &sie,
+                                    Indexer_t &&lambda = nullptr) const {
+    t_.InternalEnergyFromDensityPressure(rho, P, sie, lambda);
+  }
+  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,
                                  Indexer_t &&lambda = nullptr) const {
diff --git a/singularity-eos/eos/modifiers/scaled_eos.hpp b/singularity-eos/eos/modifiers/scaled_eos.hpp
index c7c2292a3e9..63851ded5fb 100644
--- a/singularity-eos/eos/modifiers/scaled_eos.hpp
+++ b/singularity-eos/eos/modifiers/scaled_eos.hpp
@@ -136,6 +136,12 @@ class ScaledEOS : public EosBase<ScaledEOS<T>> {
     return t_.GruneisenParamFromDensityTemperature(scale_ * rho, temperature, lambda);
   }
   template <typename Indexer_t = Real *>
+  PORTABLE_FUNCTION void
+  InternalEnergyFromDensityPressure(const Real rho, const Real P, Real &sie,
+                                    Indexer_t &&lambda = nullptr) const {
+    t_.InternalEnergyFromDensityPressure(scale_ * rho, P, sie, lambda);
+  }
+  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,
                                  Indexer_t &&lambda = nullptr) const {
diff --git a/singularity-eos/eos/modifiers/shifted_eos.hpp b/singularity-eos/eos/modifiers/shifted_eos.hpp
index ff52b2eac26..cd0bad539b6 100644
--- a/singularity-eos/eos/modifiers/shifted_eos.hpp
+++ b/singularity-eos/eos/modifiers/shifted_eos.hpp
@@ -127,6 +127,14 @@ class ShiftedEOS : public EosBase<ShiftedEOS<T>> {
     return t_.GruneisenParamFromDensityTemperature(rho, temperature, lambda);
   }
   template <typename Indexer_t = Real *>
+  PORTABLE_FUNCTION void
+  InternalEnergyFromDensityPressure(const Real rho, const Real P, Real &sie,
+                                    Indexer_t &&lambda = nullptr) const {
+    sie -= shift_;
+    t_.InternalEnergyFromDensityPressure(rho, P, sie, lambda);
+    sie += shift_;
+  }
+  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,
                                  Indexer_t &&lambda = nullptr) const {
diff --git a/singularity-eos/eos/modifiers/zsplit_eos.hpp b/singularity-eos/eos/modifiers/zsplit_eos.hpp
index ce850acd4e4..ee9f260e400 100644
--- a/singularity-eos/eos/modifiers/zsplit_eos.hpp
+++ b/singularity-eos/eos/modifiers/zsplit_eos.hpp
@@ -169,6 +169,17 @@ class ZSplit : public EosBase<ZSplit<ztype, T>> {
     return t_.GruneisenParamFromDensityInternalEnergy(rho, iscale * sie, lambda);
   }
 
+  template <typename Indexer_t = Real *>
+  PORTABLE_FUNCTION void
+  InternalEnergyFromDensityPressure(const Real rho, const Real P, Real &sie,
+                                    Indexer_t &&lambda = nullptr) const {
+    const Real scale = GetScale_(lambda);
+    const Real iscale = GetInvScale_(lambda);
+    sie *= iscale;
+    t_.InternalEnergyFromDensityPressure(rho, P, sie, lambda);
+    sie *= scale;
+  }
+
   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_gruneisen.cpp b/test/test_eos_gruneisen.cpp
index 376760c64b1..ac1df8f9e6f 100644
--- a/test/test_eos_gruneisen.cpp
+++ b/test/test_eos_gruneisen.cpp
@@ -1,5 +1,5 @@
 //------------------------------------------------------------------------------
-// © 2021-2024. Triad National Security, LLC. All rights reserved.  This
+// © 2021-2025. 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
@@ -208,6 +208,67 @@ SCENARIO("Gruneisen EOS", "[VectorEOS][GruneisenEOS]") {
   }
 }
 
+SCENARIO("Gruneisen energy from pressure and temperature",
+         "[Gruneisen][VectorEOS][SieFromRhoP]") {
+  GIVEN("An ideal gas object") {
+    // Unit conversions
+    constexpr Real cm = 1.;
+    constexpr Real us = 1e-06;
+    constexpr Real Mbcc_per_g = 1e12;
+    // Gruneisen parameters for copper
+    constexpr Real C0 = 0.394 * cm / us;
+    constexpr Real S1 = 1.489;
+    constexpr Real S2 = 0.;
+    constexpr Real S3 = 0.;
+    constexpr Real Gamma0 = 2.02;
+    constexpr Real b = 0.47;
+    constexpr Real rho0 = 8.93;
+    constexpr Real T0 = 298.;
+    constexpr Real P0 = 0.;
+    constexpr Real Cv = 0.383e-05 * Mbcc_per_g;
+    // Create the EOS
+    EOS eos_host = Gruneisen(C0, S1, S2, S3, Gamma0, b, rho0, T0, P0, Cv);
+    EOS eos = eos_host.GetOnDevice();
+
+    WHEN("We compute pressure from energy for several energies") {
+      constexpr int N = 20;
+      Real *rhos = (Real *)PORTABLE_MALLOC(sizeof(Real) * N);
+      Real *sies = (Real *)PORTABLE_MALLOC(sizeof(Real) * N);
+      Real *Ps = (Real *)PORTABLE_MALLOC(sizeof(Real) * N);
+      portableFor(
+          "Set rho, sie, and P", 0, N, PORTABLE_LAMBDA(const int i) {
+            rhos[i] = 0.5 * (eos.MinimumDensity() + eos.MaximumDensity());
+            sies[i] = Cv * (298 + i);
+            Ps[i] = eos.PressureFromDensityInternalEnergy(rhos[i], sies[i]);
+          });
+
+      THEN("We can compute the internal energy from the pressure") {
+        Real *sies_new = (Real *)PORTABLE_MALLOC(sizeof(Real) * N);
+        eos_host.InternalEnergyFromDensityPressure(rhos, Ps, sies_new, N);
+
+        AND_THEN("They agree with the previous pressures") {
+          std::size_t nwrong = 0;
+          portableReduce(
+              "Compute sie diff", 0, N,
+              PORTABLE_LAMBDA(const int i, std::size_t &nw) {
+                nw += !isClose(sies[i], sies_new[i], 1e-12);
+              },
+              nwrong);
+          REQUIRE(nwrong == 0);
+        }
+
+        PORTABLE_FREE(sies_new);
+      }
+
+      PORTABLE_FREE(Ps);
+      PORTABLE_FREE(sies);
+      PORTABLE_FREE(rhos);
+    }
+
+    eos.Finalize();
+  }
+}
+
 SCENARIO("Aluminum Gruneisen EOS", "[GruneisenEOS]") {
   GIVEN("Parameters for an aluminum Gruneisen EOS") {
     // Unit conversions
diff --git a/test/test_eos_ideal.cpp b/test/test_eos_ideal.cpp
index 6685370744c..4b90f524894 100644
--- a/test/test_eos_ideal.cpp
+++ b/test/test_eos_ideal.cpp
@@ -89,6 +89,55 @@ SCENARIO("Ideal gas entropy", "[IdealGas][Entropy][GibbsFreeEnergy]") {
   }
 }
 
+SCENARIO("Ideal gas energy from pressure and temperature",
+         "[IdealGas][VectorEOS][SieFromRhoP]") {
+  GIVEN("An ideal gas object") {
+    constexpr Real Cv = 5.0;
+    constexpr Real gm1 = 0.4;
+    EOS eos_host = IdealGas(gm1, Cv);
+    EOS eos = eos_host.GetOnDevice();
+
+    WHEN("We compute pressure from energy for several energies") {
+      constexpr int N = 20;
+      constexpr Real sie_min = 0;
+      constexpr Real dsie = 0.5;
+      Real *rhos = (Real *)PORTABLE_MALLOC(sizeof(Real) * N);
+      Real *sies = (Real *)PORTABLE_MALLOC(sizeof(Real) * N);
+      Real *Ps = (Real *)PORTABLE_MALLOC(sizeof(Real) * N);
+      portableFor(
+          "Set rho, sie, and P", 0, N, PORTABLE_LAMBDA(const int i) {
+            rhos[i] = 1;
+            sies[i] = sie_min + dsie * i;
+            Ps[i] = eos.PressureFromDensityInternalEnergy(rhos[i], sies[i]);
+          });
+
+      THEN("We can compute the internal energy from the pressure") {
+        Real *sies_new = (Real *)PORTABLE_MALLOC(sizeof(Real) * N);
+        eos_host.InternalEnergyFromDensityPressure(rhos, Ps, sies_new, N);
+
+        AND_THEN("They agree with the previous pressures") {
+          std::size_t nwrong = 0;
+          portableReduce(
+              "Compute sie diff", 0, N,
+              PORTABLE_LAMBDA(const int i, std::size_t &nw) {
+                nw += !isClose(sies[i], sies_new[i], 1e-12);
+              },
+              nwrong);
+          REQUIRE(nwrong == 0);
+        }
+
+        PORTABLE_FREE(sies_new);
+      }
+
+      PORTABLE_FREE(Ps);
+      PORTABLE_FREE(sies);
+      PORTABLE_FREE(rhos);
+    }
+
+    eos.Finalize();
+  }
+}
+
 SCENARIO("Ideal gas mean atomic properties",
          "[IdealGas][MeanAtomicMass][MeanAtomicNumber]") {
   constexpr Real Cv = 5.0;