Skip to content

Commit e21f059

Browse files
bobmyhillbobmyhill
authored
Merge pull request #618 from bobmyhill/spock_eos
Spock EOS paper
2 parents aaaf264 + b7ca85e commit e21f059

File tree

6 files changed

+786
-0
lines changed

6 files changed

+786
-0
lines changed

‎contrib/spock_eos/Au_comparison.py‎

Lines changed: 240 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,240 @@
1+
from scipy.special import gamma, gammainc, exp1
2+
import numpy as np
3+
import matplotlib.pyplot as plt
4+
from burnman.eos.reciprocal_kprime import RKprime
5+
from burnman.eos.birch_murnaghan import BM3
6+
from burnman.eos.birch_murnaghan_4th import BM4
7+
from burnman.eos.vinet import Vinet
8+
from burnman.eos.modified_tait import MT
9+
from burnman.eos.spock import SPOCK
10+
from burnman.eos.macaw import MACAW, make_params
11+
from burnman import Mineral
12+
from burnman.tools.eos import check_eos_consistency
13+
import warnings
14+
15+
"""
16+
Au_comparison
17+
-------------
18+
"""
19+
20+
# Colorblind friendly colors
21+
# https://www.nature.com/articles/nmeth.1618.pdf
22+
colors = {
23+
"SPOCK": "#000000",
24+
"BM3": "#E69F00",
25+
"none": "#56B4E9",
26+
"MT": "#009E73",
27+
"Vinet": "#F0E442",
28+
"MACAW": "#0072B2",
29+
"BM4": "#D55E00",
30+
"RK": "#CC79A7",
31+
}
32+
33+
34+
# Define the equations of state we wish to compare
35+
eoses = {
36+
"Vinet": Vinet(),
37+
"BM3": BM3(),
38+
"BM4": BM4(),
39+
"MT": MT(),
40+
"RK": RKprime(),
41+
"MACAW": MACAW(),
42+
"SPOCK": SPOCK(),
43+
}
44+
45+
# Create the parameter dictionary with the required values for
46+
# all the equations of state.
47+
params = {
48+
"E_0": 0.0,
49+
"P_0": 1.0e5,
50+
"T_0": 0.0,
51+
"V_0": 10.1959e-6,
52+
"K_0": 156.6e9,
53+
"Kprime_0": 6.76,
54+
"Kprime_prime_0": -15.6 / 156.6e9,
55+
"Kdprime_0": -15.6 / 156.6e9,
56+
"Kprime_inf": 2.97,
57+
"molar_mass": 0.06008,
58+
"equation_of_state": None,
59+
}
60+
61+
62+
# Define the range of volumes over which to
63+
# compare the equations of state
64+
volumes = np.linspace(1.0e-6, 1.4, 1001) * params["V_0"]
65+
66+
# Initialise the figure and subplots
67+
fig = plt.figure(figsize=(10, 7))
68+
ax = [fig.add_subplot(2, 2, i) for i in [4, 2, 1, 3]]
69+
70+
for name, eos in eoses.items():
71+
72+
# Set the equation of state and initialise the mineral object
73+
params["equation_of_state"] = eos
74+
m = Mineral(params)
75+
76+
# Check that the SPOCK and MACAW equations of state are internally consistent
77+
if name == "SPOCK" or name == "MACAW":
78+
assert check_eos_consistency(m, 1.0e10, 300.0, including_shear_properties=False)
79+
80+
# Define some plotting values
81+
if name == "SPOCK":
82+
linestyle = "-"
83+
linewidth = 2.0
84+
alpha = 1.0
85+
else:
86+
linestyle = "-"
87+
linewidth = 1.5
88+
alpha = 1.0
89+
90+
# Get the values of the properties to plot from the equations of state
91+
pressures = np.empty_like(volumes)
92+
KTs = np.empty_like(pressures)
93+
Vs = np.empty_like(pressures)
94+
Es = np.empty_like(pressures)
95+
for i, volume in enumerate(volumes):
96+
with warnings.catch_warnings(record=True) as w:
97+
try:
98+
Vs[i] = volumes[i]
99+
pressures[i] = eos.pressure(0.0, volumes[i], params)
100+
KTs[i] = eos.isothermal_bulk_modulus_reuss(
101+
pressures[i], 0.0, volumes[i], params
102+
)
103+
Es[i] = (
104+
eos.molar_internal_energy(pressures[i], 0.0, volumes[i], params)
105+
/ 1.0e9
106+
)
107+
except ValueError:
108+
pressures[i] = np.nan
109+
KTs[i] = np.nan
110+
Es[i] = np.nan
111+
112+
with warnings.catch_warnings(record=True) as w:
113+
Kprimes = -np.gradient(np.log(KTs), np.log(Vs), edge_order=2)
114+
115+
# Plot the properties
116+
ax[0].plot(
117+
Vs / params["V_0"],
118+
Kprimes,
119+
linestyle=linestyle,
120+
linewidth=linewidth,
121+
alpha=alpha,
122+
label=name,
123+
c=colors[name],
124+
)
125+
ax[1].plot(
126+
Vs / params["V_0"],
127+
KTs / 1.0e9,
128+
alpha=alpha,
129+
linestyle=linestyle,
130+
linewidth=linewidth,
131+
label=name,
132+
c=colors[name],
133+
)
134+
ax[2].plot(
135+
Vs / params["V_0"],
136+
pressures / 1.0e9,
137+
alpha=alpha,
138+
linestyle=linestyle,
139+
linewidth=linewidth,
140+
label=name,
141+
c=colors[name],
142+
)
143+
ax[3].plot(
144+
Vs / params["V_0"],
145+
Es * 1.0e9,
146+
alpha=alpha,
147+
linestyle=linestyle,
148+
linewidth=linewidth,
149+
label=name,
150+
c=colors[name],
151+
)
152+
153+
154+
def Kdprime_0K_0_MACAW(params):
155+
A, B, C = make_params(params["K_0"], params["Kprime_0"], params["Kprime_inf"])
156+
return (
157+
C
158+
* (
159+
-3.0 * B**4
160+
- 12.0 * B**3 * C
161+
+ 3.0 * B**3
162+
- 18.0 * B**2 * C**2
163+
+ 9.0 * B**2 * C
164+
+ 3.75 * B**2
165+
- 12.0 * B * C**3
166+
+ 9.0 * B * C**2
167+
+ 21.0 * B * C
168+
- 2.25 * B
169+
- 3.0 * C**4
170+
+ 3.0 * C**3
171+
- 3.0 * C**2
172+
)
173+
/ (
174+
2.0 * B**4
175+
+ 8.0 * B**3 * C
176+
+ 4.0 * B**3
177+
+ 12.0 * B**2 * C**2
178+
+ 6.0 * B**2 * C
179+
+ 2.0 * B**2
180+
+ 8.0 * B * C**3
181+
- 2.0 * B * C
182+
+ 2.0 * C**4
183+
- 2.0 * C**3
184+
+ 0.5 * C**2
185+
)
186+
)
187+
188+
189+
print(f"MACAW value of K''_0 * K_0: {Kdprime_0K_0_MACAW(params):.2f}")
190+
191+
# Assign useful ranges and labels to the plots
192+
ax[0].set_ylim(0.0, 15.0)
193+
ax[1].set_ylim(1.0, 1.0e6)
194+
ax[2].set_ylim(-50.0, 1000.0)
195+
ax[3].set_ylim(1.0e3, 1.0e9)
196+
ax[1].set_yscale("log")
197+
ax[3].set_yscale("log")
198+
199+
for i, param in enumerate(["Kprime_0", "K_0", "P_0", "E_0"]):
200+
p = params[param]
201+
if i == 1 or i == 2:
202+
p = p / 1.0e9
203+
ax[i].set_xlim(0.0, 1.4)
204+
ax[i].plot(ax[i].get_xlim(), [p, p], c="grey", linewidth=0.5)
205+
ax[i].set_xlabel("$V/V_0$")
206+
207+
for i, param in enumerate(["Kprime_inf"]):
208+
p = params[param]
209+
if i == 1 or i == 2:
210+
p = p / 1.0e9
211+
ax[i].plot(ax[i].get_xlim(), [p, p], c="grey", linewidth=1.0)
212+
213+
for i in range(4):
214+
ax[i].set_ylim(ax[i].get_ylim())
215+
ax[i].plot([1.0, 1.0], ax[i].get_ylim(), c="grey", linewidth=0.5)
216+
217+
ax[0].set_ylabel("$K'$")
218+
ax[1].set_ylabel("$K$ (GPa)")
219+
ax[2].set_ylabel("Pressure (GPa)")
220+
ax[3].set_ylabel("Internal energy (J/mol)")
221+
222+
bbox = dict(facecolor="white", edgecolor="none", alpha=1.0)
223+
for i, label in enumerate(["d", "b", "a", "c"]):
224+
ax[i].text(
225+
0.04,
226+
0.93,
227+
label,
228+
fontsize=12,
229+
ha="center",
230+
va="center",
231+
transform=ax[i].transAxes,
232+
bbox=bbox,
233+
)
234+
235+
ax[2].legend()
236+
fig.set_tight_layout(True)
237+
238+
# Save and show the figure
239+
fig.savefig("figures/Au_comparison.pdf")
240+
plt.show()

