Energetic, vibrational, and electronic properties of silicon using a nonorthogonal tight-binding model

,


I. INTRODUCTION
As the capabilities of materials simulations increase, so does the demand for methodologies that can capture the important physics accurately while being fast enough to simulate large systems for long periods of time.[4][5][6][7][8][9][10] The main challenge in developing these models has been the determination of the Hamiltonian matrix elements ͑and in the case of nonorthogonal bases, the overlap matrix elements͒ as a function of interatomic distance.The most common approach is to fit the results of either total-energy calculations or band-structure features to experiment or ab initio calculations, assuming a particular functional form with some free parameters for the distance dependence of the matrix elements.Most of the models in the literature, however, suffer from certain shortcomings.Many models assign a large part of the total energy of the system to a repulsive pair potential to compensate for adopting an orthogonal set of basis orbitals. 2,5,6Some models do not give an accurate description of the band structure of the ground state, 7 or the energetics of important features of bulk semiconducting systems, such as point defects. 3,4Some models were intended for use for a specific application, and were optimized and tested only for geometries relevant to that application. 9Finally, other models use a very large number of free parameters and fit a very large number of structures, leading to potentially unphysical values for the parameters and the suspicion that the good fit to the data set may not guarantee reliability. 8 the past, TB models have mainly been used as interpolative schemes assuming that the configurations that were being fit will encompass the configurational space where the model will be used.In this paper we use the ͑nonorhogonal͒ NRL-TB method, 11 an extrapolative method that uses parameters obtained from fitting ab initio calculations of a few high-symmetry structures, to compute the energies of a wide range of geometries for silicon.We find that the results compare very well to ab initio calculations for configurations that are substantially different from those included in the fitting data set.This increases our confidence that the reason for the accurate results in the tested configurations is that the physics underlying the model is correct.
This paper is organized as follows: In Sec.II we describe the functional form of our TB parametrization and the fitting data set.In Sec.III we discuss applications of the TB model to bulk properties such as the diamond lattice band structure and energies of other lattices, point defect properties such as formation and relaxation energies, and surface energies and reconstructions.We also report on two applications of the TB model that would be impractical with ab initio methods, one using molecular dynamics ͑MD͒ simulations to compute finite temperature properties of the bulk crystal, and another using the model to compute electronic properties of large amorphous systems.In the final section we summarize the results.

II. FUNCTIONAL FORM AND FITTING
In this paper we present results for two parametrizations, one using a sp 3 basis, which has already been presented in some detail, 10 and another using a sp 3 d 5 basis.Since their functional forms are nearly identical and have already been presented, we will only give a brief summary here.The total energy of the system is written as the sum of the energies of the occupied electronic eigenstates.The onsite Hamiltonian matrix elements vary with the local density, allowing the NRL-TB method be fit to linearized augmented plane wave ͑LAPW͒ eigenvalues that have been shifted so that the LAPW total energy is equal to the eigenvalue sum.Therefore all of the contributions to the total energy are accounted for in the eigenvalue sum and the addition of a repulsive pair potential, a feature common to most TB models, is not needed.
The energies of the electronic states and the corresponding eigenvectors are the solutions to a generalized eigenvalue equation with Hamiltonian and overlap matrix elements parametrized as follows: The basis used to describe the Hamiltonian and overlap matrices is a set of one s, three p, and for one of the parametrizations five d orbitals around each atom, with all interactions assumed to be in the two-center approximation. 12A local atomic density at atom i is defined as where R ជ i is the position of atom i and is a fitting param- eter.The cutoff function f (R) is given by

