From cfa3b3a3b7ed670e4ce7421f7aeb18185f2202c3 Mon Sep 17 00:00:00 2001
From: Ian Bell <ian.bell@nist.gov>
Date: Tue, 24 Jan 2023 20:25:29 -0500
Subject: [PATCH] Add L-J EOS from Johnson
And some tests
---
include/teqp/json_builder.hpp | 3 +
include/teqp/models/fwd.hpp | 4 +-
include/teqp/models/mie/lennardjones.hpp | 147 +++++++++++++++++++++++
src/tests/catch_test_LJmodels.cxx | 61 ++++++++++
4 files changed, 214 insertions(+), 1 deletion(-)
create mode 100644 src/tests/catch_test_LJmodels.cxx
diff --git a/include/teqp/json_builder.hpp b/include/teqp/json_builder.hpp
index 4e4487f..68b2f11 100644
--- a/include/teqp/json_builder.hpp
+++ b/include/teqp/json_builder.hpp
@@ -104,6 +104,9 @@ namespace teqp {
else if (kind == "LJ126_KolafaNezbeda1994"){
return LJ126KolafaNezbeda1994();
}
+ else if (kind == "LJ126_Johnson1993"){
+ return LJ126Johnson1993();
+ }
else if (kind == "2CLJF-Dipole"){
return twocenterljf::build_two_center_model_dipole(spec.at("author"), spec.at("L^*"), spec.at("(mu^*)^2"));
}
diff --git a/include/teqp/models/fwd.hpp b/include/teqp/models/fwd.hpp
index 4937934..a0e4297 100644
--- a/include/teqp/models/fwd.hpp
+++ b/include/teqp/models/fwd.hpp
@@ -41,6 +41,7 @@ namespace teqp {
using vdWEOS_t = vdWEOS<double>;
using twocenterLJF_t = decltype(twocenterljf::build_two_center_model_dipole(std::string{}, double{}, double{}));
using LJ126KolafaNezbeda1994_t = LJ126KolafaNezbeda1994;
+ using LJ126Johnson1993_t = LJ126Johnson1993;
using idealgas_t = IdealHelmholtz;
@@ -58,6 +59,7 @@ namespace teqp {
SW_EspindolaHeredia2009_t,
EXP6_Kataoka1992_t,
twocenterLJF_t,
- LJ126KolafaNezbeda1994_t
+ LJ126KolafaNezbeda1994_t,
+ LJ126Johnson1993_t
>;
}
diff --git a/include/teqp/models/mie/lennardjones.hpp b/include/teqp/models/mie/lennardjones.hpp
index 53d6434..7df7055 100644
--- a/include/teqp/models/mie/lennardjones.hpp
+++ b/include/teqp/models/mie/lennardjones.hpp
@@ -1,6 +1,7 @@
#pragma once
#include "teqp/models/multifluid.hpp"
+#include "teqp/types.hpp"
namespace teqp {
@@ -203,4 +204,150 @@ namespace teqp {
}
};
+ /**
+ J. KARL JOHNSON, JOHN A. ZOLLWEG and KEITH E. GUBBINS
+ The Lennard-Jones equation of state revisited
+ MOLECULAR PHYSICS,1993, VOL. 78, No. 3, 591-618
+ doi: 10.1080/00268979300100411
+ */
+ class LJ126Johnson1993 {
+
+ private:
+ template<typename T>
+ auto POW2(const T& x) const -> T{
+ return x*x;
+ }
+
+ template<typename T>
+ auto POW3(const T& x) const -> T{
+ return POW2(x)*x;
+ }
+
+ template<typename T>
+ auto POW4(const T& x) const -> T{
+ return POW2(x)*POW2(x);
+ }
+
+ const double gamma = 3.0;
+ const std::vector<double> x {
+ 0.0, // placeholder for i=0 term since C++ uses 0-based indexing
+ 0.8623085097507421,
+ 2.976218765822098,
+ -8.402230115796038,
+ 0.1054136629203555,
+ -0.8564583828174598,
+ 1.582759470107601,
+ 0.7639421948305453,
+ 1.753173414312048,
+ 2.798291772190376e03,
+ -4.8394220260857657e-2,
+ 0.9963265197721935,
+ -3.698000291272493e01,
+ 2.084012299434647e01,
+ 8.305402124717285e01,
+ -9.574799715203068e02,
+ -1.477746229234994e02,
+
+ 6.398607852471505e01,
+ 1.603993673294834e01,
+ 6.805916615864377e01,
+ -2.791293578795945e03,
+ -6.245128304568454,
+ -8.116836104958410e03,
+ 1.488735559561229e01,
+ -1.059346754655084e04,
+ -1.131607632802822e02,
+ -8.867771540418822e03,
+ -3.986982844450543e01,
+ -4.689270299917261e03,
+ 2.593535277438717e02,
+ -2.694523589434903e03,
+ -7.218487631550215e02,
+ 1.721802063863269e02
+ };
+
+ // Table 5
+ template<typename TTYPE>
+ auto get_ai(const int i, const TTYPE& Tstar) const -> TTYPE{
+ switch(i){
+ case 1:
+ return x[1]*Tstar + x[2]*sqrt(Tstar) + x[3] + x[4]/Tstar + x[5]/POW2(Tstar);
+ case 2:
+ return x[6]*Tstar + x[7] + x[8]/Tstar + x[9]/POW2(Tstar);
+ case 3:
+ return x[10]*Tstar + x[11] + x[12]/Tstar;
+ case 4:
+ return x[13];
+ case 5:
+ return x[14]/Tstar + x[15]/POW2(Tstar);
+ case 6:
+ return x[16]/Tstar;
+ case 7:
+ return x[17]/Tstar + x[18]/POW2(Tstar);
+ case 8:
+ return x[19]/POW2(Tstar);
+ default:
+ throw teqp::InvalidArgument("i is not possible in get_ai");
+ }
+ }
+
+ // Table 6
+ template<typename TTYPE>
+ auto get_bi(const int i, const TTYPE& Tstar) const -> TTYPE{
+ switch(i){
+ case 1:
+ return x[20]/POW2(Tstar) + x[21]/POW3(Tstar);
+ case 2:
+ return x[22]/POW2(Tstar) + x[23]/POW4(Tstar);
+ case 3:
+ return x[24]/POW2(Tstar) + x[25]/POW3(Tstar);
+ case 4:
+ return x[26]/POW2(Tstar) + x[27]/POW4(Tstar);
+ case 5:
+ return x[28]/POW2(Tstar) + x[29]/POW3(Tstar);
+ case 6:
+ return x[30]/POW2(Tstar) + x[31]/POW3(Tstar) + x[32]/POW4(Tstar);
+ default:
+ throw teqp::InvalidArgument("i is not possible in get_bi");
+ }
+ }
+
+ // Table 7
+ template<typename FTYPE, typename RHOTYPE>
+ auto get_Gi(const int i, const FTYPE& F, const RHOTYPE& rhostar) const -> RHOTYPE{
+ if (i == 1){
+ return forceeval((1.0-F)/(2.0*gamma));
+ }
+ else{
+ // Recursive definition of the other terms;
+ return forceeval(-(F*powi(rhostar, 2*(i-1)) - 2.0*(i-1)*get_Gi(i-1, F, rhostar))/(2.0*gamma));
+ }
+ }
+
+ template<typename TTYPE, typename RHOTYPE>
+ auto get_alphar(const TTYPE& Tstar, const RHOTYPE& rhostar) const{
+ std::common_type_t<TTYPE, RHOTYPE> summer = 0.0;
+ auto F = exp(-gamma*POW2(rhostar));
+ for (int i = 1; i <= 8; ++i){
+ summer += get_ai(i, Tstar)*powi(rhostar, i)/static_cast<double>(i);
+ }
+ for (int i = 1; i <= 6; ++i){
+ summer += get_bi(i, Tstar)*get_Gi(i, F, rhostar);
+ }
+ return forceeval(summer);
+ }
+
+ public:
+ // We are in "simulation units", so R is 1.0, and T and rho that
+ // go into alphar are actually T^* and rho^*
+ template<typename MoleFracType>
+ double R(const MoleFracType &) const { return 1.0; }
+
+ template<typename TTYPE, typename RHOTYPE, typename MoleFracType>
+ auto alphar(const TTYPE& Tstar, const RHOTYPE& rhostar, const MoleFracType& /*molefrac*/) const {
+ return forceeval(get_alphar(Tstar, rhostar)/Tstar);
+ }
+
+ };
+
};
diff --git a/src/tests/catch_test_LJmodels.cxx b/src/tests/catch_test_LJmodels.cxx
new file mode 100644
index 0000000..30f4caa
--- /dev/null
+++ b/src/tests/catch_test_LJmodels.cxx
@@ -0,0 +1,61 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/catch_approx.hpp>
+using Catch::Approx;
+
+#include "teqp/models/mie/lennardjones.hpp"
+#include "teqp/algorithms/critical_pure.hpp"
+#include "teqp/derivs.hpp"
+
+using namespace teqp;
+
+TEST_CASE("Test for critical point of Kolafa-Nezbeda", "[LJ126]")
+{
+ auto m = LJ126KolafaNezbeda1994();
+ auto soln = solve_pure_critical(m, 1.3, 0.3);
+ auto Tc = 1.3396, rhoc = 0.3108;
+ CHECK(std::get<0>(soln) == Approx(Tc).margin(0.001));
+ CHECK(std::get<1>(soln) == Approx(rhoc).margin(0.001));
+}
+TEST_CASE("Test for critical point of Thol", "[LJ126-Thol]")
+{
+ auto m = build_LJ126_TholJPCRD2016();
+ auto soln = solve_pure_critical(m, 1.3, 0.3);
+ auto Tc = 1.32, rhoc = 0.31;
+ CHECK(std::get<0>(soln) == Approx(Tc).margin(0.01));
+ CHECK(std::get<1>(soln) == Approx(rhoc).margin(0.01));
+}
+TEST_CASE("Test for critical point of Johnson", "[LJ126]")
+{
+ auto m = LJ126Johnson1993();
+ auto soln = solve_pure_critical(m, 1.3, 0.3);
+ auto Tc = 1.313, rhoc = 0.310;
+ CHECK(std::get<0>(soln) == Approx(Tc).margin(0.001));
+ CHECK(std::get<1>(soln) == Approx(rhoc).margin(0.001));
+}
+TEST_CASE("Test single point values for Johnson against calculated values from S. Stephan", "[LJ126]")
+{
+ auto m = LJ126Johnson1993();
+ auto z = (Eigen::ArrayXd(1) << 1.0).finished();
+
+ auto Bn = VirialDerivatives<decltype(m)>::get_B2vir(m, 0.8, z);
+ CHECK(Bn == Approx(-7.821026827));
+ auto Bnmcx = VirialDerivatives<decltype(m)>::get_Bnvir<2,ADBackends::multicomplex>(m, 0.8, z)[2];
+ CHECK(Bnmcx == Approx(-7.821026827));
+ auto Bnad = VirialDerivatives<decltype(m)>::get_Bnvir<2,ADBackends::autodiff>(m, 0.8, z)[2];
+ CHECK(Bnad == Approx(-7.821026827));
+
+ auto ar = m.alphar(1.350, 0.600, z);
+ CHECK(ar == Approx(-1.237403479));
+}
+TEST_CASE("Test single point values for K-N against calculated values from S. Stephan", "[LJ126]")
+{
+ auto m = LJ126KolafaNezbeda1994();
+ auto z = (Eigen::ArrayXd(1) << 1.0).finished();
+
+ auto Bn = VirialDerivatives<decltype(m)>::get_B2vir(m, 0.8, z);
+ CHECK(Bn == Approx(-7.821026827).margin(0.0005));
+ auto Bnmcx = VirialDerivatives<decltype(m)>::get_Bnvir<2,ADBackends::multicomplex>(m, 0.8, z)[2];
+ CHECK(Bnmcx == Approx(-7.821026827).margin(0.0005));
+ auto Bnad = VirialDerivatives<decltype(m)>::get_Bnvir<2,ADBackends::autodiff>(m, 0.8, z)[2];
+ CHECK(Bnad == Approx(-7.821026827).margin(0.0005));
+}
--
GitLab