From 332682a5f5275e4796a183c468713302be82621b Mon Sep 17 00:00:00 2001
From: Cory Kinney <cory.kinney@ucf.edu>
Date: Mon, 2 Jan 2023 14:50:58 -0500
Subject: [PATCH 01/17] Add initial SoaveRedlichKwong files

---
 include/cantera/thermo/SoaveRedlichKwong.h | 35 ++++++++++++++++++++++
 src/thermo/SoaveRedlichKwong.cpp           | 20 +++++++++++++
 2 files changed, 55 insertions(+)
 create mode 100644 include/cantera/thermo/SoaveRedlichKwong.h
 create mode 100644 src/thermo/SoaveRedlichKwong.cpp

diff --git a/include/cantera/thermo/SoaveRedlichKwong.h b/include/cantera/thermo/SoaveRedlichKwong.h
new file mode 100644
index 0000000000..270d30642e
--- /dev/null
+++ b/include/cantera/thermo/SoaveRedlichKwong.h
@@ -0,0 +1,35 @@
+//! @file SoaveRedlichKwong.h
+
+// This file is part of Cantera. See License.txt in the top-level directory or
+// at https://cantera.org/license.txt for license and copyright information.
+
+#ifndef CT_SOAVEREDLICHKWONG_H
+#define CT_SOAVEREDLICHKWONG_H
+
+#include "MixtureFugacityTP.h"
+#include "cantera/base/Array.h"
+
+namespace Cantera
+{
+/**
+ * Implementation of a multi-species Soave-Redlich-Kwong equation of state
+ *
+ * @ingroup thermoprops
+ */
+class SoaveRedlichKwong : public MixtureFugacityTP
+{
+public:
+    explicit SoaveRedlichKwong(const std::string& infile="",
+                               const std::string& id="");
+
+    virtual std::string type() const {
+        return "Soave-Redlich-Kwong";
+    }
+private:
+    static const double omega_a;
+    static const double omega_b;
+    static const double omega_vc;
+}
+}
+
+#endif
diff --git a/src/thermo/SoaveRedlichKwong.cpp b/src/thermo/SoaveRedlichKwong.cpp
new file mode 100644
index 0000000000..77042c0345
--- /dev/null
+++ b/src/thermo/SoaveRedlichKwong.cpp
@@ -0,0 +1,20 @@
+//! @file SoaveRedlichKwong.cpp
+
+// This file is part of Cantera. See License.txt in the top-level directory or
+// at https://cantera.org/license.txt for license and copyright information.
+
+#include "cantera/thermo/SoaveRedlichKwong.h"
+#include "cantera/thermo/ThermoFactory.h"
+#include "cantera/thermo/Species.h"
+#include "cantera/base/stringUtils.h"
+
+#include <boost/math/tools/roots.hpp>
+
+namespace bmt = boost::math::tools;
+
+namespace Cantera
+{
+const double SoaveRedlichKwong::omega_a = 4.27480233540E-01;
+const double SoaveRedlichKwong::omega_b = 8.66403499650E-02;
+const double SoaveRedlichKwong::omega_vc = 3.33333333333333E-01;
+}

From fd27820838966ca67f1994ac3893d4187bf5611b Mon Sep 17 00:00:00 2001
From: Cory Kinney <cory.kinney@ucf.edu>
Date: Wed, 4 Jan 2023 15:57:50 -0500
Subject: [PATCH 02/17] Add solveCubic

---
 include/cantera/thermo/SoaveRedlichKwong.h |  8 +++++++-
 src/thermo/SoaveRedlichKwong.cpp           | 19 +++++++++++++++++++
 2 files changed, 26 insertions(+), 1 deletion(-)

diff --git a/include/cantera/thermo/SoaveRedlichKwong.h b/include/cantera/thermo/SoaveRedlichKwong.h
index 270d30642e..19e59ce247 100644
--- a/include/cantera/thermo/SoaveRedlichKwong.h
+++ b/include/cantera/thermo/SoaveRedlichKwong.h
@@ -25,11 +25,17 @@ class SoaveRedlichKwong : public MixtureFugacityTP
     virtual std::string type() const {
         return "Soave-Redlich-Kwong";
     }
+
+public:
+    int SoaveRedlichKwong::solveCubic(double T, double pres, double a, double b, double aAlpha,
+                                  double Vroot[3]) const;
+
 private:
     static const double omega_a;
     static const double omega_b;
     static const double omega_vc;
-}
+
+};
 }
 
 #endif
diff --git a/src/thermo/SoaveRedlichKwong.cpp b/src/thermo/SoaveRedlichKwong.cpp
index 77042c0345..b66e7c56b9 100644
--- a/src/thermo/SoaveRedlichKwong.cpp
+++ b/src/thermo/SoaveRedlichKwong.cpp
@@ -17,4 +17,23 @@ namespace Cantera
 const double SoaveRedlichKwong::omega_a = 4.27480233540E-01;
 const double SoaveRedlichKwong::omega_b = 8.66403499650E-02;
 const double SoaveRedlichKwong::omega_vc = 3.33333333333333E-01;
+
+
+
+int SoaveRedlichKwong::solveCubic(double T, double pres, double a, double b, double aAlpha,
+                                  double Vroot[3]) const
+{
+    double an = 1.0;
+    double bn = - GasConstant * T / pres;
+    double cn = (aAlpha - b * GasConstant * T) / pres - b * b;
+    double dn = aAlpha * b / pres;
+
+    double tc = a * omega_b / (b * omega_a * GasConstant);
+    double pc = omega_b * R * tc / b;
+    double vc = omega_vc * GasConstant * tc / pc
+
+    return MixtureFugacityTP::solveCubic(T, pres, a, b, aAlpha, Vroot,
+                                         an, bn, cn, dn, tc, vc);
+}
+
 }

From 4363c810596a1a0440b416a41b5502522e6b9eef Mon Sep 17 00:00:00 2001
From: Cory Kinney <cory.kinney@ucf.edu>
Date: Thu, 5 Jan 2023 13:08:42 -0500
Subject: [PATCH 03/17] Add derivatives

---
 include/cantera/thermo/SoaveRedlichKwong.h | 62 +++++++++++++++++-
 src/thermo/SoaveRedlichKwong.cpp           | 76 +++++++++++++++++++++-
 2 files changed, 135 insertions(+), 3 deletions(-)

diff --git a/include/cantera/thermo/SoaveRedlichKwong.h b/include/cantera/thermo/SoaveRedlichKwong.h
index 19e59ce247..4b5668073d 100644
--- a/include/cantera/thermo/SoaveRedlichKwong.h
+++ b/include/cantera/thermo/SoaveRedlichKwong.h
@@ -26,9 +26,67 @@ class SoaveRedlichKwong : public MixtureFugacityTP
         return "Soave-Redlich-Kwong";
     }
 
+protected:
+
+    virtual double dpdVCalc(double T, double molarVol, double& presCalc) const;
+
+    // Special functions not inherited from MixtureFugacityTP
+
+    //! Calculate temperature derivative \f$d(a \alpha)/dT\f$
+    /*!
+     *  These are stored internally.
+     */
+    double daAlpha_dT() const;
+
+    //! Calculate second derivative \f$d^2(a \alpha)/dT^2\f$
+    /*!
+     *  These are stored internally.
+     */
+    double d2aAlpha_dT2() const;
+
 public:
-    int SoaveRedlichKwong::solveCubic(double T, double pres, double a, double b, double aAlpha,
-                                  double Vroot[3]) const;
+    //! Calculate \f$dp/dV\f$ and \f$dp/dT\f$ at the current conditions
+    /*!
+     *  These are stored internally.
+     */
+    void calculatePressureDerivatives() const;
+
+    //! Prepare variables and call the function to solve the cubic equation of state
+    int solveCubic(double T, double pres, double a, double b, double aAlpha,
+                   double Vroot[3]) const;
+
+protected:
+    //! Value of \f$b\f$ in the equation of state
+    /*!
+     *  `m_b` is a function of the mole fractions and species-specific b values.
+     */
+    double m_b;
+
+    //! Value of \f$a\f$ in the equation of state
+    /*!
+     *  `m_a` depends only on the mole fractions.
+     */
+    double m_a;
+
+    //! Value of \f$\alpha\f$ in the equation of state
+    /*!
+     *  `m_aAlpha_mix` is a function of the temperature and the mole fractions.
+     */
+    double m_aAlpha_mix;
+
+    //! The derivative of the pressure with respect to the volume
+    /*!
+     * Calculated at the current conditions. temperature and mole number kept
+     * constant
+     */
+    mutable double m_dpdV;
+
+    //! The derivative of the pressure with respect to the temperature
+    /*!
+     *  Calculated at the current conditions. Total volume and mole number kept
+     *  constant
+     */
+    mutable double m_dpdT;
 
 private:
     static const double omega_a;
diff --git a/src/thermo/SoaveRedlichKwong.cpp b/src/thermo/SoaveRedlichKwong.cpp
index b66e7c56b9..3466000c9b 100644
--- a/src/thermo/SoaveRedlichKwong.cpp
+++ b/src/thermo/SoaveRedlichKwong.cpp
@@ -20,6 +20,80 @@ const double SoaveRedlichKwong::omega_vc = 3.33333333333333E-01;
 
 
 
