分子动力学实验报告

分子动力学实验报告

实验名称平衡晶格常数和体弹模量

实验目的

1、学习Linux系统的指令

2、学习lammps脚本的形式和内容

实验原理

原子、离子或分子在三维空间做规则的排列，相同的部分具有直线周期平移的特点。为了描述晶体结构的周期性，人们提出了空间点阵的概念。为了说明点阵排列的规律和特点，可以在点阵中去除一个具有代表性的基本单元作为点阵的组成单元，称为晶胞。晶胞的大小一般是由晶格常数衡量的，它是表征晶体结构的一个重要基本参数。

在本次模拟实验中，给定Si集中典型立方晶体结构：fcc，bcc，sc，dc。根据

可判定dc结构是否能量最低，即是否最稳定

材料在弹性变形阶段，其应力和应变成正比例关系（即符合胡克定律），其比例系数称为弹性模量。弹性模量是描述物质弹性的一个物理量，是一个总称，包括杨氏模量、剪切模量、体积模量等。在弹性变形范围内,物体的体应力与相应体应变之比的绝对值称为体弹模量。表达式为

B=?

dP dV V

式中，P为体应力或物体受到的各向均匀的压强，dV V为体积的相对变化。对于立方晶胞，总能量可以表示为ε=ME，E为单个原子的结合能，M 为单位晶胞内的原子数。晶胞体积可以表示为V=a3，那么压强P为

P=?dε

dV

=?

M

3a2

dE

da

故体积模量可以表示为

根据实验第一部分算出的平衡晶格常数，以及能量与晶格间距的函数关系，可以求得对应晶格类型的体积模量。并与现有数据进行对比。

实验过程

（1）平衡晶格常数

将share文件夹中关于第一次实验的文件夹拷贝到本地，其中包含势函数文件和input文件。

$ cp□-r□share/md_1□.

$ cd□md_1

$ cd□1_lattice

通过LAMMPS执行in.diamond文件，得到输出文件，包括体系能量和cfg文件，log文件。

$ lmp□-i□in.diamond

用gnuplot画图软件利用输出数据作图，得到晶格长度与体系能量的关系，能量最低处对应的晶格长度即是晶格常数。

Si为diamond晶格结构时晶格长度与体系能量关系图如图，

由图可得能量最小处对应取a0=5.43095。

Si为fcc晶格结构时晶格长度与体系能量关系图如图，

a0=4.15。

改写后的sc、bcc脚本文件分别如图所示

Si为sc晶格结构时晶格长度与体系能量关系图如图，

a0=2.6。

Si为sc晶格结构时晶格长度与体系能量关系图如图，

a0=3.25。

（2）体弹模量

利用gnuplot的二次拟合功能，对上述步骤中得出的晶格常数与体系能量关系图进行二次拟合，

f(x)=a+b*x+c*x**2

fit□f(x)□‘data’□via□a,b,c

得到各晶格结构下的二次拟合函数表达式中的系数a，b，c，如图

Si在diamond结构下的拟合结果

Si在fcc结构下的拟合结果

Si在sc结构下的拟合结果

Si在bcc结构下的拟合结果

结果与分析

（1）平衡晶格常数

由实验过程中所得的各个关系图可知：

Si在diamond结构下，平衡晶格常数为5.43095×10-10m，对应能量约为-936.706eV；

Si在fcc结构下，平衡晶格常数为4.1×10-10m，对应能量约为-420eV；

Si在sc结构下，平衡晶格常数2.61×10-10m，对应能量约为-109.65eV；

Si在bcc结构下，平衡晶格常数为3.245×10-10m，对应能量约为-219eV；

从以上对比可得，Si能量为-936.706eV时为最小值，即对Si来说，diamond 结构最稳定，且此时平衡晶格常数为5.43095×10-10m。

（2）体弹模量

。由此公式计算得：

由体弹模量计算公式得B=Mc×160.22

9a0

Si在diamond结构下，M=8，a0=5.43095，c=417.69，得B=10953.23684GPa；

Si在fcc结构下，M=4，a0=4.1，c=2122.6，得B=36865.36282GPa；

Si在sc结构下，M=1，a0=2.61，c=201.149，得B=1371.99203GPa；

Si在bcc结构下，M=2，a0=3.245，c=359.335，得B=3942.65733GPa；

实验结论

通过本次实验，我们在计算不同结构下的单晶硅的平衡晶格常数与其对应能量，发现单晶硅在diamond结构下最为稳定。在diamond结构下，单晶硅的平衡晶格常数为5.43095×10-10m，对应能量约为-936.706eV，体弹模量为10953.23684Gpa。在Si为diamond晶格结构时晶格长度与体系能量关系图中有一条首尾相连的直线，可能为在课余时间多次计算而产生的。在不同结构下体弹模量的误差较大，可能为拟合过程中产生的了一定的误差，应该尝试更高次数的拟合方程。

