2.6 Real space functions

The "Real space functions" in this manual are referred to as the functions whose variables are coordinate of the three-dimension space of present system. Real space function analysis is one of the most important functions of Multiwfn, the supported real space functions are listed below. All wavefunctions are assumed to be real type, all units are in atomic unit (a.u.).

Notice that for speeding up calculation, especially for big system, when evaluating an exponential function (except for some real space functions, such as 12, 14 and 16), if the exponent is more negative than -40, then this evaluation will be skipped. The default cutoff value is safe enough and cannot cause detectable loss of precision even in quantitative analysis, you can also disable this treatment or adjust cutoff, see “expcutoff” in settings.ini.

1 Electron density

\[\rho(\mathbf{r})=\sum_i \eta_i\left|\varphi_i(\mathbf{r})\right|^2=\sum_i \eta_i\left|\sum_\mu C_{\mu,i}\chi_\mu(\mathbf{r})\right|^2\]

where η_i is occupation number of orbital i, \(\varphi\) is orbital wavefunction, χ is basis function. C is coefficient matrix, the element of the ith row jth column corresponds to the expansion coefficient of orbital j respect to basis function i. Atomic unit for electron density can be explicitly written as 1/Boh\(r^{3}\) (which corresponds to 6.74833/ų, since 1 Bohr = 0.529177Å).

When input file does not contain GTF information, this function will be calculated as promolecular density, which is approximate molecular electron density simply constructed by superposing built-in spherically averaged free-state atomic density of all atoms in the system. See Appendix 3 on how the built-in atomic densities were derived.

It is worth to note that distribution character of valence electron density is much more informative than electron density, this point was thoroughly discussed in my work Acta Phys. -Chim. Sin., 34, 503 (2018) DOI: 10.3866/pku.Whxb201709252.

2 Gradient norm of electron density

\[\left|\nabla\rho(\mathbf{r})\right|=\sqrt{\left(\frac{\partial\rho(\mathbf{r})}{\partial x}\right)^2+\left(\frac{\partial\rho(\mathbf{r})}{\partial y}\right)^2+\left(\frac{\partial\rho(\mathbf{r})}{\partial z}\right)^2}\]

3 Laplacian of electron density

\[\nabla^2\rho(\mathbf{r})=\frac{\partial^2\rho(\mathbf{r})}{\partial x^2}+\frac{\partial^2\rho(\mathbf{r})}{\partial y^2}+\frac{\partial^2\rho(\mathbf{r})}{\partial z^2}\]

The positive and negative values of this function correspond to electron density is locally depleted and locally concentrated, respectively. The relationships between \({\mathrm{\nabla }}^{2}\rho\) and valence shell electron pair repulsion (VSEPR) model, chemical bond type, electron localization and chemical reactivity have been built by Bader and many other researchers.

Notice that in some other codes, such as AIMALL and AIM2000, Laplacian of electron density is defined as \((-\frac{1}{4}){\mathrm{\nabla }}^{2}\rho\), which differs from the value given by Multiwfn by a factor of -1/4.

If “laplfac” in settings.ini is set to other value rather than the default one 1.0, \({\mathrm{\nabla }}^{2}\rho\) will be multiplied with this value.

The \({\mathrm{\nabla }}^{2}\rho\) outputted by Multiwfn is in atomic unit, which can be explicitly written as 1 Boh\(r^{-5}\) (corresponding to 24.09874 Å^-5, since 1 Bohr = 0.529177 Å).

4 Value of orbital wavefunction; 44 Orbital probability density

They are respectively defined as

\[\begin{gathered}\varphi_i(\mathbf{r})=\sum_\mu C_{\mu,i}\chi_\mu(\mathbf{r})\\ \rho_i(\mathbf{r})=\left|\varphi_i(\mathbf{r})\right|^2\end{gathered}\]

When you select these functions from function list, you will be prompted to input the index of the orbital of interest.

5 Electron spin density

Spin density is defined as the difference between alpha and beta density

\[\rho^s(\mathbf{r})=\rho^\alpha(\mathbf{r})-\rho^\beta(\mathbf{r})\]

If “ipolarpara” in settings.ini is set to 1, then spin polarization parameter function will be returned instead of spin density

\[\zeta(\mathbf{r})=\frac{\rho^\alpha(\mathbf{r})-\rho^\beta(\mathbf{r})}{\rho^\alpha(\mathbf{r})+\rho^\beta(\mathbf{r})}\]

The absolute value of ζ going from zero to unity corresponds to the local region going from unpolarized case to completely polarized case.

6 Hamiltonian kinetic energy density K(r)

The kinetic energy density is not uniquely defined, since the expected value of kinetic energy operator \(\left<\varphi\left|-(1/2)\nabla^2\right|\varphi\right>\) can be recovered by integrating kinetic energy density from alternative definitions. One of commonly used definitions is

\[K(\mathbf{r})=-\frac{1}{2}\sum_i \eta_i \varphi_i^*(\mathbf{r})\nabla^2\varphi_i(\mathbf{r})\]

7 Lagrangian kinetic energy density G(r)

Relative to K(r), the local kinetic energy definition given below guaranteed positiveness everywhere, hence the physical meaning is clearer and more commonly used. G(r) is also known as positive definite kinetic energy density.

\[G(\mathbf{r})=\frac{1}{2}\sum_i\eta_i\left|\nabla\varphi_i(\mathbf{r})\right|^2=\frac{1}{2}\sum_i\eta_i\left\{\left(\frac{\partial\varphi_i(\mathbf{r})}{\partial x}\right)^2+\left(\frac{\partial\varphi_i(\mathbf{r})}{\partial y}\right)^2+\left(\frac{\partial\varphi_i(\mathbf{r})}{\partial z}\right)^2\right\}\]

Since \(\nabla\rho_i=\nabla(\eta_i\varphi_i^2)=2\eta_i\varphi_i\nabla\varphi_i\), there is an equivalent form of G(r): \(\frac{1}{8}\sum_i\frac{\nabla\rho_i\cdot\nabla\rho_i}{\rho_i}\).

K(r) and G(r) are directly related by Laplacian of electron density

\[\nabla^2\rho(\mathbf{r})/4=G(\mathbf{r})-K(\mathbf{r})\]

8 Electrostatic potential from nuclear / atomic charges

