Diameters in SAFT-VR-Mie appear to be incorrect

[ianhbell]

EDIT I think I calculated diameters at the wrong temperatures; too many moving pieces

TL;DR: The diameters for SAFT-VR-Mie are calculated in a weird way and don't seem to give the right values

I have nearly finished my SAFT-VR-Mie implementation in teqp : https://github.com/usnistgov/teqp/tree/SAFTVRMie for mixtures (but no association). In so doing, I have checked all intermediate calculations (mostly good to numerical precision), and by stepping into the Julia code, I have noticed a difference in the calculation of the diameter (see #150). I don't follow the approach taken in Clapeyron.jl to calculate the diameter, and indeed, it does not match the 5-point Gauss-Legendre quadrature integration. I picked a random temperature for the coefficients of ethane and calculated values of d:

Npoints value   error
3   3.7257000E-10   1.55%
5   3.7157102E-10   1.27%
7   3.6794295E-10   0.29%
15  3.6689926E-10   0.00%
30  3.6689376E-10   0.00%

exact   3.6689376E-10   

Clapeyron.jl    3.5978393E-10   -1.94%

So Clapeyron.jl ends up with a value that is not in the convergence path of the calculation with fixed quadratures. The code for "exact" is below.

This is not just inside baseball; this has knock on effects on all calculated properties. For instance the critical points are not the same as in the paper of Lafitte.

My Julia code:

using Clapeyron
model = SAFTVRMie(["Ethane"])
display(model.params.segment)
display(model.params.sigma)
display(model.params.epsilon)
display(model.params.lambda_r)
display(model.params.lambda_a)
display(Clapeyron.data(model, 1/12000.0, 300.0, [1.0]))
# (_d,_ρ_S,ζi,_ζ_X,_ζst,σ3x,m̄)
ahs = Clapeyron.a_hs(model, 1/12000.0, 300.0, [1.0])
amono = Clapeyron.a_mono(model, 1/12000.0, 300.0, [1.0])
ares = Clapeyron.a_chain(model, 1/12000.0, 300.0, [1.0])

And python code to calculate the exact d with adaptive quadrature

# Here is a small example of using adaptive quadrature 
# to obtain the quasi-exact value of d for ethane
import numpy as np 
import scipy.integrate

epskB = 206.12
sigma_m = 3.7257e-10
lambda_r = 12.4
lambda_a = 6.0
C = lambda_r/(lambda_r-lambda_a)*(lambda_r/lambda_a)**(lambda_a/(lambda_r-lambda_a))
T = 100.0 # [K]
def integrand(r_m):
    u = C*epskB*((sigma_m/r_m)**(lambda_r) - (sigma_m/r_m)**(lambda_a))
    return 1-np.exp(-u/T)
print('(value, error estimate):')
print(scipy.integrate.quad(integrand, 0.0, sigma_m))

yielding

(value, error estimate):
(3.668937640717724e-10, 1.437139875642782e-12)

[longemen3000]

equivalent julia code

using Cubature
epskB = 206.12
sigma_m = 3.7257e-10
lambda_r = 12.4
lambda_a = 6.0
C = lambda_r/(lambda_r-lambda_a)*(lambda_r/lambda_a)^(lambda_a/(lambda_r-lambda_a))
T = 300.0 # [K]
function integrand(r_m)
    u = C*epskB*((sigma_m/r_m)^(lambda_r) - (sigma_m/r_m)^(lambda_a))
    return -expm1(-u/T) #exp(x) - 1
end
Cubature.hquadrature(integrand,0.0,sigma_m) #alternatively, Clapeyron.Solvers.integral21(integrand,0.0,sigma_m)

seems weird?

#T = 300K
julia> Clapeyron.d(model,1/12000,300,Clapeyron.SA[1.0])
1-element StaticArraysCore.SVector{1, Float64} with indices SOneTo(1):
 3.597839296564643e-10

#T = 100K
julia> Clapeyron.d(model,1/12000,100,Clapeyron.SA[1.0])
1-element StaticArraysCore.SVector{1, Float64} with indices SOneTo(1):
 3.668937610504847e-10

[ianhbell]

Wait, maybe the problem is on my side. I think one half of the test is at 300 K, and the other half at 100 K.

[ianhbell]

Sorry for the noise folks, user error. Still very interested to understand the magic that Clapeyron.jl uses to calculate d; it is surprisingly accurate.

[longemen3000]

given:

  C = (λr/(λr-λa))*(λr/λa)^(λa/(λr-λa)) 
  u(r) = C*ϵ*((σ/r)^λr - (σ/r)^λa)
  f(r) = 1-exp(-u(r)/T)
  d = integral(f,0,σ)

1 is constant, so:

  g(r) = exp(-u(r)/T)
  d = σ - integral(g,0,σ)

we can do a change of variables, j = r/σ, dj/dr = σ

  h(j) = exp(C*ϵ*((1/j)^λr - (1/j)^λa)/T)
  h(j) = exp(C*ϵ*(j^-λr - j^-λa)/T)
  g(r) = exp(-u(r)/T)
  d = σ - σ*integral(g,0,1) = σ(1 - integral(h,0,1))

got stuck here, but i saw this comment on the SAFTgammaMie paper:

The integral of Eq. (10) is obtained by means of Gauss-Legendre quadrature, a technique previously employed by Paricaud[106] who showed that a 5-point Gauss-Legendre procedure is adequate for an accurate representation of the effective hard-sphere diameter dkk.

That reference got me to the real integrand used:

[ianhbell]

What is $x_{\infty}$?

[longemen3000]

from the paper: "x∞ is the reduced distance below which the potential is considered as infinite" . that x∞ also appears in the SAFT-VRQ-Mie paper, but it is used with the vanilla integral instead https://github.com/ClapeyronThermo/Clapeyron.jl/blob/062350298afd1e3c62e003b0e2d7bbe7258f4f27/src/models/SAFT/SAFTVRMie/variants/SAFTVRQMie.jl#L181-L203 the first VRQ paper has a nice diagram explaining that

[ianhbell]

That makes sense, you want to be putting in your integration quadrature where the function is actually contributing to the integral. Now you have to solve for where the function is zero in double precision. That implies a significant complexity since solving for the $r/\sigma$ where the integrand is one is a non-trivial exercise.

But this isn't what is done in Clapeyron.jl: https://github.com/ClapeyronThermo/Clapeyron.jl/blob/master/src/models/SAFT/SAFTVRMie/SAFTVRMie.jl#L167

[ianhbell]

I see now how the Aasen method in SAFT-VRQ works. It's very nice. Some python code that roughly implements the method:

# The classical method based on adaptive quadrature
def integrand(r_m):
    u = C*epskB*((sigma_m/r_m)**(lambda_r) - (sigma_m/r_m)**(lambda_a))
    return 1-np.exp(-u/T)
print('quasi-exact; (value, error estimate):')
print(scipy.integrate.quad(integrand, 0.0, sigma_m))

r0 = 0.8*sigma_m # manual rootfinding
print('integrand at r0:', integrand(r0))
for n in [5,10,20,40]:
    print(n, r0+scipy.integrate.fixed_quad(integrand, r0, sigma_m, n=n)[0])

yielding

quasi-exact; (value, error estimate):
(3.597838581533227e-10, 1.1907314145502913e-12)
integrand at r0: 0.9999999999999846
5 3.597658405106842e-10
10 3.597838619543654e-10
20 3.59783859272095e-10
40 3.5978385927209507e-10

They used 21-point quadrature which seems a little bit overkill, although very accurate. Somehow the Clapeyron.jl method is competitive without doing the interval splitting.