A Filon-like integration strategy for calculating exact exchange in periodic boundary conditions: a plane-wave DFT implementation

An efficient and accurate approach for calculating exact exchange and other two-electron integrals has been developed for periodic electronic structure methods. Traditional approaches used for integrating over the Brillouin zone in band structure calculations, e.g. trapezoidal or Monkhorst-Pack, are not accurate enough for two-electron integrals. This is because their integrands contain multiple singularities over the double integration of the Brillouin zone, which with simple integration methods lead to very inaccurate results. A common approach to this problem has been to replace the Coulomb interaction with a screened Coulomb interaction that removes singularities from the integrands in the two-electron integrals, albeit at the inelegance of having to introduce a screening factor which must precomputed or guessed. Instead of introducing screened Coulomb interactions in an ad hoc way, the method developed in this work derives an effective screened potential using a Filon-like integration approach that is based only on the lattice parameters. This approach overcomes the limitations of traditionally defined screened Coulomb interactions for calculating two-electron integrals, and makes chemistry many-body calculations tractable in periodic boundary conditions. This method has been applied to several systems for which conventional DFT methods do not work well, including the reaction pathways for the addition of H2 to phenol and Au20−\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$_{20}^{-}$\end{document} nanoparticle, and the electron transfer of a charge trapped state in the Fe(II) containing mica, annite.


Introduction
One of the more computationally demanding and important scientific simulations for materials and chemistry is the ab initio Molecular Dynamics (AIMD) method (Car and Parrinello 1985;Remler and Madden 1990;Payne et al. 1992;Marx and Hutter 2009;Bylaska et al. 2011b;Bylaska 2017) in which the motions of the atoms are simulated using Newton's laws where the forces on the atoms are obtained using the plane-wave DFT methods. Plane-wave DFT methods are advantageous for implementing AIMD because the calculation of the forces is computationally less expensive than more traditional LCAO DFT methods. An exciting use for AIMD is to be able to directly simulate chemical reactions in condensed phases. However, these types of simulations are rarely attempted, because even with an efficient plane-wave DFT method, AIMD requires large computational resources not accessible to most researchers. This situation is even worse for the direct simulation of many standard chemical reactions, because higher and more expensive levels of DFT, which explicitly include exact exchange (i.e. Hartree-Fock exchange), are needed to accurately represent the energetics of transition states (Bylaska et al. 2011a;Bylaska 2017).
A fundamental challenge in including exact exchange in plane-wave DFT methods for condensed phase systems is that one typically wants to use periodic boundary conditions to carry out the simulation. While this is natural for plane-wave DFT methods without exact exchange, it is significantly more complicated when using exact exchange as it requires special integration strategies to handle the integration of the Brillouin zone. At first glance, periodic many-body calculations would appear to be intractable because the expansion of one-electron orbitals in terms of Bloch states leads to a large number of orbitals describing the first Brillouin zone. Simple approximations to the integration over the Brillouin zone in the exact exchange and other two-electron integrals lead to very inaccurate results, e.g. a straightforward point approximated calculation results in the two-electron integrals being infinite.
To handle this problem the condensed phase community typically introduces a screened Coulomb interaction to calculate the two-electron integrals (Stampfl et al. 2001;Heyd et al. 2003;Heyd and Scuseria 2004;Peverati and Truhlar 2012). The types of screened interactions in these approaches are highly dependent on the type of material and are calculated using an expensive coupled perturbed DFT calculation or the screened interaction is directly optimized during the calculation as in an RPA (Langreth and Perdew 1977;Niquet et al. 2003) or GW (Hedin 1965;Aryasetiawan and Gunnarsson 1998;van Schilfgaarde et al. 2006) calculation. In principle, this type of approach can be made accurate. However, the highly engineered nature of these types of approaches, which introduce material dependent screenings, makes it difficult to use them for systems in which the screenings vary across the simulation cell, e.g. catalytic reactions at metal surfaces, adsorption at solvated mineral surfaces, the structure of complex brines containing highly charged ions, extended defects in solids, and electron transfer in solids and at complex interfaces.
Another approach that can be used that is less engineered is the Wannier orbital approach (Bylaska et al. 2011a;Bylaska 2017). This approach is preferable for molecular and liquid systems, and when used in combination with advanced high-performance computing algorithms it has been used with success on a variety of non-trivial systems, including highly charged ions in solution (Atta-Fynn et al. 2011;Fulton et al. 2012), mineral surfaces (Chen et al. 2016;2017), reactions in solution (Bylaska 2017), and electron transfer (Bylaska and Rosso 2018). Despite these successes, there are still some well known drawbacks. In particular, for periodic solids it is not obvious how to extend the Wannier orbital approach to systems that can be described using small unit cells and large k-point samplings (e.g. small band gap semiconductors), as it requires the use of large unit cells and large band gaps to be accurate.
In this work, we present an efficient and accurate approach for calculating exact exchange and other two-electron integrals in periodic electronic structure methods. In contrast to prior screened interaction and Wannier orbital based methods (Bylaska et al. 2011a), our new proposed method is equally applicable to both condensed phase systems and molecular systems and can easily be generalized to work with standard Monkhorst-Pack integration techniques used in band structure calculations. The manuscript first shows, in "Derivation of exact exchange in periodic boundary conditions" section, the derivation of exact exchange in periodic boundary conditions, and "Integration strategy for exact exchange in periodic boundary conditions" section presents a Filon-like integration (Filon 1930) approach for it, in which an effective screened or filtered exchange kernel is derived solely in terms of the size and shape of the first Brillouin zone. A key aspect of our derivation is the observation that the integrands in the exchange integrals are composed of a smooth function times a stiff function, where the stiff function is independent of the one-electron orbitals and is purely a function of first Brillouin zone. In "Integration strategies for generating filtered potentials" section, an accurate and efficient strategy for generating the filtered potential is presented, and in "Applications" section the new method is applied to Hybrid DFT and electron transfer calculations. Finally, conclusions are given in "Conclusion" section.

