############################# INITIALIZATION ###############################

     units         metal
     dimension     3
     boundary      p p p
     atom_style     atomic
     atom_modify     map array sort 0 0.0
     neighbor       2.0 bin
     neigh_modify   every 1 delay 0 check yes

#================================================================
#----------------Potential & Structure---------------------------
#================================================================
lattice         fcc 3.604 origin 0.25 0.25 0.25
#region          box block 0 80 0 40 0 40     #bigger block contains Ni
#region          Inconel block 0 40 0 40 0 40   #block of inconel half along X, OTHERS SAME 
region          box block 0 20 0 10 0 10     #bigger block contain Ni
region          Inconel block 0 10 0 10 0 10   #block of inconel half along X, OTHERS SAME
create_box      3 box
create_atoms    1 box

#region  pore sphere 40 20 20 5.5 side in
#delete_atoms region pore

group inconel region Inconel

set   group inconel type/fraction 3 0.47 1234         #Fe
set   group inconel type/fraction 2  0.21 12345       #Cr

mass 1 58.6934 # Ni
mass 2 51.9961 # Cr
mass 3 55.8450 # Fe

pair_style       hybrid/overlay eam/alloy zbl  0.90 1.250
pair_coeff * *   eam/alloy FeNiCr.eam.alloy Ni Cr Fe
pair_coeff 1 1   zbl 28.0 28.0
pair_coeff 2 2   zbl 24.0 24.0
pair_coeff 3 3   zbl 26.0 26.0
pair_coeff 1 2   zbl 28.0 24.0
pair_coeff 1 3   zbl 28.0 26.0
pair_coeff 2 3   zbl 24.0 26.0

# ---------------------------------------------
# MINIMIZATION
# ---------------------------------------------
	 
     # min_style cg
     # timestep 0.001
     # minimize 1.0e-4 1.0e-6 100 1000
     # run 5000

# ---------------------------------------------
# VARIABLE DEFINITIONS
# ---------------------------------------------
variable natoms equal "count(all)"					# Total number of atoms 
print "Number of atoms = ${natoms}"
variable T equal ${natoms}
variable t equal 1
variable naccept equal 0						# Total number of accepted moves
variable increment equal v_naccept+1					# Number of increments of the accepted moves
variable niterations equal 1000000 					# Number of MC iterations	
variable kBT equal 0.086173 						# Boltzmann's constant * Temperature (ambient thermal energy) 
variable h equal 4.1357e-15						# Planck's constant
variable k0 equal 2*${kBT}/$h						# Constant for calculation of total time 
variable ktot equal 0							# Initialization for calculation of total time within each iteration 
variable ttotal equal 0 						# Initialization to calculate total time 
# ----- INITIALIZING SEEDS -----
variable seed equal 27847
variable seed1 equal 27848
variable seed2 equal 27849
variable seed3 equal 27850
variable seed4 equal 27851
variable seed5 equal 27852
variable seed6 equal 27853
# ----- INITIALIZATION OF VARIOUS COUNTERS -----
variable numtype1 equal 0
variable numtype2 equal 0
variable numtype3 equal 0
variable ntype1 equal 0
variable ntype2 equal 0
variable ntype3 equal 0
variable countnn_Ni equal 0
variable countnn_Cr equal 0
variable countnn_Fe equal 0
variable countNi equal 0
variable countCr equal 0
variable countFe equal 0
variable rateA equal 0
variable typeNi equal 1
variable typeCr equal 2
variable typeFe equal 3
# ----- RATE ARRAY -----
variable rateNiNi equal 0.156
variable rateNiCr equal 0.034
variable rateNiFe equal 0.275
variable rateCrCr equal 0.005
variable rateCrFe equal 0.069
variable rateFeFe equal 0.461
# ----- VARIABLE DEFINITION OF THE DISTANCE UP TO THE SECOND NEAREST NEIGHBOUR ATOM -----
variable latt equal 2*2.4890 						# Nearest-neighbour distance between atoms
variable l equal ${latt}						# Assigning the calculated distance to a variable 'l'
# ----- DEFINITION OF RATE TABLE ------
variable r1 equal ${rateNiNi}
variable rate2 equal ${rateNiCr}+${r1}
variable r2 equal ${rate2}
variable rate3 equal ${rateNiFe}+${r2}
variable r3 equal ${rate3}
variable rate4 equal ${rateCrCr}+${r3}
variable r4 equal ${rate4}
variable rate5 equal ${rateCrFe}+${r4}
variable r5 equal ${rate5}
variable rate6 equal ${rateFeFe}+${r5}
variable r6 equal ${rate6}

