Electronic couplings for molecular charge transfer: benchmarking CDFT, FODFT, and FODFTB against high-level ab initio calculations.

Electronic couplings for molecular charge transfer: Benchmarking CDFT, FODFT, and FODFTB against high-level ab initio calculations Adam Kubas, Felix Hoffmann, Alexander Heck, Harald Oberhofer, Marcus Elstner, and Jochen Blumberger Citation: The Journal of Chemical Physics 140, 104105 (2014); doi: 10.1063/1.4867077 View online: http://dx.doi.org/10.1063/1.4867077 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/140/10?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Ab initio benchmark study for the oxidative addition of CH 4 to Pd: Importance of basis-set flexibility and polarization J. Chem. Phys. 121, 9982 (2004); 10.1063/1.1792151 An ab initio and experimental study of vibrational effects in low energy O + +C 2 H 2 charge-transfer collisions J. Chem. Phys. 115, 3184 (2001); 10.1063/1.1385793 Density-functional study of intramolecular ferromagnetic interaction through m-phenylene coupling unit (II): Examination of functional dependence J. Chem. Phys. 113, 10486 (2000); 10.1063/1.1290008 Comparison of methods for calculating the properties of intramolecular hydrogen bonds. Excited state proton transfer J. Chem. Phys. 111, 849 (1999); 10.1063/1.479371 High level ab initio molecular orbital study of the structures and vibrational spectra of CHBr + and CBr + J. Chem. Phys. 109, 134 (1998); 10.1063/1.476530

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 209.183.183.254 On: Mon, 24 Nov 2014 03:43:53

THE JOURNAL OF CHEMICAL PHYSICS 140, 104105 (2014)

Electronic couplings for molecular charge transfer: Benchmarking CDFT, FODFT, and FODFTB against high-level ab initio calculations Adam Kubas,1,a) Felix Hoffmann,1,2,a) Alexander Heck,3,a) Harald Oberhofer,4 Marcus Elstner,3 and Jochen Blumberger1,b) 1

Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, United Kingdom 2 Lehrstuhl für Theoretische Chemie, Ruhr-Universität Bochum, Universitätsstr. 150, 44801 Bochum, Germany 3 Institute of Physical Chemistry, Karlsruhe Institute of Technology, Fritz-Haber-Weg 6, 76131 Karlsruhe, Germany 4 Department of Chemistry, Technical University of Munich, Lichtenbergstr. 4, 85747 Garching, Germany

(Received 7 November 2013; accepted 17 February 2014; published online 11 March 2014) We introduce a database (HAB11) of electronic coupling matrix elements (Hab ) for electron transfer in 11 π -conjugated organic homo-dimer cations. High-level ab inito calculations at the multireference configuration interaction MRCI+Q level of theory, n-electron valence state perturbation theory NEVPT2, and (spin-component scaled) approximate coupled cluster model (SCS)-CC2 are reported for this database to assess the performance of three DFT methods of decreasing computational cost, including constrained density functional theory (CDFT), fragment-orbital DFT (FODFT), and selfconsistent charge density functional tight-binding (FODFTB). We find that the CDFT approach in combination with a modified PBE functional containing 50% Hartree-Fock exchange gives best results for absolute Hab values (mean relative unsigned error = 5.3%) and exponential distance decay constants β (4.3%). CDFT in combination with pure PBE overestimates couplings by 38.7% due to a too diffuse excess charge distribution, whereas the economic FODFT and highly cost-effective FODFTB methods underestimate couplings by 37.6% and 42.4%, respectively, due to neglect of interaction between donor and acceptor. The errors are systematic, however, and can be significantly reduced by applying a uniform scaling factor for each method. Applications to dimers outside the database, specifically rotated thiophene dimers and larger acenes up to pentacene, suggests that the same scaling procedure significantly improves the FODFT and FODFTB results for larger π -conjugated systems relevant to organic semiconductors and DNA. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4867077] I. INTRODUCTION

Diabatic electronic states are an important concept in chemical physics. They are the natural starting point for the description of electron transfer (ET) reactions1–3 and they have been successfully used also in the framework of (empirical) valence bond theory ((E)VB) for the description of adiabatic chemical reactions4, 5 and proton transport.6 The appealing feature of diabatic electronic states is their chemically intuitive, spatially localized character. This latter property can be highly beneficial, e.g., in calculations on large systems where low scaling of the computational cost with system size is desired (leading, e.g., to tight binding methods), or where electronic structure methods struggle to correctly describe charge delocalized adiabatic states, e.g., in density functional theory (DFT) with approximate exchange-correlation functionals.7 The cost one has to pay is that diabatic states are not eigenfunctions of the electronic Hamiltonian, but require the calculation of electronic off-diagonal elements, also denoted as Hab , Hab = ψa |H|ψb ,

(1)

a) A. Kubas, F. Hoffmann, and A. Heck contributed equally to this work. b) Author to whom correspondence should be addressed. Electronic mail:

[email protected] 0021-9606/2014/140(10)/104105/21/$30.00

where H is the electronic Hamiltonian and ψ a are ψ b are the diabatic states a and b. In ET theory, where only two diabatic states are considered (initial and final ET state), the off-diagonal element Hab , also denoted as electronic coupling matrix element or charge transfer integral, takes a prominent role as the ET rate is proportional to |Hab |2 for small values of Hab (non-adiabatic limit) and the ET activation barrier is lowered by Hab for large values of Hab (adiabatic limit).1, 3 If one choses to work with diabatic states one has to make a choice with regard to the method for generation of diabatic states and the electronic structure method for the actual calculation. There are several methods available to diabatize the ET states of the system such as block diagonalization of the adiabatic electronic Hamiltonian,8–10 generalized Mulliken– Hush method (GMH),10–17 fragment charge difference,18 fragment energy difference,19 fragment orbital,20–23 projection methods,24, 25 constrained density functional theory (CDFT),26–30 and frozen density embedding31, 32 (see also Refs. 33 and 34 for a discussion). Accurate electronic structure methods have been used mostly in the framework of GMH theory, but of course, only for relatively small systems. In the pioneering study of Cave and Newton complete active space self-consistent field (CASSCF) calculations were

140, 104105-1

© 2014 AIP Publishing LLC


104105-2

Kubas et al.

carried out to obtain electronic coupling matrix elements for − 10, 11 More recently, Zn+ 2 , benzene-Cl and similar systems. benchmark-quality data were provided by Pieniazek et al.35 for He+ 2 , at the full configuration interaction (FCI) level and on small systems composed of first and second row atoms such as the ethylene dimer radical cation. Pavenello et al.32 used the equation-of-motion coupled cluster singles and doubles method (EOM-CCSD) to test their frozen density embedding scheme on the ethylene and water dimer radical cation. Moreover, Blancafort and Voityuk performed extensive CASPT2 calculations on hole36 and electron37 transfer between stacked DNA nucleobases. It is clear that computationally less expensive electronic structure methods need to be used if one would like to study larger molecular systems that are of high practical significance in materials science and biology, or if coupling calculations need to be carried out for a very large number of times as, e.g., in FODFTB non-adiabatic MD simulations38, 39 or kinetic Monte Carlo simulations of charge transport,40–43 To this end, we have previously used more economic schemes for estimation of Hab based on constrained density functional theory (CDFT),30 fragment orbital density functional theory (FODFT),30, 43–46 and density functional tight-binding (FODFTB).38, 47–53 Especially the latter two methods can be used when a very large number of calculations on large molecular systems is required, such as, e.g., in organic semiconductors, DNA and biological redox cofactors.38, 43, 45–47, 49, 51–53 However, it remains unclear how accurate these methods are, partly because of the limited availability of high-level ab initio reference data. In our previous work we have tested the CDFT and FODFT methods against high-level ab initio calculations and found surprisingly good agreement for small diatomics + − 30 But for larger systems such He+ 2 , Zn2 and benzene-Cl . as the mixed-valence tetrathiafulvalenediquinone anion30 and for hole transfer between oxygen vacancies in crystalline MgO,54, 55 we found that the results obtained with CDFT were rather sensitive on the fraction of exact exchange used to approximate the exchange-correlation functional. The FODFTB couplings were tested for hole transfer in DNA where they were found to be comparable to data obtained by CASPT2 calculations using the GMH method.47 Here we introduce the HAB11 database that consists of 11 small and medium sized π -conjugated organic dimers (see Figure 1 for dimer alignment) with increasing number of multiple bonds and varying number of heteroatoms (N, O, S atoms, see Table I). The choice of molecules reflects our interest in describing charge transport in organic semiconducting materials and DNA. The molecules are small enough so that electronic coupling matrix elements between dimer

J. Chem. Phys. 140, 104105 (2014)

structures can be obtained with high-level ab initio quantum chemistry methods at the MRCI+Q and NEVPT2 levels of theory. These values will provide new benchmarks on which the performance of our CDFT, FODFT, and FODFTB implementations as well as the methods of other groups can be evaluated (the atomic coordinates of the molecular dimers are appended in the supplementary material).56 Here we focus in particular on the performance of the above mentioned three DFT-based methods in reproducing absolute Hab values for electron transfer in cationic systems (the term “electron transfer” is used throughout the paper as there is no bridge between donor and acceptor to distinguish between electron and hole transfer) and the corresponding exponential distance decay constants β, |Hab | = A exp(−βd/2).

(2)

We find that CDFT reproduces the benchmark values very well if it is used in combination with a functional that includes 50% Hartree-Fock exchange. FODFT and FODFTB give reasonable results; the couplings are consistently too small and the decay constants are too large. However, the errors of the latter two methods are systematic within the chemical space of our database and can be simply corrected by applying a uniform scaling factor. We also present applications to larger systems with increasing number of condensed benzene rings (acenes, naphthalene to pentacene) and find that the results obtained for the molecules in the database can likely be generalized to larger π -conjugated systems relevant to organic semiconducting materials. This paper is organized is follows. In the Sec. II we describe the HAB11 database. Sections III and IV contain details of the various methods used for calculation of Hab . To provide rigorous reference values, we exploit the fact that for symmetric dimers Hab is uniquely defined as half of the first adiabatic excitation energy. The latter is obtained from MRCI+Q and NEVPT2 reference calculations with adequately tested basis sets. In addition, we will investigate the performance of more economic CC2 and SCS-CC2 ab initio calculations, which will be used as a reference for the larger acene dimers outside the database. We then briefly review the theory underlying the CDFT, FODFT, and FODFTB approaches and give all the computational details. In Sec. V, we first present and discuss the performance of CDFT, paying particular attention to the dependence on the percentage of Hartree-Fock exchange used. We then discuss the application of a simple scaling factor that considerably improves the performance of FODFT and FODFTB. After the calculations on the stacked dimers in the database, two case examples are considered to probe whether our conclusions carry over to more general situations: (i) the orientational dependence of Hab in a thiophene dimer and (ii) the electronic couplings in larger acene dimers. In Sec. VI we discuss the relation between electronic couplings obtained on the various levels of approximations. Finally, we draw our conclusions in Sec. VII. II. HAB11 DATABASE

FIG. 1. Schematic representation of the alignment of a dimer. The neutral geometry was replicated along the axis perpendicular to the plane intercepting all heavy atoms. Hab values were obtained for d = 5.0, 4.5, 4.0, 3.5 Å.

In this study, 11 homo-dimers of organic molecules56 were selected with different numbers of π -bonds and


104105-3

Kubas et al.

J. Chem. Phys. 140, 104105 (2014)

TABLE I. HAB11 database of organic molecules56 for the calculation of electron transfer electronic coupling matrix elements. For simplification only monomers are presented but the actual calculations were performed on homo-dimeric π -stacked cationic structures. For each entry the molecular point group, electronic states included in the state-averaging procedure, active space in multireference calculations, and reference method used to estimate couplings are given. The multiplicity of each electronic state is 2 (omitted for clarity) and the ground state term is given in bold. Name

Structure

Molecular point group

Electronic states

Active space

Reference

Ethylene

D2h

Ag , B2u

(3,4)

MRCI+Q

Acetylene

D2h

Ag , B3g , B1u , B2u

(7,8)

MRCI+Q

Cyclopropene

C2v

A1 , B2

(3,4)

MRCI+Q

Cyclobutadiene

D2h

B2g , B3u

(7,8)

MRCI+Q

Cyclopentadiene

C2v

A2 , B1

(7,8)

MRCI+Q

Furane

C2v

A2 , B2

(11,10)

MRCI+Q

Pyrrole

C2v

A2 , B2

(11,10)

MRCI+Q

Thiophene

C2v

A1 , A2 , B1 , B2

(11,10)

NEVPT2

Imidazole

Cs

A , A

(11,10)

NEVPT2

Benzene

D2h

B2g , B2u

(11,12)

NEVPT2

Phenol

Cs

A , A

(11,12)

NEVPT2

heteroatoms (see Table I). Three molecules have a single double or triple bond (ethylene, acetylene, and cyclopropane), six molecules have two double bonds and represent different types of aromaticity: cyclobutadiene is antiaromatic, cyclopentadiene is nonaromatic, and the five-membered heterocyclic compounds are aromatic (furane, pyrrole, thiophene, and imidazole). The high-level ab initio methods used limit the largest systems that can be studied to molecules with three conjugated double bonds, benzene and phenol. The structures of each monomer were optimized in their neutral charge states using standard DFT calculations with the BP86 functional,57–61 def2-TZVP basis set,62–64 and resolution of identity approximation65, 66 as implemented in the Turbomole package.67 Individual SCF energies were converged to 10−7 a.u. while gradient convergency criteria were set to 10−4 a.u. An enlarged DFT integration grid (m4) was used. Each structure was subsequently submitted to frequency calculations at the same level of theory to ensure that all normal modes were positive. In order to obtain the dimer of a given molecule, the optimized monomer structure was replicated along the axis perpendicular to the plane intercepting all heavy atoms of a single unit (Figure 1). Hence, due to sym-

metry, the initial and final states for electron transfer from one monomer to the other are the same.

III. THEORY A. Reference calculations

For a simple 2-state donor-acceptor system, the adiabatic ground (E1 ) and first excited potential energies (E2 ) are related to the diabatic (charge-localized) potential energies (Ea , Eb ) by E2,1 =

1 Ea + Eb ± (Ea − Eb )2 + 4|Hab |2 . 2

(3)

