低標(biāo)度XYG3雙雜化密度泛函的開(kāi)發(fā)與測(cè)評(píng)

鄭若昕，張 穎，徐 昕

（復(fù)旦大學(xué)化學(xué)系,上海市分子催化和功能材料重點(diǎn)實(shí)驗(yàn)室,物質(zhì)計(jì)算科學(xué)教育部重點(diǎn)實(shí)驗(yàn)室,上海200433）

基于密度泛函理論（DFT）和Kohn-Sham（KS）無(wú)相互作用框架下的密度泛函近似方法（DFAs）是目前使用最廣泛的第一性原理電子結(jié)構(gòu)計(jì)算方法［1］.按照Perdew等［2］提出的密度泛函發(fā)展途徑——雅各布天梯（Jacob’s Ladder），目前常用的DFAs主要隸屬于天梯的前4階近似，僅采用KS占據(jù)軌道信息來(lái)構(gòu)造近似泛函［3～25］.與之不同，第五階泛函需引入復(fù)雜的KS未占軌道的信息，被認(rèn)為是最靠近化學(xué)精度天堂的最后一階泛函近似［26～28］.

目前最成功的第五階泛函近似是雙雜化泛函（DHA）［26～28］.顧名思義，雙雜化泛函在采用Hartree-交換形式顯式地引入占據(jù)軌道信息的同時(shí)，更進(jìn)一步使用二級(jí)微擾（Second Order Energy of Perturbation Theory，PT2）相關(guān)形式引入未占軌道的信息：

式（1）為閉殼層體系的正則PT2相關(guān)能表述形式；開(kāi)殼層體系的公式相對(duì)復(fù)雜，但可由式（1）直接推廣而得.本課題組于2009年提出了XYG3雙雜化泛函［28］，具體表達(dá)式為

系統(tǒng)的測(cè)評(píng)和后續(xù)應(yīng)用證明［28～37］，XYG3型雙雜化方法可非常準(zhǔn)確地描述主族元素化學(xué)，以二階微擾波函數(shù)方法PT2的計(jì)算消耗，接近、甚至達(dá)到媲美高等級(jí)波函數(shù)方法的計(jì)算精度.

最近，我們課題組在數(shù)值基組框架下的全電子計(jì)算軟件平臺(tái)FHI-aims實(shí)現(xiàn)了低標(biāo)度XYG3方法.低標(biāo)度XYG3方法的程序基礎(chǔ)是本文作者之一——張穎博士領(lǐng)銜開(kāi)發(fā)的低標(biāo)度PT2算法［39，40］.PT2的能量計(jì)算是基于Hartree-Fock軌道和密度，且包含一份完整的PT2相關(guān)能貢獻(xiàn).與之不同，XYG3的能量計(jì)算是基于B3LYP軌道和密度，且僅包含32.11%的PT2相關(guān)能貢獻(xiàn)［見(jiàn)式（4）］.將低標(biāo)度算法應(yīng)用于XYG3方法，有可能帶來(lái)不一樣的數(shù)值收斂表現(xiàn).本文將簡(jiǎn)要介紹低標(biāo)度XYG3方法的基本原理，并開(kāi)展系統(tǒng)的測(cè)評(píng)工作.利用一系列水簇體系和ISOL22測(cè)試集證明了低標(biāo)度XYG3可以在保證計(jì)算精度的同時(shí)有效節(jié)約計(jì)算時(shí)間，減少內(nèi)存消耗并降低計(jì)算標(biāo)度，對(duì)發(fā)展適合于上千核并行的低標(biāo)度雙雜化泛函方法，并將其應(yīng)用于復(fù)雜大分子體系有重要意義.

1 公式推導(dǎo)與程序?qū)崿F(xiàn)

由式（3）可知，XYG3雙雜化泛函的計(jì)算瓶頸在于PT2部分.低標(biāo)度PT2的算法推導(dǎo)和具體實(shí)現(xiàn)見(jiàn)文獻(xiàn)［40］.簡(jiǎn)要來(lái)說(shuō)，低標(biāo)度XYG3方法繼承了H?ser等［41］提出的LT-PT2算法，通過(guò)拉普拉斯變換，將原本在分母上的特征值轉(zhuǎn)換為指數(shù)積分形式［見(jiàn)式（1）］，因此亦可稱之為L(zhǎng)T-XYG3；拉普拉斯變換需數(shù)值格點(diǎn)化，以便具體的程序?qū)崿F(xiàn).我們的LT-XYG3方法采用了Takatsuka等［42］提出的Minimax Approximation.

