11.16 Electronic Couplings for Electron- and Energy Transfer

11.16.2 Diabatic-State-Based Methods

11.16.2.1 Electronic coupling in charge transfer

A charge transfer involves a change in the electron numbers in a pair of molecular fragments. As an example, we will use the following reaction when necessary, and a generalization to other cases is straightforward:

D-ADA- (11.98)

where an extra electron is localized to the donor (D) initially, and it becomes localized to the acceptor (A) in the final state. The two-state secular equation for the initial and final electronic states can be written as

𝐇-E𝐒=(Hii-SiiEHif-SifEHif-SifEHff-SffE)=0 (11.99)

This is very close to an eigenvalue problem except for the non-orthogonality between the initial and final states. A standard eigenvalue form for Eq. (11.99) can be obtained by using the Löwdin transformation:

𝐇eff=𝐒-1/2𝐇𝐒-1/2, (11.100)

where the off-diagonal element of the effective Hamiltonian matrix represents the electronic coupling for the reaction, and it is defined by

V=Hifeff=Hif-Sif(Hii+Hff)/21-Sif2 (11.101)

In a general case where the initial and final states are not normalized, the electronic coupling is written as

V=SiiSff×Hif-Sif(Hii/Sii+Hff/Sff)/2SiiSff-Sif2 (11.102)

Thus, in principle, V can be obtained when the matrix elements for the Hamiltonian H and the overlap matrix S are calculated.

The direct coupling (DC) scheme calculates the electronic coupling values via Eq. (11.102), and it is widely used to calculate charge transfer coupling.687, 121, 257, 1047 In the DC scheme, the coupling matrix element is calculated directly using charge-localized determinants (the “diabatic states” in electron transfer literature). In electron transfer systems, it has been shown that such charge-localized states can be approximated by symmetry-broken unrestricted Hartree-Fock (UHF) solutions.687, 121, 669 The adiabatic eigenstates are assumed to be the symmetric and anti-symmetric linear combinations of the two symmetry-broken UHF solutions in a DC calculation. Therefore, DC couplings can be viewed as a result of two-configuration solutions that may recover the non-dynamical correlation.

The core of the DC method is based on the corresponding orbital transformation469 and a calculation for Slater’s determinants in Hif and Sif.257, 1047 Unfortunately, the calculation of Hif is not available for DFT method because a functional of the two densities ρi and ρf is unknown and there are no existing approximate forms for Hif.1024 To calculate charge transfer coupling with DFT, we can use the CDFT-CI method (Section 5.13.3), the frontier molecular orbital (FMO) approach939, 1038 (Section 11.16.2.6) or a hybrid scheme – DC with CDFT wave functions (Section 11.16.2.4).

11.16.2.2 Corresponding orbital transformation

Let |Ψa and |Ψb be two single Slater-determinant wave functions for the initial and final states, and 𝐚 and 𝐛 be the spin-orbital sets, respectively:

𝐚 = (a1,a2,,aN) (11.103)
𝐛 = (b1,b2,,bN) (11.104)

Since the two sets of spin-orbitals are not orthogonal, the overlap matrix 𝐒 can be defined as:

𝐒=𝐛𝐚𝑑τ. (11.105)

We note that 𝐒 is not Hermitian in general since the molecular orbitals of the initial and final states are separately determined. To calculate the matrix elements Hab and Sab, two sets of new orthogonal spin-orbitals can be used by the corresponding orbital transformation.469 In this approach, each set of spin-orbitals 𝐚 and 𝐛 are linearly transformed,

𝐚^ = 𝐚𝐕 (11.106)
𝐛^ = 𝐛𝐔 (11.107)

where 𝐕 and 𝐔 are the left-singular and right-singular matrices, respectively, in the singular value decomposition (SVD) of 𝐒:

𝐒=𝐔𝐬^𝐕 (11.108)

The overlap matrix in the new basis is now diagonal