The adiabatic states, which diagonalize the electronic Hamiltonian, are well defined and can be obtained with standard electronic structure methods. Hence, for Ea = Eb , the first excitation energy between the adiabatic states, E12 = E2 − E1 , is exactly equal to twice the electronic coupling matrix element Hab between the diabatic states, 2|Hab | = E12 .

(4)


104105-4

Kubas et al.

We have chosen to include only homo-dimers (with identical monomer geometries) in our database because for these symmetric donor-acceptor pairs Ea = Eb , i.e., the reference electronic coupling matrix element does not depend on the diabatization procedure, but is uniquely defined by the excitation energy between the adiabats, Eq. (4). For unsymmetric cases where Ea = Eb , the GMH method10–17 gives access to the diabatic states and coupling between them via the transformation of an adiabatic dipole moment matrix. It is assumed that the diabatic states localized at different sites possess zero off-diagonal dipole moment matrix elements. Furthermore, only the vector components of each dipole matrix element in the direction defined by the difference between adiabatic dipole moments of the ground and excited states are considered. For a two-state problem, the 2 × 2 dipole moment matrix μDA is diagonalized and the same transformation is applied to the adiabatic (diagonal) Hamiltonian matrix. The symmetric matrix obtained has offdiagonal elements |Hab | referred to as electronic coupling matrix elements. In this way |Hab | can be obtained from a simple formula: |μ12 |E12 , (5) |Hab | = (μ11 − μ22 )2 + 4(μ12 )2 where E12 is the energy difference between two adiabatic states 1 and 2, μ11 and μ22 denote the adiabatic state dipole moments and μ12 is the transition dipole moment between these states. Equation (5) reduces to Eq. (4) in case of μ11 = μ22 (symmetric donor and acceptor system). The GMH approach can be used in combination with any ab initio method that gives access to excited states. Wavefunction-based approaches offer a ladder of methods of increasing accuracy for excited state calculations. The solution of the full, non-relativistic, and time-independent Schrödinger equation is only known for very small chemical systems by means of the full configuration interaction (FCI) method within a given basis set. Unfortunately, the number of determinants in the FCI wavefunction increases exponentially with the system size and more approximate methods have to be applied. In addition, special attention has to be paid on a balanced description of static and dynamic correlation: the former is connected to deficiencies of a single-determinantal Hartree-Fock (HF) reference and the latter is a consequence of the mean-field approximation in the HF theory. Thus, we selected three methods of decreasing computational complexity that can serve as a reference for small- to mediumsize molecules: (i) multireference configuration interaction singles and doubles method68 with a posteriori Davidson correction as generalized from the single-reference case by Peyerimhoff 69, 70 (MRCI+Q), (ii) perturbation approach with multideterminantal reference – n-electron valence state perturbation theory71–73 (NEVPT2), and (iii) HF-based secondorder approximate coupled-cluster method without74 and with75 spin-component scaling (CC2 and SCS-CC2, respectively). In each case the reference wavefunction for an openshell cation was used, CASSCF or UHF for (i)–(ii) or (iii), respectively. All three methods were used in case of ethylene, acetylene, cyclopropene, cyclobutadiene, cyclopentadiene, furane, and pyrrole. Moreover, it was possible to

J. Chem. Phys. 140, 104105 (2014)

perform NEVPT2 and CC2 calculations for thiophene, imidazole, benzene, and phenol. It is worth noting that the NEVPT2 method was recently benchmarked by Schapiro, Sivalingam, and Neese76 for low-lying vertical excitation energies of small organic molecules. It was found that the method performs very well and gives results similar to the more established CASPT2 approach. The CC2 method was used also for the larger acenes, naphthalene to pentacene. B. Constrained DFT

In conventional DFT calculations the energy functional E[ρ] is minimized with respect to an electron density ρ. This gives the adiabatic ground state E1 . The excited state E2 could in principle be calculated with TDDFT. However, due to the electron delocalisation error of common exchange correlation functionals the adiabatic states obtained from DFT and TDDFT are bound to be inaccurate. This problem is mitigated in CDFT calculations, where charge-localized diabatic states are obtained by imposing a charge constraint during minimization of the functional.27, 29, 77 In what follows we review the basic principles of the CDFT method as originally proposed by Wu and Van Voorhis.78 While this will aid the discussion on methodology, we refer here to the literature for a more detailed treatment.26, 30 The central ingredient of a CDFT calculation is the charge constraint Nc , defined as an integral of the electron density ρ(r) multiplied with a weight function w(r), Nc = ∫ w(r)ρ(r)dr.

(6)

The weight function w selects the spatial region subjected to the charge constraint and it defines the type of charge obtained from real-space integration. In the CPMD implementation of CDFT29, 30, 79 w is chosen so that Eq. (6) gives Hirshfeld charges.80 For instance, if one is interested in describing electron transfer from a donor (D) to an acceptor (A) and one would like to constrain the Hirshfeld charge difference between these two fragments, w takes the form, i∈D ρi (r − Ri ) − i∈A ρi (r − Ri ) , (7) w(r) = N i=1 ρi (r − Ri ) where the sums in the numerator run over all atoms of D and A fragments, respectively, whereas the sum in the denominator runs over all N atoms of the system and ρ i (r − Ri ) is the unperturbed (or pro-molecular) electron density of atom i. Minimization of E[ρ] under the constraint Eq. (6) can be carried out efficiently using the method of undetermined Lagrange multipliers. According to the formulation of Wu and Van Voorhis,78 a new functional W [ρ, V ] is defined, (8) W [ρ, V ] = E[ρ] + V w(r)ρ(r)dr − Nc , and the energy of the selected charge constraint state obtained by making W stationary with respect to ρ and the Lagrange multiplier V , E = minρ maxV [W [ρ, V ]].

(9)

The function V w(r) can be interpreted as the external potential that needs to be added to the Kohn-Sham equations to


104105-5

Kubas et al.

J. Chem. Phys. 140, 104105 (2014)

enforce the constraint Eq. (6). V is the height of this potential since w(r) ≤ 1∀r. Note that according to Eq. (7), the sign of the potential changes midway between donor and acceptor. The initial and final charge localized (diabatic) states (A) and (B) of the ET reaction under consideration can now be obtained by solving Eq. (9) for values Nc = Nc,A and Nc = Nc,B that represent the charge difference of donor and acceptor in the respective state (for the simple ET reactions considered here Nc,A = −Nc,B = 1 corresponding to full localization of an electron on the donor in state A or on the acceptor in state B). This gives two diabatic state Kohn-Sham (KS) determinants, ψ A and ψ B , and energies EA and EB , respectively, minimizing W for V = VA and V = VB , respectively. To proceed with the calculation of the electronic coupling matrix elements Hab , we note that in ET theory Hab is defined as the off-diagonal element of the exact electronic Hamiltonian H Eq. (1). In CDFT H is approximated by the state dependent Kohn-Sham Hamiltonian HKS and the two diabatic wavefunctions are approximated by the Kohn-Sham determinants obtained from the charge constrained minimization as explained above, ψ A and ψ B . The two states obtained from CDFT are in general non-orthogonal, which is why we have distinguished them from the orthogonal states in Eq. (1) by capital subscripts. Importantly, within the above approximations, it is possible to rigorously calculate the electronic coupling between the two KS determinants, HAB = ψA |HBKS |ψB = FB SAB − VB WAB ,

(10)

KS |ψA = FA SBA − VA WBA , HBA = ψB |HA

(11)

KS where HA = HKS [ρA ], HBKS = HKS [ρB ], ρ A , ρ B are the electron densities corresponding to ψ A and ψ B , FA and FB are the potential energies including the interaction with the exn KS ternal potential, F A = ψA |HA + VA i=1 w(ri )|ψA , FB = ψB |HBKS + VB ni=1 w(r the i )|ψB , SAB = ψA |ψB overlap and WAB = ψA | ni=1 w(ri )|ψB . Denoting the KS |ψA = EA and HBB diagonal elements by HAA = ψA |HA KS = ψB |HB |ψB = EB , we define the matrix,

H H AA AB H = , (12) HBA HBB

but note that this is not the final matrix representation of the 2 × 2 Hamiltonian as ψ A and ψ B are not orthogonal. Here we orthogonalize ψ A and ψ B by transforming them to the eigenstates ψ a and ψ b of the 2 × 2 matrix representation of n i=1 w(ri ), W , W V = SV L, where V is the 2 × 2 matrix of generalized eigenstates and L is the diagonal matrix of generalized eigenvalues. The final diabatic Hamiltonian matrix H in the orthogonal ψ a and ψ b basis is obtained from the similarity transformation, H = V † H V .

(13)

The final electronic coupling matrix element is the offdiagonal element of H, Hab = [H]ab .

(14)

We note that for the homo-dimer systems considered in this work, H is Hermitian. This is not necessarily true for non-

symmetric systems (although usually the deviation from nonHermiticity is small), in which case the coupling matrix element is taken as the average over the two off-diagonal elements. C. Fragment orbital DFT

The FODFT method is conceptually simpler than CDFT. The main difference is that in FODFT the diabatic state wavefunctions are constructed from orbitals of the isolated donor and acceptor fragments, i.e., unlike in CDFT, the interaction between donor and acceptor fragments is neglected. In the following, we review the FODFT method for electron transfer between a donor D and an acceptor A in case of cationic complexes as implemented in the CPMD program, and we refer to Ref. 43 for the case of anionic systems. For simplicity we assume that D and A are neutral and have the same number of electrons, N, and that the electron transfer is described by D+ + A → D + A+ , i.e., in the initial state D has N − 1 electrons and A has N electrons, and vice versa for the final state. The corresponding initial and final state diabatic wavefunctons are again described by a single determinant of 2N − 1 spin-orbitals φ. 1 det φa1 . . . φa2N−1 , (2N − 1)!

(15)

1 ψb = √ det φb1 . . . φb2N−1 . (2N − 1)!

(16)

ψa = √

As mentioned above, in FODFT the determinants ψ a and ψ b are constructed from the orbitals of the neutral, isolated (noninteracting) D and A fragments, ψa ≈ ψaDA = √

1 1 N−1 1 . . . φD φA . . . φAN , det φD (2N − 1)! (17)

