# Excitation energies through Becke’s exciton model within a Cartesian-grid KS
DFT

Abhisek Ghosal, Tarun Gupta, Kishalay Mahato  Amlan K. Roy  Email:
<EMAIL_ADDRESS><EMAIL_ADDRESS> Department of Chemical Sciences
Indian Institute of Science Education and Research (IISER) Kolkata
Nadia, Mohanpur-741246, WB, India

###### Abstract

Photon-induced electronic excitations are ubiquitously observed in organic
chromophore. In this context, we present a simple, alternative time-
independent DFT procedure, for computation of single-particle excitation
energies, in particular, the lower bound excited singlet states, which are of
primary interest in photochemistry. This takes inspiration from recently
developed Becke’s exciton model, where a key step constitutes the accurate
evaluation of correlated singlet-triplet splitting energy. It introduces a
non-empirical model, both from “adiabatic connection theorem” and “virial
theorem” to analyze the role of 2e- integral in such calculations. The latter
quantity is efficiently mapped onto a real grid and computed accurately using
a purely numerical strategy. Illustrative calculations are performed on 10
$\pi$-electron organic chromophores within a Cartesian-grid implementation of
pseudopotential Kohn-Sham (KS) DFT, developed in our laboratory, taking SBKJC-
type basis functions within B3LYP approximation. The triplet and singlet
excitation energies corresponding to first singly excited configuration, are
found to be in excellent agreement with TD-B3LYP calculations. Further, we
perform the same for a set of larger molecular systems using the
asymptotically corrected LC-BLYP, in addition to B3LYP. A systematic
comparison with theoretical best estimates demonstrates the viability and
suitability of current approach in determining optical gaps, combining
predictive accuracy with moderate computational cost.

Keywords: Density functional theory, adiabatic connection theorem, virial
theorem, Becke’s exciton model, range-separated hybrid functional, singlet-
triplet splitting.

##  I Introduction

Over the past decades, Kohn-Sham density functional theory (KS-DFT) (Kohn and
Sham, 1965) has emerged as a very powerful and successful tool for ground-
state calculations of many-electron systems, such as atoms, molecules and
solids Becke (2014). Moreover, its time-dependent (TD) variant, generally
entitled linear-response (LR)-TDDFT has been developed for computing
excitation energies (Runge and Gross, 1984; Casida, 1995), mostly for discrete
spectrum. In comparison to the conventional wave function-based formalisms,
such as configuration interaction (CI), multireference CI, complete active
space self-consistent field, equation-of-motion coupled cluster etc., LR-TDDFT
has gained popularity due to its reasonable trade-off between accuracy and
efficiency. This is evident from its applicability to medium and large systems
(Casida and Huix-Rotllant, 2012). However, despite its huge success, it has
well-known difficulties regarding double excitation, charge transfer and
Rydberg excitation. This is mainly due to the adverse effects of exchange-
correlation (XC) potential within adiabatic approximation. A detailed analysis
could be found in (Maitra, 2016).

Besides TDDFT, attempts were also made to obtain excited states within a time-
independent DFT rubric–several elegant approaches are available in the
literature Theophilou (1979); Gross et al. (1988); Singh and Deb (1996); Roy
et al. (1997); Ayers et al. (2015). In particular, the $\Delta$SCF method
(Ziegler et al., 1977; Kowalczyk et al., 2011)was developed within the
standard self-consistent field iteration, by adopting non-Aufbau occupations
at each iterations, and it targets the non-Aufbau solutions using the ground-
state functional. It has favorable scaling like ground-state DFT, and hence
smaller computational resource compared to LR-TDDFT methods. Although it
already provides a good estimation of excitation energies for molecular
systems (Seidu et al., 2014a), its tendency to variationally collapse to the
ground state is a serious concern. Such features during the $\Delta$SCF
procedure has been characterized in recent times (Seidu et al., 2014b); this
is particularly severe in systems with a dense energy spectrum near the Fermi
energy levels. Several sophisticated schemes were put forth to alleviate this
issue, which include a few constrained-DFT formalisms (Ziegler et al., 2012;
Barca et al., 2014; Ramos and Pavanello, 2016) as well as gentlest ascent
dynamics and meta dynamics related methods (Li et al., 2015), etc.  These are
often somehow involved, but are quite successful where TDDFT fails to perform
well. Furthermore, a lot of recent developments have been placed in the field
of unconstrained excited-state orbital relaxation that have the same
complexity as normal ground-state or TDDFT, while using a simple $\Delta$SCF-
style ansatz (Frank et al., 1998; Kowalczyk et al., 2013; Shea and Neuscamman,
2018; Hait and Head-Gordon, 2020; Carter-Fenk and Herbert, 2020).  On the
other hand, orbital energies themselves have physical meaning as excitation
energies (Van Meer et al., 2014), but are quite sensitive to the choice of
density functionals and its implementations. In a recent article (Haiduke and
Bartlett, 2018), KS-DFT has been used successfully to extract the photo-
electron and electronic excitation spectrum for molecules from orbital
energies through range-separated (RS) XC functionals. Moreover, efforts are
also known the development of multi-reference (MR)-DFT by means of
configurations mixing such as, multiconfiguration pair DFT (Li Manni et al.,
2014) and MR-DFT with generalized auxiliary systems (Chen et al., 2017).

In a series of articles, Becke (Becke, 2016, 2018a, 2018b; Becke et al.,
2018), introduced a simple novel model for estimation of lowest single-
particle excitation energies via correlated singlet-triplet splitting (STS)
energy. It offers an accurate, economical way to compute optical gaps in large
molecules. The main focus here lies in the lowest singlet excited state, which
are of fundamental interest in photochemistry. While triplet states can be
easily obtained, _in principle_ , standard DFT cannot be used in a singlet
excited state due to its multi-determinantal nature. The present scheme
advocates two separate single-determinant DFT calculations—one for the closed-
shell ground state and another for open-shell lowest triplet excited state,
followed by a simple two-electron integral (Coulomb self energy) evaluation
corresponding to the HOMO-LUMO transition. We note that there is no concern
about standard density functionals applied on a given triplet excited state as
it is represented by a single Slater determinant and is characterized by Fermi
hole. Thus the main ingredient is the correlated STS energy which can be
approached by means of (i) adiabatic connection theorem (Becke, 2016) and (ii)
virial theorem (Becke, 2018a). In a sense, this avoids the configuration
mixing and is altogether non-empirical.

In the present work, we have adopted the above approach for lowest single-
particle excitation energy corresponding to singlet excited state in
molecules. The required computations are performed following a pseudopotential
KS-DFT implemented in Cartesian coordinate grid (CCG), as developed in our
laboratory (Roy, 2008a, b, 2009, 2010, 2011; Ghosal and Roy, 2016; Ghosal et
al., 2018; Mandal et al., 2019; Ghosal et al., 2019). The pertinent two-
electron integral (from appropriate orbitals obtained after solving the lowest
triplet excited state) is carried out numerically using a recently designed
algorithm (Ghosal et al., 2019). This employs the Fourier convolution theorem
in conjunction with a RS Coulomb interaction kernel. The latter is efficiently
mapped onto the real grid through a simple grid optimization prescription,
giving rise to some constraints in RS parameter. The “adiabatic” model (Becke,
2016) corresponding to correlated STS energy is also provided accurately,
leading to an easy route to compute the singlet excitation energy. Following
the arguments of Becke (Becke, 2018a), an elegant “virial” theorem is then
engaged to analyze the role of two-electron integral in determining the
correlated STS energy. Both procedures are applied to calculate optical gaps
arising from both $\pi\to\pi^{\star}$ and $n\to\pi^{\star}$ transitions in a
decent number of molecules adopting B3LYP (Becke, 1993) functional. For a
proper comparison, we carry out parallel calculations corresponding to lowest
triplet excited states, from GAMESS quantum chemistry package (Schmidt et al.,
1993) using TD-B3LYP method, which reveals a better accuracy. Further, to
extend the scope and applicability of this approximation, we execute the same
using “virial” theorem for a set of large molecular systems, and employing an
additional XC functional from RS hybrid (Baer et al., 2010) family. The
effectiveness of this scheme is illustrated by the respective statistical
analysis. The article is organized as follows. The underlying theorem of
single-particle excitation energy along with a brief summary of our KS-DFT
framework in CCG, is provided in Sec. II. Section III offers the necessary
computational and technical details. Finally, the feasibility, performance and
accuracy of our results are critically assessed in Sec. IV. Some concluding
remarks as well as the future prospects are given in Sec. V.

