改进的模拟退火算法用于团簇结构优化

刘凤萍

（上海海事大学信息工程学院，上海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的团簇，能求得其近似最优结构。