# 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.

Caroline Desgranges and Jerome Delhommelle  Department of Chemistry, 151
Cornell Street Stop 9024, University of North Dakota, Grand Forks ND 58202

###### Abstract

We extend Expanded Wang-Landau (EWL) simulations beyond classical systems and
develop the EWL method for systems modeled with a tight-binding Hamiltonian.
We then apply the method to determine the partition function and thus all
thermodynamic properties, including the Gibbs free energy and entropy, of the
fluid phases of Si. We compare the results from quantum many-body (QMB) tight
binding models, which explicitly calculate the overlap between the atomic
orbitals of neighboring atoms, to those obtained with classical many-body
force fields (CMB), which allow to recover the tetrahedral organization in
condensed phases of Si through e.g. a repulsive 3-body term that favors the
ideal tetrahedral angle. Along the vapor-liquid coexistence, between $3000$ K
and $6000$ K, the densities for the two coexisting phases are found to vary
significantly (by $5$ orders of magnitude for the vapor and by up to $25$ %
for the liquid) and to provide a stringent test of the models. Transitions
from vapor to liquid are predicted to occur for chemical potentials that are
$10-15$ % higher for CMB models than for QMB models, and a ranking of the
force fields is provided by comparing the predictions for the vapor pressure
to the experimental data. QMB models also reveal the formation of a gap in the
electronic density of states of the coexisting liquid at high temperatures.
Subjecting Si to a nanoscopic confinement has a dramatic effect on the phase
diagram, 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 presented here provide a full picture of
the impact of the strategy (CMB or QMB) chosen to model many-body effects on
the thermodynamic properties of the fluid phases of Si.

##  I Introduction

In recent years, the development of force fields, that are able to take into
account many-body effects Tao et al. (2012); DiStasio Jr et al. (2014);
Nasrabad et al. (2004); Moosavi and Goharshadi (2008); Song et al. (2011);
Malijevskỳ et al. (2007); Nasrabad and Laghaei (2006); Jäger et al. (2011);
del Río et al. (2013); Guzman et al. (2011); Cencek et al. (2013); Wang and
Sadus (2006a); Wiebke et al. (2012); Tang et al. (2012); Wiebke et al. (2014),
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 van der Hoef and
Madden (1999); Jakse et al. (2002); Wang and Sadus (2006b); Moosavi and
Goharshadi (2008); Malijevskỳ and Malijevskỳ (2003); Goharshadi and Abbaspour
(2006); Nasrabad and Deiters (2003); Leonhard and Deiters (2000); Wang and
Sadus (2006a); Anta et al. (1997); Wang and Sadus (2006c); Vogt et al. (2001);
Marcelli and Sadus (1999); Desgranges and Delhommelle (2015), to molecular
systems Eckl et al. (2008); Tkatchenko and Scheffler (2009); Tkatchenko et al.
(2012a); Pai and Bae (2014); Kennedy et al. (2014); Tkatchenko et al. (2012b);
Bereau and von Lilienfeld (2014); Tainter et al. (2011); Babin et al. (2012);
Mas et al. (2003); McDaniel and Schmidt (2014); Schmidt et al. (2015),
nanostructures Gobre and Tkatchenko (2013); Tkatchenko (2015) and biological
systems von Lilienfeld and Tkatchenko (2010). These effects have also been
shown to become increasingly significant in highly inhomogenous systems, such
as e.g. nanoconfined systems Lee et al. (1984); Eslami and Mehdipour (2012);
Reiter and Deb (2014); Strekalova et al. (2012); Wu et al. (2007); Fullerton
and Maranas (2005); Habenicht and Paddison (2011); de Los Santos and Franzese
(2012); Kalcher et al. (2010); Li and Zeng (2012); Stroberg et al. (2012). 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 e.g. the point-defect
aggregation in Silicon, the nucleation and growth of Silicon crystals as well
as the formation of carbon nanotubes at SiC interfaces Choudhary and Clancy
(2005); Honda and Nagasaka (1999); Makhov and Lewis (2003); Lorazo et al.
(2006); Sinno (2007); Kapur et al. (2010); Desgranges and Delhommelle (2011);
Kuwahara et al. (2012); Gehrmann et al. (2015); Ogasawara et al. (2014).
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 force fields (CMB), such
as e.g. the well-known Stillinger-Weber potential Stillinger and Weber (1985),
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\theta=1/3$, between triplets of Si atoms Stillinger
and Weber (1985)). 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 Chadi and Cohen (1975); Sankey and Allen
(1986); Porezag et al. (1995); Harrison (1980); Cohen et al. (1994); Menon and
Subbaswamy (1997); Elstner et al. (1998); Kwon et al. (1994); Lenosky et al.
(1997); Cook and Clancy (1993); Ohta et al. (2008, 2009); Berdiyorov et al.
(2014); Cui and Elstner (2014). 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 Ghasemi
et al. (2010), 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 Desgranges and Delhommelle (2012a, b); Desgranges and
Delhommelle (2014). Given the successes of tight-binding approaches Porezag et
al. (1995); Ohta et al. (2008, 2009); Berdiyorov et al. (2014); Cui and
Elstner (2014) 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 Desgranges and Delhommelle
(2012a, b); Desgranges and Delhommelle (2014). 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)
Stillinger and Weber (1985) and the Tersoff potential (CMB-T) Tersoff (1989).
For QMB force fields, we use the Kwon model (QMB-K) Kwon et al. (1994) and the
Lenosky model (QMB-L) Lenosky et al. (1997). 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 the next section, 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 the last section.

##  II Expanded Wang-Landau sampling for tight-binding models

###  II.1 Theoretical framework

In the first papers of the series Desgranges and Delhommelle (2012a, b);
Desgranges and Delhommelle (2014), 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

$\Theta(\mu,V,T)=\sum_{N=0}^{\infty}Q(N,V,T)\exp(\beta\mu N)$ (1)

where $\beta=1/k_{B}T$, $N$ the number of atoms, $\mu$ the chemical potential
of atoms and $Q(N,V,T)$ is the canonical partition function given by

$Q(N,V,T)={V^{N}\over{N!\Lambda^{3N}}}\int\exp\left(-\beta
U(\mathbf{\Gamma})\right)d{\mathbf{\Gamma}}$ (2)

where $\mathbf{\Gamma}$ denotes a specific configuration of the system and
$\Lambda$ 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 Lyubartsev et al. (1992); Escobedo and de
Pablo (1996); Muller and Paul (1994); Escobedo and Abreu (2006); Abreu and
Escobedo (2006); Singh and Errington (2006); Escobedo and Martinez-Veracoechea
(2007); Escobedo (2007); Escobedo and Martinez-Veracoechea (2008); Shi and
Maginn (2008); Hicks et al. (2012); Koenig et al. (2014); Hicks et al. (2014);
Owen et al. (2015), which consists in dividing the insertion/deletion of a
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 methods Rane et al.
(2013) or the continuous fraction component methods Yee et al. (2013); Sikora
et al. (2015), 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 Desgranges and Delhommelle (2012a); Ganzenm$\ddot{\mathrm{u}}$ller
and Camp (2007). In this approach, throughout the simulation, the system
contains $N$ full atoms and a fractional atom at stage $l$ (with $0\leq l\leq
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
at stage $l-M$ 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 at stage $l+M$ is created,
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 Desgranges and Delhommelle (2012a, b);
Desgranges and Delhommelle (2014) as

$\Theta_{SEGC}(\mu,V,T)=\sum_{N=0}^{\infty}\sum_{l=0}^{M-1}Q(N,V,T,l)\exp(\beta\mu
N)\\\ $ (3)

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

$Q(N,V,T,l)={V^{N}\over{N!\Lambda^{3N}\Lambda_{l}^{3}}\int\exp\left(-\beta
U({\mathbf{\Gamma}})\right)d{\mathbf{\Gamma}}}$ (4)

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 Trebst et al. (2004); Escobedo and Abreu (2006)). As discussed in
previous work Desgranges and Delhommelle (2012a, b); Desgranges and
Delhommelle (2014), 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 ($\Lambda_{l}=\Lambda$).

