The application of the Young–Laplace equation to a solid–liquid interface is considered. Computer simulations show that the pressure inside a solid cluster of hard spheres is smaller than the external pressure of the liquid (both for small and large clusters). This would suggest a negative value for the interfacial free energy. We show that in a Gibbsian description of the thermodynamics of a curved solid–liquid interface in equilibrium, the choice of the thermodynamic (rather than mechanical) pressure is required, as suggested by Tolman for the liquid–gas scenario. With this definition, the interfacial free energy is positive, and the values obtained are in excellent agreement with previous results from nucleation studies. Although, for a curved fluid–fluid interface, there is no distinction between mechanical and thermal pressures (for a sufficiently large inner phase), in the solid–liquid interface, they do not coincide, as hypothesized by Gibbs.

Under certain conditions (i.e., constant number of particles N, volume V, and temperature T), it is possible to have a spherical phase in equilibrium with another phase around it. There, the Helmholtz free energy F is a local/global minimum representing a metastable/stable equilibrium state.1–8 This equilibrium implies that T and chemical potential μ are homogeneous. Thus, ∇T(r) = 0 and ∇μ(r) = 0, where r is the position vector. However, the number density ρ(r) = dN(r)/dV(r) and the pressure tensor p(r) are inhomogeneous.9,10 By taking the center of mass of the cluster (COM) as the origin and using spherical coordinates,

p(r)=pN(r)[erer]+pT(r)[eθeθ+eϕeϕ],
(1)

where er, eθ, and eϕ are the unitary vectors. Then, the condition of mechanical equilibrium, ∇ · p = 0, implies9–11 

pT(r)=pN(r)+r2dpNdr.
(2)

In the late 1970s and 1980s, Rusanov, Rowlinson, and Gubbins and co-workers pioneered the application of computer simulations to study fluid–fluid spherical interfaces at equilibrium.1,11,12 Lately, there has been a revival in the study of these systems,5,6,8,13–27 including also solid–fluid curved interfaces.4,7,28–31 In fact, we have recently shown that the stable equilibrium observed in the isochoric-isothermal (NVT) ensemble is an unstable equilibrium in the isobaric-isothermal (NpT) ensemble that corresponds to a maximum in the Gibbs free energy G. Thus, nucleation can be studied via both stable and unstable equilibrium as they are two sides of the same coin.7,8

The best thermodynamic description of a system with a curved interface in equilibrium can be found in the book of Rowlinson and Widom.9,32 Following in the spirit of Gibbs, one assumes two macroscopic phases that are homogeneous up to the interface and accounts for an additional contribution due to the interface itself. Taking into account that μ is homogeneous,

F=Nμpint43πR3pextV43πR3+4πR2γ,
(3)

where γ is the interfacial free energy, R is the radius of the spherical phase, and pint and pext are the respective internal and external pressures.

At a molecular scale, there is some arbitrariness in determining R. Since F, μ, pint, and pext are fixed, changing R also changes γ. There are two popular choices for R. The first is the Gibbs dividing surface, R = Re, for which the number of excess particles is zero (meaning that particles belong either to the solid or to the liquid, but not to the interfacial region). The second is the surface of tension, R = Rs, for which γ is a minimum (γs).

Actually, by taking the notational derivative (i.e., an arbitrary change in R without any physical change in the system), one obtains9,32

pintpext=2γR+dγdR.
(4)

By definition, dγ/dR=0 when R = Rs leading to the celebrated Young–Laplace equation,

pintpext=2γsRs.
(5)

Since γs is positive, this equation shows that the pressure inside the spherical phase is higher than outside and the difference depends on γs and Rs. Simulation studies of fluid–fluid phases by Gubbins and co-workers11 and Vrabec and co-workers33 have confirmed the higher pressure of the internal phase. However, this equation has not been tested for a solid–fluid curved interface. This is the goal of this letter.

Recently, we simulated several solid clusters in equilibrium with a liquid7 via the pseudo-hard-sphere PHS continuous potential34 (hereafter, simply HS), which allows us to simulate with the standard molecular dynamics package GROMACS.35 

Here, for three selected clusters labeled IV, VII, and VIII in our previous work7 (see Ref. 7 for further details on the size of these solid clusters and the way they were obtained using NVT simulations), we launch new trajectories (in the NVT ensemble) saving configurations very often allowing us to compute the pressure tensor. Since the definition of the pressure tensor is locally arbitrary, we choose to use the Irving–Kirkwood36 convention in which the forces between two particles act in the line connecting them. Further details can be found in the supplementary material (SM). In addition, density profiles are provided.