+double SoaveRedlichKwong::dpdVCalc(double T, double molarVol, double& presCalc) const
+{
+    double denom = molarVol * (molarVol + m_b);
+    double vmb = molarVol - m_b;
+    return -GasConstant * T / (vmb * vmb) + m_aAlpha_mix * (2 * molarVol + m_b) / (denom * denom);
+}
+
+void SoaveRedlichKwong::calculatePressureDerivatives() const
+{
+    double T = temperature();
+    double mv = molarVolume();
+    double pres;
+
+    m_dpdV = dpdVCalc(T, mv, pres);
+    m_dpdT = GasConstant / (mv - m_b) - daAlpha_dT() / (mv * (mv + m_b));
+}
+
+double SoaveRedlichKwong::daAlpha_dT() const
+{
+    // Same as PengRobinson::daAlpha_dT()
+    double daAlphadT = 0.0, k, Tc, sqtTr, coeff1, coeff2;
+    for (size_t i = 0; i < m_kk; i++) {
+        // Calculate first derivative of alpha for individual species
+        Tc = speciesCritTemperature(m_a_coeffs(i,i), m_b_coeffs[i]);
+        sqtTr = sqrt(temperature() / Tc);
+        coeff1 = 1 / (Tc*sqtTr);
+        coeff2 = sqtTr - 1;
+        k = m_kappa[i];
+        m_dalphadT[i] = coeff1 * (k*k*coeff2 - k);
+    }
+    // Calculate mixture derivative
+    for (size_t i = 0; i < m_kk; i++) {
+        for (size_t j = 0; j < m_kk; j++) {
+            daAlphadT += moleFractions_[i] * moleFractions_[j] * 0.5
+                         * m_aAlpha_binary(i, j)
+                         * (m_dalphadT[i] / m_alpha[i] + m_dalphadT[j] / m_alpha[j]);
+        }
+    }
+    return daAlphadT;
+}
+
+double SoaveRedlichKwong::d2aAlpha_dT2() const
+{
+    // Same as PengRobinson::d2aAlpha_dT2()
+    for (size_t i = 0; i < m_kk; i++) {
+        double Tcrit_i = speciesCritTemperature(m_a_coeffs(i, i), m_b_coeffs[i]);
+        double sqt_Tr = sqrt(temperature() / Tcrit_i);
+        double coeff1 = 1 / (Tcrit_i*sqt_Tr);
+        double coeff2 = sqt_Tr - 1;
+        // Calculate first and second derivatives of alpha for individual species
+        double k = m_kappa[i];
+        m_dalphadT[i] = coeff1 * (k*k*coeff2 - k);
+        m_d2alphadT2[i] = (k*k + k) * coeff1 / (2*sqt_Tr*sqt_Tr*Tcrit_i);
+    }
+
+    // Calculate mixture derivative
+    double d2aAlphadT2 = 0.0;
+    for (size_t i = 0; i < m_kk; i++) {
+        double alphai = m_alpha[i];
+        for (size_t j = 0; j < m_kk; j++) {
+            double alphaj = m_alpha[j];
+            double alphaij = alphai * alphaj;
+            double term1 = m_d2alphadT2[i] / alphai + m_d2alphadT2[j] / alphaj;
+            double term2 = 2 * m_dalphadT[i] * m_dalphadT[j] / alphaij;
+            double term3 = m_dalphadT[i] / alphai + m_dalphadT[j] / alphaj;
+            d2aAlphadT2 += 0.5 * moleFractions_[i] * moleFractions_[j]
+                           * m_aAlpha_binary(i, j)
+                           * (term1 + term2 - 0.5 * term3 * term3);
+        }
+    }
+    return d2aAlphadT2;
+}
+
+
 int SoaveRedlichKwong::solveCubic(double T, double pres, double a, double b, double aAlpha,
                                   double Vroot[3]) const
 {
@@ -30,7 +104,7 @@ int SoaveRedlichKwong::solveCubic(double T, double pres, double a, double b, dou
 
     double tc = a * omega_b / (b * omega_a * GasConstant);
     double pc = omega_b * R * tc / b;
-    double vc = omega_vc * GasConstant * tc / pc
+    double vc = omega_vc * GasConstant * tc / pc;
 
     return MixtureFugacityTP::solveCubic(T, pres, a, b, aAlpha, Vroot,
                                          an, bn, cn, dn, tc, vc);

From 2f7d216a17090c9b45290d0fecbe7b48948f6ff8 Mon Sep 17 00:00:00 2001
From: Cory Kinney <cory.kinney@ucf.edu>
Date: Thu, 5 Jan 2023 15:30:35 -0500
Subject: [PATCH 04/17] Add initialization methods

---
 include/cantera/thermo/SoaveRedlichKwong.h |  30 ++++
 src/thermo/SoaveRedlichKwong.cpp           | 155 +++++++++++++++++++++
 2 files changed, 185 insertions(+)

diff --git a/include/cantera/thermo/SoaveRedlichKwong.h b/include/cantera/thermo/SoaveRedlichKwong.h
index 4b5668073d..7ff6999cbf 100644
--- a/include/cantera/thermo/SoaveRedlichKwong.h
+++ b/include/cantera/thermo/SoaveRedlichKwong.h
@@ -26,6 +26,31 @@ class SoaveRedlichKwong : public MixtureFugacityTP
         return "Soave-Redlich-Kwong";
     }
 
+
+    //! Calculate species-specific critical temperature
+    /*!
+     *  The temperature dependent parameter in SRK EoS is calculated as
+     *       \f[ T_{crit} = (0.0778 a)/(0.4572 b R) \f] TODO
+     *  Units: Kelvin
+     *
+     * @param a    species-specific coefficients used in P-R EoS
+     * @param b    species-specific coefficients used in P-R EoS
+     */
+    double speciesCritTemperature(double a, double b) const;
+
+    //! @name Initialization Methods - For Internal use
+    //!
+    //! The following methods are used in the process of constructing
+    //! the phase and setting its parameters from a specification in an
+    //! input file. They are not normally used in application programs.
+    //! To see how they are used, see importPhase().
+    //! @{
+
+    virtual bool addSpecies(shared_ptr<Species> spec);
+    virtual void initThermo();
+    virtual void getSpeciesParameters(const std::string& name,
+                                      AnyMap& speciesNode) const;
+
 protected:
 
     virtual double dpdVCalc(double T, double molarVol, double& presCalc) const;
@@ -88,6 +113,11 @@ class SoaveRedlichKwong : public MixtureFugacityTP
      */
     mutable double m_dpdT;
 
+
+    enum class CoeffSource { EoS, CritProps, Database };
+    //! For each species, specifies the source of the a, b, and omega coefficients
+    std::vector<CoeffSource> m_coeffSource;
+
 private:
     static const double omega_a;
     static const double omega_b;
diff --git a/src/thermo/SoaveRedlichKwong.cpp b/src/thermo/SoaveRedlichKwong.cpp
index 3466000c9b..ae1feca931 100644
--- a/src/thermo/SoaveRedlichKwong.cpp
+++ b/src/thermo/SoaveRedlichKwong.cpp
@@ -19,6 +19,161 @@ const double SoaveRedlichKwong::omega_b = 8.66403499650E-02;
 const double SoaveRedlichKwong::omega_vc = 3.33333333333333E-01;
 
 
+double SoaveRedlichKwong::speciesCritTemperature(double a, double b) const
+{
+    if (b <= 0.0) {
+        return 1000000.;
+    } else if (a <= 0.0) {
+        return 0.0;
+    } else {
+        return a * omega_b / (b * omega_a * GasConstant);
+    }
+}
+
+
+bool SoaveRedlichKwong::addSpecies(shared_ptr<Species> spec)
+{
+    bool added = MixtureFugacityTP::addSpecies(spec);
+    if (added) {
+        m_a_coeffs.resize(m_kk, m_kk, 0.0);
+        m_b_coeffs.push_back(0.0);
+        m_aAlpha_binary.resize(m_kk, m_kk, 0.0);
+        m_kappa.push_back(0.0);
+        m_acentric.push_back(0.0);
+        m_alpha.push_back(0.0);
+        m_dalphadT.push_back(0.0);
+        m_d2alphadT2.push_back(0.0);
+        m_pp.push_back(0.0);
+        m_partialMolarVolumes.push_back(0.0);
+        m_dpdni.push_back(0.0);
+        m_coeffSource.push_back(CoeffSource::EoS);
+    }
+    return added;
+}
+
+
+void SoaveRedlichKwong::initThermo()
+{
+    // Contents of 'critical-properties.yaml', loaded later if needed
+    AnyMap critPropsDb;
+    std::unordered_map<std::string, AnyMap*> dbSpecies;
+
+    for (auto& item : m_species) {
+        auto& data = item.second->input;
+        size_t k = speciesIndex(item.first);
+        if (m_a_coeffs(k, k) != 0.0) {
+            continue;
+        }
+        bool foundCoeffs = false;
+        if (data.hasKey("equation-of-state") &&
+            data["equation-of-state"].hasMapWhere("model", "Soave-Redlich-Kwong"))
+        {
+            // Read a and b coefficients and acentric factor w_ac from species input
+            // information, specified in a YAML input file.
+            auto eos = data["equation-of-state"].getMapWhere(
+                "model", "Soave-Redlich-Kwong");
+            if (eos.hasKey("a") && eos.hasKey("b") && eos.hasKey("acentric-factor")) {
+                double a0 = eos.convert("a", "Pa*m^6/kmol^2");
+                double b = eos.convert("b", "m^3/kmol");
+                // unitless acentric factor:
+                double w = eos["acentric-factor"].asDouble();
+                setSpeciesCoeffs(item.first, a0, b, w);
+                foundCoeffs = true;
+            }
+
+            if (eos.hasKey("binary-a")) {
+                AnyMap& binary_a = eos["binary-a"].as<AnyMap>();
+                const UnitSystem& units = binary_a.units();
+                for (auto& item2 : binary_a) {
+                    double a0 = units.convert(item2.second, "Pa*m^6/kmol^2");
+                    setBinaryCoeffs(item.first, item2.first, a0);
+                }
+            }
+            if (foundCoeffs) {
+                m_coeffSource[k] = CoeffSource::EoS;
+                continue;
+            }
+        }
+
+        // Coefficients have not been populated from model-specific input
+        double Tc = NAN, Pc = NAN, omega_ac = NAN;
+        if (data.hasKey("critical-parameters")) {
+            // Use critical state information stored in the species entry to
+            // calculate a, b, and the acentric factor.
+            auto& critProps = data["critical-parameters"].as<AnyMap>();
+            Tc = critProps.convert("critical-temperature", "K");
+            Pc = critProps.convert("critical-pressure", "Pa");
+            omega_ac = critProps["acentric-factor"].asDouble();
+            m_coeffSource[k] = CoeffSource::CritProps;
+        } else {
+            // Search 'crit-properties.yaml' to find Tc and Pc. Load data if needed.
+            if (critPropsDb.empty()) {
+                critPropsDb = AnyMap::fromYamlFile("critical-properties.yaml");
+                dbSpecies = critPropsDb["species"].asMap("name");
+            }
+
+            // All names in critical-properties.yaml are upper case
+            auto ucName = boost::algorithm::to_upper_copy(item.first);
+            if (dbSpecies.count(ucName)) {
+                auto& spec = *dbSpecies.at(ucName);
+                auto& critProps = spec["critical-parameters"].as<AnyMap>();
+                Tc = critProps.convert("critical-temperature", "K");
+                Pc = critProps.convert("critical-pressure", "Pa");
+                omega_ac = critProps["acentric-factor"].asDouble();
+                m_coeffSource[k] = CoeffSource::Database;
+            }
+        }
+
+        // Check if critical properties were found in either location
+        if (!isnan(Tc)) {
+            double a = omega_a * std::pow(GasConstant * Tc, 2) / Pc;
+            double b = omega_b * GasConstant * Tc / Pc;
+            setSpeciesCoeffs(item.first, a, b, omega_ac);
+        } else {
+            throw InputFileError("SoaveRedlichKwong::initThermo", data,
+            "No Soave-Redlich-Kwong model parameters or critical properties found for "
+            "species '{}'", item.first);
+        }
+    }
+}
+
+
+void SoaveRedlichKwong::getSpeciesParameters(const std::string& name,
+                                        AnyMap& speciesNode) const
+{
+    MixtureFugacityTP::getSpeciesParameters(name, speciesNode);
+    size_t k = speciesIndex(name);
+    checkSpeciesIndex(k);
+
+    // Pure species parameters
+    if (m_coeffSource[k] == CoeffSource::EoS) {
+        auto& eosNode = speciesNode["equation-of-state"].getMapWhere(
+            "model", "Soave-Redlich-Kwong", true);
+        eosNode["a"].setQuantity(m_a_coeffs(k, k), "Pa*m^6/kmol^2");
+        eosNode["b"].setQuantity(m_b_coeffs[k], "m^3/kmol");
+        eosNode["acentric-factor"] = m_acentric[k];
+    } else if (m_coeffSource[k] == CoeffSource::CritProps) {
+        auto& critProps = speciesNode["critical-parameters"];
+        double Tc = speciesCritTemperature(m_a_coeffs(k, k), m_b_coeffs[k]);
+        double Pc = omega_b * GasConstant * Tc / m_b_coeffs[k];
+        critProps["critical-temperature"].setQuantity(Tc, "K");
+        critProps["critical-pressure"].setQuantity(Pc, "Pa");
+        critProps["acentric-factor"] = m_acentric[k];
+    }
+    // Nothing to do in the case where the parameters are from the database
+
+    if (m_binaryParameters.count(name)) {
+        // Include binary parameters regardless of where the pure species parameters
+        // were found
+        auto& eosNode = speciesNode["equation-of-state"].getMapWhere(
+            "model", "Soave-Redlich-Kwong", true);
+        AnyMap bin_a;
+        for (const auto& item : m_binaryParameters.at(name)) {
+            bin_a[item.first].setQuantity(item.second, "Pa*m^6/kmol^2");
+        }
+        eosNode["binary-a"] = std::move(bin_a);
+    }
+}
 
 double SoaveRedlichKwong::dpdVCalc(double T, double molarVol, double& presCalc) const
 {

From 024f6f46158603c3e5365cc8b1c3c2d545f4f3d3 Mon Sep 17 00:00:00 2001
From: Cory Kinney <cory.kinney@ucf.edu>
Date: Thu, 5 Jan 2023 16:28:35 -0500
Subject: [PATCH 05/17] Add mixing functions

---
 include/cantera/thermo/SoaveRedlichKwong.h | 18 ++++++++++
 src/thermo/SoaveRedlichKwong.cpp           | 39 ++++++++++++++++++++--
 2 files changed, 54 insertions(+), 3 deletions(-)

diff --git a/include/cantera/thermo/SoaveRedlichKwong.h b/include/cantera/thermo/SoaveRedlichKwong.h
index 7ff6999cbf..25fb842d74 100644
--- a/include/cantera/thermo/SoaveRedlichKwong.h
+++ b/include/cantera/thermo/SoaveRedlichKwong.h
@@ -76,6 +76,24 @@ class SoaveRedlichKwong : public MixtureFugacityTP
      */
     void calculatePressureDerivatives() const;
 
+    //! Update the \f$a\f$, \f$b\f$, and \f$\alpha\f$ parameters
+    /*!
+     *  The \f$a\f$ and the \f$b\f$ parameters depend on the mole fraction and the
+     *  parameter \f$\alpha\f$ depends on the temperature. This function updates
+     *  the internal numbers based on the state of the object.
+     */
+    virtual void updateMixingExpressions();
+
+    //! Calculate the \f$a\f$, \f$b\f$, and \f$\a alpha\f$ parameters given the temperature
+    /*!
+     * This function doesn't change the internal state of the object, so it is a
+     * const function.  It does use the stored mole fractions in the object.
+     *
+     * @param aCalc (output)  Returns the a value
+     * @param bCalc (output)  Returns the b value.
+     * @param aAlpha (output) Returns the (a*alpha) value.
+     */
+    void calculateAB(double& aCalc, double& bCalc, double& aAlphaCalc) const
     //! Prepare variables and call the function to solve the cubic equation of state
     int solveCubic(double T, double pres, double a, double b, double aAlpha,
                    double Vroot[3]) const;
diff --git a/src/thermo/SoaveRedlichKwong.cpp b/src/thermo/SoaveRedlichKwong.cpp
index ae1feca931..4d45fa7a8a 100644
--- a/src/thermo/SoaveRedlichKwong.cpp
+++ b/src/thermo/SoaveRedlichKwong.cpp
@@ -30,7 +30,6 @@ double SoaveRedlichKwong::speciesCritTemperature(double a, double b) const
     }
 }
 