###  II.2 Expanded Wang-Landau sampling.

The Wang-Landau sampling relies on an iterative evaluation of the biased
distribution, $p_{bias}$ Wang and Landau (2001a, b); Shell et al. (2002, 2003,
2004); de Pablo et al. (2002); Luettmer-Strathmann et al. (2008);
Ganzenm$\ddot{\mathrm{u}}$ller and Camp (2007); Desgranges and Delhommelle
(2009); Aleksandrov et al. (2010); Desgranges et al. (2010a, b); Malakis et
al. (2010); Do et al. (2012); Ngale et al. (2012); Desgranges et al. (2012);
Aleksandrov et al. (2012). 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 ($\mathbf{\Gamma_{o}},N_{o},l_{o}$) to a new state
($\mathbf{\Gamma_{n}},N_{n},l_{n}$)

$acc(o\to n)=min\left[1,{p_{bias}(\mathbf{\Gamma_{n}},N_{n},l_{n})\over
p_{bias}(\mathbf{\Gamma_{o}},N_{o},l_{o})}\right]$ (5)

In the case of a single-component system in the SEGC ensemble, the joint
Boltzmann distribution $p(\mathbf{\Gamma},N,l)$ is defined as

${p(\mathbf{\Gamma},N,l)}={V^{N+1}\exp\left(-\beta\left[U(\Gamma)-\mu
N\right]\right)\over{N!\Lambda^{3(N+1)}}\Theta_{SEGC}(\mu,V,T)}\\\ $ (6)

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

$p(N,l)=\int{p(\mathbf{\Gamma},N,l)}d\Gamma={Q(N,V,T,l)\exp(\beta\mu
N)\over\Theta_{SEGC}(\mu,V,T)}\\\ $ (7)

Finally, ${p_{bias}(\mathbf{\Gamma},N,l)}=p(\mathbf{\Gamma},N,l)/p(N,l)$ is
given by

${p_{bias}(\mathbf{\Gamma},N,l)}={V^{N+1}\exp\left(-\beta\left[U(\Gamma)-\mu
N\right]\right)\over{N!\Lambda^{3(N+1)}Q(N,V,T,l)}}\\\ $ (8)

###  II.3 Extension to tight-binding schemes

In tight-binding schemes Chadi and Cohen (1975); Sankey and Allen (1986);
Porezag et al. (1995); Harrison (1980); Cohen et al. (1994); Menon and
Subbaswamy (1997); Elstner et al. (1998), 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

$\begin{array}[]{lll}U&=&U_{TB}+U_{R}\\\
&=&\sum_{n}2\left<\Psi_{n}|H_{TB}|\Psi_{n}\right>+U_{R}\\\ \end{array}$ (9)

In Eq. LABEL:TB, 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-bing 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 Slater and Koster (1954) 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\times 4N$
dimension and consists of $4\times 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\neq j)$ contain the hopping matrix elements,
calculated from the distance-dependent tight-binding overlaps $h_{a}(r_{ij})$
(with $a$ denoting either $ss\sigma$, $sp\sigma$, $pp\sigma$ and $pp\pi$)
Goodwin et al. (1989); Kwon et al. (1994).

