Implementation of Constrained DFT for Computing Charge Transfer

Article pubs.acs.org/JCTC

Implementation of Constrained DFT for Computing Charge Transfer Rates within the Projector Augmented Wave Method Marko Melander,*,† Elvar Ö . Jónsson,‡ Jens J. Mortensen,¶ Tejs Vegge,† and Juan Maria García Lastra*,† †

Department of Energy Conversion and Storage, Technical University of Denmark, DK-4000 Roskilde, Denmark COMP, Applied Physics Department, Aalto University FI-00076 Aalto, Espoo, Finland ¶ Department of Physics, Center for Atomic-scale Materials Design, Technical University of Denmark, DK-2800 Kgs. Lyngby, Denmark ‡

ABSTRACT: Combining constrained density function theory (cDFT) with Marcus theory is an efficient and promising way to address charge transfer reactions. Here, we present a general and robust implementation of cDFT within the projector augmented wave (PAW) framework. PAW pseudopotentials offer a reliable frozen-core electron description across the whole periodic table, with good transferability, as well as facilitate the extraction of all-electron quantities. The present implementation is applicable to two different wave function representations, atomic-centered basis sets (LCAO) and the finite-difference (FD) approximation utilizing real-space grids. LCAO can be used for large systems, molecular dynamics, or quick initialization, while more accurate calculations are achieved with the FD basis. Furthermore, the calculations can be performed with flexible boundary conditions, ranging from isolated molecules to periodic systems in one-, two-, or three-dimensions. As such, this implementation is relevant for a wide variety of applications. We also present how to extract the electronic coupling element and reorganization energy from the resulting diabatic cDFT-PAW wave functions for the parametrization of Marcus theory. Here, the combined method is applied to important test cases where practical implementations of DFT fail due to the self-interaction error, such as the dissociation of the helium dimer cation, and it is compared to other established cDFT codes. Moreover, for charge localization in a diamine cation, where it was recently shown that the commonly used generalized gradient and hybrid functionals of DFT failed to produce the localized state, cDFT produces qualitatively and quantitatively accurate results when benchmarked against self-interaction corrected DFT and high-level CCSD(T) calculations at a fraction of the computational cost.

■

INTRODUCTION

Electron transfer is a crucial elementary step occurring in several biologically, chemically, and industrially import reactions. Albeit being conceptually the simplest of chemical reactions, charge transfer reactions are notoriously difficult to describe within density functional theory (DFT). This is due to the spurious selfinteraction error (SIE)1 inherent to practical implementations of the theory. In many cases, the charge transfer rate can be computed using the quantum mechanical version of Marcus theory2−5 2π kct = ℏ

⎡ (λ + ΔG0)2 ⎤ exp⎢ − ⎥ 4kBTλ ⎦ ⎣ 4πkBTλ

Figure 1. Schematic view of the Marcus theory. The green and red lines define the ground and excited adiabatic states, respectively, whereas the dotted lines present the initial and final diabatic states.

|Hab|2

(1)

where the electronic coupling element Hab, reorganization free energy λ, and reaction free energy ΔG0 define the charge transfer rate. The variables are shown schematically in Figure 1. The electronic coupling element and reorganization free energy need to be computed using diabatic or charge-localized states that are not eigenfunctions of the electronic Hamiltonian. The traditional way of obtaining accurate diabatic states is to diabatize the © XXXX American Chemical Society

adiabatic states from a high-level electronic structure calculation (see section 4.1.6 of ref 6 and references therein). Instead of using such computationally expensive methods, it has recently been shown that constructing charge-localized states Received: August 18, 2016 Published: October 17, 2016 A

DOI: 10.1021/acs.jctc.6b00815 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation using cDFT7−9 offers an accurate and efficient way of computing the electronic coupling elements and the reorganization free energy in many cases.10,11 Compared to conventional DFT, which suffers from unphysical electron delocalization due to the SIE, cDFT offers a way for constructing charge-localized states by introducing an additional external potential to the Kohn−Sham (KS) equations, resulting in an effective self-interaction correction (SIC). Thus, efficient underlying DFT machinery with small modifications can be used for creating the diabatic states needed for computing the parameters in Marcus theory. cDFT was first established in 1984 as an extension to DFT in order to study lowest-energy excited states under a specified charge or magnetization constraint imposed by an additional external potential.12 In the first cDFT implementations, the external potential satisfying the constraint was not solved selfconsistently. Modern, self-consistent cDFT was formulated in 2005,7 and since then, the method has been implemented in several codes. Most of them utilize localized basis sets (NWChem,7 SIESTA,13 Q-Chem,14 deMon2k,15 ADF16), and also plane wave implementations exist (CPMD,17,18 QuantumESPRESSO,19 VASP20 for noncollinear spin only). Linear-scaling cDFT implementations rely on either numerical basis functions (CONQUEST21) or wavelets (BigDFT22). NWChem and QChem implementations are intended for molecular calculations, deMon2k focuses on QM/MM calculations using LCAO and auxiliary density fitting, ADF combines cDFT with density functional embedding,23 and periodic calculations can be done with (CPMD, SIESTA, QuantumESPRESSO, and VASP) or without k-points (CONQUEST and BigDFT). If pseudopotentials are used, they are norm-conserving Goedecker−Teter− Hutter24 (GTH) or Troullier−Martins25 (TM) pseudopotentials as in CPMD and SIESTA, respectively. A projector augmented wave (PAW) implementation is included in VASP but only for the treatment of noncollinear spin calculations (i.e., it is not a general tool of the code that would allow one to study charge transfer reactions, for instance). Our current implementation differs from its predecessors in several ways. First, we utilize modern, non-normconserving PAW pseudopotentials for arbitrary charge and spin constraints. The PAW method offers increased accuracy over the GTH and TM pseudopotentials,26,27 can cover the entire periodic table28 with good transferability,29 and offers the possibility of extracting (frozen core) all-electron quantities. Second, in GPAW, the cDFT-PAW wave functions can be presented numerically on a grid to facilitate accurate and basis-set-free calculations or using a linear combination of atomic orbitals (LCAO) for fast calculations such as initial geometry optimization or ab initio molecular dynamics, for example. Third, the periodicity can be chosen flexibly between isolated molecules and periodic systems in one-, two-, or three-dimensions with k-point sampling. Fourth, we have also implemented the possibility of spin-constrained regions, which is missing in other solid-state cDFT implementations. All in all, the combination of these features offers the possibility of treating molecules, surfaces, and the bulk on the same footing with the use of accurate and transferable pseudopotentials utilizing quick LCAO or more accurate gridbased calculations. cDFT has been verified as a useful tool for widening the scope of ground-state DFT to excitation processes, correcting for selfinteraction energy in current DFT functionals, excitation energy, and electron transfer as well as parametrizing model Hamiltonians, for example.9 Our long-term goal is to study charge transfer kinetics at surfaces, in the bulk, and in solution. At

surfaces, the interest is in modeling heterogeneous charge transfer, which is fundamental to all electrocatalysis. Thus, far, commonly used models for predicting and explaining electrocatalytic reactions are based on equilibrium thermodynamic reasoning,30 while the charge transfer kinetics at the solid−liquid interface are often neglected. In the bulk, our aim is to understand and predict charge conduction in semiconductors and insulators, where polaronic conduction via consecutive charge hopping is dominant. In solution, our interest is in understanding charge transfer reactions in oxygen evolution and reduction reactions taking place in metal−air batteries, among others. The current implementation offers a generally applicable, efficient, and accurate tool for computation of charge transfer kinetics. The article is arranged as follows. In the Theory section, cDFT is presented and our choice of weight function for the constraining potential discussed. A brief overview of the PAW method follows, including the constraining potential, as well as additional force terms due to the constraints in the PAW framework. Finally, the parametrization of Marcus theory using cDFT is presented. Next is the Computational Details section, followed by the Test Calculations section, where a variety of systems are calculated with cDFT and the results compared to other implementations thereof. A special case is reported in the Charge Localization in an N,N′-dimethylpiperazine Cation section, where the results are compared to self-interaction corrected DFT calculations and high-level quantum chemistry calculations. The implementation and results are summarized in the Conclusions section.

