理論背景 - 晶體成長與形貌調控之粗粒化蒙地卡羅模擬

在本章節中將會一一介紹此模擬系統中所使用的基本背景知識。並以此背景知識作為建構模擬晶體成長之基礎。

2-1 單位晶格

HAp 的分子式為 Ca10(PO4)6OH2，晶體結構為六角柱狀，圖 2-1 為 HAp 的單位晶格¹。將 HAp 中的氫氧基以氟取代可以得到 FHAp (Ca10(PO4)6F2)。FHAp 單位晶格之 a, b 軸長為 9.3684 Å 、c 軸長為 6.8841 Å 。

本研究根據 FHAp 單位晶格設計模擬之晶體單元大小，晶體單元長寬設置為 FHAp 延 c 軸堆疊 2 個單位晶格之長寬比 1.47 而定。將晶體單元長寬分別設置為 50 與 34 單位，每單位大小為 10^10 m。

圖 2-1：HAp 單位晶格示意圖¹。

2-2 Helmholtz 自由能

Helmholtz 自由能(A)在熱力學中可用來表示固定溫度下密閉系統受到有效功的熱力學位能。其基本定義如式 2-1。在統計熱力學中：固定粒子數(N)、體積(V)、溫度(T)情況下，可將 Helmholtz 自由能表示成式 2-2。

Q : canonical partition function k : Boltzmann factor

excess N kT

NVT N

ln _NVT^ideal _NVT^excess

ideal excess

模擬中，藉由計算系統之 Helmholtz 自由能得知晶體單元間的作用情形。模擬時，每次在系統中加入一顆晶體單元，晶體單元在系統中的位置以晶體單元中心 x、y 座標以及與垂直方向的夾角( )表示。新晶體單元在模擬系統中的 Helmholtz 自由能可設為以下形式：

 

能分別以 Debye length 與排斥體積觀念描述系統內能(E_area)與熵的相關項(E_asa) 表示。

2-2-1 Debye length

系統內能的設定中，先以 Debye length 來描述晶體間作用距離。Debye length 是指在稀釋溶液中離子運動並非理想情形，溶液中離子彼此會有作用³。溶液中帶電離子的運動可用帕松方程式(Poisson equation)描述，如式 2-8。當電解質溶液處於平衡時，帶電離子分布會符合波茲曼分布(Boltzmann

distribution)。溶液中正負離子濃度可由式 2-9、式 2-10 表示。

 

  

2 : Laplacian operator

 : electrostatic potential

 : charge density

式 2-即為 linearized Poisson–Boltzmann 方程又稱為 Debye–Hückel 方程，

為 Debye–Hückel 參數(式 2-15)，^¹即為 Debye length (_D)。

Debye length 可以表示帶電粒子的遮蔽距離，亦即帶電粒子在溶液中能感受到其他粒子並與之作用的距離。當特定兩粒子間的距離小於 Debye length 時，兩粒子彼此會有作用力；若特定兩粒子間的距離大於 Debye length 時，兩

向(orientation-dependent)差異所引起之熵的變化⁵。Onsager 理論可簡單的用排斥體積來描述。排斥體積表示：每個分子所佔據的空間不可讓其他分子進入。

此概念被廣泛利用在描述液晶與聚合物的理論上⁶。在液晶的向列型狀態中，當

液晶整體的排斥體積變小時，會使系統中平移熵(translational entropy)變大⁷。當溶液中的大分子排列良好會使其排斥體積會變小，讓溶液中其他小分子、粒子可移動的區域變大，增加系統熵、降低系統的 Helmholtz 自由能。

系統熵在平衡系統中一般可以寫成式 2-17，W 為粒子在系統中所有組態，

表示系統熵與系統中的組態數目有關。

ln S k 

