張廣迪 毛力2)? 徐紅星2)3)4)?

1) (武漢大學物理科學與技術學院,武漢 430070)

2) (武漢量子技術研究院,武漢 430205)

3) (武漢大學微電子學院,武漢 430072)

4) (河南省科學院,鄭州 450046)

通過將平面波與高斯函數或者樣條函數結合到一起,本文構建了一種新的復合基組.利用格拉姆-施密特正交化方法或者L?wdin 正交化方法,對復合基組進行正交歸一化.通過選擇平面波函數中波矢的絕對值,選擇性地求解某個能量區(qū)間內的本征態(tài),將大型哈密頓矩陣的計算轉變?yōu)槎鄠€小型矩陣的計算,以及通過減少電子勢能平緩部分展開基矢數目,極大地加快了計算速度.以一維有限深勢阱為例,通過與嚴格計算方法的對比,驗證了本文復合基組能夠在加速計算的情況下保證求解精度.同時,本文還研究了不同的參數設置對計算精度的影響,包括復合基矢的疏密度、高斯函數的寬值,以及樣條函數不同區(qū)域占函數總寬度的比值等參數.最后該復合基組可以直接應用到對大尺寸納米金屬結構的等離激元數值計算當中.

1 引言

金屬微納結構(納米金屬棒、顆?；虬?間的納米間隙具有很重要的性質,在物理、化學和生物等領域都有著廣泛的研究與應用[1,2],比如在表面等離激元方面的應用[3,4].表面等離激元是電磁場與固體中自由電子集體振蕩強耦合形成的一種元激發(fā),它可以將光場局域在突破衍射極限的亞波長范圍內.例如金屬顆粒之間的納米間隙[5-10],人們稱之為“熱點”.其中的電場強度被等離激元共振增強,增強倍數甚至可以達到102,從而使得間隙內部的分子拉曼散射信號極大增強(由于拉曼信號正比于電場的4 次方,所以最大可達108),這就是表面增強拉曼光譜[11-13].一般而言間隙間距越小電場增強越大,但是間隙小到納米尺度時量子效應會起作用[14-17],以往僅利用經典的電磁學方法計算電磁場[18-20]會失效.此時可以使用密度泛函理論進行計算[21,22],以獲得更準確的結果.

密度泛函理論(DFT)[23-26]以及含時密度泛函理論(TDDFT)[27,28]是研究多電子體系的重要方法,在物理和化學的很多領域都有著廣泛的應用.密度泛函理論基于Hohenberg-Kohn 定理: 給定一個在特定初始狀態(tài)下的電子系統(tǒng),外部的勢能和電子密度之間存在著一一對應關系.相互作用電子系統(tǒng)的基態(tài)密度和勢場分布可以由對應的無相互作用系統(tǒng)的Kohn-Sham 方程組的解確定.因此,密度泛函理論的關鍵在于快速準確地求解Kohn-Sham 方程組[29,30].然而,在實際應用中,很多器件總體尺寸在微米量級.如果使用密度泛函的方法,為了使計算結果足夠準確,需要使用足夠多的波函數進行展開,Kohn-Sham 方程組中的哈密頓量的矩陣維度將會非常大.以直徑10 nm 長1 μm 的金納米棒為例,電子數目將達到百萬級,由于需要計算哈密頓量的矩陣元,以及對矩陣進行對角化,矩陣維度過大,此時計算量將會變得非常大(對角化矩陣的計算復雜度為O(N3) ,同時需要計算N2量級的勢能積分).這樣會給計算資源帶來很大壓力,甚至可能因為矩陣維度過大而無法計算.如果直接減少波函數的數目來加快計算速度,將會導致計算結果的誤差較大.因此,急需新的算法來加快大尺度密度泛函理論的計算速度.

減小矩陣的維度有不同的方法,例如Meng等[31,32]發(fā)展了混合量子力學和經典電磁學(QM/EM)方法,將系統(tǒng)依據特征劃分為量子區(qū)間與經典區(qū)間,在經典區(qū)間用時域麥克斯韋方程組求解.在量子力學區(qū)間用量子非平衡格林函數求解,量子力學和經典區(qū)間界面處的瞬態(tài)電勢分布和電流密度分別用作量子力學和電磁模擬的邊界條件.由于此時量子區(qū)間非常小,可以有效地減少矩陣維度.除此以外還可以通過基組的選取減小矩陣維度.在DFT的計算中,選取合適的基組,是一個非常重要的問題,合適的基組可以使得計算準確與快速.很多研究都在關注基組的選取,例如何禹等[33]通過研究基函數效應,找到了DFT 理論中B972-PFD 方法的最具有性價比的基組.另外將不同的基矢復合到一起,作為新的基矢對研究的系統(tǒng)進行展開,有時可以得到比單一基矢更好的結果.段宜武等[34]發(fā)展了少體物理中的簡諧振子展開方法,引入了低能態(tài)波函數并將其與簡諧振子波函數混合作為基矢,利用初步的近似波函數,與一些沒有使用過的簡諧振子乘積波函數乘到一起,構成新的復合基矢,從而減小基矢總數.本文將使用復合基組的方法,將平面波與高斯函數或者樣條函數乘在一起來減小密度泛函理論中矩陣的維度.

在第一性原理計算中,平面波基組應用非常廣泛(例如著名的VASP 軟件)[35],它在處理周期系統(tǒng)時非常有效.然而,在開放系統(tǒng)中,體系尺寸往往較大,有的地方需要精確的密集數據,有些地方只需要少量的稀疏數據就可以得到較精確的結果.例如在求解銀納米棒的表面等離激元激發(fā)時,銀納米棒的兩端以及表面電子勢能變化快,往往需要精確計算.而納米棒“棒芯”部分由于電子勢能平緩,原則上僅需少量信息就可以描述.因此并不需要花費過多的計算資源在占比極大的中間部分.作為一種非局域基矢,平面波基組需要大量的平面波來對物理空間的所有部分進行展開,包括對我們的問題不重要的中間區(qū)域.因此選擇線性復雜度的局域基組(例如高斯基組[36]),可以根據需要調整不同區(qū)域的基矢疏密度,減少所用的基矢總數,從而加快計算.盡管如此,高斯基組間的間距也是原子尺度,在介觀尺度下需要的基組數目仍然非常龐大.本文基于這一點,發(fā)展了一種新型的高斯函數或含有多項式的樣條函數與平面波復合的基組.這種基組同時利用高斯函數或者樣條函數的局域特性,通過選擇不同平面波波矢來控制求解能量區(qū)間,實現了對哈密頓矩陣的分區(qū)求解.在每一個小區(qū)間可以顯著地減少使用的基矢數,在保持相應精確度的情況下實現最終的加速計算.

2 計算方法

2.1 初始復合基組

本文使用高斯函數或樣條函數與平面波函數相乘,得到一組新基矢:

或樣條函數,這里σ為高斯函數的寬值,xn為不同高斯函數或樣條函數中心的位置,樣條函數參數cij和u0由寬度參數al計算得到,高斯函數或樣條函數的疏密度可以方便地進行調整.使用的高斯函數與樣條函數如圖1 所示.

圖1 用來復合的高斯函數與樣條函數 (a)高斯函數;(b)樣條函數Fig.1.Gaussian wave function and spline function for composition: (a) Gaussian wave function;(b) spline function.

可以通過將平面波函數與高斯函數或樣條函數簇相乘得到一系列未正交歸一化新基矢,即(1)式.由完備性考慮平面波波矢k應該同時選取正負某個值,通過改變k的絕對值可以選擇求解的本征值區(qū)間.特定k值的平面波相當于一個在勢能平緩區(qū)間的試探解,算法自動組合各局域高斯或樣條函數來匹配試探解與嚴格解之間的k值差距.

在量子力學的相關計算中,基組的選取與構建需要遵循一定的原則.可以看出,本文的復合基組有以下特點: 首先,由基矢的傅里葉變換可知,基組在動量空間中±k附近是近似完備的;其次,基組是正交歸一化的,2.2 節(jié)將進行基組的正交歸一化處理;然后,基組中基矢的形式在數學上是便于積分的,高斯函數或樣條函數,以及它們與平面波函數的乘積,均可以求出解析的積分結果;最后,由于復合基矢的形式中含有平面波,基矢的形式與系統(tǒng)波函數有一定的類似性.在滿足上述原則的情況下,本文的復合基組可以很好地應用于具有電勢邊界變化劇烈而內部平緩這一特征的物理體系.另外,由于復合基矢同時包含空間局域化函數(高斯函數或樣條函數)與空間擴展函數(平面波函數),便于通過改變高斯函數或樣條函數的空間分布疏密度,以及改變平面波波矢的選取這兩種途徑,來加速本征值的求解.

2.2 復合基組的正交歸一化

新基矢的數目需要根據系統(tǒng)的長度和精度要求進行調整.在進行計算前還需要對它們進行正交歸一化.本文使用格拉姆-施密特正交化方法[37,38],或者L?wdin 的正交化方法[39],得到一系列新的正交歸一化基矢.

格拉姆-施密特正交化方法的基本流程是,選第一個原基矢作為第一個新基矢,不斷構建新的基矢,使每個新基矢與前面所有的新基矢正交,并且對新基矢進行歸一化.格拉姆-施密特正交化方法可以用如下公式表示:

其中ψn為正交歸一后的基矢,?n為原基矢,ψ1=?1,由表達式可知其是一個迭代方法,直接用(4)式一般不太穩(wěn)定,研究者們提出了許多改進方法,格拉姆-施密特正交化本質上等價于對矩陣做QR 分解,可以直接調用現在非常成熟的線性代數包求正交矩陣Q和系數矩陣R來計算.

L?wdin 正交化方法通過構建一個變換矩陣來把原基矢進行正交歸一化:

Si為(6)式中矩陣的第i個特征值,可以得到變換矩陣:

本文使用上述兩種正交化方法對復合基組進行正交歸一化,并進行比較,可以發(fā)現,兩種方法均可以對復合基組進行足夠精確的正交歸一化處理.然而由于格拉姆-施密特正交化方法需要不斷地迭代,當兩個基矢相距較近時,會使得分母接近于零從而導致誤差的積累,需要進行更細致繁瑣的參數調整來進行處理.因此本文使用L?wdin 正交化方法,或者調用科學計算庫中(例如GSL)的QR分解(以代替格拉姆-施密特正交化方法),來進行正交歸一化處理.

3 結果與討論

3.1 一維有限深勢阱的能級分布

利用復合基矢計算時,某些情況下有嚴格的解析結果,可以通過對比誤差來分析它的有效性.本文利用復合基矢計算一維有限深勢阱的能級分布與波函數.一維有限深勢阱為

該系統(tǒng)有解析解[40].下面用前文得到的基矢對一維有限深勢阱進行展開,求得能量值后與解析結果進行對比.

系統(tǒng)的哈密頓量為

其中m為電子質量,? 為約化普朗克常數.用本文的復合基矢可以求得哈密頓矩陣的矩陣元為Hmn=〈ψm|H|ψn〉,將(2)式或(3)式中的函數代入其中并數值對角化即可求得能量本征值.這些計算值與理論值之間的相對誤差與系統(tǒng)的長度、使用的基矢的數目、基矢結構參數均有關系.在圖2 中橫坐標為不同的本征值,縱坐標為計算值與理論值之間的相對誤差,我們僅使用少量基矢,發(fā)現即使在這種情況下面相對誤差也很小.但是圖2 可以說明本文方法精確度非常高.計算的相對誤差基本上都在10-5和 10-6的量級.可以看出,減小中間基矢的數目,從而減小總基矢的數目,仍然可以得到足夠精確的結果,同時又可以明顯加快計算速度.當系統(tǒng)的物理尺寸較大,因而需要較多的基矢總數時,3.3 節(jié)將通過調整基矢中平面波波矢的值來提高計算的效率.

圖2 僅使用少量基矢時,計算值與理論值的相對誤差(x 軸為不同的本征態(tài),y 軸為計算值和理論值之間的相對誤差,后文相同)Fig.2.Relative error between calculated and theoretical values when only a small number of base vectors are used(x-axis,different eigenstates;y-axis,relative error between calculated and theoretical values;the following image is the same as the Fig.2).

3.2 不同的函數設定對精度的影響

由于高斯函數或樣條函數均與坐標有關,并且可以改變高斯函數的寬值,或者分段函數各段的長度,因此,本文研究了不同的參數對新的復合基矢精度的影響.以高斯函數構成的復合基矢為例.由于本節(jié)使用了較多的基矢,計算的相對誤差與2.1 節(jié)相比會有略微的不同.

首先考慮不同的疏密度比率對精度的影響.如圖3 所示,橫坐標為不同的本征態(tài)計數,縱坐標為計算值與理論值之間的相對誤差,r為疏密度比率,疏密度比率是單位長度內系統(tǒng)兩端(即電子勢能變化快速部分)使用的基矢數與系統(tǒng)中間部分(電子勢能平緩部分)使用的基矢數的比值.從圖3 可以看出,即使在疏密度比率等于1 時,也就是系統(tǒng)兩端部分和中間部分的基矢密度相等時,復合基矢仍然保持著較高的精度,所求的系統(tǒng)本征值與理論本征值之間的相對誤差的數量級在 10-4的位置.而隨著疏密度比率逐漸增大,相對誤差逐漸減小,直到 10-5的位置.由于系統(tǒng)中間部分需要最低數量的基矢來展開,因此隨著疏密度比率的增大,相對誤差不會無限減小,在本文對應的程序里面可以根據不同的系統(tǒng)求出一個相對誤差接近最小的疏密度比率.綜上,通過調整基矢疏密度確實能夠提高精度,圖3 顯示大約可以提高一個量級.

圖3 不同的疏密度比率對精度的影響Fig.3.Effect of different degree of closing on precision.

接著討論不同的基矢寬度對精度的影響.圖4橫坐標為不同的本征態(tài)計數,縱坐標為利用復合基矢求得的本征值與理論本征值之間的相對誤差,s為基矢寬度與基矢之間距離的比值,圖中疏密度比率為5.可以看出,整體的相對誤差數量級在10-5的位置.相對誤差隨著基矢寬度變寬先減小然后增大,在疏密度比率為5 的條件下,最優(yōu)基矢寬度大概為0.85 相鄰基矢間距.

圖4 不同的基矢寬度對精度的影響Fig.4.Effect of different basis sets width on precision.

對比圖3 和圖4 可以看出,基矢的不同疏密度比率對精度的影響,比不同的基矢寬度對精度的影響更大,疏密度比率是更具決定性的影響精度的因素.

3.1 與3.2 節(jié)僅計算了能量最低(平面波波矢k的取值為0,即不乘平面波只用局域基)的少量本征值.要想精確求解大量本征值,傳統(tǒng)方法需要加大展開基矢數目,但是這樣會顯著增加計算量.本文方法可以通過選取平面波波矢的絕對值來計算特定能量區(qū)間的本征態(tài).同時選擇樣條函數能夠進一步提高計算精度.

3.3 基組中平面波波矢的值對計算效率的影響

以1 μm 厚度的銀納米板為例,勢阱深度大約10 eV,費米能5.5 eV,銀納米板對應在厚度方向大概3800 條能級被占據,如果用以往方法大約需要數萬數量的基矢,才可得到比較精確的計算結果.利用本文方法可以極大減小這個數值.例如求解費米面附近的電子時,在勢阱里面的區(qū)域選擇201 個樣條函數,勢阱外選擇60 個,由于平面波基矢用到了正負的k,這樣總共的基矢空間為522(為了進一步減小基矢空間,也可以在勢阱外不再乘以平面波,也就是只使用樣條函數,此時基矢空間為462).通過對角化哈密頓矩陣,得到522 個能量本征值.由于選擇的基矢數目非常少,只有能量在 ?2k2/(2m) 附近非常有限的能量本征值準確(能夠算準的個數大約為1/10 中間樣條函數個數).在上述條件下總共需要約190 次這樣的計算.于是這種方法就相當于把對大型矩陣的計算,轉變?yōu)閷υS多小型矩陣的計算,從而提高計算效率.通過計算能量在 ?2k2/(2m) 附近波函數與坐標軸交點個數,可以判斷出它們對應于哪個嚴格解,并且通過對比它們的本征值與嚴格解,也確認了該挑選方法是準確的.

