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

2024/6/29材料的材料的计算算机模机模拟方法方法2023/7/30材料的材料的计计算机模算机模拟拟方法方法内容内容材料材料计算方法以及最新的算方法以及最新的进展；展；含能材料的引含能材料的引发反反应；静高静高压下材料的行下材料的行为；材料的降解材料的降解动力学模力学模拟；气气态CO在高温高在高温高压下的歧化反下的歧化反应；HMX固体的性能理固体的性能理论模模拟内容材料内容材料计计算方法以及最新的算方法以及最新的进进展；展；1 计算方法算方法简介介在含能材料的模在含能材料的模拟计算中，目前主要算中，目前主要有三种方法：有三种方法：第一性原理（或从第一性原理（或从头算）的算）的电子子结构构计算算(ab initio)、半半经验方法方法分子力学方法分子力学方法(MM)。1 计计算方法算方法简简介在含能材料的模介在含能材料的模拟计拟计算中，目前主要有三种方法算中，目前主要有三种方法方法比方法比较方法方法优点优点缺点缺点最适宜的领域最适宜的领域Ab initioF采用量子物理方法采用量子物理方法F数学上严格：没有数学上严格：没有经验参数经验参数半经验方法半经验方法F采用量子物理方法采用量子物理方法F使用实验参数使用实验参数F使用更多的近似使用更多的近似分子力学方法分子力学方法F采用经典物理学方采用经典物理学方法法F需要有嵌入经验参需要有嵌入经验参数的力场数的力场J适用的体系广泛适用的体系广泛J不需要实验数据不需要实验数据J可以计算激发态和可以计算激发态和过渡态过渡态J与从头算方法相比与从头算方法相比对计算机的要求较对计算机的要求较低低J可计算过渡态和激可计算过渡态和激发态发态J在有限的计算资源在有限的计算资源下计算下计算“便宜便宜”，快且有用快且有用J可以用于大分子如可以用于大分子如酶的计算中酶的计算中L非常耗费机时非常耗费机时L需要有从头算方法需要有从头算方法或实验给出参数或实验给出参数L精度比从头算方法精度比从头算方法低低L不能计算电子的性不能计算电子的性质质L需要从头算或实验需要从头算或实验数据作参数数据作参数L商业软件只能用于商业软件只能用于有限范围的分子有限范围的分子Q小体系（数十个原小体系（数十个原子子几百个原子）几百个原子）Q有电子转移的体系有电子转移的体系Q不需要实验数据的不需要实验数据的体系体系Q需要高精度的体系需要高精度的体系Q中等尺寸的体系中等尺寸的体系(几百个原子几百个原子数千数千个原子个原子)Q有电子转移的体系有电子转移的体系Q大体系（数千大体系（数千数数亿个原子）亿个原子）Q体系中不能包含有体系中不能包含有键断裂的情况键断裂的情况方法比方法比较较方法方法优优点缺点最适宜的点缺点最适宜的领领域域Ab initio适用的体系适用的体系 通通过对材料的模材料的模拟计算可以得到如下的性算可以得到如下的性质：生成生成热、键能和反能和反应能、能、分子（或晶体）的能量和分子（或晶体）的能量和结构构热化学化学稳定性定性)、过渡渡态的能量和的能量和结构（活化能）、构（活化能）、反反应路径以及路径以及动力学、力学、机理、机理、分子中的分子中的电荷分布（反荷分布（反应点）、点）、取代效取代效应、振振动频率（率（红外和外和Raman光光谱）、）、电子的激子的激发跃迁（迁（UV/可可见光光谱）、）、磁屏蔽效磁屏蔽效应（NMR谱）。p 通通过对这些性些性质的的预测，在含能材料中有，在含能材料中有许多多实用用价价值，如：可以研究合成路径、反，如：可以研究合成路径、反应产物、爆物、爆轰产物状物状态方程、爆方程、爆轰引引发的机理等等。的机理等等。通通过对过对材料的模材料的模拟计拟计算可以得到如下的性算可以得到如下的性质质：1.1 从从头算算(ab initio)方法方法从从头算方法在算方法在计算中是基于量子力学基本定律，没有引入算中是基于量子力学基本定律，没有引入经验或或实验参数，参数，仅使用使用5个常数（个常数（电子子质量：量：9.10938215E-31Kg、电子子电量：量：1.602176487E-19C、光速：、光速：c=2.99792458108 m/s、Planck常数：常数：h=6.6260689610-34 Js和波和波尔兹曼常数曼常数1.3806504E-23J/K），），对化学体系利用化学体系利用满足足严格数学近似的基格数学近似的基组来解来解Schrdinger方程：方程：动能是所有粒子的动能加和动能是所有粒子的动能加和动能是所有粒子的动能加和动能是所有粒子的动能加和 势能代表粒子之间的库仑相互作用势能代表粒子之间的库仑相互作用势能代表粒子之间的库仑相互作用势能代表粒子之间的库仑相互作用（即：核与核之间、电子与电子之间（即：核与核之间、电子与电子之间（即：核与核之间、电子与电子之间（即：核与核之间、电子与电子之间的排斥以及电子与核之间的吸引）的排斥以及电子与核之间的吸引）的排斥以及电子与核之间的吸引）的排斥以及电子与核之间的吸引）对电子而言，对电子而言，对电子而言，对电子而言，对电子而言，对电子而言，q qq=-=-=-e ee，对核而言，原子序数为，对核而言，原子序数为，对核而言，原子序数为，对核而言，原子序数为，对核而言，原子序数为，对核而言，原子序数为Z ZZ：q qq=+=+=+ZeZeZe。1.1 从从头头算算(ab initio)方法从方法从头头算方法在算方法在计计算中算中通通过对Schrdinger方程的形式方程的形式观察，我察，我们知道知道对于多于两个以上粒于多于两个以上粒子的分子，子的分子，严格解格解Schrdinger方程是不可能的，必方程是不可能的，必须引入一些必要引入一些必要的合理近似，甚至一些的合理近似，甚至一些经验参数，如非相参数，如非相对论解、解、Born-Oppenheimer近似、近似、单粒子近似、原子粒子近似、原子轨道的道的线性性组合等。合等。对解出解出的的Schrdinger方程是否合适，必方程是否合适，必须满足以下判据：足以下判据：p这个解必个解必须很好地被定很好地被定义且能且能够用来用来说明分子明分子结构和分子中构和分子中电子的状子的状态。p分子的分子的势能按照原子核的位移必能按照原子核的位移必须平滑平滑连续可可变。p模型不能有偏离（即：假模型不能有偏离（即：假设化学化学键只存在于原子之只存在于原子之间）。）。p模型必模型必须与体系大小相关（即解和相关与体系大小相关（即解和相关误差必差必须与分子的大小成比例）与分子的大小成比例），p模型必模型必须是可是可变的（即：近似解必的（即：近似解必须提供体系能量真提供体系能量真值的上界，因此的上界，因此近似解近似解应具有最低能量具有最低能量值与真与真实的波函数相一致）。的波函数相一致）。通通过对过对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及其同事及其同事对基函数作出了很大基函数作出了很大贡献，他献，他们用一用一组高斯函数来描述内高斯函数来描述内层和外和外层电子。子。l如如6-21G 表示内层原子轨道用六个高斯基函数来描述，内层价壳则用两个高表示内层原子轨道用六个高斯基函数来描述，内层价壳则用两个高斯基函数表示，外层价壳用一个高斯函数表示。斯基函数表示，外层价壳用一个高斯函数表示。l类似地还有：类似地还有：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的水平，的水平，p 单取代（取代（CIS）打开一个占据）打开一个占据轨道和一个虚道和一个虚轨道，它等价于道，它等价于单电子激子激发态。p 更高更高级的的计算包括双取代（算包括双取代（CID）对应于打开两于打开两对轨道：道：CISD代表增加单取代和双取代的行列式，代表增加单取代和双取代的行列式，CISD（T）则代表单、双和三激发态。）则代表单、双和三激发态。理论极限为全组态相互作用，它是用理论极限为全组态相互作用，它是用HF行列式的线性组合表达可能取代的行列式。全行列式的线性组合表达可能取代的行列式。全CI方法提方法提供了分子体系最完全的非相对论处理，当然它也非常耗机时。供了分子体系最完全的非相对论处理，当然它也非常耗机时。组态相互作用（组态相互作用（CI）：）：由于在分子的由于在分子的HF波函数里只有一个行列式，因此它在分子中只能描波函数里只有一个行列式，因此它在分子中只能描对典型的分子而言，大部分的基典型的分子而言，大部分的基态能量都源自于能量都源自于单电子子HF的能量。的能量。MP微微扰理理论假假设电子相关的子相关的效效应很小，可以用很小的很小，可以用很小的矫正（即微正（即微扰）来得到）来得到HF的解。的解。MP方法假方法假设分子真正的分子真正的Hamiltonian量量可分成两部分：可分成两部分：（其中：（其中：代表代表单电子的子的贡献；献；代表代表组态相互作用的相互作用的贡献）献）则代表能量和分子波函数的代表能量和分子波函数的展开系数，即：展开系数，即：；代入分子的代入分子的Schrdinger方程中就可以得到：方程中就可以得到：EHF=E(0)+E(1)+E(2)+这就是所就是所谓的的MP2、MP3、MP4、MP5等等。等等。Mller-Plesset 微扰理论微扰理论对对典型的分子而言，大部分的基典型的分子而言，大部分的基态态能量都源自于能量都源自于单电单电子子HF的能量。的能量。这种方法是将种方法是将电子相互作用看作是子相互作用看作是电子密度的函数，子密度的函数，现行的行的DFT方法方法通通过Kohn-Sham方程将方程将电子能量分成子能量分成4个部分：个部分：E=ET+EV+E J+E XC，其中，其中ET 代表代表电子的子的动能能(来自于来自于电子的运子的运动)；EV 代表代表电子的子的势能（包括核与核之能（包括核与核之间的排斥能、核与的排斥能、核与电子之子之间的相互作用能）；的相互作用能）；E J 代表代表电子之子之间的排斥能；的排斥能；E XC代表代表电子相关能。子相关能。p由于由于对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 闭闭壳壳层层和开壳和开壳层计层计算大多数分子具有算大多数分子具有电电子子单线态单线态基基态态（开壳开壳层非限制性自旋非限制性自旋计算适用的体系：算适用的体系：具有奇数个具有奇数个电子的分子子的分子(例如一些离子和自由基例如一些离子和自由基)；激激发态；分子有三分子有三线态或更高自旋的基或更高自旋的基态(例如例如O2)；电子子对的分离的分离过程程(例如分解等例如分解等)；共振体系，存在共振体系，存在电子离域。子离域。开壳开壳层层非限制性自旋非限制性自旋计计算适用的体系：算适用的体系：自旋自旋污染染非限制性自旋波函数并不是自旋算符非限制性自旋波函数并不是自旋算符的本征函数，因此要的本征函数，因此要对自旋自旋态进行混行混合，从而合，从而导致自旋致自旋污染。染。在大多数开壳在大多数开壳层系系统中，自旋中，自旋污染是可染是可以忽略的。然而以忽略的。然而对于未成于未成对电子离域的子离域的分子，自旋分子，自旋污染效染效应相当明相当明显，尤其是，尤其是次高自旋次高自旋态。NONO的非限制性计算，的非限制性计算，的非限制性计算，的非限制性计算，NONO含有一个未成对电子（双线态基态含有一个未成对电子（双线态基态含有一个未成对电子（双线态基态含有一个未成对电子（双线态基态22），当用小基），当用小基），当用小基），当用小基组组组组(STO-3G(STO-3G 或或或或 3-21G)3-21G)计算时就呈现出严重的自旋污染。计算时就呈现出严重的自旋污染。计算时就呈现出严重的自旋污染。计算时就呈现出严重的自旋污染。自旋自旋污污染非限制性自旋波函数并不是自旋算符的本征函数，因此要染非限制性自旋波函数并不是自旋算符的本征函数，因此要对对p在含能材料的在含能材料的计算中，由于包含硝基算中，由于包含硝基(-NO2)，在，在计算算激激发态或有自由基参与的反或有自由基参与的反应里，也存在自旋里，也存在自旋污染染现象。象。p自旋自旋污染效染效应的存在会造成几何的存在会造成几何预测和振和振动频率的率的预测不准确，当然自旋不准确，当然自旋污染效染效应会随基会随基组的增大而降低。的增大而降低。p当占据当占据轨道和非占据道和非占据轨道非常靠近道非常靠近时，自旋，自旋污染就特染就特明明显。如当。如当键被拉伸到快要断裂被拉伸到快要断裂时，这时成成键和反和反键轨道的能量差接近于零，道的能量差接近于零，产生的波函数就包含等量的生的波函数就包含等量的单线态和三和三线态。pGaussian软件包已件包已经在其在其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自洽自洽场场方法方法(HarZerners 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 ZAM1方法方法(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)力；力；势能函数依能函数依赖经验方法如从方法如从实验或从其它或从其它计算方法中得到。算方法中得到。现行的分子力学方法是用一行的分子力学方法是用一组势能函数来表征，用以描述能函数来表征，用以描述化学力化学力力力场，这些力些力场依依赖于：原子位移（即于：原子位移（即键长）、）、原子原子类型和一些与原子型和一些与原子类型及型及键相关的其它参数。相关的其它参数。典型的分子力学方法基本假典型的分子力学方法基本假设设 每个原子都可以用一个具有特定每个原子都可以用一个具有特定质质量量pAMBER(Assisted Model Building with Energy Refinement)力力场,最初用最初用于研究生物分子如蛋白于研究生物分子如蛋白质和核苷；和核苷；pCHARMM(Chemistry at HARvard Molecular Mechanics)力力场，最初也是，最初也是用于生物分子和用于生物分子和药物研究，但也被用来研究胶束和自物研究，但也被用来研究胶束和自组装大分子；装大分子；pMMx(MM2,MM3,等等).用以研究非极性小分子的用以研究非极性小分子的结构构优化和化和热力学研究力学研究MMx 力力场包含包含对标准的准的键和和键角角势能面的二次能面的二次拟合合进行三行三级和四和四级校正，因此它校正，因此它可以可以处理分子振理分子振动中的非中的非谐效效应。不同的。不同的MMx版本其参数也不同，版本越高版本其参数也不同，版本越高的的MMx其其较前一版本而言前一版本而言针对性就越性就越强但适用性就越窄。例如但适用性就越窄。例如较新的新的MM4版版本就不能本就不能够适用于所有种适用于所有种类的分子；的分子；pOPLS(Optimised Potentials for Liquid Simulations)力力场，用于再，用于再现生物生物分子在溶液中的物理性分子在溶液中的物理性质。因因为这些力些力场方法都是方法都是针对生物分子的，生物分子的，这里我里我们不作不作过多的介多的介绍。但。但这些方法些方法对含能分子几何快速的初步含能分子几何快速的初步优化化还是有用的。是有用的。AMBER(Assisted Model Building1.4 分子分子动力学方法力学方法(MD)分子分子动力学方法是假力学方法是假设一个体系有一个体系有N个原子，而体系的状个原子，而体系的状态由由这N个个原子的位置原子的位置ri和和动量量pi或速度或速度vi来来标志，体系的能量志，体系的能量为H(ri,pi)，体系的运体系的运动方程方程为：分子分子动力学的主要目的是解上面的方程求得体系状力学的主要目的是解上面的方程求得体系状态相空相空间演化的演化的轨迹迹ri,pit0，ri,pit1，ri,pit2，ri,pit3，。进而可而可计算我算我们感感兴趣的物理量的趣的物理量的值Q(ri,pi)。在。在实际应用中，我用中，我们常将上面的哈密常将上面的哈密顿方方程化程化为下面的牛下面的牛顿方程，并且用位置方程，并且用位置ri和速度和速度vi做做为描述体系的参量：描述体系的参量：；式中，式中，V(ri)是原子是原子间相互作用相互作用势，通，通过解上面的牛解上面的牛顿方程我方程我们可以得到体系在相空可以得到体系在相空间的的轨迹，迹，进而求得物理量的平均而求得物理量的平均值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.081014sec-1，则振动周期为，则振动周期为T=1/=0.9210-14sec，那么计，那么计算时步长最大可取为约算时步长最大可取为约910-16sec下表给出了各种模拟体系的计下表给出了各种模拟体系的计算步长建议：算步长建议：体系体系体系体系运动状态运动状态运动状态运动状态步长（步长（步长（步长（secsec）简单原子体系简单原子体系简单原子体系简单原子体系移动移动移动移动1010-14-14刚性分子刚性分子刚性分子刚性分子移动、转动移动、转动移动、转动移动、转动510510-15-15软性分子、限制软性分子、限制软性分子、限制软性分子、限制键长键长键长键长移动、转动、移动、转动、移动、转动、移动、转动、扭动扭动扭动扭动210210-15-15软性分子、无键软性分子、无键软性分子、无键软性分子、无键长限制长限制长限制长限制移动、转动、移动、转动、移动、转动、移动、转动、扭动、振动扭动、振动扭动、振动扭动、振动1010-15-15或或或或510510-1616MD 的的计计算步算步长选长选取取计计算步算步长长的的选选取，在取，在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 initio2 含能材料的引含能材料的引发反反应Coffey:含能材料的引含能材料的引发反反应来自于含能材料固来自于含能材料固体内位体内位错的的贯通，当冲通，当冲击波速度足波速度足够高高时，位，位错产生的能量会直接生的能量会直接传输给含能材料分子的内在振含能材料分子的内在振动模式模式 振振动能能up-pumping模型指出：冲模型指出：冲击波会波会产生构生构成晶体的分子的最低振成晶体的分子的最低振动模式所吸收的激模式所吸收的激发声子声子热浴，增加的声子吸收以及分子浴，增加的声子吸收以及分子间振振动能的重新能的重新分配会分配会进一步激一步激发更高的振更高的振动模式，最模式，最终达到平达到平衡衡过渡渡态，从而，从而产生化学生化学键的断裂并的断裂并维持化学反持化学反应。p值得注意的是，由于爆得注意的是，由于爆轰波前波前缘的的传播速度播速度对化化学反学反应而言太快以至于不能而言太快以至于不能驱动它它们热激活。激活。2 含能材料的引含能材料的引发发反反应应Coffey:含能材料的引含能材料的引发发反反应应来自来自Gilman指出，在前指出，在前缘的的压缩会会导致未反致未反应物物质的局部金的局部金属化，属化，这种金属化由共价种金属化由共价键的弯曲或由密度达到一定的弯曲或由密度达到一定时，它它们达到达到临界点后形成。在共价性材料中，当界点后形成。在共价性材料中，当键被拉伸和被拉伸和弯曲弯曲时，就会有一个核，就会有一个核组态空空间区域，区域，这时的最高占据的最高占据轨道（道（HOMO）和最低未占）和最低未占轨道（道（LUMO）的能量差）的能量差值就会就会减小。减小。固体含能材料会固体含能材料会导致金属化致金属化 因因为爆爆轰波的前沿会引起未反波的前沿会引起未反应物物质的的压缩和剪切，特和剪切，特别是是对共价性炸共价性炸药材料而言，材料而言，键的弯曲相的弯曲相对而言比而言比键的的压缩显得更得更为重要，重要，这是因是因为与之相关的力常数会与之相关的力常数会变小。小。剪切剪切导致金属化致金属化这一一观点的主要点的主要证据是：通据是：通过对一系列不一系列不同的氨基同的氨基-硝基硝基-苯炸苯炸药的的HUMO-LUMO差差值的研究，其的研究，其“shake-up”谱与其冲与其冲击感度有直接的关感度有直接的关联 Gilman指出，在前指出，在前缘缘的的压缩压缩会会导导致未反致未反应应物物质质的局部金属化，的局部金属化，在体在体积不不变的情况下，的情况下，键的弯曲会引的弯曲会引起起HUMO-LUMO值变小从而引起在小从而引起在键合合态的的电子的离域化，子的离域化，这样会引起快会引起快速的化学反速的化学反应。对TATB而言，在体而言，在体积接近不接近不变时，临界剪切（弯曲）角与能隙的消失相界剪切（弯曲）角与能隙的消失相对应。在体在体积积不不变变的情况下，的情况下，键键的弯曲会引起的弯曲会引起HUMO-LUMO值变值变小从小从对非完美含能分子晶体非完美含能分子晶体进行的模行的模拟研究，如研究，如RDX、PETN等含有孔洞和位等含有孔洞和位错的情况，都表明当的情况，都表明当RDX晶体中存在晶体中存在单及双孔洞及双孔洞时，对其其进行行压缩它它们的的光学光学带隙会降低，隙会降低，导致它致它们的的绝缘体体金属相金属相转变所需的激所需的激发能下降。能下降。从直从直觉上上讲，非均匀分子固体在猛烈撞，非均匀分子固体在猛烈撞击下会引下会引起原子和起原子和电子性子性质的的剧烈烈变化。一个可能的机理化。一个可能的机理就是就是电子的激子的激发：电子是属于子是属于轻、快和量子化的、快和量子化的粒子，因此它首先粒子，因此它首先对外部的微外部的微扰做出响做出响应。接下。接下来就会来就会电子激子激发，几步就会引起化学反，几步就会引起化学反应的的发生。生。这包括激包括激发态表面的离解、表面的离解、辐射射过程以及非程以及非辐射射能的能的钝化。化。零温下，能带闭合也不能爆炸！零温下，能带闭合也不能爆炸！零温下，能带闭合也不能爆炸！零温下，能带闭合也不能爆炸！对对非完美含能分子晶体非完美含能分子晶体进进行的模行的模拟拟研究，如研究，如RDX、PETN等含有等含有 TATBTATB中硝基基团的弯曲振动中硝基基团的弯曲振动中硝基基团的弯曲振动中硝基基团的弯曲振动 受剪切的受剪切的受剪切的受剪切的TATBTATB在在在在FermiFermi面附近的能带结构面附近的能带结构面附近的能带结构面附近的能带结构-HMX-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共价双共