͑2͒
where R c is 12.5 a.u. and L c is 0.5 a.u.The onsite matrix elements are given in terms of the local atomic density i as where l is the orbital type index (s, p, or d) and ␣ l , ␤ l , ␥ l , and l are fitting parameters.The distance dependent parts of the two-center Hamiltonian matrix elements are given by where l and lЈ are orbital type indices, is an index for the type of interaction between orbitals (, , or ␦), and the parameters a ll Ј b ll Ј , c ll Ј , and g ll Ј are fitting parameters.The distance dependent parts of the overlap matrix elements are where t ll Ј , q ll Ј , r ll Ј , and u ll Ј are fitting parameters, and ␦ ll Ј is the Kronecker delta.Note that the symbols for some of the parameters are different from those used in Ref. 10.The overlap matrix elements have similar functional form to the Hamiltonian matrix elements, but are constrained to go to the correct value, zero or one, at zero interatomic separation.The angular dependence of the Hamiltonian and overlap matrix elements is the standard two-center Slater-Koster form. 12he 41 parameters used by the functional form for the sp 3 basis parametrization are fit to four high-symmetry crystal lattices: simple cubic ͑sc͒, face-centered cubic ͑fcc͒, bodycentered cubic ͑bcc͒, and the diamond structure.The fitting data set includes both the total energy and band structure of each lattice, as computed by LAPW ab initio densityfunctional theory ͑DFT͒ calculations in the local-density approximation ͑LDA͒ for a wide range of volumes around the energy minimum.The diamond lattice data included the widest range of volumes, from 12.2 Å 3 /atom to 22.7 Å 3 /atom.The sc lattice structures ranged from 12.6 Å 3 /atom to 18.5 Å 3 /atom, the fcc lattice from 12.7 Å 3 /atom to 15.0 Å 3 /atom, and the bcc lattice from 13.0 Å 3 /atom to 16.0 Å 3 /atom.The best-fit root-meansquare ͑rms͒ error of the valence-band energies for the diamond lattice is 0.12 eV, and the rms error for the crystal lattice total energies is 0.020 eV.
The sp 3 d 5 basis parametrization has 69 parameters, which were fit to the diamond lattice band structure at volumes ranging from 13.5 Å 3 /atom to 22.7 Å 3 /atom.This set of parameters does not allow for any Hamiltonian or overlap matrix elements between different d orbitals, but allows all interactions between s and p orbitals, as well as s Ϫd and pϪd interactions.Since this model is optimized for accuracy in the band structure, we adjusted the DFT/LDA calculations, which predict an indirect gap of 0.5 eV, by applying a uniform shift of 0.67 eV to the conduction bands of the ideal volume diamond lattice, matching the band gap to the experimental result. 13,14For the other lattice constants in the fit we shifted the conduction band so that the gap was scaled up by a factor of 1.17/0.5ϭ2.34.The shift amount increased monotonically from 0.67 eV for larger volumes and decreased monotonically for smaller volumes.The best fit rms error for the diamond lattice valence band and lower two conduction bands is 0.21 eV, and the rms error for the diamond lattice total energies is 0.004 eV.
The parameters that result from this fit for the sp 3 and sp 3 d 5 basis models are listed in Tables I and II, respectively.The sp 3 basis Hamiltonian and overlap matrix elements are plotted in Fig. 1, and the onsite matrix elements for the dia- mond structure are plotted in Fig. 2. The variation of the onsite matrix elements with nearest-neighbor distance is structure dependent because they have a nonlinear dependence on the density, which is itself a structure-dependent quantity.Note that the ss, pp, and pp Hamiltonian and overlap parameters have the expected sign, while the sp parameters are opposite in sign to the usual convention. 12owever, this sign is not physically meaningful, since it is determined by the ͑arbitrary͒ choice of sign of the s and p basis orbitals, and does not affect the eigenvalues or energies computed with the model.Cohen et al. have also used a similar method to generate parameters for silicon. 9They used a somewhat different functional form and concentrated their fitting and tests on high-pressure phases.Since we are interested in applying this TB model to complex structures, including defects and surfaces, but at or near atmospheric pressure, we have employed a different set of geometries for our fitting data set.

