用NWChem做SODFT在DFT計算中考慮旋軌耦合效應
1 前言
眾所周知，對重原子構成的體系，尤其是從第五周期開始（H,He算第一周期，下同），需要考慮相對論效應，否則結果可能嚴重錯誤。相對論效應分為標量和旋軌耦合兩部分。標量相對論效應好考慮，用相對論贗勢時即可等效體現，多數主流量化程序如Gaussian也支持標量相對論哈密頓，結合專為標量相對論計算優化的全電子基組即可（見《在贗勢下做波函數分析的一些說明》http://www.shanxitv.org/156）。
比較棘手的是旋軌耦合的考慮，本文只說單個態的旋軌耦合問題。最直接的考慮方式是做二分量相對論計算，但比較耗時，且支持它的ADF、Turbomole都是收費的，而免費的Dirac等程序則比較小眾且功能局限性大。Gaussian里有個DKHSO方法可以在單點計算時考慮旋軌耦合效應，但是只能用于HF，因此沒什么實用性。一種常用的又便宜又好的考慮旋軌耦合的做法是使用相對論贗勢，但不僅要包含我們平時常用的標量勢，還同時要考慮旋軌勢，這種方式做DFT稱為Spin-orbit DFT (SODFT)，Molpro、Turbomole、NWChem都支持。本文介紹使用其中免費的NWChem程序做SODFT的方法，并會以TlCl和HgI的解離能計算做為例子。在NWChem中，SODFT有解析梯度，因此能夠容易地優化分子結構，不過Hessian只有數值的（只考慮標量勢的話倒是有解析Hessian）。NWChem的SODFT還一個不足之處是沒法利用對稱性加速計算。
NWChem的安裝方法看此文：《NWChem 6.6編譯方法》（http://www.shanxitv.org/270）。本文使用NWChem 6.6版。
2 SODFT計算時贗勢和贗勢基組的選擇
贗勢包括標量勢、旋軌勢和核極化勢。標量勢部分是我們最常使用的，所有相對論贗勢都提供了標量勢，但不是所有贗勢都給出了旋軌勢。給出了旋軌勢的贗勢比較常見的有Stuttgart贗勢和CRENB贗勢，對周期表覆蓋得都比較全面。后者很老，精度也不比前者好，所以本文將會用前者。Stuttgart贗勢也分多種，包括HF系列（贗勢不體現相對論效應），WB系列（基于Wood-Boring準相對論計算的數據搞的）和DF系列（基于Dirac-Fock全相對論計算的數據搞的），我們一般都用M版（擬合贗勢時考慮多個價電子），如MDF。
除贗勢外還要考慮贗勢基組。基于Stuttgart贗勢的常用的幾種選擇有：
1 Stuttgart贗勢的標配基組。缺點是對于大多數元素，不像下面幾種給出了不同檔次基組的定義，選擇余地小。
2 cc-pVnZ-PP系列：結合的是Stuttgart小核MDF贗勢
3 def2系列：從第五周期開始是贗勢基組，結合的是Stuttgart小核贗勢，對s、d、鑭系元素用的是MWB，對p元素用的是MDF
4 dhf系列：和def2都是Weigend主要搞的，只定義了五、六周期的s和d族元素，p族元素和def2定義相同。此系列分為dhf和dhf-2c兩種，前者主要用于做一分量計算（即常規的贗勢下的計算），把def2對s、d元素用的標量勢從MWB改為MDF并重新調整了贗勢基組定義，由此使結果略好了一點點；后者在前者基礎上增加了額外的p,d基函數，使得帶著旋軌勢做計算時結果更好
以上贗勢基組中，做SODFT時能用dhf-2c就用dhf-2c，畢竟是給此目的專門優化的。對于本文要考察的元素，Tl、Hg、I都可以用dhf-2c，我們用其中TZVP級別的就夠了，即dhf-TZVP-2c。雖然也有更好的dhf-QZVP-2c，但對于DFT計算升到QZ級別意義并不大。對于Cl，由于相對比較輕，其引起的旋軌耦合效應不會給結果帶來明顯影響，所以我們不對它用贗勢，而對它用def2-TZVP全電子基組。
3 計算TlCl
本文我們要嘗試重復下面圖中的TlCl和HgI的數據，數據來自Theor. Chem. Acc., 130, 633-644 (2011)。圖中BDE是鍵解離能，單位是kcal/mol。鍵解離能的計算是用解離后的原子的焓減去分子的焓。圖中BDE下面的列是文中只考慮標量相對論的計算結果，BDE+SO是考慮旋軌耦合后的。
我們先算TlCl。由于文中用的是PBE0，所以我們也用PBE0，此泛函對這類體系是比較合適的。我們用的鍵長直接取實驗值，大量雙原子分子的實驗平衡距離可以在NIST網站http://webbook.nist.gov/chemistry/上很容易地查到，注意勾選Constants of diatomic molecules，之后數據頁面里r_e就是平衡距離了，對此分子是2.484826埃。
首先我們得獲取贗勢和贗勢基組定義。def2和dhf系列在Turbomole基組庫http://www.cosmologic-services.de/basis-sets/basissets.php里是最全的（EMSL上的不僅不全，而且沒有旋軌勢）。在此網站里面選上Tl，然后選dhf-TZVP-2c，輸出格式就用默認的Turbomole（沒直接提供NWChem的比較可惜）。給出的數據包括基組部分、標量勢部分和旋軌勢部分。
PS1：值得一提的是，當選擇dhf-TZVP-2c的時候，對于非第五、六周期的元素，自動給出的是def2對應檔次基組的定義
PS2：此網站輸出格式里只有選Turbomole或Molpro這兩種支持旋軌勢的程序時輸出內容才有旋軌勢，我們這里不選Molpro是因為沒有Turbomole的格式那么容易手動改成NWChem的格式
然后我們把Turbomole格式的基組、標量勢和旋軌勢的定義都改成NWChem的格式，不知道NWChem的格式怎么寫就RTFM。一個較大區別是Turbomole里把收縮系數放在第一列，指數放在最后一列，和NWChem的順序是反過來的。最終TlCl的SODFT任務的輸入文件如下
start
PRINT low
geometry noautosym
Tl 0.00000 0.00000 0.00000
Cl 0.00000 0.00000 2.484826
end
basis "ao basis" spherical
Cl library def2-TZVP
Tl s
729.65038145 0.13672829828E-03
46.665548707 0.60443439951E-02
20.970448726 -0.20022066697
14.149588677 0.40801678488
Tl s
20.730134285 -0.71861135918E-01
6.1527631309 0.98057508445
Tl s
1.5757324480 1.0000000000
Tl s
0.74980169458 1.0000000000
Tl s
0.19536816194 1.0000000000
Tl s
0.70878767298E-01 1.0000000000
Tl p
15.383852616 0.61717949180
14.814929544 -0.72859235151
6.7261253658 0.40438195364
Tl p
1.9626182149 0.43157661160
1.0331857851 0.39230403853
0.53837445996 0.14007406820
Tl p
0.24446611676 1.0000000000
Tl p
0.90785377260E-01 1.0000000000
Tl p
0.33401321498E-01 1.0000000000
Tl d
57.606819928 0.16054811114E-03
9.7368866667 0.24456562496E-01
6.9256201679 -0.69914775031E-01
2.1396230731 0.19496269490
1.0836187110 0.29731629705
0.52356298209 0.23728708102
Tl d
0.23842309375 1.0000000000
Tl d
0.95000000000E-01 1.0000000000
Tl f
0.28435 1.0000000000
Tl f
0.95 1.0000000000
Tl p
32.0 -0.0005
16.0 -1.0
Tl p
8.0 -1.0125
4.0 0.7450
Tl p
2.0 1.0
Tl d
3.0 1.0
END
ECP
Tl nelec 60
Tl ul
2 1.00000000 0.00000000
Tl s
2 12.16780500 281.28466300
2 8.29490900 62.43425100
Tl p
2 7.15149200 4.63340800
2 5.17286500 9.34175600
2 9.89107200 72.29925300
2 9.00339100 144.55803700
Tl d
2 7.13021800 35.94303900
2 6.92690600 53.90959300
2 5.41757000 10.38193900
2 5.13868100 15.58382200
Tl f
2 5.62639900 15.82548800
2 5.54895200 21.10402100
2 2.87494600 2.91512700
2 2.82145100 3.89690300
Tl g
2 6.67905700 -7.49453400
2 6.70683500 -9.54057500
2 7.20928400 -7.79799200
2 7.07096400 -9.25952400
END
SO
Tl p
2 7.151492 -9.266817
2 5.172865 9.341756
2 9.891072 -144.598506
2 9.003391 144.558037
Tl d
2 7.130218 -35.943039
2 6.926906 35.939729
2 5.417570 -10.381939
2 5.138681 10.389215
Tl f
2 5.626399 -10.550326
2 5.548952 10.552010
2 2.874946 -1.943418
2 2.821451 1.948451
Tl g
2 6.679057 3.747267
2 6.706835 -3.816230
2 7.209284 3.898996
2 7.070964 -3.703809
END
dft
xc pbe0
end
task sodft
輸入文件還是很容易理解的，照常定義基組部分，ECP段落指定標量勢部分，SO字段指定旋軌勢部分，最后task里寫上SODFT就行了。noautosym寫不寫無所謂，由于SODFT不支持對稱性，不寫它時程序也不會利用對稱性來加速。
然后我們運行此任務，輸出文件里這部分即是SODFT能量：
Total DFT energy = -632.863447528456
在筆者的Intel 4核i7-2630QM機子下，做這個計算耗時是20.8s，若只考慮標量勢是10.8s，可見考慮旋軌勢還是會增加不少耗時，不過一般完全承受得起，總比做二分量相對論計算來考慮旋軌耦合效應便宜得多。
之后我們把TlCl的文件分別改寫為Tl和Cl的文件，二者都是二重態，所以在dft字段里面要寫上mult 2。算出的結果是Tl=-172.724034492809，Cl=-459.997396831618。因此基于單點能算的鍵解離能是627.51*(-459.997396831618-172.724034492809+632.863447528456)=89.116 kcal/mol。
到此還沒完，我們還得計算焓的熱校正量對解離能的修正。我們就不用SODFT算頻率了，就直接在Gaussian里用# PBE1PBE/def2TZVP opt freq scale=0.982關鍵詞對TlCl優化和做振動分析就完了，這里0.982是ZPE校正因子（為了省事我們不單獨考慮升溫對焓影響那部分的校正因子）。實際上我們查不到PBE0/def2-TZVP級別的校正因子，但由于J. Phys. Chem. A, 111, 11683 (2007)中給出PBE0結合6-311+G(2df,p)這樣不錯級別基組的校正因子是0.9824，所以有理由認為0.982對于PBE0結合def2-TZVP是基本合適的。計算后得到的焓的熱校正量輸出為
Thermal correction to Enthalpy= 0.004388
對于Tl和Cl沒必要做振動分析，因為單個原子只有平動對焓有貢獻，且是精確已知的，即5/2*RT，常溫下對應0.002360 Hartree。所以常溫下焓的熱校正量對TlCl鍵解離能的貢獻可以計算為627.51*(0.002360+0.002360-0.004388)=0.208 kcal/mol。
最后，得到TlCl的常溫下鍵解離能是89.116+0.208=89.324 kcal/mol，這和前面貼出來的實驗值88.1 kcal/mol相符極好！（對普通DFT泛函來說能達到1kcal/mol程度誤差著實難能可貴）
4 計算HgI
HgI是二重態，解離后會生成單重態的Hg和二重態的I。HgI的鍵長在NIST上查不到，雖然某些文獻里估計能查到，但這里不深究，就直接用PBE0/def2-TZVP級別來優化鍵長，之后再做單點計算。優化后鍵長是2.75663埃。（當然，若要求更精確，可以直接用SODFT優化。筆者也試了下，i7-2630QM上耗時336s，結果是2.76379埃。可見對此體系是否考慮旋軌耦合對結構影響很小，對TlCl測試過也是如此）。
我們使用上一節的方式，獲取Hg和I的dhf-TZVP-2c贗勢基組定義以及對應的標量勢和旋軌勢，寫成HgI的PBE0的SODFT任務的單點輸入文件，如下所示：
start
PRINT low
geometry noautosym
Hg 0.00000 0.00000 0.00000
I 0.00000 0.00000 2.75663
end
basis "ao basis" spherical
I s
5899.5791533 0.24188269271E-03
898.54238765 0.15474041742E-02
200.37237912 0.42836684457E-02
31.418053840 -0.39417936275E-01
15.645987838 0.96086691992
I s
11.815741857 0.57815778954
6.4614458287 0.37374293124
I s
2.3838067579 1.0000000000
I s
1.1712089662 1.0000000000
I s
0.32115875757 1.0000000000
I s
0.12387919364 1.0000000000
I p
185.43362455 0.83127824000E-03
20.091408146 0.63991653000E-01
9.7577022390 -0.27791138000
I p
13.068307912 -0.49793790382E-01
3.5818714205 0.38212490511
2.0282441852 0.70447564804
1.0181492146 0.33781067803
I p
0.46673773115 1.0000000000
I p
0.19242597960 1.0000000000
I p
0.74508878495E-01 1.0000000000
I d
124.20341062 0.65671747209E-03
34.587311801 0.51648185674E-02
12.767328064 -0.19881371307E-01
4.7745100133 0.21386794109
2.4582209028 0.43405444707
1.1923708147 0.37850637882
I d
0.52883971906 1.0000000000
I d
0.17008164307 1.0000000000
I f
0.44141808 1.00000000
I f
2.18 1.0000000000
I p
12.0 -0.3
6.0 1.0
I d
2.0 1.0
Hg s
48.013786990 0.58613168385E-02
21.239875095 -0.17590367988
15.876100879 0.35780354753
Hg s
5.4837607070 1.0000000000
Hg s
1.5480592128 1.0000000000
Hg s
0.72425230437 1.0000000000
Hg s
0.16369906863 1.0000000000
Hg s
0.57211615334E-01 1.0000000000
Hg p
23.291760168 -0.83564430982E-02
13.028969569 0.83703058587E-01
6.5100040792 -0.31023833705
Hg p
1.8167935815 1.0000000000
Hg p
0.90079391013 1.0000000000
Hg p
0.41304090835 1.0000000000
Hg p
0.11845879331 1.0000000000
Hg p
0.36084184656E-01 1.0000000000
Hg d
15.176302343 0.60654575178E-02
6.7004896493 -0.59880306300E-01
1.9144256118 0.31411145584
0.88641552102 0.45081091161
Hg d
0.38364767725 1.0000000000
Hg d
0.14936891490 1.0000000000
Hg f
0.79569 1.0
Hg p
19.3582 0.0635
8.7992 -0.4607
Hg p
3.9996 0.8853
Hg d
2.7000000000 1.0000000000
END
ECP
I nelec 28
I ul
2 1.000000 0.000000
I s
2 40.033376 49.989649
2 17.300576 281.006556
2 8.851720 61.416739
I p
2 15.720141 67.416239
2 15.208222 134.807696
2 8.294186 14.566548
2 7.753949 28.968422
I d
2 13.817751 35.538756
2 13.587805 53.339759
2 6.947630 9.716466
2 6.960099 14.977500
I f
2 18.522950 -20.176618
2 18.251035 -26.088077
2 7.557901 -0.220434
2 7.597404 -0.221646
Hg nelec 60
Hg ul
2 1.0000000 0.0000000
Hg s
2 12.4130710 275.7747970
2 6.8979130 49.2678980
Hg p
2 11.3103200 80.5069840
2 10.2107730 161.0348240
2 5.9398040 9.0834160
2 5.0197550 18.3677730
Hg d
2 8.4078950 51.1372560
2 8.2140860 76.7074590
2 4.0126120 6.5618210
2 3.7953980 9.8180700
Hg f
2 3.2731060 9.4290010
2 3.2083210 12.4948560
Hg g
2 4.4852960 -6.3384140
2 4.5132000 -8.0998630
END
SO
I p
2 15.720141 -134.832478
2 15.208222 134.807696
2 8.294186 -29.133096
2 7.753949 28.968422
I d
2 13.817751 -35.538756
2 13.587805 35.559839
2 6.947630 -9.716466
2 6.960099 9.985000
I f
2 18.522950 13.451079
2 18.251035 -13.044039
2 7.557901 0.146956
2 7.597404 -0.110823
Hg p
2 11.31032000 -161.01396700
2 10.21077300 161.03482400
2 5.93980400 -18.16683200
2 5.01975500 18.36777300
Hg d
2 8.40789500 -51.13725600
2 8.21408600 51.13830600
2 4.01261200 -6.56182100
2 3.79539800 6.54538000
Hg f
2 3.27310600 -6.28600100
2 3.20832100 6.24742800
Hg g
2 4.48529600 3.16920700
2 4.51320000 -3.23994500
END
dft
mult 2
xc pbe0
end
task sodft
之后稍作修改，得到Hg和I的PBE0-SODFT任務輸入文件。計算后將能量求差，得到基于電子能量的BDE：
627.51*(-295.746562477964-153.684090186717+449.448365018831)=11.115 kcal/mol
如前面的圖所示，HgI的解離能的三個不同來源的實驗值為7.8、8.1、8.9 kcal/mol，當前算的和哪個都想差甚遠，因此PBE0對此體系不夠給力。對這種情況，我們可以在DFT級別計算出旋軌耦合對BDE的校正，然后加到CCSD(T)/def2-QZVP這樣只考慮標量相對論的高計算級別得到的解離能上。為此，我們把上面用的HgI、Hg、I輸入文件里的sodft改成dft（之后SO段落可以刪掉，不刪也無所謂），這樣計算就只利用標量勢部分了，旋軌耦合效應就不考慮了。結果為
627.51*(-295.697077663054-153.500765444841+449.222674994391)=15.582 kcal/mol
因此，旋軌耦合對BDE的校正量為11.115-15.582=-4.467 kcal/mol。
在CCSD(T)/def2QZVP級別算的HgI解離過程的電子能量變化：
627.51*(-153.2231386-297.2662532+450.5085777)=12.039 kcal/mol
用Gaussian在PBE0/def2-TZVP下做振動分析，得到常溫下焓的熱校正量對HgI鍵解離能的貢獻：
627.51*(0.002360+0.002360-0.004277)=0.278 kcal/mol
put together，得到HgI的常溫下的BDE：
12.039-4.467+0.278 = 7.85 kcal/mol
可見結果和實驗符合得極好！
5 SODFT方式計算BDE總結
前面HgI的例子很有代表性，說明想準確計算重原子的鍵解離能，最簡單省事的過程就是：
1 用較合適的泛函，在def2-TZVP或相似級別的結合相對論贗勢的贗勢基組（如cc-pVTZ-PP，但不劃算）下進行優化和做振動分析，得到焓的熱校正量對BDE的影響。注意應考慮ZPE校正因子
2 用較合適的泛函，在1的結構下用dhf-TZVP-2c贗勢基組與配套的MDF標量勢+旋軌勢做SODFT計算，再去掉旋軌勢做普通DFT計算，利用差值得到旋軌耦合對BDE的影響
3 基于1的結構，在盡可能高的級別下計算考慮了標量相對論效應的單點，由此得到基于電子能量算的BDE。體系很小時計算級別一般用CCSD(T)結合def2-QZVP或cc-pVQZ-PP贗勢基組，贗勢只考慮標量勢部分。可以把結果再外推到CBS。
最后把1、2、3相加，就得到了特定溫度下的鍵解離能。
上面說的是體系只含第五、六周期s,d,p族元素的情況。對于體系存在其它元素時，筆者認為應當這樣考慮：
對前三周期的元素，就用全電子基組就行，不需要考慮相對論效應。
對第四周期元素，考慮相對論效應有益。分為兩種情況：
(a)Cu之前的：Stuttgart贗勢沒有給出這些元素的旋軌勢，cc-pVnZ-PP也沒有給出這些元素的定義，def2對這些還是全電子基組。要考慮相對論效應就用Stuttgart贗勢結合其標配贗勢基組就完了。非要也考慮旋軌耦合那就用CRENB標量勢+旋軌勢結合CRENBL贗勢基組吧。
(b)Cu~Kr：用cc-pVnZ-PP贗勢基組，結合對應的MDF標量勢和旋軌勢。從EMSL上取cc-pVnZ-PP定義的時候贗勢部分只有標量勢部分，旋軌勢得從Stuttgart贗勢的官網http://www.tc.uni-koeln.de/PP/clickpse.en.html上自行取，輸出格式必須選molpro，要選小核的MDF版本（比如Br要選ECP10MDF）。輸出的贗勢里前幾行是標量勢（和EMSL上cc-pVnZ-PP直接帶的一致），后幾行是旋軌勢。
對于鑭系、錒系及一些第七周期的元素，就用Stuttgart贗勢官網上小核標量勢+旋軌勢結合標配的贗勢基組。情況較復雜。有的有MWB有的有MDF有的都有（優先用后者），有的有旋軌勢有的沒有旋軌勢（沒有的話就沒轍了）。
另外，上面說的第1步的優化+振動分析，如果求準確，可直接用SODFT來做，但是由于SODFT只有半數值頻率，體系稍大一點就會極為吃力。對于第3步高精度的考慮標量相對論的單點，不是必須用相對論贗勢，用高質量全電子基組并使用較好的標量相對論哈密頓（LUT-IOTC、DKH3等）也可以，只不過耗時會高不少而且精度也未必更好。