Charge Analysis: Ver. 1.0

1 Mulliken population

The Kohn-Sham (KS) Bloch functions ${\psi }_{\mu }$ are expanded in a form of linear combination of pseudo-atomic basis functions (LCPAO) ${\varphi }_{i\alpha }$ centered on site ${\tau }_i$ by

${\psi }_{\sigma \mu }^{(𝐤)}(𝐫)={\displaystyle \frac{1}{\sqrt{N}}}{\displaystyle \sum _{\mathrm{n}}^N}{\mathrm{e}}^{\mathrm{i}{𝐑}_{\mathrm{n}}\cdot 𝐤}{\displaystyle \sum _{i\alpha }}{c}_{\sigma \mu ,i\alpha }^{(𝐤)}{\varphi }_{i\alpha }(𝐫-{\tau }_i-{𝐑}_{\mathrm{n}}),$ (1)

where $c$ and $\varphi$ are an expansion coefficient and pseudo-atomic function, ${𝐑}_{\mathrm{n}}$ a lattice vector, $i$ a site index, $\sigma$ ($\uparrow$ or $\downarrow$) spin index, $\alpha \equiv (plm)$ an organized orbital index with a multiplicity index $p$, an angular momentum quantum number $l$, a magnetic quantum number $m$, and $N$ the number of repeated cells. The charge density operator ${\widehat{n}}_{\sigma }$ for the spin index $\sigma$ is given by

${\widehat{n}}_{\sigma }={\displaystyle \frac{1}{V_{\mathrm{B}}}}{\displaystyle {\int }_{\mathrm{B}}}𝑑k^3{\displaystyle \sum _{\mu }^{\mathrm{occ}}}|{\psi }_{\sigma \mu }^{(𝐤)}⟩⟨{\psi }_{\sigma \mu }^{(𝐤)}|,$ (2)

where ${\int }_{\mathrm{B}}$ means the integration over the first Brillouin zone of which volume is $V_{\mathrm{B}}$, and ${\sum }^{\mathrm{occ}}$ means the summation over occupied states. The charge density $n_{\sigma }(𝐫)$ with the spin index $\sigma$ is found as

$n_{\sigma }(𝐫)$ $=$ $⟨𝐫|{\widehat{n}}_{\sigma }|𝐫⟩,$ (3)

$=$ ${\displaystyle \frac{1}{V_{\mathrm{B}}}}{\displaystyle {\int }_{\mathrm{B}}}𝑑k^3{\displaystyle \sum _{\mu }^{\mathrm{occ}}}⟨𝐫|{\psi }_{\sigma \mu }^{(𝐤)}⟩⟨{\psi }_{\sigma \mu }^{(𝐤)}|𝐫⟩,$

$=$ ${\displaystyle \frac{1}{V_{\mathrm{B}}}}{\displaystyle {\int }_{\mathrm{B}}}𝑑k^3{\displaystyle \sum _{\mathrm{n}}^N}{\displaystyle \sum _{i\alpha ,j\beta }}{\displaystyle \sum _{\mu }^{\mathrm{occ}}}{\mathrm{e}}^{\mathrm{i}{𝐑}_{\mathrm{n}}\cdot 𝐤}{c}_{\sigma \mu ,i\alpha }^{(𝐤)*}{c}_{\sigma \mu ,j\beta }^{(𝐤)}{\varphi }_{j\beta }(𝐫-{\tau }_j){\varphi }_{i\alpha }(𝐫-{\tau }_i-{𝐑}_{\mathrm{n}}),$

$=$ ${\displaystyle \sum _{\mathrm{n}}^N}{\displaystyle \sum _{i\alpha ,j\beta }}{\rho }_{\sigma ,i\alpha j\beta }^{({𝐑}_{\mathrm{n}})}{\varphi }_{i\alpha }(𝐫-{\tau }_i){\varphi }_{j\beta }(𝐫-{\tau }_j-{𝐑}_{\mathrm{n}})$

with a density matrix defined by

${\rho }_{\sigma ,i\alpha j\beta }^{({𝐑}_{\mathrm{n}})}={\displaystyle \frac{1}{V_{\mathrm{B}}}}{\displaystyle {\int }_{\mathrm{B}}}𝑑k^3{\displaystyle \sum _{\mu }^{\mathrm{occ}}}{\mathrm{e}}^{\mathrm{i}{𝐑}_{\mathrm{n}}\cdot 𝐤}{c}_{\sigma \mu ,i\alpha }^{(𝐤)*}{c}_{\sigma \mu ,j\beta }^{(𝐤)}.$ (4)

Then, Mulliken populations $M_{\sigma ,i\alpha }$ are given by

$M_{\sigma ,i\alpha }={\displaystyle \sum _{\mathrm{n}}^N}{\displaystyle \sum _{j\beta }}{\rho }_{\sigma ,i\alpha j\beta }^{({𝐑}_{\mathrm{n}})}{S}_{i\alpha ,j\beta }^{({𝐑}_{\mathrm{n}})},$ (5)

