#BNGR GK
#####################################################
# Variable Definition ###############################
#####################################################
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.001 # 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
#package omp 8
package kokkos comm host
suffix kk
#suffix hybrid omp kk

# 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}
# Initialization
units           metal
boundary        p p p
atom_style      full
#processors      * * 1     # domain decomposition over x and y

# System and atom definition
# we use different molecule ids for each layer
# so that inter- and intra-layer
# interactions can be specified separately
read_data       gr_hBN_Cstack_2L_noH.data
#replicate		2 2 1
mass            1 10.8110   # boron    mass (g/mole) | membrane
mass            2 14.0067   # nitrogen mass (g/mole) | membrane
mass            3 12.0107   # carbon   mass (g/mole) | adsorbate
# Separate atom groups
#group 		hBN type 1 #1 hBN 2 C
#group 		hBN type 2
#group 		gr  type 3
group 		hBN molecule 1
group 		gr  molecule 2
# within hBN  - 11, 22, 12 within C - 33 cross hBN/C - 13, 23
######################## Potential defition ########################
pair_style  hybrid/overlay rebo tersoff ilp/graphene/hbn 16.0 coul/shield 16.0
####################################################################
pair_coeff  * * rebo               CH.rebo     NULL NULL C    # chemical
pair_coeff  * * tersoff            BNC.tersoff B    N    NULL # chemical
pair_coeff  * * ilp/graphene/hbn   BNCH.ILP    B    N    C    # long range
pair_coeff  1 1 coul/shield 0.70
pair_coeff  1 2 coul/shield 0.69498201415576216335
pair_coeff  2 2 coul/shield 0.69

#pair_coeff  1 1 coul/shield 0.70 #within hBN
#pair_coeff  1 2 coul/shield 0.70 #within hBN
#pair_coeff  2 2 coul/shield 0.70 #within hBN
#pair_coeff  2 3 coul/shield 0.69498201415576216335 #cross hBN/C
#pair_coeff  1 3 coul/shield 0.69498201415576216335 #cross hBN/C
#pair_coeff  3 3 coul/shield 0.69 #within C
####################################################################
# Neighbor update settings
neighbor        2.0 bin
neigh_modify    every 1
neigh_modify    delay 0
neigh_modify    check yes
#minimize 1e-10 1e-10  10000 100000000
#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
run_style		verlet
thermo			100000
velocity     	all create 300 12345 mom yes rot yes dist gaussian
fix				1 all npt temp 300.0 300.0 0.2 iso 0.0 0.0 0.1
thermo_style	custom step temp press
run 			2500
unfix			1
fix				1 all nvt temp 300.0 300.0 0.2
thermo_style	custom step temp press
run 			300000
unfix			1
fix				1 all nve
thermo_style	custom step temp press
run 			300000
#unfix			1
reset_timestep	0
#####################################################
# Green Kubo Method #################################
#####################################################

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 #calculates HCACF J=V/kBT^2int(J<0>.J<t>)#Nevery Nrepeat Nfreq
variable     	scale equal ${convert}/${kB}/$T/$T/$V*$s*${dt} #converts V/kBT^2
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 v_Jx v_Jy v_Jz v_k11 v_k22 v_k33
dump          	2 all custom 100 BNGr-GK.lammpstrj id element x y z

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

print		"SIMULATION DONE"