课题1笔记和操作记录

发表于 2022-10-25 更新于 2023-07-04 分类于课题阅读次数：

多看看

任务要求

计算某个简单反应的中间体和过渡态结构。其中猜测的中间体结构已经给出，需要对相关结构进行优化，并且需要找到相应结构的中间体和过渡态。

任务特征、过程与原理

任务特征与反应过程

特征

没有涉及金属有机化合物，仅仅只是对中间体进行结构优化。

反应过程

该实例中有三个中间体（A2-A4），三个过渡态（TS-1~TS-3）。

其中反应物为含碳自由基的结构和芳香醛基烯【CHO-1】结构。

反应从初始到最后分别产生了吸引（A2）——成环（A3）——氧自由基进攻后开环（A4）三个部分。现在一方面需要优化并得出三个中间体的结构。

反应原理（烯烃自由基加成）

主要涉及的原理包括：

A1+CHO-1的C自由基加成（过渡态TS1）
A2自身结构中的C自由基进攻羰基成环（过渡态TS2）
A3自身结构中的氧自由基进攻C-C键开环（过渡态TS3）

其中自由基加成的原理：

引发：催化剂的自由基传递到反应物X-Y，并使Y获得自由基；
传递：烯烃和Y自由基作用加成后，α-C获得自由基，再与另一分子X-Y作用使X连接到α-C；
终止：α-C自由基既有可能自身偶联也有可能与烯烃反应生成其他形式的新的自由基物种。

自由基的加成取向由进攻两个双键碳原子的速度相对大小来决定。实践证明：影响自由基加成进攻速度的因素主要是立体效应对进攻自由基的阻碍作用和共轭效应对生成自由基的稳定作用。

对于链端末取代的α-烯烃CH2=CH-A或CH2=CAB， A、B的空间位阻效应以及A、B的共轭效应驿自由基的稳定作用是一致的，均使自由基进攻α-C的速度远大于β-C的进攻速度，结果得到自由基加在于α-C上为主的产物。（如：(CH3)2C=CH2 + HBr → (CH3)2CHCH2Br（主产物））
对于多取代烯烃，如果A、B对两个双键碳原子的立体阻碍作用和共轭稳定作用一致，则二者共同指导自由基进攻的位置。

总之：不论取代基和自由基的性质如何，自由基对α-烯烃的加成方向总是优先选择未取代端。

输入文件相关

输入文件的参数已经事先设计好：

对于CHO-1的：

1
2
3

%mem=64GB
%nproc=26
# b3lyp/6-31g(d) 5d opt freq empiricaldispersion=gd3bj

对于过渡态结构【Radd-TS1(TS1)】：

1
2
3

%mem=64GB
%nprocshared=16
# opt=(TS,calcfc,noeigen) freq b3lyp/6-31g(d) 5d empiricaldispersion=gd3bj

任务分析与概念说明

输入文件的最前面两行并不重要，主要看相关参数的设置。

DFT泛函和基组和基组关键词

b3lyp/6-31g(d) 5d

B3LYP是DFT泛函，在大多情况下具有较高效率对优化效果。

参考：谈谈量子化学研究中什么时候用B3LYP泛函优化几何结构是适当的

值得注意的有：

用B3LYP优化主族元素构成的体系（不含大共轭特征、不显著涉及弱相互作用）一般能得到合理的结果。
如果你对含有过渡金属体系的计算关注的重点只是配体部分，比如单纯是配体部分发生的反应等，那么用B3LYP优化完全没关系。
如果你算的体系牵扯到d族过渡金属之间直接成键，此时静态相关很强，一定要用纯泛函优化。
B3LYP优化有机反应的过渡态一般来说是能给出合理结果的，即便到现在依然使用非常广泛。

部分内容之前已经摘抄到这篇文章中：量子化学计算中DFT泛函和基组的选择

**6-31g(d)**则是对应的基组。

参考高斯官网的说明：Basis Sets

Single first polarization functions can also be requested using the usual * or ** notation.
Note that (d,p) and ** are synonymous (=same meaning) —6-31G** is equivalent to 6-31G(d,p), for example—and that the 3-21G* basis set has polarization functions on second row atoms only. The + and ++ diffuse functions are available with some basis sets, as are multiple polarization functions.
The keyword syntax is best illustrated by example: 6-31+G(3df,2p) designates the 6-31G basis set supplemented by diffuse functions, 3 sets of d functions and one set of f functions on heavy atoms, and supplemented by 2 sets of p functions on hydrogens.
(which means: 6-31G†=6-31G*=6-31G(d); 6-31G‡=6-31G**=6-31G(d,p), while ‘+ and ++’ is describing diffuse functions.)
这篇帖子的回复：
6-31g**是6-31g(d,p)的等价写法，d代表给除氢、氦外的原子加一层d极化函数，p代表给氢、氦加一层p极化函数，d极化函数有两种形式：笛卡尔型高斯函数和球谐型高斯函数。不同的程序会使用不同的形式。

根据之前的文章，6-31G*应该是相关计算体系中的最低要求。因此这里用于优化简单结构也是合适的。

5d是基组修饰关键词。

根据谈谈5d、6d型d壳层基函数与它们在Gaussian中的标识和这篇帖子的回复：

在高斯等量化程序中，对于非半经验方法用的都是GTF(Gaussian type function)型实函数用来描述原子轨道，以解决双电子多中心积分的麻烦，所以与复数型轨道完全无关。
最常用的高斯型函数是笛卡尔型高斯函数，通式为N_xayb_zcexp(-qr2)，N为归一化常数，q为高斯函数的指数，r为以此基函数所在原子为中心的坐标向量。
对于描述d原子轨道，使用d型笛卡尔型GTF作为基函数，包含六种具体形式:xx,yy,zz,xy,xz,yz。例如xy型GTF，其通式中的a=b=1，c=0。然而这6种里面除了xy、yz、xz以外的其它三种笛卡尔型GTF并没有与真实的实形式的原子轨道有直接的对应关系，波函数角度部分的行为明显不一致。
还有一种球谐型高斯函数，也叫原子轨道型高斯函数，通式为N_rnexp(-q*r2)_Y(θ,φ)，其中Y(θ,φ)就是类氢原子轨道的角度部分函数，与STO的区别是指数项的r变为了r^2。球谐型高斯函数与实型原子轨道角度部分行为一致，通过调整收缩系数调整径向行为，就可以一一对应地近似描述实型原子轨道。
高斯一般默认使用6d关键字（在使用gen等关键字时往往默认为5d）。此时说明使用xx、yy、zz、yz、yx、zx这6种笛卡尔型GTF轨道，在输出的d轨道符号中也对应地以xx、yy、zz、yz、yx、zx标识。在d轨道基函数不分裂情况下，描述d轨道的可独立变分的基函数就相对于5d关键字时的5个增加到了6个，可看到输出文件中basis functions的数目对应地增加了。
6d型轨道好处是方便编程，能够直接计算各种积分。缺点很明显，由于缺乏与真实原子轨道的对应关系，结果不像用5d的结果那样方便分析。虽然这6个笛卡尔型轨道不直接对应于真实原子轨道，但在计算过程中经过变分，其结果同样会展现出真实原子轨道的行为。这是因为以它们为基函数变分的结果，等价于对它们线性变换后的基函数变分的结果。前面已提到，6个笛卡尔型d型GTF线性变换后的6个函数其中的5个就是对应真实轨道的5d型轨道，只不过多增加了一个r^2轨道。换句话说，使用6d关键字，就是在5d型的基函数基础上多添加了一个内部含有节面的s型GTF轨道，由于这个轨道与其它s轨道有不小重叠，会造成一定线性相关问题，所以并不会比5d的结果有多少改进。
对于6-31g系列的基组，Gaussian中默认是笛卡尔型高斯函数（6D），但是如果用gen关键词，然后再写6-31g**则是球谐型高斯函数（5D），如果在Gaussian中不想用笛卡尔型高斯函数（6D），可以写上独立的关键词5d。对于6-311g系列的基组，就默认是球谐型高斯函数（5D）。