𝐛^𝐚^=𝐔(𝐛𝐚)𝐕=𝐬^ (11.109)

11.16.2.3 Generalized density matrix

The Hamiltonian for electrons in molecules are a sum of one-electron and two-electron operators. In the following, we derive the expressions for the one-electron operator Ω(1) and two-electron operator Ω(2),

Ω(1) = i=1Nω(i) (11.110)
Ω(2) = 12i,j=1Nω(i,j) (11.111)

where ω(i) and ω(i,j), for the molecular Hamiltonian, are

ω(i)=h(i)=-12i2+V(i) (11.112)

and

ω(i,j)=1rij (11.113)

The evaluation of matrix elements can now proceed:

Sab=Ψb|Ψa=det(𝐔)det(𝐕)i=1Ns^ii (11.114)
Ωab(1)=Ψb|Ω(1)|Ψa=det(𝐔)det(𝐕)i=1Nb^i|ω(1)|a^ijiNs^jj (11.115)
Ωab(2)=Ψb|Ω(2)|Ψa=12det(𝐔)det(𝐕)ijNb^ib^j|ω(1,2)(1-P12)|a^ia^jki,jNs^kk (11.116)
Hab=Ωab(1)+Ωab(2) (11.117)

In an atomic orbital basis set, {χ}, we can expand the molecular spin orbitals 𝐚 and 𝐛,

𝐚=χ𝐀,𝐚^=χ𝐀𝐕=χ𝐀^ (11.118)
𝐛=χ𝐁,𝐛^=χ𝐁𝐔=χ𝐁^ (11.119)

The one-electron terms, Eq. (11.114), can be expressed as

Ωab(1) = iNλσA^λiTiiB^iσχσ|ω(1)|χλ (11.120)
= λσGλσωσλ

where Tii=Sab/s^ii and define a generalized density matrix, 𝐆:

𝐆 = 𝐀^𝐓𝐁^ (11.121)

Similarly, the two-electron terms, Eq. (11.116), are

Ωab(2) = 12ijλσμνA^λiA^σj(1s^ii)TjjB^iμB^jνχμχν|ω(1,2)|χλχσ (11.122)
= λσμνGλμLGσνRμν||λσ

where 𝐆R and 𝐆L are generalized density matrices as defined in Eq. (11.121) except Tii in 𝐆L is replaced by 1/(2sii).

The α- and β-spin orbitals are treated explicitly. In terms of the spatial orbitals, the one- and two-electron contributions can be reduced to

Ωab(1) = λσGλσαωσλ+λσGλσβωσλ (11.123)
Ωab(2) = λσμνGλμLαGσνRα(μν|λσ-μν|σλ)+λσμνGλμLβGσνRαμν|λσ (11.124)
+λσμνGλμLαGσνRβμν|λσ+λσμνGλμLβGσνRβ(μν|λσ-μν|σλ)

The resulting one- and two-electron contributions, Eqs. (11.123) and (11.124) can be easily computed in terms of generalized density matrices using standard one- and two-electron integral routines in Q-Chem.

11.16.2.4 Direct coupling method for electronic coupling

