材料的计算机模拟方法课件

2022/11/141材料的计算机模拟方法姬广富冲击波物理与爆轰物理国防科技重点实验室内容材料计算方法以及最新的进展；含能材料的引发反应；静高压下材料的行为；材料的降解动力学模拟；气态CO在高温高压下的歧化反应；

HMX固体的性能理论模拟1计算方法简介在含能材料的模拟计算中，目前主要有三种方法：第一性原理（或从头算）的电子结构计算(abinitio)、半经验方法分子力学方法(MM)。

通过对材料的模拟计算可以得到如下的性质：生成热、键能和反应能、分子（或晶体）的能量和结构热化学稳定性)、过渡态的能量和结构（活化能）、反应路径以及动力学、机理、分子中的电荷分布（反应点）、取代效应、振动频率（红外和Raman光谱）、电子的激发跃迁（UV/可见光谱）、磁屏蔽效应（NMR谱）。。。。。。通过对这些性质的预测，在含能材料中有许多实用价值，如：可以研究合成路径、反应产物、爆轰产物状态方程、爆轰引发的机理等等。1.1从头算(abinitio)方法从头算方法在计算中是基于量子力学基本定律，没有引入经验或实验参数，仅使用5个常数（电子质量：9.10938215E-31Kg、电子电量：1.602176487E-19C、光速：c=2.99792458×108m/s、Planck常数：h=6.62606896×10-34J·s和波尔兹曼常数1.3806504E-23J/K），对化学体系利用满足严格数学近似的基组来解Schrödinger方程：动能是所有粒子的动能加和势能代表粒子之间的库仑相互作用（即：核与核之间、电子与电子之间的排斥以及电子与核之间的吸引）对电子而言，q=-e，对核而言，原子序数为Z：q=+Ze。通过对Schrödinger方程的形式观察，我们知道对于多于两个以上粒子的分子，严格解Schrödinger方程是不可能的，必须引入一些必要的合理近似，甚至一些经验参数，如非相对论解、Born-Oppenheimer近似、单粒子近似、原子轨道的线性组合等。对解出的Schrödinger方程是否合适，必须满足以下判据：这个解必须很好地被定义且能够用来说明分子结构和分子中电子的状态。分子的势能按照原子核的位移必须平滑连续可变。模型不能有偏离（即：假设化学键只存在于原子之间）。模型必须与体系大小相关（即解和相关误差必须与分子的大小成比例），模型必须是可变的（即：近似解必须提供体系能量真值的上界，因此近似解应具有最低能量值与真实的波函数相一致）。1.1.1基组一般来说，基组属于数学函数的一类，用来解微分方程。基组：G、PW、G+PW量化计算，基组一般约定为高斯基组来代表原子轨道，用来优化和产生所需体系的化学性质物理计算（DFT方法）,特别是金属的计算，经常用PW；GPW标准的从头算软件包都提供基组、当分子中包含较大的原子（Z>30）时，我们经常将原子的内层电子合并作近似处理（赝势）。带极化函数的基组，在基组中加入极化函数后可以允许原子核外电荷的不均匀位移，因此它可以改进对化学键的描述。极化函数可以描述那些和孤立原子相比更高的角动量量子数的轨道（如p-型函数用于H和He，d-型函数用于Z>2的原子），并将之加入到价电子壳层中。例如，6-31G(d)基组的构建是将6个d-型简单高斯函数加入到6-31G中来描述每个非氢原子；同样地，6-31G(d,p)是6-31G(d)对重原子,但加入了一组p-型高斯函数到氢原子或氦原子。将p-轨道加入到氢原子中对那些含有以氢原子为桥原子的体系而言特别有用。带弥散函数的基组，对于那些电子密度明显远离核的中心的原子或离子（如阴离子、孤电子对以及激发态），就需要用弥散函数来解决核对最外面电子束缚较弱的问题。一般来说，建议将弥散基组用于处理阴离子的电子亲和性、质子亲和性、反演势垒以及键角等。只对非氢原子加入弥散的s-和p-型高斯函数，弥散基组用加号“+”表示，如3-21+G；若对重原子和氢原子均加入弥散函数就用两个加号“++”表示。1.1.2电子相关方法由于HF方法为单粒子近似，忽略了电子之间的相互作用，因此会过低的估计键长同时会过高地估计振动频率，因此发展出了三种考虑电子相关的方法（或称post-HF方法）：组态相互作用（configuration-interaction，CI）Møller-Plesset(MP)微扰理论密度泛函理论（densityfunctionaltheory，DFT)对典型的分子而言，大部分的基态能量都源自于单电子HF的能量。MP微扰理论假设电子相关的效应很小，可以用很小的矫正（即微扰）来得到HF的解。MP方法假设分子真正的Hamiltonian量可分成两部分：（其中：代表单电子的贡献；代表组态相互作用的贡献）λ则代表能量和分子波函数的展开系数，即：；代入分子的Schrödinger方程中就可以得到：EHF=E(0)+E(1)+E(2)+…这就是所谓的MP2、MP3、MP4、MP5等等。Møller-Plesset微扰理论这种方法是将电子相互作用看作是电子密度的函数，现行的DFT方法通过Kohn-Sham方程将电子能量分成4个部分：E=ET