III. APPLICATIONS
To test the transferability of the sp 3 parameters we computed the total energy of a range of structures important for condensed phases of silicon including bulk systems, point defects, and surfaces.First we review the diamond lattice band structure, cohesive energies of a number of bulk lattices as a function of volume, and the elastic constants of the diamond structure, as presented in Ref. 10.To address the issue of improving the diamond lattice band structure we present an electronic structure calculation using the sp 3 d 5 parametrization.We expand our analysis of the sp 3 parameters to include phonon spectra at several high-symmetry points.As a more stringent test we use this parametrization to compute the energies of some lower-symmetry configurations.For the bulk we simulate important point defects, including several low-energy interstitial configurations, the vacancy, and the concerted exchange pathway for diffusion.For the (100) and (111) surfaces we compute ideal surface energies as well as relaxation energies for a number of reconstructions.From MD simulations we compute the meansquared atomic displacement for a range of temperatures, the vibrational density-of-states, and the phonon-dispersion relations.Finally, we use the sp 3 d 5 basis model to study the electronic structure of large amorphous systems.

A. Bulk
The band structure of the diamond structure lattice, which was part of the fitting data set, is shown in Fig. 3.The valence band is in very good agreement with LAPW calculations.The conduction band is not as well described, with the  minimum indirect gap of 1.02 eV appearing at the L point rather than at approximately three fourths of the way from the ⌫ to the X point as first-principles calculations and experiments have established.The size of the gap is somewhat smaller than the experimental value of 1.17 eV, 15 although it is larger than the DFT/LDA prediction to which it was fitted.
We have addressed the issue of obtaining a better fit of the gap and the conduction band and came to the following conclusions.The addition of five d orbitals to the basis improves the fit of the conduction band, as can be seen in Fig. 3.The lowest-energy conduction band is nearly perfect when compared to a LAPW calculation with a rigid 0.67 eV shift of the conduction bands, 13,14 and the higher bands are also closer to the LAPW calculation than with the sp 3 basis model.The valence bands are very well described, although the lowest band at the ⌫ point is too flat.The density-ofstates ͑DOS͒, including its decomposition into contributions from different angular momentum states ͑which is found using TB eigenvectors that were not fitted͒, is also in very good agreement with DFT/LDA calculations.The three peaks in the valence band are clear, as is the decomposition into mainly s character in the lowest peak, mixed s and p character in the middle peak, and p character in the third peak.There is very little d character in the valence band, while the conduction band is mainly of mixed p and d character, with smaller s contribution.
To obtain such a good fit for the band structure the sp 3 d 5 model was fit only to the full diamond lattice band structure at all volumes.The lack of other structures and energy information in the fitting data set deteriorates the energetics of the model.We decided that the best compromise is to use the minimal sp 3 basis for all total-energy calculations presented in this paper, and to use the sp 3 d 5 to compute the electronic structure of amorphous silicon presented later in this section.We note in passing that in his book, Papaconstantopoulos 16 was able to obtain a good fit of the conduction band near the gap with a sp 3 basis model.However, that work differs from the present approach in two important ways: first, it utilizes three-center integrals and second, it treats the Hamiltonian and overlap matrix elements for the first three neighbor shells as independent parameters, rather than giving them as an analytical functional form that varies with distance.These differences provide the flexibility that produces a better fit to the conduction states.
The total energies as a function of volume for a wide range of structures are shown in Fig. 4, and their equilibrium structural and energy properties are listed in Table III.All of the structures have higher energy than the diamond structure, including some low-energy, rarely examined theoretical phases such as hexagonal diamond and the clathrate structures. 17The TB model reproduces the LAPW results very well for the four structures to which it was fit, as can be seen from the equilibrium energies, volumes, and bulk moduli listed in Table IV.These quantities were computed using a Birch fit 18 to the fitting data and the sp 3 TB model calculations.The elastic constants of the ground-state diamond structure are listed in Table V.Those that do not involve shear, c 11 and c 12 , are within 22% of LAPW calculations and 14% of experiment. 19The shear elastic constant computed without allowing for relaxation of the internal degrees of freedom of the two-atom unit cell c 44 * , is 34% larger than the LAPW result.Allowing for the internal relaxation brings c 44 within 19% of experiment. 19A detailed compari- son of the energetic and structural parameters of two clathrate structures, Si 34 and Si 46 , comparing results of our TB model with experiment, DFT calculations, and results of an orthogonal TB model, is shown in Table VI.The energy differences between the clathrates and the diamond structure are lower than in DFT calculations, although they are of the correct sign.The structures, including both the lattice constant and the internal structural parameters of the basis, are within 1% of the experimental values. 20,21This agreement is as good as that provided by DFT/LDA plane-wave calculations 22 or by an ab initio localized orbital method, 17 and substantially better than the orthogonal TB model of Goodwin et al. 2 tested by Kahn and Lu. 23honon frequencies at high-symmetry points in the Brillouin zone ͑BZ͒ computed with the TB model using the frozen phonon approximation are compared with experimentally measured values in Table VII. 24The agreement is quite good: the TB results are within 15% of experiment for all but three of the modes, the X 3 , L 1 , and W 2 modes, which are off by about 25%, 30%, and 60%, respectively.While this good description of the phonon spectra is a nontrivial test of the model, it represents only infinitesimal deviations of atomic positions from the diamond structure, which was part of the fitting data set.In the following sections we show that this TB model can also accurately describe the energies of configurations that are substantially different from those in the fitting data set.

