How to explain to a five year old why "DFT with local exchange–correlation functionals do(es) not describe van der Waals interactions accurately"?

Question

From How do you explain Canada's Trudeau's power-sharing agreement to a five-year-old (American)?:

The "explain to a five-year-old" question is a semi-standard format where the answer demonstrates their ability to understand both the subject matter, and the perspective and limitations of the answer's recipient.

---

Wikipedia's Bilayer graphene begins:

Bilayer graphene can exist in the AB, or Bernal-stacked form,² where half of the atoms lie directly over the center of a hexagon in the lower graphene sheet, and half of the atoms lie over an atom, or, less commonly, in the AA form, in which the layers are exactly aligned.³ In Bernal stacked graphene, twin boundaries are common; transitioning from AB to BA stacking.⁴ Twisted layers, where one layer is rotated relative to the other, have also been extensively studied.

Quantum Monte Carlo methods have been used to calculate the binding energies of AA- and AB-stacked bilayer graphene, which are 11.5(9) and 17.7(9) meV per atom, respectively.⁵ This is consistent with the observation that the AB-stacked structure is more stable than the AA-stacked structure.

²K Yan; H Peng; Y Zhou; H Li; Z Liu (2011). Formation of bilayer Bernal graphene: layer-by-layer epitaxy via chemical vapor deposition. Nano Lett. 11 (3): 1106–10

³Z Liu; K Suenaga PJF Harris; S Iijima (2009). Open and closed edges of graphene layers. Phys. Rev. Lett. 102 (1): 015501.

⁴Min, Lola; Hovden, Robert; Huang, Pinshane; Wojcik, Michal; Muller, David A.; Park, Jiwoong (2012). Twinning and Twisting of Tri- and Bilayer Graphene. Nano Letters. 12 (3): 1609–1615.

⁵E. Mostaani, N. D. Drummond and V. I. Fal'ko (2015). Quantum Monte Carlo Calculation of the Binding Energy of Bilayer Graphene. Phys. Rev. Lett. 115 (11): 115501. also in arXiv

The introduction to ref. 5 Mostaani, Drummond & Fal'ko (2015) states:

Standard first-principles approaches such as density functional theory (DFT) with local exchange–correlation functionals do not describe vdW interactions accurately...

and goes on to discuss some enhancements or alternatives to better model van der Wall interactions for vdW materials like bilayar graphene:

One technique for including vdW interactions in a first-principles framework is to add energies obtained using pairwise interatomic potentials to DFT total energies; this is the so-called DFT-D scheme [1–4]. The development of vdW density functionals (vdW-DFs) that can describe vdW interactions in a seamless fashion is another promising approach [5–8]. DFT-based random-phase approximation (RPA) calculations of the correlation energy [9, 10] provide a more sophisticated method for treating vdW interactions; however, RPA atomization energies are typically overestimated by up to 15% for solids [11, 12], and hence the accuracy of this approach is unclear. Symmetry-adapted perturbation theory based on DFT allows one to calculate the vdW interactions between molecules and hence, by extrapolation, between nanostructures [13]. Finally, empirical interatomic potentials with $r^{-6}$ tails may be used to calculate binding energies [14, 15], although such potentials give a qualitatively incorrect description of the interaction of metallic or πbonded two-dimensional (2D) materials at large separation [16].

But it does not say why "Standard first-principles approaches such as density functional theory (DFT) with local exchange–correlation functionals" fail.

Question: How do you explain to a surface science equivalent five year old why "standard first-principles approaches such as... DFT with local exchange–correlation functionals do not describe van der Waals interactions accurately"?

And if possible within the constraints of a simplified explanation, can one hint at what it is about Quantum Monte Carlo that overcomes this particular limitation?

Answer 1

Standard DFT functionals are based exclusively on local information, e.g. for a GGA functional

$$ E_{xc} = \int n({\bf r}) \epsilon_{xc}(n({\bf r}), \nabla n({\bf r})) {\rm d}^3r $$

