---
abstract: 'In this work, a method is described to extend the iterative Hirshfeld-I method, generally used for molecules, to periodic systems. The implementation makes use of precalculated pseudo-potential based charge density distributions, and it is shown that high quality results are obtained for both molecules and solids, such as ceria, diamond, and graphite. The use of such grids makes the implementation independent of the solid state or quantum chemical code used for studying the system. The extension described here allows for easy calculation of atomic charges and charge transfer in periodic and bulk systems.'
author:
- 'D. E. P. Vanpoucke'
- 'P. Bultinck'
- 'I. Van Driessche'
bibliography:
- 'danny.bib'
- 'Hirshfeld.bib'
- 'notes.bib'
title: 'Extending Hirshfeld-I to bulk and periodic materials.'
---

Introduction
============

One of the most successful concepts in chemistry is that of “atoms in molecules” (AIM). It states that the properties of a molecule can be seen as simple sums of the properties of its constituent atoms. An impressive amount of insights has been obtained from such a viewpoint, although a precise definition of an AIM remains elusive.[@AyersPW:2000JCP; @BaderRFWMattaCF:2004JPhysChemA; @BaderRFWMattaCF:2006JPhysChemA] All AIM methods have the common purpose of trying to improve our understanding of chemical concepts such as molecular similarity and transferability between molecules.[@ParrRGAyersPWNalewajskiRF:2005JPCA] Since the concept of AIM is basically about how one should divide the electrons, more specifically the charge density distribution (CDD) $\rho(\mathbf{r})$ in the molecule between the different “atoms”, this leads to two obvious categories of approaches in which most of the methods used for defining AIM can be divided. The first category of approaches is based on the wave-function/states of the system, and most of the work is performed in the Hilbert space of the basis functions used. One of the most famous examples here is (probably) the Mulliken approach.[@MullikenRS:1935JChemPhys; @MullikenRS:1955aJChemPhys] The second category of approaches is based on the division of the CDD as it exists in real-space space. In these real-space approaches the molecule is split into atomic basins, that can overlap such as in the Hirshfeld[@HirshfeldFL:1977TCA] and derived methods,[@BultinckP:2007FaradDisc_HirRef; @BultinckP:2007JCP_HirRef; @LillestolenTCWheatleyRJ:2008ChemComm; @LillestolenTCWheatleyRJ:2009JChemPhys; @ManzSholl:JChemTheoryComp2010] or that are non-overlapping such as in Bader’s approach.[@BaderRFW:1991ChemRev]\
The concept of AIM is strongly linked to the concept of transferability. Because both are central in chemistry, and chemists mainly focus on molecules, they are mostly used for molecules.[@FranciscoEPendasAM:2011CTC; @FranciscoEPendasAM:2011TCA] There is, however, no reason why these concepts should not be applicable for periodic systems such as bulk materials. Even more, if these concepts are truly valid, they should hold equally well for solids as for molecules, and should provide additional insight in the chemical properties of defects, such as dopants, interfaces and adsorption of molecules on surfaces.\
In this work, we have implemented an extension of the iterative Hirshfeld-I approach[@BultinckP:2007JCP_HirRef] to periodic systems, as a module in our HIVE-code.[@HIVE_REFERENCE] The implementation makes use of grid stored CDDs, which can easily be generated by standard solid state and quantum chemical codes.\
In Sec. \[sc:Hir\_Methods\] a short review of the parameters used in the solid state and quantum chemical codes is given. Afterwards the basic theory behind the Hirshfeld-I method is presented and extended to periodic systems. In addition, the spatial integration grid, pseudo-potentials and stored CDDs are discussed in view of the Hirshfeld-I method for periodic systems. In Sec. \[sc:Hir\_Res\_and\_Dis\] the influence of the different grids on the accuracy is discussed. We also identify a delocalization problem in the radial CDDs which originates from the plane wave approach and periodic boundary conditions (PBC) used in the solid state code, and a simple and effective solution for this issue is presented. As a last point, the atomic charges in some simple periodic systems are calculated, showing that the algorithm works correctly. Finally, in Sec. \[sc:Hir\_conclusions\] some conclusions are given.\

