lammps实例(5)

Project #5
熔化与凝固：氩，铜，铝
铜和铝的熔化转变:
对于铜和铝，LAMMPS建立8×8×5的FCC晶格体系;充分弛豫后利用Nose-Hover方法，保持压强为零，使体系从T=2.5K开始加热，直至发生熔化转变。

下面是铜熔化的输入文件：
# LAMMPS Melt_Cu or Al
units metal # 单位，指定为lammps里的金属类的单位，长度为Å，能量为eV。

boundary p p p # 周期性边界条件
atom_style atomic # 原子模式
variable x equal 2.5 # 定义变量 x为初始温度
lattice fcc 3.61 # Cu 的晶格常数3.61
#lattice fcc 4.05
region box block 0 8 0 8 0 5 # x,y,z各方向上的晶胞重复单元数，也即区域大小
create_box 1 box # 将上述区域指定为模拟的盒子
create_atoms 1 box # 将原子按晶格填满盒子
timestep 0.01 # 步长 0.005fs
thermo 1000 #每隔1000步输出热力学结果
pair_style eam/alloy # 选取 Cu 的EAM势作为模型
pair_coeff * * jin_copper_lammps.setfl Cu # EAM 势文件名称
#pair_style eam/fs
#pair_coeff * * Al_FM.eam.fs Al
neighbor 0.5 bin
neigh_modify every 5 delay 0 check yes
#velocity all create $x 825577 dist gaussian
fix 1 all nvt $x $x 1.0 drag 0.2 #保持初始温度，在NVT下弛豫
#compute 3 all pe/atom
#compute 4 all ke/atom
#compute 5 all coord/atom 3.0
#dump 1 all custom 1 dump.atom id xs ys zs c_3 c_4 c_5
run 10000 # 运行10000步
unfix 1
fix 1 all npt $x 2000 4.00 xyz 0.0 0.0 6.0 drag 0.2 #在NPT下加热至2000K
#fix 1 all npt $x 1500 4.00 xyz 0.0 0.0 6.0 drag 0.2
run 1200000 #运行1200000步
对于铜，其熔点为1357.77K，但在我们的模拟中其在1609K附近发生一级相变，比其平衡时熔点增大了18.7%。

对于铝，其熔点为934.477K, 在我们的模拟中，其在1110.0K附近发生一级相变，单位体积发生突变。

其相对于平衡熔点增大了了18.9%。

从铜和铝的熔化过热，以及后面所涉及的氩的熔化过热与凝固过冷，我们可以看出，在利用分子动力学方法模拟熔化与凝固时往往会发生过热与过冷，其值基本在10%-30%之间。

产生过热与过冷的因素基本上可以从热力学与动力学的方面阐述，例如均匀形核而导致。

铜在特定温度下的性质:
利用LAMMPS建立8×8×8的FCC格子，分别在10K，500K，1000K， 1800K和2000K下保持零外压弛豫，得到在不同温度下原子运动的情况，以及不同情况下的均方根位移。

下面是LAMMPS的输入文件in.melt_Cu_temp
# LAMMPS Melt_Cu_temp
units metal
boundary p p p
atom_style atomic
variable x index 10 500 1000 1800 2000
print "-------------------------------Temperature=$x K---------"
lattice fcc 3.62
region box block 0 8 0 8 0 8
create_box 1 box
create_atoms 1 box
timestep 0.01
thermo 1000
pair_style eam/alloy
pair_coeff * * jin_copper_lammps.setfl Cu
neighbor 0.6 bin
neigh_modify every 5 delay 0 check yes
velocity all create $x 825577 dist gaussian #初始化速度，按高斯分布
fix 1 all npt $x $x 2.0 xyz 0.0 0.0 6.0 drag 0.2
compute 3 all pe/atom
compute 4 all ke/atom
compute 5 all coord/atom 3.0
run 50000
unfix 1
fix 1 all nvt $x $x 2.0 drag 0.2
dump 1 all custom 50 dump_$x.atom id xs ys zs c_3 c_4 c_5
thermo 100
fix 2 all msd 1 msd_Cu_$x.dat #输出msd文件run 1000
clear
next x
jump in.melt_Cu_temp
原子在不同温度下的运动
T=10K T=500K
T=1000K T=2000K
均方根位移:
模拟体系中的均方根位移可以通过如下公式求得：
是指相应量的统计平均值
均方根位移的量与原子的扩散系数存在对应的关系。

当体系是固态时，即体系温度处于熔点之下时，均方根位移存在上限值；而当体系处于液态时，均方根位移呈线性关系，而且其斜率与原子的扩散系数存在如下关系：
在2维体系中上式的6应该被4所取代。

下图为T=2000K，即体系处于液态时的均方根位移图。