The simulation details including the interaction potential, GROMACS setup, and order parameter to label particles as solid or liquid are exactly the same as in our previous work.7 We shall use here reduced units. The lengths are given in units of σ (i.e., the hard sphere diameter), the time unit is σm/(kT), the densities are given as ρ = N/3, the pressures are given in units of kT/σ3, the interfacial free energies are given in units of kT/σ2, and the chemical potentials are given in units of kT. The pressure profiles are computed up to half of the simulation box L/2, whereas the density profiles, following Ref. 37, cover the whole system.

The density profile and the normal and tangential components of the pressure tensor are presented in Fig. 1. The values of the densities of the solid and the liquid when they reach a plateau and that of the pressure (far from the interface) are presented in Table I. These are obtained by averaging the data from the corresponding plateaus.

FIG. 1.

Radial density (top) and pressure (bottom) profiles from the COM for clusters IV (left) and VIII (right). For the meaning of psolμ, see the main text. ⟨p⟩ is the average pressure of the system, as obtained from the virial theorem. p¯(r) is the average pressure at a distance r, as given by (2/3)pT(r) + (1/3)pN(r). The solid black line is a fit to pN data, while the red dashed line is obtained from Eq. (2) using the pN fit. In the radial density plot, we show the value of ρsol, which would be the density of a bulk solid at psol, and ρliq, which would be the density of a bulk liquid at pliq. In the pressure profile, we show psol, which would be the pressure of a bulk solid having the density ρsol, and pliq, which would be the pressure of a bulk liquid having the density ρliq. The value labeled as psolnucleation corresponds to psolnucleation=pliq+2γs/Rs when using the value of γs and Rs from nucleation studies.7 

FIG. 1.

Radial density (top) and pressure (bottom) profiles from the COM for clusters IV (left) and VIII (right). For the meaning of psolμ, see the main text. ⟨p⟩ is the average pressure of the system, as obtained from the virial theorem. p¯(r) is the average pressure at a distance r, as given by (2/3)pT(r) + (1/3)pN(r). The solid black line is a fit to pN data, while the red dashed line is obtained from Eq. (2) using the pN fit. In the radial density plot, we show the value of ρsol, which would be the density of a bulk solid at psol, and ρliq, which would be the density of a bulk liquid at pliq. In the pressure profile, we show psol, which would be the pressure of a bulk solid having the density ρsol, and pliq, which would be the pressure of a bulk liquid having the density ρliq. The value labeled as psolnucleation corresponds to psolnucleation=pliq+2γs/Rs when using the value of γs and Rs from nucleation studies.7 

Close modal
TABLE I.

Densities and pressures determined at the respective plateaus in the density and pressure profiles. The difference in pressure is also given as Δp = psolpliq. The notation of IV, VII, and VIII refers to the clusters labeled in this way in Ref. 7.

LabelRsρsolρliqpsolpliqΔp
IV 10.791 1.0613 0.9619 12.6046 12.7437 −0.1391 
VII 15.20 1.0548 0.9560 12.3053 12.4047 −0.0994 
VIII 17.467 1.0529 0.9541 12.2199 12.3003 −0.0804 
LabelRsρsolρliqpsolpliqΔp
IV 10.791 1.0613 0.9619 12.6046 12.7437 −0.1391 
VII 15.20 1.0548 0.9560 12.3053 12.4047 −0.0994 
VIII 17.467 1.0529 0.9541 12.2199 12.3003 −0.0804 

Close to the interface, the normal and tangential components of the pressure tensor are different, albeit, far from it, both are identical. Surprisingly, the pressure inside (solid) is smaller than outside (liquid). This result, in principle, contradicts the Young–Laplace equation. Notice though that having a lower pressure for the solid phase does not violate the mechanical equilibrium condition, which only requires a certain relation between pN and pT [Eq. (2)].

This is opposite to the fluid–fluid curved interface. Actually, all previous studies on curved interfaces with fluid phases found higher pressure for the internal phase.1,11,12 Here, for a solid spherical cluster, we found lower pressure in the internal phase. One may think that this behavior is peculiar for HS, for which there are no attractive forces. However, recently, Gunawardana and Song31 have reported a similar behavior for a solid cluster of Lennard-Jones particles surrounded by liquid.

