何長春 廖繼海 楊小寶

(華南理工大學物理與光電學院,廣州 510640)

平面團簇穩(wěn)定結構的蒙特卡羅樹搜索?

何長春 廖繼海 楊小寶?

(華南理工大學物理與光電學院,廣州 510640)

(2017年4月7日收到;2017年6月18日收到修改稿)

以平面團簇為例提出了一種結合結構識別和蒙特卡羅樹技術搜索穩(wěn)定結構的新方法.體系原子之間的相互作用由兩類模型勢能函數(shù)來描述:Lennard-Jones二體勢函數(shù)與基于Lennard-Jones勢的三體勢函數(shù).考慮可能的三角晶格碎片作為候選結構,引入編號策略對結構進行快速識別,并運用蒙特卡羅樹搜索研究穩(wěn)定結構隨著原子數(shù)增大的演化過程;對于能量較低的候選結構,進一步采取局域優(yōu)化來獲得對應體系的穩(wěn)定結構.計算表明,Lennard-Jones二體勢函數(shù)對應的三角晶格團簇更穩(wěn)定;在特定的參數(shù)下,三體勢函數(shù)對應的六角晶格團簇更穩(wěn)定.結合結構識別和蒙特卡羅樹搜索可以對候選結構空間進行高效掃描,在較短時間內(nèi)更容易搜索到穩(wěn)定的團簇結構,并可以與第一原理計算結合實現(xiàn)材料的結構預測.

團簇,結構識別,蒙特卡羅樹,全局優(yōu)化

1 引 言

團簇一般由從幾個到上千個原子組成,是宏觀世界與微觀世界的連接紐帶,結構豐富并具有很多奇特的性質(zhì),在物理、化學、材料等領域具有重要的研究價值[1?4].幻數(shù)(magic number)是團簇的重要特征之一[5],而幻數(shù)結構的確定在原子團簇的研究中具有相當重要的作用[6].近些年來,各種各樣的原子或分子團簇陸續(xù)在實驗中被制備出來,也引起了理論工作者的廣泛關注[7].例如,實驗中發(fā)現(xiàn)惰性氣體的團簇往往呈現(xiàn)出高對稱結構[8];理論研究指出,金屬團簇在尺寸較小的時候是高對稱結構,而Ru55,Rh55,Pd55等的穩(wěn)定結構不是高對稱性的正二十面體結構,而是面心立方結構對應的晶格碎片[9].實驗中制備了不同尺寸的石墨烯量子點,發(fā)光頻率覆蓋所有可見光范圍,其對應的結構是平面六角結構的晶格團簇[10];在從石油中提取和分離出來的納米金剛石團簇,也是金剛石結構的晶格碎片[11,12].有趣的是,體相的硼是由正二十面體堆積起來的結構,而硼團簇傾向于形成平面三角結構,隨著原子數(shù)的增大體系將出現(xiàn)空位,出現(xiàn)由原子和空位形成的“合金”結構[13,14].

理論對團簇的結構和性質(zhì)研究也有很多相關的工作,其中包括對結構能量的評價以及不同結構的變換和篩選.一般來說,第一原理計算對能量的評價精度高,但比較耗時,同時局域優(yōu)化方法(準牛頓法、共軛梯度法等)對初始結構依賴嚴重.在做結構搜索時,我們通常需要借助勢能函數(shù)簡化能量的計算來提高搜索效率;然而,即使在給定的勢能函數(shù)的條件下,尋找原子團簇的最優(yōu)結構在理論上被證明為NP-Hard問題,通常需要利用啟發(fā)式智能算法進行優(yōu)化計算,使得候選結構中不陷入局域最小值,實現(xiàn)全局最優(yōu)的結構搜索,常見的方法包括遺傳算法[15]、粒子群算法[16]、蒙特卡羅模擬退火算法[17]、分子動力學模擬[18,19]等.邵桂芳等[20]選擇經(jīng)驗勢能來描述原子間相互作用,運用改進的Basin-Hopping Monte Carlo算法,研究了不同尺寸和不同比例下的Fe-Pt二元合金團簇.張林等[21]結合遺傳算法和基于密度泛函理論的緊束縛方法(density functional tight binding,DFTB),搜索發(fā)現(xiàn)SimGe9?m團簇存在的兩種低能原子堆積穩(wěn)定構型.潘必才等[22]運用兩種結構搜索方法(壓縮液態(tài)法、遺傳算法)與鍺的緊束縛模型勢相互結合,確定了Ge65,Ge70,Ge75的能量較低的可能結構,發(fā)現(xiàn)這三種團簇均具有兩類穩(wěn)定、能量相近的異構體:類球形和類橢球形.值得注意的是,在結構搜索過程中,如果不進行結構檢查容易導致重復計算而降低搜索效率.基于微粒群和遺傳算法等全局優(yōu)化方法,兩類常用的搜索軟件Calypso[23]和Uspex[24]都提出了表征結構相似性的算法,可以實現(xiàn)高效的結構搜索.

理論上要確定團簇的幻數(shù)結構,我們需要在給定原子數(shù)的條件下搜索能量最低結構,再根據(jù)能量隨原子數(shù)變化曲線的局域極值來實現(xiàn).Greiner等[25]選擇Lennard-Jones勢描述原子之間的相互作用,運用團簇生長模型,在 N 個原子的最穩(wěn)定團簇基礎上尋找 N+1個原子團簇的最穩(wěn)定結構,被證明是一種效率較高的搜索方法.以平面團簇的結構搜索為例,我們用Lennard-Jones勢函數(shù)描述原子之間相互作用.我們選擇三角晶格團簇作為初始結構,引入編號方法對結構進行識別檢查,避免相同結構重復的能量評價,可以很好地提高搜索效率.考慮到不同大小的晶格團簇互相聯(lián)系,我們采用蒙特卡羅樹技術[17],以低能量為目標進行搜索.該方法可以同時兼顧到搜索的深度和廣度,可以很好地增加樣本的多樣性,提高低能量結構搜索的成功率.在找到能量較低的結構后,我們進一步對結構進行局域優(yōu)化,獲得體系的最穩(wěn)定結構.

2 計算方法

我們的搜索方法主要包括結構識別和蒙特卡羅樹技術.選擇的初始結構是晶格團簇,可以根據(jù)對稱性嚴格檢查兩個結構是否全同;蒙特卡羅樹技術可以建立不同原子數(shù)構型之間的聯(lián)系,同時根據(jù)體系能量的高低對不同分枝進行重要性抽樣.

2.1 結構識別

使用坐標描述給定團簇的構型是簡單易行的方法,但在結構識別時很不方便.以平面三角晶格團簇為例,我們按照編號的方法對結構命名,可以實現(xiàn)快速的結構檢查.首先,我們在三角晶格平面上選擇一個格點為原點,根據(jù)其他格點到原點距離由小到大進行編號:對于距離相同的格點,我們從極角最小者開始,按照逆時針方向從小到大進行編號,格點的編號如圖1(a)所示.接著,對于給定的三角晶格團簇,我們可以在已經(jīng)有編號的平面上,通過各個原子的坐標分別找到對應的編號,并從小到大排序,獲得一組編號表示該結構(注意:此時等價的結構有多個可能的編號).最后,我們考慮可能的平移、反演、旋轉(zhuǎn)等平面三角晶格允許的對稱操作,并記錄所有可能的編號數(shù)組,找出最小的編號作為該結構對應的編號(我們從第一個分量開始逐次比較大小,如果第一個的分量相同就考慮下一個分量,以此類推,比如[1-2-3]小于[1-3-4]),此時獲得的編號和結構一一對應.結構檢查時只需要比較兩個編號是否完全相同即可,從而有效提高了程序的運行效率.當原子數(shù)為4時,共有七種同分異構體,其結構以及對應的編號如圖1(b)所示.

圖1 (網(wǎng)刊彩色)(a)三角晶格背景格點編號;(b)原子數(shù)為4時的七種同分異構體及其對應的編號Fig.1.(color on line)(a)The nuMbered triangu lar lattice;(b)seven isoMers With 4 atoMs and the corresponding nuMber series.

