We have developed and implemented a selfconsistent density functional method using standard norm-conserving pseudopotentials and a flexible, numerical linear combination of atomic orbitals basis set, which includes multiple-zeta and polarization orbitals. Exchange and correlation are treated with the local spin density or generalized gradient approximations. The basis functions and the electron density are projected on a real-space grid, in order to calculate the Hartree and exchange-correlation potentials and matrix elements, with a number of operations that scales linearly with the size of the system. We use a modified energy functional, whose minimization produces orthogonal wavefunctions and the same energy and density as the Kohn-Sham energy functional, without the need for an explicit orthogonalization. Additionally, using localized Wannier-like electron wavefunctions allows the computation time and memory required to minimize the energy to also scale linearly with the size of the system. Forces and stresses are also calculated efficiently and accurately, thus allowing structural relaxation and molecular dynamics simulations.

We present a method to perform fully selfconsistent density-functional calculations, which 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 selfconsistent 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 β-C3N4. We apply the method to solve the existing controversy about the faceting of large icosahedral fullerenes by performing dynamical simulations on C60, C240, and C540.A large effort has been devoted in the last few years to develop approximate methods to solve the electronic structure of large systems with a computational cost proportional to its size. 1 Several approaches are now sufficiently accurate and robust to obtain reliable results for systems with thousands of atoms. So far, however, most of these schemes have been useful only with simple Hamiltonians, like empirical tight-binding models, which provide an ideal setting for order-N calculations. First-principles order-N calculations have been performed mainly in the non-selfconsistent Harris functional version 2 of the local density approximation (LDA) for electronic exchange and correlation (XC) using minimal bases. 1,3 Linear scaling algorithms in fully selfconsistent LDA have also been tried, 4 but the results are far from the linear scaling regime, due to the relatively small number of manageable atoms in those simulations. Hernandez et al. 5 have successfully produced LDA results in large silicon systems using a real-space grid method. The computational requirements that this kind of approach demands are, however, extremely large, and calculations must be performed in massive computational platforms.We have developed a selfconsistent density-functional formulation with linear scaling, capable of producing results for very large systems, whose computational demands are not overwhelmingly large, so that systems with many hundreds of atoms can be treated in modest computational platforms like workstations, and much larger systems can be treated in massive platforms. The method is based on the linear combination of atomic orbitals (LCAO) approximation as basis of expansion of the electronic states. Non-orthogonal LCAO bases are very efficient, reducing the number of variables dramatically, compared to plane-wave (PW) or real-space-grid approaches, so that larger systems can be studied. Also, LCAO can provide up to extremely accurate base...

The performance of basis sets made of numerical atomic orbitals is explored in density-functional calculations of solids and molecules. With the aim of optimizing basis quality while maintaining strict localization of the orbitals, as needed for linear-scaling calculations, several schemes have been tried. The best performance is obtained for the basis sets generated according to a new scheme presented here, a flexibilization of previous proposals. The basis sets are tested versus converged plane-wave calculations on a significant variety of systems, including covalent, ionic and metallic. Satisfactory convergence (deviations significantly smaller than the accuracy of the underlying theory) is obtained for reasonably small basis sizes, with a clear improvement over previous schemes. The transferability of the obtained basis sets is tested in several cases and it is found to be satisfactory as well.

We have implemented a linear scaling, fully self-consistent densityfunctional method for performing first-principles calculations on systems with a large number of atoms, using standard norm-conserving pseudopotentials and flexible linear Ž . combinations of atomic orbitals LCAO basis sets. Exchange and correlation are treated within the local-spin-density or gradient-corrected approximations. The basis functions and the electron density are projected on a real-space grid in order to calculate the Hartree and exchange᎐correlation potentials and matrix elements. We substitute the customary diagonalization procedure by the minimization of a modified energy functional, which gives orthogonal wave functions and the same energy and density as the Kohn᎐Sham energy functional, without the need of an explicit orthogonalization. The additional restriction to a finite range for the electron wave functions allows the Ž . computational effort time and memory to increase only linearly with the size of the system. Forces and stresses are also calculated efficiently and accurately, allowing structural relaxation and molecular dynamics simulations. We present test calculations beginning with small molecules and ending with a piece of DNA. Using double-z, polarized bases, geometries within 1% of experiments are obtained.

A brief review of the Siesta project is presented in the context of linear-scaling density-functional methods for electronic-structure calculations and molecular-dynamics simulations of systems with a large number of atoms. Applications of the method to different systems are reviewed, including carbon nanotubes, gold nanostructures, adsorbates on silicon surfaces, and nucleic acids. Also, progress in atomic-orbital bases adapted to linear-scaling methodology is presented.