Interestingly, one can already learn the behavior of the curved interface by analyzing the behavior of the pressure tensor for the planar interface. It turns out that, in the interfacial region of a planar interface, pT < pN for fluid–fluid interfaces38 and pT > pN for solid–fluid interfaces.39 Thus, a simple analysis of the behavior of the pressure tensor for the planar interface is sufficient to know if one will have higher or lower pressure in the internal spherical phase [see Eq. (2)]. It is also interesting to point out that the pressure of the external phase, pext, is identical to the average pressure of the system, ⟨p⟩, obtained from the virial equation applied to the entire system provided that the normal and tangential components are identical at L/2, as demonstrated in the supplementary material.

Although one must accept the fact that the pressure inside a solid cluster is smaller than the pressure outside, the consequence of this appears to be dramatic as this would imply (apparently) from the Young–Laplace equation that γ is negative.40 The Young–Laplace equation is explained in any textbook of physics, and now, we have a problem about how to use it in the case of a solid cluster surrounded by liquid. How to reconcile the results of this work with the Young–Laplace equation? The key was provided by Tolman in his celebrated paper discussing the variation of γ with R in a droplet. In particular, there is a remark by Tolman,41 which we believe is highly important in this context. The remark is as follows (adapting his notation to this paper):

In applying Eqs. (2.2) and (2.3) to very small droplets, it is to be noted that pintand ρintare to be taken as the pressure and density for a large mass of internal phase in a condition at the temperature of interest to give the same value of μ as that of the vapor (cf. Gibbs, Ref. 1 , p. 253).

Notice that Tolman was describing the equilibrium between a droplet of liquid and its vapor. However, it also applies for the solid–liquid interface as we are about to show. Therefore, when using the formalism of Gibbs,42 the pressure of the internal spherical phase should not be taken as its actual value but rather from that of a bulk having the same μ as the external phase. Similar reasoning was also used by ten Wolde and Frenkel.43 

Determining the exact value of μ of inhomogeneous systems of high density is very difficult44 (notice though some recent progress45). Thus, we do not know the exact value of μ for the three systems considered in this work. However, to illustrate our main point, this is not crucial. We shall assume that the external liquid has bulk behavior so that μ in the system corresponds to that of a bulk liquid at pliq. In the inset of Fig. 2(a), μ(p) for solid and liquid bulks are presented (obtained via thermodynamic integration46 from p = 11.648 that is the coexistence pressure47 where μ is identical in both phases).

FIG. 2.

(a) Difference in pressure between the solid cluster and the liquid as a function of 1/Rs: upper curve—psolμpliq; lower curve—psolpliq. Note that there are two sets for the mechanical Δp. The solid black circles are estimated (for all clusters in Ref. 7) by using the respective densities and the equation of state (EOS) rather than by performing a more costly pressure tensor calculation, as was done for the empty circles. The former is systematically smaller than the latter, suggesting that the cluster differs slightly from a bulk. In the inset of this panel, the chemical potential of both phases is shown. (b) Snapshots of the solid cluster VIII in the initial configuration and after some time. Only solid particles within the system are shown. We followed the ten closest particles to the COM and their first coordination shells (cyan and red spheres). The remaining solid particles are shown as blue dots. In red, an example of solid particle that changed neighbor. (c) Maximum in the probability distribution function of the closest 12th neighbor considering cluster VIII and solid bulks with and without vacancies at psolVIII. The ratios 1/1000 and 1/2000 mean the proportion of vacancy per number of particles for the considered bulks. Only the particles at COM < r < 10σ were considered.

FIG. 2.

(a) Difference in pressure between the solid cluster and the liquid as a function of 1/Rs: upper curve—psolμpliq; lower curve—psolpliq. Note that there are two sets for the mechanical Δp. The solid black circles are estimated (for all clusters in Ref. 7) by using the respective densities and the equation of state (EOS) rather than by performing a more costly pressure tensor calculation, as was done for the empty circles. The former is systematically smaller than the latter, suggesting that the cluster differs slightly from a bulk. In the inset of this panel, the chemical potential of both phases is shown. (b) Snapshots of the solid cluster VIII in the initial configuration and after some time. Only solid particles within the system are shown. We followed the ten closest particles to the COM and their first coordination shells (cyan and red spheres). The remaining solid particles are shown as blue dots. In red, an example of solid particle that changed neighbor. (c) Maximum in the probability distribution function of the closest 12th neighbor considering cluster VIII and solid bulks with and without vacancies at psolVIII. The ratios 1/1000 and 1/2000 mean the proportion of vacancy per number of particles for the considered bulks. Only the particles at COM < r < 10σ were considered.

Close modal

