```python
#%matplotlib ipympl
%matplotlib inline
```

For high dpi displays.


```python
%config InlineBackend.figure_format = 'retina'
```

# 0. Note

This notebook demonstrates how different physical parameters change the equation of state.

# 1. Setup


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

In the Mie-Gruneisen approach, we divide total pressure, $P(V,T)$, into static pressure, $P(V,T_0)$, which is a reference isotherm (normally $T_0 = 300$ K), and thermal pressure, $\Delta P_\mathrm{th}(V,T)$, which is a pressure difference between 300-K and high-temperature isotherms.

\begin{equation}
P(V,T) = P(V,T_0) + \Delta P_\mathrm{th}(V,T).\label{eq:eos-totalP}
\end{equation}

# 2. Static pressure

We will use the Birch-Murnathan equation for the static part.

\begin{equation}
P = \frac{3K_0}{2} \left[ \left(\frac{V_0}{V}\right)^{7/3} -
\left(\frac{V_0}{V}\right)^{5/3} \right]\left\{
1-\xi\left[\left(\frac{V_0}{V}\right)^{2/3} -1\right] \right\},
\end{equation}
where $K_0$ is the bulk modulus (in GPa), $K'_0 = (dK/dP)$, and $\xi = \frac{3}{4}(4-K_0')$.


Generate volume strain: $\eta = V/V_0$.


```python
eta = np.linspace(1., 0.75, 21)
print(eta)
```

Suppose you have a material with the following properties measured at 300 K.


```python
v0 = 163. # A^3
k0 = ufloat(260., 3.) # GPa
k0p = ufloat(4., 0.2)  
```


```python
v = unp.uarray(eta * v0, 0.05)
```

Then you can calculate pressures for the volumes at 300 K: a 300-K isotherm.


```python
p_st = eos.bm3_p(v, v0, k0, k0p)
print(p_st)
```


```python
plt.errorbar(unp.nominal_values(p_st), unp.nominal_values(v), \
            xerr=unp.std_devs(p_st), yerr=unp.std_devs(v)) 
plt.xlabel('Pressure')
plt.ylabel('Volume')
```

## Try these...

- Increase bulk modulus, $K_0$, to 300 GPa.

- Increase pressure derivative of bulk modulus, $K'_0$, to 3 and 5.

What did you find?

Can you convert this into a density profile along an isotherm?

# 3. Thermal pressure

Now let's calculate thermal pressure.  We will use the Debye approach.

\begin{equation}
\Delta P_\mathrm{th}(V,T) = \frac{\gamma(V)}{V}\Delta E_\mathrm{th}[\theta_D(V),T],\label{eq:eos-MGD}
\end{equation}
where $\gamma$ is the Gruneisen parameter ($= V (\partial P/\partial E)_V$), and $T_0=300$ K in most experiments.  The internal energy change can be calculated from the Debye model:

\begin{equation}
E_\mathrm{th} = \frac{9nRT}{x^3}\int^x _0 \frac{\xi^3}{e^\xi -1} d\xi,
\end{equation}
where $R$ is the gas constant, $x = \dfrac{\theta_D}{T}$, and $n$ is the number of atoms per formula unit.

Also,

\begin{equation}
\gamma = \gamma_0 \left( \dfrac{V}{V_0} \right) ^q
\end{equation}



We need some parameters for the Debye thermal pressure equation.


```python
temp = 2000.
gamma0 = ufloat(1.5, 0.3)
q = ufloat(1., 0.5)
theta0 = ufloat(1000., 50.)
```


```python
n = 5.
z = 4.
```


```python
p_th = eos.constq_pth(v, temp, v0, gamma0, q, theta0, n, z)
```


```python
p_ht = p_st + p_th
print(p_ht)
```

Let's plot 300-K and 3000-K isotherms together.


```python
plt.errorbar(unp.nominal_values(p_st), unp.nominal_values(v), xerr=unp.std_devs(p_st))
plt.errorbar(unp.nominal_values(p_ht), unp.nominal_values(v), xerr=unp.std_devs(p_ht))
plt.xlabel('Pressure')
plt.ylabel('Volume')
```

## Try these...

- Increase Gruneisen parameter, $\gamma_0$, to 3.

- Increase logarithmic volume dependence of $\gamma$, $q$, to 0.1 and 2.

What did you find?

Can you convert this into a density profile along an isotherm?

# 4. Volume expansion


```python
import pytheos.scales.objs
from collections import OrderedDict
```


```python
params_st = OrderedDict([('v0', ufloat(v0, 0.01)),
                         ('k0', k0),
                         ('k0p', k0p)])
params_th = OrderedDict([('v0', ufloat(v0, 0.01)),
                         ('gamma0', gamma0),
                         ('q', q),
                         ('theta0', theta0)])
brd = pytheos.scales.objs.MGEOS(n, z, params_st=params_st, params_th=params_th,
               eqn_st='bm3', eqn_th='constq')
```


```python
p_lm = np.linspace(25., 125., 101)
```


```python
v_r = brd.cal_v(p_lm, np.ones_like(p_lm)*1500.)
v_ht = brd.cal_v(p_lm, np.ones_like(p_lm)*2000.)
```


```python
plt.plot(p_lm, (v_ht - v_r) / v_r * 100.)

plt.ylabel('$\Delta V / V$ (%)')
plt.xlabel('Pressure (GPa)')
```

e18e838461655991769d1aa1b0625cfd8be359ea
ipynb
Jupyter Notebook
0-eos.ipynb
SHDShim/CIDER2018-ex
f1d6c3b0b1f5c3b96f5e0b61ff18ca9045377752
Apache-2.0
2018-07-19T19:35:49.000Z
2020-12-02T01:06:48.000Z
0-eos.ipynb
SHDShim/CIDER2018-ex
f1d6c3b0b1f5c3b96f5e0b61ff18ca9045377752
Apache-2.0
0-eos.ipynb
SHDShim/CIDER2018-ex
f1d6c3b0b1f5c3b96f5e0b61ff18ca9045377752
Apache-2.0
Qwen/Qwen-72B
1. YES
2. YES