MM Excited-State Dynamics of Retinal Protonated Schiff

Article pubs.acs.org/JPCB

Direct QM/MM Excited-State Dynamics of Retinal Protonated Schiff Base in Isolation and Methanol Solution Chutintorn Punwong,*,† Jane Owens,#,⊥ and Todd J. Martínez*,§,¶ †

Department of Physics, Faculty of Science, and Trace Analysis and Biosensor Research Center, Prince of Songkla University, Songkhla 90112, Thailand # Department of Chemistry, University of Illinois, Urbana-Champaign, Illinois 61801, United States § Department of Chemistry and The PULSE Institute, Stanford University, Stanford, California 94305, United States ¶ SLAC National Accelerator Laboratory, 2575 Sand Hill Road, Menlo Park, California 94025, United States S Supporting Information *

ABSTRACT: We use the full multiple spawning (FMS) dynamics approach with a hybrid quantum mechanics/molecular mechanics (QM/MM) reparameterized semiempirical method to investigate the excited-state dynamics of retinal protonated Schiff base (RPSB) in isolation, in neat methanol solution, and in methanol solution with a Cl− counterion. The excited-state lifetime is dramatically affected by MeOH solvent, which slows down the photoisomerization by an order of magnitude. We show that this is due to charge migration in the RPSB chromophore and the concomitant solvent friction in polar media. Simulation results are compared to experiments where available, with good agreement for excited-state lifetimes, bond selectivity of isomerization, and the time/energy-resolved fluorescence spectrum. We find that the inclusion of a Cl− counterion in the simulations has little effect on lifetimes, mechanism, or bond selectivity. In contrast to previous studies limited to RPSB and a surrounding counterion, we find that the placement of the counterion has little effect on bond selectivity. This suggests that dielectric screening can spoil the effect of a counterion in directing excited-state reactivity.

■

INTRODUCTION Retinal protonated Schiff base (RPSB), the chromophore of the rhodopsin family of proteins, is well-known for its ability to convert light into mechanical energy. Upon photoexcitation, RPSB isomerizes at one of the double bonds in the conjugated chain. In protein environments, the isomerization is highly selective. For example, RPSB undergoes isomerization only at C11C12 in rhodopsin (Rh) and at C13C14 in bacteriorhodopsin (bR) and halorhodopsin (hR). On the other hand, in solvated environments such as methanol, several photoisomerization products are possible, including isomerization about the C9C10, C11C12, and C13C14 bonds.1,2 In addition to the bond selectivity of isomerization in the protein environment, the quantum yield of isomerization in protein is also enhanced, 67 and 65% in Rh and bR, respectively, compared to only 10−30% in methanol.1−6 Moreover, the time scale of 2−10 ps4,5,7,8 for isomerization in methanol is much slower than that in proteins, where isomerization occurs within hundreds of femtoseconds.9−11 The protein environment increases the selectivity, quantum yield, and rate of photoinduced isomerization in RPSB. The origin of this enhancement presumably lies in the specific arrangements of key residues around the chromophore. A thorough understanding of this enhancement is best obtained by comparison of the photodynamics of RBSB in isolated, solvated, and protein environ© 2014 American Chemical Society

ments. Here, we focus on the photodynamics in a solvated (methanol) environment, with comparison to the predicted dynamics of the isolated chromophore. Condensed-phase simulation of the nonadiabatic dynamics is carried out using the full multiple spawning (FMS) method with a hybrid quantum mechanics/molecular mechanics (QM/MM) scheme for the potential energy surfaces and their couplings. Direct comparison of the excited-state dynamics in vacuum and neat solvent environments provides insight into the influence of solvent on the photochemical mechanism and bond selectivity of RPSB isomerization. We also compare to simulations including a Cl− ion in order to determine its effects on photoisomerization. Most previous related work simulating photoisomerization in solvent has been restricted to small RPSB analogues, usually the three double-bond analogue known as PSB3.12−16 Studies of solvent effects for large analogues containing all isomerizable double bonds have typically been restricted to the Franck− Condon region or constrained geometries.17,18 We note that a recent study did follow the nonadiabatic dynamics of a large Special Issue: William L. Jorgensen Festschrift Received: April 21, 2014 Revised: August 29, 2014 Published: September 2, 2014 704

