I am new with LAMMPS. Kindly check below my input script and what might be th issue
===== 5 keV H → Si(10x10x260)irradiation with 100 H atoms =====
---- Initialization ----
units metal
dimension 3
boundary p p f
atom_style atomic
---- Lattice / Box Setup ----
variable a equal 5.431 # Å (Si lattice constant)
variable nx equal 10
variable ny equal 10
variable nz equal 260 # 130 cells = 706 Å depth
variable nzEXT equal ${nz}+20 # headroom above surface (cells) for injection
lattice diamond {a}
region box block 0 {nx} 0 {ny} 0 {nzEXT} units lattice
region si block 0 {nx} 0 {ny} 0 ${nz} units lattice
create_box 2 box # 1=Si, 2=H
create_atoms 1 region si
mass 1 28.0855 # Si
mass 2 1.008 # H
---- Geometry (Å) ----
variable Lx equal {nx}*{a}
variable Ly equal {ny}*{a}
variable Lz equal {nz}*{a} # Surface is at z = Lz
---- Neighbor / comm ----
neighbor 2.0 bin
neigh_modify delay 0 every 1 check yes one 20000 page 400000
comm_modify cutoff 40.0
---- Groups / Regions ----
Bottom thermostat slab
variable zbath equal 3.0 # ~3 lattice planes for fix
region bath block 0 {nx} 0 {ny} 0 ${zbath} units lattice
group Gbath region bath
group Gsi type 1 # Group for Silicon atoms
group H type 2 # Group for all Hydrogen atoms (Automatically updated)
group mobile subtract all Gbath # Mobile atoms (all atoms except the bottom fixed layer)
Injection region (5–10 Å above free surface)
variable z1 equal {Lz}+5.0
variable z2 equal {Lz}+10.0
region deposit block 0 {Lx} 0 {Ly} {z1} {z2} units box
---- Potentials ----
NOTE: Ensure Si.tersoff file is in the working directory
pair_style hybrid/overlay tersoff zbl 1.5 2.5
pair_coeff * * tersoff Si.tersoff Si NULL
pair_coeff 1 1 zbl 14 14
pair_coeff 1 2 zbl 14 1
pair_coeff 2 2 zbl 1 1
---- Integration & boundaries ----
velocity Gsi create 300.0 4928459 mom yes rot yes dist gaussian
fix fb Gbath langevin 300.0 300.0 0.1 90421
fix nve1 Gbath nve
fix nve2 mobile nve
Reflect at bottom; open at top (free surface)
fix wallbot all wall/reflect zlo EDGE
---- Timestep control ----
timestep 0.00005 # 0.05 fs
fix adaptdt all dt/reset 1 0.00002 0.0002 0.1 units box
thermo 2000
thermo_style custom step atoms temp pe ke etotal press
thermo_modify lost warn
run 0 # Initial system setup
---- Beam Velocity Calculation (Å/ps) ----
v = 9790 Å/ps for 5keV H.
variable vH equal 9790.0
---- Depth Binning (Chunk Define) ----
variable dz equal 0.2 # Fine 0.2 Å bins for resolution
compute cZ all chunk/atom bin/1d z lower ${dz} units box
---- Hydrogen Implantation Profile (Density) ----
Samples H atoms, bins them by Z, outputs density (count/volume of bin).
Accumulates over the entire simulation (100 shots).
fix Hprof H ave/chunk 10 100 1000 cZ density/number file H_profile.txt
x_modify Hprof ave running
---- Vacancy Profile (Under-coordination CN<4 on Si) ----
Vacancy profiling: under-coordination CN<4 on Si
1. Compute per-atom coordination on Si group
compute cn Gsi coord/atom cutoff 2.85 # Si-Si bond length, adjust if needed (try 3.0 if no vacancies)
2. Define per-atom vacancy variable (1 if CN < 4, 0 otherwise)
variable vdef atom (c_cn < 4)
3. Use ave/chunk on Si group to sum up the '1’s (the vacancies) in each bin.
We use ‘sum’ here to get the total number of vacancies counted per bin.
fix Vprof Gsi ave/chunk 10 100 1000 cZ v_vdef file Si_vac_profile.txt
If no vacancies are detected, try increasing cutoff to 3.0 or higher
---- Dumps ----
dump dtraj all custom 5000 dump.Si_H_5keV.lammpstrj id type x y z
dump_modify dtraj sort id
dump htraj H custom 20 H_trajectory.dump id type x y z vx vy vz # Only track H atoms
---- Bombardment Parameters (Optimized) ----
variable gap equal 1000 # Relaxation steps between shots (increased for better thermalization)
variable impact_time equal 500 # Time for high-energy impact (increased for deeper penetration)
variable n_shots equal 10 # Total H atoms to fire
variable shot_counter loop 1 ${n_shots}
---- Automated H Bombardment Sequence (Loop) ----
label BOMBARDMENT_LOOP
1. Create a new H atom in the deposit region (just above Si surface)
This ensures H atoms are injected at z = Lz+5 to Lz+10
create_atoms 2 random 1 {shot_counter} deposit variable last_atom equal atoms group newH id {last_atom}
Ensure all H atoms are always in group H (type 2)
group H type 2
2. Assign high velocity and small timestep (high-energy impact)
— Randomize impact angle to reduce channeling and transmission —
Reason: Channeling at normal incidence allows H atoms to escape through the slab. Randomizing the angle increases straggling and retention.
variable theta_deg equal random(0,10,12345) # Up to 10 degrees off normal (reduces channeling)
variable phi_deg equal random(0,360,67890) # Random azimuthal angle
variable theta_rad equal {theta_deg}*3.141592653589793/180
variable phi_rad equal {phi_deg}*3.141592653589793/180
variable vx equal {vH}*sin({theta_rad})*cos({phi_rad})
variable vy equal {vH}*sin({theta_rad})*sin({phi_rad})
variable vz equal -{vH}*cos({theta_rad})
velocity newH set {vx} {vy} ${vz} sum no units box
— Assign downward velocity (already set by vz = -vH*cos(theta_rad)) —
— Run high-energy impact and profile H atoms —
timestep 0.00001
Increased impact_time for deeper penetration and realistic straggling
run ${impact_time}
Hprof fix is active throughout simulation, so H atoms will be profiled after each run
3. Relaxation (larger timestep)
timestep 0.00002
Increased gap for better thermalization and retention
run ${gap}
4. Cleanup and iterate
group newH delete # Delete the temporary group definition
Iterate to next shot
next shot_counter
jump Si_5keV.in BOMBARDMENT_LOOP
---- Final Relaxation (Optional) ----
timestep 0.001 # Larger final timestep for long relaxation
run 50000
print “Total H atoms implanted: ${n_shots}”
print “All done.”