Density functional theory (DFT) is an effective tool in computational physics that allows for the prediction and analysis of numerous transport and thermal properties of solids1,2. The precision and reliability of DFT computations are greatly influenced by the choice of exchange–correlation functional. Correlation energy is a fundamental concept in quantum mechanics that accounts for the interactions between electrons in a many-electron system. It represents the additional energy required to describe the behavior of electrons beyond what can be explained by the mean-field approximation, such as the Hartree–Fock method3. It basically represents the effects of electron–electron correlation, which arises due to their mutual electrostatic repulsion. This correlation leads to complex quantum effects that cannot be represented by simple mathematical models. Various functionals have been developed over the years, each with its own unique advantages and disadvantages4,5. While some functionals may provide higher overall accuracy, others may be better suitable for specific kinds of systems or features. To accurately correct the positive value of Homo-energy and the orbital energy of anions, for instance, the tuned range-separated DFT approximation has been introduced by satisfying the DFT-Koopmans theory based on basis set alterations and quantum approximations6. The choice of functional depends on the particular problem being solved, the degree of required precision, and the computation time. Fitting the results of accurate Quantum Monte Carlo (QMC) computations of the uniform electron gas (UEG) has yielded multiple equations for the correlation energy in the local density approximation (LDA). These are quite complex expressions, usually denoted by their acronyms, like VWN, VWN5, and so forth. Nevertheless, when applied to molecular computation, they yield remarkably similar results. VWN functional is defined as follows7:$$ E_{c}^{VWN} = \int {d^{3} r(A\left\{ \begin{gathered} \ln \frac{{x^{2} }}{X(x)} + \frac{2b}{Q}\tan^{ – 1} \frac{Q}{2x + b} – \frac{{bx_{0} }}{{X(x_{0} )}} \hfill \\ [\ln \frac{{(x – x_{0} )^{2} }}{X(x)} + \frac{{2(b + 2x_{0} )}}{Q}\tan^{ – 1} \frac{Q}{2x + b}] \hfill \\ \end{gathered} \right\}} \,)\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, $$
(1)
where \(x = r_{s}^{1/2} \,,\,X(x) = x^{2} + bx + c\,,\,Q = (4c – b^{2} )^{1/2}\) and the parameters \({x}_{0}\) , b, and c are constants that depend on the context and the specific version of the VWN functional being used8. The error may be obtained at a relatively high level for the LDA correlation functional. An improved LDA leads to LSDA0 that is constructed to satisfy specific requirements, but remarkably agrees with favorable results for a homogeneous two-electron density in a restricted curved three-dimensional area. The correlation part of LSDA0 is expressed as follows:$$ E_{c}^{LSDA0} = \int {d^{3} rn\frac{{ – b_{1c} }}{{1 + b_{2c} r_{s}^{\frac{1}{2}} + b_{3c} r_{s} }}g_{c} (\xi )} $$
(2)
where, \(z = {{(n_{ \uparrow } – n_{ \downarrow } )} \mathord{\left/ {\vphantom {{(n_{ \uparrow } – n_{ \downarrow } )} {(n_{ \uparrow } + n_{ \downarrow } )}}} \right. \kern-0pt} {(n_{ \uparrow } + n_{ \downarrow } )}}\) is the relative spin polarization defined with the spin densities \({n}_{\uparrow }\) and \({n}_{\downarrow }\), \(r_{s}\) is the Seitz radius, b is constant and the function \({g}_{c}(z)\) is 0 for z =  ± 1 in order to satisfy the one electron self-correlation constraint and 1 for z = 09,10. The use of Generalized Gradient Approximation (GGA) can increase accuracy. Some of GGA correlations functional are LYP, PBE, and PW91. The PW91 correlation functional has achieved widespread acceptance because it covers short- and long-range electron correlation effects and accomplishes a comparatively low Mean Absolute Error (MAE) for the overall energy for a wide range of molecules11,12,13. The well-known PBE functional correlation part is given by10,14:$$ E_{c}^{PBE} = \int {n(r)\varepsilon_{c}^{PBE} (n(r))dr} $$
(3)
where \({\varepsilon }_{c}^{PBE}(n(r))\) , is the correlation energy density. The precise mathematical expression for \({\varepsilon }_{c}^{PBE}(n(r))\) is quite complex and can be found in the original publication by Pedrew, Burke, and Ernzerhof14. This functional’s error is better than LDA but still higher than other functionals. The accuracy of the LYP functional combined with Becke’s three-parameter exchange functional, or B3LYP, is also better than that of more basic functionals like LDA, especially for molecule systems and transition metals15,16.Vydrov and Van Voorhis reported that the B97M-V functional was developed for main-group thermochemistry and non-covalent interactions. These families of functions were created by Truhlar and associates, and they exhibit improved accuracy for the total energy, particularly for non-covalent interactions17,18,19. For both organic and inorganic molecules, it has demonstrated remarkable performance for the total energy with low MAE values20,21,22. The mathematical formulation of the B97M-V functional is extremely complex and comprises several parameters that are determined through experimentation22.A different correlation functional that was developed by Chachiyo et al.23 yields impressively accurate findings when combined with their exchange functional24,25,26:$$ E_{c} = \int {n\varepsilon_{c} (1 + t^{2} )^{{\frac{h}{{\varepsilon_{c} }}}} d^{3} r} $$
(4)
where t (\(t = (\frac{\pi }{3})^{1/6} \frac{1}{4}\frac{{\left| {\vec{\nabla }n} \right|}}{{n^{7/6}}}\)) is the gradient parameter, n is the electron density, \({\varepsilon }_{c}\) is the correlation energy density and h is a constant with a value of 0.06672632 Hartree23. Increasing t is intended to decrease the maximum rate of energy through the gradient suppressing factor \({(1+{t}^{2})}^{\frac{h}{{\varepsilon }_{c}}}\) .In this paper, we report a new correlation functional that minimizes the mean absolute error (MAE), exceeding existing methods in terms of accuracy. The MAE is a measure of the average difference between computed and experimental values and is an essential measure of a functional’s performance. Specifically, we perform a detailed analysis of over 64 different molecules and evaluate their key characteristics, including bond energies, total energy, dipole moments, and zero-point energy. We derive new correlation functional analytically; taking into account the functional’s dependency on both the electron density and the ionization energy. By incorporating the ionization energy into the correlation functional, we improve the accuracy. This methodology is a complementary one to our earlier research, wherein we presented a novel exchange functional that is dependent on the ionization energy27. As a result, the ionization energy is included as a significant parameter in both the correlation and exchange functional in our current study, enabling a more comprehensive description of electronic interactions27.Theoretical methodEquation (4) defines the generic form of the correlation functional23. Based on the second-order Moller–Plesset perturbation (MP2) theory, a straightforward analytical \({\varepsilon }_{c}\) is defined in the local-density approximation as follows28:$$ \varepsilon_{c} = a\ln (1 + \frac{b}{{r_{s} }} + \frac{b}{{r_{s}^{2} }})\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, $$
(5)
where a and b are constants and \({r}_{s}\) is the Wigner–Seitz density parameter.In this work, we used the ionization dependent density as29,30:$$ n(r_{s} ) \to Ar_{s}^{2\beta } e^{{ – 2(2I)^{\frac{1}{2}} r_{s} }} \, $$
(6)
where I is the ionization energy and \(\beta =\frac{1}{2}\sqrt{\frac{2}{I}} -1\) . We take A = 1 into consideration based on formula 29 in reference29. Next, we determine \({r}_{s}\) to be:$$ \begin{gathered} r_{s} = – \frac{{\sqrt{\frac{2}{I}} – 2}}{{2\sqrt {2I} }}W( – \frac{{2\sqrt {2I} }}{{\sqrt{\frac{2}{I}} – 2}}n^{{\frac{1}{{\sqrt{\frac{2}{I}} – 2}}}} ) = \,\, – \frac{{\sqrt{\frac{2}{I}} – 2}}{{2\sqrt {2I} }}( – \frac{{2\sqrt {2I} }}{{\sqrt{\frac{2}{I}} – 2}})\,n^{{\frac{1}{{\sqrt{\frac{2}{I}} – 2}}}} W(n) \hfill \\ = \,n^{{\frac{1}{{\sqrt{\frac{2}{I}} – 2}}}} W(n)\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \hfill \\ \end{gathered} $$
(7)
We obtain the following, using Eqs. (6) and (7):$$ \begin{gathered} \vec{\nabla }n = 2A\beta (n^{{\frac{1}{{\sqrt{\frac{2}{I}} – 2}}}} W(n)\,)^{2\beta – 1} e^{{ – 2(2I)^{\frac{1}{2}} n^{{\frac{1}{{\sqrt{\frac{2}{I}} – 2}}}} W(n)\,}} \hfill \\ – 2A(2I)^{\frac{1}{2}} (n^{{\frac{1}{{\sqrt{\frac{2}{I}} – 2}}}} W(n)\,)^{2\beta } e^{{ – 2(2I)^{\frac{1}{2}} n^{{\frac{1}{{\sqrt{\frac{2}{I}} – 2}}}} W(n)\,}} \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \hfill \\ \end{gathered} $$
(8)
So:$$ \begin{gathered} t = (\frac{\pi }{3})^{1/6} \frac{{\sqrt {a_{0} } }}{4}n^{{\frac{2\beta – 1}{{\sqrt{\frac{2}{I}} – 2}} – \frac{7}{6}}} (2A\beta (W(n)\,)^{2\beta – 1} e^{{ – 2(2I)^{\frac{1}{2}} n^{{\frac{1}{{\sqrt{\frac{2}{I}} – 2}}}} W(n)\,}} \hfill \\ – 2A(2I)^{\frac{1}{2}} (n^{{\frac{ – 1}{{\sqrt{\frac{2}{I}} – 2}}}} )(W(n)\,)^{2\beta } e^{{ – 2(2I)^{\frac{1}{2}} n^{{\frac{1}{{\sqrt{\frac{2}{I}} – 2}}}} W(n)\,}} )\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \hfill \\ \end{gathered} $$
(9)
Remarkably, \({\varepsilon }_{c}\) can be interpreted as an interpolation between the paramagnetic and ferromagnetic instances for spin polarization \(z = \frac{{n_{\alpha } – n_{\beta } }}{n}\)31:$$ \varepsilon_{c} (n,z) = \varepsilon_{c}^{0} + (\varepsilon_{c}^{1} – \varepsilon_{c}^{0} )f(z)\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, $$
(10)
where \(g^{3} (z) = 1 – \frac{1}{2}f(z)\quad and \quad g(z) = \frac{(1 + z)^{2/3} + (1 – z)^{2/3}}{2}\).Thus, we rewrite Eq. (5) as:$$ \varepsilon_{c} = a\ln (1 + \frac{b}{W(n)}n^{{\frac{ – 1}{{\sqrt{\frac{2}{I}} – 2}}}} + \frac{b}{{(W(n))^{2} }}n^{{\frac{ – 2}{{\sqrt{\frac{2}{I}} – 2}}}} )\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, $$
(11)
Which is clearly depends on ionization energy8. The Lambert function32 denoted as W, is related to density n. The correlation energy when the density tends to be uniform is represented by \({\varepsilon }_{c}\) in Eq. (5) that is associated with density n28. We obtain the following by substituting Eqs. (9) and (11) in Eq. (4):$$ \begin{gathered} E_{c} = \int {na\ln (1 + \frac{b}{W(n)}n^{{\frac{ – 1}{{\sqrt{\frac{2}{I}} – 2}}}} + \frac{b}{{(W(n))^{2} }}n^{{\frac{ – 2}{{\sqrt{\frac{2}{I}} – 2}}}} )} (1 + ((\frac{\pi }{3})^{1/6} \frac{{\sqrt {a_{0} } }}{4}n^{{\frac{2\beta – 1}{{\sqrt{\frac{2}{I}} – 2}} – \frac{7}{6}}} (2A\beta (W(n)\,)^{2\beta – 1} e^{{ – 2(2I)^{\frac{1}{2}} n^{{\frac{1}{{\sqrt{\frac{2}{I}} – 2}}}} W(n)\,}} \hfill \\ – 2A(2I)^{\frac{1}{2}} (n^{{\frac{ – 1}{{\sqrt{\frac{2}{I}} – 2}}}} )(W(n)\,)^{2\beta } e^{{ – 2(2I)^{\frac{1}{2}} n^{{\frac{1}{{\sqrt{\frac{2}{I}} – 2}}}} W(n)\,}} )\,\,)^{2} )^{{\frac{h}{{a\ln (1 + \frac{b}{W(n)}n^{{\frac{ – 1}{{\sqrt{\frac{2}{I}} – 2}}}} + \frac{b}{{(W(n))^{2} }}n^{{\frac{ – 2}{{\sqrt{\frac{2}{I}} – 2}}}} )}}}} d^{3} r\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \hfill \\ \end{gathered} $$
(12)
This is our new correlation energy in its final form. In order to get the exchange–correlation functional, as previously stated, this relationship was paired with the exchange functional that obtained in our prior study27, namely with:$$\Rightarrow {E}_{x}=\int \rho {\varepsilon }_{x}\frac{\begin{array}{c}\left(\sqrt{9{\pi }^{2}{\left(2-\sqrt{\frac{2}{I}} \right)}^{-2}+6}-3\pi {\left(2-\sqrt{\frac{2}{I}} \right)}^{-1}\right){x}^{2}+({\frac{3}{4\pi }\left(\sqrt{9{\pi }^{2}{\left(2-\sqrt{\frac{2}{I}} \right)}^{-2}+6}-3\pi {\left(2-\sqrt{\frac{2}{I}} \right)}^{-1}\right)}^{2}{x}^{2}\\ +\left(\sqrt{9{\pi }^{2}{\left(2-\sqrt{\frac{2}{I}} \right)}^{-2}+6}-3\pi {\left(2-\sqrt{\frac{2}{I}} \right)}^{-1}\right)x+\frac{4\pi }{3})\left(2-\sqrt{\frac{2}{I}} \right)\left(L\left(-\alpha {\left(\frac{3}{4\pi }\right)}^{\frac{-3}{\sqrt{\frac{2}{I}} -2}}{x}^{\frac{-3}{\sqrt{\frac{2}{I}} -2}}\right)+1\right)\end{array}}{\left(\left(\sqrt{9{\pi }^{2}{\left(2-\sqrt{\frac{2}{I}} \right)}^{-2}+6}-3\pi {\left(2-\sqrt{\frac{2}{I}} \right)}^{-1}\right)x+\frac{4\pi }{3}\right)\left(2-\sqrt{\frac{2}{I}} \right)\left(L\left(-\alpha {\left(\frac{3}{4\pi }\right)}^{\frac{-3}{\sqrt{\frac{2}{I}} -2}}{x}^{\frac{-3}{\sqrt{\frac{2}{I}} -2}}\right)+1\right)}{d}^{3}r$$
(13)
where, \(\varepsilon_{x} = \frac{{ – 3\sqrt {2I} }}{4\pi s}\) , the parameter \(x = \frac{4\pi }{3}s\) and \(s = \frac{{2\sqrt {2I} }}{{2(3\pi^{2} )^{\frac{1}{3}} }}A^{{\frac{ – 1}{3}}} e^{{\frac{2}{3}r\sqrt {2I} }} r^{{\frac{ – 1}{3}\sqrt{\frac{2}{I}} + \frac{2}{3}}}\).We utilized Siam quantum software to do numerical calculations and examine the exchange–correlation functional results33,34. Reference34 was utilized while determining the values of I for atoms and molecules.
Introducing a new correlation functional in density functional theory