Derivation of exact exchange in periodic boundary conditions
The standard formula for exact exchange energy of a total system is given by where ψ n (r) are the one electron spin orbitals, V is the integration volume that spans the entire system, and N σ occ are the number of occupied spin orbitals. The notation for the volume with the bar above it, e.g.V , is used to denote the range of integration in volume integrals.
The following substitutions, can be used to convert E x to periodic boundary conditions, where V BZ = 8π 3 is the volume of the first Brillouin zone, is the volume of the periodic unit cell, and ∞ is an infinite collection of supercells, with volume , separated by Bravais lattice vectors a.
Using these substitutions the exact exchange energy can be transformed into the following energy per unit cell The following relation for periodic sums can be used to simplify the exchange energy to be In this formula, ψ nk (r) and a exp(iK·(R−a)) |R−a| are periodic functions that can be represented by the following discrete Fourier sums Making these replacements and then integrating out r and r results in the following formula is the volume of the first Brillouin zone and the overlap densities in reciprocal space are given by

Integration strategy for exact exchange in periodic boundary conditions
As pointed out by Gygi and Balderechi (1985, 1986 and others Chawla and Voth (1998); Sorouri et al. (2006); Marsman et al. (2008); Görling (1996), this expression must be evaluated with some care especially for small Brillouin zone samplings and small unit cell size, because of the singularities at |G − k + l| = 0.
In previous work, we have shown that the exact exchange energy for periodic crystals can be implemented using a formalism based on localized Wannier orbitals (Marzari and Vanderbilt 1997;Silvestrelli 1999). Extensive derivations of formulae implementing exact exchange can be found in prior work by (Bylaska et al. 2011a). In this work, we show how to simplify Equation 7 so that it can be evaluated with the same level of accuracy as the localized Wannier orbital procedure (Bylaska et al. 2011a, b), and arguably with even more accuracy.
The first step is to notice that if the Brillouin zone sampling or unit cell is large enough, the overlap densities are independent of the Brillouin zone integration in the exchange energy (Bylaska 2015). That is, Equation 7 can be rearranged as follows, where V f is given by We refer to this potential as a filtered potential. Moreover this strategy can be generalized to be a 2D trapezoidal integration over k and l, where there is a filtered potential for each 2D patch in the trapezoidal integration. Our new approach is named after Filon integration because of its similarity to this style of integration, in which the integration (over k and l) of a slowly varying function that is modulated by a high frequency or stiff function is carried out. In the case of the exchange integral the slowly varying functions are the overlap densities and the stiff function is 4π |G−k+l| . The filtered potential can readily be evaluated for most of G vectors by using a simple numerical integration. The range of integration is illustrated in Fig. 1. However, for G = 0 and the G vectors just outside the first Brillouin zone, care must be taken to evaluate these integrals. The first term in Equation 11 can be evaluated by using numerical integration since the singularity at |l−k| = 0 has been removed, and the second term in this formula is approximated by expanding its range of integration to an infinite volume, and using the following formula i.e.
For simple cubic boxes, the values of the integral where the integrand is singular can be shown to be

Integration strategies for generating filtered potentials
In order to make the generation of the filtered potentials manageable, the following integrals, (for G vectors not inside the first and second Brillouin zone), and (for G vectors inside the first and second Brillouin zone), which all have positive values, have to be integrated accurately and efficiently. The integrands in the I 1 and I 2 integrals are smooth and can be numerically integrated using standard trapezoidal, Simpson, or other polynomial integration strategies. However, these integrals are expensive to evaluate, because they are six-dimensional. The I 3 integrals are even more difficult to integrate, because not only are they six-dimensional, they also have very stiff integrands (they are singular at |G+l−k| = 0) and cannot be numerically integrated with standard integration strategies accurately. Another complication in evaluating these integrals is the integration volumes (i.e. the Brillouin zone) may be trapezoidal shaped rather than cubic. This complication can be dealt with by introducing the scaled coordinates (i.e. fractional coordinates) ξ 1 , ξ 2 , ξ 3 for k, and ζ 1 , ζ 2 ζ 3 for l. These coordinates are defined by the following linear transformations where b 1 , b 2 , and b 3 are the reciprocal lattice vectors, and the Jacobian of these two trans- Using this linear transformation, the above six-dimensional integrals can be recast as and γ i = ξ i − ζ i for i = 1, 2, 3. We note that in the above and subsequent labeling of the functions I 1 , I 2 , and I 3 , the explicit inclusion of G in them has been removed to make the formulas more readable.

Reducing the six-Dimensional I 1 and I 2 integrals to three-Dimensional integrals
The integrals I 1 and I 2 can easily be reduced to three-dimensional integrals by introducing the following linear coordinate transformations where i = 1, 2, 3. The Jacobians for each component of this transformation is |J| = 1 2 and the limits of integration in the transformed variables are shown in Fig. 2. Using this transformation, the integration over each ξ i and ζ i pair is converted to and combining these transformations for each pair results in the overall integration of I 1 and I 2 being simplified to the following three-dimensional integrals The I 1 integral can be rewritten for large G, i.e. outside the first Brillouin zone, into the following generalized multipole form where P l,|m| (cos θ) is an associated Legendre polynomial (with (−1) m factor omitted) and x = (cos φ sin θ, sin φ sin θ, cos θ). This formula is straightforward to derive by using the expansion of Gegenbauer polynomials in terms of Legendre polynomials (Rainville 1960). This form can be used to simplify the computation of the I 1 integrals needed to compute the filtered potentials. Also, with this form it is easy to show that for larger unit cells (i.e. smaller reciprocal cells) the I 1 integrals essentially reduces down to be just 4π G 2 .

Changing the ranges of integration and using spherical coordinates to simplify the I 3 integrals
By shifting the center of the I 3 integrals to be at G (for G vectors inside the first and second Brillouin zone), they can be recast as double volume integrals where the volumes span different regions of space, i.e. the volumes are located on top of each other or next to each other as shown in Fig. 3. Note, this change in the center of integration produces in 27 different integrals.
For example, the ranges of integration for the 6 face sharing volumes are replaced by changing the ranges of integration for one pair of {ξ i , ζ i } to be from e.g.
Similarly, the ranges of integration for the 12 edge-sharing and 8 corner-sharing volumes are obtained by replacing two and three pairs of {ξ i , ζ i } respectively, e.g. and For these alternative ranges of integration, we again transform from ξ i and ζ i variables to γ i and η i variables using The Jacobian for this transformation is |J| = 1 2 and the limits of integration in the transformed variables will be These integrals can be simplified further by noting the decay parameter α is chosen Combining these transformations with the transformations in Eq. 26 allows one to simplify the various I 3 integrals to be three-dimensional integrals, where the 1 same integral is the 12 edge-sharing integrals are I edge,(θ=0...π,φ=0...π/2) 3,q=(1,1,0) and the 8 corner-sharing integrals are I corner,(θ=0...π/2,φ=0...π/2) 3,q=(1,1,1) Even though these integrals have been reduced to three-dimensions, they still cannot be integrated accurately in their current form, because the integrand f 3 is singular ( 1/R 2 ) at γ 1 = γ 2 = γ 3 = 0. This singularity in the integration can be removed by transforming to γ 1 , γ 2 , γ 3 to spherical coordinates using By changing the integration variables, the I 3 integrals can be written as where γ i (r, θ, φ) are now functions of r, θ and φ; R 2 (γ 1 , γ 2 , γ 3 ) is evaluated using Eq. 25 with G = 0; and the function w(γ 1 , γ 2 , γ 3 ) is used to represent the extra factors resulting from changing the integration over ξ i and ζ i to γ i . A key result of this change in integration is that the integrand is multiplied by r 2 , and since R 2 is proportional to r 2 the transformed integrand, r 2 f 3 w is no longer singular and thus it can be integrated accurately using numerical integration. To remove the 1 r 2 factor in I 3 it is convenient to introduce the following re-scaled γ i ,γ

Comparisons and approximations to filtered exchange potentials
In Figs. 5 and 6, various periodic filtered exchange potentials and their errors are shown for a L = 20 simple cubic unit cell. The V f 1 potential (solid black line) is the most accurate of the potentials shown and it was generated with formulae for the I 1 , I 2 and I 3 integrals. Even though this potential can be precomputed for use in a plane-wave code, the cost of doing these calculations can become quite expensive for large Brillouin zones (O(N 6 ), where N is number of grid points along each dimension of the 3D-FFT used in plane-wave DFT methods), making it undesirable, especially for unit cell optimizations in which the potentials are constantly being recomputed. This cost can be reduced to essentially the same as the 4π |G| 2 potential (O(N 3 )), by replacing the I 1 integrals with the more efficient generalized multipole form of Eq. 29. This potential (V f 2 in Figs. 5 and 6 (dashed red line)), is shown to produce nearly the same results as V f 1 . Making it reasonable to use in a plane-wave DFT method.
The essential character of these potentials are that they mimic a 1/|R| potential within a single unit cell, and then slightly flatten off near the unit cell boundary. For comparison, the 1/|R| potential (for a simple cubic cell V (R = 0) = 2.380077/h where h = L/N FFT with L being the side length) and traditional periodic Coulomb kernel V f 3 = 4π |G| 2 (G = 0) (dashed green line), as well as slight modifications to it, V f 4 (dashed blue line), and V f 5 (dashed indigo line), are shown in Figs. 5 and 6. The V f 4 and V f 5 potentials were modified by replacing the G vectors in the first Brillouin zone and G = 0, respectively by I 2 + I 3 . Note V f 4 is the same as V f 2 with n max = 0. As can be seen from the plots the traditional periodic form, V f 3 , has a large error and produces a kernel in real-space that is significantly below V f 1 . The V f 4 and V f 5 potentials are considerably better and appear to produce the correct long range behavior which is close to 1/|R|, however, the overall errors are still fairly significant. Finally, the V f 6 potential (dashed dark green line), i.e.
is the cutoff Coulomb kernel used in our prior exchange paper based on the Wannier orbitals. The design of this cutoff kernel is chosen to remove the interactions between redundant periodic images of Wannier orbitals, because of the long-range nature of the Coulomb potential. However, it should be noted that below cutoff radius the cutoff Coulomb kernel, while tracking the other kernels, has slightly larger values which become more pronounced for larger R. As a result the Wannier kernel will produce slightly larger exchange energies for systems in which a localization procedure produces very localized  R) are plotted in real-space along the x-axis for an FFT=32x32x32 L = 20 simple cubic unit cell. The Cartesian I 1 , I 2 and mulipole integrals (n max = 10) were computed using Simpson integration with a 228x228x228 grid in γ 1 , γ 2 and γ 3 . The radial I 3 integrals were computed using Simpson integration with a 432x432x432 grid in r, θ and φ. The parameters for V f 6 were chosen to be R cut = 8 and N = 8. In the figure, V f 1 is the most accurate definition of the kernel, V f 2 is generated with a more efficient generalized multipole form, V f 3 is the traditional periodic Coulomb kernel, V f 4 and V f 5 are further approximations to V f 2 , V f 6 is the cutoff Coulomb kernel, and 1/R is the simple Coulomb kernel for free space boundary conditions wavefunctions, i.e. systems with large band gaps. Also the V f 6 potential appears to be close to the V f 4 and V f 5 potentials for |R| < R cut .

