Lammps 石墨烯实例

一、简介
1.SiC热分解制备石墨烯
自2004年Novoselov、Geim和合作者们从石墨上剥离出世界上第一种二维材料——单层石墨：石墨烯（Graphene）以来，石墨烯就受到了科技界的广泛重视[1]。

Novoselov 和Geim两人因此在2010年获得了诺贝尔物理学奖。

因为石墨烯的独特特性，在许多技术领域例如光电子学上它都被寄予厚望。

研究石墨烯这种材料相关的物理化学特性和发展大面积、高质量生长石墨烯的技术，同时将其与器件物理学联系起来是我们研究和应用石墨烯的必由途径。

石墨烯是由碳元素组成的二维六边形材料，其在光学、电学、热学、力学等性质十分优异。

它有可能在后摩尔定律时代成为硅（Silicon）的继任者，在单分子气体传感器[2]、自旋电子学[3]、量子计算[4]、太赫兹振荡器[5]等等领域发挥重要作用。

如今，从石墨上剥离出石墨烯仍然是一种重要的石墨烯制备方方法。

然而，这种方法产生的石墨烯大小通常不超过1000 μm2，只适合实验室研究，尚不能在工业上大规模应用。

科学家发展了其他的石墨烯制备方法，包括将石墨烯视作一种薄膜来生长的化学气相沉积（Chemical Vapor Deposition, CVD）法、热分解碳化硅法（SiC thermal decomposition）、氧化石墨烯还原法（Graphene oxide reduction）等。

CVD法通过使含碳气源在有催化作用的金属表面分解或者使溶入到这些有催化作用的金属中的碳（C）发生表面偏析，使得在金属表面生成石墨烯或者多层石墨烯（Few-Layer Graphene, FLG）。

能否直接在半导体/绝缘体上生长石墨烯呢？碳化硅热分解成功的解决了这一问题。

最早试图使六方晶系的SiC晶体石墨化的研究报告见于1961年，Badami在高温和真空环境下得到了发生了一定石墨化的SiC[6]。

在一定的退火条件下，SiC晶体表面发生热分解，Si原子发生解吸附，而C原子留下来重新排列和组合可以生长成外延型的石墨烯层[7]。

更细致的研究发现用热退火的方法在六方SiC的Si面上生长的石墨烯比C面有更好的可控性，例如：可以更好的控制石墨烯的层数。

Si面上生长的石墨烯生长方向与基底晶体结构有密切关系，这样提供了在基底上均匀覆盖和特定方向生长石墨烯的可能性。

特别地，石墨烯直接生长在半导体SiC上使得我们无
须将石墨烯转移到其他绝缘体上就可直接使用。

Huang 等根据STM观察到的现象提出了SiC热分解法生长石墨烯的重要模型：Bottom-up生长模型[8]。

SiC刚开始热分解的过程中Si原子发生解吸附，剩下的C原子发生（6√3×6√3）R30º的结构重排。

虽然，这层C原子层与石墨烯有着类似的原子排布，但它并不是石墨烯层。

原因是这层C原子层与下层SiC发生较强的共价耦合作用，它已经失去了石墨烯的诸多优异特性。

当热分解的过程继续，下层的Si-C键继续断裂，Si原子继续挥发掉，下层的C原子继续发生（6√3×6√3）R30º结构重排。

此时，先前生成的与SiC结合在一起的C原子层被新形成的C原子层顶了上去，脱离了SiC基底的作用。

这层C原子层在结构和物理化学性质上与从石墨上剥离出的石墨烯较为相似，这样就生成了单层石墨烯。

如果SiC热退火的过程继续，以上描述的过程还会生成新的过渡层，从而第一层过渡层现在就变成了双层石墨烯的上层。

Bottom-up生长模型较好的解释了目前SiC热退火生长石墨烯的过程。

2.SiC介绍
与单晶Si类似，SiC是一种Si与C以共价键相键结而成的半导体化合物。

在这种晶体中，每1个Si与周围4个C以共价键连接形成四面体，在这种SiC晶体中存在两种四面体结构(图1)。

第一种SiC结构可以通过另一种SiC结构沿c轴旋转180°得到，c轴指Si-C双原子层的法向方向。

在每一层Si原子(或C原子)上，原子排列都是六方密排结构。

在垂直c轴的Si-C双层原子上，共有3种不同的占位，分别为A,B,C。

图1 SiC的两种四面体结构
以上述3种占位的双原子层分别沿c轴以一定周期排列形成特定的SiC同构异形体。

