Dear all,

I have been studying how to properly calculate the thermal conductivity of a 2d material like graphene with NEMD.

Since it’s the first time that I do this sort of calculation, I started from the example file provided and modified it accordingly. I still have some questions about the procedure LAMMPS uses and some technicalities present in the example file, that I report here in quotes:

**variables**

variable kb equal 1.3806504e-23 #J/K

variable ev2J equal 1.60217e-19

**setup**

dimension 3

boundary p p p #2D material

units metal

atom_style atomic

read_data graphene.data

mass 1 12.0107 # Carbon

velocity all create 300 423425

pair_style tersoff

pair_coeff * * BNC.tersoff C

neighbor 0.3 bin #check

neigh_modify delay 0 every 1 #check

minimize 1e-10 1e-10 10000 100000000

timestep dt

variable dt equal 0.0005

**first equilibration run #NVT for 100ps**

fix 1 all nvt temp 300 300 $(100**dt*)

thermo 100

run 2000000000

velocity all scale 300

unfix 1

**second equilibration #NVE**

compute ke all ke/atom

variable temp atom c_ke*ev2j/1.5kb #convert to SI the T

- Here we want to convert the values to SI and calculate the Temperature in Kelvin using the Boltzman constant. The equation should be E=kbT. So, where does the 1.5 come from?

fix 1 all nve

compute layers all chunk/atom bin/1d x lower 0.05 units reduced #name,on all the atoms,take part of the atoms, slab, origin, thickness

- This should select the thickness of the layer. Here I get confused by the difference with

region hot block 0 30 INF INF INF INF

region cold block lx/2 lx/2+30 INF INF INF INF

in this method one should still consider a cold and an hot region. In this line, I specifically select a slab with the chunk/atom command. Would it be sufficient? If so, why?

fix 2 all ave/chunk 10 100 1000 layers v_temp file profile.mp #Nevery, Nrepeat, Nfreq

fix 3 all thermal/conductivity 1000 x 20 #exchange kinetic energy every n steps, direction, layers in edim direction

- The last number (20), in the manual I read is the layers in the edim direction. It’s not quite clear for me what it means and how one should define these layers. Can you please comment?

variable tdiff equal f_2[11][3]-f_2[1]1[3] #T difference

thermo_style custom step temp epair etotal f_3 v_tdiff

thermo_modify colname f_3 E_delta colname v_tdiff dTemp_step

thermo 1000

run 2000000000

**thermal conductivity calculation**

**reset fix thermal/conductivity to zero energy accumulation**

fix 3 all thermal/conductivity 1000 x 20

- Why is there the need of this step in the calculation?

variable start_time equal time

variable kappa equal (f_3*ev2j/*(time-${start_time}*ps*2)/(ly*A2m*lz*A2m)/2.0)*(lx*A2m/2.0)/f_ave

- Here, since I have a 2d material, the Area on the denomitator of J should be the width times the thickness and everything should be multiplied by the length of the sheet. This means literally the x axis, where the heat flows, right? Moreover, the f_ave should be divided by 2 due to the dual direction of the flow, or this is incorporated in the average?

fix ave all ave/time 1 1 1000 v_tdiff ave running

thermo_style custom step temp epair etotal f_3 v_tdiff f_ave

thermo_modify colname f_3 E_delta colname v_tdiff dTemp_step colname f_ave dTemp

- Here, if I’m using the metal units, I guess I shouldn’t convert anything anymore since I already did in the definition of kappa. Right?

run 20000

print “Running average thermal conductivity: $(v_kappa:%.2f)”

I’m grateful for anybody willing to solve my doubts.