接下來計算這種方法的相對誤差,并分析不同的參數設置對相對誤差的影響.勢阱不同區(qū)域樣條函數疏密度的比值對精度有很大影響.由于勢阱中間部分勢能平緩,原則上講少數樣條函數足以描述;但是在勢阱壁勢能變化劇烈,需要使用較密集的函數.在基矢總數保持不變的條件下,兩端邊緣部分與中間部分的樣條函數疏密度比率對精度影響很大.圖5 為疏密度比率對精度的影響,可以看出r為11.94 時計算精度最高.

圖5 不同的疏密度比率對精度的影響Fig.5.Effect of different degree of closing on precision.

樣條函數不同區(qū)域與整個函數的寬度比值對計算精度也有一定影響.如(3)式和圖1 所示,樣條函數包括兩邊的多項式區(qū)域和中間的直線區(qū)域,為了計算方便,可以令樣條函數的多項式區(qū)域,正好與兩個樣條函數相互交疊的區(qū)域完全重合.在圖6 中,s為某個樣條函數與下一個樣條函數相互交疊的部分,與該樣條函數間距的比值.可以看到每次計算也就對應著每選取一個波矢,有少數能量本征值(例如圖5 中的綠線,大約有中間樣條函數個數除10,即20 個)可以被準確地計算.耦合部分與函數間距的比值對計算精度有明顯的影響,原則上選擇的交疊部分寬度應該比(1)式中的平面波波矢所對應的電子波長小一個量級.至于具體哪個值最優(yōu)可以通過試探、擬合的方式找到.