##  II Methodology

Let us consider an excitation of a given system, corresponding to an
electronic configuration $\varphi_{i}\varphi_{f}$ from a closed-shell ground
state. Assuming completely filled closed-shell core, this is spanned by four
Slater determinants: $|\varphi_{i}^{\alpha}\varphi_{f}^{\alpha}\rangle$,
$|\varphi_{i}^{\alpha}\varphi_{f}^{\beta}\rangle$,
$|\varphi_{i}^{\beta}\varphi_{f}^{\alpha}\rangle$ and
$|\varphi_{i}^{\beta}\varphi_{f}^{\beta}\rangle$, where $\alpha$ and $\beta$
denote up and down spin. Therefore, the coupled excited states are found by
diagonalizing the Hamiltonian matrix in the space of above four determinants.
As such, the singlet state is given by
$|\psi_{\mathrm{S}}\rangle=\frac{1}{\sqrt{2}}\\{|\varphi_{i}^{\alpha}\varphi_{f}^{\beta}\rangle-|\varphi_{i}^{\beta}\varphi_{f}^{\alpha}\rangle\\}$,
whereas the three degenerate triplet states are written as follows:
$\psi_{\mathrm{T}}=|\varphi_{i}^{\alpha}\varphi_{f}^{\alpha}\rangle\quad\mathrm{or}\quad\frac{1}{\sqrt{2}}\\{|\varphi_{i}^{\alpha}\varphi_{f}^{\beta}\rangle+|\varphi_{i}^{\beta}\varphi_{f}^{\alpha}\rangle\\}\quad\mathrm{or}\quad|\varphi_{i}^{\beta}\varphi_{f}^{\beta}\rangle$.
The corresponding singlet and triplet energies (identified by “S” and “T”
subscripts) are,

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

and

$E_{\mathrm{T}}=E^{\alpha\alpha}\quad\mathrm{or}\quad
E^{\alpha\beta}-K_{if}\quad\mathrm{or}\quad E^{\beta\beta},$ (2)

where $E^{\sigma_{1}\sigma_{2}}$ is the energy of a given determinant of form
$|\varphi_{i}^{\sigma_{1}}\varphi_{f}^{\sigma_{2}}\rangle$,
$(\sigma_{1},\sigma_{2})\\!\in\\!\\{\alpha,\beta\\}$, and $K_{if}$ is the 2e-
integral (or Coulomb self-energy of product of transition orbitals) defined
as,

$K_{if}=\int\int\frac{\varphi_{i}(\mathrm{\mathbf{r}}_{1})\varphi_{f}(\mathrm{\mathbf{r}}_{1})\varphi_{i}(\mathrm{\mathbf{r}}_{2})\varphi_{f}(\mathrm{\mathbf{r}}_{2})}{|\mathrm{\mathbf{r}}_{1}-\mathrm{\mathbf{r}}_{2}|}\
d\mathrm{\mathbf{r}}_{1}d\mathrm{\mathbf{r}}_{2}.$ (3)

Now combining Eqs. (1), (2), one can connect singlet and triplet excitation
energy as follows,

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

where $E_{0\mathrm{S}}=E_{\mathrm{S}}-E_{0},\
E_{0\mathrm{T}}=E_{\mathrm{T}}-E_{0}$, while $E_{0}$ signifies the ground-
state energy of a certain closed-shell system. But Eq. (4) is highly
inaccurate for singlet excitation energy, $E_{0\mathrm{S}}$. The problem lies
in the determination of $2K_{if}$ term, also called the zeroth-order (or
uncorrelated) STS energy. In order to tackle this issue, recently some novel
proposals (semi-empirical (Becke, 2016) as well as non-empirical (Becke, 2016,
2018a)) have appeared in the literature to derive a simplified formula for
correlated STS energy.

###  II.1 The adiabatic connection theorem

One can make use of the well-known “adiabatic connection” theorem (Harris and
Jones, 1974) to approach the correlated STS energy in single-particle
excitations, which we are discussed here. For a specific excited
configuration, all the orbitals involved in Eqs. (1), (2) are same. In that
case, singlet and triplet states have same density and non-interacting kinetic
energy. Then their energy difference can be written as ($\Delta
E_{\mathrm{STS}}^{0}=2K_{if}$),

$\displaystyle\Delta E_{\mathrm{STS}}=\Delta E_{\mathrm{STS}}^{0}+\Delta
E_{\mathrm{ST}}^{\mathrm{corr}},$ $\displaystyle\Delta
E_{\mathrm{STS}}=2K_{if}+\Delta E_{\mathrm{ST}}^{\mathrm{corr}},$ (5)

where $\Delta E_{\mathrm{ST}}^{\mathrm{corr}}$ represents the singlet-triplet
correlation energy difference. Recently, a non-empirical formula Becke (2016)
has been proposed for it, based on the inter-electronic cusp condition, and
its effect on electron correlation. Accordingly, it can be expressed as,

$\Delta E_{\mathrm{ST}}^{\mathrm{corr}}=-0.4\ \int 4\
\varphi_{i}^{2}(\mathrm{\mathbf{r}}_{1})\varphi_{j}^{2}(\mathrm{\mathbf{r}}_{1})\
z_{\mathrm{C}}^{2}\bigg{[}1-\frac{\ln(1+z_{\mathrm{C}})}{z_{\mathrm{C}}}\bigg{]}d\mathrm{\mathbf{r}}_{1}.$
(6)

The only unknown quantity, the correlation length $z_{\mathrm{C}}$, measures
the spatial extent of electron correlation in configuration,
$\varphi_{i}\varphi_{f}$. Now if one allows a “strictly correlated electrons”
limit, $z_{\mathrm{C}}$ can be written in terms of two-electron integral,
$K_{if}$, as in the following,

$0.4z_{\mathrm{C}}^{2}\int\ 4\
\varphi_{i}^{2}(\mathrm{\mathbf{r}}_{1})\varphi_{j}^{2}(\mathrm{\mathbf{r}}_{1})\
d\mathrm{\mathbf{r}}_{1}=2K_{if}.$ (7)

A detailed derivation could be found in (Becke, 2016).

###  II.2 The virial exciton model

From the above discussion, it is clear that, one may add a correlation
correction term, $\Delta E_{\mathrm{ST}}^{\mathrm{corr}}$, to $2K_{if}$ to
recover the desired correlated STS energy. As suggested in Becke (2018a), one
may also further proceed to approximate this by noting that it comprises the
kinetic and potential energy contributions,

$\Delta E_{\mathrm{STS}}^{\mathrm{corr}}=\Delta
T_{\mathrm{STS}}^{\mathrm{corr}}+\Delta V_{\mathrm{STS}}^{\mathrm{corr}}.$ (8)

