corysimon / pyiast Goto Github PK

The fitting method for the dual-site Langmuir isotherm doesn't seem to work very well. I've included an example of data where it fails.

df = pd.read_csv('./data.csv', sep=' ')

T1 = 293.15
T2 = 313.15
T3 = 333.15
T4 = 363.15

df.columns = pd.MultiIndex.from_tuples(list(itertools.product([T1, T2, T3, T4], ['P', 'N'])))

m = pyiast.ModelIsotherm(df[T1].dropna(), loading_key='N', pressure_key='P', model='DSLangmuir')
m.print_params()

DSLangmuir identified model parameters:
	M1 = 0.271962
	K1 = -0.014806
	M2 = 0.189937
	K2 = -0.021286
RMSE =  2.75210883186

data.csv.tar.gz

I think this is just because the model is hard to fit to, even using a simple scipy.curve_fit seems to give "bad" answers even with a relatively accurate initial guess, eg: (here changing one parameter guess from 0.01 to 0.01 makes the optimisation snap to the correct values)

from scipy.optimize import curve_fit

def dualsite(P, q1, b1, q2, b2):
    return q1 * b1 * P / (1 + b1 * P) + q2 * b2 * P / (1 + b2 * P)

popt, pcov = curve_fit(dualsite, df[T1]['P'].values, df[T1]['N'].values, p0=[5, 0.01, 5, 0.01])

q1: 0.3308547822386095
b1: -3453538.574244023
q2: 7.386229859865762
b2: 0.003215815376835452

popt, pcov = curve_fit(dualsite, df[T1]['P'].values, df[T1]['N'].values, p0=[5, 0.001, 5, 0.01])
q1: 7.890734382768765
b1: 0.0016454550275279897
q2: 1.468241084209643
b2: 0.024404438085929743

I've had some luck in making lmfit fit stuff for me, it implements boundaries on parameters, which seem to work even with quite coarse estimates for initial values (see below again).

import lmfit

def dualsite_lmfit(params, x, data):
    q1 = params['q1']
    q2 = params['q2']
    b1 = params['b1']
    b2 = params['b2']
    
    model = dualsite(x, q1, b1, q2, b2)
    
    return model - data

params = lmfit.Parameters()

params.add('q1', value=5.0, min=0.0, max=20)
params.add('q2', value=5.0, min=0.0, max=20)
params.add('b1', value=1, min=0.0)
params.add('b2', value=1, min=0.0)

minner = lmfit.Minimizer(dualsite_lmfit, params, fcn_args=(x, y))

result = minner.minimize()

lmfit.report_fit(result)

[[Fit Statistics]]
    # function evals   = 152
    # data points      = 29
    # variables        = 4
    chi-square         = 0.028
    reduced chi-square = 0.001
    Akaike info crit   = -193.707
    Bayesian info crit = -188.238
[[Variables]]
    q1:   1.46819959 +/- 0.285037 (19.41%) (init= 5)
    q2:   7.89071518 +/- 0.250531 (3.18%) (init= 5)
    b1:   0.02440523 +/- 0.005533 (22.67%) (init= 1)
    b2:   0.00164549 +/- 0.000264 (16.03%) (init= 1)
[[Correlations]] (unreported correlations are <  0.100)
    C(q1, b1)                    = -0.980 
    C(q1, b2)                    = -0.957 
    C(b1, b2)                    =  0.893 
    C(q2, b2)                    = -0.759 
    C(q1, q2)                    =  0.542 
    C(q2, b1)                    = -0.408 

I could implement this to the package to improve LangmuirDS (and other models too?) if there's interest?

@CorySimon @SimonEnsemble what do you think about letting users create a ModelIsotherm directly from parameters instead of fitting them? This is useful when the raw isotherm data isn't presented in the literature, but authors report the parameters of the isotherms instead. With your permission, I'd like to make a PR to add this feature.

Is this part included in the program you give?