式中：w q和t q分別為Minimax Approximation生成的積分格點(diǎn)q的權(quán)重（weight）和根（root）.本文將討論格點(diǎn)積分?jǐn)?shù)量（N q）對(duì)低標(biāo)度XYG3計(jì)算結(jié)果的影響.由此，正則PT2表述形式［式（1）］中的分子軌道耦合可被分解，從而轉(zhuǎn)化為實(shí)空間稀疏的數(shù)值原子基表述：

式中：(su|tv)代表原子基函數(shù){φs}的雙電子積分：表示占據(jù)軌道的偽密度矩陣：

由此，式（7）可進(jìn)一步簡(jiǎn)化為

考慮到3類原子基函數(shù){φs，φ-s，φ-t}均有優(yōu)良的空間局域性，圍繞式（12），有望通過(guò)引入合適的積分篩選機(jī)制以及篩選閾值，開(kāi)發(fā)出兼具高精度和效率的低標(biāo)度XYG3方法.本工作采用由Ochsenfeld等［43］提出的顯式考慮基組對(duì)距離的QQZZR4積分篩選策略，討論不同的積分篩選閾值（θ）對(duì)低標(biāo)度XYG3計(jì)算精度的影響.然而式（12）的計(jì)算需要準(zhǔn)備兩類雙電子積分增加了相關(guān)數(shù)據(jù)存儲(chǔ)和通訊方面的需求以及并行程序?qū)崿F(xiàn)的復(fù)雜度.為應(yīng)對(duì)這一挑戰(zhàn)，我們采用局域化的Resolution-of-Identity方案（RI-LVL）［39］，將四階張量的雙電子積分(st|uv)近似分解為多個(gè)二階張量的乘積：

式中：P(S，raux)定義了以φs所在原子S為中心、半徑為raux的范圍內(nèi)的原子集合.根據(jù)這一限制，{φμ}僅包含原子中心落在上述集合內(nèi)的輔助基函數(shù).系數(shù)的確定也是采用直接最小化庫(kù)侖積分誤差的方式（具體推導(dǎo)細(xì)節(jié)見(jiàn)文獻(xiàn)［40］）.式（12）～式（14）構(gòu)筑了低標(biāo)度XYG3（LT-XYG3）程序的算法基礎(chǔ).本文將討論輔助基函數(shù)構(gòu)造參數(shù)raux對(duì)LT-XYG3計(jì)算精度的影響.

為了應(yīng)對(duì)積分篩選帶來(lái)的復(fù)雜并行應(yīng)用場(chǎng)景，LT-XYG3的程序設(shè)計(jì)跳出“數(shù)據(jù)分塊-計(jì)算-通訊”的傳統(tǒng)并行模式，而采取了OpenMP［45］與MPI［46］混合并行的框架.圖1示出了LT-XYG3中PT2部分的程序偽代碼.為了更好地介紹LT-XYG3的程序結(jié)構(gòu)，將式（12）等價(jià)變換為

具體來(lái)說(shuō)，我們將計(jì)算需求分配到不同的MPI任務(wù)中，所有的MPI進(jìn)程包含在1個(gè)或多個(gè)計(jì)算域（Domain）中，每個(gè)計(jì)算域中都分配了一整套完整的所需數(shù)據(jù)，以減少全局的數(shù)據(jù)交換.而在每個(gè)計(jì)算域中包含1個(gè)或多個(gè)對(duì)應(yīng)的MPI進(jìn)程，每個(gè)進(jìn)程中使用多個(gè)OpenMP線程并行，以共享內(nèi)存的方式減少了線程間的通訊.根據(jù)原子序號(hào)分組，兩組RI系數(shù)基函數(shù)所在原子S進(jìn)行分塊，并分配至同一個(gè)MPI進(jìn)程中.后續(xù)的積分計(jì)算在相應(yīng)的進(jìn)程中進(jìn)行線程級(jí)并行計(jì)算.

Fig.1 Pseudo code of the OpenMP/MPI hybrid parallel model for the LT-XYG3 method

2 計(jì)算細(xì)節(jié)

低標(biāo)度LT-XYG3的程序開(kāi)發(fā)在FHI-aims本地開(kāi)發(fā)版本上完成.后續(xù)計(jì)算采用數(shù)值原子基組NAO-VCC-2Z［47］.配套的輔助基組采用FHI-aims缺省算法結(jié)合NAO-VCC-2Z和額外的1個(gè)g函數(shù)生成［39，47］.所有計(jì)算均在課題組自有機(jī)群上使用4個(gè)32核節(jié)點(diǎn)（共128核）完成.具體硬件配置如表1所示.

Table 1 Hardware configuration for calculations

3 結(jié)果與討論

