基於智能全局優化演算法的理論結構預測

作者：高朋越 呂健 王彥超 馬琰銘 (吉林大學物理學院 超硬材料國家重點實驗室)

基於群體智能的CALYPSO結構預測方法與軟體

1 理論結構預測難題

凝聚態物質內部的原子堆垛方式，即我們通常所說的微觀原子結構，是凝聚態物質最基本、最重要的信息，是揭開凝聚態物質的各種宏觀物理和化學性質產生根源的關鍵。一個典型例子是碳的同素異形體。如圖1 所示，石墨和金剛石都是由碳元素組成，但是由於原子排列方式的不同，石墨非常柔軟，具有良好的導電性，金剛石卻是目前已知最堅硬的材料，而且是一個寬頻隙的絕緣體。可見，即便是化學組分相同的材料，不同的微觀結構也可以導致截然不同的宏觀物理和化學性質。長期以來，確定物質的微觀原子結構一直是物理、化學、材料和生命等科學研究領域的焦點研究內容，先後有5 位科學家因在物質結構測量方法上的貢獻而獲得了諾貝爾獎(圖2)。

圖1 (a)石墨及其微觀晶體結構；(b)金剛石及其微觀晶體結構

目前確定微觀原子結構的實驗技術(如X射線衍射和中子散射等)已經相對完善，但實驗測量往往受到樣品純度、樣品所處環境和實驗信號的質量等因素的影響，特別是在某些極端條件(如高壓條件)下，僅通過實驗手段來確定結構仍面臨著諸多困難和挑戰。此時，發展只依據物質的化學組分來預測結構的理論方法尤為重要，不僅可以與實驗相輔相成最終確定物質的結構，還可以先於實驗開展材料結構的設計，根據材料的目標功能性質，逆向設計新型功能材料，指導實驗合成，節省實驗成本，縮短新材料的研發周期。

圖2 因在物質結構測量方法上的貢獻而獲得諾貝爾獎的科學家。Max van Laue 因發現晶體的X射線衍射而獲得1914 年諾貝爾物理學獎，William Henry Bragg和William Lawrence Bragg因用X射線分析晶體結構所做出的貢獻而獲得1915年諾貝爾物理學獎，HerbertA. Hauptman和Jerome Karle因在使用直接方法確定晶體結構方面的貢獻獲得了1985 年諾貝爾化學獎

理論上，在給定物質的化學組分和外界條件以後，其所有可能結構與對應的能量構成了一個高維度的勢能面(圖3)。根據能量最低原理，物質通常以其能量最低的狀態(即基態)存在，基態結構的能量處於勢能面上的全局最小值點。因此，理論結構預測的目標就是在勢能面上尋找到能量最小值點所對應的結構。物質的勢能面通常具有高度複雜性和多維度性，對於晶胞具有N個原子的晶體，其勢能面的維度是3N+3。這個高維的勢能曲面是由一系列局域的能量低谷組成的，每個能量低谷的極小值點對應著局域穩定結構，它們之間由勢壘分隔。具有相似結構特徵的能谷通常聚集在一起，形成漏斗形的大能谷，代表不同結構類型的大能谷之間通常存在著較高的勢壘。勢能面上能量的極小值點數量巨大，並隨著原子數的增加以指數方式增長，例如，具有13 個粒子的Lennard—Jones 體系，其勢能面上的能量極小值大約有103個。當粒子數增加到55 時，極小值至少有1012個。理論結構預測就是在如此龐大的結構群裡面找到全局能量最低的唯一結構，這是一個巨大的挑戰，在現有計算條件下，通過遍歷搜索勢能面上的所有結構來確定基態結構是不可能的。也正因為如此，1988 年John Maddox (《自然》雜誌社的主編)在《自然》期刊上發表社論說：「物理學的重要挑戰之一是只根據化學組分來確定物質的結構」。

圖3 二維勢能面示意圖

2 理論結構預測方法的發展

近些年來，隨著計算機計算能力的不斷提升和第一性原理計算方法的發展，只根據給定化學組分和外界條件(如壓力)，從理論上對材料結構進行預測已經成為可能。國內外科學家先後發展了多種可對物質勢能面進行智能全局搜索的結構預測方法，力爭通過較少勢能面取樣確定物質的結構，按照是否依賴初始給定結構，可把這些方法分為兩大類：一類是基於單一初始結構進行演化的躍遷勢壘方法；另一類是不依賴初始結構的群體搜索方法。