1 1 N 1 . . . φD φA . . . φAN−1 . det φD (2N − 1)! (18) This requires two Kohn-Sham (KS) DFT calculations, one for isolated D and one for isolated A in the respective geometries of the fragments in the DA complex, giving two sets of KS N 1 . . . φD ) and (φA1 . . . φAN ), respectively. Hence, in orbitals (φD order to obtain states that resemble as closely as possible the states in reference and CDFT calculations, we biorthogonalize the two sets of primed orbitals to two sets of unprimed N 1 . . . φD ) and (φA1 . . . φAN ) according to the method orbitals (φD 81 of Löwdin, which are then combined to ψaDA and ψbDA according to Eq. (17) and (18). Choosing an orthogonal representation has the additional advantage that the corresponding Hamiltonian matrix is Hermitian, which is not necessarily the case in a non-orthogonal representation.3, 44 To obtain the coupling matrix element, we build the KS-Hamiltonian of the combined system. Since the KS-Hamiltonian is state dependent, we can againform two Hamiltonians, 2N−1 KS KS KS KS HbKS = 2N−1 i=1 hb,i and Ha = i=1 ha,i , where hb,i and KS ha,i are the one-particle KS-Hamiltonians constructed from the 2N − 1 orbitals that form ψbDA and ψaDA , respectively. The off-diagonal elements of the combined Hamiltonian are ψb ≈ ψbDA = √


104105-6

Kubas et al.

J. Chem. Phys. 140, 104105 (2014)

then given by Eq. (22) below, and similarly for Hba . Hab = ψa |H|ψb

(19)

≈ ψa HbKS ψb

(20)

≈ ψaDA |HbKS |ψbDA

(21)

N = φAN |hKS b |φD .

(22)

Hence, in the FODFT approach two strong assumptions are made as was noted previously.43 First, starting from the exact expression Eq. (19), the exact Hamiltonian is replaced by the KS Hamiltonian, which is dependent on the 2N − 1 orbitals (thus the subscript b or a), leading to Eq. (20). Second, in Eq. (21) the wavefunctions ψ a and ψ b are replaced by the wavefunctions composed of the fragment orbitals, ψaDA and ψbDA . The equality Eq. (22) is exact because ψaDA , ψbDA differ only by the occupation of the HOMO orbital on each fragment (see Eqs. (17) and (18), the N − 1 orbitals below the HOMO of each fragment do not relax upon electron transfer). Therefore, if assumptions Eqs. (20) and (21) are made, the electronic coupling matrix element is simply a matrix element over orbitals instead of determinants (specifically, HOMO orbitals of neutral donor and acceptor for electron transfer in cationic systems considered in this study). This method is also denoted FODFT(2N−1) and in all subsequent sections FODFT simply refers to FODFT(2N−1). In addition to H OMO ) FODFT, where the occupation number of φAH OMO (φD is set equal to zero when constructing the Hamiltonians hKS b (hKS a ), we have also carried out calculations where the occupation number of this orbital remain unchanged, i.e., where KS hKS b and ha are constructed from 2N rather than 2N−1 orbitals. This method will be referred to as FODFT(2N). We further note that for ET between the homo-dimers studied here Hab and Hba are identical. In case of asymmetric systems, where small differences can occur, the average of the two off-diagonal elements is taken.3, 26, 30 Coupling matrix elements were also calculated with the ADF program.82 In this particular implementation of FODFT21 2N electrons are used for the construction of the Hamiltonian. We note that in the standard ADF output the off-diagonal elements of the Hamiltonian are represented in a non-orthogonal basis. To be consistent with the CPMD FODFT methods, where orthogonal orbitals are used, the 2 × 2 HOMO subspace of the obtained Hamiltonian was subsequently transformed to an orthogonal HOMO subspace according to the method of Löwdin. Hab is obtained as offdiagonal element in the orthogonal basis. This method is referred to as ADF FODFT(2N). Particular attention needs to be paid in molecules where HOMO has a symmetry equivalent orbital like in acetylene or benzene molecules. In such a case of degenerated HOMOs there exist in general several charge transfer integrals between two fragments. For two degenerated orbitals {φna , φnb } on each of two fragments n the

Hamilton matrix can be written as ⎛ 0 Haa 1 ⎜ 0 Hba 1 H=⎜ ⎝ Haa Hba 2 Hab Hbb 0

⎞ Hab Hbb ⎟ ⎟, 0 ⎠ 2

(23)

with Hab = φ1a |H|φ2b . However, there is a unitary transformation U , ⎛ ⎞ cos(α) sin(α) 0 0 ⎜ − sin(α) cos(α) 0 0 ⎟ ⎟, U =⎜ (24) ⎝ 0 0 cos(β) sin(β) ⎠ 0 0 − sin(β) cos(β) where α and β are the rotation angles for orbitals on fragment 1 and 2, so that ⎛ ⎞ 1 0 Haa Hab ⎜ 0 ⎟ 1 Hba Hbb ⎜ ⎟ H = U −1 HU = ⎜ ⎟ 2 0 ⎠ ⎝ Haa Hba Hab Hbb 0 2 ⎛ ⎞ 1 0 Haa 0 ⎜ 0 ⎟ 1 0 Hbb ⎟ ! ⎜ ⎜ ⎟ . (25) = ⎝H 0 0 ⎠ 2 aa 0 Hbb 0 2 The requirement Hab = Hba = 0 defines α and β and thus leads to a new set of orbitals where these couplings are zero due to symmetry reasons (e.g., one set σ - and one set π -interacting orbitals). In this symmetry adapted basis the orthogonalized Hamiltonian can be written as

Hortho = U −1 S−1/2 U U −1 HU U −1 S−1/2 U = U −1 S−1/2 HS−1/2 U ,

(26)

defining a unique set of couplings. D. Density functional tight binding

DFTB denotes a class of computational models derived from DFT using a Taylor series expansion of the total energy around a reference density which is written as superposition of neutral atomic densities. In the lowest order (DFTB1)83, 84 no deviation from the reference density is considered in the Hamiltonian which leads to a non-self-consistent description. The 2nd (DFTB2)85 and 3rd (DFTB3)86 order terms in the energy expansion correspond to a self-consistent representation, where the deviation of the ground state density from the reference density is represented by charge monopoles only. The DFTB method is highly efficient both due to the application of a minimal basis and a two-center approximation, which makes the tabulation of overlap and Hamilton matrix elements possible. A review of these methods can be found in Refs. 87 and 88. As mentioned above one reason for the efficiency of DFTB is founded in the application of a minimal basis set. This causes a gain in computational speed but also makes it difficult to get the correct description for two interacting atoms over a wide range of distances, since an approximative


104105-7

Kubas et al.

J. Chem. Phys. 140, 104105 (2014)

basis function can only describe a certain part of the radial atomic orbital correctly and further, e.g., diffuse, basis functions are needed to get the tail of the orbital right. The standard parameterization of DFTB focuses on the correct description of bonded atoms and therefore uses confined basis functions, which are obtained by solving modified atomic Kohn-Sham equations, 2 r 1 2 (27) φμ = εμ φμ , − ∇ + veff (ρatom ) + 2 r0 where r0 is the confinement radius and φ μ are the resulting confined atomic basis functions.83, 84 The shape of φ μ is improved by this procedure in the region of typical bond distances. Unfortunately, the tail of the basis function gets suppressed and therefore it is only a poor representation of the atomic orbital at distances beyond the bond length. However, this part of the atomic orbital is crucial for the calculation of the electronic coupling, which is therefore severely underestimated, if standard DFTB parameter is used for coupling calculations.47 With a reduced confinement during the parametrization process of FODFTB, an improved description of the atomic basis at larger distances is achieved. With this procedure one can mimic the features of a double zeta basis with the computational cost of a minimal basis. The DFTB Hamiltonian can then be used for the calculation of Hab either in charge constrained (CDFTB)89–91 or fragment-orbital calculations (FODFTB).47–53 The latter approach is used in this work. In order to calculate Hab , first the Hamilton matrix is constructed in a atomic orbital basis and then transformed to a fragment-orbital basis according to HFO = C T HAO C,

(28)

where C is a block-diagonal matrix with blocks consisting of the eigenvectors of donor and acceptor respectively. Subsequently, the matrix HFO is reduced to only HOMOs of donor and acceptor, and orthogonalized according to the method of Löwdin to give the matrix element,

H OMO

. (29) Hab = Hba = φAH OMO hKS 0 φD We note that the FODFTB method uses the same number of electrons (2N) as the FODFT(2N) and ADF FODFT(2N) methods to construct the Hamiltonian. IV. DETAILS OF CALCULATIONS A. Reference calculations

In MRCI+Q and NEVPT2 calculations a state-average complete active space self-consistent field (CASSCF) wavefunction of a cation was chosen as a reference using aug-ccpVTZ basis set for heavy atoms and the smaller cc-pVDZ for hydrogens (see a comparison of various basis sets further below).92, 93 All π and π ∗ orbitals were selected to enter the active space with corresponding electrons (see Table I for individual active space sizes) and averaging was done over ground and first excited doublet states with equal weights if not stated otherwise. For the thiophene and acetylene dimers averaging was done over the four lowest doublet states because of small energy gaps between these states. The

CASSCF energy was converged to 10−7 a.u. while the energy gradient convergence criterion was tightened to 10−5 a.u. Unfortunately, the MRCI+Q expansion explodes rapidly with the active space size (number of reference functions). Therefore, additional truncations are introduced by excluding configurations with small CASSCF weights. In the present calculations the default value of Tpre = 10−4 a.u. is small enough and does not need to be changed. The frozen core approximation was applied for 1s electrons of carbon, nitrogen and oxygen atoms. The number of double excitations is further reduced by a perturbative approximation of the interaction of doubly excited wavefunctions with the reference and subsequent elimination of weak perturbers.94 This is controlled by the parameter Tsel and it was found that with Tsel = 10−7 a.u. more than 99% of excitation energy is recovered, compared to a value of 10−10 a.u. In case of NEVPT2 calculations, only configurations with a weight equal to or larger than 10−10 were chosen as a reference. Moreover, only CI vectors with a weight larger than 10−14 entered the third order density matrix. Orbitals were requested to be canonical for each state. Because of software limitations, the frozen core approximation was not applied. All CC2 calculations presented in this study were carried out with the same basis set as the multireference calculations and were based on an unrestricted Hartree-Fock (UHF) reference, that was converged up to 10−7 a.u. with the one-electron density convergence threshold set to 10−7 . In analogy to the spin-component scaled MP2 (SCS-MP2) method, the SCSCC2 model could be derived from the original CC2 method by applying scaling factors to same and opposite spin components of the elements of a Jacobi matrix needed to obtain vertical excitation energies. Here we should only note that we used the original parameters of Grimme (cos = 6/5 and css = 1/3) although other choices are possible.95 In all calculations we took the advantage of the resolution of identity (RI) approximation65, 66 with an auxiliary basis set fitted to the atomic basis set by Weigend et al.96 Moreover, the use of symmetry led to significant computational savings in case of MRCI+Q and CC2 methods, while it gave better control over calculations in case of CASSCF wavefunction optimisation and NEVPT2 treatment. Multireference calculations were carried out with the ORCA 2.9.1 program97 and (SCS)CC2 computations with the ricc2 module98 of the Turbomole suite of programs.67 In order to obtain reliable reference data from ab initio calculations, one also has to consider the size of the oneelectron basis set. This becomes even more crucial when observables calculated with atom-centred basis sets (e.g., Gaussians) are to be compared with those obtained with plane-waves. Here we considered a systematic series of Dunning’s correlation-consistent basis sets (denoted as ccpVXZ)92, 93 with various levels of augmentation with diffuse functions (original full augmentation, denoted as augcc-pVXZ, and minimally augmented as proposed recently by Papajak et al.,99 denoted as maug-cc-pVXZ). Basis sets were retrieved from the EMSL Basis Set Exchange Database.100, 101 We chose the ethylene dimer cation in parallel orientation separated by 4.0 Å to test convergence of the MRCI+Q excitation energy with respect to the basis set used. Four


Kubas et al.

104105-8

J. Chem. Phys. 140, 104105 (2014)

TABLE II. Excitation energies E of an ethylene dimer cation in parallel orientation at the intermolecular separation of 4.0 Å calculated at CASSCF and MRCI+Q level of theory with various basis set combinations. The MRCI+Q error is reported with respect to the aug-cc-pVQZ/cc-pVQZ basis set. Basis C cc-pVDZ cc-pVTZ cc-pVQZ maug-cc-pVDZ maug-cc-pVTZ aug-cc-pVDZ aug-cc-pVTZ aug-cc-pVTZ aug-cc-pVQZ

H cc-pVDZ cc-pVTZ cc-pVQZ cc-pVDZ cc-pVTZ cc-pVDZ cc-pVTZ cc-pVDZ cc-pVQZ

E(CASSCF) E(MRCI+Q) MRCI+Q (eV) (eV) error (%) 0.561 0.576 0.579 0.583 0.581 0.585 0.581 0.580 0.581

0.518 0.527 0.533 0.534 0.535 0.548 0.544 0.542 0.540

4.3 2.6 1.6 1.3 1.2 1.1 0.4 0.2 –

π -type orbitals were selected to enter the active space together with three electrons [CAS(3,4)]. The state-averaging procedure was employed in order to simultaneously describe the doublet ground state and the first excited doublet state. The largest basis set used was a combination of aug-cc-pVQZ for C atoms and cc-pVQZ for H atoms. The excitation energy of 0.540 eV is treated as a benchmark and errors reported in Table II are with respect to this setup. The CASSCF energies are already well converged with the maug-cc-pVTZ basis set. Moreover, we found that such minimal augmentation can be used to decrease basis set error in MRCI+Q calculations to about 1%. Full augmentation of the basis set for carbon atom, however, is needed to get an error smaller than 0.5%. As expected, the small basis set can be used for H atoms as the excitation process does not involve any hydrogen orbitals. Thus, a combination of the aug-ccpVTZ basis set for heavy atoms and the cc-pVDZ basis for hydrogen atoms was found to give sufficiently accurate excitation energies for electron transfer reactions. This finding is in line with studies of Pieniazek et al.,35 where the augcc-pVTZ basis set was found to provide almost converged electronic coupling matrix elements with FCI and equationof-motion coupled cluster singles and doubles (EOM-CCSD) methods. To test the proposed set-up we calculated electronic cou+ plings, Hab , for the diatomic cations He+ 2 and Zn2 at various internuclear distances and compare the results obtained with available literature data.10, 31 For the helium dimer cation CASSCF(3,2)/MRCI+Q couplings are overestimated by 0.1–2.7 meV compared to FCI data and the exponential coupling decay constant β is underestimated by 1%. This suggests that the minimal active space is insufficient for very accurate calculations. We performed the same calculations with a larger active space where two 2s orbitals were taken into account as well. The couplings obtained are now much closer to the FCI values and the error in β is reduced to 0.01%. In case of the zinc dimer cation, we tested various active spaces and we found that the best description is obtained for an active space composed of two 4s and six 4p orbitals together with three electrons [CASSCF(3,8)]. Our CASSCF results

are close to values obtained by Cave and Newton10 in stateaverage CASSCF(3,8) calculations but note the latter calculations employed a smaller basis set and averaging was done over four states in contrast to only two lowest lying states in our case. Moreover, we did not impose any symmetry restrictions in these calculations. Interestingly, excitations and, in consequence, couplings become significantly lower with inclusion of dynamic correlation effects at the MRCI+Q level of theory. Thus, a larger decay constant β = 2.60 Å−1 was obtained with our best setup and this can be considered as a new reference. B. CDFT

In order to model electron transfer between the homodimers of organic molecules, we have used the difference in the Hirshfeld charge between the donor and acceptor monomer and included all atoms of the donor and acceptor monomers in the weight function defined in Eq. (7). The charge constraint was set to Nc,A = 1 for the initial state where the hole is localized on the donor monomer D, and to Nc,B = −1 for the final state, where the hole is localized in the acceptor monomer A. Since donor and acceptor are identical, the final Lagrange multipliers obtained for states A and B differ only in the sign, VA = −VB . CDFT calculations were carried out in the gas phase using the Perdew-Burke-Ernzerhof (PBE) functional. Modified PBE functionals were tested as well, where a given percentage of the GGA exchange was replaced by Hartree-Fock exchange (HFX). We denote these functionals with “/X” suffix, where the numeral X refers to the percentage of HFX. The dimer was placed in the centre of a rectangular box such that the smallest distance to the box edges was 5 Å. In this way, one ensures the wavefunction of the dimer to be negligibly small at the box edges (increasing the distance to 6 Å had an effect on Hab of less than 0.1 meV). The reciprocal space plane wave cutoff for KS-orbitals was set to 80 Ry, which was found to be sufficiently high. Higher cut-offs gave rise to changes in Hab of less than 1 meV. Additionally, in CDFT calculations with exact Hartree-Fock exchange the cut-off energy of the electron density was decreased by a factor of two to allow for a more efficient calculation of the electronic couplings. This could be done at virtually no loss of accuracy. Only valance electrons were considered explicitly, whereas the core electrons were described with Troullier-Martins pseudopotentials.103 The wavefunction was optimized using the method of the preconditioned conjugate gradients (PCG), where convergence was achieved if the largest element of the wavefunction gradient was below 10−6 a.u. The charge constraint was converged using the Newton scheme together with a convergence criterion of 5 × 10−5 a.u. All CDFT calculations were carried out with the CPMD program.79 To allow a direct comparison with reference data for each molecule, the electron transfer process must involve the same pairs of orbitals. As the shape of spin density isosurface resembles essentially the shape of molecular orbitals being coupled in calculations, we always inspected the superimposed spin densities of donor and acceptor fragments and a ground state spin density from reference calculations. When both


104105-9

Kubas et al.

J. Chem. Phys. 140, 104105 (2014)

TABLE III. Hab values (meV) for He+ 2 and the decay rate constant β (1/Å) calculated with the CASSCF/MRCI+Q approach as compared to full CI (FCI) numbers with aug-cc-pVTZ basis set.35 Two 1s orbitals [CASSCF(3,2)] or two 1s and two 2s orbitals [CASSCF(3,4)] along with three electrons were used to construct the reference wavefunction in MRCI+Q calculations. For the sake of completeness, FODFT(2N), FODFT(2N−1), and CDFT/030 are also reported. β values are obtained from the exponential fit to the data points given. R (Å) 2.5 3.0 4.0 5.0 β(Å−1 ) a b

FODFT(2N) a

FODFT(2N−1)

CDFT/0a

CASSCF(3,2)/MRCI+Q

CASSCF(3,4)/MRCI+Q

FCIb

170.76 56.25 6.14 0.68

153.79 45.32 3.57 0.23

207.23 67.37 6.92 0.57

198.61 62.37 5.88 0.51

195.84 60.73 5.51 0.44

195.93 60.81 5.52 0.44

4.42

5.20

4.71

4.78

4.88

4.82

Reference 30. Reference 35.

isosurfaces were qualitatively the same the calculations were accepted while when some discrepancies occurred the computations were restarted from other starting orbitals (e.g., from FODFT or other amount of HFX that was successful to give proper spin density). An example of such comparison is given in Fig. S1 in the Supporting Information (SI).56 Moreover, the value of integrated spin density in each calculations can be used as a diagnostic tool in CDFT calculations. In our case, all molecules are in a doublet spin state and should ideally show integrated absolute spin density of one. Any deviations indicate that some degree of symmetry breaking occurs, i.e., spatial part of alpha and beta orbital manifolds start to differ. Although special care was taken to couple always the same orbitals as in the reference calculation, results can still be affected by this unwanted feature to some extent. For the sake of completeness the integrated spin density values for all molecules in this study are presented in Table S1 in SI along with a deeper analysis of this rather technical problem.56

C. FODFT

FODFT and FODFT(2N) calculations were carried out for the HAB11 database with the CPMD program79 using the PBE functional and the same plane wave cutoff, box dimensions and pseudopotentials as used in CDFT calculations (see + Sec. IV B). For the FODFT calculations on He+ 2 and Zn2 , which are summarized in Tables III and IV, we have used the same system setup as in Ref. 30. ADF FODFT(2N) calculations on the database dimers were carried out with the ADF program package82 using the PBE functional and the TZ2P

all-electron basis set. The Hamiltonian obtained was transformed according to the method of Löwdin as explained in Sec. III C. The values obtained from different FODFT implementations are discussed in Sec. VI B. D. FODFTB

All calculations are carried out within the framework of FODFTB2 using our in-house program. Since most of the molecules in the regarded database are nonpolar, probably comparable couplings would have been obtained using FODFTB1 with even less computational effort. Calculations are carried out as explained in Sec. III D. In principle the full basis could have been orthogonalized, but the effect of additional orbitals on the HOMO coupling was found to be negligible. As mentioned in Sec. III D an effective double zeta description is used, which is established by using the standard mio-parameters85, 104 for the diagonal blocks of HAO which represent short-range intra-fragment interactions and a different set for the off-diagonal blocks of HAO which represent long-range inter-fragment interactions. For the latter parameter set the wave function confinement was reduced to 8 a.u. and the density confinement radius was extended to ∞ for all atoms. V. RESULTS A. High-level ab initio reference calculations on HAB11

MRCI+Q is the most accurate method used in this study. Unfortunately, the method scales unfavourably with respect to

TABLE IV. Hab values (meV) for Zn+ 2 and the decay rate constant β (1/Å) calculated with the CASSCF(3,8) wavefunction with two basis sets (aug-cc-pVTZ and Wachters+102 where “+” sign denotes some augmentation by diffuse functions performed by the authors10 ) and our best estimate CASSCF(3,8)/MRCI+Q//aug-cc-pVTZ. For the sake of completeness, FODFT(2N), FODFT(2N−1) and CDFT/030 are also reported. β values are obtained from the exponential fit to the data points given.

R (Å) 5.0 6.0 7.0 β (Å−1 ) a b

FODFT(2N)a 80.67 22.91 6.50 2.52

FODFT(2N−1) 65.10 13.82 2.31 3.34

CDFT/0a 170.32 43.63 11.30 2.71

CASSCF(3,8) aug-cc-pVTZ 173.71 56.02 15.49 2.42

CASSCF(3,8)b Wachters+ 197.56 58.78 16.95 2.46

CASSCF(3,8)/MRCI+Q aug-cc-pVTZ 149.3 43.74 11.13 2.60

Reference 30. Reference 10.


104105-10

Kubas et al.

J. Chem. Phys. 140, 104105 (2014)

TABLE V. Comparison of ab initio methods used in this study on dimers selected from the HAB11 database (Table I). MRCI+Q data serve as a reference: absolute errors of electronic coupling matrix elements Hab at various distances between planes of molecules (see Figure 1) and absolute errors of decay rate constant β for each method. Statistical evaluation of both quantities are given by means of mean unsigned error [MUE = ( n |ycalc − yref |)/n], mean relative signed error [MRSE = ( n ((ycalc − yref )/yref ))/n], mean relative unsigned error [MRUE = ( n (|ycalc − yref |/yref ))/n] and maximum unsigned error (MAX = max |ycalc − yref |). All values in meV. Dimer

NEVPT2

Ethylene

|Hab |

3.5 Å 4.0 Å 4.5 Å 5.0 Å

Acetylene

|Hab |

3.5 Å 4.0 Å 4.5 Å 5.0 Å

β

β Cyclopropene

|Hab |

3.5 Å 4.0 Å 4.5 Å 5.0 Å

− 31.8 − 18.5 − 8.4 − 2.4

Cyclobutadiene

|Hab |

3.5 Å 4.0 Å 4.5 Å 5.0 Å β

Cyclopentadiene

|Hab |

3.5 Å 4.0 Å 4.5 Å 5.0 Å

Furane

|Hab |

3.5 Å 4.0 Å 4.5 Å 5.0 Å

Pyrrole

|Hab |

3.5 Å 4.0 Å 4.5 Å 5.0 Å

β

β

β

19.5 33.3 46.4 56.7

− 0.04 − 20.6 − 12.2 − 7.5 − 5.0

− 0.75 − 9.0 7.1 20.4 30.4

− 0.06 − 50.8 − 30.9 − 17 − 7.8

− 0.59 14.1 29.3 42.5 51.1

0.03

β

CC2

− 27.8 − 14.9 − 5.2 0.7

− 0.76 33.4 49.5 62.7 72.2

− 0.10 − 25.5 − 11.2 − 2.3 3.8

− 1.04 18.9 36.4 48.6 57

− 0.18 − 27.6 − 14.1 − 4.4 2.9

− 1.01 25.1 41.2 54 62.9

− 0.15 − 31.8 − 18.5 − 8.4 − 2.4

− 0.98 19.5 33.3 46.4 56.7

− 0.04

− 0.75

|Hab |

MUE (meV) MRSE (%) MRUE (%) MAX (meV)

15.6 − 5.8 6.9 50.8

39.5 39.8 39.9 72.2

β

MUE (1/Å) MRSE (%) MRUE (%) MAX (1/Å)

0.08 − 2.1 2.8 0.18

0.86 − 30.0 30.0 1.04

the number of reference configurations and number of correlated electrons and orbitals. Thus, we tested three alternative approaches – NEVPT2, CC2, and SCS-CC2 – that can provide good excitation energies at reduced cost. For a subset of our database (see Table V) Hab values at various distances were calculated and the decay constant β was evaluated by a

SCS-CC2

MRCI+Q

− 32.6 − 22.7 − 13.3 − 6.1

519.2 270.8 137.6 68.5

0.04 − 57.6 − 45.4 − 35.4 − 28.0 0.72 − 39.7 − 26.2 − 14.4 −7 0.04 − 27.9 − 14.3 − 3.4 4.3 − 0.18 − 34.5 − 19.2 −9 − 2.2 − 0.04 − 30 − 16.3 − 5.8 1.3 − 0.12 − 32.6 − 22.7 − 13.3 − 6.1 0.04

2.70 460.7 231.8 114.8 56.6 2.80 536.6 254.0 118.4 54.0 3.06 462.7 239.1 121.7 62.2 2.67 465.8 234.4 114.3 53.4 2.89 440.3 214.9 101.8 46.0 3.01 456.3 228.6 111.3 52.2 2.89

20.5 − 10.1 10.8 57.6 0.17 2.8 6.1 0.70

linear fit to Eq. (2) with all above mentioned ab initio methods. Both NEVPT2 and SCS-CC2 approaches perform very well. The mean unsigned errors (MUE) are 15.6 meV and 20.5 meV, respectively, and the mean relative signed errors (MRSE) as well as mean relative unsigned errors (MRUE) indicate an underestimation by about 6%–10% relative to


104105-11

Kubas et al.

MRCI+Q. The exponential distance decay constants (β values) are even less sensitive to the level of theory used with MRUEs of 2.8% and 6.1% for NEVPT2 and SCS-CC2, respectively. In contrast, the unscaled CC2 method significantly overestimates couplings (MUE of 39.5 meV and MRUE of almost 40%) and gives too low decay constants (MUE larger than 0.8 Å−1 ). The reliability of the CC2 model for ground state energies could be measured by the so-called D1 diagnostic proposed by Janssen and Nielsen.105 Ideally D1(CC2) should be smaller than 0.03 but in our calculations it is rarely below 0.08. This clearly indicates high static correlation involvement that can be covered either by higher-order excitations or by more appropriate reference wavefunction, e.g., CASSCF. In the original SCS-CC2 paper75 authors argued that the method should be less sensitive (compared to CC2) to strong correlation effects caused by inclusion of single excitations in the reference due to damping of the exchange correlation effects. To the best of our knowledge, the SCS-CC2 approach was tested only for π → π ∗ and n → π ∗ excitations within neutral species.75, 106 No significant improvement was observed for vertical excitation energies. We noted, however, that our unsaturated systems are positively charged and the excitations are of π → π type. It is clear from Table V that spin-component scaling uniformly improves over the unscaled counterpart for all molecules except acetylene. This is the only molecule with a triple bond in our test set. Thus, π -conjugated systems with single and double bonds only can be well described by the SCS-CC2 method but the scope could be possibly broadened by different scaling parameters. This phenomenon deserves further investigation that is beyond the scope of this paper. Here, we should only note that the NEVPT2 method will be used as a reference in the database when MRCI+Q calculations are computationally too expensive. Since SCS-CC2 performs almost as good as the NEVPT2 approach, but at a fraction of the computational cost of the multireference calculations, it will be used as a reference for calculations on a series of acenes, reported later in this paper. B. CDFT, FODFT, and FODFTB calculations on HAB11

The high quality benchmark data obtained in Sec. V A can now be compared with the results of DFT-based approaches. The numerical results are compiled in Table VI for the HAB11 database. Interestingly, the reference values for a given separation distance are very similar, irrespective of the number of π -bonds and type of heteroatom, about 450, 200, 100, and 50 meV for 3.5, 4.0, 4.5, and 5.0 Å, respectively. We find that the CDFT/PBE method (CDFT/0) overestimates electronic couplings, MRSE = +38.7% and MRUE = 38.7% (averaged over all dimers and distances). Inclusion of Hartree-Fock exchange (HFX) leads to a systematic decrease in electronic couplings, that roughly scales with the fraction of HFX (see Figure 2). Replacing 25% of the GGA exchange in the PBE functional by HFX (designated here as CDFT/25), the MRUE decreases to 13.8%. Further increase to 50% HFX (CDFT/50) gives the best results, an MRUE of 5.3%, whereas at 100% HFX (CDFT/100) the couplings are

J. Chem. Phys. 140, 104105 (2014)

underestimated by about 19%. This trend is a consequence of the reduction of electron delocalization with increasing percentage of HFX as discussed in more detail in Sec. VI A. We would like to stress, however, that the dependence of the CDFT method on the fraction of HFX is much weaker than for standard DFT calculations. CDFT at the GGA level still gives reasonable results, whereas the error (underestimation) for the first excitation energy with TDDFT/GGA is expected to be significantly larger.107 Interestingly, the relative errors are only very weakly dependent on the donor-acceptor separation distance for any of the functionals considered. The general trend is that the errors slightly decrease with increasing separation distance (Figure 2), indicative of an improved performance of CDFT at larger distances. This is in accord with expectations, as the influence of the choice of the weight function used for charge localization becomes smaller and eventually vanishes for very large separations. (Note that DFT/TDDFT will generally not yield improved results for larger distances.) As a consequence of the weak distance dependence of the errors, the β values obtained with CDFT are very close to the reference values. In general, they are slightly overestimated with MRUEs of 2.3%, 4.3%, and 4.9% for CDFT/0, CDFT/50, and CDFT/100, respectively (averaged over all dimers). The weak but systematic increase in β values with increasing fraction of HFX is a consequence of the increasing compactness of the orbital carrying the excess charge. In the FODFT method the delocalization error of the excess charge is circumvented altogether by constructing the wavefunction from orbitals of the isolated monomers only. Because there is no direct interaction between the monomers, one can expect that this method yields smaller couplings than CDFT. This is indeed the case (see Table VI and Figure 2). The errors for FODFT are comparable in magnitude to CDFT/0 but of opposite sign, MRUE = 37.6% and MRSE = −37.3%. The MUE is slightly reduced by about 15 meV and the maximum error is reduced by 95 meV relative to CDFT/0. However, FODFT couplings drop too quickly with intermolecular distance, as reflected by an increased MRUE in β values of 13.0%, up from 2.3% for CDFT/0. The increasing relative underestimation of coupling values for larger distances is somewhat surprising and suggests that polarization effects are still important at larger distances. We expect that admixture of HFX will not further improve the FODFT results (more likely, results will worsen due to possible further contraction of orbitals), which is why we have not carried out FODFT calculations with HFX. FODFTB couplings and decay constants are similar to the results of FODFT calculations. This is not surprising as both methods rely on isolated fragment calculations. The MRUE in electronic couplings is 42.4%, only slightly higher than for FODFT, and the couplings drop again too quickly with intermolecular distance, as reflected by an overestimation of β values by 12.5%. C. Scaling factors

The results presented so far allow us to build a hierarchy of DFT-based methods for coupling matrix element


104105-12

Kubas et al.

J. Chem. Phys. 140, 104105 (2014)

TABLE VI. Absolute errors of electronic coupling matrix elements |Hab | for the HAB11 database (Table I). Values for electron transfer for the homodimer cations at various separation distances between monomers (Figure 1), and absolute errors of the decay constant β for each method are given. Statistical evaluation of both quantities are given by means of mean unsigned error [MUE = ( n |ycalc − yref |)/n], mean relative signed error [MRSE = ( n ((ycalc − yref )/yref ))/n], mean relative unsigned error [MRUE = ( n (|ycalc − yref |/yref ))/n] and maximum unsigned error (MAX = max |ycalc − yref |). CDFT/PBE%HFX Dimer Ethylene

Acetylene

Cyclopropene

Cyclobutadiene

Cyclopentadiene

0 |Hab |

|Hab |

|Hab |

|Hab |

|Hab |

Pyrrole

Thiophene

|Hab |

Hab

|Hab |

Benzene

|Hab |

|Hab |

60

75

100

FODFT

FODFTB

Ref.

102.3 44.0 20.7

41.2 9.0 1.1

37.7 7.7 0.6

− 12.7 − 20.7 − 14.8

− 28.6 − 29.6 − 19.6

− 49.4 − 41.4 − 26.1

− 151.5 − 100.9 − 61.6

− 175.5 − 99.3 − 56.5

519.2a 270.8a 137.6a

5.0 Å β 3.5 Å

9.9 0.06 103.0

− 0.3 0.11 40.1

− 0.5 0.10 − 1.6

− 8.5 0.14 − 14.1

− 11.1 0.16 − 30.1

− 14.6 0.19 − 51.0

− 35.9 0.53 − 143.8

− 32.9 0.32 − 169.2

68.5a 2.701a 460.7a

4.0 Å 4.5 Å

41.4 16.4

8.1 − 0.5

− 20.2 − 14.8

− 28.6 − 19.3

− 39.9 − 25.3

− 91.9 − 55.1

− 85.5 − 43.2

5.0 Å β 3.5 Å

5.8 0.14 279.5

− 2.5 0.17 144.2

− 7.9 0.19 49.1

− 9.7 0.21 19.7

− 11.9 0.22 − 17.6

− 15.1 0.26 − 66.4

− 32.3 0.63 − 117.8

− 21.7 0.18 − 201.1

56.6a 2.797a 536.6a

4.0 Å 4.5 Å 5.0 Å

115.6 46.4 19.9

48.5 16.3 6.3

5.5 − 2.3 − 1.8

− 7.2 − 7.9 − 4.2

− 23.5 − 15.0 − 7.4

− 44.6 − 24.2 − 11.6

− 74.8 − 42.5 − 22.3

− 121.8 − 71.2 − 41.7

254.0a 118.4a 54.0a

β 3.5 Å

0.14 183.6

0.17 28.2

0.16 − 66.8

0.16 − 93.7

0.15 − 127.1

0.15 − 169.8

0.38 − 139.4

0.71 − 182.6

3.061a 462.7a

4.0 Å 4.5 Å

87.7 42.8

4.2 − 1.6

− 44.0 − 25.9

− 57.7 − 32.7

− 74.6 − 41.3

− 96.4 − 52.5

− 90.2 − 53.8

− 103.9 − 56.7

239.1a 121.7a

5.0 Å β

25.7 − 0.01

− 4.3 0.17

− 18.5 0.25

− 22.5 0.28

− 27.4 0.33

− 33.9 0.42

− 31.7 0.47

3.5 Å

236.3

109.4

24.7

− 0.4

− 32.2

− 73.9

− 122.5

− 160.0

465.8a

4.0 Å 4.5 Å 5.0 Å

112.0 53.1 27.5

41.4 18.5 10.4

0.1 − 1.0 1.5

− 11.6 − 6.4 − 1.1

− 26.4 − 13.5 − 4.5

− 45.7 − 22.9 − 9.1

− 76.7 − 42.5 − 21.1

− 86.7 − 45.6 − 22.8

234.4a 114.3a 53.4a

0.00

0.04

− 13.6 − 11.4

0.03

0.03

3.5 Å 4.0 Å 4.5 Å

157.9 77.6 39.2

70.6 28.7 14.3

12.1 − 1.2 − 0.3

− 5.5 − 9.9 − 4.5

5.0 Å β

21.0 − 0.09

9.1 − 0.04

2.2 − 0.02

0.2 − 0.02

3.5 Å 4.0 Å

173.5 86.2

77.6 30.5

13.5 − 2.9

4.5 Å 5.0 Å β

44.2 23.5 − 0.07

14.9 9.2 − 0.01

− 1.9 1.0 0.01

3.5 Å 4.0 Å 4.5 Å

220.9 113.4 55.9

105.8 46.5 20.3

5.0 Å β Imidazole

50

3.5 Å 4.0 Å 4.5 Å

β Furane

25

3.5 Å 4.0 Å

23.0 0.06 178.5 83.5

− 27.7 − 20.8 − 9.8 − 2.4 − 0.01

0.02 − 56.8 − 34.9 − 16.8 − 5.8 0.00

− 29.4 − 24.4

− 60.0 − 39.6

− 6.6 − 1.3 0.02

− 12.5 − 4.2 0.02

− 20.2 − 8.1 0.04

27.9 6.2 0.6

4.9 − 5.1 − 4.8

− 24.3 − 19.2 − 11.8

− 61.9 − 37.3 − 20.6

− 3.5 0.16

− 6.1 0.16

89.6 35.1

31.8 6.5

15.0 − 1.5

5.8 0.14

− 5.5 − 12.5

0.02

− 9.4 0.17 − 6.3 − 11.7

− 13.7 0.18 − 33.5 − 24.6

0.26 − 124.7 − 73.3 − 39.0 − 18.5 0.24 − 127.6 − 78.6 − 43.5 − 21.9 0.29 − 107.8 − 64.8 − 37.7 − 24.0 0.40 − 100.9 − 63.4

− 29.4 0.40

0.31 − 170.1 − 97.5 − 53.8 − 28.4 0.43 − 184.2 − 85.4 − 40.4 − 23.2 0.03 − 134.1 − 75.2 − 42.9 − 26.3 0.47 − 135.3 − 72.1

231.8a 114.8a

62.2a 2.678a

2.886a 440.3a 214.9a 101.8a 46.0a 3.009a 456.3a 228.6a 111.3a 52.2a 2.890a 449.0b 218.9b 106.5b 54.4b 2.821b 411.6b 202.8b

4.5 Å 5.0 Å β

37.7 14.6 0.13

13.9 3.9 0.16

0.4 − 2.4 0.16

− 3.5 − 4.2 0.16

− 8.3 − 6.5 0.16

− 14.6 − 9.6 0.17

− 37.2 − 22.6 0.43

− 39.5 − 24.3 0.34

99.1b 49.7b 2.823b

3.5 Å 4.0 Å 4.5 Å 5.0 Å β

212.6 107.4 51.6 22.1 0.05

113.4 47.7 20.6 6.9 0.14

40.3 9.2 1.6 − 1.9 0.16

17.5 − 2.1 − 3.7 − 4.4 0.16

− 11.8 − 16.2 − 10.6 − 7.6 0.17

40.4 9.3 1.7 − 1.8 0.16

− 92.8 − 59.5 − 34.8 − 21.1 0.37

− 143.6 − 74.2 − 38.7 − 21.9 0.36

435.2b 214.3b 104.0b 51.7b 2.846b


104105-13

Kubas et al.

J. Chem. Phys. 140, 104105 (2014)

TABLE VI. (continued.) CDFT/PBE%HFX Dimer Phenol

|Hab |

β

a b

|Hab |

3.5 Å 4.0 Å 4.5 Å 5.0 Å β

MUE (meV) MRSE (%) MRUE (%) MAX (meV) MUE (1/Å) MRSE (%) MRUE (%) MAX (1/Å)

0

25

50

60

75

100

FODFT

FODFTB

Ref.

182.3 91.7 44.2 17.6 0.05

88.3 36.2 15.2 4.6 0.13

22.8 3.2 − 0.3 − 2.0 0.14

2.1 − 6.2 − 4.9 − 4.2 0.15

− 26.2 − 18.7 − 10.5 − 6.7 0.13

− 122.6 − 52.2 − 24.1 − 12.5 − 0.05

− 184.5 − 98.5 − 50.5 − 18.9 − 0.04

− 201.1 − 121.8 − 71.2 − 41.7 0.71

375.0b 179.6b 85.2b 41.3b 2.946b

83.0 38.7 38.7 279.5 0.07 1.5 2.3 0.14

32.8 13.8 13.8 144.2 0.12 3.8 3.9 0.17

11.8 − 0.8 5.3 66.8 0.13 4.4 4.3 0.25

11.7 − 6.0 7.1 93.7 0.14 4.7 4.5 0.17

21.1 − 11.7 12.2 127.1 0.14 5.0 4.7 0.17

37.2 − 18.2 19.4 169.8 0.15 5.0 4.9 0.42

69.4 − 37.3 37.6 184.5 0.37 12.8 13.0 0.63

82.4 − 42.9 42.4 201.1 0.33 11.7 12.5 0.71

MRCI+Q reference. NEVPT2 reference.

calculations. Four approaches of decreasing computational cost, CDFT/50, CDFT/0, FODFT, and FODFTB, could be applied to a wide range of molecular and extended systems ranging from small molecules to proteins and spanning several orders of magnitude in computational cost. In Figure 3 we show the coupling matrix elements obtained with each of these methods against the reference ab initio values. Particularly striking is the excellent linear correlation with R2 values of 0.9904, 0.9869, 0.9788, and 0.9853 for CDFT/50, CDFT/0, FODFT, and FODFTB, respectively. Moreover, the correlation is better for smaller couplings (larger distances).

FIG. 2. Dependence of the mean relative signed error (MRSE, top panel) and mean relative unsigned error (MRUE, bottom panel) on the intermolecular distance d (see Figure 1) for various methods. Percentage of Hartree-Fock exchange (HFX) in CDFT/X calculations is denoted by numerals in the legend (X = 0, 25, 50, 60, 75, 100).

Hence, our results suggest that the DFT-based results can be simply corrected by applying a method-specific uniform scaling factor. Here, the scaling factors are simply obtained as the inverse of the slope of the best-fit line to the data points plotted in Figure 3. The final values are 0.962, 0.721, 1.348, and 1.540 for CDFT/50, CDFT/0, FODFT, and FODFTB, respectively. The errors after scaling are presented in Table VII. The most impressive improvement is obtained for CDFT/0 where the MRUE decreases by an order of magnitude to 5.2%, followed by FODFTB and FODFT, where the error reduces by approximately a factor of three to 13.9% and two to 17.6%, respectively. The CDFT/50 scaling factor is very close to unity (0.962) and no statistical improvement is gained. This is mainly due to larger underestimation of couplings at short distances that is not compensated by smaller error at longer distances. Thus the unscaled method is recommended in general applications. The proposed scaling procedure should work well for systems where orbitals involved in charge transfer are mainly carbon-centred and possibly for other second row atoms

FIG. 3. Correlation between calculated CDFT/0, CDFT/50, FO-DFT, FODFTB |Hab | values and reference MRCI+Q/NEVPT2 data. The black dashed line represents perfect correlation (y-intercept was set to zero).


104105-14

Kubas et al.

J. Chem. Phys. 140, 104105 (2014)

TABLE VII. Scaling factors for electronic coupling matrix elements for CDFT and FODFT methods as well as for the FODFTB approach (standard error is given in brackets). Errors in |Hab | as obtained after scaling are also given.

Method

Scaling factor

MUE (meV)

MRSE (%)

MRUE (%)

MAX (meV)

0.962 (0.015) 0.721 (0.013) 1.348 (0.031) 1.540 (0.029)

11.7 10.6 24.7 20.0

− 6.6 − 1.8 − 17.3 − 12.7

7.1 5.2 17.6 13.9

81.7 71.1 142.1 69.7

TABLE VIII. The dependence of the electronic coupling matrix element |Hab | (meV) on the rotation angle α in the ethylene dimer (see Figure 4(a) for a definition) calculated with CDFT/50 and scaled CDFT/0, FO-DFT, and FODFTB approaches (note prefix “s”). Unscaled values are given in brackets and scaling factors are reported in Table VII. Reference MRCI+Q data are reported in the last column. α [deg]

CDFT/50 CDFT/0 FODFT FODFTB

involving 2p orbitals. Moreover, we find that the same scaling factor gives excellent results for FODFT couplings for the He+ 2 system, but a larger scaling factor would be needed for Zn+ 2 due to the higher polarizability of the 4s orbital (see Tables III and IV for unscaled values).

D. Performance on rotated structures

All dimers in the database are perfectly π -stacked. We found that economic methods such as CDFT/0, FODFT, or FODFTB can provide good to excellent Hab values when couplings are uniformly scaled by an empirical correction factor (see Sec. V C). When about 50% HFX is added to PBE in CDFT calculations, the method reproduces the reference data very well without scaling. We will now investigate whether these results carry over to geometries where the monomers are rotated with respect to each other so that the dimers are no longer perfectly π -stacked (see Figure 4). The calculations were performed without symmetry restrictions, although in systems shown in Figures 4(a)–4(c) the monomers are symmetry equivalent (dimers have virtually D2 , C2 , and Cs point groups, respectively). The first test case was the rotation of one ethylene molecule around the centre of mass at the intermolecular separation of 5 Å (see Figure 4(a)). Inspection of Table VIII shows clearly that the coupling matrix elements are hardly dependent on such a change of the molecular framework. The agreement with absolute reference couplings is excellent in case of CDFT/50 while it is reasonable for all scaled methods. The unscaled FODFT and FODFTB approaches have here MRUE of 53.0% and 48.9%, respectively. Significant improvement is observed after application of scaling factors –

0 30 60 90

sCDFT/0

CDFT/50

sFODFT

sFODFTB

MRCI+Q

56.5 (78.4) 56.4 (78.2) 56.1 (77.8) 55.9 (77.6)

68.0 67.7 67.5 68.2

43.9 (32.6) 43.8 (32.5) 43.5 (32.3) 43.4 (32.2)

54.9 (35.6) 54.6 (35.9) 54.1 (35.1) 53.8 (34.9)

68.5 68.9 68.8 69.8

MRUE changes to 36.7% for FODFT and to 21.3% in case of FODFTB. Scaling of CDFT/0 values leads to a change from overestimation (MRSE = +13.1% and MRUE = 13.1%) to underestimation (MRSE = −18.5% and MRUE = 18.5%) of electronic couplings. To understand the weak dependence on the rotation angle α visible at any level of theory, we recall that the electronic coupling matrix elements are approximately proportional to the overlap of HOMOs of the ethylene subunits. Since these orbitals do not feature any nodal plane perpendicular to the bond axis, no significant change in Hab should be expected upon such rotation. A greater dependence on the orientation angle is expected in case of ET in negatively charged dimers (not considered here), where the HOMO is a π ∗ orbital, which indeed features a nodal plane perpendicular to the bond plane. Larger changes in electronic coupling matrix elements are expected in case of the thiophene dimer. We considered three types of rotation: (i) one thiophene molecule is rotated by an angle δ around the centre of mass (see Figure 4(b)) as in case of ethylene, (ii) the two fragments are simultaneously rotated by an angle of +γ and –γ , respectively (Figure 4(c)), and (iii) one thiophene molecule is randomly rotated around the centre of mass (Figure 4(d)). All methods employed show correctly a minimum coupling of 0.0 meV for the angle δ = 90◦ in the sandwich rotation (i) (Figure 5(a)). Once again, these findings can be rationalized in terms of the orbital overlap of the fragment’s HOMOs. Accordingly, one would expect the maximum of the electronic coupling matrix elements for δ = 0◦ and δ = 180◦ . The configuration at δ = 90◦ leads to an almost perfect cancellation of negative and positive phases of HOMO orbitals and, in consequence, to negligible overlap and Hab values. Although a rotation about 180◦ leads to a

FIG. 4. The examined rotations: (a) one ethylene molecule was rotated around the centre of mass by α = 30◦ , 60◦ , and 90◦ ; (b) one thiophene molecule was rotated around the centre of mass by an angle δ between 0◦ and 180◦ in steps of 20◦ ; (c) simultaneous disrotation of two thiophene molecules by angle γ between 0◦ and 90◦ in steps of 10◦ ; (d) one thiophene molecule was randomly rotated around the centre of mass (black dot). The distance between molecules in case (a) and (b) is 5.0 Å. It was set to 6.57 Å in (c) to keep a minimal distance of 2.0 Å between hydrogen atoms when γ = 90◦ . The distance of 4.0 Å between the centres of mass was chosen in case (d) to span a large range of couplings.


104105-15

Kubas et al.

J. Chem. Phys. 140, 104105 (2014)

FIG. 5. Change of electronic coupling matrix element |Hab | in thiophene dimers with variation of angles (a) δ (see Figure 4(b)) and (b) γ (Figure 4(c)). Panels on the right show the curves after application of the scaling factors to CDFT/0, FO-DFT, and FODFTB methods, summarised in Table VII.

structure which is not identical to the one at 0◦ , the coupling values are very similar. In case (ii) almost all methods show a maximum for γ = 60◦ (only at CDFT/0 the maximum is shifted by 10◦ , see Figure 5(b)). At this geometry the overlap between two 2p orbitals of carbon atoms is the largest what leads to largest coupling. In general, for both rotations considered the agreement is qualitatively correct for all methods. CDFT/50 is the only approach that predicts almost quantitatively correct couplings, while as indicated in Figure 5, the relative errors of CDFT/0, FODFT, and FODFTB can be significantly reduced through scaling. In case (iii) we considered 15 random orientations of one thiophene molecule. This is the only case in this work where we deal with dimers that possess C1 symmetry, i.e., Eq. (4) cannot be applied directly. Instead, we used the GMH approach to obtain reference electronic coupling matrix elements via Eq. (5). According to Table IX, the performance of FODFT and FODFTB approaches is comparable to the full HAB11 dataset. The CDFT/0 method yields MRUE and MRSE of 124.9% which is significantly larger than 38.7% obtained for the database. Inclusion of 50% of exact exchange in the functional allows to decrease the MRUE to 23.3% which is still somehow higher than 5.3% observed for the symmetric dimers. To explain such large differences in case of

CDFT/0 calculations we first inspected the generated dimers and noted that for some dimers the shortest interatomic distances were smaller than the sum of van der Waals radii. Especially in cases of short S. . . S contacts (e.g., 3.4 Å whereas the vdW radius of sulphur is 1.8 Å108 ) strong interaction between fragments as well as a shift of electron density away from sulphur can be expected. This is observed in the CASSCF wavefunction as a strong mixing of the configuration with an unpaired electron located mainly on the sulphur 3pz orbital with the typical configuration where an unpaired electron TABLE IX. Statistical evaluation of the electronic coupling matrix elements |Hab | calculated for 15 randomly rotated thiophene dimers (see Figure 4(d)). Reference data were obtained at the NEVPT2 level using the GMH approach. The performance of the scaled methods (“s” prefix) is also given. Method CDFT/50 CDFT/0 FODFT FODFTB sCDFT/0 sFODFT sFODFTB

MUE (meV)

MRSE (%)

MRUE (%)

MAX (meV)

26.8 160.3 48.2 64.9 72.1 25.6 26.8

17.1 124.9 − 24.8 − 38.0 62.2 1.4 − 4.5

23.3 124.9 26.6 38.0 62.3 17.7 17.1

55.8 256.0 94.4 110.3 158.6 58.7 69.2


104105-16

Kubas et al.

J. Chem. Phys. 140, 104105 (2014)

occupies a carbon-based orbital. This was not the issue for symmetric cases as the transition between the two orbitals was forbidden according to the Laporte rule. Such orbital near-degeneracy resulted to some degree of symmetrybreaking in CDFT calculations which deteriorated the results, in particular for CDFT/0 method. We note, however, that addition of 50% of HFX allows for proper treatment of this situation. For larger systems the problem can be circumvent with the FO-based methods where the orbital mixing is alleviated by construction. We should also note that the scaling factors derived in Sec. V C still uniformly improve the results. E. Acenes

We now proceed to the largest molecules in this study, a series of condensed aromatic rings (acenes) with systematically increased π -conjugation: benzene, naphthalene,

FIG. 6. Homodimers of acenes (n = 0–4) were chosen to study electron transfer between fragments of systematically longer conjugation length. The total charge of the dimers was always +1. 0: Benzene, 1: naphthalene, 2: anthracene, 3: tetracene, and 4: pentacene.

anthracene, tetracene, and pentacene (Figure 6). The last three species are interesting from a practical point of view as they show photoinduced charge-transfer excitation,109, 110 whereas crystalline pentacene is widely used as an organic semiconducting material.111 The reference SCS-CC2 calculations were performed in D2h point group for all dimers. The results are summarized in Table X and illustrated in

TABLE X. Electronic coupling matrix elements |Hab | (meV) and decay constants β (1/Å) calculated with CDFT/50 and scaled CDFT/0, FO-DFT, and FODFTB approaches and compared to reference SCS-CC2 values. Statistical evaluation of both quantities are given by means of mean un signed error [MUE = ( n |ycalc − yref |)/n], mean relative signed error [MRSE = ( n ((ycalc − yref )/yref ))/n], mean relative unsigned error [MRUE = ( n (|ycalc − yref |/yref ))/n], and maximum unsigned error (MAX = max |ycalc − yref |). MUE, MRSE, MRUE, and MAX for unscaled methods are given in brackets. Dimer

sCDFT/0

Benzene

|Hab |

3.5 Å 4.0 Å 4.5 Å 5.0 Å

Naphthalene

|Hab |

3.5 Å 4.0 Å 4.5 Å 5.0 Å

Anthracene

|Hab |

3.5 Å 4.0 Å 4.5 Å 5.0 Å

β

β

β Tetracene

|Hab |

3.5 Å 4.0 Å 4.5 Å 5.0 Å β

Pentacene

|Hab |

3.5 Å 4.0 Å 4.5 Å 5.0 Å

|Hab |

MUE (meV) MRSE (%) MRUE (%) MAX (meV)

β

MUE (1/Å) MRSE (%) MRUE (%) MAX (1/Å)

β

467.0 232.0 112.2 53.2 2.90 452.8 223.7 109.3 51.4 2.90 436.0 212.8 104.3 49.0 2.91 423.6 204.7 100.3 47.5 2.91 294.8 198.0 96.9 46.0 2.52 20.2 (86.7) 11.9 (55.2) 15.0 (55.2) 87.5 (209.9) 0.32 − 7.5 10.1 0.77

CDFT/50 475.5 223.5 105.6 49.8 3.01 469.0 216.6 101.4 47.6 3.05 447.3 202.8 93.8 43.7 3.10 428.8 191.2 87.9 40.8 3.13 413.8 181.7 82.7 37.8 3.19 13.9 5.9 7.9 51.0 0.09 1.0 3.2 0.29

sFO-DFT 461.6 208.7 93.2 41.2 3.22 451.4 205.0 93.0 42.1 3.16 438.6 198.6 89.6 40.4 3.18 430.0 193.3 86.8 38.6 3.21 423.6 189.4 84.6 37.7 3.23 13.7 (42.1) 1.9 (−24.4) 7.7 (24.4) 41.3 (101.2) 0.15 4.6 5.3 0.50

sFODFTB 461.9 219.0 99.3 42.3 3.19 428.5 194.5 84.7 34.8 3.34 403.8 178.5 75.8 30.5 3.44 389.6 169.3 70.7 28.0 3.51 380.3 163.3 67.4 26.5 3.55

SCS-CC2 443.6 220.2 110.3 57.9 2.72 418.0 220.3 93.7 42.7 3.08 403.3 191.0 88.1 39.5 3.10 394.1 184.3 83.4 36.1 3.19 382.3 176.8 78.4 32.5 3.28

10.3 (69.7) − 10.1 (−41.7) 10.7 (41.7) 25.8 (143.7) 0.33 11.0 11.0 0.47


104105-17

Kubas et al.

J. Chem. Phys. 140, 104105 (2014)

FIG. 7. Correlation between calculated unscaled (left) and scaled (right) CDFT/0, CDFT/50, FO-DFT, FODFTB electronic couplings and reference SCS-CC2 data for acene homo-dimer cations. The black dashed line represents perfect correlation.

Figure 7. Examining the reference calculations first, we find that the couplings decrease with the number of rings in the system (higher n). This observation can be explained in terms of a simple particle in a box model, where the density of electronic levels increases with increasing energy levels. In consequence, the excitation energy between HOMO−1 and HOMO (CT excitation in our case), and therefore electronic coupling (see Eq. (4)), becomes smaller. Moreover, we find that the decay constant monotonically increases with n. A possible explanation is that with increasing charge delocalization over the aromatic rings, the excess charge distribution perpendicular to the molecular plane (in direction of the other monomer) decays faster. Hence, the decay constants increase with increasing n. Considering the DFT approaches, the unscaled CDFT/50 provides excellent agreement for both Hab (MRUE = 7.9%) and β values (MRUE = 3.2%). The only large deviation is the β value in case of benzene (about 0.3 Å−1 ). Overall, the performance of CDFT/50 is similarly good as for the database of small molecules, indicating that the results are likely to be transferrable to larger π -conjugated systems. The scaled variants of CDFT/0, FODFT, and FODFTB give MRUEs of 15.0%, 7.7%, and 10.7% for Hab and 10.1%, 5.3%, and 11.0% for β values, respectively. Similar errors were obtained for the database, showing again that the results can be likely generalized to larger systems. Again, the improvement through scaling is significant (see Figure 7), reducing the error in Hab by a factor of 4 and 2 in case of CDFT/0 and FODFTB while an order of magnitude for FODFT. We noticed that the sensitivity of β values with respect to Hartree-Fock exchange was larger for the acenes than for the molecules in the database. The average change for the 11 molecules in the database was 0.08 Å−1 when going from CDFT/0 to CDFT/50 while in case of acenes it is on average 0.19 Å−1 . We decided to investigate this phenomenon in more detail and plotted the density contour plot of tetracene dimer for both methods (Figure 8). In accord with expectations the spin density calculated with hybrid functionals is more compact, and hence, the decay constant is higher.

VI. DISCUSSION A. CDFT and fraction of HFX

Overall, the dependence of CDFT on the fraction of HFX is rather weak. We found that CDFT in combination with a functional including 50% HFX reproduces the HAB11 reference calculations very accurately (MRSE = –0.8%) but errors for pure GGA (+38.7%) or 100% of HFX (–18.2%) are reasonably small. An analysis of the spin density distribution shows that the excess charge on the donor or acceptor becomes somewhat more compact and less diffuse with increasing the fraction of HFX (this is especially visible for larger systems, see, for example, Figure 8 in the section discussing couplings in acenes). Since the electronic coupling is approximately proportional to the overlap between the orbitals that carry the excess charge in the two diabatic states, the coupling will be smaller for more compact spin densities, i.e., for increasing fractions of HFX. We also note that exact exchange might be important for systems with near-degenerated electronic states like in case of rotated thiophene dimers (see Sec. V D).

FIG. 8. Contour plot of the spin densities for the π -stacked tetracene dimer at an intermolecular separation of 5 Å. The contours were plotted for a plane in perpendicular orientation to the molecular plan, cutting through C-atoms 1-4 and 1 -4 . The isovalue was set to 10% of the maximal value for a given method (CDFT/0 or CDFT/50). Solid and dashed lines are used to distinguish between the two charge transfer states.


104105-18

Kubas et al.

Interestingly, for the mixed-valence tetrathiafulvalenediquinone anion29, 112 and electron hole transfer between F 0 and F + centres (oxygen vacancies) in crystalline MgO54, 55 CDFT couplings decreased much stronger with increasing fraction of HFX than for the dimers in the present study. In the two systems the medium between donor and acceptor is comprised of a “bridge” of either C atoms or Mg and O ions, respectively. This indicates that results for charge transfer through bridged donor-acceptor systems are potentially more sensitive to the fraction of HFX than for tunnelling through vacuum. Increasing the fraction of HFX is likely to open up the energy gap between donor/acceptor and bridge states (similar to the band gap problem in bulk water and semiconductors), leading to a much reduced mixing of donor/acceptor states with bridge states and a faster decay of initial and final state wavefunction, and in consequence to significantly smaller coupling values. The optimum fraction for electron tunnelling between oxygen vacancies in MgO was found to be 25%, which is smaller than the optimised value of 35% HFX used by Renz et al.113, 114 to describe intramolecular charge transfer excitations in π -conjugated organic molecules, and 50% of HFX (BHandLYP functional) was suggested in some recent studies.115, 116 Thus, the optimal amount of HFX for such bridged complexes is system-dependent and careful benchmark of the fraction of HFX is necessary, whereas for tunnelling through vacuum the amount of HFX seems to be less important.

B. FODFT and effect of electron number

In Sec. III C we showed via Eqs. (19)–(22) that the FODFT method is an approximation to the CDFT approach. In fact in the limit of non-interacting fragments both methods become equivalent. Our benchmark calculations on the HAB11 database clearly show that in practical calculations FODFT gives couplings that are a factor of about 1.5 smaller than the reference values what is most probably a collective effect of missing interactions between two fragments, like polarization or electronic stabilization. In addition, there are several different implementations of FODFT methods that differ in, e.g., the number of electrons used to construct the Hamiltonian, the type of basis set and the orthogonalisation procedure used for orbitals (see Sec. III C). We have summarized the results obtained with different FODFT schemes for the HAB11 database at vander Waals separation distance (3.5 Å) in Table XI and would like to discuss some important observations. The column “FODFT(2N−1)” refers to the standard FODFT method described in Sec. III C, using the correct number of electrons for construction of the Hamiltonian, a plane wave basis and Löwdin orthogonalisation of all 2N orbitals. Because there is no direct interaction between the monomers, this method yields generally smaller couplings than CDFT for the same density functional (CDFT/0). Column “FODFT(2N)” refers to the same calculation except that the Hamiltonian is constructed from 2N instead of 2N − 1 electrons. The consequence of the extra electron in the Hamiltonian is that the couplings increase compared to FODFT(2N−1), to val-

J. Chem. Phys. 140, 104105 (2014) TABLE XI. Hab (meV) calculated with various FODFT approaches for HAB11 dimers at intermolecular separation of 3.5 Å.

Ethylene Acetylene Cyclopropene Cyclobutadiene Cyclopentadiene Furane Pyrrole Thiophene Imidazole Benzene Phenol

FODFT (2N − 1)

FODFT (2N)

ADF FODFT(2N)

FODFTB (2N)

Ref.

367.7 316.9 418.8 323.3 343.3 315.6 328.7 341.2 310.7 342.4 190.5

389.2 345.8 443.7 346.9 360.6 334.0 347.8 357.8 328.9 353.5 211.3

388.4 345.3 439.4 345.6 358.7 333.7 347.7 356.1 328.2 354.1 279.5

343.7 212.0 367.4 261.6 283.2 280.3 286.2 264.8 277.5 299.9 231.4

519.2 460.7 536.6 462.7 465.8 440.3 456.3 449.0 411.6 435.2 375.0

ues that are slightly closer to the reference values. The same trend can be observed for small diatomic systems (see Tables III and IV). Interestingly, the ADF FODFT(2N) calculations give electronic couplings that are very close to the CPMD FODFT(2N) values; the difference is – in most cases – not more than about 1%. The only exceptions are FODFT(2N/2N−1) calculations on phenol dimer, where couplings are decreased by the same mechanism as explained for CDFT/0 in the supplementary material.56 This shows, however, that the basis set (Slater-orbitals in ADF and plane waves in CPMD) and the slightly different orthogonalization procedure have very little effect. We should also note that for both FODFT(2N) methods a slightly smaller scaling factor of 1.313 could be applied to systematically improve computed absolute values of electronic coupling matrix elements. In Table XII we summarize the beta-values obtained for some of the FODFT methods on HAB11. As one can see, the FODFT(2N) method gives beta-values that are closer to the reference values than FODFT(2N−1). The same trend is observed for small diatomics (Tables III and IV). We consider this as a fortuitous coincidence. The interaction term with the additional electron in the matrix element leads to an increase in electronic coupling and in particular to a smaller exponential distance decay in the FODFT(2N) method. Hence,

TABLE XII. Decay constant β (1/Å) calculated with various FODFT approaches for HAB11 dimers.

Ethylene Acetylene Cyclopropene Cyclobutadiene Cyclopentadiene Furane Pyrrole Thiophene Imidazole Benzene Phenol

FODFT(2N − 1)

FODFT(2N)

FODFTB(2N)

Ref.

3.23 3.42 3.44 3.15 3.15 3.25 3.18 3.22 3.25 3.22 2.91

2.75 2.86 3.10 2.75 2.83 2.92 2.83 2.92 2.91 2.96 2.74

3.00 2.96 3.23 3.37 3.27 3.30 3.29 2.82 3.28 3.16 3.28

2.70 2.80 3.06 2.68 2.89 3.01 2.89 2.82 2.82 2.85 2.95


104105-19

Kubas et al.

J. Chem. Phys. 140, 104105 (2014)

the presence of an additional electron partly replaces the effects of electronic polarization between the fragment orbitals, which is missing in all FODFT methods. The result is that FODFT(2N), which formally has one electron too much in the dimer Hamiltonian, gives better agreement with reference calculations than the FODFT(2N−1) method with the correct number of electrons. C. Model-Hamiltonian for nonadiabatic dynamics and Landauer theory vs diabatic states

In systems with more than two electron accepting/donating groups such as, e.g., DNA,49, 117–119 the molecular Hamiltonian can be used to determine the electron dynamics of the system.38, 39 The latter can give rise to rich quantum mechanical phenomena in particular when coupled to the nuclear dynamics as in non-adiabatic molecular dynamics (MD) approaches.38, 39 In principle, such a procedure can be understood in the framework of time-dependent DFT as formulated by Runge and Gross.120 Here, time-dependent equations for the electronic degrees of freedom can be derived as an exact reformulation of time-dependent quantum mechanics, leading to the time-dependent Kohn-Sham equations:121 ˙ i = Hˆ KS i . i¯

(30)

In short, this allows the usage of a simple DFT Hamiltonian for the propagation of the Kohn-Sham states. For large systems, the dynamics of the system cannot be treated with DFT, therefore often classical MD trajectories are used. These, however, are computed for the charge neutral system, and the question arises how to construct the KS-Hamiltonian used in Eq. (30). In the DFTB propagation scheme of the electron/hole wavefunction developed by some of us38, 122 we start from the DFT total energy expressed in terms of a FO Hamiltonian for the neutral system39, 120 (see for instance our recent derivations Refs. 38 and 118). We would like to point out that the off-diagonal elements of this FO-Hamiltonian are the same integrals that are used in this work (in the FODFTB(2N) method) to approximate the electronic coupling matrix element Hab between the charged diabatic states for electron transfer. However, since the FO Hamiltonian is constructed for the neutral system, the FO states in the DFTB charge propagation scheme are formally not the same as the diabatic states in Marcus ET theory (the latter carrying the total net charge of the system). Thus, for construction of the FO Hamiltonian in our DFTB charge propagation scheme (and similarly for calculations of Landauer currents), the off-diagonal elements obtained according to FODFTB(2N) should not be scaled, whereas for calculation of Hab for Marcus-type charge transfer reactions, these couplings should be scaled as recommended in Sec. V C. VII. CONCLUSIONS

In this work the HAB11 database of 11 π -stacked organic dimers was proposed and high-level ab initio reference values for electronic coupling matrix elements were reported to aid the assessment of DFT-based approaches. We found that the

NEVPT2 couplings are in very good agreement with the most accurate MRCI+Q data reported in this work. The approximate coupled cluster model (CC2) gives good couplings only after spin-component scaling (SCS-CC2). Thus, a hierarchy of wavefunction-based methods of decreasing computational complexity, yet of still sufficient accuracy is proposed in this work for calculations of electron transfer coupling matrix elements. Interestingly, similar high-level calculations for electron transfer (i.e., for the dimer anions) were not possible, as none of the dimers that could be treated at the MRCI+Q or NEVPT2 level of theory (up to benzene dimer) binds an extra electron. Among the DFT-based approaches, CDFT/50 was found to show best performance overall (–0.8%), but at short distances (3.5 Å) it was outperformed by CDFT/60. This indicates that range-separated functionals123 (which we have not investigated in this work) could further improve the overall quality of the results. However, the dependence on HFX is not very pronounced for the systems studied here. As a rough guide, we recommend a value between 25% and 50%. To get accurate results, in particular for bridged systems, careful benchmarking of this parameter is necessary. We have shown that economic approaches such as FODFT and FODFTB already provide chemically accurate estimates for electronic couplings in the HAB11 database. By this we mean that the typical error of these methods, an underestimation of Hab by a factor of about 1.5 translates to an underestimation of the non-adiabatic ET rate (proportional to |Hab |2 ) by a factor of about 2. This is well within the order of magnitude error to which chemical accuracy refers to. This should be compared to the typical error one makes for the calculations of the other quantities that enter the rate expression1, 124 such as reorganization free energy λ125–127 and free energy difference G.46, 128 For instance, an error in λ of 0.2 eV, which is easily made in particular in the condensed phase, translates into an error in the rate of one order of magnitude (at zero driving force). This is not to say that the only difficulty is to compute free energies. Indeed, there are situations where the coupling calculation is still very challenging, e.g., when there is an intervening medium (other than vacuum) and/or when donor-acceptor distance is so large that the actual coupling is smaller than the numerical accuracy of the calculations. We have also shown that the performance of FODFT and FODFTB can be further improved by application of a uniform scaling factor. In this way relative errors are reduced by at least a factor of 2. The scaled methods as well as our best unscaled CDFT/50 method were subsequently applied to some rotated isomers of selected molecules from the database as well as to a series of acenes. We noted significant improvement when applying the scaling factor for systems outside the database. These test cases suggest that our findings can be likely generalized to arbitrary orientations and larger π conjugated systems relevant to organic semiconducting materials and DNA. We conclude that when high-level ab initio calculations are impractical due to system size, CDFT with 50% of exact exchange should be a good choice to describe charge transfer in π -conjugated organic systems. For large molecules, or to


104105-20

Kubas et al.

obtain fast estimates, scaled versions of FODFT or FODFTB are recommended.

ACKNOWLEDGMENTS

A.K. was supported by EPSRC Grant No. EP/J015571/1. J.B. thanks the Royal Society for a University Research Fellowship. H.O. acknowledges the Alexander v. Humboldt Foundation and the state of Bavaria (grant “solar technologies go hybrid”) for support. This work was carried out on HECToR (Edinburgh), to which access was granted via the Materials Chemistry Consortium (EPSRC Grant Nos. EP/F067496 and EP/L000202). 1 R.

A. Marcus, J. Chem. Phys. 24, 966 (1956). Warshel, J. Phys. Chem. 86, 2218 (1982). 3 M. D. Newton, Chem. Rev. 91, 767 (1991). 4 E. Rosta and A. Warshel, J. Chem. Theory Comput. 8, 3574 (2012). 5 Q. Wu, B. Kaduk, and T. Van Voorhis, J. Chem. Phys. 130, 034109 (2009). 6 G. Voth, Acc. Chem. Res. 39, 143 (2006). 7 A. J. Cohen, P. Mori-Sanchez, and W. Yang, Science 321, 792 (2008). 8 T. Pacher and L. S. Cederbaum, J. Chem. Phys. 89, 7367 (1988). 9 W. Domcke and C. Woywood, Chem. Phys. Lett. 216, 362 (1993). 10 R. J. Cave and M. D. Newton, J. Chem. Phys. 106, 9213 (1997). 11 R. J. Cave and M. D. Newton, Chem. Phys. Lett. 249, 15 (1996). 12 R. S. Mulliken, J. Am. Chem. Soc. 74, 811 (1952). 13 R. S. Mulliken and W. B. Person, Molecular Complexes (Wiley, New York, 1969). 14 N. S. Hush, Prog. Inorg. Chem. 8, 391 (1967). 15 N. S. Hush, Electrochim. Acta 13, 1005 (1968). 16 J. R. Reimers and N. S. Hush, J. Phys. Chem. 95, 9773 (1991). 17 C. Creutz, M. D. Newton, and N. Sutin, J. Photochem. Photobiol. A 82, 47 (1994). 18 A. A. Voityuk and N. Rösch, J. Chem. Phys. 117, 5607 (2002). 19 C.-P. Hsu, Z.-Q. You, and H.-C. Chen, J. Phys. Chem. C 112, 1204 (2008). 20 A. Farazdel, M. Dupuis, E. Clementi, and A. Aviram, J. Am. Chem. Soc. 112, 4206 (1990). 21 K. Senthilkumar, F. C. Grozema, F. M. Bickelhaupt, and L. D. A. Siebbeles, J. Chem. Phys. 119, 9809 (2003). 22 X. Zeng, X. Hu, and W. Yang, J. Chem. Theory Comput. 8, 4960 (2012). 23 D. M. A. Smith, K. M. Rosso, M. Dupuis, M. Valiev, and T. P. Straatsma, J. Phys. Chem. B 110, 15582 (2006). 24 A. Migliore, S. Corni, R. Di Felice, and E. Molinari, J. Chem. Phys. 124, 064501 (2006). 25 A. Migliore, P. H.-L. Sit, and M. L. Klein, J. Chem. Theory Comput. 5, 307 (2009). 26 Q. Wu and T. Van Voorhis, J. Chem. Phys. 125, 164105 (2006). 27 Q. Wu and T. Van Voorhis, J. Chem. Theory Comput. 2, 765 (2006). 28 A. de la Lande and D. R. Salahub, J. Mol. Struct.: THEOCHEM 943, 115 (2009). 29 H. Oberhofer and J. Blumberger, J. Chem. Phys. 131, 064101 (2009). 30 H. Oberhofer and J. Blumberger, J. Chem. Phys. 133, 244105 (2010). 31 M. Pavanello and J. Neugebauer, J. Chem. Phys. 135, 234103 (2011). 32 M. Pavanello, T. Van Voorhis, L. Visscher, and J. Neugebauer, J. Chem. Phys. 138, 054101 (2013). 33 J. E. Subotnik, S. Yeganeh, R. J. Cave, and M. A. Ratner, J. Chem. Phys. 129, 244101 (2008). 34 J. E. Subotnik, R. J. Cave, R. P. Steele, and N. Shenvi, J. Chem. Phys. 130, 234102 (2009). 35 P. A. Pieniazek, S. A. Arnstein, S. E. Bradforth, A. I. Krylov, and C. D. Sherrill, J. Chem. Phys. 127, 164110 (2007). 36 L. Blancafort and A. A. Voityuk, J. Phys. Chem. A 110, 6426 (2006). 37 L. Blancafort and A. A. Voityuk, J. Phys. Chem. A 111, 4714 (2007). 38 T. Kubaˇr and M. Elstner, J. Phys. Chem. B 114, 11221 (2010). 39 T. Kubaˇr and M. Elstner, J. R. Soc. Interface 10, 20130415 (2013). 40 V. Coropceanu, J. Cornil, D. A. da Silva, Y. Olivier, R. Silbey, and J.-L. Bredas, Chem. Rev. 107, 926 (2007). 41 J. J. Kwiatkowski, J. M. Frost, and J. Nelson, Nano Lett. 9, 1085 (2009). 42 A. Troisi, D. L. Cheung, and D. Andrienko, Phys. Rev. Lett. 102, 116602 (2009). 2 A.

J. Chem. Phys. 140, 104105 (2014) 43 H.

Oberhofer and J. Blumberger, Phys. Chem. Chem. Phys. 14, 13846 (2012). 44 H. Oberhofer and J. Blumberger, Angew. Chem., Int. Ed. 49, 3631 (2010). 45 F. Gajdos, H. Oberhofer, M. Dupuis, and J. Blumberger, J. Phys. Chem. Lett. 4, 1012 (2013). 46 M. Breuer, K. M. Rosso, and J. Blumberger, Proc. Natl. Acad. Sci. U.S.A. 111, 611 (2014). 47 T. Kubaˇr, P. B. Woiczikowski, G. Cuniberti, and M. Elstner, J. Phys. Chem. B 112, 7937 (2008). 48 P. B. Woiczikowski, T. Kubaˇr, R. Gutírrez, G. Cuniberti, and M. Elstner, J. Chem. Phys. 133, 035103 (2010). 49 P. B. Woiczikowski, T. B. Steinbrecher, T. Kubaˇr, and M. Elstner, J. Phys. Chem. B 115, 9846 (2011). 50 A. Heck, P. B. Woiczikowski, T. Kubaˇr, B. Giese, M. Elstner, and T. B. Steinbrecher, J. Phys. Chem. B 116, 2284 (2012). 51 M. Wolter, P. B. Woiczikowski, M. Elstner, and T. Kubaˇr, Phys. Rev. B 85, 075101 (2012). 52 G. Lüdemann, P. B. Woiczikowski, T. Kubaˇr, M. Elstner, and T. B. Steinbrecher, J. Phys. Chem. B 117, 10769 (2013). 53 M. Wolter, M. Elstner, and T. Kubaˇr, J. Chem. Phys. 139, 125102 (2013). 54 K. P. McKenna and J. Blumberger, Phys. Rev. B 86, 245110 (2012). 55 J. Blumberger and K. P. McKenna, Phys. Chem. Chem. Phys. 15, 2184 (2013). 56 See supplementary material at http://dx.doi.org/10.1063/1.4867077 for Cartesian coordinates and isosurfaces comparison and discussion of the value of integrated spin density as a diagnostic tool in CDFT calculations. 57 P. A. M. Dirac, Proc. R. Soc. London, Ser. A 123, 714 (1929). 58 J. C. Slater, Phys. Rev. 81, 385 (1951). 59 S. H. Vosko, L. Wilk, and M. Nusair, Can. J. Phys. 58, 1200 (1980). 60 A. D. Becke, Phys. Rev. A 38, 3098 (1988). 61 J. P. Perdew, Phys. Rev. B 33, 8822 (1986). 62 K. Eichkorn, F. Weigend, O. Treutler, and R. Ahlrichs, Theor. Chem. Acc. 97, 119 (1997). 63 F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys. 7, 3297 (2005). 64 F. Weigend, M. Häser, H. Patzelt, and R. Ahlrichs, Chem. Phys. Lett. 294, 143 (1998). 65 K. Eichkorn, O. Treutler, H. Ohm, M. Häser, and R. Ahlrichs, Chem. Phys. Lett. 240, 283 (1995). 66 M. Sierka, A. Hogekamp, and R. Ahlrichs, J. Chem. Phys. 118, 9136 (2003). 67 TURBOMOLE version 6.5, a development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989–2007, TURBOMOLE GmbH, 2013, see http://www.turbomole.com. 68 P. G. Szalay, T. Müller, G. Gidofalvi, H. Lischka, and R. Shepard, Chem. Rev. 112, 108 (2012) and references therein. 69 S. R. Langhoff and E. R. Davidson, Int. J. Quantum Chem. 8, 61 (1974). 70 W. Butscher, S. Shih, R. J. Buenker, and S. D. Peyerimhoff, Chem. Phys. Lett. 52, 457 (1977). 71 C. Angeli, R. Cimiraglia, S. Evangelisti, T. Leininger, and J.-P. Malrieu, J. Chem. Phys. 114, 10252 (2001). 72 C. Angeli, R. Cimiraglia, and J.-P. Malrieu, Chem. Phys. Lett. 350, 297 (2001). 73 C. Angeli, R. Cimiraglia, and J.-P. Malrieu, J. Chem. Phys. 117, 9138 (2002). 74 O. Christiansen, H. Koch, and P. Jørgensen, Chem. Phys. Lett. 243, 409 (1995). 75 A. Hellweg, S. Grün, and C. Hättig, Phys. Chem. Chem. Phys. 10, 4119 (2008). 76 I. Schapiro, K. Sivalingam, and F. Neese, J. Chem. Theory Comput. 9, 3567 (2013). 77 P. H. Dederichs, S. Blugel, R. Zeller, and H. Akai, Phys. Rev. Lett. 53, 2512 (1984). 78 Q. Wu and T. Van Voorhis, Phys. Rev. A 72, 024502 (2005). 79 CPMD version 3.15.3, the CPMD consortium, MPI für Festkörperforschung, and the IBM Zurich Research Laboratory, 2013, see http://www.cpmd.org. 80 F. L. Hirshfeld, Theor. Chem. Acc. 44, 129 (1977). 81 P. O. Löwdin, J. Chem. Phys. 18, 365 (1950). 82 G. te Velde, F. M. Bickelhaupt, E. J. Baerends, C. Fonseca Guerra, S. J. A. van Gisbergen, J. G. Snijders, and T. Ziegler, J. Comput. Chem. 22, 931 (2001). 83 D. Porezag, T. Frauenheim, T. Köhler, G. Seifert, and R. Kaschner, Phys. Rev. B 51, 12947 (1995).


104105-21 84 G.

Kubas et al.

Seifert, D. Porezag, and T. Frauenheim, Int. J. Quantum Chem. 58, 185 (1996). 85 M. Elstner, D. Porezag, G. Jungnickel, J. Elsner, M. Haugk, T. Frauenheim, S. Suhai, and G. Seifert, Phys. Rev. B 58, 7260 (1998). 86 M. Gaus, Q. Cui, and M. Elstner, J. Chem. Theory Comput. 7, 931 (2011). 87 G. Seifert and J.-O. Joswig, WIREs Comput. Mol. Sci. 2, 456 (2012). 88 M. Gaus, Q. Cui, and M. Elstner, WIREs Comput. Mol. Sci. 4, 49 (2013). 89 M. Rapacioli and F. Spiegelman, Eur. Phys. J. D 52, 55 (2009). 90 M. Rapacioli, F. Spiegelman, A. Scemama, and A. Mirtschink, J. Chem. Theory Comput. 7, 44 (2011). 91 M. Rapacioli, A. Simon, L. Dontot, and F. Spiegelman, Phys. Status Solidi B 249, 245 (2012). 92 T. H. Dunning, Jr., J. Chem. Phys. 90, 1007 (1989). 93 D. E. Woon and T. H. Dunning, Jr., J. Chem. Phys. 98, 1358 (1993). 94 F. Neese, J. Chem. Phys. 119, 9428 (2003). 95 S. Grimme, L. Goerigk, and R. F. Fink, WIREs Comput. Mol. Sci. 2, 886 (2012). 96 F. Weigend, A. Kohn, and C. Hattig, J. Chem. Phys 116, 3175 (2002). 97 F. Neese, Comput. Mol. Sci. 2, 73 (2012). 98 C. Hättig and F. Weigend, J. Chem. Phys. 113, 5154 (2000). 99 E. Papajak, H. R. Leverentz, J. Zheng, and D. G. Truhlar, J. Chem. Theory Comput. 5, 1197 (2009). 100 D. Feller, J. Comput. Chem. 17, 1571 (1996). 101 K. L. Schuchardt, B. T. Didier, T. Elsethagen, L. Sun, V. Gurumoorthi, J. Chase, J. Li, and T. L. Windus, J. Chem. Inf. Model. 47, 1045 (2007). 102 A. J. H. Wachters, J. Chem. Phys. 52, 1033 (1970). 103 N. Troullier and J. Martins, Phys. Rev. B 43, 1993 (1991). 104 T. Niehaus, J. Mol. Struct.: THEOCHEM 541, 185 (2001). 105 C. L. Janssen and I. M. B. Nielsen, Chem. Phys. Lett. 290, 423 (1998). 106 L. Goerigk and S. Grimme, J. Chem. Phys. 132, 184103 (2010). 107 A. Dreuw and M. Head-Gordon, J. Am. Chem. Soc. 126, 4007 (2004).

J. Chem. Phys. 140, 104105 (2014) 108 A.

Bondi, J. Phys. Chem. 68, 441 (1964). Sebastian, G. Weiser, and H. Bässler, Chem. Phys. 61, 125 (1981). 110 L. Sebastian, G. Weiser, G. Peter, and H. Bässler, Chem. Phys. 75, 103 (1983). 111 Y. Y. Yamashita, Sci. Technol. Adv. Mater. 10, 024313 (2009). 112 B. Kaduk, T. Kowalczyk, and T. Van Voorhis, Chem. Rev. 112, 321 (2012). 113 M. Renz, K. Theilacker, C. Lambert, and M. Kaupp, J. Am. Chem. Soc. 131, 16292 (2009). 114 M. Renz, M. Kess, M. Diedenhofen, A. Klamt, and M. Kaupp, J. Chem. Theory Comput. 8, 4189 (2012). 115 R. J. Magyar and S. Tretiak, J. Chem. Theory Comput. 3, 976 (2007). 116 J. J. Eriksen, S. P. A. Sauer, K. V. Mikkelsen, O. Christiansen, H. J. Aa. Jensen, and J. Kongsted, Mol. Phys. 111, 1235 (2013). 117 P. B. Woiczikowski, T. Kubaˇr, R. Gutíerrez, R. A. Caetano, G. Cuniberti, and M. Elstner, J. Chem. Phys. 130, 215104 (2009). 118 T. Kubaˇr, R. Gutiérrez, U. Kleinekathöfer, G. Cuniberti, and M. Elstner, Phys. Status Solidi B 250, 2277 (2013). 119 F. C. Grozema, S. Tonzani, Y. A. Berlin, G. C. Schatz, L. D. A. Siebbeles, and M. A. Ratner, J. Am. Chem. Soc. 130, 5157 (2008). 120 E. Runge and E. K. U. Gross, Phys. Rev. Lett. 52, 997 (1984). 121 M. Marques and E. K. U. Gross, Annu. Rev. Phys. Chem. 55, 427 (2004). 122 T. Kubaˇr and M. Elstner, Phys. Chem. Chem. Phys. 15, 5794 (2013). 123 A. Savin, Theoretical and Computational Chemistry: Recent Developments and Applications of Modern Density Functional Theory (Elsevier Science B.V., Amsterdam, 1996), Vol. 4, p. 327. 124 A. Troisi, A. Nitzan, and M. A. Ratner, J. Chem. Phys. 119, 5782 (2003). 125 J. Blumberger and G. Lamoureux, Mol. Phys. 106, 1597 (2008). 126 R. Seidel, M. Faubel, B. Winter, and J. Blumberger, J. Am. Chem. Soc. 131, 16127 (2009). 127 J. Moens, R. Seidel, P. Geerlings, M. Faubel, B. Winter, and J. Blumberger, J. Phys. Chem. B 114, 9173 (2010). 128 Y. Tateyama, J. Blumberger, T. Ohno, and M. Sprik, J. Chem. Phys. 126, 204506 (2007). 109 L.


Erratum: "Electronic couplings for molecular charge transfer: Benchmarking CDFT, FODFT, and FODFTB against high-level ab initio calculations" [J. Chem. Phys. 140, 104105 (2014)].

Electronic Structure of Silicon Nanowires Matrix from Ab Initio Calculations.

Kinetic formulation of the Kohn-Sham Equations for ab initio electronic structure calculations.

Excited electronic states of thiophene: high resolution photoabsorption Fourier transform spectroscopy and ab initio calculations.

Electronic states of tetrahydrofurfuryl alcohol (THFA) as studied by VUV spectroscopy and ab initio calculations.

Benchmarking ab initio binding energies of hydrogen-bonded molecular clusters based on FTIR spectroscopy.

An ab initio electronic transport database for inorganic materials.

(4)He Thermophysical Properties: New Ab Initio Calculations.

Hydration of porphyrin and Mg-porphyrin: ab initio quantum mechanical charge field molecular dynamics simulations.

Ab initio molecular orbital calculations on furanose sugars: a study with the 6-31G basis set.

Ab initio calculations of the electronic structure of the low-lying states for the ultracold LiYb molecule.

Structural and electronic properties of organo-halide hybrid perovskites from ab initio molecular dynamics.

Hydration of the cyanide ion: an ab initio quantum mechanical charge field molecular dynamics study.

Determination of the electronic energy levels of colloidal nanocrystals using field-effect transistors and Ab-initio calculations.

The torsional potential for phosphate diesters. The effect of geometry optimization in CNDO and ab initio molecular orbital calculations.

Ab initio electronic properties of dual phosphorus monolayers in silicon.

Electronic State Spectroscopy of Halothane As Studied by ab Initio Calculations, Vacuum Ultraviolet Synchrotron Radiation, and Electron Scattering Methods.

Ab initio molecular dynamics with noisy forces: validating the quantum Monte Carlo approach with benchmark calculations of molecular vibrational properties.

Ab initio molecular simulations for proposing potent inhibitors to butyrylcholinesterases.

Ab Initio electronic properties of monolayer phosphorus nanowires in silicon.

Hydrated Electron Transfer to Nucleobases in Aqueous Solutions Revealed by Ab Initio Molecular Dynamics Simulations.

Molecular excited states: accurate calculation of relative energies and electronic coupling between charge transfer and non-charge transfer states.

Dissociation and recombination of D₂ on Cu(111): ab initio molecular dynamics calculations and improved analysis of desorption experiments.

Ab initio calculations for the far infrared collision induced absorption by N2 gas.