Self-consistent order-Ndensity-functional calculations for very large systems

We present a method to perform fully self-consistent density-functional calculations that scales linearly with the system size and which is well suited for very large systems. It uses strictly localized pseudoatomic orbitals as basis functions. The sparse Hamiltonian and overlap matrices are calculated with an $O(N)$ effort. The long-range self-consistent potential and its matrix elements are computed in a real-space grid. The other matrix elements are directly calculated and tabulated as a function of the interatomic distances. The computation of the total energy and atomic forces is also done in $O(N)$ operations using truncated, Wannier-like localized functions to describe the occupied states, and a band-energy functional which is iteratively minimized with no orthogonality constraints. We illustrate the method with several examples, including carbon and silicon supercells with up to 1000 Si atoms and supercells of $\ensuremath{\beta}$-${\mathrm{C}}_{3}$${\mathrm{N}}_{4}$. We apply the method to solve the existing controversy about the faceting of large icosahedral fullerenes by performing dynamical simulations on ${\mathrm{C}}_{60}$, ${\mathrm{C}}_{240}$, and ${\mathrm{C}}_{540}$.