|
|---|
|
|
|
|
|
q********g
发帖数: 10694 | 1 相关搜索: 热导率, 文件, 计算
作者： zhxlhdd2008 于 2010-10-28 16:23
看到有不少人在找热导率计算方面的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 |
|
|
|
|
|
