This paper presents a selfconsistent density-functional method for large systems that scales linearly with system size. The method uses strictly localized pseudoatomic orbitals as basis functions, allowing for efficient computation of the sparse Hamiltonian and overlap matrices. The long-range selfconsistent potential and its matrix elements are computed on a real-space grid, while other matrix elements are directly calculated and tabulated as functions of 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 that is iteratively minimized without orthogonality constraints. The method is illustrated with examples including carbon and silicon supercells with up to 1000 Si atoms and β-C3N4 supercells. It is applied to resolve the controversy about the faceting of large icosahedral fullerenes by performing dynamical simulations on C60, C240, and C540. The method is efficient and suitable for large systems, with computational demands not overwhelming, allowing systems with hundreds of atoms to be treated on modest computational platforms. The method is based on the linear combination of atomic orbitals (LCAO) approximation, which is efficient compared to plane-wave or real-space-grid approaches. The basis orbitals are strictly localized pseudoatomic orbitals, which give rise to sparse matrices and are efficient for large systems. The method uses standard LDA techniques for valence electrons, with core electrons replaced by pseudopotentials. The basis orbitals are defined by Sankey and Niklewski, and are strictly localized, vanishing outside a given radius. The Kohn-Sham Hamiltonian is obtained by combining techniques to handle the long-range nature of the pseudopotentials. The Hartree potential and exchange-correlation potential are computed using a real-space grid and fast Fourier transform methods. The integrals for these potentials are computed efficiently using sparse-matrix techniques, leading to linear scaling with system size. The method is used to compute the band structure energy and atomic forces, with the band energy minimized using a conjugate gradients algorithm. The total energy is computed using the selfconsistent charge and the Hartree and exchange-correlation potentials. The method is shown to be efficient for large systems, with linear scaling in CPU time and memory requirements. It is applied to study the structure of large icosahedral fullerene clusters, showing that they tend to be polyhedral rather than spherical, with results consistent with previous studies. The method is efficient and suitable for large-scale ab initio simulations.This paper presents a selfconsistent density-functional method for large systems that scales linearly with system size. The method uses strictly localized pseudoatomic orbitals as basis functions, allowing for efficient computation of the sparse Hamiltonian and overlap matrices. The long-range selfconsistent potential and its matrix elements are computed on a real-space grid, while other matrix elements are directly calculated and tabulated as functions of 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 that is iteratively minimized without orthogonality constraints. The method is illustrated with examples including carbon and silicon supercells with up to 1000 Si atoms and β-C3N4 supercells. It is applied to resolve the controversy about the faceting of large icosahedral fullerenes by performing dynamical simulations on C60, C240, and C540. The method is efficient and suitable for large systems, with computational demands not overwhelming, allowing systems with hundreds of atoms to be treated on modest computational platforms. The method is based on the linear combination of atomic orbitals (LCAO) approximation, which is efficient compared to plane-wave or real-space-grid approaches. The basis orbitals are strictly localized pseudoatomic orbitals, which give rise to sparse matrices and are efficient for large systems. The method uses standard LDA techniques for valence electrons, with core electrons replaced by pseudopotentials. The basis orbitals are defined by Sankey and Niklewski, and are strictly localized, vanishing outside a given radius. The Kohn-Sham Hamiltonian is obtained by combining techniques to handle the long-range nature of the pseudopotentials. The Hartree potential and exchange-correlation potential are computed using a real-space grid and fast Fourier transform methods. The integrals for these potentials are computed efficiently using sparse-matrix techniques, leading to linear scaling with system size. The method is used to compute the band structure energy and atomic forces, with the band energy minimized using a conjugate gradients algorithm. The total energy is computed using the selfconsistent charge and the Hartree and exchange-correlation potentials. The method is shown to be efficient for large systems, with linear scaling in CPU time and memory requirements. It is applied to study the structure of large icosahedral fullerene clusters, showing that they tend to be polyhedral rather than spherical, with results consistent with previous studies. The method is efficient and suitable for large-scale ab initio simulations.