Recent Advances in Cartesian-Grid DFT in Atoms and Molecules

In the past several decades, density functional theory (DFT) has evolved as a leading player across a dazzling variety of fields, from organic chemistry to condensed matter physics. The simple conceptual framework and computational elegance are the underlying driver for this. This article reviews some of the recent developments that have taken place in our laboratory in the past 5 years. Efforts are made to validate a viable alternative for DFT calculations for small to medium systems through a Cartesian coordinate grid- (CCG-) based pseudopotential Kohn–Sham (KS) DFT framework using LCAO-MO ansatz. In order to legitimize its suitability and efficacy, at first, electric response properties, such as dipole moment (μ), static dipole polarizability (α), and first hyperpolarizability (β), are calculated. Next, we present a purely numerical approach in CCG for proficient computation of exact exchange density contribution in certain types of orbital-dependent density functionals. A Fourier convolution theorem combined with a range-separated Coulomb interaction kernel is invoked. This takes motivation from a semi-numerical algorithm, where the rate-deciding factor is the evaluation of electrostatic potential. Its success further leads to a systematic self-consistent approach from first principles, which is desirable in the development of optimally tuned range-separated hybrid and hyper functionals. Next, we discuss a simple, alternative time-independent DFT procedure, for computation of single-particle excitation energies, by means of “adiabatic connection theorem” and virial theorem. Optical gaps in organic chromophores, dyes, linear/non-linear PAHs, and charge transfer complexes are faithfully reproduced. In short, CCG-DFT is shown to be a successful route for various practical applications in electronic systems.

1 Introduction

“The general theory of quantum mechanics is now almost complete⋯ The underlying physical laws necessary for the mathematical theory of a large part of physics and the whole of chemistry are thus completely known, and the difficulty is only that the exact application of these laws leads to equations much too complicated to be soluble. It therefore becomes desirable that approximate practical methods of applying quantum mechanics should be developed, which can lead to an explanation of the main features of complex atomic systems without too much computation” (Dirac, 1929). This famous and enlightening announcement in 1929 made physicists and chemists rise to the challenge and develop the theoretical framework needed to calculate wave functions and properties of atoms, molecules, and solids. Throughout the course of the past few decades, ab initio quantum mechanical methods have been utilized for elucidation of electronic structure, properties, and dynamics of many-electron systems. With the rapid progress of computational algorithms, in association with modern computer architecture and resources, some of the recent advancements have been extended to widely explored fields such as materials science, nanoscience, and biological science. The electronic Schrödinger equation (SE) of such an N-electron system is, in principle, a many-body problem having space, spin, and time variables, expressed as

$\mathbf{H}{\psi }_i\left({\mathbf{x}}_1,{\mathbf{x}}_2,\dots ,{\mathbf{x}}_N\right)=E_i{\psi }_i\left({\mathbf{x}}_1,{\mathbf{x}}_2,\dots ,{\mathbf{x}}_N\right),\mathbf{x}=\left\{\mathbf{r},\sigma \right\},(1)$

where H is the many-electron Hamiltonian operator consisting of M nuclei and N electrons, ψ_i is the N-electron anti-symmetric wave function corresponding to ith eigenstate of H, and E_i is its energy. However, to pursue it practically, one has to go through multiple challenges caused due to the size of many-body SE. The exact analytical solution is very hard to obtain in most cases, leaving aside a few well-known model systems. The quantum mechanical solvability of N-electron systems is exhausted by hydrogen and helium atom. That is why the question of how to deal with systems containing thousands of electrons and hundreds of nuclei has attained utmost relevance. The straightforward application of SE, either in its numerical, variational, or perturbation theory versions, is nowadays out of the reach of even the most advanced supercomputers. It is for this reason that alternative ways of handling such problems have been vigorously pursued during the last few decades. In the past few years, significant strides have been made in these aforementioned areas.

In today’s computational chemistry, it is desirable to achieve the energy of a chemical reaction within the bounds of chemical accuracy ($<$1 kcal mol⁻¹). The primary challenge is to reduce the computational cost without much compromise on accuracy for increasingly larger and complicated systems. Now, starting from less accurate, reliable empirical or semi-empirical schemes to modern ab initio methods, all fall under the roof of the current computational chemistry repertoire. Two fundamentally different routes, based on wave function approaches and functional theories, have gained popularity and credence so far.

The usual quantum-mechanical approach to SE can be summarized as follows:

$v\left(\mathbf{r}\right)\stackrel{\mathrm{SE}}{\Rightarrow }\mathrm{\Psi }\left({\mathbf{r}}_1,{\mathbf{r}}_2\dots \dots .,{\mathbf{r}}_N\right)\stackrel{<\mathrm{\Psi }|\dots |\mathrm{\Psi }>}{\Rightarrow }\text{observables}.(2)$

In other words, one specifies a system by choosing v(r), plugs it into SE, solves that equation for a wave function Ψ, and then calculates observables by taking the expectation values of operators. Many powerful approximate methods for solving SE have been developed. Thus, starting from the single determinantal uncorrelated Hartree–Fock (HF) to various correlated multi-configurational methods is available. It is well-known that the HF method is still not accurate enough for energy predictions in chemistry; for example, bond energies are significantly underestimated. Thus, post-HF methods, adding numerous other determinants involving excited or “virtual” orbitals, are generally required for a satisfactory prediction. Some notable correlated methods are diagrammatic perturbation theory (based on Feynman diagrams and Green’s functions), Möller–Plesset perturbation theory (MPn), configuration interaction (CI), coupled-cluster ansatz (comes in many flavors such as CCSD, CCSD(T), CCSDT, and CCSDTQ), and multi-reference perturbation theory (such as CASPT2), among others. These methods offer authentic and reliable results but are quite difficult to be implemented computationally for large N mainly because of their unfavorable scaling.

The above-mentioned methods based on an approximation to many-electron wave function were the only possibilities until 1964, the birth year of density functional theory (DFT). In general, the so-called functional theories rely on utilizing the limited information coming from single-particle electron density, ρ(r), reduced density matrix or Green’s function, and follow variational principle. Amongst them, however, DFT has emerged as the most versatile and outstanding candidate in electronic structure theory. This leads to the central quantity of our present article, namely, the spin-less, single-particle electron density ρ(r), which is the diagonal element of one-particle density matrix, defined as

$\rho \left(\mathbf{r}\right)=N\int \dots .\int \mathrm{\Psi }\left(\mathbf{r}\sigma ,{\mathbf{x}}_2,\dots \dots {\mathbf{x}}_N\right){\mathrm{\Psi }}^{\ast }\left(\mathbf{r}\sigma ,{\mathbf{x}}_2,\dots \dots {\mathbf{x}}_N\right)d\sigma d{\mathbf{x}}_2\dots ..d{\mathbf{x}}_N.(3)$

Walter Kohn noted in his Nobel lecture that “DFT has been most useful for systems of very many electrons where wave function methods encounter and are stopped by the exponential wall” (Kohn, 1999). Many beautiful books and in-depth reviews are available on the subject (Parr and Yang, 1989; Koch and Holthausen, 2001; Martin, 2004; Pugh and Hinchliffe, 2006; Hoffman, 2007; Champagne and Springborg, 2009; Burke, 2012; Roy, 2012; Zangwill, 2015; Yu et al., 2016; Mardirossian and Head-Gordon, 2017; Janesko, 2021).

The first genuine DFT scheme for electronic systems was given as early as in 1927 (Thomas, 1927); it was a model for calculating atomic properties, based purely on ρ(r). It assumed that electrons form a gas satisfying Fermi statistics, with electron–electron interaction energy determined from classical Coulomb potential. For kinetic energy, they adopted a local density approximation (LDA), where the contribution from a point r is determined from kinetic energy of a homogeneous electron gas with this density. The resulting Euler equation is

$\frac{5}{3}{C}_F{\rho }^{\frac{2}{3}}\left(\mathbf{r}\right)+\int d{\mathbf{r}}^{\prime }\frac{\rho \left({\mathbf{r}}^{\prime }\right)}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}+v_{\mathit{ext}}\left(\mathbf{r}\right)+\lambda =0,(4)$

where $C_F=\frac{3}{10}{(3{\pi }^2)}^{\frac{2}{3}}$, v_ext(r) is external potential, and λ is the Lagrange multiplier related to the constraint of constant particle number. The Thomas–Fermi (TF) model has severe deficiencies because of poor description of outer region of an atom. The charge density decays as r⁻⁶ far from nucleus, not exponentially as it should. Moreover, it does not bind neutral atoms or form molecules or solids. Later, Dirac (1930) noted the necessity of incorporating “exchange” phenomena, by recasting HF theory in terms of a “density function,” without reference to a single-determinantal wave function. This function leads to a correction to the energy derived from exchange energy for a homogeneous electron gas of ρ(r). The modified Thomas–Fermi–Dirac (TFD) energy density functional can be written as

$E_{\mathit{TFD}}\left[\rho \left(\mathbf{r}\right)\right]=C_F\int \rho {\left(\mathbf{r}\right)}^{\frac{5}{3}}d\mathbf{r}+\int \rho \left(\mathbf{r}\right)v_{\mathit{ext}}\left(\mathbf{r}\right)d\mathbf{r}-C_x\int \rho {\left(\mathbf{r}\right)}^{\frac{4}{3}}d\mathbf{r}+\frac{1}{2}\iint \frac{\rho \left(\mathbf{r}\right)\rho \left({\mathbf{r}}^{\prime }\right)}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}.(5)$

The overwhelming simplicity of abandoning a complicated many-electron SE for a single equation in terms of ρ(r) alone is remarkably appealing. However, underlying approximations are rather primitive, making it grossly inadequate for any practical use in quantum chemistry.

In 1964–1965, the Hohenberg–Kohn (HK) (Hohenberg and Kohn, 1964) and Kohn–Sham (KS) (Kohn and Sham, 1965) theorems, the twin pillars of modern DFT, were published. They asked a plain obvious question of whether the information contained in ρ(r) is sufficient for elucidating a many-electron system completely. In these seminal works, they founded the rigorous theory that finally legitimized the intuitive leaps of other previous works. They first proved that the external potential, v_ext(r), of a non-degenerate ground state of an N-particle system (and hence total energy), is a unique functional of ground-state electron density ρ(r). This one-to-one mapping, which is the basic preamble of this theorem, can be expressed by the following energy functional:

$E_{v_{\text{ext}}}\left[\rho \right]=\langle \mathit{\psi }\left[\rho \right]|\mathbf{T}+{\mathbf{V}}_{\text{ext}}+{\mathbf{V}}_{\text{ee}}|\mathit{\psi }\left[\rho \right]\rangle =F_{\text{HK}}\left[\rho \right]+\int v_{\text{ext}}\left(\mathbf{r}\right)\rho \left(\mathbf{r}\right)d\mathbf{r},(6)$

where F_HK[ρ] denotes the universal energy density functional containing kinetic energy and electron–electron interaction terms. Further, the total energy of a given system can be determined variationally by minimizing the functional $E_{v_{\mathit{ext}}}[\rho (\mathbf{r})]$, subject to the normalization condition, ∫ρ(r) = N, as a constant, through the following equation:

$\frac{\delta F_{\text{HK}}\left[\rho \right]}{\delta \rho \left(\mathbf{r}\right)}+v_{\text{ext}}\left(\mathbf{r}\right)=\mathit{\mu }.(7)$

Moreover, to make sure that a given density is indeed the true ground-state density, now the second HK theorem suggests

$E_{v_{\text{ext}}}\left[\rho \left(\mathbf{r}\right)\right]\ge E_{v_{\text{ext}}}\left[{\rho }_0\left(\mathbf{r}\right)\right],(8)$

where $E_{v_{\mathit{ext}}}[{\rho }_0(\mathbf{r})]$ corresponds to the ground-state energy of a Hamiltonian with v_ext(r) as external potential and ρ₀(r) is its ground-state density.

The explicit form of F_HK(ρ) is unknown as yet. Hence, this needs to be evaluated approximately. However, so far, the most successful way to implement DFT is through a method proposed in Kohn and Sham (1965). They introduced the clever concept of a fictitious, non-interacting system built from a set of KS orbitals, such that the major part of kinetic energy can be computed exactly. The remaining fairly small portion is absorbed in the non-classical contribution of electron–electron repulsion, which is also unknown. However, the advantage is that electrons now move in an effective KS single-particle potential. The mapped auxiliary system now yields the same ground-state density as the real interacting system. However, this simplifies the actual calculation tremendously, as it is operationally an independent-particle theory, simpler even than HF. Yet, it delivers, in principle, the exact density the same as ground-state density resulting from a summation of moduli of square of orbitals through

${\rho }_s\left(\mathbf{r}\right)=\sum _i^N\sum _{\sigma }{\left|{\mathrm{\Psi }}_i\left(\mathbf{r},\sigma \right)\right|}^2={\rho }_0\left(\mathbf{r}\right).(9)$

Here, σ signifies spin, and the exact total energy is expressed as

$E\left[\rho \right]=\int v_{\mathit{ext}}\left(\mathbf{r}\right)\rho \left(\mathbf{r}\right)d\mathbf{r}+J\left[\rho \right]+T_s\left[\rho \right]+E_{xc}\left[\rho \right],(10)$

where

$E_{xc}\left[\rho \right]=\left(T\left[\rho \right]-T_s\left[\rho \right]\right)+\left(E_{ee}\left[\rho \right]-J\left[\rho \right]\right)=T_c\left[\rho \right]+E_{nc}\left[\rho \right].(11)$

The associated terms have the following meanings: J[ρ] is the known classical part of E_ee[ρ], whereas E_xc[ρ] contains everything unknown, that is, non-classical electrostatic effects of E_ee[ρ] and the difference between true kinetic energy T[ρ] and T_s[ρ]. The latter represents the sum of individual kinetic energies of reference system with same density as real system as $T_s=-\frac{1}{2}{\sum }_i^N{\nabla }_i^2$ and T_c[ρ] symbolizes the correlation kinetic energy.

Thus, as apparent from above, the beauty of DFT lies in its being exact yet efficient, with a single determinant describing the ρ(r)—the whole complexity is hidden in one term, the exchange-correlation (XC) functional. Local exchange functionals in KS theory automatically include some static correlation (Cook and Karplus, 1987; Handy and Cohen, 2001). Thus, a single Slater determinant is not necessarily as poor in KS theory as in HF theory, keeping the burden at the same level as HF. Building better and improved XC functionals to treat real correlated systems within KS theory is an active area of research (Peverati and Truhlar, 2014; Dale et al., 2017; Verma et al., 2019; Wang et al., 2019; Verma and Truhlar, 2020), so much so that the story behind the success of DFT, to a large extent, is tantamount to the search for accurate reliable XC functional. The exact form remains elusive; it is necessary to use various density functional approximations (DFAs). The popular DFAs can be hierarchically characterized in the following manner. The simplest XC functionals are of LDA-type (Kohn and Sham, 1965) (containing ρ only), residing in the first rung of Jacob’s ladder. They are exact for an infinite uniform electron gas but are highly inaccurate for molecular properties because most real systems have inhomogeneous density distribution. Next, generalized gradient approximation (GGA) functionals (Becke and Dickson, 1988; Perdew et al., 1996a) (with the addition of gradient of electron density, ∇ρ) occupy the second rung of the ladder and tend to improve significantly upon LDA. In addition to combining separable exchange and correlation functionals (e.g., PBE, BP86, BLYP, PW91, rev-PBE, and RPBE), it is possible to semi-empirically parameterize GGA XC functionals (HCTH, N12, and GAM) (Boese and Handy, 2001; Peverati and Truhlar, 2012a; Yu et al., 2015). Then, the third rung belongs to meta-GGA (Tao et al., 2003; Zhao and Truhlar, 2008) functionals (addition of kinetic energy density, τ). With further inclusion of exact exchange (EEX) energy, one gets the hybrid functionals, occupying the fourth rung of the ladder (Perdew and Schmidt, 2000). We also have functionals that go beyond the fourth rung (including virtual orbitals), thus requiring an even higher computational cost.

Now with these insightful backgrounds, we proceed to investigate some aspects of the realistic solution of the KS equation. Therefore, one has to deal with several mathematically non-trivial integrals that cannot be evaluated analytically and pose a certain degree of difficulty. Without such procedures, one is left with no choice but to calculate numerically. It is well recognized that such a discrete procedure for multi-center integrals in 3D space is difficult. It becomes even more daunting when it is noted that more computation time and a larger grid size are often required to achieve a satisfactory level of chemical accuracy due to the involvement of a prohibitively huge number of operations. Cusps in the density and singular nature of Coulombic potential offer a major challenge when constructing such integration routes. Therefore, to accomplish high-accuracy calculations within a reasonable number of quadrature grids, one needs numerical integrators that are both efficient and sophisticated enough to capture the forms of density at a reasonable level. This opens up a wide range of integrators with varying degrees of performance. Of these, two distinct, well-recognized partitioning schemes have shown a great deal of promise: the Voronoi cellular approach and fuzzy cells Avenue, commonly known as the atom-centered grid (ACG). Currently, the enviable status of DFT is beholden to the ensuing basis-set calculations. As these studies are heavily dominated by ACG, the real-space grid has been invoked for fully numerical, basis-set free DFT methods. Apart from ACG, a few scattered works exist for different grids in literature, for example, an adaptive Cartesian grid with a hierarchical cubature technique, a transformed sparse-tensor product grid, and a Fourier Transform Coulomb method interpolating density from ACG to a more regular grid. This work presents a simple fruitful way to implement DFT within the basis-set framework, utilizing only CCG.

Thus, within the Born–Oppenheimer approximation and Hohenberg–Kohn–Sham framework, an implementation of DFT is offered in CCG. With the aid of a linear combination of Gaussian functions, molecular orbitals (MO) and quantities such as basis functions, ρ(r), as well as classical Hartree and non-classical XC potentials, are directly calculated on a 3D real CCG. A combination of Fast Fourier Transform (FFT) and inverse FFT is used to calculate the Coulomb potential quite accurately. In order to present the inner core electrons, analytical effective core potentials are used, whereas energy-optimized truncated Gaussian bases are used for valence electrons. The requisite work of a many-electron problem has four different proceeding directions: method development, improvement of respective new and existing theories, successful implementation of them, and finally, application of these theories in various physicochemical problems. All four genres are covered in this review. Section 2 forms the theoretical backbone for works presented throughout, followed by a systematic investigation of static linear response properties for a host of atoms and molecules in Section 3. Within the finite-field (FF) KS framework, several properties such as permanent dipole moment (μ), dipole polarizability (α), and first hyperpolarizability (β) are considered within the InDFT program (Roy et al., 2019) developed in our laboratory. A simple alternative of the FF technique, employing a rational function fit to μ with respect to the electric field, is also acquired. Next, a purely numerical, efficient computation scheme of HF exchange energy, density, and matrix elements for certain orbital-dependent and range-separated hybrid (RSH) functionals is presented in Section 4. This is inspired by a recently developed algorithm by Liu and Kong (2017), where an accurate evaluation of the electrostatic potential (ESP) integral is the rate-determining step. A combination of the Fourier convolution theorem (FCT) with an RS Coulomb interaction kernel (CIK) is introduced. The latter is efficiently mapped onto a real grid through an easy optimization procedure, leading to a constraint within the RS parameter, allowing a logarithmic scaling with respect to atomic/molecular size. Simultaneously, as an offshoot of this work, in Section 5, a self-consistent systematic optimization procedure, from first principles, is proposed for the development of optimally tuned (OT) RSH functionals through a size-dependency-based ansatz. To this end, a novel self-consistent procedure appears, which engages the characteristic length of a system with the RS parameter. Finally, Section 6 details the successful implementation of Becke’s exciton model, followed by its applications in computing the optical gap in organic chromophores and various properties of charge-transfer complexes. This is an alternative time-independent DFT procedure to compute single-particle excitation energies, which are of particular relevance in photochemistry. Correct evaluation of the correlated singlet-triplet splitting (STS) energy is the key step in this procedure. Non-empirical models from both the “adiabatic connection theorem” and “virial theorem” are introduced for such calculations. The obtained results are compared with theoretical and experimental references as appropriate. Finally, we end with a few comments in Section 7.

2 DFT in Cartesian Grid

In this segment, we briefly outline the methodology and the computational aspects. The single-particle KS equation of a many-electron system in presence of pseudopotential can be written as follows (atomic unit employed unless stated otherwise):

$\left[-\frac{1}{2}{\nabla }^2+v_{\text{ion}}^{pp}\left(\mathbf{r}\right)+v_h\left[\rho \left(\mathbf{r}\right)\right]+v_{xc}\left[\rho \left(\mathbf{r}\right)\right]\right]{\varphi }_i\left(\mathbf{r}\right)={ϵ}_i{\varphi }_i\left(\mathbf{r}\right),(12)$

where $v_{\mathit{ion}}^{\text{pp}}$ designates ionic pseudopotential, expressed as

$v_{\text{ion}}^{\text{pp}}=\sum _{R_a}{v}_{\text{ion},a}^{pp}\left(\mathbf{r}-{\mathbf{R}}_a\right).(13)$

In the above equation, $v_{\text{ion},a}^{\text{pp}}$ represents the ion-core pseudopotential associated with atom a, situated at R_a, whereas $v_h[\rho (\mathbf{r})]=\int \frac{\rho ({\mathbf{r}}^{\prime })}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}d{\mathbf{r}}^{\prime }$ describes usual classical electrostatic (Hartree) potential among valence electrons, and v_xc[ρ(r)], as usual, represents the non-classical part of many-electron Hamiltonian, H as $\frac{\delta E_{\text{xc}}[\rho (\mathbf{r})]}{\delta \rho (\mathbf{r})}$, the functional derivative of XC energy. The single-particle charge density is then given by

$\rho \left(\mathbf{r}\right)={\rho }^{\alpha }\left(\mathbf{r}\right)+{\rho }^{\beta }\left(\mathbf{r}\right)=\sum _i{f}_i^{\alpha }|{\varphi }_i^{\alpha }\left(\mathbf{r}\right){|}^2+\sum _i{f}_i^{\beta }|{\varphi }_i^{\beta }\left(\mathbf{r}\right){|}^2,(14)$

where $\{{\varphi }_i^{\sigma },\sigma =\alpha \mathrm{o}\mathrm{r}\beta \}$ corresponds to a set of N occupied orthonormal spin molecular orbitals (MO) and $f_i^{\sigma }$’s denote occupation numbers of ith spin-MO. The integral of this pseudo-valance density offers the total number of valence electrons, N_v. The benefits of pseudopotential are twofold. Firstly, as the number of KS orbitals relies exclusively on N_v only, it is particularly favorable for heavy elements. In such cases, where each atom involves several tens of electrons, N_v can be much smaller. Thus, pseudopotential can consider the non-negligible relativistic impacts in heavy elements. Secondly, as pseudo-valance orbitals are typically smoother than core orbitals, much lesser grid points suffice.

For the realistic solution of Eq. 12, the basis-set technique is by far the most convenient and pragmatic one. This is essentially motivated from the success of basis-set related methodologies in conventional wave function (such as HF and post-HF) theory through LCAO-MO ansatz. Consequently, the KS MOs are expanded in terms of K appropriately picked, known basis functions {χ_μ(r); μ = 1, 2, 3, …., K}, often called atomic orbitals (AO), in a practice homologous to that in the Roothaan-HF method, such as

