Remove overlapping using fix nve/limit

Dear all ,
I am trying to prepare LJ glass sample to apply athermal quasistatic shear deformation. I am creating 10000 particles in a 2D simulation box with density 0.976. While doing so the particles are getting highly overlapped which is generating huge forces and velocities over the particles , leading to losing atoms while running the normal dynamics.
As mentioned in the lammps manual , to overcome this issue it is better to run initial dynamics with fix nve/limit.
I am trying to implement the same but still I am facing the issue of losing particles whenever I am trying to start the equilibration of the sample at high temperature. Please help me understand if I am doing anything wrong in using fix nve/limit.
I am using lammps version of July, 2021. I have attached the code I am using and the paper I am referring for the parameters. The reference paper is using a modified version of lj/cut , I have made necessary edit in the .cpp file of the same.

units           lj   # lj or real or metal or si or cgs or electron or micro or nano (unitless
atom_style      sphere
boundary p p p
dimension       2

region      box prism 0.0 101.0 0.0 101.0 -0.5 0.5 0.0 0.0 0.0


create_box      2 box
create_atoms 1 random 5000 3000 NULL  units box
create_atoms 2 random 5000 4000 NULL  units box

set type 1 mass 1.0
set type 2 mass 1.0
set type 1 diameter 1.175
set type 2 diameter 0.618

comm_modify vel yes
neighbor  0.55 bin
pair_style lj/cut 2.5
pair_coeff 1 1 0.5 1.175
pair_coeff 2 2 0.5 0.618
pair_coeff 1 2 1.0 1.0

neigh_modify every 1 delay 0 check yes
timestep 0.0001
fix 100 all enforce2d
fix 2000 all nve/limit 0.1
thermo 100000
thermo_style custom step temp ke pe
dump    mydumpp all atom 100000 initial.dump
dump             2 all custom 100000 *.data id type diameter mass x y z vx vy fx fy
run 10000000

Any suggestion is appreciated
https:journals.aps.org:prl:pdf:10.1103:PhysRevLett.109.255502.pdf (395.0 KB)
https:journals.aps.org:prl:supplemental:10.1103:PhysRevLett.109.255502:EshelbyAppendix.pdf (291.0 KB)

Why such an old version? You are missing out on lots of useful features and bugfixes.

For example, since version 2 June 2022 there is an overlap keyword that can be used to avoid the closest contacts and thus highest potential energy for pairs of atoms caused by random placement.

Why not run a minimization first?

Dear sir
I have a question, as you suggested to run a minimisation first , that means I have to run a minimisation first then use fix nve/limit as an initegrator during the equilibration at high temperature. This is the part I am not able to understand, I will be greatful if you can help me under stand this part

My first and most important suggestion was to use a more recent LAMMPS version and then the “overlap” keyword of the create_atoms command. That will reduce the problem before it manifests itself.

You have to explain, what your concerns are. What I am suggesting is a common protocol for setting up MD simulations or all kinds. So if you want more help with that, you have to provide the information about what you are struggling with to understand. For me the whole thing is just “obvious”.