# ------------------------------------
# GROUP Definitions ATOMS BASED ON TYPE
# ------------------------------------
group Ni type 1 
variable numtype1 equal "count(Ni)" 
variable ntype1 equal ${numtype1}

group Cr type 2
variable numtype2 equal "count(Cr)" 
variable ntype2 equal ${numtype2}

group Fe type 3 
variable numtype3 equal "count(Fe)" 
variable ntype3 equal ${numtype3}

# ---------------------------------------------
# GETTING THE STARTING ENERGY OF THE STRUCTURE
# ---------------------------------------------

variable        e equal pe

run             0

variable emin equal $e

variable elast equal $e

thermo_style custom step v_emin v_elast pe
thermo_modify flush yes

run 0

variable estart equal $e
variable elast equal $e

# ----- DEFINITION OF RATE TABLE ------
#variable r1 equal "(${rateNiNi}/${rate})"
#variable r2 equal "(${rateNiCr}/${rate})"
#variable r3 equal "(${rateNiFe}/${rate})"
#variable r4 equal "(${rateCrCr}/${rate})"
#variable r5 equal "(${rateCrFe}/${rate})"
#variable r6 equal "(${rateFeFe}/${rate})"

# ---------------------------------------------
# MC LOOP
# ---------------------------------------------
label loopa 									# MC Loop label 
variable a loop ${niterations}							# Number of MC iterations			

	print "ITERATION = ${a}" # Iteration
	dump 1 all atom 10000 Inconel-Ni-KMC_Diff.$a.txt 			# Dumping the atom trajectories after every 10000 steps
	dump_modify 1 sort id

	variable randnum equal random(0,1,${seed1})				# Choosing a random number
	variable rand equal ${randnum}

	variable denom equal $t/$T
	print "Denominator = ${denom}"
	variable calc equal ${rand}/${denom}
	print "Calc = ${calc}"	
	variable i equal ${calc}+1						# Choosing the random atom based on the random number
	print "Random atom ID = $i"

	variable boltzfactor equal "exp(atoms*(v_elast - v_e) / v_kBT)" 	# Boltzmann's factor calculation

# ---------------------------------------------
# ASSIGNING THE ID, TYPE AND X,Y,Z COORDINATES OF THE CHOSEN ATOM TO VARIABLES
# ---------------------------------------------
    	variable xi equal x[v_i]						# Extracting the x-coordinate of the atom
    	variable yi equal y[v_i]						# Extracting the y-coordinate of the atom
   	variable zi equal z[v_i]						# Extracting the z-coordinate of the atom
    	variable idi equal id[v_i]						# Extracting the ID of the atom
    	variable typei equal type[v_i]						# Extracting the type of the atom

    	variable atom1 equal ${typei}						# Assigning the atom type to the variable 'atom1'
    	variable ID1 equal ${idi}						# Assigning the atom ID to the variable 'ID_1'
    	variable x0 equal ${xi}							# Assigning the x-coordinate of the atom  to the variable 'x0'
    	variable y0 equal ${yi}							# Assigning the y-coordinate of the atom  to the variable 'y0'	
   	variable z0 equal ${zi}   						# Assigning the z-coordinate of the atom  to the variable 'z0'

	region 1 sphere  ${x0} ${y0} ${z0} $l side in units box
    	group nnatoms region 1 							# Grouping the nearest-neighbour atoms of the randomly-chosen first atom
	
# ---------------------------------------------
# STEPS TO ASSIGN THE GLOBAL PROPERTIES TO THE LOCAL CALCULATED ATOM PROPERTIES
# ---------------------------------------------

	compute neighbour nnatoms property/atom id			
	compute 1 nnatoms chunk/atom c_neighbour compress yes
	fix storechunkid nnatoms store/state 0 c_1
	run 0	
		
# ---------------------------------------------------------------------	
# COUNTING THE NUMBER OF ATOMS WITHIN THE NEAREST-NEIGHBOUR RADIUS
# ---------------------------------------------------------------------	
      	variable num_nnatoms equal "count(nnatoms)"
    	variable nn equal ${num_nnatoms}

	label loopc								# Loop label 
        variable c loop ${nn} 							# Loop to go through all the nearest-neighbour atoms
	
		variable store atom "f_storechunkid==v_c"
	    	group selected variable store
            	variable id_nn equal id[v_nn]					# Extracting the ID of the atom
		
		variable type_nn equal type[v_nn]				# Extracting the type of the atom
	    	variable typ_nn equal ${type_nn} 				# Assigning the first-chosen atom type to the variable 'typ_nn'	

