J. Phys. Chem. A 2010, 114, 6757–6765
6757
Nonadiabatic Excited-State Dynamics with Hybrid ab Initio Quantum-Mechanical/ Molecular-Mechanical Methods: Solvation of the Pentadieniminium Cation in Apolar Media Matthias Ruckenbauer,† Mario Barbatti,† Thomas Mu¨ller,‡ and Hans Lischka*,†,§ Institute for Theoretical Chemistry, UniVersity of Vienna, Waehringerstrasse 17, A-1090 Vienna, Austria, Institute of AdVanced Simulation, Ju¨lich Supercomputer Centre, Research Centre Ju¨lich, D-52425 Ju¨lich, Germany, and Institute of Organic Chemistry and Biochemistry, Academy of Sciences of the Czech Republic, FlemingoVo nam. 2, CZ-16610 Prague 6, Czech Republic ReceiVed: April 6, 2010; ReVised Manuscript ReceiVed: May 7, 2010
A new implementation of nonadiabatic excited-state dynamics using hybrid methods is presented. The current approach is aimed at the simulation of photoexcited molecules in solution. The chromophore is treated at the ab initio level, and its interaction with the solvent is approximated by point charges within the electrostatic embedding approach and by a Lennard-Jones potential for the nonbonded interactions. Multireference configuration interaction (MRCI) and multiconfiguration self-consistent field (MCSCF) methods can be used. The program implementation has been performed on the basis of the Columbus and Newton-X program systems. For example, the dynamics of penta-2,4-dien-1-iminium (PSB3) and 4-methyl-penta-2,4-dien-1iminium cations (MePSB3) was investigated in gas phase and in n-hexane solution. The excited-state (S1) lifetime and temporal evolution of geometrical parameters were computed. In the case of PSB3 the n-hexane results resemble closely the gas phase data. MePSB3, however, shows a distinct extension of lifetime due to steric hindering of the torsion around the central bond because of solute-solvent interactions. 1. Introduction Processes occurring in the electronic excited states of molecules are of fundamental importance for the underlying phenomena in photochemistry and photobiology such as photostability of nucleic acids,1 photoaging,2 UV-induced immunosuppression,3 and the molecular mechanism of vision.4 In combination with experimental ultrafast (femtosecond) investiagtions,5 dynamics simulations on multiple potential energy surfaces became a valuable tool for understanding such complex phenomena in real time at the molecular level.6-11 Besides the intrinsic difficulties involved in the proper description of the excited-state electronic structure required in the dynamics simulations, the theoretical investigations of these phenomena should also take into account that the overwhelming majority of the processes mentioned above occur in condensed matter such as solution and/or biological environments. The environment naturally affects the photoexcited molecule promoting steric hindering, new channels for charge and energy flow, and energetic shifts that can substantially change the topography of the potential energy surface with respect to the gas phase.12 Despite computational advances, the purely quantum-mechanical computation of extended molecular systems, including the environment, is still not feasible even in the ground state, not to mention in excited states. The main strategy developed for treating such cases is the use of hybrid methods, where each subsystem is treated by a different computational level, choosing the accuracy/cost of each one according to the expected importance of the region. Probably the most popular hybrid method is the QM/MM approach, in which a reactive primary * To whom correspondence should be addressed. E-mail: hans.lischka@ univie.ac.at. † University of Vienna. ‡ Research Centre Ju¨lich. § Academy of Sciences of the Czech Republic.
system is treated quantum mechanically (QM), while the remaining (secondary) system is treated by molecular mechanics (MM). For a review on the large literature on QM/MM approaches see for example, refs 13-15. There are several possibilities to implement a QM/MM approach. For the present application of embedding into a solvent, the main two strategies followed are mechanical embedding and electrostatic embedding, respectively.15 While the mechanical embedding takes into account the interaction between the primary and secondary systems at MM level only, the electrostatic embedding includes the electrostatic part of this interaction at QM level by incorporating additional one-electron terms into the electronic Hamiltonian to represent the secondary system. The current implementation is designed to be flexible enough to allow both strategies, letting the final decision to the user of how to proceed for each specific application. In the case of applying QM/MM methods to photochemical investigations, the electrostatic embedding is certainly to be preferred in view of the difficulties of describing the charge distribution of an electronically excited system within the MM approach. Additionally, the electrostatic embedding has the advantage of treating the polarization of the primary system by the secondary system naturally. A number of recent investigations of photochemical nonadiabatic phenomena have been reported using the QM/MM approach and surface hopping.16-24 Many of them, however, treat the nonadiabatic event itself not in full generality since either the dynamics simulation is terminated as soon as the region of the conical intersection is reached or simplified surface hopping algorithms are employed. The inclusion of the solvent may substantially increase the complexity of the dynamics close to the crossing seam, an effect that would probably be hidden by approximate hopping algorithms. In the present implementation, the nonadiabatic coupling terms are fully computed by including the contributions from the secondary system through
10.1021/jp103101t 2010 American Chemical Society Published on Web 06/02/2010
6758
J. Phys. Chem. A, Vol. 114, No. 25, 2010
Ruckenbauer et al. the procedures for computing the analytic multireference configuration interaction (MRCI) gradients and nonadiabatic coupling vectors within electrostatic embedding, and the generation of appropriate initial conditions for the dynamics. The present purpose is to observe the effect of spatial restrictions within a solvent cage on the nonadiabatic deactivation process. The relatively small PSB3 shows, as will be discussed, no big effect in this respect. Therefore, simulations with the methylsubstituted derivative MePSB3 (3-methyl-penta-2,4-diene-1iminium cation) have been performed. The CH3-group reaching out from the molecule acts as a sterical “anchor” in the solvent, hindering the torsion of the molecule.
Figure 1. Structure of trans-PSB3 and trans-4-methyl-PSB3 in gas phase and solution. Only the first solvation shell is shown.
the one-electron terms in electronic Hamiltonian. Thus, the nonadiabatic events can be computed by means of the welltested fewest switches algorithm by Tully,25 and the contribution of the environment to the nonadiabatic events can be directly evaluated. The photophysics and photochemistry of protonated Schiff bases (CH2[CH]2n-2NH2+), here denoted as PSBn, has been subject of great interest because these systems are important models for the chromophore of the class of rhodopsins6,26-36 and because they range among computationally feasible polyatomic systems presenting nonadiabatic ultrafast behavior. In particular, the photodynamics of the PSB3 species (penta-2,4dien-1-iminium cation, see Figure 1) after UV excitation into the lowest singlet ππ* state is an interesting example, being characterized by skeletal relaxation occurring in the first few tens of femtoseconds followed by torsional isomerization around the central bond in the next hundred femtoseconds.6,34,36,37 The geometrical deformations triggered by UV excitation of PSB3, which naturally correspond to the stabilization of the excited state, also lead to destabilization of the ground state and promote the formation of a seam of conical intersections between these two states. In the proximity of this intersection seam, PSB3 has a large probability of returning to the ground state by radiationless deactivation. In spite of the fact that the conical intersection of lowest energy is located at 90° central-bond torsional angle, both minimum energy path27 and dynamics analysis9,36 show that in gas phase PSB3 returns to the ground state at conical intersections located at about 70° when the cis isomer is excited and at about 120° when the trans isomer is excited. Additionally, dynamics simulations in gas phase predict an excited-state lifetime of about 100 fs.36 Many theoretical studies have dealt with the excited states of the PSBn and related systems in other media than vacuum; most of them, however, focused on Franck-Condon properties or on description of potential energy surfaces.26,28,38-43 This work presents systematic steps for a comprehensive QM/MM approach for nonadiabatic dynamics in solution based on the program systems Columbus,44-46 Tinker,47 and Newton-X.35,48 It outlines the structure of the QM/MM program setup itself,
2. QM/MM Implementation 2.1. General Concept. The atoms of the entire system S to be treated by means of the hybrid method are divided into disjoint regions. Sets of these regions can be assembled to be treated by different methods, computing the energy, gradient, and, if possible, the nonadiabatic coupling vector. The results of these separate jobs are combined to a total result as defined by the user. At the present stage the program development is dedicated to the computation of solvent effects and of other nonbonded interactions on a single chromophore. Thus, there is no need for a special treatment of a boundary region between the different regions of atoms. For a standard QM/MM setup with electrostatic embedding these subsets are typically an inner (I) and an outer region (O). Inner and outer region are described by QM and MM, respectively. Specifically, multireference electronic structure methods are used to accurately describe multiple electronic states of the compound of interest, whereas the MM component primarily deals with secondary environmental effects. Standard parametrized force fields are employed in the MM part incorporating bonded terms (bond stretching, angle bending, proper and improper torsions), van-der-Waals interactions (Lennard-Jones type potential), and electrostatic interaction between partial point charges associated with each atom. The total energy of the entire system (S) is given by b+vdW el EQM/MM(S) ) EMM (O) + EMM (O) + Eel(I, O) + vdW (I, O) + EQM(I) (1) EMM
where the superscripts denote bonding (b), van-der-Waals (vdW) and electrostatic (el) interactions. The coupling between both vdW (I, O) comregions arises from the van-der-Waals terms EMM puted at MM level and the electrostatic coupling term Eel(I, O). In case of electrostatic embedding, the latter term is computed by including the electrostatic potential of the solvent into the Hamiltonian of the QM region via point charges as follows. ˆ (in atomic units), The standard quantum chemical Hamiltonian H given by
ˆ ) Tˆe + Vˆne + Vˆee + Vnn ) H )
∑ hˆ(i) + ∑ |rF i
i
