1，1-二氨基-2，2-二硝基乙烯晶體熱膨脹各向異性的模擬研究

錢(qián) 文，張朝陽(yáng)，宗和厚，熊 鷹，張偉斌，舒遠(yuǎn)杰

（中國(guó)工程物理研究院化工材料研究所 四川綿陽(yáng)621900）

1 引 言

1，1-二氨基-2，2-二硝基乙烯，俗稱(chēng)FOX-7，是1998年由瑞典國(guó)防研究院（Swedish Defence Research Agency，F(xiàn)OI）的一個(gè)研究小組首先報(bào)道合成的[1]，近幾年來(lái)成為高能量密度化合物的研究熱點(diǎn)之一．FOX-7分子（如圖1）屬于推－拉型烯烴，分子內(nèi)存在大的共軛體系，硝基氧原子與氨基氫原子之間還存在氫鍵，含兩個(gè)硝基因此高能而又有兩個(gè)胺基降感，是鈍感高能炸藥（IHE）的理想組分．

圖1 FOX-7的分子構(gòu)型Fig．1 Molecular structure of FOX-7

通常的結(jié)晶條件下均得到α晶型FOX-7晶體，屬于單斜晶系 P21／n空間群[2，15]．FOX-7晶體廣泛應(yīng)用于低感度炸藥配方，炸藥晶體熱膨脹性質(zhì)具有重要意義：材料內(nèi)部晶體熱膨脹的各向異性，將引發(fā)顯著的內(nèi)應(yīng)力[3]，如果內(nèi)應(yīng)力的作用超過(guò)了內(nèi)聚能，將最終導(dǎo)致內(nèi)部裂紋的出現(xiàn)；而材料內(nèi)部晶體受熱的各向異性膨脹，引起材料熱形變[4]，產(chǎn)生不可逆的尺寸或密度變化，將影響精密測(cè)量的準(zhǔn)確度．為了確保含F(xiàn)OX-7炸藥部件的精確測(cè)量和庫(kù)存可靠性，需要對(duì)FOX-7晶體特性進(jìn)行深入研究．

針對(duì)FOX-7晶體的性質(zhì)各國(guó)科學(xué)家在實(shí)驗(yàn)方面開(kāi)展了一系列工作．Bemm 等[2]和 Gilardi等[15]先后用X射線衍射研究了FOX-7的晶體結(jié)構(gòu)，發(fā)現(xiàn)了其層狀結(jié)構(gòu)；Suhithi等[5]用X射線衍射和拉曼光譜研究了在高壓下FOX-7晶體的狀態(tài)與結(jié)構(gòu)變化，發(fā)現(xiàn)壓力在1．1GPa以下時(shí)晶胞b軸的壓縮程度遠(yuǎn)大于a軸和c軸的壓縮程度；Kempa和Herrmann等[6]用 DSC、TG、TMA 等熱分析方法和X射線衍射研究了FOX-7晶體的相變行為，觀察到晶體的α、β、γ三相；Burnham 和 Weese等[7，8]用 DSC、TG、HPC、ARC 及 TG-DTA-FTIR-MS聯(lián)用技術(shù)詳細(xì)研究了FOX-7的固體相變、熱分解和活化能等熱性質(zhì)；Evers等[9]采用X射線單晶和粉末衍射詳細(xì)研究了α、β兩種晶型的FOX-7在200K-393K溫度范圍的晶體特性、熱膨脹和相變．理論研究主要集中在對(duì)FOX-7晶體的各種性能計(jì)算上，第一性原理計(jì)算和分子模擬的結(jié)果不僅與實(shí)驗(yàn)相互印證，更成為深入研究本質(zhì)機(jī)理的有效途徑．趙紀(jì)軍等[10]用密度泛函理論方法計(jì)算了高壓下FOX-7的晶胞參數(shù)和能帶結(jié)構(gòu)，所得結(jié)果與實(shí)驗(yàn)值符合；Hu，Anguang等[11]采用DFT方法分析了FOX-7晶體熱分解反應(yīng)，發(fā)現(xiàn)了在高壓下分子間氫轉(zhuǎn)移導(dǎo)致的化學(xué)分解；Dan C．Sorescu等[12]采用密度泛函理論方法計(jì)算了FOX-7的結(jié)構(gòu)性質(zhì)，結(jié)果與X射線衍射實(shí)驗(yàn)數(shù)據(jù)相符，并通過(guò)自建分子間作用勢(shì)和分子動(dòng)力學(xué)模擬研究了其各向異性性質(zhì)；DeCarlos E．Taylor等[13]以SAPT理論（symmetry adapted perturbation theory）為基礎(chǔ)建立分子力場(chǎng)（勢(shì)函數(shù)），對(duì)FOX-7晶體進(jìn)行等溫等壓下的動(dòng)力學(xué)模擬，結(jié)果與實(shí)驗(yàn)結(jié)果相符并且其隨溫度和壓力變化體現(xiàn)顯著的各向異性；邊亮等[14]采用巨正則系綜蒙特卡洛、密度泛函理論方法結(jié)合分子動(dòng)力學(xué)模擬，詳細(xì)研究了FOX-7晶體的密度、表面能、形貌和爆轟性能．以上研究為深入闡釋FOX-7晶體的熱膨脹機(jī)理及進(jìn)行相關(guān)預(yù)測(cè)提供了新的思路，本文中運(yùn)用分子動(dòng)力學(xué)模擬（MD）結(jié)合量子化學(xué)計(jì)算（QM）方法研究了FOX-7晶體的熱膨脹的各向異性，將能量變化及分子堆積模式與熱膨脹各向異性建立了關(guān)聯(lián)．

2 模擬計(jì)算方法與細(xì)節(jié)

2．1 分子模型的構(gòu)建與分析

FOX-7晶體結(jié)構(gòu)基于α-FOX-7單晶X射線衍射數(shù)據(jù)[15]構(gòu)建，晶胞參數(shù)為a＝0．6940nm，b＝0．6637nm，c＝1．1341nm，β＝90．61°，α＝γ＝90°（295K），呈波浪型π堆積．為在有限的計(jì)算資源下，盡可能獲得能夠代表實(shí)際體系的分子模型，需要以單胞結(jié)構(gòu)為基礎(chǔ)構(gòu)建超晶胞作為分子模擬的初始構(gòu)型．分別構(gòu)建了3＊3＊2超晶胞（圖3）和5＊5＊3超晶胞（圖4），與單胞一起作比較分析．

圖3 FOX-7的3＊3＊2超晶胞結(jié)構(gòu)Fig．3 3＊3＊2supercell of FOX-7

圖4 FOX-7的5＊5＊3超晶胞結(jié)構(gòu)Fig．4 5＊5＊3supercell of FOX-7

2．2 FOX-7晶體熱膨脹的分子動(dòng)力學(xué)模擬

MD模擬計(jì)算結(jié)果的優(yōu)劣很大程度上也決定于力場(chǎng)的合理選擇．FOX-7屬于由碳、氫、氧、氮四種原子組成的硝胺類(lèi)分子晶體，PCFF（The polymer consistent force field）[16－18]及 COMPASS（Condensed-phase Optimized Molecular Potentials for Atomistic Simulation Studies）[19]對(duì)其均應(yīng)有一定的適用性．通過(guò)303K下對(duì)單胞與超晶胞的模擬計(jì)算（溫度外其他參數(shù)設(shè)置如2．3），對(duì)其適用性做了比較．在NPT（等溫等壓系綜）下分子動(dòng)力學(xué)模擬細(xì)節(jié)如：各分子起始速度按Maxwell分布取樣，Velocity Verlet算法[20]進(jìn)行求解，van der waals力和coulomb作用分別用atom-based和ewald方法[21，22]計(jì)算，采用 Anderson法[23]控溫，為實(shí)現(xiàn)對(duì)晶體各向異性形變的模擬，采用Parrinello法[24]控壓，外壓力取常壓10－4GPa，總共進(jìn)行了200ps的MD模擬以保證體系達(dá)到平衡，每500fs取樣一次，共得到400幀的軌跡文件，選取平衡后的軌跡文件用于統(tǒng)計(jì)分析得到晶胞參數(shù)．303K的 MD模擬計(jì)算結(jié)果及與室溫（283K-303K）時(shí)CCDC標(biāo)準(zhǔn)晶胞參數(shù)值[15]的比較見(jiàn)表1．

表1 303K下FOX-7晶胞參數(shù)模擬計(jì)算值及其與標(biāo)準(zhǔn)值的偏差Table 1 Simulation values of cell parameters for FOX-7at 303Kand their deviation

