# ------------------------ INITIALIZATION ----------------------------
units 		real
dimension	3
boundary	p	p	p
atom_style	full
variable latparam equal 2.85

# ----------------------- ATOM DEFINITION ----------------------------
lattice		bcc ${latparam}
region		whole block 0 5 0 5 0 5
create_box	1 whole
lattice 	bcc ${latparam} orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
create_atoms	1 region whole

mass 1 55.85

# ------------------------ FORCE FIELDS ------------------------------
pair_style	reax/c NULL checkqeq yes
pair_coeff	* * ffield.reax.Fe_O_C_H Fe
fix 1 all qeq/reax 1 0.0 10.0 1.0e-6 reax/c

# ------------------------- SETTINGS ---------------------------------
compute csym all centro/atom bcc
compute peratom all pe/atom 

#from stress script
compute stress all stress/atom NULL virial
compute stress1 all reduce ave c_stress[1]
compute stress2 all reduce sum c_stress[2]
compute stress3 all reduce sum c_stress[3]
compute stress4 all reduce sum c_stress[4]
compute stress5 all reduce sum c_stress[5]
compute stress6 all reduce sum c_stress[6]

# Minimize first
#reset_timestep 0
#thermo		10
#thermo_style custom step lx ly lz press pxx pyy pzz pe c_stress1 c_stress2 c_stress3 c_stress4 c_stress5 c_stress6
#min_style cg
#fix 2 all box/relax x 0.0 z 0.0 couple none vmax 0.001 
#minimize 1.0e-25 1.0e-25 1000 10000

# unit fix 
#minimize 1.0e-25 23e-25 1000 10000

#unfix 2

####################################### EQUILIBRATION

#set timestep CHANGED TO .25
reset_timestep	0
timestep 0.25

#set temperature 
velocity all create 300 12345 mom yes rot no
fix 2 all npt temp 300 300 1 iso 0 0 1 drag 1 

# Set thermo output
thermo 1000
thermo_style custom step lx ly lz press pxx pyy pzz pe temp 

# Run for at least 10 picosecond (assuming 1 fs timestep)
run 5000
unfix 2

# Store final cell length for strain calculations
variable tmp equal "lx"
variable L0 equal ${tmp}
print "Initial Length, L0: ${L0}"


####################################### DEFORMATION CHANGED TO 1e15 to match rate
reset_timestep	0

fix		2 all npt temp 300 300 1 y 0 0 1 z 0 0 1 drag 1
variable srate equal 1.0e10
variable srate1 equal "v_srate / 1.0e15"

fix		3 all deform 1 x erate ${srate1} units box remap x

# Output strain and stress info to file
# for units metal, pressure is in [atom] = 1/9869 [GPa]
# p2, p3, p4 are in GPa
variable strain equal "(lx - v_L0)/v_L0"
variable p1 equal "v_strain"
variable p2 equal "-pxx/9869"
variable p3 equal "-pyy/9869"
variable p4 equal "-pzz/9869"
variable p5 equal "c_stress1"
variable p6 equal "c_stress2"
variable p7 equal "c_stress3"
variable p8 equal "c_stress4"
variable p9 equal "c_stress5"
variable p10 equal "c_stress6"
fix def1 all print 100 "${p1} ${p2} ${p3} ${p4} ${p5} ${p6} ${p7} ${p8} ${p9} ${p10}" file Al_SC_100.def1.txt screen no

# Use cfg for AtomEye
dump 		1 all cfg 250 dump.tensile_*.cfg mass type xs ys zs c_csym c_peratom fx fy fz
dump_modify 1 element Fe


# Display thermo CHANGED TO RUN 8000
thermo 	1000
thermo_style	custom step v_strain temp v_p2 v_p3 v_p4 ke pe press 
run		8000

######################################
# SIMULATION DONE
print "All done"