# Define looping variables
variable loopcount loop 42
variable latticeconst index 22.950000 21.412265 20.459239 19.766814 19.222652 18.774324 18.393076 18.061423 17.767934 17.504726 17.266129 17.047934 16.846923 16.660587 16.486928 16.324330 16.171469 16.027243 15.890730 15.761145 15.637820 15.520177 15.407716 15.300000 15.226669 15.149492 15.068041 14.981817 14.890225 14.792550 14.687929 14.575297 14.453324 14.320319 14.174087 14.011705 13.829155 13.620710 13.377787 13.086676 12.723410 12.240000 

# 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 15.9994

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

# 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