# 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_pressure equal 1.0
variable _u_mass equal 1.0
variable _u_force equal 1.0

# Isolated atom energy for this species in eV (computed in a separate LAMMPS calculation)
variable isolated_atom_energy equal 0

# Define looping variables
variable loopcount loop 42
variable latticeconst index 9.135420 8.523313 8.143954 7.868329 7.651721 7.473261 7.321502 7.189485 7.072660 6.967888 6.872913 6.786058 6.706044 6.631872 6.562746 6.498022 6.437175 6.379765 6.325424 6.273842 6.224752 6.177923 6.133157 6.090280 6.061090 6.030369 5.997947 5.963624 5.927165 5.888285 5.846640 5.801805 5.753253 5.700309 5.642101 5.577463 5.504797 5.421824 5.325126 5.209247 5.064646 4.872220 

# 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

# Periodic boundary conditions along all three dimensions
boundary p p p

neigh_modify one 4000

# Create a diamond lattice using a single conventional (orthogonal) unit
# cell with a lattice constant from the 'latticeconst' variable defined on line 15 above
variable latticeconst_converted equal ${latticeconst}*${_u_distance}
lattice diamond ${latticeconst_converted}
region box block 0 1 0 1 0 1 units lattice
create_box 1 box
create_atoms 1 box
mass 1 1.0 # Mass inconsequential since we're not performing time integration

# Specify which KIM Model to use
pair_style kim EAM_Dynamo_MendelevAckland_2007v3_Zr__MO_004835508849_000
pair_coeff * * Zr

# Variables used to rescale the box parameters, positions and forces so that the
# quantities in the thermo output and dumpfile are in the original metal units
# (Angstroms and eV/Angstrom) even if we're running with a Simulator Model that uses
# different units
variable     pe_metal  equal       "c_thermo_pe/v__u_energy"
variable    xlo_metal  equal              xlo/${_u_distance}
variable    xhi_metal  equal              xhi/${_u_distance}
variable    ylo_metal  equal              ylo/${_u_distance}
variable    yhi_metal  equal              yhi/${_u_distance}
variable    zlo_metal  equal              zlo/${_u_distance}
variable    zhi_metal  equal              zhi/${_u_distance}
variable     xy_metal  equal               xy/${_u_distance}
variable     xz_metal  equal               xz/${_u_distance}
variable     yz_metal  equal               yz/${_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    pxy_metal  equal              pxy/${_u_pressure}
variable    pxz_metal  equal              pxz/${_u_pressure}
variable    pyz_metal  equal              pyz/${_u_pressure}
variable      x_metal   atom                x/${_u_distance}
variable      y_metal   atom                y/${_u_distance}
variable      z_metal   atom                z/${_u_distance}
variable     fx_metal   atom                  fx/${_u_force}
variable     fy_metal   atom                  fy/${_u_force}
variable     fz_metal   atom                  fz/${_u_force}

# Set what thermodynamic information to print to log
thermo_style custom step atoms v_xlo_metal v_xhi_metal v_ylo_metal v_yhi_metal v_zlo_metal v_zhi_metal &
                    v_pe_metal press v_press_metal v_pxx_metal v_pyy_metal v_pzz_metal &
                    v_pxy_metal v_pxz_metal v_pyz_metal
thermo 10 # Print every 10 steps

# Set what information to write to dump file
dump dumpid all custom 10 output/lammps.dump id type v_x_metal v_y_metal v_z_metal &
                                            v_fx_metal v_fy_metal v_fz_metal

dump_modify dumpid format line "%d %d %16.7f %16.7f %16.7f %16.7f %16.7f %16.7f"

# Compute the energy and forces for this lattice spacing
run 0

# Calculate cohesive energy
variable natoms    equal "count(all)"
variable ecohesive equal "v_pe_metal/v_natoms - v_isolated_atom_energy"

# Output cohesive energy and equilibrium lattice constant
print "Cohesive energy: ${ecohesive} eV/atom"

# Queue next loop
clear # Clear existing atoms, variables, and allocated memory
next latticeconst # Increment latticeconst to next value
next loopcount # Increment loopcount to next value
jump SELF # Reload this input script with the new variable values