Voronoi Polyhedra

The electron density contains cusps at the nuclei and decays exponentially away from the nuclei. This suggests that the placement of quadrature points should be dependent on the nuclear positions. This can be achieved by partitioning space into regions (Voronoi Polyhedra), each containing only one atom. The total integral is the sum of single atom integrals. There are a number of such partitioning schemes [124,123]; however, the scheme used by the package produced in our research group, Q-CHEM[125], is that developed by Becke [126].

The integral
$$\begin{align}E_{xc}[\rho] = \int F({\bf r}) \, d{\bf r} \end{align}$$
is split up via weight functions $w_{A}({\bf r})$
$$\begin{align}E_{xc}[\rho] = \sum_{A} \int w_{A}({\bf r}) F({\bf r}) \, d{\bf r} \end{align}$$
where the weight functions obey
$$\begin{align}w_{A}({\bf r}) \ge 0 \qquad \text{and} \qquad \sum_{A}w_{A}({\bf r}) = 1. \end{align}$$
These weights are constructed to be almost unity if A is the closest atom, and almost zero in the vicinity of other atoms. This is achieved through the variable $\mu$, defined by
$$\begin{align}\mu_{AB} = \frac{r_{A}-r_{B}}{R_{AB}} \end{align}$$
where r_A and r_B are the distance from atoms A and B, while R_AB is the interatomic distance between atoms A and B. The weight function is then described by
$$\begin{align}w_{A}({\bf r}) &= \frac{P_{A}({\bf r})}{\sum\limits_{A\ne B}P_{B}({\bf r})} \\ P_{A}({\bf r}) &= \prod_{A\ne B}s(\mu_{AB}) \end{align}$$
where the function $s(\mu_{AB})$ is defined as
$$\begin{align}s(\mu_{AB}) = \begin{cases} 0& \text{if $0 < \mu_{AB} \le 1$ } \\ 1& \text{if $-1\le\mu_{AB}<0$ }. \end{cases}\end{align}$$
Becke then removed the discontinuity at $\mu_{AB}=0$ by redefining $s(\mu)$ as a function with
$$\begin{align}s(-1) = 1, \qquad s(+1) = 1, \qquad \text{and} \qquad \left. \frac{ds}{d\mu}\right\vert _{\pm 1} = 0. \end{align}$$

Above is the implementation in Q-CHEM. Becke's original implementation also includes a correction for atomic size. This is accomplished by a change in variable, working with $s(\nu)$ instead of $s(\mu)$, where
$$\begin{align}\nu_{AB} = \mu_{AB} + a_{AB}(1-\mu_{AB}^{2}) \end{align}$$
with a_AB defined by
$$\begin{align}a_{AB} &= \frac{u_{AB}}{u_{AB}^{2}-1} \\ u_{AB} &= \frac{\chi-1}{\chi+1} \\ \chi &= \frac{R_{A}}{R_{B}} \end{align}$$
where R_A and R_B are the Bragg-Slater radii [127,128].

Each of these single-center integrals is then calculated with a spherical polar quadrature grid (requiring the insertion of more weights, w_i), making the final expression for E_xc:
$$\begin{align}E_{xc} = \iint F(r,\theta,\phi)r^{2}\sin\theta \,dr d\theta d\phi = \sum_{i} \sum_{A} w_{A}({\bf r}_{i}) w_{i} F({\bf r}_{i}). \end{align}$$