Let us now consider a system containing $N$ full atoms and a fractional atom
at stage $l$ $(l\neq 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 $\xi_{l}=l/M$ in the
same way as in previous work Desgranges and Delhommelle (2012a, b) (potential
parameters with the dimension of an energy are scaled by $\xi_{l}^{1/3}$ and
potential parameters with the dimension of a distance are scaled by
$\xi_{l}^{1/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)\times 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

$\left(\begin{array}[]{cccc}\epsilon_{s}&0&0&0\\\ 0&\epsilon_{p}&0&0\\\
0&0&\epsilon_{p}&0\\\ 0&0&0&\epsilon_{p}\\\ \end{array}\right)$ (10)

where $\epsilon_{s}$ and $\epsilon_{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

$\left(\begin{array}[]{cccc}h^{f}_{ss\sigma}&d_{x}h^{f}_{sp\sigma}&d_{y}h^{f}_{sp\sigma}&d_{z}h^{f}_{sp\sigma}\\\
-d_{x}h^{f}_{sp\sigma}&d_{x}^{2}h^{f}_{pp\sigma}+(1-d_{x}^{2})h^{f}_{pp\pi}&d_{x}d_{y}(h^{f}_{sp\sigma}-h^{f}_{sp\pi})&d_{x}d_{z}(h^{f}_{sp\sigma}-h^{f}_{sp\pi})\\\
-d_{y}h^{f}_{sp\sigma}&d_{y}d_{x}(h^{f}_{sp\sigma}-h^{f}_{sp\pi})&d_{y}^{2}h^{f}_{pp\sigma}+(1-d_{y}^{2})h^{f}_{pp\pi}&d_{y}d_{z}(h^{f}_{sp\sigma}-h^{f}_{sp\pi})\\\
-d_{z}h^{f}_{sp\sigma}&d_{z}d_{x}(h^{f}_{sp\sigma}-h^{f}_{sp\pi})&d_{z}d_{y}(h^{f}_{sp\sigma}-h^{f}_{sp\pi})&d_{z}^{2}h^{f}_{pp\sigma}+(1-d_{z}^{2})h^{f}_{pp\pi}\\\
\end{array}\right)$ (11)

where $h^{f}_{ss\sigma}$, $h^{f}_{sp\sigma}$, $h^{f}_{pp\sigma}$ and
$h^{f}_{pp\sigma}$ are the hopping functions for the fractional-full
interactions (see the ’Simulation Models’ section for the detailed expression
for the TB schemes used in this work) and
$d_{\alpha}=(\alpha_{i}-\alpha_{N+1})/r_{i,N+1}$ (with $\alpha=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 Goedecker (1999); Voter et al. (1996);
Mauri et al. (1993); Ordejón et al. (1995); Drabold and Sankey (1993); Stechel
et al. (1994); Goedecker and Colombo (1994).

##  III Simulation Models

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) Stillinger and Weber (1985) and Tersoff (CMB-T)
Tersoff (1989) 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_{3}$. The
two-body term between two atoms $i$ and $j$ is given by

$\begin{array}[]{lll}u_{2}(r_{ij})&=&A\epsilon(B(r_{ij}/\sigma)^{-p}-(r_{ij}/\sigma)^{-q})\exp\left[{\left((r_{ij}/\sigma-a\right)^{-1}}\right],(r_{ij}/\sigma)<a\\\
&=&0,(r_{ij}/\sigma)\geq a\\\ \end{array}$ (12)

where $\epsilon$, $\sigma$, $A$, $B$, $p$ and $a$ are potential parameters
taken from previous work Stillinger and Weber (1985).

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

$u_{3}({\mathbf{r}_{i}},{\mathbf{r}_{j}},{\mathbf{r}_{k}})=\epsilon\left[h(r_{ij},r_{ik},\theta_{jik})+h(r_{ji},r_{jk},\theta_{ijk})+h(r_{ki},r_{kj},\theta_{ikj})\right]$
(13)

where the $h$ function is defined for $r<a$ as e.g. in the case of
$h(r_{ij},r_{ik},\theta_{jik})$

$h(r_{ij},r_{ik},\theta_{jik})=\lambda\exp\left[\gamma{\left(r_{ij}/\sigma-a\right)^{-1}}+\gamma{\left(r_{ik}/\sigma-a\right)^{-1}}\right]\times\left(\cos\theta_{jik}+1/3\right)^{2}$
(14)

where $\theta_{jik}$ denotes the angle between vectors ${\mathbf{r}_{ij}}$ and
${\mathbf{r}_{ik}}$, subtended by vertex $i$, and where $\lambda$ and $\gamma$
are potential parameters Stillinger and Weber (1985). 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
$\epsilon_{f}=\xi_{l}^{1/3}\epsilon$ and $\sigma_{f}=\xi_{l}^{1/4}\sigma$ for
the full-fractional interaction.

The second many-body force field used in this work is the CMB-T potential
Tersoff (1989), that is based on a bond-order potential description of the
interactions Porter et al. (1997). In this model, the interactions between two
atoms $i$ and $j$ is given by

$\begin{array}[]{lll}V(r_{ij})&=&f_{c}(r_{ij})\left[A\exp\left(-\lambda
r_{ij}\right)-B\exp\left(-\mu r_{ij}\right)b_{ij}\right]\\\
f_{c}(r_{ij})&=&{1\over 2}\left[1+\cos\left({r_{ij}-R\over
S-R}\right)\right]\\\ \end{array}$ (15)

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$, $\lambda$, $\mu$, $S$ and
$R$ are the CMB-T potential parameters. The bond order parameter directly
depends on the bond geometry according to

$\begin{array}[]{lll}b_{ij}&=&{\left(1+\beta^{n}\zeta_{ij}^{n}\right)^{-1/2n}}\\\
\zeta_{ij}&=&\sum_{k\neq i,j}f_{c}(r_{ik})g(\theta_{ijk})\\\
g(\theta_{ijk})&=&1+{c^{2}\over
d^{2}}-{c^{2}\over\left[d^{2}+\left(h-\cos\theta_{ijk}\right)\right]}\\\
\end{array}$ (16)

where $\beta$, $c$, $d$ and $h$ are potential parameters and $\theta_{jik}$
denotes the angle between vectors ${\mathbf{r}_{ij}}$ and ${\mathbf{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 $A_{f}=\xi_{l}^{1/3}A$,
$B_{f}=\xi_{l}^{1/3}B$, $\mu_{f}=\xi_{l}^{1/4}\mu$ and
$\lambda_{f}=\xi_{l}^{1/4}\lambda$ for the full-fractional interaction.

The two tight-binding schemes studied in this work are the models of Kwon
(QMB-K) Kwon et al. (1994) and Lenosky (QMB-L) Lenosky et al. (1997). 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 Ghasemi et al. (2010) as well as the
solid-liquid and solid-solid phase boundaries of silicon Kaczmarski et al.
(2005). The Kwon model uses the following short-range scaling functions for
the TB matrix elements Goodwin et al. (1989); Kwon et al. (1994) for a pair of
atoms $(i,j)$

$h_{\alpha}(r_{ij})=h_{\alpha}(r_{0})\times{\left(r_{0}\over
r_{ij}\right)}^{n}\times\exp\left[n\left(-{\left(r_{ij}\over
r_{c\alpha}\right)}^{n_{c\alpha}}+{\left(r_{0}\over
r_{c\alpha}\right)}^{n_{c\alpha}}\right)\right]$ (17)

where $\alpha$ corresponds to either $ss\sigma$, $sp\sigma$, $pp\sigma$ or
$pp\pi$ and $r_{0}$, $r_{c\alpha}$, $n_{c\alpha}$ and $n$ are potential
parameters taken from Kwon et al. Kwon et al. (1994). The QMB-K repulsive
energy is calculated as a sum of a functional of a repulsive pair potential Xu
et al. (1992) $\phi(r_{ij})$

$\begin{array}[]{lll}U_{R}&=&\sum_{i}f\left[\sum_{j}\phi(r_{ij})\right]\\\
f(x)&=&C_{1}x+C_{2}x^{2}+C_{3}x^{3}+C_{4}x^{4}\\\
\phi_{(}r_{ij})&=&{\left(r_{0}\over
r_{ij}\right)}^{m}\times\exp\left[m\left(-{\left(r_{ij}\over
d_{c}\right)}^{m_{c}}+{\left(r_{0}\over
d_{c}\right)}^{m_{c}}\right)\right]\end{array}$ (18)

where $C_{n}$ ($n=1$, $2$, $3$ or $4$), $r_{0}$, $d_{c}$, $m$ and $m_{c}$ are
potential parameters Kwon et al. (1994).

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_{\alpha}(r_{ij})$ defined in Eq. 17 and use the following scaled
parameters: $h^{f}_{\alpha}(r_{0})=\xi_{l}^{1/3}h_{\alpha}(r_{0})$ (with
$\alpha$ corresponding to either $ss\sigma$, $sp\sigma$, $pp\sigma$ or
$pp\pi$), $r^{f}_{0}=\xi_{l}^{1/4}r_{0}$ and
$r_{c\alpha}=\xi_{l}^{1/4}r_{c\alpha}$. The electronic energy is then obtained
by determining the lowest eigenvalues of the $4(N+1)\times 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: $C^{f}_{n}=\xi_{l}^{1/3}C_{n}$ ($n=1$, $2$, $3$ or $4$),
$r^{f}_{0}=\xi_{l}^{1/4}r_{0}$ and $d^{f}_{c}=\xi_{l}^{1/4}d_{c}$.

The QMB-L model Lenosky et al. (1997) 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_{\alpha}(r)=g_{\alpha}(r)/r^{2}$, where $r$ is the interatomic distance for
a pair of atoms, $g_{\alpha}$ is a cubic spline and $\alpha$ one of the 4
possible overlaps ($\alpha={ss\sigma}$, ${sp\sigma}$, ${pp\sigma}$ or
${pp\pi}$). The scaling function for a pair including a full atom and a
fractional atom is obtained by rescaling the energies by ${\xi_{l}}^{1/4}$ and
the distances by ${\xi_{l}}^{1/3}$, leading to
$h^{f}_{\alpha}(r)={\xi_{l}}^{1/4}g_{\alpha}(r/{\xi_{l}}^{1/3})/(r/{\xi_{l}}^{1/3})^{2}$,
and the scaled pair potential for the full-fractional repulsion is given by:
$\phi^{f}(r_{ij})=\xi^{1/4}\phi(r_{ij}/\xi^{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
so 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,
are 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 $\xi=0.1$) to a
larger distance (e.g. reached at $2.14$ Å for an energy of $3$ eV for
$\xi=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 Fournier et al. (1992)).

EWL simulations are performed on systems of up to 200 atoms within the
framework of Monte Carlo (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 Desgranges and Delhommelle (2012a, b); Desgranges
and Delhommelle (2014). 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 Gelb et al. (1999); Steele (1973). 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

$\phi_{wf}^{t}(z)={2\pi\rho_{w}\sigma_{w-Si}^{3}\epsilon_{w-Si}\over
3}\left[{2\over 15}\left(\sigma_{w-Si}\over
S_{z}/2-z\right)^{9}-\left(\sigma_{w-Si}\over S_{z}/2-z\right)^{3}\right]$
(19)

while the interaction with the bottom wall (located at $-S_{z}/2=-6$ Å along
the $z$-axis) is

$\phi_{wf}^{b}(z)={2\pi\rho_{w}\sigma_{w-Si}^{3}\epsilon_{w-Si}\over
3}\left[{2\over 15}\left(\sigma_{w-Si}\over
S_{z}/2+z\right)^{9}-\left(\sigma_{w-Si}\over S_{z}/2+z\right)^{3}\right]$
(20)

For the wall-fractional atom interaction, we simply use a rescaled version of
Eqs. 19 and 20 using as interaction parameters
$\epsilon^{f}_{w-Si}=\xi_{l}^{1/3}\epsilon_{w-Si}$ and
$\sigma^{f}_{w-Si}=\xi_{l}^{1/4}\sigma_{w-Si}$. We consider a graphite wall
(using the parameters for carbon determined by Steele Steele (1973), with a
number density for the wall $\rho_{w}=0.097$ atoms/Å3 taken from previous work
Porcheron et al. (1999); Padilla and Toxvaerd (1994)) and model the wall-fluid
interactions using the Silicon parameters of Murad and Puri Murad and Puri
(2007), which gives $\sigma_{w-Si}=3.71$ Å and $\epsilon_{w-Si}/k_{B}=54.15$
K.

##  IV Results and Discussion

###  IV.1 Thermodynamics of the fluid phases of Si

We plot in Fig. 2(a) the results obtained, at a temperature $T=5000$ K, for
the grand-canonical partition function $\Theta(\mu,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 ($\mu<-31000$ 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 $\Theta(\mu,V,T)$ is first
observed for the QMB-L model ($\mu>-30900$ kJ/kg), then shortly after for the
QMB-K model ($\mu>-30750$ kJ/kg) and then for the CMB-T model ($\mu>-28350$
kJ/kg) and finally for the CMB-SW model ($\mu>-28150$ 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 $\mu$ 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[\beta\mu N]$ (see Eq. 1). The
onset of the vapor-liquid transition, indicated by a sharp increase in the
grand-canonical partition function $\Theta(\mu,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 $\mu$ at coexistence shown in Fig.
2(a). This can be understood in terms of the number distribution $p(N)$ given
by

$p(N)={Q(N,V,T)\exp{\beta\mu N}\over\Theta(\mu,V,T)}\\\ $ (21)

$p(N)$ is proportional to the weighted function $Q(N,V,T)\exp{\beta\mu N}$. At
coexistence, the number distribution at low $N$ (corresponding to the vapor)
and at large $N$ (corresponding to the liquid) contribute equally. This occurs
when $Q(N,V,T)\exp{\beta\mu N}$ take similar values for both phases, or,
equivalently, when $\log Q(N,V,T)$ exhibits a slope that compensates (i.e. is
of the opposite sign) the linear term in $N$, $\beta\mu N$. This accounts for
the fact that the values for the slopes for $\log Q(N,V,T)$ are in the
opposite order of that obtained for $\mu$ 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 $\Pi_{l}$ the
probability associated with the liquid and $\Pi_{v}$ the probability
associated with the vapor, writing the condition that $\Pi_{l}=\Pi_{v}$ leads
to

$\sum_{N<N_{b}}p(N)=\sum_{N>N_{b}}p(N)\\\ $ (22)

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(\rho=N/V)$ at coexistence, in the case of the QMB-L model. For each
temperature, $p(\rho)$ 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,
$<\rho_{l}>$, and for the vapor, $<\rho_{v}>$, obtained from

$\begin{array}[]{lll}<\rho_{l}>&=&{\sum_{N>N_{b}}(N/V)p(N)\over{\sum_{N>N_{b}}p(N)}}\\\
<\rho_{v}>&=&{\sum_{N<N_{b}}(N/V)p(N)\over{\sum_{N<N_{b}}p(N)}}\\\
\end{array}$ (23)

The densities at coexistence are shown in Fig. 4(a) for the CMB and QMB models
for temperatures ranging from $3000$ K and $6500$ K. For comparison purposes,
we also include the results Honda and Nagasaka (1999) from prior work,
obtained using the CMB-SW model and the Gibbs Ensemble Monte Carlo (GEMC)
algorithm Panagiotopoulos (1987, 2000). 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^{3}$ at $3000$ K)
and the lowest for the QMB-L model (about $2.2~{}g/cm^{3})$. 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 model.

While experimental data are not available for the densities at coexistence,
the experimental data for the vapor pressure Desai (1986) 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
$\Theta(\mu,V,T)$ through the following equation (using the value at
coexistence for $\mu$)

$P={k_{B}T\ln\Theta(\mu,V,T)\over V}\\\ $ (24)

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 Desai (1986). 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 Lenosky et al. (1997). 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 obtain 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 model 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$ to $3\times
10^{3}~{}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).

###  IV.2 Thermodynamics of nanoconfined Si

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 $\Theta(\mu,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 $\mu<-22350$ 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 $\mu$
exceeds $-22350$ kJ/kg, while the CMB-SW partition function increases rapidly
once $\mu$ becomes greater than $-20100$ 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 $\mu$ by about $650$ kJ/kg at $T=5000$ K. Examining the results for
$\log Q(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 $\beta\mu N$ takes the opposite value of the
slope of $\log Q(N,V,T)$. For a given temperature (fixed $\beta$), Fig 6(b)
shows that the slope of $\log Q$ 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 $\mu$ for QMB-L
than for CMB-SW, consistently with the delay in $\mu$ for the steep increase
of $\Theta_{CMB-SW}$ when compared to $\Theta_{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(a) 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\Theta$ in Fig. 6(a) (e.g. slightly below $-20000$ 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\Theta$ 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 the two coexisting phases
from the EWL data. For this purpose, we calculate the void volume Myers and
Monson (2002); Hicks et al. (2012); Koenig et al. (2014) $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(b) 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\Theta$ 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 Desgranges and Delhommelle (2012b). 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 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$ bar 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$ bar 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$ bar. 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. 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$\to$ 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.

##  V Conclusions

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 that 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$\to$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.

Acknowledgements

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

## References

  * Tao et al. (2012) J. Tao, J. P. Perdew, and A. Ruzsinszky, Proc. Natl. Acad. Sci. 109, 18 (2012).
  * DiStasio Jr et al. (2014) R. A. DiStasio Jr, V. Gobre, and A. Tkatchenko, J. Phys. Condens. Matt. 26, 213202 (2014).
  * Nasrabad et al. (2004) A. E. Nasrabad, R. Laghaei, and U. K. Deiters, J. Chem. Phys. 121, 6423 (2004).
  * Moosavi and Goharshadi (2008) M. Moosavi and E. K. Goharshadi, Fluid Phase Equil. 274, 51 (2008).
  * Song et al. (2011) B. Song, X. Wang, J. Wu, and Z. Liu, Mol. Phys. 109, 1607 (2011).
  * Malijevskỳ et al. (2007) A. Malijevskỳ, F. Karlickỳ, R. Kalus, and A. Malijevskỳ, J. Phys. Chem. C 111, 15565 (2007).
  * Nasrabad and Laghaei (2006) A. E. Nasrabad and R. Laghaei, J. Chem. Phys. 125, 084510 (2006).
  * Jäger et al. (2011) B. Jäger, R. Hellmann, E. Bich, and E. Vogel, J. Chem. Phys. 135, 084308 (2011).
  * del Río et al. (2013) F. del Río, E. Díaz-Herrera, O. Guzmán, J. A. Moreno-Razo, and J. E. Ramos, J. Chem. Phys. 139, 184503 (2013).
  * Guzman et al. (2011) O. Guzman, F. Del Rio, and J. Eloy Ramos, Mol. Phys. 109, 955 (2011).
  * Cencek et al. (2013) W. Cencek, G. Garberoglio, A. H. Harvey, M. O. McLinden, and K. Szalewicz, J. Phys. Chem. A 117, 7542 (2013).
  * Wang and Sadus (2006a) L. Wang and R. J. Sadus, J. Chem. Phys. 123, 074503 (2006a).
  * Wiebke et al. (2012) J. Wiebke, E. Pahl, and P. Schwerdtfeger, J. Chem. Phys. 137, 064702 (2012).
  * Tang et al. (2012) L.-Y. Tang, Z.-C. Yan, T.-Y. Shi, J. F. Babb, and J. Mitroy, J. Chem. Phys. 136, 104104 (2012).
  * Wiebke et al. (2014) J. Wiebke, M. Wormit, R. Hellmann, E. Pahl, and P. Schwerdtfeger, J. Phys. Chem. B 118, 3392 (2014).
  * van der Hoef and Madden (1999) M. A. van der Hoef and P. A. Madden, J. Chem. Phys. 111, 1520 (1999).
  * Jakse et al. (2002) N. Jakse, J. Bomont, and J. Bretonnet, J. Chem. Phys. 116, 8504 (2002).
  * Wang and Sadus (2006b) L. Wang and R. J. Sadus, J. Chem. Phys. 125, 074503 (2006b).
  * Malijevskỳ and Malijevskỳ (2003) A. Malijevskỳ and A. Malijevskỳ, Mol. Phys. 101, 3335 (2003).
  * Goharshadi and Abbaspour (2006) E. K. Goharshadi and M. Abbaspour, J. Chem. Theory Comput. 2, 920 (2006).
  * Nasrabad and Deiters (2003) A. E. Nasrabad and U. K. Deiters, J. Chem. Phys. 119, 947 (2003).
  * Leonhard and Deiters (2000) K. Leonhard and U. K. Deiters, Mol. Phys. 98, 1603 (2000).
  * Anta et al. (1997) J. Anta, E. Lomba, and M. Lombardero, Phys. Rev. E 55, 2707 (1997).
  * Wang and Sadus (2006c) L. Wang and R. J. Sadus, Phys. Rev. E 74, 021202 (2006c).
  * Vogt et al. (2001) P. S. Vogt, R. Liapine, B. Kirchner, A. J. Dyson, H. Huber, G. Marcelli, and R. J. Sadus, Phys. Chem. Chem. Phys. 3, 1297 (2001).
  * Marcelli and Sadus (1999) G. Marcelli and R. J. Sadus, J. Chem. Phys. 111, 1533 (1999).
  * Desgranges and Delhommelle (2015) C. Desgranges and J. Delhommelle, J. Chem. Theory Comput. 11, 5401 (2015).
  * Eckl et al. (2008) B. Eckl, J. Vrabec, and H. Hasse, J. Phys. Chem. B 112, 12710 (2008).
  * Tkatchenko and Scheffler (2009) A. Tkatchenko and M. Scheffler, Phys. Rev. Lett. 102, 073005 (2009).
  * Tkatchenko et al. (2012a) A. Tkatchenko, R. A. DiStasio Jr, R. Car, and M. Scheffler, Phys. Rev. Lett. 108, 236402 (2012a).
  * Pai and Bae (2014) S. J. Pai and Y. C. Bae, J. Chem. Phys. 141, 064303 (2014).
  * Kennedy et al. (2014) M. R. Kennedy, A. R. McDonald, A. E. DePrince III, M. S. Marshall, R. Podeszwa, and C. D. Sherrill, J. Chem. Phys. 140, 121104 (2014).
  * Tkatchenko et al. (2012b) A. Tkatchenko, D. Alfé, and K. S. Kim, J. Chem. Theory Comput. 8, 4317 (2012b).
  * Bereau and von Lilienfeld (2014) T. Bereau and O. A. von Lilienfeld, J. Chem. Phys. 141, 034101 (2014).
  * Tainter et al. (2011) C. Tainter, P. Pieniazek, Y.-S. Lin, and J. Skinner, J. Chem. Phys. 134, 184501 (2011).
  * Babin et al. (2012) V. Babin, G. R. Medders, and F. Paesani, J. Phys. Chem. Lett. 3, 3765 (2012).
  * Mas et al. (2003) E. M. Mas, R. Bukowski, and K. Szalewicz, J. Chem. Phys. 118, 4404 (2003).
  * McDaniel and Schmidt (2014) J. G. McDaniel and J. R. Schmidt, J. Phys. Chem. B 118, 8042 (2014).
  * Schmidt et al. (2015) J. Schmidt, K. Yu, and J. G. McDaniel, Acc. Chem. Res. 48, 548 (2015).
  * Gobre and Tkatchenko (2013) V. V. Gobre and A. Tkatchenko, Nat. Commun. 4, 2341 (2013).
  * Tkatchenko (2015) A. Tkatchenko, Adv. Funct. Mater. 25, 2054 (2015).
  * von Lilienfeld and Tkatchenko (2010) O. A. von Lilienfeld and A. Tkatchenko, J. Chem. Phys. 132, 234109 (2010).
  * Lee et al. (1984) C.-Y. Lee, J. A. McCammon, and P. Rossky, J. Chem. Phys. 80, 4448 (1984).
  * Eslami and Mehdipour (2012) H. Eslami and N. Mehdipour, J. Chem. Phys. 137, 144702 (2012).
  * Reiter and Deb (2014) G. Reiter and A. Deb, J. Phys. Conf. Series 571, 012001 (2014).
  * Strekalova et al. (2012) E. Strekalova, M. Mazza, H. Stanley, and G. Franzese, J. Phys. Condens. Matt. 24, 064111 (2012).
  * Wu et al. (2007) Y.-C. Wu, J.-S. Lin, S.-P. Ju, W.-J. Lee, Y.-S. Lin, and C.-C. Hwang, Computational Materials Science 39, 359 (2007).
  * Fullerton and Maranas (2005) S. K. Fullerton and J. K. Maranas, Nano Lett. 5, 363 (2005).
  * Habenicht and Paddison (2011) B. F. Habenicht and S. J. Paddison, J. Phys. Chem. B 115, 10826 (2011).
  * de Los Santos and Franzese (2012) F. de Los Santos and G. Franzese, Phys. Rev. E 85, 010602 (2012).
  * Kalcher et al. (2010) I. Kalcher, J. C. Schulz, and J. Dzubiella, J. Chem. Phys. 133, 164511 (2010).
  * Li and Zeng (2012) H. Li and X. C. Zeng, ACS nano 6, 2401 (2012).
  * Stroberg et al. (2012) W. Stroberg, S. Keten, and W. K. Liu, Langmuir 28, 14488 (2012).
  * Choudhary and Clancy (2005) D. Choudhary and P. Clancy, J. Chem. Phys. 122, 174509 (2005).
  * Honda and Nagasaka (1999) N. Honda and Y. Nagasaka, Int. J. Thermophys. 20, 837 (1999).
  * Makhov and Lewis (2003) D. V. Makhov and L. J. Lewis, Phys. Rev. B 67, 153202 (2003).
  * Lorazo et al. (2006) P. Lorazo, L. J. Lewis, and M. Meunier, Phys. Rev. B 73, 134108 (2006).
  * Sinno (2007) T. Sinno, J. Cryst. Growth 303, 5 (2007).
  * Kapur et al. (2010) S. S. Kapur, A. M. Nieves, and T. Sinno, Phys. Rev. B 82, 045206 (2010).
  * Desgranges and Delhommelle (2011) C. Desgranges and J. Delhommelle, J. Am. Chem. Soc. 133, 2872 (2011).
  * Kuwahara et al. (2012) T. Kuwahara, H. Ito, Y. Higuchi, N. Ozawa, and M. Kubo, J. Phys. Chem. C 116, 12525 (2012).
  * Gehrmann et al. (2015) J. Gehrmann, D. Pettifor, A. Kolmogorov, M. Reese, M. Mrovec, C. Elsässer, and R. Drautz, Phys. Rev. B 91, 054109 (2015).
  * Ogasawara et al. (2014) N. Ogasawara, W. Norimatsu, S. Irle, and M. Kusunoki, Chem. Phys. Lett. 595, 266 (2014).
  * Stillinger and Weber (1985) F. H. Stillinger and T. A. Weber, Phys. Rev. B 31, 5262 (1985).
  * Chadi and Cohen (1975) D. Chadi and M. Cohen, Phys. Status Solidi B 68, 405 (1975).
  * Sankey and Allen (1986) O. F. Sankey and R. E. Allen, Phys. Rev. B 33, 7164 (1986).
  * Porezag et al. (1995) D. Porezag, T. Frauenheim, T. Köhler, G. Seifert, and R. Kaschner, Phys. Rev. B 51, 12947 (1995).
  * Harrison (1980) W. A. Harrison, _Electronic structure and the properties of solids_ (Freeman, San Francisco, 1980).
  * Cohen et al. (1994) R. E. Cohen, M. J. Mehl, and D. A. Papaconstantopoulos, Phys. Rev. B 50, 14694 (1994).
  * Menon and Subbaswamy (1997) M. Menon and K. Subbaswamy, Phys. Rev. B 55, 9231 (1997).
  * Elstner et al. (1998) M. Elstner, D. Porezag, G. Jungnickel, J. Elsner, M. Haugk, T. Frauenheim, S. Suhai, and G. Seifert, Phys. Rev. B 58, 7260 (1998).
  * Kwon et al. (1994) I. Kwon, R. Biswas, C. Wang, K. Ho, and C. Soukoulis, Phys. Rev. B 49, 7242 (1994).
  * Lenosky et al. (1997) T. J. Lenosky, J. D. Kress, I. Kwon, A. F. Voter, B. Edwards, D. F. Richards, S. Yang, and J. B. Adams, Phys. Rev. B 55, 1528 (1997).
  * Cook and Clancy (1993) S. J. Cook and P. Clancy, Phys. Rev. B 47, 7686 (1993).
  * Ohta et al. (2008) Y. Ohta, Y. Okamoto, S. Irle, and K. Morokuma, ACS Nano 2, 1437 (2008).
  * Ohta et al. (2009) Y. Ohta, Y. Okamoto, S. Irle, and K. Morokuma, Carbon 47, 1270 (2009).
  * Berdiyorov et al. (2014) G. Berdiyorov, M. Neek-Amal, F. Peeters, and A. C. van Duin, Phys. Rev. B 89, 024107 (2014).
  * Cui and Elstner (2014) Q. Cui and M. Elstner, Phys. Chem. Chem. Phys. 16, 14368 (2014).
  * Ghasemi et al. (2010) S. A. Ghasemi, M. Amsler, R. G. Hennig, S. Roy, S. Goedecker, T. J. Lenosky, C. Umrigar, L. Genovese, T. Morishita, and K. Nishio, Phys. Rev. B 81, 214107 (2010).
  * Desgranges and Delhommelle (2012a) C. Desgranges and J. Delhommelle, J. Chem. Phys. 136, 184107 (2012a).
  * Desgranges and Delhommelle (2012b) C. Desgranges and J. Delhommelle, J. Chem. Phys. 136, 184108 (2012b).
  * Desgranges and Delhommelle (2014) C. Desgranges and J. Delhommelle, J. Chem. Phys. 140, 104109 (2014).
  * Tersoff (1989) J. Tersoff, Phys. Rev. B 39, 5566 (1989).
  * Lyubartsev et al. (1992) A. P. Lyubartsev, A. A. Martsinovski, S. V. Shevkunov, and P. N. Vorontsov-Velyaminov, J. Chem. Phys. 96, 1776 (1992).
  * Escobedo and de Pablo (1996) F. Escobedo and J. J. de Pablo, J. Chem. Phys. 105, 4391 (1996).
  * Muller and Paul (1994) M. Muller and W. Paul, J. Chem. Phys. 100, 719 (1994).
  * Escobedo and Abreu (2006) F. A. Escobedo and C. R. A. Abreu, J. Chem. Phys. 124, 104110 (2006).
  * Abreu and Escobedo (2006) C. R. A. Abreu and F. A. Escobedo, J. Chem. Phys. 124, 054116 (2006).
  * Singh and Errington (2006) J. K. Singh and J. R. Errington, J. Phys. Chem. B 110, 1369 (2006).
  * Escobedo and Martinez-Veracoechea (2007) F. A. Escobedo and F. J. Martinez-Veracoechea, J. Chem. Phys. 127, 174103 (2007).
  * Escobedo (2007) F. A. Escobedo, J. Chem. Phys. 127, 174104 (2007).
  * Escobedo and Martinez-Veracoechea (2008) F. A. Escobedo and F. J. Martinez-Veracoechea, J. Chem. Phys. 129, 154107 (2008).
  * Shi and Maginn (2008) W. Shi and E. J. Maginn, J. Comp. Chem. 29, 2520 (2008).
  * Hicks et al. (2012) J. M. Hicks, C. Desgranges, and J. Delhommelle, J. Phys. Chem. C 116, 22938 (2012).
  * Koenig et al. (2014) A. R. V. Koenig, C. Desgranges, and J. Delhommelle, Molec. Simul. 40, 71 (2014).
  * Hicks et al. (2014) E. A. Hicks, C. Desgranges, and J. Delhommelle, Molec. Simul. 40, 656 (2014).
  * Owen et al. (2015) A. N. Owen, C. Desgranges, and J. Delhommelle, Fluid Phase Equil. 402, 69 (2015).
  * Rane et al. (2013) K. S. Rane, S. Murali, and J. R. Errington, J. Chem. Theory Comput. 9, 2552 (2013).
  * Sikora et al. (2015) B. J. Sikora, Y. J. Colòn, and R. Q. Snurr, Molec. Simul. 41, 1339 (2015).
  * Yee et al. (2013) P. Yee, J. K. Shah, and E. J. Maginn, J. Phys. Chem. B 117, 12556 (2013).
  * Ganzenm$\ddot{\mathrm{u}}$ller and Camp (2007) G. Ganzenm$\ddot{\mathrm{u}}$ller and P. J. Camp, J. Chem. Phys. 127, 154504 (2007).
  * Trebst et al. (2004) S. Trebst, D. A. Huse, and M. Troyer, Phys. Rev. E 70, 046701 (2004).
  * Wang and Landau (2001a) F. G. Wang and D. P. Landau, Phys. Rev. E 64, 056101 (2001a).
  * Wang and Landau (2001b) F. G. Wang and D. P. Landau, Phys. Rev. Lett. 86, 2050 (2001b).
  * Shell et al. (2002) M. S. Shell, P. G. Debenedetti, and A. Z. Panagiotopoulos, Phys. Rev. E 66, 056703 (2002).
  * Shell et al. (2003) M. S. Shell, P. G. Debenedetti, and A. Z. Panagiotopoulos, J. Chem. Phys. 119, 9406 (2003).
  * Shell et al. (2004) M. S. Shell, P. G. Debenedetti, and A. Z. Panagiotopoulos, J. Phys. Chem. B 108, 19748 (2004).
  * de Pablo et al. (2002) J. J. de Pablo, Q. Yan, and R. Faller, J. Chem. Phys. 116, 8649 (2002).
  * Luettmer-Strathmann et al. (2008) J. Luettmer-Strathmann, F. Rampf, W. Paul, and K. Binder, J. Chem. Phys. 128, 064903 (2008).
  * Desgranges and Delhommelle (2009) C. Desgranges and J. Delhommelle, J. Chem. Phys. 130, 244109 (2009).
  * Aleksandrov et al. (2010) T. Aleksandrov, C. Desgranges, and J. Delhommelle, Fluid Phase Eq. 287, 79 (2010).
  * Desgranges et al. (2010a) C. Desgranges, E. A. Kastl, T. Aleksandrov, and J. Delhommelle, Molec. Simul. 36, 544 (2010a).
  * Desgranges et al. (2010b) C. Desgranges, J. M. Hicks, A. Magness, and J. Delhommelle, Mol. Phys. 108, 151 (2010b).
  * Malakis et al. (2010) A. Malakis, A. N. Berker, I. A. Hijagapiou, N. G. Fytas, and T. Papakonstantinou, Phys. Rev. E 81, 041113 (2010).
  * Do et al. (2012) H. Do, J. D. Hirst, and R. J. Wheatley, J. Phys. Chem. B 116, 4535 (2012).
  * Ngale et al. (2012) K. N. Ngale, C. Desgranges, and J. Delhommelle, Molec. Simul. 36, 653 (2012).
  * Desgranges et al. (2012) C. Desgranges, K. N. Ngale, and J. Delhommelle, Fluid Phase Equil. 322-323, 92 (2012).
  * Aleksandrov et al. (2012) T. Aleksandrov, C. Desgranges, and J. Delhommelle, Molec. Simul. 38, 1265 (2012).
  * Slater and Koster (1954) J. C. Slater and G. F. Koster, Phys. Rev. 94, 1498 (1954).
  * Goodwin et al. (1989) L. Goodwin, A. Skinner, and D. Pettifor, Europhys. Lett. 9, 701 (1989).
  * Goedecker (1999) S. Goedecker, Rev. Mod. Phys. 71, 1085 (1999).
  * Voter et al. (1996) A. Voter, J. Kress, and R. Silver, Phys. Rev. B 53, 12733 (1996).
  * Mauri et al. (1993) F. Mauri, G. Galli, and R. Car, Phys. Rev. B 47, 9973 (1993).
  * Ordejón et al. (1995) P. Ordejón, D. A. Drabold, R. M. Martin, and M. P. Grumbach, Phys. Rev. B 51, 1456 (1995).
  * Drabold and Sankey (1993) D. A. Drabold and O. F. Sankey, Phys. Rev. Lett. 70, 3631 (1993).
  * Stechel et al. (1994) E. Stechel, A. Williams, and P. J. Feibelman, Phys. Rev. B 49, 10088 (1994).
  * Goedecker and Colombo (1994) S. Goedecker and L. Colombo, Phys. Rev. Lett. 73, 122 (1994).
  * Porter et al. (1997) L. J. Porter, J. Li, and S. Yip, J. Nucl. Mater. 246, 53 (1997).
  * Kaczmarski et al. (2005) M. Kaczmarski, O. N. Bedoya-Martinez, and E. R. Hernández, Phys. Rev. Lett. 94, 095701 (2005).
  * Xu et al. (1992) C. H. Xu, C. Z. Wang, C. T. Chan, and K. M. Ho, J. Phys. Condens. Matter 4, 4047 (1992).
  * Fournier et al. (1992) R. Fournier, S. B. Sinnott, and A. D. DePristo, J. Chem. Phys. 97, 4149 (1992).
  * Gelb et al. (1999) L. D. Gelb, K. Gubbins, R. Radhakrishnan, and M. Sliwinska-Bartkowiak, Rep. Prog. Phys. 62, 1573 (1999).
  * Steele (1973) W. A. Steele, Surf. Sci. 36, 317 (1973).
  * Porcheron et al. (1999) F. Porcheron, M. Schöen, and A. H. Fuchs, Phys. Chem. Chem. Phys. 1, 4083 (1999).
  * Padilla and Toxvaerd (1994) P. Padilla and S. Toxvaerd, J. Chem. Phys. 101, 1490 (1994).
  * Murad and Puri (2007) S. Murad and I. K. Puri, Nano Lett. 7, 707 (2007).
  * Panagiotopoulos (1987) A. Z. Panagiotopoulos, Mol. Phys. 61, 813 (1987).
  * Panagiotopoulos (2000) A. Z. Panagiotopoulos, J. Phys. Condens. Matter 12, R25 (2000).
  * Desai (1986) P. D. Desai, J. Phys. Chem. Ref. Data 15, 967 (1986).
  * Myers and Monson (2002) A. L. Myers and P. A. Monson, Langmuir 18, 10261 (2002).

List of figures and captions

FIG. 1: Comparison of fractional-full interactions for different values of the
coupling parameter $\xi$ to the full-full interaction for the Lenosky model.
(a) (Top) Repulsive pair potential and (Bottom) TB matrix elements
$h_{ss\sigma}$ (dashed lines) and $h_{sp\sigma}$ (solid lines), and (b) Dimer
potential energy.

FIG. 2: Si at $T=5000$ K. (a) Logarithm of the grand-canonical partition
function $\Theta(\mu,V,T)$ of Si for the classical force fields, with the
results for the SW model shown in black and the results for the Tersoff model
shown in red, and for the tight-binding models, with the results for the Kwon
model shown in green and the results for the Lenosky model shown in blue. (b)
Logarithm of the function $Q(N,V,T)$ (same legend as in (a)).

FIG. 3: Density probability $p(\rho)$ obtained with the EWL method for the
Lenosky model. (a) 3D-plot showing the results obtained at coexistence as a
function of the density and temperature, with the peak corresponding to the
vapor phase on the left and the peak corresponding to the liquid phase on the
right. (b) Probabilities are shown using a logarithmic scale for the density
in the case of the vapor phase (left), and using a linear scale in the case of
the liquid phase (right).

FIG. 4: Vapor-liquid equilibrium properties for Silicon. (a) Densities of the
two fluid phases at coexistence: EWL results are shown for the SW model
(circles), for the Tersoff model (diamonds), for the Kwon model (triangles up)
and for the Lenosky model (squares), and compared to previous work using GEMC
simulations for the SW model Honda and Nagasaka (1999) (crosses). (b) Vapor
pressure results obtained using the EWL method for the 4 models (same symbols
as in (a)). The inset shows a comparison in the low-temperature range to the
experimental data Desai (1986) (solid black line with stars).

FIG. 5: Electronic density of states for the Kwon model (left) and the Lenosky
model (right) for the liquid phase at coexistence. Results obtained at 3000 K
are shown in black, while results obtained at 6000 K are shown in red.

FIG. 6: (a) Logarithm of the grand-canonical partition function
$\Theta(\mu,V,T)$ for nanoconfined Si. Results for the SW model are shown in
black for $3000$ K (solid line) and for $5000$ K (dashed line). Results for
the Lenosky model are shown in red for $3000$ K (solid line with filled
squares) and for $5000$ K (dashed line with open squares). (b) Logarithm of
the function $Q(N,V,T)$ (same legend as in (a)).

FIG. 7: Adsorption of Silicon in a slit nanopore for the SW model. (Top)
Adsorption isotherms showing the variation of the number of Si atoms adsorbed
as a function of the chemical potential. (Bottom) Phase diagram for
nanoconfined Si (in red), with a comparison to the bulk (in black).

FIG. 8: Adsorption of Silicon in a slit nanopore for the Lenosky model. (Top)
Adsorption isotherms showing the variation of the number of Si atoms adsorbed
as a function of the chemical potential. (Bottom) Vapor liquid phase
equilibria for nanoconfined Si (in red), with a comparison to the bulk (in
black).

FIG. 9: Adsorption of Silicon in a slit nanopore. The left panel shows the
results for the SW model at $3000$ K (top) and $4000$ K (bottom) for the Gibbs
free energy of adsorption (black line), the entropy of adsorption (black
triangles up with a dashed line) and the enthalpy of adsorption (red circles
with a solid line). The right panel shows the results for the Lenosky model at
$3000$ K (top) and $4000$ K (bottom), with the same legend as for the left
panel.

(a) (b)

Figure 1: DESGRANGES-DELHOMMELLE

(a) (b)

Figure 2: DESGRANGES-DELHOMMELLE

(a) (b)

Figure 3: DESGRANGES-DELHOMMELLE

(a) (b)

Figure 4: DESGRANGES-DELHOMMELLE Figure 5: DESGRANGES-DELHOMMELLE

(a) (b)

Figure 6: DESGRANGES-DELHOMMELLE Figure 7: DESGRANGES-DELHOMMELLE Figure 8:
DESGRANGES-DELHOMMELLE Figure 9: DESGRANGES-DELHOMMELLE


cell data, we run DESeq with parameters test=LRT, fitType=glmGamPoi,
useT=TRUE, minmu=1e-6, minReplicatesForReplace=Inf, reduced= 1. We repeat this
for each knockdown independently in turn.

At the end of this procedure, we have $100^{2}$ unadjusted p-values. We
perform Benjamini-Hochberg correction to an FDR level of 0.05. Subsequently,
in $\text{ko}(i)$, if gene $j$ has an adjusted p-value less than .05, we call
that a true edge in downstream calculations of precision and recall.

####  7.10.2 Protein-Protein Interactions

As a second validation edge set, we use protein-protein interactions from the
STRING database [28]. We pull protein-protein interactions for Homo Sapiens
down, and use only physical evidence for interactions. Protein-protein
interactions were only used if the STRING confidence score was over 70%.

####  7.10.3 Precision Recall Tables

In this section, we present precision and recall at three edge threshold
levels ($|w|>0.2$, $|w|<0.1$, and $|w|<0.05$) for all methods. We use both
differential expression calls and protein-protein interactions as true sets,
and separate results for both. Tables 4, 5, and 6 denote results for
$|w|>0.2$, $|w|<0.1$, and $|w|<0.05$, respectively. For differential
expression calls we ignore causal direction, and instead for two given genes
denote any differential expression in either direction as a “true” edge. This
is because all methods struggle equally with inferring causal direction under
differential expression, and call an equal number of edges in the “correct”
causal direction as the “incorrect” anticausal direction (see Supplementary
Table 7).

Method | # Called Edges | % DE | % PPI | % DE | PPI | % Recall
---|---|---|---|---|---
GOLEM-EV | 18 | 44.4 | 77.8 | 83.3 | 0.813
NO TEARS | 20 | 0.5 | 0.8 | 0.85 | 0.921
GOLEM-NV | 21 | 42.9 | 76.2 | 81.0 | 0.921
DirectLiNGAM | 22 | 40.9 | 68.2 | 77.3 | 0.921
sortnregress | 31 | 35.5 | 64.5 | 67.7 | 1.12
dotears | 32 | 46.9 | 65.6 | 75.0 | 1.14
DCDI-G | 2039 | 30.0 | 4.46 | 32.6 | 44.0
GIES | 3038 | 29.2 | 15.4 | 39.2 | 64.5
IGSP | 3064 | 29.0 | 15.2 | 38.9 | 64.6
UT-IGSP | 3075 | 29.8 | 15.1 | 39.6 | 66.0
Table 4: Validation results for all methods at a threshold of $|w|<0.2$. Columns in order are: Method name (Method), number of edges called at threshold (# Called Edge), percent of called edges showing differential expression in either the causal or anticausal direction (% DE), percent of called edges with protein-protein interactions at any confidence level (% PPI), percent of edges called either differentially expressed or have protein-protein interactions (% DE | PPI), and total recall, counting both differential expression and protein-protein interactions as true positives (% Recall). Method | # Called Edges | % DE | % PPI | % DE | PPI | % Recall
---|---|---|---|---|---
NO TEARS | 69 | 40.6 | 50.7 | 66.7 | 2.49
GOLEM-EV | 71 | 39.4 | 50.7 | 66.2 | 2.55
GOLEM-NV | 72 | 38.9 | 50.0 | 63.9 | 2.49
dotears | 92 | 37.0 | 45.7 | 60.9 | 3.03
DirectLiNGAM | 121 | 37.2 | 41.3 | 61.2 | 4.01
sortnregress | 131 | 36.6 | 41.2 | 58.8 | 4.17
DCDI-G | 2039 | 30.0 | 4.46 | 32.6 | 44.0
GIES | 3038 | 29.2 | 15.4 | 39.2 | 64.5
IGSP | 3064 | 29.0 | 15.2 | 38.9 | 64.6
UT-IGSP | 3075 | 29.8 | 15.1 | 39.6 | 66.0
Table 5: Validation results for all methods at a threshold of $|w|<0.1$. Columns in order are: Method name (Method), number of edges called at threshold (# Called Edge), percent of called edges showing differential expression in either the causal or anticausal direction (% DE), percent of called edges with protein-protein interactions at any confidence level (% PPI), percent of edges called either differentially expressed or have protein-protein interactions (% DE | PPI), and total recall, counting both differential expression and protein-protein interactions as true positives (% Recall). Method | # Called Edges | % DE | % PPI | % DE | PPI | % Recall
---|---|---|---|---|---
NO TEARS | 149 | 38.3 | 40.9 | 61.7 | 4.98
GOLEM-EV | 173 | 39.9 | 40.9 | 63.1 | 5.74
GOLEM-NV | 168 | 42.2 | 39.9 | 63.1 | 5.74
dotears | 299 | 34.4 | 28.4 | 52.2 | 8.45
DirectLiNGAM | 613 | 31.6 | 25.0 | 46.5 | 15.4
sortnregress | 572 | 34.3 | 26.4 | 50.5 | 15.7
DCDI-G | 2039 | 30.0 | 4.46 | 32.6 | 44.0
GIES | 3038 | 29.2 | 15.4 | 39.2 | 64.5
IGSP | 3064 | 29.0 | 15.2 | 38.9 | 64.6
UT-IGSP | 3075 | 29.8 | 15.1 | 39.6 | 66.0
Table 6: Validation results for all methods at a threshold of $|w|<0.05$. Columns in order are: Method name (Method), number of edges called at threshold (# Called Edge), percent of called edges showing differential expression in either the causal or anticausal direction (% DE), percent of called edges with protein-protein interactions at any confidence level (% PPI), percent of edges called either differentially expressed or have protein-protein interactions (% DE | PPI), and total recall, counting both differential expression and protein-protein interactions as true positives (% Recall). Method | Threshold | # Called Edges | % Causal DE | % Anticausal DE
---|---|---|---|---
dotears | .05 | 299 | 19.7 | 18.7
dotears | .1 | 92 | 22.8 | 21.7
dotears | .2 | 32 | 31.3 | 31.3
NO TEARS | .05 | 149 | 21.5 | 22.1
NO TEARS | .1 | 69.0 | 26.1 | 23.2
NO TEARS | .2 | 20.0 | 40.0 | 30.0
sortnregress | .05 | 572 | 19.6 | 18
sortnregress | 0.1 | 131 | 22.9 | 20.6
sortnregress | 0.2 | 31 | 25.8 | 16.1
DirectLiNGAM | .05 | 613 | 17.9 | 16.8
DirectLiNGAM | .1 | 121 | 23.1 | 20.7
DirectLiNGAM | 0.2 | 22 | 22.7 | 22.3
GOLEM-EV | .05 | 173 | 20.8 | 23.7
GOLEM-EV | .1 | 71 | 22.5 | 25.4
GOLEM-EV | .2 | 18 | 27.8 | 33.3
GOLEM-NV | .05 | 168 | 19.0 | 30.0
GOLEM-NV | .1 | 72 | 20.8 | 26.4
GOLEM-NV | .2 | 21 | 28.6 | 28.6
GIES |  | 3038 | 15.3 | 15.6
IGSP |  | 3064 | 14.8 | 16.2
UT-IGSP |  | 3075 | 16.0 | 15.8
DCDI-G |  | 2982 | 18.7 | 13.7
Table 7: All methods call an equal proportion of causal and anticausal edges.
Columns in order are: Method name (Method), number of edges called at
threshold (# Called Edge), proportion of called edges that show differential
expression in the correct causal direction (% Causal DE), proportion of called
edges that show differential expression in the anticausal direction (%
Anticausal DE).

## Code Availability

Code to reproduce the experiments in this paper, as well as a working
implementation of dotears, is available at https://github.com/asxue/dotears/.

## Acknowledgments

AX was supported by the NIH Training Grant in Genomic Analysis and
Interpretation T32HG002536. This work was partially funded by HHMI Hanna H
Gray and Sloan fellows programs to HP. This work was partially funded by NIH
grant (R35GM125055) and NSF grants (CAREER-1943497, IIS-2106908) to SS. We
thank Nathan LaPierre for helpful conversations.