通過Multiwfn繪制等化學屏蔽表面(ICSS)研究芳香性

注1：經常有Multiwfn用戶閱讀本文不仔細導致操作不當，無法正確繪制出ICSS，故這里提供苯的ICSS計算涉及到的所有文件：http://www.shanxitv.org/attach/216/benzene.rar。繪制不成功者請先嘗試直接用這里的文件進行繪制，并仔細對照這里面的文件和自己算出來的文件。

注2：本文很多例子中利用了用戶自定義函數-1來基于ICSS_ZZ格點數據做三線性插值繪制平面圖、曲線圖。而后來Multiwfn支持了更好的三維B-spline插值，得到的插值函數以及相應的圖像明顯更平滑，故建議把本文中設iuserfunc為-1的地方都改為-3來使用B-spline插值。

2023-Aug-8重要補充：Multiwfn后來專門加入了繪制一維NICS曲線圖和二維NICS平面圖的功能，比起用本文介紹的ICSS格點數據插值繪制方式步驟明顯更簡單，而且更關鍵的是耗時低得多得多。見《使用Multiwfn繪制一維NICS曲線并通過積分衡量芳香性》（http://www.shanxitv.org/681）和《使用Multiwfn巨方便地繪制二維NICS平面圖考察芳香性》（http://www.shanxitv.org/682）。因此如果你只是要繪制曲線圖或平面圖，強烈推薦用這兩篇文章里的做法而不建議用本文的做法。

 

通過Multiwfn繪制等化學屏蔽表面(ICSS)研究芳香性

文/Sobereva @北京科音
First release: 2013-Dec-30  Last update: 2023-Aug-11

1 等化學屏蔽表面的含義

核獨立化學位移(NICS)，以及其改進版NICS_ZZ是研究芳香性問題的極其重要的手段，可參見這兩個博文中的介紹和討論：《衡量芳香性的方法以及在Multiwfn中的計算》（http://www.shanxitv.org/176）、《將核獨立化學位移(NICS)分解為sigma和pi軌道的貢獻》（http://www.shanxitv.org/145）。

通常NICS的研究是0維的，也就是只考察特定的點上的NICS數值，最典型的就是研究環中心的位置。也有一些文章將NICS的研究擴展到一維，討論在一條線上NICS值的變化。還有少數文章，在二維層面上討論NICS，例如取一個垂直于并穿越環的截面，考察NICS值在這個平面各個位置上是如何分布的。而本文介紹的相當于是在三維層面上研究NICS的方法，也就是計算磁屏蔽值在三維空間中的格點數據，并繪制成等值面圖，這被稱為等化學屏蔽表面（ICSS, iso-chemical shielding surfaces），這是在J. Chem. Soc., Perkin Trans. 2, 2001, 1893-1898當中被提出的。這種研究方法比起傳統的計算環中心NICS數值的研究方式直觀形象得多，并且同時展現出更豐富的信息。

原文定義的ICSS對應的是NICS的原始定義，體現的是對外磁場的各項同性(Isotropic)屏蔽。對于平面體系，Org. Lett., 8, 863 (2006)已證明NICS_ZZ比起NICS更具有優勢，更適合展現芳香性，因此筆者也定義了ICSS_ZZ，也就是相當于NICS_ZZ的等值面。這里Z方向是指垂直于環平面的方向，已經假設分子平面是在XY平面上，顯然如果分子平面在YZ方向就要研究NICS_XX和ICSS_XX。另外筆者還定義了ICSS_ani，對應的是NICS_ani等值面。NICS_ani這里指的是磁屏蔽各向異性(anisotropy)值，定義為ε_3 - (ε_1 + ε_2)/2，ε代表磁屏蔽張量的本征值，由小到大排序。各向異性值越大的地方說明這個地方對不同方向外磁場的屏蔽能力差異越大。

特別注意的是，ICSS和NICS的符號定義是恰好相反的。NICS取的是磁屏蔽值的負值，也就是說，NICS越負，在這個點上對外磁場屏蔽程度越大，NICS越正則說明這個位置的去屏蔽程度越強。而ICSS直接展現的就是不同位置的磁屏蔽值，并沒有取負號。

功能最強大的波函數分析程序Multiwfn從3.2.1版開始加入了繪制ICSS、ICSS_ZZ（或ICSS_XX、ICSS_YY）和ICSS_ani的功能，是目前唯一公開的能夠繪制ICSS的程序，可在此處免費下載http://www.shanxitv.org/multiwfn。在本文下面將通過幾個簡單的實例介紹如何進行操作，同時討論ICSS的意義。同時還介紹如何通過Multiwfn研究一條線上磁屏蔽值的變化，或研究特定平面上的磁屏蔽值的分布。本文所用的.gjf文件在Multiwfn程序的examples\ICSS目錄下都已經提供了，平面或準平面分子的平面都是處在XY平面上。

