diff --git a/doc/source/models/Association.ipynb b/doc/source/models/Association.ipynb
index 2ef6341..633318c 100644
--- a/doc/source/models/Association.ipynb
+++ b/doc/source/models/Association.ipynb
@@ -7,7 +7,7 @@
"source": [
"# Association\n",
"\n",
- "\n",
+ "\n",
"\n",
"## Site Interactions\n",
"\n",
@@ -70,6 +70,28 @@
"$$\n",
"This is the method utilized in Langenbach and Enders.\n",
"\n",
+ "## Simplified analysis for pure fluids\n",
+ "\n",
+ "In the case that there is only one non-zero interaction strength, the mathematics can be greatly simplified. This section is based on [the derivations of Pierre Walker]( https://github.com/ClapeyronThermo/Clapeyron.jl/discussions/260#discussioncomment-8572132). The general result for a pure fluid with N types of sites, which looks like\n",
+ "\n",
+ "is as shown in [Eq. A39 in Lafitte et al.](https://doi.org/10.1063/1.4819786):\n",
+ "\n",
+ "$$ \n",
+ "X_{ai} = \\dfrac{1}{1+\\rho_{\\rm N}\\displaystyle\\sum_{j=1}^nx_j\\displaystyle\\sum_{b=1}^{s_j}n_{b,j}X_{bj}\\Delta_{abij}}\n",
+ "$$\n",
+ "\n",
+ "If you have a pure fluid that has two types of sites of arbitrary multiplicity, and no site-site self association (meaning that A cannot dock with A and B cannot dock with B), you can write out the law of mass-action as\n",
+ "\n",
+ "$$X_A = \\frac{1}{1+\\rho_{\\rm N}n_BX_B\\Delta_{AB}} $$\n",
+ "$$X_B = \\frac{1}{1+\\rho_{\\rm N}n_AX_A\\Delta_{BA}} $$\n",
+ "\n",
+ "and further assuming that $\\Delta=\\Delta_{BA}=\\Delta_{AB}$ and the simplification of $\\kappa_k=\\rho_{\\rm N}n_k\\Delta$ yields\n",
+ "\n",
+ "$$X_A = \\frac{1}{1+\\kappa_BX_B} $$\n",
+ "$$X_B = \\frac{1}{1+\\kappa_AX_A} $$\n",
+ "\n",
+ "which can be solved simultaneously from a quadratic equation (see below)\n",
+ "\n",
"## Interaction strength\n",
"\n",
"The interaction site strength is a matrix with side length of the number of ``siteid``. It is a block matrix because practically speaking the interaction sites are still about molecule-molecule interactions\n",
@@ -103,6 +125,46 @@
"K. Langenbach & S. Enders (2012): Cross-association of multi-component systems, Molecular Physics, 110:11-12, 1249-1260; https://dx.doi.org/10.1080/00268976.2012.668963"
]
},
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "7768d099",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Here is the algebraic solutions to the laws of mass-action for the simplified cases\n",
+ "# covered in Huang & Radosz for the case of a pure fluid with two types \n",
+ "# of sites of arbitrary multiplicity\n",
+ "import sympy as sy\n",
+ "\n",
+ "rhoN, kappa_B, kappa_A, X_A, X_B = sy.symbols('rho_N, kappa_B, kappa_A, X_A, X_B')\n",
+ "\n",
+ "# Definitions of the equations to be solved\n",
+ "# In Eq(), the first arg is the LHS, second is RHS\n",
+ "eq1 = sy.Eq(X_B, 1/(1+kappa_A*X_A))\n",
+ "eq2 = sy.Eq(X_A, 1/(1+kappa_B*X_B))\n",
+ " \n",
+ "# The solutions\n",
+ "solns = sy.solve([eq1, eq2], [X_A, X_B])\n",
+ "for soln in solns:\n",
+ " for x in soln:\n",
+ " display(x)\n",
+ " \n",
+ "# 2B scheme; one site of type A, one site of type B, no A-A or B-B interactions\n",
+ "print('2B solutions:')\n",
+ "for soln in solns:\n",
+ " print('X_A, X_B:')\n",
+ " for x in soln:\n",
+ " display(x.subs(kappa_A,rhoN*1*Delta).subs(kappa_B,rhoN*1*Delta))\n",
+ " \n",
+ "# 3B scheme; one site of type A, two sites of type B, no A-A or B-B interactions\n",
+ "print('3B solutions:')\n",
+ "for soln in solns:\n",
+ " print('X_A, X_B:') \n",
+ " for x in soln:\n",
+ " display(x.subs(kappa_A,rhoN*1*Delta).subs(kappa_B,rhoN*2*Delta))"
+ ]
+ },
{
"cell_type": "code",
"execution_count": null,
diff --git a/doc/source/models/SingleABMolecule.drawio.svg b/doc/source/models/SingleABMolecule.drawio.svg
new file mode 100644
index 0000000..7b39355
--- /dev/null
+++ b/doc/source/models/SingleABMolecule.drawio.svg
@@ -0,0 +1,4 @@
+
+
+
+
\ No newline at end of file
diff --git a/include/teqp/models/association/association.hpp b/include/teqp/models/association/association.hpp
index 3d06be3..9ad092d 100644
--- a/include/teqp/models/association/association.hpp
+++ b/include/teqp/models/association/association.hpp
@@ -30,6 +30,7 @@ struct AssociationOptions{
association::radial_dists radial_dist;
association::Delta_rules Delta_rule = association::Delta_rules::CR1;
std::vector self_association_mask;
+ bool allow_explicit_fractions=true;
double alpha = 0.5;
double rtol = 1e-12, atol = 1e-12;
int max_iters = 100;
@@ -380,13 +381,35 @@ class Association{
}
using rDDXtype = std::decay_t>; // Type promotion, without the const-ness
+ Eigen::ArrayX> X = X_init.template cast(), Xnew;
+
Eigen::MatrixX rDDX = rhomolar*N_A*(Delta.array()*D.cast().array()).matrix();
for (auto j = 0; j < rDDX.rows(); ++j){
rDDX.row(j).array() = rDDX.row(j).array()*xj.array().template cast();
}
// rDDX.rowwise() *= xj;
- Eigen::ArrayX> X = X_init.template cast(), Xnew;
+ // Use explicit solutions in the case that there is a pure
+ // fluid with two kinds of sites, and no self-self interactions
+ // between sites
+ if (options.allow_explicit_fractions && molefracs.size() == 1 && mapper.counts.size() == 2 && (rDDX.matrix().diagonal().unaryExpr([](const auto&x){return getbaseval(x); }).array() == 0.0).all()){
+ auto Delta_ = Delta(0, 1);
+ auto kappa_A = rhomolar*N_A*static_cast(mapper.counts[0])*Delta_;
+ auto kappa_B = rhomolar*N_A*static_cast(mapper.counts[1])*Delta_;
+ // See the derivation in the docs in the association page; see also https://github.com/ClapeyronThermo/Clapeyron.jl/blob/494a75e8a2093a4b48ca54b872ff77428a780bb6/src/models/SAFT/association.jl#L463
+ auto X_A1 = (kappa_A-kappa_B-sqrt(kappa_A*kappa_A-2.0*kappa_A*kappa_B + 2.0*kappa_A + kappa_B*kappa_B + 2.0*kappa_B+1.0)-1.0)/(2.0*kappa_A);
+ auto X_A2 = (kappa_A-kappa_B+sqrt(kappa_A*kappa_A-2.0*kappa_A*kappa_B + 2.0*kappa_A + kappa_B*kappa_B + 2.0*kappa_B+1.0)-1.0)/(2.0*kappa_A);
+ // Keep the positive solution, likely to be X_A2
+ if (getbaseval(X_A1) < 0 && getbaseval(X_A2) > 0){
+ X(0) = X_A2;
+ }
+ else if (getbaseval(X_A1) > 0 && getbaseval(X_A2) < 0){
+ X(0) = X_A1;
+ }
+ auto X_B = 1.0/(1.0+kappa_A*X(0)); // From the law of mass-action
+ X(1) = X_B;
+ return X;
+ }
for (auto counter = 0; counter < options.max_iters; ++counter){
// calculate the new array of non-bonded site fractions X
diff --git a/src/tests/catch_test_association.cxx b/src/tests/catch_test_association.cxx
index 3fe4e9b..73827d7 100644
--- a/src/tests/catch_test_association.cxx
+++ b/src/tests/catch_test_association.cxx
@@ -237,3 +237,45 @@ TEST_CASE("Benchmark association evaluations", "[associationbench]"){
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;
}
+
+TEST_CASE("Check explicit solutions for association fractions from old and new code","[XA]"){
+ double T = 298.15, rhomolar = 1000/0.01813;
+ double epsABi = 16655.0, betaABi = 0.0692, bcubic = 0.0000145, RT = 8.31446261815324*T;
+ auto molefrac = (Eigen::ArrayXd(1) << 1.0).finished();
+
+ // Explicit solution from Huang & Radosz (old-school method)
+ auto X_Huang = teqp::CPA::XA_calc_pure(4, teqp::CPA::association_classes::a4C, teqp::CPA::radial_dist::CS, epsABi, betaABi, bcubic, RT, rhomolar, molefrac);
+
+ auto b_m3mol = (Eigen::ArrayXd(1) << 0.0145/1e3).finished();
+ auto beta = (Eigen::ArrayXd(1) << 69.2e-3).finished();
+ auto epsilon_Jmol = (Eigen::ArrayXd(1) << 166.55*100).finished();
+
+ std::vector> molecule_sites = {{"e", "e", "H", "H"}};
+ association::AssociationOptions opt;
+ opt.radial_dist = association::radial_dists::CS;
+ opt.max_iters = 1000;
+ opt.allow_explicit_fractions = true;
+ opt.interaction_partners = {{"e", {"H",}}, {"H", {"e",}}};
+ association::Association aexplicit(b_m3mol, beta, epsilon_Jmol, molecule_sites, opt);
+
+ opt.allow_explicit_fractions = false;
+ association::Association anotexplicit(b_m3mol, beta, epsilon_Jmol, molecule_sites, opt);
+
+ Eigen::ArrayXd X_init = Eigen::ArrayXd::Ones(2);
+
+ // Could be short-circuited solution from new derivation, and should be equal to Huang & Radosz
+ auto X_newderiv = aexplicit.successive_substitution(T, rhomolar, molefrac, X_init);
+
+ // Force the iterative routines to be used as a further sanity check
+ auto X_newderiv_iterative = anotexplicit.successive_substitution(T, rhomolar, molefrac, X_init);
+
+ CHECK_THAT(X_Huang(0), WithinRel(X_newderiv(0), 1e-10));
+ CHECK_THAT(X_Huang(0), WithinRel(X_newderiv_iterative(0), 1e-10));
+
+ BENCHMARK("Calling explicit solution"){
+ return aexplicit.successive_substitution(T, rhomolar, molefrac, X_init);
+ };
+ BENCHMARK("Calling non-explicit solution"){
+ return anotexplicit.successive_substitution(T, rhomolar, molefrac, X_init);
+ };
+}