- Original Article
- Open Access
- Published:

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

*Materials Theory*
**volume 4**, Article number: 3 (2020)

## Abstract

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 H_{2} to phenol and Au\(_{20}^{-}\) 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}^{\sigma }\) are the number of occupied spin orbitals. The notation for the volume with the bar above it, e.g. \(\bar {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} = \frac {8\pi ^{3}}{\Omega }\) 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 \(\sum _{\mathbf {a}} \frac {\exp \left (i \mathbf {K} \cdot \left (\mathbf {R} - \mathbf {a}\right)\right)}{|\mathbf {R} - \mathbf {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

where \(V_{BZ}=\frac {(2\pi)^{3}}{\Omega }\) 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, 1989) 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 \(\frac {4\pi }{|\mathbf {G}-\mathbf {k}+\mathbf {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 transformations is *J*=*b*_{1}·(*b*_{2}×*b*_{3})=*V*_{BZ}. Using this linear transformation, the above six-dimensional integrals can be recast as

where the kernel functions are defined as

with

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|=\frac {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

with

where *P*_{l,|m|}(cos*θ*) is an associated Legendre polynomial (with (−1)^{m} factor omitted) and \(\hat {\mathbf {x}} = (\cos \phi \sin \theta, \sin \phi \sin \theta, \cos \theta)\). 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 \(\frac {4\pi }{\mathbf {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|=\frac {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 such that *f*(|*γ*_{i}|≥1)=0, i.e.

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 6 face-sharing integrals are

the 12 edge-sharing integrals are

and the 8 corner-sharing integrals are

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 \(\frac {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*_{f1} 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 \(\frac {4\pi }{|\mathbf {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*_{f2} in Figs. 5 and 6 (dashed red line)), is shown to produce nearly the same results as *V*_{f1}. 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_{f3} = \frac {4\pi }{|\mathbf {G}|^{2}}\) (**G**≠0) (dashed green line), as well as slight modifications to it, *V*_{f4} (dashed blue line), and *V*_{f5} (dashed indigo line), are shown in Figs. 5 and 6. The *V*_{f4} and *V*_{f5} potentials were modified by replacing the **G** vectors in the first Brillouin zone and **G**=0, respectively by *I*_{2}+*I*_{3}. Note *V*_{f4} is the same as *V*_{f2} with *n*_{max}=0. As can be seen from the plots the traditional periodic form, *V*_{f3}, has a large error and produces a kernel in real-space that is significantly below *V*_{f1}. The *V*_{f4} and *V*_{f5} 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*_{f6} 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 wavefunctions, i.e. systems with large band gaps. Also the *V*_{f6} potential appears to be close to the *V*_{f4} and *V*_{f5} 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) and hybrid PBE0 (Adamo and Barone 1999) exchange correlation potentials were used. For the electron transfer calculations, open-shell spin-unrestricted 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 (Hermes et al. 2019; Hermes 2019) with a force convergence criterion of 10^{−4}eV/Å. For each of these transition states, the reaction pathways were then optimized using the intrinsic reaction coordinate (IRC) with a force termination criterion of 0.05eV/Å. In addition, we also performed Gaussian basis set calculations using the def2-TZVP basis set (Weigend and Ahlrichs 2005) for comparison. These results, shown in Fig. 8, not only produced 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\(_{20}^{-1}\) 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 chemi-adsorbed state, whereas the DFT method (PBE) only predicts an endothermic chemi-adsorbed 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 et al. (2003); Rosso and Ilton (2003); Rosso and Ilton (2005); Smith et al. (2006); Windus et al. (2003); Skomurski et al. (2011); Skomurski et al. (2010); Kerisit et al. (2007); 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 (\(\hbar \nu _{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 (*Δ**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 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 plane-wave 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 3*s* and 3*p* states are highly polarizable and for this element we generated the pseudo potential using a [Ne] core, where the 3*s* and 3*p* 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,\(\frac {Z_{I}}{|\mathbf {r}-\mathbf {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_{M}^{valence}(\mathbf {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, \(\hat {\mathbf {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

and the coefficient *h*_{l} is

where \(\tilde {\varphi }(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

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 \(\tilde {\phi }_{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 formation will be effected and the pseudopotential representation will produce incorrect bonding results. An example of the derived smooth pseudopotential and nodeless pseudo wavefunctions is given Fig. 12 for the Fe ^{3+} ion.

## Availability of data and materials

The codes developed are available in the NWChem software(http://www.nwchem-sw.org/). Input and output decks used to test the codes are available from the corresponding author.

## References

C. Adamo, V. Barone, Toward reliable density functional methods without adjustable parameters: The PBE0 model. J. Chem. Phys.

**110**(13), 6158–70 (1999).N. Adelstein, J. B. Neaton, M. Asta, L. C. De Jonghe, Density functional theory based calculation of small-polaron mobility in hematite. Phys. Rev. B.

**89**(24), 245115 (2014).G. Adragna, R. Del Sole, A. Marini, Ab initio calculation of the exchange-correlation kernel in extended systems. Phys. Rev. B.

**68**(16), 165108 (2003).V. I. Anisimov, J. Zaanen, O. K. Andersen, Band theory and mott insulators: Hubbard u instead of stoner i. Phys. Rev. B.

**44**(3), 943 (1991).F. Aryasetiawan, O. Gunnarsson, The gw method. Rep. Prog. Phys.

**61**(3), 237 (1998).R. Atta-Fynn, E. J. Bylaska, G. K. Schenter, W. A. De Jong, Hydration shell structure and dynamics of curium (iii) in aqueous solution: first principles and empirical studies. J. Phys. Chem. A.

**115**(18), 4665–4677 (2011).J. Andzelm, J. Baker, A. Scheiner, M. Wrinn, A density functional study of chemical reactions. Int. J. Quantum Chem.

**56**(6), 733–746 (1995).A. Barducci, G. Bussi, M. Parrinello, Well-tempered metadynamics: a smoothly converging and tunable free-energy method. Phys. Rev. Lett.

**100**(2), 020603 (2008).P. K. Behara, M. Dupuis, Electron transfer in extended systems: characterization by periodic density functional theory including the electronic coupling. Phys. Chem. Chem. Phys. (2019). https://doi.org/10.1039/C9CP05133C.

J. Behler, B. Delley, K. Reuter, M. Scheffler, Nonadiabatic potential-energy surfaces by constrained density-functional theory. Phys. Rev. B.

**75**(11), 115409 (2007).A. Bergner, M. Dolg, W. Küchle, H. Stoll, H. Preuß, Ab initio energy-adjusted pseudopotentials for elements of groups 13–17. Mol. Phys.

**80**(6), 1431–1441 (1993).G. R. Bowman, K. A. Beauchamp, G. Boxer, V. S. Pande, Progress and challenges in the automated construction of markov state models for full protein systems. J. Chem. Phys.

**131**(12), 124101 (2009).G. Bussi, A. Laio, M. Parrinello, Equilibrium free energies from nonequilibrium metadynamics. Phys. Rev. Lett.

**96**(9), 090601 (2006).E. J. Bylaska, P. R. Taylor, R. Kawai, J. H. Weare, Lda predictions of c20 isomerizations: Neutral and charged species. J. Phys. Chem.

**100**(17), 6966–6972 (1996).E. J. Bylaska, M. Valiev, R. Kawai, J. H. Weare, Parallel implementation of the projector augmented plane wave method for charged systems. Comput. Phys. Commun.

**143**(1), 11–28 (2002).E. J. Bylaska, K. Tsemekhman, S. B. Baden, J. H. Weare, H. Jonsson, Parallel implementation of

*γ*-point pseudopotential plane-wave dft with exact exchange. J. Comput. Chem.**32**(1), 54–69 (2011a).E. Bylaska, K. Tsemekhman, N. Govind, M. Valiev, in

*Computational Methods for Large Systems*. Large-scale plane-wave-based density functional theory: formalism, parallelization, and applications (Wiley, 2011b), pp. 77–116. https://doi.org/10.1002/9780470930779.ch3.E. J. Bylaska, in

*Paper presented at The 27th Annual Workshop on Recent Developments in Electronic Structure Theory, University of Washington, Seattle, WA*. Improving the performance of ab initio molecular dynamics simulations and band structure calculations for actinide and geochemical systems with new algorithms and new machines, (2015). http://monalisa.phys.washington.edu/ES2015/InvitedAbstract/BylaskaEric.pdf. Accessed 1 Nov 2019.E. J. Bylaska, in

*Annual Reports in Computational Chemistry*, 13. Plane-wave dft methods for chemistry (Elsevier, 2017), pp. 185–228. https://doi.org/10.1016/bs.arcc.2017.06.006.E. J. Bylaska, K. Rosso, Corresponding orbitals derived from periodic bloch states for electron transfer calculations of transition metal oxides. J. Chem. Theory Comput.

**14**(8), 4416–4426 (2018).R. Car, M. Parrinello, Unified approach for molecular dynamics and density-functional theory. Phys. Rev. Lett.

**55**(22), 2471 (1985).S. Chawla, G. A. Voth, Exact exchange in ab initio molecular dynamics: An efficient plane-wave based algorithm. J. Chem. Phys.

**108**(12), 4697–4700 (1998).Y. Chen, E. Bylaska, J. Weare, in

*Molecular Modeling of Geochemical Reactions: An Introduction*. First principles estimation of geochemically important transition metal oxide properties (Wiley, 2016), pp. 107–149. https://doi.org/10.1002/9781118845226.ch4.Y. Chen, E. J. Bylaska, J. H. Weare, Weakly bound water structure, bond valence saturation and water dynamics at the goethite (100) surface/aqueous interface: ab initio dynamical simulations. Geochem. Trans.

**18**(1), 3 (2017).S. Dudarev, G. Botton, S. Savrasov, C. Humphreys, A. Sutton, Electron-energy-loss spectra and the structural stability of nickel oxide: An lsda+ u study. Phys. Rev. B.

**57**(3), 1505 (1998).R. Dogonadze, A. Kuznetsov, M. Vorotyntsev, On the theory of adiabatic and non-adiabatic electrochemical reactions. J. Electroanal. Chem. Interfacial Electrochem.

**25**(2), 17–19 (1970).A. Farazdel, M. Dupuis, E. Clementi, A. Aviram, Electric-field induced intramolecular electron transfer in spiro. pi.-electron systems and their suitability as molecular electronic devices. a theoretical study. J. Am. Chem. Soc.

**112**(11), 4206–4214 (1990).L. N. G. Filon, Iii.—on a quadrature formula for trigonometric integrals. Proc. R. Soc. Edinb.

**49:**, 38–47 (1930).D. Frenkel, B. Smit,

*Understanding Molecular Simulation: from Algorithms to Applications, vol. 1*(Academic press, San Diego, 2001).K. Fukui, The path of chemical reactions - the IRC approach. Acc. Chem. Res.

**14**(12), 363–368 (1981). https://doi.org/10.1021/ar00072a001.J. L. Fulton, E. J. Bylaska, S. Bogatko, M. Balasubramanian, E. Cauët, G. K. Schenter, J. H. Weare, Near-quantitative agreement of model-free dft-md predictions with xafs observations of the hydration structure of highly charged transition-metal ions. J. Phys. Chem. Lett.

**3**(18), 2588–2593 (2012).M. Garcia-Viloca, J. Gao, M. Karplus, D. G. Truhlar, How enzymes work: analysis by modern rate theory and computer simulations. Science.

**303**(5655), 186–195 (2004).A. Görling, Exact treatment of exchange in kohn-sham band-structure schemes. Phys. Rev. B.

**53**(11), 7024 (1996).C. Gu, H. -W. Chang, L. Maibaum, V. S. Pande, G. E. Carlsson, L. J. Guibas, Building markov state models with solvent dynamics. BMC Bioinformatics.

**14**(2), 8 (2013).F. Gygi, A. Baldereschi, in

*Helvetica Physica Acta*, 58. Exact exchange calculations of electronic-properties of solids, (1985), pp. 928–928. BIRKHAUSER VERLAG AG PO BOX 133 KLOSTERBERG 23, CH-4010 BASEL, SWITZERLAND.F. Gygi, A. Baldereschi, Self-consistent hartree-fock and screened-exchange calculations in solids: Application to silicon. Phys. Rev. B.

**34**(6), 4405 (1986).F. Gygi, A. Baldereschi, Quasiparticle energies in semiconductors: Self-energy correction to the local-density approximation. Phys. Rev. Lett.

**62**(18), 2160 (1989).D. Hamann, M. Schlüter, C. Chiang, Norm-conserving pseudopotentials. Phys. Rev. Lett.

**43**(20), 1494 (1979).D. R. Hamann, Generalized norm-conserving pseudopotentials. Phys. Rev. B.

**40**(5), 2980 (1989).J. Harl, G. Kresse, Cohesive energy curves for noble gas solids calculated by adiabatic connection fluctuation-dissipation theory. Phys. Rev. B.

**77**(4), 045136 (2008).J. Harl, G. Kresse, Accurate bulk properties from approximate many-body techniques. Phys. Rev. Lett.

**103**(5), 056401 (2009).J. Harl, L. Schimka, G. Kresse, Assessing the quality of the random phase approximation for lattice constants and atomization energies of solids. Phys. Rev. B.

**81**(11), 115126 (2010).P. J. Hay, W. R. Wadt, Ab initio effective core potentials for molecular calculations. potentials for k to au including the outermost core orbitals. J. Chem. Phys.

**82**(1), 299–310 (1985).L. Hedin, New method for calculating the one-particle green’s function with application to the electron-gas problem. Phys. Rev.

**139**(3A), 796 (1965).E. D. Hermes, K. Sargsyan, H. N. Najm, J. Zádor, Accelerated Saddle Point Refinement through Full Exploitation of Partial Hessian Diagonalization. J. Chem. Theory. Comput. (2019). https://doi.org/10.1021/acs.jctc.9b00869.

E. D. Hermes, Sella (2019). https://doi.org/10.5281/zenodo.3379094.

J. Heyd, G. E. Scuseria, M. Ernzerhof, Hybrid functionals based on a screened coulomb potential. J. Chem. Phys.

**118**(18), 8207–8215 (2003).J. Heyd, G. E. Scuseria, Efficient hybrid density functional calculations in solids: Assessment of the heyd–scuseria–ernzerhof screened coulomb hybrid functional. J. Chem. Phys.

**121**(3), 1187–1192 (2004).S. Hirata, M. R. Hermes, J. Simons, J. Ortiz, General-order many-body green’s function method. J. Chem. Theory Comput.

**11**(4), 1595–1606 (2015).T. Holstein, Studies of polaron motion: Part i. the molecular-crystal model. Ann. Phys.

**8**(3), 325–342 (1959).T. Holstein, Studies of polaron motion: Part ii. the "small" polaron. Ann. Phys.

**8**(3), 343–389 (1959).N. Iordanova, M. Dupuis, K. M. Rosso, Charge transport in metal oxides: a theoretical study of hematite

*α*-fe 2 o 3. J. Chem. Phys.**122**(14), 144305 (2005).N. Iordanova, M. Dupuis, K. M. Rosso, Theoretical characterization of charge transport in chromia (

*α*-cr 2 o 3). J. Chem. Phys.**123**(7), 074710 (2005).A. W. Jasper, C. Zhu, S. Nangia, D. G. Truhlar, Introductory lecture: Nonadiabatic effects in chemical dynamics. Faraday Discuss.

**127:**, 1–22 (2004).R. A. Kendall, E. Aprà, D. E. Bernholdt, E. J. Bylaska, M. Dupuis, G. I. Fann, R. J. Harrison, J. Ju, J. A. Nichols, J. Nieplocha, T. P. Straatsma, T. L. Windus, A. T. Wong, High performance computational chemistry: An overview of NWChem a distributed parallel application. Comp. Phys. Commun.

**128**(1), 260–283 (2000). https://doi.org/10.1016/S0010-4655(00)00065-5.S. Kerisit, K. M. Rosso, Charge transfer in feo: A combined molecular-dynamics and ab initio study. J. Chem. Phys.

**123**(22), 224712 (2005).S. Kerisit, K. M. Rosso, Computer simulation of electron transfer at hematite surfaces. Geochim. Cosmochim. Acta.

**70**(8), 1888–1903 (2006).S. Kerisit, K. M. Rosso, M. Dupuis, M. Valiev, Molecular computational investigation of electron-transfer kinetics across cytochrome- iron oxide interfaces. J. Phys. Chem. C.

**111**(30), 11363–11375 (2007).H. F. King, R. E. Stanton, H. Kim, R. E. Wyatt, R. G. Parr, Corresponding orbitals and the nonorthogonality problem in molecular quantum mechanics. J. Chem. Phys.

**47**(6), 1936–1941 (1967).L. Kleinman, D. Bylander, Efficacious form for model pseudopotentials. Phys. Rev. Lett.

**48**(20), 1425 (1982).W. Kohn, L. J. Sham, Self-consistent equations including exchange and correlation effects. Phys. Rev.

**140**(4A), 1133 (1965).K. Kowalski, K. Bhaskaran-Nair, W. Shelton, Coupled-cluster representation of green function employing modified spectral resolutions of similarity transformed hamiltonians. J. Chem. Phys.

**141**(9), 094102 (2014).Y. Kurosaki, Ab initio molecular orbital study of the c2h4+ cl2 → c2h4cl2 reaction. J. Mol. Struct. THEOCHEM.

**503**(3), 231–240 (2000).A. M. Kuznetsov, J. Ulstrup,

*Electron Transfer in Chemistry and Biology: an Introduction to the Theory*(Wiley, Hoboken, 1999).A. Laio, F. L. Gervasio, Metadynamics: a method to simulate rare events and reconstruct the free energy in biophysics, chemistry and material science. Rep. Prog. Phys.

**71**(12), 126601 (2008).D. C. Langreth, J. P. Perdew, Exchange-correlation energy of a metallic surface: Wave-vector analysis. Phys. Rev. B.

**15**(6), 2884 (1977).V. Levich, R. -R. Dogonadze, E. German, A. Kuznetsov, Y. I. Kharkats, Theory of homogeneous reactions involving proton transfer. Electrochimica Acta.

**15**(2), 353–367 (1970).G. I. Likhtenshtein,

*Solar Energy Conversion: Chemical Aspects*(Wiley, Hoboken, 2012).J. Malinsky, Y. Magarshak, A green’s function approach to long-range electron transfer in macromolecules. J. Phys. Chem.

**96**(7), 2849–2851 (1992).L. Maragliano, E. Vanden-Eijnden, Single-sweep methods for free energy calculations. J. Chem. Phys.

**128**(18), 184110 (2008).R. A. Marcus, Chemical and electrochemical electron-transfer theory. Ann. Rev. Phys. Chem.

**15**(1), 155–196 (1964).R. A. Marcus, Electron transfer reactions in chemistry. theory and experiment. Rev. Modern Phys.

**65**(3), 599 (1993).R. A. Marcus, Electron transfer reactions in chemistry: theory and experiment (nobel lecture). Angew. Chem. Int. Ed.

**32**(8), 1111–1121 (1993).M. Marsman, J. Paier, A. Stroppa, G. Kresse, Hybrid functionals applied to extended systems. J. Phys. Condens. Matter.

**20**(6), 064201 (2008).D. Marx, J. Hutter,

*Ab Initio Molecular Dynamics: Basic Theory and Advanced Methods*(Cambridge University Press, 2009). https://doi.org/10.1017/cbo9780511609633.N. Marzari, D. Vanderbilt, Maximally localized generalized wannier functions for composite energy bands. Phys. Rev. B.

**56**(20), 12847 (1997).R. Morrison, R. Boyd (1992).

P. Nachtigall, K. Jordan, A. Smith, H. Jónsson, Investigation of the reliability of density functional methods: Reaction and activation energies for si–si bond cleavage and h2 elimination from silanes. J. Chem. Phys.

**104**(1), 148–58 (1996).M. D. Newton, N. Sutin, Electron transfer reactions in condensed phases. Ann. Rev. Phys. Chem.

**35**(1), 437–480 (1984).M. D. Newton, H. L. Friedman, Green function theory of charge transfer processes in solution. J. Chem. Phys.

**88**(7), 4460–4472 (1988).Y. Niquet, M. Fuchs, X. Gonze, Exchange-correlation potentials in the adiabatic connection fluctuation-dissipation framework. Phys. Rev. A.

**68**(3), 032507 (2003).V. S. Pande, K. Beauchamp, G. R. Bowman, Everything you wanted to know about markov state models but were afraid to ask. Methods.

**52**(1), 99–105 (2010).M. C. Payne, M. P. Teter, D. C. Allan, T. A. Arias, JD Joannopoulos, Iterative minimization techniques for ab initio total-energy calculations: molecular dynamics and conjugate gradients. Rev Mod Phys.

**64:**, 1045–97 (1992). https://doi.org/10.1103/RevModPhys.64.1045.J. P. Perdew, K. Burke, M. Ernzerhof, Generalized gradient approximation made simple. Phys. Rev. Lett.

**77**(18), 3865 (1996).R. Peverati, D. G. Truhlar, Screened-exchange density functionals with broad accuracy for chemistry and solid-state physics. Phys. Chem. Chem. Phys.

**14**(47), 16187–16191 (2012).W. E. Pickett, Pseudopotential methods in condensed matter applications. Comput. Phys. Rep.

**9**(3), 115–197 (1989).E. D. Rainville,

*Special Functions*(Macmillan, New York, 1960).J. Reimers, N. Hush, Electronic properties of transition-metal complexes determined from electroabsorption (stark) spectroscopy. 2. mononuclear complexes of ruthenium (ii). J. Phys. Chem.

**95**(24), 9773–9781 (1991).J. R. Reimers, L. K. McKemmish, R. H. McKenzie, N. S. Hush, Non-adiabatic effects in thermochemistry, spectroscopy and kinetics: the general importance of all three born–oppenheimer breakdown corrections. Phys. Chem. Chem. Phys.

**17**(38), 24641–24665 (2015).D. K. Remler, P. A. Madden, Molecular dynamics without effective potentials via the car-parrinello approach. Mol. Phys.

**70**(6), 921–66 (1990).H. Ren, M. R. Provorse, P. Bao, Z. Qu, J. Gao, Multistate density functional theory for effective diabatic electronic coupling. J. Phys. Chem. Lett.

**7**(12), 2286–2293 (2016).L. Reining, V. Olevano, A. Rubio, G. Onida, Excitonic effects in solids described by time-dependent density-functional theory. Phys. Rev. Lett.

**88**(6), 066404 (2002).K. M. Rosso, E. S. Ilton, Charge transport in micas: The kinetics of fe ii/iii electron transfer in the octahedral sheet. J. Chem. Phys.

**119**(17), 9207–9218 (2003).K. M. Rosso, D. M. Smith, M. Dupuis, Aspects of aqueous iron and manganese (ii/iii) self-exchange electron transfer reactions. J. Phys. Chem. A.

**108**(24), 5242–5248 (2004).K. M. Rosso, D. M. Smith, Z. Wang, C. C. Ainsworth, J. K. Fredrickson, Self-exchange electron transfer kinetics and reduction potentials for anthraquinone disulfonate. J. Phys. Chem. A.

**108**(16), 3292–3303 (2004).K. M. Rosso, D. M. Smith, M. Dupuis, An ab initio model of electron transport in hematite (

*α*-fe 2 o 3) basal planes. J. Chem. Phys.**118**(14), 6455–6466 (2003).K. M. Rosso, E. S. Ilton, Effects of compositional defects on small polaron hopping in micas. J. Chem. Phys.

**122**(24), 244709 (2005).M. Sarich, R. Banisch, C. Hartmann, C. Schütte, Markov state models for rare events in molecular dynamics. Entropy.

**16**(1), 258–286 (2013).S. S. Shaik, P. C. Hiberty,

*A Chemist’s Guide to Valence Bond Theory*(Wiley, Hoboken, 2007).M. Shishkin, M. Marsman, G. Kresse, Accurate quasiparticle spectra from self-consistent gw calculations with vertex corrections. Phys. Rev. Lett.

**99**(24), 246403 (2007).D. M. Sherman, The electronic structures of fe3+ coordination sites in iron oxides: Applications to spectra, bonding, and magnetism. Phys. Chem. Miner.

**12**(3), 161–175 (1985).D. M. Sherman, Cluster molecular orbital description of the electronic structures of mixed-valence iron oxides and silicates. Solid State Commun.

**58**(10), 719–723 (1986).D. M. Sherman, Molecular orbital (scf-x

*α*-sw) theory of metal-metal charge transfer processes in minerals. Phys. Chem. Miner.**14**(4), 355–363 (1987).P. L. Silvestrelli, Maximally localized wannier functions for simulations with supercells of general symmetry. Phys. Rev. B.

**59**(15), 9703 (1999).C. Stampfl, W. Mannstadt, R. Asahi, A. J. Freeman, Electronic structure and physical properties of early transition metal mononitrides: Density-functional theory lda, gga, and screened-exchange lda flapw calculations. Phys. Rev. B.

**63**(15), 155106 (2001).F. N. Skomurski, S. Kerisit, K. M. Rosso, Structure, charge distribution, and electron hopping dynamics in magnetite (fe 3 o 4)(100) surfaces from first principles. Geochim. Cosmochim. Acta.

**74**(15), 4234–4248 (2010).F. N. Skomurski, E. S. Ilton, M. H. Engelhard, B. W. Arey, K. M. Rosso, Heterogeneous reduction of u 6+ by structural fe 2+ from theory and experiment. Geochim. Cosmochim. Acta.

**75**(22), 7277–7290 (2011).D. M. Smith, K. M. Rosso, M. Dupuis, M. Valiev, T. Straatsma, Electronic coupling between heme electron-transfer centers and its decay with distance depends strongly on relative orientation. J. Phys. Chem. B.

**110**(31), 15582–15588 (2006).A. Sorouri, W. M. C. Foulkes, N. D. Hine, Accurate and efficient method for the treatment of exchange in a plane-wave basis. J. Chem. Phys.

**124**(6), 064105 (2006).A. G. Stack, K. M. Rosso, D. M. Smith, C. M. Eggleston, Reaction of hydroquinone with hematite: Ii. calculated electron-transfer rates and comparison to the reductive dissolution rate. J. Colloid Interf. Sci.

**274**(2), 442–450 (2004).L. Song, J. Gao, On the construction of diabatic and adiabatic potential energy surfaces based on ab initio valence bond theory. J. Phys. Chem. A.

**112**(50), 12925–12935 (2008).J. E. Subotnik, S. Yeganeh, R. J. Cave, M. A. Ratner, Constructing diabatic states from adiabatic states: Extending generalized mulliken–hush to multiple charge centers with boys localization. J. Chem. Phys.

**129**(24), 244101 (2008).D. Ter Haar,

*Collected Papers of LD Landau*(Elsevier, Amsterdam, 2013).N. Troullier, J. L. Martins, Efficient pseudopotentials for plane-wave calculations. Phys. Rev. B.

**43**(3), 1993 (1991).M. Valiev, E. J. Bylaska, N. Govind, K. Kowalski, T. P. Straatsma, H. J. Van Dam, D. Wang, J. Nieplocha, E. Aprà, T. L. Windus, et al., NWChem: a comprehensive and scalable open-source solution for large scale molecular simulations. Comput. Phys. Commun.

**181**(9), 1477–89 (2010).E. Vanden-Eijnden, Some recent techniques for free energy calculations. J. Comput. Chem.

**30**(11), 1737–1747 (2009).M. van Schilfgaarde, T. Kotani, S. Faleev, Quasiparticle self-consistent g w theory. Phys. Rev. Lett.

**96**(22), 226402 (2006).T. Van Voorhis, T. Kowalczyk, B. Kaduk, L. -P. Wang, C. -L. Cheng, Q. Wu, The diabatic picture of electron transfer, reaction barriers, and molecular dynamics. Ann. Rev. Phys. Chem.

**61:**, 149–170 (2010).Z. Wang, C. Liu, X. Wang, M. J. Marshall, J. M. Zachara, K. M. Rosso, M. Dupuis, J. K. Fredrickson, S. Heald, L. Shi, Kinetics of reduction of fe (iii) complexes by outer membrane cytochromes mtrc and omca of shewanella oneidensis mr-1. Appl. Environ. Microbiol.

**74**(21), 6746–6755 (2008).F. Weigend, R. Ahlrichs, Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for h to rn: Design and assessment of accuracy. Phys. Chem. Chem. Phys.

**7**(18), 3297–3305 (2005).T. L. Windus, E. J. Bylaska, M. Dupuis, S. Hirata, L. Pollack, D. M. Smith, T. P. Straatsma, E. Aprà, in

*International Conference on Computational Science*. Nwchem: new functionality (Springer, 2003), pp. 168–177.C. Wittig, The landau- zener formula. J. Phys. Chem. B.

**109**(17), 8428–8430 (2005).S. Yeganeh, M. A. Ratner, V. Mujica, Dynamics of charge transfer: Rate processes formulated with nonequilibrium Green’s functions. J. Chem. Phys.

**126**(16), 161103 (2007). https://doi.org/10.1063/1.2735606.C. Zener, in

*Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences*, 137. Non-adiabatic crossing of energy levels (The Royal Society, 1932), pp. 696–702.C. Zener. Dissociation of excited diatomic molecules by external perturbations, Proceedings of the Royal Society of London. Series A, vol. 140, (1933), pp. 660–668.

Y. Zhao, N. González-García, D. G. Truhlar, Benchmark database of barrier heights for heavy atom transfer, nucleophilic substitution, association, and unimolecular reactions and its use to test theoretical methods. J. Phys. Chem. A.

**109**(9), 2012–2018 (2005).

## Acknowledgements

We would like to thank Kiril Tsemekhman, Hannes Jonsson, Scott Baden, and John H. Weare for there contributions to previous work on exact exchange that inspired this manuscript.

## Funding

This work was supported by the U.S. Department of Energy (DOE), Office of Science, Office of Basic Energy Sciences, Chemical Sciences, Geosciences, and Biosciences Division (DOE BES Chem) through its Geosciences program at Pacific Northwest National Laboratory (PNNL), DE-AC06-76RLO 1830. We also would like to thank the DOE BES Chem CCS, DOE BES Chem QIS, and DOE Advanced Scientific Computing Research (ASCR) ECP NWChemEx programs for their support of software development for high-performance computers and computer time needed to carry out the work. This research used resources of the National Energy Research Scientific Computing Center (NERSC), a User Facility supported by the Office of Science of the U.S. DOE under Contract No. DE-AC02-05CH11231. We wish to thank for support from the NERSC NESAP program. A portion of this research was performed using EMSL, a U.S. DOE Office of Science User Facility sponsored by the Office of Biological and Environmental Research, DE-AC06-76RLO 1830.

## Author information

### Affiliations

### Contributions

All authors contributed to the writing of the article. EJB and KW developed the model and developed the codes. EJB, JZ, EH, and KMR contributed to application of the codes. JZ and KMR were the PIs of computing projects used in the development of the codes. All authors read and approved the final manuscript.

### Corresponding author

## Ethics declarations

### Competing interests

The authors declare that they have no competing interests.

## Additional information

### Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

## Rights and permissions

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

## About this article

### Cite this article

Bylaska, E.J., Waters, K., Hermes, E.D. *et al.* A Filon-like integration strategy for calculating exact exchange in periodic boundary conditions: a plane-wave DFT implementation.
*Mater Theory* **4, **3 (2020). https://doi.org/10.1186/s41313-020-00019-9

Received:

Accepted:

Published:

### Keywords

- NWChem
- High-performance chemistry
- Plane-wave DFT
- Pseudopotentials
- Projector augmented wave
- PAW
- PSPW
- Ab initio Molecular Dynamics
- AIMD
- Periodic exact exchange
- UHF
- DFT
- Electron transfer