■

THEORY Constrained DFT. The constrained density functional method (cDFT) as presented here follows the formulation by Van Voorhis.7−9 The key modification to normal DFT is the introduction of an auxiliary potential to force a certain region (in real space and around a molecule, molecular fragment, or atom) to carry a predefined charge. The energy functional can then be written as a sum of the usual KS functional and a constraining term F[n(r), Vc] = EKS[n] +

∑ Vc ∑ (∫ dr wcσ(r)nσ (r) − Nc) c

σ

(2)

Here wσi (r) is the weight function that defines how the charge is to be partitioned, that is, the regions where charge is to be localized, and Ni is the desired number of electrons within the constrained region. This will be discussed in more detail below. A charge density difference between a donor (D) and an acceptor (A) is achieved by defining w = wA − wD, while the local magnetization can be constrained by working with spin densities and by setting wβ = −wα. All of the above constraints are included in eq 2. The introduction of constraining terms in eq 2 leads to a new effective potential defined as σ = veff

δF[n(r)] δEKS[n(r)] = + δn(r) δn(r)

∑ Vcwcσ(r) c

(3)

Thus, the cDFT potential is just the sum of the usual KS potential and the constraining potential, which is also used in the selfconsistent calculation. The constraint is further enforced by introducing the convergence criteria B


Article

Journal of Chemical Theory and Computation C ≥ |∑

∫ dr wcσ(r)nσ (r) − Nc|, ∀ c

σ

Projector Augmented Wave Method. In the PAW method,38 the rapidly oscillating KS wave functions {ψn(r)} are transformed to smooth pseudo wave functions {ψ̃ n(r)} near the nuclear regions. This also affects the electron density and how the cDFT energy is computed. We therefore review the basic ingredients of the PAW formalism needed for the present cDFT implementation, while more complete presentations can be found in refs 38−41. The transformation from the KS wave functions to the pseudo wave functions is formally carried out using an operator ;

(4)

From now on, we will only work with the unified notation, performing the spin summation according to the weight function rules above and omitting explicit spin notation, unless necessary. Weight Function. The weight function has the only criterion that wc(r) ∈ [0,1], which offers a lot of freedom for its construction. Several alternative definitions have been discussed in ref 9, but we chose to use Hirshfeld-type31 partitioning. Then, the weight function for constraint c is written as wc(r) =

|ψn⟩ = ;|ψñ ⟩

∑j ∈ c ρj (r) ∑k ρk (r)

⎧ −(r − R )2 ⎫ k ⎬ exp⎨ 2 2π σ 2 ⎩ ⎭ k

a

(5)

i

(8)

atom a. The smooth and all-electron partial wave function are chosen to be identical beyond the cutoff radius rPAW , that is, c outside of the augmentation sphere. Therefore, the pseudo wave functions differ from the all-electron wave functions only in the vicinity of the nuclei. The projectors ⟨p̃ai | are chosen to form a local, complete dual basis to the pseudo partial waves a

a

∑ |ϕĩ ⟩⟨piã | = 1

⟨pjã |ϕĩ ⟩ = δij

i

(9)

With these definitions, the all-electron wave function can be recovered using |ψn⟩ = |ψñ ⟩ +

a

∑ ∑ (|ϕia⟩ − |ϕĩ ⟩)⟨piã |ψñ ⟩ a

ψn(r) = ψñ (r) +

i

∑ [ψna(r) − ψñ a(r)] a

(10a)

(10b)

while all practical calculations can be carried out using the smooth wave functions. Above, ψan(r) − ψ̃ an(r) = ∑i (|ϕai ⟩ − | ϕ̃ ai ⟩)⟨p̃ai |ψ̃ n⟩ is the local atom-wise correction to the smooth pseudo wave function that is nonzero only inside of the augmentation sphere a. It is also convenient to define the atomic density matrices

Nel, k σk

a

∑ ∑ (|ϕia⟩ − |ϕĩ ⟩)⟨piã |

where |ϕ̃ ai ⟩ are the smooth, atomic pseudo partial wave functions and |ϕai ⟩ are the atomic all-electron partial wave functions for

where ρj are spherically symmetric quantities that could be freeatom densities (as in the original Hirshfeld scheme) or iteratively defined atomic densities such as in the iterative Hirshfeld32−34 or the iterative Stockholder schemes.35,36 While these iterative schemes are possible in principle, they require cDFT weight functions that would depend on the electron density and would need to be updated at every SCF cycle, making these iterative schemes demanding. To our knowledge, an iterative scheme within cDFT has not yet been developed. Oberhofer17 observed that the results of cDFT calculations are relatively insensitive to how the atomic densities are defined. In particular, they found that simple Gaussian model densities gave qualitatively the same results as the spherical free-atom densities used in the original Hirshfeld partitioning. In the following, we use the Hirshfeld partition with Gaussian atomic densities. This choice for the weight function has several advantages, such as the Gaussian functions have analytical derivatives making the computation of nuclear derivatives particularly simple (see the Forces in cDFT section). Moreover, the real-space definition of the weight function offers a consistent way to access the real-space grid and LCAO basis sets available in GPAW. Therefore, the current cDFT implementation in GPAW uses the following atomic densities in the weight function ρk (r, R k) =

;=1+

Dija =

∑ ⟨piã |ψñ ⟩fn ⟨ψñ |pjã ⟩

(11)

n

(6)

With these definitions, the all-electron density can be written as

This choice for the weight function introduces one parameter σk for each atom type. This choice affects the decay of the weight function; a small σk results in rapid decrease and a highly localized weight function. Currently, scaled covalent radii are used, but the results hardly vary when σk ∈ (0.5,1.0) Å. Finally, due to the rapid decay and hence efficient screening of spatial contributions along any meaningful periodic dimension, periodic boundary conditions can be treated with the minimum image convention.37 To further increase numerical stability, the weight functions are confined within a real-space cutoff, Rc. Hence, the weight function obtains the form

n(r) = n(̃ r) +

∑ [na(r − Ra) − nã (r − Ra)] (12)

a

where n(̃ r) =

∑ fn ||ψñ ⟩|2 + ∑ nc̃a(r) n

a

(13)

and ñac is the smooth pseudo core density equal to the all-electron core density nac outside of the augmentation sphere. The atomwise local correction takes the following form

⎧ N ⎧ − (r − R )2 ⎫ el, k k ⎪ ⎬ if |r − R k| ≤ Rc exp⎨ ⎪ 2σk 2 ρk (r, R k) = ⎨ σk 2π ⎩ ⎭ ⎪ ⎪0 if |r − R k| > Rc ⎩

[na(r − Ra) − n ã (r − Ra)] =

∑ Dija(ϕja *(r)ϕia(r) − ϕjã *(r)ϕĩa(r)) ij

+ (nca(r) − nc̃a(r))

(7)

nac (r)

To avoid division by 0 in eq 5, the denominator is set to one in regions where the denominator would equal 0. This choice does not affect any calculated value as the numerator equals 0 whenever the denominator is modified.

(14)

∑i ∥ϕai,c⟩|2.

and = Instead of working with the all-electron quantities explicitly, they are obtained by manipulating the atomic corrections. Toward this, a localized compensation charge (Z̃ a(r)) is introduced to correct for the lack of normC


Article

Journal of Chemical Theory and Computation conservation and nuclear charges. The compensation charge is added to the pseudo electron density to give a neutral charge density ρ̃(r) ρ̃(r) = n(̃ r) +

c

∑∫

∑ Z̃ a(r)

a∈c

(15)

a

a

∑ Q LagL̃a (r)

(16)

L

g̃aL(r)

where = is a shape function written in terms of the spherical harmonics YL(r̂). The shape function is bound to fulfill the normalization condition

Now the requirement of identical multipoles between the pseudo and all-electron charge densities can be expressed as

=

∑

ΔaL , i , j Dia, j