Now invoking the standard virial theorem, the above equation can be simplified
as,

$\Delta E_{\mathrm{STS}}^{\mathrm{corr}}=\frac{1}{2}\Delta
V_{\mathrm{STS}}^{\mathrm{corr}}.$ (9)

One may further presume that the non-local part of $\Delta
V_{\mathrm{STS}}^{\mathrm{corr}}$ is dominated by pair-density effects. In
that scenario, correlation would then lower the potential energy of singlet
state relative to triplet state, by $\Delta
V_{\mathrm{STS}}^{\mathrm{corr}}=-2K_{if}$. Hence, one can write,

$\Delta E_{\mathrm{STS}}^{\mathrm{corr}}=-K_{if}.$ (10)

which surprisingly leads to a very simple relation as below,

$\Delta E_{\mathrm{STS}}=2K_{if}-K_{if}=K_{if}.$ (11)

It may be noted that, this model (as well as that in II(A)) is also a purely
non-empirical one. But obviously, this expression for correlated STS energy is
much simpler involving only the two-electron integral.

###  II.3 The semi-empirical approach

As implied, both the routes in Eqs. (5), (11) are purely non-empirical,
suggesting $\Delta E_{\mathrm{STS}}$ to be less than $2K_{if}$. This prompts
one to formally define a molecule-independent re-scaling parameter $f$ such
that,

$\Delta E_{\mathrm{STS}}=2fK_{if},\quad 0<f<1.$ (12)

and if determined semi-empirically, it might offer good excitation energies
overall, as long as de-localization error (J. P. Perdew, R. G. Parr, M. Levy,
and J. L. Balduz Jr., 1982) is not a serious issue. This led to a semi-
empirical technique (Becke, 2016) to obtain an optimum value of $f$. This was
accomplished in Becke (2016) by fitting the results with best estimated
theoretical excitation energy data set of (Silva-Junior et al., 2010a),
resulting in a value for $f$ as 0.486. It is surprisingly close to 0.5, as
also appears from a consideration of II(B). The same author has also reported
that the obtained mean average error (using this $f$) remains similar to that
in non-empirical calculations, in II(A).

##  III Computational details

Our desired quantity, the singlet excitation energy $E_{0\mathrm{S}}$ can be
obtained from two single-determinant calculations (both are spin polarized),
namely, one for closed-shell ground state and another for fully relaxed
triplet state. These are computed using our in-house pseudopotential KS-DFT
program in CCG A. K. Roy, A. Ghosal and T. Mandal (2019). An initial version
was first developed by the corresponding author in 2008, on the basis of works
(Roy, 2008a, b, 2009, 2010, 2011), which has been later extended by his group
Ghosal and Roy (2016); Ghosal et al. (2018); Mandal et al. (2019); Ghosal et
al. (2019). Since the details have been published in the above references,
here we give only some essential aspects. Accordingly, the single-particle KS
equation in presence of pseudopotential is written as (atomic unit employed
unless stated otherwise),

$\bigg{[}-\frac{1}{2}\nabla^{2}+v^{\mathrm{p}}_{\mathrm{ion}}(\mathrm{\mathbf{r}})+v_{\mathrm{ext}}(\mathrm{\mathbf{r}})+v_{\mathrm{h}}[\rho(\mathrm{\mathbf{r}})]+v_{\mathrm{xc}}[\rho(\mathrm{\mathbf{r}})]\bigg{]}\varphi_{i}^{\sigma}(\mathrm{\mathbf{r}})=\epsilon_{i}\varphi_{i}^{\sigma}(\mathrm{\mathbf{r}}),$
(13)

where $v^{\mathrm{p}}_{\mathrm{ion}}$ denotes the ionic pseudopotential,
written as below,

$v^{\mathrm{p}}_{\mathrm{ion}}(\mathrm{\mathbf{r}})=\sum_{\mathrm{\mathbf{R}}_{a}}v^{\mathrm{p}}_{\mathrm{ion},\mathrm{a}}(\mathrm{\mathbf{r}}-\mathrm{\mathbf{R}}_{a}).$
(14)

In the above equation, $v^{\mathrm{p}}_{\mathrm{ion},\mathrm{a}}$ signifies
ion-core pseudopotential associated with atom A, situated at
$\mathrm{\mathbf{R}}_{a}$. The classical Coulomb (Hartree) term,
$v_{\mathrm{h}}[\rho(\mathrm{\mathbf{r}})]$ describes usual electrostatic
interaction amongst valence electrons, whereas
$v_{\mathrm{xc}}[\rho(\mathrm{\mathbf{r}})]$ represents the non-classical XC
part of latter, and
$\\{\varphi^{\sigma}_{i},\sigma=\alpha\quad\mathrm{or}\quad\beta\\}$
corresponds to a set of $N$ occupied orthonormal spin molecular orbitals
(spin-MOs).

Now various quantities like localized basis function, two-electron potential,
MOs and electron density are directly set up on a 3D cubic box,

$\displaystyle r_{i}=r_{0}+(i-1)h_{r},\quad i=1,2,3,....,N_{r}~{},\quad
r_{0}=-\frac{N_{r}h_{r}}{2},\quad r\in\\{x,y,z\\},$ (15)

where $h_{r},N_{r}$ denote grid spacing and total number of points along three
directions respectively. While one-electron contributions of KS-Fock matrix
are evaluated through well-established recursion relations, the two-electron
matrix elements are computed via direct numerical integration in the grid,

$F_{\mu\nu}^{\mathrm{hxc}}=\langle\chi_{\mu}(\mathrm{\mathbf{r}}_{g})|v_{\mathrm{hxc}}(\mathrm{\mathbf{r}}_{g})|\chi_{\nu}(\mathrm{\mathbf{r}}_{g})\rangle=h_{x}h_{y}h_{z}\sum_{g}\chi_{\mu}(\mathrm{\mathbf{r}}_{g})v_{\mathrm{hxc}}(\mathrm{\mathbf{r}}_{g})\chi_{\nu}(\mathrm{\mathbf{r}}_{g}).$
(16)

where $v_{hxc}$ refers to the classical (Hartree) and non-classical (XC)
potential combined. The construction of HF exchange in KS-Fock matrix has been
documented in our earlier work (Ghosal et al., 2019); hence not repeated here.

The actual implementation of the central quantity, $K_{if}$, is quite
different from the original prescription in (Becke, 2016). Here, we use the
Fourier convolution theorem (FCT) using a range-separated technique in Coulomb
interaction kernel, instead of the multi-center numerical integral procedure
put forth in (Becke and Dickson, 1988). According to (Ghosal et al., 2019),
this two-electron integral $K_{if}$ can be rewritten as,

$\displaystyle
K_{if}=\int\int\frac{\varphi_{i}(\mathrm{\mathbf{r}}_{1})\varphi_{f}(\mathrm{\mathbf{r}}_{1})\varphi_{i}(\mathrm{\mathbf{r}}_{2})\varphi_{f}(\mathrm{\mathbf{r}}_{2})}{|\mathrm{\mathbf{r}}_{1}-\mathrm{\mathbf{r}}_{2}|}\
d\mathrm{\mathbf{r}}_{1}d\mathrm{\mathbf{r}}_{2}$
$\displaystyle=\int\varphi_{i}(\mathrm{\mathbf{r}}_{1})\varphi_{f}(\mathrm{\mathbf{r}}_{1})v_{if}(\mathrm{\mathbf{r}}_{1})d\mathrm{\mathbf{r}}_{1}.$
(17)

