# ------------------------ INITIALIZATION ----------------------------
units 		metal
dimension	3
boundary	p p p
atom_style	atomic
variable latparam equal 2.856

# ----------------------- ATOM DEFINITION ----------------------------
lattice		bcc ${latparam}
region		whole block 0 200 0 10 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

# ------------------------ regions ------------------------------

region root block 0 20 0 10 0 5
group root region root
region middle block 20 180 0 10 0 5
group middle region middle
region tip block 180 204 0 10 0 5
group tip region tip
# ------------------------ FORCE FIELDS ------------------------------

pair_style	eam/fs
pair_coeff	* * Fe_mm.eam.fs Fe

# ------------------------- pressure calculations ---------------------------------
compute csym all centro/atom bcc
#compute peratomp all pe/atom 

variable volroot index 23295.64
variable volmiddle index 186365.51
variable voltip index 23295.64

#variable volroot equal (.1*lx*ly*lz*2.856*2.856*2.856)
#variable volmiddle equal (.8*lx*ly*lz*2.856*2.856*2.856)
#variable voltip equal (.1*lx*ly*lz*2.856*2.856*2.856)


compute peratom all stress/atom NULL
compute p all reduce sum c_peratom[1] c_peratom[2] c_peratom[3]
variable press equal -(c_p[1]+c_p[2]+c_p[3])/(3*vol)
variable pressx equal -(c_p[1]/vol)

compute peratom1 root stress/atom NULL
compute p1 root reduce sum c_peratom1[1] c_peratom1[2] c_peratom1[3]
variable rootpress equal -(c_p1[1]+c_p1[2]+c_p1[3])/(3*v_volroot)
variable rootpressx equal -(c_p1[1]/v_volroot)

compute peratom2 middle stress/atom NULL
compute p2 middle reduce sum c_peratom2[1] c_peratom2[2] c_peratom2[3]
variable middlepress equal -(c_p2[1]+c_p2[2]+c_p2[3])/(3*v_volmiddle)
variable middlepressx equal -(c_p2[1]/v_volmiddle)

compute peratom3 tip stress/atom NULL
compute p3 tip reduce sum c_peratom3[1] c_peratom3[2] c_peratom3[3]
variable tippress equal -(c_p3[1]+c_p3[2]+c_p3[3])/(3*v_voltip)
variable tippressx equal -(c_p3[1]/v_voltip)

variable nroot equal count(root,root)
variable nmiddle equal count(middle,middle)
variable ntip equal count(tip,tip)

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

#set timestep
reset_timestep	0
timestep 0.001

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

# Set thermo output
thermo 1000
thermo_style custom vol step press temp v_nroot v_nmiddle v_ntip

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

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

variable tmp1 equal "ly"
variable L1 equal ${tmp1}
print "Initial Length, L1: ${L1}"

variable tmp2 equal "lz"
variable L2 equal ${tmp2}
print "Initial Length, L2: ${L2}"

####################################### DEFORMATION
reset_timestep	0

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

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 strainx equal "(lx - v_L0)/v_L0"
variable strainy equal "(ly - v_L1)/v_L1"
variable strainz equal "(lz - v_L2)/v_L2"
variable a1 equal "v_strainx"
variable a2 equal "v_strainy"
variable a3 equal "v_strainz"
variable a4 equal "-pxx/10000"
variable a5 equal "v_pressx/10000"
variable a6 equal "v_rootpressx/10000"
variable a7 equal "v_middlepressx/10000"
variable a8 equal "v_tippressx/10000"

fix def1 all print 100 "${a1} ${a2} ${a3} ${a4} ${a5} ${a6} ${a7} ${a8}" file Fe.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

dump 1 all custom 250 md.dump id type x y z

# Display thermo
thermo 	500
thermo_style	custom step v_nroot v_nmiddle v_ntip press v_press pxx v_pressx v_rootpressx v_middlepressx v_tippressx temp
run		2000

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