Inclusion of Van der Waals Interactions in DFT using Wannier Functions without empirical parameters

We describe a method for including van der Waals (vdW) interactions in Density Functional Theory (DFT) using the Maximally-Localized Wannier functions (MLWFs), which is free from empirical parameters. With respect to the previous DFT/vdW-WF2 version, in the present DFT/vdW-WF2-x approach, the empirical, short-range, damping function is replaced by an estimate of the Pauli exchange repulsion, also obtained by the MLWFs properties. Applications to systems contained in the popular S22 molecular database and to the case of adsorption of Ar on graphite, and Xe and water on graphene, indicate that the new method, besides being more physically founded, also leads to a systematic improvement in the description of systems where vdW interactions play a significant role.


Introduction
Density Functional Theory (DFT) is a well-established computational approach to study the structural and electronic properties of condensed matter systems from first principles. Although current, approximated density functionals allow a quantitative description at much lower computational cost than other first principles methods, they fail [1] to properly describe dispersion interactions. These forces originate from correlated charge oscillations in separate fragments of matter and the leading component (namely, the lowest-order perturbative contribution) is represented by the R −6 van der Waals (vdW) interaction [2], due to correlated instantaneous dipole fluctuations. vdW interactions play a fundamental role in determining the structure, stability, and function of a wide variety of systems, including molecules, clusters, proteins, nanostructured materials, molecular solids and liquids, and in adsorption processes of fragments weakly interacting with a substrate ("physisorbed").
In all these methods a certain degree of empiricism is present since the energetic vdW-correction term is multiplied by a short-range damping function, which is introduced not only to avoid the unphysical divergence of the vdW correction at small fragment separations, but also to eliminate double countings of correlation effects (in fact standard DFT approaches are able to describe short-range correlations). This damping function contains one or more empirical parameters which are typically set by a trial and error approach or/and are fitted using some reference database.
In a recent paper [24] we have presented a new method, named DFT/vdW-WF2-x, to overcome the above limitation, where the empirical, short-range, damping function is replaced by an estimate of the Pauli exchange repulsion, also obtained by the MLWFs properties. Here we use the DFT/vdW-WF2-x scheme, already successfully applied [24] to the popular S22 benchmark set [25] of weakly interacting molecules and to the case of an Ar atom interacting with graphite, also to describe the adsorption processes of a single Xe atom and a single water molecule interacting with graphene. Rare gases adsorbed on graphene are interesting since they represent prototypical adsorption systems. Moreover the study of the interactions of graphene with water is very important as a model for the characterization of the interface between water and hydrophobic substrates [26]; in addition, a single water molecule adsorbed on graphene represents a weakly interacting system involving a complex mixture of hydro-gen bonding, electrostatic, and van der Waals interaction, thus providing a significant test model for any electronic structure theory. Our results are compared with reference data and confirm that the DFT/vdW-WF2-x method leads to a systematic improvement in the description of systems where vdW interactions play a significant role.

