改進的模擬退火算法用于團簇結構優化

劉鳳萍

（上海海事大學信息工程學院，上海201306）

0 引言

團簇是由幾個乃至上千個原子、分子或離子通過物理或化學結合力組成的相對穩定的微觀或亞微觀聚集體，其物理和化學性質與其結構息息相關，不同的原子數目所能形成的穩定結構也不相同，因此研究不同原子數目的團簇的穩態結構是當今物理學和化學中的一個重要的前沿課題。應用儀器通過實驗實現團簇穩態結構的測定是十分困難而又不切實際的，目前為止，用的較為廣泛的一種辦法是：找到一種能夠表示團簇能量與其空間結構關系的能量函數，進而使用搜索算法搜索該能量函數的最小值點。該方法叫做從頭預測方法，其原理是基于熱物理理論：團簇的最穩定結構即為其能量最低時候的結構。

用于解決團簇結構優化問題的算法主要有遺傳算法[1-3]、模擬退火算法[4]和粒子群算法[5-6]等。隨原子數目增加，團簇的空間結構自由度也迅速增長，傳統搜索算法對團簇結構的優化效果也迅速減弱。對傳統的優化算法進行改進，是更好的選擇。因模擬退火算法具有計算過程簡單、魯棒性強等優點，故本文在模擬退火算法的基礎上進行算法改進。

1 能量模型

本文選擇蘭納瓊斯（Lennard-Jones）勢能公式來表示團簇能量與其空間結構之間的關系。

蘭納瓊斯勢是數學家John Lennard-Jones在1924年提出的，主要是用于模擬兩個電中性的分子或原子間相互作用勢能的數學模型，其表示見公式（1）。

上述公式中：ε表示勢能阱的深度，σ表示互相作用的勢能恰為零時的兩分子或原子間的距離，為了便于實驗結果的對比，本文的實驗中ε、σ兩個參數均參考文獻[7]，設置為1。rij表示兩體間的歐幾里得度量，其計算公式如下（2）：

含N個原子的團簇的蘭納瓊斯總勢能為：

2 搜索算法

2.1基本模擬退火算法

模擬退火算法思想受啟發于固體的退火過程：溫度很高的固體，其內部粒子都在做高速無序運動，具有較大的內能。當固體的溫度逐漸降低時，其內能也逐漸減小，內部粒子逐漸趨于有序狀態。直到固體溫度到達常溫時，其內能將至最低，這時內部粒子最為穩定。這恰似求解優化問題時，搜索目標函數最優解的過程。

模擬退火算法用于求解優化問題時，先從初始溫度出發，隨機生成初始解x，并計算其對應的目標函數值f(x)，溫度每降低一次產生一個新解，并計算其對應的目標函數值f)。當新解優于原解的時候接受新解，否則根據Metropolis準則判斷是否接受新解。Metropolis準則定義為：首先計算新解對應的目標函數值與原解目標函數值的差值?E，如下公式（4），若?E＞0，則接受新解；否則按概率P接受新解，概率P的計算見公式（5）。

算法流程如下：

Step1：設置初始溫度t，終止溫度tf，退火速度speed，隨機生成初始解x，計算目標函數f(x)；

Step2：若t＞tf，擾動產生新的解x，，計算對應目標函數f(x，)，否則終止算法；

Step3：根據Metropolis準則判斷新解是否被接受：如果接受則令變。t=t*speed，跳到Step2。

2.2基于種子策略的模擬退火算法的提出

（1）種子策略的提出

基本的SA算法具有較強的跳出局部最優解的能力，但是當搜索空間較大時，因其產生的解在多樣性上的局限性，很難求解出全局最優解。針對該問題，提出一次性產生多個初始解，分別對這些解進行適應度計算后，對產生的每個初始解進行擾動以產生一批新解，進而通過Metropolis準則判斷是否接受新解，循環直至到達終止溫度。

