ARTICLE pubs.acs.org/JPCA
Semiempirical Hamiltonian for Simulation of Azobenzene Photochemistry Teresa Cusati,† Giovanni Granucci,† Emilio Martínez-Nu~nez,‡ Francesca Martini,† Maurizio Persico,*,† and Saulo Vazquez‡ † ‡
Dipartimento di Chimica e Chimica Industriale, Universita di Pisa, v. Risorgimento 35, I-56126 Pisa, Italy Departamento de Química Física y Centro Singular de Investigacion en Química Biologica y Materiales Moleculares (CIQUS), Universidad de Santiago de Compostela, E-15782 Santiago de Compostela, Spain ABSTRACT: We present a semiempirical Hamiltonian that provides an accurate description of the first singlet and triplet potential energy surfaces of azobenzene for use in direct simulations of the excited-state dynamics. The parameterization made use of spectroscopic and thermochemical data and the best ab initio results available to date. Two-dimensional potential energy surfaces based on constrained geometry optimizations are presented for the states that are most relevant for the photochemistry of azobenzene, namely, S0, S1, and S2. In order to run simulations of the photodynamics of azobenzene in hydrocarbons or hydroxylic solvents, we determined the interactions of methane and methanol with the azo group by ab initio calculations and fitted the interactions with a QM/MM interaction Hamiltonian.
1. INTRODUCTION Azobenzene and its derivatives have been the subject of many experimental19 and theoretical4,1033 studies aimed at unraveling the excited-state dynamics and the mechanism of the photoisomerization. Interest in these compounds is stimulated by a variety of proposed applications in the fields of photoresponsive materials and supramolecular devices.3446 The computational simulation approach plays an important role in the investigation of molecular dynamics, either for isolated systems or in condensed phase. The direct or “on the fly” simulation methods, after proving successful in the field of ground-state processes,4749 have found increasingly widespread application to the excited-state nonadiabatic dynamics, either ab initio4,30,5054 or semiempirical.17,19,20,25,32,33,5557 The basic advantage of the direct methods is to avoid working out analytic expressions for the potential energy surfaces (PES), which is a particularly complex task in the case of surface crossings. However, they tend to be computationally very exacting, especially for long integration times (i.e., beyond the picosecond time scale). As a consequence, the level of ab initio theory that can be used in the simulations is normally below the state of the art. On the other hand, the direct semiempirical approach is computationally much faster and can rely on quite accurate potential energy surfaces insofar as a suitable set of parameters is adopted. In fact, standard methods, such as MNDO, AM1, PMx, or OMx,5860 are targeted to reproduce properties (e.g., entalpy differences) of molecules in the ground state and near their equilibrium geometries, and an ad hoc reparameterization is often necessary to represent correctly electronic excitations and bond breaking. Moreover, the state mixing and the distorted geometries r 2011 American Chemical Society
frequently occurring in photoprocesses are best described by configuration interaction (CI) wave functions, but this treatment conflicts with the standard parameterizations that are devised to compensate for the lack of electron correlation at the SCF level. The search for an efficient electronic structure method, suitable for nonadiabatic dynamics simulations, led us to develop a semiempirical CI approach based on SCF molecular orbitals with floating occupation numbers (FOMOCI).56,61 The FOMO CI is technically simpler than state-averaged CASSCF but is equally able to describe the lowest excited states, the homolytic dissociation of bonds, and (near-)degeneracy situations for states and/or orbitals. Being faster than CASSCF, it is ideally suited in the semiempirical context but has also been implemented in ab initio codes.6264 In this work, as well as in the nonadiabatic dynamics simulations, we used a development version of the MOPAC 2002 program.65 We applied the FOMOCI method with a reparameterization worked out about 10 years ago17 to run simulations of the photodynamics of azobenzene17,19,22 and two of its derivatives.20,25 The results were in good agreement with the experimental data, except that the computed trans f cis photoisomerization quantum yield of azobenzene was slightly overestimated and we did not reproduce the weak component with long lifetimes in the fluorescence and differential absorption transients. These minor discrepancies might be due to faults in the surface hopping nonadiabatic dynamics, to inaccuracies in the Received: September 6, 2011 Revised: November 18, 2011 Published: November 22, 2011 98
dx.doi.org/10.1021/jp208574q | J. Phys. Chem. A 2012, 116, 98–110
The Journal of Physical Chemistry A
ARTICLE
Figure 1. Azobenzene structures. For TAB and CAB the respective Cartesian frames are also shown. The bond and dihedral angles referred to in the text and formulas can be identified by means of the four labeled atoms: θ t — CNNC t — C1N1N2C2; α1 t — NNC1 t — C1N1N2; α2 t — NNC2 t — C2N2N1. The NC bond torsion angles ϕ1 and ϕ2 are defined as combinations of four dihedrals (see text after eq 4): N2N1C1C3, N2N1C1C5, N1N2C2C4, and N1N2C2C6, where C3 and C5 are bonded to C1 and C4 and C6 to C2.
semiempirical PESs, and/or to the neglect of solvent effects since we ran simulations for the isolated azobenzene molecule. All three potential sources of error have been addressed in a new series of simulations of azobenzene photodynamics in the condensed phase.66 This paper presents a new parameterization of the AM1 Hamiltonian (section 2) on which the latest simulations66 are based, taking into account the most accurate ab initio results that have been published after we began working in this field.13,15,16,23,26 The new adiabatic PESs include a state-specific correction that improves their accuracy without forsaking their consistency with the semiempirical wave functions (section 3). The PESs of the lowest lying singlet states as well as their spectroscopic properties are explored and described in section 4. We also determined the interaction potentials of azobenzene with two model compounds (methane and methanol) by MP2 calculations. From these ab initio results we extracted the parameters for a QM/MM treatment of azobenzene solvent interactions in hydrocarbons and hydroxylic compounds (section 5).
electrons. Only the subset of MOs that are active in the CI has variable occupations, the others having either ni = 2 or ni = 0. The Gaussian width w can be considered an additional semiempirical parameter and is chosen so as to ensure an easy convergence of SCF and a smooth dependence of the active MO space on the molecular geometry (see below). After a few trials, we adopted the value w = 0.1 au The active space included 7 occupied and 6 virtual MOs. Of these, the highest three occupied and the lowest virtual are one π and one π* orbital and the two nitrogen lone pairs, easily identified at planar geometries and mainly localized on the NdN azo group. We included in the CI subspace all the determinants that can be generated by excitations within these four orbitals, plus all the singly excited determinants within the 13 MO active space. In total, the CI subspace counted 94 determinants. These choices were done by trial and error, also aiming at a substantial “robustness” of the CI results, vis a vis of the changing nature of the active MOs, that depend on the molecular geometry and on the interaction with the solvent. In order to stabilize the composition of the active space, we also implemented a procedure that allows the use of different semiempirical parameters in the SCF and CI calculations. We observe that the energies and wave functions of the states depend directly on the set of parameters used in the CI calculation. However, considering that a truncated configuration subspace is used, the one-electron functions also influence the results. For this reason, also the parameters used in the SCF calculation have an effect, although indirect. In our case, we fixed the parameters βS = 38 eV and βP = 35 eV for the carbon atoms to high
2. PARAMETERIZATION OF THE AM1 HAMILTONIAN In the floating occupation SCF the occupation number ni of the ith orbital depends on its energy εi according to pffiffiffi Z 2 εf ðε εi Þ2 =2ω2 ni ¼ pffiffiffiffiffiffiffi e dε ð1Þ πω ∞ where εf indicates the Fermi level energy obtained by imposing that the sum of occupation numbers be equal to the number of 99
dx.doi.org/10.1021/jp208574q |J. Phys. Chem. A 2012, 116, 98–110
The Journal of Physical Chemistry A
ARTICLE
values in the SCF calculation only, in order to increase the separation between the occupied and virtual orbitals of the phenyl groups. In this way, we have a better control of the number of aromatic π and π* orbitals belonging to the active space (at most geometries where they are easily recognized, three aromatic π orbitals and four π* are in the active space). In Appendix A we briefly describe the procedure for computing CI energy gradients with the different sets of semiempirical parameters in the SCF and CI calculations. In searching for an optimal set of semiempirical parameters we took advantage of previous reparameterizations performed by Toniolo et al. for the excited states of the benzene molecule,67 based on the AM1 and PM3 Hamiltonians. For our purposes, we found that AM1 yields better results, so we used the C atom AM1 parameter set optimized for benzene, the standard AM1 for H (as in Toniolo’s work), and we redetermined the N parameters only, with the exclusion of those defining the core repulsion. Optimization of a set of semiempirical parameters P is based on the search of a minimum of the evaluation function !2 Vi ðSÞ ðPÞ Vi ðTÞ Wi ð2Þ FðPÞ ¼ Vi ðTÞ i
data. Among the target values we considered the excitation energies of azobenzene, some energy differences between constrained or unconstrained extrema of the potential energy surfaces, and a set of geometrical parameters. To simulate the photodynamics with excitation wavelenghts λexc g 250 nm, we are interested in the first five singlet states (S0 S4) but shall also consider the first two triplets. Tables 1 and 2 report the V(T) i values along with the best V(S) semiempirical results and the i associated weights Wi. The most complete and accurate exploration of the PES of azobenzene has been done by Orlandi’s group,15,16 and we shall privilege their results in our choice of target values for the sake of consistency. The work described in this section was concluded before publication of their last ab initio study on the conical intersections of S2 and higher states,26 so it does not take into account those results. In the text, tables, and figures we shall use the following acronyms, also illustrated in Figure 1: TAB = trans-azobenzene; CAB = cis-azobenzene; ROT = rotamer = dihedral angle CNNC near to 90; INV = invertomer = one of the two angles NNC near to 180; PlanINV = planar invertomer = both phenyl groups are on the plane formed by the atoms CNNC; PerpINV = perpendicular invertomer = the phenyl group attached to the nitrogen atom that undergoes the inversion is perpendicular to the CNNC plane while the other one is on the plane. Among the energetic target values, we have first of all the vertical transition energies ΔE(S0 Sn) of TAB and CAB. They were deduced, when possible, from the maxima of the absorption spectra in the vapor phase,70,71 because in the simulations the solvent effects are taken into account independently by a QM/ MM approach. The spectral bands are well defined for TAB S1 (1Bg), TAB S2 (1Bu), CAB S1 (1B), and CAB S4 (1B). The TAB S3 (1Ag) and TAB S4 (1Bu) states are very near in energy to S2, but they have very small oscillator strengths; thus, their bands are not evident in the spectrum. The best ab initio calculations put them a fraction of eV above S2.4,11,13 The CAB S2 (1B) and S3 (1A) states probably correspond to the shoulder extending from 4.4 to 4.8 eV, next to the more intense band that we attribute to S4: two computational works confirm this interpretation.11,13 For the triplets, we adopted as target values the averages of several experimental and computational data that cover the following ranges: 1.572.09 eV for TAB T110,11,16 and 2.102.83 eV for TAB T2.10,11,72 In the case of CAB, the ab initio values for T1 range from 1.57 to 2.20 eV,10,11,16 and for T2 one finds 3.2911 and 4.12 eV. 10 Given the large discrepancy in the T2 ab initio results, we did not include ΔE(S0 T 2) in the target set; after reparameterization, the semiempirical computed value was 4.12 eV. The energy difference between CAB and TAB was measured in n-heptane (0.51 ( 0.06 eV),73 while most of the calculated values for the isolated molecule are slightly higher (0.55 0.69 eV).10,15,7477 We have chosen the value 0.55 eV, which is compatible with the PES computed by Orlandi’s group15,16 and very close to the experimental value. Concerning the conformation of TAB there has been in the past some uncertainty, but a thorough ab initio study by Briquet et al.78 firmly established that TAB is planar, with C2h symmetry. The torsions of the NC bonds are almost independent, as shown by MP2 calculations:23,79 the barrier for rotation of one phenyl group is 0.210.23 eV, and that for simultaneous rotation of both groups is 0.48 eV. We took the last value as a target, because of the easy optimization of the double saddle point with C2h symmetry and all NNCC dihedral angles at 90. Preliminary calculations
∑
Here V(S) i (P) is the value of a molecular property (such as an energy difference or an equilibrium coordinate) as computed is the corresponding target value, and with the parameters P, V(T) i Wi is the associated weight; the index i runs over all properties included in the evaluation set. The constant denominators V(T) i are not strictly necessary, but they facilitate the choice of the weights. To locate the minima of F(P) we applied the simplex method,68,69 combined with a simulated annealing procedure. The latter consists of adding to each F(P) value considered by the simplex algorithm a term ΔF, randomly chosen according to the probability distribution exp(ΔF/T). When the “temperature” T is large in comparison with the range of values spanned by F(P), the optimization is not bound to end up in the local minimum closest to the starting point. In this way, a large portion of the parameter P space is explored and deeper minima can be located. As the optimization proceeds, T is decreased so that accurate convergence to a minimum of F(P) is achieved. As already found in the previous reparameterization,17 optimization of the parameters could not eliminate a drawback, shared by AM1 with other semiempirical Hamiltonians, that yield too low N-inversion barriers and too large NNC equilibrium angles. For this reason we added a correction term to all the electronic PESs with the following expression
# " α1 P2 α2 P2 Uinv ðθ, α1 , α2 Þ ¼ P1 cos π þ 1 cos π þ 1 π P2 π P2
P3 þ cosð2θÞ P3 þ 1
ð3Þ
Here and in the following, θ is the CNNC dihedral angle, while α1 and α2 are the NNC angles (see Figure 1 and its caption for the definitions of angles and dihedrals). The parameter P1 is an energy, in eV, P2 is an angle and corresponds to the minimum of Uinv, and P3 is a dimensionless constant. The dependence of Uinv on the dihedral angle CNNC vanishes if one of the NNC angles is near to 180, as it must. The three parameters PI have been optimized together with the semiempirical ones. Before performing the optimization one has to choose the from the available experimental and ab initio target values V(T) i 100
dx.doi.org/10.1021/jp208574q |J. Phys. Chem. A 2012, 116, 98–110
The Journal of Physical Chemistry A
ARTICLE
Table 1. Energy Target Values and Semiempirical Results Obtained with the Optimized Parametersa
a
target value V(T) i
semiemp. value V(S) i
absolute error V(S) V(T) i i
relative error (%)
weight Wi 8/11
TAB, ΔE(S0 S1)
2.82
2.83
0.01
0.3
TAB, ΔE(S0 S2)
4.12
4.37
0.25
6.1
8/11
TAB, ΔE(S0 S3)
4.20
4.71
0.51
12.2
2/11 2/11
TAB, ΔE(S0 S4)
4.20
4.73
0.53
12.7
TAB, ΔE(S0 T1)
1.70
1.61
0.09
5.0
1/11
TAB, ΔE(S0 T2)
2.20
2.66
0.46
21.1
1/11
CAB, ΔE(S0 S1)
2.92
2.89
0.03
0.9
16/21
CAB, ΔE(S0 S2) CAB, ΔE(S0 S3)
4.60 4.60
5.44 4.86
0.84 0.26
18.3 5.7
4/21 4/21
CAB, ΔE(S0 S4)
5.17
5.53
0.36
7.0
16/21
CAB, ΔE(S0 T1)
1.90
1.57
0.33
17.3
2/21
ΔE(CAB-TAB),S0
0.55
0.61
0.06
10.5
2/3
ΔE 90 phenyl rotat. TAB S0
0.48
0.40
0.08
17.4
1/3
PerpINV S0 TS
1.55
1.93
0.38
24.7
4/15
PerpINV S1 (S0 geom.)
2.75
3.09
0.34
12.2
4/15
PerpINV S2 (S0 geom.) PerpINV S3 (S0 geom.)
4.55 6.80
5.30 5.48
0.75 1.32
16.5 19.4
4/15 2/15
PlanINV S0
2.10
2.16
0.06
3.0
1/15
ROT S0 TS
1.70
2.17
0.41
23.9
4/13
ROT S1 (S0 TS)
2.53
2.50
0.03
1.0
4/13
ROT S2 (S0 TS)
2.91
3.24
0.33
11.4
2/13
ROT S3 (S0 TS)
4.48
5.90
1.42
31.8
1/13
ROT S4 (S0 TS)
4.48
6.65
2.17
48.5
2/13
TAB S1 min TAB S0 (S1 min)
2.27 1.02
2.36 0.65
0.09 0.37
4.0 35.8
4/5 2/5 4/5
PlanINV S1
3.30
2.92
0.38
11.4
ROT S0/S1 CoIn
2.27
2.24
0.03
1.5
8/9
TAB S0/S1 CoIn
3.30
2.96
0.34
10.4
8/9
ROT S2 min
2.30
3.18
0.88
38.0
2/9
ROT T1 min
1.25
1.11
0.14
11.2
1/2
ROT S0 (T1 min)
1.50
1.40
0.10
6.8
1/2
Energies, in eV, given as differences from TAB S0 (except for the vertical excitation energies of CAB).
showed that a good reproduction of the NC torsional potential using semiempirical methods is quite difficult. For this reason, we gave a small weight to this target value and introduced an ad hoc correction after optimization of the parameters (see the end of this section). Configuration interaction calculations, with geometries optimized at thr CASSCF level,10,74 show that the preferred pathway for isomerization in the ground state is the N inversion with the phenyl ring attached to the inverting nitrogen perpendicular to the rest of the molecule. The experimental values of activation energy for conversion CAB f TAB should therefore be referred to this geometry (perpendicular invertomer, PerpINV). In n-heptane, the activation ΔH of the cis f trans conversion is about 1 eV,80 i.e., 1.55 eV with respect to the TAB. Two more saddle points connecting the TAB and CAB minima have been located by geometry optimizations: one related to the torsion of the NdN double bond (rotamer, ROT),10,16,74 and one for the N inversion with all atoms on the same plane (planar invertomer, PlanINV).10,74 Target values were assigned to these geometries according to the quoted computational results: from our previous work10 we took the ground-state energy of the ROT TS and the transition energies of the PerpINV TS, while those of ROT TS were found in Gagliardi et al.15 At a semiempirical level,
for the ROT TS we performed a saddle point search, while for PerpINV and PlanINV we located two constrained minima, fixing one of the two NNC angles at 180. Our knowledge of the excited-state PESs outside the FranckCondon region is essentially based on ab initio calculations of CASSCF type, in some cases with addition of perturbative CI treatments.10,12,1416,26 Optimization of the internal coordinates in the S1 PES, with the constraint of a trans planar geometry, decreases the energy by 0.55 eV.14,16 Combining these data with the experimental vertical excitation energy the target value of 2.27 eV is obtained. From Cembran et al.16 we get the energy difference S0 S1 at the same geometry (1.25 eV), which puts the energy of S0 equal to 1.02 eV. For the invertomer optimized in the S1 state (PlanINV with one of the NNC angles fixed at 180) we have chosen a value of 0.48 eV above the vertical excitation energy as a compromise between two rather different ab initio values: 0.20 and 0.86 eV, respectively, from Diau14 and Cattaneo and Persico.10 The S0 and S1 surfaces present a crossing seam.17,22 The lowest energy point (ROT S0/S1 CoIn) coincides with the minimum of the S1 PES, which is found along the NdN torsional coordinate. Its energy has been determined by Cembran et al.16 in 2.03 eV with respect to TAB. Considering that in their calculations the 101
dx.doi.org/10.1021/jp208574q |J. Phys. Chem. A 2012, 116, 98–110
The Journal of Physical Chemistry A
ARTICLE
Table 2. Target Values of Internal Coordinates and Semiempirical Results Obtained with the Optimized Parametersa target value V(T) i
absolute error V(S) V(T) i i
relative error (%)
weight Wi
TAB, RNN
1.243
1.280
0.0369
3.0
TAB, RNC
1.422
1.434
0.0122
0.9
1/5
2.6850
2.3
3/5 10/61
TAB, — NNC
115.1
117.8
1/5
CAB, RNN
1.242
1.254
0.0117
0.9
CAB, RNC
1.435
1.452
0.0165
1.2
10/61
3.1206
2.6
30/61
CAB, — NNC CAB, — NNCC CAB, — CNNC TS ROT S0, RNN TS ROT S0, RNC TS ROT S0, — NNC TS ROT S0, — CNNC
122.4
125.5
62.0
63.2
4.2 1.304 1.370
1.2104
2.0
10/61
0.3 1.294
3.9414 0.0096
93.8 0.7
1/61 1/8
1.395
1/8
0.0250
1.8
122.2
128.4
6.1972
5.1
3/8
85.3
91.4
6.1189
7.2
3/8
PerpINV S0, RNN
1.233
1.239
0.0061
0.5
1/6
PerpINV S0, RNC1
1.425
1.438
0.0127
0.9
1/6
PerpINV S0, RNC2
1.318
1.391
0.0732
5.6
1/6
116.0 1.247
125.9 1.245
9.9521 0.0016
8.6 0.1
3/6 1/6
PlanINV S0, RNC1
1.424
1.441
0.0168
1.2
1/6
PlanINV S0, RNC2
1.354
1.394
0.0395
2.9
1/6
8.2112
7.0
3/6
PerpINV S0, — NNC1 PlanINV S0, RNN
PlanINV S0, — NNC1
118.0
126.2
TAB S1 min, RNN
1.253
1.260
0.0068
0.5
1/5
TAB S1 min, RNC
1.366
1.407
0.0412
3.0
1/5
3.5320
2.7
3/5
0.0742 0.0007
5.7 0.0
1/6 1/6
TAB S1 min, — NNC PlanINV S1, RNN PlanINV S1, RNC1 PlanINV S1, RNC2 PlanINV S1, — NNC1
128.8 1.309 1.403 1.360 118.7
132.3 1.235 1.404 1.382 141.1
0.0224
1.6
1/6
22.3831
18.9
3/6
ROT S0/S1 CoIn, RNN
1.261
1.267
0.0058
0.5
1/14
ROT S0/S1 CoIn, RNC
1.397
1.413
0.0161
1.2
1/14
1.371
1.398
0.0265
1.9
1/14
ROT S0/S1 CoIn, — NNC
117.3
128.0
10.7301
9.1
3/14
ROT S0/S1 CoIn, — CNNC
136.0 94.0
139.5 96.1
3.5095 2.1173
2.6 2.3
3/14 3/14
ROT S0/S1 CoIn, — NNCC
25.9
9.4
16.5335
63.8
1/14
18.8
9.2
9.5815
51.0
1/14
1.230
0.0352
2.9
1/5
1.378
0.0517
3.9
1/5
0.3267
0.2
3/5
TAB S0/S1 CoIn, RNN TAB S0/S1 CoIn, RNC TAB S0/S1 CoIn, — NNC a
semiemp. value V(S) i
1.195 1.326 156.6
156.9
Distances in Angstroms and angles in degrees.
who provided geometry data for TAB S0 and S1, CAB S0, ROT S0 TS, and ROT S0/S1 CoIn. The calculation level, CASSCF with a large active space, guarantees a good agreement with the experimental data for the ground state of TAB (see Tsuji et al.79 and references therein) and CAB.81 The geometries PlanINV and PerpINV S 0 were taken from our previous CASSCF calculations.10 We corrected our distances and bond angles, considering the differences found with respect to those of Cembran et al.16 for the TAB and CAB S0 geometries. For the PlanINV and PerpINV geometries, the two RNC bond distances are different: the longer one (RNC1) is relative to the sp2 nitrogen, while the shorter one (RNC2) concerns the nitrogen undergoing inversion. Also, for PlanINV S1 we can only rely on our previous work,10 but in this case we did not apply any correction because of the lack of more accurate data for similar geometries in the
vertical transition energy for TAB and CAB is understimated by 0.20.3 eV, we took a target value of 2.27 eV for the conical intersection. The crossing seam can also be reached by keeping a planar geometry of TAB type, deformed for the symmetric opening of the two angles NNC, as shown by Diau:14 in this case, the energy is slightly higher than the vertical excitation one (near 3.3 eV). The energy value of the minimum of S2 has been taken from Gagliardi et al.15 and that of T1 as well as the energy of S0 at the same geometry from Cembran et al.16 Concerning the target values of internal coordinates relative to critical points of the PESs we note that the internal consistence of such data is important. In fact, the dynamical behavior depends on the optimal geometry differences between ground and excited states or between isomers and transition states. For this reason we relied, whenever possible, on the data of Cembran et al.,16 102
dx.doi.org/10.1021/jp208574q |J. Phys. Chem. A 2012, 116, 98–110
The Journal of Physical Chemistry A
ARTICLE
Table 3. Semiempirical Parameters for Azobenzenea units
Table 4. Parameters That Define the Added Potentials Uinv and Uph
N
C
H 11.396427
parameter P1
eV
0.382809
6.173787
P2 P3
degrees
82.172261 13.416564
1.188078
P4
eV
Uss
eV
68.388062*
49.536242*
Upp
eV
56.045018*
33.722921*
βs
eV
11.018709*
13.797578*
βp
eV
19.591942*
10.113264*
ζs
bohr1
2.169870*
1.412377*
ζp
1
bohr
1.917021*
1.749266*
gss
eV
16.667653*
12.492690*
gsp gpp
eV eV
11.470530* 13.339308*
11.571701* 12.017613*
gp2
eV
10.242204*
7.923532*
hsp
eV
1.794342*
2.739211*
α
Å1
P5
value
0.0949700 0.660000
12.848000
The optimized semiempirical parameters are listed in Table 3 and those of eqs 3 and 4 in Table 4. The semiempirical results of Tables 1 and 2 were obtained with the corrections so far considered, Uinv and Uph.
2.947286
2.626199*
2.882324
K1
0.025251
0.011355
0.122796
K2
0.028953
0.045924
0.005090
K3
0.005806
0.020061
0.018336 5.000000
K4 L1
Å1
5.000000
0.001260 5.000000
L2
Å1
5.000000
5.000000
5.000000
L3
Å1
2.000000
5.000000
2.000000
L4
Å1
M1
Å
1.500000
1.600000
1.200000
M2
Å
2.100000
1.850000
1.800000
M3
Å
2.400000
2.050000
2.000000
M4
Å
3. STATE-SPECIFIC CORRECTIONS OF THE PESS The reparameterization described in the previous section was quite successful for the S1 (n f π*) excited state: its semiempirical PES reproduces very well the best ab initio data15,16,26 and in the FranckCondon region is even more accurate than them, thanks to inclusion in the target values of spectroscopic data (see the next section for details). For the upper π f π* states, S2S4, the results were slightly less satisfactory. For instance, all vertical excitation energies are overevaluated, the largest relative error being 18% for CAB S2 (notice however that the CAB S2 transition is less well characterized than the TAB one, see discussion of the target values in the previous section). On the other hand, new accurate ab initio data26 became available after our reparameterization was completed, including determination of the S1/S2 conical intersection, which is very important for decay of the S2 state.17 For these reasons, we devised a statespecific correction to the PESs, limited almost exclusively to the S2S4 states. This correction was not applied in our simulations of the photodynamics following n f π* excitation66 but only in the case of πfπ* excitation. State-specific corrections applied to the adiabatic PES, i.e., to the eigenvalues Uk of the electrostatic Hamiltonian, must not alter the ordering of the states within a given spin manifold (same eigenvalue of S2) in order to preserve the consistency between electronic energies and wave functions. For instance, at a nearcrossing situation one has large nonadiabatic couplings because of the electronic wave function mixing, and it is of paramount importance that the quasi-degeneracy of energies and the wave function switching still coincide in the nuclear coordinate space, after correcting the PESs. In direct dynamics simulations one cannot easily alter the wave functions, from which the nonadiabatic couplings are computed, so the locations of the energy crossings must also remain unchanged. In general, two kinds of corrections can be envisaged: one is to add the same potential energy term to all state energies, as already discussed in the previous section; the second consists in changing the energy difference between any two consecutive states, Uk Uk1, by a factor Fk > 0, function of the nuclear coordinates Q. In this way, the set of molecular geometries where the two states k and k 1 are degenerate remains unaltered. However, the energy of the crossing seam and the position of the minimum energy conical intersection as well as those of the local minima of the PESs can be modified. In the case of a two-state avoided crossing, changing the eigenvalue difference by a factor F without changing the wave functions amounts to multiplying by F all the elements of the
5.000000
2.650000
a
The starred values differ from the standard AM1. The names of the parameters are those used in the MOPAC 2002 documentation65 and in the literature.58,60
same state to compare with. The geometry of the S0/S1 conical intersection of TAB, obtained by symmetric opening of the NNC angles, was calculated by Diau.14 The values of RNN and RNC have been corrected by comparison between the data of Diau14 and that of Cembran et al.16 for the TAB S1 minimum. As remarked above, we could not obtain an accurate torsional potential of the phenyl rings about the NC bonds. Therefore, we applied a further additive correction to all the PES, besides that expressed by eq 3, in the form Uph ðθ, α1 , α2 , ϕ1 , ϕ2 Þ ¼ P4
unit
P5 cos θ ð1 þ cos α1 Þ P5 þ 1
ð1 þ cos α2 Þðsin2 ϕ1 þ sin2 ϕ2 Þ
ð4Þ
Here, ϕ1 and ϕ2 are the torsion angles of the two phenyl groups: ϕ1 = ( — N2N1C1C3 + — N2N1C1C5 π)/2 and ϕ2 = ( — N1N2C2C4 + — N1N2C2C5 π)/2. The parameters P4 and P5 have been determined after optimization of the semiempirical parameters. P4 is an energetic factor, and P5 modulates the dependence on the angle θ t CNNC. Before addition of the Uph potential, the phenyl groups of the CAB isomer lied closer than expected to the CNNC plane: we found equilibrium angles ϕ1 = ϕ2 = 46, while our target value was 62. Instead, for TAB the torsional potential of the phenyl groups was too flat, in comparison with ab initio calculations and experimental evidence.23,78,79 To determine P4 and P5 we considered the torsional barriers for TAB, the equilibrium geometry of CAB, and the CAB TAB energy difference. 103
dx.doi.org/10.1021/jp208574q |J. Phys. Chem. A 2012, 116, 98–110
The Journal of Physical Chemistry A
ARTICLE
Hamiltonian matrix in the diabatic representation.82,83 We feel that this conservative way of correcting the PES could be profitably applied also to direct dynamics simulations with ab initio electronic structure calculations, which cannot be as accurate as those affordable when characterizing a few critical points of the PESs. Thus, for each state within the same spin manifold, starting with the second one (k = 1), the corrected PES U0 k is recurrently computed as U 0k ðQ Þ ¼ U 0k1 ðQ Þ þ Fk ðQ Þ½Uk ðQ Þ Uk1 ðQ Þ
Table 5. State-Specific Correction Parameters (αk,left and αk,right angles in degrees) state pair
k
αk,left
αk,right
Fk,left
Fk,right
S0 S1
1
110
115
0.911
1.000
S1 S2 S2 S3
2 3
110 0
150 180
0.750 0.517
1.000 0.517
Table 6. Target Values and Semiempirical Results Obtained with and without State-Specific Corrections (SSC) of the PESa
ð5Þ
The first PES (ground state for each manifold) can be corrected by an additive term, which is automatically applied to all the other states U00 ðQ Þ ¼ U0 ðQ Þ þ Uadd ðQ Þ
semiemp.
ð6Þ
In our case, Uadd is the sum of the correction terms Uinv and Uph (see previous section). Of course, also the derivatives of the potential energies with respect to the nuclear coordinates are altered ∂Uk0
∂Q r
¼ ð1 Fk Þ
0 ∂Uk1
∂Q r
þ
∂Fk ∂Uk ðUk Uk1 Þ þ Fk ∂Q r ∂Q r
ð7Þ
Thus, when such corrections are introduced in simulations, one needs to calculate the energy gradients for all the PESs lying below the one where the trajectory is running. To correct the S2S4 PES in the FranckCondon regions of CAB and TAB and the energy of the S1/S2 conical intersection at transoid geometries, we introduced the F1, F2, and F3 factors for the S0 S1, S1 S2, and S2 S3 energy differences, respectively. The S4 state was not an explicit object of state-specific corrections, because it is relatively unimportant in the photodynamics and because accurate target data were not available. Nevertheless, the S4 PES is affected by the corrections applied to the lower states. We defined the Fk factors as functions of the NNC angles, because we aimed at a better approximation of the π f π* state energies for small NNC, where one finds the minima of S0 and S2, and the S1/S2 CoIn. For the latter, we added new target values taken from the work of Conti et al.26 Since they overestimate the vertical transition energy of the S2 state by ∼0.1 eV, we lowered their prediction for the CoIn energy (3.84 eV) to 3.75 eV. The form of the correction factors is Fk ðxÞ ¼ Fk, lef t þ ðFk, right Fk, lef t ÞSðxÞ 2 cos αk, lef t cos α1 cos α2 2ðcos αk, lef t cos αk, right Þ
and S(x) is the sigmoid function 8 > for x e 0 :1 for x g 1
without SSC
semiemp. with SSC
TAB, ΔE(S0 S2)
4.12
4.37
4.04
TAB, ΔE(S0 S3)
4.20
4.71
4.22
TAB, ΔE(S0 S4)
4.20
4.73
4.24
CAB, ΔE(S0 S2)
4.60
5.44
4.88
CAB, ΔE(S0 S3)
4.60
4.86
4.59
CAB, ΔE(S0 S4) PerpINV S2 (S0 geom.)
5.17 4.55
5.53 5.30
4.98 5.27
PerpINV S3 (S0 geom.)
6.80
5.48
5.37
ROT S2 (S0 TS)
2.91
3.24
3.16
ROT S3 (S0 TS)
4.48
5.90
4.53
ROT S4 (S0 TS)
4.48
6.65
5.28
ROT S2 min
2.30
3.18
3.11
TAB S1/S2 CoIn
3.75
4.05
3.75
TAB S1/S2 CoIn, RNN TAB S1/S2 CoIn, RNC
1.350 1.310
1.341 1.404
1.341 1.405
TAB S1/S2 CoIn, — NNC
109.0
110.4
110.2
a
Energies, in eV, given as differences from TAB S0. Distances are in Angstroms and angles in degrees.
however, we find the S1/S2 CoIn, so this part of the PESs is very important in the case of π f π* excitation. The S2 PES was modified for a broader range of angles (α < α2,right = 150), while for the S3 S2 energy difference we used a constant factor, F3,left = F3,right. In Table 6 we report the only semiempirical values that are affected by the state-specific corrections. The vertical excitation energy of TAB to S2 is lowered and so is the S1/S2 conical intersection, but the energy difference between these two important points remains about the same. Larger changes occur in the CAB FranckCondon region, which is characterized by steep slopes in the excited states. In the next section we present a fuller characterization of the excited-state PESs.
ð8Þ
where x¼
target value
ð9Þ
4. POTENTIAL ENERGY SURFACES AND TRANSITION DIPOLES Thanks to the reduced cost of the semiempirical calculations we could perform an extensive exploration of the PESs. Figures 2, 3, and 4 show the PESs of S0, S1, and S2, respectively, as functions of the CNNC dihedral and of one of the NNC angles (NNC1). For each state and for each pair of — CNNC, — NNC1 values, all other internal coordinates have been optimized, including the NNC2 angle. When more than one minimum is found, we choose the lowest energy one. As a consequence, in some cases two neighboring points in the PES correspond to rather different geometries and in particular to different rotations of the phenyl groups around the NC bonds. The minimum energy pathways
ð10Þ
As in the previous section, α1 and α2 are the NNC angles. The final parameters αk,left, αk,right, Fk,left, and Fk,right are given in Table 5. The S1 PES has been modified only for NNC angles smaller than α1,right = 115, i.e., in a region hardly accessible by trajectories starting with the n f π* excitation; at about α = 110, 104
dx.doi.org/10.1021/jp208574q |J. Phys. Chem. A 2012, 116, 98–110
The Journal of Physical Chemistry A
ARTICLE
Figure 3. Potential energy map of S1 as a function of one of the NNC angles and of the CNNC dihedral angle. The radial coordinate is 180 — NNC, and the angular one is — CNNC. All other internal coordinates have been optimized. Contour line increments: red lines, 0.1 eV; blue lines, 0.5 eV. In the yellow shaded regions the energy gap with the ground state is less than 0.3 eV, and in the orange ones it is less than 0.15 eV. Two schematic trajectories going from the TAB and CAB FranckCondon points to the region of the conical intersection with S0 are drawn by dashed black lines.
Figure 2. Potential energy map of S0 as a function of one of the NNC angles and of the CNNC dihedral angle. The radial coordinate is 180 — NNC, and the angular one is — CNNC. All other internal coordinates have been optimized. Contour line increments: red lines, 0.1 eV; blue lines, 0.5 eV.
for double-bond torsion and N inversion in the same three states are shown in Figure 5. The figures have the form of polar plots, where the radial coordinate is 180 — NNC1 and the angular one is — CNNC. In all diagrams a straight line connecting the cis and trans regions represents the pure inversion pathway and a semicircle the torsional one. In the ground-state PES the minimum corresponding to the TAB isomer is unique. On the contrary, there are two enantiomeric forms of CAB that differ mainly because of the rotation of the phenyl rings and also because of the sign of the CNNC dihedral. However, the semiempirically computed equilibrium value of — CNNC is so small ((0.3) that the two points almost coincide in Figure 2. In the same diagram we also marked the PerpINV and ROT transition states: it can be seen that they are only separated by a modest ridge topping at about 0.1 eV above the ROT TS. To the best of our knowledge, an accurate ab initio characterization of the ground-state transition states and of their relationships is still missing. In the S1 PES the TAB- and CAB-optimized geometries with — CNNC, respectively, equal to 180 and 0 are transition states, with the negative curvature mode largely coinciding with the CNNC torsional coordinate. In fact, while there is a PlanINV TS along the pure inversion pathway, the torsion leads without barriers to a set of four equivalent ROT minima, starting from either TAB or CAB. The minima at — CNNC = 96 are enantiomers of those at — CNNC = 264 (i.e., — CNNC = 96). Since every minimum has forsaken the C2 symmetry, the two minima on the same side correspond to — NNC1 > — NNC2 and to — NNC1 < — NNC2, respectively. The ROT minima coincide with the optimized S0 S1 conical intersection points, and in their proximity the crossing seam lies on the PES shown in Figure 3; in other parts of the — CNNC/ — NNC diagram the crossing seam is at higher energies with respect to the optimized S1 PES. A simpler one-dimensional view of the relationship between the S0 S1 crossing seam and the torsional minimum energy pathway in S1 is offered in Figure 5. The
Figure 4. Potential energy map of S2 as a function of one of the NNC angles and of the CNNC dihedral angle. The radial coordinate is 180 — NNC, and the angular one is — CNNC. All other internal coordinates have been optimized. Contour line increments: red lines, 0.1 eV; blue lines, 0.5 eV. In the yellow shaded regions the energy gap with S1 is less than 0.3 eV, and in the orange ones it is less than 0.15 eV.
yellow/orange shading in Figure 3 marks the regions with a small S0 S1 energy gap, resulting in a non-negligible nonadiabatic transition rate. Quite clearly the shape of the PES attracts both the TAB and the CAB isomers toward the fast decay regions. 105
dx.doi.org/10.1021/jp208574q |J. Phys. Chem. A 2012, 116, 98–110
The Journal of Physical Chemistry A
ARTICLE
Figure 5. Potential energy curves of the first three singlets, as functions of the CNNC dihedral (left) and of one of the NNC angles. All other internal coordinates have been optimized for each state. The blue dotted line in the left panel indicates the S0/S1 crossing seam. The full circles show the vertical excitation energies (FranckCondon points), and the empty one shows the S1/S2-optimized conical intersection.
Table 7. Computed and Experimental Spectroscopic Dataa TAB
CAB
S1
S2
S3
S4
S1
S2
S3
1
1
1
1
1
1
1
property method 1 Bg
1 Bu
are facilitated by the excess vibrational energy generated by the electronic excitation. The whole reparameterization process was driven by energy and molecular geometry data. The only wave function property we considered was the symmetry of the electronic states, which is important in assigning the energies of S2, S3, and S4. Among the wave function-dependent properties, the transition dipoles are most directly related with photochemistry. In Table 7 we compare the semiempirically computed spectroscopic quantities with experimental70 and ab initio data. The latter were obtained by CC213 and CASPT226 calculations. All semiempirical and theoretical results in Table 7 refer to vertical excitation starting from the equibrium ground-state geometry, and the experimental values were obtained from published absorption spectra.70 At D2h geometries the S0 f S1 (n f π*) transition dipole of TAB vanishes; thus, to evaluate the intensity of this band one must take into account the distorsions of the equilibrium geometry that arise because of the zero-point vibrations and because of thermal excitation of the lowest frequency modes. An accurate study based on a TDDFT treatment of the n f π* transition showed that the thermal effects are preponderant in this case and yielded an oscillator strength f = 0.0110.013 depending on the computed PES.23 With the present semiempirical PES and wave functions a simpler Monte Carlo integration in the internal coordinate space (680 000 steps) also yielded f = 0.013. Knowing the direction of the n f π* transition dipole is important to interpret the fluorescence anisotropy measurements, performed to shed light on the photoisomerization and decay mechanism of TAB.3,7 Figure 1 shows the coordinate frame we use: for both TAB and CAB the N atoms lie on the x axis and the z axis coincides with the C2 symmetry axis. By averaging the TDDFT transition dipole one gets a vector that lies in the xy plane and makes an angle with the x axis of 53 (the signs of μx and μy being opposite).23 The semiempirical Monte Carlo result is in good agreement, yielding an angle of 48. In both cases the S0 f S1 transition turns out to be approximately parallel to the long axis of the molecule, just as the S0 S2 one, from which it borrows its intensity (see the transition dipole components in Table 7). The CAB n f π* transition dipole also lies in the xy plane and yields an f value in good agreement with the experimental one. In the enantiomer depicted in Figure 1, to which the results of Table 7 refer, both phenyl rings are rotated out of the xz plane in the counterclockwise sense around their respective NC axes. For this enantiomer the μx and μy components have opposite
2 Bu 2 Ag
1B
2B
S4
2 A 31B
ΔE, eV
SE
2.83
4.04 4.24 4.22
2.89
4.88 4.59 4.98
ΔE, eV
CC
2.84
4.04 4.44 4.45
3.00
4.49 4.65 4.82
ΔE, eV
PT
2.53
4.23 4.46 4.53
2.72
4.49 4.54 4.72
ΔE, eV
exp
2.82
4.12
2.92
f
SE
0.000
1.08 0.15 0.00
0.019
0.28 0.02 0.27
f f
CC PT
0.000 0.000
0.85 0.02 0.00 0.62 0.07 0.00
0.030 0.015
0.06 0.02 0.07 0.04 0.01 0.19
f
exp
0.007
0.53
0.017
μx, au
SE
0.00
2.34 1.21 0.00
0.31
1.39 0.00 1.35
μy, au
SE
0.00
2.36 0.04 0.00 0.42
0.49 0.00 0.59
μz, au
SE
0.00
0.00 0.00 0.00
0.00
5.17
0.21
0.00 0.43 0.00
a
The SE, CC, PT, and exp labels mean semiempirical (this work), ab initio CC2,13 ab initio CASPT2,26 and experimental, respectively. The theoretical energy differences ΔE and transition dipoles μB are computed for the vertical excitation from the S0 minimum, and the oscillator strengths are given by f = (2ΔEμ2)/3 (in au). The experimental ΔE values are obtained from the maxima of the absorption bands and f from graphical integration of the published spectra.70
Figure 3 is probably the most accurate representation of the S1 PES obtained to date. We cannot claim the same accuracy for S2 and the higher states, not only because the semiempirical method is less apt to describe these states but also because the relevant spectroscopic data are not fully unambiguous and the ab initio results are less accurate and complete than in the case of S1. In the S2 PES we find again a real minimum for the TAB geometry and two shallow minima at — CNNC ( 15 corresponding to the enantiomeric forms of CAB. Both the TAB and the CAB minima are close to the FranckCondon points, i.e. to the equilibrium geometry of the ground state. They can be escaped through easily accessible transition states to reach the two equivalent global minima that lie close to half way between the two isomers along the torsional path. Another way to leave the TAB and CAB FranckCondon regions of S2 is by nonadiabatic transitions to S1. In a wide region around the CAB minima the energy difference S1 S2 is small (area denoted with yellow/orange shading), and the same occurs close to the TAB minimum. Of course, both processes, geometrical relaxation and decay to S1, 106
dx.doi.org/10.1021/jp208574q |J. Phys. Chem. A 2012, 116, 98–110
The Journal of Physical Chemistry A
ARTICLE
to the experimental value of 0.21 attributed to the band centered at 5.17 eV one should add about 0.05 for its long wavelength tail for a total f = 0.26. We obtain two close lying strong B transitions plus a weak A transition (notice that with respect to the ab initio results we get a reversed energy ordering of the 21B and 21A states). Summing up their oscillator strengths we get f = 0.57, which is overestimated as in the case of TAB. The ab initio CASPT2 values are instead quite good, while the CC2 ones are underestimated. The transition dipoles of the stronger transitions are approximately aligned with the NN axis and more closely if we consider their average.
5. QM/MM PARAMETERS FOR THE INTERACTION WITH SOLVENT MOLECULES We computed the interaction potential of azobenzene with two molecules at a sufficiently accurate ab initio level in order to extract the QM/MM parameters to be used in simulations of the photochemistry in solution. We chose methane and methanol as representatives of two classes of solvents, alkanes and alcohols. We assumed that the standard force-field parameters correctly describe the interaction of alkyl and hydroxyl groups with the phenyl rings; therefore, we concentrated on the interactions with the azo group. All calculations were performed with Dunning’s cc-pVDZ basis set and the GAUSSIAN03 package.84 To take into account the dispersion forces, certainly quite important, we used the MP2 method. Basis set superposition errors are not negligible; therefore, we applied the counterpoise (CP) correction of the BSSE.85 Full optimizations with the CP correction being unpractical, we first optimized the azobenzenesolvent complex at the MP2/ccpVDZ level. Then we computed the potential energy curves either with or without CP corrections, starting with the optimized complex geometry and only varying the distance RNH between the closest N and H atoms (N belonging to azobenzene and H to the solvent molecule). The rationale for this procedure is that the BSSE depends essentially on the distance between the two interacting partners, here taken as RNH, and much less on all other internal coordinates of the complex. For the TABmethane complex the MP2/cc-pVDZ optimization placed CH4 above the plane of TAB, with the line connecting the closest NH pair approximately perpendicular to the plane of the C1N1N2 group (TAB-Met-Perp geometry, see Figure 6). We also optimized a complex with the constraint that the closest H atom of methane lies in the C1N1N2 plane (TABMet-Plan geometry). The complex with CAB has CH4 approximately on the C1N1N2 plane. Methanol prefers to point its hydroxyl hydrogen toward the lone pairs of both TAB and CAB, and its OH bond lies approximately in the C1N1N2 plane (TAB-MeOH-Plan and CAB-MeOH geometries, see Figure 6). Also, for the TABmethanol complex we performed a constrained optimization, forcing the hydroxyl H atom to be on a plane that contains the N atoms, perpendicular to that of the C1N1N2 group (see again Figure 6). The well depths (see Table 8) were obtained by subtracting the potential minima from the energy of the separated molecules, optimized at the MP2/cc-pVDZ level. We see that the CP correction of the BSSE reduces all the De values by a factor larger than 2 and elongates considerably the equilibrium distances RNH. In the molecular dynamics simulations we use a QM/MM representation of the solutesolvent interaction.86 The MM potential concerning only solvent molecules is computed by the
Figure 6. Azobenzenemethane and azobenzenemethanol complexes as obtained by MP2/cc-pVDZ optimizations.
Table 8. Well Depths De and Equilibrium Distances RNH for the AzobenzeneMethane and AzobenzeneMethanol Complexesa MP2/no CP
MP2/CP
semiemp.
De
RNH
De
RNH
De
RNH
TAB-Met-Plan
1.48
2.83
0.57
3.18
0.70
3.11
TAB-Met-Perp
2.37
2.77
0.92
3.08
0.96
3.08
2.31 11.25
2.80 1.93
0.98 4.24
3.03 2.11
1.12 3.54
3.28 2.29
complex
CAB-Met TAB-MeOH-Plan TAB-MeOH-Perp CAB-MeOH
4.69
2.10
1.95
2.38
2.26
2.49
10.53
2.00
4.98
2.18
5.09
2.17
a
Ab initio (MP2 with and without CP correction) and semiempirical QM/MM results. Energies in kcal/mol and distances in Angstroms.
signs, while for the other enantiomer they agree in sign. As already shown by Ootani et al.30 and Weingart et al.33 the cis f trans isomerization can occur by twisting the NdN double bond in one of two opposite directions, i.e., by a clockwise or a counterclockwise rotation of the phenyl ring labeled 1 around the x axis, the other ring moving in the opposite sense. The counterclockwise rotation promptly moves the two rings away from each other and, assuming the z axis remains fixed, produces TAB with the orientation displayed in Figure 1. This photoisomerization pathway is the preferred one30,33 and entails a minimal change in the orientation of the transition dipole with respect to the molecular frame. On the contrary, the clockwise rotation requires a preemptive enantiomerization to avoid repulsion between the two rings, thus resulting in a substantial reorientation of the dipole. Such considerations are relevant to shed light on the mechanism of photoinduced anisotropy.35,4346 In the case of the UV absorption of TAB the oscillator strength f reported in the table refers to the whole band centered at 4.12 eV that includes at least the transitions to S2, S3, and S4, the first one being by far the strongest. As already observed, the S0 f S2 transition dipole is approximately aligned with the long axis of the molecule. All computational methods overestimate f and the semiempirical more than the ab initio ones. In the case of CAB, 107
dx.doi.org/10.1021/jp208574q |J. Phys. Chem. A 2012, 116, 98–110
The Journal of Physical Chemistry A
ARTICLE
6. CONCLUDING REMARKS To obtain a quantum chemical description of the lowest lying electronic states of azobenzene, suitable for “on the fly” excitedstate dynamics, we resorted to the already developed FOMOCI semiempirical method. Some technical advances, described in section 3 and in Appendix A, proved to be useful. The new set of semiempirical parameters yielded a very accurate description of the S1 surface, already exploited in a full set of simulations66 (trans f cis and cis f trans photochemistry in vacuo and in solution). We also obtained a fairly accurate description of the higher lying π f π* states, as far as the fragmentary available ab initio and experimental data allow us to judge. We note that probably only a little more than the FranckCondon region of the S2 S4 PESs is relevant to photochemistry, because these states decay very fast to S1 before undergoing a considerable geometrical relaxation.17 We also determined the interactions of azobenzene with two model molecules, methane and methanol, by MP2 calculations. From these ab initio data we extracted two sets of QM/MM parameters that can be used in the frame of the FOMOCI method to run condensed-phase simulations of the nonadiabatic dynamics where azobenzene interacts with alkanes or alcohols.66 Overall, the work described in this paper is preliminary to investigate photochemical and photophysical phenomena in a variety of conditions not yet tackled by computational simulations, such as the dependence of condensed-state photochemistry on the exciting wavelength, the effects of high viscosity and solid matrices and the photoinduced anisotropy, the photosensitized isomerization,89 or the ability of different spectroscopic techniques to probe the excitedstate dynamics of azobenzene and its derivatives.
Table 9. LennardJones Parameters and Charge Factor for the Semiempirical QM/MM Representation of Azobenzene Methane and AzobenzeneMethanol Interactions parameter
methane
methanol
εN, kcal/mol
0.000018
0.1490
σN, Å
5.782
2.081
εH, kcal/mol
0.0289
σH, Å
2.837
Q fac
0.4097
package TINKER87 that we interfaced to MOPAC2002. In this work we adopted the all-atom OPLS force field.88 The interaction between nuclei and electrons of the solute (QM subsystem) and the atomic charges of the solvent is included in the semiempirical Hamiltonian, so that the electrostatic field generated by the solvent affects the solute electronic wave functions and energies in a state-specific way.86 The QM atoms also interact with the MM ones through 6-12 LennardJones (LJ) potentials. The QM/MM interaction therefore depends on three kinds of parameters that are held fixed during the dynamics: the MM atomic charges and the ε and σ LJ constants for both the QM and the MM atoms. We made use of the OPLS atomic charges, but in the case of methanol we scaled them by a factor Q fac, which plays the role of an inverse dielectric constant. The LJ parameters for the aromatic (QM) and alkylic (MM) C and H atoms also kept their OPLS values, and the same holds for the hydroxyl oxygen atom. Therefore, we only optimized the N atom LJ parameters and in the case of methanol those of the hydroxyl H atom (which is not a center of LJ interactions in OPLS) and the Q fac factor. We stress that these parameters are only used in the QM/MM part of the calculation and do not affect the MM potential. The optimizations were performed with the same procedure already used for the semiempirical parameters (section 2), having as targets the three pairs of De and RNH values obtained with CP corrections for methane and methanol. The optimized parameters are shown in Table 9, while the De and RNH values obtained by the semiempirical QM/MM method are in Table 8. We see that the ab initio values are fairly well reproduced. In the case of methane εN is very small, so apparently the related LJ terms are practically of no consequence in comparison with the electrostatic interactions. The azobenzenemethanol attractive interactions, either at the ab initio or at the semiempirical QM/MM level, are weaker than those between two methanol molecules: the OPLS force field yields a well depth of 6.04 kcal/mol for the methanol dimer. This explains why hydrogen bonding does not play any special role in azobenzene solvation.66 For a sake of simplicity, we determined the interaction potentials and the QM/MM parameters for the ground state only. However, in our semiempirical treatment the electrostatic interactions between the QM and the MM subsystems depend on the QM charge distribution and are therefore state specific; thus, for instance, the spectral shifts due to solvation could in principle be reproduced. In the case of azobenzene, such effects are rather unimportant, as seen by comparing the absorption spectra in different solvents,90,91 probably because of the relative weakness of the interactions with H-bond donors already discussed. In summary, the isolated molecule potential energy surfaces are only slightly modified by solvation, and the important photophysical effects observed both experimentally3,7 and in simulations66 are essentially due to solvent caging.
7. APPENDIX A: SEMIEMPIRICAL CI ENERGY GRADIENT We extend here the Z-vector technique for evaluation of CI analytical gradients in semiempirical framework with floating occupation molecular orbitals (FOMO), outlined in a previous paper,20 to the more general case where the FOMOs are obtained with a different set of semiempirical parameters with respect to the one used in the CI. These two sets of semiempirical parameters will be called the CI set and the SCF set. According to our previous formulation,20 we assume that only the CI-active orbitals are allowed to have floating occupation numbers. In this appendix the molecular orbitals are labeled as follows a,b,... = any orbital i,j,... = active orbitals d,e,... = inactive orbitals v,w,... = virtual orbitals In the following, we will use a tilde to denote quantities evaluated with the CI set of semiempirical parameters. We limit here our description to the case where the variable semiempirical parameters are of the one-electron type. Therefore, if h is the one~ = h + T, the same electron Hamiltonian matrix, defining T as h ~ = F + T. FOMOs relation holds unaltered for the Fock matrix: F ~. We have then for the SCF energy61 are eigenvectors of F, not F ~HF ¼ EHF þ E
∑a Taa Oa
ð11Þ
where Oa is the occupation number of molecular orbital ja. As shown in our previous paper,20 the CI energy of state n becomes ~nCI ¼ EnCI þ 2 E 108
∑d Tdd þ ∑ij TijΔnij
ð12Þ
dx.doi.org/10.1021/jp208574q |J. Phys. Chem. A 2012, 116, 98–110
The Journal of Physical Chemistry A
ARTICLE
where Δn is the one-electron spinless CI density matrix for state n. The gradient of the CI energy can be partitioned in a “static” and a “response” term, the last one containing derivatives of the MO coefficients, energies, or occupation numbers. The static term is immediately obtained from the above equation; we focus then on the response term. The response derivative of T is ∂T ∂Cþ ¼ Bα T TBα , Bα ¼ C ð13Þ ∂Rα ∂Rα
(6) Satzger, H.; Root, C.; Braun, M. J. Phys. Chem. A 2004, 108, 6265–6271. (7) Chang, C.-W.; Lu, Y.-C.; Wang, T.-T.; Diau, E. W.-G. J. Am. Chem. Soc. 2004, 126, 10109–10118. (8) Lu, Y.-C.; Diau, E. W.-G.; Rau, H. J. Phys. Chem. A 2005, 109, 2090–2099. (9) Stuart, C. M.; Frontiera, R. R.; Mathies, R. A. J. Phys. Chem. A 2007, 111, 12072–12080. (10) Cattaneo, P.; Persico, M. Phys. Chem. Chem. Phys. 1999, 1, 4739–4743. (11) Åstrand, P.-O.; Ramanujam, P. S.; Hvilsted, S.; Bak, K. L.; Sauer, S. P. A. J. Am. Chem. Soc. 2000, 122, 3482–3487. (12) Ishikawa, T.; Noro, T.; Shoda, T. J. Chem. Phys. 2001, 115, 7503–7512. (13) Fliegl, H.; K€ohn, A.; H€attig, C.; Ahlrichs, R. J. Am. Chem. Soc. 2003, 125, 9821–9827. (14) Diau, E. W.-G. J. Phys. Chem. A 2004, 108, 950–956. (15) Gagliardi, L.; Orlandi, G.; Bernardi, F.; Cembran, A.; Garavelli, M. Theor. Chem. Acc. 2004, 111, 363–372. (16) Cembran, A.; Bernardi, F.; Garavelli, M.; Gagliardi, L.; Orlandi, G. J. Am. Chem. Soc. 2004, 126, 3234–3243. (17) Ciminelli, C.; Granucci, G.; Persico, M. Chem.—Eur. J. 2004, 10, 2327–2341. (18) Barada, D.; Itoh, M.; Yatagai, T. J. Appl. Phys. 2004, 96, 4204–4210. (19) Toniolo, A.; Ciminelli, C.; Persico, M.; Martínez, T. J. J. Chem. Phys. 2005, 123, 234308/1–10. (20) Ciminelli, C.; Granucci, G.; Persico, M. J. Chem. Phys. 2005, 123, 174317/1–10. (21) Crecca, C. R.; Roitberg, A. E. J. Phys. Chem. A 2006, 110, 8188–8203. (22) Granucci, G.; Persico, M. Theor. Chem. Acc. 2007, 117, 1131–1143. (23) Cusati, T.; Granucci, G.; Persico, M.; Spighi, G. J. Chem. Phys. 2008, 128, 194312/1–9. (24) Creatini, L.; Cusati, T.; Granucci, G.; Persico, M. Chem. Phys. 2008, 347, 492–502. (25) Ciminelli, C.; Granucci, G.; Persico, M. Chem. Phys. 2008, 349, 325–333. (26) Conti, I.; Garavelli, M.; Orlandi, G. J. Am. Chem. Soc. 2008, 130, 5216–5230. (27) Wang, L.; Wu, W.; Yi, C.; Wang, X. J. Mol. Graphics Modell. 2009, 27, 792–796. (28) B€ ockmann, M.; Doltsinis, N. L.; Marx, D. Phys. Rev. E 2008, 78, 036101/1–4. (29) B€ ockmann, M.; Doltsinis, N. L.; Marx, D. J. Phys. Chem. A 2010, 114, 745–754. (30) Ootani, Y.; Satoh, K.; Nakayama, A.; Noro, T.; Taketsugu, T. J. Chem. Phys. 2009, 131, 194306/1–10. (31) Tiberio, G.; Muccioli, L.; Berardi, R.; Zannoni, C. ChemPhysChem 2010, 11, 1018–1028. (32) Carstensen, O.; Sielk, J.; Shoenborn, J. B.; Granucci, G.; Hartke, B. J. Chem. Phys. 2010, 133, 124305/1–11. (33) Weingart, O.; Lan, Z.; Koslowski, A.; Thiel, W. J. Phys. Chem. Lett. 2011, 2, 1506–1509. (34) Cojocariu, C.; Rochon, P. Pure Appl. Chem. 2004, 76, 1479–1497. (35) Yager, K. G.; Barrett, C. J. J. Photochem. Photobiol. A 2006, 182, 250–261. (36) Browne, W. R.; Feringa, B. L. Nat. Nanotechnol. 2006, 1, 25–35. (37) Camorani, P.; Fontana, M. P. Phys. Rev. E 2006, 73, 011703/ 1–6. (38) Vecchi, I.; Arcioni, A.; Bacchiocchi, C.; Tiberio, G.; Zanirato, P.; Zannoni, C. J. Phys. Chem. B 2007, 111, 3355–3362. (39) Ambrosio, A.; Camposino, A.; Maddalena, P.; Patane, S.; Allegrini, M. J. Microsc. 2008, 229, 307–312. (40) Rais, D.; Zakrevskyy, Y.; Stumpe, J.; Nespøurek, S.; Sedlakova, Z. Opt. Mater. 2008, 30, 1335–1342.
response
where C is the matrix of MO coefficients and Bα is the antisymmetric matrix collecting their derivatives. We have then ~nCI ∂E ∂EnCI ¼ þ 4 ∑ Bαda Tad þ 2 ∑ Bαia Taj Δnij ∂Rα ∂Rα aij ad response
response
ð14Þ The above equation can be rewritten in the compact form ~nCI ∂E ~ ab ¼ X αab Q ð15Þ ∂Rα aeb
∑
response
where the matrix Xα, collecting the CPHF unknowns, is equal to Bα supplemented, on the diagonal, by the derivatives of the (untilded) orbital SCF energies,20,92 and ~ ii ¼ Q ii Q ~ Q ij ¼ Q ij þ 2
∑k ðTjk Δnik TikΔnjk Þ ~ di ¼ Q di þ 2 ∑ Tdj Δnij 4Tid Q j ~ iv ¼ Q iv þ 2 ∑ Tvj Δnij Q j
ð16Þ
~ dv ¼ 4Tdv Q The untilded Q terms are defined as in Ciminelli et al.20 As Xα is unaffected by the CI parameter set, the Z-vector CPHF formalism ~ already formulated20 can be applied just by replacing Q with Q without further modifications.
’ AUTHOR INFORMATION Corresponding Author
*E-mail:
[email protected].
’ ACKNOWLEDGMENT This work was supported by grants of the Italian MIUR (PRIN project 2006030944) and of the University of Pisa. ’ REFERENCES (1) Tamai, N.; Miyasaka, H. Chem. Rev. 2000, 100, 1875–1890. (2) Fujino, T.; Arzhantsev, S. Yu.; Tahara, T. Bull. Chem. Soc. Jpn. 2002, 75, 1031–1040. (3) Lu, Y.-C.; Chang, C.-W.; Diau, E. W.-G. J. Chin. Chem. Soc. 2002, 49, 693–701. (4) Schultz, T.; Quenneville, J.; Levine, B.; Toniolo, A.; Martínez, T. J.; Lochbrunner, S.; Schmitt, M.; Schaffer, J. P.; Zgierski, M. Z.; Stolow, A. J. Am. Chem. Soc. 2003, 125, 8098–8099. (5) Satzger, H.; Sp€orlein, S.; Root, C.; Wachtveitl, J.; Zinth, W.; Gilch, P. Chem. Phys. Lett. 2003, 372, 216–223. 109
dx.doi.org/10.1021/jp208574q |J. Phys. Chem. A 2012, 116, 98–110
The Journal of Physical Chemistry A
ARTICLE
(78) Briquet, L.; Vercauteren, D. P.; Perpete, E. A.; Jaquemin, D. Chem. Phys. Lett. 2006, 417, 190–195. (79) Tsuji, T.; Takashima, H.; Takeuchi, H.; Egawa, T.; Konaka, S. J. Phys. Chem. A 2001, 105, 9347–9353. (80) Brown, E. V.; Granneman, G. R. J. Am. Chem. Soc. 1975, 97, 621–627. (81) Mostad, A.; Romming, C. Acta Chem. Scand. 1971, 25, 3561–3568. (82) Pacher, T.; Cederbaum, L. S.; K€oppel, H. Adv. Chem. Phys. 1993, 84, 293–391. (83) Persico, M. Electronic diabatic states: definition, computation and applications. In Encyclopedia of Computational Chemistry; Schleyer, P. v. R., Allinger, N. L., Clark, T., Gasteiger, J., Kollman, P. A., Schaefer III, H. F., Schreiner, P. R., Eds.; Wiley: Chichester, 1998; p 852. (84) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery Jr., J. A.; Vreven, T.; Kudin, K. N.; Burant, J. C.; et al. Gaussian 03; Gaussian Inc., Wallingford CT, 2004. (85) Boys, S. F.; Bernardi, F. Mol. Phys. 1970, 19, 553–557. (86) Persico, M.; Granucci, G.; Inglese, S.; Laino, T.; Toniolo, A. J. Mol. Struct. (THEOCHEM) 2003, 621, 119–126. (87) Ponder, J. W. TINKER 4.2; Washington University School of Medicine: St. Louis, MO, 2004; http://dasher.wustl.edu/tinker. (88) Jorgensen, W. L.; Tirado-Rives, J. J. Am. Chem. Soc. 1988, 110, 1657–1666. (89) Granucci, G.; Persico, M. J. Comput. Chem. 2011, 32, 2690–2696. (90) Ronayette, J.; Arnaud, R.; Lebourgeois, P.; Lemaire, J. Can. J. Chem. 1974, 52, 1848–1857. (91) N€agele, T.; Hoche, R.; Zinth, W.; Wachtveitl, J. Chem. Phys. Lett. 1997, 272, 489–495. (92) Errata corrige in Ciminelli et al.:20 eq A17 should read Qiv = qiv.
(41) Ferri, V.; Elbing, M.; Pace, G.; Dickey, M. D.; Zharnikov, M.; Samor, P.; Rampi, M. A. Angew. Chem., Int. Ed. 2008, 47, 4307–4309. (42) Sasaki, T.; Tour, J. M. Org. Lett. 2008, 10, 897–900. (43) Oki, K.; Nagasaka, Y. Colloid Surf. A: Physicochem. Eng. Aspects 2009, 333, 182–186. (44) Domenici, V.; Ambrozic, G.; Copic, M.; Lebar, A.; Drevensek igon, M. Polymer 2009, Olenik, I.; Umek, P.; Zalar, B.; Zupancic, B.; Z 50, 4837–4844. (45) Devetak, M.; Zupancic, B.; Lebar, A.; Umek, P.; Zalar, B.; igon, M.; C opic, M.; Drevensek-Olenik, I. Domenici, V.; Ambrozic, G.; Z Phys. Rev. E 2009, 80, 050701/1–4. (46) Yoshino, T.; Kondo, M.; Mamiya, J.; Kinoshita, M.; Yu, Y.; Ikeda, T. Adv. Mater. 2010, 22, 1361–1363. (47) Leforestier, C. J. Chem. Phys. 1978, 68, 4406–4410. (48) Car, R.; Parrinello, M. Phys. Rev. Lett. 1985, 55, 2471–2474. (49) Guidon, M.; Schiffmann, F.; Hutter, J.; VandeVondele, J. J. Chem. Phys. 2008, 128, 214104/1–15. (50) Vreven, T.; Bernardi, F.; Garavelli, M.; Olivucci, M.; Robb, M. A.; Schlegel, H. B. J. Am. Chem. Soc. 1997, 119, 12687–12688. (51) Kaledin, A. L.; Morokuma, K. J. Chem. Phys. 2000, 113, 5750–5762. (52) Ben-Nun, M.; Martínez, T. J. Adv. Chem. Phys. 2002, 121, 439–512. (53) Barbatti, M.; Granucci, G.; Persico, M.; Ruckenbauer, M.; Vazdar, M.; Eckert-Maksic, M.; Lisckha, H. J. Photochem. Photobiol. A 2007, 190, 228–240. (54) Levine, B. G.; Martínez, T. J. J. Phys. Chem. A 2009, 113, 12815–12824. (55) Smith, B. R.; Bearpark, M. J.; Robb, M. A.; Bernardi, F.; Olivucci, M. Chem. Phys. Lett. 1995, 242, 27–32. (56) Granucci, G.; Persico, M.; Toniolo, A. J. Chem. Phys. 2001, 114, 10608–10615. (57) Fabiano, E.; Keal, T. W.; Thiel, W. Chem. Phys. 2008, 349, 334–347. (58) Weber, W.; Thiel, W. Theor. Chem. Acc. 1999, 103, 495–506. (59) Otte, N.; Scholten, M.; Thiel, W. J. Phys. Chem. A 2007, 111, 5751–5751. (60) Stewart, J. J. P. J. Mol. Mod. 2007, 13, 1173–1213. (61) Granucci, G.; Toniolo, A. Chem. Phys. Lett. 2000, 325, 79–85. (62) Toniolo, A.; Persico, M.; Pitea, D. J. Chem. Phys. 2000, 112, 2790–2997. (63) Collaveri, C.; Granucci, G.; Persico, M. J. Chem. Phys. 2001, 115, 1251–1263. (64) Slavícek, P.; Martínez, T. J. J. Chem. Phys. 2010, 132, 234102/ 1–10. (65) Stewart, J. J. P. MOPAC 2000 and MOPAC 2002; Fujitsu Ltd.: Tokio, Japan, 2000 and 2002. (66) Cusati, T.; Granucci, G.; Persico, M. J. Am. Chem. Soc. 2011, 133, 5109–5123. (67) Toniolo, A.; Thompson, A. L.; Martínez, T. J. Chem. Phys. 2004, 304, 133–145. (68) Nelder, J.; Mead, R. Comput, J, 1965, 7, 308–313. (69) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes in Fortran 77; Cambridge U. P, 2001. (70) Andersson, J.-Å.; Petterson, R.; Tegner, L. J. Photochem. 1982, 20, 17–32. (71) Biancalana, A.; Campani, E.; Di Domenico, G.; Gorini, G.; Iacoponi, A.; Masetti, G. Spectrochim. Acta A 1999, 55, 2883–2887. (72) Shashoua, V. E. J. Am. Chem. Soc. 1960, 82, 5505–5506. (73) Adamson, A. W.; Vogler, A.; Kunkely, H.; Wachter, R. J. Am. Chem. Soc. 1978, 100, 1298–1300. (74) Cimiraglia, R.; Hofmann, H.-J. Chem. Phys. Lett. 1994, 217, 430–435. (75) Chen, P. C.; Chieh, Y. C. J. Mol. Struct. (THEOCHEM) 2003, 624, 191–200. (76) Tiago, M. L.; Ismail-Beigi, S.; Louie, S. G. J. Chem. Phys. 2005, 122, 094311/1–7. (77) F€uchsel, G.; Klamroth, T.; Dokic, J.; Saalfrank, P. J. Phys. Chem. B 2006, 110, 16337–16345. 110
dx.doi.org/10.1021/jp208574q |J. Phys. Chem. A 2012, 116, 98–110