Taking, for instance, system VIII where pliq = 12.3003 and psol = 12.2199, it is possible to determine μ for a bulk liquid phase at this pressure and also the pressure of a solid that has the same value of μ, which we found to be psolμ=12.37. The superscript μ reminds us that this pressure is not the mechanical pressure of the solid but rather the pressure of a bulk solid that has the same μ as that of a bulk liquid at pliq. We shall denote this as the thermodynamical pressure (as opposed to the mechanical pressure psol). Notice that psolμ>pliq>psol. This finding suggests that the cluster must be different from a perfect bulk at the same pressure psol; otherwise, μ could not be homogeneous and no equilibrium could be reached. In order to find differences, we followed the evolution of the closest particles to the COM finding that the solid cluster is a “living” structure that can melt in a certain region and grow in another while keeping the size approximately constant. As can be seen in Fig. 2(b), the selected particles ended up quite close to the interface and some of them changed their neighbor. This is likely due to the presence of vacancies in the cluster that lead to a relative diffusion. By computing such diffusion for the cluster as well as for solid bulks with and without vacancies at psolVIII, we could estimate the cluster to have about one vacancy per four thousand particles (1/4000). More relevant is the distribution of distances from a given particle to its 12th closest neighbor considering only particles that fulfill COM < r < 10σ in order to avoid surface effects (setting the upper limit in 7σ did not produce any difference). As shown in Fig. 2(c), such distribution is shifted to higher distances with respect to a bulk with and without vacancies indicating a tiny expansion in the cluster lattice. Further work is needed to completely understand this. Nevertheless, it is clear that the cluster is not identical to a bulk solid with or without vacancies at the same mechanical pressure.

Therefore, for the solid–liquid interface, the Young–Laplace equation must be written as

psolμpliq=2γsRs.
(6)

In Table II, the value of psolμ for the three systems considered in this work is presented along with the difference in pressure Δpμ=psolμpliq.

TABLE II.

Thermodynamic pressure of the solid and the difference with the pressure of the liquid phase (from Table I). By using the values of γs from nucleation studies, we estimated the term 2γs/Rs and found it to be in excellent agreement with the difference in pressure obtained when using pμ.

LabelRspsolμΔpμ2γs/Rs7 
IV 10.791 12.8627 0.1190 0.1164 
VII 15.20 12.4846 0.0799 0.0793 
VIII 17.467 12.3700 0.0697 0.0694 
LabelRspsolμΔpμ2γs/Rs7 
IV 10.791 12.8627 0.1190 0.1164 
VII 15.20 12.4846 0.0799 0.0793 
VIII 17.467 12.3700 0.0697 0.0694 

As can be seen, Δpμ is positive. Thus, by using Tolman’s suggestion, one recovers a “normal” Young–Laplace equation. In previous work where the same clusters (IV, VII, and VIII) were studied, we obtained the values of γs from nucleation studies. Since, according to Eq. (6), Δpμ corresponds to 2γs/Rs, it is of interest to analyze whether our previously reported values of γs and Rs are consistent with this difference of pressures. As shown in Table II, results are fully consistent. Thus, the physical meaning of 2γs/Rs obtained for the values of γs and Rs from nucleation studies7,48 is now clear. Note that, for a fluid–fluid interface, there is no difference between pintμ and pint (for a sufficiently large inner phase), whereas, in a solid–liquid system, we could not find such agreement even for very large clusters. We have plotted the difference in pressure between the solid and liquid as a function of 1/Rs in Fig. 2(a). As can be seen, there is no evidence that this difference can become positive for a certain value of Rs.

The idea that, for the solid–liquid interface, the difference in pressure between the phases may not lead to γ was already insinuated by Gibbs. Later, Cahn,49 Cammarata,50,51 and others52 suggested that the strain, which is present in solids and not in fluids, was behind this.

In summary, we have computed the pressure tensor for a HS system at constant N, V, and T where one has a stable solid cluster in contact with a liquid away from coexistence conditions. We found that the internal pressure (solid) is lower than the external one (liquid). This would lead to a negative γ. However, as suggested by Tolman (and insinuated by Gibbs), defining a thermal pressure for the inner phase, which corresponds to that of a solid with the same chemical potential as the external liquid, allows us to recover a normal Young–Laplace equation, where the pressure of the internal phase is higher, leading to a positive γ. The values of γ from this scheme are in excellent agreement with recent results from nucleation studies. Thus, for a solid–liquid interface, one should distinguish between the mechanical and the thermodynamic pressure. This distinction is not so necessary for a fluid–fluid curved interface as they are comparable.25,51 However, it is crucial in understanding the meaning of the Young–Laplace equation for a solid–fluid interface. Computer simulations have been the key to solve this subtle issue.

