淺析軟膠體粒子形成束晶的動力學模擬范文

本站小編為你精心準備了淺析軟膠體粒子形成束晶的動力學模擬參考范文，愿這些范文能點燃您思維的火花，激發您的寫作靈感。歡迎深入閱讀并收藏。

摘要：利用廣義指數模型描述軟膠體粒子，結合分子動力學模擬研究軟膠體粒子形成束晶的動力學過程．通過等溫壓縮和等密度降溫2個不同的過程，研究了束晶形成過程中結構變化特征和動力學路徑對結構的影響規律．研究發現，與蒙特卡洛模擬結果相比，分子動力學模擬得到的結構隨著密度的變化有明顯的遲滯現象，這是由于考慮了真實的動力學因素引起的差異．此外，在相同溫度和壓力下通過不同的動力學路徑得到的相結構不完全相同，這是由于動力學形成過程會對相結構產生很大的影響．

關鍵詞：軟粒子束晶;分子動力學模擬;廣義指數模型;面心立方;體心立方

軟膠體粒子體系具有豐富而復雜的相行為，在外部條件改變時，軟膠體粒子可處于液態、固態或液固共存態．通過改變溫度和密度，軟膠體粒子不僅具有從液態到固態的一階相變，還具有重入熔融和同態結構轉變等相行為［1］．超軟的膠體粒子還可以形成束晶(Clustercrystal)，這是一種由多個膠體粒子占據同一晶格的晶體結構［2～7］．在束晶的1個粒子束中，膠體粒子的質心位置并沒有完全重疊，即使膠體粒子表現出非同尋常的軟化性質，其質心重疊度仍然是有限度的［8］．在實際生活中也存在許多類似的束晶結構，比如由顆粒團簇形成的面心立方(fcc)或體心立方(bcc)晶體［9～11］．由于束晶的獨特性質，其一經發現就引起了軟物質領域的廣泛研究興趣［12～17］．許多研究者利用軟勢能形式(如高斯勢、簡諧勢和赫茲勢等)來描述軟粒子，并用計算機模擬手段研究了其相行為［4，18～21］．Miyazaki等［20］利用分子動力學模擬研究了廣義赫茲勢相互作用的超軟粒子在低溫和不同密度下形成不同粒子束結構，以及在流體相中發生的相轉變．張凱等［7］利用蒙特卡洛模擬繪制了廣義指數模型粒子的相圖．目前的研究表明，軟粒子豐富的相行為與其勢能形式息息相關．當粒子質心重合時對于一個相互作用勢能達到最大值的軟勢，如果傅里葉變換后只有實部，體系隨著密度的增高，會出現重入熔融現象;如果傅里葉變換后有虛部，所有的溫度下都沒有重入熔融現象，而在高密度下會形成束晶［4］．目前的研究多集中在束晶的結構和相行為等方面，對于束晶形成的動力學過程，以及調控動力學過程對結構的影響規律還不是很清楚．本文利用廣義指數模型(GEM)，研究束晶形成和結構轉變的動力學過程．廣義指數模型定義為φ(r)=εexp［－(r/σ)n］，n＞2，ε和σ分別為能量和長度［2，6，7，22］．本文利用分子動力學模擬，對GEM(n=4)模型粒子體系進行等溫壓縮和等壓降溫2個不同熱力學路徑的退火，研究束晶形成過程中粒子束和晶體結構出現的先后順序;通過與蒙特卡羅模擬結果相比，明晰動力學過程對束晶結構的影響．以期揭示束晶的形成動力學過程，為軟粒子膠體的結構調控提供幫助．

1模型及研究方法