df_N2 = pd.read_csv("N2.csv")
N2_isotherm = pyiast.ModelIsotherm(df_N2,
loading_key="Loading(mmol/g)",
pressure_key="P(bar)",
model="Henry")

df_CO2 = pd.read_csv("../CO2.csv")
CO2_isotherm = pyiast.ModelIsotherm(df_CO2,re_key="P(bar)",
model="Langmuir")

df_H2O = pd.read_csv("../H2O.csv")
H2O_isotherm = pyiast.ModelIsotherm(df_H2O,
loading_key="Loading(mmol/g)",
pressure_key="P(bar)",
fill_vaule=df_H2O["Loading(mmol/g)"].max())
#list of CO2,N2 and H2O oartial pressures (bar)
partial_pressures = [0.166,0.679,0.02]
component_loadings = pyiast.iast(partial_pressures,
[CO2_isotherm,N2_isotherm, H2O_isotherm])

total_pressure = 65.0 #total bulk gas pressure (bar)
adsorbed_mole_fralctions = [0.5,0.5] #deslired adsorbed mole fractions
gas_mole_fractions, component_loadings = pyiast.reverse_iast(
adsorbed_mole_fralctions,
total_pressure,
[ch3ch3_isotherm, ch4_isotherm])

I find this in paper but don't find it in code

Hi!

Which units should I use to work with equilibrium data from liquid systems?

Using loadings in mmol/m3 and the bulk fluid phase concentration in mmol/m3 is enough to hold IAST for liquid systems?

Is this part included in the program you give?

df_N2 = pd.read_csv("N2.csv")
N2_isotherm = pyiast.ModelIsotherm(df_N2,
loading_key="Loading(mmol/g)",
pressure_key="P(bar)",
model="Henry")

df_CO2 = pd.read_csv("../CO2.csv")
CO2_isotherm = pyiast.ModelIsotherm(df_CO2,re_key="P(bar)",
model="Langmuir")

df_H2O = pd.read_csv("../H2O.csv")
H2O_isotherm = pyiast.ModelIsotherm(df_H2O,
loading_key="Loading(mmol/g)",
pressure_key="P(bar)",
fill_vaule=df_H2O["Loading(mmol/g)"].max())
#list of CO2,N2 and H2O oartial pressures (bar)
partial_pressures = [0.166,0.679,0.02]
component_loadings = pyiast.iast(partial_pressures,
[CO2_isotherm,N2_isotherm, H2O_isotherm])

total_pressure = 65.0 #total bulk gas pressure (bar)
adsorbed_mole_fralctions = [0.5,0.5] #deslired adsorbed mole fractions
gas_mole_fractions, component_loadings = pyiast.reverse_iast(
adsorbed_mole_fralctions,
total_pressure,
[ch3ch3_isotherm, ch4_isotherm])

I find this in paper but don't find it in code

Using the latest version (8ebc14) with pandas>=0.24.0 in python 2.7 the following FutureWarning shows up:

~/.virtualenvs/aiidapy/local/lib/python2.7/site-packages/pyiast/isotherms.py:434: FutureWarning: Sorting because non-concatenation axis is not aligned. A future version
of pandas will change to not sort by default.

To accept the future behavior, pass 'sort=False'.

To retain the current behavior and silence the warning, pass 'sort=True'.

  }, index=[0]), df

Can you please update to solve the issue?
Thanks!

way to do this is play naive:

adsorbed_mole_fraction_guess = np.zeros(2)
for g, gas in enumerate(gases):
    adsorbed_mole_fraction_guess[g] = ads_model[gas].loading(p[g])
adsorbed_mole_fraction_guess = adsorbed_mole_fraction_guess / np.sum(adsorbed_mole_fraction_guess)

naive b/c assumes adsorption of component g will be pure-component adsorption at partial pressure p[g] for this component (wrong, but good enough for an initial guess).

after experimenting this gives quite a robust initial guess for adsorbed_mole_fraction_guess.