Source-linked AI summary
O(N) methods in electronic structure calculations
D. R. Bowler, T. Miyazaki
TL;DR
Large-scale DFT requires methods that avoid the memory and cubic-cost limits of conventional calculations while retaining accurate electronic-structure modelling. This review explains locality, surveys real-space techniques and seven areas of linear-scaling methods, and discusses their applications and implementations. Methods capable of exact behaviour are reported as more efficient than conventional methods for systems above about a thousand atoms, while major challenges remain.
Problem
Conventional DFT calculations face memory and O(N^3) CPU-time costs, motivating methods with improved scaling for very large systems and massively parallel computers.
Method
The paper reviews electronic-structure locality, real-space methods, and linear-scaling approaches across seven areas, including their technical developments, implementations, and applications.
Results
Methods with the capacity for exact behaviour are more efficient than conventional methods for systems over about a thousand atoms and are beginning to show predictive applications.
Takeaways & Limitations
The review concludes that accurate linear-scaling DFT can address systems of unprecedented size, including emerging calculations involving millions of atoms.
Takeaways & Limitations
Recursive approaches have stable, readily evaluated diagonal Green’s-function elements but difficult and potentially unstable off-diagonal elements needed for force calculations.
Abstract
from arXiv · showhide
Linear scaling methods, or O(N) methods, have computational and memory requirements which scale linearly with the number of atoms in the system, N, in contrast to standard approaches which scale with the cube of the number of atoms. These methods, which rely on the short-ranged nature of electronic structure, will allow accurate, ab initio simulations of systems of unprecedented size. The theory behind the locality of electronic structure is described and related to physical properties of systems to be modelled, along with a survey of recent developments in real-space methods which are important for efficient use of high performance computers. The linear scaling methods proposed to date can be divided into seven different areas, and the applicability, efficiency and advantages of the methods proposed in these areas is then discussed. The applications of linear scaling methods, as well as the implementations available as computer programs, are considered. Finally, the prospects for and the challenges facing linear scaling methods are discussed.
1. Introduction
Conventional plane-wave DFT has enabled increasingly large and accurate calculations, but its computational bottlenecks limit scalability on very large systems and massively parallel computers. The review surveys real-space and O(N) alternatives designed to improve scaling and parallel efficiency.
- Motivation: DFT has expanded from small periodic solids to a broad range of materials and aperiodic systems containing hundreds of atoms.Calculated total energies and atomic forces provide data independently from experiments.
- Conventional DFT: Plane-wave bases provide systematically improvable accuracy and accurate energies and forces without Pulay corrections from atomic-position changes.Pseudopotentials smooth valence-electron wavefunctions, reducing the plane-wave burden while retaining accurate calculations.
- Conventional DFT: FFT reduces the plane-wave transformation cost to N_BN_G ln(N_G), making plane-wave DFT substantially more efficient despite this remaining bottleneck for small systems.Without FFT, the corresponding operation count is of order N_BN_G^2.
- Review scope: The review identifies massively parallel suitability and better-than-cubic, ideally linear, scaling as the two key requirements for extremely large DFT calculations.Real-space methods are considered advantageous for the first requirement, while the paper surveys O(N) methods across seven areas.
- Scaling challenges: Orthogonalisation becomes the dominant cost beyond a few hundred atoms, scaling as O(N^3) when both the number of bands and plane waves grow with N.This bottleneck makes further increases in system size difficult.
- Review scope: The paper reviews recent O(N) progress, localisation, seven classes of linear-scaling methods, implementations, applications, and challenges, with high parallel efficiency as a success criterion.The review first surveys real-space methods and later discusses technical issues for practitioners.
2. Real-space methods
Real-space electronic-structure methods exploit spatial locality to produce sparse matrices, enabling efficient solvers and massive parallelisation. The review covers finite differences, finite elements, local basis functions, and atomic-like orbitals, while highlighting accuracy–cost trade-offs and convergence limitations.
- Real-space approaches: Real-space methods gain efficiency from spatial locality, which produces sparse Hamiltonians and supports massive parallelisation.Finite-difference, finite-element, and local-orbital approaches represent electronic structure locally in real space.
- Finite differences: Finite-difference methods represent wavefunctions on fixed grids, with grid spacing and discretisation order controlling convergence.Higher-order discretisations improve representation but enlarge the operator range and affect efficiency.
- Finite differences: Finite-difference discretisation is not variational because errors can have either sign, making convergence of computational parameters harder.Variational kinetic-energy representations have been proposed to alleviate this problem.
- Finite elements: Finite-element methods use local piecewise-polynomial basis functions and lead to generalised eigenvalue equations because the basis is non-orthogonal.The simulation cell is divided into elements whose shapes need only fill the cell.
- Atomic-like orbitals: Atomic-like orbitals are spatially local and computationally efficient, but their basis-set convergence lacks a simple systematic rule.For Group IV elements, matrix-multiplication costs rise by factors of 11 from four to nine orbitals and 34 from four to thirteen orbitals.
- Electrostatics: Electrostatic treatments include FFT-based Poisson solvers with O(N log N) scaling, Ewald sums with O(N^3/2) scaling, and asymptotically O(N) fast multipole methods.Neutral-atom potentials can remove the need for the ionic Ewald step.
3. Linear Scaling Methods
Linear-scaling methods seek the electronic ground state without solving the most expensive conventional eigenvalue problem. The review motivates density-matrix approaches and notes that orbital-free DFT offers efficiency for large metallic systems at a possible accuracy cost.
- Motivation: Local basis sets can make Hamiltonian construction linear scaling, but solving for eigenstates remains the dominant poor-scaling step.This motivates methods that improve the solution stage rather than only matrix construction.
- Orbital-free DFT: Orbital-free DFT searches over charge densities instead of wavefunctions, using an approximate kinetic-energy functional.The approximation can reduce accuracy, while the approach is used for efficient calculations on large metallic systems.
- Density-matrix formulation: The density matrix is defined from the system’s eigenstates and their occupancies, and it can be expressed as a finite-temperature Fermi function of the Hamiltonian.It is the single-particle, two-point density matrix within the Kohn–Sham framework.
- Density-matrix formulation: Quantum interference makes the density matrix ranged, providing the locality exploited by linear-scaling methods.The detailed functional form of its decay is connected to Wannier-function localisation.
3.1. Density Matrix Properties and Wannier Functions
Wannier functions provide a local representation of electronic structure, and their localisation is closely tied to density-matrix range. The review surveys localisation theory, maximally localised constructions, applications, and unresolved links to disorder, polarisation, and topological materials.
- Wannier functions: Bloch states are delocalised, whereas Wannier functions use phase freedom to construct localised functions.Localised functions are useful for interpreting bonding and electronic structure.
- Physical connections: Wannier-function localisation has been connected to insulating versus metallic behaviour, polarisation, Berry phases, and localised molecular orbitals.The review focuses on aspects relevant to linear-scaling methods rather than covering the full theory.
- Wannier functions: Maximally localised Wannier functions are constructed by minimising a spread functional and have been applied to disordered and entangled-band systems.Efficient iterative approaches have also been developed.
- Applications: Localisation methods are being used to extend DFT, Hartree–Fock, and perturbative quantum-chemistry calculations to improved scaling and accuracy.Examples include linear-scaling exact exchange and localised orbitals for fast MP2 calculations.
- Connection to density matrices: Wannier localisation and density-matrix localisation are intimately related because unitary band transformations leave the density matrix unchanged.The density matrix can therefore be written in terms of Wannier functions, linking their localisation properties.
- Density-matrix range: Density-matrix range varies with band gap, temperature, and disorder, while its relationship to disorder, polarisation, topological insulators, and graphene remains an active area.The review identifies strong effects from defects in graphene as one motivation for further work.
- Density-matrix range: Density-matrix decay is supported as exponential in insulators and in metals at finite electronic temperature, with details depending strongly on the system.For one-dimensional systems, one reported decay prefactor is α = 0.5.
3.2. Solving for the density matrix
Linear-scaling ground-state methods use sparse density matrices represented in localised orbitals. They impose finite range and physical constraints such as electron number and idempotency, while non-orthogonal orbitals introduce additional formal complications and truncation error.
- Core representation: The density matrix is the fundamental quantity in linear-scaling methods, and its sparsity enables linear computational effort after sparse storage and operations are established.A localised basis makes the matrices sparse, but implementing sparse operations requires significant extra effort.
- Truncation: Linear-scaling methods impose a finite density-matrix range as a controllable approximation, which introduces extra error relative to the unconstrained problem.The truncation is an additional constraint on the system.
- Localised-orbital formulation: Ground-state searches commonly represent the density matrix with localised orbitals or support functions and assume a separable form.The approximation assumes only that the original density matrix has finitely many non-zero eigenvalues within the chosen local-orbital representation.
- Constraints: The density matrix must satisfy the correct electron number, with the overlap matrix entering the electron-number expression for non-orthogonal localised orbitals.The review discusses methods for imposing electron number separately.
- Constraints: Idempotency requires the density matrix to project onto the occupied subspace, with eigenvalues restricted to zero or one.Weak idempotency relaxes exact enforcement, while iterative purification can restore exact idempotency.
- Non-orthogonality: Non-orthogonal localised orbitals require either an inverse overlap matrix or an orthogonalisation procedure, both of which complicate sparse calculations.For the occupied subspace, the density matrix can equal the inverse overlap matrix of the local orbitals.
3.2.1. Direct and Iterative Approaches
Direct and iterative approaches construct linear-scaling density matrices or localized orbitals by exploiting idempotency, locality, and sparse matrices. Their methods differ in variationality, stability, convergence, and force-calculation implications.
- Density Matrix Minimisation: Density Matrix Minimisation uses McWeeny purification to impose idempotency during a variational search over an auxiliary density matrix.Sparsity in the density matrix and Hamiltonian enables linear scaling when electron number is handled appropriately.
- Iterative Approaches: Iterative purification methods restore idempotency through repeated transformations, including LNV, trace-resetting, and related sign-matrix approaches.These methods extend across orthogonal and non-orthogonal bases, density functional theory, and finite electronic temperatures.
- Performance and Trade-offs: Purification variants improve efficiency or accuracy in particular filling regimes, but iterative methods are non-variational and complicate force calculations.TRS4 is reported as more efficient and accurate than the original Palser-Manolopoulos method for high- and low-filling problems, with similar efficacy at mid-filling.
- Orbital Minimisation: Orbital Minimisation removes explicit orthogonalisation and can achieve linear scaling by exploiting sparse overlap and Hamiltonian matrices.The original formulation can suffer from many local minima and severe convergence problems; augmented formulations use localized kernel functions to improve stability and convergence.
- DMM and OMM: The DMM and OMM methods share a fundamental connection, while DMM is more commonly implemented because of its simplicity, stability, and single functional minimum.The augmented OMM is presented as a potential improvement, and both methods are among the most commonly used linear-scaling approaches.
3.2.2. Divide and Conquer
Divide-and-conquer methods partition systems into subsystems, solve local electronic-structure problems, and combine their results through shared constraints or density information. They are simple and attractive for large systems, but partitioning and subsystem-size errors limit full-accuracy quantitative use.
- Method: Subsystem solutions are obtained by local diagonalisation, while a system-wide Fermi level enforces the electron-number constraint and feeds charge density back into local Hamiltonians.This produces a self-consistent solution for the entire system after partitioning and choosing local basis sets.
- Solid-State Extensions: Solid-state extensions require buffer atoms around each subsystem to reduce surface effects.Cohesive energy converges to 0.1 eV with 40–50 buffer atoms, whereas the density of states converges much more slowly.
- Method Variants: Density-matrix partitioning offers a refinement in which charge density is constructed from a Mulliken-like assignment of the density matrix to subsystems.This changes the partitioning object from wavefunctions to the density matrix.
- Implementations and Related Methods: Divide-and-conquer implementations exist in OpenMX, SIESTA, hierarchical real-space approaches, and QUAMBO-based methods.The approach has also been applied to quantum-chemistry systems and related large-scale calculations.
- Limitations: The method’s main approximation is projection onto small regions, so results depend on partitioning and subsystem size unless care is taken.Slow convergence with subsystem size and the lack of a variational principle make full-DFT-accuracy quantitative calculations difficult.
3.2.3. Recursive and Stochastic Approaches
Recursive, polynomial-expansion, stochastic, and subspace methods exploit localized Hamiltonian information to construct electronic properties with linear-scaling approximations. Their advantages are balanced by convergence, stability, and variational limitations.
- Recursion methods: Local densities of states can be represented through moments related to powers of a localized-basis Hamiltonian.The pth LDOS moment corresponds to Hamiltonian powers and can be interpreted as hopping around closed lattice loops.
- Recursion methods: Recursion methods repeatedly apply the Hamiltonian, tridiagonalize it, and use a continued-fraction Green’s function to obtain densities of states and density matrices.The procedure starts from a localized state, generates a Krylov subspace, and constructs the required quantities from the resulting tridiagonal Hamiltonian.
- Recursion methods: Off-diagonal Green’s-function elements needed for forces are numerically unstable, motivating Bond Order Potential and block-BOP treatments.Restricting recursion to finite clusters and truncating moments yields linear scaling, although BOP methods are more commonly used between tight-binding and empirical models.
- Polynomial expansions: Chebyshev-based Fermi Operator Expansion constructs the finite-temperature density matrix recursively, becoming linear scaling when matrix-column elements are truncated.Without truncation the approach scales as O(N^2); truncation of localized columns gives linear scaling with atom number.
- Polynomial and stochastic methods: Related polynomial and stochastic approaches reconstruct the density of states or thermodynamic quantities using weighted polynomials, random vectors, and sparse-matrix inversions.The Kernel Polynomial Method reduces finite-order oscillations with Gibbs factors, while stochastic methods estimate inverses through random-vector linear solves.
- Renormalisation and subspace methods: Energy renormalisation expands the density matrix into progressively lower-temperature corrections, while subspace methods use Krylov vectors and shifted Green’s-function solves.The first energy-renormalisation term is short-ranged, but the approach had not demonstrated true linear scaling; subspace methods can reuse vectors across energies.
- Assessment: FOE-like methods have achieved the greatest ab initio recursion-method success, whereas convergence and limited variational properties reduce the success of other approaches.These limitations are identified as reasons recursion methods have been less successful than competing methods.
3.2.4. Penalty functionals
Penalty-functional methods optimize a density matrix while enforcing idempotency through an energetic penalty and controlling electron number with a chemical potential. Their performance depends on the penalty parameter and on avoiding minimization pathologies.
- Penalty-functional formulation: The penalty-functional formulation augments the Kohn–Sham energy with chemical-potential and idempotency-penalty terms.For an idempotent density matrix, the penalty vanishes; the chemical potential sets the electron number, leaving the penalty coefficient as the tunable parameter.
- Penalty-functional formulation: For α greater than a critical value αC, the correct idempotent ground-state density matrix is obtained, but αC cannot be predicted exactly.Too-small α can produce local minima, whereas too-large α slows convergence.
- Limitations and corrections: A square-root penalty creates a branch point at the minimum, preventing standard conjugate-gradient minimization from being effective.A corrected functional was proposed to remove this difficulty.
- Limitations and corrections: The corrected approach does not impose idempotency exactly and introduces occupancy and total-energy errors dependent on α before correction.Occupied bands can have more than one electron and unoccupied bands negative occupancies; following correction, the total energy becomes α-independent.
- Implementation and adoption: Penalty-functional methods appear in ONETEP as one stage in a sequence of density-matrix optimization methods but have not been widely adopted.The sequence proceeds from canonical purification through penalty functionals and finally LNV.
3.2.5. Orbital-free DFT
Orbital-free DFT seeks the ground-state charge density without calculating orbitals, simplifying large-system calculations. Its central challenges are constructing accurate kinetic-energy functionals and using sufficiently transferable local pseudopotentials.
- Concept and functionals: Thomas–Fermi and von Weizsäcker functionals represent limiting behaviors for uniform and rapidly varying charge densities, respectively.Early orbital-free approaches combined the two functionals into a single form.
- Limitations: Because charge density is local, orbital-free DFT requires local pseudopotentials, restricting transferability and accuracy.This pseudopotential constraint is identified as a significant problem for the method.
- Developments: Subsequent work introduced density-dependent kernels, non-local functionals, finite-element and coarse-grained schemes, and extensions to covalent systems and semiconductors.A key advance was a non-local kinetic-energy functional with parameters transferable across environments having similar coordination numbers.
- Assessment: Orbital-free DFT provides a good route for modeling large metallic systems, but its kinetic-energy functional remains inexact and its pseudopotentials remain limited.The paper identifies detailed accuracy and convergence comparisons with recursion methods and other approaches as valuable.
3.2.6. Expansion of the density matrix and tensorial approaches
Density-matrix and tensorial approaches use parameterizations, unitary transformations, and manifold-based optimization to preserve electronic constraints while seeking linear scaling. Their promise is tempered by unresolved dense-algebra costs and limited adoption.
- Exponential and unitary parameterizations: Exponential parameterizations impose density-matrix idempotency or orbital orthogonalization implicitly, providing routes to conventional optimization and possible linear scaling.The first linear-scaling proposal used an exponential density-matrix parameterization.
- Exponential and unitary parameterizations: The exponential density-matrix transformation preserves idempotency and trace, and its Baker–Campbell–Hausdorff expansion supports linear-scaling Hartree–Fock and DFT formulations.Related orthogonal-basis curvy-step methods and extended-Lagrangian applications were also reported.
- Curvy-step optimization: Curvy-step optimization uses unitary transformations to preserve idempotency and electron number while expanding the energy as a polynomial search.The formulation emphasizes tensorial correctness in non-orthogonal bases.
- Grassmann-manifold methods: Grassmann-manifold optimization directly minimizes wavefunction coefficients while retaining the orbital subspace on which the energy depends.The method represents basis-set coefficients as a matrix transformation of the bands.
- Sparsity methods: Sparsity-oriented unitary transformations maximize coefficient sparsity with an ℓ1-norm, but the initial formulation defined only steepest-descent gradients.The approach was intended for linear scaling and was described as promising.
- Sparsity methods: CP2K orbital transformation and iterative refinement replace constrained orthogonality optimization with locally unconstrained variables and matrix-function iterations.The proposed fourth-order refinement is paired with Taylor expansions rather than diagonalization to target linear scaling.
- Limitations and outlook: Dense-algebra preconditioners still impose O(N^3) scaling in the orbital-transformation method, and the broader class is not yet widely used.The paper suggests that the dense-preconditioner bottleneck may eventually be removed.
3.2.7. Quantum Chemistry
Linear-scaling quantum-chemistry methods exploit locality to reduce the cost of correlation calculations, extending approaches from MP2 toward coupled-cluster and configuration-interaction methods.
- Local correlation methods reduce scaling, with Laplace-transformed MP2 extended to atomic orbitals to achieve asymptotic N^2 scaling.
- Linear-scaling formulations have been demonstrated for CCSD using fragment-confined non-orthogonal orbitals or local divide-and-conquer bases, with dynamic localisation for error control.
- MRSD-CI has also achieved linear scaling through local correlation and integral screening.
- Divide-and-conquer extensions to MP2 and CCSD use full Hartree–Fock subsystem orbitals, while calculations on many tens of atoms are currently possible.
- Beyond ground states, related work is extending linear-scaling approaches toward responses and excitations, although the survey is not comprehensive.
3.2.8. Extensions
Extensions of linear-scaling electronic-structure methods address time-dependent response, excitation energies, polarizability, and band-edge calculations using density-matrix and recursion-based approaches.
- Real-time TDDFT can propagate the density matrix rather than wavefunctions, with time variation governed by its commutator with the Hamiltonian.
- Truncation effects on propagation accuracy have not yet been comprehensively investigated, and approximately 35 fs of propagation is needed for 0.1 eV energy resolution in medium-to-large systems.
- Linear-response TDDFT uses Liouvillean eigenvalues for vertical excitation energies, with Krylov, variational, and recursion approaches avoiding explicit virtual-orbital representations.
- Density-matrix perturbation theory recursively generates perturbed density matrices and has been applied to polarizability, including water clusters containing 150 molecules.
- Laplace-transform and projection-based reformulations improve efficiency for polarizability calculations, with applications including NMR and competitive performance against other linear-scaling methods.
- Band edges can be obtained from the density matrix by seeking extreme eigenvalues associated with the highest occupied and lowest unoccupied states.
4. Technical Details & Parallelisation
Technical developments address non-orthogonal representations, sparse overlap-matrix inversion, electron-number preservation, and parallel or domain-decomposed calculations, while structural relaxation remains a distinct scaling challenge.
- 4.1. Non-orthogonal Basis Functions: Non-orthogonal localised orbitals are more contracted and can provide computational advantages, but require careful covariant, contravariant, and mixed tensor representations.
- 4.1. Non-orthogonal Basis Functions: The inverse overlap matrix is exponentially localised when the overlap is localised, but its range depends on the overlap spectrum and its condition number affects inversion difficulty.
- 4.1. Non-orthogonal Basis Functions: Sparse linear-scaling inversion methods include iterative Hotelling–Schultz updates, conjugate-gradient formulations, divide-and-conquer clusters, and sparsity-constrained minimisation.
- 4.1. Non-orthogonal Basis Functions: Cubic-scaling overlap inversion may have a small prefactor but is not scalable, while iterative inversion can stall when truncation error matches the per-step inverse change.
- 4.2. Preserving Electron Number: The electron-number-preserving functional keeps the ground-state search from perturbing electron number, requiring only the chemical potential to be updated from a correct initial electron count.
- 4.2. Preserving Electron Number: A transformed formulation can introduce multiple minima, although reported results indicate that convergence is not adversely affected.
- 4.4. Structure Relaxation: Linear-scaling electronic ground states do not guarantee equally scalable atomic relaxation, especially for systems with low and complex curvature.
5. Implementation and Applications
O(N) DFT methods are being applied to increasingly large biological and nanoscale systems, with studies demonstrating useful agreement, convergence, and experimentally consistent structural energetics. Implementations also expose practical limitations, including functional, basis-set, and scaling constraints.
- 5.1. Implementations: The applicability of O(N) implementations is constrained by method-specific limitations: PROFESS permits only local pseudopotentials and restricted functionals, while basis sets lack systematic convergence.For the cited LNV implementation, a local-pseudopotential contribution becomes significant around 15,000 atoms despite its small prefactor.
- 5.2.1. O(N) calculations on biological systems: O(N) DFT studies target large biological systems because conventional force fields can depend strongly on interatomic-potential parameters.Linear scaling methods are expected to enable DFT calculations on entire or sufficiently large complex biomolecular regions.
- 5.2.1. O(N) calculations on biological systems: A 715-atom dry-DNA calculation found orbital cleavage at an introduced mutation, consistent with high resistivity measured for λ-DNA.The study examined eleven guanine-cytosine base pairs and compared orbitals before and after swapping one base pair.
- 5.2.1. O(N) calculations on biological systems: For a 3,439-atom hydrated DNA system, DMM accuracy was reported as surprisingly good, with total energy examined against density-matrix cutoff length.The corresponding dry system contained 643 atoms, and calculations used a minimal PAO basis in a non-self-consistent setting.
- 5.2.1. O(N) calculations on biological systems: A 24,000-atom influenza haemagglutinin calculation found that sialoside binding was not regulated by allosteric effects.The FMO approach used a polarisable continuum model and supports interaction analysis through pair interaction energy decomposition.
- 5.2.2. Order-N DFT study on nanoscale structures of Ge islands on Si(001): For Ge islands on Si(001), structural optimisation remained robust for a system of approximately 23,000 atoms, while 3D huts became more stable above 2.7 ML coverage.The calculated 2D-to-3D transition coverage agreed with the experimental value, although kinetic aspects remained important for simulating growth.
- 5.2.3. Other examples: Ab initio O(N) molecular dynamics reproduced several experimental properties of liquid ethanol, including a self-diffusion coefficient of 8.2×10^-6 cm^2s^-1 versus 1.1×10^-5 cm^2s^-1.The radial distribution function and hydrogen-bond-induced O-H vibrational red shift also compared well with experiment.
6. Conclusions
Linear-scaling methods have progressed from theoretical proposals to practical, predictive calculations for systems beyond about a thousand atoms, including the first true DFT calculations on millions of atoms. Important challenges remain in accuracy, metallic systems, functionality, applications, and adapting to emerging computer architectures.
- Linear-scaling methods introduce additional parameters and instability sources, while their convergence criteria and basis-set choices remain less established than in standard codes.
- Methods with the capacity for exact behaviour are more efficient than conventional methods above about 1,000 atoms and are demonstrating predictive applications.
- Low- or zero-gap metallic systems with extended electronic structure still lack a clear linear-scaling solution of sufficient quantitative accuracy.
- Real-space methods suit multicore and GPU architectures, but communication patterns designed for earlier systems may not transfer efficiently to novel architectures.
- Future development must improve functionality, including exact exchange and dispersion forces, while preserving linear scaling and efficient parallelisation.
- The first true DFT calculations on millions of atoms have been achieved, with fully converged highly accurate results expected soon.