分子動力學模擬及其在微晶玻璃中的應用綜述

李楊井敏武吉偉劉立強

(1.山東建筑大學 材料科學與工程學院，山東 濟南250101；2.國家包裝產品質量監(jiān)督檢驗中心，山東 濟南250100)

0 引言

微晶玻璃作為一種新型無機材料，具備優(yōu)良的化學、力學、機械、光電性能，有著廣闊的應用前景，相關研究一直是關注熱點[1]。 傳統微晶玻璃的制備通常從改變原料配比和熱處理方式等方面開展。為得到一組可行有效的配方，通常需要反復的實驗，不僅消耗了大量的人力、物力，而且成效低并存在產品易開裂、成品率低、結構和性能關系尚不明細等問題[2]。 隨著微晶玻璃研究的深入，要求從微觀和介觀角度能夠準確、定量地描述微晶玻璃結構以及析晶過程規(guī)律，如研究微晶玻璃高溫熔體液體結構和性能的關系、堿金屬氧化物添加量對高溫熔體的黏度的影響、晶核劑的形核機理等。 分子動力學模擬作為凝聚態(tài)物理學常用的計算機模擬技術，從原子、分子層次認識物質的組成，利用計算機直觀展示和量化晶體結構，通過構建物質的微觀結構和數值模擬熱力學運動過程，得到每個粒子運動規(guī)律、能量波動等信息[3-4]。 微晶玻璃體系一般性模擬集中在三元或四元體系，完全可以借助基于分子動力學開發(fā)的模擬軟件開展工作，通過調整密度大小、粒子的數目、力場參數、施加的溫度、壓力，構建需要的微晶玻璃模型，從微觀層次上對微晶玻璃進行結構優(yōu)化，從而表征宏觀性能，對于改進微晶玻璃的性能從而指導生產優(yōu)良性能的微晶玻璃具有社會價值。

1 分子動力學模擬概述

1.1 分子動力學的發(fā)展歷史

20 世紀50 年代初，ALDER 等[5]用分子動力學方法設計了一個具有周期性邊界條件的粒子系統解決有關動量與壓力的問題，第一次真正意義上從微觀角度模擬了物質宏觀性能。 20 世紀70 至80 年代，ANDERSEN[6]對傳統的動力學方法加以約束條件，開創(chuàng)了恒壓條件下分子動力學方法。 GILLAN等[7]提出非平衡狀態(tài)動力學方法。 NOSE[8]提出恒溫條件下分子動力學方法。 CAR[9]解決了半導體和金屬勢函數難以模型化的問題。 CAGIN 等[10]提出了基于巨正則系統的分子動力學方法。 這些方法的提出都極大促進了分子動力學的的發(fā)展。 進入21世紀，隨著基礎學科基本理論的完善、計算機硬件的更新換代以及考慮更多作用力的多體勢函數的開發(fā)，使緩慢發(fā)展的分子動力學得到迅速發(fā)展，已廣泛應用于材料、醫(yī)藥、機械、化學、生物等多個學科，包含了晶體、非晶體、液態(tài)溶液、復合結構等方面，尤其對于材料極限條件(超臨界、深過冷、高溫、高壓)的實驗，分子動力學模擬優(yōu)越性顯著[11-16]。

1.2 分子動力學模擬的原理

第三種科學研究手段是除理論分析和實驗觀察之外的分子動力學方法，也被稱為“計算機實驗”手段[17]。 分子動力學(Molecular Dynamics，MD)方法和蒙特卡洛(Monte Carlo，MC)方法既有相同點又有所區(qū)別，最初MD 和MC 方法的產生的原因一樣，都是為了計算積分，不同點在于MD 方法在微正則系綜算積分，而MC 方法在正則系綜算積分；與MC 方法相比，MD 方法優(yōu)勢在于計算的內容和性質更多，得到的信息更多，而且MD 方法可使用的軟件要遠多于MC 方法，適應于大多數體系。 分子動力學模擬是借助MD 模擬軟件[18](如Moldy、Materials Studio、LAMMPS、DL_POLY 等)在原子及分子水平對多體系統進行求解的計算機模擬方法，在模擬過程遵循兩個假設[19]:(1) 體系中粒子運動遵循牛頓運動定律；(2) 粒子之間的相互作用滿足疊加定理。忽略量子效應和多體作用的模擬研究對象，微觀粒子的運動定義為體系原子核的運動，因此忽略了原子內電子的分布情況，這與實際存在的物理系統有一定的差別。

分子動力學模擬體系內的所有粒子的運動都遵循牛頓運動方程，即Fi(t)＝ miai(t) ，Fi(t)為粒子i在系統中受到的力；mi為粒子i的質量；ai(t)為粒子i的加速度。 勢能函數U對坐標ri求導，可得對于體系任何一個粒子，都由式(1)和(2)表示為

