10.15 Electronic Couplings for Electron- and Energy Transfer

10.15.3 Fragment-Based Methods for Electronic Coupling

10.15.3.1 Approach based on absolutely localized molecular orbitals

One can use absolutely localized molecular orbitals (ALMOs, see Chapter 12) to construct charge-localized diabatic states directly from DFT calculations. The ALMOs on each fragment are expanded by the AO basis functions belonging to the same fragment alone, whose corresponding MO coefficient matrix is fragment block-diagonal.Khaliullin:2006 In energy decomposition analysis methods,Khaliullin:2007, Horn:2016c ALMOs are utilized to separate the effects of polarization and charge transfer in intermolecular binding, because they have the useful property that they do not allow for charge transfer between fragments under the Mulliken definition of charge population. Making use of this property, one can construct charge-localized diabats for hole and electron transfer. For example, considering the initial and final states of a hole transfer process, |D+A and |DA+, the two diabats can be represented in the following form:

|ψa =1(N-1)!det{ϕD1(a),ϕD2(a),,ϕDnD-1(a)ϕA1(a),ϕA2(a),,ϕAnA(a)} (10.104a)
|ψb =1(N-1)!det{ϕD1(b),ϕD2(b),,ϕDnD(b)ϕA1(b),ϕA2(b),,ϕAnA-1(b)} (10.104b)

For systems where the donor and acceptor moieties are well-separated, one can construct the ALMO-based diabats by simply concatenating orbitals obtained from isolated fragment calculations: D+ and A for one diabat, and D and A+ for the other. The energy of each ALMO diabat can then be variationally optimized with respect to orbital rotations on fragment, using the SCF-MI technique (see Section 12.4).Stoll:1980, Gianinetti:1996, Khaliullin:2006 These ALMO-based diabatic states are variationally optimized such that the associated nuclear forces can be easily computed.Mao:2017 The mutual polarization of donor and acceptor moieties in the presence of each other is also taken into account.

To calculate the electronic coupling between two ALMO diabats, one should first construct the diabatic Hamiltonian in the ALMO state basis

𝐇=(HaaHabHbaHbb) (10.105)

and then transform that into the Löwdin-orthogonalized basis

𝐇=𝐒-1/2𝐇𝐒-1/2 (10.106)

whose off-diagonal element, Hab, corresponds to the diabatic coupling to be evaluated. In the 2-state case, we have

Hab=11-Sab2|Hab-Haa+Hbb2Sab| (10.107)

which requires the overlap between two ALMO diabats and the diagonal and off-diagonal elements of 𝐇. The interstate overlap is given by

Sab=ψa|ψb=det[(𝐂o(a))𝐒𝐂o(b)]. (10.108)

where 𝐂o(a) and 𝐂o(b) are MO coefficients for the occupied orbitals in diabats |ψa and |ψb, respectively, and 𝐒 is the AO overlap matrix.

The elements of the diabatic Hamiltonian matrix can be evaluated using the multi-state DFT (MSDFT) approach.Cembran:2009, Ren:2016, Mao:2019 For the diagonal elements, it is straightforward to employ the KS energies of the two diabats:

Haa=EaKS[𝐏(a)],Hbb=EbKS[𝐏(b)] (10.109)

where 𝐏(a) and 𝐏(b) are the one-electron density matrices associated with two ALMO states |ψa and |ψb, respectively. The approximation for the off-diagonal element is theoretically more challenging. In the original MSDFT scheme,Cembran:2009, Ren:2016

Hab=Sab[Vnn+𝐏ab𝐡+12𝐏ab𝐈𝐈𝐏ab+12(ΔEac+ΔEbc)] (10.110)

where 𝐏ab is the one-particle transition density matrix between two ALMO states

𝐏ab=𝐂o(a)[(𝐂o(b))𝐒𝐂o(a)]-1(𝐂o(b)) (10.111)

The first three terms on the right-hand side of Eq. (10.110) correspond to the contributions from nuclear repulsion, one-electron Hamiltonian (kinetic energy and nuclei-electron attraction), and full two-electron integrals (Coulomb and full HF exchange), which can be derived as in non-orthogonal CI.Thom:2009b The last term accounts for the contribution from exchange-correlation (XC) functional as a correction to the HF coupling, which is given by the average of the difference between the KS and HF energies calculated from the same one-electron density matrix for each diabat:

ΔEac =EaKS[𝐏(a)]-EaHF[𝐏(a)] (10.112a)
ΔEbc =EbKS[𝐏(b)]-EbHF[𝐏(b)]. (10.112b)