+EV

+EJ

+EXC，其中ET

代表电子的动能(来自于电子的运动)；EV

代表电子的势能（包括核与核之间的排斥能、核与电子之间的相互作用能）；EJ

代表电子之间的排斥能；EXC代表电子相关能。由于对EXC的定义不同便产生了不同的DFT方法，如：LDA、GGA、B3LYP、B3PW91等等。DFT方法可分成三类：局域密度近似（LDA）、广义梯度近似（GGA）和DFT与HF项结合的杂化方法：LDA方法只含有与电子密度有关的项，可以用来预测块材料的性质，但对于孤立分子的性能预测效果就不行；GGA方法同时考虑了电子密度和电子密度梯度，用它预测的键离解能与实际仅相差几个kJ/mol。但GGA方法却严重地过低估计了一些反应的活化位垒，这是因为这种方法忽略了电子的库仑子相互作用的原因；通过DFT与HF杂化的方法加以解决GGA方法产生的问题。密度泛函理论（DFT）1.1.3闭壳层和开壳层计算大多数分子具有电子单线态基态（电子的总自旋数为零），这时分子内电子完全成对地在限制在轨道内，这时就要用闭壳层自旋限制计算。对于那些具有非成对电子的体系，就要采用开壳层非自旋限制方法计算，计算中考虑电子的自旋状态（分为自旋向上α电子和自旋向下β电子）。

开壳层非限制性自旋计算适用的体系：具有奇数个电子的分子(例如一些离子和自由基)；激发态；分子有三线态或更高自旋的基态(例如O2)；电子对的分离过程(例如分解等)；共振体系，存在电子离域。自旋污染非限制性自旋波函数并不是自旋算符的本征函数，因此要对自旋态进行混合，从而导致自旋污染。在大多数开壳层系统中，自旋污染是可以忽略的。然而对于未成对电子离域的分子，自旋污染效应相当明显，尤其是次高自旋态。NO的非限制性计算，NO含有一个未成对电子（双线态基态2π），当用小基组(STO-3G或3-21G)计算时就呈现出严重的自旋污染。1.2半经验方法半经验方法是对从头算方法的简化，如通过限制分子轨道的选择或仅考虑价电子等再同实验数据进行拟合（如对有机分子的结构和构象能量拟合），由此可以大大加快计算的速度。扩展Hückel计算，它忽略了电子之间的相互作用，因而计算速度很快，但精度较差。该方法提供了分子轨道的空间形状和相对能量的定性估计以及电子密度的空间分布近似。该方法对化学结构的可视化以及化学反应的前线轨道处理非常有用。忽略微分重叠(Neglectofdifferentialoverlap，NDO)方法，NDO模型忽略了一些但不是全部电子-电子之间的相互作用。Zerner’sINDO方法

(ZINDO/1andZINDO/S)，MichaelZerner’s(Florida大学)开发的INDO版本，主要用于含过渡金属的分子体系，ZINDO/1用于预测分子几何而ZINDO/S则用于预测UV光谱。忽略双原子微分重叠(Neglectofdiatomicdifferentialoverlap,NDDO)方法，该方法基于INDO方法，提供了对完整HF矩阵的最相近的处理。它仅对属于不同原子的原子轨道使用零微分重叠近似，而保留全部双中心排斥积分。由于所保留的积分数目太多，因而难以过渡到半经验水平上去。对某些有机小分子和无机络合物的非经验NDDO计算，结果并不令人满意。改进的忽略双原子微分重叠(modifiedneglectofdiatomicdifferentialoverlay;MNDO，MNDDO)。它用半经验模型处理双中心排斥积分，其主要特点是把双原子的电荷分布之间的相互作用看成是两堆电荷分布的多极距相互作用之和，然后用实验数据优化参数。其参数仅包含单中心单电子能量和轨道指数，不仅所用参数比MINDO/3少，计算结果也比MINDO/3大有改进，而且机时只增加20％。但它对原子空间排列比较紧密地分子，如四员环、氢键、超键化合物、硝基和过氧化物等的预测则较差。一般地MNDO会过高地估计化学反应的活化位垒。AM1方法(AustinMethod,version1)，对MNDO方法的参数重新处理，包括核排斥项的改变。尽管AM1方法比MNDO方法精度要高，但它对磷氧键、硝基化合物以及过氧化物的键的处理不好。PM3、PM5、PM6方法(ParameterisationModel,version3、5和6)，另一种对MNDO方法的参数重新处理方法，功能类似于AM1，但比AM1方法有明显的改进。至今还没有发现其缺点。DFTB方法（DFT+TB），是将密度泛函方法与紧束缚方法结合起来，具有处理含有化学反应的较大体系的计算能力。典型的分子力学方法基本假设每个原子都可以用一个具有特定质量的粒子来代表；化学键可用弦来表示，用以表征两个原子间相互作用势能的力常数，势能函数可以描述分子内键的弯曲、拉伸和扭转以及分子间的一些现象，如经典相互作用或范德华(vanderWaals，VDW)力；势能函数依赖经验方法如从实验或从其它计算方法中得到。现行的分子力学方法是用一组势能函数来表征，用以描述化学力—力场，这些力场依赖于：原子位移（即键长）、原子类型和一些与原子类型及键相关的其它参数。AMBER(AssistedModelBuildingwithEnergyRefinement)力场,最初用于研究生物分子如蛋白质和核苷；CHARMM(ChemistryatHARvardMolecularMechanics)力场，最初也是用于生物分子和药物研究，但也被用来研究胶束和自组装大分子；MMx(MM2,MM3,等)

