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"