It genrates a big file bonds.reaxff and an empty file species.out and species.del.
Look at the line
fix 5 all reaxff/bonds 10 bonds.reaxff
fix 6 all reaxff/species 1 5 100 species.out element Si H O F
dump 2 all custom 100 dump.reax.incident.xyz id type x y z
run 1000 # 250 fs
Here is the full input scripts.
LAMMPS Input Script for SiO2 and HF Etching Process using ReaxFF
Initialize simulation
units real
dimension 3
boundary p p f
atom_style charge
Define SiO2 lattice and create simulation box
region simbox block 0 210 0 110 0 650 units box
create_box 4 simbox
mass 1 28.085 # Si
mass 2 15.999 # O
mass 3 1.008 # H
mass 4 19.00 # F
Create SiO2 base
lattice custom 5.43 a1 5.43 0.0 0.0 a2 0.0 5.43 0.0 a3 0.0 0.0 5.43 &
basis 0.0 0.0 0.0 basis 0.5 0.5 0.5 basis 0.25 0.25 0.25 basis 0.75 0.75 0.75
region sio2 block 0 210 0 110 0 240 units box
create_atoms 1 region sio2 basis 1 1 basis 2 1 basis 3 2 basis 4 2
Define potential using ReaxFF
pair_style reaxff NULL safezone 1.6 mincap 100
pair_coeff * * ffield.reax Si O H F
fix 1 all qeq/reax 1 0.0 10.0 1e-6 reaxff
Define neighbor settings
neighbor 6.0 bin
neigh_modify every 1 delay 0 check yes
Output initial atomic coordinates for debugging
dump initial_atoms all custom 1 initial_atoms.xyz id type x y z
run 0
undump initial_atoms
Initial minimization with BFGS method
min_style fire
min_modify dmax 0.1 line quadratic
minimize 1.0e-4 1.0e-6 10000 10000
Check minimization result
dump minimized_atoms all custom 1 minimized_atoms.xyz id type x y z
run 0
undump minimized_atoms
Equilibration in NVT ensemble for 100 ps
fix 2 all nvt temp 308.15 308.15 100.0
timestep 0.25
dump 1 all custom 10000 dump.sio2.*.xyz id type x y z
run 10000 # 2.5 ps to check initial relaxation
Remove previous dump to reset timestep
undump 1
Reset timestep and continue simulation
reset_timestep 0
Define variables for HF insertion
variable E_in equal 40.0 # eV, modify for different energies
variable HF_mass equal 20.0 # amu, approximate mass of HF molecule
variable v_in equal sqrt(2*{E_in}*1.60218e-19*1000/({HF_mass}*1.66054e-27))1.0e100.25/1.0e15 # Å/fs
variable rand_x equal random(0,210,12345)
variable rand_y equal random(0,110,54321)
create_atoms 3 single {rand_x} {rand_y} 649 units box
create_atoms 4 single {rand_x} {rand_y} 649 units box
group hf type 3 4
velocity hf set 0.0 0.0 -${v_in} units box
fix 3 all nve
fix 4 all temp/berendsen 308.15 308.15 100.0
thermo_modify lost ignore
fix 5 all reaxff/bonds 10 bonds.reaxff
fix 6 all reaxff/species 1 5 100 species.out element Si H O F
dump 2 all custom 100 dump.reax.incident.xyz id type x y z
run 1000 # 250 fs
label repeat_start
variable j loop 2
Insert HF molecules every 250 fs for 50 ps
label loop_start
variable i loop 199
variable rand_x2 equal random(0,210,12345)
variable rand_y2 equal random(0,110,54321)
create_atoms 3 single {rand_x2} {rand_y2} 649 units box
create_atoms 4 single {rand_x2} {rand_y2} 649 units box
group hf2 type 3 4
velocity hf2 set 0.0 0.0 -${v_in} units box
thermo_modify lost ignore
dump 3 all custom 1000 dump.reax.incident2-200.xyz id type x y z
run 1000 # 250 fs
undump 3
next i
jump SELF loop_start
Continue simulation for additional 50 ps after HF insertion
run 200000 # 50 ps
Remove gas molecules H2, H2O, O2, and SiFx
fix 7 all reaxff/species 1 5 100 species.out delete species.del specieslist 6 H2 H2O O2 SiF4 SiF5 SiF6
Repeat the entire process two more times
next j
jump SELF repeat_start
Final output settings
unfix 2
unfix 3
unfix 4
unfix 5
Analyze bond formations and output relevant information
compute 2 all pair/local dist
compute 3 all property/atom type
write_data final.data