式中v為速度矢量，對式(1)和(2)求解，可以求得體系中任一粒子的位置和速度。

1.3 有限差分算法

對于多元體系，牛頓運動方程求解有一定的局限性，一般利用有限差分法完成對運動方程的求解，常用的方法包括Verlet 算法、Velocity-Verlet 算法、Beeman 算法、Gear 算法、Rahman 算法、蛙跳算法(Leap-frog algorithm)、Nordsieck 算法。

1.3.1 Verlet 算法

其算法計算過程較簡單，對計算機性能要求不高，每次積分只計算一次力而且時間可逆，因此運用最為廣泛，而其缺點是計算結果精度不高、軌跡與速度無關[20]。

1.3.2 Velocity-Verlet 算法

Velocity-Verlet 算法，其特點是可以同時計算出粒子一段時間內的位置、速度、加速度，并且準確性較高、計算量適中，已廣泛應用于分子模擬中。 其缺點是計算過程相較Verlet 算法復雜，用時較長[21]。

1.3.3 Beeman 算法

相較于Verlet 算法，Beeman 算法更為復雜，是一種更為精確的速度表達式，其動能是由速度計算直接得到的，并盡可能得到能量守恒結果，但是實際應用中對模擬條件要求較高，而且計算量較大[22]。

1.3.4 Gear 算法

Gear 算法由泰勒式展開，能夠預測每一個粒子新位置、速度、加速度[23]。

試驗地位于合作市卡加曼鄉(xiāng)新集村的甘南州農科所綜合試驗站，海拔2 721 m，年平均氣溫3.0℃，年降水量639.8 mm左右，無霜期93 d左右，耕種亞高山草甸草原土，旱川地，地力中等，前茬作物為油菜。播種時施尿素150 kg/hm2、磷酸二銨225 kg/hm2作基肥一次性施入，人工犁開溝溜種條播。2017年4月3日播種，6月6日中耕除草，田間管理略高于大田。

1.3.5 Rahman 算法

其算法表達式很復雜、計算量較大，但是計算結果精確，不適用于中小型體系模擬。 實際中很少采用Rahman 算法[24]。

1.3.6 蛙跳算法(Leap-frog algorithm)

蛙跳算法由Verlet 算法發(fā)展而來，相較Verlet式帶有t時刻的速度，不計算2y(t) 和y(t - δt) ，減少了計算量，提高了精確度，軌跡與速度無關。 其缺點是位置項和速度項不能同時計算[25]。 其表達式由式(3)表示為

1.3.7 Nordsieck 算法

Nordsieck 算法適用于解決大系統粒子牛頓方程的求解[26]。

1.4 勢函數

用來描述系統原子間或者分子間作用力的函數稱為勢函數[27]，也稱為力場。 勢函數不僅決定模擬能否順利進行，還將影響計算結果的精度。 精準且簡練的勢函數往往會使模擬效率提高，同時模擬結果符合物理規(guī)律。 20 世紀中后期，ALDER 等[5]借助數學方法，將粒子形象化硬球，一定距離之內會發(fā)生彈性碰撞，這種模型雖有缺陷，但是處理問題的思想促進了勢函數的發(fā)展，勢函數深刻體現了力場的關系。 在廣大科研人員共同推動下，出現了多種形式的勢函數，一般來說，根據相互作用力包含的粒子數分為二體勢、三體勢及多體勢，目前對勢、無方向性多體對泛函數、考慮角度效應多體勢是最常見的種類[28]。 常見的勢函數見表1。

表1 常見的勢函數表

1.5 系綜

2 微晶玻璃的分子動力學模擬綜述

2.1 分子動力學模擬對微觀網絡結構的分析

微晶玻璃的制備中存在高溫熔融階段，微晶玻璃高溫熔體的微觀網絡結構特征一直是模擬的重點問題。 需要結合均方位移函數(Mean Square Displacement， MSD)、 徑 向 分 布 函 數( Radial Distribution Function，RDF)、鍵長鍵角分布、配位數和橋氧數量等多個角度進行綜合分析，粒子在固定時間段內位移平方的平均值為均方位移[31]。 均方位移的量與原子的擴散系數存在著對應關系，均方位移數值越大，粒子的擴散系數越大。 對于研究微晶玻璃網絡結構，由MSD 可以推測擴散系數小的粒子很可能為網絡形成子或者網絡中間子，反之，則為網絡修飾子(不參與網絡結構的構建)。 RDF 描述的是某個原子周圍其他原子的分布特征，反映出給定某一粒子的位置坐標，在其半徑r的空間范圍內發(fā)現另一原子的概率[32]，用于表征材料的有序程度，利用這一特性，可以探究微晶玻璃體系粒子之間的相互結合能力大小。 鍵長、鍵角分布規(guī)律能反映出粒子鍵能強弱。 配位數是原子的第一近鄰原子的個數，橋氧定義為連接兩個網絡四面體且同時被兩個四面體所共用的氧粒子[33]。 橋氧數量的多少反映了硅氧網絡的完整性、致密性。 在微晶玻璃網絡結構中，四面體是構成網絡骨架的主要部分，與橋氧數量類似其也能反映網絡的致密性[34]。

