Al/MoO3体系含能材料反应力场关键参数的研究

高 超 胡丽娟 张巧凤 刘 珂 梁婉怡 谢耀平

(上海大学材料科学与工程学院材料研究所、微结构重点实验室，上海 200072)

亚稳态分子间复合物(metastable intermolecular composites, MIC)是一种由纳米尺度的燃料(通常为Al)和氧化剂(MoO3、CuO、Fe2O3等)颗粒组成的新型含能材料。未加载环境下，MIC处于亚稳态，在外界载荷(如激光、热诱导等)下会发生反应形成一种更稳定的合金或复合物。MIC具有燃烧效率高、能量释放快等特点。如Al和MoO3构成的MIC含能材料发生反应时，反应传播速率在600～1 000 m/s之间，能量释放为1.12 kcal/g[1]，这可以与常见的硝胺基爆炸物相媲美。

为了探讨MIC燃烧的反应机制，许多学者对其做了研究。Asay等[2]试验研究了MIC材料Al/MoO3体系的燃烧过程，认为对流机制是MIC材料燃烧的反应机制。Dean等[3]以Al/CuO和Al/NiO作为试验材料，推测MIC的燃烧机制为熔化- 爆裂机制。但由于MIC燃烧反应一般都在远离平衡态或者是极端条件下发生，这给试验研究带来了一定困难，从而难以确定MIC材料燃烧的反应机制。分子动力学(molecular dynamics，MD)作为一种理论研究手段可以弥补试验方法在时间和空间尺度上的局限性，能在原子尺度上直观地观测到反应过程中各个原子的运动状况，进而研究MIC燃烧的反应机制。

MD采用经验的势函数描述原子的运动，因此找到合理的势函数来准确描述MIC的反应过程是必要的。对势可以较好地描述在平衡状态的体系，嵌入原子势[4](embedded- atom method，EAM)可以对金属体系进行较好的描述。但它们都不能描述键的断裂和生成，不适用于描述MIC的燃烧反应。由Van Duin等[5]开发的势函数——反应力场(reactive force filed，ReaxFF)可以精确描述复杂化学反应过程中原子间的相互作用力，能很好地模拟键的断裂与形成过程，是研究MIC燃烧反应的有效方法。目前，运用ReaxFF已经成功描述了有机物[5]、金属氧化态[6]、Al/Al2O3界面[7]等结构的物理和化学性质。由于必须获取材料的ReaxFF参数才可以运用ReaxFF对其进行MD模拟，为了获取适用于MIC材料的ReaxFF参数，需要对ReaxFF参数进行拟合。通过将拟合获得的ReaxFF参数引入ReaxFF，才可以将ReaxFF运用于MIC体系的计算模拟。

本文选取典型的MIC材料Al/MoO3体系，拟合其ReaxFF参数，并验证ReaxFF对该体系的适用情况。通过分析ReaxFF势函数和力场拟合方法，综合运用拟合好的ReaxFF势函数计算Al/MoO3体系的一系列性质，并将其与第一性原理计算值或试验值进行对比，验证ReaxFF对Al/MoO3体系的适用情况。

1 ReaxFF反应力场模型的建立

1.1 ReaxFF反应力场

键级是ReaxFF的核心概念，它能很好地反映化学键的信息，从而能够描述键的断裂和形成过程。如式(1)～式(3)即为ReaxFF计算系统能量的表达式[5]。

Esystem=Ebond+Eangle+Etorsion+

Evdwaals+ECoulomb

(1)

(2)

(3)

1.2 ReaxFF拟合方法

ReaxFF参数决定着MIC的反应过程分子动力学模拟的准确度及精度，因此拟合出合理可靠的ReaxFF参数可以说是MIC的反应过程分子动力学模拟的一项核心，并具有一定挑战性。ReaxFF参数的拟合主要分为三部分。首先，获得合适的训练基(training sets)。训练基是应用于ReaxFF参数拟合过程中的一组输入数据，为第一性原理计算值或试验值，组成训练基的数据有键能、键角能、生成焓、晶格常数、状态方程、空位形成能等。下一步是反应力场的拟合，由于ReaxFF有大量的参数需要拟合，本文采用单参数连续优化方法(successive one- parameter search method)[9]依次优化一套力场参数中的每个参数，每次优化使得ReaxFF计算值与训练基中标准值的偏差平方和最小，此时得到的参数即为该参数的最佳值。接着通过该自洽方法依次循环优化每个参数，最终得到一套最佳的力场参数。最后，对该套力场参数进行验证，确定其合理性。

