Nonadiabatic Dynamics of Photoinduced Proton-Coupled Electron

Dec 29, 2014 - Puja Goyal, Christine A. Schwerdtfeger, Alexander V. Soudackov, and Sharon Hammes-Schiffer*. Department of Chemistry, 600 South ...
0 downloads 0 Views 2MB Size
Article pubs.acs.org/JPCB

Nonadiabatic Dynamics of Photoinduced Proton-Coupled Electron Transfer in a Solvated Phenol−Amine Complex Puja Goyal, Christine A. Schwerdtfeger, Alexander V. Soudackov, and Sharon Hammes-Schiffer* Department of Chemistry, 600 South Mathews Avenue, University of Illinois at Urbana-Champaign, Urbana, Illinois 61801, United States S Supporting Information *

ABSTRACT: Photoinduced concerted electron-proton transfer (EPT), denoted photo-EPT, is important for a wide range of energy conversion processes. Transient absorption and Raman spectroscopy experiments on the hydrogen-bonded p-nitrophenylphenol−t-butylamine complex, solvated in 1,2-dichloroethane, suggested that this complex may undergo photoEPT. The experiments probed two excited electronic states that were interpreted as an intramolecular charge transfer (ICT) state and an EPT state. Herein mixed quantum mechanical/molecular mechanical nonadiabatic surface hopping dynamics is used to investigate the relaxation pathways following photoexcitation. The potential energy surface is generated on the fly with a semiempirical floating occupation molecular orbital complete active space configuration interaction method for the solute molecule and a molecular mechanical force field for the explicit solvent molecules. The free energy curves along the proton transfer coordinate illustrate that proton transfer is thermodynamically and kinetically favorable on the lower-energy excited state but not on the higher-energy excited state, supporting the characterization of these states as EPT and ICT, respectively. The nonadiabatic dynamics simulations indicate that the population decays from the ICT state to the EPT state in ∼100 fs and from the EPT state to the ground state on the slower time scale of ∼1 ps, qualitatively consistent with the experimental measurements. For ∼54% of the trajectories, the proton transfers from the phenol to the amine in ∼400 fs on the EPT state and then transfers back to the phenol rapidly upon decay to the ground state. Thus, these calculations augment the original interpretation of the experimental data by providing evidence of proton transfer on the EPT state prior to decay to the ground state. The fundamental insights obtained from these simulations are also relevant to other photo-EPT processes.

1. INTRODUCTION

photoinduced EPT mechanisms, denoted photo-EPT, with the aim of long-term applications to photoelectrochemical cells.1,12 Organic hydrogen-bonded dyes, as well as inorganic metal− ligand complexes, have been investigated and demonstrated to undergo photo-EPT mechanisms.1,12−16 A particular system that has been studied experimentally is the hydrogen-bonded complex composed of p-nitrophenylphenol and t-butylamine, solvated in 1,2-dichloroethane (DCE).17 This complex is depicted schematically in Figure 1. On the basis of transient absorption spectra, two distinct excited

Proton-coupled electron transfer (PCET) reactions are ubiquitous in energy conversion processes throughout chemistry and biology. Examples of such energy conversion processes are the light-driven reduction of CO2 to carbohydrates, accompanied by water oxidation, in natural photosynthesis and light-driven water splitting or water reduction to methanol, methane, or other hydrocarbons in dye-sensitized photoelectrochemical cells.1−7 PCET processes in which the proton and electron transfer reactions occur sequentially often involve high-energy intermediates that are difficult to access. In contrast, concerted mechanisms, denoted electron-proton transfer (EPT), often lead to more thermodynamically favorable pathways.8−11 Recently, experimental efforts have been directed toward the design of systems that exhibit © XXXX American Chemical Society

Special Issue: Photoinduced Proton Transfer in Chemistry and Biology Symposium Received: December 19, 2014

A

DOI: 10.1021/jp5126969 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

the EPT state and, more generally, to provide atomic-level insights into the relaxation pathways following photoexcitation to either the ICT or the EPT state. Previously our group21 used time-dependent density functional theory (TDDFT)22 with long-range corrected functionals and the polarizable continuum model23 to investigate the lowlying singlet states of the hydrogen-bonded complex composed of p-nitrophenylphenol and methylamine in DCE. In this study, the S1 state at the Franck−Condon geometry corresponded to a ππ* excited state with charge transfer character and a change in electronic density on the amine. The minimum energy structure on this S1 state corresponded to the proton bonded to the nitrogen acceptor, thereby favoring proton transfer. The S2 state at the Franck−Condon geometry corresponded to an nπ* excited state with little charge transfer character and negligible change in electronic density on the amine. The minimum energy structure for this S2 state corresponded to the proton bonded to the oxygen donor and therefore did not favor proton transfer. These TDDFT calculations also provided evidence of an avoided crossing between the S1 and S2 states near the Franck−Condon geometry. Thus, the S1 and S2 states in this study are most likely associated with the states denoted EPT and ICT, respectively, in the experimental analysis. The negligible charge transfer character of this so-called ICT state may be due to limitations of TDDFT in describing charge transfer states as well as limitations of the dielectric continuum description of the solvent. While this previous study provided insight into the nature of the electronic states, dynamical studies in conjunction with a multiconfigurational electronic structure method and an explicit solvent treatment are required to obtain reliable information about relaxation pathways and time scales. In the present study, we perform nonadiabatic dynamics simulations to gain an understanding of the time evolution of this system following photoexcitation. We use a reparameterized semiempirical implementation of the floating occupation molecular orbital complete active space configuration interaction (FOMOCASCI) method24 to generate the potential energy surfaces on the fly during the dynamics. The FOMO-CASCI method includes the required multireference character, as well as dynamical correlation through a semiempirical treatment with parameters fit to multistate complete active space second order perturbation theory (CASPT2) results.25 Solvation effects are included by an explicit treatment of the solvent molecules in a quantum mechanical/molecular mechanical (QM/MM) framework.26 Finally, the molecular dynamics with quantum transitions (MDQT) surface hopping method27 is used to describe the nonadiabatic dynamics. These simulations assist in the interpretation of the experimental data, particularly the issue regarding proton transfer on the EPT state following photoexcitation. The outline of this paper is as follows. Section 2 summarizes the computational methodology. In section 3, we present and discuss the results, starting with characterization of the excited electronic states and reparameterization of the semiempirical method in the gas phase, followed by the calculation of absorption spectra, free energy profiles for proton transfer, and nonadiabatic dynamics trajectories in DCE solution. Section 4 summarizes the results and provides concluding remarks.

Figure 1. Experimentally studied hydrogen-bonded complex composed of p-nitrophenylphenol and t-butylamine. For efficient high-level computations, t-butylamine is replaced with ammonia in this study. Use of this reduced system is validated by data for the full complex provided in the Supporting Information.

Figure 2. Schematic illustration of the electronic states probed by experiments performed on the hydrogen-bonded complex shown in Figure 1. Upon vertical excitation to the intramolecular charge transfer (ICT) state, the proton remains covalently bonded to O. In contrast, vertical excitation to the electron proton transfer (EPT) state leads to electron transfer accompanied by a shift in electronic charge density from the O−H bond to the N−H bond, resulting in an elongated N− H bond.