\[V_{\mathrm{nuc}}(\mathbf{r})=\sum_A\frac{Z_A}{\left|\mathbf{r}-\mathbf{R}_A\right|}\]

where \(R_{A}\) and \(Z_{A}\) denote position vector and nuclear charge of atom A, respectively. If pseudopotential is used, then Z is the number of explicitly expressed electrons. When .chg file is used as input, Z will stand for the atomic charges recorded in the file (the fourth column), at this time \(V_{\mathrm{nu}}\) is useful for analyzing the difference between exact electrostatic potential and the electrostatic potential reproduced by atomic charges.

Notice that at nuclear positions, this function will be infinite and may cause some numerical problems in program, hence at these cases this function always returns 1000 instead of infinity.

9 Electron localization function (ELF)

The larger the electron localization is in a region, the more likely the electron motion is confined within it. If electrons are completely localized, then they can be distinguished from the ones outside. Bader found that the regions which have large electron localization must have large magnitudes of Fermi hole integration. However, the Fermi hole is a six-dimension function and thus difficult to study visually. Becke and Edgecombe noted that spherically averaged like-spin conditional pair probability has direct correlation with the Fermi hole and then suggested electron localization function (ELF) in the paper J. Chem. Phys., 92, 5397 (1990). Angew. Chem. Int. Ed., 137, e202504895 (2025) DOI: 10.1002/anie.202504895 is a nice review illustrating the great practical value of ELF in deciphering electronic structure, and Chapter 5 of Theoretical Aspects of Chemical Reactivity (2007) provides some additional information.

The ELF adopted in Multiwfn is the form generalized for spin-polarized system by me, see my paper Acta Phys. -Chim. Sin., 27, 2786 (2011) for introduction.

\[\mathrm{ELF}(\mathbf{r})=\frac{1}{1+[D(\mathbf{r})/D_0(\mathbf{r})]^2}\]

where

\[\begin{gathered}D(\mathbf{r})=\frac{1}{2}\sum_i\eta_i\left|\nabla\varphi_i(\mathbf{r})\right|^2-\frac{1}{8}\left[\frac{\left|\nabla\rho_\alpha(\mathbf{r})\right|^2}{\rho_\alpha(\mathbf{r})}+\frac{\left|\nabla\rho_\beta(\mathbf{r})\right|^2}{\rho_\beta(\mathbf{r})}\right]\\ D_0(\mathbf{r})=\frac{3}{10}(6\pi^2)^{2/3}\left[\rho_\alpha(\mathbf{r})^{5/3}+\rho_\beta(\mathbf{r})^{5/3}\right]\end{gathered}\]

For closed-shell system, since \(\rho_\alpha=\rho_\beta=(1/2)\rho\), D and \(D_{0}\) terms can be simplified as

\[\begin{gathered}D(\mathbf{r})=\frac{1}{2}\sum_i\eta_i\left|\nabla\varphi_i(\mathbf{r})\right|^2-\frac{1}{8}\frac{\left|\nabla\rho(\mathbf{r})\right|^2}{\rho(\mathbf{r})}\\ D_0(\mathbf{r})=(3/10)(3\pi^2)^{2/3}\rho(\mathbf{r})^{5/3}\end{gathered}\]

Savin et al. have reinterpreted ELF in the view of kinetic energy, see Angew. Chem. Int. Ed. Engl., 31, 187 (1992), which makes ELF also meaningful for Kohn-Sham DFT wavefunction and even multi-configuration wavefunction. The first term of D(r) can be seen as the exact kinetic energy density of the non-interacting electron system defined by KS-DFT theory, namely \(\tau_{\mathrm{S}}(\mathbf{r})=(1/2)\sum_i\eta_i\left|\nabla\varphi_i(\mathbf{r})\right|^2\), while the second term is equivalent to Weizsäcker kinetic energy density \(\tau_W(\mathbf{r})=(1/8)\left|\nabla\rho\right|^2/\rho\), therefore the \(D(\mathbf{r})=\tau_{\mathrm{S}}(\mathbf{r})-\tau_{\mathrm{W}}(\mathbf{r})\) reveals the excess kinetic energy density caused by Pauli repulsion and it is known as Pauli kinetic energy density. The \(D_{0}\)(r) can be interpreted as Thomas-Fermi kinetic energy density τ_TF, which is the exact kinetic energy density of non-interacting, uniform electron gas. Since \(D_{0}\)(r) is introduced into ELF as reference, what the ELF reveals is actually degree of relative localization.

ELF is within the range of [0,1]. A large ELF value means that electrons are greatly localized, indicating that there is a covalent bond, a lone pair or inner shells of the atom involved. ELF has been widely used for a wide variety of systems, such as organic and inorganic small molecules, atomic crystals, coordination compounds, clusters, and for different problems, such as the revealing atomic shell structure, classification of chemical bonding, verification of charge-shift bonds, studying aromaticity.

Note

BTW: Someone questioned the implementation of ELF in Multiwfn for multi-configuration wavefunctions. The following paragraphs are my response to a Multiwfn user in an E-mail, which justifies the definition of ELF employed by Multiwfn that given above:

It should be noted that ELF has different interpretations, or there are different ways of deriving ELF. Originally, Becke et al. derived ELF based on average Fermi hole for HF wavefunction, therefore the original form of ELF is incompatible with multi-configuration wavefunctions like MP2, CASSCF and CCSD (it is not even strictly applicable to DFT wavefunction!). I noticed that some people extended ELF to multi-configuration wavefunctions by following Becke's original idea, which, however, is fairly complicated and difficult to implement, and thus not adopted by Multiwfn.

In contrast, for multi-configuration wavefunctions, Multiwfn calculates ELF simply based on their natural orbitals, as mentioned in Section 2.6 of manual. I believe this is never incorrect but physically sound. In Angew. Chem. Int. Ed. Engl., 31, 187 (1992), Savin et al. reinterpreted the physical meaning of ELF in terms of kinetic energy, this paves a way of calculating ELF for any kind of wavefunction as long as kinetic energy density is available. By following this viewpoint, all three kinds of kinetic energy densities involved in the expression of ELF are calculated based on natural orbitals. Although the result of this implementation is not identical as the one mentioned above (i.e. based on Fermi hole exactly derived from multi-configuration wavefunctions), it is also strict and physically meaningful.

