Add explicit solutions for simple association in pure fluids

Additional generalizations are possible, but this is a good starting point
closes #137

doc/source/models/Association.ipynb

@@ -7,7 +7,7 @@
    "source": [
     "# Association\n",
     "\n",
-    "![Two example molecules, each with multiple assocation sites](Molecule.drawio.svg)\n",
+    "![Two example molecules, each with multiple association sites](Molecule.drawio.svg)\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",
+    "![One molecule, each with two association sites](SingleABMolecule.drawio.svg)\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,

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<bool> 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<std::common_type_t<typename decltype(Delta)::Scalar, decltype(rhomolar), decltype(molefracs[0])>>; // Type promotion, without the const-ness
+        Eigen::ArrayX<std::decay_t<rDDXtype>> X = X_init.template cast<rDDXtype>(), Xnew;
+
         Eigen::MatrixX<rDDXtype> rDDX = rhomolar*N_A*(Delta.array()*D.cast<resulttype>().array()).matrix();
         for (auto j = 0; j < rDDX.rows(); ++j){
             rDDX.row(j).array() = rDDX.row(j).array()*xj.array().template cast<rDDXtype>();
         }
 //        rDDX.rowwise() *= xj;
 
-        Eigen::ArrayX<std::decay_t<rDDXtype>> X = X_init.template cast<rDDXtype>(), 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<double>(mapper.counts[0])*Delta_;
+            auto kappa_B = rhomolar*N_A*static_cast<double>(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

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<std::vector<std::string>> 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);
+    };
+}