一硼化鉿的晶體結構預測與構效關系基因

魏曉婷 曾慶豐＊, 張 琪 馮欽顥

（1西北工業大學材料學院材料基因組國際合作研究中心，西安 710072）

（2西北工業大學材料學院超高溫結構復合材料重點實驗室，西安 710072）

最近10年，材料晶體結構預測技術取得了重大突破[1-3]。晶體結構預測技術可能使新材料發現的周期會大大縮短，也成為材料基因組的重要技術之一[4-6]。例如，通過晶體結構預測，研究人員已經發現了一系列新型的超硬材料[7-9]、超導材料[10-11]、拓撲絕緣體[12]等。除了輔助用于新材料的發現外，晶體結構預測技術的另一應用是研究人員能夠通過無偏見的理論預測來獲得任意給定材料體系的潛在晶體結構[1]。這一重要應用有助于消除研究人員對某些材料晶體結構的不確定性，并進而對材料的合成和使用提供指導。近幾年來，研究人員運用晶體結構理論預測技術針對過渡金屬硼化物體系進行了很多研究[13-16]。一方面，是因為過渡金屬硼化物往往兼具十分優異的力學、熱學和電學等性能，是一類應用廣泛和具有重大研究價值的材料[17-18]；另一方面，則是因為過渡金屬硼化物往往能夠形成多種配比的穩定晶體結構，而僅依靠于實驗研究并不能十分完全地發現和確定它們的晶體結構。晶體結構決定材料性質。因此，通過理論預測出不同配比過渡金屬硼化物的晶體結構并進而研究它們的理論性質，無疑將對過渡金屬硼化物的實驗合成和實際應用起到十分重要的指導作用。

本研究采用晶體結構預測技術確定過渡金屬硼化鉿體系（HfB）的晶體結構。根據實驗和理論研究報道[19]，硼化鉿具有2種不同的化學計量比形式，一硼化鉿（HfB）和二硼化鉿（HfB2）。 HfB2是一種常用的超高溫陶瓷材料，具有十分優異的材料綜合性能如高熔點、高硬度、高力學強度等[20]。HfB2具有AlB2型晶體結構，屬六方晶系結構，空間群為P6/mmm。相比于HfB2，HfB的確切晶體結構和性質仍存在諸多疑問。研究人員僅僅猜測HfB可能具有FeB型[19]、CrB 型[21]、MoB 型[22]或 NaCl型[23]晶體結構。基于這些HfB晶體結構進行的理論計算表明，FeB型HfB晶體結構的基態能量最低，是基態下最穩定的晶體結構[21-23]。然而，考慮到迄今為止并無直接證據表明HfB具有FeB型晶體結構，因而有必要利用晶體結構預測技術針對HfB的晶體結構進行全局搜索。特別要指出的是，研究人員通過晶體結構預測已成功糾正過很多過渡金屬硼化物的晶體結構。在本研究中，通過晶體結構預測，發現了2個比FeB型結構具有更低基態能量的HfB晶體結構。通過第一性原理計算，系統地論證了這2個新發現結構的晶格動力學和力學穩定性。進一步地，基于新發現的2個和已報道的4個HfB晶體結構，運用密度泛函微擾理論研究了溫度對HfB晶體結構穩定性的影響。作為補充，還通過第一性原理計算給出了這些HfB材料的力學性質和電子性質。研究結果可為HfB材料的后續實驗和理論研究提供依據。

1 計算方法