总结就是：对于6-31g系列的基组，Gaussian中默认是笛卡尔型高斯函数（6D），但是不想用笛卡尔型高斯函数，可以写上独立的关键词5d，这种情况下的对轨道能量对分析会更简单一些。

基准振动频率

opt freq

首先是opt和freq关键词：

freq和opt关键词一般连用。 用于计算基准振动频率，可计算结构中所有的振动模式与参数。

opt关键词单独使用，一般指通过角度和距离参数优化，找到能量极小点【即平衡结构/最适结构】；

freq用于计算力常数及其振动频率，另外还计算强度。

参考：[Numerical frequency calculations][https://sites.google.com/site/orcainputlibrary/vibrational-frequencies-thermochemistry](Google sites)

Always perform frequency calculations on optimized structures. Theory level must be identical for the structure and frequencies or you will get imaginary frequencies or meaningless results. A combined Opt+Freq job is the easiest way of making sure of this (even though you have already optimized the geometry).

在Radd-TS1中，opt后还添加了括号表明优化细节。

opt=(TS,calcfc,noeigen)

使用上述关键词一般用于计算过渡态变化。

参考：简谈Gaussian里找过渡态的关键词opt=TS和QST2、3：

TS代表向过渡态优化的关键词，其算法和优化极小点结构非常类似，只优化方向是势能面的一阶鞍点；
CalcFC代表精确计算初始结构的Hessian矩阵（能量相对于核坐标的二阶导数），是给过渡态计算提供Hessian/黑塞矩阵的关键词。（注：黑塞矩阵是一个多元函数的二阶偏导数构成的方阵，描述了函数的局部曲率。黑塞矩阵常用于牛顿法解决优化问题，利用黑塞矩阵可判定多元函数的极值问题。）
NoEigen代表优化过程中不对每一步做Hessian矩阵本征值数目的检测。换而言之就是为了提高计算速度，回避了优化过程中，Hessian 矩阵中可能会出现多个负特征值的情况。

后记：本次计算的目的是算中间体而不是算过渡态，所以上述写法其实是有问题的。

正确的写法应该是：opt freq

色散关键词

empiricaldispersion=gd3bj

empiricaldispersion代表启用经验色散矫正的关键词，在不带色散计算效果的基组（即不含D字母结尾）中作为补充进行使用（如B3LYP、M06等）。不同的关键词代表不同的经验色散矫正模型。

参考：GD3BJ

Add the D3 version of Grimme’s dispersion with Becke-Johnson damping [Grimme11]. The table below gives the list of functionals in Gaussian 16 for which GD3BJ parameters are defined. The functionals highlighted in bold include this dispersion model by default when the indicated keyword is specified (e.g., B2PLYPD3). For the rest of the functionals, dispersion is requested with EmpiricalDispersion=GD3BJ.

反之，如果基组中含D结尾，如B2PLYPD3，则无需empiricaldispersion=gd3bj。

另外，还有PFD、GD2、GD3等不同等色散模型。

理想化溶剂模型和高斯中的使用

高斯软件中，将复杂的溶剂模型做了极端化的隐式模型简化，这意味着将液相的溶剂体系简化为不同介电常数的“溶剂场”。

【类似中学物理的电场概念，可以理解为不同极性的溶剂对于溶质分子产生不同的“电场力”。】这一体系的简化尽管大幅简化了计算量，但是也无法解释一些溶剂对反应选择性产生影响的现象。这种情况下的溶剂模型被称为理想化溶剂模型。

隐式溶剂化

隐式溶剂化（有时称为连续溶剂化）是一种表示方法。溶剂作为一种连续的介质，而不是单个的“显式”溶剂分子。该方法常用于估算自由能的溶质-溶剂结构和化学过程中的相互作用。

参考：隐式溶剂化 - Implicit solvation

隐式溶剂模型就是不具体描述溶质附近的溶剂分子的具体结构和分布，而是把溶剂环境简单地当成可极化的连续介质来考虑的。这种考虑溶剂效应的好处是可以表现溶剂的平均效应而不需要像显式溶剂模型那样需要考虑各种可能的溶剂层分子的排布方式，而且不至于令计算耗时增加很高，因此被广泛用于量子化学和分子模拟领域。

参考：谈谈隐式溶剂模型下溶解自由能和体系自由能的计算

SMD（Solvation Model Based on Density，常用溶剂化模型的一种）

SMD方法是一种基于SCRF参数化的溶剂化模型，由Truhlar等人开发。SMD设计用于溶剂化自由能的预测，其所用的参数不同于SCRF模型的默认参数。对于一个给定的溶质，SMD溶剂化自由能计算的是在298.15K将同样浓度的理想气体转移到理想溶液的自由能。也就是说计算得到的溶剂化自由能是指1M气体->1M溶液的溶剂化自由能。而1M浓度的气体等于24.5 atm下的1 mole气体，而非标准状态的理想气体（在1个标准大气压下的气体）。

参考：Gaussian教程 | 溶剂化自由能的计算

SMD的使用

一般内置溶剂：scrf里写solvent=xxx即可，xxx是溶剂名字。比如如果想用SMD溶剂模型表现乙醇溶剂，就写scrf(SMD,solvent=ethanol)。
其他自定义溶剂和混合溶剂的情况我们组的计算经验好像没有考虑过，这里略过。

计算特定溶剂的单点能时，按照以下模板套用：

# b3lyp/6-31g(d) scrf(SMD,solvent=ethanol)

不标注opt一类的内容时，即计算单点能【SP】；
然后放入一个精度较高的基组&方法；
最后放scrf和smd溶剂模型即可。

计算溶剂模型中需要记录的数据

气相优化的记录数据（任何情况都需要记录）

气相（理想气体）中优化的数据是所有DFT都需要记录的数据，主要包括：

（Results-Summary）

上图中对应的内容如下：

英语名称	中文名	单位
Electronic Energy（EE）	电子自由能	Hartee
Thermal Correction to Enthalpy（Hcorrection）	焓变热力学修正	Hartee
Thermal Correction to Free Energy（Gcorrection）	自由能热力学修正	Hartee
EE + Hcorrection	-	Hartee
EE + Gcorrection	-	Hartee

注意：高斯默认单位Hartee需要换算为电子伏特/eV才能使用，换算系数为：1 hartee to eV = 27.2114 eV。

液相环境的记录和计算数据

一般在液相环境中，较少直接进行DFT运算，这是因为溶剂与一些金属有机分子的相互作用可能会导致结构难以收敛。

我们主要在此过程获得溶剂环境下的溶质自由能（溶质在气相下的自由能加上溶解自由能）：

更准确的说法是：

溶质在溶剂环境下标况（298.15K、1M）时的自由能 = 溶质在1atm气相下的自由能 + 直接通过隐式溶剂模型计算的溶解自由能（始末都是1M） + 1.89kcal/mol（标况下1atm→1M自由能变）

以在我们组为例，首先通过上面提到的内容获得两个热力学校正值，即自由能矫正Gcorr和焓变矫正Hcorr。

然后通过SMD模型下单点计算获得指定溶剂在SMD溶剂模型下的电子能量。

最终可得到：G solvent = E solvent + G correction ; H solvent = E solvent + H correction

注意点：

默认设定下，通过freq关键词或者热力学组合方法得到的气相下的自由能就是标况（298.15K、1atm）状态下的。
对绝大多数中性且局部不显离子性的分子，气相和溶剂下结构差异甚微，此时结构优化及之后的任务用气相下优化的结构就够了。但如果溶剂效应对当前溶质结构影响可能较大，或者想追求稳妥，结构优化、振动分析获得热力学校正量都应当在隐式溶剂模型下进行，并且之后算溶解自由能和高精度单点能的时候也用这个结构。

溶剂部分的写法

参考自：http://bbs.keinsci.com/thread-10624-1-1.html

Gaussian中用scrf=solvent=x定义溶剂的时候，以下每一行的写法是等价的，因此可以挑简单的、便于记忆的形式写

[acc status=”close” title=”溶剂写法”]

Water H2O
Acetonitrile CH3CN
Methanol CH3OH
Ethanol CH3CH2OH
DiethylEther Ethylether Ether
Dichloromethane MethyleneChloride CH2Cl2
DiChloroEthane 1,2-DiChloroEthane CH2ClCH2Cl
CarbonTetraChloride CCl4
Benzene C6H6
Toluene C6H5CH3
ChloroBenzene C6H5Cl
NitroMethane CH3NO2
Heptane n-Heptane
CycloHexane C6H12
Aniline C6H5NH2
Acetone CH3COCH3
TetraHydroFuran THF
DiMethylSulfoxide DMSO CH3SOCH3
n-Octanol 1-Octanol
[/acc]

区分过渡态、中间体、内坐标优化的关键词

判断过渡态结构的正确性

打开过渡态结构的频率分析窗口（右键Results-Vibrations）

第一个行存在一个虚频（即频率值小于0），然后点击下方的start Animation观察一下振动方式：

ts-1

可以看到振动方向和成键（C-C）方向是一致的，因此这个结构暂时可以认为符合要求。

区分过渡态和中间体的方法

虚频判断。过渡态：有虚频（小于0）；中间体：没有虚频，全正；
输入文件。过渡态：opt=(TS,calcfc,noeigen)；中间体：opt
断键位置的键长。过渡态：超过正常键长；中间体：在正常键长范围内。

常见的三种输入文件opt部分设置

过渡态：*# opt=(TS,calcfc,noeigen) freq*
中间体：*# opt freq*
固定键长：*# opt=z-matrix freq 文件尾也需要加内容*

限制性优化：内坐标优化介绍

上述固定键长属于限制性优化，详细可以参考：在Gaussian中做限制性优化的方法

量化计算时一般做的几何优化叫做全优化(full optimization)，与之相对的叫做限制性优化，即优化时冻结某些变量来实现特殊目的。对于单分子体系，在Gview自动建立的内坐标系下优化一般比在笛卡尔坐标系下优化收敛会快。这不难理解。比如对应化学键长度的坐标之间，尤其是相隔较远这样的坐标之间，没有明显的耦合，一个键稍微伸缩一下不至于使得另一个键长坐标上的受力受到明显影响。为了加快收敛，可以采用内坐标关键词（Opt=Z-Matrix）进行坐标转化。
内坐标不太适合有较多分子或原子的团簇体系，因为每个单体/原子之间内坐标怎么去链接比较合适往往说不清楚，很有可能会弄得很糟。内坐标尤为不适合环状体系的优化；内坐标对于高对称性体系的优化特别有用。通过巧妙地自行设定要内坐标，并结合一些虚原子作为辅助，可以大大减少要优化的坐标数目，不仅易于收敛也显著减少了计算梯度的耗时

单点能、结构扫描

单点能

定义：指对给定几何构型的分子的能量以及性质进行计算。在势能曲面上相对应的只是一个点，所以简称单点（能）计算。

用途：单点能计算可以用于：计算分子的基本信息；作为分子构型优化前对分子的检查；在由较低等级计算得到的优化结果上进行高精度的计算；在计算条件下，体系只能进行单点计算

关键词设置：设置计算的理论等级和计算类型 # HF/3-21G Pop=Reg SCF=Tight

这一部分需要出现的关键词有：

计算的理论：如HF（默认关键词，可以不写），B3PW91；
计算采用的基组，如6—31G， Lanl2DZ；
布局分析方法，如Pop=Reg；
波函数自恰方法，如SCF=Tight。

Pop=Reg只在输出文件中打印出最高的5条HOMO轨道和最低的5条LOMU轨道，而采用Pop=Full则打印出全部的分子轨道。SCF设置是指波函数的收敛计算时的设定，一般不用写，SCF=Tight设置表示采用比一般方法较严格的收敛。

内容参考自：单点能计算 ppt课件和 Gaussian软件应用——单点能计算

势能面扫描

针对此主题我还做了大量摘抄可供参考。

参考：详谈使用Gaussian做势能面扫描

刚性扫描

刚性扫描指的是让程序按照你指定的扫描设定依次改变体系几何结构，结构每变一次就算一次单点能，而不被扫描的那些几何变量还是保持在初始结构的值。

在Gaussian里做刚性扫描很简单，关键词里写上scan，然后对要扫的变量进行定义即可。比如某个键长变量原本定义是B1= 1.3，那么把这个变量设置改为B1= 1.3 10 0.1，就代表以1.3埃为初始值扫描10步，每步增加0.1埃，最终扫到1.3+10*0.1=2.3埃。由于初始结构也会被计算一次，所以总共会计算11个单点能。算完之后，可以从输出文件末尾读取信息汇总，也可以把输出文件用GaussView打开观看扫描过程每一步的结构变化，还可以在Result - Scan里图形化显示能量在扫描过程中的变化。

柔性扫描

柔性扫描中，每次按照你的要求而改变特定几何参数时，做的不是单点计算而是限制性优化，即没有被扫描的变量都会被优化以使得体系能量尽可能低，相当于允许这些变量自发地弛豫（relax）。

柔性扫描一般用于经验不足时，通过估计过渡态中断裂键的“长度范围”（或者二面角、夹角等参数等变化），按照一定的间隔（一般以0.1Å间隔）设置计算次数/steps，每次计算时都固定该断裂键的键长计算对应的单点能，然后确定其中最高的能量结构作为过渡态初猜结构。

注意：这里的过渡态结构也需要符合逻辑，一旦优化出现结构差异太大的情况，就需要手动调整初始结构中可能存在的问题。

在经验值积攒足够以后，就可以考虑缩小筛选范围了。

柔性扫描模板

这个模板是根据我们组的经验进行的调整，不一定适用于其他组。

扫描类型：

Lines in a ModRedundant input section use the following syntax:

[Type] N1 [N2 [N3 [N4]]] [[+=]value] [A | F] [[min] max]]

