Nonadiabatic Molecular Dynamics with Tight-Binding Fragment

Article pubs.acs.org/JCTC

Nonadiabatic Molecular Dynamics with Tight-Binding Fragment Molecular Orbitals Alexey V. Akimov* Department of Chemistry, University at Buffalo, The State University of New York, Buffalo, New York 14260-3000, United States ABSTRACT: This work presents a nonadiabatic molecular dynamics methodology that relies on the use of fragment molecular orbitals computed using tightbinding Hamiltonians. The approach aims to model charge and energy transfer in large systems via quantum-classical trajectory-based approaches. The technique relies on a chemically motivated fragmentation of the overall system into arbitrary fragments. Several types of fragment molecular orbitals (FMO) can be constructed and used in nonadiabatic simulations, comprising quasidiabatic, adiabatic, and Löwdin-transformed ones. The adiabatic FMOs are found to be most suitable for modeling nonadiabatic dynamics in complex molecular systems. The overall algorithm shows advantageous scaling properties, making it possible to model long time scale charge transfer processes in large systems with many hundreds of atoms. The approach is applied to study charge transfer in subphtalocyanine(SubPc)/C60 heterojunction. The computational results emphasize the importance of decoherence and details of interfacial structure for obtaining accurate charge transfer time scales in SubPc/C60 herejunctions.

1. INTRODUCTION Charge (CT) and energy (ET) transfer are the key processes that occur in many artificial and natural systems. Their outcomes determine the efficiency of photovoltaic and photocatalytic materials, light-emitting diodes (LED), and molecular electronics complexes. CT and ET play a key role in photosynthesis and light-induced biochemical processes in plants and animals. Direct simulation of CT and ET in molecular models of artificial materials and biological systems relies on the quantum-classical methods that describe coupled evolution of electronic excites states and nuclear dynamics. These methods commonly fall into a category of nonadiabatic molecular dynamics (NA-MD) techniques. NA-MD is a powerful tool that can provide comprehensive insights into the mechanisms of CT and ET processes, key to the operation of artificial and natural light-harvesting materials.1 This knowledge can help one to rationally design new materials or tune the properties of the existing ones. Many artificial and biological systems are very large for direct simulations that rely on a quantum description. As the solar energy research progresses, many novel nanostructures are synthesized and proposed for energy harvesting and conversion applications. The composites often rely on the nanostructures such as quantum dots,2−6 nanotubes,7−11 organic crystals, polymers and extended oligomers,12−19 or any combinations of the above. A direct modeling of CT and ET processes in these systems is problematic because of the rapidly growing costs of electronic structure calculations. Semiempirical20−22 or density functional theories (DFT)23 provide notable acceleration in comparison to ab initio calculations. In the past, several groups applied semiempirical24−29 or density functional tight-binding (DFTB) strategies to simulate NA-MD dynamics,30 extending © 2016 American Chemical Society

the modeling capabilities up to 1000 atoms or more. However, these methods rely on the use of all-system adiabatic states, which require expensive poorly scaling computations. Namely, these methods show a cubic scaling, arising from the eigenvalue problem and matrix multiplications for the full-rank matrices. The cubic scaling implies that the computational costs increase by almost an order of magnitude every time the system’s size is doubled. To surmount the scaling limitation, alternative formulations of the NA-MD are needed. The past three decades have witnessed an inspiring advancement in the field of linearly scaling methods for quantum chemical calculations.22,31−34 Fragment-based techniques such as divide-and-conquer (DC)35−40 or fragment molecular orbital (FMO)41−44 methods are perhaps the most popular nowadays. The methods rely heavily on the use of localized states characteristic to each subsystem (fragment) of the original system of interest. Specifically designed algorithms ensure self-consistency in determining the states of all fragments. The standard quantum chemical calculations are typically performed in the adiabatic representation. The adiabatic states are intrinsic to the whole system and are often delocalized over a notable part of the system. Thus, the computations are expensive and require full-dimensional Hamiltonian diagonalization. The fragmentation approaches break the dimensionality bottleneck by partitioning the overall system into smaller fragments, each of which can be treated much faster. Such calculations produce localized states that can be considered quasidiabatic. In the limit of vanishing interaction between the Received: September 27, 2016 Published: November 17, 2016 5719

DOI: 10.1021/acs.jctc.6b00955 J. Chem. Theory Comput. 2016, 12, 5719−5736

Article

Journal of Chemical Theory and Computation

heterojunction are essential factors that can not be ignored if accurate estimates of the CT rates are desired.

fragments, these states become diabatic (noninteracting) with respect to the states of other fragments, yet they remain adiabatic with respect to other states of the same fragment. A number of techniques can be used to generate such quasidiabatic states: block-localized wave function theory (BLWFT),45−47 block-localized DFT,48 molecular-orbital valence bond,49−51 fragment orbital DFT,52,52,53 constrained DFT,54−59 superposition of fragment states,60 dimer splitting,61,62 or fragment orbital63−65 methods, to name a few. The approaches produce charge-separated states in a self-consistent manner. The self-consistency is essential for static electronic structure calculations but is unnecessary for dynamical simulations. Instead, one can utilize an appropriately chosen basis set of physically meaningful states. Such an approach is computationally advantageous, because the electronic structure calculations can be performed on all subsystems independently, without the need for self-consistency loops that incur additional costs. In the past, localized quasidiabatic states have been utilized with some success in phenomenological models of charge transport.45−53,60,66−73 Nonetheless, such approaches lack the atomistic description of the internal dynamics. For instance, the block-localized wave function (BLWF) method45−47 had been utilized to compute the parameters of the phenomenological Marcus-type CT theories.67,68 Several other theories follow the same phenomenological approach: the diabatic coupling elements are computed, either from the charge-localized states48,53,72 or even from the delocalized adiabatic wave function.74−76 The couplings are then utilized in a simple twolevel Marcus-type expression to analyze the CT dynamics. Because the nuclear dynamics is largely neglected, the effects of molecular orientation are oversimplified,72 although the couplings may be very sensitive to molecular orientation and packing.52,77 Most importantly, the presence of multiple electronic states close in energy to those of the donor and acceptor may need to be included, because the transitions via these states may facilitate the overall CT dynamics. To account for these effects, an extension beyond the standard 2-level Marcus-type model is required. In this work, I report a fully atomistic NA-MD approach with the quasi-diabatic fragment-localized electronic states used as the basis. The fragment-localized states are constructed in independent electronic structure calculations on isolated fragments. No self-consistent adjustment of the computed states or diabatization74−76,78−80 of the adiabatic states is required, making the approach computationally efficient and applicable to extended systems. The computational scheme rests on the classical path approximation for electron−nuclear dynamics augmented with the molecular mechanics description of interatomic interactions that govern the evolution of nuclei. The electronic structure of the fragments is described using the semiempirical extended Hückel theory (EHT), which is an inexpensive approximation of a tight-binding type. Despite its simplicity, it has a great capacity of providing qualitative and semiquantitative insights into electronic structure and dynamics of diverse molecular and solid-state systems. The methodology is augmented with the trajectory surface hopping (TSH) simulation and optional decoherence correction. The methodology is validated by its application to modeling the CT dynamics in subphthalocyanine (SubPc)/fullerene (C60) heterojunction. The analysis of the computational results suggests that decoherence and surface structure of a

2. THEORY AND METHODS The starting point of the present approach is the system’s fragmentation into smaller subsystems. An exemplary partitioning scheme for a prototypical A-L-B system is presented in Figure 1. In the scheme, L can be associated with a linker

Figure 1. Possible fragmentation schemes for a prototypical A-L-B system. Scheme 1: relying on nonoverlapping fragments. Scheme 2: involving overlapping fragments.

fragment, whereas A and B can be donor/acceptor groups. Of course, more general assignments of the fragments can be envisioned. The fragments can be either nonoverlapping (not having common atoms) or overlapping. In both cases, the fragments can be defined quite arbitrarily, including the cases that require cutting covalent bonds. Multiple distinct partitioning schemes can be considered. The molecular orbitals of the resulting fragment sets can be concatenated to yield a complete basis sufficient for modeling quantum dynamical phenomena in complex systems. The present work reports a NA-MD method formulated in a basis of fragmental states computed using tight-binding Hamiltonian. The mathematical and computational details needed to implement and utilize the approach are presented. The formulation is sufficiently general and is applicable to the case of arbitrary overlapping fragments. The present work provides the proof of the principle and establishes the mathematical and computational basis of the method. Thus, in the application part, only the simplest case of nonoverlapping, noncovalently bound fragments is considered. In the present approach, electronic states of all fragments are determined in a nonself-consistent fashionvia computing orbitals of each individual fragment. The resulting orbitals of each fragment, {ψi(A)}, form an orthonormal set within (A) themselves, ⟨ψ(A) i |ψj ⟩ = δij, but they are nonorthogonal to (B) orbitals of other fragments, ⟨ψ(A) i |ψj ⟩ = SIJ. Here, the capital letters index the orbitals on different fragments in the total pool of orbitals: I = (i, A) and J = (j, B). In the future discussion, the lower-case indices will be utilized, keeping in mind that they refer to the presently capitalized letters. 2.1. Quantum-Classical Dynamics in the Nonorthogonal Basis. Coherent dynamics of electronic states is a starting point for many quantum-classical NA-MD formulations. Because the constructed fragment orbitals used as the basis for coherent dynamics are nonorthogonal, it is important to revise the relevant formulations. The approach, being a quantum-classical trajectory-based method, treats nuclei classi5720


Article


within an orthonormal basis, |ψ̃ ⟩. Also, the definition and analysis of state populations are much more transparent in orthogonal bases. Thus, the two most natural choices for selecting such bases are considered: (a) Löwdin-transformed fragment molecular orbitals (L-FMO), |ψ̃ L‑FMO⟩ and (b) adiabatic fragment molecular orbitals (A-FMO), |ψ̃ A‑FMO⟩, which diagonalize the electronic Hamiltonian in the reduced subspace of (quasi)diabatic fragment molecular orbitals (DFMO), |ψD‑FMO⟩ = |ψ⟩. Löwdin transformation is known to minimize the difference between the original and transformed orbitals in the mean square sense.82 Thus, the formulation of the surface hopping procedures in the orthogonal states, |ψ̃ L⟩, would resemble that in a chosen set of quasidiabatic states. On the contrary, the choice of the A-FMOs (yet only within the subset of all available orbitals) is the closest to the global adiabatic states of the overall system. When all generated quasidiabatic orbitals would be used for its construction, it is expected to produce the results converging to those from the full adiabatic MO simulation. The two types of orbitals are defined according to

cally whereas electrons are described quantum mechanically. The electrons are described by an overall wave function, |Ψ⟩, represented as a coherent superposition of system’s eigenstates: |Ψ⟩ =

∑ ci(t )|ψi(r ;R(t ))⟩ = |ψ ⟩c(t ) i

(1)

Here, c(t) = (c0(t) c1(t) ... cN−1(t))T are the explicitly timedependent amplitudes of the corresponding electronic states, |ψ(r;R(t))⟩ = |ψ⟩ = (|ψ0⟩ |ψ1⟩ ... |ψN−1⟩). The latter are parametrized by nuclear configurations, R(t), at which the electronic states are computed. These configurations evolve according to a predefined procedure, leading to the corresponding classical trajectories. Most often, the trajectories are found by integrating classical Newton’s equations of motion using a given potential energy surface (PES). This approach is known as a classical path approximation (CPA).81 The choice of the PES can vary. In the Ehrenfest, or the mean field (MF), method the governing PES is found as a superposition of all considered electronic PES, weighted by the corresponding populations and coherences. In the trajectory surface hopping (TSH) approaches, a single (“current”) PES is used, the one corresponding to an instantaneous electronic state of the system. A simpler version of the TSH, which is also called CPA assumes that the nuclei are always governed by the ground-state PES. The evolution of a quantum-classical system is described by the time-dependent Schrodinger equation (TD-SE):

|ψ̃ L‐FMO⟩ = |ψ ⟩S −1/2

and |ψ̃ A‐FMO⟩ = |ψ ⟩U HD‐FMOU = SH A‐FMOU

dc = H vibc dt

(7)

thus leading to the relationships between the time-evolved amplitudes:

(3a)

(3b)

c L̃ ‐FMO(t ) = S1/2c(t )

(8a)

cÃ ‐FMO(t ) = U −1c(t ) = U +Sc(t )

(8b)

2.2. Vibronic Hamiltonian. 2.2.1. Electronic Couplings. The vibronic Hamiltonian, eq 4, contains two main components: the electronic coupling of the quasi-diabatic ̂ (B) fragmental orbitals, ⟨ψ(A) i |H|ψj ⟩, and the corresponding NAC d terms, ⟨i| dt |j⟩. The computation of electronic couplings has been discussed by many authors within different levels of theory. For instance, van Voorhis developed the constrained DFT (c-DFT)55−57 approaches to compute diabatic (electronic) couplings between donor and acceptor states defined by particular spatial localization of the charge densities. Very recently, Ren et al.83 have adapted the multistate density functional theory (MSDFT)48,69 to compute the electronic couplings. Many authors considered computing electronic coupling via various diabatization schemes. For instance, Voityuk74−76 and Truhlar78−80 developed a number of methods that approximate diabatic couplings on the basis of adiabatic properties of the system. However, such approaches would compromise the overall fragmentation idea for reduced scaling:

where Hvib is the vibronic Hamiltonian, whose matrix elements are defined as

H vib, ij = Hij − iℏdij

(6)

|Ψ⟩ = |ψ ⟩c(t ) = |ψ̃ L‐FMO⟩c L̃ ‐FMO(t ) = |ψ̃ A‐FMO⟩cÃ ‐FMO(t )

or, in a matrix-vector form iℏS

H A‐FMO = diag(ε A‐FMO)

The dependence of quasidiabatic electronic states on nuclear coordinates is only parametric. The functional time dependence of the overall wave function, |Ψ⟩, on time is encoded into the evolution of electronic amplitudes of the corresponding basis states. The overall wave function, eq 1, remains the same in all representations:

∑ (Hel,ij − iℏdij)cj j

(5b)

where

d|Ψ⟩ iℏ = Ĥ |Ψ⟩ (2) dt Here, Ĥ = T̂ + Ĥ el is the quantum-mechanical Hamiltonian, 1 T̂ = ∑α − 2M ∇α 2 is the nuclear kinetic energy operator, and α Ĥ el is the electronic Hamiltonian operator. Within the general CPA framework, eq 1, the TD-SE, eq 2, reduces to the TD-SE describing coherent dynamics of electronic degrees of freedom, ci(t): ∂c (t ) iℏ ∑ Sij i = dt j

(5a)

(4)

d ⟨ψi| dt |ψj⟩

Here, dij = is the nonadiabatic coupling (NAC) scalar, Hij = ⟨ψi|Ĥ el|ψj⟩ is the matrix element of the electronic Hamiltonian, and Sij = ⟨ψi|ψj⟩ ⇔ S = ⟨ψ|ψ⟩ is the overlap of two electronic states (or an overlap matrix). Most commonly, the electronic states, |ψ⟩, form an orthonormal basis (⟨ψi|ψj⟩ = δi,j, where δi,j is a Kronecker delta symbol), so eq 3 can be simplified to drop the overlap matrix. However, in the present formulation, the electronic states that belong to different fragments are obtained in separate calculations and therefore are not orthogonal. Thus, the overlap matrix should be retained in eq 3. Nonetheless, although the TD-SE, eq 3b, can be integrated in a nonorthogonal basis of independent fragment states, |ψ⟩, the surface hopping procedures are formulated 5721


Article


considered a simplified version of a more elaborate selfconsistent procedure.45−47 Unlike the latter, the present formulation does not involve self-consistent determination of the fragment orbitals. In principle, the consistency may be important because it can stabilize charge-separated states, especially in polar environments. Due to the tight-binding and semiempirical nature of the EHT approach, the energy of electronic states does not depend on the charge density distribution (core Hamiltonian approximation) but can be calibrated to reproduce experimentally known energy levels. In a sense, the present approach extends the tight-binding AOEHT to the FMO-EHT level. As with many semiempirical theories, essential correlation and electrostatic effects can be implicitly included via a parametrization tuned to a specific system. Further extensions to incorporate charge-transfer effects within the tight-binding formulation can be envisioned along the lines of the tight-binding configuration interaction (TB-CI) method developed in the Truhlar group.85 The approximations above have been discussed keeping only the nonorthogonal basis states in mind, |ψ⟩. In the basis of the Löwdin-orthogonalized states, |ψ̃ L⟩, the corresponding electronic Hamiltonian is obtained by the similarity transformation:

the computation of the adiabatic wave functions of the overall system is necessary to determine the adiabatic properties such as charges or dipole moments, yet it represents the timelimiting step we try to avoid using the fragmentation ansatz. It is important to recall that once the diabatic states are constructed, the computation of electronic coupling is straightforward if the action of the system’s Hamiltonian on the corresponding MOs or AOs is known. Such action can be derived explicitly for the ab initio or DFT Hamiltonians, involving 1- and 2-electron integrals and leading to the most rigorous formulation. For semiempirical methods, the Hamiltonian operator is commonly defined in the AO space via the corresponding matrix elements, Hab = ⟨χa|Ĥ |χb⟩. Here, the AOs |χa⟩ and |χb⟩ can belong to the same or different fragments. Specifically, in the extended Hückel theory (EHT) we utilize in this work, the AO matrix elements are given in terms of valence-state AO energies, Ii, and the AO overlap integrals, SAO ab = ⟨χa|χb⟩. With each fragment MO expanded in the basis of the same fragment AOs: |ψi(A)⟩ =

∑ TA,ai|χa ⟩

(9)

a∈A

S −1/2HD‐FMOS −1/2 = HL‐FMO

the electronic couplings between different fragments are computed as ⟨ψi(A)|Ĥ |ψj(B)⟩

=

∑ a∈A b∈B

=

(11)

If the adiabatic fragment |ψ̃ A‑FMO⟩ basis is utilized, the electronic Hamiltonian becomes diagonal, and the states become coupled only via nonadiabatic couplings (adjusted accordingly). If the A-FMO basis is utilized, the computational scheme follows closely the standard A-MO scheme. 2.2.2. Nonadiabatic Couplings. Nonadiabatic coupling is the second ingredient needed for constructing vibronic Hamiltonian, eq 4, which governs coherent dynamics of the selected basis states. In the present work, the NAC computation is based on the numerical differentiation, following the standard technique of Hammes-Schiffer and Tully.86 The approach assumes the orthogonality of the basis states. Thus, it is directly applicable to fragment orbitals of LFMO or A-FMO types, |ψ̃ L‑FMO⟩ or |ψ̃ A‑FMO⟩. However, if the nonorthogonal states of D-FMO type, |ψ⟩, are utilized, the formulation must be generalized. Specifically,

(B) ̂ Ta(A) i ⟨χa |H |χb ⟩T bj

D‐MO AO = TATHAB TB ∑ Ta(iA)HabTb(B)j ⇔ HAB a∈A b∈B

(10)

Here, TA is the MO-LCAO coefficients matrix for the fragment A found via an independent diagonalization of the corresponding fragment Hamiltonian in the AO basis: HAO AA → TA. An alternative approach to compute electronic couplings between fragment-localized states would involve the use of adiabatic molecular orbitals of the overall system, as was defined in several recent works.84 As it has been noted above, computing the adiabatic MOs of the entire system would require the diagonalization of the overall Hamiltonian, compromising the fragmentation idea for large-scale systems. Thus, in the present work we adapt the approach based on the transformation, eq 10, as it provides a compromise between simplicity and computational efficiency. One of the key elements of the present computational approach, is the reduction of the AFMO basis size: although all localized MOs are determined for each fragment, only a subset of orbitals from each fragment is utilized as a basis for NA-MD simulations. Only the orbitals within a certain energy window, suitable for the problem of interest, are selected. Thus, the overall dimensionality of the AFMO basis is drastically reduced in comparison to the total number of the adiabatic MOs (A-MO) or D-FMOs one could have otherwise. This reduction allows minimizing the costs of the matrix transformations in eq 10. At the same time, the transformation in eq 10 involves the construction of the overall system Hamiltonian in the AO basis. This operation becomes a computational bottleneck, determining the efficiency of the overall scheme for computing electronic couplings. With the tight-bonding EHT approach adopted in the present work, a quadratic scaling with the number of orbitals is expected. In eq 10, the D-FMO Hamiltonian of the overall system is constructed by the block-diagonalization procedure. It can be

⎛ dt ⎞ d ⎛ dt ⎞ Dij⎜t + ⎟ ≡ ⟨i| |j⟩⎜t + ⎟ ⎝ ⎠ ⎝ 2 dt 2⎠ ≈

ψi(t + dt ) + ψi(t ) ψj(t + dt ) − ψj(t ) dt

2

⎛ dt ⎞ 1 dS ⎜⎛ dt ⎞ = dij⎜t + ⎟ + t+ ⎟ ⎝ ⎠ ⎝ 2 2 dt 2⎠

(12)

where ⟨i(t )|j(t + dt )⟩ − ⟨i(t + dt )|j(t )⟩ ⎛ dt ⎞ dij⎜t + ⎟ ≡ ⎝ ⎠ 2 2dt

(13)

⟨i(t + dt )|j(t + dt )⟩ − ⟨i(t )|j(t )⟩ dS ⎛⎜ dt ⎞ t+ ⎟≡ dt ⎝ 2⎠ dt

(14)

The terms, dij, are the same as in the original HammesSchiffer−Tully formulation. One can observe that the 1 dS corresponding matrix d is antisymmetric, dij = −dji. The 2 dt term, however, is a new feature that originates from the nonorthonormality of the basis states. Essentially, it describes how fast the orbital overlap changes in time (simple finite 5722


Article

Journal of Chemical Theory and Computation difference formula). Because in the orthonormal basis the 1 dS overlap matrix, S, does not change in time, the 2 dt term vanishes exactly. As a consequence, the standard HammesSchiffer−Tully formula is recovered. In the case of nonorthogonal states, the term remains. Similar to the matrix S, the d matrix dS is symmetric. Hence, the sum −iℏ⟨i| dt |j⟩ is neither dt symmetric nor antisymmetric, making the vibronic Hamiltonian, eq 5, non-Hermitian. Working with non-Hermitian Hamiltonians is disadvantageous, because the solution of the corresponding TD-SE is not bound to preserve the wave function norm. This may cause numerical instabilities and loss of accuracy in long-time simulations. In the following section, I will show how the TD-SE with such a non-Hermitian vibronic Hamiltonian can be integrated stably and accurately. 2.3. Integration of the TD-SE. The integration of the TDSE, eq 3b, is straightforward if electronic (fragment) basis states are orthogonal. In this case, the overlap matrix on the left-hand side reduces to the identity matrix and the vibronic Hamiltonian on the right-hand side is Hermitian. However, if the nonorthogonal basis is used, one needs a reliable, robust, and convenient way to solve the equation. The complications come due to (1) the presence of the overlap matrix and (2) the non-Hermiticity of the vibronic Hamiltonian (as shown in the previous section). To solve the TD-SE, eq 3b, we first perform the vector transformations: C = S1/2c

(15a)

c = S −1/2C

(15b)

iℏ

⎛ i·dt −1/2 ⎞ C(t + dt ) = exp⎜ − S H̃ vibS −1/2⎟C(t ) ⎝ ℏ ⎠

(

The exponent, exp −

(16)

⎛ 1 dc dS dC ⎞ = iℏ⎜ − S −1/2 C + S1/2 ⎟ = H vibS −1/2C ⎝ 2 dt dt dt ⎠ (17)

After multiplying both sides of eq 17 by S−1/2 from the left, one obtains ⎛ 1 dS dC ⎞⎟ iℏ⎜ − S −1 C + = S −1/2H vibS −1/2C ⎝ 2 dt dt ⎠

Rearranging the terms, one solves for iℏ

(18)

(19)

P D‐FMO(t ) = |⟨ψ |ψ̃ L‐FMO⟩c L̃ ‐FMO(t )|2 = |⟨ψ |ψ ⟩S −1/2|2 |c L̃ ‐FMO(t )|2

(20)

= |S1/2|2 |c L̃ ‐FMO(t )|2

where H̃ vib, ij = Hij − iℏdij

(24)

Here, |ψO⟩ is the basis of observable states (in the sense that their population dynamics is being tracked). For instance, if the coefficients c̃L‑FMO(t) = C are evolved (one of the most straightforward practical approaches), and we are interested in the population of the original nonorthogonal basis states, |ψ⟩, the result can be computed by

At this point, it is important to recall that the vibronic Hamiltonian, Hvib, is not Hermitian but can be represented as a iℏ dS sum of Hermitian (H̃ vib ) and anti-Hermitian (− 2 dt ) components: iℏ dS H vib = H̃ vib − 2 dt

), can be computed by

P(t ) = |⟨ψ O|Ψ⟩|2

dC : dt

dC iℏ dS = S −1/2H vibS −1/2C + S −1 C dt 2 dt

i·dt −1/2 S H̃ vibS −1/2 ℏ

(23)

first diagonalizing the argument, forming a diagonal matrix with the diagonal elements set to the exponents of the eigenvalues, and transforming the resulting matrix back to the original representation. The detailed description of the procedure can be found elsewhere.87 One can see that both H̃ vib and S−1/2H̃ vibS−1/2 are Hermitian matrices. Therefore, the integration of eq 22 preserves the norm of the wave function, as given by C+C = c+Sc. The proof is straightforward, especially keeping eq 23 in mind. Comparing the transformations in eqs 15 and 8a as well as the expressions for the effective Hamiltonian, eqs 22 and 11, one can realize that the auxiliary coefficients C are nothing, but the amplitudes of the Löwdin-orthogonalized fragment states. Thus, the two formulations are complementary: if the computations of the vibronic Hamiltonian and the overlap matrix are more convenient to perform in the original fragments basis, |ψ⟩, the evolution of the amplitudes is more convenient in the Löwdin-transformed one, |ψ̃ L‑FMO⟩. 2.4. Surface Hopping and Decoherence. As it has been noted above, surface hopping calculations are more convenient to carry out in an orthonormal rather than in a nonorthogonal basis. First, the hopping in the nonorthogonal basis leads to the conceptual difficulties in assigning populations of such states. Second, the surface hopping probabilities would need to be redefined to account for implicit time dependence of the overlap matrix. On the contrary, the diabatic fragment states (orthogonal or nonorthogonal) are well-defined, because they are always localized on particular fragments. In addition, their character does not change upon nuclear evolution, solving the trivial crossing problem.88 As a compromise between all the considerations above, the orthogonal fragment states are chosen for surface hopping dynamics. However, to keep track the fragment depopulation dynamics, suitable projections of the time-evolved wave function onto fragment states of interest are computed:

