3.18 Fuzzy atomic space analysis (15)

3.18.0 Basic concepts

Before introducing each individual function, here I first introduce some basic concepts of fuzzy atomic space.

Atomic space is the local space attributed to a specific atom in the whole three-dimension molecular space. Below we will express atomic space as weighing function w. The methods used to partition the whole space into atomic spaces can be classified into two categories:

1 Discrete partition methods: The two representative methods are Bader's partition (also known as AIM partition) and Voronoi partition. They partition molecular space discretely, so any point can be attributed to only one atom, in other words,

\[\left\{\begin{gathered}w_A(\mathbf{r})=1\qquad \mathrm{if}\ \mathbf{r}\in\Omega_{\mathrm{A}}\\ w_A(\mathbf{r})=0\qquad \mathrm{if}\ \mathbf{r}\notin\Omega_{\mathrm{A}}\end{gathered}\right.\]

where Ω_A is atomic space of atom A.

2 Fuzzy partition methods: The representative methods include Hirshfeld, Becke, Hirshfeld-I, MBIS and ISA. They partition molecular space contiguously, atomic spaces overlap with each other, any point may be simultaneously attributed to many atoms to different extent, and the weights are normalized to unity. In other words, the two following conditions hold for all atoms and any point

\[\begin{gathered}0\leq w_A(\mathbf{r})\leq 1\qquad \forall A\\ \sum_B w_B(\mathbf{r})=1\end{gathered}\]

The most significant advantage of the fuzzy partition may be that the integration of real space function in fuzzy atomic space is much easier than in discrete atomic space. By using Becke's numerical DFT integration scheme (J. Chem. Phys., 88, 2547 (1988)), high accuracy of integration in fuzzy atomic space can be achieved for most real space functions at the expense of relatively low computation effort. In the fuzzy atomic space analysis module of Multiwfn, all integrations are realized by this scheme. The more integration points are used, the higher integration accuracy can be reached, one can adjust the number of points by "radpot" and "sphpot" parameter in settings.ini.

In fuzzy atomic space analysis module of Multiwfn, one can obtain many properties that based on fuzzy atomic spaces. Currently, the most widely used definitions of fuzzy atomic spaces, namely Hirshfeld, Hirshfeld-I and Becke are supported, they are introduced below. One can choose which fuzzy atomic spaces will be used by option -1.

Hirshfeld atomic space: In Theor. Chim. Acta (Berl.), 44, 129 (1977), Hirshfeld defined the atomic space as

\[w_A^{\mathrm{Hirsh}}(\mathbf{r})=\frac{\rho_A^{\mathrm{free}}(\mathbf{r}-\mathbf{R}_A)}{\sum_B\rho_B^{\mathrm{free}}(\mathbf{r}-\mathbf{R}_A)}\]

where R is coordinate of nucleus, ρ^free denotes spherically averaged atomic electron density in free-state.

In option -1, you will find two options "Hirshfeld" and "Hirshfeld*". The former uses atomic .wfn files to calculate the weights, they must be provided yourself or let Multiwfn automatically invoke Gaussian to generate them, see Section 3.7.3 for detail. The latter evaluates the weights directly based on built-in radial atomic densities and thus is more convenient, detail can be found in Appendix 3. I strongly suggest using "Hirshfeld*" instead of "Hirshfeld".

Hirshfeld-I (HI) atomic space: This is a well-known extension of Hirshfeld method, it was proposed in J. Chem. Phys., 126, 144111 (2007). Commonly the atomic space defined by HI is more physically meaningful than that of Hirshfeld, since it can respond to the actual molecular environment. Unfortunately, HI is much more expensive than Hirshfeld due to its iterative nature. Details of Hirshfeld-I and its implementation in Multiwfn have been introduced in Section 3.9.13 and thus will not be repeated here. When you choose HI in option -1, Multiwfn will first perform regular HI iterations (If you are confused by the operations, please consult the example of computing HI charges in Section 4.7.4). After HI atomic spaces have converged, you can do subsequent analyses.

MBIS atomic space: Like HI, MBIS refine atomic spaces iteratively. See Section 3.9.18 for details. When you choose MBIS in option -1, you will enter the interface of performing MBIS iteration, you should choose option 1 to start the atomic space refinement process. After MBIS atomic spaces have converged, you can do subsequent analyses.

Becke atomic space: First, consider a function p

\[p(d)=(3/2)d-(1/2)d^3\]

which can be iterated many times

\[\begin{gathered}f_1(d)=p(d)\\ f_2(d)=p[p(d)]\\ f_3(d)=p\{p[p(d)]\}\\ \cdots\end{gathered}\]

Then define a function s

\[s_k(t)=(1/2)[1-f_k(t)]\]

The plot of \(s_{k}\) versus to t is

From the graph above it can be seen that \(s_{k}\) gradually reduces from 1 to 0 with t varying from -1 to 1. The larger the k is, the sharper the curve becomes. The weighting function of Becke atomic space is based on simple transformation of \(s_{k}\), for details please consult original paper J. Chem. Phys., 88, 2547 (1988).

\[w_A^{\mathrm{Becke}}(\mathbf{r})=\frac{P_A(\mathbf{r})}{\sum_B P_B(\mathbf{r})}\]

\[\begin{gathered}P_A(\mathbf{r})=\prod_{B\neq A}s_k(v_{AB}(\mathbf{r}))\qquad v_{AB}(\mathbf{r})=\mu_{AB}(\mathbf{r})+a_{AB}(1-\mu_{AB}(\mathbf{r})^2)\\ a_{AB}=\frac{u_{AB}}{u_{AB}^2-1}\qquad u_{AB}=\frac{\chi_{AB}-1}{\chi_{AB}+1}\qquad \chi_{AB}=\frac{R_A^{\mathrm{cov}}}{R_B^{\mathrm{cov}}}\\ \left\{\begin{gathered}a_{AB}=-0.5\qquad \mathrm{if}\ a_{AB}<-0.5\\ a_{AB}=0.5\qquad \mathrm{if}\ a_{AB}>0.5\end{gathered}\right.\\ \mu_{AB}(\mathbf{r})=\frac{r_A-r_B}{R_{AB}}\qquad R_{AB}=\left|\mathbf{R}_A-\mathbf{R}_B\right|\qquad r_A=\left|\mathbf{r}-\mathbf{R}_A\right|\qquad r_B=\left|\mathbf{r}-\mathbf{R}_B\right|\end{gathered}\]

where R stands for coordinates of nucleus. \(R^{\mathrm{cov}}\) denotes covalent radius.

The number of iterations, namely k value, can be set by option -3. The default value (3) is appropriate for most cases. The definition of the covalent radius used to generate Becke atomic space can be chosen by option -2. Through corresponding suboptions, one can directly select a set of built-in radii (CSD radii, modified CSD radii, Pyykkö radii, Suresh radii, Hugo radii), load radii information from external plain text file (the format required is described in the program prompts), or modify current radii by manual input.

The origin paper of CSD radii is Dalton Trans., 2008, 2832, these radii were deduced from statistic of Cambridge Structural Database (CSD) for the elements with atomic numbers up to 96. Pyykkö radii were defined in Chem. Eur. J., 15, 186 (2008), which covers the entire periodic table, Groups 1–18, Z=1–118. Suresh radii were proposed in J. Phys. Chem. A, 105, 5940 (2001), which is based on theoretically calculated geometries of H₃C-E\(H_{n}\), the defined radii cover most of main group and transition elements in periodic table. Hugo radii were proposed in Chem. Phys. Lett., 480, 127 (2009), which has clear physical meaning and is based on atomic ionization energy. Notice that Hugo radius for hydrogen is rather large (even larger than Kr by 0.01 Bohr).

I found it is inappropriate to directly use any covalent radii definition shown above to define Becke's atomic space. The covalent radii of metal elements in IA and IIA groups are always large, e.g. CSD radius of lithium is 1.28 Å. While covalent radii of elements in such as VIIA group are always small, e.g. CSD radius of fluorine is only 0.58 Å. For main groups, the elements with small (large) covalent radius generally have large (small) electronegativity. So, in molecule environment, the atoms with small (large) covalent radius prefer to withdraw (donate) electrons to expand (shrink) their effective size, this behavior makes actual radii of main group elements in each row equalized. In order to faithfully reflect this behavior, I defined the so-called "modified CSD radii", namely the CSD radii of all main group elements (except for the first row) are replaced by CSD radii of the IVA group element in corresponding row, while transition elements still use their original CSD radii. The modified CSD radii are the default radii definition for Becke's atomic space.

The Becke atomic space of carbon in acetamide constructed by default parameters is illustrated below

3.18.1 Integration of a real space function in fuzzy atomic spaces (1)

This function is used to integrate real space function f in atomic spaces

\[I_A=\int_A w_A(\mathbf{r})f(\mathbf{r})\mathrm{d}\mathbf{r}\]

For example, if f is chosen as electron density, then \(I_{A}\) will be the electron population number of atom A.

f may be also chosen as the real space functions involving coordinates of two electrons, such as exchange-correlation density and source function. For this case, the coordinate of reference point can be set by option -10 (this is equivalent to set "refxyz" in settings.ini). If you have carried out topology analysis, you can also use option -11 to set a critical point as reference point, this is especially convenient for studying source function (for which bond critical point is usually set as reference point).

The "% of sum" and "% of sum abs" in output are defined as \(\left(I_A/\sum_B I_B\right)\times 100\%\) and \(\left(I_A/\sum_B\left|I_B\right|\right)\times 100\%\), respectively.

By default, all atomic spaces will be integrated. If you only need integral value of certain atoms, you can use option -5 to define the atom list.

3.18.2 Integration of a real space function in overlap spaces (8)

This function is used to integrate specified real space function f in overlap spaces between atomic pairs

\[I_{AB}=\int_A w_A(\mathbf{r})w_B(\mathbf{r})f(\mathbf{r})\mathrm{d}\mathbf{r}\]

For example, if f is chosen as electron density, then \(I_{\mathrm{AB}}\) will be the number of electrons shared by atom A and B. f may be also chosen as the real space functions involving coordinates of two electrons.

Integrals of positive and negative parts of f are outputted separately. Meanwhile, sum of diagonal elements \(\sum_{A} {I}_{AA}\), sum of non-diagonal elements \(\sum_{A} \sum_{B\neq A} {I}_{AB}\) and sum of all elements \(\sum_{A} \sum_{B} {I}_{AB}\) for positive and negative parts are also outputted together.

Currently only the fuzzy atomic space defined by Becke can be employed in this function.

3.18.3 Atomic and molecular multipole moments and <r²> (2)

This function is used to evaluate atomic and molecular monopole, dipole, quadrupole moments and octopole moments as well as <\(r^{2}\)>. All units in the output are in a.u.

In below formulae, superscript A means an atom named A. x, y and z are the components of electron coordinate r relative to nuclear coordinate R.

\[x=r_x-R_x^A\qquad y=r_y-R_y^A\qquad z=r_z-R_z^A\]

and \(r^{2}\) = \(x^{2}\) + \(y^{2}\) + \(z^{2}\).

Atomic monopole moment due to electrons is just negative of electron population number

\[p_A=-\int w_A(\mathbf{r})\rho(\mathbf{r})\mathrm{d}\mathbf{r}\]

Atomic charges are outputted together, namely \(q_{A}\) = \(p_{A}\) + \(Z_{A}\), where Z denotes nuclear charge.

Atomic dipole moment is useful to measure polarization of electron distribution around the atom, which is defined as

\[\unicode{x1D6CD}^A=\begin{bmatrix}\mu_x^A\\ \mu_y^A\\ \mu_z^A\end{bmatrix}=-\int\begin{bmatrix}x\\ y\\ z\end{bmatrix}w_A(\mathbf{r})\rho(\mathbf{r})\mathrm{d}\mathbf{r}\]

Its magnitude, or say its norm, is

\[\left|\unicode{x1D6CD}^A\right|=\sqrt{(\mu_x^A)^2+(\mu_y^A)^2+(\mu_z^A)^2}\]

Multiwfn also outputs the contribution of present atom to total molecular dipole moment, which is evaluated as \(q_{A}\)R + µ^A.

Traceless Cartesian form of atomic quadrupole moment tensor is defined as (see Section 1.8.7 of book The Quantum Theory of Atoms in Molecules-From Solid State to DNA and Drug Design).

\[\unicode{x1D6AF}^A=\begin{bmatrix}\Theta_{xx}^A & \Theta_{xy}^A & \Theta_{xz}^A\\ \Theta_{yx}^A & \Theta_{yy}^A & \Theta_{yz}^A\\ \Theta_{zx}^A & \Theta_{zy}^A & \Theta_{zz}^A\end{bmatrix}=-\frac{1}{2}\int\begin{bmatrix}3x^2-r^2 & 3xy & 3xz\\ 3yx & 3y^2-r^2 & 3yz\\ 3zx & 3zy & 3z^2-r^2\end{bmatrix}w_A(\mathbf{r})\rho(\mathbf{r})\mathrm{d}\mathbf{r}\]

whose magnitude can be calculated as

\[\left|\unicode{x1D6AF}^A\right|=\sqrt{(2/3)\left[(\Theta_{xx}^A)^2+(\Theta_{yy}^A)^2+(\Theta_{zz}^A)^2\right]}\]

Atomic quadrupole moments in Cartesian form can be used to exhibit deviation of electron distribution from spherical symmetry around nuclei. Specifically, \({\Theta}_{ii}^{A}\)<0 (\({\Theta}_{ii}^{A}\)>0) indicates that the electron density of atom A is elongated (contracted) along i direction. If the atomic electron density has exact spherical symmetry, then Θ_xx = Θ_yy = Θ_zz. Noticeably, the Cartesian quadrupole moment tensor Θ given here is traceless, namely the condition Θ_xx + Θ_yy + Θ_zz = 0 holds.

Standard Cartesian form of atomic quadrupole moment tensor is defined as follows. It is not outputted by default because it is rarely useful. However, if you hope it to be outputted, you can set “ispecial” in settings.ini to 1.

\[\unicode{x1D6AF}^A=\begin{bmatrix}\Theta_{xx}^A & \Theta_{xy}^A & \Theta_{xz}^A\\ \Theta_{yx}^A & \Theta_{yy}^A & \Theta_{yz}^A\\ \Theta_{zx}^A & \Theta_{zy}^A & \Theta_{zz}^A\end{bmatrix}=-\int\begin{bmatrix}x^2 & xy & xz\\ yx & y^2 & yz\\ zx & zy & z^2\end{bmatrix}w_A(\mathbf{r})\rho(\mathbf{r})\mathrm{d}\mathbf{r}\]

Inspired by electronic spatial extent (see Section 3.300.5), I defined atomic electronic spatial extent \(\left\langle {r}_{A}^{2} \right\rangle\), it is expressed as

\[\left\langle {r}_{A}^{2} \right\rangle=\int {r}^{2}{w}_{A}\left( \mathbf{r} \right)\rho\left( \mathbf{r} \right)\mathrm{d}\mathbf{r}=\left\langle {x}_{A}^{2} \right\rangle+\left\langle {y}_{A}^{2} \right\rangle+\left\langle {z}_{A}^{2} \right\rangle\]

whose X component is expressed as follows

\[\left\langle {x}_{A}^{2} \right\rangle=\int {x}^{2}{w}_{A}\left( \mathbf{r} \right)\rho(\mathbf{r}\mathrm{)d}\mathbf{r}\]

and similarly for \(\left\langle {y}_{A}^{2} \right\rangle\) and \(\left\langle {z}_{A}^{2} \right\rangle\).

\(\left\langle {r}_{A}^{2} \right\rangle\) is a useful metric of overall spatial extent of electron distribution within a fuzzy atom, while its Cartesian component reveals electronic spatial extent in specific direction.

The atomic quadrupole and octopole moments in spherical harmonic form are also outputted. The general expression of multipole moments in spherical harmonic form is

\[Q_{l,m}^A=-\int R_{l,m}^A(\mathbf{r})w_A(\mathbf{r})\rho(\mathbf{r})\mathrm{d}\mathbf{r}\]

All of the five components of quadrupole moment in spherical harmonic form correspond to

\[\begin{gathered}R_{2,0}=(3z^2-r^2)/2\\ R_{2,-1}=\sqrt{3}yz\qquad R_{2,1}=\sqrt{3}xz\\ R_{2,-2}=\sqrt{3}xy\qquad R_{2,2}=(\sqrt{3}/2)(x^2-y^2)\end{gathered}\]

All of the 7 components of octopole moment in spherical harmonic form correspond to

\[\begin{gathered}R_{3,0}=(1/2)(5z^2-3r^2)z\\ R_{3,-1}=\sqrt{3/8}(5z^2-r^2)y\qquad R_{3,1}=\sqrt{3/8}(5z^2-r^2)x\\ R_{3,-2}=\sqrt{15}xyz\qquad R_{3,2}=(\sqrt{15}/2)(x^2-y^2)z\\ R_{3,-3}=\sqrt{5/8}(3x^2-y^2)y\qquad R_{3,3}=\sqrt{5/8}(x^2-3y^2)x\end{gathered}\]

The magnitude of multipole moments in spherical harmonic form is calculated as

\[\left|Q_l^A\right|=\sqrt{\sum_m(Q_{l,m}^A)^2}\]

At the end of the calculation, the total number of electrons, molecular dipole moment and its magnitude are outputted. Molecular dipole moment is calculated as the sum of all contributions from atomic dipole moments and atomic charges (i.e. the sum of all "Contribution to molecular dipole moment" terms in the output information)

\[\unicode{x1D6CD}^{\mathrm{mol}}=\sum_A\left(q_A\mathbf{R}^A+\unicode{x1D6CD}^A\right)\]

In addition, Multiwfn outputs molecular quadrupole and octopole moments in Cartesian form and spherical harmonic form, they can also be viewed as sum of contributions of atoms. For example, molecular quadrupole moment of Θ_xy is expressed as

\[\Theta_{xy}=\frac{3}{2}\left[\sum_A R_x^A R_y^A Z_A-\sum_A\int xyw_A(\mathbf{r})\rho(\mathbf{r})\mathrm{d}\mathbf{r}\right]\]

where x, y, z in this context are Cartesian components of r with respect to (0,0,0) position. <\(r^{2}\)> of molecule can be written as

\[\left\langle {r}^{2} \right\rangle=\sum_{A} \int {r}^{2}{w}_{A}\left( \mathbf{r} \right)\rho\left( \mathbf{r} \right)\mathrm{d}\mathbf{r}\]

where r is radial distance with respect to (0,0,0).

By default, atomic multipole moments and <\(r^{2}\)> for all atoms are evaluated, and finally, these quantities of the whole system are printed. If you only need them for specific atoms, you can use option -5 to define an atom list, in this case only the quantities of selected atoms will be calculated and outputted. In addition, via this feature you can calculate the quantities of a molecule in a molecular complex, or calculate them of a fragment in a molecule, because in this case the "Molecular dipole and multipole moments" printed at the end of output are only contributed by the atoms in the defined list. See example in Section 4.15.3 for illustration of use of this feature.

After entering the present function, you will be asked to choose destination of outputting. If you choose 2 to output result to multipole.txt, a file named atom_moment.txt will also be produced in the current folder. Based on this file, atomic electric dipole and quadrupole moments can be visualized in VMD program via a special script, see Section 4.15.5 for detail.

PS 1

If your purpose is only calculating electric dipole/multipole moments and <\(r^{2}\)> for the whole system, it is best to use the function described in Section 3.300.5, it is significantly faster and more accurate since it calculates them analytically.

PS 2

If “ispecial” in settings.ini is set to 1, then the electron density involved in this function will be replaced with user-defined function. Via this feature, it is possible to realize some special purpose, such as calculating atomic dipole moments corresponding to variation of electron density, see #10 and relevant discussions in http://sobereva.com/wfnbbs/viewtopic.php?id=650.

3.18.4 Atomic overlap matrix and fragment overlap matrix (3, 33)

Subfunction 3 in fuzzy analysis module is used to calculate atomic overlap matrix (AOM) for orbitals in atomic spaces, the AOM will be outputted to AOM.txt in current folder. The element of AOM is defined as

\[S_{ij}(A)=\int_A\varphi_i(\mathbf{r})\varphi_j(\mathbf{r})\mathrm{d}\mathbf{r}\]

where i and j are orbital indices, integration is performed within fuzzy space of atom A. For unrestricted wavefunctions, AOMs between α orbitals and between β orbitals are outputted separately for each atom.

Notice that the highest virtual orbitals will not be taken into account during calculation. For example, present system has 10 orbitals in total, in which 7, 8, 9, 10 are not occupied, and user has set the occupation number of orbital 3 to zero via option 26 in main function 6, then the dimension of each AOM outputted by Multiwfn will be (6,6), corresponding to the overlap integral between the first 6 orbitals in each atomic space. If you hope to take all orbitals into account, set “ispecial” in settings.ini to 3.

Since orbitals are orthonormal in the whole space, in principle, summing up AOMs for all atoms (corresponding to integrating in the whole space) should yield an identity matrix

\[\mathbf{SUM}=\sum_A\mathbf{s}(A)=\mathbf{I}\]

Of course, this condition is not strictly held, because the integration is performed numerically. The deviation of SUM to identity matrix is a useful metric of integration accuracy

\[\mathrm{Error}=\frac{\sum_i\sum_j\left|\mathbf{SUM}_{i,j}-I_{i,j}\right|}{N_{\mathrm{atom}}}\]

Multiwfn automatically outputs the "Error" value. If it is not small enough, e.g. >0.001, then you may want to improve the integration accuracy via following ways

• (1) Enlarge "radpot" and "sphpot" in settings.ini

• (2) Set "radcut" in settings.ini to 0

• (3) Choose option -6 to change the default atomic integration grid to the much more expensive molecular integration grid

• (4) If diffuse functions were heavily employed, remove them

Fragment overlap matrix

Fragment overlap matrix (FOM) is simply sum of AOM of the atoms in the fragment. FOM of one fragment or two fragments can be calculated by subfunction 33 in fuzzy analysis module, the result is outputted to FOM.txt in current folder. You can directly define the atoms in the fragment(s) in this subfunction.

When atomic integration grid is used to evaluate AOM, and the number of atoms involved in the one or two fragments is significantly smaller than total number of atoms, the calculation cost of FOM is significantly lower than using subfunction 3 to calculate the entire AOM, because atoms not involved in the fragment(s) will simply be skipped.

3.18.5 Localization index (LI) and delocalization index (DI) (4, 44)

3.18.5.1 Theoretical background

Definition of LI and DI

For open-shell systems, the LI (λ) and DI (δ) are calculated for each spin of electrons respectively. Below only the expression of LI and DI for α electrons is given. For β electrons, just replacing α with β, similarly hereinafter. The electrons in atomic space A that can delocalize to atomic space B is computed as

\[\delta^\alpha(A\to B)=-\int_A\int_B\Gamma_{\mathrm{XC}}^{\alpha,\mathrm{tot}}(\mathbf{r}_1,\mathbf{r}_2)\mathrm{d}\mathbf{r}_1\mathrm{d}\mathbf{r}_2\]

where Г_XC is exchange-correlation density; if you are not familiar with it, please consult the discussion in part 17 of Section 2.6. The electrons in atomic space B that can delocalize to atomic space A is

\[\delta^\alpha(B\to A)=-\int_B\int_A\Gamma_{\mathrm{XC}}^{\alpha,\mathrm{tot}}(\mathbf{r}_1,\mathbf{r}_2)\mathrm{d}\mathbf{r}_1\mathrm{d}\mathbf{r}_2\]

Clearly, above two terms are identical in value, therefore we define DI between A and B as below, it measures the total number of α electrons shared by atom A and B

\[\delta^\alpha(A,B)=\delta^\alpha(A\to B)+\delta^\alpha(B\to A)=-2\int_A\int_B\Gamma_{\mathrm{XC}}^{\alpha,\mathrm{tot}}(\mathbf{r}_1,\mathbf{r}_2)\mathrm{d}\mathbf{r}_1\mathrm{d}\mathbf{r}_2\]

\(\lambda^\alpha\) measures the number of α electrons localized in an atom. Note that this quantity is not additive.

\[\lambda^\alpha(A)=-\int_A\int_A\Gamma_{\mathrm{XC}}^{\alpha,\mathrm{tot}}(\mathbf{r}_1,\mathbf{r}_2)\mathrm{d}\mathbf{r}_1\mathrm{d}\mathbf{r}_2=\delta^\alpha(A,A)/2\]

The relationship between LI, DI and the population number of electrons in atomic space is given below, the physical meaning is that the sum of α electrons of atom A that localized in atom A and that delocalized to other regions is the total number of α electrons in space A.

\[\begin{aligned}&\lambda^\alpha(A)+(1/2)\sum_{B\neq A}\delta^\alpha(A,B)\equiv\lambda^\alpha(A)+\sum_{B\neq A}\delta^\alpha(A\to B)=\\ &=-\int_A\int\Gamma_{\mathrm{XC}}^{\alpha,\mathrm{tot}}(\mathbf{r}_1,\mathbf{r}_2)=\int_A\rho^\alpha(\mathbf{r})\mathrm{d}\mathbf{r}=N_A^\alpha\end{aligned}\]

Using the Müller approximate expression of Г_XC, the DI and LI can be explicitly written as follows (δ in this form is also known as Fulton index, see Phys. Chem. Chem. Phys., 28, 19133 (2026) for a review)

\[\begin{gathered}\delta^\alpha(A,B)=2\sum_{i\in\alpha}\sum_{j\in\alpha}\sqrt{\eta_i\eta_j S_{ij}(A)S_{ij}(B)}\\ \lambda^\alpha(A)=\sum_{i\in\alpha}\sum_{j\in\alpha}\sqrt{\eta_i\eta_j S_{ij}(A)S_{ij}(A)}\end{gathered}\]

where S is atomic overlap matrix (AOM), see Section 3.18.4 for introduction.

Total DI and LI are the summation of α part and β part

\[\begin{gathered}\delta(A,B)=\delta^\alpha(A,B)+\delta^\beta(A,B)\\ \lambda(A)=\lambda^\alpha(A)+\lambda^\beta(A)\end{gathered}\]

Also worth noting is the DI of Ángyán-Loos-Mayer (ALM) formulation (J. Phys. Chem., 98, 5244 (1994)) as shown below, where orbitals are spatial with η within [0.0,2.0]. It is not implemented in Multiwfn.

\[\delta^\alpha(A,B)=\sum_i\sum_j\eta_i\eta_j S_{ij}(A)S_{ij}(B)\]

Special form of closed-shell cases

Since in closed-shell case \(\delta^\alpha(A,B)=\delta^\beta(A,B)\), one can evaluate total DI as

\[\delta(A,B)=2\delta^\alpha(A,B)=2\times 2\sum_m\sum_n\sqrt{\frac{\eta_m}{2}\frac{\eta_n}{2}}S_{mn}(A)S_{mn}(B)=2\sum_m\sum_n\sqrt{\eta_m\eta_n}S_{mn}(A)S_{mn}(B)\]

where m and n denote closed-shell natural orbitals. Similarly, the total LI for closed-shell cases is

\[\lambda(A)=\sum_m\sum_n\sqrt{\eta_m\eta_n}S_{mn}(A)S_{mn}(A)\]

For closed-shell systems, it is argued that the value of total DI is a quantitative measure of the number of electron pairs shared between two atoms. For example, total δ(A,B)=1.0 implies a pair of electrons (an α and a β electrons) is shared between atom A and B. (In fact, this is strictly true only for nonpolar bonds such as H-H bond in H₂. In polar bonds, the DI must be lower than formal bond, because what total DI actually reflects is the effective number of electron pairs shared by two atoms and thus somewhat reflects covalency. Note that the value of DI is very sensitive to the definition of atomic space employed.

Fragment LI and interfragment DI

Interfragment DI (IFDI) between fragments F and G can be evaluated as

\[\begin{aligned}&\delta^\alpha(F,G)=2\sum_{i\in\alpha}\sum_{j\in\alpha}\sqrt{\eta_i\eta_j S_{ij}(F)S_{ij}(G)}\\ &\to 2\sum_{i\in\alpha}\sum_{j\in\alpha}\sqrt{\eta_i\eta_j}\sum_{A\in F}S_{ij}(A)\sum_{B\in G}S_{ij}(B)\\ &\to \sum_{A\in F}\sum_{B\in G}\left[2\sum_{i\in\alpha}\sum_{j\in\alpha}\sqrt{\eta_i\eta_j S_{ij}(A)S_{ij}(B)}\right]\\ &\to \sum_{A\in F}\sum_{B\in G}\delta^\alpha(A,B)\end{aligned}\]

In Phys. Chem. Chem. Phys., 24, 11486 (2022), it was demonstrated that IFDI between two terminal fragments in a globally conjugated system is very useful in characterizing extent of global delocalization, and it IFDI is found to be well positively correlated with rotational barrier between the two fragments. This is because the stronger the original π conjugation is, the more obvious the destruction of conjugation will be when the two groups rotate relative to each other, and the higher the energy will rise.

Fragment LI (FLI) of fragment F can be evaluated as

\[\begin{aligned}&\lambda^\alpha(F)=\sum_{i\in\alpha}\sum_{j\in\alpha}\sqrt{\eta_i\eta_j S_{ij}(F)S_{ij}(F)}\\ &\to \sum_{A\in F}\sum_{B\in F}\sum_{i\in\alpha}\sum_{j\in\alpha}\sqrt{\eta_i\eta_j S_{ij}(A)S_{ij}(B)}\\ &\to \sum_{A\in F}\sum_{i\in\alpha}\sum_{j\in\alpha}\sqrt{\eta_i\eta_j S_{ij}(A)S_{ij}(A)}+\sum_{(B>A)\in F}2\sum_{i\in\alpha}\sum_{j\in\alpha}\sqrt{\eta_i\eta_j S_{ij}(A)S_{ij}(B)}\\ &\to \sum_{A\in F}\lambda^\alpha(A)+\sum_{(B>A)\in F}\delta^\alpha(A,B)\end{aligned}\]

Special form for single-determinant wavefunctions

For single-determinant wavefunction, because of integer occupation number of orbitals, DI and LI can be simplified as

\[\begin{gathered}\delta^\alpha(A,B)=2\sum_{i\in\alpha}^{\mathrm{occ}}\sum_{j\in\alpha}^{\mathrm{occ}}S_{ij}(A)S_{ij}(B)\\ \delta(A,B)=4\sum_m^{\mathrm{occ}}\sum_n^{\mathrm{occ}}S_{mn}(A)S_{mn}(B)\\ \lambda^\alpha(A)=\sum_{i\in\alpha}^{\mathrm{occ}}\sum_{j\in\alpha}^{\mathrm{occ}}S_{ij}(A)S_{ij}(A)\\ \lambda(A)=2\sum_m^{\mathrm{occ}}\sum_n^{\mathrm{occ}}S_{mn}(A)S_{mn}(A)\end{gathered}\]

Relationship between DI and fuzzy bond order

Conventionally, LI and DI are calculated in AIM atomic spaces (also called AIM basins). While in fuzzy atomic space analysis module of Multiwfn, they are calculated in fuzzy atomic space, the physical nature is the same. According to the discussion presented in J. Phys. Chem. A, 109, 9904 (2005) (compare Eq. 13 and Eq. 18), the DI calculated in fuzzy atomic space is just the so-called fuzzy bond order, which was defined by Mayer in Chem. Phys. Lett., 383, 368 (2004).

For closed-shell system, atomic valence can be calculated as the sum of its fuzzy bond orders

\[V(A)=\sum_{B\neq A}\delta(A,B)\]

Separation of σ and π contributions

For strictly planar molecules, because overlap integral of σ orbital and π orbital is exactly zero in atomic space, the contributions from σ and π electrons to DI can be exactly decomposed as DI-σ and DI-π

\[\begin{gathered}\delta^\alpha(A,B)=\delta_\sigma^\alpha(A,B)+\delta_\pi^\alpha(A,B)\\ \delta_\sigma^\alpha(A,B)=2\sum_{i\in\alpha}^{\sigma}\sum_{j\in\alpha}^{\sigma}\sqrt{\eta_i\eta_j S_{ij}(A)S_{ij}(B)}\\ \delta_\pi^\alpha(A,B)=2\sum_{i\in\alpha}^{\pi}\sum_{j\in\alpha}^{\pi}\sqrt{\eta_i\eta_j S_{ij}(A)S_{ij}(B)}\end{gathered}\]

Similarly, LI can be decomposed as LI-σ and LI-π. Summing up corresponding off-diagonal elements in DI-σ and DI-π matrix gives σ-atomic valence and π-atomic valence, respectively. If you want to compute DI/LI-σ (DI/LI-π), before the DI/LI calculation, you should set the occupation numbers of all π orbitals (σ orbitals) to zero by subfunction 26 of main function 6.

Covariance and relative fluctuation parameter

In some papers, especially the ones written by Bernard Silvi, the variance of electronic fluctuation in atomic space σ²(A) and the covariance of fluctuation of electron pair between two atomic spaces cov(A,B) are discussed. They are not directly outputted by Multiwfn, because there is a very simple relationship correlates σ²(A), cov(A,B) and DI(A,B), thus you can calculate them quite easily, see Chem. Rev., 105, 3911 (2005) for derivation

\[\begin{gathered}\mathrm{cov}(A,B)=-\delta(A,B)/2\\ \sigma^2(A)=N_A-\lambda(A)=-\sum_{B\neq A}\mathrm{cov}(A,B)=\sum_{B\neq A}\delta(A,B)/2\end{gathered}\]

where \(N_{A}\) is the electron population number in A. As mentioned above, the diagonal terms of the DI matrix outputted by Multiwfn are calculated as the sum of off-diagonal elements in the corresponding row (or column), hence you can simply obtain σ² by dividing corresponding diagonal term of DI matrix by two.

A quantity closely related to σ² is the relative fluctuation parameter introduced by Bader, which indicates the electronic fluctuations for a given atomic space relative to its electron population, you can calculate it manually if you want

\[\lambda_{\mathrm{F}}(A)=\sigma^2(A)/N_A\]

Alternatively, you can calculate below value to measure the proportion of the electrons localized in the atomic space

\[l(A)=\lambda(A)/N_A\]

3.18.5.2 Usage

In Multiwfn, before calculating LI and DI, AOM is calculated first automatically, this is the most time-consuming step. For open-shell systems, the LI and DI for α and β electrons, as well as for all electrons are outputted respectively. Notice that the diagonal terms of DI matrix are calculated as the sum of corresponding off-diagonal row (or column) elements. For closed-shell system, as stated above, they correspond to atomic valence.

3.18.6 Para-delocalization index (PDI) (5)

Para-delocalization index (PDI) is a quantity used to measure aromaticity of six-membered rings. PDI was first proposed in Chem. Eur. J., 9, 400 (2003), also see Chem. Rev., 105, 3911 (2005) for more discussion. PDI is essentially the averaged para-delocalization index (para-DI) in six-membered rings.

\[\mathrm{PDI}=\frac{\delta(1,4)+\delta(2,5)+\delta(3,6)}{3}\]

The basic idea behind PDI is that Bader and coworkers reported that DI in benzene is greater for para-related than for meta-related carbon atoms. Obviously, the larger the PDI, the larger the delocalization, and the stronger the aromaticity. The main limitation of the definition of PDI is that it can only be used to study aromaticity of six-membered rings, and it was shown that PDI is inappropriate for the cases when the ring plane has an out-plane distortion.

In Multiwfn, before calculating PDI, AOM and DI are first calculated automatically. Then you will be prompted to input the indices of the atoms in the ring that you are interested in, the input order must be consistent with atom connectivity.

PDI currently is only available for closed-shell systems, although theoretically it may be possible to be extended to open-shell cases.

Note that for completely planar systems, since DI can be decomposed to α and π parts, PDI can also be separated as PDI-α and PDI-π to individually study α aromaticity and π aromaticity. In order to calculate PDI-α (PDI-π), before entering present module, you should first manually set occupation number of all MOs except for π (α) MOs to zero (or you can utilize option 22 in main function 100 to do this step, which will be much more convenient).

3.18.7 Aromatic fluctuation index (FLU) and FLU-π (6,7)

Aromatic fluctuation index (FLU) was proposed in J. Chem. Phys., 122, 014109 (2005), also see Chem. Rev., 105, 3911 (2005) for more discussion. Like PDI, FLU is an aromaticity index based on DI, but can be used to study rings with any number of atoms. The FLU index was constructed by following the HOMA philosophy (see Section 3.28.6), i.e. measuring divergences (DI differences for each single pair bonded) from aromatic molecules chosen as a reference. FLU is defined as below

\[\mathrm{FLU}=\frac{1}{n}\sum_{A-B}^{\mathrm{ring}}\left[\left(\frac{V(B)}{V(A)}\right)^\alpha\left(\frac{\delta(A,B)-\delta_{\mathrm{ref}}(A,B)}{\delta_{\mathrm{ref}}(A,B)}\right)\right]^2\]

where the summation runs over all adjacent pairs of atoms around the ring, n is equal to the number of atoms in the ring, δ_ref is the reference DI value, which is precalculated parameter. α is used to ensure the ratio of atomic valences is greater than one

\[\alpha=\left\{\begin{gathered}1\qquad V(B)>V(A)\\ -1\qquad V(B)\leq V(A)\end{gathered}\right.\]

The first factor in the formula of FLU penalizes those with highly localized electrons, while the second factor measures the relative divergence with respect to a typical aromatic system. Obviously, lower FLU corresponds to stronger aromaticity.

The dependence on reference value is one of main weaknesses of FLU. The default δ_ref in Multiwfn for C-C, C-N, B-N are 1.468, 1.566 and 1.260 respectively, they are obtained from calculation of benzene, pyridine and borazine respectively under HF/6-31G* (geometry is optimized at the same level. Becke's atomic space with modified CSD radii and with sharpness parameter k=3 is used to derive δ_ref). Users can modify or add δ_ref through option -4.

The original paper of FLU also defined FLU-π, which is based on DI-π and π-atomic valence

\[\mathrm{FLU}_\pi=\frac{1}{n}\sum_{A-B}^{\mathrm{ring}}\left[\left(\frac{V_\pi(B)}{V_\pi(A)}\right)^\alpha\left(\frac{\delta_\pi(A,B)-\delta_{\mathrm{avg}}}{\delta_{\mathrm{avg}}}\right)\right]^2\]

where δ_π is the average value of the DI-π for the bonded atomic pairs in the ring, and the other symbols denote the aforementioned quantities calculated using π-orbitals only. The advantage of FLU-π over FLU is that FLU-π does not rely on predefined reference DI value, while the disadvantage is that FLU-π can only be exactly calculated for planar molecules.

Akin to FLU, the lower the FLU-π, the stronger aromatic the ring. If FLU-π is equal to zero, that means DI-π is completely equalized in the ring. The reasonableness to measure aromaticity by FLU-π is that aromaticity for most aromatic molecules is almost purely contributed by π electrons, rather than σ electrons.

In fuzzy atomic space analysis module of Multiwfn, PDI, FLU and FLU-π are calculated in fuzzy atomic spaces. In J. Phys. Chem. A, 110, 5108 (2006), the authors showed that the correlation between the PDI, FLU and FLU-π calculated in fuzzy atomic space and the ones calculated in AIM atomic space is excellent.

In Multiwfn, before calculating FLU and FLU-π, AOM will be calculated automatically. If you are calculating FLU-π, you will be prompted to input the indices of π orbitals, you can find out their indices by checking isosurface of all orbitals by main function 0. Then DI or DI-π matrix will be generated. Next, you should input the indices of the atoms in the ring, the input order must be consistent with atom connectivity. Besides FLU or FLU-π value, the contributions from each bonded atomic pair are outputted too.

FLU and FLU-π are only available for closed-shell system in Multiwfn. It is not well known whether FLU and FLU-π are also applicable for open-shell systems.

3.18.8 Condensed linear response kernel (CLRK) (9)

Linear response kernel (LRK) is an important concept defined in DFT framework, which can be written as

\[\chi(\mathbf{r}_1,\mathbf{r}_2)=\left(\frac{\delta^2 E}{\delta v(\mathbf{r}_1)\delta v(\mathbf{r}_2)}\right)_N=\left(\frac{\delta\rho(\mathbf{r}_1)}{\delta v(\mathbf{r}_2)}\right)_N\]

This quantity reflects the impact of the perturbation of external potential at \(r_{2}\) on the electron density at \(r_{1}\), which may also be regarded as the magnitude coupling between electron at \(r_{1}\) and \(r_{2}\).

In Multiwfn, LRK is evaluated by an approximation form based on second-order perturbation theory (see Eq.3 of Phys. Chem. Chem. Phys., 14, 3960 (2012))

\[\chi(\mathbf{r}_1,\mathbf{r}_2)\approx 4\sum_{i\in\mathrm{occ}}\sum_{j\in\mathrm{vir}}\frac{\varphi_i^*(\mathbf{r}_1)\varphi_j(\mathbf{r}_1)\varphi_j^*(\mathbf{r}_2)\varphi_i(\mathbf{r}_2)}{\varepsilon_i-\varepsilon_j}\]

where φ is molecular orbital, ε stands for MO energy. Note that this approximation form is only applicable to HF/DFT closed-shell systems, therefore present function only works for HF/DFT closed-shell systems.

Condensed linear response kernel (CLRK) is calculated as

\[\chi_{A,B}=\int_A\int_B\chi(\mathbf{r}_1,\mathbf{r}_2)\mathrm{d}\mathbf{r}_1\mathrm{d}\mathbf{r}_2=4\sum_{i\in\mathrm{occ}}\sum_{j\in\mathrm{vir}}\frac{S_{ij}(A)S_{ji}(B)}{\varepsilon_i-\varepsilon_j}\]

where A and S(A) denote fuzzy atomic space and atomic overlap matrix for atom A, similar for atom B. In Phys. Chem. Chem. Phys., 15, 2882 (2013), it was shown that CLRK is useful for investigation of aromaticity and anti-aromaticity.

Present function is used to calculate CLRK between all atomic pairs in current system, and the result will be outputted as a matrix. Due to evaluation of LRK requires virtual MO information, in current version .mwfn/.fch/.molden/.gms file must be used as input file.

Note that CLRK can be decomposed to orbital contribution, e.g. for MO i

\[\chi_{A,B}^{(i)}=4\sum_{j\in\mathrm{vir}}\frac{S_{ij}(A)S_{ji}(B)}{\varepsilon_i-\varepsilon_j}\]

For instance, assume that you want to evaluate the contribution from MO 3,4,7, then before calculating CLRK, you should enter main function 6 and use option 26 to set occupation number of all MOs except for 3,4,7 to zero. (Note that the virtual MOs used to calculate LRK will automatically still be the original virtual MOs, rather than the ones after modification of MO occupation numbers.)

3.18.9 Para linear response index (PLR) (10)

The definition of para linear response index (PLR) has an analogy to PDI, the only difference is that DI is replaced by CLRK

\[PLR(A,B)=\frac{\chi_{1,4}+\chi_{2,5}+\chi_{3,6}}{3}\]

In Phys. Chem. Chem. Phys., 14, 3960 (2012), the authors argued that PLR is as useful as PDI in quantitatively measuring aromaticity, and it is found that the linear relationship between PLR and PDI is as high as \(R^{2}\)=0.96.

Present function is used to calculate PLR. Multiwfn will first calculate CLRK, and then you should input the indices of the atoms constituting the ring in question, e.g. 3,5,6,7,9,2. The input order must be consistent with atom connectivity. Then PLR will be immediately outputted on screen. PLR is only applicable to HF/DFT closed-shell systems, and currently .mwfn/.fch/.molden/.gms must be used as input file.

Note that for completely planar systems, PLR can be exactly separated as PLR-α and PLR-π to individually study α aromaticity and π aromaticity. In order to calculate PLR-α (PLR-π), before entering present module, you should first manually set occupation number of all MOs except for π (α) MOs to zero (or you can utilize option 22 in main function 100 to do this step, which will be much more convenient).

3.18.10 Multi-center delocalization index (11)

n-center multi-center DI is calculated as

\[\delta(A,B,C\cdots H)=2^{n-1}\sum_i\sum_j\sum_k\cdots\sum_qS_{ij}(A)S_{jk}(B)S_{kl}(C)\cdots S_{qi}(H)\]

where i, j, k... only cycle occupied orbitals. The normalized form of multi-center DI is defined as δ^1/n, and may be compared between rings with different numbers of members.

Currently this function is only available for single-determinant closed-shell wavefunctions, and supports up to 10 centers. Note that for relatively large size of systems, calculating multi-center DI for more than 6 centers may be quite time-consuming.

3.18.11 Information-theoretic aromaticity index (12)

In ACS Omega, 3, 18370 (2018) it is shown that arithmetic mean of some information-theoretic quantities of the atoms constituting a ring has good linear relationship with other widely accepted aromaticity indices, such as HOMA and aromatic stabilization energy (ASE). It is thus clear that the arithmetic mean may be used as index for measuring aromaticity, although this point needs to be further explored.

The information-theoretic aromaticity index, namely the above-mentioned arithmetic mean can be calculated via subfunction 12 of fuzzy analysis module. After entering this function, you should choose the way of defining atomic information-theoretic quantity, three choices are currently available:

Atomic Shannon entropy: \(s_{\mathrm{S}}(A)=\int-\rho(\mathbf{r})\ln\rho(\mathbf{r})w_A(\mathbf{r})\mathrm{d}\mathbf{r}\)

Atomic Fisher information: \(i_{\mathrm{F}}(A)=\int\left|\nabla\rho(\mathbf{r})\right|^2/\rho(\mathbf{r})w_A(\mathbf{r})\mathrm{d}\mathbf{r}\)

Atomic GBP entropy: \(s_{\mathrm{GBP}}(A)=\int(3/2)\rho(\mathbf{r})\{\lambda+\ln[t(\mathbf{r})/t_{\mathrm{TF}}(\mathbf{r})]\}w_A(\mathbf{r})\mathrm{d}\mathbf{r}\)

Essentially, the three quantities correspond to the integral of user-defined functions 50, 51 and 54 in fuzzy atomic space. In this function, you also need to input the index of the atoms in the ring. Once calculation of the selected quantity for all atoms in the ring is finished, the average will be shown, and it can be regarded as an aromaticity index.

Before using this function, you can first select the way of defining atomic space. In the original paper, Hirshfeld partition was employed, while the default partition method of the fuzzy analysis module is Becke.

Information needed by fuzzy analysis module: GTFs, atom coordinates

3.18.12 Atomic effective volume, free volume, polarizability and C₆ coefficient (13)

Theory

Atomic (effective) polarizability of zero frequency, \({\alpha}_{}^{\mathrm{eff}}(0)\), in a molecule system, is a quantity of great important; however, there is no unique way to estimate it and it is not experimentally observable in general. In the original paper of Tkatchenko-Scheffler (TS) dispersion correction method, Phys. Rev. Lett., 102, 073005 (2009), the authors suggested a simple way of estimating it by scaling known free-atom polarizability \({\alpha}_{}^{\mathrm{free}}(0)\), they assumed that atomic polarizability of an element is positively proportional to its atomic volume. Then in a review article Chem. Rev., 117, 4714 (2017), this method is explicitly expressed as

\[\alpha_A^{\mathrm{eff}}(0)=\alpha_A^{\mathrm{free}}(0)\frac{V_A^{\mathrm{eff}}}{V_A^{\mathrm{free}}}\]

where atomic effective volume \(V^{\mathrm{eff}}\) in a molecule and atomic free volume \(V^{\mathrm{free}}\) are expressed as

\[\begin{gathered}V_A^{\mathrm{eff}}=\int\rho(\mathbf{r})w_A(\mathbf{r})\left|\mathbf{r}-\mathbf{R}_A\right|^3\mathrm{d}\mathbf{r}\\ V_A^{\mathrm{free}}=\int\rho_A^{\mathrm{free}}(\mathbf{r})\left|\mathbf{r}-\mathbf{R}_A\right|^3\mathrm{d}\mathbf{r}\end{gathered}\]

in which \(w_{A}\) is atomic weighting function of atom A. \(R_{A}\) is nuclear position of atom A. ρ is molecular electron density, and \({\rho}_{A}^{\mathrm{free}}\) is electron density of atom A in its free state. Note that the calculation level used for evaluating ρ and ρ^free must be exactly the same.

Ideally, the sum of \({\alpha}_{}^{\mathrm{eff}}(0)\) of all atoms should be equal to static polarizability of the whole system. However, since the above method is not rigorous, this condition is obviously impossible to achieve. In addition, it is worth noting that the choice of atomic weighting function significantly affects calculation result, but it is unclear which weighting function is the best choice for this purpose.

In my opinion, it is useful to define percentage contribution of an atom to total polarizability as

\[\alpha_A^{\%}=\frac{\alpha_A^{\mathrm{eff}}(0)}{\sum_A\alpha_A^{\mathrm{eff}}(0)}\times 100\%\]

This quantity is evidently useful in analyzing the major source of molecular polarizability.

According to the TS method, \(C_{6}\) dispersion coefficient of an atom in practical chemical environment can be easily and approximately evaluated as

\[C_{6,AA}^{\mathrm{TS}}=\left(\frac{V_A^{\mathrm{eff}}}{V_A^{\mathrm{free}}}\right)^2 C_{6,AA}^{\mathrm{free}}\]

where \({C}_{6,AA}^{\mathrm{free}}\) is the known \(C_{6}\) dispersion coefficient of the atom in free state. Furthermore, \(C_{6}\) between different two atoms can be calculated as

\[C_{6,AB}^{\mathrm{TS}}=\frac{2C_{6,AA}^{\mathrm{TS}}C_{6,BB}^{\mathrm{TS}}}{\frac{\alpha_B^{\mathrm{eff}}(0)}{\alpha_A^{\mathrm{eff}}(0)}C_{6,AA}^{\mathrm{TS}}+\frac{\alpha_A^{\mathrm{eff}}(0)}{\alpha_B^{\mathrm{eff}}(0)}C_{6,BB}^{\mathrm{TS}}}\]

Finally, intermolecular \(C_{6}\) coefficient can be obtained as follows

\[C_6^{\mathrm{mol}}=\sum_{A\in\mathrm{mol1}}\sum_{B\in\mathrm{mol2}}C_{6,AB}^{\mathrm{TS}}\]

Usage

Multiwfn is able to calculate \(V^{\mathrm{eff}}\), \(V^{\mathrm{free}}\), \({\alpha}_{}^{\mathrm{eff}}(0)\), \({\alpha}_{}^{\%}\), and \({C}_{6,AA}^{\mathrm{TS}}\) for all atoms. In addition, the \({C}_{6}^{\mathrm{mol}}\) between two present systems, in other words, homomolecular \({C}_{6}^{\mathrm{mol}}\), is also given. To evaluate them, the steps are:

1. Manually generate wavefunction file of each kind of atom in the current system by your favourite quantum chemistry code.

2. Boot up Multiwfn and load the wavefunction file of the system.

3. Enter main function 15, choose the weighting function you want to use by option -1.

4. Choose option 13. Multiwfn will ask you to input path of wavefunction file of each kind of element involved in the present system, the corresponding density will be used to evaluate \(V^{\mathrm{free}}\). After inputting the paths, calculation will be started.

During calculation, you can see \(V^{\mathrm{eff}}\), \(V^{\mathrm{free}}\) and their ratio \(V^{\mathrm{eff}}\)/\(V^{\mathrm{free}}\) are printed for each atom. Once the calculation is completely finished, Multiwfn will print \({\alpha}^{\mathrm{eff}}(0)\), \({\alpha}^{\%}\) and \({C}_{6,AA}^{\mathrm{TS}}\) for all atoms. Note that the \({\alpha}^{\mathrm{free}}(0)\) used for their evaluations come from the recommended values in CTCP atomic polarizability table (http://ctcp.massey.ac.nz/index.php?menu=dipole&page=dipole, data of Nov 16, 2020 version was taken). The built-in \({C}_{6,AA}^{\mathrm{free}}\) comes from J. Chem. Phys., 121, 4083 (2004), the value of H comes from Table I of original paper of TS method. Only the \({C}_{6,AA}^{\mathrm{free}}\) values of elements of the first five rows except for Zr-Cd are available, other elements have zero value.

If you only need aforementioned data for certain atoms in the system, you can use option -5 to define an atom list first, then other atoms will not be involved in the calculation of option 13.

Note

It is worth noting that if you choose to construct atomic weighting functions by Hirshfeld or Hirshfeld-I method in the way that manually providing atomic wavefunction files, these atomic wavefunction files are completely irrelevant to those used to evaluate \(V^{\mathrm{free}}\). The former must correspond to sphericalized density, while this requirement is dispensable for the latter.

An example of using this function is given in Section 4.15.4.