用Gaissian16中的GIC功能实现翻转过程的势能面扫描
用Gaissian16中的GIC功能实现翻转过程的势能面扫描
一、简介
势能面扫描在计算化学领域中有很多应用场景，本公众号之前也推送过相关介绍：
势能面扫描前需要用户对扫描坐标有一个明确的定义。在Gaussian 16的广义内坐标（GIC）功能出现之前，我们只能对一些简单的结构参数，如笛卡尔坐标、键长、键角、二面角做势能面扫描。GIC的出现让我们可以定义更加复杂的扫描坐标。
翻转过程在自由基和三取代的N、P体系中很常见。以氨气分子的翻转为例，三个H原子构成一个平面，翻转过程可以看成是N原子到此平面的距离不断变化的过程：
因此一个很自然的想法是：将N原子到三个H原子形成的平面的距离设置为扫描坐标，以此坐标做势能面扫描。实际操作中我并没有将距离作为扫描坐标，详细的扫描坐标的设置会在下文中解释。不使用GIC的话，这是很难实现的，使用GIC后问题就能迎刃而解。
二、相关数学公式的推导
为了得到扫描坐标在数学上的定义，我们需要对相关公式进行推导，先求出不共线三点的平面方程，再求出第四个点到这个平面的距离。
2.1 求过（不共线）三点的平面方程
已知三个点的坐标
得到两个向量
则垂直于此平面的法向量为
其中
则点法式平面方程可写为
将平面方程化为一般式
可求得
综上，不共线三点的平面方程为
其中
2.2 点到平面的距离
设空间中第四个点的坐标为
则此点到平面
的距离为
三、氨气分子翻转过程的势能面扫描
以氨气分子为例，介绍GIC在扫描翻转过程的势能面上的应用。
利用GIC功能的输入文件如下：
%nprocshared=24
%mem=24GB
#p opt=(modredundant,GIC) 6-31g(d,p) nosymm m062x
Title Card Required
0 1
H 1.80953700 -0.74056700 0.00001400
H 1.80955100 0.66571400 0.81195400
H 1.80951800 0.66574100 -0.81197200
N 1.41988000 0.19697800 0.00000500
X1=X(1)
Y1=Y(1)
Z1=Z(1)
X2=X(2)
Y2=Y(2)
Z2=Z(2)
X3=X(3)
Y3=Y(3)
Z3=Z(3)
X4=X(4)
Y4=Y(4)
Z4=Z(4)
NA=(Y2-Y1)*(Z3-Z1)-(Y3-Y1)*(Z2-Z1)
NB=(Z2-Z1)*(X3-X1)-(Z3-Z1)*(X2-X1)
NC=(X2-X1)*(Y3-Y1)-(X3-X1)*(Y2-Y1)
ND=-1.0*(NA*X1+NB*Y1+NC*Z1)
Dist(NSteps=20,StepSize=-0.07363)=(NA*X4+NB*Y4+NC*Z4+ND)/SQRT(NA**2+NB**2+NC**2)
其中
X1=X(1)
...
Z4=Z(4)
的作用是用自定义的变量X1，Y1，Z1等表示每个原子的笛卡尔坐标。
NA=(Y2-Y1)*(Z3-Z1)-(Y3-Y1)*(Z2-Z1)
NB=(Z2-Z1)*(X3-X1)-(Z3-Z1)*(X2-X1)
NC=(X2-X1)*(Y3-Y1)-(X3-X1)*(Y2-Y1)
ND=-1.0*(NA*X1+NB*Y1+NC*Z1)
的作用是计算三个H原子构成的平面方程。
Dist(NSteps=20,StepSize=-.07363)=(NA*X4+NB*Y4+NC*Z4+ND)/SQRT(NA**2+NB**2+NC**2)
的作用有两个，一是定义扫描坐标Dist，二是给定扫描的步数和步长。这是整个流程中最关键的部分。我们在2.2中计算的是点到平面的距离，是一个非负数。而这里Dist的表达式与2.2中有一个区别是去掉了绝对值符号，这使得Dist的值可以是负数。
易得按照新的Dist的定义，N原子在3个H原子构成的平面的两侧时，Dist有相反的符号。反之亦然，Dist改变大小和符号时，N原子就可以在3个H原子构成平面的两侧移动。
那么如何得到初始结构下Dist的大小和符号，让我们对扫描的步数和步长有一个定义呢？根据右手定则可以确定Dist的符号，不借助其他方法Dist的大小只能目测。一个比较简单的方法是用半经验方法算一个单点，Gaussian就会输出各GIC的值。
用半经验方法算单点的输入文件如下：
%nprocshared=24
%mem=24GB
#p nosymm pm6 geom=addGIC
Title Card Required
0 1
H 1.80953700 -0.74056700 0.00001400
H 1.80955100 0.66571400 0.81195400
H 1.80951800 0.66574100 -0.81197200
N 1.41988000 0.19697800 0.00000500
X1=X(1)
Y1=Y(1)
Z1=Z(1)
X2=X(2)
Y2=Y(2)
Z2=Z(2)
X3=X(3)
Y3=Y(3)
Z3=Z(3)
X4=X(4)
Y4=Y(4)
Z4=Z(4)
NA=(Y2-Y1)*(Z3-Z1)-(Y3-Y1)*(Z2-Z1)
NB=(Z2-Z1)*(X3-X1)-(Z3-Z1)*(X2-X1)
NC=(X2-X1)*(Y3-Y1)-(X3-X1)*(Y2-Y1)
ND=-1.0*(NA*X1+NB*Y1+NC*Z1)
Dist=(NA*X4+NB*Y4+NC*Z4+ND)/SQRT(NA**2+NB**2+NC**2)
注意，由于只需要做单点计算，RouteSection部分有改动。并且在Dist的定义部分去掉了势能面扫描的内容。
单点计算结束后，我们可以在输出文件中找到以下内容：
NA GIC-1 -8.1553
NB GIC-2 0.0
NC GIC-3 0.0002
ND GIC-4 27.8873
Dist GIC-5 0.7363可见Gaussian默认输出所有自定义的GIC。借助这个特性，我们也可以先随意给定NSteps和StepSize的值，提交Gaussian作业后马上kill掉进程，这样也可以在输出文件中找到初始结构下GIC的值。
在得到Dist数值之后，我们就可以设置扫描步数和步长为：
NSteps=20,StepSize=-0.07363
计算完毕后，用GaussView打开输出文件，可以看到氨气分子结构随扫描过程的变化及势能曲线。
取出势能曲线最高点的结构做过渡态搜索，优化两步就能收敛到翻转过渡态：
四、三取代膦翻转过程的势能面扫描
与三取代的胺的性质相似，三取代膦也可以有翻转过程。有文献指出，膦被氧化生成自由基正离子之后，翻转过程的能垒会降低。这里我们首先完成膦的翻转过程的势能面扫描，再做翻转过程的过渡态搜索。
中性膦的势能面扫描的输入文件如下：
%nprocshared=24
%mem=48GB
#p opt=(modredundant,GIC) m062x/6-31g(d,p) scrf=(smd,solvent=acetonitrile) nosymm
title
0 1
P -0.08670900 1.26287500 -0.95222100
C -1.47074400 0.32462200 -0.16515700
C -2.64121000 0.11500200 -0.89937700
C -1.39311100 -0.16786900 1.14315500
C -3.71993900 -0.56673800 -0.33663100
H -2.70680000 0.48251800 -1.92025200
C -2.46631700 -0.85288100 1.70502000
H -0.48207200 -0.02125600 1.71892100
C -3.63247400 -1.05223900 0.96516400
H -4.62347800 -0.72323100 -0.91717100
H -2.39458900 -1.23297200 2.71923500
H -4.46832100 -1.58825600 1.40318700
C 1.36920700 0.36444500 -0.26127900
C 1.61072200 -0.96464400 -0.67382300
C 2.27789400 0.98437200 0.60369100
C 2.74461700 -1.62383700 -0.19552200
C 3.40744700 0.31243600 1.06745000
H 2.11239800 2.00712700 0.92566900
C 3.64063200 -0.99858600 0.66843900
H 2.92760800 -2.64729400 -0.51196700
H 4.09752100 0.81484900 1.73762200
H 4.51618400 -1.53337000 1.02229700
C -0.13426100 2.78566600 0.10023900
H -1.07659600 3.29649900 -0.11053200
H -0.08645200 2.55787100 1.16860700
H 0.68341700 3.45840600 -0.16864500
C 0.66647000 -1.68997800 -1.59928300
H -0.25773600 -1.96742600 -1.08128400
H 0.37874200 -1.06348500 -2.45041700
H 1.12722100 -2.60172900 -1.98447600
X1=X(2)
Y1=Y(2)
Z1=Z(2)
X2=X(23)
Y2=Y(23)
Z2=Z(23)
X3=X(13)
Y3=Y(13)
Z3=Z(13)
X4=X(1)
Y4=Y(1)
Z4=Z(1)
NA=(Y2-Y1)*(Z3-Z1)-(Y3-Y1)*(Z2-Z1)
NB=(Z2-Z1)*(X3-X1)-(Z3-Z1)*(X2-X1)
NC=(X2-X1)*(Y3-Y1)-(X3-X1)*(Y2-Y1)
ND=-1.0*(NA*X1+NB*Y1+NC*Z1)
Dist(NSteps=20,StepSize=-0.16)=(NA*X4+NB*Y4+NC*Z4+ND)/SQRT(NA**2+NB**2+NC**2)
计算完毕后，用GaussView打开势能面扫描的输出文件，可以看到结构随扫描过程的变化及势能曲线。
取出势能曲线最高点的结构做过渡态搜索，优化三步就能收敛到翻转过渡态：
自由基正离子膦的势能面扫描的输入文件如下：
%nprocshared=24
%mem=48GB
#p opt=(modredundant,GIC) um062x/6-31g(d,p) scrf=(smd,solvent=acetonitrile) nosymm
title
1 2
P 0.11694400 1.13577500 0.56859500
C 1.56222700 0.22089900 0.07077600
C 2.82938400 0.69566800 0.44800900
C 1.43525800 -0.95614100 -0.68416900
C 3.96295700 -0.00178100 0.05363800
H 2.92984100 1.59576800 1.04687900
C 2.57899900 -1.64232500 -1.07070700
H 0.45518900 -1.31362900 -0.98503900
C 3.83790800 -1.16810400 -0.70181100
H 4.94362000 0.36112400 0.33926900
H 2.48792500 -2.54443000 -1.66494700
H 4.72693800 -1.71044700 -1.00496300
C -1.42137300 0.35948300 0.11779600
C -1.79142200 -0.86909900 0.71024200
C -2.26044300 1.00942000 -0.80274800
C -3.00255300 -1.43274500 0.31044200
C -3.45971400 0.42077300 -1.17677800
H -1.97302900 1.96156900 -1.23585100
C -3.82502300 -0.80361900 -0.62103400
H -3.31110200 -2.37555300 0.75066700
H -4.10351500 0.91537000 -1.89475800
H -4.76363900 -1.26757500 -0.90466800
C 0.21246100 2.87741700 0.09391000
H 1.10895300 3.31314400 0.53803100
H 0.27677500 2.95886300 -0.99456000
H -0.66737000 3.40215200 0.46823700
C -0.94138000 -1.55987300 1.74591200
H -0.08734000 -2.07014000 1.28988800
H -0.55120600 -0.85359000 2.48523900
H -1.52990900 -2.30908900 2.27677400
X1=X(2)
Y1=Y(2)
Z1=Z(2)
X2=X(23)
Y2=Y(23)
Z2=Z(23)
X3=X(13)
Y3=Y(13)
Z3=Z(13)
X4=X(1)
Y4=Y(1)
Z4=Z(1)
NA=(Y2-Y1)*(Z3-Z1)-(Y3-Y1)*(Z2-Z1)
NB=(Z2-Z1)*(X3-X1)-(Z3-Z1)*(X2-X1)
NC=(X2-X1)*(Y3-Y1)-(X3-X1)*(Y2-Y1)
ND=-1.0*(NA*X1+NB*Y1+NC*Z1)
Dist(NSteps=20,StepSize=-0.08964)=(NA*X4+NB*Y4+NC*Z4+ND)/SQRT(NA**2+NB**2+NC**2)
计算完毕后，用GaussView打开势能面扫描的输出文件，可以看到结构随扫描过程的变化及势能曲线：
取出势能曲线最高点的结构做过渡态搜索，优化三步就能收敛到翻转过渡态：
易知膦被氧化为自由基正离子后翻转过程变得更加容易了：
五、总结
使用GIC可以让我们快速实现复杂反应坐标的势能面扫描。面对翻转过程这一类扫描坐标很复杂的情况，如果不使用GIC，势能面扫描是相当难实现的。势能面扫描给出的初猜也可以让我们很快地找到翻转过程的过渡态。
参考文献
Kyle D. R., Daniel H. E. , and Alexander T. R.* J. Am. Chem. Soc. 2013, 135, 9354−9357