# The units we assume we'll use. The kim-lammps-preprocessor
# may swap this line out if running against a Simulator Model
# whose units are not 'metal'
units metal

# Variables that can be adjusted by kim-lammps-preprocessor
# to switch unit sets for Simulator Models
variable _u_distance equal 1.0
variable _u_energy equal 1.0
variable _u_time equal 1.0
variable _u_pressure equal 1.0
variable _u_temperature equal 1.0

dimension 3

# This line may be swapped out by kim-lammps-preprocessor if
# running against a Simulator Model whose atom_style is not
# 'atomic'
atom_style atomic

boundary p p p

# Set neighbor skin
variable neigh_skin equal 2.0*${_u_distance}
neighbor ${neigh_skin} bin

# Atom definition
variable latparam_converted equal 3.8232301473617563*${_u_distance}

# 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. Start with x and y dimensions.
# Before converting to the active unit set, get their value in Angstroms
# for later normalization
variable xdim     equal (3.8232301473617563*sqrt(6)/2)*${nlat_x}
variable ydim     equal (3.8232301473617563*sqrt(2)/2)*${nlat_y}

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

# Now convert xdim and ydim to the active unit set and do z direction
variable xdim     equal ${xdim}*${_u_distance}
variable ydim     equal ${ydim}*${_u_distance}
variable zdim     equal (${latparam_converted}/sqrt(3))*${nlat_z}
variable layer_z  equal (${latparam_converted}/sqrt(3))

# Geometry definition
lattice fcc ${latparam_converted}

# Define regions for stacking fault test
# NOTE: -0.001 in the x- and y- directions only to overcome numerical
#       precision issues (no need to worry about units on these quantities)

# 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.001 &
                          0-0.001  ${ydim}+0.001 &
                          0-0.001  ${zdim}+0.001  units box
create_box   1     whole
lattice fcc ${latparam_converted} 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 * * Pd
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

# Variables used to rescale the energy and stress so that the quantities
# in the thermo output are in the original metal units (eV and bars)
# even if we're running with a Simulator Model that uses different units
variable pe_metal    equal "c_thermo_pe/v__u_energy"
variable press_metal equal "c_thermo_press/v__u_pressure"
variable pxx_metal   equal pxx/${_u_pressure}
variable pyy_metal   equal pyy/${_u_pressure}
variable pzz_metal   equal pzz/${_u_pressure}
variable temp_metal  equal temp/${_u_temperature}

# Compute initial energy and output
thermo 100
thermo_style custom step v_pe_metal press &
        v_press_metal v_pxx_metal v_pyy_metal v_pzz_metal &
        temp v_temp_metal

#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_converted}/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 "v_pe_metal"
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_converted}*sqrt(6)/2)/${n_incr_x}
variable inc_y equal (1.0*${latparam_converted}*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 ${pe_metal}
      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"