[转] 热导率计算的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