2.2 蒙特卡羅樹搜索技術

對于平面三角晶格團簇,我們根據(jù)編號方法可以很快遍歷得到原子數(shù)不大于10的全部結構(其中原子數(shù)是10時對應的結構數(shù)為30490),并根據(jù)勢能函數(shù)計算獲得體系的最優(yōu)結構.隨著原子數(shù)的增多,可能的候選結構數(shù)目急劇上升,這要求我們只能選擇有限數(shù)量結構來繼續(xù)增加原子.然而,簡單選擇一些能量較低的結構,容易讓體系陷入局域最優(yōu);過多選擇結構,在有限的時間內(nèi)很難得到較為穩(wěn)定的結構.因此,我們借鑒了用于處理連續(xù)決策優(yōu)化問題的上限信心策略(upper confidence bound).該算法最初是為了研究多臂賭徒模型(mu lti-armed bandit Model)[26,27],核心在于計算其信心索引Ij:

本文算法的流程圖如圖2(a)所示:我們根據(jù)初始的原子數(shù)為N的結構,按照能量的大小,依照玻爾茲曼分布進行重要性抽樣得到一定數(shù)目的樣本,樣本j被抽到的概率Pj表示為

這里Ej是結構樣本中第j個結構的能量;k是玻爾茲曼常數(shù);T是溫度,用來調(diào)節(jié)搜索的廣度與深度.在這些樣本上增加一個原子拓展得到N+1個原子的結構,最后對能量最低的100個候選結構進行局域優(yōu)化.圖2(b)給出原子數(shù)為3的結構通過選擇-拓展-選擇得到原子數(shù)為4的結構的示例.首先,原子數(shù)N=3時有三個不等價結構,根據(jù)(2)式進行重要性抽樣得到兩個結構,然后對這兩個結構進行完全拓展得到原子數(shù)N=4的六個結構,再根據(jù)(2)式進行重要性抽樣得到三個結構,實際操作中我們選取每次抽樣的樣本數(shù)目隨著原子數(shù)的增加而相應增加,重復進行上述操作得到原子數(shù)更多的結構.

圖2 (a)搜索算法流程圖;(b)原子數(shù)從3到4的操作舉例Fig.2.(a)The fl oWchart of the search algorithm;(b)the exaMp le for the case froMN=3 to N=4.

相鄰原子數(shù)對應的結構互相關聯(lián),我們可以得到類似“樹”一樣的結構,在10個原子以內(nèi)通過遍歷得到的是完整的“樹”.之后,我們根據(jù)能量的高低進行了剪枝,使得“樹”出現(xiàn)非對稱性隨機生長,傾向于保留能量較低的結構,可以有效縮小范圍并提高搜索效率.這里,我們在樹的每一層以一定概率選取初始的父結構,然后在這些結構中進行拓展計算出能量最低的結構,并根據(jù)概率選擇作為下一層的初始結構,并通過循環(huán)操作不斷增加體系的原子數(shù).其中概率是按照玻爾茲曼分布設置的,如(2)式所示.可以看到,能量越低的結構被選中的概率越大,但是能量較高的結構仍有被選中的概率,可以保證樣本的多樣性,避免陷入局域極小值.(2)式中的T類似于模擬退火算法中的溫度,溫度越小則全局最小的結構被搜索出來的概率越高.當T趨于零時,該算法將趨于完全的深度搜索;T趨于無窮大時,該算法將趨于完全的廣度搜索.

3 結果與討論

為了檢驗算法的有效性,我們選擇了兩類勢能函數(shù)描述原子之間的相互作用,并搜索了給定勢函數(shù)下不同原子數(shù)的穩(wěn)定結構.

3.1 Lennard-Jones勢

首先,我們考慮體系中粒子之間的相互作用符合Lennard-Jones勢函數(shù),則體系的總能量V為

其中rij為原子i,j之間的距離.圖3給出了我們搜索得到的最低的平均能量趨勢圖以及幻數(shù)所對應的穩(wěn)定結構.在該勢能函數(shù)下,原子傾向于形成盡可能多的最近鄰,在平面結構中對應的是六個最近鄰,同時要求團簇的邊界呈凸包絡,可以使得邊緣的懸掛鍵數(shù)目比例最低,增加體系穩(wěn)定性.圖4給出了原子數(shù)較多時對應的幻數(shù)結構.可以看出,幻數(shù)所對應的結構往往具有較高的對稱性,邊界對應的是六邊形;隨著原子數(shù)的增多,穩(wěn)定結構包含了所有的正六邊形的結構,例如原子數(shù)為19,37,61,91時,分別對應于邊長為2,3,4,5的正六邊形.

圖3中用×表示的即為在搜索空間下所訪問過的結構的平均原子能量,我們稱之為候選結構.這里在每一次的循環(huán)操作中將會以(2)式的概率隨機選取一定數(shù)量的結構進行拓展,以生成N+1個原子的結構,在N≤100時一共計算了8×105個結構.將給定原子數(shù)的最低能量的平均值連接起來即為圖3中的折線.

圖3 (網(wǎng)刊彩色)Lennard-Jones勢函數(shù)下原子團簇平均能量的變化趨勢Fig.3. (color on line)The average energy of clusters as a function of the nuMber of atoMs under the Lennard-Jones potential.

圖4 (網(wǎng)刊彩色)Lennard-Jones勢函數(shù)下的幻數(shù)Fig.4.(color on line)TheMagic nuMber structures.

利用蒙特卡羅樹搜索,我們可以得到能量較低的三角晶格團簇,屬于平面密堆結構,與Yang和Zhang[29]之前工作相符合.在本文的算法中,35個原子以內(nèi)的體系通常需要計算體系的勢函數(shù)1000—3000次就可以搜索到基態(tài)結構,而36—50個原子的體系則需要評價3000—5000次,但是一般的微粒群算法需要對結構進行8000次的評價才可以搜索出基態(tài)結構,可見本文的算法在結構搜索的效率上有了一定的提升.我們可以進一步對能量較低的候選結構進行局部優(yōu)化,這里我們使用的是最速下降法.由圖3可以看出,局部優(yōu)化后的結構的能量均有所下降(如圖3中的點劃線所示),但原子位置沒有很大的變動,對應的幻數(shù)也沒有變化.

3.2 Lennard-Jones三體勢函數(shù)

作為比較,我們接下來考慮的體系中,粒子之間的相互作用包括二體作用與三體作用,此時體系的總能量V為

這里θjik表示向量rij,rik之間的夾角,Vij和Vik分別表示為原子i與原子j和k之間的二體作用勢. 特別地,這里我們?nèi)(θ)=|θ? 2π/3|,所以θ=2π/3時f(θ)最小,從而使得基態(tài)體系容易形成角度為2π/3的鍵角.

我們利用類似的方法進行結構搜索,對應的幻數(shù)結構和能量變化曲線如圖5所示.該勢函數(shù)由兩部分貢獻,第一部分使得體系原子之間距離傾向于為1,第二部分使得體系的成鍵鍵角傾向于2π/3.在該勢函數(shù)下,原子團簇傾向于形成類似石墨烯的碎片結構.由圖5可見,這些較為穩(wěn)定的原子團簇結構是不斷向外生長出來的正六邊形.比如S10,S14,S16分別對應的是萘、蒽、芘的結構,S24對應于輪烯的結構.需要說明的是,基態(tài)結構和f(θ)的具體形式有關,當改變f(θ)的極值點,基態(tài)的結構會發(fā)生較大的變化.

圖5 (網(wǎng)刊彩色)Lennard-Jones三體勢函數(shù)下原子團簇的平均能量變化趨勢Fig.5. (color on line) The average energy of clusters as a function of the nuMber of atoMs under the Lennard-Jones th ree body interaction potential.

4 結 論

