J. Phys. Chem. A 2010, 114, 12585–12590
12585
Azomethane: Nonadiabatic Photodynamical Simulations in Solution Matthias Ruckenbauer,† Mario Barbatti,‡ Bernhard Sellner,† Thomas Muller,§ and Hans Lischka*,† Institute for Theoretical Chemistry, UniVersity of Vienna, Waehringerstrasse 17, 1090 Vienna, Austria, Max Plank Institute fuer Kohlenforschung, Kaiser-Wilhelm-Platz 1, 45470 Muelheim, Germany, and Institute for AdVanced Simulation, Research Center Ju¨lich, 53425 Juelich, Germany ReceiVed: September 16, 2010; ReVised Manuscript ReceiVed: October 26, 2010
The nonadiabatic deactivation of trans-azomethane starting from the nπ* state has been investigated in gas phase, water, and n-hexane using an on-the-fly surface-hopping method. A quantum mechanical/molecular mechanics (QM/MM) approach was used employing a flexible quantum chemical level for the description of electronically excited states and bond dissociation (generalized valence bond perfect-pairing complete active space). The solvent effect on the lifetime and structural parameters of azomethane was investigated in detail. The calculations show that the nonadiabatic deactivation is characterized by a CNNC torsion, mainly impeded by mechanic interaction with the solvent molecules. The similar characteristics of the dynamics in polar and nonpolar solvent indicate that solvent effects based on electrostatic interactions do not play a major role. Lifetimes increase by about 20 fs for both solvents with respect to the 113 fs found for the gas phase. The present subpicosecond dynamics also nicely show an example of the suppression of C-N dissociation by the solvent cage. Introduction Symmetric azoalkanes are known to dissociate after irradiation with ultraviolet light or after thermal activation and are commonly used as a source for clean alkyl radicals. Their reactivity has been investigated experimentally and theoretically already for a long time.1-5 Theoretical investigations have concentrated on the simplest azoalkane, azomethane, but also aromatic azo compounds have been studied intensively.6,7 Mostly, the photochemical properties of the first excited singlet state, which is a dark nπ* state, have been considered. This state can be populated by photonic excitation in the 350-300 nm range (see, e.g., refs 4 and 8-10) via vibronic coupling. It is generally agreed that this state is deactivated in an ultrafast process through a conical intersection between the S1 and S0 states. The main geometrical distortion leading from the trans isomer to this intersection is a torsion around the NdN bond in the excited state, which results in an isomerization with a high quantum yield. More details can be found in refs 4, 5, and 9. Experimental and theoretical investigations show that in the gas phase the time for the deactivation of the excited state to the ground state is in the subpicosecond range. Two basic models have been proposed for the subsequent CN dissociation process: the impulsive variant effecting a fast high-energy dissociation of the first CH3 group immediately after the transition to the ground state and the statistic process which involves an intermediate hot ground state motion and subsequent dissociation due to the large amount of kinetic energy available to the molecule from the electronic excitation. Most studies favor the statistical over the impulsive model,8,11 and recent ab initio dynamics simulations of gas-phase azomethane5 have showed fast dissociations to be a very rare event in the subpicosecond range. * To whom correspondence should be addressed. E-mail: hans.lischka@ univie.ac.at. † University of Vienna. ‡ Max Plank Institute fuer Kohlenforschung. § Research Center Ju¨lich.
In solvents, the CN dissociation is highly inhibited, and isomerized azomethane is the main product after irradiation.12,13 The suppression of dissociation seems to be widely invariant with the nature of the solvent. Cattaneo and Persico studied this effect theoretically for aqueous solvation based on semiclassical nonadiabatic dynamics on energy surfaces fitted to ab initio calculations.3,14,15 They found that the short time dynamics and the isomerization are only marginally influenced by polar solvation, and the main effect of the solvent lies in the suppression of the dissociation through a cage effect. In this paper, a comparative study of the nonadiabatic shorttime photodynamics of trans-azomethane in gas phase, polar (aqueous), and nonpolar (n-hexane) solvation is presented with the goal of studying in detail the environmental effects on the photodeactivation mechanism. A recently developed quantum mechanical/molecular mechanics (QM/MM) procedure combined with multireference configuration interaction (MRCI) methods is used.16 To the best of our knowledge, this is the first attempt to simulate the photoreaction of azoalkanes in solution using high-level multireference ab initio QM/MM onthe-fly methods with treatment of all degrees of freedom and the nonadiabatic effects in full generality. For a detailed understanding of the possible pathways leading to the nonadiabatic deactivation and the solvent’s influence on them, the employment of such a high level is of vital importance. Computational Details QM/MM Ansatz. A QM/MM approach was used17-19 separating the complete system of azomethane and solvent into two subsets of atoms: an inner (I) and an outer (O) region. The inner and outer regions are described by quantum mechanics and molecular mechanics, respectively. Specifically, multireference electronic structure methods are used to accurately describe multiple electronic states of the compound of interest, whereas the MM component deals with environmental effects. Standard parametrized force fields are employed in the MM part incorporating bonded terms (bond stretching, angle bending,
10.1021/jp108844g 2010 American Chemical Society Published on Web 11/11/2010
12586
J. Phys. Chem. A, Vol. 114, No. 48, 2010
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 vdW EQM/MM(S) ) EQM(I) + Eel(I, O) + EMM (I, O) + b+vdW el (O) + EMM (O) (1) EMM
where the superscripts denote bonding (b), van der Waals (vdW) and electrostatic (el) interactions. An electrostatic embedding scheme is used in which the effective point charges of the atoms of the solvent molecules are included in the quantum mechanical Hamiltonian. More details of the specific implementation can be found in ref 16. QM and MM Levels. All ab initio calculations were performed at a state-averaged multiconfiguration self-consistent field (SA-MCSCF) level using a complete active space (CAS) combined as direct product with a generalized valence bond-perfect pairing approach.20,21 The CAS consisted of six electrons distributed into four orbitals (2 lone pairs (n), π, π*) and the σ-orbitals and their corresponding σ*-orbitals were kept in nine pairs (σNN, 2 × σCN, 6 × σCH) wherein each orbital pair two electrons are localized. State-averaging over two states was performed, and the 6-31G(d) basis set22 was used. The carbon and nitrogen inner shell orbitals were kept doubly occupied. The overall quantum mechanical level is denoted as SA-2-GVBPP(9)-CAS(6,4)/6-31G(d) or short as GVB-CAS. By comparison with extended MRCI calculations, this level was proven to be well suited for the accurate description of azomethane’s energy surfaces as well as for use in the molecular dynamics.5 For the molecular mechanics part of the calculations the van der Waals parameters, intramolecular parameters and effective charges of the solvent molecules were taken from the OPLS/ AA force field.23 For water, the TIP3P water model was chosen. For the quantum mechanical part of the calculation, the program system COLUMBUS24-26 was used. It provides the required analytic gradients and nonadiabatic coupling vectors.27-31 The implementation of the electrostatic embedding in COLUMBUS is described in ref 16. The molecular mechanics calculations were performed using TINKER.32 The combination of the hybrid energies and gradients, the integration of the equations of motion and the time-dependent Schroedinger equation, and the surface hopping were performed using NEWTON-X.33,34 Initial Conditions. For small- to medium-sized molecules in the gas phase, a Wigner distribution for the quantum harmonic oscillator in the ground state can be used for creation of initial geometries and velocities (see, e.g., ref 35). This is also the method used for the creation of the gas phase initial conditions and for the initial conditions of the solvated azomethane molecule. In earlier simulations,5 the assignment of additional quanta to the torsional mode in the initial conditions was used to mimic the torsional bias in previous femtosecond pump-probe experiments.4 This procedure also allows, without qualitatively changing the behavior of the dynamics, the observation of at least a few dissociations within a rather short simulation time of maximum 350 fs after the transition to the ground state. We made use of this advantage by assigning five additional quanta to the central torsional mode, as described in ref 5. In solution. the large number of solvent atoms to be considered requires a different scheme for obtaining initial structures and velocities. A commonly used alternative is to perform a thermalized ground state dynamics and select uncorrelated geometry/velocity points and use those as initial
Ruckenbauer et al. conditions. This procedure is well-suited to sample the intermolecular modes of the solvent. However, it has an inherent defect when sampling the modes of the solute: because of the classical motion, at room temperature, the thermostat will not provide enough kinetic energy for the higher vibrational frequencies of the molecule(s); in fact, it will for most cases not even provide the zero-point energy of the system. Because nonadiabatic dynamics can be sensitive to the distribution of initial conditions, it seems potentially unfavorable to apply the classical procedure for the inner region where the electronic excitation occurs. To overcome this problem, a scheme to create mixed initial conditions using a gas phase Wigner distribution for the geometries and velocities of the solute and an equilibrated and thermalized solvent structure around these geometries was developed. A spherical cluster of solvent molecules with constant volume was chosen as the environment (see below for a description of the setup). This scheme is characterized by the following points: (1) The solvent of the entire cluster was equilibrated and thermalized at the MM level keeping the molecule of the QMregion at the frozen ground state gas phase equilibrium structure. (2) After thermalization, statistically independent sample geometries of the solvent were selected from a thermalized ground state trajectory in steps of 1 ps with the QM-region still frozen at the equilibrium structure. (3) For the isolated molecule of the QM region, initial conditions were calculated using a Wigner distribution of a quantum harmonic oscillator. (4) The equilibrium QM structure embedded in each of the individual solvent clusters was replaced by one of the displaced structures of the Wigner distribution. (5) The MM region of each of these samples structures was then thermalized around its frozen, displaced QM-region structure so as to adapt the cavity to the slightly changed new structure. The ground state minimum structure of azomethane was optimized at the same level of theory as used in the dynamics. The initial thermalization was performed for 100 ps at 298 K using a Nose´-Hoover thermostat36 and TINKER for the dynamics using the OPLS/AA force field. Effective charges for azomethane in the molecular mechanics were ChelpG-charges37 obtained from a calculation at the same level as the nonadiabatic dynamics. The equations of motion were integrated using the Beeman algorithm38 with a step size of 1 fs. After thermalization, this trajectory was continued to select the solvent structures. The force constant matrix for the Wigner distribution was calculated using COLUMBUS at the same level as the nonadiabatic dynamics. The relaxation of the solvent cavity was performed for 20 ps for water and 10 ps for n-hexane with a 1 fs step size using NEWTON-X. An Anderson thermostat39,40 is used in NEWTON-X (T ) 298 K; collision frequency, 0.2 fs-1). In this way, the geometries and velocities of the quantum mechanical core were obtained with respect to the quantum mechanical nature of the description of this region while still getting suitable solvent structures with a thermalized distribution of geometries and velocities. Figure 1 displays examples of initial conditions of azomethane in gas phase, water, and n-hexane. Setup of the Solvent Cluster. The azomethane solute was included in a spherical cluster of 150 n-hexane and 300 water molecules, respectively. The initial packing was performed using the PACKMOL41 program. The solvent structure was then thermalized and relaxed around individual solute structures
Azomethane
J. Phys. Chem. A, Vol. 114, No. 48, 2010 12587 TABLE 1: Lifetimes, Average Central Torsional Angles at the Point of First Hopping and Average Vertical Excitation Energies at the Initial Geometries of the Dynamics
environment
t1 (fs)
t2 (fs)
τ (fs)
hopping angle, 〈θj 〉 ( σ (°)
gas phase n-hexane aqua
46 44 40
67 95 97
113 139 137
90 ( 9 90 ( 7 90 ( 8
time constants
initial condition av, Eex,vert (eV) 3.71 3.67 3.73
a consequence of the hopping. For the determination of the possibility of a back-hopping, only the kinetic energy of the core atoms was taken in consideration. The temperature was kept constant at 298 K using an Anderson thermostat. To not interfere with the nonadiabatic treatment, the action of the thermostat was restricted to the solvent molecules. One hundred trajectories each have been computed for both solvated systems and for the gas phase, 300 trajectories in total. Results and Discussion
Figure 1. Selected structures of azomethane in (a) gas phase, (b) n-hexane, and (c) water with atomic numbering given in the gas phase structure.
obtained from a Wigner distribution with procedures described in the section Initial Conditions. To keep the gross density at a given value during the dynamics, the entire cluster was included in a spherical, rigid boundary for the whole time of the simulation (elastic collision boundary). The volume of this sphere was determined by the condition of maintaining the experimental density of 0.651 g/L for n-hexane42 and 0.997 g/L for water.42 Dynamics Details. Mixed quantum-classical dynamics was performed with on-the-fly calculation of the electronic energies, energy gradients, and nonadiabatic couplings. The details of the method have been described elsewhere.33 In brief, the nuclear coordinates were integrated with a 0.5 fs time step while the time-dependent Schro¨dinger equation was integrated with a 0.01 fs time step using interpolated quantities. To reduce the computational demands, the partial coupling approximation43 was employed computing only the nonadiabatic couplings, including the current state. Decoherence effects were taken into account by means of the model presented in ref 44 (R ) 0.1 hartree). Transition probabilities to other surfaces were computed by means of the fewest-switches algorithm.45 The maximum simulation time was 700 fs. Trajectories staying in the ground state for more than 350 fs after a hopping were terminated. The nonadiabatic couplings were restricted to the atoms treated fully quantum mechanically (core atoms) by setting the nonadiabatic coupling vector components of the other atoms to zero. This procedure ensures that the nonadiabatic hopping probability is governed only by the quantum mechanically treated region. In addition, the excess kinetic energy at the time of hopping is distributed only to these core atoms, which prevents an unphysical drain of kinetic energy to the solvent as
Solvatochromic Effects. Static solvent effects (i.e., the effects on the potential energy surface of the solute itself) can be observed by a change in the vertical excitation energies of the solute. The average excitation energy sampled from the initial geometries of the dynamics is used to characterize the solvent shifts. The method of generation of initial conditions in this study (see above) ensures that the solute geometries of the initial conditions are exactly the same in gas phase and in solvent. A change in the average vertical excitation energy of the initial conditions is therefore related only to solvent effects on the electronic structure. Table 1 shows the average vertical excitation energies of the initial conditions of all three investigated systems. The maximum difference is 0.06 eV between aqueous (3.73 eV) and n-hexane (3.67 eV) solvation; the value for gas phase is between these two (3.71 eV). The vertical excitation energy of the relaxed ground state structure in gas phase computed at the same theoretical level is 3.87 eV. The fact that the polarity of the solvent does not influence the vertical excitation significantly is not surprising in view of the fact that azomethane is neither in the ground state nor in the excited state polar. Excited State Lifetimes. In Figure 2, the evolution of the excited-state population for all three investigated systems is displayed. The population shows a latency time up to which the S1 population remains practically constant at 100%. A time of ∼50 fs is needed to reach the vicinity of the crossing seam with hoppings to the ground state. From this point on, the decay of the excited state is exponential. To obtain the lifetime, the range from 50 to 700 fs of the decay curves of Figure 2 has been fitted to the formula
Figure 2. Time evolution of the average excited-state occupation of azomethane in gas phase, n-hexane, and aqueous solution.
12588
J. Phys. Chem. A, Vol. 114, No. 48, 2010
(
y ) y0 + (1 - y0) exp -
(t - t1) t2
Ruckenbauer et al.
)
where t1 is the initial delay or latency time and t2 is the decay constant. The baseline, y0, was always set to zero, reflecting the assumption of a purely exponential decay with complete depletion of the excited state in the long time limit. The excitedstate lifetime (τ) is the sum of t1 and t2. The latency time, t1, is the time needed for the first trajectories to reach the crossing seam and perform the transition to the ground state. Reaching the crossing seam already requires a change in the molecular geometry. A change in t1 due to solvation reflects the influence of the solvent on the least hindered trajectories. It therefore shows in the absence of static solvent effects the minimal mechanic restriction on the excited state motion. In contrast, the decay constant t2 reflects the reaction rate of the nonadiabatic S1 f S0 deactivation after the first trajectories reached the conic intersection. Under the assumption that the time spent in the vicinity of the crossing seam is short when compared with the time needed to reach it, this parameter is dominated by the average velocity of the geometrical changes in the excited state. A change in this parameter, again in the absence of static solvent effects, gives information about the average mechanic hindrance of the pathway to the crossing seam. Computed lifetimes and average hopping angles can be found in Table 1. Azomethane in the gas phase shows a latency time, t1, of 46 fs and a decay constant, t2, of 67 fs, summing up to a total lifetime for the excited state of 113 fs. Using no additional torsional quanta in the initial conditions, the lifetime rises to ∼190 fs.5 These decay times are much shorter than the previous results obtained by Cattaneo and Persico (∼400 fs)3 but are much closer to the experimentally found 70 fs reported by Zewail and Diau.4 Solvation neither in water nor in n-hexane changes the latency time beyond the statistical uncertainty. In aqueous solution, the initial delay is 40 fs; in n-hexane, 44 fs, as compared with 46 fs in gas phase. The effect on the decay constant is significant, elongating t2 from 67 fs in gas phase to 97 fs in water and to 95 fs in n-hexane, giving total lifetimes of 137 fs in water (+24 fs relative to gas phase) and 139 fs in n-hexane (+26 fs relative to gas phase). The—within the limits of the statistics—identical latency times t1 between gas phase and solution confirm that the initial hindrance caused by the solvent is not very strong. The fact that the elongation of the decay constant t2 is independent of the polarity and the nature of the solvent together with the fact that no strong solvatochromic shift was observed leads to the conclusion that neither the potential energy surfaces nor the nonadiabatic couplings are largely affected by static electronic effects. This leaves mechanic effects as the major source of solvent influence on the dynamics of azomethane. Geometrical Evolution. The main geometrical distortion leading to the crossing seam in azomethane is torsion around the central NdN bond, which is a double bond in the ground state. This torsion is not characterized by a rigid torsion with widely moving -CH3 groups around a fixed N-N axis but involves coupling with other coordinates, especially with the CNN angles. Figure 3 shows illustrations of the pathway from the Franck-Condon point to the point of hopping for representative trajectories in gas phase and in aqueous and n-hexane solution by plotting the overlying atomic spheres along the trajectory. One can observe a very efficient and space-saving movement already in the gas phase, which affords only small displacements of all atoms. In water and n-hexane, the trajec-
Figure 3. Pathway from the vertical excitation to the conical intersection of typical azomethane trajectories for gas-phase, water, and n-hexane solvation. The C-N-N-C backbone of the initial condition (green lines) and point of hopping (red lines) and the path of the atoms between these two times are displayed. The initial conditions of the azomethane moiety were identical in all three trajectories.
tories displayed in Figure 3 show that some atoms of the azomethane are inhibited in their motion. However, this hindrance can be easily compensated by only small changes in the movement of the other atoms. This is an extremely flexible and adaptive system that can easily adapt to environmental restrictions. Figure 4 shows the time evolution of the torsional angles for all trajectories in the gas phase, water, and n-hexane. The points of transition to the ground state are marked by black rhombs. Torsion in azomethane is equally distributed in both possible directions. Therefore, the central torsional angle is folded in the following into the range between 0° and 180°. In the gas phase, water, and n-hexane, the trajectories switched to the ground state at an average angle of 90°, with a standard deviation from 7° to 9°, close to the minimum of the crossing seam located at a 90° torsion of the central C-NdN-C dihedral angle.5 Figure 5 shows the time evolution of the average central torsional angle. During the first 40 fs (the delay time for the first trajectories to reach the crossing seam), the torsion is not affected at all by the solvent. After that, for both solvents, the central torsion does not approach the conical intersection at approximately 90° so quickly as in the gas phase. This is consistent with the observations regarding the lifetime of the excited state. The reduced torsion lets azomethane reach the crossing seam later, on average, due to the interaction with the solvent molecules. Another interesting question is the CN dissociation leading finally to N2 and two CH3• radicals. Most of these dissociations take place on a picosecond time scale3,5 but, even though they are a rare event at the subpicosecond scale, they can be observed in the present dynamics, as well. To characterize the single and sequential double dissociations, the absolute value of the difference of the two C-N bond lengths (absolute bond length difference, aBLD ) |C3N1-C4N2|) was used. Single dissociations are characterized by a continuous increase in the aBLD, whereas the sequential double dissociations can be identified
Azomethane
Figure 4. Evolution of the central torsional angle for all trajectories in the gas phase, water, and n-hexane. Black diamonds mark the points of transition to the ground state.
Figure 5. Evolution of the average central torsion of azomethane.
by an increase in the aBLD value with a subsequent leveling off or even a decrease. Concerted dissociations will not show up in this kind of analysis. However, an investigation of all trajectories brought up no such case of a concerted dissociation. Figure 6 shows the aBLD evolution for all systems in all trajectories. In the gas phase, overall, seven dissociations occur within the simulation time; two trajectories show double dissociations; and three, a single dissociation. In water and n-hexane, the dissociation is nearly completely suppressed. In n-hexane, no dissociation can be observed at all, and in water, only one trajectory shows a single dissociation. In the latter case, the curve of the aBLD of the dissociation trajectory has a maximum as would be typical for a double dissociation in the gas phase, but inspection of the individual bond lengths indicates that only one CH3• fragment dissociates. The other C-N bond shows no sign of dissociation at all. The origin of the maximum in the aBLD curve is that further separation of the dissociated fragment is inhibited by the solvent, and the CH3 group is again driven back toward the NN-CH3 fragment. The trajectory shows a subsequent increase in the aBLD value, which will certainly end up in a damped oscillatory motion. The trajectory in the gas phase starting with the same
J. Phys. Chem. A, Vol. 114, No. 48, 2010 12589
Figure 6. Absolute bond length difference, aBLD, of the two C-N bonds in the gas phase, water, and n-hexane.
Figure 7. Comparison of the C-N bond lengths along two trajectories in the gas phase and water showing dissociation. The initial conditions of the azomethane moiety were identical in both trajectories.
initial conditions of the azomethane moiety also shows dissociation. This fact makes the two trajectories directly comparable. Figure 7 shows the evolution of the individual C-N bond lengths in both trajectories. It can be seen that the dissociation in the gas phase and in aqueous solution starts at a similar time. Only one vibrational period of the C-N bond more is needed in the gas phase before the bond is fully broken. After the dissociation in the gas phase, the fragment can drift off freely while the cage of the water molecules keeps the CH3• group near the remainder of the molecule in the solvated trajectory. Even though this observation is based only on one occurrence, in the case of hydration, it nevertheless gives a good idea of how the dissociation is inhibited in solution. This finding is also in agreement with the very low N2 formation quantum yield of solvated azomethane in comparison with the highly efficient dissociation reaction gas phase.2,12 This cage effect has also been reported for longer time scales by Cattaneo and Persico,3 and a similar environmental effect has been observed in the simulation of the photodissociation of formamide in an argon matrix.46 Conclusions A systematic investigation of the nonadiabatic short-time photodynamics of trans-azomethane in the gas phase and polar and
12590
J. Phys. Chem. A, Vol. 114, No. 48, 2010
nonpolar solvents has been performed using a QM/MM approach describing azomethane with a highly flexible generalized valence bond, perfect pairing, complete active space (GVB-CAS) method. The goal was to obtain detailed information about the influence of solvation on the excited state deactivation of the S1 (nπ*) state of azomethane. Two solvents, water and n-hexane, strongly differing in their polarity were chosen. The first excited singlet state of azomethane is deactivated in a subpicosecond time scale involving a torsion around the central NdN double bond. This process is not affected by electrostatic solvent effects but is mainly influenced by mechanic restrictions arising from interactions with the solvent. The excitation energies show no solvatochromic shifts. The average lifetime of the excited state is elongated by ∼20-30 fs in both solvents, as compared with the gas phase, independent of the polarity of the environment. The basic feature of the evolution of structural parameters (mainly the torsional motion around the NdN bond) in the excited state are similar in the gas phase and in solution. The structural distortions needed to access the crossing seam are variable and flexible enough to allow for the partial compensation of the mechanic restrictions imposed by the solvent. A few percent of the molecules in the gas phase show dissociation of one or both C-N bonds within a short time after the transition to the ground state (i.e.,