欲計算粒子在系統中的所有組態，對粒子的描述會使用統計力學中的硬球體假設。硬球體假設認為：原子、分子存在的區域不可與其他原子或分子重疊，其作用力與距離關係可用圖 2-2 表示。當兩原子或分子靠得太近時，便會產生一強大斥力使得兩原子或分子分離；當兩原子或分子距離太遠則無作用力。當兩粒子在一維空間中，其 Helmholtz 自由能計算可表示成圖 2-3a。兩粒子在一維空間中可移動的區域如圖 2-3b、配分函數可以寫成式 2-18。圖 2-3b 之紅色部分為兩粒子間的排斥體積，會影響配分函數的計算，進而影響系統中的 Helmholtz 自由能。

圖 2-2：兩相同粒子間作用力與距離關係簡化示意圖。

式 2-17

圖 2-3：兩粒子在一維空間中(a)位置示意圖，(b)位置與分佈關係圖。紅色表示 粒子之排斥體積、藍色表示兩粒子間彼此有作用力、灰色表示粒子間可移動範圍。

將此一維空間計算推廣至二維空間中，用以表示模擬系統中晶體成長的 Helmholtz 自由能計算。在本研究中，模擬的晶體單元為長方形，晶體單元間彼此排列關係會影響排斥體積大小。將晶體單元在系統中的所有組態寫成式 2-19，

Vi 表示晶體單元在系統中的所有區域、Vex表示晶體單元在系統中的排斥體積、

V0 為此模擬系統中的單位體積。晶體單元在系統中的組態數目以晶體單元在系

統中存在的全部區域扣除排斥體積估算。當晶體單元彼此排列較鬆散時，由於晶體單元間可作用的距離不變，但排斥體積變大會使得配分函數變小，系統中的 Helmholtz 自由能下降。反之，晶體單元彼此排列整齊時，系統中的 Helmholtz 自由能較高。因此，模擬時以計算晶體單元間的排斥體積表示系統中熵的相關項。

i ex

V V

   (b) (a)

式 2-18

式 2-19

2-3 晶格排列

單位晶格有序且重複地排列會形成晶體結構⁸，晶體結構與晶體形貌相關。

單晶(single crystal)為單位晶格有週期且連續排列形成；單晶中所有單位晶格以相同方式互相鍵結且具有同樣之方向。故在單位晶格排列過程中，彼此為對齊排列不易產生錯位情形。

在本研究中將此概念運用至模擬系統：當模擬系統中沒有添加物存在時，

晶體單元會因晶格關係彼此排列整齊，不易形成錯位。然而當模擬系統中有添加物存在時，我們認為添加物會影響晶體排列，因此不考慮晶體單元間是否有錯位。

2-4 Terrace ledge kink 模型

Kossel 與 Stranski 認為：新晶體在表面出現後，受到熱力學因素會使新晶

體傾向移動至與最多相鄰原、分子結合之位置⁹。以二維空間中之正方形為例

10(圖 2-4)，鍵結數(N )與原子數(n)關係如

表 2-1。式 2-20 表示鍵結數與原子數關係之通式。系統能量受鍵結數目不同而改變，系統能量變化情形即為鍵能( )與鍵結數相乘 (式 2-21)。

圖 2-4：二維空間中之正方形晶體與其相鄰晶體示意圖。

表 2-1：相鄰原子數與鍵結數關係

Number of atoms, n Number of bonds, N

1 0

4 4

9 12

16 24

2 2 N  n n

2 2 N n n

此現象在實驗的三維空間中稱為 terrace ledge kink 模型(TLK 模型)，目前晶體成長的統計力學模型常以此為基礎。假設晶體為正立方體，晶體表面各位置如圖 2-5 所示，各位置之相鄰晶體數如表 2-2。TLK 模型描述：新加入的粒子在晶體表面生成後，會迅速移動到階梯或扭折位以增加相鄰晶體數，進而使系統能量更穩定。在本研究中之模擬參考此模型作設定：模擬系統中出現扭折位時，晶體單元會優先在此位置附近生成，考慮完扭折位附近所有可能位置後才會計算其他位置之系統能量，以得到較穩定之系統能量、增加模擬效率。

