notebook.community



In [1]:

    
%cat 0Source_Citation.txt









    



Source and citation

- This notebook is a part of the `pytheos` package.
- Website: http://github.com/SHDShim/pytheos.
- How to cite: S.-H. Shim (2017) Pytheos - a python tool set for equations of state. DOI: 10.5281/zenodo.802392



In [2]:

    
%matplotlib inline
# %matplotlib notebook # for interactive

For high dpi displays.



In [3]:

    
%config InlineBackend.figure_format = 'retina'

0. General note

This example compares pressure calculated from pytheos and original publication for the gold scale by Tsuchiya 2003.

1. Global setup



In [4]:

    
import matplotlib.pyplot as plt
import numpy as np
from uncertainties import unumpy as unp
import pytheos as eos

3. Compare



In [5]:

    
eta = np.linspace(0., 0.34, 18)
print(eta)









    



[0.   0.02 0.04 0.06 0.08 0.1  0.12 0.14 0.16 0.18 0.2  0.22 0.24 0.26
 0.28 0.3  0.32 0.34]



In [6]:

    
tsuchiya_au = eos.gold.Tsuchiya2003()



In [7]:

    
help(tsuchiya_au)









    



Help on Tsuchiya2003 in module pytheos.scales.gold object:

class Tsuchiya2003(pytheos.scales.objs.MGEOS)
 |  Tsuchiya 2003 JGR 108. 2462+
 |  
 |  Original k0 = 166.7, k0p = 6.12.  However, I cannot reproduce their table
 |  exactly unless I change the values to 166.1.
 |  
 |  If reproduce_table=True, the adjusted k0 will be used for reproducing
 |  the table value down to the first number after decimal point.
 |  
 |  Method resolution order:
 |      Tsuchiya2003
 |      pytheos.scales.objs.MGEOS
 |      builtins.object
 |  
 |  Methods defined here:
 |  
 |  __init__(self, v0=67.84742110765599, reproduce_table=False)
 |      :param params_st: elastic parameters for static EOS in an OrderedDict
 |          [v0 in A^3, k0 in GPa, k0p]
 |      :param params_th: thermal parameters for thermal EOS in an OrderedDict.
 |          The component can differ depending on the equation used.
 |      :param params_anh: anharmonic parameters for anharmonic correction in
 |          an OrderedDict.  The component can differ depending on the
 |          equation used.
 |      :param params_el: electronic parameters for electronic correction in
 |          an OrderedDict. The component can differ depending on the
 |          equation used.
 |      :param eqn_st: equation type for the static EOS. 'bm3', 'vinet', or
 |          'kunc'
 |      :param eqn_th: equation type for the thermal EOS. 'constq', 'tange',
 |          'speziale', 'dorogokupets2007', 'dorogokupets2015', 'alphakt'
 |      :param eqn_anh: equation type for anharmonic correction. 'zharkov',
 |      :param eqn_el: equation type for electonic correction. 'zharkov',
 |          'tsuchiya'
 |      :param t_ref: reference temperature, 300 K
 |      :param three_r: 3 times gas constant.
 |          Jamieson modified this value to compensate for mismatches
 |      :param reference: reference for the EOS
 |  
 |  ----------------------------------------------------------------------
 |  Methods inherited from pytheos.scales.objs.MGEOS:
 |  
 |  cal_p(self, v, temp)
 |      calculate total pressure at given volume and temperature
 |      
 |      :param v: unit-cell volume in A^3
 |      :param temp: temperature in K
 |      :return: pressure in GPa
 |      :note: 2017/05/10 temp must be numpy array.  If not, such as list,
 |          create an error.
 |  
 |  cal_panh(self, v, temp)
 |      calculate pressure from anharmonic contributions
 |      
 |      :param v: unit-cell volume in A^3
 |      :param temp: temperature in K
 |      :return: pressure in GPa
 |  
 |  cal_pel(self, v, temp)
 |      calculate pressure from electronic contributions
 |      
 |      :param v: unit-cell volume in A^3
 |      :param temp: temperature in K
 |      :return: pressure in GPa
 |  
 |  cal_pst(self, v)
 |      calculate static pressure at 300 K.
 |      
 |      :param v: unit-cell volume in A^3
 |      :return: static pressure at t_ref (=300 K) in GPa
 |  
 |  cal_pth(self, v, temp)
 |      calculate thermal pressure
 |      
 |      :param v: unit-cell volume in A^3
 |      :param temp: temperature in K
 |      :return: thermal pressure in GPa
 |  
 |  cal_v(self, p, temp, min_strain=0.2, max_strain=1.0)
 |      calculate unit-cell volume at given pressure and temperature
 |      
 |      :param p: pressure in GPa
 |      :param temp: temperature in K
 |      :param min_strain: minimum strain searched for volume root
 |      :param max_strain: maximum strain searched for volume root
 |      :return: unit-cell volume in A^3
 |      :note: 2017/05/10 I found wrap function is not compatible with
 |          OrderedDict. So I convert unp array to np array.
 |  
 |  print_equations(self)
 |      show equations used for the EOS
 |  
 |  print_parameters(self)
 |      show thermoelastic parameters for the EOS
 |  
 |  print_reference(self)
 |      show reference for the EOS
 |  
 |  ----------------------------------------------------------------------
 |  Data descriptors inherited from pytheos.scales.objs.MGEOS:
 |  
 |  __dict__
 |      dictionary for instance variables (if defined)
 |  
 |  __weakref__
 |      list of weak references to the object (if defined)