Applications
All calculations in this study were performed using the pseudopotential plane-wave program (NWPW module) contained in the NWChem computational chemistry package (Valiev et al. 2010). For the DFT calculations, both the gradient corrected PBE96 (Perdew et al. 1996)  Hartree-Fock (UHF) broken-symmetry wavefunctions were used to represent the diabatic wavefunctions. The valence electron interactions with the atomic core are approximated with generalized norm-conserving Hamann (1989) or Troullier-Martins (1991) pseudopotentials modified to the separable form suggested by Kleinman and Bylander (1982). A brief description of the pseudopotentials used in this study is given in Appendix 4. The exchange integrals in the hybrid-DFT and electron transfer calculations were performed using the filtered exchange kernel described above with n max = 10, and also with the Wannier orbital based screened Coulomb kernel from the prior work of Bylaska et al. (Bylaska et al. 2011a) with parameters R cut = 8 Bohrs and N = 8.

Inclusion of exchange in organic reactions: addition of halogens to the carbon-carbon double bond
We begin our testing of the filtered exchange potential by examining the reaction pathway for the electrophillic addition of a Cl 2 molecule onto the carbon-carbon double bond in ethene. This well known reaction is extremely slow and to make it feasible it is often catalyzed with an Iron(III) chloride solid. It is also important in production of 1,2-dichloroethane, which is an intermediate for other organic compounds such as vinyl chloride. As a well known reaction it is one of the first reactions students learn about in organic chemistry. While seemingly simple, its reaction pathway is surprising to students and a challenge to model. The addition of halogens to alkenes is believed to involve two steps, where the first involves the formation of chloriran-1-ium, a cyclic chloronium cation (SMILES: C1C[Cl+]1), and a chloride, and then is subsequently followed by the addition of the chloride (Morrison and Boyd 1992). However, it should be noted that in the gas phase the reaction energy to form the chloriran-1-ium and chloride is well over 100 kcal/mol and unlikely to occur. Prior calculations by Kurosaki (2000) have shown that the lowest transition state of the gas-phase reaction is the extraction of a radical Cl from Cl 2 by C 2 H 4 , after which the radicals recombine to form C 2 H 4 Cl 2 .
The foremost challenge to modeling this type of reaction with a periodic plane-wave program has been that standard DFT methods significantly underestimate the reaction barrier for this type of reaction and are not reliable. Hybrid DFT methods are thought to produce higher and more accurate barrier heights (Nachtigall et al. 1996;Andzelm et al. 1995;Zhao et al. 2005). Moreover, this type of basic organic reaction, in which the bonding is primarily covalent, is expected to be well described by the Wannier exchange methods, and is thus a good test of the consistency of our newly developed exchange potential with prior exact exchange implementations for isolated systems.
The reaction pathway for the addition of Cl 2 to the double bond in ethene (C 2 H 4 ) is shown in Figs. 7 and 8. As expected, the hybrid-DFT PBE0 calculations produce a reaction pathway with a higher barrier by ≈10 kcal/mol than the regular DFT PBE calculations, and the filtered periodic and Wannier exact exchange kernels produced nearly identical reaction energy pathways. These calculations were carried out in a periodic FCC unit cell with L=20.109 Å at a cutoff energy of 100 Ry. Standard pseudopotentials contained in the NWChem program package (see Appendix 4) and both Wannier and filtered periodic exact exchange kernels were used in these calculations. These potential energy surfaces were computed with the PBE and PBE0 exchange-correlation potentials.
In Fig. 7, the pathway results are shown from a series of optimizations using a "bondings" penalty function, where the collective variable (or reaction coordinate) γ is defined as the difference squared between the sum of bonds broken and the sum of the bonds made, i.e.
and γ 0 are constants used to target the reaction coordinate. The target values for the reaction coordinate γ 0 were chosen to be between -11.0 Bohr 2 and 1.0 Bohr 2 , and the spring constant was set to K = 1 Hartree/Bohr 2 .
To further quantify this reaction we also determined saddle points and their corresponding reaction pathways using the intrinsic reaction coordinate (IRC). The IRC is defined as the minimum energy reaction pathway in mass-weighted Cartesian coordinates between the transition state and the reactants and products (Fukui 1981). In these calculations, the saddle points were optimized using the Sella optimization package Intrinsic reaction coordinate (IRC) calculations for the ethene + Cl 2 → 1,2-dichloroethane addition reaction from periodic plane-wave DFT calculations and non-periodic Gaussian DFT calculations. Results using Wannier and filtered periodic exact exchange kernels are shown for the plane-wave hybrid-DFT PBE0 calculations. The hybrid-DFT PBE0 calculations contain 0.25 exact exchange. One transition state (TS1) was found with the DFT PBE calculations, and two different transition states (TS1 and TS2) were found in the hybrid-DFT PBE0 calculations. The geometries are illustrated as TS1 and TS2 for the two different transition states duced nearly identical results between the filtered periodic and Wannier exact exchange kernels, but also showed the plane-wave and Gaussian basis set calculations agreed very well with one another.
Besides producing a slightly higher reaction barriers, these more extensive calculations also showed that the PBE and PBE0 calculations produced subtly different results in that only one transition state was found with the PBE calculations, whereas two transition states were found with the PBE0 calculations. The first transition state had the same structure as was found with the PBE and penalty function method calculations, where one of the chlorine atoms is preferentially bound to one carbon atom and the other is slight dissociated. The second transition state, which is slightly higher in energy, has a structure where the chlorine atoms are loosely bound to each carbon atom.

