Kohn-Sham density functional theory (KSDFT) is the most widely used electronic structure theory for condensed matter systems. The standard method for solving KSDFT requires solving N eigenvectors for an O(N) * O(N) Kohn-Sham  Hamiltonian matrix, with N being the number of electrons in the system.  The computational cost for such procedure is expensive and scales as O(N^3).  We have developed pole expansion plus selected inversion (PEXSI) method, in which KSDFT is solved by evaluating the selected elements of the inverse of a series of sparse symmetric matrices, and the overall algorithm scales at most O(N^2) for all materials including metallic and insulating systems.  The PEXSI method can be used with orthogonal or nonorthogonal basis set, and the electron density, total energy, Helmholtz free energy and atomic force are calculated simultaneously and accurately without using the eigenvalues and eigenvectors.  Combined with atomic orbital basis functions, the PEXSI method can be applied to study the electronic structure of nanotube systems with more than 10,000 atoms with minimum basis set even with on single processor.  The recently developed parallel PEXSI method further allows the accurate treatment of layered graphene-like structure with  more than 20,000 atoms, and can be efficiently parallelized to more than 10,000 processors on leadership class machines.