首先，將原胞作為初始構(gòu)型直接進(jìn)行分子動(dòng)力學(xué)模擬，所得到的晶胞參數(shù)模擬值與實(shí)驗(yàn)值的偏差很大，超晶胞由于更接近真實(shí)體系，所得的數(shù)據(jù)偏差較小，可以進(jìn)行進(jìn)一步比較；其次，總體而言PCFF與COMPASS力場(chǎng)的計(jì)算值與實(shí)驗(yàn)值偏差都較小，相對(duì)地在選擇PCFF力場(chǎng)計(jì)算時(shí)得到的晶胞參數(shù)偏差較大，而選擇COMPASS力場(chǎng)計(jì)算時(shí)，晶胞參數(shù)計(jì)算結(jié)果與實(shí)驗(yàn)值的相對(duì)偏差均在3%以內(nèi)，其中5＊5＊3超晶胞的棱邊夾角α、β、γ相對(duì)偏差很小，均在0．5%以內(nèi)，證明COMPASS力場(chǎng)和5＊5＊3超晶胞模型具有較好的適用性．

此外，晶格能被作為驗(yàn)證力場(chǎng)適用性的能量參數(shù)．晶格能可通過(guò)升華熱的文獻(xiàn)值計(jì)算，公式如

其中R是標(biāo)準(zhǔn)氣體常數(shù)，T為熱力學(xué)溫度．FOX-7晶體升華熱文獻(xiàn)值[25]為108．784kJ／mol，計(jì)算所得晶格能為 －113．739kJ／mol；晶格能亦可通過(guò)在選定力場(chǎng)下晶體總能及氣態(tài)分子總能和的差值計(jì)算，近似公式為

對(duì)FOX-7晶體Z＝4，COMPASS力場(chǎng)下所得晶格能為－109．526kJ／mol，可見(jiàn)兩個(gè)晶格能值相差小于10kJ／mol，因此可以認(rèn)為COMPASS力場(chǎng)對(duì)FOX-7晶體具有較好的適用性．

選取COMPASS力場(chǎng)，在NPT系綜下，以FOX-7的5＊5＊3超晶胞模型為初始結(jié)構(gòu)，其他參數(shù)設(shè)置不變，分別完成在303K、318K、333K、348K、363K、378K下對(duì)FOX-7晶體的分子動(dòng)力學(xué)模擬，得到平衡構(gòu)象，實(shí)現(xiàn)對(duì)一定溫

度范圍內(nèi)晶體各向異性熱膨脹的模擬．

2．3 FOX-7晶體隨各向棱長(zhǎng)膨脹的能量變化計(jì)算

分子晶體的熱膨脹性質(zhì)與晶格能相關(guān)．對(duì)于FOX-7晶體，其晶體內(nèi)相互作用越強(qiáng)，則熱膨脹系數(shù)越小，因此可以通過(guò)能量計(jì)算來(lái)進(jìn)一步闡釋熱膨脹各向異性．經(jīng)典的處理方法難以描述晶體內(nèi)部各種復(fù)雜的相互作用，需要應(yīng)用量子理論．晶體結(jié)合能為晶體總能與自由原子的總能之差，計(jì)算晶體總能是關(guān)鍵所在，通過(guò)密度泛函理論方法對(duì)FOX-7晶體進(jìn)行計(jì)算可得到其能量變化情況．

圖5 FOX-7晶體分子動(dòng)力學(xué)模擬的溫度／能量－時(shí)間平衡曲線Fig．5 Temperature／energy-time equilibrium curves for MD simulation of FOX-7crystal

采用基于密度泛函理論贗勢(shì)平面波方法的CASTEP程序[26，27]，原子核和電子之間的相互作用用超軟贗勢(shì)[28]來(lái)描述，截?cái)嗄苋?50eV；雖然LDA（local density approximation）[29，30]是一個(gè)相對(duì)簡(jiǎn)單且簡(jiǎn)化的近似，但它仍能滿足某些交換相干能在原理上就原本應(yīng)該滿足具有的求和規(guī)則，在計(jì)算晶體結(jié)構(gòu)方面獲得了出乎意料的成功，而GGA（Generalized Gradient Approximation）[31]考慮了粒子數(shù)密度的一階梯度，但實(shí)際計(jì)算時(shí)可能出現(xiàn)矯枉過(guò)正的情況，反而不如LDA的結(jié)果好，因此交換關(guān)聯(lián)項(xiàng)采取了LDA近似．通過(guò)對(duì)沿各棱邊方向隨0．1%膨脹率（100%-100．5%，該范圍接近于303K-378K熱膨脹的尺度）的各晶胞模型進(jìn)行幾何構(gòu)型全優(yōu)化和總能計(jì)算，求算膨脹時(shí)各向能量變化率，實(shí)現(xiàn)從能量角度闡述FOX-7晶體熱膨脹各向異性．

