units		    metal
dimension       3
atom_style	    atomic
boundary        p p p
timestep        0.001
neighbor	    0.3 bin
neigh_modify    every 1 delay 0 check yes


read_data    1.xyz 
# ---------------------- Induce EAM Potentials -------------------

pair_style      eam/alloy
pair_coeff * * FeNiCr_Bonny_2013_ptDef.eam.alloy Fe Ni Cr

minimize		1.0e-12 1.0e-12 10000 10000
minimize		1.0e-12 1.0e-12 10000 10000
minimize		1.0e-12 1.0e-12 10000 10000
minimize		1.0e-12 1.0e-12 10000 10000
minimize		1.0e-12 1.0e-12 10000 10000
minimize		1.0e-12 1.0e-12 10000 10000



########## relaxation ####################################
variable	l1 equal lx
variable	l2 equal ly
variable	l3 equal lz
variable	st equal step
variable	tp equal temp

reset_timestep	0
timestep	0.001

velocity	all create 100.0 6666 rot yes dist gaussian
fix			relax all npt temp 300.0 300.0 0.1 x 0 0 10 y 0 0 10 z 0 0 10
run			5000

unfix		relax
write_data      3.xyz
reset_timestep	0


dump             2  all custom 10000 dump.xyz.*  id type x y z
# ---------- PKA          -----------------
label	setloop
variable	j loop 200

label           shift   
# ---------- random shift  -----------------

variable	ran_x equal random(0,20,1234)        
variable	ran_y equal random(0,20,5678)
variable	ran_z equal random(0,20,9011)

variable	disp_x equal ${ran_x}
variable	disp_y equal ${ran_y}
variable	disp_z equal ${ran_z}

displace_atoms all move ${disp_x} ${disp_y} ${disp_z}  units box
 
variable        xx equal 0.5*${l1}
variable        yy equal 0.5*${l2}
variable        zz equal 0.5*${l3}

variable        region_radius equal  1.5   
region          pkaa sphere ${xx} ${yy} ${zz} ${region_radius} units box  


group          pka region pkaa
variable	    num equal count(pka)

if "${num} > 1" then "jump SELF selectfail"

jump            SELF select     

label           selectfail              
region	        pkaa delete
group           pka delete

jump            SELF shift        



label           select

variable        cos equal      10
variable	    xdel equal ${l1}-10
variable	    ydel equal ${l2}-10
variable	    zdel equal ${l3}-10
region         bound  block 10 ${xdel} 10 ${ydel} 10 ${zdel}  units box 
group           bound  region bound
group          bound1  subtract all bound
    
fix		    1 all nve
fix                2 bound temp/berendsen 300.0 300.0 $(100.0*dt)    
fix         3 all recenter INIT INIT INIT
# ----------set velocity -----------------

#nve Cascade collision
variable    ee    equal 1000
variable	epka  equal ${ee}*1.60217662*10^(-19)#J
variable	matom equal 58.6934*1.9927*10^(-26)/12#kg   

variable	jj equal random(10000,99999,456789)
variable	alphaa equal random(0,180,${jj})            
variable	betaa equal random(0,360,${jj})

variable	alpha equal ${alphaa}*PI/180 
variable	beta  equal ${betaa}*PI/180


#Spherical coordinate system
variable	vel equal sqrt(2*v_epka/v_matom)*0.01#A/ps   
variable	vxx equal v_vel*sin(v_alpha)*cos(v_beta)
variable	vyy equal v_vel*sin(v_alpha)*sin(v_beta)
variable	vzz equal v_vel*cos(v_alpha)
velocity	    pka set ${vxx} ${vyy} ${vzz} sum yes units box  



# ----------run -----------------
fix   adaptive    all dt/reset  1 1.0E-7  1.0E-3 0.014 units box  



run		50000

unfix	1
unfix	2
unfix   3
unfix             adaptive




region			pkaa delete
group			pka delete
region			bound delete
group			bound delete
group			bound1 delete
next			j
jump          SELF setloop
        
quit