Now, the main concern lies in the evaluation of $v_{if}$ integral which is
related to electro static potential integral evaluation in (Ghosal et al.,
2019). Again, this can be recast as,

$\displaystyle
v_{if}(\mathrm{\mathbf{r}}_{1})=\int\frac{\varphi_{i}(\mathrm{\mathbf{r}}_{2})\varphi_{f}(\mathrm{\mathbf{r}}_{2})}{|\mathrm{\mathbf{r}}_{1}-\mathrm{\mathbf{r}}_{2}|}d\mathrm{\mathbf{r}}_{2}$
$\displaystyle=\int\frac{\varphi_{if}(\mathrm{\mathbf{r}}_{2})}{|\mathrm{\mathbf{r}}_{1}-\mathrm{\mathbf{r}}_{2}|}d\mathrm{\mathbf{r}}_{2}=\varphi_{if}(\mathrm{\mathbf{r}}_{1})\star
v^{\mathrm{c}}(\mathrm{\mathbf{r}}_{1})$ (18)

The last expression is in terms of convolution integral, where $\varphi_{if}$
denotes simple multiplication of ith and fth MO from lowest triplet excited
state and $v^{\mathrm{c}}(\mathrm{\mathbf{r}})$ represents the Coulomb
interaction kernel. Now one can invoke the FCT to further simplify this
integral,

$v_{if}(\mathrm{\mathbf{r}})=\mathcal{F}^{-1}\\{v^{\mathrm{c}}(\mathrm{\mathbf{k}})\varphi_{if}(\mathrm{\mathbf{k}})\\}\quad\mathrm{where}\quad\varphi_{if}(\mathrm{\mathbf{k}})=\mathcal{F}\\{\varphi_{if}(\mathrm{\mathbf{r}})\\}$
(19)

Here $v^{c}(\mathrm{\mathbf{k}})$ and $\varphi_{if}(\mathrm{\mathbf{k}})$
stand for Fourier integrals of Coulomb kernel and MOs respectively. The main
concern lies in an accurate mapping of the former, which has singularity at
$\mathrm{\mathbf{r}}\to 0$. In order to alleviate this problem, we apply a
simple RS technique from (Ghosal et al., 2019), expanding the kernel into
long- and short-range components with a suitably chosen RS parameter
($\zeta$). This gives rise to the following expression,

$\displaystyle
v^{c}(\mathrm{\mathbf{r}})=\frac{\mathrm{erf}(\zeta\mathrm{\mathbf{r}})}{\mathrm{\mathbf{r}}}+\frac{\mathrm{erfc}(\zeta\mathrm{\mathbf{r}})}{\mathrm{\mathbf{r}}}$
$\displaystyle
v^{c}(\mathrm{\mathbf{r}}_{g})=v^{c}_{long}(\mathrm{\mathbf{r}}_{g})+v^{c}_{short}(\mathrm{\mathbf{r}}_{g}).$
(20)

In the above equation, $\mathrm{erf}(x)$ and $\mathrm{erfc}(x)$ denote error
function and its complement respectively, while the second expression is
written in real-space CCG. The Fourier integral of Coulomb kernel can be
separated out as follows: (i) FT of short-range part is treated analytically
and (ii) long-range portion is computed directly from FFT of corresponding
real-space CCG values. Then the remaining problem lies in finding an optimum
value of parameter $\zeta_{\mathrm{opt}}$ for successful mapping of Coulomb
kernel in CCG from _first principles_. In this regard, we proposed a simple
procedure which is found to be sufficiently accurate for lowest triplet
excited states (as exemplified in results that follows) as well. This prompts
us to write,

$\zeta_{\mathrm{opt}}\equiv\underset{\zeta}{\mathrm{opt}}\
E_{\mathrm{tot}}=\underset{N_{x},N_{y},N_{z}}{\mathrm{opt}}E_{\emph{tot}},\quad\mathrm{at\
fixed}\ h_{\mathrm{\mathbf{r}}},$ (21)

using a suitably defined constraint Ghosal et al. (2019) $(\zeta\times L=7)$,
where $L\ (=N_{\mathrm{\mathbf{r}}}h_{\mathrm{\mathbf{r}}};\
\mathrm{\mathbf{r}}\in\\{x,y,z\\})$ is the smallest side of simulating box. In
the same spirit, other necessary quantities such as correlation length,
$z_{\mathrm{C}}$ and singlet-triplet correlation energy difference, $\Delta
E_{\mathrm{C}}$ are directly computed in real-space grid using pseudo KS
orbitals $\varphi_{i}$ and $\varphi_{f}$.

The strategy as outlined for HF exchange component, is implemented in case of
B3LYP, containing a fixed amount of former with conventional DFT XC
functional. As such, the XC energy corresponding to this functional is defined
as follows,

$E_{\mathrm{xc}}^{\mathrm{B3LYP}}=(1-a_{0})E^{\mathrm{x}}_{\mathrm{LSDA}}+a_{0}E^{\mathrm{x}}_{\mathrm{HF}}+a_{\mathrm{x}}E^{\mathrm{x}}_{\mathrm{B88}}+a_{\mathrm{c}}E^{\mathrm{c}}_{\mathrm{LYP}}+(1-a_{\mathrm{c}})E^{\mathrm{c}}_{\mathrm{VWN}}.$
(22)

Throughout our current presentation, we use the recommended values of
$a_{0},a_{\mathrm{x}},a_{\mathrm{c}}$, as advocated in Stephens et al. (1994),
i.e., $0.2$, $0.72$, $0.81$ for B3LYP. This is computed using KS orbitals
obtained from the solution of Eq. (13) in real-space CCG. The relevant LDA-
and GGA-type functionals in connection with B3LYP are: (i) Vosko-Wilk-Nusair
(VWN)–with the homogeneous electron gas correlation proposed in
parametrization formula V (Vosko et al., 1980) (ii) B88–incorporating Becke
(Becke, 1988) semi-local exchange (iii) Lee-Yang-Parr (LYP) (Lee et al., 1988)
semi-local correlation. All these are adopted from density functional
repository program (rep, 2001) except LDA.

The present work employs SBKJC (Stevens et al., 1984) effective core potential
basis sets for species containing Group II elements. These are imported from
EMSL Basis Set Library (Feller, 1996). The norm-conserving pseudopotential
matrix elements in contracted basis are collected from GAMESS (Schmidt et al.,
1993). The triplet calculations refer to _restricted_ open-shell, so that
orbitals $\varphi_{i}$ and $\varphi_{f}$ are well defined. The convergence
criteria imposed in this communication are slightly tighter than our earlier
works (Roy, 2008a, b, 2011); this is to generate accurate orbital energies,
especially for triplet excited states. Changes in following quantities were
followed during self-consistent field process, _viz.,_ (i) orbital energy
difference between two successive iterations and (ii) absolute deviation in
density matrix elements. They both were required to remain below a prescribed
threshold set to $10^{-8}$ a.u.; this ensured that total energy maintained a
convergence of at least this much, in general. In order to perform discrete
Fourier transform, standard FFTW3 package Frigo and Johnson (2005) was
invoked. The resulting generalized matrix-eigenvalue problem is solved through
standard LAPACK routine Anderson et al. (1999) easily. All molecular
calculations are performed with calculated geometry (B3LYP XC functional and
cc-PVTZ basis set), taken from NIST database (Johnson, 2016) (otherwise stated
below). Other details including scaling properties may be found in references
(Roy, 2008a, b, 2009, 2010, 2011; Ghosal and Roy, 2016; Ghosal et al., 2018;
Mandal et al., 2019; Ghosal et al., 2019).

