Fragment Molecular Orbital Nonadiabatic Molecular Dynamics for

Article pubs.acs.org/JPCA

Fragment Molecular Orbital Nonadiabatic Molecular Dynamics for Condensed Phase Systems Ben Nebgen and Oleg V. Prezhdo* Department of Chemistry, University of Southern California, Los Angeles, California 90037, United States ABSTRACT: A method for efficiently simulating nonadiabatic molecular dynamics (NAMD) of nanoscale and condensed phase systems is developed and tested. The electronic structure, including force and nonadiabatic coupling, are obtained with the fragment molecular orbital (FMO) approximation, which provides significant computational savings by splitting the system into fragments and computing electronic properties of each fragment subject to the external field due to other all other fragments. The efficiency of the developed technique is demonstrated by studying the effect of explicit solvent molecules on excited state relaxation in the Fe(CO)4 complex. The relaxation in the gas phase occurs on a 50 fs time scale, which is in excellent agreement with previously recorded femtosecond pump−probe spectroscopy. Adding a solvation shell of ethanol molecules to the simulation results in an increase in the excited state lifetime to 100 fs, in agreement with recent femtosecond X-ray spectroscopy measurements.

1. INTRODUCTION The surface hopping (SH) formulation of quantum-classical nonadiabatic (NA) dynamics, as first explored by Tully in 1971,1 has become an important tool in the modern computational chemistry toolbox. It is seeing increased use in describing electron motion in systems such as small molecules,2−5 quantum dots,6−8 conjugated π systems,9−16 bulk-heterojunctions,17,18 and metal complexes.19 These systems exhibit interesting chemical and physical properties, including singlet fission,12,14,20,21 electron transfer,22,23 proton transfer,4,24 and other excited state reactions that will prove useful in solving the future energy crisis. However, NA molecular dynamics (NAMD) methods, such as Ehrenfest dynamics25−28 and surface hopping,1,29 can be very computationally demanding, especially with ab initio description of electronic structure. Although full ab initio NAMD calculations have been performed on solution phase systems in the past,30 approximations such as frozen density embedding31 prove critical when large systems or systems with significant environmental contributions are simulated.2,32 Current solutions to this problem include parametrizing classical or semiclassical force fields to ab initio calculations,33−35 or performing QM/MM calculations with the environment being simulated classically.36 A more direct solution, developed in the current work, employs the FMO molecular dynamics (FMOMD) method37−41 to simplify the computational complexity of ab initio NAMD. The FMO technique is used to compute both the MD trajectory and the NA coupling constants required for NAMD calculations. The FMO methodology was developed in 1999 by Kitaura and co-workers to treat polymers, molecules in solution, and other systems that can be easily divided into fragments.42 Since © 2016 American Chemical Society

then, the method has been used to successfully study a variety of systems, ranging from proteins,43 such as ubiquitin44 and human estrogen receptor,45 to condensed phase reactions.39,46 Both the time saving nature of fragmentation due to the nonlinear scaling of quantum simulations and the ease of parallelization47 have put the FMO method in a good position to handle the quantum chemistry problems of the future. In the FMO methodology, a large system is initially divided into fragments consisting of a small number of heavy atoms. After a Hückel guess is made on each fragment, ab initio calculations are performed on individual fragments in the presence of the electric field from all other fragments. The process is repeated until the electron densities of each fragment are self-consistent with the densities on other fragments. This results in an efficient way to obtain the gradient for MD. Here, the FMO methodology is extended to computation of NA coupling needed for NAMD. To illustrate and test the developed FMO-NAMD approach, we consider excited state dynamics of the Fe(CO)4 complex in the gas phase and in solution. Gas phase photodegradation of Fe(CO)5 to form Fe(CO)4 has been extensively studied over the last few decades to gain a better understanding of photodegradation of metal complexes.48−53 Briefly, Fe(CO)5 is excited by a 267 nm photon to the first singlet excited state (S1).50 Then, in less than 100 fs, a CO ligand is eliminated, while the molecule still remains in the excited state.50 Following the elimination of a ligand, the complex decays into the much more stable S0 state, though this state is still higher than the Received: June 3, 2016 Revised: July 21, 2016 Published: August 18, 2016 7205

DOI: 10.1021/acs.jpca.6b05607 J. Phys. Chem. A 2016, 120, 7205−7212

Article

The Journal of Physical Chemistry A triplet ground state.50,54 The lifetime of the S1 state is dependent on environmental effects. Recent advances in Xray laser femtosecond spectroscopy54,55 has allowed further characterization of the photodegradation of Fe(CO)5 in ethanol. Important differences between the gas phase and solution phase processes have been observed. In particular, it is observed that the lifetime of the S1 state of Fe(CO)4 in solution is increased by a factor of 4 over that of the same molecule in the gas phase, namely, from 50 to 200 fs. Further, the ordering between the decay of the S1 state and the attachment of the ethanol molecule to the iron core is unknown.54 The NAMD FMO interface is used to investigate this process.

SijIJ = ⟨φi I|φjJ⟩ S̅ ijIJ = ⟨φ(t )iI |φ(t −1) jJ ⟩

