{'fraction_non_alphanumeric': 0.03731293201798537, 'fraction_numerical': 0.011878397422991745, 'mean_word_length': 4.789432789432789, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Van der Waals Interactions in DFT using Wannier Functions without empirical parameters

A new implementation is proposed 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 method, 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 an Ar atom interacting with graphite, and comparison with reference data, indicate that the new method, besides being more physically founded, also leads to a systematic improvement in the description of vdW-bonded systems.


I. 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. Dispersion forces originate from correlated charge oscillations in separate fragments of matter and the most important component is represented by the R −6 van der Waals (vdW) interaction, 2 due to correlated instantaneous dipole fluctuations. These interactions play a fundamental role in determining the structure, stability, and function of a wide variety of systems, including molecules, clusters, proteins, nanostructered materials, molecular solids and liquids, and in adsorption processes of fragments weakly interacting with a substrate ("physisorbed").
In the last few years a variety of practical methods have been proposed to make DFT calculations able to accurately describe vdW effects (for a review, see, for instance, refs. [3][4][5][6][7][8]. In this respect, a family of such methods, all based on the Maximally Localized Wannier Functions (MLWFs), 9 has been developed, namely the original DFT/vdW-WF scheme, [10][11][12] DFT/vdW-WF2 13 (based on the London expression and taking into account the intrafragment overlap of the MLWFs), DFT/vdW-WF2s 14 (including metal-screening corrections), and DFT/vdW-QHO-WF 15 (adopting the coupled Quantum Harmonic Oscillator model), successfully applied to a variety of systems: [10][11][12][13][14][15][16][17][18][19][20][21][22] small molecules, water clusters, graphite and graphene, water layers interacting with graphite, interfacial water on semiconducting substrates, hydrogenated carbon nanotubes, molecular solids, the interaction of rare gases and small molecules with metal surfaces,... In all these methods a certain degree of empiricism is present since the energetic vdWcorrection 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 the present paper we overcome the above limitation by presenting a new method, here-after referred to as DFT/vdW-WF2-x, where the empirical, short-range, damping function is replaced by an estimate of the Pauli exchange repulsion, also obtained by the MLWFs properties. The new approach is successfully applied to the popular S22 benchmark set 23 of weakly interacting molecules and also to the case of an Ar atom interacting with graphite.
The results are compared with reference data and indicate that the new method leads to a systematic improvement in the description of vdW-bonded systems.

II. METHOD
Here we describe and apply a new implementation of the DFT/vdW-WF2 method, introduced in ref. 13 and summarized below. Basically, the electronic charge partitioning is achieved using the Maximally-Localized Wannier Functions (MLWFs), which are obtained from a unitary transformation in the space of the occupied Bloch states, by minimizing the total spread functional: 9 The localization properties of the MLWFs are of particular interest for the implementation of an efficient vdW correction scheme: in fact, the MLWFs represent a suitable basis set to evaluate orbital-orbital vdW interaction terms. In particular, 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 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. This appears to be a physically sound choice since, in the H case, one knows the exact analytical spread, S i = S H = √ 3 a.u.
In order to achieve a better accuracy, one must properly deal with intrafragment MLWF overlap (we refer here to charge overlap, not to be confused with wave functions 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. In particular, since in the DFT/vdW-WF2 method the i-th MLWF is approximated with a homogeneous charged sphere of radius S i , then the overlap among neighboring MLWFs can be evaluated as the geometrical overlap among neighboring spheres. To derive the correct volume reweighting factor for dealing with overlap effects, we first consider the limiting case of two pairs (one for each fragment) of completely overlapping MLWFs, which would be, for instance, applicable to two interacting He atoms if each MLWF just describes the density distribution of a single electron; then we can evaluate a single C 6 coefficient using eq. (4) with Z A,B = 2, so that: Alternatively, the same expression can be obtained by considering the sum of 4 identical pairwise contributions (with Z = 1), by introducing a modification of the effective volume in such a way to take the overlap into account and make the global interfragment C 6 coefficient equivalent to that in eq. (5). This is clearly accomplished by replacing S 3 i in eq. (4) with ξS 3 i , where ξ = 1/2. This procedure can be easily generalized to multiple overlaps, by weighting the overlapping volume with the factor n −1 , where n is the number of overlapping MLWFs. Finally, by extending the approach to partial overlaps, we define the free volume of a set of MLWFs belonging to a given fragment (in practice three-dimensional integrals are evaluated by numerical sums introducing a suitable mesh in real space) as: where w f ree (r l ) is equal to 1 if |r l − r i | < S i for at least one of the fragment MLWFs, and is 0 otherwise.
The corresponding effective volume is instead given by where the new weighting function is defined as that is equal to the number of MLWFs contemporarily satisfying the relation |r l − r i | < S i .
Therefore, the non overlapping portions of the spheres (in practice the corresponding mesh points) will be associated to a weight factor 1, those belonging to two spheres to a 1/2 factor, and, in general, those belonging to n spheres to a 1/n factor. The average ratio between the effective volume and the free volume (V ef f /V f ree ) is then assigned to the factor ξ, appearing in eq. (8). Although in principle the correction factor ξ must be evaluated for each MLWF and the calculations must be repeated at different fragment-fragment separations, our tests show that, in practice, if the fragments are rather homogeneous all the ξ factors are very similar, and if the spreads of the MLWFs do not change significantly in the range of the interfragment distances of interest, the ξ's remain essentially constant; clearly, exploiting this behavior leads to a significant reduction in the computational cost of accounting for the intrafragment overlap. 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. The need for a proper treatment of overlap effects has been also pointed out by Andrinopoulos et al., 12 who however applied a correction only to very closely centred WFCs.
Finally, in the original DFT/vdW-WF2 method, the vdW interaction energy was computed as: where f (R ij ) 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 In order to overcome this limitation, we propose to improve the approach by replacing the somehow artificial, short-range 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. 27, 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 : Therefore, in this new DFT/vdW-WF2-x version of the method, the vdW interaction energy is computed as: In this way the vdW energy correction is evaluated as the sum of two terms, both ex-   39 and PBE+MBD. 40 For the sake of completeness we also report data relative to the semiempirical PBE-D3 26 approach and to the bare, non-vdW-corrected, PBE and BLYP functionals. In Table II is comparable, as far as the S22 database is concerned to PBE-D3, and its performances are only inferior to those of the rVV10, VV10, PBE+TS-vdW, and PBE+MBD schemes, which are among the most accurate vdW-corrected DFT approaches for noncovalently bound complexes. 38,40 As can be seen, looking at Table II, DFT/vdW-WF2-x(PBE) turns out to be better than DFT/vdW-WF2-x(BLYP) for both dispersion-dominated and mixed complexes, while instead the opposite is true for Hydrogen-bonded complexes.
In order to test the applicability of the present DFT/vdW-WF2-x method also to a representative of extended systems, which of course is of particular interest because, in this case, high-quality chemistry methods are typically too computationally demanding, we considered the adsorption of a single Ar atom on graphite, that is a typical physisorption process. Calculations have been performed using the same DFT approach followed in ref. 41 and considering the adsorption on the more favored hollow site only. Table III , is compared to that (see eq. (16)) of DFT/vdW-WF2- x(PBE), i<j J ij ex (R ij ). Clearly, the repulsive term vanishes for Ar-graphite distances larger than 4Å, however, around the equilibrium distance (and at shorter distances), it is much more substantial in the DFT/vdW-WF2-x(PBE) approach than in DFT/vdW-WF2(PBE).
The DFT/vdW-WF2-x(PBE) binding energy (-115 meV) is very close to the values reported by Tkatchenko et al. 42 (-116 meV) and Bichoutskaia and Pyper 43 (-111 meV). Moreover, this value is also compatible with one of the few experimental estimates, represented by the measurement of the latent heat of condensation relative to the adsorption of an Ar monolayer on graphite: -119 ± 2 meV/atom. 44 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., 45 that is -99 ± 4 meV. Our DFT/vdW-WF2-x(PBE) computed Ar-graphite equilibrium distance is instead somehow larger than that reported in other theoretical studies 42,43 (3.3Å), and also than the "best estimate" by Vidali et al. 45 of 3.1 ± 0.1Å (and an old experimental measurement of 3.2 ± 0.1Å, see ref. 45). 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.

IV. CONCLUSIONS
In summary, we have presented a new method for including van der Waals (vdW) interactions in Density Functional Theory using the MLWFs, which is free from empirical parameters. With respect to the previous DFT/vdW-WF2 method, in the present DFT/vdW-WF2x 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 an Ar atom interacting with graphite, and comparison with reference data, indicate that the new method, besides being more physically founded, also leads to a systematic improvement in the description of vdW-bonded systems.      Table III).