##  IV Result and Discussion

A set of 10 selective organic chromophores ($\pi$ containing molecules) from
our previous work (Ghosal et al., 2019) has been chosen, excluding those which
have degenerate frontier orbitals, as the method in its present form is not
applicable for these. These results pertain to the pseudopotential
approximation neglecting the effects of core electrons. This may, by and
large, be justifiable on the ground that only valence excitations are taken
into account. By doing so, one can make a balance between accuracy and cost.
To establish the effectiveness of this approach, we have chosen “SBKJC” type
of pseudopotential basis set which does not contain any diffuse function, and
this suffices in our present scenario. It may be emphasized that, in this
exploratory study, only lowest singlet and triplet excited state corresponding
to first single excitation, is pursued for each molecule.

Each calculation begins by optimizing the non-uniform grid for a given
molecule, both in its closed-shell ground and lowest excited triplet state
separately, employing a simple procedure as reported earlier (Ghosal and Roy,
2016; Ghosal et al., 2018). Individual singlet excitation energies (in eV)
through B3LYP XC functional, are displayed in Table I, along with mean
absolute error (MAE) and mean error (ME) statistics. It contains both non-
empirical and empirical excitation values as following: (i) PR1 represents Eq.
(12) using semi-empirical re-scaling parameter $f=0.486$ (ii) PR2 presents Eq.
(5) using non-empirical model from “adiabatic connection” (iii) PR3 refers to
results from Eq. (11) using non-empirical model from virial theorem. To put
things in perspective, column 6 reports the respective TD-B3LYP energies
employing same functional and basis set, computed from GAMESSSchmidt et al.
(1993). The correlation lengths $z_{\mathrm{C}}$ manifestly fall into two
distinct groups with $\pi\to\pi^{\star}$ and $n\to\pi^{\star}$ transition–a
feature that has been observed in (Becke, 2016). Furthermore, in consonance
with their finding, we also witness a similar kind of pattern in the
calculated $z_{\mathrm{C}}$ values, having smaller $\varphi_{i}\varphi_{f}$
overlap for the $n\to\pi^{\star}$ transition for current set of molecules.
Overall, the obtained singlet excitation energies are within fraction of eV
compared to the TD-B3LYP results. As expected from (Becke, 2018a), the
correlated virial theorem (PR3) confers remarkable improvement over the
adiabatic connection formula (PR2). Here, the MAE value is about two times
smaller compared to the latter one. On the other hand, the semi-empirical
approach (PR2) performs very close to PR3 results, as evidenced from the
proximity in their corresponding MAE values. The worst case candidate is
propene, with an excitation energy $0.55$ eV too low, in comparison to
TD-B3LYP result. We note that the MAE and ME remain same in magnitude and
sign, indicating a systematic error in calculated values. Out of three
methods, the virial PR3 route provides the best estimate in exciting energies;
quite competitive with that of TD-B3LYP, with an MAE of 0.34 eV.

