30
$\begingroup$

I am evaluating the atomic charge of a system using inter-atomic potentials and comparing with using DFT. I know about the following types of partial charge: Mulliken, Bader, Qeq. I wonder what are the differences between the methods, such as pros and cons.

Perhaps we can get an answer explaining each of the following types of partial charge:

Here is a good starting point for anyone that wants to explain one of the above types of partial charge in the format of my (Nike Dattani's) answer.

| cite | improve this question | |
$\endgroup$
  • $\begingroup$ Don't use Mulliken for charge analysis. It is basically a random number generator, since it lacks a basis set limit. By choosing a different basis set representation, you can basically freely move the electrons around; e.g. in a one-center expansion all the electrons are counted for the expansion center whereas all the other nuclei in the system become bare. The same problem also exists in the Löwdin method (which somehow many people think is better!) but in a much worse way. Bader charges are often exaggerated. I've never heard of Qeq charges. $\endgroup$ – Susi Lehtola Jul 1 at 14:59
  • 1
    $\begingroup$ Might be a basic question, but what is "an IP system"? The most common abbreviation I've see for IP is ionization potential, though that doesn't fit in this context. $\endgroup$ – Tyberius♦ Jul 1 at 17:08
  • 3
    $\begingroup$ I would guess IP is interatomic potential and Qeq to stand for charge equilibration that is used for reactive force fields such as reax ff or comb. $\endgroup$ – Fabian Jul 1 at 17:16
  • 1
    $\begingroup$ Missing from the list is also NBO charges. $\endgroup$ – Andrew Rosen Jul 1 at 22:42
  • 2
    $\begingroup$ Thanks @AndrewRosen, I've added it now. You can also feel free to edit the question yourself next time (I wouldn't mind, and if the original user does, they can always click "rollback "easily). I guess NBO charges are also missing from this list: en.wikipedia.org/wiki/Partial_charge along with analogs based on many other types of bonds or orbitals. $\endgroup$ – Nike Dattani Jul 1 at 23:12

10 Answers 10

Active Oldest Votes
19
$\begingroup$

Just to add to the discussion:

  1. Mulliken charges are flawed in many aspects, but we know how and why, and therefore we accept its use, since it is simple and easily computed. But very dependent on the size of the basis set.

  2. Mulliken charges do not reproduce the dipole (or higher) moment, but can be made to do so easily: Thole, van Duijnen, "A general population analysis preserving the dipole moment” Theoret. Chim. Acta 1983, 63, 209–221 www.dx.doi.org/10.1007/BF00569246

  3. Based on a multipole expansion (used in ADF for the Coulomb potential) we have extended this to quadrupoles etc. M. Swart, P.Th. van Duijnen and J.G. Snijders "A charge analysis derived > from an atomic multipole expansion" J. Comput. Chem. 2001, 22, 79-88 http://www.dx.doi.org/10.1002/1096-987X(20010115)22:1%3C79::AID-JCC8%3E3.0.CO;2-B

Note that the multipoles result directly from the charge density, no fitting needed to the electrostatic potential at some grid outside the molecule (as done by other electrostatic potential fitted charge analyses). For each atom, its multipoles are represented by redistributed fractional atomic charges (with a weight function based on distance to keep these as close to the original atom as possible), the sum of these fractional atomic charges then add up to e.g. MDC-m (when only monopoles are redistributed), MDC-d (both monopoles and dipoles redistributed), MDC-q (monopoles, dipoles, quadrupoles redistributed).

For the N:+:C60 quartet state as mentioned by Tom, the MDC-m works best.

MDC-m, charge N -0.017, spin-dens. charge N 2.874; charge C 0.0003, spin-dens. charge C 0.002 MDC-d, charge N 0.136, spin-dens. charge N 0.722; charge C -0.002, spin-dens. charge C 0.038 MDC-q, charge N 0.062, spin-dens. charge N 0.729; charge C ranging from +0.05 to -0.05, spin-dens. charge C 0.038 In this case, the only other places where to put fractional charges are at the cage, hence not a representative sample for the method an sich.

I must add here that the Mulliken analysis works excellently: charge N -0.062, spin-dens. charge N 2.970. Hirshfeld (0.138) and Voronoi (0.259) charges for N are larger (no spin-density equivalent available within ADF).

  1. Finally, see also the following papers:

M. Cho, N. Sylvetsky, S. Eshafi, G. Santra, I. Efremenko, J.M.L. Martin The Atomic Partial Charges Arboretum: Trying to See the Forest for the Trees ChemPhysChem 2020, 21, 688-696 www.dx.doi.org/10.1002/cphc.202000040

G. Aullón, S Alvarez Oxidation states, atomic charges and orbital populations in transition metal complexes. Theor. Chem. Acc. 2009, 123, 67-73 www.dx.doi.org/10.1007/s00214-009-0537-9

G. Knizia Intrinsic Atomic Orbitals: An Unbiased Bridge between Quantum Theory and Chemical Concepts J. Chem. Theory Comp. 2013, 9, 4834-4843 www.dx.doi.org/10.1021/ct400687b

| cite | improve this answer | |
$\endgroup$
16
$\begingroup$

Bader charge analysis

In Bader's theory of Atoms in Molecules, we partition a molecule into "atoms" which are separated from each other by surfaces of minimum charge density:

                                         

You can then calculate the partial charges of the "atoms" in the molecule, e.g. H$_2$O might yield:

\begin{array}{cc} \rm{Atom} & \rm{Charge}\\ \hline \rm{O} & -1.150\\ \rm{H} & +0.425\\ \rm{H} & +0.425 \end{array}

Meaning that each hydrogen has "given away" 0.575 of an electron.

Pros:

Cons:

  • I will need help here. I only found out what "charge analysis" is after doing research to try to answer this question. I've never used (or learned) Mulliken, Bader, or any other type of charge analysis before.
  • Susi mentioned that Bader charges can be counterintuitive, which based on my hour of research on the topic today, seems to be a fundament of Bader's AIM concept for interpreting chemical bonds in general.
| cite | improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ Bader has a basis set limit. $\endgroup$ – Susi Lehtola Jul 1 at 17:16
  • $\begingroup$ @SusiLehtola I had no idea what Mulliken or Bader analysis was, before today, so I've taken your advice and added the basis set limit as a "pro" to this method. Feel free to edit my answer to make it correct if needed! Thank you. $\endgroup$ – Nike Dattani Jul 1 at 18:19
  • 1
    $\begingroup$ Bader himself allegedly never liked it, people calling the derived charges "Bader charges". Also, Multiwfn should be among the software list. This old question of mine is quite related: chemistry.stackexchange.com/q/48831/4945 $\endgroup$ – Martin - マーチン Jul 1 at 21:56
  • $\begingroup$ Thanks @Martin-マーチン I've added Multiwfn now. Feel free to edit my answer since I knew nothing about AIM or Bader's work before seeing this question, other than the fact that he wrote a book called QTAIM and is credited as inventing QTAIM (whatever that may be; I guess it has something to do with partitioning molecules into sections based on the surfaces of minimum charge density). $\endgroup$ – Nike Dattani Jul 1 at 23:06
16
$\begingroup$

Electrostatic potential (ESP) derived charges

Note that ESP1 derived charges include ChelpG (CHarges from ELectrostatic Potentials using a Grid-based method), the Merz–Kollman (MK) 2, and the RESP (restrained electrostatic potential) [3] scheme. While there are differences between the approaches the general idea is similar between the different methods. The main difference between the approaches is how the "grids" are selected. The points are selected in a regularly spaced cubic grid for CHELPG while the MK and the Resp schemes use points located on nested Connolly surfaces.

Pros:

  • Basis set completeness: The charges calculated with the CHELPG method are more systematic and predictable than charge methods based on the wave function or the electron density topology [4]

  • These types of charges are used in molecular mechanics. AMBER developers use RESP/MK, GLYCAM developers use RESP/CHELPG and CHARMM developers ESP/CHELPG and RESP/MK ChelpG, MK.

  • For floppy/flexible molecules ESP charges can be fit to multiple conformers to provide a better overall fit.

  • The charges of atoms can be constrained to certain values like unit charge. This is essential for making the building blocks of complex systems like proteins

    • Amino acids are parameterized with "caps" e.g. ACE or NME caps which complete the valence, these caps are constrained to unit charge so when they are removed the remaining "lego" piece has unit charge and can be easily combined with other amino acids to create a protein.

  • Many different programs can calculate RESP and CHELPG.