如果你的文章使用了ICSS，應引用其原文J. Chem. Soc., Perkin Trans. 2, 2001, 1893。如果你的文章使用了筆者提出的ICSS_ZZ、ICSS_ani，應引用ICSS_ZZ研究18碳環的文章Carbon, 165, 468 (2020)。也推薦引用筆者的Chem. Eur. J., 28, e202300348 (2023) DOI: 10.1002/chem.202300348、J. Phys. Chem. C, 123, 18593 (2019)、Carbon, 165, 468 (2020)、Chem. Eur. J., 28, e202103815 (2022)，都是非常好的ICSS的應用范例。當然，也別忘了必須同時引用Multiwfn程序的原文。

 

2 苯

Multiwfn自身沒法做NMR計算得到空間中各個位置上的磁屏蔽張量。Gaussian雖然能做NMR計算，只要通過Bq原子就可以指定要計算磁屏蔽張量的位置，但是沒法通過這種方式直接產生磁屏蔽值的格點數據。另外，Gaussian程序自身有限制，或者某種程度說存在bug，如果要計算的格點數據包含幾萬或者更多的點，那么就無法一次性把這些點都作為Bq寫進一個Gaussian輸入文件里，否則程序會報錯，而必須分成多次處理，比如每個輸入文件計算8000個點的屏蔽張量的。用Multiwfn研究ICSS時，程序會讀取分子結構并由用戶對格點進行設定，然后會生成一批Gaussian的NMR任務的輸入文件，在使用Gaussian對它們計算后，Multiwfn將把Gaussian輸出文件里的磁屏蔽張量信息進行匯總并轉換成格點數據，然后就可以觀看ICSS了。

先看最簡單的例子，苯分子。首先準備一個當前體系的標準的Gaussian輸入文件，形式如下所示
%chk=./benzene.chk
# b3lyp/6-31+g(d)
[空行]
test
[空行]
0 1
 C                  0.00000000    1.39078800    0.00000000
 C                  1.20445700    0.69539400    0.00000000
 C                  1.20445700   -0.69539400    0.00000000
 C                  0.00000000   -1.39078800    0.00000000
 C                 -1.20445700   -0.69539400    0.00000000
 C                 -1.20445700    0.69539400    0.00000000
 H                  0.00000000    2.47288200    0.00000000
 H                  2.14157800    1.23644100    0.00000000
 H                  2.14157800   -1.23644100    0.00000000
 H                  0.00000000   -2.47288200    0.00000000
 H                 -2.14157800   -1.23644100    0.00000000
 H                 -2.14157800    1.23644100    0.00000000

注意必須根據實際情況明確定義%chk。體系應當事先已經優化好，坐標以笛卡爾坐標表示。關鍵詞里只需指定理論方法和基組（也可以加上一些幫助SCF收斂的選項之類），不要寫opt、freq之類涉及任務的關鍵詞。此文件里這些關鍵詞將成為之后Multiwfn產生的做NMR任務的Gaussian輸入文件里的關鍵詞。

值得一提的是在使用Bq后無法開啟對稱性，Gaussian也因此不會自動調整體系朝向并由此帶來坐標不對應的問題，所以不需要加上nosymm來手動關閉對稱性，因為nosymm此時已經是默認的。計算NMR用B3LYP/6-31+G*是個比較不錯的選擇，不昂貴，精度也還不錯。如果體系較大，算ICSS格點數據難以算得動，可以降到6-31G*，但不宜進一步降低。

