The paper presents a method for performing fully self-consistent density-functional calculations that scale linearly with the system size, making it suitable for very large systems. The method uses strictly localized pseudatomic orbitals as basis functions, with sparse Hamiltonian and overlap matrices calculated in \(O(N)\) effort. The long-range self-consistent potential and its matrix elements are computed in real-space grids, while other matrix elements are directly calculated and tabulated based on interatomic distances. The total energy and atomic forces are computed in \(O(N)\) operations using truncated, Wannier-like localized functions and a band-energy functional iteratively minimized without orthogonality constraints. The method is illustrated with examples including carbon and silicon supercells with up to 1000 Si atoms and \(\beta\)-C\(_{3}\)N\(_{4}\) supercells. The authors also apply the method to resolve the controversy about the faceting of large icosahedral fullerenes by performing dynamical simulations on C\(_{60}\), C\(_{240}\), and C\(_{540}\). The method demonstrates linear scaling in CPU time and memory requirements, making it feasible for large-scale ab initio simulations.The paper presents a method for performing fully self-consistent density-functional calculations that scale linearly with the system size, making it suitable for very large systems. The method uses strictly localized pseudatomic orbitals as basis functions, with sparse Hamiltonian and overlap matrices calculated in \(O(N)\) effort. The long-range self-consistent potential and its matrix elements are computed in real-space grids, while other matrix elements are directly calculated and tabulated based on interatomic distances. The total energy and atomic forces are computed in \(O(N)\) operations using truncated, Wannier-like localized functions and a band-energy functional iteratively minimized without orthogonality constraints. The method is illustrated with examples including carbon and silicon supercells with up to 1000 Si atoms and \(\beta\)-C\(_{3}\)N\(_{4}\) supercells. The authors also apply the method to resolve the controversy about the faceting of large icosahedral fullerenes by performing dynamical simulations on C\(_{60}\), C\(_{240}\), and C\(_{540}\). The method demonstrates linear scaling in CPU time and memory requirements, making it feasible for large-scale ab initio simulations.