lammps实例(4)

1  Project#2金属中的点缺陷：空位和间隙原子 一、空位从晶体中移去一个原子，即可形成空位。本例将运用LAMMPS计算空位形成能,Ev.LAMMPS输入文件为in.vacancy1)在fcc结构的完整Cu晶体中引入一个空位沿100方向构造一个4×N×N×N的晶体。N为input文件中lattice命令指定的个方向上的晶胞重复单元数。2)弛豫当一个原子从晶体中移走之后，周围的原子将相应地调整位置以降低体系势能。为得到稳定的构型，需要对体系进行弛豫，relaxation.LAMMPS提供两种能量最小化方式，cg和sd。本例中选用sd方式进行能量最小化。如下是输入文件，in.vacancy:unitsmetal#单位为lammps中的metel类型boundaryppp#周期性边界条件atom_styleatomic#原子模式latticefcc3.61#Cu的晶格常数3.61regionboxblock060606#x,y,z各方向上的晶胞重复单元数，也即区域大小create_box1box#将上述区域指定为模拟的盒子create_atoms1box#将原子按晶格填满盒子pair_styleeam/alloy#选取Cu的EAM势作为模型pair_coeff**jin_copper_lammps.setflCu#EAM势文件名称run0#运行0步，仅为启动lammps的热力学数据计算variableEequalpe#定义变量E为系统总势能variableNequalatoms#定义变量N为系统总原子数printthenumberofatoms&systemenergynoware$N$E#打印信息regioncenterpointblock33.0533.0533.05#指定一个原子大小的区域delete_atomsregioncenterpoint#删除这个区域的原子min_stylesd#能量最小化模式，sdminimize1.0e-121.0e-1210001000#能量最小化参数，指数越大最小化程度越深printthenumberofatoms&systemenergynoware$N$E#打印信息compute3allpe/atom#计算每个原子的势能compute4allke/atom#计算每个原子的动能compute5allcoord/atom3.0#计算每个原子的近邻原子数dump1allcustom1dump.atomidxsyszsc_3c_4c_5#将指定的各原子信息写入dump.atom.timestep0.005#步长0.005fsrun1#运行1步2  3)运行lammps4)计算空位形成能空位浓度由下式给出：[n]=exp(−Fv/kBT).其中Fv=Ev−TSv为形成一个空位所需要的Helmholtz自由能.忽略熵Sv,空位浓度公式简化为[n]=exp(−Ev/kBT).设E1为完整晶体能量，含N个原子；E2为弛豫后的晶体能量，含N–1个原子。空位形成能Ev为：211vNEEEN或21vcohEENE,其中Ecoh=E1/N,为完整晶体的内聚能。本例中以EAM模型计算4×(20×20×20)=32000个原子的体系，得到空位形成能Ev~1.26eV，文献中的实验值为~1.28eV，符合较好。另由上式计算得到，300K温度下的空位浓度为~7.59×10-22，1350K(Tm)时的空位浓度~2.2×10-5（文献中的实验值为~2×10-4）。换算时注意(1eV/kB=1.1604×10−4K)图1.空位处于4×(6×6×6)的FCC晶体中心，106ca,206ca,306ca.颜色依据原子势能标注。3  二、间隙原子向完整晶体中插入一个原子，即形成间隙原子。如果新插入的原子和晶体原子相同，则为自间隙原子，self-interstitial。与空位计算类似，我们用如下式子计算金属Cu中的自间隙原子形成能Ei211iNEEENEi可能取决于间隙原子引入的初始位置。但对应最小的Ei值的，应该是唯一的位置。在Cu,Ni和Pt等FCC金属中，最稳定的自间隙原子构型均为[100]方向的哑铃型。如下图：图2.FCC金属中的哑铃状自间隙原子构型。右图为LAMMPS得到的Cu的构型。晶体中引入间隙原子后，周围原子将做相应的位置调整以期达到最低能量状态。为了得到弛豫后的构型，我们采用LAMMPS里的cg和sd的能量最小化方法。本例中采用的是sd方法。相对于空位，间隙原子的引入需要更大程度的弛豫。结合能量最小化方法，我们采用NVT或NVE系统的热力学平衡方法。给体系升温，让原子充分动起来，找到最稳定的位置，得到最稳定的构型。然后淬火quench到0K。最后再运用能量最小化。Theinputfilein.interstitiallookslike4  计算32000个原子的体系，得到Cu的自间隙原子形成能为Ei~3.1eV.类似上述计算，Cu在T=300K和1350K(Tm)时的间隙原子浓度分别为~8.4×10-53和~2.7×10-12.讨论如果体系未得到充分弛豫，可以得到各种不同的间隙原子构型，如图3和图4.1.观察Ei与模拟体系大小的关系。改变盒子大小。2.改变间隙原子的引入位置，计算可能的间隙原子构型，并指出最稳定的间隙原子构型和形成能。unitsmetal#单位为lammps中的metel类型boundaryppp#周期性边界条件atom_styleatomic#原子模式latticefcc3.61#Cu的晶格常数3.61regionboxblock040404#x,y,z各方向上的晶胞重复单元数，也即区域大小create_box1box#将上述区域指定为模拟的盒子create_atoms1box#将原子按晶格填满盒子pair_styleeam/alloy#选取Cu的EAM势作为模型pair_coeff**jin_copper_lammps.setflCu#EAM势文件名称run0#运行0步，仅为启动lammps的热力学数据计算variableEequalpe#定义变量E为系统总势能variableNequalatoms#定义变量N为系统总原子数printthenumberofatoms&systemenergynoware$N$E#打印信息create_atoms1single2.452.052.05#在该位置插入一个原子min_stylesd#能量最小化模式，sdminimize1.0e-121.0e-1210001000#能量最小化参数，指数越大最小化程度越深printinterstitialintroduced,minimized:$Natoms,energyis$Efix1allnvt100100100drag0.2#nvt系综，原子数、体积和温度保持不变;T=100Ktimestep0.005#步长0.005fsrun1000#运行1000步printnvtperformed,temperatureup:$Natoms,totalenergyis$Efix1allnvt1000.0001100drag0.2#nvt系综，温度由100K到0.0001Krun1000#运行1000步printnvtperformed,temperaturedown:$Natoms,totalenergyis$Ecompute3allpe/atom#计算每个原子的势能compute4allke/atom#计算每个原子的动能compute5allcoord/atom3.0#计算每个原子的近邻原子数dump1allcustom1dump.atomidxsyszsc_3c_4c_5#将信息写入dump.atommin_stylesdminimize1.0e-121.0e-121000010000#再次能量最小化printthefinalstate:$Natoms,totalenergyis$E#打印信息5  (a)(b)(c)图3.(a)FCC晶体中的八面体（红色）和四面体（蓝色）间隙位置；(b)LAMMPS计算所得最稳定的间隙原子位置，恰为八面体中心；(c)LAMMPS计算所得的另一构型，为四面体中心，体系能量比较高。6  图4.LAMMPS计算得到的五种可能的间隙原子构型。