10.15.2 Diabatic-State-Based Methods

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}⟶{\mathrm{DA}}^{-}$$

(10.81)

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

(10.82)

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

$${𝐇}_{\text{eff}}={𝐒}^{-1/2}{\mathrm{𝐇𝐒}}^{-1/2},$$

(10.83)

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

$$V=H_{\mathrm{if}}^{\text{eff}}=\frac{H_{\mathrm{if}}-S_{\mathrm{if}}(H_{\mathrm{ii}}+H_{\mathrm{ff}})/2}{1-S_{\mathrm{if}}^2}$$

(10.84)

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_{\mathrm{if}}^2}$$

(10.85)

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. (10.85), and it is widely used to calculate charge transfer coupling.^{707, 104, 251, 1098} 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.^{707, 104, 688} 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 transformation⁴⁷⁶ and a calculation for Slater’s determinants in $H_{\mathrm{if}}$ and $S_{\mathrm{if}}$.^{251, 1098} 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}}$.¹⁰⁷³ To calculate charge transfer coupling with DFT, we can use the CDFT-CI method (Section 5.13.3), the frontier molecular orbital (FMO) approach^{980, 1089} (Section 10.15.2.5) or a hybrid scheme – DC with CDFT wave functions (Section 10.15.2.4).

Let $|{\mathrm{\Psi }}_a⟩$ and $|{\mathrm{\Psi }}_b⟩$ be two single Slater-determinant wave functions for the initial and final states, and $𝐚$ and $𝐛$ be the spin-orbital sets, respectively:

	$𝐚$	$=$	$(a_1,a_2,\mathrm{\cdots },a_N)$		(10.86)
	$𝐛$	$=$	$(b_1,b_2,\mathrm{\cdots },b_N)$		(10.87)

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

$$𝐒=\int {𝐛}^{†}𝐚𝑑\tau .$$

(10.88)

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 $H_{ab}$ and $S_{ab}$, two sets of new orthogonal spin-orbitals can be used by the corresponding orbital transformation.⁴⁷⁶ In this approach, each set of spin-orbitals $𝐚$ and $𝐛$ are linearly transformed,

	$\widehat{𝐚}$	$=$	$\mathrm{𝐚𝐕}$		(10.89)
	$\widehat{𝐛}$	$=$	$\mathrm{𝐛𝐔}$		(10.90)

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

$$𝐒=𝐔\widehat{𝐬}{𝐕}^{†}$$

(10.91)

The overlap matrix in the new basis is now diagonal

$$\int {\widehat{𝐛}}^{†}\widehat{𝐚}={𝐔}^{†}\left(\int {𝐛}^{†}𝐚\right)𝐕=\widehat{𝐬}$$

(10.92)

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

	${\mathrm{\Omega }}^{(1)}$	$=$	${\displaystyle \sum _{i=1}^N}\omega (i)$		(10.93)
	${\mathrm{\Omega }}^{(2)}$	$=$	${\displaystyle \frac{1}{2}}{\displaystyle \sum _{i,j=1}^N}\omega (i,j)$		(10.94)

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

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

(10.95)

and

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

(10.96)

The evaluation of matrix elements can now proceed:

$$S_{ab}=⟨{\mathrm{\Psi }}_b|{\mathrm{\Psi }}_a⟩=det(𝐔)det({𝐕}^{†})\prod _{i=1}^N{\widehat{s}}_{ii}$$

(10.97)

$${\mathrm{\Omega }}_{ab}^{(1)}=⟨{\mathrm{\Psi }}_b|{\mathrm{\Omega }}^{(1)}|{\mathrm{\Psi }}_a⟩=det(𝐔)det({𝐕}^{†})\sum _{i=1}^N⟨{\widehat{b}}_i|\omega (1)|{\widehat{a}}_i⟩\cdot \prod _{j\ne i}^N{\widehat{s}}_{jj}$$

(10.98)

$${\mathrm{\Omega }}_{ab}^{(2)}=⟨{\mathrm{\Psi }}_b|{\mathrm{\Omega }}^{(2)}|{\mathrm{\Psi }}_a⟩=\frac{1}{2}det(𝐔)det({𝐕}^{†})\sum _{ij}^N⟨{\widehat{b}}_i{\widehat{b}}_j|\omega (1,2)(1-P_{12})|{\widehat{a}}_i{\widehat{a}}_j⟩\cdot \prod _{k\ne i,j}^N{\widehat{s}}_{kk}$$

(10.99)

$$H_{ab}={\mathrm{\Omega }}_{ab}^{(1)}+{\mathrm{\Omega }}_{ab}^{(2)}.$$

(10.100)

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

	$𝐚=\chi 𝐀,\widehat{𝐚}=\chi \mathrm{𝐀𝐕}=\chi \widehat{𝐀}$		(10.101)
	$𝐛=\chi 𝐁,\widehat{𝐛}=\chi \mathrm{𝐁𝐔}=\chi \widehat{𝐁}$		(10.102)

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

	${\mathrm{\Omega }}_{ab}^{(1)}$	$=$	${\displaystyle \sum _i^N}{\displaystyle \sum _{\lambda \sigma }}{\widehat{A}}_{\lambda i}{T}_{ii}{\widehat{B}}_{i\sigma }^{†}⟨{\chi }_{\sigma }|\omega (1)|{\chi }_{\lambda }⟩$		(10.103)
		$=$	${\displaystyle \sum _{\lambda \sigma }}{G}_{\lambda \sigma }{\omega }_{\sigma \lambda }$

where $T_{ii}=S_{ab}/{\widehat{s}}_{ii}$ and define a generalized density matrix, $𝐆$:

$$𝐆=\widehat{𝐀}𝐓{\widehat{𝐁}}^{†}.$$

(10.104)

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

	${\mathrm{\Omega }}_{ab}^{(2)}$	$=$	${\displaystyle \frac{1}{2}}{\displaystyle \sum _{ij}}{\displaystyle \sum _{\lambda \sigma }}{\displaystyle \sum _{\mu \nu }}{\widehat{A}}_{\lambda i}{\widehat{A}}_{\sigma j}\left({\displaystyle \frac{1}{{\widehat{s}}_{ii}}}\right)T_{jj}{\widehat{B}}_{i\mu }^{†}{\widehat{B}}_{j\nu }^{†}⟨{\chi }_{\mu }{\chi }_{\nu }|\omega (1,2)|{\chi }_{\lambda }{\chi }_{\sigma }⟩$		(10.105)
		$=$	${\displaystyle \sum _{\lambda \sigma \mu \nu }}{G}_{\lambda \mu }^L{G}_{\sigma \nu }^R⟨\mu \nu ||\lambda \sigma ⟩$

where ${𝐆}^R$ and ${𝐆}^L$ are generalized density matrices as defined in Eq. (10.104) except $T_{ii}$ in ${𝐆}^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

	${\mathrm{\Omega }}_{ab}^{(1)}$	$=$	${\displaystyle \sum _{\lambda \sigma }}{G}_{\lambda \sigma }^{\alpha }{\omega }_{\sigma \lambda }+{\displaystyle \sum _{\lambda \sigma }}{G}_{\lambda \sigma }^{\beta }{\omega }_{\sigma \lambda }$		(10.106)
	${\mathrm{\Omega }}_{ab}^{(2)}$	$=$	${\displaystyle \sum _{\lambda \sigma \mu \nu }}{G}_{\lambda \mu }^{L\alpha }{G}_{\sigma \nu }^{R\alpha }(⟨\mu \nu |\lambda \sigma ⟩-⟨\mu \nu |\sigma \lambda ⟩)+{\displaystyle \sum _{\lambda \sigma \mu \nu }}{G}_{\lambda \mu }^{L\beta }{G}_{\sigma \nu }^{R\alpha }⟨\mu \nu |\lambda \sigma ⟩$		(10.107)
			$+{\displaystyle \sum _{\lambda \sigma \mu \nu }}{G}_{\lambda \mu }^{L\alpha }{G}_{\sigma \nu }^{R\beta }⟨\mu \nu |\lambda \sigma ⟩+{\displaystyle \sum _{\lambda \sigma \mu \nu }}{G}_{\lambda \mu }^{L\beta }{G}_{\sigma \nu }^{R\beta }(⟨\mu \nu |\lambda \sigma ⟩-⟨\mu \nu |\sigma \lambda ⟩)$

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

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 (10.81), 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. (10.85).

  Example 10.40  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 10.41  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
  

The frontier molecular orbital (FMO) approach is often used with DFT to calculate ET coupling.^{980, 1089} 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}$$

(10.108)

where $f_{\mathrm{DA}}=⟨{\varphi }_{\mathrm{FMO}}^{\mathrm{D}}|\widehat{f}|{\varphi }_{\mathrm{FMO}}^{\mathrm{A}}⟩$, with $\widehat{f}$ being the Kohn–Sham operator of the donor-acceptor system. ${\varphi }_{\mathrm{FMO}}^{\mathrm{D}(\mathrm{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 has a setup that is similar to the “1+1” 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.