在Gaussian16中同时扫描两个反应坐标
本公众号之前推送过在高斯中的两种常见势能面扫描：
用高斯做势能面扫描（一）：刚性扫描
用高斯做势能面扫描（二）：柔性扫描
可能大家都熟知，在柔性扫描中如果写了两个扫描坐标，如
B 1 5 S 7 0.1
B 1 6 S 7 0.1
是依次扫描两个坐标，无法做到同时，因此得到的是一张二维势能面，总扫描点数是两个坐标扫描点数的乘积，计算量较大。然而有时候我们只想同时扫描两个反应坐标，即两个坐标同时改变，得到一条曲线。
例如，找[2+2]环加成反应的过渡态经常会碰到这种问题，对于复杂的分子结构，手动调整过渡态初猜很难合适，此时使用opt=ts找到过渡态的成功率自然也不高，这时候我们可能就想，取柔性扫描势能曲线（面）上的突跃点作为过渡态初猜，然而去扫两个坐标得到一张二维势能面未免过于耗时。
对于这种问题，笔者以往采用了两种做法：（1）写了一个小程序来产生调整键长后新的结构（不仅仅是拉进/远两个原子），然而产生的结构我并不满意，算法还需改进或者仍有bug。（2）若仅算一两步反应，那么就手动在GaussView里调整好两个键长，每次算完下载下来再调键长，这样扫描5个点就要下载、调整5次，甚是麻烦。
由于G16推出了广义内坐标（GIC）功能，于是笔者便研究了一下官网的说明
http://gaussian.com/gic/
琢磨出了一个同时扫描两个键长的文件模板。可能（大概率）这不是最优的方案，欢迎有更简洁做法的小伙伴们留言分享你的技巧，我们可以日后再出一期更新推送。
为简洁起见，本文以甲醛和水的加成反应（真空中）为例，示范如何写输入文件。
在这个反应中有两个主反应坐标（C−O键和O−H键）同时在动，单独去扫描C−O键或者O−H键能量都会一直升高，并不会有突跃点。为防止混淆，短横线−符号左边始终表示甲醛分子中的原子，符号右边则表示水分子中的原子。
当然，化学直觉较强的同学可以直接构造出这个简单反应的合理过渡态初始结构，但这招对复杂结构很难见效，因此有必要掌握同时扫描两个反应坐标的技巧。假设我们希望C−O键和O−H键按如下距离逐渐减小
C−O: 1.8, 1.7, 1.6, 1.5
O−H: 1.7, 1.5, 1.3, 1.1
则利用GIC功能的输入文件如下：
%chk=methanal_water.chk
%mem=24GB
%nprocshared=24
#p M062X/6-31G(d) nosymm geom=addgic
Title Card Required
0 1
C -1.13256115 -0.04504504 0.00000000
H -0.59939740 -0.97274996 0.00000000
H -2.20256115 -0.04504504 0.00000000
O -0.50552072 1.04600530 0.00000000
O -0.71081810 -0.28742611 1.83667749
H -0.42043150 -1.07235298 2.30695978
H -0.23246558 0.54490727 1.83667749
R15=R(1,5)
R47=R(4,7)
R15(value=1.8)
R47(value=1.7)
--Link1--
%chk=methanal_water.chk
%mem=24GB
%nprocshared=24
#p opt M062X/6-31G(d) nosymm geom(allcheck,addGIC) guess=read
R15(freeze)
R47(freeze)
--Link1--
%chk=methanal_water.chk
%mem=24GB
%nprocshared=24
#p M062X/6-31G(d) nosymm geom(allcheck,addGIC) guess=read
R15(value=1.7)
R47(value=1.5)
--Link1--
%chk=methanal_water.chk
%mem=24GB
%nprocshared=24
#p opt M062X/6-31G(d) nosymm geom(allcheck,addGIC) guess=read
R15(freeze)
R47(freeze)
--Link1--
%chk=methanal_water.chk
%mem=24GB
%nprocshared=24
#p M062X/6-31G(d) nosymm geom(allcheck,addGIC) guess=read
R15(value=1.6)
R47(value=1.3)
--Link1--
%chk=methanal_water.chk
%mem=24GB
%nprocshared=24
#p opt M062X/6-31G(d) nosymm geom(allcheck,addGIC) guess=read
R15(freeze)
R47(freeze)
--Link1--
%chk=methanal_water.chk
%mem=24GB
%nprocshared=24
#p M062X/6-31G(d) nosymm geom(allcheck,addGIC) guess=read
R15(value=1.5)
R47(value=1.1)
--Link1--
%chk=methanal_water.chk
%mem=24GB
%nprocshared=24
#p opt M062X/6-31G(d) nosymm geom(allcheck,addGIC) guess=read
R15(freeze)
R47(freeze)
其中R15和R47是新起的变量名，可以改成其他名称。而R(1,5)则表示1号和5号原子之间的距离。注意value并不能限制键长严格等于给定的值，只能是十分接近，这是高斯官网提到的（不过这点对普通用户并不重要）。
第一步是先用GIC调整两根化学键至目标键长，调整完成后会做单点计算，下一个任务--Link1--读取之，将调整好的键长冻结、做限制性优化；完成后再用GIC调整两根键至新的键长，算单点，再--Link1--读取之，将调整好的键长冻结、做限制性优化。。。如此几步下来就能同时让两根键缩短了。从始至终只用了一个gjf文件，中途无需人为干预，可以认为是达到了目的。完成后用GV 6.0打开log文件
会发现有很多任务，这是用了--Link1--造成的，只能打开一个个Optimization查看结构和能量，看在哪个任务处发生了突跃。这点不是很方便，后续如果没有发现更优的方案，笔者会考虑写个脚本从log里提取信息，避免这一步操作。
回到这个问题，经过查看后发现第三个任务是能量最高点，下一步能量就会下降了。将第三个任务中最后一帧结构另存为、提取出来作为过渡态初猜，经过10步可以收敛到正确的过渡态上
总结：本文用一个简单的反应展示了如何同时扫描两个反应坐标，该反应若仅扫描其中任何一个反应坐标都是得不到突跃点的。对于更密的扫描步长、及扫描键角等等，读者可根据文中提供的示例文件自己举一反三。
PS1：作为练习，感兴趣的读者可以自己试试找溴化氢HBr与乙烯或乙炔加成的过渡态，考验一下自己的化学直觉。若找不到的话，可试试本文的技巧。
PS2: GIC功能仅在>= G16 A版本才有，G09无法使用该功能。