Finding Thermal conductivity of solid argon

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

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.


Dear Sina

I can’t understand your approach to determining the heat/flux
but i think that you can time average the global variable beforehand:

delE=E - 1.5(N)(Kb)(T)

where E is the total energy of the system contains N atoms.
note that you have to determine delE from the above formula each time step
using appropriate T before the thermostat is applied.

about the heat/flux command why can’t you use that?


Your script is complex and you are trying to convert LJ units
to metric units which could have several subtle issues.
I suggest you first try to compute kapps in LJ units
and verify it is correct wrt literature values.

It is better to use fix heat for your hot/cold regions than
fix nvt. Using fix nvt on regions is conceptually problematic
since the number of atoms in the regions is changing and
the thermostat stores time history.

Using compute heat/flux for kappa is an equilibrium method,
not a non-equil method like you are attempting here. See
Section 6.20 of the manual for a discussion.