# 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:

 $\mathrm{D}^{-}\mathrm{A}\longrightarrow\mathrm{D}\mathrm{A}^{-}$ (11.97)

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

 ${\bf H}-E{\bf S}=\left(\begin{array}[]{cc}H_{\mathrm{ii}}-S_{\mathrm{ii}}E&H_{% \mathrm{if}}-S_{\mathrm{if}}E\\ H_{\mathrm{if}}-S_{\mathrm{if}}E&H_{\mathrm{ff}}-S_{\mathrm{ff}}E\end{array}% \right)=0$ (11.98)

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.98) can be obtained by using the Löwdin transformation:

 $\mathbf{H}_{\textrm{eff}}=\mathbf{S}^{-1/2}\mathbf{H}\mathbf{S}^{-1/2},$ (11.99)

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

 $V=H^{\textrm{eff}}_{\mathrm{if}}=\frac{H_{\mathrm{if}}-S_{\mathrm{if}}(H_{% \mathrm{ii}}+H_{\mathrm{ff}})/2}{1-S^{2}_{\mathrm{if}}}$ (11.100)

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

 $V=\sqrt{S_{\mathrm{ii}}S_{\mathrm{ff}}}\times\frac{H_{\mathrm{if}}-S_{\mathrm{% if}}(H_{\mathrm{ii}}/S_{\mathrm{ii}}+H_{\mathrm{ff}}/S_{\mathrm{ff}})/2}{S_{% \mathrm{ii}}S_{\mathrm{ff}}-S^{2}_{\mathrm{if}}}$ (11.101)

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.101), and it is widely used to calculate charge transfer coupling.681, 119, 256, 1036 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.681, 119, 663 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 transformation466 and a calculation for Slater’s determinants in $H_{\mathrm{if}}$ and $S_{\mathrm{if}}$.256, 1036 Unfortunately, the calculation of $H_{\mathrm{if}}$ is not available for DFT method because a functional of the two densities $\rho_{\mathrm{i}}$ and $\rho\mathrm{f}$ is unknown and there are no existing approximate forms for $H_{\mathrm{if}}$.1013 To calculate charge transfer coupling with DFT, we can use the CDFT-CI method (Section 5.13.3), the frontier molecular orbital (FMO) approach929, 1027 (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 $|\Psi_{a}\rangle$ and $|\Psi_{b}\rangle$ be two single Slater-determinant wave functions for the initial and final states, and $\mathbf{a}$ and $\mathbf{b}$ be the spin-orbital sets, respectively:

 $\displaystyle\mathbf{a}$ $\displaystyle=$ $\displaystyle\left(a_{1},a_{2},\cdots,a_{N}\right)$ (11.102) $\displaystyle\mathbf{b}$ $\displaystyle=$ $\displaystyle\left(b_{1},b_{2},\cdots,b_{N}\right)$ (11.103)

Since the two sets of spin-orbitals are not orthogonal, the overlap matrix $\mathbf{S}$ can be defined as:

 $\mathbf{S}=\int\mathbf{b}^{\dagger}\mathbf{a}\ d\tau.$ (11.104)

We note that $\mathbf{S}$ is not Hermitian in general since the molecular orbitals of the initial and final states are separately determined. To calculate the matrix elements $H_{ab}$ and $S_{ab}$, two sets of new orthogonal spin-orbitals can be used by the corresponding orbital transformation.466 In this approach, each set of spin-orbitals $\mathbf{a}$ and $\mathbf{b}$ are linearly transformed,

 $\displaystyle\hat{\mathbf{a}}$ $\displaystyle=$ $\displaystyle\mathbf{a}\mathbf{V}$ (11.105) $\displaystyle\hat{\mathbf{b}}$ $\displaystyle=$ $\displaystyle\mathbf{b}\mathbf{U}$ (11.106)

where $\mathbf{V}$ and $\mathbf{U}$ are the left-singular and right-singular matrices, respectively, in the singular value decomposition (SVD) of $\mathbf{S}$:

 $\mathbf{S}=\mathbf{U}\hat{\mathbf{s}}\mathbf{V}^{\dagger}$ (11.107)

The overlap matrix in the new basis is now diagonal

 $\int\hat{\mathbf{b}}^{\dagger}\hat{\mathbf{a}}=\mathbf{U}^{\dagger}\left(\int% \mathbf{b}^{\dagger}\mathbf{a}\right)\mathbf{V}=\hat{\mathbf{s}}$ (11.108)

## 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 $\Omega^{(1)}$ and two-electron operator $\Omega^{(2)}$,

 $\displaystyle\Omega^{(1)}$ $\displaystyle=$ $\displaystyle\sum_{i=1}^{N}\omega(i)$ (11.109) $\displaystyle\Omega^{(2)}$ $\displaystyle=$ $\displaystyle\frac{1}{2}\sum^{N}_{i,j=1}\omega(i,j)$ (11.110)

where $\omega(i)$ and $\omega(i,j)$, for the molecular Hamiltonian, are

 $\omega(i)=h(i)=-\frac{1}{2}\nabla_{i}^{2}+V(i)$ (11.111)

and

 $\omega(i,j)=\frac{1}{r_{ij}}$ (11.112)

The evaluation of matrix elements can now proceed:

 $S_{ab}=\langle\Psi_{b}|\Psi_{a}\rangle=\det(\mathbf{U})\det(\mathbf{V}^{% \dagger})\prod_{i=1}^{N}\hat{s}_{ii}$ (11.113)
 $\Omega^{(1)}_{ab}=\langle\Psi_{b}|\Omega^{(1)}|\Psi_{a}\rangle=\det(\mathbf{U}% )\det(\mathbf{V}^{\dagger})\sum^{N}_{i=1}\langle\hat{b}_{i}|\omega(1)|\hat{a}_% {i}\rangle\cdot\prod^{N}_{j\neq i}\hat{s}_{jj}$ (11.114)
 $\Omega^{(2)}_{ab}=\langle\Psi_{b}|\Omega^{(2)}|\Psi_{a}\rangle=\frac{1}{2}\det% (\mathbf{U})\det(\mathbf{V}^{\dagger})\sum^{N}_{ij}\langle\hat{b}_{i}\hat{b}_{% j}|\omega(1,2)(1-P_{12})|\hat{a}_{i}\hat{a}_{j}\rangle\cdot\prod^{N}_{k\neq{i,% j}}\hat{s}_{kk}$ (11.115)
 $H_{ab}=\Omega^{(1)}_{ab}+\Omega^{(2)}_{ab}$ (11.116)

In an atomic orbital basis set, $\{\chi\}$, we can expand the molecular spin orbitals $\mathbf{a}$ and $\mathbf{b}$,

 $\displaystyle\mathbf{a}=\chi\mathbf{A},\qquad\hat{\mathbf{a}}=\chi\mathbf{A}% \mathbf{V}=\chi\hat{\mathbf{A}}$ (11.117) $\displaystyle\mathbf{b}=\chi\mathbf{B},\qquad\hat{\mathbf{b}}=\chi\mathbf{B}% \mathbf{U}=\chi\hat{\mathbf{B}}$ (11.118)

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

 $\displaystyle\Omega^{(1)}_{ab}$ $\displaystyle=$ $\displaystyle\sum^{N}_{i}\sum_{\lambda\sigma}\hat{A}_{\lambda i}T_{ii}\hat{B}^% {\dagger}_{i\sigma}\langle\chi_{\sigma}|\omega(1)|\chi_{\lambda}\rangle$ (11.119) $\displaystyle=$ $\displaystyle\sum_{\lambda\sigma}G_{\lambda\sigma}\omega_{\sigma\lambda}$

where $\displaystyle T_{ii}=S_{ab}/\hat{s}_{ii}$ and define a generalized density matrix, $\mathbf{G}$:

 $\displaystyle\mathbf{G}$ $\displaystyle=$ $\displaystyle\hat{\mathbf{A}}\mathbf{T}\hat{\mathbf{B}}^{\dagger}$ (11.120)

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

 $\displaystyle\Omega^{(2)}_{ab}$ $\displaystyle=$ $\displaystyle\frac{1}{2}\sum_{ij}\sum_{\lambda\sigma}\sum_{\mu\nu}\hat{A}_{% \lambda i}\hat{A}_{\sigma j}\left(\frac{1}{\hat{s}_{ii}}\right)T_{jj}\hat{B}^{% \dagger}_{i\mu}\hat{B}^{\dagger}_{j\nu}\langle\chi_{\mu}\chi_{\nu}|\omega(1,2)% |\chi_{\lambda}\chi_{\sigma}\rangle$ (11.121) $\displaystyle=$ $\displaystyle\sum_{\lambda\sigma\mu\nu}G^{L}_{\lambda\mu}G^{R}_{\sigma\nu}% \langle\mu\nu||\lambda\sigma\rangle$

where $\mathbf{G}^{R}$ and $\mathbf{G}^{L}$ are generalized density matrices as defined in Eq. (11.120) except $T_{ii}$ in $\mathbf{G}^{L}$ is replaced by $1/(2s_{ii})$.

The $\alpha$- and $\beta$-spin orbitals are treated explicitly. In terms of the spatial orbitals, the one- and two-electron contributions can be reduced to

 $\displaystyle\Omega^{(1)}_{ab}$ $\displaystyle=$ $\displaystyle\sum_{\lambda\sigma}G^{\alpha}_{\lambda\sigma}\omega_{\sigma% \lambda}+\sum_{\lambda\sigma}G^{\beta}_{\lambda\sigma}\omega_{\sigma\lambda}$ (11.122) $\displaystyle\Omega^{(2)}_{ab}$ $\displaystyle=$ $\displaystyle\sum_{\lambda\sigma\mu\nu}G^{L\alpha}_{\lambda\mu}G^{R\alpha}_{% \sigma\nu}(\langle\mu\nu|\lambda\sigma\rangle-\langle\mu\nu|\sigma\lambda% \rangle)+\sum_{\lambda\sigma\mu\nu}G^{L\beta}_{\lambda\mu}G^{R\alpha}_{\sigma% \nu}\langle\mu\nu|\lambda\sigma\rangle$ (11.123) $\displaystyle+\sum_{\lambda\sigma\mu\nu}G^{L\alpha}_{\lambda\mu}G^{R\beta}_{% \sigma\nu}\langle\mu\nu|\lambda\sigma\rangle+\sum_{\lambda\sigma\mu\nu}G^{L% \beta}_{\lambda\mu}G^{R\beta}_{\sigma\nu}(\langle\mu\nu|\lambda\sigma\rangle-% \langle\mu\nu|\sigma\lambda\rangle)$

The resulting one- and two-electron contributions, Eqs. (11.122) and (11.123) 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.97), 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.101).