Table 1:  Singlet excitation energies (in eV) using B3LYP functional. See text for details. Molecule | State |  PR1 |  PR2 |  PR3 |  TD-B3LYP (Schmidt et al., 1993) | $z_{\mathrm{C}}$
---|---|---|---|---|---|---
Ethylene |  B1u($\pi\rightarrow\pi^{\star}$) | 7.78 | 7.63 | 7.87 | 8.09 | 2.97
Propene |  A′($\pi\rightarrow\pi^{\star}$) | 7.18 | 7.05 | 7.26 | 7.81 | 2.99
1,3-Butadiene (E) |  B($\pi\rightarrow\pi^{\star}$) | 5.63 | 5.42 | 5.70 | 6.02 | 3.29
1,3-Pentadiene (E) |  A′($\pi\rightarrow\pi^{\star}$) | 5.47 | 5.28 | 5.54 | 5.88 | 3.29
1,3,5-Hexatriene (E) |  Bu($\pi\rightarrow\pi^{\star}$) | 4.38 | 4.14 | 4.44 | 4.79 | 3.53
2,4-Hexadiene (E,E) |  Bu($\pi\rightarrow\pi^{\star}$) | 5.39 | 5.20 | 5.45 | 5.79 | 3.27
1,3-Cyclo-pentadiene |  A′($\pi\rightarrow\pi^{\star}$) | 5.12 | 5.03 | 5.17 | 5.28 | 2.98
Thiophene |  B2($\pi\rightarrow\pi^{\star}$) | 5.61 | 5.31 | 5.66 | 6.02 | 3.99
Formaldehyde |  A2($n\rightarrow\pi^{\star}$ | 3.52 | 3.63 | 3.53 | 3.98 | 1.30
Acetaldehyde |  A′′($n\rightarrow\pi^{\star}$) | 4.67 | 4.75 | 4.68 | 5.07 | 1.44
MAE | – | 0.40 | 0.53 | 0.34 | – | –
ME | – | 0.40 | 0.53 | 0.34 | – | –

Next we analyze the effects of triplet excitations and correlated STS term
separately for B3LYP functional in Table II, giving a side-by-side comparison
with TD-B3LYP in both cases. As singlet excitation energies are obtained from
these two above quantities, a detailed assessment of them is worthwhile
considering. Due to the superior performance, henceforth in future
proceedings, we engage only the PR3. The $E_{0\mathrm{T}}$ error metrics are:
MAE = 0.40 eV and ME = $-$0.39 eV, having very similar magnitude, but
different sign. The triplet MAE remains only marginally larger than singlet
MAE values of PR3 results of previous table. On the other hand, STS measures
read: MAE = 0.70 eV and ME = 0.69 eV, showing very close agreement in
magnitude and also importantly without changing sign. However magnitude of
these errors (MAE and ME) in STS are relatively larger than the $E_{0T}$
errors. A careful look on respective ME in these two quantities reveals that,
ME in STS is about 75% larger compared to that in E0T, but with an opposite
sign. And that makes the singlet excitation energies more close to TD-B3LYP
results, occurring through a systematic error cancellation. This leads us to
conclude that, _STS term and not_ $E_{0\mathrm{T}}$, is the major source of
error. It is worthwhile mentioning that, the success of this approach relies
on an accurate estimation of the two-electron integral, which in turn depends
on the accuracy of triplet excited state. Thus it is evident that our current
prescription can be reliably used for determining optical gaps using the non-
empirical time-independent methods of Sec. II, at least for small molecular
systems.

Table 2:  Triplet excitation energies and correlated STS energies (in eV) using B3LYP functional. Molecule | State | $E_{0\mathrm{T}}$ | $\Delta E_{\mathrm{STS}}$
---|---|---|---
|  |  Ref. (A. K. Roy, A. Ghosal and T. Mandal, 2019) | TD-B3LYP |  PR3 | TD-B3LYP
Ethylene |  B1u($\pi\rightarrow\pi^{\star}$) | 4.47 | 4.03 | 3.40 | 4.06
Propene |  A′($\pi\rightarrow\pi^{\star}$) | 4.44 | 4.03 | 2.82 | 3.78
1,3-Butadiene (E) |  B($\pi\rightarrow\pi^{\star}$) | 3.26 | 2.71 | 2.44 | 3.31
1,3-Pentadiene (E) |  A′($\pi\rightarrow\pi^{\star}$) | 3.24 | 2.71 | 2.30 | 3.17
1,3,5-Hexatriene (E) |  Bu($\pi\rightarrow\pi^{\star}$) | 2.42 | 1.85 | 2.02 | 2.94
2,4-Hexadiene (E,E) |  Bu($\pi\rightarrow\pi^{\star}$) | 3.22 | 2.70 | 2.23 | 3.09
1,3-Cyclo-pentadiene |  A′($\pi\rightarrow\pi^{\star}$) | 3.21 | 2.70 | 1.96 | 2.58
Thiophene |  B2($\pi\rightarrow\pi^{\star}$) | 3.88 | 3.47 | 1.78 | 2.55
Formaldehyde |  A2($n\rightarrow\pi^{\star}$) | 3.19 | 3.20 | 0.34 | 0.78
Acetaldehyde |  A′′($n\rightarrow\pi^{\star}$) | 4.39 | 4.44 | 0.29 | 0.24
MAE | – | 0.40 | – | 0.70 | –
ME | – |  $-$0.39 | – | 0.69 | –
Table 3:  Excitation energies (in eV) of organic chromophores from “virial
theorem”.

Molecule | State | $E_{0\mathrm{T}}$ |  $E_{0\mathrm{S}}$ (PR3) |  Lit.‡
---|---|---|---|---
|  | B3LYP | LC-BLYP | B3LYP | LC-BLYP |
Cyclopropene |  B2($\pi\rightarrow\pi^{\star}$) | 4.03 | 4.05 | 7.04 | 7.07 | 7.01
Norbornadiene |  A2($\pi\rightarrow\pi^{\star}$) | 4.62 | 4.23 | 5.77 | 5.45 | 4.91
Naphthalene |  B2u($\pi\rightarrow\pi^{\star}$) | 3.22 | 3.53 | 4.74 | 5.20 | 4.64
Furan |  B2($\pi\rightarrow\pi^{\star}$) | 4.15 | 4.25 | 6.65 | 6.85 | 6.57
Pyrrole |  B2($\pi\rightarrow\pi^{\star}$) | 4.45 | 4.52 | 6.88 | 7.07 | 6.85
Pyridine |  B1($n\rightarrow\pi^{\star}$) | NC | 4.22 | NC | 4.65 | 4.74
Pyridazine |  B3u($n\rightarrow\pi^{\star}$) | 2.76 | 2.78 | 3.35 | 3.34 | 3.57
s-tetrazine |  B2($\pi\rightarrow\pi^{\star}$) | 1.39 | 1.63 | 1.83 | 2.07 | 2.15
Acetone |  A2($n\rightarrow\pi^{\star}$) | 3.71 | 3.80 | 3.98 | 4.08 | 4.11
Formamide |  A′′($n\rightarrow\pi^{\star}$) | NC | 5.03 | NC | 5.25 | 5.42
Acetamide |  A′′($n\rightarrow\pi^{\star}$) | 5.20 | 5.15 | 5.45 | 5.37 | 5.46
Propanamide |  A′′($n\rightarrow\pi^{\star}$) | 5.64 | 5.19 | 6.07 | 5.40 | 5.49
Thymine |  A′($\pi\rightarrow\pi^{\star}$) | 3.51 | 3.59 | 5.50 | 5.73 | 5.48
Uracil |  A′′($n\rightarrow\pi^{\star}$) | 4.28 | 3.76 | 5.66 | 5.86 | 4.66
Imidazole |  A′($\Pi\rightarrow\pi^{\star}$) | NC | 4.75 | NC | 7.15 | 6.89
E-Octatetraene |  Bu($\pi\rightarrow\pi^{\star}$) | 1.93 | 1.92 | 4.32 | 3.69 | 4.22
MAE from Lit. results(Becke, 2018a) . |  |  |  | 0.27 | 0.23 |
ME from Lit. results (Becke, 2018a). |  |  |  |  $-$0.16 |  $-$0.07 |
MAE from TBE-2 results (Silva-Junior et al., 2010a) . |  |  |  | 0.38 | 0.45 |
ME from TBE-2 results (Silva-Junior et al., 2010a). |  |  |  |  $-$0.05 |  $-$0.09 |

aMAE and ME values from RO-PBE0 calculations using 6-31G⋆ are 0.37 and $-$0.22
respectively (Kowalczyk et al., 2013).
---
bMAE and ME values from RO-LC$\omega$PBE0 calculations using 6-31G⋆ are 0.31
and $-$0.05 respectively (Kowalczyk et al., 2013).
†NC denotes “not converged”. ‡ This corresponds to $E_{0\mathrm{S}}$, from
(Becke, 2018a).

At this stage, in order to assess the validity and usefulness, an additional
set of larger molecular systems are taken up in Tables III and IV. We consider
a representative set of $16$ organic chromophores and $5$ linear acene
molecules with 2–6 rings for which accurate theoretical reference data is
available (Silva-Junior et al., 2010b). All the molecular geometry of Table
III are taken from (Silva-Junior et al., 2010b) whereas for Table IV, these
are optimized with B3LYP and cc-pVDZ basis set using GAMESS. Since the present
scheme has earlier been verified to deliver practically identical results as
that of GAMESS, for brevity, we have employed the same for calculation of
$E_{0}$ and $E_{T}$ in the same basis. In addition to B3LYP, another XC
functional from RS hybrid family, has been invoked. Accordingly, the inter-
electronic space is separated through the use of a RS-operator
$g(\gamma,\mathrm{\mathbf{r}})$ such that,

$\frac{1}{\mathrm{\mathbf{r}}}=\frac{\tilde{g}(\gamma,\mathrm{\mathbf{r}})}{\mathrm{\mathbf{r}}}+\frac{g(\gamma,\mathrm{\mathbf{r}})}{\mathrm{\mathbf{r}}},$
(23)

where $\tilde{g}(\gamma,\mathrm{\mathbf{r}})$ represents the complementary RS-
operator. Here, $\gamma$ is an important RS parameter which has a pivotal role
to adjust the contribution of HF exchange between short-range (SR) and long-
range (LR) region, for a specific $g(\gamma,\mathrm{\mathbf{r}})$. In
particular, we consider the long-range correction (LC) scheme presented in
(Iikura et al., 2001), given by,

$\left.\begin{aligned}
E_{\mathrm{xc}}^{\mathrm{LC}}=E^{\mathrm{dfa,sr}}(\gamma)+E^{\mathrm{x,lr}}(\gamma)+E^{\mathrm{c,dfa}},\\\
g(\gamma,\mathrm{\mathbf{r}})=\mathrm{erf}(\gamma\mathrm{\mathbf{r}})\quad\mathrm{and}\quad\tilde{g}(\gamma,\mathrm{\mathbf{r}})=\mathrm{erfc}(\gamma,\mathrm{\mathbf{r}}).\end{aligned}\right\\}$
(24)

It uses full HF exchange at LR region ($E^{\mathrm{x,lr}}$) in contrast to the
conventional hybrid functionals, and successfully mitigates the lack of having
correct asymptotic behaviour of exact exchange at LR region present in hybrid
family. But at the same time, it fails to achieve the standard energy
calculations at the level of hybrid ladder. To be consistent with B3LYP, we
utilize the particular LC-BLYP (with $\gamma=0.33$), where interelectronic
distance dependent B88 exchange (Iikura et al., 2001) (in place of original
B88 exchange (Becke, 1988)) and LYP (Lee et al., 1988) semi-local
correlational functionals were used as $E^{\mathrm{dfa,sr}}(\gamma)$ and
$E^{\mathrm{c,dfa}}$ respectively.

Table 4:  Excitation energies (in eV) of linear acenes from “virial theorem”.