.用以研究非极性小分子的结构优化和热力学研究MMx力场包含对标准的键和键角势能面的二次拟合进行三级和四级校正，因此它可以处理分子振动中的非谐效应。不同的MMx版本其参数也不同，版本越高的MMx其较前一版本而言针对性就越强但适用性就越窄。例如较新的MM4版本就不能够适用于所有种类的分子；OPLS(OptimisedPotentialsforLiquidSimulations)力场，用于再现生物分子在溶液中的物理性质。因为这些力场方法都是针对生物分子的，这里我们不作过多的介绍。但这些方法对含能分子几何快速的初步优化还是有用的。1.4分子动力学方法(MD)分子动力学方法是假设一个体系有N个原子，而体系的状态由这N个原子的位置{ri}和动量{pi}或速度{vi}来标志，体系的能量为H({ri,pi})，体系的运动方程为：

分子动力学的主要目的是解上面的方程求得体系状态相空间演化的轨迹{ri,pi}t0，{ri,pi}t1，{ri,pi}t2，{ri,pi}t3，……。进而可计算我们感兴趣的物理量的值Q({ri,pi})。在实际应用中，我们常将上面的哈密顿方程化为下面的牛顿方程，并且用位置ri和速度vi做为描述体系的参量：；式中，V({ri})是原子间相互作用势，通过解上面的牛顿方程我们可以得到体系在相空间的轨迹，进而求得物理量的平均值[t(1),t(2),t(3),…t(M)]：，（m=1，2，……M）解牛顿方程的方法有很多种，通常有Verlet法则、Verlet速度法则、蛙跳法则、Gear提出的校正预测法等，最常用的为Verlet速度法则。分子动力学方法中最关键的是粒子间相互作用。系综选取根据研究的体系和方法还可将MD方法分成许多系综：微正则系综（NVT、NVE）等压系综（NPT、NPH）。分子动力学计算一般采用周期边界方法来满足实际模拟的需要，计算系统中分子间作用力时，采取最近镜象法，如下图所示：由于使用了最近镜象法，在计算时就必须采用截断半径(cut-offradius)方法计算远程的作用力，否则就会因重复计算粒子所受的力而导致不正确的结果。在MD模拟中截断半径不能大于周期盒边长的一半。粒子的最近镜像MD的计算步长选取计算步长的选取，在MD模拟中最重要的工作就是选取合适的步长。选取的标准是既能加快模拟速度又不能失去精度，通常的原则是取小于系统最快运动周期的1/10。例如水分子，分子内的运动为键长和键角的变化，分子间的运动质心的运动和分子的转动。分子间的运动源自范德华作用，一般较慢，而分子内运动则较快。由红外谱可知它的最大振动频率为1.08×1014sec-1，则振动周期为T=1/ν=0.92×10-14sec，那么计算时步长最大可取为约9×10-16sec下表给出了各种模拟体系的计算步长建议：体系运动状态步长（sec）简单原子体系移动10-14刚性分子移动、转动5×10-15软性分子、限制键长移动、转动、扭动2×10-15软性分子、无键长限制移动、转动、扭动、振动10-15或5×10-16MD方法的注意点经典分子动力学方法不能对含有化学反应的体系进行模拟，最近发展的反应力场方法(ReaxFF，A.vanDuin,W.A.Goddard,III)和从头算的CPMD（CarandParrinelloMD）以及半经验的MD方法可以对含化学反应的体系进行模拟。ReaxFF方法ReaxFF方法是嫁接量子力学(QM)和经典力场(FF)方法的桥梁，它属于经验方法而更接近于FF方法而不是QM方法。ReaxFF反应力场模型可以比QM方法花费低得多的计算机能力来模拟化学反应，其与传统FF方法相比主要的修改在于：引入了键级/键距方法，给出了共价键级连续解离的描述；加入了键级校正方案，可以允许键级方法区分自由基（离域）和基态（定域）键；在所有的键相互作用中都引入了键级（包括基于价键角和二面角相互作用相）；在所有原子间都进行非键相互作用（vanderWaals和Coulomb）计算，因此摒弃了不相容原理；加入了可极化电荷模型。AbinitioMD从头算分子动力学（abinitioMD）的计算方法同传统的MD模拟一样，都是基于牛顿定律，只不过abinitioMD方法通过解Schrödinger方程得到原子的力常数，但abinitioMD方法可以同时得到原子的电子结构，正因为如此，abinitioMD方法非常耗机时。——BOMD、CPMD半经验分子动力学的计算方法同从头算MD模拟一样，只不过半经验

