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:
kappa=-0.5dq/(areatime)*dx/dT
I considered the cross section area taking the simulation box lengths given by LAMMPS in the directions perpendicular to the heat flux:
area=ly*lz
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 lx3.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]"