基于Gaussian的NMR=CSGT任务得到各个轨道对NICS贡献的方法
基于Gaussian的NMR=CSGT任务得到各个轨道对NICS贡献的方法
The method to obtain the contribution of each orbital to NICS based on the Gaussian NMR=CSGT task
文/Sobereva@北京科音 2023-Jun-3
1 前言
NICS是最流行的考察芳香性的方法之一，《衡量芳香性的方法以及在Multiwfn中的计算》（http://sobereva.com/176）中我做过介绍，如果不会算的话先仔细看此文。对于某些电子结构特殊的体系，比如18碳环（http://sobereva.com/carbon_ring.html）等双芳香性体系，将NICS分解为不同类型轨道进行研究对于考察芳香性的本质特征是很有好处的。做这种分解的方法之一是利用NBO的NCS分析功能，但利用它将NICS分解为分子轨道的贡献需要花钱买NBO，而Gaussian自带的NBO 3.1没法实现。不花钱的方法是借助Gaussian的AICD程序的界面，之前我在《使用AICD 2.0绘制磁感应电流图》（http://sobereva.com/294）里简单提过一句，在本文将更具体说明。下文会首先举例对苯计算pi和其它轨道对NICS产生的贡献，然后介绍我写的一个可以自动计算每个轨道对NICS的贡献的脚本的使用。本文使用G16 C.02进行计算，读者使用的版本必须是G09 D.01及之后。
注意，由于NCS是基于GIAO做的，本文的方法是基于CSGT做的，因此即便计算所有电子贡献的总的NICS也会存在差异，但一般差异很小，可忽略不计。这两种方法产生的每个轨道对NICS的贡献可能相差非常大，但每类轨道贡献的总和差异通常很小，例如苯的所有pi轨道的贡献。
下文涉及的所有文件都可以在http://sobereva.com/attach/670/file.zip里找到。
本文的这种NICS分解方法在Multiwfn绘制NICS曲线图中也支持，由此可以考察不同轨道对特定路径上磁屏蔽的贡献，见《使用Multiwfn绘制一维NICS曲线并通过积分衡量芳香性》（http://sobereva.com/681）里的介绍。
《详谈分子轨道对磁感生电流的贡献》（http://sobereva.com/703）是非常重要的一篇博文，一定要仔细阅读，否则可能会对轨道对NICS的贡献的意义有错误的认识。
2 计算苯的pi电子和其它电子对NICS_ZZ的贡献
首先我们基于CSGT计算一下苯的NICS(1)ZZ，顺带观看分子轨道。输入文件如下，对应本文文件包里的benzene_NICS.gjf，结构是在B3LYP/6-31G*下优化过的
%chk=benzene.chk
# b3lyp/6-31g(d) NMR=CSGT
[空行]
b3lyp/6-31G* opted
[空行]
0 1
C 0.00000000 1.39661199 0.00000000
C 1.20950146 0.69830599 0.00000000
C 1.20950146 -0.69830600 0.00000000
C -0.00000000 -1.39661199 0.00000000
C -1.20950146 -0.69830599 0.00000000
C -1.20950146 0.69830600 0.00000000
H 0.00000000 2.48362808 0.00000000
H 2.15088501 1.24181404 0.00000000
H 2.15088501 -1.24181404 0.00000000
H -0.00000000 -2.48362808 0.00000000
H -2.15088501 -1.24181404 0.00000000
H -2.15088501 1.24181404 0.00000000
Bq 0. 0. 1.0
[空行]
[空行]
计算完毕后，看到NICS(1)ZZ是-27.3147：
13 Bq Isotropic = 13.6602 Anisotropy = 20.4819
XX= 6.8585 YX= 0.0000 ZX= 0.0000
XY= -0.0000 YY= 6.8072 ZY= 0.0000
XZ= 0.0000 YZ= -0.0000 ZZ= 27.3147
将benzene.chk转成fch文件，然后用Multiwfn按照http://sobereva.com/269观看轨道，或者让Multiwfn按照《在Multiwfn中单独考察pi电子结构特征》（http://sobereva.com/432）自动识别pi轨道，可知17、20、21号轨道是pi轨道。下面就计算这三个pi轨道对NICS的贡献。创建如下输入文件，对应本文文件包里的benzene_NICS-pi.gjf。其中guess=read读取之前的chk文件是为了省SCF过程的时间。IOp(10/93=2)必须写，这样才能利用Gaussian的AICD界面对考虑的轨道进行指定。此界面产生的给AICD程序用的文件跟我们当前没关系，故产生的文件随便起个名字叫nouse.txt。
%chk=benzene.chk
# b3lyp/6-31g(d) NMR=CSGT IOp(10/93=2) guess=read
[空行]
b3lyp/6-31G* opted
[空行]
0 1
C 0.00000000 1.39661199 0.00000000
C 1.20950146 0.69830599 0.00000000
C 1.20950146 -0.69830600 0.00000000
C -0.00000000 -1.39661199 0.00000000
C -1.20950146 -0.69830599 0.00000000
C -1.20950146 0.69830600 0.00000000
H 0.00000000 2.48362808 0.00000000
H 2.15088501 1.24181404 0.00000000
H 2.15088501 -1.24181404 0.00000000
H -0.00000000 -2.48362808 0.00000000
H -2.15088501 -1.24181404 0.00000000
H -2.15088501 1.24181404 0.00000000
Bq 0. 0. 1.0
[空行]
nouse.txt
[空行]
17 20 21
[空行]
[空行]
计算结果如下，即曰pi电子对NICS(1)ZZ产生的贡献为-29.7306
13 Bq Isotropic = 14.7855 Anisotropy = 22.4176
XX= 7.3287 YX= -0.0000 ZX= 0.0000
XY= -0.0000 YY= 7.2973 ZY= -0.0000
XZ= -0.0000 YZ= -0.0000 ZZ= 29.7306
其它电子对NICS(1)ZZ产生的贡献可以用NICS(1)ZZ总值减去pi电子的贡献，即-27.3147-(-29.7306) = 2.4159。如果你想直接算也可以，结果是一样的，即把上面输入文件里的17 20 21改为其它占据轨道的序号1-16 18 19，相应的输入文件是本文文件包里的benzene_NICS-sigma.gjf，结果如下，可见其它电子对NICS(1)ZZ的贡献确实为2.4158
13 Bq Isotropic = -1.1254 Anisotropy = 0.9827
XX= -0.4702 YX= 0.0000 ZX= 0.0000
XY= -0.0000 YY= -0.4900 ZY= 0.0000
XZ= 0.0000 YZ= -0.0000 ZZ= -2.4158
顺带提醒一下，对于非限制性开壳层体系，Gaussian的AICD界面在指定轨道序号的时候的规则比较特殊，注意看《使用AICD 2.0绘制磁感应电流图》（http://sobereva.com/294）里面的说明，千万别搞错了。
3 计算各个占据轨道对NICS的贡献
看完上一节，读者自然知道怎么计算各个分子轨道对NICS的贡献了，设好序号后计算即可。但如果你想对几十甚至上百个轨道产生的贡献都一一计算，不借助脚本而人力操作就太像原始人了。我特意写了个Bash shell脚本自动调用Gaussian干这个事，是本文文件包里的orbNICS.sh。用法是：
(1) 对当前体系做个单点任务，令所得的chk文件叫做NICS.chk
(2) 创建模板文件NICS.gjf，内容是按照上一节的方法计算特定轨道对NICS贡献的输入文件，必须有且只有一个Bq，轨道序号那行写为orbidx，会自动被脚本替换为要算的轨道序号。并且用guess=read要求从%oldchk指定的NICS.chk中读取初猜波函数
(3) 修改orbNICS.sh，令第6行的循环的起始和终止序号对应要考察的轨道编号范围。并且恰当设置里面的Gaussian运行命令，默认是g16。
(4) 运行orbNICS.sh，脚本会调用Gaussian依次计算各个轨道对NICS的贡献，对Bq位置的各向同性磁屏蔽值和ZZ分量的贡献会分别写入到当前目录下的iso.txt和ZZ.txt中
下面就看具体计算例子，计算苯的所有占据轨道（1-21号）的贡献。首先运行以下任务产生NICS.chk，此输入文件是本文文件包里的SP.gjf。
%chk=NICS.chk
#p b3lyp/6-31G*
[空行]
test
[空行]
0 1
C 0.00000000 1.39661199 0.00000000
C 1.20950146 0.69830599 0.00000000
C 1.20950146 -0.69830600 0.00000000
C -0.00000000 -1.39661199 0.00000000
C -1.20950146 -0.69830599 0.00000000
C -1.20950146 0.69830600 0.00000000
H 0.00000000 2.48362808 0.00000000
H 2.15088501 1.24181404 0.00000000
H 2.15088501 -1.24181404 0.00000000
H -0.00000000 -2.48362808 0.00000000
H -2.15088501 -1.24181404 0.00000000
H -2.15088501 1.24181404 0.00000000
[空行]
[空行]
创建以下模板文件，是本文文件包里的NICS.gjf。
%oldchk=NICS.chk
# b3lyp/6-31G* nmr=csgt IOp(10/93=2) guess=read
[空行]
test
[空行]
0 1
C 0.00000000 1.39661199 0.00000000
C 1.20950146 0.69830599 0.00000000
C 1.20950146 -0.69830600 0.00000000
C -0.00000000 -1.39661199 0.00000000
C -1.20950146 -0.69830599 0.00000000
C -1.20950146 0.69830600 0.00000000
H 0.00000000 2.48362808 0.00000000
H 2.15088501 1.24181404 0.00000000
H 2.15088501 -1.24181404 0.00000000
H -0.00000000 -2.48362808 0.00000000
H -2.15088501 -1.24181404 0.00000000
H -2.15088501 1.24181404 0.00000000
Bq 0. 0. 1.
[空行]
nouse.txt
[空行]
orbidx
[空行]
[空行]
把orbNICS.sh里轨道序号循环范围设成1到21。令NICS.chk、NICS.gjf和orbNICS.sh都在当前目录，运行chmod +x orbNICS.sh增加可执行权限，然后运行orbNICS.sh，之后会看到
Calculating NICS contributed by orbital 1...
Calculating NICS contributed by orbital 2...
Calculating NICS contributed by orbital 3...
...略
Calculating NICS contributed by orbital 19...
Calculating NICS contributed by orbital 20...
Calculating NICS contributed by orbital 21...
算完后，当前目录下的iso.tx如下所示。第1列是轨道序号，第2列是各个轨道对Bq位置的isotropic磁屏蔽值的贡献，取负值就是对NICS的贡献。末尾输出的Total是所有轨道数值的加和，其13.6602和第2节直接用NMR=CSGT算出来的总磁屏蔽完全一致。
1 -0.0061
2 -0.2657
3 -0.2657
...略
17 2.8926
18 -2.7037
19 -2.7032
20 5.9459
21 5.9470
Total: 13.6602
当前目录下的ZZ.txt是各个轨道对Bq位置磁屏蔽张量ZZ分量的贡献，取负值就是对NICS_ZZ的贡献：
1 0.0162
2 -0.1581
3 -0.1581
...略
17 1.5800
18 -5.4334
19 -5.4334
20 14.0752
21 14.0754
Total: 27.3147
大家可以把NICS.chk转成fch后用Multiwfn一边看轨道图形一边照着以上列出的贡献值进行分析。
当前目录下还有NICS1.out、NICS2.out ... NICS21.out，是计算各个轨道贡献时留下的Gaussian输出文件，可以删掉。每次脚本启动时都会自动删当前目录下所有的.out文件和iso.txt、ZZ.txt。