Question about computing thermal conductivity of graphene using fix heat

Dear LAMMPS users,

I am trying to compute the thermal conductivity (kappa) of graphene using LAMMPS using fix heat by changing the corresponding input in examples/KAPPA. Atomic interactions are described by a Tersoff potential (Physical Review B 81, 205441 2010) optimized for phonon thermal transport calculations. From several sources (experimental and simulations), thermal conductivity of graphene is between 2000-4000 W/mK. I’ve obtained ~19 W/mK. After several tests (running the simulations for up to 1 ns, changing the amount of energy added to/removed from the hot/cold layer, increasing the graphene sheet area, etc), the results do not change too much. This makes me guess (if there is nothing wrong with the input below) that the problem is related to the cross section area in the thermal conductivity calculation equation:


I considered the cross section area taking the simulation box lengths given by LAMMPS in the directions perpendicular to the heat flux:


but it is not clear to me if lz can be taken as the thickness of a 2D material such as graphene. In papers reporting this kind of MD simulation this point is overlooked as if it was obvious. Can anyone with previous experience in this kind of simulation help me please?Thanks in advance,Roberto------------------------------------------------------------------------------------------------Input:# ----------------------------- Initialization
dimension 3
units metal
boundary p p s
atom_style atomic

----------------------------- Conversion variables

variable ev2j equal 1.602e-19
variable a2m equal 1.0e-10
variable ps2s equal 1.0e-12
variable convert equal {ev2j}/{ps2s}/${a2m} # [W/m]

----------------------------- Work variables

variable kb equal 8.617e-5
variable dt equal 0.001
variable dk equal 10.0
variable nswap equal 1
variable nbins equal 40
variable nequil equal 100000
variable nrun equal 100000

----------------------------- Atom definition

read_data newdata.read_data
pair_style tersoff
pair_coeff * * BNC.tersoff C# ----------------------------- Define the heat layers
variable width equal lx/{nbins} # Width of the layers variable hotini equal lx/4.0 # Begin: hot layer variable hotend equal lx/4.0+{width} # End: hot layer
variable coldini equal lx3.0/4.0 # Begin: the cold layer
variable coldend equal lx
3.0/4.0+{width} # End: the cold layer region hot block {hotini} {hotend} INF INF INF INF units box region cold block {coldini} ${coldend} INF INF INF INF units box
compute thot all temp/region hot
compute tcold all temp/region cold

----------------------------- Settings

neighbor 2.0 nsq
neigh_modify every 1 delay 10 check yes
mass 1 12.0107

----------------------------- NVE equilibration

thermo 1000
reset_timestep 0
timestep {dt} fix nve all nve fix hot all heat {nswap} {dk} region hot fix cold all heat {nswap} -{dk} region cold thermo_style custom step c_thot c_tcold run {nequil}

----------------------------- Production run

reset_timestep 0
compute ke all ke/atom
variable temp atom c_ke/(1.5*{kb}) compute layers all chunk/atom bin/1d x lower {width} units box
fix 2 all ave/chunk 10 100 1000 layers v_temp ave running file profile.heat
variable tdiff equal f_2[31][3]-f_2[11][3]
fix 3 all ave/time 1 1 1000 v_tdiff ave running start 10000
thermo_style custom step temp c_thot c_tcold v_tdiff f_3
run ${nrun}

----------------------------- Thermal conductivity calculation

variable dT equal f_3
variable dx equal lx/2.0
variable tgrad equal {dT}/{dx}
variable area equal lylz
variable trun equal {nrun}*{dt}
variable dq equal {dk}*{nrun}/{nswap}*{dt}
variable kappa equal -0.5
{dq}/({area}{trun})/{tgrad}{convert} print ">>>>>>> Thermal conductivity: {kappa} [W/mK]"

1 Like