Optimal atomic structure of amorphous silicon obtained from density functional theory calculations

Atomic structure of amorphous silicon consistent with several reported experimental measurements has been obtained from annealing simulations using electron density functional theory calculations and a systematic removal of weakly bound atoms. The excess energy and density with respect to the crystal are well reproduced in addition to radial distribution function, angular distribution functions, and vibrational density of states. No atom in the optimal configuration is locally in a crystalline environment as deduced by ring analysis and common neighbor analysis, but coordination defects are present at a level of 1%–2%. The simulated samples provide structural models of this archetypal disordered covalent material without preconceived notion of the atomic ordering or fitting to experimental data.

0.07eV atom −1 [6]. Available a-Si models often correspond to higher excess energy [7,14,[18][19][20][21], suggesting that they might not be representative of the experimental a-Si structure. Another important property is the density, ρ, which has been found experimentally to be 1%-2% lower than the density of c-Si in several measurements [12,22,23]. As shown below, the density is useful for assessing the quality of a structural model, and the convergence of annealing simulations.
The generation of a-Si structures without preconceived notion of the atomic ordering or input from experimental observations has proven to be challenging. One possibility is to simulate a rapid quench of the liquid. But, due to computational limitations, empirical potential functions are often used to describe atomic interactions. Their limited accuracy has resulted in a high density of structural defects [24,25], and an overall poor agreement with available experimental information on a-Si. In order to improve on this approach, we have carried out annealing simulations with a density functional theory (DFT) description of the electronic degrees of freedom so the final structure is dictated by atomic interactions determined from first principles calculations. In order to speed up the annealing process and optimize the density as well as the energy, a systematic procedure for removing weakly bound atoms was developed. This turns out to be an efficient procedure for generating an a-Si structure with calculated properties in excellent agreement with reported experimental measurements of ion implantation samples, including RDF, angular distributions and vibrational density of states, as well as the excess energy and relative density with respect to c-Si.
The starting configuration was generated by rapid cooling of liquid Si represented by a 216 atom system subject to periodic boundary conditions with volume chosen to correspond to that of the crystalline phase of silicon. The temperature was reduced at a rate of 2×10 12 K s −1 , typical of laser experiments [24], down to a temperature of 10K. The empirical potential function of Tersoff [26] was used in this initial phase of the simulation to reduce the computational effort. Following the cooling, the energy was minimized with respect to an isotropic scaling of the volume of the simulation box and with respect to atomic coordinates before any DFT calculations were done. The structural characteristics of the configuration obtained in this way is in some ways already in relatively good agreement with experiments although there are also clear discrepancies. The RDF is close to the experimental one, but the 1st peak is shifted 2% compared with experiments towards larger pair distances [23]. Also, while the angular distribution function has a maximum at 108.1°, in agreement with measurements giving 108.5° [2], and 107.8° [23], the dispersion is too large, 15.4°, compared with experimental values ranging from 9.0°to 11° [2,6,23].
The sample was then annealed by carrying out simulations where the energy and atomic forces were evaluated with DFT. The PBE generalized gradient approximation to the exchange and correlation functional [27], was used with valence electrons represented in a plane wave basis extended to 18Ry and only Γ-point in the k-point sampling while inner electrons were represented with projector augmented-wave [28]. This gives a lattice constant of 5.468Å for the Si crystal. The VASP software was used for the calculations [29]. A classical dynamics simulation was carried out at a temperature of 600 K followed by short, 1 ps dynamics at 300K and finally conjugate gradient energy minimization with respect to atomic coordinates at fixed volume. Figure 1 shows how the excess energy of the sample depends on the length of the time interval simulated at 600K. While the first 10 ps had a large effect, an extension to 20 ps gave only insignificant further reduction in the energy.
The density at this point is too high compared with the density of the Si crystal at the DFT/PBE level of theory since the volume had been fixed at the volume of the crystal described by the empirical potential [26]. A minimization of the DFT/PBE energy with respect to the volume box was then carried out by an isotropic scaling of the atom coordinates, followed by minimization of the energy with respect to atom coordinates. The excess energy of this configuration, labeled e in figure 1, is 0.167 eV atom −1 . This is significantly lower than previously published theoretical values [10,14,18,20,30], but higher than the experimental estimates [1,6,17]. The density, is also too high, close to the DFT/PBE crystal density, and the position of the 1st peak in the RDF is almost the same as in c-Si. The dispersion of the 1st RDF peak is 0.045Å. This is similar to experimental values [23]. The angular distribution has a peak at 108.4°with a dispersion of 10.5°. This indicates that the DFT/PBE description of the atomic forces is significantly better than the empirical potential function [26] used to generate the initial configuration. In fact, we have found that much longer annealing simulations with the empirical potential function were not able to lower the excess energy to the extent that has been done at this point in the DFT/PBE calculations. Using a cutoff bond distance of 2.8Å, configuration e is found to include 7 fivefold and 1 threefold coordinated atoms, corresponding to 4% coordination defects, and the average coordination number is 4.03.
In order to optimize the atomic structure and density further, an annealing scheme was developed that is more efficient than direct dynamical simulations. This is important since the time interval that can be simulated with reasonable computational effort is very limited. Two different calculations were carried out, one starting from configuration e which has high density, and the other from a low density configuration. In the end similar density was obtained for the optimal, low energy configurations, confirming the efficiency of the approach. The presence of coordination defects suggests that further optimization of the atomic structure and density could be achieved by removing atoms with low binding energy. Previous investigations of point defects in a-Si models have shown that in some cases the vacancy formation energy is negative, i.e. the energy of the system could be lowered by removing an atom [31]. An analysis of the structure of bond defects in our samples revealed that fivefold coordinated atoms tend to cluster together, as had been noted in some previous investigations [24]. We performed a systematic search for such 'high energy atoms' using a procedure which in the following will be referred to as point defect annihilation (PDA): (1) An atom is removed followed by a minimization of the energy with respect to the coordinates of the rest of the atoms. This is repeated for all the atoms in the system and the final configuration with the lowest energy selected for further processing. (2) An annealing of the selected configuration is then carried out for 1ps at 600K and for1 ps at 300K followed by minimization of the energy with respect to coordinates and rescaling of the volume. Steps (1) and (2) were repeated until the PDA procedure led to an increase in the final energy. We note that due to the annealing and relaxation with respect to the volume, the PDA procedure can lead to an increase in the density despite the fact that atoms are removed from the system.
The results of the PDA calculations starting with configuration e are shown in figure 1. The excess energy dropped from 0.167 to 0.156 eV atom −1 after four iterations (removing 4 atoms, configuration i in figure 1).
Since the PDA procedure described above started from a configuration with high density, close to that of the crystal, it is possible that the configuration obtained is biased towards high density. Therefore, a different initial configuration for the PDA procedure was generated by increasing the volume of the simulation box by ca 2% followed by annealing at 900, 600 and 300K for 10ps at each temperature setting, and then energy minimization with respect to the atom coordinates and an isotropic scaling of the simulation box (configuration x in figure 1). The density of this configuration is 1.8% lower than the density of the crystal and the excess energy is 0.160 eV atom −1 . When the PDA procedure is applied to this configuration, the density increases to a similar optimal value (sample y in figure 1) as the optimal configuration obtained in the previous PDA calculation, ca 99% of the crystal density. This optimal density has, therefore, been obtained with the PDA procedure starting both from values above and from below.
This second PDA annealing simulation generates, however, a configuration with even lower excess energy than obtained before, 0.149eV atom −1 after removal of just one atom. This value of the excess energy is much lower than what has been reported previously from model calculations and is in good agreement with experimental estimates for samples created by ion implantation [1,17]. We note that the total simulated annealing time used to generate this optimal configuration is quite short, 53 ps.
At the optimum for the PDA annealing, the energy of the system is not lowered further by removing an atom, i.e. no atom gives rise to a negative vacancy formation energy as shown in figure 2. As expected, the PDA procedure indeed results in a reduction of the number of coordination defects. The configuration with the , , : after simulating dynamics at 600K using DFT (for 1, 10 and 20 ps, resp.) followed by energy minimization with respect to atom coordinates at fixed volume. Two PDA annealing calculations are then carried out. I (black): energy was minimized with respect to the volume of the box (black dashed line) prior the PDA treatment, which is an iterative removal of atoms with the smallest binding energy and volume optimization. After the first iteration configuration f is obtained, g results after removing two atoms, etc, and the PDA annealing leads to the optimal, low energy configuration labeled i containing 212 atoms. II (red): the volume is increased to reduce the density to ca 2% below that of c-Si and dynamics simulated at 900, 600 and 300K for 10ps in each case, followed by minimization of the energy with respect to both atomic coordinates and volume (red dashed line), leading to a lower density sample, x. From the following PDA calculations, the optimal low energy configuration labeled y with 215 atoms results. The optimal configurations obtained from the PDA annealing I and II have nearly the same density, 99% of the c-Si density.
lowest energy per atom (sample y in figure 1) only includes one threefold and three fivefold coordinated atoms. Sample i in figure 1 which has larger excess energy actually contains fewer coordination defects, only one fivefold and one threefold coordinated atom. The excess energy is therefore not necessarily related to the number of coordination defects, in agreement with an earlier study by Bernstein et al [20]. This indicates that it can be energetically advantageous to introduce coordination defects in a-Si. The abundance of coordination defects in these samples corresponds to 1%-2% of the atoms, close to estimates for samples prepared by ion implantation [23].
The optimal density obtained in these simulations, 99% of the crystal density, is within the range of values determined experimentally for a-Si samples prepared by ion implantation [6,22,23]. Some measurements have given a density as low as 98% but, as stated by the authors, it is likely that those samples include a small concentration of vacancies even after annealing [23]. This is consistent with the low value of the coordination number determined for those samples, 3.88. In all our optimized samples, the coordination number is higher, 3.99-4.02. Another explanation could be that the experiments are carried out at non-zero temperature, and therefore contain an equilibrium concentration of point defects such as vacancies. While thermal defect concentration is usually small in c-Si due to the high formation energy of about 3.6eV, vacancies in a-Si can be created at a much lower cost. In fact, the computed distribution of the vacancy formation energy in our optimal a-Si sample, shown in figure 2, includes a value as low as 0.13eV.
Other structural characteristics of the lowest energy configuration are also in close agreement with available experimental data, as shown in figure 3. In order to compare with room temperature measurements, the RDF was calculated as an average over a large number of configurations generated by DFT/PBE classical dynamics after heating the lowest energy configuration, sample y, to 300K. The agreement is excellent, as shown in figure 3. The 1st peak position in the RDF is similar to the 1st neighbor distance in c-Si as has been found in experiments [23]. The bond angle distribution is also in close agreement with experiment as shown in figure 3 [23]. The dispersion in the bond angle distribution is also in excellent agrement with experimental estimates [2,6,23].
It has been postulated that medium-range order in a-Si is due to the presence of crystalline regions [12], and models of a-Si have been proposed that include crystalline grains [14]. These models can reproduce the variance of the measured diffraction signal in fluctuation electron microscopy experiments, but are characterized by large excess energy compared to c-Si, larger than the experimental estimates. One way to characterize medium-range order is the dihedral angle distribution. The samples generated here have distributions with well defined peaks at 60°and 180°, a clear indication of medium-range order, see figure 3.
A more detailed analysis of the medium range order can be carried out using ring analysis, i.e. count the number of atoms in shortest rings formed by Si-Si bonds [32]. The results of such analysis is shown in table 1. In Figure 2. Distribution of the vacancy formation energy obtained from the density functional theory calculations. The upper panel is for configurations labeled b and d in figure 1 and shows the effect of extending the annealing at 600 K from 1 ps (dashed line) to 20 ps (solid line and gray-shaded area). The total energy of the sample decreases during the extended annealing period but a few of the atoms become energetically less stable as indicated by larger negative values of the vacancy formation energy. The lower panel shows results for the optimal configuration after the PDA procedure (y in figure 1) where the weakly bound atoms have been removed followed by short time annealing simulations. There, all values of the vacancy formation energy are positive, the smallest one being 0.13 eV (the curve has a tail into the region of negative values because of Gaussian smearing).
c-Si, each atom forms 6 rings with 6 atoms. Such analysis has been used to identify atoms locally in crystalline environment in models of a-Si [19]. Our optimal configuration contains no atom with the ring characteristics of c-Si.
Another structural analysis method that has been used in studies of amorphous packings of atoms is based on counting the number of common neighbors and analyzing the number and arrangement of bonds between common neighbors [33][34][35]. This is referred to as common neighbor analysis (CNA). The method has recently  been extended to systems where the atoms ideally have fourfold coordination [36]. There, the analysis is carried out for second nearest neighbors rather than first nearest neighbors. An atom in c-Si forms 12 pairs of the (421) type with its second neighbors, while an atom in a hexagonal diamond lattices forms 6 pairs of (421) type and 6 pairs of (422) type. The results of CNA for our optimal configuration are shown in table 1. None of the atoms in our optimal configuration forms the set of pairs characteristic of crystalline ordering. Analogous to the results of the ring analysis, no atom is found to be in a local environment characteristic of the crystal. This shows that medium-range order could be present in well relaxed a-Si models without the presence of crystal grains. The vibrational density of states, computed within the harmonic approximation for our optimal configuration is also shown in figure 3 along with experimental data obtained by neutron scattering of a sputter deposited sample of a-Si [37]. The force constant matrix was evaluated from the DFT/PBE calculations. The agreement between the calculation and the measurement is excellent, further supporting the high quality of the atomic configuration.
Finally, we have made detailed comparison of our optimal configuration with a high-quality CRN sample [10] containing 200 atoms. The generation of such samples is not based on the energy but is subject to the constraint that all atoms are fourfold coordinated, i.e. no coordination defects. After minimization of the DFT/ PBE energy with respect to atom coordinates and volume starting with the coordinates of the CRN, the excess energy with respect to the crystal is found to be practically the same as that of our lowest energy sample, 0.149eV atom −1 . This is quite low, comparable to our optimal configuration. However, the density of the relaxed CRN configuration is significantly lower, only 96.7% of the crystal density, well below experimental estimates. The CRN configuration is also quite different from our optimal configuration in many respects. It contains no coordination defects even after relaxation with DFT/PBE while our lowest energy sample has 2% coordination defects. This again raises the question whether some small density of coordination defects is in fact an essential characteristic of a-Si. The results of ring and CNA of the CRN configuration is given in table 1. No atom is classified as being in a crystalline environment by either analysis method, analogous to our optimal configuration. But, there is a clear difference in the abundance of rings, the CRN having fewer sixfold rings and more fiveand sevenfold rings. The CNA results show a larger abundance of pairs corresponding to 3 common second neighbors while the abundance of pairs corresponding to 4 and 5 is lower. The abundance of (555) and (544) pairs, which are elements of an (distorted) icosahedron, is only ca 2/3 of that in our optimal configuration.
The electronic density of states calculated with DFT/PBE for the lowest energy configuration gives a band gap that is a factor of 1.48 larger than that calculated for the crystal at the same level of theory. While the calculated band gaps are known to be too small in DFT/PBE, this ratio is in excellent agreement with the ratio of experimentally determined band gaps, 1.55 (see [38]). Two defect peaks are found in the gap, 0.09 and 0.15eV above the valence band edge, presumably corresponding to the coordination defects. The DFT/PBE calculations of the relaxed CRN sample, however, showed no defect states within the electronic band gap consistent with the fact that no coordination defects are present. Further study of the electronic properties of the various a-Si samples will be reported elsewhere.
In summary, we have used an annealing procedure based on electron DFT to generate an atomic configuration for amorphous silicon that is in excellent agreement with the various experimental measurements, including excess energy with respect to the crystal, density, coordination defect density, RDF, angular distribution functions and vibrational density of states. Our approach is in part based on a systematic procedure for removing weakly bound atoms. The optimal configuration obtained has no atoms in a local environment characteristic of crystalline order as deduced from ring analysis and CNA.