Number of rings | $E_{0\mathrm{T}}$ |  $E_{0\mathrm{S}}$ (PR3) |  Lit.†
---|---|---|---
| B3LYP | LC-BLYP | B3LYP | LC-BLYP |
2 | 3.24 | 3.56 | 4.75 | 5.22 | 4.65
3 | 2.22 | 2.46 | 3.71 | 4.43 | 3.58
4 | 1.55 | 1.81 | 2.85 | 3.41 | 2.75
5 | 1.08 | 1.32 | 2.30 | 2.96 | 2.22
6 | 0.75 | 0.98 | 1.89 | 2.50 | 1.82
MAE from Lit. results (Becke, 2018a) |  |  | 0.09 | 0.70 |
ME from Lit. results (Becke, 2018a) |  |  |  $-$0.10 |  $-$0.70 |
MAE from experimental results (Grimme and Parac, 2003) |  |  | 0.09 | 0.60 |
ME from experimental results (Grimme and Parac, 2003) |  |  |  $-$0.16 |  $-$0.60 |

† This corresponds to $E_{0\mathrm{S}}$, from (Becke, 2018a).
---

The calculated $E_{0\mathrm{S}}$, as computed above, are compared with the
all-electron results from Becke’s “virial theorem” method (Becke, 2018a) along
with respective MAE and ME values for two concerned functionals at the bottom.
As in Table II, here also, due to the proximity of PR3 results to that of
TD-B3LYP, it suffices to report only PR3 excitation energies. For sake of
completeness, this also accompanies the lowest triplet excitation energies.
Note that, for three molecules (pyridine, formamide, imidazole), the spin-
restricted triplet calculations did not lead to convergence in the default
options, for B3LYP. Hence they are not included in MAE and ME evaluations. It
is evident that the performance of B3LYP (excitation energies are
overestimated) is usually much more consistent than that of LC-BLYP. A
systematic error cancellation may be involved in these calculation. On the
other hand, the MAE value of LC-BLYP (0.23) is slightly lower than that of
B3LYP (0.27), but did not outperform extensively as was expected from the
foregoing discussion. Initially, one may draw the conclusion that the effect
of full HF exchange at long-range does not have significant impact on
excitation energies, despite the fact that it enriches the behavior of
frontier orbitals involved in $K_{if}$ calculations. Such discrepancy has
occurred due to the fact that we assume $\gamma$ is independent of system
size. However, it has been shown that it is possible to attain a higher level
of performance provided one treats $\gamma$ as a system-dependent parameter
(functional of $\rho$) tuned from _first principles_ (Baer et al., 2010). It
is our belief that, an optimally tuned $\gamma$ (perhaps in the spirit of size
dependency principle) will provide better results than the conventional hybrid
and RS hybrid functionals. Note that, the less dramatic performance of LC-BLYP
in both ground and triplet excited state has negligible effect on singlet
excitation energy, relative to the STS term, which is in harmony with the
observation made in (Becke, 2018b).  Additionally, we have also estimated MAE
and ME values from TBE-2 results (Silva-Junior et al., 2010a). It is found
that these in Table III (with respect to TBE-2), do not alter significantly
from their counterparts in Becke (Becke, 2018a), which was of course,
expected. So the basic conclusions will remain same. Furthermore, we have also
quoted restricted open-shell KS (ROKS) error statistics (MAE and ME) for
different XC functionals in footnote of Table III from (Kowalczyk et al.,
2013). They are all-electron ROKS lowest vertical excitation results computed
using 6-31G⋆ basis set. Our performance (MAE and ME from TBE-2) are in harmony
with ROKS MAE and ME values.  A similar qualitative trend is also found for
linear acenes, given in Table IV, where again both $E_{0\mathrm{S}}$ as well
as $E_{0\mathrm{T}}$, are given. In this occasion, the better compatibility of
MAE and ME values signify that the error is systematic in nature in case of
LC-BLYP. But these are always overestimated and the extent is larger than
B3LYP. This may arise due to a greater involvement of $\gamma$ on _system
size_ ; hence possibly requires an optimally tuned $\gamma$ in theses cases.
This demands further elaborate investigation, and may be considered in future
works.

##  V Future and Outlook

We have examined the feasibility and practicability of a simple yet accurate
time-independent DFT approach for realistic calculation of single excitation
energies through DFT in CCG. This was applied to a host of organic
chromophores and linear acenes having $\pi$ network. The excitation energies
were estimated using two representative set of functionals–each from
respective hybrid and RS hybrid within a pseudopotential approximation. The
obtained results for all these species, from virial theorem are in quite good
agreement with the reference results. It can be easily extended to the
excitation energies in charge transfer complexes (Becke et al., 2018; Feng et
al., 2018) as well. Its application with RS hybrid, hyper, as well as local
hybrid XC functionals would further enhance the scope in a wide range of
systems. It may also be desirable to examine its performance in various
excited configurations other than the lowest excited state. But this may pose
some computational challenge regarding the state-specific SCF convergence. In
such occasions, some solution, as suggested in (Gilbert et al., 2008; Barca et
al., 2014; Shea and Neuscamman, 2018; Hait and Head-Gordon, 2020; Carter-Fenk
and Herbert, 2020), may offer some valuable guidelines. To conclude, the
present work has demonstrated the usefulness of a simple scheme in predicting
optical gap in organic chromophores within a CCG-DFT.

##  VI Acknowledgement

AG is grateful to UGC for a senior research fellowship. TG acknowledges
INSPIRE program for financial support. AKR thankfully acknowledges funding
from DST SERB, New Delhi, India (sanction order: EMR/2014/000838).