采用基于進化算法開發的晶體結構預測軟件USPEX[1-2]對HfB的晶體結構進行全局搜索。僅需指定化學元素種類、原子個數和目標性質等基本參數信息，并結合外部接口軟件對目標性質進行計算，USPEX軟件便可以利用遺傳和變異等進化操作手段發現該指定化學體系下具有預期優選性質的晶體結構。當用戶選擇能量最低作為優化條件時，USPEX軟件搜索得到的便是最穩定的晶體結構。在當前的晶體結構預測中，考慮每一晶胞可包含1～4倍化學計量比的HfB。在密度泛函理論（DFT）[24]的框架下，針對晶體結構預測過程中由USPEX軟件產生的每一晶體結構，運用VASP軟件包[25]對其進行充分的幾何結構優化和總能計算，具體參數設置如下：采用全電子投影綴加波（PAW）[26]方法，廣義梯度近似（GGA）Perdew-Burke-Emzerhof（PBE）[27]形式，平面波截斷能設置為500 eV，倒易空間k點網格分辨率為2π×0.6 nm-1。以上這些參數設置既可以保證晶體結構的能量得到較好收斂又能確保USPEX軟件進行結構預測的準確性和效率。對于預測發現的2個新型低能量HfB晶體結構，連同文獻已報道的4個晶體結構一起，利用VASP軟件并采用更小的倒易空間 k點網格分辨率（2π×0.3 nm-1）對它們的結構、力學和電子結構等性質進行計算。需要說明，考慮到某些HfB晶體結構能量間的微小差異，對以上相關計算參數的選取已進行了嚴格和全面的測試，測試結果表明選用的計算參數是合理的。為了驗證新發現的2個HfB晶體結構的晶格動力學穩定性，采用 Phonopy 軟件[28]基于密度泛函微擾理論（DFPT）[29]計算了這些晶體結構的聲子色散曲線。此外，采用VESTA軟件[30]對HfB晶體結構進行可視化展示。

2 研究結果

2.1 HfB晶體結構的預測

為了探究HfB的確切晶體結構，研究人員已經提出過4種可能的晶體結構類型[19,21-23]：FeB型（空間群：Pnma）、CrB 型（空間群：Cmcm）、MoB 型（空間群：I41/amd）、以及 NaCl型（Fm3m）。 基于這 4 個 HfB 晶體結構進行的理論研究表明：Pnma結構具有最低的基態能量，因而是基態下最穩定的HfB晶體結構。通過對HfB的基態晶體結構進行全面預測，發現了2個比Pnma結構具有更低基態能量的晶體結構：P6m2和R3m。在不考慮零點能的情況下，相比于Pnma結構，P6m2結構的平均每原子能量相比Pnma結構要低12 meV；R3m結構的能量略低于Pnma結構，二者能量相差不到2 meV·atom-1。

圖1 HfB晶體結構的模擬X射線衍射圖Fig.1 Simulated X-ray diffraction patterns for HfB compounds

對新發現的2個以及已報道的4個HfB晶體結構進行了充分的幾何結構優化，這些晶體結構優化后的幾何結構信息列于表1中。對于已報道的HfB晶體結構，表1也列出了它們幾何結構信息的文獻報道值[23]。通過對比發現，本次計算結果與前人的報道結果基本一致，表明了當前計算結果的正確性和可靠性。表1也列出了以上6個HfB晶體結構在基態下的每原子平均能量。在這6個HfB晶體結構中，P6m2結構的能量最低，R3m結構的能量次低，然后依次是Pnma結構，Cmcm結構，I41/amd結構和Fm3m結構。如果僅從能量角度考慮不同結構間的穩定性差異，以上6個HfB晶體結構在基態下的穩定性由高到低的次序為：P6m2→R3m→Pnma→Cmcm→I41/amd→Fm3m。需要指出的是，前4個HfB晶體結構（P6m2，R3m，Pnma 和 Cmcm 結構）間的能量差異較小 （≤20 meV·atom-1），意味著這4個HfB晶體結構實驗存在的可能性較大并需要被重點關注。為了便于后續的實驗驗證，給出了這6個HfB晶體結構的模擬X射線衍射圖，如圖1所示。

2.2 HfB晶體結構的穩定性

理論上，判斷一個晶體結構在基態下的穩定性，除了需要考慮晶體結構的基態能量外，還需要考慮晶體結構的晶格動力學和力學穩定性。晶體結構的晶格動力學穩定性可以由其聲子色散曲線反映出來，力學穩定性則可以通過驗證其彈性常數是否符合Born-Huang穩定準則[31]進行判斷。

表1 基態條件下HfB晶體結構的空間群、晶格參數、原子Wyckoff位置和能量Table 1 Space group,lattice parameter,Wyckoff position,and energy of HfB compounds at 0 K

對于新發現的P6m2和R3m結構，首先分別構建了3×3×1和 2×3×3超胞結構，然后利用這2個超胞結構計算了它們的聲子色散曲線。如圖 2（a）、（b）所示，在整個第一布里淵區內，P6m2和R3m結構的聲子振動頻率均無虛頻，說明這2個新發現的HfB晶體結構在基態下均具有晶格動力學穩定性。本研究也補充驗證了已報道的4個HfB晶體結構在基態下的晶格動力學穩定性，如圖 2（c）～（f）所示。

