# Define looping variables
variable loopcount loop 42
variable latticeconst index 16.331500 15.237232 14.559048 14.066312 13.679080 13.360045 13.088745 12.852738 12.643888 12.456586 12.286798 12.131528 11.988487 11.855888 11.732310 11.616604 11.507826 11.405194 11.308049 11.215836 11.128076 11.044360 10.964332 10.887680 10.835497 10.780576 10.722615 10.661257 10.596078 10.526572 10.452122 10.371971 10.285173 10.190525 10.086464 9.970911 9.841006 9.692673 9.519805 9.312647 9.054141 8.710140 

# Define unit set and class of atomic model
units metal
atom_style atomic

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

# Create a bcc lattice using a single conventional (orthogonal)
# unit cell with a lattice constant from the 'latticeconst' variable defined on line 3 above
lattice bcc ${latticeconst}
region box block 0 1 0 1 0 1 units lattice
create_box 1 box
create_atoms 1 box
mass 1 91.224

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

# Set what thermodynamic information to print to log
thermo_style custom step atoms xlo xhi ylo yhi zlo zhi pe press pxx pyy pzz pxy pxz pyz
thermo 10 # Print every 10 steps

# Set what information to write to dump file
dump id all custom 10 output/lammps.dump id type x y z fx fy fz
dump_modify id 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

# Define auxiliary variables to contain cohesive energy and equilibrium lattice constant
variable poteng    equal "c_thermo_pe"
variable natoms    equal "count(all)"
variable ecohesive equal "v_poteng/v_natoms"

# 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