LDA Code: Background and Technical Details
We investigate the microscopic properties of semiconductors in both bulk and surface conditions through the use of a sophisticated LDA computer code. This code, originated by Dr. GP Srivastava, has undergone over fifteen years of development and is continually being refined to extend its accuracy and range of applicability. Essentially the code implements the hugely successful density functional theory (DFT) of Hohenberg and Kohn (1964) and Kohn and Sham (1965). In this approach the ground state electron density of a system is considered to be its fundamental variable (rather than the much more complicated many-body electron wavefunction). All other ground state properties of the system are then expressible (in principle) as functionals of the ground state electron density.
The Kohn-Sham equation (derived from the variational principle with respect to fluctuations in the electron density distribution) provides a means of determining the electron density distribution in practice. The equation is Schrodinger-like in form, with a Hamiltonian composed of an electron kinetic energy term, an electron-ion interaction potential, the Hartree potential (this is the potential due to the electron density distribution neglecting exchange and correlation effects) and an unknown exchange-correlation potential. This last potential accounts for the exchange interaction between electrons (due to the Pauli exclusion principle) and the electrostatic screening of each electron by the correlated motion of all the other electrons. In truth the evaluation of this potential is a subtle many-body problem of many years standing and is still a matter for current debate in the literature, however it has long been known that the so-called local density approximation (LDA) performs remarkably well in this role. The LDA is based upon the simple idea of using the known exchange-correlation potential of the homogeneous electron gas with density equal to the local density of the true inhomogeneous electron gas in question.
Most quantities appearing in the Kohn-Sham Hamiltonian are dependant upon the electron density distribution, which, of course, cannot be known until the Kohn-Sham equation has been solved. Thus the Kohn-Sham equation is commonly solved by making an initial guess for the electron density distribution, allowing the Hamiltonian to be constructed and a better approximation to the true electron density distribution to be calculated. The process can then be repeated until an electron density distribution of the desired accuracy has been achieved.
A popular alternative method is to solve the Car-Parrinello (1985) equation instead of the Kohn-Sham equation. This is simply a Lagrangian reformulation of the Kohn-Sham equation in which the electronic wavefunctions are treated as co-ordinates for fictitious particles with fictitious mass. Of course an iterative method of solution is still required but it the advantage is that the ionic degrees of freedom and the electronic degrees of freedom are now treated on the same footing and the ground state with respect to both can be found simultaneously (i.e. the atoms relax into their minimum energy configuration while the electron density is being refined). In the Kohn-Sham scheme it is necessary to solve the electronic degrees of freedom iteratively for a given atomic configuration and then move the atoms towards their minimum energy configuration guided by the calculated forces acting upon them.
The foregoing discussion is entirely general and common to many DFT/LDA computer codes. In the following sections we shall discuss some of the more technical aspects of our code (which is of the Kohn-Sham type), some of which it shares with others and some of which are unique....
There are advantages and disadvantages to both methods. In the real space approach it is natural to describe isolated clusters of atoms and so the method is suited to investigating properties of large molecules, or of isolated defects in bulk materials (i.e. if the cluster is large enough and surface dangling bonds are passivated with Hydrogen atoms then an extended solid can be modelled). On the other hand the momentum space approach copes with extended systems (e.g. surfaces) in a much more natural way because of its use of a unit cell (or large supercell) and periodic boundary conditions.
There are two LDA codes in use at Exeter: our own plane wave code described here, which we use primarily for surface calculations, and the cluster code used by Dr. Jones' group to study defects.
In either case the electron density which is calculated at each iterative step may simply be found by integrating the squares of the magnitudes of the Kohn-Sham eigenfunctions over the Brillouin zone. The integration is done in practice by summing over a relatively small number of k-points in the irreducible segment of the Brillouin zone and then using the lattice symmetry to extend this result to the rest of the zone. By choosing the k-points used with some care and weighting them appropriately it is possible to reduce the number of points needed to an absolute minimum. Typically we use less than ten special k-points depending on circumstances (see Evarestov and Smirnov 1983).
Once again there are two common approaches used to overcome this problem. In the so-called all-electron methods the region immediately surrounding each ion is dealt with analytically using some simplifying approximation (e.g. assuming spherical symmetry). A basis set is then chosen composed of functions which match the core region solutions within the relevant regions but which reduce to simple forms in interstitial space.
The other popular approach is the pseudopotential approach in which a specially constructed non-singular electron-ion potential is used in place of the true potential. This pseudopotential is non-singular because it includes the screening of the nucleus by the tightly bound core electrons in a special way: it is constructed to project out the rapidly oscillatory part of the Kohn-Sham solutions. It can be formally shown that a properly constructed pseudopotential yields the same results as the true electron-ion potential, but with much better convergence properties. Our code employs ab initio, non-local, norm conserving pseudopotentials (that is they are constructed from first principles, they include the effect of a variety of different orbital angular momentum states of the core electrons, and they give the correct wavefunctions in the interstital regions). We usually use the Bachelet-Hamann-Schluter pseudopotentials (1982), but when necessary we can use more modern smoother pseudopotentials instead.