Why does orbital-free DFT scale linearly with system size?

Question

In answering another question here (Is there any relevant DFT formalism apart from the Kohn-Sham approach?), I came across numerous statements that Orbital-Free DFT should scale linearly with system size. From the implementation details in GPAW, this seems to be a result of using the same machinery as for a typical Kohn-Sham DFT calculation, but with only a single "orbital" (actually the square root of the density). But from the other statements, there seems to be a more fundamental physical reason why we would expect orbital free DFT to scale linearly with system size. Why should orbital free DFT be linearly scaling and by what metric for system size (since KS DFT is usually said to scale cubically with the number of orbitals)?

Answer 1

In orbital-free DFT, the key quantity is the particle density (not the density matrix, the actual density). The particle density is a scalar 3D field; when we increase the particle number in our simulation, the values of the particle density change, subject to the constraint that for $N$ particles in a volume $V$ of space,

$$ \iiint_V \rho({\bf r})d^3{\bf r} = N. $$

Although the values change if $N$ changes, the density remains a 3D scalar field and takes up the same amount of computer storage. The size of the density does not depend on the number of particles and so scales as O(1) in particle count -- better than linear scaling!

If we increase the simulation volume instead, then the density now extends over a greater volume and clearly the 3D field must be extended over the additional space. The size of the density scales linearly with the simulation volume, since the volume directly controls how much space we have to represent with our density.

Hopefully you can now see why the size of the density object scales linearly with system size. The size of a data object gives a lower bound for the computational time, assuming that every element of the object has to be checked at least once. However, there is no guarantee that the computational time does scale linearly with system size. It is easy to construct functionals whose computational cost scales far worse than the size of the data object, for example multiple nested integrals over all spatial dimensions. I think the argument for linear-scaling computational time is based on using a semi-local energy functional; even including a Hartree term makes the calculation worse than linear-scaling, strictly speaking, although most people wouldn't complain about the extra factor of $\log(N)$.

Note that whether orbital-free methods scale linearly or not is unrelated to whether the material being simulated is an insulator or a metal. This is in contrast to most "linear-scaling DFT" approaches, which rely on the density matrix decaying exponentially in real-space in order to truncate it safely, which is only valid for systems with a band-gap. Orbital-free methods do not use the density-matrix, so are not subject to the same constraints; in fact the most successful applications of orbital-free DFT are mostly to high-temperature metals (e.g. liquid sodium); even Thomas-Fermi(-Dirac) can perform reasonably in these situations, yet they are amongst the worst cases for density-matrix-based linear-scaling DFT.

Answer 2

There is linear scaling DFT, where you use density matrix as your object throughout SCF iterations instead of the wavefunctions. The linear scaling comes from the sparsity assumption/truncation on the density matrix. Usually, the wavefunctions are dense vectors, and any eigensolver results in cubic scaling with respect to the number of electrons/bases. However, if the density matrix is truncated to a sparse matrix with linear scaling sparsity (w.r.t. number of electrons), then the linear eigensolver in SCF iteration can be implemented through a density filtering method. The main operation in the density filtering method is multiplying the Hamiltonian operator to the density matrix, which scales linearly. Therefore, in the end, the overall cost for the DFT scales linearly.

Regarding the sparsity of density matrix, it can be related to Wannier function, which is localized molecular orbital of crystalline system. The Wannier function is much sparser than the wavefunction (theoretical understanding is available for gaped systems). In DFT calculation, Wannier function can be constructed through rotating (transforming) your wavefunctions. Since Wannier function decays exponentially, it can be well approximated by a truncated function with constant size support. On the other hand, the density matrix (the product of wavefucntions) equals to the product of Wannier functions (two rotation matrices in the middle cancel with each other). Hence the density matrix is sparse and the sparsity scales linearly with respect to the number of electrons.

I am not familiar with OFDFT. Since OFDFT also uses the density matrix as the object, I guess the linear scaling argument is due to the similar reason above.

Answer 3

For linear scaling both Hamiltonian matrix ($H$) and density matrix ($D$) need to be sparse. If either $D$ or $H$ are dense, then the resulting scaling is quadratic $O(N^2)$. When both $H$ and $D$ are dense, then the scaling is cubic $O(N^3)$. As mentioned above the underlying computational cost is determined by matrix-matrix multiplication. For example, to calculate $C = A * B$ we note that one has to find $N^2$ elements of $C$ and each element of $C$ requires $2N$ operations:

$$C(i,j) = \sum_{k=1}^N A(i,k) * B(k,j)$$

The reduction from cubic to quadratic or linear comes from the assumption that $A$ and/or $B$ are not dense and so the summation is reduced to some "fixed" number of elements smaller than N or perhaps the structure of C fixes the number of elements. Hence, not all $N^2$ elements are required to be calculated. For interesting problems of materials science, physics, chemistry, it is very difficult to achieve linear scaling because, depending on the choice of basis set representation, one matrix can be dense which kills linear scaling. Typically, localized orbital basis sets are most promising because they lead to dense Hamiltonian, but for metallic systems the Density matrix becomes delocalized and dense. You can switch to delocalized orbitals in hope to get sparse density matrix but you will get delocalized, dense hamiltonian matrix. In practice the linear scaling is only achieved in the so called "biological systems" (such as a box of waters) in which wide band gap leads to sparse density and hamiltonian matrices. There is a good paper on this topic: Stefan Godecker, Rev. Mod. Phys. 1998.