3 結(jié)果與討論

3．1 體系平衡判定

體系達(dá)到平衡的標(biāo)志是能量和溫度均趨于平衡，以在303K時(shí)為例，平衡時(shí)的溫度－時(shí)間曲線及能量－時(shí)間曲線分別如圖5所示，后100ps時(shí)體系處于平衡，選取相應(yīng)軌跡文件作為統(tǒng)計(jì)分析的原始數(shù)據(jù)．

3．2 平衡構(gòu)象分析

圖6 303K下FOX-7超晶胞MD模擬后的平衡構(gòu)象Fig．6 Equilibrium structure of FOX-7supercell after MD simulation at 303K

比較各溫度下的平衡構(gòu)象發(fā)現(xiàn)分子在晶體中的取向及堆積方式未見(jiàn)變化．分子內(nèi)鍵角的微小扭轉(zhuǎn)將會(huì)一定程度上影響分子間相互作用，但相對(duì)于溫度的影響，這種偏差微乎其微，下面的數(shù)據(jù)也表明模擬值與實(shí)驗(yàn)值的偏差較小，證明了模擬方法較為準(zhǔn)確可靠．

3．3 MD計(jì)算結(jié)果分析

與303K的MD模擬類(lèi)似，分別對(duì)FOX-7晶體在318K、333K、348K、363K、378K下 MD模擬得到的平衡軌跡文件進(jìn)行統(tǒng)計(jì)分析，最后得到晶體在一定溫度范圍內(nèi)（303K-378K未發(fā)生相變）隨溫度升高的晶胞參數(shù)變化，如圖7和表2所示．

表2 不同溫度下FOX-7晶體的晶胞參數(shù)計(jì)算值Table 2 Simulation values of cell parameters for FOX-7 crystal at different temperatures

圖7 FOX-7的晶胞棱長(zhǎng)與溫度的線性關(guān)系（303K-378K）Fig．7 Linear relation between cell length and temperature for FOX-7（303K-378K）

圖6為303K時(shí)MD模擬平衡后FOX-7超晶胞結(jié)構(gòu)，可以看出FOX-7分子的硝基和碳碳雙鍵平面間的扭轉(zhuǎn)角有微小變化但分子堆積方式不變，

可見(jiàn)，模擬得到的晶胞，隨溫度上升，其棱邊夾角α、β、γ偏差和變化不大，棱長(zhǎng)a、b、c與溫度成良好的線性關(guān)系．

3．4 各向熱膨脹系數(shù)的計(jì)算

為研究晶體熱膨脹的各向異性，根據(jù)微分線膨脹系數(shù)計(jì)算公式

求算各方向的線膨脹系數(shù)，計(jì)算結(jié)果如表3所示：

可見(jiàn)，隨溫度升高FOX-7晶體沿三個(gè)方向均膨脹；線膨脹系數(shù)值與文獻(xiàn)值[9，12，13]相近，雖然本文計(jì)算值與文獻(xiàn)[13]相比較大，但無(wú)論文獻(xiàn)[9]的實(shí)驗(yàn)值還是文獻(xiàn)[12]的理論計(jì)算值都能與本文結(jié)果相互印證；沿b軸方向的線膨脹系數(shù)值αb明顯較大，與文獻(xiàn)[9]中實(shí)驗(yàn)值 12．28×10－5K－1及文獻(xiàn)[12，13]理論計(jì)算值符合很好；而αc與αa的值接近但不相等，能與文獻(xiàn)[9]實(shí)驗(yàn)值和文獻(xiàn)[13]計(jì)算值相互符合．由此可見(jiàn)，本文的FOX-7晶體熱膨脹值和各向異性特點(diǎn)是可信且較佳的結(jié)果．

