The determination of molecular and crystal structure:
The Fundamentals of AIMPRO

To determine the structure of an assembly of atoms, we need to solve the Schrödinger equation for both the ions and electrons. It is impossible to give an exact solution: progress can only be made if a set of approximations are made.

The first is due to Born and Oppenheimer who argue that, since the nuclear masses are so much larger than those of the electrons, then we can treat the nuclear motion classically and reduce the Schrödinger equation to one only involving the electrons moving in a potential of fixed nuclear sites denoted by R. The solution of this equation is the structural potential energy E(R). To determine dynamical quantities like vibrational frequencies and diffusion barriers, we need then to solve the classical problem of the nuclei moving in the force field tex2html_wrap_inline99 .

Let us now consider the Schrödinger equation of the electrons in a field of frozen nuclei. This equation still cannot be solved exactly because of the inter-electronic Coulomb terms, tex2html_wrap_inline101 . If these terms were neglected, then the Schrödinger equation is separable and the exact wavefunction can be written down as a sum of product orbitals tex2html_wrap_inline103 where tex2html_wrap_inline105 are the eigenfunctions of the one-electron Hamiltonian. Of course, we must require that the total wavefunction is antisymmetric and this can be achieved by a suitable sum over the product orbitals. For, example for the H tex2html_wrap_inline107 molecule, the antisymmetric wavefunction is

displaymath109

This type of wavefunction is called a Slater determinant and is specified by the electronic configuration, ie. a listing of the orbitals tex2html_wrap_inline111 occupied by the N electrons.

The neglect of the electron-electron interaction term is too extreme to make this theory useful and various attempts have been made to take this interaction into account. There are two principal schemes: Hartree-Fock (HF) and density functional theory (DFT). Both of these replace the electrostatic potential acting on electron i by the other electrons, ie.  tex2html_wrap_inline117 by an average over the positions of the other electrons. Thus this potential becomes tex2html_wrap_inline119 . Clearly the Schrödinger equation is then separable and we can write the wavefunction as a Slater determinant. In these theories, the effective potential is composed of two terms: the Hartree potential tex2html_wrap_inline121 and the exchange-correlation potential tex2html_wrap_inline123 . The former is the electrostatic potential due to the electron density of charge n(r) at point r. The density is determined by the probability of finding a electron at the point, or

displaymath129

Here, the sum is over the occupied orbitals. The exchange-correlation potential differs in the two theories. In HF, it is a complicated function determined by all the orbitals tex2html_wrap_inline105 . In DFT, it is rigorously determined only by the electron density. This is an important result but the precise dependence on density is not known except for a particular problem: the homogeneous electron gas. In this case, the equations have been solved numerically and the dependence of tex2html_wrap_inline133 on n is known exactly. For any other problem where the electron density varies throughout space, we assume that tex2html_wrap_inline123 at point r is given by the homogeneous electron gas value involving the density n(r) at the same point r. This is called local density functional (LDF) theory. For high electron densities, it can be shown that tex2html_wrap_inline133 is proportional to tex2html_wrap_inline147 and thus in LDF theory tex2html_wrap_inline123 is also proportional to tex2html_wrap_inline151 . For lower densities, tex2html_wrap_inline133 can still be written down but is given by a more complicated expression

AIMPRO uses LDF theory and thus the electron density is a quantity of fundamental significance.

The problem has now been reduced to a one-electron one involving the interaction of the electron with all the ions together with tex2html_wrap_inline155 . Not all the electrons need be considered. Electrons are divided into two groups: core and valence electrons. The former are bound close to the ions and do not play a part in bonding. It is highly advantageous to remove these from the theory. This can be done through the use of the pseudopotential. For example, in Si the 1s, 2s, and 2p electrons are regarded as part of the core while the 3s and 3p electrons are part of the valence shell. The four bonds associated with each Si atom in the bulk arise from these four valence electrons. The pseudopotential is constructed by insisting that it has exactly the same valence energy levels as the true atomic potential, ie. the 3s and 3p and even the 3d levels are the same (the 3d levels are required because these orbitals are occupied to some extent in the solid). Moreover, the corresponding pseudo-wavefunctions are exactly equal to the true wavefunctions outside a core whose size depends on the identity of the atom. Inside the core, the pseudo-wavefunctions are not equal to the true valence ones but are very smooth nodeless functions which are easy to represent mathematically. The true 3s and 3p wavefunctions of Si oscillate inside the core and are rather difficult to represent mathematically. Fig. 1 shows these functions for Ni.

  figure56