∂⟨ψñ σ |

(18)

σ = Vcwcσ f nσ |ψñ σ ⟩ + hatomic |ψñ σ ⟩

where is the spin-dependent atomic cDFT Hamiltonian, which is obtained by expanding the atomic local correction according to eq 14

+ Δ δl 0 (19)

where L is the combined angular (l) and magnetic (m) quantum number. The constants ΔaL,i,j and Δaδl0 are written in terms of the spherical harmonics Y L (r̂) and atomic nuclear charges Z a(r) = −A aδ(r ) with A a being the atomic number. Then the quantities in eq 19 have the following form ΔaL , i , j =

∫ dr YL(r)̂ r l[ϕja*(r)ϕia(r) − ϕjã *(r)ϕĩa(r)]

Δaδl 0 =

∫ dr Y00(r)[̂ −Aaδ(r) + nac(r) − nãc(r)]

σ hatomic |ψñ σ ⟩ =

=

a∈c

dr Vcwcσ (r) ∑a [na(r) − nã (r)] ∂⟨ψñ σ |

a a dr Vcwcσ (r) ∑a ∑i , j Dia, j, σ [ϕiaϕja − ϕi ̃ ϕj̃ ] ∂Dia, j, σ

∂Dia, j, σ

⟨ψñ σ |

∫a∈c dr wcσ(r) ∑ ∑ [ϕiaϕja − ϕĩaϕjã ]f nσ |piã ⟩⟨pjã |ψñ σ⟩ i,j

(25)

This atomic cDFT Hamiltonian is added to the usual atomic KS Hamiltonian, which has an analogous form.38−41 We further note that [ϕai ϕaj − ϕ̃ ai ϕ̃ aj ] is nonzero only inside of the augmentation spheres. Because the augmentation spheres are much smaller than the weight function cutoffs and similar to the Gaussian widths, wσc ≈ ±1 inside of augmentation spheres of the atoms included in wc and 0 otherwise. This allows us to approximate hσatomic as

(20)

σ hatomic ≈ ±Vc

∑ ∑ |piã ⟩ΔOija⟨pjã | a ∈ wc i , j

a

(26)

with ΔOaij = [ϕai ϕaj − ϕ̃ ai ϕ̃ aj ]. This quantity is closely related to the PAW overlap operator Ô , which is defined as Ô = 1 + ∑a ∑i,j | p̃ai ⟩ΔOaij⟨p̃aj |.38−41 All of the quantities used in eq 26 are readily available as atomic PAW parameters, making computation of the atomic Hamiltonian facile. The approximation made for eq 26 holds well even at short atomic distances and overlapping augmentation spheres, as noted in the Test Calculations section for N2. Updating the Lagrange Multiplier. All quantities in eq 2 are now defined in terms of PAW quantities, but an efficient way for determining the constraining potential Vc is needed. Toward this, eq 2 can be written in terms of Vc only because the density is computed from Vc. Van Voorhis has also shown that the functional F[Vc] is strictly concave and the minimum energy of

a a 4π Dia, jY00[ϕja *(r)ϕia(r) − ϕj̃ *(r)ϕĩ (r)]

i,j

+

a∈c

∂∫

a

∫ dr n(r) = ∫ dr n(̃ r) + ∑ ∫ dr Δna(r) ∑

∂∫

= Vc

By noting that the spherical harmonic Y00 = 1/ 4π , the above equations can be used for constructing the integral over the allelectron density. This is very useful as the cDFT energy in eq 2 and stationary condition in eq 27 are written in terms of such integrals. To be explicit, the integral over the all-electron density is

Δna(r) =

(24)

hσatomic

a

i,j

(23)

∂[∑n fn ⟨ψñ σ |Vcwcσ |ψñ σ ⟩ + ∫ dr Vcwcσ (r) ∑a Δna(r) − VcNc]

Inserting eqs 14 and 16 into the above equation, the multipole moments take the following form39 Q La

(22)

Here, only the part related to the cDFT external potential is presented, while other terms in the Hamiltonian can be found in refs 38−41. Within the PAW, the cDFT potential is obtained by using eq 23 on a single constrained region

(17)

∫ dr r l[nã (r) + Z̃a(r) − na(r) − Z a(r)]YL(r)̂ = 0

dr wc(r)Δna(r) − Nc)

δE = fn H̃̂ |ψñ ⟩ δ⟨ψñ |

g̃al (r)YL(r̂)

∫ dr r lgl̃ a (r)YL(r)̂ = 1

Ωa

where the latter integration is performed only inside of the augmentation spheres Ωa belonging to region c. For a self-consistent calculation, the cDFT Hamiltonian must also be presented within the PAW method. The PAW Hamiltonian can written as H̃̂ = ; †Ĥ ; ,38−40 but it can also be obtained from the relation39,41

The compensation charges are constructed by requiring that the pseudo and all-electron charge densities have the same multipole moments. To that end, the compensation charge is defined in terms of multipole moments QaL Z̃ (r) =

∑ Vc(∫ dr wc(r)n(̃ r)+

F[n(r), Vc] = EKS[n(r)] +

4π Y00[−A aδ(r ) + nca(r) − nc̃a(r)] + A a (21)

We again note that Δna(r) can be obtained from atomic coefficients and Di,ja. The pseudo electron density is computed self-consistently and Dai,j is obtained from the resulting pseudoorbitals, {ψ̃ n}. The KS wave functions and energies (with an external potential) are solved using the pseudo wave functions and densities. Therefore, part of the cDFT energy is included in the KS calculation, but we need to correct for the lack of norm conservation in the PAW wave functions. Thus, the following formula is used for computing the cDFT energy D


Article

Journal of Chemical Theory and Computation

discontinuities in cDFT if the underlying DFT energy landscape exhibits transitions between different stable solutions,42 but these were not observed in the present work. Forces in cDFT. Introduction of the constraining potential in eq 2 gives rise to an additional term in the Hellman−Feynman forces

the system under the constraint in eq 4 is found by carrying out a two-stage optimization.7−9 F[Vc] = min max[EKS[n(r)] + n

{Vc}

∑ Vc(∫ dr wc(r)n(r) − Nc)] c

(27)

The outer optimization is done using a regular KS-DFT minimization, for which efficient algorithms exist. The inner optimization requires that F[Vc] is stationary with respect to each Vc7−9 dF = dVc

fA = −∇A F = −∇A EKS −

c

= f AKS + f AcDFT

∫ dr wc(r)n(r) − Nc = 0

dFcDFT[n(r, Vc)] dRA

⎡ = ∑ ⎢ −Vc ⎢ i ⎣

Thus, the total cDFT force is a sum of the usual KS force and cDFT contribution from the weight function and the all-electron density. The cDFT force is divided into pseudo and atomic components based on the cDFT energy expression of eq 22. This is combined with eqs 24 and 26 with the force definition within the PAW method given in refs 38, 39, and 41] to yield the cDFT forces as

∫

∂ ∫ dr Vcwc(r)n(̃ r) ∂RA

|ψñ ⟩

∑∫

∂wc(r, R) − ∂RA

∑∫

∑ ⎢⎢−Vc ∫ dr n(̃ r) i

=−

∂w (r, R) − dr n(̃ r) c ∂RA

⎡

≈

⎣

i,j

i,j

where the same approximation used to obtain eq 26 is made. We stress that the pseudo wave functions are kept fixed and Pulay forces38 are added as done in GPAW (see eq 50 of ref 39). The latter cDFT force component, that is, the atomic correction, is computed with the rest of the atomic Hamiltonian contributions as done in GPAW (see again eq 50 of ref 39) in order to maintain orthogonality between orbitals. However, the first contribution to the forces needs to be handled separately. Using the Hirschfeld weight function defined in eq 7, the partial derivative is obtained as

i,j

Vcwc(r)ΔOia, j

a ∂Vcwc(r)Δna ∂Dij ∂Dija ∂RA

∂Dia, j ⎤ ⎥ ∂RA ⎥⎦

∂Dia, j ⎤ ⎥ ∂RA ⎥⎦