氩的熔化与凝固：
LAMMPS中对氩取LJ约化单位，其与国际单位制的转换如下:
• mass = mass or m
• distance = sigma, where x* = x / sigma
• time = tau, where tau = t* = t (Kb T / m / sigma^2)^1/2
• energy = epsilon, where E* = E / epsilon
• velocity = sigma/tau, where v* = v tau / sigma
• force = epsilon/sigma, where f* = f sigma / epsilon
• temperature = reduced LJ temperature, where T* = T Kb / epsilon
• pressure = reduced LJ pressure, where P* = P sigma^3 / epsilon
for argon: 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)
reduced LJ temperature=119.799(K)
reduced LJ velocity=157.8967679(m/s)
reduced LJ pressure=41909784.02(Pa)
氩的熔化转变:
对于氩，LAMMPS建立8×8×5的FCC晶格体系;充分弛豫后利用Nose-Hover方法，保持压强为零，使体系从T=0.01开始加热，直至发生熔化转变。

如下是输入文件 in.melt_Ar
# 2d Lennard-Jones melt
units lj
atom_style atomic
boundary p p p
#processors 2 2 1
1.073
lattice fcc
region box block 0 8 0 8 0 5
box
create_box 1
create_atoms 1 box
1.0
mass 1
velocity all create 0.01 872877
timestep 0.01
dump 1 all xyz 1000 melt.xyz
pair_style l j/cut 2.5
pair_coeff 1 1 1.0 1.0 2.5
neighbor 0.3 bin
neigh_modify every 10 delay 0 check yes
thermo 1000
fix 1 all npt 0.01 0.01 1.00 xyz 0.0 0.0 1.0 drag 0.2
run 50000
unfix 1
thermo 1000
fix 1 all npt 0.01 0.85 2.0 xyz 0.0 0.0 1.0 drag 0.2
run 1000000
在Tm=0.775左右发生一级相变，原子单位体积发生跃变；而氩的平衡熔化温度约为Tm=0.66，故其在分子动力学方法模拟下过热约20%。

氩的凝固转变：
对于氩，LAMMPS建立8×8×5的FCC晶格体系;充分弛豫后利用Nose-Hover方法，保持压强为零，使体系从T=0.85开始降温，发生凝固转变，直至温度降至0附近。

其下为输入文件in.quench
# 2d Lennard-Jones quench
units lj
atom_style atomic
boundary p p p
0.851
lattice fcc
region box block 0 8 0 8 0 5
box
create_box 1
create_atoms 1 box
1.0
mass 1
velocity all create 0.85 872877
timestep 0.01
pair_style l j/cut 2.5
pair_coeff 1 1 1.0 1.0 2.5
neighbor 0.3 bin
neigh_modify every 10 delay 0 check yes
thermo 1000
fix 1 all npt 0.85 0.85 2.0 xyz 0.0 0.0 1.00
run 50000
unfix 1
dump 1 all xyz 1000 quench.xyz
fix 1 all npt 0.85 0.01 2.0 xyz 0.0 0.0 1.00
run 1000000
如下图所示，在分子动力学模拟下，氩在T=0.44附近发生一级相变，单位原子体积突然下降：相对与其平衡凝固温度Tm=0.44，约有30%左右的过冷。

特定温度下的性质:
利用LAMMPS建立8×8×8的FCC格子，分别在T=0.1，0.4，0.6，0.8下保持零外压弛豫，得到在不同温度下原子运动的情况, 以及不同径向分布函数和速度自相关函数。

下面是LAMMPS的输入文件in.melt_Ar_temp
# 2d Lennard-Jones melt_temp
units lj
atom_style atomic
boundary p p p
variable x index 0.80 0.60 0.40 0.10
0.888
lattice fcc
region box block 0 8 0 8 0 8
create_box 1
box
create_atoms 1 box
1.0
mass 1
print "----------------------------Temperature=$x--------------------"
velocity all create $x 872877
timestep 0.01
pair_style l j/cut 2.5
pair_coeff 1 1 1.0 1.0 2.5
neighbor 0.3 bin
neigh_modify every 20 delay 0 check no
thermo 1000
fix 1 all npt $x $x 2.0 xyz 0.0 0.0 1.0
run 50000
unfix 1
fix 1 all nvt $x $x 1.0
compute 3 all pe/atom
compute 4 all ke/atom
compute 5 all coord/atom 3.0
dump 1 all custom 50 dump_$x.atom id xs ys zs c_3 c_4 c_5 run 1000
clear
next x
jump in.melt_temp
径向分布函数即是原子径向上的原子密度与体系总密度的比值函数（具体参见Project_5.1）。

下图为四个温度下的氩体系的径向分布函数。

可以发现在温度较低的情况下，径向分布函数的峰比较尖锐，各个峰所对应的径向值，分别对应最近邻、次近邻等配位的位置；随着温度的升高，径向分布函数的峰变宽，一些位置上的峰消失；直至到达液相，此时分布函数的值不再表示配位情况，而是反映了此时其它原子相对于中心原子的位置的概率分布。

速度自相关函数的物理意义为原子某一时候的速度与其在初始时速度的函数关系，如下（具体参见Project_5.1）
C v(t) = < v i(0) . v i(t) >.
下图为4个不同温度下的速度自相关函数。

可以发现如下规律：随着时间的变大，原子的速度基本上与初始值没有相关关系（即Cv(t)‐>0）。

另外，随着温度的提高，原子的振动幅度与频率上升，自相关函数的第一个波谷深度变小，而且到达无关速度所花费的时间减少。