SAFT gamma mie : How to add second order parameters

[YouHyin]

Hello Gustavo.
I'm trying to put second order parameters (2nd order NH – H2O) in the database in the following paper, do I just add values to the secondmie and secondasso tabs in saftgamma_database.xlsx?

Paper : Perdomo, F.A., et al., A predictive group-contribution framework for the thermodynamic modelling of CO2 absorption in cyclic amines, alkyl polyamines, alkanolamines and phase-change amines: New data and SAFT-γ Mie parameters. Fluid Phase Equilibria, 2023. 566: p. 113635.

from sgtpy import database
# setting an unlike mie interaction
group_k = 'NH'
group_l = 'H2O'
eps_kl = 488.28
lr_kl = 49.901
database.new_interaction_mie(group_k, group_l, eps_kl, lr_kl, overwrite=True)
# setting an unlike association interaction
group_k = 'NH'
group_l = 'H2O'
site_k = 'e1'
site_l = 'H'
epsAB_kl = 305.80
kAB_kl = 0.011099
database.new_interaction_asso(group_k, group_l, site_k, site_l, epsAB_kl, kAB_kl, overwrite=True)
group_k = 'NH'
group_l = 'H2O'
site_k = 'H'
site_l = 'e1'
epsAB_kl = 1756.2
kAB_kl = 24.260
database.new_interaction_asso(group_k, group_l, site_k, site_l, epsAB_kl, kAB_kl, overwrite=True)

Or I wrote this code based on this example(https://github.com/gustavochm/sgtpy/blob/master/examples/GC%20SAFT-Gamma-Mie%20Examples/1.%20Database.ipynb), but I'm not sure how to add the second order parameter because it only shows how to add the unlike Mie potential interaction, association parameters, not the second order parameter.

[YouHyin]

First, I added a parameter like this and ran the following code as you gave me like this.

import numpy as np
import matplotlib.pyplot as plt
from sgtpy import component, mixture, saftgammamie
from sgtpy.equilibrium import bubblePy


co2 = component(GC={'CO2': 1})

water = component(GC={'H2O': 1})

diethanolamine = component(GC={'NH': 1, 'CH2': 2, 'CH2OH': 2})

mixture = diethanolamine + co2 + water
mixture.saftgammamie()

eos = saftgammamie(mixture)


n = 2000
alpha = np.linspace(1e-5, 1, n)

mass_base = 1. 

mass_amine = 0.35 
moles_amine = mass_amine / eos.Mw[0] #

moles_co2 = alpha * moles_amine 
mass_co2 = moles_co2 * eos.Mw[1] 

mass_water = mass_base - mass_amine - mass_co2
moles_water = mass_water / eos.Mw[2]

moles_array = np.stack([moles_amine * np.ones(n), moles_co2, moles_water])

# mole fractions for alpha from 0 to 1
x = moles_array / np.sum(moles_array, axis=0)

T323 = 323.15 # K

y323 = np.zeros_like(x)
p323 = np.zeros(n)
vl323 =  np.zeros(n)
vv323 = np.zeros(n)

i = 0
y0 = np.array([0, 1, 0])
p0 = 1e5
sol = bubblePy(y0, p0, x[:, i], T348, eos, full_output=True)
y323[:, i] = sol.Y
p323[i] = sol.P
vl323[i] = sol.v1
vv323[i] = sol.v2


for i in range(1, n):
    y0 = y323[:, i-1]
    p0 =  p323[i-1]
    v0 = [vl323[i-1], vv323[i-1]]
    sol = bubblePy(y0, p0, x[:, i], T323, eos, full_output=True, v0=v0, good_initial=True)
    y323[:, i] = sol.Y
    p323[i] = sol.P
    vl323[i] = sol.v1
    vv323[i] = sol.v2
    
    
y_CO2_323 = y323[1, :]  
P_CO2_323 = y_CO2_323 * p323

but I get this error. How can I fix it?
C:\Users\user\anaconda3\lib\site-packages\sgtpy\gammamie_mixtures\saftgammamie_mixture.py:1238: RuntimeWarning: invalid value encountered in log
lnphi = mures - np.log(Z)

Thanks as always for your help!
I'm getting a lot of help.
Best whises!

[gustavochm]

Hi YouHyin,

This is something I didn't consider before. when I implemented this SAFT I only included the second-order parameters of the following paper:

Haslam, A. J., González-Pérez, A., Di Lecce, S., Khalit, S. H., Perdomo, F. A., Kournopoulos, S., Kohns, M., Lindeboom, T., Wehbe, M., Febra, S., Jackson, G., Adjiman, C. S., & Galindo, A. (2020). Expanding the Applications of the SAFT-γMie Group-Contribution Equation of State: Prediction of Thermodynamic Properties and Phase Behavior of Mixtures. Journal of Chemical and Engineering Data, 65(12), 5862–5890. https://doi.org/10.1021/acs.jced.0c00746

As you mentioned, the parameters are stored in the "secondmie" and "secondasso" of the saftgamma_database.xslx file. However, these parameters are only used in very specific molecular environments. In the end, I made custom functions to identify the molecular environment, but this is very manually done, see this file, for example.
If you modify the "secondmie" and "secondasso" database, the parameters will be there, but the specific molecular environment won't be checked, and hence, the custom parameters won't be used :/ I clearly will need to work on a more flexible API for this. Additionally, I'll definitely add Felipe's parameter in the next version of the package.

The other option that you show of replacing the "normal" parameters with the second-order parameters of the entire database should work only if the "normal" parameters were not used elsewhere in your mixture. I think this is the case with the mixture of diethanolamine + co2 + water. The warning you get is because probably the density solver probably failed at some point (hence a negative pressure and negative Z). This could be the case of starting from a very small value of CO2 (alpha=1e-5) and then even with a large value of n (n=2000 in this case), the previous solution as initial guess is not good enough (the volume/composition change is too big even with a small change of CO2).
If you just change alpha to alpha = np.linspace(1e-3, 1, n) does the trick without any warning :) (even works with n=100). I got the following plot

