用高斯做势能面扫描（一）：刚性扫描

转自公众号“量子化学”

势能面扫描可用于研究一个区域内的势能面，通常是在不同结构上执行一系列的单点能计算，得到能量或其他性质随几何结构变化的情况。势能面扫描分以下两种类型：

刚性扫描 (rigid scan)只改变被扫描的变量，而分子内其他变量固定不变，对分子进行一系列的单点能计算。柔性扫描 (relaxed scan)改变被扫描的变量，在每个结构处固定当前变量，并优化分子结构（实际是一个限制性优化）。

本篇先介绍刚性扫描的做法。刚性扫描使用 scan 关键词，与一般计算不同的是必须使用内坐标 (internal coordinate) 或称 Z-matrix 来描述分子结构。首先介绍一下内坐标的书写方法。

内坐标书写

在内坐标中，使用原子间的相对位置来定义原子坐标。具体来说使用的是键长、键角和二面角三个量。

以丁烷分子为例：

用 GaussView 构建完分子后，在保存 gjf 文件时，可以将 Write Cartesians 选项前的勾去掉，就可以保存体系的内坐标。如下所示：

 C  
 H    1   B1 
 H    1   B2    2   A1 
 H    1   B3    2   A2    3   D1    0 
 C    1   B4    2   A3    3   D2    0 
 H    5   B5    1   A4    2   D3    0 
 H    5   B6    1   A5    2   D4    0 
 C    5   B7    1   A6    2   D5    0 
 H    8   B8    5   A7    1   D6    0 
 H    8   B9    5   A8    1   D7    0 
 C    8  B10    5   A9    1   D8    0 
 H   11  B11    8  A10    5   D9    0 
 H   11  B12    8  A11    5  D10    0 
 H   11  B13    8  A12    5  D11    0   
 
   B1            1.06999912   
   B2            1.06999995   
   B3            1.07000079   
   B4            1.54000051   
   B5            1.07000043   
   B6            1.06999979   
   B7            1.54000000   
   B8            1.06999991   
   B9            1.06999952   
   B10           1.53999941   
   B11           1.07000000   
   B12           1.07000073   
   B13           1.06999983   
   A1          109.47125498   
   A2          109.47121447   
   A3          109.47117182   
   A4          109.47119590   
   A5          109.47120631   
   A6          109.47119171   
   A7          109.47118682   
   A8          109.47122479   
   A9          109.47122429   
   A10         109.47121247   
   A11         109.47122973   
   A12         109.47123833   
   D1          120.00006239   
   D2         -120.00003390   
   D3          179.88897367   
   D4          -60.11103467   
   D5           59.88896396   
   D6         -120.00003791   
   D7          120.00000474   
   D8           -0.00000000   
   D9          180.00000000   
   D10          60.00005254   
   D11         -60.00471913   
   

内坐标的书写格式为：

元素  原子1  键长  原子2  键角   原子3  二面角

表示当前原子与原子 1 的键长，当前原子与原子 1、原子 2 形成的键角和当前原子与原子 1、原子 2、原子 3 形成的二面角。这三个“参考原子”必须在当前原子之前已经定义好。它们的选取具有一定的随机性，因此，一个体系的 Z-matrix 可以有很多种。在上述例子中，使用的是变量的表示方式，将数据先用符号表示，再在最后一个原子之后空一行给出各变量的值。也可以直接将数值写在 Z-matrix 中。在用 GaussView 自动生成的内坐标中第四行开始的最后还有一个数字 0，这个数字在 ONIOM 计算中会遇到，而在一般的计算中可以不写，或写 0。

计算实例

以下我们以一个实例来说明如何使用 scan 关键词进行势能面扫描。下图是邢其毅《基础有机化学（第三版）》上册中的“正丁烷各种构象的势能关系图”。我们尝试绘制这幅图。

其输入文件如下：

%nproc=48
%mem=20GB
%chk=C4.chk
#p m062x/6-311G** scan nosymm

