Zhi
Li
*a,
Christophe
Winisdoerffer
b,
François
Soubiran
ac and
Razvan
Caracas
ad
aCNRS, École Normale Supérieure de Lyon, Laboratoire de Géologie de Lyon LGLTPE UMR 5276, 46 allée d'Italie, Lyon 69364, France. E-mail: zhi.li@ens-lyon.fr
bUniversité de Lyon, École Normale Supérieure de Lyon, Centre de Recherche Astrophysique de Lyon UMR 5574, 46 allée d'Italie, Lyon 69364, France
cCEA, DAM, DIF, 91297 Arpajon, France
dThe Centre for Earth Evolution and Dynamics (CEED), University of Oslo, Postbox 1028 Blindern, N-0315 Oslo, Norway
First published on 21st December 2020
The ab initio (ai) Gibbs ensemble (GE) Monte Carlo (MC) method coupled with Kohn–Sham density functional theory is successful in predicting the liquid–vapour equilibrium of insulating systems. Here we show that the aiGEMC method can be used to study also metallic systems, where the excited electronic states play an important role and cannot be neglected. For this we include the electronic free energy in the formulation of the effective energy of the system to be used in the acceptance criteria for the MC moves. The application of this aiGEMC method to sodium yields a good agreement with available experimental data on the liquid–vapour equilibrium densities. We predict a critical point for sodium at 2338 ± 108 K and 0.24 ± 0.03 g cm−3. The liquid structure stemming from aiGEMC simulations is very similar to the one from ab initio molecular dynamics. Since this method can determine phase transition without computing the Gibbs free energy, it may offer a new possibility to study other materials with a reasonable computational cost.
Given the available experimental apparatus, the equilibria between condensed phases and gases are much more difficult to investigate than the transitions in condensed matters. And difficulties are only increasing when addressing the critical points, i.e. the points where differences between two phases, like liquids and gases, vanish. Critical points have been extensively studied experimentally especially under near-ambient conditions.7 Under higher pressure, multi-anvil press and diamond anvil cell have been used for decades.8 Unfortunately, critical points lying above 1500 K and at low densities are extremely challenging to obtain experimentally. For these particular conditions, computational studies can be extremely useful, and sometimes are the only way forward for the time being.
The numerical study of phase equilibria is a relatively straightforward problem since it means equating the pressures, temperatures and chemical potentials of two different systems.9 If the first two are relatively easy to achieve numerically, the difficulty comes from the last equality: computing chemical potentials is extremely difficult since it requires the knowledge of the entropy of each phase, which can only be derived from the absolute partition function. Techniques such as molecular dynamics (MD) or Monte Carlo (MC) simulations only provide average quantities but not the absolute partition function.10 Numerous methods have been proposed to circumvent this issue. The most evident one consists in computing the free energy of the system using a reference point and the thermodynamic integration (TDI).11 This is however computationally expensive and a reference system is not always easy to find. A more direct approach consists in simulating layers of two different phases and observe the evolution of the interfaces to determine the most stable phase.12,13 This is more efficient than the TDI but is challenging to set up and highly prone to finite-size effects. Finally, in the Gibbs ensemble Monte Carlo (GEMC) method, one simulates two nearly independent phases at the same time through a MC simulation and ensures the equilibrium by performing specific Metropolis moves allowing to exchange volume and particles.14 This directly mimics the macroscopic equilibrium at a microscopic level, but avoids the interface issues allowing for a more reliable determination of the phase equilibrium.
The success of all the methods outlined in the previous paragraph rely on the quality of the description of the interatomic interaction. In classical simulations these interactions are reproduced by simplified functions that disregard the actual complexity of the electronic interactions. A viable alternative is to fully consider the quantum effects, which are especially important for molecular systems or for metals where the electronic states and occupations change with temperature and pressure. The state-of-the-art method to tackle this issue is the Kohn–Sham formulation (KS) of density functional theory (DFT),15,16 and its finite temperature extension (FT) by Mermin,17 which can account for the quantum behavior of the electrons in a computationally cost-effective manner. TDI18–22 and two-phase simulations23,24 have been linked to DFT calculations to significantly improve the accuracy of the predictions of phase transitions for condensed matter.
In the classical Gibbs ensemble scheme,25 the potential energy evaluated from a force field is used in the acceptance criterion for a trial move. This can be extended to include ab initio simulations, where the internal energy of the electron–ion system plays a similar role as the potential energy in the classical scheme. In the pioneering work of Siepmann et al.,26 the GEMC method has been combined with KS-DFT – often referred as ab initio (aiGEMC) – to improve its accuracy. Ever since different groups have applied the aiGEMC technique to calculate the liquid–vapor equilibrium line of water,27–29 methanol, methane30 and argon.31 All above materials are insulators and have a critical point that is below 700 K. Due to the relatively large band gaps, the thermal excitation of electrons in these systems is not significant, while it may play an important role in determining the phase equilibrium of metallic systems at high temperature.
In metals, in order to capture the thermal excitations of electrons, it is necessary to make full use of FT-KS-DFT which explicitly includes the contribution of excited electronic states. Because of the extra electronic entropy term in FT-KS-DFT framework, at finite temperature, a great care should be given to the electronic contribution as the excited states also play a role in the acceptance ratio. A first attempt of ab initio MC simulations on a metal was performed on liquid lithium in the isothermal–isobaric and canonical ensemble32 with encouraging success.
In the present paper, we present the formalism of the Gibbs ensemble coupled to FT-KS-DFT, followed by its practical implementation. We then apply the finite temperature aiGEMC method to compute the liquid–vapor equilibrium of sodium for which several experimental data sets are available.33–35 The critical densities and temperatures range from 0.175 to 0.3 g cm−3 in density and from 2485 to 2573 K in temperature. The sodium vapor is mainly an atomic gas without large clusters,36 which thus greatly simplifies the structural constraints. Moreover, as FT-KS-DFT has proven to be reliable in the study of the condensed phases of sodium,37,38 it may serve as an excellent testing ground to extend this application of FT-KS-DFT to cover the gaseous phase of a metal as well. We demonstrate that the aiGEMC offers a great way of computing a fully ab initio liquid–vapor phase diagram and of positioning the critical point.
(1) |
(2) |
(3) |
The first term is the kinetic nuclear term and the second is the effective nuclear potential term; the latter includes the free energy of the electrons of the system for a given nuclei configuration {R}N:
(4) |
The first term is the Coulomb interaction of the nuclei of charge Z, and the second term is the electronic free energy at fixed nuclear configuration:
(5) |
(6) |
The electrons are thus treated for a specific nuclear configuration and then an average must be performed over the different configurations. Note that in eqn (6), the electronic positions and momenta are operators because the electrons must be treated in the quantum mechanical framework.
At this point, we need to stress that eqn (2) is the standard expression of the classical partition function for the canonical ensemble. It means that the nuclei are treated according to the classical mechanics in an effective potential Ueff [eqn (4)] which includes all electronic degrees of freedom [eqn (5)] plus the coulombic interaction between nuclei. This effective nuclear potential Ueff encompasses the quantum behaviour of the electrons in a consistent way through the computation of the corresponding free energy contribution.
In practice, we split the computation of the partition function eqn (2) in two distinct problems. The first deals with the evaluation of the multi-dimensional (6N) integral, which is sampled using the MC Gibbs ensemble approach. The second corresponds to the computation of the effective nuclear potential energy including the free energy of electrons, which is done using FT-KS-DFT.17 In the following sub-sections we detail the methodology that we have developed in order to perform such a computation.
(7) |
Ensemble averages can be estimated by the famous MC importance sampling algorithm introduced by Metropolis et al.43 which is at the core of the GEMC approach. For obvious reasons, full monographs were dedicated to the presentation and discussion of this method and its possible improvements.47–51 Three types of trial move o(ld) → n(ew) are considered in our GEMC simulations: (i) random displacement of a randomly chosen particle, (ii) random volume rearrangement, (iii) random transfer of a randomly chosen particle between the two sub-systems. The Metropolis scheme defines the following acceptance rules as:
(8) |
• (i) particle (belonging to box j) displacement:
(9) |
• (ii) volume exchange V1(o) → V1(n) and V2(o) → V2(n):
(10) |
• (iii) particle transfer from box j = 1 to box j = 2:
(11) |
In particular for sodium we used as initial conditions a liquid box of 17 Å with 80 atoms and an empty vapor box of 25 Å for the simulations below 2000 K, and a liquid box of 18 Å with 80 atoms and an empty vapor box of 18 Å at 2000 K. One accepted configuration at 2000 K served as starting point for the simulations at 2100 K, 2200 K, 2300 K and 2400 K. At 2500 K, we start with two 40-atoms boxes of 17 Å. In addition, we have performed an extra simulation at 2000 K with 120 atoms to estimate the finite-size effect, which is found to be small at this condition (see Fig. 1). All the simulation boxes are cubic and the dimensions above define the edge of the box. The equilibrium values and their uncertainty were calculated using the autocorrelation technique.10 The error bars reported in the following sections are one-Sigma error bars.
When compared to single-box simulations, GEMC is rightfully considered to be affected by surface finite-size effects. It means that if a liquid–gas phase transition should happen, interface issues are avoided. However one should keep in mind that other finite-size effects exist, such as divergence of the correlation length and critical slowing down, and may play an important role, especially in the vicinity of a critical point. These effects are at the heart of the brilliant refinements of the “basic algorithm” that are implemented in codes designed to run GEMC simulations such as MCCCS Towhee,53 VMMC54 and CP2K.26 The limited number of atoms we have in our cells does not preclude such finite-size effects. However, as we will show in Section 3, the relatively good agreement with the experimental data is reassuring. We also chose on purpose to not include pre-biased moves or more complex moves to our algorithm in order to carefully check our implementation using VASP and the electron free energy. We thus limit the possible flaws in our sampling. It is clear that in order to have a faster algorithm one needs to consider such more intricate moves.
To compare with the Gibbs ensemble method, we also performed molecular dynamics simulations of the liquid phase in the DFT framework at several densities along the 2000 K isotherm. In order to stay consistent with the MC calculations we used the same DFT parameters for the MD simulations. The temperature was kept constant thanks to a Nosé thermostat.61 We used a fixed volume cell containing 80 atoms. The time step was set to 2 fs for a total duration of 20 ps.
The acceptance ratios reported in Table 2 are very satisfactory with typical values between 10 and 75% for each move. The 1200 K simulation shows a very low acceptance rate for the particle exchange because the temperature is very low compared to the energy barrier. Since the acceptance rate is non zero and since we ensured a long enough run we anticipate that the results are still reliable. Fig. 1 shows the evolution of different quantities as a function of the MC steps at 2000 K. Table 1 lists all the values of the thermodynamic quantities.
Temperature (K) | N vap | N liq | ρ vap (g cm−3) | ρ liq (g cm−3) | P vap (kbar) | P liq (kbar) | U vap (eV) | U liq (eV) |
---|---|---|---|---|---|---|---|---|
a We perform an extra simulation at 2000 K with 120 atoms to estimate the finite-size effect on the liquid–vapour equilibrium density. b At 2400 and 2500 K, only one phase is present. The density, pressure, effective energy of gas and liquid phase indicate the average properties in each box. | ||||||||
1200 | 0.04 ± 0.02 | 79.96 ± 0.02 | 9 × 10−5 ± 4 × 10−5 | 0.60 ± 0.01 | 2 × 10−4 ± 1 × 10−4 | −0.36 ± 0.35 | −0.011 ± 0.008 | −83.47 ± 0.16 |
1500 | 0.52 ± 0.32 | 79.48 ± 2.35 | 0.001 ± 0.0008 | 0.58 ± 0.01 | 0.003 ± 0.001 | −0.37 ± 0.34 | −0.23 ± 0.04 | −81.55 ± 0.32 |
1800 | 2.26 ± 0.67 | 78.74 ± 1.52 | 0.014 ± 0.004 | 0.54 ± 0.01 | 0.06 ± 0.03 | −0.21 ± 0.28 | −1.01 ± 0.33 | −78.94 ± 0.45 |
2000 | 4.22 ± 1.86 | 75.78 ± 1.88 | 0.03 ± 0.01 | 0.50 ± 0.02 | 0.08 ± 0.03 | 0.0 ± 0.2 | −2.33 ± 0.37 | −74.19 ± 0.67 |
2000a | 7.48 ± 2.90 | 112.52 ± 3.22 | 0.04 ± 0.01 | 0.51 ± 0.02 | 0.10 ± 0.01 | 0.20 ± 0.34 | −4.34 ± 1.51 | −111.38 ± 1.03 |
2100 | 7.20 ± 1.59 | 72.80 ± 2.42 | 0.05 ± 0.01 | 0.46 ± 0.02 | 0.12 ± 0.05 | 0.20 ± 0.23 | −69.44 ± 1.45 | −5.11 ± 1.93 |
2200 | 12.27 ± 2.35 | 59.20 ± 6.29 | 0.074 ± 0.02 | 0.43 ± ± 0.05 | 0.15 ± 0.08 | 0.44 ± 0.50 | −8.06 ± 1.54 | −55.99 ± 6.52 |
2300 | 14.47 ± 3.84 | 57.07 ± 12.22 | 0.15 ± 0.04 | 0.42 ± 0.07 | 0.13 ± 0.52 | 0.36 ± 0.67 | −10.35 ± 3.28 | −52.97 ± 11.89 |
2400b | 38.37 ± 5.44 | 41.63 ± 5.32 | 0.28 ± 0.04 | 0.26 ± 0.03 | 0.33 ± 0.13 | 0.43 ± 0.11 | −32.25 ± 0.98 | −34.63 ± 0.74 |
2500b | 38.54 ± 6.60 | 41.46 ± 7.98 | 0.31 ± 0.03 | 0.32 ± 0.04 | 0.44 ± 0.24 | 0.46 ± 0.20 | −31.42 ± 1.71 | −38.27 ± 0.84 |
Temperature [K] | Particle displacement in Box1 | Particle displacement in Box2 | Volume exchange | Atom swap acceptance ratio | |||
---|---|---|---|---|---|---|---|
Acceptance ratio | δsmax | Acceptance ratio | δsmax | Acceptance ratio | δVmax [Å3] | ||
a We perform an extra simulation at 2000 K with 120 atoms to estimate the finite-size effect on the liquid–vapour equilibrium density. | |||||||
1200 | 0.50 | 0.035 | 0.78 | 0.23 | 0.55 | 300 | 0.006 |
1500 | 0.41 | 0.052 | 0.73 | 0.16 | 0.50 | 420 | 0.07 |
1800 | 0.58 | 0.05 | 0.41 | 0.10 | 0.45 | 700 | 0.13 |
2000 | 0.54 | 0.055 | 0.50 | 0.13 | 0.60 | 550 | 0.18 |
2000a | 0.48 | 0.04 | 0.39 | 0.1 | 0.46 | 700 | 0.14 |
2100 | 0.52 | 0.07 | 0.54 | 0.1 | 0.56 | 700 | 0.21 |
2200 | 0.49 | 0.072 | 0.58 | 0.08 | 0.57 | 710 | 0.26 |
2300 | 0.61 | 0.07 | 0.45 | 0.10 | 0.52 | 800 | 0.29 |
2400 | 0.53 | 0.078 | 0.55 | 0.08 | 0.59 | 650 | 0.39 |
2500 | 0.50 | 0.078 | 0.53 | 0.08 | 0.52 | 650 | 0.37 |
As temperature increases the distributions become less and less separated, eventually preventing for a clear difference between the two phases. As we approach the critical temperature, the simulations have a higher probability of switching identity or even having two phases at once in the same simulation box. The latter is due to a comparable magnitude of the surface tension effect and the entropy contribution as already observed by Smit et al. (1989).46 It results in the appearance of three peaks in the density distribution plot. In order to better quantify the density of the three phases (gas, liquid and the mixed phase), we fitted the three peaks by three Gaussian functions (solid black line in Fig. 2) as suggested by Smit et al. (1989).46 The center of the Gaussian is assumed to be the average density of each phase and its width is the standard deviation entering in the determination of the uncertainty. The low density peak corresponds to the gaseous phase and the high density peak to the liquid phase. The middle peak, close to the average density of the two boxes is the mixed phase, and is disregarded in the liquid–vapor equilibrium analysis. At 2300 K, the fluctuations become extremely large and we thus decided to only show the results in the density plot for reference but not to use them for the fit since they offer too loose a constraint on the critical point. For 2400 and 2500 K both boxes reach a very similar equilibrium and seem identical. This means that these conditions are above the critical point.
Based on the equilibrium densities, we can plot a vapor–liquid coexistence curve as shown in Fig. 3. In general we obtain a good agreement compared to experimental data available in the literature.33 At 1200 K, the saturated liquid density is 20% lower than that of experiments, which may be due to the choice of PBE as an exchange correlation functional as already noted in previous studies of argon.31 A precise direct determination of the critical point is made difficult by the finite size of our system since the correlation length is expected to tend to infinity at the critical point; this cannot be captured within our small simulation cells. The critical point may however be approximated using the law of rectilinear diameter:62
(12) |
(13) |
Fig. 3 Liquid–vapor equilibrium of Na obtained from the ab initio Gibbs ensemble simulations (blue circles) and its comparison with the experiment33 (black stars). The black solid line is a fit for the scaling law to all data points above 2000 K, while the black dotted line denotes the extrapolation to lower temperature. The red solid and dotted line is the law of rectilinear diameter with A = 0.87 ± 0.1 and B = 0.19 ± 0.03 parameters (see text for more details). |
We obtained the critical point by applying the scaling law to all data points above 2000 K. Using our data we obtain our best fit with A = 0.87 ± 0.1 and B = 0.19 ± 0.03. The critical point lies at 2338 ± 108 K and 0.24 ± 0.03 g cm−3. We stress here that the density values at 1800 K are compatible with the extrapolation of the scaling law, providing confidence in our fit. Our theoretical critical temperature is slightly lower than the experimental value at 2573 ± 171 K,33 and the critical density is similar to the experimental value of 0.21 ± 0.02 g cm−3. We want to underline that both our values and the experimental ones are the result of extrapolations, as in both cases it is too challenging to obtain equilibrium data in the very vicinity of the critical point. With this in mind, the good agreement that we obtain between our calculations and experiments confirms the suitability of the ab initio Gibbs ensemble method for the determination of accurate coexistence curves. We also include the Clausius–Clapeyron plot of the saturated vapor pressure and density as a function of the inverse temperature in Fig. 5. We obtain a nice affine behavior for both the logarithm of the density and of the pressure on these plots. We also have a relatively good agreement with the experimental density data.
In order to check our Gibbs ensemble results, we also performed a set of ab initio molecular dynamics simulations in the canonical ensemble at 2000 K for different densities. This allows us to analyse the structure of the liquid and to determine the corresponding spinodal point.65Fig. 4 shows the variation of pressure as a function of the density along the 2000 K isotherm. This curve exhibits a clear minimum close to 0.48 g cm−3; this is the liquid spinodal point. This is the smallest density at which the liquid is metastable. Above this density the fluid is homogeneous, as shown for example in the insets of Fig. 4. At lower densities, in the unstable branch, bubbles form proving that the liquid becomes unstable. The density of the spinodal point is close yet lower than the equilibrium density of 0.50 ± 0.02 g cm−3 predicted by the Gibbs ensemble method. We thus have a full consistency between these two completely different methods ensuring the reliability of the Gibbs ensemble method.
Fig. 4 Pressure evolution during the isothermal volume expansion for sodium at 2000 K. The insets show the snapshots at 0.52 and 0.38 g cm−3, respectively. |
Fig. 5 Clausius–Clapeyron plot of the logarithm of the saturated vapor pressure and density as a function of the inverse temperature. The solid black squares are experimental date.33 The dashed line on the right graph denotes the pressures from fitting vapor density to the ideal gas law. As expected, a good agreement with calculated pressure is achieved in the low density range. In the high density region, a deviation is observed. |
(14) |
Fig. 6 shows the RDF at several temperatures as extracted from our Gibbs ensemble simulations. The main peak lies around 3.5 Å. Along the vapor–liquid equilibrium line, the position of the peak changes slightly and broadens due mostly to the temperature effects. The spherical integration of the RDF from 0 to its first minimum gives the coordination number.
At 2000 K and 0.52 g cm−3 the agreement of the RDF as obtained in the MC and in the MD simulations is excellent. As the two methods start from different initial configurations and use different paths, they give a remarkable consistent outcome as they explore the configurational space. While this is not an absolute proof that we achieved ergodicity, it strongly suggests that our simulations are satisfying it.
This method offers a new possibility to study phase transitions without computing the complete Gibbs free energy and therefore paves the way for the outstanding questions regarding phase equilibria. It is clear that the high computer cost of such an ab initio Gibbs ensemble method precludes the use of large simulation cells, at least for now. However, considering a good agreement between our results and experimental data, we can expect to obtain very satisfactory liquid–vapor equilibrium curves for other materials.
This journal is © the Owner Societies 2021 |