lammps实例.pdf
- 文档编号:14657204
- 上传时间:2023-06-25
- 格式:PDF
- 页数:31
- 大小:5.30MB
lammps实例.pdf
《lammps实例.pdf》由会员分享,可在线阅读,更多相关《lammps实例.pdf(31页珍藏版)》请在冰点文库上搜索。
Project#5熔化与凝固:
氩,铜,铝熔化与凝固:
氩,铜,铝铜和铝的熔化转变:
铜和铝的熔化转变:
对于铜和铝,LAMMPS建立885的FCC晶格体系;充分弛豫后利用Nose-Hover方法,保持压强为零,使体系从T=2.5K开始加热,直至发生熔化转变。
下面是铜熔化的输入文件:
#LAMMPSMelt_CuorAlunitsmetal#单位,指定为lammps里的金属类的单位,长度为,能量为eV。
boundaryppp#周期性边界条件atom_styleatomic#原子模式variablexequal2.5#定义变量x为初始温度latticefcc3.61#Cu的晶格常数3.61#latticefcc4.05regionboxblock080805#x,y,z各方向上的晶胞重复单元数,也即区域大小create_box1box#将上述区域指定为模拟的盒子create_atoms1box#将原子按晶格填满盒子timestep0.01#步长0.005fsthermo1000#每隔1000步输出热力学结果pair_styleeam/alloy#选取Cu的EAM势作为模型pair_coeff*jin_copper_lammps.setflCu#EAM势文件名称#pair_styleeam/fs#pair_coeff*Al_FM.eam.fsAlneighbor0.5binneigh_modifyevery5delay0checkyes#velocityallcreate$x825577distgaussianfix1allnvt$x$x1.0drag0.2#保持初始温度,在NVT下弛豫#compute3allpe/atom#compute4allke/atom#compute5allcoord/atom3.0#dump1allcustom1dump.atomidxsyszsc_3c_4c_5run10000#运行10000步unfix1fix1allnpt$x20004.00xyz0.00.06.0drag0.2#在NPT下加热至2000K#fix1allnpt$x15004.00xyz0.00.06.0drag0.2run1200000#运行1200000步对于铜,其熔点为1357.77K,但在我们的模拟中其在1609K附近发生一级相变,比其平衡时熔点增大了18.7%。
对于铝,其熔点为934.477K,在我们的模拟中,其在1110.0K附近发生一级相变,单位体积发生突变。
其相对于平衡熔点增大了了18.9%。
从铜和铝的熔化过热,以及后面所涉及的氩的熔化过热与凝固过冷,我们可以看出,在利用分子动力学方法模拟熔化与凝固时往往会发生过热与过冷,其值基本在10%-30%之间。
产生过热与过冷的因素基本上可以从热力学与动力学的方面阐述,例如均匀形核而导致。
铜在特定温度下的性质铜在特定温度下的性质:
利用LAMMPS建立888的FCC格子,分别在10K,500K,1000K,1800K和2000K下保持零外压弛豫,得到在不同温度下原子运动的情况,以及不同情况下的均方根位移。
下面是LAMMPS的输入文件in.melt_Cu_temp#LAMMPSMelt_Cu_tempunitsmetalboundarypppatom_styleatomicvariablexindex10500100018002000print-Temperature=$xK-latticefcc3.62regionboxblock080808create_box1boxcreate_atoms1boxtimestep0.01thermo1000pair_styleeam/alloypair_coeff*jin_copper_lammps.setflCuneighbor0.6binneigh_modifyevery5delay0checkyesvelocityallcreate$x825577distgaussian#初始化速度,按高斯分布fix1allnpt$x$x2.0xyz0.00.06.0drag0.2compute3allpe/atomcompute4allke/atomcompute5allcoord/atom3.0run50000unfix1fix1allnvt$x$x2.0drag0.2dump1allcustom50dump_$x.atomidxsyszsc_3c_4c_5thermo100fix2allmsd1msd_Cu_$x.dat#输出msd文件run1000clearnextxjumpin.melt_Cu_temp原子在不同温度下的运动原子在不同温度下的运动T=10KT=500KT=1000KT=2000K均方根位移均方根位移:
模拟体系中的均方根位移可以通过如下公式求得:
是指相应量的统计平均值均方根位移的量与原子的扩散系数存在对应的关系。
当体系是固态时,即体系温度处于熔点之下时,均方根位移存在上限值;而当体系处于液态时,均方根位移呈线性关系,而且其斜率与原子的扩散系数存在如下关系:
在2维体系中上式的6应该被4所取代。
下图为T=2000K,即体系处于液态时的均方根位移图。
氩的熔化与凝固:
氩的熔化与凝固:
LAMMPS中对氩取LJ约化单位,其与国际单位制的转换如下:
mass=massormdistance=sigma,wherex*=x/sigmatime=tau,wheretau=t*=t(KbT/m/sigma2)1/2energy=epsilon,whereE*=E/epsilonvelocity=sigma/tau,wherev*=vtau/sigmaforce=epsilon/sigma,wheref*=fsigma/epsilontemperature=reducedLJtemperature,whereT*=TKb/epsilonpressure=reducedLJpressure,whereP*=Psigma3/epsilonforargon:
mass=6.6362126e-26(kg)sigma=3.405e-10(m);epsilon=1.6545e-21(J);Kb=8.314/(6.02e23);tau=2.156472211e-12(s)reducedLJtemperature=119.799(K)reducedLJvelocity=157.8967679(m/s)reducedLJpressure=41909784.02(Pa)氩的熔化转变氩的熔化转变:
对于氩,LAMMPS建立885的FCC晶格体系;充分弛豫后利用Nose-Hover方法,保持压强为零,使体系从T=0.01开始加热,直至发生熔化转变。
如下是输入文件in.melt_Ar#2dLennard-Jonesmeltunitsljatom_styleatomicboundaryppp#processors221latticefcc1.073regionboxblock080805create_box1boxcreate_atoms1boxmass11.0velocityallcreate0.01872877timestep0.01dump1allxyz1000melt.xyzpair_stylelj/cut2.5pair_coeff111.01.02.5neighbor0.3binneigh_modifyevery10delay0checkyesthermo1000fix1allnpt0.010.011.00xyz0.00.01.0drag0.2run50000unfix1thermo1000fix1allnpt0.010.852.0xyz0.00.01.0drag0.2run1000000在Tm=0.775左右发生一级相变,原子单位体积发生跃变;而氩的平衡熔化温度约为Tm=0.66,故其在分子动力学方法模拟下过热约20%。
氩的凝固转变:
氩的凝固转变:
对于氩,LAMMPS建立885的FCC晶格体系;充分弛豫后利用Nose-Hover方法,保持压强为零,使体系从T=0.85开始降温,发生凝固转变,直至温度降至0附近。
其下为输入文件in.quench#2dLennard-Jonesquenchunitsljatom_styleatomicboundaryppplatticefcc0.851regionboxblock080805create_box1boxcreate_atoms1boxmass11.0velocityallcreate0.85872877timestep0.01pair_stylelj/cut2.5pair_coeff111.01.02.5neighbor0.3binneigh_modifyevery10delay0checkyesthermo1000fix1allnpt0.850.852.0xyz0.00.01.00run50000unfix1dump1allxyz1000quench.xyzfix1allnpt0.850.012.0xyz0.00.01.00run1000000如下图所示,在分子动力学模拟下,氩在T=0.44附近发生一级相变,单位原子体积突然下降:
相对与其平衡凝固温度Tm=0.44,约有30%左右的过冷。
特定温度下的性质特定温度下的性质:
利用LAMMPS建立888的FCC格子,分别在T=0.1,0.4,0.6,0.8下保持零外压弛豫,得到在不同温度下原子运动的情况,以及不同径向分布函数和速度自相关函数。
下面是LAMMPS的输入文件in.melt_Ar_temp#2dLennard-Jonesmelt_tempunitsljatom_styleatomicboundarypppvariablexindex0.800.600.400.10latticefcc0.888regionboxblock080808create_box1boxcreate_atoms1boxmass11.0print-Temperature=$x-velocityallcreate$x872877timestep0.01pair_stylelj/cut2.5pair_coeff111.01.02.5neighbor0.3binneigh_modifyevery20delay0checknothermo1000fix1allnpt$x$x2.0xyz0.00.01.0run50000unfix1fix1allnvt$x$x1.0compute3allpe/atomcompute4allke/atomcompute5allcoord/atom3.0dump1allcustom50dump_$x.atomidxsyszsc_3c_4c_5run1000clearnextxjumpin.melt_temp径向分布函数:
径向分布函数:
径向分布函数即是原子径向上的原子密度与体系总密度的比值函数(具体参见Project_5.1)。
下图为四个温度下的氩体系的径向分布函数。
可以发现在温度较低的情况下,径向分布函数的峰比较尖锐,各个峰所对应的径向值,分别对应最近邻、次近邻等配位的位置;随着温度的升高,径向分布函数的峰变宽,一些位置上的峰消失;直至到达液相,此时分布函数的值不再表示配位情况,而是反映了此时其它原子相对于中心原子的位置的概率分布。
速度自相关函数:
速度自相关函数:
速度自相关函数的物理意义为原子某一时候的速度与其在初始时速度的函数关系,如下(具体参见Project_5.1)Cv(t)=.下图为4个不同温度下的速度自相关函数。
可以发现如下规律:
随着时间的变大,原子的速度基本上与初始值没有相关关系(即Cv(t)0)。
另外,随着温度的提高,原子的振动幅度与频率上升,自相关函数的第一个波谷深度变小,而且到达无关速度所花费的时间减少。
1Project#1硅的晶格常数和体弹模量的计算硅的晶格常数和体弹模量的计算一、平衡晶格常数和内聚能一、平衡晶格常数和内聚能自然条件下硅为金刚石结构(dc)。
计算模拟时,我们可以假定它为各种结构,fcc,bcc,sc,dc.可以预测,模拟的dc结构的硅的体系能量最低,也即最稳定。
下面我们将运用LAMMPS来对硅的各种结构进行模拟。
定义晶格能量为,数密度为:
potENNV其中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$grepdc.log$lmp_serialin.Siliconlogdc.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#划定区域,x0,3,y0,3z0,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.43153=0.049926,也即我们直接通过grep命令得到的第二项值;纵坐标为每个原子的能量,为第三项值。
金刚石之外,还需计算其他结构。
只需对in.Silicon做稍微改动:
首先,将in.Silicon复制成in.fcc:
然后编辑in.fcc改动如下几项:
然后如下命令执行:
相应的,如下命令查看log文件中的数据:
$grepfcc.log$lmp_serialin.fcclogfcc.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().52345E104799.5791127112.82866a15604.99601a7592.99351a1134.25198a57.82091a二、体弹模量二、体弹模量我们同时可以从晶格能曲线在最低处得到体弹模量的信息。
体弹模量定义为:
/dPBdVVV和P分别为晶胞的体积和体系的压强。
我们已经得到了内聚能与晶格参数的函数关系,对于立方晶胞而言,23dMdEPdVada由此,02209aMdEBada其中a0为平衡晶格常数,M为体积为3Va的晶胞中的原子数目。
从多项式拟合,可以得到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$lammps2cfg1Project#2金属中的点缺陷:
空位和间隙原子金属中的点缺陷:
空位和间隙原子一、空位一、空位从晶体中移去一个原子,即可形成空位。
本例将运用LAMMPS计算空位形成能,Ev.LAMMPS输入文件为in.vacancy1)在在fcc结构的完整结构的完整Cu晶体中引入一个空位晶体中引入一个空位沿方向构造一个4NNN的晶体。
N为input文件中lattice命令指定的个方向上的晶胞重复单元数。
2)弛豫弛豫当一个原子从晶体中移走之后,周围的原子将相应地调整位置以降低体系势能。
为得到稳定的构型,需要对体系进行弛豫,relaxation.LAMMPS提供两种能量最小化方式,cg和sd。
本例中选用sd方式进行能量最小化。
如下是输入文件,in.vacancy:
unitsmetal#单位为lammps中的metel类型boundaryppp#周期性边界条件atom_styleat
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- lammps 实例