2 结果与讨论

为了拟合出适用于MIC材料Al/MoO3体系的ReaxFF参数，以前期研究开发的Al/Mo/O ReaxFF参数[10]为基础，开发适用于Al/MoO3体系的Mo- O二元与Al- Mo- O三元参数，从而完善该ReaxFF参数，使其可以描述MIC材料Al/MoO3体系。

以Mo- O二元结构MoO2、MoO3、Mo4O11、Mo8O23、Mo9O26和Al- Mo- O三元结构Al8Mo12O48、Al16Mo24O96为研究对象，运用第一性原理(density functional theory，DFT)计算这些结构的几何参数和能量作为训练基来拟合ReaxFF参数。DFT计算采用了VASP[11]程序包，利用投影缀加平面波方法(PAW)来描述电子与原子间作用，交换关联泛函采用GGA- PW91处理交换关联能，布里渊区内则采用MP方法自动产生K点。分子动力学模拟则采用大型分子动力软件LAMMPS[12]进行。

2.1 能量变化

由于DFT方法可以较精确地描述结构的能量，可以视其状态方程(equation of states，EOS)曲线为真实值。因此采用拟合完毕的Al/Mo/O ReaxFF分别计算了Mo- O二元结构和Al- Mo- O三元结构能量随晶格常数变化的EOS，将其与DFT计算结果作对比，以验证ReaxFF对Al/MoO3体系的适用情况，如图1所示。由图1可以看出，ReaxFF势对Al/MoO3体系的各种结构的描述较好，能够准确区分出Al/MoO3体系的不同结构，并较准确地描述各个结构的能量变化。需要指出的是，图1(d)中，Al8Mo12O48和Al16Mo24O96的EOS曲线与DFT计算的结果有些许偏差，这是由于ReaxFF可以精确描述复杂化学反应过程中原子间的相互作用力，但对金属键的描述能力不如DFT那样准确。比如宋文雄等[13]使用ReaxFF对金属Al以及Narayanan等[14]使用ReaxFF对Li/Al/Si/O体系含有金属键结构的描述都出现了一定的偏差。但这种情况并不影响力场的应用，因为力场可以相对准确描述出Al8Mo12O48和Al16- Mo24O96结构的整体能量变化。

图1 运用DFT和ReaxFF势函数计算的Al/MoO3体系结构的状态方程(EOS)曲线Fig.1 Equation of states of Al/MoO3 structures calculated by using DFT and ReaxFF

为了验证拟合好的ReaxFF对Al/MoO3体系表面的描述能力，以Al/MoO3体系晶体结构中原子排布清晰的表面为研究对象，分别用DFT和ReaxFF计算了这些表面的表面能，如表1所示。由表1可以看出，ReaxFF计算结果与DFT计算结果一致，说明拟合后的ReaxFF能够描述Al/MoO3体系结构的表面能。

2.2 晶体结构

表1 运用DFT与ReaxFF方法计算得到的材料表面能Table 1 Surface energy of some structures calculated by using DFT and ReaxFF eV

运用DFT和ReaxFF计算了Mo- O二元结构和Al- Mo- O三元结构的晶格常数与结合能，如表2所示。运用ReaxFF计算的不同结构的晶格常数和结合能值与DFT计算值或试验值相差均小于5%，说明拟合得到的力场参数较好地重现了Al/MoO3体系各个结构的晶格参数和结合能，可以较准确地描述MIC材料Al/MoO3体系的晶体结构。ReaxFF计算得到的Al8Mo12O48和Al16- Mo24O96结构的结合能分别为3.61和3.60 eV，这与DFT计算得到的2.89和2.87 eV相比有一定偏差，这是由于ReaxFF对金属键的描述能力较弱造成的。另外，可以发现，DFT计算的晶格常数比试验值要大一些，这是由于DFT计算采用的是GGA势函数，GGA势函数在描述氧化物时会引起晶格常数偏大。Scanlon等[15]使用GGA势函数对MoO2和MoO3结构的晶格常数进行了DFT计算，发现这些结构的晶格常数均大于试验值，与本文结果一致。

