劉鳳萍

（上海海事大學(xué)信息工程學(xué)院，上海201306）

0 引言

團(tuán)簇是由幾個(gè)乃至上千個(gè)原子、分子或離子通過(guò)物理或化學(xué)結(jié)合力組成的相對(duì)穩(wěn)定的微觀或亞微觀聚集體，其物理和化學(xué)性質(zhì)與其結(jié)構(gòu)息息相關(guān)，不同的原子數(shù)目所能形成的穩(wěn)定結(jié)構(gòu)也不相同，因此研究不同原子數(shù)目的團(tuán)簇的穩(wěn)態(tài)結(jié)構(gòu)是當(dāng)今物理學(xué)和化學(xué)中的一個(gè)重要的前沿課題。應(yīng)用儀器通過(guò)實(shí)驗(yàn)實(shí)現(xiàn)團(tuán)簇穩(wěn)態(tài)結(jié)構(gòu)的測(cè)定是十分困難而又不切實(shí)際的，目前為止，用的較為廣泛的一種辦法是：找到一種能夠表示團(tuán)簇能量與其空間結(jié)構(gòu)關(guān)系的能量函數(shù)，進(jìn)而使用搜索算法搜索該能量函數(shù)的最小值點(diǎn)。該方法叫做從頭預(yù)測(cè)方法，其原理是基于熱物理理論：團(tuán)簇的最穩(wěn)定結(jié)構(gòu)即為其能量最低時(shí)候的結(jié)構(gòu)。

用于解決團(tuán)簇結(jié)構(gòu)優(yōu)化問(wèn)題的算法主要有遺傳算法[1-3]、模擬退火算法[4]和粒子群算法[5-6]等。隨原子數(shù)目增加，團(tuán)簇的空間結(jié)構(gòu)自由度也迅速增長(zhǎng)，傳統(tǒng)搜索算法對(duì)團(tuán)簇結(jié)構(gòu)的優(yōu)化效果也迅速減弱。對(duì)傳統(tǒng)的優(yōu)化算法進(jìn)行改進(jìn)，是更好的選擇。因模擬退火算法具有計(jì)算過(guò)程簡(jiǎn)單、魯棒性強(qiáng)等優(yōu)點(diǎn)，故本文在模擬退火算法的基礎(chǔ)上進(jìn)行算法改進(jìn)。

1 能量模型

本文選擇蘭納瓊斯（Lennard-Jones）勢(shì)能公式來(lái)表示團(tuán)簇能量與其空間結(jié)構(gòu)之間的關(guān)系。

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

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

含N個(gè)原子的團(tuán)簇的蘭納瓊斯總勢(shì)能為：

2 搜索算法

2.1基本模擬退火算法

模擬退火算法思想受啟發(fā)于固體的退火過(guò)程：溫度很高的固體，其內(nèi)部粒子都在做高速無(wú)序運(yùn)動(dòng)，具有較大的內(nèi)能。當(dāng)固體的溫度逐漸降低時(shí)，其內(nèi)能也逐漸減小，內(nèi)部粒子逐漸趨于有序狀態(tài)。直到固體溫度到達(dá)常溫時(shí)，其內(nèi)能將至最低，這時(shí)內(nèi)部粒子最為穩(wěn)定。這恰似求解優(yōu)化問(wèn)題時(shí)，搜索目標(biāo)函數(shù)最優(yōu)解的過(guò)程。

模擬退火算法用于求解優(yōu)化問(wèn)題時(shí)，先從初始溫度出發(fā)，隨機(jī)生成初始解x，并計(jì)算其對(duì)應(yīng)的目標(biāo)函數(shù)值f(x)，溫度每降低一次產(chǎn)生一個(gè)新解，并計(jì)算其對(duì)應(yīng)的目標(biāo)函數(shù)值f)。當(dāng)新解優(yōu)于原解的時(shí)候接受新解，否則根據(jù)Metropolis準(zhǔn)則判斷是否接受新解。Metropolis準(zhǔn)則定義為：首先計(jì)算新解對(duì)應(yīng)的目標(biāo)函數(shù)值與原解目標(biāo)函數(shù)值的差值?E，如下公式（4），若?E＞0，則接受新解；否則按概率P接受新解，概率P的計(jì)算見(jiàn)公式（5）。

