Gaussian软件应用——激发态计算

Gaussian软件应用——激发态计算 第九章 激发态计算   本章探讨激发态的计算.激发态是指分子体系的稳定的,高能级的电子结构状态. 比如在高能紫外可见光谱中的激发态的分子.激发态在化学的很多领域都有用, 包括光化学,电子光谱.   对于激发态分子的模拟是困难的,因为很难避免体系向基态变化.模拟激发态的   理论很少,在Gaussian中,提供的是结构相关(Configuration Interaction  approach),采用Hartree-Fock单置换,所以也被称为CI-Singles.这一方法被认为  是很多分子激发态的零级处理的有效方法.在定性上有比较好的结果,但在定量 方面,并不总是好的.   和Hartree-Fock方法一样,CI-Singles方法也不昂贵.   9.1 运行激发态计算   下面是激发态计算的关键词   CIS 进行CI-Silgles激发态计算,这是方法关键词,在其前面可以   加上R和U作为闭壳层分子和开壳层分子的标识.   CIS=(Root=n) 确定研究哪一个激发态,默认为1   CIS=(NStates=n) 确定研究的激发态的数量.默认为3   CIS=50-50 同时计算单重态和三重态.Triplets选项为计算单重态,默认   为单重态.当选择Triplets时,一定要设置NStates数量.   CIS=Read 从临时文件中读入初始猜测.用于利用前一个工作中的结果   进行计算,与Guess=Read和Geom=Check一起使用.   Density=Current 要求计算激发态的布局分析   Pop=Reg 要求较详细的布局分析结果,包括分子轨道系数   例9.1 文件e9_01 乙烯的激发态   本例计算乙烯的最低的四个激发态.   乙烯的激发态已经有了很好的研究,最低的四个如下   多重度 Modern标记 Mulliken标记 跃迁   1 3 B1u T pi --> pi*   2 3 B3u TR pi --> 3s   3 1 B3u R pi --> 3s   4 1 B1u V pi --> pi*   计算设置如下   # RCIS=(NStates=2, 50-50)/6-31+G(d) Test   由于即需要单重态,也需要三重态,所以设置50-50,每个状态计算两个.   注意这里使用了弥散基组,弥散基组在计算过渡态方面有好的结果.   输出文件中有   Excited State 1: Triplet-B1U 3.7768eV 328.27 nm f=0.000   多重度 对称性 激发能 振荡强度   8 --> 11 0.52952   跃迁轨道数 激发态波函数系数   8 --> 17 -0.45942   This State for optimization and/or second-order correction:   Total Energy, E(Cis) = -77.8969983928   激发态能量   Copying the Cisingles density for this state as the   1-particle RhoCI density.   计算结果是   1 2 3 4   对称性 3B1u 3B3u 1B3u 1B1u   能量 3.78 7.43 7.83 7.98   实验值 4.36 7.66 6.98 7.15   总的说计算的结果与实验值符合的不错,更加精确的结果就需要大的基组计算了.   9.2 激发态优化和频率分析   在Gaussian中可以进行激发态的优化和频率分析.步骤是,首先进行能量计算得到   激发态,然后在此结构基础上进行结构优化和频率分析   例9.2 文件e9_02 甲醛的激发态优化   下面是输入文件   %Chk=es_form   #T RCIS/6-31+G(d) Test   formaldehyde Excited States   0 1   Ground sate molecule specification   --Link1--   %Chk=se_form   %NoSave   #T RCIS(Root=1,Read)/6-31+G(d) Opt Freq Geom=Check Guess=Read Test   第一个步骤计算最低的三个激发态,第二步骤是利用其结果进行结构优化和   频率分析   在本例的计算中,得到的激发态出现了虚频,其结构显示碳原子要离开原子平面.   由于原来的分子结构中,所有原子都在同一平面上,所以,也得到一个共平面的   激发态结构.这样,就要设法产生一个正确的几个结构.(激发态的结构优化目前   没有包含内坐标冗余,所以在初始结构中定义一个平面结构,得到的激发态结构   就也是共平面的).我们将在练习9.2中继续讨论.   练习   练习9.1 文件9_01 亚甲基环丙烯的激发态   亚甲基环丙烯在二十世纪八十年代中期发现的,其紫外光谱有三个峰,   位置(nm) 对称性 能量(ev) 相对面积   309 1B2 4.01 0.2   242 1B1 5.12 0.01   206 1A1 6.02 1.5   该分析结构如下   H H   \ /   C1   ||   C2   / \   C3==C4   / \   H H   半方法的计算结果显示在1A1下面另一个1B1的激发态,最初的解释为其可能   被206峰掩盖了.   在MP2/6-31G(d)水平进行结构优化,然后进行激发态计算   # RCIS(NStates=5)/6-31+G(d) Density=All Test   得到了五个激发态,如下   状态 对称性 能量(eV) 强度   1 1B2 5.48 4.01 0.05   2 1B1 5.91 5.12 0.02   3 1A2 6.30 0.0   4 1B1 6.38 0.03   5 1A1 6.41 6.02 0.37   得到的一,二,五三个激发态都与实验符合.   对于第三个和第四个,第三个激发态的强度为零,自然观测不到,第四个激发态的   强度很弱,而且与第五个激发态位置相近,也被其掩盖了.   下面得到的是偶极矩和原子电荷   方法 偶极矩 C1 C2 C3 C4   SCF -2.39z -0.5 -0.2 -0.1 -0.1   CI 1-Particle 4.75z -0.004 -0.06 -0.4 -0.4   CI 2.56z -0.2 -0.2 -0.3 -0.3   两个激发态偶极矩计算的结果都表明偶极矩的符号变化,其中CI 1-particle方法   夸大了这一变化,这表明了在激发态电子密度的迁移.   这一计算偶极矩的方法仍然利用了传统的方法,由于其根据波函数的平方来进行   计算,就倾向于产生这样的趋势.据此,我们强烈推荐采用波函数的解析导数来进   行偶极矩计算.设置为Density=CI或在CI-Silgles中的Density=Current   练习9.2 文件9_02 甲醛的激发态优化   确定甲醛第一激发态的结构,比较红外光谱.   下面是分子的Z-矩阵   C   O 1 RCO   X 1 1.0 2 ACO   H 1 RCH 3 ACH 2 90.   H 1 RCH 3 ACH 2 -90.   RCO=1.25   RCH=1.08   ACO=145.0   ACH=60.   这里使用了虚原子,使得对于分子的描述更加简单,这个分子也不再是共平面的   结构.计算采用校正因子0.8929,结果如下   计算值 校正值 实验值   平面外弯曲 495 442 683   CH2摇摆 978 873 898   CH2剪式 1426 1273 1290   CO伸缩 1647 1471 1173   对称CH伸缩 3200 2857 2847   反对称CH伸缩 3295 2942 2968   基本与实验值相符   练习9.3 文件9_03 丙烯醛的激发态优化   丙烯醛的激发态结构仍然是平面结构   练习9.4 文件 9_04a~b 苯的激发态优化   在CI-Singlets方法的研究过程中,有这样的文字,   进一步,我们可以得出结论,CI-Singlets方法的成功很大程度上依赖所选用的   基组,弥撒的激发态往往需要在分裂的价键轨道基组上增加弥散函数.   苯是这一结论的好的例子.   苯的激发态研究也显示了如下两点   * 理论预测必须与高精度的实验比较.允许跃迁可以与普通的单光子光谱结果   比较,但禁阻跃迁要有多光子实验证明.   * 当需要确定预测的和观测的状态时,必须应用激发态的对称性.一般的,   Gaussian提供激发态的对称性,但当不能提供时,需要检验跃迁的波函数.   本例中,没有使用弥散基组得到的激发态与实验结果相差较大,缺少了三个低能   量的激发态,而使用弥散基组的计算结果与实验吻合较好.   练习9.5 文件 9_05a~b 使用CASSCF方法研究激发态   练习使用CASSCF(Conplete Active Space Multiconfiguration SCF)方法.   研究体系为练习9.3研究的丙烯醛.这练习9.3中,第一激发态的能量相差较大,   本练习试图寻找减小误差的.   CASSCF计算在SCF计算基础上增加了完全电子相关计算.由于SCF计算有倾向基态   的趋势,CASSCF是其的改进.   关键词就是CASSCF,需要两个参数,活性空间的(Active Space)电子数量和轨道   数量.活性空间指CI计算中的轨道.活性空间的电子从HOMO轨道中取得,轨道由   活性空间电子的HOMO轨道和一些LUMO轨道组成.   例如一个单重态体系的CASSCF(4,6),表示活性空间由两个HOMO轨道(4个电子),   和4个LUMO轨道组成.   Guess=Alter用于指定组成活性空间的轨道.一个快速的Guess=Only和/或   Pop=Full计算用于生成具体的轨道信息,来确定参与活性空间的轨道.   采用CASSCF方法确定丙烯醛的第一激发态能量,需要如下步骤,   * 运行UHF/STO-3G Pop=NaturalOrbitals计算三重态丙烯醛的初始轨道和对   称性,选择将要组成活性空间的轨道.由于第一激发态电子由氧原子的孤对   电子跃迁到CO键的pi轨道,所需要的组成活性空间的轨道就要包括氧的孤对   电子以及羰基的pi轨道.   * 在6-31G(d)基组确定基态和激发态的能量.由于CASSCF波函数很难收敛,我们   需要做如下工作   * 运行CASSCF(6,5,UNO)/STO-3G Guess=(Read, Alter)计算三重态丙烯醛,   改进临时文件中的轨道.之所以采用三重态,是因为它比较容易收敛.   * 从STO-3G的波函数,6-31G(d)基组进行第二个计算.   * 采用6-31G(d)基组计算能量,然后,在前一个计算得到的波函数基础上,   设置多重度为1,Geom=Check,进行计算,这里需要给出NRoot=2来确定激   发态.   * 运行另一个CASSCF(6,5)/6-31G(d)计算得到基态能量,方法与得到激发态结果   的方法相同.   CASSCF方法不再是个黑箱子,它需要使用者的参与,而且必须很小心,很有耐心.   注意从默认的猜测进行的CASSCF计算几乎肯定是要失败的.选择合适的活性空间   轨道需要对所研究的体系有深刻的认识.   下面是计算的具体方法   首先是第一组计算,检验轨道,确定活性空间组成.   %Chk=acro_cas   #T UHF/STO-3G Test Pop=NaturalOrbotals   UHF on triplet acrolein at CAS(6,5) 6-31G(d) geometry   0 3   molecule specification   --Link1--   %Chk=acro_cas   %Chk=acro_cas   #T CAS(6,5,UNO) Guess=(Read,Only) Test Geom=Check   ...   下面是从第二次得到的收敛的波函数的对称性.   Orbital Symmetries:   Occupied (A') (A') (A') (A') (A') (A") (A') (A') (A')   (A') (A') (A') (A") (A') (A")   Virtual (A") (A') (A') (A') (A') (A') (A') (A')   活性空间将要由4个最高占据轨道和一个最低空轨道组成.我们需要四个A"和   一个A'轨道,这样,我们就必须讲第6号轨道与13或15号轨道置换.这样,我们就   必须检验这些轨道究竟代表什么,在分子的什么地方.这里,15号轨道在氧原子   的2px,2py轨道上有很大的系数,而其他几乎为零,所以需要保留这条轨道.   这样,就要置换6和13号轨道,输入文件如下   %Chk=acro_cas   #T CAS(6,5,UNO)/STO-3G Test Geom=Check Guess=(Read,Alter)   CAS 6,5 Using triplet UNO orbitals   0 3   ! Bring A" into the active space   6 13   这个计算的输出文件就是新的三重态的CASSCF表述.我们要用这个结果进行   进一步的计算,首先是更大基组的CASSCF计算,然后是单重态的CASSCF计算.   --Link1--   %Chk=acro_cas   %Chk=acro_cas   #T CAS(6,5)/6-31G(d) Test Geom=Check Guess=Read   CAS 6,5 in extended basis set   0 3   --Link1--   %Chk=acro_cas   #T CAS(6,5,NRoot=2)/6-31G(d) Test Geom=Check Guess=Read)   Singlet n-pi state (should have similar orbitals)   0 1   最后是基态的能量计算   --Link1--   %Chk=acro_cas   %NoSave   #T CAS(6,5)/6-31G(d) Test Geom=Check Guess=Read   Ground state starting from excited state orbitals   0 1   最终得到的能量差是4.035eV,与实验值3.71eV比较,已经很接近了.   练习9.6 文件 9_06a~d CASSCF方法研究丁二烯光化学   当光照射丁二烯分子时,分子由基态升至激发态,然后激发态又要通过一个非   辐射(不发射光子)衰减回到基态势能面.   这个非辐射衰减有两个方法解释.传统观点认为分子从激发态的势能面通过其   和基态势能面的最接近的点转到基态的势能面,这一点上,激发态的势能面仅   比基态势能面高一点.另外的观点认为,存在一个区域,两个势能面相交,跃迁   在相交的区域完成.   具体描述可以是,丁二烯(包括顺式和反式,我们这里讨论的是顺式)首先受激   跃迁至第一激发态,这一激发态是单重态的1Bu状态,表示的是pi体系的HOMO到   LUMO的跃迁.第二激发态是1Ag状态,表示的是HOMO到LUMO的双重跃迁,这两个   激发态的顺序并不清楚,但是1Bu在飞秒的时间内很快转换到1Ag.然后反应从   1Ag状态的势能面变化到接近,穿过基态的势能面,在接近或穿过基态势能面的   时候,分子降到基态.   本练习只研究这个体系的一小部分,在基态结构下研究三个状态的相对能量,   并确定两个势能面相交的点.   对于基态,完成如下步骤,   * 在RHF/3-21G水平检验轨道,来确定活性空间.这里将采用4电子空间,6条活性   轨道.轨道要包括pi轨道,并且具有A2,B1对称性(对应两个激发态的对称性)   * 使用选区的活性空间,在基态进行一系列态平均(State Average)4,4 CAS计算,   一般的,CAS计算优化轨道和所研究状态的波函数,而态平均计算提供符合所   需要的群的所有轨道的最佳描述.   在这个计算中,需要定义NRoot选项,在其所指定的状态和基态之间的所有状   态都将要被平均化.同时,各个状态的权重也需要提供.   这里我们关心最低的三个状态,所以,NRoot=3,设权重=0.333333.我们将在   3-21G, 4-31G, 6-31G(d),6-31+G(d)四个基组水平上运行平均化.   在最后一个计算的设置行设定#P选项来输出各个状态的表述.表述如下   VERTOR EIGENvalueS CORRESPONDING EIGENVECTOR   1 -154.87477374 0.95038483 -0.32993992E-01 .....   状态序号 能量 各个电子结构的系数   下面是4,4 CAS计算中的前两个激发态结构   PRIMARY BASIS FUNCTION=1 2 1 2   2 SYMMETRY TYPE=0   结构序号   1 2   1 3   活性空间中的电子位置   3 SYMMETRY TYPE=0   1 3   1 3   结构2是单重激发态,一个电子从活性空间的第二条轨道跃迁到第三条轨   道,结构3是双重激发态,活性空间第二轨道中的两个电子都跃迁到第三   轨道   * 计算完成后,表征两个激发态,比较能量   * 可以进行第二系列的4,6CAS计算,考察结果的变化   研究相交点,步骤如下   * 以已经得到的结构为起点,设置   UHF/STO-3G Guess=Mix Pop=NaturalOrbitals NoSymm   计算,检验轨道   * 建立好的CAS体系,运行两次4,4CAS态平均计算,在第一次运算中采用STO-3G   基组,设置UNO, NoFullDiog和CASSCF关键词;第二个计算中采用4-31G基组.   两个计算都有NoSymm关键词   * 运行CAS Opt=Conical, Nosymm, iop(1/8=5)确定相交点   * 在优化的相交点的结构上运行态平均计算,采用4-31G基组,   部分结果如下   方法 激发态1 激发态2 能量差   CAS(4,4)/6-31+G(d,p) 2 1Ag 1 1Bu 0.03702   CAS(4,6)/6-31+G(d,p) 2 1Ag 1 1Bu 0.00540   两个结果中,1Ag的能量都比1Bu低,随着方法的改善,能差减小.增加MP2计算   可以进一步缩小这个差距.   这个寻找相交点的计算还是成功的,这一结构上得到的基态和激发态的差距   仅有0.00014Hartree.