diff --git a/CHANGELOG.md b/CHANGELOG.md
index 3d938dc22c..0ead0a9bd6 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -3,6 +3,7 @@
 ## Current develop
 
 ### Added (new features/APIs/variables/...)
+- [[PR416]](https://github.com/lanl/singularity-eos/pull/416) Gibbs free energy
 - [[PR361]](https://github.com/lanl/singularity-eos/pull/361) Added tests for PTEsolver and added a strawman kinetic phase transition framework
 - [[PR410]](https://github.com/lanl/singularity-eos/pull/410) Enable serialization and de-serialization
 - [[PR330]](https://github.com/lanl/singularity-eos/pull/330) Piecewise grids for Spiner EOS.
diff --git a/doc/sphinx/src/using-eos.rst b/doc/sphinx/src/using-eos.rst
index fe570ad625..c2fcc8c3e3 100644
--- a/doc/sphinx/src/using-eos.rst
+++ b/doc/sphinx/src/using-eos.rst
@@ -981,6 +981,50 @@ 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 EntropyFromDensityTemperature(const Real rho, const Real temperature,
+                                      Indexer_t &&lambda = nullptr) const;
+
+returns the entropy as a function of density in :math:`g/cm^3` and
+temperature in Kelvin.
+
+.. code-block:: cpp
+
+   template <typename Indexer_t = Real*>
+   Real EntropyFromDensityInternalEnergy(const Real rho, const Real sie,
+                                         Indexer_t &&lambda = nullptr) const;
+
+returns the entropy as a function of density in :math:`g/cm^3` and
+specific internal energy in :math:`erg/g`.
+
+.. code-block:: cpp
+
+   template <typename Indexer_t = Real*>
+   Real GibssFreeEnergyFromDensityTemperature(const Real rho, const Real temperature,
+                                              Indexer_t &&lambda = nullptr) const;
+
+returns the Gibbs free energy as a function of density in :math:`g/cm^3` and
+temperature in Kelvin.
+
+.. code-block:: cpp
+
+   template <typename Indexer_t = Real*>
+   Real GibbsFreeEnergyFromDensityInternalEnergy(const Real rho, const Real sie,
+                                                 Indexer_t &&lambda = nullptr) const;
+
+
+returns the Gibbs free energy as a function of density in :math:`g/cm^3` and
+specific internal energy in :math:`erg/g`.
+
+.. warning::
+
+  Not all equations of state provide entropy and Gibbs free
+  energy. These are coupled, however, so if one is provided, the other
+  will be too. If you call an entropy for a model that does not
+  provide it, ``singularity-eos`` will return an error.
+
 The function
 
 .. code-block:: cpp
diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp
index 0d627e75fe..8b305dba54 100644
--- a/singularity-eos/eos/eos_base.hpp
+++ b/singularity-eos/eos/eos_base.hpp
@@ -79,7 +79,9 @@ char *StrCat(char *destination, const char *source) {
   using EosBase<EOSDERIVED>::GruneisenParamFromDensityInternalEnergy;                    \
   using EosBase<EOSDERIVED>::FillEos;                                                    \
   using EosBase<EOSDERIVED>::EntropyFromDensityTemperature;                              \
-  using EosBase<EOSDERIVED>::EntropyFromDensityInternalEnergy;
+  using EosBase<EOSDERIVED>::EntropyFromDensityInternalEnergy;                           \
+  using EosBase<EOSDERIVED>::GibbsFreeEnergyFromDensityTemperature;                      \
+  using EosBase<EOSDERIVED>::GibbsFreeEnergyFromDensityInternalEnergy;
 
 // This macro adds several methods that most modifiers will
 // want. Not ALL modifiers will want these methods as written here,
@@ -184,6 +186,28 @@ class EosBase {
                                     std::forward<Args>(args)...);
   }
 
+  // Scalar member functions that get shared
+  template <typename Indexer_t = Real *>
+  PORTABLE_INLINE_FUNCTION Real GibbsFreeEnergyFromDensityTemperature(
+      const Real rho, const Real T,
+      Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
+    const CRTP copy = *(static_cast<CRTP const *>(this));
+    Real sie = copy.InternalEnergyFromDensityTemperature(rho, T, lambda);
+    Real P = copy.PressureFromDensityTemperature(rho, T, lambda);
+    Real S = copy.EntropyFromDensityTemperature(rho, T, lambda);
+    return sie + (P / rho) - T * S;
+  }
+  template <typename Indexer_t = Real *>
+  PORTABLE_INLINE_FUNCTION Real GibbsFreeEnergyFromDensityInternalEnergy(
+      const Real rho, const Real sie,
+      Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
+    const CRTP copy = *(static_cast<CRTP const *>(this));
+    Real T = copy.TemperatureFromDensityInternalEnergy(rho, sie, lambda);
+    Real P = copy.PressureFromDensityTemperature(rho, T, lambda);
+    Real S = copy.EntropyFromDensityTemperature(rho, T, lambda);
+    return sie + (P / rho) - T * S;
+  }
+
   // Vector member functions
   template <typename RealIndexer, typename ConstRealIndexer, typename LambdaIndexer>
   inline void
@@ -618,6 +642,73 @@ class EosBase {
     GruneisenParamFromDensityInternalEnergy(rhos, sies, gm1s, num,
                                             std::forward<LambdaIndexer>(lambdas));
   }
+  template <typename RealIndexer, typename ConstRealIndexer, typename LambdaIndexer>
+  inline void GibbsFreeEnergyFromDensityTemperature(ConstRealIndexer &&rhos,
+                                                    ConstRealIndexer &&Ts,
+                                                    RealIndexer &&Gs, 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) {
+          Gs[i] = copy.GibbsFreeEnergyFromDensityTemperature(rhos[i], Ts[i], lambdas[i]);
+        });
+  }
+  template <typename RealIndexer, typename ConstRealIndexer, typename LambdaIndexer,
+            typename = std::enable_if_t<!is_raw_pointer<RealIndexer, Real>::value>>
+  inline void
+  GibbsFreeEnergyFromDensityTemperature(ConstRealIndexer &&rhos, ConstRealIndexer &&Ts,
+                                        RealIndexer &&Gs, Real * /*scratch*/,
+                                        const int num, LambdaIndexer &&lambdas) const {
+    GibbsFreeEnergyFromDensityTemperature(
+        std::forward<ConstRealIndexer>(rhos), std::forward<ConstRealIndexer>(Ts),
+        std::forward<RealIndexer>(Gs), num, std::forward<LambdaIndexer>(lambdas));
+  }
+  template <typename LambdaIndexer>
+  inline void GibbsFreeEnergyFromDensityTemperature(const Real *rhos, const Real *Ts,
+                                                    Real *Gs, Real * /*scratch*/,
+                                                    const int num,
+                                                    LambdaIndexer &&lambdas,
+                                                    Transform && = Transform()) const {
+    GibbsFreeEnergyFromDensityTemperature(rhos, Ts, Gs, num,
+                                          std::forward<LambdaIndexer>(lambdas));
+  }
+  template <typename RealIndexer, typename ConstRealIndexer, typename LambdaIndexer>
+  inline void GibbsFreeEnergyFromDensityInternalEnergy(ConstRealIndexer &&rhos,
+                                                       ConstRealIndexer &&sies,
+                                                       RealIndexer &&Gs, 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) {
+          Gs[i] =
+              copy.GibbsFreeEnergyFromDensityInternalEnergy(rhos[i], sies[i], lambdas[i]);
+        });
+  }
+  template <typename RealIndexer, typename ConstRealIndexer, typename LambdaIndexer,
+            typename = std::enable_if_t<!is_raw_pointer<RealIndexer, Real>::value>>
+  inline void GibbsFreeEnergyFromDensityInternalEnergy(ConstRealIndexer &&rhos,
+                                                       ConstRealIndexer &&sies,
+                                                       RealIndexer &&Gs,
+                                                       Real * /*scratch*/, const int num,
+                                                       LambdaIndexer &&lambdas) const {
+    GibbsFreeEnergyFromDensityInternalEnergy(
+        std::forward<ConstRealIndexer>(rhos), std::forward<ConstRealIndexer>(sies),
+        std::forward<RealIndexer>(Gs), num, std::forward<LambdaIndexer>(lambdas));
+  }
+  template <typename LambdaIndexer>
+  inline void GibbsFreeEnergyFromDensityInternalEnergy(const Real *rhos, const Real *sies,
+                                                       Real *Gs, Real * /*scratch*/,
+                                                       const int num,
+                                                       LambdaIndexer &&lambdas,
+                                                       Transform && = Transform()) const {
+    GibbsFreeEnergyFromDensityInternalEnergy(rhos, sies, Gs, num,
+                                             std::forward<LambdaIndexer>(lambdas));
+  }
+
   template <typename RealIndexer, typename LambdaIndexer>
   inline void FillEos(RealIndexer &&rhos, RealIndexer &&temps, RealIndexer &&energies,
                       RealIndexer &&presses, RealIndexer &&cvs, RealIndexer &&bmods,
@@ -632,6 +723,7 @@ class EosBase {
                        output, lambdas[i]);
         });
   }