躍遷勢壘方法主要包括盆地跳躍、模擬退火和極小值跳躍方法。它們都是從一個人為給定的初始結構出發，通過結構在不同能谷之間的躍遷來尋找基態結構。盆地跳躍方法是蒙特卡羅技術與局域優化技術的結合。該方法通過蒙特卡羅模擬來進行結構的演化，並根據Metropolis 準則來決定是否接受新產生的結構。局域優化技術的引入可以把原始的勢能面投影到階梯型的勢能面上，從而提高了蒙特卡羅模擬的勢壘躍遷能力。模擬退火方法是源於對熱力學中退火過程的模擬，在蒙特卡羅或分子動力學的模擬過程中，通過緩慢降低溫度參數使體系最終達到基態結構。極小值跳躍方法同樣採用了高效的分子動力學模擬來進行結構演化。其主要思想是在分子動力學的模擬過程中，通過局域優化使體系達到能量極小值點，並記錄已經探索過的區域，通過不斷調整溫度參數，改變躍遷勢壘的能力，不斷探索未知區域。

群體搜索方法主要包括隨機尋找方法和基因遺傳演算法。隨機尋找方法完全隨機產生大量結構，然後對結構進行局域優化來尋找基態結構。由於其原理簡單並且實現方便，在中小體系的結構預測方面具有較為廣泛的應用。基因遺傳演算法是一種基於群體搜索的通用問題求解方法，它通過模擬達爾文的進化理論，利用「適者生存」和隨機信息交換的思想，通過交叉、變異和選擇等操作來尋找探索空間的最優解。1995年，Deaven 等人將遺傳演算法引入到團簇結構預測領域，從隨機結構出發，成功得到了C60的基態結構。隨後，該方法在結構預測領域得到了較多應用，國內外科學家先後對遺傳演算法的基本操作進行了不同程度的發展和改進，1999 年Woodley首次將遺傳演算法應用到了三維晶體結構預測領域，隨後遺傳演算法在晶體結構預測領域得到了更多推廣。這些方法在解決中小體系結構問題上獲得了巨大的成功，如預言了高壓下透明絕緣體鈉和離子硼高壓新相的結構等。

上述結構預測方法都是由國外科學家率先發展起來的，國內結構預測研究起步相對較晚，但發展極為迅速，目前已經在國際上佔據了重要一席。復旦大學劉智攀教授課題組發展了隨機勢能面搜索的結構預測方法。該方法通過給定初始結構在勢能面上的隨機行走，通過勢壘的跨越實現勢能面的探索，被成功地應用於SiO2、SrTiO3晶體和硼富勒烯等體系的結構研究。復旦大學龔新高教授研究組基於多目標差分進化演算法，發展了以材料的功能性質為導向的逆向設計方法，並編寫了相應的結構預測程序(IM2ODE)，該方法成功預言了新型窄帶隙TiO2半導體。2010年，我們課題組將基於群體智能理論的粒子群優化演算法與多種結構處理方法相結合，提出並發展了CALYPSO 結構預測方法，在此基礎上開發了擁有自主知識產權的CALYPSO 結構預測軟體包。它可以僅依據材料的化學組分和外界條件(如壓力)開展三維晶體、二維層狀材料、二維表面重構和零維團簇結構預測，並可以根據功能需求進行功能導向的材料結構設計(如超硬材料等)。目前，CALYPSO 已成為國際結構預測領域的主流軟體之一，在56 個國家和地區得到了推廣，被1700餘位同行簽訂版權協議來使用。

綜上所述，結構預測方法的發展在國內外已經取得了可喜的進展，不同方法對物質勢能面全局探索的策略也不盡相同、各具特色，它們在物質結構的研究中發揮了重要作用，解決了大量科學難題。下面將主要介紹我們課題組所發展的基於群體智能理論的CALYPSO 結構預測方法。

3 基於粒子群優化演算法的CALYPSO 結構預測方法

圖4 為CALYPSO 結構預測方法的流程，主要包括以下幾步：(1)在對稱性限制的條件下隨機產生一定數目的初始結構，作為群體的第一代；(2)對所產生結構進行相似性判斷，從而排除相似結構；(3)利用第一性原理或經驗勢方法對所產生的結構進行局域優化和能量計算；(4)構建下一代結構，其中群體中能量較低的部分結構(例如整個群體的60%)通過群智演算法(粒子群優化演算法)演化為新一代結構，能量較高的部分結構(例如整個群體的40%)被隨機產生的結構替換。CALYPSO 方法通過上述過程的不斷循環，實現整個群體結構的不斷演化和對勢能面的全局探索，主要包括以下4種關鍵的結構處理方法：

圖4 CALYPSO方法流程圖