針對以上6個HfB晶體結構，進一步計算了它們的彈性常數并依據Born-Huang穩定準則對它們的力學穩定性進行了判斷。通過判斷，這6個HfB晶體結構在基態下均具有力學穩定性。通過計算獲得的HfB晶體結構的獨立彈性常數列于表2中。

圖2 HfB晶體結構在基態下的聲子色散曲線Fig.2 Calculated phonon dispersion curves for HfB compounds

表2 基態條件下HfB晶體結構的獨立彈性常數CijTable 2 Computed elastic constants Cijfor HfB compounds at 0 K

2.3 溫度對HfB晶體結構能量的影響

上文中已經比較了不同HfB晶體結構在基態下的能量。這里進一步考慮溫度對HfB晶體結構能量的影響。當不考慮非諧效應時，晶體結構在零壓強以及不同溫度條件下的自由能可通過如下公式進行計算式中F（T）即晶體結構在不同溫度T下的自由能，Etot是晶體結構在基態下的能量（表1），Ezpe是晶體結構的零點能（Ezpe= ∫f（ω）hωdω），k 為玻爾茲曼常數，h 為普朗克常數，f（ω）是聲子態密度,ω是聲子振動頻率。基于密度泛函微擾理論，通過計算得到了6個HfB晶體結構的聲子態密度，并由此獲得了這些HfB晶體結構在不同溫度下的自由能，如表3所列。需要指出，通過對聲子態密度進行積分可以獲得晶體結構的零點能 （常用于修正晶體結構在基態下的能量）。對比表1和表3可以發現，考慮零點能修正后，基態下6個HfB晶體結構間的能量大小次序會有所變化。具體是，通過預測發現的R3m結構在0 K時的能量會高于已報道的Pnma結構，而其他結構的能量大小順序沒有改變。零點能修正的結果表明：在不考慮零點能的情況下進行晶體結構預測時，研究人員不能僅僅只關注能量最低的結構，還需要關注其他一些能量較低的結構。

表3 HfB晶體結構在不同溫度下的自由能Table 3 Computed temperature-dependent free energy for HfB compounds

利用表3的數據，以Pnma結構的自由能作為零基準，進一步獲得了其它5個HfB晶體結構的自由能隨溫度的變化曲線，如圖3所示。由于Fm3m結構與 Pnma結構的能量相差過大（＞200 meV·atom-1），因此圖3中并沒有顯示出Fm3m結構的信息。不難發現，溫度可以顯著改變不同HfB晶體結構間的能量差異甚至是相對大小。例如，在0 K下，新發現的P6m2結構比Pnma結構具有更低的能量。但是，當T＞777 K后，Pnma結構的能量卻比P6m2結構的能量低。類似的現象同樣可見于R3m結構和Cmcm結構：當T＜585 K時，R3m結構的能量比Cmcm結構的能量小；而當T＞585 K時，Cmcm結構比R3m結構具有更低的能量。

圖3 HfB晶體結構在0～1 000 K溫度范圍內的自由能（以Pnma結構作為零基準）Fig.3 Computed free energy for HfB compounds at 0～1 000 K

由圖3可以發現在0～1 000 K范圍內，雖然P6m2、R3m、Pnma和Cmcm這4個HfB晶體結構間的能量大小次序發生了變化，但是這些結構間的能量差異卻始終維持在20 meV·atom-1以內。這說明，即便在高溫條件下，這4個HfB晶體結構也均有較大可能被實驗合成和發現。而對于I41/amd和Fm3m結構來說，由于它們具有相對較高能量，推測它們實際存在的可能性要低得多。

2.4 HfB晶體結構的結構特征

圖4為6個HfB晶體結構的示意圖。如圖4（a）、（b）所示，在P6m2和 R3m結構中，每個B原子與6個Hf原子相連形成Hf6B正三棱柱，每2個Hf6B三棱柱通過共面形成Hf8B2四棱柱，B原子相互連接形成二維類石墨烯B層，B-B鍵長約為0.182 nm。如圖 4（c）～（e）所示，在 Pnma、Cmcm 和 I41/amd 結構中，每個B原子與6個Hf原子相連形成Hf6B三棱柱（近正三棱柱，略有扭曲），每2個Hf6B三棱柱通過共面形成Hf8B2四棱柱 （近棱形四棱柱，略有扭曲），B原子相互連接形成zig-zag鏈而不組成二維類石墨烯B層，B-B鍵長范圍為0.189～0.199 nm。在Fm3m結構中，每個B原子與6個Hf原子相連形成Hf6B正八面體，B原子之間不相鄰，如圖4（f）所示。