Cons:

  • ESP charges are not good for atoms which are located far from the points at which the electrostatic potential is calculated (see ChelpG). The RESP is apparently slightly better because it restrains certain charges with a hyperbolic restraint function. However when dealing with more complex systems the RESP charges are always calculated, when possible, with caps to limit the size and therefore to obtain more realistic charges [5].

References:

  1. The electrostatic potential is the amount of work necessary to move a charge from infinity to that point.
  2. Besler, B. H.; Merz, K. M.; Kollman, P. A. J Comput Chem 1990, 11, 431
  3. C. I. Bayly et. al. J. Phys. Chem. 97, 10269 (1993)
  4. Charge distribution in the water molecule -A comparison of methods, F. MARTIN, H. ZIPSE, Journal of Computational Chemistry, 2005, 10.1002/jcc.20157
  5. http://ambermd.org/tutorials/advanced/tutorial1/section1.htm
| cite | improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ +1. Two questions: How large is large (in the cons section), and, is there a reference for the statement that CHELPG is more systematic and predictable than density-based methods like Bader (written in the pros section)? $\endgroup$ – Nike Dattani Jul 2 at 3:08
  • $\begingroup$ Okay let me make that edit $\endgroup$ – Cody Aldaz Jul 2 at 3:15
16
$\begingroup$

Hirshfeld / CM5 (Charge Model 5)

My main interest is in explaining CM5 charges, but to do that it's necessary to briefly explain what Hirshfeld charges are.

Hirshfeld charges are obtained as: $$q_X=Z_X-\int\frac{\rho^0_X(\mathbf{r})}{\sum_Y\rho^0_Y(\mathbf{r})}\rho(r)d\mathbf{r}$$ where $Z_X$ is the atomic number of element $\ce{X}$, $\rho$ is the molecular density and $\rho_X^0$ is density of $\ce{X}$ as an isolated atom. Essentially, the density, and thus the charge, is partitioned in proportion to the atomic density. This approach has been found to be less basis set dependent than similar population analysis methods (e.g. Mulliken, Lowdin).