MD方法中部分解来自于实验或其它低水平方法，由于半经验MD方法也可以同时得到原子的电子结构，而且需要的资源也远少于从头算MD方法，因此在含有化学反应的爆轰模拟上有独到的优势。——PM6/MD、DFTB/MD2含能材料的引发反应Coffey:含能材料的引发反应来自于含能材料固体内位错的贯通，当冲击波速度足够高时，位错产生的能量会直接传输给含能材料分子的内在振动模式振动能up-pumping模型指出：冲击波会产生构成晶体的分子的最低振动模式所吸收的激发声子热浴，增加的声子吸收以及分子间振动能的重新分配会进一步激发更高的振动模式，最终达到平衡过渡态，从而产生化学键的断裂并维持化学反应。值得注意的是，由于爆轰波前缘的传播速度对化学反应而言太快以至于不能驱动它们热激活。Gilman指出，在前缘的压缩会导致未反应物质的局部金属化，这种金属化由共价键的弯曲或由密度达到一定时，它们达到临界点后形成。在共价性材料中，当键被拉伸和弯曲时，就会有一个核组态空间区域，这时的最高占据轨道（HOMO）和最低未占轨道（LUMO）的能量差值就会减小。——固体含能材料会导致金属化因为爆轰波的前沿会引起未反应物质的压缩和剪切，特别是对共价性炸药材料而言，键的弯曲相对而言比键的压缩显得更为重要，这是因为与之相关的力常数会变小。

剪切导致金属化这一观点的主要证据是：通过对一系列不同的氨基-硝基-苯炸药的HUMO-LUMO差值的研究，其“shake-up”谱与其冲击感度有直接的关联在体积不变的情况下，键的弯曲会引起HUMO-LUMO值变小从而引起在键合态的电子的离域化，这样会引起快速的化学反应。对TATB而言，在体积接近不变时，临界剪切（弯曲）角与能隙的消失相对应。对非完美含能分子晶体进行的模拟研究，如RDX、PETN等含有孔洞和位错的情况，都表明当RDX晶体中存在单及双孔洞时，对其进行压缩它们的光学带隙会降低，导致它们的绝缘体—金属相转变所需的激发能下降。从直觉上讲，非均匀分子固体在猛烈撞击下会引起原子和电子性质的剧烈变化。一个可能的机理就是电子的激发：电子是属于轻、快和量子化的粒子，因此它首先对外部的微扰做出响应。接下来就会电子激发，几步就会引起化学反应的发生。这包括激发态表面的离解、辐射过程以及非辐射能的钝化。零温下，能带闭合也不能爆炸！！