IJ E = ∑ E + ∑ Ẽ + I

IJ

Δt

|φIi ⟩

Here is a single basis function located on monomer I, t represents the current time step, and (t − 1) represents the previous time step. The derivative overlap matrices as well as the Fock matrices are essential in computing the nonadiabatic couplings for the FSSH calculations as discussed below. 2.2. Calculation of the Nonadiabatic Couplings. To perform nonadiabatic coupling calculations on the full system, the dimer Fock and overlap matrices are combined to form an approximate Fock and overlap matrix for the whole system. Because the overlap matrices will not depend on the order of terms kept in the FMO calculation, they can be obtained exactly from the dimer calculation by combining the dimer matrices in block fashion.

2. THEORY 2.1. Fragment Molecular Orbital Method. In this section the fragment molecular orbital method will be briefly discussed. The underlying principal of FMO is the many-body expansion, which is a prescription for computing properties of a large system using only calculations performed on small subsystems, or fragments. For instance, the energy manybody expansion is56 I

⟨φ(t )iI |φ(t ) jJ ⟩ − ⟨φ(t )iI |φ(t −1) jJ ⟩

S′ijIJ =

∑ ẼIJK + ··· + Ẽ IJK

IJ

E ̃ = EIJ − EI − E J Here, Ẽ IJ is the interaction energy between two monomers and can be computed by computing the dimer energy and subtracting the two monomer energies. Similar equations exist for the three and higher body energy terms. One of the main features of the many-body expansion is that if it is not truncated it will simplify to the exact energy, thus offering a clear path to improve the quality of calculations. Additionally, this many-body expansion holds for many properties of the system, facilitating the calculation of things like the energy gradient for molecular dynamics calculations. A brief summary of running an FMO calculation will be given here with a more detailed description to be found in the review by Fedorov and Kitaura.46 The calculation begins by dividing the system into fragments to use as pieces of the many-body expansion. Ab-initio calculations are then performed on each of these fragments, which generate an initial electron density on each fragment. With these densities, single electron electrostatic interactions are added to each monomer Hamiltonian on the basis of the effects of the other nearby fragments and the ab initio calculation is then repeated on each monomer including these effects. These monomer calculations are then repeated until the interfragment electrostatic interactions are selfconsistent. After the monomer calculations are completed, higher order cluster calculations are performed. Usually these higher order clusters include either just dimer calculations (FMO2) or dimer and trimer calculations (FMO3). This paper remains restricted to FMO2 for computational expediency. To interface the FMO method to FSSH, we will need the following three quantities computed from the FMO-MD run, namely, the Fock matrix for each fragment and fragment pair (FI and FIJ), the atomic orbital overlap matrix for each fragment pair (SIJ), and the overlap between atomic orbitals on subsequent time steps (S̅IJ). The first two of these quantities are normally computed during a FMO-MD calculation, whereas the code to compute the third quantity was added to an in house copy of Gamess. The two overlap matrices are used to compute the derivative overlap matrix (S′IJij ) as follows:

⎡ S(t )1 ⋯ S(t )1N ⎤ ⎢ ⎥ ⋱ ⋮ ⎥ S(t ) = ⎢ ⋮ ⎢ ⎥ ⎢⎣ S(t )N1 ⋯ S(t )N ⎥⎦

(4)

⎡ S′(t )1 ⋯ S′(t )1N ⎤ ⎥ ⎢ ⋱ ⋮ ⎥ S′(t ) = ⎢ ⋮ ⎥ ⎢ ⎢⎣ S′(t )N1 ⋯ S′(t )N ⎥⎦

(5)

In the event of the FMO calculation not performing a dimer calculation between to fragments because it is deemed to be too insignificant of a contribution, the overlap matrix between these two monomers is assumed to be zero. To obtain the Fock matrix for the full system, we apply the many-body expansion in much the same way that other quantities are computed in the FMO method. F=

∑ FI + ∑ F ̃IJ + ∑ F ̃IJK + ... + F I

IJ

(6)

IJK

̃IJ

I

where F is the Fock matrix for monomer I, F is the two-body energy for dimer IJ, etc. The many-body expansion of the Fock matrix is truncated at the order of the underlying FMO simulation, which is 2 in this paper:

F≈

∑ FI + ∑ F ̃IJ I

IJ

(7)

where the two-body contribution to the Fock matrix can be written IJ

F ̃ = F IJ − F I − F J

(8)

Once the approximate Fock matrix is obtained, it is transformed into an orthonormal basis using the full system overlap matrix and diagonalized, in a method analogous to a single iteration of the Hartree−Fock procedure. The transformed Fock matrix will be denoted F̅ to avoid confusion with the time derivatives: F̅ = (S−1/2)T F(S−1/2)

(9)

Diagonalizaing this matrix results in the transformed orbital coefficients C̅ and the orbital energies ε, where 7206


Article

The Journal of Physical Chemistry A C = S−1/2C̅

(10)

Finally, the full system molecular orbitals can be expressed in terms of these obtained coefficients: |ψi ⟩ =

F2 ⋮ 0