圖4 HfB晶體結構的示意圖Fig.4 Crystal structures of HfB

依據B原子的分布形式，將這6個HfB晶體結構大致分為3大類；第Ⅰ類：含二維類石墨烯B層結構（P6m2和R3m結構）；第Ⅱ類：含zig-zag形B鏈結構（Pnma、Cmcm 和 I41/amd 結構）；第Ⅲ類：孤立 B原子結構（Fm3m結構）。

因為B原子的分布不同，這3種不同類型的HfB晶體結構的性質將大不一樣。例如，這3類HfB晶體結構的基態能量有區別（見第2.1小節）。第Ⅰ類HfB晶體結構的基態能量最低；第Ⅱ類HfB晶體結構的基態能量稍高；第Ⅲ類HfB晶體結構的基態能量最高。下面將進一步描述這3類HfB晶體結構力學性質和電子結構的差異。

2.5 HfB的力學性質

HfB材料的一大應用是作為結構陶瓷材料，因而關注它的力學性質是十分必要的。利用計算的彈性常數，基于Voigt-Reuss-Hill近似[32-34]首先估算了多晶HfB材料的體模量B和剪切模量G。基于估算的體模量B和剪切模量G，進一步獲得了多晶HfB材料的楊氏模量E、泊松比ν、Pugh比k以及維氏硬度 Hv等。 相關計算公式為[35-36]：E=9BG/（3B+G）；ν=0.5（3B-2G）/（3B+G）；k=G/B；Hv=2（G3/B2）0.585-3。

表4所列為6種不同晶體結構類型HfB材料在基態下的力學性質。這6種HfB材料均具有較高的體模量（179～206 GPa），表明這些 HfB材料具有良好的抗壓縮性。在這6種HfB材料中，第Ⅱ類HfB材料具有最高的體模量（203～206 GPa），第Ⅰ類 HfB材料和第Ⅲ類HfB材料的體模量比第Ⅱ類的略低大約20 GPa。不同于體模量，這3類HfB材料的剪切模量存在著較大的差別。如表4所列，第Ⅱ類HfB材料具有最大的剪切模量，達到152～166 GPa；第Ⅰ類HfB材料的剪切模量次高，為115～127 GPa；第Ⅲ類HfB材料的剪切模量最低，僅為77 GPa。剪切模量間的明顯差異造成這6種HfB材料在其它力學性質上的明顯差異，尤其是楊氏模量和維氏硬度。第Ⅱ類HfB材料的楊氏模量最高（365～393 GPa），維氏硬度也最大（23.9～28.1 GPa）；第Ⅰ類HfB材料的楊氏模量較高（300 GPa左右），維氏硬度也較大（16.1～19.1 GPa）；第Ⅲ類 HfB 材料的楊氏模量最低（202 GPa），維氏硬度也最小（6.4 GPa）。材料的Pugh比常用于衡量材料的韌脆性[37]：若Pugh比大于0.57，材料為脆性材料，反之則為韌性材料。如表4所列，第Ⅰ類和第Ⅱ類HfB材料的Pugh比均大于0.57，表明它們為脆性材料；而第Ⅲ類HfB材料則為韌性材料。泊松比表示材料中共價鍵的方向性程度[38]：典型共價材料的泊松比要小于0.1，而金屬材料的泊松比一般大于0.33。以上3類HfB材料泊松比的相對大小次序為：第Ⅱ類＜第Ⅰ類＜第Ⅲ類，表明第Ⅱ類HfB材料中共價鍵的方向性最強，因而力學強度也最高。

表4 基態條件下HfB材料的力學性質:體模量B、剪切模量G、楊氏模量E、Pugh比k、泊松比ν和維氏硬度HvTable 4 Computed bulk moduli B,shear moduli G,Young′s moduli E,Pugh′s ratio,Poisson′s ratio,and Vickers hardness for HfB compounds at 0 K

2.6 HfB晶體結構的電子性質

