Mulliken charges and bond populations are calculated according to the formalism described by Segall et al. (1996a and 1996b).
A disadvantage of the use of a plane wave (PW) basis set is that, due to the delocalized nature of the basis states, it provides no information regarding the localization of the electrons in the system. In contrast, a Linear Combination of Atomic Orbitals (LCAO) basis set provides a natural way of specifying quantities such as atomic charge, bond population, charge transfer etc.
Population analysis in CASTEP is performed using a projection of the PW states onto a localized basis using a technique described by Sanchez-Portal et al. (1995). Population analysis of the resulting projected states is then performed using the Mulliken formalism (Mulliken, 1955). This technique is widely used in the analysis of electronic structure calculations performed with LCAO basis sets.
The eigenstates |ψα(k)>, obtained from the PW calculation when sampling at a given wavevector k, are projected onto Bloch functions formed from a LCAO basis set |ϕμ(k)>. In general, such a localized basis set will be neither orthonormal nor complete. Therefore, care must be taken when performing the projection.
The overlap matrix of the localized basis set, S(k), is defined as:
Eq. CASTEP 63
The quality of the projection may be assessed by calculating a spilling parameter:
Eq. CASTEP 64
where Nα is the number of PW states, wk are the weights associated with the calculated k-points in the Brillouin zone and p(k) is the projection operator of Bloch functions with wavevector k generated by the atomic basis set:
Eq. CASTEP 65
where |ϕμ(k)>are the duals of the LCAO basis, such that:
Eq. CASTEP 66
The spilling parameter varies between one, when the LCAO basis is orthogonal to the PW states and zero when the projected wavefunctions perfectly represent the PW states.
The density operator may be defined:
Eq. CASTEP 67
where nα are the occupancies of the PW states, |χα(k)> are the projected PW states p(k)|ψα(k)> and |χα(k)> are the duals of these states. From this operator the density matrix for the atomic states may be calculated as follows:
Eq. CASTEP 68
The density matrix P(k) and the overlap matrix S(k) are sufficient to perform population analysis of the electronic distribution. In Mulliken analysis (Mulliken, 1955) the charge associated with a given atom, A, is determined by:
Eq. CASTEP 69
and the overlap population between two atoms, A and B, is:
Eq. CASTEP 70
The weight of a band on a given orbital may be calculated by simply projecting the band onto the selected orbital. So, the weight of band α on orbital μ is given by:
Eq. CASTEP 71
The contribution of each band α to the density of states is multiplied by the weight Wαμ(k) to get the projected density of states for orbital μ.
In this case, the natural choice of basis set is that of pseudo-atomic orbitals, generated from the pseudopotentials used in the electronic structure calculation. These are generated automatically when population analysis is performed. The orbitals are calculated by solving for the lowest energy eigenstates of the pseudopotential in a sphere, using a spherical Bessel basis set.
By default, the basis set used on an atom of a given species is that of the orbitals in the closed valence shell of that species. This is usually a sufficiently good basis set for the purposes of population analysis and partial density of states projection. If the spilling parameter indicates that the basis set does not represent the PW states with sufficient accuracy, the number of orbitals in the atomic basis set may be increased. A spilling parameter in the region of a few percent or less is sufficiently good for population analysis and partial density of states projection. Care should be exercised when increasing the basis set dramatically, as this can lead to artificial transfers of charge.
The elements of the overlap matrix S(k) are computed in reciprocal space, as the imposition of periodic boundary conditions implies that the orbitals need only be evaluated on a discrete reciprocal space grid. The overlaps between orbitals on different atomic sites are calculated on the same grid with the application of a phase factor. The overlaps between the plane wave states and the basis functions:
Eq. CASTEP 72
are also calculated in reciprocal space as this is the natural representation of the plane wave states.
The duals of the orbitals are constructed as follows:
Eq. CASTEP 73
The duals of the projected wavefunctions are generated in a similar manner:
Eq. CASTEP 74
where R(k) is the overlap matrix between projected states:
Eq. CASTEP 75
Two auxiliary matrices can be defined for efficiency, namely:
Eq. CASTEP 76
Eq. CASTEP 77
These matrices are used to calculate the spilling parameter:
Eq. CASTEP 78
Similarly, the density matrix may be straightforwardly calculated as:
Eq. CASTEP 79
In the case of spin dependent systems, the populations of the two spin components are calculated separately and the sum and difference of the atomic and overlap populations are calculated to find the net charge and spin, respectively.
It is widely accepted that the absolute magnitude of the atomic charges yielded by population analysis have little physical meaning, since they display a high degree of sensitivity to the atomic basis set with which they were calculated (Davidson and Chakravorty, 1992). However, consideration of their relative values can yield useful information (Segall et al., 1996a; Segall et al., 1996b; Winkler et al., 2001), provided a consistent basis set is used for their calculation.
The spilling parameter should be in the region of a few percent or less in order to ensure that the results are reliable. If this is not the case, the number of orbitals in the basis set should be increased with care. If this does not significantly improve the spilling parameter, it may be that atomic orbitals do not provide a good representation of some of the states in the system.
In addition to providing an objective criterion for bonding between atoms, the overlap population may be used to assess the covalent or ionic nature of a bond. A high value of the bond population indicates a covalent bond, while a low value indicates an ionic interaction. A further measure of ionic character may be obtained from the effective ionic valence, which is defined as the difference between the formal ionic charge and the Mulliken charge on the anion species. A value of zero indicates a perfectly ionic bond, while values greater than zero indicate increasing levels of covalency. See Segall et al. (1996b) for more details.
Care should be taken when interpreting the bond order values for small unit cells. The results will be incorrect whenever an atom is bonded to a number of periodic images of itself. For example, bond orders for the primitive cell of silicon are reported as 3.0, while the same calculation for the conventional cell gives 0.75. The factor of four appears because CASTEP sums up the bond orders for the four bonds that link the Si atom at the origin with the periodic images of itself, in the primitive cell. A manual analysis of the multiplicities is required in such cases, specifically: determine the number of bonds to periodic images of the same atom and divide the bond order by that number.
The physical meaning of Mulliken charges or bond populations for metallic systems is unclear. This calculation is allowed, but the scientific interpretation of the results will be left to the user.