B. Point defects
The ground-state structure for silicon is the covalently bonded, open network of the diamond structure.Since strong covalent bonds in the ideal lattice allow for little atomic motion at temperatures below the melting point, processes such as diffusion are dominated by point defects, which are far more mobile than perfectly bonded atoms. 25The formation energy of such defects strongly influences their concentrations and is therefore an important material property.In Equilibrium energies and structural features of hypothetical crystal lattices for Si computed with the sp 3 TB model.E is equilibrium energy relative to the diamond structure in eV per atom, V is the volume in Å 3 per atom, and c/a is the unit cell aspect ratio.The internal structural parameter x for the H-Dia structure is the position ( 13 , 2 3 ,x) of the atom at site 4 f , and for the BC8 structure is the position (x,x,x) of the atom at site 16c.Notation for the lattice types is the same as in Fig. 4 Lattice  46 .⌬E is given in eV/atom, the lattice constant a 0 is given in Å , and the internal parameters are given in terms of the lattice constant.Using the notation of Ref. 17, the parameters of the Si 34 structure are the position (x e ,x e ,x e ) of the atom at site 32a, and the position (x g ,x g ,z g ) of the atom at site 96g.The parameters of the Si 46 structure are the position (x i ,x i ,x i ) of the atom at site 16i and the position (0,y k ,z k ) of the atom at site 24k.PW denotes plane-wave basis DFT/LDA calculations from Ref. 22 and of the vacancy.All defect energies were computed using a 16.29 Å cube cell with 216Ϯ1 atoms, and sampling the BZ at the ⌫ point.To compute the relaxed configurations a conjugate-gradient algorithm was used, 26 with the atomic positions relaxed until the force on each atom was less than 3 meV/Å .In agreement with ab initio calculations the ͗110͘ split is the lowest-energy interstitial configuration, followed by the tetrahedral and hexagonal configurations. 27,288][29] The formation energy of the ideal vacancy is also accurately predicted, although its relaxation energy is about twice as large as DFT/LDA calculations predict. 27,30he relaxed geometries are in approximate agreement with ab initio results, 27,28,30 although there are some differ-ences.The tetrahedral interstitial reduces its symmetry, with three of the four atoms that surround it relaxing outward and the fourth relaxing inward, while the interstitial atom itself moves parallel to the fourth atom.A diagram of this configuration, shown in Fig. 5, makes clear that the relaxed interstitial atom has moved part way ͑0.34Å͒ towards a hexagonal site without significantly stretching (Ͻ2.5%) any of its bonds.This geometry differs from the DFT/LDA work where the relaxed interstitial had the full symmetry of the initial configuration, 27 although whether those calculations were restricted to maintain tetrahedral symmetry or simply found no energy gain from breaking the symmetry is not stated.The hexagonal interstitial retains its ideal hexagonal symmetry with the six ring atoms moving outward by 0.1 Å, in perfect agreement with the DFT/LDA results. 27The four atoms around vacancy relax inward by about 0.28 Å , in good agreement with DFT/LDA calculations, which find an inward relaxation of 0.25 Å .However the structure keeps the full tetrahedral symmetry rather than reducing to the symmetry of the tetragonal structure that ab initio calculations predict. 27,30o test the accuracy of the TB model in describing the breaking of bonds within a relatively undisturbed solid we computed the energy for the concerted exchange pathway for diffusion. 31,32The energy along the path for the ideal and relaxed configurations, as well as DFT/LDA results, 31,32 are plotted in Fig. 6.The ideal energy is within 22% of the DFT/LDA calculation throughout the path, with no spurious minima.The relaxation energy is substantially too high, although the relaxed length of the bond between the diffusing atoms of 2.15 Å is nearly identical to the DFT/LDA result.
Point defect configurations include substantial deviations from the ideal lattice geometry and several inequivalent atomic sites.In such a situation it is possible for charge to be transfered between atoms.If this charge transfer is substantial, the applicability of a model with no Coulombic interaction or charge self-consistency may be in doubt.The variation of the onsite energies in our TB model could potentially exacerbate this effect.To measure the amount of charge transfer we performed a Mulliken population analysis on the  vacancy, the tetrahedral interstitial, and the hexagonal interstitial.The deviation of the charge from neutrality ͑four electrons per atom͒ was modest, less than half an electron in every case, and as low or lower for the relaxed configurations as for the ideal ones.By comparison, a nonorthogonal TB model 7 with constant onsite energies predicted somewhat smaller charge transfers ͑by 30% to 50%͒, except for the ideal tetrahedral interstitial where it predicted Ϫ0.65e as compared with Ϫ0.35e on the interstitial atom.This comparison indicates that the variation of the onsite elements does not substantially increase the charge transfer, which are neglected by most TB models.