In short, for multi-configuration wavefunctions, the realization of ELF is not unique, which depends on how one looks at ELF.

Notice that there is a deficiency of ELF, sometimes with r going beyond from molecular boundary, D(r) decreases faster than \(D_{0}\)(r) and then ELF reaches 1 (completely localized). To overcome the problem, Multiwfn automatically adds a minimal value 10^-5 to D(r), this treatment almost does not affect the ELF value in interesting regions. You can also disable this treatment by modifying “ELF_addminimal” in settings.ini to 0.

Tsirelson and Stash put forward an approximate version of ELF in Chem. Phys. Lett., 351, 142 (2002), in which the actual kinetic energy term in D(r) is replaced by Kirzhnits type second-order gradient expansion, that is

\[(1/2)\sum_i\eta_i\left|\nabla\varphi_i(\mathbf{r})\right|^2\approx\tau_{\mathrm{TA}}(\mathbf{r})+(1/72)\left|\nabla\rho(\mathbf{r})\right|^2/\rho(\mathbf{r})+(1/6)\nabla^2\rho(\mathbf{r})\]

so that ELF is totally independent from wavefunction, and then can be used to analyze electron density from X-ray diffraction data. Of course, Tsirelson’s ELF can also be used to analyze electron density from quantum chemistry calculation, but is not as good as the ELF defined by Becke owing to the approximation introduced in kinetic energy term, however, qualitative conclusions can still be recovered in general. If you want to use Tsirelson’s definition of ELF, change “ELFLOL_type” in settings.ini from 0 to 1.

If “ELFLOL_type” is set to 2, another formalism will be used:

\[\frac{1}{1+D(\mathbf{r})/D_0(\mathbf{r})}\]

A real space function closely related to ELF is SCI index, it is very useful for identifying strong covalent bonds, see introduction of user-defined function 37 in Section 2.7.

10 Localized orbital locator (LOL)

This is another function for locating high localization regions likewise ELF, defined by Schmider and Becke in the paper J. Mol. Struct. (THEOCHEM), 527, 51 (2000):

\[\mathrm{LOL}(\mathbf{r})=\frac{\tau(\mathbf{r})}{1+\tau(\mathbf{r})}\]

where

\[\tau(\mathbf{r})=\frac{D_0(\mathbf{r})}{(1/2)\sum_i\eta_i\left|\nabla\varphi_i(\mathbf{r})\right|^2}\]

\(D_{0}\)(r) for spin-polarized system and closed-shell system are defined in the same way as in ELF.

LOL has a similar expression compared to ELF. Actually, the chemically significant regions that are highlighted by LOL and ELF are generally qualitatively comparable. It is highly recommended to read my review Angew. Chem. Int. Ed., 137, e202504895 (2025) DOI: 10.1002/anie.202504895, which exemplified the important practical value of LOL. Obviously, LOL can be interpreted in kinetic energy way as for ELF, however, LOL can also be interpreted in view of localized orbital. Small (large) LOL value usually appears in boundary (inner) region of localized orbitals because the gradient of orbital wavefunction is large (small) in this area. The value range of LOL is identical to ELF, namely [0,1].

Multiwfn also supports the approximate version of LOL defined by Tsirelson and Stash (Acta. Cryst., B58, 780 (2002)), namely the actual kinetic energy term in LOL is replaced by second-order gradient expansion, as what they do for ELF. This Tsirelson’s version of LOL can be activated by setting “ELFLOL_type” to 1.

For special reason, if “ELFLOL_type” in settings.ini is changed from 0 to 2, another formalism will be used: \(\mathrm{LOL}(\mathbf{r})=\frac{1}{1+[1/\tau(\mathbf{r})]^2}\).

11 Local information entropy

Information entropy is a quantification of information, this theory was proposed by Shannon in his study of information transmission in noise channel, nowadays its application has been largely widened to other areas, including theoretical chemistry. For example, Aslangul and coworkers attempted to decompose diatomic and triatomic molecules into mutually exclusive space by minimizing information entropy (Adv. Quantum Chem., 6, 93 (1972)), Parr et al. discussed the relationship between information entropy and atom partition as well as molecular similarity (J. Phys. Chem. A, 109, 3957 (2005)), Noorizadeh and Shakerzadeh suggested using information entropy to study aromaticity (Phys. Chem. Chem. Phys., 12, 4742 (2010)). The formula of Shannon’s information entropy for normalized and continuous probability function is

\[S(x)=-\int P(x)\ln P(x)\,\mathrm{d}x\]

For chemical system, if P(x) is replaced with ρ(r)/N, then the integrand may be called as local information entropy of electrons

\[S(\mathbf{r})=-\frac{\rho(\mathbf{r})}{N}\ln\frac{\rho(\mathbf{r})}{N}\]

where N is the total number of electrons in current system. Integrating this function over whole space yields information entropy.

12 Total electrostatic potential (ESP)