electronic states have been characterized for this system. Figure 2 provides a schematic illustration of these two electronic states. Both states were interpreted as intramolecular charge transfer (ICT) states involving electron transfer within the p-nitrophenylphenol, but one of these states was also interpreted to involve proton transfer and was therefore denoted an EPT state. In the conventional ICT state, the proton remains bonded to the donor oxygen, whereas in the EPT state, the proton becomes bonded to the acceptor nitrogen in that the electronic density has shifted from the O−H bond to the N−H bond. Thus, photoexcitation to the EPT state at the equilibrium ground state geometry (i.e., the Franck−Condon geometry) involves concerted electron and proton transfer upon optical excitation. According to the Franck−Condon principle,18−20 the hydrogen nucleus remains fixed during the optical excitation, forming a highly elongated N−H bond that subsequently relaxes to equilibrium. Although both singlet and triplet states were observed in the transient absorption spectra, the ultrafast dynamics is dominated by the singlet states because relaxation to the ground state via triplet states occurs on a time scale of hundreds of picoseconds.1,17 Spectra observed ∼1 ps after photoexcitation indicated that both the ICT and EPT states are populated by 355 nm excitation but only the EPT state is populated by 388 nm excitation. This observation implies that the EPT state is lower in energy than the ICT state at the Franck−Condon geometry and that the population decays from the ICT state to the EPT state on a time scale less than 1 ps. Moreover, the population decay from the EPT state was found to occur on a slower time scale of 4.5 ps. Coherent Raman experiments performed on this system shortly after photoexcitation to the EPT state provided evidence of an elongated N−H bond at the ground state equilibrium geometry. These experiments were interpreted to imply that proton transfer from the oxygen to the nitrogen occurs on the EPT state prior to relaxation to the ground state, but this proton transfer reaction could not be detected experimentally. Herein we perform nonadiabatic molecular dynamics simulations to examine the hypothesis that proton transfer occurs on

2. METHODS This section presents the methods implemented for the on-thefly QM/MM nonadiabatic dynamics simulations. The first two B

DOI: 10.1021/jp5126969 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B subsections describe the ab initio methods used to characterize the electronic states of this system and the reparameterization of the semiempirical PM3 potential for use with FOMOCASCI. The third subsection provides the details of the solvation and equilibration procedures, as well as the QM/MM potential used for the dynamical simulations. The last three subsections describe the methods used to generate the potential of mean force for each electronic state, the nonadiabatic dynamics surface hopping trajectories, and the ground state electronic absorption spectrum. The experimentally studied system is the hydrogen-bonded complex comprised of p-nitrophenylphenol and t-butylamine17 (Figure 1). For more efficient high-level computations, the amine is replaced by an ammonia molecule in this study because the alkyl substituent on the amine is not expected to be crucial to the PT process. To verify this assumption, we also performed some of the PMF calculations and nonadiabatic dynamics simulations for the complex containing t-butylamine. The corresponding results are similar to those for the reduced system and are reported in the Supporting Information (Figures S8 and S9). On the basis of this agreement, the converged nonadiabatic dynamics simulations and complete analysis were performed on the reduced system. 2.1. Ab Initio Methods for Characterization of Electronic States. The optimized geometry on the ground state, denoted the Franck−Condon (FC) geometry, of the hydrogen-bonded complex was obtained with second-order Møller−Plesset perturbation theory (MP2).28 The 6-31+G(d) basis set, with additional diffuse and polarization basis functions placed on the hydrogen atom at the hydrogen-bonding interface, was used for all ab initio calculations in this work. A single point B3LYP29,30 calculation at the geometry optimized with MP2 provided the molecular orbitals used as the initial guess for the complete active space self-consistent field31−33 (CASSCF) calculations at this geometry. For all other geometries, the initial guess for the CASSCF calculations used orbitals from CASSCF calculations at the FC geometry or similar geometries. All CASSCF and CASPT234 calculations were performed using MOLCAS 7.8,35−37 and the MP2 calculations were performed using GAUSSIAN 09.38 Several different active spaces and state-averaging procedures were examined for the CASSCF and CASPT2 calculations. The orbitals constituting an active space of 10 electrons in 10 orbitals, denoted (10,10), were found to be independent of the initial ordering of the B3LYP orbitals and the number of singlet electronic states (4 or 6) involved in the state-averaging39,40 (SA) procedure. For the SA-CASSCF calculations, we used static weighting with equal weights for all of the states. The two lowest-energy bright singlet excited electronic states from 4SACAS(10,10)SCF were found to be the same as those from 5SACAS(14,14)SCF. The multistate CASPT225 (MS-CASPT2) vertical excitation energies for these electronic states were also in very good agreement for the two sets of calculations. During potential energy scans along various internal coordinates, the (10,10) active space was found to be stable for all geometries explored. Hence, all reported scans were performed using this active space. A level shift of 0.25 Hartree was used in the MSCASPT2 calculations. Several different active spaces were also examined for the semiempirical FOMO-CASCI calculations. We performed a variety of test calculations using the FOMO-CASCI method24 with PM341 parameters developed specifically for benzene excited state dynamics, denoted RPM3.42 On the basis of these

test calculations, an active space of 6 electrons in 4 orbitals (6,4) was chosen. Note that FOMO-CASCI calculations in conjunction with an appropriate Gaussian width do not require as large of an active space as do CASSCF calculations because the FOMO procedure centers the most important orbitals involved in electronic transitions around the Fermi level. The two excited electronic states of interest are S2 and S3 with 4SACAS(10,10)SCF but are S1 and S2 with FOMO-CAS(6,4)CI. For simplicity, we refer to the three singlet electronic states of interest as S0, S1, and S2 throughout this work. All FOMOCASCI calculations were carried out using a modified version of MOPAC 2000, version 1.0.0.43 2.2. PM3 FOMO-CASCI Reparameterization. We found that the RPM3 parameters, which were developed for the excited state dynamics of benzene,42 were not sufficiently accurate for the molecule in the present study. As a result, we performed a substantial reparameterization of PM3 to ensure that the semiempirical FOMO-CASCI method reproduces the ab initio MS-CASPT2 and MP2 results for the important regions of the potential energy surface. The following properties were considered in this reparameterization procedure: (1) the optimized geometry on the ground state, (2) the vertical excitation energies at the FC geometry, (3) the potential energy as a function of the N−O distance altered by shifting the rigid p-nitrophenylphenol and ammonia fragments, (4) the proton potential energy curves along the N−O axis at different N−O distances, and (5) the potential energy curves as a function of the dihedral angle between the planes of the two benzene rings. To properly describe the balance between N−O repulsion in hydrogen-bonded and non-hydrogen-bonded geometries, we also included some points along the potential energy curve as a function of the N−O distance for the nonhydrogen-bonded complex. More details about the protocol followed for the generation of geometries corresponding to rigid scans along the internal coordinates, as well as the twostep fitting procedure, are provided in the Supporting Information. The target function that was minimized is the root-meansquare deviation (RMSD) of the semiempirical FOMO-CASCI data from the ab initio target data N

f (X ⃗ ) =

∑i = 1 [(yt, i − ys, i )wi]2 N

(1)

where X⃗ is a vector that contains all semiempirical parameters to be optimized, yt,i and ys,i are target and semiempirical values, respectively, wi is a weight, and N is the number of data points involved in the optimization. A weight of 4.0 was assigned to each of the vertical excitation energies, while a weight of 1.0 was assigned to each of the other quantities. The final set of parameters obtained is denoted R2PM3. 2.3. Solvation, Equilibration, and QM/MM Methods. The solvent was described by explicit DCE molecules using the general AMBER force field (GAFF).44,45 The solvent force field is an all-atom, flexible model. The hydrogen-bonded complex was optimized on the ground state, S0, in the gas phase using R2PM3. This optimized complex was placed at the center of a box of pre-equilibrated DCE molecules. Solvent molecules with any heavy atom within 2.8 Å of any solute heavy atom, as well as solvent molecules with no atoms within 20 Å from the center, were deleted. This procedure led to a total of 258 DCE molecules around the solute. The GAFF46 Lennard-Jones parameters were used for the solute atoms in all of the QM/ C

