使用Multiwfn巨方便地繪制二維NICS平面圖考察芳香性

使用Multiwfn巨方便地繪制二維NICS平面圖考察芳香性

文/Sobereva@北京科音   2023-Aug-9

1 前言

NICS在衡量芳香性上用得很普遍，在《衡量芳香性的方法以及在Multiwfn中的計算》（http://www.shanxitv.org/176）里有充分介紹，不了解者務必先閱讀。強大的波函數分析程序Multiwfn具有很好用的繪制一維NICS圖的功能，見《使用Multiwfn繪制一維NICS曲線并通過積分衡量芳香性》（http://www.shanxitv.org/681）；此外，Multiwfn支持的ICSS分析在本質上是考察三維空間的NICS（ICSS和NICS僅相差正負號），見《通過Multiwfn繪制等化學屏蔽表面(ICSS)研究芳香性》（http://www.shanxitv.org/216）。雖然如此文所示例的，利用Multiwfn的ICSS模塊結合Gaussian程序產生三維ICSS格點數據后，可以利用Multiwfn以格點數據插值方式得到任意平面上的二維ICSS/NICS圖，但如果你的目的僅僅是繪制二維NICS圖，那么為此很耗時地計算三維ICSS格點數據明顯是不劃算的。本文就演示Multiwfn里新加入的專門繪制二維NICS圖的功能，僅需要借助Gaussian計算分布在作圖平面上的格點即可，耗時比計算三維ICSS格點數據低一、兩個數量級。

讀者請務必使用2023-Aug-7及以后更新的Multiwfn版本，可以在官網http://www.shanxitv.org/multiwfn免費下載。對Multiwfn不了解者建議閱讀《Multiwfn入門tips》（http://www.shanxitv.org/167）和《Multiwfn FAQ》（http://www.shanxitv.org/452）。本文用的Gaussian是G16 C.02。

2 例1：繪制暈苯平面上方1埃處NICS_ZZ平面圖

Multiwfn程序包里的examples\NICS_scan\coronene.pdb是B3LYP/6-31G*級別優化的暈苯的結構文件，如下所示，分子在Z=0的XY平面上。此例繪制分子上方1埃處的平面的NICS_ZZ的填色圖。注意本文里的NICS_ZZ值是指垂直于被考察的環平面的磁屏蔽張量分量的負值，而非特指笛卡爾Z軸方向的分量，后同。

啟動Multiwfn，然后輸入
examples\NICS_scan\coronene.pdb
25   //離域性與芳香性分析
14   //繪制NICS二維平面圖
1   //填色圖
[回車]  //用默認的格點數，即兩個方向都是100個點，因此要計算100*100=10000個Bq
0   //設置延展距離
1  //1 Bohr。對當前體系用比默認小得多的延展距離便足以，免得圖中有過多的體系外圍的不感興趣的區域。不知道什么叫延展距離的話看《Multiwfn FAQ》（http://www.shanxitv.org/452）的4.2節
1   //XY平面
1a   //Z=1埃
1   //產生Gaussian的NICS二維掃描的輸入文件
examples\NICS_scan\template_NMR.gjf   //這是Gaussian做NMR計算的模板文件，里面原子坐標部分用[geometry]代替，會被自動替換

現在當前目錄下產生了NICS_2D.gjf。可以根據實際情況進行修改，比如修改凈電荷、自旋多重度、基組等等，不懂的地方別亂改。這里把基組設為6-31G*，其它的不變。之后用Gaussian計算之，在2*7R32機子96核并行情況下花了1m50s算完。相應的輸入和輸出文件是examples\NICS_scan\目錄下的coronene_NICS_2D.gjf和coronene_NICS_2D.out。

然后在Multiwfn界面里輸入2選擇載入Gaussian輸出文件，然后輸入其路徑examples\NICS_scan\coronene_NICS_2D.out。程序問你要獲得哪種NICS，可以選各項同性值、各向異性值、平行于笛卡爾X或Y或Z方向的分量值、順著特定矢量的分量值。這里選擇5，即平行于Z方向的值。由于當前體系平行于XY平面，因此這么選得到的對應于一般意義的NICS_ZZ。

然后圖像馬上就蹦出來了，同時文本窗口顯示了平面數據最負和最正數值，如下所示，可見作圖平面上NICS_ZZ最負值為-43.3 ppm。可以根據此值恰當地設置色彩刻度上下限范圍
The minimum of data:  -43.3034000000000
The maximum of data:   11.3487000000000