⋯ 0⎤ ⎥ ⋯ 0⎥ ⎥ ⋱ 0 ⎥ ⎥ 0 FN ⎦

dij =

Ψi

d Ψj dt

ψp

i ,m

m

i ,m

d ψ dt pj ,m (15)

d ψ dt pj ,m

∑ Cp

dij =

(16)

i ,m

, kφk

k

(12)

d dt

∑ Cp

j,m

, lφl

(17)

l

Applying the chain rule and recalling that both C and φ are time-dependent quantities can simplify this expression to one dependent on only the orbital coefficients and previously computed overlap matrices. (C p l(t ) − C p l(t −1))

M

dij =

∑ Cp

j,m

j,m

k i ,m

Δt

∑ Cp

kC pj , m lS′kl (t )

l ,k

Skl(t )

M

+

l ,k

i ,m

(18)

Finally, these nonadiabatic couplings and the orbital energies obtained from the diagonalization of the approximate Fock matrix are then fed into the standard surface hopping equations, discussed below. 2.3. Overview of Nonadiabatic Molecular Dynamics. NAMD dynamics is a widely used method,8,9,57−63 and fewest switches SH (FSSH) is the most popular NAMD technique. Introduced by Tully in 1990,1 FSSH minimizes the number of hops, which is a desirable feature in condensed phase simulations, and satisfies several important criteria, including trajectory branching, internal consistency and detailed balance between transitions upward and downward in energy.58 Later, FSSH was advanced to include decoherence effects,64 to account for the superexchange mechanism of population transfer,63,65 and to resolve the numerical problem with FSSH probabilities in cases of small NA coupling.66 Due to the short time scales of the excited state decay in the current system, the original FSSH algorithm is employed. The NAMD methodology is rooted in the time-dependent Schrödinger equation:

(13)

pi,j is a list determining which orbitals are occupied in the i’th state, and Â is the antisymmetric operator. Nonadiabatic coupling matrix elements are computed between these states in the same way they are computed for other surface hopping methods:

dij =

ψp

From here the expansion of the full system orbitals into the basis functions can be applied to compute the value of di,j:

̂ |Ψ⟩ i = A ∏ |ψp ⟩ i ,j

∑

Because orbitals are orthogonal, di,j will only be nonzero if the two states differ by only a single orbital. Define pi,m and pj,m as the indices of the single different orbital for states i and j.

This Fock matrix represents uncoupled monomers and, because it is block diagonal, will have eigenvectors that are the monomer orbitals assuming the overlaps between basis functions on different monomers are small. This simple result effectively represents the zeroth-order case in the many-body expansion and will prevent almost all energy and electron transfer between monomers. The other limit is where the Fock matrix is obtained from an untruncated many-body expansion, in which case the Fock matrix is identical to that obtained from a full system ab initio calculation. Because the overlap matrix, S, is also identical to that of the full system, the eigenvectors of this matrix are the orbitals obtained from the full system ab initio calculation. When the FSSH method is applied to these orbitals, the result will be identical to that obtained using a FSSH calculation on the whole system without FMO. Thus, this method maintains the core principal of FMO that untruncated FMO calculation will be identical to a full system calculation performed without FMO. Once the orbitals |ψi⟩ are obtained, they can be used to construct Slater determinants for the ground and excited states, where the excited states are expressed as single Slater determinant wave functions with one or more electrons excited from occupied orbitals to unoccupied orbitals.

j

l

=

i ,l

i ,l j ,l

l≠m

This method of obtaining the full system orbitals from the fragment and cluster simulations has two important limits. The first is the limit of using only the one-body terms, where the Fock matrix becomes block diagonal with the monomer Fock matrices making up the diagonal blocks: 0

∏ ψp

∏ δ p ,p

(11)

⎡ F1 ⎢ ⎢0 I ∑F = ⎢ ⎢ ⋮ I ⎢ ⎣ 0

i ,k

k

∑ Cij|φi⟩ j

d dt

∏ ψp

dij =

i∇

∂Ψ(r ,R ,t ) = H(r ,R ,t ) Ψ(r ,R ,t ) ∂t

(19)

where Ψ is the nuclear/electronic wave function, H is the system Hamiltonian, r represents the electronic coordinantes, and R is the nuclear coordinates. The wave function can be expressed as a sum over (adiabatic) electronic states multiplied by a weighting function.

(14)

where {ψj} are the set of all basis functions for the system. Applying these two equations together, we get the following quantity for the nonadiabatic coupling. It is important to note that due to the orthogonality of the orbitals, the antisymmetrizing operator has no effect on this derivation and can be ignored.

Ψ(r ,R ,t ) =

∑ χi (t ,R(t )) Φi(r ,R(t )) j

7207

(20) DOI: 10.1021/acs.jpca.6b05607 J. Phys. Chem. A 2016, 120, 7205−7212

Article

The Journal of Physical Chemistry A Here, χi(t,R(t)) is a time-dependent coefficient determining the contribution of a given electronic state to the evolving wave function. The Hamiltonian for the system is the sum of nuclear and electronic kinetic energies, as well as the electronic potential energy determined from the interactions between the nuclei and electrons. H = Tnucl + Tel + Vel