(30)

∂(|r − RA|) r − xA =− x ∂xA |r − RA|

(33)

dr Vc ΔOia, j

∂Θ(Rc − |r − RA|)ρ(|r − RA|) ∂(|r − RA|) = −Θ(r − RA)

ρA (|r − RA|)|r − RA| σA 2

+ ρA (|r − RA|)δ(R c − |r − RA|)

(34)

where Rc is the cutoff radius. The δ terms, so-called surface terms in ref 17, are due to the cutoff radius. To ensure the correctness of the analytical forces, they were compared to forces computed using finite differences described below. The forces can also be easily computed using finite differences. This is especially important if one uses another weight function specification. In this case, the derivate is simply

(31)

where δA∈c = 1 if atom A belongs to weight function i and 0 otherwise. Then, we are left with the calculation of the derivative of the atomic weight function ρA. For a Gaussian weight function, this can be done analytically. The Gaussian weight functions depend on the nuclear coordinates, and the applied cutoff takes the form of a Heaviside function Θ, giving the Gaussian derivative of atom A in the x direction as ∂Θ(Rc − |r − RA|)ρ(|r − RA|) ∂(|r − RA|) ∂(|r − RA|) ∂xA

A∈c

∑ ∫ dr ∑ a

dr

A∈c

−

and the first partial derivative in eq 32 in analytical form is

∑k ρk δA ∈ c − ∑j ∈ c ρj δ − wc ∂wc(r, R) = ρA′ = A ∈ c ρ′ 2 ∑k ρk A ∂RA (∑k ρk )

ρA′ , x =

(29)

(28)

By virtue of eq 27, this stationary state is a maximum. Note that here the all-electron density is used as presented in eq 21. As the first derivatives are easily computed, efficient optimization algorithms can be employed. Also, second derivatives can be used, in which case cDFT response functions should also be evaluated.42 It must be noted that it is possible to observe multiple (self-consistent) solutions, hysteresis, and energy f AcDFT = −

∑ Vc∇A ∫ dr n(r)wc(r)

ρA′ , x ≈

ρA (r, RA + Δx) − ρA (r, RA − Δx) 2Δx

(35)

where Δx is either a small number or an adjacent grid point making the computation of finite difference derivatives efficient. Parametrization of Marcus Theory Using cDFT. Instead of working with the usual Marcus theory, in eq 1, which is valid at the nonadiabatic limit only, it is more convenient to use the Landau−Zener charge transfer equation,5 which combines the nonadiabatic and adiabatic limits

(32)

using the chain rule. The latter partial derivative is E


Article

Journal of Chemical Theory and Computation ⎡ −G‡ ⎤ ka → b = κelνn exp⎢ adia ⎥ ⎢⎣ kBT ⎥⎦

⎛ H′AA H′AB ⎞ ⎟⎟ H′ = ⎜⎜ ′ HBB ′ ⎠ ⎝ HBA

(36)

where the diagonal element is H′AA = = E |Ψ=ΨcDFT and similarly for HBB ′ . The off-diagonal elements can be computed after noting that ψA is an eigenstate of the cDFT-KS Hamiltonian HKS A + VA ∑i w(r)i. Thus, adding and subtracting VA ∑i w(r)i to HKS A yields an expression for the off-diagonal matrix element18

This function depends parametrically on the adiabatic activation free energy (G‡adia), the coupling element between the initial and final diabatic states (Hab), and the reorganization free energy λ. νn is an effective nuclear frequency along the reaction coordinate, and κel is the electron transmission coefficient computed using the Landau−Zener equation5,6,43 ⎧ 2PLZ if ΔG 0 ≥ −λ ⎪ κel = ⎨ 1 + PLZ ⎪ 0 ⎩ 2PLZ(1 − PLZ) if ΔG < −λ

(37)

PLZ = 1 − exp[−2πγ ]

(38)

2πγ =

H′AB = ⟨ψA|HBKS|ψB⟩ = Q BSAB − VBWAB

(39)

where h is Planck’s constant, ΔG is the reaction free energy, and PLZ is the Landau−Zener transition probability for a crossing between two diabatic states that can give both the adiabatic (2πγ ≫ 1, κ ≈ 1) and nonadiabatic (2πγ < 1) limits of electron transfer. It is also worth noting that the transmission coefficient in eq 37 takes different forms for the normal (ΔE ≥ −λ) and inverted (ΔE < −λ) Marcus regions. λ is the reorganization free energy obtained using cDFT 0

λ = ⟨FA(R B)⟩ − ⟨FA(RA)⟩

(40)

H = S−1/2H′S−1/2

where ⟨FA(RA)⟩ is the cDFT free energy computed at the initial state (RA) with the charge constrained at the initial state while ⟨FA(RB)⟩ is computed at the geometry of the final state (RB) but charge-constrained at the initial state. The free energies can be computed by ab initio molecular dynamics utilizing the cDFT machinery.6 The adiabatic transition state energy can be directly computed using normal, unconstrained DFT or approximated with cDFT from the intersection of Marcus parabola, including a correction for adiabaticity44

±

(45)