In [8]:

    
tsuchiya_au.print_equations()









    



P_static:  vinet
P_thermal:  constq
P_anharmonic:  None
P_electronic:  None



In [9]:

    
tsuchiya_au.print_equations()









    



P_static:  vinet
P_thermal:  constq
P_anharmonic:  None
P_electronic:  None



In [10]:

    
tsuchiya_au.print_parameters()









    



Static:  OrderedDict([('v0', 67.84742110765599+/-0.001), ('k0', 166.7+/-0), ('k0p', 6.12+/-0)])
Thermal:  OrderedDict([('v0', 67.84742110765599+/-0.001), ('gamma0', 3.16+/-0), ('q', 2.15+/-0), ('theta0', 180.0+/-0)])
Anharmonic:  None
Electronic:  None



In [11]:

    
v0 = 67.84742110765599



In [12]:

    
tsuchiya_au.three_r









    Out[12]:





24.943379399999998



In [13]:

    
v = v0 * (1.-eta) 
temp = 2500.



In [14]:

    
p = tsuchiya_au.cal_p(v, temp * np.ones_like(v))



In [15]:

    
print('for T = ', temp)
for eta_i, p_i in zip(eta, p):
    print("{0: .3f} {1: .2f}".format(eta_i, p_i))









    



for T =  2500.0
 0.000  16.94+/-0.00
 0.020  20.13+/-0.00
 0.040  23.86+/-0.00
 0.060  28.21+/-0.00
 0.080  33.26+/-0.00
 0.100  39.11+/-0.00
 0.120  45.89+/-0.00
 0.140  53.72+/-0.01
 0.160  62.75+/-0.01
 0.180  73.17+/-0.01
 0.200  85.19+/-0.01
 0.220  99.04+/-0.01
 0.240  115.01+/-0.01
 0.260  133.45+/-0.01
 0.280  154.73+/-0.01
 0.300  179.34+/-0.01
 0.320  207.81+/-0.02
 0.340  240.81+/-0.02



In [16]:

    
v = tsuchiya_au.cal_v(p, temp * np.ones_like(p), min_strain=0.6)
print(1.-(v/v0))









    



[0.   0.02 0.04 0.06 0.08 0.1  0.12 0.14 0.16 0.18 0.2  0.22 0.24 0.26
 0.28 0.3  0.32 0.34]

I cannot quite reproduce the table values. The mismatch is about 1 GPa at 2500 K and 240 GPa.
This means his parameters may have been rounded.
Therefore, I readjusted the eos parameters from Tsuchiya to match their table values better. Users have a choice if they use the table values or the parameter values. If reproduce_table sets to True, the difference reduces to 0.1 GPa.



In [17]:

    
tsuchiya_au = eos.gold.Tsuchiya2003(reproduce_table=True)



In [18]:

    
p = tsuchiya_au.cal_p(v, temp * np.ones_like(v))



In [19]:

    
print('for T = ', temp)
for eta_i, p_i in zip(eta, p):
    print("{0: .3f} {1: .2f}".format(eta_i, p_i))









    



for T =  2500.0
 0.000  16.94+/-0.00
 0.020  20.12+/-0.00
 0.040  23.83+/-0.00
 0.060  28.16+/-0.00
 0.080  33.19+/-0.00
 0.100  39.03+/-0.00
 0.120  45.78+/-0.00
 0.140  53.57+/-0.01
 0.160  62.58+/-0.01
 0.180  72.96+/-0.01
 0.200  84.93+/-0.01
 0.220  98.73+/-0.01
 0.240  114.65+/-0.01
 0.260  133.01+/-0.01
 0.280  154.22+/-0.01
 0.300  178.73+/-0.01
 0.320  207.10+/-0.02
 0.340  239.98+/-0.02



In [ ]:



In [ ]: