NMR化學位移和耦合常數計算簡介

本文簡單介紹一下怎麼用高斯的DFT計算NMR化學位移和耦合常數，以及如何在GaussView裡觀看計算生成的圖譜。計算和可視化所用的軟體分別是G16 A.03, GV 6.0.16。

用量子化學方法計算NMR化學位移現在已經有了常用套路，DFT泛函一般選取B97-2（註：高斯關鍵詞寫作B972），基組常用pcSseg-1或pcSseg-2，更換其他泛函（如B3LYP, M062X）或基組（如6-311G**），結果精度在統計意義上是不會比上述常用搭配要好的。當然，如果體系很小或土豪機器配置很高，用得起cc-pVQZ基組，那自然是用QZ更好。

 

與實驗上的1H-NMR化學位移需要用到標準物質四甲基矽烷（TMS）一樣，理論計算上也需要分別計算TMS與目標分子中氫的核磁屏蔽數值，然後相減就是目標分子中氫的核磁化學位移。下面展示的是TMS分子的高斯輸入文件

%chk=TMS_def2TZVP_opt.chk%mem=16GB%nprocshared=28[空行]geometry opt of TMS at Td point group [標題行，可瞎寫][空行]0 1Si    0.00000000   0.00000000   0.00000000C     0.00000000   0.00000000   1.94000006H    -0.87364980  -0.50440195   2.29666647H     0.87364980  -0.50440195   2.29666647H     0.00000000   1.00880390   2.29666647C    -0.00000000  -1.82904960  -0.64666669H    -0.87364980  -2.33345189  -0.29000077H    -0.00000000  -1.82904994  -1.71666493H     0.87364980  -2.33345189  -0.29000077C     1.58400341   0.91452480  -0.64666669H     1.58400372   0.91452497  -1.71666493H     1.58400372   1.92332887  -0.29000077H     2.45765352   0.41012302  -0.29000077C    -1.58400341   0.91452480  -0.64666669H    -1.58400372   1.92332887  -0.29000077H    -1.58400372   0.91452497  -1.71666493H    -2.45765352   0.41012302  -0.29000077[空行]--Link1--%chk=TMS_def2TZVP_opt.chk%mem=16GB%nprocshared=28[空行][自行去EMSL網站上拷貝pcSseg-2基組數據，貼於此]

  在GV裡構建TMS分子結構後，可以利用右鍵Tools -> Point Group工具，然後勾選Enable Point Group Symmetry，點擊右中區域Symmetrize，使結構符合Td點群對稱性，這樣所有氫是核磁化學等價的，算出來的核磁屏蔽數值一模一樣。
從上面的輸入文件可以看出，先用M062X/def2TZVP進行結構優化，後用B972/pcSseg-2算NMR化學位移。計算核磁的關鍵詞為NMR，相當於NMR=GIAO，表示用默認的GIAO方法計算。兩步任務所使用的泛函和基組不要求一致。這有兩點原因：（1）幾何結構的描述上B97-2泛函未必是最好的；（2）pcSseg-2基組挺大的（比cc-pVTZ還大），拿它算單點還行，做結構優化時間就比較長了。如果將來研究的體系比較大，無論是TMS還是目標分子，基組建議統一降為pcSseg-1。隱式溶劑化模型的溶劑用的是氯仿，相當於實驗上的氘代氯仿。若不用隱式溶劑化模型，氣相中的結果一般也還行。
算完以後用文本編輯器（如Notepad++）打開log（或out）文件，搜索Isotropic，會發現每個原子有一處，每處等號後面的小數就是對應原子的核磁屏蔽數值，氫、碳、矽全都有（見下圖），這裡氫的核磁屏蔽各向同性數值為31.6753。Linux下可以直接用命令grep 'H   Isotropic =' xx.log搜索、列印到屏幕上。
 
接下來我們以丙酮分子為例，計算其1H-NMR化學位移。丙酮分子的輸入文件與TMS類似，這裡就不貼出了。同樣，在構造丙酮分子之後可以選擇點群工具使結構對稱化為C2v點群；不做對稱化、直接優化也行。要注意的是，對稱化這一步並不能保證最後的優化結果被高斯識別為C2v點群，可能會有微小的偏差；嚴格的C2v點群需要使用內坐標定義分子結構。上面的TMS分子優化後仍保持Td點群是湊巧的。下圖展示的是從丙酮輸出文件grep出的結果
3  H    Isotropic =    29.3619   Anisotropy =     4.82534  H    Isotropic =    29.2212   Anisotropy =     4.20415  H    Isotropic =    29.8594   Anisotropy =     6.64857  H    Isotropic =    29.3618   Anisotropy =     4.82498  H    Isotropic =    29.2211   Anisotropy =     4.20449  H    Isotropic =    29.8594   Anisotropy =     6.6487

可以看出左側3個氫（3,4,5）與右側3個氫（7,8,9）有十分微小的差異，這是略微偏離C2v點群造成的，影響可以忽略。但是3,4,5三個氫的數值怎麼會不相等呢？實驗上只有一個甲基氫的峰。這是由於我們計算的只是一個靜態結構，室溫（至0 ℃範圍均可）時丙酮分子中的甲基是可以自由旋轉的，因此實驗測得的是大量不同構象（符合玻爾茲曼分布）下的統計平均結果。如果是凝聚相下大分子中化學位移的計算，是要通過跑MD從軌跡中選取幀數進行統計平均的；不過對於常見有機小分子，當精度要求不是很高的時候，我們只要簡單地把優化得到的結構中的三個氫平均一下再與TMS進行比較即可：