基本SA算法用于求解原子數目較少的團簇具有較好的效果，但是隨著原子數目增多，原子的空間結構自由度呈指數增長，基本SA算法的效果大大降低。考慮這樣一個現象：當某個含有N(其中N＞1)個原子的團簇W處于其能量最低狀態，也就是其穩定結構時，假設此時W中各原子的空間坐標為(X1,X2,…,Xi,…XN)，其中Xi表示第i個原子的三維坐標，即可表示為此時我們再往該團簇W中增加一個原子，該原子的坐標在區域范圍內隨機選取，設為產生一個新的團簇W，，此時新的團簇W，未必能達到其穩定狀態，但是一般情況下，這樣產生的團簇的能量會比隨機生成N+1個原子的三維坐標組合在一起的能量小，因此，本文采用了參考已求得其穩定狀態的原子數相對少的團簇的空間結構，初始化所求的原子數相對多的團簇的空間結構，進而再通過模擬退火算法進行結構優化的策略，此種策略稱為種子策略。

（2）算法流程

Step1：設置初始溫度t，終止溫度tf，退火速度speed，根據種子策略生成p個初始解x1,x2,…,xp，并計算各個解對應的目標函數

Step2：若t＞tf，擾動產生新的解，計算對應目標函數，否則終止算法；

Step3：根據Metropolis準則判斷新解是否被接受，t=t*speed，跳到 Step2。

3 實驗結果與分析

3.1參數設置

為檢驗提出的SS-SA算法對解決團簇結構優化問題的效果，選擇原子數在2到14之間的團簇的Len?nard-Jones勢能進行測試，并與基本模擬退火算法進行效果對比，進行試驗的兩個算法參數設置見下表1，其中t表示初始溫度，tf表示終止溫度，speed表示降溫速度。

表1 SA算法和SS-SA算法的參數設置

3.2實驗結果

SS-SA算法和SA算法對原子數為3到14的團簇結構優化的結果如表2，對于原子數相同的團簇，分別列出了SA算法和SS-SA算法對其Lennard-Jones勢能函數上獨立運行20次所得優化后勢能的最優值、平均值和標準差。表2中所展示的團簇對應的Lennard-Jones勢能的理想值參考于文獻[7]。

表2 SA算法和SS-SA算法對團簇結構的優化結果

SS-SA算法和SA算法分別對包含不同原子數的團簇單獨進行20次優化，可求得的團簇最低勢能與理想勢能之間存在一定誤差，誤差曲線見圖1；20次單獨優化過程中，SS-SA算法和SA算法對原子數目不同的團簇優化的結果的平均值與理想勢能也存在一定誤差，誤差曲線見圖2。

3.3結果分析

根據表2可知，當團簇的原子數目為3到5時，SA算法和SS-SA算法運行20次都可求出團簇的最穩定結構，而參考平均值和標準差，可知SS-SA算法穩定性強于SA算法；當原子數為6到10的時候，SA算法所能搜索到的對應團簇的最穩定結構已和理想值有所差距，但是SS-SA算法可求得團簇的最穩定結構；當原子數為11到14時，SA算法所求得的其最優結構已和其理想結構相去甚遠，而此時SS-SA所能求得的最優值與理想最優值誤差相對較小，認為此時可求得團簇的近似穩定結構，再參考均值和方差可知其較為穩定。綜上，對原子數在3到14之間的團簇進行結構優化時，本文提出的SS-SA算法，優化能力強于SA算法。

圖1 SA和SS-SA算法所求團簇能量最優值與理想值的誤差曲線

圖2 SA和SS-SA算法所求團簇能量平均值與理想值的誤差曲線

4 結語

為了解決團簇結構優化問題，本文所提出的在基本模擬退火算法的基礎上，增加了種子策略的算法——基于種子策略的模擬退火算法（SS-SA），提高了基本SA的算法性能。經過實驗驗證：用于解決團簇結構優化問題時，SS-SA算法能求得原子數目為3到10的團簇的穩定結構，對原子數目為11到14的團簇，能求得其近似最優結構。