A question about defining a new component

[kaige-web]

Hi! I'm trying to do some calculations about drugs using clapeyron.jl. I referred to the information in the user guide (using your own parameters). I have written:

model = PCSAFT(["griseofulvin","water"];userlocations=["C:/Users/GK/Desktop/Clapeyron.jl-master/PCSAFT_like.csv","C:/Users/GK/Desktop/Clapeyron.jl-master/PCSAFT_unlike.csv","C:/Users/GK/Desktop/Clapeyron.jl-master/PCSAFT_assoc.csv"])
T = 310.15
p = 1.01325e5
z=[0.00024993751562109475, 0.9997500624843789]
γ1=activity_coefficient(model,p,T,z)

However, when I calculated the activity coefficient of griseofulvin in water, the results were wrong.
The activity coefficient of griseofulvin, 3.740122415836513e8, The activity coefficient of water 1.0000354641606668

In fact, the correct result should be 47384.9 for griseofulvin and 1.04428e+10 for water.

I feel like I have a problem defining new parameters. Could you provide an example of defining a new component?

Also, binary interactions are usually linear functions of temperature, which the current code doesn't seem to take into account.

[pw0908]

Hi there; I'm currently away from my laptop and either Andrés or Paul may be able to check the code for you. But, just looking at your set up, you seem to have done everything correctly.

Thinking about the thermodynamics though, I'm a bit confused about the magnitude of the experimental data points? Water is almost pure in this system, thus, the activity coefficient should approach 1. Similarly, the interactions between water and griseofulvin are unfavourable so the activity coefficient should be greater than 1.

Did you mean to calculate the fugacity coefficient? Or perhaps did the experimentalist use another definition for the activity coefficients?

With regards to temperature dependent k parameters, you are correct; this has been a very long-standing problem in Clapeyron and we're trying to address it now that we are trying to implement parameter estimation.

[kaige-web]

The activity coefficient of griseofulvin is 47384.9 obtained from the ratio of fugacity coefficients.

I got the activity coefficient of water wrong. It should be 1. But I confirmed the activity coefficient of griseofulvin.
This paper gives the activity coefficient of griseofulvin (4.99e4) at a certain solubility (7.54e-7).

https://pubs.acs.org/doi/10.1021/mp500824d

model = PCSAFT(["griseofulvin","water"];userlocations=["C:/Users/GK/Desktop/Clapeyron.jl-master/PCSAFT_like.csv","C:/Users/GK/Desktop/Clapeyron.jl-master/PCSAFT_unlike.csv","C:/Users/GK/Desktop/Clapeyron.jl-master/PCSAFT_assoc.csv"])
T = 310.15
p = 1.01325e5
z=[7.54e-7, 1-7.54e-7]
γ1=activity_coefficient(model,p,T,z)
print(γ1)
[4.981194806210843e8, 1.000000000340072]

4.981194806210843e8 4.99e4. It seems that the difference in quantity

I also tried the itraconazole in the above figure. The activity cofficient should be 1.13e4.

model = PCSAFT(["itraconazole","water"];userlocations=["C:/Users/GK/Desktop/Clapeyron.jl-master/PCSAFT_like.csv","C:/Users/GK/Desktop/Clapeyron.jl-master/PCSAFT_unlike.csv","C:/Users/GK/Desktop/Clapeyron.jl-master/PCSAFT_assoc.csv"])
T = 310.15
p = 1.01325e5
z=[1.28e-7, 1-1.28e-7]
γ1=activity_coefficient(model,p,T,z)
print(γ1)
[205095.67903461048, 1.0000000000075153]

[pw0908]

These results make more sense; one potential mistake is the water model used. We use an older one whilst this paper uses the latest one which has a long correlation for the segment size. Looking at the parameters, the fact that this water model is slightly "longer" might be able to account for the difference in activity coefficient.

Ive implemented a version of this water model in the past; @longemen3000 , should we implement this water model? It involved having a few if statements in our PC-SAFT code...

[longemen3000]

i recently added the capability of T-dependent kij for PCSAFT. We are discussing if we add all the changes to the original PCSAFT model or we make a new WaterPCSAFT (tentative name) with all the modifications

[longemen3000]

