# 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}^{1}1In the case of
vdW-DF2Lee et al. (2010) by a formulation that reflects an exchange-scaling to the
high-density limit. this internal function comprises LDA with gradient corrections defined by analysis of screened exchange.Langreth and Vosko (1987); Dion 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}^{2}2In discussions of DFT, is sometimes called the single-particle kinetic energy. We prefer the Kohn-Sham label as is always a
single-particle operator. from single-particle expectation values

(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}^{3}3The wording ‘binding contribution’ is used to describe any component of the molecular binding even if,
for example, is negative. Similarly, we
use the word ‘spatially resolved binding contributions’ to describe binding-induced changes in the spatial variation of, for example,
XC energy terms, like ; Again this
term is used without regards to the sign of the integrated values.

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 . However, the signatures of the nonlocal kinetic-correlation part dominate in the spare-density regions for dispersion-bonded systems. The nonlocal part also remains a non-vanishing part of full kinetic-correlation energy in the hydrogen-bonded and mixed binding cases.

Figure 6 shows the computed binding contributions in the stacked uracil dimer and in the stacking of adenine and thymine. The top row shows the investigated geometries from the S22 benchmark set.Jurecka et al. (2006) Supplementary materials Fig. S.4 shows the variation in the kinetic-energy binding contributions for these systems.

The panels in the bottom two rows of Figure 6 contrast the variation in the and accounts of the nonlocal-correlation binding. The contributions are computed for the cuts indicated by the two dashed lines in the top panels. We find that including the nonlocal part of the kinetic-correlation energy enhances the binding signatures found in the variation in our mapping of the total nonlocal correlation binding, shown in the pair of lower panels. The same is true for the wider set of S22 cases. Comparing the fourth and fifth column of Fig. 5 we see that they are essentially negative prints of each other, i.e., naturally leading to an enhancement of signatures in .

For both inter- and intra-molecular interactions the variation, can effectively be mapped using either or .

We also find that the signatures are, in effect, channelled into pockets of binding, Figs. 5 and 6. The dominant contributions are located in the intermolecular (trough) regions but concentrated in areas that resemble orbitals. In other words, we can, in principle, use a similar form of bond-type characterization for a qualitative discussion of the nonlocal-correlation binding, drawing on an analogy with discussions of chemical bonds.

Before using this vdW-bond mapping analysis, below, we emphasize that this binding is much weaker and that the concentration of binding (in such inter-molecular pathways) arises for a different physical reasonThonhauser et al. (2007) than in chemical bonds. Chemical binding can exist in the limit (in a Hartree-Fock description) by orbital hybridization, but that cannot happen for the nonlocal-correlation (or vdW) binding. While the inclusion of the energy term in DFT calculations causes density changes, and therefore an electrostatic signature,Thonhauser et al. (2007) the total nonlocal-correlation binding reflects, instead, an energy gainThonhauser et al. (2007) produced by collective exitations, i.e., by plasmons described by the screening properties.Mahan (1965); Maggs and Ashcroft (1987); Andersson et al. (1996); Rydberg et al. (2003); Dion et al. (2004); Langreth et al. (2005); Hyldgaard et al. (2014)

(Å) | (meV) | (meV) | |||
---|---|---|---|---|---|

Ne dimer | 3.09 | (3.09) | 10.2 | (3.64) | 14.0 |

Ne trimer | 3.09 | 30.2 | 42.8 | ||

Ar dimer | 3.99 | (3.76) | 18.9 | (12.3) | 25.7 |

Ar trimer | 4.02 | 55.1 | 75.4 | ||

Kr dimer | 4.33 | (4.01) | 22.1 | (17.3) | 30.7 |

Kr trimer | 4.35 | 64.6 | 90.8 |