材料的電子結構能夠幫助研究人員更為充分地了解材料的性質。通過理論計算，得到了6個HfB晶體結構的電子態密度 （包括總態密度和分波態密度）。如圖5所示，所有這6個HfB晶體結構在費米面處均具有非零的電子態密度，表明它們均具有金屬性。從分波態密度圖可以進一步得出費米面處的非零電子態主要來源于Hf的5d電子。由圖5還可看到，在費米面以下某些區間內，Hf的5d電子與B的2p電子的態密度曲線具有類似的形狀，意味著Hf的5d電子軌道與B的2p電子軌道存在雜化。這說明在這6個HfB晶體結構中可能均存在較強的Hf-B共價鍵。

為了進一步了解6個HfB晶體結構中的化學鍵，計算了它們的電子局域函數（ELF[39]）。ELF取值0～1，其數值大小可用于區分金屬鍵、共價鍵和離子鍵[40]：ELF=1表示純共價鍵或孤對電子，ELF=0.5則對應均勻的電子氣。圖6所示為6個HfB晶體結構在某些特定平面上的二維ELF圖。在第Ⅰ類HfB晶體結構和第Ⅱ類HfB晶體結構中，B原子間的區域均具有較大的ELF值，表明在這些結構中具有強的B-B共價鍵。而在第Ⅲ類HfB晶體結構中，B原子之間由于相距較遠而不成鍵。在所有3類HfB晶體結構中，Hf原子之間的電子局域分布表明Hf原子之間形成金屬鍵。此外，在這些HfB晶體結構中，Hf原子和B原子之間具有較大的ELF值，并且在偏向于B原子的區域顯示最大的ELF值，這意味著Hf-B共價鍵具有一定的離子鍵特征。

圖5 HfB晶體結構的電子態密度圖Fig.5 Density of states of HfB compounds

圖6 HfB晶體結構的二維ELF圖Fig.6 Computed 2D ELF for HfB compounds

正是由于晶體結構中具有強的B-B和Hf-B的共價鍵，第Ⅰ類和第Ⅱ類HfB材料具有高的彈性模量和維氏硬度。而由于缺少B-B共價鍵，相比于第Ⅰ類和第Ⅱ類HfB材料來說，第Ⅲ類HfB材料彈性模量和維氏硬度明顯偏低。

3 結 論

利用晶體結構預測軟件USPEX，系統搜索了HfB在基態下的穩定晶體結構。新發現2個HfB晶體結構（空間群：P6m2和R3m）。相比于文獻已報道的 4個 HfB 晶體結構（空間群：Pnma、Cmcm、I41/amd和Fm3m），P6m2結構具有最低的基態能量。計算了這6個HfB晶體結構在0～1 000 K溫度范圍內的自由能。計算結果表明：溫度能夠改變不同HfB晶體結構間能量的大小次序。特別地，當T＞777 K時，Pnma結構會取代P6m2結構成為自由能最低的HfB晶體結構。這是由于HfB晶體結構中B原子連接方式的差異性所引起的。在P6m2結構中，B原子間相互連接形成的二維類石墨烯B層造成高的聲子振動頻率，因而P6m2結構的自由能隨溫度升高而增加的量更多。而在Pnma結構中，B原子僅通過相互連接形成zig-zag形B鏈，因而Pnma結構的聲子振動頻率較低，進而造成其自由能隨溫度升高而增加的量相對更少。基于6個不同的HfB晶體結構，計算了它們的力學性質。除具有Fm3m結構的HfB材料外，其它5種HfB材料均具有較好的力學強度。考慮到Fm3m結構穩定存在的可能性較低，預測實際的HfB材料極有可能具有良好的力學強度。進一步的電子結構分析表明：造成HfB材料具有良好力學強度的原因是由于其晶體結構中具有強的B-B和Hf-B的共價鍵。本研究是材料基因組技術用于材料發現與功能性質基因分析的典型案例，對加速新型功能材料的研發與應用具有重要的借鑒價值。

[1]Oganov A R.Modern Methods of Crystal Structure Prediction.Berlin：Wiley-VCH,2010.

[2]Oganov A R,Lyakhov A O,Valle M.Acc.Chem.Res.,2011,44（3）：227-237

[3]Wang Y C,Lv J,Li Z,et al.Comput.Phys.Commun.,2012,183（10）：2063-2070

[4]Jain A,Ong S P,Hautier G,et al.APL Mater.,2013,1（1）：011002

