# Signatures of van der Waals binding: A coupling-constant scaling analysis

## Abstract

The van der Waals (vdW) density functional (vdW-DF) method [ROPP 78, 066501 (2015)] describes dispersion or vdW binding by tracking the effects of an electrodynamic coupling among pairs of electrons and their associated exchange-correlation holes. This is done in a nonlocal-correlation energy term , which permits density functional theory calculation in the Kohn-Sham scheme. However, to map the nature of vdW forces in a fully interacting materials system, it is necessary to also account for associated kinetic-correlation energy effects. Here, we present a coupling-constant scaling analysis which permits us to compute the kinetic-correlation energy that is specific to the vdW-DF account of nonlocal correlations. We thus provide a more complete spatially-resolved analysis of the electrodynamical-coupling nature of nonlocal-correlation binding, including vdW attraction, in both covalently and non-covalently bonded systems. We find that kinetic-correlation energy effects play a significant role in the account of vdW or dispersion interactions among molecules. Furthermore, our mapping shows that the total nonlocal-correlation binding is concentrated to pockets in the sparse electron distribution located between the material fragments.

## I Introduction

Many-body effects are essential for an accurate description of materials bonds.
Such correlation effects must be accurately reflected in density functional
theory (DFT) as we seek approximate evaluations of ground-state expectation values,
and , of operators for the
kinetic energy and for the electron-electron interaction . This is
clear, for example, because dispersion or van der Waals (vdW) interactions arise from an
electrodynamical coupling among collective excitations.Mahan (1965); Langreth and Perdew (1977); Rapcewicz and Ashcroft (1991); Maggs and Ashcroft (1987); Langreth and Vosko (1987); Hyldgaard *et al.* (2014)
In the Kohn-Sham (KS) schemeKohn and Sham (1965) for efficient, in principle, exact DFT calculations,
we handle all many-body effects by a trick. We focus on an independent-particle approximation,
the KS kinetic-energy term , while embedding the difference, termed the
kinetic-correlation energy,

(1) |

in the exchange-correlation (XC) energy, . However, the
kinetic-correlation energy, Eq. (1), can still
be unmasked as a functional of the electron density
using a formally exact scaling analysis.Levy and Perdew (1985); Levy (1991); Görling and Levy (1993); Burke *et al.* (1997)
An evaluation of is equivalent to correcting
beyond the Hartree approximationLevy and Perdew (1985) and thus allows an exploration
of many-electron interaction effects.

The van der Waals (vdW) density functional (vdW-DF) method for general-purpose
DFT calculations relies on truly nonlocal formulations , Refs. Andersson *et al.*, 1996; Rydberg *et al.*, 2000, 2003; Dion *et al.*, 2004; Thonhauser *et al.*, 2007; Lee *et al.*, 2010; Berland and Hyldgaard, 2014; Thonhauser *et al.*, 2015; Berland *et al.*, 2015. The vdW-DF functional design can be seen as a systematic
extension of the local density approximations (LDA) and of the generalized gradient
approximation (GGA). In its original and most commonly used form,Dion *et al.* (2004); Berland and Hyldgaard (2014)
it relies on the same many-body perturbation theory analysisLangreth and Vosko (1987); Thonhauser *et al.* (2007)
that underpins the formulations of PBEPerdew *et al.* (1996a) and PBEsol,Perdew *et al.* (2008) and it
adheres to the same fundamental principle, that physics constraints, including
chargeDion *et al.* (2004) and currentDion *et al.* (2004); Hyldgaard *et al.* (2014) conservation, should guide the XC
functional design.Rydberg *et al.* (2003); Rydberg (2001); Berland *et al.* (2014)
However, as part of a rationale for constraint-based GGA, Langreth
and Vosko showed that a gradient-corrected formulation of correlation cannot
naturally account for vdW interactions.Langreth and Vosko (1987) The vdW-DF method overcomes
that limitation, noting that electrons and their associated GGA-type XC holes
themselves form dipole systems with internal dynamics.Maggs and Ashcroft (1987); Rapcewicz and Ashcroft (1991); Andersson *et al.* (1996); Becke and Johnson (2005, 2007); Hyldgaard *et al.* (2014)
The vdW-DF method tracks screening effects produced by the mutual electrodynamical
coupling of such virtual dipoles. It thus extends GGA within the vdW-DF framework,
capturing screened dispersion bindingMaggs and Ashcroft (1987); Rapcewicz and Ashcroft (1991) by summing coupling-induced shifts in the
collective plasmon excitations.Mahan (1965); Hyldgaard *et al.* (2014)

The vdW-DF method is computationally efficient since the dispersion-energy
gainsHyldgaard *et al.* (2014) are evaluated in a truly nonlocal-correlation energy
term that is an explicit functional of the density.Dion *et al.* (2004); Román-Pérez and Soler (2009)
This is done by using the adiabatic-connection formulaLangreth and Perdew (1975); Gunnarsson and Lundqvist (1976); Langreth and Perdew (1977)
(ACF) for the exact XC functional to define an effective dielectric
function ,Dion *et al.* (2004); Hyldgaard *et al.* (2014); Berland *et al.* (2015) and by
expanding in terms of a plasmon-pole approximation that reflects
the response corresponding to an internal semilocal functionalRydberg (2001); Dion *et al.* (2004); Lee *et al.* (2010); Hyldgaard *et al.* (2014)
. In the original general-geometry vdW-DFDion *et al.* (2004)
and in the recent consistent-exchange vdW-DF-cxBerland and Hyldgaard (2014)
formulations,^{1}*et al.* (2004); Thonhauser *et al.* (2007, 2015)
The total functional specificationDion *et al.* (2004); Berland and Hyldgaard (2014)

(2) |

generally also contains a cross-over term that contains nothing but
gradient-corrected exchange.Berland *et al.* (2014, 2015) The
total exchange functional is semilocal; the correlation part of the functional,
, comprises LDA correlation (from ) and
.

An elegant illustration of the many-body physics nature of vdW binding
can be obtained by computing the spatially resolved componentCallsen *et al.* (2012); Lazić *et al.* (2012)

(3) |

of the total vdW-DF nonlocal-correlation energy

(4) |

The spatial resolution Eq. (4) is a natural extension
of how we normally analyze total-energy contributions arising from the
semilocal components of the XC energy.Langreth and Perdew (1977); Perdew *et al.* (1996a); Burke *et al.* (1997); Dion *et al.* (2004)
The spatially resolved energy Eq. (3) is
given by the vdW-DF
kernelDion *et al.* (2004, 2005) for which there exist both formal
analysisThonhauser *et al.* (2007) and an efficient evaluation scheme.Román-Pérez and Soler (2009)
With Eq. (3) one can track and understand
binding-induced changes ,
for example, for benzene adsorption on graphene.Lazić *et al.* (2012)
The mapping confirms that the dominant contributions to the vdW binding arise
in the regions of sparseLangreth *et al.* (2009) (but not vanishingly low) electron
density between molecules and surfaces.Rydberg *et al.* (2003); Kleis *et al.* (2008); Berland and Hyldgaard (2010, 2013); Hyldgaard *et al.* (2014)

In this paper, we seek a characterization of electron-electron interaction effects
that underpin vdW attraction between molecules. Many computational descriptions of
the vdW attraction build on a discussion of the electron response and dielectric
function in the physical, fully interacting system,Lifshitz (1955); Mahan (1965); Zaremba and Kohn (1976, 1977); Nordlander and Harris (1984); Hult *et al.* (2001); Rapcewicz and Ashcroft (1991); Maggs and Ashcroft (1987); Andersson *et al.* (1996); Rydberg *et al.* (2000, 2003); Kleis *et al.* (2007); Grimme (2004, 2006); Grimme *et al.* (2010); Becke and Johnson (2005, 2007); Silvestrelli (2008); Tkatchenko and Scheffler (2009); Vydrov and Van
Voorhis (2010); Ruiz *et al.* (2012); Tkatchenko *et al.* (2012) although the actual
response behavior is sometimes approximated by an independent-particle description. The ACF specifies the
exact XC functional as an average over the electron response, denoted , that reflects a
ramping () of assumed electron-electron interaction strengths, .Langreth and Perdew (1975); Gunnarsson and Lundqvist (1976); Langreth and Perdew (1977); Perdew *et al.* (1996b); Burke *et al.* (1997) This ACF view is explicitly maintained in vdW
density functionalsLundqvist *et al.* (1995); Dobson and Dinte (1996); Dion *et al.* (2004); Berland *et al.* (2014, 2015); Thonhauser *et al.* (2015)
which track the vdW binding produced by plasmon-energy shiftsMahan (1965); Rapcewicz and Ashcroft (1991) in a KS framework.Langreth and Perdew (1977); Hyldgaard *et al.* (2014)
However, the electrodynamical-coupling mechanism for vdW attractionMaggs and Ashcroft (1987); Rapcewicz and Ashcroft (1991) is at work in the physical
system, i.e., at full coupling-constant strength . For a more complete mapping of the
nature of vdW attraction,Rapcewicz and Ashcroft (1991) we therefore seek to (a) compute an XC energy, denoted
, that instead reflects the physical response , and
(b) extract and study the component, denoted ,
that corresponds to nonlocal-correlation effects in .

Our central observation is that such information is directly available from the vdW-DF
functional form, Eq. (2), by applying the formally exact coupling-constant scaling
analysisLevy and Perdew (1985); Levy (1991, 1995) on the vdW-DF method. The formal analysis rests on density scaling, which provides a complete specification of the would-be XC energy that reflects the response function assuming only that the -averaged response defines the specific
form; the analysis can be made for a given problem once we know the self-consistent solution density .
We present details of how to extend the scaling analysis from semilocal functionalsGörling and Levy (1993); Perdew *et al.* (1996c); Burke *et al.* (1997); Ernzerhof *et al.* (1997) to the truly nonlocal-correlation
term of the vdW-DF method.

We note that the formal scaling analysis permits calculations of the kinetic-correlation energy,Levy and Perdew (1985) Eq. (1). For practical calculations, we present a code, termed ppACF, that computes the component

(5) |

which is specific to . Eq. (5) is also combined with the known coupling-constant scaling analysis for LDA correlation,Levy and Perdew (1985); Levy (1991); Görling and Levy (1993); Levy (1995) for a full specification of the kinetic-correlation energy

(6) |

Finally we rely on the formal equivalenceLevy and Perdew (1985)

(7) |

to extract a representation

(8) |

of the mutual plasmon electrodynamical coupling in the physical
system.Mahan (1965); Rapcewicz and Ashcroft (1991); Maggs and Ashcroft (1987); Hyldgaard *et al.* (2014)

As implied in Eqs. (5) and (6), the code also gives us
access to spatially resolved kinetic-correlation energies,
and , that are consistent with Eqs. (3) and (4)
and with the standard resolution of XC energy contributions. Using ppACF, we can thus
compute and discuss the nature of binding-induced changes
, ,
and ,
in the spatially resolved descriptions. Our ppACF code can provide
this analysis for most versions or variants of the vdW-DF method.Dion *et al.* (2004); Lee *et al.* (2010); Cooper (2010); Klimeš *et al.* (2010, 2011); Berland and Hyldgaard (2014); Hamada (2014); Berland *et al.* (2015)
Here we work with the consistent-exchange vdW-DF-cx formulation,Berland and Hyldgaard (2014); Berland *et al.* (2014)
which can effectively be seen as a mean-value evaluation of the ACF.Hyldgaard *et al.* (2014)

We find that and both
contain signatures of directed binding: the dominant binding contributions are
channeled into pockets. Also, the signatures in and
in typically mirror each other,
up to a sign. This means that the concentration of vdW bonding is
further enhanced in the contribution that characterizes the electrodynamical
coupling mechanism behind the vdW attraction.Mahan (1965); Maggs and Ashcroft (1987); Rapcewicz and Ashcroft (1991); Hyldgaard *et al.* (2014)

Overall, our results show that there is an important kinetic-energy nature of
vdW binding and confirm that the density tails, rather than the atomic centers,
play the decisive role in setting dispersion forces at binding separations.Rydberg *et al.* (2003); Kleis *et al.* (2008); Berland and Hyldgaard (2010); Berland *et al.* (2011); Lee *et al.* (2011, 2012a); Lazić *et al.* (2012); Berland and Hyldgaard (2013, 2014); Hyldgaard *et al.* (2014) Our results also
suggest that there exists an orbital-like structure of dispersion binding, although
much weaker than in chemical bonds and originating in different mechanisms.Maggs and Ashcroft (1987); Thonhauser *et al.* (2007); Hyldgaard *et al.* (2014)
This observation could be useful for qualitative discussions of the nature and
variation in vdW forces in materials.

The rest of this paper is organized as follows. Section II details the coupling constant analysis of the vdW-DF method. Section III provides computational details. In Sec. IV we document signatures of the vdW attraction in both noncovalent and covalent molecular binding. Section V contains a summary and discussion. The paper has one appendix.

## Ii Theory

A systematic theory characterization of the screened response in a homogeneous and weakly perturbed electron gasLangreth and Perdew (1975); Gunnarsson and Lundqvist (1976); Langreth and Perdew (1977, 1980); Langreth and Mehl (1981); Langreth and Vosko (1987, 1990) has led to the definition of a range of successful constraint-based functionals for the XC energy and broad use of DFT. We use to denote the full electron-electron interaction. We consider the density changes produced by an external field , and compute the electron-gas density response as a function of the assumed coupling constant for an adiabatic turn on of the many particle interaction, . The exact XC energy is given by the ACF,

(9) |

which links , the (complex) frequency , and spatial variations in the response function to the XC energy. We use to denote the density operator, and the last term of Eq. (9) is the electron self energy

The exact XC energy can be recast as an electrostatic interactionLangreth and Perdew (1975); Gunnarsson and Lundqvist (1976); Langreth and Perdew (1977)

(10) |

between the electrons and associated, so-called, XC holes .
The XC hole reflects a average of the response . An emphasis
on the assumed plasmon-nature of the electron-response, a reliance on formal many-body perturbation
theory, and the imposing of additional physics constraints, such as charge conservation
of the XC hole, has led to formulations of LDA,Hedin and Lundqvist (1971); Perdew and Wang (1992)
of the PBE and PBEsol versions of GGAs,Perdew *et al.* (1996a, 2008) and of
the vdW-DF method.Dion *et al.* (2004); Thonhauser *et al.* (2007); Lee *et al.* (2010); Berland *et al.* (2014); Hyldgaard *et al.* (2014); Berland *et al.* (2015)

At any given coupling constant , the response function defines an approximation for the exchange-correlation hole

(11) |

The actual XC hole then emerges simply as an average,

(12) |

Using Eq. (11) it is meaningful to define and discuss also the coupling-constant dependence of the XC functional:

(13) |

Same as for the holes, the actual functional, Eq. (10), is given by an average over ,

(14) |

