lammps实例(3)

1  Project#1硅的晶格常数和体弹模量的计算 一、平衡晶格常数和内聚能自然条件下硅为金刚石结构(dc)。计算模拟时，我们可以假定它为各种结构，fcc,bcc,sc,dc.可以预测，模拟的dc结构的硅的体系能量最低，也即最稳定。下面我们将运用LAMMPS来对硅的各种结构进行模拟。定义晶格能量为Φ，数密度为ρ：potENNV其中Epot为势能，N为体系总原子数，V为体系的体积。选取Stillinger-Weber(SW)，以下面命令执行lammps运算：其中，lmp_serial为lammps命令；””符号为读取符；in.Silicon为输入文件，里面包含运算所需要的各种数据和命令；-log指定输出文件的名称。可以看到屏幕上显示出lammps运行的信息。这个计算量很小，所以很快就结束。接下来以如下命令来查看计算得到的数据：grep是linux中一个很重要的命令，用来搜索文本，读取匹配的行并打印出来。这里是搜索dc.log文件，将@开头的行打印出来。如下：晶格参数为5.4305埃,数密度为0.0499540303，每个原子的能量为-4.336599609eV.@@@@(latticeparameter,rho,energyperatom):5.43050.0499540303-4.336599609@@@@(latticeparameter,rho,energyperatom):5.43060.04995127076-4.336599763@@@@(latticeparameter,rho,energyperatom):5.43070.04994851143-4.336599879$grep^@dc.log$lmp_serialin.Silicon–logdc.log2  下面具体来看刚才给的输入文件，in.Silicon.dc.log文件中有原子总数的信息，每个金刚石晶胞中有8个原子，383216，所以是216个原子。如下给出各种结构下的体系的原子数：^^^^^^^^^^^^^^^^^Thisisanewloop^^^^^^^^^^^^^^#注释行，标志新循环开始Latticespacinginx,y,z=5.43155.43155.4315#晶胞大小，也即为指定的晶格参数xCreatedorthogonalbox=(000)to(16.294516.294516.2945)#构造的正交盒子的坐标1by1by1processorgrid#一个核在运算，串行Created216atoms#构造了216个原子#bulkSiliconlattice#注释行，随便给unitsmetal#单位，指定为lammps里的金属类的单位，长度为Å，能量为eV。boundaryppp#周期性边界条件atom_styleatomic#原子模式variablexindex5.43055.43065.43075.43085.43095.43105.43115.43125.43135.43145.4315#定义变量x，在运行中x逐一取这些值。本例中为各个晶格常数。latticediamond$x#晶格，指定金刚石结构的晶格，晶格常数为x的值#latticediamond5.431#latticefcc3.615#如果计算fcc结构的晶格，则将晶格常数取在3.615附近#latticebcc3.28#同上#latticesc2.60#同上regionboxblock030303#划定区域，x∈[0,3],y∈[0,3]z∈[0,3]，单位为晶胞create_box1box#在上面这个区域里创建一个模拟的盒子create_atoms1box#将这个盒子按晶格填满一种原子pair_stylesw#选取sw势pair_coeff**Si.swSi#势文件名为Si.swmass128#给定硅的质量，此处与势对应neighbor1.0binneigh_modifyevery1delay5checkyesvariablePequalpe/216#定义P为每个原子的势能。216为原子数，pe为体系总势能variablerequal216/($x*3)^3#定义r为数密度，单位体积里的原子数timestep0.005#步长为0.005飞秒thermo10#每10步在屏幕上打印一次热力学状态信息min_stylesd#能量最小化模式，lammps提供sd和cg可选minimize1.0e-121.0e-1210001000#能量最小化参数设置，第一项和第二项为能量和力的判据，指数越高，最小化程度越高，效果越好，但计算量也多。后两项为限制的最多迭代次数。compute3allpe/atom#计算每个原子的势能compute4allke/atom#计算每个原子的动能compute5allcoord/atom3.0#计算每个原子的近邻原子数dump1allcustom1dump.atomidxsyszsc_3c_4c_5#输出到dump文件中print@@@@(latticeparameter,rho,energyperatom):$x$r$P#打印clear#清除该次循环里的数据信息nextx#跳转到下一个xjumpin.Silicon#跳到in.Silicon文件，再从头读起。也即循环。3  晶体结构类型晶胞中的原子数总原子数简单立方SC127体心立方BCC254面心立方FCC4108金刚石DC8216表1.不同晶体结构中的原子数下图是计算模拟得出的各种结构下的数密度与每个原子能量的关系图。横坐标为数密度，以金刚石为例，ρ=8/5.4315^3=0.049926，也即我们直接通过grep命令得到的第二项值；纵坐标为每个原子的能量，为第三项值。金刚石之外，还需计算其他结构。只需对in.Silicon做稍微改动：首先，将in.Silicon复制成in.fcc:然后编辑in.fcc改动如下几项：然后如下命令执行：相应的，如下命令查看log文件中的数据：$grep^@fcc.log$lmp_serialin.fcc–logfcc.logvariablexindex3.6103.6113.6123.6133.6143.6153.6163.6173.618#设置x值，为平衡晶格常数3.615附近的值。#latticediamond$x#将金刚石这项注释掉#latticediamond5.431latticefcc3.615#启用此项，fcc结构#latticebcc3.28#latticesc2.60variablePequalpe/108#此时原子总数为108variablerequal108/($x*3)^3#相应修改为108$geditin.fcc$cpin.Siliconin.fcc4  以同样方法编辑in.bcc,in.sc，计算不同晶格参数时的体系能量值，并绘制下图：图1.不同结构下的硅的晶格能。可以看出金刚石结构对应最低能量，最为稳定下图更为细致地画出金刚石结构中，不同晶格参数所对应的内聚能。内聚能(cohesiveenergyEcoh)的定义是，最小的晶格能。由图可以得到，平衡晶格常数为a0=5.431Å，内聚能为Ecoh=−4.3365eV.图2.金刚石结构中的晶格能VS晶格参数。五阶拟合得到平衡晶格常数5.43095(Å).5  2345E104799.5791127112.82866a15604.99601a7592.99351a1134.25198a57.82091a二、体弹模量我们同时可以从晶格能曲线在最低处得到体弹模量的信息。体弹模量定义为：/dPBdVVV和P分别为晶胞的体积和体系的压强。我们已经得到了内聚能与晶格参数的函数关系，对于立方晶胞而言，23dMdEPdVada由此，02209aMdEBada其中a0为平衡晶格常数，M为体积为3Va的晶胞中的原子数目。从多项式拟合，可以得到022adEda,此例中为3.87081eV/Å2(a0=5.431Å).由上面公式，计算得到硅的体弹模量为B=101.366GPa.文献中的实验数据为99GPa.可视化每一次lammps运行后，会生成一个dump.atom文件。可以通过如下命令转换为Atomeye可读取的.cfg文件：这时当前目录中就会生成001.cfg之类的.cfg文件，然后通过Atomeye查看：Tab键可以改变观察方位上下左右箭头键可以转动原子PgUp和PgDn改变原子大小Alt+1和Alt+2以及Alt+3改变原子颜色先按一次9这个键，转动时就是每次以90度转动$A.i686001.cfg$lammps2cfg