24小时热门版块排行榜

返回列表

当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

zxzj05

荣誉版主 (著名写手)

应助: 0 (幼儿园)
贵宾: 1.242
金币: 5912.5
红花: 2
帖子: 1528
在线: 144.5小时
虫号: 406147
注册: 2007-06-18
性别: GG
专业: 金属功能材料
管辖: 第一性原理

[交流] 【转帖】热导率计算的in文件已有12人参与

看到有不少人在找热导率计算方面的in文件，我就贡献三个in文件吧，仅供参考。

同时，附件里贴出了我的计算结果。EMD的输出结果（compute heat/flux command+compute tc command的计算结果）中， “ac.dat”（见附件中的"ac.wmf"）是热流自相关函数（我已经修改了compute_tc.cpp，目前输出的是normalized HCACF，但结果中给出的还是没有归一化的热流自相关函数，但形状和归一化的是一样的，请注意！）随m的变化，"tc.dat"(见附件中的"tc.wmf"

是热导率随m的变化（m的涵义请参看热导率计算的Green Kubo离散化公式,见附件"Comparison of atomic-level simulation methods for computing thermal conductivity”中的(9)式），"tc_time.dat"(见附件中的"tc_time.wmf"

是热导率随时间的变化；NEMD的结果中，"temp.profile"(见附件中的"temp_distribution.wmf"

是z方向上的温度分布，"thermal_conductivity.dat“(见附件中的"NEMD.wmf"

表示热导率随时间的变化，两者都包含了fix heat command 和fix thermal/conductivity command的计算结果。

1.用compute heat/flux command+compute tc command得到热流自相关函数和热导率（EMD方法）

# MD simulation of Ar thermal conductivity
# Initialization
units          lj
dimension       3
newton          on
boundary       p p p
atom_style       atomic
neighbor       0.3 bin
neigh_modify    check  yes
lattice          fcc 0.844
region          box  block -4  4  -4  4  -4  4  units lattice
create_box    1  box
create_atoms    1  box
mass    1  1.0
velocity       all  create  0.71 458127641 mom yes  rot yes dist gaussian units box
# LJ potential *********************************************************
pair_style       lj/cut 2.8
pair_coeff       1  1 1.0                1.0                   #  LJ parameters for Ar-Ar
fix             temp all  temp/berendsen 0.71 0.71 0.000466
fix             nve  all  nve
thermo_style    custom step temp etotal vol
thermo_modify    lost warn
thermo          100
# Run

timestep       0.000466
run             200000
reset_timestep 0
# -------------- Flux calculation in nve ---------------
compute    myKE all ke/atom
compute    myPE all pe/atom
compute    myStress all stress/atom virial
variable       factor_ac equal 1.0
variable       factor_tc equal 1.3806504e-23*sqrt(1.67e-21/6.633e-26)/3.405e-10^2
compute    jflux all heat/flux myKE myPE myStress
compute          tc all tc c_thermo_temp c_jflux v_factor_ac v_factor_tc iso first 10000 900000 100000
fix             tc_out  all  ave/time  1  1  1  c_tc file  tc_time.dat
thermo_style    custom  step  temp
restart          100000 restart.*
run          1000000

2. 用fix thermal/conductivity command得到温度梯度，进而得到热导率(NEMD方法)

# MD simulation of Ar thermal conductivity
# Initialization
units          lj
dimension       3
newton          on
boundary       p p p
atom_style       atomic
neighbor       0.3 bin
neigh_modify    check  yes
lattice          fcc 0.844
region          box  block -4  4  -4  4  -4  4  units lattice
create_box    1  box
create_atoms    1  box
region          up1 block  INF INF  INF INF  -0.5  -0.25  units lattice
region          up2 block  INF INF  INF INF 0.5 0.75  units lattice
region          up  union 2 up1 up2
region          down1  block  INF INF  INF INF  -3.5  -3.25  units lattice
region          down2  block  INF INF  INF INF 3.5 3.75  units lattice
region          down union 2 down1 down2
mass    1  1.0

velocity       all  create  0.71 458127641 mom yes  rot yes dist gaussian units box
# Tersoff potential *********************************************************
pair_style       lj/cut 2.8
pair_coeff       1  1 1.0                1.0                   #  LJ parameters for Ar-Ar
fix             temp all  temp/berendsen 0.71 0.71 0.0466
fix             nve  all  nve
compute          ke  all  ke/atom
variable       temp atom  c_ke/(1.5*1.0)
fix             temp_profile all ave/spatial  1  100000  100000  z  lower  0.25    v_temp  file  temp.profile units  lattice
compute          up_temp all  temp/region up
compute          down_temp  all  temp/region down
variable       delta_temp equal  c_up_temp-c_down_temp
fix             delta_out  all  ave/time  1  100000  100000  v_delta_temp file  delta_temp.dat
thermo_style    custom step temp etotal vol
thermo_modify    lost warn
thermo          100
# Run

timestep       0.000466
run             100001
unfix          temp
fix             heat_swap all  thermal/conductivity  10  z  32
fix             e_exchange  all  ave/time  10  10000  100000  f_heat_swap  file  e_exchange.dat
variable       thermal_conductivity equal f_e_exchange/(0.000466*10.0*4.0*f_delta_out)*1.3806504e-23/3.405e-10/3.405e-10*sqrt(1.67e-21/6.633e-26)*6.0/8.0

# 以上variable命令需要特别注意，因为我所模拟的系统，盒子边长Lx=Ly=Lz，热导率计算公式经过推导变成为e_exchange/(4.0*t*L*delta_T)，
# 为了不在in文件里给L赋值，我修改了fix_thermal_conductivity.cpp文件(见附件)，将e_exchange修改成了 e_exchange += force->mvv2e * (all[0].value - all[1].value) / (domain->zprd); 同时在end_of_step()
里添加了一句 “ e_exchange = 0.0;“,详见附件中的fix_thermal_conductivity.cpp，这样所得的e_exchange曲线基本上是一条水平曲线，而不是用原来的fix thermal/conductivity command所得到的斜向上的曲线，请注意！！！
# 所以才出现以上variable的表达式。
# 请看明白后再做计算，免得算出错误的结果！！！

fix             thermal_conductivity_out  all  ave/time  100000  1 100000  v_thermal_conductivity file  thermal_conductivity.dat

# Run
run             10000000

3. 用fix heat command建立温度梯度，进而得到热导率（NEMD方法）

# MD simulation of Ar thermal conductivity
# Initialization
units          lj
dimension       3
newton          on
boundary       p p p
atom_style       atomic
neighbor       0.3 bin
neigh_modify    check  yes
lattice          fcc 0.844
region          box  block -4  4  -4  4  -4  4  units lattice
create_box    1  box
create_atoms    1  box
region          up1 block  INF INF  INF INF  -0.5  -0.25  units lattice
region          up2 block  INF INF  INF INF 0.5 0.75  units lattice
region          up  union 2 up1 up2
region          down1  block  INF INF  INF INF  -3.5  -3.25  units lattice
region          down2  block  INF INF  INF INF 3.5 3.75  units lattice
region          down union 2 down1 down2
region          hot  block  INF INF  INF INF    0.0 0.25  units lattice
group          hot  region  hot
region          cold  block  INF INF  INF INF -4.0  -3.75  units lattice
group          cold  region  cold
mass    1  1.0
mass0          6.633e-26
epsilon0       1.67e-21
sigma0          3.405e-10
velocity       all  create  0.71 458127641 mom yes  rot yes dist gaussian units box
# Tersoff potential *********************************************************
pair_style       lj/cut 2.8
pair_coeff       1  1 1.0                1.0                   #  LJ parameters for Ar-Ar
fix             temp all  temp/berendsen 0.71 0.71 0.0466
fix             nve  all  nve
compute          ke  all  ke/atom
variable       temp atom  c_ke/(1.5*1.0)
fix             temp_profile all ave/spatial  1  100000  100000  z  lower  0.25    v_temp  file  temp.profile units  lattice
compute          up_temp all  temp/region up
compute          down_temp  all  temp/region down
variable       delta_temp equal  c_up_temp-c_down_temp
fix             delta_out  all  ave/time  1  100000  100000  v_delta_temp file  delta_temp.dat
thermo_style    custom step temp etotal vol
thermo_modify    lost warn
thermo          100
# Run

timestep       0.000466
run             100001
unfix          temp
fix             hot all  heat  1 50 region hot
fix             cold  all  heat  1  -50 region cold
variable       thermal_conductivity equal 50.0*0.5*1.67e-21/3.405e-10/sqrt(6.633e-26/1.67e-21)/((4.0*8.0*8.0*8.0/0.844)^(1.0/3.0)*3.405e-10*2.0*f_delta_out*1.67e-21/1.3806504e-23)*6.0/8.0
fix             thermal_conductivity_out  all  ave/time  100000  1 100000  v_thermal_conductivity file  thermal_conductivity.dat

# Run
run             10000000

【转自】http://www.mdbbs.org/thread-26451-1-1.html

回复此楼