Improved proton-transfer barriers with van der Waals density functionals: Role of repulsive non-local correlation

Proton-transfer (PT) between organic complexes is a common and important biochemical process. Unfortunately, PT energy barriers are difficult to accurately predict using density functional theory (DFT); in particular, the generalized gradient approximation (GGA) tends to underestimate PT barriers. Moreover, PT typically occurs in environments where dispersion forces contribute to the cohesion of the system; thus, a suitable exchange-correlation functional should accurately describe both dispersion forces and PT barriers. This paper provides benchmark results for the PT barriers of several density functionals including several variants of the van der Waals density functional (vdWDF). The benchmark set comprises small organic molecules with inter- and intra-molecular PT. The results show that replacing GGA correlation with a fully non-local vdW-DF correlation increases the PT barriers, making it closer to the quantum chemical reference values. In contrast, including nonlocal correlations with the Vydrov-Voorhis (VV) method or dispersion-corrections at the DFT-D3 or the Tkatchenko-Scheffler (TS) level has barely any impact on the PT barriers. Hybrid functionals also increase and improve the energies and the best performance is provided by a hybrid version of the consistent-exchange van der Waals density functional vdW-DF-cx. For the formic acid dimer PT system, we analyzed the GGA exchange and non-local correlation contributions. The analysis shows that the repulsive part of the non-local correlation kernel plays a key role in the PT energy barriers predicted with vdW-DF.


INTRODUCTION
Proton transfer (PT) is an ubiquitous chemical reaction and many biochemical reactions involve PT. For instance, proton-coupled charge transfer is an established mechanism in enzymology 3 and collective proton transfer in DNA base pairs can give rise to rare tautomers which may lead to mutations. 4 In organic solid-state systems, PT can change the nature of the bonding and, consequently, the properties of a crystal. 5 PT is also one mechanism for electric polarization switching in the organic ferroelectrics. [6][7][8] High level quantum chemical methods, in particular, coupled-cluster with single and double and perturbative triple excitations (CCSD(T)) can provide accurate reference data for PT energy barriers, but its high computational cost makes it ill-suited for complex PT systems. The organic complexes in which PT occurs are often held together by van der Waals forces, making it important to assess the accuracy of predicted PT barriers with functionals that include dispersion forces.
Earlier benchmark studies have found that density functional theory (DFT) in the local density approximation (LDA) 9 and generalized-gradient approximation (GGA) 10,11 tend to underestimate PT barriers. 1,12 Hybrid functionals, which mix in a fraction of Fock exchange, [13][14][15][16] improve performance 1 which, in part, can be linked to reduced self-interaction error. 17 Patchkovskii et al. 18 found that Perdew-Zunger self-interaction corrected DFT improves reaction and activation energy barriers for 11 selected "difficult" reactions compared to * E-mail: kristian.berland@nmbu.no LDA and GGA. Although self interaction is one source of the PT barrier underestimation, lack of non-local correlation effects has also been suggested as a possible source of inaccuracy. 19 Several different methods for including dispersion forces in DFT have been developed, but most of these use "dispersion correction" scheme in which van der Waals forces are reintroduced with atomistic force-field correction. [20][21][22][23][24][25] One class of functionals [26][27][28][29][30][31][32][33][34][35] with a nonlocal account of the correlation is the Chalmers-Rutgers van der Waals (vdW-DF) density functional 26,34,36 which is a popular method to describe materials bonded by dispersion forces. vdW-DF is derived from exact criteria using a plasmon-model of the response properties of the electron-gas. 26,30,[35][36][37][38] In the theory, the non-local correlation energy takes the form, in which a kernel function φ(r 1 , r 2 ) connects two density regions n(r 1 ) and n(r 2 ). Unlike other dispersioncorrection methods, vdW-DF is designed so that E nl c vanishes seamlessly into a homogeneous electron gas limit without the inclusion of damping terms. Therefore, vdW-DF does not include gradient components of the GGA correlation in the total exchange-correlation en- . Over the years, several variants of vdW-DF [26][27][28][30][31][32][33]39 have been developed including hybrid variants. 40,41 vdW-DF has also inspired other non-local correlation functionals; in particular the Vydrov-Voorhis (VV10) 42 functional and its revision for planewave codes (rVV10). 43 While vdW-DF is foremost developed for describing dispersion-bonded systems, its explicit non-local corre- lation has also been found to improve various material properties 28,30,37,44 including image plane states on graphene. 45 This paper provides benchmark results for 22 different functionals including a number of vdW-DFs. The benchmark set labelled PT14 is based on the 9 charge-neutral intra-molecular PT systems by Karton et al. 2 and 5 interand intra-molecular PT systems by Mangiatordi et al. 1 computed with coupled-cluster method at the CCSD(T) level of theory. The set of systems are displayed in Fig. 1. In addition to vdW-DFs, we tested several GGAs, 46,47 two standard hybrid functionals, 15,16 and the strongly constrained and appropriately normed (SCAN) meta-GGA functional. 48 We also tested the rVV10 variant 42,43 as well as the SCAN-rVV10 functional. 49 Moreover, for several of the GGAs, we tested the effect of adding forcefield dispersion corrections at the D3 23 and Tkatchenko-Scheffler (TS) 20 level of theory.
We found, as detailed in Sec. III, that vdW-DF tends to increase PT barriers compared to GGA and therefore reduce the deviation with the reference data, in particular if the exchange functional is kept fixed. This trend was found to arise from a reduction in the negative correlation contribution to the PT barriers. To gain better understanding of this result, we performed an in-depth analysis of the case of formic acid dimer (system 2 in Fig. 1), as detailed in Sec. IV. The non-local correlation contribution was analyzed and found to be linked to repulsive short-range non-local correlation effects. This effect is similar to that of GGA-type correlation but with the local geometry-sensitivity inherited to vdW-DF.

