[转] 热导率计算的in文件 (分子模拟论坛)
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
附件为计算结果:
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的计算结果。
加载中,请稍候......