############################################################################################
##
##
##	sample LAMMPS script for EMD (Green-Kubo method) [not fully developed but works]
##	written by Afnan Mostafa
##
##
############################################################################################

dimension 3
units metal
processors * * 2
boundary p p p

neighbor 0.3 bin
neigh_modify every 1 delay 0 check yes

#define variables

variable T equal 300.0
variable V equal lx*ly*10.05  # 3.35 (1st layer) + gap (0) + 3.35 (2nd layer)
variable dt equal 0.0001

variable Nevery equal 3500

variable Nrepeat equal 100

variable Nfreq equal 350000


#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 atomic
read_data graphene_3_layers_ab_10x10.data

#define interatomic potentials

pair_style			airebo 3 1 0
pair_coeff 			* * CH.airebo C C C
#group hydro type 4
#group carbon type 1 2 3
compute 1 all stress/atom NULL
compute 2 all displace/atom
#compute 3 carbon coord/atom cutoff 1.73 group all
#compute 4 hydro coord/atom cutoff 1.3 group carbon 
#variable coords equal "c_3+c_4"

dump mini all custom 5000 minimize.xyz id type x y z c_1[1] c_1[2] c_1[3] c_1[4] c_1[5] c_1[6]
dump_modify mini sort id

#define thermo variable settings

timestep ${dt}

#### BOND FORMATION ####

fix noCOM all recenter INIT INIT INIT units box
thermo_style    	custom step time temp press pxx pyy pzz pe ke etotal vol lx ly lz f_noCOM
thermo          	1000
thermo_modify   	lost error

# ======= minimization =======
min_style		cg
minimize		1e-8 1e-8 400000 400000

#initialize velocities

velocity all create $T 23121 mom yes rot yes dist gaussian units box

# ======= thermal equilibration =======

fix				NVE all nve
fix				ts all temp/rescale 1 $T $T $(100.0*dt) 1.0
run				120000
unfix 		ts
unfix			NVE

fix 1 all nvt temp $T $T $(100.0*dt) drag 2.0
run 120000
unfix 1

fix 1 all nve
run 50000
unfix 1
unfix noCOM
undump mini

dump emd all xyz 100000 emd.xyz
dump_modify emd sort id

fix 1 all nve

#thermal conductivity calculation

reset_timestep 0

thermo ${Nfreq}

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

variable t equal step*${dt}/1000

fix JJ all ave/correlate ${Nevery} ${Nrepeat} ${Nfreq} &
	c_flux[1] c_flux[2] c_flux[3] type auto file J0Jt.data 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 press v_Jx v_Jy v_Jz v_k11 v_k22 v_k33 v_k vol pe ke fnorm lx ly 

fix print_tc all print ${Nfreq} "$t ${k11} ${k22} ${k33} $k" file "./tc_on_the_fly.txt" screen no title "#time_in_ns kxx kyy kzz k"

restart         8750000 tmp.restart

run 87500000

print "average thermal conductivity: $k [W/mK] @ $T K"

write_data post_emd.data