${\varphi }_i\left(\mathbf{r}\right)=\sum _{\mu =i}^K{C}_{\mu i}{\chi }_{\mu }\left(\mathbf{r}\right),i=\mathrm{1,2,3},\dots .K.(15)$

The electron density then takes the following expression in this basis

$\rho \left(\mathbf{r}\right)=\sum _{i=1}^{N_v}\sum _{\mu =1}^K\sum _{\nu =1}^K{C}_{\mu i}{C}_{\nu i}{\chi }_{\mu }\left(\mathbf{r}\right){\chi }_{\nu }\left(\mathbf{r}\right).(16)$

In principle, one requires a complete basis set (K = ∞) to get exact expansion of MOs, but it is not achievable in reality. Subsequently, appropriate truncation is needed for practical computational purposes; it suffices to work with a finite basis set.

Now embedding Eq. 15 in Eq. 12, multiplying left side of resulting equation with ${\chi }_{\nu }^{*}(\mathbf{r})$, and then integrating over entire space, trailed by some algebraic manipulation gives rise to the following KS matrix equation, in parallel to the HF case:

${\mathbf{F}}^{\text{KS}}\mathbf{C}=\mathbf{S}\mathbf{C}\mathit{ϵ},(17)$

where F and S stand for K × K real, symmetric total KS and overlap matrices, respectively. The eigenvector matrix C contains basis-set expansion coefficients C_μi and diagonal matrix ϵ holds orbital energies ϵ_i. It could be promptly solved by standard numerical techniques of linear algebra. Individual components of KS matrix can be expressed as

$\begin{array}{ccc}\hfill F_{\mu \nu }^{\text{KS}}& =\hfill & \int {\chi }_{\mu }\left(\mathbf{r}\right)\left[h^{\text{core}}+v_{\text{hxc}}\left(\mathbf{r}\right)\right]{\chi }_{\nu }\left(\mathbf{r}\right)\mathrm{d}\mathbf{r}\hfill \\ \hfill & =\hfill & H_{\mu \nu }^{\text{core}}+\langle {\chi }_{\mu }\left(\mathbf{r}\right)|v_{\text{hxc}}\left(\mathbf{r}\right)|{\chi }_{\nu }\left(\mathbf{r}\right)\rangle =H_{\mu \nu }^{\text{core}}+J_{\mu \nu }+V_{\mu \nu }^{\text{xc}}\hfill \\ \hfill v_{\text{hxc}}\left(\mathbf{r}\right)& =\hfill & v_{\text{h}}\left(\mathbf{r}\right)+v_{\text{x}}\left(\mathbf{r}\right)+v_{\text{c}}\left(\mathbf{r}\right),\hfill \end{array}(18)$

where $H_{\mu \nu }^{\text{core}}$ includes the core bare-nucleus Hamiltonian matrix element comprising of kinetic energy and nuclear-electron attraction, representing 1e⁻ contributions. These can be evaluated analytically with the assistance of well-established recursion relations (Obara and Saika, 1986) for Gaussian type orbital (GTO) bases, which is employed here (see later). The second term v_hxc(r) contains all 2e⁻ interactions including classical Coulomb repulsion and non-classical XC potential. J_μν signifies matrix element of v_h(r), whereas the remaining term, $V_{\mu \nu }^{\text{xc}}$ supplies XC contribution into 2e⁻ matrix element, whose development remains one of the fundamental steps in the entire KS-DFT process. In absence of any analytical method, the latter can be either computed numerically or fitted by an auxiliary set of Gaussian functions (Sambe and Felton, 1975; Dunlap et al., 1979a; Dunlap et al., 1979b). The current procedure does not engage any fitting. For gradient-corrected functionals, non-local XC contribution of KS matrix is implemented by a finite-orbital basis expansion, without requiring an evaluation of the density Hessians. Thus, in such cases, XC contribution is written in a convenient working form, given in (Pople et al., 1992)

$F_{\mu \nu }^{\text{xc}\alpha }=\int \left[\frac{\partial f}{\partial {\rho }_{\alpha }}{\chi }_{\mu }{\chi }_{\nu }+\left(2\frac{\partial f}{\partial {\gamma }_{\alpha \alpha }}\nabla {\rho }_{\alpha }+\frac{\partial f}{\partial {\gamma }_{\alpha \beta }}\nabla {\rho }_{\beta }\right)\cdot \nabla \left({\chi }_{\mu }{\chi }_{\nu }\right)\right]d\mathbf{r},(19)$

where γ_αα = |∇ρ_α|², γ_αβ = ∇ ρ_α ⋅∇ ρ_β, γ_ββ = |∇ρ_β|², and f is a function only of local quantities ρ_α, ρ_β, and their gradients.

To continue further, all relevant quantities, such as basis function, electron density, MO, and all 2e⁻ potentials are straightforwardly set up on a real 3D CCG:

$r_i=r_0+\left(i-1\right)h_r,i=\mathrm{1,2,3},\dots .,N_r;\mathrm{f}\mathrm{o}\mathrm{r}r\in \left\{x,y,z\right\},(20)$

where h_r, N_r imply grid spacing and total number of grid points, respectively, $(r_0=-\frac{N_r{h}_r}{2})$. As necessary, the single-particle density in the active grid can be expressed as

$\rho \left({\mathbf{r}}_g\right)=\sum _{\mu \nu }{P}_{\mu \nu }{\chi }_{\mu }\left({\mathbf{r}}_g\right){\chi }_{\nu }\left({\mathbf{r}}_g\right),(21)$

where r_g is the gth grid point in the simulation box. Similarly, any multi-centre integration involved in KS-DFT, such as classical Hartree energy and XC energy, in principle, can be directly set up in 3D real-space grid utilizing small 3D unit volume:

$I=\int F\left(\mathbf{r}\right)\mathrm{d}\mathbf{r}\approx \sum _{i}F\left({\mathbf{r}}_i\right)h_x{h}_y{h}_z.(22)$

The 2e⁻ matrix elements can be acquired by direct numerical integration in CCG:

$\langle {\chi }_{\mu }\left(\mathbf{r}\right)|v_{\mathit{hxc}}\left(\mathbf{r}\right)|{\chi }_{\nu }\left(\mathbf{r}\right)\rangle \approx h_x{h}_y{h}_z\sum _{\text{grid}}{\chi }_{\mu }\left(\mathbf{r}\right)v_{\mathit{hxc}}\left(\mathbf{r}\right){\chi }_{\nu }\left(\mathbf{r}\right).(23)$

The pragmatic implementation of Eq. 22 is a lot less involved than that in ACG. We note that the presence of a cusp in density and singularity in Coulomb potential would create some computational burden if F(r) is not smooth enough throughout the simulation box.

A pressing issue in the grid-based approach comprises an accurate estimation of v_h(r). For finite systems, the simplest and crudest route for computing this is through direct numerical integration on the grid. For smaller systems, this is a feasible option; in any remaining cases, it is generally tedious and cumbersome. However, the most rewarding and widespread approach is through a solution of the corresponding Poisson equation:

${\nabla }^2{v}_h\left(\mathbf{r}\right)=-4\pi \rho \left(\mathbf{r}\right).(24)$

The standard method for tackling this is by conjugate gradients (Saad, 2003) or multi-grid solvers (Brandt, 1977). As another option, a conventional FCT was originally suggested by Martyna and Tuckerman (1999), Mináry et al. (2002), Skylaris et al. (2002) and adapted in the context of molecular modeling (Hine et al., 2011; Chang et al., 2012). To give a layout of this theorem, let us start with two functions f(r) and F(k) in r and k spaces, which are Fourier transforms (FT) (denoted by $\mathcal{F}$) of one another as follows:

$\begin{array}{ccc}\hfill f\left(\mathbf{r}\right)& =\hfill & \mathcal{F}\left(\mathbf{k}\right)={\int }_{-\infty }^{\infty }F\left(\mathbf{k}\right)e^{2\pi i\mathbf{r}\mathbf{k}}d\mathbf{k}\hfill \\ \hfill F\left(\mathbf{k}\right)& =\hfill & {\mathcal{F}}^{-1}\left(\mathbf{r}\right)={\int }_{-\infty }^{\infty }f\left(\mathbf{r}\right)e^{-2\pi i\mathbf{r}\mathbf{k}}d\mathbf{r}.\hfill \end{array}(25)$

With two functions f(r), g(r), one can frame the convolution, characterized by

$f\star g={\int }_{-\infty }^{\infty }f\left({\mathbf{r}}^{\prime }\right)g\left(\mathbf{r}-{\mathbf{r}}^{\prime }\right)d{\mathbf{r}}^{\prime }.(26)$

The FCT directs that FT of convolution is just the product of individual FTs:

$\mathcal{F}\left(f\star g\right)=F\left(\mathbf{k}\right).G\left(\mathbf{k}\right),f\star g={\mathcal{F}}^{-1}\left(F\left(\mathbf{k}\right).G\left(\mathbf{k}\right)\right).(27)$

The above theorem can be utilized to construct v_h(r) efficiently as follows:

$v_h\left(\mathbf{r}\right)=FF{T}^{-1}\left\{v_h^c\left(\mathbf{k}\right)\rho \left(\mathbf{k}\right)\right\}\mathrm{a}\mathrm{n}\mathrm{d}\rho \left(\mathbf{k}\right)=FFT\left\{\rho \left(\mathbf{r}\right)\right\},(28)$

where $v_h^c(\mathbf{k})$ and ρ(k) represent Fourier integrals of the CIK and density. The quantity ρ(k) can be handily calculated using discrete FT of its real-space value. Accordingly, our essential worry here lies in the calculation of CIK, which has a singularity in real space. For this, we use an Ewald summation-type approach (Chang et al., 2012), expanding the Hartree kernel into long-range (LR) and short-range (SR) parts:

$v_h^c\left(\mathbf{r}\right)=\frac{\mathrm{e}\mathrm{r}\mathrm{f}\left(\zeta r\right)}{r}+\frac{\mathrm{e}\mathrm{r}\mathrm{f}\mathrm{c}\left(\zeta r\right)}{r}\equiv v_{h_{\text{lr}}}^c\left(\mathbf{r}\right)+v_{h_{\text{sr}}}^c\left(\mathbf{r}\right),(29)$

where erf(x) and erfc(x) denote error function and its complementary function, respectively. The FT of SR part can be dealt with analytically, whereas the LR segment necessitates to be computed directly from FFT of real-space values. A convergence parameter ζ is utilized to adjust the range of $v_{h_{\text{sr}}}^c(\mathbf{r})$, so that the error is minimized. Following the conjecture of Martyna and Tuckerman (1999), here we employ ζ × L = 7 (L indicates smallest side of the box), which produces quite precise outcomes (Ghosal et al., 2016; Ghosal et al., 2018; Ghosal et al., 2019). A few other courses that evaluate LR part proficiently are fast multi-pole (Beatson and Greengard, 1997), multi-level summation (Skeel et al., 2002), and fast Fourier-Poisson (York and Yang, 1994) method, among others.

Next, we follow a simple grid-optimization technique as follows:

${\left\{N_x,N_y,N_z\right\}}_i\iff {\rho }_i^{\text{SCF}}\iff E_{\text{tot},i}^{\text{SCF}},(30)$

for a fixed grid spacing h_r. Here, “i” denotes the ith combination of N_x, N_y, N_z for the active box. With this grid parameter, the corresponding self-consistent field (SCF) density and total energy can be labelled as ${\rho }_i^{\text{SCF}}$ and $E_{\text{tot},i}^{\text{SCF}}$, respectively. Now, one can systematically increase the value of {N_x, N_y, N_z} from an appropriate starting point and eventually reach the optimal value of {N_x, N_y, N_z} such that

${\left\{N_x,N_y,N_z\right\}}_i={\left\{N_x,N_y,N_z\right\}}_{\text{opt}},\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{n}\mathrm{\Delta }E=\left(E_{\text{tot},i}^{\text{SCF}}-E_{\text{tot},i-1}^{\text{SCF}}\right)<\mathrm{t}\mathrm{h}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{h},(31)$

where thresh is the grid accuracy for total energy convergence, that is, the energy difference between two successive calculations with different grid parameters. A detailed demonstration of this simple optimization strategy has been well documented (Ghosal et al., 2016; Ghosal et al., 2018; Mandal et al., 2019).

For practical, useful electronic structure calculation, it is of foremost significance to choose suitable functions that imitate KS MOs as precisely as possible. Numerical accuracy of KS-DFT is very delicate to the choice and design of a basis set for a particular problem, as an incomplete basis set inducts certain restrictions on the relaxation of density through KS orbitals. A sizeable number of elegant, flexible, versatile basis sets have been proposed from various perspectives, of which GTOs stand out as our most appealing choice. Moreover, for a practical purpose, rather than involving individual GTOs as the basis, it is customary to utilize a fixed linear combination of GTOs, called contracted GTOs. It is cordial in terms of ease and efficiency of computation of essential integrals. While man attractive choices exist for full calculations, which contain higher angular momentum orbitals, the option is much restricted for pseudopotential approximations. We have employed the following effective core potential basis sets: SBKJC (Stevens et al., 1984) for species containing Group-II elements, LANL2DZ (Wadt and Hay, 1985) for Group-III or higher group elements, and Labello–Ferreira–Kurtz (LFK) basis as proposed in Labello et al. (2005), based in the light of a method to incorporate diffuse and polarization functions in familiar Sadlej basis set (Sadlej, 1992). These are adopted from EMSL Basis Set Library (Feller, 1996). All 1e⁻ integrals are generated by standard recursion relations (Obara and Saika, 1986) utilizing GTOs as primitive basis functions. The norm-conserving pseudopotential matrix elements on a contracted basis are imported from the GAMESS (Schmidt et al., 1993) program package. The discrete Fourier transforms are incorporated from the FFTW3 package (Frigo and Johnson, 2005). The above features are implemented in the InDFT (Roy et al., 2019) program developed in our laboratory over the years, which has been employed in several practical applications in a series of articles (Roy, 2008a; Roy, 2008b; Roy, 2010a; Roy, 2010b; Roy, 2011; Ghosal et al., 2016; Ghosal et al., 2018; Ghosal et al., 2019; Mandal et al., 2019; Ghosal et al., 2021; Roy et al., 2021; Ghosal and Roy, 2022a; Ghosal and Roy, 2022b; Roy et al., 2022). At this point, this includes some of the frequently used, prominent XC functionals mentioned at relevant places in the article. The following sections summarize some of the applications of InDFT that have taken place in our laboratory in almost the last 5 years.

3 Electric Response Properties

A salient feature of atoms, molecules, and clusters is the electric dipole polarizability; in other words, their ability to respond to an external electric field (George, 2006; Champagne and Springborg, 2009). An accurate description of this has a prominent role in exploring various interesting phenomena of field-matter interaction and inter-particle collision, such as Rayleigh and Raman scattering, second-order Stark effect, and electron detachment process (El Ghazaly et al., 2005). The linear and nonlinear electric properties, such as μ, α, β are highly relevant in many applications, for example, the development of nonlinear optical materials, structural identification of atomic clusters, Raman and infrared spectroscopy, and separations of molecular isomers. Note that these symbols have also been used for basis set indices. However, there should be no confusion, as their meanings would be apparent from the context of their usage. From a technological point of view, it is interesting to synthesize and design novel optical materials and molecular assemblies with large non-linear optical coefficients (Kümmel and Kronik, 2006). Several distinct theoretical routes were put forward in the literature to obtain these properties within the KS-DFT rubric. Some of the noteworthy ones are the coupled-perturbed Kohn–Sham (CPKS) method (Fournier, 1990; Colwell et al., 1993), linear-response time-dependent DFT (Jansik et al., 2005; Helgaker et al., 2012), perturbative sum-over state expression over all dipole-allowed electronic transitions (Orr and Ward, 1971; Bishop, 1994), the numerical method using the Sternheimer approach (Talman, 2012), auxiliary DFT (Flores-Moreno and Köster, 2008; Carmon-Espíndola et al., 2012), non-iterative CPKS (Shophy et al., 2008) and the fully numerical FF method (Bishop and Pipin, 1987; Maroulis and Thakkar, 1988). The least expensive method, from a computational standpoint, is FF. This approach does not require any analytical derivatives or information about the excited state; implementation is also quite favorable because only the one-body Hamiltonian is perturbed by the applied field (Kurtz et al., 1990). These are the reasons for its success and popularity over other methods (Bulat et al., 2005; de Wergifosse et al., 2014; Wouters et al., 2016).

3.1 FF KS Method

The response of a many-electron system can be represented by expanding field-dependent μ, computed from the field-instigated charge distribution, as a power series in external electric field F (provided the field strength remains small), as

${\mathit{\mu }}_i\left(\mathbf{F}\right)={\mathit{\mu }}_i\left(0\right)+\sum _j{\mathit{\alpha }}_{ij}{F}_j+\frac{1}{2}\sum _{j,k}{\mathit{\beta }}_{\mathit{ijk}}{F}_j{F}_k+\cdots .(32)$

In this equation, three consecutive terms on the right-hand side designate static dipole moment μ_i(0), dipole polarizability ${\mathit{\alpha }}_{ij}=\frac{\partial {\mathit{\mu }}_i}{\partial F_j}$, and first-hyperpolarizability ${\mathit{\beta }}_{\mathit{ijk}}=\frac{{\partial }^2{\mathit{\mu }}_i}{\partial F_j\partial F_k}$ (Mclean and Yoshimine, 1967), respectively. An alternative representation is also available in terms of field-induced energy; however, both are equivalent according to the Hellmann–Feynman theorem (Feynman, 1939).

The components of α, β can be obtained using well-known finite-difference formulas (Smith, 1978):

$\left.\begin{array}{c}\hfill {\mathit{\alpha }}_{ii}{F}_i=\frac{2}{3}\left[{\mathit{\mu }}_i\left(F_i\right)-{\mathit{\mu }}_i\left(-F_i\right)\right]-\frac{1}{2}\left[{\mathit{\mu }}_i\left(2F_i\right)-{\mathit{\mu }}_i\left(-2F_i\right)\right]\hfill \\ \hfill {\mathit{\alpha }}_{ij}{F}_j=\frac{2}{3}\left[{\mathit{\mu }}_i\left(F_j\right)-{\mathit{\mu }}_i\left(-F_j\right)\right]-\frac{1}{2}\left[{\mathit{\mu }}_i\left(2F_j\right)-{\mathit{\mu }}_i\left(-2F_j\right)\right]\hfill \\ \hfill \mathit{\beta }iii{F}_i^2=\frac{1}{3}\left[{\mathit{\mu }}_i\left(2F_i\right)+{\mathit{\mu }}_i\left(-2F_i\right)\right]-\frac{1}{3}\left[{\mathit{\mu }}_i\left(F_i\right)+{\mathit{\mu }}_i\left(-F_i\right)\right]\hfill \\ \hfill \mathit{\beta }ijj{F}_j^2=\frac{1}{3}\left[{\mathit{\mu }}_i\left(2F_j\right)+{\mathit{\mu }}_i\left(-2F_j\right)\right]-\frac{1}{3}\left[{\mathit{\mu }}_i\left(F_j\right)+{\mathit{\mu }}_i\left(-F_j\right)\right]\hfill \end{array}\right\}.(33)$

Furthermore, in addition to α, β tensors, for a given species, one can also compute the so-called experimentally determined quantity, the average polarizability $(\bar{\alpha })$, in the form of

$\bar{\alpha }=\frac{1}{3}\left({\mathit{\alpha }}_{zz}+{\mathit{\alpha }}_{xx}+{\mathit{\alpha }}_{yy}\right).(34)$

In order to get these tensors from μ of the system (expressed as a function of F), one needs the perturbed density matrix at various field strengths. This can be obtained from the SCF solution of Eq. 13. Hence, the core part of the Hamiltonian (symbolized by a prime) will now be altered by a relevant/appropriate field-dependent term conventionally expressed as

$H_{\mathit{\mu }\nu }^{\prime \mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}}=H_{\mathit{\mu }\nu }^{\text{core}}+F_i\langle \mu |\mathbf{r}|\nu \rangle ,i\in \left\{x,y,z\right\}.(35)$

Here, $H_{\mathit{\mu }\nu }^{\text{core}}$ denotes the unperturbed core Hamiltonian mentioned above, F_i refers to ith component of F, and ⟨μ|r|ν⟩ provides the dipole moment integral related to length vector r. The two-body matrix elements do not change during FF calculations. Eventually, μ of a molecule can be described as follows:

$\mathit{\mu }\equiv {\mathit{\mu }}_{el}=\sum _{\mu \nu }{P}_{\mu \nu }\langle \mu |\mathbf{r}|\nu \rangle +\sum _a{Z}_a{\mathbf{R}}_a,(36)$

where Z_a and R_a are nuclear charge and position of atom “a”, respectively.

In the FF method, numerical accuracy is a crucial factor. The system’s dipole moment is computed in the presence of F, and the respective finite differences are used to approximate the derivatives. Hence, it is very sensitive with respect to F, and the field needs to be chosen with utmost care such that 1) it is adequately large to subdue finite-precision artifacts for a meaningful estimation of essential finite differences, particularly in the nonlinear domain for β and 2) it must also be small enough to be able to neglect the higher-order derivatives for one particular coefficient. Selection of the appropriate field strengths initiates with picking up an initial field strength (F₀), around which other field strengths (F) are distributed. This is achieved by selecting the field distribution according to the following relation:

$F_n=F_0{2}^{\frac{3n}{100}},(37)$

where n ranges from 0 to 160, and F₀ = 0.0005 a.u. It gives a maximum field strength larger than 0.01 a.u. For a given F_n, these properties are calculated for a fixed grid and basis set.

3.2 Field Sensitivity and Its Optimization

To address the delicate nature of electric field on these properties, in addition to the above procedure, we also adopted a recently proposed (Patel et al., 2017) technique, whereby the energy is fitted with respect to electric-field coefficients in the form of a rational function. This is examined by a fitting strategy for induced dipole moment in terms of the electric field as follows:

$\mathit{\mu }\left(\mathbf{F}\right)=\frac{a+b\mathbf{F}+c{\mathbf{F}}^2+d{\mathbf{F}}^3\cdots }{1+B\mathbf{F}+C{\mathbf{F}}^2+D{\mathbf{F}}^3\cdots }.(38)$

where a,b,c,d,⋯ and B,C,D,⋯ are fitting coefficients. If the denominator coefficients are set to zero, this gives rise to a generalized form of Taylor expansion. Such a recipe has the advantage that it provides a less sensitive (thus more effective) dependence on the electric field, as it enlarges the range of the field. In the FF technique, μ needs to be computed at certain field strengths. That requires a proper selection of initial field strength (F₀), which is achieved here via a proposal put forth by Patel et al. (2017):

$F_n=x^n{F}_0,x=2^{\frac{p}{100}}.(39)$

The recommended value of p is 50, corresponding to a geometric progression. This was arrived on the basis of a systematic and detailed analysis of α and γ for a set of 121 and β for 91 molecules. Following Ghosal et al. (2018), an initial value of 10^−2.5 was found to be quite appropriate. The optimal form of rational function is adopted (Patel et al., 2017) as