C. Surfaces
The properties of the silicon surface are critical for processes such as surface growth and etching.Since atoms on a surface have an asymmetric environment and lower coordination than in the bulk, the resulting configurations are qualitatively different from the types of geometries that this TB parametrization was fitted to.Therefore, surface structures provide an important test of the transferability of the param-eters.We compute the surface energies of the ideal ͑100͒ and ͑111͒ surfaces, as well as the relaxation energies for some simple but representative reconstructions of these surfaces.For both surfaces we use a symmetric slab configuration ͓24 layers for the ͑100͒ surface, 6 bilayers for the ͑111͒ surface͔ and a set of k points equivalent to a 4ϫ4 mesh in the full planar BZ of the 1ϫ1 surface unit cell.For the ͑100͒ surface we examined the 2ϫ1 buckled dimer reconstruction, and for the ͑111͒ surface we examined 2ϫ2 reconstructions with adatoms at three inequivalent sites, the T 4 , H 3 , and B 2 .The results of these calculations are listed in Table IX.
The energetics of the ideal surfaces are in agreement with DFT/LDA calculations, 33,34 with the ͑111͒ surface about 0.8 eV lower in energy than the ͑100͒ surface.Both the relaxation energy and the structure of the buckled dimer reconstruction of the ͑100͒ surface are in good agreement with DFT/LDA calculations. 33,35,36One of the three adatom reconstructions of the ͑111͒ surface, with the adatom at the T 4 site, is related to the dominant features of the ground-state 7ϫ7 reconstruction. 37The H 3 adatom site is a metastable position that is involved in the diffusion of adatoms.The energy of the B 2 site, which lies halfway along the path connecting the T 4 and H 3 sites, determines the barrier for migration between them.
The TB model predicts the same ordering in energy as DFT/LDA calculations for the three adatom sites, 34,38,39 an improvement over previous nonorthogonal TB models. 7The agreement is not quantitative, however, as the TB model overestimates the binding energy of the adatom at the T 4 site by about 50%.In this position the adatom forms a bond of length of 2.34 Å to the atom underneath, and bonds of length 2.42 Å to the three next-nearest-neighbors.DFT/LDA calculations predict corresponding bond lengths of 2.43 Å and 2.47 Å . 34The surface atom that does not form bonds to the adatom, called the rest atom, makes three bonds of length 2.38 Å with angles of 100.5°, in very good agreement with the DFT/LDA result of 2.34 Å long bonds with angles of 99.9°.The binding energy of the adatom in the H 3 site is somewhat underestimated as compared with DFT/LDA calculations. 34,39In this site the adatom makes three 2.39 Å