這相當於氫原子之間交換序號得到的local minima的構象平均。實驗上丙酮在氘代氯仿溶劑中的化學位移是多少呢？筆者不會做實驗，只好去找資料庫了，比如AIST的SDBS有機分子譜學資料庫
https://sdbs.db.aist.go.jp/sdbs/cgi-bin/cre_index.cgi
常見小分子各種譜的數據挺詳細的。同意網站搜索協議之後，搜索acetone



接下來選擇對應的細分類、溶劑、氫譜，可以看到實驗值為2.162 ppm。說明上述計算結果還是挺接近的。按照相應的流程，也可以計算13C-NMR碳譜、19F-NMR氟譜等等，只要注意對應的標準物分別是什麼就行了。不過非氫元素的化學位移值本身往往很大，因此計算與實驗相差2 ppm左右也是正常的，DFT結合隱式溶劑化模型不能苛求更高的精度了。
 
接下來介紹一下怎麼在GV裡看NMR圖譜（這純粹是一種直觀的圖示功能，有的同學直接把化學位移值標在相應的鍵線式原子旁也很直觀）。進入GV安裝目錄的data文件夾，找到nmr.data文件用文本編輯器打開（在修改前可以先備份之），可以看到裡面已經預存了一些標準物的NMR屏蔽值，但都是過時的方法或基組，不符合現今需要。我們可以在底下添加一些常用的標準物質


雙引號裡面是注釋，記錄是在什麼水平下算的，以防日後忘記。修改之後保存，重啟GV，打開丙酮的log文件，面板上選擇Results -> NMR，會彈出NMR譜圖對話框，依次選擇Element: H, Reference: TMS B972/pcSseg-2 ultrafineGIAO （對應剛添加的nmr.data文件內容）。


橫坐標表示NMR化學位移，縱坐標是簡併度，即表示有幾個相同的氫。上圖可以看出丙酮左右兩邊甲基的氫是對稱的；同時由於我們算的是靜態結構，同一個甲基上的三個氫不等價。如果想顯示等價的圖示效果，可以用平均值將log文件內6個氫的不同數值替換掉，再打開就會顯示等價了。GV的NMR譜圖對話框Plots菜單中還有導出圖片、導出數據及在對應原子上顯示化學位移等各種選項，這裡不就一一介紹了。
 
至於耦合常數的計算，只需要加上spinspin子關鍵詞，即NMR=spinspin。不過要注意的是，耦合常數計算時間是頻率計算時間的大概兩倍[1]，計算量較大。我們再試一個1,1-二氯乙烷的NMR化學位移和耦合常數計算，結果總結如下：
實驗值同樣從SDBS資料庫查得，耦合常數實驗值為6.0。我們在log文件內搜索Total nuclear spin-spin coupling J，找到矩陣元(6,2),(6,3), (6,4)，將三個矩陣元取平均算得5.6，這就是甲基氫與alpha位氫的耦合常數JAB，與實驗值吻合得還行。注1：筆者構建的1,1-二氯乙烷分子中2,3,4是甲基氫，6號是alpha位氫，讀者需根據自己計算時構建的序號來找對應的矩陣元。注2：截圖中D是一種Fortran雙精度浮點數的格式，相當於Python和C中的E，表示乘以10的多少次方的意思。
 
 
在結構優化之後，使用常見套路B972/pcSseg-2計算NMR化學位移和耦合常數，結果與實驗值吻合不錯。
 
PS1：標準物質的NMR屏蔽數值顯然只需計算一次即可，以後對同一計算級別的目標分子可以重複使用。

PS2：如果是量化新手，不妨練習一下如何找到丙酮中甲基轉動的過渡態（在優化好的結構上旋轉一個甲基60°即可作為過渡態初猜，筆者算的Gibbs自由能壘是0.07 kcal/mol）。
 
PS3：計算核磁時需要分別引用泛函B97-2（見高斯官網http://gaussian.com/dft）、基組pcSseg-2（在EMSL上打開基組數據時就有引用文獻），及此搭配計算NMR化學位移的文章（任一篇高水平文章，用於支持這麼做的理由）。
 
PS4：ChemDraw中也有預測氫譜和碳譜化學位移的功能，畫好鍵線式之後選中，菜單欄上點擊Structure -> predict 1H-NMR shifts，即點即出，對於定性、半定量地與實驗值對比，還是很方便的。
 
PS5：推薦一篇計算NMR化學位移的實例應用paper，DOI: 10.1021/acs.jctc.7b00380。
 
PS6：幾個常見譜圖資料庫https://mason.gmu.edu/~sslayden/Lab/spec-db.htm
 
參考資料
1. http://gaussian.com/nmr/
2. Sobereva. 《談談如何又好又快地計算NMR化學位移》
http://bbs.keinsci.com/thread-4434-1-1.html