關閉圖像，輸入以下命令修改作圖效果
4   //顯示原子標簽
1  //紅色
8    //顯示化學鍵
14  //棕色
17   //設置顯示標簽的距離閾值
5   //距離作圖平面5 Bohr以內的原子標簽才顯示出來
y   //更遠的原子用細體字顯示標簽
1   //修改色彩刻度范圍
-45,45
-8   //坐標軸改為以埃為單位
-2   //修改坐標軸刻度
2,2,10
2   //顯示出等值線
3   //修改等值線設置
8   //按等差數列生成等值線數值
-50,5,21  //起始值，步長，步數
y  //替換原有的等值線數值。之后等值線數值為-50,-45,-40...略...40,45,50
1   //保存并返回
-1  //重新作圖

現在看到下圖

此圖顏色越深說明NICS_ZZ越負，即對垂直于體系平面方向的磁場屏蔽效應越強。此圖中外圍六元環的顏色明顯比中間的六元環更深，體現出外圍六元環有更強的芳香性。用其它芳香性指標也能體現這一點，比如用Multiwfn計算多中心鍵級的話，會發現中間的環的六中心鍵級是0.022，而外圍的是0.034，也展現出外圍的六元環的芳香性相對更強。

在后處理菜單還可以選-7用來給平面數據乘上一個因子，如果設-1，則數據的符號就會反轉，之后繪制的相當于ICSS平面圖。

3 例2：繪制三苯胺的苯環上方1埃處的NICS_ZZ平面圖

在B3LYP/6-31G*下優化的三苯胺的結構文件是examples\NICS_scan\N(phenyl)3.pdb，如下所示，可見環都是傾斜著的。此例繪制23,24,26,30,28,25原子組成的環（已高亮顯示）的上方1埃平面處的NICS_ZZ平面圖，取垂直于環方向的磁屏蔽張量分量。

啟動Multiwfn，然后輸入
examples\NICS_scan\N(phenyl)3.pdb
25   //離域性與芳香性分析
14   //繪制NICS二維平面圖
1   //填色圖
[回車]  //用默認的格點數100*100
8   //繪制某些原子構成的擬合平面的上方或下方的平面
23,24,26,30,28,25   //環中的原子
注意此時在屏幕上顯示了這6個原子擬合的平面的單位法矢量：
The unit normal vector is    0.33076524    0.57300118    0.74984265

接著在Multiwfn里輸入
1   //繪制的是與擬合平面相平行而在它上方1埃的平面。如果輸入負值代表在下方
6   //繪圖平面的邊長為6埃

此時在Multiwfn窗口中顯示了在VMD程序里繪制作圖平面區域的命令：
draw triangle {   2.495   2.581  -1.739} {  -1.211  -0.235   2.047} {   6.776  -1.451  -0.547}
draw triangle {  -1.211  -0.235   2.047} {   6.776  -1.451  -0.547} {   3.070  -4.266   3.239}
draw material Transparent
如果將N(phenyl)3.pdb載入VMD，然后在文本窗口里輸入以上命令，恰當修改顯示方式后就可以看到下圖，藍色透明區域對應實際作圖平面，由此可以檢驗作圖平面設定得是否合理，當前顯然是合理的。被選擇的六個原子的幾何中心在作圖平面上的投影位置對應于作圖平面的正中心。

接著在Multiwfn里輸入
1   //產生Gaussian的NICS二維掃描的輸入文件
examples\NICS_scan\template_NMR.gjf   //Gaussian做NMR計算的模板文件

當前目錄下產生了NICS_2D.gjf。將基組改為6-31G*，然后讓Gaussian運行，在2*7R32 96核機子上2m5s算完。相應的輸入和輸出文件分別是examples\NICS_scan\目錄下的N(phenyl)3_NICS_2D.gjf和N(phenyl)3_NICS_2D.out。

接著在Multiwfn里輸入
2  //載入Gaussian的NICS掃描的輸出文件
examples\NICS_scan\N(phenyl)3_NICS_2D.out
0   //對各個點，取NICS在特定矢量上的投影
0.33076524    0.57300118    0.74984265   //之前Multiwfn顯示的擬合平面的單位法矢量

之后可見當前的作圖平面上的數據的最負和最正值：
The minimum of data:  -34.0268462149550
The maximum of data:   4.05938599881831

關閉圖像，然后輸入
4   //顯示原子標簽
1  //紅色
8    //顯示化學鍵
14  //棕色
17   //設置顯示標簽的距離閾值
5   //距離作圖平面5 Bohr以內的原子標簽才顯示出來
y   //更遠的原子用細體字顯示標簽
1   //修改色彩刻度范圍
-34,0
-8   //坐標軸改為以埃為單位
-2   //修改坐標軸刻度
1,1,5
19   //修改色彩變化方式
19   //黃-橙-黑
-1  //重新作圖