Multiplying both sides of eq 16 by iℏS from the left, and recalling the TD-SE, eq 3b, one further gets iℏS

(22)

The formal solution of eq 22 is given by the matrix exponential:

From the definition, eq 15b, one obtains dc 1 dS dC = − S −3/2 C + S −1/2 dt 2 dt dt

dC = S −1/2H̃ vibS −1/2C dt

= S ·|c L̃ ‐FMO(t )|2 (21)

(25a)

If the propagation is performed in the |ψ̃ A‑FMO⟩ basis, and the |ψ̃ ⟩ states are being tracked, the matrix is defined as L‑FMO

With these definitions, eq 19 can be simplified: 5723


Article


Tretiak group.95 Analogous approaches to introduce decoherence were also discussed by Hammes-Schiffer and Tully86 and Landry and Subotnik.91 According to the technique, named “instantaneous decoherence at attempted hops” (ID-A), periodic adjustments of the coherently evolving electronicstate amplitudes are performed to mimic stochastic collapse of the wave function and to achieve a consistency with the TSH process. According to the technique, every attempted hop causes decoherence. If the hop is successful, the coherent superposition is collapsed onto the final state by setting cf = 1 (up to a complex phase factor) and ca = 0, ∀ a ≠ f. If the hop is unsuccessful, the coherent superposition is collapsed onto the currently running (“initial”) state: ci = 1 and ca = 0, ∀ a ≠ i. In the present work, the Boltzmann scaling technique is combined with the ID-A scheme, to account for the detailed balance and decoherence, respectively. Specifically, a random number, ξ1, is taken from a uniform distribution on the [0, 1] interval and is compared with the hopping probabilities, eq 26. If ξ1 > Pi→j, no transition occurs. The wave function, eq 1, continues to evolve coherently. If ξ1 < Pi→j, the corresponding transition is regarded as an attempted hop. Another random number, ξ2, is drawn from the same distribution and is compared with the Boltzmann factor, f b. If ξ2 < f b, the hop is regarded successful: the electronic transition is accepted and the coherent superposition is collapsed onto the final state. Alternatively, the hop is regarded unsuccessful (frustrated): the electronic transition is not accepted and the coherent superposition is collapsed back to the initial (currently running) state.

P L‐FMO(t ) = |⟨ψ |ψ̃ A‐FMO⟩cÃ ‐FMO(t )|2 = |⟨ψ |ψ ⟩U |2 |cÃ ‐FMO(t )|2 = U +S2U |cÃ ‐FMO(t )|2

(25b)

In analogy, when analyzing the statistics of the hopping trajectories, one should use fractional numbers based on the matrix PO(t) for assigning weights the trajectory contributes to different observed states. The time-evolved amplitudes of electronic states are used to compute the surface hopping dynamics. Unlike the commonly used fewest switches surface hopping (FSSH) probabilities,81 the present work utilizes the Markov-state surface hopping technique (MSSH).89 According to it, the probability to hop into a target state j from any source state i is determined by the populations of the target state: Pi → j = |ci|2

(26)

The scheme allows for several simultaneous transitions to occur such that the total result is consistent with the overall change of the quantum populations. Such a scheme has been shown to be particularly useful for modeling nonadiabatic dynamics in a diabatic basis.89 The electronic transitions that rely solely on eq 26 would not account for the detailed balance in the combined quantumclassical system. To ensure the balance, the following hop rejection scheme is applied. Assume the potential energies of the system before and after an attempted (resulting from the fortunate stochastic process realization) hop are Ei and Ef, respectively. The total energy must be conserved during the hop, resulting in the condition, Ei + Ti = Ef + Tf, where Ti and Tf are the initial and final kinetic energies of the system. In the situation when Ef > Ei + Ti, the final kinetic energy would be negative, Tf = Ei + Ti < 0. Therefore, such hops are not physical and are not classically allowed. No change of the discrete electronic state occurs. If the final kinetic energy is positive, the hops are accepted, but the final velocities are rescaled to accommodate the difference in the potential energies of the initial and final states. A simplified approach to the hop rejection has been utilized by the Prezhdo group.90 The hop rejection is decided on the basis of the Boltzmann factor consideration. Effectively, the hopping probabilities in eq 26 are scaled by the factor

(

(

fb = min 1, exp −

Ef − E i kBT

)),

3. RESULTS AND DISCUSSION 3.1. Molecular Structure and Nuclear Dynamics. To validate the developed methodology, it is applied to calculate CT rates in the recently studied SubPc/C60 heterojunction.96 The measurements suggested the CT from SubPc to C60 takes place on a 10 ps time scale, whereas the computations based on the phenomenological Fermi Golden Rule (FGR) predicted an overestimated rate, with the time scale on the order of 500 fs. The system is convenient for the present purposes because it consists of noncovalently bound molecular components, making the fragmentation scheme straightforward: SubPc and C60 molecules are treated as separate fragments. This partitioning also precludes the inclusion of many high-energy states that might appear when covalent bonds are broken. To demonstrate the applicability of the developed method to large systems and to examine the role of the model size choice in simulating interfacial CT dynamics, we consider two types of molecular models. The first type, further referred to as 1, consists of a single acceptor molecule, C60, and a single donor SubPc molecule (Figure 2a). This is a minimal computational model, also used by Wilcox et al.96 The second type of models, referred to as 2, contains one donor SubPc and ten C60 molecules: seven in the surface layer, and three in the bottom layer (Figure 2b). In the models of this type, all atoms are involved in nuclear dynamics driven by classical molecular mechanics interactions. However, only a subset of atoms (defined by the number of electronically active C60 units) is included in electronic structure and NA-MD calculations. In the SubPc/(C60)n system, the SubPc and n C60 fragments are considered electronically active, whereas all remaining fragments are needed only for structural (including nuclear dynamic) purposes. The models of type 2 can be regarded realistic representations of the interfacial heterostructure. The

to disfavor hops to higher

energy states. Although the total energy is not conserved in this treatment, the approach is computationally efficient, because it does not perturb the nuclear dynamics. Therefore, multiple realizations of the stochastic surface hopping can be computed for each single nuclear trajectory. Specifically, the present calculations employ 1000 realizations of the stochastic TSH process for each nuclear trajectory. To account for thermal averaging, 10 independent nuclear trajectories are computed for each system. Thus, the total number of the stochastic-thermal sampling trajectories is 10 000 for each system. In addition to the proper description of thermal equilibrium, it is important to account for electronic coherences. Both the TSH and the underlying MF dynamics are overcoherent, which may translate into the largely overestimated electronic transition rates.91−94 A simple way to account for electronic decoherence was comprehensively examined recently by the 5724


Article


(NVT) at 300 K. A total of 5000 steps are computed with the integration time step of 0.5 fs, resulting in a 2.5 ps thermalization period. The constant temperature is maintained by the Nose−Hoover chain thermostat.113−117 The thermalized systems are then used to continue NVE dynamics with the onthe-fly fragment NA-MD calculations. In the present calculations, 5000 and 10 000 0.5 fs steps (including NA-MD part) have been computed for the systems of type 2 and 1, respectively. All computations have been performed on a single CPU node with a limited memory requirement. In the present work, the neglect of back-reaction approximation has been adopted. According to it, the nuclear dynamics of the overall system governs the electronic dynamics via parametric dependence of the vibronic Hamiltonian on nuclear coordinates, whereas the electronic dynamics does not affect the nuclear one. This version of the classical path approximation (CPA)118 had been successfully utilized in the past.119−125 Its validity is justified by the fact that the molecular systems are sufficiently large and rigid. Specifically, nuclear forces are determined by the electronic structure of the overall system. Exciting a single electron in such a large system is not expected to strongly perturb nuclear evolution. Thus, as long as the nuclear dynamics is governed by the overall system and not by the individual fragments, the CPA is expected to be valid, especially as the size of the system increases. In the present work the nuclear dynamics is governed by a classical potential of the overall system. A similar approach to nuclear dynamics within the CPA NA-MD has been successfully utilized in a number of recent studies.126−128 Even though the fragmentation ansatz is introduced, it is not used to determine nuclear forces. Thus, the nuclear dynamics remains unaffected by the fragmentation scheme. As another consequence, the nonorthogonality of quasi-diabatic states does not introduce additional complications of the present computational procedure. The molecular mechanics-based CPA is an important acceleration approach, because it helps avoiding expensive computations of forces, especially those on excited states. 3.2. Static Electronic Structure and the Choice of Electronic Basis. The electronic structure calculations on each fragment are performed using the semiempirical EHT.129−131 EHT offers a chemically transparent and computationally efficient description of electronic states and their energies and was utilized by a number of authors.24,26,28,29,126−128 At the same time, a judicious parametrization is often required. In the present work, the Calzaferri131 formulation of the EHT has been utilized, according to which the electronic Hamiltonian in the AO basis can be expressed via:

Figure 2. Molecular models of the SubPc/C60 interface. (a) System of type 1, a single C60 is used to model the surface. (b) System of type 2, a cluster of ten C60 molecules is used to represent the surface. The color code for the atoms: light-blue, C; dark-blue, N; ochre, B; green, Cl; white, H. The CPK rendering model is used in (a) to represent the SubPc molecule whereas the vdW rendering model is used to represent the C60. The vdW rendering model is used for all molecules in (b).

initial molecular structures have been constructed using the IQmol software.97 The visualization and analysis of the MOs, molecular structures, and MD trajectories have been performed using the VMD package.98 Further optimization of the initially prepared structure has been performed using the Libra code99 as will be explained below. The nuclear dynamics of the systems is computed using classical MD simulations with the universal force field (UFF),100 which represents intermolecular interactions in the ground electronic state. Due to its generic nature, UFF is known to be less accurate than specialized force field. Nonetheless, it has been successfully utilized to study various molecular systems, also those involving C60 and other organic molecules.101−109 Importantly, because molecular dynamics is a method to sample a phase space rather than a potential energy surface, the intrinsic inaccuracies of the UFF may be comparable to the uncertainties in entropic terms that also present in specialized force fields. The nuclear dynamics computations are based on the system-wide interactions, thus including the interfragmental interactions. This procedure is separate from the electronic structure calculations, which depend on the type of fragments included in electronic dynamics. In principle, generalizations utilizing nuclear gradients computed with self-consistent FMOs fully quantummechanically are possible. Ground-state molecular dynamics and geometry optimization calculations of large systems based on such approach have been reported.38,40,43,110−112 However, such approaches extend beyond the scope of the present method. To mimic strong dispersion interactions due to π−π stacking, the default vdW interaction strength parameter of carbon atoms has been scaled by the factor of 3. If the original parameter is used, the SubPc unit desorbs from the surface of the (C60)n cluster, whereas it stays adsorbed for the whole duration of MD simulations, if the scaled parameter is utilized. To ensure a stable and meaningful integration, the molecular structures have been optimized using 20 cycles of simulated annealing. In each cycle, 500 microcanonical (NVE) MD steps are carried out with the 0.5 fs integration time step. In the end of each cycle, the atomic velocities are reset to zero, removing generated kinetic energy from the system, and moving the system to a lower total energy level. The annealing is followed by the isothermal MD

Hij =

1 KSij[Hii + Hjj] 2

(27a)

K = 1 + (0.75 + Δ2 − 0.75Δ4 ) exp(− 0.13· (rij − rij0)) (27b)

Δ=

Hii − Hjj Hii + Hjj

(27c)

Here, Sij is the overlap of two AOs i and j, Hii is the energy of the orbital i, typically taken to be negative of the valence-state ionization potential, rij is the actual separation between two AO centers, whereas r0ij is the reference value, computed using atomic radii. The latter are computed using properties of valence shell orbitals, as prescribed in the original paper.131 The 5725


Article


was explained in the Theory and Methods section, the use of such diabatic states leads to some inconveniences, prompting for alternative representations. Specifically, as has been discussed above, the adiabatic (Figure 4, left column) and

approach introduces a geometry-dependent EHT proportionality constant, eq 27b. In the past, the model was successfully applied to compute the excited states of metal−organic complexes.132 The EHT model parameters (Table 1) have been adjusted to reproduce the experimentally measured highest occupied MO Table 1. Adjusted Parameters of the EHT Model orbital

IP, eV

ζ, au−1

H, 1s B, 2s B, 2p C, 2s C, 2p N, 2s N, 2p Cl, 3s Cl, 3p

−13.6 −13.462 −8.432 −19.376 −14.072 −26.223 −14.841 −29.196 −13.780

1.3000 1.2650 1.1340 1.5760 1.4350 1.6760 1.5350 2.2500 1.9000

