微观理论计算和模拟方法

1.3.2  分子动力学模拟及其基本原理

本论文采用的主要是分子动力学模拟(Molecular Dynamics Simulation，MD Simulation)，它是时下最广泛采用的计算庞大复杂系统的计算方法。自1970年起，基于分子力学的迅速发展，人们系统地建立了许多适用于聚合物、生化分子体系、金属与非金属材料的力场，使得计算复杂体系的结构、热力学与光谱性质的能力和精确性大大提升；MD模拟是应用这些力场及根据牛顿运动力学原理所发展的计算方法，此方法的优点为：⑴精确度高，⑵可同时获得系统的热力学与动态统计资料，⑶可以广泛地应用于探讨各种系统的各类特性；经过许多改进，MD的计算方法现已日趋成熟，由于其计算能力强，能满足各类问题的需求[16,17]。

1957年，Alder等在硬球模型下首次运用MD方法研究气、液状态方程, 开创了运用MD方法研究物质宏观性质的先河；1972年，Less等将此方法的运用领域推广到具有速度梯度的非平衡体系；1980年，Andersen等在前人研究基础上开创了恒压分子动力学方法；1983年，Gillan等将MD方法扩展到存在温度梯度的非平衡系统, 从而建立了非平衡系统MD方法体系；1984年，Nose等创建了恒温MD方法；1985年，Car等针对半导体和金属势函数模型化比较困难的问题，提出了将电子论与MD方法有机统一起来的设想，创建了第一性原理MD方法；1991年，Cagin等进一步提出了适用于处理吸附类问题的巨正则系综MD方法；20世纪末，伴随着计算机技术的迅猛发展,加之多体势函数的提出与发展,使得MD模拟技术有了长足的进步[18]。

在经典MD模拟中，假设原子运动遵循牛顿运动方程，也就是说原子的运动按特定轨迹进行；假定原子核运动的量子效应可以忽略，并且绝热近似严格成立，即设任一时刻下电子均处于相应原子结构的基态[16,17]。为了进行MD模拟需要知道各个原子间正确的相互作用势，在经典MD模拟中，用经验势能函数表示原子间的相互作用势。

通常势能可分为分子内部势能(imt)与分子间(或分子内)原子间的非键结范德瓦耳斯作用(vdW)两部分，即：

U=UvdW+Uint                                         (1-1)

在经典力学中，系统中任一原子i所受的力等于势能的梯度：

                (1-2)

依此，依据牛顿运动定律可求得i原子的加速度为：

                              (1-3)

将牛顿运动定律方程式(1-4)对时间进行积分，可以预测出i原子经过时间t后的位置、速度与加速度。

                        (1-4)

式(1-4)中，粒子的位置和速度分别用 和 来表示，上标“0”为各物理量的初始值[16]。

在一定力场下进行MD模拟，首先由系统中各原子的位置计算系统势能，再由式(1-2)和式(1-3)计算出系统中各原子所受的力与加速度，然后令式(1-4)中t=Δt，Δt表示极短的时间间隔。重复以上步骤。这样，随着时间的推移，原子位置不断向前移动，从而可以得到体系中原子运动时间与原子位置的关系，即运动轨迹[17]；基于牛顿力学基本方程的MD模拟既能得到各原子的运动轨迹，又能进一步基于轨迹计算所需的各项性质，还能如同做实验一样进行各种观察，成为理论与实验的有效补充[19]。

1.3.3  分子动力学模拟在含能材料领域的应用

1993年，在第二十五届ICT国际会议上，Cumming[20]发表了关于HMX和PNMO相互作用的MM和MD模拟计算工作，开启了近代计算模拟化学用于含能材料分子间相互作用理论研究的先河。该文为PBX配方设计提供了一种崭新的研究工具。随后，关于运用MD模拟研究HMX单体炸药的报道越来越多。Manaa等[21]应用MD模拟研究了近似爆炸条件下HMX的热分解，其分解产物为CO2、CO、N2和H2O，气体产物的反应速率与采用热力学方法计算的最大浓度出现时间相吻合。Sewell[22]等运用MD模拟计算了不同HMX晶型的弹性系数和工程模量，模拟计算结果表明β-HMX晶体接近于各向同性的弹性体，刚性较强、属于脆性物质，温度变化对其力学性能的影响不大。Y. Kohno[23]通过设定不同的模拟条件，研究晶型对HMX初始热分解的影响，发现各晶型HMX初始热分解历程不同的原因在于其分子中N-N键键长的差异。马秀芳等[24]应用MD模拟研究了晶体缺陷对β-HMX晶体力学性能的影响，发现相对于HMX完美晶体，缺陷晶体的弹性系数和工程模量均有所下降。朱伟等[25,26]用MD模拟，在NVT系综下对不同配比的高氯酸铵和HMX所组成的二元混合体系进行模拟计算，研究得到PBX结构与感度的关系。