Hi,
I want to find thermal conductivity of solid argon, and I wrote this code for it:
boundary p p p
atom_style atomic
units lj
#Defining the region and pair styles
neighbor 0.3 bin
neigh_modify delay 1
lattice fcc 0.6
region box block 0 10 0 3 0 3
create_box 1 box
create_atoms 1 region box
mass 1 1.0
pair_style lj/cut 2.5
pair_coeff * * 1.0 1.0 2.5
main groups
region a block INF 1 INF INF INF INF
group heatsource region a
region b block 1 5 INF INF INF INF
group mid1 region b
region c block 5 6 INF INF INF INF
group heatsink region c
region d block 1 2 INF INF INF INF
group 1 region d
region e block 2 3 INF INF INF INF
group 2 region e
region f block 3 4 INF INF INF INF
group 3 region f
region g block 4 5 INF INF INF INF
group 4 region g
region h block 6 INF INF INF INF INF
group mid2 region h
region i block 6 7 INF INF INF INF
group 5 region i
region j block 7 8 INF INF INF INF
group 6 region j
region k block 8 9 INF INF INF INF
group 7 region k
region l block 9 INF INF INF INF INF
group 8 region l
group temp union heatsource heatsink
group mid union mid1 mid2
velocity all create 1 100000 dist gaussian
#Thermalization---------------------------------------------
fix 4 all npt temp .248 .248 .1 iso 0.00238 0.00238 1
thermo_style custom step temp pe etotal press vol
thermo 1000
calculating the temp gradient
compute Tsource heatsource temp
compute T1 1 temp
compute T2 2 temp
compute T3 3 temp
compute T4 4 temp
compute Tsink heatsink temp
compute T5 5 temp
compute T6 6 temp
compute T7 7 temp
compute T8 8 temp
#Averaging and printing Temperature of different groups
fix 8 all ave/time 1 49990 50000 c_Tsource c_T1 c_T2 c_T3 c_T4 c_Tsink c_T5 c_T6 c_T7 c_T8 file Temp.txt
timestep 0.0001
run 600000
unfix 4
#Continue to relax the system with NVE ensemble
fix nve all nve
timestep 0.0001
run 600000
unfix nve
#Creating a hot region and a cold one
fix nve2 mid nve
fix 5 heatsource nvt temp 0.33068 0.33068 0.01
fix 6 heatsink nvt temp 0.16534 0.16534 0.01
fix_modify 5 energy yes
fix_modify 6 energy yes
thermo_style custom step temp pe etotal press vol f_5 f_6
timestep .0001
run 2000000
It’s obvious that first I thermalized the system with NPT ensemble and then I continue to relax the system with NVE. After these steps I made one hot region and one cold region via NVT ensemble and the rest of atoms where integrated with NVE.
The cumulative energy that was added or subtracted in hot and cold region is printed using “f_5” and “f_6”. For finding heat flux I find the slope of cumulative energy vs timestep and then I convert it to energy change vs time.
The slope that I found is 1.510^(-6) that when I want to convert it to energy change vs time gives me 1.510^(-2). For finding heat flux we should divide it to the area that is 9 units. So the heat flux is 0.00167.
The slope of temp gradient that I made is 0.034, which is quite linear. So Delta T/Delta X is equal to 0.034. Now we can find thermal conductivity in dimensionless units.
The thermal conductivity is heat flux/(DeltaT/DeltaX). So thermal conductivity is equal to 0.00167/0.034= 0.049. which is in dimensionless units. Now if we want to convert it to metric units we should divide it by 52.91. which yields 0.000926! The thermal conductivity should be 1~2!( Because the pressure in the last step is about 0.182 in dimensionless units which is about 75 atm, when I have one hot region and one cold region)
Now I have three questions:
1)Why the thermal conductivity that I get is much far different from what there is in references?
2) Why the pressure of the system in the last step ( in which I have one hot region and one cold region becomes higher than I had in NPT ensemble( it becomes about 0.18). It doesn’t make any problem if I use a berendsen barostat in the last step?
3) Can I use compute heat/flux command here? Because the manual has used it just for Green-Kubo method.
Best,
Sina