ARTICLE pubs.acs.org/JPCC
Extending the Density Functional Tight Binding Method to Carbon Under Extreme Conditions Nir Goldman* and Laurence E. Fried Physical and Life Sciences Directorate, Lawrence Livermore National Laboratory, Livermore, California 94550, United States ABSTRACT:
We report herein on simulations of carbon under pressures up to 2000 GPa and 30 000 K using the density functional tight binding method (DFTB) with a parameter set we have specifically designed for these conditions. The DFTB method can provide a high throughput simulation capability compared to KohnSham density functional theory while retaining most of its accuracy. We fit the DFTB repulsive energy to measured and computed diamond isothermal compression data and show that this yields accurate compression curves for diamond, graphite, and the BC8 phase, as well as material properties for all three phases. We then show that our new repulsive energy yields predictions of the Hugoniot of diamond shock compressed to the conducting liquid that are within the range of different experimental measurements. 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.
’ INTRODUCTION The equation of state (EOS) of materials under extreme pressures and temperatures is of great importance for understanding planetary interiors and the chemical reactivity that occurs under strong dynamic compression. Knowledge of the high pressuretemperature 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 Diamond anvil cell experiments have successfully accessed high pressure, low temperature states of matter,4 as well as the lower pressure, high temperature melting lines of compressed materials.5 Thermodynamic states that have been inaccessible with diamond anvil cells have traditionally been achieved through shock compression. Shock compression dynamically strains the sample in one spatial dimension while simultaneously heating the sample.6 Laser-driven shock compression experiments have elucidated many high pressure properties of carbon including the transition to a high pressure conducting liquid79 and possible observation of exotic high pressure phases such as BC8 carbon.10 (The BC8 phase of carbon is a tetrahedrally coordinated solid with eight atoms per unit cell that has been predicted to exist above 1075 GPa.11). However, experiments tend to rely on EOS models for interpretation, which have been shown to be inaccurate for r 2011 American Chemical Society
some systems12,13 and often neglect chemical kinetics that can affect experimental measurements.13 Molecular dynamics (MD) simulations provide an independent route to EOS and chemical reactivity determination, where material properties such as pressures, temperatures, and the shock Hugoniot states are readily computed. Empirical bond-order potentials have been used to conduct MD simulations of solid diamond under shock compression14 and the diamond melt line up to 400 GPa and 12 000 K.15 However, accurate modeling of the breaking and forming of chemical bonds and the thermal electronic excitations that can occur at higher pressures and temperatures usually require the use of quantum theories such as density functional theory (DFT), e.g., ref 16. DFT has been shown to accurately reproduce the high pressuretemperature phase boundaries5,17,18 and shock Hugoniot properties of many materials.1922 DFT-MD Special Issue: Chemistry and Materials Science at High Pressures Symposium Received: July 15, 2011 Revised: September 29, 2011 Published: November 16, 2011 2198
dx.doi.org/10.1021/jp206768x | J. Phys. Chem. C 2012, 116, 2198–2204
The Journal of Physical Chemistry C simulations, though, scale poorly with computational effort and thus are generally limited to system sizes of hundreds of atoms and time-scales of tens of picoseconds.2125 By contrast, laser driven dynamic compression experiments span nanosecond time scales.7,9,10,13,26 As a result, DFT-MD simulations can 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 (with self-consistent atomic charges) 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.27 The DFTB method uses a minimal atomcentered basis set and an approximate Hamiltonian based on KohnSham DFT.28 Thermal populations of excited electronic states are computed through the Mermin functional.29 DFTB thus has the potential to achieve time scales relevant to experiments while providing an accurate picture of chemical reactivity under those conditions. In particular, the DFTB total energy expression contains a pairwise repulsive energy VIJrep (discussed below) that can be fit to experimental and DFT computed data in order to extend the technique to a wide variety of conditions.30,31 DFTB has been shown to yield accurate results for energetic materials at relatively low pressures and temperatures (P < 200 GPa, T < 4000 K).20,3234 However, to date there has been little validation of DFTB at the extreme thermodynamic conditions that pertain to shock compression experiments on carbon (P > 500 GPa, T > 12000 K). To this end, we report on the determination of the DFTB repulsive energy for carbon at pressures up to 2000 GPa and temperatures up to 30 000 K. We first fit the repulsive energy to diamond cold curve data from both experiment and theory. We compute the density and bulk modulus of diamond from our fit and compare to experiment and previous calculations. We then test the transferability of our repulsive energy by computing cold curve data to high pressure and material properties for graphite and the BC8 phase of carbon. Our results show that our fit yields significantly improved agreement with experiment and converged DFT calculations. Finally, in order to examine high-temperature conditions with our new repulsive energy, we conducted MD simulations of diamond shock compressed to the conducting liquid. We observe that our new fit yields pressures, temperatures, and densities that lie within the range of experimental measurements.
’ COMPUTATIONAL DETAILS AND METHODS We now briefly introduce the DFTB method and motivate our fitting of VIJrep. The DFTB Hamiltonian with self-consistent charges (SCC)28 is derived by assuming a spherical charge distribution n(r) and expanding the KohnSham DFT expression for the system energy E[n] to second order in charge fluctuation δn(r) (see ref 35 for a review of different tight binding charge equilibration methods). 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 ionion repulsive interactions as well as exchange-correlation interactions. Erep can be expressed as a sum of pairwise interactions Erep = ∑I < dn ðRc RIJ Þn , IJ Vrep ¼ n¼2 > :0
∑
ðRIJ < Rc Þ otherwise
Here, dn corresponds to the polynomial coefficients and Rc is the cutoff radius. The repulsive energy is thus guaranteed to be zero at R > Rc and to have smooth behavior at Rc. VIJrep can be expressed by any number of functional forms, although for the purpose of this work we choose to retain the above polynomial expression. VIJrep is usually obtained by fitting to high-level DFT of dimer interactions and other small systems.30 DFTB calculations were driven by the LAMMPS molecular software simulation suite,36 with the DFTB+ code27 used to compute forces and the cell stress tensor. All of our calculations only include s- and p-orbital interactions. SCC charge equilibration was utilized due to the possibility of pressure-driven ionization under the extreme thermodynamic conditions studied here. Geometry optimizations were all performed using a force convergence criteria of 106 au or smaller. The bulk moduli B0 for diamond and BC8 carbon were determined by using finite differences to compute the elastic constants with a fractional volumetric strain of 5 106 and then applying the relation B0 = (C11 + 2C12)/3. Convergence of our results was confirmed by performing additional calculations with strains from 106 to 104. Values of B0 for graphite were computed by fitting the results of a series of constant pressure optimizations (i. e., all degrees of freedom including cell lattice vectors were optimized at a specified pressure) to the Vinet EOS.37 DFT calculations were performed using the VASP planewave basis set DFT code.3840 Fractional electron occupancies were calculated using the Mermin functional29 with the smearing set to 0.2 eV, although calculations were found to be converged with lower smearing values. For our cold curve calculations, the wave function convergence criteria was set to 105 eV and the force convergence set to 102 eV/Å. All of our calculations were performed using spin polarization with projector-augmented wave (PAW) pseudopotentials41,42 and the PerdewBurke Ernzerhof generalized gradient approximation (GGA) functional.43 System sizes of 64 atoms were used for diamond and graphite, and 128 atoms for BC8 carbon. Sampling of the Brillouin zone for diamond, graphite, and the BC8 phase were performed with a MonkhorstPack grid44 of up to 2 2 2. Previous calculations under similar thermodynamic conditions have shown the pressures and energies of 64 atom system sizes are converged to within 1% with sampling of the Γ-point, only.70 Our bulk moduli from DFT reported below were estimated by fitting cold curve results to the Vinet EOS, similar to previous work.45 Although previous DFT data exist for diamond, graphite and BC8,4550 here we recompute their cold curves due to the lack of published tabulated data and in order to provide a consistent level of theory for the calculated results. Shock compression simulations were conducted with the multiscale shock technique (MSST).20,21,5153 MSST is a simulation methodology based on the NavierStokes equations for compressible flow. Instead of simulating a shock wave within a large computational cell with many atoms, the MSST computational cell follows a Lagrangian point through the shock wave. This is accomplished by time-evolving equations of motion for the atoms and volume of the computational cell to constrain the stress in the propagation direction to the Rayleigh line and the energy of the system to the Hugoniot energy condition.51,53 (The Hugoniot is 2199
dx.doi.org/10.1021/jp206768x |J. Phys. Chem. C 2012, 116, 2198–2204
The Journal of Physical Chemistry C
ARTICLE
Table 1. Diamond PressureVolume Data Points Used for a the VIJ ref Fit volume (Å3/atom)
P (GPa)
DFTB tuned
3.323
13.193
14.787
3.225 3.129
28.470 45.106
28.071 43.501
3.036
63.181
61.116
2.946
82.774
81.056
2.858
103.971
103.541
2.773
126.861
128.720
2.732
138.938
142.408
2.312
349.567
344.158
2.180 1.993
438.713 625.360
443.749 623.997
1.873
767.057
768.294
1.819
850.923
848.954
1.764
933.280
936.072
1.710
1032.055
1030.420
a
Pressures below 140 GPa were computed from the experimental EOS,63 and those above were computed from DFT.
Table 2. Repulsion Energy Parameters for Standard and Tuned DFTBa
a
parameter
DFTB std.
DFTB tuned
d2
0.0180339
0.0646869
d3
0.0947251
0.0439292
d4 d5
0.128899 0.116864
0.137595 0.0462101
d6
0.301533
0.339320
d7
0.172412
0.202761
d8
0.032959
0.0118606
d9
0.0
0.0262189
Rc
3.83847
3.486079
All values are displayed in atomic units.
the locus of thermodynamic endstates achieved by a specific shock velocity and given set of initial thermodynamic conditions.) For a given shock speed, these two relations describe a steady planar shock wave within continuum theory. MSST thus enables simulation of the shock wave with small system sizes,20 making it possible to simulate with DFT or other computationally intensive force fields such as DFTB. MSST has been used in conjunction with DFT to accurately reproduce the shock Hugoniot of a number of systems.2125 MSST has also been shown to accurately reproduce the sequence of thermodynamic states throughout the reaction zone of shock compressed explosives with analytical equations of state.53 New MSST equations of motion used in this study allow for a self-consistent dynamic electron temperature, where the ionic and electronic temperatures are kept equal every time step.54 We have used the extended Lagrangian BornOppenheimer molecular dynamics (XL-BOMD) approach for propagation of the electronic degrees of freedom5560 in our MD simulations. XL-BOMD uses an extended Lagrangian framework in order to create Verlet-like equations of motion for the initial guess for the charges on the system for each MD time step. Here, we have employed a fifth-order dissipative Verlet-like integrator.57 The
Figure 1. Diamond cold compression curve. The “+” symbols correspond to the experimental EOS,63 open squares correspond to our computed DFT results, open circles correspond to results from standard DFTB, and closed circles correspond to results from our tuned DFTB parameter set.
approximate time-reversibility of the integrator improves the stability of the simulation, thus allowing for incomplete convergence of the SCC charges. Our MSST results reported below all used a maximum of four SCC steps per MD time step, although simulations with as few as two SCC steps were found to allow for stable conserved quantities in our simulations.
’ RESULTS AND DISCUSSION Fitting of the repulsive energy was performed using a Levenberg Marquardt nonlinear least-squares fitting routine.61 We use a previous fitting of the repulsive energy62 from carboncarbon interaction parameters available for download (pbc-0-3 parameter set from http://www.dftb.org) as the starting point for our fit. The pressure-volume data points used in our fit are shown in Table 1, and our new repulsive energy parameters are listed in Table 2. Pressure-volume data points were computed up to 140 GPa using the experimentally determined EOS63 (the high pressure limit of the experiments). We included in our fit DFT data points up to pressures of 1032.1 GPa. For the purpose of this work, our fitted repulsive energy is called the tuned parameter set (DFTB tuned) and the repulsive energy that we used as the starting point for our fit is denoted as the standard parameter set (DFTB std.) The computed diamond cold curve (Figure 1) indicates that the tuned parameter set accurately reproduces both experiment and DFT cold curve points. We observe an overall softening of the diamond cold curve at higher pressure as the result of our fitting procedure. The next-nearest-neighbor distance in diamond under ambient conditions is ∼2.52 Å, and the standard parameter set has a relatively long cutoff radius of 2.02 Å. As a result, starting at a density of 6.70 g/cm3 (V = 2.97 Å3/atom), next-nearest-neighbor interactions will be included in the repulsive energy, which will cause the resulting EOS to stiffen further at pressures greater than ∼1300 GPa. By contrast, our tuned parameter set has a much shorter cutoff radius of 1.84 Å, which excludes next-nearest-neighbor interactions until a density of 8.95 g/cm3 (V = 2.23 Å3/atom) and a pressure of ∼2200 GPa are reached. 2200
dx.doi.org/10.1021/jp206768x |J. Phys. Chem. C 2012, 116, 2198–2204
The Journal of Physical Chemistry C
ARTICLE
Table 3. Density, Lattice Parameters, and Bulk Moduli for Diamond, the BC8 Phase and Graphitea phase
F0
a
c
(g/cm3)
(Å)
(Å)
3.49 3.53
3.58 3.56
352.5 527.2
DFT
3.48
3.58
421.0
expt.63
3.51
3.57
446
DFTB tuned
3.37
4.56
343.8
DFTB std.
3.50
4.50
360.7
DFT
3.50
4.50
DFTB tuned
1.74
2.48
9.14
3.69
0.7
DFTB std. DFT50
1.66 1.71
2.48 2.47
9.04 8.84
3.65 3.58
9.8 1.0
expt.66
2.26
2.462 6.707 2.724
31
method
diamond DFTB tuned DFTB std.
BC8
graphite
a
bulk modulus c/a
(GPa)
347.0
DFT results for graphite are taken from ref 50.
Material properties of diamond are reported in Table 3. Our computed density for the tuned parameter set at P = 0 (F0) is 3.49 g/cm3, which agrees within 1% with experiment63 and our computed DFT result. The bulk modulus from the tuned parameter set (352.5 GPa) is lower than the reported experimental value (446 GPa).63 Our DFT estimate of the bulk modulus (421 GPa) compares well with results from computed elastic constants (430 GPa),64 and is significantly higher than a previous fit of the Vinet EOS model to DFT data (368.2 GPa) from Correa et al.45 Given that we have fit the repulsive energy over a wide range of pressures, we do not expect our results to necessarily be optimal right in the neighborhood of V0. Rather, our repulsive energy provides an accurate representation of diamond over a broad range of volumes. Similar to diamond, we observe that our tuned parameter set provides improved agreement with the cold compression curve for the BC8 phase computed from DFT (Figure 2a). Values of B0 and F0 are shown in Table 3. The value of F0 from the tuned parameter set (3.37 g/cm3) is ∼4% lower than that from DFT (3.50 g/cm3), while the bulk modulus (343.8 GPa) agrees to within 1% with our DFT estimate (347.0 GPa). Our DFT result for F0 agrees within ∼1% with previous DFT calculations,49 although it is somewhat lower than calculations using the local density approximation (LDA) (3.65 g/cm346 and 3.70 g/cm347) as well the Correa EOS (3.92 g/cm3).45 However, the improved accuracy due to use of GGA techniques have been known for some time.65 Similarly, the computed bulk modulus for BC8 carbon from the LDA calculations (411 GPa46 and 454 GPa47) is somewhat higher than our result, as is that from Correa et al. (539.7 GPa).45 For the cold compression curve for graphite (Figure 2b), we show the compression in the basal plane only due to the fact that our calculations did not include any dispersion corrections, which are needed in order to accurately model the interplanar van der Waals interactions.50 We observe that in addition to yielding a cold curve that is too stiff at lower pressure, the standard parameter set exhibits a sudden drop in pressure at a relative volume of ∼0.81 that is absent in both the DFT and tuned parameter set results. We attribute this pressure drop to a buckling of the graphitic planes when bonding between planes begins to occur. DFTB calculations with both parameter sets and our DFT results all yield values for F0 that are lower and the c lattice parameters that are higher than
Figure 2. BC8 phase (a) and graphite basal plane (b) cold curve compression results. The thick solid line corresponds to results from DFT, open circles correspond to standard DFTB, and closed circles correspond to the tuned DFTB parameter set. Cold curves are plot as a function of relative volume (V/V0) due to differing values of V0. The thin dashed line in both panels is drawn as a guide for the standard DFTB results.
experimental values (Table 3). The F0 (1.74 g/cm3) value from the tuned parameter set shows slightly closer agreement with that from DFT (1.71 g/cm3),50 compared to the result from the standard parameter set (1.66 g/cm3). By contrast, we observe close agreement between both parameter sets, DFT and experiment for the a lattice parameter. The value of B0 from the tuned parameter set (0.7 GPa) agrees more closely with that from DFT50 (1.4 GPa), although the experimental result is reported between 31 and 34 GPa.66,67 These discrepancies with experiment are due to the strong underestimation of dispersion forces between graphitic planes in these calculations. Comparison of the computed relative energies of graphite, diamond, and BC8 phase (Table 4) shows that the tuned parameter set yields the correct energetic ordering of the phases. 2201
dx.doi.org/10.1021/jp206768x |J. Phys. Chem. C 2012, 116, 2198–2204
The Journal of Physical Chemistry C
ARTICLE
Table 4. Relative Energies of Graphite, Diamond and BC8 Carbon (eV/atom) phase
DFTB tuned
DFTB std.
DFT
graphite
0.0
0.0
0.0
diamond BC8
0.11 0.70
0.046 0.69
0.0120.01247,48,50 0.72
Figure 4. Shock velocity versus particle velocity Hugoniot results. Symbols are the same as defined in Figue 3.
Figure 5. Pressure versus temperature Hugoniot results. The dashed line corresponds to results from the Correa EOS,45 the “+” symbols correspond to DFT,70 the solid line correspond to experiment (with error bars),26 the closed black circles correspond to our tuned DFTB parameter set, and the closed red circles correspond to standard DFTB.
Figure 3. Pressure versus density Hugoniot results for both standard and tuned parameter sets, only (a) and compared to experimental results (b). The open red circles and green triangles correspond to experimental results from Hicks et al.9 and Brygoo et al.,69 respectively, using a quartz standard, and the blue squares correspond to experimental results from Nagao8 using an aluminum standard. The closed black circles correspond to results from the tuned DFTB parameter set and the closed red circles from standard DFTB. The blue “+” symbols in (b) correspond to results from DFT70 and the dashed line from the Correa EOS.45
However, the predicted energetic difference between graphite and diamond of 0.11 eV/atom is higher than the experimentally
measured value of 0.026 eV/atom.68 The standard parameter set yields a diamond energy that is 0.046 eV/atom lower than that of graphite, and several different DFT calculations yield results between 0.01248 to 0.01250 eV/atom. Both standard and tuned parameter sets yield an energy difference of ∼0.7 eV/atom between graphite and the BC8 phase, in agreement with our DFT result and previous DFT calculations.64 Finally, we present results of MD simulations of diamond shock compressed to the conducting liquid. For the pressure versus density Hugoniot (Figure 3), our tuned parameter set provides Hugoniot pressures that are lower by 200300 GPa for a given density (Figure 3a) than standard DFTB. The tuned results show improved agreement with experimental results using an aluminum standard8 (Figure 3b). However, results from both parameter sets fall within the rather broad range of the experimental uncertainty.8,69,9 Both parameter sets also yield 2202
dx.doi.org/10.1021/jp206768x |J. Phys. Chem. C 2012, 116, 2198–2204
The Journal of Physical Chemistry C Hugoniot pressures that are higher than previously computed results from DFT70 and the Correa EOS.45 The high coordination numbers (4.56.4)70 from DFT at these Hugoniot points could indicate that d-orbitals play a role in the bonding that occurs. By contrast, both the standard and tuned parameter sets only contain s- and p-orbital bonding interactions. We have also computed the shock velocity Us versus particle velocity Up relation from our simulations, using the relation Up = Us(1 F1/F2), where F1 corresponds to the density of the preshocked state and F2 the density of the shocked state. Results show that both parameter sets yield results within the experimental uncertainty (Figure 4), with the tuned parameter yielding slightly improved agreement with experiment. The temperature versus pressure Hugoniot from the tuned parameter set (Figure 5) yields good agreement with recent experimental measurements26 and results from DFT in the region of overlap. However, we note that the DFT Hugoniot temperatures extrapolate to higher values than DFTB and experiment at pressures greater than ∼1400 GPa, in rough agreement with the Correa EOS.45 This is consistent with softer Hugoniot pressures from DFT, which again could be the result of excluding d-orbital interactions from the DFTB simulations.
’ CONCLUSIONS We have shown that by fitting the repulsive energy to the diamond cold curve, DFTB can be extended to the extreme pressures and temperatures of laser-driven dynamic compression experiments. Our new repulsive energy has been used to make accurate predictions of material properties of diamond, graphite, and the BC8 phase over pressures up to 2000 GPa, including initial densities, bulk moduli, and compression curves. Agreement with DFT and experimental Hugoniot data for liquid carbon could possibly be improved by including d-orbital interactions in the DFTB parameter sets, which is the subject of future work. The approach we have used for carbon could be extended to compute the EOS and kinetics of any number of materials related to geology and planetary science, including silicon, SiO2, and hydrocarbon systems. The increased computational efficiency of DFTB could yield accurate EOS and chemical data of these materials at extreme thermodynamic conditions on experimental time scales for the first time. ’ AUTHOR INFORMATION Corresponding Author
*E-mail:
[email protected].
’ ACKNOWLEDGMENT This work was performed under the auspices of the U.S. Department of Energy by the Lawrence Livermore National Laboratory (LLNL) under Contract DE-AC52-07NA27344. Computations were performed at LLNL using the Aztec massively parallel computer. We thank Alfredo Correa for tabulated results for the EOS-determined PF Hugoniot curve shown in Figure 3. ’ REFERENCES (1) Hubbard, W. B. Science 1981, 214, 145. (2) Segretain, L.; Chabrier, G.; Hernanz, M.; Garcia-Berro, E.; Isern, J.; Mochkovitch, R. Astrophys. J. 1994, 434, 641. (3) Seager, S.; Kuchner, M.; Hier-Majumder, C. A.; Militzer, B. Astrophys. J. 2007, 669, 1279.
ARTICLE
(4) Goncharov, A. F.; Struzkhin, V. V.; Somayazulu, M. S.; Hemley, R. J.; Mao, H. K. Science 1996, 273, 218. (5) Goncharov, A. F.; Goldman, N.; Fried, L. E.; Crowhurst, J. C.; Kuo, I.-F. W.; Mundy, C. J.; Zaug, J. M. Phys. Rev. Lett. 2005, 94, 125508. (6) Zel’dovitch, Y. B.; Raizer, Y. P. Physics of Shock Waves and HighTemperature Hydrodynamic Phenomena; Dover Publications: Mineola, NY, 2002. (7) Bradley, D. K.; Eggert, J. H.; Hicks, D.; Celliers, P.; Moon, S. J.; Cauble, R. C.; Collins, G. Phys. Rev. Lett. 2004, 93, 195506. (8) Nagao, H.; et al. Phys. Plasmas 2006, 13, 052705. (9) Hicks, D. G.; Boehly, T. R.; Celliers, P. M.; Bradley, D. K.; Eggert, J. H.; McWilliams, R. S.; Jeanloz, R.; Collins, G. W. Phys. Rev. Lett. 2008, 78, 174102. (10) Knudson, M. D.; Desjarlais, M. P.; Dolan, D. H. Science 2008, 322, 1822. (11) Correa, A. A.; Bonev, S. A.; Galli, G. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 1204. (12) Lyzenga, G. A.; Ahrens, T. J.; Nellis, W. J.; Mitchell, A. C. J. Chem. Phys. 1982, 76, 6282. (13) Barrios, M. A.; Hicks, D. G.; Boehly, T. R.; Fratanduono, D. E.; Eggert, J. H.; Celliers, P. M.; Collins, G. W.; Meyerhofer, D. D. Phys. Plasmas 2010, 17, 056307. (14) Zybin, S. V.; Elert, M. L.; White, C. T. Phys. Rev. B 2002, 66, 220102. (15) Ghiringhelli, L. M.; Los, J. H.; Meijer, E. J.; Fasolino, A.; Frenkel, D. Phys. Rev. Lett. 2005, 94, 145701. (16) Gygi, F.; Galli, G. Phys. Rev. B 2002, 65, 220102. (17) Goldman, N.; Fried, L. E.; Kuo, I.-F. W.; Mundy, C. J. Phys. Rev. Lett. 2005, 94, 217801. (18) Schwegler, E.; Sharma, M.; Gygi, F.; Galli, G. Proc. Natl. Acad. Sci. U.S.A. 2008, 105, 14779. (19) Kress, J. D.; Mazevet, S.; Collins, L. A.; Wood, W. W. Phys. Rev. B 2000, 63, 024203. (20) Reed, E. J.; Manaa, M. R.; Fried, L. E.; Glaesemann, K. R.; Joannopoulos, J. D. Nat. Phys. 2008, 4, 72–76. (21) Mundy, C. J.; Curioni, A.; Goldman, N.; Kuo, I.-F.; Reed, E.; Fried, L. E.; Ianuzzi, M. J. Chem. Phys. 2008, 128, 184701. (22) Goldman, N.; Reed, E. J.; Kuo, I.-F. W.; Fried, L. E.; Mundy, C. J.; Curioni, A. J. Chem. Phys. 2009, 130, 124517. (23) Goldman, N.; Reed, E. J.; Fried, L. E. J. Chem. Phys. 2009, 131, 204103. (24) Wu, C.; Fried, L. E.; Yang, L. H.; Goldman, N.; Bastea, S. Nat. Chem. 2009, 1, 57. (25) Goldman, N.; Reed, E. J.; Fried, L. E.; Kuo, I.-F. W.; Maiti, A. Nat. Chem. 2010, 2, 949. (26) Eggert, J. H.; Hicks, D. G.; Celliers, P. M.; Bradley, D. K.; McWilliams, R. S.; Jeanloz, R.; Miller, J. E.; Boehly, T. R.; Collins, G. Nat. Phys. 2010, 6, 40. (27) Aradi, B.; Hourahine, B.; Frauenheim, T. J. Phys. Chem. A 2007, 111, 5678(http://www.dftb-plus.info). (28) Elstner, M.; Porezag, D.; Jungnickel, G.; Elsner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Phys. Rev. B 1998, 58, 7260. (29) Mermin, N. D. Phys. Rev. 1965, 137, 1441. (30) Porezag, D.; Frauenheim, T.; K€ohler, T.; Seifert, G.; Kaschner, R. K. Phys. Rev. B 1995, 51, 12947. (31) Koskinen, P.; M€akinen, V. Comput. Mater. Sci. 2009, 47, 237. (32) Manaa, M. R.; Fried, L. E.; Melius, C. F.; Elstner, M.; Frauenheim, T. J. Phys. Chem. A 2002, 106, 9024. (33) Margetis, D.; Kaxiras, E.; Elstner, M.; Frauenheim, T.; Manaa, M. R. J. Chem. Phys. 2002, 117, 788. (34) Manaa, M. R.; Reed, E. J.; Fried, L. E.; Galli, G.; Gygi, F. J. Chem. Phys. 2004, 120, 10146. (35) Gaus, M.; Cui, Q.; Elstner, M. J. Chem. Theory Comput. 2011, 7, 931. (36) Plimpton, S. J. Comput. Phys. 1995, 117, 1(http://lammps. sandia.gov). (37) Vinet, P.; Smith, J. R.; Ferrante, J.; Rose, J. H. Phys. Rev. B 1987, 35, 1945. 2203
dx.doi.org/10.1021/jp206768x |J. Phys. Chem. C 2012, 116, 2198–2204
The Journal of Physical Chemistry C
ARTICLE
(38) Kresse, G.; Hafner, J. Phys. Rev. B 1993, 47, R558. (39) Kresse, G.; Hafner, J. Phys. Rev. B 1994, 49, 14251. (40) Kresse, G.; Furthm€uller, J. Phys. Rev. B 1996, 54, 11169. (41) Bl€ochl, P. E. Phys. Rev. B 1994, 50, 17953. (42) Kresse, G.; Joubert, D. Phys. Rev. B 1999, 59, 1758. (43) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1997, 78, 1396–1396. (44) Monkhorst, H. J.; Pack, J. D. Phys. Rev. B 1976, 13, 5188. (45) Correa, A. A.; Benedict, L. X.; Schwegler, D. A. Y. E.; Bonev, S. Phys. Rev. B 2008, 78, 024101. (46) Fahy, S.; Louie, S. G. Phys. Rev. B 1987, 36, 3373. (47) Mailhiot, C.; McMahan, A. K. Phys. Rev. B 1991, 44, 11578. (48) Ribeiro, F. J.; Tangney, P.; Louie, S. G.; Cohen, M. L. Phys. Rev. B 2005, 72, 214109. (49) Sun, J.; Klug, D. D.; Martonak, R. J. Chem. Phys. 2009, 130, 194512. ngyan, J. G. J. Phys. Chem. A (50) Bucko, T.; Hafner, J.; Lebegue, S.; A 2010, 114, 11814. (51) Reed, E. J.; Fried, L. E.; Joannopoulos, J. D. Phys. Rev. Lett. 2003, 90, 235503. (52) Reed, E.; Fried, L. E.; Manaa, M. R.; Joannopoulos, J. D. A multi-scale approach to molecular dynamics simulations of shock waves. In Chemistry at Exteme Conditions; Manaa, M., Ed.; Elsevier: Amsterdam, 2005; pp 297326. (53) Reed, E. J.; Fried, L. E.; Henshaw, W. D.; Tarver, C. M. Phys. Rev. E 2006, 74, 056706. (54) Reed, E. J. to be published in this issue. (55) Niklasson, A. M. N.; Tymczak, C. J.; Challacombe, M. Phys. Rev. Lett. 2006, 97, 123001. (56) Niklasson, A. M. N. Phys. Rev. Lett. 2008, 100, 123004. (57) Niklasson, A. M. N.; Steneteg, P.; Odell, A.; Bock, N.; Challacombe, M.; Tymczak, C. J.; Holmstr€om, E.; Zheng, G.; Weber, V. J. Chem. Phys. 2009, 130, 214109. (58) Odell, A.; Delin, A.; Johansson, B.; Bock, N.; Challacombe, M.; Niklasson, A. M. N. J. Chem. Phys. 2010, 131, 244106. (59) Zheng, G.; Niklasson, A. M. N.; Karplus, M. J. Chem. Phys. 2011, 135, 044122. (60) Sanville, E. J.; Bock, N.; Niklasson, A. M. N.; Cawkwell, M. J.; Sewell, T. D.; Dattelbaum, D. M.; Sheffield, S. A. Proceedings of the 14th International Symposium on Detonation 2011, in press. (61) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes in Fortran 77; Cambridge University Press: Cambridge, 1992. (62) K€ohler, C.; Frauenheim, T. Surf. Sci. 2006, 600, 453. (63) Occelli, F.; Loubeyre, P.; Letoullec, R. Nat. Mater. 2003, 2, 151. (64) Wen, B.; Takami, S.; Kawazoe, Y.; Adschiri, T. Physica B 2011, 406, 2654. (65) Perdew, J. P.; Chevary, J. A.; Vosko, S. H.; Jackson, K. A.; Pederson, M. R.; Singh, D. J.; Fiolhais, C. Phys. Rev. B 1992, 46, 6671. (66) Zhao, Y. X.; Spain, I. L. Phys. Rev. B 1989, 40, 994. (67) Hanfland, M.; Beister, H.; Syassen, K. Phys. Rev. B 1989, 39, 12598. (68) Fahy, S.; Wang, X. W.; Louie, S. G. Phys. Rev. B 1990, 42, 3503. (69) Brygoo, S.; Henry, E.; Loubeyre, P.; Eggert, J.; Koenig, M.; Loupias, B.; Benuzzi-Mounaix, A.; Gloahec, M. R. L. Nat. Mater. 2007, 6, 274. (70) Romero, N. A.; Mattson, W. D. Phys. Rev. B 2007, 76, 214113.
2204
dx.doi.org/10.1021/jp206768x |J. Phys. Chem. C 2012, 116, 2198–2204