DOI: 10.1021/jp5126969 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

avoid unphysical initial geometries resulting from the very lowfrequency normal modes, the seven normal modes with frequencies lower than 100 cm−1 were not included in this procedure. Excluding these low-frequency modes is reasonable because they contribute negligibly to the zero point energy. Following the approach in ref 59, a 5 ps equilibration trajectory of the solvent with frozen solute was performed using Brownian dynamics for each of the sampled geometries of the solute. These trajectories were propagated at 296 K on the ground electronic state generated at the QM/MM R2PM3 FOMOCASCI level. The geometries and associated solvent configurations obtained after solvent equilibration were used as initial configurations for surface hopping trajectories starting on both the S1 and S2 states. All degrees of freedom in the solute and solvent were propagated in the surface hopping trajectories. The solute initial velocities were sampled from the quantum harmonic oscillator distributions, while the solvent initial velocities were sampled from a Maxwell−Boltzmann distribution at 296 K. The MDQT trajectories were propagated in the microcanonical (NVE) ensemble with a time step of 0.1 fs. We propagated 230 trajectories initiated on each state, leading to a total of 460 trajectories, and confirmed that the population decay times were converged. 2.6. Electronic Absorption Spectrum. For each of the solute geometries obtained from normal mode sampling on the ground state, vertical excitation energies and electronic transition dipole moments for S1 and S2 were calculated. These calculations were performed both in the gas phase and in the presence of MM solvent. For each solute geometry, the solvent coordinates were obtained at the end of a 5 ps trajectory around the frozen solute, as used in the equilibration procedure described above. For each equilibrated configuration, a Gaussian function centered at the vertical excitation energy was constructed. The height of the Gaussian function was set equal to the Einstein coefficient for absorption, which is proportional to the square of the magnitude of the electronic transition dipole moment, and was assigned a full width at halfmaximum of 20 nm. The sum of all of the Gaussian functions plotted as a function of the transition energy provides a semiclassical approximation to the inhomogeneous absorption spectrum.61

MM simulations. To restrain the solvent molecules within a sphere of radius 20 Å about the center, the CAP facility in AMBER was used with a force constant of 2 kcal mol−1 Å−2. All calculations that included solvent were carried out using an interface between MOPAC and AMBER 6.47 The detailed equilibration procedure for the solvent molecules is described in the Supporting Information. 2.4. Potential of Mean Force Calculations. The potential of mean force (PMF) for PT between N and O on each electronic state of interest was calculated using the umbrella sampling48 method. The potential energy surface was generated with the QM/MM R2PM3 FOMO-CASCI method. The reaction coordinate for these umbrella sampling calculations was the difference between the N−H and O−H distances (i.e., rNH−rOH). We employed 13−18 windows spaced by 0.1 or 0.2 Å. The force constant used for the restraint was 100 kcal mol−1 Å−2. We performed ∼10 ps of equilibration followed by ∼50 ps of production sampling per window. Each PMF calculation was performed on a single electronic adiabatic surface without any nonadiabatic transitions. The weighted histogram analysis method49,50 (WHAM) was used to obtain the final PMF for each state. 2.5. Surface Hopping Simulations. The nonadiabatic dynamics simulations were performed with the molecular dynamics with quantum transitions (MDQT) surface hopping method.27 In this approach, an ensemble of trajectories is initiated on an excited electronic state. Each trajectory moves classically on a single electronic state except for instantaneous transitions that are incorporated with the probabilistic fewestswitches algorithm.51 The time-dependent wave function is expanded as a function of the adiabatic electronic states, and the time-dependent Schrödinger equation is propagated in parallel with the classical equations of motion to determine the quantum amplitudes associated with each electronic state at each classical time step. This procedure relies on the nonadiabatic couplings between the adiabatic states, which were calculated numerically as prescribed elsewhere.27 In our implementation, the overlaps between pairs of determinants composed of nonorthogonal molecular orbitals were calculated using previously described methods.52 If a nonadiabatic transition occurred, the nonadiabatic coupling vector was calculated analytically,53,54 comprised of both a configuration interaction and a molecular orbital contribution, to determine the velocity scaling vector required for energy conservation. The probabilistic algorithm used to determine when the nonadiabatic transitions occur is designed to ensure that the population of each state at each time is the quantum probability determined from the quantum amplitudes. Note that decoherence effects55−57 are not expected to be significant in this study because of the fast time scale of these trajectories and the observation that very few trajectories exhibit multiple nonadiabatic transitions back and forth between two states. The initial conditions for the surface hopping trajectories were chosen to represent the FC geometry of the solute with the solvent equilibrated to the ground state. The zero point energy contributions from the high-frequency vibrational modes, such as the O−H stretch, are expected to be important for the dynamics of this system. Thus, zero point energy contributions were included in the initial conditions by sampling coordinates and velocities for the solute from quantum harmonic oscillator distributions corresponding to the R2PM3 normal modes calculated at the ground state geometry optimized with R2PM3 in the gas phase.58−60 To

3. RESULTS AND DISCUSSION 3.1. Characterization of Electronic States with CASSCF/CASPT2. We analyzed the two lowest-energy bright singlet states obtained at various levels of theory. For consistency throughout the paper, we refer to these two states as S1 and S2, even though the specific numbering of the states differs for each level of theory due to the presence of dark states. Figure 3 illustrates the nature of these states obtained from 4SA-CAS(10,10)SCF calculations at the FC geometry. Table 1 reveals that the S1 state is characterized by a much larger dipole moment compared to the S0 and S2 states, thereby implying a significant amount of intramolecular charge transfer in the S1 state. The S2 state also has a larger dipole moment than the ground state, indicating a degree of intramolecular charge transfer in this excited state also. The 5SA-CASSCF calculations with a larger active space of (14,14) provide a qualitatively similar picture. The vertical excitation energies obtained from MS-CASPT2 calculations are also very consistent for the two active spaces. D

DOI: 10.1021/jp5126969 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Table 2. Values of Key Quantities Used in the Reparameterization to Obtain the R2PM3 Parameter Set targeta

R2PM3b

Vertical Excitation Energy (eV) Eexc S1 4.46 4.35 Eexc S2 5.22 5.21 At rN−O = 2.808 Å, Energy Relative to rN−H = 1.8 Å (eV)c rN−H = 1.0 Å 0.97 0.94 rN−H = 1.4 Å 0.82 0.82 rN−H = 1.95 Å 0.53 0.58 At rN−O = 2.7 Å, Energy Relative to rN−H = 1.7 Å (eV)c rN−H = 1.0 Å 0.96 0.84 rN−H = 1.35 Å 0.60 0.70 rN−H = 1.85 Å 0.63 0.64 At rN−O = 2.6 Å, Energy Relative to rN−H = 1.6 Å (eV)c rN−H = 1.0 Å 0.89 0.74 rN−H = 1.3 Å 0.42 0.56 rN−H = 1.75 Å 0.65 0.60 With a H-Bond Present, Energy Relative to rN−O = 2.8 Å (eV)c rN−O = 2.6 Å 0.03 0.09 rN−O = 3.0 Å 0.05 0.06 Energy Relative to a Twist Angle of 45° (in eV)c 5° −0.06 0.00 85° 0.66 0.30 RMSD of Optimized Geometry on S0 with Respect to the Targetd (in Å) S0 0.00 0.296