This approach was denoted as ALMO(MSDFT) in Ref. Mao:2019 and it was found to overestimate the electronic couplings for the tested hole and electron transfer systems. A modified approach, denoted as ALMO(MSDFT2), was proposed in Ref. Mao:2019, which evaluates the XC contribution using the XC energy of the symmetrized transition density matrix

Hab=Sab[Vnn+𝐏ab𝐡+12𝐏ab𝐈𝐈𝐏ab+Exc[𝐏~ab]] (10.113)

where

𝐏~ab=12(𝐏ab+𝐏ba). (10.114)

Note that in Eq. (10.113), 𝐈𝐈 includes only Coulomb integrals and a fraction of exact exchange if hybrid functionals are employed.

According to the benchmark results in Ref. Mao:2019, ALMO(MSDFT2) shows better accuracy than the original MSDFT method for hole and electron transfer, and thus it is implemented as the default approach to compute electronic couplings between ALMO diabats in Q-Chem. We note that the results given by Eq. (10.113) may become inaccurate when the overlap between two states becomes near-singular, as

𝝈ba=(𝐂o(b))𝐒𝐂o(a) (10.115)

is inverted when constructing the transition density [Eq. (10.111)]. To circumvent this numerical issue, one can replace the inverse in Eq. (10.111) with the Penrose pseudo-inverse, which was suggested for a similar objective in Ref. Pavanello:2013.

10.15.3.2 Other fragment based diabatization methods

Besides ALMO-based diabatization method, we also implemented two other fragment-based approaches: the projection operator diabatization (POD) and the fragment orbital DFT (FODFT) approach. The POD method Kondov:2007 starts from a standard KS-DFT calculation of the system and post-processes the converged Fock matrix. It first transforms the Fock matrix into the Löwdin-orthogonalized AO basis and then partitions that into the donor and acceptor blocks, assuming that these orthogonalized AO basis functions still retain their original fragment tags:

𝐅~=𝐒-1/2𝐅𝐒-1/2=(𝐅~dd𝐅~da𝐅~ad𝐅~aa) (10.116)

One then diagonalizes 𝐅~dd and 𝐅~aa separately

ϵd=𝐃d𝐅~dd𝐃d,ϵa=𝐃a𝐅~aa𝐃a, (10.117)

where the eigenvectors 𝐃d and 𝐃a define the single-particle “diabatic states”:

|φ¯p(d)=μ|χ~μ(d)(Dd)pμ|φ¯p(a)=μ|χ~μ(a)(Da)pμ, (10.118)

and transforms the off-diagonal block of the Fock matrix into this diabatic basis

𝐅¯da=𝐃d𝐅~da𝐃a (10.119)

yielding

𝐅¯=(ϵd𝐅¯da𝐅¯adϵa) (10.120)

The couplings between these single-particle diabats can then be directly read off from the elements of 𝐅¯da.

The Q-Chem implementation of the POD method follows the description in Refs. Kondov:2007 and Yang:2018b, where a closed-shell reference system is used to generate the Fock matrix to be processed, i.e., 𝐅 in Eq. (10.116).

