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

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,*,单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,*,单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,*,单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,*,单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,*,2024/8/23,材料的计算,机模拟方法,2023/8/30材料的计算机模拟方法,内容,材料计算方法以及最新的进展；,含能材料的引发反应；,静高压下材料的行为；,材料的降解动力学模拟；,气态CO在高温高压下的歧化反应；,HMX,固体的性能理论模拟,内容材料计算方法以及最新的进展；,1,计算方法简介,在含能材料的模拟计算中，目前主要有三种方法：,第一性原理（或从头算）的电子结构计算,(,ab initio,),、,半经验方法,分子力学方法,(MM),。,1 计算方法简介在含能材料的模拟计算中，目前主要有三种方法,方法比较,方法,优点,缺点,最适宜的领域,Ab initio,采用量子物理方法,数学上严格：没有经验参数,半经验方法,采用量子物理方法,使用实验参数,使用更多的近似,分子力学方法,采用经典物理学方法,需要有嵌入经验参数的力场,适用的体系广泛,不需要实验数据,可以计算激发态和过渡态,与从头算方法相比对计算机的要求较低,可计算过渡态和激发态,在有限的计算资源下计算,“,便宜,”,，快且有用,可以用于大分子如酶的计算中,非常耗费机时,需要有从头算方法或实验给出参数,精度比从头算方法低,不能计算电子的性质,需要从头算或实验数据作参数,商业软件只能用于有限范围的分子,小体系（数十个原子,几百个原子）,有电子转移的体系,不需要实验数据的体系,需要高精度的体系,中等尺寸的体系,(,几百个原子,数千个原子,),有电子转移的体系,大体系（数千,数亿个原子）,体系中不能包含有键断裂的情况,方法比较方法优点缺点最适宜的领域Ab initio适用的体系,通过对材料的模拟计算可以得到如下的性质：,生成热、,键能和反应能、,分子（或晶体）的能量和结构热化学稳定性,),、,过渡态的能量和结构（活化能）、,反应路径以及动力学、,机理、,分子中的电荷分布（反应点）、,取代效应、,振动频率（红外和,Raman,光谱）、,电子的激发跃迁（,UV/,可见光谱）、,磁屏蔽效应（,NMR,谱）,。,通过,对这些性质的预测，在含能材料中有许多实用价值，如：可以研究合成路径、反应产物、爆轰产物状态方程、爆轰引发的机理等等。,通过对材料的模拟计算可以得到如下的性质：,1.1,从头算,(,ab initio,),方法,从头算方法在计算中是基于量子力学基本定律，没有引入经验或实验参数，仅使用,5,个常数（电子质量： 9.10938215E,-31,Kg、电子电量：1.602176487E,-19,C、光速：,c,= 2.99,792458,10,8,m/s,、,Planck,常数：,h,= 6.6260689610,-34,Js,和波尔兹曼常数1.3806504E,-23,J/K），对化学体系利用满足严格数学近似的基组来解,Schrdinger,方程：,动能是所有粒子的动能加和,势能代表粒子之间的库仑相互作用,（即：核与核之间、电子与电子之间,的排斥以及电子与核之间的吸引）,对电子而言，,q,=-,e,，对核而言，原子序数为,Z,：,q,=+,Ze,。,1.1 从头算(ab initio) 方法从头算方法在计算中,通过对,Schrdinger,方程的形式观察，我们知道对于多于两个以上粒子的分子，严格解,Schrdinger,方程是不可能的，必须引入一些必要的合理近似，甚至一些经验参数，如非相对论解、,Born-Oppenheimer,近似、单粒子近似、原子轨道的线性组合等。对解出的,Schrdinger,方程是否合适，必须满足以下判据：,这个解必须很好地被定义且能够用来说明分子结构和分子中电子的状态。,分子的势能按照原子核的位移必须平滑连续可变。,模型不能有偏离（即：假设化学键只存在于原子之间）。,模型必须与体系大小相关（即解和相关误差必须与分子的大小成比例），,模型必须是可变的（即：近似解必须提供体系能量真值的上界，因此近似解应具有最低能量值与真实的波函数相一致）。,通过对Schrdinger方程的形式观察，我们知道对于多于,1.1.1,基组,一般来说，基组属于数学函数的一类，用来解微分方程。,基组：,G,、,PW,、,G+PW,量化计算，基组一般约定为高斯基组来代表原子轨道，用来优化和产生所需体系的化学性质,物理计算（,DFT,方法）,特别是金属的计算，经常用,PW,；,GPW,标准的从头算软件包都提供基组、当分子中包含较大的原子（,Z30,）时，我们经常将原子的内层电子合并作近似处理（赝势）。,1.1.1 基组一般来说，基组属于数学函数的一类，用来解微分,极小基包含最少的原子轨道基函数，如对,H,和,He,仅包含,1s,电子、对,Li,到,Ne,只含,1s, 2s, 2px, 2py, 2pz,电子。极小基组的一个例子是,STO-3G,每个基函数采用,3,个,gaussian,型函数,(3G),用来表示,Slater,型原子轨道近似。尽管极小基不能精确预测分子能量，但是用它来做分子和化学键的可视化却非常好。,劈裂价基使用约定的高斯函数加上一些高斯基函数的线性组合来表示每个原子价轨道， 劈裂价基常用一些指定的价轨道函数来表征。,双,zeta,基就使用两个基函数来描述价电子、,三,zeta,基用三个基函数等等。,Pople,及其同事对基函数作出了很大贡献，他们用一组高斯函数来描述内层和外层电子。,如,6-21G,表示内层原子轨道用六个高斯基函数来描述，内层价壳则用两个高斯基函数表示，外层价壳用一个高斯函数表示。,类似地还有：,3-21G, 4-31G,以及,6-311G,等等。,极小基包含最少的原子轨道基函数，如对H和He仅包含1s电子、,带极化函数的基组，在基组中加入极化函数后可以允许原子核外电荷的不均匀位移，因此它可以改进对化学键的描述。极化函数可以描述那些和孤立原子相比更高的角动量量子数的轨道（如,p-,型函数用于,H,和,He,，,d-,型函数用于,Z2,的原子），并将之加入到价电子壳层中。例如，,6-31G(d),基组的构建是将,6,个,d,-,型简单高斯函数加入到,6-31G,中来描述每个非氢原子； 同样地，,6-31G(d,p),是,6-31G(d),对重原子,但加入了一组,p,-,型高斯函数到氢原子或氦原子。 将,p,-,轨道加入到氢原子中对那些含有以氢原子为桥原子的体系而言特别有用。,带弥散函数的基组，对于那些电子密度明显远离核的中心的原子或离子（如阴离子、孤电子对以及激发态），就需要用弥散函数来解决核对最外面电子束缚较弱的问题。一般来说，建议将弥散基组用于处理阴离子的电子亲和性、质子亲和性、反演势垒以及键角等。只对非氢原子加入弥散的,s,-,和,p,-,型高斯函数，弥散基组用加号“,+”,表示，如,3-21+G,；若对重原子和氢原子均加入弥散函数就用两个加号“,+”,表示。,带极化函数的基组，在基组中加入极化函数后可以允许原子核外电荷,高角动量基组（一般的,Hartree-Fock,（以后简称,HF,）计算中不需要）,带有多个极化函数的基组，已经被用于很多体系的计算中了。这种基组对于描述有电子之间相互作用的电子相关方法很有用。高角动量基组的例子有：,6-31G(2d) ,两个,d,-,函数加到重原子中；,6-311G(2df, pd) .,除了,311,价函数外，两个,d,函数和一个,f,函数加到重原子中，而对氢原子则考虑了,p,和,d,轨道函数；,6-311G(3df, 2df, p) ,三个,d,函数和一个,f,函数加到,Z, 11,的重原子中，两个,d,函数和一个,f,函数加到第一行原子,(Li,到,Ne),，而对氢原子则考虑了,p,轨道函数。,高角动量基组（一般的Hartree-Fock（以后简称HF）,1.1.2,电子相关方法,由于,HF,方法为单粒子近似，忽略了电子之间的相互作用，因此会过低的估计键长同时会过高地估计振动频率，因此发展出了三种考虑电子相关的方法（或称,post-HF,方法）：,组态相互作用（,configuration-interaction,，,CI,）,Mller-Plesset (MP),微扰理论,密度泛函理论（,density functional theory,，,DFT),1.1.2 电子相关方法由于HF方法为单粒子近似，忽略了电子,由于在分子的,HF,波函数里只有一个行列式，因此它在分子中只能描述一个单电子组态。这就限制了,HF,方法,它对最简单的分子,氢的二聚体的离解能都不可以精确预测。,组态相互作用方法是在原来的,HF,波函数的基础上再构建一个行列式包含更高能量（激发态）的空轨道，这个新的行列式通过用更高能量的未占轨道（虚轨道）取代一个或更多个已占轨道来建立，行列式中取代的数目代表了,CI,的水平，,单取代（,CIS,）打开一个占据轨道和一个虚轨道，它等价于单电子激发态。,更高级的计算包括双取代（,CID,）对应于打开两对轨道：,CISD,代表增加单取代和双取代的行列式，,CISD,（,T,）则代表单、双和三激发态。,理论极限为全组态相互作用，它是用,HF,行列式的线性组合表达可能取代的行列式。全,CI,方法提供了分子体系最完全的非相对论处理，当然它也非常耗机时。,组态相互作用（,CI,）：,由于在分子的HF波函数里只有一个行列式，因此它在分子中只能描,对典型的分子而言，大部分的基态能量都源自于单电子,HF,的能量。,MP,微扰理论假设电子相关的效应很小，可以用很小的矫正（即微扰）来得到,HF,的解。,MP,方法假设分子真正的,Hamiltonian,量可分成两部分：,（其中： 代表单电子的贡献； 代表组态相互作用的贡献）,则代表能量和分子波函数的展开系数，即： ； 代入分子的,Schrdinger,方程中就可以得到：,E,HF,=,E,(0),+,E,(1),+,E,(2),+,这就是所谓的,MP2,、,MP3,、,MP4,、,MP5,等等。,Mller-Plesset,微扰理论,对典型的分子而言，大部分的基态能量都源自于单电子HF的能量。,这种方法是将电子相互作用看作是电子密度的函数，现行的,DFT,方法通过,Kohn-Sham,方程将电子能量分成,4,个部分：,E,=,E,T,+,E,V,+,E,J,+,E,XC,，其中,E,T,代表电子的动能,(,来自于电子的运动,),；,E,V,代表电子的势能（包括核与核之间的排斥能、核与电子之间的相互作用能）；,E,J,代表电子之间的排斥能；,E,XC,代表电子相关能。,由于对,E,XC,的定义不同便产生了不同的,DFT,方法，如：,LDA,、,GGA,、,B3LYP,、,B3PW91,等等。,DFT,方法可分成三类：局域密度近似（,LDA,）、广义梯度近似（,GGA,）和,DFT,与,HF,项结合的杂化方法：,LDA,方法只含有与电子密度有关的项，可以用来预测块材料的性质，但对于孤立分子的性能预测效果就不行；,GGA,方法同时考虑了电子密度和电子密度梯度，用它预测的键离解能与实际仅相差几个,kJ/mol,。但,GGA,方法却严重地过低估计了一些反应的活化位垒，这是因为这种方法忽略了电子的库仑子相互作用的原因；,通过,DFT,与,HF,杂化的方法加以解决GGA方法产生的问题。,密度泛函理论（,DFT,）,这种方法是将电子相互作用看作是电子密度的函数，现行的DFT方,1.1.3,闭壳层和开壳层计算,大多数分子具有电子单线态基态（电子的总自旋数为零），这时分子内电子完全成对地在限制在轨道内，这时就要用闭壳层自旋限制计算。,对于那些具有非成对电子的体系，就要采用开壳层非自旋限制方法计算，计算中考虑电子的自旋状态（分为自旋向上,电子和自旋向下,电子）。,1.1.3 闭壳层和开壳层计算大多数分子具有电子单线态基态（,开壳层非限制性自旋计算适用的体系：,具有奇数个电子的分子,(,例如一些离子和自由基,),；,激发态；,分子有三线态或更高自旋的基态,(,例如,O,2,),；,电子对的分离过程,(,例如分解等,),；,共振体系，存在电子离域。,开壳层非限制性自旋计算适用的体系：,自旋污染,非限制性自旋波函数并不是自旋算符的本征函数，因此要对自旋态进行混合，从而导致自旋污染。,在大多数开壳层系统中，自旋污染是可以忽略的。然而对于未成对电子离域的分子，自旋污染效应相当明显，尤其是次高自旋态。,NO,的非限制性计算，,NO,含有一个未成对电子（双线态基态,2,），当用小基组,(STO-3G,或,3-21G),计算时就呈现出严重的自旋污染。,自旋污染非限制性自旋波函数并不是自旋算符的本征函数，因此要对,在含能材料的计算中，由于包含硝基,(-NO,2,),，在计算激发态或有自由基参与的反应里，也存在自旋污染现象。,自旋污染效应的存在会造成几何预测和振动频率的预测不准确，当然自旋污染效应会随基组的增大而降低。,当占据轨道和非占据轨道非常靠近时，自旋污染就特明显。如当键被拉伸到快要断裂时，这时成键和反键轨道的能量差接近于零，产生的波函数就包含等量的单线态和三线态。,Gaussian,软件包已经在其,UHF,计算里自动加入了自旋污染校正。当然密度泛函方法与,HF,方法相比有较低的自旋污染效应，计算时可以采用，但最好在存在自旋污染的体系里不要采用,DFT-HF,杂化方法。,在含能材料的计算中，由于包含硝基(-NO2)，在计算激发态或,1.2,半经验方法,半经验方法是对从头算方法的简化，如通过限制分子轨道的选择或仅考虑价电子等再同实验数据进行拟合（如对有机分子的结构和构象能量拟合），由此可以大大加快计算的速度。,1.2 半经验方法半经验方法是对从头算方法的简化，如通过限制,扩展,Hckel,计算，它忽略了电子之间的相互作用，因而计算速度很快，但精度较差。 该方法提供了分子轨道的空间形状和相对能量的定性估计以及电子密度的空间分布近似。该方法对化学结构的可视化以及化学反应的前线轨道处理非常有用。,忽略微分重叠,(Neglect of differential overlap,，,NDO),方法，,NDO,模型忽略了一些但不是全部电子,-,电子之间的相互作用。,扩展Hckel计算，它忽略了电子之间的相互作用，因而计算速,各种近似解,Schrdinger,方程的,HF,自洽场方法,(Hartree-Fock Self-Consistent Field,，,HF-SCF),全略微分重叠,NDO (Complete-NDO,，,CNDO),方法，不同原子的两个原子轨道之间的乘积都为零，这时不同轨道间电子的排斥仅依赖于所含原子的性质而与具体的轨道无关。因为,CNDO,几乎忽略了对所有电子交换性质的描述 ，它不能区分具有相同电子组态而自旋值不同的那些电子态。,非全略微分重叠,(Intermediate NDO,，,INDO),方法，在描述电子排斥方面。考虑了相同原子之间的轨道重叠微分，但忽略了不同原子之间的轨道重叠微分,改进的,INDO,方法,(,Modified INDO,版本,3,，,MINDO/3),，,INDO,的参数化版本，可以很好地预测一些化学体系的生成焓以及更合理的分子几何，特别是含硫的化合物、碳阳离子以及多硝基化合物。,各种近似解Schrdinger方程的HF自洽场方法(Har,Zerners INDO,方法,(ZINDO/1 and ZINDO/S),，,Michael Zerners (Florida,大学,),开发的,INDO,版本，主要用于含过渡金属的分子体系，,ZINDO/1,用于预测分子几何而,ZINDO/S,则用于预测,UV,光谱。,忽略双原子微分重叠,(Neglect of diatomic differential overlap, NDDO),方法，该方法基于,INDO,方法，提供了对完整,HF,矩阵的最相近的处理。它仅对属于不同原子的原子轨道使用零微分重叠近似，而保留全部双中心排斥积分。由于所保留的积分数目太多，因而难以过渡到半经验水平上去。对某些有机小分子和无机络合物的非经验,NDDO,计算，结果并不令人满意。,改进的忽略双原子微分重叠,( modified neglect of diatomic differential overlay;MNDO,，,MNDDO),。它用半经验模型处理双中心排斥积分，其主要特点是把双原子的电荷分布之间的相互作用看成是两堆电荷分布的多极距相互作用之和，然后用实验数据优化参数。其参数仅包含单中心单电子能量和轨道指数，不仅所用参数比,MINDO/3,少，计算结果也比,MINDO/3,大有改进，而且机时只增加,20,。但它对原子空间排列比较紧密地分子，如四员环、氢键、超键化合物、硝基和过氧化物等的预测则较差。一般地,MNDO,会过高地估计化学反应的活化位垒。,Zerners INDO方法 (ZINDO/1 and Z,AM1,方法,(,Austin Method, version 1,),，对,MNDO,方法的参数重新处理，包括核排斥项的改变。尽管,AM1,方法比,MNDO,方法精度要高，但它对磷氧键、硝基化合物以及过氧化物的键的处理不好。,PM3,、PM5、PM6方法,(Parameterisation Model, version3,、5和6,),，另一种对,MNDO,方法的参数重新处理方法，功能类似于,AM1,，但比,AM1,方法有明显的改进。至今还没有发现其缺点。,DFTB,方法（,DFT+TB,），是将密度泛函方法与紧束缚方法结合起来，具有处理含有化学反应的较大体系的计算能力。,AM1方法 (Austin Method, version,1.3,分子力学方法,(Molecular mechanics,，,MM),分子力学方法对于一些大体系以及非对称的化学体系如蛋白质或聚合物的预测非常有效的手段，,MM,方法是一个纯粹的经验方法，它仅依赖于经典物理定律而忽略了电子的处理，因此它不能预测分子的化学性质，结果是,MM,方法不能涉及键的断裂和形成等有电子或量子效应参与的体系，,MM,方法预测的能量是以无意义的绝对量，它只能作相对比较时才有意义。,但它确是连接量子力学和连续力学的桥梁，而具有特殊的应用价值。它已被广泛地应用在含能材料的细观效应研究里，包括反应模型和经典势能面的分裂、晶体的平衡性质（即：密度、填充、比热等）研究、晶体及其缺陷在动态加载下的相互作用以及分子晶体中模拟爆轰等。,用作分子的可视化,1.3 分子力学方法(Molecular mechanics,典型的分子力学方法基本假设,每个原子都可以用一个具有特定质量的粒子来代表；,化学键可用弦来表示，用以表征两个原子间相互作用势能的力常数，势能函数可以描述分子内键的弯曲、拉伸和扭转以及分子间的一些现象，如经典相互作用或范德华,(van der Waals,，,VDW),力；,势能函数依赖经验方法如从实验或从其它计算方法中得到。,现行的分子力学方法是用一组势能函数来表征，用以描述化学力,力场，这些力场依赖于：原子位移（即键长）、原子类型和一些与原子类型及键相关的其它参数。,典型的分子力学方法基本假设 每个原子都可以用一个具有特定质量,AMBER,(,Assisted Model Building with Energy Refinement,),力场,最初用于研究生物分子如蛋白质和核苷；,CHARMM,(,Chemistry at HARvard Molecular Mechanics,),力场，最初也是用于生物分子和药物研究，但也被用来研究胶束和自组装大分子；,MMx,(,MM2, MM3,等,),.,用以研究非极性小分子的结构优化和热力学研究,MMx,力场包含对标准的键和键角势能面的二次拟合进行三级和四级校正，因此它可以处理分子振动中的非谐效应。不同的,MMx,版本其参数也不同，版本越高的,MMx,其较前一版本而言针对性就越强但适用性就越窄。例如较新的,MM4,版本就不能够适用于所有种类的分子；,OPLS,(,Optimised Potentials for Liquid Simulations,),力场，用于再现生物分子在溶液中的物理性质。,因为这些力场方法都是针对生物分子的，这里我们不作过多的介绍。但这些方法对含能分子几何快速的初步优化还是有用的。,AMBER (Assisted Model Building,1.4,分子动力学方法,(MD),分子动力学方法是假设一个体系有,N,个原子，而体系的状态由这,N,个原子的位置,r,i,和动量,p,i,或速度,v,i,来标志，体系的能量为,H,(,r,i,p,i,),，体系的运动方程为：,分子动力学的主要目的是解上面的方程求得体系状态相空间演化的轨迹,r,i,p,i,t,0,，,r,i,p,i,t,1,，,r,i,p,i,t,2,，,r,i,p,i,t,3,，,。进而可计算我们感兴趣的物理量的值,Q(,r,i,p,i,),。在实际应用中，我们常将上面的哈密顿方程化为下面的牛顿方程，并且用位置,r,i,和速度,v,i,做为描述体系的参量：,； 式中，,V,(,r,i,),是原子间相互作用势，通过解上面的牛顿方程我们可以得到体系在相空间的轨迹，进而求得物理量的平均值,t(1),t(2),t(3),t(M),：,，（,m,=1,，,2,，,M,）,解牛顿方程的方法有很多种，通常有,Verlet,法则、,Verlet,速度法则、蛙跳法则、,Gear,提出的校正预测法等，最常用的为,Verlet,速度法则。分子动力学方法中最关键的是粒子间相互作用。,1.4 分子动力学方法(MD)分子动力学方法是假设一个体系有,常用的势函数和应用范围,两体势,(Lennard-Jones,势和,Morse,势,),， 适用于惰性气体、简单金属；,原子内嵌势,(Embed Atom Methods,，,EAM,势,),，适用于简单金属，过渡金属；,类原子内嵌势,(Johnson, Meipotential, Glue Potential, Finnis-Sutton,势,),，适用于简单金属，过渡金属；,紧束缚势,(Tight Binding),，适用于合金等；,Tersoff, Brennerd, Stillinger-Weber,势，适用于半导体；,从头算势函数，适用范围宽。,常用的势函数和应用范围 两体势(Lennard-Jones势,系综选取,根据研究的体系和方法还可将,MD,方法分成许多系综：,微正则系综（,NVT,、,NVE,）,等压系综（,NPT,、,NPH,）。,系综选取根据研究的体系和方法还可将MD方法分成许多系综：,分子动力学计算一般采用周期边界方法来满足实际模拟的需要，计算系统中分子间作用力时，采取最近镜象法，如下图所示：,由于使用了最近镜象法，在计算时就必须采用截断半径,(cut-off radius),方法计算远程的作用力，否则就会因重复计算粒子所受的力而导致不正确的结果。在,MD,模拟中截断半径不能大于周期盒边长的一半。,粒子的最近镜像,分子动力学计算一般采用周期边界方法来满足实际模拟的需要，计算,MD,的计算步长选取,计算步长的选取，在,MD,模拟中最重要的工作就是选取合适的步长。选取的标准是既能加快模拟速度又不能失去精度，通常的原则是取小于系统最快运动周期的,1/10,。例如水分子，分子内的运动为键长和键角的变化，分子间的运动质心的运动和分子的转动。分子间的运动源自范德华作用，一般较慢，而分子内运动则较快。由红外谱可知它的最大振动频率为,1.0810,14,sec,-1,，则振动周期为,T=1/=0.9210,-14,sec,，那么计算时步长最大可取为约,910,-16,sec,下表给出了各种模拟体系的计算步长建议：,体系,运动状态,步长（,sec,）,简单原子体系,移动,10,-14,刚性分子,移动、转动,510,-15,软性分子、限制键长,移动、转动、扭动,210,-15,软性分子、无键长限制,移动、转动、扭动、振动,10,-15,或,510,-16,MD 的计算步长选取计算步长的选取，在MD模拟中最重要的工作,MD,方法的注意点,经典分子动力学方法不能对含有化学反应的体系进行模拟，,最近发展的反应力场方法,(ReaxFF,，,A. van Duin, W. A. Goddard, III ),和从头算的,CPMD,（,Car and Parrinello MD,）以及半经验的MD方法可以对含化学反应的体系进行模拟。,MD方法的注意点经典分子动力学方法不能对含有化学反应的体系进,ReaxFF,方法,ReaxFF,方法是嫁接量子力学,(QM),和经典力场,(FF),方法的桥梁，它属于经验方法而更接近于,FF,方法而不是,QM,方法。,ReaxFF,反应力场模型可以比,QM,方法花费低得多的计算机能力来模拟化学反应，其与传统,FF,方法相比主要的修改在于：,引入了键级,/,键距方法，给出了共价键级连续解离的描述；,加入了键级校正方案，可以允许键级方法区分自由基（离域）和基态（定域）键；,在所有的键相互作用中都引入了键级（包括基于价键角和二面角相互作用相）；,在所有原子间都进行非键相互作用（,van der Waals,和,Coulomb,）计算，因此摒弃了不相容原理；,加入了可极化电荷模型。,ReaxFF方法 ReaxFF方法是嫁接量子力学(QM)和经,Ab initio MD,从头算分子动力学（,ab initio,MD,）的计算方法同传统的,MD,模拟一样，都是基于牛顿定律，只不过,ab initio,MD,方法通过解,Schrdinger,方程得到原子的力常数，但,ab initio,MD,方法可以同时得到原子的电子结构，正因为如此，,ab initio,MD,方法非常耗机时。BOMD、CPMD,半经验分子动力学的计算方法同从头算,MD,模拟一样，只不过半经验,MD,方法中部分解来自于实验或其它低水平方法，由于半经验,MD,方法也可以同时得到原子的电子结构，而且需要的资源也远少于从头算MD方法，因此在含有化学反应的爆轰模拟上有独到的优势。PM6/MD、DFTB/MD,Ab initio MD从头算分子动力学（ab initio,2,含能材料的引发反应,Coffey :,含能材料的引发反应来自于含能材料固体内位错的贯通，当冲击波速度足够高时，位错产生的能量会直接传输给含能材料分子的内在振动模式,振动能,up-pumping,模型指出：冲击波会产生构成晶体的分子的最低振动模式所吸收的激发声子热浴，增加的声子吸收以及分子间振动能的重新分配会进一步激发更高的振动模式，最终达到平衡过渡态，从而产生化学键的断裂并维持化学反应。,值得注意的是，由于爆轰波前缘的传播速度对化学反应而言太快以至于不能驱动它们热激活。,2 含能材料的引发反应Coffey :含能材料的引发反应来自,Gilman,指出，在前缘的压缩会导致未反应物质的局部金属化，这种金属化由共价键的弯曲或由密度达到一定时，它们达到临界点后形成。在共价性材料中，当键被拉伸和弯曲时，就会有一个核组态空间区域，这时的最高占据轨道（,HOMO,）和最低未占轨道（,LUMO,）的能量差值就会减小。,固体含能材料会导致金属化,因为爆轰波的前沿会引起未反应物质的压缩和剪切，特别是对共价性炸药材料而言，键的弯曲相对而言比键的压缩显得更为重要，这是因为与之相关的力常数会变小。,剪切导致金属化这一观点的主要证据是：通过对一系列不同的氨基,-,硝基,-,苯炸药的,HUMO-LUMO,差值的研究，其“,shake-up”,谱与其冲击感度有直接的关联,Gilman指出，在前缘的压缩会导致未反应物质的局部金属化，,在体积不变的情况下，键的弯曲会引起,HUMO-LUMO,值变小从而引起在键合态的电子的离域化，这样会引起快速的化学反应。,对,TATB,而言，在体积接近不变时，临界剪切（弯曲）角与能隙的消失相对应。,在体积不变的情况下，键的弯曲会引起HUMO-LUMO值变小从,对非完美含能分子晶体进行的模拟研究，如,RDX,、,PETN,等含有孔洞和位错的情况，都表明当,RDX,晶体中存在单及双孔洞时，对其进行压缩它们的光学带隙会降低，导致它们的绝缘体,金属相转变所需的激发能下降。,从直觉上讲，非均匀分子固体在猛烈撞击下会引起原子和电子性质的剧烈变化。一个可能的机理就是电子的激发：电子是属于轻、快和量子化的粒子，因此它首先对外部的微扰做出响应。接下来就会电子激发，几步就会引起化学反应的发生。这包括激发态表面的离解、辐射过程以及非辐射能的钝化。,零温下，能带闭合也不能爆炸！,对非完美含能分子晶体进行的模拟研究，如RDX、PETN等含有,TATB,中硝基基团的弯曲振动 受剪切的,TATB,在,Fermi,面附近的能带结构,-HMX,的点电荷静电势,TATB中硝基基团的弯曲振动,3,静高压下含能材料的行为,现在认为,TATB,晶体和硝基甲烷中硝基的弯曲会导致能带的闭合，考察分子内硝基的弯曲在现在实验所能达到的压力范围内能带是否能够闭合可以通过静高压、均匀压力和非轴向加压下对硝基甲烷的电子结构计算，得到其光学带隙来确定。,电子结构的计算可以采用自洽电荷密度泛函,紧束缚理论（,SCCDFTB,）得到。这是基于密度泛函理论下对标准紧束缚模型的拓展。自洽部分描述的是总能、原子力和电荷的输运。实际应用中对于每个固定的晶胞参数集，所有的原子位置均由能量最小化确定，即原子是全弛豫的。,SCCDFTB,的紧束缚特点是可以处理大的超晶胞，可包含更多个分子。早期的第一性原理计算则只能计算较小的超晶胞和较低的精度，但都给出了类似的带隙结果。,3 静高压下含能材料的行为现在认为TATB晶体和硝基甲烷中硝,以硝基甲烷为例，固定比值连续降低初始晶胞参数，压力,P,利用低温公式来估计。对静水压而言，当体积为原来的,50%,时，压力为约,50GPa,，而,HOMO-LUMO,减小了约,0.6eV,，即为原来值的,13%,。,HOMO,和,LUMO,同时增加而其差值却几乎随体积的变化单调下降（参见下左图）。在非常高的静水压下带隙的变化值示下面右图中。压缩率最高达,70%,。,HOMO-LUMO,值的减小，从理论上可以认为感度在增加。,以硝基甲烷为例，固定比值连续降低初始晶胞参数，压力P利用低温,完美硝基甲烷晶体的原胞的,HOMO-LUMO,在静水压下硝基甲烷,HOMO-LUMO,值在静水压和非静水压下随比容的变化的变化值，突然向下的拐点对应高的,C-H,键弯曲,完美硝基甲烷晶体的原胞的HOMO-LUMO 在静水压下硝基,C-H键形变的机理如下：,随着初始,C-N共价单键的收缩，电荷被传输到C原子和N原子之间的区域，随后C-N键变得更强，使得该键变为C=N共价双键。比邻的C-H键作为电子的供源，至少这些H原子中有一个的部分电子被拉走，导致该键变弱。另一方面，局域在硝基附近的轨道并没有被压缩的明显变形。,C-H键形变的机理如下：随着初始C-N共价单键的收缩，电荷被,由,Mulliken电荷布居可以看出，H原子总的Mulliken电荷丢失随每个分子的形变而单调增加。这些发现意味着化学反应和共价键相离子键的转化会导致硝基甲烷的爆轰会始于带隙的闭合之前。最近的计算研究预测TATB金属化的压力在120GPa，远高于爆轰压力。对于硝基甲烷的情况，引起爆轰的原因应是质子的离解。这就有点意思了，在研究含能材料时，应该更注重它在单轴压缩而不是静水压缩下的化学行为，这与Dick对硝基甲烷的研究结果相一致。,V,/,V,0,C,1,H,1,a,H,1,b,H,1,c,C,2,H,2,a,H,2,b,H,2,c,0.79,0.012,-0.010,0.016,0.016,0.012,-0.010,0.016,0.020,0.60,0.126,-0.016,-0.026,-0.026,0.108,-0.019,-0.008,-0.155,0.51,0.294,0.014,-0.074,-0.074,0.168,-0.051,-0.030,-0.245,硝基甲烷分子,1,和,2,中的,C,原子和,H,原子的,Mulliken,电荷（在单晶中沿,b,方向压缩）,由Mulliken电荷布居可以看出，H原子总的Mullike,4 硝胺炸药的降解动力学模拟,对凝聚炸药在高密度和高温下的化学反应机理的详细描述对于理解这些材料在燃烧和爆轰条件下反应面上所发生的事件有重要意义。如在冲击条件下，含能材料会被快速升温到几千度并被压缩到几百个千巴 ，体积的减少可达,30%。这时会引发复杂的化学反应，释放出的大量能量会导致反应持续下去。很明显弄清极端条件下各种化学反应机理对于构建宏观水平上连续爆轰模型有重要意义。,4 硝胺炸药的降解动力学模拟对凝聚炸药在高密度和高温下的化,静态模拟结果,基于自洽密度泛函紧束缚算法（,SCC-DFTB）,产物,H,2,O,、,N,2,、,CO,2,和,CO,的浓度,C(t,)随时间,t,的演化,H,2,O、N,2,、CO,2,和CO的有效反应生成速率分别为：0.48、0.08、0.05和0.11ps,-1,（,a,） （,b,） （,c,） （,d,）,静态模拟结果基于自洽密度泛函紧束缚算法（SCC-DFTB）产,-HMX,撞击加载下的初始分解动力学,图,1,一个分子分解成,C,3,H,6,N,6,O,6,和,CH,2,N,2,O,2,两部分，另一个分子分解成,C,3,H,6,N,8,O,6,；,CH2,和氧原子；,图,2,显示,C,3,H,6,N,6,O,6,进一步分解为,C,2,N,4,O,4,H,4,，然后又和邻近的,HMX,生成的片断,CN,2,O,2,H,2,反应生成,C,3,H,6,N,6,O,6,；,图,3,显示生成物中有,C,2,N,4,O,4,H,4,、,CH,2,、,N,2,O,、,NO,和,CH,2,N,2,O,2,等产物，,图,4,显示生成物中有,C,3,H,6,N,6,O,6,、,C,2,N,3,O,2,H,4,、,C,2,N,2,O,2,H,4,、,NO,2,和,N,2,O,2,；,图,5,显示生成物中有,C,3,H,5,O,4,H,6,、,N,2,O,2,、,CH,2,N,2,O,2,、,CH,2,、,N,2,O,2,和,NO,；,图,6,显示到了反应后期生成物基本是一些小分子产物，而且里边还有一些游离的,O,原子、,H,原子和,N,原子等，这些原子很容易在当系统温度冷却下来的时候，生成,H,2, O,2, N,2,及,H,2,O,等小分子量产物。我们在此模拟中尚未发现有游离碳存在，这可能与我们模拟的时间较短有关，这个问题将在后续的工作中进一步研究。,-HMX撞击加载下的初始分解动力学 图1 一个分子分解成C,材料的计算机模拟方法课件,问题,:,冲击压缩模拟在以下方面非常必要,:,诠释高温高压实验、炸药的冲击感度,沿着冲击路径研究动力学和亚稳态,解决办法,:,新的冲击模拟技术,(MSST),使得在计算速度上大大加快,原子尺度下的量子力学模拟为这种研究提供了工具,得到化学反应的直观图像,冲击条件下的从头算模拟,问题: 冲击压缩模拟在以下方面非常必要:冲击条件下的从头算模,计算原理,基于,Navier-Stokes,方程的,NEMD,方法结合冲击动力学的三个守恒条件，在冲击条件下粒子的速度,(,u,),与冲击波速度,(,v,s,),之间的关系可以写作：,u,=,v,s,(1-,0,/),， （,1,）,Rayleigh,线为冲击波方向上连接初态与终态（,Hugoniot,）的热力学路径：,p,-,p,0,=,v,s,2,0,(1-,0,/),， （,2,）,在冲击条件下遵循能量守恒：,e,-,e,0,=,p,0,(1/,0,- 1/)+,v,s,2,(1-,0,/)/2,， （,3,）,对分子动力学模拟而言，我们采用拉格朗日量：,（,4,）,这里,T,和,V,分别代表单位质量下的动能和势能,Q,为模拟晶胞尺寸下的质量参数。我们可以发现当粒子速度的一阶微分为零时从哈密顿量的形式上方程,(3),和,(4),是等价的，这是因为,T,+,V,=,e,。由,（,5,）,可以看出，当粒子速度的二阶微分为零时，方程,(5),就约化为方程,(2),。因此通过方程,(5),，我们就可以将冲击,Hugoniot,关系结合到非平衡分子动力学中。,计算原理基于Navier-Stokes方程的NEMD方法结合,冲击模拟可以诠释与化学反应相关的热力学途径,只能模拟静态的热力学单态,整个依赖于热力学路径的化学过程可以得到,反应只与冲击波速度和初始条件有关,通过一次冲击模拟可以获得大量的热力学态,静态计算,vs.,冲击模拟,仅为给定温度和压力下的平衡态,根据产物得到数据,冲击模拟可以诠释与化学反应相关的热力学途径只能模拟静态的热力,冲击波速度在,9km/s,时的初始分解与产物演化,反应始于C-H键及N-O键的首先断裂；,随着反应时间的进行，这些碎片会进一步分解成一个,N,2,O,2,、,C,2,N,2,O,2,H,4,、,NO,2,和,C,2,N,3,O,2,H,4,；,进而,CH,2,和,NO,2,从各个碎片中分离出来，分别形成,CN,2,O,2,H,2,、,C,2,NH,4,、,CNH,2,、,N,2,O,、,N,2,O,2,和,CN,4,O,4,H,2,显示生成物中有,N,2,、,O,2,、,NO,2,、,N,2,O,和游离的,O,、,H,、,CH,、,N,、,C,等这些离子在最终会形成,H,2,O,、,CO,、,CO,2,和凝聚碳等。,冲击波速度在9km/s时的初始分解与产物演化 反应始于C-H,材料的计算机模拟方法课件,5 CO在高温高压下的歧化反应,CO气体在自然界大量存在，也是炸药爆轰的主要产物之一。,CO作为炸药的降解产物之一，并且一经产生便存在于高温高压的环境之中，在这样的极端条件下很容易形成更复杂的化合物，因此能否准确描述在特定条件下产物的化学组分以及化学反应过程对于理解炸药爆轰过程以及爆轰产物状态方程具有很重要的意义。,5 CO在高温高压下的歧化反应CO气体在自然界大量存在，也,CO在冲击状态下的聚合分解反应,CO在冲击状态下的聚合分解反应,CO在冲击状态下的聚合分解反应,t=0.16ps,t=0.18ps,CO在冲击状态下的聚合分解反应t=0.16pst=0.18p,气体,CO,在,5000K,25GPa,条件下的一些产物,气体CO在5000K,25GPa条件下的一些产物,6,极端条件下含能固体HMX结构与性能的理论研究,6 极端条件下含能固体HMX结构与性能的理论研究,研究背景：,感度,分子结构,晶体结构,晶格聚集态,缺陷,空穴,起爆机理,起爆过程,效率设计,HMX,研究背景：感度分子结构晶体结构晶格聚集态缺陷空穴起爆机理起爆,研究背景：,Bowden,1,和,Kanllet,2,等人很早就对感度与分子结构之间的关系做了大量的工作，认为认识含能材料的感度不仅要从分子内基团的特征和数目，还必须要深入到分子的电子层结构层次，才能触及问题的本质。一般认为,C-NO,2,键的断裂常常被认为是硝基化合物分子分解过程中非常关键的一步，但有的硝基化合物分子的分解是首先发生消环反应（如黑索金,(RDX),分子）。因而，确定与,C-NO,2,键的强度相关的因素就具有十分重要的意义。,-NO,2,：,致敏,-NH,3,：,致钝,1 Bowden F.P., Yoffe A.D., Initiation and Growth of Explosions in Liquids and Solids,Cambridge University Press, Cambridge, 1952.,2,Kamlet M.J., Proceedings 6th Symposium (International) on Detonation, San Diego,California, ADAO 59120, 1976 .,研究背景：Bowden 1和Kanllet 2等人很,但是仅从分子角度出发得出结论往往并不完全正确，因为不同的晶体构型其感度相差很大，比如,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即奥克托金或环,-HMX,的晶体及椅子型分子构型,研究背景：,包括,2,个,C,4,H,8,N,8,O,8,分子，以,b,为对称轴，属于单斜晶系的,P2,1,/,c,准椅状的,C,4,H,8,O,8,N,8,-HMX 的晶体及椅子型分子构型研究背景：包括2个C4H8,HMX,的晶体结构,包括,8,个,C,4,H,8,N,8,O,8,分子，属于斜方晶系的,F,dd,2,空间群,准船状的,C,4,H,8,O,8,N,8,研究背景：,HMX的晶体结构包括8个C4H8N8O8分子，属于斜方晶,HMX,的晶体结构,包括,6,个,C,4,H,8,N,8,O,8,分子，以,c,轴为对称轴，属于,P,6,1,22,空间结构,研究背景：,HMX的晶体结构包括6个C4H8N8O8分子，以c 轴为,HMX,的晶体结构,包括两个,C,4,H,8,N,8,O,8,0.5H,2,O,分子，具有,Pn,空间对称群,H,2,O,研究背景：,HMX的晶体结构包括两个C4H8N8O80.5H2O分,研究背景：,在实验方面：,由于,HMX,的特殊性，实验条件上受到很大限制。多集中在对,-HMX,的高压压砧实验，温度和压强对,-HMX,的,光谱振动的影响；,还有大量的实验是对,-HMX,研究样品的宏观性质的研究。,在理论方面：,理论方面的研究多集中对在,-HMX,在一定条件下的分子间及分子内力场、振动、动能模型的,研究。而对,HMX,四种构相的系统研究，尤其是其在极端条件下的系统研究就更少。,研究背景：在实验方面： 由于HMX的特殊性，实验,研究内容及方法：,在微观尺度上对四种结构的,HMX,在极端条件（高温高压）下结构的、电子的、光学的及热力学性质进行系统研究，对三种纯,HMX,晶体在高温高压下的相变、力学性质以及,常温常压下,-HMX,基,PBXs,在,不同晶面上的力学性质进行探讨，尝试在介观尺度上对添加粘结剂条件下,-HMX,的相关性质进行讨论。,研究内容：,研究内容及方法： 在微观尺度上对四种结构的HMX在极端,研究内容及方法：,研究方法：,密度泛函理论,经典分子力学方法,耗散粒子动力学方法,使用的软件：,Material Studio,研究内容及方法：研究方法：密度泛函理论经典分子力学方法耗散粒,采用了基于密度泛函理论的平面波赝势法，使用,MS,中的,Castep,模块；,电子之间的交换相互作用尝试使用了广义梯度近似,（,GGA,）,中的,PBE,和局域密度近似,（,LDA,）,中的,CA,交换关联势；,离子实和价电子之间的相互作用采用的是,Vanderbilt,超软赝势；,布里渊区积分采用,Monkhorst-pack,方法，,几何优化中采用,Broyden-Fletcher-Goldfarb-Shanno (BFGS),的算法，收敛标准为总能量偏差小于,110,-6,eV/atom,，原子的最大,Hellmann-Feynman,力小于,0.05eV/,，最大应力偏差小于,0.01GPa,，最大位移偏差小于,0.002 ,。,结果与讨论：,采用了基于密度泛函理论的平面波赝势法，使用MS中的Caste,对,-HMX,采用,423,的,K,点对布里渊区进行采样（网格原点位于,点），平面波基函数的截断能取,500eV,；对,-HMX,采用,214,的,K,点，,500eV,的截断能；对,-HMX,采用,331,的,K,点，,500eV,的截断能；对,-HMX,采用,121,的,K,点，,380eV,的截断能。,在进行深入的模拟之前，对同一体系采用更高的截断能和更密的,K,点采样进行试验，结果表明，这些改变使总能量的变化总小于,0.002eV/atom,。因此，我们所选取的参数对于对,HMX,来说是足够的。,结果与讨论：,对-HMX采用423的K点对布里渊区进行采样（网格原点,结果与讨论：,-HMX,的体积与压强的关系,整个范围内，用,GGA,方法得到的数据比实验值偏高，而用,LDA,方法得到的数据比实验值偏低，相当而言,LDA,方法的结果要好一点。 低压下的差别都很大。,由于压强的增加，分子间距不断减小，分子间相互作用力不断增强，而使得,GGA,方法能更好的适合此时的体系，在,15,到,25GPa,的范围内,GGA,的方法比,LDA,方法的结果好,。,27GPa,时发生相变，我们的结果与实验结果不符。,结果与讨论：-HMX的体积与压强的关系整个范围内，用GGA,结果与讨论：,各个方向晶格参数的变化率和压强间的关系，,b,的压缩性最大，,c,的压缩性最小。表明了在各个方向的可压缩性有很大的差异。,结果与讨论：各个方向晶格参数的变化率和压强间的关系，b的压缩,-HMX,各项异性,结果与讨论,：,-HMX各项异性结果与讨论：,结果与讨论：,N-N,键的压缩率最大,N-O,键的压缩率最小,结果与讨论：N-N键的压缩率最大N-O键的压缩率最小,结果与讨论：,N,7,-O,3,-O,4,几乎不随压强的变化而变化,H,13,-C,2,-N,4,的变化最大,结果与讨论：N7-O3-O4几乎不随压强的变化而变化H13-,结果与讨论：,H,14,-C,2,