首先研究了拉普拉斯變換積分點(diǎn)的數(shù)量N q［見(jiàn)式（6）］、QQZZR4積分篩選閾值θ［見(jiàn)式（12）及相關(guān)討論］以及廣義RI-LVL近似下局域輔助基組的選擇半徑raux［見(jiàn)式（14）］對(duì)LT-XYG3計(jì)算精度的影響.圖2示出了20個(gè)水分子團(tuán)簇計(jì)算中LT-XYG3總能量的參數(shù)收斂行為；圖3示出了LT-XYG3計(jì)算中由于低標(biāo)度算法所帶來(lái)的近似誤差隨水簇體系尺寸的變化.我們考察了10～120個(gè)不同大小水分子簇［43］.所使用的坐標(biāo)數(shù)據(jù)由文獻(xiàn)［48］中附錄得到.

采用標(biāo)準(zhǔn)RI-V近似下的正則XYG3結(jié)果為參考標(biāo)準(zhǔn)，僅探討了LT-XYG3方法中低標(biāo)度算法引入的上述3個(gè)參數(shù)的誤差收斂行為.由圖2不難看出，LT-XYG3的計(jì)算結(jié)果對(duì)篩選閾值比較敏感.為了達(dá)到平均每個(gè)原子1 meV以下［1 meV/atom，約為0.02 kcal/mol（1 kcal=4186 J）］的計(jì)算精度，需要使用θ=10-8或者10-9的積分篩選閾值.與之相反，計(jì)算結(jié)果對(duì)格點(diǎn)積分Nq和輔助基函數(shù)選擇半徑raux收斂很快.Nq=6，raux=6a0（波爾）已經(jīng)能獲得足夠精確的LT-XYG3計(jì)算結(jié)果.

我們隨后固定3個(gè)參數(shù)的設(shè)定（raux=6a0，N q=6，θ=10-9），研究了不同大小的水簇體系.由圖3可以看出，LT-XYG3的計(jì)算精度并沒(méi)有隨著體系尺寸的增大而降低，與基于RI-V的正則XYG3方法相比，誤差始終維持在2 meV/atom之內(nèi).這表明我們基于局域RI-LVL方法分別生成的原子四中心積分和變換四中心積分在數(shù)值上是穩(wěn)定的.Neese等［49］提出局域相關(guān)方法準(zhǔn)確性的判據(jù)是要包含99.9%的總相關(guān)能.本文結(jié)果顯示，LT-XYG3方法的計(jì)算精度已經(jīng)遠(yuǎn)低于該判標(biāo).

Fig.2 Accuracy of the LT-XYG3 for a water cluster with 20 molecules

Fig.3 Absolute errors per atom in the LT-XYG3 total energies compared with the canonical XYG3 results

圖4示出了（LT）-XYG3計(jì)算中，與PT2計(jì)算瓶頸相關(guān)的總時(shí)間隨體系尺寸變化的趨勢(shì).這里考察的計(jì)算消耗主要包含RI-V或RI-LVL的系數(shù)構(gòu)建和（LT）-PT2的能量計(jì)算所花費(fèi)的時(shí)間.從圖4中可以看出，基于局域RI-LVL的LT-XYG3方法展現(xiàn)出了巨大的時(shí)間優(yōu)勢(shì).在100個(gè)水分子的體系下（這是在所使用的硬件條件下，標(biāo)準(zhǔn)XYG3計(jì)算的極限體系），低標(biāo)度算法LT-XYG3可以帶來(lái)5倍以上的加速.圖4顯示，這一加速優(yōu)勢(shì)隨著體系尺寸的增大將更加明顯，展示了LT-XYG3可應(yīng)用于大體系計(jì)算的巨大潛力.

Fig.4 Total time used for the calculation of the XYG3 correlation energy with both the canonical XYG3 and LT-XYG3 implementations

Fig.5 Memory consumption for the generation of RI coefficients in both canonical XYG3(green)and LT-XYG3(blue)

內(nèi)存消耗同樣是電子結(jié)構(gòu)算法開(kāi)發(fā)中值得關(guān)注的部分.圖5示出了2種方法在最關(guān)鍵的部分，即生成RI系數(shù)上的內(nèi)存消耗.可以看出，基于RI-V方法的標(biāo)準(zhǔn)XYG3計(jì)算中，RI系數(shù)的內(nèi)存消耗約是體系尺寸3次方的標(biāo)度；而基于RI-LVL方法的LT-XYG3由于進(jìn)一步利用了局域性質(zhì)，RI系數(shù)的內(nèi)存需求僅為體系的2次方標(biāo)度.隨著體系尺寸的進(jìn)一步增長(zhǎng)，低標(biāo)度LT-XYG3方法在數(shù)據(jù)存儲(chǔ)上的優(yōu)勢(shì)將更加明顯.

