基于Gaussian的NMR=CSGT任務得到各個軌道對NICS貢獻的方法
基于Gaussian的NMR=CSGT任務得到各個軌道對NICS貢獻的方法
文/Sobereva@北京科音 2023-Jun-3
1 前言
NICS是最流行的考察芳香性的方法之一，《衡量芳香性的方法以及在Multiwfn中的計算》（http://www.shanxitv.org/176）中我做過介紹，如果不會算的話先仔細看此文。對于某些電子結構特殊的體系，比如18碳環（http://www.shanxitv.org/carbon_ring.html）等雙芳香性體系，將NICS分解為不同類型軌道進行研究對于考察芳香性的本質特征是很有好處的。做這種分解的方法之一是利用NBO的NCS分析功能，但利用它將NICS分解為分子軌道的貢獻需要花錢買NBO，而Gaussian自帶的NBO 3.1沒法實現。不花錢的方法是借助Gaussian的AICD程序的界面，之前我在《使用AICD 2.0繪制磁感應電流圖》（http://www.shanxitv.org/294）里簡單提過一句，在本文將更具體說明。下文會首先舉例對苯計算pi和其它軌道對NICS產生的貢獻，然后介紹我寫的一個可以自動計算每個軌道對NICS的貢獻的腳本的使用。本文使用G16 C.02進行計算，讀者使用的版本必須是G09 D.01及之后。
注意，由于NCS是基于GIAO做的，本文的方法是基于CSGT做的，因此即便計算所有電子貢獻的總的NICS也會存在差異，但一般差異很小，可忽略不計。這兩種方法產生的每個軌道對NICS的貢獻可能相差非常大，但每類軌道貢獻的總和差異通常很小，例如苯的所有pi軌道的貢獻。
下文涉及的所有文件都可以在http://www.shanxitv.org/attach/670/file.zip里找到。
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://www.shanxitv.org/269觀看軌道，或者讓Multiwfn按照《在Multiwfn中單獨考察pi電子結構特征》（http://www.shanxitv.org/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://www.shanxitv.org/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。