鋁相變臨界壓力的分子動(dòng)力學(xué)計(jì)算

鋁相變臨界壓力的分子動(dòng)力學(xué)計(jì)算 內(nèi)容提要： 本文用分子動(dòng)力學(xué)模擬的方法，采用了EAM勢(shì)函數(shù), 模擬計(jì)算了fcc、hcp、bcc結(jié)構(gòu)鋁在壓縮過程中單原子能量隨單原子體積的變化規(guī)律，通過此規(guī)律分析繪制了fcc、hcp、bcc結(jié)構(gòu)鋁的冷能曲線，并由冷能曲線分析判斷了鋁在高壓下可能出現(xiàn)的相變規(guī)律，計(jì)算了各相變所對(duì)應(yīng)的相變臨界壓力點(diǎn)。計(jì)算結(jié)果表明，鋁在高壓下存在fcc-hcp和hcp-bcc相變，對(duì)應(yīng)的相變臨界壓力分別為61 Gpa和210 Gpa。 關(guān)鍵詞： 鋁，相變，分子動(dòng)力學(xué)，勢(shì)函數(shù)，冷能曲線。 Molecular Dynamics Calculation of Phase Transition Critical Pressure of Aluminum Abstract: By molecular dynamics simulations employing an embedded atom method potential, we investigate structural transformations in single crystal Al. In our simulation, the atomic energy as the function of the single atomic volume of the fcc, hcp, bcc structure of aluminum in the compression process were calculated. The cold energy curve of the fcc, hcp, bcc structure of aluminum were plotted. The results show that the fcc-hcp and hcp-bcc transformation will occur in the compression process, and the transition pressures are 61 Gpa and 210 Gpa respectivly. Keywords: Aluminum , phase transition , Molecular dynamics , potential function , cold energy curve . 引言： 鋁,化學(xué)符號(hào)是Al，它的原子序數(shù)是13。鋁元素在地殼中的含量?jī)H次于氧和硅，居第三位，是地殼中含量最豐富的金屬元素。在金屬品種中，僅次于鋼鐵，為第二大類金屬。至19世紀(jì)末，鋁才嶄露頭角，成為在工程應(yīng)用中具有競(jìng)爭(zhēng)力的金屬，且風(fēng)行一時(shí)。航空、建筑、汽車三大重要工業(yè)的發(fā)展，要求材料特性具有鋁及其合金的獨(dú)特性質(zhì)，這就大大有利于這種新金屬鋁的生產(chǎn)和應(yīng)用。 鋁的應(yīng)用極為廣泛。 鋁是一種非常年輕的金屬，從1944年到1988年的44年間，世界鋁年產(chǎn)量增加了19倍，鋁發(fā)展如此之快，說明鋁本身有很多優(yōu)點(diǎn)，因此鋁不斷在各個(gè)領(lǐng)域中代鋼代銅代木也代陶土，鋁就發(fā)展的一個(gè)重要形式就是“取而代之”。 在鋁短暫的發(fā)展過程中，需求起了巨大的推動(dòng)作用，如二戰(zhàn)時(shí)期由于制造飛機(jī)，推動(dòng)鋁急劇發(fā)展，1940年到1943年鋁的產(chǎn)量由23萬噸猛增到146萬噸。而1960～1970年間鋁由446萬噸增加至1020萬噸，其原因是建筑業(yè)廣泛以鋁代鋼引起的。汽車工業(yè)在發(fā)達(dá)國(guó)家是重要的支柱工業(yè)，當(dāng)前汽車工業(yè)激烈競(jìng)爭(zhēng)中一個(gè)新的方向即汽車“輕量化”－“以鋁代鋼”是必然趨勢(shì)。 隨著鋁在現(xiàn)代生產(chǎn)生活中的大量應(yīng)用，對(duì)鋁材料性能的研究也逐漸深入，鋁在各種特殊環(huán)境下的性能特征成了研究的重點(diǎn)，尤其是鋁在高壓下的相變研究得到了廣泛的關(guān)注。 相變是指當(dāng)外界約束（溫度或壓力）作連續(xù)變化時(shí)，在特定條件（溫度或壓力達(dá)到某定值）下，物相突然發(fā)生改變。在通常條件下物質(zhì)有固、液、氣三種狀態(tài)，分別稱為固相、液相和氣相，在極高溫高壓下還存在等離子態(tài)。從廣義上講，構(gòu)成物質(zhì)的原子（或分子）的聚合狀態(tài)（相狀態(tài)）發(fā)生變化的過程均稱為相變，如從液相到固相的凝固過程、從液相到氣相的汽化過程。金屬和陶瓷等固態(tài)材料在溫度和壓力改變時(shí)，其內(nèi)部組織或結(jié)構(gòu)會(huì)發(fā)生變化，即發(fā)生從一種相狀態(tài)到另一種相狀態(tài)的轉(zhuǎn)變，這種轉(zhuǎn)變稱為固態(tài)相變[1]。 對(duì)于相變的研究，一般情況下可以從熱力學(xué)、結(jié)構(gòu)變化以及相變動(dòng)力學(xué)幾方面進(jìn)行分類。其中熱力學(xué)方面主要從相變臨界條件方面考慮；結(jié)構(gòu)變化方面主要研究相變的產(chǎn)物；而相變動(dòng)力學(xué)主要指的是相變的機(jī)制等。本文研究的鋁的相變屬于固態(tài)相變的范圍。在許多情況下，材料的性能決定于材料中的某些特定結(jié)構(gòu)和微觀組織結(jié)構(gòu)和分布。這種組織結(jié)構(gòu)包括晶體的形狀和分布（晶粒和金相組織）、構(gòu)成的晶體結(jié)構(gòu)及其中的晶體缺陷、電子組態(tài)等，當(dāng)然也包括它們當(dāng)中的組織缺陷。因此研究固態(tài)相變對(duì)控制金屬、合金以及某些非金屬材料性能有極為重要的理論和實(shí)踐意義。 引起固態(tài)相變的原因非常多，其中很重要的一種就是高壓引起的相變。在理論研究方面，人們認(rèn)為鋁在高壓下存在結(jié)構(gòu)相變。在20世紀(jì)80年代Mcmahan和Moriarity等[2]用GPT(The Generalized Pseudopotential Theory)和LMTO (The Linear Muffin-Tin-Orbitals)兩種方法分別計(jì)算了鋁的面心晶格結(jié)構(gòu)(fcc)、六角密堆積結(jié)構(gòu)(hcp)和體心晶格結(jié)構(gòu)(bcc)，得到fcc-hcp以及hcp-bcc高壓結(jié)構(gòu)相變時(shí)的臨界壓力值分別是360 GPa和560 GPa(GPT)，120 GPa和200 GPa(LMTO)；Boettger和Trickey[3]用LCGTO(The Linear Combination of Gaussian-Type Orbitals)方法采用LDA近似后得到fcc-hcp以及hcp-bcc高壓結(jié)構(gòu)相變時(shí)的臨界壓力值分別是(205±20)GPa和(565±60)GPa，同時(shí)得到fcc-bcc高壓結(jié)構(gòu)相變時(shí)臨界壓力值是330 GPa。Lam和Cohen[4]利用從頭計(jì)算贗勢(shì)平面波(PW)的方法得到fcc-hcp及hcp-bcc高壓結(jié)構(gòu)相變時(shí)的臨界壓力值分別是200 GPa和400 GPa，同時(shí)得到fcc-bcc高壓結(jié)構(gòu)相變時(shí)的臨界壓力值為300 GPa。Mishrab 和 Chaturvedi[5]經(jīng)過計(jì)算后認(rèn)為fcc-bcc高壓結(jié)構(gòu)相變時(shí)的臨界壓力值為275 Gpa?？梢?，由不同的理論方法計(jì)算而來的相變點(diǎn)的差異還是比較大的。實(shí)驗(yàn)上對(duì)于鋁相變的研究也有很多。1988年，Greene[6]等人用DAC(Diamond-Anvil—Cel1)方法研究了鋁在高壓下的相變，在當(dāng)時(shí)的實(shí)驗(yàn)條件，壓力加載最大只能達(dá)到220 Gpa左右。在這種實(shí)驗(yàn)條件下，Greene等人發(fā)現(xiàn)在小于219 GPa(V／V0<0．5)時(shí)fcc結(jié)構(gòu)的鋁并沒有發(fā)生相變。隨著實(shí)驗(yàn)手段的不斷進(jìn)步，在2006年，Akahama[7]等人用DAC(Diamond-Anvil—Cel1)方法將壓力值提高到 300Gpa附近并考察了鋁在這一加載過程的結(jié)構(gòu)穩(wěn)定性。在這一工作中，Akahama等人發(fā)現(xiàn)鋁在 217±10 GPa附近發(fā)生了fcc-hcp的結(jié)構(gòu)相變。這是人們首次在實(shí)驗(yàn)中觀察到了鋁在高壓下的結(jié)構(gòu)相變現(xiàn)象，并且相變臨界壓力與Lam 和 Boettger的理論計(jì)算工作符合的非常好，這在一定程度上激發(fā)了人們對(duì)鋁高壓相變的研究熱情。 事實(shí)上，要充分了解相變的微觀機(jī)制和建立相應(yīng)的理論，就必須從微觀層面上研究結(jié)構(gòu)相變，對(duì)結(jié)構(gòu)相變做微觀觀測(cè)。幾十年來，雖然結(jié)構(gòu)相變過程及其相變機(jī)制受到了廣泛重視，但由于結(jié)構(gòu)相變過程是一個(gè)突變過程，轉(zhuǎn)變往往在ps量級(jí)完成，因而在實(shí)驗(yàn)上很難詳細(xì)觀察到相關(guān)的相變機(jī)制，尤其是微觀（原子）尺度的變化機(jī)制。相變機(jī)制往往是通過相變時(shí)一些宏觀量的變化特征推測(cè)得到的。這一情況隨著近年來一些實(shí)驗(yàn)技術(shù)的飛速發(fā)展有所改觀。比如，利用超快激光技術(shù)來研究結(jié)構(gòu)相變。利用飛秒激光產(chǎn)生的納秒以下脈寬的X光可以在材料發(fā)生結(jié)構(gòu)相變時(shí)實(shí)時(shí)測(cè)量晶體結(jié)構(gòu)的變化，這一技術(shù)在研究鐵的沖擊相變中得到了成功應(yīng)用。研究結(jié)構(gòu)相變機(jī)制的另外一個(gè)非常重要的工具就是分子動(dòng)力學(xué)計(jì)算機(jī)模擬。分子動(dòng)力學(xué)計(jì)算機(jī)模擬是研究分子、原子系統(tǒng)的有力工具。它具有溝通宏觀與微觀的作用，其模擬結(jié)果對(duì)理論分析和實(shí)驗(yàn)都是有力的補(bǔ)充，因此它在物理學(xué)、化學(xué)、生物學(xué)和材料學(xué)等許多領(lǐng)域中得到廣泛應(yīng)用[8-11]。當(dāng)實(shí)驗(yàn)研究方法不能完全滿足研究工作的需求時(shí)，用分子動(dòng)力學(xué)計(jì)算機(jī)模擬可以提供實(shí)驗(yàn)上尚無法獲得或很難獲得的重要信息；雖然分子動(dòng)力學(xué)計(jì)算機(jī)模擬不能取代實(shí)驗(yàn)，但可以用來指導(dǎo)實(shí)驗(yàn)，并驗(yàn)證某些理論上的假設(shè)，以此促進(jìn)理論和實(shí)驗(yàn)的發(fā)展。特別是結(jié)構(gòu)相變中有許多與原子有關(guān)的微觀細(xì)節(jié)，在實(shí)驗(yàn)中不容易獲得，而在計(jì)算機(jī)模擬中卻可以方便地得到。這種優(yōu)點(diǎn)使分子動(dòng)力學(xué)計(jì)算機(jī)模擬在研究結(jié)構(gòu)相變時(shí)顯得非常有吸引力。 隨著計(jì)算機(jī)計(jì)算速度與計(jì)算容量的不斷提高，現(xiàn)有的大規(guī)模并行分子動(dòng)力學(xué)模擬可以處理多達(dá)百萬的原子系統(tǒng)，在模擬結(jié)構(gòu)相變的過程中，能夠獲取原子位置變化的詳細(xì)信息，模擬結(jié)果可以和實(shí)驗(yàn)結(jié)果相校驗(yàn)，因此成為在微觀尺度研究結(jié)構(gòu)相變的有力工具。本文即是用分子動(dòng)力學(xué)計(jì)算機(jī)模擬方法研究了鋁的相變臨界壓力等高壓相變規(guī)律。 一．研究方法簡(jiǎn)介 1.分子動(dòng)力學(xué) （1）分子動(dòng)力學(xué)簡(jiǎn)介 研究材料的性質(zhì)需要更深入的認(rèn)識(shí)材料的物理特性和物理規(guī)律，這就要求研究領(lǐng)域必須從傳統(tǒng)的宏觀層次逐步轉(zhuǎn)入細(xì)觀乃至微觀層次，以便考察材料的微觀作用機(jī)理。受目前實(shí)驗(yàn)手段在時(shí)間和空間分辨率上的限制，現(xiàn)在還無法通過傳統(tǒng)的實(shí)驗(yàn)方法直接獲取原子尺度上的信息。分子動(dòng)力學(xué)較好地解決了這個(gè)問題。 分子動(dòng)力學(xué)(Molecular Dynamics，MD)是一門結(jié)合物理，數(shù)學(xué)和化學(xué)的綜合技術(shù)。分子動(dòng)力學(xué)是一套分子模擬方法，該方法主要是依靠牛頓力學(xué)來模擬分子體系的運(yùn)動(dòng)，以在由分子體系的不同狀態(tài)構(gòu)成的系綜中抽取樣本，從而計(jì)算體系的構(gòu)型積分，并以構(gòu)型積分的結(jié)果為基礎(chǔ)進(jìn)一步計(jì)算體系的熱力學(xué)量和其他宏觀性質(zhì)。 分子動(dòng)力學(xué)(Molecular Dynamics，MD) 模擬是指對(duì)于原子核和電子所構(gòu)成的多體系統(tǒng)，把其中每一個(gè)原子核視為在全部其它原子核和電子作用下運(yùn)動(dòng)，在此基礎(chǔ)上求解運(yùn)動(dòng)方程。通過分析各粒子的受力情況，用經(jīng)典或量子的方法求解每個(gè)粒子在某一時(shí)刻的位置和速度，確定系統(tǒng)中各粒子的運(yùn)動(dòng)狀態(tài)，然后使用一定的統(tǒng)計(jì)方法計(jì)算出系統(tǒng)的力學(xué)、熱力學(xué)、動(dòng)力學(xué)等性質(zhì)。 分子動(dòng)力學(xué)計(jì)算機(jī)模擬是研究分子、原子系統(tǒng)的有力工具。分子動(dòng)力學(xué)模擬具有溝通宏觀特性與微觀的作用，對(duì)于許多在理論分析和實(shí)驗(yàn)觀察上都難以觀察的現(xiàn)象做出一定的解釋，因此它在物理學(xué)、化學(xué)、生物學(xué)和材料學(xué)等許多領(lǐng)域中得到廣泛應(yīng)用。分子動(dòng)力學(xué)是由Alder和Wainwrigh[12]在1957年首創(chuàng)的，他們首先在硬球模型下，采用經(jīng)典分子動(dòng)力學(xué)模擬方法研究了氣體和液體的狀態(tài)方程，取得了嘗試性的成功。但是由于計(jì)算機(jī)速度和內(nèi)存的限制，早期模擬的空間尺度和時(shí)間尺度都受到了很大的限制。計(jì)算機(jī)技術(shù)的飛速發(fā)展和多體勢(shì)函數(shù)的提出與發(fā)展，為分子動(dòng)力學(xué)模擬技術(shù)注入了新的活力。1972年，Lees和Edwards[13]首先將分子動(dòng)力學(xué)模擬應(yīng)用到了非平衡的研究，進(jìn)一步擴(kuò)展了分子動(dòng)力學(xué)模擬方法的應(yīng)用范圍。它提供了一種研究高溫高壓下的液體和固體平衡特性的切實(shí)可行的方法，它不僅有助于理論研究，而且有助于解釋實(shí)驗(yàn)中出現(xiàn)的現(xiàn)象，它是研究復(fù)雜的凝聚態(tài)系統(tǒng)的有力工具，是對(duì)物理系統(tǒng)的確定的微觀描述，這個(gè)系統(tǒng)可以是少體系統(tǒng)，也可以是多體系統(tǒng)，并可通過哈密頓量、拉格朗日量或牛頓運(yùn)動(dòng)方程來描述。由于分子動(dòng)力學(xué)模擬是用運(yùn)動(dòng)方程來計(jì)算系統(tǒng)的性質(zhì)，這一技術(shù)既能得到原子的運(yùn)動(dòng)軌跡，又能像做實(shí)驗(yàn)一樣進(jìn)行各種觀察，其結(jié)果既有系統(tǒng)的靜態(tài)特性，也有系統(tǒng)的動(dòng)態(tài)特性。對(duì)于平衡系統(tǒng)，可以在一個(gè)分子動(dòng)力學(xué)觀察時(shí)間內(nèi)作時(shí)間平均來計(jì)算物理量的統(tǒng)計(jì)平均值；對(duì)于非平衡系統(tǒng)，發(fā)生在一個(gè)分子動(dòng)力學(xué)觀察時(shí)間內(nèi)（一般為幾十個(gè)ps）的物理現(xiàn)象都可以用分子動(dòng)力學(xué)進(jìn)行仿真。特別是許多與原子有關(guān)的微觀細(xì)節(jié)，在實(shí)際實(shí)驗(yàn)中無法獲得，而在分子動(dòng)力學(xué)模擬中可以方便清晰地得到，這對(duì)理論和實(shí)驗(yàn)無疑是有力的補(bǔ)充。 （2）分子動(dòng)力學(xué)原理 分子動(dòng)力學(xué)(Molecular Dynamics，MD) 模擬是指對(duì)于原子核和電子所構(gòu)成的多體系統(tǒng)，通過經(jīng)典力學(xué)計(jì)算系統(tǒng)中各個(gè)原子的運(yùn)動(dòng)軌跡，然后使用統(tǒng)計(jì)方法計(jì)算出系統(tǒng)的力學(xué)、熱力學(xué)、動(dòng)力學(xué)性質(zhì)的方法。 分子動(dòng)力學(xué)方法中有兩個(gè)基本假設(shè)，所有的分子動(dòng)力學(xué)的模擬計(jì)算都是在這兩個(gè)基本假設(shè)下進(jìn)行的。這兩個(gè)基本假設(shè)是： 1、所有粒子的運(yùn)動(dòng)都遵循經(jīng)典牛頓運(yùn)動(dòng)定律； 2、所有粒子間的相互作用都滿足疊加原理。 這兩個(gè)假設(shè)意味著分子動(dòng)力學(xué)雖然是在原子層次研究問題，但它忽略了量子效應(yīng)，是一種近似計(jì)算模型。 在分子動(dòng)力學(xué)中，首先將由N個(gè)粒子構(gòu)成的系統(tǒng)抽象為N個(gè)相互作用的質(zhì)點(diǎn)，每個(gè)質(zhì)點(diǎn)具有相應(yīng)的坐標(biāo)、質(zhì)量、及成鍵方式，按體系溫度根據(jù)Boltzmann分布隨機(jī)給定各質(zhì)點(diǎn)的初始速度，然后根據(jù)所選用的勢(shì)函數(shù)來描述質(zhì)點(diǎn)間的相互作用。接著依據(jù)牛頓力學(xué)、根據(jù)所選用的勢(shì)函數(shù)來計(jì)算各質(zhì)點(diǎn)的加速度及速度，得到在系統(tǒng)演化時(shí)各時(shí)刻各質(zhì)點(diǎn)的坐標(biāo)和速度，設(shè)定時(shí)間間隔對(duì)這些計(jì)算結(jié)果進(jìn)行保存。最后可以對(duì)這些結(jié)果進(jìn)行結(jié)構(gòu)、能量、熱力學(xué)、動(dòng)力學(xué)等屬性特征的分析，得到感興趣的計(jì)算結(jié)果，從而得到系統(tǒng)在演化過程中的結(jié)構(gòu)、能量、溫度等屬性的變化特征。其優(yōu)點(diǎn)在于系統(tǒng)中粒子的運(yùn)動(dòng)有正確的物理依據(jù)，準(zhǔn)確性高，并且可同時(shí)獲得系統(tǒng)的動(dòng)態(tài)與熱力學(xué)統(tǒng)計(jì)信息，可廣泛地用于各種系統(tǒng)及各類特性的探討。缺點(diǎn)是模擬移動(dòng)時(shí)間間隔不能過長(zhǎng)，通常為1飛秒。按時(shí)間步長(zhǎng)為1 飛秒計(jì)算，粒子移動(dòng)103次模擬時(shí)間僅為1納秒，盡管如此，對(duì)納米尺度的系統(tǒng)進(jìn)行模擬，分子動(dòng)力學(xué)有著無可比擬的優(yōu)勢(shì)。 分子動(dòng)力學(xué)模擬的實(shí)際步驟可以劃分為四步： 1、設(shè)定模擬所采用的模型，給定初始條件；使得體系粒子位置、速度要與模擬系統(tǒng)相容。 2、計(jì)算作用到每個(gè)原子上的力。這一步是模擬的關(guān)鍵，在這其中要看所用到的勢(shì)函數(shù)是否合理，能否精確地描述體系的性質(zhì)。 3、趨于平衡的計(jì)算過程；在這一步中我們要用一些系統(tǒng)控制方法使得體系符合模擬所要求的壓力、溫度等。并且在程序的設(shè)計(jì)中算法要求精確、高效、簡(jiǎn)潔。 4、宏觀物理量的得出及感興趣量的提??；分子動(dòng)力學(xué)計(jì)算的結(jié)果是體系中粒子的位置、速度、加速度等量，我們需要從結(jié)果中分析全部或局部體系的狀態(tài)，統(tǒng)計(jì)出相應(yīng)的宏觀量和體系特征。 （3）分子動(dòng)力學(xué)模擬的基本方程 分子動(dòng)力學(xué)假定原子的運(yùn)動(dòng)服從確定的描述，原子的運(yùn)動(dòng)和確定的軌跡聯(lián)系在一起。分子動(dòng)力學(xué)方法中描述原子的運(yùn)動(dòng)方程主要有兩種：拉格朗日運(yùn)動(dòng)方程和哈密頓運(yùn)動(dòng)方程。 1．拉格朗日（Lagrangian）運(yùn)動(dòng)方程 對(duì)于由N個(gè)自由度的相互作用的質(zhì)點(diǎn)構(gòu)成的運(yùn)動(dòng)系統(tǒng)，拉格朗日方程為： ， i=1,2,…,n (1-1) 其中q 為廣義坐標(biāo)，代表質(zhì)點(diǎn)的空間位置，為質(zhì)點(diǎn)位置對(duì)時(shí)間的導(dǎo)數(shù) 由n 個(gè)原子構(gòu)成的模擬系統(tǒng)的拉格朗日方程在笛卡兒坐標(biāo)系下可寫為： ，i=1,2,…,n (1-2) 對(duì)于互不影響的粒子系統(tǒng)，可取拉格朗日量： (1-3) 式中為第i個(gè)粒子的質(zhì)量，L為系統(tǒng)的動(dòng)能，即系統(tǒng)中所有原子的動(dòng)能之和。 如果考慮到原子間的相互作用，L可改記為: (1-4) 式(2-5)中等式右端的兩項(xiàng)分別代表系統(tǒng)的動(dòng)能和勢(shì)能。代入拉格朗日方程中可得系統(tǒng)的牛頓運(yùn)動(dòng)方程： ， (1-5) 式中即為原子i所受的內(nèi)力，即由于系統(tǒng)中其它原子的作用而在原子i上體現(xiàn)出的合力。 在分子動(dòng)力學(xué)模擬中，經(jīng)典拉格朗日方程常被用來計(jì)算原子的運(yùn)動(dòng)，更適合于求解系統(tǒng)運(yùn)動(dòng)的過程。在計(jì)算中能夠施加外部載荷如外力、約束、邊界條件等，在這些條件下通過求解原子的速度和位置，得到單個(gè)原子的運(yùn)動(dòng)軌跡，描述系統(tǒng)的運(yùn)動(dòng)過程，從而反映在系統(tǒng)整體特征上的變形、缺陷等。 2．哈密頓（Hamiltonian）運(yùn)動(dòng)方程 采用廣義坐標(biāo)和動(dòng)量的形式來描述粒子的運(yùn)動(dòng)，則可以得到哈密頓形式的運(yùn)動(dòng)方程，求解多粒子系統(tǒng)的狀態(tài)和演化過程。 含n個(gè)粒子的保守系統(tǒng)的系統(tǒng)哈密頓量在笛卡兒坐標(biāo)系中可以寫為： (1-6) 式中粒子的動(dòng)量，哈密頓正則方程為： (1-7) (1-8) 給定系統(tǒng)的初始狀態(tài)（初始位置和速度），則可以求解上述方程，得到系統(tǒng)中原子在任意時(shí)刻的位置和動(dòng)量，得到運(yùn)動(dòng)軌跡，由統(tǒng)計(jì)平均得到系統(tǒng)的熱力學(xué)表征。保守系統(tǒng)的哈密頓正則方程適合于求解系統(tǒng)的熱力學(xué)狀態(tài)，如溫度、熱流動(dòng)等。另外，在實(shí)際模擬中，哈密頓運(yùn)動(dòng)方程可以在保守系統(tǒng)的基礎(chǔ)上，通過改變系統(tǒng)狀態(tài)變量，與外界進(jìn)行能量交換，構(gòu)成不同系綜來進(jìn)行模擬計(jì)算。 （4）分子動(dòng)力學(xué)方法——程序流程 分子動(dòng)力學(xué)模擬過程與真實(shí)情況下的實(shí)驗(yàn)過程非常相似。在真實(shí)實(shí)驗(yàn)的進(jìn)行過程中，首先要做的工作就是準(zhǔn)備實(shí)驗(yàn)研究需要的材料，然后把實(shí)驗(yàn)儀器和測(cè)量?jī)x器與被研究材料連接起來進(jìn)行實(shí)驗(yàn)，在實(shí)驗(yàn)進(jìn)行過程的某一時(shí)刻或時(shí)間段，測(cè)量我們感興趣的物理量，進(jìn)來得到被研究材料的某些特性。與這一過程相似，在分子動(dòng)力學(xué)模擬中也遵循同樣的方法。一個(gè)典型的分子動(dòng)力學(xué)模擬的步驟如下： 1）讀入指定運(yùn)算條件的參數(shù)（如初始溫度、粒子數(shù)、密度、時(shí)間等）。 2）體系初始化（即選定初始位置和速度）。 3）計(jì)算作用于所有粒子上的力。 4）解牛頓運(yùn)動(dòng)方程。這一步和上一步構(gòu)成了模擬的核心。重復(fù)這兩步直 至我們計(jì)算體系的演化到指定的時(shí)間長(zhǎng)度。 5）中心循環(huán)完成之后，計(jì)算并打印測(cè)定量的平均值，模擬結(jié)束。 具體流程圖如圖1所示。 圖1 分子動(dòng)力學(xué)模擬流程圖 2.勢(shì)函數(shù) 在分子動(dòng)力學(xué)模擬計(jì)算中，原子間的相互作用由勢(shì)函數(shù)來具體描述。因此，勢(shì)函數(shù)的選取是分子動(dòng)力學(xué)模擬的核心問題，它的準(zhǔn)確性對(duì)計(jì)算結(jié)果的精度影響極大，直接決定了模擬結(jié)果正確與否。勢(shì)函數(shù)一般表示為粒子坐標(biāo)r的函數(shù)形式，它給出體系能量與粒子坐標(biāo)的關(guān)系。如果這個(gè)函數(shù)對(duì)體系的描述越準(zhǔn)確，那么用它來計(jì)算所得到的體系的性質(zhì)就與真實(shí)情況越接近。 （1）對(duì)勢(shì) 在分子動(dòng)力學(xué)模擬計(jì)算的早期，分子動(dòng)力學(xué)模擬所采用的勢(shì)函數(shù)一般是對(duì)勢(shì)模型( Pair potential)。對(duì)勢(shì)模型認(rèn)為原子之間的相互作用是兩兩之間的作用，并認(rèn)為系統(tǒng)能量為各粒子能量總和。對(duì)勢(shì)中一些參數(shù)主要通過宏觀實(shí)驗(yàn)參數(shù)來確定，其中有平衡點(diǎn)陣常數(shù)、彈性常數(shù)以及空位形成能、內(nèi)聚能和層錯(cuò)能等。它們?cè)谔囟ǖ膯栴}中均有各自的優(yōu)越性，可以比較好地描述除金屬和半導(dǎo)體以外的幾乎所有無機(jī)化合物。對(duì)勢(shì)雖然形式簡(jiǎn)單，并且不能很好的描述金屬及半導(dǎo)體的特性，但它的提出對(duì)勢(shì)函數(shù)的發(fā)展及分子動(dòng)力學(xué)的發(fā)展起到了不可替代的作用。 比較經(jīng)典的對(duì)勢(shì)模型有Lennard-Jones勢(shì)、Buckingham勢(shì)、Morse勢(shì)及附加勢(shì)等。對(duì)勢(shì)模型僅考慮原子對(duì)之間的作用力，無法考慮非金屬材料化學(xué)鍵的方位變化，鑒于此，人們開始考慮粒子間的多體作用，構(gòu)造出多體勢(shì)結(jié)構(gòu)，用它來精確描述更為復(fù)雜的系統(tǒng)。 （2）多體勢(shì) 多體勢(shì)于20 世紀(jì)80 年代初期開始出現(xiàn)，1984年，Daw和Baskes[14]首次對(duì)金屬提出了唯像近似的嵌入原子(embedded-atom method，EAM)。該理論的基本思想是：在一系列相互作用的原子中，一個(gè)原子對(duì)能量的貢獻(xiàn)是該原子所在處電子密度的函數(shù)。具體的操作方式是把晶體的總勢(shì)能分為兩部分：一部分是位于晶格點(diǎn)陣上的原子核之間的相互作用對(duì)勢(shì)，另一部分是原子核鑲嵌在電子云背景中的嵌入能。Baskes教授的領(lǐng)導(dǎo)下研究小組已經(jīng)對(duì)元素周期表中常用的42種元素測(cè)出了EAM所必須的物理參數(shù)，該項(xiàng)研究為原子多體勢(shì)理論的應(yīng)用做出了極大的貢獻(xiàn)。 在EAM模型中，系統(tǒng)總能量可表示為： (1-9) 其中F是嵌入能，第二項(xiàng)是對(duì)勢(shì)項(xiàng)。為系統(tǒng)中所有其它原子的核外電子在i處產(chǎn)生的局域背景電子密度之和，表示為： 具體到用于模擬計(jì)算某確定晶格模型的時(shí)候，上式中的各參數(shù)必須根據(jù)晶格模型具體包含的原子及原子間的可能結(jié)合方式來確定，在確定這些函數(shù)的具體形式和相關(guān)參數(shù)時(shí)，通常不再過分關(guān)注他們所代表的實(shí)際物理意義。實(shí)際應(yīng)用中在這些函數(shù)形式中將包含一些可調(diào)參數(shù)，通過調(diào)整這些參數(shù)，使所得模擬計(jì)算結(jié)果和實(shí)驗(yàn)結(jié)果盡可能的一致。 在Daw和Baskes公布了原型EAM理論不久，1984年，F(xiàn)innis和Sinclair[15] 基于束縛能帶理論發(fā)展了一種在數(shù)學(xué)上等同于EAM的勢(shì)函數(shù)。即將嵌入能函數(shù)設(shè)為平方根形式。Acldand[16]等在此基礎(chǔ)上通過擬合金屬的點(diǎn)陣常數(shù)、彈性常數(shù)、壓強(qiáng)體積關(guān)系、空位形成能、聚合能給出了適用于Cu，Al，Ni，Ag的多體勢(shì)函數(shù)。其中(4-7)式中的多體項(xiàng)及對(duì)勢(shì)項(xiàng)分別為 (1-10) 在上式中，當(dāng)時(shí)，；當(dāng)時(shí)，。,,,為常數(shù)，且有，。具體的取值因物質(zhì)不同而不同。 基于EAM勢(shì)的勢(shì)函數(shù)還有很多。由于這種勢(shì)函數(shù)很好地描述了金屬原子間的相互作用，因而成為描述金屬體系最常用的一種勢(shì)函數(shù)?？傮w來說，分子動(dòng)力學(xué)模擬中要求所采用的勢(shì)函數(shù)能夠盡可能詳細(xì)地、真實(shí)地描述原子間的作用力，但是目前MD 所采用的多體勢(shì)是用驗(yàn)的或半經(jīng)驗(yàn)的方法建立起來的，用其計(jì)算所得的理論數(shù)據(jù)與實(shí)驗(yàn)數(shù)據(jù)比較起來總會(huì)有一些差距，這在一定程度上影響和制約了分子動(dòng)力學(xué)的發(fā)展。因此，要使分子動(dòng)力學(xué)技術(shù)在材料性質(zhì)及規(guī)律的模擬方面有進(jìn)一步的應(yīng)用和發(fā)展，需要更新更優(yōu)良的數(shù)值計(jì)算方法，同時(shí)，更有效的模擬理論也亟待出現(xiàn)。 本文所采用的勢(shì)函數(shù)的具體形式及參數(shù)由Adams[17]等人給出,在Adams等人的工作中，對(duì)EAM勢(shì)參數(shù)進(jìn)行了擬合，擬合結(jié)果很好地滿足了鋁的晶格常數(shù)、結(jié)合能、體模量、彈性系數(shù)、聲子譜、熔點(diǎn)、熔化潛熱等，擬合結(jié)果表示在表1中，從表1中可以看出，該工作的擬合結(jié)果與已有的大量實(shí)驗(yàn)數(shù)據(jù)吻合的很好。隨后，在用分子動(dòng)力學(xué)方法模擬鋁的特性的工作中，該勢(shì)函數(shù)被廣泛地應(yīng)用。因此，在本文中采用這一勢(shì)函數(shù)來描述鋁在高壓下的相變特征。 二．計(jì)算過程 1．搭建樣品 分別利用分子動(dòng)力學(xué)程序搭建樣品: (1) 晶格常數(shù)為 4.032 Å，每個(gè)坐標(biāo)方向有20個(gè)原子的面心立方（fcc）結(jié)構(gòu)的樣品； （2）晶格常數(shù)為 3.147 Å， 每個(gè)坐標(biāo)方向有30個(gè)原子的體心立方（bcc）結(jié)構(gòu)的樣品； （3）晶格常數(shù) a=2.854 Å，c=4.652 Å，每個(gè)坐標(biāo)方向有20個(gè)原子的六角密排（hcp）結(jié)構(gòu)的樣品[18]。 表1 Adams擬合的勢(shì)函數(shù)計(jì)算的數(shù)據(jù)與已有實(shí)驗(yàn)數(shù)據(jù)的比較 2．用分子動(dòng)力學(xué)程序計(jì)算冷能曲線 利用分子動(dòng)力學(xué)程序，將三個(gè)樣品沿三個(gè)方向等比例壓縮，壓縮至原來體積的0.42倍，即壓縮至平衡晶格常數(shù)的0.75倍。通過修改分子動(dòng)力學(xué)程序，控制樣品壓縮過程，在壓縮過程中記錄樣品的能量，壓力等隨體積的變化，繪制冷能曲線，冷能曲線即單原子能量和單原子體積的關(guān)系曲線。 圖2 fcc結(jié)構(gòu)鋁、bcc結(jié)構(gòu)鋁、hcp結(jié)構(gòu)鋁的冷能曲線 圖3 fcc結(jié)構(gòu)鋁、bcc結(jié)構(gòu)鋁、hcp結(jié)構(gòu)鋁的冷能曲線 圖 2、圖 3 是本文計(jì)算的fcc結(jié)構(gòu)鋁、bcc結(jié)構(gòu)鋁、hcp結(jié)構(gòu)鋁在高壓下的冷能曲線。隨著壓強(qiáng)的增大，從圖 2 所表示的冷能曲線可以看到，隨壓力的增大，鋁在壓縮過程中會(huì)發(fā)生有fcc-hcp的相變。而從圖 3 所表示的冷能曲線可以看到，隨壓力的進(jìn)一步增大，在壓縮過程中鋁還會(huì)發(fā)生hcp-bcc的相變。 3.由冷能曲線判定鋁的高壓相變 從冷能曲線上可以判斷鋁在高壓下有無相變發(fā)生及發(fā)生相變時(shí)相變點(diǎn)所對(duì)應(yīng)的壓力。在相變點(diǎn)，以下兩個(gè)條件同時(shí)滿足： 1) 兩相的壓力相等， 2) 相變點(diǎn)的焓相等。 由條件1)得到 (2-1) 而通過條件2)得到 (2-2) 變換后得到 (2-3) (2-4) 因此在冷能曲線上，兩線切線的斜率就是所對(duì)應(yīng)的相變壓力點(diǎn)，因?yàn)榇藭r(shí)所得到的壓力正好滿足上述的兩個(gè)條件。 四．計(jì)算結(jié)果與討論 1．由圖 2可以看出，隨著壓強(qiáng)增大，原子體積的減小，fcc結(jié)構(gòu)鋁的冷能曲線與hcp結(jié)構(gòu)鋁的冷能曲線有一交點(diǎn)M?？梢?，在壓力加載體積減小的過程中，鋁會(huì)發(fā)生fcc-hcp的結(jié)構(gòu)相變。通過求兩冷能曲線的斜率可得到這一相變的相變臨界壓力為 61 Gpa。由圖 2 可知，隨著壓力進(jìn)一步升高，原子體積繼續(xù)減小，hcp結(jié)構(gòu)的冷能曲線和bcc結(jié)構(gòu)的冷能曲線出現(xiàn)交點(diǎn)N，在點(diǎn)N鋁發(fā)生fcc-hcp的結(jié)構(gòu)相變，在高壓下鋁的bcc結(jié)構(gòu)為穩(wěn)定相，由兩冷能曲線的斜率可得到這一轉(zhuǎn)變壓力為 210 Gpa。 2. 用分子動(dòng)力學(xué)模擬方法得到鋁fcc-hcp結(jié)構(gòu)的相變臨界壓力為 61 Gpa，由hcp-bcc結(jié)構(gòu)的相變臨界壓力為 210 Gpa。分子動(dòng)力學(xué)模擬很好地證明了鋁在高壓下會(huì)發(fā)生結(jié)構(gòu)相變，并計(jì)算出了相變臨界壓力分別為61 Gpa和210 Gpa，但由表 2 可以看出，鋁相變的第一性原理計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)明顯高于分子動(dòng)力學(xué)模擬得到的相變臨界壓力，例如：文獻(xiàn)[22]中利用LMTO-ASA的方法計(jì)算所得的fcc-hcp及hcp-bcc的相變臨界壓力分別為120 Gpa和200 Gpa，在文獻(xiàn)[19]于2010年的模擬結(jié)果中，fcc-hcp相變壓力為90 Gpa，同樣高于本文分子動(dòng)力學(xué)模擬結(jié)果。 Fcc hcp Hcp bcc Fcc bcc P V/V300 P V/V300 P V/V300 LCGTO-FF 205±20 0.510 565±60 0.364 340±15 0.436 GPT 360 560 LMTO-ASA 120 200 AP 220 0.50 380 0.40 300 0.45 LCGTO 330 0.446 Expt ＞219 ＜0.50 表2 鋁相變的第一性原理計(jì)算結(jié)果和實(shí)驗(yàn)數(shù)據(jù)[4,20,21] 導(dǎo)致本文分子動(dòng)力學(xué)模擬明顯低于理論計(jì)算結(jié)果的原因主要有兩點(diǎn)：（1）本文分子動(dòng)力學(xué)模擬結(jié)果是在溫度為0 K下的模擬結(jié)果，沒有考慮到溫度的影響，因此，相變壓力比實(shí)際相變壓力以及考慮到溫度影響的理論計(jì)算壓力要低的多。例如：J. C. Boettger 等人在計(jì)算中考慮到溫度的影響認(rèn)為300 K下鋁出現(xiàn)fcc-hcp相變的相變臨界壓力應(yīng)該在220 Gpa附近。（2）在本文模擬中是通過判斷在一定壓縮體積下哪一種相具有較低自由能來得到相變規(guī)律的，而實(shí)際相變需要一個(gè)形核長(zhǎng)大的過程，因此實(shí)際中出現(xiàn)相變的壓力點(diǎn)要比模擬計(jì)算所得到的結(jié)果高得多。Y.Akahama等人的研究中，在靜壓力加載下觀測(cè)到fcc-hcp結(jié)構(gòu)相變的臨界壓力為209 Gpa；文獻(xiàn)[19]于2010年的研究工作中，沿[001]方向單軸加載下，在90 Gpa附近觀測(cè)到fcc-hcp結(jié)構(gòu)相變。而在本文模擬結(jié)果中，fcc-hcp的相變臨界壓力為60 Gpa，可見相變過程中形核過程的存在使得實(shí)際相變臨界壓力滯后于本文計(jì)算的結(jié)果。 五．結(jié)論 本文采用分子動(dòng)力學(xué)模擬方法，用EAM勢(shì)函數(shù)描述鋁分子間的作用力，通過計(jì)算fcc、bcc和hcp結(jié)構(gòu)的鋁在壓縮過程中單原子能量和單原子體積間的關(guān)系，繪制出了鋁的冷能曲線。并利用冷能曲線判斷鋁在壓縮過程中的相變點(diǎn)，從而計(jì)算出了鋁在壓縮過程中的相變臨界壓力。計(jì)算結(jié)果表明：fcc-hcp結(jié)構(gòu)相變壓力為61 Gpa，hcp-bcc結(jié)構(gòu)相變壓力為210 Gpa。分子動(dòng)力學(xué)很好地證明了鋁在高壓下存在結(jié)構(gòu)相變。但由于本文的模擬中未考慮溫度及形核過程對(duì)相變壓力的影響，導(dǎo)致模擬計(jì)算結(jié)果低于理論計(jì)算和實(shí)驗(yàn)結(jié)果。 參考文獻(xiàn)： [1] 趙品，謝輔州，孫振國(guó)﹒材料科學(xué)基礎(chǔ)教程﹒第三版﹒哈爾濱：哈爾濱工業(yè)大學(xué)出版社，2007﹒170-178 [2] J. A. Moriarty and A. K. McMahan Phys. Rev. Lett (1982) 48, 809 [3] J. C. Boettger and S. B. Trickey Phys. Rev. B (1996) 29, 6434 [4] P. K. Lam and M. L .Cohen Phys. Rev. B(1983) 27, 5986 [5] V Mishra, S Chaturvedi Physica.B(2007) 393, 278 [6] R .G .Green, H. Luo and A L Ruoff Phys. Rev. Lett (1994) 73, 2075 [7] Y Akahama, M Nishimura, K Kinoshita, H Kawamura Phys. Rev. Lett, (2006) 96, 045505 [8] B Yaakobi, T R Boehly, D D Meyerhofer, et al Phys of plasmas, (2005) 12,092703 [9] S J Zhou, D M Beazey, P S Lomdahl,et al. Phys. Rev. Lett,(1996) 78, 01980 [10] B L Holian, P S Lomdahl, Science, (1998) 280, 285 [11] V V Zhakhovskii, V ZybinS, K Nishihara, et al. Suppl Prog Theor Phys（2000）138,223 [12] B J Alder, T E Wainwright J Chem Phys,（1957）27,1208 [13] A W Lees, S F Edwards. J Phys C: Solid St Phys, (1975) 5,1921 [14] M S Daw, M I Baskes, Phys. Rev. B, (1984) 29(12), 6443 [15] M W Finnis, J E Sinclair, Philosophic Magazine A, (1984) 50(1), 45 [16] G J Ackland, V Vitek, Physical Review B, (1990) 41(15),10324 [17] D L Adams, B Nielsen, J N Andersen, Phys.Scr, (1983) T4 22 [18] 黃昆，韓汝琦﹒固體物理學(xué)﹒北京：高等教育出版社，1985.1-6 [19] L Li, J L Shao, S Q Duan and J Q Liang New Journal of Physics (2010)12, 033011 [20] J. C. Boettger and S. B. Trickey Phys. Rev. B(1996) 6,53 [21] A. K. McMahan and J. A. Moriarty, Phys. Rev. B (1983)27, 3235