See the supplementary material for a description of the system and its interaction potential, details on the pressure tensor calculation, values for the fitting parameters, and the demonstration that the average pressure in the system equals the external pressure. Additional figures and information can also be found.

This work was funded by Grant Nos. FIS2016-78117-P and PID2019-105898GB-C21 and FIS2017-89361-C3-2-P of the MEC, by Project No. UCM-GR17-910570 from UCM, and by the USA National Science Foundation under Award No. CBET-1855465. P.M.d.H. acknowledges financial support from the FPI Grant No. BES-2017-080074.

The data that support the findings of this work are available within the article and its supplementary material.

1.
D. J.
Lee
,
M. M.
Telo da Gama
, and
K. E.
Gubbins
,
J. Chem. Phys.
85
,
490
(
1986
).
2.
A. I.
Rusanov
,
Russ. Chem. Rev.
33
,
385
(
1964
).
3.
A. K.
Shchekin
,
K.
Koga
, and
N. A.
Volkov
,
J. Chem. Phys.
151
,
244903
(
2019
).
4.
P.
Koß
,
A.
Statt
,
P.
Virnau
, and
K.
Binder
,
Mol. Phys.
116
,
2977
(
2018
).
5.
V. G.
Baidakov
,
Chem. Phys.
525
,
110407
(
2019
).
6.
V.
Baidakov
and
S.
Protsenko
,
Doklady Physics
(
Springer
,
2004
), Vol. 49, pp.
69
72
.
7.
P.
Montero de Hijes
,
J. R.
Espinosa
,
V.
Bianco
,
E.
Sanz
, and
C.
Vega
,
J. Phys. Chem. C
124
,
8795
(
2020
).
8.
P.
Rosales-Pelaez
,
I.
Sanchez-Burgos
,
C.
Valeriani
,
C.
Vega
, and
E.
Sanz
,
Phys. Rev. E
101
,
022611
(
2020
).
9.
J. S.
Rowlinson
and
B.
Widom
,
Molecular Theory of Capillarity
(
Oxford University Press
,
1982
).
10.
J. S.
Rowlinson
,
Pure Appl. Chem.
65
,
873
(
1993
).
11.
S. M.
Thompson
,
K. E.
Gubbins
,
J. P. R. B.
Walton
,
R. A. R.
Chantry
, and
J. S.
Rowlinson
,
J. Chem. Phys.
81
,
530
(
1984
).
12.
A. I.
Rusanov
and
E. N.
Brodskaya
,
J. Colloid Interface Sci.
62
,
542
(
1977
).
13.
A.
Tröster
,
M.
Oettel
,
B.
Block
,
P.
Virnau
, and
K.
Binder
,
J. Chem. Phys.
136
,
064709
(
2012
).
14.
B. J.
Block
,
S. K.
Das
,
M.
Oettel
,
P.
Virnau
, and
K.
Binder
,
J. Chem. Phys.
133
,
154702
(
2010
).
15.
M.
Schrader
,
P.
Virnau
, and
K.
Binder
,
Phys. Rev. E
79
,
061104
(
2009
).
16.
M.
Schrader
,
P.
Virnau
,
D.
Winter
,
T.
Zykova-Timan
, and
K.
Binder
,
Eur. Phys. J.: Spec. Top.
177
,
103
(
2009
).
17.
K.
Binder
,
B. J.
Block
,
P.
Virnau
, and
A.
Tröster
,
Am. J. Phys.
80
,
1099
(
2012
).
18.
A.
Tröster
and
K.
Binder
,
Phys. Rev. Lett.
107
,
265701
(
2011
).
19.
P.
Virnau
,
F.
Schmitz
, and
K.
Binder
,
Mol. Simul.
42
,
549
(
2016
).
20.
A.
Tröster
,
F.
Schmitz
,
P.
Virnau
, and
K.
Binder
,
J. Phys. Chem. B
122
,
3407
(
2017
).
21.
L. G.
MacDowell
,
P.
Virnau
,
M.
Müller
, and
K.
Binder
,
J. Chem. Phys.
120
,
5293
(
2004
).
22.
L. G.
MacDowell
,
V. K.
Shen
, and
J. R.
Errington
,
J. Chem. Phys.
125
,
034705
(
2006
).
23.
Ø.
Wilhelmsen
,
D.
Bedeaux
,
S.
Kjelstrup
, and
D.
Reguera
,
J. Chem. Phys.
140
,
024704
(
2014
).
24.
J. G.
Sampayo
,
A.
Malijevský
,
E. A.
Müller
,
E.
de Miguel
, and
G.
Jackson
,
J. Chem. Phys.
132
,
141101
(
2010
).
25.
A.
Malijevský
and
G.
Jackson
,
J. Phys.: Condens. Matter
24
,
464121
(
2012
).
26.
G. V.
Lau
,
I. J.
Ford
,
P. A.
Hunt
,
E. A.
Müller
, and
G.
Jackson
,
J. Chem. Phys.
142
,
114701
(
2015
).
27.
S. H.
Min
and
M. L.
Berkowitz
,
J. Chem. Phys.
150
,
054501
(
2019
).
28.
A.
Statt
,
P.
Virnau
, and
K.
Binder
,
Phys. Rev. Lett.
114
,
026101
(
2015
).
29.
P.
Koß
,
A.
Statt
,
P.
Virnau
, and
K.
Binder
,
Phys. Rev. E
96
,
042609
(
2017
).
30.
D.
Richard
and
T.
Speck
,
J. Chem. Phys.
148
,
224102
(
2018
).
31.
K. G. S. H.
Gunawardana
and
X.
Song
,
J. Chem. Phys.
148
,
204506
(
2018
).
32.
S.
Kondo
,
J. Chem. Phys.
25
,
662
(
1956
).
33.
J.
Vrabec
,
G. K.
Kedia
,
G.
Fuchs
, and
H.
Hasse
,
Mol. Phys.
104
,
1509
(
2006
).
34.
J.
Jover
,
A. J.
Haslam
,
A.
Galindo
,
G.
Jackson
, and
E. A.
Müller
,
J. Chem. Phys.
137
,
144505
(
2012
).
35.
B.
Hess
,
C.
Kutzner
,
D.
van der Spoel
, and
E.
Lindahl
,
J. Chem. Theory Comput.
4
,
435
(
2008
).
36.
J. H.
Irving
and
J. G.
Kirkwood
,
J. Chem. Phys.
18
,
817
(
1950
).
37.
M.
Deserno
, “
How to calculate a three-dimensional g(r) under periodic boundary conditions
,” www.cmu.edu/biolphys/deserno/pdf/gr_periodic.pdf,
2004
.
38.
A.
Trokhymchuk
and
J.
Alejandre
,
J. Chem. Phys.
111
,
8510
(
1999
).
39.
R. L.
Davidchack
and
B. B.
Laird
,
J. Chem. Phys.
108
,
9452
(
1998
).
40.
A. O.
Tipeev
,
J. Chem. Phys.
151
,
017101
(
2019
).
41.
R. C.
Tolman
,
J. Chem. Phys.
17
,
333
(
1949
).
42.
J.
Gibbs
,
The Collected Works
, Vol. 1. Thermodynamics (
Yale University Press
,
1948
).
43.
P. R.
ten Wolde
and
D.
Frenkel
,
J. Chem. Phys.
109
,
9901
(
1998
).
44.
J. G.
Powles
,
S. E.
Baker
, and
W. A. B.
Evans
,
J. Chem. Phys.
101
,
4098
(
1994
).
45.
C.
Perego
,
O.
Valsson
, and
M.
Parrinello
,
J. Chem. Phys.
149
,
072305
(
2018
).
46.
C.
Vega
,
E.
Sanz
,
J. L. F.
Abascal
, and
E. G.
Noya
,
J. Phys.: Condens. Matter
20
,
153101
(
2008
).
47.
J. R.
Espinosa
,
E.
Sanz
,
C.
Valeriani
, and
C.
Vega
,
J. Chem. Phys.
139
,
144502
(
2013
).
48.
P.
Montero de Hijes
,
J. R.
Espinosa
,
E.
Sanz
, and
C.
Vega
,
J. Chem. Phys.
151
,
144501
(
2019
).
50.
R. C.
Cammarata
and
K.
Sieradzki
,
Annu. Rev. Mater. Sci.
24
,
215
(
1994
).
51.
R. C.
Cammarata
,
Solid State Phys.
61
,
1
(
2009
).
52.
B. B.
Laird
,
R. L.
Davidchack
,
Y.
Yang
, and
M.
Asta
,
J. Chem. Phys.
131
,
114110
(
2009
).

Supplementary Material