iron carbon bond from 1.71 to 1.81 Å, though the carbon oxygen bond remained the same at 1.17 Å. Because solvent substrate interactions have been successfully described in the past with ground state calculations,78,79 the geometric factors were deemed to be small enough to use the classical path approximation to FSSH.67 Dynamics were performed on Fe(CO)4 in the gas phase and Fe(CO)4 surrounded by up to 12 ethanol molecules to determine the effect of solvation on the complex NA dynamics. Ten different initial configuration of atoms for Fe(CO)4 surrounded with 12 ethanol molecules were generated with the Packmol program.80 Dynamics were performed using velocity Verlet integrator with a Nose/Hoover thermostat set to 300 K. The initial velocities for the atoms were chosen randomly from a Boltzmann distribution. The solvated clusters were held together using a harmonic potential with a force constant of 0.75 kcal/(mol Å2) that began 6.5 Å from the center of the box. Similar procedures for containing a box of molecules have been used in other applications.39,81 Recent work on FMO with periodic boundary conditions may allow a more accurate treatment of the boundary in the future.40 2.5. Timing. Timings for FMO calculations have been well reported on in the past.41,47,82 A brief overview of the timings associated with the current NAMD calculations is presented in Table 1. Simulations on Fe(CO)4 with 12 ethanol molecules

(21)

This form for the Hamiltonian can be substituted into the time-dependent Schrödinger equation, along with the assumed form for each electronic state, to produce a series of differential equations governing the evolution of the expansion coefficient i∇

∂χi (t ,R(t ))

= ∂t ⎛ ⎞ ⎯⇀ ⎯ (1) (2) ⎯⇀ ⎜ ⎟ d d ij ij ∇ − ∇2 ∑ ⎜Tnucl + Ei(t ,R) − ∇2 ⎟χ (t ,R(t )) M 2M ⎟ i j ⎜ ⎝ ⎠ (22)

This equation is simplified further by combining the diagonal part of the second-order derivative term with electronic energy and ignoring the off-diagonal part. In a fully quantum description of NA dynamics, χi(t,R(t)) represents a nuclear wave packet associated with the ith electronic state. A quantum-classical approximation employs a classical description of nuclear motion, turning χi(t,R(t)) into an expansion coefficient, ci(t), Ψ(r ,R ,t ) =

∑ ci(t ) φi(r ,R(t )) i

Table 1. Comparison between the Time (s) Required To Perform a Single Time Step on the Given System Both Using the FMO Method and Running the Whole System in a Single Calculationa

(23) using FMO without FMO

Equation 4 transforms into i∇

dci(t ) = ∂t

12 ethanol molecules (3 nodes)

64 ethanol molecules (8 nodes)

809 975

960

a

A calculation involving 64 ethanol molecules could not be performed without FMO due to memory limitations.

∑ (εiδij + i∇dij)cj j

(24)

were run both with the FMO method and as a large single ab initio calculation. FMO offered a 17% speedup in the calculation, which is promising as the system is small by FMO standards. As more ethanol molecules are added, calculations run without FMO will quickly become infeasible. Indeed, dynamics were also performed on a system containing a single Fe(CO)4 and 64 ethanol molecules. Though data sufficient for statistical sampling needed to converge NAMD calculations were not collected, due to the computational expense and success of the smaller simulations, the available results do illustrate how FMO can be used to efficiently simulate larger systems. Even though the larger system contains 4.5 times more electrons, it only required 3 times as much CPU time. Note the difference in the number of nodes used for the 12 and 64 ethanol simulations, Table 1. This better than linear scaling is due to the dynamic load balancing implemented in the FMO method, which works better with many small jobs to distribute between the nodes. A calculation with 64 ethanol molecules could not be performed as a single piece for comparison with the FMO calculation due to memory constraints.

The only pieces of information that must be obtained from quantum calculations, then, are the state energies, εi, and the nonadiabatic coupling matrix elements, dij, the calculation of which has been discussed above. Once the MD trajectory and the electronic data are gathered using GAMESS, 1000 stochastic FSSH calculations are performed using the PYXAID computational suite of programs.67,68 PYXAID was released under the GNU license on the basis of the NAMD methodologies developed in the Prezhdo group and applied to a variety of nanoscale systems.69−71 2.4. Computational Details. Density functional theory (DFT) and time-dependent DFT (TDDFT) calculations were performed in Gamess72 using the PBE073 functional and the 631G basis74,75 for all non-iron atoms. For the iron atom, the model core potential (MCP) method at the TZP level was used to reduce the computational complexity by removing core electrons.76 Model core potentials are designed to obtain the correct correlation and nodal structure of valence orbitals.77 These core potentials have previously been interfaced with the FMO method to compute polypeptide fragments in solution.37 Fe(CO)4 with no environmental effects was optimized in both the ground and first excited state to determine the effect of electronic excitation on the geometry of the molecule. This revealed that excitation resulted in a slight elongation in the

3. RESULTS 3.1. Validation. To validate our underlying DFT calculations, a simple single-point comparison between DFT 7208


Article

