基于智能全局优化算法的理论结构预测

基于智能全局优化算法的理论结构预测 | 视频

2017 年 10 月 11 日 中国物理学会期刊网 高朋越等

作者：高朋越吕健王彦超马琰铭 (吉林大学物理学院超硬材料国家重点实验室)

基于群体智能的CALYPSO结构预测方法与软件

1 理论结构预测难题

凝聚态物质内部的原子堆垛方式，即我们通常所说的微观原子结构，是凝聚态物质最基本、最重要的信息，是揭开凝聚态物质的各种宏观物理和化学性质产生根源的关键。一个典型例子是碳的同素异形体。如图1 所示，石墨和金刚石都是由碳元素组成，但是由于原子排列方式的不同，石墨非常柔软，具有良好的导电性，金刚石却是目前已知最坚硬的材料，而且是一个宽带隙的绝缘体。可见，即便是化学组分相同的材料，不同的微观结构也可以导致截然不同的宏观物理和化学性质。长期以来，确定物质的微观原子结构一直是物理、化学、材料和生命等科学研究领域的焦点研究内容，先后有5 位科学家因在物质结构测量方法上的贡献而获得了诺贝尔奖(图2)。

图1 (a)石墨及其微观晶体结构；(b)金刚石及其微观晶体结构

目前确定微观原子结构的实验技术(如X射线衍射和中子散射等)已经相对完善，但实验测量往往受到样品纯度、样品所处环境和实验信号的质量等因素的影响，特别是在某些极端条件(如高压条件)下，仅通过实验手段来确定结构仍面临着诸多困难和挑战。此时，发展只依据物质的化学组分来预测结构的理论方法尤为重要，不仅可以与实验相辅相成最终确定物质的结构，还可以先于实验开展材料结构的设计，根据材料的目标功能性质，逆向设计新型功能材料，指导实验合成，节省实验成本，缩短新材料的研发周期。

图2 因在物质结构测量方法上的贡献而获得诺贝尔奖的科学家。Max van Laue 因发现晶体的X射线衍射而获得1914 年诺贝尔物理学奖，William Henry Bragg和William Lawrence Bragg因用X射线分析晶体结构所做出的贡献而获得1915年诺贝尔物理学奖，Herbert A. Hauptman和Jerome Karle因在使用直接方法确定晶体结构方面的贡献获得了1985 年诺贝尔化学奖

理论上，在给定物质的化学组分和外界条件以后，其所有可能结构与对应的能量构成了一个高维度的势能面(图3)。根据能量最低原理，物质通常以其能量最低的状态(即基态)存在，基态结构的能量处于势能面上的全局最小值点。因此，理论结构预测的目标就是在势能面上寻找到能量最小值点所对应的结构。物质的势能面通常具有高度复杂性和多维度性，对于晶胞具有N个原子的晶体，其势能面的维度是3N+3。这个高维的势能曲面是由一系列局域的能量低谷组成的，每个能量低谷的极小值点对应着局域稳定结构，它们之间由势垒分隔。具有相似结构特征的能谷通常聚集在一起，形成漏斗形的大能谷，代表不同结构类型的大能谷之间通常存在着较高的势垒。势能面上能量的极小值点数量巨大，并随着原子数的增加以指数方式增长，例如，具有13 个粒子的Lennard—Jones 体系，其势能面上的能量极小值大约有10³个。当粒子数增加到55 时，极小值至少有10¹²个。理论结构预测就是在如此庞大的结构群里面找到全局能量最低的唯一结构，这是一个巨大的挑战，在现有计算条件下，通过遍历搜索势能面上的所有结构来确定基态结构是不可能的。也正因为如此，1988 年John Maddox (《自然》杂志社的主编)在《自然》期刊上发表社论说：“物理学的重要挑战之一是只根据化学组分来确定物质的结构”。

图3 二维势能面示意图

2 理论结构预测方法的发展

近些年来，随着计算机计算能力的不断提升和第一性原理计算方法的发展，只根据给定化学组分和外界条件(如压力)，从理论上对材料结构进行预测已经成为可能。国内外科学家先后发展了多种可对物质势能面进行智能全局搜索的结构预测方法，力争通过较少势能面取样确定物质的结构，按照是否依赖初始给定结构，可把这些方法分为两大类：一类是基于单一初始结构进行演化的跃迁势垒方法；另一类是不依赖初始结构的群体搜索方法。

跃迁势垒方法主要包括盆地跳跃、模拟退火和极小值跳跃方法。它们都是从一个人为给定的初始结构出发，通过结构在不同能谷之间的跃迁来寻找基态结构。盆地跳跃方法是蒙特卡罗技术与局域优化技术的结合。该方法通过蒙特卡罗模拟来进行结构的演化，并根据Metropolis 准则来决定是否接受新产生的结构。局域优化技术的引入可以把原始的势能面投影到阶梯型的势能面上，从而提高了蒙特卡罗模拟的势垒跃迁能力。模拟退火方法是源于对热力学中退火过程的模拟，在蒙特卡罗或分子动力学的模拟过程中，通过缓慢降低温度参数使体系最终达到基态结构。极小值跳跃方法同样采用了高效的分子动力学模拟来进行结构演化。其主要思想是在分子动力学的模拟过程中，通过局域优化使体系达到能量极小值点，并记录已经探索过的区域，通过不断调整温度参数，改变跃迁势垒的能力，不断探索未知区域。

群体搜索方法主要包括随机寻找方法和基因遗传算法。随机寻找方法完全随机产生大量结构，然后对结构进行局域优化来寻找基态结构。由于其原理简单并且实现方便，在中小体系的结构预测方面具有较为广泛的应用。基因遗传算法是一种基于群体搜索的通用问题求解方法，它通过模拟达尔文的进化理论，利用“适者生存”和随机信息交换的思想，通过交叉、变异和选择等操作来寻找探索空间的最优解。1995年，Deaven 等人将遗传算法引入到团簇结构预测领域，从随机结构出发，成功得到了C₆₀的基态结构。随后，该方法在结构预测领域得到了较多应用，国内外科学家先后对遗传算法的基本操作进行了不同程度的发展和改进，1999 年Woodley首次将遗传算法应用到了三维晶体结构预测领域，随后遗传算法在晶体结构预测领域得到了更多推广。这些方法在解决中小体系结构问题上获得了巨大的成功，如预言了高压下透明绝缘体钠和离子硼高压新相的结构等。

上述结构预测方法都是由国外科学家率先发展起来的，国内结构预测研究起步相对较晚，但发展极为迅速，目前已经在国际上占据了重要一席。复旦大学刘智攀教授课题组发展了随机势能面搜索的结构预测方法。该方法通过给定初始结构在势能面上的随机行走，通过势垒的跨越实现势能面的探索，被成功地应用于SiO₂、SrTiO₃晶体和硼富勒烯等体系的结构研究。复旦大学龚新高教授研究组基于多目标差分进化算法，发展了以材料的功能性质为导向的逆向设计方法，并编写了相应的结构预测程序(IM²ODE)，该方法成功预言了新型窄带隙TiO₂半导体。2010年，我们课题组将基于群体智能理论的粒子群优化算法与多种结构处理方法相结合，提出并发展了CALYPSO 结构预测方法，在此基础上开发了拥有自主知识产权的CALYPSO 结构预测软件包。它可以仅依据材料的化学组分和外界条件(如压力)开展三维晶体、二维层状材料、二维表面重构和零维团簇结构预测，并可以根据功能需求进行功能导向的材料结构设计(如超硬材料等)。目前，CALYPSO 已成为国际结构预测领域的主流软件之一，在56 个国家和地区得到了推广，被1700余位同行签订版权协议来使用。

综上所述，结构预测方法的发展在国内外已经取得了可喜的进展，不同方法对物质势能面全局探索的策略也不尽相同、各具特色，它们在物质结构的研究中发挥了重要作用，解决了大量科学难题。下面将主要介绍我们课题组所发展的基于群体智能理论的CALYPSO 结构预测方法。

3 基于粒子群优化算法的CALYPSO 结构预测方法

图4 为CALYPSO 结构预测方法的流程，主要包括以下几步：(1)在对称性限制的条件下随机产生一定数目的初始结构，作为群体的第一代；(2)对所产生结构进行相似性判断，从而排除相似结构；(3)利用第一性原理或经验势方法对所产生的结构进行局域优化和能量计算；(4)构建下一代结构，其中群体中能量较低的部分结构(例如整个群体的60%)通过群智算法(粒子群优化算法)演化为新一代结构，能量较高的部分结构(例如整个群体的40%)被随机产生的结构替换。CALYPSO 方法通过上述过程的不断循环，实现整个群体结构的不断演化和对势能面的全局探索，主要包括以下4种关键的结构处理方法：

图4 CALYPSO方法流程图

(1)基于对称性限制的随机结构产生方法。为了保证在势能面取样的广泛性和群体的多样性，基于群体搜索的方法通常会完全随机地产生初始结构。但随着粒子数的增加，随机产生的结构会越来越接近无序的液态结构，从而导致结构多样性的迅速降低。我们研究发现，在随机产生结构时，如果引入对称性的限制，不但可以减小结构局域优化的变量，还可以有效地增加群体的多样性，为结构预测提供一个较好的出发点。众所周知，三维晶体结构的对称性由230 种空间群来表征，二维平面结构或者表面是由17 种平面空间群来描述，而零维的孤立分子和团簇具有点群对称性。因此预测不同维度体系的结构，我们采用不同的对称性规则限制。例如，对于三维晶体结构预测，我们在230 个空间群中随机地选择一个空间群，晶格参数按照这一空间群所对应的布拉菲晶格产生，原子位置通过这一空间群的Wyckoff的占位组合而成。

(2)基于成键特征矩阵的结构表征方法。结构预测需要随机生成大量的结构，并对它们进行局域结构优化和能量的计算。其中一定会存在着许多类似的结构，对应于势能面上相同的能量低谷。它们经过局域结构优化以后会变成类似或同一个结构。对相似甚至相同结构的重复计算会严重影响结构预测的效率。如果我们能够设计一个函数，通过该函数来唯一地表征物质的结构，并能够定量地表征结构之间的相似程度，那么就可以利用该函数来排除相似结构，提高结构预测的效率。这样的函数还可以用来监控结构预测过程中群体多样性的变化情况，并以此来引导结构的演化方向。为此，我们发展了一套名为成键特征矩阵的结构表征方法。该方法通过改进Steinhardt发展的键取向序参数技术来实现对键角信息的量化表征，并且通过e 指数函数来量化键长信息。具体地说，当结构中两个原子间距离小于给定的截断长度时，记录该成键的种类δ_AB，键长和键角信息。对于每个种类的成键，成键特征矩阵表示为

其中r_ij表示由原子i 指向j 矢量的长度。θ_ij，φ_ij是r_ij在极坐标下的方位角，A(B)表示第i( j)个原子种类。Y_lm 是球谐函数， N_δAB表示A与B元素成键的数目，b_AB表示每种成键的最短长度，α是一个可调参数，为了消除坐标依赖，将这个矩阵用转动不变的形式进行表示：

每一个结构可以通过这样的矩阵来表征。两个结构的相似度通过它们对应的成键特征矩阵之间的欧氏距离来判断：

其中u和v分别表示两个结构。

(3)结构的局域优化。在结构预测过程中，对产生的结构进行局域优化可以使其迅速达到对应能谷的极小值点，产生物理上更加合理的结构。这些极小值点对应的亚稳定结构将为产生下一代结构提供更合理的结构信息，显著提高结构全局搜索的效率。因此，尽管局域优化增加了评估结构适应值(能量)所用的时间，但几乎所有结构预测方法都采用了局域优化技术。目前，CALYPSO方法集成了多种基于第一性原理或经验势方法的软件来实现结构的局域优化。

图5 (a)一维势能面上粒子群优化算法示意图；(b) 全局粒子群优化算法示意图；(c)局域粒子群优化算法示意图

(4)基于粒子群优化算法的结构演化方法。粒子群优化算法是由Eberhart 和Kennedy 于1995 年提出的一种基于群体搜索策略的全局优化算法。它源于对鸟类捕食行为的模拟，把求解问题中的一个可能解看成是一个粒子，每个粒子在搜索空间中追寻最优粒子来进行探索，作为一种高效的多目标优化算法，前期已经应用于系统辨识和神经网络训练等领域。在2010 年，我们课题组首次将该算法引入到晶体结构预测领域，发展了CALYPSO 结构预测方法。具体来说，从第二代开始，一定数目的新结构(默认是种群的60％)是通过粒子群优化算法产生的，每一个结构被看成是搜索空间中的一个粒子，一系列结构组成一个群体。图5(a)是一维势能面上某一粒子运动的示意图，在演化过程中粒子的位置通过下面的公式进行更新：

其中x_i^t 代表第i 个粒子在第t 代的位置，其在t+1代的速度v_it+1可以通过粒子前一代的位置( x_i^t )、当前最优的位置( pbest_i^t )、整个种群中最好的粒子的位置( gbest_i^t )以及粒子上一代的速度(v_i^t )计算出来，具体的表达式如下：

其中c₁是个体学习因子，表示粒子对过去自身经验的依赖，而c₂为全局学习因子，表示粒子对整个种群的信赖程度。早期研究表明c₁和c₂的值为2时，可以给出较好的优化结果。r₁和r₂是两个分布在0 和1 之间的随机数。这两个随机数可以保证我们的方法收敛到全局最小值点而不会陷入某个极小值点。显然粒子在搜索空间中的移动受每个粒子过去经验和种群经验的影响，追随能量最低的粒子向全局最小点移动。ω是惯性权重，当其取较大值的时候能够提高算法全局搜索的能力，而当取较小值的时候可以提高算法搜索的精度。在CALYPSO 方法中，随着迭代次数(iter)的增加，按照下面的公式，ω动态地从0.9 变到0.4：

其中ω_max 为0.9， ω_min 为0.4， iter_max为最大演化代数。

目前CALYPSO 方法中采用了全局和局域两个版本的粒子群优化算法。如图5(b)，(c)所示，在全局粒子群优化算法中，所有粒子追寻目前最优的粒子运动，具有快速收敛的特点，适用于势能面简单的体系(原子数约小于30)；局域粒子群优化算法将势能面分成不同的区域，粒子追寻所在区域中的最优粒子进行探索。局域粒子群优化算法可以看成是多个信息可以共享的全局粒子群优化算法的组合，具有更强的抗过早熟能力，适用于势能面复杂体系(模拟晶胞中的原子数大于30) 。

4 CALYPSO方法的应用

CALYPSO 方法的有效性已经通过对近百种已知实验结构的测试得到了证实，并且在科研实践中得到了进一步的检验，自2010 年软件发布以来，7 年的时间里，国内外同行利用CALYPSO软件已经在Nature Chem.，Nature Commun.，PRL，PNAS，JACS 等国际顶级期刊发表了400 多篇论文，在物理、化学、材料和地球科学等领域解决了若干长期无法解决的科学难题，典型工作包括：

(1)破解了单质锂的高压半导体相的结构难题。作为高压研究的模型体系，单质锂在高压下的相变极为复杂，其高压下的结构问题一直是高压研究的焦点。自从2002 年实验发现单质锂在60万大气压以上存在新的高压相以来，科学家们就开展了深入的理论和实验研究。但由于理论技术和实验条件的限制，新的高压相的结构一直没有得到解决。2009 年3 月，日本科学家在《自然》期刊上发表了单质锂的高压电学测量结果，令人意想不到的是，单质锂的新的高压相竟然是半导体，这更加激发了人们对其结构的研究热情。随后，国外几个研究小组分别利用近几年发展起来的基因算法和随机搜索算法等晶体结构预测技术模拟了锂的半导体相，但研究结果争议很大。针对上述问题，我们利用CALYPSO 方法对锂的高压半导体相结构进行了系统的探索，预言了一个晶体学单胞内含有40 个原子的复杂底心正交Aba2-40 结构(Pearson 符号，oC40，图6(a)) 。该结构的能量远远低于前期基因算法和随机搜索算法所获得的结构。oC40 结构中，锂的价电子由于受到芯电子的排斥完全局域到晶格间隙之中，失去了自由电子的特性，金属锂变成了半导体。在我们的工作发表不久后，英国爱丁堡大学的Guillaume 等人独立报道了锂的高压单晶X射线衍射实验数据，实验发现锂的高压半导体相的确具有oC40 结构。英国爱丁堡大学的Margues等人随后通过实验和理论相结合的方法进一步确定oC40 结构就是锂的半导体相结构，支持了我们的理论预言。

图6 (a)单质锂的高压半导体相oC40 结构；(b)地核环境下XeFe₃的晶体结构；(c)高压下H₂S金属相的晶体结构；(d)高压下CaH₆的晶体结构