$\mathit{\mu }\left(\mathbf{F}\right)=\frac{a+b\mathbf{F}+c{\mathbf{F}}^2+d{\mathbf{F}}^3}{1+B\mathbf{F}+C{\mathbf{F}}^2},(40)$

containing four and three terms in numerator and denominator. Now putting value of μ at F = 0 in the above equation leads to a trivial relation, μ(0) = a. The remaining unknown coefficients can be determined employing different F_n. For five unknown coefficients, six minimum equations can be constructed (as both +F_n and −F_n used), each F_n giving two equations. Employing μ(0) for a, one may write the following matrix equation of form, Ax = b:

$\left[\begin{array}{ccccc}\hfill -F_0\mathit{\mu }\left(F_0\right)\hfill & \hfill -F_0^2\mathit{\mu }\left(F_0\right)\hfill & \hfill F_0\hfill & \hfill F_0^2\hfill & \hfill F_0^3\hfill \\ \hfill F_0\mathit{\mu }\left(-F_0\right)\hfill & \hfill -F_0^2\mathit{\mu }\left(-F_0\right)\hfill & \hfill -F_0\hfill & \hfill F_0^2\hfill & \hfill -F_0^3\hfill \\ \hfill -x{F}_0\mathit{\mu }\left(x{F}_0\right)\hfill & \hfill -x^2{F}_0^2\mathit{\mu }\left(x{F}_0\right)\hfill & \hfill x{F}_0\hfill & \hfill x^2{F}_0^2\hfill & \hfill x^3{F}_0^3\hfill \\ \hfill x{F}_0\mathit{\mu }\left(-x{F}_0\right)\hfill & \hfill -x^2{F}_0^2\mathit{\mu }\left(-x{F}_0\right)\hfill & \hfill -x{F}_0\hfill & \hfill x^2{F}_0^2\hfill & \hfill -x^3{F}_0^3\hfill \\ \hfill -x^2{F}_0^2\mathit{\mu }\left(x^2{F}_0^2\right)\hfill & \hfill -x^4{F}_0^4\mathit{\mu }\left(x^2{F}_0^2\right)\hfill & \hfill x^2{F}_0^2\hfill & \hfill x^4{F}_0^4\hfill & \hfill x^6{F}_0^6\hfill \\ \hfill x^2{F}_0^2\mathit{\mu }\left(-x^2{F}_0^2\right)\hfill & \hfill -x^4{F}_0^4\mathit{\mu }\left(-x^2{F}_0^2\right)\hfill & \hfill -x^2{F}_0^2\hfill & \hfill x^4{F}_0^4\hfill & \hfill -x^6{F}_0^6\hfill \end{array}\right]\left[\begin{array}{c}\hfill B\hfill \\ \hfill C\hfill \\ \hfill b\hfill \\ \hfill c\hfill \\ \hfill d\hfill \end{array}\right]=\left[\begin{array}{c}\hfill \mathit{\mu }\left(F_0\right)-\mathit{\mu }\left(0\right)\hfill \\ \hfill \mathit{\mu }\left(-F_0\right)-\mathit{\mu }\left(0\right)\hfill \\ \hfill \mathit{\mu }\left(x{F}_0\right)-\mathit{\mu }\left(0\right)\hfill \\ \hfill \mathit{\mu }\left(-x{F}_0\right)-\mathit{\mu }\left(0\right)\hfill \\ \hfill \mathit{\mu }\left(x^2{F}_0^2\right)-\mathit{\mu }\left(0\right)\hfill \\ \hfill \mathit{\mu }\left(-x^2{F}_0^2\right)-\mathit{\mu }\left(0\right)\hfill \end{array}\right].(41)$

The solution of this equation is overdetermined, as both (+)ve and (−)ve fields are used. A least-squares method can be convenient; else, one may disregard any one equation. The current work invokes the latter, where one of the six equations is arbitrarily eliminated. The required properties can be determined from derivatives of Eq. 38 at F = 0:

$\mathit{\alpha }={\mathit{\mu }}^{\prime }\left(0\right)=b-aB,\mathit{\beta }={\mathit{\mu }}^{″}\left(0\right)=2c-2bB-2aC-2a{B}^2.(42)$

Further details of the method and its validation can be found in Mandal et al. (2019).

3.3 Numerical Tests and Convergence

This section presents some sample results to demonstrate the applicability of the above-described method. At first, a few practical points may be pointed out. All computations are executed involving norm-conserving pseudopotential at the experimental geometries taken from the NIST database (Johnson, 2016). A simple grid-optimization technique has been followed, ensuring a grid precision of at least 5 × 10⁻⁶ a.u., all through, at a fixed h_r = 0.3. It was noticed that the optimal non-uniform grid marginally differs from functional to functional. We employ the LFK basis set for this study. The properties are inspected for four representative XC functionals, namely, LDA, BLYP, PBE, and LBVWN. The SCF convergence criteria imposed in this calculation to generate the perturbed density matrix are as follows: 1) orbital energy difference between two consecutive iterations and 2) absolute deviation in a density matrix element. They both must stay below a specific predetermined threshold, which is set to 1 × 10⁻⁸ a.u. This assured that total energy maintains a convergence to at least this level. In order to facilitate the convergence, an unperturbed (field-free) density matrix was employed as trial input. The convergence was carefully examined with respect to all parameters, such as grid and field optimization, both in the absence and presence of the electric field.

Let us now examine the α and β tensors. Following Buckingham and Hirschfelder (1967), we have two independent components (α_xx=yy, α_zz) associated with α and β (β_xxz=yyz, β_zzz), for a hetero-nuclear diatomic molecule having C_∞v symmetry. The maximum field response toward the electron density is then found along the z direction as it is the molecular axis. Now, as a check, we have performed the GAMESS calculations (Schmidt et al., 1993) with default grid options, that is, 96 radial and 302 angular points for the spatial grid and 0.001 for field strength. A recent study of grid effects (based on ACG), reported in Castet and Champagne (2012), suggested a spatial grid of 99 radial and 974 angular points to be an optimal solution. It has been observed that the default option delivers results that are practically coincident with that from the finer grid; we have verified this for three diatomic molecules (HCl, HBr, and HI). Thus, for all practical purposes, the default grid suffices the current purpose.

3.3.1 Comparison With Standard Packages

Now with this preamble, at first, we report the non-zero components of α, β in addition to $\bar{\alpha }$ and μ, of a representative test molecule, HCl in Table 1. These are supplied for all four XC functionals. We quote the reference values (except for LBVWN) acquired from GAMESS software (Schmidt et al., 1993). Comparative components for a few other selective diatomic molecules are additionally presented by Ghosal et al. (2018), which are excluded here to save space. The largest mean absolute deviation (MAD) in μ_z in HCl is around 5 × 10⁻⁵ a.u., for PBE, whereas the LDA and BLYP results impeccably coincide with reference (in all the digits quoted) for LDA and BLYP. The α, β tensors of our calculations are also similarly consistent with the reference data. For comparison, a few relevant theoretical outcomes are cited in the footnote, alongside the methods (such as higher-order perturbation theory, MCSCF, and CCSD(T)) and basis set. Also, experimental values for $\bar{\alpha }$ are additionally recorded from two different kinds of experimental strategies (Olney et al., 1997; Hohm, 2013). These values contain just the electronic part, and neither geometry relaxation in the presence of the electric field nor vibrational contribution is considered. It uncovers a fascinating fact that all three traditional functionals (LDA, BLYP, and PBE) overestimate both experimental outcomes; however, LBVWN underestimates. These conclusions are in accordance with the behavioral pattern of these functionals for $\bar{\alpha }$ in Cohen et al. (1999) and Vasiliev and Chelikowsky (2010). This verification for the diatomic hydrides goes about as a test-bed for the following arrangement of atoms and molecules. Note that because the converged properties reproduce standard GAMESS results (Schmidt et al., 1993) for all XC functionals (verified for other systems as well), these reference values are discarded hereafter.

TABLE 1

TABLE 1. Static dipole moment μ_z and FF $\mathit{\alpha },\bar{\alpha },\mathit{\beta }$ values (in a.u.) of HCl for different XC functionals. PR implies present result. More details could be found in Ghosal et al. (2018).

3.3.2 Results on $\mathit{\mu },\bar{\alpha },\mathit{\beta }$

Now, in Table 2, we present μ for some linear and non-linear molecules covering both close and open shells of systems extending from diatomic to the hexa-atomic ones at equilibrium geometry. These correspond to the electronic part only; geometry relaxation or vibration contribution has not been incorporated. The total energies are accurately reproduced by the present calculation and not reported here. These can be found in Mandal et al. (2019). The computed zero values of components of μ in the case of non-polar molecules have been correctly produced and henceforth not discussed. The polar molecules show good overall concurrence with experimental results. For a bunch of 29 molecules, the MAD from respective experimental outcomes is 13% (for PBE, LDA) and 10% (for LBVWN, BLYP), separately.

TABLE 2

TABLE 2. Permanent dipole moment of molecules for four XC functionals. All results in a.u. and taken from Mandal et al. (2019).

Next, we advance toward $\bar{\alpha }$ of a cross-section of atoms and molecules, in Table 3. For the sake of comparison, accessible theoretical values from the NIST database and experimental values of zero frequency (containing electronic part only) (Johnson, 2016; Miller and Bederson, 1997) are quoted for comparison. It reveals that each of the three customary functionals (LDA, BLYP, and PBE) overestimate the experimental references. However, LBVWN offers underestimation in all cases (except for Be) and thus differentiates from the other three mentioned functionals.

TABLE 3

TABLE 3. Average static polarizability, $\bar{\alpha }$ for some atoms and molecules using FF KS method, for four XC functionals. All results in a.u. For details, see Mandal et al. (2019).

Now, in this segment, we proceed to higher-order coefficients; Table 4 reports non-zero components of β (using T convention) for some representative molecules. It is evident that the components of a given molecule differ fundamentally from functional to functional—in some cases, including even the sign changes. One such candidate is HI, where β_xxz, β_yyz signs for LDA, PBE functionals are opposite from those of BLYP, PBE. For a comparative understanding, a few selected high-level all-electron calculations (such as CCSD, CAS, and CCSD(T)) in elaborate basis sets (NLO-II, Sadlej, qaug-sadlej, and taug-cc-pVTZ) are provided, along with certain experiments. For clear reasons, our outcomes differ from extended calculations rather significantly. However, this is not the primary objective of this work.

TABLE 4

TABLE 4. The components of β for some selected molecules for four XC functionals. All results are in a.u. See Mandal et al. (2019) for details.

3.3.3 Distorted Geometries

We now offer some illustrative results to explore the efficacy of CCG in determining non-zero components of μ and α, β tensors, as functions of R, in Table 5. As an example, HCl is chosen with R ranging from 1.5 to 3.0 a.u. In general, beyond equilibrium geometry, the static correlation becomes predominant; subsequently, the role of XC functional is of utmost significance. Moreover, the role of the basis set is likewise a major factor. The computed values are in excellent concurrence with theoretical references for all XC functionals throughout the entire region. Upon closer investigation, there is an adjustment in sign in β_xxz=yyz on shifting R from 2.5 to 3 a.u., which is quite satisfactorily captured in InDFT (Roy et al., 2019). Besides that, a comparison with all-electron calculations is additionally performed and portrayed in Figure 1, which are done using the Sadlej basis (Sadlej, 1992) and standard B3LYP functional through the GAMESS program. All functionals reproduce the qualitative shape of α_xx=yy and α_zz very well for the whole range. In both panels, PBE is the nearest to Sadlej-B3LYP results around R_eq. While in panel (a), all plots stay well separated, a distinct crossover is recorded in panel (b) as one moves farther past R_eq. The PBE plot in (b) tends to deviate maximum from all-electron results in the case of all functionals. Consequently, InDFT (Roy et al., 2019) can produce competitive results for α_xx=yy, α_zz, with more elaborate full calculations, both around and away from equilibrium. More point-by-point results and thorough analysis could be found in Mandal et al. (2019).

TABLE 5

TABLE 5. Static dipole moment μ_z, along with FF $\mathit{\alpha },\bar{\alpha },\mathit{\beta }$ (in a.u.) of HCl molecule at various distorted geometries. All quantities are in a.u. More details are available in Ghosal et al. (2018).

FIGURE 1

FIGURE 1. Impact of R on (A) α_xx=yy and (B) α_zz of HCl molecule, taken from Ghosal et al. (2018).

4 HF Exchange Through FCT

While DFT has witnessed an overwhelming number of successful applications, in general, the DFA arouses certain discomfitures. These are 1) piece-wise linearity (PWL) of total energy in fractional particle numbers (Perdew and Wang, 1992; Yang et al., 2000), 2) non-cancellation of fictitious Coulomb self-repulsion energy, often called self-interaction error (SIE) (Perdew and Zunger, 1981; Bao et al., 2018), and 3) asymptotically correct XC potential behavior at LR (Levy et al., 1984). The above three issues are not equivalent, but, to a certain extent, they are interconnected (Kronik and Kümmel, 2020). They provide crucial guidelines in the development of advanced density functionals. A prominent route through which these problems can be addressed is via the introduction of EEX into the picture, which may be combined empirically (Becke, 1993) or non-empirically (Perdew et al., 1996b; Guido et al., 2013), with semi-local functionals. This can improve the asymptotic nature, and, as a consequence, the SIE may get reduced significantly. Following this, the global hybrid functionals (Becke, 1993) (a classic case being B3LYP) were subsequently proposed (Adamo and Barone, 1999; Ernzerhof and Scuseria, 1999; Zhao and Truhlar, 2008); this has tremendously enhanced the chemical applicability of DFT in a large range of chemical systems (Silva-Junior et al., 2008; Mangiatordi et al., 2012). Recently, static correlation in covalent bonds has been treated with general single-determinant model functionals up to the dissociation limit (Becke, 2013; Kong and Proynov, 2015). These emerging hyper-GGA functionals involve exact exchange energy density, e^x, as a fundamental variable, requiring a higher computational cost than the global hybrid ones. Another promising route through which the above conditions can be controlled is optimal tuning of RSH functionals (Baer et al., 2010; Kronik and Kümmel, 2018).

Within the scope of basis-set expansion of MOs, HF exchange energy and related matrix elements can be calculated analytically through four-center electron repulsion integrals (ERI), when GTOs are used. This provides the following contribution to the KS-Fock matrix:

$F_{\mu \nu }^{\text{x}\sigma }=\sum _{\lambda \eta }{P}_{\lambda \eta }^{\sigma }\left(\mu \lambda |\eta \nu \right)=\sum _{\lambda \eta }{P}_{\lambda \eta }^{\sigma }\iint \frac{{\chi }_{\mu ,\sigma }\left(\mathbf{r}\right){\chi }_{\lambda ,\sigma }\left(\mathbf{r}\right){\chi }_{\nu ,\sigma }\left({\mathbf{r}}^{\prime }\right){\chi }_{\eta ,\sigma }\left({\mathbf{r}}^{\prime }\right)}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}d\mathbf{r}d{\mathbf{r}}^{\prime }.(43)$

Here, ERI is represented by (μλ|ην) with μ, ν, λ, η denoting the contracted AO basis, and $P_{\lambda \eta }^{\sigma }$ represents an element of single-particle spin density matrix, P^σ with spin σ. The corresponding HF exchange energy density can be expressed as

$\begin{array}{ccc}\hfill e_{\sigma }^{\text{x}}\left(\mathbf{r}\right)=-\sum _i^{\mathit{occ}}\sum _j^{\mathit{occ}}\int \frac{{\varphi }_{i,\sigma }^{*}\left(\mathbf{r}\right){\varphi }_{j,\sigma }^{*}\left(\mathbf{r}\right){\varphi }_{i,\sigma }\left({\mathbf{r}}^{\prime }\right){\varphi }_{j,\sigma }\left({\mathbf{r}}^{\prime }\right)}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}d{\mathbf{r}}^{\prime }& \hfill & \hfill \\ \hfill =-\sum _{\mu \nu }\sum _{\lambda \eta }{P}_{\mu \nu }^{\sigma }{P}_{\lambda \eta }^{\sigma }\int \frac{{\chi }_{\mu ,\sigma }\left(\mathbf{r}\right){\chi }_{\lambda ,\sigma }\left(\mathbf{r}\right){\chi }_{\nu ,\sigma }\left({\mathbf{r}}^{\prime }\right){\chi }_{\eta ,\sigma }\left({\mathbf{r}}^{\prime }\right)}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}d{\mathbf{r}}^{\prime }.& \hfill & \hfill \end{array}(44)$

The first and second mathematical forms are written in terms of KS-occupied MO (ϕ_σ) and AO (χ_σ). The real form of density matrix and basis gives us the liberty to omit the complex conjugate sign. At first glance, its computational cost appears to be higher than the regular exchange energy calculation (as it requires to be computed at each grid point with four AO indices). Of late, a few proposals have been reported showing considerable computational cost lightening through a pair-atomic RI (Proynov et al., 2010; Liu et al., 2012; Proynov et al., 2012) or a semi-numerical (SNR) (Kossmann and Neese, 2009; Bahmann and Kaupp, 2015; Liu and Kong, 2017; Laqua et al., 2018) approximation.

In what follows, we present a simple, novel strategy for calculating HF exchange density, energy, and necessary matrix elements in CCG, which are significant components for some XC functionals (especially orbital-dependent ones). This takes inspiration from Liu and Kong (2017), where evaluating an intermediate quantity, such as the two-center ESP integral, is a vital step, given as follows:

$v_{\nu \eta }\left(\mathbf{r}\right)=\int \frac{{\chi }_{\nu }\left({\mathbf{r}}^{\prime }\right){\chi }_{\eta }\left({\mathbf{r}}^{\prime }\right)}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}d{\mathbf{r}}^{\prime }=\sum _p\sum _q\int \frac{{\phi }_{\nu }^p\left({\mathbf{r}}^{\prime }\right){\phi }_{\eta }^q\left({\mathbf{r}}^{\prime }\right)}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}d{\mathbf{r}}^{\prime }.(45)$

The two expressions are based on AO basis and primitive functions, ${\phi }_{\nu }^p$ for a particular χ_ν (spin indices omitted for simplicity). We offer a direct, efficient numerical (NR) strategy for the accurate, reliable calculation of ESP integral. This is founded on FCT and employs an RS technique, leading to an LR and SR interaction for CIK. A critical point is the characterization of optimal RS parameter for successful mapping of CIK in CCG from first principles and not empirically. Here, it is achieved through a grid-optimization technique (Ghosal et al., 2016; Ghosal et al., 2018) with respect to the total energy in CCG through a well-defined constraint. This is helpful in the additional development of so-called RSH functionals in combination with the generalized KS theorem (Seidl et al., 1996).

4.1 HF Exchange Energy, Density, and Matrix Elements

This section delineates a numerical methodology for EEX energy and potential, where the former is evaluated by integrating the respective density, $e_{\sigma }^{\text{x}}(\mathbf{r})$, given as

$E_{\sigma }^{\text{x}}=\frac{1}{2}\int e_{\sigma }^{\text{x}}\left(\mathbf{r}\right)d\mathbf{r}.(46)$

Now, Eq. 44 can be rewritten as follows:

$e_{\sigma }^{\text{x}}\left(\mathbf{r}\right)=-\sum _{\nu }{Q}_{\nu }^{\sigma }\left(\mathbf{r}\right)M_{\nu }^{\sigma }\left(\mathbf{r}\right).(47)$

One may anticipate the computation of $e_{\sigma }^{\text{x}}(\mathbf{r})$ in three phases. The first component, $Q_{\nu }^{\sigma }(\mathbf{r})$ can be written as follows:

$Q_{\nu }^{\sigma }\left(\mathbf{r}\right)=\sum _{\mu }{\chi }_{\mu ,\sigma }\left(\mathbf{r}\right)P_{\mu \nu }^{\sigma },(48)$

in which a simple matrix multiplication is used to combine density matrix with the AOs. The step computationally scales as $\mathcal{O}(N_g{N}_B^2)$, where N_B, N_g refer to total number of basis functions and grid points, respectively. The next vital step pertains to the evaluation of ESP integral (contained in $M_{\nu }^{\sigma }(\mathbf{r})$), which as per Eq. 45 also scales the same way. This integral is usually performed analytically (second expression of Eq. 45) using various types of recursion relations such as Obara–Saika (OS) (Obara and Saika, 1986; Obara and Saika, 1988), Head-Gordon–Pople (Head-Gordon and Pople, 1988), or some other combinations (Liu et al., 2016). Here, we have performed this integral through the OS scheme in terms of expensive incomplete Gamma function; this is referred to as the SNR-OS method. Evidently, each ESP integral scales as $\mathcal{O}(N_g{N}_P^2)$, where N_P denotes the average number of primitive functions. The final step requires the computation of $M_{\nu }^{\sigma }(\mathbf{r})$ as follows:

$M_{\nu }^{\sigma }\left(\mathbf{r}\right)=\sum _{\eta }{Q}_{\eta }^{\sigma }\left(\mathbf{r}\right)v_{\nu \eta }^{\sigma }\left(\mathbf{r}\right).(49)$

This segment also scales in the same way as ESP integral calculation, but with lesser steps than the latter; in the innermost loop, only one multiplication and addition are required. This provides a NR route to evaluate the HF exchange matrix, which can be further modified as

$\frac{\partial E_{\sigma }^x}{\partial P_{\lambda \eta }^{\sigma }}=F_{\mu \nu }^{\text{x}\sigma }=-{\int }_{\mathbf{r}}{\chi }_{\mu }\left(\mathbf{r}\right)M_{\nu }^{\sigma }\left(\mathbf{r}\right)d\mathbf{r}.(50)$

Then, the HF exchange energy and KS-Fock matrix with its contribution, can be evaluated accurately in CCG in a purely NR way via the following equations:

$\begin{array}{ccc}\hfill E_{\sigma }^{\text{x}}=\frac{1}{2}{h}_x{h}_y{h}_z\sum _g{e}_{\sigma }^{\text{x}}\left({\mathbf{r}}_g\right),& \hfill & \hfill \\ \hfill F_{\mu \nu }^{\text{x}\sigma }=-h_x{h}_y{h}_z\sum _g{\chi }_{\mu }\left({\mathbf{r}}_g\right)M_{\nu }^{\sigma }\left({\mathbf{r}}_g\right).& \hfill & \hfill \end{array}(51)$

Now, to evaluate ESP integral in CCG, one can rewrite Eq. 45 in the form of

$v_{\nu \eta }\left(\mathbf{r}\right)=\int \frac{{\chi }_{\nu }\left({\mathbf{r}}^{\prime }\right){\chi }_{\eta }\left({\mathbf{r}}^{\prime }\right)}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}d{\mathbf{r}}^{\prime }=\int \frac{{\chi }_{\mu \nu }\left({\mathbf{r}}^{\prime }\right)}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}d{\mathbf{r}}^{\prime }={\chi }_{\nu \eta }\left(\mathbf{r}\right)\star v^c\left(\mathbf{r}\right).(52)$