(HOMO)−lowest unoccupied MO (LUMO) gaps of isolated SubPc and C60 molecules, as well as the alignment of the LU orbitals of the two molecules with respect to each other.96 The projected densities of states (pDOS) for the two molecules computed using the optimized parameters are presented in Figure 3a. One can observe a favorable band alignment of the

Figure 4. Three sets of electronic basis states, corresponding to different representations: A-FMO (left column), A-MO (middle column), and L-FMO (right column). Each set includes LUMO (top) to LUMO+4 (bottom) orbitals.

Figure 3. (a) Projected density of states of the isolated molecules. Positive values on the y-axis correspond to the SubPc states, and negative, to the C60 states. The latter are reduced by the factor of 2, for clarity. (b) Orbital isosurfaces, corresponding to 0.015 au−3 isovalue. The upper row shows two closely spaced donor orbitals: LUMO and LUMO+1 of SubPc. The lower row shows tree degenerate P-type acceptor orbitals: LUMO, LUMO+1, and LUMO+2 of C60.

Löwdin-transformed (Figure 4, right column) FMOs are considered. In addition, reference NA-MD simulations are performed within the basis of adiabatic MOs of the overall system (Figure 4, middle column). Figure 4 shows the orbitals included in the basis of electronic states within which NA-MD calculations discussed below are carried out. The basis sets feature three degenerate LUMO states (in the internal counting convention denoted as LUMO, LUMO+1, and LUMO+2), corresponding to the acceptor states localized on C60. Two higher states (LUMO+3 and LUMO+4) represent the donor states and are mostly localized on the SubPc unit. As one can observe, all sets agree with each other very well in a qualitative sense. Because the wave function is defined up to its sign, one can observe occasional sign alternations, but the charge density distributions are very similar. A careful examination reviles a very close similarity of the A-FMO and L-FMO sets, showing definite localization of the orbitals of both types on the corresponding fragments, even though these are not just the orbitals of the isolated fragments. On the contrary, A-MOs LUMO+3 and LUMO+4 are notably delocalized across the C60 and SubPc fragments. Figure 4 suggests that the difference between A-FMO and LFMO might be relatively small. However, there is a significant conceptual difference. The Löwdin-transformed orbitals, due to

SubPc and C60 LU levels. The closely spaced LUMO and LUMO+1 levels of SubPc can act as donor states, whereas the degenerate LUMO, LUMO+1, and LUMO+2 orbitals of C60 constitute the acceptor states. The computed orbital isosurfaces are shown in Figure 3b. One can observe a P-like character of the triply degenerate acceptor state. The donor states show certain symmetry of the wave function localization on the lobes of the SubPc molecule. No significant wave function localization on B or Cl atoms is present in the LUMO and LUMO +1 states, unlike the HOMO (not shown), in which the charge density is localized closely around the center of the molecule, including on B and Cl atoms. The localization of the LUMO/ LUMO+1 donor states on the lobes of the molecule, and the cup-like geometry of SubPc maximize the overlaps of these states with the states of C60 molecule, to enable efficient charge transfer. The orbitals illustrated in Figure 3 are the noninteracting (quasidiabatic) orbitals of the SubPc and C60 subsystems. As 5726


Article

Journal of Chemical Theory and Computation the construction, resemble most closely the “raw” quasidiabatic states of the noninteracting fragments. Therefore, such states will remain localized mostly on the corresponding fragments. On the contrary, A-FMO are defined as the fragment orbitals that diagonalize the electronic Hamiltonian in the D-FMO basis. Thus, one can expect the A-FMOs to be composed of linear combinations of original quasidiabatic FMOs. The coefficients in this combination are related to energies of the corresponding quasidiabatic orbitals: the degeneracy of the latter may lead to A-FMOs that are well-delocalized over several fragments. As an illustration, Figure 5, shows several randomly selected A-FMOs and L-FMOs computed for the system of type 2 with

interaction (CI)-like character of the A-FMOs. Similar to standard CI wave functions, the A-FMOs may be considered a way to account for nondynamic correlation. One of the consequences is the possibility to describe extended charge carriers as well as long-range charge transport and tunneling. Within the present methodology, the originally constructed D-FMOs are approximated by the one-electron orbitals. As has been shown in earlier works,133,134 one can think in terms of excited Slater determinants constructed using such 1-electron orbitals. For the purposes of NA-MD simulations, the two descriptions would be nearly equivalent providing same recipes for computing NACs and, under additional approximations, same energetics. The latter may be different, however, if the many-body energy terms are included.99 The independent orbital (1-electron) approximation has been in use for a long time, mostly due to its efficiency.24,26,119,122,135−138 It is an inevitable approximation for large systems, where manyelectron wave functions are presently beyond the computational feasibility. The many-electron wave functions are often dominated by the contributions from one of two excited Slater determinants. However, even if this is not the case, and the many-electron wave function is a mixture of many terms, ΨI = ∑i CiI Φi , the couplings between two many-electron states, are bound by the couplings of the singly excited determinants:

Figure 5. Randomly selected A-FMOs (top row) and L-FMOs (bottom row) computed for the system 2 with nine electronically active C60 fragments.

||⟨Ψ|I Ĥ |ΨJ ⟩|| < || max⟨ΦI |Ĥ |ΦJ ⟩|| i

(28)

with a similar inequality applicable to the NAC magnitudes. Thus, the time scales computed using 1-electron wave functions are expected to be the lower bound for the time scales on the basis of the many-electron wave functions. It is interesting to note that the introduced A-FMOs not only constitute a convenient basis for NA-MD simulations but also introduce some many-body effects via a tight-binding configuration interaction. Indeed, the A-FMOs are constructed as linear combinations of the D-FMOs, with the latter describing various locally excited Slater determinants. The AFMO Hamiltonian (e.g., eqs 6 and 10) is therefore a CI Hamiltonian computed within a tight-binding approximation. The corresponding A-FMO states are linear combinations of different excitations on different fragments. The localization properties of the considered orbitals are important for understanding the corresponding NACs and electronic couplings. Another parameter playing a crucial role in the NA-MD is the energy of electronic state. Selecting different representations may, in principle, change the

9 C 60 fragments included in the electronic structure calculations. One can clearly see the that A-FMOs can be delocalized over 2 or 3 C60 fragments, with any pattern of localization (Figure 5, upper row). As the energy of the corresponding A-FMOs increases, the delocalization is expected to extended onto more that 2 or 3 fragments (not shown). On the contrary to A-FMOs, all L-MOs are primarily localized on a single fragment (Figure 5, bottom row). The extended character of the A-FMOs may be important for modeling the dynamics of extended charge carriers, those whose size may extend onto many molecules (e.g., in molecular solids like C60). As Figure 5 (top row) suggests, A-FMOs can delocalize onto several, even distant, fragments. Importantly, this localization is determined by the energies of all fragments and does not depend on their relative positions. If D-FMOs of two or more fragments are energetically close or resonant, the resulting A-FMO would be well extended over a system of interest. This result simply follows from the configuration

Figure 6. Evolution of the donor and acceptor energy levels corresponding to several types of basis states. 5727


Article


Figure 7. Kinetics of charge transfer in different bases and methods. (a, b, c) Coherent evolution of populations with the propagation performed in one of the 3 bases considered: (a) A-FMO; (b) A-MO; (c) L-FMO. For the FMO-based orbitals, the evolution of related FMO states is also shown. (d, e, f) Evolution of the SH populations computed with different leading representation: (d) A-FMO; (e) A-MO; (f) L-FMO. (g, h, i) Evolution of the populations of the donor states, computed with the inclusion of decoherence: (g) A-FMO; (h) A-MO; (i) L-FMO.

energetics of different transitions quite dramatically. To understand this factor, the evolution of energies of the donor and acceptor orbitals has been tracked (Figure 6). As one can observe, the differences are relatively minor, which can be attributed to good localization of the orbitals and their relatively weak interaction with each other. One notable feature is that the dispersion of the A-MO donor levels is somewhat larger than that of the A-FMO and L-FMO. On the contrary, the dispersion of the donor levels in the A-MO representation is smaller than that in the A-FMO or L-FMO representations. These differences are unlikely to introduce any unusual features in the dynamics of excited states. Importantly, the donor− acceptor energy gaps are comparable to each other in all three representations. Thus, from the energetics point of view, the three representations are equivalent. This property is likely explained by a suitable choice of the fragmentation scheme, which does not introduce any bond cuts and follows a natural partitioning of the system onto individual molecular fragments. 3.3. Electronic Dynamics in the SubPc/C60 Dimer. The electronic dynamics calculations combine on-the-fly construction of the quasi-diabatic fragment states and corresponding

vibronic Hamiltonian, the propagation of electronic amplitudes, and the surface hopping calculations with or without decoherence correction. First, the dynamics of excited states in the system type 1 is considered. The relatively small size of the system allows for a detailed investigation of several factors that may affect the dynamics. The NA-MD in this system can be performed in the A-MO basis, providing a reference for the fragment-based NA-MD calculations. For the present purpose of computing electron-transfer rates, the kinetics curves showing the total population on the donor group (SubPc) have been computed. All simulations resulted in a reasonable Gaussian or exponential interfragmental CT kinetics, as will be discussed further. The upper row of Figure 7 summarizes the total population of donor states evolving according to a fully coherent TD-SE. Both the A-FMO (Figure 7a) and A-MO (Figure 7b) populations decay very slowly, with comparable rates. The DFMO and L-FMO populations computed according to transformations eqs 24 and 25 from the A-FMO wave functions are shown. They are superimposed on the A-FMO populations, which represent actual “working” variables used in the surface 5728


Article

Journal of Chemical Theory and Computation hopping procedure. The D-FMO and L-FMO are shown here only to provide a sense of how the populations in different representations evolve along the evolution of the working variable. As one can observe, the coherent TD-SE populations in these three representations are very close to each other. In contrast to the decay of the A-FMO and A-MO states, the L-FMO states decay much faster. This acceleration is dictated by the nonvanishing electronic couplings between such states. Recall that, unlike the A-FMO or A-MO, the L-FMOs do not diagonalize the electronic Hamiltonian. Therefore, the electron is transferred by both adiabatic and nonadiabatic mechanisms. Moreover, because of such strong coupling, the donor states are not depopulated completely. Instead, an equilibrium between the donor and acceptor states is established. One should also note that because the L-FMO states are not the eigenstates of the electronic Hamiltonian, the analysis of the dynamics based on such states is not straightforward, because each L-FMO can correspond to a mixture of pure donor and acceptor states. Therefore, the “relaxation” time scales in this basis are ultrafast and correspond to the process of transforming one mixture of pure states into another. The second row of Figure 7 shows the evolution of the total population of the donor states, as computed by the MSSH method without the inclusion of decoherence. The relaxation

( t)

single exponential kinetics, P(t ) = exp − τ . The decrease is more drastic for the originally slow processes: if the A-MO time scales increase from 205 fs to 3.8 ps, the A-FMO time scales change from 380 fs to 23.4 ps. It is important to note that the differences in the A-FMO and A-MO time scales that are revealed in the MSSH + ID-A calculations point mostly to the differences in the decoherence treatment within the two approaches, and not to the differences in the intrinsic electronic dynamics. For the latter, the comparison of the corresponding coherent (Figure 7a vs Figure 7b) and overcoherent SH (Figure 7d vs Figure 7e) simulations is more suitable. These results indicate reasonable agreement of the dynamics produced using the A-FMO and A-MO basis sets. At the same time, the comparison of Figure 7g and Figure 7h, indicates that the treatment of decoherence may be quite different in the two bases. Despite the differences in the CT dynamics computed with the A-FMO and A-MO states, it is clear that decoherence plays an essential role in determining correct range of CT time scales. Namely, for the system type 1 the present calculations with the inclusion of decoherence suggest a 3.8−23.4 ps time scale, on the same order of magnitude as those inferred from the experiment.96 On the contrary, when decoherence is not included, the computed time scales are underestimated to be on the order of 205−380 fs. The FGR theory does partially account for decoherence effects, via the incorporated dependence on the spectral density. At the same time, the FGRcomputed time scales of approximately 500 fs are strongly underestimated with respect to the experimental values. This may indicate that the phenomenological decoherence treatment included in the FGR theory may be incomplete. Specifically, the parameters computed for the SubPc/C60 dimer and utilized in the calculations may be not suitable for the actual interfacial system. As it has been mentioned above, the time scales for CT in the A-MO basis may be somewhat underestimated with respect to CT in the A-FMO basis, because the initial A-MOs are notably delocalized onto the C60 fragment. The delocalized states are likely to possess longer decoherence times, because they may resemble the coherent superposition much closer than a pure localized state. Extended decoherence times facilitate faster CT because the intervals of relatively fast coherent dynamics are increased. As a consequence, the CT dynamics becomes much faster in the A-MO basis than in the A-FMO one, even when the decoherence correction is included. 3.4. Electronic Dynamics in SubPc/(C60)n Complexes. To demonstrate the real usability of the developed methodology, large molecular complexes of type 2 have been considered. Due to their size, these systems are expected to better mimic a realistic interfacial environment. Thus, it is expected that the experimental time scales should be compared to those computed for this type of molecular models rather than to those found for the simplistic system of type 1. Simulations based on the adiabatic MOs become prohibitively expensive for the systems of type 2, especially when the number of electronically active C60 groups is large. Thus, only the AFMO results are reported. Both the overcoherent MSSH and the decoherence-corrected MSSH + ID-A types of NA-MD are considered. Unlike the dimer model 1, the CT dynamics in the systems 2 computed using only MSSH follows the Gaussian kinetics,