-
 bool SoaveRedlichKwong::addSpecies(shared_ptr<Species> spec)
 {
     bool added = MixtureFugacityTP::addSpecies(spec);
@@ -51,7 +50,6 @@ bool SoaveRedlichKwong::addSpecies(shared_ptr<Species> spec)
     return added;
 }
 
-
 void SoaveRedlichKwong::initThermo()
 {
     // Contents of 'critical-properties.yaml', loaded later if needed
@@ -137,7 +135,6 @@ void SoaveRedlichKwong::initThermo()
     }
 }
 
-
 void SoaveRedlichKwong::getSpeciesParameters(const std::string& name,
                                         AnyMap& speciesNode) const
 {
@@ -192,6 +189,42 @@ void SoaveRedlichKwong::calculatePressureDerivatives() const
     m_dpdT = GasConstant / (mv - m_b) - daAlpha_dT() / (mv * (mv + m_b));
 }
 
+void SoaveRedlichKwong::updateMixingExpressions()
+{
+    // Same as PengRobinson::updateMixingExpressions()
+    double temp = temperature();
+
+    // Update individual alpha
+    for (size_t j = 0; j < m_kk; j++) {
+        double critTemp_j = speciesCritTemperature(m_a_coeffs(j,j), m_b_coeffs[j]);
+        double sqt_alpha = 1 + m_kappa[j] * (1 - sqrt(temp / critTemp_j));
+        m_alpha[j] = sqt_alpha*sqt_alpha;
+    }
+
+    // Update aAlpha_i, j
+    for (size_t i = 0; i < m_kk; i++) {
+        for (size_t j = 0; j < m_kk; j++) {
+            m_aAlpha_binary(i, j) = sqrt(m_alpha[i] * m_alpha[j]) * m_a_coeffs(i,j);
+        }
+    }
+    calculateAB(m_a,m_b,m_aAlpha_mix);
+}
+
+void SoaveRedlichKwong::calculateAB(double& aCalc, double& bCalc, double& aAlphaCalc) const
+{
+    // Same as PengRobinson::calculateAB()
+    bCalc = 0.0;
+    aCalc = 0.0;
+    aAlphaCalc = 0.0;
+    for (size_t i = 0; i < m_kk; i++) {
+        bCalc += moleFractions_[i] * m_b_coeffs[i];
+        for (size_t j = 0; j < m_kk; j++) {
+            aCalc += m_a_coeffs(i, j) * moleFractions_[i] * moleFractions_[j];
+            aAlphaCalc += m_aAlpha_binary(i, j) * moleFractions_[i] * moleFractions_[j];
+        }
+    }
+}
+
 double SoaveRedlichKwong::daAlpha_dT() const
 {
     // Same as PengRobinson::daAlpha_dT()

From 54d7857e4ad0661739d3a80a892943d2ff68827b Mon Sep 17 00:00:00 2001
From: Cory Kinney <cory.kinney@ucf.edu>
Date: Thu, 5 Jan 2023 16:54:33 -0500
Subject: [PATCH 06/17] Add density functions

Copied from `PengRobinson`/`RedlichKwongMFTP`
Also copied over other variable declarations to the header file
---
 include/cantera/thermo/SoaveRedlichKwong.h |  51 +++++++-
 src/thermo/SoaveRedlichKwong.cpp           | 133 +++++++++++++++++++++
 2 files changed, 183 insertions(+), 1 deletion(-)

diff --git a/include/cantera/thermo/SoaveRedlichKwong.h b/include/cantera/thermo/SoaveRedlichKwong.h
index 25fb842d74..6285ee8b0d 100644
--- a/include/cantera/thermo/SoaveRedlichKwong.h
+++ b/include/cantera/thermo/SoaveRedlichKwong.h
@@ -52,7 +52,13 @@ class SoaveRedlichKwong : public MixtureFugacityTP
                                       AnyMap& speciesNode) const;
 
 protected:
+    // Special functions inherited from MixtureFugacityTP
 
+    virtual double liquidVolEst(double T, double& pres) const;
+    virtual double densityCalc(double T, double pressure, int phase, double rhoguess);
+
+    virtual double densSpinodalLiquid() const;
+    virtual double densSpinodalGas() const;
     virtual double dpdVCalc(double T, double molarVol, double& presCalc) const;
 
     // Special functions not inherited from MixtureFugacityTP
@@ -117,6 +123,33 @@ class SoaveRedlichKwong : public MixtureFugacityTP
      */
     double m_aAlpha_mix;
 
+    // Vectors required to store species-specific a_coeff, b_coeff, alpha, kappa
+    // and other derivatives. Length = m_kk.
+    vector_fp m_b_coeffs;
+    vector_fp m_kappa;
+    vector_fp m_acentric; //!< acentric factor for each species, length #m_kk
+    mutable vector_fp m_dalphadT;
+    mutable vector_fp m_d2alphadT2;
+    vector_fp m_alpha;
+
+    // Matrices for Binary coefficients a_{i,j} and {a*alpha}_{i.j} are saved in an
+    // array form. Size = (m_kk, m_kk).
+    Array2D m_a_coeffs;
+    Array2D m_aAlpha_binary;
+
+    //! Explicitly-specified binary interaction parameters, to enable serialization
+    std::map<std::string, std::map<std::string, double>> m_binaryParameters;
+
+    int m_NSolns;
+
+    double m_Vroot[3];
+
+    //! Temporary storage - length = m_kk.
+    mutable vector_fp m_pp;
+
+    // Partial molar volumes of the species
+    mutable vector_fp m_partialMolarVolumes;
+
     //! The derivative of the pressure with respect to the volume
     /*!
      * Calculated at the current conditions. temperature and mole number kept
@@ -131,16 +164,32 @@ class SoaveRedlichKwong : public MixtureFugacityTP
      */
     mutable double m_dpdT;
 
+    //! Vector of derivatives of pressure with respect to mole number
+    /*!
+     *  Calculated at the current conditions. Total volume, temperature and
+     *  other mole number kept constant
+     */
+    mutable vector_fp m_dpdni;
 
     enum class CoeffSource { EoS, CritProps, Database };
     //! For each species, specifies the source of the a, b, and omega coefficients
     std::vector<CoeffSource> m_coeffSource;
 
 private:
+    //! Omega constant: omega_a used in Soave-Redlich-Kwong equation of state
+    /*!
+     *  This value is calculated by solving SRK cubic equation at the critical point.
+     */
     static const double omega_a;
