University of Exeter Shield

Plane Wave LDA Code


We investigate the microscopic properties of semiconductors in both bulk and surface conditions through the use of a sophisticated LDA computer code, known as EKSETER (standing for Evaluation of the Kohn-Sham Equation for Total Energy and Relaxation). 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 the EKSETER code (which is of the Kohn-Sham type), some of which it shares with others and some of which are unique....

Basis Set Expansion in Plane Waves and Brillouin Zone Integration

In order accurately to describe the electron wavefunctions of the system being studied it is necessary to choose a suitable set of simple basis functions in terms of which one can write a series expansion of the (much more complicated) electronic wavefunctions. There have been several approaches to this matter over the years. One is to consider the most natural basis functions from a real space viewpoint, namely atomic-like basis functions (e.g. Gaussians). Alternatively one could employ a basis set more suitable for a momentum space description of the material, namely plane waves. Our code follows the latter approach.

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 EKSETER plane wave code described here, which we use primarily for surface calculations, and the cluster code (AIMPRO) 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).

Pseudopotential Description of the Electron-Ion Interaction

Whether one chooses a real space or momentum space representation for the wavefunctions there remains a problem associated with the electron-ion interaction. The inverse square law of electrostatic interactions dictates that there is a singularity at each ionic site (where the electron-ion potential tends to minus infinity). This singularity means that simply increasing the size of one's basis set will not necessarily lead to satisfactory numerical convergence in the electron density and total energy derived from the Kohn-Sham equation.

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.

Diagonalising the Hamiltonian and Updating the Potential

The code supports three alternative methods for diagonalisation of the Hamiltonian matrix. The first is the standard Householder diagonalisation scheme, but for large Hamiltonians this becomes slow and memory intensive. A better scheme supported by the code is to iteratively diagonalise the Hamiltonian by the conjugate gradient method. Finally there is the RM-DIIS (residual minimisation by direct inversion in iterative subspace) which is applied in a semiblock form so that only the lowest eigenvectors (i.e. those which are occupied and contribute to the electron density) are calculated. After the Hamiltonian is diagonalised to find the eigenvectors the electron density distribution may be calculated, but this is not used directly as the input for the next iterative solution of the Kohn-Sham equation since this would lead to unacceptably slow convergence. Instead the next input potential is found by updating the previous potential with the new potential corresponding to the newly calculated density according to a quasi-Newton-Raphson (Broyden's Jacobian update) scheme.

Relaxation of Atoms

Having solved the Kohn-Sham equation for a given ionic configuration it is possible to calculate the forces on the ions using the Hellmann-Feynman theorem. Using these forces as a guide the code then moves the atoms to a second "test" configuration and solves the Kohn-Sham equation once again. At this point we know the total energy of both configurations and all the forces in both configurations. This information allows the code to minimise the total energy by finding some third configuration somewhere between the first two. Subsequently the code refines the ionic configuration by making successive adjustments in a simliar manner. At each step the "test" configuration is found using both the Hellmann-Feynman forces and the constraint that each relaxation must be conjugate to all previous relaxations in the configuration space. Thus the iterative relaxation scheme is an example of a conjugate gradient minimisation scheme. In this way we are able to start with a plausible guess at the microscopic atomic structure of a system and confidently find at least a local energy minimum. Of course, we cannot guarantee that this is a global energy minimum (i.e. the ground state), but if we start with a variety of plausible configurations and show that one specific end-configuration is lower than all others that we find then we can be quite confident in assigning this to be the ground state. The emphasis must thus be on a good understanding of the physical processes which dictate the structure and the discipline to carefully narrow down a range of possibilities to find the correct one.

Limitations of the Method

The DFT/LDA approach has an excellent track record in calculations of the ground states for a wide variety of materials. Atomic structures can be calculated with great accuracy, especially if one employs one of the improved LDA exchange-correlation potentials such as that of Ceperley and Alder (1980). However, there remains one intrinsic limitation of the DFT method per se: its validity derives entirely from the variational principle and as such it is unsuitable for describing excited states. The band structure of a solid is thus beyond the range of DFT to describe, since the valence and conduction bands must be considered to be excited states either of holes (valence bands) or electrons (conduction bands). In practice the problem only manifests itself in materials possessing a band gap, so metal band structures are generally reasonable within the LDA while semiconductor and insulator band gaps are underestimated by 50-100%. To overcome this problem it is necessary to calculate corrections to DFT band stuctures by a so-called quasiparticle method.

References

 B

A
Up. Directory Tree.
Stephen Jenkins
Last modified: Fri May 5 11:58:14 1995