( t)

follows the exponential decay kinetics, P(t ) = exp − τ , with

the time scales, τ, notably decreased in comparison to those of the coherent evolution. The depopulation in the basis of AFMOs proceeds on the time scale of 380 fs. Although this time scale is close to the 500 fs found according to FGR,96 there is no rigorous ground to compare the two values directly because the latter does account for some decoherence, whereas the former does not. The reference calculations using A-MO basis suggest somewhat accelerated CT, with the time scale of 205 fs. This increase can be attributed to the significant delocalization of the “donor” A-MOs onto both SubPc and C60 units, hence leading to stronger couplings between the donor and acceptor states due to significantly increased overlap of the C60-localized orbitals. One should be careful when comparing the two time scales. If the 380 fs corresponds to the transfer of an entire electron from the SubPc to the C60 fragment, the 205 fs corresponds to the transfer of only a fraction of electron, because the other fraction is already localized on the C60 (and hence the acceptor). A visual examination of the LUMO+3 and LUMO+4 of the overall system 1, suggests that these orbitals contain approximately a half of electron on the SubPc unit. If the transfer of a half electron takes 205 fs, the transfer of a full electron would require about 410 fs, the number close to the CT time scale obtained with the A-FMO basis. Thus, the computations with the A-FMO and A-MO bases yield comparable results, suggesting also that the fragmentation approach with the adiabatic fragment orbitals is a reasonable way to model NA-MD in large systems. As it is expected from the coherent evolution of the populations, the L-FMO produces much shorter time scales for CT, which can not be easily interpreted. Therefore, the following discussion will mostly focus on the A-FMO orbitals. The use of L-FMOs is not recommended. The third row in Figure 7 illustrates the dynamics of CT with decoherence correction included. As it is generally expected, the rates are decreased significantly with respect to those computed with MSSH. The relaxation is well described by the 5729


Article


(

P(t ) = exp −

numbers on the number of electronically active C60 molecules (Figure 8b). One can clearly see a very large change when one goes from the complex 1 to all other complexes 2. Interestingly, the CT time scales decrease dramatically, even though the ovecoherent time scales, dictated by the couplings, increase (e.g., Figure 8a). Thus, the effect can be attributed to decoherence, which may occur much slower in the systems of type 2 than in the dimer 1. This leads to a larger acceleration of the CT rates in 2 than in 1. Thus, even though the time scales computed for 1 using MSSH + ID-A are quite distinct from the experimentally measured times of approximately 10 ps, the time scales in the interfacial system 2 are in much better agreement with the experiment. This result highlights the importance of considering bigger systems, beyond the size of the minimalistic models. Another interesting trend observed in the MSSH + ID-A calculations is that, unlike the MSSH case, the CT rates slowly increase as the size of the substrate increases. Thus, the model corrected for decoherence effects can capture the dependence of CT rates on the density of acceptor states, at it is expected from the FGR. Unlike in the overcoherent description, in the decoherence-corrected treatment the system’s wave function is periodically collapsed onto a pure state randomly selected from a manifold of available states. The wave function collapse explicitly depends on the density of acceptor states, because the outcome of this stochastic process is determined by the probability to pick a given pure state. On a qualitative level, each stochastic collapse transforms the coherent wave function of the system into a particular state, with a well-defined identity. The identities are determined by the availability of various fragmental orbitals, making the wave function collapse and hence CT dependent on the density of acceptor states. In principle, the FGR accounts for some decoherence effects, in an effective way. For instance, spectral density, reorganization energy, and density of the acceptor states are all encoded into the FGR expression. Nonetheless, these factors are introduced via effective parameters, corresponding to static or ensemble-averaged structures, which may in addition be overly minimalistic. For instance, the effective parameters used by Wilcox et al.96 are computed for the dimer complex 1, in which the relative coordination of the donor and acceptor groups may be different from that in a larger interfacial system. The present calculations (Figure 8) indicate a dramatic difference between the minimalistic complex 1 and the larger interfacial systems 2, which may be more accurate models of the realistic material. This result suggests that the minimalistic model 1 may be insufficient for accurate modeling of CT in SubPc/C60 system, leading to the discrepancies between the FGR calculations and the experimental results, discussed by Wilcox et al.96 In addition, as has been mentioned above, the treatment of decoherence via the model FGR parameters and via the explicit atomistic simulations may lead to different results. As it is found in this study, the explicit atomistic models lead to the time scales ranging from 5 to 10 ps, in much closer agreement with the experimental value of 10 ps than the corresponding static FGR prediction. It is expected that using the FGR approach with the parameters derived for a larger model of interfacial system may lead to an improved description of CT rates in SubPc/C60 system, in analogy with the present real-time NAMD simulations. 3.5. Implementation and Scaling Characteristics. All calculations have been performed using the recently developed Libra package,99 which provides basic building blocks for the

2

( τt ) ), not the exponential. The change of the

regime from the single exponential to the Gaussian can be attributed to the presence of different types of acceptor sites (C60 molecules), those at different distances from the donor SubPc group. The inferred time scales are therefore almost twice as large as the exponential decay time scale for the dimer 1 (Figure 8a). Remarkably, the time scales vary around 500 fs,

Figure 8. Dependence of CT time scales in SubPc/(C60)n systems on the number of electronically active C60 fragments, n, computed with (a) MSSH, no decoherence; (b) MSSH + ID-A, with decoherence.

coincidentally approaching the reported FGR results.96 One can conclude that the time scales computed without decoherence correction are relatively insensitive to the number of C60 fragments in the substrate. The decreased rate for CT in the interfacial complexes 2 in comparison to the dimer 1 can be attributed to the details of geometry. In the latter, the SubPc and C60 units are tightly fit to each other, leading to strong intermolecular interactions, large orbital overlaps and large electronic and nonadiabatic couplings. On the contrary, additional C60 fragments present in the systems 2 induce a steric repulsion of the SubPc fragment from the surface, decreasing the couplings and hence the CT rates. At the same time, these results emphasize the importance of the interfacial structure of the SubPc/C60. Note that the nuclear dynamics of SubPc/(C60)n includes interaction of SubPc with all C60 units. Thus, the details of interfacial nuclear dynamics do not depend on the number of electronically active acceptor molecules, n. The insensitivity of the CT rates to the number of electronically active acceptor states in the overcoherent MSSH simulations is an interesting and somewhat unexpected observation. It is expected that the density of acceptor states should affect the CT rates directly. The density of states increases with the number of electronically active C60 molecules included in the NA-MD simulations. Yet, at the MSSH level, this factor does not affect the time scales (Figure 8a). This can be attributed to the overcoherence of the underlying SH scheme. According to the overcoherent picture, the electrons always “see” the same set of acceptor statesthe one defined by the coherent superposition. No distinction between the fragment-specific states is implied. The states originating from different fragments lose their identity, making the density of acceptor states to be effectively independent of the number of electronically active C60 fragments. The CT dynamics in the systems of type 2 computed using MSSH + ID-A follows a single exponential kinetics,

( t)

P(t ) = exp − τ . As it has been discussed already, the inclusion of the decoherence correction increases the CT time scales. Here, we focus on the dependence of these 5730


Article


linearized using distance cutoff criteria. In general, the electronic Hamiltonian includes long-ranged Coulombic and exchange terms. The expenses due to electrostatic terms can be overcome using methods like fast multipole moments (FMM),140−144 with specific generalizations to quantum distributions.142,145,146 Additional techniques such as the resolution of identity,147−150 combined with the singular value151 or Cholesky152 decompositions can help reduce the costs for multicenter integral evaluation. The integrals screening based on the Schwartz inequality153 can be another way of reducing the amount of necessary calculations.

algorithm development: functions for molecular system construction, atom type perception, molecular mechanics and classical MD, Nose−Hoover thermostat used for thermal equlibration, electronic structure calculations with EHT (in a minimal STO-3G basis), the routines for computing pDOS and orbital isosurfaces, the modules for coherent electronic propagation and for the trajectory surface hopping calculations. The resulting implementation of the entire fragment-based NAMD simulation protocol constitutes a moderate-size Python script, which can be utilized for further developments. The script is distributed as a part of the Libra open-source package.139 To test the scaling properties of the developed algorithm, the times required for electronic structure calculations on all fragments and for computation of nonadiabatic and electronic couplings at each nuclear iteration are estimated. These times are obtained for all SubPc/(C60)n systems, with n = 1, 3, 5, 7, and 9. As it is expected from the fragmentation ansatz, the total time for electronic structure calculations is found to scale linearly with the number of C60 fragments included in NA-MD simulations (Figure 9a). This number is a good descriptor of

4. CONCLUSIONS In this work, I have presented the methodology for NA-MD simulations with fragment molecular orbitals computed using tight-binding EHT Hamiltonian. The mathematical formulation for computing NA-MD in a nonorthogonal basis has been elaborated. The presence of the orbital overlap and the nonHermiticity of effective vibronic Hamiltonian can be dealt with by a suitable basis transformation. Specifically, Lö wdin orthogonalization helps casting the TD-SE equations into a convenient form. At the same time, Löwdin-transformed quasidiabatic states have been found to be relatively poor choice for running quantum-classical trajectory surface hopping simulations. The primary difficulty arises due to the mixed nature of such orbitals leading to significantly increased couplings and unclear physical nature of such states. The coherent dynamics in this basis does not lead to full relaxation and results in a fast equilibration instead. The use of another basis set (adiabatic fragment MOs) has been tested. The basis has been found to be most suitable and convenient for formulating fragment-based NA-MD protocol. By the construction, the A-FMO basis has a capacity to introduce many-body effects, providing more accurate wave functions, also convenient for NA-MD simulations in large systems. The comparison with the all-system adiabatic MO calculations suggested that both the mean-field and overcoherent MSSH descriptions in these bases yield comparable results, with the A-MO showing somewhat accelerated transitions. This discrepancy has been attributed to a delocalized and mixed nature of the A-MO states. Because the A-FMO and A-MO may have different properties, a certain care must be exercised when comparing the results obtained in the two types of simulations. The inclusion of decoherence correction at the ID-A level has been found to be able to induce significant differences in the time scales computed with the AFMO and A-MO states. This effect is attributed to difference in decoherence times of pure and mixed states. The developed approach has been implemented as a relatively compact yet generic Python script driven by the methodology development Libra backbone. The developed algorithm allows simulating several picoseconds of NA-MD in systems with over 600 atoms requiring only a single node and limited memory and time resources. The CPU time required for electronic structure has been found to scale linearly with the size of the system. The CPU time for computing vibronic Hamiltonian (electronic and nonadiabatic couplings) needed for NA-MD scales quadratically with the size of the system but can be further improved. The developed tight-binding fragment-based NA-MD technique has been applied to study CT in the extensive SubPc/(C60)n systems with n varying up to 10. The present calculations suggest the importance of the careful selection of

Figure 9. Scaling of the CPU time required for (a) electronic structure calculations on all fragments and for (b) computing electronic couplings and overlaps (blue) and nonadiabatic couplings (green) with the size of the system. Note that the former is linear, whereas the latter is quadratic in the number of electronically active C60 fragments. The R2 values for the above relationships are 0.9962 (a), 0.9777 (b, blue), and 0.9993 (b, green).

the system’s size, because the C60 is the largest fragment and the number of molecular orbitals included in the calculations is proportional the number of atoms. A slight deviation from the linearity is observed for n = 1, because for such system the size of C60 and SubPc units are comparable. Nonetheless, high values of the R2 factor indicate a good degree of linearity. Despite the linear scaling of the time for electronic structure calculations, the time for computing nonadiabatic and electronic couplings scales quadratically with the size of the system (Figure 9b), as expected. Such a scaling follows simply from the notion that the number of orbital pairs scales quadratically with the number of orbitals. In principle, this scaling can be improved further by additional screening of the orbital pairs for which these quantities should be computed. The NA couplings are related to orbital overlaps, which decay exponentially with the interatomic separation. Thus, NACs can be computed only for the orbital pairs whose centers are apart from each other no more than a predefined cutoff distance. The reduction of the scaling for electronic couplings can be more elaborate than that for computing NACs. In the present formulation, the tight-binding Hamiltonian is utilized. As such, all the electronic couplings decay quite rapidly as a function of interatomic distance. Thus, the present scheme can be 5731


Article


