I. INTRODUCTION Within density functional theory (DFT), ground-state energy of a many-electron system is conveniently divided into specific components as follows: Here, from left to right, the terms in right-hand side denote kinetic energy of the non- interacting electrons, nuclear-electron attraction, classical electron-electron repulsion, cor- rection to kinetic energy arising out of the interacting nature of electrons, and all non- classical corrections to the electron repulsion energy (including exchange and correlation), respectively. Now, using an orbital expression for density, above equation may be recast as, where N, M denote number of electrons and nuclei respectively, whereas density of a Slater determinantal wave function (exact for a non-interacting system) is simply ρ = ∑ N i 〈 ψ$_{i}$ | ψ$_{i}$ 〉 . As in the Hartree-Fock (HF) case, Kohn-Sham (KS) DFT equations also cannot be solved exactly and recourse must be taken to approximations. However, the required iterative process of solving KS equations is conceptually very similar to that encountered for solving HF equations. Two major routes have been explored for practical solution of a molecular KS equation. The real-space method [1–4] involves an iterative mechanism for a discretized KS equation on a real mesh using either of finite-difference, finite-element or wavelet technique. Typically this whole molecular grid belongs to either uniform or refined uniform grids. Some important advantages of this method are that (a) grid-based matrix representation offers highly structured banded matrices (b) potential operator is diagonal in coordinate space whereas Laplacian operator is nearly local (c) they are easily amenable to the so-called linear scaling methods. Usually these schemes require exceedingly large grid points to capture the complete physical system of interest, which is essential to deliver physically and chemically acceptable results. However with the introduction of higher order and multigrid techniques, grid points could be considerably cut down without sacrificing much accuracy. The other alternative, which is more often used these days, relies on an expansion of eigenfunctions in terms of some suitable basis functions such as Slater or Gaussian type functions (GTF), plane waves, numerical functions, augmented plane waves, linear muffin- 2 E [ ρ ( r )] = N ∑ i ( 〈 ψ$_{i}$ ∣ ∣$^{1}$ 2 ∇ 2 i ∣ ψ$_{i}$ 〉 − 〈 ψ$_{i}$ ∣ M ∑ k Z$_{k}$ | r$_{i}$ − r$_{k}$ | ∣ ∣$^{ψ$_{i}$}$ 〉 + 〈 ψ$_{i}$ ∣ ∣$^{1}$ 2 ∫ ρ ( r $^{′}$) | r$_{i}$ − r $^{′}$| d r ′ ∣ ∣$^{ψ$_{i}$}$ 〉) + E$_{xc}$ [ ρ ( r )] (2) E [ ρ ( r )] = T$_{ni}$ [ ρ ( r )] + V$_{ne}$ [ ρ ( r )] + V$_{ee}$ [ ρ ( r )] + Δ T [ ρ ( r )] + Δ V$_{ee}$ [ ρ ( r )] (1)