while van der Waals interactions are induced dipole - induced dipole interactions which are by definition non-local and have a $r^{-6}$ behavior.

However, van der Waals effects can often be successfully modeled with empirical atom-pair corrections, as in e.g. Grimme's DFT-D3 or DFT-D4 methods.

Answer 2

I'll turn my comment as an answer, as requested by the OP and Nike.

Susi Lehtola's answer has already pointed out that dispersion (which may be a more accurate term than van der Waals interaction here) is by definition a non-local phenomenon, and semilocal functionals are at least unsuitable for describing dispersion. While it's not easy to prove that semilocal functionals can never describe dispersion accurately for Coulombic systems, one can at least show that for some (not necessarily realistic) external potential, semilocal functionals fail to describe dispersion at all. For example consider two helium atoms, each being placed at the center of an infinitely deep potential well. Then move the two atoms (together with their wells) from infinitely far away to a finite distance from each other. Any semilocal density functional will predict a zero energy change, because there is no mechanism for any of the helium atoms to sense whether the other helium atom is near or not: inside the wells none of the density derivatives change due to the incoming helium atom, while outside the wells all density derivatives vanish. However the true energy of course decreases as the helium atoms approach each other, since even strictly confined densities can interact with each other via dispersion. This is a case where semilocal functionals not only fail to describe dispersion accurately, but actually fail to describe it at all.

In order for the functional to describe at least some dispersion in the two strictly confined helium atoms case, one must at least use a double integral over two spatial coordinates, such that every spatial point can know about the density behavior of any other spatial point. For example, a GGA-level dispersion functional (like VV10) has the following form $$ E_c = \frac{1}{2}\iint n(\mathbf{r}) \Phi(\mathbf{r}-\mathbf{r}', n(\mathbf{r}), \nabla n(\mathbf{r}), n(\mathbf{r}'), \nabla n(\mathbf{r}')) n(\mathbf{r}') d\mathbf{r} d\mathbf{r}' \tag{1} $$ where, unlike the single-integral case, $\Phi$ is allowed to depend on the spatial coordinates directly, but only through $\mathbf{r}-\mathbf{r}'$, otherwise one will break translational invariance. In VV10 the functional form is further restricted so that the $\mathbf{r}-\mathbf{r}'$ dependence is only through $|\mathbf{r}-\mathbf{r}'|$, and there are no dot products like $(\mathbf{r}-\mathbf{r}')\cdot\nabla n(\mathbf{r})$ or $\nabla n(\mathbf{r})\cdot\nabla n(\mathbf{r}')$, but one can also imagine functionals without such restrictions. (By the way, the use of such a functional form to describe the whole exchange-correlation energy, instead of just the dispersion energy, seems to be an under-researched direction, and I strongly believe that one can get much better functionals by using a VV10-like functional form to describe exchange and/or correlation.) Interestingly, one can show that even Eq. (1) cannot describe dispersion exactly, because it predicts that the dispersion energy of three non-overlapping densities is exactly the sum of the three pairwise dispersion energies, while in reality this is not the case due to the Axilrod-Teller-Muto three-body dispersion term.

However, there are still a bunch of semilocal functionals that try to describe dispersion without even a form like Eq. (1). While they certainly cannot describe dispersion when the densities of the interacting molecules do not overlap, they may still yield some attractive interaction when the densities overlap. The trick is to notice that if a system is dispersion-bound, this usually means that there are regions where the density is small and the density gradient is even smaller. Only directly contacting but nonbonded atom pairs can create such regions. Therefore, any semilocal functional that gives a sufficiently negative XC energy density ($\epsilon_{xc}$) when $n$ is small and $\nabla n$ is even smaller, but a less negative $\epsilon_{xc}$ when $n$ is small and $\nabla n$ is not too small, can potentially give a dispersion-like attraction; the former requirement stabilizes dispersion-bounded dimers while the latter destabilizes isolated molecules. By tuning the functional behavior in the small-$n$ region, one may even make the functional at least semi-quantitatively reproduce experimental vdW binding energies.

