Possible improvements on association implementation

[ianhbell]

With a colleague we are working on adding association to teqp. As you have implemented this into Clapeyron.jl already, we plan to mirror your implementation quite closely since we can easily test against it in julia and jupyter.

Before we begin, are there any places there where you wish you had done something different from the beginning so that we can avoid needing to learn the same hard lessons as you? Is it sufficient to have only electron donors and acceptors in the full generality? I guess perhaps one could consider to do something similar with halogen bonding, at least in principle?

[longemen3000]

I have some notes here: #70

Some notes:

• the sparse solver was deprecated, accessing an multiplying a dense matrix is faster
  -there are some cases where sigma is T-dependent. (e-PCSAFT) In those cases, we can't do the combining beforehand.
• there are some papers that use the association as a way to specify for chemical reactions (SAFTγMie), I don't know if just having acceptors and donors would suffice to reproduce that?
• Until really recently (around 2 versions ago) we did allocate the Vector of non bonded fractions for a single component. We now, skip that and evaluate the results directly, resulting in a zero-allocation PCSAFT with assoc

In things I wish to implement someday:

• a way to do the Newton iterations without additional allocations (reusing the caches maybe?) it can converge faster, but it is unstable if a good initial point is not provided. feos does this, but the work there is easier, as they only have acceptor, donor and extra (A, B, C) site types. They also don't create the Vector of non bonded fractions (Output the fraction of non-bonded association sites feos-org/feos#210). We do, and use it for some electrolyte EOS, but I still have to check if that makes a difference.
• the exact solution of the problem when there are three site pairs (that covers a binary with mixing)
• a chain rule that implements implicit differenciation to solve the derivatives

[ianhbell]

Thank very much @longemen3000 for your insights, it is really helpful!

[longemen3000]

there is also this paper, https://pubs.acs.org/doi/10.1021/acs.iecr.2c02723 that has an alternative way (and by the looks of it, a faster one) to solve association, but that requires having the association strength defined for each site instead of defined for each site pair. if there was a way to convert between those, that would speed up the calculations a lot

[ianhbell]

@longemen3000 We wrote a little python script to do the D matrix construction of Langenbach. It took us a lot of time, it is not so obvious, but it also in the end is quite straightforward to understand once you know exactly what you are doing. And you can do it quite generically, with as many sites as you like. I think this is what we will implement in teqp.

[longemen3000]

i was just reading what i wrote some time ago, and i did, indeed cite the wrong paper 😓 , we don't use the michaelsen matrix, but a variation of the Langenbach one (i say variation because i did the work on compressing the indices before reading that paper). if you look at the example of 2B + 3B, we produce a permutation of the columns (1,2,3,4) -> (3,4,1,2)
Langenbach:

Clapeyron:

julia> model = SAFTVRMie(["ethanol","methanol"],assoc_options = AssocOptions(combining = :elliott_runtime))
SAFTVRMie{BasicIdeal, Float64} with 2 components:
 "ethanol"
 "methanol"
Contains parameters: Mw, segment, sigma, lambda_a, lambda_r, epsilon, epsilon_assoc, bondvol

julia> model.sites
SiteParam with 2 components:
 "ethanol": "e" => 1, "H" => 1
 "methanol": "e" => 2, "H" => 1

julia> Clapeyron.assoc_site_matrix(model,0.03,373.15,[0.4,0.6])
4×4 Matrix{Float64}:
 0.0         1.13795e26  0.0         1.08907e12
 1.13795e26  0.0         2.17814e12  0.0
 0.0         7.26046e11  0.0         0.00694862
 7.26046e11  0.0         0.0138972   0.0

i remember trying the pure initial point, but i saw no improvement on multicomponent systems with strong cross-interactions

[ianhbell]

Ok, good, I think we were intending to try that with the pure fluid association starting values, but you think we should just start with association fractions of zero?

[longemen3000]