For simplicity, the spin indices are omitted here. The final mathematical form involves a convolution integral, with χ_νη representing a straightforward multiplication of two AO basis functions, and v^c(r) denoting the usual CIK. This is further simplified by invoking FCT as follows:

$v_{\nu \eta }\left(\mathbf{r}\right)={\mathcal{F}}^{-1}\left\{v^c\left(\mathbf{k}\right){\chi }_{\nu \eta }\left(\mathbf{k}\right)\right\}\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}{\chi }_{\nu \eta }\left(\mathbf{k}\right)=\mathcal{F}\left\{{\chi }_{\nu \eta }\left(\mathbf{r}\right)\right\}.(53)$

Here, v^c(k) and χ_νη(k) signify Fourier integrals of CIK and AO basis functions. The key issue is obtaining a precise mapping of the former, which involves a singularity at r → 0. To deal with this concern, we use a simple RS strategy based on the works of Gill and Adamson (1996) and Gill et al. (1996), expanding the CIK into LR and SR parts using a suitable RS parameter (ζ) as follows:

$\begin{array}{ccc}\hfill v^c\left(\mathbf{r}\right)=\frac{\mathrm{e}\mathrm{r}\mathrm{f}\left(\zeta \mathbf{r}\right)}{\mathbf{r}}+\frac{\mathrm{e}\mathrm{r}\mathrm{f}\mathrm{c}\left(\zeta \mathbf{r}\right)}{\mathbf{r}}& \hfill & \hfill \\ \hfill v^c\left({\mathbf{r}}_g\right)=v_{\mathit{long}}^c\left({\mathbf{r}}_g\right)+v_{\mathit{short}}^c\left({\mathbf{r}}_g\right).& \hfill & \hfill \end{array}(54)$

The last issue is determining an optimum value of parameter ζ_OT, from first principles. This is achieved through the following relation:

${\zeta }_{\text{opt}}\equiv \underset{\zeta }{\mathrm{o}\mathrm{p}\mathrm{t}}{E}_{\sigma }^{\text{x}}\equiv \underset{\zeta }{\mathrm{o}\mathrm{p}\mathrm{t}}{E}_{\text{tot}}=\underset{N_x,N_y,N_z}{\mathrm{o}\mathrm{p}\mathrm{t}}{E}_{\mathit{\text{tot}}},\mathrm{a}\mathrm{t}\mathrm{a}\mathrm{fi}\mathrm{x}\mathrm{e}\mathrm{d}{h}_{\mathbf{r}},(55)$

which is reminiscent of the grid-optimization strategy employed in Section 2.

It is worth noting that, every ESP integral is computed using only a collection of FFTs (two forward and one backward transformation) resulting in a $\mathcal{O}(N_g\log N_g)$ scaling. This contrasts with quadratic scaling with respect to N_P (apart from N_g), in the SNR-OS scheme. Hence, the computational cost is unaffected by the degree of contraction; but in SNR-OS, the cost grows quadratically with the degree of contraction. It is favorable for basis sets with substantial degrees of contraction, needed for a system with decent grid size.

4.1.1 Computational Time of SNR-OS Versus NR

Now, we venture into a comparative discussion on the NR and SNR-OS schemes. Toward this pursuit, a real-time comparison of the performance in terms of the average effective CPU time for an SCF iteration of these two approaches is quoted for a representative set of molecules from Ghosal et al. (2019) in Table 6. Calculations are carried out on a system with Intel Core i7-7700 CPU (3.6 GHz) using the identical optimized grid. A study of the ratio $\left(\frac{SNR-OS}{NR}\right)$ suggests that it hovers in the range of 3.70 (Si₂H₆)–10.83 (CH₄), implying that the NR route offers an advantageous scaling over SNR-OS. This is consistent with the scaling relations given previously. As demonstrated, the ease of implementation of this method is quite encouraging and can be intuited to have a fruitful application in future extensions.

TABLE 6

TABLE 6. Timing (in s) comparison between NR and SNR-OS schemes for one SCF iteration for some representative systems, adopted from Ghosal et al. (2019).

4.1.2 Orbital-Dependent Hybrid Functionals via RS-FCT

The strategy described above is applied in constructing three global hybrid functionals, namely, B3LYP, PBE0, and BHLYP, with a variable amount of former, as well as the traditional XC functional. The XC energy corresponding to each functional is expressed as follows:

$E_{\text{xc}}^{\text{B}3\text{LYP}}=\left(1-a_0\right)E_{\text{LSDA}}^{\text{x}}+a_0{E}_{\text{HF}}^{\text{x}}+a_{\text{x}}{E}_{\text{B}88}^{\text{x}}+a_{\text{c}}{E}_{\text{LYP}}^{\text{c}}+\left(1-a_{\text{c}}\right)E_{\text{VWN}}^{\text{c}},(56)$

$E_{\text{xc}}^{\text{PBE0}}=b_0{E}_{\text{HF}}^{\text{x}}+\left(1-b_0\right)E_{\text{PBE}}^{\text{x}}+E_{\text{PBE}}^{\text{c}},(57)$

$E_{\text{xc}}^{\text{BHLYP}}=\left(1-c_0\right)E_{\text{LSDA}}^{\text{x}}+c_0{E}_{\text{HF}}^{\text{x}}+E_{\text{LYP}}^{\text{c}}.(58)$

Following (Stephens et al., 1994), a₀, a_x, a_c are 0.2, 0.72, 0.81 for B3LYP, whereas in case of PBE0 (Perdew et al., 1996b), b₀ = 0.25. Note that, in PBE0, the contribution of HF exchange is slightly higher than B3LYP, but a higher proportion (c₀ = 0.5) is assigned in BHLYP.

The pertinent LDA- and GGA-type functionals related to B3LYP, BHLYP, and PBE0 are as follows: 1) Vosko–Wilk–Nusair (VWN), the homogeneous electron gas correlation proposed in parametrization formula V (Vosko et al., 1980); 2) B88–incorporating Becke (Becke, 1988) semi-local exchange; 3) Lee–Yang–Parr (LYP) (Lee et al., 1988) semi-local correlation; and 4) Perdew–Burke–Ernzerhof (PBE) (Perdew et al., 1996a) functional for semi-local XC. Other computational details and scaling properties are available in Roy (2008a), Roy (2008b), Roy (2010b), Roy (2011), Ghosal et al. (2016), Ghosal et al. (2018), Mandal et al. (2019).

4.2 Analysis of Hybrid Functionals

Let us begin with the total energies for a few representative atoms and molecules in Table 7. The NR and SNR-OS results are reported for three sets of computations: HF, B3LYP, and PBE0. In all cases, the same pseudopotential, basis set, and convergence (both grid and SCF) criteria of the previous section were engaged. The highest absolute difference in energy (labeled E_diff) between NR and SNR-OS is displayed side by side for simple comparison in all instances. With the exception of O and CH₄, where absolute deviations are far below 0.0004 and 0.001 a.u., the overall agreement between these two sets of results is excellent, showing that the two energies are practically indistinguishable for all species. Needless to say, these energies are in close agreement with those from the standard package GAMESS (Schmidt et al., 1993). For a more detailed analysis, the interested reader may consult (Ghosal et al., 2019).

TABLE 7

TABLE 7. HF, B3LYP, and PBE0 energies (a.u.) of atoms and molecules. E_diff = |E_NR − E_SNR-OS|. These are adopted from Ghosal et al. (2019).

Now, the highest occupied molecular orbital (HOMO) energies are investigated in atoms and molecules for some selected cases in Tables 8 and 9. These are collected from Ghosal et al. (2019), where a broader set of results are available. In addition to the three functionals of the (Table 7), here we also include BHLYP and experimental results. The latter table contains some π-electron molecules (simple, aromatic, and conjugated), where the fundamental gap (difference in energy between HOMO and LUMO) plays an important role. Accurate knowledge of such orbital energies is required for a satisfactory estimation of such gaps. As the outcomes of the NR and SNR-OS schemes are almost identical, we only proceed with the former. A comparison with available theoretical and experimental results reveals that HF HOMO energies (which do not incorporate correlation or relaxation effects) are better than any of the four DFT functionals investigated in terms of agreement with the experiment. This is generally true for a larger data set (Ghosal et al., 2019). Furthermore, it is interesting to note that, with the increase in the fractional contribution of HF exchange (which plays a key role in determining accurate asymptotic behavior at LR) in the hybrid functionals, the deviation falls (e.g., from B3LYP to BHLYP). This could be beneficial in larger systems that require a highly contracted basis, although just a moderate size of grid suffices for the purpose. The precision and ease with which it can be implemented augurs well for its future use in the development of RSH functionals, which may eventually lead to a comprehensive view of HF exchange in the asymptotic limit, bridging theoretical and experimental results.

TABLE 8

TABLE 8. Negative HOMO energies, − ϵ_HOMO (in a.u.) for atoms and molecules using HF, B3LYP, PBE0 XC functionals. For details, consult Ghosal et al. (2019).

TABLE 9

TABLE 9. Negative HOMO energies, − ϵ_HOMO (in a.u.) for selected π-electron molecules using HF, B3LYP, PBE0, and BHLYP XC functionals. Further details are available in Ghosal et al. (2019).

5 OT-RSH Functionals

This section presents an outgrowth of the prior work described earlier. The OT-RSH functionals perform remarkably well in resolving some of the important issues in connection with DFAs (detailed in Section 6). Generally, this is based on a partitioning of CIK into SR and LR parts, using an RS operator, g(γ, r), and an RS parameter, γ, as follows:

$\frac{1}{\mathbf{r}}=\frac{\tilde{g}\left(\gamma ,\mathbf{r}\right)}{\mathbf{r}}+\frac{g\left(\gamma ,\mathbf{r}\right)}{\mathbf{r}},(59)$

where $\tilde{g}(\gamma ,\mathbf{r})$ denotes the complementary RS operator. Historically, this was introduced for the first time (Leininger et al., 1997) in context of multi-resolution CI, where dynamical correlation hardly impacts LR interactions due to fast decaying features. In this scenario, γ plays a central role in adjusting the EEX contribution from SR to LR region for a certain g(γ, r). For a system, usually these two regions are treated separately: the SR region is represented using a revised inter-electronic distance-dependent local/semi-local DFA, while the LR sector, by EEX with g(γ, r)/r rectification. Based on the partitioning scheme, the XC energy can be obtained as follows:

$E_{\text{xc}}=a_{\text{eex}}^{\text{sr}}{E}_{\text{eex}}^{\text{sr}}\left(\gamma \right)+\left(1-a_{\text{eex}}^{\text{sr}}\right)E_{\text{dfa}}^{\text{x,\hspace{0.17em}sr}}\left(\gamma \right)+b_{\text{eex}}^{\text{lr}}{E}_{\text{eex}}^{\text{lr}}\left(\gamma \right)+\left(1-b_{\text{eex}}^{\text{lr}}\right)E_{\text{dfa}}^{\text{x,\hspace{0.17em}lr}}\left(\gamma \right)+E_{\text{dfa}}^{\text{c}},(60)$

where $E_{\text{eex}}^{\text{sr}}$, $E_{\text{eex}}^{\text{lr}}$ signify EEX energy contribution, whereas $E_{\text{dfa}}^{\text{x,\hspace{0.17em}sr}}$, $E_{\text{dfa}}^{\text{x,\hspace{0.17em}lr}}$ represent DFA exchange, at SR and LR segments. A certain set of ($a_{\text{eex}}^{\text{sr}}$, $b_{\text{eex}}^{\text{lr}}$) characterizes a particular mode of partitioning for a given g(γ, r). RSH functionals with $b_{\text{eex}}^{\text{lr}}=1$ offer an asymptotically correct XC potential at LR region. Furthermore, choosing an optimal $a_{\text{eex}}^{\text{sr}}$ strikes a fine balance between EEX and dynamical correlation. As a result, these functionals are not fully free from SIE and also, if not tuned optimally for a desired system, do not follow the PWL condition. This occurs mainly due to a pre-defined default γ, generally obtained semi-empirically by fitting some reference data (Iikura et al., 2001; Tawada et al., 2004; Yanai et al., 2004; Lange et al., 2008).

In OT parlance, γ is usually determined from first principles by imposing Koopmans’ theorem (Salzner and Baer, 2009). It helps satisfy PWL conditions, makes XC potential asymptotically correct at the LR region, and preserves the size-dependency of γ on ρ. Consequently, these functionals improve properties that are rooted in orbital energies, such as vertical ionization energy (IE), fundamental gap, electron affinity (EA), charge-transfer (CT) excitation, optical gap, and Rydberg excitation (Livshits and Baer, 2007; Stein et al., 2010). However, these are hard to maintain with a universal γ (Baer et al., 2010). In recent years, techniques based on electron localization function and localized orbital locator have been attempted, which necessitates only one single SCF calculation (Borpuzari and Kar, 2017; Wang and Zhang, 2018). Also, a self-consistent OT-RSH approach (Tamblyn et al., 2014), based on a minimization of inter-atomic forces, has been reported as well; it produces better geometries and vibrational modes.

Following the general framework of RSH functionals presented in Eq. 60, three distinct types are considered in this rubric which obey a certain well-established mode of partitioning. The first one is the long-range correction (LC) scheme (Iikura et al., 2001) which looks like this

$\left.\begin{array}{c}\hfill E_{\text{xc}}^{\text{LC}}=E_{\text{dfa}}^{\text{x,\hspace{0.17em}sr}}\left(\gamma \right)+E_{\text{eex}}^{\text{lr}}\left(\gamma \right)+E_{\text{dfa}}^{\text{c}},\hfill \\ \hfill g\left(\gamma ,\mathbf{r}\right)=\mathrm{e}\mathrm{r}\mathrm{f}\left(\gamma \mathbf{r}\right)\mathrm{a}\mathrm{n}\mathrm{d}\tilde{g}\left(\gamma ,\mathbf{r}\right)=\mathrm{e}\mathrm{r}\mathrm{f}\mathrm{c}\left(\gamma ,\mathbf{r}\right).\hfill \end{array}\right\}(61)$

The second one is termed as Coulomb-attenuating method (CAM) approach (Yanai et al., 2004), originally introduced utilizing a more general form of g(γ, r) as follows:

$\left.\begin{array}{c}\hfill g_{\alpha ,\beta }\left(\gamma ,\mathbf{r}\right)=\alpha +\beta \mathrm{e}\mathrm{r}\mathrm{f}\left(\gamma \mathbf{r}\right)\mathrm{a}\mathrm{n}\mathrm{d}{\tilde{g}}_{\alpha ,\beta }\left(\gamma ,\mathbf{r}\right)=1-\left[\alpha +\beta \mathrm{e}\mathrm{r}\mathrm{f}\left(\gamma ,\mathbf{r}\right)\right],\hfill \\ \hfill 0\le \alpha +\beta \le 1,0\le \alpha \le 1,\mathrm{a}\mathrm{n}\mathrm{d}0\le \beta \le 1.\hfill \end{array}\right\}(62)$

The α parameter ensures the EEX contribution over the whole range by a factor of α, whereas the β parameter is responsible for the incorporation of DFA throughout the complete range by a factor of 1 − (α + β). In the particular scenario of α = 0, β = 1, the CAM approach gives rise to the previously mentioned LC scheme. These two parameters are related to $a_{\text{eex}}^{\text{sr}}$ in a rather difficult way. The third one, considered here, is called the long-range-corrected (LRC) (Rohrdanz et al., 2009) one, including an additional parameter in $E_{\text{eex}}^{\text{sr}}$ as

$\left.\begin{array}{c}\hfill E_{\text{xc}}^{\text{LRC}}=a_{\text{eex}}^{\text{sr}}{E}_{\text{eex}}^{\text{sr}}\left(\gamma \right)+\left(1-a_{\text{eex}}^{\text{sr}}\right)E_{\text{dfa}}^{\text{x,\hspace{0.17em}sr}}\left(\gamma \right)+E_{\text{eex}}^{\text{lr}}\left(\gamma \right)+E_{\text{dfa}}^{\text{c}},\hfill \\ \hfill g\left(\gamma ,\mathbf{r}\right)=\mathrm{e}\mathrm{r}\mathrm{f}\left(\gamma \mathbf{r}\right)\mathrm{a}\mathrm{n}\mathrm{d}\tilde{g}\left(\gamma ,\mathbf{r}\right)=\mathrm{e}\mathrm{r}\mathrm{f}\mathrm{c}\left(\gamma ,\mathbf{r}\right).\hfill \end{array}\right\}(63)$

This extra parameter $a_{\text{eex}}^{\text{sr}}$ accounts for the incorporation of a certain particular amount of $E_{\text{eex}}^{\text{sr}}$ by a factor $a_{\text{eex}}^{\text{sr}}$. In a special case, when $a_{\text{eex}}^{\text{sr}}=0$, LRC leads to LC. Several other partitioning schemes and different g(γ, r) have been reported in the literature, primarily to deal with accurate thermochemistry and reaction height (Chai and Head-Gordon, 2008a; Chai and Head-Gordon, 2008b; Peverati and Truhlar, 2012b; Vikramaditya et al., 2018; Chan et al., 2019).

To properly incorporate $E_{\text{dfa}}^{\text{x,\hspace{0.17em}sr}}$ in RSH functionals, a number of schemes were proposed, such as adiabatic connection theorem (Baer and Neuhauser, 2005; Cohen et al., 2007), model exchange hole (Iikura et al., 2001; Henderson et al., 2008), and exchange energy density (Chai and Head-Gordon, 2008b; Lin et al., 2012). The present work invokes the framework of Iikura et al. (2001), involving a modified Fermi wave vector in exchange enhancement factor, applicable to any LDA or GGA type DFAs. Later, this was also utilized for CAM-B3LYP (Yanai et al., 2004) through a general g(γ, r), defined in Eq. 62. In this way, the SR GGA-exchange energy can be put forth as

$\left.\begin{array}{c}\hfill E_{\text{gga}}^{\text{x,\hspace{0.17em}sr}}=-\frac{1}{2}\sum _{\sigma }\int {\rho }_{\sigma }^{\frac{4}{3}}{K}_{\text{gga},\sigma }^{\text{x,\hspace{0.17em}sr}}d\mathbf{r},\hfill \\ \hfill K_{\text{gga},\sigma }^{\text{x,\hspace{0.17em}sr}}=K_{\text{gga},\sigma }^{\text{x}}\left[\left(1-\alpha \right)-\beta \left\{\frac{8}{3}{a}_{\sigma }\left[\sqrt{\pi }\mathrm{e}\mathrm{r}\mathrm{f}\left(\frac{1}{2a_{\sigma }}+2a_{\sigma }\left(b_{\sigma }-c_{\sigma }\right)\right)\right]\right\}\right],\hfill \\ \hfill a_{\sigma }=\frac{\gamma }{2K_{\text{gga},\sigma }^{\text{f}}},b_{\sigma }=\exp \left(-\frac{1}{4a_{\sigma }^2}\right),c_{\sigma }=2a_{\sigma }^2{b}_{\sigma }+\frac{1}{2},K_{\text{gga},\sigma }^{\text{f}}={\left(\frac{9\pi }{K_{\text{gga},\sigma }^{\text{x}}}\right)}^{\frac{1}{2}}{\rho }_{\sigma }^{\frac{1}{3}},\hfill \end{array}\right\}(64)$

where $K_{\text{gga},\sigma }^{\text{x}}$ signifies the enhancement factor. The average relative momentum for GGA, $K_{\text{gga},\sigma }^{\text{f}}$, is used to define the modified GGA-enhancement factor, $K_{\text{gga},\sigma }^{\text{x,\hspace{0.17em}sr}}$. One can see that Eq. 64 reproduces the original GGA DFA, when γ = α = 0. The respective potential is evaluated using $K_{\text{gga},\sigma }^{\text{x,\hspace{0.17em}sr}}$, as it was employed for standard GGA DFAs (Johnson et al., 1993). More detailed discussion on SR DFAs can be found in the literature (Iikura et al., 2001; Baer and Neuhauser, 2005; Cohen et al., 2007; Chai and Head-Gordon, 2008b; Henderson et al., 2008; Lin et al., 2012).

Now, we proceed for the discussion on OT-RSH functionals and properties derived from them. Three different kinds of RSH functionals (LC, CAM, and LRC) are employed in our calculations. As in the original articles, the segmentation mode and RS operator remain unaltered. Here, however, γ_OT is determined following the strategy expressed as

${\gamma }_{\text{OT}}\equiv \underset{N_x,N_y,N_z}{\mathrm{o}\mathrm{p}\mathrm{t}}{E}_{\text{tot},\gamma },\mathrm{a}\mathrm{t}\mathrm{fi}\mathrm{x}\mathrm{e}\mathrm{d}{h}_{\mathbf{r}}.(65)$

This optimization technique does not require any fitting scheme. Through the characteristic length of a system, this procedure satisfies the size-dependency principle. In InDFT (Roy et al., 2019), these are implemented for five representative sets of functionals having a varying amount of SR/LR EEX with SR DFA exchange and traditional correlation functional. We consider the LC-BLYP (Tawada et al., 2004) and LC-PBE (Perdew et al., 1996a; Iikura et al., 2001) functionals from the LC-hybrid group with γ = 0.33 and γ = 0.30, respectively. Furthermore, for the CAM-hybrid group, CAM-B3LYP (Yanai et al., 2004) with α = 0.19, β = 0.46, γ = 0.33 and CAM-PBE0 (Lange et al., 2008) with α = 0.25, β = 0.75, γ = 0.30 are utilized. With slight modifications, the original LRC-ωPBEh functional (Rohrdanz et al., 2009) with a^{x, sr} = 0.2, γ = 0.2 is employed for the LRC-hybrid group. Here, it is designated as LRC-ωPBEh^⋆. To distinguish it from the original, it is superscripted with ⋆. The sole difference is about the construction of the SR DFA exchange. All the parameters are left unchanged as in the original paper, except for γ, which is represented by the subscript “ot.” B3LYP and PBE0, the two global hybrid functionals with a configurable quantity of EEX energy and a conventional DFA, are also compared side by side (Stephens et al., 1994; Perdew et al., 1996b). For ease of discussion, the five functionals (LC-BLYP, CAM-B3LYP, LC-PBE0, CAM-PBE0, and LRC-ωPBEh^⋆) are categorized into two distinct blocks: B3LYP (containing B88 exchange and LYP correlation) and PBE0 (PBE exchange and correlation) types.

5.1 Ionization Energy

It is well-known that even if we have the EEX potential, the physical interpretation of KS frontier orbitals is not straightforward; the only exception is the HOMO. The ionization energy of a system can be assigned utilizing a KS analog to Koopmans’ theorem in HF theory and can be written in the form of

$\mathrm{I}\mathrm{E}=-{ϵ}_{\text{HOMO}}.(66)$

When it comes to LDA or GGA-type DFAs, Eq. 66 underestimates HOMO energy. Moreover, this will not work for other functionals outside the KS regime, especially those that interest us. The RSH functionals have, in principle, correct asymptotic behavior in the LR region, but the (G)KS version of Koopmans’ theorem is required to fully capture the essence of HOMO and its energy. It is proved that, for a selected case of an EEX operator, it is still possible to spot the (G)KS HOMO energy with − IE(M) (Görling and Levy, 1997), and, accordingly,

