Nbin in fix thermal/conductivity

Dear Arthur,

Thank you for your informative response.

Best regards,

Nima Pirouzmand

Also, I should mention that I have used Green-Kubo method for bilayer graphene and got equal values for components of kappa in x,y and z directions!


Well, then there’s certainly an error in your setup. If you want I can check your input file.


Thank you very much for your favour. Here is my input file:

units metal
boundary p p p
atom_style atomic
#velocity all create 1.35 87287
pair_style hybrid tersoff tersoff lj/cut 6

#read_restart 1.2_relaxed_structure.100000
read_data 1.2_relaxed_structure.dat
#pair_coeff 1 2 lj/cut 0.0046616 3.35
pair_coeff 1 2 lj/cut 0.0048 3.851
pair_coeff * * tersoff 1 SiC.tersoff C NULL
pair_coeff * * tersoff 2 SiC.tersoff NULL C
newton on
timestep 0.001
group one type 1
group two type 2
#variable R equal 0.00198722
variable R equal 8.617343e-5

1st equilibration run

fix 1 all nvt temp 200 200 100 drag 1.0
thermo 100
run 100000

#velocity all scale 1.35

unfix 1

2nd equilibration run

compute ke all ke/atom
variable temp atom c_ke/(1.5*$R)

fix 1 all nve
fix 2 all ave/spatial 10 100 1000 z lower 0.05 v_temp file profile.mp units reduced
fix 3 one ave/spatial 10 100 1000 z lower 0.05 v_temp file one_temperature.dat units reduced
fix 4 two ave/spatial 10 100 1000 z lower 0.05 v_temp file two_temperature.dat units reduced
fix 5 all thermal/conductivity 10 z 20

variable tdiff equal f_4[11][3]-f_3[1][3]
thermo_style custom step temp epair etotal f_5 v_tdiff

thermo 1000
run 100000

thermal conductivity calculation

reset fix thermal/conductivity to zero energy accumulation

fix 5 all thermal/conductivity 10 z 20

fix ave all ave/time 1 1 1000 v_tdiff ave running
thermo_style custom step temp epair etotal f_5 v_tdiff f_ave

run 2000000

Best regards,

Nima Pirouzmand

May be you were talking about the Green-Kubo. Here is my Green-Kubo input to calculate thermal conductivity:
units real
variable T equal 70
variable V equal vol
variable dt equal 1.0
variable p equal 2000 # correlation length
variable s equal 4 # sample interval
variable d equal $p*$s # dump interval

convert from LAMMPS real units to SI

variable kB equal 1.3806504e-23 # [J/K] Boltzmann
variable kCal2J equal 4186.0/6.02214e23
variable A2m equal 1.0e-10
variable fs2s equal 1.0e-15
variable convert equal {kCal2J}*{kCal2J}/{fs2s}/{A2m}

setup problem

dimension 3
boundary p p p
atom_style full
pair_style tersoff
read_data relaxed_structure.dat
pair_coeff * * coeffs.dat C
timestep ${dt}
thermo $d

equilibration and thermalization

velocity all create $T 102486 mom yes rot yes dist gaussian
fix NVT all nvt temp $T $T 10 drag 0.2
run 20000

reset_timestep 0

compute myKE all ke/atom
compute myPE all pe/atom
compute myStress all stress/atom NULL virial

compute flux all heat/flux myKE myPE myStress

variable Jx equal c_flux[1]/vol
variable Jy equal c_flux[2]/vol
variable Jz equal c_flux[3]/vol

fix JJ all ave/correlate $s $p $d c_flux[1] c_flux[2] c_flux[3] type auto file J0Jt.dat ave running

variable scale equal {convert}/{kB}/$T/$T/$V*s*{dt}
variable k11 equal trap(f_JJ[3]){scale} variable k22 equal trap(f_JJ[4])*{scale}
variable k33 equal trap(f_JJ[5])

thermo_style custom step temp v_Jx v_Jy v_Jz v_k11 v_k22 v_k33
run 400000
variable k equal (v_k11+v_k22+v_k33)/3.0
variable kx equal v_k11
variable ky equal v_k22
variable kz equal v_k33
variable ndens equal count(all)/vol
print “Thermal_conductivity in x direction: $kx[W/mK] @ $T K”
print “Thermal_conductivity in y direction: $ky[W/mK] @ $T K”
print “Thermal_conductivity in z direction: $kz[W/mK] @ $T K”
print “Average thermal conductivity: $k[W/mK] @ T K, {ndens} /A^3”

Best regards,

Nima Pirouzmand

I was indeed talking about the Green-Kubo approach. I can see a number of issues:

  • how large is your cell? It has to be pretty large - in fact, large enough so that you see a convergence in the thermal conductivity you obtain
  • the equilibration time (20ps NH nvt) is much too small. Instead of using a Nosé-Hoover thermostat, try using a Berendsen thermostat (+ time integration), which will get you much more quickly to thermal equilibrium
  • did you check total energy conservation (you can also look at the temperature) during the sampling in NVE?
  • the sampling time (400ps) is also too short. I’m not sure for bilayer graphene, but for a single layer at RT, you can expect a correlation time of your heat flux density autocorrelation function of a couple of ns. It’s in general not even enough, and you have to fit in a post-processing step the tail of the autocorrelation function with say a power law, and extrapolate up to the point where your thermal conductivity is converged (to a certain tolerance).


Dear Arthur,

I will work on the comments and let you know the results.Thank you.

Best regards,

Nima Pirouzmand

Another comment: in your script, it looks like you’re only using a Tersoff potential (which is suitable for intra-layer interactions, ie. covalently bonded C atoms) - but then you’re missing inter-layer interactions. If you have bilayer graphene, then if the minimum distance between the two sheets is greater than the C-C Tersoff cutoff distance (and it certainly is), then these two sheets do not interact. You need to model dispersion interaction, for example with a LJ type potential. You should perhaps take a look at REBO like potentials, or COMB3.

I did a very quick and dirty - absolutely unconverged - calculation using bulk graphite, and just a Tersoff potential (so no van der Waals interactions), and here is what I’m getting:

lambda x   300.59 +/- 0.89
lambda y      285 +/- 1.2
lambda z   0.0087 +/- 0.0017

Nothing is converged here, hence the low x and y conductivities. But anyway, you can still see that lambda z is very close to zero, since there are no interactions between the different sheets. 


Thank you Arthur for help and comments.

Best regards,
Nima Pirouzmand