Evaluation of the grand-canonical partition function using expanded Wang-Landau simulations. IV. Performance of many-body force fields and tight-binding schemes for the fluid phases of silicon

In recent years, the development of force fields, that are able to take into account many-body effects,^1–15 has been the focus of intense research. This is especially crucial since the addition of many-body terms has been shown to improve greatly the accuracy of the predictions from molecular simulation calculations for a wide range of systems, from simple systems of rare gases and their mixtures,^4,12,16–27 to molecular systems,^28–39 nanostructures,^40,41 and biological systems.⁴² These effects have also been shown to become increasingly significant in highly inhomogenous systems, such as nanoconfined systems.^43–53 In the case of silicon, simulations using many-body force fields have led to a new understanding of a wide range of phenomena such as the point-defect aggregation in silicon, the nucleation and growth of silicon crystals, as well as the formation of carbon nanotubes at SiC interfaces.^54–63 Different strategies have been proposed to model many-body interactions in silicon, either relying on a purely classical approach or on a quantum approach. The idea underlying the classical many-body (CMB) force fields, such as the well-known Stillinger-Weber potential,⁶⁴ consists in using a combination of a two-body potential with an effective many-body potential (e.g., a repulsive 3-body term that favors the ideal tetrahedral angle, cosθ = 1/3, between triplets of Si atoms⁶⁴). This allows to recover the tetrahedral arrangement of Si atoms found in the condensed phases of Si. The alternative approach, used in quantum many-body (QMB) force fields, consists in evaluating the overlap between the atomic orbitals of neighboring Si atoms, as, e.g., calculated in the tight-binding Hamiltonian matrix.^65–78 In this case, the tetrahedral ordering in the condensed phases of Si directly results from the overlap between the 4 valence orbitals of Si atoms. While a comparison of the CMB and QMB approaches has been made recently on the crystalline phases of Si and on Si clusters,⁷⁹ a full assessment of the relative performance of these two classes of force fields for the fluid phases of Si and for nanoconfined Si has yet to be carried out. In this work, in order to carry out this assessment, we extend the recently developed Expanded-Wang Landau (EWL) simulations beyond classical systems.^80–82 Given the successes of tight-binding approaches^67,75–78 in computational materials science, we develop the EWL formalism to study systems modeled within tight-binding schemes. The EWL approach is an accurate and versatile scheme that allows to determine the grand-canonical partition function of systems.^80–82 This, in turn, gives a direct access to all thermodynamic properties, including the Gibbs free energy and the entropy, from the partition functions, through the application of the formalism of statistical mechanics. Using EWL simulations, we determine the thermodynamics properties of the fluid phases of Si under a wide range of conditions, i.e., at the vapor-liquid coexistence, in compressed liquids, and under a nanoscopic confinement. We consider both CMB and QMB types of models. For CMB force fields, we consider the Stillinger-Weber potential (CMB-SW)⁶⁴ and the Tersoff potential (CMB-T).⁸³ For QMB force fields, we use the Kwon model (QMB-K)⁷² and the Lenosky model (QMB-L).⁷³ Applying the EWL approach to Si, modeled with CMB or QMB force fields, provides a full picture of the impact of the two types of strategies (CMB or QMB) on the thermodynamic properties of the fluid phases of Si in a wide range of conditions and settings. The paper is organized as follows. In Sec. II, we discuss how we extend the EWL approach for QMB tight-binding systems. We also detail how the EWL approach is used in conjunction with CMB force fields. Then, we present the EWL results obtained, using both classes of model, for the grand-canonical partition function of Si in the bulk and under a nanoscopic confinement. In particular, we assess the relative performance of each model and carry out a comparison of the EWL results to the experimental data. We finally draw the main conclusions from this work in Sec. V.

In the first papers of the series,^80–82 we have developed the expanded Wang-Landau approach to determine the grand-canonical partition function of single-component systems and mixtures modeled with classical force fields. Here we extend this approach to the case of systems modeled within a tight binding scheme. The grand-canonical partition function for such a system is given by

where β = 1/k_BT, N the number of atoms, μ the chemical potential of atoms, and Q(N, V, T) is the canonical partition function given by

where Γ denotes a specific configuration of the system and Λ is the de Broglie wavelength.