TATB中硝基基团的弯曲振动受剪切的TATB在Fermi面附近的能带结构β-HMX的点电荷静电势3静高压下含能材料的行为现在认为TATB晶体和硝基甲烷中硝基的弯曲会导致能带的闭合，考察分子内硝基的弯曲在现在实验所能达到的压力范围内能带是否能够闭合可以通过静高压、均匀压力和非轴向加压下对硝基甲烷的电子结构计算，得到其光学带隙来确定。电子结构的计算可以采用自洽电荷密度泛函—紧束缚理论（SCC—DFTB）得到。这是基于密度泛函理论下对标准紧束缚模型的拓展。自洽部分描述的是总能、原子力和电荷的输运。实际应用中对于每个固定的晶胞参数集，所有的原子位置均由能量最小化确定，即原子是全弛豫的。SCC—DFTB的紧束缚特点是可以处理大的超晶胞，可包含更多个分子。早期的第一性原理计算则只能计算较小的超晶胞和较低的精度，但都给出了类似的带隙结果。以硝基甲烷为例，固定比值连续降低初始晶胞参数，压力P利用低温公式来估计。对静水压而言，当体积为原来的50%时，压力为约50GPa，而HOMO-LUMO减小了约0.6eV，即为原来值的13%。HOMO和LUMO同时增加而其差值却几乎随体积的变化单调下降（参见下左图）。在非常高的静水压下带隙的变化值示下面右图中。压缩率最高达70%。HOMO-LUMO值的减小，从理论上可以认为感度在增加。完美硝基甲烷晶体的原胞的HOMO-LUMO在静水压下硝基甲烷HOMO-LUMO值在静水压和非静水压下随比容的变化的变化值，突然向下的拐点对应高的C-H键弯曲C-H键形变的机理如下：随着初始C-N共价单键的收缩，电荷被传输到C原子和N原子之间的区域，随后C-N键变得更强，使得该键变为C=N共价双键。比邻的C-H键作为电子的供源，至少这些H原子中有一个的部分电子被拉走，导致该键变弱。另一方面，局域在硝基附近的轨道并没有被压缩的明显变形。由Mulliken电荷布居可以看出，H原子总的Mulliken电荷丢失随每个分子的形变而单调增加。这些发现意味着化学反应和共价键相离子键的转化会导致硝基甲烷的爆轰会始于带隙的闭合之前。最近的计算研究预测TATB金属化的压力在120GPa，远高于爆轰压力。对于硝基甲烷的情况，引起爆轰的原因应是质子的离解。这就有点意思了，在研究含能材料时，应该更注重它在单轴压缩而不是静水压缩下的化学行为，这与Dick对硝基甲烷的研究结果相一致。V/V0C1H1aH1bH1cC2H2aH2bH2c0.790.012-0.0100.0160.0160.012-0.0100.0160.0200.600.126-0.016-0.026-0.0260.108-0.019-0.008-0.1550.510.2940.014-0.074-0.0740.168-0.051-0.030-0.245硝基甲烷分子1和2中的C原子和H原子的Mulliken电荷（在单晶中沿b方向压缩） 4硝胺炸药的降解动力学模拟对凝聚炸药在高密度和高温下的化学反应机理的详细描述对于理解这些材料在燃烧和爆轰条件下反应面上所发生的事件有重要意义。如在冲击条件下，含能材料会被快速升温到几千度并被压缩到几百个千巴，体积的减少可达30%。这时会引发复杂的化学反应，释放出的大量能量会导致反应持续下去。很明显弄清极端条件下各种化学反应机理对于构建宏观水平上连续爆轰模型有重要意义。静态模拟结果基于自洽密度泛函紧束缚算法（SCC-DFTB）产物H2O、N2、CO2和CO的浓度C(t)随时间t的演化H2O、N2、CO2和CO的有效反应生成速率分别为：0.48、0.08、0.05和0.11ps-1（a）（b）（c）（d）β-HMX撞击加载下的初始分解动力学图1一个分子分解成C3H6N6O6

和CH2N2O2

两部分，另一个分子分解成C3H6N8O6；CH2和氧原子；图2显示C3H6N6O6

进一步分解为C2N4O4H4

，然后又和邻近的HMX生成的片断CN2O2H2

反应生成C3H6N6O6；图3显示生成物中有C2N4O4H4、CH2、

N2O、NO和CH2N2O2

等产物，图4显示生成物中有C3H6N6O6、C2N3O2H4、C2N2O2H4、NO2和N2O2；图5显示生成物中有C3H5O4H6、N2O2、CH2N2O2、CH2、N2O2和NO；图6显示到了反应后期生成物基本是一些小分子产物，而且里边还有一些游离的O原子、H原子和N原子等，这些原子很容易在当系统温度冷却下来的时候，生成H2,O2,N2,及H2O等小分子量产物。我们在此模拟中尚未发现有游离碳存在，这可能与我们模拟的时间较短有关，这个问题将在后续的工作中进一步研究。问题:

冲击压缩模拟在以下方面非常必要:诠释高温高压实验、炸药的冲击感度沿着冲击路径研究动力学和亚稳态解决办法:新的冲击模拟技术(MSST)使得在计算速度上大大加快原子尺度下的量子力学模拟为这种研究提供了工具得到化学反应的直观图像冲击条件下的从头算模拟计算原理基于Navier-Stokes方程的NEMD方法结合冲击动力学的三个守恒条件，在冲击条件下粒子的速度(u)与冲击波速度(vs)之间的关系可以写作：

u=vs(1-ρ0/ρ)，（1）Rayleigh线为冲击波方向上连接初态与终态（Hugoniot）的热力学路径：

