diff --git a/include/teqp/cpp/teqpcpp.hpp b/include/teqp/cpp/teqpcpp.hpp
index ab6ddd5a483dde8d0a9d6eb253a533a89e6ce809..7ce5bee4476ad924dbcc4d83ede38cf0d4b4d5bb 100644
--- a/include/teqp/cpp/teqpcpp.hpp
+++ b/include/teqp/cpp/teqpcpp.hpp
@@ -121,7 +121,7 @@ namespace teqp {
virtual EArray33d get_deriv_mat2(const double T, double rho, const EArrayd& z ) const = 0;
- virtual std::tuple<double, double> solve_pure_critical(const double T, const double rho, const std::optional<nlohmann::json>&) const = 0;
+ virtual std::tuple<double, double> solve_pure_critical(const double T, const double rho, const std::optional<nlohmann::json>& = std::nullopt) const = 0;
virtual std::tuple<double, double> extrapolate_from_critical(const double Tc, const double rhoc, const double Tgiven) const = 0;
virtual std::tuple<EArrayd, EMatrixd> get_pure_critical_conditions_Jacobian(const double T, const double rho, int alternative_pure_index=-1, int alternative_length=2) const = 0;
virtual EArray2 pure_VLE_T(const double T, const double rhoL, const double rhoV, int maxiter) const = 0;
diff --git a/include/teqp/json_builder.hpp b/include/teqp/json_builder.hpp
index f7026ec5498311533d0ec7ba01855506569b3d10..4b795d2f76e74adc4379d9ac75cd736d499be4a1 100644
--- a/include/teqp/json_builder.hpp
+++ b/include/teqp/json_builder.hpp
@@ -11,34 +11,6 @@
namespace teqp {
inline AllowedModels build_model(const nlohmann::json& json) {
-
- auto build_square_matrix = [](const nlohmann::json& j){
- if (j.is_null() || (j.is_array() && j.size() == 0)){
- return Eigen::ArrayXXd(0, 0);
- }
- try{
- const std::valarray<std::valarray<double>> m = j;
- // First assume that the matrix is square, resize
- Eigen::ArrayXXd mat(m.size(), m.size());
- if (m.size() == 0){
- return mat;
- }
- // Then copy elements over
- for (auto i = 0; i < m.size(); ++i){
- auto row = m[i];
- if (row.size() != mat.rows()){
- throw std::invalid_argument("provided matrix is not square");
- }
- for (auto j = 0; j < row.size(); ++j){
- mat(i, j) = row[j];
- }
- }
- return mat;
- }
- catch(const nlohmann::json::exception&){
- throw teqp::InvalidArgument("Unable to convert this kmat to a 2x2 matrix of doubles:" + j.dump(2));
- }
- };
// Extract the name of the model and the model parameters
std::string kind = json.at("kind");
@@ -112,48 +84,7 @@ namespace teqp {
}
}
else if (kind == "SAFT-VR-Mie") {
- using namespace SAFTVRMie;
- std::optional<Eigen::ArrayXXd> kmat;
- if (spec.contains("kmat") && spec.at("kmat").is_array() && spec.at("kmat").size() > 0){
- kmat = build_square_matrix(spec["kmat"]);
- }
-
- if (spec.contains("names")){
- std::vector<std::string> names = spec["names"];
- if (kmat && kmat.value().rows() != names.size()){
- throw teqp::InvalidArgument("Provided length of names of " + std::to_string(names.size()) + " does not match the dimension of the kmat of " + std::to_string(kmat.value().rows()));
- }
- return SAFTVRMieMixture(names, kmat);
- }
- else if (spec.contains("coeffs")){
- std::vector<SAFTVRMieCoeffs> coeffs;
- for (auto j : spec["coeffs"]) {
- SAFTVRMieCoeffs c;
- c.name = j.at("name");
- c.m = j.at("m");
- c.sigma_m = (j.contains("sigma_m")) ? j.at("sigma_m").get<double>() : j.at("sigma_Angstrom").get<double>()/1e10;
- c.epsilon_over_k = j.at("epsilon_over_k");
- c.lambda_r = j.at("lambda_r");
- c.lambda_a = j.at("lambda_a");
- c.BibTeXKey = j.at("BibTeXKey");
- if (j.contains("(mu^*)^2") && j.contains("nmu")){
- c.mustar2 = j.at("(mu^*)^2");
- c.nmu = j.at("nmu");
- }
- if (j.contains("(Q^*)^2") && j.contains("nQ")){
- c.Qstar2 = j.at("(Q^*)^2");
- c.nQ = j.at("nQ");
- }
- coeffs.push_back(c);
- }
- if (kmat && kmat.value().rows() != coeffs.size()){
- throw teqp::InvalidArgument("Provided length of coeffs of " + std::to_string(coeffs.size()) + " does not match the dimension of the kmat of " + std::to_string(kmat.value().rows()));
- }
- return SAFTVRMieMixture(coeffs, kmat);
- }
- else{
- throw std::invalid_argument("you must provide names or coeffs, but not both");
- }
+ return SAFTVRMie::SAFTVRMiefactory(spec);
}
else if (kind == "multifluid") {
return multifluidfactory(spec);
diff --git a/include/teqp/json_tools.hpp b/include/teqp/json_tools.hpp
index 2974feab6dc88302bb8fa9a8ae5530785014fa3d..952b3edb9dadc34a8f4e9ab4cd6e0c84cf37f7cd 100644
--- a/include/teqp/json_tools.hpp
+++ b/include/teqp/json_tools.hpp
@@ -4,6 +4,7 @@
#include <set>
#include <filesystem>
#include <fstream>
+#include "teqp/exceptions.hpp"
namespace teqp{
@@ -29,4 +30,32 @@ namespace teqp{
for (auto k : ks) { lengths.insert(j.at(k).size()); }
return lengths.size() == 1;
}
+
+ auto build_square_matrix = [](const nlohmann::json& j){
+ if (j.is_null() || (j.is_array() && j.size() == 0)){
+ return Eigen::ArrayXXd(0, 0);
+ }
+ try{
+ const std::valarray<std::valarray<double>> m = j;
+ // First assume that the matrix is square, resize
+ Eigen::ArrayXXd mat(m.size(), m.size());
+ if (m.size() == 0){
+ return mat;
+ }
+ // Then copy elements over
+ for (auto i = 0; i < m.size(); ++i){
+ auto row = m[i];
+ if (row.size() != mat.rows()){
+ throw std::invalid_argument("provided matrix is not square");
+ }
+ for (auto j = 0; j < row.size(); ++j){
+ mat(i, j) = row[j];
+ }
+ }
+ return mat;
+ }
+ catch(const nlohmann::json::exception&){
+ throw teqp::InvalidArgument("Unable to convert this kmat to a 2x2 matrix of doubles:" + j.dump(2));
+ }
+ };
}
diff --git a/include/teqp/models/saft/polar_terms.hpp b/include/teqp/models/saft/polar_terms.hpp
index c558a7aebf2ac425ea1097d7a07a2c6b0332b4dd..4cd3d7d2909f501eeff5c735831927233346d1a7 100644
--- a/include/teqp/models/saft/polar_terms.hpp
+++ b/include/teqp/models/saft/polar_terms.hpp
@@ -379,7 +379,11 @@ public:
const auto& sigma = sigma_m; // concision
const auto N = mole_fractions.size();
- std::common_type_t<TTYPE, RhoType, decltype(mole_fractions[0])> summer112 = 0.0, summer123 = 0.0, summer224 = 0.0;
+ std::common_type_t<TTYPE, RhoType, decltype(mole_fractions[0])> alpha2_112 = 0.0, alpha2_123 = 0.0, alpha2_224 = 0.0;
+
+ const auto factor_112 = forceeval(-2.0*PI_*rhoN/3.0); //*POW2(4*PI_*epsilon_0)
+ const auto factor_123 = forceeval(-PI_*rhoN/3.0);
+ const auto factor_224 = forceeval(-14.0*PI_*rhoN/5.0);
for (auto i = 0; i < N; ++i){
for (auto j = 0; j < N; ++j){
@@ -391,16 +395,12 @@ public:
auto leading = x[i]*x[j]/(Tstari*Tstarj); // common for all alpha_2 terms
auto Tstarij = forceeval(T/epskij);
- summer112 += leading*POW3(sigma[i]*sigma[j])/POW3(sigmaij)*mubar2[i]*mubar2[j]*J6.get_J(Tstarij, rhostar);
- summer123 += leading*POW3(sigma[i])*POW5(sigma[j])/POW5(sigmaij)*mubar2[i]*Qbar2[j]*J8.get_J(Tstarij, rhostar);
- summer224 += leading*POW5(sigma[i]*sigma[j])/POW7(sigmaij)*Qbar2[i]*Qbar2[j]*J10.get_J(Tstarij, rhostar);
+ alpha2_112 += factor_112*leading*POW3(sigma[i]*sigma[j])/POW3(sigmaij)*mubar2[i]*mubar2[j]*J6.get_J(Tstarij, rhostar);
+ alpha2_123 += factor_123*leading*POW3(sigma[i])*POW5(sigma[j])/POW5(sigmaij)*mubar2[i]*Qbar2[j]*J8.get_J(Tstarij, rhostar);
+ alpha2_224 += factor_224*leading*POW5(sigma[i]*sigma[j])/POW7(sigmaij)*Qbar2[i]*Qbar2[j]*J10.get_J(Tstarij, rhostar);
}
}
- auto alpha2_112 = -2.0*PI_*rhoN*POW2(4*PI_*epsilon_0)/3.0*summer112;
- auto alpha2_123 = -PI_*rhoN*POW2(4*PI_*epsilon_0)/3.0*summer123;
- auto alpha2_224 = -14.0*PI_*rhoN*POW2(4*PI_*epsilon_0)/5.0*summer224;
-
return forceeval(alpha2_112 + 2.0*alpha2_123 + alpha2_224);
}
@@ -440,31 +440,43 @@ public:
auto sigmaik = (sigma[i]+sigma[k])/2;
auto sigmajk = (sigma[j]+sigma[k])/2;
auto leadingijk = x[i]*x[j]*x[k]/(Tstari*Tstarj*Tstark);
+
+ auto get_Kijk = [&](const auto& Kint){
+ return pow(Kint.get_K(Tstarij, rhostar)*Kint.get_K(Tstarik, rhostar)*Kint.get_K(Tstarjk, rhostar), 1.0/3.0);
+ };
- auto K222333 = pow(K222_333.get_K(Tstarij, rhostar)*K222_333.get_K(Tstarik, rhostar)*K222_333.get_K(Tstarjk, rhostar), 1.0/3.0);
- summerB_112_112_112 += leadingijk*POW3(sigma[i]*sigma[j]*sigma[k])/(sigmaij*sigmaik*sigmajk)*mubar2[i]*mubar2[j]*mubar2[k]*K222333;
- auto K233344 = pow(K233_344.get_K(Tstarij, rhostar)*K233_344.get_K(Tstarik, rhostar)*K233_344.get_K(Tstarjk, rhostar), 1.0/3.0);
- summerB_112_123_123 += leadingijk*POW3(sigma[i]*sigma[j])*POW5(sigma[k])/(sigmaij*POW2(sigmaik*sigmajk))*mubar2[i]*mubar2[j]*Qbar2[k]*K233344;
- auto K334445 = pow(K334_445.get_K(Tstarij, rhostar)*K334_445.get_K(Tstarik, rhostar)*K334_445.get_K(Tstarjk, rhostar), 1.0/3.0);
- summerB_123_123_224 += leadingijk*POW3(sigma[i])*POW5(sigma[j]*sigma[k])/(POW2(sigmaij*sigmaik)*POW3(sigmajk))*mubar2[i]*Qbar2[j]*Qbar2[k]*K334445;
- auto K444555 = pow(K444_555.get_K(Tstarij, rhostar)*K444_555.get_K(Tstarik, rhostar)*K444_555.get_K(Tstarjk, rhostar), 1.0/3.0);
- summerB_224_224_224 += leadingijk*POW5(sigma[i]*sigma[j]*sigma[k])/(POW3(sigmaij*sigmaik*sigmajk))*Qbar2[i]*Qbar2[j]*Qbar2[k]*K444555;
+ if (std::abs(mubar2[i]*mubar2[j]*mubar2[k]) > 0){
+ auto K222333 = get_Kijk(K222_333);
+ summerB_112_112_112 += leadingijk*POW3(sigma[i]*sigma[j]*sigma[k])/(sigmaij*sigmaik*sigmajk)*mubar2[i]*mubar2[j]*mubar2[k]*K222333;
+ }
+ if (std::abs(mubar2[i]*mubar2[j]*Qbar2[k]) > 0){
+ auto K233344 = get_Kijk(K233_344);
+ summerB_112_123_123 += leadingijk*POW3(sigma[i]*sigma[j])*POW5(sigma[k])/(sigmaij*POW2(sigmaik*sigmajk))*mubar2[i]*mubar2[j]*Qbar2[k]*K233344;
+ }
+ if (std::abs(mubar2[i]*Qbar2[j]*Qbar2[k]) > 0){
+ auto K334445 = get_Kijk(K334_445);
+ summerB_123_123_224 += leadingijk*POW3(sigma[i])*POW5(sigma[j]*sigma[k])/(POW2(sigmaij*sigmaik)*POW3(sigmajk))*mubar2[i]*Qbar2[j]*Qbar2[k]*K334445;
+ }
+ if (std::abs(Qbar2[i]*Qbar2[j]*Qbar2[k]) > 0){
+ auto K444555 = get_Kijk(K444_555);
+ summerB_224_224_224 += leadingijk*POW5(sigma[i]*sigma[j]*sigma[k])/(POW3(sigmaij*sigmaik*sigmajk))*Qbar2[i]*Qbar2[j]*Qbar2[k]*K444555;
+ }
}
}
}
- auto alpha3A_112_112_224 = 8.0*PI_*rhoN*POW3(4*PI_*epsilon_0)/25.0*summerA_112_112_224;
- auto alpha3A_112_123_213 = 8.0*PI_*rhoN*POW3(4*PI_*epsilon_0)/75.0*summerA_112_123_213;
- auto alpha3A_123_123_224 = 8.0*PI_*rhoN*POW3(4*PI_*epsilon_0)/35.0*summerA_123_123_224;
- auto alpha3A_224_224_224 = 144.0*PI_*rhoN*POW3(4*PI_*epsilon_0)/245.0*summerA_224_224_224;
+ auto alpha3A_112_112_224 = 8.0*PI_*rhoN/25.0*summerA_112_112_224;
+ auto alpha3A_112_123_213 = 8.0*PI_*rhoN/75.0*summerA_112_123_213;
+ auto alpha3A_123_123_224 = 8.0*PI_*rhoN/35.0*summerA_123_123_224;
+ auto alpha3A_224_224_224 = 144.0*PI_*rhoN/245.0*summerA_224_224_224;
- auto alpha3A = 3.0*alpha3A_112_112_224 + 6.0*alpha3A_112_123_213 + 6.0*alpha3A_123_123_224 + alpha3A_224_224_224;
+ auto alpha3A = forceeval(3.0*alpha3A_112_112_224 + 6.0*alpha3A_112_123_213 + 6.0*alpha3A_123_123_224 + alpha3A_224_224_224);
- auto alpha3B_112_112_112 = 32.0*POW3(PI_)*POW2(rhoN)*POW3(4*PI_*epsilon_0)/135.0*sqrt(14*PI_/5.0)*summerB_112_112_112;
- auto alpha3B_112_123_123 = 64.0*POW3(PI_)*POW2(rhoN)*POW3(4*PI_*epsilon_0)/315.0*sqrt(3.0*PI_)*summerB_112_123_123;
- auto alpha3B_123_123_224 = -32.0*POW3(PI_)*POW2(rhoN)*POW3(4*PI_*epsilon_0)/45.0*sqrt(22.0*PI_/63.0)*summerB_123_123_224;
- auto alpha3B_224_224_224 = 32.0*POW3(PI_)*POW2(rhoN)*POW3(4*PI_*epsilon_0)/2025.0*sqrt(2002.0*PI_)*summerB_224_224_224;
+ auto alpha3B_112_112_112 = 32.0*POW3(PI_)*POW2(rhoN)/135.0*sqrt(14*PI_/5.0)*summerB_112_112_112;
+ auto alpha3B_112_123_123 = 64.0*POW3(PI_)*POW2(rhoN)/315.0*sqrt(3.0*PI_)*summerB_112_123_123;
+ auto alpha3B_123_123_224 = -32.0*POW3(PI_)*POW2(rhoN)/45.0*sqrt(22.0*PI_/63.0)*summerB_123_123_224;
+ auto alpha3B_224_224_224 = 32.0*POW3(PI_)*POW2(rhoN)/2025.0*sqrt(2002.0*PI_)*summerB_224_224_224;
- auto alpha3B = alpha3B_112_112_112 + 3.0*alpha3B_112_123_123 + 3.0*alpha3B_123_123_224 + alpha3B_224_224_224;
+ auto alpha3B = forceeval(alpha3B_112_112_112 + 3.0*alpha3B_112_123_123 + 3.0*alpha3B_123_123_224 + alpha3B_224_224_224);
return forceeval(alpha3A + alpha3B);
}
@@ -503,7 +515,7 @@ public:
}
};
-/// The variant including the multipolar terms that can be provided
+/// The variant containing the multipolar types that can be provided
using multipolar_contributions_variant = std::variant<
MultipolarContributionGrossVrabec,
MultipolarContributionGubbinsTwu<LuckasJIntegral, LuckasKIntegral>,
diff --git a/include/teqp/models/saftvrmie.hpp b/include/teqp/models/saftvrmie.hpp
index aa9beb5c40788bda5f67b9fcf0aa19ef3d69e167..d77c6ca61548601cfec79944c973219b483c2e48 100644
--- a/include/teqp/models/saftvrmie.hpp
+++ b/include/teqp/models/saftvrmie.hpp
@@ -8,6 +8,7 @@
#include "nlohmann/json.hpp"
#include "teqp/types.hpp"
+#include "teqp/json_tools.hpp"
#include "teqp/exceptions.hpp"
#include "teqp/constants.hpp"
#include "teqp/math/quadrature.hpp"
@@ -708,9 +709,9 @@ private:
std::vector<std::string> names;
const SAFTVRMieChainContributionTerms terms;
- std::optional<SAFTpolar::multipolar_contributions_variant> polar; // Can be present or not
+ const std::optional<SAFTpolar::multipolar_contributions_variant> polar; // Can be present or not
- void check_kmat(const Eigen::ArrayXXd& kmat, std::size_t N) {
+ static void check_kmat(const Eigen::ArrayXXd& kmat, std::size_t N) {
if (kmat.size() == 0){
return;
}
@@ -721,11 +722,12 @@ private:
throw teqp::InvalidArgument("kmat needs to be a square matrix the same size as the number of components");
}
};
- auto get_coeffs_from_names(const std::vector<std::string> &names){
+ static auto get_coeffs_from_names(const std::vector<std::string> &names){
SAFTVRMieLibrary library;
return library.get_coeffs(names);
}
- auto build_chain(const std::vector<SAFTVRMieCoeffs> &coeffs, const std::optional<Eigen::ArrayXXd>& kmat){
+public:
+ static auto build_chain(const std::vector<SAFTVRMieCoeffs> &coeffs, const std::optional<Eigen::ArrayXXd>& kmat){
if (kmat){
check_kmat(kmat.value(), coeffs.size());
}
@@ -748,7 +750,6 @@ private:
return SAFTVRMieChainContributionTerms(m, epsilon_over_k, sigma_m, lambda_r, lambda_a, std::move(mat));
}
}
-
auto build_polar(const std::vector<SAFTVRMieCoeffs> &coeffs) -> decltype(this->polar){
Eigen::ArrayXd mustar2(coeffs.size()), nmu(coeffs.size()), Qstar2(coeffs.size()), nQ(coeffs.size());
auto i = 0;
@@ -773,10 +774,18 @@ private:
public:
SAFTVRMieMixture(const std::vector<std::string> &names, const std::optional<Eigen::ArrayXXd>& kmat = std::nullopt) : SAFTVRMieMixture(get_coeffs_from_names(names), kmat){};
SAFTVRMieMixture(const std::vector<SAFTVRMieCoeffs> &coeffs, const std::optional<Eigen::ArrayXXd> &kmat = std::nullopt) : terms(build_chain(coeffs, kmat)), polar(build_polar(coeffs)) {};
+ SAFTVRMieMixture(SAFTVRMieChainContributionTerms&& terms, std::optional<SAFTpolar::multipolar_contributions_variant> &&polar = std::nullopt) : terms(std::move(terms)), polar(std::move(polar)) {};
+
// PCSAFTMixture( const PCSAFTMixture& ) = delete; // non construction-copyable
SAFTVRMieMixture& operator=( const SAFTVRMieMixture& ) = delete; // non copyable
+ auto chain_factory(const std::vector<SAFTVRMieCoeffs> &coeffs, const std::optional<Eigen::ArrayXXd>& kmat){
+ SAFTVRMieMixture::build_chain(coeffs, kmat);
+ }
+
+ const auto& get_polar(){ return polar; }
+
const auto& get_terms() const { return terms; }
auto get_core_calcs(double T, double rhomolar, const Eigen::ArrayXd& mole_fractions) const {
auto val = terms.get_core_calcs(T, rhomolar, mole_fractions);
@@ -845,7 +854,9 @@ public:
auto visitor = [&](const auto& contrib){
if constexpr(std::decay_t<decltype(contrib)>::arg_spec == mas::TK_rhoNA3_packingfraction_molefractions){
auto rho_A3 = forceeval(rhomolar*N_A*1e-30);
- alphar += contrib.eval(T, rho_A3, vals.zeta[3], mole_fractions).alpha;
+ auto packing_fraction = vals.zeta[3];
+ auto alpha = contrib.eval(T, rho_A3, packing_fraction, mole_fractions).alpha;
+ alphar += alpha;
}
else if constexpr(std::decay_t<decltype(contrib)>::arg_spec == mas::TK_rhoNm3_molefractions){
auto rhoN_m3 = forceeval(rhomolar*N_A);
@@ -861,6 +872,145 @@ public:
return forceeval(alphar);
}
};
+
+static auto SAFTVRMiefactory(const nlohmann::json & spec){
+
+ std::optional<Eigen::ArrayXXd> kmat;
+ if (spec.contains("kmat") && spec.at("kmat").is_array() && spec.at("kmat").size() > 0){
+ kmat = build_square_matrix(spec["kmat"]);
+ }
+
+ if (spec.contains("names")){
+ std::vector<std::string> names = spec["names"];
+ if (kmat && kmat.value().rows() != names.size()){
+ throw teqp::InvalidArgument("Provided length of names of " + std::to_string(names.size()) + " does not match the dimension of the kmat of " + std::to_string(kmat.value().rows()));
+ }
+ return SAFTVRMieMixture(names, kmat);
+ }
+ else if (spec.contains("coeffs")){
+ bool something_polar = false;
+ std::vector<SAFTVRMieCoeffs> coeffs;
+ for (auto j : spec["coeffs"]) {
+ SAFTVRMieCoeffs c;
+ c.name = j.at("name");
+ c.m = j.at("m");
+ c.sigma_m = (j.contains("sigma_m")) ? j.at("sigma_m").get<double>() : j.at("sigma_Angstrom").get<double>()/1e10;
+ c.epsilon_over_k = j.at("epsilon_over_k");
+ c.lambda_r = j.at("lambda_r");
+ c.lambda_a = j.at("lambda_a");
+ c.BibTeXKey = j.at("BibTeXKey");
+
+ // These are legacy definitions of the polar moments
+ if (j.contains("(mu^*)^2") && j.contains("nmu")){
+ c.mustar2 = j.at("(mu^*)^2");
+ c.nmu = j.at("nmu");
+ something_polar = true;
+ }
+ if (j.contains("(Q^*)^2") && j.contains("nQ")){
+ c.Qstar2 = j.at("(Q^*)^2");
+ c.nQ = j.at("nQ");
+ something_polar = true;
+ }
+ if (j.contains("Q_Cm2") || j.contains("Q_DA") || j.contains("mu_Cm") || j.contains("mu_D")){
+ something_polar = true;
+ }
+ coeffs.push_back(c);
+ }
+ if (kmat && kmat.value().rows() != coeffs.size()){
+ throw teqp::InvalidArgument("Provided length of coeffs of " + std::to_string(coeffs.size()) + " does not match the dimension of the kmat of " + std::to_string(kmat.value().rows()));
+ }
+
+ if (!something_polar){
+ // Nonpolar, just m, epsilon, sigma and possibly a kmat matrix with kij coefficients
+ return SAFTVRMieMixture(SAFTVRMieMixture::build_chain(coeffs, kmat));
+ }
+ else{
+ // Polar term is also provided, along with the chain terms
+ std::string polar_model = "GrossVrabec"; // This is the default, as it was the first one implemented
+ if (spec.contains("polar_model")){
+ polar_model = spec["polar_model"];
+ }
+
+ // Go back and extract the dipolar and quadrupolar terms from
+ // the JSON, in base SI units
+ const double D_to_Cm = 3.33564e-30; // C m/D
+ const double mustar2factor = 1.0/(4*static_cast<double>(EIGEN_PI)*8.8541878128e-12*1.380649e-23); // factor=1/(4*pi*epsilon_0*k_B), such that (mu^*)^2 := factor*mu[Cm]^2/((epsilon/kB)[K]*sigma[m]^3)
+ const double Qstar2factor = 1.0/(4*static_cast<double>(EIGEN_PI)*8.8541878128e-12*1.380649e-23); // same as mustar2factor
+ auto N = coeffs.size();
+ Eigen::ArrayXd ms(N), epsks(N), sigma_ms(N), mu_Cm(N), Q_Cm2(N), nQ(N), nmu(N);
+ Eigen::Index i = 0;
+ for (auto j : spec["coeffs"]) {
+ double m = j.at("m");
+ double sigma_m = (j.contains("sigma_m")) ? j.at("sigma_m").get<double>() : j.at("sigma_Angstrom").get<double>()/1e10;
+ double epsilon_over_k = j.at("epsilon_over_k");
+ auto get_dipole_Cm = [&]() -> double {
+ if (j.contains("(mu^*)^2") && j.contains("nmu")){
+ // Terms defined like in Gross&Vrabec; backwards-compatibility
+ double mustar2 = j.at("(mu^*)^2");
+ return sqrt(mustar2*(m*epsilon_over_k*pow(sigma_m, 3))/mustar2factor);
+ }
+ else if (j.contains("mu_Cm")){
+ return j.at("mu_Cm");
+ }
+ else if (j.contains("mu_D")){
+ return j.at("mu_D").get<double>()*D_to_Cm;
+ }
+ else{
+ return 0.0;
+ }
+ };
+ auto get_quadrupole_Cm2 = [&]() -> double{
+ if (j.contains("(Q^*)^2") && j.contains("nQ")){
+ // Terms defined like in Gross&Vrabec; backwards-compatibility
+ double Qstar2 = j.at("(Q^*)^2");
+ return sqrt(Qstar2*(m*epsilon_over_k*pow(sigma_m, 5))/Qstar2factor);
+ }
+ else if (j.contains("Q_Cm2")){
+ return j.at("Q_Cm2");
+ }
+ else if (j.contains("Q_DA")){
+ return j.at("Q_Cm2").get<double>()*D_to_Cm*1e10;
+ }
+ else{
+ return 0.0;
+ }
+ };
+ ms(i) = m; sigma_ms(i) = sigma_m; epsks(i) = epsilon_over_k; mu_Cm(i) = get_dipole_Cm(); Q_Cm2(i) = get_quadrupole_Cm2();
+ nmu(i) = (j.contains("nmu") ? j["nmu"].get<double>() : 0.0);
+ nQ(i) = (j.contains("nQ") ? j["nQ"].get<double>() : 0.0);
+ i++;
+ };
+
+ using namespace SAFTpolar;
+ if (polar_model == "GrossVrabec"){
+ auto mustar2 = (mustar2factor*mu_Cm.pow(2)/(ms*epsks*sigma_ms.pow(3))).eval();
+ auto Qstar2 = (Qstar2factor*Q_Cm2.pow(2)/(ms*epsks*sigma_ms.pow(5))).eval();
+ auto polar = MultipolarContributionGrossVrabec(ms, sigma_ms*1e10, epsks, mustar2, nmu, Qstar2, nQ);
+ return SAFTVRMieMixture(SAFTVRMieMixture::build_chain(coeffs, kmat), std::move(polar));
+ }
+ else if (polar_model == "GubbinsTwu+Luckas"){
+ using MCGTL = MultipolarContributionGubbinsTwu<LuckasJIntegral, LuckasKIntegral>;
+ auto mubar2 = (mustar2factor*mu_Cm.pow(2)/(epsks*sigma_ms.pow(3))).eval();
+ auto Qbar2 = (Qstar2factor*Q_Cm2.pow(2)/(epsks*sigma_ms.pow(5))).eval();
+ auto polar = MCGTL(sigma_ms, epsks, mubar2, Qbar2);
+ return SAFTVRMieMixture(SAFTVRMieMixture::build_chain(coeffs, kmat), std::move(polar));
+ }
+ else if (polar_model == "GubbinsTwu+GubbinsTwu"){
+ using MCGG = MultipolarContributionGubbinsTwu<GubbinsTwuJIntegral, GubbinsTwuKIntegral>;
+ auto mubar2 = (mustar2factor*mu_Cm.pow(2)/(epsks*sigma_ms.pow(3))).eval();
+ auto Qbar2 = (Qstar2factor*Q_Cm2.pow(2)/(epsks*sigma_ms.pow(5))).eval();
+ auto polar = MCGG(sigma_ms, epsks, mubar2, Qbar2);
+ return SAFTVRMieMixture(SAFTVRMieMixture::build_chain(coeffs, kmat), std::move(polar));
+ }
+ else{
+ throw teqp::InvalidArgument("didn't understand this polar_model:"+polar_model);
+ }
+ }
+ }
+ else{
+ throw std::invalid_argument("you must provide names or coeffs, but not both");
+ }
+}
} /* namespace SAFTVRMie */
}; // namespace teqp
diff --git a/src/tests/catch_test_SAFTpolar.cxx b/src/tests/catch_test_SAFTpolar.cxx
index 6f84f68b3a8f3915629cd4272469e63c98898f6f..1a02bd81c3770330ccfa93097ddf4d1ed3499170 100644
--- a/src/tests/catch_test_SAFTpolar.cxx
+++ b/src/tests/catch_test_SAFTpolar.cxx
@@ -1,8 +1,10 @@
#include <catch2/catch_test_macros.hpp>
#include <catch2/catch_approx.hpp>
-
using Catch::Approx;
+#include "nlohmann/json.hpp"
+#include "teqp/cpp/teqpcpp.hpp"
+
#include "teqp/models/saft/correlation_integrals.hpp"
#include "teqp/models/saft/polar_terms.hpp"
#include "teqp/types.hpp"
@@ -25,6 +27,35 @@ TEST_CASE("Evaluation of K(xxx, yyy)", "[LuckasKnn]")
CHECK(Kval45 == Approx(0.01541).margin(0.02));
}
+TEST_CASE("Evaluation of J^{(n)}", "[checkJvals]")
+{
+ double Tstar = 2.0, rhostar = 0.4;
+ SECTION("Luckas v. G&T"){
+ for (auto n = 4; n < 37; ++n){
+ CAPTURE(n);
+ auto L = LuckasJIntegral(n).get_J(Tstar, rhostar);
+ auto G = GubbinsTwuJIntegral(n).get_J(Tstar, rhostar);
+ CHECK(L == Approx(G));
+ }
+ }
+}
+
+TEST_CASE("Evaluation of K^{(n,m)}", "[checkKvals]")
+{
+ double Tstar = 2.0, rhostar = 0.4;
+ SECTION("Luckas v. G&T"){
+ std::vector<std::tuple<int, int>> nm = {{222,333},{233,344},{334,445},{444,555}};
+ for (auto [n,m] : nm){
+ auto L = LuckasKIntegral(n,m).get_K(Tstar, rhostar);
+ auto G = GubbinsTwuKIntegral(n,m).get_K(Tstar, rhostar);
+ CAPTURE(n);
+ CAPTURE(m);
+ CHECK(L == Approx(G));
+ }
+ }
+}
+
+
using MCGTL = MultipolarContributionGubbinsTwu<LuckasJIntegral, LuckasKIntegral>;
using MCGG = MultipolarContributionGubbinsTwu<GubbinsTwuJIntegral, GubbinsTwuKIntegral>;
@@ -69,3 +100,78 @@ TEST_CASE("Check derivative of |x|", "[diffabs]")
CHECK(safe_abs2(std::complex<double>(-3.1, h)).imag()/h == -1);
}
}
+
+TEST_CASE("Check critical points with polarity terms", "[SAFTVRMiepolar]"){
+ double rhostar_guess_init = 0.27;
+ double Tstar_guess_init = 1.5;
+ double ek = 100;
+ double sigma_m = 3e-10;
+ auto j = nlohmann::json::parse(R"({
+ "kind": "SAFT-VR-Mie",
+ "model": {
+ "coeffs": [
+ {
+ "name": "Stockmayer126",
+ "m": 1.0,
+ "sigma_m": 3e-10,
+ "epsilon_over_k": 100,
+ "lambda_r": 12,
+ "lambda_a": 6,
+ "BibTeXKey": "me"
+ }
+ ]
+ }
+ })");
+ CHECK_NOTHROW(teqp::cppinterface::make_model(j));
+ auto nonpolar_model = teqp::cppinterface::make_model(j);
+
+// SECTION("Backwards-compat Gross&Vrabec"){
+// const bool print = false;
+// double Tstar_guess = Tstar_guess_init, rhostar_guess = rhostar_guess_init;
+// if (print) std::cout << "(mu^*)^2, T^*, rho^*" << std::endl;
+//
+// auto critpure = nonpolar_model->solve_pure_critical(Tstar_guess*ek, rhostar_guess/(N_A*pow(sigma_m, 3)));
+// Tstar_guess = std::get<0>(critpure)/ek;
+// rhostar_guess = std::get<1>(critpure)*N_A*pow(sigma_m, 3);
+// if (print) std::cout << "0, " << Tstar_guess << ", " << rhostar_guess << std::endl;
+//
+// for (double mustar2 = 0.1; mustar2 < 5; mustar2 += 0.1){
+// j["model"]["coeffs"][0]["(mu^*)^2"] = mustar2;
+// j["model"]["coeffs"][0]["nmu"] = 1;
+// auto model = teqp::cppinterface::make_model(j);
+// auto pure = model->solve_pure_critical(Tstar_guess*ek, rhostar_guess/(N_A*pow(sigma_m, 3)));
+// if (print) std::cout << mustar2 << ", " << std::get<0>(pure)/ek << ", " << std::get<1>(pure)*N_A*pow(sigma_m, 3) << std::endl;
+// Tstar_guess = std::get<0>(pure)/ek;
+// rhostar_guess = std::get<1>(pure)*N_A*pow(sigma_m, 3);
+// }
+// CHECK(Tstar_guess == Approx(2.29743));
+// CHECK(rhostar_guess == Approx(0.221054));
+// }
+
+ SECTION("With Gross&Vrabec"){
+ for (std::string polar_model : {"GrossVrabec", "GubbinsTwu+Luckas", "GubbinsTwu+GubbinsTwu"}){
+ const bool print = true;
+ double Tstar_guess = Tstar_guess_init, rhostar_guess = rhostar_guess_init;
+ if (print) std::cout << "(mu^*)^2, T^*, rho^*" << std::endl;
+
+ auto critpure = nonpolar_model->solve_pure_critical(Tstar_guess*ek, rhostar_guess/(N_A*pow(sigma_m, 3)));
+ Tstar_guess = std::get<0>(critpure)/ek;
+ rhostar_guess = std::get<1>(critpure)*N_A*pow(sigma_m, 3);
+ if (print) std::cout << "0, " << Tstar_guess << ", " << rhostar_guess << std::endl;
+
+ double mustar2factor = 1.0/(4*static_cast<double>(EIGEN_PI)*8.8541878128e-12*1.380649e-23);
+ j["model"]["polar_model"] = polar_model;
+ for (double mustar2 = 0.1; mustar2 < 5; mustar2 += 0.1){
+ j["model"]["coeffs"][0]["mu_Cm"] = sqrt(mustar2/mustar2factor*(ek*pow(sigma_m, 3)));
+ j["model"]["coeffs"][0]["nmu"] = 1;
+ auto model = teqp::cppinterface::make_model(j);
+ auto pure = model->solve_pure_critical(Tstar_guess*ek, rhostar_guess/(N_A*pow(sigma_m, 3)));
+ if (print) std::cout << mustar2 << ", " << std::get<0>(pure)/ek << ", " << std::get<1>(pure)*N_A*pow(sigma_m, 3) << std::endl;
+ Tstar_guess = std::get<0>(pure)/ek;
+ rhostar_guess = std::get<1>(pure)*N_A*pow(sigma_m, 3);
+ }
+// CHECK(Tstar_guess == Approx(2.29743));
+// CHECK(rhostar_guess == Approx(0.221054));
+ }
+ }
+}