D. Finite temperature properties
To examine some finite temperature properties of the diamond structure crystal we used the NRL TB-MD moleculardynamics package 40 developed by Kirchhoff to evolve a 512atom unit cell at constant energy for 2000 moleculardynamics time steps ͑each step corresponds to 2.0 fs͒, varying the initial kinetic energy of the atoms.From the resulting positions we computed the mean-squared displacement.This quantity, plotted against the temperature measured in the sample, is shown in Fig. 7.A linear fit of the mean-squared displacement gives a slope of 1.72 ϫ10 5 Å 2 /K.In comparison, the Debye temperature, extracted from experimental measurements of the temperature dependence of x-ray diffraction peak broadening, 41 corresponds to a mean-squared atomic displacement temperature coefficient of 1.75ϫ10 5 Å 2 /K. 42rom the Fourier transforms of velocity autocorrelation functions calculated during MD simulation we obtained the vibrational density-of-states. 43The phonon-dispersion curves were extracted from the Fourier transform of the velocity and position dependent autocorrelation function for a given polarization and k vector. 43The vibrational densities-of-states from two MD simulations, one at 300 K and one at 1500 K, are plotted in Fig. 8.The overall shapes are quite similar, although the high-temperature curve is smoother.The peaks are shifted to lower frequencies at 1500 K, indicating a softening of the vibrational modes.The corresponding phonondispersion relations are plotted in Fig. 8, and compared against experimental values extracted from the graphs of Ref. 44.The dispersion relations from the 300 K MD simulation are in good agreement with experiment.The only significant errors are an underestimate of the frequencies of the second branch between the ⌫ and L points ͓⌳ 1 (A)͔ and an overestimate of the frequencies of the highest two branches near the X point ͓⌬ 5 (O), ⌬ 2 Ј(O), ⌺ 2 (O), and ⌺ 1 (O)].Some, although not all, of the vibrational frequencies in the 1500 K MD simulation are lower than at 300 K, a result that is consistent with the differences between the vibrational densities-of-states in Fig. 8.The high-frequency branches have the largest shifts, with many of them shifting down by 0.5 THz, about 5%.The low-frequency branches, with the exception of the ⌳ direction, also shift down, with a particularly prominent shift of the second branch at the K point ͓⌺ 3 (A)͔.