+
   // Report minimum values of density and temperature
   PORTABLE_FORCEINLINE_FUNCTION
   Real MinimumDensity() const { return 0; }
diff --git a/singularity-eos/eos/eos_eospac.hpp b/singularity-eos/eos/eos_eospac.hpp
index e65c71d95c..00450a2a9c 100644
--- a/singularity-eos/eos/eos_eospac.hpp
+++ b/singularity-eos/eos/eos_eospac.hpp
@@ -472,6 +472,7 @@ class EOSPAC : public EosBase<EOSPAC> {
                                                              lambdas);
   }
 
+  // TODO(JMM): Add performant entropy and Gibbs Free Energy
   using EosBase<EOSPAC>::FillEos;
   using EosBase<EOSPAC>::EntropyIsNotEnabled;
 
diff --git a/singularity-eos/eos/eos_variant.hpp b/singularity-eos/eos/eos_variant.hpp
index 507e888423..938bb1da40 100644
--- a/singularity-eos/eos/eos_variant.hpp
+++ b/singularity-eos/eos/eos_variant.hpp
@@ -178,6 +178,28 @@ class Variant {
         eos_);
   }
 
+  template <typename Indexer_t = Real *>
+  PORTABLE_INLINE_FUNCTION Real GibbsFreeEnergyFromDensityTemperature(
+      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.GibbsFreeEnergyFromDensityTemperature(rho, T, lambda);
+        },
+        eos_);
+  }
+
+  template <typename Indexer_t = Real *>
+  PORTABLE_INLINE_FUNCTION Real GibbsFreeEnergyFromDensityInternalEnergy(
+      const Real rho, const Real sie,
+      Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
+    return mpark::visit(
+        [&rho, &sie, &lambda](const auto &eos) {
+          return eos.GibbsFreeEnergyFromDensityInternalEnergy(rho, sie, lambda);
+        },
+        eos_);
+  }
+
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityTemperature(
       const Real rho, const Real temperature,
@@ -672,6 +694,113 @@ class Variant {
         eos_);
   }
 
