#Parameters of simulation and atoms 
units metal 
dimension 3 
boundary p p p
atom_style atomic
newton on off 

#Atom definition 
read_data  graphene_1008atoms

region up block INF INF INF INF -5 5 units box
region down block INF INF INF INF -50 -40 units box
 
pair_style tersoff 
pair_coeff * * SiC.tersoff C 

neighbor 3.5 bin 
#neigh_modify delay 0 every 1 check yes 

timestep 0.0005
 
velocity all create 300.0 31654 mom yes rot yes dist gaussian units box
 
# take some potential energy out of the system 

minimize 1.0e-8 1.0e-6 10000 100000 
#reset_timestep 0 
 
fix thermostat all nvt temp 300.0 300.0 100

compute alltemp all temp 

#perform time-averaging to obtain temperature gradient 

compute ke all ke/atom 

variable temp atom c_ke/(1.5*1.0)/8.617343*100000.0 

fix temp_profile all ave/spatial 1 1000 1000 z -50 50 v_temp file temp.profile units box

fix temp_atom all ave/atom 1 100 100 v_temp 
compute up_temp all temp/region up 
compute down_temp all temp/region down 
variable delta_t equal c_up_temp-c_down_temp

fix delta_t all ave/time 10 10000 100000 v_delta_t file delta_t.dat 

# calculate total kinetic energy 

unfix thermostat

fix ensemble all nve 

fix heat_swap all thermal/conductivity 40 z 10
fix e_exchange all ave/time 40 250 10000 f_heat_swap file e_exchange.dat 


variable thermal_conductivity equal f_heat_swap*1.60217e-19/(2*0.0005*($step+1)*1.0e-12*50*51.77*(f_delta_t+1)*1.0e-10)*100

#2: To address heat flux in periodic boundary conditions
#0.0005 – Timestep
#100000 – Total number of steps in simulation
#1.0e-12: ps to second conversion
#50: Half the length of simulation box along z (i.e. heat flow direction)
#51.77: Total length of simulation box along y direction
#1.0e-10: Angstrom to metre conversion 
#100: Length of simulation box in z direction i.e. in direction of heat flow


fix thermal_conductivity all ave/time 100000 1 100000 v_thermal_conductivity file thermal_conductivity.dat 
thermo_style custom step atoms temp press pe ke etotal xlo xhi ylo yhi zlo zhi vol f_delta_t f_thermal_conductivity
run 900000