由圖可見在環平面上方的六元環內區域的NICS_ZZ很負，體現出這個環的顯著芳香性。更具體來說，環中心正上方不是最負的，靠近C-C鍵的內側區域是稍微更負的。這樣的NICS分布是普遍現象。

4 例3：繪制C16的垂直于環平面的NICS_ZZ平面圖

筆者之前對18碳環體系（http://www.shanxitv.org/carbon_ring.html）做過大量研究，其中包括其芳香性。其類似物16碳環（C16）最近也被合成了出來。16碳環的in-plane和out-of-plane兩類pi軌道都有16個電子，滿足Huckel反芳香性規則，因此可以預期C16具有雙反芳香性特征。這里使用Multiwfn對垂直于它環平面的截面繪制NICS_ZZ圖，以考察分子由內到外不同區域對垂直于環方向磁場的屏蔽效果。wB97XD/def2-TZVP優化過的C16的結構文件是examples\NICS_scan\C16.pdb，如下所示，可見繪制X=0的YZ平面或者Y=0的XZ平面都可以滿足當前目的。

啟動Multiwfn，然后輸入
examples\NICS_scan\C16.pdb
25   //離域性與芳香性分析
14   //繪制NICS二維平面圖
1   //填色圖
[回車]  //用默認的格點數，即兩個方向都是100個點
0   //設置延展距離
1  //9 Bohr。把延展距離設得比默認更大以展現更廣闊區域的磁屏蔽情況
2   //XY平面
0   //Y=0
1   //產生Gaussian的NICS二維掃描的輸入文件
examples\NICS_scan\template_NMR.gjf   //Gaussian做NMR計算的模板文件

之后把計算級別改成wB97XD/def2-TZVP，因為筆者在Carbon, 165, 468 (2020)里證明了這個級別可以描述好碳環的結構，而常用的B3LYP則絕對不能用，《談談量子化學研究中什么時候用B3LYP泛函優化幾何結構是適當的》（http://www.shanxitv.org/557）中也提到了這點。而且《我對一篇存在大量錯誤的J.Mol.Model.期刊上的18碳環研究文章的comment》（http://www.shanxitv.org/584）里也說了6-31G*無法合理描述碳環類體系的電子結構，因此這里用的基組質量比前面的例子更好。

此例的NICS掃描的輸入輸出文件是examples\NICS_scan\目錄下的C16_NICS_2D.gjf和C16_NICS_2D.out。在Multiwfn當前的界面里選2并載入此輸出文件，關閉圖像后用以下命令調節作圖設置
4   //顯示原子標簽
12  //深綠色
8   //顯示化學鍵
14   //棕色
17   //設置顯示標簽的距離閾值
50   //距離作圖平面50 Bohr以內的原子標簽才顯示出來
y   //更遠的原子用細體字顯示標簽
1   //修改色彩刻度范圍
-50,50
-8   //坐標軸改為以埃為單位
-2   //修改坐標軸刻度
2,2,10
2   //顯示出等值線
3   //修改等值線設置
8   //按等差數列生成等值線數值
-80,5,33  //起始值，步長，步數
y  //替換原有的等值線數值
1   //保存并返回
19  //設置色彩變化
8    //藍-白-紅
-1  //重新作圖

看到的圖像如下。可見在環中央及其上下方區域NICS_ZZ都為非常顯著的正值，Multiwfn顯示的平面上最正值都達到了72.3 ppm，說明C16有極強的反芳香性。

上圖可以跟《通過Multiwfn繪制等化學屏蔽表面(ICSS)研究芳香性》（http://www.shanxitv.org/216）里列舉的Carbon, 165, 468 (2020)文中的18碳環截面的ICSS二維圖進行對比，可見特征一樣但符號相反（再次注意ICSS和NICS相差正負號），也即C16和C18的芳香性特征正好是相反的。

5 總結

本文介紹了Multiwfn的非常簡單易用的繪制二維NICS平面圖的功能，定義作圖平面的方式非常靈活，作圖設置豐富，效果極好，在像樣的機子上用Gaussian計算繪制這種圖所需的磁屏蔽張量數據的耗時一般也比較低（在較好的服務器上算不很大的體系屬于立等可取的耗時），非常推薦在研究芳香性問題時恰當地使用。

使用本文的功能繪制NICS平面圖用于論文的情況請按照Multiwfn啟動時的提示恰當引用Multiwfn的原文。