IOP PUBLISHING JOURNAL OF PHYSICS: CONDENSED MATTER
J. Phys.: Condens. Matter 21 (2009) 084204 (7pp) doi:10.1088/0953-8984/21/8/084204
A grid-based Bader analysis algorithm
without lattice bias
WTang
1
, E Sanville
2
and G Henkelman
1
1
Department of Chemistry and Biochemistry, The University of Texas at Austin,
Austin, TX 78712-0165, USA
2
Department of Mathematical Sciences, Loughborough University,
Loughborough LE11 3TU, UK
E-mail: henkelman@mail.utexas.edu
Received 19 March 2008, in final form 6 June 2008
Published 30 January 2009
Online at stacks.iop.org/JPhysCM/21/084204
Abstract
A computational method for partitioning a charge density grid into Bader volumes is presented
which is efficient, robust, and scales linearly with the number of grid points. The partitioning
algorithm follows the steepest ascent paths along the charge density gradient from grid point to
grid point until a charge density maximum is reached. In this paper, we describe how accurate
off-lattice ascent paths can be represented with respect to the grid points. This improvement
maintains the efficient linear scaling of an earlier version of the algorithm, and eliminates a
tendency for the Bader surfaces to be aligned along the grid directions. As the algorithm assigns
grid points to charge density maxima, subsequent paths are terminated when they reach
previously assigned grid points. It is this grid-based approach which gives the algorithm its
efficiency, and allows for the analysis of the large grids generated from plane-wave-based
density functional theory calculations.
(Some figures in this article are in colour only in the electronic version)
1. Introduction
First principles methods, and especially density functional
theory (DFT), are commonly used to calculate the many-body
electronic interactions between atoms in molecules and in the
solid state. Accurately describing these complex interactions
is difficult, but it can also be challenging to rationalize the
calculated energetics. One powerful technique for doing this
is to decompose properties of the molecule or material into
contributions from the individual atoms. Bader suggested
an elegant way to do this partitioning [1]. His idea was
to use the charge density to divide space within molecular
systems into atomic (Bader) volumes. Each Bader volume
contains a single charge density maximum, and is separated
from other volumes by surfaces on which the charge density
is a minimum normal to the surface. Typically, there is one
charge density maximum at each atomic center and one Bader
volume for each atom, but this is not required; there are cases
in which Bader volumes do not contain a nucleus [2]. The
dividing surfaces (also called zero-flux surfaces) separating
these volumes lie in the bonding regions between atoms. This
Bader partitioning has an advantage over other partitioning
schemes (e.g. Mulliken population analysis) in that it is based
upon the charge density, which is an observable quantity that
can be measured experimentally or calculated. Furthermore, in
a converged electronic structure calculation, the charge density
is insensitive to the basis set used. In this regard, the Bader
analysis is more robust than wavefunction-based population
methods [3–5].
There are several different approaches to calculating
Bader volumes. Early algorithms were implemented for
quantum chemistry calculations of small molecules, in which
the gradient of the charge density can be calculated from
derivatives of an analytic wavefunction [1, 6, 7]. These
methods first find stationary points in the charge density and
then follow trajectories along the density gradient from these
points to map out their connectivity and the zero-flux dividing
surfaces. With the dividing surfaces represented in this way,
the charge in each Bader volume can be integrated radially
from the charge density maximum to the surface. While this
approach works for small molecules, a high density of descent
trajectories is needed to accurately represent the surface away
from the critical points, and the method has been criticized
as being computationally expensive for large systems [8, 5].
0953-8984/09/084204+07$30.00 © 2009 IOP Publishing Ltd Printed in the UK1
J. Phys.: Condens. Matter 21 (2009) 084204 WTanget al
It was also found to fail for complex bonding geometries and
when the radial integration rays have multiple intersection
points with the dividing surface [9–11].
Several improvements to this original approach have
been proposed. Popelier developed a method to more
accurately integrate the charge density within each Bader
volume using a Fourier–Chebyshev fit [12] and a bisection
method to analytically represent the dividing surfaces [13].
Other approaches suggested by Popelier include the use
of the divergence theorem to replace the three-dimensional
integration over the Bader volumes into a two-dimensional
integral over the dividing surfaces [8], and the use of a tree
search to treat complex bonding topologies between critical
points [14, 15].
Most current implementations of Bader’s analysis are bas-
ed upon a grid of charge density values [6, 16, 17, 11, 18, 19].
This is particularly important for plane-wave-based DFT
calculations, because it allows for the analysis of condensed
phase systems with many atoms. The algorithm described
in this work is an improvement upon such a grid-based
method [18]. In that original work, an algorithm was
introduced in which ascent trajectories along the charge density
were followed between grid points to determine the Bader
volumes. The constraint of ascent trajectories to the grid
means that each point need be considered only once. The
method is highly efficient, scales linearly with system size,
and is robust to complex bonding topology found in condensed
systems. Recently, however, Sanville et al showed that
this algorithm introduces a bias, causing the Bader surfaces
to artificially follow the orientation of the lattice [19]. A
modified algorithm was implemented by these authors in which
ascent trajectories along interpolated gradients of the charge
density were not constrained to the grid. This algorithm
removes the bias and retains the linear scaling of the original
method. Here, we present a modified version of our original
algorithm which also removes the lattice bias in a somewhat
different way. Trajectories are constrained to grid points at
which charge density gradients are evaluated, and a correction
vector is calculated which points from the nearest grid point
to the unbiased (off-lattice) trajectory. In this way, lattice
bias is removed and we have an algorithm which retains the
efficiency, linear scaling, and robustness of the original grid-
based method.
2. Methodology
This method for finding Bader volumes without lattice bias is
an improvement to an earlier grid-based method [18]. We first
describe the original ‘on-grid’ method, the bias that it has, and
then this new ‘near-grid’ modification which removes the bias.
The input to these methods is a grid of charge density
values defined on an orthogonal lattice. The generalization
to non-orthogonal lattices and specific boundary conditions
follows without significant complication, and has been done
in the software available from [20]. The charge density grid
can be obtained from ab initio calculations or from experiment.
The method described here is particularly well suited to DFT
calculations of large molecules or materials, with many atoms
and complex bonding geometries.
2.1. On-grid method
The following steps summarize the on-grid Bader analysis
method. Additional details are provided in [18].
(i) First, an initial grid point is chosen. To associate this point
with a Bader volume, a path of steepest ascent is followed
between neighboring grid points along the charge density
gradient.
(ii) From each grid point along the path,
(i, j, k),the
projection of the charge density gradient is calculated
along the direction to each of the 26 neighboring grid
points,
∇ρ(i, j, k) ·ˆr(di, dj, dk) =
ρ
|ˆr|
.
(1)
Here,
(di, dj, dk) is a vector of integers describing the
step along the grid to the neighbor. The integers
di, dj,
and
dk can each take the values {−1, 0, 1}, excluding
di = dj = dk = 0. There are 26 neighbors to each grid
point. The vector
ˆr(di, dj, dk) is the normalized direction
to the neighbor reached by the grid step
(di, dj, dk).The
gradient projections in equation (1) are calculated using
finite difference, where the change in charge density to the
neighbor is
ρ = ρ(i + di, j +dj, k + dk) −ρ(i, j, k), (2)
and the distance to the neighbor is
|r|=|r(i + di, j + dj, k + dk) −r(i, j, k)|. (3)
(iii) One of the 26 neighbors is then determined as the next
point along the ascent path. This neighbor,
(i + di, j +
dj, k + dk)
, is the one which maximizes the gradient
projection from equation (1). This gradient must be a
positive value in order to make the step. If there are
no positive values, the point
(i, j, k) is a charge density
maximum.
(iv) The steepest ascent path is followed until a charge density
maximum is found. Each new maximum found is assigned
an integer value corresponding to the order in which it was
found. In an array with the same dimensions as the charge
density grid (initialized to zero), we take the number of the
maximum to which the trajectory terminated, and assign
that integer to all the grid points along the trajectory. We
can do this because an ascent trajectory from each point
along the path terminates at that same maximum. In this
way, all points along the path are associated with a Bader
volume around the charge density maximum.
(v) Each grid point is analyzed in this way. The order in which
grid points are analyzed is not particularly important, and
it is easiest to cycle through them in a loop over all values
of
i, j,andk. Grid points that are already assigned to a
Bader volume are skipped, so that ascent trajectories are
only followed from unassigned grid points. Trajectories
are terminated if they either reach a grid point which
has already been assigned, or when a new charge density
maximum is found. Figure 1 illustrates both kinds of
trajectories. If the trajectory terminates at an assigned grid
2
J. Phys.: Condens. Matter 21 (2009) 084204 WTanget al
Figure 1. An illustration of the steepest ascent paths (a) on a charge
density grid to find the Bader volumes using the on-grid analysis
method. These ascent trajectories are constrained to the grid points,
moving at each step to the neighboring grid point towards which the
charge density gradient is maximized. Each trajectory either
terminates at a new charge density maximum,
m
i
, or at a grid point
which has already been assigned. After all grid points are assigned
(b), the set of points which terminate at each maximum (green to
m
1
and blue to m
2
) constitute that Bader volume. The Bader surfaces
(red curved line) separate the volumes.
point, the initial point and all points along the trajectory
are assigned to that same Bader volume. If the trajectory
terminates at a new charge density maximum, a new entry
is made in the list of known charge density maxima, and all
points along the trajectory are assigned to that new Bader
volume number.
(vi) After analyzing all grid points in this way, each is assigned
to a Bader volume. The total charge in each Bader region
is found by integrating the charge density over the grid
points assigned to that region. The surface around a Bader
volume can be visualized by plotting the charge density of
that individual Bader volume.
2.2. Lattice bias in the on-grid method
The on-grid method is simple to implement, efficient, and
robust. Recently, however, Sanville et al showed that the
method results in some systematic error of the Bader surfaces
and population analysis [19]. This error is caused by the
fact that the ascent trajectories in the on-grid method are
constrained to the grid points. With trajectories following grid
directions, the Bader surfaces become artificially angular with
facets parallel to the grid. This bias remains in the limit of a
fine charge density grid.
The lattice bias in the on-grid method is illustrated in
figure 2. Ascent trajectories step between grid points in the
direction that is most aligned with the charge density gradient
(see equation (1)). Since there are only 26 such directions (8
in two-dimensions), the trajectories accumulate error when the
gradients lines do not run parallel to the lattice. In figure 2,an
on-grid trajectory is shown crossing the true dividing surface,
resulting in points (e.g. the light-blue point) that are assigned to
the incorrect Bader volume, and a surface which is artificially
aligned along the lattice.
2.3. Near-grid method without lattice bias
Here, we describe a method which fixes the lattice bias error
in the on-grid method. This new method is described as the
Figure 2. Illustration of lattice bias in the on-grid method. The true
dividing surface (red) runs parallel to the gradient lines, but the
on-grid ascent trajectories follow the lattice direction along which the
projection of the charge density gradient is maximized. Starting from
the initial point, this direction is along
+x and the trajectory moves
from grid point to grid point in this direction. The error can be seen
as the on-grid trajectory (straight blue arrows through the light-blue
point) deviates from the true trajectory (solid blue curved arrow).
The resulting dividing surface follows the
x lattice direction instead
of the true dividing surface.
near-grid method; trajectories still follow from grid point to
grid point, but a correction vector is remembered, which points
to the location of the actual trajectory as it passes near each
grid point. The near-grid method eliminates the lattice bias
because trajectories are not constrained to grid points. It also
retains the simplicity and linear scaling with system size of the
on-grid method.
The following steps describe the near-grid algorithm.
(i) Starting from an initial grid point,
(i, j, k), the trajectory
of steepest ascent is followed up the charge density. Each
step is made to one of the 26 neighboring grid points.
(ii) In order to make a step from point
(i, j, k), the charge
density gradient,
∇ρ, is calculated from the charge density
at the six closest neighbors using a central finite difference
scheme. The components of the charge density in the x,
y,andz directions are:
∇ρ
x
=
ρ(i +
1, j, k) −ρ(i −1, j, k)
|r(i + 1, j, k) −r(i −1, j, k)|
,
(4a)
∇ρ
y
=
ρ(i, j +
1, k) −ρ(i, j − 1, k)
|r(i, j +1, k) −r(i, j −1, k)|
,
(4b)
∇ρ
z
=
ρ(i, j, k +
1) −ρ(i, j, k − 1)
|r(i, j, k + 1) −r(i, j, k − 1)|
.
(4c)
(iii) To follow the gradient, a step along the gradient vector,
r
grad
= c(∇ρ
x
, ∇ρ
y
, ∇ρ
z
), (5)
should be taken. For this step to advance the trajectory by
one grid point, and no more in any direction, the constant
is chosen as
c = min(dx/|∇ρ
x
|, dy/|∇ρ
y
|, dz/|∇ρ
z
|), (6)
3
J. Phys.: Condens. Matter 21 (2009) 084204 WTanget al
Figure 3. An illustration of the near-grid step and how a correction
vector points to the true trajectory. Steps are taken between grid
points. The first step between grid points,
r
grid
and along the gradient,
r
grad
, differ by the correction vector r
1
. After the second step, the
difference is
r
2
, and the total correction vector is r = r
1
+r
2
.
Now, the
y component of the correction vector is larger than half of
the grid spacing so that a correction step is taken in the
−y direction,
and the correction vector is recalculated as
r
2
.
where dx,dy,anddz are the grid spacings along the x,
y,andz Cartesian directions, respectively. The step r
grad
takes the trajectory (in general) off the grid. To retain
the efficiency of the on-grid method, we then jump to the
nearest grid point, so that the trajectory can be described
as a hopping process between grid points. If we label this
step between grid points as
r
grid
, we can then describe a
correction vector (initially zero) from the final grid point
to the location of the true trajectory,
r = r +
r
grad
−r
grid
.
(7)
In figure 3 the first
r
grid
step is along +x, and the correction
vector,
r
1
, points along −y to the true trajectory.
(iv) From each new point along the ascent trajectory, the
vectors
r
grad
and r
grid
are calculated, and the correction
vector is accumulated so that it always points to the true
trajectory. When the length of any component of
r
is larger than half of the grid spacing, a correction step
is taken in that direction. The correction vector is then
recalculated by subtracting the correction step. In this
way, the true trajectory is never more than half a grid point
from the current grid point in any direction.
(v) The ascent trajectory is terminated when one of two
criteria is met: (i) when it reaches a charge density
maximum, or (ii) when it reaches a point for which the
point itself and all of its neighbors are assigned to the same
Bader region. In either case, all the grid points along the
ascent trajectory are assigned to the same Bader region as
the end points of the trajectory.
(vi) An ascent trajectory is calculated from each grid point
until all are assigned to a Bader volume. From each new
trajectory, the correction vector
r is reset to (0, 0, 0).
dividing surfaces
(a) on-grid
method
(b) near-grid, before
edge-refinement
(c) near-grid, after
edge-refinement
error regions
Figure 4. Comparison of the on-grid and near-grid methods for a
two-dimensional charge density surface constructed from three
Gaussian functions. By construction, the true dividing surface (gray
lines in (a)) are at a small angle to the lattice. The on-grid method
results in biased dividing surfaces that follow the grid; error regions,
in which the Bader volume is incorrectly assigned, are colored white.
In a single iteration of the near-grid method (b), the error regions are
confined to the Bader surfaces. A subsequent refinement iteration (c)
corrects all but two grid points, which are mis-assigned because of
the low resolution of the grid. As the grid density is increased, the
near-grid method finds the exact Bader volumes.
(vii) When all the points are assigned, a final refinement of
the grid points adjacent to the Bader surfaces is required.
Depending upon the order in which grid points are
analyzed, the grid point adjacent to the Bader surface can
be assigned to one of the two volumes on either side of
the dividing surface. This ambiguity is due to the fact
that the trajectory between grid points deviates from the
true trajectory by up to half a grid step. The refinement
step corrects this by first identifying all grid points on
the boundary of a Bader volume. An ascent trajectory
is followed from each of these boundary grid points to
determine in which volume it belongs. In this case, only
the initial point is assigned, so that (in the case of smooth
charge density grids) the refinement process need not be
repeated.
3. Results
3.1. Two-dimensional model
The lattice bias in the on-grid method is particularly severe
when the true dividing surface runs at a small angle to the
grid over long distances. In figure 4 we have constructed
such a two-dimensional charge density grid from the sum of
three Gaussian functions. The on-grid method finds vertical
dividing surfaces that deviate from the true (gray) dividing
surface, resulting in the (white) error regions. In these
regions the method has assigned charge to the wrong Bader
region.
The near-grid method corrects this lattice bias. After
ascent trajectories are followed from each grid point, error
regions are localized to grid points adjacent to the true
dividing surface. A subsequent refinement step eliminates
these errors, and gives the correct Bader volumes in the limit
of a sufficiently fine grid.
4
J. Phys.: Condens. Matter 21 (2009) 084204 WTanget al
O
H
H
O
H
H
(a) on-grid method (b) near-grid method
Figure 5. Comparison of Bader dividing surfaces found with the
on-grid method (left) and the near-grid method (right). The on-grid
surfaces are angular with facets oriented along the grid directions.
This bias is removed in the near-grid method.
3.2. A water molecule
Finding the Bader volumes in a water molecule is a more
realistic test. The geometry and charge density of the molecule
was calculated with Gaussian 98 [21] using the aug-cc-pVDZ
basis set at the MP2 level of perturbation theory. The charge
density was written on a 257
3
orthogonal grid. The dividing
surfaces found with the on-grid and near-grid methods are
shown in figure 5. The lattice bias is apparent in the on-
grid method from the angular shape of the surfaces (see
also [18]). This bias is removed with the near-grid method and
the surfaces become smooth. The ripples in the surface shown
in figure 5 are due to the finite resolution of the grid.
3.3. Ionic charge in a NaCl crystal
These grid-based analysis methods can be used for solid-
state systems as well as molecular ones. In order to show
convergence of the method with respect to grid density, we
have calculated the Bader charges in a NaCl salt crystal,
using the eight-atom unit cell illustrated in figure 6 (inset).
The charge density was generated using the VASP plane-
wave-based DFT code [22], using the PW91 generalized
gradient functional [23]. Frozen core charges were represented
with pseudopotentials of the Vanderbilt form [24]usingthe
projector augmented wave (PAW) framework [25]. Including
this frozen core charge in the charge density grid is important
for calculating accurate Bader charges. A plane wave energy
cutoff of 262.5 eV was used, and the Brillouin zone was
sampled with a 3
×3 ×3 Monkhorst–Pack [26] k-point mesh.
An optimal lattice constant of 5.86
˚
A was determined for the
NaCl crystal.
A set of charge density grids ranging from 60
3
points to
350
3
points were calculated. The grid spacing is 0.095
˚
Ainthe
former case and 0.016
˚
A in the latter. The calculated valance
charge on the Na ions is shown in figure 6. The near-grid
method shows monotonic and smooth convergence to a valance
charge of 0.828
e, whereas the on-grid method has systematic
errors in the limit of a fine grid. Although the systematic error
of less than 0.01 e in the on-grid method is acceptable for
most calculations, the near-grid method should be used for its
improved accuracy and systematic convergence.
Figure 6. Convergence with respect to grid density of the Bader
charges for ions in a NaCl crystal. A fully converged charge on each
Na ion is 0.828
e (blue dashed line). The (old) on-grid method
deviates from this value by ca 0.01
e for a fine grid with 40 million
points. The (new) near-grid method rapidly converges to the
correct value.
Figure 7. The Bader surface in H
2
O is strongly dependent upon the
orientation of the molecule with the on-grid method. This orientation
dependence is due to the bias in the method which tends to orient the
Bader surfaces along the grid directions. In the near-grid method,
this bias is removed and the Bader surfaces are insensitive to the
orientation of the molecule.
3.4. Variation of charge with molecular orientation
The lattice bias in the on-grid method can be seen by
comparing the Bader volumes and charges for a molecule
aligned at different orientations with respect to the charge
density lattice. A converged calculation should not depend
upon this orientation. For this test, we have calculated the
Bader volumes and charges around an H
2
O molecule using the
VASP code. The details of these calculations are the same as
for the NaCl calculations with the exception of a 250 eV energy
cutoff and a single
-point for the isolated molecule.
Figure 7 shows a change in the shape of the Bader
volumes for the on-grid method as the H
2
O molecule is
rotated by 45
◦
in the plane of the figure. With the
5
