The one-electron charge density, usually written as

$$\rho (\mathbf{r})=\sum _{\mu \nu}{P}_{\mu \nu}{\varphi}_{\mu}(\mathbf{r}){\varphi}_{\nu}(\mathbf{r})$$ | (11.1) |

represents the probability of finding an electron at the point $\mathbf{r}$, 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}(\mathbf{r})\right|}^{2}$$ | (11.2) |

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

$$N=\sum _{\mu \nu}{P}_{\mu \nu}{S}_{\mu \nu}=\sum _{\mu}{\left(\mathrm{\mathbf{P}\mathbf{S}}\right)}_{\mu \mu}=\mathrm{Tr}(\mathrm{\mathbf{P}\mathbf{S}})$$ | (11.3) |

where we interpret ${(\mathrm{\mathbf{P}\mathbf{S}})}_{\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{\mathbf{P}\mathbf{S}}\right)}_{\mu \mu}$$ | (11.4) |

where ${Z}_{A}$ is the atom’s nuclear charge. This is called a Mulliken population
analysis.^{896} Q-Chem performs a Mulliken population analysis by default.

POP_MULLIKEN

Controls running of Mulliken population analysis.

TYPE:

LOGICAL/INTEGER

DEFAULT:

TRUE (or 1)

OPTIONS:

FALSE (or 0)
Do not calculate Mulliken populations.
TRUE (or 1)
Calculate Mulliken populations.
2
Also calculate shell populations for each occupied orbital.
$-1$
Calculate Mulliken charges for both the ground state and any CIS,
RPA, or TDDFT excited states.

RECOMMENDATION:

Leave as TRUE, unless excited-state charges are desired.
Mulliken analysis is a trivial additional calculation, for ground or
excited states.

LOWDIN_POPULATION

Run Löwdin population analysis.

TYPE:

LOGICAL

DEFAULT:

FALSE

OPTIONS:

FALSE
Do not calculate Löwdin populations.
TRUE
Run Löwdin population analysis.

RECOMMENDATION:

None

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,^{591} 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”).^{118} 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,^{118} 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.^{369} (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.^{387})
For that reason, an alternative procedure based on atom-centered Lebedev grids
is also available,^{387} 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.

CHELPG

Controls the calculation of CHELPG charges.

TYPE:

LOGICAL

DEFAULT:

FALSE

OPTIONS:

FALSE
Do not calculate ChElPG charges.
TRUE
Compute ChElPG charges.

RECOMMENDATION:

Set to TRUE if desired. For large molecules, there is some overhead
associated with computing
ChElPG charges, especially if the number of grid points is large.

CHELPG_HEAD

Sets the “head space”^{118} (radial extent) of the ChElPG grid.

TYPE:

INTEGER

DEFAULT:

30

OPTIONS:

$N$
Corresponding to a head space of $N/10$, in Å.

RECOMMENDATION:

Use the default, which is the value recommended by Breneman and Wiberg.^{118}

CHELPG_DX

Sets the rectangular grid spacing for the traditional Cartesian ChElPG grid or
the spacing between concentric Lebedev shells (when the variables CHELPG_HA
and CHELPG_H are specified as well).

TYPE:

INTEGER

DEFAULT:

6

OPTIONS:

$N$
Corresponding to a grid space of $N/20$, in Å.

RECOMMENDATION:

Use the default, which corresponds to the “dense grid” of Breneman and
Wiberg,^{118}, unless the cost is prohibitive, in which case a
larger value can be selected. Note that this default value is set with the
Cartesian grid in mind and not the Lebedev grid. In the Lebedev case, a larger
value can typically be used.

CHELPG_H

Sets the Lebedev grid to use for hydrogen atoms.

TYPE:

INTEGER

DEFAULT:

NONE

OPTIONS:

$N$
Corresponding to a number of points in a Lebedev grid.

RECOMMENDATION:

CHELPG_H must always be less than or equal to CHELPG_HA. If it is greater,
it will automatically be set to the value of CHELPG_HA.

CHELPG_HA

Sets the Lebedev grid to use for non-hydrogen atoms.

TYPE:

INTEGER

DEFAULT:

NONE

OPTIONS:

