tin orbitals, etc. Of these, GTFs are the most favorites, for they provide easy analytic routes towards relevant multi-center integrals, A central problem common to almost all DFT approaches is that of obtaining the classical Coulomb ( V$_{c}$ ) and XC ( V$_{xc}$ ) potentials from electron density. In general, these cannot be obtained in any analytic form, and hence numerical treatments are necessary for evaluation of subsequent matrix elements. This is also true for energy integrals associated with XC energy; E$_{xc}$ [ ρ ( r )] = ∫ ϵ$_{xc}$ [ ρ ( r )] d r . In order to circumvent this problem, in some molecular DFT implementations (see, [5, 6], for example), electron density and XC potentials are also expanded in some auxiliary Gaussian bases viz., ρ ( r ) ≃ ˜ ρ ( r ) = ∑ N$_{ρ}$ i C ρ i f ρ i ( r ) , V$_{xc}$ ( r ) ≃ ˜ V$_{xc}$ ( r ) = ∑ N$_{xc}$ j C xc j f xc j ( r ), following some fitting procedure (so called discrete variational method) [7, 8], in addition to an MO expansion. This facilitates an N 3 scaling. A vast majority of modern DFT implementations employ the so-called atom-centered grid (ACG), pioneered by Becke [9], where a molecular grid is conveniently described in terms of some suitable 3D quadratures. The basic idea is to decompose real molecular space into overlapping atomic regions which are described by fuzzy Voronoi polyhedra. These discrete mono-centric atomic integrals can be individually approximated using some stan- dard numerical techniques. Finally, summing these contributions with appropriate weight functions leads to the desired molecular integration result. The atomic grid constitutes of a tensor product between radial part, defined in terms of some quadrature formulas such as Gauss-Chebyshev, Gaussian, Euler-McLaurin, multi-exponential numerical, etc., [10–16] and Lebedev angular quadratures (order as high as 131 has been reported, although usually much lower orders suffice; 59th order is the one most frequently used) [17–19]. Many variants of original Becke integration scheme have been proposed thereafter, mainly to prune away some extraneous grid points which could be quite beneficial. Attempts were also made to perform such integrations by dividing whole space and invoking product Gauss rule to com- plete the subsequent integrations [20]. In a variational integration scheme, on the other hand, molecular space has been categorized into three different regions such as atomic spheres, ex- cluded cubic region and interstitial parallelepiped [21]. In the Fourier transform Coulomb and multiresolution technique, Cartesian coordinate grid (CCG) was used in addition to ACG [22, 23]; former divides Gaussian shell pairs into “smooth” and “sharp” categories on 3 ψ$_{i}$ ( r ) = K ∑ µ =1 C$_{µi}$χ$_{µ}$ ( r ) , i = 1 , 2 , · · · , K (3)