啟動Multiwfn后輸入
examples\ICSS\benzene.gjf   //這個自帶的文件就是前面給出的那個gjf文件
25  //離域性和芳香性分析
3   //ICSS分析
-10   //設定格點數據計算區域的延展距離
12    //延展距離設為12 Bohr，也就是相對于邊界原子向四周擴12 Bohr。延展距離不能太小，否則繪制出來的ICSS等值面會被截斷。通常12 Bohr夠大了
1     //低質量格點。Gaussian接下來將總共計算約125000個位置的磁屏蔽張量。如果用更高質量的格點計算耗時將會很高
n    //不跳過Gaussian輸入文件生成的步驟
現在Multiwfn以benzene.gjf文件作為模板，在當前目錄下生成了NICS0001.gjf、NICS0002.gjf...共17個Gaussian NMR任務的輸入文件，大家有興趣可以自行打開看看了解一下內容，可見里面除了原子坐標外就是以Bq表示的每個格點的坐標。除了最后一個gjf文件外，每個文件中體系的原子數和Bq數加起來都為8000個。然后大家自行通過Gaussian執行它們，其中NICS0001.gjf必須作為第一個執行，因為它會產生chk文件，之后的gjf里都有guess=read關鍵詞從這個chk里讀取波函數作為初猜來節約時間。為了批量運行省事，在examples目錄下提供了兩個腳本runall.sh和runall.bat，前者和后者分別對應Linux和Windows系統，執行它們之后就可以調用g09批量執行當前目錄下所有.gjf文件，并在當前目錄下產生名字相同的.out文件。注意Windows版用戶別忘了設定Gaussian的環境變量，否則調用Gaussian時會出錯，設置方法參見Multiwfn手冊附錄1。

注：如果Gaussian運行時馬上報錯，可能是由于每個輸入文件里原子數（包括Bq）太多所致。每個輸入文件里的原子數的上限和內存大小以及Gaussian的版本密切相關。為解決此問題，可以將Multiwfn的settings.ini文件中的NICSnptlim參數改小些，比如改為5000，然后重啟Multiwfn并重新生成gjf文件，此時每個gjf文件里實際原子和Bq數之和將為5000。由于總的點數是固定的，NICSnptlim越小則產生的gjf文件數目越多，總計算耗時也會越長（但也不是說NICSnptlim越大總耗時就一定越低，比如設成10萬個反倒總運行時間還更長）。如果是G09 D.01版，一定要加上guess=huckel，否則NICSnptlim必須設得特別小才行。如果是G09 E.01版，.gjf直接運行不成功，用guess=huckel還是出錯，則改用guess=core再試。G16不需要添加guess關鍵詞。

用Gaussian運算后，假設生成的NICS0001.out、NICS0002.out...都放在了C:\ltwd\benzene目錄下，那么在Multiwfn窗口中接著輸入C:\ltwd\benzene\NICS，然后選擇你感興趣的量，Isotropic、Anisotropic分別代表將要產生各項同性和各向異性磁屏蔽值的格點數據，XX component、YY component和ZZ component分別代表將要產生X、Y、Z方向屏蔽值的格點數據。然后Multiwfn就會讀取這些.out文件并把里面的磁屏蔽張量信息收集起來，產生你指定的量的格點數據。接下來，可以選-1重新選擇并生成感興趣的量的格點數據，可以選0退回出菜單，可以選1來觀看你感興趣的量的等值面，或者選2來把格點數據導出成cube文件以便用其它可視化工具來觀看等值面。這里假設我們選的量是Isotropic，然后又選了1觀看等值面，就會進入到一個圖形界面，把等值面設為0.25后，就會看到下圖，即ICSS

旋轉一下體系，并從窗口上方的菜單選Isosurface style - Use transparent surface，可以得到俯視圖

綠色等值面代表磁屏蔽值為0.25ppm的區域。藍色等值面代表屏蔽值為-0.25ppm，或者說去屏蔽值為0.25ppm的區域。可以看出，在垂直于苯環的方向，等值面凸出來一大塊，說明在苯的內部區域磁屏蔽較強，這是由于苯具有極強的芳香性，pi電子可以自由運動，外磁場令pi電子產生的磁感應環電流會在苯環內部區域顯著抵消外磁場。同時從俯視圖也可以清楚看到，即便是C-H鍵的區域，和pi電子沒什么關系，外磁場卻也被明顯屏蔽了，因而也被綠色等值面所包圍。其原因是構成C-H鍵的sigma電子在C-H鍵區域形成了局部感應環流，所以C-H鍵區域的外磁場也被很大程度屏蔽了，如下圖所示。下圖左邊是分子平面上的情況，對應的是sigma電子的感應電流密度，每個sigma鍵都有局部環流；而右圖是分子平面上方一定距離的圖，體現的主要是高度離域的pi電子的感應電流密度。

注：若想作電流密度圖，可以參見《使用AICD程序研究電子離域性和磁感應電流密度》（http://www.shanxitv.org/147）和《使用GIMIC計算和分析磁感應電流密度》（http://www.shanxitv.org/155）。

為何ICSS圖中在苯環四周出現了一圈明顯的去屏蔽區域？從下面這張外磁場、感應電流和感生磁場的關系就很容易理解。