To perform an accurate sampling of the grand-canonical ensemble, it is necessary to implement a very efficient scheme for the insertion/deletion of atoms. For this purpose, we have developed an approach based on the expanded grand-canonical ensemble approach,^84–97 which consists in dividing the insertion/deletion of full atoms into M stages. Recent work has shown that the implementation of efficient schemes for the insertions/deletion steps, e.g., the expanded ensemble approach within transition matrix Monte Carlo (MC) methods⁹⁸ or the continuous fraction component methods,^99,100 greatly improves the accuracy of the simulation results. Here, the combination of a Wang-Landau sampling with the expanded grand-canonical approach yields much more accurate results for the thermodynamic properties in the low temperature high density regime, most notably for the chemical potential.^80,101 In this approach, throughout the simulation, the system contains N full atoms and a fractional atom at stage l (with 0 ≤ l ≤ M − 1). The coupling (or interaction) of the fractional atom with the N full atoms depends on the value of l and will be discussed in detail in Section II C. A fractional atom at stage l = 0 is considered as void and does not interact with the N full atoms. If l is increased and its new value exceeds M, the fractional atom becomes a full atom and a new fractional atom is created, leading to a system which now has N + 1 full atoms and a new fractional atom at stage l − M. Similarly, if l is decreased and its new value is less than 0, the fractional atom is deleted and a randomly chosen full atom becomes a new fractional atom, leading to a system which now has N − 1 full atoms and a new fractional atom at stage l + M. For this system, we define a simplified expanded grand-canonical (SEGC) partition function^80–82 as

in which Q(N, V, T, l) is the canonical partition function for a system of N full atoms and a fractional atom at stage l > 0, given by

The SEGC grand-canonical partition function differs from the conventional expanded grand-canonical partition, since it does not require the use of a weighting function (usually optimized numerically for given sets of (N, l) value^87,102). As discussed in previous work,^80–82 the fact that we use a Wang-Landau sampling scheme ensures a uniform sampling of all possible (N, l) values, thereby alleviating the need for a weighting function. Finally, the mass of a fractional atom (for any value of l other than 0) is chosen to be the same as that of a full atom, so that the de Broglie wavelength of the fractional atom is the same as for a full atom (Λ_l = Λ).

The Wang-Landau sampling relies on an iterative evaluation of the biased distribution, p_bias.^{101,103–118} In the context of simulations carried out in the Monte Carlo framework, we have the following Metropolis criterion for a move from an old state (Γ_o, N_o, l_o) to a new state (Γ_n, N_n, l_n):

In the case of a single-component system in the SEGC ensemble, the joint Boltzmann distribution p(Γ, N, l) is defined as

Eq. (6) is written above for l > 0. For a void fractional particle, l = 0, the (N + 1) terms are replaced by N.

The number distribution p(N, l) can be calculated from Eq. (6) as

Finally, p_bias(Γ, N, l) = p(Γ, N, l)/p(N, l) is given by

In tight-binding (TB) schemes,^65–71 the energy includes an electronic part (calculated as the sum of single-electron energy eigenvalues of the Schrödinger equation) and a phenomenological short-ranged repulsive part (corresponding to the repulsion between the atomic core electrons and nuclei). It is given by

In Eq. (9), the repulsive energy U_R is function of the position of atoms only and is assumed to be independent from their electronic states. U_TB denotes the electronic energy, obtained from the lowest eigenvalues of the tight-binding Hamiltonian H_TB (the factor of 2 accounting for spin). The matrix elements of the Hamiltonian H_TB are obtained using the Slater-Koster formalism¹¹⁹ within the context of the two-center approximation of the tight-binding theory, with bonding occurring as a result of the coupling between pairs of neighboring atoms. In the case of silicon, a basis set of 4 atomic orbitals (s, p_x, p_y, p_z) is assigned to each atom, and these orbitals overlap with the orbitals of neighboring atoms. For a system of N atoms, the H_TB matrix has a 4N × 4N dimension and consists of 4 × 4 blocks—one block per atomic pair (i, j). Diagonal blocks (i = j) are diagonal themselves, with matrix elements along the diagonal taken to be equal to the on-site energy for the s and p orbitals. Off-diagonal blocks (i≠j) contain the hopping matrix elements, calculated from the distance-dependent tight-binding overlaps h_a(r_ij) (with a denoting either ssσ, spσ, ppσ, or ppπ).^72,120