本文提出了一種結合結構識別和蒙特卡羅樹技術搜索穩(wěn)定結構的新方法.以平面團簇為例,我們考慮了包含二體和三體相互作用的Lennard-Jones勢函數(shù),并搜索確定了不同原子數(shù)的穩(wěn)定結構.我們選擇三角晶格團簇作為候選結構,引入編號策略對結構進行快速識別,并選擇能量較低的候選結構進行局域優(yōu)化來獲得對應體系的穩(wěn)定結構.我們的算法避免了相同結構的重復計算,建立了相鄰原子數(shù)的晶格團簇的互相聯(lián)系,采用蒙特卡羅樹技術兼顧搜索的深度和廣度.該方法可以對候選結構空間進行高效掃描,在較短時間內(nèi)更容易搜索到穩(wěn)定的團簇結構,并可以與第一原理計算方法結合,進行真實材料體系的結構搜索.

[1]Baletto F,Ferrando R 2005 Rev.Mod.Phys.77 371

[2]Gong X F,Wang Y,N ing X J 2008 Chin.Phys.Lett.25 468

[3]Liu T D,Zheng JW,Shao G F,Fan T E,Wen Y H 2015 Chin.Phys.B 24 33601

[4]Zhang M,Gao Y,Fang H P 2016 Chin.Phys.B 25 13602

[5]de Heer WA 1993 Rev.Mod.Phys.65 611

[6]K night WD,C leMenger K,Heer WA D,Saunders WA,Chou MY,Cohen ML 1984 Phys.Rev.Lett.52 2141

[7]Honeycutt J D,Andersen H C 1987 J.Phys.Chem.91 4950

[8]Liu G Q,Zhang Y L,Wang Z X,Wang Y Z,Zhang X X,Zhang WX 2012 CoMput.Theor.Chem.993 118

[9]Li S F,Zhao X J,Xu X S,Gao Y F,Zhang Z Y 2013 Phys.Rev.Lett.111 115501

[10]K iMS,Hwang SW,K iMMK,Shin D Y,Shin D H,K iMC O,Yang S B,Park J H,Hwang E,Choi SH,Ko G,SiMS,Sone C,ChoiH J,Bae S,Hong B H 2012 ACS Nano 6 8203

[11]Dah l J E,Liu SG,Carlson R MK 2003 Science 299 96

[12]Yang X B,Zhao Y J,Xu H,Yakobson B I 2011 Phys.Rev.B 83 205314

[13]Sergeeva A P,Popov I A,Piazza Z A,Li WL,RoManescu C,Wang L S,Boldy rev A I 2014 Acc.Chem.Res.47 1349

[14]Xu SG,Zhao Y J,Liao J H,Yang X B 2015 J.Chem.Phys.142 214307

[15]Hartke B 1993 J.Phys.Chem.97 9973

[16]Wang Y C,LüJ,Zhu L,Ma Y M2010 Phys.Rev.B 82 094116

[17]Frontera C,V ives E,Castan T,P lanes A 1995 Phys.Rev.B 51 11369

[18]K resse G,Jurgen H 1993 Phy.Rev.B 47 558

[19]Zhang Y J,X iao X Y,Li Y Q,Yan Y H 2012 Acta Phys.Sin.61 093602(in Chinese)[張英杰,肖緒洋,李永強,顏云輝2012物理學報61 093602]

[20]Liu T D,Li Z P,Ji Q S,Shao G F,Fan T E,Wen Y H 2017 Acta Phys.Sin.66 053601(in Chinese)[劉暾東,李澤鵬,季清爽,邵桂芳,范天娥,文玉華 2017物理學報 66 053601]

[21]Wu L J,Sui Q T,Zhang D,Zhang L,Qi Y 2015 Acta Phys.Sin.64 042102(in Chinese)[吳麗君,隨強濤,張多,張林,祁陽2015 物理學報64 042102]

[22]Li P F,Zhang Y G,Lei X L,Pan B C 2013 Acta Phys.Sin.62 143602(in Chinese)[李鵬飛,張艷革,雷雪玲,潘必才2013物理學報62 143602]

[23]LüJ,Wang Y C,Zhu L,Ma Y M2012 J.Chem.Phys.137 084104