It is important to obtain proper charge-localized initial and final states for the DC scheme, and this step determines the quality of the coupling values. Q-Chem provides three approaches to construct charge-localized states:

  • The “1+1” approach
    Since the system consists of donor and acceptor molecules or fragments, with a charge being localized either donor or acceptor, it is intuitive to combine wave functions of individual donor and acceptor fragments to form a charge-localized wave function. We call this approach “1+1” since the zeroth order wave functions are composed of the HF wave functions of the two fragments.

    For example, for the case shown in Example (11.98), we can use Q-Chem to calculate two HF wave functions: those of anionic donor and of neutral acceptor and they jointly form the initial state. For the final state, wave functions of neutral donor and anionic acceptor are used. Then the coupling value is calculated via Eq. (11.102).

    Example 11.38  To calculate the electron-transfer coupling for a pair of stacked-ethylene with “1+1” charge-localized states

    $molecule
      -1  2
    --
      -1  2,  0  1
       C    0.662489    0.000000    0.000000
       H    1.227637    0.917083    0.000000
       H    1.227637   -0.917083    0.000000
       C   -0.662489    0.000000    0.000000
       H   -1.227637   -0.917083    0.000000
       H   -1.227637    0.917083    0.000000
    --
       0  1, -1  2
       C    0.720595    0.000000    4.5
       H    1.288664    0.921368    4.5
       H    1.288664   -0.921368    4.5
       C   -0.720595    0.000000    4.5
       H   -1.288664   -0.921368    4.5
       H   -1.288664    0.921368    4.5
    $end
    
    $rem
       JOBTYPE           SP
       METHOD            HF
       BASIS             6-31G(d)
       SCF_PRINT_FRGM    FALSE
       SYM_IGNORE        TRUE
       SCF_GUESS         FRAGMO
       STS_DC            TRUE
    $end
    

    In the $molecule subsection, the first line is for the charge and multiplicity of the whole system. The following blocks are two inputs for the two molecular fragments (donor and acceptor). In each block the first line consists of the charge and spin multiplicity in the initial state of the corresponding fragment, a comma, then the charge and multiplicity in the final state. Next lines are nuclear species and their positions of the fragment. For example, in the above example, the first block indicates that the electron donor is a doublet ethylene anion initially, and it becomes a singlet neutral species in the final state. The second block is for another ethylene going from a singlet neutral molecule to a doublet anion.

    Note that the last three $rem variables in this example, SYM_IGNORE, SCF_GUESS and STS_DC must be set to be the values as in the example in order to perform DC calculation with “1+1” charge-localized states. An additional $rem variable, SCF_PRINT_FRGM is included. When it is TRUE a detailed output for the fragment HF self-consistent field calculation is given.

  • The “relaxed” approach
    In “1+1” approach, the intermolecular interaction is neglected in the initial and final states, and so the final electronic coupling can be underestimated. As a second approach, Q-Chem can use “1+1” wave function as an initial guess to look for the charge-localized wave function by further HF self-consistent field calculation. This approach would ‘relax’ the wave function constructed by “1+1” method and include the intermolecular interaction effects in the initial and final wave functions. However, this method may sometimes fail, leading to either convergence problems or a resulting HF wave function that cannot represent the desired charge-localized states. This is more likely to be a problem when calculations are performed with diffusive basis functions, or when the donor and acceptor molecules are very close to each other.

    To perform relaxed DC calculation, set STS_DC = RELAX.

  • A hybrid scheme – constrained DFT charge-localized states
    Constrained DFT (Section 5.13) can be used to obtain charge-localized states. It is recommended to set both charge and spin constraints in order to generate proper charge localization. To perform DC calculation with CDFT states, set SAVE_SUBSYSTEM = 10 and SAVE_SUBSYSTEM = 20 to save CDFT molecular orbitals in the first two jobs of a batch jobs, and then in the third job of the batch job, set SCF_GUESS = READ and STS_DC = TRUE to compute electronic coupling values.

    Example 11.39  To calculate the electron-transfer coupling for a pair of stacked-ethylene with CDFT charge-localized states

    $molecule
     -1  2
     C  0.662489  0.000000  0.000000
     H  1.227637  0.917083  0.000000
     H  1.227637 -0.917083  0.000000
     C -0.662489  0.000000  0.000000
     H -1.227637 -0.917083  0.000000
     H -1.227637  0.917083  0.000000
     C  0.720595  0.000000  4.500000
     H  1.288664  0.921368  4.500000
     H  1.288664 -0.921368  4.500000
     C -0.720595  0.000000  4.500000
     H -1.288664 -0.921368  4.500000
     H -1.288664  0.921368  4.500000
    $end
    
    $rem
       METHOD          =  wb97xd
       BASIS           =  cc-pvdz
       CDFT            =  true
       SYM_IGNORE      =  true
       SAVE_SUBSYSTEM  =  10
    $end
    
    $cdft
     1.0
     1.0  1 6
     1.0
     1.0  1 6  s
    $end
    
    @@@
    
    $molecule
     read
    $end
    
    $rem
       method          =  wb97xd
       basis           =  cc-pvdz
       cdft            =  true
       sym_ignore      =  true
       save_subsystem  =  20
    $end
    
    $cdft
     1.0
     1.0  7 12
     1.0
     1.0  7 12  s
    $end
    
    @@@
    
    $molecule
     read
    $end
    
    $rem
       method      =  hf
       basis       =  cc-pvdz
       sym_ignore  =  true
       scf_guess   =  read
       sts_dc      =  true
    $end
    