#---------------------------------------------------------------------	
# If the type of the atom is Ni (type 1)/Cr (type 2)/Fe (type 3), count
# the number of each type of atoms in the neighbour group of atoms.
# ---------------------------------------------------------------------	
            	if "${typ_nn}==1" then &					
                	"variable countnn_Ni equal ${countnn_Ni}+1" &		
                	"variable countNi equal ${countnn_Ni}" & 		
            	elif "${typ_nn}==2" &						
                	"variable countnn_Cr equal ${countnn_Cr}+1" &		
                	"variable countCr equal ${countnn_Cr}" &		
            	else &								
                	"variable countnn_Fe equal ${countnn_Fe}+1" &		
                	"variable countFe equal ${countnn_Fe}"			
        next c  
        jump KMC_LATEST.in loopc

# ---------------------------------------------
# Depending on the type of atom chosen and the random numebr chosen, the second atom type is decided.
# A random atom is chosen from the grouped neighbour atoms of the chosen type and the activation energy is chosen.
# The atom type, ID and coordinates of the chosen random atom are extracted and assigned to variables.
# ---------------------------------------------
	
	if "(${atom1}==1) && (${rand}>0) && (${rand}<=${r1}) && (${countnn_Ni}!=0)" then &
		"variable chosen_atomtype equal 1" &
		"print 'First atom = Ni'" &
		"print 'Second atom = Ni'" &
		"variable jrandom equal random(1,${countNi},${seed2})" &
		"variable j equal ${jrandom}" & 
                "variable E equal 1.03369" & 
		"variable E_a equal $E" &
		"variable xj equal x[v_j]" &								
    		"variable yj equal y[v_j]" &				
    		"variable zj equal z[v_j]" &				
    		"variable idj equal id[v_j]" &				
    		"variable typej equal type[v_j]" &				
	    	"variable atom2 equal ${typej}" & 				
    		"variable ID2 equal ${idj}" &				
    		"variable x1 equal ${xj}" &					
    		"variable y1 equal ${yj}" &					
    		"variable z1 equal ${zj}" &
	elif "(${atom1}==1) && (${rand}>${r1}) && (${rand}<=${r2}) && (${countnn_Cr}!=0)" &
		"variable chosen_atomtype equal 2" &
		"print 'First atom = Ni'" &
		"print 'Second atom = Cr'" &
		"variable jrandom equal random(1,${countCr},${seed2})" &
		"variable j equal ${jrandom}" & 
                "variable E equal 1.16538" & 
		"variable E_a equal $E" &
		"variable xj equal x[v_j]" &								
    		"variable yj equal y[v_j]" &				
    		"variable zj equal z[v_j]" &				
    		"variable idj equal id[v_j]" &				
    		"variable typej equal type[v_j]" &				
	    	"variable atom2 equal ${typej}" & 				
    		"variable ID2 equal ${idj}" &				
    		"variable x1 equal ${xj}" &					
    		"variable y1 equal ${yj}" &					
    		"variable z1 equal ${zj}" &
	elif "(${atom1}==1) && (${rand}>${r2}) && (${rand}<=${r3}) && (${countnn_Fe}!=0)" &
		"variable chosen_atomtype equal 3" &
		"print 'First atom = Ni'" &
		"print 'Second atom = Fe'" &
		"variable jrandom equal random(1,${countFe},${seed2})" &
		"variable j equal ${jrandom}" & 
                "variable E equal 0.98484" & 
		"variable E_a equal $E" &
		"variable xj equal x[v_j]" &								
    		"variable yj equal y[v_j]" &				
    		"variable zj equal z[v_j]" &				
    		"variable idj equal id[v_j]" &				
    		"variable typej equal type[v_j]" &				
	    	"variable atom2 equal ${typej}" & 				
    		"variable ID2 equal ${idj}" &				
    		"variable x1 equal ${xj}" &					
    		"variable y1 equal ${yj}" &					
    		"variable z1 equal ${zj}" &
	elif "(${atom1}==2) && (${rand}>${r1}) && (${rand}<=${r2}) && (${countnn_Ni}!=0)" &
		"variable chosen_atomtype equal 1" &
		"print 'First atom = Cr'" &
		"print 'Second atom = Ni'" &
		"variable jrandom equal random(1,${countNi},${seed2})" &
		"variable j equal ${jrandom}" & 
                "variable E equal 1.16538" & 
		"variable E_a equal $E" &
		"variable xj equal x[v_j]" &								
    		"variable yj equal y[v_j]" &				
    		"variable zj equal z[v_j]" &				
    		"variable idj equal id[v_j]" &				
    		"variable typej equal type[v_j]" &				
	    	"variable atom2 equal ${typej}" & 				
    		"variable ID2 equal ${idj}" &				
    		"variable x1 equal ${xj}" &					
    		"variable y1 equal ${yj}" &					
    		"variable z1 equal ${zj}" &
	elif "(${atom1}==2) && (${rand}>${r3}) && (${rand}<=${r4}) && (${countnn_Cr}!=0)" &
		"variable chosen_atomtype equal 2" &
		"print 'First atom = Cr'" &
		"print 'Second atom = Cr'" &
		"variable jrandom equal random(1,${countCr},${seed2})" &
		"variable j equal ${jrandom}" & 
                "variable E equal 1.33551" & 
		"variable E_a equal $E" &
		"variable xj equal x[v_j]" &								
    		"variable yj equal y[v_j]" &				
    		"variable zj equal z[v_j]" &				
    		"variable idj equal id[v_j]" &				
    		"variable typej equal type[v_j]" &				
	    	"variable atom2 equal ${typej}" & 				
    		"variable ID2 equal ${idj}" &				
    		"variable x1 equal ${xj}" &					
    		"variable y1 equal ${yj}" &					
    		"variable z1 equal ${zj}" &
	elif "(${atom1}==2) && (${rand}>${r4}) && (${rand}<=${r5}) && (${countnn_Fe}!=0)" &
		"variable chosen_atomtype equal 3" &
		"print 'First atom = Cr'" &
		"print 'Second atom = Fe'" &
		"variable jrandom equal random(1,${countFe},${seed2})" &
		"variable j equal ${jrandom}" & 
                "variable E equal 1.10366" & 
		"variable E_a equal $E" &
		"variable xj equal x[v_j]" &								
    		"variable yj equal y[v_j]" &				
    		"variable zj equal z[v_j]" &				
    		"variable idj equal id[v_j]" &				
    		"variable typej equal type[v_j]" &				
	    	"variable atom2 equal ${typej}" & 				
    		"variable ID2 equal ${idj}" &				
    		"variable x1 equal ${xj}" &					
    		"variable y1 equal ${yj}" &					
    		"variable z1 equal ${zj}" &
	elif "(${atom1}==3) && (${rand}>${r2}) && (${rand}<=${r3}) && (${countnn_Ni}!=0)" &
		"variable chosen_atomtype equal 1" &
		"print 'First atom = Fe'" &
		"print 'Second atom = Ni'" &
		"variable jrandom equal random(1,${countNi},${seed2})" &
		"variable j equal ${jrandom}" & 
                "variable E equal 0.98484" & 
		"variable E_a equal $E" &
		"variable xj equal x[v_j]" &								
    		"variable yj equal y[v_j]" &				
    		"variable zj equal z[v_j]" &				
    		"variable idj equal id[v_j]" &				
    		"variable typej equal type[v_j]" &				
	    	"variable atom2 equal ${typej}" & 				
    		"variable ID2 equal ${idj}" &				
    		"variable x1 equal ${xj}" &					
    		"variable y1 equal ${yj}" &					
    		"variable z1 equal ${zj}" &
	elif "(${atom1}==3) && (${rand}>${r4}) && (${rand}<=${r5}) && (${countnn_Cr}!=0)" &
		"variable chosen_atomtype equal 2" &
		"print 'First atom = Fe'" &
		"print 'Second atom = Cr'" &
		"variable jrandom equal random(1,${countCr},${seed2})" &
		"variable j equal ${jrandom}" & 
                "variable E equal 1.10366" & 
		"variable E_a equal $E" &
		"variable xj equal x[v_j]" &								
    		"variable yj equal y[v_j]" &				
    		"variable zj equal z[v_j]" &				
    		"variable idj equal id[v_j]" &				
    		"variable typej equal type[v_j]" &				
	    	"variable atom2 equal ${typej}" & 				
    		"variable ID2 equal ${idj}" &				
    		"variable x1 equal ${xj}" &					
    		"variable y1 equal ${yj}" &					
    		"variable z1 equal ${zj}" &
	elif "(${atom1}==3) && (${rand}>${r5}) && (${rand}<=${r6}) && (${countnn_Fe}!=0)" &
		"variable chosen_atomtype equal 2" &
		"print 'First atom = Fe'" &
		"print 'Second atom = Fe'" &
		"variable jrandom equal random(1,${countFe},${seed2})" &
		"variable j equal ${jrandom}" & 
                "variable E equal 0.9404" & 
		"variable E_a equal $E" &
		"variable xj equal x[v_j]" &								
    		"variable yj equal y[v_j]" &				
    		"variable zj equal z[v_j]" &				
    		"variable idj equal id[v_j]" &				
    		"variable typej equal type[v_j]" &				
	    	"variable atom2 equal ${typej}" & 				
    		"variable ID2 equal ${idj}" &				
    		"variable x1 equal ${xj}" &					
    		"variable y1 equal ${yj}" &					
    		"variable z1 equal ${zj}" &
	else &
		"variable E equal 0" & 
		"variable E_a equal $E"	

	print "Activation energy = ${E_a}"

