File:Drude Model.gif

(*WARNING: Highly unoptimized code!!!*)

Clear[evo, pos, posold, posoldold, accelerations, oldaccelerations, velocities, halfvelocities, tableForces];
Lx = 10; Ly = 3; (*(Half) Size of the box*)
\[Sigma] = 0.5; (*particle radius*) \[Sigma]plot = 0.55 \[Sigma];(*particle radius for graphing*)
dt = 0.002;
nelectrons = 300;
LJforce[p1_, p2_] := 10*2 (12 \[Sigma]^12/Max[Norm[p1 - p2], \[Sigma]]^13(*-6 \[Sigma]^6/Max[Norm[p1-p2],\[Sigma]]^7*)) Normalize[p1 - p2];
externalForce = {10, 0};
poselectrons = RandomVariate[ UniformDistribution[{{-Lx - 100, -Lx + 1}, {-Ly, Ly}}], nelectrons];
posscatterers = Select[Flatten[ Table[N[i {1, 0} + j {Cos[\[Pi]/3], Sin[\[Pi]/3]}], {i, -Lx + 2, Lx + 2, Lx/6}, {j, -Ly, Ly, Ly/2}], 1], (#[[1]] < Lx + 1 && #[[1]] > -Lx + 2) &];
 nscatterers = Dimensions[posscatterers][[1]];
posscatterers = posscatterers + RandomReal[{-\[Sigma], \[Sigma]}/3, {nscatterers, 2}];
initialvelocities = Table[{5, 0}, {j, 1, nelectrons}];
tableForces = Chop@Table[ LJforce[poselectrons[[i]], posscatterers[[j]]], {i, 1, nelectrons}, {j, 1, nscatterers}]; (*electrons do not interact with each other*)
accelerations = Total@Transpose@tableForces + Threaded@externalForce;
posold = poselectrons;
oldaccelerations = accelerations;
poselectrons = posold + initialvelocities*dt + accelerations/2 dt^2;
velocities = initialvelocities + (accelerations + oldaccelerations)/2 dt;

nsteps = 2000;
evo1 = Reap[Do[
      tableForces = Chop@Table[ LJforce[poselectrons[[i]], posscatterers[[j]]], {i, 1, nelectrons}, {j, 1, nscatterers}]; 
      oldaccelerations = accelerations;
      posold = poselectrons;
      halfvelocities = velocities + accelerations/2 dt;
      
      For[j = 1, j <= nelectrons, j++,
       If[poselectrons[[j, 2]] >= Ly, 
        halfvelocities[[j, 2]] = -Abs[halfvelocities[[j, 2]] ] ];
       If[poselectrons[[j, 2]] <= -Ly, 
        halfvelocities[[j, 2]] = Abs[halfvelocities[[j, 2]] ] ];
       ];
      
      poselectrons = posold + halfvelocities dt;
      accelerations = Total@Transpose@tableForces + Threaded@externalForce;
      velocities = halfvelocities + (accelerations + oldaccelerations)/2 dt;
      Sow[poselectrons];
      , nsteps];][[2, 1]];

current = Reverse@GaussianFilter[ Dimensions[Select[#, (#[[1]] > Lx + 1 && #[[1]] < Lx + 3) &]][[1]] & /@ evo1, 10];
start = 800; length = 700;
frames = Table[
   Grid[{{Graphics[{
        Gray, Disk[#, \[Sigma]plot] & /@ posscatterers,
        Black, Disk[#, \[Sigma]plot] & /@ evo1[[j]]
        ,
        Line[{{-Lx - 5, -Ly - \[Sigma]plot}, {Lx + 5, -Ly - \[Sigma]plot}}], 
        Line[{{-Lx - 5, Ly + \[Sigma]plot}, {Lx + 5, Ly + \[Sigma]plot}}]
        }, PlotRange -> {{-Lx - 1, Lx + 1}, {-Ly - 1, Ly + 1}}]
      ,
      ListPlot[current, Joined -> True , PlotRange -> {{0, length} + {nsteps - j, nsteps - j}, {0, 15}}, AspectRatio -> 1/3, Ticks -> None, Epilog -> {Text["Current", Scaled[{0.5, 0.8}] ]} ]
      }}]
   , {j, start, start + length}];
ListAnimate[frames]