Example 11.35  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.36 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.929, 1027 FMO coupling value is essentially an off-diagonal Kohn–Sham matrix element with the overlap effect accounted  $V^{\mathrm{FMO}}=\frac{f_{\mathrm{DA}}-S\left(f_{\mathrm{DD}}+f_{\mathrm{AA}}% \right)/2}{1-S^{2}}$ (11.124) where $f_{\mathrm{DA}}=\langle\phi^{\mathrm{D}}_{\mathrm{FMO}}|\hat{f}|\phi^{\mathrm{% A}}_{\mathrm{FMO}}\rangle$, with $\hat{f}$ being the Kohn–Sham operator of the donor-acceptor system. $\phi^{\mathrm{D}(\mathrm{A})}_{\mathrm{FMO}}$ 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.929, 1027 FMO coupling value is essentially an off-diagonal Kohn–Sham matrix element with the overlap effect accounted  $V^{\mathrm{FMO}}=\frac{f_{\mathrm{DA}}-S\left(f_{\mathrm{DD}}+f_{\mathrm{AA}}% \right)/2}{1-S^{2}}$ (11.125) where $f_{\mathrm{DA}}=\langle\phi^{\mathrm{D}}_{\mathrm{FMO}}|\hat{f}|\phi^{\mathrm{% A}}_{\mathrm{FMO}}\rangle$, with $\hat{f}$ being the Kohn–Sham operator of the donor-acceptor system. $\phi^{\mathrm{D}(\mathrm{A})}_{\mathrm{FMO}}$ 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.37 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.