基于模拟退火优化算法的二元团簇结构优化方法

摘要

本发明公开了一种基于模拟退火优化算法的二元团簇结构优化方法。该方法通过模拟固体退火过程逐步降温以达到最低能量状态而获得二元团簇全局最稳定结构。基于退火进度表的扰动操作有效地拓展了算法的搜索空间，针对二元团簇中同时存在几何结构异构体以及由于不同类型的原子相对位置上的差异而引起的位置同分异构问题，算法中采用二元动态格点搜索方法加快了模拟退火中一个状态到另一个状态变化的速度，得到更为稳定的几何结构，还利用局部迭代搜索算法来解决位置同分异构问题。以二元Lennard‑Jones团簇势函数为优化算法的评价函数，通过模拟退火算法全局搜索与高效能的格点搜索及局部迭代搜索相结合来提高算法的效率，以获得最稳定的二元团簇结构。

说明书

技术领域

本发明涉及一种基于模拟退火优化算法的二元团簇结构优化方法。

背景技术

在催化、光学、电学等诸多领域，二元合金团簇在基础科学与应用方面吸引了广泛的研究热情。而这些性质与二元团簇的几何结构有着密切的联系，因此，确定其稳定几何构型是研究特殊性质的首要条件。但是，确定二元团簇最稳定结构是一项艰巨的任务。目前，人们已经设计了如遗传算法(GA)、basin-hopping算法(BH)、自适应免疫优化算法(AIOA)等算法并用于单一原子类型团簇的结构优化研究。而与单一原子类型团簇中仅存在几何构型同分异构现象相比，二元团簇中存在由于不同类型原子相对位置差异引起的位置同分异构问题。两者的同时存在使得寻找最优结构极为困难，易陷入不同构型的局部极优搜索空间，算法普遍存在效率低下的缺点。

发明内容

本发明要解决的技术问题是提供一种基于模拟退火优化算法的二元团簇结构优化方法。

为了解决上述技术问题，本发明采用的技术方案是，基于模拟退火优化算法的二元团簇结构优化方法，包括以下步骤：

步骤1：设置初始温度T_max，设计退火进度表：T_k＝T_max*exp[-1.5*(k/K)^0.25]；式中k为迭代次数，K为常数；

步骤2：产生初始结构库：确定二元A_nB_m团簇原子总数为n+m，其中n为A类型原子数目，m为B类型原子数目；随机地产生M个初始结构；根据二元Lennard-Jones团簇原子间相互作用计算结构库中各个个体的能量值，计算公式为：

其中，r_ij为个体中原子i和原子j间的距离；σ_ij代表原子i和原子j核间的距离；ε是势能井深度；并对这些初始结构采用限制内存的拟牛顿法进行局部优化操作，构成了初始结构库X₀；此时，k＝1；

步骤3：扰动操作：对结构库中的个体，在[1/3*r_max，2/3*r_max]范围内随机地选择I个原子，并对这些原子坐标x₀进行随机扰动，其中r_max为最大结构半径：x＝x₀±σ_AB*T_k，其中σ_AB为原子A和B核间距离，再进行局部优化操作，得到新解X；

步骤4：采用基于能量的原子移动操作、表面优化操作或二元动态格点搜索方法对新解X中所有的个体进行表面搜索，得到新解X′；

步骤5：位置异构体迭代操作：对步骤4得到的个体采用首次改进的局部搜索算法寻找最优的位置同分异构体；因此，基于解X′得到新解X“；

步骤6：判断能量是否降低：如果解X“的势能量值低于解X的值，则用新解X“替换解X，并返回步骤4；否则，当前能量最低的解为X₁，并转步骤7；

步骤7：根据Metropolis准则判断是否接受当前最优解X₁；若E(X₁)-E(X₀)≤0，则令X₀＝X₁，否则若exp(-(E(X₁)-E(X₀))/T_k)＞η，其中η为(0，1)上均匀分布的随机数，则令X₀＝X₁；此时，k＝k+1；

步骤8：如果满足k＜k_max，则转步骤9；否则，转到步骤3；

步骤9：输出全局最优值，即为二元团簇的最稳定结构。

作为优选，步骤1中，初始温度T_max和K值需要根据优化规模进行设定，T_max一般取值3.0，K取值为20。

作为优选，步骤3中，在[1/3*r_max，2/3*r_max]范围内随机地选择I个原子，I的一般取值范围为[36]，其中r_max为最大结构半径，并对这些原子坐标x₀进行与退火进度表有关的随机扰动，温度越高，扰动范围越大；随着温度降低，扰动范围变小，且此扰动操作提升了在团簇结构全局空间进行搜索的能力。

作为优选，步骤4中进行团簇结构表面操作时采用二元动态格点搜索方法。

本发明的有益效果是：

(1)基于退火进度表设计的团簇内部原子移动策略保障了模拟退火算法全局空间搜索能力；

(2)二元格点搜索方案在局部极小搜索空间中快速地找到最稳定几何构型；

(3)迭代的局部搜索方法降低了针对特定的几何构型寻找最优位置同分异构体的时间复杂度；(4)综合考虑了搜索广度、个体多样性、快速收敛性、全局搜索能力等方面的性能，以获取最稳定的二元团簇结构。

附图说明

下面结合附图和具体实施方式对本发明作进一步详细的说明。

图1是本发明实施例的流程图。

具体实施方式

本实施例是一种基于模拟退火优化算法用于二元团簇结构优化的方法，工作流程参阅图1所示。

以二元Lennard-Jones团簇A_nB_m(n+m＜30)的结构优化为例，具体实施步骤为：

(1)设置初始温度T_max＝3.0，设计退火进度表：T_k＝T_max*exp[-1.5*(k/K)^0.25]；式中k为迭代次数，K＝20。

(2)产生初始结构库。在一个半径为R＝σ_AB[(3*(n+m))/(4π)]^1/3(σ_AB为原子A和B核间距离)球形中随机地生成初始结构的坐标(n为A类型原子数目，m为B类型原子数目)，并从这n+m个坐标中随机地选择n个A类型原子，则剩余的为B类型原子。随机地产生M个初始结构。根据二元Lennard-Jones团簇原子间相互作用计算结构库中各个个体的能量值，计算公式为：

其中，r_ij为个体中原子i和原子j间的距离；σ_ij代表原子i和原子j核间的距离，σ_AA＝1.0，σ_BB＝1.3，σ_AB＝1.15；ε是势能井深度，取值为1.0。并对这些初始结构采用限制内存的拟牛顿法进行局部优化操作，构成了初始结构库X₀；此时，k＝1。

(3)扰动操作。对结构库中的个体，在[1/3*r_max，2/3*r_max](r_max为最大结构半径)范围内随机地选择5个原子，并对这些原子坐标x₀进行随机扰动：x＝x₀±σ_AB*T_k(σ_AB为原子A和B核间距离)，再进行局部优化操作，得到新解X。

(4)二元动态格点搜索方法(binary dynamic lattice searching，BDLS)。对新解X中所有的个体采用二元动态格点搜索方法进行表面搜索，得到新解X′。

(5)位置异构体迭代操作。对(4)得到的个体采用首次改进的局部搜索算法寻找最优的位置同分异构体。因此，基于解X′得到新解X“。

(6)判断能量是否降低。如果解X“的势能量值低于解X的值，则用新解X“替换解X，并返回(4)；否则，当前能量最低的解为X₁，并转(7)。

(7)根据Metropolis准则判断是否接受当前最优解X₁。若E(X₁)-E(X₀)≤0，则令X₀＝X₁，否则若exp(-(E(X₁)-E(X₀))/T_k)＞η，其中η为(0，1)上均匀分布的随机数，则令X₀＝X₁。此时，k＝k+1。

(8)如果满足k＜k_max(k_max＝10)，则转步骤9；否则，转到步骤3。

(9)输出全局最优值，即为二元团簇的最稳定结构。

作为优选，步骤1中，初始温度T_max和K值需要根据优化规模进行设定，T_max一般取值3.0，K取值为20。

步骤4中进行团簇结构表面操作时还可以采用基于能量的原子移动操作或表面优化操作进行表面搜索。

各个尺寸二元团簇结构的优化结果见表1。

表1

表1显示了本实施例与传统遗传算法(GA)、basin-hopping(BH)算法在优化部分二元Lennard-Jones团簇稳定结构的优化结果。GA与BH算法运行时的初始库数量均为20，最大迭代数目为100。各种算法独立运行100次，统计最优结果。如表1所示，本方法能够优化到所有最稳定结构，并且对于除尺寸为19的二元团簇之外，其余团簇结构优化的成功率均为100％，且明显地优于其他算法，表现出优异的优化性能。此外，从优化时间上看，该方法寻找到最稳定结构所需时间远少于其他算法。由此可见，该算法具有优异的全局搜索能力以及快速收敛速度。

以上所述的本发明实施方式，并不构成对本发明保护范围的限定。任何在本发明的精神和原则之内所作的修改、等同替换和改进等，均应包含在本发明的权利要求保护范围之内。