实验名称晶体点缺陷

实验目的

利用分子动力学方法构建和研究fcc 晶体Cu 的点缺陷，空位和间隙原子，并利用以上公式计算空位形成能和间隙原子形成能。

实验原理

晶体中的缺陷包括从原子、电子水平的微观缺陷到显微缺陷。按几何形态来分类，点缺陷是零维缺陷，这种缺陷在各个方向上的延伸都很小，仅发生在晶格中的一个原子尺寸范围，如空位、间隙原子、置换原子。点缺陷在晶体中呈随机、无序的分布状态。

点缺陷普遍存在于晶体材料中，它是晶体中最基本的结构缺陷，对材料的物理和化学性质影响很大。

根据点缺陷相对于理想晶格位置可能出现的几种主要偏差状态，可将其命名如下：

（1）空位：正常节点位置上出现的原子空缺。

（2）间隙原子（离子）：指原子（离子）进入正常格点位置之间的间隙位置。

（3）杂质原子（离子）：晶体组分意外的原子进入晶格中即为杂质，杂质原子若取代晶体中正常格点位置上的原子（离子）即为置换原子（离子），也可进入正常格点位置之间的间隙位置而成为填隙的杂质原子（离子）。

常见的空位、间隙原子都是构成晶体的原子或离子偏离原有格点所形成的热缺陷。在一定温度下，晶体中各原子的热振动状态和能量并不同，遵循麦克斯韦（Maxwell）分布规律。热振动的原子某一瞬间可能获得较大的能量，这些较高

能量的原子可以挣脱周围质点的作用而离开平衡位置，进入到晶格内的其他位置，于是在原来的平衡格点位置上留下空位。根据原子进入晶格内的不同位置，可以将缺陷分为弗伦克尔（Frenkel）缺陷和肖特基（Schottky）缺陷。

空位和间隙原子都只有一个原子大小的尺度，因此很难通过实验对其进行直接的观察。通过场离子显微镜可分辨金属表面上的原子排列而直接观察到金属表层中的空位位置。利用电子显微镜薄膜投射法可观察到空位片或间隙原子片，但实验方法研究缺陷时利用较多的还是缺陷对晶体性质的影响。例如，通过测量晶体的膨胀率和电阻率的变化规律，即可对点缺陷的存在、运动和相互作用等方面展开间接的研究。分子动力学方法对金属材料原子尺度物理和化学过程的研究具有实验法无法比拟的优势，可直观的模拟和分析晶体中的点缺陷。若我们搭建完整晶体的原子个数为N，能量为E1，通过删除和增加一个原子得到空位和间隙原子，充分弛豫后体系能量为E2，则空位形成能E v和间隙原子形成能E i分别为：

实验过程

进入实验2的目录，运行引入空位原子的执行文件$ lmp□-in□in.vacancy

查看实验输出的cfg文件

$ A.i686□v20.cfg

$ A.i686□v26.cfg

所得空位原子如图，

查看引入空位原子后的输出data，

$ cat□data.v

引入空位后所输出的能量数据如图，

同理，间隙原子与引入间隙原子后所输出的能量数据如图，

尝试修改脚本文件，研究空位浓度对其形成能的影响，改变模拟体系的大小。修改后的脚本文件如图，

所得结果如图，

尝试修改脚本文件，研究间隙原子浓度对其形成能的影响，改变模拟体系的大小。修改后的脚本文件如图（仅显示修改部分），

所得结果如图，

结果与分析

E1，计算得空位的由原始的空位原子能量数据，根据公式E v=E2?N?1

N

形成能为1.257eV；而在修改空位浓度后，空位形成能为1.256eV。

E1，计算得间隙原由原始的间隙原子能量数据，根据公式E i=E2?N+1

N

子的形成能为3.298eV；而在修改间隙原子浓度后，间隙原子形成能为

3.287eV。

实验结论

通过本次实验，在比较不同空位浓度下的空位形成能可以发现，空位浓度越高，空位形成能越高。

同理，在比较不同间隙原子浓度下的形成能可以发现，间隙原子浓度越高，间隙原子的形成能也越高。

实验名称：刃位错和螺位错

实验目的

1、观察刃位错和螺位错的应力应变场的分布

2、学会计算单位位错的位错能。

实验原理

晶体的线缺陷表现为各种类型的位错。位错的概念最早是由Taylor、Orowan 和Polanyi 在研究晶体滑移过程时提出来的，他们认为，晶体实际滑移过程并不是滑移面两边的所有原子都同时做整体刚性移动，而是通过在晶体存在着的称为位错的线缺陷来进行的，位错在较低应力的状态下就开始移动，使滑移区逐渐扩大，直至整个滑移面上的原子都发生相对位移。位错有几种基本的类型，首先是螺型位错。

