二苯乙烯光致顺反异构化反应.docx

二苯乙烯光致顺反异构化反应.docx

《二苯乙烯光致顺反异构化反应.docx》由会员分享，可在线阅读，更多相关《二苯乙烯光致顺反异构化反应.docx（19页珍藏版）》请在冰豆网上搜索。

二苯乙烯光致顺反异构化反应

研究论文

Article

*E-mail:

cyzhu@mail.nctu.edu.tw

ReceivedApril21,2012;publishedAugust6,2012.

ProjectsupportedbytheNationalNaturalScienceFoundationofChina（Nos.21003100,21033001,21103136,21173166.项目受国家自然科学基金（Nos.21003100,21033001,21103136,21173166资助.

化学学报

ACTACHIMICASINICA

改进的半经典动力学模拟二苯乙烯光致顺反异构化反应

雷依波a,b朱超原*,b文振翼a,b林聖聖b

（a合成与天然功能分子化学教育部重点实验室西北大学化学与材料科学学院西安710069

（b西北大学现代物理研究所西安710069

摘要发展了一种改进的半经典动力学模拟方法,并将其程序化用于气相二苯乙烯光致顺反异构化反应的机理研究.新的方法不仅采用e指数模型改进了原有Zhu-Nakamura理论中计算电子非绝热跃迁几率的计算方法,而且将约束哈密顿方法用于限制性分子动力学模拟过程中.计算结果表明,采用此方法得到的统计平均的量子产率及反应机理与以前的实验与理论结果吻合较好,从而可以应用于全量子动力学方法无法进行的大分子体系的动力学研究.关键词改进的半经典动力学模拟;约束哈密顿系统;Zhu-Nakamura理论;二苯乙烯顺反异构化;二维解析势能面

NewImplementationofSemi-classicalDynamicSimulationonthePhotoisomerizationofcis-andtrans-IsomersofFreeStilbene

Lei,Yiboa,bZhu,Chaoyuan*,bWen,Zhenyia,bLin,Sheng-Hsienb

（aKeyLaboratoryofSyntheticandNaturalFunctionalMoleculeChemistryofMinistryofEducation,TheCollegeofChemis-try&MaterialsScience,ShaanxiKeyLaboratoryofPhysico-InorganicChemistry,NorthwestUniversity,Xi'an710069

（bInstituteofModernPhysics,NorthwestUniversity,Xi'an710069AbstractNewimplementationofsemi-classicaltrajectorysurfacehoppingdynamicsimulationhasbeendevelopedandappliedtothephotoisomerizationofcis-andtrans-isomersonthegasphase.Thismethodnotonlyusestheexponentialmodeltothemodificationoftheoriginallyanalyticalnon-adiabatictransitionprobabilityformula,butalsoinvolvesthecon-strainedHamiltoniansystemintotheconstrainedmoleculardynamicsimulation.Two-dimensionalpotentialenergysurfacesofgroundS0andexcitedS1statesareconstructedanalyticallyfittingtoabinitiocalculationsintermsoftorsionangleandonedihedralanglearoundthecentralethylenicC=Cbondasvariables,andtheotherinternalcoordinatesareallfixedatconfigurationofone-bondflipconicalintersection.TheanalyticalPESsarequiteaccurateandthemeanabsoluteerrorislessthan2.4kcal•mol-1,andmuchlessthan1.0kcal•mol-1aroundconicalintersectionregion.Astraightseamlineisfoundonpotentialenergysurfacesthatsimplyseparatesthecis-areawiththetrans-area.TheconstrainedHamiltoniansystemisem-ployedtoruntrajectoriesintheCartesiancoordinatesystemandsurfacehoppingintermsofthetwointernaldihedralangles.Typicaltrajectoriesarefoundinwhichthetorsionanglechangesmonotonicallyforbothcis-totrans-andtrans-tocis-isomerizations.Thisisanexactpictureofone-bondflipmechanismofphotoisomerizationaroundtheconicalintersection.Quantumyieldfortrans-tocis-isomerizationissimulatedas60.45%inverygoodagreementwithexperimentalvalue55.0%,whilequantumyieldforcis-totrans-isomerizationissimulatedas42.3%incomparisonwithexperimentalvalue35.0%.AstheS1energyinlocalminimumofcis-areaishigherthanthatintrans-area,andthuscis-totrans-isomerizationisquitepossibletoaccesstoanotherHula-Twistconicalintersection.Thesesimulationresultsdemonstratethatthecomputedcumulativequantumyieldandreactionmechanismareconsistentwiththepreviouslyexperimentalandtheoreticalresults.ThismeansthatthepresenttrajectorysurfacehoppingmethodwouldbegoodatthedynamicsimulationonthelargesystemwithorwithoutconstraintHamiltonianincomparisonwiththequantummoleculardynamics.

