basis functions; therefore a price is to be paid for efficiency in terms of loss of generality (in contrast to, for example, a plane wave method, where “one basis fits all” philosophy works). Several important factors must be considered for building and choosing basis functions for a particular problem; two most important of them being (i) reduction of number of functions and (ii) ease of computation of relevant integrals. The electron density is described in terms of basis functions and corresponding one-body density matrix P as, where P$_{µν}$ denotes an element of the density matrix. In a spin-unrestricted formalism, ρ ( r ) = ρ $^{α}$( r ) + ρ $^{β}$( r ) , P = P α + P $^{β}$, and KS SCF wave function satisfies the following set of equations, which are reminiscent of Pople-Nesbet equations in HF theory, with the orthonormality conditions, Here C $^{α}$, C β are matrices containing MO coefficients, S is the atomic overlap matrix, and ϵ $^{α}$, ϵβ are diagonal matrices of orbital eigenvalues. F $^{α}$, F β are KS matrices corresponding to α, β spins respectively, having matrix elements as, Here H core µν represents the bare-nucleus Hamiltonian matrix that accounts for one-electron energies, including contributions from kinetic energy plus nuclear-electron attraction. J$_{µν}$ term refers to matrices from classical Coulomb repulsion whereas the third term arises from non-classical XC effects. Last one remains the most difficult and challenging part of the whole SCF process. All one-electron integrals including overlap, kinetic-energy, nuclear-electron attraction as well as pseudopotential matrix elements are identical to those found in HF theory in Cartesian Gaussian functions and are generated by standard recursion algorithms. Here we employ angular-momentum dependent pseudopotentials as those from [29, 30]. Classical Hartree potential is computed by means of a Fourier convolution technique [27, 28], shown to be quite accurate and efficient for molecular modeling. This relies on a Ewald summation 6 F α µν = ∂E$_{KS}$ ∂P α µν = H core µν + J$_{µν}$ + F XCα µν , and F β µν = ∂E$_{KS}$ ∂P β µν = H core µν + J$_{µν}$ + F XCβ µν (8) ( C $^{α}$) $^{†}$SC α = I, and ( C $^{β}$) $^{†}$SC β = I (7) F $^{α}$C α = SC $^{α}$ϵ$^{α}$, and F $^{β}$C β = SC $^{β}$ϵβ (6) ρ ( r ) = N ∑ i =1 K ∑ µ =1 K ∑ ν =1 C$_{µi}$C$_{νi}$χ$_{µ}$ ( r ) χ$_{ν}$ ( r ) = ∑ µ ∑ ν P$_{µν}$ χ$_{µ}$ ( r ) χ$_{ν}$ ( r ) (5)