(**Velocity Verlet**)

Clear[evo, pos, posold, posoldold, accelerations, oldaccelerations, 
  velocities, halfvelocities, tableForces];
L = 10; (*(Half) Size of the box*)
\[Sigma] = 0.8; (*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) Normalize[p1 - p2];
pos = Flatten[Table[N[{i, j}], {i, -L + 1, 0, L/6}, {j, 0, L - 1, L/6}], 1];
nparticles = Dimensions[pos][[1]]
m = Table[1, {nparticles}];
initialvelocities = Table[0, {j, 1, nparticles}]; initialvelocities[[15]] = {1, 10}; initialvelocities[[4]] = {5, 5};
tableForces = Chop@Table[ LJforce[pos[[i]], pos[[j]]], {i, 1, nparticles}, {j, 1, nparticles}];
MatrixForm[tableForces];
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 = 1000;
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]];

dt = -0.005;
pos = Flatten[Table[N[{i, j}], {i, -L + 1, 0, L/6}, {j, 0, L - 1, L/6}], 1];
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;
evo2 = 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]];

frames3 = Table[
   Graphics[{
     Disk[#, \[Sigma]plot] & /@ evo1[[j]], 
     Disk[#, \[Sigma]plot] & /@ (({2.5 L, 0} + #) & /@ evo2[[j]])
     ,
     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}}],
     Line[({2.5 L, 0} + # &) /@ {{-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}}]
     ,
     Text[ Style["Time flows forward", Bold, Black, FontSize -> 15], {0, 1.2 L}],
     Text[ Style["Time flows backward", Bold, Black, FontSize -> 15], {2.5 L, 1.2 L}],
     Text[ Style["Identical initial conditions", Bold, Black, FontSize -> 16], {1.2 L, L + 4.5}]
     }, PlotRange -> {{-L - 1, 3.5 L + 1}, {-L - 1, L + 5}}, 
    ImageSize -> Large]
   , {j, 1, nsteps}];
ListAnimate[frames3]