Keywordsnewimplementationofsemi-classicaldynamicsimulation;constrainedHamiltoniansystem;Zhu-Nakamuratheory;isomerizationofcis-andtrans-stilbene;twodimensionalanalyticalpotentialenergysurfaces

1引言

众所周知,二苯乙烯在光照下可以从其顺式构型转变为反式构型,也可以从反式结构异构化到顺式结构,其中反式二苯乙烯分子构型示于图1.反应过程中电子首先受光照从基态激发到激发态,再经过无辐射跃迁回到分子基态.一般情况下,电子跃迁主要集中在第一激

发态S1和基态S0之间[1],其中S1对应电子从S0的最高占据分子轨道（HOMO到最低未占据轨道（LUMO的跃迁.

大量的实验数据表明光照下二苯乙烯从顺式转换到反式和从反式到顺式的量子产率分别约为35%和55%（或52%[2~5],其中顺式到反式产率较低的原因是有

DOI:

10.6023/A12040139

10%的顺式二苯乙烯经过环化反应生成一个副产物4a,4b-二氢菲（DHP[6~8].理论分析认为反应过程很有可能经过一个能量相对较低的圆锥相交点（标记为OBF-CI[1],当受激反应物靠近此相交点时,电子就从激发态跃迁到达产物或者返回到基态.此前的理论和实验报道倾向于反应坐标主要由苯环及氢原子绕着中心乙烯双键的旋转所决定[9],具体对应图1中的两个二面角D1和D2.而Fuβ等认为此顺反异构化过程中还应该考虑苯环自身转动的影响,此结论得到了一些实验的证实[10~12].以上这些反应机理大多数只是依靠对圆锥相交点的分析.为了更好地研究反应的动力学过程,豆育升及其合作者通过实时激光诱导的动力学模拟对反应机理进行了研究,模拟结果证实反应过程中描述苯环或氢原子绕中心烯键扭转的二面角的变化较大,而苯环自身转动的二面角的变化较小[13].此动力学模拟过程中只考察了一个典型反应轨迹的运动,没有考虑动力学中

的相效应,因此无法得到反应的统计平均的量子产率.另外一维和两维的势能面也被构建用于研究反应的机理,但其中有些势能面只考虑苯环绕中心烯键扭转[14];有些通过实验数据拟合得到[15];有些势能面虽然通过从头算方法计算得到,但所选反应坐标非独立坐标,无法用于动力学研究[1].

图1反式二苯乙烯的几何构型图

Figure1Structureandstandardnumberingoftrans-stilbene

鉴于此,本文首先采用态平均多组态的自洽场方法

（SA-CASSCF在6-31G的基组水平上构建新的二维势能面[16,17],所选独立变量为图1所述二面角D1和D2的线性组合二面角DD1和DD2.其中DD1为D1和D2的平均值（D2+D1/2;DD2表示D2相对D1的扭转大小,即（D2-D1/2.为了简化计算量,我们拟合了势能面的解析表达式,并用于提供经典轨迹中原子核运动所需的力.由于此分子体系较大,无法构建3N-6维的全维势能面（其中N为原子数,因此本文采用约束Hamilton系统限制部分核坐标的运动,例如在动力学模拟过程中每个苯环将始终作为一个刚体运动.在此基础上将最新改进的TrajectorySurfaceHopping（TSH方法用于半经典动力学模拟二苯乙烯的光致顺反异构化过程中.TSH方法由Tully和Preston[18]首先提出,其基本思想是反应轨迹总是在单一势能面上演化,电子的非绝热跃迁几率可以由Laudau-Zener公式或数值求解含时耦合方程得到.20世纪90年代,朱超原与Nakamura共同提出的