(8) Eder, D. Carbon Nanotube−Inorganic Hybrids. Chem. Rev. 2010, 110, 1348−1385. (9) Lee, K.; Mazare, A.; Schmuki, P. One-Dimensional Titanium Dioxide Nanomaterials: Nanotubes. Chem. Rev. 2014, 114, 9385− 9454. (10) Li, Z.; Nieto-Pescador, J.; Carson, A. J.; Blake, J. C.; Gundlach, L. Efficient Z-Scheme Charge Separation in Novel Vertically Aligned ZnO/CdSSe Nanotrees. Nanotechnology 2016, 27, 135401. (11) Blake, J. C.; Eldridge, P. S.; Gundlach, L. Spatial Variation in Carrier Dynamics along a Single CdSSe Nanowire. Chem. Phys. 2014, 442, 128−131. (12) Chan, W.-L.; Tritsch, J. R.; Zhu, X.-Y. Harvesting Singlet Fission for Solar Energy Conversion: One- versus Two-Electron Transfer from the Quantum Mechanical Superposition. J. Am. Chem. Soc. 2012, 134, 18295−18302. (13) Chan, W.-L.; Ligges, M.; Jailaubekov, A.; Kaake, L.; Miaja-Avila, L.; Zhu, X.-Y. Observing the Multiexciton State in Singlet Fission and Ensuing Ultrafast Multielectron Transfer. Science 2011, 334, 1541− 1545. (14) Pandit, B.; Jackson, N. E.; Zheng, T.; Fauvell, T. J.; Manley, E. F.; Orr, M.; Brown-Xu, S.; Yu, L.; Chen, L. X. Molecular Structure Controlled Transitions between Free-Charge Generation and Trap Formation in a Conjugated Copolymer Series. J. Phys. Chem. C 2016, 120, 4189−4198. (15) Schneider, A. M.; Lu, L.; Manley, E. F.; Zheng, T.; Sharapov, V.; Xu, T.; Marks, T. J.; Chen, L. X.; Yu, L. Wide Bandgap OPV Polymers Based on Pyridinonedithiophene Unit with Efficiency > 5%. Chem. Sci. 2015, 6, 4860−4866. (16) Spencer, S. D.; Bougher, C.; Heaphy, P. J.; Murcia, V. M.; Gallivan, C. P.; Monfette, A.; Andersen, J. D.; Cody, J. A.; Conrad, B. R.; Collison, C. J. The Effect of Controllable Thin Film Crystal Growth on the Aggregation of a Novel High Panchromaticity Squaraine Viable for Organic Solar Cells. Sol. Energy Mater. Sol. Cells 2013, 112, 202−208. (17) Spencer, S.; Cody, J.; Misture, S.; Cona, B.; Heaphy, P.; Rumbles, G.; Andersen, J.; Collison, C. Critical Electron Transfer Rates for Exciton Dissociation Governed by Extent of Crystallinity in Small Molecule Organic Photovoltaics. J. Phys. Chem. C 2014, 118, 14840−14847. (18) Bhosale, S.; Sisson, A. L.; Talukdar, P.; Furstenberg, A.; Banerji, N.; Vauthey, E.; Bollot, G.; Mareda, J.; Roger, C.; Wurthner, F.; Sakai, N.; Matile, S. Photoproduction of Proton Gradients with Pi-Stacked Fluorophore Scaffolds in Lipid Bilayers. Science 2006, 313, 84−86. (19) Bhosale, R.; Míšek, J.; Sakai, N.; Matile, S. Supramolecular N/pHeterojunction Photosystems with Oriented Multicolored Antiparallel Redox Gradients (OMARG-SHJs). Chem. Soc. Rev. 2010, 39, 138− 149. (20) Christensen, A. S.; Kubař, T.; Cui, Q.; Elstner, M. Semiempirical Quantum Mechanical Methods for Noncovalent Interactions for Chemical and Biochemical Applications. Chem. Rev. 2016, 116, 5301− 5337. (21) Bredow, T.; Jug, K. Theory and Range of Modern Semiempirical Molecular Orbital Methods. Theor. Chem. Acc. 2005, 113, 1−14. (22) Akimov, A. V.; Prezhdo, O. V. Large-Scale Computations in Chemistry: A Bird’s Eye View of a Vibrant Field. Chem. Rev. 2015, 115, 5797−5890. (23) Cohen, A. J.; Mori-Sánchez, P.; Yang, W. Challenges for Density Functional Theory. Chem. Rev. 2012, 112, 289−320. (24) Rego, L. G. C.; Batista, V. S. Quantum Dynamics Simulations of Interfacial Electron Transfer in Sensitized TiO2 Semiconductors. J. Am. Chem. Soc. 2003, 125, 7989−7997. (25) Abuabara, S. G.; Rego, L. G. C.; Batista, V. S. Influence of Thermal Fluctuations on Interfacial Electron Transfer in Functionalized TiO2 Semiconductors. J. Am. Chem. Soc. 2005, 127, 18234− 18242. (26) da Silva, R.; Rego, L. G. C.; Freire, J. A.; Rodriguez, J.; Laria, D.; Batista, V. S. Study of Redox Species and Oxygen Vacancy Defects at TiO2-Electrolyte Interfaces. J. Phys. Chem. C 2010, 114, 19433−19442.

the atomistic model. A simple dimer model, also used by other authors,96 has been found insufficient. It leads to increased CT rates in comparison to the larger interfacial system, if no decoherence treatment is included. The neglect of decoherence has been found to be responsible for a counterintuitive insensitivity of the CT time scales on the density of acceptor states in such calculations. On the contrary, the inclusion of decoherence is critical for the simulations to account for the density of acceptor states. Under such conditions, the minimalistic dimer model leads to underestimated CT rates, because the density of donor states is different from that in the interfacial heterostructure. Only the use of the extended model, most closely mimicking the structure and dynamics of the realistic interface, and the inclusion of decoherence allowed predicting the CT time scales of 5−10 ps, in reasonably good agreement with the experiment (10 ps), especially considering all the approximations made in the model. Thus, both the details of the interfacial structure and the atomistic timedomain treatment of decoherence and nonadiabatic CT are essential for obtaining accurate time scales. It is also suggested that the FGR model used earlier by Wilcox et al.96 may need to be refined by using the effective parameters computed for an extended system as opposed to the minimalistic dimer model.

■

AUTHOR INFORMATION

Corresponding Author

*Email: alexeyak@buffalo.edu. ORCID

Alexey V. Akimov: 0000-0002-7815-3731 Notes

The author declares no competing financial interest.

■

ACKNOWLEDGMENTS The author acknowledges the financial support from the University at Buffalo, The State University of New York startup package. Support of computations is provided by the Center for Computational Research at the University at Buffalo.

■

REFERENCES

(1) Akimov, A. V.; Neukirch, A. J.; Prezhdo, O. V. Theoretical Insights into Photoinduced Charge Transfer and Catalysis at Oxide Interfaces. Chem. Rev. 2013, 113, 4496−4565. (2) Carey, G. H.; Abdelhady, A. L.; Ning, Z.; Thon, S. M.; Bakr, O. M.; Sargent, E. H. Colloidal Quantum Dot Solar Cells. Chem. Rev. 2015, 115, 12732−12763. (3) Beard, M. C. Multiple Exciton Generation in Semiconductor Quantum Dots. J. Phys. Chem. Lett. 2011, 2, 1282−1288. (4) Nozik, A. J.; Beard, M. C.; Luther, J. M.; Law, M.; Ellingson, R. J.; Johnson, J. C. Semiconductor Quantum Dots and Quantum Dot Arrays and Applications of Multiple Exciton Generation to ThirdGeneration Photovoltaic Solar Cells. Chem. Rev. 2010, 110, 6873− 6890. (5) Kern, M. E.; Watson, D. F. Linker-Assisted Attachment of CdSe Quantum Dots to TiO2: Time- and Concentration-Dependent Adsorption, Agglomeration, and Sensitized Photocurrent. Langmuir 2014, 30, 13293−13300. (6) Dibbell, R. S.; Youker, D. G.; Watson, D. F. Excited-State Electron Transfer from CdS Quantum Dots to TiO2 Nanoparticles via Molecular Linkers with Phenylene Bridges. J. Phys. Chem. C 2009, 113, 18643−18651. (7) Wang, S.; Khafizov, M.; Tu, X.; Zheng, M.; Krauss, T. D. Multiple Exciton Generation in Single-Walled Carbon Nanotubes. Nano Lett. 2010, 10, 2381−2386. 5732


Article


Approximate Electrostatic Potential. Chem. Phys. Lett. 2002, 351, 475− 480. (45) Mo, Y.; Peyerimhoff, S. D. Theoretical Analysis of Electronic Delocalization. J. Chem. Phys. 1998, 109, 1687−1697. (46) Mo, Y.; Zhang, Y.; Gao, J. A Simple Electrostatic Model for Trisilylamine: Theoretical Examinations of the N→σ* Negative Hyperconjugation, Pπ →dπ Bonding, and Stereoelectronic Interaction. J. Am. Chem. Soc. 1999, 121, 5737−5742. (47) Mo, Y.; Gao, J.; Peyerimhoff, S. D. Energy Decomposition Analysis of Intermolecular Interactions Using a Block-Localized Wave Function Approach. J. Chem. Phys. 2000, 112, 5530−5538. (48) Cembran, A.; Song, L.; Mo, Y.; Gao, J. Block-Localized Density Functional Theory (BLDFT), Diabatic Coupling, and Their Use in Valence Bond Theory for Representing Reactive Potential Energy Surfaces. J. Chem. Theory Comput. 2009, 5, 2702−2716. (49) Mo, Y.; Gao, J. An Ab Initio Molecular Orbital-Valence Bond (MOVB) Method for Simulating Chemical Reactions in Solution. J. Phys. Chem. A 2000, 104, 3012−3020. (50) Mo, Y.; Gao, J. Ab Initio QM/MM Simulations with a Molecular Orbital-Valence Bond (MOVB) Method: Application to an SN2 Reaction in Water. J. Comput. Chem. 2000, 21, 1458−1469. (51) Song, L.; Gao, J. On the Construction of Diabatic and Adiabatic Potential Energy Surfaces Based on Ab Initio Valence Bond Theory †. J. Phys. Chem. A 2008, 112, 12925−12935. (52) Oberhofer, H.; Blumberger, J. Revisiting Electronic Couplings and Incoherent Hopping Models for Electron Transport in Crystalline C60 at Ambient Temperatures. Phys. Chem. Chem. Phys. 2012, 14, 13846−13852. (53) Gajdos, F.; Valner, S.; Hoffmann, F.; Spencer, J.; Breuer, M.; Kubas, A.; Dupuis, M.; Blumberger, J. Ultrafast Estimation of Electronic Couplings for Electron Transfer between π-Conjugated Organic Molecules. J. Chem. Theory Comput. 2014, 10, 4653−4660. (54) Wu, Q.; Van Voorhis, T. Direct Optimization Method to Study Constrained Systems within Density-Functional Theory. Phys. Rev. A: At., Mol., Opt. Phys. 2005, 72, 024502. (55) Wu, Q.; Van Voorhis, T. Constrained Density Functional Theory and Its Application in Long-Range Electron Transfer. J. Chem. Theory Comput. 2006, 2, 765−774. (56) Wu, Q.; Kaduk, B.; Van Voorhis, T. Constrained Density Functional Theory Based Configuration Interaction Improves the Prediction of Reaction Barrier Heights. J. Chem. Phys. 2009, 130, 034109. (57) van Voorhis, T.; Kowalczyk, T.; Kaduk, B.; Wang, L.-P.; Cheng, C.-L.; Wu, Q. The Diabatic Picture of Electron Transfer, Reaction Barriers, and Molecular Dynamics. Annu. Rev. Phys. Chem. 2010, 61, 149−170. (58) Yost, S. R.; Wang, L.-P.; Van Voorhis, T. Molecular Insight Into the Energy Levels at the Organic Donor/Acceptor Interface: A Quantum Mechanics/Molecular Mechanics Study. J. Phys. Chem. C 2011, 115, 14431−14436. (59) Sena, A. M. P.; Miyazaki, T.; Bowler, D. R. Linear Scaling Constrained Density Functional Theory in CONQUEST. J. Chem. Theory Comput. 2011, 7, 884−889. (60) Valeev, E. F.; Coropceanu, V.; da Silva Filho, D. A.; Salman, S.; Brédas, J.-L. Effect of Electronic Polarization on Charge-Transport Parameters in Molecular Organic Semiconductors. J. Am. Chem. Soc. 2006, 128, 9882−9886. (61) Lemaur, V.; da Silva Filho, D. A.; Coropceanu, V.; Lehmann, M.; Geerts, Y.; Piris, J.; Debije, M. G.; van de Craats, A. M.; Senthilkumar, K.; Siebbeles, L. D. A.; Warman, J. M.; Brédas, J.-L.; Cornil, J. Charge Transport Properties in Discotic Liquid Crystals: A Quantum-Chemical Insight into Structure−Property Relationships. J. Am. Chem. Soc. 2004, 126, 3271−3279. (62) Hutchison, G. R.; Ratner, M. A.; Marks, T. J. Intermolecular Charge Transfer between Heterocyclic Oligomers. Effects of Heteroatom and Molecular Packing on Hopping Transport in Organic Semiconductors. J. Am. Chem. Soc. 2005, 127, 16866−16881.