2.2 分子動力學在微晶玻璃研究中的應用進展

分子動力學模擬由于獨特的實驗操作環(huán)境，在國內發(fā)展較快，受到越來越多的科研人員的青睞，為了更好地了解微晶玻璃微觀結構性質，指導特定性能微晶玻璃的生產，開展分子動力學模擬研究微晶玻璃的性能十分必要，但是目前在微晶玻璃領域的應用還較少。 ROSSANO 等[35]利用分子動力學模擬了CaO-FeO-2SiO2系玻璃，得到了鐵周圍氧原子的徑向分布情況，測量了四配位和五配位鐵原子的氧原子平均距離分別為1.99 和2.15 nm，這些基礎的有關原子距離的模擬為以后的諸多模擬正確性評價提供了一定的參考價值。 FUMIYA[36]選取NPT 系綜，采用Ewald 求和法模擬了Na2O·3SiO2熔體，發(fā)現Na—O 和Na—Na 的距離隨壓力的增加而縮短，Si—O—Si 鍵角隨壓力的增加而減小，說明了—Si—O—網絡的崩潰，表明了—Si—O—網絡的變形自由度增加，網絡在高壓下趨于崩潰，正由于這些結構的變化，Na2O·3SiO2熔體產生了致密化。 相較大多數模擬以單一組分對其他基礎玻璃成分的影響，研究壓力變化與微觀網絡結構的關系打開了另一種思路，基礎組分固然影響宏觀性能，但是某些特定溫度、壓力變化等外界因素也應被重視。 朱才鎮(zhèn)等[37-38]采用了一種二體和三體的多體相互作用勢探究了CaO-Al2O3-SiO2系微晶玻璃成分比例組成、微觀結構之間的關系，發(fā)現Ca/Al＝1/2 時，CaO-Al2O3-SiO2系微晶玻璃并不是傳統理論中認為的完整的網絡結構，而是存在一定的非橋氧，Ca/Al＜1/2 時，Si比Al 更容易形成網絡中間體。 通過模擬了不同Ca、Al 元素比情況下網絡結構的特征變化，發(fā)現了橋氧的存在，這是對傳統硅酸鹽網絡結構完整的觀點的一次否定，體現了計算模擬用來解決現實問題、提供理論創(chuàng)新的作用。 吳永全等[39-40]采用了BMH勢函數的基礎上，利用經典分子動力學模擬xCaO-(1-x)Al2O3(x為成分的變化)高溫熔體，重點探討了有關Al 的配位數及其網絡中的構成、微結構單元分布等結構性質，發(fā)現Al 的微觀結構單元是四面體的形式，類似硅酸鹽網絡中Si—O 鍵的結構特點，證明了Al 的配位數為4，Al 在網絡架構中起到了網絡形成子的作用。 可以看出，大多數的模擬是以基礎玻璃、三元體系為主，這是因為元素越多，需要考慮的作用力越多，開發(fā)的勢函數越復雜，因此分子動力學研究領域中，開發(fā)簡潔而且精確的勢函數是重要的方向。 王艷偉等[41]采用EAM 原子勢，在NPT 系綜條件下模擬了Ni 的含量變化對Ni-Zr 非晶玻璃微觀結構的影響，發(fā)現ico＋other(ico 為體系的十二面結構，other 為未知的配位結構)結構隨著Ni 數量的增加先增加后減少，ico＋other 結構反映了非晶的形成能力，當Ni 含量為50%時，結構數量最多，即說明了一定數量的Ni 有利于促進玻璃的形成。 趙亞賢等[42]利用經典分子動力學模擬對比了Na-Al-Si系微晶玻璃和高堿高鋁酸鹽玻璃中堿金屬的擴散行為，統計了兩種玻璃的Si—O、Al—O 鍵長，擴散系數，發(fā)現高堿高鋁酸鹽玻璃中堿金屬離子的擴散系數更大，這是因為[AlO4]四面體的存在，相較于[SiO4]四面體體積更大，會使網絡結構疏松，有利于堿離子擴散，證明了高含量的Al 為堿金屬離子的擴散提供了更多的擴散通道。 王海龍等[43]采用Mishin 嵌入原子勢模擬了金屬玻璃中的Cu 在急速冷條件下的形成過程，利用Verlet 算法求解運動方程，證明了隨著應變率的提高，Cu 的塑性流動應力增大，晶化程度增加，應力效應和溫度效應共同作用導致了金屬玻璃的晶化，對于解釋金屬的晶化起到了一定的促進作用。 肖成[44]構建了CaO-SiO2-Al2O3-Na2O四元體系的微觀玻璃模型，建立了1 873 K下體系中溫度和黏度的函數關系，發(fā)現Al2O3濃度高的體系傾向于形成更加復雜的網絡結構。 BINGHUI 等[45]巧妙設計了一種自上而下的模擬方法，探究影響Al2O3-SiO2系微晶玻璃斷裂韌性的因素，發(fā)現并證明了增強程度隨著納米晶體的尺寸、長徑比和排列角度的變化而變化。 王亞文[46]模擬了形核劑在熔融態(tài)高爐渣中的擴散行為，高爐渣成分與微晶玻璃成分相似，非均勻形核過程需要加入晶核劑降低析晶活化能。 通過數值與物理模擬發(fā)現形核劑的混勻時間與形核劑密度、等效直徑成正比，與熔體黏度成反比。 可以看出，有關微晶玻璃的模擬已經取得了一些進展，從鍵長、鍵角的分布，溫度、壓力等外界因素對網絡結構的影響，以及各元素在網絡中的位置、作用等都做出了一定的解釋，這些研究為制備性能優(yōu)異的微晶玻璃提供了一定的幫助。