E. Amorphous structures
As an illustration of the expanded modeling capabilities offered by the present TB parametrizations for Si, we study the electronic properties of bulk and surface structures in amorphous Si ͑a-Si͒.We emphasize at the outset that the following discussion is concerned mostly with demonstrating the capabilities of the approach, rather than the physics of a-Si, the latter being a broader problem beyond the scope of the present paper.
Based on many experimental and theoretical studies, 45,46 it is widely accepted that a-Si has the basic structure of a continuous random network ͑CRN͒ of tetrahedrally bonded atoms, [47][48][49] but the question of defects has been the subject of considerable debate in recent years 50,51 and remains controversial. 52Experimental results appear to favor undercoordinated, three-fold bonded atoms as the dominant defects ͑so-called ''dangling bonds''͒.On the other hand, theoretical simulations, using a wide variety of ab initio, 52,53 semiempirical [54][55][56] and empirical [57][58][59][60][61][62][63] methods, consistently produce both undercoordinated as well as overcoordinated ͑fivefold bonded͒ defects, with a significant preference for the latter.The type of bonding arrangements at the surface of a-Si is even less clear than in the case of bulk defects, since surface-specific measurements are not readily available.It has been reported that ab initio relaxation of a bulkterminated CRN model produces a surface with roughly equal numbers of threefold and fivefold bonded atoms. 644,65 In order to study the electronic properties of bulk and surface a-Si, we first prepare a 32.9 Å ϫ65.8 Å ϫ16.5 Å bulk amorphous sample with 1728 atoms.Because such a large sample is impractical to generate using any electronic structure based simulation method, we simulate the quenching of the liquid with an interatomic potential, 66 following a procedure similar to the one used in Ref. 62.The resulting structure has over 96% tetrahedral coordination, with only fivefold coordinated defects.
To model the surfaces of a-Si, we considered two qualitatively different 1728 atom slabs.The ''cleaved sample'' is created from the bulk structure by turning off the periodic boundary conditions in the third direction.The ''quenched sample'' is created directly from the liquid phase by turning off periodic boundary conditions in the third direction and quenching the resulting liquid slab with the interatomic potential.Not surprisingly, the quenched surfaces are slightly rougher ͑by about 1 Å͒ than the cleaved surfaces.The cleaved surfaces regions contain mostly ͑65%͒ fourfold coordinated atoms, with a predominance of threefold ͑29%͒ over fivefold ͑6%͒ coordinated atoms.On the other hand, the quenched surfaces have somewhat higher fourfold coordination ͑72%͒, with many ͑27%͒ fivefold and almost no ͑1%͒ threefold coordinated atoms.Note that these surfaces are fully reconstructed, and therefore ''bulk'' concepts of defects do not apply; many of the fourfold coordinated atoms do not have tetrahedral bond angles, and the fivefold coordinated atoms tend to appear in clusters near the top layer of the surface, rather than as isolated floating bonds below the top layer.
While the size of these samples makes calculating their electronic structure with first-principles methods impractical, the sp 3 d 5 basis TB parametrization described in Sec.II makes such a study feasible.In Fig. 9 we compare the electronic DOS of the three amorphous samples computed with the TB model ͑see also Fig. 3 for the diamond structure crystal electronic DOS͒.The bulk amorphous sample DOS is much smoother than the crystal, in agreement with other simulation results. 53,67There are only two peaks in the valence band, corresponding to the highest and lowest peaks of the DOS of the crystal.Despite the structural differences between the two surface models, their overall DOS is quite similar, showing that the electronic signatures of undercoordinated and overcoordinated atoms at the amorphous surface are difficult to distinguish.This result is consistent with the arguments of Pantelides for bulk defects. 50The DOS of the surface samples differ from that of the bulk sample essentially only in the gap region.This indicates that all surfacerelated defects produce gap states, consistent with analogous results for bulk defects. 52,54,65It is also interesting that the surface DOS above the gap region is depleted relative to the bulk one.A more detailed analysis of these results will be given elsewhere. 68Here, we wish to point out only the effi-ciency of the approach in generating reliable electronic structure information for large systems.

IV. SUMMARY
We have applied the NRL-TB method to generate TB models for Si that were fit to LAPW results of a small number of high-symmetry crystal structures.We found that the resulting Hamiltonians are transferable to a much wider range of geometries.A model with a nonorthogonal sp 3 basis reproduces DFT/LDA and experimental measurements for a wide range of material properties, including elastic constants and phonon frequencies, point defect formation energies, and surface energies and reconstructions.In fact, this TB model is as good or better at describing the energetics of point defects than some models that included such structures in their fitting process.It is also the only nonorthogonal TB model we are aware of that correctly describes the energy sequence of different adatom configurations on the ͑111͒ surface of silicon.The ability of the model to accurately describe such diverse systems, despite having been fitted to a small number of high-symmetry crystal structures, increases our confidence that the model captures the essential physics of bonding in solid-state silicon systems.
The efficiency of the model makes it possible to study finite temperature properties of large silicon systems through molecular-dynamics simulation, an application that would be impractical with DFT/LDA methods.We have shown that the model reproduces experimental results for atomic meansquared displacements as a function of temperature in bulk silicon.Phonon densities-of-states and dispersion curves extracted from MD simulations at different temperatures show good agreement with experiment at low temperatures, and a substantial softening of many modes at higher temperatures.By adding d orbitals and modifying the fitting data set, we obtained a model that accurately reproduces both the valence and conduction bands of silicon in the diamond structure, at the price of deterioration in the accuracy of the energetics.This sp 3 d 5 parametrization makes possible the study of the electronic structure of amorphous systems with nearly 2000 atoms.

