將核獨立化學位移(NICS)分解為sigma和pi軌道的貢獻
2023-Jun-3注：筆者后來又寫了《基于Gaussian的NMR=CSGT任務得到各個軌道對NICS貢獻的方法》（http://www.shanxitv.org/670）介紹另一種分解NICS的方法，可以得到各個分子軌道的貢獻，只需要用Gaussian而不依賴于NBO
將核獨立化學位移(NICS)分解為sigma和pi軌道的貢獻
文/Sobereva @北京科音
First release: 2012-May-31 Last update: 2014-Aug-1
1 前言
衡量芳香性的指標很多，也沒有統一的定義。現今存在幾十種指標，從磁性質、電子離域性、能量、結構等諸多方面對芳香性進行描述。NICS是被使用得最廣泛的衡量芳香性的指標。它的含義是在某個人為設定的不在原子核位置上的磁屏蔽值的負值，越負（對磁場屏蔽越強）則芳香性越強。這個位置設定得有一定任意性。原文(JACS,118,6317)是取在共軛環的幾何中心，為了明確起見，后來被特稱為NICS(0)。有人認為取在環中心會將sigma和pi軌道的貢獻疊加在一起而說不清楚，于是提出取在平面上方或下方1埃的位置，稱為NICS(1)，這個指標體現的主要是pi電子的貢獻。后來有人更進一步認為，NICS值不該取各項同性值（磁屏蔽張量XX、YY、ZZ之和的平均），而應該只考察ZZ值（Z垂直于環平面），這樣最適合體現pi芳香性，這被稱為NICS(1)_ZZ。雖然NICS(1)_ZZ很適合討論pi芳香性，但是還是稍微摻入了一點（盡管可忽略）的sigma軌道的影響，而且也不可能通過選取某個空間位置的來單獨討論sigma芳香性。
磁屏蔽值可以是分解為軌道的貢獻的，因此，直接研究各個sigma軌道和各個pi軌道對NICS值的貢獻就可以直接研究sigma和pi芳香性。然而一般量子化學程序都不支持分解為軌道的貢獻，雖然一些文獻中使用DeMon-Master程序去做，但是此程序貌似是私下傳播，難以獲取（注意這并不是公開的DeMon2k，DeMon2k是做不到的）。相對明朗一些的方法是使用NBO程序的NCS（自然化學屏蔽）功能來做，它既可以將磁屏蔽張量分解為定域化軌道的貢獻，也可以分解為正則分子軌道(CMO)的貢獻。
實際上，NBO程序3.1版及后續版本的NCS功能有很大不同。NBO 3.1能做的是將磁屏蔽張量分解為NLMO的貢獻，而不能分解為NBO或MO的貢獻。而從后續版本，比如常用的5.0版，能做的是將磁屏蔽張量分解為NBO或MO的貢獻，據我所知，反倒不能再分解為NLMO的貢獻了。
能做NCS的NBO程序必須是內嵌到Gaussian作為l607模塊的版本，獨立運行的NBO不行。在第二節，將介紹如何用NBO5.0版l607做NCS，這也被稱為NBO5.G。NBO5.G是收費的，但網上能找到源碼，可以自行編譯，方法見《編譯NBO5.0獨立運行版和嵌入Gaussian03 C02版的方法》（http://www.shanxitv.org/144）。在第三節將介紹如何直接用Gaussian自帶的NBO3.1版的L607做NCS。
2 使用NBO 5.0分解磁屏蔽張量
NBO5.0和G03的C及以后版本之間有點兼容性問題。也就是，直接通過Bq虛原子定義NICS計算的位置的話，NBO5.0發現Bq原子沒有S基函數，就會報錯終止。為解決這個問題，要么用G03的A或B版本，要么用更新版本的NBO5.x程序。雖說可以把Bq寫成比如H-Bq，這樣的話這個虛原子上就有了當前基組中H的基函數，可以避免報錯，但這額外的基函數會明顯影響NICS值。更好的辦法是給這個虛原子只賦予一個彌散范圍很大的S基函數，這個基函數基本不會對體系能量、NICS值帶來什么影響，卻能夠由此繞過這個不兼容問題。
本文都用B3LYP/6-31+G*優化的苯作為例子，NICS、NCS計算也是在這個級別下進行。輸入文件如下，Z方向垂直于苯環平面
#P b3lyp/gen NMR pop=nboread //雖然相關資料說還需要寫IOp(10/46=1)，但筆者發現沒必要
b3lyp/6-31+g(d) opted
0 1
C 0.00000000 1.39864200 0.00000000
C 1.21125900 0.69932100 0.00000000
C 1.21125900 -0.69932100 0.00000000
C 0.00000000 -1.39864200 0.00000000
C -1.21125900 -0.69932100 0.00000000
C -1.21125900 0.69932100 0.00000000
H 0.00000000 2.48606800 0.00000000
H 2.15299800 1.24303400 0.00000000
H 2.15299800 -1.24303400 0.00000000
H 0.00000000 -2.48606800 0.00000000
H -2.15299800 -1.24303400 0.00000000
H -2.15299800 1.24303400 0.00000000
bq 0. 0. 0. //環中心
bq 0. 0. 1. //環平面中心上方1埃
C H 0
6-31+g(d)
****
13 14 0 //設定13、14號原子，也就是虛原子的基函數
S 1 1.0 //都只含有一個收縮度為1的S型基函數
0.005 1. //第一個數值越小，表明這個基函數越彌散，它的存在對結果影響也會越小，0.005已經是十分彌散了。如果設得太小會導致報錯。第二個數是收縮系數，對于當前基函數收縮度為1的情況這個值是隨意的，一般就設1.0。
****
$NBO NCS <XYZ MO> $END // $NBO和$END夾著的是傳遞給NBO模塊的關鍵詞。被<>括住的關鍵詞是NCS分析的具體選項
Gaussian輸出的兩個虛原子的磁屏蔽張量信息為
13 Bq Isotropic = 7.8997 Anisotropy = 7.6878
XX= 5.3444 YX= 0.0000 ZX= 0.0000
XY= 0.0000 YY= 5.3298 ZY= 0.0000
XZ= 0.0000 YZ= 0.0000 ZZ= 13.0249
Eigenvalues: 5.3298 5.3444 13.0249
14 Bq Isotropic = 10.0608 Anisotropy = 27.9171
XX= 0.7644 YX= 0.0000 ZX= 0.0000
XY= 0.0000 YY= 0.7457 ZY= 0.0000
XZ= 0.0000 YZ= 0.0000 ZZ= 28.6722
Eigenvalues: 0.7457 0.7644 28.6722
這表明NICS(0)為-7.8997，NICS(0)_ZZ=-13.0249。NICS(1)為-10.0608，NICS(1)_ZZ=-28.6722。可見，無論是在環中心還是其上方1埃處，環上離域電子對垂直于環平面的磁場屏蔽強度比起平行于環平面的要明顯大得多。
Isotropic（各向同性屏蔽值）就是一般所說的絕對化學位移，單位是ppm，和溶液中的普通核磁共振實驗數據是對應的（相對于混亂運動的分子，外磁場相當于沒有特定方向，故視為各向同性），這個值是磁屏蔽張量矩陣的跡的平均，比如7.8997=(5.3444+5.3298+13.0249)/3。Anisotropy用來衡量各項異性的程度，計算方法是Eig(3) - (Eig(2) + Eig(1))*0.5，Eig(i)代表第i個本征值，從小到大排序。如果屏蔽張量在各方向相同，此值即為0。
為了具體了解各個NBO及其離域效應對屏蔽值的影響，需要在NCS關鍵詞后的<>內寫上XYZ關鍵詞。這樣就會依次輸出每個原子的磁屏蔽張量的來自各個NBO的貢獻。比如對于環中心的虛原子：
NBO
L NL XX YX ZX XY YY ZY XZ YZ ZZ
===============================================================================
1. -4.16 2.65 0.00 2.65 -1.10 0.00 0.00 0.00 -3.73
50. -0.95 0.30 0.00 -0.43 0.14 0.00 0.00 0.00 -1.52
51. 0.66 -0.22 0.00 0.12 -0.17 0.00 0.00 0.00 0.88
54. 0.33 -0.14 0.00 0.40 -0.21 0.00 0.00 0.00 0.23
55. 0.15 -0.05 0.00 0.03 -0.03 0.00 0.00 0.00 0.13
56. 0.39 -0.21 0.00 0.57 -0.32 0.00 0.00 0.00 0.15
......略
21. -0.08 -0.51 0.00 -0.51 -0.67 0.00 0.00 0.00 -1.76
24. 0.00 0.00 0.00 0.02 0.04 0.00 0.00 0.00 0.15
28. -0.01 -0.03 0.00 0.03 0.05 0.00 0.00 0.00 0.11
76. 0.02 0.04 0.00 0.01 0.02 0.00 0.00 0.00 0.15
80. 0.03 0.05 0.00 0.00 0.00 0.00 0.00 0.00 0.11
-------------------------------------------------------------------------------
L -3.84 -0.02 0.00 -0.01 -3.85 0.00 0.00 0.00 -22.66
NL 9.25 -0.02 0.00 -0.03 9.24 0.00 0.00 0.00 35.69
-------------------------------------------------------------------------------
Total 5.41 -0.04 0.00 -0.04 5.40 0.00 0.00 0.00 13.03
-------------------------------------------------------------------------------
第一列是NBO的編號，如果靠左邊寫就說明它是占據數較高的NBO軌道，也被稱為Lewis型(L)。靠右寫的是非Lewis型(NL)，絕大多數都是占據數很低的NBO，但也可以是有一定占據但不那么接近2的NBO。比如XX標簽下面的第三個數0.66，可以解釋為1號NBO向51號NBO離域的微量電子對屏蔽張量XX分量的影響為0.66ppm（更確切來說，是51號NBO對應的NLMO由1號NBO引入的“離域性尾巴”的貢獻）。最后Total顯示的總XX/YY/ZZ值5.41/5.40/13.03和Gaussian前面給出的XX/YY/ZZ方向屏蔽值5.3444/5.3298/13.0249是十分吻合的。Total當中L就是所有L項的加和，NL就是所有NL項的加和。注意只有數值大于特定閾值的項才列在表里，閾值可以通過比如NCS=0.02關鍵詞來設定為0.02。仔細觀看ZZ列，并結合軌道圖形觀察，會發現三個C-C pi型NBO對磁場有很強屏蔽作用，每個貢獻都高達9.88ppm。而C-C和C-H的sigma鍵都會產生一定的去屏蔽作用，每個鍵都貢獻負三點幾。
緊接著輸出的是這個原子屏蔽張量主軸方向的值，33對應整體屏蔽磁場最強的方向，11對應最弱的方向。（如果表中有任何一項的12、13或23方向的值不很小，那么非對角元也會被輸出出來，不過此例所有項的非對角元都很小所以沒顯示）。主軸與笛卡爾坐標的變換矩陣也會輸出。ISO就是各向同性屏蔽值，CSA就是各向異性值，和前面給出的計算式子一樣，比如第一行，-1.09=-3.73-0.5*(0.24-5.50)。由于列數比較少，所以每行開頭都留出了空間直接給出了NBO編號對應的含義。如NL 54. C 3(ry*)就代表這是54號NBO，是3號碳上的里德堡軌道。
NBO 11 22 33 CSA ISO
============================================================
1. C 1- C 2 0.24 -5.50 -3.73 -1.09 -3.00
NL 50. C 3(ry*) -0.37 -0.43 -1.52 -1.11 -0.77
NL 51. C 3(ry*) 0.12 0.36 0.88 0.64 0.46
NL 54. C 3(ry*) 0.14 -0.02 0.23 0.17 0.11
NL 55. C 3(ry*) 0.03 0.09 0.13 0.07 0.08
NL 56. C 3(ry*) 0.15 -0.08 0.15 0.11 0.07
......略
20. C 5(cr) 0.08 -0.83 -1.76 -1.38 -0.84
NL 63. C 4(ry*) 0.01 0.03 0.15 0.12 0.06
NL 67. C 4(ry*) 0.02 0.01 0.11 0.10 0.05
NL 89. C 6(ry*) 0.00 0.05 0.15 0.12 0.06
NL 93. C 6(ry*) -0.01 0.04 0.11 0.10 0.05
21. C 6(cr) -0.94 0.18 -1.76 -1.38 -0.84
NL 24. C 1(ry*) 0.04 0.01 0.15 0.12 0.06
NL 28. C 1(ry*) 0.02 0.01 0.11 0.10 0.05
NL 76. C 5(ry*) 0.05 0.00 0.15 0.12 0.06
NL 80. C 5(ry*) 0.04 -0.01 0.11 0.10 0.05
------------------------------------------------------------
Lewis -3.86 -3.83 -22.66 -18.82 -10.12
non-Lewis 9.23 9.27 35.69 26.44 18.06
------------------------------------------------------------
Total 5.37 5.44 13.03 7.62 7.95
------------------------------------------------------------
最后NCS分析會輸出以下匯總內容，如果不寫<XYZ>的話NCS分析也只會輸出這些而不會輸出上述每個原子的成分細節
Summary of isotropic NMR chemical shielding
Total Lewis (L) and non-Lewis (NL) contributions: (ppm)
NBO C( 1) C( 2) C( 3) C( 4) C( 5) C( 6) H( 7)
------------ ------- ------- ------- ------- ------- ------- -------
1. C 1- C 2 L -49.21 -49.22 -4.75 -0.06 -0.06 -4.75 -1.67
NL -2.91 -2.90 3.56 -0.07 -0.08 3.56 -0.26
2. C 1- C 2 L 0.27 0.28 3.24 1.61 1.61 3.25 0.58
NL 1.47 1.47 -3.04 2.87 2.88 -3.07 0.18
3. C 1- C 6 L -49.21 -4.75 -0.06 -0.06 -4.75 -49.22 -1.67
NL -2.91 3.56 -0.08 -0.07 3.56 -2.90 -0.26
4. C 1- H 7 L -35.94 -4.61 -0.82 -0.45 -0.82 -4.61 28.45
NL 7.33 2.22 0.49 0.19 0.49 2.22 -1.23
......略
20. C 5(cr) L -0.26 -0.21 -0.26 -0.41 203.56 -0.41 -0.10
NL 0.08 0.05 0.08 0.05 0.11 0.05 0.02
21. C 6(cr) L -0.41 -0.26 -0.21 -0.26 -0.41 203.56 -0.02
NL 0.05 0.08 0.05 0.08 0.05 0.11 0.00
------------ ------- ------- ------- ------- ------- ------- -------
Lewis 51.84 51.86 51.91 51.87 51.91 51.93 26.16
non-Lewis 15.86 15.89 15.87 15.84 15.84 15.84 -1.38
------------ ------- ------- ------- ------- ------- ------- -------
Total 67.69 67.75 67.78 67.71 67.76 67.77 24.77
這里C1-C2出現兩次是因為它們之間被認為是雙鍵。
如果在NCS后面的<>里寫上MO，則會輸出每個原子屏蔽張量來自于每個占據正則分子軌道的貢獻，如下所示。還會緊跟著輸出在主軸方向下的值，最后也會輸出匯總信息。
Full Cartesian NMR shielding tensor (ppm) for atom gh(13):
Canonical MO contributions
MO XX YZ ZX XY YY ZY XZ YZ ZZ
===============================================================================
1. 0.12 0.00 0.00 0.00 0.12 0.00 0.00 0.00 0.22
2. 0.12 0.00 0.00 0.00 0.12 0.00 0.00 0.00 -1.20
3. 0.12 0.00 0.00 0.00 0.12 0.00 0.00 0.00 -1.20
4. 0.07 0.00 0.00 0.00 0.22 0.00 0.00 0.00 -2.99
5. 0.22 0.00 0.00 0.00 0.07 0.00 0.00 0.00 -2.99
......略
19. 1.27 0.00 0.00 0.00 8.03 0.00 0.00 0.00 21.90
20. -0.44 0.00 0.00 0.00 4.87 0.00 0.00 0.00 -5.07
21. 4.85 0.00 0.00 0.00 -0.45 0.00 0.00 0.00 -5.08
-------------------------------------------------------------------------------
Total 41.91 0.00 0.00 0.00 41.93 0.00 0.00 0.00 70.75
-------------------------------------------------------------------------------
遺憾的是NBO5.0結合G03 C.02版給出的結果我發現是明顯錯誤的，根本沒法用，總的XX、YY、ZZ值明顯和Gaussian算出來的以及分解為NBO貢獻時的總值不符。結合G98得到的結果肯定是正確的，如果是G03/09的用戶，可能得用更新的NBO版本才行，但官方主頁上宣稱只支持G9X。
3 使用NBO 3.1分解磁屏蔽張量
在G03/09中只要用了NMR關鍵詞，并且同時用了比如pop=nbo關鍵詞調用了NBO 3.1分析模塊，則NBO 3.1模塊就會自動進行NCS分析，例如：
#P b3lyp/6-31+G* NMR pop=nbo
b3lyp/6-31+g(d) opted
0 1
C 0.00000000 1.39864200 0.00000000
C 1.21125900 0.69932100 0.00000000
C 1.21125900 -0.69932100 0.00000000
C 0.00000000 -1.39864200 0.00000000
C -1.21125900 -0.69932100 0.00000000
C -1.21125900 0.69932100 0.00000000
H 0.00000000 2.48606800 0.00000000
H 2.15299800 1.24303400 0.00000000
H 2.15299800 -1.24303400 0.00000000
H 0.00000000 -2.48606800 0.00000000
H -2.15299800 -1.24303400 0.00000000
H -2.15299800 1.24303400 0.00000000
bq 0. 0. 0.
bq 0. 0. 1.
這里不必像NBO5.0那樣用混合基組的技巧解決與Gaussian兼容性問題，NBO3.1對虛原子的處理和Gaussian是兼容的。
例如，在NBO模塊輸出信息末尾可以找到各個NLMO對苯環中心的虛原子的磁屏蔽張量：
Contributions to Isotropic shielding for atom 13 ( 7.9768 ppm):
Standard Origin Atomic Origin
NLMO Total Paramagnetic Diamagnetic Paramagnetic Diamagnetic
1 -1.3247 -3.4073 2.0825 -1.3247 0.0000
2 -1.3248 -3.4072 2.0825 -1.3248 0.0000
3 8.1896 3.0323 5.1573 8.1896 0.0000
4 -0.8489 -1.4446 0.5957 -0.8489 0.0000
5 -1.3247 -3.4080 2.0833 -1.3247 0.0000
6 8.1936 3.0366 5.1570 8.1936 0.0000
7 -0.8492 -1.4451 0.5959 -0.8492 0.0000
8 -1.3247 -3.4073 2.0825 -1.3247 0.0000
9 -0.8492 -1.4451 0.5959 -0.8492 0.0000
10 -1.3248 -3.4072 2.0825 -1.3248 0.0000
11 8.1896 3.0323 5.1573 8.1896 0.0000
12 -0.8489 -1.4446 0.5957 -0.8489 0.0000
13 -1.3246 -3.4080 2.0835 -1.3246 0.0000
14 -0.8492 -1.4451 0.5959 -0.8492 0.0000
15 -0.8492 -1.4451 0.5959 -0.8492 0.0000
Sum 11.5298 -20.0134 31.5433 11.5298 0.0000
順磁和抗磁成份這里也單獨給了出來，但是一般只需要關注它們的總效應，即Total那列就行了。
可見這里Sum的值11.5298ppm和Gaussian算出來的各向同性磁屏蔽值7.9768對應不上，這是因為默認情況下，貢獻較小的NLMO不會被輸出出來，而Sum值僅僅是所有全部顯示出來的貢獻值的加和，這里只顯示出了全部21個NLMO中的15個的貢獻。由于沒有算進去那些沒顯示出來的，所以總值自然對應不上。這個默認設定非常傻。為了讓所有NLMO的貢獻都顯示出來，需要加上IOp(6/65=-1)。
默認情況下輸出的只是各向同性磁屏蔽值來自各NLMO的貢獻，如果想要將NLMO對磁屏蔽張量的各個元素的貢獻也輸出出來，就要用pop=NCSall關鍵詞，而pop=nbo就不必寫了。例如上面例子的route section寫成
#P b3lyp/6-31+G* NMR pop=NCSall IOp(6/65=-1)
在輸出文件中就能找到諸如下面這樣的信息，全部NLMO對環中心的磁屏蔽張量ZZ分量的貢獻都給出了。各個NLMO對NICS(0)_ZZ的貢獻也就是Total那列的負值。
Contributions to Atom 13 Sigma-ZZ ( 13.1863 ppm):
Standard Origin Atomic Origin
NLMO Total Paramagnetic Diamagnetic Paramagnetic Diamagnetic
1 -0.8999 -5.9092 5.0093 -0.8999 0.0000
2 12.0045 5.7848 6.2197 12.0045 0.0000
3 -0.8999 -5.9092 5.0093 -0.8999 0.0000
4 -1.6607 -1.3298 -0.3309 -1.6607 0.0000
5 -0.9002 -5.9111 5.0108 -0.9002 0.0000
6 -1.6618 -1.3315 -0.3303 -1.6618 0.0000
7 -0.8999 -5.9092 5.0093 -0.8999 0.0000
8 12.0045 5.7848 6.2197 12.0045 0.0000
9 -1.6618 -1.3315 -0.3303 -1.6618 0.0000
10 -0.8999 -5.9092 5.0093 -0.8999 0.0000
11 -1.6607 -1.3298 -0.3309 -1.6607 0.0000
12 -0.9003 -5.9111 5.0108 -0.9003 0.0000
13 12.0057 5.7863 6.2194 12.0057 0.0000
14 -1.6618 -1.3315 -0.3303 -1.6618 0.0000
15 -1.6618 -1.3315 -0.3303 -1.6618 0.0000
16 -1.2438 -0.7266 -0.5172 -1.2438 0.0000
17 -1.2431 -0.7258 -0.5172 -1.2431 0.0000
18 -1.2431 -0.7258 -0.5172 -1.2431 0.0000
19 -1.2438 -0.7266 -0.5172 -1.2438 0.0000
20 -1.2431 -0.7258 -0.5172 -1.2431 0.0000
21 -1.2431 -0.7258 -0.5172 -1.2431 0.0000
Sum 13.1863 -30.4450 43.6313 13.1863 0.0000
可以按照《使用Multiwfn繪制NBO及相關軌道》（http://www.shanxitv.org/134）來觀看這些NLMO的圖形，找出sigma和pi軌道。最好不要用pop=saveNLMO將NLMO保存到check文件中去觀看，而建議用plot文件方式來觀看，以免自動對軌道重新排序而不便與NCS輸出信息對應。通過圖形觀察可知2、8、13號NLMO都是C-C pi軌道，磁屏蔽作用很強，每個都達到12ppm，因此NICS(0)_ZZ_pi值即是3*-12=-36ppm。有興趣的讀者可以參看Org. Lett., 8, 863文章的補充材料的末尾，作者是利用NBO5.G結合G98將磁屏蔽張量分解為CMO的貢獻并計算的NICS(0)_ZZ_pi，他的結果和本文的值相符很好。
總的來說，用Gaussian自帶的NBO3.1來做NCS比起用NBO5.0來做明顯更方便，而且也不用買NBO5.0，也不涉及Gaussian源代碼版權問題。