#LAMMPS input script
units           metal
atom_style      atomic
boundary        p p p
atom_modify     map array

variable        joules2eV equal 1.6E-19
variable        mole equal 6.022E23
variable        g2kg equal 1E-3
variable        dt_anneal equal 0.5e-3

read_data       /ccs/home/mitwood/SetupLAMMPSRuns/DataFiles/W/Pressurize_1000K_1atm_NPT_150Cell.data
mass            1 183.84

compute         lavg all ave/sphere/atom cutoff 5.488
compute         pea all pe/atom
compute         kea all pe/atom
comm_modify     cutoff 5.488
variable        tdamp equal ${tswitch}
variable        twidth equal ${tdamp}/14.0
variable        psratio1 atom 1/(1+exp(-(c_lavg[2]-${tdamp})/v_twidth))         #High T surface
variable        psratio2 atom 1-v_psratio1                                      #Low T surface
compute         sumpsratio all reduce sum v_psratio1
compute         maxpsratio1 all reduce max v_psratio1
compute         minpsratio2 all reduce min v_psratio2

pair_style      hybrid/scaled 1.0 zbl 3.927000 4.488000 v_psratio2 pace product v_psratio1 pace product #Dynamic
pair_coeff      1 1 zbl 74 74
pair_coeff      * * pace 1 W_pot_T0.0.yace W
pair_coeff      * * pace 2 W_pot_T0.5.yace W

thermo_style    custom step etotal pe ke temp press pxx pyy pzz pxy pxz pyz lx ly lz xy xz yz vol density spcpu
thermo          2000
#thermo_modify   norm yes
thermo_modify   lost warn flush yes
timestep        0.5e-3
neighbor        1.0 bin
neigh_modify    once no every 10 delay 0 check yes

variable        buff_therm equal lx*0.025
region          pka_region0 block $(xlo + v_buff_therm) $(xhi - v_buff_therm) $(ylo + v_buff_therm) $(yhi - v_buff_therm) $(zlo + v_buff_therm) $(zhi - v_buff_therm)
group           pka_region0 region pka_region0
group           thermo_region subtract all pka_region0
region          pka_region block $(xlo + v_buff_therm*2) $(xhi - v_buff_therm*2) $(ylo + v_buff_therm*2) $(yhi - v_buff_therm*2) $(zlo + v_buff_therm) $(zhi - v_buff_therm)

group           pka_region region pka_region
variable        atomstrike universe 3391458 #THIS NEEDS TO MATCH THE CELL SIZE USED
variable        PKAEnergy equal random(${pkae},${pkae},${seed1})
group           pka_group id ${atomstrike}
# Set Velocity
variable        pka_mass equal mass[${atomstrike}]
variable        mass_si equal ${pka_mass}*${g2kg}/(${mole})
variable        inner equal 2*${PKAEnergy}*(${joules2eV}/${mass_si})
variable        velocity_si equal sqrt(${inner})
variable        pka_velocity equal ${velocity_si}*1E-2
# Pick Angle
# theta: 0-180          phi: 0-360
variable        u equal random(0,1,${seed3})
variable        v equal random(0,1,${seed4})
variable        thetaInner equal 1-2*${v}
variable        thetaRad equal acos(${thetaInner})
variable        theta equal ${thetaRad}*180/PI
variable        phi equal 360*${u}
# PKA velocity components 
variable        v_x equal ${pka_velocity}*sin(${theta})*cos(${phi})
variable        v_y equal ${pka_velocity}*sin(${theta})*sin(${phi})
variable        v_z equal ${pka_velocity}*cos(${theta})
velocity        pka_group set ${v_x} ${v_y} ${v_z} units box

compute         pe all pe/atom
compute         ke all ke/atom
compute         disp all displace/atom

fix             hist all ave/histo 10 1 100 0.0 100000.0 2000 c_lavg[2] file temperature_adaptive.hist mode vector
fix             1 all dt/reset 100 1.0E-9 5.0E-4 0.05 units box
fix             2 pka_region0 nve
fix             3 thermo_region nvt temp 1000.0 1000.0 0.1 #Just Right
fix             4 all print 100 "$(step) $(dt) $(time)" file out.pka_dt title "step dt time"
variable        pkaSpeed equal (vcm(pka_group,x)^2+vcm(pka_group,y)^2+vcm(pka_group,z)^2)^0.5
thermo_style    custom step dt time etotal pe ke temp press pxx pyy pzz pxy pxz pyz v_pkaSpeed c_sumpsratio spcpu
thermo 100

run 0

#dump           1 all custom 500 DuringAdaptiveTime.dump id x y z c_pe c_disp[*] c_lavg[2]
variable current_dt equal dt # get first dt
label loop_cascade # Loop to run cascade + anneal, break on time = total_cascade_time
variable a loop ${loopmax}
        run 100
        print "current dt ${current_dt}"
        print "dt_anneal ${dt_anneal}"
        if "${current_dt} == ${dt_anneal}" then "print 'Running final steps of cascade loop'" "run 100" "jump SELF break_cascade"
        variable current_dt delete
        variable current_dt equal dt
        next a
jump SELF loop_cascade

label break_cascade

unfix           hist
fix             hist all ave/histo 10 1 1000 0.0 100000.0 2000 c_lavg[2] file temperature_anneal.hist mode vector

timestep        0.5e-3

unfix           1
unfix           2
unfix           3

write_data      postpka.data
#undump         1
dump            1 all custom 20000 PostPKA.*.dump id type x y z c_pe c_disp[*] c_lavg[2]

fix             2 all nvt temp 1000.0 1000.0 0.1 
run             60000 
unfix           2