FIG. 1 .
FIG. 1. Hamiltonian matrix elements ͑upper panel͒ and overlap matrix elements ͑lower panel͒ for the sp 3 parametrization plotted as a function of interatomic distance.

FIG. 3 .
FIG. 3. Band structure of Si in the diamond lattice computed with the sp 3 ͑upper panel͒ and sp 3 d 5 ͑middle panel͒ bases, and the density-of-states for the sp 3 d 5 model ͑bottom panel͒.Dashed lines are TB results, solid lines are DFT/LDA results ͑a rigid shift has been applied to the DFT/LDA conduction-band results used for the sp 3 d 5 basis model͒.All energies are referred to as the valence-band maximum.

FIG. 5 .
FIG. 5. Structure of the tetrahedral interstitial in ideal ͑white spheres͒ and relaxed ͑black spheres͒ configurations, with bonds connecting atoms within 2.5 Å .The geometry is viewed along the ͗110͘ direction with the ͗11 ¯1͘ direction pointing up.

FIG. 6 .
FIG. 6.Energy along the concerted exchange pathway for the diffusion of atoms without vacancies or interstitials.DFT/LDA calculations of the ideal configuration from Ref. 31 are plotted in a solid line, TB calculations of the ideal configurations in a dashed line, and TB calculations of the relaxed configurations in a dot-dash line.The relaxed DFT/LDA value is only available for the saddle point.

FIG. 7 .
FIG.7.Mean-squared displacement of atoms in the diamond lattice as a function of temperature.The points are computed from MD simulations, the solid line is a linear fit going through the origin, and the dashed line is a line with a slope corresponding to the experimental measurement as discussed in the text.

FIG. 9 .
FIG. 9.The electronic TB DOS computed with the sp 3 d 5 parametrization for the bulk and surface a-Si samples.The lower curve ͑dashed line͒ is the DOS for the bulk amorphous sample.The two upper curves show the excess DOS associated with the surfaces, plotted as the difference between the cleaved sample DOS and bulk amorphous DOS ͑dotted line͒, and the difference between the quenched sample DOS and the bulk amorphous DOS ͑solid line͒.

TABLE I .
Parameters for the sp 3 basis tight-binding model.

TABLE II .
Parameters for the sp 3 d 5 basis tight-binding model.

TABLE IV .
Equilibrium energy relative to the diamond structure ͑E͒ in eV per atom, volume ͑V͒ in Å 3 per atom, and bulk modulus ͑B͒ in GPa, for the lattice structures in the fitting data set, computed with the sp 3 TB model and with LAPW DFT/LDA.

TABLE V .
Elastic constants in GPa for the diamond structure, computed with the sp 3 TB model, LAPW DFT/LDA calculations, and experiment.c 44 * is the shear elastic modulus computed without allowing for relaxation of the internal degrees of freedom in the two-atom unit cell.

TABLE VI .
Energy differences ⌬E relative to the diamond structure and structural parameters for the two optimized clathrate structures, Si 34 and Si , LO denotes local orbital DFT/LDA calculations from Ref. 17, OTB denotes the orthogonal tight-binding results from Ref. 23, and NRL-TB denotes the present paper.Note that the experimental samples are actually of Na x Si 34 with xϽ11, and Na 8 Si 46 .

TABLE IX .
39rface energies ␥ ͑in eV per 1ϫ1 surface unit cell͒, for the (100) and (111) surfaces of Si, relaxation energies ⌬␥ ͑in eV per 1ϫ1 surface unit cell͒ relative to the ideal surface, and selected structural features, computed using the TB model and compared with DFT/LDA results.PRB 62bonds with its nearest neighbors.No DFT/LDA results for adatom-surface bond lengths are available.The rest atom makes three 2.40 Å long bonds with angles of 101.1°, as compared with 2.34 Å long bonds with angles of 104.9°inDFT/LDA calculations.In the B 2 site the TB model adatom makes two bonds of length 2.33 Å with the surface.The binding energy in this configuration is underestimated as compared to the DFT/LDA result.39 a From Refs.33-36 and 39.