+
+    //! Omega constant: omega_b used in Soave-Redlich-Kwong equation of state
+    /*!
+     *  This value is calculated by solving SRK cubic equation at the critical point.
+     */
     static const double omega_b;
-    static const double omega_vc;
 
+    //! Omega constant for the critical molar volume
+    static const double omega_vc;
 };
 }
 
diff --git a/src/thermo/SoaveRedlichKwong.cpp b/src/thermo/SoaveRedlichKwong.cpp
index 4d45fa7a8a..db8a17ed08 100644
--- a/src/thermo/SoaveRedlichKwong.cpp
+++ b/src/thermo/SoaveRedlichKwong.cpp
@@ -172,6 +172,139 @@ void SoaveRedlichKwong::getSpeciesParameters(const std::string& name,
     }
 }
 
+double SoaveRedlichKwong::liquidVolEst(double T, double& presGuess) const
+{
+    double v = m_b * 1.1;
+    double atmp, btmp, aAlphatmp;
+    calculateAB(atmp, btmp, aAlphatmp);
+    double pres = std::max(psatEst(T), presGuess);
+    double Vroot[3];
+    bool foundLiq = false;
+    int m = 0;
+    while (m < 100 && !foundLiq) {
+        int nsol = solveCubic(T, pres, atmp, btmp, aAlphatmp, Vroot);
+        if (nsol == 1 || nsol == 2) {
+            double pc = critPressure();
+            if (pres > pc) {
+                foundLiq = true;
+            }
+            pres *= 1.04;
+        } else {
+            foundLiq = true;
+        }
+    }
+
+    if (foundLiq) {
+        v = Vroot[0];
+        presGuess = pres;
+    } else {
+        v = -1.0;
+    }
+    return v;
+}
+
+double SoaveRedlichKwong::densityCalc(double T, double presPa, int phaseRequested,
+                                 double rhoGuess)
+{
+    // It's necessary to set the temperature so that m_aAlpha_mix is set correctly.
+    setTemperature(T);
+    double tcrit = critTemperature();
+    double mmw = meanMolecularWeight();
+    if (rhoGuess == -1.0) {
+        if (phaseRequested >= FLUID_LIQUID_0) {
+            double lqvol = liquidVolEst(T, presPa);
+            rhoGuess = mmw / lqvol;
+        }
+    } else {
+        // Assume the Gas phase initial guess, if nothing is specified to the routine
+        rhoGuess = presPa * mmw / (GasConstant * T);
+    }
+
+    double volGuess = mmw / rhoGuess;
+    m_NSolns = solveCubic(T, presPa, m_a, m_b, m_aAlpha_mix, m_Vroot);
+
+    double molarVolLast = m_Vroot[0];
+    if (m_NSolns >= 2) {
+        if (phaseRequested >= FLUID_LIQUID_0) {
+            molarVolLast = m_Vroot[0];
+        } else if (phaseRequested == FLUID_GAS || phaseRequested == FLUID_SUPERCRIT) {
+            molarVolLast = m_Vroot[2];
+        } else {
+            if (volGuess > m_Vroot[1]) {
+                molarVolLast = m_Vroot[2];
+            } else {
+                molarVolLast = m_Vroot[0];
+            }
+        }
+    } else if (m_NSolns == 1) {
+        if (phaseRequested == FLUID_GAS || phaseRequested == FLUID_SUPERCRIT
+            || phaseRequested == FLUID_UNDEFINED)
+        {
+            molarVolLast = m_Vroot[0];
+        } else {
+            return -2.0;
+        }
+    } else if (m_NSolns == -1) {
+        if (phaseRequested >= FLUID_LIQUID_0 || phaseRequested == FLUID_UNDEFINED
+            || phaseRequested == FLUID_SUPERCRIT)
+        {
+            molarVolLast = m_Vroot[0];
+        } else if (T > tcrit) {
+            molarVolLast = m_Vroot[0];
+        } else {
+            return -2.0;
+        }
+    } else {
+        molarVolLast = m_Vroot[0];
+        return -1.0;
+    }
+    return mmw / molarVolLast;
+}
+
+double SoaveRedlichKwong::densSpinodalLiquid() const
+{
+    double Vroot[3];
+    double T = temperature();
+    int nsol = solveCubic(T, pressure(), m_a, m_b, m_aAlpha_mix, Vroot);
+    if (nsol != 3) {
+        return critDensity();
+    }
+
+    auto resid = [this, T](double v) {
+        double pp;
+        return dpdVCalc(T, v, pp);
+    };
+
+    boost::uintmax_t maxiter = 100;
+    std::pair<double, double> vv = bmt::toms748_solve(
+        resid, Vroot[0], Vroot[1], bmt::eps_tolerance<double>(48), maxiter);
+
+    double mmw = meanMolecularWeight();
+    return mmw / (0.5 * (vv.first + vv.second));
+}
+
+double SoaveRedlichKwong::densSpinodalGas() const
+{
+    double Vroot[3];
+    double T = temperature();
+    int nsol = solveCubic(T, pressure(), m_a, m_b, m_aAlpha_mix, Vroot);
+    if (nsol != 3) {
+        return critDensity();
+    }
+
+    auto resid = [this, T](double v) {
+        double pp;
+        return dpdVCalc(T, v, pp);
+    };
+
+    boost::uintmax_t maxiter = 100;
+    std::pair<double, double> vv = bmt::toms748_solve(
+        resid, Vroot[1], Vroot[2], bmt::eps_tolerance<double>(48), maxiter);
+
+    double mmw = meanMolecularWeight();
+    return mmw / (0.5 * (vv.first + vv.second));
+}
+
 double SoaveRedlichKwong::dpdVCalc(double T, double molarVol, double& presCalc) const
 {
     double denom = molarVol * (molarVol + m_b);

From d2ad5eabaea389ad00d97db98e3d388297453fef Mon Sep 17 00:00:00 2001
From: corykinney <cory.kinney@ucf.edu>
Date: Mon, 9 Jan 2023 17:04:38 -0500
Subject: [PATCH 07/17] Add isothermal compressibility and thermal expansion
 coefficient

---
 include/cantera/thermo/SoaveRedlichKwong.h |  6 ++++
 src/thermo/SoaveRedlichKwong.cpp           | 36 ++++++++++++++++++++++
 2 files changed, 42 insertions(+)

diff --git a/include/cantera/thermo/SoaveRedlichKwong.h b/include/cantera/thermo/SoaveRedlichKwong.h
index 6285ee8b0d..15b168b7c3 100644
--- a/include/cantera/thermo/SoaveRedlichKwong.h
+++ b/include/cantera/thermo/SoaveRedlichKwong.h
@@ -76,6 +76,12 @@ class SoaveRedlichKwong : public MixtureFugacityTP
     double d2aAlpha_dT2() const;
 
 public:
+    //! Returns the isothermal compressibility. Units: 1/Pa.
+    double isothermalCompressibility() const;
+
+    //! Return the volumetric thermal expansion coefficient. Units: 1/K.
+    double thermalExpansionCoeff() const;
+
     //! Calculate \f$dp/dV\f$ and \f$dp/dT\f$ at the current conditions
     /*!
      *  These are stored internally.
diff --git a/src/thermo/SoaveRedlichKwong.cpp b/src/thermo/SoaveRedlichKwong.cpp
index db8a17ed08..bb0190a916 100644
--- a/src/thermo/SoaveRedlichKwong.cpp
+++ b/src/thermo/SoaveRedlichKwong.cpp
@@ -312,6 +312,42 @@ double SoaveRedlichKwong::dpdVCalc(double T, double molarVol, double& presCalc)
     return -GasConstant * T / (vmb * vmb) + m_aAlpha_mix * (2 * molarVol + m_b) / (denom * denom);
 }
 
+double SoaveRedlichKwong::isothermalCompressibility() const 
+{
+    double P = pressure();
+    double mv = molarVolume();
+    double T = temperature();
+    double RT_ = RT();
+    double Z = P * mv / RT_;
+
+    double A = m_aAlpha_mix * P / (RT_ * RT_);
+    double B = m_b * P / RT_;
+
+    double dAdP = A / P;
+    double dBdP = B / P;
+    double dZdP = ((B - Z) * dAdP + (A + B + 2 * B * Z) * dBdP) / (3 * Z * Z - 2 * Z + A - B - B * B);
+
+    return 1 / P - dZdP / Z;
+}
+
+double SoaveRedlichKwong::thermalExpansionCoeff() const 
+{
+    double P = pressure();
+    double mv = molarVolume();
+    double T = temperature();
+    double RT_ = RT();
+    double Z = P * mv / RT_;
+
+    double A = m_aAlpha_mix * P / (RT_ * RT_);
+    double B = m_b * P / RT_;
+
+    double dAdT = P / (RT_ * RT_) * (daAlpha_dT() - 2 * m_aAlpha_mix / T);
+    double dBdT = -B / T;
+    double dZdP = ((B - Z) * dAdT + (A + B + 2 * B * Z) * dBdT) / (3 * Z * Z - 2 * Z + A - B - B * B);
+
+    return 1 / T + dZdT / Z;
+}
+
 void SoaveRedlichKwong::calculatePressureDerivatives() const
 {
     double T = temperature();

From f64df51722c102a93db9707d75ace7a0e7dfed48 Mon Sep 17 00:00:00 2001
From: corykinney <cory.kinney@ucf.edu>
Date: Tue, 10 Jan 2023 09:46:28 -0500
Subject: [PATCH 08/17] Add pressure

---
 include/cantera/thermo/SoaveRedlichKwong.h | 44 ++++++++++++++++++++++
 src/thermo/SoaveRedlichKwong.cpp           |  9 +++++
 2 files changed, 53 insertions(+)

diff --git a/include/cantera/thermo/SoaveRedlichKwong.h b/include/cantera/thermo/SoaveRedlichKwong.h
index 15b168b7c3..b7fbf94ebb 100644
--- a/include/cantera/thermo/SoaveRedlichKwong.h
+++ b/include/cantera/thermo/SoaveRedlichKwong.h
@@ -26,6 +26,50 @@ class SoaveRedlichKwong : public MixtureFugacityTP
         return "Soave-Redlich-Kwong";
     }
 
+    //! @}
+    //! @name Mechanical Properties
+    //! @{
+
+    //! Return the thermodynamic pressure (Pa).
+    /*!
+     * Since the mass density, temperature, and mass fractions are stored,
+     * this method uses these values to implement the
+     * mechanical equation of state \f$ P(T, \rho, Y_1, \dots, Y_K) \f$.
+     *
+     * \f[
+     *    P = \frac{RT}{v-b_{mix}}
+     *        - \frac{\left(a\alpha\right)_{mix}}{v\left(v + b_{mix}\right)}
+     * \f]
+     *
+     * where:
+     *
+     * \f[
+     *    \alpha = \left[ 1 + \kappa \left(1-T_r^{0.5}\right)\right]^2
+     * \f]
+     *
+     *  and
+     *
+     * \f[
+     *    \kappa = \left(0.48508 + 1.55171\omega - 0.15613\omega^2\right)
+     * \f]
+     *
+     * Coefficients \f$ a_{mix}, b_{mix} \f$ and \f$(a \alpha)_{mix}\f$ are calculated as
+     *
+     * \f[
+     *    a_{mix} = \sum_i \sum_j X_i X_j a_{i, j} = \sum_i \sum_j X_i X_j \sqrt{a_i a_j}
+     * \f]
+     *
+     * \f[
+     *    b_{mix} = \sum_i X_i b_i
+     * \f]
+     *
+     * \f[
+     *   {a \alpha}_{mix} = \sum_i \sum_j X_i X_j {a \alpha}_{i, j}
+     *       = \sum_i \sum_j X_i X_j \sqrt{a_i a_j} \sqrt{\alpha_i \alpha_j}
+     * \f]
+     */
+    virtual double pressure() const;
+
 
     //! Calculate species-specific critical temperature
     /*!
diff --git a/src/thermo/SoaveRedlichKwong.cpp b/src/thermo/SoaveRedlichKwong.cpp
index bb0190a916..ab58597e38 100644
--- a/src/thermo/SoaveRedlichKwong.cpp
+++ b/src/thermo/SoaveRedlichKwong.cpp
@@ -19,6 +19,15 @@ const double SoaveRedlichKwong::omega_b = 8.66403499650E-02;
 const double SoaveRedlichKwong::omega_vc = 3.33333333333333E-01;
 
 
+
+double SoaveRedlichKwong::pressure() const
+{
+    _updateReferenceStateThermo();
+    // Get a copy of the private variables stored in the State object
+    double T = temperature();
+    double mv = molarVolume();
+    return GasConstant * T / (mv - m_b) - m_aAlpha_mix / (mv * (mv - m_b))
+}
 double SoaveRedlichKwong::speciesCritTemperature(double a, double b) const
 {
     if (b <= 0.0) {

From d986d398caae3505aded5f54f4f74e672018e9d4 Mon Sep 17 00:00:00 2001
From: corykinney <cory.kinney@ucf.edu>
Date: Tue, 10 Jan 2023 09:47:44 -0500
Subject: [PATCH 09/17] Add constructor and coefficient functions

---
 src/thermo/SoaveRedlichKwong.cpp | 75 ++++++++++++++++++++++++++++++++
 1 file changed, 75 insertions(+)

diff --git a/src/thermo/SoaveRedlichKwong.cpp b/src/thermo/SoaveRedlichKwong.cpp
index ab58597e38..b12ee14e66 100644
--- a/src/thermo/SoaveRedlichKwong.cpp
+++ b/src/thermo/SoaveRedlichKwong.cpp
@@ -19,6 +19,81 @@ const double SoaveRedlichKwong::omega_b = 8.66403499650E-02;
 const double SoaveRedlichKwong::omega_vc = 3.33333333333333E-01;
 
 
+SoaveRedlichKwong::SoaveRedlichKwong(const std::string& infile="", const std::string& id="")  :
+    m_b(0.0),
+    m_a(0.0),
+    m_aAlpha_mix(0.0),
+    m_NSolns(0),
+    m_dpdV(0.0),
+    m_dpdT(0.0)
+{
+    std::fill_n(m_Vroot, 3, 0.0);
+    initThermoFile(infile, id_);
+}
+
+void SoaveRedlichKwong::setSpeciesCoeffs(const std::string& species, double a, double b,
+                                         double w)
+{
+    size_t k = speciesIndex(species);
+    if (k == npos) {
+        throw CanteraError("SoaveRedlichKwong::setSpeciesCoeffs",
+            "Unknown species '{}'.", species);
+    }
+
+    // Calculate value of kappa (independent of temperature)
+    // w is an acentric factor of species
+    m_kappa[k] = 0.48508 + 1.55171*w - 0.15613*w*w;
+    m_acentric[k] = w; // store the original acentric factor to enable serialization
+
+    // Calculate alpha (temperature dependent interaction parameter)
+    double critTemp = speciesCritTemperature(a, b);
+    double sqt_T_r = sqrt(temperature() / critTemp);
+    double sqt_alpha = 1 + m_kappa[k] * (1 - sqt_T_r);
+    m_alpha[k] = sqt_alpha*sqt_alpha;
+
+    m_a_coeffs(k,k) = a;
+    double aAlpha_k = a*m_alpha[k];
+    m_aAlpha_binary(k,k) = aAlpha_k;
+
+    // standard mixing rule for cross-species interaction term
+    for (size_t j = 0; j < m_kk; j++) {
+        if (k == j) {
+            continue;
+        }
+        double a0kj = sqrt(m_a_coeffs(j,j) * a);
+        double aAlpha_j = a*m_alpha[j];
+        double a_Alpha = sqrt(aAlpha_j*aAlpha_k);
+        if (m_a_coeffs(j, k) == 0) {
+            m_a_coeffs(j, k) = a0kj;
+            m_aAlpha_binary(j, k) = a_Alpha;
+            m_a_coeffs(k, j) = a0kj;
+            m_aAlpha_binary(k, j) = a_Alpha;
+        }
+    }
+    m_b_coeffs[k] = b;
+}
+
+void SoaveRedlichKwong::setBinaryCoeffs(const std::string& species_i,
+        const std::string& species_j, double a0)
+{
+    size_t ki = speciesIndex(species_i);
+    if (ki == npos) {
+        throw CanteraError("SoaveRedlichKwong::setBinaryCoeffs",
+            "Unknown species '{}'.", species_i);
+    }
+    size_t kj = speciesIndex(species_j);
+    if (kj == npos) {
+        throw CanteraError("SoaveRedlichKwong::setBinaryCoeffs",
+            "Unknown species '{}'.", species_j);
+    }
+
+    m_a_coeffs(ki, kj) = m_a_coeffs(kj, ki) = a0;
+    m_binaryParameters[species_i][species_j] = a0;
+    m_binaryParameters[species_j][species_i] = a0;
+    // Calculate alpha_ij
+    double alpha_ij = m_alpha[ki] * m_alpha[kj];
+    m_aAlpha_binary(ki, kj) = m_aAlpha_binary(kj, ki) = a0*alpha_ij;
+}
 
 double SoaveRedlichKwong::pressure() const
 {

From f4faa34a3d8c1f6fa4f1fe1ca11ce7a295c352c0 Mon Sep 17 00:00:00 2001
From: corykinney <cory.kinney@ucf.edu>
Date: Tue, 10 Jan 2023 09:58:34 -0500
Subject: [PATCH 10/17] Add calculation of critical conditions

---
 include/cantera/thermo/SoaveRedlichKwong.h |  3 +++
 src/thermo/SoaveRedlichKwong.cpp           | 19 +++++++++++++++++++
 2 files changed, 22 insertions(+)

diff --git a/include/cantera/thermo/SoaveRedlichKwong.h b/include/cantera/thermo/SoaveRedlichKwong.h
index b7fbf94ebb..ce60029045 100644
--- a/include/cantera/thermo/SoaveRedlichKwong.h
+++ b/include/cantera/thermo/SoaveRedlichKwong.h
@@ -150,6 +150,9 @@ class SoaveRedlichKwong : public MixtureFugacityTP
      * @param aAlpha (output) Returns the (a*alpha) value.
      */
     void calculateAB(double& aCalc, double& bCalc, double& aAlphaCalc) const
+
+    void calcCriticalConditions(double& pc, double& tc, double& vc) const;
+
     //! Prepare variables and call the function to solve the cubic equation of state
     int solveCubic(double T, double pres, double a, double b, double aAlpha,
                    double Vroot[3]) const;
diff --git a/src/thermo/SoaveRedlichKwong.cpp b/src/thermo/SoaveRedlichKwong.cpp
index b12ee14e66..202292d6f4 100644
--- a/src/thermo/SoaveRedlichKwong.cpp
+++ b/src/thermo/SoaveRedlichKwong.cpp
@@ -103,6 +103,7 @@ double SoaveRedlichKwong::pressure() const
     double mv = molarVolume();
     return GasConstant * T / (mv - m_b) - m_aAlpha_mix / (mv * (mv - m_b))
 }