The behavior of is exclusively set by exchange effects at . This follows because exchange reflects
an independent-particle behavior and, unlike correlation, it is
independent of . At the other physical limit, the plasmon character can be expected to dominate in the response. One therefore also expects that
becomes accurate at if Eq. (13) reflects a plasmon-based analysis of electron response, for example, as used the early LDA formulations Hedin and Lundqvist (1971); Gunnarsson and Lundqvist (1976), in the constraint-based GGAs Perdew *et al.* (1996a, 2008), and in vdW-DF-cx Berland and Hyldgaard (2014).

It is instructive to split the XC hole into exchange and correlation components

(15) |

and to define (at every ) a spatially resolved correlation term

(16) |

This term provides a mapping of the total correlation effects at :

(17) |

Also, there exists a coupling-constant scaling
analysis for LDA correlationLevy (1995); Görling and Levy (1993); Burke *et al.* (1997)
with spatial resolution

(18) |

Accordingly we isolate a spatially resolved nonlocal-correlation energy

(19) |

corresponding the coupling-constant scaling of the total nonlocal-correlation energy

(20) |

Equation (3) is the coupling-constant integral of .

To map the electrodynamical-coupling nature of vdW attraction, we seek to compute binding-induced changes .

### ii.1 Density scaling in the exact XC energy

Coupling-constant scaling analysisLevy and Perdew (1985) is a natural
tool for exploring the nature of both exchange-based GGAsBurke *et al.* (1997); Ernzerhof *et al.* (1997)
and of vdW-DF-cx. For any given solution density , we define
a rescaled density

(21) |

and resolve Eq. (14) into -specific contributions using
the exact resultLevy and Perdew (1985); Burke *et al.* (1997); Ernzerhof *et al.* (1997)

(22) |

Since there is no dependence for exchange, we can recast Eq. (22) using the correlation-energy density:

(23) |

The scaling results for and can be directly applied to individual components of the XC functional (as they are linear in the functional expression).

The scaling results, Eqs. (22) and (23), reflect properties
of the approximations that are implicitly made in crafting the PBE and vdW-DF-cx functionals. The existence of a well-understood coupling-constant scaling has been used to rationalize the formulation of the PBE0 hybridAdamo and Barone (1999) based on PBE.Burke *et al.* (1997); Ernzerhof *et al.* (1997)
Noting that a similar rationale exists for the coupling constant scaling of vdW-DF-cx, some of us have recently motivated the introduction of correspondingly defined vdW-DF hybrids,
including vdW-DF-cx0, which replace the
vdW-DF-cx exchange component with a fraction
of Fock exchange.Berland *et al.* (2017)

The scaling results, Eqs. (22) and (23), follow from an analysis of the many-particle wavefunction ground-state solution corresponding to a specific density and a specific strength of the electron-electron interaction. The detailed arguments are given elsewhere; For completeness, we include a renormalization-type argument for this observation in the appendix. Here we simply note that the wavefunctions solving the Hamiltonian themselves scale according to

(24) |

and that this formal equivalence is sufficient to establish the
scaling.Levy and Perdew (1985); Burke *et al.* (1997)

### ii.2 Access to the kinetic-correlation energy

Below we drop explicit references to the density functional nature when working with spatially-resolved energy contributions such as and , except when specifically needed for the discussion.

In the KS scheme Kohn and Sham (1965), we nominally focus on computing the so-called KS kinetic energy^{2}

(25) | |||||

(26) |

for occupied orbitals . As in the Quantum-Espresso
package,Giannozzi *et al.* (2009) we compute the KS kinetic energy as a spatial integration

(27) |

over positive definite contributions

(28) |

defined by the set of occupied orbitals. This representation of the KS kinetic energy is simply related to the summation over single-particle contributions, Eq. (26). The descriptions differ only in the inclusion of an Poisson-type term

(29) |

and give the same total KS kinetic energy, Eq. (27), upon spatial integration.

We typically compute DFT energies of combined systems ‘AB’ and of the relevant fragments, ‘A’ or ‘B’, to understand binding (with suitable adjustments in the case of related problems such as material cohesion). The mean-field electrostatic energy among electrons, that is, the Hartree term

(30) |

is one important contribution as it approximates . For analysis, we track binding-induced changes like

(31) | |||||

(32) | |||||

(33) |

We also track corresponding expressions for binding-induced
changes in, for example, the total nonlocal-correlation term .
In our discussion, we call such differences binding contributions.^{3}

Computational results for the binding-induced changes in the KS kinetic energy
and in the mean-field electrostatic energy
often suffice for a characterization of covalent bonds in
molecules and materials.von Hopffgarten and Frenking (2012) This is because
the combination allows us to characterize and understand orbital
hybridization.von Hopffgarten and Frenking (2012); Schmidt *et al.* (2014); Rohrer and Hyldgaard (2011)
However, for noncovalent bonds we have to look further than changes
in . One can generally sort chemical bonds from knowledge of
the average orbital energy.Rahm and Hoffmann (2016) The average orbital energy
is a measure that will, in principle, reflect all correlation effects,
including those that are manifested in the kinetic energy.