(27) Hoff, D. A.; da Silva, R.; Rego, L. G. C. Coupled Electron−Hole Quantum Dynamics on D−π−A Dye-Sensitized TiO2 Semiconductors. J. Phys. Chem. C 2012, 116, 21169−21178. (28) Akimov, A. V.; Prezhdo, O. V. Persistent Electronic Coherence Despite Rapid Loss of Electron−Nuclear Correlation. J. Phys. Chem. Lett. 2013, 4, 3857−3864. (29) Rego, L. G. C.; Hames, B. C.; Mazon, K. T.; Joswig, J.-O. Intramolecular Polarization Induces Electron−Hole Charge Separation in Light-Harvesting Molecular Triads. J. Phys. Chem. C 2014, 118, 126−134. (30) Pal, S.; Trivedi, D. J.; Akimov, A. V.; Aradi, B.; Frauenheim, T.; Prezhdo, O. V. Nonadiabatic Molecular Dynamics for Thousand Atom Systems: A Tight-Binding Approach towards PYXAID. J. Chem. Theory Comput. 2016, 12, 1436−1448. (31) Shao, Y.; Molnar, L. F.; Jung, Y.; Kussmann, J.; Ochsenfeld, C.; Brown, S. T.; Gilbert, A. T. B.; Slipchenko, L. V.; Levchenko, S. V.; O’Neill, D. P.; DiStasio, R. A., Jr; Lochan, R. C.; Wang, T.; Beran, G. J. O.; Besley, N. A.; Herbert, J. M.; Yeh Lin, C.; Van Voorhis, T.; Hung Chien, S.; Sodt, A.; Steele, R. P.; Rassolov, V. A.; Maslen, P. E.; Korambath, P. P.; Adamson, R. D.; Austin, B.; Baker, J.; Byrd, E. F. C.; Dachsel, H.; Doerksen, R. J.; Dreuw, A.; Dunietz, B. D.; Dutoi, A. D.; Furlani, T. R.; Gwaltney, S. R.; Heyden, A.; Hirata, S.; Hsu, C.-P.; Kedziora, G.; Khalliulin, R. Z.; Klunzinger, P.; Lee, A. M.; Lee, M. S.; Liang, W.; Lotan, I.; Nair, N.; Peters, B.; Proynov, E. I.; Pieniazek, P. A.; Min Rhee, Y.; Ritchie, J.; Rosta, E.; David Sherrill, C.; Simmonett, A. C.; Subotnik, J. E.; Lee Woodcock, H., III; Zhang, W.; Bell, A. T.; Chakraborty, A. K.; Chipman, D. M.; Keil, F. J.; Warshel, A.; Hehre, W. J.; Schaefer, H. F., III; Kong, J.; Krylov, A. I.; Gill, P. M. W.; HeadGordon, M. Advances in Methods and Algorithms in a Modern Quantum Chemistry Program Package. Phys. Chem. Chem. Phys. 2006, 8, 3172−3191. (32) Beran, G. J. O.; Hirata, S. Fragment and Localized Orbital Methods in Electronic Structure Theory. Phys. Chem. Chem. Phys. 2012, 14, 7559−7561. (33) 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. (34) Pruitt, S. R.; Bertoni, C.; Brorsen, K. R.; Gordon, M. S. Efficient and Accurate Fragmentation Methods. Acc. Chem. Res. 2014, 47, 2786−2794. (35) Yang, W. Direct Calculation of Electron Density in DensityFunctional Theory. Phys. Rev. Lett. 1991, 66, 1438−1441. (36) Yang, W. Direct Calculation of Electron Density in DensityFunctional Theory: Implementation for Benzene and a Tetrapeptide. Phys. Rev. A: At., Mol., Opt. Phys. 1991, 44, 7823−7826. (37) Yang, W.; Lee, T.-S. A Density-Matrix Divide-and-Conquer Approach for Electronic Structure Calculations of Large Molecules. J. Chem. Phys. 1995, 103, 5674−5678. (38) Lee, T.-S.; York, D. M.; Yang, W. Linear-Scaling Semiempirical Quantum Calculations for Macromolecules. J. Chem. Phys. 1996, 105, 2744−2750. (39) Dixon, S. L.; Merz, K. M. Semiempirical Molecular Orbital Calculations with Linear System Size Scaling. J. Chem. Phys. 1996, 104, 6643−6649. (40) Dixon, S. L.; Merz, K. M. Fast, Accurate Semiempirical Molecular Orbital Calculations for Macromolecules. J. Chem. Phys. 1997, 107, 879−893. (41) 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. (42) Fedorov, D. G.; Ishida, T.; Kitaura, K. Multilayer Formulation of the Fragment Molecular Orbital Method (FMO). J. Phys. Chem. A 2005, 109, 2638−2646. (43) Fedorov, D. G.; Asada, N.; Nakanishi, I.; Kitaura, K. The Use of Many-Body Expansions and Geometry Optimizations in FragmentBased Methods. Acc. Chem. Res. 2014, 47, 2846−2856. (44) Nakano, T.; Kaminuma, T.; Sato, T.; Fukuzawa, K.; Akiyama, Y.; Uebayasi, M.; Kitaura, K. Fragment Molecular Orbital Method: Use of 5733


Article

Journal of Chemical Theory and Computation (63) Orlandi, G.; Troisi, A.; Zerbetto, F. Simulation of STM Images from Commercially Available Software. J. Am. Chem. Soc. 1999, 121, 5392−5395. (64) Troisi, A.; Orlandi, G. Band Structure of the Four Pentacene Polymorphs and Effect on the Hole Mobility at Low Temperature. J. Phys. Chem. B 2005, 109, 1849−1856. (65) Migliore, A. Nonorthogonality Problem and Effective Electronic Coupling Calculation: Application to Charge Transfer in π-Stacks Relevant to Biochemistry and Molecular Electronics. J. Chem. Theory Comput. 2011, 7, 1712−1725. (66) Kondov, I.; Cizek, M.; Benesch, C.; Wang, H.; Thoss, M. Quantum Dynamics of Photoinduced Electron-Transfer Reactions in Dye-Semiconductor Systems: First-Principles Description and Application to Coumarin 343-TiO2. J. Phys. Chem. C 2007, 111, 11970− 11981. (67) Mo, Y. Two-State Model Based on the Block-Localized Wave Function Method. J. Chem. Phys. 2007, 126, 224104. (68) Mo, Y.; Song, L.; Lin, Y.; Liu, M.; Cao, Z.; Wu, W. BlockLocalized Wavefunction (BLW) Based Two-State Approach for Charge Transfers between Phenyl Rings. J. Chem. Theory Comput. 2012, 8, 800−805. (69) Cembran, A.; Payaka, A.; Lin, Y.; Xie, W.; Mo, Y.; Song, L.; Gao, J. A Non-Orthogonal Block-Localized Effective Hamiltonian Approach for Chemical and Enzymatic Reactions. J. Chem. Theory Comput. 2010, 6, 2242−2251. (70) Oberhofer, H.; Blumberger, J. Electronic Coupling Matrix Elements from Charge Constrained Density Functional Theory Calculations Using a Plane Wave Basis Set. J. Chem. Phys. 2010, 133, 244105. (71) Gajdos, F.; Oberhofer, H.; Dupuis, M.; Blumberger, J. On the Inapplicability of Electron-Hopping Models for the Organic Semiconductor Phenyl-C61-Butyric Acid Methyl Ester (PCBM). J. Phys. Chem. Lett. 2013, 4, 1012−1017. (72) Kubas, A.; Gajdos, F.; Heck, A.; Oberhofer, H.; Elstner, M.; Blumberger, J. Electronic Couplings for Molecular Charge Transfer: Benchmarking CDFT, FODFT and FODFTB against High-Level Ab Initio Calculations. II. Phys. Chem. Chem. Phys. 2015, 17, 14342− 14354. (73) Zhang, J.; Valeev, E. F. Hybrid One-Electron/many-Electron Methods for Ionized States of Molecular Clusters. Phys. Chem. Chem. Phys. 2012, 14, 7863−7871. (74) Voityuk, A. A. Estimation of Electronic Coupling for Singlet Excitation Energy Transfer. J. Phys. Chem. C 2014, 118, 1478−1483. (75) Voityuk, A. A. Fragment Transition Density Method to Calculate Electronic Coupling for Excitation Energy Transfer. J. Chem. Phys. 2014, 140, 244117. (76) Voityuk, A. A. Electronic Coupling for Charge Transfer in Donor−bridge−acceptor Systems. Performance of the Two-State FCD Model. Phys. Chem. Chem. Phys. 2012, 14, 13789. (77) Wang, L.; Olivier, Y.; Prezhdo, O. V.; Beljonne, D. Maximizing Singlet Fission by Intermolecular Packing. J. Phys. Chem. Lett. 2014, 5, 3345−3353. (78) Xu, X.; Yang, K. R.; Truhlar, D. G. Diabatic Molecular Orbitals, Potential Energies, and Potential Energy Surface Couplings by the 4Fold Way for Photodissociation of Phenol. J. Chem. Theory Comput. 2013, 9, 3612−3625. (79) Yang, K. R.; Xu, X.; Truhlar, D. G. Direct Diabatization of Electronic States by the Fourfold-Way: Including Dynamical Correlation by Multi-Configuration Quasidegenerate Perturbation Theory with Complete Active Space Self-Consistent-Field Diabatic Molecular Orbitals. Chem. Phys. Lett. 2013, 573, 84−89. (80) Hoyer, C. E.; Xu, X.; Ma, D.; Gagliardi, L.; Truhlar, D. G. Diabatization Based on the Dipole and Quadrupole: The DQ Method. J. Chem. Phys. 2014, 141, 114104. (81) Tully, J. C. Molecular Dynamics with Electronic Transitions. J. Chem. Phys. 1990, 93, 1061−1071. (82) Löwdin, P.-O. Some Properties of Inner Projections. Int. J. Quantum Chem. 1970, 5, 231−237.

(83) Ren, H.; Provorse, M. R.; Bao, P.; Qu, Z.; Gao, J. Multistate Density Functional Theory for Effective Diabatic Electronic Coupling. J. Phys. Chem. Lett. 2016, 7, 2286−2293. (84) Azpiroz, J. M.; Infante, I.; De Angelis, F. First-Principles Modeling of Core/Shell Quantum Dot Sensitized Solar Cells. J. Phys. Chem. C 2015, 119, 12739−12748. (85) Iron, M. A.; Heyden, A.; Staszewska, G.; Truhlar, D. G. TightBinding Configuration Interaction (TBCI): A Noniterative Approach to Incorporating Electrostatics into Tight Binding. J. Chem. Theory Comput. 2008, 4, 804−818. (86) Hammes-Schiffer, S.; Tully, J. C. Proton Transfer in Solution: Molecular Dynamics with Quantum Transitions. J. Chem. Phys. 1994, 101, 4657−4667. (87) Akimov, A. V.; Prezhdo, O. V. Advanced Capabilities of the PYXAID Program: Integration Schemes, Decoherence Effects, Multiexcitonic States, and Field-Matter Interaction. J. Chem. Theory Comput. 2014, 10, 789−804. (88) Wang, L.; Prezhdo, O. V. A Simple Solution to the Trivial Crossing Problem in Surface Hopping. J. Phys. Chem. Lett. 2014, 5, 713−719. (89) Akimov, A. V.; Trivedi, D.; Wang, L.; Prezhdo, O. V. Analysis of the Trajectory Surface Hopping Method from the Markov State Model Perspective. J. Phys. Soc. Jpn. 2015, 84, 094002. (90) Duncan, W. R.; Craig, C. F.; Prezhdo, O. V. Time-Domain Ab Initio Study of Charge Relaxation and Recombination in DyeSensitized TiO2. J. Am. Chem. Soc. 2007, 129, 8528−8543. (91) Landry, B. R.; Subotnik, J. E. Communication: Standard Surface Hopping Predicts Incorrect Scaling for Marcus’ Golden-Rule Rate: The Decoherence Problem Cannot Be Ignored. J. Chem. Phys. 2011, 135, 191101−191101. (92) Jaeger, H. M.; Fischer, S.; Prezhdo, O. V. Decoherence-Induced Surface Hopping. J. Chem. Phys. 2012, 137, 22A545. (93) Akimov, A. V.; Long, R.; Prezhdo, O. V. Coherence Penalty Functional: A Simple Method for Adding Decoherence in Ehrenfest Dynamics. J. Chem. Phys. 2014, 140, 194107. (94) Akimov, A. V.; Prezhdo, O. V. Nonradiative Relaxation of Charge Carriers in GaN-InN Alloys: Insights from Nonadiabatic Molecular Dynamics. ACS Symp. Ser. 2015, 1196, 189−200. (95) Nelson, T.; Fernandez-Alberti, S.; Roitberg, A. E.; Tretiak, S. Nonadiabatic Excited-State Molecular Dynamics: Treatment of Electronic Decoherence. J. Chem. Phys. 2013, 138, 224111. (96) Wilcox, D. E.; Lee, M. H.; Sykes, M. E.; Niedringhaus, A.; Geva, E.; Dunietz, B. D.; Shtein, M.; Ogilvie, J. P. Ultrafast Charge-Transfer Dynamics at the Boron Subphthalocyanine Chloride/C60 Heterojunction: Comparison between Experiment and Theory. J. Phys. Chem. Lett. 2015, 6, 569−575. (97) IQmol Molecular Viewer http://iqmol.org/ (accessed May 18, 2016). (98) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33−38. (99) Akimov, A. V. Libra: An Open-Source “methodology Discovery” Library for Quantum and Classical Dynamics Simulations. J. Comput. Chem. 2016, 37, 1626−1649. (100) Rappe, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard, W. A., III; Skiff, W. M. UFF, a Full Periodic Table Force Field for Molecular Mechanics and Molecular Dynamics Simulations. J. Am. Chem. Soc. 1992, 114, 10024−10035. (101) Dzubak, A. L.; Lin, L.-C.; Kim, J.; Swisher, J. A.; Poloni, R.; Maximoff, S. N.; Smit, B.; Gagliardi, L. Ab Initio Carbon Capture in Open-Site Metal-Organic Frameworks. Nat. Chem. 2012, 4, 810−816. (102) Höft, N.; Horbach, J. Condensation of Methane in the Metal− Organic Framework IRMOF-1: Evidence for Two Critical Points. J. Am. Chem. Soc. 2015, 137, 10199−10204. (103) Bury, W.; Fairen-Jimenez, D.; Lalonde, M. B.; Snurr, R. Q.; Farha, O. K.; Hupp, J. T. Control over Catenation in Pillared Paddlewheel Metal−Organic Framework Materials via SolventAssisted Linker Exchange. Chem. Mater. 2013, 25, 739−744. 5734


