计算化学-10.计算类型及高斯实现方法.ppt

袁淑萍袁淑萍 国家重点实验室培育基地国家重点实验室培育基地408408室室 E-mail: E-mail: 电话：电话：8595069085950690，1358920823513589208235 计算化学计算化学 计算作业类型 势能面扫描 单点能 几何优化 反应过渡态 寻找 反应路径 振动频率分 析 几何优化的目的 寻找势能面上的极小 点, 确定分子可能的稳定结构 几何优化得到的仅仅是势能面上的局部极小点! 能量 OPT的更多选项 Maxcycle=n: 设定优化的最大步数 NoFreeze: 活 化所有冻结的坐标 (Geom=Check) Expert: 在力常 数和步长上放松限制 Tight, VeryTight: 提高收敛标准(力和 步长) Loose: 在工作 初期使用 MaxStep=m: 最大步长 = 0.01 * m； IOP(1/8=m): 改变步长 0.01 * m Restart: 优化 作业的重启 GEOM的选项 Checkpoint: 从 Chekpoint文件中读入分子数据( 通常还用到Guess=Read). Modify: 读入并 修改Checkpoint文件中的分子数 据. NoCrowd: 允许 原子间距小于0.5A. NoKeep: 丢弃 关于冻结的信息. Step=n: 在前一 次优化失败后, 把它的第n步作 为起点进行再次优化（ Geom=(Check, Step=n) ） Handling Difficult Optimization Cases 复杂体系的优化 q对于大体系，先用简单方法作一粗略优化、后用高级方 法和高级基组优化，有时可能更快地达到目的。 q后续作业可以从前步作业的结果中取得几何构型和波函 数作为继续优化的起点，此时，后续作业再不需要输入 分子几何（但仍需要分子电荷和自旋多重度)。 q也可先用分子力学优化，作为SP的输入构型或进一步作 量子化学优化的初始构型。分子力学优化相当快，机时 远少于量子化学优化。 几何优化的策略 初学者有时急于求成，希望一步成功，一开始就 采用高级方法和高级基组而忘了欲速则不达的道 理。 q 优化中途变量出错。最常见的现象是键角超越了 最大限度180o： q可以减小步长OPT(stepsize=N)或在进程行增加 Opt=calcfc q先用较低水平的方法和小基组作频率计算，再在后续的 优化作业中使用Opt=Readfc。但这种计算往往是费时的 。 q 初始构型或优化过程中的某种构型在势能面上处 于较为平坦的区域，容易造成振荡, 甚至永远不 收敛。 q解决的办法是加一个哑原子, 重新定义键角。 q步数超出 q检查非常容 易改变的坐标和/或强烈耦合的坐 标 q增加几何优 化的最大步数 OPT (Restart, Maxcyc=N) q在优化过程中点 群改变了 q检查结构和 /或不使用对称性(NoSymm) 批处理任务 任务比较多，互相之间没有关联，不能用多 步处理的方法写在同一个输入文件里。 编辑批处理文件，程序按次序调用并逐步计 算 计算作业类型 势能面扫描 单点能 几何优化 反应过渡态 搜索 反应路径 振动频率分 析 搜索过渡态的目的 寻找势能面上的一级鞍点，确定分子可能的稳定 结构之间的变换可能经过的态 一级鞍点满足的条件: 过渡态优化算法 过渡态优化 的算法: TS，QST2， QST3 寻找过渡态 的实用建议 关键词: Opt=QST2 Optimization of Transition States OPT = TS -输入初始 猜测的过渡态结构 -确认初始 Hessian矩阵有一个本征向 量有负的本征值 -可能的话 计算力常数矩阵(CALCFC 或READFC) 寻找过渡态 用OPT=QST2 或 OPT=QST3 在冗余内坐标下 寻找 QST2: 输入反应 物附近和产物附近的两个结构(通 过对冗余坐标的线性插值估计过 渡态) QST3: 输入反应 物、产物和估计的过渡态 Opt=QST2 的输入 #B3LYP/6-31G(d,p) OPT=QST2 H3CO-Title 1 0 2 C1 0.0 0.0 0.0 02 0.0 0.0 1.3 H3 0.0 0.9 -.3 H4 0.8 -.2 -.6 H5 -.8 -.2 -.6 CH2OH - Title 2 0 2 C1 0.0 0.0 0.0 02 0.0 0.0 1.4 H3 0.0 0.92 1.7 H4 0.7 -.1 -.7 H5 -.7 -.1 -.7 每个结构找到原子 顺序应该一致 输入结构不需要是 优化的稳定结构 QST3 增加了第3个 标题和过渡态结构 的估计 为过渡态估算Hessian矩阵 初始的Hessian 矩阵必须有一个负的本征值, 及 其合适的本征向量 从低等级 计算得到近似的Hessian矩阵 从低等级 的完整Hessian矩阵来计算 (READFC : 从频率计算中读入) 在同等级 计算完整的Hessian矩阵 (CALCFC) 在优化的 每一步都重新计算完整的Hessian 矩阵 (CALCALL)，非常耗时 过渡结构算例 # T U H F / 6 - 3 1 G ( d ) O p t = Q S T 2 H 3 C O - - H 2 C O H R e a c t a n t s 0 , 2 C O 1 1 . 4 8 H 1 R 2 A H 1 1 . 0 8 2 1 1 0 . 3 1 2 0 . H 1 1 . 0 8 2 1 1 0 . 3 - 1 2 0 . R = 1 . 0 8 A = 1 1 0 . H 3 C O - - H 2 C O H p r o d u c t s 0 , 2 C O 1 1 . 4 8 H 1 R 2 A H 1 1 . 0 8 2 1 1 0 . 3 1 2 0 . H 1 1 . 0 8 2 1 1 0 . 3 - 1 2 0 . R = 1 . 9 A = 3 0 . 过渡结构的确认 计算完整的 Hessian矩阵，确保有且仅有一个负 本征值 通过振动分析计算完 成，可以使用连续作业 # B3LYP /6 -31g opt=(ts,calcfc) freq 确认本征矢量的 方向正确连接到反应物和产物 通过计算反应途径进 行确认（IRC） 寻找过渡态失败怎么办 1)有多 余的负本征值 沿着不 是过渡向量的那些负本征值 向量的方向寻找 2) 没有负本征值 在反应坐标上扩大扫描, 寻找最高能量对应的 坐标点，沿最小本征值方向优化 需要注意的问题 一般过渡态 都涉及到键的断裂与生成 ，因此即使是自旋多重度 为1的体系，也应该用非限 制性的理论方法来寻找。 设想一个键 的均裂反应，原来成对的 两个电子在反应过程中是 要分别属于两个自由基的 ，不会占据同一个空间轨 道了。 计算作业类型 势能 面扫描 单点 能 几何 优化 反应 过渡态搜索 振动 频率分析 反应 路径 振动频率分析计算可以作为一个单 独的作业，但更常见的是附在分子构型优化之后 成为多步作业，以验证优化结果是否稳定构型或 过渡态。 频率计算与构型优化用同样的计算 方法和基组进行，关键词是Freq。下面是频率计 算作为后继作业的输入实例： % chk=H2O.chk # HF/6-31G(d) Freq Geom=Checkpoint Test H2O Frequencies 0 1 频率计算的目的 预测分子的红外和拉曼光谱(频率和强度) 在HF水平计算时，由于忽略相关能，计算的频率高于实 验值。通常在计算前或计算后利用校正因子来校正。 影响频率计算准确性的因素还有：确定力常数矩阵所需 的微分技术的不准确性和平衡几何构型的不确定性等。 几何构型优化时为难进行的优化计算力常 数矩阵 力常数是能量对坐标的导数。从数学的角度，导数显示 着函数的变化趋势（可以根据导数值确定函数在某个范 围内是变大还是变小）。 优化的目的是寻找对应能量最低的构型。所以对难进行 的优化可以通过提供力常数值来指引优化的方向，使优 化能够顺利完成。 确定势能面上驻点性质（极小点？鞍点 ？） 虚频率的有无及其数目可以帮助我们判断优化结果是鞍 点或能量极小构型，分子稳定构型不应有虚频率，而过 渡态则应该且只有一个虚频率。 计算零点能和热化学数据，对体系总能对体系总能 量、量、 焓、熵进行校正焓、熵进行校正 极化度与超极化度。 Hessian矩阵本征值与振动模式的关系 1，Hessian矩阵的 本征值对应于简正振动的力常数 2，简正振动的频 率等于力常数的平方根 3，简正振动的振 动模式是相应本征值的本征矢 振动强度 振动强度用于 光谱指认 IR光谱的振动 带的强度由偶极矩对简正模 式的微商确定 Raman光谱的 振动带强度由极化率对简正 模式的微商的平方确定 振动分析算例 %ch k=freqhessian.chk #p UHF/6-31G(d) freq=(noraman) H3 CO H2COH TRANSITION STATE 0 2 C 0.03104 0.63055 0.0 O 0.03104 -0.73696 0.0 H - 0.99138 -0.13555 0.0 H 0.27839 1.12398 0.92614 H 0.27839 1.12398 - 0.92614 注意:频率计算只能在优化 好的结构上进行，只能使用 与优化方法同样的方法来做 频率计算。 振动频率的校正因子 计算得到的简正频率比实验值一般高约10% 这是由于谐振子近似和理论的近似而产生的 热化学数据 所有频率计算都包含有体系的热化学分析。默认是在 298.15K一个大气压条件下计算得到。对每一种元素都使 用含量最高的同位素。下面是热化学的输出部分: - Thermochemistry - Temperature 298.150 Kelvin. Pressure 1.00000 Atm(计算 的具体条件：温度，气压) Atom 1 has atomic number 35 and mass 78.91834 (计算所 用元素的确切质量) Atom 2 has atomic number 16 and mass 31.97207 。 可以通过参数ReadIsotopes根据自己的具体需要设定计算条 件 Freq(ReadIsotopes)选项： Specify alternate temperature, pressure, and/or isotopes This information appears in a separate input section having the format: temp pressure scale Must be real numbers isotope mass for atom 1 isotope mass for atom 2 . isotope mass for atom n where temp, pressure, and scale are the desired temperature, pressure, and an optional scale factor for frequency data when used for thermochemical analysis (the default is unscaled). %Chk=freq # HF/6-31G(d,p) Freq Test Frequencies at STP molecule specification -Link1 %Chk=freq %NoSave # HF/6-31G(d,p) Freq(ReadIso, ReadFC) Geom=Check Test Repeat at 300 K 0,1 300.0 1.0 16 2 3 . The isotope masses for the various atoms in the molecule, arranged in the same order as they appeared in the molecule specification section. If integers are used to specify the atomic masses, the program will automatically use the corresponding actual exact isotopic mass (e.g., 18 specifies O18, and Gaussian uses the value 17.99916). 后面四行的四个能量分别为E0，E，H，G。计算关系为 E0=E(elec)+ZPE (表示零K时体系的能量等于Eelec（就是优化计算中所 得的SCF能量）和零点能之和) E=E0+E(vib)+E(rot)+E(trans) (在某一特定温度和压力下体系的能 量等于零K时体系的能量加上内能值) H=E+RT (在某一特定温度和压力下体系的焓等于在这一条件下的能量 和RT之和) G=H-TS (在某一特定温度和压力下体系的自由能等于在这一条件下的 焓和TS的差。S表示熵 ) 零点能是对分子的电子能量的矫正，表明了在0K温度下的分子的振动能 量对电子能量的影响。需要根据所用的算法采用不同的因子进行矫正。 在较高的温度下预测体系能量时，要给总能加上内 能项。 内能包括特定温度下的平动能，转动能和振动能， 这些因素都受到具体条件的影响。 E (Thermal) Contributions to the thermal en