Figure 3. Isosurfaces of electronic density differences between the bright singlet excited states and the ground state of the pnitrophenylphenol−ammonia complex in the gas phase at the MP2 optimized ground state geometry. The density differences correspond to the (a) S0 → S1 (EPT) and (b) S0 → S2 (ICT) transitions. The electronic densities were obtained from 4SA-CAS(10,10)SCF calculations. Isosurfaces of 0.0001 and −0.0001 bohr−3 are depicted in green and pink, respectively. The EPT transition exhibits greater changes in electronic density near the hydrogen-bonding interface.

For further comparison, we also performed TDDFT calculations with the long-range corrected functional CAMB3LYP.62 The TDDFT results lead to excited electronic states with similar dipole moments, but the vertical excitation energies are underestimated, possibly due to an inadequate treatment of the charge transfer character of these states.63 These results are qualitatively consistent with our previous TDDFT calculations on p-nitrophenylphenol in the gas phase.21 Only the CASSCF and CASPT2 results were used for the reparameterization procedure described in the next subsection. 3.2. Optimization of PM3 Parameters. We reparameterized the PM3 potential for application of the semiempirical FOMO-CASCI method to the p-nitrophenylphenol−ammonia hydrogen-bonded complex. The new parameters are given in Table S2 of the Supporting Information. As depicted in Figure S6 (Supporting Information), the ground state geometries optimized with MP2 and R2PM3 are in good agreement, with a root-mean-square deviation (RMSD) of 0.296 Å. The key energetic properties of this system obtained using the new PM3 parameters, denoted R2PM3, are given in Table 2. The vertical excitation energies obtained with R2PM3 and MS-CASPT2 are in very good agreement. The accuracy of the proton potential energy curves is critical for studying the PT reaction. As illustrated in Figure 4, the proton potential energy curves generated with R2PM3 FOMOCASCI are in reasonable agreement with the MS-CASPT2

a

Target values were obtained from MS-4SA-CAS(10,10)PT2 unless otherwise specified. The rigid scans along various coordinates were based on the geometry optimized with MP2 on S0. bSemiempirical results were obtained using FOMO-CAS(6,4)CI unless otherwise specified. The vertical excitation energies were calculated for the geometry optimized with R2PM3 on S0. The other energy calculations were performed on the same structures as used to obtain the target values. The Gaussian width parameter ω = 0.18 au was used for all FOMO-CASCI calculations. cRelative energies were calculated on the S1 state. dThe MP2 geometry was used as a target.

results for the ground state equilibrium N−O distance of 2.808 Å. The deviations on the S0 and S2 states occur in the higherenergy regions that are unlikely to be encountered in dynamical trajectories. The agreement for the S1 state is better because several energy values along this curve were used in the reparameterization procedure. At N−O distances of 2.6 and 2.7 Å, which are shorter than the N−O distance at the FC geometry, the R2PM3 FOMO-CASCI PT barrier on S1 is higher than the MS-CASPT2 barrier by ∼0.1 eV (Figure S1, Supporting Information), but the overall agreement is still reasonable. The potential energy dependence on the N−O distance is also relevant to the PT reaction. This property is particularly

Table 1. MS-CASPT2 and TD-CAM-B3LYP Results for the Three Singlet Electronic States of Interest for the HydrogenBonded p-Nitrophenylphenol−Ammonia Complex in the Gas Phase MS-4SA-CAS(10,10)PT2d,e S0 S1 S2

MS-5SA-CAS(14,14)PT2d,f

TD-CAM-B3LYPg

fa

Eexcb (eV)

|μ|c (D)

fa

Eexcb (eV)

|μ|c (D)

fa

Eexcb (eV)

|μ|c (D)

0.3677 0.0294

4.46 5.22

9.0 22.2 15.6

0.6715 0.0328

4.46 5.19

8.3 18.5 13.3

0.5349 0.0117

4.04 4.73

9.9 20.5 13.6

a Oscillator strength. bVertical excitation energy. cDipole moment; values reported under CASPT2 were obtained from CASSCF calculations. dSingle point calculation at the MP2 optimized geometry. eS1 and S2 are actually S2 and S3 in the calculations. fS1 and S2 are actually S3 and S4 in the calculations. gSingle point calculation at the CAM-B3LYP optimized geometry. S1 and S2 are actually S2 and S4 in the calculations.

E

DOI: 10.1021/jp5126969 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 5. Calculated electronic absorption spectra of the pnitrophenylphenol−ammonia complex. The solid curves represent the spectra in 1,2-dichloroethane (DCE) solution, and the dashed curves represent the spectra in the gas phase. The wavelengths corresponding to the maximum absorption for S1 in DCE, S1 in the gas phase, S2 in DCE, and S2 in the gas phase are 308, 294, 252, and 246 nm, respectively.

Figure 4. Proton potential energy curves in which the proton is moved along the N−O axis for the p-nitrophenylphenol−ammonia complex at a fixed N−O distance of 2.808 Å in the gas phase. The MS-CASPT2 calculations were performed on ground state structures based on the geometry optimized with MP2, while the R2PM3 FOMO-CASCI calculations were performed on ground state structures based on the geometry optimized with R2PM3. The dashed green curve coincides with the dashed red curve and therefore is not visible. The corresponding figures at N−O distances of 2.7 and 2.6 Å are provided in the Supporting Information (Figure S1).

with the experimental interpretation, as discussed further below. Experimental measurements indicated that the formation of the hydrogen-bonded complex with t-butylamine caused a red shift of ∼20 nm in the excitation wavelength of p-nitrophenylphenol.17 A similar trend is observed in our simulations, as depicted in Figure S7 (Supporting Information). Using data from ground state trajectories in which the hydrogen bond between p-nitrophenylphenol and ammonia breaks and reforms sporadically, the formation of the hydrogen bond increases the excitation wavelength by ∼20−30 nm for the S1 state, as observed experimentally. 3.4. Proton Transfer Free Energy Profiles in Solution. While the gas phase proton potential energy curves (Figure 4 and Figure S1, Supporting Information) show that PT from O to N is not favored on any of the electronic states, the original interpretation of the experimental data suggests that PT may occur on S1. Moreover, our previous work,21 which utilized TDDFT with implicit solvent, predicted PT from O to N to be favorable on S1. Here we investigate the feasibility of PT from O to N on S1 in an explicit solvent environment by calculating the free energy profile for proton transfer. For this purpose, we use the QM/MM R2PM3 FOMO-CASCI method to calculate the PMF as a function of the PT coordinate rNH−rOH. To gain further insight into how the S1 state differs from the S0 and S2 states, we also calculate the PMFs for the S0 and S2 states. The PMF curves depicted in Figure 6 illustrate that the proton is thermodynamically favored to bond to O on the S0 and S2 states but is thermodynamically favored to bond to N on the S1 state. Moreover, the free energy barrier for PT from O to N is only ∼0.17 eV (∼4 kcal/mol) on the S1 state and is 0.95 eV (∼22 kcal/mol) on the S0 and S2 states. As a result, the proton is expected to remain bonded to O after photoexcitation to the S2 state but can transfer to N after photoexcitation to the S1 state. The differences between the gas phase proton potential energy curve (Figure 4) and the PMF in solution (Figure 6) for the S1 state illustrate the importance of the solvent for enabling the PT reaction. Solvent effects stabilize the configuration with the transferred proton more than the initial configuration because the former has ionic character, with a negative charge on the phenol oxygen and a positive charge on the ammonium. This ionic configuration is more stabilized by solvent effects than is the initial configuration. The S1 state has a smaller energy difference between the two proton positions than do the S0 and S2 states, leading to the PMFs depicted in Figure 6.

