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.5*10^(-6) that when I want to convert it to energy change vs time gives me 1.5*10^(-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