Article pubs.acs.org/JPCC
Determination of a Density Functional Tight Binding Model with an Extended Basis Set and Three-Body Repulsion for Carbon Under Extreme Pressures and Temperatures Nir Goldman,*,† Sriram Goverapet Srinivasan,‡ Sebastien Hamel,† Laurence E. Fried,† Michael Gaus,§ and Marcus Elstner∥ †
Physical and Life Sciences Directorate, Lawrence Livermore National Laboratory, Livermore, California 94550, United States Department of Mechanical and Nuclear Engineering, Pennsylvania State University, State College, Pennsylvania § Department of Chemistry and Theoretical Chemistry Institute, University of Wisconsin, Madison, 1101 University Avenue, Madison, Wisconsin 53706, United States ∥ Institute of Physical Chemistry, Karlsruhe Institute of Technology, Kaiserstrasse 12, 76131 Karlsruhe, Germany ‡
ABSTRACT: We report here on development of a density functional tight binding (DFTB) simulation approach for carbon under extreme pressures and temperatures that includes an expanded basis set and an environmentally dependent repulsive energy. We find that including d-orbital interactions in the DFTB Hamiltonian improves determination of the electronic states at high pressure−temperature conditions, compared to standard DFTB implementations that utilize s- and p-orbitals only for carbon. We then determine a three-body repulsive energy through fitting to diamond, BC8, and simple cubic cold compression curve data, as well pressures from metallic liquid configurations from density functional theory (DFT) simulations. Our new model (DFTB-p3b) yields approximately 2 orders of magnitude increase in computational efficiency over standard DFT while retaining its accuracy for condensed phases of carbon under a wide range of conditions, including the metallic liquid phase at conditions up to 2000 GPa and 30 000 K. Our results provide a straightforward method by which DFTB can be extended to studies of covalently bonded materials under extremely high pressures and temperatures such as the interiors of planets and other large celestial bodies.
■
GPa.12) Experiments tend to rely on equation of state (EOS) models for interpretation, which have been shown to be inaccurate for some systems15,16 and often neglect chemical kinetics that can affect experimental measurements.16 Determination of experimental temperatures generally remains an unresolved issue, partially due to large uncertainties in the calibration of pyrometric measurements.13,14 In addition, carbon has been computed to have a number of high pressure metastable states with high energetic barriers.17 These states could be accessed during dynamic compression experiments, whereas EOS models predict equilibrium states of matter over essentially infinite time scales.18 Molecular dynamics (MD) simulations provide an independent route to time-resolved phase transformations and chemical reactivity determination, where material properties such as compression curves and phase diagrams are readily computed. Empirical potentials have been used successfully to conduct MD simulations of several reactive systems under elevated pressure and temperature.19,20 However, accurate
INTRODUCTION The equation of state of materials under extreme pressures and temperatures (P > 100 GPa, T > 1000 K; 1 GPa = 10 kbar) is of great importance for understanding planetary interiors and the chemical reactivity that occurs under strong dynamic compression. Knowledge of the high pressure−temperature properties of simple compounds such as solid phases of carbon is needed to devise models of Neptune and Uranus,1 brown dwarfs,2 and extra-solar carbon planets.3 In addition, carbon is major component of candidate materials for design of fusion capsules for the National Ignition Facility (NIF).4 Diamond anvil cell experiments have successfully accessed high-pressure, low-temperature states of matter,5 as well as the lower-pressure, high-temperature melting lines of compressed materials.6 Thermodynamic states that have been inaccessible with diamond anvil cells can be achieved through dynamic compression experiments. Laser-driven ramp compression7 and shock compression experiments8−10 have elucidated many high-pressure properties of carbon, including the transition to a conducting liquid and indirect observation of exotic solid phases such as BC8 carbon.11 (The BC8 phase of carbon exhibits a body-centered cubic lattice with eight atoms per unit cell, and has been predicted to exist above 1075 © 2013 American Chemical Society
Received: December 26, 2012 Revised: March 13, 2013 Published: March 20, 2013 7885
dx.doi.org/10.1021/jp312759j | J. Phys. Chem. C 2013, 117, 7885−7894
The Journal of Physical Chemistry C
Article
system pressures and RDFs with converged DFT calculations. Finally, we have conducted MD simulations of diamond under shock compression conditions over pressures up to 2000 GPa and temperatures up to 30 000 K. We observe that DFTB-p3b yields results that compare closely to those from DFT and experiments for both solid diamond and metallic liquid end states.
modeling of the breaking and forming of chemical bonds and transitions from insulating to metallic phases usually requires the use of quantum theories such as density functional theory (DFT), for example, ref 21. DFT has been shown to accurately reproduce the high pressure−temperature phase boundaries6,22,23 and shock Hugoniot properties of many materials.24−27 DFT-MD simulations, however, scale poorly with computational effort and thus are generally limited to system sizes of hundreds of atoms and time-scales of tens of picoseconds.26−30 In contrast, laser driven dynamic compression experiments span nanosecond time scales.7,8,10,11,14,16 DFT-MD simulations can consequently be intractable for many high pressure studies, where chemical kinetic effects can span the entire time-scale of experiments. The density functional tight binding method (DFTB) holds promise as an alternative approach for simulations of materials at high pressures and temperatures. DFTB is an approximate quantum simulation technique that allows for several orders of magnitude increase in computational efficiency while retaining most of the accuracy of standard DFT.31 DFTB methods generally use a minimal atom-centered basis set (e.g., s- and porbitals for carbon, only) and an approximate Hamiltonian based on Kohn−Sham DFT.32 It thus has the potential to achieve time scales relevant to experiments while providing an accurate picture of chemical reactivity under those conditions.33,34 DFTB has been shown to yield accurate results for energetic materials at relatively low conditions (P < 200 GPa, T < 4000 K).28,35−38 We recently determined a pairwise repulsive energy for carbon that yielded accurate results for diamond and the BC8 phase up to ∼2000 GPa.39 However, our model yielded pressures for the metallic liquid phase that were significantly higher than many experimental results and produced radial distribution functions (RDFs) that were overstructured relative to previous DFT results. We speculated that this discrepancy was partially due to the inability of the porbital interactions and pairwise repulsive energy to account for the up to 6-fold coordination of carbon under these hot, compressed conditions. To this end, we report on the determination of a new highpressure DFTB model for carbon that includes both d-orbital interactions, as well as a screened, three-body empirical function for the repulsive energy. We call our new model DFTB-p3b, where “p” stands for polarizable functions (i.e., dorbitals) and “3b” stands for the three-body repulsive energy. We first analyze the effect of including d-orbital interactions on the computed electronic states and compute the many-center effects on the DFT total energy for diamond at high density. We then determine the three-body repulsive energy by fitting to experimental and computed cold curve data for diamond, as well as computed cold curve data for the simple cubic and BC8 phases and pressures from snapshots of metallic liquid configurations. MD simulations under metallic liquid conditions show that DFTB-p3b yields excellent agreement for
■
COMPUTATIONAL DETAILS AND METHODS We now briefly introduce the DFTB method and motivate our determination of the repulsive energy. The DFTB Hamiltonian with self-consistent charges (SCC) is derived by assuming a neutral, spherical charge distribution n0(r) and expanding the Kohn−Sham DFT expression for the system energy E[n] to second order in charge fluctuation δn(r). The resulting equation for the system energy is Etot = EBS[n0] + Ecoul[δn] + Erep(R). EBS is the band structure energy (i.e., Hamiltonian excluding charge transfer), Ecoul is the energy from charge fluctuations, and Erep is the repulsive energy, which contains both ion−ion repulsive interactions as well as exchange− correlation interactions. For standard tight binding approaches, EBS is computed using a two-center approximation, and Erep is taken to be pairwise and short-ranged.32−34 Erep can be expressed by any number of functional forms and has previously been obtained by fitting to converged DFT-GGA calculations of dimer interactions and other small systems.33 For conditions where electron delocalization is significant for such as metallization, the two-center pairwise approximations can yield inaccurate bonding energies and band structures.40,41 Attempts have been made to improve upon this through both explicit calculation of three-center Hamiltonian terms40,42,43 and inclusion of an environmental dependence within the twocenter approximation.41,44,45 The latter has been approximated in tight binding approaches through the use of screening functions,41,45 which are designed to mimic the electronic screening effects in solids. In this case, the interaction energy between two atomic sites approaches zero as a third site approaches the midpoint of the line connecting them. Here, we create an environmentally dependent repulsive energy based on the screening principle where Erep is expressed as a sum of pairwise interactions multiplied by three-body screening terms, Erep = ΣiΣj>iΣkVij(1 − Sijk). Vij is expressed as a ninth-order polynomial ⎧ 9 ⎪ ∑ Cn(R cij − R ij)n Vij = ⎨ n = 2 ⎪ ⎩ 0
, (R ij < R cij) otherwise
The Cn correspond to the polynomial coefficients and Rijc is the pairwise cutoff radius. Sijk is expressed as the multiplication of a sixth-order polynomial and hyperbolic tangent function
⎧ ⎡ 6 ⎡ ⎛ ⎛ R ik + R jk ⎞a3⎞⎤ ⎛ ijk R ik + R jk ⎞m⎤ 1⎢ ⎪ ⎢ ⎟⎟ ⎟⎥ , ⎪ ∑ b ⎜R − ⎟ ⎥tanh⎢a1 exp⎜⎜ −a 2⎜⎜ ⎟ Sijk = ⎨ 3 ⎢⎣ m = 2 m⎝ c R 2 ⎠ ⎥⎦ ⎝ ⎠ ⎠⎥⎦ ij ⎝ ⎣ ⎪ ⎪ ⎩ 0
In this case, bm are the polynomial coefficients and a1, a2, and a3 are also parameters to be fit. Rijk c is the three-body interaction
⎡ (R ik + R jk) ⎤ ⎢ < R cijk ⎥ ⎢⎣ ⎥⎦ 2 otherwise
cutoff radius. Both pairwise and three-body cutoff radii guarantee that Vij and Sijk will have smooth behavior around 7886
dx.doi.org/10.1021/jp312759j | J. Phys. Chem. C 2013, 117, 7885−7894
The Journal of Physical Chemistry C
Article
Rijc and Rijk c and will be equal to zero beyond. The bm coefficients were constrained to yield values of Sijk between 0 and 1, where the function is smoothly varying. As a result, the repulsive energy is completely screened and equal to zero when a third atomic site exists at the midpoint of the line connecting two atomic sites and recovers a pairwise functional form when the third atomic site is beyond Rijk c (Figure 1).
and p-orbitals only). We have analyzed results for both ambient diamond (3.51 g/cm3 and 300 K) and one metallic liquid configuration (7.668 g/cm3 and 15 890 K), taken from a DFTMD simulation (Figure 2). The DOS from our tight-binding calculations are independent of Erep, allowing for a direct comparison to density of states results from DFT without any empirical function fitting. DFT results at ambient conditions show that the occupied states are largely composed of s- and porbitals with d-orbitals making more significant contributions for the unoccupied states (Figure 2a). Correspondingly, we observe that both DFTB-p3b and standard DFTB yield accurate results for the density of occupied states with DFTB-p3b yielding somewhat closer results to DFT for both the band gap and the density of unoccupied states. However, upon forming the metallic liquid we observe significant dorbital character for both occupied and unoccupied states, which increases monotonically with energy (Figrue 2b). This is consistent with results from both our own DFT-MD simulations (discussed below) and previous results,61 where the observed ∼6-fold coordination of the first solvation shell can necessitate incorporating d-orbital interactions into bonding. Our results show that DFTB-p3b yields significantly improved agreement with DFT over standard DFTB for the electron states at extreme thermodynamic conditions. The DFTB-p3b DOS exhibits a large peak at ∼20 eV that is not present in the results from DFT, which we attribute to errors due to the two-center Hamiltonian approximation. However, this peak occurs at sufficiently high energy that these electronic states are inaccessible even at the extreme temperature simulated here. Improving calculation of the electronic DOS over the entire energy spectrum is the subject of future study. The many-center contributions to the DFT total energy can be computed using the equation Enb = EDFT − Σi