important for the S1 state, on which the PT reaction occurs. The R2PM3 FOMO-CASCI potential energy scan along the N−O distance in the hydrogen-bonded configuration on S1 (Figure S3, Supporting Information) agrees well with the MSCASPT2 results near the minimum and is more repulsive by only ∼0.1 eV at short N−O distances of ∼2.6 Å. At shorter N− O distances of ∼2.5 Å, the deviation is significantly greater, but such short N−O distances are not encountered often in the trajectories. The R2PM3 FOMO-CASCI potential energy scans along the N−O distance in the non-hydrogen-bonded configuration are also in qualitative agreement with the MSCASPT2 results for the relevant N−O distances (Figure S4, Supporting Information). The potential energy scans along the dihedral angle between the planes of the two benzene rings (i.e., the “twist” angle) are not appreciably affected by reparameterization (Figure S5, Supporting Information). The barriers for rotation are underestimated by R2PM3 FOMO-CASCI on all three electronic states. Because the magnitude of the deviation is different for each electronic state, the potential energy scans cannot be corrected by the addition of a single correction term to the Hamiltonian.64 The decay from the S2 to the S1 state is fast enough that the twist angle does not change significantly on the S2 state. The system spends more time on the S1 state, where the twist angle changes but does not play an important role in the relaxation to the S0 state. Therefore, the errors in the potential energy along this coordinate are not expected to impact the qualitative relaxation behavior. 3.3. Absorption Spectra in Solution. Figure 5 depicts the ground state absorption spectra of the p-nitrophenylphenol− ammonia hydrogen-bonded complex in the gas phase and in DCE. According to experimental measurements,17 excitation at 388 nm populates only the S1 state, and excitation at 355 nm populates both the S1 and S2 states. The calculated excitation energies from the S0 state to the S1 and S2 states are 4.0 and 4.9 eV, respectively. Moreover, relative to the gas phase, the calculated vertical excitation energies for the S1 and S2 states in solution are lowered by ∼0.2 and ∼0.1 eV, respectively, according to the absorption spectra in Figure 5. Although the absolute excitation energies are overestimated, as typically found for these types of calculations,65−67 the nature and relative ordering of these two excited states are in agreement F

DOI: 10.1021/jp5126969 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

the same time scale as did the trajectories that were initiated on the S1 state. Specifically, the S2 state population decayed to the S1 state within ∼100 fs, while the S1 state population decayed to the S0 state on a slower time scale of ∼0.9 ps. Transient absorption experiments indicated a relaxation time from the S1 state to the ground state of ∼4.5 ps.1,17 Hence, the calculated time scale of relaxation to the ground state is in qualitative agreement with the experimentally measured time scale. The quantitative difference between the calculated and experimental time scales could arise from inaccuracies in the semiempirical potential energy surfaces and vertical excitation energies, as well as the classical treatment of the transferring proton. Note that the zero point energy that was included in the initial conditions for the solute does not leak into the solvent significantly because of the short time scale of the decay, thereby maintaining the temperature of the solvent bath during the relaxation process (Figure S10, Supporting Information). In addition, the transient absorption experiments indicated that vertical excitation to the ICT state, S2, leads to population of the lower-energy EPT state within 1 ps. The observation of very fast population decay from the S2 to the S1 state in the MDQT simulations is in agreement with this experimental observation. Moreover, the calculations suggest that the fast decay from the S2 state to the S1 state occurs by way of a conical intersection (CI) between these two states. A gas phase optimization of the S2/S1 minimum energy CI at the R2PM3 FOMO-CASCI level revealed a “pyramidalization” of the NO2bearing phenyl ring (Figure 8). Moreover, the MDQT

Figure 6. Potential of mean force (PMF) for proton transfer between p-nitrophenylphenol and ammonia in DCE on the three electronic states of interest. These curves were generated with the QM/MM R2PM3 FOMO-CASCI method. The minima on the right side correspond to the proton bonded to the hydroxyl O atom of pnitrophenylphenol, and the minima on the left side correspond to the proton bonded to the N atom of ammonia.

These PMF curves indicate that S2 is an ICT state, which exhibits intramolecular charge transfer but not PT, and S1 is an EPT state, which exhibits intramolecular charge transfer as well as PT. Upon photoexcitation to the S1 state at the FC geometry, in which the proton is bound to O, the zero point energy will be sufficient to surmount the free energy barrier of ∼0.17 eV for PT from the O to the N. Although the proton is treated classically in the nonadiabatic dynamics trajectories, the zero point energies of the solute modes are included in the initial conditions. Note that the free energy differences in the PMF curves for the three states are significantly greater than the deviations observed between the R2PM3 FOMO-CASCI and CASPT2 results. Thus, the trends in the free energy curves and the characterization of the S1 and S2 states as EPT and ICT, respectively, are reliable. Moreover, this characterization of the two electronic states is consistent with the previous interpretation of the experimental data17 and the previous TDDFT calculations in implicit solvent.21 3.5. Excited State Decay Time Scales and Pathways. Figure 7 depicts the population decays from the S1 and S2 states down to the ground state, as obtained from an ensemble of surface hopping trajectories initiated on each state. The trajectories initiated on the S2 state quickly hopped down to the S1 state and then eventually hopped down to the S0 state on

Figure 8. (a) Overlay of the geometry optimized with R2PM3 on the S0 state (orange) and the geometry at the S2/S1 minimum energy conical intersection determined with R2PM3 FOMO-CASCI (blue) in the gas phase. (b) Schematic picture used to depict the values of key dihedral angles at the S2/S1 minimum energy conical intersection geometry. The thick black lines represent the central bonds in these dihedral angles, which are approximately zero at the ground state minimum energy geometry. The dihedral angle value in a particular color corresponds to the dihedral angle involving the atoms at the ends of the bonds in the same color.

trajectories in solution were also found to exhibit this “pyramidalization” during nonadiabatic transitions between the S2 and S1 states. Note that this type of pyramidalization has also been observed in gas phase CI geometries of ethylene and benzene,42,68 although the minimum energy CI has a half boat-like structure in benzene but a chair-like structure for the NO2-bearing phenyl ring in the present study. The experimentally measured Raman spectrum following photoexcitation to the S1 state indicated an elongated N−H bond suggestive of PT,17 thereby justifying the characterization of the S1 state as EPT. As discussed above, the PMF for the S1

Figure 7. Population decay of trajectories initiated on the S1 (upper panel) and S2 (lower panel) electronic states, as obtained from the MDQT simulations. The inset in the lower panel shows the fast population decay from S2 that occurs within ∼100 fs. G

DOI: 10.1021/jp5126969 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B state (Figure 6) is consistent with this EPT characterization. In addition, transient absorption experiments have been interpreted to indicate a time scale of ∼200−250 fs for PT on the EPT state. In the MDQT simulations, ∼54% of the trajectories exhibited excited state PT on the S1 state with an average time scale of ∼400 fs, where proton transfer was defined as the N−H distance becoming shorter than 1.15 Å. Hence, in a majority of the trajectories, the highly elongated N−H bond relaxed to its equilibrium length before the electronic transition to the ground state. As illustrated by the histogram in Figure 9,

Figure 11. Representative MDQT trajectory initiated on S2. The decay from the S2 to the S1 state at ∼164 fs is followed by proton transfer from O to N on the S1 state at ∼600 fs. Upon decay to the S0 state at ∼745 fs, the proton transfers back to O. The O−H and N−H distances are calculated for the proton that is involved in the hydrogen bond. Figure 9. Histogram corresponding to the calculated time for the first proton transfer between p-nitrophenylphenol and ammonia in DCE after photoexcitation to state S1. This histogram was constructed from the trajectories that exhibited proton transfer on S1, which constituted ∼57% of all MDQT trajectories initiated on state S1. Proton transfer was defined to occur when the proton involved in the hydrogen bond was within a distance of 1.15 Å from the nitrogen atom of ammonia.

