lammps实例(3)

Project #1

硅的晶格常数和体弹模量的计算

一、平衡晶格常数和内聚能

自然条件下硅为金刚石结构(dc )。计算模拟时，我们可以假定它为各种结构，f cc, bcc, sc, dc. 可以预测，模拟的dc 结构的硅的体系能量最低，也即最稳定。下面我们将运用LAMMPS 来对硅的各种结构进行模拟。

定义晶格能量为Φ， 数密度为 ρ：

pot

E N Φ=

N V

ρ= 其中E pot 为势能，

N 为体系总原子数，V 为体系的体积。选取 Stillinger-Weber (SW)，以下面命令执行 lammps 运算：

其中，lmp_serial 为 lammps 命令；”<” 符号为读取符；in.Silicon 为输入文件，里面包含运算所需要的各种数据和命令；-log 指定输出文件的名称。

可以看到屏幕上显示出lammps 运行的信息。这个计算量很小，所以很快就结束。接下来以如下命令来查看计算得到的数据：

grep 是linux 中一个很重要的命令，用来搜索文本，读取匹配的行并打印出来。这里是搜索 dc.log 文件，将 @ 开头的行打印出来。如下：

晶格参数为5.4305埃,数密度为0.0499540303，每个原子的能量为-4.336599609eV.

下面具体来看刚才给的输入文件，in.Silicon .

dc.log 文件中有原子总数的信息，

每个金刚石晶胞中有8个原子，383216⨯=，所以是216个原子。如下给出各种结构下的体系的原子数：

晶体结构类型

晶胞中的原子数 总原子数 简单立方SC

1 27 体心立方BCC

2 54 面心立方FCC

4 108 金刚石DC 8 216

表1.

不同晶体结构中的原子数

下图是计算模拟得出的各种结构下的数密度与每个原子能量的关系图。

横坐标为数密度， 以金刚石为例，ρ= 8/5.4315^3=0.049926，也即我们直接通过 grep 命令得到的第二项值；纵坐标为每个原子的能量，为第三项值。

金刚石之外，还需计算其他结构。只需对 in.Silicon 做稍微改动：

首先，将in.Silicon 复制成in.fcc :

然后编辑 in. fcc

改动如下几项：

然后如下命令执行：

相应的，如下命令查看log 文件中的数据：

以同样方法编辑in.bcc, in.sc，计算不同晶格参数时的体系能量值，并绘制下图：

图1. 不同结构下的硅的晶格能。可以看出金刚石结构对应最低能量，最为稳定

下图更为细致地画出金刚石结构中，不同晶格参数所对应的内聚能。内聚能(cohesive energy E coh)的定义是，最小的晶格能。由图可以得到，平衡晶格常数为a0 = 5.431 Å，内聚能为E coh = −4.3365 eV.

图2. 金刚石结构中的晶格能 VS 晶格参数。五阶拟合得到平衡晶格常数5.43095(Å).

2345E 104799.5791127112.82866a 15604.99601a 7592.99351a 1134.25198a 57.82091a =-++-+-

二、体弹模量

我们同时可以从晶格能曲线在最低处得到体弹模量的信息。体弹模量定义为：

/≡-dP B dV V

V 和 P 分别为晶胞的体积和体系的压强。

我们已经得到了内聚能与晶格参数的函数关系，对于立方晶胞而言，

23ε=-

=-d M dE P dV

a da 由此， 022

09=a M d E B a da

其中 a 0 为平衡晶格常数，M 为体积为3=V a 的晶胞中的原子数目。

从多项式拟合，可以得到 0

22a d E da , 此例中为 3.87081 eV/ Å2 (a 0 =5.431 Å). 由上面

公式，计算得到硅的体弹模量为 B = 101.366 GPa. 文献中的实验数据为 99 GPa.

可视化

每一次 lammps 运行后，会生成一个 dump.atom 文件。可以通过如下命令转换为 Atomeye 可读取的 .cfg 文件：

这时当前目录中就会生成 001.cfg 之类的 .cfg 文件，然后通过 Atomeye 查看： Tab 键可以改变观察方位

上下左右箭头键可以转动原子

PgUp 和 PgDn 改变原子大小

Alt + 1 和 Alt + 2 以及 Alt + 3 改变原子颜色

先按一次9这个键，转动时就是每次以90度转动