Band gaps of insulators
For our second application, we calculate the band gaps of two well known insulators using DFT and hybrid-DFT. Several studies have shown that traditional DFT methods, or density only exchange-correlations methods such as PBE96, significantly underestimate the band gaps for semiconductors and insulators by up to 60%. This large error has been attributed to the fact that density only exchange-correlation functionals do not have a discontinuity at the Fermi level, because the effective potentials for the occupied and unoccupied states in them are the same rather than different.
In Table 1, the band gaps of Al 2 O 3 and SiO 2 at the PBE96 DFT and hybrid PBE0 DFT levels are reported. In these calculations a cutoff energy of 100 Ry was used, and the sizes of the supercells were 80 and 72 atoms respectively for the two crystals. The default pseudopotentials of NWChem were used (see Appendix 4). These results show that band gaps at the hybrid-DFT level are considerably better than those at the conventional DFT level, and that the filtered periodic kernel produced slightly lower band gaps then with the Wannier exact exchange kernel. This is not unexpected since the Wannier kernel will produce a slightly lower HOMO (highest occupied molecular orbital) energy and hence a slightly larger band gap. The band gaps were estimated by taking the difference between the HOMO and LUMO (lowest unoccupied orbital) single particle energies at the point using large supercells. With this approach only the minimal gaps at the point or at any point folded into the point were obtained. However, since large unit cells were used, many of the high-symmetry points in the Brillouin zone were covered.