minimum energy CI between S1 and S0 in solvent exhibits the proton bonded to the N atom of ammonia and pyramidalization about the N atom of the nitro group. The observation of the proton bonded to the ammonia at the minimum energy CI in contrast to the proton bonded to the phenol at the equilibrium geometry on the S0 state indicates that proton transfer plays an important role in the S1 to S0 decay. Other points on the conical intersection seam could also be important for the decay process. Figure 11 illustrates the time evolution of the occupied state and proton transfer distances for a representative trajectory started on S2. In this trajectory, fast decay from the S2 state to the S1 state is followed by proton transfer from p-nitrophenylphenol to ammonia on the S1 state. When the proton is bonded to the nitrogen, the S1 state is stabilized while the S0 state is destabilized, as shown by the PMFs in Figure 6. This phenomenon leads to a smaller splitting between these two states, thereby facilitating the decay from the S1 state to the S0 state. Following this decay, the proton rapidly transfers back to the phenol. This trajectory exemplifies the behavior that was suggested by the transient absorption and Raman experiments on this system.

however, the spread of time scales for PT is quite broad. Furthermore, analysis of the trajectories indicated a coupling between PT and decay to the ground electronic state, suggesting that the configuration in which the proton is bonded to the nitrogen is low in energy on S1 but high in energy on S0. In other words, the splitting between these two states is relatively small at this configuration, thereby facilitating a nonadiabatic transition down to the ground state. This observation is consistent with the PMFs associated with PT on these states (Figure 6), namely, that the S1 state is stabilized and the S0 state is destabilized by PT to the N. In addition, population decay from the S1 state to the S0 state was found to be accompanied by PT back to the phenol, as shown to be thermodynamically favorable in Figure 6. To further investigate this PT reaction, we optimized the conical intersection starting from the solvent configuration obtained at the end of the QM/MM solvent equilibration around the frozen solute at the FC geometry. As shown in Figure 10, the

4. CONCLUSIONS In this paper, we characterized the excited electronic states of the p-nitrophenylphenol−ammonia system with the CASSCF and CASPT2 methods and used these results to parameterize a semiempirical FOMO-CASCI method to enable QM/MM molecular dynamics simulations of this molecule in explicit solvent. Our calculations identified two charge transfer excited states with properties that are in qualitative agreement with experimental observations. Both the S1 and S2 states exhibit charge-transfer character within the phenyl rings, but only the S1 state exhibits changes in electronic charge density near the proton transfer interface. In the gas phase, both of these states have minima corresponding to the proton bonded to the oxygen of the phenol and do not favor proton transfer to the nitrogen. This situation is qualitatively different in solution, however, as indicated by the calculated potentials of mean force along the proton transfer coordinate. For the S1 state, the

Figure 10. Minimum energy conical intersection between S0 and S1 in the presence of solvent. The optimization was started from a solvent configuration obtained at the end of the QM/MM solvent equilibration around the frozen solute at the FC geometry. H

DOI: 10.1021/jp5126969 J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B



proton is thermodynamically favored to bond to the nitrogen by ∼0.30 eV (∼7 kcal/mol), with a free energy barrier of only ∼0.17 eV (∼4 kcal/mol) for proton transfer from the oxygen to the nitrogen. For the S2 state, the proton is thermodynamically and kinetically favored to remain on the oxygen. These analyses support the characterization of the S1 state as EPT and the S2 state as ICT and are consistent with the experimental coherent Raman measurements indicating an elongated N−H bond upon excitation to the S1 state. We also performed MDQT surface hopping simulations to investigate the ultrafast dynamics following excitation to either the S1 or the S2 state. Our calculations indicated that decay from the S2 to the S1 state occurs in ∼100 fs and decay from the S1 to the S0 state occurs on the slower time scale of ∼1 ps. These time scales are qualitatively consistent with the experimental measurements, which provided evidence of decay from the S2 state to the S1 state in less than 1 ps and decay from the S1 state on a time scale of ∼4.5 ps. The simulations suggested that the decay from the S2 to the S1 state occurs via a conical intersection associated with pyramidalization of the NO2-bearing phenyl ring. Moreover, for ∼54% of the trajectories, the proton transfers from the oxygen to the nitrogen in ∼400 fs on the S1 state and then transfers back to the oxygen quickly upon decay to the S0 state. The experiments were unable to detect the actual transfer of the proton from the oxygen to the nitrogen. Complementing these experiments, the simulations provide compelling evidence that the proton transfer reaction occurs on the S1 state prior to decay to the ground state. Thus, the calculations clarify and augment the original interpretation of the experiments. Photo-EPT processes encompass both intramolecular electron transfer and proton transfer on a single excited electronic state upon photoexcitation. Understanding the underlying physical principles of these processes has important implications for a wide range of energy conversion systems. QM/MM nonadiabatic dynamics simulations provide atomiclevel insights that can assist in the interpretation of experimental data and generate experimentally testable predictions. The simulations described herein treat the transferring hydrogen nucleus classically in the nonadiabatic dynamics trajectories, although the zero point energies of the solute modes are included in the initial conditions. Future simulations will treat the transferring hydrogen nucleus quantum mechanically during the nonadiabatic dynamics trajectories69−71 and will be utilized to calculate hydrogen/ deuterium isotope effects72 to stimulate further experiments. This feedback between experiments and simulations will play a vital role in elucidating fundamental aspects of photo-EPT processes.



Article

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS This material is based upon work supported by the Air Force Office of Scientific Research under AFOSR Award No. FA9550-14-1-0295. We are grateful to Todd Martinez and James Stewart for providing a modified MOPAC code that includes the semiempirical FOMO-CASCI QM/MM method.



REFERENCES

