用高斯做势能面扫描（二）：柔性扫描

Jun 28, 2020 增加保存数据的操作方法及其他细节
Sep 5, 2019 初版

柔性扫描简介

势能面扫描可以用来研究势能随少数几个坐标变化的情况，分为刚性扫描和柔性扫描两种。刚性扫描是指被扫描的坐标变化后，不改变未被扫描的坐标，相当于做一系列单点能计算。对于刚性扫描，本公众号之前已经有过介绍，见《用高斯做势能面扫描（一）：刚性扫描》。柔性扫描是指被扫描的坐标变化后，固定被扫描的坐标优化分子结构，相当于做一系列限制性优化。刚性扫描所需的计算资源更少，但是由于刚性扫描只考虑了一个自由度，其计算结果对势能面的描述不如柔性扫描更真实。

做柔性扫描时高斯输入文件的书写

柔性扫描可以使用任意坐标描述分子结构（内坐标或笛卡尔坐标），由于要做一系列限制性优化，routine section中要写入opt=modredundant，一般还要写入nosymm，最后在输入文件末尾写入被扫描的坐标及其步数和步长。

柔性扫描时文件末尾内坐标的描述方式为

[Type] N1 N2 [N3 [N4]] S steps step-size

中括号代表可选项。Type代表内坐标的类型，写入B代表键长，A代表键角，D代表二面角。NX代表第X个原子在这个分子结构中的序号，当只有N1和N2存在时代表扫描N1-N2间的键长，N3存在时代表扫描N1-N2-N3间的键角，N3和N4都存在时代表扫描N1-N2-N3-N4间的二面角。[Type]为可选项，因为内坐标的类型可以通过给出原子的个数确定。S为保留字，表明这是个柔性扫描作业。steps代表步数，step-size代表步长（需写浮点数），这两项与描述刚性扫描时的输入项相同。被扫描变量的初始值即为当前输入文件中该变量的值，因此需要注意将初始结构调整为想要的结构。

例如，给出一个柔性扫描作业的输入文件，文件尾写着

2 3 S 3 0.05

代表扫描原子2和原子3间的键长，扫描3步（加上起点一共得到4个结构），每步键长增加0.02埃。

笛卡尔坐标同样也可以用于柔性扫描，不过此时扫描的不再是内坐标，而是笛卡尔坐标。当用笛卡尔坐标描述柔性扫描时，输入文件结尾应写入

X atom S x-steps x-size y-steps y-size z-steps z-size

其中X为保留字代表笛卡尔坐标，atom代表原子在分子中的号，x-steps和x-size代表原子在x方向上的步数和步长，y方向和z方向上同理。

例如，给出一个笛卡尔坐标柔性扫描的作业，文件尾写着

X 1 S 5 0.2 0 0.0 3 0.1

这代表相对起始结构，原子1在x方向移动5步，每步0.2埃，在z方向上移动3步，每步0.1埃。

注意，柔性扫描时，扫描内坐标并不要求分子结构用内坐标表示，扫描笛卡尔坐标不要求分子结构用笛卡尔坐标表示。与刚性扫描要求分子结构用内坐标描述不同，柔性扫描对分子结构的描述方式没有要求。

下面用两个实例说明柔性扫描的输入文件的书写和计算结果的分析。第一个例子为一维势能面柔性扫描，第二个例子为二维势能面柔性扫描。

一维势能面柔性扫描

为了与刚性扫描一文对比，本文同样用丁烷分子绕C-C键旋转过程为例。丁烷分子的结构如下所示，由丁烷分子的三维结构可知，被扫描的内坐标应该选取1-5-8-11二面角。

丁烷分子的三维结构

其输入文件如下：

%nproc=48
%mem=20GB
%chk=C4_relaxed.chk
#p m062x/6-311G** opt=modredundant nosymm

n-butane relaxed scan

0 1
C      0.00000000    0.00000000    0.00000000
H      0.00000000    0.00000000    1.07000000
H      1.00880000    0.00000000   -0.35670000
H     -0.50440000   -0.87370000   -0.35670000
C     -0.72600000    1.25740000   -0.51330000
H     -0.72430000    1.25840000   -1.58330000
H     -1.73530000    1.25640000   -0.15830000
C     -0.00160000    2.51480000    0.00230000
H     -0.50600000    3.38850000   -0.35440000
H      1.00770000    2.51580000   -0.35280000
C     -0.00410000    2.51340000    1.54230000
H      0.49920000    3.38700000    1.90060000
H      0.50030000    1.63970000    1.89900000
H     -1.01340000    2.51250000    1.89740000

D 1 5 8 11 S 36 10.00

关键词opt=modredundant表明每一步都需要做限制性优化，写入nosymm一方面是为了方便观察轨迹，另一方面防止扫描过程中对称性改变带来的问题。D 1 5 8 11 S 36 10.00表明扫描原子1-5-8-11构成的二面角，扫描步长为36步，每一步二面角增大10度。因为步长是36，再加上初始结构，我们最后一共会得到37个结构。

得到的扫描轨迹如下：

柔性扫描轨迹

计算结束后，输出文件中有每一步能量和结构的记录：

Summary of Optimized Potential Surface Scan (add -158.0 to energies):
                   1        2        3        4        5
Eigenvalues -- -0.40297 -0.40237 -0.40105 -0.39920 -0.39725
……
                   6        7        8        9        10
Eigenvalues -- -0.39576 -0.39541 -0.39641 -0.39817 -0.40010
……
……

用GaussView打开输出文件，点击Results->Scan可以得到GaussView绘制的势能曲线，如下图所示：

柔性扫描势能曲线

一般在论文中不使用GaussView作的图，而是自行使用Origin等程序作图。柔性扫描不像刚性扫描一样在文件末尾以两列的形式给出扫描变量和能量，若要取出能量，可以借助GaussView，在上述图势能曲线处右击，选择Save Data，即可将扫描变量与能量写入txt文件中。

柔性扫描势能曲线与刚性扫描得到的势能曲线形状相同，但是由于进行过限制性优化总体能量更低，而且几个势能面上的极大值相对能量也有差别。将势能曲线中的最高点作为初猜做过渡态搜索，依然可以搜索到相应过渡态。读者可以自行将柔性扫描与刚性扫描的结果做比较。

二维势能面柔性扫描

Gaussian除了支持一维势能面扫描外，还支持多维的势能面扫描，这里以二维势能面扫描为例。二维势能面扫描的输入文件与一维势能面扫描的唯一不同就是文件尾有两行描述被扫描的两个内坐标。

[Type] N1 N2 [N3 [N4]] S steps step-size

这里我们称第一行描述的坐标为1st coordinate，第二行描述的坐标为2nd coordinate。

做一维势能面扫描时，程序只需要沿着势能面上的一个方向移动就可以了，做二维势能面扫描时就会出现移动方向的选择问题。Gaussian做二维势能面扫描时移动方向如下图所示：

二维势能面扫描示意图

从起点O开始，首先固定1st coordinate不变，依次改变2nd coordinate直至最大步长。在2nd coordinate的最大步长处，固定2nd coordinate，让1st coordinate移动一步。继续固定1st coordinate不变，依次改变2nd coordinate直至最小步长。在2nd coordinate的最小步长处，固定2nd coordinate，让1st coordinate移动一步。以此类推。因此二维势能面扫描产生的结构的总数为(steps_1+1)*(steps_2+1)，更多维势能面扫描产生的结构总数同理。在做多维势能面扫描时，变量之间最好不要有关联，否则会使扫描结果混乱或出错。

本文二维势能面柔性扫描所研究的体系为SN2反应，反应式如下。

Cl−+CH3Cl → ClCH3+Cl−

定义的两个内坐标分别为下图中两个C-Cl键的键长。扫描的起点时两个键长均为2.0埃，两个坐标的步数均为7，步长均为0.1埃。

输入文件如下：

%chk=sn2-scan.chk
%mem=8GB
%nprocshared=8
#p opt=modredundant nosymm b3lyp/6-31g(d,p)
 
Title Card Required
 
-1 1
 C    -0.00000100    0.00001000   -0.00009800
 H     0.00000500   -0.90861200    0.57116500
 H     0.00000500    0.94905800    0.50114900
 H     0.00000500   -0.04041300   -1.07262700
 Cl    1.99999900   -0.00000100    0.00000777
 Cl   -2.00000100   -0.00000100    0.00000777
 
 B 1 5 S 7 0.1
 B 1 6 S 7 0.1

计算结束后，在输出文件中同样有每一步的能量和结构信息：

Summary of Optimized Potential Surface Scan (add-960.0 to energies):
                   1        2        3        4        5
Eigenvalues -- -0.32917 -0.34305 -0.35429 -0.36323 -0.37009
……
……

用和上文打开一维势能面扫描输出文件时一样的操作，打开二维势能面扫描的输出文件时，GaussView会输出三维的势能图，如下所示。可以看到有鞍点存在。

二维势能面

推荐阅读：

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

2. http://gaussian.com/opt/

3. Exploring Chemistry with Electronic Structure Methods, Third Edition, Chapter 6, P238