Method
Here we briefly describe the DFT/vdW-WF2-x method, introduced in ref. [24], where further details can be found. Basically, the electronic charge partitioning is achieved using the Maximally-Localized Wannier Functions (ML-WFs), which are obtained from a unitary transformation in the space of the occupied Bloch states, by minimizing the total spread functional [10]: If two interacting atoms, A and B, are approximated[2] by coupled harmonic oscillators, the vdW energy correction can be taken to be the change of the zero-point energy of the coupled oscillations as the atoms approach; if only a single excitation frequency is associated to each atom, ω A , ω B , then where Z A,B is the total charge of A and B, and R AB is the distance between the two atoms (e and m are the electronic charge and mass). Now, adopting a simple classical theory of the atomic polarizability, the polarizability of an electronic shell of charge eZ i and mass mZ i , tied to a heavy undeformable ion can be written as Then, given the direct relation between polarizability and atomic volume [27], we assume that α i = γS 3 i , where γ is a proportionality constant, so that the atomic volume is expressed in terms of the MLWF spread, S i . Rewriting eq. (2) in terms of the quantities defined above, one obtains an explicit, simple expression for the C 6 vdW coefficient: The constant γ can then be set up by imposing that the exact value for the H atom polarizability (α H =4.5 a.u.) is obtained.
In order to achieve a better accuracy, one must properly deal with intrafragment MLWF overlap: in fact, the method is strictly valid for nonoverlapping fragments only; now, while the overlap between the MLWFs relative to separated fragments is usually negligible for all the fragment separation distances of interest, the same is not true for the MLWFs belonging to the same fragment, which are often characterized by a significant overlap. This overlap affects the effective orbital volume, the polarizability, and the excitation frequency (see eq. (3)), thus leading to a quantitative effect on the value of the C 6 coefficient. We take into account the effective change in volume due to intrafragment MLWF overlap by introducing a suitable reduction factor ξ obtained by interpolating between the limiting cases of fully overlapping and non-overlapping MLWFs [24]. We therefore arrive at the following expression for the C 6 coefficient: where ξ A,B represents the ratio between the effective and the free volume associated to the A-th and B-th MLWF [24].
Finally, in the previous DFT/vdW-WF2 method, the vdW interaction energy was computed as: where f (R i j ) is a short-range damping function, which is introduced not only to avoid the unphysical divergence of the vdW correction at small fragment separations, but also to eliminate double countings of correlation effects (in fact standard DFT approaches are able to describe short-range correlations); it is defined as: The parameter R s represents the sum of the vdW radii R s = R vdW i + R vdW j , with (by adopting the same criterion chosen above for the γ parameter) where R vdW H is the literature [28] (1.20 Å) vdW radius of the H atom, and, following Grimme et al. [29], a 20. Although a is the only ad-hoc parameter of the method, while all the others are only determined by the basic information given by the MLWFs (that is from first principles calculations) and in many applications the results are only mildly dependent on the particular value of a, nonetheless, this parameter, together with the choice of a specific form of the above damping function, clearly imply a certain degree of empiricism.
In order to overcome this limitation, in the DFT/vdW-WF2-x approach we replace the somehow artificial, shortrange damping function by a term that directly measures the quantum mechanical Pauli exchange repulsion between electronic orbitals and can be entirely expressed in terms of the MLWFs properties, without the need of introducing empirical parameters. Following ref. [30], using the dipole approximation for the Coulomb interaction, the exchange integral, for two closed electronic shells with total zero spin, is simply given by: where q indicates the electronic charge of each electronic shell and O is the overlap integral between the electronic shells, separated by R. In our specific case, assuming that an electronic orbital is described by the wave function relative to a quantum harmonic oscillator: where S A is the spread of the corresponding MLWF, then one can easily obtain that : Then, the exchange integral can be expressed in terms of the MLWFs spreads as : (12) Therefore, in the DFT/vdW-WF2-x method, the vdW interaction energy is computed as: In this way the vdW energy correction is evaluated as the sum of two terms, both expressed in terms of the MLWfs spreads, thus making explicit the direct connection between attractive and repulsive parts of the vdW interaction [30].
Of course there are very weakly bonded systems, entirely dominated by vdW effects, where the repulsive term is not relevant for determining the equilibrium complex configuration. For instance, in the Ar-dimer case, DFT/vdW-WF2 and DFT/vdW-WF2-x predict the same equilibrium Ar-Ar distance and the same binding energy (within 0.1 meV). However, in most cases, a proper treatment of short-range repulsion is crucial to correctly describe the minimum, equilibrium configuration.
Applications of DFT/vdW-WF2-x to the case of Xe and water adsorbed on graphene have been performed with the Quantum-ESPRESSO ab initio package [31] and the MLWFs have been generated as a post-processing calculation using the WanT program [32], with ultrasoft pseudopotentials to describe the electron-ion interactions and taking PBE [33] as the reference, Generalized Gradient Approximation (GGA) DFT functional. PBE is chosen because it probably represents the most popular GGA functionals for standard DFT simulations of condensed-matter systems.
As expected, considering the whole S22 database, pure PBE performs poorly and a substantial improvement can be obtained by vdW-corrected approaches. More importantly, the performances of the new DFT/vdW-WF2-x scheme are clearly better than those of the previous DFT/vdW-WF2 method. In particular, the general tendency of DFT/vdW-WF2 to overbind is considerably reduced by DFT/vdW-WF2-x. Interestingly, with DFT/vdW-WF2-x the mean absolute error (MAE), 0.78 kcal/mol, is well below the so-called "chemical accuracy" threshold of 1 kcal/mol, required to attribute a genuine quantitative character to the predictions of an ab initio scheme. Moreover, DFT/vdW-WF2-x performs better than the more sophisticated vdW-DF and vdW-DF2 methods, based on the use of a nonlocal expression for the correlation energy-term, is comparable, as far as the S22 database is concerned to PBE-D3, and its performances are not too inferior to those of the optB88-vdW, optB86b-vdW, vdW-DF-C09, vdW-DF-cx, rev-vdW-DF2, rVV10, VV10, PBE+TS-vdW, and PBE+MBD schemes, which are among the most accurate, presently available, vdW-corrected DFT approaches for noncovalently bound complexes [38-42, 44, 46].
It is certainly important to test the applicability of the DFT/vdW-WF2-x method also to extended systems, which of course are of particular interest because, in this case, higher-quality chemistry methods are typically too computationally demanding. We have therefore considered [24] the adsorption of a single Ar atom on graphite, that is a typical physisorption process, and also the cases of a Xe atom and a water molecule interacting with graphene, which are adsorption processes characterized by a little stronger interactions with the substrate. Calculations have been performed using the same DFT approach followed in ref. [20] and considering the adsorption on the most favored hollow sites only, illustrated in Figure 1 for the case of water on graphene (with the two OH bonds pointing down to graphene).  Tables 2 and 3 report the binding energy, E b , and equilibrium distance, R, for Ar-graphite, Xe-graphene, and water-graphene. Data obtained at the PBE, DFT/vdW-WF2, and DFT/vdW-WF2-x level are compared to reference theoretical and experimental estimates (the latter relative to adsorptions on graphite, in the absence of specific experimental data on graphene).
As can be seen, in all the cases the pure PBE functional largely underbinds while vdW-corrected schemes predict much larger binding energies and shorter equilibrium distances. Moreover, DFT/vdW-WF2-x represents an evident improvement compared to DFT/vdW-WF2, in line with what previously observed in the application to the S22 database, showing again that, at equilibrium distances, the effect of the repulsive term is significant. In particular, for Ar-graphite the DFT/vdW-WF2-x binding energy (-115 meV) is very close to the values reported by Tkatchenko et al. [48] (-116 meV) and Bichoutskaia and Pyper [49] (-111 meV). Moreover, this value is also compatible with one of the few experimental estimates, represented by the Table 2. Binding energy, E b for an Ar atom adsorbed on graphite, and Xe and water on graphene, compared to available theoretical (for Ar-graphite refs. [48,49], for Xe-graphene refs. [50][51][52], for Xe-graphite refs. [53,54], for water-graphene refs. [55,56]) and experimental (for Ar-graphite refs. [50,57], for Xe-graphite refs. [50,54], for water-graphite refs. [58]) reference data.   [59]. It is also close to the "best estimate" (obtained from a combination of experimental and theoretical, mainly semiempirically-based, data) reported in the milestone review paper by Vidali et al. [50], that is -99 ± 4 meV. Our Ar-graphite equilibrium distance computed by DFT/vdW-WF2-x is instead somehow larger than that reported in other theoretical studies [48,49] (3.3 Å), and also than the "best estimate" by Vidali et al. [50] of 3.1 ± 0.1 Å (and an old experimental measurement of 3.2 ± 0.1 Å, see ref. [57]). However one sould observe that an accurate estimate of this quantity is more difficult due to the relatively shallow potential energy curve of this system. Concerning the interaction of Xe with graphene, previous theoretical calculations exist (see Tables 2 and 3); in particular, the MP2 method [51] gives an estimate of -143 meV for the binding energy at a separation of the Xe atom from the substrate of 3.6 Å. Our DFT/vdW-WF2-x result of -170 meV, is not far from the MP2 value also considering that the MP2 results were not directly obtained from a calculation on an infinite graphene sheet, as in the present approach, but were instead derived by an extrapolation procedure using calculations on Xe interacting with polycylic aromatic hydrocarbons such as coronene. Moreover, the DFT/vdW-WF2-x binding energy agrees well with the "best estimate" (-162 ±4 meV) reported by Vidali et al. [50] for Xe on graphite. Instead the binding energy obtained by the previous DFT/vdW-WF2 scheme turns out to be evidently overestimated, and the same conclusion applies by referring to the experimental reference range too (see Table 2). With respect to the reference, theoretical and experimental estimates for the equilibrium distance, our DFT/vdW-WF2-x value of 4.0 Å appears to be significantly larger, but, for this quantity, the same considerations made above for Ar-graphite apply.
Finally, as far as the adsorption of a water molecule on graphene/graphite is concerned, the experimental estimates of the binding energy are scarce [26]: for watergraphite an experimental estimate of -156 meV exist, which is very close to our DFT/vdW-WF2-x value, while interaction energies in the range from -65 to -97 meV were obtained from empirical force-field simulations [60], aimed at reproducing the experimentally observed contact angles of water droplets on graphite. From the theoretical point of view, highly accurate ab initio data (coupled clusters at the complete basis set limit) are currently limited to the water-benzene system [61] and other quantum-chemistry estimates are quite scattered and based on extrapolations of the interaction energies calculated for polycyclic aromatic hydrocarbons of increasing size. A DFT/coupled cluster method (DFT/CC) has been applied to evaluate the electronic binding energy of a water molecule with graphite and graphene, giving as a result, -155 meV at 3.24 Å and -138 meV at 3.26 Å, respectively [26]. Also in this case DFT/vdW-WF2-x performs well: not only the binding energy but even the equilibrium distance are close to the reference values, although the improvement with respect to the previous DFT/vdW-WF2 method is smaller than for Ar-graphite and Xegraphene. This is probably because the water-graphene interaction is not entirely due to vdW forces but results from a complex interplay of hydrogen bonding, electrostatic, and vdW effects.

Conclusions
In summary, we have described the DFT/vdW-WF2-x method for including van der Waals (vdW) interactions in Density Functional Theory using the MLWFs, which is free from empirical parameters, and, with respect to the previous DFT/vdW-WF2 method, replaces the empirical, short-range, damping function by an estimate of the Pauli exchange repulsion, also obtained by the MLWFs properties. Applications to systems contained in the popular S22 molecular database and to the case of adsorption of Ar on graphite, and Xe and water on graphene, indicate that the new method, besides being more physically founded, represents a systematic improvement in the description of systems where vdW interactions play a significant role.