Let us now consider a system containing N full atoms and a fractional atom at stage l (l≠0). We now define the coupling between the fractional atom and the N full atoms. The repulsive energy between a full atom and a fractional atom is rescaled with a coupling parameter ξ_l = l/M in the same way as in previous work^80,81 (potential parameters with the dimension of an energy are scaled by $ξl1/3$ and potential parameters with the dimension of a distance are scaled by $ξl1/4$⁠). The electronic energy for a system of N full atoms and a fractional atom is calculated as follows. The tight-binding matrix H_TB for such a system is now chosen as a matrix of 4(N + 1) × 4(N + 1) dimension, with the 4 additional dimensions (beyond 4N) corresponding to the fractional atom, labeled as the (N + 1)th particle in the system. The diagonal block for the fractional particle is given by

where ϵ_s and ϵ_p are the on-site energies for the s and p orbitals of the fractional atom (chosen to be the same as that of a full atom). The off-diagonal blocks for the atomic pairs involving a full atom i and the fractional atom N + 1 are given by

where $hssσf$⁠, $hspσf$⁠, $hppσf$⁠, and $hppσf$ are the hopping functions for the fractional-full interactions (see Sec. III for the detailed expression for the TB schemes used in this work) and d_α = (α_i − α_N+1)/r_i,N+1 (with α = x, y or z). Once the tight-binding matrix is defined, the electronic energy can be obtained by diagonalizing the H_TB matrix and by taking the lowest eigenvalues of H_TB, considering that both full and fractional atoms have 4 valence electrons. Finally, we add that, instead of the direct diagonalization of the H_TB matrix, linear-scaling methods can also be used to determine the electronic energy.^121–127 

In this work, we use classical many-body force fields and tight-binding models for Si. The many-body force fields studied here are the widely used Stillinger-Weber (CMB-SW)⁶⁴ and Tersoff (CMB-T)⁸³ potentials. The Stillinger-Weber potential U_CMB−SW is defined as the sum of a two-body term and of a three-body term u₃. The two-body term between two atoms i and j is given by

where ϵ, σ, A, B, p, and a are potential parameters taken from previous work.⁶⁴ 

The three-body term between 3 atoms i, j, and k is given by

where the h function is defined for r < a as, e.g., in the case of h(r_ij, r_ik, θ_jik),

where θ_jik denotes the angle between vectors r_ij and r_ik, subtended by vertex i, and where λ and γ are potential parameters.⁶⁴ The interaction between a full atom and the fractional atom is defined by rescaling the two-body and three-body energy terms between full atoms. For the CMB-SW potential, we use Eqs. (12)-(14) with the following parameters $ϵf=ξl1/3ϵ$ and $σf=ξl1/4σ$ for the full-fractional interaction.

The second many-body force field used in this work is the CMB-T potential,⁸³ which is based on a bond-order potential description of the interactions.¹²⁸ In this model, the interactions between two atoms i and j are given by

where b_ij is the bond order parameter, which is a many-body term that depends on the strength of the interaction between atoms i and j, f_c(r_ij) is a cutoff function, and A, B, λ, μ, S, and R are the CMB-T potential parameters. The bond order parameter directly depends on the bond geometry according to

where β, c, d, and h are potential parameters and θ_jik denotes the angle between vectors r_ij and r_ik. We finally define the interaction between a full atom and a fractional atom for the CMB-T potential by rescaling the repulsive and attractive terms for V(r_ij) in Eq. (15), using the following parameters: $Af=ξl1/3A$⁠, $Bf=ξl1/3B$⁠, $μf=ξl1/4μ$⁠, and $λf=ξl1/4λ$ for the full-fractional interaction.

The two tight-binding schemes studied in this work are the models of Kwon (QMB-K)⁷² and Lenosky (QMB-L).⁷³ Both are orthogonal TB models, with a minimal (s, p) basis and a repulsive potential, and have been shown to be highly transferable as they model accurately the properties of crystal phases, clusters,⁷⁹ as well as the solid-liquid and solid-solid phase boundaries of silicon.¹²⁹ The Kwon model uses the following short-range scaling functions for the TB matrix elements^72,120 for a pair of atoms (i, j),

where α corresponds to either ssσ, spσ, ppσ, or ppπ and r₀, r_cα, n_cα, and n are potential parameters taken from Kwon et al.⁷² The QMB-K repulsive energy is calculated as a sum of a functional of a repulsive pair potential¹³⁰ ϕ(r_ij),

where C_n (n = 1, 2, 3 or 4), r₀, d_c, m, and m_c are potential parameters.⁷² 