+  template <typename RealIndexer, typename ConstRealIndexer>
+  inline void GibbsFreeEnergyFromDensityTemperature(ConstRealIndexer &&rhos,
+                                                    ConstRealIndexer &&temperatures,
+                                                    RealIndexer &&Gs,
+                                                    const int num) const {
+    NullIndexer lambdas{}; // Returns null pointer for every index
+    return GibbsFreeEnergyFromDensityTemperature(
+        std::forward<ConstRealIndexer>(rhos),
+        std::forward<ConstRealIndexer>(temperatures), std::forward<RealIndexer>(Gs), num,
+        lambdas);
+  }
+
+  template <typename RealIndexer, typename ConstRealIndexer, typename LambdaIndexer>
+  inline void GibbsFreeEnergyFromDensityTemperature(ConstRealIndexer &&rhos,
+                                                    ConstRealIndexer &&temperatures,
+                                                    RealIndexer &&Gs, const int num,
+                                                    LambdaIndexer &&lambdas) const {
+    return mpark::visit(
+        [&rhos, &temperatures, &Gs, &num, &lambdas](const auto &eos) {
+          return eos.GibbsFreeEnergyFromDensityTemperature(
+              std::forward<ConstRealIndexer>(rhos),
+              std::forward<ConstRealIndexer>(temperatures), std::forward<RealIndexer>(Gs),
+              num, std::forward<LambdaIndexer>(lambdas));
+        },
+        eos_);
+  }
+
+  template <typename RealIndexer, typename ConstRealIndexer>
+  inline void GibbsFreeEnergyFromDensityTemperature(ConstRealIndexer &&rhos,
+                                                    ConstRealIndexer &&temperatures,
+                                                    RealIndexer &&Gs, Real *scratch,
+                                                    const int num) const {
+    NullIndexer lambdas{}; // Returns null pointer for every index
+    return GibbsFreeEnergyFromDensityTemperature(
+        std::forward<ConstRealIndexer>(rhos),
+        std::forward<ConstRealIndexer>(temperatures), std::forward<RealIndexer>(Gs),
+        scratch, num, lambdas);
+  }
+
+  template <typename RealIndexer, typename ConstRealIndexer, typename LambdaIndexer>
+  inline void GibbsFreeEnergyFromDensityTemperature(ConstRealIndexer &&rhos,
+                                                    ConstRealIndexer &&temperatures,
+                                                    RealIndexer &&Gs, Real *scratch,
+                                                    const int num,
+                                                    LambdaIndexer &&lambdas) const {
+    return mpark::visit(
+        [&rhos, &temperatures, &Gs, &scratch, &num, &lambdas](const auto &eos) {
+          return eos.GibbsFreeEnergyFromDensityTemperature(
+              std::forward<ConstRealIndexer>(rhos),
+              std::forward<ConstRealIndexer>(temperatures), std::forward<RealIndexer>(Gs),
+              scratch, num, std::forward<LambdaIndexer>(lambdas));
+        },
+        eos_);
+  }
+
+  template <typename RealIndexer, typename ConstRealIndexer>
+  inline void GibbsFreeEnergyFromDensityInternalEnergy(ConstRealIndexer &&rhos,
+                                                       ConstRealIndexer &&sies,
+                                                       RealIndexer &&Gs,
+                                                       const int num) const {
+    NullIndexer lambdas{}; // Returns null pointer for every index
+    return GibbsFreeEnergyFromDensityInternalEnergy(
+        std::forward<ConstRealIndexer>(rhos), std::forward<ConstRealIndexer>(sies),
+        std::forward<RealIndexer>(Gs), num, lambdas);
+  }
+
+  template <typename RealIndexer, typename ConstRealIndexer, typename LambdaIndexer>
+  inline void GibbsFreeEnergyFromDensityInternalEnergy(ConstRealIndexer &&rhos,
+                                                       ConstRealIndexer &&sies,
+                                                       RealIndexer &&Gs, const int num,
+                                                       LambdaIndexer &&lambdas) const {
+    return mpark::visit(
+        [&rhos, &sies, &Gs, &num, &lambdas](const auto &eos) {
+          return eos.GibbsFreeEnergyFromDensityInternalEnergy(
+              std::forward<ConstRealIndexer>(rhos), std::forward<ConstRealIndexer>(sies),
+              std::forward<RealIndexer>(Gs), num, std::forward<LambdaIndexer>(lambdas));
+        },
+        eos_);
+  }
+
+  template <typename RealIndexer, typename ConstRealIndexer, typename LambdaIndexer>
+  inline void GibbsFreeEnergyFromDensityInternalEnergy(ConstRealIndexer &&rhos,
+                                                       ConstRealIndexer &&sies,
+                                                       RealIndexer &&Gs, Real *scratch,
+                                                       const int num) const {
+    NullIndexer lambdas{}; // Returns null pointer for every index
+    return GibbsFreeEnergyFromDensityInternalEnergy(
+        std::forward<ConstRealIndexer>(rhos), std::forward<ConstRealIndexer>(sies),
+        std::forward<RealIndexer>(Gs), scratch, num, lambdas);
+  }
+
+  template <typename RealIndexer, typename ConstRealIndexer, typename LambdaIndexer>
+  inline void GibbsFreeEnergyFromDensityInternalEnergy(ConstRealIndexer &&rhos,
+                                                       ConstRealIndexer &&sies,
+                                                       RealIndexer &&Gs, Real *scratch,
+                                                       const int num,
+                                                       LambdaIndexer &&lambdas) const {
+    return mpark::visit(
+        [&rhos, &sies, &Gs, &scratch, &num, &lambdas](const auto &eos) {
+          return eos.GibbsFreeEnergyFromDensityInternalEnergy(
+              std::forward<ConstRealIndexer>(rhos), std::forward<ConstRealIndexer>(sies),
+              std::forward<RealIndexer>(Gs), scratch, num,
+              std::forward<LambdaIndexer>(lambdas));
+        },
+        eos_);
+  }
+
   template <typename RealIndexer, typename ConstRealIndexer>
   inline void SpecificHeatFromDensityTemperature(ConstRealIndexer &&rhos,
                                                  ConstRealIndexer &&temperatures,
diff --git a/test/test_eos_ideal.cpp b/test/test_eos_ideal.cpp
index d33c883ea0..b97f3de4de 100644
--- a/test/test_eos_ideal.cpp
+++ b/test/test_eos_ideal.cpp
@@ -32,7 +32,7 @@
 using singularity::IdealGas;
 using EOS = singularity::Variant<IdealGas>;
 
-SCENARIO("Ideal gas entropy", "[IdealGas][Entropy]") {
+SCENARIO("Ideal gas entropy", "[IdealGas][Entropy][GibbsFreeEnergy]") {
   GIVEN("Parameters for an ideal gas with entropy reference states") {
     // Create ideal gas EOS ojbect
     constexpr Real Cv = 5.0;
@@ -54,6 +54,16 @@ SCENARIO("Ideal gas entropy", "[IdealGas][Entropy]") {
         auto entropy = host_eos.EntropyFromDensityTemperature(rho, T);
         INFO("Entropy: " << entropy << "  True entropy: " << entropy_true);
         CHECK(isClose(entropy, entropy_true, 1e-12));
+
+        AND_THEN("The free energy agrees") {
+          const Real sie = host_eos.InternalEnergyFromDensityTemperature(rho, T);
+          const Real P = host_eos.PressureFromDensityTemperature(rho, T);
+          const Real G_true = sie + (P / rho) - T * entropy;
+          const Real GT = host_eos.GibbsFreeEnergyFromDensityTemperature(rho, T);
+          CHECK(isClose(GT, G_true, 1e-12));
+          const Real Gsie = host_eos.GibbsFreeEnergyFromDensityInternalEnergy(rho, sie);
+          CHECK(isClose(Gsie, G_true, 1e-12));
+        }
       }
     }
     GIVEN("A state at the reference density and a temperature whose square is the "