Development of the Fragment Molecular Orbital Method for

Article Cite This: J. Phys. Chem. A XXXX, XXX, XXX−XXX

pubs.acs.org/JPCA

Development of the Fragment Molecular Orbital Method for Calculating Nonlocal Excitations in Large Molecular Systems Takatoshi Fujita*,† and Yuji Mochizuki‡,§ †

Institute for Molecular Science, Okazaki, Aichi 444-0865, Japan Department of Chemistry and Research Center for Smart Molecules, Faculty of Science, Rikkyo University, 3-34-1 Nishi-Ikebukuro, Toshima-ku, Tokyo 171-8501, Japan § Institute for Industrial Science, The University of Tokyo, 4-6-1 Komaba, Meguro-ku, Tokyo 153-8505, Japan ‡

S Supporting Information *

ABSTRACT: We developed the fragment-based method for calculating nonlocal excitations in large molecular systems. This method is based on the multilayer fragment molecular orbital method and the configuration interaction single (CIS) wave function using localized molecular orbitals. The excited-state wave function for the whole system is described as a superposition of configuration state functions (CSFs) for intrafragment excitations and for interfragment charge-transfer excitations. The formulation and calculations of singlet excited-state Hamiltonian matrix elements in the fragment CSFs are presented in detail. The efficient approximation schemes for calculating the matrix elements are also presented. The computational efficiency and the accuracy were evaluated using the molecular dimers and molecular aggregates. We confirmed that absolute errors of 50 meV (relative to the conventional calculations) are achievable for the molecular systems in their equilibrium geometries. The perturbative electron correlation correction to the CIS excitation energies is also demonstrated. The present theory can compute a large number of excited states in large molecular systems; in addition, it allows for the systematic derivation of a model exciton Hamiltonian. These features are useful for studying excitedstate dynamics in condensed molecular systems based on the ab initio electronic structure theory.

1. INTRODUCTION Theoretical studies of electronically excited states are essential for understanding many photophysical and photochemical processes. However, predicting excited states of large systems with reasonable accuracy is still a challenging issue in quantum chemistry. A sophisticated electronic structure theory that can treat an excited state with high accuracy requires considerable computational time, which scales steeply with system size. The highly nonlinear scaling of the electronic structure theory places a restriction on a method that is practically useful for large systems. On the one hand, the time-dependent density functional theory (TDDFT) is widely used for treating large molecules1,2 because of its efficient treatment of dynamical correlations. On the other hand, a starting point of the wave function-based theory for excited states is the configuration interaction single (CIS).1,3,4 The wave function-based theory enables the systematic improvement of the accuracy, such as with the perturbative corrections,5,6 but at additional computational costs. Thanks to the Davidson subspace iterative diagonalization algorithm,7 TDDFT and CIS can compute several excited states with comparable computational times to the ground state. However, obtaining all excited states by TDDFT or CIS requires the computational time that increases asymptotically with the sixth power of system size. For example, © XXXX American Chemical Society

to study an absorption spectrum or excitation energy transfer (EET) in an organic molecular crystal, one must compute the excited state manifold that contains as many states as the number of molecules present. The development of an excitedstate theory suitable for large systems is necessary for efficient computations of such an excited-state manifold. Numerous efforts have focused on the development of a linear-scaling electronic structure theory for treating large systems. One such successful linear-scaling approach is to exploit the locality of the wave function and to adopt a variety of fragmentation, incrementation, or multilayer schemes (see ref 8 and references therein). Some of the fragment-based approaches have been successfully applied to local electronic excitations.9,10 However, it is not straightforward to calculate nonlocal excited states by the fragment-based methods. One type of nonlocal excited state is a collective excited state that is a superposition of local excitations. Other types of nonlocal excited states include charge-transfer (CT) excitation between different subunits and electronic transitions between delocalized molecular orbitals (MOs). Several methods have been Received: January 15, 2018 Revised: March 12, 2018 Published: March 28, 2018 A

DOI: 10.1021/acs.jpca.8b00446 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

(hole) (H(elec) IK ) and hole (HIK ) Hamiltonians and the two-body part (HIJ,KL) that includes excitonic couplings and electron−hole Coulomb interactions. This model includes the CT excitations by its construction. Moreover, the electronic transition between delocalized orbitals can be treated in this model. Because diagonalizing the one-body electron or hole Hamiltonian with translational symmetry results in conduction or valence band, the electronic transition between delocalized orbitals is considered. The two-body electron−hole Hamiltonian was adopted to study the reaction center in photosynthetic systems,38 organic molecular crystals,39,40 and organic/organic interfaces.41,42 The implementation of the electron−hole Hamiltonian requires several parameters including orbital energies, transfer integrals, and electron−hole Coulomb interactions. The TDFI-TI method proposed by Fujimoto and co-workers 15,16 can be regarded as an ab initio implementation of the electron−hole Hamiltonian, where the excited states of the organic crystal are written as the superposition of local excitations and intermolecular CT states. On the basis of the TDFI-TI method, Fujita et al. parametrized the electron−hole Hamiltonian using the fragment-based electronic structure theory and applied it to the optical spectroscopy and quantum dynamics of a p-type organic semiconductor.43 As a similar approach, Li et al. have extended their ab initio exciton model to include CT states.44 In this study, we develop the excited-state method for calculating nonlocal excitations from the fragment molecular orbital (FMO) method.45,46 The FMO method is one of fragment-based molecular orbital theories suitable for treating large systems.8 In the FMO method, a system is divided into several fragments, and the energy of the whole system is approximated by the self-consistent field (SCF) calculations of fragment monomers and dimers embedded in the electrostatic potential from other fragments.47 The excited-state methods have been developed within the multilayer FMO (MLFMO) scheme,48 which calculates excited states of a target fragment after a ground-state FMO calculation. The local excitation of a target fragment monomer was treated using variants of CIS theory49,50 or TDDFT.51 The many-body correction to the local excitation energy was developed,9,52,53 utilizing the pairwise correction scheme proposed by Hirata and coworkers.54 However, the method mentioned above was restricted for treating local excitations, and the nonlocal excitations cannot be obtained. Very recently, Wen and Ma proposed the MLFMO-based excited-state theory for calculating collective excitations in conjugated polymers.55 Their method parametrizes the Frenkel exciton model from the excited-state calculations of fragment monomers and dimers, in combination with the diabatization scheme of Arago and Troisi.56 Here, we present the FMO-based excited-state theory for treating nonlocal excitations; the method presented here is an extension of the author’s previous effort43 that combines MLFMO48,49 and TDFI-TI.15,16 The method is based on the CIS theory and was designed to reproduce low-energy excited states for a whole system with small dimensions. To describe nonlocal excited states for the whole system, configuration state functions (CSFs) were constructed from localized MOs and the CIS amplitudes of intrafragment excitation via ground-state FMO and MLFMO−CIS calculations. The excited states of the whole system were written as the superposition of fragment CSFs for local excitations and interfragment CT excitations. Exploiting fragment CSFs allows for efficient truncation of the

proposed to calculate nonlocal excited states from localized MOs or fragment-based methods. Neugebauer and co-workers developed the subsystem TDDFT based on the frozen-density embedding.11,12 Mata and Stoll have developed an improved incremental approach for describing electronically excited state, with the inclusion of dominant natural transition orbitals.13 Fujimoto and co-workers developed the transition-density fragment interaction (TDFI) method to calculate electronic couplings for EET and treated the collective excited states within the quantum mechanics/molecular mechanics (QM/ MM) approach.14 Their method was later extended to include CT excitations,15,16 namely, the TDFI-transfer integral (TDFITI) method, wherein the interactions between local excitations and CT states were considered with the transfer integrals. Recently, Liu and Hybertshon developed the TDDFT based on localized nonorthogonal MOs.17 The equivalent wave function formulation, referred to as the absolutely localized molecular orbital-configuration interaction single (ALMO−CIS), was also proposed by Head-Gordon and co-workers.18,19 More recently, Nakai and Yoshikawa developed the linear-scaling TDDFT based on the divide-and-conquer method and dynamic polarizability calculation.20 Parker and Shiozaki proposed the sophisticated fragment-based theory,21 adopting an idea from the density matrix renormalization group.22 Their active space decomposition method21 describes a wave function of the whole system as a sum of products of fragment basis states, with monomer states being determined by the complete active space self-consistent field method. The aforementioned ab initio studies were restricted to the single point calculation of excited states. However, an optical transition to the excited state triggers various dynamic processes, such as EET and photoinduced electron transfer. Such dynamic processes have been studied with the help of a model Hamiltonian that originates from the exciton concept.23,24 The exciton model also adopts a similar idea to the fragment-based methods for constructing the excited-state wave functions. Frenkel first introduced the exciton concept, where the excited states of solids are described as a superposition of excitation waves.25 The idea was later extended to write an excited-state Hamiltonian in the following tight-binding form: ⟨I |H |J ⟩ = δIJEI + (1 − δIJ )VIJ

(1)

