Copyright 2018 Oliver Lord
This file is part of Xtools.
Xtools is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
Xtools is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License along with Xtools. If not, see http://www.gnu.org/licenses/.
In [8]:
%%capture
!jupyter nbextension enable --py --sys-prefix widgetsnbextension
In [9]:
import ipywidgets as widgets
import numpy as np
import sys
import array
import matplotlib.pyplot as plt
from ipywidgets import Layout, HBox, VBox, Dropdown, interactive, interact, Label
from IPython.core.display import display
from scipy.optimize import fsolve
from scipy.integrate import quad
In [10]:
global P_max, phase, EoS
P, P_max, V0, EoS, phase, ref, K, Kp, a0, b0, c0 = 0, 200, 162.373, 1, 1, 1, 256.7, 4.09, 4.7769, 4.9284, 6.8972
D0, g0, g_inf, q, apfu, apuc = 950, 1.54, 1, 1.5, 5, 4
VT = V0
def phases(phase,EoS,ref):
global V0, K, Kp, a0, b0, c0, D0, g0, g_inf, q, apfu, apuc
if phase == 'MgSiO3 perovskite':
#Data from: Tange et al. (2012) 10.1029/2011JB008988
dropdown3.options, dropdown3.disabled = ['Tange et al. (2012)'], True
dropdown2.options, dropdown2.disabled = ['3rd Order Birch-Murnaghan','Vinet'], False
if EoS == '3rd Order Birch-Murnaghan':
V0, K, Kp, a0, b0, c0 = 162.373, 256.7, 4.09, 4.7769, 4.9284, 6.8972
D0, g0, g_inf, q, apfu, apuc = 950, 1.54, 1, 1.5, 5, 4
if EoS =='Vinet':
V0, K, Kp, a0, b0, c0 = 162.373, 258.4, 4.1, 4.7769, 4.9284, 6.8972
D0, g0, g_inf, q, apfu, apuc = 940, 1.55, 1, 1.1, 5, 4
if phase == 'MgCO3 magnesite':
#Data from: Fiquet et al. (2002) American Mineralogist 87:1261-1265
dropdown2.value, dropdown2.disabled = '3rd Order Birch-Murnaghan', True
dropdown3.options, dropdown3.disabled = ['Fiquet et al. (2002)'], True
V0, K, Kp, a0, b0, c0 = 279.2, 108, 5, 4.6377, 4.6377, 15.007
D0, g0, g_inf, q, apfu, apuc = 1, 0, 0, 1, 1, 1
if phase == 'SiO2 stishovite':
#Data from: Wang et al. (2012) 10.1029/2011JB009100
dropdown3.options, dropdown3.disabled = ['Wang et al. (2012)'], True
dropdown2.options, dropdown2.disabled = ['3rd Order Birch-Murnaghan','Vinet'], False
if EoS == '3rd Order Birch-Murnaghan':
V0, K, Kp, a0, b0, c0 = 46.55, 294, 4.85, 4.66, 4.66, 2.8403
D0, g0, g_inf, q, apfu, apuc = 1, 0, 0, 1, 1, 1
if EoS == 'Vinet':
V0, K, Kp, a0, b0, c0 = 46.55, 292, 5.01, 4.66, 4.66, 2.8403
D0, g0, g_inf, q, apfu, apuc = 1, 0, 0, 1, 1, 1
if phase == 'Pt':
#Data from: Dorogokupets & Dewaele (2007) 10.1080/08957950701659700
dropdown2.value, dropdown2.disabled = '3rd Order Birch-Murnaghan', True
dropdown3.options, dropdown3.disabled = ['Dorogokupets & Dewaele (2007)'], True
V0, K, Kp, a0, b0, c0 = 60.38, 277.3, 5.12, 3.9231, 3.9231, 3.9231
D0, g0, g_inf, q, apfu, apuc = 1, 0, 0, 1, 1, 1
if phase == 'Au':
#Data from: Fei et al. (2012) 10.1073/pnas.0609013104
dropdown3.options, dropdown3.disabled = ['Fei et al. (2012)'], True
dropdown2.options, dropdown2.disabled = ['3rd Order Birch-Murnaghan','Vinet'], False
if EoS == '3rd Order Birch-Murnaghan':
V0, K, Kp, a0, b0, c0 = 67.85, 167, 5.77, 4.0786, 4.0786, 4.0786
D0, g0, g_inf, q, apfu, apuc = 1, 0, 0, 1, 1, 1
if EoS == 'Vinet':
V0, K, Kp, a0, b0, c0 = 67.85, 167, 6.00, 4.0786, 4.0786, 4.0786
D0, g0, g_inf, q, apfu, apuc = 1, 0, 0, 1, 1, 1
if phase == 'NaCl B1':
dropdown3.options, dropdown3.disabled = ['Dorogokupets & Dewaele (2007)', 'Decker (1971)'], False
dropdown2.options, dropdown2.disabled = ['3rd Order Birch-Murnaghan','Vinet'], True
#Data from: Dorogokupets & Dewaele (2007) 10.1080/08957950701659700
if ref == 'Dorogokupets & Dewaele (2007)':
dropdown2.value = 'Vinet'
V0, K, Kp, a0, b0, c0 = 179.44, 23.83, 5.09, 5.6404, 5.6404, 5.6404
D0, g0, g_inf, q, apfu, apuc = 1, 0, 0, 1, 1, 1
#Data from: Decker (1971) 10.1080/08957950701659700
if ref == 'Decker (1971)':
dropdown2.value = '3rd Order Birch-Murnaghan'
V0, K, Kp, a0, b0, c0 = 179.43, 24.02, 4.74, 5.6402, 5.6402, 5.6402
D0, g0, g_inf, q, apfu, apuc = 1, 0, 0, 1, 1, 1
if phase == 'NaCl B2':
#Data from: Dorogokupets & Dewaele (2007) 10.1080/08957950701659700
dropdown2.value, dropdown2.disabled = 'Vinet', True
dropdown3.options, dropdown3.disabled = ['Dorogokupets & Dewaele (2007)'], True
V0, K, Kp, a0, b0, c0 = 40.73, 29.72, 5.14, 3.4407, 3.4407, 3.4407
D0, g0, g_inf, q, apfu, apuc = 1, 0, 0, 1, 1, 1
if phase == 'CO2-V':
#Data from: Datchi et al. (2012) 10.1103/PhysRevLett.108.125701
dropdown2.value, dropdown2.disabled = 'Vinet', True
dropdown3.options, dropdown3.disabled = ['Datchi et al. (2012)'], True
V0, K, Kp, a0, b0, c0 = 90.99, 136.0, 4, 3.4407, 3.4407, 3.4407
#V0, K, Kp, a0, b0, c0 = 90.99, 136.0, 3.7, 3.7939, 3.7939, 6.3221
D0, g0, g_inf, q, apfu, apuc = 1, 0, 0, 1, 1, 1
initialise(EoS)
def phase_selector(phase):
dropdown2.value, dropdown2.disabled = 'MgSiO3 perovskite', False
dropdown3.value, dropdown3.disabled = '3rd Order Birch-Murnaghan', False
phases
def EoS_selector(EoS):
phases
def ref_selector(ref):
phases
def initialise(EoS):
K_box.value = ("K = " + str(K) + " GPa")
Kp_box.value = ("K' = " + str(Kp))
V0_box.value = ("V₀ = " + str(V0) + " ų")
def BM3(V,P):
return P-((3/2)*(K*(((V0/V)**(7/3))-((V0/V)**(5/3)))*(1-(((3/4)*(4-Kp))*(((V0/V)**(2/3))-1)))))
def Vinet(V,P):
return P-(3*K*((1-((V/V0)**(1/3)))/((V/V0)**(2/3)))*np.exp((3/2)*(Kp-1)*(1-((V/V0)**(1/3)))))
phase_items = ['MgSiO3 perovskite', 'MgCO3 magnesite', 'SiO2 stishovite', 'Pt', 'Au','NaCl B1','NaCl B2', 'CO2-V']
EoS_items = ['3rd Order Birch-Murnaghan', 'Vinet']
dropdown1 = widgets.Dropdown(options = phase_items,description='Phase:',layout=Layout(width='75.4%'))
dropdown2 = widgets.Dropdown(options = EoS_items,description='EoS:',layout=Layout(width='75.4%'),disabled=False)
dropdown3 = widgets.Dropdown(options = ['Tange et al. (2012)'],description='Ref:',layout=Layout(width='75.4%'),disabled=True)
w1=widgets.interactive(phase_selector,phase=dropdown1,EoS=dropdown2,ref=dropdown3)
w2=widgets.interactive(phases,phase=dropdown1,EoS=dropdown2,ref=dropdown3)
w3=widgets.interactive(phases,phase=dropdown1,EoS=dropdown2,ref=dropdown3)
K_box = widgets.Text(disabled=True,description='Parameters',layout=Layout(width='35%'))
Kp_box = widgets.Text(disabled=True,layout=Layout(width='19.8%'))
V0_box = widgets.Text(disabled=True,layout=Layout(width='19.8%'))
initialise('3rd Order Birch-Murnaghan')
VBox([HBox([dropdown1]), HBox([dropdown2]), HBox([dropdown3]), HBox([K_box, Kp_box, V0_box])])
In [11]:
global P_max, phase, EoS, V0, K, Kp, a0, b0, c0
def V_calc(P,T):
if g0 == 0:
T_slider.value, T_slider.disabled = 300, True
else:
T_slider.disabled = False
VVT_box.description = (str(round(T_slider.value)) + 'K:')
P_range, V_range, VT_range = list(range(0,P_max+1,4)), list(range(0,P_max+1,4)), list(range(0,P_max+1,4))
if EoS == 1:
Vp = fsolve(lambda V,P: BM3(V,P), 1, args=(P), xtol=1e-5)
VT = VT_calc(Vp,P,T)
for i, j in enumerate(P_range):
V_range[i] = fsolve(lambda V,P: BM3(V,P), 1, args=(j), xtol=1e-5)
if T != 300:
VT_range[i] = VT_calc(Vp,j,T)
else:
VT_range[i] = V_range[i]
if EoS == 2:
Vp = fsolve(lambda V,P: Vinet(V,P), 1, args=(P), xtol=1e-5)
VT = VT_calc(Vp,P,T)
for i, j in enumerate(P_range):
V_range[i] = fsolve(lambda V,P: Vinet(V,P), 1, args=(j), xtol=1e-5)
if T != 300:
VT_range[i] = VT_calc(Vp,j,T)
else:
VT_range[i] = V_range[i]
VV_box.value = ("V = " + str(np.around(Vp,3))[1:-1] + " ų")
Va_box.value = ("a = " + str(np.around((Vp/V0)**(1/3)*a0,3))[1:-1] + " Å")
Vb_box.value = ("b = " + str(np.around((Vp/V0)**(1/3)*b0,3))[1:-1] + " Å")
Vc_box.value = ("c = " + str(np.around((Vp/V0)**(1/3)*c0,3))[1:-1] + " Å")
VVT_box.value = ("V = " + str(np.around(VT,3))[1:-1] + " ų")
VaT_box.value = ("a = " + str(np.around((Vp/V0)**(1/3)*a0,3))[1:-1] + " Å")
VbT_box.value = ("b = " + str(np.around((Vp/V0)**(1/3)*b0,3))[1:-1] + " Å")
VcT_box.value = ("c = " + str(np.around((Vp/V0)**(1/3)*c0,3))[1:-1] + " Å")
PVT_plot(P_max,P_range,V_range,P,Vp,VT,VT_range)
def PVT_plot(P_max,P_range,V_range,P,Vp,VT,VT_range):
plt.figure(num=1, figsize=(10.72, 6), dpi=80, facecolor='w', edgecolor='k')
line1, line2, line3, line4 = plt.plot(P_range,VT_range,'r-',P_range,V_range,'b-',P,VT,'ro',P,Vp,'bo')
plt.xlim(-P_max*0.02, P_max*1.02)
plt.ylim(V_range[-1]*.98, VT_range[0]*1.02)
plt.ylabel('Unit Cell Volume (ų)', fontsize=18)
plt.xlabel('Pressure (GPa)', fontsize=18)
def VT_calc(V,P,T):
VT = fsolve(lambda V,P,T: MGD(V,P,T), 1, args=(P,T), xtol=1e-5)
return(VT)
def MGD(Vp,P,T):
if EoS == 1:
P300 = ((3/2)*(K*(((V0/Vp)**(7/3))-((V0/Vp)**(5/3)))*(1-(((3/4)*(4-Kp))*(((V0/Vp)**(2/3))-1)))))
if EoS == 2:
P300 = (3*K*((1-((Vp/V0)**(1/3)))/((Vp/V0)**(2/3)))*np.exp((3/2)*(Kp-1)*(1-((Vp/V0)**(1/3)))))
P300=((3/2)*(K*(((V0/Vp)**(7/3))-((V0/Vp)**(5/3)))*(1-(((3/4)*(4-Kp))*(((V0/Vp)**(2/3))-1)))))
g = g_inf+(g0-g_inf)*(Vp/V0)**q
D = D0*(Vp/V0)**-g
int_T = quad(debye_integrand, 0, D/T)
int_300 = quad(debye_integrand, 0, D/300)
E_T = 9*apfu*8.314472*T*((T/D)**3)*int_T[0]
E_300 = 9*apfu*8.314472*300*((300/D)**3)*int_300[0]
PTH = (E_T-E_300)*(g/(((Vp/apuc)*6.022E+23)/1E+21))
PTOTAL = P300 + PTH
return PTOTAL-P
def debye_integrand(x):
return x**3/((np.exp(x) - 1))
P_slider = widgets.FloatSlider(min=0, max=200, step=0.1, value=0, description='P in GPa:',layout=Layout(width='76.9%'),continuous_update=False)
P_slider.style.handle_color = 'lightblue'
T_slider = widgets.FloatSlider(min=300, max=3000, step=1, value=300, description='T in K:',layout=Layout(width='76.9%'),continuous_update=False)
T_slider.style.handle_color = 'lightblue'
widgets.interactive(V_calc,P=P_slider, T=T_slider)
VV_box = widgets.Text(disabled=True,layout=Layout(width='18.55%'),description=('300K:'))
Va_box = widgets.Text(disabled=True,layout=Layout(width='18.55%'))
Vb_box = widgets.Text(disabled=True,layout=Layout(width='18.55%'))
Vc_box = widgets.Text(disabled=True,layout=Layout(width='18.55%'))
VVT_box = widgets.Text(disabled=True,layout=Layout(width='18.55%'),description=(str(round(T_slider.value)) + 'K:'))
VaT_box = widgets.Text(disabled=True,layout=Layout(width='18.55%'))
VbT_box = widgets.Text(disabled=True,layout=Layout(width='18.55%'))
VcT_box = widgets.Text(disabled=True,layout=Layout(width='18.55%'))
intwidget = interactive(V_calc, P=P_slider, T=T_slider)
intoutput = intwidget.children[-1]
display(VBox([HBox([P_slider]), HBox([T_slider]), HBox([VV_box, Va_box, Vb_box, Vc_box]), HBox([VVT_box, VaT_box, VbT_box, VcT_box]), HBox([intoutput])]))
intwidget.update()
In [12]:
A, B = 1904, 9.500
def ruby_scale(scale1):
global A, B
if scale1 == 1:
#Data from: Dewaele et al. (2004) 10.1103/PhysRevB.70.094112
A, B = 1904, 9.500
if scale1 == 2:
#Data from: Mao et al. (1986)
A, B = 1904, 7.665
initialise2()
def raman_scale(scale2):
return scale2
def ruby(R1):
R1_0,V1_0=get_boxes()
delta=R1-R1_0
ruby_P=A/B*((1+delta/R1_0)**B-1)
P_slider2.value=ruby_P
raman_slider.value=(ruby_P*(1/0.501))+V1_0
def raman(V1):
R1_0,V1_0=get_boxes()
print(R1_0,V1_0)
raman_P=0.501*(V1-V1_0)
P_slider2.value=raman_P
ruby_slider.value=R1_0*(((B*((A/B)+raman_P))/A)**(1/B))
def cal_plot(i):
R1_0,V1_0=get_boxes()
dV = raman_slider.value-V1_0
dV_range = np.arange(0., (raman_slider.max-V1_0), 0.2)
plt.figure(num=1, figsize=(10.72, 6), dpi=80, facecolor='w', edgecolor='k')
line1, line2, = plt.plot(dV_range, dV_range*0.501, 'r-', dV, dV*0.501, 'bo')
plt.ylim(0, P_slider2.max)
plt.xlim(0, raman_slider.max-V1_0)
plt.ylabel('Ruby Pressure (GPa)', fontsize=18)
plt.xlabel('Δν (cm⁻¹)', fontsize=18)
def R1_0(R1_0_value):
R1_0 = (float(R1_0_box.value[0:-3]))
initialise2()
def V1_0(V1_0_value):
V1_0 = (float(V1_0_box.value[0:-5]))
initialise2()
def get_boxes():
R1_0 = (float(R1_0_box.value[0:-3]))
V1_0 = (float(V1_0_box.value[0:-5]))
return (R1_0,V1_0)
def initialise2():
R1_0,V1_0=get_boxes()
ruby_slider.max = 750
ruby_slider.value = R1_0
ruby_slider.min = R1_0
P_slider2.max = A/B*((1+(ruby_slider.max-R1_0)/R1_0)**B-1)
P_slider2.value = 0
raman_slider.max = ((P_slider2.max*(1/0.501))+V1_0)
raman_slider.value = V1_0
raman_slider.min = V1_0
ruby_dropdown = widgets.Dropdown(options = {'Hydrostatic - Dewaele et al. 2004': 1, 'Non-hydrostatic - Mao et al. 1986': 2},description='RUBY:',layout=Layout(width='35%'))
raman_dropdown = widgets.Dropdown(options = {'Walter et al. 2014': 1},description='RAMAN:',layout=Layout(width='35%'),disabled=True)
R1_0_box = widgets.Text(value='694.2 nm',disabled=False,layout=Layout(width='18%'),description='R1 λ₀:')
V1_0_box = widgets.Text(value='1332 cm⁻¹',disabled=False,layout=Layout(width='18%'),description='D ν₀:')
ruby_slider = widgets.FloatSlider(min=694.2, max=750, step=0.01, value=694.2, description='R1 λ:',layout=Layout(width='52.5%'),continuous_update=False)
ruby_slider.style.handle_color = 'lightblue'
raman_slider = widgets.FloatSlider(min=1332, max=1600, step=0.1, value=1332, description='001 ν:',layout=Layout(width='52.5%'),continuous_update=False)
raman_slider.style.handle_color = 'lightblue'
P_slider2 = widgets.FloatSlider(min=0, max=200, step=0.1, value=0, description='P in GPa:',layout=Layout(width='100%'),disabled=True,continuous_update=False)
P_slider2.style.handle_color = 'lightblue'
w5=widgets.interactive(ruby_scale,scale1=ruby_dropdown)
w6=widgets.interactive(raman_scale,scale2=raman_dropdown)
w7=widgets.interactive(R1_0,R1_0_value=R1_0_box)
w8=widgets.interactive(V1_0,V1_0_value=V1_0_box)
w9=widgets.interactive(ruby,R1=ruby_slider)
w10=widgets.interactive(raman,V1=raman_slider)
initialise2()
intwidget = interactive(cal_plot, i=P_slider2)
intoutput = intwidget.children[-1]
display(VBox([HBox([ruby_dropdown,R1_0_box,ruby_slider]),HBox([raman_dropdown,V1_0_box,raman_slider]),HBox([P_slider2]),HBox([intoutput])]))
intwidget.update()
In [14]:
h = 4.1357e-15
c = 299792458
def energy(e):
wavelength_slider.value=((h*c)/(e*1000))*1e10
def wavelength(w):
energy_slider.value=((h*c)/w)*1e7
energy_slider = widgets.FloatSlider(min=10, max=100, step=0.01, value=1, description='E (KeV)',layout=Layout(width='100%'),continuous_update=False,readout_format='.2f')
energy_slider.style.handle_color = 'lightblue'
wavelength_slider = widgets.FloatSlider(min=0.123985, max=1.23985, step=0.001, value=1.23985, description='λ (Å)',layout=Layout(width='100%'),continuous_update=False,readout_format='.4f')
wavelength_slider.style.handle_color = 'lightblue'
w11=widgets.interactive(energy,e=energy_slider)
w12=widgets.interactive(wavelength,w=wavelength_slider)
display(VBox([HBox([energy_slider]),HBox([wavelength_slider])]))
In [ ]: