calculate thermal conductivity of graphene

Dear All,

I am trying to calculate the thermal conductivity of graphene with LAMMPS. I am using the Muller-Plathe method.
But the value I got is 2 orders larger than the reported value( 4.8410^3 ~ 5.3 10^3 W/mK, Nano Lett., 2008, 8(3), pp902-907). I checked the temperature distribution on my system
and it looks quite linear. I really appreciate if anyone could give me some advice on this. Pleae refer to my code
posted here:

units metal
boundary p p p
atom_style molecular
lattice custom 2.4595 a1 1 0 0 a2 0 1.73203 0 a3 0.0 0.0 10 basis 0 0 0 basis 0.5 0.16666666666666666 0 basis 0.5 0.5 0 basis 0 0.6666666666666666 0
#YL heat flow is in the y direction; z is perpendicular to the single layer graphene; graphene is in the x-y plane
region box block 0 80 0 120 -0.5 0.5
create_box 1 box
create_atoms 1 box
region coldest block INF INF 0 1.5 INF INF
region hottest block INF INF 60 61.5 INF INF
mass 1 12

pair_style tersoff
pair_coeff * * SiC_1994.tersoff C
neighbor 2.0 bin
neigh_modify delay 5

dump snap all cfg 10000 dump.config.*.cfg id type xs ys zs vx vy vz


velocity all create 600 458127641 mom yes rot yes dist gaussian
fix NPT all npt temp 300 300 1 iso 0 0 0.1 drag 1
thermo_style custom step temp etotal vol lx ly lz press pxx pyy pzz
thermo_modify lost warn
thermo 1000

timestep 0.001
run 50000
write_restart restart.NPT50k.equil
reset_timestep 0
unfix NPT

fix temp all temp/berendsen 300 300 0.1
fix nve all nve
variable kb equal 1.3806504e-23
variable ev2J equal 1.60217e-19
compute ke all ke/atom
variable temp atom c_ke*{ev2J}/(1.5*{kb})
fix temp_profile all ave/spatial 1 100000 100000 y lower 1.5 v_temp file temp.profile units lattice #YL for test
#calculate temperature difference***************************
compute hottest_temp all temp/region hottest
compute coldest_temp all temp/region coldest
variable delta_temp equal c_hottest_temp-c_coldest_temp
fix delta_out all ave/time 1 100000 100000 v_delta_temp file delta_temp.dat

fix hottesttemp_out all ave/time 1 100000 100000 c_hottest_temp file hottesttemp.dat
fix coldesttemp_out all ave/time 1 100000 100000 c_coldest_temp file coldesttemp.dat
run 100001
unfix temp

fix heat_swap all thermal/conductivity 10 y 80
fix e_exchange all ave/time 10 1000 100000 f_heat_swap file e_exchange.dat
variable ps2s equal 1.0e-12
variable An2m equal 1.0e-10

variable accumu_thermalK equal f_e_exchange*{ev2J}/(4*0.001*10*{ps2s}196.7624.595*${An2m}*f_delta_out)*511.191
#0.001 is the timestep
#5 is the swap frequency
#196.76 is the width of the graphene layer
#1.42 is the thickness of the graphene layer
#511.191 is the length of the graphene layer (also the heat flux direction)

fix accumu_thermalK_out all ave/time 100000 1 100000 v_accumu_thermalK file thermal_conductivity.dat

the thermal conductivity calculated on step i should be

accumu_thermalK_out (i) minus accumu_thermalK_out (i-1)

run 2000000

Thank you very much.

Best Regards,

Dear Leon

Are you modeling isotopically pure graphene, i.e. do all your carbon atoms have the same mass? The thermal conductivity is determined by the mean path length of the phonons, which in turn is determined by impurities in the material – even the difference in mass between C12 and C13 has a large effect, particularly at low temperatures. While it may not be all of the difference, I think at room temperature, using isotopically pure C might enhance the thermal conductivity my a factor of 10. Also, the correctly calculated thermal conductivity will always tend to be higher than experimental measurements due to defects in the real material, for example there may be some N atoms, and perhaps 5 and 7 membered rings in graphene.


[email protected]…1795…

Phone: +1.760.495.4924 x 102
Cell: +1.505.603.8182

This email is covered by the Electronic Communications Privacy Act, 18 U.S.C. Section 2510-2521. This message and any attachments hereto may contain confidential information intended only for the use of the individual or entity named above. If you are not the intended recipient(s), or the employee or agent responsible for delivery of this message to the intended recipient(s), you are hereby notified that any dissemination, distribution or copying of this email message is strictly prohibited. If you have received this message in error, please immediately notify the sender and delete this email from your computer. The sender does not waive any privilege in the event this message was inadvertently disseminated.

MedeA® and Materials Design® are registered trademarks of Materials Design, Inc.

I would like to amplify on what Paul said, and I have to say that I have been wondering when this was going to happen. Transport properties usually depend on spatial dimensionality. Thermal conductivity is among those that are dimensionally dependent. If you had

boundary p p f # perfect 2d system

then Alder and Wainwright (Phys Rev 1970) predict infinite conductivities. The prediction has been tested with hard spheres and discs. Below 3 spatial dimensions, our normal concepts of transport go haywire.

Graphene sheets are essentially two-dimensional. As it stands, your “3d” simulation only has a little thermal noise to produce scattering that gives the full 3d dynamics needed to calculate a finite thermal conductivity. The experimental guys likely have grain boundaries and other defects to help them. (Speculating here, but Paul seems to agree).

That’s how things look to me.

Of course, it could be LAMMPS-centric things like a terrible potential (or an error), scale-factor in your formula for accumu_thermalK, but I doubt any of these are to blame.

Steve Valone


Your and Paul’s comments help me a lot.

Best Regards,


Thank you very much.

Best Regards,