Hybrid DFT calculations for the adsorption reaction of Au − 20 + H 2
For our next application, we examine the reaction pathway for the adsorption of H 2 onto the top site of tetrahedral Au −1 20 nanoparticle. Neutral and negatively charged gold nanoparticles have gained interest in the catalysis community, because of their potential catalytic activity and selectivity for a variety of chemical reactions important to industry and environmental remediation, including anaerobic oxidation of methanol, CO oxidation, NO x reduction, and selective hydrogenation of unsaturated hydrocarbons including acrolein and nitroaromatic compounds. The selectivity seen in hydrogenation reactions appears to be unsystematic. For example, the critical first step of H 2 adsorption produces different trends as a function of nanoparticle size for the elementary reaction energies and reaction barriers; where the reaction of Au 13 +H 2 is exothermic with a modest barrier, Au 20 +H 2 is slightly endothermic with a high barrier, and Au 55 +H 2 is very endothermic with a barrier in between Au 13 and Au 20 . The reaction pathway for the adsorption of H 2 unto the top atom of the tetrahedral Au − 20 nanoparticle is shown in Fig. 9. These calculations were carried out in a periodic FCC unit cell with L=22.695 Å at a cutoff energy of 100 Ry. Both Wannier and filtered periodic exact exchange kernels were used in these calculations. The potential energy surface was computed with the PBE and PBE0 exchange-correlation potentials using a "bondings" penalty function, where the collective variable (or reaction coordinate) γ is defined as the difference squared between the sum of bonds broken and the sum of the bonds made, i.e.
and γ 0 are constants used to target the reaction coordinate. The target values for the reaction coordinate γ 0 were chosen to be between -8.0 Bohr 2 and -1.50 Bohr 2 , and the spring constant was set to K = 1 Hartree/Bohr 2 . Standard DFT methods are known to be unreliable for this type of reaction, because it involves a neutral closed shell molecule interacting with a radical. Not surprisingly, the hybrid-DFT PBE0 calculations produce different reaction pathways than the regular DFT PBE calculations. In addition, both the filtered periodic and Wannier exact exchange kernels produced nearly identical reaction energy pathways, which suggests that both kernels are equally valid for computing relative energies in large unit cell calculations.
However, somewhat surprising is the dramatic nature of the change of the curves at large distances, which show that the hybrid DFT method with exact exchange (PBE0) predicts an exothermic physi-adsorbed state and a near thermodynamic neutral chemiadsorbed state, whereas the DFT method (PBE) only predicts an endothermic chemiadsorbed state. At this point, it is uncertain which curve better represents reality. Future work will focus on using other approaches for estimating the correlation energy such the random phase approximation (RPA) and Muliti-Configurational Self-Consistent Field (MCSCF) methods. We note, the filtered kernels presented in this work are currently being used in the development periodic RPA and MCSCF methods.

