在贗勢下做波函數分析的一些說明

 
經常有人問在贗勢下能不能做AIM分析、Multiwfn (http://sobereva.com/multiwfn)支持不支持贗勢等等問題。此文專門就此問題進行探討。
 

1 全電子相對論計算

對于重原子，從第四周期開始（K~Kr）相對論效應就已經體現出來了，第四周期之后就必須考慮相對論效應了。對于使用全電子基組時，一般通過DKH等方法做標量相對論計算，比如在Gaussian中用DKH2方法時只需寫上int=dkh2就可以了。

平時常用的基組的收縮方式并不是為了相對論計算而優化的，用于相對論計算并不適合，建議用相對論計算專用基組。比如(aug)-cc-pVnZ-DK系列比(aug)-cc-pVnZ系列明顯更適合DKH2計算。有些人也將已有的非相對論基組重新調整收縮方式使之適合相對論計算。完全未收縮的基組也可以直接用在相對論計算上，但是計算量往往過大。

適合相對論計算的全電子基組有這些：

• (aug)-cc-pVnZ-DK系列：(aug)-cc-pVnZ的專用于DKH計算的重收縮版本。主要適合后HF計算，由于高角動量函數多，廣義收縮程度又高，DFT計算用它的話虧。支持前六周期元素（K/Rb/Cs、Ca/Sr/Ba除外）和錒系的Th/Pa/U。
• def2系列重收縮版：def2系列基組從H到Xe都有全電子基組版本，但沒有對相對論計算而優化。ORCA的一伙人對這系列基組進行了重收縮使之適合DKH2、ZORA計算。這也是ORCA對于H~Xe的DKH/ZORA計算時標配的基組。很適合DFT計算。基組定義可以令ORCA程序用PrintBasis關鍵詞打印出來放到別的程序中使用。
• SARC（Segmented all-electron relativistically contracted）：ORCA程序的主要開發者提出的，是很適合DKH2、ZORA結合DFT計算的全電子基組，基組不太大結果也挺好。目前公開發表的只有鑭系、錒系、第三周期過渡金屬和Tl~Rn。很適合搭配def2系列重收縮版，這就可以涵蓋周期表絕大多數元素了。后來又提出了SARC2，比SARC精度更高，但覆蓋的元素目前還較少。
• ANO-RCC：元素覆蓋全面，很適合考慮電子相關作用的相對論計算。但是基組巨大，收縮度特別高，高角動量函數又多，對于HF/DFT這樣不需要很高角動量的方法來說就非常浪費了，稍大點的體系很難算得動。
• Turbomole的x2c-基組：def2的作者對標量和二分量X2C計算專門優化了全電子2-zeta和3-zeta基組，最高支持到Rn。包括x2c-SV(P)all、x2c-SVPall、x2c-TZVPall、x2c-TZVPPall，質量與def2系列同檔次版本相仿佛，并且也定義了/J輔助基組。質量很理想
  
• WTBS：元素涵蓋全面，是極小基，收縮程度很大，是給非相對論HF計算專用的，能接近HF能量極限。由于不是為相對論而優化的收縮的形式，用于相對論時收縮形式必須得改改才行。
• UGBS：元素涵蓋廣，由于完全沒有收縮所以雖然不是特意給相對論計算提出的但也毫無問題，GTF雖不多但是基函數很多，計算費時，在Gaussian中還很容易出現積分精度不夠的錯誤（需要用較高數目徑向積分格點來解決，如int=grid=400434關鍵詞）。此基組可以加上極化和彌散函數（見Gaussian手冊）。
• Sapporo全電子相對論基組：札幌的一批人提出了DZ、TZ、QZ級別的用于非相對論和DKH3相對論計算的片段收縮基組，可在此獲得：http://sapporo.center.ims.ac.jp/sapporo/Order.do。
• Hirao等人在JCP,115,4463(2001)提出了專供DKH3相對論計算的基組，參數通過最小化原子DKH3標量相對論計算得到的能量得到，覆蓋1-103號元素，定義可在此獲得：http://www.riken.jp/qcl/publications/dk3bs/periodic_table.html。
• Dyall搞了專門適合四分量、二分量計算的基組，是完全非收縮的，從第四及之后周期有定義，DZ、TZ、QZ級別都有。可在此獲得：http://dirac.chem.sdu.dk/basisarchives/dyall/。
• RPF-4Z是一種性價比很高的4-zeta級別用于相對論計算的基組。RAGBS（Relativistic Adapted Gaussian Basis Sets）號稱是相對論計算時沒有變分塌陷問題的基組中最小的一種，可在此獲得：http://basis-sets.iqsc.usp.br/basis-sets/。
• Jorge課題組開發了一套DZ到6Z級別的基組，涵蓋了周期表大部分，有用于普通計算和專為DKH2計算優化的兩類基組。可在此獲得：http://qcgv.ufes.br/downloadbasis.html。
•  DGauss系列基組是專做DFT的DGauss程序御用的基組。定義了從H~Xe，對非相對論LSDA計算優化得到，包含不同檔次。只有B~Ne公開發表了，其它的內置于DGauss和Gaussian程序中（在后者中叫DGDZVP、DGDZVP2、DGTZVP，DG=DGauss）。不適合直接用于相對論計算。
  

有興趣者可參考《從H到Lr所有元素的全電子波函數文件》（http://sobereva.com/235），其中用全電子基組結合DKH2計算了H到Lr所有元素，輸入輸出文件都提供了。

2 大核贗勢和小核贗勢

量子化學中通常指的贗勢(ECP)就是用有效勢來表達內核電子的勢場，這樣就只需要計算外層的，也即與化學過程密切相關的電子就可以了，大大節約了計算量。使用了等效體現相對論效應的贗勢(RECP)還避免了處理重元素時需要做標量相對論計算，量化研究中常用的贗勢都是RECP。簡單來說，如果被贗勢代替的內核電子（贗化的電子）較少就叫小核贗勢，用這樣的贗勢計算時，價層、次價層乃至更內層的電子都需要明確計算。如果贗化的電子數較多，就叫大核贗勢，用這類贗勢時只有價層電子是要明確地計算的。

顯然，小核贗勢比大核贗勢耗時，原則上說要比大核贗勢好，但是實際上也并非總如此，比如對于鑭系錒系，記得Xiaoyan Cao就指出用DFT時小核贗勢（及全電子基組時）結果還不如用大核贗勢，可能是DFT都沒對這么重的元素優化。如果能確定價態的話，DFT最好使用相應價態的大核贗勢（有些贗勢對于原子在不同價態時是不同的）。

對于一個不熟悉的贗勢，想知道是大核贗勢還是小核贗勢，一種方法是看原文，另一種方法就是到BSE（https://www.basissetexchange.org），如果能查到相應贗勢，從中查看有多少電子被贗化。比如對于SBKJC贗勢的碘，用Gaussian格式查看贗勢定義，就會看到這一行
I-ECP     3     46
這代表46個內核電子被贗化了，由于碘是53個電子，所以明確考慮的有7個電子，顯然就是5s 5p價層上的7個電子，因此，這是大核贗勢。而對于cc-pVDZ-PP，對于碘被贗化的是28個電子，因此比SBKJC時要多考慮18個電子，自然這就是4s 4p 4d上的2+8+10=18個電子（不考慮4f，因為其能量較5p更高）。對于碘，由于cc-pVDZ-PP不僅考慮了價層的電子還考慮了次價層的電子，所以是小核贗勢（具體來說是小核贗勢基組，見下文）。

Lanl系列贗勢中，對于過渡金屬，Lanl1贗勢就是大核贗勢，比如對于Fe，只考慮了4s 4d。Lanl2贗勢就是也考慮了次價層電子的小核贗勢，比如對于Fe，考慮了3s 3p 4s 4d。特別要注意的是，lanl1、lanl2對于主族元素（不包括第一和第二主族，下同）的定義是完全一樣的，且是大核贗勢！順帶一提，和多數贗勢不同的是，Lanl贗勢對第四周期元素并沒有考慮相對論效應，Kr之后的才考慮了。

SBKJC贗勢對主族是大核，對過渡金屬是小核。

Stuttgart (SDD)贗勢對多數過渡金屬是小核，對一些過渡金屬尤其是稀土金屬元素也有大核版本。在Gaussian中使用時，可以用XYn形式來指定，n代表多少內核電子被贗化，可用的組合方式見手冊Pseudo關鍵詞末尾部分。Stuttgart贗勢對于主族也有大核和小核贗勢之分，不過Gaussian內置的只有大核的。

贗勢和贗勢基組是兩個概念，只有二者結合才能用于實際計算。lanl系列、SBKJC系列、Stuttgart系列等各種贗勢都有官方標配的贗勢基組，也有一些第三方的贗勢基組，比如(aug)-cc-pVnZ-PP實際上不是贗勢的名字，而是基于Stuttgart贗勢優化出的適合于相關計算的贗勢基組。(aug)-cc-pVnZ-PP對于所有元素，無論過渡金屬還是主族，結合的都是小核版本的Stuttgart贗勢。def2系列基組，諸如def2-TZVP，對于前四周期都是全電子基組，但是從Rb開始，盡管還是叫做原先的名字，但本質上就成了贗勢基組了，結合的贗勢都是Stuttgart小核版本的贗勢。

由于絕大多數贗勢對主族元素都是大核贗勢，如果想用小核贗勢算主族元素，推薦用(aug)-cc-pVnZ-PP或def2系列，它對過渡金屬和主族元素都是小核贗勢。用Stuttgart小核贗勢結合標配基組算主族也可以，但是Gaussian里沒自帶，需要去Stuttgart贗勢的主頁上拷貝下來：http://www.theochem.uni-stuttgart.de/pseudopotentials/clickpse.en.html。

3 贗勢下能做波函數分析乎？

原理上說，只要嚴格滿足了兩個前提，贗勢下做波函數分析就完全沒問題：1 關心的只是價層或更外層的性質 2：只有很內層的電子被贗化，價層的電子基本不受影響。

舉幾個例子，假設前述條件完全滿足：
1 分析靜電勢時，由于實際體系中原子內核的密度基本不受化學環境的影響而是球對稱分布的，因此它對分子價層及以外區域的靜電勢的貢獻會與等電量的核電荷產生的貢獻精確抵消，所以把內核電子贗化不會影響分子價層及以外的靜電勢的分析。
2 分析鍵級時，由于內層軌道不參與成鍵，因此顯然可以放心地把內層電子贗化。
3 計算原子間離域性指數時，由于內核電子定域性很強，幾乎不會離域到其它原子上，因此贗化內層電子不影響結果。
4 通過等值面分析電子定域化函數(ELF)時，如果贗化內核電子，盡管會導致看不到內核電子的定域化域，但由于價層波函數仍保持實際情況，原子間共價鍵和孤對電子的定域化域依然照常存在。
5 分析AIM臨界點、鍵徑時，由于贗化了內核電子，使得原子核位置的電子密度的峰消失，會導致找不到核臨界點(NCP)（就算找到也只是數值原因產生的巧合）。而且由于密度不再是簡單地從核中心單調下降而變得略復雜，因此在內核附近區域還可能會產生一大堆額外的臨界點混淆視聽。但是對應成鍵的鍵臨界點(BCP)，以及環和籠的臨界點(RCP和CCP)的位置和性質并不會因此發生變化，完全可以照常分析。而尋找鍵徑時自然不可能再找到連接核中心NCP和對應成鍵的BCP的鍵徑，不過諸如BCP-RCP這樣原本只在價層出現的鍵徑則不會受到影響。

上述這些都是在滿足那兩個前提的理想條件下的結果。但實際情況更復雜，因為內核電子并不只分布在內核區域，它們還與價層區域電子存在重疊，這是因為用作基函數的高斯函數是平滑衰減的，內核區域的高斯函數因此會在價層區域露著“尾巴”。因此，對于基于實空間的波函數分析，比如分析密度、ELF等，贗化內層電子會多少影響到價層性質的定量乃至定性結果。即便是基于希爾伯特空間的分析，比如計算Mayer鍵級，一些價層以內的電子，尤其是次價層電子，也可能一定程度地表現出價層電子的性質，涉及到原子間相互作用，因此將它們贗化會影響分析結果。

由于次價層電子最接近價層電子，所以為了盡可能少地避免贗勢對波函數分析結果的影響，次價層電子總應當保留，也就是總應當使用小核贗勢，否則可能定性正確的結果也得不到。但是也不代表用小核贗勢就保險了，贗化比次價層更內層的電子依然可能對價層的分析帶來一定定量的擾動，但起碼此時定性正確的結果可以獲得。在JCC,18,416中，作者分析了贗勢對AIM分析的影響，發現使用小核贗勢時雖然價層的動能密度、拉普拉斯值與實際有些偏差，但結果還是可信的，該有的臨界點也都有，而使用大核贗勢時，有的BCP就不見了。

4 wfx文件與EDF信息

wfn文件格式是最早Bader小組開發的AIMPAC程序引入的，用于記錄分子坐標、高斯函數和軌道波函數信息，已被Gaussian, GAMESS-US/UK, Firefly, Q-Chem等眾多從頭算軟件所支持。這種格式的三大弊端是：
1 高斯函數最高角動量只支持到f。雖然其格式實際上也可以用于記錄更高角動量，但是更高角動量的標號沒有統一定義（不過從G09 B.01開始輸出的wfn文件已經支持更高角動量了，Multiwfn也可以正確載入）
2 沒有明確記錄軌道自旋類型，是alpha還是beta軌道有時得靠猜
3 數據的記錄精度有限

Gaussian從Gaussian09 B.01開始支持了wfx格式，用out=wfx即可輸出，它是對wfn文件的擴展。wfx格式搞得很繁復，我個人認為格式定義得不好，程序讀取起來也麻煩，但是帶來了幾大好處：
1 高斯函數角動量可以支持到無限高，標號有統一的定義
2 軌道自旋類型有了明確的記錄
3 數據的記錄精度得到了提高

除了在上述三點的改進，wfx格式最關鍵的是加入了一個EDF(electron density function)字段，這對于在使用贗勢時進行基于密度的分析十分有用。這個字段記錄了贗勢計算時內核電子的分布，將這部分密度和實際計算得到的價層的密度疊加起來和全電子計算時的密度是幾乎一致的。各元素內核部分的密度是內置在Gaussian程序中的，是事先通過多個s型高斯函數對全電子相對論獲得的內核密度擬合得到的，會在Gaussian計算完畢后輸出wfx文件時直接寫入EDF字段。關于擬合有關的更詳細信息見JPCA,115,12879。EDF字段定義的內核密度和贗勢種類無關，是看元素和被贗化的電子數。比如對于Cu使用lanl2DZ和SDD，由于都是小核，所以.wfx文件里的EDF字段的信息是相同的。

注意：對于極個別體系，Gaussian某些版本生成的wfx文件的EDF字段信息不對。比如Ti原子在lanl2下本該有10個內核電子，但是在個別含Ti體系中，G09給出的wfx文件里對應的EDF字段只描述了4個電子，導致此錯誤的原因目前不明（在C.01版仍存在）。如果不放心，可以在使用Multiwfn分析wfx文件時，先進入主功能100里的功能4，然后選1，這會對整個空間的電子密度進行積分。如果積分的結果和總的電子數（內核電子+價電子）比較接近，那么說明內核電子的密度都在EDF字段中正確描述了。

Multiwfn支持wfx格式，并且會自動利用其中EDF字段的信息把內層電子密度在計算總密度（包括各種基于密度的函數，但靜電勢除外）的時考慮進去。不過從頭算程序對wfx支持得還很不普遍。Firefly程序也支持產生.wfx文件，但是并不會把EDF字段寫進去。

欲更詳細了解wfn格式見《高斯fch文件與wfn波函數文件的介紹及轉換方法》（http://sobereva.com/55）。欲了解wfx格式的具體定義參見http://aim.tkgristmill.com/wfxformat.html（可能需要fan wall）。

Multiwfn從3.4版開始有了一個極為重要改進，就是內置了Wenli Zou等人構建的EDF內核電子密度數據庫。如果載入的.wfx文件里有EDF字段，那么Multiwfn照常從.wfx文件里讀取EDF字段；而當.wfx文件里沒有EDF字段，或者用的是其它包含波函數信息的輸入文件（如.wfn、.fch、.molden等），默認情況下Multiwfn對于使用了贗勢基組的原子會自動從自帶的EDF數據庫中選取合適的EDF信息載入，從而提供對內核電子的描述。而且這個Multiwfn內置的EDF數據庫比Gaussian自帶的EDF數據庫的質量更好，支持的元素更多。因此，對于Multiwfn來說，使用wfx文件在描述內層密度上并沒有任何優勢。

如果不想讓Multiwfn從.wfx文件中讀取EDF字段，或者不想從自帶的EDF庫中自動讀取EDF信息，可以用settings.ini文件里的readEDF和isupplyEDF參數修改規則，詳見此文件中的注釋。另外，Multiwfn中也允許從原子的.wfx文件中載入EDF信息。詳見Multiwfn手冊附錄4，不過這種做法意義不大，因為靠內置的EDF庫就足夠了。

注意EDF信息提供的只是內層電子的密度，在Multiwfn中也只有直接依賴于電子密度的實空間函數，比如電子密度、電子密度梯度模、電子密度拉普拉斯值、約化密度梯度(RDG)等的分析會受到EDF的影響，而其它的實空間函數，比如ELF、自旋密度、靜電勢，以及不涉及實空間函數的波函數分析，如Mayer鍵級、Mulliken分析，都不會受到EDF的任何影響。

雖然EDF可以讓贗勢下的與電子密度有關的分析和全電子基組時符合較好，但是，筆者依然鼓勵用小核贗勢，因為在考慮EDF下，小核贗勢的結果還是會比大核贗勢更接近全電子基組的結果。

5 總結

綜上所述，如果在贗勢下做波函數分析，一定要用小核贗勢。如果用了小核贗勢，但還是還想獲得更高精度結果，或者懷疑結果可能有誤，那么不如干脆直接上全電子標量相對論計算（對第四周期元素，忽略相對論效應用普通基組計算也可以）。如果是在贗勢下分析和密度相關的問題，如根據AIM理論對電子密度做拓撲分析、分析電子密度拉普拉斯值，用Multiwfn可以非常放心，因為無論哪種格式的波函數文件，程序都會自動讀入EDF信息，從而把內核密度如實地描述出來，結果會和全電子基組時相符很好。