(1) Gagliardi, C. J.; Westlake, B. C.; Kent, C. A.; Paul, J. J.; Papanikolas, J. M.; Meyer, T. J. Integrating proton coupled electron transfer (PCET) and excited states. Coord. Chem. Rev. 2010, 254, 2459−2471. (2) Gust, D.; Moore, T. A.; Moore, A. L. Solar Fuels via Artificial Photosynthesis. Acc. Chem. Res. 2009, 42, 1890−1898. (3) Magnuson, A.; Anderlund, M.; Johansson, O.; Lindblad, P.; Lomoth, R.; Polivka, T.; Ott, S.; Stensjo, K.; Styring, S.; Sundstrom, V.; Hammarstrom, L. Biomimetic and Microbial Approaches to Solar Fuel Generation. Acc. Chem. Res. 2009, 42, 1899−1909. (4) Prezhdo, O. V.; Duncan, W. R.; Prezhdo, V. V. Dynamics of the photoexcited electron at the chromophore-semiconductor interface. Acc. Chem. Res. 2008, 41, 339−348. (5) Sproviero, E. M.; Gascon, J. A.; McEvoy, J. P.; Brudvig, G. W.; Batista, V. S. Quantum mechanics/molecular mechanics structural models of the oxygen-evolving complex of photosystem II. Curr. Opin. Struct. Biol. 2007, 17, 173−180. (6) Alstrum-Acevedo, J. H.; Brennaman, M. K.; Meyer, T. J. Chemical approaches to artificial photosynthesis. 2. Inorg. Chem. 2005, 44, 6802−6827. (7) Gratzel, M. Photoelectrochemical cells. Nature 2001, 414, 338− 344. (8) Huynh, M. H. V.; Meyer, T. J. Proton-coupled electron transfer. Chem. Rev. 2007, 107, 5004−5064. (9) Hammes-Schiffer, S.; Soudackov, A. V. Proton-Coupled Electron Transfer in Solution, Proteins, and Electrochemistry. J. Phys. Chem. B 2008, 112, 14108−14123. (10) Hammes-Schiffer, S.; Stuchebrukhov, A. A. Theory of Coupled Electron and Proton Transfer Reactions. Chem. Rev. 2010, 110, 6939− 6960. (11) Warren, J. J.; Tronic, T. A.; Mayer, J. M. Thermochemistry of Proton-Coupled Electron Transfer Reagents and its Implications. Chem. Rev. 2010, 110, 6961−7001. (12) Concepcion, J. J.; Brennaman, M. K.; Deyton, J. R.; Lebedeva, N. V.; Forbes, M. D. E.; Papanikolas, J. M.; Meyer, T. J. Excited-state quenching by proton-coupled electron transfer. J. Am. Chem. Soc. 2007, 129, 6968−6969. (13) Eisenhart, T. T.; Dempsey, J. L. Photo-induced Proton-Coupled Electron Transfer Reactions of Acridine Orange: Comprehensive Spectral and Kinetics Analysis. J. Am. Chem. Soc. 2014, 136, 12221− 12224. (14) Chen, J.; Kuss-Petermann, M.; Wenger, O. S. Distance Dependence of Bidirectional Concerted Proton-Electron Transfer in Phenol-Ru(2,2′-bipyridine)3(2+) Dyads. Chem.Eur. J. 2014, 20, 4098−4104. (15) Wenger, O. S. Proton-Coupled Electron Transfer with Photoexcited Metal Complexes. Acc. Chem. Res. 2013, 46, 1517−1526. (16) Wenger, O. S. Proton-coupled electron transfer with photoexcited ruthenium(II), rhenium(I), and iridium(III) complexes. Coord. Chem. Rev. 2015, 282−283, 150−158. (17) Westlake, B. C.; Brennaman, M. K.; Concepcion, J. J.; Paul, J. J.; Bettis, S. E.; Hampton, S. D.; Miller, S. A.; Lebedeva, N. V.; Forbes, M. D. E.; Moran, A. M.; Meyer, T. J.; Papanikolas, J. M. Concerted

ASSOCIATED CONTENT

S Supporting Information *

Additional details of the reparameterization process, including potential energy scans along various coordinates and the R2PM3 parameter set; details of the solvent equilibration procedure; dependence of excitation wavelengths on the hydrogen-bonding interaction between p-nitrophenylphenol and ammonia; and simulation results for the hydrogen-bonded complex between p-nitrophenylphenol and t-butylamine. This material is available free of charge via the Internet at http:// pubs.acs.org. I

DOI: 10.1021/jp5126969 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B electron-proton transfer in the optical excitation of hydrogen-bonded dyes. Proc. Natl. Acad. Sci. U.S.A. 2011, 108, 8554−8558. (18) Franck, J.; Dymond, E. G. Elementary processes of photochemical reactions. Trans. Faraday Soc. 1926, 21, 536−542. (19) Condon, E. A theory of intensity distribution in band systems. Phys. Rev. 1926, 28, 1182−1201. (20) Birge, R. T. The band spectra of carbon monoxide. Phys. Rev. 1926, 28, 1157−1181. (21) Ko, C.; Solis, B. H.; Soudackov, A. V.; Hammes-Schiffer, S. Photoinduced Proton-Coupled Electron Transfer of HydrogenBonded p-Nitrophenylphenol-Methylamine Complex in Solution. J. Phys. Chem. B 2013, 117, 316−325. (22) Casida, M. E. Time-Dependent Density Functional Response Theory for Molecules. In Recent Advances in Density Functional Methods; Chong, D. P., Ed.; World Scientific: Singapore, 1995; Vol. 1, p 155. (23) Tomasi, J.; Mennucci, B.; Cammi, R. Quantum mechanical continuum solvation models. Chem. Rev. 2005, 105, 2999−3094. (24) Granucci, G.; Toniolo, A. Molecular gradients for semiempirical CI wavefunctions with floating occupation molecular orbitals. Chem. Phys. Lett. 2000, 325, 79−85. (25) Finley, J.; Malmqvist, P. A.; Roos, B. O.; Serrano-Andres, L. The multi-state CASPT2 method. Chem. Phys. Lett. 1998, 288, 299−306. (26) Senn, H. M.; Thiel, W. QM/MM Methods for Biomolecular Systems. Angew. Chem., Int. Ed. 2009, 48, 1198−1229. (27) Hammes-Schiffer, S.; Tully, J. C. Proton Transfer in Solution: Molecular Dynamics with Quantum Transitions. J. Chem. Phys. 1994, 101, 4657−4667. (28) Møller, C.; Plesset, M. S. Note on an approximation treatment for many-electron systems. Phys. Rev. 1934, 46, 0618−0622. (29) Becke, A. D. Density-Functional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 1993, 98, 5648−5652. (30) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. Ab Initio Calculation of Vibrational Absorption and Circular Dichroism Spectra Using Density Functional Force Fields. J. Phys. Chem. 1994, 98, 11623−11627. (31) Roos, B. O.; Taylor, P. R.; Siegbahn, P. E. M. A Complete Active Space SCF Method (CASSCF) Using a Density Matrix Formulated Super-CI Approach. Chem. Phys. 1980, 48, 157−173. (32) Siegbahn, P.; Heiberg, A.; Roos, B.; Levy, B. A comparison of the Super-Cl and the Newton-Raphson Scheme in the Complete Active Space SCF Method. Phys. Scr. 1980, 21, 323−327. (33) Siegbahn, P. E. M.; Almlof, J.; Heiberg, A.; Roos, B. O. The Complete Active Space SCF (CASSCF) Method in a NewtonRaphson Formulation with Application to the HNO Molecule. J. Chem. Phys. 1981, 74, 2384−2396. (34) Andersson, K.; Malmqvist, P. A.; Roos, B. O. Second-Order Perturbation Theory with a Complete Active Space Self-Consistent Field Reference Function. J. Chem. Phys. 1992, 96, 1218−1226. (35) Aquilante, F.; De Vico, L.; Ferre, N.; Ghigo, G.; Malmqvist, P. A.; Neogrady, P.; Pedersen, T. B.; Pitonak, M.; Reiher, M.; Roos, B. O.; Serrano-Andres, L.; Urban, M.; Veryazov, V.; Lindh, R. MOLCAS 7: The Next Generation. J. Comput. Chem. 2010, 31, 224−247. (36) Karlstrom, G.; Lindh, R.; Malmqvist, P. A.; Roos, B. O.; Ryde, U.; Veryazov, V.; Widmark, P. O.; Cossi, M.; Schimmelpfennig, B.; Neogrady, P.; Seijo, L. MOLCAS: a program package for computational chemistry. Comput. Mater. Sci. 2003, 28, 222−239. (37) Veryazov, V.; Widmark, P. O.; Serrano-Andres, L.; Lindh, R.; Roos, B. O. 2MOLCAS as a development platform for quantum chemistry software. Int. J. Quantum Chem. 2004, 100, 626−635. (38) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega,