在加了外磁場（紅色箭頭）后，苯環整體形成了黃圈所示的pi磁感應環流，在苯環里側，紫色所示的感生磁場方向和外磁場相反，產生了屏蔽。而在苯環外側，pi環流產生的磁場和外磁場方向一致，對它產生了增強效果，因此苯環外側對磁場是去屏蔽的。

我們看看ICSS_ZZ和ICSS圖有什么區別。接著之前的分析，關閉GUI窗口后在Multiwfn窗口里選-1 Load another ICSS form，然后選5: ZZ component，Multiwfn會重新讀取那些Gaussian的輸出文件并產生ICSS_ZZ的格點數據，選1查看等值面并把等值面數值改為2，結果如下所示

可見，和ICSS的圖大同小異。但是此時的等值面綠色和藍色表現的是對Z方向的外磁場產生2ppm的屏蔽和去屏蔽的區域，物理意思和ICSS有一定區別。對于平面體系的芳香性研究，筆者建議用ICSS_ZZ代替ICSS，道理和NICS_ZZ優于NICS一樣。因為平面體系的芳香性應當靠對于垂直于分子平面的磁場的屏蔽能力來表現，而其它分量的影響理應是被去除的。雖然在此例中ICSS和ICSS_ZZ沒太大區別，但是換成其它體系，差異有可能是定性的。

假設用戶已經有了Gaussian的輸出文件，想直接做ICSS分析，對于當前體系可以按以下步驟操作
examples\ICSS\benzene.gjf
25  //離域性和芳香性分析
3   //ICSS分析
-10   //設定格點數據計算區域的延展距離
12    //延展距離設為12 Bohr
1     //低質量格點
y    //跳過Gaussian輸入文件生成的步驟
C:\ltwd\benzene\NICS      //將載入C:\ltwd\benzene\NICS0001.out、NICS0002.out...
1   //選感興趣的量，例如選Isotropic
可見，重做ICSS分析的時候不需要再重新生成Gaussian輸入文件了，但是在格點設定時必須和之前生成Gaussian輸入文件時的格點設定完全一致，否則在載入之前得到的Gaussian輸出文件時就會不對應，導致Multiwfn崩潰或者結果無意義。

3 環丁二烯

做法和苯的例子一樣，請大家基于examples\ICSS\下的cyclobutadiene.gjf自行完成繪制。所得ICSS=0.5的圖如下。為了清楚起見，用了兩種顯示方式

環丁二烯的情況和苯截然相反。環的內側和上下部分都是去屏蔽區域，而環的周圍都是屏蔽區域。這是因為環丁二烯是典型的4n電子反芳香性體系，pi環電流的方向和苯的情況是相反的。

可見，ICSS方式的研究比起計算特定點的NICS值生動得多，信息也更為豐富，可以掌握整個空間各個位置上對磁場的屏蔽或去屏蔽情況。而且也能避免對于一些結構復雜的體系選取計算NICS點的位置的困難。

4 環庚三烯(cycloheptatriene)

環庚三烯這個體系結構如下

雖然也是4n+2電子，但是表面上看共軛路徑上隔著一個sp3的碳，并且是非平面的，于是Wiki英文百科上說此分子是非芳香性的，這種說法靠譜否？如果有芳香性，對外磁場又是如何屏蔽的？按照前述步驟算一下這個體系的ICSS，結果如下，ICSS=0.15和ICSS=0.25的都給出了

0.15的圖和苯的情況很相似，也是環內側有明顯屏蔽，周圍有去屏蔽區域。但是，如果把等值面數值增加到和苯的圖一樣，即0.25，那么去屏蔽的環狀等值面就明顯收縮并且斷開了。這樣的觀測表明環庚三烯確實有芳香性，能夠形成整體的環流而對內部區域的磁場造成明顯屏蔽，但是程度要比苯更弱。我們也可以進一步通過Multiwfn計算多中心鍵級來考察一下。苯在B3LYP/6-31G*下的多中心鍵級為0.086，如果我們取環庚三烯的六個含有pi電子的碳（即1-2-3-5-6-4）在同樣條件下計算多中心鍵級，結果為0.03，所以結論和ICSS的一致，環庚三烯有芳香性但比苯要弱。

5 丙烷

這回來看一個典型的非芳香性分子，它的ICSS=0.12、0.08和0.02的圖如下

我們先看0.08的等值面，從中看到，分子主體被綠色等值面包圍，這是sigma鍵的電子形成的局部環流產生的磁屏蔽效應。有屏蔽區域就伴隨著去屏蔽區域，因為環電流不可能使得任何位置都處于屏蔽區域，所以在旁邊還有藍色去屏蔽區域。但是去屏蔽效應很弱，等值面數值在0.08時才看到一點而已，若把等值面數值調大到0.12就沒有藍色區域了。如果把等值面數值降低到0.02，則可以看到去屏蔽區域涉及空間范圍其實很大，不過這些區域內去屏蔽程度非常弱。

6 pyracylene

這個體系光從下面的結構式來看很不好說這個分子是芳香性的還是反芳香性的，亦或是局部芳香性/反芳香性。我們通過ICSS_ZZ來研究下

計算過程為
examples\ICSS\pyracylene.gjf
25
3
-10
12
4   //自定義每個方向的格點間距或格點數
0.5,0.5,0.5   //每個方向格點間距為0.5 Bohr
n
之后和前面的例子一樣，也是生成Gaussian輸入文件（共35個），然后自行計算它們，再把它們載入到Multiwfn里。這個例子和之前唯一不同的是我們沒用low quality grid，而是自定義了格點間距。因為low quality grid是指總點數約為125000個，但當前這個體系略大，因此計算的空間區域比前幾例也都更大，因此如果還用low quality grid的話格點密度就偏低了，等值面顯示效果可能就差了。所以這里我們直接指定格點間距，0.5bohr間距基本算是達到精度和計算量的平衡點。如果是更大的體系，諸如富勒烯，那么為了省時可以把格點間距設得更大一些。

這個體系的ICSS_ZZ圖如下

可見五元環是反芳香性的，而六元環有一定芳香性。定量討論不同環的芳香性往往我們需要確定一個平面，例如分子平面上方1埃處，看不同環的NICS(1)ZZ值。而光從等值面圖上我們難以比較某個平面上的屏蔽值分布。Multiwfn是十分靈活的程序，繪制指定平面上的屏蔽值，包括填色圖、等值線圖、地形圖等等在Multiwfn中都可以做到。首先我們選2 Export the grid data to ICSSZZ.cub current folder把當前格點數據導出到ICSSZZ.cub文件中。然后把settings.ini里的iuserfunc改為-1，這樣自定義函數就變為了通過格點數據插值獲得的數值。啟動新的Multiwfn，然后輸入
ICSSZZ.cub
4    //繪制平面圖
100     //自定義函數。此時這個函數在任意一點的數值即是基于ICSSZZ.cub插值得到的值
1    //填色圖
[回車]
1    //XY平面
1.89    //Z=1.89 Bohr，相當于Z=1埃
在剛出現的圖上點右鍵關閉之，并輸入
1    //更改色彩刻度
-40,30
4    //顯示原子標號
1    //文字為紅色
17    //設定顯示文字的距離閾值
2     //只要原子核距離作圖平面小于2埃則這個原子的符號就顯示在圖上
n
2    //顯示等值線
-1    //重新作圖
此時會看到下圖

從這張圖上我們可以清楚看出，在五元環上方1埃處，去屏蔽非常嚴重，達到-40ppm左右。實際上精確計算的五元環的NICS(1)ZZ值就是40.3476。而六元環上方1埃處，對應的顏色是綠色的，對照色彩刻度條就知道屏蔽值很低，也就比0大不了多少，實際上精確計算的六元環的NICS(1)ZZ值僅為-1.5543。我們還可以通過圖形間接地看出六元環的pi電子分布很不均勻，其中邊緣的C-C鍵上屏蔽值較高，體系中間的C-C間屏蔽值也較高，而邊緣與中間的碳之間的屏蔽值較低，這表明pi電子在每個六元環上沒有充分地離域化，而是分離成了兩個定域性相對較高的區域，電子容易在這兩個區域內形成環流。也正因為pi電子缺乏在六元環上的充分離域，沒有使得六元環中心的上方出現較高的屏蔽值。

順便作為對比，這里給出苯的分子平面上方1埃處的Z方向磁屏蔽值的填色圖以茲對比。可見苯上方磁場被屏蔽得非常顯著，相比之下，pyracylene的六元環的芳香性就很不顯著了

請大家自行研究pyracylene的ICSS，與上面討論的ICSS_ZZ進行對比，在此例它們的差異還是挺明顯的。另外請自行計算分子平面(Z=0)和分子平面上方2埃(Z=2)的XY平面的磁屏蔽值的填色圖。