表2 运用DFT与ReaxFF方法计算得到的晶格常数和结合能及其试验值Table 2 Lattice constants and binding energies calculated by DFT and ReaxFF, and their experimental data

2.3 弹性性质

晶体的状态方程、热容、德拜温度、熔点等固态性质都与其弹性性质密切相关。Yun等[17]开发了Fe/Al/Ni合金的ReaxFF势函数，并使用其计算了Fe/Al/Ni合金的弹性常数以验证ReaxFF适用性。Zhang等[7]运用了ReaxFF计算了α- Al2O3的弹性常数，发现ReaxFF能够准确描述α- Al2O3的弹性性能。可见，弹性性能已经成为验证ReaxFF适用性的有效手段。通过计算弹性常数，可以验证拟合好的ReaxFF参数对晶体弹性性能的描述情况。根据弹性应变能与应变的关系来计算弹性常数，如式(4)所示。

(4)

对晶体施加不同的应变ei、ej，根据晶体施加应变前后应变能的变化△E及晶体体积V即可求得晶体的弹性常数Cij。分别使用DFT与ReaxFF计算正交晶系MoO3、Mo4O11和Al8Mo12O48结构的弹性常数，以验证拟合后的力场参数对晶体力学稳定性的描述情况。正交晶系的弹性常数为9个，分别为C11、C12、C13、C23、C22、C33、C44、C55、C66。表3给出了用DFT和ReaxFF方法计算得到各个结构的弹性常数。由表3可以看出，对于不同结构，ReaxFF计算的弹性常数与DFT计算结果有一定偏差。由于ReaxFF添加了Morse函数来描述范德瓦尔斯力，而Morse函数考虑两个原子间的相互作用，但不会考虑键角能、二面角能等多体相互作用。另一方面，ReaxFF又通过计算键级来描述多体相互作用。在拟合ReaxFF参数时，可能并未处理好Morse函数表示的范德华作用与键级描述的多体相互作用，导致ReaxFF计算的弹性常数偏大。但是ReaxFF计算的弹性常数与DFT计算结果的数量级一致，而且不同弹性常数间大小关系也几乎与DFT计算结果一致。式(5)～式(7)给出了正交晶系的力学稳定性判据，采用DFT和ReaxFF计算得到的弹性常数都满足此判据，表明这两种方法得到的晶体结构是稳定的，这与事实相符。因此ReaxFF预测Al/MoO3体系的弹性常数值虽然大于DFT计算结果，但是ReaxFF可以准确地描述Al/MoO3体系的力学稳定性。

C11>0;C22>0;C33>0;

C44>0;C55>0;C66>0

(5)

C11+C22-2C12>0;C11+C33-2C13>0;

C22+C33-2C23>0

(6)

C11+C22+C33+2(C12+C13+C23)>0

(7)

表3 运用DFT与ReaxFF方法计算得到的弹性常数Table 3 Elastic constants of some structures calculated by DFT and ReaxFF kbar

3 结论

为了使用ReaxFF反应力场对MIC体系燃烧反应机制进行理论研究，以常见MIC材料Al/MoO3体系为研究对象，对其ReaxFF参数的Mo- O二元与Al- Mo- O三元参数进行了拟合优化，并得到了一套较为完善的力场参数。运用该力场参数计算了Mo- O和Al- Mo- O结构的物理性质，如状态方程、晶格常数、表面能和弹性常数。并将计算结果与第一性原理(DFT)计算结果进行对比，以验证该力场对Al/MoO3体系的适用性。