[5]WANG Hong（汪洪）,XIANG Yong（向勇）,XIANG Xiao-Dong（項曉東）,et al.Science and Technology Review（科技導報）,2015,33（10）：13-19

[6]YE Shao-Qing（王紹青）,YE Heng-Qiang（葉恒強）.Chin.Sci.Bull.（科學通報）,2013,58（35）：3623-3632

[7]Oganov A R,Chen J H,Gatti C,et al.Nature,2009,457（7231）：863-867

[8]Li Y W,Hao J,Liu H Y,et al.Phys.Rev.Lett.,2015,115（10）：105502

[9]Zhang M,Liu H Y,Li Q,et al.Phys.Rev.Lett.,2015,114（1）：015502

[10]Duan D F,Liu Y X,Tian F B,et al.Sci.Rep.,2014,4：6968

[11]Errea I,Calandra M,Pickard C J,et al.Phys.Rev.Lett.,2015,114（15）：157004

[12]Li R H,Xie Q,Cheng X Y,et al.Phys.Rev.B,2015,92（20）：205130

[13]Zhang M G,Yan H Y,Wei Q,et al.Comp.Mater.Sci.,2013,68（2）：371-378

[14]Niu H Y,Chen X Q,Ren W J,et al.Phys.Chem.Chem.Phys.,2014,16（30）：15866-15873

[15]ZHAO Li-Ka（趙立凱）,ZHAO Er-Jun（趙二俊）,WU Zhi-Jian（武志堅）.Acta Phys.Sin.（物理學報）,2013,62（4）：383-391

[16]Ma T,Li H,Zheng X,et al.Adv.Mater.,2017,29（3）：-

[17]TAO Qiang（陶強）,MA Shuai-Ling（馬帥領）,CUI Tian（崔田）,et al.Acta Phys.Sin.（物理學報）,2017,66（3）：79-93

[18]WU Chuan（吳川）,WANG Xin（王鑫）,WU Feng（吳鋒）,et al.Modern Chemical Industry（現代化工）,2007,27（s1）：155-158

[19]Rogl P,Potter P E.Calphad,1988,12（2）：191-204

[20]ZOUWu-Zhuang（鄒武裝）.Zirconium and Hafnium Handbook（鋯·鉿手冊）.Beijing：Chemical Industry Press,2012.

[21]Xu X W,Fu K,Li L L,et al.Physica B,2013,419：105-111

[22]Failamani F,Gschl K,Reisinger G,et al.J.Phase Equilib.Diff.,2015,36（6）：620-631

[23]Huang B,Duan Y H,Hu W C,et al.Ceram.Int.,2015,41（5）：6831-6843

[24]Kohn W,Sham L J.Phys.Rev.,1965,140（4A）：A1133

[25]Kresse G,Furthmüller J.Phys.Rev.B,1996,4（16）：11169-11186

[26]Bl?chl P E.Phys.Rev.B,1994,50（24）：17953-17979

[27]Perdew J P,Burke K,Ernzerhof M.Phys.Rev.Lett.,1996,77（18）：3865-3868

[28]Togo A,Oba F,Tanaka I.Phys.Rev.B,2008,78（13）：134106

[29]Baroni S,de Gironcoli S,Dal Corso A,et al.Rev.Mod.Phys.,2001,73（2）：515-562

[30]Momma K,Izumi F.J.Appl.Crystallogr.,2011,44（6）：1272-1276

[31]Born M,Huang K.Dynamical Theory of Crystal Lattices.Oxford：Oxford University Press,1998.

[32]Voigt W.Lehrburch der Kristallphysik.Germany Leipzig：Teubner,1928.

[33]Reuss A.Z.Angew.Math.Mech.,1929,9（1）：49-58

[34]Hill R.J.Mech.Phys.Solids,1963,11（5）：357-372

[35]Hill R.Proc.Phys.Soc.London,Sect.B,1952,65（389）：396

[36]Chen X Q,Niu H Y,Li D Z,et al.Intermetallics,2011,19（9）：1275-1281

[37]Pugh S F.Philos.Mag.,1954,45（367）：823-843

[38]Haines J,Leger J M,Bocquillon G.Annu.Rev.Mater.Res.,2001,31（1）：1-23

[39]Becke A D,Edgecombe K E.J.Chem.Phys.,1990,92（9）：5397-5403

[40]Silvi B,Savin A.Nature,1994,371（6499）：683-686