11.16.2.5 The frontier molecular orbital approach

The frontier molecular orbital (FMO) approach is often used with DFT to calculate ET coupling.939, 1038 FMO coupling value is essentially an off-diagonal Kohn–Sham matrix element with the overlap effect accounted

VFMO=fDA-S(fDD+fAA)/21-S2 (11.125)

where fDA=ϕFMOD|f^|ϕFMOA, with f^ being the Kohn–Sham operator of the donor-acceptor system. ϕFMOD(A) is the Kohn–Sham frontier molecular orbital for the donor (acceptor) fragment, which represents one-particle scheme of a charge transfer process.

11.16.2.6 The frontier molecular orbital approach

The frontier molecular orbital (FMO) approach is often used with DFT to calculate ET coupling.939, 1038 FMO coupling value is essentially an off-diagonal Kohn–Sham matrix element with the overlap effect accounted

VFMO=fDA-S(fDD+fAA)/21-S2 (11.126)

where fDA=ϕFMOD|f^|ϕFMOA, with f^ being the Kohn–Sham operator of the donor-acceptor system. ϕFMOD(A) is the Kohn–Sham frontier molecular orbital for the donor (acceptor) fragment, which represents one-particle scheme of a charge transfer process.

In this approach, computations are often performed separately in the two fragments, and the off-diagonal Kohn–Sham operator (and the overlap matrix) in the FMOs is subsequently calculated. To compute FMO couplings, Q-Chem a setup similar to the “1+1” approach

Example 11.40  To calculate the electron-transfer coupling for a pair of stacked-ethylene with the FMO approach

$molecule
 0  1
--
 0  1
 C  0.662489  0.000000  0.000000
 H  1.227637  0.917083  0.000000
 H  1.227637 -0.917083  0.000000
 C -0.662489  0.000000  0.000000
 H -1.227637 -0.917083  0.000000
 H -1.227637  0.917083  0.000000
--
 0  1
 C  0.720595  0.000000  4.5
 H  1.288664  0.921368  4.5
 H  1.288664 -0.921368  4.5
 C -0.720595  0.000000  4.5
 H -1.288664 -0.921368  4.5
 H -1.288664  0.921368  4.5
$end

$rem
   method              =  lrcwpbe
   omega               =  370
   basis               =  dz*
   scf_print_frgm      =  true
   sym_ignore          =  true
   scf_guess           =  fragmo
   sts_dc              =  fock
   sts_trans_donor     =  2-3    ! use HOMO, HOMO-1 and LUMO, LUMO+1, LUMO+2 of donor
   sts_trans_acceptor  =  1-2 ! use HOMO and LUMO, LUMO+1 of acceptor
$end

$rem_frgm
 print_orbitals = 5
$end

Note:  The FMOs are not always HOMO or LUMO of fragments. We can use STS_TRANS_DONOR (and STS_TRANS_ACCEPTOR) to select a range of occupied and virtual orbitals for FMO coupling calculations.