本文最后一項(xiàng)測(cè)試采用ISOL22測(cè)試集.該測(cè)試集包含22個(gè)有機(jī)分子異構(gòu)化反應(yīng)能［50］.所涉及的有機(jī)分子尺寸不一，包含24～52個(gè)不等的原子數(shù)目.以往的測(cè)試證明XYG3可以精確地描述此類化學(xué)相互作用［50］.在NAO-VCC-2Z基組下，我們采用基于標(biāo)準(zhǔn)RI-V方法的XYG3反應(yīng)能量作為參考值，計(jì)算了raux=6a0，N q=6，θ=10-9參數(shù)設(shè)置下的LT-XYG3結(jié)果.表2示出了XYG3與LT-XYG3低標(biāo)度算法在異構(gòu)化反應(yīng)能計(jì)算上的數(shù)值誤差.可以看出，所有的反應(yīng)能量誤差均低于1.7 kJ/mol，而絕對(duì)平均誤差僅為0.45 kJ/mol.可見(jiàn)，對(duì)于重要的化學(xué)反應(yīng)能的描述，引入基于RI-LVL的LT-XYG3低標(biāo)度近似算法展現(xiàn)出了令人滿意的精度.當(dāng)然，值得注意的是，LT-XYG3的數(shù)值誤差有一定的系統(tǒng)依賴性，有可能出現(xiàn)LT-XYG3的上述精度表現(xiàn)無(wú)法滿足特定模擬需求的情況.如表2所示，反應(yīng)14和16具有非常低的異構(gòu)化能絕對(duì)值.標(biāo)準(zhǔn)XYG3方法給出的異構(gòu)化能分別為-0.59 kJ/mol和1.85 kJ/mol.LT-XYG3低標(biāo)度近似引入的誤差在反應(yīng)14上僅有-0.01 kJ/mol（相對(duì)誤差也僅為1.7%）.然而在反應(yīng)16上，雖然LT-XYG3的絕對(duì)誤差只有0.86 kJ/mol，但是其相對(duì)誤差超過(guò)了總異構(gòu)化能的45%.在實(shí)際研究中若碰到類似情況，需要通過(guò)提高上述介紹的3個(gè)與低標(biāo)度算法精度相關(guān)的參數(shù)設(shè)定來(lái)達(dá)到更高的精度需求.

Table 2 XYG3 isomerization energy for the ISOL22 test set

4 結(jié) 論

XYG3雙雜化密度泛函方法通過(guò)引入PT2形式的相關(guān)能信息，在計(jì)算分子能量性質(zhì)上表現(xiàn)優(yōu)秀，然而在計(jì)算精度提高的同時(shí)也帶來(lái)了與PT2相當(dāng)?shù)挠?jì)算消耗［O(N5)］.為了進(jìn)一步增強(qiáng)雙雜化泛函方法的實(shí)用性，使其能夠應(yīng)對(duì)復(fù)雜大體系的研究需求，我們程序化實(shí)現(xiàn)了低標(biāo)度XYG3方法（LT-XYG3），通過(guò)結(jié)合局域RI方法（RI-LVL）和拉普拉斯變換PT2（LT-PT2）有效降低了雙雜化泛函的計(jì)算消耗.在該方法中，拉普拉斯變換將PT2的標(biāo)準(zhǔn)離域分子軌道表述轉(zhuǎn)化為更加定域的原子基表述，使得引入積分篩選來(lái)降低計(jì)算量的低標(biāo)度算法成為可能.另一方面，我們采用RI-LVL以及輔助基函數(shù)的局域性質(zhì)進(jìn)一步有效降低了計(jì)算需求與內(nèi)存消耗.然而這一低標(biāo)度算法也同時(shí)帶來(lái)了不規(guī)則的復(fù)雜并行應(yīng)用場(chǎng)景，因此無(wú)法通過(guò)正常分塊的方式對(duì)算法做并行處理.LT-XYG3程序采用MPI與OpenMP相結(jié)合的混合并行機(jī)制，其在被應(yīng)用于不同尺寸的水簇以及ISOL22測(cè)試集的計(jì)算時(shí)，可以在保持極高的計(jì)算精度的同時(shí)顯著降低內(nèi)存消耗并提升計(jì)算效率，證明了可兼顧高精度和高效率的低標(biāo)度雙雜化泛函方法在有機(jī)大分子計(jì)算上的可靠性和巨大應(yīng)用前景.