units metal
newton on on
dimension 3
boundary s p s
echo both

variable ht equal 10
variable avol equal 17.167

timestep 0.0005
thermo 100

atom_style atomic

neighbor 5.0 bin
neigh_modify every 5 delay 10 check yes


restart 1000 restart.A restart.B


# simulation cell definition

variable 	lat equal 3.25

lattice custom ${lat} a1 1 0 0 a2 0 1 0 a3 0 0 1 basis 0.0 0.0 0.0 basis 0.5 0.5 0.5 
region 	box cylinder y 0 0 9 0 50  

create_box	2 box

region		workpiece_lower cylinder y 0 0 9 0 1
region		workpiece_upper cylinder y 0 0 9 49 50

region 		boundw union 2 workpiece_lower workpiece_upper

region		thermw2	cylinder y 0 0 9 48 49
region		thermw3	cylinder y 0 0 9 1 2

region		thermw union 2 thermw2 thermw3

region		newtonw cylinder y 0 0 9 2 48

region 	workpiece union 3 newtonw thermw boundw 

create_atoms 2 region workpiece basis 1 1  basis 2 2 

pair_style	meam
pair_coeff	* * library.meam Zr Cu ZrCu.meam Zr Cu

mass		* 1



# some group definitions
group lower region workpiece_lower

group upper region workpiece_upper

group thermw region thermw

group nvew region newtonw

group fixed region boundw

group moving union nvew thermw

group workpiece region workpiece

# initialize temperature, except for fixed atoms
velocity workpiece create 20 58142275 units box 


# monitor temperature of groups
compute t1 workpiece temp/com

dump  1  all custom  10000 v_ht.*.lammpstrj id type xu yu zu vx vy vz

#minimize 1.0-e4  1.0e-4 1000 100000
fix 		md workpiece nve


variable    kB equal 1.3806504e-23    # [J/K] Boltzmann
variable    ev2J equal 1.602e-19
variable    A2m equal 1.0e-10
variable    ps2s equal 1.0e-12
variable    convert equal ${ev2J}*${ev2J}/${ps2s}/${A2m} 


compute         myKE all ke/atom
compute         myPE all pe/atom
compute         myStress all stress/atom NULL virial


fix tempnve workpiece langevin 0 ${ht} 10 92830

run 5000

unfix  tempnve

undump 1

reset_timestep	0


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


variable srate equal 5.0e8
variable srate1 equal "v_srate / 1.0e12" #put negative for compression
fix		2 all deform 1 y erate ${srate1} units box remap x

# http://lammps.sandia.gov/threads/msg28919.html

# Output strain and stress info to file
# for units metal, pressure is in [bars] = 100 [kPa] = 1/10000 [GPa]
# p2, p3, p4 are in GPa
variable strain equal "(ly - v_L0)/v_L0"
variable p1 equal "v_strain"
variable p2 equal "-pxx/10000"
variable stress equal "-pyy/10000"
variable p4 equal "-pzz/10000"


 compute  tt workpiece stress/atom NULL virial
variable hot1 atom  c_tt[1]/(v_avol*10000)
variable hot2 atom c_tt[2]/(v_avol*10000)
variable hot3 atom c_tt[3]/(v_avol*10000)
variable hot4 atom c_tt[4]/(v_avol*10000)
variable hot5 atom c_tt[5]/(v_avol*10000)
variable hot6 atom c_tt[6]/(v_avol*10000)

variable hot7 atom (0.7071*sqrt((v_hot1-v_hot2)*(v_hot1-v_hot2)+(v_hot2-v_hot3)*(v_hot2-v_hot3)+(v_hot3-v_hot1)*(v_hot3-v_hot1)+6*(v_hot4*v_hot4+v_hot5*v_hot5+v_hot6*v_hot6))) # mechanics of materials

fix k1 workpiece ave/atom 10 100 1000 v_hot1
fix k2 workpiece ave/atom 10 100 1000 v_hot2
fix k3 workpiece ave/atom 10 100 1000 v_hot3
fix k4 workpiece ave/atom 10 100 1000 v_hot4
fix k5 workpiece ave/atom 10 100 1000 v_hot5
fix k6 workpiece ave/atom 10 100 1000 v_hot6
fix k7 workpiece ave/atom 10 100 1000 v_hot7


thermo_style custom step temp press c_t1 v_strain v_stress vol density enthalpy atoms pe ke
dump  1  all custom  5000 v_ht.tensile.*.lammpstrj id type xu yu zu f_k1 f_k2 f_k3 f_k4 f_k5 f_k6 f_k7

run 800000