Fe(II)/Fe(II) electron transfer calculations in annite
For our last application, we examine the electron transfer (ET) in annite, an Fe(II) containing trioctahedral mica. The structure of this mica, shown in Fig. 10, contains Fe(II) octahedral layers sandwiched between two (Si,Al) tetrahedral layers. During oxidation, electrons may be removed from these Fe(II) ions, most likely resulting in localized Fe(III) charge states that are expected to be mobile in which the hole on the Fe(III) hops to a neighboring Fe(II) site (Rosso and Ilton 2003;2005). This type of ET is conceptually quite simple as it does not involve any bond breaking, complicated changes in electronic states, or any long-range rearrangement of atoms typically found in ET reactions occurring in a polar solvent. However, even for such a simple ET reaction, it still requires a "beyond DFT method" to describe the localization of positive charge in the d-bands of the material and its rate of transfer between Fe(II) ions. It should be noted that a major failing of commonly used DFT functionals, such as PBE96, is that they cannot describe these types of localized Fe(III) charge states. This is a well known failure, and it is a by-product of these functionals containing self-interaction, which creates an artificial Coulomb barrier to charge localization, even for systems where significant charge localization is expected.
One of the more computationally accessible (beyond DFT) ET methods available today is the approach developed by Farazdel et al. (1990), and used extensively by ; Rosso and Ilton (2003); Rosso and Ilton (2005) Kerisit and Rosso (2006); Kerisit and Rosso (2005); Stack et al. (2004); Wang et al. (2008). This ET method uses key ideas from semiclassical Landau-Zener theory for ET along with valence bond theory to calculate the V AB electronic coupling matrix element. This method has been extended to periodic systems Bylaska and Rosso (2018). A caveat of this method is that it requires an accurate calculation of exact exchange with periodic boundary conditions. In this previous work we used a Wannier orbital based method to calculate periodic exact change. We note that very recently, Behara and Dupuis have also been exploring the use of approximate forms of this theory based on DFT+U method (Dudarev et al. 1998;Anisimov et al. 1991), which does not need accurate exchange calculations, with some success (Behara and Dupuis 2019). In this section, we demonstrate the use of the filtered exchange potential with this ET method for hole transport in annite.
For the transition metal oxide systems, where the ET typically occurs by electron or hole small polaron hopping, the Landau-Zener equation can be written (Holstein 1959a; Holstein 1959b) in terms of V AB , temperature (T), transition state energy ( G † ), and a typical longitudinal phonon energy ( ν 0 ) as where P A→B (Q C ) is the probability of the system not transitioning (i.e. staying on the adiabatic surface) from the ground-state to the first excited state at the transition state geometry, Q C , between the reactant state A and the product state B. Typically, this formula is extended to account for multiple crossings and re-crossings through the intersection region to properly describe this electron transmission. This modified transmission factor κ el , which is sometimes called the electron factor, is the probability of ET after the system has evolved to Q C , and is given by Newton and Sutin (1984) A nice feature of the semiclassical Landau-Zener approach for estimating non-adiabatic ET rates k non−adiabatic is that these rates relate to the adiabatic rates k adiabatic as This approximate relation, which becomes more correct when nuclear tunneling and the initial non-equilibrium effects are small (Garcia-Viloca et al. 2004;Reimers et al. 2015), allows one to make use of standard molecular modeling techniques appropriate for the adiabatic ET and then estimate the non-adiabatic rate by scaling it down by κ el . Figure 11 and Table 2 show the results for the ET between two neighboring cis-Fe in annite along with comparisons to prior ET calculations with small cluster models by Rosso and Ilton (Rosso and Ilton 2003). For these calculations, internal energy differences R is the distance between the iron atom with a localized electron and its neighbor. λ/4 is the Marcus estimation of the diabatic activation energy. E † diabatic is diabatic activation energy. E * adiabatic is the adiabatic activation energy. V AB is the electronic coupling. P A→B is the adiabaticity computed using ν 0 = 0.012 eV. κ elc is the steady-state transmission factor.
( E) rather than free energies differences ( G) were used. The differences of the enthalpy and entropy corrections due to vibrations orthogonal to the generalized reaction path are expected to be small in this system.
In these calculations, the supercell contained 88 atoms (K 4 Fe 12 Si 12 Al 4 O 48 H 8 ) with lattice parameters determined using the PBE96 exchange-correlation functional. The plane-wave calculations used a cutoff energy of 100 Ry and the ion-electron interactions were described using a small core relativistic pseudopotential for Fe. It was shown in the previous work by Bylaska and Rosso (2018) that a small core pseudopotential was needed to produce reliable results, as a large core pseudopotential for Fe predicted overlaps and energy splittings that were too large (i.e., more adiabatic). We used open-shell spin-unrestricted Hartree-Fock (UHF) broken-symmetry wavefunctions to represent the Fig. 11 Calculated diabatic and adiabatic potential energy surfaces for electron transfer between two nearest neighbor cis-Fe octahedra in annite diabatic wave functions in the calculation of ET matrix elements H AB and S AB . We note that MCSCF and other approaches such as constrained DFT and localization methods could also be used to define diabatic wave functions, but they were not used in these calculations. To localize the hole on any of the Fe atoms in the cell a spin-penalty function was used.
A plot of the potential energies for the diabatic and adiabatic states of A (Q) and B (Q) (i.e., H AA and H BB ) as a function of the ET reaction coordinate, Q, are shown in Figure 11. To calculate the valence bond diabatic states as a function of Q, the geometry at Q A was taken to be the optimized geometry in which the hole was localized on a cis-Fe, and the geometry at the other end point Q B was defined by localizing the hole on a neighboring cis-Fe. As expected, the energies of Q A and Q B were nearly the same. The geometries between Q A and Q B were then defined using a simple linearized reaction coordinate (the slight error the optimization of Q A and Q B is not expected to affect the reaction coordinate at the transition state). The adiabatic curves were then determined by calculating the matrix elements H AA , H BB , H AB , and S AB at each value of Q and then diagonalizing the 2 × 2 generalized eigenvalue problem (See Bylaska and Rosso (2018)) for details on the computational procedure). The Landau-Zener equation was used to calculate the probability P A→B (Q C ) of the system not transitioning (i.e., staying on the adiabatic surface) from the ground-state to the first excited state at the transition state geometry, Q C , between the reactant state A and the product state B at which the electron is transferred.
As can be seen from the results, the nearest-neighbor ET between cis-Fe produced a large splitting, and as a result, the ET is predicted to be adiabatic. Our work agrees with the assertion by Sherman (1987) and the previous results of Rosso and Ilton based on cluster models that have predicted this type of ET to be adiabatic (Rosso and Ilton 2003). However, the adiabatic splitting for the ET at the transition state was found in our calculations to be 0.102 eV, whereas Rosso and Ilton found it to be 0.064 eV. The difference in results is somewhat larger than seen in previous comparisons between these models for polaron hopping in edge-sharing Fe octahedral chains (Bylaska and Rosso 2018). Despite these differences, the electronic coupling, V AB , appears to be consistently (inversely) correlated to the distance between Fe atoms involved in the ET. The cluster model predicts a larger distance (3.21 Å) and smaller V AB than the periodic model. The distance seen in the periodic model (3.11 Å) is closer to that found in the original crystal structure, which suggests that the underlying lattice constrains the amount of distortion coming from the electron hole, a strength of the periodic planewave approach over use of small cluster models. However, another possibility is that the size of the supercell, which contains a 4 × 3 2D-grid of Fe octahedra, needs to be increased, since using a larger supercell will loosen up the structure to allow for larger distortions.
Finally, these results further support the use of plane-wave electronic structure methods containing accurate implementations of periodic exchange for modeling charge transfer process in solids. With large scale high-performance computers these methods are capable of performing ab initio molecular dynamics efficiently enough to carry out standard rare event simulations (e.g., potential of mean force, metadynamics, TAMD, or Markov state modeling with ab initio molecular dynamics) with ground-state adiabatic surfaces to find adiabatic ET rates, while at the same time being accurate enough to model ET rates and transmission factors, κ el