For the QMB-K model, we model the interaction between a full atom and a fractional atom as follows. We use the same functional form h_α(r_ij) defined in Eq. (17) and use the following scaled parameters: $hαf(r0)=ξl1/3hα(r0)$ (with α corresponding to either ssσ, spσ, ppσ, or ppπ), $r0f=ξl1/4r0$⁠, and $rcα=ξl1/4rcα$⁠. The electronic energy is then obtained by determining the lowest eigenvalues of the 4(N + 1) × 4(N + 1) matrix (for a system of N full atoms + 1 fractional particle), keeping in mind that the full atoms as well as the fractional atom all have 4 valence electrons. For the repulsive part, we use the same set of equations (Eq. (18)) with parameters scaled as follows: $Cnf=ξl1/3Cn$ (n = 1, 2, 3 or 4), $r0f=ξl1/4r0$⁠, and $dcf=ξl1/4dc$⁠.

The QMB-L model⁷³ is also an orthogonal TB model that uses cubic splines to represent the 4 functions defining the TB matrix elements as well as the pair repulsive potential. To model the interaction between a full atom and a fractional atom, we apply a rescaled version of the cubic splines defined by Lenosky et al. In the QMB-L model, the TB matrix elements are defined from the scaling functions h_α(r) = g_α(r)/r², where r is the interatomic distance for a pair of atoms, g_α is a cubic spline, and α one of the 4 possible overlaps (α = ssσ, spσ, ppσ, or ppπ). The scaling function for a pair including a full atom and a fractional atom is obtained by rescaling the energies by ξ_l^1/4 and the distances by ξ_l^1/3, leading to $hαf(r)=ξl1/4gα(r/ξl1/3)/(r/ξl1/3)2$⁠, and the scaled pair potential for the full-fractional repulsion is given by ϕ^f(r_ij) = ξ^1/4ϕ(r_ij/ξ^1/3). As for the QMB-K model, the potential energy is obtained by adding the repulsive energy to the electronic energy, obtained from the lowest eigenvalues of the TB matrix. Fig. 1(a) summarizes the dependence of the full-fractional pair potential on the coupling parameter. As the coupling parameter decreases, the repulsive pair potential smoothly decreases, ensuring that the insertion of the fractional atom is facilitated. Similarly, the TB matrix elements decrease with the coupling parameter, which leads to a smooth transition from a fractional atom to a full atom as it is grown during the EWL simulations. The resulting potential energy (calculated as the sum of the repulsive and electronic contributions) for dimers, composed of a full atom and of a fractional atom, is shown in Fig. 1(b) for different values of the coupling parameters. As the coupling parameter increases, the minimum is smoothly shifted from a shorter distance (reached at 1.37 Å for an energy of 1.3 eV for ξ = 0.1) to a larger distance (e.g., reached at 2.14 Å for an energy of 3 eV for ξ = 0.9), until the full Lenosky dimer energy is recovered (with a minimum of 3.4 eV reached for a distance of 2.28 Å, in very good agreement with the ab initio results¹³¹).

EWL simulations are performed on systems of up to 200 atoms within the framework of MC simulations, with the following two types of MC moves: (i) a translation of a single, full or fractional, atom (75% of the MC steps) and (ii) a change in (N, l) (25% of the MC steps). The technical details regarding the Wang-Landau scheme are the same as described in the first papers of the series.^80–82 To study the properties of the fluid phases of silicon under nanoscopic confinement, we use a slit pore geometry and model the interactions between the fluid atoms and the two confining walls with the well-known Steele 9-3 potential.^132,133 More specifically, the fluid is confined between 2 planar walls separated by a distance of 12 Å along the z-axis. The effective interaction between the atoms of the fluid with the top wall (located at S_z/2 = 6 Å along the z-axis) is given by

while the interaction with the bottom wall (located at −S_z/2 = − 6 Å along the z-axis) is

For the wall-fractional atom interaction, we simply use a rescaled version of Eqs. (19) and (20) using as interaction parameters $ϵw−Sif=ξl1/3ϵw−Si$ and $σw−Sif=ξl1/4σw−Si$⁠. We consider a graphite wall (using the parameters for carbon determined by Steele,¹³³ with a number density for the wall ρ_w = 0.097 atoms/ų taken from previous work^134,135) and model the wall-fluid interactions using the silicon parameters of Murad and Puri,¹³⁶ which gives σ_w−Si = 3.71 Å and ϵ_w−Si/k_B = 54.15 K.