算法流程如下：

Step1：設(shè)置初始溫度t，終止溫度tf，退火速度speed，隨機(jī)生成初始解x，計(jì)算目標(biāo)函數(shù)f(x)；

Step2：若t＞tf，擾動(dòng)產(chǎn)生新的解x，，計(jì)算對(duì)應(yīng)目標(biāo)函數(shù)f(x，)，否則終止算法；

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

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

（1）種子策略的提出

基本的SA算法具有較強(qiáng)的跳出局部最優(yōu)解的能力，但是當(dāng)搜索空間較大時(shí)，因其產(chǎn)生的解在多樣性上的局限性，很難求解出全局最優(yōu)解。針對(duì)該問(wèn)題，提出一次性產(chǎn)生多個(gè)初始解，分別對(duì)這些解進(jìn)行適應(yīng)度計(jì)算后，對(duì)產(chǎn)生的每個(gè)初始解進(jìn)行擾動(dòng)以產(chǎn)生一批新解，進(jìn)而通過(guò)Metropolis準(zhǔn)則判斷是否接受新解，循環(huán)直至到達(dá)終止溫度。

基本SA算法用于求解原子數(shù)目較少的團(tuán)簇具有較好的效果，但是隨著原子數(shù)目增多，原子的空間結(jié)構(gòu)自由度呈指數(shù)增長(zhǎng)，基本SA算法的效果大大降低。考慮這樣一個(gè)現(xiàn)象：當(dāng)某個(gè)含有N(其中N＞1)個(gè)原子的團(tuán)簇W處于其能量最低狀態(tài)，也就是其穩(wěn)定結(jié)構(gòu)時(shí)，假設(shè)此時(shí)W中各原子的空間坐標(biāo)為(X1,X2,…,Xi,…XN)，其中Xi表示第i個(gè)原子的三維坐標(biāo)，即可表示為此時(shí)我們?cè)偻搱F(tuán)簇W中增加一個(gè)原子，該原子的坐標(biāo)在區(qū)域范圍內(nèi)隨機(jī)選取，設(shè)為產(chǎn)生一個(gè)新的團(tuán)簇W，，此時(shí)新的團(tuán)簇W，未必能達(dá)到其穩(wěn)定狀態(tài)，但是一般情況下，這樣產(chǎn)生的團(tuán)簇的能量會(huì)比隨機(jī)生成N+1個(gè)原子的三維坐標(biāo)組合在一起的能量小，因此，本文采用了參考已求得其穩(wěn)定狀態(tài)的原子數(shù)相對(duì)少的團(tuán)簇的空間結(jié)構(gòu)，初始化所求的原子數(shù)相對(duì)多的團(tuán)簇的空間結(jié)構(gòu)，進(jìn)而再通過(guò)模擬退火算法進(jìn)行結(jié)構(gòu)優(yōu)化的策略，此種策略稱為種子策略。

（2）算法流程

Step1：設(shè)置初始溫度t，終止溫度tf，退火速度speed，根據(jù)種子策略生成p個(gè)初始解x1,x2,…,xp，并計(jì)算各個(gè)解對(duì)應(yīng)的目標(biāo)函數(shù)

Step2：若t＞tf，擾動(dòng)產(chǎn)生新的解，計(jì)算對(duì)應(yīng)目標(biāo)函數(shù)，否則終止算法；

Step3：根據(jù)Metropolis準(zhǔn)則判斷新解是否被接受，t=t*speed，跳到 Step2。

3 實(shí)驗(yàn)結(jié)果與分析

3.1參數(shù)設(shè)置

為檢驗(yàn)提出的SS-SA算法對(duì)解決團(tuán)簇結(jié)構(gòu)優(yōu)化問(wèn)題的效果，選擇原子數(shù)在2到14之間的團(tuán)簇的Len?nard-Jones勢(shì)能進(jìn)行測(cè)試，并與基本模擬退火算法進(jìn)行效果對(duì)比，進(jìn)行試驗(yàn)的兩個(gè)算法參數(shù)設(shè)置見(jiàn)下表1，其中t表示初始溫度，tf表示終止溫度，speed表示降溫速度。