Article

Journal of Chemical Theory and Computation (104) Wilmer, C. E.; Leaf, M.; Lee, C. Y.; Farha, O. K.; Hauser, B. G.; Hupp, J. T.; Snurr, R. Q. Large-Scale Screening of Hypothetical Metal−organic Frameworks. Nat. Chem. 2011, 4, 83−89. (105) Skoulidas, A. I.; Sholl, D. S. Self-Diffusion and Transport Diffusion of Light Gases in Metal-Organic Framework Materials Assessed Using Molecular Dynamics Simulations. J. Phys. Chem. B 2005, 109, 15760−15768. (106) Akimov, A. V.; Kolomeisky, A. B. Unidirectional Rolling Motion of Nanocars Induced by Electric Field. J. Phys. Chem. C 2012, 116, 22595−22601. (107) Akimov, A. V.; Kolomeisky, A. B. Molecular Dynamics Study of Crystalline Molecular Gyroscopes. J. Phys. Chem. C 2011, 115, 13584−13591. (108) Konyukhov, S. S.; Kupchenko, I. V.; Moskovsky, A. A.; Nemukhin, A. V.; Akimov, A. V.; Kolomeisky, A. B. Rigid-Body Molecular Dynamics of Fullerene-Based Nanocars on Metallic Surfaces. J. Chem. Theory Comput. 2010, 6, 2581−2590. (109) Akimov, A. V.; Williams, C.; Kolomeisky, A. B. Charge Transfer and Chemisorption of Fullerene Molecules on Metal Surfaces: Application to Dynamics of Nanocars. J. Phys. Chem. C 2012, 116, 13816−13826. (110) Ganesh, V.; Dongare, R. K.; Balanarayan, P.; Gadre, S. R. Molecular Tailoring Approach for Geometry Optimization of Large Molecules: Energy Evaluation and Parallelization Strategies. J. Chem. Phys. 2006, 125, 104109. (111) Hu, X.; Jin, Y.; Zeng, X.; Hu, H.; Yang, W. Liquid Water Simulations with the Density Fragment Interaction Approach. Phys. Chem. Chem. Phys. 2012, 14, 7700. (112) Chiba, M.; Fedorov, D. G.; Nagata, T.; Kitaura, K. Excited State Geometry Optimizations by Time-Dependent Density Functional Theory Based on the Fragment Molecular Orbital Method. Chem. Phys. Lett. 2009, 474, 227−232. (113) Nose, S. A Unified Formulation of the Constant Temperature Molecular Dynamics Methods. J. Chem. Phys. 1984, 81, 511−519. (114) Hoover, W. G. Generalization of Nose’s Isothermal Molecular Dynamics: Non-Hamiltionian Dynamics for the Canonical Ensemble. Phys. Rev. A: At., Mol., Opt. Phys. 1989, 40, 2814−2815. (115) Nose, S.; Klein, M. L. Constant-Temperature-ConstantPressure Molecular-Dynamics Calculations for Molecular Solids: Application to Solid Nitrogen at High Pressure. Phys. Rev. B: Condens. Matter Mater. Phys. 1986, 33, 339−342. (116) Kleinerman, D. S.; Czaplewski, C.; Liwo, A.; Scheraga, H. A. Implementations of Nosé−Hoover and Nosé−Poincaré Thermostats in Mesoscopic Dynamic Simulations with the United-Residue Model of a Polypeptide Chain. J. Chem. Phys. 2008, 128, 245103. (117) Kamberaj, H.; Low, R. J.; Neal, M. P. Time Reversible and Symplectic Integrators for Molecular Dynamics Simulations of Rigid Molecules. J. Chem. Phys. 2005, 122, 224114. (118) Augustin, S. D.; Rabitz, H. The Classical Path Approximation in Time-Dependent Quantum Collision Theory. J. Chem. Phys. 1978, 69, 4195. (119) Kilina, S. V.; Craig, C. F.; Kilin, D. S.; Prezhdo, O. V. Ab Initio Time-Domain Study of Phonon-Assisted Relaxation of Charge Carriers in a PbSe Quantum Dot. J. Phys. Chem. C 2007, 111, 4871−4878. (120) Isborn, C. M.; Kilina, S. V.; Li, X.; Prezhdo, O. V. Generation of Multiple Excitons in PbSe and CdSe Quantum Dots by Direct Photoexcitation: First-Principles Calculations on Small PbSe and CdSe Clusters. J. Phys. Chem. C 2008, 112, 18291−18294. (121) Kilin, D. S.; Tsemekhman, K. L.; Kilina, S. V.; Balatsky, A. V.; Prezhdo, O. V. Photoinduced Conductivity of a Porphyrin−Gold Composite Nanowire. J. Phys. Chem. A 2009, 113, 4549−4556. (122) Kilina, S. V.; Neukirch, A. J.; Habenicht, B. F.; Kilin, D. S.; Prezhdo, O. V. Quantum Zeno Effect Rationalizes the Phonon Bottleneck in Semiconductor Quantum Dots. Phys. Rev. Lett. 2013, 110, 180404. (123) Chen, J.; Schmitz, A.; Inerbaev, T.; Meng, Q.; Kilina, S.; Tretiak, S.; Kilin, D. S. First-Principles Study of P-N-Doped Silicon

Quantum Dots: Charge Transfer, Energy Dissipation, and TimeResolved Emission. J. Phys. Chem. Lett. 2013, 4, 2906−2913. (124) Long, R.; Fang, W.; Akimov, A. V. Nonradiative Electron− Hole Recombination Rate Is Greatly Reduced by Defects in Monolayer Black Phosphorus: Ab Initio Time Domain Study. J. Phys. Chem. Lett. 2016, 7, 653−659. (125) Madjet, M. E.-A.; Akimov, A. V.; El-Mellouhi, F.; Berdiyorov, G. R.; Ashhab, S.; Tabet, N.; Kais, S. Enhancing the Carrier Thermalization Time in Organometallic Perovskites by Halide Mixing. Phys. Chem. Chem. Phys. 2016, 18, 5219−5231. (126) Torres, A.; Oliboni, R. S.; Rego, L. G. C. Vibronic and Coherent Effects on Interfacial Electron Transfer Dynamics. J. Phys. Chem. Lett. 2015, 6, 4927−4935. (127) da Silva, R.; Hoff, D. A.; Rego, L. G. C. Coupled QuantumClassical Method for Long Range Charge Transfer: Relevance of the Nuclear Motion to the Quantum Electron Dynamics. J. Phys.: Condens. Matter 2015, 27, 134206. (128) Akimov, A. V.; Asahi, R.; Jinnouchi, R.; Prezhdo, O. V. What Makes the Photocatalytic CO2 Reduction on N-Doped Ta 2O5 Efficient: Insights from Nonadiabatic Molecular Dynamics. J. Am. Chem. Soc. 2015, 137, 11517−11525. (129) Hoffmann, R. An Extended Hückel Theory. I. Hydrocarbons. J. Chem. Phys. 1963, 39, 1397−1412. (130) Ammeter, J. H.; Bürgi, H. B.; Thibeault, J. C.; Hoffmann, R. Counterintuitive Orbital Mixing in Semiempirical and Ab Initio Molecular Orbital Calculations. J. Am. Chem. Soc. 1978, 100, 3686− 3692. (131) Calzaferri, G.; Forss, L.; Kamber, I. Molecular Geometries by the Extended Hueckel Molecular Orbital (EHMO) Method. J. Phys. Chem. 1989, 93, 5366−5371. (132) Amouyal, E.; Mouallem-Bahout, M.; Calzaferri, G. Excited States of M(II,d6)-4′-Phenylterpyridine Complexes: Electron Localization. J. Phys. Chem. 1991, 95, 7641−7649. (133) Akimov, A. V.; Prezhdo, O. V. The PYXAID Program for NonAdiabatic Molecular Dynamics in Condensed Matter Systems. J. Chem. Theory Comput. 2013, 9, 4959−4972. (134) Ryabinkin, I. G.; Nagesh, J.; Izmaylov, A. F. Fast Numerical Evaluation of Time-Derivative Non-Adiabatic Couplings for Mixed Quantum-Classical Methods. J. Phys. Chem. Lett. 2015, 6, 4200−4203. (135) Shenvi, N.; Roy, S.; Tully, J. C. Dynamical Steering and Electronic Excitation in NO Scattering from a Gold Surface. Science 2009, 326, 829−832. (136) Krishna, V.; Tully, J. C. Vibrational Lifetimes of Molecular Adsorbates on Metal Surfaces. J. Chem. Phys. 2006, 125, 054706. (137) Hyeon-Deuk, K.; Prezhdo, O. V. Time-Domain Ab Initio Study of Auger and Phonon-Assisted Auger Processes in a Semiconductor Quantum Dot. Nano Lett. 2011, 11, 1845−1850. (138) Fischer, S. A.; Habenicht, B. F.; Madrid, A. B.; Duncan, W. R.; Prezhdo, O. V.; et al. Regarding the Validity of the Time-Dependent Kohn-Sham Approach for Electron-Nuclear Dynamics via Trajectory Surface Hopping. J. Chem. Phys. 2011, 134, 024102. (139) alexvakimov/libra-code https://github.com/alexvakimov/libracode (accessed Nov. 4, 2016). (140) Ochsenfeld, C.; Kussmann, J.; Lambrecht, D. S. Linear-Scaling Methods in Quantum Chemistry. Rev. Comput. Chem. 2007, 23, 1−82. (141) Greengard, L.; Rokhlin, V. A Fast Algorithm for Particle Simulations. J. Comput. Phys. 1987, 73, 325−348. (142) White, C. A.; Johnson, B. G.; Gill, P. M. W.; Head-Gordon, M. The Continuous Fast Multipole Method. Chem. Phys. Lett. 1994, 230, 8−16. (143) White, C. A.; Johnson, B. G.; Gill, P. M. W.; Head-Gordon, M. Linear Scaling Density Functional Calculations via the Continuous Fast Multipole Method. Chem. Phys. Lett. 1996, 253, 268−278. (144) Sun, X.; Pitsianis, N. P. A Matrix Version of the Fast Multipole Method. SIAM Rev. 2001, 43, 289−300. (145) Kutteh, R.; Apra, E.; Nichols, J. A Generalized Fast Multipole Approach for Hartree-Fock and Density Functional Computations. Chem. Phys. Lett. 1995, 238, 173−179. 5735


Article

Journal of Chemical Theory and Computation (146) Burant, J. C.; Strain, M. C.; Scuseria, G. E.; Frisch, M. J. Analytic Energy Gradients for the Gaussian Very Fast Multipole Method (GvFMM). Chem. Phys. Lett. 1996, 248, 43−49. (147) Dunlap, B. I.; Connolly, J. W. D.; Sabin, J. R. On First-Row Diatomic Molecules and Local Density Models. J. Chem. Phys. 1979, 71, 4993−4999. (148) Dunlap, B. I. Accurate Density Functional Calculation on Large Systems. Int. J. Quantum Chem. 1996, 58, 123−132. (149) Ren, X.; Rinke, P.; Blum, V.; Wieferink, J.; Tkatchenko, A.; Sanfilippo, A.; Reuter, K.; Scheffler, M. Resolution-of-Identity Approach to Hartree−Fock, Hybrid Density Functionals, RPA, MP2 and GW with Numeric Atom-Centered Orbital Basis Functions. New J. Phys. 2012, 14, 053020. (150) Sodt, A.; Subotnik, J. E.; Head-Gordon, M. Linear Scaling Density Fitting. J. Chem. Phys. 2006, 125, 194109. (151) Foerster, D. Elimination, in Electronic Structure Calculations, of Redundant Orbital Products. J. Chem. Phys. 2008, 128, 034108. (152) Beebe, N. H.; Linderberg, J. Simplifications in the Generation and Transformation of Two-Electron Integrals in Molecular Calculations. Int. J. Quantum Chem. 1977, 12, 683−705. (153) Aquilante, F.; Pedersen, T. B.; Lindh, R. Low-Cost Evaluation of the Exchange Fock Matrix from Cholesky and Density Fitting Representations of the Electron Repulsion Integrals. J. Chem. Phys. 2007, 126, 194106.

5736


Nonadiabatic Molecular Dynamics with Tight-Binding Fragment

Recommend Documents