第35卷第1期 人 工 晶 体 学 报 Vol . 35 No . 1 2006年2月 JOURNAL OF SY NTHETI C CRYST ALS February, 2006
晶体生长机制和生长动力学的蒙特卡罗模拟研究朱 阁, 卢贵武, 李英峰, 蓝建慧, 张 军, 郑庆彬, 黄乔松, 孙 洵, 夏海瑞
(1. 中国石油大学物理科学与技术学院, 东营257061; 2. 中国石油大学数理系, 北京102249;
3. 山东大学晶体材料国家重点实验室, 济南250100; 4. 山东大学物理与微电子学院, 济南250100) 121111134
摘要:采用动力学蒙特卡罗方法, 对在完整光滑界面上低过饱和度溶液中的晶体生长机制和动态过程进行计算机模拟, 得到了晶体生长速率与溶液过饱和度之间的关系以及晶体生长的表面形态的动力学生长规律进行分析, 值, 讨论了热粗糙度、表面扩散、关键词:蒙特卡罗模拟; 生长机制; 晶体表面形态; 中图分类号:O361 文献标识码:A :2985X (2006) 0120024208
M on te S on Crysta l
is m and K i n eti cs
LU Gui 2w u , L I Ying 2feng , LAN J ian 2hui , ZHAN G Jun ,
ZHEN G Q ing 2bin , HUAN G Q iao 2song , SUN X un , X I A Hai 2rui 11342111
(1. College of Physics Science and Technol ogy, China University of Petr oleum, Dongying 257061, China;
2. Depart m ent of Math and Physics, China University of Petr oleum, Beijing 102249, China;
3. State Key Laborat ory of CrystalMaterials, Shandong University, J inan 250100, China;
4. College of Physics and M icr o 2electr onics, Shandong University, Jinan 250100, China )
(Received 1A ugust 2005, accepted 18O ctober 2005)
Abstract:An i m p r oved kinetic Monte Carl o method was adop ted t o si m ulate the dyna m ic p r ocess of crystal gr owth on the perfect and s mooth surface in order t o penetrate the m icr ocos m ic gr owth p r ocess of crystal, es pecially the crystal gr owth mechanis m in the l ow supersaturated s oluti on . A s a result, the relati onshi p bet w een the crystal gr owth rate and the supersaturati on of the s oluti on was given, as well as the surface mor phol ogy of crystal gr owth . The kinetic la w of crystal gr owth was analyzed, which indicated the t w o 2di m ensi onal nucleati on being the main gr owth mechanis m. The gr owth dead z one of t w o 2di m ensi onal nucleati on gr owth and a critical point transfor m ing the mononuclear gr owth int o a multi p le nucleati on gr o wth were discovered . Ther mal r oughness, surface diffusi on, average height of the step s and the i m pact of the surface size on the average gr owth rate of the crystal are als o discussed in this paper .
Key words:Monte Carl o si m ulati on; gr o wth mechanis m; crystal surface mor phol ogy; crystal gr owth kinetics 1 引 言
目前, 对于晶体生长机制和生长动力学, 国内外已开展了大量的实验研究。例如A. McPhers on 等收稿日期:2005208201; 修订日期:2005210218
基金项目:国家自然科学基金(No . 10274043) ; 山东省自然科学基金(No
. Y2003A01)
作者简介:朱阁(19782) , 男, 山东省人, 硕士研究生。
通讯作者:卢贵武, E 2mail:lugui w u@s ohu . com [1]采
第1期 朱 阁等:
晶体生长机制和生长动力学的蒙特卡罗模拟研究25用原子力显微镜(AF M ) 观察了法向生长、位错生长、二维成核生长以及三维成核生长机制下的晶体生长形
[2]态图, 并在分子水平上进一步研究了缺陷对生长机制的影响; 于锡玲等采用相衬2微分干涉显微术连接全
[3]息图像分析系统, 实时研究了有机金属配合物晶体的截面动力学过程; 卢贵武等采用迈克尔逊干涉实验
装置实时观测K DP 晶体的生长动力学规律。这些实验结果一方面给我们提供了更直接和直观的晶体生长信息, 另一方面也验证和丰富了晶体生长理论(如三维成核生长机制的发现) 。计算机晶体生长是建立在晶体生长理论和实验基础之上, 可以比较真实地反映晶体生长微观过程, 深入研究晶体生长机理和进一步指导
[4]晶体生长实验的一项重要研究方法。自上世纪70年代计算机模拟方法被应用到晶体生长的研究以来, 计
[527]算机模拟晶体生长已经从最初的无机晶体发展到有机晶体和蛋白质等生物大分子晶体, 计算机晶体生长
[8210]的尺度也开始从介观尺度发展到原子尺度。传统的蒙特卡罗(MC ) 模拟方法在当今计算机运算性能显
著提高的情况下其优越性也得到了越来越好的发挥。
相比国外计算机晶体生长的研究, 国内这方面的文献还很鲜见。本文将结合晶体生长动力学理论采用MC 方法对低过饱和度溶液下的晶体生长过程进行计算机模拟, 研究晶体生长表面形态和生长机制, 并着重讨论晶体生长条件的变化对晶体生长速率的影响。
2 蒙特卡罗模拟
模拟在简立方(001) 晶面进行, , Kossel 模型, 生长基元设为正立方体, 且只考虑第一近邻间相互作用。, 排除基元空位和附着悬挂。2晶相之间和晶相2液相之间的相互作用, 模拟中采用周期性边界条件。
晶相2液:生长基元被吸附到晶体表面, 生长基元脱附到溶液, 。晶体表面每个座点配位数的不同, 晶体和溶液的热力学性质的差异() 以及界面层内生长基元数量(与溶液过饱和度成比例) 均影响三种基本事件发生的概率。基元吸附频率k i 、脱附频率k i 和扩散频率k ij 具体可用公式表述为
k i
k i +-+-[8, 11]:(1) (2) =f t ・exp [-α(2-i ) /4+β/2]=f t ・exp [α(2-i ) /4-β/2]
2k ij =f t ・(X S /a ) ・exp [α(2+j -i ) /4-β/2](3)
其中i (=0, 1, 2, 3, 4) 是晶体表面生长基元的最近邻
配位数, 在S OS 假定下不考虑基元的垂直键合。扩散
频率k ij 中的i 代表该基元在扩散前的最近邻配位数, j
代表扩散到另一座点时的配位数。图1给出了晶体
生长过程中三种基本事件的示意图。从(a ) →(b )
和从(c ) →(b ) 均为脱附过程, (b ) →(c ) 及(b ) →
(a ) 为简单吸附过程, (a ) 和(c ) 之间为一可逆的扩
散过程。事件发生的概率要依据座点周围配位数来
+计算。例如在计算(b ) →(a ) 的吸附频率k i 时, 吸
附点处在有两个最近邻生长基元的扭折位置, 因此配
位数i 取值为2。而在计算事件(c ) →(a ) 的扩散频
率k ij 时, 配位数需要分别计算它的脱附配位数和吸附
配位数, i 和j 的取值应分别为4和2。图中可以看到
扩散事件可以被看作是脱附和吸附两个事件的组合。
依据表面自由能最小原理, 图1中(b ) →(c ) 的吸附
概率要比从(b ) →(a ) 的吸附概率大, 也就是说(b )
→(c ) 的吸附过程更易发生。同样从(a ) →(c ) 的
扩散过程比其逆过程更易发生。图1 平衡体系下晶体表面存在的三种可能构型以及三种基本事件示意图Fig . 1 Possible surface configurati ons in the equilibriu m ense mble of syste m s . A si m p le creati on converts (b ) int o (a ) and an annihilati on converts (c ) int o (b ) . An ele mentary p r ocess bet w een (a ) and (c ) is a surface diffusi on ju mp. Those p r obabilities depend on the nu mbers of their lateral s olid neighbors
26人工晶体学报 第35卷
式(1) 2(3) 中f t 代表单位时间出现的尝试数, 它由粘性流体的活化自由能ΔG 决定, 即:
f t =exp (-ΔG /kT ) h [8](4) 其中h 是普朗克常数, k 是玻耳兹曼常数, T 是温度。参数α用于描述晶体表面的热粗糙度, 其定义可表述为:
(5) α=[2(
体基元对的键能。其中
参数β用于描述结晶成核的驱动力(溶液过饱和度) , 有时称为动力学粗糙系数, 它与流体和固体间化学势μ的差值成正比, 即:
β=Δμ/kT
取值范围0～2。
在Gil m er 和Benne ma 的模型中, 似计算。这里我们采用所有座点事件可能发生频率总和K , 并假设事件的发生是顺序的且中间没有时间间隔, 即[11](6) (3) 式中的X S 为基元的平均扩散距离, a 为正立方体生长基元的边长X S /a的:
/K
5
-
i (x, y ) K =x y k x y +k +j (x, y ) +
q =1∑k ) i (x, y ) , j (x ’, y ’(7)
其中(x, y ) 。表面扩散就是基元方块在某一表面座点(有i 个近邻方块) 离开(有j 个近邻方块) 进入晶体表面两种事件的组合, 所以可能的扩散方向不仅包括最近邻的四个方向, 而且还应包括本身位置不变的情况。扩散方向的选取根据该座点周围环境的不同, 即扩散至不同方向去的扩散概率的不同而随机选取。这一算法的改进与文献[11]中提出的扩散模型(即只考虑四
) 表示扩散的终止位置。个方向, 且四个扩散方向的几率相等) 相比显然更具合理性。(x ’, y ’
因此晶面的法向生长速率可表示为:
r =t N (8)
其中ΔN 是在t 时间内晶体表面固体方块数目的增量, d 是一个生长层的厚度(即生长单元的高度) , N 是在一个生长层中所有生长基元立方块的总数量。
本模拟我们着重研究参数α, β, X S 对生长速率的影响, 因此引入平均生长速率:
r ’=df t (9)
该物理量为一无量纲参数。
模拟中得到的生长速率和下文中没有特别说明的生长速率均是无量纲的法向平均生长速率。
本文的计算机蒙特卡罗模拟程序运用Fortran 语言自行编写, 程序的具体算法和流程参见文后附录A 。3 结果和讨论
本文对不同条件下的晶体生长进行了大量的模拟计算, 研究了生长速率与过饱和度、热粗糙度、扩散距离、表面尺寸以及平均生长高度等参量之间的关系。
3. 1 晶体生长机制
3. 1. 1 法向生长机制
图2给出了在热粗糙度α较小, 晶体生长到平均台阶高度h =3. 3时, 过饱和度β以及平均扩散距离X s
对平均生长速率r ’的影响。从图中可以看出, 生长速率r ’随过饱和度β的增加其曲线基本上均呈线性上升趋势。在α=3时, 曲线的斜率较小, 生长速率的增长较平缓。图3给出了在不同热粗糙度下模拟得到的晶体生长表面形态, 从图3(a ) 和3(b ) 可以看出, 在热粗糙度比较小时(α
。
3. 1. 2 二维成核生长机制
当α>4时晶体的生长机制由法向生长转变为二维成核生长, 表面生长形态如图3(c ) 和3(d ) 所示。两种生长机制造成的晶体表面形态上的差异在图3中可以清晰地看到, 相比较二维成核生长机制下的晶体表面更光滑和有规律。Vol m er 认为:结晶质点从溶液到达晶体表面上作横向移动, 并形成一个吸附层。由于热力学涨落, 层中质点做频繁的非弹性碰撞, 可形成二维晶体小岛即二维晶核。这种二维核一旦形成, 新质点结合到小岛的台阶上, 就不会遇到任何能量上的困难, 于是新层开始生长。20世纪50年代初, Burt on, Cabrera 和Frank (BCF ) 对二维成核机制作了普遍化的处理[14], 认为形成一个超临界半径的晶核需要大的活化能。所以在极低过饱和度时, 通过二维成核生长不大可能[15]。在我们的模拟当中也发现了文献[16]中所指出的这种所谓" 生长死区" 。当矩阵大小为20×20, α=4、5、6时, 分别存在β=0. 05、0. 15、0. 35的晶体生长临界值(见图4) , 在小于此临界值的过饱和度区域内无论是否考虑扩散, 晶体都无法成核生长。但是如果事先在光滑表面上设置一定大小的晶核, 生长基元将围绕该晶核铺满全层, 然后生长停止。如果表面存在固定的和永不消失的位错露头, 即使在临界值以下晶体也会按螺旋位错的生长机制持续生长。
二维成核生长机制根据成核率和台阶扩展速度的关系还可以细分为单核生长和多核生长[12]。单核生长是指成核速率较小, 在相当长时间内不可能形成二维临界核, 一旦出现一个二维临界核, 就很快地形成一新的结晶层, 相对来说台阶扩展速度要远远大于成核率。相反, 如果成核率比较大, 而台阶横向扩展的速度
却相对比较慢时, 生长界面会同时存在许多稳定的二维晶核
生长, 而且若相邻台阶相遇, 它们可以不留痕迹地合并成新
的结晶层。另外还可能在比较大的晶核之上又有新的晶核
形成, 这种二维生长机制被称为多核生长。单核生长机制在
实际晶体生长中很少能观测到, 原因是实际晶体表面相对我
们模拟的60×60的晶体表面要大的多。但在过饱和度较低
和晶体表面足够小时, 这种单核生长机制是完全有可能的,
在我们的模拟中也的确观察到了这种生长机制。图5为模
拟结果中截取的两张分别代表单核和多核生长机制的生长
表面形态图。
整体而言, 二维成核生长的动力学规律即生长速率与生
长驱动力(溶液过饱和度) 之间的函数关系为指数关系。其中图4 热粗糙度较大时生长速率与过饱和度之间的关系Fig . 4 The relati onshi p bet w een the gr owth rate and upersaturati on with the larger ther mal r oughening coefficient 单核生长速率随过饱和度的增大基本上成线性增长, , 即呈指数
4增长。因此在两种生长机制之间有一个临界值, 通过图456) 时, 二维成核
生长机制从单核生长过渡到多核生长时的临界值大致分别在0. 3、0. [11]中得到的结果基本一致。我们注意到α=4, 生长速率比较快。α=6时的生长曲线在β=0. 35～0. 75, 3(d ) 和5(a ) 所示, 晶体生长质量较高,
图5 以二维成核机制在完整光滑平面上的晶体生长表面形态图
Fig . 5 The surface configurati on of crystal gr o wth on the perfect and s mooth p lane under the t w o 2di m ensi onal nucleati on
mechanis m (surface matrix is 60×60) :(a ) T wo 2di m ensi onal single nucleati on gr o wth mechanis m ;
(b ) T wo 2di m ensi onal multi p le nucleati on gr o wth mechanis m
3. 2 影响晶体生长速率的其它因素
3. 2. 1 表面扩散对晶体生长速率的影响
扩散是晶体表面存在的三种基本事件之一, 也是本文模拟中重点考查的对象。图2给出了平均扩散距离X S 对生长速率的影响。结果发现考虑扩散时的生长速率比不考虑扩散情况时要大, 扩散距离大小与生长速率成正比关系, 而且随着溶液过饱和度的增大扩散的影响更为显著。但整体而言, 表面扩散的影响并不引起生长机制的变化, 从图2中的生长曲线来看依然是线性关系。图6给出了在不同扩散距离下生长速率与热粗糙度之间的关系。可以清楚的看到, 在相同的溶液过饱和度下随着热粗糙度α的不断增大, 扩散对生长速率的影响逐渐减小。也就是说表面扩散对法向生长机制的影响更大, 而对二维成核机制尤其是单核生长机制的影响比较小。从图中还可以看到当α>4时生长速率趋缓, 可说明生长机制由法向生长过渡到了二维成核生长机制。
3. 2. 2 表面尺寸对晶体生长速率的影响
, 图7给整体来说在小的晶体表面上模拟得到尤其是在以单核机制生长时这种差别是非常明显的, 随着过饱和度的增大, 尺寸较小的晶体表面生长速率起伏开始增大。晶3. 2. 3 平均生长高度对晶体生长速率的影响
本文计算的晶体生长速率为平均生长速率, 我们发现相
同条件下晶体生长到不同台阶高度所对应的平均生长速率是
有差别的, 如图8所示, 生长机制以单核生长为主时, 这种因
生长到不同台阶高度而产生的平均生长速率的差别是非常小
的; 生长机制过渡到多核生长后, 随着过饱和度的增大, 平均
生长速率的差别开始越来越大。例如h =1. 1与h =2. 1时的
差别明显比h =2. 1与h =3. 1的差别要大。我们发现, 最初
在台阶高度为1的光滑平面上成核生长相对所花时间较长,
而在多核生长机制下台阶的扩展和新晶核的形成是同时的,
生长点急剧增多, 因此晶体生长速率会逐渐增大, 其平均生长
速率渐渐趋于一个稳定值。所以在模拟时我们尽量令晶体生图8 生长台阶的平均高度对生长速率的影响Fig . 8 The relati onshi p bet w een the gr o wth rate and mean gr o wth height
长到较高的同一台阶高度(例如h =3. 3) , 以便得到比较稳定的和具有可比性的平均生长速率。
图8中的生长曲线在β≥0. 5后以锯齿形上升, 原因就在于表面矩阵没有达到足够大而引起的生长速率的不稳定。这种现象也说明了计算机晶体生长本身具有的实验属性。因此要在考虑计算机运算速度和实际模拟所用时间的基础上尽量选取较大尺寸的晶体表面进行模拟。
4 结 论
(1) 作为在完整光滑晶面上进行的晶体生长模拟, 主要的晶体生长机制是二维成核生长。当然在热粗糙度比较低时也可以观察到法向生长机制, 而在过饱和度非常低时也发现了二维成核生长存在的生长死区;
(2) 通过模拟出的生长曲线结合晶体生长的表面形态可以进一步把二维成核生长机制细分为单核生长和多核生长。证实在热粗糙度比较大、过饱和度又比较低时, 线性的单核生长是主要生长机制;
(3) 表面扩散对晶体生长速率有一定影响。扩散优化了生长基元在晶体表面位置的排放, 使得晶体生长的表面能更小, 表面形态更趋稳定, 提高了晶体的生长速率。需要说明的是这种影响并不引起生长机制的变化, 而且随着热粗糙度的增大这种影响也越来越小;
(4) 计算机模拟的晶体表面尺寸对生长速率存在一定影响, 尤其是在单核生长时, 晶体生长速率随表面尺寸的增大明显增大。在多核生长机制下表面尺寸大时晶体生长速率趋缓, 相反在晶体表面尺寸小时生长曲线局部起伏较大。但表面大小对晶体生长曲线的整体走势没有影响;
(5) 晶体生长到不同平均高度时得到的平均生长速率是不同的, 在计算平均生长速率时, 需选定一个可比较的标准, 例如一定的生长时间或者本文采取的同一台阶平均高度, 而不是通常情况下简单的蒙特卡罗模拟步数的递加。
计算机晶体生长实质上是对晶体复杂生长过程的理解和反映。由于实际上不可能对微观的生长基元进行跟踪和定位, 因此注重模拟结果的真实具有重要意义。主要生长机制的晶体生长过程, 的应用前景, 可以进一步结合对电解质溶液的认识, , 采用E wald 求和等方法考虑长在第2, 建立了晶体生长的理论模型, 在此基础上我们给出:
(1) 初始化参量
h (x , y ) 表示晶体表面不同位置的高度, 初始值设为1, {x , y }表示晶体表面矩阵, x ∈(1, x max ) , y ∈(1, y max ) 。依据各点高度分别计算表面各位置的吸附配位数和脱附配位数, 以便得到k i , k i , k ij , K xy , K x , K 值。并预先设置基元生长个数h N, 基元个数n =0, 其中K x =+-∑K
y xy ,N 为每一生长层内的基元个数。
(2) 产生一个在(0, K ) 区间均匀分布的随机数ξ
(3) 随机选取某一点(x , y ) :若x
若y
(4) 选择事件:若ξ
否则若ξ
其它情形, 基元扩散, 取q =ξ-k i (x, y ) -k j (x, y )
(5) 选择扩散方向
若q
否则若q
令q =q -k i (x, y ) , j (x +1, y ) ;
否则若q
令q =q -k i (x, y ) , j (x, y -1) ;
否则若q
其它则h (x , y ) =h (x , y ) 。
以上寻找坐标点(x, y ) 以及选择事件的过程可以参见图9, 随机数按照标注1, 2, 3的顺序依次进行选择。+-+-
图9 一个蒙特卡罗循环步的流程图
Fig . 9 Fl ow chart of a Monte Carl o circulati on:(a ) fix a site (x, y ) ; (b ) choose a certain event
(6) 增加时间记数:t =t +1/K
(7) 刷新(x , y ) 及周围12个点的配位数, 12个点的配位数, 以便得到新的K xy , K x , K
(8) 若n
(9) 输出h (x , y ) , 以及生长速率等结果。
参[1]
[2]McPhers on A lexander, G, omolecular Crystal Gr owth as Revealed byA t om ic ForceM icr oscopy [J ].Journal of S tructural B iology , :32. Yu X L, W et al . I nterface Kinetics and Crystal Gr owth Mechanis m of a Ne w O rganometallic Coordinati on Compound:
Triallylthi ourea B r om ide [J ].Journal of C rystal Gro w th , 1997, 182:4282433.
[3]卢贵武, 李春喜, 等. K DP 晶体台阶生长动力学的激光干涉实验研究[J ].化学物理学报, 2003, 4:2892292.
[4]
[5]
[6]
[7]
[8]
[9]
[10]
[11]
[12]Andrea C L, M ir oslav Kotrla . Theory and Si m ulati on of Crystal Gr owth [J ].J. Phys . :Condens M atte . , 1997, 9:2992344. Cheng V K W et al . The Kinetic A sy mmetry bet w een Nucleati on Gr owth and D iss oluti on:a Monte Carl o Study [J ].Journal of C rystal Gro w th , 1989, 99(1) :2932299. Van EnckevortW J P, et al . I m purity B l ocking of Crystal Gr owth:a Monte Carl o Study [J ].Journal of C rystal Gro w th , 1998, 183(2) :4412448. Cheng V K W. A Monte Carl o Study ofMoving Step s during Crystal Gr owth and D iss oluti on [J ].Journal of C rystal Gro w th , 1993, 134(2) :3562365. Gil m er G H, Bennema P . Si m ulati on of Crystal Gr owth with Surface D iffusi on [J ].J. Appl . Phys . , 1972, 43:134721360. Lee Jae 2Wook, Hwang Nong M , Ki m Doh 2Yeon . Gr owth Mor phol ogy of Perfect and T winned Face 2centered 2cubic Crystal by Monte Carl o Si m ulati on [J ].Journal of C rystal Gro w th , 2003, 250:5382545. W ats on Grae me W , O liver PeterM , Parker Stephen C . A t om istic Si m ulati on of Crystal Gr owth at the a 〈100〉Scre w D isl ocati on Ter m inating at the {100}Surface of Mg O [J ].Surface Science , 2001, 474:1852190. Rak M ir osla wa, I zdebskiMarek, B r ozi Andrzej . Kinetic Monte Carl o Study of Crystal Gr owth fr om Soluti on [J ].Co m p . Phys . Co mm . , 2001, 138:2502263. Malkin A J, Kuznets ov Yu G, Land T A, et al . Mechanis m s of Gr owth f or Pr otein and V irus Crystals [J ].N ature S tructural biology , 1995, 2
(11) :9562959.
[13]张克从, 张乐惠, 主编. 晶体生长科学与技术(第二版) 上册[M].北京:科学出版社, 1997.
[14]
[15]
[16]Burt on W K, Cabrera N, Frank F C . The Crystal of Crystals [J ].Phil. Trans . , 1951, A243:2992358. De Yoreo, Land T A, Dair B. Gr owth Mor phol ogy of V icinal H ill ocks on the {101}Face of KH 2P O 4:fr om Step 2fl ow t o Layer 2by 2layer Gr owth [J ].Physical Revie w L etters , 1994, 73:8382841. Malkin A J, Kuznets ov Y G et al . I n Situ A t om ic Force M icr oscopy Studies of Surface Mor phol ogy, Gr owth Kinetics, Defect Structure and
D iss oluti on in Macr omolecular Crystallizati on [J ].Journal of C rystal Gro w th , 1999, 196:4712488.
第35卷第1期 人 工 晶 体 学 报 Vol . 35 No . 1 2006年2月 JOURNAL OF SY NTHETI C CRYST ALS February, 2006
晶体生长机制和生长动力学的蒙特卡罗模拟研究朱 阁, 卢贵武, 李英峰, 蓝建慧, 张 军, 郑庆彬, 黄乔松, 孙 洵, 夏海瑞
(1. 中国石油大学物理科学与技术学院, 东营257061; 2. 中国石油大学数理系, 北京102249;
3. 山东大学晶体材料国家重点实验室, 济南250100; 4. 山东大学物理与微电子学院, 济南250100) 121111134
摘要:采用动力学蒙特卡罗方法, 对在完整光滑界面上低过饱和度溶液中的晶体生长机制和动态过程进行计算机模拟, 得到了晶体生长速率与溶液过饱和度之间的关系以及晶体生长的表面形态的动力学生长规律进行分析, 值, 讨论了热粗糙度、表面扩散、关键词:蒙特卡罗模拟; 生长机制; 晶体表面形态; 中图分类号:O361 文献标识码:A :2985X (2006) 0120024208
M on te S on Crysta l
is m and K i n eti cs
LU Gui 2w u , L I Ying 2feng , LAN J ian 2hui , ZHAN G Jun ,
ZHEN G Q ing 2bin , HUAN G Q iao 2song , SUN X un , X I A Hai 2rui 11342111
(1. College of Physics Science and Technol ogy, China University of Petr oleum, Dongying 257061, China;
2. Depart m ent of Math and Physics, China University of Petr oleum, Beijing 102249, China;
3. State Key Laborat ory of CrystalMaterials, Shandong University, J inan 250100, China;
4. College of Physics and M icr o 2electr onics, Shandong University, Jinan 250100, China )
(Received 1A ugust 2005, accepted 18O ctober 2005)
Abstract:An i m p r oved kinetic Monte Carl o method was adop ted t o si m ulate the dyna m ic p r ocess of crystal gr owth on the perfect and s mooth surface in order t o penetrate the m icr ocos m ic gr owth p r ocess of crystal, es pecially the crystal gr owth mechanis m in the l ow supersaturated s oluti on . A s a result, the relati onshi p bet w een the crystal gr owth rate and the supersaturati on of the s oluti on was given, as well as the surface mor phol ogy of crystal gr owth . The kinetic la w of crystal gr owth was analyzed, which indicated the t w o 2di m ensi onal nucleati on being the main gr owth mechanis m. The gr owth dead z one of t w o 2di m ensi onal nucleati on gr owth and a critical point transfor m ing the mononuclear gr owth int o a multi p le nucleati on gr o wth were discovered . Ther mal r oughness, surface diffusi on, average height of the step s and the i m pact of the surface size on the average gr owth rate of the crystal are als o discussed in this paper .
Key words:Monte Carl o si m ulati on; gr o wth mechanis m; crystal surface mor phol ogy; crystal gr owth kinetics 1 引 言
目前, 对于晶体生长机制和生长动力学, 国内外已开展了大量的实验研究。例如A. McPhers on 等收稿日期:2005208201; 修订日期:2005210218
基金项目:国家自然科学基金(No . 10274043) ; 山东省自然科学基金(No
. Y2003A01)
作者简介:朱阁(19782) , 男, 山东省人, 硕士研究生。
通讯作者:卢贵武, E 2mail:lugui w u@s ohu . com [1]采
第1期 朱 阁等:
晶体生长机制和生长动力学的蒙特卡罗模拟研究25用原子力显微镜(AF M ) 观察了法向生长、位错生长、二维成核生长以及三维成核生长机制下的晶体生长形
[2]态图, 并在分子水平上进一步研究了缺陷对生长机制的影响; 于锡玲等采用相衬2微分干涉显微术连接全
[3]息图像分析系统, 实时研究了有机金属配合物晶体的截面动力学过程; 卢贵武等采用迈克尔逊干涉实验
装置实时观测K DP 晶体的生长动力学规律。这些实验结果一方面给我们提供了更直接和直观的晶体生长信息, 另一方面也验证和丰富了晶体生长理论(如三维成核生长机制的发现) 。计算机晶体生长是建立在晶体生长理论和实验基础之上, 可以比较真实地反映晶体生长微观过程, 深入研究晶体生长机理和进一步指导
[4]晶体生长实验的一项重要研究方法。自上世纪70年代计算机模拟方法被应用到晶体生长的研究以来, 计
[527]算机模拟晶体生长已经从最初的无机晶体发展到有机晶体和蛋白质等生物大分子晶体, 计算机晶体生长
[8210]的尺度也开始从介观尺度发展到原子尺度。传统的蒙特卡罗(MC ) 模拟方法在当今计算机运算性能显
著提高的情况下其优越性也得到了越来越好的发挥。
相比国外计算机晶体生长的研究, 国内这方面的文献还很鲜见。本文将结合晶体生长动力学理论采用MC 方法对低过饱和度溶液下的晶体生长过程进行计算机模拟, 研究晶体生长表面形态和生长机制, 并着重讨论晶体生长条件的变化对晶体生长速率的影响。
2 蒙特卡罗模拟
模拟在简立方(001) 晶面进行, , Kossel 模型, 生长基元设为正立方体, 且只考虑第一近邻间相互作用。, 排除基元空位和附着悬挂。2晶相之间和晶相2液相之间的相互作用, 模拟中采用周期性边界条件。
晶相2液:生长基元被吸附到晶体表面, 生长基元脱附到溶液, 。晶体表面每个座点配位数的不同, 晶体和溶液的热力学性质的差异() 以及界面层内生长基元数量(与溶液过饱和度成比例) 均影响三种基本事件发生的概率。基元吸附频率k i 、脱附频率k i 和扩散频率k ij 具体可用公式表述为
k i
k i +-+-[8, 11]:(1) (2) =f t ・exp [-α(2-i ) /4+β/2]=f t ・exp [α(2-i ) /4-β/2]
2k ij =f t ・(X S /a ) ・exp [α(2+j -i ) /4-β/2](3)
其中i (=0, 1, 2, 3, 4) 是晶体表面生长基元的最近邻
配位数, 在S OS 假定下不考虑基元的垂直键合。扩散
频率k ij 中的i 代表该基元在扩散前的最近邻配位数, j
代表扩散到另一座点时的配位数。图1给出了晶体
生长过程中三种基本事件的示意图。从(a ) →(b )
和从(c ) →(b ) 均为脱附过程, (b ) →(c ) 及(b ) →
(a ) 为简单吸附过程, (a ) 和(c ) 之间为一可逆的扩
散过程。事件发生的概率要依据座点周围配位数来
+计算。例如在计算(b ) →(a ) 的吸附频率k i 时, 吸
附点处在有两个最近邻生长基元的扭折位置, 因此配
位数i 取值为2。而在计算事件(c ) →(a ) 的扩散频
率k ij 时, 配位数需要分别计算它的脱附配位数和吸附
配位数, i 和j 的取值应分别为4和2。图中可以看到
扩散事件可以被看作是脱附和吸附两个事件的组合。
依据表面自由能最小原理, 图1中(b ) →(c ) 的吸附
概率要比从(b ) →(a ) 的吸附概率大, 也就是说(b )
→(c ) 的吸附过程更易发生。同样从(a ) →(c ) 的
扩散过程比其逆过程更易发生。图1 平衡体系下晶体表面存在的三种可能构型以及三种基本事件示意图Fig . 1 Possible surface configurati ons in the equilibriu m ense mble of syste m s . A si m p le creati on converts (b ) int o (a ) and an annihilati on converts (c ) int o (b ) . An ele mentary p r ocess bet w een (a ) and (c ) is a surface diffusi on ju mp. Those p r obabilities depend on the nu mbers of their lateral s olid neighbors
26人工晶体学报 第35卷
式(1) 2(3) 中f t 代表单位时间出现的尝试数, 它由粘性流体的活化自由能ΔG 决定, 即:
f t =exp (-ΔG /kT ) h [8](4) 其中h 是普朗克常数, k 是玻耳兹曼常数, T 是温度。参数α用于描述晶体表面的热粗糙度, 其定义可表述为:
(5) α=[2(
体基元对的键能。其中
参数β用于描述结晶成核的驱动力(溶液过饱和度) , 有时称为动力学粗糙系数, 它与流体和固体间化学势μ的差值成正比, 即:
β=Δμ/kT
取值范围0～2。
在Gil m er 和Benne ma 的模型中, 似计算。这里我们采用所有座点事件可能发生频率总和K , 并假设事件的发生是顺序的且中间没有时间间隔, 即[11](6) (3) 式中的X S 为基元的平均扩散距离, a 为正立方体生长基元的边长X S /a的:
/K
5
-
i (x, y ) K =x y k x y +k +j (x, y ) +
q =1∑k ) i (x, y ) , j (x ’, y ’(7)
其中(x, y ) 。表面扩散就是基元方块在某一表面座点(有i 个近邻方块) 离开(有j 个近邻方块) 进入晶体表面两种事件的组合, 所以可能的扩散方向不仅包括最近邻的四个方向, 而且还应包括本身位置不变的情况。扩散方向的选取根据该座点周围环境的不同, 即扩散至不同方向去的扩散概率的不同而随机选取。这一算法的改进与文献[11]中提出的扩散模型(即只考虑四
) 表示扩散的终止位置。个方向, 且四个扩散方向的几率相等) 相比显然更具合理性。(x ’, y ’
因此晶面的法向生长速率可表示为:
r =t N (8)
其中ΔN 是在t 时间内晶体表面固体方块数目的增量, d 是一个生长层的厚度(即生长单元的高度) , N 是在一个生长层中所有生长基元立方块的总数量。
本模拟我们着重研究参数α, β, X S 对生长速率的影响, 因此引入平均生长速率:
r ’=df t (9)
该物理量为一无量纲参数。
模拟中得到的生长速率和下文中没有特别说明的生长速率均是无量纲的法向平均生长速率。
本文的计算机蒙特卡罗模拟程序运用Fortran 语言自行编写, 程序的具体算法和流程参见文后附录A 。3 结果和讨论
本文对不同条件下的晶体生长进行了大量的模拟计算, 研究了生长速率与过饱和度、热粗糙度、扩散距离、表面尺寸以及平均生长高度等参量之间的关系。
3. 1 晶体生长机制
3. 1. 1 法向生长机制
图2给出了在热粗糙度α较小, 晶体生长到平均台阶高度h =3. 3时, 过饱和度β以及平均扩散距离X s
对平均生长速率r ’的影响。从图中可以看出, 生长速率r ’随过饱和度β的增加其曲线基本上均呈线性上升趋势。在α=3时, 曲线的斜率较小, 生长速率的增长较平缓。图3给出了在不同热粗糙度下模拟得到的晶体生长表面形态, 从图3(a ) 和3(b ) 可以看出, 在热粗糙度比较小时(α
。
3. 1. 2 二维成核生长机制
当α>4时晶体的生长机制由法向生长转变为二维成核生长, 表面生长形态如图3(c ) 和3(d ) 所示。两种生长机制造成的晶体表面形态上的差异在图3中可以清晰地看到, 相比较二维成核生长机制下的晶体表面更光滑和有规律。Vol m er 认为:结晶质点从溶液到达晶体表面上作横向移动, 并形成一个吸附层。由于热力学涨落, 层中质点做频繁的非弹性碰撞, 可形成二维晶体小岛即二维晶核。这种二维核一旦形成, 新质点结合到小岛的台阶上, 就不会遇到任何能量上的困难, 于是新层开始生长。20世纪50年代初, Burt on, Cabrera 和Frank (BCF ) 对二维成核机制作了普遍化的处理[14], 认为形成一个超临界半径的晶核需要大的活化能。所以在极低过饱和度时, 通过二维成核生长不大可能[15]。在我们的模拟当中也发现了文献[16]中所指出的这种所谓" 生长死区" 。当矩阵大小为20×20, α=4、5、6时, 分别存在β=0. 05、0. 15、0. 35的晶体生长临界值(见图4) , 在小于此临界值的过饱和度区域内无论是否考虑扩散, 晶体都无法成核生长。但是如果事先在光滑表面上设置一定大小的晶核, 生长基元将围绕该晶核铺满全层, 然后生长停止。如果表面存在固定的和永不消失的位错露头, 即使在临界值以下晶体也会按螺旋位错的生长机制持续生长。
二维成核生长机制根据成核率和台阶扩展速度的关系还可以细分为单核生长和多核生长[12]。单核生长是指成核速率较小, 在相当长时间内不可能形成二维临界核, 一旦出现一个二维临界核, 就很快地形成一新的结晶层, 相对来说台阶扩展速度要远远大于成核率。相反, 如果成核率比较大, 而台阶横向扩展的速度
却相对比较慢时, 生长界面会同时存在许多稳定的二维晶核
生长, 而且若相邻台阶相遇, 它们可以不留痕迹地合并成新
的结晶层。另外还可能在比较大的晶核之上又有新的晶核
形成, 这种二维生长机制被称为多核生长。单核生长机制在
实际晶体生长中很少能观测到, 原因是实际晶体表面相对我
们模拟的60×60的晶体表面要大的多。但在过饱和度较低
和晶体表面足够小时, 这种单核生长机制是完全有可能的,
在我们的模拟中也的确观察到了这种生长机制。图5为模
拟结果中截取的两张分别代表单核和多核生长机制的生长
表面形态图。
整体而言, 二维成核生长的动力学规律即生长速率与生
长驱动力(溶液过饱和度) 之间的函数关系为指数关系。其中图4 热粗糙度较大时生长速率与过饱和度之间的关系Fig . 4 The relati onshi p bet w een the gr owth rate and upersaturati on with the larger ther mal r oughening coefficient 单核生长速率随过饱和度的增大基本上成线性增长, , 即呈指数
4增长。因此在两种生长机制之间有一个临界值, 通过图456) 时, 二维成核
生长机制从单核生长过渡到多核生长时的临界值大致分别在0. 3、0. [11]中得到的结果基本一致。我们注意到α=4, 生长速率比较快。α=6时的生长曲线在β=0. 35～0. 75, 3(d ) 和5(a ) 所示, 晶体生长质量较高,
图5 以二维成核机制在完整光滑平面上的晶体生长表面形态图
Fig . 5 The surface configurati on of crystal gr o wth on the perfect and s mooth p lane under the t w o 2di m ensi onal nucleati on
mechanis m (surface matrix is 60×60) :(a ) T wo 2di m ensi onal single nucleati on gr o wth mechanis m ;
(b ) T wo 2di m ensi onal multi p le nucleati on gr o wth mechanis m
3. 2 影响晶体生长速率的其它因素
3. 2. 1 表面扩散对晶体生长速率的影响
扩散是晶体表面存在的三种基本事件之一, 也是本文模拟中重点考查的对象。图2给出了平均扩散距离X S 对生长速率的影响。结果发现考虑扩散时的生长速率比不考虑扩散情况时要大, 扩散距离大小与生长速率成正比关系, 而且随着溶液过饱和度的增大扩散的影响更为显著。但整体而言, 表面扩散的影响并不引起生长机制的变化, 从图2中的生长曲线来看依然是线性关系。图6给出了在不同扩散距离下生长速率与热粗糙度之间的关系。可以清楚的看到, 在相同的溶液过饱和度下随着热粗糙度α的不断增大, 扩散对生长速率的影响逐渐减小。也就是说表面扩散对法向生长机制的影响更大, 而对二维成核机制尤其是单核生长机制的影响比较小。从图中还可以看到当α>4时生长速率趋缓, 可说明生长机制由法向生长过渡到了二维成核生长机制。
3. 2. 2 表面尺寸对晶体生长速率的影响
, 图7给整体来说在小的晶体表面上模拟得到尤其是在以单核机制生长时这种差别是非常明显的, 随着过饱和度的增大, 尺寸较小的晶体表面生长速率起伏开始增大。晶3. 2. 3 平均生长高度对晶体生长速率的影响
本文计算的晶体生长速率为平均生长速率, 我们发现相
同条件下晶体生长到不同台阶高度所对应的平均生长速率是
有差别的, 如图8所示, 生长机制以单核生长为主时, 这种因
生长到不同台阶高度而产生的平均生长速率的差别是非常小
的; 生长机制过渡到多核生长后, 随着过饱和度的增大, 平均
生长速率的差别开始越来越大。例如h =1. 1与h =2. 1时的
差别明显比h =2. 1与h =3. 1的差别要大。我们发现, 最初
在台阶高度为1的光滑平面上成核生长相对所花时间较长,
而在多核生长机制下台阶的扩展和新晶核的形成是同时的,
生长点急剧增多, 因此晶体生长速率会逐渐增大, 其平均生长
速率渐渐趋于一个稳定值。所以在模拟时我们尽量令晶体生图8 生长台阶的平均高度对生长速率的影响Fig . 8 The relati onshi p bet w een the gr o wth rate and mean gr o wth height
长到较高的同一台阶高度(例如h =3. 3) , 以便得到比较稳定的和具有可比性的平均生长速率。
图8中的生长曲线在β≥0. 5后以锯齿形上升, 原因就在于表面矩阵没有达到足够大而引起的生长速率的不稳定。这种现象也说明了计算机晶体生长本身具有的实验属性。因此要在考虑计算机运算速度和实际模拟所用时间的基础上尽量选取较大尺寸的晶体表面进行模拟。
4 结 论
(1) 作为在完整光滑晶面上进行的晶体生长模拟, 主要的晶体生长机制是二维成核生长。当然在热粗糙度比较低时也可以观察到法向生长机制, 而在过饱和度非常低时也发现了二维成核生长存在的生长死区;
(2) 通过模拟出的生长曲线结合晶体生长的表面形态可以进一步把二维成核生长机制细分为单核生长和多核生长。证实在热粗糙度比较大、过饱和度又比较低时, 线性的单核生长是主要生长机制;
(3) 表面扩散对晶体生长速率有一定影响。扩散优化了生长基元在晶体表面位置的排放, 使得晶体生长的表面能更小, 表面形态更趋稳定, 提高了晶体的生长速率。需要说明的是这种影响并不引起生长机制的变化, 而且随着热粗糙度的增大这种影响也越来越小;
(4) 计算机模拟的晶体表面尺寸对生长速率存在一定影响, 尤其是在单核生长时, 晶体生长速率随表面尺寸的增大明显增大。在多核生长机制下表面尺寸大时晶体生长速率趋缓, 相反在晶体表面尺寸小时生长曲线局部起伏较大。但表面大小对晶体生长曲线的整体走势没有影响;
(5) 晶体生长到不同平均高度时得到的平均生长速率是不同的, 在计算平均生长速率时, 需选定一个可比较的标准, 例如一定的生长时间或者本文采取的同一台阶平均高度, 而不是通常情况下简单的蒙特卡罗模拟步数的递加。
计算机晶体生长实质上是对晶体复杂生长过程的理解和反映。由于实际上不可能对微观的生长基元进行跟踪和定位, 因此注重模拟结果的真实具有重要意义。主要生长机制的晶体生长过程, 的应用前景, 可以进一步结合对电解质溶液的认识, , 采用E wald 求和等方法考虑长在第2, 建立了晶体生长的理论模型, 在此基础上我们给出:
(1) 初始化参量
h (x , y ) 表示晶体表面不同位置的高度, 初始值设为1, {x , y }表示晶体表面矩阵, x ∈(1, x max ) , y ∈(1, y max ) 。依据各点高度分别计算表面各位置的吸附配位数和脱附配位数, 以便得到k i , k i , k ij , K xy , K x , K 值。并预先设置基元生长个数h N, 基元个数n =0, 其中K x =+-∑K
y xy ,N 为每一生长层内的基元个数。
(2) 产生一个在(0, K ) 区间均匀分布的随机数ξ
(3) 随机选取某一点(x , y ) :若x
若y
(4) 选择事件:若ξ
否则若ξ
其它情形, 基元扩散, 取q =ξ-k i (x, y ) -k j (x, y )
(5) 选择扩散方向
若q
否则若q
令q =q -k i (x, y ) , j (x +1, y ) ;
否则若q
令q =q -k i (x, y ) , j (x, y -1) ;
否则若q
其它则h (x , y ) =h (x , y ) 。
以上寻找坐标点(x, y ) 以及选择事件的过程可以参见图9, 随机数按照标注1, 2, 3的顺序依次进行选择。+-+-
图9 一个蒙特卡罗循环步的流程图
Fig . 9 Fl ow chart of a Monte Carl o circulati on:(a ) fix a site (x, y ) ; (b ) choose a certain event
(6) 增加时间记数:t =t +1/K
(7) 刷新(x , y ) 及周围12个点的配位数, 12个点的配位数, 以便得到新的K xy , K x , K
(8) 若n
(9) 输出h (x , y ) , 以及生长速率等结果。
参[1]
[2]McPhers on A lexander, G, omolecular Crystal Gr owth as Revealed byA t om ic ForceM icr oscopy [J ].Journal of S tructural B iology , :32. Yu X L, W et al . I nterface Kinetics and Crystal Gr owth Mechanis m of a Ne w O rganometallic Coordinati on Compound:
Triallylthi ourea B r om ide [J ].Journal of C rystal Gro w th , 1997, 182:4282433.
[3]卢贵武, 李春喜, 等. K DP 晶体台阶生长动力学的激光干涉实验研究[J ].化学物理学报, 2003, 4:2892292.
[4]
[5]
[6]
[7]
[8]
[9]
[10]
[11]
[12]Andrea C L, M ir oslav Kotrla . Theory and Si m ulati on of Crystal Gr owth [J ].J. Phys . :Condens M atte . , 1997, 9:2992344. Cheng V K W et al . The Kinetic A sy mmetry bet w een Nucleati on Gr owth and D iss oluti on:a Monte Carl o Study [J ].Journal of C rystal Gro w th , 1989, 99(1) :2932299. Van EnckevortW J P, et al . I m purity B l ocking of Crystal Gr owth:a Monte Carl o Study [J ].Journal of C rystal Gro w th , 1998, 183(2) :4412448. Cheng V K W. A Monte Carl o Study ofMoving Step s during Crystal Gr owth and D iss oluti on [J ].Journal of C rystal Gro w th , 1993, 134(2) :3562365. Gil m er G H, Bennema P . Si m ulati on of Crystal Gr owth with Surface D iffusi on [J ].J. Appl . Phys . , 1972, 43:134721360. Lee Jae 2Wook, Hwang Nong M , Ki m Doh 2Yeon . Gr owth Mor phol ogy of Perfect and T winned Face 2centered 2cubic Crystal by Monte Carl o Si m ulati on [J ].Journal of C rystal Gro w th , 2003, 250:5382545. W ats on Grae me W , O liver PeterM , Parker Stephen C . A t om istic Si m ulati on of Crystal Gr owth at the a 〈100〉Scre w D isl ocati on Ter m inating at the {100}Surface of Mg O [J ].Surface Science , 2001, 474:1852190. Rak M ir osla wa, I zdebskiMarek, B r ozi Andrzej . Kinetic Monte Carl o Study of Crystal Gr owth fr om Soluti on [J ].Co m p . Phys . Co mm . , 2001, 138:2502263. Malkin A J, Kuznets ov Yu G, Land T A, et al . Mechanis m s of Gr owth f or Pr otein and V irus Crystals [J ].N ature S tructural biology , 1995, 2
(11) :9562959.
[13]张克从, 张乐惠, 主编. 晶体生长科学与技术(第二版) 上册[M].北京:科学出版社, 1997.
[14]
[15]
[16]Burt on W K, Cabrera N, Frank F C . The Crystal of Crystals [J ].Phil. Trans . , 1951, A243:2992358. De Yoreo, Land T A, Dair B. Gr owth Mor phol ogy of V icinal H ill ocks on the {101}Face of KH 2P O 4:fr om Step 2fl ow t o Layer 2by 2layer Gr owth [J ].Physical Revie w L etters , 1994, 73:8382841. Malkin A J, Kuznets ov Y G et al . I n Situ A t om ic Force M icr oscopy Studies of Surface Mor phol ogy, Gr owth Kinetics, Defect Structure and
D iss oluti on in Macr omolecular Crystallizati on [J ].Journal of C rystal Gro w th , 1999, 196:4712488.