A formal analysis of the DFT variational schemeLevy and Perdew (1985) shows that

(34) |

where , Refs. Levy and Perdew, 1985; Burke *et al.*, 1997.
Using the density-scaling analysis, it immediately follows that

(35) |

A similar equation connects and . For any given system (solution density ), we use numerical differentiation to compute and from Eq. (22), and and from Eq. (34).

The electron-electron interaction effects in the physical systems are
now formally available for computation (in the approximations that define ).
In particular, we can study the electrodynamical coupling among plasmonsMahan (1965); Maggs and Ashcroft (1987); Rapcewicz and Ashcroft (1991); Hyldgaard *et al.* (2014)
at since is available via
Eq. (35). This value is the nonlocal-correlation
part of which, by definition, is given by a contour integral of
the response evaluated at full electron-electron interaction strength,
Eqs. (11) and (13).

We note in passing that is also the nonlocal-correlation part of the electron-electron interaction expectation value

(36) |

Since the XC energy functional is defined we can use Eq. (35) for the formal identification

(37) |

The formal equivalence Eq. (36) follows by subtracting the LDA and gradient-corrected exchange components.

Figure 1 shows (computed results for) the coupling constant scaling
for the XC contribution (solid red curve) to the total energy of the N molecules.
The specific scaling results are here provided for vdW-DF-cx (using the formal
derivation of the scaling for detailed in the following
subsection): However, the behavior is generic and thus similar to what has previously been reported and discussed for PBE.Burke *et al.* (1997); Ernzerhof *et al.* (1997)

We note that the exchange and correlation components, and , used for DFT calculations in the KS scheme, are integrals of the
indicated variations. The
exchange value traces a horizontal line (dotted curve separating red and green areas) in
Fig. 1. In contrast, the correlation begins at zero but changes to a significant magnitude at . It is straightforward to verifyBurke *et al.* (1997) that the area of the green region is minus the functional
approximation for . Importantly, we can immediately extract the corresponding kinetic-correlation energy using Eq. (35), that is, as the area of the blue region below the variation but above the value of the
limit.

### ii.3 Coupling-constant scaling and kinetic-correlation energy in vdW-DF-cx

To compute the kinetic-energy component of vdW binding, we need only consider the density
scaling for the correlation parts, namely . Moreover,
the coupling constant scaling for the LDA part, , has previously been discussed, as it is part of the GGA characterization.Levy and Perdew (1985); Levy *et al.* (1985); Levy (1991); Görling and Levy (1993); Levy (1995); Levy *et al.* (1996); Perdew *et al.* (1996c); Ernzerhof *et al.* (1997); Burke *et al.* (1997)

To explore the coupling-constant scaling of we first summarize the vdW-DF formulation of this nonlocal-correlation energy. Any semilocal XC density functional can be characterized by a local energy-per-particle density

(38) |

where is a function of just the
local density and the scaled density gradient .
We further split into
exchange and correlation components, and
. The local variation in the
inverse length scale for the plasmon-pole description can
then be expressed Dion *et al.* (2004); Thonhauser *et al.* (2007)

(39) | |||||

(40) |

Here denotes the local value of the Fermi wave-vector
and is the energy-per-particle density in LDA exchange.
The nonlocal-correlation energy is computedDion *et al.* (2004)

(41) |

using a universal-kernel formulation,Dion *et al.* (2004, 2005); Thonhauser *et al.* (2007); Román-Pérez and Soler (2009)

(42) | |||||

(43) | |||||

(44) |

The universal kernel is tabulated and permits
an efficient numerical evaluation through fast-Fourier transforms.Román-Pérez and Soler (2009)
The generalization to scaling in spin-polarized cases is also completely
specified, since it amounts to a simple rescaling of the inverse length scale ,
Ref. Thonhauser *et al.*, 2015.

For scaling analysis (and coding), it is convenient to introduce as a short hand for the coordinate scaling and to represent the density variation in terms of . The local values of the scaled density gradient are . The density scaling leaves unchanged and the effect amounts to computing the changes in Fermi vector and in the LDA correlation components. This is done in terms of the corresponding scaling .

Using to denote the exchange-enhancement factor of , the overall scaling of the inverse length scale can be expressed , where

(45) |

Introducing also

(46) | |||||

(47) |

we can compute

(48) |

We complete the scaling analysis of via Eq. (22) and of via Eq. (23). Specifically, we express the scaling of the nonlocal-correlation energy density

(49) |

For the numerical evaluation we adopted the scheme proposed by Román-Pérez and Soler Román-Pérez and Soler (2009) (as implemented in quantum-espresso) to calculate ; the
calculation is similar to the calculation of in Ref. Callsen *et al.*, 2012. We note that
provides a spatial
mapping of all nonlocal-correlation effects that exist in
the fully interacting system, as a direct implication of
Eq. (35).

Moreover, as part of this characterization, we can now compute the spatially resolved kinetic-correlation energy and nonlocal-kinetic-correlation energy . First, we simply add the knownLevy and Perdew (1985); Levy (1991); Görling and Levy (1993); Levy (1995) coupling constant scaling of the LDA correlation-energy density, , entering in Eq. (18):

