我在用NEMD计算Cu/diamond界面热导时，遇到了用Langevin热浴控温的情况下冷热源热流混乱的问题，具体情况是统计结果中热源的f_hot值时正时负其绝对值的大小在整体上也不是增大的，f_cold也出现了类似的情况，不知具体原因出在哪里，望解惑。
#Cu-diamond界面热阻
#---------initialization---------
units metal
dimension 3
boundary p p s
atom_style atomic
#构建模型
lattice fcc 3.61
region box block 0 36 0 36 0 140 units box
create_box 2 box
create_atoms 1 box
region cu_del block INF INF INF INF 70 140 units box
delete_atoms region cu_del
lattice diamond 3.57
region dia block INF INF INF INF 70 140 units box
create_atoms 2 region dia
#设定原子质量
mass 1 63.5
mass 2 12.011
#--------force-field------------
neighbor 0.3 bin
neigh_modify delay 0
pair_style hybrid/overlay morse 10.0 eam tersoff
pair_coeff 1 1 eam Cu_u3.eam
pair_coeff * * tersoff c.tersoff NULL C
pair_coeff 1 2 morse 0.087 51.4 0.205
##### setting some important variables
timestep 0.0001
variable Time equal step*dt/1000 #time in picoseconds
#### 定义热源和热汇
velocity all create 300 123456
# ----------------- Minimization-----------------
min_style cg
min_modify dmax 0.4
minimize 1e-15 1e-15 5000 5000
reset_timestep 0
#------------------NVT -equilibrium-run-----------
fix 1 all nvt temp 300 300 0.01
thermo_style custom step v_Time temp press vol etotal ke pe evdwl
thermo 10000
dump 1 all custom 10000 nvt.xyz id type x y z vx vy vz
run 500000
unfix 1
undump 1
reset_timestep 0
write_data nvt.data
write_restart nvt.restart
velocity all scale 300
#------------------set -equilibrium-run-----------
region setl block INF INF INF INF 0 5 units box
region setr block INF INF INF INF 135 140 units box
group setl region setl
group setr region setr
group set2 union setl setr
velocity set2 set 0.0 0.0 0.0
fix 2 set2 nve
fix set set2 setforce 0.0 0.0 0.0
region hot block INF INF INF INF 5 15 units box
region cold block INF INF INF INF 125 135 units box
group hot region hot
group cold region cold
compute Thot all temp/region hot
compute Tcold all temp/region cold
#------------------NVE ---消除内应力---run-----------
fix 1 all nve
fix hot all langevin 320 320 0.01 59804 tally yes
fix cold all langevin 280 280 0.01 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 press vol etotal ke pe evdwl c_Thot c_Tcold f_hot f_cold v_tdiff
thermo 10000
dump 1 all custom 10000 nve.xyz id type x y z vx vy vz
run 500000
undump 1
reset_timestep 0
write_data nve.data
write_restart nve.restart
velocity all scale 300
#---------------set
region cu block INF INF INF INF 20 70 units box
region c block INF INF INF INF 70 120 units box
group cu region cu
group c region c
region set3 block INF INF INF INF 0 5 units box
region setr1 block INF INF INF INF 135 140 units box
group set3 region set3
group setr1 region setr1
group set4 union set3 setr1
region hot1 block INF INF INF INF 5 15 units box
region cold1 block INF INF INF INF 125 135 units box
group hot1 region hot1
group cold1 region cold1
compute Thot1 all temp/region hot1
compute Tcold1 all temp/region cold1
velocity set4 set 0.0 0.0 0.0
fix 6 set4 nve
fix set4 set4 setforce 0.0 0.0 0.0
dump 222222 cu custom 100 cu.txt id type vx vy vz
dump 2222222 c custom 100 c.txt id type vx vy vz
dump_modify 222222 sort id
dump_modify 2222222 sort id
# thermal conductivity calculation
# reset langevin thermostats to zero energy accumulation
variable kB equal 8.625e-5
compute ke all ke/atom
variable temp atom c_ke/1.5/${kB}
fix hot1 all langevin 320 320 0.01 59804 tally yes
fix cold1 all langevin 280 280 0.01 287859 tally yes
fix_modify hot1 temp Thot1
fix_modify cold1 temp Tcold1
dump 3 all custom 10000 all.xyz id type x y z v_temp
variable tdiff equal c_Thot1-c_Tcold1
fix ave all ave/time 10 1000 10000 v_tdiff ave running
thermo_style custom step temp press vol etotal ke pe ebond eangle edihed evdwl c_Thot1 c_Tcold1 f_hot1 f_cold1 v_tdiff f_ave
compute layers all chunk/atom bin/1d z lower 0.025 units reduced
fix 7 all ave/chunk 10 100000 1000000 layers v_temp file Neverybin_Temperature.langevin
发表于 2024-4-25 10:08:48
收藏
分享
淘帖