(* Initialization *)
n=36;
dt=0.004;
ts=1000;
m=1.12246;
dr=0.2;
v0=1.0;
rcut=3.5;

L=Sqrt[n]*m;

(*Initial positions and velocities*)
Do[
Do[
gx[i+IntegerPart[Sqrt[n]]*(j-1)]=(i-1)*m;
gy[i+IntegerPart[Sqrt[n]]*(j-1)]=(j-1)*m,
{i,1,IntegerPart[Sqrt[n]]}],
{j,1,IntegerPart[Sqrt[n]]}];
Do[
x[i,1]=gx[i]+2*(Random[]-0.5)*dr;
y[i,1]=gy[i]+2*(Random[]-0.5)*dr;
vx[i,1]=2*(RandomInteger[]-0.5)*v0;
vy[i,1]=2*(RandomInteger[]-0.5)*v0;
x[i,0]=x[i,1]-vx[i,1]*dt;
y[i,0]=y[i,1]-vy[i,1]*dt,
{i,1,n}];

(*Visualize the generated initial positions*)
pos1=Table[{x[i,1],y[i,1]},{i,n}];
Do[
p1[1,j]={Part[pos1,j,1],Part[pos1,j,2]};
g1=Array[p1,{1,n}],
{j,1,n}];
Manipulate[Graphics[{PointSize[Large],Point[Part[g1]],PlotRange->{-dr,bl}}],{i,1}]

(* Calculation *)
bl=Sqrt[n]*m+dr;
Do[
Do[
ax=0;
ay=0;
u[i,t]=0;
Do[
If[j≠i,
dx[1]=x[i,t]-x[j,t];
dy[1]=y[i,t]-y[j,t];
(*check periodic boundary condition*)
dx[2]=x[i,t]-x[j,t]+bl;
dy[2]=y[i,t]-y[j,t]+bl;
dx[3]=x[i,t]-x[j,t]-bl;
dy[3]=y[i,t]-y[j,t]-bl;
Δx=Min[Abs[dx[1]],Abs[dx[2]],Abs[dx[3]]];
Δy=Min[Abs[dy[1]],Abs[dy[2]],Abs[dy[3]]];
Do[
If[Δx==-dx[k],Δx=-Δx];
If[Δy==-dy[k],Δy=-Δy],
{k,1,3}];
(*Force calculation*)
rij=Sqrt[((Δx)^2)+((Δy)^2)];
(*If r_ij > r_cut, skip j (vice versa)*)
If[rij<=rcut,
ax=ax+(48(rij^(-14))-24(rij^(-8)))*Δx;
ay=ay+(48(rij^(-14))-24(rij^(-8)))*Δy;
u[i,t]=u[i,t]+4((rij^(-12))-(rij^(-6)))]
],
{j,1,n}];
(*Ordinary Verlet Algorithm*)
x[i,t+1]=2x[i,t]-x[i,t-1]+ax*(dt^2);
y[i,t+1]=2y[i,t]-y[i,t-1]+ay*(dt^2);
vx[i,t]=(x[i,t+1]-x[i,t-1])/(2dt);
vy[i,t]=(y[i,t+1]-y[i,t-1])/(2dt);
(*Periodic Boundary Condition*)
If[x[i,t+1]>  bl,x[i,t+1]=x[i,t+1]-bl;x[i,t]=x[i,t]-bl];
If[x[i,t+1]<-dr,x[i,t+1]=x[i,t+1]+bl;x[i,t]=x[i,t]+bl];
If[y[i,t+1]>  bl,y[i,t+1]=y[i,t+1]-bl;y[i,t]=y[i,t]-bl];
If[y[i,t+1]<-dr,y[i,t+1]=y[i,t+1]+bl;y[i,t]=y[i,t]+bl],
{i,1,n}];
step[t+1]=step[t]+dt,
{t,1,ts}];
(*Visualize*)
pos=Table[{x[i,t],y[i,t]},{t,1,ts},{i,1,n}];
Do[
Do[
p[t,j]=Part[pos,t,j];
g=Array[p,{t,j}],
{j,1,n}],
{t,1,ts}];
Manipulate[Graphics[{PointSize[Large],Point[Part[g,t]],PlotRange->{-dr,bl}}],{t,1,ts,1}]

