【LAMMPS如何系列】计算平衡晶格常数

平衡晶格常数（equilibrium lattice constant）对应的体系能量是最低的，因此只需要计算一些列不同晶格常数下体系的能量，那么体系能量最小时对应的晶格常数就是平衡晶格常数。

在初次计算时，建议将晶格常数增加的步长取得大一点，比如0.1A，这样通过结果数据就可以大致确定晶格常数的范围。

然后再将晶格常数的范围缩小，减小步长，比如0.01A，这样做就可以获得精度比较高的结果。

下面以金刚石为例介绍具体如何计算平衡晶格常数。

输入文件： in.diamond

# This input script is used to calculate
# the lattice constant of diamond
# Powered by Xianbao Duan
# Email: xianbao.d@gmail.com
# Website: http://www.52souji.net/

units			metal
boundary        p p p
atom_style		atomic

variable		i loop 20
#variable		x equal 2.5+0.1*$i
variable        x equal 3.5+0.01*$i

# build the model
lattice			diamond $x
region			box block 0 10 0 10 0 10
create_box		1 box
create_atoms	1 box

# specify the potential
pair_style		tersoff
pair_coeff		* * SiC.tersoff C
mass			1 12

variable		n equal count(all)
variable        P equal pe/$n
#variable		v equal vol

timestep		0.005
thermo			10

# minimize the total energy
min_style		cg
minimize		1.0e-12 1.0e-12 1000 1000

print			"@ $x $P"

# loop
clear
next 			i
jump			in.diamond

脚本很简单，只要一个简单的能量最小化就可以了，然后输出晶格常数和结合能log文件中。

然后从log文件中将晶格常数和结合能的数据提取出来。

$ grep ^@ log.lammps > lat.vs.Ecoh.step1

grep是一个linux命令，上面的意思就是将log.lammps文件中以@开始的行输出到 lat.vs.Ecoh.step1 文件中。

得到的文件中，每一行开头会有一个@，通过列编辑器一起删掉，就可以获得有如下内容的文件。

3.51 -7.35079243
3.52 -7.356562846
3.53 -7.361141466
3.54 -7.364556691
3.55 -7.366836396
3.56 -7.368007936
3.57 -7.368098161
3.58 -7.367133418
3.59 -7.365139565
3.6 -7.362141974
3.61 -7.358165545
3.62 -7.353234707
3.63 -7.347373432
3.64 -7.340605242
3.65 -7.332953213
3.66 -7.324439984
3.67 -7.315087768
3.68 -7.304918353
3.69 -7.293953116
3.7 -7.282213024

上面的数据实际上对应的就是结合能曲线，其能量最低点对应的就是平衡晶格常数。当然，能量最低点的能量就是晶体的结合能。

具体如何从上面的数据中获得最终结果，可以参考我之前写的一篇文章：MATLAB计算平衡晶格常数。

晶格常数表可以参考：所有元素的晶格常数查询表。