例如，在3C-SiC中，双原子层的堆垛次序为ABCABCABC…其中字母“C”指的是立方晶体结构(Cubic Crystal Structure)，数字“3”指的是在一个重复单元中双原子层的数目(ABC)。

3C-SiC是SiC中唯一的一种立方晶体结构，又称之为β-SiC。

SiC晶
体同构异形体中大多数属于六方密排堆积结构。

最常见的六方密排结构是2H,4H和6H 堆垛方式，其堆垛次序分别为ABAB…, ABCBABCB…和ABCACBABCACB…类似的字母“H”指的是六方密排结构(hexagonal crystal structure)，数字“6”指的是在一个重复单元中双原子层的数目(ABCACB)。

同时也存在其他SiC晶体构型，例如15R-SiC 中，字母“R”指的是菱形晶体结构(rhombohedral crystal structure)。

值得注意的是，在堆垛次序中每一个Si原子层(或C原子层)都存在一个局部立方晶格(k)或六方密排(h)类型环境。

例如在3C-SiC中仅存有k类型的原子，2H-SiC仅存有h类型的原子。

在4H-SiC和6H-SiC中，k和h原子数比分别为1:1和2:1[9]。

(图2和图3)
图2 沿着c轴堆垛方向的三种(A,B,C)占位示意图。

c轴垂直纸面方向。

图3 在最常见的三种SiC构型中Si–C 双原子层堆垛方式
MMPS介绍
LAMMPS是一个经典的分子动力学代码，可以模拟液体中的粒子，固体和气体的系综，也可以采用不同的力场和边界条件来模拟全原子，聚合物，生物，金属，粒状和粗料化体系。

LAMMPS可以计算的体系小至几个粒子，大到上百万甚至是上亿个粒子。

LAMMPS可以在单个处理器的台式机和笔记本本上运行且有较高的计算效率，但是它是专门为并行计算机设计的。

他可以在任何一个按装了C++编译器和MPI的平台上运算，这其中当然包括分布式和共享式并行机和Beowulf型的集群机。

LAMMPS是一个可以修改和扩展的计算程序，比如，可以加上一些新的力场，原子模型，边界条件和诊断功能等[10]。

二、模型和模拟方法
1.经验势的选取
为了使用lammps分子动力学模拟SiC衬底，我们首先需要描述C和Si之间的相互作用，在本例中我们采取了最为流行的半经验Tersoff势。

共价键分子中原子间的相互作用势不仅取决于原子间的距离, 而且与原子间的成键方向有密切联系, 因此在计算共
价键原子间的相互作用时必须考虑周围原子的影响。

在Tersoff势或TEA势中经验势的一般表达方式为
E=
X
i E i=
1
2
X
j=i
V(r ij)
其中E i为在i位置处一个原子的总势能，
V(R ij)=f C(r ij)[a ij V R(r ij)+b ij V A(r ij)]
上式f c是一个从1到0的连续变化的截断函数。

它使两原子间的相互作用在距离较远的时候趋于零。

我们由Tersoff势选择截断函数f C为：
f C(r ij)=8
><
>:
1;r ij<R ij
1
2
+1
2
cos(¼(r ij¡R ij)=(S ij¡R ij));R ij<r ij<S ij 0;r ij>R ij
上式i和j分别表示C或Si, r ij是ij的键长，S ij=(S i S j)1=2和R ij=(R i R j)1=2分别为原子范围的最大和最小。

2.构建6H-SiC衬底
为了模拟6H-SiC衬底外延石墨烯，我们首先需要知道SiC衬底上C原子的位置。

另外为了实现SiC在真空环境中Si的蒸发，我们还需要在SiC表面设计一个富碳的环境。

根据6H-SiC晶体资料，我们知道6H-SiC属于六方密排结构，晶格参数和键角的基本信息为a=b=c；®=¯=90±；°=120±。

每个6H-SiC单元有12个原子(6个Si，6个C)，我们通过Avogadro软件导出单个6H-SiC每个原子的坐标，原子类型信息。

如图4。

图4 单个6H-SiC晶胞各原子位置(来自Avogadro)
由于单个6H-SiC晶胞属于三斜晶系，我们由Avogadro导出的数据可以计算出相应的晶格常数：
a=b=3:095c=15:17
®=¯=90±°=120±
在lammps软件中，默认是使用正交模拟盒子，通过定义原点和三组棱向量得到
a=(xhi-xlo;0;0)
b=(0;yhi-ylo;0)
c=(0;0;zhi-zlo)
lammps也提供了一组三斜(非正交)模拟盒子来定义平行六面体结构，定义原点和
三组棱向量[11]：
a=(xhi-xlo;0;0)
b=(xy;yhi-ylo;0)
c=(xz;yz;zhi-zlo)
三斜晶系结构定义使用了晶格常数a b c和三个角度®¯°。