(1)基於對稱性限制的隨機結構產生方法。為了保證在勢能面取樣的廣泛性和群體的多樣性，基於群體搜索的方法通常會完全隨機地產生初始結構。但隨著粒子數的增加，隨機產生的結構會越來越接近無序的液態結構，從而導致結構多樣性的迅速降低。我們研究發現，在隨機產生結構時，如果引入對稱性的限制，不但可以減小結構局域優化的變數，還可以有效地增加群體的多樣性，為結構預測提供一個較好的出發點。眾所周知，三維晶體結構的對稱性由230 種空間群來表徵，二維平面結構或者表面是由17 種平面空間群來描述，而零維的孤立分子和團簇具有點群對稱性。因此預測不同維度體系的結構，我們採用不同的對稱性規則限制。例如，對於三維晶體結構預測，我們在230 個空間群中隨機地選擇一個空間群，晶格參數按照這一空間群所對應的布拉菲晶格產生，原子位置通過這一空間群的Wyckoff的佔位組合而成。

(2)基於成鍵特徵矩陣的結構表徵方法。結構預測需要隨機生成大量的結構，並對它們進行局域結構優化和能量的計算。其中一定會存在著許多類似的結構，對應於勢能面上相同的能量低谷。它們經過局域結構優化以後會變成類似或同一個結構。對相似甚至相同結構的重複計算會嚴重影響結構預測的效率。如果我們能夠設計一個函數，通過該函數來唯一地表徵物質的結構，並能夠定量地表徵結構之間的相似程度，那麼就可以利用該函數來排除相似結構，提高結構預測的效率。這樣的函數還可以用來監控結構預測過程中群體多樣性的變化情況，並以此來引導結構的演化方向。為此，我們發展了一套名為成鍵特徵矩陣的結構表徵方法。該方法通過改進Steinhardt發展的鍵取向序參數技術來實現對鍵角信息的量化表徵，並且通過e 指數函數來量化鍵長信息。具體地說，當結構中兩個原子間距離小於給定的截斷長度時，記錄該成鍵的種類δAB，鍵長和鍵角信息。對於每個種類的成鍵，成鍵特徵矩陣表示為

其中rij表示由原子i 指向j 矢量的長度。θij，φij是rij在極坐標下的方位角，A(B)表示第i( j)個原子種類。Ylm是球諧函數， NδAB表示A與B元素成鍵的數目，bAB表示每種成鍵的最短長度，α是一個可調參數，為了消除坐標依賴，將這個矩陣用轉動不變的形式進行表示：

每一個結構可以通過這樣的矩陣來表徵。兩個結構的相似度通過它們對應的成鍵特徵矩陣之間的歐氏距離來判斷：

其中u和v分別表示兩個結構。

(3)結構的局域優化。在結構預測過程中，對產生的結構進行局域優化可以使其迅速達到對應能谷的極小值點，產生物理上更加合理的結構。這些極小值點對應的亞穩定結構將為產生下一代結構提供更合理的結構信息，顯著提高結構全局搜索的效率。因此，儘管局域優化增加了評估結構適應值(能量)所用的時間，但幾乎所有結構預測方法都採用了局域優化技術。目前，CALYPSO方法集成了多種基於第一性原理或經驗勢方法的軟體來實現結構的局域優化。

圖5 (a)一維勢能面上粒子群優化演算法示意圖；(b) 全局粒子群優化演算法示意圖；(c)局域粒子群優化演算法示意圖

(4)基於粒子群優化演算法的結構演化方法。粒子群優化演算法是由Eberhart 和Kennedy 於1995 年提出的一種基於群體搜索策略的全局優化演算法。它源於對鳥類捕食行為的模擬，把求解問題中的一個可能解看成是一個粒子，每個粒子在搜索空間中追尋最優粒子來進行探索，作為一種高效的多目標優化演算法，前期已經應用於系統辨識和神經網路訓練等領域。在2010 年，我們課題組首次將該演算法引入到晶體結構預測領域，發展了CALYPSO 結構預測方法。具體來說，從第二代開始，一定數目的新結構(默認是種群的60％)是通過粒子群優化演算法產生的，每一個結構被看成是搜索空間中的一個粒子，一系列結構組成一個群體。圖5(a)是一維勢能面上某一粒子運動的示意圖，在演化過程中粒子的位置通過下面的公式進行更新：

其中xit代表第i 個粒子在第t 代的位置，其在t+1代的速度vit+1可以通過粒子前一代的位置( xit)、當前最優的位置( pbestit)、整個種群中最好的粒子的位置( gbestit)以及粒子上一代的速度(vit)計算出來，具體的表達式如下：

