A fast and robust algorithm for Bader decomposition of charge density

Abstract

An algorithm is presented for carrying out decomposition of electronic charge density into atomic contributions. As suggested by Bader [R. Bader, Atoms in Molecules: A Quantum Theory, Oxford University Press, New York, 1990], space is divided up into atomic regions where the dividing surfaces are at a minimum in the charge density, i.e. the gradient of the charge density is zero along the surface normal. Instead of explicitly finding and representing the dividing surfaces, which is a challenging task, our algorithm assigns each point on a regular (x, y, z) grid to one of the regions by following a steepest ascent path on the grid. The computational work required to analyze a given charge density grid is approximately 50 arithmetic operations per grid point. The work scales linearly with the number of grid points and is essentially independent of the number of atoms in the system. The algorithm is robust and insensitive to the topology of molecular bonding. In addition to two test problems involving a water molecule and NaCl crystal, the algorithm has been used to estimate the electrical activity of a cluster of boron atoms in a silicon crystal. The highly stable three-atom boron cluster, B₃I is found to have a charge of −1.5 e, which suggests approximately 50% reduction in electrical activity as compared with three substitutional boron atoms.

Introduction

The properties of chemicals and materials are often described in terms of charge transfer between atoms and the presence of ionic charges or electric multipoles on atoms or molecules. Theoretical calculations producing estimates of the electronic charge distribution in the system can, in principle, provide this type of information but it is not clear how to extract it. Atomic charges in molecules or solids are not observables and, therefore, not defined by quantum mechanical theory. The output of quantum mechanical calculations is a continuous electronic charge density and it is not clear how one should partition electrons amongst fragments of the system such as atoms or molecules. Many different schemes have been proposed, some based on electronic orbitals and others based on only the charge density (see for example Ref. [1]).

The most commonly used orbital based method is the so-called Mulliken analysis. It can be applied when basis functions centered on atoms are used in the calculation of the electronic wavefunction of the system. The charge associated with the basis functions centered on a particular atom is then assigned to that atom. This can be a fast and useful way of determining partial charges on atoms but it has the major drawback that the analysis is sensitive to the choice of basis set. Furthermore, in the limit of an infinite basis, the assignment becomes arbitrary. As an extreme case, a complete basis set can be formed where all the basis functions are centered on just one atom in the system. The Mulliken analysis would then clearly be misleading as it would assign all the electrons in the system to that one atom. Moreover, many calculations are carried out with plane wave basis functions which are not associated with any particular atom in the system. Mulliken analysis is not applicable in such cases.

A different approach is to focus entirely on the charge density as has been proposed by Bader [2]. Space is then divided into regions by surfaces that run through minima in the charge density. More precisely, at a point on a dividing surface the gradient of the electron density has no component normal to the surface. We will refer to regions bounded by such dividing surfaces as Bader regions. Because this analysis is based solely on the charge density, it is rather insensitive to the basis set used in the electron wavefunction calculation and can be used to analyze plane wave based calculations as well as calculations using atomic basis functions. Furthermore, it can be used to analyze electronic density obtained by other means, for example from X-ray crystallography. Each Bader region often contains one nucleus, but this is not necessarily so, sometimes no nucleus is found within a Bader region. By integrating the electronic density within the Bader region where an atom’s nucleus is located, and possibly adding the electronic charge in nearby regions that do not include a nucleus, the total charge on an atom can be estimated.

A common complaint about Bader analysis is the computational effort and complexity of the algorithms that have been developed [3], [1]. Commonly used implementations [2], [4], [5] involve finding the critical points of the charge density where ∇ρ = 0, followed by the construction of the zero-flux surfaces which intersect these points and then integration of the electronic density within each region. Several refinements have been made since the initial implementation of the method. Popelier improved the basin integration method by using Chebyshev polynomials [6], along with a bisection method [7] for finding dividing surfaces and increased integration accuracy around critical points [2]. Popelier has also suggested integrating the electric field on the dividing surface around atoms in order to find the enclosed charge through the divergence theorem [3]. Most recently Malcolm and Popelier [8], [9] have used dynamic grid techniques introduced by Silvi and coworkers [10], [11] in order to effectively treat complicated bonding topology. A very different approach was presented by Uberuaga et al. where multiple replicas of the system connected by springs form a discrete representation of an ‘elastic sheet’ that stretches along the dividing surface surrounding a given Bader region [12].

In this paper, a simple and fast method for carrying out Bader decomposition of electron density is presented. It is different from the methods mentioned above in that an explicit representation of the dividing surfaces is not used and no attempt is made to locate stationary points of the charge density. This makes the algorithm more robust. Algorithms that are commonly used today are known to have convergence problems in some cases (an example involving clusters of water molecules is given in Ref. [12]). As in several previous implementations of Bader analysis, such as TopMoD [10], InteGriTy [13], Extreme 94 [14], and MORPHY [4], it is assumed that the electron density is represented by its values on a regular grid in three-dimensions, but our method is different from previous methods in that only steepest ascent trajectories confined to the grid points are used to identify the Bader regions. The algorithm is robust since dividing surfaces and critical points do not need to be found explicitly. This makes the method applicable to systems with complicated bonding topology where the representation of the dividing surfaces is difficult. Although the dividing surfaces are not found explicitly, it is easy to render them for visualization after the analysis is complete. Bader regions which do not contain a nucleus [15] do not pose a problem because the method does not use the nuclear coordinates, neither in finding Bader regions nor in the integration of the electron density. Furthermore, no assumption about the shape of the Bader region is made in the integration of the electron density. Such an assumption is, for example, made when radial integration away from nuclei is used (see for example [13]). Such a procedure is valid when radial rays do not cross the dividing surface more than once but this is not always the case [16]. The number of arithmetic operations per grid point in our algorithm is insensitive to the size of the grid and the number of atoms in the system so the computational effort scales linearly with the number of grid points used in the representation of the electron density. A description of the algorithm will now be presented.