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

  • There exist different formulations for the Gruneisen parameter and thermal pressure. In this notebook we will compare the different formulations for MgO.

1. General setup


In [4]:
import uncertainties as uct
import numpy as np
from uncertainties import unumpy as unp
import matplotlib.pyplot as plt
import pytheos as eos

2. Calculate Gruneisen parameter

We test three different ways to get Gruneisen parameter at high pressure-temperature.


In [5]:
help(eos.constq_grun)


Help on function constq_grun in module pytheos.eqn_therm_constq:

constq_grun(v, v0, gamma0, q)
    calculate Gruneisen parameter for constant q
    
    :param v: unit-cell volume in A^3
    :param v0: unit-cell volume in A^3 at 1 bar
    :param gamma0: Gruneisen parameter at 1 bar
    :param q: logarithmic derivative of Grunseinen parameter
    :return: Gruneisen parameter at a given volume


In [6]:
help(eos.speziale_grun)


Help on function speziale_grun in module pytheos.eqn_therm_Speziale:

speziale_grun(v, v0, gamma0, q0, q1)
    calculate Gruneisen parameter for the Speziale equation
    
    :param v: unit-cell volume in A^3
    :param v0: unit-cell volume in A^3 at 1 bar
    :param gamma0: Gruneisen parameter at 1 bar
    :param q0: logarithmic derivative of Gruneisen parameter
    :param q1: logarithmic derivative of Gruneisen parameter
    :return: Gruneisen parameter


In [7]:
help(eos.altshuler_grun)


Help on function altshuler_grun in module pytheos.eqn_therm_Dorogokupets2007:

altshuler_grun(v, v0, gamma0, gamma_inf, beta)
    calculate Gruneisen parameter for Altshuler equation
    
    :param v: unit-cell volume in A^3
    :param v0: unit-cell volume in A^3 at 1 bar
    :param gamma0: Gruneisen parameter at 1 bar
    :param gamma_inf: Gruneisen parameter at infinite pressure
    :param beta: volume dependence of Gruneisen parameter
    :return: Gruneisen parameter

We will test for MgO.


In [8]:
v0 = 74.698
v = np.linspace(v0, v0*0.8, 10)
  • From Dewaele et al. (2000, JGR):
$$\dfrac{\gamma}{\gamma_0} = \left(\dfrac{V}{V_0}\right)^q$$

where $\gamma_0 = 1.45$ and $q = 0.8 \pm 0.5$.

  • From Speziale et al. (2001, JGR):
$$\gamma = \gamma_0 \exp\left\{ \dfrac{q_0}{q_1} \left[ \left(\dfrac{V}{V_0}\right)^{q_1} -1 \right] \right\}$$

where $\gamma_0 = 1.49 \pm 0.03$, $q_0 = 1.65 \pm 0.4$, and $q_1 = 11.8 \pm 0.2$.

  • From Dorogokupets and Dewaele (2007, HPR):
$$q = \beta x^\beta \dfrac{\gamma_0 - \gamma_\inf}{\gamma} $$$$\gamma = \gamma_0 x^\beta$$

where $\gamma_0 = 1.50$, $\gamma_\inf = 0.75$, and $\beta = 2.96$.


In [9]:
gamma_de = eos.constq_grun(v, v0, 1.45, 0.8)
gamma_sp = eos.speziale_grun(v, v0, uct.ufloat(1.49, 0.03), uct.ufloat(1.65, 0.4), uct.ufloat(11.8, 0.2))
gamma_do = eos.altshuler_grun(v, v0, 1.50, 0.75, 2.96)

In [10]:
plt.plot(v, gamma_de, label='Dewaele2000')
plt.errorbar(v, unp.nominal_values(gamma_sp), yerr=unp.std_devs(gamma_sp), label='Speziale2001')
plt.plot(v, gamma_do, label='Dorogokupets2007')
plt.xlabel('Unit-cell volume ($\mathrm{\AA}^3$)')
plt.ylabel(r"$\mathrm{Gr{\"u}neisen parameter}$")
plt.legend();


3. Calculate Debye temperature


In [11]:
help(eos.constq_debyetemp)


Help on function constq_debyetemp in module pytheos.eqn_therm_constq:

constq_debyetemp(v, v0, gamma0, q, theta0)
    calculate Debye temperature for constant q
    
    :param v: unit-cell volume in A^3
    :param v0: unit-cell volume in A^3 at 1 bar
    :param gamma0: Gruneisen parameter at 1 bar
    :param q: logarithmic derivative of Gruneisen parameter
    :param theta0: Debye temperature at 1 bar
    :return: Debye temperature in K


In [12]:
help(eos.speziale_debyetemp)


Help on function speziale_debyetemp in module pytheos.eqn_therm_Speziale:

speziale_debyetemp(v, v0, gamma0, q0, q1, theta0)
    calculate Debye temperature for the Speziale equation
    
    :param v: unit-cell volume in A^3
    :param v0: unit-cell volume in A^3 at 1 bar
    :param gamma0: Gruneisen parameter at 1 bar
    :param q0: logarithmic derivative of Gruneisen parameter
    :param q1: logarithmic derivative of Gruneisen parameter
    :param theta0: Debye temperature at 1 bar in K
    :return: Debye temperature in K