Here, |I⟩ is the localized excited state of the Ith molecule with its excitation energy EI, and VIJ is the excitonic couplingsthe Coulomb interaction between transition densities of local excitations. The use of this simplified model Hamiltonian allows to study exciton dynamics and linear and time-resolved optical spectroscopy.26−29 The excitonic couplings were determined by using a variety of approximations,29,30 while the ab initio parametrization of the Frenkel exciton model is also possible.31−35 Wannier36 and Mott37 introduced another model that describes the exciton state as a free electron and a free hole with their mutual Coulomb interaction. Knox developed the general exciton model that includes the Frenkel and Wannnier−Mott exciton models as the limiting case.23 The approach writes the excited-state Hamiltonian with the electron−hole basis: (e) ⟨eI hJ |H |eK hL⟩ = δJLHIK + δIK H (JLh) + HIJ ,KL

(2)

Here, |eIhJ⟩ denotes a state with electron and hole being localized in Ith and Jth molecules, respectively. The exciton Hamiltonian is composed of the one-body parts of electron B


Article

The Journal of Physical Chemistry A dimension of the excited-state Hamiltonian. The excited-state Hamiltonian represented with fragment CSFs is calculated, and the resulting eigenvalue problem is solved to determine the excited states of the whole system. The calculations of the singlet Hamiltonian matrix elements in the FMO framework is presented in detail: The one-electron part of the excited-state Hamiltonian was treated by the FMO-linear combination of molecular orbital (FMO-LCMO) method,57,58 while the twoelectron parts were efficiently treated within the two-body expansion of the FMO. We also introduce the approximation scheme to calculate the two-electron part, in accordance with the electrostatic approximation for separated fragment dimers.59 The proposed theory is incorporated into a developer version of the ABINIT-MP program.46 The accuracy and efficiency of the presented theory were tested for illustrative examples of molecular dimers and molecular aggregates. Dynamical correlation corrections will also be demonstrated.

For example, the matrix element between two monomer highest occupied molecular orbitals (HOMOs) yields electronic coupling for hole transfer, while that between two monomer lowest unoccupied molecular orbitals (LUMOs) provides electronic coupling for electron transfer.61−63 In terms of the model exciton Hamiltonian, the occupied−occupied and virtual−virtual blocks of the Fock matrix correspond to the one-body Hamiltonian for a hole and an electron in eq 2, respectively. We note that the nonorthogonality of interfragment MOs implies that the virtual MOs of a monomer are not strictly virtual, because they can overlap with the occupied MOs of other fragments. Within the ALMO−CIS method, Closser et al.18 have proposed to project occupied space out from virtual MOs, to avoid a possible mixing between ground and excited states. However, it is expected that the interfragment MOs overlap are small for weakly interacting molecular aggregates. Therefore, according to the original FMO-LCMO method, the monomer MOs were adopted as is to represent the Fock operator and the excited-state Hamiltonian in the following section. 2C. Excited-State Wavefuncton from Fragment CSFs. The singlet excited-state wave function for the whole system was constructed from fragment CSFs. To clarify the connection to the previous study,43 the following notation was adopted for singly substituted CSFs.

2. THEORY 2A. Notation. The following notation is used throughout the paper. We will use i,j as occupied MO indices, a,b as virtual MO indices, and p,q,r as the general MO indices. μ,ν,λ,σ refer to the atomic orbitals (AOs). I,J,K,L indicate fragment indices, and a fragment dimer is denoted by IJ. Hereafter, we will adopt “monomer” or “dimer” to refer to fragment monomer or fragment dimer. The MO of a monomer is represented by |ψ(I) p ⟩ (I) with its orbital energy as ϵ(I) and the MO coefficient as C . m,n i pi denote indices for intrafragment excited states, and bia(Im) indicates a singlet CIS amplitude for the mth excited state of Ith fragment. 2B. FMO-LCMO. We begin with the FMO-LCMO method, by which the one-electron parts of the excited-state Hamiltonian are calculated. In the FMO-LCMO method,57,58 the Fock operator for the whole system is expanded up to an Nbody correction, analogous to the FMO total energy. Within the two-body expansion of the FMO (FMO2), the Fock operator for the whole system, F, is approximated from monomer and dimer Fock operators: F=

|eIahJi⟩ =

I>J

|LE Im⟩ =

∑ bia(Im)|eIahIi⟩

(6)

ia

(3)

Intermolecular (CT) states, |eIahJi⟩ (I ≠ J) were also considered. The excited-state wave function for the whole system was described as the superposition of those LE and CT CSFs as follows:

where the summation signs should be taken in the tensor sense, so that each block, I or IJ, are added to the supermatrix F in the appropriate location. The monomer Fock operator, F(I), is (I) diagonal with respect to the monomer MOs, F(I) pq = ϵp δpq. The dimer Fock operator is represented in the monomer MOs and (IJ) (I) (IJ) IJ (I) is calculated according to57 F(IJ) pq = ∑rϵr ⟨ψp |ψr ⟩⟨ψr |ψq ⟩. Because the MOs between different fragments are not orthogonal, the overlap matrix between monomer MOs, Spq = (J) ⟨ψ(I) p |ψq ⟩, must be considered. Diagonalizing the FMO-LCMO Fock matrix as a generalized eigenvalue problem provides an approximate solution for the canonical orbitals of the whole system. The FMO-LCMO method was originally proposed to approximate canonical MOs of the whole system, while it can be used to compute electronic couplings for CT. In general, the electronic couplings are defined as the interaction between two diabatic states associated with EET or CT.60 In the FMOLCMO method, the electronic couplings for CT are calculated as the off-diagonal elements of the Fock matrix61−63 that is transformed by the Löwdin orthogonalization.64 F′ = S−1/2FS−1/2

(5)

† where aIa,α (aIa,β) is the creation (annihilation) operator of the α-spin (β-spin) electron in the ath orbital of the Ith fragment, and |G⟩ denotes the ground-state wave function for the whole system. The following fragment CSFs were considered to describe the excited-state wave function efficiently. The CSFs of intrafragment local excitations (LE) are described from the CIS amplitude (b(Im) ia ) and the intrafragment singly substituted CSFs:

∑ ⊕ FI + ∑ ⊕ (FIJ − FI ⊕ F J ) I

1 † † (aIa, αaJi, α + aIa, β aJi, β )|G⟩ 2

|Ψ⟩ =

CT |eIahJi⟩ ∑ c ImLE|LEIm⟩ + ∑ cIaJi

(7)

where |Ψ⟩ is the excited state wave function for the whole CT system, and cLE Im and cIaJi are the coefficients for fragment LE and CT CSFs, respectively. The coefficient of fragment CSFs were obtained by solving the eigenvalue problem Hc = cE. 2D. Excited-State Hamiltonian. This subsection derives and presents the matrix elements of the excited-state Hamiltonian represented in the fragment CSFs. The starting point of the derivation is the singlet CIS Hamiltonian with localized MOs.3,4 For convenience, the excited-state Hamiltonian is decomposed into one-electron and two-electron parts. ⟨eIahJi|H − E|eKbhLj⟩ = ⟨eIahJi|H1e|eKbhLj⟩ + ⟨eIahJi|H2e|eKbhLj⟩ (8)

We approximately rely on the conventional Slater−Condon rule to derive the Hamiltonian matrix elements. Here, the one-

(4) C


Article


process between LUMOs or HOMOs. Correction of the CIS amplitude can describe the interfragment electron or hole transfer following electron excitation from occupied to virtual MOs within a monomer.15 The CT-CT part has already been described in eq 9; it vanishes unless I = K or J = L. The matrix elements of the CT-CT blocks are responsible for delocalization of electron or hole wave function and thus formation of a conduction or valence band. The off-diagonal elements of the two-electron part contain four-center electron repulsion integrals (i(J)a(I)|j(L)b(K)) that can be involved in up to a fragment tetramer. In the present study, the four-center electron integrals that are involved in fragment trimers and tetramers were neglected in accordance with the FMO2. The following terms were considered for the offdiagonal elements of the two-electron part. The off-diagonal elements of the LE-LE block are given by

electron part is composed of the FMO-LCMO Fock operator in the virtual and occupied spaces. ′ − δIK δabFij′ ⟨eIahJi|H1e|eKbhLj⟩ = δJLδijFab

(9)

The two-electron part describes the electron−electron interactions represented in the monomer MOs. H2e = 2(i(J )a(I )|j(L) b(K )) − (i(J )j(L) |a(I )b(K ))

(10)

Here, the four-center electron repulsion integrals are (i(J )a(I )|j(L) b(K )) =

∫ dr1 ∫ dr2ψi(J)(r1)ψa(I)(r1) 1 ψ (L)(r2)ψb(K )(r2) |r1 − r2| j

(11)

