一种用Gaussian 16中的GIC功能实现同时扫描多个坐标的方法

一、简介

势能面扫描是我们用Gaussian常做的计算，一般可以分为刚性扫描和柔性扫描。如果在柔性扫描中给定两个坐标，那么我们将会得到二维势能面。但是有时候我们只希望两个坐标同时变化得到一条势能曲线，这可以通过使用Gaussian中的GIC（广义内坐标）实现。本公众号之前也给出了一个可行的解决方案，见《在Gaussian16中同时扫描两个反应坐标》。但是之前方案的缺点是使用了Link1，在用GaussView打开输出文件时不能很方便地显示能量的变化趋势，这在找能量极大，极小点时会带来困难。因此这里给出了一个新的方案，不使用Link1，让势能曲线可以直观地显示出来。

这个新的方案同样需要使用GIC，没有接触过GIC的同学可以在Gaussian官网上学习相关的资料 (http://gaussian.com/gic)。GIC使Gaussian用户可以自定义一些结构参数，如键长、键角、二面角等。在基础的结构参数上还可以用数学运算，如加减乘除、平方开方、三角函数等，定义更复杂的结构参数。根据定义出的结构参数，我们可以做限制性优化，势能面扫描等等。

二、甲醛与水的加成反应

还是以甲醛与水的加成反应为例，介绍GIC输入文件的写法。

我们需要扫描的是O4−H7和C1−O5间的键长，因此需要对相关结构参数有一个定义：

RCO=R(1,5)

ROH=R(4,7)

RCO=R(1,5)定义了RCO为1号原子和5号原子之间的距离。ROH=R(4,7)定义了ROH为4号原子和7号原子之间的距离。注意新建变量名时，不建议用R12、R34这种变量名，因为可能会被高斯程序内部占用，引起冲突报错。与之前公众号推送的内容相同，我们希望RCO和ROH按照如下距离逐渐缩小：

RCO=1.8, 1.7, 1.6, 1.5, 1.4

ROH=1.7, 1.5, 1.3, 1.1, 0.9

正如前面提到的，我们不希望做二维势能面扫描，因此我们只能选择一个扫描坐标。这里选择RCO，我们只需要对RCO=R(1,5)做一些改动即可：

RCO(NSteps=4,StepSize=-0.1)=R(1,5)

NSteps=4指扫描步数为4，StepSize=−0.1指扫描步长为−0.1。我们先将初始结构调整为RCO=1.8，ROH=1.7，上面的命令可以使RCO按每一步0.1的步长减小，共减小5步，这正是我们所希望的。

在完成了RCO扫描的设置后，我们设置ROH，让ROH能随着RCO变化而变化。这可以用GIC中的Frozen功能完成。使用Frozen功能的前提是找出我们需要固定的量。我们将RCO和ROH的值输入到Excel中，并作出趋势线：

得到RCO和ROH间满足的关系为

ROH=2*RCO-1.9

这个等式说明2.0*RCO-ROH在整个扫描过程中可以作为一个不变量。

我们定义一个量F=2.0*RCO-ROH，并将其固定：

F(Frozen)=2.0*RCO-ROH

F(Frozen)可以让Gaussian在优化过程中固定F这个值不变。在扫描过程中，每当我们有一个新的RCO，由于F的限制，我们总会得到相应的ROH。这就让ROH可以随RCO变化而变化。

使用GIC功能的完整的输入文件如下：

%mem=24GB%nprocshared=24#p opt=(modredundant,GIC) M062X/6-31G(d) nosymm
Title Card Required
0 1 C       -1.12689502   -0.06090992    0.04839780 H       -0.59373127   -0.98861484    0.04839780 H       -2.19689502   -0.06090992    0.04839780 O       -0.49985459    1.03014042    0.04839780 O       -0.71648423   -0.27156123    1.78827969 H       -0.42609763   -1.05648810    2.25856198 H       -0.25552070    0.59195719    1.67268115
RCO(NSteps=4,StepSize=-0.1)=R(1,5)ROH=R(4,7)F(Frozen)=2.0*RCO-ROH

由于没有使用Link1，用GaussView打开输出文件不会显示多步任务。查看Scan的结果，可以很快地找到最高点，如下图所示：

将最高点作为初猜做过渡态搜索，即可找到甲醛和水加成的过渡态。另外，若把opt=(modredundant,GIC)替换成两个关键词opt geom=addGIC，效果一样。

简单总结一下，写同时扫描多个坐标所需Gaussian输入文件的通用步骤为：

1. 首先指定第一个扫描坐标，例如

RCO(NSteps=4,StepSize=-0.1)=R(1,5)

2. 定义下一个坐标，并根据新坐标和旧坐标的关系定义出一个不变量，如

ROH=R(4,7)

F(Frozen)=2.0*RCO-ROH

这一步可以通过函数的拟合来完成，相当于找到一个函数F(RCO,ROH)=0。

3. 重复2，直到所有坐标都定义完毕。

三、三个水分子间的质子转移反应

下面以三个水分子间的质子转移反应为例，演示如何同时扫描多个坐标。

我们首先对三个水分子的团簇做结构优化，得到稳定结构。

得到O−H间键长为0.97819，不成键的O−H原子间距离为1.83887。

为了得到3个质子同时转移的过渡态，需要同时缩短H2−O4，H6−O7，H9−O1间距离。扫描步数预设为10步，可得到步长约为−0.08。将H2−O4间距离定为扫描坐标，为了让H6−O7，H9−O1间距离随着H2−O4间距离减小而减小，需要引入两个不变量F1、F2：

F1(Frozen)=RH2O4-RH6O7

F2(Frozen)=RH2O4-RH9O1

根据前文的介绍，很容易写出GIC计算的输入文件：

%nprocshared=24%mem=24GB#p opt=(modredundant,GIC) m062x/6-31g(d,p) nosymm
Title Card Required
0 1 O      -2.12109900    0.15275700    0.36236700 H      -1.26182000   -0.15720700    0.01249000 H      -2.74456300    0.04217400   -0.36151100 O       0.46851300    0.12482800   -0.54235200 H       1.12745300   -0.25541100    0.04672800 H       0.37017800    1.05352900   -0.25131900 O      -0.56877500    2.41551400    0.54801000 H      -0.86226500    3.12772800   -0.02774600 H      -1.33253600    1.81105900    0.63231100
RH2O4(NSteps=10,StepSize=-0.08)=R(2,4)RH6O7=R(6,7)F1(Frozen)=RH2O4-RH6O7RH9O1=R(9,1)F2(Frozen)=RH2O4-RH9O1

用GaussiView打开输出文件，可以看到结构随扫描过程的变化及势能曲线。

取出能量最高点的结构做过渡态搜索，优化3步之后就能收敛到相应过渡态。可见柔性扫描可以给出相当好的过渡态初猜。

综上，GIC是一个Gaussian中很有用的功能，GIC结合柔性扫描对过渡态搜索有很大的帮助。