[Type] N1 [N2 [N3 [N4]]] [[+=]value] B [[min] max]]

[Type] N1 [N2 [N3 [N4]]] K | R [[min] max]]

[Type] N1 [N2 [N3 [N4]]] [[+=]value] D [[min] max]]

[Type] N1 [N2 [N3 [N4]]] [[+=]value] H diag-elem [[min] max]]

[Type] N1 [N2 [N3 [N4]]] [[+=]value] S nsteps stepsize [[min] max]]

Code	Meaning
A	Activate the coordinate for optimization if it has been frozen.
F	Freeze the coordinate in the optimization.
B	Add the coordinate and Build all related coordinates.
K	Remove the coordinate and Kill all related coordinates containing this coordinate.
R	Remove the coordinate from the definition list (but not the related coordinates).
D	Calculate numerical second Derivatives for the row and column of the initial Hessian for this coordinate.
H	Change the diagonal element for this coordinate in the initial Hessian to diag-elem.
S	Perform a relaxed potential energy surface Scan. Set the initial value of this coordinate to value (or its current value), and increment the coordinate by stepsize a total of nsteps times, performing an optimization from each resulting starting geometry.

其中S用于逐步扫描 F用于固定键长（冻结相关坐标）~~A用于反向激活（冻结体系，激活部分区域）~~

