Skip to content

Commit

Permalink
added SPOCK paper scripts
Browse files Browse the repository at this point in the history
  • Loading branch information
bobmyhill committed Nov 28, 2024
1 parent d26c6bb commit e535b7e
Show file tree
Hide file tree
Showing 5 changed files with 712 additions and 0 deletions.
221 changes: 221 additions & 0 deletions contrib/spock_eos/Au_comparison.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,221 @@
from scipy.special import gamma, gammainc, exp1
import numpy as np
import matplotlib.pyplot as plt
from burnman.eos.reciprocal_kprime import RKprime
from burnman.eos.birch_murnaghan import BM3
from burnman.eos.birch_murnaghan_4th import BM4
from burnman.eos.vinet import Vinet
from burnman.eos.modified_tait import MT
from burnman.eos.spock import SPOCK
from burnman.eos.macaw import MACAW, make_params
from burnman import Mineral
from burnman.tools.eos import check_eos_consistency
import warnings

"""
Au_comparison
-------------
"""

# Define the equations of state we wish to compare
eoses = {
"Vinet": Vinet(),
"BM3": BM3(),
"BM4": BM4(),
"MT": MT(),
"RK": RKprime(),
"MACAW": MACAW(),
"SPOCK": SPOCK(),
}

# Create the parameter dictionary with the required values for
# all the equations of state.
params = {
"E_0": 0.0,
"P_0": 1.0e5,
"T_0": 0.0,
"V_0": 10.1959e-6,
"K_0": 156.6e9,
"Kprime_0": 6.76,
"Kprime_prime_0": -15.6 / 156.6e9,
"Kdprime_0": -15.6 / 156.6e9,
"Kprime_inf": 2.97,
"molar_mass": 0.06008,
"equation_of_state": None,
}


# Define the range of volumes over which to
# compare the equations of state
volumes = np.linspace(1.0e-6, 1.4, 1001) * params["V_0"]

# Initialise the figure and subplots
fig = plt.figure(figsize=(10, 7))
ax = [fig.add_subplot(2, 2, i) for i in [4, 2, 1, 3]]

for name, eos in eoses.items():

# Set the equation of state and initialise the mineral object
params["equation_of_state"] = eos
m = Mineral(params)

# Check that the SPOCK and MACAW equations of state are internally consistent
if name == "SPOCK" or name == "MACAW":
assert check_eos_consistency(m, 1.0e10, 300.0, including_shear_properties=False)

# Define some plotting values
if name == "SPOCK":
linestyle = "-"
linewidth = 2.0
alpha = 1.0
else:
linestyle = "--"
linewidth = 1.0
alpha = 1.0

# Get the values of the properties to plot from the equations of state
pressures = np.empty_like(volumes)
KTs = np.empty_like(pressures)
Vs = np.empty_like(pressures)
Es = np.empty_like(pressures)
for i, volume in enumerate(volumes):
with warnings.catch_warnings(record=True) as w:
try:
Vs[i] = volumes[i]
pressures[i] = eos.pressure(0.0, volumes[i], params)
KTs[i] = eos.isothermal_bulk_modulus_reuss(
pressures[i], 0.0, volumes[i], params
)
Es[i] = (
eos.molar_internal_energy(pressures[i], 0.0, volumes[i], params) / 1.0e9
)
except ValueError:
pressures[i] = np.nan
KTs[i] = np.nan
Es[i] = np.nan

with warnings.catch_warnings(record=True) as w:
Kprimes = -np.gradient(np.log(KTs), np.log(Vs), edge_order=2)

# Plot the properties
ax[0].plot(
Vs / params["V_0"],
Kprimes,
linestyle=linestyle,
linewidth=linewidth,
alpha=alpha,
label=name,
)
ax[1].plot(
Vs / params["V_0"],
KTs / 1.0e9,
alpha=alpha,
linestyle=linestyle,
linewidth=linewidth,
label=name,
)
ax[2].plot(
Vs / params["V_0"],
pressures / 1.0e9,
alpha=alpha,
linestyle=linestyle,
linewidth=linewidth,
label=name,
)
ax[3].plot(
Vs / params["V_0"],
Es * 1.0e9,
alpha=alpha,
linestyle=linestyle,
linewidth=linewidth,
label=name,
)