# --------------------------------------------------------------------------------------------	
# UPDATING THE POSITIONS OF THE SWITCHED ATOMS    
# --------------------------------------------------------------------------------------------	
	
	if "${E_a}!=0" then &
		"set atom $j x ${x0}" & 
    		"set atom $j y ${y0}" &
  		"set atom $j z ${z0}" &
		"set atom $i x ${x1}" &
    		"set atom $i y ${y1}" &
    		"set atom $i z ${z1}"
	
	run 1 pre no post no						# Run once to get the energy of the present structure

# ---------------------------------------------
# METROPOLIS CRITERION TO CALCULATE TIME ADVANCED
# ---------------------------------------------

	if "${E_a}==0" then &
		"variable k equal ${k0}*exp(-${E_a}/${kBT})" &				
		"variable Rtot equal ${ktot}+${k}"	&				
		"variable timerandno equal random(0.0,1.0,${seed4})" &			
        	"variable deltat equal -log(${timerandno})/${Rtot}" &				
		"variable ttotal equal ${ttotal}+${deltat}" &
	elif "(${E_a}!=0) && ($e<=${elast})" &
		"variable elast equal $e" &
		"variable k equal ${k0}*exp(-${E_a}/${kBT})" &				
		"variable Rtot equal ${ktot}+${k}" &				
		"variable timerandno equal random(0.0,1.0,${seed4})" &			
        	"variable deltat equal -log(${timerandno})/${Rtot}" &				
		"variable ttotal equal ${ttotal}+${deltat}" &																		
	elif "(${E_a}!=0) && (${rand} <= ${boltzfactor})" &														
             	"variable elast equal $e" &
		"variable k equal ${k0}*exp(-${E_a}/${kBT})" &				
		"variable Rtot equal ${ktot}+${k}"&				
		"variable timerandno equal random(0.0,1.0,${seed4})" &			
        	"variable deltat equal -log(${timerandno})/${Rtot}" &				
		"variable ttotal equal ${ttotal}+${deltat}" &         
        else & 									
             	"set atom $i x ${x0}" &
             	"set atom $i y ${y0}" &
             	"set atom $i z ${z0}" &
             	"set atom $j x ${x1}" &
             	"set atom $j y ${y1}" &
             	"set atom $j z ${z1}" &
		"variable k equal ${k0}*exp(-${E_a}/${kBT})" &				
		"variable Rtot equal ${ktot}+${k}" &				
		"variable timerandno equal random(0.0,1.0,${seed4})" &			
        	"variable deltat equal -log(${timerandno})/${Rtot}" &				
		"variable ttotal equal ${ttotal}+${deltat}" 

# ----- Print the calculated advanced time (deltat) and the total time of the process after every MC iteration (ttotal) -----
      
	print "delta_t = ${deltat}"
	print "Total time = ${ttotal}"

	run 100000
	undump 1
		
	uncompute neighbour
	uncompute 1
	unfix storechunkid
	region 1 delete	
	group nnatoms delete
	group selected delete

# ----- Next iteration in the Monte Carlo loop -----
	
	next a
	jump KMC_LATEST.in loopa

# ----------------------------------------------------------------------------------
    
print "Total time of simulation = ${ttotal}" 					# Display total time of simulation