贗勢的函數形式以及在量子化學程序中定義的方式
1 前言
贗勢在量子化學計算中使用得極為普遍，有時需要用的贗勢在量子化學程序中沒有內置，就得通過自定義的方式添加，例如Stuttgart贗勢對于Au按照Gaussian程序格式的寫法是
AU 0
AU-ECP 4 60
g-ul potential
1
2 1.000000000 0.000000000
s-ul potential
2
2 13.205100000 426.709840000
2 6.602550000 35.938824000
p-ul potential
2
2 10.452020000 261.161023000
2 5.226010000 26.626284000
[...略]
雖然定義的格式對于各個量化程序都大同小異，大家也經常接觸到，但是其含義很多人不很明白。一方面是主流的量化書籍中往往對贗勢都不怎么詳談，另一方面網上也沒有什么資料專門討論這個問題。本文主要目的就是以不很嚴謹但是淺顯易懂的語言說明白通常用的贗勢的函數形式，特別是在量化程序中是怎么靠類似上述內容定義出來的。在此之前先談一下贗勢的基本特點。
2 贗勢的基本特點
本文說的贗勢是量化研究中最常用的贗勢，如Hay和Wadt提出的Los Alamos贗勢（或稱Lanl贗勢）、Stuttgart/Dresden贗勢（或稱SDD贗勢）、SBKJC贗勢（也稱CEP贗勢）等等。
贗勢的主要用處在于：(1)將化學上不感興趣的內層電子以等效的勢場描述，從而不需要將內層電子顯式地表達出來，大大節約了計算量 (2)可以等效地體現相對論效應。從第四周期相對論效應開始顯現，但不考慮也問題不大；而對第五周期及之后相對論效應是不可忽視的，不考慮甚至結果定性錯誤。
使用贗勢后，所得的軌道可以稱為贗軌道。最低的一批贗軌道的能量和全電子計算時各個價層軌道能量接近或嚴格相一致；軌道形狀上，在內層區域，這些贗軌道沒有了實際價層軌道的節點，而在價層及靠外的區域（即大于cutoff半徑rc的區域）贗軌道和實際軌道的形狀是很接近或嚴格一致的。只要贗勢設定得合適，靠贗軌道一般可以合理表現出實際價層軌道所展現的性質，如成鍵方式、價層區域電子密度分布等。
贗勢一般是對于元素來定義的，有的時候還根據特定的價態、考慮相對論的方式、被贗化的電子數等因素來具體分為不同版本。對于分子體系，贗勢就是體系中各個原子對應的贗勢的加和。所以這也要求贗勢要有較好可移植性，適用于不同化學環境。
贗勢的定義是依賴于被作用的波函數的角動量的（即依賴于角量子數l）。而對于主量子數n和磁量子數m的依賴在一般的研究中是不考慮的。（2、4分量相對論計算時所用的贗勢還依賴于總角量子數j，本文不涉及）
構建贗勢的方法各有不同。例如Los Alamos贗勢是形狀一致型贗勢的典型，其構建的方法是：(1)對指定的原子，選定好哪些軌道作為內核軌道，哪些作為價軌道，并且選定原子所處的組態 (2)對原子做數值Hartree-Fock全電子計算得到價層軌道波函數 (3)把價層軌道波函數的r<rc部分的節點抹掉而成為平滑函數（以多項式表達），r>rc區域的軌道波函數則不作改動。這樣就構建出了贗軌道 (4)對于不同角動量的贗軌道，根據贗軌道波函數及其能量反推出一個勢函數（暫以列表方式數值表達），使得贗軌道正好是HF方程外加這個勢函數下的本征函數，且相應的本征值恰為實際價軌道的能量 (5)將數值描述的勢函數以高斯函數擬合來表達。由于勢函數是平滑的，所以這容易做到 (6)選取一定數量的原始高斯函數(GTF)，然后優化它們的指數和收縮系數，使得它們能盡可能好地表現贗軌道，這就定義出了與此贗勢對應的贗勢基組。
而以Stuttgart贗勢為代表的能量一致型贗勢的構建思路與形狀一致型贗勢有明顯差別。它將目標數據選取為原子及其離子態的成百上千個電子態的總價電子能量（而非一個個價軌道能），通過不斷調整贗勢參數，使贗勢與相對論全電子計算得到的這些態的總價電子能量的平均差異盡可能小。盡管此類贗勢沒有嚴格要求r>rc范圍內贗軌道和實際軌道完全一致，但結果上看卻也非常接近。
顯然，用贗勢就得用價電子基組。贗勢在提出時都有標配的贗勢基組。同一種贗勢也可以有多種不同質量或不同研究者定義的贗勢基組，如Los Alamos贗勢有Lanl2DZ、Lanl2DZdp、Lanl08、LANL2TZ、LAV或LACV系列基組等等。對于Stuttgart贗勢，除了它標配的基組外，(aug)-cc-pVnZ-PP系列、較重元素的def2-系列基組用的也都是Stuttgart贗勢。如果沒有經驗，對于特定贗勢，一般不建議使用未專門對它優化過的價層基組，不同贗勢各自適用的贗勢基組也不建議隨意互換，這是因為不同贗勢對應的贗軌道在內層區域徑向形狀是有差異的，亂用贗勢基組會造成不能合理表現這部分區域。
3 贗勢的函數形式與角動量問題
使用贗勢時，體系的哈密頓為H = ∑[i]h'_i + ∑∑[i>j]1/r_ij，第一項是價電子的有效單電子哈密頓的加和，第二項是價電子間的靜電相互作用。有效單電子哈密頓為h'=h+∑[a]U_a，其中h是常規的單電子哈密頓，即動能+核吸引勢。U_a代表a原子對應的贗勢。
有了H，就可以基于贗勢求解化學問題了。贗勢對于HF/DFT/post-HF計算都可以用。對于HF（或DFT）而言，HF方程只考慮價層電子，fock算符中包含只由價層電子構建的庫侖項和交換項、∑[a]U_a以及核吸引勢。這樣求解HF方程得到的就是贗價層軌道。
下面來看看贗勢的函數形式。由于贗勢都是球形平均的，所以只需考慮徑向坐標。各種贗勢都寫為這樣的通式：
U(r)=U_L(r)+∑[U_l(r)-U_L(r)]*P_l （加和是對于角量子數l，從0到L-1）
式中L往往取被贗化的內核軌道中的最高角動量+1。P_l=|l><l|代表l角動量投影算符，可以進一步寫為對應的磁角動量投影算符的加和P_l=∑[m=-l→l]P_l,m。比如P_2投影算符投影出的是d成份（l=2），若將它作用在p軌道（l=1）上，結果就為0；若作用在d軌道上，則軌道保持原樣；若作用在一個復雜的軌道上，則得到徑向、角度方向可分離變量形式R(r)*∑c_lm*Y_lm(θ,φ)，其中R是此軌道的徑向部分函數，Y_lm是l,m型球諧函數，c_lm則是此軌道中Y_lm球諧函數的分量。
上述贗勢U(r)表達式表現的含義是對于任何軌道都施加U_L勢。而對于軌道中的角量子數l在0到L-1之間的成份，則進一步施加相應的校正勢U_l-U_L。換句話說，軌道中角動量l>=L的部分感受到的是相同的U_L勢，而角動量l<L的部分感受到的則是依賴于l而因此各自不同U_l勢。
前面已提到贗勢是依賴于角動量的。假設內核軌道中角動量最高的為t，那么，對于價層軌道l=0,1,2...t角動量成份都必須單獨定義相應的贗勢。這是由于不同角動量的內層軌道與價軌道的相應角動量成份之間有著明顯不同的相互作用（具體來說，是正交以及交換勢）。對于l>t的情況，隨著l增大，l與l+1的勢的差異很快就變得十分微小，這是為什么在通常的贗勢中，會選取一個L，令l>=L的角動量成份都使用相同的勢U_L，而不再對它們單獨定義勢。L的選取使情況而定。例如，對于第一周期過渡金屬，Hay和Wadt發現d和f軌道對應的勢差異比較顯著，而f和g的勢就沒什么差別，所以對于第一周期過渡金屬Los Alamos贗勢中對于f及以上角動量所用的贗勢一致，即L=3，而對s、p、d都單獨定義了贗勢。而對于第三周期過渡金屬則需要明確考慮更高角動量，Los Alamos贗勢將L設到了4 (g)，Stuttgart贗勢甚至設到了5 (h)。
注意贗勢基組的最高角動量和贗勢定義中的L沒必要一致，一般L是大于或等于贗勢基組中最高角動量+1的。比如對于Cu，Los Alamos贗勢的L=3，而Lanl2DZ贗勢基組中Cu的基函數的最高角動量是d，用以描述價層d電子。之所以基函數沒有f但是贗勢中卻定義了角動量大于等于f成份所感受到的勢，是因為實際計算的是分子體系而非原子體系，贗勢的作用對象不是當前原子的基函數而是贗軌道，其它原子的基函數都會參與贗軌道的構成，因此贗軌道的形狀是復雜的，投影算符對它們是可以投影出大于等于f的成份的，這部分也應感受到贗勢，所以贗勢不能只給s,p,d成份來定義。但是，假設計算的只是Cu這一個原子，那么此時的贗軌道，換句話說也就是沒有內部節點的原子軌道，最高角動量只能到d，而投影不出f及以上的成份，因此此時只需要定義s、p、d受到的贗勢就夠了，f及以上成份感受的贗勢即使定義了也沒用。
為了實際求解時方便，數值描述的贗勢一般用多項式與高斯函數的乘積來擬合來得到解析表達式。贗勢一般的展開形式都是統一的：
r^2*U(r)=∑[k]d_k*r^n_k*exp(-ζ_k*r^2)，等價于U(r)=∑[k]d_k*r^(n_k -2)*exp(-ζ_k*r^2)。
其中r是徑向坐標，d_k是展開系數，n_k可以為0、1或2，ζ是高斯函數的指數。正是由于贗勢的實際表達形式如同高斯型基組一樣遵循統一的規范，所以各種贗勢原則上都可以在各種支持贗勢的量化程序中使用。通常贗勢表達方法是先對U_L(r)用上式進行擬合，得出對于l>=L角動量成份的贗勢的d、n、ζ數據，然后對于l<L的角動量部分，分別按上式擬合出對U_L贗勢的修正量U_l(r)-U_L(r)的d、n、ζ數據。
當量子化學程序使用某贗勢時，必須有能力做到L角動量基函數的積分。例如GAMESS-US的積分代碼支持的基函數最高只到g，而對于Hg，Stuttgart贗勢的L=5 (h)，所以盡管其贗勢基組的最高角動量只是d，在GAMESS-US中也沒法用。注意輸出的波函數中最高角動量函數與贗勢的L無關。比如標準的wfn文件，最高只能記錄到f角動量高斯函數，對于含Hg體系并且用了Stuttgart贗勢的情況，由于贗勢基組自身沒有高于f角動量的基函數，若其它原子所用的基組也沒有高于f角動量的基函數，那么wfn文件就可以正常記錄此體系的波函數信息。對于波函數分析程序，諸如Multiwfn，支持的高斯函數最高角動量是h，用于這樣的體系的分析也是完全沒有問題的。
4 贗勢的定義格式
贗勢在不同文獻中的表達方式、在不同量子化學程序中定義的格式都是大同小異的。下面按照Gaussian程序中定義贗勢的格式來進行說明。贗勢基組的定義方式和普通的基組定義方式沒有任何區別，所以這里不提及。
下面這是Zn的Los Alamos贗勢的定義，//后面是注釋：
ZN 0 //ZN是被定義贗勢的元素名，如果只對某個原子定義，則可以寫原子序號。0是必須寫的
ZN-ECP 3 18 //贗勢名字；贗勢的L值，比如此例對于f及以上角動量用相同的勢，故L=3；用贗勢描述的內核電子數，這里為18。贗軌道中不同角動量成份所受的勢分別寫在接下來的一個個block里。第一個block定義的是l>=L成份所受的勢（即U_L），接下來依次定義l=0、l=1、l=2...直到l=L-1對應的勢相對于U_L的修正量
f-ul potential //從這里開始是U_L勢的block。這一行是block的名字，可以隨意設，不影響結果
5 //U_L展開為前述的∑[k]d_k*r^(n_k -2)*exp(-ζ_k*r^2)形式，k=1、2、3、4、5
1 386.7379660 -18.0000000 //n_1、ζ_1和d_1
2 72.8587359 -124.3527403 //n_2、ζ_2和d_2
2 15.9066170 -30.6601822 //n_3、ζ_3和d_3，以此類推
2 4.3502340 -10.6358989
2 1.2842199 -0.7683623
s-ul potential //這個block定義的是s角動量成份所受的勢相對于U_L的修正，即U_0-U_L。這行是block的名字，也可以隨意設
5 //U_0-U_L也是作如上形式的展開，展開為以下5項。下面每一行都是相應項的n、d和ζ。
0 19.0867858 3.0000000
1 5.0231080 22.5234225
2 1.2701744 48.4465942
2 1.0671287 -44.5560119
2 0.9264190 12.9983958
p-ul potential //和上一個block類似，此block定義的是p角動量部分對于U_L的修正，即U_1-U_L
5
0 43.4927750 5.0000000
1 20.8692669 20.7435589
2 21.7118378 90.3027158
2 6.3616915 74.6610316
2 1.2291195 9.8894424
d-ul potential //和上一個block類似，此block定義的是d角動量部分對于U_L的修正，即U_2-U_L
3
2 13.5851800 -4.8490359
2 9.8373050 3.6913379
2 0.8373113 -0.5037319
由此例可見，每一個block的名字都是隨便起的，比如d-ul也可以寫為d and up，s-ul可以寫成s-d，或者瞎寫一個名字，諸如寫寫Mio_akiyama也無妨。關鍵的是，block的順序不能錯，一定是先定義U_L，再按角動量由低到高依次定義U_l-U_L，l=0,1,...L-1。
5 修改贗勢的寫法
前面已經說明了贗勢的形式和定義方法，為了加深理解，這一節我們對贗勢的寫法進行各種修改來其考察影響。
下面是Pt的Stuttgart贗勢
PT 0
PT-ECP 5 60
h-ul potential
1
2 1.000000000 0.000000000
s-ul potential
2
2 13.428651000 579.223861000
2 6.714326000 29.669491000
p-ul potential
2
2 10.365944000 280.860774000
2 5.182972000 26.745382000
d-ul potential
2
2 7.600479000 120.396444000
2 3.800240000 15.810921000
f-ul potential
1
2 3.309569000 24.314376000
g-ul potential
1
2 5.277289000 -24.218675000
值得一提的是這個勢對于大于等于h角動量的成份其實是沒有勢的作用的，因為其block只有一項，而且系數d為0，這一項其實只是走個形式。之所以對>=h角動量沒有定義贗勢，這和前面提到的Stuttgart贗勢生成方式有關，因為定義>=h的贗勢對計算原子各種組態的總價電子能量的影響是可忽略的，于是就不去弄相應的勢了。
如果我們計算的是單個Pt原子，由于價軌道只由自身的s、p、d基函數來描述，沒有f及以上角動量成份，所以贗勢定義簡化為下面這樣，即去掉f、g以及>=h的定義，對Pt單原子計算結果是沒有任何影響的：
PT 0
PT-ECP 3 60 //L=3，以下依次定義>=f的勢以及s、p、d對勢的修正
f-ul potential
0 //含0項，>=f的成份將不受贗勢作用。盡管計算Pt單原子考慮f是多余的，但按照基組定義格式，必須先定義>=L所受的勢，所以這里走個形式
s-ul potential
2
2 13.428651000 579.223861000
2 6.714326000 29.669491000
p-ul potential
2
2 10.365944000 280.860774000
2 5.182972000 26.745382000
d-ul potential
2
2 7.600479000 120.396444000
2 3.800240000 15.810921000
如果我們計算的是含Pt化合物，做將原始的贗勢改為如上這樣就會一定程度影響結果，因為贗軌道比原子軌道復雜得多，因此可以投影出>=f的成份。對于那些最高角動量有比較大限制的程序，比如GAMESS-US只能做到g，若非要用這樣L達到h（或更高）的贗勢，就不得不人為地修改贗勢來對角動量定義進行截斷（也就是讓贗勢定義只包含>=g,s,p,d,f的定義，且讓>=g部分不受贗勢作用），盡管計算分子時這會一定程度上影響結果。并且非要這么做的話需要在文章中注明用的是自行對角動量進行截斷過的贗勢。另一種辦法就是換其它的贗勢，比如對于Pt，Los Alamos贗勢和SBKJC贗勢的L=4，比Stuttgart贗勢低一階，所以GAMESS-US中就能夠直接用。
再舉一個例子，對于前面的Zn的Los Alamos贗勢，我們進行改寫，不設U_L勢（即令之定義為空），且將原先U_L勢和各個角動量的修正勢U_l-U_L直接合并，得到下面的形式：
ZN 0
ZN-ECP 3 18
f-ul potential
0
s-ul potential
10
1 386.7379660 -18.0000000
2 72.8587359 -124.3527403
2 15.9066170 -30.6601822
2 4.3502340 -10.6358989
2 1.2842199 -0.7683623
0 19.0867858 3.0000000
1 5.0231080 22.5234225
2 1.2701744 48.4465942
2 1.0671287 -44.5560119
2 0.9264190 12.9983958
p-ul potential
10
1 386.7379660 -18.0000000
2 72.8587359 -124.3527403
2 15.9066170 -30.6601822
2 4.3502340 -10.6358989
2 1.2842199 -0.7683623
0 43.4927750 5.0000000
1 20.8692669 20.7435589
2 21.7118378 90.3027158
2 6.3616915 74.6610316
2 1.2291195 9.8894424
d-ul potential
8
1 386.7379660 -18.0000000
2 72.8587359 -124.3527403
2 15.9066170 -30.6601822
2 4.3502340 -10.6358989
2 1.2842199 -0.7683623
2 13.5851800 -4.8490359
2 9.8373050 3.6913379
2 0.8373113 -0.5037319
對于Zn單個原子的計算，這種贗勢寫法和原先的贗勢的寫法計算結果是完全一致的，因為s、p、d成份感受到的贗勢與原先相同。但是，對于分子體系計算，由于贗軌道存在f及以上角動量成份，我們又沒設這部分的贗勢，而原先定義中這部分贗勢是非0的，所以計算結果將會有差異。