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}\u27f6{\mathrm{DA}}^{-}$$ | (10.76) |
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
$$\mathbf{H}-E\mathbf{S}=\left(\begin{array}{cc}\hfill {H}_{\mathrm{ii}}-{S}_{\mathrm{ii}}E\hfill & \hfill {H}_{\mathrm{if}}-{S}_{\mathrm{if}}E\hfill \\ \hfill {H}_{\mathrm{if}}-{S}_{\mathrm{if}}E\hfill & \hfill {H}_{\mathrm{ff}}-{S}_{\mathrm{ff}}E\hfill \end{array}\right)=0$$ | (10.77) |
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.77) can be obtained by using the Löwdin transformation:
$${\mathbf{H}}_{\text{eff}}={\mathbf{S}}^{-1/2}{\mathrm{\mathbf{H}\mathbf{S}}}^{-1/2},$$ | (10.78) |
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.79) |
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.80) |
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.80), and it is widely used to calculate charge transfer coupling.^{Ohta:1986, Broo:1990, Farazdel:1990, Zhang:1997} 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.^{Ohta:1986, Broo:1990, Newton:1991} 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^{King:1967} and a calculation for Slater’s determinants in ${H}_{\mathrm{if}}$ and ${S}_{\mathrm{if}}$.^{Farazdel:1990, Zhang:1997} 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}}$.^{Wu:2006c} To calculate charge transfer coupling with DFT, we can use the CDFT-CI method (Section 5.13.3), the frontier molecular orbital (FMO) approach^{Valeev:2006, You:2015a} (Section 10.15.2.5) or a hybrid scheme – DC with CDFT wave functions (Section 10.15.2.4).
Let $|{\mathrm{\Psi}}_{a}\u27e9$ and $|{\mathrm{\Psi}}_{b}\u27e9$ 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:
$\mathbf{a}$ | $=$ | $({a}_{1},{a}_{2},\mathrm{\cdots},{a}_{N})$ | (10.81) | ||
$\mathbf{b}$ | $=$ | $({b}_{1},{b}_{2},\mathrm{\cdots},{b}_{N})$ | (10.82) |
Since the two sets of spin-orbitals are not orthogonal, the overlap matrix $\mathbf{S}$ can be defined as:
$$\mathbf{S}=\int {\mathbf{b}}^{\u2020}\mathbf{a}\mathit{d}\tau .$$ | (10.83) |
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.^{King:1967} In this approach, each set of spin-orbitals $\mathbf{a}$ and $\mathbf{b}$ are linearly transformed,
$\widehat{\mathbf{a}}$ | $=$ | $\mathrm{\mathbf{a}\mathbf{V}}$ | (10.84) | ||
$\widehat{\mathbf{b}}$ | $=$ | $\mathrm{\mathbf{b}\mathbf{U}}$ | (10.85) |
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}\widehat{\mathbf{s}}{\mathbf{V}}^{\u2020}$$ | (10.86) |
The overlap matrix in the new basis is now diagonal
$$\int {\widehat{\mathbf{b}}}^{\u2020}\widehat{\mathbf{a}}={\mathbf{U}}^{\u2020}\left(\int {\mathbf{b}}^{\u2020}\mathbf{a}\right)\mathbf{V}=\widehat{\mathbf{s}}$$ | (10.87) |
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)}$ | $=$ | $\sum _{i=1}^{N}}\omega (i)$ | (10.88) | ||
${\mathrm{\Omega}}^{(2)}$ | $=$ | $\frac{1}{2}}{\displaystyle \sum _{i,j=1}^{N}}\omega (i,j)$ | (10.89) |
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.90) |
and
$$\omega (i,j)=\frac{1}{{r}_{ij}}.$$ | (10.91) |
The evaluation of matrix elements can now proceed:
$${S}_{ab}=\u27e8{\mathrm{\Psi}}_{b}|{\mathrm{\Psi}}_{a}\u27e9=det(\mathbf{U})det({\mathbf{V}}^{\u2020})\prod _{i=1}^{N}{\widehat{s}}_{ii}$$ | (10.92) |
$${\mathrm{\Omega}}_{ab}^{(1)}=\u27e8{\mathrm{\Psi}}_{b}|{\mathrm{\Omega}}^{(1)}|{\mathrm{\Psi}}_{a}\u27e9=det(\mathbf{U})det({\mathbf{V}}^{\u2020})\sum _{i=1}^{N}\u27e8{\widehat{b}}_{i}|\omega (1)|{\widehat{a}}_{i}\u27e9\cdot \prod _{j\ne i}^{N}{\widehat{s}}_{jj}$$ | (10.93) |
$${\mathrm{\Omega}}_{ab}^{(2)}=\u27e8{\mathrm{\Psi}}_{b}|{\mathrm{\Omega}}^{(2)}|{\mathrm{\Psi}}_{a}\u27e9=\frac{1}{2}det(\mathbf{U})det({\mathbf{V}}^{\u2020})\sum _{ij}^{N}\u27e8{\widehat{b}}_{i}{\widehat{b}}_{j}|\omega (1,2)(1-{P}_{12})|{\widehat{a}}_{i}{\widehat{a}}_{j}\u27e9\cdot \prod _{k\ne i,j}^{N}{\widehat{s}}_{kk}$$ | (10.94) |
$${H}_{ab}={\mathrm{\Omega}}_{ab}^{(1)}+{\mathrm{\Omega}}_{ab}^{(2)}.$$ | (10.95) |
In an atomic orbital basis set, $\{\chi \}$, we can expand the molecular spin orbitals $\mathbf{a}$ and $\mathbf{b}$,
$\mathbf{a}=\chi \mathbf{A},\widehat{\mathbf{a}}=\chi \mathrm{\mathbf{A}\mathbf{V}}=\chi \widehat{\mathbf{A}}$ | (10.96) | ||
$\mathbf{b}=\chi \mathbf{B},\widehat{\mathbf{b}}=\chi \mathrm{\mathbf{B}\mathbf{U}}=\chi \widehat{\mathbf{B}}$ | (10.97) |
The one-electron terms, Eq. (10.92), can be expressed as
${\mathrm{\Omega}}_{ab}^{(1)}$ | $=$ | $\sum _{i}^{N}}{\displaystyle \sum _{\lambda \sigma}}{\widehat{A}}_{\lambda i}{T}_{ii}{\widehat{B}}_{i\sigma}^{\u2020}\u27e8{\chi}_{\sigma}|\omega (1)|{\chi}_{\lambda}\u27e9$ | (10.98) | ||
$=$ | $\sum _{\lambda \sigma}}{G}_{\lambda \sigma}{\omega}_{\sigma \lambda$ |
where ${T}_{ii}={S}_{ab}/{\widehat{s}}_{ii}$ and define a generalized density matrix, $\mathbf{G}$:
$$\mathbf{G}=\widehat{\mathbf{A}}\mathbf{T}{\widehat{\mathbf{B}}}^{\u2020}.$$ | (10.99) |
Similarly, the two-electron terms, Eq. (10.94), are
${\mathrm{\Omega}}_{ab}^{(2)}$ | $=$ | $\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}^{\u2020}{\widehat{B}}_{j\nu}^{\u2020}\u27e8{\chi}_{\mu}{\chi}_{\nu}|\omega (1,2)|{\chi}_{\lambda}{\chi}_{\sigma}\u27e9$ | (10.100) | ||
$=$ | $\sum _{\lambda \sigma \mu \nu}}{G}_{\lambda \mu}^{L}{G}_{\sigma \nu}^{R}\u27e8\mu \nu ||\lambda \sigma \u27e9$ |
where ${\mathbf{G}}^{R}$ and ${\mathbf{G}}^{L}$ are generalized density matrices as defined in Eq. (10.99) except ${T}_{ii}$ in ${\mathbf{G}}^{L}$ is replaced by $1/(2{s}_{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)}$ | $=$ | $\sum _{\lambda \sigma}}{G}_{\lambda \sigma}^{\alpha}{\omega}_{\sigma \lambda}+{\displaystyle \sum _{\lambda \sigma}}{G}_{\lambda \sigma}^{\beta}{\omega}_{\sigma \lambda$ | (10.101) | ||
${\mathrm{\Omega}}_{ab}^{(2)}$ | $=$ | $\sum _{\lambda \sigma \mu \nu}}{G}_{\lambda \mu}^{L\alpha}{G}_{\sigma \nu}^{R\alpha}(\u27e8\mu \nu |\lambda \sigma \u27e9-\u27e8\mu \nu |\sigma \lambda \u27e9)+{\displaystyle \sum _{\lambda \sigma \mu \nu}}{G}_{\lambda \mu}^{L\beta}{G}_{\sigma \nu}^{R\alpha}\u27e8\mu \nu |\lambda \sigma \u27e9$ | (10.102) | ||
$+{\displaystyle \sum _{\lambda \sigma \mu \nu}}{G}_{\lambda \mu}^{L\alpha}{G}_{\sigma \nu}^{R\beta}\u27e8\mu \nu |\lambda \sigma \u27e9+{\displaystyle \sum _{\lambda \sigma \mu \nu}}{G}_{\lambda \mu}^{L\beta}{G}_{\sigma \nu}^{R\beta}(\u27e8\mu \nu |\lambda \sigma \u27e9-\u27e8\mu \nu |\sigma \lambda \u27e9)$ |
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.76), 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.80).
$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.
$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.^{Valeev:2006, You:2015a} 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.103) |
where ${f}_{\mathrm{DA}}=\u27e8{\varphi}_{\mathrm{FMO}}^{\mathrm{D}}|\widehat{f}|{\varphi}_{\mathrm{FMO}}^{\mathrm{A}}\u27e9$, 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.