Methods {#sc:Hir_Methods}
=======

Atomic and molecular calculations {#ssc:Hir_Methods_solstatecodes}
---------------------------------

Hirshfeld-I calculations require CDDs as input. These can be obtained from electronic structure calculations using standard solid state or quantum chemical codes. In this work we have chosen to perform these calculations within the DFT framework using the projector augmented wave (PAW) approach for the core-valence interaction and the local density approximation (LDA) for the exchange-correlation functional as implemented in the <span style="font-variant:small-caps;">VASP</span> code.[@Blochl:prb94; @Kresse:prb99] The kinetic energy cut-off is set at $500$ eV and the $k$-point set is reduced to the $\Gamma$-point for molecular and atomic calculations. For the bulk materials sufficiently converged $k$-point sets were used. To optimize the geometry of the molecules and periodic materials a conjugate gradient algorithm is applied. For molecular calculations only the atom positions are optimized, for bulk materials the cell parameters are optimized simultaneously. All molecules are placed in periodic cells of $20.0\times 21.0\times 20.5$Å$^3$, which provide a sufficiently large vacuum region between periodic copies of the molecules, to prevent interaction.\
Hirshfeld-I data computed using the approach detailed in the present study, are compared to those obtained using more common molecular calculations of AIM properties. For the set of $168$ neutral molecules previously studied by Bultinck *et al.*,[@BultinckP:2007JCP_HirRef; @VanDamme_Bultinck:JCTC2009] geometry optimization and Hirshfeld-I charge calculations are performed at the Local Spin Density Approximation[@HohenbergKohn:PR1964; @KohnSham:PR1965] level with the Slater exchange functional[@SlaterJC:1974] and the VWN$5$ correlation functional[@VoskoWilkNusai:CanJP1980] as implemented in Gaussian-03.[@Gaussian03] Numerical integrations are carried out using Becke’s integration grid with $170$ angular points in the Lebedev-Laikov grid.[@BeckeAD:1987JCP; @LebedevVI_grid:1999DokladyMath] The Hirshfeld-I charges are considered converged if the largest change in charge of any atom is below $0.0005$e.

Hirshfeld methods
-----------------

The basic idea behind the Hirshfeld method,[@HirshfeldFL:1977TCA] also known as the stockholder method, is that the AIM share the charge density in each point of space. This means that the AIM CDD becomes a weighted partition of the molecular CDD. Formally this can be written as: $$\begin{aligned}
\rho_{mol}(\mathbf{r})&=&\sum_{A}{\rho_A^{AIM}(\mathbf{r})}\nonumber\\
\rho_A^{AIM}(\mathbf{r})&=&w_{A}^{H}(\mathbf{r})\rho_{mol}(\mathbf{r}) \qquad \forall \mathbf{r}\\
\mathrm{with\ }\sum_{A}{w_{A}^{H}(\mathbf{r})} &\equiv& 1 \nonumber,\end{aligned}$$ with $\rho_{mol}(\mathbf{r})$ and $\rho_{A}^{AIM}(\mathbf{r})$ the CDDs for the molecule and the AIM. All sums are taken over the entire set of AIM, and  is the Hirshfeld weight function for atom $A$. From these equations the Hirshfeld weight can be written as: $$w_{A}^{H}(\mathbf{r})=\frac{\rho^{AIM}_{A}(\mathbf{r})}{\sum_{B}{\rho_{B}^{AIM}(\mathbf{r})}}.\label{eq:Hir_w_mol}$$ Since the $\rho^{AIM}(\mathbf{r})$ are the CDDs sought, they cannot be used as input. The Hirshfeld method circumvents this problem by using spherically averaged reference state atomic CDDs $\rho_{X}^{\circ}(\mathbf{r})$. In the original paper by Hirshfeld the neutral atomic ground state is used as reference state.[@HirshfeldFL:1977TCA] When summing these isolated atomic CDDs over all AIM, one gets the so-called ‘*promolecular*’ CDD instead of the actual molecular CDD: $$\label{eq:promoldef}
\rho_{promol}(\mathbf{r})=\sum_{B}{\rho^{\circ}_{B}(\mathbf{r})}.$$ It is then assumed that the difference between this promolecular CDD and the actual molecular CDD has only little influence on the Hirshfeld weight . As a result one can write the CDD of an AIM $A$ as: $$\label{eq:Hir_chg_AIM}
\rho^{AIM}_{A}(\mathbf{r})=\frac{\rho^{\circ}_{A}(n_{A}^{\circ},\mathbf{r})}{\rho_{promol}(\mathbf{r})}\rho_{mol}(\mathbf{r})=w^{H}_{A}(\mathbf{r})\rho_{mol}(\mathbf{r}),$$ where the population of the atom $A$ is given by $n_{A}^{\circ}$. In the original Hirshfeld approach, neutral atoms were used as reference. This, however, has been identified by several authors to be a major weakness of the method as changing the choice of the promolecular atom charges can have a highly significant effect on the resulting AIM[@BultinckP:2007JCP_HirRef; @BaderRFWMattaCF:2004JPhysChemA; @BaderRFWMattaCF:2006JPhysChemA; @DavidsonERChakravortyS:1992TheorChimActa; @FranciscoEPendasAM:2007JPhysChemA] From eq. it is easy to understand that the resulting $\rho^{AIM}_{A}(\mathbf{r})$ will tend to be as similar to $\rho^{\circ}_{A}(n_{A}^{\circ},\mathbf{r})$ as possible[@AyersPW:2000JCP; @AyersPW:2002JCP; @BultinckP:2007JCP_HirRef], explaining why the Hirshfeld populations strongly depend on the choice of reference atomic CDDs. Fortunately, this problem can be resolved by using the iterative Hirshfeld-I scheme.[@BultinckP:2007FaradDisc_HirRef; @BultinckP:2007JCP_HirRef] For each iteration $i$, the obtained $\rho^{AIM}_{A}(\mathbf{r})$ are used to calculate the population $n^{i}_{A}$ of each atom $A$. The (spherically symmetric) CDD $\rho^{i}_{A}(n_{A}^{i},\mathbf{r})$ of a free atom $A$ with population $n^{i}_{A}$ is then used as atomic CDD in . For each iteration the new promolecular density $\rho^{i}_{promol}(\mathbf{r})$ is obtained by summing the density distributions $\rho^{i}_{A}(n_{A}^{i},\mathbf{r})$ for all atoms of the molecule. This setup is independent of the initial choice of atomic CDDs and the convergence of the iterative scheme is determined by the convergence of the populations of the AIM.[@BultinckPAyersPWFiasS:2007ChemPhysLett; @GhillemijnDBultinckP:2011JCompChem] Note that the first step in this scheme usually corresponds to the standard Hirshfeld method with $n_{X}^{\circ}=Z_{X}$. At this point it is important to note that the CDDs $\rho^{AIM}_{A}(\mathbf{r})$ and $\rho_{A}^{i}(n_{A}^{i},\mathbf{r})$ will generally be different, despite having the same population. $\rho_{A}^{i}(n_{A}^{i},\mathbf{r})$ is constructed as a spherically symmetric CDD, whereas $\rho_{A}^{AIM}(\mathbf{r})$ is a weighted part of the molecular CDD. The resulting CDD will generally be not spherically symmetric but show protrusions along the directions bonds are formed.\
The extension of the Hirshfeld and Hirshfeld-I methods from molecules to bulk and other periodic materials is quite trivial from the formal perspective. The main problem lies in the fact that a bulk system is considered to consist of an infinite number of “*atoms in the system*”(AIS). Calculating the atomic charge densities for all AIS can, for a periodic system, be reduced to only the atoms in a single unit cell since all periodic copies should yield the same results.\
In addition to this, also the summation limits in Eqs. ,  and change. Where for molecules the sum over $B$ is a finite sum over all AIM, it becomes an infinite sum over all AIS. Because atomic CDDs drop exponentially, the density contribution to the “prosystem” CDD $\rho_{prosys}(\mathbf{r})=\sum^{AIS}_{B}{\rho_B^{\circ}(n_{B}^{\circ},\mathbf{r})}$ of atoms at larger distances becomes negligible. This allows us to truncate the infinite sum to include only the atoms within a certain ‘*sphere of influence*’ (SoI), *i.e.* all atoms of which the contribution to the prosystem CDD is not negligible. Within the iterative Hirshfeld-I scheme we then get: $$\label{eq:Hir_w_iterativeHirshfeldDef}
w^{H,i}_{A}(\mathbf{r}) = \frac{\rho^{i-1}_{A}(n_{A}^{i-1},\mathbf{r})}{\sum_{B}^{SoI}\rho_{B}^{i-1}(n_{B}^{i-1},\mathbf{r})} \quad \forall A\ \epsilon\ \mathrm{unit\ cell}$$ where $$\sum_{A}^{SoI}w^{H,i}_{A}(\mathbf{r})\leq \sum_{A}^{AIS}w^{H,i}_{A}(\mathbf{r})\equiv 1 \qquad \forall \mathbf{r},$$ with $i$ indicating the iteration step, and $\rho^{i}_{A}(n_{A}^{i},\mathbf{r})$ the atomic CDD for an atom $A$ with a population $n_{A}^{i}$ given by $$\label{eq:Hir_pop_int}
n_{A}^{i} = \int{w^{H,i}_A(\mathbf{r})\rho_{system}(\mathbf{r})\mathrm{d}\mathbf{r}},$$ where $\rho_{system}(\mathbf{r})$ is the CDD of the periodic system.

Spatial integration of the population {#ssc:spatInteg}
-------------------------------------

In chemistry, due to the exponential decay of the charge density of atoms and molecules, and due to the sharp cusps present in CDDs at the atomic nuclei, atom centered grids are widely and successfully used. This makes them ideally suited for integrations such as given in Eq. . The multicenter integration scheme proposed by Becke splits up the full space integration into a set of overlapping atom centered spherical integrations.[@BeckeAD:1987JCP] To solve the problem of double counting in the overlapping regions a weight $h_{A}(\mathbf{r})$ is given to each point in space for every atom A in the system, such that $$\sum_{A}^{AIS}{h_{A}(\mathbf{r})}\equiv 1\qquad \forall\mathbf{r}.$$ This weight function indicates how much a point ‘belongs’ to a certain atom $A$. The weight function can be binary, when the space is split up in Voronoi or Wigner Seitz cells,[@VoronoiGeorges:JRAM1907; @WignerSeitzCell:PhysRev1933] or smoothly varying, as is the case in the Becke scheme.[@BeckeAD:1987JCP] As a result, an integrand $F(\mathbf{r})$ can be decomposed as $F(\mathbf{r})=\sum_{A}^{AIS}{h_{A}(\mathbf{r})F_{A}(\mathbf{r})}$ and the full integration becomes $$\label{eq:Becke_IntSplit}
I=\int{F(\mathbf{r})\mathrm{d}\mathbf{r}}=\sum_{A}^{AIS}{\int{h_{A}(\mathbf{r})F_{A}(\mathbf{r})\mathrm{d}\mathbf{r}}},$$ where the sum over all AIS is again an infinite sum. However, in numerical implementations for periodic systems, the exponential decay of the atomic CDD allows us to truncate both the infinite sum and integration region of Eq., without significant loss of numerical accuracy. The sum can be reduced to contain only the atoms included in the SoI of atom $A$ (orange circles in Fig. \[fig:Hir\_SphereOfInfluence\]), because only these atoms contribute significantly to the density in the integration region around atom $A$. In addition, the integration region for all atoms in the SoI can be reduced even further, without loss of accuracy, to only the region that overlaps with the spherical integration region of atom $A$ (blue shaded disc in Fig. \[fig:Hir\_SphereOfInfluence\]).

![(color online) Schematic representation of the integration scheme used for a linear periodic system of atoms (green discs). The black rectangle indicates the unit cell, and the circles indicate the boundaries of the spherical integration regions. The spherical integration regions for the atoms in the sphere of influence (SoI) of atom A are shown in orange, those outside are shown in grey. All integrations can be limited to the blue region (or sections thereof) around atom A (see text). []{data-label="fig:Hir_SphereOfInfluence"}](Hir_SphereOfInfluence.png "fig:"){width="8.0cm"}\

Grid stored charge densities and frozen core pseudo-potentials.
---------------------------------------------------------------

In periodic systems, the use of PBC allows one to reduce the system size dramatically. For bulk materials this even allows simple systems, such as face-centered cubic Cu or Ni, to be represented using single atom unit cells. A useful side effect of such reduced cells is that it is easily possible and relatively cheap to store the CDD $\rho_{system}(\mathbf{r})$ on a three dimensional grid covering the unit cell, and thus fully describing the entire infinite system. The use of such precalculated charge density grids speeds up the Hirshfeld method significantly, since there is no more need to calculate the charge density at any given grid point starting from the wave function of the system. This also makes the implementation independent of the code used to generate the system. The drawback, however, is a slightly reduced accuracy. Since the charge density grid has a finite resolution, interpolation between the stored grid points is needed. This effect is discussed in the following section.\
In chemistry, it is usually sufficient to consider only the valence electrons to describe the interactions between atoms. The core electrons are often considered inert as a result, and are kept frozen during calculations, greatly reducing the computational cost for heavy atoms. This has two small but interesting side effects on the calculated radial charge distributions. Firstly, the integrated charge equals the number of valence electrons only, since the core electrons are not explicitly treated. Secondly, the resulting radial profile is not necessarily monotonically decreasing, showing a minimum or even negative values at the core of the atom (*e.g.* Fig. \[fig:radialdens\]). The origin of this behavior lies in the practical implementation of the PAW pseudo-potentials.[@Blochl:prb94; @Kresse:prb99] By using the same pseudo-potentials for generating the atomic charge distributions and the charge distribution of the system under study, however, no errors are introduced since both systems contain the same frozen core contribution.

Results and Discussion {#sc:Hir_Res_and_Dis}
======================

To this date we are aware of few implementations for periodic systems of the Hirshfeld and approach, but their number is steadily rising. Pendás *et al.*[@PendasAM:2002JChemPhys] investigated Hirshfeld surfaces as approximations for interatomic surfaces for LiF and CS$_2$ crystals. The Cut3D plugin of the <span style="font-variant:small-caps;">ABINIT</span> code can be used to calculate Hirshfeld charges,[@ABINIT_REFERENCE:ComputerPhysCommun2009; @SmithMD:CrystGrowthDes2011; @GlorEC:JSolStateChem2011; @SchierJ:2011ApplMaterInter] and recently Leenaerts *et al.* implemented a “subsystem” based Hirshfeld-I method to study graphane, graphene fluoride and paramagnetic adsorbates on graphene.[@LeenaertsOPeetersF:2010PhysRevB; @LeenaertsO:AplPhysLett2008; @LeenaertsO:PRB2009] In a recent publication Watanabe *et al.*[@WatanabeManzSholl:JPhysChemC2011] presented Hirshfeld results for metal–organic frameworks using the DDEC-code of Manz and Sholl.[@ManzSholl:JChemTheoryComp2010] The same code was probably also used in the investigation of charge injection in graphene layers by Rogers and Liu.[@RogersGWLiuJZ:JACS] These authors present Hirshfeld-I charges, though they do not mention explicitly how these were obtained.\
It is clear that there is a growing interest in codes that can provide atomic populations, however, only few true bulk systems have been investigated using a purely atom based Hirshfeld-I method. For this reason most numerical tests in this work are performed on molecules, though we will investigate the behavior of periodic systems at the end of this section. As a first test system we have chosen the CO molecule. Its small size makes it easily suitable for quick test calculations, and its heteronuclear structure should result in a non-zero charge transfer, at least at equilibrium distance.\
Before proceeding, the different grids used in our current setup of the Hirshfeld-I scheme for periodic systems are introduced. In this setup there are two ‘types’ of CDDs: that of the system and that of the free atoms/ions, which are indicated in the following as the subscripts ‘$sys$’ and ‘$atom$’, respectively. In this, ‘system’ refers to the object of which we want to obtain the atomic charges, and can thus refer to bulk materials, wires, molecules or even single ions. The free atoms/ions on the other hand refer to the single atoms which are used for the generation of the reference radial CDD $\rho_{X}^{i}$ of eqn.. For both types there are two kinds of 3D grids involved:

1.  Linear grids: Instead of using the analytical expression for the underlying wave function, we use the CDDs stored by <span style="font-variant:small-caps;">VASP</span> on a finite numerical grid. These grids span a single unit cell and use uniformly spaced grid points in direct coordinates.[@fn:CHGCAR] (*cf.* Sec.\[ssc:chargedensitygrid\]) In the remainder, the following notation is used:

    -   V$_{atom}$ : The linear grid for the reference atom density distributions as obtained from the atomic calculations.

    -   V$_{sys}$ : The linear grid for the (poly)atomic system under study

2.  Spherical grids: These are atom centered grids which are not limited to a single unit cell. The spherical grids decompose into a radial and a shell grid. In our current setup, a logarithmic grid is used as radial grid, such that closer to the core the grid is sufficiently dense to describe this region accurately. The number of radial points was chosen to equal the numbers suggested by Becke.[@BeckeAD:1987JCP] At each point in this radial grid grid, a shell is located on which grid points are distributed according to a Lebedev-Laikov grid.[@LebedevVI_grid:1999DokladyMath] (*cf.* Sec. \[ssc:sphereIntgrid\]) The total number of grid points $S$ equals the sum over all atoms of the number of radial points ($R_{A}$) used for that atom times the number of points on each shell ($\sigma$): $S=\sum_{A}{R_{A}\cdot\sigma}$. In the remainder only $\sigma$ is varied to study the stability of the integrations. Two three-dimensional spherical grids are distinguished:

    -   S$_{atom}$ : Total spherical grid used to generate the reference spherically averaged radial density distribution for the atoms.

    -   S$_{sys}$ : The multi-center spherical grid for the system under study.

Note that the results of the Hirshfeld-I AIM analysis depend directly on $S_{sys}$ but also indirectly on $S_{atom}$ as this determines the quality of the isolated atomic CCDs.

![The Hirshfeld-I charge for the C atom in a CO molecule as function of the V$_{sys}$ grid resolution. The different curves show the results for the use of different resolutions in the V$_{atom}$ grids, used for generating the atomic radial densities. In all molecular calculations we used spherical integration grids of $1202$ grid points per shell.[]{data-label="fig:HIR_GridConv"}](atdensmodifications3.png "fig:"){width="8cm"}\

[c|ccc]{} &\
& $10$ & $20$ & $30$\
\
$\Delta^{10}_{20}$ & 0.00738 & 0.00743 & 0.00744\
$\Delta^{20}_{30}$ & 0.00126 & 0.00128 & 0.00128\
$\Delta^{30}_{40}$ & 0.00049 & 0.00048 & 0.00049\
&\
& $\Delta^{10}_{20}$ & $\Delta^{20}_{30}$ &\
\
$10$ & -0.00617 & -0.00112 &\
$20$ & -0.00612 & -0.00111 &\
$30$ & -0.00610 & -0.00111 &\
$40$ & -0.00611 & -0.00110 &\

Charge density grids V$_{sys}$ and V$_{atom}$ {#ssc:chargedensitygrid}
---------------------------------------------

As was mentioned in the previous section, the use of stored grid-based CDDs introduces small inaccuracies due to the need for interpolation between the existing grid points. The charge of the C atom in a CO molecule is shown in Fig. \[fig:HIR\_GridConv\] as a function of the grid spacing used in the V$_{sys}$ grid. The different curves are for different grid spacings used in the V$_{atom}$ grids, from which the atomic radial CDD $\rho_{C}^{i}(n_{C}^{i},\mathbf{r})$ and $\rho_{O}^{i}(n_{O}^{i},\mathbf{r})$ are generated. It clearly shows the influence of both grids to be independent, since all curves have the same shape. Looking in detail at the exact numbers reveals that for both grids the same accuracy is obtained (*cf.* Table \[table:gridConv\]). This means that the change in population of the C atom in CO is the same when the same changes are made in either the V$_{atom}$ or the V$_{sys}$ grid; *i.e.* the change in the population (in absolute value) of the C atom is comparable when going from the black to the red curve and when going from a point at $10$ grid points per Å to a point at $20$ grid points per Å on the same curve in Fig. \[fig:HIR\_GridConv\]. Figure \[fig:HIR\_GridConv\] and Table \[table:gridConv\] also show that quite a dense mesh is needed to obtain very accurate results. Though this is not a big problem for periodic systems with small unit cells, it could become problematic for molecules which require big unit cells to accommodate the vacuum required to prevent interaction between the periodic copies. The same is true for high accuracy V$_{atom}$ grids, which are required for high accuracy atomic radial densities. Fortunately, these must only be generated once, and the high resolution radial densities can then be stored in a library.\

![(a) Convergence behavior of the population/charge as function of the number of grid points per spherical integration shell. For the C(black line) and O(red line) atom of the CO molecule the value shown is the absolute value of the difference between the calculated population and the calculated population using the most dense grid. In case of the C atom in the diamond system (green line) the charge should be zero, so the absolute value of the calculated charge is presented. (b) The total number of grid points as function of the number of Lebedev-Laikov grid points per spherical shell.[@LebedevVI_grid:1999DokladyMath] The large number of grid points for the diamond system, is due to its large SoI (see text).[]{data-label="fig:HirI_nlebConv"}](nleb_conv3.png "fig:"){width="8cm"}\

Spherical integration grids S$_{sys}$ and S$_{atom}$ {#ssc:sphereIntgrid}
----------------------------------------------------

Because the populations and charges in the Hirshfeld(-I) approach are obtained by an integration of the CDD of the system, attention needs to be paid to the integration grid used. The current implementation uses a multi-centered grid, as was noted in Sec. \[ssc:spatInteg\]. In this setup we used Lebedev-Laikov grids on the spherical shells of the integration grids.[@LebedevVI_grid:1999DokladyMath] Figure \[fig:HirI\_nlebConv\]a shows the influence of the number of grid points per spherical shell on the accuracy; note the log-log scale used. In general, a denser grid results in a more accurate value for the population. However, different atoms in the same system, and the same calculation can show different convergence, as is shown by the curves of the C and O atoms of a CO molecule. To have the population of the system atoms converged to within $0.001$ electron a few hundred grid points per shell are required. For a molecular calculation the multi-center grid, S$_{sys}$, only requires a few tens of thousands of grid points, as is seen in Fig. \[fig:HirI\_nlebConv\]b. However, a Hirshfeld-I calculation for a bulk material such as diamond, which also has only two atoms in its unit cell, requires a multi-center grid S$_{sys}$ with several tens of millions of grid points. This difference by a factor thousand originates from the fact that the SoI for the diamond unit cell contains a few thousand atoms, which all contribute to the total number of grid points that need to be evaluated. This makes it very important to reduce the SoI to an as small as possible size without significant loss of accuracy.

![image](radialdens_art.png){width="18cm"}\

Atomic radial CDDs
------------------

The atomic CDDs are calculated in a periodic cell under PBC, and stored on $3$D charge density grids V$_{atom}$ (*cf.* Sec. \[ssc:chargedensitygrid\]). The atomic radial CDDs are obtained from these through spherical averaging of the density distributions. Spherical averaging is done using the S$_{atom}$ grid, with spherical shells containing Lebedev-Laikov grids of $5810$ grid points.[@LebedevVI_grid:1999DokladyMath] The resulting distributions for different C ions are shown in Fig. \[fig:radialdens\]a. The populations obtained by spherical integration over these distributions shows that the correct populations are obtained for the neutral and the positively charged ions. The negatively charged ions show a population which is too small. Moreover, the curves in Fig. \[fig:radialdens\]a increase again for longer ranges. This increase is not (only) due to overlapping tails of periodic copies in neighboring cells; extrapolating the decreasing part of the curve, and multiplying with the number of periodic nearest neighbors would give a much smaller value at a distance of $10$Å than is now the case. The origin of the increase lies in the fact that the plane wave approach used for the atom calculations can only bind a limited amount of extra electrons to a given atom. This amount varies from atom to atom. As a result, it tries to place the excess electrons as far from the atom as possible. Due to strong delocalization inherent to plane waves these electrons are spread out over the vacuum between the atoms, with the highest electron density at the center of the unit cell; *i.e.* as far from the atoms as possible.\
It is interesting to note that this artifact is purely due to the use of plane waves, which try to smear out the unbound electrons over the entire (empty) space. If one uses Gaussians instead, the additional electrons are bound by definition through the basis set used, even if these electrons should not be bound to the atom anymore. From the modeling point of view one might prefer this latter type of artifact over the former, since we are interested in the CDD of the electrons for ‘free ions’, irrespective of their bound or unbound nature. Later in this section we show how this delocalization problem can be solved in a simple way.\
Figure \[fig:radialdens\]b shows that this artifact is independent of the resolution used for the V$_{atom}$ grid, as the different curves nicely overlap. Figure \[fig:radialdens\]c shows the influence of the periodic cell size on this artefact. In case of the presented hydrogen ions, it shows that using a cubic periodic cell with a side of $20$Å, the curves for H$^-$ and H$^{2-}$ coincide in the short range region. Moreover, they don’t show the expected exponential decay, and increase sharply at longer range. Using a larger unit cell appears to solve these problems: firstly, the radial distributions for the two ions become distinguishable, and show the expected exponential decay. Secondly, the point where the excess electrons start to interfere noticeably is pushed back to a larger distance.\
This type of behavior is seen for all atom types investigated, positive up to neutral ions give the expected radial distributions, while the negative ions seem hampered by the fact that only a fraction of the additional electrons can be attached to the atom. Fluorine and chlorine are in this respect exceptional since for these atoms also the F$^-$ and Cl$^-$ ions give good distributions and populations. This could be considered a result of the high electron affinity of these elements.\
We find that for most negatively charged ions the populations are too small. As a result, the calculated Hirshfeld weight  for a negatively charged AIM A is underestimated. The easiest way to compensate for this discrepancy is by scaling these specific distributions such that the correct population is found after integration. This way the shape of the curve is maintained, but the resulting weights  increase. To investigate the effects of such a scaling and the erroneous tails shown in Fig. \[fig:radialdens\], we compare the results of three types of atomic radial density distributions. To that end, the reference set of molecules previously used in Hirshfeld-I studies by Bultinck *et al.*[@BultinckP:2007JCP_HirRef; @VanDamme_Bultinck:JCTC2009] is used. This set consists of $168$ neutral molecules containing only H, C, N, O, F, and Cl.\
The first set of atomic radial CDDs, called R1, contains the density distributions as shown in Fig. \[fig:radialdens\]a, where the radial distribution is obtained from a periodic cell of $20\times 20\times 20$Å$^3$.\
For the second set, which is referred to as R2, we have combined these results with the results from a periodic cell of $40\times 40\times 40$Å$^3$ but with a lower grid resolution. In this case, the core part of the radial distribution is taken from the $20\times 20\times 20$Å$^3$ unit cell with the high resolution grid and connected to the tail part obtained from the $40\times 40\times 40$Å$^3$ unit cell. As a result, the high accuracy for the core part of the distribution is maintained, and the tail is corrected through the removal of the delocalized electron contribution. Note that for these distributions, the curves are limited to a distance of $10$Å from the core, i.e. the same maximum radius as is available from the $20\times 20\times 20$Å$^3$ periodic cells. Because the excess tail electrons are not included anymore, the spherically integrated populations of the negative ions are slightly smaller than they are for the R1 set.\
For the third set, R3, the same procedure as for R2 is used, but this time the curves are normalized for the negative ions such that the correct population is given on spherical integration of these radial CDDs.\
To test the accuracy of the results obtained in our periodic implementation using these three different atomic radial density distribution sets, we compare the results for a large benchmark set of molecules with those obtained by a Hirshfeld-I implementation based on a molecular program (*cf.* Sec. \[ssc:Hir\_Methods\_solstatecodes\]).\

![Correlation plots of the charges obtained using the periodic Hirshfeld-I implementation and the standard molecular implementation. (a) using the atomic radial density distributions of set R2 with the corrected tails, and (b) those of set R3 with the corrected tails and correctly normalized distributions (see text). In both cases the VASP optimized geometry was used for the VASP based Hirshfeld-I calculations, while Gaussian optimized geometries were used for the reference calculations.[]{data-label="fig:HirI_CompGvsV"}](VASPvsGaussian2.png "fig:"){width="8cm"}\

First, two sets of calculations are performed to check the influence of the geometry on the results. The first set used the optimized geometries obtained from Gaussian calculations, and the second set used optimized geometries obtained through <span style="font-variant:small-caps;">VASP</span> calculations. For both sets the R1 atomic radial distributions are used. Table \[table:VASPvsGaussianCorr\] shows that the resulting correlations are nearly identical, despite the small differences in geometries for the two sets. Although the correlation coefficients $R^2$ can be considered reasonable, the spread on the C and N data points is quite large. Furthermore, the slope of the linear fit of the C and N data is much too big. In both cases the intercepts are acceptably small.\
Looking at the effect of fixing the tail of the atomic distributions, through the use of the R2 set of atomic radial density distributions, we see a slight increase in the slope for the C and N data, but in general the obtained values for the correlation, intercept and slope remain comparable to the previous sets of calculations. For the results obtained with the R2 radial distributions, Table \[table:VASPvsGaussianCorr\] also shows the standard deviation of the difference between our calculated results and those of the reference data. The values are similar for the first two sets, which are therefore not shown. Table \[table:VASPvsGaussianCorr\] shows the deviation for the C and N data sets to be one order of magnitude larger than for the other elements. It is unclear to the authors why specifically these two elements show such a bad behavior. Looking at the underestimation of the atomic population for the monovalent anions we find an error of roughly $0.5$e for H and N, and $0.3$e for C and O. On the other hand for F and Cl an underestimation of $0.1$ and $0.0$e, respectively, is found. The only aspect in which C and N differ from the other elements is that both can show large positive *and* negative charges.\
The correlation results of the R3 atomic radial density distributions show a clear improvement over the previous results. For each data set the slope is closer to unity, though the intercepts remain as before. Especially the C and N results show a large improvement. Their standard deviation gets halved which is clearly visible in the correlation plots in Fig. \[fig:HirI\_CompGvsV\]a and b. This immediately shows that the simple scaling used for fixing the underestimated population due to the delocalization problem does not introduce large artifacts. In contrast, it actually shows that the shape of the radial CDD of the negatively charged ions is correct.\

  ---- --------- ---------- --------- --------- ---------- --------- --------- ---------- --------- ---------- --------- ---------- --------- ----------
                                                                                                                                              
                                                                                                                                              
          $a$                 $R^2$      $a$                 $R^2$      $a$                 $R^2$    $\sigma$     $a$                 $R^2$    $\sigma$
  H     $1.187$    $-0.014$  $0.940$   $1.208$    $-0.019$  $0.944$   $1.185$    $-0.008$  $0.933$   $0.026$    $1.109$    $-0.021$  $0.979$   $0.016$
  C     $1.353$    $ 0.020$  $0.966$   $1.344$    $ 0.025$  $0.970$   $1.439$    $ 0.020$  $0.963$   $0.181$    $1.131$    $ 0.013$  $0.989$   $0.086$
  N     $1.454$    $ 0.026$  $0.978$   $1.455$    $ 0.018$  $0.979$   $1.483$    $ 0.025$  $0.975$   $0.154$    $1.161$    $ 0.032$  $0.995$   $0.072$
  O     $1.131$    $ 0.066$  $0.992$   $1.123$    $ 0.057$  $0.991$   $1.152$    $ 0.064$  $0.992$   $0.031$    $1.025$    $ 0.048$  $0.993$   $0.019$
  F     $1.110$    $ 0.064$  $0.982$   $1.150$    $ 0.071$  $0.977$   $1.197$    $ 0.070$  $0.980$   $0.011$    $1.162$    $ 0.071$  $0.987$   $0.009$
  Cl    $1.117$    $-0.009$  $0.979$   $1.124$    $-0.018$  $0.980$   $1.110$    $-0.013$  $0.966$   $0.015$    $1.072$    $-0.018$  $0.993$   $0.007$
  ---- --------- ---------- --------- --------- ---------- --------- --------- ---------- --------- ---------- --------- ---------- --------- ----------

![Ball-and-stick representations of (a) the cubic diamond super cell, (b) the graphene sheet (unit cell indicated with black parallelogram), (c) graphite (unit cell indicated) in black, (d) the cubic CeO$_2$ super cell, and (e) the Ce$_2$O$_3$ unit cell. The black, red and yellow spheres indicate the positions of the carbon, oxygen and cerium atoms,respectively. The inequivalent C atoms in the graphite structure are indicated as $A$ and $B$. The $B$ atoms are always located at the center of the hexagons of the neighboring sheets. The inequivalent O positions in Ce$_2$O$_3$ are indicated as $1$ and $2$.\
The ball-and-stick representations are generated using the <span style="font-variant:small-caps;">VESTA</span> visualization tool.[@VESTA:JApplCryst2008][]{data-label="fig:solidsGeom"}](solids_geom.png "fig:"){width="8cm"}\

Periodic systems
----------------

In this final section some actual periodic systems are considered. Their structures are presented in Fig. \[fig:solidsGeom\]: the systems of choice are diamond, graphene, graphite, CeO$_2$, and Ce$_2$O$_3$. The diamond and graphene systems, have the property that all C atoms are equivalent, which should result in zero charges on all atoms. Graphite is quite similar to graphene, however, two inequivalent C positions are present (*cf.* Fig. \[fig:solidsGeom\]c). The ceria systems on the other hand are chosen for the presence of Ce atoms with different valency; tetravalent Ce in CeO$_2$ and trivalent Ce in Ce$_2$O$_3$. Note that for both systems the same Ce pseudo-potential is used. Furthermore, in CeO$_2$, the O ions are all equivalent, while in Ce$_2$O$_3$ the Ce ions are equivalent but the O ions are not, only two of them are equivalent. Furthermore, for Ce$_2$O$_3$ we consider both the ferromagnetic (FM) and anti-ferromagnetic (AF) configurations, allowing us to check how strongly different spin-configurations influence the results in this system.\
For all these systems we use the radial atomic CDDs of the R3 set, with a spherical integration grids S$_{atom}$ containing $5810$ grid points per shell. Table \[table:SolidGrids\] shows the $k$-point sets used for the periodic systems, the number of atoms per unit cell and the number of atoms included in the SoI. The grid point separation for the V$_{sys}$ grid for each of these systems was set to $\leq0.01$Å. The Hirshfeld(-I) populations are calculated using spherical integration grids S$_{sys}$ with $1202$ Lebedev-Laikov grid points per shell.\

  ---------------- ------ ----------- ---------- ----------------- -------- -------
                            atoms per   atoms in    grid points             
                            unit cell    the SoI  ($\times 10^6$)           
  diamond            $21$        $21$       $21$        $2$         $6374$   $128$
  graphene           $21$        $21$        $1$        $2$         $1276$   $22$
  graphite           $21$        $21$       $11$        $4$         $4618$   $104$
  CeO$_2$             $8$         $8$        $8$        $3$         $3063$   $69$
  Ce$_2$O$_3$ FM     $10$        $10$        $5$        $5$         $3009$   $77$
  Ce$_2$O$_3$ AF     $10$        $10$        $5$        $5$         $3025$   $78$
  ---------------- ------ ----------- ---------- ----------------- -------- -------

  : $k$-point sets and the number of atoms in the unit cell and SoI for the periodic systems under investigation. In addition, also the total number of grid points used for the spherical integration grid are given. \[table:SolidGrids\]

  ------------------ ------- ------------ ------------ --
                                                       
                                                       
  diamond            C         $-$0.00007   $-$0.00007 
  \[2mm\] graphene   C            0.00000      0.00000 
  \[2mm\]            C$_A$        0.00113      0.00705 
                     C$_B$     $-$0.00115   $-$0.00707 
  \[2mm\]            Ce           0.59463      2.79393 
                     O         $-$0.30091   $-$1.40056 
  \[2mm\]            Ce           0.49081      2.32119 
                     O$_1$     $-$0.31318   $-$1.61119 
                     O$_2$     $-$0.33576   $-$1.51663 
  \[2mm\]            Ce           0.48575      2.32755 
                     O$_1$     $-$0.30825   $-$1.63062 
                     O$_2$     $-$0.33317   $-$1.51378 
  ------------------ ------- ------------ ------------ --

  : Hirshfeld (H$^0$) and Hirshfeld-I (H$^i$) charges calculated using LDA generated CDDs. The geometries of the periodic systems are shown in Fig. \[fig:solidsGeom\] where the labels for the inequivalent atoms are given.\[table:HirshfeldResults\]

The resulting Hirshfeld (H$^0$) and Hirshfeld-I (H$^i$) charges for the inequivalent atoms are shown in Table \[table:HirshfeldResults\]. It clearly shows the H$^0$ values are closer to zero (i.e. the charge at which the atoms are initialized) than the H$^i$ ones. This is the expected behavior and its origin was discussed earlier by Ayers[@AyersPW:2000JCP] and Bultinck *et al.*[@BultinckP:2007JCP_HirRef] The diamond and graphene charges are (nearly) zero as one would expect based on symmetry arguments. This shows there are no significant artifacts which introduce spurious charges due to the PBC. The results for graphite are somewhat remarkable. Table \[table:HirshfeldResults\] shows there is a small charge transfer going from the A to the B sites. This could be understood as a consequence of the very weak bonding between the A sites in different sheets. For each C atom, three electrons are placed in hybridized $sp^2$ orbitals, where the fourth electron delocalizes in distributed $\pi$ bonds. For the A site C atoms, the contribution to the AIS charge of these $\pi$ bond electrons is shared between the A sites of neighboring sheets. Since the C atom at the B site has no direct neighbor on the neighboring sheet the contribution goes entirely to this atom, resulting in the slightly negatively charged C atom. Charge neutrality results in a slightly positively charged C atom at the A site. Similar behavior was observed by Baranov and Kohout[@KohoutMBaranovA:2011JComputChem] using the Bader approach. These authors, however, find a larger and opposite charge transfer, resulting in a charge of $+0.08$e and $-0.08$e on the C$_{A}$ and C$_{B}$ atoms respectively. This difference could originate from the different methods used.\
The ceria compounds show the behavior expected with regard to equivalent/inequivalent atoms. The Hirshfeld-I values presented in Table \[table:HirshfeldResults\] are comparable with Bader charges presented in literature. Castleton *et al.* found for CeO$_2$ Bader charges of $+2.3$e and $-1.15$e for Ce and O, respectively.[@castleton:jcp2007] The Mulliken atomic charges for Ce$_2$O$_3$ presented in literature appear strongly dependent on the functional used, varying from $+1.29$e for PBE up to $+2.157$e for Hartree-Fock.[@hay:jcp2006; @YangJDolgM:TheorChemAcc2005] This makes it difficult to make a qualitative assessment of our obtained results.\
Table \[table:HirshfeldResults\] also shows there is a clear difference between the tri-and tetra-valent Ce ions, also the different configurations for the O ions show distinctly different charges. Looking at the relative atomic charges of Ce and O atoms in CeO$_2$ and Ce$_2$O$_3$ we find the same relative order as was found by Hay *et al.*[@hay:jcp2006] and a difference in atomic charge for the Ce ions of comparable size. The different charges for the tri-and tetra-valent Ce ions might tempt one to consider these charges as indicators of the oxidation state if not the actual oxidation state of the atoms involved. As a result one could then assume that the same charge in a different configuration would be the result of the same oxidation state (*cf.* concept of transferability). Looking at the charges of the O atoms in both CeO$_2$ and Ce$_2$O$_3$ shows this is clearly not the case, since all O atoms formally have the same oxidation state, while the calculated Hirshfeld-I charges vary $0.2$ electron. The Hirshfeld charges on the other hand show a much smaller variation of only $0.03$ electron. At this point, it is important to stress that atomic charges do not, as opposed to what is often assumed, directly reveal the oxidation state, nor the valence of an atom. A Hirshfeld(-I) analysis, just like a Bader analysis, can only reveal atomic domains. The actual valence of an atom can be derived from the localization indices, which correspond in the simplest form to integrating twice the exchange-correlation density over the same atomic domain.[@BaderRFWStephens:JAmChemSoc1975; @PonecRCooperD:FaradDisc2007] (De)localization indices are obtained from double integration of the exchange correlation density, which requires density matrices, in the simplest case only first order density matrices.[@BultinckPatrick:2010JPCA] Such matrices have been less thoroughly explored in solid state calculations.[@KohoutMBaranovA:2011JComputChem]\
Another interesting point to note is that different spin-configurations have little to no influence on the obtained charges. This is seen when comparing the FM and AF configurations of Ce$_2$O$_3$. This means that for generating the required system CDDs generally non-spin-polarized calculations suffice for the studied systems. Note, however, that the single atom calculations are spin-polarized.

Conclusion {#sc:Hir_conclusions}
==========

We have presented an implementation of the Hirshfeld-I method specifically aimed at periodic systems, such as wires, surfaces, and bulk materials. Instead of calculating the charge densities at each point in space on the fly using the precalculated wave function of the system, we interpolate the charge density from a precalculated CDD on a dense spatial grid, speeding up the calculation of the density significantly. The use of such grids is possible because PBC allow for the use of a relatively small grid to describe the entire system accurately.\
Unlike total energy calculations, the number of atoms involved can not be fully reduced to only those in the unit cell. Although, the populations only need to be calculated for the atoms in the unit cell, the Hirshfeld-I calculations require a large ‘*sphere of influence*’ containing a few thousand of atoms. Using multi-center integration grids, the computational cost can be seriously reduced by cleverly selecting only the grid points involved in the calculations.\
We have shown that the uniform grids used to store both the atom and the system CDDs have an equal influence on the accuracy of the final Hirshfeld-I calculated populations, leading to the suggestion of building the library of atomic radial CDDs using as dense as possible grids. In addition, we have shown that both different atomic types and different chemical environments give rise to a different convergence behavior as function of the spherical integration grid.\
The problems observed for the atomic radial CDDs of negatively charged ions are solved in a simple way, and we show that the introduced scaling of the distributions significantly improves the obtained results for the Hirshfeld-I charges. The resulting values for a benchmark set of $168$ neutral molecules show very good agreement with the values obtained by a previous implementation of the Hirshfeld-I method aimed solely at molecular systems.\
In the final section we have investigated some periodic systems to show the validity of our implementation. For each of these systems the expected behavior of the charges is observed. Because of their simplicity these systems are ideal test cases for Hirshfeld-I implementations for periodic systems.\

Acknowledgement
===============

The research was financially supported by FWO-Vlaanderen, project n$^{\circ}$ $3$G$080209$. This work was carried out using the Stevin Supercomputer Infrastructure at Ghent University, funded by Ghent University, the Hercules Foundation and the Flemish Government–department EWI.
