Accurate structure factors from pseudopotential methods
Highly accurate experimental structure factors of silicon are available in the literature, and these provide the ideal test for any ab initio method for the construction of the all-electron charge density. In a recent paper [J. R. Trail and D. M. Bird, Phys. Rev. B 60, 7863 (1999)] a method has been developed for obtaining an accurate all-electron charge density from a first principles pseudopotential calculation by reconstructing the core region of an atom of choice. Here this method is applied to bulk silicon, and structure factors are derived and compared with experimental and Full-potential Linear Augmented Plane Wave results (FLAPW). We also compare with the result of assuming the core region is spherically symmetric, and with the result of constructing a charge density from the pseudo-valence density + frozen core electrons. Neither of these approximations provide accurate charge densities. The aspherical reconstruction is found to be as accurate as FLAPW results, and reproduces the residual error between the FLAPW and experimental results.
pacs:78.70.Ck, 71.15.Hx, 71.15.Ap
Pseudopotential methods, particularly within the framework of total energy plane-wave calculations, are extremely powerful for the ab initio description of large system of atoms due to their computational efficiency and suitability for structural optimisationpayne92 (). However, they do not yield the correct charge density of the system studied, but a ‘pseudo’ charge density that does not include core electrons and is incorrect close to atomic nuclei. This precludes the direct application of these methods to the prediction of any properties of the material that depend directly on the charge density, such as hyperfine couplings or X-ray structure factors. In this paper we reconstruct the all-electron charge density for bulk silicon from a pseudopotential calculation using the method described by Trail and Birdtrail99 () (hereafter referred to as I), and from this we derive the X-ray structure factors. These are then compared with experimental results and the results of other theoretical approximations and methods to assess the importance of various assumptions often made in the calculation of structure factors, and to evaluate the success of the reconstruction method.
Previous methods gardner86 (); vackar94 (); kuzmiak91 (); meyer95 () for solving this reconstruction problem have relied on the assumption that the potential in the core region is spherically symmetric in order to decouple the differential equations that must be solved, and in many cases the charge density itself is assumed to be spherical. The method used here does not require this to be the case, and we compare the structure factors resulting from assuming spherical symmetry to justify the extra effort necessary to develop an aspherical reconstruction procedure. The reconstruction method itself is based around the embedding approach of Inglesfield inglesfield81 (). Results from this localised calculation are used to replace the pseudo-charge density where this is incorrect, leading to the required structure factors.
In section II a brief summary of the reconstruction approach is given - a full description can be found in I. We describe how we obtain the structure factors from the reconstruction in section III, and compare the reconstruction results with accurate experimental and theoretical structure factors. Rydberg atomic units are used throughout the paper.
Ii Reconstruction Method
The first step in the reconstruction procedure is to obtain an accurate approximation for the real space single particle Green function of the substrate system, in this case bulk silicon. We begin with a total energy pseudopotential calculation performed with a plane-wave basis and using the Local Density Approximation (LDA) for exchange and correlation. A plane-wave energy cut-off of eV is used and Monkhorst-Pack -pointsmonkhorst76 () are included in the irreducible wedge of the FCC Brillouin zone. These values are more than sufficient to obtain essentially perfect convergence of the self-consistent density and potential, which allows us to attribute any errors in our results to the reconstruction procedure. A norm-conserving Kerkerkerker80 () pseudopotential is used, with a maximum core radius of au. As explained in I the method requires to be less than half the nearest-neighbour atomic separation in the crystal. The resulting self-consistent potential is used to obtain a set of eigenstates by direct matrix diagonalisation, at -points in the irreducible wedge of the Brillouin zone and with a eV plane-wave energy cut-off. Careful tests have been carried out to confirm that these values are sufficient to provide structure factors with a precision of order milli-electrons/atom. The spectral representation morse53 () is used to construct a Green function from the set of plane-wave states. This Green function is then used to obtain an embedding potential, which is a term that is added to the Kohn-Sham Hamiltonian for the localised core region of an atom of interest. The effect of the embedding potential is to take into account the lattice of pseudo-atoms surrounding the chosen atom. The localised embedded Hamiltonian is then solved self-consistently (again using the LDA) to obtain the Green function of the embedded system, from which the charge density in the core region can be obtained (see paper I).
The reconstruction is performed using the same parameters as in I, with the embedding radius chosen as the ‘touching spheres’ radius (in this case au). Again, convergence with respect to all parameters has been thoroughly checked. The final result we arrive at is for the charge density of a single all-electron atom embedded in a lattice of pseudo-atoms. This information can be used together with the original pseudo-charge density to construct an accurate all-electron charge density and hence the structure factors for the crystal.
Iii Structure Factors
Extremely accurate experimental structure factors for silicon have been available in the literature for some time cumming88 (); aldred73 (); teworte84 (); saka86 (). These results have been used by a number of workers to assess the accuracy of parameterised models deutsch92 (), FLAPW and other ab initio methods lu93 () and Generalised Gradient Approximations to the exchange-correlation potential zuo97 (). In view of the accuracy and range of data available, both experimental and theoretical, the reconstructed silicon charge densities are here used to construct structure factors for comparison with experimental data and the results of FLAPW calculations.
iii.1 Structure Factors from Reconstructed Charge Densities
To obtain the structure factors we require the Fourier coefficients of the charge density,
where denotes the unit cell volume and is the total real space charge density. consists of the original pseudo-charge density between atoms, and the reconstructed total charge density within the embedding sphere surrounding each atom. Since this integral is a linear operation on the charge density it is possible to subtract the contribution to the pseudo-density from the embedding regions around each atom and add on the contributions from a reconstruction calculation. This gives the expression
where are the position vectors of the atoms in the unit cell. The quantities and are the Fourier integrals of the reconstructed and pseudo-densities respectively, carried out over the reconstruction sphere surrounding each atom, and are given by
where is the radius of the reconstruction sphere and is the appropriate charge density, reconstructed or pseudo. The original pseudo-charge densities are available in reciprocal space directly from the plane-wave calculation, and the reconstructed charge densities are given as an expansion in spherical harmonics trail99 (), which allows Eq. (2) and (3) to be evaluated.
For an atom situated at the origin Eq. (3) takes the form
where the charge density has been explicitly written as an expansion in spherical harmonics, and the identity
has been used. The radial integral in Eq. (4) is carried out numerically. In our calculations for silicon we choose the primitive unit cell, and reconstruct the core region of one atom chosen to be at the origin, hence the integral in Eq. (4) is carried out over a sphere centred on this atom. Other atoms within the unit cell must also be taken into account, and in the case of silicon there is another atom at related to the origin by an inversion symmetry at . The contribution to Eq. (2) from this atom can be derived from the symmetry of the unit cell. If the atom at the origin is related to an atom at site by the space group operator (defined by ) altmann91 (), where is a unitary transformation, then the integral, is
By transforming coordinates this reduces to
For silicon the atom at is related to the atom at the origin by the operator , an inversion followed by a translation. In this case the above expression, together with the expansion around the origin in spherical harmonics, yields
where . The transformation results in the phase factor, and the inversion results in the power of present in the sum.
Eq. (4) and (8) are applied to both the reconstructed charge density and the pseudo-density (expanded in spherical harmonics) and are then substituted into Eq. (2) to yield the structure factor as a function of the reciprocal lattice vector, . At first it seems a roundabout route to calculate the radial expansion of the pseudo-density only to convert this back to a reciprocal space representation, but this is the most straightforward way of replacing the pseudo-charge density with the reconstructed charge density in the sphere around each atom. One final point is the position of the origin. The coordinate system used for the reconstruction has the origin on one of the silicon atoms in the unit cell (at ), whereas the system normally chosen for crystallographic studies has the origin at the inversion centre, () hahn95 (). Placing the origin at the inversion centre gives real structure factors, and the origin can easily be shifted to this point by introducing an appropriate phase factor into Eq. (1), or simply by taking the magnitude of the complex structure factors.
iii.2 Comparison with Experimental Results and FLAPW Calculations
Before comparison can be made between the theoretical and experimental results two further factors must be considered. First, the experimentally measured quantity (normally given in the literature) is not the Fourier coefficient of the charge density, , but the form factor , which takes into account the lattice structure. This is defined as zuo97 ()
where are the indices of the reciprocal lattice vector. For values that satisfy the criteria for integer, the denominator on the right hand side is zero, and the structure factor is given.
The second effect that must be taken into account when correlating the theoretical and experimental results is the thermal motion of the lattice. The majority of experimental data for structure factors are taken at room temperature, and the thermal energy ‘smears out’ the charge density, reducing the amplitude of the higher order structure factors. This can be described by a convolution integral in real space, which corresponds to a further correction factor in reciprocal space to give the dynamic structure factor
Structure factors obtained from the core reconstruction are here compared with those obtained from three sources: from the simple addition of free atom core states to the original pseudo-charge density lu93 (), structure factors obtained using the FLAPW method by Zuo et al zuo97 () and structure factors determined experimentally by Cumming and Hart cumming88 () and Saka and Kato saka86 (), as given by Zuo et al zuo97 (). The pseudo+core structure factors are obtained from the charge density of the original pseudopotential calculation together with the core charge densities of the original atomic calculations used to create the pseudopotential. The contribution from the atomic core charge density is included at the atomic sites in the same manner as described above, ie
where is the contribution from the core states at site . This structure factor is expected to show significant error, since the valence charge density will be incorrect close to the atomic sites.
Zuo et al zuo97 () calculated Si structure factors using the FLAPW method, and they give results using the LDA, and two different GGAs. Since the core reconstruction calculation carried out here employs the LDA, the reconstruction results are only compared with the LDA FLAPW results. For a successful reconstruction scheme we would expect to accurately reproduce these results, since the same physical approximations have been made even though the algorithmic implementation of the two methods is entirely different.
We begin by comparing the theoretical form factors, uncorrected for temperature, and only for the values for which experimental data are available (experimental data are given in Table 1). Fig. 1 shows the difference between the reconstructed and FLAPW form factors, and the difference between the pseudo+core and FLAPW form factors. It is apparent that the reconstruction agrees very well with the FLAPW results - the average absolute difference for the reconstructed results is only 3 milli-electrons/atom whereas for the pseudo+core result the average absolute difference is over 25 times greater at 76 milli-electrons/atom.
iii.3 Experimental, FLAPW and Reconstructed Structure Factors
In order to compare the static structure factors given above with experimental data a value for the Debye-Waller parameter in Eq. (10) is required. This is commonly taken to be a free parameter and varied to minimise the error between the experimental and theoretical results. The value of used here is that employed by Zuo et alzuo97 (), Å . In their paper values of are obtained by minimising the error of high values only, for a number of different ab initio methods. These high structure factors depend largely on the core states of the atoms that make up the lattice, so the best values should result from methods that accurately describe the core states. Zuo et al found that a calculation of these high order structure factors using the Multi Configuration Dirac-Fock (MCDF) method grant80 () gives the best fit at high , and therefore took the associated parameter to be the best estimate. It should be noted that a better fit between experiment and theory can be obtained for a different value, but this would effectively use the description of a physical effect, the thermal smearing, to adjust for deficiencies in the theory, such as the LDA.
Table 1 gives the experimental data, with reconstructed, FLAPW and pseudo+core dynamic form factors. The quality of the theoretical data is assessed by two statistics - the -factor and GOF parameter. The -factor is given by
and the goodness of fit parameter by
where is the sample variance of the form factor. The variance is taken to be the average of the estimated error for all data points in line with the approach of Zuo et al, and takes the value .
|1 1 1||10.6025(29)||10.6020||10.5995||10.5824|
|2 2 0||8.3881(22)||8.3955||8.3952||8.3531|
|3 1 1||7.6814(19)||7.6879||7.6909||7.6373|
|2 2 2||0.1820(10)||0.1695||0.1615||0.1650|
|4 0 0||6.9958(12)||6.9924||6.9933||6.9287|
|3 3 1||6.7264(20)||6.7081||6.7031||6.6365|
|4 2 2||6.1123(22)||6.0890||6.0897||6.0145|
|3 3 3||5.7806(21)||5.7456||5.7552||5.6732|
|5 1 1||5.7906(27)||5.7754||5.7761||5.6984|
|4 4 0||5.3324(20)||5.3119||5.3136||5.2339|
|5 3 1||5.0655(17)||5.0447||5.0490||4.9670|
|6 2 0||4.6707(9)||4.6542||4.6561||4.5748|
|5 3 3||4.4552(11)||4.4485||4.4444||4.3661|
|4 4 4||4.1239(18)||4.1069||4.1085||4.0285|
|7 1 1||3.9282(22)||3.9213||3.9229||3.8449|
|5 5 1||3.9349(34)||3.9255||3.9248||3.8482|
|6 4 2||3.6558(54)||3.6413||3.6427||3.5671|
|7 3 1||3.4919(11)||3.4868||3.4869||3.4135|
|5 5 3||3.5055(14)||3.4805||3.4883||3.4108|
|8 0 0||3.2485(34)||3.2458||3.2470||3.1766|
|7 3 3||3.1270(14)||3.1112||3.1154||3.0453|
|8 2 2||2.9111(15)||2.9096||2.9105||2.8456|
|6 6 0||2.9143(16)||2.9095||2.9105||2.8458|
|5 5 5||2.8009(21)||2.8008||2.7947||2.7361|
|7 5 1||2.8006(25)||2.7951||2.7976||2.7341|
|8 4 0||2.6200(7)||2.6216||2.6219||2.5631|
|9 1 1||2.5325(8)||2.5232||2.5242||2.4678|
|7 5 3||2.5274(29)||2.5264||2.5229||2.4688|
|6 6 4||2.3677(9)||2.3724||2.3733||2.3208|
|8 4 4||2.1506(24)||2.1572||2.1581||2.1115|
|8 8 0||1.5325(26)||1.5365||1.5370||1.5095|
From the data in Table 1 it can be seen that the reconstruction calculation describes the experimental data as well as the FLAPW results. For both sets of data the -factor is %, and the GOF is , with the GOF for the reconstruction slightly greater than that for FLAPW. The average absolute error is 10 milli-electrons/atom for both the FLAPW and reconstruction calculations and 70 milli-electrons/atom for the pseudo+core results. The maximum error is roughly milli-electrons/atom for the FLAPW and reconstruction results, and milli-electrons/atom for the pseudo+core results. Fig. 2a shows the residual error () of the FLAPW results together with the error bars of the experimental data, and Fig. 2b the residual error for the reconstruction calculation. The errors are very similar, even to the point of a significant correlation existing between the two. This suggests that the errors present are largely due to the theory shared by the calculations, specifically the LDA. It should also be noted that the data presented by Zuo et al is calculated for a lattice constant of Å, whereas the reconstruction calculations are carried out for Å.
Finally we give the -factor and GOF parameter comparing the pseudo+core and reconstructed results with the FLAPW results. The pseudo+core form factors give a -factor and GOF of % and respectively, while the reconstruction gives % and .
iii.4 Spherical Symmetry
One of the strengths of our reconstruction method is that it does not require spherical symmetry of the charge density in the reconstruction region near the cores of the atom. To assess the importance of the aspherical components of the charge density we replace the original pseudo-density in the reconstruction sphere with the spherical part only of the reconstructed density. Fig. 3 gives the residual error of the reconstructed form factors from the experimental data for this case. The -factor is % with a GOF of - considerably worse than for the full aspherical reconstruction. From this data it is apparent that the aspherical components of the charge density are essential for the calculation of accurate form factors. However, it is interesting to note that if we replace the spherical part of the pseudo-density with the spherical part of the reconstructed density (that is the charge density components for within the embedding sphere are given by the pseudo-charge density) we obtain form factors that are almost as accurate as the FLAPW and fully aspherical results. In this case the -factor is %, the GOF is and the mean absolute error is milli-electrons/atom).
In this paper all-electron states have been reconstructed successfully from a total energy pseudopotential calculation, giving an accurate charge density in the region near atomic sites. This reconstruction is carried out using the embedding method described in a recent paper trail99 (). The reconstruction calculation itself uses a scalar relativistic description for the valence states, in a fully aspherical potential and using LAPW basis functions. The core states are calculated fully relativistically by direct solution of the Dirac equation in a spherical average of the self consistent potential. It is apparent that the reconstruction method itself has a lot in common with FLAPW methods.
Structure factors have been derived from the reconstructed silicon charge density for comparison with accurate experimental data and FLAPW calculations. Agreement is excellent with both the FLAPW and reconstructed form factors agreeing with experimental results with an average absolute error of milli-electrons/atom while the experimental data itself is accurate to milli-electrons/atom. The FLAPW and reconstructed form factors agree extremely well with each other, with an average absolute difference of milli-electrons/atom. In addition to this the residual errors for both methods of calculation show significant correlation, indicating that they arise from the physical approximations common to both methods.
Acknowledgements.This work has been supported by the United Kingdom Engineering and Physical Sciences Research Council. We thank S. Crampin and J. E. Inglesfield for helpful discussions.
- (1) M. C. Payne, M. P. Teter, D. C. Allan, T. A. Arias and J. D. Joannopoulos, Rev. Mod. Phys. 64, 4 1045 (1992).
- (2) J. R. Trail and D. M. Bird, Phys. Rev. B 60, 7863 (1999).
- (3) J. R. Gardner and N. A. W. Holzwarth, Phys. Rev. B 33, 10 7139 (1986).
- (4) J. Vackář and A. Šimu̇nek, J.Phys.: Condens. Matter 6, 3025 (1994).
- (5) V. Kuzmiak, J.Zavadil and K. Žďánský, Phys. Stat. Sol.(b) 168, 547 (1991).
- (6) B. Meyer, K. Hummler, C. Elsässer and M. Fähnle, J.Phys.: Condens. Matter 7, 9201 (1995).
- (7) J. E. Inglesfield, J.Phys. C 14, 3795 (1981).
- (8) H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 12 5188 (1976).
- (9) G. P. Kerker, J.Phys. C 13, L189 (1980).
- (10) P. M. Morse and H. Feshbach, Methods of Theoretical Physics (Plenum Press, New York, 1953).
- (11) S. Cumming and M. Hart, Aust. J.Phys. 41, 423 (1988).
- (12) P. J. E. Aldred and M. Hart, Proc. R. Soc. A 332, 223 (1973).
- (13) R. Teworte and U. Bonse, Phys. Rev. B 29, 2102 (1984).
- (14) T. Saka and N. Kato, Acta Crystallogr. A 42, 469 (1986).
- (15) M. Deutsch, Phys. Rev. B 45, 2 646 (1992).
- (16) Z. W. Lu. A. Zunger and M. Deutsch, Phys. Rev. B 47, 15 9385 (1993).
- (17) J. M. Zuo, P. Blaha and K. Schwarz, J.Phys.: Condes. Matter 9, 7541 (1997).
- (18) S. L. Altmann, Band Theory of Solids: An Introduction from the Point of View of Symmetry (Oxford University Press, Oxford, 1991).
- (19) T. Hahn (Ed.), International Tables for Crystallography: Vol. A Space-group symmetry (Kluwer Academic Publishers, Dordrecht, 1995)
- (20) I. P. Grant, B. J. McKenzie, P. H. Norrington, D. F. Mayers and N. C. Pyper, Comput. Phys. Commun. 21, 207 (1980).