Dear Arthur,

Thank you for your informative response.

Best regards,

Nima Pirouzmand

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!

Nima

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

Arthur

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

fix 1 all nvt temp 200 200 100 drag 1.0

thermo 100

run 100000

#velocity all scale 1.35

unfix 1

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

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

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}

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

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])*${scale}

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).

Arthur

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.
```

Arthur

Thank you Arthur for help and comments.

Best regards,

Nima Pirouzmand