+
 double SoaveRedlichKwong::speciesCritTemperature(double a, double b) const
 {
     if (b <= 0.0) {
@@ -534,6 +535,24 @@ double SoaveRedlichKwong::d2aAlpha_dT2() const
     return d2aAlphadT2;
 }
 
+void SoaveRedlichKwong::calcCriticalConditions(double& pc, double& tc, double& vc) const
+{
+    if (m_b <= 0.0) {
+        tc = 1000000.;
+        pc = 1.0E13;
+        vc = omega_vc * GasConstant * tc / pc;
+        return;
+    }
+    if (m_a <= 0.0) {
+        tc = 0.0;
+        pc = 0.0;
+        vc = 2.0 * m_b;
+        return;
+    }
+    tc = m_a * omega_b / (m_b * omega_a * GasConstant);
+    pc = omega_b * GasConstant * tc / m_b;
+    vc = omega_vc * GasConstant * tc / pc;
+}
 
 int SoaveRedlichKwong::solveCubic(double T, double pres, double a, double b, double aAlpha,
                                   double Vroot[3]) const

From 7f77f153066e297b472a74a7fc5fd540e56ee724 Mon Sep 17 00:00:00 2001
From: corykinney <cory.kinney@ucf.edu>
Date: Tue, 10 Jan 2023 09:59:08 -0500
Subject: [PATCH 11/17] Fix missing declarations and comments

---
 include/cantera/thermo/SoaveRedlichKwong.h | 33 ++++++++++++++++++++--
 1 file changed, 31 insertions(+), 2 deletions(-)

diff --git a/include/cantera/thermo/SoaveRedlichKwong.h b/include/cantera/thermo/SoaveRedlichKwong.h
index ce60029045..d623bbe704 100644
--- a/include/cantera/thermo/SoaveRedlichKwong.h
+++ b/include/cantera/thermo/SoaveRedlichKwong.h
@@ -19,6 +19,15 @@ namespace Cantera
 class SoaveRedlichKwong : public MixtureFugacityTP
 {
 public:
+
+    //! Construct and initialize a SoaveRedlichKwong object directly from an
+    //! input file
+    /*!
+     * @param infile    Name of the input file containing the phase YAML data.
+     *                  If blank, an empty phase will be created.
+     * @param id        ID of the phase in the input file. If empty, the
+     *                  first phase definition in the input file will be used.
+     */
     explicit SoaveRedlichKwong(const std::string& infile="",
                                const std::string& id="");
 
@@ -77,8 +86,8 @@ class SoaveRedlichKwong : public MixtureFugacityTP
      *       \f[ T_{crit} = (0.0778 a)/(0.4572 b R) \f] TODO
      *  Units: Kelvin
      *
-     * @param a    species-specific coefficients used in P-R EoS
-     * @param b    species-specific coefficients used in P-R EoS
+     * @param a    species-specific coefficients used in SRK EoS
+     * @param b    species-specific coefficients used in SRK EoS
      */
     double speciesCritTemperature(double a, double b) const;
 
@@ -95,6 +104,26 @@ class SoaveRedlichKwong : public MixtureFugacityTP
     virtual void getSpeciesParameters(const std::string& name,
                                       AnyMap& speciesNode) const;
 
+    //! Set the pure fluid interaction parameters for a species
+    /*!
+     *  @param species   Name of the species
+     *  @param a         \f$a\f$ parameter in the Soave-Redlich-Kwong model [Pa-m^6/kmol^2]
+     *  @param b         \f$a\f$ parameter in the Soave-Redlich-Kwong model [m^3/kmol]
+     *  @param w         acentric factor
+     */
+    void setSpeciesCoeffs(const std::string& species, double a, double b,
+                          double w);
+
+    //! Set values for the interaction parameter between two species
+    /*!
+     *  @param species_i   Name of one species
+     *  @param species_j   Name of the other species
+     *  @param a           \f$a\f$ parameter in the Soave-Redlich-Kwong model [Pa-m^6/kmol^2]
+     */
+    void setBinaryCoeffs(const std::string& species_i,
+                         const std::string& species_j, double a);
+    //! @}
+
 protected:
     // Special functions inherited from MixtureFugacityTP
 

From c918d10b4264310a90aef68da5022d8130a53766 Mon Sep 17 00:00:00 2001
From: corykinney <cory.kinney@ucf.edu>
Date: Tue, 10 Jan 2023 10:01:05 -0500
Subject: [PATCH 12/17] Fix pressure equation typo

---
 src/thermo/SoaveRedlichKwong.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/thermo/SoaveRedlichKwong.cpp b/src/thermo/SoaveRedlichKwong.cpp
index 202292d6f4..4a45d82917 100644
--- a/src/thermo/SoaveRedlichKwong.cpp
+++ b/src/thermo/SoaveRedlichKwong.cpp
@@ -101,7 +101,7 @@ double SoaveRedlichKwong::pressure() const
     // Get a copy of the private variables stored in the State object
     double T = temperature();
     double mv = molarVolume();
-    return GasConstant * T / (mv - m_b) - m_aAlpha_mix / (mv * (mv - m_b))
+    return GasConstant * T / (mv - m_b) - m_aAlpha_mix / (mv * (mv + m_b))
 }
 
 double SoaveRedlichKwong::speciesCritTemperature(double a, double b) const

From cb9b292c40a370b7d46848e6abb4527c602794ee Mon Sep 17 00:00:00 2001
From: corykinney <cory.kinney@ucf.edu>
Date: Tue, 10 Jan 2023 17:14:53 -0500
Subject: [PATCH 13/17] Cleanup and fix errors

Add `SoaveRedlichKwong` to `ThermoFactory`
Change coefficients to original SRK values
Fix build issues
Fix typo in `isothermal_compressibility` and `thermal_expansion_coeff` and remove unused variables
---
 src/thermo/SoaveRedlichKwong.cpp | 21 +++++++++------------
 src/thermo/ThermoFactory.cpp     |  2 ++
 2 files changed, 11 insertions(+), 12 deletions(-)

diff --git a/src/thermo/SoaveRedlichKwong.cpp b/src/thermo/SoaveRedlichKwong.cpp
index 4a45d82917..4eed854a1c 100644
--- a/src/thermo/SoaveRedlichKwong.cpp
+++ b/src/thermo/SoaveRedlichKwong.cpp
@@ -14,12 +14,12 @@ namespace bmt = boost::math::tools;
 
 namespace Cantera
 {
+
 const double SoaveRedlichKwong::omega_a = 4.27480233540E-01;
 const double SoaveRedlichKwong::omega_b = 8.66403499650E-02;
 const double SoaveRedlichKwong::omega_vc = 3.33333333333333E-01;
 
-
-SoaveRedlichKwong::SoaveRedlichKwong(const std::string& infile="", const std::string& id="")  :
+SoaveRedlichKwong::SoaveRedlichKwong(const std::string& infile, const std::string& id_)  :
     m_b(0.0),
     m_a(0.0),
     m_aAlpha_mix(0.0),
@@ -42,7 +42,7 @@ void SoaveRedlichKwong::setSpeciesCoeffs(const std::string& species, double a, d
 
     // Calculate value of kappa (independent of temperature)
     // w is an acentric factor of species
-    m_kappa[k] = 0.48508 + 1.55171*w - 0.15613*w*w;
+    m_kappa[k] = 0.480 + 1.574*w - 0.176*w*w;
     m_acentric[k] = w; // store the original acentric factor to enable serialization
 
     // Calculate alpha (temperature dependent interaction parameter)
@@ -400,17 +400,15 @@ double SoaveRedlichKwong::dpdVCalc(double T, double molarVol, double& presCalc)
 double SoaveRedlichKwong::isothermalCompressibility() const 
 {
     double P = pressure();
-    double mv = molarVolume();
-    double T = temperature();
     double RT_ = RT();
-    double Z = P * mv / RT_;
+    double Z = z();
 
     double A = m_aAlpha_mix * P / (RT_ * RT_);
     double B = m_b * P / RT_;
 
     double dAdP = A / P;
     double dBdP = B / P;
-    double dZdP = ((B - Z) * dAdP + (A + B + 2 * B * Z) * dBdP) / (3 * Z * Z - 2 * Z + A - B - B * B);
+    double dZdP = ((B - Z) * dAdP + (A + Z + 2 * B * Z) * dBdP) / (3 * Z * Z - 2 * Z + A - B - B * B);
 
     return 1 / P - dZdP / Z;
 }
@@ -418,17 +416,16 @@ double SoaveRedlichKwong::isothermalCompressibility() const
 double SoaveRedlichKwong::thermalExpansionCoeff() const 
 {
     double P = pressure();
-    double mv = molarVolume();
     double T = temperature();
     double RT_ = RT();
-    double Z = P * mv / RT_;
+    double Z = z();
 
     double A = m_aAlpha_mix * P / (RT_ * RT_);
     double B = m_b * P / RT_;
 
     double dAdT = P / (RT_ * RT_) * (daAlpha_dT() - 2 * m_aAlpha_mix / T);
     double dBdT = -B / T;
-    double dZdP = ((B - Z) * dAdT + (A + B + 2 * B * Z) * dBdT) / (3 * Z * Z - 2 * Z + A - B - B * B);
+    double dZdT = ((B - Z) * dAdT + (A + Z + 2 * B * Z) * dBdT) / (3 * Z * Z - 2 * Z + A - B - B * B);
 
     return 1 / T + dZdT / Z;
 }
@@ -560,10 +557,10 @@ int SoaveRedlichKwong::solveCubic(double T, double pres, double a, double b, dou
     double an = 1.0;
     double bn = - GasConstant * T / pres;
     double cn = (aAlpha - b * GasConstant * T) / pres - b * b;
-    double dn = aAlpha * b / pres;
+    double dn = -aAlpha * b / pres;
 
     double tc = a * omega_b / (b * omega_a * GasConstant);
-    double pc = omega_b * R * tc / b;
+    double pc = omega_b * GasConstant * tc / b;
     double vc = omega_vc * GasConstant * tc / pc;
 
     return MixtureFugacityTP::solveCubic(T, pres, a, b, aAlpha, Vroot,
diff --git a/src/thermo/ThermoFactory.cpp b/src/thermo/ThermoFactory.cpp
index fd6edb4c92..7a477d368e 100644
--- a/src/thermo/ThermoFactory.cpp
+++ b/src/thermo/ThermoFactory.cpp
@@ -25,6 +25,7 @@
 #include "cantera/thermo/PureFluidPhase.h"
 #include "cantera/thermo/RedlichKwongMFTP.h"
 #include "cantera/thermo/PengRobinson.h"
+#include "cantera/thermo/SoaveRedlichKwong.h"
 #include "cantera/thermo/SurfPhase.h"
 #include "cantera/thermo/EdgePhase.h"
 #include "cantera/thermo/MetalPhase.h"
@@ -104,6 +105,7 @@ ThermoFactory::ThermoFactory()
     reg("binary-solution-tabulated", []() { return new BinarySolutionTabulatedThermo(); });
     addAlias("binary-solution-tabulated", "BinarySolutionTabulatedThermo");
     reg("Peng-Robinson", []() { return new PengRobinson(); });
+    reg("Soave-Redlich-Kwong", []() { return new SoaveRedlichKwong(); });
 }
 
 ThermoPhase* ThermoFactory::newThermoPhase(const std::string& model)

From 74d7624ed0313647fbd2641b572d56bd0f4f4341 Mon Sep 17 00:00:00 2001
From: corykinney <cory.kinney@ucf.edu>
Date: Wed, 11 Jan 2023 08:21:21 -0500
Subject: [PATCH 14/17] Temporary

---
 include/cantera/thermo/SoaveRedlichKwong.h | 57 +++++++++++++-
 src/thermo/SoaveRedlichKwong.cpp           | 86 +++++++++++++++++++++-
 2 files changed, 139 insertions(+), 4 deletions(-)

diff --git a/include/cantera/thermo/SoaveRedlichKwong.h b/include/cantera/thermo/SoaveRedlichKwong.h
index d623bbe704..e1d56a8524 100644
--- a/include/cantera/thermo/SoaveRedlichKwong.h
+++ b/include/cantera/thermo/SoaveRedlichKwong.h
@@ -35,6 +35,12 @@ class SoaveRedlichKwong : public MixtureFugacityTP
         return "Soave-Redlich-Kwong";
     }
 
+    //! @name Molar Thermodynamic properties
+    //! @{
+
+    virtual double cp_mole() const;
+    virtual double cv_mole() const;
+
     //! @}
     //! @name Mechanical Properties
     //! @{
@@ -79,6 +85,49 @@ class SoaveRedlichKwong : public MixtureFugacityTP
      */
     virtual double pressure() const;
 
+    //! @}
+
+    //! Returns the standard concentration \f$ C^0_k \f$, which is used to
+    //! normalize the generalized concentration.
+    /*!
+     * This is defined as the concentration by which the generalized
+     * concentration is normalized to produce the activity.
+     * The ideal gas mixture is considered as the standard or reference state here.
+     * Since the activity for an ideal gas mixture is simply the mole fraction,
+     * for an ideal gas,  \f$ C^0_k = P/\hat R T \f$.
+     *
+     * @param k Optional parameter indicating the species. The default is to
+     *          assume this refers to species 0.
+     * @return
+     *   Returns the standard Concentration in units of m^3 / kmol.
+     */
+    virtual double standardConcentration(size_t k=0) const;
+
+    //! Get the array of non-dimensional activity coefficients at the current
+    //! solution temperature, pressure, and solution concentration.
+    /*!
+     * For all objects with the Mixture Fugacity approximation, we define the
+     * standard state as an ideal gas at the current temperature and pressure of
+     * the solution. The activities are based on this standard state.
+     *
+     * @param ac Output vector of activity coefficients. Length: m_kk.
+     */
+    virtual void getActivityCoefficients(double* ac) const;
+
+    //! @name  Partial Molar Properties of the Solution
+    //! @{
+
+    virtual void getChemPotentials(double* mu) const;
+    virtual void getPartialMolarEnthalpies(double* hbar) const;
+    virtual void getPartialMolarEntropies(double* sbar) const;
+    virtual void getPartialMolarIntEnergies(double* ubar) const;
+    //! Calculate species-specific molar specific heats
+    /*!
+     *  This function is currently not implemented for Soave-Redlich-Kwong phase.
+     */
+    virtual void getPartialMolarCp(double* cpbar) const;
+    virtual void getPartialMolarVolumes(double* vbar) const;
+    //! @}
 
     //! Calculate species-specific critical temperature
     /*!
@@ -89,7 +138,7 @@ class SoaveRedlichKwong : public MixtureFugacityTP
      * @param a    species-specific coefficients used in SRK EoS
      * @param b    species-specific coefficients used in SRK EoS
      */
-    double speciesCritTemperature(double a, double b) const;
+    virtual double speciesCritTemperature(double a, double b) const;
 
     //! @name Initialization Methods - For Internal use
     //!
@@ -126,6 +175,8 @@ class SoaveRedlichKwong : public MixtureFugacityTP
 
 protected:
     // Special functions inherited from MixtureFugacityTP
+    virtual double sresid() const;
+    virtual double hresid() const;
 
     virtual double liquidVolEst(double T, double& pres) const;
     virtual double densityCalc(double T, double pressure, int phase, double rhoguess);
@@ -169,7 +220,7 @@ class SoaveRedlichKwong : public MixtureFugacityTP
      */
     virtual void updateMixingExpressions();
 
-    //! Calculate the \f$a\f$, \f$b\f$, and \f$\a alpha\f$ parameters given the temperature
+    //! Calculate the \f$a\f$, \f$b\f$, and \f$a \alpha\f$ parameters given the temperature
     /*!
      * This function doesn't change the internal state of the object, so it is a
      * const function.  It does use the stored mole fractions in the object.
@@ -178,7 +229,7 @@ class SoaveRedlichKwong : public MixtureFugacityTP
      * @param bCalc (output)  Returns the b value.
      * @param aAlpha (output) Returns the (a*alpha) value.
      */
-    void calculateAB(double& aCalc, double& bCalc, double& aAlphaCalc) const
+    void calculateAB(double& aCalc, double& bCalc, double& aAlphaCalc) const;
 
     void calcCriticalConditions(double& pc, double& tc, double& vc) const;
 
diff --git a/src/thermo/SoaveRedlichKwong.cpp b/src/thermo/SoaveRedlichKwong.cpp
index 4eed854a1c..643d48a2f4 100644
--- a/src/thermo/SoaveRedlichKwong.cpp
+++ b/src/thermo/SoaveRedlichKwong.cpp
@@ -95,13 +95,71 @@ void SoaveRedlichKwong::setBinaryCoeffs(const std::string& species_i,
     m_aAlpha_binary(ki, kj) = m_aAlpha_binary(kj, ki) = a0*alpha_ij;
 }
 
+// ------------Molar Thermodynamic Properties -------------------------
+
+double SoaveRedlichKwong::cp_mole() const
+{
+    return 1.0;
+}
+
+double SoaveRedlichKwong::cv_mole() const
+{
+    _updateReferenceStateThermo();
+    double T = temperature();
+    calculatePressureDerivatives();
+    return (cp_mole() + T * m_dpdT * m_dpdT / m_dpdV);
+}
+
 double SoaveRedlichKwong::pressure() const
 {
     _updateReferenceStateThermo();
     // Get a copy of the private variables stored in the State object
     double T = temperature();
     double mv = molarVolume();
-    return GasConstant * T / (mv - m_b) - m_aAlpha_mix / (mv * (mv + m_b))
+    return GasConstant * T / (mv - m_b) - m_aAlpha_mix / (mv * (mv + m_b));
+}
+
+double SoaveRedlichKwong::standardConcentration(size_t k) const
+{
+    getStandardVolumes(m_tmpV.data());
+    return 1.0 / m_tmpV[k];
+}
+
+void SoaveRedlichKwong::getActivityCoefficients(double* ac) const
+{
+    
+}
+
+// ---- Partial Molar Properties of the Solution -----------------
+
+void SoaveRedlichKwong::getChemPotentials(double* mu) const
+{
+
+}
+
+void SoaveRedlichKwong::getPartialMolarEnthalpies(double* hbar) const
+{
+
+}
+
+void SoaveRedlichKwong::getPartialMolarEntropies(double* sbar) const
+{
+
+}
+
+void SoaveRedlichKwong::getPartialMolarIntEnergies(double* ubar) const
+{
+
+}
+
+void SoaveRedlichKwong::getPartialMolarCp(double* cpbar) const
+{
+
+}
+
+void SoaveRedlichKwong::getPartialMolarVolumes(double* vbar) const
+{
+
 }
 
 double SoaveRedlichKwong::speciesCritTemperature(double a, double b) const
@@ -257,6 +315,32 @@ void SoaveRedlichKwong::getSpeciesParameters(const std::string& name,
     }
 }
 
+double SoaveRedlichKwong::sresid() const
+{ // Replace
+    double molarV = molarVolume();
+    double hh = m_b / molarV;
+    double zz = z();
+    double alpha_1 = daAlpha_dT();
+    double vpb = molarV + (1.0 + Sqrt2) * m_b;
+    double vmb = molarV + (1.0 - Sqrt2) * m_b;
+    double fac = alpha_1 / (2.0 * Sqrt2 * m_b);
+    double sresid_mol_R = log(zz*(1.0 - hh)) + fac * log(vpb / vmb) / GasConstant;
+    return GasConstant * sresid_mol_R;
+}
+
+double SoaveRedlichKwong::hresid() const
+{ // Replace
+    double molarV = molarVolume();
+    double zz = z();
+    double aAlpha_1 = daAlpha_dT();
+    double T = temperature();
+    double vpb = molarV + (1 + Sqrt2) * m_b;
+    double vmb = molarV + (1 - Sqrt2) * m_b;
+    double fac = 1 / (2.0 * Sqrt2 * m_b);
+    return GasConstant * T * (zz - 1.0)
+        + fac * log(vpb / vmb) * (T * aAlpha_1 - m_aAlpha_mix);
+}
+
 double SoaveRedlichKwong::liquidVolEst(double T, double& presGuess) const
 {
     double v = m_b * 1.1;

From f323fa5c3f3d5aacb846a5bce5ee2967e4c52fec Mon Sep 17 00:00:00 2001
From: Cory Kinney <cory.kinney@ucf.edu>
Date: Fri, 13 Jan 2023 10:39:53 -0500
Subject: [PATCH 15/17] [Test] Switch real gas derivative tests to thermo
 consistency tests

Removed `getIsothermalCompressibility` and `getThermalExpansionCoeff` from `PengRobinson_Test` and `RedlichKwongMFTP_Test`.
Add finite difference tests to `TestConsistency` as `betaT_eq_minus_dmv_dP_const_T_div_mv` and `alphaV_eq_dmv_dT_const_P_div_mv`.
---
 test/thermo/consistency.cpp | 50 +++++++++++++++++++++++++++++++++++++
 1 file changed, 50 insertions(+)

diff --git a/test/thermo/consistency.cpp b/test/thermo/consistency.cpp
index 0004d24212..78a8d2602a 100644
--- a/test/thermo/consistency.cpp
+++ b/test/thermo/consistency.cpp
@@ -345,6 +345,56 @@ TEST_P(TestConsistency, dSdv_const_T_eq_dPdT_const_V) {
     }
 }
 
+TEST_P(TestConsistency, betaT_eq_minus_dmv_dP_const_T_div_mv)
+{
+    double betaT1;
+    try {
+        betaT1 = phase->isothermalCompressibility();
+    } catch (NotImplementedError& err) {
+        GTEST_SKIP() << err.getMethod() << " threw NotImplementedError";
+    }
+
+    double T = phase->temperature();
+    double P1 = phase->pressure();
+    double mv1 = phase->molarVolume();
+    
+    double P2 = P1 * (1 + 1e-7);
+    phase->setState_TP(T, P2);
+    double betaT2 = phase->isothermalCompressibility();
+    double mv2 = phase->molarVolume();
+
+    double betaT_mid = 0.5 * (betaT1 + betaT2);
+    double mv_mid = 0.5 * (mv1 + mv2);
+    double betaT_fd = -1 / mv_mid * (mv2 - mv1) / (P2 - P1);
+
+    EXPECT_NEAR(betaT_fd, betaT_mid, max({rtol_fd * betaT_mid, rtol_fd * betaT_fd, atol}));
+}
+
+TEST_P(TestConsistency, alphaV_eq_dmv_dT_const_P_div_mv)
+{
+    double alphaV1;
+    try {
+        alphaV1 = phase->thermalExpansionCoeff();
+    } catch (NotImplementedError& err) {
+        GTEST_SKIP() << err.getMethod() << " threw NotImplementedError";
+    }
+
+    double P = phase->pressure();
+    double T1 = phase->temperature();
+    double mv1 = phase->molarVolume();
+
+    double T2 = T1 * (1 + 1e-7);
+    phase->setState_TP(T2, P);
+    double alphaV2 = phase->thermalExpansionCoeff();
+    double mv2 = phase->molarVolume();
+
+    double alphaV_mid = 0.5 * (alphaV1 + alphaV2);
+    double mv_mid = 0.5 * (mv1 + mv2);
+    double alphaV_fd = 1 / mv_mid * (mv2 - mv1) / (T2 - T1);
+
+    EXPECT_NEAR(alphaV_fd, alphaV_mid, max({rtol_fd * alphaV_mid, rtol_fd * alphaV_fd, atol}));
+}
+
 // ---------- Tests for consistency of standard state properties ---------------
 
 TEST_P(TestConsistency, hk0_eq_uk0_plus_p_vk0)

From b02fa1a820e7a9e4ff9cc705db14b02be54da9c2 Mon Sep 17 00:00:00 2001
From: corykinney <cory.kinney@ucf.edu>
Date: Fri, 13 Jan 2023 13:08:38 -0500
Subject: [PATCH 16/17] [Test] Add SRK example input file and set up thermo
 consistency tests

---
 test/data/co2_SRK_example.yaml   | 129 +++++++++++++++++++++++++++++++
 test/data/consistency-cases.yaml |   9 +++
 test/thermo/consistency.cpp      |   6 ++
 3 files changed, 144 insertions(+)
 create mode 100644 test/data/co2_SRK_example.yaml

diff --git a/test/data/co2_SRK_example.yaml b/test/data/co2_SRK_example.yaml
new file mode 100644
index 0000000000..9a1ce6ee34
--- /dev/null
+++ b/test/data/co2_SRK_example.yaml
@@ -0,0 +1,129 @@
+units: {length: cm, quantity: mol, activation-energy: cal/mol}
+
+phases:
+- name: CO2-PR
+  species: [CO2, H2O, H2, CO, CH4, O2, N2]
+  thermo: Soave-Redlich-Kwong
+  kinetics: gas
+  reactions: all
+  state: {T: 300, P: 1 atm, mole-fractions: {CO2: 0.99, H2: 0.01}}
+
+
+species:
+- name: CO2
+  composition: {C: 1, O: 2}
+  thermo:
+    model: NASA7
+    temperature-ranges: [200.0, 1000.0, 3500.0]
+    data:
+    - [2.35677352, 8.98459677e-03, -7.12356269e-06, 2.45919022e-09, -1.43699548e-13,
+      -4.83719697e+04, 9.90105222]
+    - [3.85746029, 4.41437026e-03, -2.21481404e-06, 5.23490188e-10, -4.72084164e-14,
+      -4.8759166e+04, 2.27163806]
+    note: L7/88
+  equation-of-state:
+    model: Soave-Redlich-Kwong
+    a: 3.700552E+11
+    b: 29.6547
+    acentric-factor: 0.228
+- name: H2O
+  composition: {H: 2, O: 1}
+  thermo:
+    model: NASA7
+    temperature-ranges: [200.0, 1000.0, 3500.0]
+    data:
+    - [4.19864056, -2.0364341e-03, 6.52040211e-06, -5.48797062e-09, 1.77197817e-12,
+      -3.02937267e+04, -0.849032208]
+    - [3.03399249, 2.17691804e-03, -1.64072518e-07, -9.7041987e-11, 1.68200992e-14,
+      -3.00042971e+04, 4.9667701]
+    note: L8/89
+  equation-of-state:
+    model: Soave-Redlich-Kwong
+    a: 5.608487E+11
+    b: 21.1282
+    acentric-factor: 0.344
+- name: H2
+  composition: {H: 2}
+  thermo:
+    model: NASA7
+    temperature-ranges: [200.0, 1000.0, 3500.0]
+    data:
+    - [2.34433112, 7.98052075e-03, -1.9478151e-05, 2.01572094e-08, -7.37611761e-12,
+      -917.935173, 0.683010238]
+    - [3.3372792, -4.94024731e-05, 4.99456778e-07, -1.79566394e-10, 2.00255376e-14,
+      -950.158922, -3.20502331]
+    note: TPIS78
+  equation-of-state:
+    model: Soave-Redlich-Kwong
+    a: 2.494771E+10
+    b: 18.4290
+    acentric-factor: -0.22
+- name: CO
+  composition: {C: 1, O: 1}
+  thermo:
+    model: NASA7
+    temperature-ranges: [200.0, 1000.0, 3500.0]
+    data:
+    - [3.57953347, -6.1035368e-04, 1.01681433e-06, 9.07005884e-10, -9.04424499e-13,
+      -1.4344086e+04, 3.50840928]
+    - [2.71518561, 2.06252743e-03, -9.98825771e-07, 2.30053008e-10, -2.03647716e-14,
+      -1.41518724e+04, 7.81868772]
+    note: TPIS79
+  equation-of-state:
+    model: Soave-Redlich-Kwong
+    a: 1.502575E+11
+    b: 27.4578
+    acentric-factor: 0.049
+- name: CH4
+  composition: {C: 1, H: 4}
+  thermo:
+    model: NASA7
+    temperature-ranges: [200.0, 1000.0, 3500.0]
+    data:
+    - [5.14987613, -0.0136709788, 4.91800599e-05, -4.84743026e-08, 1.66693956e-11,
+      -1.02466476e+04, -4.64130376]
+    - [0.074851495, 0.0133909467, -5.73285809e-06, 1.22292535e-09, -1.0181523e-13,
+      -9468.34459, 18.437318]
+    note: L8/88
+  equation-of-state:
+    model: Soave-Redlich-Kwong
+    a: 2.333891E+11
+    b: 29.8499
+    acentric-factor: 0.01
+- name: O2
+  composition: {O: 2}
+  thermo:
+    model: NASA7
+    temperature-ranges: [200.0, 1000.0, 3500.0]
+    data:
+    - [3.78245636, -2.99673416e-03, 9.84730201e-06, -9.68129509e-09, 3.24372837e-12,
+      -1063.94356, 3.65767573]
+    - [3.28253784, 1.48308754e-03, -7.57966669e-07, 2.09470555e-10, -2.16717794e-14,
+      -1088.45772, 5.45323129]
+    note: TPIS89
+  equation-of-state:
+    model: Soave-Redlich-Kwong
+    a: 1.400265E+11
+    b: 22.0823
+    acentric-factor: 0.022
+- name: N2
+  composition: {N: 2}
+  thermo:
+    model: NASA7
+    temperature-ranges: [300.0, 1000.0, 5000.0]
+    data:
+    - [3.298677, 1.4082404e-03, -3.963222e-06, 5.641515e-09, -2.444854e-12,
+      -1020.8999, 3.950372]
+    - [2.92664, 1.4879768e-03, -5.68476e-07, 1.0097038e-10, -6.753351e-15,
+      -922.7977, 5.980528]
+    note: '121286'
+  equation-of-state:
+    model: Soave-Redlich-Kwong
+    a: 1.388390E+11
+    b: 31.2734
+    acentric-factor: 0.04
+
+
+reactions:
+- equation: CO2 + H2 <=> CO + H2O  # Reaction 1
+  rate-constant: {A: 1.2E+3, b: 0, Ea: 0}
diff --git a/test/data/consistency-cases.yaml b/test/data/consistency-cases.yaml
index dd06e62482..2e2761a16c 100644
--- a/test/data/consistency-cases.yaml
+++ b/test/data/consistency-cases.yaml
@@ -55,6 +55,15 @@ peng-robinson:
   - {T: 320, P: 200 bar, X: {CO2: 0.3, CH4: 0.5, H2O: 0.2}}
   - {T: 600, P: 200 bar, X: {CO2: 0.4, CH4: 0.2, H2O: 0.4}}
 
+soave-redlich-kwong:
+  setup:
+    file: co2_SRK_example.yaml
+    rtol_fd: 1e-5
+  states:
+  - {T: 300, P: 101325, X: {CO2: 0.7, CH4: 0.2, H2O: 0.1}}
+  - {T: 320, P: 200 bar, X: {CO2: 0.3, CH4: 0.5, H2O: 0.2}}
+  - {T: 600, P: 200 bar, X: {CO2: 0.4, CH4: 0.2, H2O: 0.4}}
+
 ideal-molal-solution:
   setup:
     file: thermo-models.yaml
diff --git a/test/thermo/consistency.cpp b/test/thermo/consistency.cpp
index 78a8d2602a..59e0c0fc93 100644
--- a/test/thermo/consistency.cpp
+++ b/test/thermo/consistency.cpp
@@ -606,6 +606,12 @@ INSTANTIATE_TEST_SUITE_P(PengRobinson, TestConsistency,
         testing::ValuesIn(getStates("peng-robinson")))
 );
 
+INSTANTIATE_TEST_SUITE_P(SoaveRedlichKwong, TestConsistency,
+    testing::Combine(
+        testing::Values(getSetup("soave-redlich-kwong")),
+        testing::ValuesIn(getStates("soave-redlich-kwong")))
+);
+
 INSTANTIATE_TEST_SUITE_P(IdealMolalSolution, TestConsistency,
     testing::Combine(
         testing::Values(getSetup("ideal-molal-solution")),

From d9e783a8daf9b5d98d4f2ea337c02621f0ed3d2c Mon Sep 17 00:00:00 2001
From: Cory Kinney <cory.kinney@ucf.edu>
Date: Thu, 27 Apr 2023 16:53:53 -0400
Subject: [PATCH 17/17] Implement simplified compressibility properties

---
 src/thermo/SoaveRedlichKwong.cpp | 33 ++++++--------------------------
 1 file changed, 6 insertions(+), 27 deletions(-)

diff --git a/src/thermo/SoaveRedlichKwong.cpp b/src/thermo/SoaveRedlichKwong.cpp
index 643d48a2f4..62b2c43e6b 100644
--- a/src/thermo/SoaveRedlichKwong.cpp
+++ b/src/thermo/SoaveRedlichKwong.cpp
@@ -481,37 +481,16 @@ double SoaveRedlichKwong::dpdVCalc(double T, double molarVol, double& presCalc)
     return -GasConstant * T / (vmb * vmb) + m_aAlpha_mix * (2 * molarVol + m_b) / (denom * denom);
 }
 
-double SoaveRedlichKwong::isothermalCompressibility() const 
+double SoaveRedlichKwong::isothermalCompressibility() const
 {
-    double P = pressure();
-    double RT_ = RT();
-    double Z = z();
-
-    double A = m_aAlpha_mix * P / (RT_ * RT_);
-    double B = m_b * P / RT_;
-
-    double dAdP = A / P;
-    double dBdP = B / P;
-    double dZdP = ((B - Z) * dAdP + (A + Z + 2 * B * Z) * dBdP) / (3 * Z * Z - 2 * Z + A - B - B * B);
-
-    return 1 / P - dZdP / Z;
+    calculatePressureDerivatives();
+    return -1 / (molarVolume() * m_dpdV);
 }
 
-double SoaveRedlichKwong::thermalExpansionCoeff() const 
+double SoaveRedlichKwong::thermalExpansionCoeff() const
 {
-    double P = pressure();
-    double T = temperature();
-    double RT_ = RT();
-    double Z = z();
-
-    double A = m_aAlpha_mix * P / (RT_ * RT_);
-    double B = m_b * P / RT_;
-
-    double dAdT = P / (RT_ * RT_) * (daAlpha_dT() - 2 * m_aAlpha_mix / T);
-    double dBdT = -B / T;
-    double dZdT = ((B - Z) * dAdT + (A + Z + 2 * B * Z) * dBdT) / (3 * Z * Z - 2 * Z + A - B - B * B);
-
-    return 1 / T + dZdT / Z;
+    calculatePressureDerivatives();
+    return -m_dpdT / (molarVolume() * m_dpdV);
 }
 
 void SoaveRedlichKwong::calculatePressureDerivatives() const