Integration with Symbolics.jl

[longemen3000]

at the moment, due to how we define functions and the fact that our models need to support ForwardDiff.jl, we have the ability, in theory, to generate fully symbolic versions of an EoS. for a proof of concept, this can be done (if we remove the NaN-safe protections):

using Clapeyron, Symbolics
model = UNIFAC(["water","ethanol"])
@variables v0 T0, x1, x2
x = [x1,x2]
Ge = Clapeyron.excess_gibbs_free_energy(model,v0,T0,x)
dGe = Symbolics.gradient(Ge,x)
act = Symbolics.simplify(exp.(dGe) ./(Clapeyron.R̄ * T0))

that gives the expression for the activity coefficient in terms of T0 and x (v0 is ignored).

 (0.12027235504272604exp(8.31446261815324T0*((-2.4561x1*(2.4561 / (2.4561x1 + 2.6616x2) + (-6.032427210000001x1) / ((2.4561x1 + 2.6616x2)^2) + (-2.19256047x2*exp((801.9 + 0.007514(T0^2) - 3.824T0) / T0)) / ((2.4561x1 + 2.6616x2)^2) + (-1.73916441x2*exp((3.6156T0 - 1391.3 - 0.001144(T0^2)) / T0)) / ((2.4561x1 + 2.6616x2)^2) + (-2.60543088x2*exp((3.6156T0 - 1391.3 - 0.001144(T0^2)) / T0)) / ((2.4561x1 + 2.6616x2)^2))) / ((2.4561x1) / (2.4561x1 + 2.6616x2) + (0.7081x2*exp((3.6156T0 - 1391.3 - 0.001144(T0^2)) / T0)) / (2.4561x1 + 2.6616x2) + (1.0608x2*exp((3.6156T0 - 1391.3 - 0.001144(T0^2)) / T0)) / (2.4561x1 +
2.6616x2) + (0.8927x2*exp((801.9 + 0.007514(T0^2) - 3.824T0) / T0)) / (2.4561x1 + 2.6616x2)) - 2.4561log((2.4561x1) / (2.4561x1 + 2.6616x2) + (0.7081x2*exp((3.6156T0 - 1391.3 - 0.001144(T0^2)) / T0)) / (2.4561x1 + 2.6616x2) + (1.0608x2*exp((3.6156T0 - 1391.3 - 0.001144(T0^2)) / T0)) / (2.4561x1 + 2.6616x2)
+ (0.8927x2*exp((801.9 + 0.007514(T0^2) - 3.824T0) / T0)) / (2.4561x1 + 2.6616x2)) - x2*((1.7689((-1.73916441x2) / ((2.4561x1 + 2.6616x2)^2) + (2.4561exp((17.253 - 0.8389T0 - 0.0009021(T0^2)) / T0)) / (2.4561x1 + 2.6616x2) + (-2.60543088x2) / ((2.4561x1 + 2.6616x2)^2) + (-6.032427210000001x1*exp((17.253 - 0.8389T0 - 0.0009021(T0^2)) / T0)) / ((2.4561x1 + 2.6616x2)^2) + (-2.19256047x2*exp((4.746T0 - 1606.0 - 0.0009181(T0^2)) / T0)) / ((2.4561x1 + 2.6616x2)^2))) / ((0.7081x2) / (2.4561x1 + 2.6616x2) + (1.0608x2) / (2.4561x1 + 2.6616x2) + (0.8927x2*exp((4.746T0 - 1606.0 - 0.0009181(T0^2)) / T0)) / (2.4561x1 + 2.6616x2) + (2.4561x1*exp((17.253 - 0.8389T0 - 0.0009021(T0^2)) / T0)) / (2.4561x1 + 2.6616x2)) + (0.8927((-2.19256047x2) / ((2.4561x1 + 2.6616x2)^2) + (2.4561exp((8.673T0 - 1460.0 - 0.01641(T0^2)) / T0)) / (2.4561x1 + 2.6616x2) + (-6.032427210000001x1*exp((8.673T0 - 1460.0 - 0.01641(T0^2)) / T0)) / ((2.4561x1 + 2.6616x2)^2) + (-1.73916441x2*exp((4.674T0 - 2777.0 - 0.001551(T0^2)) / T0)) / ((2.4561x1 + 2.6616x2)^2) + (-2.60543088x2*exp((4.674T0 - 2777.0 - 0.001551(T0^2)) / T0)) / ((2.4561x1 + 2.6616x2)^2))) / ((0.8927x2) / (2.4561x1 + 2.6616x2) + (0.7081x2*exp((4.674T0 - 2777.0 - 0.001551(T0^2)) / T0)) / (2.4561x1 + 2.6616x2) + (1.0608x2*exp((4.674T0 - 2777.0 - 0.001551(T0^2)) / T0)) / (2.4561x1 + 2.6616x2) + (2.4561x1*exp((8.673T0 - 1460.0 - 0.01641(T0^2)) / T0)) / (2.4561x1 + 2.6616x2)))) + 8.31446261815324T0*((x1*(1.5106856162723488x1 + 1.9853131596298699x2)*(1.5106856162723488 / (1.5106856162723488x1 + 1.9853131596298699x2) + (-2.282171031212166(x1 + x2)) / ((1.5106856162723488x1 + 1.9853131596298699x2)^2))) / (1.5106856162723488(x1 + x2)) + (x2*(1.5106856162723488x1 + 1.9853131596298699x2)*(1.9853131596298699 / (1.5106856162723488x1 + 1.9853131596298699x2) + (-2.999184034049054(x1 + x2)) / ((1.5106856162723488x1 + 1.9853131596298699x2)^2))) / (1.9853131596298699(x1 + x2)) + (21.2870187x1*(2.4561x1 + 2.6616x2)*(4.25740374 / (1.7334(2.4561x1 + 2.6616x2)) + (-10.456609325814(1.7334x1 + 2.4951999999999996x2)) / (3.0046755600000004((2.4561x1 + 2.6616x2)^2)))) / (2.4561(1.7334x1 + 2.4951999999999996x2)) + (33.206121599999996x2*(2.4561x1 + 2.6616x2)*(4.6136174400000005 / (2.4951999999999996(2.4561x1 + 2.6616x2)) + (-16.311511052352(1.7334x1 + 2.4951999999999996x2)) / (6.2260230399999985((2.4561x1 + 2.6616x2)^2)))) / (2.6616(1.7334x1 + 2.4951999999999996x2)) + 12.2805log((2.4561(1.7334x1 + 2.4951999999999996x2)) / (1.7334(2.4561x1 + 2.6616x2))) + log((1.5106856162723488(x1 + x2)) / (1.5106856162723488x1 + 1.9853131596298699x2))))) / T0
 (0.12027235504272604exp(8.31446261815324T0*((-2.4561x1*((0.7081exp((3.6156T0 - 1391.3 - 0.001144(T0^2)) / T0)) / (2.4561x1 + 2.6616x2) + (1.0608exp((3.6156T0 - 1391.3 - 0.001144(T0^2)) / T0)) / (2.4561x1 + 2.6616x2) + (0.8927exp((801.9 + 0.007514(T0^2) - 3.824T0) / T0)) / (2.4561x1 + 2.6616x2) + (-6.53715576x1) / ((2.4561x1 + 2.6616x2)^2) + (-2.3760103200000002x2*exp((801.9 + 0.007514(T0^2) - 3.824T0) / T0)) / ((2.4561x1 + 2.6616x2)^2) + (-1.8846789599999998x2*exp((3.6156T0 - 1391.3 - 0.001144(T0^2)) / T0)) / ((2.4561x1 + 2.6616x2)^2) + (-2.82342528x2*exp((3.6156T0 - 1391.3 - 0.001144(T0^2)) / T0)) / ((2.4561x1 + 2.6616x2)^2))) / ((2.4561x1) / (2.4561x1 + 2.6616x2) + (0.7081x2*exp((3.6156T0 - 1391.3 - 0.001144(T0^2)) / T0)) / (2.4561x1 + 2.6616x2) + (1.0608x2*exp((3.6156T0 - 1391.3 - 0.001144(T0^2)) / T0)) / (2.4561x1 + 2.6616x2) + (0.8927x2*exp((801.9 + 0.007514(T0^2) - 3.824T0) / T0)) / (2.4561x1 + 2.6616x2)) - 0.8927(log((0.8927x2) / (2.4561x1 + 2.6616x2) + (0.7081x2*exp((4.674T0 - 2777.0 - 0.001551(T0^2)) / T0)) / (2.4561x1 + 2.6616x2) + (1.0608x2*exp((4.674T0 - 2777.0 - 0.001551(T0^2)) / T0)) / (2.4561x1 + 2.6616x2) + (2.4561x1*exp((8.673T0 - 1460.0 - 0.01641(T0^2)) / T0)) / (2.4561x1 + 2.6616x2)) - log(0.37571385632702137(0.8927 + 1.7689exp((4.674T0 - 2777.0 - 0.001551(T0^2)) / T0)))) - 1.7689(log((0.7081x2) / (2.4561x1 + 2.6616x2) + (1.0608x2) / (2.4561x1
+ 2.6616x2) + (0.8927x2*exp((4.746T0 - 1606.0 - 0.0009181(T0^2)) / T0)) / (2.4561x1 + 2.6616x2) + (2.4561x1*exp((17.253 - 0.8389T0 - 0.0009021(T0^2)) / T0)) / (2.4561x1 + 2.6616x2)) - log(0.37571385632702137(1.7689 + 0.8927exp((4.746T0 - 1606.0 - 0.0009181(T0^2)) / T0)))) - x2*((1.7689(0.7081 / (2.4561x1 +
2.6616x2) + 1.0608 / (2.4561x1 + 2.6616x2) + (-1.8846789599999998x2) / ((2.4561x1 + 2.6616x2)^2) + (0.8927exp((4.746T0 - 1606.0 - 0.0009181(T0^2)) / T0))
/ (2.4561x1 + 2.6616x2) + (-2.82342528x2) / ((2.4561x1 + 2.6616x2)^2) + (-6.53715576x1*exp((17.253 - 0.8389T0 - 0.0009021(T0^2)) / T0)) / ((2.4561x1 + 2.6616x2)^2) + (-2.3760103200000002x2*exp((4.746T0 - 1606.0 - 0.0009181(T0^2)) / T0)) / ((2.4561x1 + 2.6616x2)^2))) / ((0.7081x2) / (2.4561x1 + 2.6616x2) + (1.0608x2) / (2.4561x1 + 2.6616x2) + (0.8927x2*exp((4.746T0 - 1606.0 - 0.0009181(T0^2)) / T0)) / (2.4561x1 + 2.6616x2) + (2.4561x1*exp((17.253 - 0.8389T0 - 0.0009021(T0^2)) / T0)) / (2.4561x1 + 2.6616x2)) + (0.8927(0.8927 / (2.4561x1 + 2.6616x2) + (-2.3760103200000002x2) / ((2.4561x1 + 2.6616x2)^2) + (0.7081exp((4.674T0 - 2777.0 - 0.001551(T0^2)) / T0)) / (2.4561x1 + 2.6616x2) + (1.0608exp((4.674T0 - 2777.0 - 0.001551(T0^2)) / T0)) / (2.4561x1 + 2.6616x2) + (-6.53715576x1*exp((8.673T0 - 1460.0 - 0.01641(T0^2)) / T0)) / ((2.4561x1 + 2.6616x2)^2) + (-1.8846789599999998x2*exp((4.674T0 - 2777.0 - 0.001551(T0^2)) / T0)) / ((2.4561x1 + 2.6616x2)^2) + (-2.82342528x2*exp((4.674T0 - 2777.0 - 0.001551(T0^2)) / T0)) / ((2.4561x1 + 2.6616x2)^2))) / ((0.8927x2) / (2.4561x1
+ 2.6616x2) + (0.7081x2*exp((4.674T0 - 2777.0 - 0.001551(T0^2)) / T0)) / (2.4561x1 + 2.6616x2) + (1.0608x2*exp((4.674T0 - 2777.0 - 0.001551(T0^2)) / T0))
/ (2.4561x1 + 2.6616x2) + (2.4561x1*exp((8.673T0 - 1460.0 - 0.01641(T0^2)) / T0)) / (2.4561x1 + 2.6616x2)))) + 8.31446261815324T0*((x1*(1.5106856162723488x1 + 1.9853131596298699x2)*(1.5106856162723488 / (1.5106856162723488x1 + 1.9853131596298699x2) + (-2.999184034049054(x1 + x2)) / ((1.5106856162723488x1 +
1.9853131596298699x2)^2))) / (1.5106856162723488(x1 + x2)) + (x2*(1.5106856162723488x1 + 1.9853131596298699x2)*(1.9853131596298699 / (1.5106856162723488x1 + 1.9853131596298699x2) + (-3.941468341799537(x1 + x2)) / ((1.5106856162723488x1 + 1.9853131596298699x2)^2))) / (1.9853131596298699(x1 + x2)) + (21.2870187x1*(2.4561x1 + 2.6616x2)*(6.12846072 / (1.7334(2.4561x1 + 2.6616x2)) + (-11.331505794384002(1.7334x1 + 2.4951999999999996x2)) / (3.0046755600000004((2.4561x1 + 2.6616x2)^2)))) / (2.4561(1.7334x1 + 2.4951999999999996x2)) + (33.206121599999996x2*(2.4561x1 + 2.6616x2)*(6.641224319999999 / (2.4951999999999996(2.4561x1 + 2.6616x2)) + (-17.676282650111997(1.7334x1 + 2.4951999999999996x2)) / (6.2260230399999985((2.4561x1 + 2.6616x2)^2)))) / (2.6616(1.7334x1 + 2.4951999999999996x2)) + 13.308log((2.6616(1.7334x1 + 2.4951999999999996x2)) / (2.4951999999999996(2.4561x1 + 2.6616x2))) + log((1.9853131596298699(x1 + x2)) / (1.5106856162723488x1 + 1.9853131596298699x2))))) / T0

In practice, we need additional support to make that a seamless experience. my vision is that:
Clapeyron.activity_coefficient(model,v,T,z) just works and returns the expression above. for that, we would need:

• Register nan-safe functions (and its derivatives) (tier 0: support primal functions)
• Register the differentials so they return symbolic derivatives instead of ForwardDiff ones (tier 1:support bulk properties)

i don't know if we can/should support higher tiers. any commentary about that it is appreciated