拟合好的Al/Mo/O ReaxFF势对Mo- O和Al- Mo- O结构的状态方程(EOS)和晶格常数描述较好，能够较准确地描述各个结构的能量变化、晶体结构和表面能。在描述晶体弹性性能方面，拟合好的ReaxFF计算的Al/MoO3体系的弹性常数值，虽与DFT计算结果有一定偏差，但可以准确描述Al/MoO3体系的力学稳定性。这说明拟合好的ReaxFF可以较准确地描述Al/MoO3体系的物理性质，可作为研究MIC材料Al/MoO3燃烧反应机制的分子动力学势函数。

[1] BOCKMON B S, PANTOYA M L, SON S F, et al. Combustion velocities and propagation mechanisms of metastable interstitial composites [J]. Journal of Applied Physics, 2005,98(6): 064903.

[2] ASAY B W, SON S F, BUSSE J R, et al. Ignition characteristics of metastable intermolecular composites [J]. Propellants Explosives Pyrotechnics, 2004, 29(4): 216- 219.

[3] DEAN S W, PANTOYA M L, GASH A E, et al. Enhanced convective heat transfer in nongas generating nanoparticle thermites [J]. Journal of Heat Transfer, 2010, 132(11): 2201- 2210.

[4] DAW M S, BASKES M I. Embedded- atom method: Derivation and application to impurities, surfaces, and other defects in metals [J]. Physical Review B, 1984, 29(12): 6443- 6453.

[5] VAN DUIN A C T, DASGUPTA S, LORANT F, et al. ReaxFF: a reactive force field for hydrocarbons [J]. The Journal of Physical Chemistry A, 2001, 105(41): 9396- 9409.

[6] ARYANPOUR M, VAN DUIN A C T, KUBICKI J D. Development of a reactive force field for iron- oxyhydroxide systems[J]. The Journal of Physical Chemistry A,2010,114(21): 6298- 6307.

[7] ZHANG Q, CAGIN T, DUIN A V, et al. Adhesion and nonwetting- wetting transition in the Al/α- Al2O3interface [J]. Physical Review B, 2004, 69(4): 1129- 1133.

[8] MORTIER W J, GHOSH S K, SHANKAR S. Electronegativity- equalization method for the calculation of atomic charges in molecules [J]. Journal of the American Chemical Society, 1986, 108(15): 4315- 4320.

[9] VAN DUIN A C T, JAN M, DE GRAAF B. Delft molecular mechanics: a new approach to hydrocarbon force fields. Inclusion of a geometry- dependent charge calculation [J]. Journal of the Chemical Society, Faraday Transactions, 1994, 90(19): 2881- 2895.

[10]SONG W X, ZHAO S J. Development of the ReaxFF reactive force field for aluminum- molybdenum alloy [J]. Journal of Materials Research, 2013, 28(9): 1155- 1164.

[11] KRESSE G, FURTHMULLER J. Efficiency of ab- initio total energy calculations for metals and semiconductors using a plane- wave basis set [J]. Computational Materials Science, 1996, 6(1): 15- 50.

[12] PLIMPTON S. Fast parallel algorithms for short- range molecular dynamics [J]. Journal of Computational Physics, 1995, 117(1): 1- 19.

[13] 宋文雄,赵世金.ReaxFF反应力场在金属铝中的适用性[J].含能材料,2012,20(5): 571- 574.

[14] NARAYANAN B, VAN DUIN A C T, KAPPES B B, et al. A reactive force field for lithium- aluminum silicates with applications to eucryptite phases [J]. Modelling and Simulation in Materials Science and Engineering, 2011, 20(1):386- 401.

[15] SCANLON D O, WATSON G W, PAYNE D, et al. Theoretical and experimental study of the electronic structures of MoO3and MoO2[J]. The Journal of Physical Chemistry C, 2010, 114(10): 4636- 4645.

[16] BREWER L, LAMOREAUX R H,FERRO R. The Mo- O System (Molybdenum- Oxygen) [J]. Journal of Phase Equilibria, 1980, 1(2): 85- 89.

[17] YUN K S, KWAK H, ZOU C, et al. Development and Validation of a ReaxFF Reactive Force Field for Fe/Al/Ni Alloys: Molecular Dynamics Study of Elastic Constants, Diffusion, and Segregation [J]. Journal of Physical Chemistry A, 2012, 116(49): 12163- 12174.