# ---------- Initialize Simulation --------------------- 
clear 
units real
dimension 3 
boundary p p p 
#atom_style atomic 
#atom_modify map array
atom_style charge

# ---------- Create Atoms --------------------- 
lattice         fcc 4.080
region  box block 0 1 0 1 0 1 units lattice
create_box      1 box

lattice fcc 4.080 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1  

create_atoms 1 box
replicate 1 1 1

mass 1 196.966
set type 1 charge 0.0
# ---------- Define Interatomic Potential --------------------- 
pair_style reax/c NULL
pair_coeff * * ffield.reax.AuO Au
#neighbor 2.0 bin 
#neigh_modify delay 10 check yes 
 
# ---------- Define Settings --------------------- 
compute eng all pe/atom 
compute eatoms all reduce sum c_eng 

# ---------- Run Minimization --------------------- 
#reset_timestep 0 
fix 1 all qeq/reax 1 0.0 10.0 1e-6 reax/c
run 0

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

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

print "All done!"