We plot in Fig. 2(a) the results obtained, at a temperature T = 5000 K, for the grand-canonical partition function Θ(μ, V, T) of Si modeled with the two classical many-body force fields (CMB-SW and CMB-T) and the two quantum many-body models (QMB-K and QMB-L) considered in this work. For low chemical potentials (μ < − 31 000 kJ/kg), the values for the partition function predicted by the different models are in excellent agreement with each other. As the chemical potential increases, Fig. 2(a) shows that the grand-canonical partition sharply increases, first for the two QMB models, and then for the two CMB models. The sudden increase in Θ(μ, V, T) is first observed for the QMB-L model (μ > − 30 900 kJ/kg), then shortly after for the QMB-K model (μ > − 30 750 kJ/kg) and then for the CMB-T model (μ > − 28 350 kJ/kg) and finally for the CMB-SW model (μ > − 28 150 kJ/kg). This sharp increase in the partition function indicates the onset of a phase transition from a low density phase towards a phase of higher density (in our case here, the vapor-liquid transition for silicon). This means that the predicted value for the chemical potential at coexistence will be the lowest for the QMB-L model, followed by the QMB-K model, the CMB-T model, and finally the CMB-SW model. The predictions by the two QMB models for μ at coexistence are close (within 0.5% of each other) and around 8% below the CMB predictions. As the chemical potential further increases beyond its value at coexistence, Fig. 2(a) shows that, for all models, the grand-canonical partition function increases steadily with the chemical potential.

The grand-canonical partition function is obtained by summing up the functions Q(N, V, T), with a weighting factor of exp[βμN] (see Eq. (1)). The onset of the vapor-liquid transition, indicated by a sharp increase in the grand-canonical partition function Θ(μ, V, T), can also be identified on a plot of Q(N, V, T) against N. Fig. 2(b) shows that the values for the slopes of the logarithm of Q(N, V, T) are in the following order: QMB-L > QMB-K > CMB-T > CMB-SW. This means that the order observed for the slopes is the opposite of the order obtained for μ at coexistence shown in Fig. 2(a). This can be understood in terms of the number distribution p(N) given by

where p(N) is proportional to the weighted function Q(N, V, T)expβμN. At coexistence, the number distribution at low N (corresponding to the vapor) and at large N (corresponding to the liquid) contributes equally. This occurs when Q(N, V, T)expβμN take similar values for both phases, or, equivalently, when logQ(N, V, T) exhibits a slope that compensates (i.e., is of the opposite sign) the linear term in N, βμN. This accounts for the fact that the values for the slopes for logQ(N, V, T) are in the opposite order of that obtained for μ at coexistence.

The next step consists in using the results obtained for the partition functions to determine the thermodynamic properties of the fluid phases of Si, as predicted by the different CMB and QMB models studied here. The number distribution p(N) can be used to determine numerically the value of the chemical potential at the vapor-liquid coexistence for each model. At coexistence, the two phases have equal probabilities. Noting Π_l the probability associated with the liquid and Π_v the probability associated with the vapor, writing the condition that Π_l = Π_v leads to

with N_b denoting the boundary value (the value of N for which p(N) reaches a minimum between the two peaks corresponding to the liquid and vapor phases).

We show in Fig. 3 the results obtained for the density distribution p(ρ = N/V) at coexistence, in the case of the QMB-L model. For each temperature, p(ρ) exhibits two peaks, associated with either the vapor phase for low densities or the liquid phase for high densities. Graphing the number distribution against the density and the temperature allows us to obtain the 3-D plot shown in Fig. 3(a). This 3-D plot is especially interesting since it provides a 3-D sketch of the phase diagram, drawn on the basis of the phase probabilities, with the two series of peaks underlying the phase envelope. To give a better account of the variation, as a function of temperature, of the densities of the two phases at coexistence, we plot in Fig. 3 the probabilities obtained for each phase, either on a logarithmic scale for the vapor (see the left-hand-side of Fig. 3(b)) or on a linear scale for the liquid (see the right-hand-side of Fig. 3(b)). These plots show the expected behavior for the densities of the two phases at coexistence, with a shift of the maximum probability for the vapor towards larger values for the density as the temperature increases (with a density increase of 5 orders of magnitude as the temperature goes from 3000 K to 6500 K) and a shift of the maximum probability for the liquid towards lower densities as the temperature increases (with a decrease in density by 27% as the temperature goes from 3000 K to 6500 K).