$\mathrm{I}\mathrm{E}=-{ϵ}_{\text{HOMO}}^{\gamma }.(67)$

Like KS mapping, the (G)KS map is not unique. When the RSH functional is considered, any choice of γ can generate a viable approximation of the (G)KS map. The apparent question is whether the (G)KS HOMO energy with − IE(M) can be accurately approximated by the RSH functional with a fixed value of γ, regardless of the system we are interested in. Hence, the comparison of (G)KS HOMO energy with experimental − IE(M) is a good check in determining γ_OT via Eq. 55. For that, the calculated negative (G)KS HOMO energies of a few atoms and molecules for 12 functionals are presented in Table 10. More detailed results are offered in Ghosal and Roy (2022a). It suggests that, for all the functionals, these are close to each other, but a comparison with experimental values indicates an underestimation of obtained results. The fruitfulness of OT-RSH functionals can be probed through a quantity designed as, ϒ = MAE(RSH)/MAE(OT-RSH); this ratio ϒ suggests the reduction in error relative to its unoptimized counterpart (fixed γ). Here, MAE signifies the mean average error. Based on this measure, the five OT functionals from both blocks, for atoms, can be arranged in the following descending order: ${\mathrm{\Upsilon }}_{LR{C}_{\omega }PBE{h}_{ot}},{\mathrm{\Upsilon }}_{CAM-PBE{0}_{ot}},{\mathrm{\Upsilon }}_{LC-PB{E}_{ot}},{\mathrm{\Upsilon }}_{LC-BLY{P}_{ot}},{\mathrm{\Upsilon }}_{CAM-B3LY{P}_{ot}}$. However, for molecules, the arrangement is bit different and the descending order of performance is ${\mathrm{\Upsilon }}_{LR{C}_{\omega }PBE{h}_{ot}},{\mathrm{\Upsilon }}_{LC-PB{E}_{ot}},{\mathrm{\Upsilon }}_{CAM-PBE{0}_{ot}},{\mathrm{\Upsilon }}_{LC-BLY{P}_{ot}},{\mathrm{\Upsilon }}_{CAM-B3LY{P}_{ot}}$. For easy understanding, a subscript “ot” has been added to identify the respective OT functionals. It is seen that, generally, the PBE0 block performs better than B3LYP. Also, within a given block, LC functionals perform better than CAM. Perhaps this is because the auxiliary parameters (α, β, a^x,sr) may have some sensitivity during the self-consistent tuning process, which are kept unaltered. Note that, during optimization of γ_OT, its compatibility with other auxiliary parameters should be taken care of. This is not examined yet and remains a matter for future investigation. In any case, however, the accuracy of RSH functionals is always improved by OT-RSH functionals irrespective of the system or block under consideration.

TABLE 10

TABLE 10. Ionization energies, − ϵ_HOMO for selected atoms and molecules in eV, adopted from Ghosal and Roy (2022a).

5.2 Fundamental Gap

We now proceed to report some properties, which are challenging within DFT, mostly due to the inaccurate description of employed functionals. With that in mind, the estimated (G)KS gap along with the experiment (Johnson, 2016) is tabulated in Table 11 for all functionals for a few illustrative atoms. These are collected from a detailed report available in Ghosal and Roy (2022a). It indicates that, for all functionals, these are comparable to each other except for LC-PBE_ot. Taking the same measure as earlier, the five functionals in descending order of performance are as follows: LRC-ωPBEh^⋆_ot, LC-PBE_ot, LC-BLYP_ot, CAM-PBE0_ot, and CAM-B3LYP_ot. On the contrary, if we compare the MAE, then OT-RSH (LC) functionals seem to do better than OT-RSH (LRC) and OT-RSH (CAM). As found in the previous case, here also, CAM-PBE0 tends to be more accurate than CAM-B3LYP. Thus, these conclusions are in accordance with the earlier findings. For all species, the overall performance of OT-RSH functionals is better than that of RSH. Amongst them, LRC-ωPBEh^⋆_ot and LC-BLYP_ot exhibit excellent performance. Note that the above calculations are done with a pseudo-valence basis and SR LDA/GGA exchange. As a result, there is a substantial prospect of additional improvement, employing an all-electron basis set, including modern SR exchange functionals.

TABLE 11

TABLE 11. (G)KS gap vs. experimental fundamental gap for selected atoms in eV. Results are adopted from Ghosal and Roy (2022a).

5.3 Fractional Charge

Another important measure of success of DFT is a proper description of fractional charge systems. According to the PWL condition (Perdew et al., 1982), for the ground-state energy of systems (of M electrons) with fractional number of electrons (δ), the energy versus δ curve should be a straight line connecting the values at integer. It can be expressed as

$\left.\begin{array}{c}\hfill E_{\text{frac}}\left(N\right)=E\left(N\right)-E_{\text{PWL}}\left(N\right),N=M+\delta ,\hfill \\ \hfill E_{\text{PWL}}\left(N\right)=\left(1-\delta \right)E\left(M\right)+\delta E\left(M+1\right),0\le \delta \le 1,\hfill \\ \hfill E_{\text{PWL}}\left(N\right)=\left(1+\delta \right)E\left(M\right)-\delta E\left(M-1\right),-1\le \delta \le 0,\hfill \end{array}\right\}(68)$

where E(N) and E_PWL(N) define the energy and PWL interpolation of energy for fractional number of electrons, respectively. Here, E_frac is a measure of deviation from PWL behavior. Depending on the choice of range, two different cases arise; when E_frac < 0, the curve is convex and it is concave for E_frac > 0. The well-known DFAs meet certain challenges, resulting in a smooth convex curve, whereas EEX demonstrates reverse trend. RSH functionals which comprise these two ingredients have shown improvement in this direction (Mori-Sánchez et al., 2006).

According to Janak’s theorem (Janak, 1978), for (G)KS HOMO, the change in ϵ_HOMO as a function of fractional occupation number should be straight line. This is then followed by a finite jump at integer point (due to derivative discontinuity). After that again, the variation should be a straight line till the following integer point. With this realization, fractional occupation number n_i, is introduced in density (ignoring degeneracy in HOMO and spin indices) as

$\rho \left(\mathbf{r}\right)=\sum _i{n}_i|{\varphi }_i\left(\mathbf{r}\right){|}^2,n_i\left\{\begin{array}{cc}0,\hfill & i>i_{\mathit{max}}\hfill \\ \delta ,\hfill & i=i_{\mathit{max}}\hfill \\ 1,\hfill & i<i_{\mathit{max}}\hfill \end{array}\right.,(69)$

where i_max corresponds to the HOMO level.

In Figure 2, we illustrate the relative performance of OT-RSH for the C atom as a test case. The pattern in the ranges −1 ≤ δ ≤ 0(0 ≤ N ≤ 1) and 0 ≤ δ ≤ 1(1 ≤ N ≤ 2) is obtained from experimental IE and EA, respectively. Three respective regions in the upper panel are 0, ≤ N ≤ 2 (a), 0, ≤ N ≤ 1 (c), and 1 ≤ N ≤ 2 (e), containing functionals from the B3LYP block, whereas the lower block corresponds to PBE0 functionals, in panels (b), (d), and (f). Once again, in all cases, OT-RSH shows its superiority over the respective RSH functionals. A deeper analysis of panels (c) and (e) reveals that LC-BLYP_ot is quite similar to the straight-line behavior in both regions, and in 1 ≤ N ≤ 2, it performs significantly better. From panels (d) and (f), it follows that LRC-ωPBEh^⋆_ot and CAM-PBE0_ot behave similarly. In region 1 ≤ N ≤ 2, their performance is very close to the experiment with a tiny overall positive shift in energy. When two blocks are compared, OT-RSHs (PBE0) seem to be substantially better than OT-RSH (B3LYP); more precisely, the CAM-PBE0_ot outperforms the CAM-B3LYP_ot. Based on all these facts, the five OT functionals can now be organized in declining order of performance as follows: LRC-ωPBEh^⋆_ot ≈ CAM-PBE0_ot > LC-PBE_ot ≈ LC-BLYP_ot > CAM-B3LYP_ot.

FIGURE 2

FIGURE 2. Performance of various functionals on fractional occupation in C atom. The left panel shows HOMO energy as a fraction of occupied p-electron number for (A) 0, ≤ N ≤ 2, (C) 0, ≤ N ≤ 1, (E) 1 ≤ N ≤ 2, for B3LYP block functionals. Right panels (B,D,F) refer to the PBE0 block, taken from Ghosal and Roy (2022a).

6 Excitation Energies

DFT has feathers in its cap as regards its application to a vast array of electronic systems in ground states. However, it suffers in the case of excited states due to several non-trivial reasons. In today’s time, the most popular approach with reasonable negotiation between accuracy and efficiency is the so-called time-dependent (TD) DFT within the linear response framework, which is a TD variant of KS-DFT. While it is, in principle, an exact method, its success mainly lies in the DFA employed and, specifically, the magnitude to which XC energies are impacted. Despite its numerous successful applications, it faces difficulties in characterizing phenomena such as double excitation, charge transfer, and Rydberg excitation. Here, we offer a simple route for accurate calculation of excited states within a time-independent approach; specifically, we are interested in the ΔSCF method (Ziegler et al., 1977; Kowalczyk et al., 2011). This uses a standard SCF iteration with the non-Aufbau occupation at each iteration using the ground-state functional. From a computational perspective, it is more convenient than linear response TDDFT due to its favorable ground-state-like scaling. Even though it provides a fair enough estimation of excitation energies, it has a tendency to variational collapse. Several sophisticated schemes, such as constrained-DFT (Ziegler et al., 2012; Barca et al., 2014; Ramos and Pavanello, 2016), methods connected to meta dynamics, and gentlest ascent dynamics (Li et al., 2015), have been quite popular in alleviating these issues.

In photo-induced electronic excitation, the lowest singlet excited energies play a crucial role. In principle, due to its multi-determinantal nature, the standard DFT cannot be used on such occasions. However, the calculation of triplet states is quite straightforward. An economical way to compute optical gaps within time-independent DFT in large molecules was introduced lately (Becke, 2016; Becke, 2018a; Becke, 2018b; Becke et al., 2018). For estimating the lowest single-particle excitation energies, the basic model employs a term known as the correlated STS energy. In essence, it involves two independent single-determinant DFT calculations: one for a closed-shell ground state and the other for the lowest triplet state with an open shell. This also requires evaluating a simple 2e⁻ integral (Coulomb self-energy) related to the HOMO-LUMO transition. Interestingly, this strategy is free from concerns involved in standard DFA applied to a triplet excited state, defined by a single Slater determinant and represented by a Fermi hole. As a result, the correlated STS energy ΔE_STS is the major element that could possibly be dealt with using 1) the adiabatic connection theorem (Becke, 2016) and 2) the virial theorem (Becke, 2018a). In a way, this is altogether a non-empirical approach that circumvents the configuration mixing.

We employ these aforementioned approaches to find the excitation energies in molecules by the above approach. The idea is to apply the FCT in conjunction with CIK to deal with the relevant 2e⁻ integral numerically. Also, to analyze its suitability and efficacy, the scheme is adopted to characterize a few properties in molecules of two different genres. First, we consider the organic chromophores, which are extremely significant in nature, showing a photo-luminescence property; a few prominent examples are photosynthesis (Murata, 1969), vision (Palczewski, 2012), and bio-luminescence (Navizet et al., 2011). These materials have wide applications in the development of unique technologies, such as organic light-emitting devices (Mitschke and Bäuerle, 2000; Forrest, 2004), fluorescent sensors (Hou et al., 2015; Li et al., 2016), organic solar cells (Taouali et al., 2018; Khan et al., 2020), medical imaging (Kundu et al., 2009; Li et al., 2016), and laser (Chénais and Forget, 2012). Another class of molecules is charge-transfer (CT) complexes, which form a distinct class of inter- and intra-molecular compounds. These molecules are characterized by the presence of a certain low-energy transition with a relatively strong oscillator strength. According to Mulliken’s quantum theory (Mulliken and Person, 1962; Mulliken and Person, 1969), the ground state of these complexes is typically a resonance hybrid wave function composed of an interacting donor (D) and acceptor (A) that can be expressed as the total of the terms, namely, neutral (DA) and dative (D⁺ A⁻, D⁻ A⁺) states as given in Eq. 70. A partial ground state charge transfer occurs when an electron is transported from donor to acceptor (D⁺ A⁻) and acceptor to donor (D⁻ A⁺). Such complexes have a wave function that looks like (Winget and Brédas, 2011)

$\mathrm{\Psi }\left(\mathrm{D},\mathrm{A}\right)=a\mathrm{\Psi }\left(\mathrm{D}\mathrm{A}\right)+\sum _i{b}_i\mathrm{\Psi }\left({\mathrm{D}}^{+},{\mathrm{A}}^{-}\right)+\sum _i{c}_i\mathrm{\Psi }\left({\mathrm{D}}^{-},{\mathrm{A}}^{+}\right).(70)$

Many processes involving electron-transfers mechanism and molecular conductance rely on these especially characterized excited states.

6.1 Virial Theorem and Adiabatic Connection Theorem for Singlet-Triplet Splitting

Let us consider an excitation of a given system from a closed-shell ground state with an electronic configuration φ_iφ_f. With the assumption of completely filled core with closed shell, this is made up of four Slater determinants: $|{\phi }_i^{\alpha }{\phi }_f^{\alpha }\rangle$, $|{\phi }_i^{\alpha }{\phi }_f^{\beta }\rangle$, $|{\phi }_i^{\beta }{\phi }_f^{\alpha }\rangle$, and $|{\phi }_i^{\beta }{\phi }_f^{\beta }\rangle$, where α, β denote up and down spins. Therefore, diagonalization of Hamiltonian matrix in the vicinity of the above four determinants gives the coupled excited states. The singlet state is given by $|{\psi }_{\text{S}}\rangle =\frac{1}{\sqrt{2}}\{|{\phi }_i^{\alpha }{\phi }_f^{\beta }\rangle -|{\phi }_i^{\beta }{\phi }_f^{\alpha }\rangle \}$, and the three degenerate triplet states, on the contrary, is represented as follows:${\psi }_{\text{T}}=|{\phi }_i^{\alpha }{\phi }_f^{\alpha }\rangle \mathrm{o}\mathrm{r}\frac{1}{\sqrt{2}}\{|{\phi }_i^{\alpha }{\phi }_f^{\beta }\rangle +|{\phi }_i^{\beta }{\phi }_f^{\alpha }\rangle \}\mathrm{o}\mathrm{r}|{\phi }_i^{\beta }{\phi }_f^{\beta }\rangle$. The energies of respective singlets and triplets (denoted by “S” and “T” subscripts) have the form of

$E_{\text{S}}=E^{\alpha \beta }+K_{if}(71)$

and

$E_{\text{T}}=E^{\alpha \alpha }\mathrm{o}\mathrm{r}{E}^{\alpha \beta }-K_{if}\mathrm{o}\mathrm{r}{E}^{\beta \beta },(72)$

where $E^{{\sigma }_1{\sigma }_2}$ is the energy of a given determinant of form $|{\phi }_i^{{\sigma }_1}{\phi }_f^{{\sigma }_2}\rangle$, (σ₁, σ₂) ∈ {α, β}, and K_if is the 2e⁻ integral (Coulomb self-energy of product of transition orbitals) defined as

$K_{if}=\iint \frac{{\phi }_i\left({\mathbf{r}}_1\right){\phi }_f\left({\mathbf{r}}_1\right){\phi }_i\left({\mathbf{r}}_2\right){\phi }_f\left({\mathbf{r}}_2\right)}{|{\mathbf{r}}_1-{\mathbf{r}}_2|}d{\mathbf{r}}_{1}d{\mathbf{r}}_2.(73)$

A combination of Eqs. 71 and 72 gives singlet and triplet excitation energies as follows:

$E_{0\mathrm{S}}=E_{0\mathrm{T}}+2K_{if},(74)$

where E_0S = E_S − E₀, E_0T = E_T − E₀ and ground-state energy of the closed-shell system is denoted by E₀. However, the problem in determining correlated STS energy makes Eq. 74 highly inaccurate for calculation of E_0S. One way to deal with this is to use the well-known “adiabatic connection” theorem (Harris and Jones, 1974), which suggests that, in case of single-particle excitation, the single-triplet energy difference can be expressed in the form of

$\mathrm{\Delta }{E}_{\text{STS}}=\mathrm{\Delta }{E}_{\text{STS}}^0+\mathrm{\Delta }{E}_{\text{ST}}^{\text{corr}},\mathrm{\Delta }{E}_{\text{STS}}=2K_{if}+\mathrm{\Delta }{E}_{\text{ST}}^{\text{corr}}.(75)$

Here, $\mathrm{\Delta }{E}_{\text{STS}}^0$ is the uncorrelated STS energy and $\mathrm{\Delta }{E}_{\text{ST}}^{\text{corr}}$ represents the difference between singlet-triplet correlation energies.

Recently, a non-empirical formula (Becke, 2016) has been proposed to tackle the $\mathrm{\Delta }{E}_{\text{ST}}^{\text{corr}}$ term. This is derived from the inter-electronic cusp condition and the effect it causes to electron correlation. Consequently, it can be written as

$\mathrm{\Delta }{E}_{\text{ST}}^{\text{corr}}=-0.4\int 4{\phi }_i^2\left({\mathbf{r}}_1\right){\phi }_j^2\left({\mathbf{r}}_1\right)z_{\text{C}}^2\left[1-\frac{\ln \left(1+z_{\text{C}}\right)}{z_{\text{C}}}\right]d{\mathbf{r}}_1.(76)$

In this prescription, the only unknown quantity is correlation length z_C. This is nothing but the measure of spatial extent of electron correlation in configuration, φ_iφ_f. In the limit of “strictly correlated electrons,” z_C can be expressed in terms of 2e⁻ integral, K_if, as

$0.4z_{\text{C}}^2\int 4{\phi }_i^2\left({\mathbf{r}}_1\right){\phi }_j^2\left({\mathbf{r}}_1\right)d{\mathbf{r}}_1=2K_{if}.(77)$

Moreover, invocation of standard virial theorem to it (Becke, 2018a) allows one to write

$\mathrm{\Delta }{E}_{\text{ST}}^{\text{corr}}=-K_{if}.(78)$

This surprisingly leads to a further simplification, which reduces the relation to

$\mathrm{\Delta }{E}_{\text{STS}}=2K_{if}-K_{if}=K_{if}.(79)$

It is to be pointed out here that this is a completely non-empirical expression of correlated STS energy but a more simplified one, involving only the 2e⁻ integral.

Therefore, both routes prescribed in Eqs 76 and 79 are entirely non-empirical, implying ΔE_STS to be lower than 2K_if. This leads to a formal definition of a molecule-independent re-scaling parameter f such that

$\mathrm{\Delta }{E}_{\text{STS}}=2f{K}_{if},0<f<1.(80)$

As long as the de-localization error (Perdew et al., 1982) is not a serious concern, determining this parameter in a semi-empirical way might offer overall good quality excitation energies. Optimization of f through a semi-empirical technique (Becke, 2016) gives the value as 0.486 when the results are fitted using the best approximated theoretical excitation energy data set of Silva-Junior et al. (2010). Surprisingly, this is close to the value (0.5) obtained from a consideration of the “virial theorem.”

Here, the central quantity, K_if, is implemented in a manner considerably different from the original prescription (Becke, 2016). In place of the multi-center numerical integral procedure used in Becke and Dickson (1988), we employ an FCT scheme using an RS technique in CIK. As per the description in Ghosal et al. (2019), K_if can be recast as

$K_{if}=\iint \frac{{\phi }_i\left({\mathbf{r}}_1\right){\phi }_f\left({\mathbf{r}}_1\right){\phi }_i\left({\mathbf{r}}_2\right){\phi }_f\left({\mathbf{r}}_2\right)}{|{\mathbf{r}}_1-{\mathbf{r}}_2|}d{\mathbf{r}}_{1}d{\mathbf{r}}_2=\int {\phi }_i\left({\mathbf{r}}_1\right){\phi }_f\left({\mathbf{r}}_1\right)v_{if}\left({\mathbf{r}}_1\right)d{\mathbf{r}}_1.(81)$

Now, the tricky job is to evaluate the v_if integral, for which we employ the RS-FCT procedure once again. This can be further manipulated as

$v_{if}\left({\mathbf{r}}_1\right)=\int \frac{{\phi }_i\left({\mathbf{r}}_2\right){\phi }_f\left({\mathbf{r}}_2\right)}{|{\mathbf{r}}_1-{\mathbf{r}}_2|}d{\mathbf{r}}_2=\int \frac{{\phi }_{if}\left({\mathbf{r}}_2\right)}{|{\mathbf{r}}_1-{\mathbf{r}}_2|}d{\mathbf{r}}_2={\phi }_{if}\left({\mathbf{r}}_1\right)\star v^{\text{c}}\left({\mathbf{r}}_1\right).(82)$

The last expression involves the convolution integral, where φ_if indicates a simple multiplication of ith and fth MO from lowest triplet excited state, whereas v^c(r) signifies the CIK. Further simplification of this integral can be made using FCT as

$v_{if}\left(\mathbf{r}\right)={\mathcal{F}}^{-1}\left\{v^{\text{c}}\left(\mathbf{k}\right){\phi }_{if}\left(\mathbf{k}\right)\right\}\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}{\phi }_{if}\left(\mathbf{k}\right)=\mathcal{F}\left\{{\phi }_{if}\left(\mathbf{r}\right)\right\}.(83)$

Here, v^c(k), φ_if(k) represent Fourier integrals of Coulomb kernel and MOs. Other necessary quantities such as correlation length, z_C, and difference in singlet-triplet correlation energy, are simply evaluated in the same essence, in real-space grid using pseudo KS orbitals φ_i and φ_f.

6.2 Excitation Energies From Becke’s Exciton Model

The effectiveness of the above-mentioned approach is presented through an “SBKJC” type pseudopotential basis set, which is devoid of any diffused function. Furthermore, the lowest singlet and triplet excited states in this exploratory study correspond to the first single excitation for every molecule. Calculations are pursued with an optimized grid using a similar technique to that mentioned in previous sections. In this background, E_0S and E_0T (in eV), as well as correlated STS terms, are tabulated separately in Table 12, along with correlation length (z_C), for B3LYP functional. A cross-section of molecules is presented from a larger set provided in Ghosal et al. (2021); current conclusions are drawn based on that set. Note that E_0S presented here is of both non-empirical and empirical nature: 1) PR₁ uses Eq. 79, which is a semi-empirical approach with re-scaling parameter f = 0.486, 2) PR₂ refers to a non-empirical model from “adiabatic connection” defined in Eq. 75, and 3) PR₃ presents results obtained from Eq. 78 employing a non-empirical model from the virial theorem. In columns 6 and 8, corresponding TD-B3LYP energies (E_0S, E_0T) calculated from GAMESS (Schmidt et al., 1993) using the same functional and basis set are presented for side-by-side comparison. An analysis of E_0S suggests that PR₃ improves results from PR₂. The performance of PR₁ is in close proximity to PR₃. Therefore, out of three methods, PR₃ proves to be the best estimate as those are quite competitive with TD-B3LYP. The next columns provide E_0T and ΔE_STS for an understanding of the contribution of these terms in calculating E_0S. Though the overall agreement is satisfactory for both E_0S and E_0T, the worst performance is observed for propene (relative to TD-B3LYP). A careful analysis suggests that the major source of inaccuracy in E_0S is the STS term, not E_0T. This leads to the understanding that the success of our method depends on an accurate estimation of 2e⁻ integrals or, in other words, the accuracy of triplet states.

