###########################################################################
# __             ___       __        __   __        __   ______    ______ #
#|  |           / _ \     |  \      /  | |  \      /  | |  __  |  |  ____|#
#|  |          / /_\ \    |   \    /   | |   \    /   | | |__|  | | |____ #
#|  |         /  ___  \   |    \  /    | |    \  /    | |   ___|  |____  |# 
#|  |_____   /  /   \  \  |  |\ \/ /|  | |  |\ \/ /|  | |  |       ____| |#
#|________| /__/     \__\ |__| \__/ |__| |__| \__/ |__| |__|      |______|#
###########################################################################

###########################################################################
################SYMMETRIC TILT GRAIN BOUNDARY ENERGY SCRIPT################
###############This script requires an input tilt direction,###############
##############normal direction, and other varibales in order ##############
#############to create a simulation with two grains stacked on#############
############top of each other with a boundary in in the x-z    ############
###########direction. The lower boundary is then moved relative ###########
##########to the upper boundary in the x and z directions, and   ##########
#########delete overlap distance is increased, all to find the    #########
########lowest energy and thus most stable structure. LAMMPS dumps ########
#######the GB Energy for every configuration it tries, with the     #######
######counter number and GBE as the two numbers in the file name to  ######
#####simplify post processing.                                        #####
###########################################################################

###########################################################################
######################Required Input-Orientations##########################
###########################################################################

#Tilt, normal, and periodic directions are required,
#and must be (at least very close to) mutally orthogonal
#in all three dimensions in order to produce good results


#NOTE: This script will only work with the [0 0 0 1] tilt direction, and
#assumes an initial periodic direction of [2 -1 -1 0] (that is, the periodic
#direction if the misorientation angle is equal to 0).

#Tilt direction, in the form of [t1 t2 t3 t4]

variable t1 equal 0
variable t2 equal 0
variable t3 equal -(${t1}+${t2})
variable t4 equal 1

#Normal direction, in the form of [n1 n2 n3 n4]

variable n1 equal 7
variable n2 equal 10
variable n3 equal -(${n1}+${n2})
variable n4 equal 0

#Periodic direction, in the form of [p1 p2 p3 p4]

variable p1 equal 9
variable p2 equal -8
variable p3 equal -(${p1}+${p2})
variable p4 equal 0

###########################################################################
#########################Required Input-Variables##########################
###########################################################################

variable latticeconst equal 3.181269601
variable caratio equal 1.632993162
variable minimumenergy equal -1.528662571
variable xtranslations equal 5
variable ztranslations equal 5
variable pairstyle index "eam/fs"
variable paircoeff index "* * Mg_mm.eam.fs Mg Mg"

#Desired filename for LAMMPS to create and dump atoms to
variable filename index "minimum_energy"

###########################################################################
#########################Direction transformations#########################
###########################################################################

#Convert directions into xyz coordinates
variable f1 equal ${p1}-0.5*${p2}-0.5*${p3}
variable f2 equal sqrt(3)/2*${p2}-sqrt(3)/2*${p3}
variable f3 equal ${p4}*${caratio}

variable g1 equal ${n1}-0.5*${n2}-0.5*${n3}
variable g2 equal sqrt(3)/2*${n2}-sqrt(3)/2*${n3}
variable g3 equal ${n4}*${caratio}

variable h1 equal ${t1}-0.5*${t2}-0.5*${t3}
variable h2 equal sqrt(3)/2*${t2}-sqrt(3)/2*${t3}
variable h3 equal ${t4}*${caratio}

variable f4 equal ${f1}
variable f5 equal -1*${f2}
variable f6 equal ${f3}

variable g4 equal -1*${g1}
variable g5 equal ${g2}
variable g6 equal ${g3}

variable h4 equal ${h1}
variable h5 equal ${h2}
variable h6 equal ${h3}

###########################################################################
####################Calculate a1 a2 a3 components##########################
###########################################################################

#This will be used for inducing the rotation whenever the custom lattice command
#is implemented. It has to be done in this method instead of "orients" due to 
#magnesium's non-repeating, non-terminating c/a ratio causing LAMMPS to crash
#on its othogonality checks. This method is more forgiving, and any round-off
#errors are negated when the initial minimization is run later.

# Calculate lattice parameters  

