3.25 Conceptual density functional theory (CDFT) analysis (22)

The conceptual density functional theory (CDFT) originally developed by Robert Parr and extended by numerous researchers is a theory framework aiming for unraveling reactivity of chemical systems. CDFT contains numerous concepts and quantities, some of them can be used to predict favorable reactive sites and reactive character, and some of them can compare reactivity among different chemical species. Due to the high popularity and important role of CDFT in quantum chemistry, as well as there are so many meaningful relevant quantities, I believe it is very useful to develop a module to calculate all commonly employed quantities involved in CDFT with minimal steps.

In Part 1 of this section, I will briefly describe the definition of all quantities that can be studied via this module; then in Part 2, I will show how to use this module. This module is able to calculate the so-called orbital-weighted quantities, which will be specifically described in Part 3.

Note that aside from CDFT, Multiwfn also supports many other methods for revealing reactive sites, see Section 4.A.4 for overview.

If the functions described in this section are used in your research, please NOT ONLY cite original papers of Multiwfn, BUT ALSO cite the following book chapter, which comprehensively introduces feature and implementation of this module:

Tian Lu, Qinxue Chen, Realization of Conceptual Density Functional Theory and Information-Theoretic Approach in Multiwfn Program. In Conceptual Density Functional Theory, WILEY-VCH GmbH: Weinheim (2022); pp 631-647. DOI: 10.1002/9783527829941.ch31

3.25.1 Theory

To yield all below quantities, electronic energy (\(E\)) and electron density of \(N\), \(N+1\) and \(N-1\) electronic states must be available. Commonly \(N\) refers to the number of electrons carried by a chemical system at its most stable status. Geometry optimized for N-electrons state is employed for all calculations.

• Global indices

  • First vertical ionization potential (I₁): \(E(N-1)\) − \(E(N)\)

  • First vertical electron affinity (\(A\)): \(E(N)\) − \(E(N+1)\)

  • Mulliken electronegativity (χ): (I₁+\(A\))/2

  • Chemical potential (µ): −χ

  • Hardness (η): I₁−\(A\), which is also equivalent to fundamental gap. See J. Am. Chem. Soc., 105, 7512 (1983). Note that according to the convention employed by many CDFT papers, the prefix of 1/2 in original definition of η is dropped.

  • Softness (\(S\)): 1/η. See Proc. Nati. Acad. Sci., 82, 6723 (1985)

  • Electrophilicity index (ω): µ²/(2η). See J. Am. Chem. Soc., 121, 1922 (1999)

  • Nucleophilicity index (\(N_{\mathrm{Nu}}\)): \(E_{\mathrm{HOMO}}\)(Nu) − \(E_{\mathrm{HOMO}}\)(TCE), where Nu denotes nucleophile, TCE denotes tetracyanoethylene, whose HOMO energy is almost the lowest one among all organic molecules and therefore it is chosen as reference system. See J. Org. Chem., 73, 4615 (2008).

• Real space functions

  • Fukui function \(f(\mathbf{r})\) and dual descriptor Δ\(f(\mathbf{r})\): See Section 4.5.4 for detailed introduction

  • Local softness: \(s^+(\mathbf{r})\) = \(Sf^+(\mathbf{r})\), \(s^-(\mathbf{r})\) = \(Sf^-(\mathbf{r})\), \(s^0(\mathbf{r})\) = \(Sf^0(\mathbf{r})\) for nucleophilic, electrophilic, radical attacks, respectively, where \(f(\mathbf{r})\) is Fukui function of corresponding type. See Proc. Nati. Acad. Sci., 82, 6723 (1985)

  • Local hyper-softness: \(s^{(2)}\) ≈ \(S^2\)Δ\(f(\mathbf{r})\), see J. Math. Chem., 62, 461 (2024)

  • Local electrophilicity index: \({\omega}_{}^{\mathrm{loc}}(\mathbf{r})=\omega{f}^{+}(\mathbf{r})\)

  • Local nucleophilicity index: \({N}_{\mathrm{Nu}}^{\mathrm{loc}}(\mathbf{r})={N}_{\mathrm{Nu}}^{}{f}^{-}(\mathbf{r})\)

• Atom indices

  • Condensed Fukui function (\(f\)_A) and dual descriptor (Δ\(f\)_A): See Section 4.7.3 for detailed introduction

  • Condensed local softness

    • For nucleophilic attack: \({s}_{A}^{+}=S{f}_{A}^{+}\)

    • For electrophilic attack: \({s}_{A}^{-}=S{f}_{A}^{-}\)

    • For radical attack attack: \({s}_{A}^{0}=S{f}_{A}^{0}\)

  • Relative electrophilicity index: \({s}_{A}^{+}/{s}_{A}^{-}\), see J. Phys. Chem. A, 102, 3746 (1998)

  • Relative nucleophilicity index: \({s}_{A}^{-}/{s}_{A}^{+}\), see J. Phys. Chem. A, 102, 3746 (1998)

  • Condensed local electrophilicity index: \({\omega}^{A}=\omega{f}_{A}^{+}\)

  • Condensed local nucleophilicity index: \({N}_{\mathrm{Nu}}^{A}={N}_{\mathrm{Nu}}^{}{f}_{A}^{-}\)

  • Condensed local hyper-softness: \({s}_{A}^{(2)}≈{S}^{2}∆{f}_{A}\)

• Bond dual descriptor (BDD)

  BDD is an index defined for each bond, negative (positive) BDD implies the bond is nucleophilic (electrophilic), and the more negative (positive), the stronger the nucleophilicity (electrophilicity)

  BDD essentially is the first derivative of bond order (\(P\)) with respect to number of electrons under constant external potential. According to J. Math. Chem., 63, 1588 (2025), using finite-difference approximation, BDD for bond \(A\)-\(B\) can be calculated as \({BDD}_{A,B}={f}_{A,B}^{+}-{f}_{A,B}^{-}\), where \({f}_{A,B}^{+}\) and \({f}_{A,B}^{-}\) are two types of bond Fukui function and evaluated as follows

  \[{f}_{A,B}^{+}={P}_{A,B}\left( N+1 \right)-{P}_{A,B}\left( N \right)\]
  \[{f}_{A,B}^{-}={P}_{A,B}\left( N \right)-{P}_{A,B}(N-1)\]

  In J. Math. Chem., 63, 1588 (2025) the author employed Wiberg bond order based on atomic natural orbitals as \(P\), but in the current implementation of Multiwfn, fuzzy bond order based on Hirshfeld partition (see Section 3.11.6 for details) is used as \(P\) for evaluating BDD, because this bond order is insensitive to basis set (much better than Mayer bond order in this regards), works reasonably (according to my test), and relatively easy to realize.