分子動力學(MD)模擬是原子、分子水平上的計算機模擬方法，該方法基于牛頓力學，通過求解運動方程來精確得到原子或分子在不同時刻的位置和速度，進而獲得計算體系熱力學量和其它宏觀量［23，24］．假設1個含有N個粒子的體系，粒子之間相互作用勢是位置的函數U(r1，r2，…，rN)，根據經典力學可知體系中任一粒子所受的力為勢能的梯度:求解運動方程獲得粒子不同時刻的速度和位置．分子動力學模擬過程如下:首先由體系粒子的位置計算體系的作用勢，再由式(1)和式(2)計算各粒子的加速度，然后選取1個小的δt并采取合適的算法解運動方程，獲得經過δt后各粒子的位置和速度．重復以上步驟獲得經過2δt后粒子的位置和速度，如此反復，可得到各時間下體系中粒子的運動軌跡及速度信息，由此獲得相關表征量．采用分子動力學模擬方法，在NVT系綜下，用速度Verlet法積分牛頓運動方程．溫度由Nosé-Hoover熱浴控制［24］，模擬中采用的時間步長為0.001．為了獲得可靠的結果，每106時間步記錄一次體系的物理性質．分子動力學模擬使用大尺度的GPU加速的GALAMOST軟件包［25］．該體系包含4096個膠體粒子，使用GEM(n=4)模型描述膠體粒子之間的相互作用．約化溫度定義為T=kBT/ε，ε為能量的基本單位．模型中σ=1.0，ε=1.0，粒子之間相互作用的截斷半徑rc=2.0．在高斯勢下，通過改變體系的密度和溫度考察膠體粒子形成束晶的動力學過程．等溫壓縮和等密度降溫的熱力學路徑由基于蒙特卡羅模擬繪制的GEM-4的相圖(圖1)中給出．為了對模擬結果進行有效分析，計算了體系的有序度和排列結構．傳統意義上的共近鄰方法(CAN)對于分析fcc、bcc和密排六方(HCP)等結構是非常有效的手段．該方法通過設定截斷半徑能清楚地提供局部粒子排列的信息，從不同晶型粒子體系中分辨出結晶缺陷．但是，在多相體系中，由于全局定義的截斷半徑的選擇性不是特別明確，使得這類傳統CAN方法并不適用于多相體系．因此，在本文模擬體系中，采用OVITO軟件自帶的自適應方法(a-CAN)分析體系粒子排列情況．自適應a-CNA方法是傳統CNA方法的衍生物，主要用于分析多相體系．通過事先指定體系中每一個粒子的截斷半徑，并通過與參照結構比較來自動調整每個粒子的截斷半徑，進而分析體系局部粒子的排列信息［26］．利用鍵取向序參量(BOO)表征粒子與其近鄰之間相對位置的取向特性，BOO也被稱為Steinhardt有序參數，通常用于確定分子模擬中的晶體結構．尤其是Q4和Q6常被用作區分立方和六方結構的最優選擇［27］．

2結果與討論

2．1等溫壓縮過程模擬將模擬盒子的大小設定為27×27×27，盒子中包含的粒子數N=4096．在體系溫度恒定為0.01時，對盒子進行壓縮，即增大體系的密度，觀察膠體粒子在超軟勢作用下形成結構及結構演化的動力學過程．圖2是最終形成的結構．隨著密度的增大，膠體粒子的結構由最初的無序液態逐漸變得有序并且進一步形成規則的晶體．經結構分析發現，這是由單粒子組成的面心立方(fcc1)晶體結構．繼續增加密度，單粒子晶體逐漸轉變成二聚體的晶體．當密度增加到1.0左右時，體系完全演化成二聚體的面心立方(fcc2)晶體結構．該結果與文獻［7，18］的研究結果非常吻合．如圖3所示，此動力學過程中束晶的類型隨著密度增大而改變．尤其是當ρ＞0.9時，軟粒子束晶的聚集狀態發生了突變，單粒子數減少，二聚體數目增加．ρ=0.94是二聚體和單粒子數量之間轉換的臨界值，當ρ＞0.94時二聚體的數目超過單粒子的個數，體系以二聚體為主體．隨著ρ增加，Q4和Q6值發生變化(Q4和Q6判斷晶體結構的依據參考文獻［27］)．在ρ=0.57～0.70范圍內，體系屬于fcc晶體．從ρ=0.70增加到ρ=0.95，體系經歷了兩相共存區．在這一過程中fcc晶體經歷了1個不連續的同構相變［28～30］．在ρ=0.7～1.0之間，聚集狀態由單粒子轉換成二聚體，且fcc晶體的同構轉變為部分bcc晶體類型．當ρ=1時，體系主體為二聚體的fcc2晶型．與蒙特卡洛模擬得到的結果相比，分子動力學過程中發生的一系列結構的轉變具有明顯的遲滯現象．這是由于動力學因素對體系的影響［31］．由于在共存相及其附近成核的自由能壘較高，因此在較短的動力學模擬時間內很難觀察到相轉變現象，因此導致了遲滯效應［24］．當溫度恒定為0.01時，對體系進行等溫壓縮，使ρ從0.1變化到1.5，體系結構隨著密度的變化過程如圖4所示．紅色粒子是被標記的粒子，利用示蹤粒子的方法，給出不同密度下束晶形成過程中粒子聚集狀態以及體系的結構演化過程(圖5)，在此動力學過程中，體系密度在ρ=1.1～1.45之間，單粒子、二聚體以及三聚體的數量均發生了轉變．密度在ρ=0.7～0.85之間結構屬于fcc晶體．在ρ=0.85～1.4之間，Q4和Q6出現了2次降低，分別對應聚集體類型的變化．此時的晶體結構主要為fcc和bcc兩相共存區，且fcc晶體逐漸向bcc晶體的結構轉變．從圖5可以看出，體系的晶體最終不能完全轉化成三聚體的面心立方(fcc3)結構，而是形成了二聚體的體心立方(bcc2)晶體結構．這一結果與隨機蒙特卡羅模擬方法得到的結果明顯不同．這可能是由于fcc2與fcc3結構之間的自由能壘較高，因而采用動力學方法難以獲得熱力學平衡態結構．此外，還觀察到三聚體的形成有3種過程:(1)由3個單粒子聚集形成;(2)由1個二聚體和另一個二聚體融化后分離出的單粒子組成;(3)由1個單粒子和1個二聚體聚集形成．經數據分析發現，第3種的形成方式占主導。

