(**Velocity Verlet**)
Clear[evo, pos, posold, posoldold, accelerations, oldaccelerations, velocities, halfvelocities, tableForces];
L = 10.25; (*(Half) Size of the box*)
\[Sigma] = 0.55; (*particle radius*) \[Sigma]plot = 0.55 \[Sigma];(*particle radius for graphing*)
dt = 0.005;
LJforce[p1_, p2_] :=  4*2 (12 \[Sigma]^12/Max[Norm[p1 - p2], \[Sigma]]^13 - 6 \[Sigma]^6/Max[Norm[p1 - p2], \[Sigma]]^7) Normalize[p1 - p2];
pos1 = Flatten[ Table[N[{i, j} + RandomReal[{-L/30, L/20}, {2}]], {i, -L + \[Sigma]plot, 0, L/6}, {j, -L + \[Sigma]plot, L, L/6}], 1];
pos2 = Flatten[ Table[N[{i, j} + RandomReal[{-L/20, L/20}, {2}]], {i, 0, L - \[Sigma]plot, L/6}, {j, -L + \[Sigma]plot, L, L/6}], 1];
pos = Join[pos1, pos2];
nparticles = Dimensions[pos][[1]]
m = Table[1, {nparticles}];
initialvelocities = RandomReal[{-10, 10}, {nparticles, 2}];
tableForces = Chop@Table[ LJforce[pos[[i]], pos[[j]]], {i, 1, nparticles}, {j, 1, nparticles}];
accelerations = Table[Sum[m[[i]]*tableForces[[i, j]], {j, 1, nparticles}], {i, 1, nparticles}];
posold = pos;
oldaccelerations = accelerations;
pos = posold + initialvelocities*dt + accelerations/2 dt^2;
velocities = initialvelocities + (accelerations + oldaccelerations)/2 dt;
nsteps = 300;
evo1 = Reap[Do[
      tableForces = Chop@Table[ LJforce[pos[[i]], pos[[j]]], {i, 1, nparticles}, {j, 1, nparticles}];
      oldaccelerations = accelerations;
      posold = pos;
      halfvelocities = velocities + accelerations/2 dt;
      
      For[j = 1, j <= nparticles, j++,
       If[pos[[j, 2]] >= L, halfvelocities[[j, 2]] = -Abs[halfvelocities[[j, 2]] ] ];
       If[pos[[j, 2]] <= -L, halfvelocities[[j, 2]] = Abs[halfvelocities[[j, 2]] ] ];
       If[pos[[j, 1]] >= L, halfvelocities[[j, 1]] = -Abs[halfvelocities[[j, 1]] ] ];
       If[pos[[j, 1]] <= -L, halfvelocities[[j, 1]] = Abs[halfvelocities[[j, 1]] ] ];
       ];
      
      pos = posold + halfvelocities dt;
      accelerations = Chop@Table[ Sum[m[[i]]*tableForces[[i, j]], {j, 1, nparticles}], {i, 1, nparticles}];
      velocities = halfvelocities + (accelerations + oldaccelerations)/2 dt;
      Sow[pos];
      , nsteps];][[2, 1]];
frames1 = Table[Graphics[{
     Disk[#, \[Sigma]plot] & /@ evo1[[j, 1 ;; nparticles/2]]
     ,
     Line[{{-L - \[Sigma]plot, -L - \[Sigma]plot}, {-L - \[Sigma]plot, L + \[Sigma]plot}, {L + \[Sigma]plot, L + \[Sigma]plot}, {L + \[Sigma]plot, -L - \[Sigma]plot}, {-L - \[Sigma]plot, -L - \[Sigma]plot}}]
     ,
     Orange, 
     Disk[#, \[Sigma]plot] & /@ evo1[[j, nparticles/2 + 1 ;; -1]]
     }, PlotRange -> {{-L - 1, L + 1}, {-L - 1, L + 1}}]
   , {j, 1, nsteps}];
ListAnimate[frames1]