‎contrib/spock_eos/Readme.md‎

Lines changed: 28 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,28 @@
1+
## Companion code to
2+
# The SPOCK equation of state for condensed phases under arbitrary compression
3+
## R. Myhill
4+
5+
The python scripts contained in this directory accompany the paper
6+
submitted by Myhill (2025).
7+
There are comments in each of the scripts, and it is hoped that by reading
8+
these in tandem with the original paper, the reader will get a feel
9+
for how to use the code for their own purposes.
10+
11+
The scripts contained in this directory are as follows:
12+
13+
fit_spock.py
14+
------------
15+
This script fits various equations of state to the Au and Pt data of
16+
Fratanduono et al., 2021. Plots Figures 1 and 2 from the paper, as well
17+
as corner plots demonstrating the parameter uncertainties.
18+
19+
Au_comparison.py
20+
----------------
21+
This script compares various equations of state given the standard state
22+
and infinite pressure parameters optimised for the Au SPOCK equation of state.
23+
Plots Figure 3.
24+
25+
SPOCK_EoS.xlsx
26+
--------------
27+
This Excel spreadsheet contains an implementation of the SPOCK equation
28+
of state.

‎contrib/spock_eos/SPOCK_EoS.xlsx‎

46.2 KB
Binary file not shown.

‎contrib/spock_eos/data/Au.dat‎

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

0 commit comments

Comments
 (0)