# Initialize the problem
units metal
dimension 3
boundary p p p
atom_style atomic
variable latparam equal 4.12579125166

# Declare the dimensions of the model in lattice cells
variable nlat_x equal 1
variable nlat_y equal 1
variable nlat_z equal 58

# Define the box/various dimensions
variable xdim      equal (${latparam}*sqrt(6)/2)*${nlat_x}
variable ydim      equal (${latparam}*sqrt(2)/2)*${nlat_y}
variable zdim      equal (${latparam}/sqrt(3))*${nlat_z}
variable layer_z   equal (${latparam}/sqrt(3))

# Area of the stacking fault plane in Angstrom^2
variable Area equal 2*${xdim}*${ydim}

# Geometry definition
lattice fcc ${latparam}

# Define regions for stacking fault test
# NOTE: -0.001 in the x- and y- directions only to overcome numerical precision issues

# Define a region for each  of the top N_layers/2
variable ntwin_layers equal ${nlat_z}/2

variable i loop ${ntwin_layers}
label start_of_loop_region

   variable j equal ${i}+${ntwin_layers}
   variable zmin equal (${j}-1.0)*${layer_z}-0.001
   variable zmax equal ${j}*${layer_z}-0.001
   region  ${i}  block  -0.001 ${xdim} -0.001 ${ydim}  ${zmin} ${zmax} units box
   next i

jump SELF start_of_loop_region

variable i delete
variable j delete

# Declare group for the rigid blocks for intrinsic stacking fault
variable zmin equal (15-1.0)*${layer_z}-0.001
variable zmax equal (45-1.0)*${layer_z}+0.001
region  stack_region block  -0.001 ${xdim} -0.001 ${ydim}  ${zmin} ${zmax} units box

# Declare group for the rigid blocks for the extrisic stacking fault
variable zmin equal (16-1.0)*${layer_z}-0.001
variable zmax equal (44-1.0)*${layer_z}+0.001
region  twin_region block  -0.001 ${xdim} -0.001 ${ydim}  ${zmin} ${zmax} units box

# Create simulation box and atoms
region       whole block  0-0.001  ${xdim}+0.01    0-0.01 ${ydim}+0.01    0-0.01  ${zdim}+0.01   units box
create_box   1     whole
lattice fcc ${latparam} orient x 1 1 2 orient y -1 1 0 orient z -1 -1 1
create_atoms 1 region whole

# Define the interatomic potentials
pair_style kim LJ_ElliottAkerson_2015_Universal__MO_959249795837_003
pair_coeff * * Rn
mass * 1.0

# Rigid holding settings for each layer as a group
variable i loop ${ntwin_layers}
label start_of_loop_group

   group ${i} region ${i}
   next i

jump SELF start_of_loop_group
variable i delete

# Rigid holding settings for the stacking fault group
group stack_group region stack_region

# Rigid holding settings for the twinning fault group
group twin_group region twin_region

# Compute initial energy and output
compute peratom all pe/atom
compute eatoms all reduce sum c_peratom
thermo 100
thermo_style custom step pe c_eatoms press pxx pyy pzz pe temp
#dump config all atom 1000 dump.Stack

# Perform initial relaxation
reset_timestep 0
fix 1 all box/relax x 0.0 y 0.0 z 0.0 vmax 0.01
min_style cg
minimize 1e-25 1e-25 10000 10000

# Twin top N_Twin_Layers
variable twin_move equal -(1.0*${latparam}/sqrt(6))
variable i loop 29
label start_of_loop_twin

   variable k equal ${i}*${twin_move}
   displace_atoms ${i} move ${k} 0.0 0.0 units box
   next i

jump SELF start_of_loop_twin
variable i delete
variable k delete
minimize 1e-25 1e-25 100000 100000
variable E equal "c_eatoms"
variable Eini equal ${E} 
# Create the output file and write header
variable outfile string "./output/stack.dat"
variable nx_points equal 50
variable ny_points equal 50
variable n_incr_x equal ${nx_points}-1
variable n_incr_y equal ${ny_points}-1

variable inc_x equal (-1.0*${latparam}*sqrt(6)/2)/${n_incr_x}
variable inc_y equal (1.0*${latparam}*sqrt(2)/2)/${n_incr_y}

print "Header: Printing Gamma Surface" file ${outfile} screen no

#dump config all atom 1000 dump.Stack
fix 2 all setforce 0 0 NULL

# Outer loop: Apply the y - displacement
variable j loop ${ny_points}
label start_of_loop_y

   # Apply y-displacement from 2nd step
   if "$j > 1" then &
      "displace_atoms stack_group move 0.0 ${inc_y} 0.0 units box"

   # Inner loop: Apply the x- displacement
   variable i loop ${nx_points}
   label start_of_loop_x

      # Apply x-displacement from 2nd step
      if "$i > 1" then &
         "displace_atoms stack_group move ${inc_x} 0.0 0.0 units box" &

      # Relax in z direction
      minimize 1e-25 1e-25 10000 10000

      variable Ecur equal "c_eatoms"
      variable SFED equal (${Ecur}-${Eini})/${Area}
      variable totdisp_x equal (${i}-1)/${n_incr_x}
      variable totdisp_y equal (${j}-1)/${n_incr_y}
      print "${totdisp_x} ${totdisp_y} ${SFED}" append ${outfile} screen no

      next i

   jump SELF start_of_loop_x
   next j

jump SELF start_of_loop_y
variable i delete
variable j delete

# SIMULATION DONE
print "All done"