An early example is the mPW family of functionals, where the exchange energy density has a minimum w.r.t. $x=|\nabla n|/n^{4/3}$ when $x$ is large, and becomes less and less negative as $x$ increases beyond this minimum. Because $n^{4/3}$ decays faster than $|\nabla n|$ as one goes into the density tail, the large $x$ behavior of the functional is responsible for the description of weak interactions, and in the case of mPW, this means that molecules can reduce their energy by forming non-covalent complexes, because this reduces the total volume in which $x$ is large: in the middle of these two molecules, $|\nabla n|$ is reduced and $n^{4/3}$ is increased due to the superposition of the molecular densities. The later development of the Minnesota functionals, from the earliest M05/M05-2X and the extremely popular M06-2X to the later M11L, MN15-L etc., incorporate the kinetic energy in their functional form, which provided greater flexibility for fine-tuning the dispersion energy of different systems.

However, this kind of functionals have two inherent problems:

1. The computed dispersion energy depends almost exclusively on the density tails, while the actual dispersion interaction is equally, if not more, dependent on the bulk of the density. This is because dispersion is, to put it loosely, the attraction of instant, fluctuating dipoles, and the dipoles do no need to localize at the periphery of the molecules, but can as well be located deep inside molecules, even near atomic nuclei. Thus, semilocal functionals that try to get dispersion right have to use the information in the density tails to guess the location and intensity of the fluctuating dipoles inside the molecules. This may be possible (due to some holographic principle, for example), but it completely misses the physics.
2. The functional likely does not describe dispersion between next nearest neighbors correctly. This is because the density overlap of a molecule with its next nearest neighbors is obstructed by the density of its nearest neighbors. There is no easy mechanism for a molecule to know whether it lives in a small cluster with only a single solvation layer, or if it lives in bulk solvent.

Finally, I'd like to discuss the question of whether semilocal functionals can reproduce the $r^{-6}$ behavior of the dispersion energy. Although no common functional seems to have demonstrated this, it appears that this is in fact possible. Note that the density decays as $n(r) = \exp(-\alpha r)$ for some $\alpha$ (see my answer to another question). So if we have $\epsilon_{xc} = O(|\ln n|^\beta/n)$ in the small $n$ limit for some number $\beta$, the contribution of the density to the total energy will decay polynomially w.r.t. $r$. This will yield a polynomially decaying interaction energy between molecules, and it is conceivable that we can make the decay behavior be like $r^{-6}$. However we would still suffer from the two problems mentioned above, and now there is a third problem: the $r^{-6}$ behavior will be lost if we use Gaussian basis sets, or any other basis set that does not decay exponentially, and if we use a basis set with a finite spatial extent (like a numerical AO basis set) we lose the dispersion interaction completely.

TL;DR: (1) There exist model systems where all semilocal functionals necessarily yield zero dispersion energy; (2) semilocal functionals can yield an attraction between non-polar molecules, and by empirical tuning one can make it resemble dispersion quantitatively, but even so they get the right answer for the wrong reason; (3) semilocal functionals fail to describe non-nearest-neighbor dispersion, and will be less and less accurate for large systems where non-nearest-neighbor interactions are abundant; (4) basically none of the existing semilocal functionals reproduce the $r^{-6}$ behavior of dispersion (so they fail for molecules that are much farther than their vdW equilibrium distance), and while one may design a new semilocal functional that does so, the functional will be extremely sensitive to the basis set; (5) to describe dispersion accurately while not departing from the Kohn-Sham framework (i.e. neither leaning to the more empirical side and use DFT-D, nor leaning to the more wavefunction-theory side and use double hybrid DFT), one should use non-local functionals such as the VV10 family, which describe two-body dispersion in the correct way, but even they fail to capture at least some of the three-body dispersion effects.