圖 2-5：正立方晶體表面各位置示意圖¹¹。

式 2-20

式 2-21

表 2-2：正立方晶體不同位置相鄰晶體數¹⁰

Atom Number of nearest neighbors

adatom 1

step adatom 2

kink atom 3

step atom 4

surface atom 5

2-5 蒙地卡羅方法

蒙地卡羅方法(Monte Carlo method)以隨機亂數取樣，使計算所得之結果符合統計熱力學¹²。蒙地卡羅方法使用亂數產生器(random number generator)產生隨機且平均的亂數；以此得到一個均勻分布的亂數數列。當產生的亂數足夠多時，能使平均值趨近於理論值；亦即需要大量重複抽樣才能得到結果。

蒙地卡羅方法在取樣中已考慮能量的波茲曼分布情形；亦即在取樣過程中已考慮各狀態存在機率。使得取樣的有效性提升，也能減少計算資源的浪費。

在取樣時，先計算系統初始能量(E )並將之隨機移動得到新系統能量(₁ E )。若₂ 新系統能量較低時接受新狀態，否則會計算新系統能量之波茲曼權重( f )，並 和一隨機數R比較。當 f 值較大，接受新狀態；否則保留原狀態。

本研究中以蒙地卡羅法搜尋晶體構型最佳化位置，在模擬系統中的計算流程如下：

1. 計算系統中新加入的晶體單元之系統能量(E )。 ₁

2. 將新晶體單元在系統中隨機移動、轉動，計算系統能量(E )。 ₂ 3. 若E₂ E₁，接受新狀態。回到步驟 1. 。

4. 若E₂ E₁，計算

E2 E1

f e kT

 

 。

5. 產生一個隨機數R介在。

6. 若 f R，接受新狀態。回到步驟 1.。

若 f R，拒絕新狀態，保留已存在的狀態。回到步驟 1. 。

2-6 參考資料

(1) Menéndez-Proupin, E.; Cervantes-Rodríguez, S.; Osorio-Pulgar, R.; Franco-Cisterna, M.; Camacho-Montes, H.; Fuentes, M. E. J. Mech. Behav. Biomed.

Mater. 2011, 4 (7), 1011.

(2) Chandler, D. Introduction to Modern Statistical Mechanics, 1st Edition edition.;

OUP USA: New York, 1987.

(3) Dill, K. A.; Bromberg, S. Molecular Driving Forces: Statistical

Thermodynamics in Biology, Chemistry, Physics, and Nanoscience; Garland Science, 2011.

(4) Muthukumar, M. J. Chem. Phys. 2009, 130, 161101.

(5) Onsager, L. Ann. N. Y. Acad. Sci. 1949, 51 (4), 627.

(6) Hill, T. L. An Introduction to Statistical Thermodynamics; Courier Corporation, 1960.

(7) Billeter, J. L. PhD Thesis 1999, 226.

(8) Shmueli, U. Theories and Techniques of Crystal Structure Determination;

Oxford University Press, UK: Oxford, GBR, 2007.

(9) Kossel, W. Nachrichten Von Ges. Wiss. Zu Gött. Math.-Phys. Kl. 1927, 1927, 135.

(10) Oura, K.; Katayama, M.; Zotov, A. V.; Lifshits, V. G.; Saranin, A. A. Surface Science; Advanced Texts in Physics; Springer Berlin Heidelberg: Berlin, Heidelberg, 2003.

(11) Vekilov, P. G. Cryst. Growth Des. 2007, 7 (12), 2796.

(12) Frenkel, D.; Smit, B. Understanding Molecular Simulation, Second Edition:

From Algorithms to Applications, 2 edition.; Academic Press: San Diego, 2001.

在文檔中晶體成長與形貌調控之粗粒化蒙地卡羅模擬 (頁 37-50)