# ---------- Initialize Simulation --------------------- 
units metal 
dimension 3 
boundary p p p 
atom_style atomic 

# ---------- Create Atoms --------------------- 
lattice     bcc 1.5 # giving wrong lattice constant on purpose. Should be ~2.85
region  box block 0 1 0 1 0 1 units lattice
create_box  1 box
create_atoms 1 box
replicate 4 4 4 

# ---------- Define Interatomic Potential --------------------- 
pair_style eam/alloy 
pair_coeff * * ./Fe_mm.eam.fs Fe
neighbor 2.0 bin 
neigh_modify delay 10 check yes 
 
# ---------- Run Minimization --------------------- 
fix 1 all box/relax iso 0.0 vmax 0.001
thermo 10 
thermo_style custom step pe lx ly lz press pxx pyy pzz pe
dump 1 all atom 100 dump.lammpstrj
min_style cg 
minimize 1e-25 1e-25 5000 10000 

variable natoms equal "count(all)" 
variable teng equal "pe"
variable length equal "lx"
variable ecoh equal "v_teng/v_natoms"

print ""
print "-------------------------------"
print "Total energy (eV) = ${teng};"
print "Number of atoms = ${natoms};"
print "Lattice constant (Angstoms) = ${length};"
print "Cohesive energy (eV) = ${ecoh};"

print "All done!" 
print "-------------------------------"
print ""

# ---------- Run Dynamics --------------------- 
unfix 1
velocity all create 300.0 87979
fix  2 all nve
run  20000