A我暂时还没用上，就不瞎扯了，等用上了再具体写。

from：https://flex.phys.tohoku.ac.jp/texi/gview/modred.htm

模版：

%mem=64GB
%nprocshared=16
# b3lyp/3-21G opt=modredundant nosymm 5d empiricaldispersion=gd3bj
# 基组可大可小，色散、5d等关键词该加还是要加 [空行]
Title Card Required
# [空行]
# [原子坐标]
# [原子坐标结束]
# [空行] 键长B：
B [No1] [No2] S [Steps] [Stepsize]
# 或者 键角A：
A [No1] [No2] [No3] S [Steps] [Stepsize]
# 或者 二面角D：
D [No1] [No2] [No3] [No4] S [Steps] [Stepsize]
# 固定键长的写法：
B [No1] [No2] [Bondlength] F

举例

单点能计算输入/输出文件

输入文件：

输出文件：

刚性扫描（扫描甲醇的羟基的二面角）

做扫描之前，我们先在B3LYP/6-31G*下对甲醇做一下几何优化，优化后是下图的结构：

在B3LYP/6-31G*下做扫描，想让C1-O5键转180度，使得上图的右图中H6从H2的反方向转到与H2重合。

通过gjf文件可知，D3对应H6-O5-C1-H2二面角，因此如果让D3变化，即绕着C1-O5转，就可以达到目的了。我们让这个变量每10度变化一次，因此总步数应当是180/10=18步。

因此对应的输出文件如下：

# b3lyp/6-31g(d) scan nosymm

Title Card Required

0 1
 C             
 H                  1            B1
 H                  1            B2    2            A1
 H                  1            B3    2            A2    3            D1    0
 O                  1            B4    2            A3    3            D2    0
 H                  5            B5    1            A4    2            D3    0

   B1             1.09337920
   B2             1.10125884
   B3             1.10125884
   B4             1.41863900
   B5             0.96872260
   A1           108.05956655
   A2           108.05956655
   A3           106.69693538
   A4           107.66742719
   D1           117.09085412
   D2          -121.45457294
   D3           180.00000000 18 10.

注意，凡是程序需要读入浮点数的地方，步长必须写为浮点数，如10.或者10.0而不能写整数10。

柔性扫描（乙醇脱水的势能面扫描）

假定H2与O8的距离发生变化

其中输入文件片段对应如下：

1
2
3

# PM7 opt=modredundant nosymm
...
2 8 S 17 -0.1

解释：最下行表示H2-O8距离缩短。一开始的结构H2-O8距离2.667埃，而水的O-H键长约0.96埃，因此我们可以从初始结构开始扫17步，每步减少0.1埃，最后一步对应2.667-17*0.1=0.967埃。

计算结果：

由于被扫描的变量是逐渐减小的，所以图应该从右往左看。随着扫描的进行，由于C-H和C-O键逐渐被拉长，能量逐渐升高，直到到达图中的最高峰。再往下走一步，能量突然降低，这是因为O-H距离此时已经比较短了，对其它变量进行优化后产生的结构已很像水+乙烯复合物了。凭基本化学直觉就能看出扫描曲线上最高点的样子差不多就是乙醇脱水的过渡态结构，即处于乙醇结构和水+乙烯复合物结构之间，因此可以用这个点作为找过渡态的初猜结构。

过渡态初猜结构需要看一下最高点对应的具体结构，如果结构不像过渡态，显然不能当做找过渡态的初猜结构。

对称性关键词

参考自：谈谈Gaussian中的对称性与nosymm关键词的使用

启用对称性的作用

没对称性的分子所属点群就是C1，有对称性的分子则属于比C1更高阶的点群。量化计算中电子积分的计算是很耗时的部分，对于有对称性的分子，要算的电子积分中会有很多是等价的，利用对称性的话就可以避免计算大量重复的积分，从而大大节约时间。一些计算在其它部分也能利用对称性显著节约时间，如积分变换、导数计算等，而且不光能节约时间，很多时候还能显著降低对内存/硬盘的占用。对称性越高的体系，利用对称性可以带来越多的加速。

启动对称性还有个好处是Gaussian可以显示出轨道、电子态、振动模式的不可约表示。如果不启用对称性，即总是当成C1对称性来算，即便体系本身有对称性，不可约表示也会都显示为A，此时轨道的不可约表示就只能基于等值面图判断、电子态的不可约表示需要自己再对占据轨道的不可约表示做直积得到、振动模式的不可约表示也只能基于正则振动坐标的箭头图来判断了，这就麻烦多了。另外，允许程序利用对称性，且正确判断出了实际点群的时候，还可以确保不同不可约表示表示的分子轨道在SCF过程中不发生混合，这对于考察一些问题以及SCF收敛有时会有帮助。

由于利用对称性时有上述优点，所以Gaussian默认是启动对称性的。计算一开始时会自动根据一定算法和阈值判断输入文件里分子的点群。

输入文件里的坐标称为输入朝向（Input orientation）。Gaussian为了判断和利用对称性，会先把输入文件里的坐标调整到标准朝向（Standard orientation）下，然后再进行点群的判断及后续的各种计算。这个调整过程会把体系进行平移，使原子核电荷中心处于笛卡尔坐标原点，并且根据一定规则对体系进行旋转。

将体系弄到标准朝向会带来一些麻烦的问题。有些时候，我们必须要求体系的位置和朝向是和输入文件里一致的。比如，我们要沿着某个键的键轴加用field关键词外电场，我们设好了外场矢量，但是如果体系朝向被Gaussian旋转到标准朝向了，那么外电场的方向就不合我们的意愿了（除非输入文件里的坐标已经是标准朝向的坐标），显然结果无意义。

令对称性算法失效的关键词：`nosymm`

理解了上面关于对称性的讨论，就很容易理解nosymm的带来的好处和坏处了。nosymm的含义是：让Gaussian在计算中完全不利用对称性。

关键词的用途：