通過Multiwfn我們還可以研究從五元環和六元環中心開始，沿著分子平面垂直向上（即沿著Z方向）的磁屏蔽值的變化。首先我們先求出五元環和六元環的中心。返回主菜單后，選100，再選21，然后輸入五元環的原子編號4,5,20,19,9，程序就給出了幾何中心
Geometry center (X/Y/Z):    0.00000000   -1.84880000    0.00000000 Angstrom
（以Bohr為單位相當于0,-3.49,0）
然后輸入
q   //返回
0   //返回主菜單
3   //繪制函數在一條直線上的變化曲線
100   //自定義函數
2    //自定義這條線的兩端坐標
0,-3.49,0,0,-3.49,12    //曲線的始端為五元環中心(0,-3.49,0），末端為向Z方向移動10Bohr的位置(0,-3.49,12)
馬上出現下面的圖像，虛線代表y=0的位置

從圖可見，五元環中心的Z方向磁場去屏蔽值極大，隨著沿著垂直于分子平面方向移動，去屏蔽能力迅速下降，在5Bohr左右屏蔽值恰好為0。然后繼續上升，屏蔽值成為微小的正值，再之后則十分緩慢地下降。成為正值絕非表明五元環呈現芳香性，而是由于其余區域電子造成的磁屏蔽效應恰延伸到了五元環上方較遠的地方。大家作一下等值面數值較小的ICSS_ZZ圖，比如ICSS_ZZ=2的圖就能看到，此時藍色的去屏蔽等值面被綠色的等值面給淹沒了。之所以后來曲線又開始緩慢下降，是因為離體系越遠，那些電子產生的磁屏蔽能力逐漸削弱。Multiwfn的文本窗口中可以看到一些有用的信息，比如曲線上最大、最小值
Minimal/Maximum value: -0.71408290D+02  0.13046734D+01
關閉圖像后，如果用選項7，我們可以得到Y為指定值時X的位置，例如我們選7并輸入0，結果為4.92Bohr，即屏蔽值恰為0的位置是五元環上方4.92Bohr處。如果用選項6，程序可以給出曲線上極大極小點，我們選擇6之后程序發現此曲線上有一個極大點，位置是7Bohr處，屏蔽值為1.3ppm。

我們選2把曲線數據導出到當前目錄下的line.txt，并改名為five.txt。然后再效仿上面的步驟，把六元環在垂直于環平面上的磁屏蔽值曲線作出來，同樣導出到line.txt。然后在origin里利用這兩個數據文件，作成同一張曲線圖以便于對比，如下所示。注意在導出曲線數據時Multiwfn窗口上提示導出的數據文件中有5列，前三列是曲線上每個點的實際的X,Y,Z坐標，這對作圖沒用，在origin里作曲線圖時要將最后兩列分別用作X和Y，坐標單位為埃。

黑色曲線是六元環的結果，紅色是五元環的，剛才已經見過了。綠色虛線把曲線與Y=0的相交的位置標明了。我們發現，六元環在分子平面處實際上磁屏蔽值也是負值，但是到了分子平面上方1埃處屏蔽值已經升到了正值（NICS(1)ZZ為-1.5543），隨后在1.6埃處達到了最高值，約6點幾ppm。隨后隨著距離平面漸遠，屏蔽能力一直下降。所以，六元環有一些芳香性，且是pi電子引起。不過倒是不能武斷地說六元環有sigma反芳香性，因為這可能是Paratropic電流太弱而不敵Diatropic電流所致的（詳見PCCP,13,20500關于這兩類感生電流的討論）。

另外，由于我們已經有了ICSS_ZZ的格點數據，在Multiwfn中我們還可以直接通過插值計算指定點的數值。例如我們想計算3,-2.2,0.4處的值，就返回主菜單，選擇主功能1，輸入3,-2.2,0.4，然后選擇輸入的坐標單位，程序輸出的User-defined real space function值就是這個點的Z方向磁屏蔽值了。當然受制于格點數據精度，這么算的結果和直接在Gaussian里用Bq指定位置來精確計算的會有一點偏差，但可以忽略。顯然，只要我們已經有了ICSS_ZZ格點數據，想獲得NICS(0)ZZ、NICS(1)ZZ就根本沒必要再單獨做一次NMR計算了，因為你只要在主功能1里輸入環中心、環中心上方1埃處的坐標，將輸出的User-defined real space function取負值后直接就分別是NICS(0)ZZ和NICS(1)ZZ了。

上面的曲線圖，以及填色圖的等值線上，都可以看到在某些部位不是很光滑，這是因為ICSS_ZZ的格點數據的格點密度不是很高，所以插值出來的結果準確度有限。想改進圖像質量和精確度就得在ICSS計算過程中設定格點時用更密的格點。另外注意，插值不能超過格點數據的空間范圍。例如這一節的例子，計算的ICSS_ZZ.cub的延展距離只有12 Bohr，所以想計算從環中心向垂直于環平面上延伸到20 Bohr的曲線圖是不可能的，到了12 Bohr時曲線就會自動成為0，沒有意義。