Zhu-Nakamura理论简化了非绝热跃迁几率的计算[19~25].此理论分别采用Landau-Zener和nonadiabatictunneling两种模型推出两种类型的几率计算公式[19,20].采用哪一种模型取决于非绝热区两个电子态的势能的梯度的符号是否相同,前者对应同向,后者对应反向,因此在计算跃迁几率时需要首先确定势能面的类型.最近,我们基于此理论进一步采用e指数模型推出了非绝热跃迁几率的新公式,避免了预先确定势能面类型的工作.下面将简要介绍新的计算方法和程序化过程,并详细报道基于此方法所做的动力学模拟的结果与讨论.

2理论背景

如前所述,TSH方法中核运动采用约束Hamilton系统进行经典计算,而电子是否发生无辐射跃迁取决于非绝热跃迁几率的大小.下面我们简要介绍一下这两部分内容的理论背景.2.1约束系统Hamilton

约束系统Hamilton量:

2

3（,,（（2C

NN

ikkik

ipHqpVqgqmλλ∑∑=++

（1

应该满足约束条件

（0kgq=

（2

从而保持体系能量守恒.此Hamilton系统所对应的正则方程可以描述为[26~29]:

（,,iii

ipHqpqpmλ∂∂==（3

（（,,（CNucc

k

iiikk

iiigqHqpVqpffqqqλλ∂∂∂∂∂∂∑=-=+=--

（4其中mi,qi和pi分别对应原子i的质量、坐标和动量.gk（q和λk（t分别对应第k个约束条件（包括键长、键角和二面角约束及所对应的拉格朗日乘子.式（4中ucif所对应的是体系势场所提供的力,而cif是用来约束核运动的力,其贡献来自于约束方程关于核坐标的微分（/kigqq∂∂.我们采用拉格朗日乘子法求解满足约束方程式（2的λk（t.首先确定t时刻的核坐标及动量,然后数值积分式（3得到t+Δt时刻的坐标[27]

1212

（（（（（（（C

ucciiiiNuckiik

k

i

qttqttmtftgqqttmtqλ−−ΔΔΔ∂ΔΔ∂∑+=++=+-（5

其中（uciqttΔ+是t+Δt时刻未受约束的核坐标.将qi（t+Δt带入式（2从而构造包括键长、键角和二面角约束的非线性方程组.求解这些关于λk（t的非线性方程组需要

采用Powell’sDogLeg方法[30].此方法通过如下公式

00（||kkgJgλλλ===-

（6

逐步更新λk,直至2-范数kg小于预先设定的阈值τ.其中Jg为约束方程gk（q所对应的雅克比矩阵.此方法已经收录在免费程序包MINPACK[31,32]中,只需将预先得到的gk（q和Jg输入此程序中即可得到所需求解的λk（t.将其代入式（3和（4应用数值积分方法四阶龙格-库塔方法（RK4[33,34],即可求解t+Δt时刻体系的约束坐标及动量.将此作为起始可以得到经过下一个时间步长Δt后体系的坐标及动量.这样沿着固定时间步长数值运算核坐标及其动量时间的过程对应体系的动力学模拟过程.2.2非绝热跃迁几率

在一维透热模型的势能曲线上,体系的无辐射跃迁与两个参数有关,分别是有效耦合常数a2和有效碰撞能b2[24].它们分别表示为:

22

213

（28x

FFFamV=-=（7a

221（

（

2xx

FFbEEFV-=-

（7b

其中m是体系的约化质量,Fi（i=1,2是透热势能的斜率,F.Vx代表透热势能耦合项,Ex是透热势能.一般情况下,电子跃迁是在绝热表象下进行.基于e指数模型,参数a2和b2在此表象可以表示为:

02

22

3

01[（（]16（xxaWxWxmxE+−=∂

⎡⎤

⎢⎥∂Δ⎣

⎦==

+（8a2

2x

EEbEΔ-=

（8b

其中W+（x和W―（x分别是能量较高和较低的电子态的能量,此时的（（

2

xWxWxE+-+=

而0EΔ=

（（

2

WxWx

+--.其中式（8b中E为体系能量.基于此模

型,非绝热跃迁几率可以表示为:

0ZN1

exp4exp（（]xxPabWxWxπ−+−=⎛

⎞⎜⎟⎝⎠

⎡⎤⎢⎥

⎢⎥⎣⎦

=-=

+

（9

众所周知,电子的无辐射跃迁沿着非绝热耦合向量最强的方向进行.为了计算电子非绝热跃迁几率PZN需要首先确定此向量的方向.对于大分子体系来说,计算此向量非常耗时,而实践证明势能面上的接缝线与非绝热耦合向量基本垂直[35],所以找到此接缝线就意味着

