From d8a1337933e2a013e16fb82db3b2753129304036 Mon Sep 17 00:00:00 2001
From: Ray Speth <speth@mit.edu>
Date: Thu, 17 Jan 2019 17:04:14 -0500
Subject: [PATCH] [Input] Create RedlichKister objects from YAML definitions

---
 src/thermo/RedlichKisterVPSSTP.cpp      | 10 +++++++++
 test/data/RedlichKisterVPSSTP_valid.xml |  7 +++++++
 test/data/thermo-models.yaml            | 28 +++++++++++++++++++++++++
 test/thermo/thermoFromYaml.cpp          | 22 +++++++++++++++++++
 4 files changed, 67 insertions(+)

diff --git a/src/thermo/RedlichKisterVPSSTP.cpp b/src/thermo/RedlichKisterVPSSTP.cpp
index 137048a5d..f368f3c78 100644
--- a/src/thermo/RedlichKisterVPSSTP.cpp
+++ b/src/thermo/RedlichKisterVPSSTP.cpp
@@ -179,6 +179,16 @@ void RedlichKisterVPSSTP::getPartialMolarVolumes(doublereal* vbar) const
 
 void RedlichKisterVPSSTP::initThermo()
 {
+    if (m_input.hasKey("interactions")) {
+        for (const auto& item : m_input["interactions"].asVector<AnyMap>()) {
+            auto& species = item["species"].asVector<string>(2);
+            vector_fp h_excess = item.convertVector("excess-enthalpy", "J/kmol");
+            vector_fp s_excess = item.convertVector("excess-entropy", "J/kmol/K");
+            addBinaryInteraction(species[0], species[1],
+                                 h_excess.data(), h_excess.size(),
+                                 s_excess.data(), s_excess.size());
+        }
+    }
     initLengths();
     GibbsExcessVPSSTP::initThermo();
 }
diff --git a/test/data/RedlichKisterVPSSTP_valid.xml b/test/data/RedlichKisterVPSSTP_valid.xml
index 572d1439c..a0617a1bd 100644
--- a/test/data/RedlichKisterVPSSTP_valid.xml
+++ b/test/data/RedlichKisterVPSSTP_valid.xml
@@ -61,6 +61,10 @@
            <cp0 units="J/mol/K">0.0</cp0>
         </const_cp>
       </thermo>
+      <standardState model="constant_incompressible">
+         <molarVolume> 0.036 </molarVolume>
+      </standardState>
+
     </species>
 
     <!-- species V(C6)    -->
@@ -74,6 +78,9 @@
            <cp0 units="J/mol/K">0.0</cp0>
         </const_cp>
       </thermo>
+      <standardState model="constant_incompressible">
+         <molarVolume> 0.036 </molarVolume>
+      </standardState>
     </species>
   </speciesData>
   <reactionData id="reaction_data"/>
diff --git a/test/data/thermo-models.yaml b/test/data/thermo-models.yaml
index 4f60f5c36..48c739edc 100644
--- a/test/data/thermo-models.yaml
+++ b/test/data/thermo-models.yaml
@@ -102,6 +102,16 @@ phases:
   standard-concentration-basis: unity
   species: [Li(l)]
 
+- name: Redlich-Kister-LiC6
+  thermo: Redlich-Kister
+  species: [Li(C6), V(C6)]
+  interactions:
+  - species: [Li(C6), V(C6)]
+    excess-enthalpy: [-3.268e+06, 3.955e+06, -4.573e+06, 6.147e+06, -3.339e+06,
+                      1.117e+07, 2.997e+05, -4.866e+07, 1.362e+05, 1.373e+08,
+                      -2.129e+07, -1.722e+08, 3.956e+07, 9.302e+07, -3.280e+07]
+    excess-entropy: [0.0]
+
 species:
 - name: NaCl(s)
   composition: {Na: 1, Cl: 1}
@@ -169,6 +179,24 @@ species:
     units: {mass: g, length: cm}
     data: [0.536504, -1.04279e-4, 3.84825e-9, -5.2853e-12]
 
+- name: Li(C6)
+  composition: {C: 6, Li: 1}
+  thermo:
+    model: constant-cp
+    h0: -11.65 kJ/mol
+  equation-of-state:
+    model: constant-volume
+    molar-volume: 0.036 m^3/kmol
+
+- name: V(C6)
+  composition: {C: 6}
+  thermo:
+    model: constant-cp
+  equation-of-state:
+    model: constant-volume
+    molar-volume: 0.036 m^3/kmol
+
+
 ideal-molal-fake-species:
 # Fake thermo data (GRI 3.0 coefficients for H2)
 - name: H2O(l)
diff --git a/test/thermo/thermoFromYaml.cpp b/test/thermo/thermoFromYaml.cpp
index bf32c4a03..f7c71b185 100644
--- a/test/thermo/thermoFromYaml.cpp
+++ b/test/thermo/thermoFromYaml.cpp
@@ -238,3 +238,25 @@ TEST(ThermoFromYaml, IdealSolnGas_liquid)
     EXPECT_NEAR(thermo->enthalpy_mass(), 1236522.9439646902, 2e-8);
     EXPECT_NEAR(thermo->entropy_mole(), 49848.48843237689, 2e-8);
 }
+
+TEST(ThermoFromYaml, RedlichKister)
+{
+    AnyMap infile = AnyMap::fromYamlFile("thermo-models.yaml");
+    auto phaseNodes = infile["phases"].asMap("name");
+    auto thermo = newPhase(*phaseNodes.at("Redlich-Kister-LiC6"), infile);
+
+    vector_fp chemPotentials(2);
+    vector_fp dlnActCoeffdx(2);
+    thermo->setState_TP(298.15, OneAtm);
+    thermo->setMoleFractionsByName("Li(C6): 0.6375, V(C6): 0.3625");
+    thermo->getChemPotentials(chemPotentials.data());
+    thermo->getdlnActCoeffdlnX_diag(dlnActCoeffdx.data());
+    EXPECT_NEAR(chemPotentials[0], -1.2618554504124604e+007, 1e-6);
+    EXPECT_NEAR(dlnActCoeffdx[0], 0.200612, 1e-6);
+
+    thermo->setMoleFractionsByName("Li(C6): 0.8625, V(C6): 0.1375");
+    thermo->getChemPotentials(chemPotentials.data());
+    thermo->getdlnActCoeffdlnX_diag(dlnActCoeffdx.data());
+    EXPECT_NEAR(chemPotentials[0], -1.1792994839484975e+07, 1e-6);
+    EXPECT_NEAR(dlnActCoeffdx[0], -0.309379, 1e-6);
+}