其中c1是個體學習因子，表示粒子對過去自身經驗的依賴，而c2為全局學習因子，表示粒子對整個種群的信賴程度。早期研究表明c1和c2的值為2時，可以給出較好的優化結果。r1和r2是兩個分布在0 和1 之間的隨機數。這兩個隨機數可以保證我們的方法收斂到全局最小值點而不會陷入某個極小值點。顯然粒子在搜索空間中的移動受每個粒子過去經驗和種群經驗的影響，追隨能量最低的粒子向全局最小點移動。ω是慣性權重，當其取較大值的時候能夠提高演算法全局搜索的能力，而當取較小值的時候可以提高演算法搜索的精度。在CALYPSO 方法中，隨著迭代次數(iter)的增加，按照下面的公式，ω動態地從0.9 變到0.4：

其中ωmax為0.9， ωmin為0.4， itermax為最大演化代數。

目前CALYPSO 方法中採用了全局和局域兩個版本的粒子群優化演算法。如圖5(b)，(c)所示，在全局粒子群優化演算法中，所有粒子追尋目前最優的粒子運動，具有快速收斂的特點，適用於勢能面簡單的體系(原子數約小於30)；局域粒子群優化演算法將勢能面分成不同的區域，粒子追尋所在區域中的最優粒子進行探索。局域粒子群優化演算法可以看成是多個信息可以共享的全局粒子群優化演算法的組合，具有更強的抗過早熟能力，適用於勢能面複雜體系(模擬晶胞中的原子數大於30) 。

4 CALYPSO方法的應用

CALYPSO 方法的有效性已經通過對近百種已知實驗結構的測試得到了證實，並且在科研實踐中得到了進一步的檢驗，自2010 年軟體發布以來，7 年的時間裡，國內外同行利用CALYPSO軟體已經在Nature Chem.，Nature Commun.，PRL，PNAS，JACS 等國際頂級期刊發表了400 多篇論文，在物理、化學、材料和地球科學等領域解決了若干長期無法解決的科學難題，典型工作包括：

(1)破解了單質鋰的高壓半導體相的結構難題。作為高壓研究的模型體系，單質鋰在高壓下的相變極為複雜，其高壓下的結構問題一直是高壓研究的焦點。自從2002 年實驗發現單質鋰在60萬大氣壓以上存在新的高壓相以來，科學家們就開展了深入的理論和實驗研究。但由於理論技術和實驗條件的限制，新的高壓相的結構一直沒有得到解決。2009 年3 月，日本科學家在《自然》期刊上發表了單質鋰的高壓電學測量結果，令人意想不到的是，單質鋰的新的高壓相竟然是半導體，這更加激發了人們對其結構的研究熱情。隨後，國外幾個研究小組分別利用近幾年發展起來的基因演算法和隨機搜索演算法等晶體結構預測技術模擬了鋰的半導體相，但研究結果爭議很大。針對上述問題，我們利用CALYPSO 方法對鋰的高壓半導體相結構進行了系統的探索，預言了一個晶體學單胞內含有40 個原子的複雜底心正交Aba2-40 結構(Pearson 符號，oC40，圖6(a)) 。該結構的能量遠遠低於前期基因演算法和隨機搜索演算法所獲得的結構。oC40 結構中，鋰的價電子由於受到芯電子的排斥完全局域到晶格間隙之中，失去了自由電子的特性，金屬鋰變成了半導體。在我們的工作發表不久後，英國愛丁堡大學的Guillaume 等人獨立報道了鋰的高壓單晶X射線衍射實驗數據，實驗發現鋰的高壓半導體相的確具有oC40 結構。英國愛丁堡大學的Margues等人隨後通過實驗和理論相結合的方法進一步確定oC40 結構就是鋰的半導體相結構，支持了我們的理論預言。

圖6 (a)單質鋰的高壓半導體相oC40 結構；(b)地核環境下XeFe3的晶體結構；(c)高壓下H2S金屬相的晶體結構；(d)高壓下CaH6的晶體結構