Conclusion
In summary, we have developed a new method for accurately calculating exact change integrals in periodic systems that is based on a Filon-like integration approach in which an effective screened potential is defined using a narrow integration that accounts for the singularity in the two-electron Coulomb integrand. The technique developed can readily be employed in plane-wave DFT programs, and it is applicable to both confined and extended systems, as well as AIMD simulations and periodic electron transfer calculations.
This new development has been demonstrated on several non-trivial applications, including hybrid-DFT calculations for the chlorination of ethene, adsorption of H 2 onto an Au − 20 nanoparticle, and electron hole transfer calculations in an Fe(II) containing mica, annite. These results show it to be accurate, and computationally this technique for evaluating periodic exact exchange has the same cost and at least the same amount of accuracy for large point calculations as our previous implementation based on based on maximally localized orbitals (Bylaska et al. 2011a). In addition, this new method has more flexibility in that several other approximate filtered exchange potentials can be derived with varying degrees of accuracy, including ones that can be used with Monkhorst-Pack Brillouin zone sampling. However, a current drawback of the method is that it requires the I 2 and I 3 integrals in Eqs.17-18 to be calculated with high accuracy. Future work will focus on improving the numerical convergence of these integrals with the singular integrands, as well as generalizing the method to evaluate the two-electron integrals between two determinant states (used by MCSCF and CASSCF methods), and evaluate two-electron integrals containing 4 molecular orbitals (used by many-body methods based on second quantization such as CI, CCSD, etc.). We are optimistic that this new development activity will be able to overcome the limitations of screened Coulomb interactions, and other highly-engineered approaches such as DFT+U (Dudarev et al. 1998;Anisimov et al. 1991), and make it easier to extend accurate many-body chemistry calculations to use periodic boundary conditions in the near future.

Appendix A: Generalized norm conserving pseudopotentials
The valence electron interactions with the atomic core are approximated with a generalized norm-conserving (Hamann 1989) or Troullier and Martins (1991) pseudopotentials modified to a separable form suggested by Kleinman and Bylander (1982). The input decks used to generate the pseudopotentials in this study can be obtained from the following directory on Github, https://github.com/nwchemgit/nwchem/tree/master/src/ nwpw/libraryps. Hamann pseudopotentials were used for hydrogen, carbon, oxygen, and silicon, and Troullier-Martins pseudopotentials were used for potassium, gold, and iron. The pseudopotentials were constructed using the following core radii, H: r cs =0.8 a.u and r cp =0.8 a.u.; C: r cs =0.8 a.u, r cp =0.85 a.u, and r cd =0.85 a.u.; O: r cs =0.7 a.u, r cp =0.7 a.u and r cd =0.7 a.u; Si: r cs =1.059 a.u, r cp =1.286 a.u and r cd =1.286 a.u.; K: r cs =4.2 a.u, r cp =4.05 a.u and r cd =3.75 a.u., r semicore =2.0 a.u.; Au: r cs =2.0761953 a.u, r cp =3.0973349 a.u, r cd =1.9261744 a.u. and r semicore = 1.2 a.u.; and Fe: r cs =1.8 a.u, r cp =1.8 a.u, r cd =1.8 a.u. For Fe, the 3s and 3p states are highly polarizable and for this element we generated the pseudo potential using a [Ne] core, where the 3s and 3p orbitals are placed in the valence space. These potentials were generated using the PBE96 exchange correlation functional.
In carrying out plane-wave DFT calculations, the stiffness of the atom center electron potential, Z I |r−R I | , needs be reduced for Kohn-Sham wave functions, ψ i (r), to be expanded in a plane-wave basis set. Pseudopotentials are one approach to do this that are widely used (Pickett 1989). An important feature of pseudopotentials is that they can be defined to varying degrees of accuracy, and it is typical to modify them in order that they have the required accuracy for a particular calculation.
A significant advance in the development of the pseudopotential method was made by Hamann et al. (1979) with the introduction of norm-conserving pseudopotentials. While there are differences between various pseudopotential approaches, all popular pseudopotentials adopt the basic prescription outlined the original work. These methods have been highly developed in the condensed matter community and are well explained and reviewed (Pickett 1989). The basic idea of pseudopotential is that the core region of the atomic potential is replaced by a much slower varying function designed to specifically reproduce the behavior of the valence wave functions in regions outside the core (presumed to be the bonding region). The smoothed potential has a node-less solution that can be expanded by a smaller plane-wave basis. It can be shown that with proper care, replacing the atomic potential with a pseudopotential will produce the same solutions beyond the region of replacement, while also maintaining the normalization of the orbital function.
In the norm-conserving pseudopotential, a pseudopotential for each total angular momentum l is found from the direct inversion of the Schrödinger equation (with a selected DFT functional) for an atom (Hamann 1989). This produces a non-local pseudopotential of the form, where V valence M (r) is the Coulomb and exchange potential due to the (non active) valence electrons, is the spherical harmonic defined by the angular momentum, l, and magnetic quantum, m, numbers,r is a unit vector in the r direction, and V l (r) is the radial potential found from the inversion of the DFT solution to the radial Schrödinger equation for the equivalent atomic problem (see reference Hamann et al. (1979)). The non-local pseudopoential operator V non−local acts on a function ψ i (r) by, This operator has a semi-local form, neither just local (radial) or fully separable (see KB [Bylander, 1984 ]). In this semi-local form, the pseudopotential is computationally difficult to calculate with a plane-wave basis set, since the kernel integration is not separable in r and r (see reference Kleinman and Bylander (1982)). To produce a more efficient calculation while retaining as much of the atomic form as possible Kleinman and Bylander approximated the form by, where the atom-centered projectors P lm (r) are of the form whereφ(r) are the zero-radial node pseudo-wavefunctions of the potentials, V l (r), calculated in the atomic environment. Note that i.e., that the fully non-local KB form preserves the form of the potential in the atomic problem. The choice of the local potential V local (r) is somewhat arbitrary, but for transition metals it is often chosen to be the V l=0 potential. A larger series expansion in pseudowavefunctions can be used to improve the fully local description of the semilocal form. This leads to the general form V non−local = V local (r) + l,m n,n P lm (r)h n,n l P * l,m (r ) Pseudopotentials are developed entirely from fitting atomic calculations and, therefore, should not be considered as part of the data fitting process. Nevertheless, there are questions about accuracy of the representation, and in particular how manyφ nl (r) are required to accurately represent the valence structure of the condensed system and how much of the unscreened atomic potential is assigned as the core region (roughly speaking the region removed). We have found for many 1st row transition metal systems that including the 3s and 3p functions in the active space of the pseudopotential considerably improved the agreement with the scattering data. The default pseudopotential included only the 3d orbitals. Additional issues that have to be considered in the pseudopotential representation include the functional form used for V l (r) and its parameterization that is obtained by comparison to atomic calculations.
A key parameter in these fittings is the radius of the replacement region (aka core region). The parameter to a large extent the smoothness of the pseudo potential (the larger the region the smoother the pseudopotential). However, if this region is too large the bond