## References

  * Kohn and Sham (1965) W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
  * Becke (2014) A. D. Becke, J. Chem. Phys. 140, 18A301 (2014).
  * Runge and Gross (1984) E. Runge and E. K. U. Gross, Phys. Rev. Lett. 52, 997 (1984).
  * Casida (1995) M. E. Casida, in _Recent Advances in Computational Chemistry; D. P. Chong (Ed.)_ (World Scientific, NJ, 1995), pp. 155–192.
  * Casida and Huix-Rotllant (2012) M. Casida and M. Huix-Rotllant, Annu. Rev. Phys. Chem. 63, 287 (2012).
  * Maitra (2016) N. T. Maitra, J. Chem. Phys. 144, 220901 (2016).
  * Theophilou (1979) A. K. Theophilou, J. Phys. C 12, 5419 (1979).
  * Gross et al. (1988) E. K. U. Gross, L. N. Oliveira, and W. Kohn, Phys. Rev. A 37, 2809 (1988).
  * Singh and Deb (1996) R. Singh and B. M. Deb, J. Chem. Phys. 104, 5892 (1996).
  * Roy et al. (1997) A. K. Roy, R. Singh, and B. M. Deb, J. Phys. B 30, 4763 (1997).
  * Ayers et al. (2015) P. Ayers, M. Levy, and Á. Nagy, J. Chem. Phys. 143, 191101 (2015).
  * Ziegler et al. (1977) T. Ziegler, A. Rauk, and E. J. Baerends, Theor. Chem. Acc. 43, 261 (1977).
  * Kowalczyk et al. (2011) T. Kowalczyk, S. R. Yost, and T. V. Voorhis, J. Chem. Phys. 134, 054128 (2011).
  * Seidu et al. (2014a) I. Seidu, M. Krykunov, and T. Ziegler, J. Phys. Chem. A 119, 5107 (2014a).
  * Seidu et al. (2014b) I. Seidu, M. Krykunov, and T. Ziegler, Mol. Phys. 112, 661 (2014b).
  * Ziegler et al. (2012) T. Ziegler, M. Krykunov, and J. Cullen, J. Chem. Phys. 136, 124107 (2012).
  * Barca et al. (2014) G. M. J. Barca, A. T. B. Gilbert, and P. M. W. Gill, J. Chem. Phys. 141, 111104 (2014).
  * Ramos and Pavanello (2016) P. Ramos and M. Pavanello, Phys. Chem. Chem. Phys. 18, 21172 (2016).
  * Li et al. (2015) C. Li, J. Lu, and W. Yang, J. Chem. Phys. 143, 224110 (2015).
  * Frank et al. (1998) I. Frank, J. Hutter, D. Marx, and M. Parrinello, J. Chem. Phys. 108, 4060 (1998).
  * Kowalczyk et al. (2013) T. Kowalczyk, T. Tsuchimochi, P.-T. Chen, L. Top, and T. Van Voorhis, J. Chem. Phys. 138, 164101 (2013).
  * Shea and Neuscamman (2018) J. A. Shea and E. Neuscamman, J. Chem. Phys. 149, 081101 (2018).
  * Hait and Head-Gordon (2020) D. Hait and M. Head-Gordon, J. Chem. Theory Comput. 16, 1699 (2020).
  * Carter-Fenk and Herbert (2020) K. Carter-Fenk and J. M. Herbert, J. Chem. Theory Comput. (2020).
  * Van Meer et al. (2014) R. Van Meer, O. V. Gritsenko, and E. J. Baerends, J. Chem. Theory Comput. 10, 4432 (2014).
  * Haiduke and Bartlett (2018) R. L. A. Haiduke and R. J. Bartlett, J. Chem. Phys. 149, 131101 (2018).
  * Li Manni et al. (2014) G. Li Manni, R. K. Carlson, S. Luo, D. Ma, J. Olsen, D. G. Truhlar, and L. Gagliardi, J. Chem. Theory Comput. 10, 3669 (2014).
  * Chen et al. (2017) Z. Chen, D. Zhang, Y. Jin, Y. Yang, N. Q. Su, and W. Yang, J. Phys. Chem. Lett. 8, 4479 (2017).
  * Becke (2016) A. D. Becke, J. Chem. Phys. 145, 194107 (2016).
  * Becke (2018a) A. D. Becke, J. Chem. Phys. 148, 044112 (2018a).
  * Becke (2018b) A. D. Becke, J. Chem. Phys. 149, 081102 (2018b).
  * Becke et al. (2018) A. D. Becke, S. G. Dale, and E. R. Johnson, J. Chem. Phys. 148, 211101 (2018).
  * Roy (2008a) A. K. Roy, Int. J. Quant. Chem. 108, 837 (2008a).
  * Roy (2008b) A. K. Roy, Chem. Phys. Lett. 461, 142 (2008b).
  * Roy (2009) A. K. Roy, in _Handbook of Computational Chemistry Research_ , edited by C. T. Collett and C. D. Robson (Nova publishers, New York, 2009).
  * Roy (2010) A. K. Roy, Trends in Phys. Chem. 14, 27 (2010).
  * Roy (2011) A. K. Roy, J. Math. Chem. 49, 1687 (2011).
  * Ghosal and Roy (2016) A. Ghosal and A. K. Roy, _In Specialist Periodical Reports: Chemical Modelling, Applications and Theory; M. Springborg and J.-O. Joswig (Eds.) Vol. 13_ (Royal Society of Chemistry, London, 2016).
  * Ghosal et al. (2018) A. Ghosal, T. Mandal, and A. K. Roy, Int. J. Quant. Chem. 118, e25708 (2018).
  * Mandal et al. (2019) T. Mandal, A. Ghosal, and A. K. Roy, Theor. Chem. Acc. 138, 10 (2019).
  * Ghosal et al. (2019) A. Ghosal, T. Mandal, and A. K. Roy, J. Chem. Phys. 150, 064104 (2019).
  * Becke (1993) A. D. Becke, J. Chem. Phys. 98, 5648 (1993).
  * Schmidt et al. (1993) M. W. Schmidt, K. K. Baldridge, J. A. Boatz, S. T. Elbert, M. S. Gordon, J. H. Hensen, S. Koseki, N. Matsunaga, K. A. Nguyen, S. J. Su, et al., J. Comput. Chem. 14, 1347 (1993).
  * Baer et al. (2010) R. Baer, E. Livshits, and U. Salzner, Annu. Rev. Phys. Chem. 61, 85 (2010).
  * Harris and Jones (1974) J. Harris and R. Jones, J. Phys. F 4, 1170 (1974).
  * J. P. Perdew, R. G. Parr, M. Levy, and J. L. Balduz Jr. (1982) J. P. Perdew, R. G. Parr, M. Levy, and J. L. Balduz Jr., Phys. Rev. Lett. 49, 1691 (1982).
  * Silva-Junior et al. (2010a) M. R. Silva-Junior, M. Schreiber, S. P. Sauer, and W. Thiel, J. Chem. Phys. 133, 174318 (2010a).
  * A. K. Roy, A. Ghosal and T. Mandal (2019) A. K. Roy, A. Ghosal and T. Mandal, _InDFT: A DFT Program for Atoms and Molecules in CCG_ (Theoretical Chemistry Laboratory, IISER Kolkata, India, 2019), This is based on the extension of an initial version of the code, established by A. K. Roy, in 2008, whose results were published in refs. [26-30].
  * Becke and Dickson (1988) A. D. Becke and R. M. Dickson, J. Chem. Phys. 89, 2993 (1988).
  * Stephens et al. (1994) P. J. Stephens, F. J. Devlin, C. F. Chabalowski, and M. J. Frisch, J. Chem. Phys. 98, 11623 (1994).
  * Vosko et al. (1980) S. H. Vosko, L. Wilk, and M. Nusair, Can. J. Phys. 58, 1200 (1980).
  * Becke (1988) A. D. Becke, Phys. Rev. A 38, 3098 (1988).
  * Lee et al. (1988) C. Lee, W. Yang, and R. G. Parr, Phys. Rev. B 37, 785 (1988).
  * rep (2001) _Density Functional Repository, Quantum Chemistry Group_ (CCLRC Daresbury Laboratory, Daresbury, Cheshire, UK, 2001).
  * Stevens et al. (1984) W. J. Stevens, H. Basch, and M. Krauss, J. Chem. Phys. 81, 6026 (1984).
  * Feller (1996) D. Feller, J. Comp. Chem. 17, 1571 (1996).
  * Frigo and Johnson (2005) M. Frigo and S. G. Johnson, Proceedings of the IEEE 93, 216 (2005).
  * Anderson et al. (1999) E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Dongarra, J. D. C. A. Greenbaum, S. Hammarling, A. A. McKenney, and D. Sorensen, _LAPACK Users’ Guide_ , vol. 9 (SIAM, 1999).
  * Johnson (2016) R. D. Johnson, _III (Ed.) NIST Computational Chemistry Comparisons and Benchmark Database, NIST Standard Reference Database, Number, Release 18_ (NIST, Gaithersburg, MD, 2016).
  * Silva-Junior et al. (2010b) M. R. Silva-Junior, M. Schreiber, S. P. Sauer, and W. Thiel, J. Chem. Phys. 133, 174318 (2010b).
  * Iikura et al. (2001) H. Iikura, T. Tsuneda, and K. Hirao, J. Chem. Phys. 115, 3540 (2001).
  * Grimme and Parac (2003) S. Grimme and M. Parac, Chem. Phys. Chem 4, 292 (2003).
  * Feng et al. (2018) X. Feng, A. D. Becke, and E. R. Johnson, J. Chem. Phys. 149, 231101 (2018).
  * Gilbert et al. (2008) A. T. B. Gilbert, N. A. Besley, and P. M. W. Gill, J. Phys. Chem. A 112, 13164 (2008).