I hope that helps!
Gustavo

[YouHyin]

Hello Gustavo,
As you suggested, I wrote the following functions in this file https://github.com/gustavochm/sgtpy/blob/master/sgtpy/secondorder.py (and also added the parameters in the "secondmie" and "secondasso" of the saftgamma_database.xslx file).

if bool622 and bool62 and bool619:

            index22 = cond6_in22[0] + index0
            
            index14 = np.where(subgroups == 'H2O')[0]
            
            boolNH = secondorder.group_k == 'NH'
            boolH2O = secondorder.group_l == 'H2O'
            eps2214 = secondorder[boolNH & boolH2O]['eps_kl'].values
            lr2214 = secondorder[boolNH & boolH2O]['lr_kl'].values

            index17 = np.where(subgroups == 'CO2')[0]

            boolNH = secondorder.group_k == 'NH'
            boolCO2 = secondorder.group_l == 'CO2'
            eps2217 = secondorder[boolNH & boolCO2]['eps_kl'].values
            lr2217 = secondorder[boolNH & boolCO2]['lr_kl'].values

            eps_kl[index17, index22] = eps2217
            eps_kl[index22, index17] = eps2217

            lr_kl[index17, index22] = lr2217
            lr_kl[index22, index17] = lr2217

        del bool622, bool62, bool619

    bool22 = np.hstack(bool22)
    cond_asso22_0 = np.any(bool22)
    cond_asso22_14 = np.any(subgroup_id_asso == 'H2O')
    cond_asso22_17 = np.any(subgroup_id_asso == 'CO2')
    cond_asso22_19 = np.any(subgroup_id_asso == 'CH2OH')

    where22all = np.where(subgroup_id_asso == 'NH')
    where22all_id = molecule_id_index_asso[where22all]
    where22second = bool22[where22all_id]
    where22 = where22all[0][where22second]
    ncond22 = len(where22)

    if cond_asso22_0 and cond_asso22_14:

        where14 = np.where(subgroup_id_asso == 'H2O')

        boolk = secondasso.group_k == 'NH'
        booll = secondasso.group_l == 'H2O'
        df = secondasso[boolk & booll]
        len1 = df.shape[0]

        for k in range(ncond22):
            for j in range(len1):

                values = df.iloc[j].values
                groupK2, siteK, groupL2, siteL, epsAB, kAB = values

                if siteK == 'e1' and siteL == 'H':
                    indexH = sites_cumsum[where14]
                    indexNH = sites_cumsum[where22[k]]

                    epsAB_kl[indexH, indexNH] = epsAB
                    epsAB_kl[indexNH, indexH] = epsAB
                    kAB_kl[indexH, indexNH] = kAB
                    kAB_kl[indexNH, indexH] = kAB

