#SLBN GK
#####################################################
# Variable Definition ###############################
#####################################################

#package gpu 0

variable		damp equal 0.5
variable    	seed equal 102486 #just a random number
variable    	T equal 300 # Temperature
variable    	V equal vol #volume
variable    	P equal 1.01325 # 1 atm in bar 
variable    	dt equal 0.0005 # 1 fs in metal (timestep)
variable    	r equal 100000 #??
variable    	p equal 10000    # correlation length should be 10000
variable    	s equal 1     	 # sample interval should be 1
variable    	d equal $p*$s  	 # dump interval
variable    	kB equal 1.3806504e-23    # [J/K] Boltzmann constant

# Convert from LAMMPS metal units to SI
variable    	eV2J equal  1.60218e-19 # eV to Joule
variable    	ps2s equal 1.0e-12 # picosec to sec
variable    	A2m equal 1.0e-10 #angstrom to m
variable    	convert equal ${eV2J}*${eV2J}/${ps2s}/${A2m} #

#####################################################
# Initializing simulation box #######################
#####################################################
dimension		3
boundary		p p p
units           metal
atom_style      full
read_data		SLBN50x50.data
replicate		2 2 1
mass			1 10.811  # B
mass			2 14.0067  # N
set				type 1 charge -0.42
set				type 2 charge 0.42
pair_style		tersoff
pair_coeff  	* *		BNC.tersoff B N # chemical
neighbor        2.0 bin
neigh_modify    every 1
neigh_modify    delay 0
neigh_modify    check yes
minimize 1e-10 1e-10  10000 100000000
#delete_atoms overlap 0.1 all all

#####################################################
# 1st equilibration run #############################
######################################################
#At the beginning of each simulation, a system is optimized and then equilibrated
#with NPT(isothermal-isobaric) ensemble under 300K temperature and zero pressure.
#All of the simulations are performed by time step of 1fs with periodic boundary condition(PBC)
#NPT 2.5ns 2500st, NVT 1ns 1000st, NVE 1ns 1000st
timestep  		0.0005
velocity  		all create 300.0 12345 dist gaussian mom yes rot yes
thermo			10000
fix				1 all nvt temp 300 300 0.3 
run 200000
unfix 1
fix				1 all npt temp 300 300 0.3 aniso 1.0 1.0 5
run 500000
unfix 1
fix				1 all nve
run 1000000
#fix				1 all nve
#fix 			2 all press/berendsen aniso 1.0 1.0 500.0 #recommended 1000 timesteps = 1000*0.5fs = 500fs = 0.5
#fix 			3 all temp/berendsen 300.0 300.0 5 #recommended 100 timesteps = 100*0.5fs = 50fs = 0.05
thermo_style	custom step temp press
#run 			100000
#unfix 2
#unfix 3
reset_timestep	0
#####################################################
# Green Kubo Method #################################
#####################################################

thermo       	10000
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]/vol
variable     	Jy equal c_flux[2]/vol
variable     	Jz equal c_flux[3]/vol
fix          	JJ all ave/correlate $s $p $d & 
             	c_flux[1] c_flux[2] c_flux[3] type auto file J0Jt.dat ave running #Nevery Nrepeat Nfreq
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}
variable     	k11r equal trap(f_JJ[3])
variable     	k22r equal trap(f_JJ[4])
variable     	k33r equal trap(f_JJ[5])
thermo_style 	custom step temp v_V v_convert v_k11r v_k22r v_k33r v_scale v_Jx v_Jy v_Jz v_k11 v_k22 v_k33 
dump myDump1 all atom 100 SLG-GK.lammpstrj

run				500000 #should be 100000

variable   	k equal (v_k11+v_k22+v_k33)/3
variable   	ndens equal count(all)/vol
print        	"average conductivity: $k[W/mK] @ $T K, ${ndens} /A\^3"

print		"SIMULATION DONE"



#heatflux = W/m2=J/sm2=x eV/psA2

#V/kBT2= m3K/Js = A3K/eVK2

#conversion factor = A2m / ps2s*K