10.2.1 Population Analysis

The one-electron charge density, usually written as

$$\rho (𝐫)=\sum _{\mu \nu }{P}_{\mu \nu }{\varphi }_{\mu }(𝐫){\varphi }_{\nu }(𝐫)$$

(10.1)

represents the probability of finding an electron at the point $𝐫$, but implies little regarding the number of electrons associated with a given nucleus in a molecule. However, since the number of electrons $N$ is related to the occupied orbitals ${\psi }_i$ by

$$N=2\sum _a^{N/2}{\left|{\psi }_a(𝐫)\right|}^2$$

(10.2)

We can substitute the atomic orbital (AO) basis expansion of ${\psi }_a$ into Eq. (10.2) to obtain

$$N=\sum _{\mu \nu }{P}_{\mu \nu }{S}_{\mu \nu }=\sum _{\mu }{\left(\mathrm{𝐏𝐒}\right)}_{\mu \mu }=\mathrm{Tr}(\mathrm{𝐏𝐒})$$

(10.3)

where we interpret ${(\mathrm{𝐏𝐒})}_{\mu \mu }$ as the number of electrons associated with ${\varphi }_{\mu }$. If the basis functions are atom-centered, the number of electrons associated with a given atom can be obtained by summing over all the basis functions. This leads to the Mulliken formula for the net charge of the atom $A$:

$$q_A=Z_A-\sum _{\mu \in A}{\left(\mathrm{𝐏𝐒}\right)}_{\mu \mu }$$

(10.4)

where $Z_A$ is the atom’s nuclear charge. This is called a Mulliken population analysis.⁹⁴² Q-Chem performs a Mulliken population analysis by default.

Although conceptually simple, Mulliken population analyses suffer from a heavy dependence on the basis set used, as well as the possibility of producing unphysical negative numbers of electrons. An alternative is that of Löwdin population analysis,⁶⁰⁸ which uses the Löwdin symmetrically orthogonalized basis set (which is still atom-tagged) to assign the electron density. This shows a reduced basis set dependence, but maintains the same essential features.

While Mulliken and Löwdin population analyses are commonly employed, and can be used to produce information about changes in electron density and also localized spin polarizations, they should not be interpreted as oxidation states of the atoms in the system. For such information we would recommend a bonding analysis technique (LOBA or NBO).

A more stable alternative to Mulliken or Löwdin charges are the so-called “ChElPG” charges (“Charges from the Electrostatic Potential on a Grid”).¹⁰³ These are the atom-centered charges that provide the best fit to the molecular electrostatic potential, evaluated on a real-space grid outside of the van der Waals region and subject to the constraint that the sum of the ChElPG charges must equal the molecular charge. Q-Chem’s implementation of the ChElPG algorithm differs slightly from the one originally algorithm described by Breneman and Wiberg,¹⁰³ in that Q-Chem weights the grid points with a smoothing function to ensure that the ChElPG charges vary continuously as the nuclei are displaced.³⁷⁵ (For any particular geometry, however, numerical values of the charges are quite similar to those obtained using the original algorithm.) Note, however, that the Breneman-Wiberg approach uses a Cartesian grid and becomes expensive for large systems. (It is especially expensive when ChElPG charges are used in QM/MM-Ewald calculations.³⁹⁴) For that reason, an alternative procedure based on atom-centered Lebedev grids is also available,³⁹⁴ which provides very similar charges using far fewer grid points. In order to use the Lebedev grid implementation the $rem variables CHELPG_H and CHELPG_HA must be set, which specify the number of Lebedev grid points for the hydrogen atoms and the heavy atoms, respectively.

Hirshfeld population analysis³⁸⁸ provides yet another definition of atomic charges in molecules via a Stockholder prescription. The charge on atom $A$, $q_A$, is defined by

$$q_A=Z_A-\int 𝑑𝐫\frac{{\rho }_A^0(𝐫)}{{\sum }_B{\rho }_B^0(𝐫)}\rho (𝐫),$$

(10.5)

where $Z_A$ is the nuclear charge of $A$, ${\rho }_B^0$ is the isolated ground-state atomic density of atom $B$, and $\rho$ is the molecular density. The sum goes over all atoms in the molecule. Thus computing Hirshfeld charges requires a self-consistent calculation of the isolated atomic densities (the promolecule) as well as the total molecule. Prior to the SCF calculation, the Hirshfeld atomic density matrix is constructed. After SCF convergence, numerical quadrature is used to evaluate the integral in Eq. (10.5). Neutral ground-state atoms are used, as the choice of appropriate reference for a charged molecule is ambiguous (such jobs will crash). As numerical integration (with default quadrature grid) is used, charges may not sum precisely to zero. A larger XC_GRID may be used to improve the accuracy of the integration, but the magnitude of the Hirshfeld charges should be largely independent of grid choice.

The charges (and corresponding molecular dipole moments) obtained using Hirshfeld charges are typically underestimated as compared to other charge schemes or experimental data. To correct this, Marenich et al. introduced “Charge Model 5” (CM5),⁶³⁴ which employs a single set of parameters to map the Hirshfeld charges onto a more reasonable representation of the electrostatic potential. CM5 charges generally lead to more accurate dipole moments as compared to the original Hirshfeld charges, at negligible additional cost. CM5 is available for molecules composed of elements H–Ca, Zn, Ge–Br, and I.