n-butane scan

 0 1 
 C    
 H    1     B1  
 H    1     B2    2      A1  
 H    1     B3    3      A2    2      D1    
 C    1     B4    4      A3    3      D2   
 H    5     B5    1      A4    4      D3   
 H    5     B6    1      A5    4      D4  
 C    5     B7    1      A6    4      D5   
 H    8     B8    5      A7    1      D6  
 H    8     B9    5      A8    7      D6   
 C    8    B10    5      A9    6      D6  
 H   11    B11    8     A10    9      D7   
 H   11    B12    8     A11    9      D8   
 H   11    B13    8     A12    9      D9       
 
 B1            1.07000000   
 B2            1.07000000   
 B3            1.07000000   
 B4            1.54000000   
 B5            1.07000000   
 B6            1.07000000   
 B7            1.54000000   
 B8            1.07000000   
 B9            1.07000000   
 B10           1.54000000   
 B11           1.07000000   
 B12           1.07000000   
 B13           1.07000000   
 A1          109.47120255   
 A2          109.47121829   
 A3          109.47121829   
 A4          109.47120255   
 A5          109.47123134   
 A6          109.47120255   
 A7          109.47120255   
 A8          109.47123134   
 A9          109.47120255   
 A10         109.47120255   
 A11         109.47120255   
 A12         109.47123134   
 D1         -120.00000060   
 D2          120.00003407   
 D3          -60.11108461   
 D4           59.88890799   
 D5          179.88890060   
 D6         -180.0 36 10.00   
 D7          -60.00000740   
 D8          180.00000000   
 D9           59.99527604   
 

使用的关键词为 scan ,对需要扫描的变量，将其写为

变量名  初始值  扫描步数  步长

其中步数是指改变的步数（需要写成整数），实际计算的结构要多一个初始结构。也可以对多个变量进行扫描，则结构的数目为每个变量扫描步数的乘积。

一方面为了方便观察轨迹，另一方面防止扫描过程中对称性改变带来的问题，一般在势能面扫描时加上 nosymm。

在本例中，输入文件中的内坐标与高斯自动生成的坐标不同。实际需要做的是旋转C5-C8单键，对应的则是改变 1、5、8、9 四个原子形成的二面角，但并不能只改变这一个二面角，在 C5-C8 单键旋转时，C8 上的 H9、H10 和 C11 需要同时进行旋转，因此需要修改这三个原子的二面角的定义，也设置为 D6。

最后，C11 上的三个氢原子也要注意随着 C11 一起运动，也要对它们的坐标进行修改。一个简单的修改方法是判断好该原子需要与哪些原子保持相对位置的不变，则用它们作为参考原子。最终内坐标的写法如上所示，得到的扫描轨迹如下：

计算结束后，在输出文件中会有各步能量的汇总：

 Summary of the potential surface scan:   
    N      D6          SCF  
....   ........   ............
   1  -180.0000   -158.39616  
   2  -170.0000   -158.39444   
   3  -160.0000   -158.39212   
   4  -150.0000   -158.38918   
   5  -140.0000   -158.38552   
   6  -130.0000   -158.38200   
   7  -120.0000   -158.38047   
   8  -110.0000   -158.38200   
   9  -100.0000   -158.38554  
  10   -90.0000   -158.38921   
  ...

可以将数据提取到作图软件中绘制势能面，以可以在 GaussView 中观察。点击 GaussView 中的 Results -> Scan 即可得到如下势能曲线：

各结构出现的顺序与教材略有不同，但图象上几个关键点是一致的。(1)为能量最高点——全重叠构象，(4)为能量最低点——对交叉构象。在笔者所得的结果中，(1)、(2)、(3)与(4)的能量差分别为 47.4、3.9、14.2 kJ/mol，后两个数与教材中的值吻合，而第一个能垒偏大。于是笔者将(1)的结构作为初始猜测，进行了 opt=ts 的优化，得到的结构与(4)的能量差为 21.8  kJ/mol，与教材中的值吻合。由此也可以看出，势能面扫描可以用于辅助寻找过渡态，势能面上的极大点往往是过渡态的一个很好的初始猜测。

---

学习更多固体与表面理论计算相关内容，欢迎参加首届固体与表面理论计算（固体场）中级培训班（北京，20191017-20191021）；详情请见：第一届固体与表面理论计算（固体场）中级培训班（北京，10月17日-10月21日）。咨询请联系科老师，手机：[phone]（微信同号）

相关内容链接：
1、从「电子结构」开始，教你如何解析理论催化计算结果2、研之成理理论催化实战教学 —— 电子结构分析二3、研之成理理论催化实战教学 —— 电子结构分析三4、研之成理理论催化实战教学 —— 电子结构分析四5、研之成理理论催化实战教学 —— 电子结构分析六COOP6、研之成理理论催化实战教学 —— 电子结构分析七COHP

更多理论计算相关内容，请在公众号内搜索 “理论”查看全部内容。