variable xlattice equal ${latticeconst}
variable ylattice equal sqrt(3)
variable zlattice equal ${caratio}

variable fmag equal sqrt((${f1})^2+(${f2})^2+(${f3})^2)
variable gmag equal sqrt((${g1})^2+(${g2})^2+(${g3})^2)
variable hmag equal sqrt((${h1})^2+(${h2})^2+(${h3})^2)

#For upper Grain Boundary
variable fx equal ${f1}/${fmag}
variable fy equal ${g1}/${gmag}
variable fz equal ${h1}/${hmag}
variable gx equal ${f2}*${ylattice}/${fmag}
variable gy equal ${g2}*${ylattice}/${gmag}
variable gz equal ${h2}*${ylattice}/${hmag}
variable hx equal ${f3}*${zlattice}/${fmag}
variable hy equal ${g3}*${zlattice}/${gmag}
variable hz equal ${h3}*${zlattice}/${hmag}

#For Lower Grain Boundary
variable fx2 equal ${f4}/${fmag}
variable fy2 equal ${g4}/${gmag}
variable fz2 equal ${h1}/${hmag}
variable gx2 equal ${f5}*${ylattice}/${fmag}
variable gy2 equal ${g5}*${ylattice}/${gmag}
variable gz2 equal ${h2}*${ylattice}/${hmag}
variable hx2 equal ${f6}*${zlattice}/${fmag}
variable hy2 equal ${g6}*${zlattice}/${gmag}
variable hz2 equal ${h3}*${zlattice}/${hmag}

###########################################################################
############################Simulation Parameters##########################
###########################################################################

variable etol equal 1.0e-25 
variable ftol equal 1.0e-25 
variable maxiter equal 10000
variable maxeval equal 100000
variable overlapboth equal 2 
variable counter equal 0 

#Calculate simulation box size

variable xlen equal ${xlattice}*sqrt((${f1})^2+(${f2})^2+(${f3})^2)
variable ylen equal ${xlattice}*sqrt((${g1})^2+(${g2})^2+(${g3})^2)
variable zlen equal ${xlattice}*sqrt((${h1})^2+(${h2})^2+(${h3})^2)

#Find y period length to mitigate the affect of the the 
#complimentary GB caused by y periodic conditions
#Simulations run to determine required length converged 
#on 160 nm or greater

variable mult equal "ceil(160/v_ylen)"
variable ylengb equal "v_ylen * v_mult"
variable negylengb equal "-1*v_ylengb"
print "xlen = ${xlen}, ylen = ${ylengb}, zlen = ${zlen}"

# Implement overlap criterion
variable overlapinc equal 531

#Create directory to dump atoms
shell mkdir ${filename}

###########################################################################
####################Begin Loop Structure###################################
###########################################################################
#This will allow multiple configurations as defined by the
#Xtranslations,Ztranslations, and Overlap distance

label loopa 
variable a loop ${xtranslations} 
variable tx equal "(v_a-1)/v_xtranslations*v_xlen" 
    label loopb 
    variable b loop ${ztranslations} 
    variable tz equal "(v_b-1)/v_ztranslations*v_zlen" 
        label loopd 
        variable d loop ${overlapboth} 
            label loopc 
            variable c loop ${overlapinc} 
            variable overlapdist equal ".500 + 0.005 * (v_c-1)" 
 
#Calculate counter for each loop
            variable ctemp equal ${counter}+1 
            variable counter equal ${ctemp} 
            variable ctemp delete 
            print "Counter: ${counter}" 
            
###########################################################################
#####################Initiallize Simulation################################
###########################################################################

            clear
            units metal
            dimension	3
            boundary	p	p	p
            atom_style	atomic
            atom_modify map array