• ω_cubic electrophilicity index

  The electrophilicity index ω_cubic introduced in J. Phys. Chem. A, 124, 2090 (2020) is somewhat special, it also relies on \(N-2\) electronic states. It includes higher-order terms than the aforementioned electrophilicity index ω. Its definition is

  \[\omega_{\mathrm{cubic}}=\omega\left(1+\frac{\mu}{3\eta^2}\gamma\right)\]

  In practice, it is calculated as

  \[\omega_{\mathrm{cubic}}=\frac{(\mu_{\mathrm{cubic}})^2}{2\eta_{\mathrm{cubic}}}\left[1+\frac{\mu_{\mathrm{cubic}}}{3(\eta_{\mathrm{cubic}})^2}\gamma_{\mathrm{cubic}}\right]\]

  where

  \[\begin{gathered}\mu_{\mathrm{cubic}}=(1/6)(-2A-5I_1+I_2)\\\eta_{\mathrm{cubic}}=I_1-A\\\gamma_{\mathrm{cubic}}=2I_1-I_2-A\end{gathered}\]

  where I₂ is the second vertical ionization potential and defined as \(E(N-2)\) − \(E(N-1)\). Correspondingly, there is a cubic form of condensed local electrophilicity index \({\omega}_{\mathrm{cubic}}^{A}={\omega}_{\mathrm{cubic}}{f}_{A}^{+}\). In J. Phys. Chem. A, 124, 2090 (2020) is shown that \({\omega}_{\mathrm{cubic}}^{A}\) value of halogen atom (which behaves as Lewis acid due to its σ-hole) in halogen-bond dimers R-X···NH₃ has excellent correlation with calculated binding energies (however, note that they employed AIM partition for atomic spaces rather than the Hirshfeld partition utilized in the present module).

• Electrophilic descriptor (ε)

  The electrophilic descriptor (ε) was introduced in Int. J. Quantum Chem., 124, e27366 (2024), it was shown to correlate with Mayr’s electrophilic parameter (\(E\)) significantly better than the electrophilicity index (ω) using a test set consisting of 35 organic molecules. In contrast to ω, whose derivation is only based on second-order Taylor expansion of system energy, derivation of ε is based on third-order expansion, which makes ε explicitly involve hyperhardness. Like ω_cubic, calculation of ε also relies on \(N-2\) electronic state.

  The working equation for computing ε is as follows

  \[\varepsilon=\chi\left(\frac{\phi}{\gamma}\right)-\left(\frac{\phi}{\gamma}\right)^2\left(\frac{\eta}{2}+\frac{\phi}{6}\right)\]

  with

  \[\phi=\sqrt{\eta^2-2\gamma\mu}-\eta\]

  It is important to note that the η, γ, χ, µ involved in above equations should be calculated in a different way than those described earlier, namely

  \[\begin{gathered}\mu=a\\\chi=-a\\\eta=2(b-ac)\\\gamma=-3c(b-ac)\end{gathered}\]

  where

  \[\begin{gathered}c=\frac{I_2-2I_1+A'}{2I_2-I_1-A'}\\b=\frac{I_1-A'}{2}-\frac{I_1+A'}{2}c\\a=-\frac{I_1+A'}{2}+\frac{I_1-A'}{2}c\end{gathered}\]

  with \(A\)’ = \(E(N+1)\) − \(E(N)\), which is electron affinity but differs with standard definition by the sign.

