of 11
Supplementary materials
1
0.1. The Matsui potential for silicate systems
2
The Matsui (1998) potential model in its original form can be represented by the following
3
equation:
4
V
(
r
ij
) =
q
i
q
j
r
ij
C
j
C
j
r
6
ij
+
f
(
B
i
+
B
j
) exp
(
A
i
+
A
j
r
ij
B
i
+
B
j
)
(1)
From left to right, the three terms describes Coulomb interaction, van der Waals attraction
5
and Born-Mayer-Huggins repulsion, respectively. Born-Mayer-Huggins repulsion is a soft-
6
sphere model for closed-shell atoms and ions published by Gilbert (1968). The physical
7
meaning of
f
in this term is an arbitrary standard force by which the atom pairs are pushed
8
apart until the distance between them is the sum of their soft-core radius. Matsui set his
9
standard force to be 4.184 kJ Å
1
mol
1
. The soft-core radius of atom
i
and
j
are denoted
10
by
A
i
and
A
j
, respectively.
B
i
and
B
j
are called softness (or hardness) parameters.
11
To simplify the form of this pairwise potential, we can perform the following transformations
12
to equation 1.
13
C
ij
=
C
i
C
j
(2)
B
ij
=
B
i
+
B
j
(3)
A
ij
=
fB
ij
exp
(
A
i
+
A
j
B
ij
)
(4)
And we can get a new form as follows.
14
V
(
r
ij
) =
q
i
q
j
e
2
4
π
0
r
ij
+
A
ij
exp
(
r
ij
B
ij
)
C
ij
r
6
ij
(5)
In this new form,
A
ij
and
C
ij
are called energy parameters for the pair
ij
and they describe
15
the repulsive and van der Waals attractive forces, respectively.
B
ij
is an e-folding length
16
characterizing the radially symmetric decay of electron repulsion energy between atom pair
17
ij
. The second and third terms together are known as a Buckingham potential.
18
In LAMMPS, this pairwise potential can be described by the built-in buck/coul/long pair
19
style.
20
E
Buck
=
Ae
r/η
C
r
6
, r < r
c
(6)
E
Coul
=
Cq
i
q
j
r
, r < r
c
(7)
E
total
=
E
Buck
+
E
Coul
(8)
In our force field optimizations, we started from the parameters given by Matsui (1998) and
21
modified the van der Waals part (
E
Buck
) only.
22
In all simulations, we used PPPM long-range interactions for the Coulomb potential and set
23
the cut-off radius for the Buckingham potential to be
r
c
=
11 Å.
24
1
0.2. Transformation of Buckingham potential
25
To understand the physical meanings of the optimization parameters and make more
26
precise adjustments, we transformed the van der Waals part of the pairwise potential
27
according to the exponential-6 (X6) form (Mayo, Olafson and Goddard, 1990).
28
29
E
X6
vdw
=
D
0
[
6
ζ
6
e
ζ
(1
ρ
)
ζ
ζ
6
ρ
6
]
(9)
where
D
0
is the van der Waals well depth,
ρ
=
R/R
0
is relative bond length where
R
0
is
30
equilibrium bond length and
ζ
is a dimensionless scaling parameter.
ζ
= 12
.
0
means the X6
31
form and the LJ form have the same long range interactions and
ζ
= 13
.
772
means they have
32
the same force constant. The X6 form has advantages in describing short range interactions
33
and is suitable for our force field where covalent and ionic bonds are dominant.
34
By comparing Eqn. 6 and Eqn. 9, we can find a set of solutions that represents
D
0
,
R
0
and
35
ζ
using
A
,
η
and
C
. The transformation of force field parameters is automated by Python
36
and done in Mathematica (Wolfram Research, Inc., 2020).
37
0.3. Optimization details
38
The parameters that are tuned in our optimization are the transformed Buckingham pa-
39
rameters
D
0
,
R
0
and
ζ
for all Cat-O pairs (tuned simutaneously) and for O-O pairs. There
40
are six parameters in total.
41
To avoid cancelling of different parts, we used the sum of log (e-based) losses as our opti-
42
mization goal.
43
L
= log(
α
P
L
P
+
α
C
V
L
C
V
+
α
γ
L
γ
)
(10)
where
L
is the total loss,
α
’s are the scaling parameters for each kind of loss and
L
P
,
L
C
V
44
and
L
γ
are the losses corresponding to pressure, heat capacity and Grüneisen parameter,
45
respectively. For a simulation containing 4 reference points (which yields 4 pressures, 2
46
heat capacities and 2 Grüneisen parameters because the latter two are calculated by finite
47
difference), each part of the loss function is defined as:
48
L
P
=
1
4
4
i
=1
(
P
i
P
i,ref
)
2
(11)
L
C
V
=
1
2
2
i
=1
(
C
V i
C
V i,ref
)
2
(12)
L
γ
=
1
2
2
i
=1
(
γ
i
γ
i,ref
)
2
(13)
2
The choice of scaling parameters considers both the relative magnitude of each physical
49
quantity and their degree of fluctuation. For our final optimizations, the choices are:
50
α
P
= 0
.
01
(14)
α
C
V
= 100
(15)
α
γ
= 70000000
(16)
Both BFGS and gradient descent (GD) algorithms were used in the optimization and al-
51
though BFGS reduced the loss function faster at the beginning of the optimization, the best
52
result was produced by GD, which reduced the loss function from 3.90 to 3.59 and achieved
53
satisfying accuracy for
C
V
and
γ
.
54
0.4. The force field parameters
55
Compared to the original Matsui force field, the optimized force field scales each param-
56
eter differently (Table S1).
57
The
van
der
Waals
well
depth
parameter
(
D
0
)
has
the
largest
change
of
all,
which
can
Table
S
1:
Final
scaling
factors
of
the
optimized
force
field
Optimization paramter
Scaling=new/original
D
0
(Cat-O)
0
.
478192
R
0
(Cat-O)
0
.
771891
ζ
(Cat-O)
1
.
044594
D
0
(O-O)
0
.
684991
R
0
(O-O)
0
.
990607
ζ
(O-O)
0
.
804164
58
be attributed to the fact that high temperature and pressure conditions reduce the bond
59
strength. For all of the oxygen-containing ion pairs, the change of equilibrium bond length
60
(
R
0
) is very small, indicating that the structure of silicates can be relatively accurately
61
described by the original force field. For Cation-Oxygen and Oxygen-Oxygen pairs, the cur-
62
vature parameter (
ζ
) changes in opposite directions, which could be related to the difference
63
in long range attraction between cation-anion pairs and anion-anion pairs.
64
These parameters are then transformed to input parameters for LAMMPS and can be found
65
in section 0.9.
66
0.5. Calculating Hugoniot from Equation of State
67
The original Hugoniot data shown in Figure 5 were presented by Asimow and Ahrens
68
(2010). The measured quantities in each experiment are shock velocity
U
S
, particle velocity
69
u
p
, and initial sample density
ρ
o
. These results are converted to the Pressure (
P
)-Density
70
(
ρ
) plane using the first two Rankine-Hugoniot equations (i.e., conservation of mass and
71
3
momentum) for a single shock driven into a sample initially at zero pressure and at rest in
72
the laboratory reference frame:
73
ρ
o
U
S
=
ρ
(
U
S
u
p
)
(17)
P
=
ρ
o
U
S
u
p
.
(18)
In most of the experiments plotted, the initial state of the sample is liquid at 0.1 MPa (i.e.,
74
approximately zero) pressure at 1673 K, which is also the reference state for the equation of
75
state adopted. One experiment begins from the initially solid state at 300 K. Interpreting
76
this experiment requires knowing the initial density of the isochemical mixture of solids at
77
300 K (
ρ
OS
), as well the internal energy difference between the 1673 K liquid and 300 K
78
solid states,
E
tr
, which is the obtained from the integral of the solid heat capacity from
79
300 K to the eutectic temperature (1547 K), the enthalpy of fusion, and the integral of the
80
liquid heat capacity from 1547 K to 1673 K.
81
These data are compared to Hugoniots computed from the thermal equation of state as
82
follows. The independent variable is
ρ
. At each value of
ρ
, we obtain the Eulerian finite
83
strain parameter
f
relative to the reference volume
ρ
o
:
84
f
=
{
(
ρ
ρ
o
)
2
3
1
}
.
(19)
and then the pressure on the reference isentrope is calculated from the
3
rd
order Birch-
85
Murnaghan equation of state:
86
P
S
= 3
K
oS
f
(2
f
+ 1)
5
2
[
1 +
3
2
(
K
S
4)
f
]
,
(20)
In which
K
oS
is the isentropic bulk modulus at the reference state (which is known precisely
87
from ultrasonic measurements and therefore is held constant during fitting) and
K
S
is the
88
first derivative of the isentropic bulk modulus along the isentrope, evaluated at the reference
89
state (which is not known from ambient pressure data and is therefore a fitted parameter).
90
To obtain the pressure on the Hugoniot at shock density
ρ
, we take advantage of the fact
91
that internal energy
E
is a path-independent state variable and therefore the difference in
92
E
between the initial reference state and the shock state is the same whether the change of
93
state is accomplished by passage of the shock (
E
H
) or by decomposition of the path into
94
a series of three processes: transformation from the initial state to the reference state at 1
95
bar (
E
tr
, if necessary), isentropic compression from the reference state at
ρ
o
to the shock
96
density
ρ
(
E
S
) and finally isochoric heating from the isentrope to the corresponding state
97
on the Hugoniot (
E
V
). That is:
98
E
H
= ∆
E
tr
+ ∆
E
S
+ ∆
E
V
.
(21)
The increase in internal energy imparted by passage of a shock is obtained from the
3
rd
99
Rankine-Hugoniot jump condition:
100
E
H
=
1
2
(
P
+
P
o
)
(
1
ρ
oo
1
ρ
)
,
(22)
4
in which ambient pressure
P
o
is negligible compared to shock pressure
P
and
ρ
oo
is the
101
density of the initial state, which may be equal to
ρ
o
if the sample begins in the reference
102
state,
ρ
os
for the case where the initial state is a cold solid assemblage, or any other value
103
as appropriate (e.g. for a porous target).
104
From the first and second laws of thermodynamics, the total derivative of specific internal
105
energy d
E
=
T
d
S
P
d
V
, where
T
is temperature,
S
is specific entropy, and
V
is specific
106
volume
(
V
= 1
)
. Hence, the increase in internal energy along an isentropic path (for
107
which d
S
= 0
) is obtained by integrating equation (20):
108
E
S
=
9
2
K
oS
ρ
o
[
f
2
+ (
K
S
4)
f
3
]
.
(23)
The increase in energy along the isochore from the isentrope to the Hugoniot is obtained
109
from the selected formulation of the Grüneisen parameter, since
γ
V
(
∂P/∂E
)
V
. For the
110
conventional Mie-Grüneisen approximation, in which
γ
is constant along an isochore, this is
111
a simple finite difference:
112
E
V
=
V
γ
(
V
)
(
P
P
S
)
.
(24)
In this work, when using the Mie-Grüneisen approximation we have assumed that
γ
(
V
)
is
113
given by
114
γ
(
V
)
γ
o
=
(
V
V
o
)
q
,
(25)
In which the Grüneisen parameter at the reference state
γ
o
is known precisely from ultra-
115
sonic, thermal expansion, and heat capacity data and is a fixed parameter during fitting,
116
whereas
q
is unknown and is a fitted parameter.
117
The solution of these equations for Hugoniot pressure at a given density, for the Mie-
118
Grüneisen case, can be compactly written
119
P
H
=
P
S
ργ
(∆
E
S
+ ∆
E
tr
)
1
γ
2
(
ρ
ρ
oo
1
)
,
(26)
If the sample begins in the reference state,
E
tr
= 0
. Hence equation (26) is used to
120
compute both the hot liquid Hugoniot and the cold solid Hugoniot curves for the Mie-
121
Grüneisen model.
122
Note that an alternative expression for
γ
is
γ
= (
ln
T/∂
ln
V
)
S
. Hence the temperature
123
along an isentrope can be obtained by integrating this expression. For the Mie-Grüneisen
124
case and the form of equation (23), we have
125
T
S
=
T
o
exp
{
γ
o
q
(
1
(
ρ
o
ρ
)
q
)}
.
(27)
5
For the model introduced in this paper,
γ
(
E,V
)
is not a constant along the isochore from
126
the isentrope to the Hugoniot and so we cannot use equation (24). In this model, we have
127
adopted the form
128
γ
(
E,V
) =
(
a
V
+
b
)
(
E
E
) +
γ
,
(28)
where
a
,
b
,
E
, and
γ
are the adjustable parameters. Note that the reference datum
129
γ
(
E
o
,V
o
) =
γ
o
(29)
is added to the fit and assigned an extremely high weight, in order to enforce the known
130
reference value of the Grüneisen parameter.
131
The relationship between pressure increase and internal energy increase along the isochore
132
from the isentrope to the Hugoniot state is given in this case by
133
P
P
S
=
1
V
E
H
E
S
+∆
E
tr
γ
(
E
,V
)
d
E
.
(30)
In combination with equations (21) and (22), integrating the linear dependence on
E
from
134
(28) yields a quadratic expression for
P
, of which one root is physical and one root is non-
135
physical. With appropriate choice of
ρ
oo
and
E
tr
, this quadratic solution yields either the
136
Hugoniot from the hot liquid initial state or from the cold solid initial state, for the
γ
(
E,V
)
137
model.
138
Since the changes in both
E
and
P
with changes in
V
along the isentrope are known from
139
equations (20) and (23), it is straightforward to obtain
γ
(
E,V
)
along the isentrope and
140
hence to compute the temperature increase along the isentrope in this model as well.
141
The final equation that is relevant to fitting the model form to the experimental dataset for
142
AnDi liquid incorporates the measurement of bulk sound velocity
C
B
in the shock state. As
143
with any bulk sound wave,
C
B
=
V K
S
, where again
K
S
is the isentropic bulk modulus.
144
If the slope of the Hugoniot
(
∂P/∂V
)
H
is obtained from the above equations by finite
145
difference, then it can be shown that
146
K
S
=
1
2
(
γP
(
∂P
∂V
)
H
[
2
V
γ
(
1
ρ
oo
V
)])
.
(31)
The inversion for fitted parameters
{
K
S
,a,b,E
}
is therefore based on
γ
o
, all
147
twelve available preheated liquid Hugoniot shots, one cold solid aggregate Hugo-
148
niot shot, and one sound speed datum. All the constants
{
ρ
o
os
,K
oS
o
,
E
tr
}
149
are given by Asimow and Ahrens (2010) but are repeated here for convenience:
150
{
2618
kg m
–3
,
3078
kg m
–3
,
22
.
98
GPa
,
0
.
356
,
1
.
954106
MJ kg
–1
}
.
151
0.6. Ergodicity check by mean squared displacement (MSD)
152
In liquid phase, particles follow the laws of Brownian mothion and their MSD is linear
153
in time (namely, MSD
= 6
Dt
where
D
is the diffusion coefficient.) However, in glass phase,
154
6