使用后Gaussian就不会再试图去判断和利用对称性，因此也就不会把体系搞到标准朝向下，这就确保了实际计算过程中的坐标和输入文件里的坐标一致，从而可以避免前面说的诸如外加电场、作密度差图、CDA、ETS-NOCV分析时遇到的问题，所以做这些任务的时候总要加上nosymm。
做几何优化、柔性扫描的时候，Gaussian会把每一步的坐标都调整到标准朝向，这往往导致笛卡尔坐标变化不连续（但内坐标变化是连续的），导致在GaussView里观看优化、扫描轨迹过程中体系朝向会发生剧烈跳变，从而难以考察结构变化过程。解决方法有二，其一就是用nosymm，这样优化的每一步的结构都不会被弄到标准朝向下，而只会输出Input orientation坐标，gview也只会读这些坐标，优化/扫描过程的轨迹看起来就连续了。
nosymm还有个情况也应该加，也就是做对称破缺计算的时候（详见《谈谈片段组合波函数与自旋极化单重态》http://sobereva.com/82），如果体系本身是有对称性的，那么加上nosymm往往才能得到对称破缺态。比如将乙烯扭转到90度让pi键被破坏，此时是双自由基状态，如果直接用UB3LYP/6-31G* guess=mix可能因为波函数对称性的约束还是得到闭壳层波函数，而加了nosymm就可以得到双自由基状态了。

nosymm只有以上三个用处，如果你的情况不属于这三个，就别用nosymm！！！对于有对称性的体系，用nosymm带来的坏处是很明显的，首先是不能利用对称性来加速计算，其次是没法显示不可约表示。

几何优化中维持/破坏对称性

首先需要明白以下两个问题：

（1）在默认的允许利用对称性时，即便优化一开始把体系判断成了高对称性，优化过程中也可能对称性发生降低甚至变为完全没有对称性。

（2）在使用nosymm禁止使用对称性时，一开始的高对称性结构也依然可能继续保持下去。

所以，对称性是否会在优化中改变，或者说对称性是否会被限制在初始对称性下，和是否使用nosymm没必然的关系。

如果你想在优化过程中100%保证对称性肯定不发生改变，唯一做法是用内坐标描述结构，而且根据对称性来写内坐标，使得对称等价的坐标共享相同的变量（比如苯，让六个C-C键都共享一个键长变量），且优化的时候用opt=z-matrix让Gaussian在内坐标下优化。

对于有对称性的体系的计算，一般的考虑方式是：先让初始结构在gview里做对称化成为你期望的点群，然后用Gaussian优化，查看输出文件确保任务最开始时判断对了实际点群（如果没判断对，用symm=loose）。如果优化过程中点群发生了不希望的降低，且你用了弥散函数，把弥散函数去掉重新优化看看；对于DFT，可以进一步提升积分格点精度再试（如G09的用户可以加int=ultrafine试试，G16的用户可以加int=superfine试试）。如果重新优化后对称性还是出现了下降，换成其它你觉得对此体系也靠谱的泛函，也可以尝试不依赖于积分格点的MP2（前提是适合当前体系）或CCSD（量力而行）。

IRC

参考自：http://bbs.keinsci.com/thread-1906-1-1.html

介绍

因为TS-3在挣扎了好几次之后还有点问题，同学建议我跑一下IRC，所以查一下相关资料。

下面是一个最简单、典型的找TS+走IRC的过程：

在gview里凭借化学直觉将体系摆成尽可能接近于过渡态的结构，保存输入文件为TS.gjf。
将TS.gjf里的关键词设为# B3LYP/6-31G* opt(calcfc,noeigen,TS)，用Gaussian执行之
打开上一步任务的输出文件，保存成输入文件IRC.gjf
将IRC.gjf里的关键词设为# B3LYP/6-31G* IRC=calcfc，用Gaussian运行即可得到IRC。此任务出现报错几率很大，仔细耐心阅读完本文就知道怎么办了。

参考自：http://sobereva.com/400

呃，也就是说把opt部分关键词改成“IRC=”的关键词就行。

maxpoints和stepsize应该怎么设合适？这要看你跑IRC到底想干嘛，以下是笔者的建议：

验证过渡态找没找对：对于这个目的，IRC只要跑一小段就足够看出趋势了，也不需要IRC跑得质量多高。maxpoints用默认即可，stepsize可以设大到15或20，从而在不增加计算量的前提下能比默认时稍微跑长一些。建议加上LQA关键词避免HPC方法因校正步不收敛而中断。
获得近似的反应物和产物结构：一切同上，但是把maxpoints设50甚至更大，从而使IRC尽可能跑完整。
获得高质量、完整的IRC曲线：这主要是用于发文章目的，粗糙、不光滑的IRC曲线图放到文章中肯定会遭人嫌，我们也希望IRC尽量跑完整从而能够充分描述整个反应历程。为了让IRC质量高，可以把stepsize设小到5。由于设小了步长，又想要IRC完整，maxpoints必须设很大，比如200。如果此时跑出来的IRC还是不平滑，或者因HPC方法的校正步不收敛而中途报错中断，应再加上calcall重新跑。如果计算能力不足用不起calcall，也可以用比如recalc=3或5。

Q&A

Q：IRC 任务报错啦！末尾提示 Maximum number of corrector steps exceded 咋办？

A：这就是前面反复说的 HPC 方法校正步不收敛。有以下办法：

(1) 用 LQA，比如 IRC(calcfc,LQA)，由于此时不涉及校正步了，因此 100% 解决问题。但这容易造成 IRC 不够准确、不光滑。

(2) 如果你不希望用 LQA 在避免报错的同时牺牲 IRC 精度，则尝试减小步长（越小避免报错的几率越高），比如用 IRC(calcfc,stepsize=5)。或者加上 calcall，若嫌太昂贵就改用 recalc=x（x 越小避免报错的几率越高，但也越昂贵）。

(3) 改用 GS2 算法，即 IRC(calcfc,GS2)，可完全避免以上报错。耗时比 LQA 高得多，但精度也比 LQA 好。

(4) IRC 里加上 ReCorrect=never，这使得 HPC 方法不做校正步，故也完全避免以上报错。此时耗时和 LQA 相同，但所得 IRC 精度不如 LQA，因此强烈不建议用。

(5) IRC 里加上 maxcyc=N（N 应大于默认的 20）来加大 HPC 校正步迭代次数上限。笔者时常在菜鸟的 IRC 输入文件里看到这关键词，笔者强烈不建议用。虽然此方法不是说解决问题的可能性精确为 0%，但可能性实在甚微，没有试的必要。“不收敛就直接加大循环次数上限” 是菜鸟最常见的思维方式。

顺带一提，有时笔者看到有人的 IRC 输入文件里在 IRC 里用了 tight，这是用来把校正步收敛限设得更严的，用这个完全是莫名其妙，也不知道从哪里学来的。本来默认的收敛限下就容易不收敛，居然还给设得更严，明显会导致出现上述报错的几率大增。

Q：怎么 IRC 刚走了几步就正常结束了？怎么 IRC 走出来的两侧的曲线是相同的？

A：此问题是继上一个问题在网上被问得最多的与 IRC 有关的问题。出现这种问题都是因为优化过渡态时定位准确度不够。看下图，当优化出的过渡态位置不准确时，结构就不是在 IRC 的极大点了，而是稍微偏离一些的红球的地方

出现这种情况时，往右边产生 IRC 能正常产生，但是从红球位置往左边产生 IRC 时，还没怎么走，程序就发现能量升高了，误以为 IRC 已经走到了离极小点很近的位置，于是就不再继续走了，就正常结束了。还有一种情况，是刚往左边走 IRC，由于体系受力是冲着右边的，导致马上转了个弯就往右边走了，就呈现了 IRC 左右两边曲线都一样的结果。

