过渡态算例
QST2法寻找反应的过渡态
25 Nov
QST2法需要输入分子的起始结构和产物结构,两个分子的原子序号要求一致,否则不能运算。下面以H3CO 和H2COH 的制作说明之。
先在Gview 中制作H3CO ,保存为H3CO.gjf ;再制作H2COH ,保存为H2COH.gjf 。分别用Edit –>Atom List…打开两个制作的分子模型,Tag 那一列就是原子的序号,比较两个分子,发现4和5代表的原子在H2COH 中颠倒了,若直接用这样的结构,提交的时候,会给出“Invalid QST2 job! Atom list ordering must be the same for all input geometries! ”错误。这时,点击Tag 列中的4和5,分别改成5和4,然后把这一列按序号重新排一下(Rows –>Sort Selected),保存。
用ultraedit 打开两个分子,就可以看到所有原子的坐标(默认是Cartesian 坐标,在Gview 中保存的时候,去掉左下角的Write Cartesians的对勾,就可以保存在Z-matrix 坐标),拷贝到高斯输入文件中即可。制作好的输入文件如下:
============================================== #T UHF/6-31G(d) opt=QST2 Test
H3CO –>H2COH Reactions
0,2
C 0.00000000 0.00000000 0.00000000
H 0.00000000 0.00000000 1.07000001
H 1.00880579 0.00000000 -0.35666635
H -0.50440269 -0.87365135 -0.35666685
O -0.67410884 1.16759007 -0.47666623
H3CO –>H2COH Products
0,2
C 0.0000000 0.0000000 0.0000000
H 0.0000000 0.0000000 1.0700000
H 0.9277049 0.0000000 -0.5331638
H -1.3684675 0.9048960 -0.7874965
O -1.0873056 0.0000000 -0.6335117
============================================== 反应物
过渡态
产物
用Gaussian 研究过度态(简明教程,适合零基础,转载+错误更正+重新排版分节) Exploring Chemistry with Electronic Structure Method
===================================================
输入输出文件介绍
---------------------------------------------------
输入文件
#T RHF/6-31G(d) Test
My first Gaussian job :water single point energy
0 1
O -0.464 0.177 0.0
H -0.464 1.137 0.0
H 0.441 -0.143 0.0
第一行以#开头, 是运行的说明行,#T表示指打印重要的输出部分,#P表示打印更多的信息. 后面的RHF 表示限制性Hartree-Fock 方法, 这里要输入计算所选用的理论方法. 6-31G(d)是计算所采用的基组, 就是用什么样的函数组合来描述轨道.
Test 是指不记入Gaussian 工作档案, 对于单机版没有用处.
第三行是对于这个工作的描述, 写什么都行, 自己看懂就是了.
第二行是空行, 这个空行以及第四行的空行都是必须的.
第五行的两个数字, 分别是分子电荷和自旋多重度.
第六行以后是对于分子几何构性的描述. 这个例子采用的是迪卡尔坐标.
分子结构输入完成后要有一个空行.
输出文件
输出文件一般很长, 对于上面的输入文件, 其输出文件中, 首先是版权说明, 然后是作者,Pople 的名字在最后一个.
然后是Gaussian 读入输入文件的说明, 再将输入的分子坐标转换为标准内坐标, 这些东西都不用去管. 当然, 验证自己的分子构性对不对就要看这个地方.
关键的是有SCF Done 的一行, 后面的能量可是重要的, 单位是原子单位,Hartree, 1Hartree = 4.3597482E-18 Joules 或 =2625.500 kJ/mol=27.2116 eV; 再后面是布居分析, 有分子轨道情况, 各个轨道的本征值(能量), 各个原子的电荷, 偶极距. 然后是整个计算结果的一个总结, 各小节之间用\分开, 所要的东西基本在里面了.
然后是一句格言, 随机有Gaussian 程序从它的格言库里选出的(在l9999.exe 中, 想看的可以用文本格式打开这个文件, 自己去找, 学英语的好机会).[1]
然后是CPU 时间, 注意这不是真正的运行时间, 是CPU 运行的时间, 真正的时间要长一些. 如果几个工作一起做的话(Window 下好像不可能,Unix/Linix 下可以同时做多个工作), 实际计算时间就长很多了.
最后一句话,"Normal termination of Gaussian 94" 很关键, 如果没有这句话, 说明工作是失败的, 肯定在什么地方出错误了. 这是这里应该有出错信息.
===================================================
计算模型
---------------------------------------------------
计算化学的方法主要有分子力学理论(Molecular Mechanics) 和电子结构理论(ElectronicStructure Theory).
两者的共同点是:
计算分子的能量, 分子的性质可以根据能量按照一定的方法得到.
进行几何优化, 在起始结构的附近寻找具有最低的能量的结构. 几何优化是根据能量的一阶导数进行的.
计算分子内运动的频率. 计算依据是能量的二阶导数.
===================================================
分子力学理论
---------------------------------------------------
分子理论采用经典物理对分子进行处理, 可以在MM3,HyperChem,Quanta,Sybyl,Alchemy 等软件中看到. 根据所采用的力场的不同, 分子理论又分为很多种.
分子力学理论方法很便宜(做量化的经常用贵和便宜来描述计算, 实际上就是计算时间的长短, 因为对于要花钱上机的而言, 时间就是金钱; 对于自己有机器的, 要想算的快, 也要多在机器上花钱), 可以计算多达几千个原子的体系. 其缺点是
每一系列参数都是针对特定原子得出的. 没有对于原子各个状态的统一参数.
计算中忽略了电子, 只考虑键和原子, 自然就不能处理有很强电子效应的体系, 比如不能描述键的断裂.
电子结构理论
这一理论基于薛定鄂方程, 采用量子化学方法对分子进行处理. 主要有两类:
1. 半经验方法
包括AM1,MINDO/3,PM3,常见的软件包有MOPAC,AMPAC,HyperChem, 以及Gaussian. 半经验方法采用了一些实验得来的参数, 来帮助对薛定鄂方程的求解.
2. 从头算
从头算, 在解薛定鄂方程的过程中, 只采用了几个物理常数, 包括光速, 电子和核的质量, 普朗克常数, 在求解薛定鄂方程的过程中采用一系列的数学近似, 不同的近似也就导致了不同的方法. 最经典的是Hartree-Fock 方法, 缩写为HF.
从头算能够在很广泛的领域提供比较精确的信息, 当然计算量要比前面讲的方法大的多, 就是贵得多了.
密度泛函(Density Functional Methods)
密度泛函是最近几年兴起的第三类电子结构理论方法. 它采用泛函(以函数为变量的函数) 对薛定鄂方程进行求解, 由于密度泛函包涵了电子相关, 它的计算结果要比HF 方法好, 计算速度也快.
===================================================
一个理论模型, 必须适用于任何种类和大小体系, 它的应用限制只应该来自于计算条件的限制. 这里包括两点:
一个理论模型应该对于任何给定的核和电子有唯一的定义, 就是说, 对于解薛定鄂方程来讲, 分子结构本身就可以提供充分的信息.
一个理论模型是没有偏见的, 指不依靠于任何的化学结构和化学过程.
这样的理论可以被认为是化学理论模型(theoretical-model chemistry), 简称化学模型(modelchemistry).
===================================================
高斯计算中结构优化所用的原理
结构优化的方法一般是先计算出解析能量梯度, 即总能量对原子核坐标的微商. 再利用共厄梯度法计算步长, 其中的关键在于Hessian 矩阵的逆矩阵的构造, 我不清楚Gaussian 里具体用什么方法, 但是常用的是BFGS 方法, 再结合GDIIS 方法加速收敛.
===================================================
分子的几何构型(Molecular Geometry) 及优化
分子的平衡构型(molecular equilibrium geometry) 是分子电子能量和核间排斥能量最小时分子的核排列.
---------------------------------------------------
分子势能
一个含有N 个原子核的非线性分子的几何构型可以用3N-6 个独立的核坐标决定, 分子的电
子能量,U(q1,q2,…,q3N-6) 是这些坐标的函数.U = Ee +VNN
注意到3 个平移和3 个转动自由度(线性分子的转动自由度为2) 对U 是没有贡献的, 因此对一个双原子分子,U 的表达式中仅仅保护一个变量, 即两个核之间的距离,U(R).
对一个多原子分子,U 是每两个原子核之间距离的函数, 是分子势能面(potential energy sur face, PES) 的一部分. 对某一特定的分子核排列下U 的计算被成为单点(single-point)计算, 因为这一计算仅仅涉及到分子PES 上的一个点.
一个大分子可能在其PES 上有多个极小点, 对应于不同的平衡构象和鞍点. 分子构象(molecular conformation) 可以通过指定围绕单键的二面角的指得到. 在能量极小点处的分子构象称为构型(conformer).
---------------------------------------------------
几何构型优化
从初始几何构型出发寻找U 的极小值的过程称几何构型优化(geometry optimization) 或者能量极小化(energy minimization). 极小化的算法同时计算U 和U 梯度.
在一个局部最小点,U 的3N-6 个偏微分都是0.PES 上▽U=0 的点称为稳定点(stationary point) 或者判据点(critical point), 它可以是极小点, 极大点或者鞍点.
除了▽U 之外, 一些最小化方法使用到U 的二阶偏微分, 从而生成Hessian 矩阵, 又称为力常数(force constant) 矩阵, 因为d^2U/Qi^2 = fi 为力常数.
如果一个稳定点是电子能量面上的一个极小点, 其力常数矩阵的所有特征值都是正值. 然而, 若一个稳定点是过渡态(transition state, TS), 其中一个特征值是负值.
---------------------------------------------------
Newton-Rapson
Newton-Rapson 方法是一种非常有效的寻找多变量函数的局部极小点的算法, 它将函数用Taylor 展开到二次项, 包括函数的一次和二次微分, 并以此作为函数的近似.
Quasi-Newton-Rapson
计算自洽场(self consistent field, SCF) 能量的二阶微分是非常耗时的, 因此在优化时经常使用一种修正的方法, 即quasi-Newton(或quasi-Newton-Rapson) 方法. 这种方法在每一步优化中通过计算梯度对Hessian 值进行初始估算.
---------------------------------------------------
优化方法
为了优化几何构型, 要先对平衡构型做一个估算, 通常使用键长和键角的经验值. 此外, 我们还要选择好适当的方法和基组, 然后就可以在所估算的平衡构型附近进行极小点的搜索了. 软件对电子Schrodinger 方程进行求解并得到U 及其在初始构型处的梯度. 通过这些数值对H essian 矩阵进行估算, 并调整3N-6 个核坐标以得到在初始构型附近但能量更低的新的分子构型. 对新构型的U 和▽U 进行计算以继续改进分子坐标使分子构型更接近于极小点.SCF 计算对新的分子构型不断重复, 直到▽U 和0 之间的差足够小能满足对极小点的判据(critia).
一般而言, 需要进行3N-6 至两倍的SCF 循环次数以找到一个极小点.(Gaussian 默认的循环次数为两倍要优化的变量, 有时候对OPT=Tight/Vtight 需要增大循环次数).
---------------------------------------------------
中间体
局域最小点代表反应物, 产物或者一个中间体(Intermediate),该中间体在多步反应中是既是一个反应的产物, 同时又是另一个反应的反应物. 因为反应中间体通常很快被反应掉或者其寿命很短而不能被光谱仪器检测到, 中间体结构和能量的计算是计算化学中一个重要的应用. 反应表面
为了得到一个完整的反应表面U, 假设我们需要在解电子Schrodinger 方程时在反应面上对3N-6 个变量各取约10 个点, 这样一共要进行103N-6 个计算. 实际上并非如此, 我们利用反应面别的重要性质来得到局部最小点和过渡态. 对于利用U 表面特征的分析算法可以参照S zabo & Ostlund, Modern Quantum Chemistry, pp. 437 - 458.
---------------------------------------------------
TS 计算算法
为寻找过渡态, 在Gaussian 程序中有几种方法可供选择.
***************************************************
TS
在Gaussian 中的OPT=(TS,CalcFC)算法, 需要输入一个估计的过渡态结构, 该过渡态
是反应物和产物之间的一个中间体.CalcFC 指定对初始结构进行Hessian 数值计算(用Ne wton-Rapson 方法), 而不是默认地对其二阶微分进行估算(用quasi-Newton-Rapson 方法). 这一步骤增加了计算时间, 但增大了TS 寻找地成功率.
***************************************************
QST2
另一个Gaussian TS 搜索算法, 由OPT=QST2 指定, 让程序自己通过反应物和产物的结构计算一个TS 的初始结构. 然后Gaussian 用这一初始TS 作为TS 搜索的起始点. ***************************************************
QST3
第三种Gaussian TS 搜索方法, 用OPT=QST3 指定, 在反应物, 产物和用户给定的初始TS 基础上进行搜索. 判断一个鞍点和TS 构型, 我们需要进行频率计算. 如果我们想知道一个鞍点是否是TS 构型, 并要研究反应机理, 需要从鞍点出发, 沿U 表面的反应路径分别得到反应物和
产物. 对一个TS, 其虚频所对应的简振模式下的原子位移方向分别指象反应物和产物. 当计算出TS 构型和振动频率以及能垒高度后, 就可以通过过渡态理论(Transition State Theory) 对反应速率进行计算了.
***************************************************
一阶鞍点
如果一个稳定点, 即电子能量面上的一个极小点, 是一个过渡态TS, 其Hessian 矩阵有且仅有一个负特征值. 这个一阶鞍点条件意味着TS 的能量在一个方向(负特征值) 取得最大值, 这一方向即反应坐标方向; 而在其它方向上能量都是最小值(正特征值).
***************************************************
TS 振动模式
一个TS 有3N-7(线性分子为3N-6) 个标准振动, 比正常的分子少一个. 负的力常数的物理意义在于, 与Hooke 定律所定义的力不同, 与TS 虚频v #相关的力沿反应坐标, 指向相反的, 朝着产物的方向. 频率v #, 沿着反应坐标, 是一个虚数.
关于TS 搜索和优化
一般来说直接用OPT=QST2/3或OPT=TS优化过渡态是比较困难的对于到底什么反应应该用QST, 什么反应应该用TS 我也没有把握, 一般做过渡态优化可以用这样的方法.
QST2必须输入反应物和产物坐标,Gaussian 猜测一个初始的TS 坐标, 即将二者取平均 我觉得对你这种反应不太适合比如你说的这个硅烷脱氢反应, 最好先把反应机理弄明白. ***************************************************
如果假设其机理是通过下面的TS 完成的:反应 SiH4-->TS-->SiH2+H2,那你就要先选定反
应坐标(Reaction Coordinate,RC) 反应坐标是在反应中发生变化的坐标, 一般是用内坐标例如这个反应中反应物键角A(H1-Si-H2)从109.5度变到TS 处大约60度, 产物则更小键长R(H1-H2) 也发生变化, 因此这两个坐标都可以设置为RC 选定RC 后, 可以通过做一potential e nergy profile 来找到TS 的初始结构 即固定该RC, 使之慢慢从反应物变化到产物并优化 T S 在所得到的势能曲线最高点附近. 这时候在势能曲线最高点两边分别取点作为反应物和产物结构做QST2或QST3, 或者直接以所得到的TS 初始结构做OPT=TS成功率就大了. 假设一个合理的反应机理对TS 搜索非常重要, 有的反应可能不止有一个TS, 有可能是通过A -->TS1-->A1-->TS2-->...-->P这样的路线完成的, 其中A1可以叫中间体这时候每个TS 都要分别寻找, 并通过反应速率的计算判断哪一步是反应决速步 TS 搜索和优化是一个非常复杂的问题, 不是某个关键词就能解决的.
OPT=TS最好和CalcFC 联用, 前着是用Quasi-Newton-Rhapson 方法用U 的一阶微分来估算力常数的, 可能往往不收敛,CalcFC 是用Newton-Rhapson 方法直接对二阶微分进行计算来求力常数, 加大了计算量, 也加大了成功率
***************************************************
输入
QST2需要反应物和产物结构,QST3需要反应物, 产物和TS 初始结构, TS 只需要TS 初始结构(通过内禀反应坐标搜索), 如:
# RHF/6-31G(d) OPT=QST2 Freq Test
Untitled
0 1
the configuration of SiH4
Untitled
0 1
the configuration of SiH2+H2