the theory of quadratic matrix equations says that if you do succesive substitution, you should start at 1, the initialization i do is here:

Clapeyron.jl/src/models/SAFT/association.jl

Lines 291 to 295 in 2a3d6ea

	Kmin,Kmax = nonzero_extrema(K) #look for 0 < Amin < Amax
	if Kmax > 1
	f = true/Kmin
	else
	f = true-Kmin

where K (the entries of our matrix) are:

# V = volume (m3)
# z = mol amount (mol)
# Δ = association energy at the calculated index
# n = amount of sites of molecule j in position b
Kij = (N_A/V)*njb*z[j]*Δ[idx]

[ianhbell]

So you think the simple SS method in Langenbach (Eq. 19) is too simple? That was going to be my first thing to try

[longemen3000]

no, i don't think so (we also use SS with dampening), just a different starting point. in theory (quadratic matrix equation theory), any sucessive substitution method to solve this system is globally convergent if x0 > xsol, so the "best" starting point is 1. the function above is just a compromise to get closer to the real solution.

[ianhbell]

Do you have a reference I can cite on the global stability of damped SS from quadratic matrix equation theory?

[longemen3000]

https://doi.org/10.1016/j.laa.2011.05.036

[ianhbell]

Where are the routines in Clapeyron for association calculations? I poked around and had a hard time finding them.

[longemen3000]

• the association calculation occurs here:
  https://github.com/ClapeyronThermo/Clapeyron.jl/blob/master/src/models/SAFT/association.jl
  Relevant structures:
• association parameters: https://github.com/ClapeyronThermo/Clapeyron.jl/blob/master/src/database/params/AssocParam.jl
• amount of sites parameter: https://github.com/ClapeyronThermo/Clapeyron.jl/blob/master/src/database/params/SiteParam.jl
• struct that holds the association data: https://github.com/ClapeyronThermo/Clapeyron.jl/blob/master/src/database/params/paramvectors.jl

[ianhbell]

Our python code for the Langenbach method (the other methods not included) looks like:

interaction_partners = {
    'B': ('N', 'P', 'B'),
    'N': ('P', 'B'),
    'P': ('N', 'B')
}
        
def get_DIJ(I, J):
    """ 
    Return the value of an entry in the D_{IJ} matrix 
    """
    i, typei = inv_mapping[I]
    j, typej = inv_mapping[J]
    if typej in interaction_partners[typei]:
        return counts[(j, typej)]
    return 0

Operating with strings as site keys is no downside since this only happens at model instantiation and thereafter everything is with integers

[longemen3000]

Just an idea i'm working on. technically solving association is solving a system of polynomial equations. An homotopy continuation method could be used to solve the problem, allowing a fast newton solver to be used while having the good convergence properties of the damped SS method. a proof of concept:

julia> using Clapeyron, HomotopyContinuation

julia> model = PCSAFT(["water","ethanol"],assoc_options = AssocOptions(combining = :cr1))      
PCSAFT{BasicIdeal, Float64} with 2 components:
 "water"
 "ethanol"
Contains parameters: Mw, segment, sigma, epsilon, epsilon_assoc, bondvol

julia> A = Clapeyron.assoc_site_matrix(model,0.03,373.15,[0.4,0.6])
4×4 Matrix{Float64}:
 0.0         0.0061563   0.0        0.0118999
 0.0061563   0.0         0.0118999  0.0
 0.0         0.00793324  0.0        0.015334
 0.00793324  0.0         0.015334   0.0

julia> X = Clapeyron.X(model,0.03,373.15,[0.4,0.6]).v #solves association, returns raw vector
4-element Vector{Float64}:
 0.982623238409625
 0.982623238409625
 0.9777200008431199
 0.9777200008431199

julia> @var x1 x2 x3 x4 #define variables
(x1, x2, x3, x4)

julia> x = [x1,x2,x3,x4]
4-element Vector{Variable}:
 x1
 x2
 x3
 x4