[24]Oganov A R,G lass C W2006 J.Chem.Phys.124 244704

[25]Solovyov I A,Solovyov A V,G reiner W,Koshelev A,Shutovich A 2003 Phys.Rev.Lett.90 053401

[26]SWiechoWskiM,Mandziuk J,Ong Y S 2016 IEEE Trans.CoMp.Intel.AI.8 218

[27]V illar S S,BoWden J,Wason J 2015 Stat.Sci.30 199

[28]Sasaki Y,de Garis H 2004 Proceedings of the 2003 Congress on Evolutionary CoMputation Canberra,ACT,Australia,December 8–12,2003 p886

[29]Yang J,Zhang WQ 2007 Acta Phys.Sin.56 4017(in Chinese)[楊炯,張文清 2007物理學報 56 4017]

PACS:36.40.–c,31.15.–p,71.15.PdDOI:10.7498/aps.66.163601

*Project supported by the National Natural Science Foundation of China(G rant No.11474100)and the FundaMental Research Fund for the Central Universities,China(G rant No.2017MS119).

?Corresponding author.E-Mail:scxbyang@scut.edu.cn

Monte-Carlo tree search for stab le structu res of p lanar clusters?

He Chang-Chun Liao Ji-Hai Yang Xiao-Bao?

(DepartMent of Physics,South China University of Technology,Guangzhou 510640,China)

7 Ap ril 2017;revised Manuscrip t

18 June 2017)

Illustrated by the case of the p lanar clusters,we propose a neWMethod to search the possible stable structures by combining the structural identification and Monte-Carlo tree algorithm.We adopt two kinds of Model-potential to describe the interaction between atoMs:the pair interaction of Lennard-Jones potential and three-body interaction based on the Lennard-Jones potential.Taking the possible triangular lattice fragMent as candidates,we introduce a neWnomenclature to distinguish the structures,which can be used for the rapid congruence check.1)We label the atoMs on the triangu lar lattice according to the distances and the polar angles.where a given triangular structure has a corresponding serial number in the numbered p lane.Note that the congruent structures can have a group of possible serial numbers.2)We consider all the possible symmetrical operations including translation,inversion and rotation,and obtain the sMallest one for the unique noMenclature of the structure.In conventional search ofMagic clusters,the global optiMizations are perforMed for the structures With given number of atoMs.Herein,we perforMthe Monte-Carlo tree search to study the evolution of stab le structures With various numbers of atoMs.FroMthe structures of given number of atoMs,we saMp le the structures according to their energy With the iMportance saMp ling,and then expand the structures to the structures With oneMore atom,where the congruence check With the noMenclature is adop ted to avoid numerous repeated evaluations of candidates.Since the structures various numbers of atoMs are correlated With each other,a searching tree Will be obtained.In order to prevent the over-expansion of branches,we prove the“tree”according to energy to Make the tree asymMetric groWth to retain the loWenergy structure.The Width and dep th of search is balanced by the control of teMperature in the Monte-Carlo tree search.For the candidatesWith lower energies,we further perforMthe localop tiMization to obtain theMore stable structures.Our calculations shoWthat the triangular lattice fragmentsWill bemore stab le under the pair interaction of Lennard-Jones potential,which are in agreement With the previous studies.Under the three body interaction With the specifi c paraMeter,the hexagonal lattice fragMentsWill be More stable,which are siMilar to the configurations of graphene nano-flakes.Combining the congruence check and Monte-Carlo tree search,we provide an eff ective avenue to screen the possible candidatesand obtain the stab le structures in a shorter period of tiMe coMpared With the comMon global op tiMizationsWithout the structural identification,which can be extended to search the stable structure for Materials by the fi rst-princip les calculations.

cluster,structural identification,Monte-Carlo tree,global optiMization

10.7498/aps.66.163601

?國家自然科學基金(批準號：11474100)和中央高?；究蒲袠I(yè)務費(批準號：2017MS119)資助的課題．

?通信作者.E-Mail:scxbyang@scut.edu.cn

?2017中國物理學會C h inese P hysica l Society

http://Wu lixb.iphy.ac.cn