# input script for SWNT simulation in LAMMPS
#-------10--------20--------30--------40--------50--------60--------70--------80

# Variables
variable input index coordi-cnt.txt

# INITIALIZATION
units		    metal
boundary            p p p 
atom_style          atomic

variable	T equal 300.0
variable 	dt equal 0.0002 		# ps
variable    p equal 40000     		# correlation length
variable    s equal 5      		# sample interval
variable	d equal $p*$s 			# dump interval

thermo ${d}
timestep ${dt}            # ps

# ATOM DEFINITION

read_data           ${input}

# FORCE FIELD

pair_style     airebo 3.0 1 1
pair_coeff     * * CH.airebo C
 
# SETTINGS

neighbor	    1.0 bin
neigh_modify	    delay 1


# Calculating Volume

variable 	diameter equal (sqrt(3)*1.42*sqrt(10^2+10*10+10^2))/PI			#(10,10) armchair CNT; assuming C-C bond distance of 1.42
variable 	radius equal ${diameter}/2
variable 	height equal "lz"												#Height is equal to z-dimension length
variable 	V equal 2*PI*${radius}*3.4*${height}							

# Convert from metal units to SI

variable    kB equal 1.3806504e-23    # [J/K] Boltzmann
variable    eV2J equal 1.602176565e-19
variable    A2m equal 1.0e-10
variable    ps2s equal 1.0e-12
variable    convert equal ${eV2J}*${eV2J}/${ps2s}/${A2m}

group Carbon type 1

# ------------------------------------------------------------------------------
variable temp equal temp

# output for VMD
dump 		1 all custom 10000 dump.uhmwpe-coord id type x y z
dump 		2 all xyz 10000 dump.uhmwpe-xyz-format
dump_modify 	2 element C 

## --- Relaxation ---
fix 51 all nvt temp $T $T 0.1
run 200000

## --- Thermal Conductivity Calculation ---
unfix 51
fix 51 all nve

reset_timestep 0
compute      myKE all ke/atom
compute      myPE all pe/atom
compute      myStress all stress/atom NULL 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 $s $p $d &
             c_flux[1] c_flux[2] c_flux[3] type auto file J0Jt.dat ave running
variable     scale equal ${convert}/${kB}/$T/$T/$V*$s*${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}
thermo_style custom step temp cpu v_Jx v_Jy v_Jz v_k11 v_k22 v_k33

run          2000000

variable     k equal (v_k11+v_k22+v_k33)/3.0
print        "average conductivity: $k[W/mK] @ $T K"
unfix 51

print "Run is complete!"