II. METHODS
The benchmark calculations were carried out with the VASP software package, 50-52 except for the vdW-DF3-opt1, vdW-DF3-opt2, 33 and B3LYP(-D3) 16 calculations for which Quantum Espresso 53,54 was used, as these functionals are only implemented in the latter. In the supercells, 15Å vacuum padding was used to isolate the molecular systems employing the dipole correction scheme of Neugebauer et al. 55 The electronic self-consistency criteria was set to 10 −6 eV. The plane-wave energy cutoff was set to 1000 eV based on our convergence study as detailed in the next section. The VASP calculations used hard projected augmented waves (PAW) pseudopotentials (PPs) while the Quantum Espresso used the recently developed optimized norm-conserving Vanderbilt (ONCV). 56 For the exchange and correlation analysis, we obtained input data from the recently developed ppACF post-processing tool 57 which is distributed as part of the Quantum Espresso software package.

A. Pseudopotential choice and convergence
In the energy cutoff convergence study and PP selection, we compared the ultra-soft Rappe-Rabe-Kaxiras-Joannopoulos (RRKJUS) pseudopotential, 58 the hard norm-conserving PP Hartwigsen-Goedecker-Hutter (HGH), 59 and the recently developed ONCV as implemented in Quantum Espresso as well as the standard and hard PAW PP 60,61 as implemented in VASP. Fig. 2 displays the results for the case of PT barriers of the formic acid dimer. It shows that while ONCV, HGH, and hard PAW converge to a similar value of 0.12 eV, this value differs from the converged value of the ultrasoft RRKJUS and standard PAW, by 0.07 eV and 0.03 eV, respectively. In comparison, the literature value obtained with 6-311+G(3df,3pd) orbital basis set using Gaussian program package 62 is 0.13 eV. 1 The similarity of the ONCV, HGH, and hard PAW results demonstrates the reliability of these PPs. Our results are in line with the fact that ONCV has been shown to perform well compared to all-electron results for a selected set of solids including covalent, ionic, and metallic bonding. 56 In contrast to ONCV and hard PW which converge smoothly to within 1 meV at 1000 eV, the PT barriers obtained with HGH fluctuates significantly with the cutoff until the energy differences fall within 1 meV at an energy cutoff of 2800 eV.