def Kdprime_0K_0_MACAW(params):
A, B, C = make_params(params["K_0"], params["Kprime_0"], params["Kprime_inf"])
return (
C
* (
-3.0 * B**4
- 12.0 * B**3 * C
+ 3.0 * B**3
- 18.0 * B**2 * C**2
+ 9.0 * B**2 * C
+ 3.75 * B**2
- 12.0 * B * C**3
+ 9.0 * B * C**2
+ 21.0 * B * C
- 2.25 * B
- 3.0 * C**4
+ 3.0 * C**3
- 3.0 * C**2
)
/ (
2.0 * B**4
+ 8.0 * B**3 * C
+ 4.0 * B**3
+ 12.0 * B**2 * C**2
+ 6.0 * B**2 * C
+ 2.0 * B**2
+ 8.0 * B * C**3
- 2.0 * B * C
+ 2.0 * C**4
- 2.0 * C**3
+ 0.5 * C**2
)
)


print(f"MACAW value of K''_0 * K_0: {Kdprime_0K_0_MACAW(params):.2f}")

# Assign useful ranges and labels to the plots
ax[0].set_ylim(0.0, 15.0)
ax[1].set_ylim(1.0, 1.0e6)
ax[2].set_ylim(-50.0, 1000.0)
ax[3].set_ylim(1.0e3, 1.0e9)
ax[1].set_yscale("log")
ax[3].set_yscale("log")

for i, param in enumerate(["Kprime_0", "K_0", "P_0", "E_0"]):
p = params[param]
if i == 1 or i == 2:
p = p / 1.0e9
ax[i].set_xlim(0.0, 1.4)
ax[i].plot(ax[i].get_xlim(), [p, p], c="grey", linewidth=0.5)
ax[i].set_xlabel("$V/V_0$")

for i, param in enumerate(["Kprime_inf"]):
p = params[param]
if i == 1 or i == 2:
p = p / 1.0e9
ax[i].plot(ax[i].get_xlim(), [p, p], c="grey", linewidth=1.0)

for i in range(4):
ax[i].set_ylim(ax[i].get_ylim())
ax[i].plot([1.0, 1.0], ax[i].get_ylim(), c="grey", linewidth=0.5)

ax[0].set_ylabel("$K'$")
ax[1].set_ylabel("$K$ (GPa)")
ax[2].set_ylabel("Pressure (GPa)")
ax[3].set_ylabel("Internal energy (J/mol)")

bbox = dict(facecolor="white", edgecolor="none", alpha=1.0)
for i, label in enumerate(["d", "b", "a", "c"]):
ax[i].text(
0.04,
0.93,
label,
fontsize=12,
ha="center",
va="center",
transform=ax[i].transAxes,
bbox=bbox,
)

ax[2].legend()
fig.set_tight_layout(True)

# Save and show the figure
fig.savefig("figures/Au_comparison.pdf")
plt.show()
23 changes: 23 additions & 0 deletions contrib/spock_eos/Readme.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,23 @@
## Companion code to
# The SPOCK equation of state for condensed phases under arbitrary compression
## R. Myhill

The python scripts contained in this directory accompany the paper
submitted by Myhill (2025).
There are comments in each of the scripts, and it is hoped that by reading
these in tandem with the original paper, the reader will get a feel
for how to use the code for their own purposes.

The scripts contained in this directory are as follows:

fit_spock.py
------------
This script fits various equations of state to the Au and Pt data of
Fratanduono et al., 2021. Plots Figures 1 and 2 from the paper, as well
as corner plots demonstrating the parameter uncertainties.