确定了电子跃迁的方向.朱超原等报道称此接缝线可以在运行经典轨迹之前预先确定,并将体系在此线所处的区域分为三部分,分别为透热区、绝热区和非绝热区.划分的标准由参数a2决定,当a2>103时,体系处于透热区,此时确定跃迁发生;当a2<10-2,体系处于绝热区,此时无跃迁发生;当103≥a2≥10-2,体系处于非绝热区,此时电子跃迁发生与否取决于此时PZN是否大于一个0和1之间随机数[35].

由式（8a可知,计算a2需要势能面的信息.如图2所示,我们在SA-CASSCF/6-31G水平上[16,17]

构建了二维势能面,其中所选活性空间包括2个活性电子及2个活性轨道HOMO和LUMO,势能变化沿着二面角DD1和DD2进行.采用线性最小二乘法[36]拟合得到了此势能面的解析表达式（10.

图2二苯乙烯围绕圆锥相交点OBF-CI的基态与第一激发态的二维解析势能面,其中变量为组合的二面角DD1和DD2

Figure2Analyticaltwo-dimensionalPESsaroundOBF-CIforthegroundstateandthefirstexcitedstatewithrespecttothecombinedinter-nalcoordinatesDD1andDD2

22010101234567891011121314（,exp[（x（]cos（/2cos（/2cos（cos（cos（2cos（2cos（3cos（3sin[（/2]sin[（2/2]sin[（2/2]sin（sin（2sin（2Wxycaxbyycxcycxcycxcycxcycxycxycxycxycxycxy=----+++++

++++

+++++++++++1516171819202122232425262728sin（22sin（3sin（3sin（32sin（23sin（33cos[（/2]cos[（2/2]cos[（2/2]cos（cos（2cos（2cos（22cos（3cxycxycxycxycxycxycxycxycxycxycxycxycxycxyc+++++++++++++++++++++++++++++29303132cos（3cos（32cos（23cos（33

xycxycxycxy+++++++

（10

其中,x0=-103.3°和y0=43.3°为相交点OBF-CI的DD1

和DD2值,而其他的系数a1,b1和c0~c32在表1中.解析势能面与计算势能面的对比表明,无论S0还是S1的平均偏差都小于2.4kcal•mol-1,而且最小偏差出现在相交点OBF-CI附近,此处平均偏差小于1.0kcal•mol-1,即最重要的相交区域的计算结果也最为精确,这为此后的动力学计算打下良好基础.

表1二苯乙烯围绕圆锥相交点OBF-CI关于DD1和DD2的基态与第一激发态的二维解析势能面的拟合参数（其中c0~c32单位为eV

Table1Thefittedcoefficientsoftwo-dimensionalS0andS1PESsaroundOBF-CIwithrespecttothecombinedinternalanglesDD1andDD2（c0~c32inunitsofeVIndexa1b1

c0

c1

c2

c3c4

c5

c6

S0353.50.495231364.2436166.2537-63.2233-228.305-12.777325.79321S1393.9-0.52219256.2847-18.8715

-74.255-85.6398-4.5113116.61385Indexc7c8

c9

c10c11c12c13c14

c15S0-3.78713-3.972889.467412106.3829-36.167-72.94725.205437-0.08421-4.58729S1-1.65933-2.44485-219.158171.352998.32809-110.789

9.156298-4.96308-1.1299Indexc16c17c18

c19c20c21c22

c23c24S05.1348253.857498-6.209160.1711332.090849-538.156-56.4412254.7531157.3086S10.0595234.654567-1.20414-1.73863

0.658395-277.529

87.56788112.793447.53832Indexc25c26c27c28c29c30c31c32S017.19558-58.19-8.1617.2671314.132838-3.086071.7343210.088946S1-7.42468

-41.5302

7.965267

1.25399

3.375719

-0.42815

-0.71119

-0.09718

结合解析势能面表达式（10与式（8a,得到了如图3所示的a2值.由此图可知,沿着DD1约等于-103.35°的直线（即接缝线方向上,a2都相对较大,最大值在圆锥相交点OBF-CI处.

当反应轨迹接近此线时,通过此时的a2判断体系所处的区域,可以判断电子是否发生跃迁,即电子跃迁只可能发生在此接缝线上,而势能面上的其他区域无需