The vapor-liquid equilibria for each model can be plotted in the temperature-density plane by reporting the densities at coexistence for the liquid, 〈ρ_l〉, and for the vapor, 〈ρ_v〉, obtained from

The densities at coexistence are shown in Fig. 4(a) for the CMB and QMB models for temperatures ranging from 3000 K to 6500 K. For comparison purposes, we also include the results⁵⁵ from prior work, obtained using the CMB-SW model and the Gibbs Ensemble Monte Carlo (GEMC) algorithm.^137,138 The GEMC results for the CMB-SW are found to be in excellent agreement with the EWL results obtained in this work. While the chemical potentials at coexistence were found to be similar within the same class of model (either CMB or QMB—see Fig. 2), a similar type of classification does not strictly apply to the results for the densities. Fig. 4(a) shows that the EWL densities are very sensitive to the model parameters. For instance, the density of the liquid phase at coexistence is found to be the highest for the QMB-K model (close to, e.g., 2.7 g/cm³ at 3000 K) and the lowest for the QMB-L model (about 2.2 g/cm³). On the other hand, the density of the vapor is found to be lower for the QMB models (with the QMB-K density being roughly an order of magnitude less than for the QMB-L model over the temperature range considered here) than for the CMB models (with the CMB-SW density being an order of magnitude less than for the CMB-T model over the temperature range). Overall, on the basis of the VLE densities, a critical comparison of the various models leads to the following conclusion. The two models exhibiting extreme behaviors are observed for the QMB-K model (highest liquid density-lowest vapor density) and for the CMB-T model (second lowest liquid density-highest vapor density), while intermediate behaviors are observed by the CMB-SW and QMB-L models.

While experimental data are not available for the densities at coexistence, the experimental data for the vapor pressure¹³⁹ provide a way to assess the performance of the force fields to model the fluid phases of Si. In EWL simulations, the pressure at coexistence can be determined from Θ(μ, V, T) through the following equation (using the value at coexistence for μ):

The EWL results for the pressure at coexistence are shown in Fig. 4(b) for the 4 models and compared to the available experimental data.¹³⁹ The results confirm the conclusions drawn from the results for the densities at coexistence, with the results for QMB-K and CMB-T appearing to be outliers. The QMB-K pressure underestimates by almost 2 orders of magnitude the experimental data at low temperatures, which is consistent with the very low density predicted for the vapor. This behavior could be attributed to the definition of the QMB-K model, which gives an incorrect energy for isolated atoms.⁷³ Such cases occur frequently and thus become of great significance at low temperatures, when the vapor density becomes very low. The CMB-T pressure consistently overestimates the experimental data as well as the pressures predicted by all other models, especially at high temperatures where the vapor density predicted by the CMB-T model appears to be too large. In line with the results for the densities, the CMB-T and QMB-K results bracket the predictions by the CMB-SW model (which are remarkably close to the experimental data) and by the QMB-L model. The QMB models also provide access to additional information that cannot be obtained when classical many-body force fields are used. These include the electronic density of state plotted in Fig. 5 for the two QMB models studied in this work. Fig. 5 shows the electronic density of states for the liquid, along the coexistence line, at T = 3000 K and T = 6000 K. The QMB-K model predicts that the electronic density of states retains the same qualitative features as temperature increases, consistently with the very moderate change in liquid density with the temperature (the liquid density at coexistence decreases by only 13.3% for the QMB-K model as the temperature increases from 3000 K to 6000 K). On the other hand, increasing the temperature results in a qualitatively different electronic density of states for the QMB-L model, with the formation of a gap around 1.9 eV at high temperature (see the right panel of Fig. 5) as a result of the larger change in liquid density observed for this model (the liquid density decreases by 24.2% for the QMB-L model compared to 13.3% for the QMB-K model, over the same temperature interval).