对这个问题，应按照以下方式排查和尝试解决

(1) 先确保初始结构是之前优化 TS 得到的结构，而且过渡态优化和走 IRC 都是在严格相同级别下进行的。

(2) 提高过渡态定位精度。在找过渡态时候用 tight，对于 DFT 再同时结合 int=ultrafine（此时产生 IRC 也必须用 int=ultrafine）。如果还不行，优化过渡态时用 calcall（或者用诸如 recalc=3）。

(3)如果反复尝试了 (2) 的方法还是不行，或者你不想尝试(2)，毕竟会增加很多耗时，那也可以尝试增大 IRC 步长，比如 20 乃至 30。由于步长大了，从上图红球的位置往左走的时候可能一下子就越过了 TS，之后就能正常继续往左产生 IRC 了。不过步长大了容易导致 HPC 校正步不收敛、IRC 不准确不平滑等问题，怎么考虑和处理前面已经说了。

另外，出现这种问题还有一种可能是在 IRC 任务中，基于自动初猜的波函数做 SCF 后收敛到的波函数与找过渡态任务最终得到的波函数不同，此时相当于 IRC 任务所在的势能面和过渡态搜索任务所在的势能面不同，这也会导致 IRC 异常，因为类似于违背了前述的走 IRC 的 “任何影响势能面的设定必须严格相同” 的这个前提。出现这种情况时，你会发现 IRC 任务第一次输出的 SCF Done 能量和找过渡态最后一步的 SCF Done 能量明显不同。为解决此问题，走 IRC 的时候可以用 guess=read 关键词，从优化过渡态的 chk 文件中读取最后的波函数（并且最好用 forward 和 reverse 关键词通过两个任务分别跑正向和逆向 IRC），这样通常可以确保 IRC 任务所在的势能面和优化过渡态时相同。

用 SMD 溶剂模型时，也可能个别时候由于数值噪音问题出现 IRC 走几步就停了的现象。可将优化和走 IRC 用的溶剂模型都改为 IEFPCM 再试，说不定能解决。

Q：我的 IRC 怎么成这样了？怎么有个点跑到下头去了？

A：即便 IRC 任务还在算着，也可以用 gview 打开输出文件看看走到第多少步了，检查是否正常。当 IRC 还没有走完的时候就打开输出文件，或者看的是 IRC 失败的任务的输出文件，往往就会看到中间有一个点突然掉下去的现象。任务正常结束后再看就不会有这个现象。

Q：走出来的 IRC 最两端的点是反应物和产物结构么？

A：否。IRC 和几何优化不一样，几何优化的步长是随机应变的，而 IRC 的步长是固定的。从 TS 开始走，即便把 maxpoints 设得非常大，stepsize 设得很小，来让 IRC 跑得又完整质量又高，也注定不可能恰好有一个 IRC 点正好落在与反应物或产物对应的极小点的位置。因此，必须对 IRC 两端的点进一步做几何优化才能得到准确的反应物或产物结构。

Q：为什么我看 IRC 趋势是对的，我对 IRC 两边的点优化，优化得到的结构却不是我想要的极小点结构？

A：有两种可能：

(1) 与这个 TS 直接连通的极小点本来就是那样。由于你的 IRC 跑得还太短，在加上你的化学直觉不准，导致并没有估计对真正极小点结构。可以让 IRC 尽量跑完整一些来检验。

(2) 优化的时候由于走得步子太大，或者 Hessian 矩阵不准确，导致优化到了其它邻近的极小点去。此时可以在优化时用较小的步长上限并结合 notrust（相关信息看《量子化学计算中帮助几何优化收敛的常用方法》），或者结合 calcall/recalc=x 提供精确 Hessian 矩阵，或者先把 IRC 跑得尽量长一些，使两端结构尽可能接近极小点时再做几何优化。

Q：柔性扫描和 IRC 有什么区别？柔性扫描能代替 IRC 么？

A：二者概念完全不同，不能代替。IRC 是从 TS 开始顺着虚频方向，沿着梯度负方向移动结构走出来的轨迹，所有坐标是同步变化的。而柔性扫描，是令某个（或某些）坐标以特定方式不断改变，每次都优化其余所有变量所得到的（限制性优化），故所有变量间的变化不是协同的，得到的能量变化曲线与 IRC 也是大相径庭的。柔性扫描曲线及其轨迹经常会出现突跃，而不像 IRC 那样总是平滑变化。不过，当某个化学过程主要只涉及一个变量发生变化的话，柔性扫描这个变量和走 IRC 产生的曲线倒是比较相近，此时如果你用的计算级别没有解析 Hessian 而不容易走 IRC 的话，用柔性扫描也未尝不可（比如当年没有 TDDFT 解析 Hessian 的时代，由于不方便找 TS 走 IRC，对激发态质子转移过程的研究经常是靠柔性扫描来代替）。

Q：我的 IRC 曲线在过渡态位置怎么特别尖 / 不平滑，貌似不对劲？

A：应当确认初始结构是优化过的过渡态结构，并且找过渡态和走 IRC 用的计算级别（包括任何能够影响单点能的因素）严格相同。如果你就是这样做的，那么尝试增加过渡态的定位精度，比如优化过渡态时用 tight 收敛限、calcall/recalc=x、增加积分格点精度等。如果这些也都试了，走 IRC 时用 GS2 算法，或加上 calcall 试试。另外也要确保 IRC 任务第一步收敛到的波函数与过渡态搜索最后一步相同，这在前面已经说过了。

Q：IRC 任务如何续算？

在 IRC 里添加 restart 关键词即可，例如 IRC(restart,maxpoint=15) b3lyp/6-31g(d)。此时 %chk 应当和之前中断的 IRC 任务相同。

Q：我发现之前 IRC 跑得不够长，如何在不完全重算的前提下在两侧再增加一些点？

A：比如 maxpoints=15 的设定下 IRC 任务已经正常跑完，现在还想两个方向各增加 5 个，执行 IRC(restart,maxpoint=20)，%chk 应当和之前的 IRC 任务相同。

注意 gview 观看 restart 得到的 IRC 轨迹要打开 chk/fch 文件，而不要打开 out/log 文件，否则看到的只是续跑出来的点。

Q：IRC 任务出现 L502 错误怎么办？

A：这是走 IRC 过程中出现了 SCF 不收敛所导致，按照常规解决 SCF 不收敛的方法尝试添加适当的关键词解决，见《解决 SCF 不收敛问题的方法》（http://sobereva.com/61）。另外读者也可以尝试降低 IRC 步长，因为每一步初猜波函数沿用上一步收敛的波函数，显然如果两帧结构相差越小，那么上一步收敛的波函数就是下一步越好的初猜波函数。

Q：我的 IRC 看起来比较怪，合理么？怎么边上还凸起来一块？（下图为笔者随便找的例子）

A：很正常。只要确信自己计算流程合理，而且轨迹看着无异常，TS 位置曲线也平滑，就没有问题。出现诸如肩峰往往暗示这个反应涉及的电子结构变化复杂，可能整个反应虽然是一步，却呈现出两个阶段，应当结合结构变化特征、波函数分析手段（比如用 Multiwfn 分析反应过程中键级、ELF、密度差、原子电荷变化等）予以考察。

Q：我的理论方法支持解析梯度但不支持解析 Hessian（比如 CCSD），怎么跑找 TS 并且跑 IRC？