dx.doi.org/10.1021/jp5038798 | J. Phys. Chem. B 2015, 119, 704−714

The Journal of Physical Chemistry B

Article

Figure 1. Structure and numbering system of all-trans RPSB (left) and the truncated analogue (right) studied in this paper.

molecules with less than 20 atoms due to computational expense.12,30−32 The QM region in our simulations includes all isomerizable double bonds in the RPSB chromophore. This RPSB analogue includes 31 atoms and excludes the β-ionone ring and butylamine linker at the two ends of the full RPSB chromophore, as depicted in Figure 1. The broadening parameter for the FOMO-SCF procedure is set to 0.2 atomic units, and the broadening is restricted to six active orbitals. The FOMO-CASCI calculations are performed with an active space of six electrons in six orbitals, denoted as CASCI(6/6). We have reoptimized the semiempirical parameters (based on the PM3 Hamiltonian and parameterization33) for RPSB by comparison to ab initio calculations, following the strategy discussed in our previous work.21,34 As shown in Figures S2−S4 of the Supporting Information (through comparison to ab initio CAS(10/10)-PT2 calculationd), the new parameter set provides a good description of the ground- and excited-state potential energy surfaces for the RPSB analogue. This validates both the reparameterization and the choice of active space. The MM region (85 methanol solvent molecules) is described using the OPLS-AA force field. As usual, the MM Hamiltonian, Ĥ MM, contains bonded and nonbonded interactions, with the general form

analogue of RPSB in MeOH solution using a QM/MM method but was restricted to few initial conditions and ∼2 ps because of computational limitations.19 The simulations reported here use a large analogue of RPSB including all isomerizable double bonds, follow the nonadiabatic dynamics for up to 10 ps, and encompass up to 100 trajectories (including spawning) for each of the initial conditions examined.

■

THEORY AND METHODS 1. Hybrid QM/MM Method. Realistic treatment of the potential energy surfaces and nonadiabatic couplings involved in excited-state dynamics requires a QM model for the chromophore. In order to describe conical intersections correctly, the QM method should allow for multireference character in the electronic wave function (because two or more electronic states will be exactly degenerate at a conical intersection). Because the electronic excitation is normally localized to the chromophore, it may be sufficient to model the surrounding environment using a MM force field. Coupling the QM treatment of the chromophore with a MM description of the surroundings gives rise to the QM/MM method, first developed by Warshel and Levitt.20 In this method, the Hamiltonian is partitioned as Ĥ total = Ĥ QM + ĤMM + Ĥ QM/MM

(1)

EMM = E bonded + Enonbonded

where Ĥ QM is the Hamiltonian of the QM part, Ĥ MM is the Hamiltonian of the MM part, and Ĥ QM/MM is the coupling between QM and MM regions. In this work, we use a multireference reparameterized semiempirical method21,22 to describe the QM region (Ĥ QM). The self-consistent field problem is solved using a floating occupation molecular orbital (FOMO) method,23,24 which can be viewed as a finite temperature ensemble technique. This ensures that degenerate molecular orbitals are assigned the same (possibly noninteger) electronic occupation. After the molecular orbitals have been determined, a full configuration interaction wave function is formed for an active space with a specified number of electrons and orbitals. We refer to this as FOMO-CASCI(ne/no) for an active space consisting of ne electrons in no orbitals. The FOMO-CASCI method has been implemented in a modified version22,23,25,26 of MOPAC2000.27 Ab initio implementation of FOMO-CASCI has been shown24 to provide similar results to the ab initio state-averaged complete active space (SA-CASSCF) method,28,29 but in the present work, a semiempirical representation of the electronic integrals is used. Thus, some treatment of dynamic electron correlation is implicit in the method, and one could hope for accuracy intermediate between CASSCF and CASPT2. Its computational cost is significantly less than ab initio methods such as SA-CASSCF because of the simplifications from the semiempirical approximations. This is critical for the present work because excited-state dynamics calculations with ab initio methods including both multireference character and dynamic electron correlation (e.g., CASPT2) are generally limited to

= [E bond + Eangle + Edihedral] + [E ES + EvdW ] ⎡ 1 =⎢ ⎢⎣ 2