1．螺型位错

一个螺型位错标志着晶体中滑移区与未滑移区之间的分界线，这一边界与滑移方向平行。可以想象用刀子将晶体切开一部分，然后使之平行于切割边界切移一个原子间距，这样造成一个螺型位错，螺型位错使相继的原子平面变成一个螺旋曲面，这就是这类位错名称的由来。

要准确地对晶体中位错周围的弹性应力场进行定量计算，是复杂而困难的。为简化起见，通常可以采用弹性连续介质模型来进行计算。该模型首先假设晶体是完全弹性体，服从胡克定律；其次，把晶体看作是各向同性的；

第三，近似地认为晶体内部是由连续介质组成，因此晶体中的位移、应力和应变等量是连续的，可用连续函数表示。设想有一各向同性材料的空心圆柱体，先把圆柱体沿XZ面切开，然后使两个切开面沿Z 方向做相对位移b，再把这两个面胶合起来，这样就相当于形成了一个伯氏矢量为b 的螺型位错。圆柱体只有Z方向的位移，即Ux=Uy=0。Uz 在切边处不连续，因此

假设在各向同性的介质中，UZ随着θ 角均匀的增大，可以得到位移UZ

与θ 和r 的关系如下

此即为螺位错的位移场公式。也是在本次实验中，我们用来构造螺位错的依据。首先，我们搭建一个完整晶体，以中心处为位错的核心，然后根据位移场公式相继移动体系内的每个原子，使其符合螺位错的位移场分布。经过能量最小化后，我们可以得到一个稳定的含有螺位错的构型。由于圆柱体只有沿z 方向的位移，因此只有切应变。相应的，各应力分量为

2．刃型位错

滑移区与未滑移区之间的边界就称为位错。这个位错的位置由挤入上半部分晶体的额外垂直半原子面的边缘标志。在位错附近，晶体的形变可以看作是由于在晶体上半部分插入了一片额外的原子面所产生。这个原子面的插入使上半部分晶体中的原子受到挤压，而使下半部分晶体中的原子受到拉伸。

按照弹性力学理论可以得到，刃型位错的位移场诸分量为：

相应的，各应力分量为：

3．位错的应变能

位错的能量可以分为两部分：位错中心畸变能Ec 和位错应力场引起的弹性应变能Ee。位错中心区域由于点阵畸变很大，不能用胡克定律，而需借助于点

阵模型直接考虑晶体结构和原子间的相互作用。据估算，这部分能量大约为总应变能的1/10～1/15 左右，故常予以忽略，而以中心区域以外的弹性应变能代表位错的应变能，此项能量可以采取连续介质弹性模型根据单位长度位错所作的功求得。对于刃型位错，

同理，对于螺型位错有

实验过程

（1）刃位错

进入刃位错对应的文件夹中，运行刃位错的执行文件，

$ lmp_serial□-i□in.edge

查看运行完毕后输出的cfg文件，

$ A.i686□198.cfg

对所得结构图用切应力染色，得到图，

进行柏氏矢量标定后如图所示，

分别采用势能及stress_1到stress_6染色，观察刃位错的能量及各力场分布，得到图

从运行刃位错的结果得知刃位错能量为E=-220806.385。运行

= -220827.838。利用gedit查看edge.data文件，in.perfect文件，得到E

得到刃位错长度为L=16.152×10-10m。

（2）螺位错

进入螺位错对应的文件夹中，运行刃位错的执行文件，

$ lmp_serial□-i□in.screw

查看运行完毕后输出的cfg文件，

$ A.i686□28.cfg

对所得结构图分别用势能、切应力及stress_1到stress_6染色，得到

图

从运行螺位错的结果得知螺位错能量为E=-355522.615。运行in.perfect = -355536.927。利用gedit查看screw.data文件，得到螺位错文件，得到E

长度为L=9.891×10-10m。

结果与分析

中，计算可得刃位错的单位长度形成能为：将实验所得各数据带入公式γ=γe

L

γ=1.328eV/?=2.127×10-9J/m

同理计算可得螺位错的单位长度形成能为：

γ=1.447eV/?=2.318×10-9J/m

由实验教程所给的公式计算得刃位错形成能理论值为：1.14×10-9J/m，螺位错形成能理论值为：1.67×10-9J/m

实验结论

通过本次实验，我们可以直观的看出刃位错与螺位错结构不同，刃位错柏氏矢量垂直于滑移线，而螺位错的柏氏矢量则平行于滑移线。由应力场的分布可看出刃位错在xy、xz方向的应力较小，而螺位错在xx、xy方向上应力较小。在能量方面，螺位错的形成能比刃位错的形成能高，但实验值与理论值之间相差较大。