Development of accurate efficient methods for the electronic structure of molecules, solids, clusters has been a topic of increasingly intense interest in chemistry, physics, materials sci- ence, etc. In the past two decades, density functional theory (DFT) has been applied with remarkable success to investigate the structural (such as geometries, vibrational frequen- cies, energetics of chemical reactions, etc.) as well as dynamical (such as photoabsorption, photoemission, multi-photon ionization, high-order harmonic generation, etc.) character- istics of such systems. In essence, the ground-state electronic energy is partitioned as: E [ ρ ( r )] = T$_{s}$ [ ρ ( r )] + U [ ρ ( r )] + E$_{xc}$ [ ρ ( r )], where the three terms in the right-hand side denote kinetic energy of a system of non-interacting particles of same density ρ ( r ), classical electro- static energy (including electron-nuclear attraction, electron-electron and nuclear-nuclear repulsion), and exchange-correlation (XC) energy respectively. The single-particle wave functions are obtained from the Kohn-Sham (KS) equation as: where v$_{eff}$ ( r ) contains the relevant one- and two-body potentials, while the electron density is given by the squared sum of the occupied orbitals, ρ ( r ) = ∑ i f$_{i}$ | ψ$_{i}$ ( r ) | $^{2}$. Broadly speaking, two approaches have gained popularity for the practical solution of above KS equation. In the real-space method [1–5] the discretized KS equation is solved iteratively on a mesh using either finite-difference, finite-element or wavelets. The grid-based matrix presentation produces structured, highly banded matrices. Moreover, the near locality (the potential operator is diagonal in coordinate space whereas the Laplacian operator is nearly local) can potentially remove the omnipresent nagging problem arising due to the basis-set incompleteness in the contrasting basis-set approaches. Furthermore, they are easily amenable to the so-called linear-scaling methods. While such schemes usually require a larger number of grids to achieve physically meaningful results, through introduction of higher-order and multigrid techniques, the effective grid can be reduced significantly. In basis-set approaches, on the other hand, the single-particle wave functions are represented by a variety of functions such as Slater-type orbitals, gaussian type functions (GTF), numerical functions, plane waves (PW), augmented plane waves, linear muffin-tin orbitals etc. While currently, there is a preponderance of the atom-centered localized GTFs in the field of molecular quantum chemistry [6] (chiefly due their ability to deal with some of the important 2 I. INTRODUCTION [ 1 2 ∇ 2 + v$_{eff}$ ( r ) ] ψ$_{i}$ ( r ) = ϵ$_{i}$ψ$_{i}$ ( r ) (1)