梁雪梅,黃吳夕,馮子儼,陳恒杰
(重慶科技學院 數理與大數據學院,重慶 401331)
團簇是物質從微觀到宏觀的過渡狀態(tài),是物質存在的重要形式之一。分子團簇的結構對其性質有著重要的影響,因此,對分子團簇研究的首要任務為確定其結構。然而,現有的實驗很難分離得到不同尺度的團簇,對團簇結構的測定則更加困難,特別是當團簇尺度較大時,這項任務變的幾乎不可能。因此,在團簇結構的研究方面往往需要結合理論計算幫助獲得團簇的結構[1,2]。目前針對團簇結構主要有兩個研究方向:一、利用基于量子力學的各類從頭計算軟件得到團簇在某一結構下的單點能,利用有偏或無偏優(yōu)化算法來更新結構,通過多次迭代得到最低能結構,限于當前計算機資源的限制,這種方式也僅在中小尺度的團簇下適用,對中大尺度團簇,由于耗費的計算資源巨量而很難開展。二、利用實驗或基于準確的量子理論計算下的團簇特性,構造出一個合理的勢能函數,然后用該勢能函數結合各類優(yōu)化算法對包括中大尺度在內的分子團簇結構進行計算或預測[3-5],其數學本質為一能量最小化的最優(yōu)化問題。兩種方法都需要研究者相當熟悉相應的優(yōu)化算法,強大的編程能力,因此學習曲線相對較長[6-8]。
1stOpt被認為是當前全世界最好的全局優(yōu)化程序[9],其強大的通用/穩(wěn)健全局優(yōu)化算法(UGO),接近數學描述且簡單的編程語言和靈活多變的編程方式。我們想能不能通過二次開發(fā),利用1stOpt的優(yōu)化功能,在團簇優(yōu)化時基于該軟件進行結構更新,這樣研究者可節(jié)約大量研究算法本身及其算法編寫的時間和精力,使其能更多地關注到團簇問題本身。同時也可以通過該問題地解決來測試1stOpt的優(yōu)化能力,給我們解決其他問題奠定基礎。
本文共測試了三種勢函數,rij表示各原子間距,本質是對角線為零的對稱方陣,為自變量。ULJ,UG和USC分別表示函數的勢能,為因變量,N為構成團簇的原子總數,三種函數的具體形式如下[10,11]:
Lennard-Jones勢:
(1)
Gupta勢:
(2)
(3)
Sutton-Chen勢:
(4)
(5)
圖1為程序設計思想及流程圖,其中結構更新由1stOpt黑箱實現,是團簇優(yōu)化中尋找最低能結構的核心問題,也是1stOpt的優(yōu)勢所在。收斂判斷部分由1stOpt的圖形界面設定,其余由1stOpt結合內嵌語言或動態(tài)庫編程實現。原則上自變量可選擇笛卡爾坐標(xyz),分子間距(rij)或分子內坐標(鍵長,鍵角,二面角等),為了計算的便利,本文選擇xyz型坐標,可通過簡單處理或借助于量化工具即可得到分子間距或分子內坐標。由于分子是可平移,旋轉的,所以優(yōu)化的xyz坐標不唯一,但與之對應的分子間距、內坐標和最低能是唯一的。
圖1 程序流程圖
結合1stOpt編程,分別用內嵌程序語言(如Pascal,Fortran等)以及調用外部編譯器編寫和編譯的動態(tài)庫(dll),實現了上述三種勢函數的單原子分子團簇結構優(yōu)化程序,探索了雙原子合金團簇結構優(yōu)化程序,計算中采用全局優(yōu)化算法或穩(wěn)健全局優(yōu)化算法,自動選擇,根據團簇尺度大小,依次選擇并行數100~1 000,比率0.3~0.5,控制數30~50,收斂判斷迭代50,變異率0.3~0.50。為了檢驗編寫程序的正確性,本文所有工作同時用ABCluster軟件進行了計算。
2.1.1 Gupta勢
計算時,取參數A=0.006 946 731,ξ=0.437 81,d=3.5,p=24.943 17,q=0.292 646。表1為基于Gupta勢函數,結合內嵌Pascal語言編寫的1stOpt程序,用戶只需改變其中的原子個數num,執(zhí)行程序后便可獲得相應尺度團簇中每個原子的xyz坐標和最低能量?;谠摮绦?,我們測試了num=2~20時的團簇,結果見表2,可以看出,本文完美地再現了Darby利用遺傳算法[12]和王全武等利用人工蜂群算法[13]得到的結果,說明利用1stOpt編程實現團簇結構優(yōu)化是可行的,同時也驗證了當前程序是可靠的。編寫該程序時,研究者并不需要任何優(yōu)化算法理論知識,也不需要將優(yōu)化算法轉化為程序的經驗,只需簡單了解團簇問題涉及的勢函數,簡單學習1stOpt編程即可完成以往需要研究者花費大量時間,前期研究算法,后期大量編程才能完成的工作,將大大降低該研究問題的知識儲備和門檻。
表1 Gupta勢函數計算程序
表2 Gupta勢函數結果
2.1.2 Sutton-Chen勢
為進一步說明程序的可靠性,我們也針對文獻中經常使用的Sutton-Chen勢函數進行編寫,計算時取參數ε=1,c=39.432,a=1.0,n=9,m=6,程序見表3,計算結果見表4。
表3 Sutton-Chen勢函數計算程序
表4 Sutton-Chen勢函數結果
為了比較各團簇優(yōu)化算法的有效性,劍橋大學的Doye等搜集了多種勢函數下已出版的最低能結構和最低能量值,本文后續(xù)均標記為Best。Doye沒有給出n=2時的能量,我們優(yōu)化的結果為-269.566 745 129 083,另外,同Gupta函數結果一樣,在不考慮精度的前提下,當前獲得的最優(yōu)結果與Doye給出的最低能完全一致,再次檢驗了1stOpt優(yōu)化原子團簇的可行性與可靠性。
2.1.3 LJ勢
在一些情況下,若研究者手上已有計算勢函數的子程序,或熟悉的其他計算機語言非內嵌在1stOpt中,且研究者也并不想過多的學習其他內嵌語言。這時,應用熟悉的計算機程序語言編寫子程序并編譯成動態(tài)庫,再利用1stOpt調用該動態(tài)庫則顯得非常合適,這樣能大大縮短程序開發(fā)周期,節(jié)約研究時間。表5展示了用Compaq Visual Fortran編寫的基于LJ勢函數的動態(tài)庫源文件,接著,在Windows操作系統下,打開cmd窗口,并cd到存放該源文件的當前目錄,用命令DF.exe/all“l(fā)j.f90”進行編譯,之后將生成lj.dll文件及其他幾個支撐文件,再用表6所示的1stOpt程序調用lj.dll,由于所有文件都在當前路徑的同一個文件夾下,1stOpt程序調用時沒有采用絕對路徑。
為便于比較,LJ中的參數采用ε=1,σ=1,計算結果見表7??梢?,利用1stOpt調用動態(tài)庫結果和ABCluster計算的結果完全匹配,也與Doye給出的最優(yōu)結果一致,不僅實踐了動態(tài)庫源程序的編寫,編譯和調用,也進一步檢驗了1stOpt優(yōu)化團簇結構的可行性。
表5 LJ勢函數下的動態(tài)庫源程序
表6 1stOpt調用動態(tài)庫
表7 Lennard-Jones勢函數結果
續(xù)表
原子團簇是最簡單的團簇形式,而很多情況下,團簇可能由多種原子組合而成,或者大量的幾種原子加少量某種原子摻雜而成,如AgXCuy,AgXCuyAuz等。接著我們以二元合金AgXCuy團簇為例,探索了1stOpt在解決復雜團簇結構方面的可行性,程序見表8,計算結果見表9。毫無意外,1stOpt完全重復了ABCluster結果,兩者完美的自洽,基于我們對ABCluster計算的經驗和單原子分子團簇計算結果的比較,有理由相信兩種計算結果得到的都是最低能。
表8 二元合金團簇計算程序
續(xù)表
表9 二元合金團簇AgXCuy能量優(yōu)化結果
基于三種常見的勢能函數,結合1stOpt及其內嵌的程序語言和外部調用動態(tài)庫,多種程序,多種手段,初步探索了1stOpt在優(yōu)化團簇結構方面的可行性,與Doye搜集的最優(yōu)結果和專業(yè)團簇優(yōu)化程序ABCluster計算結果比較,我們發(fā)現基于1stOpt的程序完全能夠重現上述結果,因此利用1stOpt優(yōu)化團簇結構是可行的。同時在執(zhí)行過程中也發(fā)現了一些不足。在原子數目較少時,基于通用/穩(wěn)健全局優(yōu)化算法的1stOpt效率非常高,很多情況下一次或者幾次迭代就能找到最低能,然而,當n增大即變量增多時,盡管最終也能找到最低能,但效率會大打折扣,有時還會陷入局部最優(yōu),我們猜測可能與以下幾方面有關:一、我們的探索目前僅在停留在利用1stOpt實現團簇優(yōu)化功能的可行性,因此編寫的程序未進行任何優(yōu)化和通用化;二、1stOpt目前僅支持Windows操作系統,而ABCluster計算是在Linux服務器上完成,兩種操作系統在效率上天然的優(yōu)劣,是造成效率差異的最最重要原因之一;三、ABCluster專門為團簇優(yōu)化設計,為有偏算法,具有加速功能,而1stOpt為通用的優(yōu)化程序,具有無偏性,因此在處理該類特定問題時效率沒能大幅提升?;谏鲜鰡栴},提出幾點展望:首先可通過文本等操作、對自編的程序進行優(yōu)化,解決同一問題盡量減少變量的同時,使得程序通用化(如多元合金不需要再寫更復雜的程序,僅給出原子種類和相互間參數即可);其次,通過動態(tài)庫/內嵌語言,參數設置等實現有偏操作,縮小搜索區(qū)間,減少優(yōu)化時間,提高優(yōu)化效率;最后,也是研究者不可控的一點,期望七維高科開發(fā)出基于Linux操作系統的版本,開發(fā)出更加靈活的動態(tài)庫接口。