(50) |

Next, we adapt Eq. (34) to descriptions of energy densities

(51) | |||||

(52) |

and evaluate the derivatives numerically.

Figure 2 documents the coupling-constant scaling of and explains an approximately linear variation. The top panel shows the coupling-constant scaling for the ratio for typical contributions to the binding of an N molecule. Specifically, starting from the known and values, the panel traces the variation in the XC components of for conditions that roughly correspond to the binding region of N (red curves) and to electron density tails of atoms and molecules (blue curves).

The middle panel of Fig. 2 shows the coupling-constant scaling of (solid curve) for the density of the N molecule. We note that the scaling in is exactly offset by the scaling of coordinates in Eqs. (46) and (47). Thus if we assume that the scaling of is set by the scaling in , there would be no -dependence in the kernel arguments, and . In this type of approximations there is then no scaling in the corresponding approximations for .

The middle panel furthermore shows two potentially relevant such approximations motivated by the analysis of the typical variations in the ratio. The first assumes that we can ignore the influence of the correlation part completely (giving the green dashed-dotted line); the second assumes that the ratio can at any given point be taken as fixed at the value (giving the blue dashed curve). The second choice effectively amounts to simply setting . Neither of them is a good description for Eq. (48). On the other hand, we add a weight on when computing , Eq. (22). We label the second approximation as ‘linear’ since it leads to and this is sometimes an acceptable approximation.

The bottom panel of Fig. 2 contrasts the resulting scaling of the ACF integrand for the nonlocal correlation contribution against that of the LDA correlation . The panel also shows (green and blue dashed curves) the scaling that results by inserting either of the approximations discussed in the middle panel into Eq. (22). Interestingly, the scaling relevant for the N-molecule total energy is found well approximated by using . The nonlocal-correlation part of the kinetic energy is just minus the total nonlocal-correlation binding contribution in such special cases.

## Iii Computational details

We focus our discussion and mapping of many-body physics effects in results obtained
using the vdW-DF-cx version.Berland and Hyldgaard (2014) In vdW-DF-cx, the total exchange component in Eq. (2) is picked so that binding contributions from can generally be ignored.Berland and Hyldgaard (2014); Hyldgaard *et al.* (2014)

The vdW-DF-cx version performs well, on par with or better than
PBE, for characterizations of many bulk, surface, and interface properties.Berland *et al.* (2014, 2015); Rangel *et al.* (2016); Brown-Altvater *et al.* (2016); Ambrosetti and Silvestrelli (2016); Gharaee *et al.* (2017); Petuya and Arnau (2017); Loncaric *et al.* (2017); Olsson *et al.* (2017); Wang *et al.* (2017a, b) The vdW-DF-cx version
has proven itself useful also in the description of binding and
function of layered materials, at surfaces, and of
molecules.Björkman (2014); Erhart *et al.* (2015); Sadeghi *et al.* (2015); Clay *et al.* (2016); Löfgren *et al.* (2016a); Fritz *et al.* (2016); Javaid and Akhtar (2016); Borca *et al.* (2015); Hellström *et al.* (2016); Mehboudi *et al.* (2016a, b); Löfgren *et al.* (2016b); Kuisma *et al.* (2016); Cechal *et al.* (2016); Zhong *et al.* (2017); Borca *et al.* (2017); Brihuega and Yndurain (2017); Kebede *et al.* (2017); Borck and Schröder (2017); Berland *et al.* (2017)

Our calculations are based on the plane-wave Quantum Espresso
package,Giannozzi *et al.* (2009) which already has the consistent exchange vdW-DF-cx
version Berland and Hyldgaard (2014) as well as the rigorous spin extension
of the vdW-DF method.Thonhauser *et al.* (2007) Core electrons are
represented by Troullier-Martins type normal-conserving
pseudo potentials using a 80 Ry wavefunction cutoff.

This paper also introduces a post-processing ACF-analysis
code, termed ppACF, which tracks the
system-specific coupling constant variation in
(for standard GGA and vdW-DF versions).
The code adapts the post-processing components of the
quantum espresso package,Giannozzi *et al.* (2009) into which ppACF
will also be released.

The ppACF code takes as input the set of quantum-espresso solution files (available after completion of the DFT calculations). It outputs the coupling-constant scaling analysis and the spatial variation in the kinetic-correlation energy density. For convenience, it also outputs the spatial variation in the set of XC components.

Our numerical analysis is based on comparing binding-energy contributions for the various components of the total DFT description. To discuss the binding ‘AB’ of fragments ‘A’ and ‘B’, the ppACF code outputs the spatial variation in all XC components. The code furthermore uses the coupling-constant scaling in the spatially resolved correlation terms, and , to numerically determine and output (for any given fragment and for the combined system) the spatial variation in and as well as in and . We then obtain, for the (spatially-resolved) binding-energy contributions

(53) | |||||

(54) | |||||

(55) | |||||

(56) |

from simple numerical subtractions.

For completeness, the ppACF code outputs the spatial variation in KS kinetic energy and in a spatially resolved measure of the full kinetic energy

(57) |

Again by numerical subtractions we can then define spatially-resolved kinetic binding energy contributions

(58) | |||||

(59) |

A mapping of the total kinetic energy binding contributions Eq. (59) will, in principle, always change if we base the definition, Eq. (57), on instead of on , using Eq. (29). This is true even if the integral values and remain the same. Qualitative differences in the resulting total-kinetic-energy mappings are visible for covalent bonding, but not for the cases of noncovalent inter-molecular interactions that we have investigated.

The set of top panels of Fig. 3 compares the coupling constant variation in the contributions and to the H, N, and O atomization energies, as computed in vdW-DF-cx. The scaling and the total kinetic-correlation energy contributions vary significantly between these traditional molecular binding examples. The total kinetic-correlation energy contribution to binding is given by the light blue area under the scaling curve. The value of is dominated by the part that originates from the LDA correlation energy.

The set of bottom panels focus on the coupling-constant scaling of the nonlocal-correlation contribution to the molecular cohesion. The coupling-constant variation in can be either upward or downward concave because it is only a part of the kinetic-correlation energy. The upwards and downwards concave behavior corresponds to positive and negative values of binding contributions, respectively. The dark-blue areas indicate the magnitude of this binding contribution. The supplementary materials includes a broad listing and comparisons of molecular-binding contributions , and . The comparison also lists KS binding contributions , making it clear that the kinetic-correlation energy can only play a significant role in the case of inter-molecular binding.

In the case of binding in the H, N, and O molecules, Fig. 3, we observe that the nonlocal-correlation contribution to binding is offset by a contribution to the nonlocal part of the correlation-kinetic energy. As shown in the supplementary materials, the same is true for many intra-molecular bonds, Tables S.I and S.II, and for all investigated inter-molecular interaction cases, Table S.III.

The binding in the total correlation term, , will be offset by a negative kinetic-correlation energy contribution , as suggested by the virial theorem. However, this need not hold generally for the nonlocal part of the kinetic-correlation energy contribution for intra-molecular binding, as further documented in the supplementary materials, Tables S.I and S.II. On the other hand, the compensation can be expected when the nonlocal part of the correlation-kinetic energy is a significant component, such as in most inter-molecular interactions.

## Iv Kinetic-energy mappings of molecular binding

We analyze and discuss the nature of binding both in H, N, and O molecules (having traditional chemical bonds) and in non-covalently bonded systems (where, in contrast, there is no pronounced orbital hybridization).

### iv.1 Intra-molecular interactions

Figure 4 shows that the kinetic-correlation energy is important in characterizations of intra-molecular binding. The figure details the spatial variation in the kinetic-correlation binding energy contribution , in the nonlocal-correlation-kinetic energy contribution , and in the vdW-DF-cx nonlocal correlation energy binding contribution for the H, N, and O molecules.

We note that the magnitude of the variation in is about an order-of-magnitude smaller than the KS kinetic-energy binding contribution (not shown) for these covalently bonded systems. Nevertheless, there is clear structure in both and and a directed nature or signature of vdW interactions even in these strongly bonded dimer molecules.

Supplementary materials Tables S.I
and S.II, supported by Fig. S1, provides a broader analysis of such
intra-molecular bindings. This is done both for the set of molecules
for which there already exists a PBE-based coupling-constant
analysis,Burke *et al.* (1997); Ernzerhof *et al.* (1997) and for the G2-1 benchmark set
of molecular atomization energies. Most of these systems are
covalently bonded, meaning that orbital hybridization plays the decisive role.
However, there are also some G2-1 cases, for example, alkali dimers, where the
nonlocal-correlation energy and the nonlocal part of kinetic-correlation
energy are important. Here, in the main text we concentrate on
characterizing the binding in H, N, and O.

The O kinetic-correlation energy binding contribution deserves a special discussion. The first O panel of Fig. 4 shows the variation of in a plane that contains the binding axis in the dimer. This plot has areas of opposite signs and implies a compensation. However, the overall kinetic-correlation energy contribution is still negative, , because the negative regions, away from the axis, have greater weight as we perform the spatial integration. The total, negative kinetic-correlation energy binding contribution is given by the light blue area shown in the left column in Fig. 3.

The second column of Fig. 4 shows the spatial variation in the nonlocal-correlation part of the kinetic-energy binding contribution, , for the three molecules. The variation in (first column) is generally dominated by the LDA contribution but adjusted by the variation in . The integrated binding contribution from the nonlocal part of the kinetic-correlation energy can be both negative or positive (as exemplified by the H and N cases). In such positive- cases, the binding contribution from is negative, i.e., the nonlocal-correlation energy is actually causing a repulsion in these intra-molecular binding cases. The supplementary materials Tables S.I and S.II provide a broader overview of the variation in nonlocal-correlation energy and kinetic-correlation energy effects that we document for intra-molecular binding.