TABLE 12

TABLE 12. E_0S, E_0T, and ΔE_STS (in eV) using B3LYP functional. For details, see Ghosal et al. (2021).

In order to validate the usefulness of the proposed method and enlarge the scope of applicability, some larger molecular systems are approached in Table 13. Thus, E_0S, E_0T for some representative organic chromophores and linear acenes (Ghosal et al., 2021) are presented. Geometries of these systems are taken from Silva-Junior et al. (2010), whereas the same for linear acenes are obtained from GAMESS calculation using the B3LYP functional and CC-pVDZ basis. In addition to B3LYP, LC-BLYP from the family of RSH functionals is also invoked. All the results henceforth refer to the PR₃ calculation. Apparently, the performance of B3LYP is much more consistent than that of LC-BLYP; the former shows an overestimation in excitation energies. Also, instead of a dramatic one, only a subtle betterment of results is observed using LC-BLYP. The effect of full HF exchange at LR has no dramatic effect on excitation energies, although it enhances the behavior of frontier orbitals used in K_if computations. This discordance has occurred as γ is assumed to be independent of system size. It is possible to achieve a greater level of performance by treating γ as a system-dependent parameter (functional of ρ) estimated from the first principles (Baer et al., 2010). It is believed that an optimally tuned (in the spirit of the size-dependency principle) γ will outperform the conventional hybrid and RSH functionals in terms of results. An analogous qualitative trend is also depicted by linear acenes.

TABLE 13

TABLE 13. Excitation energies of organic chromophores and linear acenes from “virial theorem.” These are taken from Ghosal et al. (2021).

6.2.1 Basis-Set Dependence

This section produces optical gaps generated from n → π* and π → π* transitions in selected molecules in Table 14. In order to probe the dependence on the basis set, singlet excitation energies are reported employing two basis sets, namely, 6-311G (B1) and 6-311 + +G* (B2), both with B3LYP functional. The corresponding geometries are taken from supplementary materials of Silva-Junior et al. (2010). All calculations are done with all-electron orbitals. The energies E₀ and E_T are calculated using the GAMESS program package. The triplet calculations correspond to restricted open-shell. The K_if integrals are evaluated numerically with our InDFT program (Roy et al., 2019), taking orbitals from GAMESS. The results are compared with the “theoretical best estimate” TBE-2 (Silva-Junior et al., 2010) benchmark values tabulated in column 3 and are comparable. A detailed analysis in terms of MAE and ME values for a larger set is available in Roy et al. (2021), of which Table 15 is a subset. Now, similar to our previous conclusion, these results are also comparable with TD-B3LYP energies; in fact, the ones with B2 are in better agreement (Roy et al., 2021). This reflects the basis set dependency of these quantities. It is also delineated by Roy et al. (2021) that TD-B3LYP significantly underestimates the excitation energies from the current procedure.

TABLE 14

TABLE 14. E_0S (in eV) in organic dyes, using B3LYP functional. See Roy et al. (2021) for details.

6.2.2 Dependence on XC Functional

Next, in Figure 3, we investigate the functional reliance of the present approach. In this regard, in left and right panels, E_0S for polyenes and linear acenes with three different functionals, namely, B3LYP, ωB97X, and CAM-B3LYP, are plotted with respect to polyene length and the number of rings. These three functionals have a variable amount of EEX contribution throughout the inter-electronic distance. This study aims to demonstrate the role of EEX in determining the optical gap in systems with reduced HOMO-LUMO gap. When the chain length (or rings) increases, the gap squeezes. As one moves from panels (a) to (b), the optical gap exhibits a similar pattern. This feature is well reproduced by all three functionals, but when compared to TBE-2 (Silva-Junior et al., 2010), B3LYP outperforms all of them. Specifically, ωB97X deviates most from reference, followed by CAM-B3LYP. The impact of EEX in the LR region is less important; rather, an appropriate balance between EEX and correlation is required throughout the inter-electronic distances. This is consistent with the fact that, unlike Rydberg and CT excitations, the optical gap has no LR feature. According to Becke (2018b), the optimal contribution of EEX is around 21%; thus, the success of B3LYP over CAM-B3LYP and ωB97X is easily understandable. When comparing CAM-B3LYP with ωB97X, the former wins as EEX’s contribution in SR is still minor when compared to ωB97X. However, the system-independent value of the RS parameter in ωB97X and CAM-B3LYP fails to induce the size dependency of the optical gap problem, implying that this size-dependency component is necessary in the parameter. Overall, the present method shows sensitivity toward EEX contributions in SR (rather LR).

FIGURE 3

FIGURE 3. E_0S with different functionals against (A) polyene length and (B) the number of rings. Both panels employ a B1 basis. More details are given by Roy et al. (2021).

6.2.3 Optical Gaps in Organic Chromophores

Now, we explore a few applications pertaining to the photoluminescence effects in organic chromophores. This requires a detailed account of low-lying excited states. In this regard, the HOMO-LUMO gap (E_L−H) and E_0S for a few representative linear and nonlinear poly-cyclic aromatic hydrocarbons (PAH) and organic dyes (comparatively difficult systems) from Roy et al. (2021) are tabulated in Table 15. Geometries for the organic dyes are supplied by Kowalczyk et al. (2011). These correspond to B3LYP/cc-PVTZ calculations. Along with the calculated excitation energies, experimental and TD-B3LYP results are also presented here. The consistency in overall performance is quite encouraging. A careful analysis (Roy et al., 2021) shows that PR is comparable with TD-B3LYP for both PAHs and organic dyes. While excitation energies calculated with TD-B3LYP shows systematic underestimation, the present scheme, in contrast, shows an overestimation consistently, which follows Becke (2018a). For organic dyes, while the efficiency of TD-B3LYP is superior to the current approach, the inaccuracy is more consistent in the latter. Note that the energies obtained for organic dyes have shown sensitivity toward the chosen basis set.

TABLE 15

TABLE 15. HOMO-LUMO gap (E_L-H), HOMO-LUMO singlet excitation energy in organic chromophores. PR ≡ present result. Details are available in Roy et al. (2021).

6.2.4 Charge-Transfer Excitation Within a Hybrid (G)KS Framework

The working equation on this occasion is Eq. 79. The respective single-point energies (E₀ and E_T) are obtained from conventional DFT calculations using Gaussian 09 package (Frisch et al., 2016). Full calculations are performed with the cc-PVTZ basis set employing BLYP, B3LYP, and LC-BLYP functionals. Similar to the previous case, K_if integrals are accomplished using InDFT (Roy et al., 2019). Two different kinds of CT complexes are chosen, which are relatively hard to deal with numerically. First, we consider weakly bound CT complexes (first four molecules of CT-A of Table 16). They are bounded with non-covalent bonds of length in the range of 3.2–3.6 Å, making them quite attractive and challenging. A set of intra-molecular CT complexes is also studied (last three of CT-A). Finally, we look at some harder organic compounds (presented in column CT-B) that display thermally activated delayed fluorescence (TADF). The geometry of these molecules is sourced from the supplementary materials of Hait et al. (2016). In Table 16, excitation energies, PR_n, calculated from the current scheme and TD_n referring to the same with TDDFT, are presented. Here, subscript “n” symbolizes three XC functionals; n = 1, 2, 3 stands for BLYP, B3LYP, and LC-BLYP. For comparison, literature/experimental results are also provided side by side. Now, energies from restricted open shell or unrestricted calculations are reported as per convergence; we have taken this liberty as the difference between these two calculations is not significant enough to bring any perceptible change in excitation energy (Roy et al., 2022). As we shift from BLYP to LC-BLYP via B3LYP for both PR and TD, a general trend of increment in excitation energy is noticed in both CT complexes. An in-depth examination indicates that, for molecules in CT-A, the impact of K_if gets dominated by E_0T for all the XC functionals considered here. Thus, in the framework, to a large extent, E_0T determines the observed trend in E_0S. However, for systems in CT-B, there is no way to determine a priori which will be a dominant factor between E_0T and K_if for a given case. Therefore, for systems with E_0T dominance over K_if, B3LYP performs better than RSH functionals, but when their contributions are comparable, the mutual cancellation of these determines the excitation energy.

TABLE 16

TABLE 16. E_0S from the present result (PR) and TDDFT (TD) in some CT complexes. NC implies “not converged.” More details are available in Roy et al. (2022).

Finally, we proceed to investigate the asymptotic limit of CT excitation in some weakly interacting systems, characterized by R⁻¹ energy decay (R is inter-molecular separation). This study also intends to exploit the perseverance of the present scheme throughout R. With this in mind, in Figure 4, the excitation energy of an inter-molecular dimer C₂H₄ − C₂F₄ is depicted as a function of R for functionals that have already been mentioned. Results are compared with configuration interaction singles (CIS) (Foresman et al., 1992) which may be a benchmark. For the entire region, all the PR_n(n = 1–3) energies are in admissible accordance with CIS, without any substantial difference between them. However, this is not true for TD_n(n = 1–3) results. They remain distinctly separated. TD₃ offers overestimation from CIS, whereas the other two record underestimation. Note that the B3LYP result (Feng et al., 2018) within the “virial theorem” model slightly overestimates CIS for the entire region of R, whereas PR₂ energies are underestimated. This lean disparity perhaps is introduced from two separate numerical strategies used for the calculation of K_if. In contrast to the present approach, TDDFT appears to be more sensitive to R and shows CT breakdown. Surprisingly, TD₃ appears to provide the poorest result. This is mainly due to the system independence of γ. The invocation of a size-dependent γ possibly can partly alleviate this problem.

FIGURE 4

FIGURE 4. E_0S for C₂ H₄ − C₂ F₄ versus R. PR₁, PR₂, and PR₃ denote the present result with BLYP, B3LYP, and LC-BLYP functional, and TD₁, TD₂, and TD₃ represent the same within TDDFT. Adapted from Roy et al. (2022).

7 Future and Outlook

We have exhibited the legitimacy and viability of the Gaussian-based LCAO-MO approach to DFT involving CCG in the context of atomic and molecular properties. A wide variety of atoms and molecules are used as test beds to examine the efficacy of CCG in the context of μ, α, β with an optimized FF procedure. Comparison with existing theoretical and experimental data vouch for its suitability and effectiveness. The feasibility and practicability of a direct NR approach for coupling accurate exchange energy density, which is a key component of hyper-type of functionals, energy, and matrix in real-space CCG, is discussed. This was done for a variety of atoms and molecules, with properties including total energy and orbital energies. These were also shown for B3LYP, PBE0, and BHLYP hybrid functionals. The effectiveness of this strategy is reliant on the precise estimation of ESP integral, which in turn, depends on the optimization of the RS parameter in CIK. The scaling suggests that this approach could be very useful in massive large-scale DFT calculations incorporating orbital-dependent functionals. Most importantly, within the CCG, the NR scheme appeared to be more proficient than the alternative SNR-OS technique. Carrying forward the success of this approach, in the following segment, the suitability and performance of a self-consistent systematic optimization procedure for OT-RSH functionals are presented. Their performance was assessed by probing properties such as frontier orbital energies, fractional occupation of electron on HOMO, and fundamental gap. For finite systems, the predominance of OT-RSH functionals is observed over the respective RSH functionals. The success of this method relies on the precise estimation of γ_OT based on the size-dependency principle.

We then move to the realm of excited states (within a time-independent scheme) and detail the practicability and convenience of a simple yet accurate TIKS-DFT approach, to calculate single excitation energies in a realistic manner. This was tested for a host of linear acenes having π network, organic chromophores, linear acenes, and charge-transfer complexes. The derived results from the virial theorem are in appreciable accordance with reference results for all species. This simple scheme has been exhibited as a feasible choice for predicting optical gaps in organic chromophores. The above-mentioned outcome of the CCG method encourages us to use more effective core potential, more elaborate and sophisticated basis sets, and superior quality density functionals (e.g., RSH, hyper, and local hybrid XC functionals) to study different physicochemical aspects of many-electron systems. It might likewise be alluring to inspect its performance in a variety of exciting configurations apart from the lowest excited state. This approach could be highly beneficial in real-time dynamical approach, especially laser-molecules interactions in the intense domain within the TDDFT framework. In this pursuit, early results have very recently got published in Ghosal and Roy (2022b). In this article, the application of an intense laser field on electron dynamics of H₂ and N₂ molecules has been performed using real-time TDDFT. Moreover, other than single-point calculations, it would also be interesting to assess the merit and suitability of the current method for geometry optimization of molecules in CCG. A significant concern would be to reduce the computation cost by invoking a linear scaling technique. Some of these works are presently being scrutinized.

Author Contributions

The main idea of this work was initiated by AR. The core program was written by AR. The manuscript is written by SM and AR.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We are grateful to the Editors TH and Sam P De Visser for the kind invitation for this special issue. We also sincerely thank the Editorial office for keeping our request to extend the deadline on more than one occasion. Financial support from DST SERB, New Delhi, India (sanction order CRG/2019/000293) is gratefully acknowledged. Partial funding from MATRICS, DST, New Delhi (sanction order MTR/2019/000012) and CSIR, New Delhi [sanction order: 01(3027)/21/EMR-II] is also deeply appreciated.

References

Adamo, C., and Barone, V. (1999). Toward Reliable Density Functional Methods without Adjustable Parameters: The PBE0 Model. J. Chem. Phys. 110, 6158–6170. doi:10.1063/1.478522

CrossRef Full Text | Google Scholar

A. K. Roy (Editor) (2012). Theoretical and Computational Developments in Modern Density Functional Theory (New York: Nova Science Publishers).

Google Scholar

Baer, R., Livshits, E., and Salzner, U. (2010). Tuned Range-Separated Hybrids in Density Functional Theory. Annu. Rev. Phys. Chem. 61, 85–109. doi:10.1146/annurev.physchem.012809.103321

PubMed Abstract | CrossRef Full Text | Google Scholar

Baer, R., and Neuhauser, D. (2005). Density Functional Theory with Correct Long-Range Asymptotic Behavior. Phys. Rev. Lett. 94, 043002. doi:10.1103/physrevlett.94.043002

PubMed Abstract | CrossRef Full Text | Google Scholar

Bahmann, H., and Kaupp, M. (2015). Efficient Self-Consistent Implementation of Local Hybrid Functionals. J. Chem. Theory Comput. 11, 1540–1548. doi:10.1021/ct501137x

PubMed Abstract | CrossRef Full Text | Google Scholar

Bao, J. L., Gagliardi, L., and Truhlar, D. G. (2018). Self-Interaction Error in Density Functional Theory: An Appraisal. J. Phys. Chem. Lett. 9, 2353–2358. doi:10.1021/acs.jpclett.8b00242

PubMed Abstract | CrossRef Full Text | Google Scholar

Barca, G. M. J., Gilbert, A. T. B., and Gill, P. M. W. (2014). Communication: Hartree-Fock Description of Excited States of H2. J. Chem. Phys. 141, 111104. doi:10.1063/1.4896182

PubMed Abstract | CrossRef Full Text | Google Scholar

Beatson, R., and Greengard, L. (1997). “A Short Course on Fast Multipole Methods,” in Numerical Mathematics and Scientific Computation (Oxford, United Kingdom: Oxford University Press), 1–37.

Google Scholar

Becke, A. D. (1993). A New Mixing of Hartree-Fock and Local Density‐functional Theories. J. Chem. Phys. 98, 1372–1377. doi:10.1063/1.464304

CrossRef Full Text | Google Scholar

Becke, A. D. (2018b). Communication: Optical Gap in Polyacetylene from a Simple Quantum Chemistry Exciton Model. J. Chem. Phys. 149, 081102. doi:10.1063/1.5050615

CrossRef Full Text | Google Scholar

Becke, A. D., Dale, S. G., and Johnson, E. R. (2018). Communication: Correct Charge Transfer in CT Complexes from the Becke'05 Density Functional. J. Chem. Phys. 148, 211101. doi:10.1063/1.5039742

PubMed Abstract | CrossRef Full Text | Google Scholar

Becke, A. D. (2013). Density Functionals for Static, Dynamical, and Strong Correlation. J. Chem. Phys. 138, 074109. doi:10.1063/1.4790598

PubMed Abstract | CrossRef Full Text | Google Scholar

Becke, A. D. (1988). Density-functional Exchange-Energy Approximation with Correct Asymptotic Behavior. Phys. Rev. A 38, 3098–3100. doi:10.1103/physreva.38.3098

PubMed Abstract | CrossRef Full Text | Google Scholar

Becke, A. D., and Dickson, R. M. (1988). Numerical Solution of Poisson's Equation in Polyatomic Molecules. J. Chem. Phys. 89, 2993–2997. doi:10.1063/1.455005

CrossRef Full Text | Google Scholar

Becke, A. D. (2018a). Singlet-triplet Splittings from the Virial Theorem and Single-Particle Excitation Energies. J. Chem. Phys. 148, 044112. doi:10.1063/1.5012033

PubMed Abstract | CrossRef Full Text | Google Scholar

Becke, A. D. (2016). Vertical Excitation Energies from the Adiabatic Connection. J. Chem. Phys. 145, 194107. doi:10.1063/1.4967813

PubMed Abstract | CrossRef Full Text | Google Scholar

Bishop, D. M. (1994). Explicit Nondivergent Formulas for Atomic and Molecular Dynamic Hyperpolarizabilities. J. Chem. Phys. 100, 6535–6542. doi:10.1063/1.467062

CrossRef Full Text | Google Scholar

Bishop, D. M., and Norman, P. (1999). Effects of Vibration on the Polarizability and the First and Second Hyperpolarizabilities of HF, HCl, and HBr. J. Chem. Phys. 111, 3042–3050. doi:10.1063/1.479661

CrossRef Full Text | Google Scholar

Bishop, D. M., and Pipin, J. (1987). Field and Field-Gradient Polarizabilities of H2O. Theor. Chim. Acta 71, 247–253. doi:10.1007/bf00529096

CrossRef Full Text | Google Scholar

Boese, A. D., and Handy, N. C. (2001). A New Parametrization of Exchange-Correlation Generalized Gradient Approximation Functionals. J. Chem. Phys. 114, 5497–5503. doi:10.1063/1.1347371

CrossRef Full Text | Google Scholar

Borpuzari, M. P., and Kar, R. (2017). A New Nonempirical Tuning Scheme with Single Self-Consistent Field Calculation: Comparison with Global and IP-Tuned Range-Separated Functional. J. Comput. Chem. 38, 2258–2267. doi:10.1002/jcc.24876

PubMed Abstract | CrossRef Full Text | Google Scholar

Brandt, A. (1977). Multi-level Adaptive Solutions to Boundary-Value Problems. Math. Comp. 31, 333–390. doi:10.1090/s0025-5718-1977-0431719-x

CrossRef Full Text | Google Scholar

Brémond, É., Ottochian, A., Pérez‐Jiménez, Á. J., Ciofini, I., Scalmani, G., Frisch, M. J., et al. (2021). Assessing Challenging Intra‐ and Inter‐molecular Charge‐transfer Excitations Energies with Double‐hybrid Density Functionals. J. Comput. Chem. 42, 970–981. doi:10.1002/jcc.26517

PubMed Abstract | CrossRef Full Text | Google Scholar

Buckingham, A. D. (1967). in Advances in Chemical Physics. Editor J. O. Hirschfelder (New York: Interscience).

Google Scholar

Bulat, F. A., Toro-Labbé, A., Champagne, B., Kirtman, B., and Yang, W. (2005). Density-functional Theory (Hyper)polarizabilities of Push-Pull π-conjugated Systems: Treatment of Exact Exchange and Role of Correlation. J. Chem. Phys. 123, 014319. doi:10.1063/1.1926275

PubMed Abstract | CrossRef Full Text | Google Scholar

Burke, K. (2012). Perspective on Density Functional Theory. J. Chem. Phys. 136, 150901. doi:10.1063/1.4704546

PubMed Abstract | CrossRef Full Text | Google Scholar

Carmon-Espíndola, J., Flores-Moreno, R., and Köster, A. M. (2012). Static and Dynamic First Hyperpolarizabilities from Time-dependent Auxiliary Density Perturbation Theory. Int. J. Quant. Chem. 112, 3461–3471.

Google Scholar

Castet, F., and Champagne, B. (2012). Assessment of DFT Exchange-Correlation Functionals for Evaluating the Multipolar Contributions to the Quadratic Nonlinear Optical Responses of Small Reference Molecules. J. Chem. Theory Comput. 8, 2044–2052. doi:10.1021/ct300174z

PubMed Abstract | CrossRef Full Text | Google Scholar

Chai, J.-D., and Head-Gordon, M. (2008). Optimal Operators for Hartree-Fock Exchange from Long-Range Corrected Hybrid Density Functionals. Chem. Phys. Lett. 467, 176–178. doi:10.1016/j.cplett.2008.10.070

CrossRef Full Text | Google Scholar

Chai, J.-D., and Head-Gordon, M. (2008). Systematic Optimization of Long-Range Corrected Hybrid Density Functionals. J. Chem. Phys. 128, 084106. doi:10.1063/1.2834918

PubMed Abstract | CrossRef Full Text | Google Scholar

Champagne, B. (2009). in Specialist Periodical Reports: Chemical Modelling, Applications and Theory. Editor M. Springborg (London: Royal Society of Chemistry), Vol. 6.

Google Scholar

Chan, B., Kawashima, Y., and Hirao, K. (2019). The reHISS Three-Range Exchange Functional with an Optimal Variation of Hartree-Fock and its Use in the reHISSB-D Density Functional Theory Method. J. Comput. Chem. 40, 29–38. doi:10.1002/jcc.25383

PubMed Abstract | CrossRef Full Text | Google Scholar

Chang, C.-M., Shao, Y., and Kong, J. (2012). Ewald Mesh Method for Quantum Mechanical Calculations. J. Chem. Phys. 136, 114112. doi:10.1063/1.3694829

PubMed Abstract | CrossRef Full Text | Google Scholar

Chénais, S., and Forget, S. (2012). Recent Advances in Solid-State Organic Lasers. Polym. Int. 61, 390–406. doi:10.1002/pi.3173

CrossRef Full Text | Google Scholar

Cohen, A. J., Handy, N. C., and Tozer, D. J. (1999). Density Functional Calculations of the Hyperpolarisabilities of Small Molecules. Chem. Phys. Lett. 303, 391–398. doi:10.1016/s0009-2614(99)00248-1

CrossRef Full Text | Google Scholar

Cohen, A. J., Mori-Sánchez, P., and Yang, W. (2007). Development of Exchange-Correlation Functionals with Minimal Many-Electron Self-Interaction Error. J. Chem. Phys. 126, 191109. doi:10.1063/1.2741248

PubMed Abstract | CrossRef Full Text | Google Scholar

