|
|
本帖最后由 xwjjie 于 2024-6-27 11:55 编辑
各位老师，我最近在做，使用dp势函数在lammps中模拟计算晶体的热导率。但是我无论是用平衡法[size=9.9626pt]Green-Kubo法，还是非平[size=9.9626pt]衡法langevin控温法，得到的结果都和预期不相符。我查看了我的in文件，似乎并没有问题，希望各位老师能够提点一下我，给出点意见，谢谢大家。
[size=13.2835px]
[size=13.2835px]Green-Kubo法in文件参数：
# 模拟参数
units metal
dimension 3
boundary p p p
atom_style atomic
neighbor 0.3 bin
neigh_modify delay 0 every 1
# 读取模型体系模型
read_data element.data
# 创建x*x*x的box
replicate 5 5 5
# 力场设置
pair_style deepmd element-compress.pb
pair_coeff * *
# 初始化参数
variable T equal 300 # temperature
#热力学信息输出
#----- initialization -----
velocity all create $T 98989 dist gaussian
#-----relax process-----
#温度初始化
fix 1 all npt temp $T $T 0.1 iso 0 0 1 drag 0.5
thermo 10000
thermo_style custom step temp vol press
dump 1 all custom 10000 element-npt.xyz id type x y z
run 100000
unfix 1
undump 1
# 1st equilibration run
fix 1 all nvt temp $T $T 0.5
# 热力学信息输出
thermo 10000
thermo_style custom step temp vol press
dump 1 all custom 10000 element-nvt.xyz id type x y z
run 100000
unfix 1
undump 1
reset_timestep 0
variable s equal 20 # sample interval
variable p equal 3000 # correlation length
variable d equal $p*$s # dump interval
variable w equal 200
variable r equal $w*$d
#单位转换公式
variable kB equal 1.3806504e-23 # [J/K] Boltzmann
variable eV2J equal 1.60218E-19 # eV to Joule
variable A2m equal 1.0e-10 # Angstrom to metre
variable ps2s equal 1.0e-12 # picoseconds to s
variable convert equal ${eV2J}*${eV2J}/${ps2s}/${A2m}
variable dt equal 0.001
timestep ${dt}
# thermal conductivity calculation
compute myKE all ke/atom # 计算原子动能
compute myPE all pe/atom #计算原子势能
compute myStress all centroid/stress/atom NULL virial # 计算原子应力
compute flux all heat/flux myKE myPE myStress
variable Jx equal c_flux[1]/vol
variable Jy equal c_flux[2]/vol
variable Jz equal c_flux[3]/vol
fix 1 all nve
fix JJ all ave/correlate $s $p $d &
c_flux[1] c_flux[2] c_flux[3] type auto &
file heatflux.txt ave running
#
variable Temp equal temp
variable Vol equal vol
variable scale equal ${convert}*$s*${dt}/${Temp}/${Temp}/${Vol}/${kB}
variable k11 equal trap(f_JJ[3])*${scale}
variable k22 equal trap(f_JJ[4])*${scale}
variable k33 equal trap(f_JJ[5])*${scale}
variable kappa equal (v_k11+v_k22+v_k33)/3.0
variable custom_step equal step
variable custom_temp equal temp
variable kappa_value equal v_kappa
# 启动运算
thermo $d
thermo_style custom step temp v_Jx v_Jy v_Jz v_k11 v_k22 v_k33 v_kappa
fix myprint all print $d "${custom_step} ${custom_temp} ${kappa_value} " file kappa_heatflux.txt screen no
run $r
variable N equal count(all)
variable NUM equal $N
print "Totel number is ${NUM}"
print "Running average thermal conductivity: $(v_kappa)"
非平[size=9.9626pt]衡法langevin控温法in文件参数：
# 模拟参数
units metal
dimension 3
boundary p p p
atom_style atomic
neighbor 0.3 bin
neigh_modify delay 0 every 1
# 读取模型体系模型
read_data element.data
# 创建x*x*x的box
replicate 4 4 8
# 力场设置
pair_style deepmd element-compress.pb
pair_coeff * *
variable T equal 300
#热力学信息输出
# velocity all create $T 98989 dist gaussian
thermo 1000
thermo_style custom step temp vol press
#温度初始化
fix 1 all npt temp ${T} ${T} 0.1 iso 0 0 1 drag 0.5
dump 1 all custom 1000 element-npt.xyz id type x y z
run 100000
unfix 1
undump 1
reset_timestep 0
# variable dt equal 0.0001
# timestep ${dt}
variable thi equal ${T}+50
variable tlo equal ${T}-50
#单位转换公式
variable kB equal 1.3806504e-23 # [J/K] Boltzmann
variable eV2J equal 1.60218E-19 # eV to Joule
variable A2m equal 1.0e-10 # Angstrom to metre
variable ps2s equal 1.0e-12 # picoseconds to s
variable convert equal ${eV2J}/${ps2s}/${A2m}
# heat layers
region hot block INF INF INF INF 0 1
region cold block INF INF INF INF 10 11
compute Thot all temp/region hot
compute Tcold all temp/region cold
# 1st equilibration run
fix 1 all nvt temp ${T} ${T} 0.5
thermo 1000
thermo_style custom step temp vol press
dump 1 all custom 1000 element-nvt.xyz id type x y z
run 100000
unfix 1
undump 1
reset_timestep 0
# 2nd equilibration run
thermo 1000
fix 1 all nve
fix hot all langevin ${thi} ${thi} 1.0 59804 tally yes
fix cold all langevin ${tlo} ${tlo} 1.0 287859 tally yes
fix_modify hot temp Thot
fix_modify cold temp Tcold
variable tdiff equal c_Thot-c_Tcold
thermo_style custom step temp c_Thot c_Tcold f_hot f_cold v_tdiff
thermo_modify colname c_Thot Temp_hot colname c_Tcold Temp_cold &
colname f_hot E_hot colname f_cold E_cold &
colname v_tdiff dTemp_step
run 100000
reset_timestep 0
# thermal conductivity calculation
# reset langevin thermostats to zero energy accumulation
thermo 1000
compute ke all ke/atom
variable temp atom c_ke/${kB}/1.5*${eV2J}
fix hot all langevin ${thi} ${thi} 1.0 59804 tally yes
fix cold all langevin ${tlo} ${tlo} 1.0 287859 tally yes
fix_modify hot temp Thot
fix_modify cold temp Tcold
fix ave all ave/time 10 100 1000 v_tdiff ave running
thermo_style custom step temp c_Thot c_Tcold f_hot f_cold v_tdiff f_ave
thermo_modify colname c_Thot Temp_hot colname c_Tcold Temp_cold &
colname f_hot E_hot colname f_cold E_cold &
colname v_tdiff dTemp_step colname f_ave dTemp
compute layers all chunk/atom bin/1d z lower 0.05 units reduced
fix 2 all ave/chunk 10 100 1000 layers v_temp file langevin.txt
variable start_time equal time
variable kappa equal ${convert}*(0.5*(abs(f_hot)+abs(f_cold))/(time-${start_time})/(lx*ly)/2)*(lz/2.0)/f_ave
variable custom_step equal step
variable Energy_hot equal abs(f_hot)
variable Energy_cold equal abs(f_cold)
variable dtemp equal f_ave
variable axis_X equal lx
variable axis_Y equal ly
variable axis_Z equal lz
variable dtime equal time-${start_time}
fix myprint all print 1000 &
"${custom_step} ${convert} ${axis_X} ${axis_Y} ${axis_Z} ${dtime} ${Energy_hot} ${Energy_cold} ${dtemp}" &
file kappa_langevin.txt screen no
run 2000000
|
|
发表于 Post on 2024-6-26 11:32:07
收藏 Add to favorites