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