###########################################################################
#######################Create Lattice Structure############################
###########################################################################

            region	box block 0 ${xlen} ${negylengb} ${ylengb} 0 ${zlen} units box
            create_box	2 box
            lattice custom ${xlattice} a1 ${fx} ${fy} ${fz} a2 ${gx} ${gy} ${gz} a3 ${hx} ${hy} ${hz} basis 0.0 0.0 0.0 basis 0.5 0.5 0 basis 0.5 0.83333333 0.5 basis 0 0.33333333 0.5 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
            region	top block INF INF 0 INF INF INF units box
            create_atoms 1 region top basis 1 1 basis 2 1 basis 3 1 basis 4 1  

            lattice custom ${xlattice} a1 ${fx2} ${fy2} ${fz2} a2 ${gx2} ${gy2} ${gz2} a3 ${hx2} ${hy2} ${hz2} basis 0.0 0.0 0.0 basis 0.5 0.5 0 basis 0.5 0.83333333 0.5 basis 0 0.33333333 0.5 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1 
            region  bottom block INF INF INF 0 INF INF units box
            create_atoms 2 region bottom basis 1 2 basis 2 2 basis 3 2 basis 4 2  
            group upper type 1
            group lower type 2

###########################################################################
#########################Induce Pair Potentials############################
###########################################################################

            pair_style	${pairstyle}
            pair_coeff	${paircoeff}
            neighbor 2.0 bin
            neigh_modify delay 10 check yes

###########################################################################
###############Translate atoms and delete overlapping atoms################
###########################################################################

            displace_atoms upper move ${tx} 0 ${tz} units lattice 
            if $c == 1 then "variable atomprev equal 1" 
            if $d == 1 then "delete_atoms overlap ${overlapdist} lower upper"
            if $d == 2 then "delete_atoms overlap ${overlapdist} upper lower" 

#Ensure that equivalent configuration has not been tested
#If the same number of atoms is present as before, keep looping
#until overlap distance increases enough to delete more atoms

            variable natoms equal "count(all)" 
            print "Previous: ${atomprev}, Present: ${natoms}" 
            if ${atomprev} == ${natoms} then "jump grain_boundary_simplified.in loopend" 
 
###########################################################################
#############################Computes######################################
###########################################################################

            compute csym all centro/atom 
            compute eng all pe/atom 
            compute eatoms all reduce sum c_eng 

###########################################################################
########################Initial Minimization###############################
###########################################################################

            reset_timestep 0 
            thermo 100
            thermo_style custom step pe lx ly lz press pxx pyy pzz c_eatoms
            min_style cg
            minimize 1e-15 1e-15 ${maxiter} ${maxeval}
            
 
#Calculate GB Energy after first minimization

            variable gbarea equal "lx * lz * 2" 
            variable gbe equal "(c_eatoms - (v_minimumenergy * count(all)))/v_gbarea" 
            variable gbemJm2 equal ${gbe}*16021.7733 
            variable gbernd equal round(${gbemJm2}) 
            print "After first minimization:" 
            print "GB energy is ${gbemJm2} mJ/m^2" 

###########################################################################
##############################Second Minimization##########################
###########################################################################

#This time allow box bound to expand and relax while minimizing
#and force convergence tolerance to decrease

            reset_timestep 0
            thermo 100 
            thermo_style custom step pe lx ly lz press pxx pyy pzz c_eatoms 
            min_style cg
            fix 1 all box/relax aniso NULL 0.0 NULL vmax 0.0001
            minimize ${etol} ${ftol} ${maxiter} ${maxeval}
 
 
#Calculate GB Energy after first minimization 

            variable gbarea equal "lx * lz * 2" 
            variable gbe equal "(c_eatoms - (v_minimumenergy * count(all)))/v_gbarea" 
            variable gbemJm2 equal ${gbe}*16021.7733 
            variable gbernd equal round(${gbemJm2}) 
            print "After second minimization:" 
            print "GB energy is ${gbemJm2} mJ/m^2" 

# Store number of atoms for overlap criterion for use in comparing next configurations 
            variable atomprev equal ${natoms}

###########################################################################
#####################Dump data into Data file##############################
###########################################################################

            shell cd ${filename}
            reset_timestep 0 
            timestep 0.001 
            dump 1 all custom 5000 dump.${counter}_${gbernd} id type x y z c_csym c_eng 
            run 0             
            shell cd ..

###########################################################################
#####################Finish Loop Structure#################################
###########################################################################

            label loopend 
            next c 
            jump grain_boundary_simplified.in loopc 
        variable c delete 
        next d 
        jump grain_boundary_simplified.in loopd 
    variable d delete 
    next b 
    jump grain_boundary_simplified.in loopb 
variable b delete 
next a 
jump grain_boundary_simplified.in loopa 
print "All done"