堆垛层错能的计算

# 模型初始设置

units       metal

dimension   3

boundary    p p s

atom_style  atomic

variable latparam1 equal 3.52

variable x_displace equal -1*(${latparam1}/sqrt(6))

variable z_displace equal -1*(${latparam1}/sqrt(3))

variable xdim equal ${latparam1}*sqrt(6)/2*10

variable ydim equal ${latparam1}*sqrt(2)/2*10

#Ni的内聚能

variable Ecoh equal -4.45026

#晶格常数

lattice     fcc ${latparam1}

#设置三个区域

region      1 block -.001 ${xdim} -.001 ${ydim} -.001 15.4839 units box

region      2 block -.001 ${xdim} -.001 ${ydim} 15.4839 17.5161 units box

region      3 block -.001 ${xdim} -.001 ${ydim} 17.5161 36.5808 units box

region      whole block 0 ${xdim} 0 ${ydim} 0 36.5808 units box

#创建box

create_box  3 whole

#设置晶格取向，并在三个区域内填充原子

lattice fcc ${latparam1} orient x 1 1 2 orient y -1 1 0 orient z -1 -1 1

create_atoms    1 region 1

lattice fcc ${latparam1} orient x 1 1 2 orient y -1 1 0 orient z -1 -1 1

create_atoms    2 region 2 

lattice fcc ${latparam1} orient x 1 1 2 orient y -1 1 0 orient z -1 -1 1

create_atoms    3 region 3 

#力场设置

pair_style  eam/alloy  

pair_coeff * * FeCuNi.eam.alloy Ni Ni Ni

#计算势能

compute peratom all pe/atom

compute eatoms all reduce sum c_peratom

#设置热力学信息输出

thermo 1

thermo_style custom step pe c_eatoms

#原子轨迹输出

dump 1 all atom 1 dump.xyz

run 0

#计算原子个数

variable N equal count(all)

#初始原子个数

variable No equal $N

#初始能量

variable E equal "c_eatoms"

variable Eo equal $E

#对原子进行分组

group bot region 1

group middle region 2

group top region 3

#删除一层原子

delete_atoms group middle

#记录当前原子个数

variable Nf equal $N

#上层原子向下移动

displace_atoms top move 0.0 0.0 ${z_displace} units box

#在x和y方向固定原子

fix 1 all setforce 0 0 NULL

#能量最小化

min_style cg

minimize 1e-10 1e-10 1000 1000

#计算删除原子后的能量

variable Ef equal "c_eatoms"

#界面面积

variable A equal (${xdim}*${ydim})*1e-20

#单位转换系数

variable Cf equal 1.60217657e-16

#删除原子的能量

variable Er equal (${No}-${Nf})*${Ecoh}

#计算堆垛层错能

variable SFE equal ((${Ef}-(${Eo}-${Er}))*${Cf})/${A}

#输出模拟结果

print "The Stacking-fault energy = ${SFE} mJ/m^2"