clear
kim init Tersoff_LAMMPS_ByggmastarNagelAlbe_2019_FeO__MO_608695023236_000 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 3.647307395935058*${_u_distance}

variable xmin_converted equal 0.0*${_u_distance}
variable xmax_converted equal 22.483512784634385*${_u_distance}
variable ymin_converted equal -47.694733068375555*${_u_distance}
variable ymax_converted equal 47.694733068375555*${_u_distance}
variable zmin_converted equal 0.0*${_u_distance}
variable zmax_converted equal 5.158071585475056*${_u_distance}

# Simulation domain
lattice fcc ${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            fcc ${lattice_constant_converted} &
                                       orient x -1 1 -6 &
                                       orient y -3 3 1 &
                                       orient z 1 1 0

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 5.926874518394469 0 0.0 units box

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

create_atoms       2 region lower
group              lower type 2

# Move atoms in lower region
displace_atoms     lower move -5.926874518394469 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 fcc
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_153.4746/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_153.4746/energy.out

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

# Write out dump file
write_dump         all cfg output/dump_153.4746/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_153.4746/success.out