Fragment orbital DFT (FODFT) Oberhofer:2012, Senthilkumar:2003, Schober:2016 is an approach to compute the diabatic couplings for hole and electron transfer between fragments. There have been several different flavors of FODFT approaches developed in literature, and here we introduce the most recent variant by Schober et al.Schober:2016 Considering a hole transfer process D++AD+A+ or an electron transfer process D-+AD+A-, where the donor (D) and acceptor (A) fragments have nD and nA electrons, respectively, the procedure is as follows:

  • Perform KS-DFT calculations for isolated donor and acceptor fragments; collect the converged fragment orbitals:
    {ϕD1,ϕD2,,ϕDnD±1} and {ϕA1,ϕA2,,ϕAnA}

  • Löwdin-orthogonalize the occupied orbitals on two fragments. The reactant diabat (D+A or D-A) can be represented as

    |ψ¯a=1(N-1)!det{ϕ¯D1,ϕ¯D2,,ϕ¯DnD±1ϕ¯A1,ϕ¯A2,,ϕ¯AnA} (10.121)

    where “ϕ¯” denotes Löwdin-orthogonalized orbitals, and N=nD+nA. Note that the lowest unoccupied orbital where the electron is transferring to, ϕDnD in the case of HT or ϕAnA+1 in the case of ET, also needs to be made orthogonal to the space spanned by all occupied orbitals.

  • Construct the product diabat (DA+ or DA-), simply by moving the hole from ϕ¯DnD to ϕ¯AnA (HT), or the excess electron from ϕ¯DnD+1 to ϕ¯AnA+1 (ET)

    |ψ¯b=1(N-1)!det{ϕ¯D1,ϕ¯D2,,ϕ¯DnDϕ¯A1,ϕ¯A2,,ϕ¯AnA±1} (10.122)
  • Compute the electronic coupling between |ψ¯a and |ψ¯b, which is approximated by the coupling of the orthogonalized fragment orbitals through the Kohn-Sham Fock operator (built from the reactant diabat)

    ψ¯a|H^|ψ¯b{ϕDnD|f^KS|ϕAnA,HTϕDnD+1|f^KS|ϕAnA+1,ET (10.123)

The approach described above is denoted as FODFT(2n-1)@D+A (HT) / FODFT(2n+1)@D-A (ET) Schober:2016 as the charged fragment is explicitly taken into account when preparing the fragment orbitals and the KS Fock matrix is built from 2n1 occupied orbitals. Besides this, there are two other variants of FODFT:

  1. 1.

    FODFT(2n)@DA:Senthilkumar:2003 fragment orbitals prepared with D and A both closed-shell; KS Fock operator constructed from 2n occupied orbitals

  2. 2.

    FODFT(2n-1)@DA (HT) / FODFT(2n+1)@D-A- (ET): Oberhofer:2012 fragment orbitals prepared with the system having one excess electron (DA for HT and D-A- for ET), while one occupied orbital is removed when building the KS Fock operator

According to the benchmark results,Schober:2016, Mao:2019 FODFT(2n-1)@D+A (HT) / FODFT(2n+1)@D-A (ET) is the best-performing method, possibly because of its explicit account for charged fragments and consistent electron count in the preparation of fragment orbitals and in the construction of Fock matrix.

One issue associated with the FODFT methods is that for asymmetric systems, the results would depend on how one chooses the initial and final states for an electron or hole transfer process (e.g. D+A vs. DA+), especially for the two variants that build the Fock matrix with 2n±1 occupied orbitals. Mao:2019 The Q-Chem implementation of FODFT(2n-1)@DA / FODFT(2n+1)@D-A- automatically computes Hab in both ways and then reports the average, as it only requires an extra Fock matrix build. This, however, is not automatically done for FODFT(2n-1)@D+A / FODFT(2n+1)@D-A.

10.15.3.3 Job control of fragment based diabatization methods

POD, FODFT, and ALMO(MSDFT) calculations in Q-Chem require specification of fragments in the $molecule section (see Sec. 12.2). For ALMO(MSDFT) calculations, one also needs to specify the charge and multiplicity of each fragment in each diabatic state in the $almo_coupling section, where two hyphens indicate the separation of different diabats:

¯$almo_coupling
¯¯¯charge_frag_1     mult_frag_1        !diabat 1
¯¯¯charge_frag_2     mult_frag_2
¯¯¯--
¯¯¯charge_frag_1     mult_frag_1        !diabat 2
¯¯¯charge_frag_2     mult_frag_2
¯$end
¯

The current implementation of FODFT is limited to hole transfer between the HOMOs of two fragments or electron transfer between the LUMOs, and the implementation of ALMO(MSDFT) is limited to ground state electron/hole transfer involving two states. We will further generalize the implementation of these methods in future releases.

FRAG_DIABAT_METHOD
       Specify fragment based diabatization method
TYPE:
       STRING
DEFAULT:
       NONE
OPTIONS:
       ALMO_MSDFT Perform ALMO(MSDFT) diabatization POD Perform projection operator diabatization ESID The energy-split-in-dimer method,Valeev:2006 which is equivalent to the FMO approach introduced in Sec. 10.15.2.5 FODFT Calculate electronic coupling using fragment orbital DFT
RECOMMENDATION:
       NONE

MSDFT_METHOD
       Specify the scheme for ALMO(MSDFT)
TYPE:
       INTEGER
DEFAULT:
       2
OPTIONS:
       1 The original MSDFT scheme [Eq. (10.110)] 2 The ALMO(MSDFT2) approach [Eq. (10.113)]
RECOMMENDATION:
       Use the default method. Note that the method will be automatically reset to 1 if a meta-GGA functional is requested.

MSDFT_PINV_THRESH
       Set the threshold for pseudo-inverse of the interstate overlap
TYPE:
       INTEGER
DEFAULT:
       4
OPTIONS:
       n Set the threshold to 10-n
RECOMMENDATION:
       Use the default value

FRAG_DIABAT_DOHT
       Specify whether hole or electron transfer is considered
TYPE:
       BOOLEAN
DEFAULT:
       TRUE
OPTIONS:
       TRUE Do hole transfer FALSE Do electron transfer
RECOMMENDATION:
       Need to be specified for POD and FODFT calculations

FRAG_DIABAT_PRINT
       Specify the print level for fragment based diabatization calculations
TYPE:
       INTEGER
DEFAULT:
       0
OPTIONS:
       0 No additional prints 1 Currently it can be used to print out the entire 𝐅¯da in POD
RECOMMENDATION:
       Use 1 if electron/hole transfer between multiple orbital pairs needs to considered in POD

FODFT_METHOD
       Specify the flavor of FODFT method
TYPE:
       INTEGER
DEFAULT:
       1
OPTIONS:
       1 FODFT(2n-1)@D+A (HT) / FODFT(2n+1)@D-A (ET) 2 FODFT(2n)@DA 3 FODFT(2n-1)@DA (HT) / FODFT(2n+1)@D-A- (ET)
RECOMMENDATION:
       The default approach shows the best overall performance

FODFT_DONOR
       Specify the donor fragment in FODFT calculation
TYPE:
       INTEGER
DEFAULT:
       1
OPTIONS:
       1 First fragment as donor 2 Second fragment as donor
RECOMMENDATION:
       With FODFT_METHOD = 1, the charged fragment needs to be the donor fragment

Example 10.41  ALMO(MSDFT2) calculation for hole transfer in ethylene dimer

$molecule
1 2
--
1 2
  C      0.000000    0.000000    0.000000
  C      1.332000    0.000000    0.000000
  H     -0.574301    0.000000   -0.928785
  H     -0.574301    0.000000    0.928785
  H      1.906301    0.000000    0.928785
  H      1.906301    0.000000   -0.928785
--
0 1
  C     -0.000000    4.000000    0.000000
  C      1.332000    4.000000   -0.000000
  H     -0.574301    4.000000    0.928785
  H     -0.574301    4.000000   -0.928785
  H      1.906301    4.000000   -0.928785
  H      1.906301    4.000000    0.928785
$end

$rem
  JOBTYPE  SP
  METHOD   PBE0
  BASIS    6-31+G(D)
  UNRESTRICTED  TRUE
  THRESH  14
  SCF_CONVERGENCE  8
  SYMMETRY   FALSE
  SYM_IGNORE  TRUE
  SCFMI_MODE   1
  FRGM_METHOD  STOLL
  FRAG_DIABAT_METHOD ALMO_MSDFT
$end

$almo_coupling
  1  2
  0  1
  --
  0  1
  1  2
$end

Example 10.42  POD diabatization method for hole transfer in ethylene dimer

$molecule
0 1
--
0 1
  C      0.000000    0.000000    0.000000
  C      1.332000    0.000000    0.000000
  H     -0.574301    0.000000   -0.928785
  H     -0.574301    0.000000    0.928785
  H      1.906301    0.000000    0.928785
  H      1.906301    0.000000   -0.928785
--
0 1
  C     -0.000000    4.000000    0.000000
  C      1.332000    4.000000   -0.000000
  H     -0.574301    4.000000    0.928785
  H     -0.574301    4.000000   -0.928785
  H      1.906301    4.000000   -0.928785
  H      1.906301    4.000000    0.928785
$end

$rem
  jobtype sp
  frag_diabat_method pod
  method  lrc-wpbeh
  basis   6-31+g(d)
  scf_convergence 8
  thresh 14
  symmetry false
  sym_ignore true
$end

Example 10.43  FODFT(2n-1)@D+A calculation for hole transfer in ethylene dimer

$molecule
1 2
--
1 2
  C      0.000000    0.000000    0.000000
  C      1.332000    0.000000    0.000000
  H     -0.574301    0.000000   -0.928785
  H     -0.574301    0.000000    0.928785
  H      1.906301    0.000000    0.928785
  H      1.906301    0.000000   -0.928785
--
0 1
  C     -0.000000    4.000000    0.000000
  C      1.332000    4.000000   -0.000000
  H     -0.574301    4.000000    0.928785
  H     -0.574301    4.000000   -0.928785
  H      1.906301    4.000000   -0.928785
  H      1.906301    4.000000    0.928785
$end

$rem
  jobtype sp
  method  wb97x-d
  basis   6-31+g(d)
  unrestricted true
  scf_convergence 8
  thresh 14
  symmetry false
  sym_ignore true
  frag_diabat_method fodft
  fodft_method 1
$end