基于分子動力學(xué)原理的性質(zhì)計(jì)算模擬

胡永波 聶德志

【摘 要】使用SCIGRESS軟件，應(yīng)用分子動力學(xué)方法，選擇NPT系綜，模擬從2200K降溫到1000K的過程中各參數(shù)的變化情況，得到了的勢能、密度、雙體分布函數(shù)、均方位移、多面體和多邊形隨溫度變化的計(jì)算結(jié)果。

【關(guān)鍵詞】SCIGRESS;分子動力學(xué);雙體分布函數(shù);擴(kuò)散系數(shù)

一、分子動力學(xué)原理

分子動力學(xué)是在原子分子水平上求解多體問題的重要計(jì)算機(jī)模擬方法，可以預(yù)測納米尺度上的材料動力學(xué)特征。通過求解所有粒子的運(yùn)動方程，分子動力學(xué)方法可以用于模擬與原子運(yùn)動路徑相關(guān)的基本過程，在分子動力學(xué)中，粒子的行為是通過經(jīng)典的牛頓運(yùn)動方程，所描述的分子動力學(xué)是確定性方法，一旦出使構(gòu)型和速度確定了，分子隨時(shí)間產(chǎn)生的運(yùn)動軌跡也就確定了。分子動力學(xué)的標(biāo)準(zhǔn)流程如下：

1.Initialize——初始化粒子的空間位置和速度

Boltemann能量均分定理：對于處在溫度為T的平衡狀態(tài)的經(jīng)典系統(tǒng)，粒子能量中每一個(gè)平方項(xiàng)的平均值等于kBT/2，下列公式表達(dá)了微觀系統(tǒng)總動能和宏觀溫度的關(guān)系：

2.Forces——計(jì)算粒子的受力

在保守力場中，保守力等于勢能函數(shù)的梯度，即：

其中，U1表示外力場或邊界條件（如容器壁）;U2表示兩體作用勢，一對原子之間的作用力，不包含其他原子的影響;U3表示三體作用勢，描述由于第三個(gè)原子的存在對原子對之間作用力的影響。

在具體問題當(dāng)中，需要選取合適的勢函數(shù)以及周期性邊界條件來通過有限的模擬區(qū)間，得到具有宏觀尺寸的結(jié)構(gòu)材料的特性。

3.Motion——使用牛頓運(yùn)動方程得到系統(tǒng)隨時(shí)間的演進(jìn)

在該情形下，需假設(shè)粒子受力不變且粒子運(yùn)動符合牛頓運(yùn)動方程。根據(jù)Verlet算法原理，粒子位置的Taylor展開式如下：

據(jù)此求解運(yùn)動方程，求解得到系統(tǒng)的粒子位置、速度、加速度隨時(shí)間的演進(jìn)。這一步驟的實(shí)現(xiàn)可以依賴于多種成熟的算法。除此之外，還可使用Leap-frog算法、Velocity Verlet算法等。

4.Analysis——計(jì)算系統(tǒng)的密度、溫度、勢能等參數(shù)

在本次實(shí)驗(yàn)中，除了勢能、溫度之外，我們還得到了雙體分布函數(shù)及多面體、多邊形的結(jié)果。

5.Summarize——對計(jì)算數(shù)據(jù)進(jìn)行整理分析

根據(jù)得到的數(shù)據(jù)對演化過程的物理學(xué)意義進(jìn)行分析。

二、模擬計(jì)算步驟

1.選取區(qū)域大小（5*5*5）、晶體結(jié)構(gòu)、NPT系綜、步長、溫度梯度（2200-1000k），確定勢函數(shù)

2.設(shè)定初始條件，原子分子初始速度服從波爾茲曼隨機(jī)分布

3.求解運(yùn)動方程，等待系統(tǒng)達(dá)到平衡態(tài)，輸出穩(wěn)定的參數(shù)值

4.用時(shí)間平均代替系綜平均，對宏觀物理量進(jìn)行統(tǒng)計(jì)

三、計(jì)算結(jié)果和分析

3.1 體系的T溫度、D密度、P勢能隨步長的變化情況

定性分析程序的運(yùn)行結(jié)果，該體系在溫度從2200K達(dá)到1000K的過程中，密度隨時(shí)間不斷增大。勢能隨時(shí)間逐漸減小。說明體系性質(zhì)正在由液體逐步趨近于固體。

密度是金屬的一個(gè)重要的熱物理性質(zhì)，對于理解液-固相轉(zhuǎn)變、凝固過程、傳質(zhì)、熱對流，探索液態(tài)結(jié)構(gòu)和亞穩(wěn)態(tài)相變等有一定意義。為了更加清晰的表示密度的變化情況進(jìn)行了縮放操作，實(shí)驗(yàn)所獲得的密度逐漸增長。

3.2 、均方位移變化趨勢

程序運(yùn)行得出鎳和鋯的均方位移變化趨勢如上圖所示。通過分段讀數(shù)，顯示Plot Information，讀到擴(kuò)散系數(shù)：

以上表格所列出的是和在不同溫度下的擴(kuò)散系數(shù)以及整個(gè)系統(tǒng)演變過程中的擴(kuò)散系數(shù)。結(jié)果表明擴(kuò)散系數(shù)逐步降低。這是由于溫度的降低減緩了擴(kuò)散的發(fā)生。

3.3計(jì)算所得的雙體分布函數(shù)

實(shí)驗(yàn)得到了計(jì)算所得的2200K，2000K，1800K，1600K，1400K，1200K下、、和總體系的雙體分布函數(shù)，以1600K的為例，結(jié)果如下：

根據(jù)以上數(shù)據(jù)，繪制、和總體系的雙體分布函數(shù)第一峰的峰值變化圖。結(jié)果表明，隨著溫度的降低雙體分布函數(shù)第一峰、第二峰的峰值在不斷上升。

雙體分布函數(shù)顯示，液態(tài)金屬的結(jié)構(gòu)介于固體和液體之間，兼具短程有序和長程無序的特征。原子可以在確定的范圍內(nèi)運(yùn)動，因此雙體分布函數(shù)在短距離有類似于固體的峰，峰的寬化由原子運(yùn)動引起。在長距離時(shí)，雙體分布函數(shù)的值趨近于1，表明在長程范圍內(nèi)液態(tài)金屬具有氣體的特征，原子的分布完全無序。

3.4多面體、多邊形結(jié)果

選取了三個(gè)代表點(diǎn)，展示了演化初期（液態(tài)）、中期、后期（固態(tài)）的多面體結(jié)果和多邊形結(jié)果：

結(jié)果表明，大多數(shù)多面體的數(shù)量隨時(shí)間未發(fā)生明顯變化，主要是四面體和六面體的數(shù)量變化較大。其中，四面體逐漸增多，六面體先減少，后增多。變化原因主要是固體和液體的微觀結(jié)構(gòu)不同。

大多數(shù)種類多邊形的數(shù)量隨時(shí)間未發(fā)生明顯變化，存在量最大的是四邊形、三邊形次之，接下來是五邊形，最少的始終是六邊形。變化原因主要是固體和液體的微觀結(jié)構(gòu)不同。

四、總結(jié)

隨著降溫過程的進(jìn)行，液態(tài)金屬由液體變?yōu)楣腆w。溫度越低，體積越小，密度越大，擴(kuò)散系數(shù)也逐步降低。雙體分布函數(shù)計(jì)算結(jié)果顯示，第一峰的峰值不斷上升。最后計(jì)算并分析了液態(tài)金屬的多面體、多邊形結(jié)果，體系中多面體和多邊形請款隨溫度無明顯變化。

參考文獻(xiàn)：

[1]王海鵬.液態(tài)金屬結(jié)構(gòu)與性質(zhì).西北工業(yè)大學(xué).2020.8.

[2]劉金遠(yuǎn)，段萍，鄂鵬.計(jì)算物理學(xué).科學(xué)出版社.2016.

[3]武玉琴.金屬凝固過程中過渡態(tài)現(xiàn)象研究.山東大學(xué).2008.

作者簡介：

胡永波，男，2000，陜西渭南，主要研究材料及物理，西北工業(yè)大學(xué)物理科學(xué)與技術(shù)學(xué)院，西安，中國陜西省西安市長安區(qū)東祥路一號西北工業(yè)大學(xué)長安校區(qū)。

聶德志，男，2000，江蘇鹽城，主要研究能量采集及微機(jī)電系統(tǒng)，西安，中國西北工業(yè)大學(xué)長安校區(qū)。

（作者單位：西北工業(yè)大學(xué)）