I have constructed a function of dispersion for the contribution to the activity coefficient in the script of COSMO-PurePython.ipynb, the detail of code is below, I have written the code according to the section "3.3. Model of Hsieh et al. COSMO-SAC-dsp" in the article to add the dispersion term. But the predicted VLE plot for system 'CHLORODIFLUOROMETHANE+HEXAFLUOROPROPYLENE' is not consistent with experimental VLE data trend. I am wondering if I can add the dispersion term of activity coefficient into COSMO-PurePython.ipynb.
import pandas, numpy as np
import matplotlib.pyplot as plt
from collections import namedtuple
import json
import scipy.optimize
q0 = 79.53 # [A^2]
r0 = 66.69 # [A^3]
z_coordination = 10
c_hb = 85580.0 # kcal A^4 / mol/e^2
R = 8.3144598/4184 # 0.001987 # but really: 8.3144598/4184
sigma_hb = 0.0084
EPS = 3.667
AEFFPRIME = 7.5
EO = 2.395e-4
FPOL = (EPS-1.0)/(EPS+0.5)
ALPHA = (0.3*AEFFPRIME**(1.5))/(EO)
alpha_prime = FPOL*ALPHA
sigma_tabulated = np.linspace(-0.025, 0.025, 51)
sigma_m = np.tile(sigma_tabulated,(len(sigma_tabulated),1))
sigma_n = np.tile(np.array(sigma_tabulated,ndmin=2).T,(1,len(sigma_tabulated)))
sigma_acc = np.tril(sigma_n) + np.triu(sigma_m,1)
sigma_don = np.tril(sigma_m) + np.triu(sigma_n,1)
DELTAW = (alpha_prime/2)*(sigma_m + sigma_n)**2+ c_hb*np.maximum(0, sigma_acc - sigma_hb)*np.minimum(0, sigma_don+sigma_hb)
ΔW(σm,σn)=(α'/2)(σm+σn)^2+chb.max[0,σacc-σhb]min[0,σdon+σhb]
def get_sigma_profile(name):
df = pandas.read_csv('profiles/VT2005/Sigma_Profile_Database_Index_v2.txt', sep = '\t')
mask = df['Compound Name'] == name
assert(sum(mask)==1)
index = int(df[mask]['Index No.'].iloc[0])
V_COSMO = float(df[mask]['Vcosmo, A3'].iloc[0])
with open('profiles/VT2005/Sigma_Profiles_v2/VT2005-'+'{0:04d}'.format(index)+'-PROF.txt') as fp:
dd = pandas.read_csv(fp,names=['sigma [e/A^2]','p(sigma)*A [A^2]'],sep='\s+')
dd['A'] = dd['p(sigma)*A [A^2]'].sum()
dd['p(sigma)'] = dd['p(sigma)*A [A^2]']/dd['A']
return dd, V_COSMO
dispersive_parameter_lib = {
'C(sp3)': 115.7023,
'C(sp2)': 117.4650,
'C(sp)': 66.0691,
'-O-': 95.6184,
'=O': -11.0549,
'N(sp3)': 15.4901,
'N(sp2)': 84.6268,
'N(sp)': 109.6621,
'F': 52.9318,
'Cl': 104.2534,
'H(OH)': 19.3477,
'H(NH)': 141.1709,
'H(water)': 58.3301,
'H(other)': 0
}
n1=1 #C
n2=2 #F
n3=1 #H
n4=1 #Cl
n=n1+n2+n3+n4
disp_mole_1=1/n*(n1*dispersive_parameter_lib['C(sp3)']+n2*dispersive_parameter_lib['F']+n3*dispersive_parameter_lib['H(other)']+n4*dispersive_parameter_lib['Cl'])
N1=1 #C(sp3)
N2=2 #C(sp2)
N3=6 #F
N=N1+N2+N3
disp_mole_2=1/N*(N1*dispersive_parameter_lib['C(sp3)']+N2*dispersive_parameter_lib['N(sp2)']+N3*dispersive_parameter_lib['F'])
w=0.27027
A1=w*(0.5*(disp_mole_1+disp_mole_2)-(disp_mole_1*disp_mole_2)**0.5)**
def get_Gamma(T, psigma):
"""
Get the value of Γ (capital gamma) for the given sigma profile
"""
Gamma = np.ones_like(psigma)
AA = np.exp(-DELTAW/(R*T))*psigma
for i in range(50):
Gammanew = 1/np.sum(AA*Gamma, axis=1)
difference = np.abs((Gamma-Gammanew)/Gamma)
Gamma = (Gammanew + Gamma)/2
if np.max(difference) < 1e-8:
break
else:
pass
return Gamma
def get_lngamma_resid(T, i, psigma_mix, prof, lnGamma_mix = None):
"""
The residual contribution to ln(γ_i)
"""
# For the mixture
if lnGamma_mix is None:
lnGamma_mix = np.log(get_Gamma(T, np.array(psigma_mix)))
# For this component
psigma = np.array(prof['p(sigma)'])
A_i = prof['A'].iloc[0]
lnGammai = np.log(get_Gamma(T, psigma))
lngammai = A_i/AEFFPRIME*np.sum(psigma*(lnGamma_mix - lnGammai))
return lngammai
def get_lngamma_comb(T, x, i, profs, V_COSMO_A3):
"""
The combinatorial part of ln(γ_i)
"""
A = np.array([prof['p(sigma)*A [A^2]'].sum() for prof in profs])
q = A/q0
r = V_COSMO_A3/r0
theta_i = x[i]*q[i]/np.dot(x,q)
phi_i = x[i]*r[i]/np.dot(x,r)
l = z_coordination/2*(r-q) - (r-1)
return (np.log(phi_i/x[i])+z_coordination/2*q[i]*np.log(theta_i/phi_i)
+l[i]-phi_i/x[i]*np.dot(x,l))
def lngamma_dsp(A1,x,i):
Lngamma_Dsp=A1*x[i]
return Lngamma_Dsp
def get_lngamma(T, x, i, psigma_mix, profs, V_COSMO_A3, lnGamma_mix = None):
"""
Sum of the contributions to ln(γ_i)
"""
return (get_lngamma_resid(T, i, psigma_mix, profs[i], lnGamma_mix=lnGamma_mix)
+ get_lngamma_comb(T, x, i, profs, V_COSMO_A3) + lngamma_dsp(A1,x,i))
profs = [get_sigma_profile(n)[0] for n in ["CHLORODIFLUOROMETHANE",'HEXAFLUOROPROPYLENE']]
As = [prof['p(sigma)*A [A^2]'].sum() for prof in profs]
profs,V_COSMO_A3 = zip(*[get_sigma_profile(name) for name in ["CHLORODIFLUOROMETHANE",'HEXAFLUOROPROPYLENE']])
V_COSMO_A3 = np.array(V_COSMO_A3)
Gamma1,Gamma2=[],[]
T=293.15
x1=np.linspace(1e-6,1-1e-6,51)
for x0L in x1:
x=[x0L,1-x0L]
# 然后计算psigma_mix
psigma_mix = sum([x[i]*profs[i]['p(sigma)*A [A^2]'] for i in range(2)])/sum([x[i]*As[i] for i in range(2)])
lnGamma_mix = np.log(get_Gamma(T, np.array(psigma_mix)))
gamma1=np.exp(get_lngamma(T, x, 0, np.array(psigma_mix), profs, V_COSMO_A3, lnGamma_mix = lnGamma_mix))
gamma2=np.exp(get_lngamma(T, x, 1, np.array(psigma_mix), profs, V_COSMO_A3, lnGamma_mix = lnGamma_mix))
Gamma1.append(gamma1)
Gamma2.append(gamma2)
ps1=918.00
ps2=645.00
x2=1-x1
p=ps1*np.array(Gamma1)*x1+ps2*np.array(Gamma2)*x2
y1=ps1*np.array(Gamma1)*x1/p
import matplotlib.pyplot as plt
import pandas as pd
df=pd.read_excel(r'C:\Users\mi\OneDrive\桌面\撰写论文的数据\论文数据\Supporting Information\COSMO-SAC P-xy Data(已自动还原).xlsx',sheet_name='CHLORODIFLUOROMETHANE')
P1=df.iloc[3:18,2]
X1=df.iloc[3:18,3]
Y1=df.iloc[3:18,4]
P2=p
x1=x1
y1=y1
plt.figure(dpi=600)
plt.rc('font',family='Times New Roman')
plt.rcParams['xtick.direction'] = 'in'
plt.rcParams['ytick.direction'] = 'in'
plt.plot(X1,P1,'.',markersize=3,color='red',label='Liquid (Exp.)')
plt.plot(Y1,P1,'>',markersize=3,color='red',label='Vapor (Exp.)')
plt.plot(x1,P2,marker='.',markersize=3,color='black',label='Liquid (COSMO-SAC)')
plt.plot(y1,P2,marker='>',markersize=3,color='black',label='Vapor (COSMO-SAC)')
plt.xlabel('Liquid/Vapor Composition of CHLORODIFLUOROMETHANE')
plt.ylabel('Total Pressure(KPa)')
plt.title('CHLORODIFLUOROMETHANE+HEXAFLUOROPROPYLENE(293.15K)')
plt.legend()
plt.show()
I have constructed a function of dispersion for the contribution to the activity coefficient in the script of COSMO-PurePython.ipynb, the detail of code is below, I have written the code according to the section "3.3. Model of Hsieh et al. COSMO-SAC-dsp" in the article to add the dispersion term. But the predicted VLE plot for system 'CHLORODIFLUOROMETHANE+HEXAFLUOROPROPYLENE' is not consistent with experimental VLE data trend. I am wondering if I can add the dispersion term of activity coefficient into COSMO-PurePython.ipynb.