∑

kr(r − req)2 +

bonds

1 2

∑

kθ(θ − θeq)2

angles

⎤ Vn (1 − ( −1)n cos(nϕ)⎥ ⎥⎦ dihedrals n = 1 2 ⎡qq ⎛ σ 12 σij6 ⎞⎤ i j ij ⎢ ⎜ +∑ + 4εij⎜ 12 − 6 ⎟⎟⎥fij ⎢ R ij ⎠⎥⎦ ⎝ R ij i < j ⎣ R ij 3

+

∑ ∑

(2)

Force constants, equilibrium bond lengths and angles, and torsional potential parameters are denoted as kr/kθ, req/θeq, and Vn, respectively. The charge (fixed throughout the simulations) on the ith atom is given as qi, and the van der Waals well depth and separation for the ij atomic pair are given as εij and σij, respectively. The nonbonded terms include electrostatic and van der Waals interactions, and these are excluded for atoms that are within three atoms of each other (as judged by bond and angle terms in the potential). As usual, in the OPLS/AA force field,35 1−4 nonbonded interactions are scaled, that is, f ij is unity except for 1−4 interactions where f ij = 0.5. The coupling between the QM and MM regions, Ĥ QM/MM, consists of electrostatic and van der Waals interactions ES vdW Ĥ QM/MM = Ĥ QM/MM + Ĥ QM/MM

705

(3)



Article

Figure 2. RPSB analogue (in blue stick representation) solvated by methanol molecules. The RBSB analogue is modeled with a reparameterized multireference semiempirical method, and the surrounding methanol molecules are modeled using the OPLS/AA force field.

methanol molecules are explicitly represented using the OPLS/ AA force field35 implemented in a development version of MOPAC2000.27 After random placement of the methanol molecules around the chromophore, the coordinates of all atoms (QM and MM) are relaxed using the QM/MM potential energy function. This is a microsolvated cluster, that is, we do not impose periodic boundary conditions. All atoms in the simulation (85 MM MeOH molecules and the QM RPSB analogue) are allowed to move without constraints during both optimization and dynamics. There is negligible (at most, three MeOH molecules over all of the simulations reported here) evaporation over the course of the simulations, which is partially due to the fact that the excited-state lifetime is quite short and only up to 10 ps of simulation time is required. 2. FMS Dynamics. We model the nonadiabatic dynamics after photoexcitation using the FMS method, which can describe the breakdown of the Born−Oppenheimer approximation36 and population transfer between electronic states, for example, near-conical intersections. The FMS37−40 method is designed to treat problems involving dynamics on multiple electronic states by employing the quantum mechanical nonadiabatic couplings that trigger the transfer of population between electronic states while maintaining some connection to classical dynamics in the description of the nuclear wave function. Specifically, the nuclear wave function is expanded in a basis set of complex Gaussian functions, called trajectory basis functions (TBFs). The phase space center of each of these TBFs evolves according to classical mechanics. Potential energy surfaces (PESs), governing the propagation of the nuclear basis function, and nonadiabatic coupling are calculated every time step “on-the-fly” using the QM/MM method detailed above. During nonadiabatic events, the nuclear basis set expands as needed (“spawning”), and throughout the dynamics, the amplitudes of the nuclear basis functions are determined by solving the nuclear Schrödinger equation in the time-evolving basis set. The FMS method has been discussed in detail previously,37−40 and the reader is referred to previous work for more information. The initial conditions (center of position and momentum) for each independent TBF (after relaxation of the solvated

The electrostatic interaction is given by ES Ĥ QM/MM =

−qm

∑ i ∈ QM electrons m ∈ MM atoms

|ri ̂ − R m|

+

∑ j ∈ QM nuclei m ∈ MM atoms

qjqm R jm (4)

where the first term includes the Coulomb interactions between QM electrons and point charges centered on MM atoms and the second term represents the interactions between nuclear core charges of the QM atoms and the MM atomic point charges. Note that the nuclear charge of the QM atoms generally is not the same as the atomic number because the innermost core electrons (e.g., the 1s electrons of carbon) are often not modeled explicitly in semiempirical methods. Thus, the nuclear charge of the H atom (which has no innermost core electrons) is +1, and that of the C atom is +4 (both values given in atomic units, which are used throughout unless otherwise specified). The distance between the QM nucleus j and MM atom m is denoted as Rjm. The first term of QM/MM electrostatic interactions directly involves electronic coordinates and therefore requires the electronic structure calculations, FOMO-CASCI,22,23,25,26 as described above. The van der Waals interaction is included to introduce quantum repulsion and dispersion and is given as vdW Ĥ QM/MM =

∑ j ∈ QM atoms m ∈ MM atoms

6 ⎤ ⎡ σ 12 σjm jm 4εjm⎢ 12 − 6 ⎥ ⎢⎣ R jm R jm ⎥⎦

(5)

where εjm and σjm are Lennard-Jones parameters for the interaction between atoms j and m. In this case, the QM atoms need to be specified according to the corresponding MM force field so that the Lennard-Jones parameters (εj and σj) can be defined. These two parameters are then used to determine εjm and σjm by the Lennard-Jones combination rules for the vdW interactions of intersystem QM/MM coupling (we use geometric mean combination rules for both of these, as is also done in the OPLS-AA force field). The QM/MM system is prepared by solvation of the QM RPSB analogue with 85 MM methanol molecules in an approximately 6000 Å3 box, as shown in Figure 2. The MM 706



Article

Table 1. Averaged Results from FMS Dynamics of Photoexcited RPSB in Isolation (Gas), Neat Solution (MeOH), and Solution with a Counterion (MeOH/Cl−)a number of TBFs

fraction cis-S0

isomerization time (ps)

bond

gas

MeOH

MeOH/Cl−

C7C8 C9C10 C11C12 C13C14 C15Nb total expt1−8

0 8 17 5 0 30

0 11 22 3 4 40

0 7 17 7 9 40

gas − 0.27 0.30 0.27 − 0.28

± 0.04 ± 0.04 ± 0.05 ± 0.03

MeOH − 2.5 ± 3.2 ± 4.1 ± 3.9 ± 3.2 ± 2−10

0.4 0.4 2.0 1.0 0.3

MeOH/Cl− − 3.7 ± 4.2 ± 5.4 ± 2.4 ± 3.9 ± 2−10

gas 0 0.52 0.31 0.40 0 0.38

0.5 0.7 0.9 0.4 0.4

± 0.1 ± 0.1 ± 0.2 ± 0.1

MeOH

MeOH/Cl−

0 0.43 ± 0.1 0.40 ± 0.1 0.41 ± 0.1 0.49 ± 0.2 0.37 ± 0.1 0.1−0.3

0 0.37 ± 0.1 0.47 ± 0.1 0.42 ± 0.1 0.51 ± 0.1 0.34 ± 0.1 0.1−0.3

a The column labeled “Number of TBFs” indicates the number of initial conditions that decay to the ground state by isomerization about a specific bond. The average fraction of the spawned TBFs that relax to a cis conformation after nonadiabatic relaxation due to isomerization about a specific bond is indicated in the column labeled “Fraction cis-S0”. Error bars are determined by bootstrapping.68 bThe C15N isomerized product is experimentally indistinguishable from reactant. Thus, this is not included in the reported isomerization quantum yield.

chromophore, as discussed above) are equilibrated on the ground state (S0) for 1 ps using Brownian motion at 298 K (with a friction coefficient of 0.1 fs−1). After this equilibration, the TBF is placed on the excited state (S1), and FMS dynamics begins with a single Gaussian wavepacket. A fourth-order Runga−Kutta scheme41 is used with a time step of 0.5 fs to integrate the equations of motion for the corresponding classical trajectories and complex amplitudes of the wave functions. This procedure is repeated with 30/40 distinct initial conditions for simulations in vacuum/methanol, respectively. When the RPSB chromophore is solvated by MeOH, the S1 and S2 states become energetically quite close in the FOMOCASCI calculations (in contrast to the isolated case, where these are well-separated). This has also been observed previously with ab initio methods, both in MeOH17,18 and also when surrounded by a Cl− ion.42 We ran test calculations for the solvated case including S2, and these showed very little (

MM Excited-State Dynamics of Retinal Protonated Schiff

Recommend Documents