7 薁(azulene)

azulene的結構如下

這也是個不做計算很難估算芳香性的體系。我們就來算一算ICSS_ZZ。作圖過程和上一節一樣。ICSS_ZZ=2的等值面圖如下

可見這個分子顯示出明顯的芳香性，無論是五元環還是七元環都被明顯屏蔽，周圍出現了和苯一樣的一圈明顯的去屏蔽區域。如果我們想研究五元環和七元環誰的芳香性更強，我們可以逐漸增加等值面數值，并以各種角度進行觀察。ICSS_ZZ=30時的側視圖如下，左半邊是五元環右半邊是七元環，可以看出左半邊的等值面比右半邊鼓出一塊，這表明五元環的芳香性要強過七元環，因為被屏蔽得更厲害。如果繼續增加ICSS_ZZ值到30，得到下圖的右圖，可見此時在七元環上已經出現了窟窿，而五元環上還是被等值面覆蓋著，也說明五元環芳香性顯著更強。

效仿上一節，作分子平面上方1埃處的填色圖，如下所示，結果更明確了，七元環的對Z方向磁場的屏蔽能力比五元環差得遠

我們可以進而再作Z方向磁屏蔽值的橫截面圖，作圖平面垂直于體系，并且穿過C1和C15原子，也就是X=0的YZ平面的填色圖，結果如下。圖中左半部分，也就是五元環部分，在環的中心區域雖然也有些對Z方向磁場的屏蔽能力，但是明顯不如環上下區域屏蔽能力高，這是表現出這個環呈現pi芳香性而沒有sigma芳香性。而圖的右半邊，顯著體現出七元環的屏蔽能力很弱，環中心區域屏蔽值甚低，環的上下方只是比它微微強一點點。

在ICSS的文獻中，總是看到那種同時顯示多個不同數值的ICSS等值面的圖，并且做出了截面效果。每層截面其實本質上相當于上面的填色圖上的等值線。不過通過多個等值面+截面方式顯示，比起上面給出的填色圖有時會顯得更形象生動一些。這里就介紹一下怎么結合VMD 1.9.1作這種圖。

首先，在獲得當前體系的ICSS_ZZ格點數據后將之導出到ICSSZZ.cub文件里，然后啟動VMD，將這個文件拖入至VMD主窗口。然后選Graphics-representations，然后點Create Rep新建顯示方式，Drawing Method改為Isosurface，Isovalue改為2，Draw改為Solid surface，Show改為Isosurface，Coloring Method改為ColorID，并且選一種喜歡的顏色比如Red。之后，再次點擊Create Rep并進行同樣的操作，但是Isovalue改為4，ColorID改為另一種顏色。然后繼續不斷創建新顯示方式，每次Isovalue數值增加一倍直到32，ColorID隨意設定使之與之前的不同。最后再加一個等值面顯示方式，令數值為-1，專用來表現去屏蔽的區域。然后選Extensions-Visualization-Clipping plane tool，點Active按鈕，取消Normal follow view復選框。由于我們這次作的是平行于YZ平面的截面，所以把Normal（法矢量）改為1 0 0。雖然效果已經有了，但是分子自身也被截斷了。因此，我們再次把ICSSZZ.cub拖進VMD主窗口，可以看到VMD主窗口中就又多了一個ID，并且開頭寫著T（Top）字樣。切換回Clip Tool窗口，把Active按鈕取消掉，此時分子結構沒被截斷了，我們再修改分子結構顯示方式。進入Graphics-representations，確認第一欄已經切換到了ID=1的體系，把默認的Lines顯示方式改為CPK。然后控制臺窗口運行color Display Background white把背景改為白色，便得到下圖

可見，如果只是研究一個數值的等值面，局限性是很大的。比如這張圖上紅色對應屏蔽值為2ppm。如果只考察這個等值面，結論只是體系整體有芳香性，卻無法得知五元環和七元環各自芳香性的強弱。而到了棕綠色，即16ppm的ICSS_ZZ等值面上，我們才能發現五元環芳香性更強，從32ppm的ICSS_ZZ的等值面上這點體現得更為明顯。所以，在通過ICSS討論一個體系的磁屏蔽問題時，應當逐漸調整等值面數值來考察，以免有失偏頗。而將圖片用于文獻時，建議給出上圖這樣的多層ICSS等值面圖，或者給出關鍵性截面上的填色圖，這樣才能將體系的磁屏蔽能力準確完整地展現出來。當我們再回首極為常見的NICS研究方法，即只是計算幾個點上的NICS數據，這和本文的研究方式相比無疑顯得太狹隘了。