p-p0=vs2ρ0(1-ρ0/ρ)，（2）在冲击条件下遵循能量守恒：e-e0=p0(1/ρ0-1/ρ)+vs2(1-ρ0/ρ)/2，（3）对分子动力学模拟而言，我们采用拉格朗日量：（4）这里T和V分别代表单位质量下的动能和势能,Q为模拟晶胞尺寸下的质量参数。我们可以发现当粒子速度的一阶微分为零时从哈密顿量的形式上方程(3)和(4)是等价的，这是因为T+V=e。由（5）可以看出，当粒子速度的二阶微分为零时，方程(5)就约化为方程(2)。因此通过方程(5)，我们就可以将冲击Hugoniot关系结合到非平衡分子动力学中。冲击模拟可以诠释与化学反应相关的热力学途径只能模拟静态的热力学单态整个依赖于热力学路径的化学过程可以得到反应只与冲击波速度和初始条件有关通过一次冲击模拟可以获得大量的热力学态静态计算vs.冲击模拟仅为给定温度和压力下的平衡态根据产物得到数据冲击波速度在9km/s时的初始分解与产物演化反应始于C-H键及N-O键的首先断裂；随着反应时间的进行，这些碎片会进一步分解成一个N2O2、C2N2O2H4、NO2

和C2N3O2H4；进而CH2和NO2

从各个碎片中分离出来，分别形成CN2O2H2、C2NH4、CNH2

、N2O、N2O2和CN4O4H2

显示生成物中有N2、O2、NO2、N2O和游离的O、H、CH、N、C等这些离子在最终会形成H2O、CO、CO2和凝聚碳等。5CO在高温高压下的歧化反应CO气体在自然界大量存在，也是炸药爆轰的主要产物之一。CO作为炸药的降解产物之一，并且一经产生便存在于高温高压的环境之中，在这样的极端条件下很容易形成更复杂的化合物，因此能否准确描述在特定条件下产物的化学组分以及化学反应过程对于理解炸药爆轰过程以及爆轰产物状态方程具有很重要的意义。CO在冲击状态下的聚合分解反应CO在冲击状态下的聚合分解反应t=0.16pst=0.18ps

气体CO在5000K,25GPa条件下的一些产物6极端条件下含能固体HMX结构

与性能的理论研究研究背景：感度分子结构晶体结构晶格聚集态缺陷空穴起爆机理起爆过程效率设计HMX研究背景：Bowden[1]和Kanllet[2]等人很早就对感度与分子结构之间的关系做了大量的工作，认为认识含能材料的感度不仅要从分子内基团的特征和数目，还必须要深入到分子的电子层结构层次，才能触及问题的本质。一般认为C-NO2键的断裂常常被认为是硝基化合物分子分解过程中非常关键的一步，但有的硝基化合物分子的分解是首先发生消环反应（如黑索金(RDX)分子）。因而，确定与C-NO2键的强度相关的因素就具有十分重要的意义。-NO2："致敏"-NH3：

"致钝"[1]BowdenF.P.,YoffeA.D.,InitiationandGrowthofExplosionsinLiquidsandSolids,CambridgeUniversityPress,Cambridge,1952.[2]KamletM.J.,Proceedings6thSymposium(International)onDetonation,SanDiego,California,ADAO59120,1976.

但是仅从分子角度出发得出结论往往并不完全正确，因为不同的晶体构型其感度相差很大，比如HMX炸药有四种不同的晶体结构，其中以β-HMX的感度最低。因此我们要寻求感度和其他因素间的关系。研究背景：研究背景：

HMX即奥克托金或环四次甲基四硝基胺(octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine)首先是由W.E.Bachmann和J.Z.Sheehan于20世纪30年代合成，是生产RDX的副产物，但比RDX有着更高的密度和熔点，自其合成以来以其最优的综合性能、极为广泛的使用率成为最炙手可热的猛炸药，许多含能材料的综合性能均将其作为参照比较的标准。

-HMX

的晶体及椅子型分子构型研究背景：包括2个C4H8N8O8分子，以b为对称轴，属于单斜晶系的P21/c准椅状的C4H8O8N8－HMX的晶体结构包括8个C4H8N8O8分子，属于斜方晶系的Fdd2空间群准船状的C4H8O8N8研究背景：－HMX的晶体结构包括6个C4H8N8O8分子，以c

轴为对称轴，属于P6122空间结构研究背景：－HMX的晶体结构包括两个C4H8N8O8·0.5H2O分子，具有Pn空间对称群H2O研究背景：研究背景：在实验方面：