The Journal of Physical Chemistry A using the 6-31G/MCP basis as used in the dynamics simulations, DFT with a CC-PVTZ83 basis, and EOMCCSD84 using the 6-31G/MCP basis (Table 2). This analysis Table 2. Table of Transition Energy, Dipole Moment, and Orbital Transition for the First Excited State of Fe(CO)4 Obtained Using Various Methodsa

method

excitation energy (eV)

ground state transition dipole moment magnitude

orbital transition

DFT, 6-31G/MCP DFT, CC-PVTZ EOM-CCSD, 6-31G/MCP

1.577 1.522 1.679

0.26 0.25 0.46

35, 36 → 37 41 → 42 36 → 37

Figure 2. Isolated Fe(CO)4 with no additional ethanol molecules. The computed lifetime of the S1 state is 50.8 fs, in excellent agreement with the experimental value of 47 fs.

T0 state. This is a good indication that the simulations are working as expected. Following the isolated gas phase simulations, Fe(CO)4 bound to a single ethanol molecule was investigated. The nonradiative decay dynamics obtained from these simulations are reported in Figure 3. The lifetime of the first excited state in

a

These data are provided to verify the accuracy of DFT used with the MCP basis set as applied throughout the remainder of the paper.

indicated that the orbital ordering as well as transition strength and excitation energy were agreed upon by all methods. These calculations were again performed in GAMESS. The oscillator strength is provided merely as a confirmation that the methods are returning the same states. It is not important that the ground to first excited state transition is dark because the state is being populated by a transition from the Fe(CO)5 excited state. 3.2. Gas Phase Dynamics. Initially, an isolated Fe(CO)4 was simulated in the gas phase so that the de-excitation rate could be compared to experiment and the solution phase simulations. Figure 1 shows the key orbitals of Fe(CO)4 as

Figure 3. Fe(CO)4 bound to a single ethanol molecule. The computed lifetime of the S1 state is 65.5 fs, which is longer than the gas phase result, Figure 2.

the presence of a single ethanol molecule grows to 65.5 fs from 50.8 fs in the gas phase. The increase in the lifetime is a positive reflection of what happens in solution, where experimentally the lifetime is increased to 200 fs, though the full effect of the condensed phase environment is not seen yet. 3.3. Condensed Phase Dynamics. To model the Fe(CO)4 complex in solution, dynamics were carried out on Fe(CO)4 in a box of 12 ethanol molecules. Twelve ethanol molecules were used to simulate the first layer of solvation, while also keeping the simulation relatively affordable. Condensed phase dynamics involves binding of a solvent molecule to the complex. Representative MD frames showing the system before and after the binding event are shown in Figure 4. The initial conditions of the system were such that the ethanol molecule was not bound to the iron center at the outset of the calculation but became bound at some subsequent time, usually around 400−600 fs after the start of the simulation. This allowed us to compute the nonradiative decay rate of the first excited state both before binding and after binding. The 200 fs for each run before the binding event was averaged together to produce Figure 5. The data indicate that, prior to ethanol binding to the iron atom, the lifetime of the excited state is near 100 fs. Although this is not in perfect agreement with the experimental results54 obtained with femtosecond X-ray spectroscopy and indicating a 200 fs lifetime in solution, it does demonstrate a significant increase in lifetime from the 50 fs in the gas phase calculations. The trend demonstrated with Figures 2−4 leads us to expect that adding more solvent molecules will increase the simulated time-

Figure 1. Highest four occupied and lowest four unoccupied gas phase orbitals of Fe(CO)4. The orbitals were computed using DFT and the combined 6-31G/MCP basis.

computed by DFT. The first excited state of the complex was assumed to be a simple HOMO to LUMO transition, which agrees well with previous characterization by Wernet et al. that the S1 state is π → σ* transition.54 The results of the NAMD simulations on the gas phase complex are reported in Figure 2. A decay rate of 50.8 fs for the S1 → S0 transitions was obtained, showing excellent agreement with the gas phase experimental value of 47 fs obtained by Schmid and co-workers through femtosecond pump probe spectroscopy.50 Although the S0 state is not the ground state of FeCO4, the decay to the lowest singlet state is the primary decay channel of the S1 state due to spin change to the ground 7209


Article

The Journal of Physical Chemistry A

Figure 6. Fe(CO)4 surrounded by ethanol molecules after one binds to the iron atom. The computed lifetime of the S1 state is 11.4 fs. This indicates that once an ethanol molecule binds to the iron center, the decay of the excited state will happen soon after.

