本文格式為Word版，下載可任意編輯——分子動力學與原子多體勢

一．分子動力學簡介

隨著納米科技的到來，大量新的學科產生了，例如納米電子學、納米生物學、納米材料學、納米機械學等。人們的注意力逐漸從宏觀物體轉向小尺度及相應的器件，其中微機械系統(mieromachine)或稱微型機電系統(mieroe一eetro一meeh耐ealsystem，MEMs)特別取得了成功，并正被拓展應用于各種工業過程。由圖1可知，分子動力學正是處于nm尺度下的研究方法。

圖1.不同模擬方法所對應的空間和時間尺度

1957年Alder和Wainwright[1]開創了分子動力學(Moleeularnynamies，MD)方法，之后經過多位科學家的努力，拓展了分子動力學方法的理論、技術及應用領域，特別是在20世紀80年代由Andersen等[2]先后完成的恒溫、定壓分子動力學方法，標志著分子動力學方法的科研應用進入了一個新階段。

分子動力學方法是研究納米尺度物理現象的重要手段。隨著越來越多的材料原子間作用勢函數被確切描述并經過試驗驗證、計算機硬件水平的快速更新以及高效率新算法的提出，分子動力學模擬被廣泛應用于納米尺度力學行為和納米材料力學性能的研究。

在納米尺度下，材料由離散的原子排列而成，由于比表面積大、表面效應明顯，材料的力學性能和力學行為將與宏觀材料迥異。基于連續性假設的宏觀連續介質理論在研究材料的損傷演化、失效過程時，往往在時間和空間上將原子尺度的缺陷進行平均化處理，但這種處理僅適用于大量缺陷分布在材料中計算區域的情形，而對大量微弱觀材料和力學試驗觀測到的現象都無法解釋，如疲乏與蠕變過程中的位錯模式、塑性變形的不均勻性、脆性斷裂的統計本質、尺寸效應等。因此，連續介質理論顯然難以確鑿求解納米尺度的力學問題。同時，假使直接從第一原理出發進行計算，除了類氫原子以外其他材料的薛定愕方程求解難度都太大，而且局域密度泛函近似理論并不是總能滿足實際問題的需要。另一方面，材料本身在空間、時間和能量等方面存在藕合和脫禍現象[3,4]，直接從頭開始的量子力學計算難以很好地應用到幾百個原子以下的計算規模中，無法達到一般納米材料和器件的模擬要求。此外，由于試驗條件控制的困難和合成、制備方式不同，各種納米材料力學性能的有關試驗結果分散性較大甚至相反，以至于目前難以通過納米力學試驗得到普遍適用的定量力學規律。

鑒于理論和試驗上的困難，通過分子動力學方法模擬納米尺度的力學性能和行為來摸索納米尺度的一般規律，是進行納米力學研究的有效方法。分子動力學最早用于熱動力學和物理化學，計算不同物理系統如固體、氣體、液體的整體或平均熱化學性能。1957年Alder首次提出并采用分子動力學方法分析剛性球系統的固液相變問題取得成功，此后，分子動力

1

學開始逐漸應用于材料領域。隨著上世紀80年代計算機硬件水平的提高和各種描述原子間作用的勢函數的提出，分子動力學模擬日益活躍。通過分子動力學模擬不僅可以深入了解繁雜的機制，發現本質上嶄新的現象，而且可以定量模擬真實固體中所發生的過程，是諸如表面結構和擴散中的動力學和穩定性等大量結果的唯一來源[5-9]。在EAM理論逐漸成熟和Baskes試驗室[10-13]、Ackland試驗室[14-18]等確切測出大量常用材料的EAM參數以后，分子動力學方法在模擬材料的物理性能和現象方面逐漸顯示出強大的計算能力和較高的精度，大量的模擬特別是固體結構、位錯運動、表面界面現象、力學性能、變形機制和流體中的電泳、電滲透流等都得到了理想的結果。只要能將基于物理的模型建立起來，通過分子動力學計算就可以透露出物理現象的本質，逐漸被廣泛應用于凝聚態物理、材料學、力學、生物學、微電子學和微納米加工等領域。

晶體是由大量的原子有序排列而成，材料的強度來源于原子間的相互作用，塑性來源于原子間的相互運動。因此，直接從原子尺度上對材料的微觀力學行為進行研究顯得十分有必要。分子動力學模擬技術既能得到原子的運動軌跡，還能像做試驗一樣觀測。對于平衡系統，可以在一個分子動力學觀測時間(ObservationTime)內做時間平均來計算一個物理量的統計平均值，對于一個非平衡系統過程，只要發生在一個分子動力學觀測時間內(一般為1一10ps)的物理現象也可以用分子動力學計算進行直接模擬。特別是大量與原子有關的微觀細節，在實際試驗中無法獲得，而在計算機模擬中可以便利的得到。這種優點使得分子動力學方法廣泛運用于材料科學與工程中，如材料設計、斷裂分析等。

近年來，利用計算機模擬技術研究材料的力學性能日益成為人們感興趣的課題。由于計算機處理速度的迅速提高，計算機模擬已經和試驗觀測、理論分析并列成為本世紀科學研究的三種方法[19]。計算機模擬的數據(從模型中得來)可以用來比較、驗證各種近似理論;同時，計算機模擬方法還可用來對試驗和模型進行比較，從而提供了評估一個模型正確與否的手段。計算機模擬方法還有一個優點就是可以溝通理論和試驗。某些量或行為可能是無法或難以在試驗中測量的，而用計算機模擬方法，這些量可以被確切的計算出來[20]。分子動力學模擬方法更以其建模簡單、模擬結果確鑿的特征而倍受研究者的關注。

分子動力學模擬(MolecularDynamicsSimulation)用于計算以固體、液體、氣體為模型的單個分子的運動狀況，是一種聯系微觀世界與宏觀世界的強有力的計算機模擬方法側。

二．分子動力學基本原理

分子動力學方法是納米計算力學的主要手段。它對多經典粒子系統的運動方程組進行數值積分，得到相空間軌道，進而研究系統的平衡熱力學性質、結構動力學性質和非平衡傳輸性質等。在分子動力學模擬中，將原子視為質點或忽略內部結構的固體球。首先建立粒子系統的幾何模型，通過描述原子間相互作用的勢函數，根據給定的邊界條件、初始條件求出系統中每一時刻單個粒子或原子的能量和所受到的力，代入牛頓動力學方程組求解原子的位置和速度，得到系統在相空間的運動軌跡。對足夠大數量的粒子在足夠長時間的結果進行統平均，則可以得到類似于宏觀意義上的物理量和力學量。

早期研究孤立系統的保守系綜分子動力學基于兩個基本假設[21]:

a.粒子是相互作用的質點，運動由位置矢量和速度矢量來描述。粒子間的相互作用取決于彼此的空間位置。

b.系統無質量交換。即模擬過程中系統的原子數不變。

保守系綜分子動力學常用于模擬能量守恒的孤立絕熱系統。但是要確鑿模擬實際納米材料和器件的表面、界面等現象，真實反映粒子系統受外界約束狀況下的物理行為，以及進行跨尺度的禍合模擬，保守系綜分子動力學模擬明顯是遠遠不夠的。在此基礎上，研究者提出了耗散系統與周邊介質進行能量交換的不同理論和算法，對粒子系統進行溫度、壓力、粒子數、體積等不同控制，適合于不同系綜的分子動力學模擬。

2

分子動力學模擬中假設系統所有粒子的運動遵循經典牛頓運動定律，且忽略原子背景電子云的量子效應，原子間的相互作用滿足疊加原理。可見，分子動力學方法是一種對廣義牛頓運動方程的近似的數值積分方法。

1.下面分別描述分子動力學方法的基本方程。

(1)拉格朗日（Lagrangian）運動方程

對于N自由度的由相互作用的質點構成的運動系統，拉格朗日方程為：

?為質點位置對時間的導數。在笛式中q為廣義坐標，指定質點的空間位置。q卡兒坐標系下，由N個原子構成的模擬系統的拉格朗日方程可寫為：

式中

子系統，可取拉格朗日量：

，為原子i的位置矢量，在三維空間n=N/3。對于互不影響的粒

式中mi，為粒子i的質量，L即為系統的動能，即系統所有原子的動能之和。假使考慮

到原子間的相互作用，L可改一記為

式右端的兩項分別表示系統的動能和勢能。代入拉格朗日方程可得系統牛頓運動方程:

式中Fi即為原了i所受的內力，即由于系統中其他原子的作用而在原子i上表達出的合力。由牛頓運動方程建立線性微分方程組，給定初始條件(初始位置、初始速度)，求解該封閉方程，可得到確定解，求出原子運動的軌跡，即任意時刻原子的位置ri(t)和速度ri(t)。分子動力學模擬中的經典拉格朗日方程常用來計算原子的運動，得到單個原子的運動軌跡，描述原子系統的運動過程，以及反映在原子系統整體特征下的位移、變形、缺陷等。(2)哈密頓(Hamiltonian)運動方程：

假使采用廣義坐標系和動量的形式來描述粒子系統的運動，則可以得到哈密頓形式的運動方程，求出多粒子系統的狀態和演化過程。

在笛卡兒坐標系中，對含n個粒子的保守系統，取系統哈密頓量為：

3

.

式中Pi，為粒子i的動量，Pi=mii，哈密頓正則方程為:

假使給定系統的初始狀態(初始位置和速度)，則可以求解線性微分方程組，得到系統原子的運動軌跡，即任意時刻的原子位置和動量，由統計平均得到系統的熱力學表征。

在分子動力學模擬的基本方程中，拉格朗日運動方程更適合于求解原子系統運動的過程(求解原子速度和位置)，并能施加外部荷載如外力、約束、邊界條件等，而保守統的哈密頓正則方程更適合于求解系統的動力學演化過程和熱力學狀態，如溫度、熱滾動。在實際的分子動力學模擬中，拉格朗日方程尋常與外部約束和邊界條件一起，構成特定原子系統模型的初值問題;而Hamiltonian運動方程往往在保守系統的基礎上，通過改變系統狀態變量，與外部環境進行能量交換，構成不同系綜的分子動力學模擬。

2.分子動力學原子間作用勢函數分子動力學方法是通過原子間的相互作用勢，依照經典牛頓運動定律求出原子運動軌跡及其演化過程。分子動力學計算的關鍵是原子間勢函數的選取，它決定著計算的工作量以及計算模型與真實系統的近似程度，直接影響到模擬結果的成功與否。由于物質系統的繁雜性以及原子間相互作用類型的不同，很難得到滿足各種不同體系和物質的一般性而又精度較高的勢函數。所以針對不同的物質體系人們不斷發展了大量的經驗和半經驗的勢函數[22]。

從分子動力學模擬的基本方程可以看出，分子動力學模擬即求解偏微分方程，求解的精度關鍵在于勢函數U(r1,r2,?rn)描述原子間相互作用的確鑿性。

多原子系統的勢函數可以表示成:

r4

式中Um為m體勢，即原子勢能受m階多體效應的影響，其中U1為系統原子受外力場(如重力)的影響項。在尋常的計算中，由于系統內原子受原子間相互作用的影響遠大于單個原子本身受外力場的影響，同時為了減少計算量，一般忽略外力場的影響U1和三階以上的多體效應Um(m>3)，這就是應用于分子動力學模擬的對勢模型(PairPotential)。對不同的固體材料，將其原子三階以上多體效應的影響作為修正項記入對勢模型的二階作用項中，就形成了各種不同的多體勢函數。（1）對偶式

對偶勢理論認為任何兩個原子結合鍵的強度不由于它周邊其它原子的存在而受到影響。1.Lennard一Jones(L一J)勢函數

在對勢模型的典型代表、應用廣泛的Lennard一Jones模型[23-24]中，取原子間作用的勢函數方程為:

式中σ—原子與零點勢的距離;ε—最小勢能處的能量;

m,n—調整系數，一般(m,n)取值為(12,6),(10,5),C8,4)

方程右端的第一項描述原子間的排斥作用，其次項描述原子間的吸引作用。|ry|=|ri-rj|，即原子間的距離。當原了間距離rij=σ時，原子勢能Ф(rij)=0;rij>σ時，Ф(rij)0。ε表征原子間的吸引/排斥強度，對原子中的一個從平衡位置r0移到無窮遠處所作的功，對于不同物質的原子,ε各不一致，尋常為10-19J或ev量級。能量最低時刻的r值，即r0，表示原子間的平衡位置，此時原子間作用力為零。對上述方程求極值，可得r0=62??1.22σ

將上式代入式模型中原了間作用力為：

，可得Lconard一JoneS

2.Morse勢函數

1929年，Morse注意到雙原子分子的振動譜的量子力學問題可用指數形式的勢函數解析地解決，并發現計算結果