$N$
Corresponding to a number of points in a Lebedev grid (see Section 5.5.1.

RECOMMENDATION:

None.

Hirshfeld population analysis^{381} 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 \mathit{d}\mathbf{r}\frac{{\rho}_{A}^{0}(\mathbf{r})}{{\sum}_{B}{\rho}_{B}^{0}(\mathbf{r})}\rho (\mathbf{r}),$$ | (11.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. (11.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),^{612} 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.*,^{123, 939} 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}(\mathbf{r})$ with electronic population ${N}_{i}^{0}=\int \mathit{d}\mathbf{r}{\rho}_{i}^{0}(\mathbf{r})={Z}_{i}$. This affords charges

$${q}_{i}^{1}={Z}_{i}-\int \mathit{d}\mathbf{r}\left(\frac{{\rho}_{i}^{0}(\mathbf{r})}{{\sum}_{i}^{A}{\rho}_{i}^{0}(\mathbf{r})}\right)\rho (\mathbf{r})={Z}_{i}-{N}_{i}^{1}$$ | (11.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 \mathit{d}\mathbf{r}{\rho}_{i}^{1}({\mathbf{r}}_{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}(\mathbf{r})=\frac{{\rho}_{i}^{k-1}(\mathbf{r})}{\sum _{i\in A}{\rho}_{i}^{k-1}(\mathbf{r})}.$$ | (11.7) |

The atomic densities ${\rho}_{i}^{k}(\mathbf{r})$ with corresponding fractional
electron numbers ${N}_{i}^{k}$ are obtained by linear interpolation between
${\rho}_{i}^{0,\lfloor {N}_{i}^{k}\rfloor}(\mathbf{r})$ and
${\rho}_{i}^{0,\lceil {N}_{i}^{k}\rceil}(\mathbf{r})$ of the same atom:^{123, 246}

$${\rho}_{i}^{k}(\mathbf{r})=\left(\lceil {N}_{i}^{k}\rceil -{N}_{i}^{k}\right){\rho}_{i}^{0,\lfloor {N}_{i}^{k}\rfloor}(\mathbf{r})+\left({N}_{i}^{k}-\lfloor {N}_{i}^{k}\rfloor \right){\rho}_{i}^{0,\lceil {N}_{i}^{k}\rceil}(\mathbf{r}),$$ | (11.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. (11.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
$$.^{123, 872}

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. (11.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.

HIRSHFELD

Controls running of Hirshfeld population analysis.

TYPE:

LOGICAL

DEFAULT:

FALSE

OPTIONS:

TRUE
Calculate Hirshfeld populations.
FALSE
Do not calculate Hirshfeld populations.

RECOMMENDATION:

None

HIRSHFELD_READ

Switch to force reading in of isolated atomic densities.

TYPE:

LOGICAL

DEFAULT:

FALSE

OPTIONS:

TRUE
Read in isolated atomic densities from previous Hirshfeld calculation from disk.
FALSE
Generate new isolated atomic densities.

RECOMMENDATION:

Use the default unless system is large. Note, atoms should be in the same
order with same basis set used as in the previous Hirshfeld calculation
(although coordinates can change). The previous calculation should be run with
the -save switch.

HIRSHFELD_SPHAVG

Controls whether atomic densities should be spherically averaged in pro-molecule.

TYPE:

LOGICAL

DEFAULT:

TRUE

OPTIONS:

TRUE
Spherically average atomic densities.
FALSE
Do not spherically average.

RECOMMENDATION:

Use the default.

CM5

Controls running of CM5 population analysis.

TYPE:

LOGICAL

DEFAULT:

FALSE

OPTIONS:

TRUE
Calculate CM5 populations.
FALSE
Do not calculate CM5 populations.

RECOMMENDATION:

None

HIRSHITER

Controls running of iterative Hirshfeld population analysis.

TYPE:

LOGICAL

DEFAULT:

FALSE

OPTIONS:

TRUE
Calculate iterative Hirshfeld populations.
FALSE
Do not calculate iterative Hirshfeld populations.

RECOMMENDATION:

None

HIRSHITER_THRESH

Controls the convergence criterion of iterative Hirshfeld population analysis.

TYPE:

INTEGER

DEFAULT:

5

OPTIONS:

$N$
Corresponding to the convergence criterion of $N/10000$, in $e$.

RECOMMENDATION:

Use the default, which is the value recommended in Ref. 123

$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.^{66} This atomic quadrature
scheme is predicated on the definition of atomic cell functions ${P}_{a}(\text{\mathbf{r}})$,
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}(\mathbf{r})=\frac{{P}_{n}(\mathbf{r})}{\sum _{m}{P}_{m}(\mathbf{r})}$$ | (11.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.

POP_BECKE

Controls the printing of atomic Becke populations.

TYPE:

LOGICAL

DEFAULT:

FALSE

OPTIONS:

TRUE
Print atomic Becke populations.
FALSE
Do not print atomic Becke populations.

RECOMMENDATION:

None

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.^{66} 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^{843} and the ab initio based weights of
Pacios.^{693}

BECKE_SHIFT

Controls atomic cell shifting in determination of Becke weights.

TYPE:

STRING

DEFAULT:

UNSHIFTED

OPTIONS:

UNSHIFTED
Use the original weighting scheme of Becke (bisection point).
BRAGG_SLATER
Use the empirically derived Bragg-Slater radii.
UNIVERSAL_DENSITY
Use the ab initio derived Pacios radii.

RECOMMENDATION:

If interested in the partitioning of the default atomic quadrature, use UNSHIFTED.
If using for physical interpretation, choose BRAGG_SLATER or UNIVERSAL_DENSITY.