In the rest of the paper, we focus on the results obtained for the CMB-SW and QMB-L models, since these have been shown to provide predictions that are closest to the experiment for the vapor-liquid coexistence. We compare the results obtained for compressed liquid for the CMB-SW and QMB-L models at T = 5000 K and for pressures ranging from 0.2 GPa to 1 GPa. Comparing the two sets of results, we find that the predictions for the properties of compressed liquids exhibit similar deviations to those observed for the vapor-liquid phase diagram. We find a shift of the order of 2–3 × 10³ kJ/kg for the Gibbs free energy (in line with the results for the chemical potentials at coexistence) and larger densities for the CMB-SW model, by up to about 9% (consistently with the findings for the liquid densities at coexistence). Both models are found to capture the essential features of the thermodynamic response of Si to the increase in pressure. In particular, the effect of pressure on the density is of the same order for both models. The dependence upon pressure of the thermodynamic properties is also found to be very similar for both models, with a similar increase in Gibbs free energy over the pressure range (about 400 kJ/kg for both models) and in the enthalpy (close to 250 kJ/kg for both models) and a similar decrease in entropy over the pressure range (0.02 kJ/kg for both models).

We now apply the EWL method to determine the partition functions for nanoconfined fluid phases of Si. We begin our analysis with the results obtained for the grand-canonical partition function Θ(μ, V, T). Fig. 6(a) shows the partition functions obtained for nanoconfined Si with the two models for T = 3000 K and T = 5000 K. The grand-canonical partition functions are in very good agreement for the two models at low chemical potentials (e.g., for μ < − 22 350 kJ/kg at T = 3000 K). As the chemical potential increases, the partition functions for the models start to depart strongly from each other. The QMB-L partition function exhibits a steep increase when μ exceeds −22 350 kJ/kg, while the CMB-SW partition function increases rapidly once μ becomes greater than −20 100 kJ/kg. A similar behavior is observed at higher temperatures, as shown for T = 5000 K on the left of Fig. 6(a), with the QMB-L grand-canonical partition function increasing steeply for a chemical potential about 2900 kJ/kg less than for the CMB-SW partition function. The steep increase in the partition function is associated with a transition from a nanoconfined phase of low density to a nanoconfined phase of high density. Comparing the results for nanoconfined Si to those obtained for the bulk, we find that, for a given model (either CMB-SW or QMB-L), confining Si results in a shift of the transition, from the low density phase to the high density phase, in μ by about 650 kJ/kg at T = 5000 K. Examining the results for logQ(N, V, T) provides a direct way to locate the chemical potentials for which the phase transition takes place. At coexistence, the chemical potential is such that the linear term βμN takes the opposite value of the slope of logQ(N, V, T). For a given temperature (fixed β), Fig. 6(b) shows that the slope of logQ is greater for the QMB-L model than for the CMB-SW model. This means that the transition from the low-density nanoconfined phase to the high-density nanoconfined phase occurs at lower μ for QMB-L than for CMB-SW, consistently with the delay in μ for the steep increase of Θ_CMB−SW when compared to Θ_QMB−L (see Fig. 6(a)).

We now turn to the thermodynamics of nanoconfined fluid phases of Si. Starting with the predictions from the CMB-SW model, we plot in Fig. 7 (top) the adsorption isotherms for Si for temperatures ranging from 3000 K to 6500 K. The adsorption isotherms exhibit a step, corresponding to the transition from the low-density to the high-density phase, for a value of the chemical potential that coincides with the chemical potential marking the steep increase in logΘ in Fig. 6(a) (e.g., slightly below −20 000 kJ/kg at T = 3000 K). As the temperature increases, the adsorption step is shifted towards the lower end of the range for the chemical potential, in agreement with the behavior of logΘ as a function of T in Fig. 6(a). To analyze further the EWL results, we determine the phase diagram for nanoconfined Si by evaluating the properties of the two coexisting phases from the EWL data. For this purpose, we calculate the void volume^94,95,140 V_void, i.e., the volume of the nanopore accessible to the fluid and find a value of V_void = 0.973V_pore (V_pore being the total volume of the nanopore considered here). The density of the two nanoconfined phases is then obtained from the number distribution p(N) through Eq. (23). The phase diagram so obtained for nanoconfined Si is plotted in Fig. 7 (bottom) and compared to the EWL results for the bulk. Nanoconfinement strongly impacts the phase coexistence. For a given temperature, subjecting Si to a nanoscopic confinement leads to an increase in the density of the vapor at coexistence (e.g., by 93% at T = 6000 K) and to a decrease in the density of the liquid at coexistence (by 47% at T = 6000 K). The adsorption isotherms for Si predicted by the QMB-L model are plotted in Fig. 8(a). Fig. 8(a) shows that the adsorption step occurs concomitantly with the steep increase in logΘ seen in Fig. 6(a). This further confirms the connection between the steep increase in the partition function and the phase transition taking place in the nanoconfined fluid. Moreover, as temperature increases, the adsorption step occurs for lower values of the chemical potential, consistently with the findings for the CMB-SW model. The phase diagram for nanoconfined Si predicted by the QMB-L model is compared to the vapor-liquid equilibrium for the bulk in Fig. 8(b). As for the CMB-SW model, the QMB-L model predicts that the coexistence curve for nanoconfined Si is located inside the phase envelope of the bulk. Nanoconfinement results in sharp changes in the densities of the two coexisting phases, with an increase in the density of the vapor at coexistence (e.g., by 173% at T = 6000 K) and to a decrease in the density of the liquid at coexistence (by 51% at T = 6000 K).