Colwell, S. M., Murray, C. W., Handy, N. C., and Amos, R. D. (1993). The Determination of Hyperpolarisabilities Using Density Functional Theory. Chem. Phys. Lett. 210, 261–268. doi:10.1016/0009-2614(93)89131-z

CrossRef Full Text | Google Scholar

Cook, M., and Karplus, M. (1987). Electron Correlation and Density-Functional Methods. J. Phys. Chem. 91, 31–37. doi:10.1021/j100285a010

CrossRef Full Text | Google Scholar

Dale, S. G., Johnson, E. R., and Becke, A. D. (2017). Interrogating the Becke'05 Density Functional for Non-locality Information. J. Chem. Phys. 147, 154103. doi:10.1063/1.5000909

PubMed Abstract | CrossRef Full Text | Google Scholar

de Wergifosse, M., Castet, F., and Champagne, B. (2015). Frequency Dispersion of the First Hyperpolarizabilities of Reference Molecules for Nonlinear Optics. J. Chem. Phys. 142, 194102. doi:10.1063/1.4920977

PubMed Abstract | CrossRef Full Text | Google Scholar

de Wergifosse, M., Liégeois, V., and Champagne, B. (2014). Evaluation of the Molecular Static and Dynamic First Hyperpolarizabilities. Int. J. Quantum Chem. 114, 900–910. doi:10.1002/qua.24685

CrossRef Full Text | Google Scholar

Dirac, P. A. M. (1930). Note on Exchange Phenomena in the Thomas Atom. Math. Proc. Camb. Phil. Soc. 26, 376–385. doi:10.1017/s0305004100016108

CrossRef Full Text | Google Scholar

Dirac, P. A. M. (1929). Quantum Mechanics of Many-Electron Systems. Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 123, 714. doi:10.1098/rspa.1929.0094

CrossRef Full Text | Google Scholar

Dudley, J. W., and Ward, J. F. (1985). Measurements of Second‐ and Third‐order Nonlinear Polarizabilities for HF and HCl. J. Chem. Phys. 82, 4673–4677. doi:10.1063/1.448726

CrossRef Full Text | Google Scholar

Dunlap, B. I., Connolly, J. W. D., and Sabin, J. R. (1979b). On First-Row Diatomic Molecules and Local Density Models. J. Chem. Phys. 71, 4993. doi:10.1063/1.438313

CrossRef Full Text | Google Scholar

Dunlap, B. I., Connolly, J. W. D., and Sabin, J. R. (1979a). On Some Approximations in Applications ofXα Theory. J. Chem. Phys. 71, 3396–3402. doi:10.1063/1.438728

CrossRef Full Text | Google Scholar

El Ghazaly, M. O. A., Svendsen, A., Bluhme, H., Nielsen, S. B., and Andersen, L. H. (2005). Electron Scattering on P-Benzoquinone Anions. Chem. Phys. Lett. 405, 278–281. doi:10.1016/j.cplett.2005.02.048

CrossRef Full Text | Google Scholar

Ernzerhof, M., and Scuseria, G. E. (1999). Assessment of the Perdew-Burke-Ernzerhof Exchange-Correlation Functional. J. Chem. Phys. 110, 5029–5036. doi:10.1063/1.478401

CrossRef Full Text | Google Scholar

Feller, D. (1996). The Role of Databases in Support of Computational Chemistry Calculations. J. Comput. Chem. 17, 1571–1586. doi:10.1002/(sici)1096-987x(199610)17:13<1571::aid-jcc9>3.0.co;2-p

CrossRef Full Text | Google Scholar

Feng, X., Becke, A. D., and Johnson, E. R. (2018). Communication: Becke's Virial Exciton Model Gives Accurate Charge-Transfer Excitation Energies. J. Chem. Phys. 149, 231101. doi:10.1063/1.5078515

PubMed Abstract | CrossRef Full Text | Google Scholar

Fernández, B., Coriani, S., and Rizzo, A. (1998). MCSCF Polarizability and Hyperpolarizabilities of HCl and HBr. Chem. Phys. Lett. 288, 677–688. doi:10.1016/s0009-2614(98)00355-8

CrossRef Full Text | Google Scholar

Feynman, R. P. (1939). Forces in Molecules. Phys. Rev. 56, 340–343. doi:10.1103/physrev.56.340

CrossRef Full Text | Google Scholar

Flores-Moreno, R., and Köster, A. M. (2008). Auxiliary Density Perturbation Theory. J. Chem. Phys. 128, 134105. doi:10.1063/1.2842103

PubMed Abstract | CrossRef Full Text | Google Scholar

Foresman, J. B., Head-Gordon, M., Pople, J. A., and Frisch, M. J. (1992). Toward a Systematic Molecular Orbital Theory for Excited States. J. Phys. Chem. 96, 135–149. doi:10.1021/j100180a030

CrossRef Full Text | Google Scholar

Forrest, S. R. (2004). The Path to Ubiquitous and Low-Cost Organic Electronic Appliances on Plastic. Nature 428, 911–918. doi:10.1038/nature02498

PubMed Abstract | CrossRef Full Text | Google Scholar

Fournier, R. (1990). Second and Third Derivatives of the Linear Combination of Gaussian Type Orbitals-Local Spin Density Energy. J. Chem. Phys. 92, 5422–5429. doi:10.1063/1.458520

CrossRef Full Text | Google Scholar

Frigo, M., and Johnson, S. G. (2005). The Design and Implementation of FFTW3. Proc. IEEE 93, 216–231. doi:10.1109/jproc.2004.840301

CrossRef Full Text | Google Scholar

Frisch, M. J., Trucks, G. W., Schlegel, H. B., Scuseria, G. E., Robb, M. A., Cheeseman, J. R., et al. (2016). Gaussian 16 Revision C.01. Wallingford CT: gaussian Inc.

Google Scholar

George, M. (2006). Atoms, Molecules and Clusters in Electric Fields: Theoretical Approaches to the Calculation of Electric Polarizability, 1. World Scientific. doi:10.1142/p464

CrossRef Full Text | Google Scholar

Ghosal, A., Gupta, T., Mahato, K., and Roy, A. K. (2021). Excitation Energies through Becke's Exciton Model within a Cartesian-Grid KS DFT. Theor. Chem. Acc. 140, 2. doi:10.1007/s00214-020-02699-5

CrossRef Full Text | Google Scholar

Ghosal, A., Mandal, T., and Roy, A. K. (2018). Density Functional Electric Response Properties of Molecules in Cartesian Grid. Int. J. Quantum Chem. 118, e25708. doi:10.1002/qua.25708

CrossRef Full Text | Google Scholar

Ghosal, A., Mandal, T., and Roy, A. K. (2019). Efficient HF Exchange Evaluation through Fourier Convolution in Cartesian Grid for Orbital-dependent Density Functionals. J. Chem. Phys. 150, 064104. doi:10.1063/1.5082393

PubMed Abstract | CrossRef Full Text | Google Scholar

Ghosal, A., and Roy, A. K. (2016). in Specialist Periodical Reports: Chemical Modelling, Applications and Theory. Editors M. Springborg, and J.-O. Joswig (London: Royal Society of Chemistry), Vol. 13.

Google Scholar

Ghosal, A., and Roy, A. K. (2022b). A Real-Time TDDFT Scheme for Strong-Field Interaction in Cartesian Coordinate Grid. Chem. Phys. Lett. 796, 139562. doi:10.1016/j.cplett.2022.139562

CrossRef Full Text | Google Scholar

Ghosal, A., and Roy, A. K. (2022a). A Self-Consistent Systematic Optimisation of Range-Separated Hybrid Functionals from First Principles. Mol. Phys. 120, e1983056. doi:10.1080/00268976.2021.1983056

CrossRef Full Text | Google Scholar

Gill, P. M. W., and Adamson, R. D. (1996). A Family of Attenuated Coulomb Operators. Chem. Phys. Lett. 261, 105–110. doi:10.1016/0009-2614(96)00931-1

CrossRef Full Text | Google Scholar

Gill, P. M. W., Adamson, R. D., and Pople, J. A. (1996). Coulomb-attenuated Exchange Energy Density Functionals. Mol. Phys. 88, 1005–1009. doi:10.1080/00268979609484488

CrossRef Full Text | Google Scholar

Görling, A., and Levy, M. (1997). Hybrid Schemes Combining the Hartree-Fock Method and Density-Functional Theory: Underlying Formalism and Properties of Correlation Functionals. J. Chem. Phys. 106, 2675–2680. doi:10.1063/1.473369

CrossRef Full Text | Google Scholar

Grimme, S., and Parac, M. (2003). Substantial Errors from Time-dependent Density Functional Theory for the Calculation of Excited States of Large π Systems. ChemPhysChem 4, 292–295. doi:10.1002/cphc.200390047

PubMed Abstract | CrossRef Full Text | Google Scholar

Guido, C. A., Brémond, E., Adamo, C., and Cortona, P. (2013). Communication: One Third: A New Recipe for the PBE0 Paradigm. J. Chem. Phys. 138, 021104. doi:10.1063/1.4775591

PubMed Abstract | CrossRef Full Text | Google Scholar

Hait, D., Zhu, T., McMahon, D. P., and Van Voorhis, T. (2016). Prediction of Excited-State Energies and Singlet-Triplet Gaps of Charge-Transfer States Using a Restricted Open-Shell Kohn-Sham Approach. J. Chem. Theory Comput. 12, 3353–3359. doi:10.1021/acs.jctc.6b00426

PubMed Abstract | CrossRef Full Text | Google Scholar

Handy, N. C., and Cohen, A. J. (2001). Left-right Correlation Energy. Mol. Phys. 99, 403–412. doi:10.1080/00268970010018431

CrossRef Full Text | Google Scholar

Harris, J., and Jones, R. O. (1974). The Surface Energy of a Bounded Electron Gas. J. Phys. F. Met. Phys. 4, 1170–1186. doi:10.1088/0305-4608/4/8/013

CrossRef Full Text | Google Scholar

Head-Gordon, M., and Pople, J. A. (1988). A Method for Two‐electron Gaussian Integral and Integral Derivative Evaluation Using Recurrence Relations. J. Chem. Phys. 89, 5777.

Google Scholar

Helgaker, T., Coriani, S., Jørgensen, P., Kristensen, K., Olsen, J., and Ruud, K. (2012). Recent Advances in Wave Function-Based Methods of Molecular-Property Calculations. Chem. Rev. 112, 543–631. doi:10.1021/cr2002239

PubMed Abstract | CrossRef Full Text | Google Scholar

Henderson, T. M., Janesko, B. G., and Scuseria, G. E. (2008). Generalized Gradient Approximation Model Exchange Holes for Range-Separated Hybrids. J. Chem. Phys. 128, 194105. doi:10.1063/1.2921797

PubMed Abstract | CrossRef Full Text | Google Scholar

Hine, N. D. M., Dziedzic, J., Haynes, P. D., and Skylaris, C.-K. (2011). Electrostatic Interactions in Finite Systems Treated with Periodic Boundary Conditions: Application to Linear-Scaling Density Functional Theory. J. Chem. Phys. 135, 204103. doi:10.1063/1.3662863

PubMed Abstract | CrossRef Full Text | Google Scholar

Hoffman, E. O. (2007). Progress in Quantum Chemistry Research. New York: Nova Science Publishers.

Google Scholar

Hohenberg, P., and Kohn, W. (1964). Inhomogeneous Electron Gas. Phys. Rev. 136, B864–B871. doi:10.1103/physrev.136.b864

CrossRef Full Text | Google Scholar

Hohm, U. (2013). Experimental Static Dipole-Dipole Polarizabilities of Molecules. J. Mol. Struct. 1054-1055, 282–292. doi:10.1016/j.molstruc.2013.10.003

CrossRef Full Text | Google Scholar

Hou, X., Ke, C., Bruns, C. J., McGonigal, P. R., Pettman, R. B., and Stoddart, J. F. (2015). Tunable Solid-State Fluorescent Materials for Supramolecular Encryption. Nat. Commun. 6, 6884. doi:10.1038/ncomms7884

PubMed Abstract | CrossRef Full Text | Google Scholar

Iikura, H., Tsuneda, T., Yanai, T., and Hirao, K. (2001). A Long-Range Correction Scheme for Generalized-Gradient-Approximation Exchange Functionals. J. Chem. Phys. 115, 3540–3544. doi:10.1063/1.1383587

CrossRef Full Text | Google Scholar

Janak, J. F. (1978). Proof that∂E∂ni=εin Density-Functional Theory. Phys. Rev. B 18, 7165–7168. doi:10.1103/physrevb.18.7165

CrossRef Full Text | Google Scholar

Janesko, B. G. (2021). Replacing Hybrid Density Functional Theory: Motivation and Recent Advances. Chem. Soc. Rev. 50, 8470–8495. doi:10.1039/d0cs01074j

PubMed Abstract | CrossRef Full Text | Google Scholar

Jansik, B., Sałek, P., Jonsson, D., Vahtras, O., and Ågren, H. (2005). Cubic Response Functions in Time-dependent Density Functional Theory. J. Chem. Phys. 122, 054107. doi:10.1063/1.1811605

CrossRef Full Text | Google Scholar

Johnson, B. G., Gill, P. M. W., and Pople, J. A. (1993). The Performance of a Family of Density Functional Methods. J. Chem. Phys. 98, 5612–5626. doi:10.1063/1.464906

CrossRef Full Text | Google Scholar

Khan, M. U., Hussain, R., Yasir Mehboob, M., Khalid, M., Shafiq, Z., Aslam, M., et al. (2020). In Silico Modeling of New "Y-Series"-Based Near-Infrared Sensitive Non-fullerene Acceptors for Efficient Organic Solar Cells. ACS Omega 5, 24125–24137. doi:10.1021/acsomega.0c03796

PubMed Abstract | CrossRef Full Text | Google Scholar

Koch, W., and Holthausen, M. C. (2001). A Chemist’s Guide to Density Functional Theory. New York: John Wiley.

Google Scholar

Kohn, W. (1999). Nobel Lecture: Electronic Structure of Matter-Wave Functions and Density Functionals. Rev. Mod. Phys. 71, 1253–1266. doi:10.1103/revmodphys.71.1253

CrossRef Full Text | Google Scholar

Kohn, W., and Sham, L. J. (1965). Self-Consistent Equations Including Exchange and Correlation Effects. Phys. Rev. 140, A1133–A1138. doi:10.1103/physrev.140.a1133

CrossRef Full Text | Google Scholar

Kong, J., and Proynov, E. (2015). Density Functional Model for Nondynamic and Strong Correlation. J. Chem. Theory Comput. 12, 133–143. doi:10.1021/acs.jctc.5b00801

PubMed Abstract | CrossRef Full Text | Google Scholar

Kossmann, S., and Neese, F. (2009). Comparison of Two Efficient Approximate Hartee-Fock Approaches. Chem. Phys. Lett. 481, 240–243. doi:10.1016/j.cplett.2009.09.073

CrossRef Full Text | Google Scholar

Kowalczyk, T., Yost, S. R., and Voorhis, T. V. (2011). Assessment of the ΔSCF Density Functional Theory Approach for Electronic Excitations in Organic Dyes. J. Chem. Phys. 134, 054128. doi:10.1063/1.3530801

PubMed Abstract | CrossRef Full Text | Google Scholar

Kronik, L., and Kümmel, S. (2018). Dielectric Screening Meets Optimally Tuned Density Functionals. Adv. Mat. 30, 1706560. doi:10.1002/adma.201706560

PubMed Abstract | CrossRef Full Text | Google Scholar

Kronik, L., and Kümmel, S. (2020). Piecewise Linearity, Freedom from Self-Interaction, and a Coulomb Asymptotic Potential: Three Related yet Inequivalent Properties of the Exact Density Functional. Phys. Chem. Chem. Phys. 22, 16467–16481. doi:10.1039/d0cp02564j

PubMed Abstract | CrossRef Full Text | Google Scholar

Kümmel, S., and Kronik, L. (2006). Hyperpolarizabilities of Molecular Chains: A Real-Space Approach. Comput. Mater. Sci. 35, 321–326. doi:10.1016/j.commatsci.2004.09.057

CrossRef Full Text | Google Scholar

Kundu, K., Knight, S. F., Willett, N., Lee, S., Taylor, W. R., and Murthy, N. (2009). Hydrocyanines: A Class of Fluorescent Sensors that Can Image Reactive Oxygen Species in Cell Culture, Tissue, and In Vivo. Angew. Chem. Int. Ed. 48, 299–303. doi:10.1002/anie.200804851

CrossRef Full Text | Google Scholar

Kurtz, H. A., Stewart, J. J. P., and Dieter, K. M. (1990). Calculation of the Nonlinear Optical Properties of Molecules. J. Comput. Chem. 11, 82–87. doi:10.1002/jcc.540110110

CrossRef Full Text | Google Scholar

Labello, N. P., Ferreira, A. M., and Kurtz, H. A. (2005). An Augmented Effective Core Potential Basis Set for the Calculation of Molecular Polarizabilities. J. Comput. Chem. 26, 1464–1471. doi:10.1002/jcc.20282

PubMed Abstract | CrossRef Full Text | Google Scholar

Lange, A. W., Rohrdanz, M. A., and Herbert, J. M. (2008). Charge-Transfer Excited States in a π-Stacked Adenine Dimer, as Predicted Using Long-Range-Corrected Time-dependent Density Functional Theory. J. Phys. Chem. B 112, 6304–6308. doi:10.1021/jp802058k

PubMed Abstract | CrossRef Full Text | Google Scholar

Laqua, H., Kussmann, J., and Ochsenfeld, C. (2018). Efficient and Linear-Scaling Seminumerical Method for Local Hybrid Density Functionals. J. Chem. Theory Comput. 14, 3451–3458. doi:10.1021/acs.jctc.8b00062

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, C., Yang, W., and Parr, R. G. (1988). Development of the Colle-Salvetti Correlation-Energy Formula into a Functional of the Electron Density. Phys. Rev. B 37, 785–789. doi:10.1103/physrevb.37.785

PubMed Abstract | CrossRef Full Text | Google Scholar

Leininger, T., Stoll, H., Werner, H.-J., and Savin, A. (1997). Combining Long-Range Configuration Interaction with Short-Range Density Functionals. Chem. Phys. Lett. 275, 151–160. doi:10.1016/s0009-2614(97)00758-6

CrossRef Full Text | Google Scholar

Levy, M., Perdew, J. P., and Sahni, V. (1984). Exact Differential Equation for the Density and Ionization Energy of a Many-Particle System. Phys. Rev. A 30, 2745–2748. doi:10.1103/physreva.30.2745

CrossRef Full Text | Google Scholar

Li, C., Lu, J., and Yang, W. (2015). Gentlest Ascent Dynamics for Calculating First Excited State and Exploring Energy Landscape of Kohn-Sham Density Functionals. J. Chem. Phys. 143, 224110. doi:10.1063/1.4936411

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, L.-L., Li, K., Liu, Y.-H., Xu, H.-R., and Yu, X.-Q. (2016). Red Emission Fluorescent Probes for Visualization of Monoamine Oxidase in Living Cells. Sci. Rep. 6, 31217. doi:10.1038/srep31217

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, Y.-S., Tsai, C.-W., Li, G.-D., and Chai, J.-D. (2012). Long-range Corrected Hybrid Meta-Generalized-Gradient Approximations with Dispersion Corrections. J. Chem. Phys. 136, 154109. doi:10.1063/1.4704370

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, F., Furlani, T., and Kong, J. (2016). Optimal Path Search for Recurrence Relation in Cartesian Gaussian Integrals. J. Phys. Chem. A 120, 10264–10272. doi:10.1021/acs.jpca.6b10468

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, F., and Kong, J. (2017). Efficient Computation of Exchange Energy Density with Gaussian Basis Functions. J. Chem. Theory Comput. 13, 2571–2580. doi:10.1021/acs.jctc.7b00055

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, F., Proynov, E., Yu, J.-G., Furlani, T. R., and Kong, J. (2012). Comparison of the Performance of Exact-Exchange-Based Density Functional Methods. J. Chem. Phys. 137, 114104. doi:10.1063/1.4752396

PubMed Abstract | CrossRef Full Text | Google Scholar

Livshits, E., and Baer, R. (2007). A Well-Tempered Density Functional Theory of Electrons in Molecules. Phys. Chem. Chem. Phys. 9, 2932. doi:10.1039/b617919c

PubMed Abstract | CrossRef Full Text | Google Scholar

Mandal, T., Ghosal, A., and Roy, A. K. (2019). Static Polarizability and Hyperpolarizability in Atoms and Molecules through a Cartesian-Grid DFT. Theor. Chem. Acc. 138, 10. doi:10.1007/s00214-018-2397-7

CrossRef Full Text | Google Scholar

Mangiatordi, G. F., Brémond, E., and Adamo, C. (2012). DFT and Proton Transfer Reactions: A Benchmark Study on Structure and Kinetics. J. Chem. Theory Comput. 8, 3082–3088. doi:10.1021/ct300338y

PubMed Abstract | CrossRef Full Text | Google Scholar

Mardirossian, N., and Head-Gordon, M. (2017). Thirty Years of Density Functional Theory in Computational Chemistry: an Overview and Extensive Assessment of 200 Density Functionals. Mol. Phys. 115, 2315–2372. doi:10.1080/00268976.2017.1333644

CrossRef Full Text | Google Scholar

Maroulis, G. (1998). A Systematic Study of Basis Set, Electron Correlation, and Geometry Effects on the Electric Multipole Moments, Polarizability, and Hyperpolarizability of HCl. J. Chem. Phys. 108, 5432–5448. doi:10.1063/1.475932

CrossRef Full Text | Google Scholar

Maroulis, G., and Thakkar, A. J. (1988). Multipole Moments, Polarizabilities, and Hyperpolarizabilities for N2from Fourth‐order Many‐body Perturbation Theory Calculations. J. Chem. Phys. 88, 7623–7632. doi:10.1063/1.454327

CrossRef Full Text | Google Scholar

Martin, R. M. (2004). Electronic Structure: Basic Theory and Practical Methods. Cambridge, UK: Cambridge University Press.

Google Scholar

Martyna, G. J., and Tuckerman, M. E. (1999). A Reciprocal Space Based Method for Treating Long Range Interactions Inab Initioand Force-Field-Based Calculations in Clusters. J. Chem. Phys. 110, 2810–2821. doi:10.1063/1.477923

CrossRef Full Text | Google Scholar

Mclean, A. D., and Yoshimine, M. (1967). Theory of Molecular Polarizabilities. J. Chem. Phys. 47, 1927–1935. doi:10.1063/1.1712220

CrossRef Full Text | Google Scholar

Miller, T. M., and Bederson, B. (1997). Advances in Atomic and Molecular Physics, 13, 1. doi:10.1016/S0065-2199(08)60054-8

CrossRef Full Text | Google Scholar

Mináry, P., Tuckerman, M. E., Pihakari, K. A., and Martyna, G. J. (2002). A New Reciprocal Space Based Treatment of Long Range Interactions on Surfaces. J. Chem. Phys. 116, 5351–5362. doi:10.1063/1.1453397

CrossRef Full Text | Google Scholar