A：先说此时怎么找 TS。有两类做法：

(1) 纯依赖于梯度的方法

用 CCSD/cc-pVDZ opt(TS,modRedundant,noeigen) 或者 opt(TS,gediis,noeigen) 或 QST2/3，这些过渡态搜索关键词都不需要提供初始 Hessian

(2) 基于半数值 Hessian 的方法

把初猜 TS 结构的 Hessian 矩阵存到 chk 文件里：# CCSD/cc-pVDZ freq（这步耗时很高）

从 chk 里读 Hessian 矩阵并优化 TS：# CCSD/cc-pVDZ opt(TS,noeigen,readfc)

找到 TS 结构后，使用这个关键词即可生成 IRC：CCSD/cc-pVDZ IRC(gradientonly,euler)。这里 gradientonly 代表使用只依赖于梯度的 IRC 算法，默认是 DVV 方法，由于结果烂到爆，所以改成相对好点的 Euler 方法，但质量也不怎么样，经常震荡得很厉害。

Q：我怎么对 IRC 过程做波函数分析，考察电子结构随反应过程的变化？

A：参看《产生 Gaussian 的 IRC 和 SCAN 任务每个点的波函数文件的工具》、《通过键级曲线和 ELF/LOL/RDG 等值面动画研究化学反应过程》。

Q：我之前对过渡态已经做过振动分析了，chk 里已经有了过渡态结构下的 Hessian 矩阵，能否在 IRC 计算时直接读取之，避免用 calcfc 计算 Hessian 以节约时间？

A：完全可以，对较大体系也推荐这么做。在 IRC 里加上 rcfc 关键词即可，此时就不需要写 calcfc 了，程序会从 %chk 文件中直接读取 Hessian 矩阵来用。

计算绘图建议

所有结构都需要自己画

首先提出的第一个意见是不要偷懒，就算ppt中有所有物质的结构，也需要自己画，这样有助于理解过程中自由基到底进攻了哪些结构，也可以避免初始结构出现问题。

建立合理的重命名习惯

接下来的意见是需要建立合理的重命名习惯。

比如可以命文件名为：任务类型-结构名称-序号（是否符合要求）这样的思路命名。

如老师的命名思路：ts-1-opt就是没有问题的最终优化的过渡态结构；

再如扫描中使用的输入文件可以命名为：scan-ts-1这样命名。

在区分不同结构时，可以按照abcd的方式添加序号，如ts-1a这样命名。

通过下一步中间体反推节约手动调整结构的时间

在中间体结构手动优化过程中，可以通过下一步的优化结构反向考虑需要从哪些键进行结构调整，从而避免计算过渡态时初始结构出现问题的情况。

随时需要监控计算步数和时间

对于小规模，尤其是本课题中这种没有涉及金属的结构，理论计算时间和步数都不会太长，因此计算超过3小时的时候就有必要从服务器调取out文件查看当前计算中的结构的情况和前面step的情况，观察结构是否处于震荡无法收敛，从而及时止损提高计算效率。

out文件中关注内容

搜索SCF done关键词，在这一关键词附近存在一些重要数据，可以反映当前优化结构是否符合常理，另外也可以查看当前状态优化的电子能。
查看优化步数是否合理，查看最后一次和倒数第二次优化时间确定实际优化时间。
查看结构，观察目前获得的结构中是否存在H-H等相互作用的可能性（H-H键长<2.4则存在相互作用），存在则需要调整再算。

Chemdraw中应使用Arial字体和ACS Document 1996格式

前面画的图都没有按照格式画，所以又重新传了一个版本，以后套用相关格式（在File中）

所有相关说明和字体都应该使用Arial字体。

另外关于单C上面连接两个F的情况，可以先接一个三元环，然后把环多余的键去掉处理：

绘图建议

以下内容基于ACS Document 1996格式：

自由基单点距离要稍微远一点，大小应放大到120%～130%
化学式一般是一对一的，没有必要框在一起
页边距不要太小，布局应该稍微宽敞一些
化合物结合时势能面可加一个小弧线；反之物质生成时则加一个带箭头的弧线
箭头使用最右侧的小箭头即可，大小不必太大
箭头右侧标记的方式如图所示
物质编号放在化学式下方，离势能面较近时可省略
过渡态加框和双线号（Double Dagger）
保持化学式和结构的居中，可利用Group成组功能便于对齐
断键单键和双键虚线应保持“连接”状态
数据的数字部分离横线距离可以稍微远一点，不用挤在一起

任务执行思路

确定中间体结构和将要发生变化的结构（本课题一开始就确定了）；
利用小基组（如3-21G）优化，进行逐步扫描（以0.1Å间隔设置计算次数）；
在计算过程中提取文件确定上述结构是否合理；
确定合适的结构后改成过渡态输入文件投入计算；
验证输出结构的构象是否存在其他相互作用，继续改进结构直到没有可以改的地方。

任务执行情况

3月22日

~~掰了一下CHO-1的结构，结果没有太大区别~~

~~首先看了一下大小，就觉得优化的没啥用了。~~

~~打开看过以后发现形状没啥差别，分子轨道能量也基本上一样。~~
2. 掰了一下Radd-TS1的结构。

3月23日

分析新算的Radd-TS1-2的out文件。

~~首先还是看了一下大小，感觉掰的好像有点用：~~

~~之前的结构中，CHO-1与另一部分的距离稍微有些大，所以主要调整了图中原子的距离和受其影响不够合理的二面角：~~

~~CYLview的优化曲线：~~

~~看看调整以后：~~

~~可以看到两个结构已经完全连在一起了。~~

当然这也有可能是因为两个部分本身距离就调的太近，使GaussView认为两者已经成键导致的。

因此，Gauss View一个非常关键的点就是：不能把键标的太当真。

~~CYLview的优化曲线：~~

~~第一次优化时，因为两个分子距离较远，所以优化到最后也没有成键，而第二次优化，因为一开始就默认两者成键，所以显得平稳了不少。~~

~~这个结构算的方法到底对不对，我得问一下同学再更新了。~~

~~画并且掰一下成环中间体2（rts-2）和开环中间体3（rts-3）的结构，丢进去算一下。~~

~~就懒得调了，先让软件优化一次试试。~~

~~非常推荐用snippaste的固定功能把结构式放在下面，然后上面用画好的结构对照一下。~~

~~记得看看中间体的多重态状况，一般是doublet。~~

3月26日

やられた。

记录一下刚才同学教我的计算上存在的问题。

首先，Radd-TS1是过渡态，而我目前需要计算的是中间体。

相关的知识介绍我放在上一节介绍了，总之错误添加关键词并且拉近C-C键后的情况是：

频率分析：

可以看到还是存在虚频。

这是对应的振动方式：

很明显，中间体不能存在这种形式的振动。

因此，接下来需要重新按照# opt freq b3lyp/6-31g(d) 5d empiricaldispersion=gd3bj投入三个中间体的任务。

4月3日

这个星期主要在折腾这个云写作，所以一直到周日上午才有时间继续折腾。

今天的主要任务是掰结构和尝试算一下三种中间体的过渡态（RTST1、RTST2、RTST3）。

~~中间体RTS-1过渡态RTST1-1：C+自由基处的C-C键断键可能性较高，所以改一下键长。~~

~~键长：1.49→1.95（标准1.54）。~~