圖6 耦合部分與函數寬度的比值對精度的影響Fig.6.Effect of the ratio of coupling part to function width on precision.

對比圖5 和圖6 可以發(fā)現,與3.2 節(jié)類似,與耦合部分占函數寬帶的比值相比,基矢的疏密度比率是更具決定性影響精度的因素.

最后估算該方法的計算效率相的提升.本文的等離激元數值計算中,最費時的地方一般是哈密頓矩陣的矩陣元與對角化計算.對于矩陣元的計算,需要求解N2個積分.對于矩陣的對角化來說,其計算復雜度為O(N3) .因此基矢總數的減少對計算速度的影響非常大,如果使用通常的平面波基矢對上述系統(tǒng)進行展開,大約需要數萬的基矢,則計算量達到了約 1013的量級.而使用本小節(jié)的方法矩陣對角化的計算量約為 190×5223,大約在 1010的量級,在保持較好的計算精度的同時,計算效率大大提高.

4 結論

本文利用平面波函數與高斯函數或者樣條函數相乘,構建了一種新的基組,并利用格拉姆-施密特正交化方法和L?wdin 正交化方法,對構建的基矢進行正交歸一化.通過選擇平面波函數中波矢的絕對值,選擇性地求解某個能量區(qū)間內的本征態(tài),將大型哈密頓矩陣的計算轉變?yōu)槎鄠€小型矩陣的計算,以及通過減少電子勢能平緩部分展開基矢數目,極大地加快計算速度.為了方便研究算法精度本文選擇一維有限深勢阱,發(fā)現新的復合基組在保持精度的前提下,可以顯著提高計算效率.同時研究了基矢寬度和疏密度比率對精確度的影響,探討了怎樣選擇最優(yōu)的參數.

一維有限深勢阱勢能在邊界跳變,在其他地方為常數,它和金屬電勢具有類似的變化規(guī)律.以納米金屬平板為例,電勢在一個費米波矢內從零變化到最小,而后在幾個費米波矢內邊振蕩邊衰減到平均值[41].上述計算方法在一維有限深勢阱效果很好,因此可以推論該方法能夠應用到銀納米板或者其他金屬納米結構等離激元DFT 計算中.給定一個合理的試探性的初始態(tài),利用DFT 理論,通過自洽求解,可以求解系統(tǒng)的基態(tài)電子密度分布,之后利用TDDFT 理論求解系統(tǒng)電磁場分布以及光學等性質.