Mitschke, U., and Bäuerle, P. (2000). The Electroluminescence of Organic Materials. J. Mat. Chem. 10, 1471–1507. doi:10.1039/a908713c

CrossRef Full Text | Google Scholar

Mori-Sánchez, P., Cohen, A. J., and Yang, W. (2006). Self-interaction-free Exchange-Correlation Functional for Thermochemistry and Kinetics. J. Chem. Phys. 124, 091102. doi:10.1063/1.2179072

CrossRef Full Text | Google Scholar

Mulliken, R. S., and Person, W. B. (1962). Donor-Acceptor Complexes. Annu. Rev. Phys. Chem. 13, 107–126. doi:10.1146/annurev.pc.13.100162.000543

CrossRef Full Text | Google Scholar

Mulliken, R. S., and Person, W. B. (1969). Molecular Compounds and Their Spectra. XXI. Some General Considerations. J. Am. Chem. Soc. 91, 3409–3413. doi:10.1021/ja01041a001

CrossRef Full Text | Google Scholar

Murata, N. (1969). Control of Excitation Transfer in Photosynthesis I. Light-Induced Change of Chlorophyll a Fluoresence in Porphyridium Cruentum. Biochimica Biophysica Acta (BBA) - Bioenergetics 172, 242–251. doi:10.1016/0005-2728(69)90067-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Navizet, I., Liu, Y.-J., Ferré, N., Roca-Sanjuán, D., and Lindh, R. (2011). The Chemistry of Bioluminescence: An Analysis of Chemical Functionalities. ChemPhysChem 12, 3064–3076. doi:10.1002/cphc.201100504

PubMed Abstract | CrossRef Full Text | Google Scholar

Nelson, R. D., Lide, D. R., and Maryott, A. A. (1967). Selected Values of Electric Dipole Moments for Molecules in the Gas Phase. New York: National Standard Reference Data System.

Google Scholar

Obara, S., and Saika, A. (1986). Efficient Recursive Computation of Molecular Integrals over Cartesian Gaussian Functions. J. Chem. Phys. 84, 3963–3974. doi:10.1063/1.450106

CrossRef Full Text | Google Scholar

Obara, S., and Saika, A. (1988). General Recurrence Formulas for Molecular Integrals over Cartesian Gaussian Functions. J. Chem. Phys. 89, 1540–1559. doi:10.1063/1.455717

CrossRef Full Text | Google Scholar

Olney, T. N., Cann, N. M., Cooper, G., and Brion, C. E. (1997). Absolute Scale Determination for Photoabsorption Spectra and the Calculation of Molecular Properties Using Dipole Sum-Rules. Chem. Phys. 223, 59–98. doi:10.1016/s0301-0104(97)00145-6

CrossRef Full Text | Google Scholar

Orr, B. J., and Ward, J. F. (1971). Perturbation Theory of the Non-linear Optical Polarization of an Isolated System. Mol. Phys. 20, 513–526. doi:10.1080/00268977100100481

CrossRef Full Text | Google Scholar

Palczewski, K. (2012). Chemistry and Biology of Vision. J. Biol. Chem. 287, 1612–1619. doi:10.1074/jbc.r111.301150

PubMed Abstract | CrossRef Full Text | Google Scholar

Parac, M., and Grimme, S. (2003). A TDDFT Study of the Lowest Excitation Energies of Polycyclic Aromatic Hydrocarbons. Chem. Phys. 292, 11–21. doi:10.1016/s0301-0104(03)00250-7

CrossRef Full Text | Google Scholar

Parr, R. G., and Yang, W. (1989). Density Functional Theory of Atoms and Molecules. New York: Oxford University Press.

Google Scholar

Pascola, D., Costa, M. F., and Santos, H. F. D. (2014). International Journal of Quantum Chemistry, 114, 796. doi:10.1002/qua.24678

CrossRef Full Text | Google Scholar

Patel, A. H. G., Mohammed, A. A. K., Limacher, P. A., and Ayers, P. W. (2017). Finite Field Method for Nonlinear Optical Property Prediction Using Rational Function Approximants. J. Phys. Chem. A 121, 5313–5323. doi:10.1021/acs.jpca.7b04049

PubMed Abstract | CrossRef Full Text | Google Scholar

Perdew, J. P., Burke, K., and Ernzerhof, M. (1996a). Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 77, 3865–3868. doi:10.1103/physrevlett.77.3865

PubMed Abstract | CrossRef Full Text | Google Scholar

Perdew, J. P., Ernzerhof, M., and Burke, K. (1996b). Rationale for Mixing Exact Exchange with Density Functional Approximations. J. Chem. Phys. 105, 9982–9985. doi:10.1063/1.472933

CrossRef Full Text | Google Scholar

Perdew, J. P., Parr, R. G., Levy, M., and Balduz, J. L. (1982). Density-Functional Theory for Fractional Particle Number: Derivative Discontinuities of the Energy. Phys. Rev. Lett. 49, 1691–1694. doi:10.1103/physrevlett.49.1691

CrossRef Full Text | Google Scholar

Perdew, J. P., and Schmidt, K. (2000). Jacob's Ladder of Density Functional Approximations for the Exchange-Correlation Energy. AIP Conf. Proc. 577, 1. doi:10.1063/1.1390175

CrossRef Full Text | Google Scholar

Perdew, J. P., and Wang, Y. (1992). Accurate and Simple Analytic Representation of the Electron-Gas Correlation Energy. Phys. Rev. B 45, 13244–13249. doi:10.1103/physrevb.45.13244

PubMed Abstract | CrossRef Full Text | Google Scholar

Perdew, J. P., and Zunger, A. (1981). Self-interaction Correction to Density-Functional Approximations for Many-Electron Systems. Phys. Rev. B 23, 5048–5079. doi:10.1103/physrevb.23.5048

CrossRef Full Text | Google Scholar

Peverati, R., and Truhlar, D. G. (2012a). Exchange-Correlation Functional with Good Accuracy for Both Structural and Energetic Properties while Depending Only on the Density and its Gradient. J. Chem. Theory Comput. 8, 2310–2319. doi:10.1021/ct3002656

PubMed Abstract | CrossRef Full Text | Google Scholar

Peverati, R., and Truhlar, D. G. (2014). Quest for a Universal Density Functional: the Accuracy of Density Functionals across a Broad Spectrum of Databases in Chemistry and Physics. Phil. Trans. R. Soc. A 372, 20120476. doi:10.1098/rsta.2012.0476

PubMed Abstract | CrossRef Full Text | Google Scholar

Peverati, R., and Truhlar, D. G. (2012b). Screened-exchange Density Functionals with Broad Accuracy for Chemistry and Solid-State Physics. Phys. Chem. Chem. Phys. 14, 16187. doi:10.1039/c2cp42576a

PubMed Abstract | CrossRef Full Text | Google Scholar

Pople, J. A., Gill, P. M. W., and Johnson, B. G. (1992). Kohn-Sham Density-Functional Theory within a Finite Basis Set. Chem. Phys. Lett. 199, 557–560. doi:10.1016/0009-2614(92)85009-y

CrossRef Full Text | Google Scholar

Proynov, E., Liu, F., Shao, Y., and Kong, J. (2012). Improved Self-Consistent and Resolution-Of-Identity Approximated Becke'05 Density Functional Model of Nondynamic Electron Correlation. J. Chem. Phys. 136, 034102. doi:10.1063/1.3676726

PubMed Abstract | CrossRef Full Text | Google Scholar

Proynov, E., Shao, Y., and Kong, J. (2010). Efficient Self-Consistent DFT Calculation of Nondynamic Correlation Based on the B05 Method. Chem. Phys. Lett. 493, 381–385. doi:10.1016/j.cplett.2010.05.029

PubMed Abstract | CrossRef Full Text | Google Scholar

Pugh, D. (2006). in Specialist Periodical Reports: Chemical Modelling, Applications and Theory. Editor A. Hinchliffe (London: Royal Society of Chemistry), Vol. 4.

Google Scholar

Ramos, P., and Pavanello, M. (2016). Constrained Subsystem Density Functional Theory. Phys. Chem. Chem. Phys. 18, 21172–21178. doi:10.1039/c6cp00528d

PubMed Abstract | CrossRef Full Text | Google Scholar

R. D. Johnson (Editor) (2016). NIST Computational Chemistry Comparisons and Benchmark Database, NIST Standard Reference Database, Number, Release 18 (Gaithersburg, MD: NIST).

Google Scholar

Risthaus, T., Hansen, A., and Grimme, S. (2014). Excited States Using the Simplified Tamm-Dancoff-Approach for Range-Separated Hybrid Density Functionals: Development and Application. Phys. Chem. Chem. Phys. 16, 14408–14419. doi:10.1039/c3cp54517b

PubMed Abstract | CrossRef Full Text | Google Scholar

Rohrdanz, M. A., Martins, K. M., and Herbert, J. M. (2009). A Long-Range-Corrected Density Functional that Performs Well for Both Ground-State Properties and Time-dependent Density Functional Theory Excitation Energies, Including Charge-Transfer Excited States. J. Chem. Phys. 130, 054112. doi:10.1063/1.3073302

PubMed Abstract | CrossRef Full Text | Google Scholar

Roy, A. K. (2008a). Grid-based Density Functional Calculations of Many-Electron Systems. Int. J. Quantum Chem. 108, 837–847. doi:10.1002/qua.21570

CrossRef Full Text | Google Scholar

Roy, A. K. (2008b). Pseudopotential Density Functional Treatment of Atoms and Molecules in Cartesian Coordinate Grid. Chem. Phys. Lett. 461, 142–149. doi:10.1016/j.cplett.2008.06.076

CrossRef Full Text | Google Scholar

Roy, A. K. (2010a). “A New Density Functional Method for Electronic Structure Calculation of Atoms and Molecules,” in Handbook of Computational Chemistry Research. Editors C. T. Collett, and C. D. Robson (New York: Nova Science Publishers).

Google Scholar

Roy, A. K. (2010b). Trends in Physical Chemistry, 14, 27.

Google Scholar

Roy, A. K. (2011). Density Functional Calculation of Many-Electron Systems in Cartesian Coordinate Grid. J. Math. Chem. 49, 1687–1699. doi:10.1007/s10910-011-9851-2

CrossRef Full Text | Google Scholar

Roy, A. K., Ghosal, A., and Mandal, T. (2019). InDFT: A DFT Program for Atoms and Molecules in CCG. IISER Kolkata, India: Theoretical Chemistry Laboratory.

Google Scholar

Roy, R., Ghosal, A., and Roy, A. K. (2021). A Simple Effective Δ SCF Method for Computing Optical Gaps in Organic Chromophores. Chem. Asian J. 16, 2729–2739. doi:10.1002/asia.202100692

PubMed Abstract | CrossRef Full Text | Google Scholar

Roy, R., Ghosal, A., and Roy, A. K. (2022). Charge-Transfer Excitation within a Hybrid-(G)KS Framework through Cartesian Grid DFT. J. Phys. Chem. A 126, 1448–1457. doi:10.1021/acs.jpca.1c10593

PubMed Abstract | CrossRef Full Text | Google Scholar

Saad, Y. (2003). Iterative Methods for Sparse Linear Systems. Philadelphia: SIAM. doi:10.1137/1.9780898718003

CrossRef Full Text | Google Scholar

Sadlej, A. J. (1992). Medium-size Polarized Basis Sets for High-Level-Correlated Calculations of Molecular Electric Properties. Theor. Chim. Acta 81, 339–354. doi:10.1007/bf01118573

CrossRef Full Text | Google Scholar

Salzner, U., and Baer, R. (2009). Koopmans' Springs to Life. J. Chem. Phys. 131, 231101. doi:10.1063/1.3269030

PubMed Abstract | CrossRef Full Text | Google Scholar

Sambe, H., and Felton, R. H. (1975). A New Computational Approach to Slater's SCF-Xα Equation. J. Chem. Phys. 62, 1122–1126. doi:10.1063/1.430555

CrossRef Full Text | Google Scholar

Schmidt, M. W., Baldridge, K. K., Boatz, J. A., Elbert, S. T., Gordon, M. S., Jensen, J. H., et al. (1993). General Atomic and Molecular Electronic Structure System. J. Comput. Chem. 14, 1347–1363. doi:10.1002/jcc.540141112

CrossRef Full Text | Google Scholar

Seidl, A., Görling, A., Vogl, P., Majewski, J. A., and Levy, M. (1996). Generalized Kohn-Sham Schemes and the Band-Gap Problem. Phys. Rev. B 53, 3764–3774. doi:10.1103/physrevb.53.3764

CrossRef Full Text | Google Scholar

Sekino, H., and Bartlett, R. J. (1993). Molecular Hyperpolarizabilities. J. Chem. Phys. 98, 3022–3037. doi:10.1063/1.464129

CrossRef Full Text | Google Scholar

Shelton, D. P., and Rice, J. E. (1994). Measurements and Calculations of the Hyperpolarizabilities of Atoms and Small Molecules in the Gas Phase. Chem. Rev. 94, 3–29. doi:10.1021/cr00025a001

CrossRef Full Text | Google Scholar

Shophy, K. B., Shedge, S. V., and Pal, S. (2008). The Journal of Physical Chemistry A, 112, 11266. doi:10.1021/jp806204q

PubMed Abstract | CrossRef Full Text | Google Scholar

Silva-Junior, M. R., Schreiber, M., Sauer, S. P. A., and Thiel, W. (2008). Benchmarks for Electronically Excited States: Time-dependent Density Functional Theory and Density Functional Theory Based Multireference Configuration Interaction. J. Chem. Phys. 129, 104103. doi:10.1063/1.2973541

PubMed Abstract | CrossRef Full Text | Google Scholar

Silva-Junior, M. R., Schreiber, M., Sauer, S. P. A., and Thiel, W. (2010). Benchmarks of Electronically Excited States: Basis Set Effects on CASPT2 Results. J. Chem. Phys. 133, 174318. doi:10.1063/1.3499598

PubMed Abstract | CrossRef Full Text | Google Scholar

Skeel, R. D., Tezcan, I., and Hardy, D. J. (2002). Multiple Grid Methods for Classical Molecular Dynamics. J. Comput. Chem. 23, 673–684. doi:10.1002/jcc.10072

PubMed Abstract | CrossRef Full Text | Google Scholar

Skylaris, C.-K., Mostofi, A. A., Haynes, P. D., Diéguez, O., and Payne, M. C. (2002). Nonorthogonal Generalized Wannier Function Pseudopotential Plane-Wave Method. Phys. Rev. B 66, 035119. doi:10.1103/physrevb.66.035119

CrossRef Full Text | Google Scholar

Smith, G. D. (1978). Numerical Solutions of Partial Differential Equation: Finite Difference Methods. New York: Oxford University Press.

Google Scholar

Stein, T., Eisenberg, H., Kronik, L., and Baer, R. (2010). Fundamental Gaps in Finite Systems from Eigenvalues of a Generalized Kohn-Sham Method. Phys. Rev. Lett. 105, 266802. doi:10.1103/physrevlett.105.266802

PubMed Abstract | CrossRef Full Text | Google Scholar

Stephens, P. J., Devlin, F. J., Chabalowski, C. F., and Frisch, M. J. (1994). Ab Initio Calculation of Vibrational Absorption and Circular Dichroism Spectra Using Density Functional Force Fields. J. Phys. Chem. 98, 11623–11627. doi:10.1021/j100096a001

CrossRef Full Text | Google Scholar

Stevens, W. J., Basch, H., and Krauss, M. (1984). Compact Effective Potentials and Efficient Shared‐exponent Basis Sets for the First‐ and Second‐row Atoms. J. Chem. Phys. 81, 6026–6033. doi:10.1063/1.447604

CrossRef Full Text | Google Scholar

Talman, J. D. (2012). Basis Set Independent Calculation of Molecular Polarizabilities. Phys. Rev. A 86, 022519. doi:10.1103/physreva.86.022519

CrossRef Full Text | Google Scholar

Tamblyn, I., Refaely-Abramson, S., Neaton, J. B., and Kronik, L. (2014). Simultaneous Determination of Structures, Vibrations, and Frontier Orbital Energies from a Self-Consistent Range-Separated Hybrid Functional. J. Phys. Chem. Lett. 5, 2734–2741. doi:10.1021/jz5010939

PubMed Abstract | CrossRef Full Text | Google Scholar

Tao, J., Perdew, J. P., Staroverov, V. N., and Scuseria, G. E. (2003). Climbing the Density Functional Ladder: Nonempirical Meta–Generalized Gradient Approximation Designed for Molecules and Solids. Phys. Rev. Lett. 91, 1464. doi:10.1103/physrevlett.91.146401

CrossRef Full Text | Google Scholar

Taouali, W., Casida, M. E., Darghouth, A. A. M. H. M., and Alimi, K. (2018). Theoretical Design of New Small Molecules with a Low Band-Gap for Organic Solar Cell Applications: DFT and TD-DFT Study. Comput. Mater. Sci. 150, 54–61. doi:10.1016/j.commatsci.2018.03.038

CrossRef Full Text | Google Scholar

Tawada, Y., Tsuneda, T., Yanagisawa, S., Yanai, T., and Hirao, K. (2004). A Long-Range-Corrected Time-dependent Density Functional Theory. J. Chem. Phys. 120, 8425–8433. doi:10.1063/1.1688752

PubMed Abstract | CrossRef Full Text | Google Scholar

Thomas, L. H. (1927). The Calculation of Atomic Fields. Math. Proc. Camb. Phil. Soc. 23, 542–548. doi:10.1017/s0305004100011683

CrossRef Full Text | Google Scholar

Vasiliev, I., and Chelikowsky, J. R. (2010). Real-space Calculations of Atomic and Molecular Polarizabilities Using Asymptotically Correct Exchange-Correlation Potentials. Phys. Rev. A 82, 012502. doi:10.1103/physreva.82.012502

CrossRef Full Text | Google Scholar

Verma, P., Janesko, B. G., Wang, Y., He, X., Scalmani, G., Frisch, M. J., et al. (2019). M11plus: A Range-Separated Hybrid Meta Functional with Both Local and Rung-3.5 Correlation Terms and High Across-The-Board Accuracy for Chemical Applications. J. Chem. Theory Comput. 15, 4804–4815. doi:10.1021/acs.jctc.9b00411

PubMed Abstract | CrossRef Full Text | Google Scholar

Verma, P., and Truhlar, D. G. (2020). Status and Challenges of Density Functional Theory. Trends Chem. 2, 302–318. doi:10.1016/j.trechm.2020.02.005

CrossRef Full Text | Google Scholar

Vikramaditya, T., Chai, J.-D., and Lin, S.-T. (2018). Impact of Non-empirically Tuning the Range-Separation Parameter of Long-Range Corrected Hybrid Functionals on Ionization Potentials, Electron Affinities, and Fundamental Gaps. J. Comput. Chem. 39, 2378–2384. doi:10.1002/jcc.25575

PubMed Abstract | CrossRef Full Text | Google Scholar

Vosko, S. H., Wilk, L., and Nusair, M. (1980). Accurate Spin-dependent Electron Liquid Correlation Energies for Local Spin Density Calculations: a Critical Analysis. Can. J. Phys. 58, 1200–1211. doi:10.1139/p80-159

CrossRef Full Text | Google Scholar

Wadt, W. R., and Hay, P. J. (1985). Ab Initio effective Core Potentials for Molecular Calculations. Potentials for Main Group Elements Na to Bi. J. Chem. Phys. 82, 284–298. doi:10.1063/1.448800

CrossRef Full Text | Google Scholar

Wang, C., and Zhang, Q. (2018). Understanding Solid-State Solvation-Enhanced Thermally Activated Delayed Fluorescence Using a Descriptor-Tuned Screened Range-Separated Functional. J. Phys. Chem. C 123, 4407–4416. doi:10.1021/acs.jpcc.8b08228

CrossRef Full Text | Google Scholar

Wang, M., John, D., Yu, J., Proynov, E., Liu, F., Janesko, B. G., et al. (2019). Performance of New Density Functionals of Nondynamic Correlation on Chemical Properties. J. Chem. Phys. 150, 204101. doi:10.1063/1.5082745

PubMed Abstract | CrossRef Full Text | Google Scholar

Winget, P., and Brédas, J.-L. (2011). Ground-State Electronic Structure in Charge-Transfer Complexes Based on Carbazole and Diarylamine Donors. J. Phys. Chem. C 115, 10823–10835. doi:10.1021/jp200666p

CrossRef Full Text | Google Scholar

Wouters, S., Van Speybroeck, V., and Van Neck, D. (2016). DMRG-CASPT2 Study of the Longitudinal Static Second Hyperpolarizability of All-Trans Polyenes. J. Chem. Phys. 145, 054120. doi:10.1063/1.4959817

PubMed Abstract | CrossRef Full Text | Google Scholar

Yanai, T., Tew, D. P., and Handy, N. C. (2004). A New Hybrid Exchange-Correlation Functional Using the Coulomb-Attenuating Method (CAM-B3lyp). Chem. Phys. Lett. 393, 51–57. doi:10.1016/j.cplett.2004.06.011

CrossRef Full Text | Google Scholar

Yang, W., Zhang, Y., and Ayers, P. W. (2000). Degenerate Ground States and a Fractional Number of Electrons in Density and Reduced Density Matrix Functional Theory. Phys. Rev. Lett. 84, 5172–5175. doi:10.1103/physrevlett.84.5172

PubMed Abstract | CrossRef Full Text | Google Scholar

York, D., and Yang, W. (1994). The Fast Fourier Poisson Method for Calculating Ewald Sums. J. Chem. Phys. 101, 3298–3300. doi:10.1063/1.467576

CrossRef Full Text | Google Scholar

Yu, H. S., Li, S. L., and Truhlar, D. G. (2016). Perspective: Kohn-Sham Density Functional Theory Descending a Staircase. J. Chem. Phys. 145, 130901. doi:10.1063/1.4963168

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, H. S., Zhang, W., Verma, P., He, X., and Truhlar, D. G. (2015). Nonseparable Exchange-Correlation Functional for Molecules, Including Homogeneous Catalysis Involving Transition Metals. Phys. Chem. Chem. Phys. 17, 12146–12160. doi:10.1039/c5cp01425e

PubMed Abstract | CrossRef Full Text | Google Scholar

Zangwill, A. (2015). A Half Century of Density Functional Theory. Phys. Today 68, 34–39. doi:10.1063/pt.3.2846

CrossRef Full Text | Google Scholar

Zhao, Y., and Truhlar, D. G. (2008). Density Functionals with Broad Applicability in Chemistry. Acc. Chem. Res. 41, 157–167. doi:10.1021/ar700111a

PubMed Abstract | CrossRef Full Text | Google Scholar

Ziegler, T., Krykunov, M., and Cullen, J. (2012). The Implementation of a Self-Consistent Constricted Variational Density Functional Theory for the Description of Excited States. J. Chem. Phys. 136, 124107. doi:10.1063/1.3696967

PubMed Abstract | CrossRef Full Text | Google Scholar

Ziegler, T., Rauk, A., and Baerends, E. J. (1977). On the Calculation of Multiplet Energies by the Hartree-Fock-Slater Method. Theor. Chim. Acta 43, 261–271. doi:10.1007/bf00551551

CrossRef Full Text | Google Scholar