The partition functions, provided by the EWL approach for the two models, can be used to predict all thermodynamic properties of adsorption, including the Gibbs free energy of adsorption, enthalpy of adsorption, and entropy of adsorption.⁸¹ We show on the left panel of Fig. 9 the predictions for these properties by the CMB-SW model at T = 3000 K and T = 4000 K. Fig. 9 is a plot of the opposite of the Gibbs free energy of adsorption, −G, against P, with −G shown as the sum of the entropic term TS and of the enthalpic term −H. For both temperatures, we observe that −G starts to decrease as the pressure increases (for P up to 0.25 bars at T = 3000 K). This is mainly due to the slow decrease in entropy, as a result of the slow density increase of the nanoconfined fluid and, to a lesser extent, to the slow increase in −H, which also increases with the density of nanoconfined Si. As P further increases (beyond 0.25 bars at T = 3000 K), the phase transition from the low-density phase to the high-density phase takes place and results in a steep drop in the entropic term and a steep rise in the enthalpic term. The two terms then reach a plateau, leading to an almost constant value for −G beyond 0.25 bars. A similar behavior is observed at higher temperature, as shown at the bottom of the left panel for T = 4000 K. However, while the entropic and enthalpic terms take similar values at T = 3000 K in the high pressure regime, the entropic term increases with T and becomes greater than the enthalpic term in the high pressure regime at T = 4000 K. The predictions obtained with the QMB-L model are shown on the right panel of Fig. 9. The Gibbs free energy of adsorption for the QMB-L model is shifted by about 10% with respect to the CMB-SW model and the transition low-density → high-density occurs at a lower pressure for the QMB-L model than for the CMB-SW model. However, there is a very good agreement between the predictions from the two models for the relative magnitude of the entropic and enthalpic contributions both at T = 3000 K and T = 4000 K. For instance, at T = 3000 K, −H and TS are found to be the same for both the CMB-SW and the QMB-L force fields in the plateau region. These sets of results confirm that the CMB-SW and QMB-L force fields, which rely on dramatically different approaches to model many-body effects, both manage to capture the essential features of the thermodynamics of nanoconfined Si.

We develop the expanded Wang-Landau method for systems modeled with a tight-binding Hamiltonian and apply the resulting method to determine the partition function and thus all thermodynamic properties of the fluid phases of Si. We consider different strategies, either classical or quantum, to take into account many-body effects in the fluid phases of silicon. These include classical many-body force fields, which mimic the tetrahedral organization in crystalline Si through a repulsive 3-body term favoring the ideal tetrahedral angle in the CMB-SW model, as well as quantum many-body tight binding models, which explicitly calculate the overlap between the 4 valence orbitals of neighboring atoms as in the QMB-T and QMB-L models. The EWL results show that the grand-canonical partition function is very sensitive to the strategy (CMB or QMB) chosen to model many-body effects. In particular, the chemical potential for which the system undergoes a vapor → liquid transition, both in the bulk and under nanoconfinement, is predicted to be 10%–15% higher for CMB models than for QMB models. When subjected to nanoconfinement, the phase diagram was shown to undergo a dramatic change, with, e.g., at 6000 K, a decrease in liquid densities by about 50% for both CMB and QMB models and an increase in vapor densities between 90% (CMB) and 170% (QMB). The results obtained in this work also allow us to rank the performance of the various models on the basis of their ability to predict the available experimental data. In particular, the CMB-SW and the QMB-L models are found to yield the results that are the closest to the experimental data for the vapor pressure. With respect to the CMB-SW model, the QMB-L model has the additional advantage of providing an insight into the dependence, upon the thermodynamic conditions, of the electronic properties of the fluid phases of silicon.

Partial funding for this research was provided by NSF through CAREER Award No. DMR-1052808.