2.7 User-defined real space function¶
In real space function selection menu, you can find a term named "User-defined real space function". In order to avoid lengthy list of real space functions, numerous uncommonly used real space functions are not explicitly presented in the list. However, if you want to use them, you can set "iuserfunc" parameter in settings.ini to one of the indices (see below), then the user-defined function will be pinned to corresponding function. For example, before running Multiwfn, if you set "iuserfunc" to 2, then the user-defined real space function will be equivalent to density of beta electrons. An alternative way of setting user-defined function is inputting iu in the main menu, you can input the index of the user-defined function.
In fact, the user-defined function corresponds to “userfunc” function in source file function.f90. By filling in proper code yourself, the functions supported by Multiwfn can be easily extended. For example, after filling the code "userfunc=fgrad(x,y,z,'t')**2/8/fdens(x,y,z)" into proper place of “function userfunc” in function.f90 and recompile Multiwfn, the integrand of Weizsäcker kinetic energy functional, namely \(\tau_\mathrm{W}[\rho]=\int\left|\nabla\rho(\mathbf{r})\right|^2/[8\rho(\mathbf{r})]\mathrm{d}\mathbf{r}\), will be ready for use. When you write your own code you can refer to existing codes, a list of built-in functions is given in Appendix 2 of this manual.
Prebuilt user-defined functions¶
-2 Promolecular density calculated based on built-in sphericalized atomic densities. You can check Appendix 3 on how these atomic densities are produced. Options -3 and -4 of main function 6 could be used to exclude contribution of some atoms to promolecular density.
Note
Please note that this promolecular density is different to the promolecular density used in real space functions 14 and 16, because in which the atomic densities from H to Ar are directly taken from the original paper of NCI method rather than calculated according to the method described in Appendix 3.
-1 Value evaluated by trilinear interpolation from grid data. The grid data can be generated by main function 5, or loaded from .cub/.grd/.vti file when Multiwfn boot up. This function is quite useful. For example, by making use of this function you can plot the grid data as curve map and plane map via main functions 3 and 4, respectively; you can also gain atomic contribution to the total value via subfunction 1 of fuzzy atomic space analysis module (main function 15).
-3 The same as -1, but use 3D cubic B-spline interpolation instead. The cost is higher than trilinear interpolation but smoother and usually more accurate especially when grid spacing is relatively large.
Note
However, the function interpolated via cubic B-spline shows unwanted feature where the variation of actual function cannot be well approximated by polynomial and in this case the result is even worse than trilinear interpolation. For example, cusp character of electron density at nucleus cannot be faithfully represented. In addition, the function shows fluctuating feature when the actual function shows a steep slope close to a region of small slope, or when the grid points inadequately represent the actual function.
Note
About consideration of periodicity for user-defined functions -1 and -3: Note that if the grid data represent a real space function in the entire crystal cell, and you want to consider periodic boundary condition during the interpolation, you should select option 18 in main function 1000 (a hidden function) to set box information of the grid data as cell information before performing calculation. Alternatively, add a string box2cell into the first line of the .cub file to be loaded to inform Multiwfn automatically doing this.
In the case that periodicity is not taken into account, the position to be evaluated should be fully within the spatial range of the grid data, otherwise the interpolation cannot be carried out and thus the returned function value will simply be zero.
0 This function corresponds to a constant value of 1.0
1 Alpha density: \(\rho^\alpha(\mathbf{r})=[\rho(\mathbf{r})+\rho^s(\mathbf{r})]/2\)
2 Beta density: \(\rho^\beta(\mathbf{r})=[\rho(\mathbf{r})-\rho^s(\mathbf{r})]/2\)
3 Integrand of electronic spatial extent <\(r^{2}\)>: \((x^2+y^2+z^2)\rho(\mathbf{r})\)
4 Weizsäcker potential (closed-shell form): \(V_\mathrm{W}(\mathbf{r})=\frac{1}{8}\frac{\left|\nabla\rho(\mathbf{r})\right|^2}{\rho^2(\mathbf{r})}-\frac{1}{4}\frac{\nabla^2\rho(\mathbf{r})}{\rho(\mathbf{r})}\)
5 Integrand of Weizsäcker functional (closed-shell form): \(\tau_\mathrm{W}(\mathbf{r})=\left|\nabla\rho(\mathbf{r})\right|^2/[8\rho(\mathbf{r})]\), which is the exact kinetic energy density of any one-orbital system (one or two electrons, or any number of bosons). If you need spin polarized form, use user-defined function 1200.
6 Radial distribution function of electron density: \(4\pi\times\rho(\mathbf{r})\times(x^2+y^2+z^2)\). Clearly, spherical symmetry of electron density is assumed.
7 Local temperature, in Hartree/\(k_\mathrm{B}\) (PNAS, 81, 8028): \(T(\mathbf{r})=[2G(\mathbf{r})]/[3\rho(\mathbf{r})]\). Note that if ρ is less than the “uservar” parameter in settings.ini, then T will be treated as 0.
8 Average local electrostatic potential (J. Chem. Phys., 72, 3027 (1980)): \(V_{\mathrm{ESP}}(\mathbf{r})/\rho(\mathbf{r})\)
9 Shape function: ρ(r)/N, where N is the total number of electrons
10 Potential energy density (Virial field): \(V(\mathbf{r})=-K(\mathbf{r})-G(\mathbf{r})=(1/4)\nabla^2\rho(\mathbf{r})-2G(\mathbf{r})\)
11 Electron energy density: \(E(\mathbf{r})=G(\mathbf{r})+V(\mathbf{r})=-K(\mathbf{r})\) (Sometimes E is written as H). Note that integral of E(r) over the whole space is not exactly identical to the electronic energy printed by quantum chemistry code, because in practice the virial ratio is always more or less violated.
-11 Scaled electron energy density: \(E_{\mathrm{scl}}(\mathbf{r})=-K(\mathbf{r})\times(R-1)\), where R is the virial ratio loaded from input file (beware that not all kinds of input file contain this information! fch/mwfn/wfx/wfn format have specific field to record it). Integral of \(E_{\mathrm{scl}}\)(r) over the whole space is exactly identical to the electronic energy printed by quantum chemistry code. Section 4.17.9 illustrates the usefulness of this function.
12 Local nuclear attraction potential energy: \(-\rho(\mathbf{r})\times V_{\mathrm{nuc}}(\mathbf{r})\)
13 Kinetic energy density per electron: \(G(\mathbf{r})/\rho(\mathbf{r})\) This quantity at bond critical point is able to discriminate covalent bonding and closed-shell interaction, see J. Am. Soc. Chem., 120, 13429 (1998)
14 Electrostatic potential from electrons: \(V_{\mathrm{ele}}(\mathbf{r})=V_{\mathrm{ESP}}(\mathbf{r})-V_{\mathrm{nuc}}(\mathbf{r})=-\int\frac{\rho(\mathbf{r}')}{\left|\mathbf{r}-\mathbf{r}'\right|}\mathrm{d}\mathbf{r}'\). The negative of this function is also known as Hartree potential, which represents classical Coulomb potential on an electron at r caused by all electrons.
15 Bond metallicity: \(\xi_\mathrm{J}(\mathbf{r})=\rho(\mathbf{r})/\nabla^2\rho(\mathbf{r})\) At bond critical point, \(\xi_\mathrm{J}\)>1 indicates metallic interaction, see J. Phys.: Condens. Matter, 14, 10251 (2002).
16 Dimensionless bond metallicity: \(\xi_\mathrm{m}(\mathbf{r})=\frac{36(3\pi^2)^{2/3}}{5}\frac{\rho(\mathbf{r})^{5/3}}{\nabla^2\rho(\mathbf{r})}\) At bond critical point, larger value corresponds to stronger metallicity of the bond, see Chem. Phys. Lett., 471, 174 (2009).
17 Energy density per electron: E(r)/ρ(r), this value at BCP is called as bond degree parameter (BD), see J. Chem. Phys., 117, 5529 (2002) for detail, in which the authors advocate that for covalent bonds (viz. \(E_{\mathrm{BCP}}\)<0), the BD renders covalence degree, and the stronger the interaction the greater the BD magnitude; while for closed-shell interactions (viz. \(E_{\mathrm{BCP}}\)>0), the BD can be viewed as softening degree, and the weaker the interaction the larger the BD magnitude.
18 Region of Slow Electrons (RoSE), which is defined in Chem. Phys. Lett., 582, 144 (2013) and has a pattern very similar to ELF, with value space of [-1,1]: \(\nu_\pm=\frac{D_0(\mathbf{r})-G(\mathbf{r})}{D_0(\mathbf{r})+G(\mathbf{r})}\)
19 Single exponential decay detector (SEDD), which is highly analogous to ELF. Its updated definition in J. Chem. Theory Comput., 10, 3745 (2014) is implemented: \(\mathrm{SEDD}(\mathbf{r})=\ln\left\{1+\left[\frac{\nabla(\nabla\rho(\mathbf{r})/\rho(\mathbf{r}))^2}{\rho(\mathbf{r})}\right]^2\right\}\)
20 Density overlap regions indicator (DORI), defined in J. Chem. Theory Comput., 10, 3745 (2014): \(\mathrm{DORI}(\mathbf{r})=\theta(\mathbf{r})/[1+\theta(\mathbf{r})]\), where \(\theta(\mathbf{r})=[\nabla(\nabla\rho(\mathbf{r})/\rho(\mathbf{r}))^2]^2/[\nabla\rho(\mathbf{r})/\rho(\mathbf{r})]^6\). DORI is mainly used to reveal interatomic interaction regions, see Section 3.23.3 for more information. IRI works much better than DORI.
21 Integrand of X component of electric dipole moment: −x×ρ(r)
22 Integrand of Y component of electric dipole moment: −y×ρ(r)
23 Integrand of Z component of electric dipole moment: −z×ρ(r)
24 Approximate form of DFT linear response kernel for closed-shell (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}\)
25 Magnitude of fluctuation of electronic momentum: \(\tilde{P}(\mathbf{r})=\left|\nabla\rho(\mathbf{r})\right|/[2\rho(\mathbf{r})]\), also known as localized electrons detector (LED), is useful to discuss bonding and very similar to reduced density gradient, see Theor. Chem. Acc., 127, 393 (2010)
26 Integrand of Thomas-Fermi kinetic energy functional (closed-shell form): τTF(r) = \(C_{\mathrm{TF}}\)ρ(r)5/3, where \(C_{\mathrm{TF}}\)=(3/10)(3π2)2/3=2.871234. This is the exact kinetic energy density of non-interacting, uniform electron gas. If you need spin-polarized form, use user-defined function 1200.
27 Local electron affinity (LEA): \(EA_\mathrm{L}(\mathbf{r})=\frac{-\sum_{i\in\mathrm{vir}}\left|\varphi_i(\mathbf{r})\right|^2\varepsilon_i}{\sum_{i\in\mathrm{vir}}\left|\varphi_i(\mathbf{r})\right|^2}\), which is very similar to average local ionization energy, but i cycles all unoccupied orbitals. See J. Mol. Model., 9, 342 (2003). Illustration of applying this function for a practical molecule is given in Section 4.12.13.
-27 Local electron attachment energy: \(E_{\mathrm{att}}(\mathbf{r})=\frac{n\sum_{i=\mathrm{LUMO}}^{\varepsilon_i<0}\left|\varphi_i(\mathbf{r})\right|^2\varepsilon_i}{\rho(\mathbf{r})}\). i loops over all unoccupied orbitals with negative energy. For restricted and unrestricted wavefunctions, n equals 2 and 1, respectively. See J. Phys. Chem. A., 120, 10023 (2016). This function has similar usefulness to LEA but is more robust. Illustration of applying this function for a practical molecule is given in Section 4.12.13.
28 Local Mulliken electronegativity: \(\chi_\mathrm{L}(\mathbf{r})=[\bar{I}(\mathbf{r})+EA_\mathrm{L}(\mathbf{r})]/2\), see J. Mol. Model., 9, 342 (2003)
29 Local hardness: \(\eta_\mathrm{L}(\mathbf{r})=[\bar{I}(\mathbf{r})-EA_\mathrm{L}(\mathbf{r})]/2\), see J. Mol. Model., 9, 342 (2003)
Note: In order to use \(EA_\mathrm{L}\), \(E_{\mathrm{att}}\), χL and ηL, the input file must contain both occupied and unoccupied orbitals of a single determinant wavefunction (but restricted open-shell is not supported), formats such as .mwfn, .fch, .molden and .gms should be used as input file.
Generally, \(EA_\mathrm{L}\) (and thus χL and ηL) is not compatible with basis sets with diffuse functions, while any basis set may be used for evaluating \(E_{\mathrm{att}}\). The much better compatibility with augmented basis set is a notable advantage of \(E_{\mathrm{att}}\) over \(EA_\mathrm{L}\). However, in order to use \(E_{\mathrm{att}}\), at least LUMO should has a negative energy, while often this condition is not satisfied under commonly used levels. It is found in the original paper of \(E_{\mathrm{att}}\) that this function works reasonably with B3LYP/6-31+G(d,p) orbitals.
30 Ellipticity of electron density: \(\varepsilon(\mathbf{r})=[\lambda_1(\mathbf{r})/\lambda_2(\mathbf{r})]-1\), where λ1 and λ2 are the lowest and the second lowest eigenvalues of Hessian matrix of ρ, respectively. At bond critical point (BCP), λ1 and λ2 are both negative and exhibit the curvatures of electron density in the two orthogonal directions that perpendicular to the bond. The ε at BCP is often viewed as an indicator of asymmetric distribution of electron density around the bond, the higher deviation to axisymmetric distribution, the larger the ε value at BCP.
31 eta index: \(\eta(\mathbf{r})=\left|\lambda_1(\mathbf{r})\right|/\lambda_3(\mathbf{r})\), where λ1 and λ3 are the lowest and the highest eigenvalues of Hessian matrix of ρ, respectively. It was argued that the value of η at bond critical point is less than unity for closed shell interactions and increases with increasing covalent character, see Angew. Chem. Int. Ed., 53, 2766 (2014), as well as J. Phys. Chem. A, 114, 552 (2010) for discussions.
32 Modified eta index by Tian Lu: \(\eta'(\mathbf{r})=\left|\lambda_1(\mathbf{r})\right|/\lambda_3(\mathbf{r})-1\). Similar to η, but negative value of η' corresponds to closed shell interactions.
33 Potential acting on one electron in a molecule (PAEM), see J. Comput. Chem., 35, 965 (2014): \(V_{\mathrm{PAEM}}(\mathbf{r})=-\sum_A\frac{Z_A}{\left|\mathbf{r}-\mathbf{R}_A\right|}+\frac{1}{\rho(\mathbf{r})}\int\frac{\pi(\mathbf{r},\mathbf{r}')}{\left|\mathbf{r}-\mathbf{r}'\right|}\mathrm{d}\mathbf{r}'\), π is pair density. It can be further written as \(V_{\mathrm{PAEM}}(\mathbf{r})=-V_{\mathrm{ESP}}(\mathbf{r})+V_{\mathrm{XC}}(\mathbf{r})=-V_{\mathrm{ESP}}(\mathbf{r})+\frac{1}{\rho(\mathbf{r})}\int\frac{\Gamma_{\mathrm{XC}}^{\alpha,\mathrm{tot}}(\mathbf{r},\mathbf{r}')+\Gamma_{\mathrm{XC}}^{\beta,\mathrm{tot}}(\mathbf{r},\mathbf{r}')}{\left|\mathbf{r}-\mathbf{r}'\right|}\mathrm{d}\mathbf{r}'\), where \(V_{\mathrm{ESP}}\) is the total molecular electrostatic potential and \(V_{\mathrm{XC}}\) is exchange-correction potential. In Multiwfn, the exchange-correlation density Γ is evaluated in terms of Müller approximation, see part 17 of Section 2.6 for detail. It was shown in the original paper that PAEM may be useful to distinguish covalent and noncovalent interactions. An example of application of PAEM can be found in Section 4.3.3.
34 The same as 33, but now \(V_{\mathrm{XC}}\) directly corresponds to DFT exchange-correlation potential. Its specific form can be chosen via "iDFTxcsel" parameter, see the end of this section for detail. This form of PAEM is significantly cheaper than 33 but only supports closed-shell wavefunction. For very large systems, cost of 33 is computationally forbidden, so 34 is the only choice.
35 \(\left| V\left( \mathbf{r} \right) \right|/G(\mathbf{r})\). In J. Chem. Phys., 117, 5529 (2002) it was proposed that this quantity at BCP can be used to discriminate interaction types. <1 corresponds to closed-shell interaction; >2 corresponds to covalent interaction; while the case that >1 and < 2 corresponds to intermediate interaction.
36 On-top pair density, namely the two positions of the pair density are identical: π(r,r). See such as Int. J. Quantum Chem., 61, 197 (1995) for discussion. "paircorrtype" parameter in settings.ini mentioned earlier affects the result.
37 The strong covalent interaction index (SCI) defined in J. Phys. Chem. A, 122, 3087 (2018) and further examined in J. Mol. Model., 24, 213 (2018). This function was shown to be very useful for identifying very strong covalent bonds. SCI is expressed as SCI(r) = 1/ζ(r), where ζ(r) = [τS(r) − τW(r)]/τTF. The meaning of τS, τW and τTF are described in Section 2.6 where ELF is introduced. It is easy to find SCI index is closely related to ELF, which can be written as ELF(r) = 1/[1+ζ 2(r)].
38 The angle between the second eigenvector of Hessian of electron density and the vector perpendicular to a given plane, which can be defined by option 4 of main function 1000; the unit vector normal to the plane will be shown on screen, assume that you use three points A, B, C to define the plane and you get vector u, but what you really want is -u, you can then input the points again but in reverse sequence, i.e. C, B, A. In J. Phys. Chem. A, 115, 12512 (2011) this quantity along bond paths was used to reveal π interaction.
39 Electrostatic potential without contribution of a specific nucleus:
\(V_\mathrm{n}(\mathbf{r})=\sum_{A\neq K}\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 the nucleus K can be set by option 3 of main function 1000 (which is hidden in main interface but can be chosen). \(V_\mathrm{n}\) is a useful quantity, for example if K is chosen as index of a hydrogen, then the value correlates with its pKa, because in this case \(V_\mathrm{n}\) approximately reflects the binding energy of a proton at the position of K and rest of the system; In addition, J. Phys. Chem. A, 118, 1697 (2014) showed that \(V_\mathrm{n}\) can be used to quantitatively predict interaction energy of the weak interactions dominated by electrostatic effect (viz. H-bonds, halogen-bonds, dihydrogen bonds), see Section 4.1.2 for introduction and example.
Note
It is worth to mention that in main function 1, when you request Multiwfn to print properties at nuclear position of an atom, the electrostatic potential without contribution of nuclear charge of this atom is automatically printed.
40 Steric energy density: \(\left|\nabla\rho(\mathbf{r})\right|^2/[8\rho(\mathbf{r})]\), which is equivalent to integrand of Weizsäcker functional.
41 Steric potential: \(\upsilon_\mathrm{s}(\mathbf{r})=\frac{1}{8}\frac{\left|\nabla\rho(\mathbf{r})\right|^2}{[\rho(\mathbf{r})+\delta]^2}-\frac{1}{4}\frac{\nabla^2\rho(\mathbf{r})}{\rho(\mathbf{r})+\delta}\). The negative of this quantity is also known as one-electron potential (OEP). Notice that the δ is a very small term artificially introduced to avoid the denominator converges to zero faster than nominator. The value of δ can be determined by "steric_addminimal" in settings.ini. If δ is set to 0, then the original expression of steric potential is recovered.
42 Steric charge: \(q_\mathrm{s}(\mathbf{r})=\nabla^2\upsilon_\mathrm{s}(\mathbf{r})/(-4\pi)\). Evaluated by finite difference of steric force.
43 Magnitude of steric force: \(F_\mathrm{s}(\mathbf{r})=\left|-\nabla\upsilon_\mathrm{s}(\mathbf{r})\right|\). Evaluated analytically.
Notice that the δ term, which is mentioned above, also affects steric force and steric charge. Discussions of steric energy/potential/force/charge can be found in J. Chem. Phys., 126, 244103 (2007).
44,45,46 Damped Steric potential, steric force based on damped steric potential, directly damped Steric force: Documented privately
47 Damped steric charge: Documented privately
49 Relative Shannon entropy density, also known as information gain density: \({i}_{\mathrm{G}}=\rho(\mathbf{r}\mathrm{)ln}\frac{\rho(\mathbf{r})}{{\rho}_{0}(\mathbf{r})}\), where ρ0(r) is promolecular density. Before performing any analysis or visualization for this function, you must enter main function 1000 (a hidden function) and then select subfunction 17 to construct promolecular wavefunction and store it into a special space in memory; this promolecular wavefunction is generated by preparing wavefunction files for all atoms (see Section 3.7.3 for details) and then combining them together. Note that in the subfunction 17, when Multiwfn asks you “Do you want to make wavefunction information in memory correspond to the just generated promolecular wavefunction? (y/n)”, you should input n to avoid replacing the present wavefunction with the promolecular wavefunction.
50 Shannon entropy density: \(s_s(\mathbf{r})=-\rho(\mathbf{r})\ln\rho(\mathbf{r})\).
51 Fisher information density: \(i_\mathrm{F}(\mathbf{r})=\left|\nabla\rho(\mathbf{r})\right|^2/\rho(\mathbf{r})\). Note that this quantity is only different from Weizsäcker functional by a constant factor of 1/8.
52 Second Fisher information density: \(i'_\mathrm{F}(\mathbf{r})=-\nabla^2\rho(\mathbf{r})\ln\rho(\mathbf{r})\). The relationship between Shannon entropy and Fisher information can be found in J. Chem. Phys., 126, 191107 (2007). The integration of this function over the whole space is exactly identical to that of Fisher information density.
53 Ghosh entropy density or Ghosh-Berkowitz-Parr (GBP) entropy density, in \(k_\mathrm{B}\) (PNAS, 81, 8028 (1984)): \(s(\mathbf{r})=(3/2)\rho(\mathbf{r})\{\lambda+\ln[t(\mathbf{r})/t_{\mathrm{TF}}(\mathbf{r})]\}\), where \(\lambda=5/3+\ln(4\pi c_{\mathrm{k}}/3)\), \(c_{k}\) is Thomas-Fermi constant (3/10)(3π2)2/3=2.871234, Lagrangian kinetic energy density G(r) is chosen to be the kinetic energy density term t(r).
54 The same as 53, but G(r) − ∇2ρ(r)/8, which is the kinetic energy density exactly corresponding to Eq. 22 of PNAS, 81, 8028 (1984), is employed as the t(r). In rare cases this definition of kinetic energy density leads to a very small negative value, since at these points \(t_{\mathrm{TF}}\)(r) is also very close to zero, in order to normally get result, ln[t(r)/\(t_{\mathrm{TF}}\)(r)] in this case is simply set to zero.
55 Integrand of the integral part of quadratic form of Rényi entropy: ρ2(r). The integral of this function over the whole space corresponds to second-order Onicescu information.
56 Integrand of the integral part of cubic form of Rényi entropy: ρ3(r) . The integral of this function over the whole space and then divided by 2 corresponds to third-order Onicescu information.
57 \(g_1(\mathbf{r})=\nabla^2\rho(\mathbf{r})\ln\frac{\rho(\mathbf{r})}{\rho_0(\mathbf{r})}\)
58 \(g_2(\mathbf{r})=\rho(\mathbf{r})\left[\frac{\nabla^2\rho(\mathbf{r})}{\rho(\mathbf{r})}-\frac{\nabla^2\rho_0(\mathbf{r})}{\rho_0(\mathbf{r})}\right]\)
59 \(g_3(\mathbf{r})=\rho(\mathbf{r})\left[\nabla\ln\frac{\rho(\mathbf{r})}{\rho_0(\mathbf{r})}\right]^2\)
Note that user-defined functions 57, 58, 59 can only be studied by plotting as map in main functions 3, 4, 5, and meantime deformation and promolecular maps are not available.
60 Pauli potential: \(V_\mathrm{θ}\). Multiwfn supports two ways to evaluate this function, depending on "ispecial" in settings.ini:
-
ispecial=0: According to Eq. 17 of Comput. Theor. Chem., 1006, 92 (2013), \(V_\mathrm{θ}\) = µ + \(V_{\mathrm{ESP}}\) − \(V_{\mathrm{XC}}\) − \(V_\mathrm{W}\), where µ is chemical potential (assumed to be zero by Multiwfn), \(V_{\mathrm{ESP}}\) is electrostatic potential (as described in Section 2.6) and \(V_{\mathrm{XC}}\) is exchange-correlation potential (its form is determined by "iDFTxcsel" in settings.ini, see later).
-
ispecial=1: Use original definition, namely \({V}_\mathrm{θ}={V}_{\mathrm{S}}-{V}_{\mathrm{W}}=\delta{\tau}_{\mathrm{S}}\left[ \rho \right]/\delta\rho-\delta{\tau}_{\mathrm{W}}\left[ \rho \right]/\delta\rho\), where \(V_\mathrm{S}\) is potential of non-interacting kinetic energy functional (τS), and \(V_\mathrm{W}\) is potential of Weizsäcker kinetic energy functional (τW). The form of τS can be chosen by "iKEDsel" in settings.ini, see later; currently only iKEDsel of 3, 5, 7 are supported.
Since Pauli force and charge are calculated by finite difference based on Pauli potential, and quantum potential/force/charge are defined based on Pauli potential/force/charge, the user-defined functions 61~65 are also affected by "ispecial".
61 Pauli force magnitude: \(F_\mathrm{θ}(\mathbf{r})=\left|-\nabla V_\mathrm{θ}(\mathbf{r})\right|\)
62 Pauli charge: \(q_\mathrm{θ}(\mathbf{r})=\nabla^2 V_\mathrm{θ}(\mathbf{r})/(-4\pi)\)
63 Quantum potential: \(V_{\mathrm{q}}=V_\mathrm{θ}+V_{\mathrm{XC}}\). When ispecial=0, clearly in Multiwfn it simply corresponds to \(V_{\mathrm{ESP}}\) − \(V_\mathrm{W}\)
64 Quantum force magnitude: \(F_{\mathrm{q}}(\mathbf{r})=\left|-\nabla V_{\mathrm{q}}(\mathbf{r})\right|\)
65 Quantum charge: \(q_{\mathrm{q}}(\mathbf{r})=\nabla^2 V_{\mathrm{q}}(\mathbf{r})/(-4\pi)\)
66 Magnitude of electrostatic force: \(F_{\mathrm{ESP}}(\mathbf{r})=\left|-\nabla[-V_{\mathrm{ESP}}(\mathbf{r})]\right|\). As discussed in Phys. Chem. Chem. Phys., 19, 1496 (2017), there is a very close relationship between steric force, quantum force and electrostatic force, which are user-defined functions 43, 64 and 66, respectively.
67 Electrostatic charge: \(q_{\mathrm{ESP}}(\mathbf{r})=\nabla^2[-V_{\mathrm{ESP}}(\mathbf{r})]/(-4\pi)\)
68 Energy density of electronic part of electrostatic term of Shubin Liu's energy decomposition: \(\upsilon_e(\mathbf{r})=\rho(\mathbf{r})\left[\int\frac{\rho(\mathbf{r}')}{\left|\mathbf{r}-\mathbf{r}'\right|}\mathrm{d}\mathbf{r}'-\sum_A\frac{Z_A}{\left|\mathbf{r}-\mathbf{R}_A\right|}\right]=-V_{\mathrm{ESP}}(\mathbf{r})\rho(\mathbf{r})\). To learn more about Shubin Liu's energy decomposition, see Section 3.24.2.
69, -69 Energy density of quantum part of Shubin Liu's energy decomposition: τs(r) − τW(r) + εXC(r), where τs is Hamiltonian kinetic energy density (if userfunc = 69) or Lagrangian kinetic energy density (if userfunc = -69). τW is integrand of Weizsäcker functional (the same as userfunc=5), εXC is exchange-correlation energy density (the same as userfunc=1000), whose form can be chosen by "iDFTxcsel" parameter in settings.ini, see later for detail.
70 Phase-space-defined Fisher information density (PS-FID): \(i_{\mathrm{f}}(\mathbf{r})=\frac{3\rho(\mathbf{r})}{k_{\mathrm{B}}T(\mathbf{r})}=\frac{9}{2}\frac{\rho^2(\mathbf{r})}{G(\mathbf{r})}\), where T(r) is local temperature as shown above. This function has very similar characters to ELF and LOL, the spatial localization of electron pairs can be clearly revealed. See Chem. Phys., 435, 49 (2014) for introduction and illustrative applications.
71, 72, 73, 74 Electron linear momentum density (EMD) in 3D representation. The electron linear momentum operator is \(-i\nabla\), for an orbital the expectation value of linear momentum in X is \(-i\left<\varphi\left|\nabla\right|\varphi\right>\), so X component of EMD for a wavefunction can be defined as (the imaginary sign is ignored) \(p_x(\mathbf{r})=-\sum_i\eta_i\varphi_i^*(\mathbf{r})\partial\varphi_i(\mathbf{r})/\partial x\), similar for Y and Z components. The 71, 72, 73th user-defined functions correspond to X, Y, Z component, respectively. The magnitude (74) of EMD is defined as \(p_{\mathrm{tot}}(\mathbf{r})=\sqrt{p_x^2(\mathbf{r})+p_y^2(\mathbf{r})+p_z^2(\mathbf{r})}\).
75, 76, 77, 78 Magnetic dipole moment density (MDMD). The operator for magnetic dipole moment is the angular momentum operator (see Theor. Chim. Acta, 6, 341 (1966))
where i, j, k are unity vectors in X, Y and Z directions, respectively. Therefore, X, Y and Z component of MDMD can be defined as follows (the imaginary sign is ignored)
The 75, 76, 77th user-defined functions correspond to X, Y, Z components, respectively. The 78th user-defined function is magnitude of MDMD: \(m_{\mathrm{tot}}(\mathbf{r})=\sqrt{m_x^2(\mathbf{r})+m_y^2(\mathbf{r})+m_z^2(\mathbf{r})}\).
79 Gradient norm of electron energy density: |∇E(r)|
80 Laplacian of electron energy density: ∇2E(r)
81, 82, 83 X, Y, Z component of Hamiltonian kinetic energy density, respectively.
84, 85, 86 X, Y, Z component of Lagrangian kinetic energy density, respectively.
87, 88, 89 Local total, dynamic, nondynamic electron correlation functions, respectively. See introduction and calculation example of these functions in Section 4.A.7.2.
90 Fractional Occupation Number Weighted Electron Density (FOD) proposed by Grimme in Angew. Chem. Int. Ed., 54, 1 (2015). See introduction and calculation example of FOD in Section 4.A.7.1.
91 δ\(g^{\mathrm{inter}}\) between two fragments defined in independent gradient model based on Hirshfeld partition (IGMH) method, see Sections 3.23.5 and 3.23.6 for detail. To use this function, you should first enter main function 1000 and select option 16 to define two fragments.
92, 93, 94 Van der Waals potential, repulsion potential and dispersion potential based on UFF forcefield parameters, see Section 3.23.7 for detailed introduction. The unit is in kcal/mol, the probe atom can be set by "ivdwprobe" in settings.ini. Note that 92 is identical to the 25th function described in Section 2.6.
95, 96, 97, 98 Orbital-weighted f +, f −, f 0 Fukui functions as well as orbital-weighted dual descriptor, respectively. They were originally proposed in J. Comput. Chem., 38, 481 (2017) and J. Phys. Chem. A, 123, 10556 (2019). They are useful in studying local reactivity for systems whose frontier molecular orbitals are (quasi-)degenerate, see Section 3.25.3 for brief introduction and Section 4.22.2 for illustrative examples. The Δ parameter in these functions can be set by option 6 of main function 1000 (a hidden option). Since these functions involve virtual orbitals, you should use mwfn/fch/molden/gms as input file. Only closed-shell single-determinant wavefunction is supported.
99 Interaction region indicator (IRI): The same as function 24 described in Section 2.6.
100 Disequilibrium (also known as semi-similarity): \({D}_{\mathrm{r}}\left( \mathbf{r} \right)={\rho}^{2}\left( \mathbf{r} \right)\). See illustrative applications in Int. J. Quantum Chem., 113, 2589 (2013).
101 Positive part of ESP:\({V}_{\mathrm{ESP}}^{+}\). In the region where \(V_{\mathrm{ESP}}\) is positive,\({V}_{\mathrm{ESP}}^{+}={V}_{\mathrm{ESP}}\); where ESP is negative,\({V}_{\mathrm{ESP}}^{+}=0\).
102 Negative part of ESP:\({V}_{\mathrm{ESP}}^{-}\), defined similarly as \({V}_{\mathrm{ESP}}^{+}\).
103 Magnitude of electric field |F|. Since electric field vector is simply negative gradient vector of ESP, therefore this quantity corresponds to norm of gradient of ESP.
110 Total energy density of the steric, electrostatic and quantum components defined in Shubin Liu's energy decomposition analysis (see Section 3.24.2 for introduction), corresponding to sum of user-defined functions 40, 68 and 69
111 Sum of steric potential, electrostatic potential and quantum potential. Note that "ispecial", "iDFTxcsel" and "iKEDsel" in settings.ini affect this and next two functions, see description of user-defined function 60 for detail.
112 Magnitude of vector sum of steric force, electrostatic force and quantum force.
113 Sum of steric charge, electrostatic charge and quantum charge.
114 Pauli kinetic energy density: \({\tau}_\mathrm{θ}(\mathbf{r})={\tau}_{\mathrm{S}}(\mathbf{r})-{\tau}_{\mathrm{W}}(\mathbf{r})\), where τS(r) is the kinetic energy density that can be chosen by "iKEDsel" in settings.ini (see later), and τW(r) is Weizsäcker kinetic energy density.
115 Stiffness \(S=|{\lambda}_{1}|/|{\lambda}_{3}|\), where λ1 < λ2 < λ3 are the three eigenvalues of Hessian matrix of electron density. See J. Comput. Chem., 37, 2722 (2016).
116, 117, 118: Stress tensor stiffness \({S}_{\sigma}=|{\lambda}_{1\sigma}|/|{\lambda}_{3\sigma}|\), stress tensor polarizability \(P=|{\lambda}_{3\sigma}|/|{\lambda}_{1\sigma}|\), stress tensor ellipticity \({\epsilon}_{\sigma H}=|{\lambda}_{1\sigma}|/|{\lambda}_{2\sigma}|-1\), respectively. λ1σ < λ2σ < λ3σ are the three eigenvalues of stress tensor σ. See J. Comput. Chem., 37, 2722 (2016) for \({S}_{\sigma}\) and P, and http://sobereva.com/wfnbbs/viewtopic.php?pid=4667 for \({\epsilon}_{\sigma H}\).
For interested users, the stress tensor σ is briefly outlined here. It was defined by Bader in J. Chem. Phys., 73, 2871 (1980) as follows (expressed in a.u. here)
explicitly, its element can be written as
where \(r_{1}\)=x, \(r_{2}\)=y, \(r_{3}\)=z, and the one-electron reduced density matrix in real space form can be evaluated based on orbitals as \({\Gamma}^{\left( 1 \right)}(\mathbf{r},{\mathbf{r}}^{'})=\sum_{i} {\eta}_{i}{\varphi}_{i}^{*}(\mathbf{r}){\varphi}_{i}(\mathbf{r'})\). Currently only real wavefunctions are considered, so elements of stress tensor are readily evaluated in Multiwfn as
If you need to obtain the entire stress tensor, you can input coordinate in main function 1, stress tensor will be outputted along with other quantities.
The path-packets analysis, which is very closely related to stress tensor, can be realized using the shell script contributed by Asdrubal Lozada based on patched Multiwfn code, check http://sobereva.com/wfnbbs/viewtopic.php?pid=4698 for relevant information and https://github.com/aslozada/Stress_tensor.
200 Random number of [0,1).
819 Ultrastrong interaction (USI): \(\mathrm{USI}(\mathbf{r})={\mathrm{\nabla }}^{2}\rho(\mathbf{r})/{\rho}^{5/3}(\mathbf{r})\). See J. Phys. Chem. A, 126, 2437 (2022) for detail.
820 Bonding and noncovalent interaction (BNI): \(\mathrm{BNI(}\mathbf{r})=[G\left( \mathbf{r} \right)-{\tau}_{\mathrm{W}}\left( \mathbf{r} \right)]/{\tau}_{\mathrm{W}}(\mathbf{r})\), where G(r) is Lagrangian kinetic energy density. See J. Phys. Chem. A, 126, 2437 (2022) for detail.
900, 901, 902 X, Y and Z coordinate variables, respectively.
Functions 910~914 shown below calculate atomic weighting function. The atom to be calculated is determined by “uservar” in settings.ini. Periodic system is supported except for Becke atomic weighting function.
910 Hirshfeld atomic weighting function.
911,912 Becke atomic weighting function using Tian Lu covalent radii and CSD covalent radii, respectively. Sharpness parameter of the boundary is 3
913 Tian Lu error function type atomic weighting function
914 Tian Lu Gaussian function type atomic weighting function
999 Local Hartree-Fock exchange energy (or Hartree-Fock exchange energy density). Its integral over the whole space corresponds to Hartree-Fock exchange energy.
Closed-shell situation:
where
in which χ is basis function, P is total density matrix.
Open-shell situation: \(e_{\mathrm{HFX}}\) is the sum of its alpha and beta parts
where
1000 Integrand of DFT exchange-correlation functionals, which is also known as exchange-correlation energy density.
1100 DFT exchange-correlation potential, only available for closed-shell systems: \(V_{\mathrm{XC}}(\mathbf{r})=\delta E_{\mathrm{XC}}/\delta\rho(\mathbf{r})\), where \(E_{\mathrm{XC}}\) is exchange-correlation energy
1101 DFT exchange-correlation potential of α electrons, only available for open-shell systems:
1102 DFT exchange-correlation potential of β electrons, only available for open-shell systems:
For user-defined functions 1000 and 1100/1101/1102, "iDFTxcsel" in settings.ini is used to select XC functional. 0-29 are X part only, 30-69 are C part only, 70-99 are entire XC. For example, if iDFTxcsel=32, then iuserfunc=1000 corresponds to integrand of LYP correlation functional, and iuserfunc=1100 corresponds to LYP potential.
[See below for details]
Available options of "iDFTxcsel" parameter
Currently only LSDA and some GGA functionals are supported
0 LSDA exchange: \(-(3/2)[3/(4\pi)]^{1/3}\left[\rho_\alpha(\mathbf{r})^{4/3}+\rho_\beta(\mathbf{r})^{4/3}\right]\). For closed-shell cases the equivalent form is \(-(3/4)(3/\pi)^{1/3}\rho(\mathbf{r})^{4/3}\).
1 Becke 88 (B88) exchange, see Phys. Rev. A, 38, 3098 (1988).
2 Perdew-Burke-Ernzerhof (PBE) exchange, see Phys. Rev. Lett., 77, 3865 (1996).
3 Perdew-Wang 91 (PW91) exchange, see Electronic Structure of Solids '91; Ziesche, P., Eschig, H., Eds.; Akademie Verlag: Berlin, 1991; p. 11.
30 Vosko-Wilk-Nusair V (VWN5) correlation, see Can. J. Phys., 58, 1200 (1980).
31 Perdew 86 (P86) correlation, see Phys. Rev. B, 33, 8822 (1986).
32 Lee-Yang-Parr (LYP) correlation, see Phys. Rev. B, 37, 785 (1988).
33 Perdew-Wang 91 (PW91) correlation, see the reference of its exchange counterpart.
34 Perdew-Burke-Ernzerhof (PBE) correlation, see the reference of its exchange counterpart.
70 Becke 97 (B97) exchange-correlation, see J. Chem. Phys., 107, 8554 (1997).
71 Hamprecht-Cohen-Tozer-Handy with 407 training molecules (HCTH407) exchange-correlation, see J. Chem. Phys., 114, 5497 (2001).
80 SVWN5 exchange-correlation.
81 BP86 exchange-correlation.
82 BLYP exchange-correlation.
83 BPW91 exchange-correlation.
84 PBEPBE exchange-correlation.
85 PW91PW91 exchange-correlation.
About iuserfunc=1200¶
The kinetic energy density (KED) denotes integrand of kinetic energy functional. User-defined function 1200 is a collection of all KEDs supported by Multiwfn, which were systematically examined in J. Chem. Phys., 150, 204106 (2019). The "iKEDsel" in settings.ini is used to choose the form of KED.
Note that if "iKEDsel" is not at its default value (0), then the τS term (Lagrangian kinetic energy density) of ELF and strong covalent interaction index (SCI) will be replaced with corresponding KED.
There is a subfunction 92 in the main function 1000 (hidden in the main menu), it can calculate integral of all KEDs over the whole space, namely evaluating kinetic energy based on the kinetic energy functionals. The GEA4 KED is calculated only when "iKEDsel" is set to 24 since it needs Laplacian of electron density, while other KEDs only need electron density and its gradient. Integral of Hamiltonian KED is not calculated by this function since the result is identical to that of Lagrangian KED.
-
iKEDsel=1: Hamiltonian KED, identical to real space function 6
-
iKEDsel=2: Lagrangian KED, identical to real space function 7
-
iKEDsel=3: Thomas-Fermi KED: \(\tau_{\mathrm{TF}}(\mathbf{r})=\sum_{\sigma=\alpha,\beta}\tau_{\mathrm{TF}}^\sigma=\sum_{\sigma=\alpha,\beta}C_{\mathrm{TF}}[\rho_\sigma(\mathbf{r})]^{5/3}\), where \(C_{\mathrm{TF}}=\frac{3}{10}(6\pi^2)^{2/3}=4.557799872\) is Thomas-Fermi constant in spin polarized case.
-
iKEDsel=4: Weizsäcker KED: \(\tau_\mathrm{W}(\mathbf{r})=\sum_{\sigma=\alpha,\beta}\frac{\left|\nabla\rho_\sigma(\mathbf{r})\right|^2}{8\rho_\sigma(\mathbf{r})}\)
Most of below-mentioned KEDs can be represented using a general form \(\tau(\mathbf{r})=\sum_{\sigma=\alpha,\beta}\tau_{\mathrm{TF}}^\sigma(\mathbf{r})F^\sigma[s^\sigma(\mathbf{r})]\), where \(s^\sigma(\mathbf{r})=\left|\nabla\rho^\sigma(\mathbf{r})\right|/\rho^\sigma(\mathbf{r})^{4/3}\), F is called as enhancement factor, see J. Chem. Phys., 127, 144109 (2007) for detail. This paper systematically introduced and compared a variety of existing KEDs. Note that a lot of expressions given in this paper are wrong, while the formulae given below are absolutely correct, and all of them are explicitly written as spin polarized form, if they will be involved in your works, please cite my paper J. Chem. Phys., 150, 204106 (2019); citation of each KED is also given in this paper. More information and comparison about KEDs can be found in Phys. Rev. A, 46, 6920 (1992) and J. Chem. Phys., 100, 4446 (1994).
-
iKEDsel=5: Second order gradient expansion approximation, \(F_{\mathrm{GEA2}}^\sigma=\left[1+\frac{1}{72C_{\mathrm{TF}}}s^\sigma(\mathbf{r})^2\right]\)
-
iKEDsel=6: Thomas-Fermi + 1/5 Weizsäcker KED, \(F_{\mathrm{TF5W}}^\sigma=\left[1+\frac{1}{40C_{\mathrm{TF}}}s^\sigma(\mathbf{r})^2\right]\)
-
iKEDsel=7: Thomas-Fermi + Weizsäcker KED, \(F_{\mathrm{TFvW}}^\sigma=\left[1+\frac{1}{8C_{\mathrm{TF}}}s^\sigma(\mathbf{r})^2\right]\)
-
iKEDsel=8: Thomas-Fermi + b/9 Weizsäcker KED, \(F_{\mathrm{TF9W}}^\sigma=\left[1+\frac{1.067}{72C_{\mathrm{TF}}}s^\sigma(\mathbf{r})^2\right]\)
-
iKEDsel=9: N-dependent Thomas-Fermi KED, \(F_{\mathrm{TF-N}}^\sigma=1+\frac{0.313}{N^{1/3}}-\frac{0.187}{N^{2/3}}\), where N is the number of total electrons in the system
-
iKEDsel=10: Pearson KED, \(\tau_{\mathrm{Pear}}^\sigma=\tau_{\mathrm{TF}}^\sigma+\frac{1}{1+[s_r^\sigma(\mathbf{r})/\zeta]^6}\frac{1}{72}\frac{\left|\nabla\rho^\sigma(\mathbf{r})\right|^2}{\rho^\sigma(\mathbf{r})}\), where ζ is fixed to 1, \(s_r^\sigma(\mathbf{r})=s^\sigma(\mathbf{r})/[2(6\pi^2)^{1/3}]\), similarly hereinafter
-
iKEDsel=11: DePristo-Kress Pade KED, \(F_{\mathrm{DK\ Pade}}^\sigma=\frac{9b_3x^4+a_3x^3+a_2x^2+a_1x+1}{b_3x^3+b_2x^2+b_1x+1}\), where \(x=[s^\sigma(\mathbf{r})]^2/(72C_{\mathrm{TF}})\), \(a_{1}\) = 0.95, \(a_{2}\) = 14.28111, \(a_{3}\) = −19.57962, \(b_{1}\) = −0.05, \(b_{2}\) = 9.99802 and \(b_{3}\)=2.96085
-
iKEDsel=12: Lee-Lee-Parr KED, \(F_{\mathrm{LLP}}^\sigma=1+\frac{0.0044188[s^\sigma(\mathbf{r})]^2}{1+0.0253s^\sigma(\mathbf{r})\mathrm{arcsinh}[s^\sigma(\mathbf{r})]}\)
-
iKEDsel=13: Ou-Yang–Levy 1 KED, \(F_{\mathrm{OL1}}^\sigma=1+\frac{1}{72C_{\mathrm{TF}}}[s^\sigma(\mathbf{r})]^2+0.00187s^\sigma(\mathbf{r})\)
-
iKEDsel=14: Ou-Yang–Levy 2 KED, \(F_{\mathrm{OL2}}^\sigma=1+\frac{1}{72C_{\mathrm{TF}}}[s^\sigma(\mathbf{r})]^2+\frac{0.0245s^\sigma(\mathbf{r})}{1+2^{5/3}s^\sigma(\mathbf{r})}\)
-
iKEDsel=15: Thakkar KED, \(F_{\mathrm{thak}}^\sigma=1+\frac{0.0055[s^\sigma(\mathbf{r})]^2}{1+0.0253s^\sigma(\mathbf{r})\mathrm{arcsinh}[s^\sigma(\mathbf{r})]}-\frac{0.072s^\sigma(\mathbf{r})}{1+2^{5/3}s^\sigma(\mathbf{r})}\)
-
iKEDsel=16: Becke 86A KED, \(F_{\mathrm{B86A}}^\sigma=1+0.0039\frac{[s^\sigma(\mathbf{r})]^2}{1+0.004[s^\sigma(\mathbf{r})]^2}\)
-
iKEDsel=17: Becke 86B KED, \(F_{\mathrm{B86B}}^\sigma=1+0.00403\frac{[s^\sigma(\mathbf{r})]^2}{\{1+0.007[s^\sigma(\mathbf{r})]^2\}^{4/5}}\)
-
iKEDsel=18: DePristo-Kress 87 KED, \(F_{\mathrm{DK87}}^\sigma=1+0.00132327[s^\sigma(\mathbf{r})]^2\frac{1+0.861504s^\sigma(\mathbf{r})}{1+0.044286[s^\sigma(\mathbf{r})]^2}\)
-
iKEDsel=19: Perdew-Wang 86 KED, \(F_{\mathrm{PW86}}^\sigma=\{1+1.296[(s_r^\sigma(\mathbf{r}))^2+14(s_r^\sigma(\mathbf{r}))^4+0.2(s_r^\sigma(\mathbf{r}))^6]\}^{1/15}\)
-
iKEDsel=20: Perdew-Wang 91 KED, \(F_{\mathrm{PW91}}^\sigma=\frac{1+a_1s_r^\sigma(\mathbf{r})\mathrm{arcsinh}[b\times s_r^\sigma(\mathbf{r})]+\{a_2-a_3e^{-100[s_r^\sigma(\mathbf{r})]^2}\}[s_r^\sigma(\mathbf{r})]^2}{1+a_1s_r^\sigma(\mathbf{r})\mathrm{arcsinh}[b\times s_r^\sigma(\mathbf{r})]+a_4[s_r^\sigma(\mathbf{r})]^4}\), where \(a_{1}\)=0.19645, \(a_{2}\)=0.2743, \(a_{3}\)=0.1508, \(a_{4}\)=0.004 and b=7.7956
-
iKEDsel=21: Lacks-Gordon 94 KED \(F_{\mathrm{CSM}}^\sigma=\frac{(1+a_1s_r^\sigma(\mathbf{r})^2+a_2[s_r^\sigma(\mathbf{r})]^4+a_3[s_r^\sigma(\mathbf{r})]^6+a_4[s_r^\sigma(\mathbf{r})]^8+a_5[s_r^\sigma(\mathbf{r})]^{10}+a_6[s_r^\sigma(\mathbf{r})]^{12})^2}{1+10^{-4}[s_r^\sigma(\mathbf{r})]^2}\), where \(a_{2}\)=(10-8+0.1234)/0.024974, \(a_{4}\)=29.790, \(a_{6}\)=22.417, \(a_{8}\)=12.119, \(a_{10}\)=1570.1, \(a_{12}\)=55.944 and b=0.024974
-
iKEDsel=22: Acharya-Bartolotti-Sears-Parr KED, \(F_{\mathrm{ABSP}}^\sigma=\frac{1}{8C_{\mathrm{TF}}}[s^\sigma(\mathbf{r})]^2+1-\frac{1.412}{N^{1/3}}\)
-
iKEDsel=23: Gázquez-Robles KED, \(F_{\mathrm{GR}}^\sigma=\frac{1}{8C_{\mathrm{TF}}}[s^\sigma(\mathbf{r})]^2+\left(1-\frac{2}{N}\right)\left(1-\frac{1.303}{N^{1/3}}+\frac{0.029}{N^{2/3}}\right)\)
-
iKEDsel=24: Fourth order gradient expansion approximation,
\[\tau_{\mathrm{GE4}}^\sigma=\tau_{\mathrm{GEA2}}^\sigma+\frac{(6\pi^2)^{2/3}}{540}[\rho^\sigma(\mathbf{r})]^{1/3}\left\{\left[\frac{\nabla^2\rho^\sigma(\mathbf{r})}{\rho^\sigma(\mathbf{r})}\right]^2-\frac{9}{8}\left(\frac{\nabla^2\rho^\sigma(\mathbf{r})}{\rho^\sigma(\mathbf{r})}\right)\left(\frac{\nabla\rho^\sigma(\mathbf{r})}{\rho^\sigma(\mathbf{r})}\right)^2+\frac{1}{3}\left|\frac{\nabla\rho^\sigma(\mathbf{r})}{\rho^\sigma(\mathbf{r})}\right|^4\right\}\]
Note that if the "uservar" in settings.ini is not equal to zero, all above-mentioned KEDs will be added by "∇2ρ /uservar" term. For example, when iuserfunc=1200, iKEDsel=5 and uservar=6, the user-defined function will correspond to τGEA2+∇2ρ /6 (which corresponds to the KED employed by Tsirelson type of ELF and LOL, as shown in Section 2.6).
1201 Difference between KED selected by "iKEDsel" and Weizsäcker KED.
1202 Difference between KED selected by "iKEDsel" and Lagrangian KED.
1203 Absolute difference between KED selected by "iKEDsel" and Lagrangian KED.
1204 Local temperature T(r)=[2τ(r)]/[3ρ(r)]. The kinetic energy density τ(r) used in this function is determined by "iKEDsel". In contrast, the user-defined function 7 uses Lagrangian kinetic energy density. Note that if ρ is less than the “uservar” parameter in settings.ini, then T will be treated as 0.
1210 Potential of kinetic energy functional, whose form is determined by "iKEDsel" in settings.ini. Only valid for closed-shell case and available for the following KEDs:
iKEDsel=3: Thomas-Fermi
iKEDsel=5: Second order GEA
iKEDsel=7: Thomas-Fermi + Weizsäcker
Main real space function as user-defined function¶
For flexibility consideration, main real space functions can also be invoked as user-defined function. If you set “iuserfunc” in settings.ini to 10000+i, then the ith real space function will be chosen as the user-defined function. For example, if “iuserfunc” is set to 10009, then the 9th real space function (ELF) will be chosen as the user-defined function.