FIG. 3: Mean (absolute) deviation (M(A)D) for PT barrier energies of the PT14 benchmarking set. The bar
colors indicate the nature of the respective exchange functionals, with GGAs indicated in gray, hybrids in blue, and meta-GGAs in pink. Stripes indicate the inclusion of dispersion forces, either with a dispersion correction (red stripes), or using full non-local correlation within vdW-DF (indigo) or at the rVV10 (light blue) level. ("vdW-" is dropped in the functional names. Fig. 3 displays computed statistical data for PT barrier energies for the PT14 set for a subset of the bench-marked functionals. The full set of results are provided in Tab. I of the appendix. Comparing the functionals, we find that the two pure GGAs, revPBE 47 and PBE, 46 (indicated by plain grey bars) significantly underestimate the PT barriers. For systems with low PT barriers, underestimating barriers can give qualitatively incorrect results. For instance, in the case of the malonaldehyde (system 1 in Fig. 1), PBE predicts 8 meV far less than the reference value of 168 meV. PBE0, which mixes in 25% Fock exchange with the PBE exchange, 15 improves the PT barriers. The meta-GGA SCAN is overall more accurate than PBE, but less accurate than PBE0. Despite the fact that dispersion forces contribute to hydrogen bonding, adding dispersion corrections at the DFT-D3 23 or at the TS 20 level to PBE has almost no impact on the predicted PT barrier. Interestingly, a similar insensitivity is also exhibited with the inclusion of rVV10 non-local correlation corrections to SCAN. 49 In contrast, several vdW-DFs, in particular vdW-DF-cx, 29,30 vdW-DF, 26 and vdW-DF2 27 have significantly smaller deviations from the reference than the GGAs. The most accurate non-hybrid functional is vdW-DF2 with a mean absolute deviation (MAD) of 0.073 eV, only a quarter of that of PBE of 0.279 eV. In comparison, the MAD of PBE0 is 0.094 eV. The reduced deviation of vdW-DF can be traced to the GGA gradient correlations having a larger negative contribution to the barrier than the fully non-local vdW-DF correlation. For instance, going from revPBE to vdW-DF, which keeps the exchange fixed, causes MAD to drop from 0.226 eV to 0.104 eV. As both of the non-local correlation and hybrid exchange increase PT barriers, the most accurate functional is obtained with the hybrid vdW-DF-cx with 20% Fock exchange (vdW-DF-cx0-20), 40,41 with a MAD of 0.037 eV, while the variant using 25% Fock exchange (vdW-DF-cx0) gives a MAD of 0.053 eV. Finally, we note the empirical B3LYP functional, 16 which also uses 20% Fock exchange and a reduced GGA correlation (0.81 LYP 63 ) also provides accurate PT barriers, with a MAD of 0.041 eV. Fig. 4 shows the individual correlation components for eight different functionals for the PT barrier of the formic acid dimer. It shows how the non-local correlation contribution of vdW-DF reduces the PT barrier less than the GGA correlation, thus making vdW-DF and vdW-DF2 the most accurate functionals for this PT barrier. As noted before, the comparison is most evident when comparing revPBE and vdW-DF as the exchange is kept fixed. The D3 correction in PBE-D3 and the rVV10 nonlocal correlation contribution in SCAN-rVV10 are barely visible.

IV. ANALYSIS OF ENERGETIC CONTRIBUTIONS OF VDW-DF
The dramatically improved performance of vdW-DF and vdW-DF2 compared to the GGAs, and the excellent performance of hybrid vdW-DF-cx prompted us to perform an analysis of the improved PT barriers. We retain the double-PT formic acid dimer (system 2 in Fig. 1) as a case study.

A. Exchange contribution
The degree of underestimation with GGAs as well as the performance of individual vdW-DF variants depends on the specific exchange functional used. Jenkins et al. 65 recently highlighted how the gradient component of the sresolved exchange energy can be a useful tool to analyze why the choice of exchange enhancement factor F x (s) causes the performance of different vdW-DFs to differ for different types of systems. The GGA exchange itself is expressed as modulation of the LDA exchange, as follows with the reduced gradient given by s(r) = |∇n(r)|/2(3π 2 ) 1/3 n(r) 4/3 .
In turn, the s-resolved The contribution to the PT energy is given by the difference between the transition state (ts) and ground state (gs) as follows: The s-integrated exchange energy density is given by Numerically, we computed Eq. 5 using the Savitzky-Golay 66 filter fitted with a third-degree polynomial to remove the noise imposed due to gridpoint integration and then take the derivative to obtain Eq. 4.
The result is shown in Fig. 5. Note that this analysis only takes the explicit energetic contributions of Eq. 3 into account, neglecting contributions from selfconsistency which can be considerable. 33,65 We can roughly discriminate two main s regions: 0.4 < s < 0.75 which reduces the PT barrier, and s > 0.75 which increases the barrier. Comparing revPBE to PBE exchange, we can trace the larger PT barrier of the former, to the increasing difference between the values of F x (s) as s increases, and that most of the difference comes from regions with s > 1.0. It is also interesting that vdW-DF and vdW-DF2 end up with quite similar energies, even though the exchange functionals, revPBE and PW86r, differ considerably in shape. This result is due to a partial cancellation of positive and negative exchange contributions to the PT barrier. Comparing the enhancement factors, F x (s), (lower panel Fig. 5) with ∆E x (s) reveals that the larger value of the PW86r F x (s) coincides with a larger slope of ∆E x (s) in most of the negatively and positively-contributing s regions (top panel). The final difference also partially cancels with the small difference in non-local correlation contributions. With a MAD of 0.144 eV, vdW-DF-cx performs more accurately than vdW-DF2-B86R with a MAD of 0.254 eV. This result can be related both to slightly smaller values of F x (s) in the regime that contributes to lowering the barrier, as well as considerably larger values beyond s > 1.0. The full effect of this is partly counteracted by self consistency which influences the kinetic energy considerably in Fig. 5. Tab. I shows that while several vdW-DFs perform better than GGAs, not all vdW-DFs are equally accurate. For instance, PBE performs better than the recently developed vdW-DF3 variants; however, this is not due to the non-local correlation, but rather their "soft" exchange functionals, i.e. the F x (s) shape has overall lower values. In fact, truly "soft" GGAs such as PBEsol 67 are even less accurate with a MAD of 0.4 eV. The "soft" form in the small-to-medium sregime is crucial for making vdW-DF accurate for other classes of systems such as coinage metals and layered systems. 30,33,36,68 In this sense, it is encouraging that vdW-DF-cx, which can be labeled "soft", performs relatively well. Moreover, the most accurate functionals are the hybrid variants of vdW-DF-cx, in line with accurate results found for other chemical reactions. 40,41

B. Analysis of the role of non-local correlation
For analyzing vdW-DF, the local responsivity q 0 (r) is the most natural variable as it enters into the non-local correlation kernel φ(r 1 , r 2 ) = φ(D, δ) through dimensionless parameters δ and D. In Sec. IV C, we will also project non-local correlation density onto s(r), for sake of comparison with GGA correlation. The parameter δ is the relative difference in responsivity q 0 (r) (inverse length scale) of two density regions δ = |q 0 (r 1 ) − q 0 (r 2 )| q 0 (r 1 ) + q 0 (r 2 ) , and the effective dimensionless separation D is is given by The parameter q 0 (r) is, within vdW-DF, given by where k F (r) is the local Fermi vector. Z ab is equal to −0.849 for vdW-DF1 and vdW-DF3-opt1 and −1.887 for vdW-DF2 and vdW-DF3-opt2. The spatial non-local correlation density is given by and, subsequently, the q 0 -resolved non-local correlation can be defined as The q 0 -resolved non-local correlation contribution to the PT barrier is then given as and the corresponding q 0 -integrated non-local correlation barrier is For the analysis of the non-local correlation contributions to PT barriers, Fig. 6 shows four interlinked panels. Panel (a) shows three selected q 0 isosurfaces for the ground state (gs) and transition state (ts). The overlaid contours indicate the non-local correlation density as given by Eq. 9. The three isosurfaces (q 0 ≈ 1. transition state, while the ground state isosurface has been disconnected. These isosurfaces differences explain the increase in the positive non-local correlation energy of the transition state, thus increasing the magnitude of the PT barrier. Finally, the upper right isosurface corresponds to q 0 = 2.2 in which the isosurface lobes around the hydrogen in the transition state has vanished, but the isosurfaces around the hydrogen in the ground state cause a significant reduction of the PT barrier. At this value of q 0 , the effective separation D has contracted significantly causing intra-molecular correlation to be dominating contribution.

C. Reduced-gradient resolved correlation comparison
The previous subsection highlighted the utility of analyzing vdW-DF in terms of the spatial distribution of q 0 (r) which enters into the kernel φ(D, δ) of vdW-DF. In GGA correlation, 46 a key variable is the reduced gradient t = |∇n(r)|/2k s (r)n(r), which is defined using Thomas-Fermi screening length k s = 4k F /π. In terms of t, the total PBE correlation energy can be expressed where H is the gradient contribution function. In order to analyze the effect of correlation and exchange on an equal footing, whether at the GGA or vdW-DF level, we benefit from resolving these quantities onto the same variable. We here choose the reduced gradient s which was used for exchange in Fig. 5. The corresponding sresolved differential PBE correlation energy is given by e PBE c (s) = d 3 r n(r) H(r, t(s)) δ(s − s(r)) , with the corresponding integrated quantity given by The projection of vdW-DF onto s is similar to Eq. 10. Fig. 7 shows the result for the PT energy barrier contribution ∆E PBE c (s) and the non-local correlation. Comparing the curves with the exchange part of Fig. 5 shows that the curve bear some resemblance, but with opposite prefactor and earlier onset of correlation contribution than exchange with increasing s. Comparing GGA and vdW-DF correlation shows that while the shapes are similar, the former has much larger positive and negative contributions beyond s ≈ 0.4. Moreover beyond s ≈ 1.5, vdW-DF correlation flattens, while the magnitude of GGA correlation continues to grow. Inspecting the isosurfaces, we found that the large s > 1 values corresponds to large isosurfaces similar to the low q 0 values.

V. CONCLUSIONS AND OUTLOOK
We have compared and analyzed the performance of different functionals for proton-transfer energy barriers and found that using non-local correlation in vdW-DF rather than GGA correlation causes a non-intuitive lowering of the energy barriers, typically improving accuracy compared to GGA, which tends to underestimate energy barriers. The best performance was provided by a hybrid version of the consistent-exchange van der Waals density functional vdW-DF-cx with 20% Fock exchange.
The improved PT transfer barriers of vdW-DF is highly encouraging because these functional can account for dispersion forces responsible for the cohesion of many proton-transfer systems. Beyond this, our study points to the possibility that the usage of a non-local kernel reflects a true geometry-sensitive repulsive short-range correlation contribution. In this context, we emphasize that while vdW-DF was designed with dispersion forces in mind, the theory is rooted in exact constraints and many-body theory and therefore is not limited to long range dispersion forces as such. vdW-DF is a true nonlocal correlation functional with contributions both at short and long ranges. At the same time, the performance of approximative exchange-correlation functionals are generally contingent on the various simplifications and parametrizations used in their construction. The improved performance may also be a mostly "fortuitous" effect arising due to typically smaller semi-local correlation-type effects in vdW-DF. Exploring this question merits further theoretical and computational investigations and may pave the way for more accurate methods using non-local functionals for improving chemical reaction energies.

ACKNOWLEDGMENTS
The computations of this work were carried out on UNINETT Sigma2 high performance computing resources (grant NN9650K). This work is supported by the Research Council of Norway as a part of the Young Research Talent project FOX (302362). We thank Rasmus Tranås and Elin Dypvik Sødahl for comments. KB acknowledges discussions with T. Jenkins and T. Thonhauser on visualization.
Appendix A: Tabulated data for PT14 set Table I shows the computed data for all studied functionals. The upper part shows results for all but the vdW-DFs, which is shown in the lower part.