Hello. I am trying to calculate the thermal conductivity along the length and width of graphene nanoribbons. I have been using a system size of 240 atoms (1 X 5 nm with armchair edge), time step of 1 fs, and applying periodic boundary conditions along the length. I essentially took the sample script to calculate the thermal conductivity of argon posted on LAMMPS and modified it for a graphene nanoribbon lattice. I equilibrate the lattice in NVT at 300K for 20,000 times steps, then let it run in NVE for 20,000 time steps before starting the heat flux calculations in NVE. However, it doesn't appear that my simulation converges after 10,000,000 time steps. Do you have any advice on how to make this calculation converge or do you see anything wrong with me code? Here is my code:

#LAMMPS input script to calculate thermal conductivity of graphene

#basic setup

dimension 3

units metal

processors * * 1

boundary p s s

neighbor 3.0 nsq

print "hello"

#define variables

variable T equal 300.0

variable V equal lx*ly*3.35

variable dt equal 0.001

variable Nevery equal 1000

variable Nrepeat equal 10

variable Nfreq equal 1000

#convert from LAMMPS metal units to SI units

variable kB equal 1.3806504e-23

variable eV2J equal 1.602177e-19

variable A2m equal 1.0e-10

variable ps2s equal 1.0e-12

variable convert equal \{eV2J\}\*{eV2J}/\{ps2s\}/{A2m}

#define lattice

atom_style charge

read_data latticetrans.txt

#define interatomic potentials

pair_style rebo

pair_coeff * * /share/apps/lammps-28Mar12/potentials/CH.airebo C

#initialize velocities

velocity all create $T 2502752 mom yes rot yes dist gaussian units box

#define thermo variable settings

timestep ${dt}

thermo ${Nfreq}

#fix ensembles

fix 1 all nvt temp $T $T 0.01 drag 2.0

dump 1 all xyz 100 equil1.xyz

run 20000

undump 1

unfix 1

fix 1 all nve

reset_timestep 0

run 20000

#thermal conductivity calculation

reset_timestep 0

compute myKE all ke/atom

compute myPE all pe/atom

compute myStress all stress/atom virial

compute flux all heat/flux myKE myPE myStress

variable Jx equal c_flux[1]/$V

variable Jy equal c_flux[2]/$V

variable Jz equal c_flux[3]/$V

fix JJ all ave/correlate \{Nevery\} {Nrepeat} ${Nfreq} &

c_flux[1] c_flux[2] c_flux[3] type auto file J0Jt.dat ave running

variable scale equal \{convert\}/{kB}/$T/$T/V\*{Nevery}*${dt}

variable k11 equal trap(f_JJ[3])*${scale}

variable k22 equal trap(f_JJ[4])*${scale}

variable k33 equal trap(f_JJ[5])*${scale}

variable k equal (v_k11+v_k22)/2.0

thermo_style custom step temp v_Jx v_Jy v_k11 v_k22 v_k vol pe ke fnorm

run 50000000

print "average thermal conductivity: $k [W/mK] @ $T K"