compiling the differences between Clapeyron "vanilla" PC-SAFT (ver. 0.3.3) and the PC-SAFT described in the paper

1. special correlation for water's σ https://doi.org/10.1016/j.cep.2007.02.034
2. kij with temperature dependence
3. combining rule for association strengh: Δ(i,j,a,b) = sqrt(Δ(i,i,a,a)*Δ(j,j,a,b)) ("Elliot Rule")
4. polar contributions (PCPSAFT) (in acetonitrile and dimethyl sulfoxide)

The first 3 options could be added to vanilla PCSAFT. (i'm working on a separate branch to add those,whereas 1 was already added on the master branch.)
For adding polar contributions, your best option would be to make a derivative PCSAFT model that includes said contribution. you can copy a (WIP) implementation from the PCPSAFT implementation in the ElectrolyteModels branch.

[longemen3000]

@kaige-web using the new pharmaPCSAFT:

griseofulvin_like = ParamTable(:like,(
  species = ["griseofulvin"],
  Mw = [352.77],
  m = [14.174],
  sigma = [3.372],
  epsilon = [221.261],
  n_H = [2],
  n_e = [2],
))

griseofulvin_unlike = ParamTable(:unlike,(
  species1 = ["griseofulvin"],
  species2 = ["water08"],
  k = [-0.155],
  kT = [2.641e-4],
))

griseofulvin_assoc  = ParamTable(:assoc,(
  species1 = ["griseofulvin"],
  site1 = ["H"],
  species2 = ["griseofulvin"],
  site2 = ["e"],
  epsilon_assoc = [1985.49],
  bondvol = [0.02],
))

model = pharmaPCSAFT(["griseofulvin","water"];userlocations=[
  griseofulvin_like,
  griseofulvin_unlike,
  griseofulvin_assoc
])
T = 310.15
p = 1.01325e5
z=[7.54e-7, 1-7.54e-7]

Results in:

julia> γ1=activity_coefficient(model,p,T,z)
2-element Vector{Float64}:
 52897.70400612142
     0.9994112472028018

[pw0908]

The differences in activity coefficients are now probably accounted for by the combining rule used by the association parameters. Unfortunately, this is not something we can support yet but you should be able to calculate manually and add to your csv files.

[kaige-web]

Wow, thank you so much for your hard work! Such swift action leads me to believe that this project must be comprehensive and wonderful in the future.

I used the python version of PC-SAFT (https://github.com/zmeri/PC-SAFT) and got the activity coefficient of griseofulvin (52857.1). Now, the results of the python version are almost indistinguishable from the results of your project (52897.7).

4.99e4 is the result of the Fortran version from professor Gabriele Sadowski. Differences may be due to different procedures, but overall results are acceptable

Thanks again for your help!

[pw0908]

The differences will most likely come down to the parameters used. The difference between zmeri's implementation probably has a different combining rule implemented for association which would explain the minor difference.

As for Sadowski's code, we haven't taken a look at it but it's possible they forgot to mention something about their implementation in the paper as the difference is significant enough that it couldn't really come down to the combining rule? Might be worth you contacting Sadowski's directly and asking for the code used to generate that table.

Out of curiosity, what are you planning to use Clapeyron for?

[kaige-web]

My advisor had a postdoc experience in Sadowski's group. So, we have her Fortran code. However, this Fortran codes perform poorly in calculating liquid-liquid equilibria for drug systems, such as LLE of drug and water to calculate solubility of amorphous drugs (https://pubs.acs.org/doi/10.1021/mp500824d) as well as LLE of drug and polymer to calculate phase separation of amorphous solid dispersions (https://www.sciencedirect.com/science/article/pii/S2590156720300323?via%3Dihub). So tp_flash also helped me a lot, avoiding my own coding.

The code of zmeri restricts the scheme of the association site (such as 2B, 3B, 4A, 4B 4C). Polymers may typically have hundreds or thousands of association sites. Here are the polymers I use.

<style> </style>

PEG 6000	6000	303.6	2.899	204.6	1799.8	0.02	2	2
PVP K30	58000	2360.6	2.71	205.599	0	0.02	521	521
HPMCE3	18000	944.2974	2.5049	142	1282.36	0.02	129	129
PVPVA 46	50000	1777.3	3.1222	205.03	0	0.02	528	528

I focus on thermodynamic modeling of solid dispersions, primarily modeling to determine reasonable drug loadings as well as stability under different relative humidity storage conditions.

I found a problem while using Clapeyron. The temperature entered must have a decimal point such as 310.15K, and an error will be reported if 310K is entered.

[pw0908]

In response to the issue you're having, try typing 310. instead of 310 (I.e. use float instead of integer values). That might help.

The work youre doing sounds really interesting! You'll be happy to know we're working on supporting SLE calculations in Clapeyron as well as more advanced flash algorithms which should be able to handle the LLE of systems you're planning to look at!

[kaige-web]

I am so looking forward to these! Thanks for your contributions!

[kaige-web]

Due to the current difficulties, could I ask what advanced flash algorithms are being developed? Although the LLE of drug and water was successful, the LLE of drug and polymer both failed. I tried genetic algorithm, dual annealing algorithm, and differential evolution algorithm, and the results were all disappointing.

[pw0908]

We are currently working on implementing the Michelsen algorithm and HELD algorithm. If you take a look at the papers related to the latter, you'll see that the benchmarks given involve examples of LLE in solvent+polymer systems.

For your case, if you just want to use the built-in methods, I would suggest using our RRTPFlash method where you specify the initial guess for the equilibrium constants and specify that you are looking for LLE. It's not perfect but it should converge to the correct equilibrium.

[kaige-web]

Thank you for your information! I will try these algorithms.

[kaige-web]

Hi, is RRTPFlash used like this?

KS=[50,128]
method = RRTPFlash(K0=KS)
(x_1PC, nvals_1PC, G_1PC) = tp_flash(model, p, T, n, method)

I used this but got a error.

MethodError: no method matching _fixpoint(::Clapeyron.var"#f#248"{pharmaPCSAFT{BasicIdeal}, Float64, Float64, Vector{Float64}, RRTPFlash{Vector{Int64}}, Vector{Int64}, Vector{Int64}, Vector{Int64}, Vector{Float64}, Vector{Float64}}, ::Vector{Int64}, ::Clapeyron.Solvers.SSFixPoint{Int64}, ::Float64, ::Float64, ::Int64, ::Bool)

[pw0908]

There are two problems:

1. You've defined the Ks as integers; Julia is very sensitive to types so you should always add a decimal place for your inputs
2. Since we're expecting a near-pure water phase, the values for Ks should be very different (orders of magnitude)
   Here's what I managed to get working:

KS=[0.001,128.]
method = RRTPFlash(K0=KS,spec=:lle)
(x_1PC, nvals_1PC, G_1PC) = tp_flash(model, p, T, n, method)

Output:

([0.5973268935514967 0.4026731065172074; 3.951321120179818e-12 0.9999999999960485], [7.539960486576025e-7 5.082877307088487e-7; 3.9513161324912614e-12 0.9999987377122692], -7.692474070520835)

Hopefully this matches your expectations?

[kaige-web]

Thanks for your answer. According to the K0 you provided, the system of drug and water I am calculating is indeed very effective.
But it doesn't seem to work with LLE for drugs and polymers. I tried to repeat this paper. https://www.sciencedirect.com/science/article/pii/S2590156720300323?via%3Dihub

At about 433.15K, the solubility of ITZ in HPC is 26wt% or 0.908 in mole fraction. I'm guessing that in the other phase, ITZ has a higher mole fraction (0.99?). This is because the molecular weights of ITZ （705.63） and HPC （20000） are so different!

For drugs and polymers, RRTPlfash results are unstable and vary with K0.

using Clapeyron
using BlackBoxOptim


ITZ_like = ParamTable(:like,(
  species = ["ITZ"],
  Mw = [705.63],
  m = [26.109],
  sigma = [2.166],
  epsilon = [252.346],
  n_H = [2],
  n_e = [2],
))

ITZ_unlike = ParamTable(:unlike,(
  species1 = ["ITZ"],
  species2 = ["HPC-UL"],
  k = [-0.037],
  kT = [0.],
))

ITZ_assoc  = ParamTable(:assoc,(
  species1 = ["ITZ"],
  site1 = ["H"],
  species2 = ["ITZ"],
  site2 = ["e"],
  epsilon_assoc = [1204.88],
  bondvol = [0.02],
))

model = pharmaPCSAFT(["ITZ","HPC-UL"];userlocations=[
  ITZ_like,
  ITZ_unlike,
  ITZ_assoc
])

T = 433.15
p = 1.01325e5


n10=0.99
n = [n10, 1-n10]


KS=[0.01,128.1]
method = RRTPFlash(K0=KS,spec=:lle,max_iters=10000)
(x_1PC, nvals_1PC, G_1PC) = tp_flash(model, p, T, n, method)

print(x_1PC)

I guess I can only look forward to your HELD algorithm. I read the original paper about HELD algorithm, but it's too hard for me who has never done algorithm training.

[pw0908]

Unfortunately, I have to agree that HELD will be much more well-suited for this problem... it will be quite a while until we fully release it Im sorry to say.

My only other piece of advice would be to plot the gibbs free energy of mixing for this system (using mixing(model,p,T,x,property)), see where the minima lie approximately and use that as the initial guess for K0. Given the size asymmetry here, I would expect values closer to KS=[1e-3,1e3].

[kaige-web]

Thanks for your advice. I am looking forward to the implementation of the HELD algorithm in the future.Thank you for your patience during this time. I believe that clapeyron.jl can become an important tool for thermodynamic research worldwide in the future！

[longemen3000]

Just an update here, with Clapeyron 0.3.7:

griseofulvin_like = ParamTable(:like,(
  species = ["griseofulvin"],
  Mw = [352.77],
  m = [14.174],
  sigma = [3.372],
  epsilon = [221.261],
  n_H = [2],
  n_e = [2],
))

griseofulvin_unlike = ParamTable(:unlike,(
  species1 = ["griseofulvin"],
  species2 = ["water08"],
  k = [-0.155],
  kT = [2.641e-4],
))

griseofulvin_assoc  = ParamTable(:assoc,(
  species1 = ["griseofulvin"],
  site1 = ["H"],
  species2 = ["griseofulvin"],
  site2 = ["e"],
  epsilon_assoc = [1985.49],
  bondvol = [0.02],
))

model = pharmaPCSAFT(["griseofulvin","water"];userlocations=[
  griseofulvin_like,
  griseofulvin_unlike,
  griseofulvin_assoc
])
T = 310.15
p = 1.01325e5
z=[7.54e-7, 1-7.54e-7]

this is new:

julia> Clapeyron.lle_init(model,p,T,z) #it takes some time
([[0.9999750008088736, 2.49991911265059e-5]], [-11.618639398827186])

where the first value is a good initial point for lle equilibria. it obtains it via michelsen tpd optimization. we also have MichelsenTPFlash that performs gibbs optimization after RR

[kaige-web]

Hi，Thanks for your many efforts in updating Clapeyron，including the new LLE algorithm （MichelsenTPFlash）. However, when updating to the latest version of v0.5.6, there seems to be a problem with the above code. When creating the model, an error was reported. It looks like there is an issue with parameter reading. Is there something new that causes this bug？
Missing values exist ∈ single parameter m: ["water08"].
in eval at base\boot.jl:368
in top-level scope at untitled:29
in at Clapeyron\3UbS1\src\utils\macros.jl:292
in #pharmaPCSAFT#579 at Clapeyron\3UbS1\src\utils\macros.jl:299
in build_eosmodel at Clapeyron\3UbS1\src\utils\macros.jl:463
in getparams at Clapeyron\3UbS1\src\database\database.jl:225
in is_valid_param at Clapeyron\3UbS1\src\database\database_rawparam.jl:308
in error at base\error.jl:44

`griseofulvin_like = ParamTable(:like,(
species = ["griseofulvin"],
Mw = [352.77],
m = [14.174],
sigma = [3.372],
epsilon = [221.261],
n_H = [2],
n_e = [2],
))

griseofulvin_unlike = ParamTable(:unlike,(
species1 = ["griseofulvin"],
species2 = ["water08"],
k = [-0.155],
kT = [2.641e-4],
))

griseofulvin_assoc = ParamTable(:assoc,(
species1 = ["griseofulvin"],
site1 = ["H"],
species2 = ["griseofulvin"],
site2 = ["e"],
epsilon_assoc = [1985.49],
bondvol = [0.02],
))

model = pharmaPCSAFT(["griseofulvin","water"];userlocations=[
griseofulvin_like,
griseofulvin_unlike,
griseofulvin_assoc
])
T = 310.15
p = 1.01325e5
z=[7.54e-7, 1-7.54e-7]`

[kaige-web]

In addition, in sle_slle example, I used the same code, but got the different results. That's weird.

[pw0908]

On the issue regarding the m parameter, there was an update a few versions back where we standardized the parameter names in Clapeyron. m has now been renamed segment:

griseofulvin_like = ParamTable(:like,(
       species = ["griseofulvin"],
       Mw = [352.77],
       segment = [14.174],
       sigma = [3.372],
       epsilon = [221.261],
       n_H = [2],
       n_e = [2],
       ))

model = pharmaPCSAFT(["griseofulvin","water"];userlocations=[
       griseofulvin_like,
       griseofulvin_unlike,
       griseofulvin_assoc
       ])
pharmaPCSAFT{BasicIdeal} with 2 components:
 "griseofulvin"
 "water"
Contains parameters: Mw, segment, sigma, epsilon, k, kT, epsilon_assoc, bondvol, water

For the SLE results: very weird, I've now re-run it on my laptop and it seems I had an older version of the paracetamol parameters where I accounted for the cross association between water and paracetamol. The values I get now are the same as you. I'll update the notebook.

[kaige-web]

Not only for the SLE results, but also for the activity of griseofulvin.
The previous version gives julia> γ1=activity_coefficient(model,p,T,z) 2-element Vector{Float64}: 52897.70400612142 0.9994112472028018 . The v0.5.6 gives γ1=activity_coefficient(model,p,T,z) 2-element Vector{Float64}: 1.4768315074772342e11 1.0000000003186411

I use the following code and only change m to segment. Also, suddenly discovered why the kij of griseofulvin and water 08 were defined, but the model was defined using griseofulvin and water?

(model = pharmaPCSAFT(["griseofulvin","water"];userlocations=[
 griseofulvin_like,
 griseofulvin_unlike,
 griseofulvin_assoc
]))

using Clapeyron

griseofulvin_like = ParamTable(:like,(
  species = ["griseofulvin"],
  Mw = [352.77],
  segment = [14.174],
  sigma = [3.372],
  epsilon = [221.261],
  n_H = [2],
  n_e = [2],
))

griseofulvin_unlike = ParamTable(:unlike,(
  species1 = ["griseofulvin"],
  species2 = ["water08"],
  k = [-0.155],
  kT = [2.641e-4],
))

griseofulvin_assoc  = ParamTable(:assoc,(
  species1 = ["griseofulvin"],
  site1 = ["H"],
  species2 = ["griseofulvin"],
  site2 = ["e"],
  epsilon_assoc = [1985.49],
  bondvol = [0.02],
))

model = pharmaPCSAFT(["griseofulvin","water"];userlocations=[
  griseofulvin_like,
  griseofulvin_unlike,
  griseofulvin_assoc
])
T = 310.15
p = 1.01325e5
z=[7.54e-7, 1-7.54e-7]

γ1=activity_coefficient(model,p,T,z)

[pw0908]

Do you remember which version you were running?

On the water vs water08 issue, that's just a typo; use water08 for everything.

[pw0908]

If you're comparing the values to the ones that Andrés shared ages ago, that might not be the best benchmark?

@longemen3000 can you reproduce those values for the activity coefficients

[longemen3000]

found the error in pharmaPCSAFT, it was introduced in 0.5.0. we were adding k twice into params.epsilon. i will add a test for this specific case and release a new version. With the fixing mixing rules:

julia> model = pharmaPCSAFT(["griseofulvin","water08"];userlocations=[
         griseofulvin_like,
           griseofulvin_unlike,
             griseofulvin_assoc
             ])
pharmaPCSAFT{BasicIdeal} with 2 components:
 "griseofulvin"
 "water08"
Contains parameters: Mw, segment, sigma, epsilon, k, kT, epsilon_assoc, bondvol, water      

julia> γ1=activity_coefficient(model,p,T,z)
2-element Vector{Float64}:
 55334.605821130834
     1.0000000001262213