Figure: The 4s (full) and pseudo- (dashed) radial wavefunction (a.u.) for the Ni atom.

The crucial point is that the pseudopotentials have no core states at all! The lowest bound state solutions are the valence energy eigenvalues and, as stated above, these are the same as the true valence energy levels. The pseudopotential has the great virtue of making the treatment of Ge say no more difficult than carbon. AIMPRO uses the pseudopotentials due to Bachelet, Hamann and Schlüter (BHS) who have given tables for all the elements from H to Pu.

Despite these great simplifications, we are still confronted with a formidable problem. We need to solve the eigenvalue problem:

displaymath179

Here H includes the kinetic energy, the sum of the pseudopotentials for all the ions present, and the effective potential tex2html_wrap_inline183 which depends on the density of electrons n(r). To proceed we expand the wavefunction in terms of a basis tex2html_wrap_inline187 . Thus

displaymath189

Two choices of tex2html_wrap_inline191 are often made: plane waves and Cartesian Gaussian orbitals. Using the first is equivalent to making a Fourier transform of the wavefunction. This has problems with certain elements, eg. first row, transition and rare-earth ones where the wavefunction varies rapidly and a great number of plane waves are required. AIMPRO uses Cartesian Gaussian orbitals of the form

displaymath193

when centred at the origin. In practice, they are centred at the nuclei and sometimes at bond centres sites between the nuclei. If tex2html_wrap_inline195 then the orbitals are simply spherically symmetric Gaussian functions sometimes called s-Gaussian orbitals. If one of the tex2html_wrap_inline199 is unity, and the others zero, then the functions are p-Gaussian orbitals while if the sum of the tex2html_wrap_inline199 is 2, the set of six orbitals generate five d-orbitals and one s orbital. The great advantage of Gaussian orbitals is that analytical expressions can be given for all the integrals.

We now have to solve the equations

displaymath209

The next step is to multiple this equation through by tex2html_wrap_inline211 and integrate over r. This reduces the differential equation to a set of matrix equations:

displaymath215

which can be solved by matrix methods. tex2html_wrap_inline217 and tex2html_wrap_inline219 are matrix elements of the Hamiltonian and the overlap respectively.

Once this equation has been solved we can generate the output charge density

displaymath221

In general this is not equal to the charge density used to generate the effective potential tex2html_wrap_inline183 . Of course, this is wrong and we need to devise a scheme whereby the input charge density, which is used to construct tex2html_wrap_inline183 , is equal to the output density. This is usually carried out by an iterative procedure and the process of achieving equality in the densities is called the self-consistent cycle.

AIMPRO uses an internal representation of the charge density by expanding it in a basis of s-Gaussian functions tex2html_wrap_inline229 . Thus

displaymath231

The particular exponents of the Gaussian tex2html_wrap_inline233 together with the exponents on the basis functions tex2html_wrap_inline187 have been calculated for atoms and given in a single data file called bhs.lib.

Once the self-consistent density has been found, the structural energy E(R) can be evaluated. We then compute he forces on each atom from tex2html_wrap_inline239 and move the atoms until this energy is a minimum and these forces vanish. This is called relaxing the assembly of atoms.

The user of AIMPRO needs to provide 1) the positions of the atoms in the cluster (in atomic units); 2) the basis to be used for the wavefunction and charge density; 3) the electronic configuration; 4) the identities of the atoms which are to be relaxed. The output gives the energy, electronic levels tex2html_wrap_inline241 , and positions of the atoms. More specialised output is available.