3．5 熱膨脹各向異性的能量分析

在接近于303K-378K溫度范圍的熱膨脹尺度100%-100．5%，采用DFT方法計(jì)算得到沿各棱邊方向隨0．1%膨脹率的晶胞總能、膨脹時(shí)各向能量變化及對(duì)應(yīng)的各晶胞棱長(zhǎng)，能量值與棱長(zhǎng)的關(guān)系如圖8所示．

隨各棱邊的不斷膨脹，晶胞能量也逐漸增大，膨脹時(shí)的能量差值具有物理意義，如表4所示．

表3 FOX-7晶體熱膨脹系數(shù)值及其與文獻(xiàn)值比較Table 3 Comparisons of CTE values for FOX-7crystal between our result and the literature

圖8 晶胞總能計(jì)算值與棱長(zhǎng)的關(guān)系（棱長(zhǎng)膨脹率0．1%）Fig．8 Correlation between the calculated total energy and the cell length

表4 FOX-7晶體隨各向膨脹的總能變化Table 4 Total energy changes for the expansion along each direction for FOX-7crystal

3．6 晶體熱膨脹各向異性與分子堆積的關(guān)聯(lián)

首先，F(xiàn)OX-7晶體各向總能量變化與不同膨脹率的棱長(zhǎng)成線性關(guān)系，即E∝L，而各棱長(zhǎng)變化與溫度成良好的線性關(guān)系，即L∝T，可以推知E∝T，即各向能量變化與溫度成線性關(guān)系，從而將各向能量變化與溫度進(jìn)行關(guān)聯(lián)．

由沿三個(gè)棱邊方向的能量（單位：kJ／mol）與棱長(zhǎng)（單位：nm）的線性關(guān)系，以及棱長(zhǎng)與熱力學(xué)溫度（單位：K）的線性關(guān)系可以聯(lián)立得到：

將各組數(shù)據(jù)代入經(jīng)整理可得：

再將具體的溫度和能量值代回方程（3．2）（3．3）（3．4）便可得到FOX-7晶體熱膨脹系數(shù)與沿各棱邊方向的能量變化的關(guān)系式，如303K下

從而建立FOX-7晶體能量變化與熱膨脹系數(shù)的關(guān)聯(lián)．

4 結(jié) 論

和c兩方向相互作用不同，因此受熱膨脹過(guò)程沿各棱邊方向的能量變化各不相同，導(dǎo)致b方向熱膨脹最為顯著，而a和c方向熱膨脹系數(shù)不同．

[1]Latypov N V，Bergman J，Langlet A，etal．Synthesis and reactions of 1，1-diamino2，2-dinitroeth-ylene[J]．Tetrahedron，1998，54：11525．

[2]Bemm U，?stmark H．1，1-Diamino-2，2-dinitroethylene：a novel energetic material with infinite layers in two dimensions[J]．ActaCrystallogr．，1998，C54：1997．

[3]Hsueh C H，Becher PF．Thermal stresses due to thermal expansion anisotropy in materials with preferred orientation[J]．JournalofMaterialsScience Letters，1991，19（10）：1165．

[4]Boas W，Honeycombe R W K．Theanisotropy of thermal expansion as a cause of deformation in metals and alloys[A]．ProceedingsoftheRoyalSocietyofLondon．SeriesA，MathematicalandPhysicalSciences[C]，（Feb．25），1947，188（1015）：427．

[5]Suhithi M，Chak P，Maija M，etal．Equation of state and structural changes in diaminodinitroethylene from experimental studies and ab initio quantum calculations[A]．12thDetonation（Int．）Symposium[C]．，SanDiego，California，2002．

[6]Kempa P B，Herrmann M，MolinaMetzger F J，et al．Phase transitions of FOX-7studied by X-ray diffraction and thermal analysis[A]．35thInt．Annu．Conf．ICT[C]，2004．

[7]Burnham A K，Weese R K，Wang R，etal．Thermalproperties of FOX-7[A]，36thInternationalAnnualConferenceofICT＆32ndInternationalPyrotechnicsSeminar[C]，Karlsruhe，Germany，June 28，2005through July 1，2005．

[8]Burnham A K，Weese R K，Wang R，etal．Solidsolid phase transition kinetics of FOX-7[A]，2005 NATASAnnualConference[C]．UniversalCity，CA，UnitedStates，September 18-21，2005．

