# cumulative H bombardment on W GB
# By Aws Al-Shalash
# YOU MUST USE lmp_g++ -in thisfilename to run. Do not use < to replace -in.

# Recent Edits: 1- Rearranged the Variables, Compute, and Thermo Style outside the loop eleminating Warning
# 2- Got the log file outside the loop which is recording all the bombardments and not the last one only! 
# 3- Changed defined 2 dump variables one will produce a string and the other a solid file 
# 4- Got the rn_seed and the rest of the variables to the loop because previously the code did not increment the The random seed is ${rn_seed}"
# 5- Added new conditions for timestep using if-every command  
# 6- Added a fix to the bottom of the substrate 
# 7- Defined a new region for H atom & "region instead of type" in group g_H & delete the group by the end of each bombardment
# 8- Added a fix to the bottom of the substrate since the restart file will lose the fix when read
# 9- Optimized the timestep by adding a loop making sure the simulation time is the same for all bombardments & at the same time decreasing the Wall time
# 10- Added a H counter to the sample that will tell the number of H trapped after each bombardment
# 11- Added Berdensten Thermostat 
variable result_file string result_file.txt
shell echo " " >> ${result_file}
shell echo "===============================" >> ${result_file}
shell echo "$(date) $(ls -1 | wc -l)" >> ${result_file}

read_restart Wbase_GB_sigma5_Z460.restart.2000

variable H_idx loop 3000      # control the number of H bombardment

# ---- Geo variable ----
variable maxx equal 50.0
variable maxy equal 50.0
variable bottomz equal -460.0
variable minz equal -440.0
variable midz equal -220.0
variable topz equal  140.0
variable fixed_W_topline equal -420

variable massW equal 183.84
variable massH equal 1.00794

variable berendsen_xmin equal 1
variable berendsen_xmax equal ${maxx}-1
variable berendsen_ymin equal 1
variable berendsen_ymax equal ${maxy}-1
variable berendsen_zmin equal ${fixed_W_topline}+3.1
variable berendsen_zmax equal 1

# Steps & Loops
variable ts equal 0.0005  # timestep for run
variable ts1 equal 0.0009 # timestep for run1
variable L equal 20       # No. of loops for run 
variable run equal 500    # Steps at 0.0005

# variable for H atom #
variable Heng equal 80 # H energy in eV
variable Hvelocity equal (-1.0)*sqrt(${Heng}*1.602e-19*2.0/(${massH}*1.661e-27))*1.e10/1.e12 # do not change this line

variable W_temperature equal 600     # kelvin

########### berendsen group, must define before the rescalling ###########

region rgn_1 block ${berendsen_xmin} ${berendsen_xmax} &
${berendsen_ymin} ${berendsen_ymax} &
${berendsen_zmin} ${berendsen_zmax} units lattice
group grp_1 region rgn_1
group grp_berendsen_W subtract all grp_1 group_fixed_W
group grp_1 delete

region A block 0 50.0 0 50.0 14.9 15.2
fix 1 all nve
#fix 2 group_mobile_W temp/rescale 10 ${W_temperature} ${W_temperature} 0.1 1.0
#unfix 2 # Maybe Doesnt work so remove it
fix 2 grp_berendsen_W temp/berendsen ${W_temperature} ${W_temperature} 0.01
fix zeroforce group_fixed_W setforce 0 0 0 
compute Wtemp group_mobile_W temp 

dump 3 all image 10000 W_${Heng}eV_${H_idx}-*.jpg type type size 1024 1280 zoom 2 box yes 0.005 axes yes 0.5 0.05 view 60 -30
log W_GB_bom_${Heng}eV_cont.log 
################################################################################
label beginning 

variable rn_seed equal ${H_idx}*20+17
variable rx equal random(0,${maxx},${rn_seed})
variable ry equal random(0,${maxy},${rn_seed})
variable rz equal random(15,15.1,${rn_seed})

#variable Wmobileline equal ${minz}+0.5 #What Doesn this do ???

reset_timestep 0

create_atoms 3 single ${rx} ${ry} ${rz} units box
group g_H region A
velocity g_H set 0 0 ${Hvelocity} # perpendicular downward

variable H_x equal xcm(g_H,x)
variable H_y equal xcm(g_H,y)
variable H_z equal xcm(g_H,z)
variable H_ke equal ke(g_H)
variable H_n equal count(g_H)

thermo_style	custom step time temp c_Wtemp v_H_x v_H_y spcpu v_H_ke v_H_z
thermo_modify lost warn
thermo_modify format 2 %4.1f
thermo_modify format 3 %4.1f
thermo_modify format 4 %4.1f 
thermo_modify format 5 %4.1f 
thermo_modify format 6 %4.1f
thermo_modify format 7 %5.2f
thermo_modify format 8 %5.2f
thermo_modify format 9 %4.1f

thermo 20

pair_style tersoff
pair_coeff * * WH6.tersoff W W H
neighbor 2.0 bin 
neigh_modify delay 0 every 1 check yes

print "==============================================================="
print "This is No. ${H_idx} bombardment. The random seed is ${rn_seed}"
print "==============================================================="

timestep ${ts}

variable a loop ${L} 
################################################################################
label run
run ${run}
if "${H_z}>6.0 || ${H_ke}<0.25" then "jump SELF run1" 
print "$(step) is the current step."
if "${a}==${L}" then "print 'Warning: max step reached.'" "jump SELF next"
next a
print "${a} is the current bomb. loop variable."
jump SELF run
################################################################################
label run1

variable run1 equal round((5-($(step)*${ts}))/${ts1}) #Maintaining the simulation time while increasing timestep
print "${run1} the remianing steps @ 0.0009 timestep"

timestep ${ts1}
run ${run1} 
################################################################################
label next 
# prepare the next bombardment

group H type 3 
variable H_NUM equal count(H)  #To count the number of H trapped at current bombardment
print "${H_NUM} is the current H atoms trapped in the W sample @ ${H_idx} ."

dump 1 all atom $(step) Hbom_${Heng}ev_${H_idx}.lammpstrj.gz
dump 2 all image $(step) W_${Heng}eV_${H_idx}-*.jpg type type size 1024 1280 zoom 2 box yes 0.005 axes yes 0.5 0.05 view 60 -30
dump_modify 2 adiam 1 0.5 adiam 2 0.5
run 0
undump 1
undump 2
variable H_NUM delete
group H delete  #delete group before next bombardment!!
variable a delete
group g_H delete  #delete group before next bombardment!!
next H_idx
jump SELF beginning