And I solved the error (mixture of diethanolamine + co2 + water) as you advised.
Thank you so much for always being very detailed and kind to me.
Your advice is very helpful.
Thank you.
Best whises

[gustavochm]

Hi YouHyin,

Thank you for coding this. I'm afraid I haven't had the time to go through it yet.
Once the functions are available on that file they should be called in when reading the database (lines 180 and 809). I know this is a bit convoluted, but I didn't manage to come up with a more simple way to include these second-order modifications...

Potentially I will include all of these improvements in the next update.

Edit: I have already uploaded a new version of the database and updated the function for the second-order interactions. I still need to work on some other improvements before publishing the final update.

Thanks!
Gustavo

[gustavochm]

Hi YouHyin,

Just an update here.
I pushed an updated version of sgtpy into PyPI. The new version includes molecular parameter obtained in the following articles: Febra et al. (2021), Wehbe et al. (2022), Perdomo et al. (2023), Valsecchi et al. (2024) and Bernet et al. (2024).

I hope that helps in your research!

Regards,
Gustavo

[YouHyin]

Hi Gustavo,
Thank you so much. I think it will really help me a lot in my research, but I have a question here.
I tried the updated version and database, but without this code, especially in the case of 3-(methylamino)propylamine (MAPA) solution, it doesn't still fit the experimental data very well, do you know why?

from sgtpy import database
# setting an unlike mie interaction
group_k = 'NH'
group_l = 'H2O'
eps_kl = 488.28
lr_kl = 49.901
database.new_interaction_mie(group_k, group_l, eps_kl, lr_kl, overwrite=True)
# setting an unlike association interaction
group_k = 'NH'
group_l = 'H2O'
site_k = 'e1'
site_l = 'H'
epsAB_kl = 305.80
kAB_kl = 0.011099
database.new_interaction_asso(group_k, group_l, site_k, site_l, epsAB_kl, kAB_kl, overwrite=True)
group_k = 'NH'
group_l = 'H2O'
site_k = 'H'
site_l = 'e1'
epsAB_kl = 1756.2
kAB_kl = 24.260
database.new_interaction_asso(group_k, group_l, site_k, site_l, epsAB_kl, kAB_kl, overwrite=True)

Before entering the code

After entering the code

MAPA : CH3 :1, NH:1, CH2:3, NH2:1

Thank you
Best whises

[gustavochm]

Have you updated to the most recent version?

When updating the code, I realized there was a bug when adding new association interactions to the new database. Basically, there were some cases in which the reference name of the site ('e1', 'e2', or 'H') was incorrectly set. I corrected that bug in the updated version.

[YouHyin]

Yes I have updated to the most recent version.

1. After updating to the latest version, for diethanolamine (DEA) (SAFT-γ Mie group)

• NH : 1, CH2 : 2, CH2OH : 2) I confirmed that it matched the experimental data well without manually entering the above code.
  

2. However, in the case of 3-(methylamino)propylamine (MAPA), the updated version does not match the experimental data, and it was confirmed that the experimental data was well matched by entering the above code. (See the picture above)

Considering the cause, second-order parameters were introduced to model unlike interactions between NH2, NH, and N groups and H2O more precisely when the amine group is in close proximity to the hydroxyl group within the alkanolamine structure. Since DEA is an alkanolamine, it was well applied to the updated second-order parameters and fits well with the experimental data, but for MAPA, it is not an alkanolamine. Therefore, it does not apply to the updated second-order parameters. I think that exceptionally, when MAPA exists like H2O, there is a remarkable polarization effect in the -NH group, so I think it should replace the normal parameters of NH-H2O with the second-order parameters to match the experimental data.

Thank you always for your response.
Best whises

[gustavochm]

Hi YouHyin,

I'm glad to see that for DEA, it is working as expected now.

I understand your point about MAPA; the thing is that I coded the second-order interactions to be applied only in the described molecular environment, which MAPA doesn't meet for the NH/NH2 groups. You could, in principle, modify those interactions by manually changing the database (as you have done above).

As I mentioned, there was a bug when introducing the new association interactions, but now it should be fixed in the recent update.

Gustavo

[YouHyin]

Thank you so much for the update.
I'm really glad that the bug has been resolved!