The use of neutral ground-state atoms to define the promolecular density in Hirshfeld scheme has no strict theoretical basis and there is no unique way to construct the promolecular densities. For example, Li${}^0$F${}^0$, Li${}^{+}$F${}^{-}$, or Li${}^{-}$F${}^{+}$ could each be used to construct the promolecular densities for LiF. Furthermore, the choice of appropriate reference for a charged molecule is ambiguous, and for this reason Hirshfeld analysis is disabled in Q-Chem for any molecule with a net charge. A solution for charged molecules is to use the iterative “Hirshfeld-I” partitioning scheme proposed by Bultinck et al.,^{110, 990} in which the reference state is not predefined but rather determined self-consistently, thus eliminating the arbitrariness. The final self-consistent reference state for Hirshfeld-I partitioning usually consists of non-integer atomic populations.

In the first iteration, the Hirshfeld-I method uses neutral atomic densities (as in the original Hirshfeld scheme), ${\rho }_i^0(𝐫)$ with electronic population $N_i^0=\int 𝑑𝐫{\rho }_i^0(𝐫)=Z_i$. This affords charges

$$q_i^1=Z_i-\int 𝑑𝐫\left(\frac{{\rho }_i^0(𝐫)}{{\sum }_i^A{\rho }_i^0(𝐫)}\right)\rho (𝐫)=Z_i-N_i^1$$

(10.6)

on the first iteration. The new electronic population (number of electrons) for atom $i$ is $N_i^1$, and is derived from the promolecular populations $N_i^0$. One then computes new isolated atomic densities with $N_i^1=\int 𝑑𝐫{\rho }_i^1({𝐫}_1)$ and uses them to construct the promolecular densities in the next iteration. In general, the new weighting function for atom $i$ in the $k$th iteration is

$$w_{i,\mathrm{HI}}^k(𝐫)=\frac{{\rho }_i^{k-1}(𝐫)}{\sum _{i\in A}{\rho }_i^{k-1}(𝐫)}.$$

(10.7)

The atomic densities ${\rho }_i^k(𝐫)$ with corresponding fractional electron numbers $N_i^k$ are obtained by linear interpolation between ${\rho }_i^{0,\lfloor N_i^k\rfloor }(𝐫)$ and ${\rho }_i^{0,\lceil N_i^k\rceil }(𝐫)$ of the same atom:^{110, 241}

$${\rho }_i^k(𝐫)=\left(\lceil N_i^k\rceil -N_i^k\right){\rho }_i^{0,\lfloor N_i^k\rfloor }(𝐫)+\left(N_i^k-\lfloor N_i^k\rfloor \right){\rho }_i^{0,\lceil N_i^k\rceil }(𝐫),$$

(10.8)

where $\lfloor N_i^k\rfloor$ and $\lceil N_i^k\rceil$ denote the integers that bracket $N_i^k.$ The two atomic densities on the right side of Eq. (10.8) are obtained from densities ${\rho }_i^{0,Z_A-2},{\rho }_i^{0,Z_A-1},\mathrm{\dots },{\rho }_i^{0,Z_A+2}$ that are computed in advance. (That is, the method uses the neutral atomic density along with the densities for the singly- and doubly-charged cations and anions of the element in equation.) The Hirshfeld-I iterations are converged once the atomic populations change insignificantly between iterations, say .^{110, 916}

The iterative Hirshfeld scheme generally affords more reasonable charges as compared to the original Hirshfeld scheme. In LiF, for example, the original Hirshfeld scheme predicts atomic charges of $\pm$0.57 while the iterative scheme increases these charges to $\pm$0.93. The integral in Eq. (10.6) is evaluated by numerical quadrature, and the cost of each iteration of Hirshfeld-I is equal to the cost of computing the original Hirshfeld charges. The $rem variable SYM_IGNORE must be set to TRUE for Hirshfeld-I analysis.

$molecule
 -1  1
  O     1.197566    -0.108087     0.000000
  H     1.415397     0.827014     0.000000
  H     0.134830    -0.084378     0.000000
  F    -1.236389     0.012239     0.000000
$end

$rem
   SYM_IGNORE  true
   METHOD      B3LYP
   BASIS       6-31G*
   HIRSHITER   true
$end

In density functional theory calculations, the integration over the total density is evaluated on a molecular grid that is systematically broken up into interlocking multi-center atomic quadrature grids.⁵² This atomic quadrature scheme is predicated on the definition of atomic cell functions $P_a(\text{𝐫})$, that define smoothed Voronoi polyhedra centered about each atom. These cell functions are products of switching functions that define the atomic cell of atom $a$, and fall rapidly from $\approx 1$ near the nucleus of $a$, to $\approx 0$ near any other nucleus. The integration weights provided by this scheme are multiplied into the Lebedev quadrature weights in any practical DFT calculation:

$$w_n(𝐫)=\frac{P_n(𝐫)}{\sum _m{P}_m(𝐫)}$$

(10.9)

In some cases, it may be useful to print out the atomic Becke populations that are defined by these atomic cell functions. Becke population analysis may be requested by setting POP_BECKE to TRUE in the input file.

The default quadrature scheme uses atomic cell functions that intersect precisely at bond midpoints. Consequently, the default atomic cell functions will yield physically meaningless atomic populations. However, it is possible to shift the intersect of the atomic cell functions using an atomic radius criterion.⁵² In shifting the intersect of neighboring atomic cell functions, the point at which the Becke weights begin to fall from $\approx 1$ to $\approx 0$ changes depending on the atomic radius of each atom. While the choice of atomic radius is arbitrary, these atomic cell shifts introduce a physical basis for the partitioning of the underlying atomic quadrature. Two choices for atomic radii exist in Q-Chem for use with Becke weights, namely the empirically derived radii introduced by Bragg and Slater⁸⁸⁶ and the ab initio based weights of Pacios.⁷²²