表1 SA算法和SS-SA算法的參數(shù)設(shè)置

3.2實(shí)驗(yàn)結(jié)果

SS-SA算法和SA算法對(duì)原子數(shù)為3到14的團(tuán)簇結(jié)構(gòu)優(yōu)化的結(jié)果如表2，對(duì)于原子數(shù)相同的團(tuán)簇，分別列出了SA算法和SS-SA算法對(duì)其Lennard-Jones勢(shì)能函數(shù)上獨(dú)立運(yùn)行20次所得優(yōu)化后勢(shì)能的最優(yōu)值、平均值和標(biāo)準(zhǔn)差。表2中所展示的團(tuán)簇對(duì)應(yīng)的Lennard-Jones勢(shì)能的理想值參考于文獻(xiàn)[7]。

表2 SA算法和SS-SA算法對(duì)團(tuán)簇結(jié)構(gòu)的優(yōu)化結(jié)果

SS-SA算法和SA算法分別對(duì)包含不同原子數(shù)的團(tuán)簇單獨(dú)進(jìn)行20次優(yōu)化，可求得的團(tuán)簇最低勢(shì)能與理想勢(shì)能之間存在一定誤差，誤差曲線見(jiàn)圖1；20次單獨(dú)優(yōu)化過(guò)程中，SS-SA算法和SA算法對(duì)原子數(shù)目不同的團(tuán)簇優(yōu)化的結(jié)果的平均值與理想勢(shì)能也存在一定誤差，誤差曲線見(jiàn)圖2。

3.3結(jié)果分析

根據(jù)表2可知，當(dāng)團(tuán)簇的原子數(shù)目為3到5時(shí)，SA算法和SS-SA算法運(yùn)行20次都可求出團(tuán)簇的最穩(wěn)定結(jié)構(gòu)，而參考平均值和標(biāo)準(zhǔn)差，可知SS-SA算法穩(wěn)定性強(qiáng)于SA算法；當(dāng)原子數(shù)為6到10的時(shí)候，SA算法所能搜索到的對(duì)應(yīng)團(tuán)簇的最穩(wěn)定結(jié)構(gòu)已和理想值有所差距，但是SS-SA算法可求得團(tuán)簇的最穩(wěn)定結(jié)構(gòu)；當(dāng)原子數(shù)為11到14時(shí)，SA算法所求得的其最優(yōu)結(jié)構(gòu)已和其理想結(jié)構(gòu)相去甚遠(yuǎn)，而此時(shí)SS-SA所能求得的最優(yōu)值與理想最優(yōu)值誤差相對(duì)較小，認(rèn)為此時(shí)可求得團(tuán)簇的近似穩(wěn)定結(jié)構(gòu)，再參考均值和方差可知其較為穩(wěn)定。綜上，對(duì)原子數(shù)在3到14之間的團(tuán)簇進(jìn)行結(jié)構(gòu)優(yōu)化時(shí)，本文提出的SS-SA算法，優(yōu)化能力強(qiáng)于SA算法。

圖1 SA和SS-SA算法所求團(tuán)簇能量最優(yōu)值與理想值的誤差曲線

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

4 結(jié)語(yǔ)

為了解決團(tuán)簇結(jié)構(gòu)優(yōu)化問(wèn)題，本文所提出的在基本模擬退火算法的基礎(chǔ)上，增加了種子策略的算法——基于種子策略的模擬退火算法（SS-SA），提高了基本SA的算法性能。經(jīng)過(guò)實(shí)驗(yàn)驗(yàn)證：用于解決團(tuán)簇結(jié)構(gòu)優(yōu)化問(wèn)題時(shí)，SS-SA算法能求得原子數(shù)目為3到10的團(tuán)簇的穩(wěn)定結(jié)構(gòu)，對(duì)原子數(shù)目為11到14的團(tuán)簇，能求得其近似最優(yōu)結(jié)構(gòu)。