在lammps中定义盒子大小(lx,ly,lz) = (xhi-xlo,yhi-ylo,zhi-zlo)和倾斜因子(xy,xz,yz)，这些参数的具体关系如下：
lammps构建样品分为内部搭建模型和内部搭建模型。

本文使用外部搭建模型，并
使用replicate命令搭建6H-SiC样品，
replicate 12 12 1
构建了12×12×1超晶格结构。

接着采用周期性边界条件，
boundary p p p
构建样品的一个关键步骤就是在样品表面创建一个10nm厚的真空层，防止重复结构对其造成影响。

3.设计富碳层
为了模拟6H-SiC高温Si原子蒸发，剩下的C原子重新自组装形成石墨烯。

我们在Avogadro软件导出原子的坐标和类型的数据进行修改。

首先沿着Z轴方向，把最上层的Si原子删掉，然后修改为C原子，xyz坐标不变。

图5 12 x 12 x 1 的6H-SiC超晶格结构(俯视图Yellow: Carbon Blue: Silicon)
图6 12 x 12 x 1 的6H-SiC超晶格结构(侧视图Yellow: Carbon Blue: Silicon )
图7 12 x 12 x 1 的6H-SiC超晶格结构(倾斜角度Yellow: Carbon Blue: Silicon) 4.模拟退火的流程
在构建样品完毕后实施退火程序，按照参考文献退火方式先从1K升温到300K并进行样品弛豫，然后以一定的升温速率达到退火温度，并在退火温度弛豫60000步，最后以一定的降温速率降到1K的温度。

图8 基于文献的退火程序
图9 基于文献的能量最小化模式
三、结果和讨论
1.单层石墨烯的生长
图10 模拟退火实际升温曲线
图11 6H-SiC表面C原子在1400K后变化图
图12 6H-SiC表面C原子在1600K后变化图
图13 6H-SiC表面C原子在1800K后变化图
图13 6H-SiC表面C原子在2000K后变化图
图13 6H-SiC表面C原子在2500K后变化图
2.结论
由上面lammps模拟结果，我们可以发现随着退火温度的不断提高，由1400K变化到2500K在6H-SiC上表面，Si原子不断蒸发后，形成一个富碳的环境，在高温的驱动下C原子重新自组装聚集。

在本次模拟中，虽然没有看到类似文献中大量成片石墨烯的形成是从温度的不断升高[12]，我们可以清楚的发现温度越高C原子聚集程度越高，而且离石墨烯结构的形成越接近，最后在1800K以上发现存在石墨烯六元环的趋势，并2500K发现存在石墨烯六元环存在。

从文献可知若修改相互作用势为更加精细的TEA势后，石墨烯趋势会更加明显。

参考文献
[1]K.S. Novoselov, A.K. Geim, S.V. Morozov, D. Jiang, Y. Zhang, S.V.Dubonos, et al.,
Electric Field Effect in Atomically Thin Carbon Films, Science, 306 (5696) (2004) 666.
[2] F. Schedin, et al., Detection of individual gas molecules adsorbed on graphene, Nature
Mater, 6 (2007) 652.
[3]Y.W. Son, M.L. Cohen, S.G. Louie, Half-metallic graphene nanoribbon, Nature, 444
(2006) 347.
[4]V.I. Fal'ko, Graphene Quantum information on chicken wire, Nature Phys, 3 (2007) 151.
[5] F. Rana, Graphene terahertz plasmon oscillators, IEEE Trans. Nanotech., 7 (2008) 91.
[6] D.V. Badami, Graphitization of α-Silicon Carbide, Nature, 193 (1962) 569.
[7]van A.J. Bommel, Van J. Crombeen and A. Tooren, LEED and Auger electron
observations of the SiC(0001) surface, Surf Sci, 48 (1975) 463.
[8]H. Huang, et al., Bottom-up growth of epitaxial graphene on 6H-SiC(0001), ACS Nano.
2 (2008) 2513.
[9]Fan, Jiyang, and P. K. Chu. General Properties of Bulk SiC. Silicon Carbide
Nanostructures. 2014:7-114.
[10]/doc/Section_intro.html
[11]/doc/Section_howto.html#triclinic-non-orthogonal-
simulation-boxes
[12]Yoon, T. L., Lim, T. L., Min, T. K., Hung, S. H., Jakse, N., & Lai, S. K. (2013). Epitaxial
growth of graphene on 6h-silicon carbide substrate by simulated annealing method.
Journal of Chemical Physics, 139(20).
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69。