Note that in eq 9 the transformed Fock matrix (F′) was adopted so that the interfragment MO overlaps can be effectively considered. Equations 9 and 10 were derived on the basis of the Slater−Condon rule, which assumes orthogonal Slater determinants. Nevertheless, the fragment CSFs defined by eq 5 are not orthogonal because of the nonzero overlaps between interfragment monomer MOs. The derivation of the matrix elements among nonorthogonal Slater determinants should follow the corresponding orbital transformation.65,66 For example, Morrison et al. have applied the corresponding orbital transformation to the ab initio exciton model.33 However, the computation of the Hamiltonian matrix elements from the corresponding orbital transformation requires the four-center AO integral of the whole system, which is computationally too demanding. Therefore, we decided to follow the conventional Slater−Condon by neglecting the overlaps among the CSFs but with the interfragment MO overlaps being effectively considered as the transformed Fock matrix in the one-electron part. The diagonal elements of the excited-state Hamiltonian in the fragment CSFs are excitation energies of the LE or CT states. The MLFMO−CIS calculations were performed for monomers after a ground-state calculation49 to obtain LE energies. When solving the CIS of the monomers, the orbital energies of the monomer (ϵ(I) p ) were replaced by corresponding ′ ). The diagonal elements of the FMO-LCMO Fock matrix (Fpp diagonal elements of the CT block were calculated according to

∑ ∑

′ − Fii′ + ⟨eIahJi|H |eIahJi⟩ = Faa

⟨LE Im|H2e|LEJn⟩ =

(Im) (Jn) T μν Tλσ [2(μν|λσ )

μν ∈ I λσ ∈ J

− (μλ|νσ )]

(15) 3

Here, the transition density matrix of the mth excited state of the Ith fragment is (Im) T μν =

∑ bia(Im)Cμ(Ii)Cν(aI)

(16)

ia

The first term in eq 15 describes Forster couplinga Coulomb interaction between transition densities. The second term describes the exchange coupling that involves simultaneous two-electron exchange. The two-electron part of the LE-CT block is given by ⟨LE Im|H2e|eIbhJj⟩ =

∑ ∑

(Im) (J ) (I ) T μν Cλj Cσb [2(μν|λσ )


− (μλ|νσ )] ⟨LE Im|H2e|eJbhIj⟩ =

∑ ∑

(17)

(Im) (I ) (J ) T μν Cλj Cσb [2(μν|λσ )


− (μλ|νσ )]

(18)

These terms can be considered as effects of the electron− electron interactions on the electron transfer or hole transfer. For the CT-CT blocks, we consider ⟨eIhJ|H2e|eIhJ⟩ and ⟨eJhI|H2e| eIhJ⟩ as follows:

Cλ(Ji )Cμ(Ji)Cν(aI)Cσ(Ia)


[2(μν|λσ ) − (μλ|νσ )]

∑ ∑

⟨eIahJi|H2e|eIbhJj⟩ =

(12)

∑ ∑

Cμ(Ji)Cν(aI)Cλ(Jj )Cσ(Ib)[2(μν|λσ )


The second and third terms on the right-hand side are electron−hole exchange and Coulomb interactions in ath and ith MOs, respectively. The off-diagonal elements of the one-electron part describe the electron or hole transfer processes. The elements of the LELE block vanish according to the Slater−Condon rule. As for the LE-CT block, the one-electron part was calculated by the CIS amplitude and the Fock matrix. ⟨LEIm|H1e|eIbhLj⟩ = −∑ bib(Im)Fij′ i

⟨LEIm|H1e|eJbhIj⟩ =

∑ b(jaIm)Fab′ a

− (μλ|νσ )]

(19)

This term contains short-range electron−hole exchange and long-range Coulomb interactions in the first and second terms, respectively. ⟨eIahJi|H2e|eJbhIj⟩ represents the simultaneous electron and hole exchange process: ⟨eIahJi|H2e|eJbhIj⟩ =

∑ ∑

Cμ(Ji)Cν(aI)Cλ(Ij )Cσ(Jb)[2(μν|λσ )


− (μλ|νσ )]

(13)

(20)

2E. Approximations for Separated Fragments. Here, the approximation methods for one-electron and two-electron parts in the excited-state Hamiltonian are presented. The FMO method adopts the electrostatic potential approximation for largely separated monomer pairs,59,67 wherein the dimer SCF

(14)

If the local excitation is composed of only the HOMO−LUMO transition, they reduce to the single electron or hole transfer D


Article


EESP using the one-electron integral and Mullikan atomic charges.59 2F. Parameters for Excited-State Wave Function. The number of fragment CSFs to describe the excited-state wave function is determined by the following parameters. One is the number of intrafragment excited states per fragment included as the LE states, nLE. nLE = M indicates that first M excited states of each fragment are included. The number of interfragment CT states was controlled by occupied and virtual MOs to define CT states, nocc and nvir. nocc = M indicates that HOMO (H) to H-M+1 orbitals are used as occupied MOs for interfragment CT states, while nvir = M indicates that LUMO (L) to L+M-1 orbitals used as virtual MOs. For an example of fragment dimer, (nocc, nvir) = (1,2) results in four interfragment CT states, |e1Lh2H⟩, |e1L+1h2H⟩, |e2Lh1H⟩, and |e2L+1h1H⟩. Different values can be used for different fragment (nLE(I),nocc(I),nvir(I)). The dimension of the wave function becomes Ndim = ∑InLE(I) + ∑I≠Jnocc(I)nvir(J), and the excited-state Hamiltonian contains Ndim2 matrix elements. (nLE,nocc,nvir) are treated as control parameters and must be specified in advance of an FMO excited-state calculation. For the case of (nLE,nocc,nvir) = (1,0,0), the Frenkel exciton model, which is of the same form as eq 1, can be derived. The two-body electron−hole Hamiltonian, in the form of eq 2, can be obtained by setting (nLE,nocc,nvir) = (1,1,1).43 FMO-LCMO Fock matrix constructions require other parameter sets: the numbers of occupied and virtual monomer MOs per fragment (nLCMO , nLCMO ). The dimensions of Fij and occ vir LCMO (I), respectively; they are Fab are ∑Inocc (I) and ∑InLCMO vir identical to the number of occupied and virtual MOs in the whole system if all occupied and virtual monomer MOs are used. However, the dimension of the FMO-LCMO Fock matrix may be reduced by considering a smaller number of monomer MOs. It is suggested from the earlier study that the canonical MOs of the whole system near the HOMO−LUMO level can be reproduced by a reduced FMO-LCMO Fock matrix.57 When calculating the excited-state Hamiltonian, the summa(IJ) (Im) (IJ) tions in eqs 13 and 14, −∑ib(Im) ib Fij and ∑abja Fab , run over the occupied and virtual MOs of Ith fragment included as the and FMO-LCMO Fock matrix, respectively. Therefore, nLCMO occ must be larger than nocc and nvir, respectively. We expect nLCMO vir that the summation in eq 13 or 14 can be well-approximated by considering a smaller number of occupied or virtual MOs, because a low-energy LE state is dominantly composed of a few CSFs (i.e., most of the CIS amplitude, b(Im) ia , are close to zero). ,nLCMO ) does not change the Note that the increase in (nLCMO occ vir dimension of the wave function but increases the computational times for constructing FMO-LCMO Fock matrix. 2G. Workflow. Finally, the computational procedure is presented briefly. Figure 1 shows the workflow of the present theory: 1. Monomer SCF calculations are performed according to the self-consistent charge (SCC) procedure. 2. Dimer SCF calculations are performed. 3. FMO-LCMO Fock matrix is constructed from the results of monomer and dimer SCF calculations. 4. MLFMO−CIS calculations are performed for each monomer, with the monomer orbital energies being replaced by the corresponding diagonal elements of the FMOLCMO Fock matrix. Resulting excitation energies are used as the diagonal elements of the LE states. 5. The excited-state Hamiltonian matrix elements are computed from the monomer MOs, the FMO-LCMO Fock matrix, and the CIS amplitudes of monomers. The two-electron terms that appear in Section 2D (e.g., eq 15) are computed as the Fock-like contraction.3,49

calculations for them are not performed. Adopting the electrostatic approximation for the separated pairs (dimer-es approximation) leads to neglecting the dimer contribution to the FMO-LCMO Fock operator, that is, F(IJ) = 0 for RIJ > Ldimer−es. Here, the interframement distance RIJ is defined as the shortest interatomic distance scaled by their van der Waals units, and Ldimer−es is the threshold value for the electrostatic approximation. Correspondingly, the one-electron part in the excited-state Hamiltonian (eqs 9, 13, and 14) was neglected. As well as the dimer-es approximation for the one-electron Hamiltonian, we introduced the additional approximation to neglect the short-range contributions in the two-electron part. The matrix elements of H2e (eqs 15−20) contain the longrange Coulomb and short-range exchange interactions. We neglect the short-range contributions and keep the long-range Coulomb interaction only, for the separated fragment monomers whose interfragment distance is larger than a threshold value. It follows from this approximation that eqs 15 and 19 are replaced by the following eqs 21 and 22, respectively. ⟨LE Im|H2e|LEJn⟩ = 2

∑ ∑

(Im) (Jn) T μν Tλσ (μν|λσ )

(21)


⟨eIahJi|H2e|eIbhJj⟩ = − ∑

∑

Cμ(Ji)Cν(aI)Cλ(Jj )Cσ(Ib)(μλ|νσ )


(22)

Equations 17, 18, and 20 are neglected if RIJ is larger than the threshold value, because they do not contain the long-range Coulomb interactions. Neglecting the short-range contributions were denoted as the dimer-2e approximation. The threshold value for the short-range two-electron part (Ldimer−2e) can be smaller than Ldimer−es, because the short-range contribution of the two-electron part is smaller than that of the one-electron part. On the one hand, it is known that the off-diagonal terms in the one-electron Hamiltonian, electron or hole transfer couplings, is exponentially dependent on the overlap between the interfagment MO, Sab or Sij. On the other hand, the exchange interaction (e.g., Dexter coupling in eq 15) is exponentially dependent on SabSij, which is smaller than the one-electron part by 2 or 3 orders of magnitude. The validity of this approximation will be discussed in the following section. The density−density Coulomb interactions in eqs 21 and 22 may be further approximated as the sums of interatomic Coulomb interactions. The atomic partial charges have been used to approximate the excitonic Coulomb couplings30 and electron−hole Coulomb interactions.68 Sets of atomic partial charges, which reproduce a transition density, HOMO density, or LUMO density, are determined by electrostatic charge procedures.69 Although we did not adopt this point-charge description in this study, it is useful for studying large systems43,70 and fluctuations of the excitonic couplings.71 We also combined the developed theory with the approximations of environmental electrostatic potentials.59,67 When solving an SCF equation of an n-mer in the FMO, the electrostatic potential acting on the n-mer from other fragments (environmental electrostatic potential (EESP)) may be approximated. We tested two approximation schemes to calculate the EESP.59 One was the Mulliken approximation (ESP-AOC), where the four-center AO integrals are approximated by three-center integrals.59. The other was the fractional point charge approximation (ESP-PTC) that approximates the E


Article


Figure 1. Workflow of the FMO excited−state calculation.

Alternatively, they may be approximated by neglecting the short-range contributions as discussed in Section 2E. Diagonalization of the excited-state Hamiltonian yields approximate solutions for the excited states of the whole system.

Figure 2. Geometries of the molecular dimers considered in this study: (a) ethylene dimer, (b) benzene dimer, and (c) ANS−BN dimer.

3. APPLICATIONS The present fragment-based theory was implemented in the developer version of ABINIT-MP program package.46 The FMO code in the ABINIT-MP program was first developed by Nakano et al.59 and later extended by Mochizuki et al. to implement the CIS49 and its modifications.6,50 In the following subsections, the accuracy of the present fragment-based theory was investigated by comparing to the conventional CIS results. All of the FMO and conventional CIS calculations were performed using the ABINIT-MP program, where the 6-31G* basis set was employed. Each molecule was assigned as one fragment in all FMO calculations. 3A. Molecular Dimer. Here, the accuracy of the present theory is evaluated for the systems of molecular dimers. To assess the accuracy, excitation energies for molecular dimers were calculated as a function of intermolecular separation (R) and were compared with the conventional CIS results. Because it is expected from earlier studies11,15,18 that the fragment-based theory works well for weakly interacting molecular systems, it is important to evaluate to what extent this expectation holds as the interaction becomes stronger at shorter intermolecular distances. As illustrative examples, molecular dimers of ethylene, benzene, and benzene derivatives were considered for π−π stacked systems. The molecular dimer of benzene derivatives is composed of anisole (ANS) and benzonitrile (BN), which have electron-donating and electron-withdrawing groups, respectively. The accuracy of the dimer-2e approximation was also tested. The other approximations in the FMO (e.g., dimer-es approximation) were not used in this subsection. Although different values can be used for nLCMO and nLCMO , the occ vir LCMO same number of MOs (n ) were adopted for simplicity to calculate the occupied−occupied and virtual−virtual blocks of the Fock matrix. The structures of the molecular dimers were prepared as follows: The structure of each molecule was first optimized at the MP2 level using the Gaussian09 software.72 Those of benzene derivatives were optimized, with carbon atoms constituting the benzene rings being fixed during the optimizations. The face-to-face dimer conformations, as seen in Figure 2, were prepared using the optimized monomer structures. For the ANS-BN dimer, R was defined as the

distance between the average positions of carbon atoms constituting the benzene rings. Because the present theory was developed to reproduce lowlying excited states, the mean absolute errors (MAEs) of several excited states were computed. The MAEs were calculated for six excited states in the ethylene dimer and benzene dimer and 10 excited states for the ANS-BN dimer. Table 1 shows the Table 1. MAEs of the FMO−CIS Excitation Energies for the Molecular Dimers at R = 3.0, 3.5, and 4.0 Å, with the Parameters for Calculating the Excited−State Hamiltonian System

nLCMO

(nLE,nocc,nvir)

MAE (meV)

Ethylene

6

(3,2,2)

Benzene

15

(4,3,3)

ANS−BN

21,b 19c

(4,5,5),b (7,5,5)c

R = 3.0 Å

3.5 Å

4.0 Å

160 108a 104 26a 34 30a

30 26a 45 23a 26 20a

7 6a 15 11a 24 11a

a

Ldimer−2e = 0. bParameters for anisole (ANS). cParameters for benzonitrile (BN).

MAEs of the molecular dimers calculated at R = 3.0, 3.5, and 4.0 Å, as well as the parameters for calculating the excited-state Hamiltonian. As expected, the results from the present theory are in quantitative agreement with the conventional CIS results at the R = 4.0 Å, and the MAEs are as small as 6−24 meV. However, the MAEs increase as the intermolecular interactions become greater with decreasing intermolecular distances. The MAEs are 23−45 meV at 3.5 Å, which are still marginal. Considerable MAEs lager than 100 meV were obtained at R = 3.0 Å for the ethylene dimer and benzene dimer. Table 1 also shows the MAEs calculated in combination with the dimer-2e approximation with Ldimer−2e = 0.0, which neglects all of the short-range contributions in H2e. The neglect of the short-range contributions tends to slightly improve the MAEs at R = 3.5 and 4.0 Å. For the ethylene dimer and benzene dimer at R = 3.0 Å, the results are rather improved by neglecting the short-range contributions, which is counterintuitive. It follows F


Article

The Journal of Physical Chemistry A that the inclusion of the short-range contributions in H2e lead to some imbalance in the excited-state Hamiltonian. The imbalance may be ascribed to the inconsistent treatments of one-electron and two-electron parts. We adopt the transformed Fock matrix, the monomer MOs of which are orthogonalized, to the one-electron parts. However, we use the nonorthogonalized monomer MOs to calculate the electron− electron interaction terms in the two-electron parts. The inconsistent treatment of the off-diagonal elements of oneelectron and two-electron parts may degrade the accuracy. Further discussion regarding the use of Löwdin-transformed Fock matrix for the excited-state Hamiltonian is presented in the Supporting Information, as well as the MAE results calculated with the nontransformed Fock matrix. To clarify the errors in more details, first to sixth excitation energies of the benzene dimer as a function of R are shown in Figure 3. Overall, the results of the conventional CIS are well-

Figure 4. Excitation energies of the first 10 excited states of the ANS− BN dimer. The results of conventional CIS and FMO−CIS are shown in solid lines and open circles, respectively. The first to tenth states are shown in red (S1), green (S2), blue (S3), purple (S4), cyan (S5), brown (S6), black (S7), orange (S8), navy (S9), and magenta (S10).

First, the nonzero intermolecular overlaps between occupied and virtual MOs imply that the ground-state component can be mixed with the excited-state wave function. Second, when being represented in the monomer MOs, the FMO-LCMO Fock matrix yields nonzero values for matrix elements in the occupied-virtual block. The nonzero occupied-virtual elements of the Fock matrix indicate that the Brillouin condition is not satisfied. Therefore, the ground-state wave function can be mixed with the excited-state wave function at the shorter intermolecular separations, which would decrease the accuracy. The results may be improved by projecting occupied MOs out from the virtual MOs, as proposed in the ALMO−CIS theory.18 Another possibility is to explicitly consider the interaction between the ground and excited states.73 Although the MAEs at R = 3.0 Å are considerable, we argue that the present theory would work with reasonable accuracy for practical applications. The face-to-face conformations of the molecular dimers are considerably short and unstable at R = 3.0 Å. For example, the minimum position of the face-to-face benzene dimer was calculated to be 3.8 Å.74 Considering that the MAEs at R = 4.0 Å are 7−24 meV, the present theory can work with satisfactory accuracy for realistic molecular assemblies in their equilibrium geometries. However, when dealing with conjugated polymers, the theory should be improved because of strong intramolecular interactions. The extension of the present theory to excited states of conjugated polymers is left for future studies. 3B. Molecular Array. This subsection briefly presents the dependence of the accuracy on total number of fragments. As described in Section 2D, the present theory approximates the two-electron terms by neglecting the four-center AO integrals of fragment trimers and tetramers. Thus, the systems of fragment dimer are not appropriate to validate the approximation. The accuracy may become worse with increasing the number of fragments due to the accumulation of the errors introduced by the AO integrals. Here, one-dimensional molecular arrays of ethylene molecules were employed to systematically change the number of molecules. The ethylene molecules were arranged into the faceto-face array with their spacing (d) of 3.0 or 4.0 Å. The geometry of each ethylene molecule was same as in Section 3A. MAEs relative to the conventional CIS were calculated for first N excited states, where N is total number of molecules (fragments). The excited-state Hamiltonians were calculated with (nLE,nocc,nvir) = (1,1,1) and nLCMO = 6. As well as in the Section 3A, the accuracy of the dimer-2e approximation was

Figure 3. Excitation energies of the first six excited states of the benzene dimer. The results of conventional CIS and FMO−CIS are shown in solid lines and open circles, respectively. The first to sixth states are shown in red (S1), green (S2), blue (S3), purple (S4), cyan (S5), and brown (S6).

reproduced by the FMO. For R > 4.0 Å, the results from the present theory are in quantitative agreement with the conventional CIS. However, the FMO fail to quantitatively reproduce the energy difference among excited states at shorter intermolecular distance. For example, the energy difference obtained by the FMO between S2 and S3 at R = 3.0 Å is 0.882 eV. The corresponding value by the conventional CIS is 0.732 eV. Because energy differences govern the photophysical processes such as EET and internal conversion, it is necessary to improve the present theory to treat strongly interacting systems. The first to tenth excitation energies of the ANS-BN dimer as a function of R are also shown in Figure 4. The MAEs of the ANS-BN dimer are smaller than those of the benzene dimer. In the ANS-BN dimer, the electron-withdrawing and -removal groups lead to more localized MOs than those of the benzene dimer, resulting in the smaller MAEs. The MAEs of the ANSBN dimer were also reduced, because the localized excited states appeared in the low-lying excited states. For example, the S7 and S8 states at R = 3.0 Å are localized states within the BN, and their excitation energies by the FMO−CIS are 6.527 and 6.919 eV, respectively. The corresponding values from the conventional CIS theory are 6.532 and 6.920 eV, respectively. From its construction for the excited-state wave function, the present theory can reproduce the intramolecular excitations. We mention that the present theory does not satisfy the Brillouin condition, which may lead to the mixing of ground and excited states at the shorter intermolecular separations. G


Article

The Journal of Physical Chemistry A tested, while the other approximations in the FMO were not used. Figure 5 presents the MAEs of the ethylene arrays of d = 3.0 or 4.0 Å. For d = 4.0 Å, the MAEs slightly increase with

Figure 6. (a) Chemical and (b) crystal structures of coumarin 4.

The excited-state Hamiltonian was calculated with (nLE,nocc,nvir) = (1,1,1) and nLCMO = 3. The S1 wave function of the isolated coumarine at the CIS level is dominantly composed of the 87% HOMO to LUMO, 5% HOMO to LUMO+1, and 3% HOMO−2 to LUMO+2 transitions. Therefore, three occupied and virtual MOs (nLCMO = 3) were included for constructing the FMO-LCMO Fock matrix. MAEs relative to the reference FMO calculation were calculated for the first 32 excited states. We confirmed that the inclusion of the second excited states of coumarine molecules by setting (nLE,nocc,nvir) = (2,1,1) have almost negligible effects on the S1state manifold. The total number of AO in this system was 6752. Table 2 shows the MAEs for the various calculation conditions. We found that the EESP approximations have

Figure 5. MAEs of the face-to-face ehylene arrays (N = 2, 4, 8, or 16) with spacing of 3.0 or 4.0 Å. Here, the MAEs were calculated for first N states.

increasing N, which is associated with the neglects of the fourcenter AO integrals. However, the MAEs of the system of N = 16 and d = 4.0 Å are still marginal (26−29 meV). It follows that the present theory can give satisfactory results for weakly interacting molecular aggregates. A similar but less systematic trend was observed for the systems of d = 3.0 Å. The MAE is largest for N = 2 and increases with increasing N for N = 4−16. With the dimer-2e approximation, the MAEs are almost independent of N for N = 4−16. Following the two-body formulation of the FMO-LCMO,57 the present theory approximates the Fock matrix of the total system from results of fragment monomers and dimers. However, for strongly interacting molecular aggregates, three-body terms in the Fock matrix58 can become significant. The results of d = 3.0 Å likely imply the importance of three-body effects both in the oneelectron (Fock matrix) and the two-electron (four-center AO integrals) parts of the excited-state Hamiltonian. 3C. Molecular Aggregate. In this subsection, we examine the accuracy and efficiency of approximations for calculating the excited-state Hamiltonian, for a realistic molecular aggregates. As well as the dimer-2e approximation, we also tested the dimer-es and EESP (ESP-AOC and ESP-PTC) approximations. Because these approximations are widely used in the applications of FMO to realistic systems, it is essential to see how the dimer-es and EESP approximations have influence on the excitation energies. To investigate the accuracy of dimer-es, dimer-2e, and EESP approximations, the results of approximate FMO will be compared with corresponding FMO calculations without the approximations. Although the ESP-AOC and ESPPTC approximations can be used together by adopting threshold values for switching them,59,67 we did not combine the EESP approximations in the present benchmark calculations. The molecular crystal of coumarin 4 was considered for the benchmark system. The optoelectronic properties of coumarin dyes are of broad interest, as they are utilized in a wide range of applications.75 The experimental crystal structure of coumarin 476 was taken from the Cambridge structure database. The 2 × 2 × 2 supercell including 32 molecules was prepared using the Mercury software 77 and was used for the benchmark calculations in this subsection. The chemical and crystal structures are shown in Figure 6.

Table 2. MAEs of the FMO−CIS Excitations Energies for the 32-Coumarin Cluster Obtained with Various Calculation Conditionsa No.

Ldimer−es

Ldimer−2e

EESP

MAE (meV)

1 2 3 4 5 6f

nonb nonb nonb 2.0 2.0 0.0

nonb nonb nonb 2.0 0.0 0.0

nonb AOCd PTCe nonb nonb nonb

N/Ac 7 × 10−4 0.4 6 8 36

a

Here, MAEs were calculated with respect to No. 1. bThe approximations were not used. cNot Assigned. dMulliken atomic population approximation. eMulliken point-charge approximation. f MLFMO−CIS with (nLE,nocc,nvir) = (1,0,0).

almost negligible effects on the excitation energies. Although the accuracy of monomer SCF calculations depend on the EESP approximations, the monomer orbital energies can be corrected by the FMO-LCMO method. The diagonal energy of (IJ) (I) the FMO-LCMO is Fpp = ϵ(I) p + ∑J≠I(Fpp − ϵp ); the second term indicates the shift of the monomer pth orbital energy by the interfragment exchange interactions. The second term also corrects the approximate treatment of EESP in the monomer SCF calculations, because the exact EESP can be recovered by (I) solving the dimer SCF and by correcting (F(IJ) pp − ϵp ) to the orbital energy. The errors from the EESP approximation result from the wave function of monomer MOs, which lead to errors in off-diagonal elements in the Fock matrix and electron− electron interactions in the H2e. The dimer-es and dimer-2e approximations have also negligible influences on the excitation energies. This can be expected from the previous subsection that the present theory H


Article


including the T1 correction, namely, PR-CIS(Ds).6 Here, the excitation energies of the LE states were corrected using CIS(D) or PR-CIS(Ds) theory but with other matrix elements of the excited-state Hamiltonian maintained at the CIS level. Here, the 1 × 1 × 1 supercell including four coumarin molecules was employed so that conventional CIS(D) or PRCIS(Ds) calculations are possible. The excited-state Hamiltonian was calculated with (nLE,nocc,nvir) = (1,1,1) and nLCMO = 3, where the ESP-AOC and separated dimer approximations with (Ldimer−es, Ldimer−2e) = (2.0, 0.0) were used. Table 3 shows

can provide satisfactory results at face-to-face molecular dimer R = 4.0 Å, where the wave function overlap is larger. For comparison, the results with (nLE,nocc,nvir) = (1,0,0) (i.e., without dimer SCF calculations) are also shown. The MAEs relative to the reference FMO calculation are still marginal (36 meV). We confirmed that the approximation schemes developed in the FMO method can be utilized in the present theory without reducing the accuracy. Next, we consider the computational times. The computational times reported here result from the calculations using 96 cores on the cluster comprising six Intel Xeon E5-2697 processors connected by InfiniBand. Figure 7 shows the

Table 3. First Four Excitation Energies of Coumarine Tetramer Using the FMO with CIS, CIS(D), or PRCIS(Ds)a FMO−CIS CIS FMO−CIS(D) CIS(D) FMO−PR-CIS(Ds) PR-CIS(Ds)

S1

S2

S3

S4

Time (min)

5.305 5.301 4.762 4.779 4.675 4.689

5.321 5.313 4.862 4.843 4.773 4.749

5.373 5.373 4.879 4.871 4.792 4.782

5.466 5.457 4.996 4.945 4.907 4.855

1.4 35.6 1.7 473.4 2.1 1175.1

a

Corresponding results from the conventional CIS, CIS(D), or PRCIS(Ds) are also shown. Results are given in electronvolts.

Figure 7. Computational times in minutes calculated with various calculation conditions. The index of calculation conditions corresponds to that of Table 2.

the first four excitation energies at the CIS, CIS(D), and PRCIS(Ds) levels for the coumarin tetramer. The computational times presented in Table 3 were obtained using 128 cores on the cluster comprising six Intel Xeon E5-2697 processors. Fortunately, the CIS(D) and PR-CIS(Ds) excitation energies can be well-reproduced by the present theory, and the errors relative to the conventional results are less than 50 meV. This agreement can be ascribed to the following reasons: the electronic couplings are well-approximated by the CIS theory,78 and the low-lying excited states of this system can be welldescribed as the superposition of the LE states. Although the results are very encouraging, the agreement observed here should not be generalized to other systems. For example, the energy difference between the LE and CT states is artificially overestimated by the correction to the LE energies, with CT energies being uncorrected. Consistent corrections to the LE and CT energies remain to be developed.

computational times decomposed into monomer SCC (Monomer), dimer SCF calculations and subsequent LCMO Fock constructions (Dimer+LCMO), and the MFMO−CIS and calculating the excited-state Hamiltonian (MLFMO). The calculation conditions correspond to those of Table 2. As shown in Figure 7, the computational times are dominated by the ground-state FMO calculations. Without the EESP approximation, the ground-state FMO calculations take considerable time, because the evaluation of the exact EESP requires four-center AO integrals, which formally scale as fourth order of the system size. The EESP approximations significantly reduced the timings of the ground-state calculations. Overall, the times for the MLFMO are comparable or less than those of ground-state calculations. For example, the times for the MLFMO part are 26.7 and 12.9 min for No. 4 and No. 5, ((Ldimer−es, Ldimer−2e) = (2.0, 2.0) and (2.0, 0.0)), respectively. Although we separately tested the separated dimer and EESP approximations, those approximations can be combined. If we use the ESP-AOC in combination with the dimer-es and dimer2e approximations with (Ldimer−es, Ldimer−2e) = (2.0, 0.0), the times are 5.9, 9.5, and 3.4 min for Monomer, Dimer+LCMO, and MLFMO, respectively. It follows that the present theory can compute a large number of excited states, with an additional computational time comparable to the ground-state calculation. Finally, we demonstrated the correlation correction to the excited-state Hamiltonian. It is well-known that the typical absolute error of CIS is as large as 0.5−1.0 eV.3 Thus, taking into account the dynamical correlations beyond the CIS level becomes essential for practical applications. We considered the perturbative treatment of dynamical correlation. Head-Gordon et al. developed a size-consistent second-order perturbative correction to the CIS, CIS(D).5 Mochizuki et al. extended the CIS(D) by introducing the partial renormalization scheme and

4. DISCUSSION Here, we discuss the possible future development of the present theory. Although we performed most of the benchmark calculations at the CIS level, it is essential to introduce dynamical correlation for practical applications. We tested the perturbative correction to the LE energies, but with CT energies being kept at the HF/CIS level. This correlation correction overestimates the energy difference between the LE and CT energies and, thus, underestimates the mixing between them. To avoid such an inconsistent treatment, both LE and CT excitation energies should be corrected at the same level of theory. The electron propagator theory79 may be useful for correcting dynamical correlation to the orbital energy for the CT energy, because the correlation energy can be computed at the same level as the perturbative CIS(D) level. Another possibility is to combine the present method with TDDFT within the Tamm−Dancoff approximation (TDA).80 The wellknown failure of TDDFT to underestimate CT excitation energies has been remedied by the long-range correction.81 Combined with the long-range corrected DFT and TDA, the I


Article


CIS result. We confirmed that the absolute errors of the excitation energies are less than 50 meV for the molecular assemblies in their equilibrium geometries, showing that the present theory can provide satisfactory results for practical applications. We have also evaluated the accuracy of various approximation methods developed in earlier FMO studies. As well as the dimer-es and EESP approximations, we have introduced the dimer-2e approximation, which neglects shortrange contributions of the two-electron part of the excited-state Hamiltonian. We have found only minor contributions from the exchange contributions and confirmed the efficiency of the dimer-2e approximation. In addition, the computational times can be significantly reduced by using the approximations, without affecting the accuracy. Although this article treats the singlet excited states, the calculation of the triplet excited states within the present theory is straightforward. Future work will be focused on improving the accuracy, which will extend the applicability of the present theory to treating larger and complex molecular systems. Specific features of our theory is the systematic derivation of an exciton model Hamiltonian and thus to bridge the gap between the ab initio electronic Hamiltonian and the model exciton Hamiltonian. Another advantage of our theory is the efficient computation of a large number of excited states. These features enable us to study optical properties and photophysical dynamics in extended molecular systems, like organic optoelectronic materials.95 Indeed, we have already applied this FMO-based theory to exciton dynamics in organic semiconductor thin films.43 Further applications to more complex systems will be considered in future works.

present theory will become a promising approach for computing excited states in large systems. We followed the original treatment of the FMO-LCMO and adopted monomer MOs to represent the excited-state Hamiltonian. Introducing a new set of occupied or virtual MOs may improve the accuracy of the present theory. As we discussed, the interfragment overlap between occupied and virtual MOs may result in the mixing of the ground and excited states. The accuracy of the present theory may be improved by defining new virtual MOs, with their occupied components being projected out. Another idea for defining new virtual MOs is a description of excitation energy or correlation energy with compact MO space. For example, virtual MOs included in the CIS calculation can be efficiently truncated by introducing the improved virtual orbitals82 or the modified virtual orbitals.83 Feller and Davidson introduced the K-orbital, which allows for efficient treatment of correlation energy.84 Flores-Moreno and Oritz adopted a similar idea within the electron propagator theory.85 Utilization of new MOs86 that allow for an efficient treatment of electron correlation will improve the accuracy of the present theory without increasing computational times. The remaining issue regarding the accuracy is to treat polarizable environments. In the FMO method, the electronic many-body polarization of the ground state is treated by the monomer SCC procedure.47 However, electronic excitations that have a large electronic density difference, such as a CT excitation, will result in additional electronic polarization of surrounding fragments. For explicit considerations of such electronic polarization, one may employ the GW approximation with the Bethe−Salpeter equation,87,88 where the dielectric function is computed for self-energy and screened Coulomb interactions. The additional electronic polarization of surrounding fragments may be treated by the response theory89 developed in the FMO. Another approach is to use nonrelaxed electron density in combination with polarizability calculations, in such a way that the polarization energy is calculated from the density difference of an excited fragment and induced dipole of surrounding fragments. We discussed the accuracy of the present theory in terms of the excitation energy. However, it is also essential to consider the accuracy of the excited-state wave functions. Although we expect that the present theory can also reproduce the excitedstate wave functions, a quantitative comparison should be made by analyzing the excited-state wave functions in detail. The excited-state wavefunctions have often been investigated by defining the natural transition orbitals.90−92 Classification of excited states inspired by the exciton concept has also been proposed.93,94 Implementation of the wave function analysis methods within the present theory will be useful for investigating excited states in large molecular systems.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpca.8b00446.

■

Further discussion on the use of the orthogonalized Fock matrix for excited-state Hamiltonian (PDF)

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: +81 564557260. Fax: +81 564534660. ORCID

Takatoshi Fujita: 0000-0003-1504-2249 Yuji Mochizuki: 0000-0002-7310-5183 Notes

The authors declare no competing financial interest.

5. CONCLUSION In this study, we have developed the FMO-based theory for calculating the nonlocal excitations of large systems. The theory writes the excited-state wave function of the whole system from fragment CSFs for intrafragment excitations and interfragment CT excitations, with their CSFs being determined by the FMO and MLFMO schemes. Computations of the excited-state Hamiltonian were presented in detail. The proposed theory was implemented into the developer version of the ABINIT-MP program.46 Using the face-to-face molecular dimers and the molecular aggregates as the benchmark systems, we evaluated the accuracy of the present theory relative to the conventional

■

ACKNOWLEDGMENTS T.F. thanks the financial support by Ministry of Education, Culture, Sports, Science and Technology (MEXT) as a Building of Consortia for the Development of Human Resources in Science and Technology. Y.M. acknowledges a partial support by MEXT as a social and scientific priority issue No. 6 (Accelerated Development of Innovative Clean Energy Systems) to be tackled by using post-K computer. Y.M. also acknowledges the support by MEXT as grant-in-aid (Kaken-hi) No. 16H04635. J


Article


■

(22) Chan, G. K.-L.; Sharma, S. The density matrix renormalization group in quantum chemistry. Annu. Rev. Phys. Chem. 2011, 62, 465− 481. (23) Knox, S., Ed. Theory of Excitons; Academic Press: New York, 1963. (24) Davidov, A. S. Theory of Molecular Excitons; Springer: New York, 1971. (25) Frenkel, J. On the transformation of light into heat in solids. I. Phys. Rev. 1931, 37, 17−44. (26) Abramavicius, D.; Palmieri, B.; Voronine, D. V.; Šanda, F.; Mukamel, S. Coherent multidimensional optical spectroscopy of excitons in molecular aggregates; quasiparticle versus supermolecule perspectives. Chem. Rev. 2009, 109, 2350−2408. (27) Spano, F. C. The spectral signatures of frenkel polarons in HAnd J-aggregates. Acc. Chem. Res. 2010, 43, 429−439. (28) Kühn, O.; Lochbrunner, S. In Quantum Efficiency in Complex Systems, Part II From Molecular Aggregates to Organic Solar Cells; Uli Wüerfel, M. T., Weber, T., Eds.; Semiconductors and Semimetals; Elsevier, 2011; Vol. 85, p 47. (29) Renger, T.; Muh, F. Understanding photosynthetic lightharvesting: a bottom up theoretical approach. Phys. Chem. Chem. Phys. 2013, 15, 3348−3371. (30) Madjet, M. E.; Abdurahman, A.; Renger, T. Intermolecular Coulomb couplings from ab initio electrostatic potentials: Application to optical transitions of strongly coupled pigments in photosynthetic antennae and reaction centers. J. Phys. Chem. B 2006, 110, 17268− 17281. (31) Ma, Y.; Liu, Y.; Ma, H. A new fragment-based approach for calculating electronic excitation energies of large systems. J. Chem. Phys. 2012, 136, 024113. (32) Sisto, A.; Glowacki, D. R.; Martinez, T. J. Ab initio nonadiabatic dynamics of multichromophore complexes: A scalable graphicalprocessing-unit-accelerated exciton framework. Acc. Chem. Res. 2014, 47, 2857−2866. (33) Morrison, A. F.; You, Z.-Q.; Herbert, J. M. Ab initio implementation of the Frenkel-Davydov exciton model: A naturally parallelizable approach to computing collective excitations in crystals and aggregates. J. Chem. Theory Comput. 2014, 10, 5366−5376. (34) Jurinovich, S.; Pescitelli, G.; Di Bari, L.; Mennucci, B. A TDDFT/MMPol/PCM model for the simulation of exciton-coupled circular dichroism spectra. Phys. Chem. Chem. Phys. 2014, 16, 16407− 16418. (35) Cupellini, L.; Jurinovich, S.; Campetella, M.; Caprasecca, S.; Guido, C. A.; Kelly, S. M.; Gardiner, A. T.; Cogdell, R.; Mennucci, B. An ab initio description of the excitonic properties of LH2 and their temperature dependence. J. Phys. Chem. B 2016, 120, 11348−11359. (36) Wannier, G. H. The structure of electronic excitation levels in insulating crystals. Phys. Rev. 1937, 52, 191−197. (37) Mott, N. F. Conduction in polar crystals. II. The conduction band and ultra-violet absorption of alkali-halide crystals. Trans. Faraday Soc. 1938, 34, 500−506. (38) Abramavicius, D.; Mukamel, S. Energy-transfer and chargeseparation pathways in the reaction center of photosystem II revealed by coherent two-dimensional optical spectroscopy. J. Chem. Phys. 2010, 133, 184501. (39) Yamagata, H.; Norton, J.; Hontz, E.; Olivier, Y.; Beljonne, D.; Brédas, J. L.; Silbey, R. J.; Spano, F. C. The nature of singlet excitons in oligoacene molecular crystals. J. Chem. Phys. 2011, 134, 204703. (40) Beljonne, D.; Yamagata, H.; Brédas, J. L.; Spano, F. C.; Olivier, Y. Charge-transfer excitations steer the Davydov splitting and mediate singlet exciton fission in pentacene. Phys. Rev. Lett. 2013, 110, 226402. (41) Tamura, H.; Burghardt, I. Potential barrier and excess energy for electron-hole separation from the charge-transfer exciton at donoracceptor heterojunctions of organic solar cells. J. Phys. Chem. C 2013, 117, 15020−15025. (42) D’Avino, G.; Muccioli, L.; Olivier, Y.; Beljonne, D. Charge separation and recombination at polymer-fullerene heterojunctions: Delocalization and hybridization effects. J. Phys. Chem. Lett. 2016, 7, 536−540.

REFERENCES

(1) Dreuw, A.; Head-Gordon, M. Single-reference ab initio methods for the calculation of excited states of large molecules. Chem. Rev. 2005, 105, 4009−4037. (2) Casida, M. E.; Huix-Rotllant, M. Progress in time-dependent density-functional theory. Annu. Rev. Phys. Chem. 2012, 63, 287−323. (3) Foresman, J. B.; Head-Gordon, M.; Pople, J. A.; Frisch, M. J. Toward a systematic molecular orbital theory for excited states. J. Phys. Chem. 1992, 96, 135−149. (4) Maurice, D.; Head-Gordon, M. Configuration interaction with single substitutions for excited states of open-shell molecules. Int. J. Quantum Chem. 1995, 56, 361−370. (5) Head-Gordon, M.; Rico, R. J.; Oumi, M.; Lee, T. J. A doubles correction to electronic excited states from configuration interaction in the space of single substitutions. Chem. Phys. Lett. 1994, 219, 21−29. (6) Mochizuki, Y.; Tanaka, K. Modification for spin-adapted version of configuration interaction singles with perturbative doubles. Chem. Phys. Lett. 2007, 443, 389−397. (7) Davidson, E. R. The iterative calculation of a few of the lowest eigenvalues and corresponding eigenvectors of large real-symmetric matrices. J. Comput. Phys. 1975, 17, 87−94. (8) Gordon, M. S.; Fedorov, D. G.; Pruitt, S. R.; Slipchenko, L. V. Fragmentation methods: A route to accurate calculations on large systems. Chem. Rev. 2012, 112, 632−672. (9) Chiba, M.; Fedorov, D. G.; Kitaura, K. Time-dependent density functional theory based upon the fragment molecular orbital method. J. Chem. Phys. 2007, 127, 104108. (10) Yoshikawa, T.; Kobayashi, M.; Fujii, A.; Nakai, H. Novel approach to excited-State calculations of large molecules based on divide-and-conquer method: Application to photoactive yellow protein. J. Phys. Chem. B 2013, 117, 5565−5573. (11) Neugebauer, J. Couplings between electronic transitions in a subsystem formulation of time-dependent density functional theory. J. Chem. Phys. 2007, 126, 134116. (12) König, C.; Schlüter, N.; Neugebauer, J. Direct determination of exciton couplings from subsystem time-dependent density-functional theory within the Tamm-Dancoff approximation. J. Chem. Phys. 2013, 138, 034104. (13) Mata, R. A.; Stoll, H. An incremental correlation approach to excited state energies based on natural transition/localized orbitals. J. Chem. Phys. 2011, 134, 034122. (14) Fujimoto, K. J. Transition-density-fragment interaction approach for exciton-coupled circular dichroism spectra. J. Chem. Phys. 2010, 133, 124101. (15) Fujimoto, K. J. Transition-density-fragment interaction combined with transfer integral approach for excitation-energy transfer via charge-transfer states. J. Chem. Phys. 2012, 137, 034101. (16) Fujimoto, K. J.; Kitamura, C. A theoretical study of crystallochromy: Spectral tuning of solid-state tetracenes. J. Chem. Phys. 2013, 139, 084511. (17) Liu, J.; Herbert, J. M. An efficient and accurate approximation to time-dependent density functional theory for systems of weakly coupled monomers. J. Chem. Phys. 2015, 143, 034106. (18) Closser, K. D.; Ge, Q.; Mao, Y.; Shao, Y.; Head-Gordon, M. Superposition of fragment excitations for excited states of large clusters with application to helium clusters. J. Chem. Theory Comput. 2015, 11, 5791−5803. (19) Ge, Q.; Mao, Y.; White, A. F.; Epifanovsky, E.; Closser, K. D.; Head-Gordon, M. Simulating the absorption spectra of helium clusters (N= 70, 150, 231, 300) using a charge transfer correction to superposition of fragment single excitations. J. Chem. Phys. 2017, 146, 044111. (20) Nakai, H.; Yoshikawa, T. Development of an excited-state calculation method for large systems using dynamical polarizability: A divide-and-conquer approach at the time-dependent density functional level. J. Chem. Phys. 2017, 146, 124123. (21) Parker, S. M.; Shiozaki, T. Quasi-diabatic states from active space decomposition. J. Chem. Theory Comput. 2014, 10, 3738−3744. K


Article

The Journal of Physical Chemistry A (43) Fujita, T.; Atahan-Evrenk, S.; Sawaya, N. P. D.; Aspuru-Guzik, A. Coherent dynamics of mixed Frenkel and charge-transfer excitons in dinaphtho[2,3-b:2′3′-f ]thieno[3,2-b]-thiophene thin films: The importance of hole delocalization. J. Phys. Chem. Lett. 2016, 7, 1374−1380. (44) Li, X.; Parrish, R. M.; Liu, F.; Kokkila Schumacher, S. I. L.; Martínez, T. J. An ab initio exciton model including charge-transfer excited states. J. Chem. Theory Comput. 2017, 13, 3493−3504. (45) Fedorov, D. G., Kitaura, K., Eds. The Fragment Molecular Orbital Method: Practical Applications to Large Molecular Systems; CRC Press: Boca Raton, FL, 2009. (46) Tanaka, S.; Mochizuki, Y.; Komeiji, Y.; Okiyama, Y.; Fukuzawa, K. Electron-correlated fragment-molecular-orbital calculations for biomolecular and nano systems. Phys. Chem. Chem. Phys. 2014, 16, 10310−10344. (47) Kitaura, K.; Ikeo, E.; Asada, T.; Nakano, T.; Uebayasi, M. Fragment molecular orbital method: An approximate computational method for large molecules. Chem. Phys. Lett. 1999, 313, 701−706. (48) Fedorov, D. G.; Ishida, T.; Kitaura, K. Multilayer formulation of the fragment molecular orbital method (FMO). J. Phys. Chem. A 2005, 109, 2638−2646. (49) Mochizuki, Y.; Koikegami, S.; Amari, S.; Segawa, K.; Kitaura, K.; Nakano, T. Configuration interaction singles method with multilayer fragment molecular orbital scheme. Chem. Phys. Lett. 2005, 406, 283− 288. (50) Mochizuki, Y.; Tanaka, K.; Yamashita, K.; Ishikawa, T.; Nakano, T.; Amari, S.; Segawa, K.; Murase, T.; Tokiwa, H.; Sakurai, M. Parallelized integral-direct CIS(D) calculations with multilayer fragment molecular orbital scheme. Theor. Chem. Acc. 2007, 117, 541− 553. (51) Chiba, M.; Fedorov, D. G.; Kitaura, K. Time-dependent density functional theory with the multilayer fragment molecular orbital method. Chem. Phys. Lett. 2007, 444, 346−350. (52) Chiba, M.; Koido, T. Electronic excitation energy calculation by the fragment molecular orbital method with three-body effects. J. Chem. Phys. 2010, 133, 044113. (53) Taguchi, N.; Mochizuki, Y.; Nakano, T. Fragment molecular orbital calculations for excitation energies of blue- and yellowfluorescent proteins. Chem. Phys. Lett. 2011, 504, 76−82. (54) Hirata, S.; Valiev, M.; Dupuis, M.; Xantheas, S. S.; Sugiki, S.; Sekino, H. Fast electron correlation methods for molecular clusters in the ground and excited states. Mol. Phys. 2005, 103, 2255−2265. (55) Wen, J.; Ma, H. A fragmentation-based approach for evaluating the intra-chain excitonic couplings in conjugated polymers. Chem. Phys. Lett. 2017, 679, 152−157. (56) Aragó, J.; Troisi, A. Dynamics of the excitonic coupling in organic crystals. Phys. Rev. Lett. 2015, 114, 026402. (57) Tsuneyuki, S.; Kobori, T.; Akagi, K.; Sodeyama, K.; Terakura, K.; Fukuyama, H. Molecular orbital calculation of biomolecules with fragment molecular orbitals. Chem. Phys. Lett. 2009, 476, 104−108. (58) Kobori, T.; Sodeyama, K.; Otsuka, T.; Tateyama, Y.; Tsuneyuki, S. Trimer effects in fragment molecular orbital-linear combination of molecular orbitals calculation of one-electron orbitals for biomolecules. J. Chem. Phys. 2013, 139, 094113. (59) Nakano, T.; Kaminuma, T.; Sato, T.; Akiyama, Y.; Uebayasi, M.; Kitaura, K.; et al. Fragment molecular orbital method: Use of approximate electrostatic potential. Chem. Phys. Lett. 2002, 351, 475− 480. (60) Hsu, C.-P. The electronic couplings in electron transfer and excitation energy transfer. Acc. Chem. Res. 2009, 42, 509−518. (61) Kitoh-Nishioka, H.; Ando, K. Charge-transfer matrix elements by FMO-LCMO approach: Hole transfer in DNA with parameter tuned range-separated DFT. Chem. Phys. Lett. 2015, 621, 96−101. (62) Kitoh-Nishioka, H.; Ando, K. FMO3-LCMO study of electron transfer coupling matrix element and pathway: Application to hole transfer between two tryptophans through cis- and trans-polyprolinelinker systems. J. Chem. Phys. 2016, 145, 114103.

(63) Fujita, T.; Haketa, Y.; Maeda, H.; Yamamoto, T. Relating stacking structures and charge transport in crystal polymorphs of the pyrrole-based π-conjugated molecule. Org. Electron. 2017, 49, 53−63. (64) Löwdin, P. On the non-orthogonality problem connected with the use of atomic wave functions in the theory of molecules and crystals. J. Chem. Phys. 1950, 18, 365−375. (65) Amos, A. T.; Hall, G. G. Single determinant wave functions. Proc. R. Soc. London, Ser. A 1961, 263, 483−493. (66) King, H. F.; Stanton, R. E.; Kim, H.; Wyatt, R. E.; Parr, R. G. Corresponding orbitals and the nonorthogonality problem in molecular quantum mechanics. J. Chem. Phys. 1967, 47, 1936−1941. (67) Fedorov, D. G.; Kitaura, K. The three-body fragment molecular orbital method for accurate calculations of large systems. Chem. Phys. Lett. 2006, 433, 182−187. (68) Yi, Y.; Coropceanu, V.; Brédas, J.-L. Exciton-dissociation and charge-recombination processes in pentacene/C60 solar cells: Theoretical insight into the impact of interface geometry. J. Am. Chem. Soc. 2009, 131, 15777−15783. (69) Singh, U. C.; Kollman, P. A. An approach to computing electrostatic charges for molecules. J. Comput. Chem. 1984, 5, 129− 145. (70) Fujita, T.; Huh, J.; Saikin, S. K.; Brookes, J. C.; Aspuru-Guzik, A. Theoretical characterization of excitation energy transfer in chlorosome light-harvesting antennae from green sulfur bacteria. Photosynth. Res. 2014, 120, 273−289. (71) Olbrich, C.; Strümpfer, J.; Schulten, K.; Kleinekathöfer, U. From atomistic modeling to excitation transfer and two-dimensional spectra of the FMO light-harvesting complex. J. Phys. Chem. B 2011, 115, 758−764. (72) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; et al. Gaussian 09, Revision E.01; Gaussian Inc.: Wallingford, CT, 2009. (73) Mochizuki, Y. A size-extensive modification of super-CI for orbital relaxation. Chem. Phys. Lett. 2005, 410, 165−171. (74) Tsuzuki, S.; Honda, K.; Uchimaru, T.; Mikami, M.; Tanabe, K. Origin of attraction and directionality of the π/π interaction: Model chemistry calculations of benzene dimer interaction. J. Am. Chem. Soc. 2002, 124, 104−112. (75) Liu, X.; Cole, J. M.; Waddell, P. G.; Lin, T. − C.; Radia, J.; Zeidler, A. Molecular origins of optoelectronic properties in coumarin dyes: Toward designer solar cell and laser applications. J. Phys. Chem. A 2012, 116, 727−737. (76) Yang, S. − P.; Han, L. − J.; Wang, D. − Q.; Xia, H. − T. Redetermination of 7−hydr−oxy−4−methyl−2H−1−benzopyran− 2−one. Acta Crystallogr., Sect. E: Struct. Rep. Online 2007, 63, o4643. (77) Macrae, C. F.; Edgington, P. R.; McCabe, P.; Pidcock, E.; Shields, G. P.; Taylor, R.; Towler, M.; van de Streek, J. Mercury: visualization and analysis of crystal structures. J. Appl. Crystallogr. 2006, 39, 453−457. (78) Yang, C. − H.; Hsu, C. − P. The dynamical correlation in spacer−mediated electron transfer couplings. J. Chem. Phys. 2006, 124, 244507. (79) Ortiz, J. V. Electron propagator theory: An approach to prediction and interpretation in quantum chemistry. WIREs Comput. Mol. Sci. 2013, 3, 123−142. (80) Hirata, S.; Head-Gordon, M. Time−dependent density functional theory within the Tamm−Dancoff approximation. Chem. Phys. Lett. 1999, 314, 291−299. (81) Tawada, Y.; Tsuneda, T.; Yanagisawa, S.; Yanai, T.; Hirao, K. A long−range−corrected time−dependent density functional theory. J. Chem. Phys. 2004, 120, 8425−8433. (82) Huzinaga, S.; Hirao, K. Improved virtual orbitals in the extended basis function space. J. Chem. Phys. 1977, 66, 2157−2160. (83) Bauschlicher, C. W. The construction of modified virtual orbitals (MVO’s) which are suited for configuration interaction calculations. J. Chem. Phys. 1980, 72, 880−885. L


Article

The Journal of Physical Chemistry A (84) Feller, D.; Davidson, E. R. An approximation to frozen natural orbitals through the use of the Hartree−Fock exchange potential. J. Chem. Phys. 1981, 74, 3977−3979. (85) Flores-Moreno, R.; Ortiz, J. V. Quasiparticle virtual orbitals in electron propagator calculations. J. Chem. Phys. 2008, 128, 164105. (86) Neogrády, P.; PitoňaK ́ , M.; Urban, M. Optimized virtual orbitals for correlated calculations: An alternative approach. Mol. Phys. 2005, 103, 2141−2157. (87) Hybertsen, M. S.; Louie, S. G. Electron correlation in semiconductors and insulators: Band gaps and quasiparticle energies. Phys. Rev. B: Condens. Matter Mater. Phys. 1986, 34, 5390−5413. (88) Rohlfing, M.; Louie, S. G. Electron−hole excitations and optical spectra from first principles. Phys. Rev. B: Condens. Matter Mater. Phys. 2000, 62, 4927−4944. (89) Nagata, T.; Brorsen, K.; Fedorov, D. G.; Kitaura, K.; Gordon, M. S. Fully analytic energy gradient in the fragment molecular orbital method. J. Chem. Phys. 2011, 134, 124115. (90) Martin, R. L. Natural transition orbitals. J. Chem. Phys. 2003, 118, 4775−4777. (91) Helmich, B.; Hàttig, C. Local pair natural orbitals for excited states. J. Chem. Phys. 2011, 135, 214106. (92) Baudin, P.; Kristensen, K. Correlated natural transition orbital framework for low−scaling excitation energy calculations (CorNFLEx). J. Chem. Phys. 2017, 146, 214114. (93) Mewes, S. A.; Plasser, F.; Dreuw, A. Communication: Exciton analysis in time−dependent density functional theory: How functionals shape excited−state characters. J. Chem. Phys. 2015, 143, 171101. (94) Hirose, D.; Noguchi, Y.; Sugino, O. Quantitative characterization of exciton from GW+Bethe−Salpeter calculation. J. Chem. Phys. 2017, 146, 044303. (95) Ostroverkhova, O. Organic optoelectronic materials: Mechanisms and applications. Chem. Rev. 2016, 116, 13279−13412.

M


Development of the Fragment Molecular Orbital Method for

Recommend Documents