′ )2 + 4(H′AB + E±SAB)2 ] (H′AA − HBB

(46)

From these eigenvalues, the vertical excitation energy ΔEv can be computed ΔEv = E+ − E− = (41)

′ )2 (H′AA − HBB + 4(H′AB)2 1 − SAB

(47)

One possibility is to compute HAB directly from the above equation as all other quantities can be obtained with cDFT. This “energy gap” approach is less accurate9 than the orthogonalization of the procedure above. Other expressions for the coupling constant based on the explicit form of the coupling constant are obtained by combining the two equations above

where ΔG0 is the free energy difference between the final and initial diabatic states. When Hab is small, the adiabatic barrier in eq 41 reduces to the nonadiabatic barrier. The coupling constant between two orthogonal diabatic states is defined as the offdiagonal element of the electronic Hamiltonian H

Hab = ⟨Ψa|H |Ψb⟩

(44)

The off diagonal elements of H are now the sought coupling constants. Using orthogonalization methods can sometimes significantly distort the eigenstates. In order to circumvent this issue, the coupling constant can be extracted directly from the nonorthogonal 2 × 2 Hamiltonian. In this approach, the secular equation for the nonorthogonal Hamiltonian is solved giving two eigenvalues9,46,47 1 ′ E± = [H′AA + HBB 2

‡ ‡ ≈ GMarcus = ΔG‡ − Δ Gadia

⎡ (λ + ΔG 0)2 λ + ΔG 0 = − ⎢|Hab| + ⎢⎣ 4λ 2 ⎤ (λ + ΔG 0)2 − + |Hab|2 ⎥ ⎥⎦ 4

KS

and analogously for HBA ′ . QA is the cDFT free energy from ⟨ψA| A HKS A + ∑i Vi wi|ψA⟩ for state A, SAB = ⟨ψA|ψB⟩ is the off-diagonal element of the overlap matrix, and WAB = ⟨ψA|∑i wBi (r)|ψB⟩ is the off-diagonal element of the weight matrix. Note that here we are dealing with the all-electron wave functions and PAW factorization of the Hamiltonian and the cDFT external potential is not necessary. The Hamiltonian matrix can be orthogonalized using, for example, Löwdin orthogonalization,9 diagonalizing the dipole moment matrix as in generalized Mulliken−Hush (GMH) treatment45 or rotating the cDFT eigenstates to eigenstates of the weight matrix W.18 Because both the weight and overlap matrices are Hermitian, we use the Lö wdin symmetric orthogonalization, where the nonorthogonal basis is made orthogonal by multiplying by S−1/2. Then, the Hamiltonian in this basis is obtained after a similarity transformation, giving the orthogonal Hamiltonian as

π 3/2 |Hab|2 hνn λkBT

(43)

⟨ψA|HKS A |ψA⟩

(42)

HAB =

In cDFT, one uses the KS Hamiltonian as the electronic Hamiltonian and the two diabatic cDFT-KS wave functions obtained as solutions of eq 2. These states are in general nonorthogonal (primed variables) but can be made orthogonal, as shown later. Within cDFT, the two-state Hamiltonian matrix is written as

′ H′ + HBB 1 H′AB − SAB AA 2 2 1 − SAB

(48)

On the basis of this equation, there are at least two feasible approaches toward computing the coupling constant. First, in the so-called “broken symmetry”9,48,49 approach,HAB ′ is calculated directly from the cDFT diabatic wave functions. In practice, this F


Article

Journal of Chemical Theory and Computation can be done using the corresponding orbital transformation49 or nonorthogonal configuration interaction.50 A second approach47 based on eq 48 is to expand the ground-state wave function (from a DFT calculation) using two cDFT states: Ψg = aΨA + bΨB. Then, the coupling constant from can be expressed in terms of the coefficients found from a = ⟨ΨA|Ψg⟩ and b = ⟨ΨB|Ψg⟩ and the energy gap ΔE = EA − EB HAB =

⎛ ab a 2 + b2 ⎞ 1 E 1 Δ + ⎟ ⎜ 2ab ⎠ 1 − SAB2 a 2 − b2 ⎝

an alternative to the density matrix metric can be based on orbital localization. The localization equals unitary transformation of the cDFT wave functions, ϕA = UψA. As unitary matrices have the property U · U† = U† · U = 1, they can be inserted in the overlap and weight matrix calculations, yielding SAB = ⟨ψA|ψB⟩ = ⟨ψAUA† ·UA|UB†·UBψB⟩ = ⟨ϕA|UA·UB†|ϕB⟩

(49)

Coupling constants computed using either the broken symmetry or the adiabatic wave function approach avoid orthogonalization, can diminish SIE, and can reduce functional dependency of the coupling constant.9 On the other hand, orthogonalization methods have been used with success in previous cDFT implementations.9,18 In the present work, we will utilize the method of Löwdin orthogonalization to enable direct comparison with previous cDFT results. It is also worth discussing that recently cDFT coupling constants have been found51 inaccurate when donor and acceptor wave functions are highly delocalized and their overlap does not decay to 0 at large distances. Such a situation can occur when strong static correlation and (nearly) degenerate orbitals are present. Furthermore, cDFT constrains the density not the orbitals involved in the charge transfer process. The chemical picture behind the observed failures is that, instead of a single electron transfer between two orbitals, fractions of multiple electrons are transferred between several orbitals. To detect when cDFT couplings are expected to be inaccurate, a simple metric based on density matrices of the donor and acceptor was developed in ref 51. However, this metric can only be applied if a LCAO basis set is used as the density matrices are not available in grid or plane-wave bases. If the source of inaccuracy in cDFT coupling originates from delocalized donor/acceptor wave functions, localizing the orbitals could offer a cure to improve the coupling terms. Also,

(50)

and similarly for the weight matrix. The cDFT wave functions are then expanded as Slater determinants, as shown below. The result from < = UA·UB† can be used as a measure of how different the localized donor and acceptor orbitals are. If ψA = ψB, UA = UB and < = I. In this case, all of the eigenvalues of < would equal unity, and the wave functions would be localized very well. If only one orbital from both states changed as a result of the charge transfer, only one of the eigenvalues of < would differ significantly from unity. In this case, < would mix the localized orbital between ϕA and ϕB, reducing the quality of localization. If multiple orbitals and fractional electrons participated in the charge transfer, more eigenvalues of < would differ significantly from unity. The number of eigenvalues differing significantly from unity then acts as a measure for the fractional charge transfer between multiple orbitals. This approach will be tested in future work. Calculation of SAB and WAB. The overlap matrix S can be computed after expanding the KS wave functions as Slater determinants using the antisymmetrizer (̂ , to give |ψB⟩ = (̂ [|ψB1⟩, |ψB2⟩...] = 1/ N! det[B]. In the following, we †

also use the equality (̂ ·(̂ = N! (̂ and that a determinant multiplied by a function is equal to one of the rows being multiplied by this function. Using these rules, the overlap matrix becomes

SAB = ⟨ψA|ψB⟩ =

†

∑ wAkwBk[⟨ψAk1,α(1)|⟨ψAk1,β(2)|⟨ψAk2,α(3)|...](̂ ·(̂ [|ψBk1,α(1)⟩|ψBk1,β(2)⟩|ψBk2,α(3)⟩...] k

N!

=

∑ wAkwBk[⟨ψAk1,α(1)|⟨ψAk1,β(2)|⟨ψAk2,α(3)|...]·(̂ [|ψBk1,α(1)⟩|ψBk1,β(2)⟩|ψBk2,α(3)⟩...] k

N!

=

∑ wAkwBk[⟨ψAk1,α(1)|⟨ψAk1,β(2)|⟨ψAk2,α(3)|...]· k

=

∑

wAkwBkdet[⟨ψAk1,α(1)|ψBk1,α(1)⟩,

1 det[|ψBk1,α(1)⟩|ψBk1,β(2)⟩|ψBk2,β(3)⟩...] N!

⟨ψAk1,β(2)|ψBk2,β(2)⟩...]

k

=

∑ wAkwBkdet[∫ dr nAα1,αB,1k(r), ∫ dr nAβ1,βB,2k(r), ...] k

where |ψk,α An ⟩ is the KS eigenstate of the diabatic state A with a specified k-point and spin dependency. The pair density is k,α k,β defined as nk,α,β nm = ⟨ψAn |ψBm⟩, and the determinants are written using the convention of presenting the diagonal values. The pair density is evaluated using the all-electron wave functions defined in eq 10a. In practice, the all-electron wave functions are obtained by adding the atomic contributions to the pseudo wave functions, which are then interpolated to a dense real-space grid to ensure

(51)

proper description of the cusp at the nuclei. An example of such an all-electron wave function is presented in Figure 2. The initial weight matrix elements Wab and Wba are computed using the weight operators Ŵ ab = ∑c VBc (r)wBc (r) and Ŵ ba = ∑c VAc (r)wAc (r), respectively. To ensure the hermicity of the weight matrix, the final off-diagonal elements are taken as WAB = WBA * = (Wab + Wba)/2. The initial elements of the weight matrix are computed using the pair density defined in eq 51. G


Article


of 0.16 Å. In the case of the diamine cation, the default LCAO basis sets54 were used for initial geometry optimizations. All calculations were carried out using an unrestricted DFT formalism without fixing the spin state. The convergence criteria for the density was 10−4 and 10−5 eV/electron for the energy and 10−6 eV/electron for the eigenstates, unless otherwise stated. The forces were computed analytically using eq 34. Adequacy of the simulation cell was ensured by increasing the cell until the energy was converged. For computation of the weight and overlap integrals, the pseudo wave functions were transformed to (frozen core) all-electron wave functions and presented on a very fine grid with 0.05 Å spacing to capture the rapidly oscillating allelectron wave functions near the nuclei. The coupling constants were computed using Löwdin orthogonalization of eq 45. Furthermore, instead of using free energies, single-point 0 K energies were used in all calculations. The PBE functional55 was used in all calculations. The first test was the charge localization in a positively charged helium dimer, which presents a serious challenge for current DFT functionals. We also computed the coupling constants as a function of distance and compared the obtained values to the literature. Then, we studied the effect of the weight function and the accuracy of the approximation made in eq 26 by studying the charge localization in N2. As a more complicated test, we computed the coupling constants in a Zn2 cation as a function of distance and compared the results to literature values. As a final test for the coupling constants, we considered a benzene−Cl anion. Finally, we studied the charge localization in an N,N′dimethylpiperazine (DMP) cation, where it is shown that cDFT can be used to construct charge-localized states, whereas commonly used generalized gradient and hybrid functionals of DFT fail to reproduce the localized state.

Figure 2. All-electron (AE) and PAW wave functions of the frontier orbitals of He2+ obtained using cDFT with the positive charge localized on the left nucleus at 3 Å. The other is at 5.5 Å, and the x-axis is along the molecules axis. ̂ |ψ ⟩ Wab = ⟨ψA|Wab B =

∑ VcB ∑ wAkwBkdet[∫ dr wcBnAα1, αB,1k(r), ∫ dr wcBnAβ1, βB,2k(r)...] c

k

(52)

Here, it must be noted that for spin constraints, the contribution from β-spin orbitals needs to be multiplied by −1 to ensure consistent treatment with the definition of weight functions in eq 2. Implementation. All of the above methods have been implemented in the GPAW software.39,40 In GPAW, the wave functions can be presented on a real-space grid or LCAO, but the potentials are always presented on a real-space grid. To be consistent with the underlying machinery, we compute all of the quantities on a dense grid. This allows all of the operations to be efficiently parallelized as usually done in GPAW. The constraining potential is added to a PAW-KS calculation as an external potential. Thus, the pseudo wave functions are computed using the modified external potential. As GPAW is coded in Python and relies on Numpy52 and Scipy,53 optimization of the Lagrangian multipliers is carried out using the optimizers included in Scipy’s optimize package. This gives access to a variety of methods, including Conjugate-Gradient, dogleg, and Quasi-Newton, among others. Using the Gaussian weight function enables the computation of the forces analytically and efficiently. Both the finite-difference and the analytical approach yield very similar forces, but the analytical computation is about 10 times faster. While a cDFT calculation can be carried out with LCAO or grid basis sets, the coupling constant is currently available only in the grid-mode.

■

■

TEST CALCULATIONS (He2+ Dissociation and Coupling Constant. The first test case is a positively charged helium dimer, which, despite its seeming simplicity, is troublesome for DFT; this is due to SIE and related charge overdelocalization.1 Current functionals predict the charge-delocalized state to be more stable than the localized state, which is, simply put, incorrect.56 In Figure 3, we show the dissociation curve for He2+ where the +1 charge is constrained on one of the He atoms and we report the energy relative to the isolated fragments. It can be seen that the current cDFT implementation correctly captures the interaction at all distances as the interaction energy is very close to 0. The discrepancies at shorter distances (not studied here) are due to formation of a chemical bond between the fragments, and constraining a unit charge on one of the He atoms becomes unphysical. In such cases, the promolecular approach by Van Voorhis9 would give a greatly improved charge constraint. We also note that constraining a +1 spin on one of the He atoms gives identical results to the charge constraint, as it should be. Also, the coupling constants from eqs 44 and 45 are shown in Figure 3. For this system, very accurate coupling constants from GMH treatment at the FCI/aug-cc-pV5Z level57 and previous cDFT implementation18 in the CPMD code are available. In general, our results are in very good agreement with both previous studies for all distances spanning 5 orders of magnitude in Hab. As reported previously,18 the coupling constant exhibits exponential decay characteristic for tunnelling: Hab ≈ exp[−rβ/2]. From our calculations, β = 4.19, which compares well

COMPUTATIONAL DETAILS

We tested our implementation on several relevant systems that have been previously investigated with other cDFT implementations available in the literature and high-level wave function methods. All calculations were carried out using the GPAW software.39,40 The final results were obtained using a grid spacing H


Article


as the N−N bond length is only 1.12 Å and the details of the weight functions are significant, and as a result, charge localization becomes rather unphysical. By varying the width of the Gaussian in the weight function by 30%, the cDFT energies can be varied by up to 5.0 eV. When the distance between the nitrogen atoms is increased to 2.0 Å, a similar change in the weight function width changes the energies by only 0.05 eV. A similar test for the charge reorganization in CO− as studied by the CPMD implementation17 also shows the sensitivity of the cDFT energies on weight function details, which decrease as atomic distances are increased. Zn2+ Coupling Constants. As another, more stringent test on the accuracy of the current cDFT implementation for computation of the coupling constants, we computed the couplings of a positively charged Zn dimer as a function of the interatomic distance. The results are given in Table 1 and Figure 3. (Upper) Dissociation curve for He2+ using DFT and cDFT. (Lower) Coupling constant of He2+ from cDFT.

Table 1. Comparison of Coupling Constants in Zn2+ as a Function of Distance to Different cDFT Implementations with Various Weight Functions and to High-Level GMH Dataa

with of values β = 4.82 and 4.98 from GMH calculations57 and the CPMD implementation of cDFT,18 respectively. In general, our values are slightly larger (0.01−5 mHa) than those obtained in previous studies,18,57 but these very small difference are most likely due to different choices for the weight function. N2. Next, the charge reorganization in N2 is considered. In this case, we vary the charge on one of the nitrogen atoms from +1 to −1, and the results are shown in Figure 4. The energy exhibits

r[Å] Becke8 Löwdin8 Hirschfeld18 GMH45 this work

5

6

7

8

9

β

5.87 15.0

1.52 6.51

0.400 2.29

0.105 0.537

0.026 0.101

7.26 16.71

2.16 4.83

0.623 1.35

0.171 0.39

0.044 0.11

2.70 2.50 2.66 2.55 2.51

a

Becke, Löwdin, and Hirschfeld refer to different weight function definitions used in different cDFT work. β is the exponential decay constant, and the unit of Hab is mHa.

compared to values from other cDFT implementations and highlevel GMH data. We again note that our results are in excellent agreement with previous work, and especially, the exponential decay constant β agrees perfectly with GMH and cDFT data. The coupling constants are slightly larger than those obtained with GMH and Becke cDFT but agree very well with the ones obtained using the Löwdin weight function. This also shows the dependency on the details of the weight function even at large atomic distances. Coupling Constants in (Benzene−Cl)−. The final test for the coupling constants is an electron transfer reaction between Cl and benzene in two different configurations. The geometries were taken from ref 45. The Cl atom is positioned 3 Å above the plane spanned by the benzene molecule and 0.604 or 1.208 Å offset from the six-fold axis of benzene. The results are given in Table 2, and it can be seen that the coupling constants are larger than those with other cDFT implementations and GMH. We believe that these discrepancies are due to weight function choices and difficulties observed during converge of the SCF,

Figure 4. (Upper) Energy of N2 with the specified charge placed at a nitrogen atom. The curve is a second-order polynomial fit to the cDFT energies. The red star is the DFT energy. (Lower) Value of the Lagrange multiplier. The line is a linear fit to the points. The N−N bond length is 1.12 Å.

parabolic dependence on the constrained charge, while the Lagrange multiplier depends linearly on the constrained charge. When the constrained charge is set to 0, the DFT energy is recovered, which also shows that the assumption made in eq 26 is valid even at short atomic distances and when the augmentation spheres overlap. The same system has been considered by Van Voorhis9 as an example of a system where details of the weight function are important. The energies of the system deviate from those presented by Van Voorhis, who also noticed that different weight functions give very different energies. However, this is expected

Table 2. Coupling Constant in Two Different (Benzene−Cl)− Configurations As Explained in the Texta r = 0.604 Å r = 1.208 Å

Hirschfeld18

FO18

Becke8

GMH45

this work

55.9 52.3

17.2/14.0 12.4/10.5

48.8 56.1

51.0 51.9

65.1 80.2

a

Becke and Hirschfeld refer to different weight function definitions used in different cDFT work, while FO is the fragment orbital method with two degenerate benzene HOMO orbitals. The unit of Hab is mHa. I


Article

Journal of Chemical Theory and Computation forcing us to use rather loose convergence criteria for the electronic structure and cDFT (density: 10−4; energy: 10−3 eV/ electron; and eigenstates: 10−3 eV/electron). Still, the results are in the same range as those from other methods, and we conclude that the current implementation exhibits satisfactory behavior also in this case.

Table 3. Comparison of the DMP-D+ and DMP-L+ Geometriesa

CHARGE LOCALIZATION IN AN N,N′-DIMETHYLPIPERAZINE CATION Recently, the relative energies of two distinct charge states of an DMP cation were experimentally measured and characterized:58 one state with the charge delocalized between the two nitrogens (DMP-D+) and another with the charge localized primarily on one of the nitrogens (DMP-L+). This system serves as a stringent test case for the SIE1 inherent to practical implementations of KS-DFT. While CCSD(T) can produce the experimentally observed states and predict accurate relative energies, all commonly used KS-DFT functionals (including hybrids) fail to capture the charge-localized state. Due to the SIE, practical implementations of KS-DFT have problems reproducing various chemical and physical properties due to overemphasis on delocalization. Apart from the case in point, this includes defect states and band gaps in crystals, the energy of the highest occupied molecular orbital not corresponding to the ionization potential, and the incorrect description of weakly bound electrons (Rydberg states). One way to address the SIE is the Perdew−Zunger SIC (PZ-SIC),59 which has been successful in describing challenging systems where KS-DFT fails. These systems include defect states in crystals60 and the binding of electrons in Rydberg states.61−63 Here, we test whether the cDFT implementation can capture the localized state and relative energies of the DMP cation. To this end, we compare unconstrained DFT with cDFT, where 0.7e are localized on one of the nitrogens and neighboring carbon atoms. The charge constraint of 0.7e is based on a Bader analysis of the PZ-SIC calculation on DMP-L+ to enable direct comparison with PZ-SIC and CCSD(T) calculations in ref 58. Taking the charge constraint from a previous calculation demonstrates one caveat of cDFT; the constraint requires knowledge of the electronic structure, or the constraint needs to be chosen based on chemical and/or physical intuition. Another possibility would be to use the promolecular densities as suggested in ref 9. In CCSD(T) and other high-level wave function methods, such a decision is not needed. In SIC functionals, a system-dependent prefactor to scale the SICs is often applied (see, e.g., ref 60 eq 12). We first optimized the structure with DFT to obtain DMP-D+ and, starting from this geometry, optimized the structure with cDFT to get the DMP-L+ geometry. The optimized geometries are given in Figure 5, and some characteristic geometric features are tabulated and compared in Table 3. DMP-D+ has both methyl groups symmetrically at pseudoaxial positions, and bond angles

a CNC is the dihedral angle in CH2−CH2−N−CH3, and CN is the distance of the N−CH3 bond (in Å). loc refers to the nitrogen where the charge is localized in DMP-L+.

DMP-D PZ-SIC58 current work

■

DMP-L

CCNC

CN

CCNC

CCNC(loc)

CN

CN(loc)

277 274

1.44 1.46

190 217

258 268

1.45 1.45

1.44 1.45

and lengths agree very well with the PZ-SIC results. In DMP-L+, the symmetry is broken and the methyl group next to the chargelocalized nitrogen (top nitrogen in Figure 5) takes a pseudoaxial position, and the dihedral angle and length again agree very well with PZ-SIC. The only notable difference between the PZ-SIC and cDFT geometries is the angle between the methyl group and the more neutral nitrogen atom. cDFT predicts the methyl group to be in a more axial position than PZ-SIC. Overall, the cDFT geometries compare well with PZ-SIC. According to cDFT, the energy of DMP-L+ is 0.31 eV higher than DMP-D+, while 0.33, 0.34, 0.38 eV differences are obtained by experiments, SIC, and CCSD(T), respectively. This shows that cDFT is a powerful tool to reduce the SIE and to obtain charge-localized states with physically sound energies and geometries. It is also worth noting that the cDFT is only 5−10 times more computationally intensive than a regular GGA calculation, while SIC with its current implementation in GPAW64 is 100 times more expensive and CCSD(T) is 1.000 times more expensive. Thus, cDFT offers a computationally feasible and appealing alternative for studying charge-localized states with good accuracy in a pragmatic way. cDFT offers one significant advantage over DFT and standard wave function methods: it can produce diabatic states. Here, we use this property for parametrizing the Marcus theory of charge transfer, shown in eq 36, to compute the rate for DMP-L+ → DMP-D+ transition. The needed parameters are easily extracted: λ = 1.12 eV, ΔE = 0.31 eV, and Hab = 3.9 × 10−5 eV. Using these values, the adiabatic reaction barrier from eq 41 gives E‡adia = 0.145 eV, which is in very good agreement with the PZ-SIC barrier58 of 0.2 eV. To compute the rate constant, the effective frequency is taken as the CCNC-dihedral bending, which is the reaction coordinate. This gives νn = 5.68 × 1012 s−1. The reaction temperature can be inferred from time-resolved photoelectron spectroscopy experiments of DMP in ref 58, where the effective vibrational temperature was evaluated to be 565−980 K. We chose an intermediate temperature of 750 K. On the basis of these values, the Landau−Zener factor and electronic transmission coeficienct clearly indicate that the reaction is diabatic as PLZ = 9.6 × 10−8 and κ = 9.6 × 10−8. This is in agreement with the minimumenergy pathway (Figure 3 of ref 58), which shows a cusp at the transition state, indicative of small coupling and diabacity. Finally, the rate can be evaluated giving kdiabatic = 1.2 × 105 s−1 for the reaction rate constant. Without taking the diabaticity into account, the rate constant would be kadiabatic = 6.1 × 1011 s−1, demonstrating the relevance of including diabatic effects in charge transfer.

■

CONCLUSIONS We have presented the implementation of constrained density functional theory in the GPAW code utilizing modern, nonnorm-conserving PAW pseudopotentials. The implementation

Figure 5. Optimized geometries of DMP-D+ (left) and DMP-L+ (right). In DMP-L+, 0.7 electrons are localized on the lower nitrogen atom. J


Article

Journal of Chemical Theory and Computation has flexible boundary conditions and can use LCAO and finitedifference grids to present the wave functions. We have also implemented the tools for computing the reorganization energy and coupling constant needed for parametrizing the Marcus model of charge transfer kinetics. The implemented methods allow us to study charge transfer kinetics in a very large variety of systems ranging from molecules to the solid state in a flexible way. The implementation is benchmarked against previous cDFT implementations and high-level quantum chemical calculations. In general, the energies and coupling constants are in very good to satisfactory agreement with previous studies, and the slight differences are most likely due to details of the weight function. The advantages of the current implementation are (1) the PAW pseudopotentials, which are very accurate and allow the extraction of the all-electron quantities, (2) the use of both LCAO and grid to present the wave functions, for quick and accurate computations, respectively, and (3) the flexible boundary conditions, which allow the treatment of molecules, surfaces, and the bulk in a consistent way. The method is applied to study charge localization in the DMP cation where a delocalized and a localized states are present. Whereas currently used DFT functionals fail to reproduce the localized state, cDFT is shown to give the localized state energy, geometry, and reaction energy barrier in excellent agreement with the computationally intensive CCSD(T) and a selfinteraction-corrected PZ-SIC functional. This shows that cDFT offers a computationally feasible method to reduce SIE and to study charge-localized states as well as charge transfer kinetics.

■

(12) Podloucky, R.; Zeller, R.; Dederichs, P. H. Phys. Rev. B: Condens. Matter Mater. Phys. 1980, 22, 5777−5790. (13) Souza, A. M.; Rungger, I.; Pemmaraju, C. D.; Schwingenschloegl, U.; Sanvito, S. Phys. Rev. B: Condens. Matter Mater. Phys. 2013, 88, 165112. (14) Wu, Q.; Kaduk, B.; Van Voorhis, T. J. Chem. Phys. 2009, 130, 034109. (15) Ř ezác,̌ J.; Lévy, B.; Demachy, I.; de la Lande, A. J. Chem. Theory Comput. 2012, 8, 418−427. (16) Ramos, P.; Pavanello, M. Phys. Chem. Chem. Phys. 2016, 18, 21172−21178. (17) Oberhofer, H.; Blumberger, J. J. Chem. Phys. 2009, 131, 064101. (18) Oberhofer, H.; Blumberger, J. J. Chem. Phys. 2010, 133, 244105. (19) Ghosh, P.; Gebauer, R. J. Chem. Phys. 2010, 132, 104102. (20) Ma, P.-W.; Dudarev, S. L. Phys. Rev. B: Condens. Matter Mater. Phys. 2015, 91, 054420. (21) Sena, A. M. P.; Miyazaki, T.; Bowler, D. R. J. Chem. Theory Comput. 2011, 7, 884−889. (22) Ratcliff, L. E.; Grisanti, L.; Genovese, L.; Deutsch, T.; Neumann, T.; Danilov, D.; Wenzel, W.; Beljonne, D.; Cornil, J. J. Chem. Theory Comput. 2015, 11, 2077−2086. (23) Wesolowski, T. A.; Warshel, A. J. Phys. Chem. 1993, 97, 8050− 8053. (24) Goedecker, S.; Teter, M.; Hutter, J. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54, 1703−1710. (25) Troullier, N.; Martins, J. L. Phys. Rev. B: Condens. Matter Mater. Phys. 1991, 43, 1993−2006. (26) Lejaeghere, K.; et al. Science 2016, 351, 1415. (27) Comparing Solid State DFT Codes, Basis Sets and Potentials. https://molmod.ugent.be/deltacodesdft, (accessed July 13,2016). (28) Dal Corso, A. Comput. Mater. Sci. 2014, 95, 337−350. (29) Blöchl, P. E.; Först, C. J.; Schimpl, J. Bull. Mater. Sci. 2003, 26, 33− 41. (30) Calle-Vallejo, F.; Koper, M. T. Electrochim. Acta 2012, 84, 3−11. (31) Hirshfeld, F. L. Theor. Chim. Acta 1977, 44, 129−138. (32) Bultinck, P.; Van Alsenoy, C.; Ayers, P. W.; Carbo-Dorca, R. J. Chem. Phys. 2007, 126, 144111. (33) Bultinck, P.; Ayers, P. W.; Fias, S.; Tiels, K.; Van Alsenoy, C. Chem. Phys. Lett. 2007, 444, 205. (34) Verstraelen, T.; Ayers, P. W.; Van Speybroeck, V.; Waroquier, M. J. Chem. Theory Comput. 2013, 9, 2221. (35) Lillestolen, T. C.; Wheatley, R. J. Chem. Commun. 2008, 2008, 5909. (36) Lillestolen, T. C.; Wheatley, R. J. J. Chem. Phys. 2009, 131, 144101. (37) Allen, M.; Tildesley, D. Computer Simulation of Liquids; Clarendon Press: Oxford, U.K., 1987. (38) Blöchl, P. E. Phys. Rev. B: Condens. Matter Mater. Phys. 1994, 50, 17953−17979. (39) Mortensen, J. J.; Hansen, L. B.; Jacobsen, K. W. Phys. Rev. B: Condens. Matter Mater. Phys. 2005, 71, 035109. (40) Enkovaara, J.; et al. J. Phys.: Condens. Matter 2010, 22, 253202. (41) Rostgaard, C. The Projector Augmented-wave Method. arXiv, https://arxiv.org/abs/0910.1921 (accessed July 15, 2016). (42) O’Regan, D. D.; Teobaldi, G. Phys. Rev. B: Condens. Matter Mater. Phys. 2016, 94, 035159. (43) Newton, M. D.; Sutin, N. Annu. Rev. Phys. Chem. 1984, 35, 437− 480. (44) Oberhofer, H.; Blumberger, J. Phys. Chem. Chem. Phys. 2012, 14, 13846−13852. (45) Cave, R. J.; Newton, M. D. Chem. Phys. Lett. 1996, 249, 15−19. (46) Newton, M. D. Chem. Rev. 1991, 91, 767−792. (47) Migliore, A. J. Chem. Theory Comput. 2011, 7, 1712−1725. (48) Farazdel, A.; Dupuis, M.; Clementi, E.; Aviram, A. J. Am. Chem. Soc. 1990, 112, 4206−4214. (49) King, H. F.; Stanton, R. E.; Kim, H.; Wyatt, R. E.; Parr, R. G. J. Chem. Phys. 1967, 47, 1936−1941. (50) Thom, A. J. W.; Head-Gordon, M. J. Chem. Phys. 2009, 131, 124113.

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected] (M.M.). *E-mail: [email protected] (J.M.G.L.). Funding

We acknowledge support from the Villum Foundation’s Young Investigator Programme (fourth round, project: In silico design of efficient materials for next generation batteries; Grant Number: 10096). Notes

The authors declare no competing financial interest.

■

REFERENCES

(1) Cohen, A. J.; Mori-Sánchez, P.; Yang, W. Science 2008, 321, 792− 794. (2) Marcus, R. A. J. Chem. Phys. 1956, 24, 966−978. (3) Marcus, R.; Sutin, N. Biochim. Biophys. Acta, Rev. Bioenerg. 1985, 811, 265−322. (4) Dogonadze, R. R.; Kuznetsov, A. M.; Chernenko, A. A. Russ. Chem. Rev. 1965, 34, 759. (5) Kuznetsov, A. M.; Ulstrup, J. Electron Transfer in Chemistry and Biology: An Introduction to the Theory; John Wiley and Sons, 1999. (6) Blumberger, J. Chem. Rev. 2015, 115, 11191−11238. (7) Wu, Q.; Van Voorhis, T. Phys. Rev. A: At., Mol., Opt. Phys. 2005, 72, 024502. (8) Wu, Q.; Van Voorhis, T. J. Chem. Phys. 2006, 125, 164105. (9) Kaduk, B.; Kowalczyk, T.; Van Voorhis, T. Chem. Rev. 2012, 112, 321−370. (10) Kubas, A.; Hoffmann, F.; Heck, A.; Oberhofer, H.; Elstner, M.; Blumberger, J. J. Chem. Phys. 2014, 140, 104105. (11) Kubas, A.; Gajdos, F.; Heck, A.; Oberhofer, H.; Elstner, M.; Blumberger, J. Phys. Chem. Chem. Phys. 2015, 17, 14342−14354. K


Article

Journal of Chemical Theory and Computation (51) Mavros, M. G.; Van Voorhis, T. J. Chem. Phys. 2015, 143.23110210.1063/1.4938103 (52) van der Walt, S.; Colbert, S. C.; Varoquaux, G. Comput. Sci. Eng. 2011, 13, 22−30. (53) Jones, E.; Oliphant, T.; Peterson, P. SciPy: Open source scientific tools for Python. http://www.scipy.org/, (accessed May 2, 2016). (54) Larsen, A. H.; Vanin, M.; Mortensen, J. J.; Thygesen, K. S.; Jacobsen, K. W. Phys. Rev. B: Condens. Matter Mater. Phys. 2009, 80, 195112. (55) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 77, 3865−3868. (56) Helgaker, T.; Jorgensen, P.; Olsen, J. Molecular electronic structure theory; Wiley, 2012; pp 164−174. (57) Pieniazek, P. A.; Arnstein, S. A.; Bradforth, S. E.; Krylov, A. I.; Sherrill, C. D. J. Chem. Phys. 2007, 127, 164110. (58) Cheng, X.; Zhang, Y.; Jónsson, E.; Jónsson, H.; Weber, P. M. Nat. Commun. 2016, 7, 11013. (59) Perdew, J. P.; Zunger, A. Phys. Rev. B: Condens. Matter Mater. Phys. 1981, 23, 5048. (60) Gudmundsdóttir, H.; Jónsson, E. O.; Jónsson, H. New J. Phys. 2015, 17, 083006. (61) Gudmundsdóttir, H.; Zhang, Y.; Weber, P. M.; Jónsson, H. J. Chem. Phys. 2014, 141, 234308. (62) Gudmundsdóttir, H.; Zhang, Y.; Weber, P. M.; Jónsson, H. J. Chem. Phys. 2013, 139, 194102. (63) Cheng, X.; Zhang, Y.; Deb, S.; Minitti, M. P.; Gao, Y.; Jónsson, H.; Weber, P. M. Chem. Sci. 2014, 5, 4394. (64) Gudmundsdóttir, H.; Jónsson, E. Ö .; Jónsson, H. New J. Phys. 2015, 17, 083006.

L


Implementation of Constrained DFT for Computing Charge Transfer

Recommend Documents