• Fukui potential and dual descriptor potential

  Definition and physical meaning of Fukui potential were carefully discussed in J. Phys. Chem. A, 115, 2325 (2011) and Int. J. Quantum Chem., 101, 520 (2005). Fukui potential is complementary to ESP in understanding energy variation in the early stage of chemical reaction when electron transfer is nonnegligible. For example, when an electrophile is attacked by a nucleophile, the external potential generated by the nucleophile felt by the electrophile is

  \[{v}_{N\mathrm{-phile}}\left( \mathbf{r} \right)≈{V}_{N-\mathrm{phile}}^{\mathrm{ESP}}\left( \mathbf{r} \right)-∆N\int \frac{{f}_{N-\mathrm{phile}}^{-}(\mathbf{r'})}{|\mathbf{r}-\mathbf{r}'|}\mathrm{d}\mathbf{r'}\]

  where Δ\(N\) is number of transferred electrons from the electrophile to nucleophile (negative value in this case), \({f}_{N-\mathrm{phile}}^{-}\) is Fukui function \(f^-\) of the nucleophile. Three kinds of Fukui potential are defined as follows

  \[{V}_{{f}^{-/+/0}}(\mathbf{r})=\int \frac{{f}^{-/+/0}(\mathbf{r'})}{|\mathbf{r}-\mathbf{r}'|}\mathrm{d}\mathbf{r}'\]

  ESP is only able to predict regioselectivity dominated by electrostatics effect with assumption that there is no electron transfer, however, evidently this assumption is far from true for general chemical reactions (but basically true for many noncovalent interactions). Clearly, Fukui potential should be more focused on in studying general reactions, especially for those with significant electron transfer.

  Fukui potential is more rigorous than Fukui function in predicting regioselectivity, as emphasized in J. Phys. Chem. A, 115, 2325 (2011) that “It is the value of the Fukui potential, more than the value of the Fukui function itself, that determines the reactive site”. However, Fukui potential is not so popular as Fukui function, because their distribution characteristics usually coincide with each other, while evaluation of Fukui potential needs calculating ESP twice, which is considerably more expensive than evaluating electron density twice. According to finite difference definition of Fukui function, for example, \({V}_{{f}^{-}}\) can be evaluated as

  \[{V}_{{f}^{-}}(\mathbf{r})=\int \frac{{f}^{-}(\mathbf{r'})}{|\mathbf{r}-\mathbf{r}'|}\mathrm{d}\mathbf{r'}=\int \frac{{\rho}_{N}\left( \mathbf{r'} \right)-{\rho}_{N-1}(\mathbf{r'})}{|\mathbf{r}-\mathbf{r}'|}\mathrm{d}\mathbf{r'=}{V}_{N-1}^{\mathrm{ESP}}(\mathbf{r})-{V}_{N}^{\mathrm{ESP}}(\mathbf{r})\]

  Furthermore, dual descriptor potential (DDP) was introduced in J. Math. Chem., 62, 1094 (2024), which is defined as

  \[DDP\left( \mathbf{r} \right)=\int \frac{∆f\left( {\mathbf{r}}^{'} \right)}{\left| \mathbf{r}-{\mathbf{r}}^{'} \right|}\mathrm{d}{\mathbf{r}}^{'}=\int \frac{{f}^{+}\left( {\mathbf{r}}^{'} \right)-{f}^{-}\left( {\mathbf{r}}^{'} \right)}{\left| \mathbf{r}-{\mathbf{r}}^{'} \right|}\mathrm{d}{\mathbf{r}}^{'}={V}_{{f}^{+}}\left( \mathbf{r} \right)-{V}_{{f}^{-}}(\mathbf{r})\]

  Unlike dual descriptor, which often has many nodal planes hindering discussion, distribution of DDP is much smoother, enabling researchers to identify preferential reactive sites easier.

  Like the Fukui function, the more positive the \({V}_{{f}^{+}}\) is in a region, the more likely it is for nucleophilic reaction to occur, and the region where the \({V}_{{f}^{-}}\) is more positive, the more likely it is for electrophilic reaction to occur. Similar to dual descriptor, the more positive and negative the DPP is in a region, the more likely this region is susceptible to undergo nucleophilic and electrophilic attacks, respectively.

3.25.2 Usage

For using the present module, namely main function 22, the file loaded after booting up Multiwfn is relatively arbitrary, the only requirement is that the atomic information in this file is identical to the system under study.

After entering the present module, you will see a menu, in which options 2, 3 and 9 are used to calculate above quantities. Before using them, generally you should provide N.wfn, N-1.wfn and N+1.wfn in current folder, which contain wavefunction and electronic energy of \(N\), \(N-1\) and \(N+1\) states respectively for present system, the geometries must be the same and correspond to the optimized geometry of \(N\) state. The calculation level of the three files must also be the same. If any of the three .wfn files is missing, Multiwfn will ask you to manually input path of .wfn file for corresponding state (.wfx, .fch and .mwfn files are also allowed, since they also carry wavefunction and electronic energy information).

Option 2: Used to calculate all aforementioned global indices and atomic indices, the result will be exported to CDFT.txt in current folder. Because as mentioned in Section 4.7.3, Hirshfeld method is an ideal choice for calculating condensed Fukui functions and may be other relevant atomic indices, therefore Hirshfeld charges are automatically calculated and used for evaluation of all atomic indices. Nucleophilicity index as well as its local version are dependent of HOMO energy of TCE, which should be calculated using the same level for present system, notice that these indices printed in present module simply employ the \(E_{\mathrm{HOMO}}\)(TCE) = -0.335198 Hartree calculated at the commonly used B3LYP/6-31G* level (clearly, if you want to get more reliable result and your current calculation level is not B3LYP/6-31G*, you should calculate \(E_{\mathrm{HOMO}}\)(TCE) yourself and then manually evaluate these indices).

Option 3: Used to calculate grid data of Fukui function, dual descriptor and functions related to them, then their isosurface maps can be directly visualized, the grid data can be exported to cube files in current folder. In this option you can set the scale factor to be multiplied to the calculated grid data. For example, if you set the scale factor to the global softness outputted by option 2, then the scaled \(f^-(\mathbf{r})\) Fukui function will correspond to \(s^-(\mathbf{r})\).

Option 9: Similar to option 3, but used to calculate grid data of Fukui potential and dual descriptor potential. Because grid data of ESP needs to be calculated for each charged state, which is often expensive, you should choose a proper grid setting to avoid being too time consuming.

Option 10: Used to calculate bond dual descriptor. Multiwfn will load .wfn file of different electronic states and calculate fuzzy bond orders in turn, and finally evaluate and print BDD values. The cost is generally low.

Generation of .wfn files

You can manually prepare the .wfn files used by options 2 and 3, alternatively, you can use option 1 of present module to automatically realize the preparation work.

After choosing option 1, you will be prompted to input Gaussian keywords for single point task, then Multiwfn asks you to input charge and spin multiplicity for \(N\), \(N+1\) and \(N-1\) states in turn, then Gaussian single point input files N.gjf, N+1.gjf and N-1.gjf will be generated in current folder (the geometry in these files correspond to the geometry in the input file of Multiwfn). Then, you can manually use Gaussian to run them to obtain N.wfn, N+1.wfn and N-1.wfn, or if Gaussian has been installed on your computer, you can directly let Multiwfn to invoke Gaussian to calculate them (in this case the "gaupath" parameter in settings.ini must have been set to actual path of Gaussian executable file), after calculations the three .wfn files will appear in current folder.

Sometimes we need to use mixed basis set, in this case you should prepare a file named basis.txt in current folder, which records the definition of basis set (may be also accompanied by pseudopotential definition). If the inputted keyword contains "gen" or "genecp", then the content of basis.txt will be automatically appended to the end of the generated .gjf files.

Multiwfn is also able to generate input files of ORCA for producing the three .wfn files. You should choose option -2 and select ORCA, then choosing option 1 will generate N.inp, N-1.inp and N+1.inp. If you have set “orcapath” in settings.ini to actual path of ORCA executable file, you can directly let Multiwfn to invoke ORCA to run them to yield N.wfn, N-1.wfn and N+1.wfn; alternatively, you can run them by ORCA manually, and then put the resulting N.wfn, N-1.wfn and N+1.wfn in current folder.

To use present module to study large organic systems, commonly I suggest using B3LYP/6-31G* level, because this level is inexpensive, while the quality of the yielded quantities is already satisfactory.

Note on calculating ω_cubic and ε

By default, Multiwfn does not calculate ω_cubic and ε, because they rely on \(N-2\) electronic state. If you need them, you should first select option -1 to switch the status to "Yes". Then you can use option 1 to help you to prepare .wfn file for \(N\), \(N+1\), \(N-1\), \(N-2\) electronic states, or you manually provide them. Then after selecting option 2, the resulting CDFT.txt file will contain condensed local ω_cubic, global ω_cubic, ε, as well as I₂.

An example of using present module to calculate various CDFT quantities for phenol is provided in Section 4.22.1. Example of calculating Fukui potential and dual descriptor potential for maleic anhydride is given in Section 4.22.4. Examples of calculating bond dual descriptor are given in Section 4.22.5.

3.25.3 Special topic 1: Orbital-weighted Fukui function and dual descriptor

Theory

The originally defined Fukui function and dual descriptor do not work well when frontier molecular orbitals are (quasi-)degenerate. For example, when HOMO and HOMO-1 have very similar or exactly identical energies, the Fukui function \(f^-\) may be unable to give meaningful result or the result is fully misleading; in addition, when the system shows point group symmetry, such as C₆₀ fullerene, the distribution of \(f^-\) is usually not in consistency with molecular symmetry, this is an apparently unexpected observation.

In order to address these problems, in J. Comput. Chem., 38, 481 (2017), the authors proposed orbital-weighted Fukui function, and in J. Phys. Chem. A, 123, 10556 (2019), they further proposed orbital-weighted dual descriptor, they are summarized below (the \({f}_{w}^{0}\) is defined by me)

\[\begin{gathered}f_w^+(\mathbf{r})=\sum_{i=\mathrm{LUMO}}^{\infty}w_i\left|\varphi_i(\mathbf{r})\right|^2\qquad w_i=\frac{\exp\left[-\left(\frac{\mu-\varepsilon_i}{\Delta}\right)^2\right]}{\sum_{i=\mathrm{LUMO}}^{\infty}\exp\left[-\left(\frac{\mu-\varepsilon_i}{\Delta}\right)^2\right]}\\f_w^-(\mathbf{r})=\sum_{i}^{\mathrm{HOMO}}w_i\left|\varphi_i(\mathbf{r})\right|^2\qquad w_i=\frac{\exp\left[-\left(\frac{\mu-\varepsilon_i}{\Delta}\right)^2\right]}{\sum_{i}^{\mathrm{HOMO}}\exp\left[-\left(\frac{\mu-\varepsilon_i}{\Delta}\right)^2\right]}\\f_w^0(\mathbf{r})=\left[f_w^+(\mathbf{r})+f_w^-(\mathbf{r})\right]/2\\\Delta f_w(\mathbf{r})=f_w^+(\mathbf{r})-f_w^-(\mathbf{r})\end{gathered}\]

where ε_i and \(\varphi\)_i are energy and wavefunction of orbital \(i\); µ is chemical potential and approximately calculated as (\(E_{\mathrm{HOMO}}\)+\(E_{\mathrm{LUMO}}\))/2 in the above formulae. The Δ is an adjustable parameter, in principle its best value is the one able to make the functions have ideal predictability of local reactivity. Clearly the most appropriate Δ is dependent on practical system, usually 0.1 Hartree is a worth-trying guess. If you find the orbital-weighted functions under this value do not work well, you can try to properly change it and redo calculations.

Compared to the frozen orbital approximation form of \(f^-\), namely \(f^-(\mathbf{r})\)=|\(\varphi\)_HOMO(r)|², the advantage of \({f}_{w}^{-}\) is that it takes all orbitals into account with different weights. From the expression it can be seen that the closer a low-lying orbital energy is to the HOMO energy, the greater its weight. Evidently degenerate orbitals share the same weight. The Gaussian function involved in the formula behaves as a decay function, the larger the Δ, the higher the contribution of low-lying orbitals to the \({f}_{w}^{-}\). When energy difference between HOMO and HOMO-1 is significant, there will be no reason to use \({f}_{w}^{-}\) instead of \(f^-\). The situation is similar for \({f}_{w}^{+}\), \({f}_{w}^{0}\) and \(∆{f}_{w}\).

Usage

Since orbital-weighted functions involve virtual orbitals, you should use .mwfn, .fch, .molden or .gms file as input file. Commonly, the geometry in the input file should correspond to the optimized geometry of \(N\)-electronic state; however, it is also possible to study them for nonequilibrium structure, such as a point in intrinsic reaction coordinate (IRC).

Only closed-shell single-determinant wavefunction is acceptable. Diffuse functions should not be used if you intend to calculate \({f}_{w}^{+}\), \({f}_{w}^{0}\) and \(∆{f}_{w}\), since they utilize virtural orbitals, whose chemical meaning may be severely broken when diffuse functions are employed.

In main function 22, four options are related to the orbital-weighted calculation:

• Option 4: Set the Δ parameter used in the subsequent orbital-weighted calculations

• Option 5: Print the highest 10 weights (i.e. the {\(w\)} in the aforementioned formulae) involved in the orbital-weighted calculations. This option is useful to check if current Δ parameter is reasonable and helps users to better understand how the orbital-weighted method works

• Option 6: Calculating condensed \({f}_{w}^{+}\), \({f}_{w}^{-}\), \({f}_{w}^{0}\) and \(\Delta{f}_{w}^{}\) values, in other words, calculating integration of these functions in Hirshfeld atomic spaces. The result is useful in quantitatively examining net amount of these functions at various atoms. The default radial and angular integration points are usually fine enough, if you find the sum of condensed \({f}_{w}^{+}\) or \({f}_{w}^{-}\) deviates from 1.0 evidently, you should set "iautointgrid" parameter in settings.ini to 0 and then properly enlarge "radpot" and "sphpot" parameters.

• Option 7: Calculating grid data of \({f}_{w}^{+}\), \({f}_{w}^{-}\), \({f}_{w}^{0}\) and \(\Delta{f}_{w}^{}\) functions, then you can directly visualize their isosurfaces or export them as cube files so that you can render them via third-part softwares such as VMD and ChimeraX.

Examples of using this module to calculate orbital-weighted Fukui function and orbital-weighted dual descriptor are provided in Section 4.22.2.

3.25.4 Special topic 2: (Quasi-)degenerate Fukui function and dual descriptor based on electron density

3.25.4.1 Closed-shell case

Theory

In the last section, I have introduced orbital-weighted Fukui function and dual descriptor, which are suitable when frontier molecular orbitals are (quasi-)degenerate. However, they are defined based on orbital approximation, namely the orbital relaxation effect is not taken into account, while this effect cannot be always safely overlooked. In J. Comput. Chem., 37, 2279 (2016), an alternative form of Fukui function and dual descriptor that work for (quasi-)degenerate HOMO/LUMO case was proposed, and this form is defined directly based on electron density, that means orbital relaxation effect is fully taken into account as the original Fukui function and dual descriptor. This (quasi-)degenerate Fukui function and dual descriptor based on electron density will be referred to as \(f_{\mathrm{Q}}\) and Δ\(f_{\mathrm{Q}}\), respectively.

The idea of \(f_{\mathrm{Q}}\) is very simple. If at electronic state \(N\) the degree of degeneracy of LUMO and HOMO is \(p\) and \(q\), respectively, then three forms of \(f_{\mathrm{Q}}\) are evaluated as follows

\[\begin{gathered}f_{\mathrm{Q}}^+(\mathbf{r})=\frac{\rho_{N+p}(\mathbf{r})-\rho_N(\mathbf{r})}{p}\\f_{\mathrm{Q}}^-(\mathbf{r})=\frac{\rho_N(\mathbf{r})-\rho_{N-q}(\mathbf{r})}{q}\\f_{\mathrm{Q}}^0(\mathbf{r})=\left[f_{\mathrm{Q}}^+(\mathbf{r})+f_{\mathrm{Q}}^-(\mathbf{r})\right]/2\end{gathered}\]

Δ\(f_{\mathrm{Q}}\) can be evaluated based on \({f}_{\mathrm{Q}}^{+}\) and \({f}_{\mathrm{Q}}^{-}\) as usual

\[\Delta f_{\mathrm{Q}}(\mathbf{r})=f_{\mathrm{Q}}^+(\mathbf{r})-f_{\mathrm{Q}}^-(\mathbf{r})\]

Clearly, if both HOMO and LUMO are nondegenerate, then \(f_{\mathrm{Q}}\) and Δ\(f_{\mathrm{Q}}\) will be equivalent to the original form of Fukui function and dual descriptor, \(f\) and Δ\(f\), respectively.

To reasonably calculate \(f_{\mathrm{Q}}\) and Δ\(f_{\mathrm{Q}}\), it is crucial to properly determine \(p\) and \(q\). Commonly they can be assigned by inspecting energies of several lowest unoccupied MOs and several highest occupied MOs, respectively. If energy difference between an occupied (unoccupied) MO and HOMO (LUMO) is very small, e.g. less than 0.01 eV, then they may be regarded as degenerate. Obvious, there is no strict energy threshold for judging orbital degeneracy, and in some cases you may need to judge by considering various factors, e.g. HOMO-LUMO gap, reasonableness of actual calculation result, orbital shape, etc.

Note that \(f_{\mathrm{Q}}\) and Δ\(f_{\mathrm{Q}}\) are defined only for closed-shell case.

Multiwfn is not only able to calculate \(f_{\mathrm{Q}}\) and Δ\(f_{\mathrm{Q}}\), but also able to calculate local properties described in Section 3.25.1 (except for ω_cubic) based on them. For example, local softness with consideration of degeneracy is product of global softness and \(f_{\mathrm{Q}}\). Note that the involved first VEA and VIP are still evaluated as usual, namely VEA = \(E(N)\) − \(E\)(\(N\)+1) and VIP = \(E\)(\(N\)-1) − \(E(N)\).

Spin multiplicity must be properly chosen for calculating wavefunction file of \(N\)+\(p\) and \(N\)-\(q\) states, usually they should be set to \(p+1\) and \(q+1\), see J. Comput. Chem., 37, 2279 (2016) for detailed discussion. This setting commonly is able to guarantee that the attached (detached) electrons equally enter (leave from) all degenerate LUMOs (HOMOs).

Usage

In main function 22, both \(p\) and \(q\) are defaulted to be 1, namely degeneracy of frontier molecular orbitals (FMOs) is not taken into account. To manually set \(p\) and \(q\) and thus consider the (quasi-)degeneracy effect in the subsequent calculations, you should choose option “-3 Set degree of FMO degeneracy” first. Then information of 10 lowest unoccupied MOs and that of 10 highest occupied MOs will be listed on screen (in the case that the input file contains wavefunction information), you should properly input \(p\) and \(q\) according to the listed MO energies.

After setting \(p\) and \(q\), you can use option 1 to use Multiwfn to help you generate input files of single point task of Gaussian or ORCA code for \(N\), \(N\)+\(p\) and \(N\)-\(q\) states, you will be asked to input net charge and spin multiplicity for these states. In addition, if either \(p\) or \(q\) is not equal to 1, then Multiwfn will also ask you if also generating input file for \(N+1\) and/or \(N-1\) states, because \(E(N+1)\) and \(E(N-1)\) are needed to evaluate first VIP, VEA and related quantities in option 2. After running these input files, .wfn files of the states will be yielded. Of course, you can also manually generate wavefunction files for the states via any of your favourite quantum chemistry codes.

Finally, you can use corresponding options to calculate various quantities or functions that defined in the CDFT framework, the output content is exactly identical to the default nondegenerate case, though \(p\) and \(q\) have been properly considered in the current situation. Note that if \(p\) or \(q\) is not equal to 1, and at the same time either N+1.wfn or N-1.wfn is not available in current folder, then quantities related to first VIP and VEA will not be calculated and outputted by option 2.

An example is given in Section 4.22.3.

3.25.4.2 Open-shell case

Theory

For open-shell cases, Martínez Araya proposed working equations of evaluating dual descriptor in Chem. Phys. Lett., 506, 104 (2011) based on the formalism of spin-polarized conceptual DFT. The equations were then further generalized (private communication) and summarized below, where \(∆{N}_{S}\) corresponds to change in spin number (\({N}_{S}={N}_{\alpha}-{N}_{\beta}\)), B denotes external magnetic field; \({p}_{\alpha}\) and \({p}_{\beta}\) are degeneracy of LUMO(α) and LUMO(β), respectively; \({q}_{\alpha}\) and \({q}_{\beta}\) are degeneracy of HOMO(α) and HOMO(β), respectively.

Different forms to evaluate \({f}^{+}={\left( \frac{∂\rho\left( \mathbf{r} \right)}{∂N} \right)}_{v\left( \mathbf{r} \right),\mathbf{B}\left( \mathbf{r} \right)}^{+}\), namely spin multiplicity increases, decreases, and approximately unchanged:

\[ \begin{gathered} {f}_{\Delta{N}_{S}<0}^{+}\left( \mathbf{r} \right)=\frac{{\rho}_{N+{p}_{\beta}}\left( \mathbf{r} \right)-{\rho}_{N}\left( \mathbf{r} \right)}{{p}_{\beta}}\\ {f}_{\Delta{N}_{S}>0}^{+}\left( \mathbf{r} \right)=\frac{{\rho}_{N+{p}_{\alpha}}\left( \mathbf{r} \right)-{\rho}_{N}\left( \mathbf{r} \right)}{{p}_{\alpha}}\\ {f}_{\Delta{N}_{S}≈0}^{+}\left( \mathbf{r} \right)=\frac{{f}_{\Delta{N}_{S}>0}^{+}\left( \mathbf{r} \right)+{f}_{\Delta{N}_{S}<0}^{+}\left( \mathbf{r} \right)}{2} \end{gathered} \]

Different forms to evaluate \({f}^{-}={\left( \frac{∂\rho\left( \mathbf{r} \right)}{∂N} \right)}_{v\left( \mathbf{r} \right),\mathbf{B}\left( \mathbf{r} \right)}^{-}\):

\[ \begin{gathered} {f}_{\Delta{N}_{S}<0}^{-}\left( \mathbf{r} \right)=\frac{{\rho}_{N}\left( \mathbf{r} \right)-{\rho}_{N-{q}_{\alpha}}\left( \mathbf{r} \right)}{{q}_{\alpha}}\\ {f}_{\Delta{N}_{S}>0}^{-}\left( \mathbf{r} \right)=\frac{{\rho}_{N}\left( \mathbf{r} \right)-{\rho}_{N-{q}_{\beta}}\left( \mathbf{r} \right)}{{q}_{\beta}}\\ {f}_{\Delta{N}_{S}≈0}^{-}\left( \mathbf{r} \right)=\frac{{f}_{\Delta{N}_{S}>0}^{-}\left( \mathbf{r} \right)+{f}_{\Delta{N}_{S}<0}^{-}\left( \mathbf{r} \right)}{2} \end{gathered} \]

Different forms of evaluating dual descriptor:

\[ \begin{gathered} \Delta{f}_{\Delta{N}_{S}<0}\left( \mathbf{r} \right)={f}_{\Delta{N}_{S}<0}^{+}\left( \mathbf{r} \right)-{f}_{\Delta{N}_{S}<0}^{-}\left( \mathbf{r} \right)\\ \Delta{f}_{\Delta{N}_{S}>0}\left( \mathbf{r} \right)={f}_{\Delta{N}_{S}>0}^{+}\left( \mathbf{r} \right)-{f}_{\Delta{N}_{S}>0}^{-}\left( \mathbf{r} \right)\\ \Delta{f}_{\Delta{N}_{S}≈0}\left( \mathbf{r} \right)={f}_{\Delta{N}_{S}≈0}^{+}\left( \mathbf{r} \right)-{f}_{\Delta{N}_{S}≈0}^{-}\left( \mathbf{r} \right) \end{gathered} \]

In principle, \(\Delta{f}_{\Delta{N}_{S}≈0}\) should be the most ideal choice, because it follows the restriction of keeping the spin-number constant. However, there are two ways to also keep the spin-number constant, which only involve change in α and β electrons, respectively, but they are not examined in literature. \(\Delta{f}_{\alpha}\left( \mathbf{r} \right)={f}_{\Delta{N}_{S}>0}^{+}\left( \mathbf{r} \right)-{f}_{\Delta{N}_{S}<0}^{-}\left( \mathbf{r} \right)\) \(\Delta{f}_{\beta}\left( \mathbf{r} \right)={f}_{\Delta{N}_{S}<0}^{+}\left( \mathbf{r} \right)-{f}_{\Delta{N}_{S}>0}^{-}\left( \mathbf{r} \right)\)

Usage

Subfunction 12 of the CDFT module is dedicated to evaluating all aforementioned functions. After booting up Multiwfn, usually you should:

• (1) Load a wavefunction file containing basis function information and all orbitals

• (2) Enter subfunction 12 of main function 22

• (3) If this system has degeneracy in frontier MOs, choose option 1 to set up degeneracy

• (4) Choose option 2 to generate input files of single point task of Gaussian for all involved electronic states. Then you can ask Multiwfn to invoke Gaussian to directly calculate them, or manually calculate them. Then put the resulting .wfn files (with the same names as the .gjf files) into the current folder.

• (5) Choose option 3 to generate grid data of electron density of different electronic states. In the post-processing menu, you can directly visualize isosurface map of Fukui function or dual descriptor of selected form, or export corresponding grid data as .cub file.

An example is given in Section 4.22.3.3.

3.25.5 Special topic 3: Nucleophilic and electrophilic superdelocalizabilities

Theory

Nucleophilic and electrophilic delocalizabilities are also known as nucleophilic and electrophilic superdelocalizabilities, they were proposed by Schüürmann in Environ. Toxicof. Chem., 9, 417 (1990) and Quant. Struct.-Act. Relat., 9, 326 (1990), and have been employed as molecular descriptors for building quantitative structure-activity relationship (QSAR) equations. In the Schüürmann’s work, nucleophilic superdelocalizability (\(D^{\mathrm{N}}\)) and electrophilic superdelocalizability (\(D^{\mathrm{E}}\)) of atom A are defined as follows, respectively

\[\begin{gathered}D^{\mathrm{N}}(A)=2\sum_{i}^{\mathrm{unocc}}\sum_{\mu\in A}\frac{C_{\mu,i}^2}{\alpha-\varepsilon_i}\\D^{\mathrm{E}}(A)=2\sum_{i}^{\mathrm{occ}}\sum_{\mu\in A}\frac{C_{\mu,i}^2}{\varepsilon_i-\alpha}\end{gathered}\]

where ε_i is energy of molecular orbital \(i\), α = (\(E_{\mathrm{HOMO}}\) + \(E_{\mathrm{LUMO}}\))/2, µ∈\(A\) stands for basis function µ of atom A, and C is coefficient matrix. Evidently, both \(D^{\mathrm{N}}\) and \(D^{\mathrm{E}}\) are negative.

However, the expression of superdelocalizabilities by Schüürmann is only suitable for semi-empirical calculation, which employs orthonormal basis functions. In Sci. Rep., 5, 13695 (2015), a different version of electrophilic superdelocalizability was proposed and it is compatible with non-orthonormal basis functions. In this work, it was shown that electrophilic superdelocalizability of an atom is closely related to its atomic polarizability.

The above definitions of superdelocalizability are based on molecular orbital expansion coefficients; in contrast, in Multiwfn, superdelocalizabilities are evaluated based on Hirshfeld partition of atomic spaces, this form is more robust and fully compatible with diffuse functions. Specifically, in Multiwfn, the nucleophilic and electrophilic superdelocalizabilities are calculated as follows

\[\begin{gathered}D^{\mathrm{N}}(A)=2\sum_{i}^{\mathrm{unocc}}\frac{\Theta_{A,i}}{\alpha-\varepsilon_i}\\D^{\mathrm{E}}(A)=2\sum_{i}^{\mathrm{occ}}\frac{\Theta_{A,i}}{\varepsilon_i-\alpha}\end{gathered}\]

where Θ_A,i is composition of atom A in orbital \(i\) calculated by Hirshfeld method, see Section 3.10.5 for detail. Multiwfn also calculates the superdelocalizabilities without the α shift parameter, namely

\[\begin{gathered}D^{\mathrm{N0}}(A)=2\sum_{i}^{\mathrm{unocc}}\frac{\Theta_{A,i}}{-\varepsilon_i}\\D^{\mathrm{E0}}(A)=2\sum_{i}^{\mathrm{occ}}\frac{\Theta_{A,i}}{\varepsilon_i}\end{gathered}\]

Usage

Since superdelocalizabilities involve virtual orbitals, you should use .mwfn, .fch, .molden or .gms file as input file. Only closed-shell single-determinant wavefunction is acceptable. Diffuse functions should not be used if you intend to calculate \(D^{\mathrm{N}}\) and \(D^{\mathrm{N0}}\), since it utilizes virtual orbitals, whose chemical meaning may be severely broken when diffuse functions are employed.

After loading input file, enter main function 22, then choose option 8, you will obtain \(D^{\mathrm{N}}\), \(D^{\mathrm{E}}\), \(D^{\mathrm{N0}}\), and \(D^{\mathrm{E0}}\) for all atoms. Example of output (examples\oxirane.fchk):

    Atom             D_N                D_E              D_N_0             D_E_0
    1(C )        -26.56182        -10.85916        -35.92863         -8.75591
    2(C )        -26.56182        -10.85916        -35.92863         -8.75591
    3(O )        -20.09086        -19.14131        -25.62598        -14.71454
    4(H )        -10.19410         -2.66387        -14.52037         -2.13024
    5(H )        -10.19410         -2.66387        -14.52037         -2.13024
    6(H )        -10.19410         -2.66387        -14.52037         -2.13024
    7(H )        -10.19410         -2.66387        -14.52037         -2.13024

Sum of D_N:        -113.99092 /Hartree
Sum of D_E:         -51.51513 /Hartree
Sum of D_N_0:      -155.56471 /Hartree
Sum of D_E_0:       -40.74734 /Hartree

3.25.6 Special topic 4: Dual delocalization descriptor

Theory

The theory of dual delocalization descriptor (DDD) was proposed in Phys. Chem. Chem. Phys., 28, 19133 (2026) by Samir Kenouche. This method quantifies responses of the interatomic delocalization index (DI) to variations in the total electron number. This theory is closely related to the bond dual descriptor (BDD) introduced in Section 3.25.1, essentially both of them exhibit how degree of electron sharing changes in the course of electron attachment or removal. However, BDD is calculated based on finite difference of bond orders between two electronic states (\(N\) and \(N-1\), or \(N\) and \(N+1\)), while DDD has an analytical form and its calculation only relies on wavefunction of \(N\) state, therefore it has advantage in computational cost. Another advantage is that the components of DDD enable one to examine the role played by frontier molecular orbitals on the response, see original paper for discussions. Note that DDD theory was derived based on frozen orbital approximation (FOA), and only single-determinant closed-shell wavefunctions are supported.

Here I describe key ideas and ingredients of DDD. This theory was derived under the Ángyán-Loos-Mayer (ALM) formulation of DI (δ):

\[{\delta}_{A,B}^{\mathrm{ALM}}=\sum_{i,j} {\eta}_{i}{\eta}_{j}{S}_{ij}(A){S}_{ij}(B)\]

where η is orbital occupation number, S is atomic overlap matrix (AOM), so \({S}_{ij}(A)\) corresponds to integral of product of MO \(i\) and MO \(j\) in the space of atom A. \(i\) and \(j\) loop over all spatial orbitals.

Derivatives of DI have discontinuity at integer number of electrons. With FOA and assuming that electrons can only leave from HOMO, it is shown that left-sided 1st derivative of DI can be expressed as follows, which reflects HOMO-driven response of DI with respect to decrease in \(N\)

\[{\left( \frac{∂{\delta}_{A,B}}{∂N} \right)}^{-}=4\sum_{i<\mathrm{H}} {S}_{\mathrm{H}i}(A){S}_{\mathrm{H}i}(B)+4{S}_{\mathrm{HH}}(A){S}_{\mathrm{HH}}(B)\]

in which the second term represents HOMO self-contribution part, while the first term is referred to as \({f}_{AB}^{-}\), namely \({f}_{AB}^{-}=4\sum_{i<\mathrm{H}} {S}_{\mathrm{H}i}(A){S}_{\mathrm{H}i}(B)\).

Right-sided 1st derivative of DI is (\(i\) loops over all doubly occupied orbitals)

\[{\left( \frac{∂{\delta}_{A,B}}{∂N} \right)}^{+}={f}_{AB}^{+}=4\sum_{i∈\mathrm{occ}} {S}_{\mathrm{L}i}(A){S}_{\mathrm{L}i}(B)\]

First-order dual delocalization descriptor is defined as \({f}_{AB}^{(1)}={f}_{AB}^{+}-{f}_{AB}^{-}\), which characterizes the asymmetry in the response of DI to variations in the total electron number, >0 and <0 clearly indicate LUMO-dominated response and HOMO-dominated response, respectively.

\(g\) terms are related to quadratic contributions to the DI response from HOMO or LUMO:

\[{g}_{AB}^{+}={S}_{\mathrm{LL}}(A){S}_{\mathrm{LL}}(B)\]

\[{g}_{AB}^{-}={2S}_{\mathrm{HH}}(A){S}_{\mathrm{HH}}(B)\]

\[{\left( \frac{{∂}^{2}{\delta}_{A,B}}{∂{N}^{2}} \right)}^{+}=2{g}_{AB}^{+}\]

\[{\left( \frac{{∂}^{2}{\delta}_{A,B}}{∂{N}^{2}} \right)}^{-}={g}_{AB}^{-}\]

When the total number of electrons changes by an integer, variation of DI can be expressed as

\[{∆{\delta}_{AB}^{+}=f}_{AB}^{+}+{g}_{AB}^{+}\]

\[{∆{\delta}_{AB}^{-}=f}_{AB}^{-}+\frac{3}{2}{g}_{AB}^{-}\]

It should be emphasized that \(∆{\delta}_{AB}^{+}\) and \(∆{\delta}_{AB}^{-}\) are just approximation to accurate variations of DI, which can be respectively evaluated using finite-difference as \({\delta}^{+}=\delta\left( N+1 \right)-\delta(N)\) and \({\delta}^{-}=\delta\left( N \right)-\delta(N-1)\).

Second-order dual delocalization descriptor is defined as \({f}_{AB}^{(2)}=∆{\delta}_{AB}^{+}-∆{\delta}_{AB}^{-}\), which corresponds to the asymmetry of the response of DI resulting from the addition and removal of an electron.

When HOMO and/or LUMO are degenerated, the aforementioned terms are calculated as follows to take the degeneracy into account:

\[{f}_{AB}^{+}=\frac{1}{{n}_{\mathrm{L}}}\sum_{l∈\mathrm{L}} {f}_{AB}^{+(l)} {f}_{AB}^{-}=\frac{1}{{n}_{\mathrm{H}}}\sum_{h∈\mathrm{H}} {f}_{AB}^{-(h)}\]

\[{g}_{AB}^{+}=\frac{1}{{n}_{\mathrm{L}}}\sum_{l∈\mathrm{L}} {g}_{AB}^{+(l)} {g}_{AB}^{-}=\frac{1}{{n}_{\mathrm{H}}}\sum_{h∈\mathrm{H}} {g}_{AB}^{-(h)}\]

where \(n_{\mathrm{L}}\) and \(n_{\mathrm{H}}\) are degeneracy of LUMO and HOMO, \(l\) and \(h\) loop over all degenerated LUMOs and HOMOs, respectively. The terms with (\(l\)) or (\(h\)) superscript refer to those calculated with only viewing \(l\) as LUMO and \(h\) and HOMO, respectively. Note that when calculating \({f}_{AB}^{-(h)}\), orbital index only loops over the occupied orbitals that do not belong to the degenerated HOMOs.

Usage

Multiwfn is able to calculate \({f}_{AB}^{+}\), \({f}_{AB}^{-}\), \({f}_{AB}^{(1)}\), \({g}_{AB}^{+}\), \({g}_{AB}^{-}\), \({f}_{AB}^{(2)}\) for any bond in the system. The system should be closed-shell, and a wavefunction file containing basis function information (e.g. .fch, .molden) should be used as input when Multiwfn boots up.

This function corresponds to subfunction 11 of main function 22. After entering it, usually you should choose option 1 or option 2 to calculate AOM using AIM partition or Hirshfeld partition of atomic spaces, respectively, then you will be asked to input the bond to be studied, then results will be immediately shown on screen. It is worth noting that if Hirshfeld partition is used to generate AOM, the resulting DI also corresponds to fuzzy bond order (Section 3.11.6), and hence response of fuzzy bond order to \(N\) is studied.

If you already have a .txt file containing AOM of all MOs of the present system exported by fuzzy analysis module or basin analysis module, in this function you can also directly choose option 3 to load and use it, so that AOM will not be calculated here. Note that before generating AOM using fuzzy or basin analysis module, you should set “ispecial” in settings.ini to 3, otherwise only occupied orbitals will be involved in the AOM.

An example is given in Section 4.22.6.