2．2等密度降溫過程模擬為了全面理解膠體粒子在高斯軟勢中的重入熔融以及同構轉變的現象，控制模擬體系的密度為ρ=1，對體系進行降溫處理(T=1.0→0.01)，觀察模擬體系形成的晶型結構，結果如圖6所示．隨著溫度的降低，單粒子數目迅速降低，二聚體數目迅速增加，且在T≈0.8時二聚體數目大于單粒子數目，表明此時體系中粒子的主要存在狀態為二聚體．當T=0.65時，Q4和Q6值迅速增加，且體系由無序的聚集體轉變成有序的fcc晶型．當溫度進一步降低到T≈0.08時，Q4和Q6略微降低，表明體系中的晶體有序度下降．結合OVITO的a-CAN分析可知，此時體系主要以fcc晶型填充．根據圖6給出的在相同密度條件下Q4和Q6以及聚集體形態隨溫度降低的變化過程可知，fcc2晶體的形成過程為先發生從單粒子到二聚體的轉變，進而二聚體逐漸聚集形成有序的fcc2晶體．當ρ=1.5時，對體系進行緩慢降溫，溫度從1.0降低到0.01，觀察體系的變化，結果如圖7所示。紅色粒子是被標記的粒子，利用示蹤粒子的方法，粒子束在不同溫度下的聚集狀態以及體系的結構變化過程如圖8所示．體系經歷了無序的聚集體到有序的晶體的變化過程．由聚集體形態以及Q4，Q6隨著溫度降低的變化過程發現，在溫度為0.8時體系中的二聚體和三聚體類型的數量發生轉變，且三聚體的數目迅速超過了二聚體和單粒子的數目．在相同溫度下，Q4和Q6也迅速增大，體系的聚集態結構為fcc晶型．結合OIVTO的a-CAN分析，可以判斷出體系中的三聚體所屬類型為fcc晶體以及少部分其它晶型的混合物．溫度進一步降低至0.65時，體系基本由fcc3晶體構成．比較等密度降溫和等溫壓縮得到的結果發現，在相同溫度和壓力下通過不同的熱力學路徑得到的相結構不完全相同，這可能是由于不同有序結構之間的自由能壘較高，采用動力學方法難以逾越較高的自由能壘，因而無法達到最終的熱力學平衡態結構．因此，需要采用增強抽樣等方法才能獲得熱力學平衡態的結構。綜上所述，采用分子動力學模擬方法研究了軟膠體粒子體系形成束晶的動力學過程，探討了等溫壓縮和等密度降溫2個不同的熱力學路徑下，膠體粒子形成有序束晶的動力學差異．隨著密度的改變，體系結構經歷“無序→fcc1→fcc2和bcc共混物→fcc2→bcc2”的變化過程，粒子形態從單粒子向二聚體轉變．與蒙特卡洛模擬結果相比，分子動力學模擬得到的結構隨著密度的變化有明顯的遲滯現象．模擬等密度降溫過程發現，該體系的結構變化為“無序的聚集體→fcc2→fcc3”．從等溫壓縮和等密度降溫2個不同的熱力學路徑得到的結果不完全相同．在相同密度和溫度下，通過等溫壓縮得到的結構是bcc2，而通過等密度降溫或者蒙特卡羅模擬得到的結構是fcc3．因此，動力學過程會對模擬結果產生很大的影響，甚至得到不同的結構．這是由于形成fcc3所需的自由能壘較高，而分子動力學模擬使得體系無法克服較高的自由能壘，因而無法獲得熱力學更加穩定的fcc3結構．

作者：馬蘭 容婧婧 朱有亮 黃以能 孫昭艷 單位：伊犁師范學院物理科學與技術學院