(* Output *)
(*Energy*)
Do[
ke[t]=0;
pe[t]=0;
Do[
ke[t]=ke[t]+0.5*(((vx[i,t])^2)+((vy[i,t])^2));
pe[t]=pe[t]+u[i,t];
etotal[t]=ke[t]+pe[t],
{i,1,n}],
{t,1,ts}]
ListPlot[Table[ke[t],{t,1,ts}],PlotStyle->Blue]
ListPlot[Table[pe[t],{t,1,ts}],PlotStyle->Green]
ListPlot[Table[etotal[t],{t,1,ts}],PlotStyle->Red]
ListPlot[{Table[ke[t],{t,1,ts}],Table[pe[t],{t,1,ts}],Table[etotal[t],{t,1,ts}]},PlotStyle->{Blue,Green,Red}]

(*Calculate Temperature from Equipartition Theorem*)
T1=0.5Mean[Table[((vx[RandomInteger[{1,n}],t])^2)+((vy[RandomInteger[{1,n}],t])^2),{t,1,ts}]]
T2=0.5Mean[Table[((vx[i,RandomInteger[{1,ts}]])^2)+((vy[i,RandomInteger[{1,ts}]])^2),{i,1,n}]]
T3=0.5Mean[Table[Mean[Table[((vx[i,t])^2)+((vy[i,t])^2),{i,1,n}]],{t,1,ts}]]

(*Speed distribution*)
width=0.1;
vmin=0;
vmax=4;
tmin=100;
tmax=ts;
count=BinCounts[Table[Mean[Table[Sqrt[((vx[i,t])^2)+((vy[i,t])^2)],{t,tmin,tmax}]],{i,1,n}],{vmin,vmax,width}]/n;
listcount=Table[{(i-1)*width+vmin,Part[count,i]},{i,1,((vmax-vmin)/width)}];
ListPlot[listcount,PlotRange->All,Filling->Axis]

(*Temperature Obtain from Curve Fitting to Maxwell distribution*)
p[v_]:=(c*(v^2)/T)Exp[-(v^2)/(2T)]
Clear[c,T]
FindFit[listcount,p[v],{c,T},v]
c=c/.%;
T=T/.%%;
Show[Plot[p[v],{v,vmin,vmax}],ListPlot[listcount,PlotRange->All,Filling->Axis],Plot[(c*(v^2)/T3)Exp[-(v^2)/(2T3)],{v,vmin,vmax},PlotStyle->Red],PlotRange->All]

(*Velocity distribution*)
widthvx=0.1;
vxmin=-3;
vxmax=3;
tmin=100;
tmax=ts;
countvx=BinCounts[Table[Mean[Table[vx[i,t],{t,tmin,tmax}]],{i,1,n}],{vxmin,vxmax,widthvx}]/n;
listcountvx=Table[{(i-1)*widthvx+vxmin,Part[countvx,i]},{i,1,(vxmax-vxmin)/widthvx}];
ListPlot[listcountvx,PlotRange->All,Filling->Axis]

(*
(*Temperature Obtain from Curve Fitting to Maxwell distribution*)
px[vx_]:=(cx/(Tx^(1/2)))Exp[-(vx^2)/(2Tx)]
Clear[cx,Tx]
FindFit[listcountvx,px[vx],{cx,Tx},vx]
cx=cx/.%
Tx=Tx/.%%;
Show[Plot[px[vx],{vx,vxmin,vxmax}],ListPlot[listcountvx,Filling->Axis],Plot[(cx/(T3^(1/2)))Exp[-(vx^2)/(2T3)],{vx,vxmin,vxmax},PlotStyle->Red],PlotRange->All]
*)