1、二苯乙烯光致顺反异构化反应研究论文Article* E-mail: cyzhumail.nctu.edu.twReceived April 21, 2012; published August 6, 2012.Project supported by the National Natural Science Foundation of China (Nos. 21003100, 21033001, 21103136, 21173166. 项目受国家自然科学基金(Nos. 21003100, 21033001, 21103136, 21173166资助.化 学 学 报ACTA CHIMICA S
2、INICA改进的半经典动力学模拟二苯乙烯光致顺反异构化反应雷依波a ,b 朱超原*,b 文振翼a,b 林聖聖b(a 合成与天然功能分子化学教育部重点实验室 西北大学化学与材料科学学院 西安 710069(b 西北大学现代物理研究所 西安 710069摘要 发展了一种改进的半经典动力学模拟方法, 并将其程序化用于气相二苯乙烯光致顺反异构化反应的机理研究. 新的方法不仅采用e 指数模型改进了原有Zhu-Nakamura 理论中计算电子非绝热跃迁几率的计算方法, 而且将约束哈密顿方法用于限制性分子动力学模拟过程中. 计算结果表明, 采用此方法得到的统计平均的量子产率及反应机理与以前的实验与理论结果吻
3、合较好, 从而可以应用于全量子动力学方法无法进行的大分子体系的动力学研究. 关键词 改进的半经典动力学模拟; 约束哈密顿系统; Zhu-Nakamura 理论; 二苯乙烯顺反异构化; 二维解析势能面New Implementation of Semi-classical Dynamic Simulation on the Photoisomerization of cis- and trans-Isomers of Free StilbeneLei, Yibo a ,b Zhu, Chaoyuan *,b Wen, Zhenyi a ,b Lin, Sheng-Hsien b(a Key La
4、boratory of Synthetic and Natural Functional Molecule Chemistry of Ministry of Education , The College of Chemis-try & Materials Science , Shaanxi Key Laboratory of Physico-Inorganic Chemistry , Northwest University , Xian 710069(b Institute of Modern Physics , Northwest University , Xian 710069 Abs
5、tract New implementation of semi-classical trajectory surface hopping dynamic simulation has been developed and applied to the photoisomerization of cis- and trans- isomers on the gas phase. This method not only uses the exponential model to the modification of the originally analytical non-adiabati
6、c transition probability formula, but also involves the con-strained Hamiltonian system into the constrained molecular dynamic simulation. Two-dimensional potential energy surfaces of ground S 0 and excited S 1 states are constructed analytically fitting to ab initio calculations in terms of torsion
7、 angle and one dihedral angle around the central ethylenic C =C bond as variables, and the other internal coordinates are all fixed at configuration of one-bond flip conical intersection. The analytical PESs are quite accurate and the mean absolute error is less than 2.4 kcal mol -1, and much less t
8、han 1.0 kcal mol -1 around conical intersection region. A straight seam line is found on potential energy surfaces that simply separates the cis-area with the trans-area. The constrained Hamiltonian system is em-ployed to run trajectories in the Cartesian coordinate system and surface hopping in ter
9、ms of the two internal dihedral angles. Typical trajectories are found in which the torsion angle changes monotonically for both cis- to trans- and trans- to cis- isomerizations. This is an exact picture of one-bond flip mechanism of photoisomerization around the conical intersection. Quantum yield
10、for trans- to cis- isomerization is simulated as 60.45% in very good agreement with experimental value 55.0%, while quantum yield for cis- to trans- isomerization is simulated as 42.3% in comparison with experimental value 35.0%. As the S 1 energy in local minimum of cis-area is higher than that in
11、trans-area, and thus cis- to trans- isomerization is quite possible to access to another Hula-Twist conical intersection. These simulation results demonstrate that the computed cumulative quantum yield and reaction mechanism are consistent with the previously experimental and theoretical results. Th
12、is means that the present trajectory surface hopping method would be good at the dynamic simulation on the large system with or without constraint Hamiltonian in comparison with the quantum molecular dynamics.Keywords new implementation of semi-classical dynamic simulation; constrained Hamiltonian s
13、ystem; Zhu-Nakamura theory; isomerization of cis- and trans-stilbene; two dimensional analytical potential energy surfaces1 引言众所周知, 二苯乙烯在光照下可以从其顺式构型转变为反式构型, 也可以从反式结构异构化到顺式结构, 其中反式二苯乙烯分子构型示于图 1. 反应过程中电子首先受光照从基态激发到激发态, 再经过无辐射跃迁回到分子基态. 一般情况下, 电子跃迁主要集中在第一激发态S 1和基态S 0之间1, 其中S 1对应电子从S 0的最高占据分子轨道(HOMO到最低未占
14、据轨道(LUMO的跃迁.大量的实验数据表明光照下二苯乙烯从顺式转换到反式和从反式到顺式的量子产率分别约为35%和55%(或52%25, 其中顺式到反式产率较低的原因是有DOI: 10.6023/A1204013910%的顺式二苯乙烯经过环化反应生成一个副产物4a, 4b-二氢菲(DHP68. 理论分析认为反应过程很有可能经过一个能量相对较低的圆锥相交点(标记为OBF-CI1, 当受激反应物靠近此相交点时, 电子就从激发态跃迁到达产物或者返回到基态. 此前的理论和实验报道倾向于反应坐标主要由苯环及氢原子绕着中心乙烯双键的旋转所决定9, 具体对应图1中的两个二面角D1和D2. 而Fu 等认为此顺反
15、异构化过程中还应该考虑苯环自身转动的影响, 此结论得到了一些实验的证实1012. 以上这些反应机理大多数只是依靠对圆锥相交点的分析. 为了更好地研究反应的动力学过程, 豆育升及其合作者通过实时激光诱导的动力学模拟对反应机理进行了研究, 模拟结果证实反应过程中描述苯环或氢原子绕中心烯键扭转的二面角的变化较大, 而苯环自身转动的二面角的变化较小13. 此动力学模拟过程中只考察了一个典型反应轨迹的运动, 没有考虑动力学中 的相效应, 因此无法得到反应的统计平均的量子产率. 另外一维和两维的势能面也被构建用于研究反应的机理, 但其中有些势能面只考虑苯环绕中心烯键扭转14; 有些通过实验数据拟合得到15
16、; 有些势能面虽然通过从头算方法计算得到, 但所选反应坐标非独立坐标, 无法用于动力学研究1.图1 反式二苯乙烯的几何构型图Figure 1 Structure and standard numbering of trans-stilbene鉴于此, 本文首先采用态平均多组态的自洽场方法(SA-CASSCF在6-31G 的基组水平上构建新的二维势能面16,17, 所选独立变量为图1所述二面角D1和D2的线性组合二面角DD1和DD2. 其中DD1为D1和D2的平均值(D2+D1/2; DD2表示D2相对D1的扭转大小, 即(D2-D1/2. 为了简化计算量, 我们拟合了势能面的解析表达式, 并用
17、于提供经典轨迹中原子核运动所需的力. 由于此分子体系较大, 无法构建3N -6维的全维势能面(其中N 为原子数, 因此本文采用约束Hamilton 系统限制部分核坐标的运动, 例如在动力学模拟过程中每个苯环将始终作为一个刚体运动. 在此基础上将最新改进的Trajectory Surface Hopping (TSH方法用于半经典动力学模拟二苯乙烯的光致顺反异构化过程中. TSH 方法由Tully 和Preston 18首先提出, 其基本思想是反应轨迹总是在单一势能面上演化, 电子的非绝热跃迁几率可以由Laudau-Zener 公式或数值求解含时耦合方程得到. 20世纪90年代, 朱超原与Nak
18、amura 共同提出的Zhu-Nakamura 理论简化了非绝热跃迁几率的计算1925. 此理论分别采用Landau-Zener 和nonadiabatic tunneling 两种模型推出两种类型的几率计算公式19,20. 采用哪一种模型取决于非绝热区两个电子态的势能的梯度的符号是否相同, 前者对应同向, 后者对应反向, 因此在计算跃迁几率时需要首先确定势能面的类型. 最近, 我们基于此理论进一步采用e 指数模型推出了非绝热跃迁几率的新公式, 避免了预先确定势能面类型的工作. 下面将简要介绍新的计算方法和程序化过程, 并详细报道基于此方法所做的动力学模拟的结果与讨论.2 理论背景如前所述,
19、TSH 方法中核运动采用约束Hamilton 系统进行经典计算, 而电子是否发生无辐射跃迁取决于非绝热跃迁几率的大小. 下面我们简要介绍一下这两部分内容的理论背景. 2.1 约束系统Hamilton约束系统Hamilton 量:23(,(2CN Ni k k i ki p H q p V q g q m =+(1应该满足约束条件(0k g q =(2从而保持体系能量守恒. 此Hamilton 系统所对应的正则方程可以描述为2629:(,i i ii p H q p q p m = (3(,(C N uc cki i i k ki i i g q H q p V q p f f q q q =-
20、=+=-(4 其中m i , q i 和p i 分别对应原子i 的质量、坐标和动量. g k (q 和k (t 分别对应第k 个约束条件(包括键长、键角和二面角约束及所对应的拉格朗日乘子. 式(4中uc i f 所对应的是体系势场所提供的力, 而c i f 是用来约束核运动的力, 其贡献来自于约束方程关于核坐标的微分(/k i g q q . 我们采用拉格朗日乘子法求解满足约束方程式(2的k (t . 首先确定t 时刻的核坐标及动量, 然后数值积分式(3得到t +t 时刻的坐标271212(Cuc c i i i i N uc k i i kkiq t t q t t m t f t g q
21、q t t m t q +=+=+- (5其中(uc i q t t +是t +t 时刻未受约束的核坐标. 将q i (t +t 带入式(2从而构造包括键长、键角和二面角约束的非线性方程组. 求解这些关于k (t 的非线性方程组需要采用Powells Dog Leg 方法30. 此方法通过如下公式00(|k k g J g =-(6逐步更新k , 直至2-范数k g 小于预先设定的阈值. 其中J g 为约束方程g k (q 所对应的雅克比矩阵. 此方法已经收录在免费程序包MINPACK 31,32中, 只需将预先得到的g k (q 和J g 输入此程序中即可得到所需求解的k (t . 将其代入
22、式(3和(4应用数值积分方法四阶龙格-库塔方法(RK433,34, 即可求解t +t 时刻体系的约束坐标及动量. 将此作为起始可以得到经过下一个时间步长t 后体系的坐标及动量. 这样沿着固定时间步长数值运算核坐标及其动量时间的过程对应体系的动力学模拟过程. 2.2 非绝热跃迁几率在一维透热模型的势能曲线上, 体系的无辐射跃迁与两个参数有关, 分别是有效耦合常数a 2和有效碰撞能b 2 24. 它们分别表示为:22213(28xF F F a m V =-= (7a221(2x xF F b E E FV -=- (7b其中m 是体系的约化质量, F i (i =1, 2是透热势能的斜率, F
23、. V x 代表透热势能耦合项, E x 是透热势能. 一般情况下, 电子跃迁是在绝热表象下进行. 基于e 指数模型, 参数a 2和b 2在此表象可以表示为:0222301(16(x x a W x W x m x E +=+ (8a 22xE E b E -=(8b其中W +(x 和W (x 分别是能量较高和较低的电子态的能量, 此时的(2x W x W x E +-+=, 而0E =(2W x W x +-. 其中式(8b中E 为体系能量. 基于此模型, 非绝热跃迁几率可以表示为:0ZN 1exp 4exp (x x P ab W x W x +=-=+(9众所周知, 电子的无辐射跃迁沿着
24、非绝热耦合向量最强的方向进行. 为了计算电子非绝热跃迁几率P ZN 需要首先确定此向量的方向. 对于大分子体系来说, 计算此向量非常耗时, 而实践证明势能面上的接缝线与非绝热耦合向量基本垂直35, 所以找到此接缝线就意味着确定了电子跃迁的方向. 朱超原等报道称此接缝线可以在运行经典轨迹之前预先确定, 并将体系在此线所处的区域分为三部分, 分别为透热区、绝热区和非绝热区. 划分的标准由参数a 2决定, 当a 2103时, 体系处于透热区, 此时确定跃迁发生; 当a 210-2, 体系处于绝热区, 此时无跃迁发生; 当103a 210-2, 体系处于非绝热区, 此时电子跃迁发生与否取决于此时P Z
25、N 是否大于一个0和1之间随机数35.由式(8a可知, 计算a 2需要势能面的信息. 如图2所示, 我们在SA-CASSCF/6-31G 水平上16,17 构建了二维势能面, 其中所选活性空间包括2个活性电子及2个活性轨道HOMO 和LUMO, 势能变化沿着二面角DD1和DD2进行. 采用线性最小二乘法36拟合得到了此势能面的解析表达式(10.图2 二苯乙烯围绕圆锥相交点OBF-CI 的基态与第一激发态的二维解析势能面,其中变量为组合的二面角DD1和DD2Figure 2 Analytical two-dimensional PESs around OBF-CI for the ground
26、state and the first excited state with respect to the combined inter-nal coordinates DD1 and DD222010101234567891011121314(,exp(x (cos(/2cos(/2cos(cos(cos(2cos(2cos(3cos(3sin(/2sin(2/2sin(2/2sin(sin(2sin(2W x y c a x b y y c x c y c x c y c x c y c x c y c x y c x y c x y c x y c x y c x y =-+151617
27、1819202122232425262728sin(22sin(3sin(3sin(32sin(23sin(33cos(/2cos(2/2cos(2/2cos(cos(2cos(2cos(22cos(3c x y c x y c x y c x y c x y c x y c x y c x y c x y c x y c x y c x y c x y c x y c +29303132cos(3cos(32cos(23cos(33x y c x y c x y c x y +(10其中, x 0=-103.3和y 0=43.3为相交点OBF-CI 的DD1和DD2值, 而其他的系数a 1,
28、 b 1和c 0c 32在表1中. 解析势能面与计算势能面的对比表明, 无论S 0还是S 1的平均偏差都小于 2.4 kcal mol -1, 而且最小偏差出现在相交点OBF-CI 附近, 此处平均偏差小于1.0 kcal mol -1, 即最重要的相交区域的计算结果也最为精确, 这为此后的动力学计算打下良好基础.表1 二苯乙烯围绕圆锥相交点OBF-CI 关于DD1和DD2的基态与第一激发态的二维解析势能面的拟合参数(其中c 0c 32单位为eVTable 1 The fitted coefficients of two-dimensional S 0 and S 1 PESs around
29、OBF-CI with respect to the combined internal angles DD1 and DD2 (c 0c 32 in units of eV Index a 1 b 1c 0c 1c 2c 3 c 4c 5c 6S 0 35 3.5 0.495231 364.2436 166.2537 -63.2233-228.305 -12.7773 25.79321 S 1 39 3.9 -0.52219 256.2847 -18.8715-74.255 -85.6398 -4.51131 16.61385 Index c 7 c 8c 9c 10 c 11 c 12 c
30、 13 c 14c 15 S 0 -3.78713 -3.97288 9.467412 106.3829 -36.167 -72.9472 5.205437 -0.08421 -4.58729S 1 -1.65933 -2.44485 -219.158 171.3529 98.32809 -110.7899.156298 -4.96308 -1.1299 Index c 16 c 17 c 18c 19 c 20 c 21 c 22c 23 c 24 S 0 5.134825 3.857498 -6.20916 0.171133 2.090849 -538.156-56.4412 254.75
31、31 157.3086 S 1 0.059523 4.654567 -1.20414 -1.738630.658395 -277.52987.56788 112.7934 47.53832 Index c 25 c 26 c 27 c 28 c 29 c 30 c 31 c 32 S 0 17.19558 -58.19 -8.161 7.267131 4.132838 -3.08607 1.734321 0.088946 S 1 -7.42468-41.53027.9652671.253993.375719-0.42815-0.71119-0.09718结合解析势能面表达式(10与式(8a, 得到了如图3所示的a 2值. 由此图可知, 沿着DD1约等于-103.35的直线(即接缝线方向上, a 2都相对较大, 最大值在圆锥相交点OBF-CI 处. 当反应轨迹接近此线时, 通过此时的a 2判断体系所处的区域, 可以判断电子是否发生跃迁, 即电子跃迁只可能发生在此接缝线上, 而势能面上的其他区域无需