由于HMX的特殊性，实验条件上受到很大限制。多集中在对-HMX的高压压砧实验，温度和压强对-HMX的光谱振动的影响；还有大量的实验是对-HMX研究样品的宏观性质的研究。在理论方面：

理论方面的研究多集中对在-HMX在一定条件下的分子间及分子内力场、振动、动能模型的研究。而对HMX四种构相的系统研究，尤其是其在极端条件下的系统研究就更少。研究内容及方法：

在微观尺度上对四种结构的HMX在极端条件（高温高压）下结构的、电子的、光学的及热力学性质进行系统研究，对三种纯HMX晶体在高温高压下的相变、力学性质以及常温常压下-HMX基PBXs在不同晶面上的力学性质进行探讨，尝试在介观尺度上对添加粘结剂条件下-HMX的相关性质进行讨论。

研究内容：研究内容及方法：研究方法：密度泛函理论经典分子力学方法耗散粒子动力学方法使用的软件：MaterialStudio采用了基于密度泛函理论的平面波赝势法，使用MS中的Castep模块；电子之间的交换相互作用尝试使用了广义梯度近似（GGA）中的PBE和局域密度近似（LDA）中的CA交换关联势；离子实和价电子之间的相互作用采用的是Vanderbilt超软赝势；布里渊区积分采用Monkhorst-pack方法，几何优化中采用Broyden-Fletcher-Goldfarb-Shanno(BFGS)

的算法，收敛标准为总能量偏差小于1×10-6eV/atom，原子的最大Hellmann-Feynman力小于0.05eV/Å，最大应力偏差小于0.01GPa，最大位移偏差小于0.002Å。结果与讨论：对-HMX采用4×2×3的K点对布里渊区进行采样（网格原点位于Γ点），平面波基函数的截断能取500eV；对-HMX采用2×1×4的K点，500eV的截断能；对-HMX采用3×3×1的K点，500eV的截断能；对-HMX采用1×2×1的K点，380eV的截断能。在进行深入的模拟之前，对同一体系采用更高的截断能和更密的K点采样进行试验，结果表明，这些改变使总能量的变化总小于0.002eV/atom。因此，我们所选取的参数对于对HMX来说是足够的。结果与讨论：结果与讨论：β-HMX的体积与压强的关系整个范围内，用GGA方法得到的数据比实验值偏高，而用LDA方法得到的数据比实验值偏低，相当而言LDA方法的结果要好一点。低压下的差别都很大。由于压强的增加，分子间距不断减小，分子间相互作用力不断增强，而使得GGA方法能更好的适合此时的体系，在15到25GPa的范围内GGA的方法比LDA方法的结果好。27GPa时发生相变，我们的结果与实验结果不符。结果与讨论：各个方向晶格参数的变化率和压强间的关系，b的压缩性最大，c的压缩性最小。表明了在各个方向的可压缩性有很大的差异。β-HMX各项异性结果与讨论：结果与讨论：N-N键的压缩率最大N-O键的压缩率最小结果与讨论：N7-O3-O4几乎不随压强的变化而变化H13-C2-N4的变化最大结果与讨论：H14-C2-N5-C1的变化最大，17.5%结果与讨论：β-HMX整个布里渊区价带和传导带的分布都很平坦，这是由-HMX的分子特点决定的。价带和导带之间的电子带隙为3.061eV标志着β-HMX是一绝缘体，在40GPa时带隙值为2.436eV。在0GPa时为传导带的G点和价带的B点的间接带隙，而压强在20GPa时为G-G的直接带隙。随着压强的增加，压强对带宽的增加的影响越来越小。

结果与讨论：随压强的升高带隙在不断的减小，使得β-HMX变得越来越不稳定，传导性在不断加强，表明通过加压可能使其成为半导体或导体。结果与讨论：0GPa时的部分态密度。可以判断β-HMX的基本性质：在高价带，C原子的态密度值比H、O和N原子的态密度值要低很多，表明了H、O和N原子在高价带比C原子的贡献多很多。在-25到-13eV的能量范围内,H、O和N原子存在着明显的峰值。在传导带区域，峰值主要源于C、H和N原子的p态。

结果与讨论：和0GPa时相比，25GPa时的态密度范围明显加宽，而且存在很大范围的负频，表明β-HMX晶体的不稳定性随压强的升高而增强。说明β-HMX随压强的增加越来越敏感。

