The melting temperature of bulk silicon from ab initio molecular dynamics simulations

doi:10.1016/j.cplett.2009.09.075

Chemical Physics Letters

Volume 481, Issues 1–3, 19 October 2009, Pages 88-90

https://doi.org/10.1016/j.cplett.2009.09.075 Get rights and content

Abstract

We estimated a melting temperature of T_m ∼ 1540 ± 50 K at zero pressure for silicon from constant enthalpy and constant pressure (NPH) Born–Oppenheimer Molecular Dynamics (BOMD) simulations of a coexisting crystalline–liquid phase. The computed T_m is below the experimental melting point of 1685 K, but it is consistent with a previously predicted first-order liquid–liquid phase transition (LLPT) at a critical point T_c ∼ 1232 K and P_c ∼ −12 kB [P. Ganesh, M. Widom, Phys. Rev. Lett. 102 (2009) 075701], which is in a highly supercooled state.

Graphical abstract

The melting temperature of crystalline silicon calculated from ab initio molecular dynamics simulations at the constant enthalpy.

Introduction

Bulk silicon and water both belong to the so-called tetrahedral liquid family [1], and thus both share certain structural similarities. For example, both liquid silicon and water expand upon freezing, and both assume diamond cubic structures in the solid state. Moreover, when deeply supercooled, both liquid silicon and water may undergo a first-order high-density liquid (HDL) to low-density liquid (LDL) phase transition [2], [3], [4], [5], [6], [7], [8]. The HDL and LDL phases are metastable compared to the bulk crystalline phase since the temperature of either HDL or LDL is below the melting point. Indeed, the liquid–liquid phase transition (LLPT) of silicon at the supercooled state has been previously studied by both classical molecular dynamics (MD) [4], [5] and ab initio MD simulations [6], [7], [8]. These computational studies of supercooled silicon may also shed light on the physical behavior of possible LLPTs in other tetrahedral liquids. The melting temperature T_m is an important thermodynamic property that is needed to obtain information about the metastable supercooled (T < T_m) and the stable liquid (T > T_m) regions of the phase diagram. The melting temperature of bulk silicon at ambient conditions has been previously computed from classical MD simulations [9], [10], [11] and from an ab initio MD simulation [12]. In the latter study, a thermodynamic integration of the free energy method was used to compute the melting temperature of silicon [12]. It will be of interest to compare the ab initio MD simulation of the coexisting system with the thermodynamic integration method for computing the melting temperature of silicon. Here, we present the results of calculating T_m from a constant enthalpy and constant pressure (NPH) Born–Oppenheimer molecular dynamics (BOMD) simulation. In our approach the energies and forces are obtained from electronic structure calculations with the Perdew–Burke–Ernzerhof (PBE) exchange-correlation density functional [13].

Section snippets

Computational method

Two commonly used computational approaches have been previously introduced to obtain T_m: (1) the direct simulation of the solid–liquid coexistence [10], [11], [14], [15], [16], [17], [18], [19], [20] system and (2) the thermodynamic integration (TI) of the free energy of the solid and liquid phases [12], [21], [22]. In the second method, T_m is determined by the condition of equality of the Gibbs free energies of the liquid and solid [12], [21], [22], viz. $G_{liq} (P, T)_{T = T_{m}} = G_{solid} (P, T)_{T = T_{m}}$ . Using

Results and discussion

The lower and upper limits of the initial temperature in the simulations were chosen with the following criteria: (i) the lower limit of the initial temperature of T = 1500 K was chosen below the melting temperature of bulk silicon (measured value T_m ∼ 1685 K) and (ii) the upper limit of the initial temperature of T = 1800 K is chosen higher than T_m ∼ 1691 K of silicon with the SW potential [9], [10], [11]. As mentioned above, the melting temperature of DFT-GGA-based Si was predicted to be T_m = 1492 ± 50 K

Acknowledgements

We acknowledge support from the Chemical Sciences, Geosciences and Biosciences Division, and the Materials Science and Engineering Division (DE-FG02-04ER46164), Office of Basic Energy Sciences, U.S. Department of Energy, and from the Nebraska Research Initiative. Battelle operates the Pacific Northwest National Laboratory for the U.S. Department of Energy. This research was performed in part using the Molecular Science Computing Facility (MSCF) in the Environmental Molecular Sciences

References (29)

C.J. Wu et al.
Phys. Rev. Lett.
(2002)
C.A. Angell et al.
Phys. Chem. Chem. Phys.
(2000)
O. Mishima et al.
Nature
(1998)
P.H. Poole et al.
Science
(1997)
S. Sastry et al.
Nat. Mater.
(2003)
P. Beaucage et al.
J. Phys.: Condens. Matter
(2005)
N. Jakse et al.
Phys. Rev. Lett.
(2007)
T. Morishita
Phys. Rev. Lett.
(2006)
P. Ganesh et al.
Phys. Rev. Lett.
(2009)
J.Q. Broughton et al.
Phys. Rev. B
(1987)

U. Landman et al.

Phys. Rev. Lett.

(1986)

S. Yoo et al.

J. Chem. Phys.

(2004)

O. Sugino et al.

Phys. Rev. Lett.

(1995)

D. Alfè et al.

Phys. Rev. B

(2003)

J.P. Perdew et al.

Phys. Rev. Lett.

(1996)

Cited by (22)