Au_comparison.py
----------------
This script compares various equations of state given the standard state
and infinite pressure parameters optimised for the Au SPOCK equation of state.
Plots Figure 3.
120 changes: 120 additions & 0 deletions contrib/spock_eos/data/Au.dat
Original file line number Diff line number Diff line change
@@ -0,0 +1,120 @@
# Density (g/cm^3) Volume (A^3/at) Stress (and err, GPa) Isentrope (and err, GPa) Isotherm (and err, GPa)
19.3200 16.9290 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
19.5200 16.7560 1.9000 0.1000 1.8000 0.1000 1.7000 0.1000
19.7200 16.5860 3.9000 0.2000 3.8000 0.2000 3.5000 0.2000
19.9200 16.4190 6.0000 0.2000 5.8000 0.2000 5.3000 0.2000
20.1200 16.2560 8.1000 0.2000 8.0000 0.2000 7.3000 0.3000
20.3200 16.0960 10.4000 0.2000 10.2000 0.3000 9.4000 0.3000
20.5200 15.9390 12.7000 0.3000 12.5000 0.3000 11.5000 0.4000
20.7200 15.7850 15.2000 0.3000 15.0000 0.3000 13.8000 0.4000
20.9200 15.6340 17.8000 0.3000 17.5000 0.4000 16.3000 0.4000
21.1200 15.4860 20.4000 0.3000 20.2000 0.4000 18.8000 0.5000
21.3200 15.3410 23.2000 0.4000 22.9000 0.4000 21.4000 0.5000
21.5200 15.1980 26.0000 0.4000 25.7000 0.4000 24.1000 0.5000
21.7200 15.0590 29.0000 0.4000 28.7000 0.5000 26.9000 0.6000
21.9200 14.9210 32.1000 0.5000 31.8000 0.5000 29.9000 0.6000
22.1200 14.7860 35.4000 0.5000 35.0000 0.5000 33.1000 0.7000
22.3200 14.6540 38.7000 0.5000 38.2000 0.6000 36.2000 0.7000
22.5200 14.5240 42.1000 0.6000 41.6000 0.6000 39.5000 0.8000
22.7200 14.3960 45.6000 0.7000 45.1000 0.7000 42.9000 0.8000
22.9200 14.2700 49.3000 0.7000 48.8000 0.8000 46.4000 0.9000
23.1200 14.1470 53.1000 0.8000 52.5000 0.8000 50.1000 1.0000
23.3200 14.0250 57.0000 0.9000 56.4000 0.9000 54.0000 1.1000
23.5200 13.9060 61.1000 1.0000 60.4000 1.0000 57.9000 1.2000
23.7200 13.7890 65.2000 1.1000 64.6000 1.2000 61.9000 1.3000
23.9200 13.6740 69.5000 1.2000 68.8000 1.3000 66.1000 1.4000
24.1200 13.5600 74.0000 1.4000 73.2000 1.5000 70.4000 1.6000
24.3200 13.4490 78.6000 1.5000 77.8000 1.6000 74.9000 1.7000
24.5200 13.3390 83.2000 1.7000 82.4000 1.8000 79.4000 1.9000
24.7200 13.2310 88.1000 1.8000 87.1000 1.9000 84.1000 2.0000
24.9200 13.1250 93.0000 1.9000 92.0000 2.0000 88.9000 2.1000
25.1200 13.0200 98.0000 2.0000 97.0000 2.2000 93.9000 2.3000
25.3200 12.9170 103.2000 2.1000 102.2000 2.3000 98.9000 2.4000
25.5200 12.8160 108.5000 2.3000 107.4000 2.5000 104.1000 2.5000
25.7200 12.7170 114.0000 2.5000 112.8000 2.6000 109.4000 2.7000
25.9200 12.6180 119.6000 2.6000 118.4000 2.8000 114.9000 2.9000
26.1200 12.5220 125.4000 2.8000 124.1000 3.0000 120.5000 3.1000
26.3200 12.4270 131.3000 2.9000 129.9000 3.1000 126.3000 3.2000
26.5200 12.3330 137.2000 3.0000 135.9000 3.3000 132.2000 3.3000
26.7200 12.2410 143.3000 3.2000 141.9000 3.4000 138.1000 3.5000
26.9200 12.1500 149.6000 3.3000 148.1000 3.5000 144.2000 3.6000
27.1200 12.0600 155.9000 3.5000 154.4000 3.7000 150.4000 3.8000
27.3200 11.9720 162.4000 3.6000 160.8000 3.9000 156.8000 4.0000
27.5200 11.8850 169.0000 3.8000 167.3000 4.1000 163.3000 4.1000
27.7200 11.7990 175.8000 4.0000 174.0000 4.3000 169.9000 4.3000
27.9200 11.7150 182.7000 4.2000 180.9000 4.5000 176.7000 4.5000
28.1200 11.6310 189.8000 4.4000 187.9000 4.7000 183.6000 4.8000
28.3200 11.5490 197.0000 4.6000 195.0000 5.0000 190.7000 5.0000
28.5200 11.4680 204.4000 4.9000 202.4000 5.2000 198.0000 5.3000
28.7200 11.3880 211.9000 5.2000 209.8000 5.5000 205.4000 5.6000
28.9200 11.3090 219.6000 5.4000 217.4000 5.8000 212.9000 5.9000
29.1200 11.2320 227.5000 5.7000 225.2000 6.1000 220.6000 6.2000
29.3200 11.1550 235.5000 6.0000 233.1000 6.4000 228.5000 6.5000
29.5200 11.0800 243.6000 6.4000 241.2000 6.8000 236.5000 6.9000
29.7200 11.0050 251.9000 6.7000 249.4000 7.2000 244.6000 7.2000
29.9200 10.9320 260.4000 7.1000 257.8000 7.5000 252.9000 7.6000
30.1200 10.8590 269.0000 7.4000 266.3000 8.0000 261.4000 8.0000
30.3200 10.7870 277.7000 7.8000 274.9000 8.3000 269.9000 8.4000
30.5200 10.7170 286.6000 8.2000 283.7000 8.7000 278.6000 8.8000
30.7200 10.6470 295.6000 8.6000 292.7000 9.1000 287.5000 9.2000
30.9200 10.5780 304.7000 8.9000 301.7000 9.5000 296.5000 9.6000
31.1200 10.5100 314.0000 9.4000 310.9000 10.0000 305.7000 10.0000
31.3200 10.4430 323.5000 9.8000 320.3000 10.4000 314.9000 10.5000
31.5200 10.3770 333.0000 10.2000 329.7000 10.8000 324.3000 10.9000
31.7200 10.3110 342.7000 10.6000 339.3000 11.3000 333.9000 11.4000
31.9200 10.2470 352.5000 11.1000 349.0000 11.8000 343.5000 11.9000
32.1200 10.1830 362.5000 11.6000 358.9000 12.3000 353.3000 12.4000
32.3200 10.1200 372.6000 12.0000 368.9000 12.9000 363.3000 12.9000
32.5200 10.0580 382.9000 12.6000 379.1000 13.5000 373.3000 13.5000
32.7200 9.9960 393.3000 13.3000 389.4000 14.2000 383.6000 14.2000
32.9200 9.9350 403.9000 13.9000 399.8000 14.8000 394.0000 14.8000
33.1200 9.8750 414.6000 14.6000 410.5000 15.6000 404.6000 15.6000
33.3200 9.8160 425.5000 15.3000 421.3000 16.4000 415.3000 16.4000
33.5200 9.7570 436.5000 16.0000 432.1000 17.0000 426.1000 17.0000
33.7200 9.7000 447.5000 16.8000 443.1000 18.0000 437.0000 18.0000
33.9200 9.6420 458.7000 17.8000 454.2000 19.0000 448.0000 19.0000
34.1200 9.5860 470.2000 19.0000 465.5000 20.2000 459.3000 20.2000
34.3200 9.5300 482.1000 19.9000 477.3000 21.2000 471.0000 21.2000
34.5200 9.4750 493.8000 20.2000 488.9000 21.5000 482.6000 21.5000
34.7200 9.4200 505.7000 20.7000 500.7000 22.1000 494.3000 22.1000
34.9200 9.3660 517.7000 21.2000 512.5000 22.6000 506.1000 22.6000
35.1200 9.3130 529.8000 21.8000 524.5000 23.2000 518.0000 23.2000
35.3200 9.2600 542.0000 22.3000 536.6000 23.8000 530.0000 23.8000
35.5200 9.2080 554.4000 23.0000 548.9000 24.5000 542.2000 24.5000
35.7200 9.1570 566.8000 23.6000 561.2000 25.1000 554.5000 25.2000
35.9200 9.1060 579.5000 24.3000 573.7000 25.9000 566.9000 25.9000
36.1200 9.0550 592.3000 25.1000 586.5000 26.7000 579.6000 26.8000
36.3200 9.0050 605.4000 25.8000 599.4000 27.4000 592.5000 27.5000
36.5200 8.9560 618.5000 26.4000 612.4000 28.1000 605.4000 28.1000
36.7200 8.9070 631.8000 26.9000 625.5000 28.6000 618.5000 28.6000
36.9200 8.8590 645.2000 27.6000 638.8000 29.4000 631.7000 29.4000
37.1200 8.8110 658.7000 28.2000 652.2000 30.0000 645.0000 30.0000
37.3200 8.7640 672.3000 28.8000 665.7000 30.7000 658.5000 30.7000
37.5200 8.7170 686.5000 30.7000 679.7000 32.7000 672.4000 32.8000
37.7200 8.6710 701.4000 32.1000 694.4000 34.2000 687.1000 34.2000
37.9200 8.6250 716.4000 32.9000 709.3000 35.0000 701.9000 35.0000
38.1200 8.5800 731.6000 33.6000 724.4000 35.8000 716.9000 35.8000
38.3200 8.5350 746.9000 34.4000 739.5000 36.6000 732.0000 36.6000
38.5200 8.4910 762.4000 35.1000 754.8000 37.3000 747.3000 37.4000
38.7200 8.4470 777.9000 35.7000 770.2000 38.1000 762.6000 38.1000
38.9200 8.4040 793.6000 36.5000 785.7000 38.8000 778.1000 38.9000
39.1200 8.3610 809.4000 37.3000 801.4000 39.7000 793.7000 39.7000
39.3200 8.3180 825.4000 38.1000 817.3000 40.6000 809.5000 40.6000
39.5200 8.2760 841.6000 38.9000 833.3000 41.5000 825.4000 41.5000
39.7200 8.2340 857.9000 39.8000 849.4000 42.4000 841.5000 42.4000
39.9200 8.1930 874.6000 41.5000 866.0000 44.2000 858.0000 44.2000
40.1200 8.1520 891.9000 42.9000 883.1000 45.7000 875.0000 45.7000
40.3200 8.1120 909.3000 44.2000 900.4000 47.0000 892.3000 47.0000
40.5200 8.0720 927.0000 45.5000 917.9000 48.4000 909.7000 48.4000
40.7200 8.0320 944.9000 46.7000 935.6000 49.7000 927.4000 49.7000
40.9200 7.9930 963.0000 47.9000 953.6000 51.0000 945.3000 51.0000
41.1200 7.9540 981.4000 49.2000 971.7000 52.4000 963.4000 52.4000
41.3200 7.9160 999.9000 50.4000 990.0000 53.6000 981.7000 53.6000
41.5200 7.8770 1018.5000 51.6000 1008.5000 55.0000 1000.1000 55.0000
41.7200 7.8400 1037.5000 53.2000 1027.3000 56.7000 1018.8000 56.7000
41.9200 7.8020 1056.8000 55.0000 1046.4000 58.6000 1037.9000 58.6000
42.1200 7.7650 1076.4000 56.4000 1065.8000 60.1000 1057.2000 60.1000
42.3200 7.7290 1096.1000 57.8000 1085.4000 61.5000 1076.7000 61.5000
42.5200 7.6920 1116.1000 59.3000 1105.2000 63.1000 1096.5000 63.1000
42.7200 7.6560 1136.6000 61.9000 1125.5000 65.8000 1116.7000 65.9000
42.9200 7.6200 1157.7000 64.8000 1146.4000 68.9000 1137.6000 69.0000
Loading

0 comments on commit e535b7e

Please sign in to comment.