[9]Evers J，Klap?tke T M，Mayer P，etal．α-andβ－FOX-7，polymorphs of a high energy density material，studied by X-ray single crystal and powder investigations in the temperature range from 200to 423K[J]．Inorganicchemistry，2006，45（13）：4996．

[10]Zhao J J，Liu H．High-pressure behavior of crystalline FOX-7by density functional theory calculations[J]．ComputationalMaterialsScience，2008，42：698．

[11]Hu A G，Larade B，Hakima A R，etal．A first principles density functional study of crystalline FOX-7chemical decomposition process under external pressure[J]．Propellants，Explosives，Pyrotechnics，2006，31（5）：355．

[12]Sorescu D C，Boatz J A，Thompson D L．Classical and quantum-mechanical studies of crystalline FOX-7（1，1-diamino-2，2-dinitroethylene）[J]．Journal ofPhysicalChemistryA，2001，105（20）：5010．

[13]Taylor D E，Rob F R，Betsy M，etal．A molecular dynamics study of 1，1-diamino-2，2-dinitroethylene（FOX-7）crystal using a symmetry adapted perturbation theory-based intermolecular force field[J]．PhysicalChemistryChemicalPhysics，2011，13（37）：16629．

[14]Bian L，Shu Y J，Li H R．Computational investigation on the surface electronic density，morphology and detonation property of 1，1-diamino-2，2-dinitroethylene（FOX-7）crystal[J]．ChineseJ．Struct．Chem．，2012，31（12）：1736．

[15]Gilardi R．Cambridge crystallographic data centre，CCDC127539，1999．

[16]Sun H，Mumby S J，Maple J R，etal．An ab initio CFF93all-atom forcefield for polycarbonates[J]．J．Am．Chem．Soc．，1994，116：2978．

[17]Sun H．Ab initio calculations and forcefield development for computer simulation of polysilanes[J]．Macromolecules，1995，28：701．

[18]Sun H，Rigby D．Polysiloxanes：ab initio forcefield and structural，conformational and thermophysical properties[J]．Spectrochim．Acta，PartA，1997，53：1301．

[19]Sun H．COMPASS：AnabInitioforce-field optimized for condensed-phase applications[J]．J．Phys．Chem．B．，1998，102（38）：7338．

[20]Verlet L．Computer experiments on classical fluids．I．thermodynamical properties of Lennard-Jones molecules[J]，Phys．Rev．，1967，159：98．

[21]Ewald P P．Die Berechnung optischer und elektrostatischer Gitterpotentiale[J]．Ann．Phys．Leipzig，1921，64：253．

[22]Karasawa N，Goddard W A．Acceleration of convergence for lattice sums[J]．J．Phys．Chem．，1989，93：7320．

[23]Andersen H C．Molecular dynamics simulations at constant pressure and／or temperature[J]．J．Phys．Chem．，1980，72：2384．

[24]Parrinello M，Rahman A．Polymorphic transitions in single crystals：a new molecular dynamics method[J]．J．Appl．Phys．，1981，52：7182．

[25]Politzer P，Concha M C，Grice M E，etal．Computational investigation of the structures and relative stabilities of amino／nitro derivatives of ethylene[J]．J．Mol．Struct．（THEOCHEM），1998，452：75．

[26]Segall M D，Lindan P J D，Probert M J，etal．First-principles simulation：ideas，illustrations，and the CASTEP code[J]．J．Phys．：CondensedMater，2002，14：2717．

[27]Clark S J，Segall M D，Pickard C J，etal．First principles methods using CASTEP[J]．Zeitschrift fuerKristallographie，2005，220（5-6）：567．

[28]Vanderbilt D．Soft self-consistent pseudopotentials in a generalized eigenvalue formalism[J]．Phys．Rev．B，1990，41：7892．

[29]Ceperley D M，Alder B J．Ground state of the electron gas by a stochastic method[J]．Phys．Rev．Lett．，1980，45：566．

[30]Perdew J P，Zunger A．Self-interaction correction to density-functional approximations for many-electron systems[J]．Phys．Rev．B，1981，23：5048．

[31]Perdew J P，Burke K，Ernzerhof M．Generalized gradient approximation made simple[J]，Phys．Rev．Lett．，1996，77：3865．