8 卟啉

這是最后一個例子。利用examples\ICSS\里的porphyrin.gjf作ICSS_ZZ圖，等值面數值設為25，結果如下所示

ICSS_ZZ顯示體系中間一大塊區域對Z方向磁場有顯著的屏蔽，體現了卟啉有整體的大共軛效應故而有顯著的芳香性。卟啉有四個五元環，從圖中可以看到氮上不帶氫的兩個五元環有一大半露在外頭，而氮上帶著氫的兩個五元環上方的等值則是和卟啉中間區域等值面直接相連的。這說明不帶氫的五元環的外側原子不是整體共軛路徑的組成部分，而帶著氫的五元環的外側原子必定是整體共軛路徑的一部分。這一點也可以利用Multiwfn繪制分子平面上方1.2Bohr處的LOL-pi圖形來考察，LOL-pi可以清晰地把pi電子離域路徑展現出來。繪制過程是
porphyrin.fch  //此文件未提供，請自行生成。假設分子平面在YZ平面上
100
22    //設定所有sigma軌道占據數為0以去除其對LOL的貢獻
0
2
0   //回主菜單
4
10
1
[Enter]
3
1.2
然后再按照前文的過程讓原子標簽顯示出來，調整色彩刻度為0~0.66，得到下圖（白色代表超過0.66處）

從圖中可以看出一圈鮮明的pi電子共軛路徑，經過含氫的五元環時傾向于走外側，而經過不含氫的五元環時傾向于走內側。之所以這么說是因為在這樣的路徑上LOL-pi值一直都比較大。假設隨意選取一條路徑，其中經過了LOL-pi很小的區域，那么電子在這里就不容易離域過去，這樣的路徑就肯定不是理想的共軛路徑。雖然研究這類問題人們更多地使用ELF-pi，但至少對于此例LOL-pi比起ELF-pi的圖像更為清晰。LOL和ELF雖然出發點不同，但其實體現的物理意義是相仿佛的。如果繪制電流密度圖，也可以得到如上的分析結論。關于Multiwfn分析pi電子更多的介紹看《在Multiwfn中單獨考察pi電子結構特征》（http://www.shanxitv.org/432）。

9 總結&其它

本文主要介紹了ICSS，同時還結合磁屏蔽值在直線和平面上的分布以及LOL討論了一些體系的芳香性問題。這些研究方法都非常實用而且直觀形象，通過Multiwfn繪制起來也十分方便，完全可以作為像NICS一樣的研究芳香性的常規手段，而且比NICS明顯深入得多。實際上磁屏蔽信息和化學鍵還有著密切關聯，特別是筆者發現數值很大的ICSS等值面和ELF/LOL、AICD等值面有著一些內在關聯和相似性，筆者以后若有機會將會另撰文討論此問題。

值得一提的是，筆者使用ICSS_ZZ方法充分考察了電子結構非常特殊的、具有雙芳香性的18碳環體系，圖像非常漂亮（如下所示），傳遞出的信息十分有價值。非常推薦大家仔細閱讀相應的論文Carbon, 165, 468 (2020)中的討論，這也可以作為ICSS的很好的例子進行引用。順帶一提，筆者對18碳環已開展了廣泛研究，發表了諸多成果，匯總見http://www.shanxitv.org/carbon_ring.html。

上文給出的平面圖的色彩都是以彩虹方式變化的，實際上完全可以在Multiwfn主功能4繪制平面圖之后通過相應選項將色彩刻度改成其它模式，并且加上等值線等效果，這樣會漂亮得多得多，上面的18碳環研究文章里的圖就是極好的范例。另外，如果想要得到更漂亮的等值面圖，建議按此文的做法結合VMD來實現：《在VMD里將cube文件瞬間繪制成效果極佳的等值面圖的方法》（http://www.shanxitv.org/483），上面的ICSS_ZZ等值面圖就是這么繪制的。

筆者在Chem. Eur. J., 28, e202103815 (2022)中利用ICSS_ZZ研究了18碳環衍生物C18-(CO)2、C18-(CO)4、C18-(CO)6的芳香性，通過ICSS圖非常直觀地證明了它們的芳香性并對它們的芳香性強弱進行了區分，是ICSS_ZZ很好的應用范例，建議閱讀和引用。筆者還專門在《深入揭示18碳環的重要衍生物C18-(CO)n的電子結構和光學特性》（http://www.shanxitv.org/640）里對此論文做了深入淺出的介紹，十分推薦一看。