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