\[V_{\mathrm{ESP}}(\mathbf{r})=V_{\mathrm{nuc}}(\mathbf{r})+V_{\mathrm{ele}}(\mathbf{r})=\sum_A\frac{Z_A}{\left|\mathbf{r}-\mathbf{R}_A\right|}-\int\frac{\rho(\mathbf{r}')}{\left|\mathbf{r}-\mathbf{r}'\right|}\mathrm{d}\mathbf{r}'\]

where Z is nuclear charge, the index A cycles over all atoms. If pseudopotential is used, then \(Z_{A}\) is the number of explicitly represented electrons. The atomic unit (a.u.) of ESP is Hartree/e, where e is elementary charge. Other commonly used units are eV/e and kcal/(mol⋅e), however for simplicity, in Multiwfn they are printed as eV and kcal/mol, respectively. The value in eV/e unit is equal to the value in SI unit (J/C).

This function measures the electrostatic interaction energy between a unit point charge placed at r and the present without consideration of charge polarization and transfer. A positive (negative) value implies that current position is dominated by nuclear (electronic) charges. Molecular electrostatic potential (ESP) has been widely used for prediction of nucleophilic and electrophilic sites for a long time. It is also valuable in studying hydrogen bonds, halogen bonds, molecular recognition and the intermolecular interaction of aromatics. Moreover, based on statistical analysis, Murray and coworkers found a set of functions called GIPF, see J. Mol. Struct. (THEOCHEM), 307, 55, which connects ESP in molecular surface and macroscopic properties. There are a lot of reviews on ESP, interested readers are suggested to consult WIREs Comput. Mol. Sci., 1, 153 (2011), Theor. Chem. Acc., 108, 134 (2002), Chapter 17 of the book Chemical Reactivity Theory-A Density Functional View, the entry "Electrostatic Potentials: Chemical Applications" (page 912) in the book Encyclopedia of Computational Chemistry and Chapter 7 of the book Reviews in Computational Chemistry vol.2.

By the way, if you only want to obtain the electrostatic potential contributed by electrons, namely \(V_{\mathrm{ele}}\)(r), you can use the 14th user-defined function. If you want to omit contribution of specific nucleus during evaluating \(V_{\mathrm{ESP}}\)(r), use 39th user-defined function. See corresponding entries in Section 2.7 for details.

Note 1

Currently there are two sets of internal code in Multiwfn for evaluating ESP, they can be selected by "iESPcode" in settings.ini. iESPcode=1 corresponds to the old and slow code written by Tian Lu, iESPcode=2 (default) corresponds to the newer and fast code provided by Jun Zhang based on his LIBRETA electron repulsion integral library. In addition, in many kinds of analyses, if you are using .fch/fchk as input file and Gaussian has been installed on your machine, you can allow Multiwfn to invoke the cubegen utility in Gaussian package to calculate ESP, the speed is even faster than "iESPcode=2" if number of your CPU cores is less than 10, see Section 5.7 for detail.

Note 2

Some functions in Multiwfn modify orbital occupations. For example, users can use subfunction 26 of main function 6 to manually modify occupation number of specific orbitals. Only the ESP calculated by internal code with "iESPcode=1" directly reflects the effect of the modification on occupation. In the case of using "iESPcode=2", in order to represent the effect of modification in the ESP calculation, you should first export the modified wavefunction to .wfn/wfx/molden/mwfn (via subfunction2 of main function 100), then reboot Multiwfn and load the file.

13 Reduced density gradient (RDG)

RDG and sign(λ₂)ρ are a pair of very important functions for revealing weak interaction regions, they are collectively employed in NCI method, see J. Am. Chem. Soc., 132, 6498 (2010) for detail, and application examples in Section 3.23.1. RDG is defined as

\[\mathrm{RDG}(\mathbf{r})=\frac{1}{2(3\pi^2)^{1/3}}\frac{\left|\nabla\rho(\mathbf{r})\right|}{\rho(\mathbf{r})^{4/3}}\]

Notice that there is a parameter “RDG_maxrho” in settings.ini, if the value is set to x, then RDG function will be set to an arbitrary big value (100.0) where the electron density is larger than x. This mechanism allows uninteresting regions to be shielded when viewing isosurfaces of weak interaction regions. By default x is 0.05, you can nullify this treatment by setting the parameter to zero.

14 Reduced density gradient (RDG) with promolecular approximation

Weak interaction has significant influence on conformation of macromolecules, binding mode of proteins and ligands; however, reproduction of electron density by ab initio and grid data calculation of RDG for such huge systems are always too time-consuming. Fortunately, it is found that weak interaction analysis under promolecular density is still reasonable. Promolecular density is simply constructed by superposing electron densities of free-state atoms and hence can be evaluated extremely rapidly

\[\rho^{\mathrm{pro}}(\mathbf{r})=\sum_A\rho_A^{\mathrm{free,fit}}(\mathbf{r}-\mathbf{R}_A)\]

where \({\rho}_{A}^{\mathrm{free,fit}}(\mathbf{r})\) is pre-fitted spherically averaged electron density of atom A. The atomic densities for H~Lr are built-in data of Multiwfn, among which the data for H~Ar are taken from supplemental material of J. Am. Chem. Soc., 132, 6498 (2010), while those for other elements are evaluated according to the description in Appendix 3. For elements heavier than Lr the promolecular approximation is not currently available.

Note

For efficiency consideration, if contribution from H, C, N or O atom to the function value at a specific point is less than 0.00001, then the contribution will not be calculated, for huge system this treatment improves efficiency several times and the result is almost unperturbed. You can also disable this treatment by setting “atomdenscut” in settings.ini to 0.

The parameter “RDGprodens_maxrho” in settings.ini is the counterpart of “RDG_maxrho” in the case of promolecular approximation.

15 sign(λ₂)ρ

\[\Omega(\mathbf{r})=sign[\lambda_2(\mathbf{r})]\rho(\mathbf{r})\]

where sign[λ₂(r)] means the sign of the second largest eigenvalue of electron density Hessian matrix at position r.

16 sign(λ₂)ρ with promolecular approximation

The actual electron density used to evaluate sign[λ₂(r)] is approximated by promolecular electron density, see the description in real space function 14.

17 Exchange-correlation density, correlation hole and correlation factor

These functions involve advanced topics, in order to clarify their physical meanings and to avoid confusion of the symbols used in a wide variety of papers, I think it is worth to use many more texts to introduce theoretical background. For more detailed discussions, please consult Chapter 2 of A Chemist's Guide to Density Functional Theory 2ed and Chapter 5 of Methods of molecular quantum mechanics (2ed, McWeeny).

======= Theoretical basis ======

Pair density is defined as

\[\pi(\mathbf{r}_1,\mathbf{r}_2)=N(N-1)\int\cdots\int\left|\Psi(\mathbf{x}_1,\mathbf{x}_2,...\mathbf{x}_N)\right|^2\mathrm{d}\sigma_1\mathrm{d}\sigma_2\mathrm{d}\mathbf{x}_3\mathrm{d}\mathbf{x}_4...\mathrm{d}\mathbf{x}_N\]

where Ψ is system wavefunction, r is space coordinate, σ is spin coordinate, x is space-spin coordinate. Pair density denotes the probability that finding an electron at \(r_{1}\) and another electron at \(r_{2,}\) regardless of the spin type. If we perform double-integration for pair density over the whole space, we will get N(N-1), reflecting the nature that there are N(N-1) electron pairs in present system. Obviously, pair density can be decomposed to contributions from different spin types of electron pairs

\[\pi(\mathbf{r}_1,\mathbf{r}_2)=\pi^{\alpha\alpha}(\mathbf{r}_1,\mathbf{r}_2)+\pi^{\alpha\beta}(\mathbf{r}_1,\mathbf{r}_2)+\pi^{\beta\alpha}(\mathbf{r}_1,\mathbf{r}_2)+\pi^{\beta\beta}(\mathbf{r}_1,\mathbf{r}_2)\]

If the electron motions are completely independent with each other, then the probability density of finding two electrons with spin σ₁ at \(r_{1}\) and with σ₂ at \(r_{2}\) respectively should simply be \(\rho^{\sigma_1}(\mathbf{r}_1)\rho^{\sigma_2}(\mathbf{r}_2)\). Of course, in real world electrons are always interacting with each other, so their motions are correlated. The pair density thereby should be corrected by exchange-correlation density Г

\[\pi^{\sigma_1\sigma_2}(\mathbf{r}_1,\mathbf{r}_2)=\rho^{\sigma_1}(\mathbf{r}_1)\rho^{\sigma_2}(\mathbf{r}_2)+\Gamma_{\mathrm{XC}}^{\sigma_1\sigma_2}(\mathbf{r}_1,\mathbf{r}_2)\]

If we have already known that an electron with spin σ₁ presents at \(r_{1}\), then the probability of finding another electron with spin σ₂ at \(r_{2}\) is known as conditional probability (This function is also known as Lennard-Jones function)

\[\Omega^{\sigma_1\sigma_2}(\mathbf{r}_1,\mathbf{r}_2)=\frac{\pi^{\sigma_1\sigma_2}(\mathbf{r}_1,\mathbf{r}_2)}{\rho^{\sigma_1}(\mathbf{r}_1)}=\rho^{\sigma_2}(\mathbf{r}_2)+\frac{\Gamma_{\mathrm{XC}}^{\sigma_1\sigma_2}(\mathbf{r}_1,\mathbf{r}_2)}{\rho^{\sigma_1}(\mathbf{r}_1)}\]

Correlation hole reveals the decrease of probability of finding another electron with spin σ₂ at \(r_{2}\) when an electron with spin σ₁ presents at \(r_{1}\) owing to electron correlation effect

\[h_{\mathrm{XC}}^{\sigma_1\sigma_2}(\mathbf{r}_1,\mathbf{r}_2)=\frac{\Gamma_{\mathrm{XC}}^{\sigma_1\sigma_2}(\mathbf{r}_1,\mathbf{r}_2)}{\rho^{\sigma_1}(\mathbf{r}_1)}=\Omega^{\sigma_1\sigma_2}(\mathbf{r}_1,\mathbf{r}_2)-\rho^{\sigma_2}(\mathbf{r}_2)=\frac{\pi^{\sigma_1\sigma_2}(\mathbf{r}_1,\mathbf{r}_2)}{\rho^{\sigma_1}(\mathbf{r}_1)}-\rho^{\sigma_2}(\mathbf{r}_2)\]

Correlation factor is a function closely related to correlation hole

\[f_{\mathrm{XC}}^{\sigma_1\sigma_2}(\mathbf{r}_1,\mathbf{r}_2)=\frac{h_{\mathrm{XC}}^{\sigma_1\sigma_2}(\mathbf{r}_1,\mathbf{r}_2)}{\rho^{\sigma_2}(\mathbf{r}_2)}=\frac{\Gamma_{\mathrm{XC}}^{\sigma_1\sigma_2}(\mathbf{r}_1,\mathbf{r}_2)}{\rho^{\sigma_1}(\mathbf{r}_1)\rho^{\sigma_2}(\mathbf{r}_2)}\]

Collectively, one can write out

\[\pi^{\sigma_1\sigma_2}(\mathbf{r}_1,\mathbf{r}_2)=\rho^{\sigma_1}(\mathbf{r}_1)\rho^{\sigma_2}(\mathbf{r}_2)+\rho^{\sigma_1}(\mathbf{r}_1)h_{\mathrm{XC}}^{\sigma_1\sigma_2}(\mathbf{r}_1,\mathbf{r}_2)=\rho^{\sigma_1}(\mathbf{r}_1)\rho^{\sigma_2}(\mathbf{r}_2)[1+f_{\mathrm{XC}}^{\sigma_1\sigma_2}(\mathbf{r}_1,\mathbf{r}_2)]\]

Г_XC can be decomposed to the sum of exchange correlation (also called as Fermi correlation) part Г_X and Coulomb correlation part Г_C, therefore \(h_{\mathrm{XC}}\) can be straightforwardly decomposed to exchange hole \(h_{X}\) (also called as Fermi hole) and Coulomb hole \(h_{C}\) as follows. Likewise, \(f_{\mathrm{XC}}\) can be decomposed to \(f_{X}\) and \(f_{C}\)

\[\begin{gathered}\Gamma_{\mathrm{XC}}^{\sigma_1\sigma_2}(\mathbf{r}_1,\mathbf{r}_2)=\Gamma_{\mathrm{X}}^{\sigma_1\sigma_2}(\mathbf{r}_1,\mathbf{r}_2)+\Gamma_{\mathrm{C}}^{\sigma_1\sigma_2}(\mathbf{r}_1,\mathbf{r}_2)\\ h_{\mathrm{XC}}^{\sigma_1\sigma_2}(\mathbf{r}_1,\mathbf{r}_2)=h_{\mathrm{X}}^{\sigma_1\sigma_2}(\mathbf{r}_1,\mathbf{r}_2)+h_{\mathrm{C}}^{\sigma_1\sigma_2}(\mathbf{r}_1,\mathbf{r}_2)\\ f_{\mathrm{XC}}^{\sigma_1\sigma_2}(\mathbf{r}_1,\mathbf{r}_2)=f_{\mathrm{X}}^{\sigma_1\sigma_2}(\mathbf{r}_1,\mathbf{r}_2)+f_{\mathrm{C}}^{\sigma_1\sigma_2}(\mathbf{r}_1,\mathbf{r}_2)\end{gathered}\]

Fermi correlation only presents between like-spin electrons; while Coulomb correlation occurs between any two electrons. Fermi correlation is much more important than Coulomb correlation, even at Hartree-Fock level, Fermi correlation is always well represented due to the anti-symmetry requirement of Slater determinant, while Coulomb is completely omitted. Only post-HF wavefunction is capable of simultaneously exhibiting Fermi and Coulomb correlation effects. Commonly we only focus on Fermi hole while neglecting Coulomb hole. One can easily show that integration of \(h_{X}\) over whole space is exactly equal to -1, hence Fermi correlation perfectly avoided self-pairing problem, which may cause significant rise in system energy; while the integration for \(h_{C}\) is zero, this is mainly why Coulomb correlation has less influence on system energy.

It is also rather straightforward to obtain the pair density and conditional probability when only exchange correlation or Coulomb correlation is taken into account.

The total π, Ω, Г_XC (or Г_X, Г_C) and \(h_{\mathrm{XC}}\) (or \(h_{X}\), \(h_{C}\)) for an electron with spin σ, regardless the spin of another electron can be defined as

\[\begin{gathered}\pi^{\sigma,\mathrm{tot}}(\mathbf{r}_1,\mathbf{r}_2)=\pi^{\alpha\alpha}(\mathbf{r}_1,\mathbf{r}_2)+\pi^{\alpha\beta}(\mathbf{r}_1,\mathbf{r}_2)\\ \Omega^{\sigma,\mathrm{tot}}(\mathbf{r}_1,\mathbf{r}_2)=\Omega^{\alpha\alpha}(\mathbf{r}_1,\mathbf{r}_2)+\Omega^{\alpha\beta}(\mathbf{r}_1,\mathbf{r}_2)\\ \Gamma_{\mathrm{XC}}^{\sigma,\mathrm{tot}}(\mathbf{r}_1,\mathbf{r}_2)=\Gamma_{\mathrm{XC}}^{\alpha\alpha}(\mathbf{r}_1,\mathbf{r}_2)+\Gamma_{\mathrm{XC}}^{\alpha\beta}(\mathbf{r}_1,\mathbf{r}_2)\\ h_{\mathrm{XC}}^{\sigma,\mathrm{tot}}(\mathbf{r}_1,\mathbf{r}_2)=h_{\mathrm{XC}}^{\alpha\alpha}(\mathbf{r}_1,\mathbf{r}_2)+h_{\mathrm{XC}}^{\alpha\beta}(\mathbf{r}_1,\mathbf{r}_2)\end{gathered}\]

======= Technical aspects ======

For single-determinant wavefunctions, exchange-correlation density for an α electron can be explicitly written as

\[\Gamma_{\mathrm{XC}}^{\alpha,\mathrm{tot}}(\mathbf{r}_1,\mathbf{r}_2)=\Gamma_{\mathrm{XC}}^{\alpha\alpha}(\mathbf{r}_1,\mathbf{r}_2)=-\sum_{i\in\alpha}^{\mathrm{occ}}\sum_{j\in\alpha}^{\mathrm{occ}}\varphi_i^*(\mathbf{r}_1)\varphi_j^*(\mathbf{r}_2)\varphi_j(\mathbf{r}_1)\varphi_i(\mathbf{r}_2)\]

Note that \(\Gamma_{\mathrm{XC}}^{\alpha\beta}(\mathbf{r}_1,\mathbf{r}_2)=\Gamma_{\mathrm{XC}}^{\beta\alpha}(\mathbf{r}_1,\mathbf{r}_2)=0\) for this type of wavefunction. To obtain the expression for β electron, just replace α by β, similarly hereinafter.

For post-HF wavefunction, exact evaluation of pair density requires two-particle density matrix (2PDM). Unfortunately, 2PDM is very difficult to obtain, mainstream quantum chemistry packages including Gaussian cannot directly output it. In Multiwfn, exchange-correlation density for post-HF wavefunction is approximately evaluated by natural orbital formalism. Note that the approximate method is not unique, see J. Chem. Theory Comput., 6, 2736 (2010) for discussion. The most popular form among them, which is firstly derived by Müller, is currently implemented in Multiwfn as below. See Phys. Lett. A, 105, 446 (1984), also see Mol. Phys., 100, 401 (2002) for extensive discussion (especially equation 32).

\[\Gamma_{\mathrm{XC,approx}}^{\alpha,\mathrm{tot}}(\mathbf{r}_1,\mathbf{r}_2)=-\sum_{i\in\alpha}\sum_{j\in\alpha}\sqrt{\eta_i\eta_j}\varphi_i^*(\mathbf{r}_1)\varphi_j^*(\mathbf{r}_2)\varphi_j(\mathbf{r}_1)\varphi_i(\mathbf{r}_2)\]

Obviously, if occupation numbers of natural spin orbitals are integer (0 or 1), then \({\Gamma}_{\mathrm{XC,approx}}^{\alpha\mathrm{,tot}}\) reduces to single-determinant form. so \({\Gamma}_{\mathrm{XC,approx}}^{\alpha\mathrm{,tot}}\) can be regarded as a general form to evaluate exchange-correlation density. Note that post-HF wavefunction has taken Coulomb correlation between unlike-spin electrons into account, however there is no way to separate \({\Gamma}_{\mathrm{XC,approx}}^{\alpha\alpha}\) and \({\Gamma}_{\mathrm{XC,approx}}^{\alpha\beta}\) from \({\Gamma}_{\mathrm{XC,approx}}^{\alpha\mathrm{,tot}}\).

The exchange-only part of Г for post-HF wavefunction can be approximately evaluated as below (of course, there is no exchange correlation between αβ electron pair)

\[\Gamma_{\mathrm{X,approx}}^{\alpha,\mathrm{tot}}(\mathbf{r}_1,\mathbf{r}_2)=\Gamma_{\mathrm{X,approx}}^{\alpha\alpha}(\mathbf{r}_1,\mathbf{r}_2)=-\sum_{i\in\alpha}\sum_{j\in\alpha}\eta_i\eta_j\varphi_i^*(\mathbf{r}_1)\varphi_j^*(\mathbf{r}_2)\varphi_j(\mathbf{r}_1)\varphi_i(\mathbf{r}_2)\]

So Coulomb-only part of Г can be evaluated as (including both αα and αβ pair contributions)

\[\begin{gathered}\Gamma_{\mathrm{C,approx}}^{\alpha,\mathrm{tot}}(\mathbf{r}_1,\mathbf{r}_2)=\Gamma_{\mathrm{XC,approx}}^{\alpha,\mathrm{tot}}(\mathbf{r}_1,\mathbf{r}_2)-\Gamma_{\mathrm{X,approx}}^{\alpha\alpha}(\mathbf{r}_1,\mathbf{r}_2)\\ \quad\quad=\sum_{i\in\alpha}\sum_{j\in\alpha}\left[(\eta_i\eta_j-\sqrt{\eta_i\eta_j})\varphi_i^*(\mathbf{r}_1)\varphi_j^*(\mathbf{r}_2)\varphi_j(\mathbf{r}_1)\varphi_i(\mathbf{r}_2)\right]\end{gathered}\]

Since we already have explicit expression to calculate Г_XC term, other quantities introduced earlier can be easily computed according to the relationships between them and Г_XC. Recall that \(\rho^\sigma(\mathbf{r})=\sum_{i\in\sigma}\eta_i\left|\varphi_i(\mathbf{r})\right|^2\).

Postscript: One can show that \({\Gamma}_{\mathrm{XC,approx}}^{\alpha\mathrm{,tot}}({\mathbf{r}}_{1},{\mathbf{r}}_{2})\) also exactly holds the requirement that integration of \(r_{2}\) over the whole space is equal to \(-{\rho}^{\alpha}({\mathbf{r}}_{1})\). However, in common, integrating \(r_{2}\) over the whole space for \({\Gamma}_{\mathrm{X,approx}}^{\alpha\mathrm{,tot}}({\mathbf{r}}_{1},{\mathbf{r}}_{2})\) and \({\Gamma}_{\mathrm{C,approx}}^{\alpha\mathrm{,tot}}({\mathbf{r}}_{1},{\mathbf{r}}_{2})\) deviate from \(-{\rho}^{\alpha}({\mathbf{r}}_{1})\) and zero, respectively, which are basic properties of exact form of \({\Gamma}_{\mathrm{X}}^{\alpha\mathrm{,tot}}\) and \({\Gamma}_{\mathrm{C}}^{\alpha\mathrm{,tot}}\).

======= Usage ======

In Multiwfn, \(r_{1}\) is seen as reference point and \(r_{2}\) is seen as variable, to define the coordinate of reference point, just modifying “refxyz” in settings.ini before booting up.

“paircorrtype” parameter in settings.ini controls which type of correlation effect will be taken into consideration in calculation of Г. Since correlation hole and correlation factor are calculated based on Г, this setting also affects them. For single-determinant wavefunction, =1 and =3 are equivalent and =2 is meaningless, because Coulomb correlation is completely omitted.

=1: Only consider exchange correlation

=2: Only consider Coulomb correlation

=3: Consider both exchange and Coulomb correlation

“pairfunctype” parameter in settings.ini controls which function and which spin will be calculated by real space function 17, see below, those enclosed by parentheses are for single-determinant wavefunction cases. Of course, for closed-shell system, the results for α spin are exactly identical to those for β spin. Note that correlation factor for post-HF wavefunction case is undefined.

=1: \(h^{α,tot}\) (\(h^{α,α}\))

=2: \(h^{β,tot}\) (\(h^{β,β}\))

=4: Undefined (f ^α,α)

=5: Undefined (f ^β,β)

=7: Γ^α,tot (Γ^α,α)

=8: Γ^β,tot (Γ^β,β)

=10: π^α,α when paircorrtype =1 (π^α,α)

=11: π^β,β when paircorrtype =1 (π^β,β)

=12: π^any,any (π^any,any)

For example, if paircorrtype=1 and pairfunctype=2, for post-HF wavefunction, what will be calculated is \(h_X^{\beta,\mathrm{tot}}(\mathbf{r}_1,\mathbf{r}_2)\), which is equivalent to \(h_{\mathrm{X}}^{\beta\beta}(\mathbf{r}_1,\mathbf{r}_2)\) since βα electron pairs have no exchange correlation. This quantity can be interpreted as Fermi hole at \(r_{2}\) caused by an β electron present at \(r_{1}\).

18 Average local ionization energy

Average local ionization energy is written as (Can. J. Chem., 68, 1440 (1990))

\[\bar{I}(\mathbf{r})=\frac{\sum_i\rho_i(\mathbf{r})\left|\varepsilon_i\right|}{\rho(\mathbf{r})}\]

where ρ_i(r) and ε_i are the electron density function and orbital energy of the ith molecular orbital, respectively. Hartree-Fock and typical DFT functionals such as B3LYP are both suitable for computing \(\bar{I}\). Lower value of \(\bar{I}\) indicates that the electrons at this point are more weakly bounded. \(\bar{I}\) has widespread uses, for example revealing atomic shell structures, measuring electronegativity, predicting p\(K_{a}\), quantifying local polarizability and hardness, but the most important one may be predicting reactive site of electrophilic or radical attack. It was proved that the minima of ALIE on vdW surface are good indicators to reveal which atoms are more likely to be the preferential site of electrophilic or radical attack. There are also many potential uses of \(\bar{I}\) waiting for further investigation. Excellent reviews of \(\bar{I}\) have been given by Politzer et al, see J. Mol. Model., 16, 1731 (2010) and Chapter 8 of the book Theoretical Aspects of Chemical Reactivity (2007).

Since \(\bar{I}\) is dependent upon orbital energies, while orbital energy for post-HF wavefunction is undefined, therefore when post-HF wavefunction is used, \(\bar{I}\) will be simply outputted as zero everywhere.

If the parameter "iALIEdecomp" in settings.ini is set to 1, in main function 1, not only \(\bar{I}\) will be outputted, the contribution from each occupied MOs will also be outputted, the contribution due to MO i is defined as

\[\bar{I}_i(\mathbf{r})=\frac{\rho_i(\mathbf{r})\left|\varepsilon_i\right|}{\rho(\mathbf{r})}\]

19 Source function

Source function was proposed by Bader and Gatti, see Chem. Phys. Lett., 287, 233 (1998).

\[SF(\mathbf{r},\mathbf{r}')=-\frac{\nabla^2\rho(\mathbf{r}')}{4\pi\left|\mathbf{r}-\mathbf{r}'\right|}\]

It can be shown that

\[\rho(\mathbf{r})=\int SF(\mathbf{r},\mathbf{r}')\mathrm{d}\mathbf{r}'\]

where r' ranges entire space. This equation suggests that SF(r,r') represents the effect of electronic Laplacian at r' on electron density at r. If at r' the electron is concentrated (namely Laplacian is negative, also suggesting potential energy dominates kinetic energy), then r' will be a source for the electron density at r; conversely, if at r' the electron is depleted, then r' diminish the electron density at r. If the range of integration in above formula is restricted to a local region Ω and we get a value S(r,Ω), then S(r,Ω)/ρ(r)×100% can be regarded as the contribution from region Ω to the electron density at r. Source function has many uses, when it is used to discuss bonding problems, usually bond critical points are taken as r. A very comprehensive review of theoretical background and applications of source function is given by Gatti in Struct. & Bond., 147, 193 (2010).

In Multiwfn, source function has two modes: (1) If "srcfuncmode" in settings.ini is set to 1, then r' is regarded as variable, while r is regarded as fixed reference point, whose coordinate is determined by “refxyz” in settings.ini. This is default mode, useful to study effects of electronic Laplacian at everywhere on specific point (2) If "srcfuncmode" is set to 2, then r becomes variable and r' becomes reference point, this is useful to study effect of electronic Laplacian at specific point on everywhere.

When r=r', this function will return −∇²ρ(r')/0.001 to avoid numerical problems.

20,21 Electron delocalization range function EDR(r;d) and orbital overlap distance function D(r)

Note

Content of this section and all analysis code of EDR(r;d) and D(r) was kindly contributed by Arshad Mehmood and then slightly adapted by Tian Lu.

The electron delocalization range function EDR(r;d) (J. Chem. Phys., 141, 144104 (2014); J. Chem. Theory Comput., 12, 3185 (2016); Angew. Chem. Int. Ed., 56, 6878 (2017)) quantifies the extent to which electrons at point r in a wave function occupy orbital lobes of size d. EDR(r;d) is built from the nonlocal one-particle reduced density matrix (1-RDM) \(\gamma(\mathbf{r},\mathbf{r}')=\sum_i\eta_i\varphi(\mathbf{r})\varphi(\mathbf{r}')\) as

\[\begin{gathered}EDR(\mathbf{r};d)=\int g_d(\mathbf{r},\mathbf{r}')\gamma(\mathbf{r},\mathbf{r}')\mathrm{d}\mathbf{r}'\\ g_d(\mathbf{r},\mathbf{r}')=\left(\frac{2}{\pi d^2}\right)^{3/4}\rho^{-1/2}(\mathbf{r})\exp\left(-\frac{\left|\mathbf{r}-\mathbf{r}'\right|^2}{d^2}\right)\end{gathered}\]

Here ρ(r) is the electron density at point r. The prefactor ensures that the EDR is between -1 and +1. The Multiwfn implementation evaluates the EDR on grids, for a single global input value of distance d. Section 4.5.6 illustrates an example.

At each point, the orbital overlap distance function D(r)=argma\(x_{d}\)EDR(r;d) corresponds to the distance d that maximizes EDR(r;d). Compact, chemically "hard" regions of small D(r) are distinguished from diffuse, chemically "soft" regions of large D(r). Atomic averages of valence-electron D(r) complement the information obtained from atomic partial charges. Plots of D(r) on density isosurfaces, and quantitative analysis of such surfaces, complement molecular electrostatic potentials. The Multiwfn implementation evaluates EDR(r;\(d_{i}\)) on a grid of distances \(d_{i}\), then uses a three-point numerical fit to find the maximum. Sections 4.5.7 and 4.12.8 illustrate example calculations.

22, 23 δg function defined in Independent Gradient Model (IGM) method

The δg function is defined as below in the original paper of Independent Gradient Model (IGM) method (Phys. Chem. Chem. Phys., 19, 17928 (2017)):

\[\delta g(\mathbf{r})=g^{\mathrm{IGM}}(\mathbf{r})-g(\mathbf{r})=\left|\sum_A \mathrm{abs}[\nabla\rho_A(\mathbf{r})]\right|-\left|\sum_A\nabla\rho_A(\mathbf{r})\right|\]

where ρ_A stands for atomic density of atom A, the abs() operator makes each of the three gradient components to be its absolute value. Specifically, Multiwfn supports two ways to calculate the δg:

• Real space function 22: Corresponding to the δg calculated under promolecular approximation, namely the ρ_A corresponds to spherically density of atom A in its isolated state, in this case only geometry information is needed to be supplied by input file, since the atomic densities are built-in data in Multiwfn (see Appendix 3 for detail).

• Real space function 23: Corresponding to the δg calculated based on Hirshfeld partition of actual molecular electron density, thus wavefunction must be supplied by input file. This definition was proposed by me in J. Comput. Chem., 43, 539 (2022). In this case the ρ_A is defined as \(w_{A}\)ρ, where ρ is molecular electron density calculated in usual way, \(w_{A}\) is Hirshfeld weighting function of atom A (see Section 3.9.1 for its definition).

The δg function at bond critical point in weak interaction region is shown to be closely related to interaction strength. This function can also be plotted as plane map or isosurface map to reveal all bonding regions. Very detailed introduction of IGM method is given in Section 3.23.5, and if you want to study weak interaction due to specific two or more fragments, you should use the function described in this section.

24 Interaction region indicator (IRI)

IRI was proposed by me in Chemistry—Methods, 1, 231 (2021), which is extremely useful in revealing all kinds of interaction regions of chemical system. See Section 3.23.8 for detailed information. IRI is defined as

\[IRI(\mathbf{r})=\frac{\left|\nabla\rho(\mathbf{r})\right|}{[\rho(\mathbf{r})]^a}\]

where a corresponds to "uservar" in settings.ini. When it is set to 0, then the recommended value 1.1 is employed. See Section 3.23.8 for detailed introduction of IRI. Note that IRI is set to an arbitrarily large value (5.0) if ρ is equal or smaller than “IRI_rhocut” in settings.ini, so that IRI isosurfaces in uninterested extremely low ρ regions will not occur. If “IRI_rhocut” is set to 0 then this treatment is not applied. For basin analysis and topology analysis for IRI or IRI-π, it should be set to 0 to avoid occurrence of artificial extrema due to this treatment (this is automatically done by Multiwfn).

25 van der Waals potential

Van der Waals (vdW) potential is very important for studying intermolecular interactions dominated by vdW effect, it has comparable role of electrostatic potential for intermolecular interactions dominated by electrostatic effect. See Section 3.23.7 for introduction of this function. The unit of this function is in kcal/mol, the probe atom can be set by "ivdwprobe" in settings.ini.