import os
import sys
import numpy as np
import scipy
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.ticker as mtick
import pandas as pd
from scipy import integrate
path='output'
exe_path=os.getcwd()
if os.path.isfile('path_file.dat'):
path_file=open('path_file.dat', 'r')
path=path_file.readline().rstrip()
path_file.close()
[docs]class latex_class():
"""
Setup for the use of LaTeX for axis labels and titles; sets of parameters
for graphics output.
"""
def __init__(self):
self.flag=False
self.dpi=300
self.font_size=14
self.tick_size=12
self.ext='jpg'
mpl.rc('text', usetex=False)
[docs] def on(self):
self.flag=True
mpl.rc('text', usetex=True)
[docs] def off(self):
self.flag=False
mpl.rc('text', usetex=False)
[docs] def set_param(self, dpi=300, fsize=14, tsize=12, ext='jpg'):
"""
Args:
dpi: resolution of the graphics file (default 300)
fsize: size of the labels of the axes in points (default 14)
tsize: size of the ticks in points (default 12)
ext: extension of the graphics file (default 'jpg'); this argument
is only used in those routines where the name of the file is
automatically produced by the program (e.g. check_poly or
check_spline functions). In other cases, the extension is
directly part of the name of the file given as argument to
the function itself, and 'ext' is ignored.
"""
self.dpi=dpi
self.font_size=fsize
self.tick_size=tsize
self.ext=ext
[docs] def get_dpi(self):
return self.dpi
[docs] def get_fontsize(self):
return self.font_size
[docs] def get_ext(self):
return self.ext
[docs] def get_tsize(self):
return self.tick_size
[docs]class name_data:
def __init__(self):
self.mineral_names=[]
[docs] def add(self,nlist):
self.mineral_names.extend(nlist)
[docs]class mineral:
def __init__(self,name,nick):
self.name=name
self.nick=nick
self.eos='m'
self.cp=[[0, 0]]
self.al=[[0, 0]]
self.k0=0.
self.kp=0.
self.dkt=0.
self.v0=0.
self.g0=0.
self.s0=0.
[docs] def info(self):
print("Mineral: %s\n" % self.name)
print("K0: %4.2f GPa, Kp: %4.2f, dK0/dT: %4.4f GPa/K, V0: %6.4f J/bar" \
% (self.k0, self.kp, self.dkt, self.v0))
print("G0: %8.2f J/mol, S0: %6.2f J/mol K\n" % (self.g0, self.s0))
print("EoS type: %s " % self.eos)
print("Cp coefficients and powers:")
for ci in self.cp:
print('{:>+8.4e}{:>+6.1f}'.format(ci[0],ci[1]))
print("\nAlpha coefficients and powers:")
for ai in self.al:
print('{:>+8.4e}{:>+6.1f}'.format(ai[0],ai[1]))
[docs] def load_ref(self,v0,g0,s0):
self.v0=v0
self.g0=g0
self.s0=s0
[docs] def load_bulk(self,k0, kp, dkt):
self.k0=k0
self.kp=kp
self.dkt=dkt
[docs] def load_eos_type(self,eos):
if eos==0:
self.eos='m'
else:
self.eos='bm'
[docs] def load_cp(self,cpc,cpp):
cl=list(zip(cpc,cpp))
item=0
self.cp=np.array([])
for ic in cl:
self.cp=np.append(self.cp,[cl[item][0], cl[item][1]])
item=item+1
self.cp=self.cp.reshape(item,2)
[docs] def load_alpha(self, alc, alp):
cl=list(zip(alc,alp))
item=0
self.al=np.array([])
for ic in cl:
self.al=np.append(self.al,[cl[item][0], cl[item][1]])
item=item+1
self.al=self.al.reshape(item,2)
[docs] def cp_t(self,tt):
cpt=0.
iterm=0
for it in self.cp:
cpt=cpt+self.cp[iterm,0]*(tt**self.cp[iterm,1])
iterm=iterm+1
return cpt
[docs] def alpha_t(self,tt):
alt=0.
iterm=0
for it in self.al:
alt=alt+self.al[iterm,0]*(tt**self.al[iterm,1])
iterm=iterm+1
return alt
[docs] def kt(self,tt):
return self.k0+(tt-298.15)*self.dkt
[docs] def entropy(self,tt):
fc=lambda ti: (self.cp_t(ti))/ti
integ, err=scipy.integrate.quad(fc,298.15,tt)
return integ+self.s0
[docs] def volume_t(self,tt):
fc= lambda ti: self.alpha_t(ti)
integ,err=scipy.integrate.quad(fc,298.15,tt)
# return (self.v0)*(1.+integ)
return (self.v0)*np.exp(integ)
[docs] def volume_p(self,tt,pp):
k0t=self.kt(tt)
vt=self.volume_t(tt)
if self.eos=='m':
fact=(1.-(pp*self.kp)/(pp*self.kp+k0t))**(1./self.kp)
return fact*vt
elif self.eos=='bm':
pf=lambda f: (3*k0t*f*(1+2*f)**(5/2)*(1+3*f*(self.kp-4)/3)-pp)**2
ff=scipy.optimize.minimize(pf,1,tol=0.00001)
return vt/((2*ff.x[0]+1)**(3/2))
[docs] def volume_fix_t_p(self,tt):
return lambda pp: self.volume_p(tt,pp)
[docs] def vdp(self,tt,pp):
fv=self.volume_fix_t_p(tt)
integ,err=scipy.integrate.quad(fv,0.000001, pp)
return integ*1e4
[docs] def g_t(self,tt):
integ,err=scipy.integrate.quad(self.entropy, 298.15, tt)
return integ
[docs] def g_tp(self,tt,pp):
return self.g0+self.vdp(tt,pp)-self.g_t(tt)
[docs] def alpha_p(self, tt, pp):
v=self.volume_p(tt,pp)
t_list=np.linspace(tt-10, tt+10, 5)
vt_list=np.array([])
for ti in t_list:
vi=self.volume_p(ti,pp)
vt_list=np.append(vt_list,vi)
fitpv=np.polyfit(t_list,vt_list,2)
fitder1=np.polyder(fitpv,1)
altp=np.polyval(fitder1,tt)
return 1*altp/v
[docs] def s_tp(self,tt,pp):
gtp=lambda tf: self.g_tp(tf,pp)
t_list=np.linspace(tt-5, tt+5, 5)
g_list=np.array([])
for ti in t_list:
gi=self.g_tp(ti,pp)
g_list=np.append(g_list,gi)
fit=np.polyfit(t_list,g_list,2)
fitder=np.polyder(fit,1)
return -1*np.polyval(fitder,tt)
[docs] def h_tp(self,tt,pp):
g=self.g_tp(tt,pp)
s=self.s_tp(tt,pp)
return g+tt*s
[docs] def pressure_vt(self, temp, volume, pmin=0., pmax=30., npp=40, p_del=1.):
p_list=np.linspace(pmin, pmax, npp)
v_list=np.array([self.volume_p(temp, pi) for pi in p_list])
v_diff=np.abs(v_list-volume)
ipos=np.argmin(v_diff)
p_approx=p_list[ipos]
p_min=p_approx-p_del/2
p_max=p_approx+p_del/2
if p_min <= 0.:
p_min=0.
p_list=np.linspace(p_min, p_max, np.int(npp/2))
v_list=np.array([self.volume_p(temp, pi) for pi in p_list])
v_diff=(v_list-volume)**2
fit_p=np.polyfit(p_list, v_diff, 2)
fit_der=np.polyder(fit_p, 1)
press=-1*fit_der[1]/fit_der[0]
return press
[docs] def helm_vt(self, temp, volume, pmin=0., pmax=30., npp=40, p_del=2):
pressure=self.pressure_vt(temp, volume, pmin, pmax, npp, p_del)
gibbs=self.g_tp(temp, pressure)
pv=1e4*pressure*volume
helm=gibbs-pv
return helm, gibbs, pressure
latex=latex_class()
name_list=name_data()
ens=mineral("enstatite","en")
cor=mineral("corindone","cor")
py=mineral("pyrope","py")
coe=mineral("coesite", "coe")
q=mineral("quartz","q")
fo=mineral("forsterite", "fo")
ky=mineral("kyanite","ky")
sill=mineral("sillimanite","sill")
andal=mineral("andalusite","and")
per=mineral("periclase","per")
sp=mineral("spinel","sp")
mao=mineral("maohokite","mao")
fmao=mineral("fe-maohokite","fmao")
stv=mineral("stishovite","stv")
cc=mineral("calcite", "cc")
arag=mineral("aragonite", "arag")
jeff=mineral("jeffbenite", "jeff")
jeff_fe=mineral("Fe-jeffbenite", "jeff_fe")
jeff_fe3p=mineral("Fe3p-jeffbenite", "jeff_fe3p")
jeff_feb=mineral("Feb-jeffbenite", "jeff_feb")
zrc=mineral("Zircon", "zrc")
diam=mineral("Diamond", "diam")
tui=mineral("Tuite", "tui")
flo=mineral("Florencite", "flo")
[docs]def export(mineral_nick):
al_power_dic={
'b1': 0,
'b2': 1,
'b3': -1,
'b4': -2,
'b5': -0.5
}
cp_power_dic={
'c1': 0,
'c2': 1,
'c3': -2,
'c4': 2,
'c5': -0.5,
'c6': -1,
'c7': -3,
'c8': 3
}
f=open('data_perplex.dat', 'a+')
nick=mineral_nick+'.'
alpha=eval(nick+'al')
cp=eval(nick+'cp')
g0=eval(nick+'g0')
v0=eval(nick+'v0')
s0=eval(nick+'s0')
name=eval(nick+'nick')
string_h=name+' '+'EoS = 2\n'
string_ini='G0 = '+str(g0)+' S0 = '+str(round(s0,4))+' V0 = '+str(round(v0,4))+'\n'
string_cp=''
for cpi in cp:
pi=cpi[1]
for ki,vi in cp_power_dic.items():
if vi==pi:
string_cp=string_cp+ki+' = '+str(round(cpi[0],6))+' '
string_al=''
for ali in alpha:
pi=ali[1]
for ki,vi in al_power_dic.items():
if vi==pi:
string_al=string_al+ki+' = '+str("%.6g" % ali[0])+' '
string_cp=string_cp+'\n'
string_al=string_al+'b6 = '+str(round(eval(nick+'k0')*1e4,6))+' '
string_al=string_al+'b7 = '+str(round(eval(nick+'dkt')*1e4,6))+' '
if eval(nick+'eos')=='m':
eos='b8 = '+str(round(eval(nick+'kp'),4))
elif eval(nick+'eos')=='bm':
eos='b8 = '+str(round(-1*eval(nick+'kp'),4))
string_al=string_al+eos+'\n\n'
f.write(string_h)
f.write(string_ini)
f.write(string_cp)
f.write(string_al)
f.close()
print("Exported data set in data_perplex.dat, with the following content:\n")
print(string_h+string_ini+string_cp+string_al)
[docs]def load_database():
ens.load_ref(6.262,-2915760, 132.5)
ens.load_bulk(107,4,-0.01605)
ens.load_cp([356.2, -.299E-2, -596900, -3185.3],[0, 1, -2, -1./2])
ens.load_alpha([.505E-4, -.505E-3],[0, -1./2])
ens.eos='m'
cor.load_ref(2.558, -1581710, 50.9)
cor.load_bulk(252,4,-0.0347)
cor.load_cp([139.5, .589E-2, -2460600,-589.2 ],\
[0, 1, -2, -1./2])
cor.load_alpha([.419E-4,-.419E-3],[0, -1./2])
cor.eos='m'
py.load_ref(11.318,-5934105,266.30)
py.load_bulk(173.7,4.00,-0.026055)
py.load_cp([633.5, -5196100, -4315.2],\
[0, -2, -0.5])
py.load_alpha([0.436E-4, -0.436E-3],\
[0, -0.5])
py.eos='bm'
load_database()
# ----------- Reactions ------------------
[docs]def equilib(tini,tfin,npoint,pini=1,prod=['py',1], rea=['ens',1.5,'cor', 1],\
out=False, tex=False, title=True, save=''):
"""
Computes the equilibrium pressure for a reaction involving a
given set of minerals, in a range of temperatures.
Args:
tini: minimum temperature in the range
tfin: maximum temperature in the range
npoint: number of points in the T range
pini (optional): initial guess for the pressure
prod: list of products of the reaction in the form
[name_1, c_name_1, name_2, c_name_2, ...]
where name_i is the name of the i^th mineral, as stored
in the database, and c_name_i is the corresponding
stoichiometric coefficient
rea: list of reactants; same syntax as the "prod" list.
Example:
equilib(300, 500, 12, prod=['py',1], rea=['ens', 1.5, 'cor', 1])
"""
if os.path.isfile("path_file.dat"):
path_file=open("path_file.dat", "r")
path=path_file.read()
path=path.rstrip()
lprod=len(prod)
lrea=len(rea)
prod_spec=prod[0:lprod:2]
prod_coef=prod[1:lprod:2]
rea_spec=rea[0:lrea:2]
rea_coef=rea[1:lrea:2]
lastr=rea_spec[-1]
lastp=prod_spec[-1]
prod_string=''
for pri in prod_spec:
prod_string=prod_string + pri
if pri != lastp:
prod_string=prod_string+' + '
rea_string=''
for ri in rea_spec:
rea_string = rea_string + ri
if ri != lastr:
rea_string=rea_string+' + '
t_list=np.linspace(tini,tfin,npoint)
p_list=np.array([])
h_list=np.array([])
s_list=np.array([])
v_list=np.array([])
cs_list=np.array([])
for ti in t_list:
pi=pressure_react(ti,pini, prod_spec, prod_coef, rea_spec, rea_coef)
hprod=0.
sprod=0.
vprod=0.
for pri, pci in zip(prod_spec, prod_coef):
hprod=hprod+(eval(pri+'.h_tp(ti,pi)'))*pci
sprod=sprod+(eval(pri+'.s_tp(ti,pi)'))*pci
vprod=vprod+(eval(pri+'.volume_p(ti,pi)'))*pci
hrea=0.
srea=0.
vrea=0.
for ri,rci in zip(rea_spec, rea_coef):
hrea=hrea+(eval(ri+'.h_tp(ti,pi)'))*rci
srea=srea+(eval(ri+'.s_tp(ti,pi)'))*rci
vrea=vrea+(eval(ri+'.volume_p(ti,pi)'))*rci
hi=hprod-hrea
si=sprod-srea
vi=vprod-vrea
dsdv_i=si/vi
p_list=np.append(p_list,pi)
h_list=np.append(h_list,hi)
s_list=np.append(s_list,si)
v_list=np.append(v_list,vi)
cs_list=np.append(cs_list, dsdv_i)
serie=(t_list.round(1),p_list.round(2),h_list.round(3), s_list.round(3), \
v_list.round(4), cs_list.round(2))
pd.set_option('colheader_justify', 'center')
df=pd.DataFrame(serie, index=['T (K)','P (GPa)','DH(J/mol)', \
'DS (J/mol K)', 'DV (J/bar)','Slope (bar/K)'])
df=df.T
df2=df.round(3)
print("")
print(df2.to_string(index=False))
ymax=max(p_list)+0.1*(max(p_list)-min(p_list))
ymin=min(p_list)-0.1*(max(p_list)-min(p_list))
xloc_py, yloc_py, xloc_en, yloc_en=field(tini,tfin, ymin, ymax, \
prod_spec, prod_coef, rea_spec, rea_coef)
dpi=80
if tex:
latex.on()
dpi=latex.get_dpi()
fontsize=latex.get_fontsize()
ticksize=latex.get_tsize()
print("\n")
fig=plt.figure()
ax=fig.add_subplot(111)
if title:
if latex.flag:
ax.title.set_text("Reaction "+ rea_string + ' $\leftrightarrow$ ' + prod_string + "\n" )
else:
ax.title.set_text("Reaction "+ rea_string + " <--> " + prod_string + "\n" )
ax.text(xloc_en, yloc_en, rea_string)
ax.text(xloc_py,yloc_py, prod_string)
ax.plot(t_list,p_list,"k-")
ax.axis([tini,tfin,ymin,ymax])
ax.yaxis.set_major_locator(plt.MaxNLocator(5))
ax.yaxis.set_major_formatter(mtick.FormatStrFormatter('%.2f'))
ax.xaxis.set_major_locator(plt.MaxNLocator(8))
ax.xaxis.set_major_formatter(mtick.FormatStrFormatter('%.1f'))
if latex.flag:
ax.set_xlabel("Temperature (K)", fontsize=fontsize)
ax.set_ylabel("Pressure (GPa)", fontsize=fontsize)
plt.xticks(fontsize=ticksize)
plt.yticks(fontsize=ticksize)
else:
ax.set_xlabel("Temperature (K)")
ax.set_ylabel("Pressure (GPa)")
if save != '':
plt.savefig(fname=path+'/'+ save,dpi=dpi, bbox_inches='tight')
plt.show()
latex.off()
clap=np.polyfit(t_list,p_list,1)
cl_s=clap[0]*1.e4
print("\nAverage Clapeyron Slope (from Delta S/Delta V): %6.2f bar/K" \
% cs_list.mean())
print("Clapeyron slope (from a linear fit of the P/T curve): %6.2f bar/K"\
% cl_s)
if out:
return t_list, p_list
[docs]def reaction(tt,pp, prod_spec, prod_coef, rea_spec, rea_coef):
"""
Computes the Gibbs free energy of reaction at given
temperature (tt) and pressure (pp), involving specified
minerals.
"""
gprod=0.
for pri, pci in zip(prod_spec, prod_coef):
gprod=gprod+(eval(pri+'.g_tp(tt,pp)'))*pci
grea=0.
for ri,rci in zip(rea_spec, rea_coef):
grea=grea+(eval(ri+'.g_tp(tt,pp)'))*rci
return gprod-grea
[docs]def pressure_react(tt,pini, prod_spec, prod_coef, rea_spec, rea_coef):
"""
Computes the pressure at which a given set of minerals
is at the equilibrium for a specified temperature.
"pini" is a initial guess for the pressure.
Output in GPa.
"pressure_react" calls "reactions, and it is invoked
by "equilib".
"""
fpr=lambda pp: (reaction(tt,pp, prod_spec, prod_coef, rea_spec, rea_coef))**2
pres=scipy.optimize.minimize(fpr,pini,tol=1)
return pres.x
[docs]def field(tmin,tmax,pmin,pmax,\
prod_spec, prod_coef, rea_spec, rea_coef, nx=6, ny=6):
t_range=np.linspace(tmin,tmax,nx)
p_range=np.linspace(pmin,pmax,ny)
fld=np.array([])
for ti in t_range:
for pi in p_range:
de=reaction(ti,pi,prod_spec, prod_coef, rea_spec, rea_coef)
fld=np.append(fld,[ti,pi,de])
fld=fld.reshape(nx*ny,3)
prodx=np.array([])
prody=np.array([])
reax=np.array([])
reay=np.array([])
for fi in fld:
if fi[2]>0:
reax=np.append(reax,fi[0])
reay=np.append(reay,fi[1])
else:
prodx=np.append(prodx,fi[0])
prody=np.append(prody,fi[1])
return prodx.mean(), prody.mean(), reax.mean(), reay.mean()
[docs]def import_database():
name_l=[]
al_power_dic={
'b1': 0,
'b2': 1,
'b3': -1,
'b4': -2,
'b5': -0.5,
}
cp_power_dic={
'c1': 0,
'c2': 1,
'c3': -2,
'c4': 2,
'c5': -0.5,
'c6': -1,
'c7': -3,
'c8': 3
}
list_cpc=[]
list_cp=[]
for ki,vi in cp_power_dic.items():
list_cpc.append(ki)
list_cp.append(vi)
list_cal=[]
list_al=[]
for ki,vi in al_power_dic.items():
list_cal.append(ki)
list_al.append(vi)
line=''
with open(exe_path+'/perplex_db.dat') as f:
jc=0
l0=['']
while True:
line=f.readline().rstrip()
if line=='': continue
if line == 'END': break
jc=jc+1
line_s=line.split()
if line_s != []:
l0=line_s[0].rstrip()
if l0=='#': continue
name=l0
name_l.append(l0)
l1=f.readline()
l2=f.readline().rstrip()
l2_s=l2.split()
g0=float(l2_s[2])
s0=float(l2_s[5])
v0=float(l2_s[8])
l3=f.readline().rstrip()
l3_s=l3.split()
l3n=len(l3_s)
coeff_cp=l3_s[2:l3n:3]
coeff_cp=[float(ci) for ci in coeff_cp]
power_ccp=l3_s[0:l3n:3]
power=[]
for cci in power_ccp:
power.append(cp_power_dic.get(cci))
l4=f.readline().rstrip()
l4_s=l4.split()
l4n=len(l4_s)
l4n_alpha=l4n-9
coeff_alpha=l4_s[2:l4n_alpha:3]
coeff_alpha=[float(ai) for ai in coeff_alpha]
power_ac=l4_s[0:l4n_alpha:3]
power_a=[]
for ai in power_ac:
power_a.append(al_power_dic.get(ai))
k0=float(l4_s[-7])/1.e4
dkt=float(l4_s[-4])/1.e4
kp=float(l4_s[-1])
eos_flag=0
eos='m'
if kp < 0.:
eos='bm'
kp=-1*kp
eos_flag=1
eval(name+'.load_ref(v0,g0,s0)')
eval(name+'.load_bulk(k0,kp,dkt)')
eval(name+'.load_cp(coeff_cp,power)')
eval(name+'.load_alpha(coeff_alpha,power_a)')
eval(name+'.load_eos_type(eos_flag)')
f.readline()
line=f.readline().rstrip()
f.close()
name_list.add(name_l)