clear
kim init EAM_MagneticCubic_MendelevHanSrolovitz_2003_Fe__MO_856295952425_002 metal unit_conversion_mode

dimension          3
boundary           p p p

# Adjust input variables for possible units changes with
# Simulator Models
variable lattice_constant_converted equal 2.855324655771256*${_u_distance}

variable xmin_converted equal 0.0*${_u_distance}
variable xmax_converted equal 9.029329371473409*${_u_distance}
variable ymin_converted equal -36.11731748874896*${_u_distance}
variable ymax_converted equal 36.11731748874896*${_u_distance}
variable zmin_converted equal 0.0*${_u_distance}
variable zmax_converted equal 2.855324655771256*${_u_distance}

# Simulation domain
lattice bcc ${lattice_constant_converted}
region             whole block ${xmin_converted} ${xmax_converted} &
                               ${ymin_converted} ${ymax_converted} &
                               ${zmin_converted} ${zmax_converted} &
                               units box
create_box         2 whole

# Upper region
region             upper block INF INF 0.0 ${ymax_converted} INF INF units box

lattice            bcc ${lattice_constant_converted} &
                                       orient x 3 -1 0 &
                                       orient y 1 3 0 &
                                       orient z 0 0 1

create_atoms       1 region upper
group              upper type 1
variable           mass1_converted equal 55.845*${_u_mass}
mass               1 ${mass1_converted}

# Move atoms in upper region
displace_atoms     upper move 0.0 0 0.0 units box

# Lower region
region             lower block INF INF ${ymin_converted} 0.0 INF INF units box
lattice            bcc ${lattice_constant_converted} &
                                orient x 3 1 0 &
                                orient y -1 3 0 &
                                orient z 0 0 1

create_atoms       2 region lower
group              lower type 2

# Move atoms in lower region
displace_atoms     lower move -0.0 0 -0.0 units box
variable           mass2_converted equal 55.845*${_u_mass}
mass               2 ${mass2_converted}

# Interatomic potential and neighbor settings
kim interactions Fe Fe

# Set up neighbor list method
variable neigh_skin_converted equal 2.0*${_u_distance}
neighbor           ${neigh_skin_converted} bin
neigh_modify       delay 10 check yes

# Delete overlapping atoms
variable delete_distance equal 0.25*${lattice_constant_converted}
delete_atoms       overlap ${delete_distance} all all

# Solver
min_style          cg
minimize           1e-6 1e-6 5000 10000
fix                1 all box/relax x 0.0 y 0.0 z 0.0 couple none vmax 0.001
minimize           1e-6 1e-6 1000 10000
unfix              1

# Variables used to rescale the positions and energies so that
# the quantities in the dumpfile are in the original metal units
# (angstrom and eV) even if we're running with a Simulator Model
# that uses different units
variable     pe_metal  equal       "c_thermo_pe/v__u_energy"
variable     lx_metal  equal             lx/${_u_distance}
variable     ly_metal  equal             ly/${_u_distance}
variable     lz_metal  equal             lz/${_u_distance}
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   mass_metal   atom               mass/${_u_mass}

# Compute centrosymmetry, particle energy, and min dist variables
# Note that if the Model does not define particle energy, this quantity will
# be zero.
compute           csym  all  centro/atom bcc
compute     particle_eng     all  pe/atom
compute  particle_engsum     all  reduce sum c_particle_eng
compute          csymsum     all  reduce sum c_csym
compute         distance     all  pair/local dist
compute          mindist     all  reduce min c_distance

# Isolated atom energy in eV (no unit conversion necessary)
variable          isolated_atom_energy equal 0.0

variable           particle_eng_metal atom "c_particle_eng/v__u_energy - v_isolated_atom_energy"
variable           csym_metal atom c_csym/(${_u_distance}*${_u_distance})
variable           csymsum_metal equal c_csymsum/(${_u_distance}*${_u_distance})
variable           mindist_metal equal c_mindist/${_u_distance}

reset_timestep     0
run                0
thermo             0
thermo_style       custom step v_pe_metal c_particle_engsum &
                          v_lx_metal v_ly_metal v_lz_metal &
                          press v_press_metal v_pxx_metal v_pyy_metal v_pzz_metal &
                          v_mindist_metal v_csymsum_metal

run                0

# Spit out number of atoms
variable           num_atoms equal "count(all)"
print              "${num_atoms}" file output/dump_036.8699/numatoms.out

# Spit out total energy relative to sum of isolated energies, converted back to our
# original metal units
variable  adjusted_pe_metal  equal ${pe_metal}-${num_atoms}*${isolated_atom_energy}
print              "${adjusted_pe_metal} eV" file output/dump_036.8699/energy.out

# Spit out distances, converted back to our original metal units
print              "${mindist_metal} Angstroms" file output/dump_036.8699/mindistance.out

# Write out dump file
write_dump         all cfg output/dump_036.8699/dump_post.cfg &
                   mass type xs ys zs id v_mass_metal v_csym_metal v_particle_eng_metal &
                   modify element Fe Fe

print              "This indicates that LAMMPS ran successfully" file output/dump_036.8699/success.out