结果与讨论：（a）通过德拜模型方法（b）通过声子谱方法。所有的焓、熵和温度乘积、随温度的升高而不断增大，自由能随温度升高而不断减小。结果与讨论：通过两种方法得到的自由能、焓、熵和温度乘积、热容随温度的变化趋势类似，但数值上相差很大。在文献[1]中的实验值是在从310到440K温度范围内β-HMX的热容变化为499.43J/mol·K增大到588.37J/mol·K，在相应的温度范围内，用声子谱法得到的热容值从589.97J/mol·K增加到775.93J/mol·K，而用德拜模型方法得出的热容值从1363.59J/mol·K增大到1380.35J/mol·K。可以判断通过声子谱方法得到的相应的数据更可靠。德拜模型的准谐近似对低对称性的有机分子晶体可能不太适合。结果与讨论：β-HMX的德拜温度(a)在0GPa时通过声子谱方法(b)在0GPa时通过德拜模型方法(c)在0K时通过德拜模型方法。变化趋势的差距很大，也进一步表现出两种方法的差异。结果与讨论：在不同压强下的吸收系数曲线的形状类似但对应的能量值不同，且随着压强的增大发生了平移。吸收系数的数值量级很大(1014nm-1)，表明β-HMX在光的作用下会很不稳定且随着压强的增加不稳定性不断增强。结果与讨论：结果与讨论：-HMX在0和30GPa时的态密度结果与讨论：-HMX在0和30GPa时的态密度结果与讨论：δ-HMX在0和30GPa时的态密度 对-、-、-HMX的态密度在费米面附近都有峰值，表明了其对应的高价带的带结构比较平缓。三种结构HMX中的C原子的态密度有很大不同，是由它们各自的分子对称性不同决定的；而苯环上的N和NO2上的N的态密度的变化却不大。在传导带范围内，峰值主要由N原子的p态决定；在能量较高的位置范围内，主要由C和N原子态密度决定；在费米面附近的总的态密度的峰值随压强的增大而增多。结果与讨论：而-HMX的吸收系数的峰值随压强的升高而减小，且-HMX的吸收系数在数值上比-、-HMX的小的多，-、-HMX的吸收系数峰值为105数量级，而-HMX的仅为104数量级，说明在高压下，-HMX对光的敏感度降低，-HMX比β-、-、-HMX更稳定些。-、-HMX的吸收系数的峰值都随压强的升高而增大，对应的能量值在减小。结果与讨论：小结：在低压下N-N键的变化对-HMX稳定性的影响最大，而高压下N-O键及O7-N3-N4键角对-HMX稳定性的影响很大。通过准谐德拜模型和声子谱两种方法研究了-HMX晶体的热力学性质表明用声子谱方法得出的结果比较合理，而准谐德拜模型方法对低对称性的-HMX的计算不够准确，原因可能是德拜模型采用的近似方法对-HMX这样的低对称性的晶体不合适。在0-40GPa内并未观察到HMX的相变，可能是我们的模拟体系不足以在空间群上反映其相变。-HMX对光的敏感程度最小，且随着压强的增大光敏感的在不断下降，而-、-HMX的光感度随压强的增大而增强，可以推断在高压下-HMX比其他几种结构相晶体都稳定。

在MS程序中的Discover模块下的恒温恒压（NPT）的分子动力学方法模拟-HMX、-HMX和-HMX在不同条件下进行研究，采用COMPASS(condensed-phaseoptimizedmolecularpotentialsforatomisticsimulationstudies)

力场，采用NPT系综中的Andersen的调温方法及Berendsen调压方法。

三种结构的HMX的初始结构如前。模拟的超晶胞包括96个HMX分子，2688个原子，对-HMX、-HMX和-HMX相应的盒子大小分别为由4×4×3、2×2×3、4×2×2的超晶胞组成，保证分子间作用势的截断半径在10.5Å。结果与讨论：HMX的温度和压力相变从左到右依次为-HMX、-HMX和-HMX晶体的超晶胞结构示意图

模拟条件为：-HMX在298K时0到31GPa及0GPa时温度从5到400K；-HMX在298K时压强从22到40GPa及0GPa时温度从350到480K；-HMX在298K时压强从20到40GPa及0GPa时温度从430到555K。所有模拟均采用2×105时间步，每个时间步为1fs，其中前1×105时间步用来平衡结构，后1×105时间步用来收集结果。在平衡结构时，为了保证HMX晶体的温度和压强能与模拟外部的温度及压强保持一致，每间隔1000时间步标定一次。相互作用是由原子核在截断范围内的最邻近原子决定的。所有的计算考虑周期性边界条件，非键相互作用分别采用Ewald和基于原子的方法来控制的库仑力和长程范德瓦尔斯力。结果与讨论：a/Åb/Åc/Åα(º)β(º)(º)V(Å3)(g/cm3)β-HMXExpt.[20]6.5411.058.7090124.390259.69351.89369Thiswork6.5811.128.7690124.390264.91351.85638α-HMXExpt.[5]15.1423.895.91909090267.20251.84048Thiswork15.3124.175.98909090276.60751.77790δ-HMXExpt.[8]7.717.7132.559090