Determination of the accuracy and reliability of molecular dynamics simulations in estimating the melting point of iron: Roles of interaction potentials and initial system configurations
2019, Journal of Molecular Liquids
Citation Excerpt :
These techniques are quite rigorous but somewhat difficult as several simulations may be required to compute and map Gibbs free energy as a function of temperature and pressure. Alternatively, another set of techniques referred to as ‘direct methods’ involve the direct dynamic simulation of the melting process based on approaches such as the hysteresis method, the voids method and solid-liquid interface-based methods [37–39]. The most direct route to simulate the melting process using MD simulations involves the gradual elevation of simulation temperature until a perfect crystal breaks down and melts; however melting points are generally overestimated due to the hysteresis effects [40,41].
Molecular dynamics (MD) simulations were carried out to investigate the melting transition of iron with an aim to determine the accuracy and reliability of simulated melting points (MP), computing efficiency, and the role played by initial configurations. Two simulation techniques with equilibrium and non-equilibrium heating, and two well-known potentials: Mendelev (EAM) and Etesami (MEAM) were used to investigate a variety of initial configurations containing perfect crystal solid lattice and solid-liquid (S-L) coexistence regions with a number of S-L interfaces and relative S-L proportions (PC, SLSS-2, SLSS-4, SLSS-2 (20%L) and SLSS-2(80%L)). The simulated magnitude of MPs using the non-equilibrium heating method was found to be ~2020 K, which is about 200 K higher than experimental MP of iron (1811 K). Similar values were obtained for both potentials, PC and SLSS-2 configurations; there was a strong evidence for hysteresis effects as well. Simulations with equilibrium heating technique showed that MPs were found to be identical for the wide range of S-L configurations investigated: Mendelev potential (1769 ± 1 K) and Etesami potential (1811 ± 1 K). Corresponding results from the purely solid PC configurations were found to be much higher (1981 ± 1 K and 1951 ± 1 K respectively). Results on other system properties including pair distribution functions, density, local melting and solidification are also reported. These studies have helped identify optimal simulation parameters, potentials and computational approaches to probe the melting region of iron towards applications to problems of industrial importance.
Affect of the graphene layers on the melting temperature of silicon by molecular dynamics simulations
2016, Computational Materials Science
Citation Excerpt :
Melting temperature of silicon in this case is Tm,pure = 1540 K observed by the curve in Fig. 6. This result is in an excellent agreement with the simulated crystalline silicon melting temperature at atmospheric pressure obtained by the two-phase method (1540 K) [23], and the free energy method (1492 K) [24], simulations of using the Environment-Dependent Interatomic potential (1520 K) [25,26]. By adding the empty intervals at the boundaries of the system along the Z direction, the silicon system has the two free surfaces during the simulation procedure.
In the paper, molecular dynamics simulations have been used to detect the melting temperature of silicon. The two models have been considered as the graphene/silicon/graphene (GSG) sandwich and the only silicon system. Atoms in the models interact with each other via a 3-body Tersoff potential with a modified part based on a Coulomb potential and the Ziegler–Biersack–Littmark universal screening function. We find that C–Si interaction prevents the melting of the two silicon surfaces adjoining the graphene layers. The melting temperature of the silicon part in the GSG sandwich (T_m = 2450 K) is 1.6 times higher than that of the only silicon system (T_m_,_pure = 1540 K). Difference of these melting temperatures has original from interaction of C–Si pairs, which causes decreasing of the energy of the silicon part in the GSG sandwich during heating process leading to an increase of the melting temperature of this silicon part.
Theoretical evaluation of corrosion inhibition performance of three antipyrine compounds
2015, Computational and Theoretical Chemistry
Citation Excerpt :
While the study of corrosion mechanisms is supported by quantum chemistry, it is not enough to only consider the isolated inhibitor molecular, the metal interface must be also considered. Since the 1990s, researchers turned to focus on the study of inhibitor molecular–metal interface system to mitigate the above problems [13–15]. From the beginning of the 21th century to the present, so many researches using molecular dynamics simulation have been emerged to investigate the inhibitor mechanism shifted from micro-scale to meso-scale [15–18].
Corrosion inhibition performance of three antipyrine derivatives namely, 4-amino-antipyrine (AAP), 4-[(4-hydroxy benzylidene)-amino]-antipyrine (SAAP) and 4-[(4-methoxy benzylidene)-amino]-antipyrine (AAAP), for metal surface corrosion was studied by using quantum chemistry, molecular mechanics and molecular dynamics simulation methods. Correlations between the quantum chemical parameters (E_HOMO, E_LUMO, the energy gap, global hardness, global softness and Fukui index) and experimentally determined inhibition efficiency have been carried out and discussed. The adsorption behavior of AAP, SAAP and AAAP on metal surface was studied using molecular dynamics simulation methods, through which the binding energy was obtained. The order is as follows: AAAP > SAAP > AAP, which is well accorded with the result (inhibition efficiencies) obtained from the experiments.
Fluid-mineral Equilibrium Under Nonhydrostatic Stress: Insight From Molecular Dynamics
2024, American Journal of Science
Determination of the phase transition of solutions of lithium salts in sulfolane by the molecular dynamics method
2023, Russian Chemical Bulletin
Modified Born method for modeling melting temperature using ab initio molecular dynamics
2023, Journal of Physics Condensed Matter

View all citing articles on Scopus

View full text

The melting temperature of bulk silicon from ab initio molecular dynamics simulations

Abstract

Graphical abstract

Introduction

Section snippets

Computational method

Results and discussion

Acknowledgements

Phys. Rev. Lett.

Phys. Chem. Chem. Phys.

Nature

Science

Nat. Mater.

J. Phys.: Condens. Matter

Phys. Rev. Lett.

Phys. Rev. Lett.

Phys. Rev. Lett.

Phys. Rev. B

Phys. Rev. Lett.

J. Chem. Phys.

Phys. Rev. Lett.

Phys. Rev. B

Phys. Rev. Lett.