One drawback is that Hirshfeld charges alone don't do a great job of reproducing experimental observables like the molecular dipole moment, suggesting they may not be physically reasonable. This is where CM5 (Charge Model 5) charges come into play [1]. These are obtained as: $$q_k^\text{CM5}=q_k^{\text{Hirsh}}+\sum_{k'\neq k}T_{kk'}B_{kk'}$$ $$B_{kk'}=\exp\big[-\alpha(r_{kk'}-R_{Z_k}-R_{Z_{k'}})\big]$$ $$T_{k,k'}=\begin{cases}D_{Z_k,Z_{k'}} & Z_k,Z_{k'}=1,6,7,8 (\ce{H,C,N,O})\\ D_{Z_k}-D_{Z_{k'}} & \text{other elements}\end{cases}$$

The basis idea is that the Hirshfeld charge are corrected based on the bond order between two atoms $B_{kk'}$. The bond order itself is parameterized by $\alpha$ and its effect on the charge of a given element is parameterized by $T_{kk'}$, more specifically the underlying $D_k$ or $D_{k,k'}$ ($k$ and $k'$ index all the atoms in the molecule, but $T_{kk'}$ only depends on which two elements are involved). This approach gives physically reasonable atomic charges that can be used to accurately calculate molecular dipole moments. Since Hirshfeld charges are fairly simple to compute for any SCF calculation, CM5 charges can easily be added on top.

The main drawback is that CM5 was only parameterized for some elements ($\ce{H}$-$\ce{Ca}$, $\ce{Zn}$-$\ce{Br}$, $\ce{I}$, along with some special parameters for common organic pairs like $\ce{O-H}$, $\ce{C-H}$, etc.). However there are still values for all elements. $D_{Z_K}$ is set equal to 0 for all transitions metals, Lanthanides, and Actinides and the elements in the same column of the periodic table are constrained to satisfy $D_{Z_k}=CD_{\bar{Z_k}}$, where $C$ is another parameter $\bar{Z_k}$ refers to the next element up in the column.

CM$x$ charges with $x<5$ employ a similar approach, but were parameterized on small training sets to correct Lowdin charges. Due to the basis set dependence of Lowdin charges, they are only suitable for use with certain small basis sets, the largest being 6-31+G(d,p).

  1. Aleksandr V. Marenich, Steven V. Jerome, Christopher J. Cramer, and Donald G. Truhlar Journal of Chemical Theory and Computation 2012 8 (2), 527-541 DOI: 10.1021/ct200866d
| cite | improve this answer | |
$\endgroup$
15
$\begingroup$

To give a example of the "random number generator" behavior of Mulliken charge analysis pointed by @SusiLehtola, below I'm showing the results of some test runs I did before on the CO2 molecule using Psi4 version 1:1.1-5 (default version on Ubuntu 18.04 repository as of july 2020).

For reproducibility purposes, first my input files:

user@machine:~/Documentos/stackexchange$ more *.in | cat
::::::::::::::
CO2_dipole_631plusGd.in
::::::::::::::
memory 4 Gb

set basis  6-31+G(d)
molecule {
0 1
   C       -3.47367        0.73246        0.22361
   O       -2.43476        1.12414       -0.22175
   O       -4.51237        0.34053        0.66926
}
optimize('B3LYP-D')
E, wfn = energy('B3LYP-D', return_wfn=True)
oeprop(wfn, "MULLIKEN_CHARGES", "DIPOLE", title = "CO2 B3LYP-D")
::::::::::::::
CO2_dipole_631plusplusGdp.in
::::::::::::::
memory 4 Gb

set basis  6-31++G(d_p)
molecule {
0 1
   C       -3.47367        0.73246        0.22361
   O       -2.43476        1.12414       -0.22175
   O       -4.51237        0.34053        0.66926
}
optimize('B3LYP-D')
E, wfn = energy('B3LYP-D', return_wfn=True)
oeprop(wfn, "MULLIKEN_CHARGES", "DIPOLE", title = "CO2 B3LYP-D")
::::::::::::::
CO2_dipole_augccpVDZ.in
::::::::::::::
memory 4 Gb

set basis  aug-cc-pVDZ
molecule {
0 1
   C       -3.47367        0.73246        0.22361
   O       -2.43476        1.12414       -0.22175
   O       -4.51237        0.34053        0.66926
}
optimize('B3LYP-D')
E, wfn = energy('B3LYP-D', return_wfn=True)
oeprop(wfn, "MULLIKEN_CHARGES", "DIPOLE", title = "CO2 B3LYP-D")
::::::::::::::
CO2_dipole_augccpVTZ.in
::::::::::::::
memory 4 Gb

set basis  aug-cc-pVTZ
molecule {
0 1
   C       -3.47367        0.73246        0.22361
   O       -2.43476        1.12414       -0.22175
   O       -4.51237        0.34053        0.66926
}
optimize('B3LYP-D')
E, wfn = energy('B3LYP-D', return_wfn=True)
oeprop(wfn, "MULLIKEN_CHARGES", "DIPOLE", title = "CO2 B3LYP-D")
::::::::::::::
CO2_dipole_augpcseg1.in
::::::::::::::
memory 4 Gb

set basis  aug-pcseg-1
molecule {
0 1
   C       -3.47367        0.73246        0.22361
   O       -2.43476        1.12414       -0.22175
   O       -4.51237        0.34053        0.66926
}
optimize('B3LYP-D')
E, wfn = energy('B3LYP-D', return_wfn=True)
oeprop(wfn, "MULLIKEN_CHARGES", "DIPOLE", title = "CO2 B3LYP-D")
::::::::::::::
CO2_dipole_augpcseg2.in
::::::::::::::
memory 4 Gb

set basis  aug-pcseg-2
molecule {
0 1
   C       -3.47367        0.73246        0.22361
   O       -2.43476        1.12414       -0.22175
   O       -4.51237        0.34053        0.66926
}
optimize('B3LYP-D')
E, wfn = energy('B3LYP-D', return_wfn=True)
oeprop(wfn, "MULLIKEN_CHARGES", "DIPOLE", title = "CO2 B3LYP-D")
::::::::::::::
CO2_dipole_pcseg1.in
::::::::::::::
memory 4 Gb

set basis  pcseg-1
molecule {
0 1
   C       -3.47367        0.73246        0.22361
   O       -2.43476        1.12414       -0.22175
   O       -4.51237        0.34053        0.66926
}
optimize('B3LYP-D')
E, wfn = energy('B3LYP-D', return_wfn=True)
oeprop(wfn, "MULLIKEN_CHARGES", "DIPOLE", title = "CO2 B3LYP-D")
::::::::::::::
CO2_dipole_pcseg2.in
::::::::::::::
memory 4 Gb

set basis  pcseg-2
molecule {
0 1
   C       -3.47367        0.73246        0.22361
   O       -2.43476        1.12414       -0.22175
   O       -4.51237        0.34053        0.66926
}
optimize('B3LYP-D')
E, wfn = energy('B3LYP-D', return_wfn=True)
oeprop(wfn, "MULLIKEN_CHARGES", "DIPOLE", title = "CO2 B3LYP-D")

Now the results I got:

user@machine:~/Documentos/stackexchange$ grep -A 4 'Mulliken Charges: (a.u.)' *.out
CO2_dipole_631plusGd.out:  Mulliken Charges: (a.u.)
CO2_dipole_631plusGd.out-   Center  Symbol    Alpha    Beta     Spin     Total
CO2_dipole_631plusGd.out-       1     C     2.62564  2.62564  0.00000  0.74871
CO2_dipole_631plusGd.out-       2     O     4.18718  4.18718  0.00000 -0.37436
CO2_dipole_631plusGd.out-       3     O     4.18718  4.18718  0.00000 -0.37436
--
CO2_dipole_631plusplusGdp.out:  Mulliken Charges: (a.u.)
CO2_dipole_631plusplusGdp.out-   Center  Symbol    Alpha    Beta     Spin     Total
CO2_dipole_631plusplusGdp.out-       1     C     2.62564  2.62564  0.00000  0.74871
CO2_dipole_631plusplusGdp.out-       2     O     4.18718  4.18718  0.00000 -0.37436
CO2_dipole_631plusplusGdp.out-       3     O     4.18718  4.18718  0.00000 -0.37436
--
CO2_dipole_augccpVDZ.out:  Mulliken Charges: (a.u.)
CO2_dipole_augccpVDZ.out-   Center  Symbol    Alpha    Beta     Spin     Total
CO2_dipole_augccpVDZ.out-       1     C     2.82315  2.82315  0.00000  0.35370
CO2_dipole_augccpVDZ.out-       2     O     4.08842  4.08842  0.00000 -0.17685
CO2_dipole_augccpVDZ.out-       3     O     4.08843  4.08843  0.00000 -0.17686
--
CO2_dipole_augccpVTZ.out:  Mulliken Charges: (a.u.)
CO2_dipole_augccpVTZ.out-   Center  Symbol    Alpha    Beta     Spin     Total
CO2_dipole_augccpVTZ.out-       1     C     2.80993  2.80993  0.00000  0.38014
CO2_dipole_augccpVTZ.out-       2     O     4.09503  4.09503  0.00000 -0.19007
CO2_dipole_augccpVTZ.out-       3     O     4.09504  4.09504  0.00000 -0.19007
--
CO2_dipole_augpcseg1.out:  Mulliken Charges: (a.u.)
CO2_dipole_augpcseg1.out-   Center  Symbol    Alpha    Beta     Spin     Total
CO2_dipole_augpcseg1.out-       1     C     2.35311  2.35311  0.00000  1.29377
CO2_dipole_augpcseg1.out-       2     O     4.32345  4.32345  0.00000 -0.64689
CO2_dipole_augpcseg1.out-       3     O     4.32344  4.32344  0.00000 -0.64688
--
CO2_dipole_augpcseg2.out:  Mulliken Charges: (a.u.)
CO2_dipole_augpcseg2.out-   Center  Symbol    Alpha    Beta     Spin     Total
CO2_dipole_augpcseg2.out-       1     C     2.51884  2.51884  0.00000  0.96233
CO2_dipole_augpcseg2.out-       2     O     4.24057  4.24057  0.00000 -0.48114
CO2_dipole_augpcseg2.out-       3     O     4.24059  4.24059  0.00000 -0.48119
--
CO2_dipole_pcseg1.out:  Mulliken Charges: (a.u.)
CO2_dipole_pcseg1.out-   Center  Symbol    Alpha    Beta     Spin     Total
CO2_dipole_pcseg1.out-       1     C     2.71634  2.71634  0.00000  0.56732
CO2_dipole_pcseg1.out-       2     O     4.14183  4.14183  0.00000 -0.28366
CO2_dipole_pcseg1.out-       3     O     4.14183  4.14183  0.00000 -0.28366
--
CO2_dipole_pcseg2.out:  Mulliken Charges: (a.u.)
CO2_dipole_pcseg2.out-   Center  Symbol    Alpha    Beta     Spin     Total
CO2_dipole_pcseg2.out-       1     C     2.70233  2.70233  0.00000  0.59534
CO2_dipole_pcseg2.out-       2     O     4.14883  4.14883  0.00000 -0.29767
CO2_dipole_pcseg2.out-       3     O     4.14883  4.14883  0.00000 -0.29767

As you can see, the results on Mulliken charge analysis can change a lot, depending on the basis set you use to run the calculation. About Bader and Qeq charges, I haven't much to say, as I lack experience dealing with them.

| cite | improve this answer | |
$\endgroup$
14
$\begingroup$

Mulliken population analysis

The Mulliken charge scheme is based on the Linear Combination of Atomic Orbitals (LCAO), so, it is based on the system wave function and was described in a series of papers by R. S. Mulliken1,2,3,4.

The idea is that the normalized Molecular Orbital (MO), $\phi_i$, of a diatomic molecule is written as a linear combination of normalized Atomic Orbitals (AO), $\chi_j$ and $\chi_k$:

$$\phi_i = c_{ij} \chi_j + c_{ik} \chi_k$$

Assuming that the MO is occupied by $N$ electrons, these $N$ electrons can be distributes as:

$$N {\phi_i}^2 = N {c_{ij}}^2 {\chi_j}^2 + N {c_{ik}}^2 {\chi_k}^2 + 2 N c_{ik} \chi_i \chi_j$$

Integrating over all electronic coordinates and as the MO and AO are normalized::

$$N = N {c_{ij}}^2 + N {c_{ik}}^2 + 2 N c_{ij} c_{ik} S_{jk}$$

$$1 = {c_{ij}}^2 + {c_{ik}}^2 + 2 c_{ij} c_{ik} S_{jk}$$

where $S_{jk}$ is the overlap integral of the two Atomic Orbitals.

According to Mulliken interpretation, the subpopulations $N {c_{ij}}^2$ and $N {c_{ik}}^2$ are called the net atomic populations on atoms $j$ and $k$ and $2 N c_{ij} c_{ik} S_{jk}$ is called the overlap population.

A convenient way to re-write the previous equation is in matrix form:

$${P_i} = \left( {\begin{array}{*{20}{c}} {c_{ij}^2}&{2{c_{ij}}c{}_{ik}{S_{jk}}}\\ {2{c_{ij}}c{}_{ik}{S_{jk}}}&{c_{ik}^2} \end{array}} \right)$$

To take into account the populations from all the electrons in all the molecular orbitals, the net population matrix can be defined as

$${\rm{Net Population}} = \sum\limits_{i = occupied} {{P_i}}. $$

As pros, we have that these populations are easily calculated (almost any software can calculate them). As cons, they are heavily dependent on system wave functions and then, on the chosen basis sets (not random!).

References:

[1] Mulliken, R.S. Electronic Population Analysis on LCAO-MO. Molecular Wave Functions. I, J. Chem. Phys. (1955), 23, 1833-1840.

[2] Mulliken, R.S. Electronic Population Analysis on LCAO-MO. Molecular Wave Functions. II. Overlap Populations, Bond Orders, and Covalent Bond Energies, J. Chem. Phys. (1955), 23, 1841-1846.

[3] Mulliken, R.S. Electronic Population Analysis on LCAO-MO. Molecular Wave Functions. III. Effects of Hybridization on Overlap and Gross AO Populations, J. Chem. Phys. (1955), 23, 2338-2342.

[4] Mulliken, R.S. Electronic Population Analysis on LCAO-MO. Molecular Wave Functions. IV. Bonding and Antibonding in LCAO and Valence-Bond Theories, J. Chem. Phys. (1955), 23, 2343-2346.

| cite | improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ +1. I suppose you could say something about the generalization beyond diatomics? the way your answer is written so far, could accidentally give the impression that it's only about diatomics. $\endgroup$ – Nike Dattani Jul 1 at 23:15
13
$\begingroup$

Mulliken and Löwdin population analysis

In the atomic orbital basis set (enumerated in Greek indices), one finds that the number of electrons $N$ is equal to the trace of the product $\mathbf{PS}$ $$ N = \sum_\mu \left(\mathbf{PS}\right)_{\mu\mu} = \mathrm{Tr} \ \mathbf{PS} $$ where $\mathbf{P}$ is the density matrix, $\mathbf{S}$ is the AO overlap matrix, and the sum runs over all basis functions. One can then decide to partition the electron population by associating non-intersecting subsets of the basis set to atoms, typically by taking those centered on atom $A$ as belonging to $A$. We will denote this as $\mu \in A$ and define the Mulliken charge on $A$ as $$ q_A^\text{Mulliken} = Z_A \ - \sum_{\mu \ \in \ A} \left(\mathbf{PS}\right)_{\mu\mu} $$ where $Z_A$ is the nuclear charge of $A$.

Non-uniqueness of $\mathbf{PS}$:

The trace has the property of cylic permutability: $$ \mathrm{Tr} \ \mathbf{ABC} = \mathrm{Tr} \ \mathbf{CAB} = \mathrm{Tr} \ \mathbf{BCA} $$ which can be applied to $\mathbf{PS}$ as follows: $$ \mathrm{Tr} \ \mathbf{PS} = \mathrm{Tr} \ \mathbf{P}\mathbf{S}^{1-x}\mathbf{S}^{x} = \mathrm{Tr} \ \mathbf{S}^{x}\mathbf{P}\mathbf{S}^{1-x} $$ which holds at least for $x \in \mathbb{Q}$.

One can then set $x = \frac{1}{2}$ and obtain the Löwdin charge $$ q_A^\text{Löwdin} = Z_A - \sum_{\mu \ \in \ A} \left(\mathbf{S}^\frac{1}{2}\mathbf{P}\mathbf{S}^{\frac{1}{2}}\right)_{\mu\mu} $$

Discussion/Cons: As pointed out by others, theses types of analysis are particularly susceptible to the basis set. Intramolecular basis set superposition error (BSSE) is also an important factor. I would further note that any scheme for atomic charge is flawed in the sense that it does not represent an observable. Therefore, values should only be considered relative to prototypical systems.

Source: A Szabo, NS Ostlund: Modern Quantum Chemistry, Dover Publications, 1996.

| cite | improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ +1. I just made some small changes so that people don't have to scroll through as much whitespace. $\endgroup$ – Nike Dattani Jul 4 at 0:08
8
$\begingroup$

Don't use Mulliken for charge analysis. It is basically a random number generator, since it lacks a basis set limit. By choosing a different basis set representation, you can basically freely move the electrons around; e.g. in a one-center expansion all the electrons are counted for the expansion center whereas all the other nuclei in the system become bare. The same problem also exists in the Löwdin method (which somehow many people think is better!) but in a much worse way.

Bader charges are often exaggerated.

I've never heard of Qeq charges.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Generally, Qeq, charge equilibration method is used with other interatomic potentials (ReaxFF, COMB3). It iteratively calculate the charge of atoms in a system to describe local chemical environment. The QEq method explicitly includes shielding between overlapping charge clouds. (ref: A.K. Rappe, W.A. Goddard III, Journal of Physical Chemistry 95 (8) (1991) 3358.) $\endgroup$ – Dong gi Kang Jul 1 at 18:05
8
+100
$\begingroup$

Voronoi charges

Voronoi charges (here called VC) are based on the partition of the real space in a system into Voronoi polyhedra.$^1$ A given point in space belongs to the polyhedron of some atom if the point is closer to that atom than to any other atom. This allows to partition the space and thus assign the point charge expectation value of a given point to the appropriate atom. After summing over all points/integrating and adding the nuclear charge, one obtains the VC.

The algorithm requires a numerical integration and thus, use of a grid. While grids are common in DFT, it is not a priori clear what their performance is for VC determination. The grids might not give good accuracy in the cell border regions.

VC would seem to be sensitive to the geometry. Consider the case of the molecule $\ce{HF}$: As one (artificially) shortens the bond distance, hydrogen will leech electron density from fluorine, though most chemists would happily assign the same $\delta\pm$ to the atoms at any (bonded) distance.

$^1$ There are several different names for this concept due to their (re)discovery by many different people. See the Wikipedia article on the Voronoi diagram for some of them.

Source: F Jensen: Introduction to Computational Chemistry, 2nd ed., Wiley, 2007.

| cite | improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ Reason I awarded the bounty here: Cody's answer was actually the only one that followed my format as the bounty & my edit to the question requested. But I tend to prefer to award bounties to people that answer after I offer the bounty, so it either had to go to TAR86's Mulliken+Lowdin or Voronoi answer. This one followed the format more closely, as recommended here. Granted your summary of Lowdin needed Mulliken too, & I do like your answers. Overall, thank you for your contributions to the site! $\endgroup$ – Nike Dattani Jul 11 at 16:51
4
$\begingroup$

Density-Derived Electrostatic and Chemical approaches

Net atomic charges have two primary applications (dual use): (1) to quantify charge transfer between atoms in materials; this identifies cations and anions and (2) to provide an electrostatic model in classical force fields using atomistic simulations (e.g., classical molecular dynamics or Monte Carlo simulations).

The older charge partitioning methods were not optimized for this dual use. For example, charges fit specifically to the electrostatic potential (CHELP, CHELPG, Merz-Kollman, etc.) often did not give any kind of reasonable chemical description for buried atoms.

The density-derived electrostatic and chemical (DDEC) methods are optimized to assign net atomic charges that give a good approximation both to the electrostatic potential surrounding the material as well as to the chemical charge states of atoms in materials. In other words, they are optimized for dual use.

A key design consideration in the DDEC family of methods is to create methods that work across an extremely broad range of material types including molecules, ions, nano-structured materials, metals, insulators, dense and porous solids, organometallics, and polymers for all chemical elements of atomic number 1 to 109.

Another key design consideration is that iterative processes to compute DDEC net atomic charges, atomic spin moments, and other atom-in-material properties should be rapid, robust, and converge to unique solutions. Several generations of improvements of the DDEC methods have been published. Unfortunately, there were some problems with the earliest DDEC approaches (e.g., DDEC1, DDEC2, and DDEC3) that caused non-unique convergence ('runaway charges') in some materials. The latest generation (DDEC6) fixes these convergence problems and is described in the following publications:

T. A. Manz and N. Gabaldon Limas, “Introducing DDEC6 atomic population analysis: part 1. Charge partitioning theory and methodology,” RSC Advances, 6 (2016) 47771-47801 DOI:10.1039/c6ra04656h

N. Gabaldon Limas and T. A. Manz, “Introducing DDEC6 atomic population analysis: part 2. Computed results for a wide range of periodic and nonperiodic materials,” RSC Advances, 6 (2016) 45727-45747 DOI:10.1039/c6ra05507a

T. A. Manz, “Introducing DDEC6 atomic population analysis: part 3. Comprehensive method to compute bond orders,” RSC Advances, 7 (2017) 45552-45581 (open access) DOI:10.1039/c7ra07400j

N. Gabaldon Limas and T. A. Manz, “Introducing DDEC6 atomic population analysis: part 4. Efficient parallel computation of net atomic charges, atomic spin moments, bond orders, and more,” RSC Advances, 8 (2018) 2678-2707 (open access) DOI:10.1039/c7ra11829e

The "density-derived" refers to atom-in-material properties (e.g., net atomic charges, atomic spin moments, bond orders, atomic multipoles, etc.) that are computed as functionals of the electron and spin density distributions. One can also imagine additional properties that are computed from the first-order density matrix or molecular orbitals that are not functionals of the electron and spin density distributions. These "orbital-derived" properties include spdfg populations of atoms in materials, projected density of states plots, bond order components assigned to individual orbitals, etc. Together, these "density-derived" and "orbital-derived" properties form the Standard Atoms in Materials Method (SAMM). In other words, recent generation DDEC (e.g., DDEC6) method is the "density-derived" part of the SAMM.

A key design consideration of the SAMM approach is that all of the various component methods should work together to provide a chemically consistent description of atoms in materials. This means the net atomic charges, atomic spin moments, atomic multipoles, bond orders, bond order components, spdfg populations, polarizabilities, dispersion coefficients, and projected density of states plots should be chemically compatible with each other. For example, when summing the individual populations of the spdfg subshell populations for spin-up and spin-down electrons in a magnetic material, these yield the prior-computed net atomic charges and atomic spin moments and have chemical consistency with the computed bond orders. Moreover, the "orbital-derived" properties from the SAMM method are designed to be chemically consistent with "density-derived" properties. For example, integrating the projected density of states (PDOS) curves for a particular atom in a material regenerates the prior-computed net atomic charges, atomic spin moments, and bond orders.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ This is a great answer, can't beat having a method explained by its developer. Welcome to the site! $\endgroup$ – Tyberius♦ Sep 16 at 22:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.