(2)提出地核壓力與溫度條件下氙和鐵/鎳的化學反應。氙氣作為惰性氣體家族中的一員在大氣層中幾乎絕跡，與氬氣和氪氣等其他惰性氣體相比，90%以上的氙氣都不知所蹤，這在科學上被稱為「氙氣的消失之謎」。氙氣是否會儲存在佔據地球總質量三分之一的地核內部一直備受關注。地核的主要成分是鐵和鎳，如果氙氣儲存於地核中，就必須與鐵或者鎳在地核的壓強和溫度環境下(即360 GPa 和6000 K)發生化學反應，形成穩定的化合物。1997 年《科學》期刊發表的理論和實驗合作論文否定了氙氣和鐵發生反應的可能性。研究工作發表後的17 年間，科學界據此公認氙氣不可能儲存在地核內部。然而，《科學》文章的結論是以人為假定的鐵—氙化合物的結構為基礎，若假定錯誤，氙氣和鐵不反應的結論就不再成立。我們利用CALYPSO 方法對氙和鐵/鎳在地核環境下的結構進行了系統研究，提出了全新的鐵/鎳—氙化合物的結構形式(圖6(b))，構築了鐵/鎳—氙化合物的高溫—高壓相圖，首次給出了氙氣和鐵/鎳在地核環境下發生化學反應的證據，提出了氙氣被捕捉在地核內部的可能性。常壓條件下，惰性氣體氙難以和其他元素反應形成化合物；而鐵/鎳容易氧化失去電子，一般是帶正電的還原劑。令人意想不到的是，高壓下氙不僅和鐵/鎳發生了反應，而且電子從氙原子轉移到鐵/鎳原子上，從而竟使鐵/鎳成為了罕見的帶負電的氧化劑。這種高壓下非常規的電子轉移現象，誘導了氙氣與鐵/鎳之間的化學反應。

(3)H2S 和CaH6的高壓相結構與超導電性。超導領域的一個共識是BCS傳統超導體的超導溫度不可能超過40 K(McMillan 極限)。美國科學家Ashcroft 在2004 年提出，富含氫化合物一旦在高壓下金屬化就可能具有較高的超導溫度。然而富含氫化合物種類繁多，實驗科學家無法尋找到合適的富含氫化合物來開展超導實驗研究。此時，通過理論結構預測方法，快速、經濟地尋找最有希望的高溫超導體對實驗研究具有重要的指導意義。我們利用CALYPSO 方法，對H2S 在高壓下的結構進行了系統的探索，首次提出高壓金屬相H2S (圖6(c))是潛在的高溫超導體，在160 GPa壓力下，超導轉變溫度可以達到80 K。受到此項結果的啟發，德國馬克思—普朗克研究所Eremets 等人開展了H2S 的高壓超導實驗，證實了我們的理論預言，並創造了203 K的超導溫度新紀錄。此外，我們還利用CALYPSO 方法，以突破常規化學計量比的新型鹼土金屬氫化物為研究對象，預言了高含氫量的CaH6結構(圖6(d)) ，其中氫形成了極其獨特的二十面體籠型單元，籠型單元的H—H原子間形成弱的共價鍵，改變了含氫化合物中的只能形成氫分子和單原子氫的傳統認識，豐富了化學成鍵理論。基於BCS 理論，我們預言150 萬大氣壓下CaH6的理論超導轉變溫度達到220 K。

5 總結和展望

近年來，隨著智能全局優化演算法的引入，國內外科學家先後發展了多種理論結構預測方法，這些方法原理不盡相同、功能各具特色，已成為物質結構研究不可或缺的工具，被廣泛應用於結構現象豐富的研究領域，在物理、化學、材料、地球等科學領域，解決了大量長期無法解決的科學難題。

然而，理論結構預測方法的發展還遠沒有結束。不難發現，現有結構預測方法還只適用於原子數較少(百原子以內)的簡單體系，鮮有成功應用於大尺寸(百原子以上)複雜體系的算例。複雜大體系結構預測的困難主要體現在以下兩個方面：(1)勢能面的複雜程度隨著體系原子數目的增加而迅速增加，其維度線性增加，極小值點個數以指數形式增長，現有結構預測方法在處理複雜體系勢能面時均面臨取樣效率嚴重下降的問題；(2)計算資源的限制。現有較為準確的能量評估方法(如基於密度泛函理論的第一性原理方法)的計算量隨體系電子數目的增加呈3 次指數變化關係，複雜體系結構預測往往需要對大量候選結構進行能量評估，因此計算成本極其昂貴。目前發展針對複雜體系的結構預測方法還面臨著困難和挑戰，但隨著全局優化演算法的不斷發展、計算機計算能力的提升和基本物理化學理論的逐漸完善，實現複雜體系的結構預測指日可待。我們相信結構預測方法必將在更大、更複雜物質結構的研究中發揮重要作用。

本文選自《物理》2017年第9期

喜歡這篇文章嗎？立刻分享出去讓更多人知道吧！

本站內容充實豐富，博大精深，小編精選每日熱門資訊，隨時更新，點擊「搶先收到最新資訊」瀏覽吧！

請您繼續閱讀更多來自 中國物理學會期刊網 的精彩文章: