王丽莉,熊鹰,谢炜宇,牛亮亮,张朝阳 分子动力学的核心在于求解正则运动方程,分子发展了许多晶体结构预测的方法,为研究含能晶体的 内和分子间的作用力都包含在分子力场中,因此计算结构及其性质注入活力。2016年,多个单位参加了对 的准确性依赖所使用的力场4-。如果缺乏完全适有机晶体结构的第六次盲测2,分子结构均来自第 用的力场参数,无法准确描述分子间的相互作用,致使原理计算和剑桥结构数据库,晶体结构通过随机搜索 晶体结构预测存在较大偏差,要提高晶体密度计算的法、遗传算法或模拟退火法产生。 准确性需自研发经验势。 3.2.1模拟退火法 313半经验方法 模拟退火法是S. Kirkpatrick等于1983年基于 半经验方法通过忽略分子积分计算或引进经验参金属退火机理而建立的全局最优化方法,通过对给定 量对 Hartree-Fock方程的近似求解,如电荷自恰密度初始状态进行模拟加温,生成一系列关联的新解并评 泛函紧束缚( Self-consistent Charge Density Function-估与初始解的能量差,再根据 Metropolis接受准则进 al Tight-binding, SCC-DFTB)方法是在DFT基础上采行结构筛选,得到一系列给定温度下允许存在新 用非正交基的非正交化方式来对紧束缚( Tight Bind-解。例如在温度T时,由当前状态i产生新状态能量 nTB)理论做修正,计算速度比传统DFT快两三个分别为E和E,若E>E则接受新状态为当前状态,否 数量级,核心是对哈密顿矩阵H的求解和对电子波函则,以概率P来接受状态,其中k为 Boltzmann常数。 数的计算,能量表达式如下所示 p=exp[ -(E-E)/k] Em=∑{的HP)+∑E (6) 在模拟退火的开始阶段通常以升高温度来实现初 式中,Oc表示电子占据数,i表示具体的价电子,EA始解在能量曲面上的上坡游走,克服势垒走向周围的 表示原子核A和B之间的排斥能,k。电子波函数主 局域最小值;这些走离初始解位置的新解通常都处于 要通过基函数的线性组合来求解: 能量面中的坡上,需要进一步通过降温处理让它们再 v(x,t)=∑va,t)(x) (7)次在能量曲面上进行下陂游走,走向周围的局域态 升温和降温时,温度变化遵循如下关系式 式中,l是电子的空间坐标,a=x,y,z表示空间维度 Tnew= Told x(1.0+T) 大大降低了哈密顿矩阵的求解难度,可用于计算晶体m=mx(0-7) (10) 结构和密度0,但该方法没有针对硼元素的经验参式中,Tm和Ta分别是变温前后的温度,K;T和T分 数,无法模拟极端高温高压条件下的物理化学变化。别是自定义的加热因子和冷却因子。 Holden用半经验AM1( Austin Model number1) Yin153和G.Z.Zhao5分别用模拟退火技术预 方法对堆积结构进行优化,预测晶胞参数。 测了1′-氨基-1H-1,5′双四唑(ATT)及其溴氰衍生物 半经验方法虽然可将计算规模增大,但还不能有和2,4,6三硝基-1,3,5-三氮1,3,5-三氧化六元氮杂 效以及长时间模拟原子数量非常庞大的系统或者结构环( TNTATO)的晶体结构。刘英哲等”基于B3LY 更加灵活的体系,比如复杂分子表面势能以及弱键能方法和6-31G(d)基组建立了适用于全氮结构的力场, 等;计算的可靠性与参数获得的途径密切相关,参数拟用 Monte-Carlo方法进行空间堆积,预测晶体结构并计 合大多依赖于DFT计算;由于采用了一系列近似,应算N4、N4、Na、N0、N12五种笼型全氮分子的晶体密度分 用于含有强吸电子基团(如硝基)的体系时,由于电荷别是1.81,2.08,2.47,2.46,2.57gcm。.F 密度不对称,计算误差会大大增加 Moxnes用 Polymorph Predictor方法和 COMPASS 3.2势能面搜索预测晶体结构 力场计算晶体密度并与实验数据比较(见表1),对大 晶体结构预测是在固定化学组分情况下确定最有多数空间群,模拟退火法高估了密度,平均绝对误差 可能的原子排列方式和晶胞参数,这一过程中用到自0.07~0.08gcm3。 然界中占据统治地位的能量最低原理,即在固定组分 模拟退火法是早期含能晶体结构预测常用方法, 的状态子空间中搜索处于能量曲面极小值的排列方这种基于初始结构演化的跃迁势垒方法对人为选择的 式。从理论上预测晶体结构一直是化学、物理和材料初始结构具有依赖性,考虑空间群不够宽泛,难以克服 科学中的一大难题,随着理论方法和计算技术的完善,退火过程的高能量势垒,需要进行多次迭代直到所得 Chinese Joumal of Energetic Materials. Vol28, No 1, 2020(1-12) 含能材料 www.energetic-materials.orgcnChinese Journal of Energetic Materials，Vol.28, No.1, 2020（1-12） 含能材料 www.energetic-materials.org.cn 王丽莉，熊鹰，谢炜宇，牛亮亮，张朝阳 分子动力学的核心在于求解正则运动方程，分子 内和分子间的作用力都包含在分子力场中，因此计算 的准确性依赖所使用的力场［48-49］。如果缺乏完全适 用的力场参数，无法准确描述分子间的相互作用，致使 晶体结构预测存在较大偏差，要提高晶体密度计算的 准确性需自研发经验势。 3.1.3 半经验方法 半经验方法通过忽略分子积分计算或引进经验参 量对 Hartree⁃Fock 方程的近似求解，如电荷自恰密度 泛函紧束缚（Self⁃consistent Charge Density Function⁃ al Tight⁃binding，SCC⁃DFTB）方法是在 DFT 基础上采 用非正交基的非正交化方式来对紧束缚（Tight Bind⁃ ing，TB）理论做修正，计算速度比传统 DFT 快两三个 数量级，核心是对哈密顿矩阵 H 的求解和对电子波函 数 λ的计算，能量表达式如下所示： EDFTB =∑i occ ψi | H | 0 ψi + ∑A ≠ B atoms E AB rep （6） 式中，occ 表示电子占据数，i 表示具体的价电子，E AB rep 表示原子核 A 和 B 之间的排斥能，kJ。电子波函数主 要通过基函数的线性组合来求解： Ψi (x，t ) =∑la Ψi (la，t )Φla (x ) （7） 式中，l 是电子的空间坐标，a = x，y，z 表示空间维度。 DFTB 方法采用经验数值代替原本复杂的积分方程， 大大降低了哈密顿矩阵的求解难度，可用于计算晶体 结构和密度［50］，但该方法没有针对硼元素的经验参 数，无法模拟极端高温高压条件下的物理化学变化。 J. R. Holden［51］ 用半经验AM1（Austin Model Number 1） 方法对堆积结构进行优化，预测晶胞参数。 半经验方法虽然可将计算规模增大，但还不能有 效以及长时间模拟原子数量非常庞大的系统或者结构 更加灵活的体系，比如复杂分子表面势能以及弱键能 等；计算的可靠性与参数获得的途径密切相关，参数拟 合大多依赖于 DFT 计算；由于采用了一系列近似，应 用于含有强吸电子基团（如硝基）的体系时，由于电荷 密度不对称，计算误差会大大增加。 3.2 势能面搜索预测晶体结构 晶体结构预测是在固定化学组分情况下确定最有 可能的原子排列方式和晶胞参数，这一过程中用到自 然界中占据统治地位的能量最低原理，即在固定组分 的状态子空间中搜索处于能量曲面极小值的排列方 式。从理论上预测晶体结构一直是化学、物理和材料 科学中的一大难题，随着理论方法和计算技术的完善， 发展了许多晶体结构预测的方法，为研究含能晶体的 结构及其性质注入活力。2016 年，多个单位参加了对 有机晶体结构的第六次盲测［52］ ，分子结构均来自第一 原理计算和剑桥结构数据库，晶体结构通过随机搜索 法、遗传算法或模拟退火法产生。 3.2.1 模拟退火法 模拟退火法是 S. Kirkpatrick 等［53］ 于 1983 年基于 金属退火机理而建立的全局最优化方法，通过对给定 初始状态进行模拟加温，生成一系列关联的新解并评 估与初始解的能量差，再根据 Metropolis 接受准则进 行结构筛选［54］，得到一系列给定温度下允许存在新 解。例如在温度 T 时，由当前状态 i 产生新状态 j，能量 分别为 Ei 和 Ej ，若 Ei > Ej 则接受新状态为当前状态，否 则，以概率 p 来接受状态，其中 k 为 Boltzmann 常数。 p = exp[ - (Ej - Ei) /kt] （8） 在模拟退火的开始阶段通常以升高温度来实现初 始解在能量曲面上的上坡游走，克服势垒走向周围的 局域最小值；这些走离初始解位置的新解通常都处于 能量面中的坡上，需要进一步通过降温处理让它们再 次在能量曲面上进行下陂游走，走向周围的局域态。 升温和降温时，温度变化遵循如下关系式： Tnew = Tol d × (1.0 + Th ) （9） Tnew = Tol d × (1.0 - Tc ) （10） 式中，Tnew 和 Tol d 分别是变温前后的温度，K；Th 和 Tc 分 别是自定义的加热因子和冷却因子。 J. Yin［55］和 G. Z. Zhao［56］分别用模拟退火技术预 测了 1′⁃氨基⁃1′H⁃1，5′⁃双四唑（ATT）及其溴氰衍生物 和 2，4，6⁃三硝基⁃1，3，5⁃三氮⁃1，3，5⁃三氧化六元氮杂 环（TNTATO）的晶体结构。刘英哲等［57］基于 B3LYP 方法和 6⁃31G（d）基组建立了适用于全氮结构的力场， 用 Monte⁃Carlo 方法进行空间堆积，预测晶体结构并计 算 N4、N6、N8、N10、N12五种笼型全氮分子的晶体密度分 别 是 1.81，2.08，2.47，2.46，2.57 g·cm-3 。 J. F. Moxnes［58］用 Polymorph Predictor 方 法 和 COMPASS 力场计算晶体密度并与实验数据比较（见表 1），对大 多数空间群，模拟退火法高估了密度，平均绝对误差 0.07~0.08 g·cm-3 。 模拟退火法是早期含能晶体结构预测常用方法， 这种基于初始结构演化的跃迁势垒方法对人为选择的 初始结构具有依赖性，考虑空间群不够宽泛，难以克服 退火过程的高能量势垒，需要进行多次迭代直到所得 4