The third column of Fig. 4 shows a mapping of the nonlocal-correlation energy contribution to binding, , allowing a contrast with the variation documented for and . We find that and are here essentially negative prints of each others. It follows that the full () nonlocal-correlation energy contribution, given by remains qualitatively described by the variation in in these cases.

Supplementary materials Tables S.I and S.II show that the nonlocal-correlation energy binding contribution can take either sign for intra-molecular binding. These tables also show that will typically then have the opposite sign. We therefore generally expect the and contributions to mirror each other, as in Fig. 4. The implication is that provides a qualitatively correct mapping of the vdW interaction in most covalently bonded cases.

Finally, we note that there are exceptions to this general trend, i.e, cases where and are both negative. Figures S1 and S2 in the supplementary materials, provide additional analysis for one of these cases, namely P. In such cases it is in principle necessary to compute both the and the variation to obtain a complete mapping, , Fig. S2. However, even in these cases it is still so that and the are approximately mirrors of each others. That is, even in P, there is no change in the qualitative observations presented above.

### iv.2 Inter-molecular interactions

Figure 5 compares kinetic-energy binding contributions , , , , in non-covalently bonded systems, i.e., in cases where there is no pronounced orbital hybridization. The top row shows our vdW-DF-cx based results for a benzene dimer, a case which is expected to have an essentially pure vdW (or dispersion) interaction. The bottom row shows results for a water dimer, a case that is predominantly hydrogen bonded, while the middle row explores the mixed-binding benzene-water case.

The supplementary materials Table S.III, supported by Fig. S3, provides a broader characterization of binding in the S22 benchmark set of such weakly bonded molecular complexes. In all of theses cases, the nonlocal-correlation energy contribution is positive while is negative. As such, the following discussion is generic.

For reference, the last column of Fig. 5 shows the spatial variation in the nonlocal-correlation energy binding, . This binding plays a decisive role in all cases. It is a core component of our vdW-DF-cx characterization and it is important for an accurate description of these molecular complexes.Berland and Hyldgaard (2014) In the case of dispersion-bonded systems, the contributions are the only sources of cohesion; In the case of the water dimer and the benzene-water complexes, there are also significant electrostatic components in the inter-molecular interactions.

The first, second, and third columns of Fig. 5 contrast the spatial variation in the binding contributions from the total kinetic energy, from the KS kinetic energy, and from the kinetic-correlation energy. There are no orbital hybridization effects in play, but it is important to note that the binding also pushes densities and thus orbitals around. Smaller signatures are therefore retained in the binding contribution from the KS kinetic energy, .

Contrasting the panels in the first and second column of Fig. 5,
we find that the KS kinetic-energy effects are still the major source of the
variation that we compute for .
Nevertheless, in the case of dispersion-bonded systems, we find
that there are also important contributions from the kinetic-correlation
energy . These contributions, shown in isolation in
the third column, arise primarily in the inter-molecular region, in areas that
have a sparseLangreth *et al.* (2009) (but not vanishing) electron density and
small-to-moderate density gradient. We sometimes refer to these binding
parts as a troughBerland and Hyldgaard (2014); Hyldgaard *et al.* (2014) but we are then emphasizing the presence
of important internal surfaces within such
sparse intermolecular regions.Rydberg *et al.* (2003); Kleis *et al.* (2008); Langreth *et al.* (2009)

The enhanced binding contributions from internal surfaces reflect the many-electron
nature of the vdW problem. The amplitudes of collective (plasmon) excitations are
themselves enhanced in the sparse surface-like intermolecular region and we should
then expect larger contributions to the systematic tracking of the electrodynamical coupling among plasmons.Mahan (1965); Langreth and Perdew (1977); Maggs and Ashcroft (1987); Rapcewicz and Ashcroft (1991); Langreth and Vosko (1987); Rydberg *et al.* (2003); Dion *et al.* (2004); Hyldgaard *et al.* (2014); Berland *et al.* (2015)
The vdW enhancement can also be interpreted as reflecting image-plane effects at (internal
or external) surfaces,Zaremba and Kohn (1976, 1977); Nordlander and Harris (1984); Rydberg *et al.* (2003); See, for example, M. Persson and S. Andersson,
Chapter 4 (2008); Lee *et al.* (2011, 2012b); Tao and Rappe (2014) or as multipole response-effects effects when arising outside
molecules.Becke and Johnson (2005, 2007); Kleis *et al.* (2008); Berland and Hyldgaard (2010); Tao and Perdew (2014); Tao *et al.* (2015, 2017)
In any case, the vdW-DF-cx handling of screeningMahan (1965); Rapcewicz and Ashcroft (1991); Dion *et al.* (2004); Berland *et al.* (2014); Hyldgaard *et al.* (2014)
provides mechanisms to track the expected vdW enhancement in the sparse intermolecular regions,
at important internal surfaces.Rydberg *et al.* (2003); Lazić *et al.* (2012); Berland and Hyldgaard (2013); Hyldgaard *et al.* (2014)

The fourth column of Fig. 5 shows the nonlocal part of the kinetic-correlation energy, . Contrasting the third and fourth columns makes it clear that the LDA component of generally masks the variation in