井敏課題組在微晶玻璃的制備方面經過多年的積累，已有一套成熟制備微晶玻璃的工藝。 對于微晶玻璃微觀結構的探究，引進了先進的材料工作室(Materials Studio，MS)軟件開展工作，取得了一系列成果。 張國瑩等[47]基于MS 軟件利用分子動力學模擬了Na2O-Al2O3-SiO2微晶玻璃體系熔體結構，發(fā)現了作為網絡游離子的Na＋能夠促進Al3＋從網絡修飾子變成網絡形成子，且在一定范圍內，隨著Na＋的增加，Si—O—Si 和Al—O—Al 鍵角的分布范圍變小。 王健健等[48-49]探究了Fe3＋在CaO-Al2O3-SiO2系微晶玻璃中的行為，發(fā)現隨著Fe3＋含量的增加Al—O 鍵長分布更廣，峰型有單峰向雙峰轉化的趨勢。 王正[50]針對RO/R2O-Al2O3-SiO2熔體黏度特性進行了相關模擬工作，采用BMH 勢，選用Shear模塊計算了Na-Al-Si 系和Ca-Al-Si 系基礎玻璃的剪切黏度。 發(fā)現隨著堿金屬元素的減少，Na-Al-Si系黏度增加要大于Ca-Al-Si 系，整體結構上，Na-Al-Si 系相較Ca-Al-Si 系網絡結構較松散。 張永豪等[51]結合BMH 勢，建立了MgO-Al2O3-SiO2-TiO2系微晶玻璃的微觀結構模型，分別探究了TiO2晶核劑含量和SiO2含量對微觀結構的影響，發(fā)現隨著SiO2添加量的增大，體系中Mg2＋的均方位移逐漸減小，當SiO2的含量為57.1%時，橋氧總量最多，硅酸鹽網絡結構最完善。 TiO2含量增加會使Al—O 鍵長的分布更廣，峰形變寬。

3 展望

目前，采用分子動力學模擬技術對于微晶玻璃體系的模擬已經取得很大的進展。 從微觀角度出發(fā)，通過原子間相互作用力、邊界條件的選取、熱力學性質、運動方程的計算對微晶玻璃宏觀的性質(黏度、致密度、硬度等)做出了解釋，改進了微晶玻璃的性能。 但是仍存在一些問題，如初始模型過于理想化，而實際晶體材料存在各種缺陷，應當深入研究構建有缺陷的建模方法；MD 方法有觀測時間和系統大小的限制，仿真的規(guī)模較小，原子數較少，缺乏普遍性，應當發(fā)展較大規(guī)模的模擬；微晶玻璃通常使用經典分子動力學模擬，忽略了電子極化效應，對于電荷的相關信息無法獲取，勢函數過于依賴力場參數，數量較少，更加精確和簡潔的勢函數是發(fā)展分子動力學的重點；微晶玻璃模擬的模型也較為簡單，原子種類和數量較少，大多數模擬忽略了含量稀少的元素，或許會影響宏觀性能。 因此，多元體系的微晶玻璃的分子動力學模擬是發(fā)展的方向之一。