Doubts on the LAMMPS example file

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:

variable kb equal 1.3806504e-23 #J/K
variable ev2J equal 1.60217e-19
dimension 3
boundary p p p #2D material
units metal
atom_style atomic
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

  1. 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

  1. 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 #Nevery, Nrepeat, Nfreq
fix 3 all thermal/conductivity 1000 x 20 #exchange kinetic energy every n steps, direction, layers in edim direction

  1. 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

  1. Why is there the need of this step in the calculation?

variable start_time equal time
variable kappa equal (f_3ev2j/(time-${start_time}ps2)/(lyA2mlzA2m)/2.0)(lx*A2m/2.0)/f_ave

  1. 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

  1. 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.

The average (translational) kinetic energy of an atom at temperature T is 3kT/2 (equipartition theorem)

Thank you! Do you have comments on the other questions?

Sorry I’m not an expert on NEMD so I can’t help too much on that. Others on the forum may have a better idea on it though.