~~中间体RTS-2过渡态RTST2-1：连接苯基和氧自由基的C-C键发生断裂，改一下这里的键长。~~

~~键长：1.55→1.85（标准1.54）。~~

~~中间体RST-3过渡态RTST3-1：虽然没有写断键，但是含C+自由基的C-C键断键可能性高，就搞这个了。~~

~~键长：1.49→1.85（标准1.54）。~~

~~然后改输入文件，改成之前同学给我的那个版本。~~

~~算完以后查一下单点计算的算法。~~

~~暂报：~~

~~RTST1-1 结构会散开，不知道对不对，很怪；~~

~~RTST2-1 看上去好像是对的，（环的振动更加强烈）但是不知道怎么表达；~~

~~RTST3-1 整个结构都在振动，不知道是不是C-C键改的不够长，加到2.0看一下。~~

老板亲自找我了，我才知道这三个结构基本上全都算的有问题。

4月8日的部分我会把相关内容全部重新盘点一下。

4月8日

前面算的rtst系列结构全部都有问题，其中只有2-1一个结构稍微没那么大问题。

这是因为这两个结构的初始结构就有问题。

其中2-1的结构按照老师的要求重命名为TS-3a（本质上是开环过程也就是TS-3的过渡态）。

重画全过程，重新命名

图上面有，复读一下：

基础结构错误需要重新开始的（TS-1、TS-2）【扫描】

TS-1的情况比较简单，就是CHO-1的双键和A1偶联后，也就是用产物A2直接进行调整就可以获得：

如上图，加粗的键通过设置扫描范围为【1.6-2.5Å，0.1Å步距，9步】进行柔性扫描并确定TS-1（初猜）结构。

（右键-view-labels查看对应的原子编号）

通过编号需要调整的C-C键编号为C28-C22

输入文件和结构：

# 文件头
%mem=64GB
%nprocshared=16
# b3lyp/3-21G opt=modredundant nosymm 5d empiricaldispersion=gd3bj
scan_ts1
# 略

B 22 28 S 9 0.1

# 文件尾

文件名：scan-ts1a；

TS-2则正常通过该过程的键变化来确定：

如上图所示的A2结构中，TS-2是a的C自由基进攻e上C的羰基（C24-C20）的结果。

因此，通过结构参数，使C24-C20距离在一定范围以下，然后再通过【-0.1Å步距】反向进行扫描，并尝试确认其中的结构。

TS-2的骨架还是可以利用刚才的scan用结构直接改造。

此处的2个碳（C24-C20）距离还略远，因此调整：

调节（12-14-17-20）的二面角（-180～-30左右）键长从4.4缩短到2.63

调节（17-20-22-28）的二面角（～88左右）使过短的分子内C-F距离拉长

这里能隐隐约约看到一些五元环的雏形，可见反向推理也是可以解决问题的。

另外，输入文件代码文件头与TS-1相同，只改一下步数和步距：

B 20 24 S 10 -0.1

文件名：scan-ts2a；

提交，下午看看结果。

计算正确但部分细节需要优化的（TS-3）

老师指出，这里的TS-2a的结构中还是存在一个H-H距离较近的情况（>2.4时才符合需求）。

因此通过旋转5元环的朝向（图中4个原子的二面角），使得不再出现距离较近的情况。

将这个文件保存为ts-3b，再投入计算。

4月19日

反应过程重画固定键长重算

最近几天没更新，写一下情况：

这里的化学式还需要调整成和势能面一样的形式。

首先，基本完成了TS-2、TS-3的能量优化，TS-1的C-C键通过多次优化，确定了键长应该在2.268左右，于是进行了结构微调，并在最后添加了冻结键长的关键词：

# b3lyp/3-21G opt=modredundant nosymm 5d empiricaldispersion=gd3bj

ts1fixed

[Describtion of molecule]

22 28 F
_
_

前面进行了几次失败的尝试，所以最后命名是TS-1d4和TS-1d5。

单点计算

本实验的溶剂是DCE（1,2-Dichloroethane 1,2-二氯乙烷）
因此使用SCRF=(SMD,SOLVENT=Dichloroethane)关键词

用# b3lyp/6-31g(d) scrf=(smd, solvent=dichloroethane) 5d empiricaldispersion=gd3bj分别计算各个产物的液相单点能。

计算完并记录好数据如下：

因为按照模板来处理了，所以效率高了不少：

前面说了，高斯默认单位Hartee需要换算为电子伏特/eV才能使用，换算系数为：1 hartee to eV = 27.2114 eV。

不过，在制作势能面的时候，单位用的是kcal/mol，因此还需要换算一下：

1 hartee = 627.5 kcal/mol
1 hartee = 27.2114 eV
1 hartee = 2625.5 kJ/mol

目前做出来的图是这样的：

现在主要问题是TS-1这里吸热却熵减，其他从能量角度来讲没啥矛盾，就看反应方面有没有啥问题了。

势能面

势能面图中的能量是各个结构在溶剂中计算的单点值。

因为整个过程的实际起点是TS-1，基础构架可以这样画（能量暂时没填进去）：

经过老师的建议，上述势能面图调整并填写能量后如下：

不过，上图中TS-2和A3的能量老师指出有问题（稳定态A3比TS-2能量高），~~准备重算。~~

结果并不是需要重算的问题，而只是A3的数据拷贝错了。

另外老师给了一个更加科学的能够处理较大数据量的数据填写表，示例如下：

这样最后得到的势能面：

总之，虽然又犯了一堆傻逼错误，最后还是把东西都算出来了。相关经验我也会更新回模板中。

5月2日

主要补充一些溶剂化计算数据问题和IRC复算问题。

高精度SMD复算

因为需要换成更高精度的基组，所以换成：

# m06/6-311+g(d,p) scrf=(smd, solvent=dichloroethane) 5d

因为m06泛函只能使用0阻尼色散，所以后面改成empiricaldispersion=gd3

（不过后面老师又说不需要加色散，所以又改了一个版本，参考上方）

不过又出现新问题，这种情况算的时候TS-3和A3的能量有些问题：

高精度复算的图就不写进来了，正好当作保密吧。

TS-3照理说应该比A3能量高，所以有必要稍微算算TS-3的IRC看看有没有什么问题（验证？）

IRC验证

TS-3有必要复算，所以参考上面摘抄的玩意写一下输入文件：

# irc=(calcfc, LQA, stepsize=20) b3lyp/6-31g(d) 5d empiricaldispersion=gd3bj

试试看有没有报错，以及报错情况。

首先，IRC不能和freq一并计算，所以上面的关键词不能带freq

去掉这个以后，报错确实没有，曲线不怎么光滑。

可见IRC的输入还是：

# irc=calcfc b3lyp/6-31g(d) 5d empiricaldispersion=gd3bj

这样写比较好。

不过从IRC的曲线趋势判断的话，照理说TS-3目前的结构应该已经算是最高点了，也就是说这个能量关系应该不会有很大问题。

5月25日

在同学的建议下，需要优化中间体和过渡态的结构的构象。

从这个例子可以看出来构象优化时比较重视的特点：

官能团交叉，优化前五元环和苯环基本上在同一平面（四面角很小）
C-C键一类的结构需要互相交叉（优化前有些偏向，优化后变为交叉）

通过优化，过渡态能量变得比较合理了：

现在结构加了网格，因为可以更好的确定能量变化情况，另外物质名称按照同学组会PPT的内容做了点改变。