4. CONCLUSIONS Here, the interface between an interface between the FMO method in GAMESS and the PYXAID program for computing NAMD dynamics in condensed phase and nanoscale systems has been detailed. The efficient treatment of many-body interaction with the FMO approximation provides great computational savings in systems that can be easily separated into weakly interacting fragments. The computational cost of the additional calculations needed for the interface was determined to be negligible, because the NA coupling is a one-electron operator that is straightforward to calculate. The nonradiative decay processes in the Fe(CO)4 complex in the gas phase and in an ethanol solution is used as a test of the FMO-NAMD method. The calculations have shown excellent agreement with the time-resolved optical and X-ray spectroscopies. The decay rates obtained using our FMO-NAMD method give us insight into the order of CO dissociation and ethanol binding. Because the lifetime of the excited state is so short once the ethanol molecule binds to the iron center in solution, we can say with confidence that the CO dissociation from excited Fe(CO)5 and ethanol attachment is not a concerted process. Our calculations suggest that a CO ligand dissociates from Fe(CO)5 and forms Fe(CO)4 in the excited state. This species then exists in solution for around 100−200 fs before rebinding to an ethanol molecule. What still remains unclear is whether binding of the ethanol causes the relaxation of the S1 state, or whether the relaxation happens before an interaction with the ethanol molecule. Both mechanisms are consistent with the current simulation, and more thorough investigation is required. The advantage of the developed method resides in both a reduction in computational complexity and the ability to model interactions between localized excitations. Future developments of the FMO-NAMD approach will include modeling of energy transfer between fragments using data obtained from the dimer calculations performed by the FMO method, the inclusion of spin orbit coupling values from monomer excited state calculations to compute intersystem crossing rates, and implementation and testing of more advanced NAMD techniques, beyond the most popular FSSH approach.

Figure 4. Frames of Fe(CO)4 in a box of ethanol molecules. (A) is a frame pulled before the binding event occurs, and (B) is after the binding event.

Figure 5. Fe(CO)4 surrounded by ethanol molecules before any binding to the iron center. The computed lifetime of the S1 state is 100 fs, representing a significant increase from the gas phase value, which is consistent with experimental data.

scale further, bringing the NAMD results in close agreement with experiment. After the binding event, the lifetime of the excited state reduces to 11.4 fs, as seen in Figure 6. This is very telling and gives hints as to whether the CO detachment and ethanol attachment are concerted or sequential. The order of magnitude acceleration of the nonradiative relaxation dynamics indicates that the binding of the ethanol group puts a hard upper limit on the lifetime of the excited state.

■ ■

AUTHOR INFORMATION

Notes

The authors declare no competing financial interest.

ACKNOWLEDGMENTS This work was supported as part of the Computational Materials Sciences Program funded by the U.S. Department 7210


Article

The Journal of Physical Chemistry A

(35) Kazaryan, A.; Lan, Z.; Schäfer, L. V.; Thiel, W.; Filatov, M. J. Chem. Theory Comput. 2011, 7, 2189. (36) Mendieta-Moreno, J. I.; Walker, R. C.; Lewis, J. P.; GómezPuertas, P.; Mendieta, J.; Ortega, J. J. Chem. Theory Comput. 2014, 10, 2185. (37) Nagata, T.; Fedorov, D. G.; Kitaura, K. Theor. Chem. Acc. 2012, 131, 1136. (38) Nakata, H.; Schmidt, M. W.; Fedorov, D. G.; Kitaura, K.; Nakamura, S.; Gordon, M. S. J. Phys. Chem. A 2014, 118, 9762. (39) Komeiji, Y.; Ishikawa, T.; Mochizuki, Y.; Yamataka, H.; Nakano, T. J. Comput. Chem. 2009, 30, 40. (40) Fujita, T.; Nakano, T.; Tanaka, S. Chem. Phys. Lett. 2011, 506, 112. (41) Nagata, T.; Brorsen, K.; Fedorov, D. G.; Kitaura, K.; Gordon, M. S. J. Chem. Phys. 2011, 134, 124115. (42) Kitaura, K.; Ikeo, E.; Asada, T.; Nakano, T.; Uebayasi, M. Chem. Phys. Lett. 1999, 313, 701. (43) Mochizuki, Y.; Yamashita, K.; Fukuzawa, K.; Takematsu, K.; Watanabe, H.; Taguchi, N.; Okiyama, Y.; Tsuboi, M.; Nakano, T.; Tanaka, S. Chem. Phys. Lett. 2010, 493, 346. (44) Komeiji, Y.; Ishida, T.; Fedorov, D. G.; Kitaura, K. J. Comput. Chem. 2007, 28, 1750. (45) Fukuzawa, K.; Kitaura, K.; Uebayasi, M.; Nakata, K.; Kaminuma, T.; Nakano, T. J. Comput. Chem. 2005, 26, 1. (46) Fedorov, D. G.; Kitaura, K. J. Phys. Chem. A 2007, 111, 6904. (47) Fletcher, G. D.; Fedorov, D. G.; Pruitt, S. R.; Windus, T. L.; Gordon, M. S. J. Chem. Theory Comput. 2012, 8, 75. (48) Karny, Z.; Naaman, R.; Zare, R. N. Chem. Phys. Lett. 1978, 59, 33. (49) Ryther, R. J.; Weitz, E. J. Phys. Chem. 1992, 96, 2561. (50) Trushin, S. A.; Fuss, W.; Kompa, K. L.; Schmid, W. E. J. Phys. Chem. A 2000, 104, 1997. (51) Lee, K.; Yoo, H. S.; Ku, J. K. Chem. Phys. Lett. 1996, 262, 610. (52) Ryther, R. J.; Weitz, E. J. Phys. Chem. 1991, 95, 9841. (53) Bañares, L.; Baumert, T.; Bergt, M.; Kiefer, B.; Gerber, G. Chem. Phys. Lett. 1997, 267, 141. (54) Wernet, P.; Kunnus, K.; Josefsson, I.; Rajkovic, I.; Quevedo, W.; Beye, M.; Schreck, S.; Grübel, S.; Scholz, M.; Nordlund, D.; Zhang, W.; Hartsock, R. W.; Schlotter, W. F.; Turner, J. J.; Kennedy, B.; Hennies, F.; de Groot, F. M. F.; Gaffney, K. J.; Techert, S.; Odelius, M.; Föhlisch, A. Nature 2015, 520, 78. (55) Huse, N.; Cho, H.; Hong, K.; Jamula, L.; de Groot, F. M. F.; Kim, T. K.; McCusker, J. K.; Schoenlein, R. W. J. Phys. Chem. Lett. 2011, 2, 880. (56) Fedorov, D. G.; Kitaura, K. J. Chem. Phys. 2004, 120, 6832. (57) Werner, U.; Mitrić, R.; Suzuki, T.; Bonačić-Koutecký, V. Chem. Phys. 2008, 349, 319. (58) Parandekar, P. V.; Tully, J. C. J. Chem. Phys. 2005, 122, 094102. (59) Schmidt, J. R.; Parandekar, P. V.; Tully, J. C. J. Chem. Phys. 2008, 129, 044104. (60) Subotnik, J. E.; Shenvi, N. J. Chem. Phys. 2011, 134, 024105. (61) Subotnik, J. E.; Vura-Weis, J.; Sodt, A. J.; Ratner, M. A. J. Phys. Chem. A 2010, 114, 8665. (62) Wang, L.; Trivedi, D.; Prezhdo, O. V. J. Chem. Theory Comput. 2014, 10, 3598. (63) Sifain, A. E.; Wang, L.; Prezhdo, O. V. J. Chem. Phys. 2015, 142, 224102. (64) Bittner, E. R.; Rossky, P. J. J. Chem. Phys. 1995, 103, 8130. (65) Akimov, A. V.; Prezhdo, O. V. Phys. Rev. Lett. 2014, 113, 153003. (66) Wang, L.; Prezhdo, O. V. J. Phys. Chem. Lett. 2014, 5, 713. (67) Akimov, A. V.; Prezhdo, O. V. J. Chem. Theory Comput. 2013, 9, 4959. (68) Akimov, A. V.; Prezhdo, O. V. J. Chem. Theory Comput. 2014, 10, 789. (69) Neukirch, A. J.; Hyeon-Deuk, K.; Prezhdo, O. V. Coord. Chem. Rev. 2014, 263−264, 161. (70) Trivedi, D. J.; Wang, L.; Prezhdo, O. V. Nano Lett. 2015, 15, 2086.

