From 16b74ffb62a8d2f9ebcba4c6197b825eeb4a02a3 Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Sat, 22 Jun 2024 20:58:13 -0400 Subject: [PATCH 01/13] First draft of Dufal association model --- CMakeLists.txt | 3 +- .../teqp/models/association/association.hpp | 153 ++++++++++++++---- .../models/association/association_types.hpp | 53 +++++- include/teqp/models/saft/genericsaft.hpp | 8 +- src/data/Dufal_assoc.cpp | 106 ++++++++++++ src/tests/catch_test_SAFTgeneric.cxx | 60 ++++++- 6 files changed, 342 insertions(+), 41 deletions(-) create mode 100644 src/data/Dufal_assoc.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index e2a6c5d2..2e32d968 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -174,8 +174,9 @@ if (NOT TEQP_NO_TEQPCPP) # Add the schema file target_sources(teqpcpp PRIVATE "${CMAKE_CURRENT_SOURCE_DIR}/interface/CPP/model_schemas.cpp") - # The FEANN data file too + # The data files too target_sources(teqpcpp PRIVATE "${CMAKE_CURRENT_SOURCE_DIR}/src/data/FEANN.cpp") + target_sources(teqpcpp PRIVATE "${CMAKE_CURRENT_SOURCE_DIR}/src/data/Dufal_assoc.cpp") target_sources(teqpcpp PRIVATE "${CMAKE_CURRENT_SOURCE_DIR}/src/data/PCSAFT.cpp") # target_compile_options(teqpcpp PRIVATE diff --git a/include/teqp/models/association/association.hpp b/include/teqp/models/association/association.hpp index 210152d1..63bfd71c 100644 --- a/include/teqp/models/association/association.hpp +++ b/include/teqp/models/association/association.hpp @@ -17,8 +17,9 @@ The implementation follows the approach of Langenbach for the index compression, #include "teqp/exceptions.hpp" #include - +#include "teqp/math/pow_templates.hpp" #include "teqp/models/association/association_types.hpp" +#include "teqp/json_tools.hpp" namespace teqp{ namespace association{ @@ -27,6 +28,7 @@ struct AssociationOptions{ std::map> interaction_partners; std::vector site_order; association::radial_dist radial_dist; + association::Delta_rules Delta_rule = association::Delta_rules::CR1; std::vector self_association_mask; double alpha = 0.5; double rtol = 1e-12, atol = 1e-12; @@ -39,6 +41,11 @@ inline void from_json(const nlohmann::json& j, AssociationOptions& o) { if (j.contains("max_iters")){ j.at("max_iters").get_to(o.max_iters); } } + +namespace DufalMatrices{ + extern const std::unordered_map bcoeffs; +} + /*** A general class for calculating association fractions and association energy for mixtures @@ -143,15 +150,11 @@ class Association{ } return D; } - auto toEig(const nlohmann::json& j, const std::string& k) const -> Eigen::ArrayXd{ + static auto toEig(const nlohmann::json& j, const std::string& k) -> Eigen::ArrayXd{ std::vector vec = j.at(k); return Eigen::Map(&vec[0], vec.size()); }; -// auto get_molecule_sites(const nlohmann::json&j){ -// std::vector> molecule_sites; -// return molecule_sites; -// } - auto get_association_options(const nlohmann::json&j){ + static auto get_association_options(const nlohmann::json&j){ AssociationOptions opt; if (j.contains("options")){ opt = j.at("options").get(); @@ -159,17 +162,50 @@ class Association{ return opt; } public: - const Eigen::ArrayXd b_m3mol, ///< The covolume b, in m^3/mol - beta, ///< The volume factor, dimensionless - epsilon_Jmol; ///< The association energy of each molecule, in J/mol const AssociationOptions options; const IndexMapper mapper; const Eigen::ArrayXXi D; - const radial_dist m_radial_dist; + const Delta_rules m_Delta_rule; + const std::variant datasidecar; - Association(const nlohmann::json&j) : Association(toEig(j, "b / m^3/mol"), toEig(j, "betaAB"), toEig(j, "epsAB/kB / J/mol"), j.at("molecule_sites"), get_association_options(j)){} - - Association(const Eigen::ArrayXd& b_m3mol, const Eigen::ArrayXd& beta, const Eigen::ArrayXd& epsilon_Jmol, const std::vector>& molecule_sites, const AssociationOptions& options) : b_m3mol(b_m3mol), beta(beta), epsilon_Jmol(epsilon_Jmol), options(options), mapper(make_mapper(molecule_sites, options)), D(make_D(mapper, options)), m_radial_dist(options.radial_dist){ + Association(const Eigen::ArrayXd& b_m3mol, const Eigen::ArrayXd& beta, const Eigen::ArrayXd& epsilon_Jmol, const std::vector>& molecule_sites, const AssociationOptions& options) : options(options), mapper(make_mapper(molecule_sites, options)), D(make_D(mapper, options)), m_Delta_rule(options.Delta_rule), datasidecar(CanonicalData{b_m3mol, beta, epsilon_Jmol, options.radial_dist}){ + } + Association(const decltype(datasidecar)& data, const std::vector>& molecule_sites, const AssociationOptions& options) : options(options), mapper(make_mapper(molecule_sites, options)), D(make_D(mapper, options)), m_Delta_rule(options.Delta_rule), datasidecar(data) { + } + static Association factory(const nlohmann::json& j){ + if (j.contains("Delta_rule")){ + std::string Delta_rule = j.at("Delta_rule"); + if (Delta_rule == "CR1"){ + CanonicalData data; + data.b_m3mol = toEig(j, "b / m^3/mol"); + data.beta = toEig(j, "beta"); + data.epsilon_Jmol = toEig(j, "epsilon / J/mol"); + auto options = get_association_options(j); + options.Delta_rule = Delta_rules::CR1; + return {data, j.at("molecule_sites"), options}; + } + else if(Delta_rule == "Dufal"){ + DufalData data; + + // Parameters for the dispersive part + data.sigma_m = toEig(j, "sigma / m"); + data.epsilon_Jmol = toEig(j, "epsilon / J/mol"); + data.lambda_r = toEig(j, "lambda_r"); + data.kmat = build_square_matrix(j.at("kmat")); + // Parameters for the associating part + data.epsilon_HB_Jmol = toEig(j, "epsilon_HB / J/mol"); + data.K_HB_m3 = toEig(j, "K_HB / m^3"); + + auto options = get_association_options(j); + options.Delta_rule = Delta_rules::Dufal; + data.apply_mixing_rules(); + return {data, j.at("molecule_sites"), options}; + } + else{ + + } + } + throw std::invalid_argument("The Delta_rule has not been specified"); } /** @@ -184,17 +220,62 @@ class Association{ using resulttype = std::common_type_t; // Type promotion, without the const-ness using Mat = Eigen::Array; auto Ngroups = mapper.to_CompSite.size(); - auto bmix = (molefracs*b_m3mol).sum(); - auto eta = bmix*rhomolar/4.0; - decltype(forceeval(eta)) g; - switch(m_radial_dist){ - case radial_dist::CS: - g = (2.0-eta)/(2.0*(1.0-eta)*(1.0-eta)*(1.0-eta)); break; - case radial_dist::KG: - g = 1.0 / (1.0 - 1.9*eta); break; - default: - throw std::invalid_argument("Bad radial distribution"); + // Calculate the radial_dist if it is meaningful + using eta_t = std::common_type_t; // Type promotion, without the const-ness + std::optional g; + if (m_Delta_rule == Delta_rules::CR1){ + const CanonicalData& d = std::get(datasidecar); + auto bmix = (molefracs*d.b_m3mol).sum(); + auto eta = bmix*rhomolar/4.0; + switch(d.radial_dist){ + case radial_dist::CS: + g = (2.0-eta)/(2.0*(1.0-eta)*(1.0-eta)*(1.0-eta)); break; + case radial_dist::KG: + g = 1.0 / (1.0 - 1.9*eta); break; + default: + throw std::invalid_argument("Bad radial distribution"); + } + } + + /// A helper function for the I integral representation of Dufal (http://dx.doi.org/10.1080/00268976.2015.1029027) + auto get_I_Dufal = [](const auto& Tstar, const auto& rhostar, const auto& lambda_r){ + + using result_t = std::decay_t>; + result_t summer = 0.0; + + std::decay_t rhostar_to_i = 1.0; + for (auto i = 0U; i <= 10; ++i){ + std::decay_t Tstar_to_j = 1.0; + for (auto j = 0U; i + j <= 10; ++j){ + double aij = 0.0, lambdar_to_k = 1.0; + for (auto k = 0; k <= 6; ++k){ + double bijk = DufalMatrices::bcoeffs.at(k)(i,j); + aij += bijk*lambdar_to_k; + lambdar_to_k *= lambda_r; + } + summer += aij*rhostar_to_i*Tstar_to_j; + Tstar_to_j *= Tstar; + } + rhostar_to_i *= rhostar; + } + return summer; + }; + + using rhostar_t = std::common_type_t; // Type promotion, without the const-ness + std::optional rhostar; + if (m_Delta_rule == Delta_rules::Dufal){ + const DufalData& d = std::get(datasidecar); + // Calculate the one-fluid vdW1 sigma from Eq. 40 of Dufal + std::decay_t sigma3_vdW1 = 0.0; + auto N = molefracs.size(); + for (auto i = 0U; i < N; ++i){ + for (auto j = 0U; j < N; ++j){ + double sigma3_ij_m3 = POW3((d.sigma_m[i] + d.sigma_m[j])/2.0); + sigma3_vdW1 += molefracs[i]*molefracs[j]*sigma3_ij_m3; + } + } + rhostar = rhomolar*N_A*sigma3_vdW1; } Mat Delta = Mat::Zero(Ngroups, Ngroups); @@ -204,12 +285,24 @@ class Association{ auto j = std::get<0>(mapper.to_CompSite.at(J)); using namespace teqp::constants; - // The CR1 rule is used to calculate the cross contributions - if (true){ // Is CR1 // TODO: also allow the other combining rule - auto b_ij = (b_m3mol[i] + b_m3mol[j])/2.0; - auto epsilon_ij_Jmol = (epsilon_Jmol[i] + epsilon_Jmol[j])/2.0; - auto beta_ij = sqrt(beta[i]*beta[j]); - Delta(I, J) = g*b_ij*beta_ij*(exp(epsilon_ij_Jmol/(R_CODATA2017*T))-1.0)/N_A; + + if (m_Delta_rule == Delta_rules::CR1){ + // The CR1 rule is used to calculate the cross contributions + const CanonicalData& d = std::get(datasidecar); + auto b_ij = (d.b_m3mol[i] + d.b_m3mol[j])/2.0; + auto epsilon_ij_Jmol = (d.epsilon_Jmol[i] + d.epsilon_Jmol[j])/2.0; + auto beta_ij = sqrt(d.beta[i]*d.beta[j]); + Delta(I, J) = g.value()*b_ij*beta_ij*(exp(epsilon_ij_Jmol/(R_CODATA2017*T))-1.0)/N_A; // m^3 + } + else if (m_Delta_rule == Delta_rules::Dufal){ + const DufalData& d = std::get(datasidecar); + auto Tstar = forceeval(T/d.EPSILONOVERKij_K(i,j)); + auto _I = get_I_Dufal(Tstar, rhostar.value(), d.LAMBDA_Rij(i, j)); + auto F_Meyer = exp(d.EPSILONOVERK_HBij_K(i, j)/T)-1.0; + Delta(I, J) = F_Meyer*d.KHBij_m3(i,j)*_I; + } + else{ + throw std::invalid_argument("Don't know what to do with this Delta rule"); } } } diff --git a/include/teqp/models/association/association_types.hpp b/include/teqp/models/association/association_types.hpp index ab49a302..ae08a9b4 100644 --- a/include/teqp/models/association/association_types.hpp +++ b/include/teqp/models/association/association_types.hpp @@ -1,4 +1,5 @@ #pragma once +#include "teqp/constants.hpp" namespace teqp { namespace association{ @@ -12,7 +13,7 @@ inline auto get_association_classes(const std::string& s) { else if (s == "3B") { return association_classes::a3B; } else if (s == "4C") { return association_classes::a4C; } else { - throw std::invalid_argument("bad association flag:" + s); + throw std::invalid_argument("bad association flag: " + s); } } @@ -22,9 +23,57 @@ inline auto get_radial_dist(const std::string& s) { if (s == "CS") { return radial_dist::CS; } else if (s == "KG") { return radial_dist::KG; } else { - throw std::invalid_argument("bad association flag:" + s); + throw std::invalid_argument("bad radial_dist flag: " + s); } } +enum class Delta_rules {not_set, CR1, Dufal}; + +inline auto get_Delta_rule(const std::string& s) { + if (s == "CR1") { return Delta_rules::CR1; } + else if (s == "Dufal") { return Delta_rules::Dufal; } + else { + throw std::invalid_argument("bad Delta_rule flag: " + s); + } +} + +struct CanonicalData{ + Eigen::ArrayXd b_m3mol, ///< The covolume b, in m^3/mol, one per component + beta, ///< The volume factor, dimensionless, one per component + epsilon_Jmol; ///< The association energy of each molecule, in J/mol, one per component + radial_dist radial_dist; +}; + +struct DufalData{ + // Parameters coming from the non-associating part, one per component + Eigen::ArrayXd sigma_m, epsilon_Jmol, lambda_r; + Eigen::ArrayXXd kmat; ///< Matrix of k_{ij} values + + // Parameters from the associating part, one per component + Eigen::ArrayXd epsilon_HB_Jmol, K_HB_m3; + + Eigen::ArrayXd SIGMA3ij_m3, EPSILONOVERKij_K, LAMBDA_Rij, EPSILONOVERK_HBij_K, KHBij_m3; + + void apply_mixing_rules(){ + std::size_t N = sigma_m.size(); + SIGMA3ij_m3.resize(N,N); + EPSILONOVERKij_K.resize(N,N); + LAMBDA_Rij.resize(N,N); + EPSILONOVERK_HBij_K.resize(N,N); + KHBij_m3.resize(N,N); + + for (auto i = 0U; i < N; ++i){ + for (auto j = 0U; j < N; ++j){ + SIGMA3ij_m3(i,j) = POW3((sigma_m[i] + sigma_m[j])/2.0); + EPSILONOVERKij_K(i, j) = (1-kmat(i,j))*sqrt(POW3(sigma_m[i]*sigma_m[j]))/SIGMA3ij_m3(i,j)*sqrt(epsilon_Jmol[i]*epsilon_Jmol[j])/constants::R_CODATA2017; + LAMBDA_Rij(i, j) = 3 + sqrt((lambda_r[i]-3)*(lambda_r[j]-3)); + EPSILONOVERK_HBij_K(i, j) = sqrt(epsilon_HB_Jmol[i]*epsilon_HB_Jmol[j])/constants::R_CODATA2017; + KHBij_m3(i,j) = POW3((3.0*K_HB_m3[i] + 3.0*K_HB_m3[j])/2.0); + } + } + } +}; + + } } diff --git a/include/teqp/models/saft/genericsaft.hpp b/include/teqp/models/saft/genericsaft.hpp index 9f6ccc25..11c0fdaf 100644 --- a/include/teqp/models/saft/genericsaft.hpp +++ b/include/teqp/models/saft/genericsaft.hpp @@ -33,13 +33,9 @@ struct GenericSAFT{ }; auto make_association(const nlohmann::json &j) -> AssociationTerms{ std::string kind = j.at("kind"); - if (kind == "canonical"){ - return association::Association(j.at("model")); + if (kind == "canonical" || kind == "Dufal"){ + return association::Association::factory(j.at("model")); } - // TODO: Add the new Dufal association term -// else if (kind == "Dufal-Mie"){ -// return std::make_unique(j.at("model")); -// } else{ throw std::invalid_argument("Not valid association kind:" + kind); } diff --git a/src/data/Dufal_assoc.cpp b/src/data/Dufal_assoc.cpp new file mode 100644 index 00000000..714d557a --- /dev/null +++ b/src/data/Dufal_assoc.cpp @@ -0,0 +1,106 @@ +#include + +namespace teqp::association::DufalMatrices{ + +auto b_0 = (Eigen::Matrix() << +0.0132970702182068, -0.0177199122935443, 0.0293736747694974, -0.0205527304404423, 0.00861683420907605, -0.00228505275303600, 0.000390171133200072, -4.26035888869942e-05, 2.86246920519487e-06, -1.07315320963937e-07, 1.70912976772329e-09, +-0.0465504528847432, 0.332597325549352, -0.326575316241193, 0.144653671541451, -0.0363193315289496, 0.00569934220115537, -0.000581966173216051, 3.83608167089024e-05, -1.50305409953983e-06, 2.66749257811143e-08, 0, +0.164972499633366, -0.974898725377830, 0.919082550772666, -0.367978443660284, 0.0788054156983951, -0.00981102799831725, 0.000717901835772044, -2.91191052989125e-05, 5.17207032026779e-07, 0, 0, +-0.553080054304108, 1.07124021914524, -0.914915281734471, 0.328132391309741, -0.0577909160509185, 0.00539941935098940, -0.000249522541631115, 4.36324500072586e-06, 0, 0, 0, +0.697481173735912, -0.127802319197572, 0.238559496985344, -0.124364954974360, 0.0183583032370354, -0.00133946361388127, 3.68888649118614e-05, 0, 0, 0, 0, +-0.00682258598593205, -0.296768597044265, 0.346077751701231, 0.0122496582163678, -0.000805951611068984, 5.01775524378700e-05, 0, 0, 0, 0, 0, +-2.46482416179796, -1.45370322973416, -0.448827080921154, -0.00566229179136722, -0.000180870029200998, 0, 0, 0, 0, 0, 0, +8.78388694047369, 3.23807384513205, 0.281816142695178, 0.00340169105539079, 0, 0, 0, 0, 0, 0, 0, +-13.5178089781880, -2.48217551606281, -0.0783334040233511, 0, 0, 0, 0, 0, 0, 0, 0, +9.42415649943917, 0.687363680163044, 0, 0, 0, 0, 0, 0, 0, 0, 0, +-2.46151453173016, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0).finished(); + +auto b_1 = (Eigen::Matrix() << +0.000556479463564548, -0.000282932524693843, -0.00212156862728691, 0.00224197388058698, -0.00115385075260200, 0.000345261305021541, -6.36772702454557e-05, 7.32538316171651e-06, -5.10881831259873e-07, 1.97049247692837e-08, -3.21285622252695e-10, +0.0492332122696642, -0.134456183191719, 0.108517185632299, -0.0425401431513956, 0.00954409805564329, -0.00131479859839035, 0.000113485041044446, -5.97590027689909e-06, 1.73112888670222e-07, -2.01831274082934e-09, 0, +-0.158447251265296, 0.430576656790814, -0.342386658999542, 0.126318819450003, -0.0257089028169037, 0.00310740070593876, -0.000226247702797060, 9.40337471999922e-06, -1.75927011626452e-07, 0, 0, +0.130961541597078, -0.358137806097004, 0.323521260711548, -0.113820803061658, 0.0192702466504444, -0.00169264438876774, 7.20470394807155e-05, -1.12895786020720e-06, 0, 0, 0, +-0.0678784049001058, -0.159716521155965, -0.0408260244006866, 0.0397261699377944, -0.00668853458452261, 0.000490515978818942, -1.26640904345952e-05, 0, 0, 0, 0, +0.611672146147809, 0.257963755040360, -0.150206376568020, 0.00387683716563604, 0.000113887448052591, -2.67212361539355e-05, 0, 0, 0, 0, 0, +-1.27512696874714, 0.421048118713917, 0.122792263121120, -0.00158386412295998, 0.000198370689434891, 0, 0, 0, 0, 0, 0, +0.269957785723510, -0.899401835359301, -0.0593934567946635, -0.00110977335239395, 0, 0, 0, 0, 0, 0, 0, +1.49379616940916, 0.628459498670159, 0.0177032259427776, 0, 0, 0, 0, 0, 0, 0, 0, +-1.50913781396606, -0.166624388301135, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0.444942467424175, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0).finished(); + +auto b_2 = (Eigen::Matrix() << +4.68753836985661e-05, -0.000201534029276569, 0.000408971640196116, -0.000334428221392103, 0.000154951071291061, -4.39272547633533e-05, 7.86877804626315e-06, -8.92538365355001e-07, 6.20541126317064e-08, -2.40904019797940e-09, 3.99225753392296e-11, +-0.00636142804034125, 0.0148870088793584, -0.0112930088940554, 0.00423537018032177, -0.000909633136025822, 0.000118406467168496, -9.31205117458790e-06, 4.07198179416416e-07, -7.24002263850399e-09, -2.85054817786792e-11, 0, +0.0211014984424621, -0.0495382445219106, 0.0371909493534914, -0.0132662359974178, 0.00264252348906371, -0.000316720398760111, 2.32682269109226e-05, -9.93969228161767e-07, 1.93247731973246e-08, 0, 0, +-0.0204522546881708, 0.0456627841733661, -0.0375241588712114, 0.0125888572116524, -0.00204736822270412, 0.000171672501041056, -6.91918528995779e-06, 1.02234151256019e-07, 0, 0, 0, +0.0267157884350682, 0.00604022182900132, 0.00822450412680841, -0.00502124889994980, 0.000803630265796884, -5.55280746686634e-05, 1.30402695735287e-06, 0, 0, 0, 0, +-0.131004401410042, -0.00985163509226622, 0.0128544254818812, -0.000216912930642026, -4.55721397631024e-05, 4.78472396442476e-06, 0, 0, 0, 0, 0, +0.260822991454304, -0.0626432527852684, -0.00959839254636274, 0.000223341727285368, -2.78886526475732e-05, 0, 0, 0, 0, 0, 0, +-0.179706752593973, 0.0989212484451187, 0.00388197187439871, 0.000151511550199890, 0, 0, 0, 0, 0, 0, 0, +-0.0315900588187112, -0.0607681438566115, -0.00141908354521815, 0, 0, 0, 0, 0, 0, 0, 0, +0.0926908832419059, 0.0151506077043458, 0, 0, 0, 0, 0, 0, 0, 0, 0, + -0.0312717011022248, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0).finished(); + +auto b_3 = (Eigen::Matrix() << +-1.52750755540612e-06, 6.65734616244750e-06, -1.49163457856162e-05, 1.28251715836463e-05, -6.13363502620892e-06, 1.77374828289499e-06, -3.21719278338969e-07, 3.67816571475073e-08, -2.57070860278200e-09, 1.00187848111417e-10, -1.66615681224878e-12, +0.000352707112753263, -0.000761494022700993, 0.000551126114938543, -0.000197194099229223, 3.97623327839760e-05, -4.65836276873807e-06, 2.92850782905453e-07, -5.67255537711216e-09, -3.50155347187566e-10, 1.64319946681537e-11, 0, +-0.00127835416513710, 0.00272631887421459, -0.00195876152374924, 0.000680170219619246, -0.000133171449466011, 1.58508191731109e-05, -1.17119945008334e-06, 5.09014996522337e-08, -1.01109171728229e-09, 0, 0, +0.00164258140843764, -0.00294360803306568, 0.00215659810736663, -0.000685370482291315, 0.000107391347801278, -8.69254977310800e-06, 3.38632787324714e-07, -4.90436065438310e-09, 0, 0, 0, +-0.00296574649489361, 0.000737340818584580, -0.000762065788926240, 0.000315097078455367, -4.64683975120556e-05, 3.00942802764189e-06, -6.47808066443911e-08, 0, 0, 0, 0, +0.0104319873447675, -0.000775429459116355, -0.000379064014723981, -1.60783833898331e-05, 4.71880061366660e-06, -3.40226871507608e-07, 0, 0, 0, 0, 0, +-0.0195358884282327, 0.00420143444237083, 0.000297766652427410, -8.36498131244475e-06, 1.51044087021239e-06, 0, 0, 0, 0, 0, 0, +0.0164362361737192, -0.00527606737824625, -9.27089410962249e-05, -9.32930119768768e-06, 0, 0, 0, 0, 0, 0, 0, +-0.00439666318131193, 0.00287829011566385, 5.24103094975034e-05, 0, 0, 0, 0, 0, 0, 0, 0, +-0.00159566365617165, -0.000667764807450442, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0.000861407443430205, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0).finished(); + +auto b_4 = (Eigen::Matrix() << +-3.67201230932920e-08, 6.65433123492297e-08, 8.09753253026481e-08, -1.29349636796981e-07, 7.75034942364271e-08, -2.50442449553800e-08, 4.81969764886444e-09, -5.68335864091982e-10, 4.02544652605731e-11, -1.57057135789188e-12, 2.59018065150187e-14, +-1.00533098186312e-05, 2.04973785455268e-05, -1.42364525555262e-05, 4.84859626533717e-06, -9.03287394521535e-07, 8.95403985256519e-08, -3.05175477078061e-09, -2.51206227550094e-10, 2.85888337544568e-11, -8.35713691699930e-13, 0, +3.98289826660923e-05, -7.92066147384541e-05, 5.49066841157859e-05, -1.86402735397850e-05, 3.59640611239294e-06, -4.25136583233798e-07, 3.14701885589044e-08, -1.37942010763254e-09, 2.76625894833152e-11, 0, 0, +-6.22606152603800e-05, 9.78661859355705e-05, -6.54711127544114e-05, 1.98394792303955e-05, -3.01458972861098e-06, 2.38076581975360e-07, -9.12137235540096e-09, 1.33298781401136e-10, 0, 0, 0, +0.000124810638661335, -5.13612052054272e-05, 3.09432853737792e-05, -1.02009751236812e-05, 1.39714413579629e-06, -8.55034278983589e-08, 1.71150044927963e-09, 0, 0, 0, 0, +-0.000374044099050907, 5.60796050391817e-05, 1.84626127086989e-06, 1.34288528575647e-06, -1.94614655230730e-07, 1.14157863743549e-08, 0, 0, 0, 0, 0, +0.000672078359176232, -0.000138031320636187, -3.14894396008512e-06, 5.76785692862702e-08, -3.88710150387342e-08, 0, 0, 0, 0, 0, 0, +-0.000609145837554261, 0.000148275594790438, 2.61353677755675e-07, 2.91087061127825e-07, 0, 0, 0, 0, 0, 0, 0, +0.000252724106819134, -7.36076396864788e-05, -1.00642768627400e-06, 0, 0, 0, 0, 0, 0, 0, 0, +-2.49255867787548e-05, 1.59720037817341e-05, 0, 0, 0, 0, 0, 0, 0, 0, 0, + -7.70341617155424e-06, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0).finished(); + +auto b_5 = (Eigen::Matrix() << +1.88048156944327e-09, -5.20041750295709e-09, 4.35881692094647e-09, -1.90967917025464e-09, 4.36098784939288e-10, -4.66145781784596e-11, 2.56821311477203e-13, 4.53836639240824e-13, -4.35569028016719e-14, 1.51660179879606e-15, -1.34231785680549e-17, +1.44492319539037e-07, -2.81833237588783e-07, 1.88795822312826e-07, -6.12284776821388e-08, 1.04054640162243e-08, -7.91894115995732e-10, -1.78097006314277e-11, 8.46173387258907e-12, -6.29468789021339e-13, 1.62914975397206e-14, 0, +-6.16808609401052e-07, 1.16624214027143e-06, -7.85736422030061e-07, 2.61841021616088e-07, -4.99004695268396e-08, 5.85970807128971e-09, -4.33266245043971e-10, 1.90360900947593e-11, -3.82506468133384e-13, 0, 0, +1.08977129458469e-06, -1.59381377689636e-06, 9.99496177147637e-07, -2.91595588554532e-07, 4.32652379690724e-08, -3.36277868350366e-09, 1.28216601435668e-10, -1.91702973428801e-12, 0, 0, 0, +-2.25971486750789e-06, 1.11917575440765e-06, -5.65159033070467e-07, 1.62627589914010e-07, -2.09918759575187e-08, 1.22756291392471e-09, -2.32414275161676e-11, 0, 0, 0, 0, +6.17883043382806e-06, -1.20280382926209e-06, 8.66499935530469e-08, -3.13263066251902e-08, 3.53265909236495e-09, -1.81505995831496e-10, 0, 0, 0, 0, 0, +-1.08048264112286e-05, 2.17366900821365e-06, -1.53022515538303e-08, 2.11500280833261e-09, 4.74260725560735e-10, 0, 0, 0, 0, 0, 0, +1.01466386059898e-05, -2.10287523493265e-06, 1.91184417519603e-08, -4.48828024315313e-09, 0, 0, 0, 0, 0, 0, 0, +-4.89026741776877e-06, 9.75457041464072e-07, 1.00074660431270e-08, 0, 0, 0, 0, 0, 0, 0, 0, +1.04333009405285e-06, -2.00783047336916e-07, 0, 0, 0, 0, 0, 0, 0, 0, 0, + -5.01517451094617e-08, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0).finished(); + +auto b_6 = (Eigen::Matrix() << +-1.84421844661105e-11, 5.39073389353030e-11, -5.61963272868663e-11, 3.23187813505929e-11, -1.13265016166587e-11, 2.59326201844053e-12, -4.01407631594698e-13, 4.20546124661511e-14, -2.87018830339552e-15, 1.15159946184350e-16, -2.05485407533207e-18, +-8.27665331374217e-10, 1.55896078845997e-09, -1.01231143749898e-09, 3.13409590234525e-10, -4.81050397623846e-11, 2.30261130144263e-12, 3.82442150265529e-13, -7.05793035532524e-14, 4.59807633591802e-15, -1.12579371971295e-16, 0, +3.74547734805679e-09, -6.83326218348459e-09, 4.50304220343962e-09, -1.47842817249421e-09, 2.78924412091172e-10, -3.25532731205188e-11, 2.40008362810406e-12, -1.05309962007498e-13, 2.11144278824276e-15, 0, 0, +-7.12756897329521e-09, 1.00393010742688e-08, -6.02106444139650e-09, 1.70643459526364e-09, -2.48716300671280e-10, 1.91426474551801e-11, -7.31095625646731e-13, 1.12112566023581e-14, 0, 0, 0, +1.48400265810797e-08, -8.18214743502133e-09, 3.81858691213582e-09, -1.00814014221601e-09, 1.24324450815098e-10, -7.02051970613448e-12, 1.27771307749498e-13, 0, 0, 0, 0, +-3.82128399646057e-08, 8.61165773602971e-09, -1.03558713701492e-09, 2.39230607884171e-10, -2.35824243963634e-11, 1.10391700857033e-12, 0, 0, 0, 0, 0, +6.55633333711198e-08, -1.31376750286488e-08, 3.55287199097211e-10, -2.85768231341124e-11, -2.19531090111075e-12, 0, 0, 0, 0, 0, 0, +-6.26595167233416e-08, 1.18345081383843e-08, -1.95221464285268e-10, 2.70989222978676e-11, 0, 0, 0, 0, 0, 0, 0, +3.22964489272539e-08, -5.24025196193952e-09, -4.14196554462777e-11, 0, 0, 0, 0, 0, 0, 0, 0, +-8.35760085684382e-09, 1.03868339224250e-09, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 8.47595203890549e-10, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0).finished(); + +extern const std::unordered_map bcoeffs = { + {0, b_0}, + {1, b_1}, + {2, b_2}, + {3, b_3}, + {4, b_4}, + {5, b_5}, + {6, b_6} +}; + +} diff --git a/src/tests/catch_test_SAFTgeneric.cxx b/src/tests/catch_test_SAFTgeneric.cxx index bf13653f..34a31570 100644 --- a/src/tests/catch_test_SAFTgeneric.cxx +++ b/src/tests/catch_test_SAFTgeneric.cxx @@ -29,8 +29,9 @@ TEST_CASE("Benchmark generic PC-SAFT+Association model", "[SAFTgeneric]"){ "kind": "canonical", "model": { "b / m^3/mol": [0.0000145], - "betaAB": [0.0692], - "epsAB/kB / J/mol": [2500.7], + "beta": [0.0692], + "Delta_rule": "CR1", + "epsilon / J/mol": [16655.0], "molecule_sites": [["e","e","H","H"]], "options": {"radial_dist": "CS"} } @@ -94,3 +95,58 @@ TEST_CASE("Benchmark generic PC-SAFT+Association model", "[SAFTgeneric]"){ }; } +TEST_CASE("Benchmark Dufal water model", "[SAFTgeneric]"){ + auto contents = R"( + { + "nonpolar": { + "kind": "SAFT-VR-Mie", + "model": { + "coeffs": [ + { + "name": "Water", + "BibTeXKey": "Dufal-2015", + "m": 1.0, + "sigma_Angstrom": 3.0555, + "epsilon_over_k": 418.00, + "lambda_r": 35.823, + "lambda_a": 6.0 + } + ] + } + }, + "association": { + "kind": "Dufal", + "model": { + "sigma / m": [3.0555e-10], + "epsilon / J/mol": [3475.445374388054], + "lambda_r": [35.823], + "epsilon_HB / J/mol": [13303.140189045183], + "K_HB / m^3": [496.66e-30], + "kmat": [[1.0]], + "Delta_rule": "Dufal", + "molecule_sites": [["e","e","H","H"]] + } + } + } + )"_json; + + nlohmann::json maincontents = { + {"kind", "genericSAFT"}, + {"model", contents} + }; + auto model = teqp::cppinterface::make_model(maincontents); + auto z = (Eigen::ArrayXd(1) << 1.0).finished(); + + BENCHMARK("Ar00"){ + return model->get_Ar00(300, 10000, z); + }; + BENCHMARK("Ar11"){ + return model->get_Ar11(300, 10000, z); + }; + BENCHMARK("Ar02"){ + return model->get_Ar02(300, 10000, z); + }; + BENCHMARK("Ar20"){ + return model->get_Ar20(300, 10000, z); + }; +} From 42c94d5ebf8ae91e73660ff545a531c7e351aafc Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Wed, 3 Jul 2024 15:02:25 -0400 Subject: [PATCH 02/13] Avoid the Ar00n function in both calling orders --- include/teqp/cpp/teqpcpp.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/include/teqp/cpp/teqpcpp.hpp b/include/teqp/cpp/teqpcpp.hpp index f64000ed..39fc2d4d 100644 --- a/include/teqp/cpp/teqpcpp.hpp +++ b/include/teqp/cpp/teqpcpp.hpp @@ -37,7 +37,6 @@ using REMatrixd = Eigen::Ref Date: Wed, 3 Jul 2024 15:03:04 -0400 Subject: [PATCH 03/13] Rename radial_dist to radial_dists --- include/teqp/models/CPA.hpp | 2 +- include/teqp/models/association/association.hpp | 6 +++--- include/teqp/models/association/association_types.hpp | 8 ++++---- 3 files changed, 8 insertions(+), 8 deletions(-) diff --git a/include/teqp/models/CPA.hpp b/include/teqp/models/CPA.hpp index 0500efb2..6f318988 100644 --- a/include/teqp/models/CPA.hpp +++ b/include/teqp/models/CPA.hpp @@ -17,7 +17,7 @@ namespace CPA { template auto POW2(X x) { return x * x; }; template auto POW3(X x) { return x * POW2(x); }; -using association::radial_dist; +using radial_dist = association::radial_dists; using association::association_classes; using association::get_radial_dist; using association::get_association_classes; diff --git a/include/teqp/models/association/association.hpp b/include/teqp/models/association/association.hpp index 63bfd71c..77e4b8b2 100644 --- a/include/teqp/models/association/association.hpp +++ b/include/teqp/models/association/association.hpp @@ -27,7 +27,7 @@ namespace association{ struct AssociationOptions{ std::map> interaction_partners; std::vector site_order; - association::radial_dist radial_dist; + association::radial_dists radial_dist; association::Delta_rules Delta_rule = association::Delta_rules::CR1; std::vector self_association_mask; double alpha = 0.5; @@ -229,9 +229,9 @@ class Association{ auto bmix = (molefracs*d.b_m3mol).sum(); auto eta = bmix*rhomolar/4.0; switch(d.radial_dist){ - case radial_dist::CS: + case radial_dists::CS: g = (2.0-eta)/(2.0*(1.0-eta)*(1.0-eta)*(1.0-eta)); break; - case radial_dist::KG: + case radial_dists::KG: g = 1.0 / (1.0 - 1.9*eta); break; default: throw std::invalid_argument("Bad radial distribution"); diff --git a/include/teqp/models/association/association_types.hpp b/include/teqp/models/association/association_types.hpp index ae08a9b4..ec9f923e 100644 --- a/include/teqp/models/association/association_types.hpp +++ b/include/teqp/models/association/association_types.hpp @@ -17,11 +17,11 @@ inline auto get_association_classes(const std::string& s) { } } -enum class radial_dist { CS, KG }; +enum class radial_dists { CS, KG }; inline auto get_radial_dist(const std::string& s) { - if (s == "CS") { return radial_dist::CS; } - else if (s == "KG") { return radial_dist::KG; } + if (s == "CS") { return radial_dists::CS; } + else if (s == "KG") { return radial_dists::KG; } else { throw std::invalid_argument("bad radial_dist flag: " + s); } @@ -41,7 +41,7 @@ struct CanonicalData{ Eigen::ArrayXd b_m3mol, ///< The covolume b, in m^3/mol, one per component beta, ///< The volume factor, dimensionless, one per component epsilon_Jmol; ///< The association energy of each molecule, in J/mol, one per component - radial_dist radial_dist; + radial_dists radial_dist; }; struct DufalData{ From 2b8d507891d07e3198eea3ecadaa00c1bf7a291e Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Wed, 3 Jul 2024 15:20:05 -0400 Subject: [PATCH 04/13] One more radial_dists --- src/tests/catch_test_association.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tests/catch_test_association.cxx b/src/tests/catch_test_association.cxx index 1ce6af37..15853aff 100644 --- a/src/tests/catch_test_association.cxx +++ b/src/tests/catch_test_association.cxx @@ -29,7 +29,7 @@ TEST_CASE("Test ethanol + water association", "[association]"){ std::vector> molecule_sites = {{"e", "H"}, {"e", "e", "H", "H"}}; association::AssociationOptions opt; - opt.radial_dist = association::radial_dist::CS; + opt.radial_dist = association::radial_dists::CS; opt.max_iters = 1000; opt.interaction_partners = {{"e", {"H",}}, {"H", {"e",}}}; From 6570b3346dbed589f8e7e9e945a3e385d3be15e0 Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Wed, 3 Jul 2024 16:22:26 -0400 Subject: [PATCH 05/13] Fix erratum in evaluation of KHB --- include/teqp/models/association/association_types.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/teqp/models/association/association_types.hpp b/include/teqp/models/association/association_types.hpp index ec9f923e..78c124a2 100644 --- a/include/teqp/models/association/association_types.hpp +++ b/include/teqp/models/association/association_types.hpp @@ -68,7 +68,7 @@ struct DufalData{ EPSILONOVERKij_K(i, j) = (1-kmat(i,j))*sqrt(POW3(sigma_m[i]*sigma_m[j]))/SIGMA3ij_m3(i,j)*sqrt(epsilon_Jmol[i]*epsilon_Jmol[j])/constants::R_CODATA2017; LAMBDA_Rij(i, j) = 3 + sqrt((lambda_r[i]-3)*(lambda_r[j]-3)); EPSILONOVERK_HBij_K(i, j) = sqrt(epsilon_HB_Jmol[i]*epsilon_HB_Jmol[j])/constants::R_CODATA2017; - KHBij_m3(i,j) = POW3((3.0*K_HB_m3[i] + 3.0*K_HB_m3[j])/2.0); + KHBij_m3(i,j) = POW3((cbrt(K_HB_m3[i]) + cbrt(K_HB_m3[j]))/2.0); // Published erratum in Dufal: https://doi.org/10.1080/00268976.2017.1402604 } } } From 275a1e2351f1c11d8453a3cc1df4edf28fd1a381 Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Wed, 3 Jul 2024 16:24:06 -0400 Subject: [PATCH 06/13] Benchmark calculations of Delta with Dufal and normal association --- src/tests/catch_test_association.cxx | 61 ++++++++++++++++++++++++++++ 1 file changed, 61 insertions(+) diff --git a/src/tests/catch_test_association.cxx b/src/tests/catch_test_association.cxx index 15853aff..38216621 100644 --- a/src/tests/catch_test_association.cxx +++ b/src/tests/catch_test_association.cxx @@ -7,6 +7,8 @@ using Catch::Matchers::WithinRel; #include "teqp/cpp/teqpcpp.hpp" #include "teqp/models/association/association.hpp" +#include "teqp/models/CPA.hpp" +#include "teqp/cpp/deriv_adapter.hpp" using namespace teqp; @@ -190,3 +192,62 @@ TEST_CASE("Trace ethanol + water isobaric VLE with CPA", "[associationVLE]"){ auto o = model->trace_VLE_isobar_binary(p, T, rhovecL, rhovecV, opt); CHECK(o.size() > 10); } + + +TEST_CASE("Benchmark assocation evaluations", "[associationbench]"){ + nlohmann::json water = { + {"a0i / Pa m^6/mol^2", 0.12277}, {"bi / m^3/mol", 0.0000145}, {"c1", 0.6736}, {"Tc / K", 647.13}, + {"epsABi / J/mol",}, {"betaABi", 0.0692}, {"sites", {"e","e","H","H"}} + }; + nlohmann::json jCPA = { + {"cubic", "SRK"}, + {"radial_dist", "CS"}, + {"pures", {water}}, + {"R_gas / J/mol/K", 8.31446261815324} + }; + nlohmann::json j = { + {"kind", "CPA"}, + {"validate", false}, + {"model", jCPA} + }; + + std::vector> molecule_sites = {{"e","e","H","H"}}; + + auto get_canon = [&](){ + auto b_m3mol = (Eigen::ArrayXd(1) << 0.0000145).finished(); + auto beta = (Eigen::ArrayXd(1) << 0.0692).finished(); + auto epsilon_Jmol = (Eigen::ArrayXd(1) << 16655.0).finished(); + association::AssociationOptions options; + options.Delta_rule = association::Delta_rules::CR1; + return association::Association(b_m3mol, beta, epsilon_Jmol, molecule_sites, options); + }; + auto get_Dufal = [&](){ + association::AssociationOptions options; + options.Delta_rule = association::Delta_rules::Dufal; + + teqp::association::DufalData data; + auto oneeig = [](double x){ return (Eigen::ArrayXd(1) << x).finished(); }; + double R = constants::R_CODATA2017; + data.sigma_m = oneeig(3.055e-10); + data.epsilon_Jmol = oneeig(418.0*R); + data.lambda_r = oneeig(35.823); + data.kmat = build_square_matrix(R"([[0.0]])"_json); + // Parameters for the associating part + data.epsilon_HB_Jmol = oneeig(1600.0*R); + data.K_HB_m3 = oneeig(496.66e-30); + data.apply_mixing_rules(); + return association::Association(data, molecule_sites, options); + }; + auto canon = get_canon(); + auto Dufal = get_Dufal(); + + auto z = (Eigen::ArrayXd(1) << 1.0).finished(); + BENCHMARK("time Delta with canonical"){ + return canon.get_Delta(300.0, 1/3e-5, z); + }; + BENCHMARK("time Delta with Dufal"){ + return Dufal.get_Delta(300.0, 1/3e-5, z); + }; + std::cout << canon.get_Delta(300.0, 1/3e-5, z) << std::endl; + std::cout << Dufal.get_Delta(300.0, 1/3e-5, z) << std::endl; +} From d5eb4a284f690027cc1dadb0612d9a5132d3b8f9 Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Wed, 3 Jul 2024 16:24:27 -0400 Subject: [PATCH 07/13] One more fix of radial_dist --- include/teqp/models/CPA.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/teqp/models/CPA.hpp b/include/teqp/models/CPA.hpp index 6f318988..7fe0cfc2 100644 --- a/include/teqp/models/CPA.hpp +++ b/include/teqp/models/CPA.hpp @@ -338,7 +338,7 @@ inline auto CPAfactory(const nlohmann::json &j){ // This is the backwards compatible approach // with the old style of defining the association class {1,2B...} std::vector classes; - radial_dist dist = get_radial_dist(j.at("radial_dist")); + auto dist = get_radial_dist(j.at("radial_dist")); std::valarray epsABi(N), betaABi(N), bi(N); std::size_t i = 0; for (auto p : j.at("pures")) { From 516b5200220f2ba3e070ebd79634dab9f607075f Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Wed, 3 Jul 2024 16:36:25 -0400 Subject: [PATCH 08/13] Match Dufal (and Clapeyron) exactly --- src/tests/catch_test_association.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tests/catch_test_association.cxx b/src/tests/catch_test_association.cxx index 38216621..c1477a62 100644 --- a/src/tests/catch_test_association.cxx +++ b/src/tests/catch_test_association.cxx @@ -228,7 +228,7 @@ TEST_CASE("Benchmark assocation evaluations", "[associationbench]"){ teqp::association::DufalData data; auto oneeig = [](double x){ return (Eigen::ArrayXd(1) << x).finished(); }; double R = constants::R_CODATA2017; - data.sigma_m = oneeig(3.055e-10); + data.sigma_m = oneeig(3.0555e-10); data.epsilon_Jmol = oneeig(418.0*R); data.lambda_r = oneeig(35.823); data.kmat = build_square_matrix(R"([[0.0]])"_json); From a242926b4e4f54843338d4f780fc5e97fa34bcfc Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Wed, 3 Jul 2024 16:36:41 -0400 Subject: [PATCH 09/13] Chunk is no longer needed --- src/tests/catch_test_association.cxx | 17 +---------------- 1 file changed, 1 insertion(+), 16 deletions(-) diff --git a/src/tests/catch_test_association.cxx b/src/tests/catch_test_association.cxx index c1477a62..666ff5e7 100644 --- a/src/tests/catch_test_association.cxx +++ b/src/tests/catch_test_association.cxx @@ -194,22 +194,7 @@ TEST_CASE("Trace ethanol + water isobaric VLE with CPA", "[associationVLE]"){ } -TEST_CASE("Benchmark assocation evaluations", "[associationbench]"){ - nlohmann::json water = { - {"a0i / Pa m^6/mol^2", 0.12277}, {"bi / m^3/mol", 0.0000145}, {"c1", 0.6736}, {"Tc / K", 647.13}, - {"epsABi / J/mol",}, {"betaABi", 0.0692}, {"sites", {"e","e","H","H"}} - }; - nlohmann::json jCPA = { - {"cubic", "SRK"}, - {"radial_dist", "CS"}, - {"pures", {water}}, - {"R_gas / J/mol/K", 8.31446261815324} - }; - nlohmann::json j = { - {"kind", "CPA"}, - {"validate", false}, - {"model", jCPA} - }; +TEST_CASE("Benchmark association evaluations", "[associationbench]"){ std::vector> molecule_sites = {{"e","e","H","H"}}; From ea25baf05bc274b88c904df1e693a2ffc4d82937 Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Thu, 4 Jul 2024 07:36:44 -0400 Subject: [PATCH 10/13] Two fixes to use of radial_dist --- include/teqp/models/association/association.hpp | 1 + src/tests/catch_test_association.cxx | 1 + 2 files changed, 2 insertions(+) diff --git a/include/teqp/models/association/association.hpp b/include/teqp/models/association/association.hpp index 77e4b8b2..c0b11c77 100644 --- a/include/teqp/models/association/association.hpp +++ b/include/teqp/models/association/association.hpp @@ -182,6 +182,7 @@ class Association{ data.epsilon_Jmol = toEig(j, "epsilon / J/mol"); auto options = get_association_options(j); options.Delta_rule = Delta_rules::CR1; + data.radial_dist = options.radial_dist; return {data, j.at("molecule_sites"), options}; } else if(Delta_rule == "Dufal"){ diff --git a/src/tests/catch_test_association.cxx b/src/tests/catch_test_association.cxx index 666ff5e7..3fe4e9b8 100644 --- a/src/tests/catch_test_association.cxx +++ b/src/tests/catch_test_association.cxx @@ -204,6 +204,7 @@ TEST_CASE("Benchmark association evaluations", "[associationbench]"){ auto epsilon_Jmol = (Eigen::ArrayXd(1) << 16655.0).finished(); association::AssociationOptions options; options.Delta_rule = association::Delta_rules::CR1; + options.radial_dist = association::radial_dists::CS; return association::Association(b_m3mol, beta, epsilon_Jmol, molecule_sites, options); }; auto get_Dufal = [&](){ From 30a38bcfe3e94dc6631db692041aa77fbd877607 Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Thu, 4 Jul 2024 07:37:02 -0400 Subject: [PATCH 11/13] Check for invalid Delta rule --- include/teqp/models/association/association.hpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/include/teqp/models/association/association.hpp b/include/teqp/models/association/association.hpp index c0b11c77..d38b134d 100644 --- a/include/teqp/models/association/association.hpp +++ b/include/teqp/models/association/association.hpp @@ -169,6 +169,9 @@ class Association{ const std::variant datasidecar; Association(const Eigen::ArrayXd& b_m3mol, const Eigen::ArrayXd& beta, const Eigen::ArrayXd& epsilon_Jmol, const std::vector>& molecule_sites, const AssociationOptions& options) : options(options), mapper(make_mapper(molecule_sites, options)), D(make_D(mapper, options)), m_Delta_rule(options.Delta_rule), datasidecar(CanonicalData{b_m3mol, beta, epsilon_Jmol, options.radial_dist}){ + if (options.Delta_rule != Delta_rules::CR1){ + throw std::invalid_argument("Delta rule is invalid; options are: {CR1}"); + } } Association(const decltype(datasidecar)& data, const std::vector>& molecule_sites, const AssociationOptions& options) : options(options), mapper(make_mapper(molecule_sites, options)), D(make_D(mapper, options)), m_Delta_rule(options.Delta_rule), datasidecar(data) { } From cad49193e9ffb09d9fe6d7648233520fbc7ffff2 Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Thu, 4 Jul 2024 11:30:44 -0400 Subject: [PATCH 12/13] Fix interaction partners in generalized association factory --- .../teqp/models/association/association.hpp | 47 ++++++++- src/tests/catch_test_SAFTgeneric.cxx | 96 ++++++++++++------- 2 files changed, 105 insertions(+), 38 deletions(-) diff --git a/include/teqp/models/association/association.hpp b/include/teqp/models/association/association.hpp index d38b134d..c4686fe6 100644 --- a/include/teqp/models/association/association.hpp +++ b/include/teqp/models/association/association.hpp @@ -176,6 +176,46 @@ class Association{ Association(const decltype(datasidecar)& data, const std::vector>& molecule_sites, const AssociationOptions& options) : options(options), mapper(make_mapper(molecule_sites, options)), D(make_D(mapper, options)), m_Delta_rule(options.Delta_rule), datasidecar(data) { } static Association factory(const nlohmann::json& j){ + + // Collect the set of unique site types among all the molecules + std::set unique_site_types; + for (auto molsite : j.at("molecule_sites")) { + for (auto & s : molsite){ + unique_site_types.insert(s); + } + } + + auto get_interaction_partners = [&](const nlohmann::json& j){ + std::map> interaction_partners; + + if (j.contains("options") && j.at("options").contains("interaction_partners")){ + interaction_partners = j.at("options").at("interaction_partners"); + for (auto [k,partners] : interaction_partners){ + if (unique_site_types.count(k) == 0){ + throw teqp::InvalidArgument("Site is invalid in interaction_partners: " + k); + } + for (auto& partner : partners){ + if (unique_site_types.count(partner) == 0){ + throw teqp::InvalidArgument("Partner " + partner + " is invalid for site " + k); + } + } + } + } + else{ + // Every site type is assumed to interact with every other site type, except for itself + for (auto& site1 : unique_site_types){ + std::vector partners; + for (auto& site2: unique_site_types){ + if (site1 != site2){ + partners.push_back(site2); + } + } + interaction_partners[site1] = partners; + } + } + return interaction_partners; + }; + if (j.contains("Delta_rule")){ std::string Delta_rule = j.at("Delta_rule"); if (Delta_rule == "CR1"){ @@ -186,6 +226,7 @@ class Association{ auto options = get_association_options(j); options.Delta_rule = Delta_rules::CR1; data.radial_dist = options.radial_dist; + options.interaction_partners = get_interaction_partners(j); return {data, j.at("molecule_sites"), options}; } else if(Delta_rule == "Dufal"){ @@ -199,15 +240,13 @@ class Association{ // Parameters for the associating part data.epsilon_HB_Jmol = toEig(j, "epsilon_HB / J/mol"); data.K_HB_m3 = toEig(j, "K_HB / m^3"); + data.apply_mixing_rules(); auto options = get_association_options(j); options.Delta_rule = Delta_rules::Dufal; - data.apply_mixing_rules(); + options.interaction_partners = get_interaction_partners(j); return {data, j.at("molecule_sites"), options}; } - else{ - - } } throw std::invalid_argument("The Delta_rule has not been specified"); } diff --git a/src/tests/catch_test_SAFTgeneric.cxx b/src/tests/catch_test_SAFTgeneric.cxx index 34a31570..8e829156 100644 --- a/src/tests/catch_test_SAFTgeneric.cxx +++ b/src/tests/catch_test_SAFTgeneric.cxx @@ -5,7 +5,7 @@ using Catch::Approx; #include "teqp/cpp/teqpcpp.hpp" #include "teqp/models/saft/genericsaft.hpp" - +#include "teqp/cpp/deriv_adapter.hpp" using namespace teqp; TEST_CASE("Benchmark generic PC-SAFT+Association model", "[SAFTgeneric]"){ @@ -95,44 +95,46 @@ TEST_CASE("Benchmark generic PC-SAFT+Association model", "[SAFTgeneric]"){ }; } -TEST_CASE("Benchmark Dufal water model", "[SAFTgeneric]"){ - auto contents = R"( - { - "nonpolar": { - "kind": "SAFT-VR-Mie", - "model": { - "coeffs": [ - { - "name": "Water", - "BibTeXKey": "Dufal-2015", - "m": 1.0, - "sigma_Angstrom": 3.0555, - "epsilon_over_k": 418.00, - "lambda_r": 35.823, - "lambda_a": 6.0 - } - ] - } - }, - "association": { - "kind": "Dufal", - "model": { - "sigma / m": [3.0555e-10], - "epsilon / J/mol": [3475.445374388054], - "lambda_r": [35.823], - "epsilon_HB / J/mol": [13303.140189045183], - "K_HB / m^3": [496.66e-30], - "kmat": [[1.0]], - "Delta_rule": "Dufal", - "molecule_sites": [["e","e","H","H"]] + +auto Dufal_contents = R"( +{ + "nonpolar": { + "kind": "SAFT-VR-Mie", + "model": { + "coeffs": [ + { + "name": "Water", + "BibTeXKey": "Dufal-2015", + "m": 1.0, + "sigma_Angstrom": 3.0555, + "epsilon_over_k": 418.00, + "lambda_r": 35.823, + "lambda_a": 6.0 } - } + ] } - )"_json; + }, + "association": { + "kind": "Dufal", + "model": { + "sigma / m": [3.0555e-10], + "epsilon / J/mol": [3475.445374388054], + "lambda_r": [35.823], + "epsilon_HB / J/mol": [13303.140189045183], + "K_HB / m^3": [496.66e-30], + "kmat": [[0.0]], + "Delta_rule": "Dufal", + "molecule_sites": [["e","e","H","H"]] + } + } +} +)"_json; +TEST_CASE("Benchmark Dufal water model", "[SAFTgeneric]"){ + nlohmann::json maincontents = { {"kind", "genericSAFT"}, - {"model", contents} + {"model", Dufal_contents} }; auto model = teqp::cppinterface::make_model(maincontents); auto z = (Eigen::ArrayXd(1) << 1.0).finished(); @@ -150,3 +152,29 @@ TEST_CASE("Benchmark Dufal water model", "[SAFTgeneric]"){ return model->get_Ar20(300, 10000, z); }; } + +TEST_CASE("Check Dufal water model", "[Dufaldetail]"){ + + nlohmann::json maincontents = { + {"kind", "genericSAFT"}, + {"model", Dufal_contents} + }; + auto model = teqp::cppinterface::make_model(maincontents); + auto z = (Eigen::ArrayXd(1) << 1.0).finished(); + + using namespace teqp::cppinterface::adapter; + const auto &gen = get_model_cref(model.get()); + const auto& a = std::get(gen.association.value()); + + // Check a result from Dufal for the non-bonded site-fraction + double rhomolar = 1.011*100*100*100/1000/0.018015; + Eigen::ArrayXd X_init_bad(40); X_init_bad = 1.0; + CHECK_THROWS(a.successive_substitution(290.0, rhomolar, z, X_init_bad)); + +// std::cout << a.D << std::endl; +// std::cout << a.get_Delta(300.0, 1/3e-5, z) << std::endl; + + Eigen::ArrayXd X_init(2); X_init = 1.0; + auto X = a.successive_substitution(290.0, rhomolar, z, X_init); + CHECK(0.08922902098124778 == Approx(X[0])); // Check value from Clapeyron +} From 29d8f4b166b21f2d59dd7ba0b1e4aef1f1342c34 Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Thu, 4 Jul 2024 11:31:05 -0400 Subject: [PATCH 13/13] Check size of X_init just to be sure --- include/teqp/models/association/association.hpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/include/teqp/models/association/association.hpp b/include/teqp/models/association/association.hpp index c4686fe6..b3a5942c 100644 --- a/include/teqp/models/association/association.hpp +++ b/include/teqp/models/association/association.hpp @@ -355,6 +355,10 @@ class Association{ template auto successive_substitution(const TType& T, const RhoType& rhomolar, const MoleFracsType& molefracs, const XType& X_init) const { + if (X_init.size() != static_cast(mapper.to_siteid.size())){ + throw teqp::InvalidArgument("Wrong size of X_init; should be "+ std::to_string(mapper.to_siteid.size())); + } + using resulttype = std::common_type_t; // Type promotion, without the const-ness using Mat = Eigen::Array;