julia> system = System(A*x .* x + x - ones(Int,4)) #quadratic polynomial system to solve
System of length 4
 4 variables: x1, x2, x3, x4

 -1 + x1 + x1*(0.0061563003744746*x2 + 0.0118998581866071*x4)
 -1 + x2 + x2*(0.0061563003744746*x1 + 0.0118998581866071*x3)
 -1 + x3 + x3*(0.00793323879107139*x2 + 0.0153339652426488*x4)
 -1 + x4 + x4*(0.00793323879107139*x1 + 0.0153339652426488*x3)

julia> result = solve(system) #solve system
Result with 6 solutions
=======================
• 6 paths tracked
• 6 non-singular solutions (6 real)
• random_seed: 0x8fc22b3a
• start_system: :polyhedral

julia> real_solutions(result)
6-element Vector{Vector{Float64}}:
 [4.649456560860215, -121.03902836969094, -87.13425260084185, -3.3419293138077375]
 [0.9826232384090261, 0.9826232384090261, 0.9777200008421453, 0.9777200008421454]
 [4.70277979831464, 4.702779798314639, -68.59842529822113, -68.59842529822114]
 [875900.535060796, 875900.535060796, -453224.4591725971, -453224.4591725971]
 [-156.9584505920934, -156.9584505920934, -3.368758186430234, -3.3687581864302336]
 [-121.03902836969094, 4.649456560860218, -3.3419293138077384, -87.13425260084185]

julia> real_solutions(result)[2] - X
4-element Vector{Float64}:
 -5.988542994828094e-13
 -5.988542994828094e-13
 -9.74664793318425e-13
 -9.745537710159624e-13

[longemen3000]

298 K, 1 atm? it does find it:

julia> model = SAFTgammaMie(["water","ethanol","methanol","carbon dioxide","MEA" => ["CH3" => 1, "CH2" => 1, "NH2" => 1]],assoc_options = AssocOptions(max_iters = 10000)) #fails with 1000 iters
SAFTgammaMie{BasicIdeal, SAFTVRMie{BasicIdeal, Float64}} with 5 components:
 "water": "H2O" => 1
 "ethanol": "CH2OH" => 1, "CH3" => 1
 "methanol": "CH3OH" => 1
 "carbon dioxide": "CO2" => 1
 "MEA": "CH3" => 1, "CH2" => 1, "NH2" => 1
Group Type: SAFTgammaMie
Contains parameters: segment, shapefactor, lambda_a, lambda_r, sigma, epsilon, epsilon_assoc, bondvol

julia> p = 1.0e5
100000.0

julia> T = 298.15
298.15

julia> z = Clapeyron.Fractions.zeros(5)
5-element Vector{Float64}:
 0.2
 0.2
 0.2
 0.2
 0.2

julia> v = volume(model,p,T,z)
4.2610813888613514e-5

julia> X = Clapeyron.X(model.vrmodel,v,T,z).v
10-element Vector{Float64}:
 0.05329649748564787
 0.27035697875231646
 0.2657359341269344
 0.21998827772804902
 0.08006041473790101
 0.3895574718713385
 0.008677421270164014
 6.560123096776866e-5
 0.15184119295721737
 0.0022790254927238735

julia> A = Clapeyron.assoc_site_matrix(model.vrmodel,v,T,z)
10×10 Matrix{Float64}:
  0.0      14.7556    0.0       35.7803    0.0      …      0.0       0.0            2.21181    
 14.7556    0.0       0.265658   0.0       7.56935         0.959134  8.13831        0.0        
  0.0       0.531317  0.0       11.8324    0.0             0.0       0.0            7.23964    
 35.7803    0.0       5.9162     0.0       0.0             0.0       0.438504       0.0        
  0.0      15.1387    0.0        0.0       0.0             0.0       0.0            0.0        
 15.1387    0.0       0.0        0.0       9.49502  …      0.0       0.0            0.0        
  0.0       0.0       0.0        0.0       0.0             0.0       0.0        50127.4        
  0.0       1.91827   0.0        0.0       0.0             0.0       1.00255e5   8464.57
  0.0       8.13831   0.0        0.438504  0.0         50127.4       0.0            0.305033   
  4.42361   0.0       7.23964    0.0       0.0          8464.57      0.610067       0.0        