of Energy, Office of Science, Basic Energy Sciences, under Award Number DE-SC00014607.

■

REFERENCES

(1) Tully, J. C. J. Chem. Phys. 1990, 93, 1061. (2) Gámez, J. A.; Weingart, O.; Koslowski, A.; Thiel, W. J. Chem. Theory Comput. 2012, 8, 2352. (3) Cantatore, V.; Granucci, G.; Persico, M. Comput. Theor. Chem. 2014, 1040−1041, 126. (4) Goyal, P.; Schwerdtfeger, C. A.; Soudackov, A. V.; HammesSchiffer, S. J. Phys. Chem. B 2015, 119, 2758. (5) Schwerdtfeger, C. A.; Soudackov, A. V.; Hammes-Schiffer, S. J. Chem. Phys. 2014, 140, 034113. (6) Liu, J.; Neukirch, A. J.; Prezhdo, O. V. J. Phys. Chem. C 2014, 118, 20702. (7) Hyeon-Deuk, K.; Prezhdo, O. V. J. Phys.: Condens. Matter 2012, 24, 363201. (8) Fischer, S. A.; Lingerfelt, D. B.; May, J. W.; Li, X. Phys. Chem. Chem. Phys. 2014, 16, 17507. (9) Jaeger, H. M.; Fischer, S.; Prezhdo, O. V. J. Chem. Phys. 2012, 137, 22A545. (10) Habenicht, B. F.; Craig, C. F.; Prezhdo, O. V. Phys. Rev. Lett. 2006, 96, xxx. (11) Habenicht, B. F.; Prezhdo, O. V. Phys. Rev. Lett. 2008, 100, 197402. (12) Wang, L.; Olivier, Y.; Prezhdo, O. V.; Beljonne, D. J. Phys. Chem. Lett. 2014, 5, 3345. (13) Liu, J.; Li, B.-W.; Tan, Y.-Z.; Giannakopoulos, A.; SanchezSanchez, C.; Beljonne, D.; Ruffieux, P.; Fasel, R.; Feng, X.; Müllen, K. J. Am. Chem. Soc. 2015, 137, 6097. (14) Tamura, H.; Huix-Rotllant, M.; Burghardt, I.; Olivier, Y.; Beljonne, D. Phys. Rev. Lett. 2015, 115, 107401. (15) Aragó, J.; Troisi, A. J. Chem. Phys. 2015, 142, 164107. (16) Knox, J. E.; Halls, M. D.; Schlegel, H. B. J. Comput. Theor. Nanosci. 2006, 3, 398. (17) Oldani, N.; Tretiak, S.; Bazan, G.; Fernandez-Alberti, S. Energy Environ. Sci. 2014, 7, 1175. (18) Jailaubekov, A. E.; Willard, A. P.; Tritsch, J. R.; Chan, W.-L.; Sai, N.; Gearba, R.; Kaake, L. G.; Williams, K. J.; Leung, K.; Rossky, P. J.; Zhu, X.-Y. Nat. Mater. 2012, 12, 66. (19) Tavernelli, I.; Curchod, B. F. E.; Rothlisberger, U. Chem. Phys. 2011, 391, 101. (20) Johnson, J. C.; Nozik, A. J.; Michl, J. Acc. Chem. Res. 2013, 46, 1290. (21) Smith, M. B.; Michl, J. Chem. Rev. 2010, 110, 6891. (22) Fornari, R. P.; Aragó, J.; Troisi, A. J. Chem. Phys. 2015, 142, 184105. (23) Lee, M. H.; Aragó, J.; Troisi, A. J. Phys. Chem. C 2015, 119, 14989. (24) Soudackov, A. V.; Hammes-Schiffer, S. J. Phys. Chem. Lett. 2014, 5, 3274. (25) Subotnik, J. E. J. Chem. Phys. 2010, 132, 134112. (26) Fischer, S. A.; Chapman, C. T.; Li, X. J. Chem. Phys. 2011, 135, 144102. (27) Li, X.; Tully, J. C.; Schlegel, H. B.; Frisch, M. J. J. Chem. Phys. 2005, 123, 084106. (28) Ding, F.; Goings, J. J.; Liu, H.; Lingerfelt, D. B.; Li, X. J. Chem. Phys. 2015, 143, 114105. (29) Fingerhut, B. P.; Dorfman, K. E.; Mukamel, S. J. Chem. Theory Comput. 2014, 10, 1172. (30) Vogel, D. J.; Kilin, D. S. Mol. Phys. 2015, 113, 397. (31) Ramos, P.; Pavanello, M. J. Chem. Theory Comput. 2014, 10, 2546. (32) Silva-Junior, M. R.; Thiel, W. J. Chem. Theory Comput. 2010, 6, 1546. (33) Zhu, C.; Jasper, A. W.; Truhlar, D. G. J. Chem. Theory Comput. 2005, 1, 527. (34) Kazaryan, A.; Kistemaker, J. C. M.; Schäfer, L. V.; Browne, W. R.; Feringa, B. L.; Filatov, M. J. Phys. Chem. A 2010, 114, 5058. 7211