(2)提出地核压力与温度条件下氙和铁/镍的化学反应。氙气作为惰性气体家族中的一员在大气层中几乎绝迹，与氩气和氪气等其他惰性气体相比，90%以上的氙气都不知所踪，这在科学上被称为“氙气的消失之谜”。氙气是否会储存在占据地球总质量三分之一的地核内部一直备受关注。地核的主要成分是铁和镍，如果氙气储存于地核中，就必须与铁或者镍在地核的压强和温度环境下(即360 GPa 和6000 K)发生化学反应，形成稳定的化合物。1997 年《科学》期刊发表的理论和实验合作论文否定了氙气和铁发生反应的可能性。研究工作发表后的17 年间，科学界据此公认氙气不可能储存在地核内部。然而，《科学》文章的结论是以人为假定的铁—氙化合物的结构为基础，若假定错误，氙气和铁不反应的结论就不再成立。我们利用CALYPSO 方法对氙和铁/镍在地核环境下的结构进行了系统研究，提出了全新的铁/镍—氙化合物的结构形式(图6(b))，构筑了铁/镍—氙化合物的高温—高压相图，首次给出了氙气和铁/镍在地核环境下发生化学反应的证据，提出了氙气被捕捉在地核内部的可能性。常压条件下，惰性气体氙难以和其他元素反应形成化合物；而铁/镍容易氧化失去电子，一般是带正电的还原剂。令人意想不到的是，高压下氙不仅和铁/镍发生了反应，而且电子从氙原子转移到铁/镍原子上，从而竟使铁/镍成为了罕见的带负电的氧化剂。这种高压下非常规的电子转移现象，诱导了氙气与铁/镍之间的化学反应。

(3)H₂S 和CaH₆的高压相结构与超导电性。超导领域的一个共识是BCS传统超导体的超导温度不可能超过40 K(McMillan 极限)。美国科学家Ashcroft 在2004 年提出，富含氢化合物一旦在高压下金属化就可能具有较高的超导温度。然而富含氢化合物种类繁多，实验科学家无法寻找到合适的富含氢化合物来开展超导实验研究。此时，通过理论结构预测方法，快速、经济地寻找最有希望的高温超导体对实验研究具有重要的指导意义。我们利用CALYPSO 方法，对H₂S 在高压下的结构进行了系统的探索，首次提出高压金属相H₂S (图6(c))是潜在的高温超导体，在160 GPa压力下，超导转变温度可以达到80 K。受到此项结果的启发，德国马克思—普朗克研究所Eremets 等人开展了H₂S 的高压超导实验，证实了我们的理论预言，并创造了203 K的超导温度新纪录。此外，我们还利用CALYPSO 方法，以突破常规化学计量比的新型碱土金属氢化物为研究对象，预言了高含氢量的CaH₆结构(图6(d)) ，其中氢形成了极其独特的二十面体笼型单元，笼型单元的H—H原子间形成弱的共价键，改变了含氢化合物中的只能形成氢分子和单原子氢的传统认识，丰富了化学成键理论。基于BCS 理论，我们预言150 万大气压下CaH6的理论超导转变温度达到220 K。

5 总结和展望

近年来，随着智能全局优化算法的引入，国内外科学家先后发展了多种理论结构预测方法，这些方法原理不尽相同、功能各具特色，已成为物质结构研究不可或缺的工具，被广泛应用于结构现象丰富的研究领域，在物理、化学、材料、地球等科学领域，解决了大量长期无法解决的科学难题。

然而，理论结构预测方法的发展还远没有结束。不难发现，现有结构预测方法还只适用于原子数较少(百原子以内)的简单体系，鲜有成功应用于大尺寸(百原子以上)复杂体系的算例。复杂大体系结构预测的困难主要体现在以下两个方面：(1)势能面的复杂程度随着体系原子数目的增加而迅速增加，其维度线性增加，极小值点个数以指数形式增长，现有结构预测方法在处理复杂体系势能面时均面临取样效率严重下降的问题；(2)计算资源的限制。现有较为准确的能量评估方法(如基于密度泛函理论的第一性原理方法)的计算量随体系电子数目的增加呈3 次指数变化关系，复杂体系结构预测往往需要对大量候选结构进行能量评估，因此计算成本极其昂贵。目前发展针对复杂体系的结构预测方法还面临着困难和挑战，但随着全局优化算法的不断发展、计算机计算能力的提升和基本物理化学理论的逐渐完善，实现复杂体系的结构预测指日可待。我们相信结构预测方法必将在更大、更复杂物质结构的研究中发挥重要作用。