julia> @var xx[1:10]
(Variable[xx₁, xx₂, xx₃, xx₄, xx₅, xx₆, xx₇, xx₈, xx₉, xx₁₀],)

julia> system = System(A*xx .* xx + xx - ones(Int,10))
System of length 10
 10 variables: xx₁, xx₂, xx₃, xx₄, xx₅, xx₆, xx₇, xx₈, xx₉, xx₁₀

 -1 + xx₁ + xx₁*(2.21180577905624*xx₁₀ + 14.755555412613*xx₂ + 35.7803360391896*xx₄ + 15.138706954092*xx₆)
 -1 + xx₂ + xx₂*(14.755555412613*xx₁ + 0.265658354440691*xx₃ + 7.56935347704601*xx₅ + 0.959133909338137*xx₈ + 8.1383083573105*xx₉)
 -1 + xx₃ + xx₃*(7.23963548945751*xx₁₀ + 0.531316708881382*xx₂ + 11.8324021672672*xx₄)
 -1 + xx₄ + xx₄*(35.7803360391896*xx₁ + 5.91620108363358*xx₃ + 0.438503818216114*xx₉)
 -1 + xx₅ + (15.138706954092*xx₂ + 18.9900407011838*xx₆)*xx₅
 -1 + xx₆ + (15.138706954092*xx₁ + 9.49502035059192*xx₅)*xx₆
 -1 + xx₇ + 50127.3946214237*xx₇*xx₁₀
 -1 + xx₈ + xx₈*(8464.57443725259*xx₁₀ + 1.91826781867627*xx₂ + 100254.789242847*xx₉)
 -1 + xx₉ + xx₉*(0.305033489335745*xx₁₀ + 8.1383083573105*xx₂ + 0.438503818216114*xx₄ + 50127.3946214237*xx₈)
 -1 + xx₁₀ + (4.42361155811247*xx₁ + 7.23963548945751*xx₃ + 50127.3946214237*xx₇ + 8464.57443725259*xx₈ + 0.61006697867149*xx₉)*xx₁₀

julia> result = solve(system)
Result with 224 solutions
=========================
• 224 paths tracked
• 224 non-singular solutions (50 real)
• random_seed: 0xd0d25c3a
• start_system: :polyhedral


julia> function is_valid_assoc(x)
           x0,x1 = extrema(x)
               return (x0 >= 0) & (x1 <= 1)
               end
is_valid_assoc (generic function with 1 method)

julia> findall(is_valid_assoc,real_solutions(result))
1-element Vector{Int64}:
 38

julia> real_solutions(result)[38] - X
10-element Vector{Float64}:
  2.5571281203617957e-12
 -2.6569857425329246e-12
  4.387179508569261e-11
 -1.7078310987628242e-11
  1.2022466355787742e-12
 -7.673472968150463e-12
  2.1689158455351354e-10
  3.2940501454997717e-15
 -2.7590985052228234e-12
 -5.696502341015486e-11

[pw0908]

Wow, that's impressive

[ianhbell]

I want to loop in @gustavochm on this thread since he is using SS plus Michelsen's second order method in sgtpy and might be interested in following this discussion

[ianhbell]

How long did the homotopy solution take?

[longemen3000]

At them moment Im not able to make a adequate comparison. For the last example, the homotopy solver found 224 solutions (so at least 224 times slower) , whereas our association solver only returns one (the one that we need). I will report back if I can handcode a homotopy solver that only traverses to the expected result, instead of all of them.

[longemen3000]

Some updates:

Michelsen was right all along, the newton minimization method is faster, even with the increased allocations (1 matrix, one vector of integers to perform LU, two support vectors). LU factorization is necessary, gauss-seidel or jacobi iterations to solve the linear system eat away the speed-up.

on successive substitution, one way to accelerate the iterations is to reuse the results in each iteration, instead of:

#SS
function f_ss!(out,K,x,w)
  Kx = K*x
  out .= w ./ (1 .+ Kx) .+ (1 - w) .* x
  return out
end

we can do something like:

#SS with acceleration
function f_ss2!(out,K,x,w)
    n = length(x)
    i_solved = 0
    for i in 1:n
        kxi = zero(eltype(x))
        #strategy 3
        Ki = @view K[i,:]
        for j in 1:i_solved
            kxi += out[j]*Ki[j] #we reuse the already calculated kxi
        end
        for j in (i_solved+1):n
            kxi += x[j]*Ki[j]
        end
        out[i] = w / (1 + kxi ) + (1 - w) * x[i]
        i_solved += 1
    end
end

This is similar to one technique used to accelerate COSMOSAC iterative calculations.

on the AD part, there are two ways to propagate derivative information.

• if you store the vector of association strengths (including derivative information), you can directly use michelsen's Q function to get derivatives of a_assoc.
• if you really need derivatives of the nonbonded fractions, you can perform an additional newton step ("perfect newton step"), but with a hessian and gradient that propagates derivative information. Some electrolyte EoS need nonbonded fractions with derivative propagation.

[longemen3000]

no reductions are done in the association solver, but it seems to give results even if one fraction is zero:

#association: Δijab = sqrt(Δiiab * Δjjab)
julia> model = PCSAFT(["water","ethanol"],assoc_options = AssocOptions(combining = :elliott_runtime))
PCSAFT{BasicIdeal, Float64} with 2 components:
 "water"
 "ethanol"
Contains parameters: Mw, segment, sigma, epsilon, epsilon_assoc, bondvol

#association: pure ethanol, 2 component model
julia> Clapeyron.X(model,v,T,[0.0,1.0])
2-element pack(::Vector{Vector{Float64}}):
 [0.04748238983368363, 0.04748238983368367]
 [0.03477405963556611, 0.03477405963556614]

#pure ethanol model
julia> model1 = split_model(model)[2]
PCSAFT{BasicIdeal, Float64} with 1 component:
 "ethanol"
Contains parameters: Mw, segment, sigma, epsilon, epsilon_assoc, bondvol

#association, pure ethanol
julia> Clapeyron.X(model1,v,T,[1.0])
1-element pack(::Vector{Vector{Float64}}):
 [0.03477405963556612, 0.034774059635566124]

maybe the Langenbach does not suffer that specific problem?

[gustavochm]

Yeah, I'm not quite sure 😅; at least in sgtpy, I had to change to solving directly using the normal Newton method (not the modified Michelsen optimization) and did the trick when one of the molar fractions was zero.

[longemen3000]

Ohhh, that seems to be the case then 😅, Im using your minimization method, so it handles those zero fractions in the same way.
On handling the zero fractions, I think that you can still do the reduction, resizing X and K, and expand again at the end, but it is a niche case

You could also deactivate the hessian and gradient at the zero fraction indices (H[i, i] = 1, H[:, i] .= 0,H[i, :] .= 0, g[i] = 0)

[ianhbell]

Just a warning: having been burned by this before in other contexts, a binary mixture with a mole fraction of 1.0 of one component is NOT the same thing as the pure component. This is because the infinite dilution derivatives are meaningless in the pure fluid case, but quite relevant for the mixture. I'm not sure that causes trouble in this case, but something to watch out for.

[longemen3000]

haha, i know about that problem, i was facing the same at the start when using mixing rules (the derivatives of sqrt(0), for example). on this specific formulation, the michelsen term has a problem on association has a problem, but i think i can handle this specific case if perform L'hopital (via AD)