N.; Millam, M. J.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, Ö .; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09; Gaussian, Inc.: Wallingford, CT, 2009. (39) Knowles, P. J.; Werner, H. J. An Efficient Second-Order MC SCF Method for Long Configuration Expansions. Chem. Phys. Lett. 1985, 115, 259−267. (40) Werner, H. J.; Knowles, P. J. A Second Order Multiconfiguration SCF Procedure with Optimum Convergence. J. Chem. Phys. 1985, 82, 5053−5063. (41) Stewart, J. J. P. Optimization of Parameters for Semiempirical Methods I. Method. J. Comput. Chem. 1989, 10, 209−220. (42) Toniolo, A.; Thompson, A. L.; Martinez, T. J. Excited state direct dynamics of benzene with reparameterized multi-reference semiempirical configuration interaction methods. Chem. Phys. 2004, 304, 133−145. (43) Stewart, J. J. P. MOPAC 2000, version 1.0.0; Fujitsu Limited: Tokyo, Japan, 2000. (44) Caleman, C.; van Maaren, P. J.; Hong, M.; Hub, J. S.; Costa, L. T.; van der Spoel, D. Force Field Benchmark of Organic Liquids: Density, Enthalpy of Vaporization, Heat Capacities, Surface Tension, Isothermal Compressibility, Volumetric Expansion Coefficient, and Dielectric Constant. J. Chem. Theory Comput. 2012, 8, 61−74. (45) van der Spoel, D. http://virtualchemistry.org/molecules/10706-2/index.php. (46) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and testing of a general amber force field. J. Comput. Chem. 2004, 25, 1157−1174. (47) Case, D. A.; Pearlman, D. A.; Caldwell, J. W.; Cheatham, T. E., III; Ross, W. S.; Simmerling, C. L.; Darden, T. A.; Merz, K. M.; Stanton, R. V.; Cheng, A. L.; Vincent, J. J.; Crowley, M.; Tsui, V.; Radmer, R. J.; Duan, Y.; Pitera, J.; Massova, I.; Seibel, G. L.; Singh, U. C.; Weiner, P. K.; Kollman, P. A. AMBER 6; University of California, San Francisco: San Francisco, CA, 1999. (48) Torrie, G. M.; Valleau, J. P. Nonphysical Sampling Distributions in Monte Carlo Free-Energy Estimation: Umbrella Sampling. J. Comput. Phys. 1977, 23, 187−199. (49) Grossfield, A. WHAM: the weighted histogram analysis method, version 2.0.2; http://membrane.urmc.rochester.edu/content/wham. (50) Kumar, S.; Rosenberg, J. M.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A. The Weighted Histogram Analysis Method for FreeEnergy Calculations on Biomolecules. I. The Method. J. Comput. Chem. 1992, 13, 1011−1021. (51) Tully, J. C. Molecular Dynamics with Electronic Transitions. J. Chem. Phys. 1990, 93, 1061−1071. (52) Granucci, G.; Persico, M.; Toniolo, A. Direct semiclassical simulation of photochemical processes with semiempirical wave functions. J. Chem. Phys. 2001, 114, 10608−10615. (53) Fabiano, E.; Keal, T. W.; Thiel, W. Implementation of surface hopping molecular dynamics using semiempirical methods. Chem. Phys. 2008, 349, 334−347. (54) Lengsfield, B. H.; Yarkony, D. R. Nonadiabatic Interactions between Potential Energy Surfaces: Theory and Applications. Adv. Chem. Phys. 2007, 82, 1−71. (55) Webster, F.; Rossky, P. J.; Friesner, R. A. Nonadiabatic Processes in Condensed Matter: Semi-classical Theory and Implementation. Comput. Phys. Commun. 1991, 63, 494−522. (56) Subotnik, J. E.; Ouyang, W.; Landry, B. R. Can we derive Tully’s surface-hopping algorithm from the semiclassical quantum Liouville equation? Almost, but only with decoherence. J. Chem. Phys. 2013, 139, 214107. (57) Granucci, G.; Persico, M.; Zoccante, A. Including quantum decoherence in surface hopping. J. Chem. Phys. 2010, 133, 134111. (58) Robinett, R. W. Quantum and Classical Probability Distributions for Position and Momentum. Am. J. Phys. 1995, 63, 823−832. J

DOI: 10.1021/jp5126969 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B (59) Ruckenbauer, M.; Barbatti, M.; Muller, T.; Lischka, H. Nonadiabatic Excited-State Dynamics with Hybrid ab Initio Quantum-Mechanical/Molecular-Mechanical Methods: Solvation of the Pentadieniminium Cation in Apolar Media. J. Phys. Chem. A 2010, 114, 6757−6765. (60) Virshup, A. M.; Punwong, C.; Pogorelov, T. V.; Lindquist, B. A.; Ko, C.; Martinez, T. J. Photodynamics in Complex Environments: Ab Initio Multiple Spawning Quantum Mechanical/Molecular Mechanical Dynamics. J. Phys. Chem. B 2009, 113, 3280−3291. (61) Lukes, V.; Solc, R.; Barbatti, M.; Elstner, M.; Lischka, H.; Kauffmann, H. F. Torsional potentials and full-dimensional simulation of electronic absorption and fluorescence spectra of para-phenylene oligomers using the semiempirical self-consistent charge densityfunctional tight binding approach. J. Chem. Phys. 2008, 129, 164905. (62) Yanai, T.; Tew, D. P.; Handy, N. C. A new hybrid exchangecorrelation functional using the Coulomb-attenuating method (CAMB3LYP). Chem. Phys. Lett. 2004, 393, 51−57. (63) Computational Photochemistry, 1st ed.; Elsevier: Amsterdam 2005; Vol. 16. (64) Cusati, T.; Granucci, G.; Martinez-Nunez, E.; Martini, F.; Persico, M.; Vazquez, S. Semiempirical Hamiltonian for Simulation of Azobenzene Photochemistry. J. Phys. Chem. A 2012, 116, 98−110. (65) Isborn, C. M.; Gotz, A. W.; Clark, M. A.; Walker, R. C.; Martinez, T. J. Electronic Absorption Spectra from MM and ab Initio QM/MM Molecular Dynamics: Environmental Effects on the Absorption Spectrum of Photoactive Yellow Protein. J. Chem. Theory Comput. 2012, 8, 5092−5106. (66) Zuev, D.; Bravaya, K. B.; Crawford, T. D.; Lindh, R.; Krylov, A. I. Electronic structure of the two isomers of the anionic form of pcoumaric acid chromophore. J. Chem. Phys. 2011, 134, 034310. (67) Silva-Junior, M. R.; Schreiber, M.; Sauer, S. P. A.; Thiel, W. Benchmarks of electronically excited states: Basis set effects on CASPT2 results. J. Chem. Phys. 2010, 133, 174318. (68) Ben-Nun, M.; Martinez, T. J. Photodynamics of ethylene: ab initio studies of conical intersections. Chem. Phys. 2000, 259, 237−248. (69) Hazra, A.; Soudackov, A. V.; Hammes-Schiffer, S. Role of Solvent Dynamics in Ultrafast Photoinduced Proton-Coupled Electron Transfer Reactions in Solution. J. Phys. Chem. B 2010, 114, 12319− 12332. (70) Soudackov, A. V.; Hazra, A.; Hammes-Schiffer, S. Multidimensional treatment of stochastic solvent dynamics in photoinduced proton-coupled electron transfer processes: Sequential, concerted, and complex branching mechanisms. J. Chem. Phys. 2011, 135, 144115. (71) Auer, B.; Soudackov, A. V.; Hammes-Schiffer, S. Nonadiabatic Dynamics of Photoinduced Proton-Coupled Electron Transfer: Comparison of Explicit and Implicit Solvent Simulations. J. Phys. Chem. B 2012, 116, 7695−7708. (72) Hazra, A.; Soudackov, A. V.; Hammes-Schiffer, S. Isotope Effects on the Nonequilibrium Dynamics of Ultrafast Photoinduced Proton-Coupled Electron Transfer Reactions in Solution. J. Phys. Chem. Lett. 2011, 2, 36−40.

K

DOI: 10.1021/jp5126969 J. Phys. Chem. B XXXX, XXX, XXX−XXX