Article

The Journal of Physical Chemistry A (71) Postupna, O.; Long, R.; Prezhdo, O. V. J. Phys. Chem. C 2015, 119, 12088. (72) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; Windus, T. L.; Dupuis, M.; Montgomery, J. A. J. Comput. Chem. 1993, 14, 1347. (73) Adamo, C.; Barone, V. J. Chem. Phys. 1999, 110, 6158. (74) Ditchfield, R.; Hehre, W. J.; Pople, J. A. J. Chem. Phys. 1971, 54, 724. (75) Hehre, W. J.; Ditchfield, R.; Pople, J. A. J. Chem. Phys. 1972, 56, 2257. (76) Ishikawa, T.; Mochizuki, Y.; Nakano, T.; Amari, S.; Mori, H.; Honda, H.; Fujita, T.; Tokiwa, H.; Tanaka, S.; Komeiji, Y.; Fukuzawa, K.; Tanaka, K.; Miyoshi, E. Chem. Phys. Lett. 2006, 427, 159. (77) Miyoshi, E.; Mori, H.; Hirayama, R.; Osanai, Y.; Noro, T.; Honda, H.; Klobukowski, M. J. Chem. Phys. 2005, 122, 074104. (78) Liu, L.-M.; Car, R.; Selloni, A.; Dabbs, D. M.; Aksay, I. A.; Yetter, R. A. J. Am. Chem. Soc. 2012, 134, 19011. (79) Wu, H.-Z.; Liu, L.-M.; Zhao, S.-J. Phys. Chem. Chem. Phys. 2014, 16, 3299. (80) Martínez, L.; Andrade, R.; Birgin, E. G.; Martínez, J. M. J. Comput. Chem. 2009, 30, 2157. (81) Fujiwara, T.; Mochizuki, Y.; Komeiji, Y.; Okiyama, Y.; Mori, H.; Nakano, T.; Miyoshi, E. Chem. Phys. Lett. 2010, 490, 41. (82) Fedorov, D. G.; Kitaura, K. Chem. Phys. Lett. 2006, 433, 182. (83) Dunning, J. T. H. J. Chem. Phys. 1989, 90, 1007. (84) Piecuch, P.; Kucharski, S. A.; Kowalski, K.; Musiał, M. Comput. Phys. Commun. 2002, 149, 71.

7212


Fragment Molecular Orbital Nonadiabatic Molecular Dynamics for

Recommend Documents