In [13]:
help(eos.altshuler_debyetemp)


Help on function altshuler_debyetemp in module pytheos.eqn_therm_Dorogokupets2007:

altshuler_debyetemp(v, v0, gamma0, gamma_inf, beta, theta0)
    calculate Debye temperature for Altshuler equation
    
    :param v: unit-cell volume in A^3
    :param v0: unit-cell volume in A^3 at 1 bar
    :param gamma0: Gruneisen parameter at 1 bar
    :param gamma_inf: Gruneisen parameter at infinite pressure
    :param beta: volume dependence of Gruneisen parameter
    :param theta0: Debye temperature at 1 bar in K
    :return: Debye temperature in K


In [14]:
theta_de = eos.constq_debyetemp(v, v0, 1.45, 0.8, 800)
theta_sp = eos.speziale_debyetemp(v, v0, uct.ufloat(1.49, 0.03), uct.ufloat(1.65, 0.4), uct.ufloat(11.8, 0.2), 773)
theta_do = eos.altshuler_debyetemp(v, v0, 1.50, 0.75, 2.96, 760)

In [15]:
plt.plot(v, theta_de, label='Dewaele2000')
plt.errorbar(v, unp.nominal_values(theta_sp), yerr=unp.std_devs(theta_sp), label='Speziale2001')
plt.plot(v, theta_do, label='Dorogokupets2007')
plt.xlabel('Unit-cell volume ($\mathrm{\AA}^3$)')
plt.ylabel("Debye Temperature (K)")
plt.legend();


4. Calculate thermal pressure


In [16]:
help(eos.constq_pth)


Help on function constq_pth in module pytheos.eqn_therm_constq:

constq_pth(v, temp, v0, gamma0, q, theta0, n, z, t_ref=300.0, three_r=24.943379399999998)
    calculate thermal pressure for constant q
    
    :param v: unit-cell volume in A^3
    :param temp: temperature
    :param v0: unit-cell volume in A^3 at 1 bar
    :param gamma0: Gruneisen parameter at 1 bar
    :param q: logarithmic derivative of Gruneisen parameter
    :param theta0: Debye temperature in K
    :param n: number of atoms in a formula unit
    :param z: number of formula unit in a unit cell
    :param t_ref: reference temperature
    :param three_r: 3R in case adjustment is needed
    :return: thermal pressure in GPa


In [17]:
help(eos.speziale_pth)


Help on function speziale_pth in module pytheos.eqn_therm_Speziale:

speziale_pth(v, temp, v0, gamma0, q0, q1, theta0, n, z, t_ref=300.0, three_r=24.943379399999998)
    calculate thermal pressure for the Speziale equation
    
    :param v: unit-cell volume in A^3
    :param temp: temperature in K
    :param v0: unit-cell volume in A^3 at 1 bar
    :param gamma0: Gruneisen parameter at 1 bar
    :param q0: logarithmic derivative of Gruneisen parameter
    :param q1: logarithmic derivative of Gruneisen parameter
    :param theta0: Debye temperature at 1 bar in K
    :param n: number of atoms in a formula unit
    :param z: number of formula unit in a unit cell
    :param t_ref: reference temperature
    :param three_r: 3R in case adjustment is needed
    :return: thermal pressure in GPa


In [18]:
help(eos.dorogokupets2007_pth)


Help on function dorogokupets2007_pth in module pytheos.eqn_therm_Dorogokupets2007:

dorogokupets2007_pth(v, temp, v0, gamma0, gamma_inf, beta, theta0, n, z, three_r=24.943379399999998, t_ref=300.0)
    calculate thermal pressure for Dorogokupets 2007 EOS
    
    :param v: unit-cell volume in A^3
    :param temp: temperature in K
    :param v0: unit-cell volume in A^3 at 1 bar
    :param gamma0: Gruneisen parameter at 1 bar
    :param gamma_inf: Gruneisen parameter at infinite pressure
    :param beta: volume dependence of Gruneisen parameter
    :param theta0: Debye temperature at 1 bar in K
    :param n: number of elements in a chemical formula
    :param z: number of formula unit in a unit cell
    :param three_r: 3 times gas constant.
        Jamieson modified this value to compensate for mismatches
    :param t_ref: reference temperature, 300 K
    :return: thermal pressure in GPa


In [19]:
n = 2; z = 4
temp = np.ones_like(v) * 2000.

In [20]:
pth_de = eos.constq_pth(v, temp, v0, 1.45, 0.8, 800, n, z)
pth_sp = eos.speziale_pth(v, temp, v0, uct.ufloat(1.49, 0.03), uct.ufloat(1.65, 0.4), uct.ufloat(11.8, 0.2), 773, n, z)
pth_do = eos.dorogokupets2007_pth(v, temp, v0, 1.50, 0.75, 2.96, 760, n, z)

In [21]:
plt.plot(v, pth_de, label='Dewaele2000')
plt.errorbar(v, unp.nominal_values(pth_sp), yerr=unp.std_devs(pth_sp), label='Speziale2001')
plt.plot(v, pth_do, label='Dorogokupets2007')
plt.xlabel('Unit-cell volume ($\mathrm{\AA}^3$)')
plt.ylabel("Thermal pressure (GPa)")
plt.legend();