where $S_{i\alpha ,j\beta }^{({𝐑}_{\mathrm{n}})}$ is an overlap integral. Since the Mulliken population can be obtained by integrating Eq. (3) over real space, and by decomposing it into each contribution specified with $\sigma$ and $i\alpha$, it can be confirmed that the sum of $M_{\sigma ,i\alpha }$ gives the number of electron $N_{\mathrm{ele}}$ per unit cell as follows:

$N_{\mathrm{ele}}={\displaystyle \sum _{\sigma }}{\displaystyle \sum _{i\alpha }}{M}_{\sigma ,i\alpha }.$ (6)

2 Voronoi charge

The fuzzy cell method decomposes real space into smeared Voronoi cells, called the fuzzy cell [2]. The fuzzy cell at the site $i$ is determined by a weighting function $w_i(𝐫)$:

$w_i(𝐫)={\displaystyle \frac{p_i(𝐫)}{{\sum }_j{p}_j(𝐫)}}$ (7)

with $p_i$ defined by

$p_i(𝐫)={\displaystyle \prod _{j\ne i}}{s}_k({\mu }_{ij}),$ (8)

${\mu }_{ij}={\displaystyle \frac{r_i-r_j}{{\tau }_{ij}}},$ (9)

$r_i=|𝐫-{\tau }_i|,r_j=|𝐫-{\tau }_j|,{\tau }_{ij}=|{\tau }_i-{\tau }_j|,$ (10)

$s_k(x)={\displaystyle \frac{1}{2}}\left\{1-f_k(x)\right\},$ (11)

$f_k(x)=f_0(f_{k-1}(x)),f_0={\displaystyle \frac{3}{2}}x-{\displaystyle \frac{1}{2}}{x}^3,$ (12)

where $k=3$ is chosen in OpenMX. As $k$ increases the fuzzy cells defined by $w$ approach to Voronoi cells (Wigner-Seitz cells). From the definition Eq. (7) it is clear that

${\displaystyle \sum _i}{w}_i(𝐫)=1.$ (13)

Thus, the integration of the charge density Eq. (3) over real space can be decomposed by employing the weighting functions as follows:

${\displaystyle \int 𝑑r^3{n}_{\sigma }(𝐫)}$ $=$ ${\displaystyle \int 𝑑r^3[\sum _i{w}_i(𝐫)]n_{\sigma }(𝐫)},$ (14)

$=$ ${\displaystyle \sum _i}{\displaystyle \int 𝑑r^3{w}_i(𝐫)n_{\sigma }(𝐫)}.$

Thus, the Voronoi charge $N_{\sigma ,i}$ at the site $i$ can be defined by

$N_{\sigma ,i}={\displaystyle \int 𝑑r^3{w}_i(𝐫)n_{\sigma }(𝐫)}.$ (15)

From Eq. (14), it is confirmed that

$N_{\mathrm{ele}}={\displaystyle \sum _{\sigma }}{\displaystyle \sum _i}{N}_{\sigma ,i}.$ (16)

3 Electro-static potential fitting (ESP) charge

Let us consider to express the Hartree potential in a system by the sum of Coulomb potentials with an effective point charge $Q_i$ located on each atomic site ${\tau }_i$ as follows:

$V^{(\mathrm{ESP})}(𝐫)={\displaystyle \sum _{i=1}^{N_{\mathrm{atom}}}}{\displaystyle \frac{Q_i}{|𝐫-{\tau }_i|}},$ (17)

where $N_{\mathrm{atom}}$ is the number of atoms in the system. The $Q_i$ can be found by a least square fitting with a constraint ${\sum }_i{Q}_i=Q_{\mathrm{tot}}$ [3, 4, 5], where $Q_{\mathrm{tot}}$ is the total charge in the system. The Lagrange multiplier method casts this to a minimization problem of the following function $F$:

$F={\displaystyle \sum _{p=1}^{N_p}}{\left(V^{(\mathrm{DFT})}({𝐫}_p)-V^{(\mathrm{ESP})}({𝐫}_p)\right)}^2-\lambda \left(Q_{\mathrm{tot}}-{\displaystyle \sum _i}{Q}_i\right),$ (18)

where $V_p^{(\mathrm{DFT})}$ and $V_p^{(\mathrm{ESP})}$ are the Hartree potential calculated by the DFT calculation and Eq. (17), respectively, $\{{𝐫}_p\}$ is a set of sampling points, and $N_p$ is the number of the sampling points. The sampling points are given by the grids in the real space between two shells of the first and second scale factors times van der Waals radii [6]. The conditions $\frac{\partial F}{\partial Q_i}=0$ and $\frac{\partial F}{\partial \lambda }=0$ lead to

(19)

with

$a_{ij}={\displaystyle \sum _{p=1}^{N_p}}{\displaystyle \frac{1}{|{𝐫}_p-{\tau }_i||{𝐫}_p-{\tau }_j|}}$ (20)

and

$b_i={\displaystyle \sum _{p=1}^{N_p}}{\displaystyle \frac{V^{(\mathrm{DFT})}({𝐫}_p)}{|{𝐫}_p-{\tau }_i|}}.$ (21)

By solving the linear equation Eq. (19), we can find the electro-static potential fitting (ESP) charges. It is noted that the ESP charge is an effective charge on each atom including the contribution of the core charge compared to the Mulliken and Voronoi charges.