一、简介
势能面扫描是我们用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结合柔性扫描对过渡态搜索有很大的帮助。
