Systems of many ordinary differential equations: Molecular dynamics

Many particles with interaction in 2 dimensions

Theory

Statistical distribution at equilibrium

If the gas of particles is rarified, collisions between particles serve to establish thermodynamic equilibrium in which potential energy of particle's interaction does not play a role. At equilibrium, the components of the particles' velocities satisfy the Maxwell distribution

dN=N(Sqrt[(m/(2π(Subscript[k, B])T))]) Exp[-((Subscript[mv, x]^2)/(2(Subscript[k, B])T))]d(Subscript[v, x])

etc. This distribution function is normalized by

∫dN=N

and the average kinetic energy of the motion along the x direction is given by

(1/N)∫((Subscript[mv, x]^2)/2)dN=(1/2)(Subscript[k, B])T.

Indeed,

Clear[T]
Integrate[(Sqrt[(m/(2π(Subscript[k, B])T))]) Exp[-((m (Subscript[v, x]^2))/(2(Subscript[k, B])T))],{(Subscript[v, x]),-∞,∞},GenerateConditions->False]
Integrate[((m (Subscript[v, x]^2))/2)(Sqrt[(m/(2π(Subscript[k, B])T))]) Exp[-((m (Subscript[v, x]^2))/(2(Subscript[k, B])T))],{(Subscript[v, x]),-∞,∞},GenerateConditions->False]

Let us consider now the distribution over kinetic energies for the particles moving in two dimensions. Replacing d(Subscript[v, x])d(Subscript[v, y])=2πvdv, one obtains

dN=N(m/(2π(Subscript[k, B])T)) Exp[-((m((Subscript[v, x]^2)+(Subscript[v, y]^2)))/(2(Subscript[k, B])T))]d(Subscript[v, x])d(Subscript[v, y])=N(m/((Subscript[k, B])T)) Exp[-((mv^2)/(2(Subscript[k, B])T))]vdv

The average kinetic energy reads

EkAvr=Integrate[((m (v^2))/2)(m/((Subscript[k, B])T)) Exp[-((m (v^2))/(2(Subscript[k, B])T))]v,{v,0,∞},GenerateConditions->False]

that is, twice as much as before because we take into account motions along both x and y axes. The quadratic dispersion of the kinetic energy at equilibrium (Subscript[ΔE, k]) is given by

((((Subscript[ΔE, k])))^2)=((∫_0)^∞)((((Subscript[E, k])-(Subscript[(OverBar[E]), k])))^2)(m/((Subscript[k, B])T)) Exp[-((mv^2)/(2(Subscript[k, B])T))]vdv

that is,

ΔEk=√Integrate[(((((m (v^2))/2)-EkAvr))^2)(m/((Subscript[k, B])T)) Exp[-((m (v^2))/(2(Subscript[k, B])T))]v,{v,0,∞},GenerateConditions->False]

the same as the average kinetic energy.

Approaching equilibrium

Starting from any initial state of a gas, collisions between particles slowly make the system relax to the equilibrium described by the Maxwell distribution above. To illustrate the process of approaching the equilibrium, we start with, say, a state in which all particles in a container move to the right with the kinetic energy (Subscript[k, B])T per particle and run a molecular dynamic (MD) simulation solving Newtonean equations of motion for each of N weakly interacting particles in a container described by a confining potential.

After several elastic collisions with the walls of the container, the particles will be equidistributed over directions of their velocities. This process is called chaotization. However, in the absence of collisions, each particle will retain its kinetic energy, except for short moments of hitting the walls of the container. Thus, after chaotization, the particles will have the correct value of the average kinetic energy (Subscript[(OverBar[E]), k])=(Subscript[k, B])T but the energy dispersion (Subscript[ΔE, k]) will be zero. Then, because of rarer collisions, (Subscript[ΔE, k]) will slowly increase to its asymptotic value (Subscript[k, B])T that corresponds to the equilibrium. Thus the energy dispersion (Subscript[ΔE, k]) can be used as an indicator of relaxation.

Parameters and definitions

NPart=50;    (* Number of particles *)
tMax=1000;   (* Maximal time of the calculation *)
m=1;             (* Mass of a particle *)
kB=1;          (* Boltzmann constant *)

(* Interaction between the particles *)
rInt=.01;     (* Set to zero to kill interaction *) 
PowInt=4;      UInt[r_]:=((((rInt/r)))^PowInt);

(* Confinement of the particles *)
rConf=10;    PowConf=50;       UConf[r_]:=((((r/rConf)))^PowConf);

(* Total potential energy *)
U:=(∑_(i=1)^NPart)UConf[(Sqrt[((x[i][t])^2)+((y[i][t])^2)])]+(∑_(i=1)^NPart)(∑_(j=1)^(i-1))UInt[(Sqrt[(((x[i][t]-x[j][t]))^2)+(((y[i][t]-y[j][t]))^2)])]  

(* Initial condition for particle's speed *)
T=1;     (* Temperature of the system *)
v0x=(Sqrt[ ((2kB T)/m)]);      (* The whole thermal energy in the x-motion *) 
v0y=0;    (* and no energy in the y motion *)

tConf=(rConf/v0x);        (* Time to cross the confinement region *)

(* Initial localization for particle's positions *)
aIni=0.5rConf;           (* In the initial state particles are randomly localized within a square of side aIni *)


Plot[{UInt[Abs[r-0.5rConf]],UConf[r]},{r,0,1.1rConf},PlotRange->{0,2}]

Equations of motions and their solution

Timing[
Equations=Flatten[Table[{
m vx[i]'[t]==-(∂_(x[i][t]))U,          x[i]'[t]==vx[i][t],
m vy[i]'[t]==-(∂_(y[i][t]))U,          y[i]'[t]==vy[i][t]    
},{i,1,NPart}]];
]

IniConds=Join[
Flatten[Table[{
x[i][0]==RandomReal[{-aIni,aIni}],
y[i][0]==RandomReal[{-aIni,aIni}]
},{i,1,NPart}]],
Flatten[Table[{vx[i][0]==v0x,vy[i][0]==v0y},{i,1,NPart}]]
];

Vars=Flatten[Table[{x[i],vx[i],y[i],vy[i]},{i,1,NPart}]];

Timing[
Sol=NDSolve[Join[Equations,IniConds],Vars,{t,0,tMax},MaxSteps->10000000];
]

xt[i_][t_]:=x[i][t]/.Sol[[1]];            vxt[i_][t_]:=vx[i][t]/.Sol[[1]];
yt[i_][t_]:=y[i][t]/.Sol[[1]];            vyt[i_][t_]:=vy[i][t]/.Sol[[1]];

Plotting results

Without interaction

Trajectory of a particle in the absence of interactions, rInt==0, is regular.

ii=1;
Show[
ParametricPlot[{xt[ii][t],yt[ii][t]},{t,0,tMax},PlotRange->{-1.01rConf,1.01rConf},PlotPoints->300],
Graphics[{
{Thick,Circle[{0,0},rConf]}
}]
]

Kinetic energy of a particle has many narrow dips due to collisions with the walls, but on average kinetic energy is conserved

Plot[((m )/2)(((vxt[ii][t])^2)+((vyt[ii][t])^2)),{t,0,tMax},PlotPoints->300,PlotRange->All]
Plot[((m )/2)(((vxt[ii][t])^2)+((vyt[ii][t])^2)),{t,0,0.02tMax},PlotRange->All]

With interaction

Trajectory of a particle in the presence of interactions, rInt > 0, shows collisions with other particles and looks random, although it is deterministic.

ii=1;
Show[
ParametricPlot[{xt[ii][t],yt[ii][t]},{t,0,tMax},PlotRange->{-1.01rConf,1.01rConf},PlotPoints->300],
Graphics[{
{Thick,Circle[{0,0},rConf]}
}]
]

Kinetic energy of a particle has many narrow dips due to collisions with the walls and has jumps because or rare collisions with other particles

Plot[((m )/2)(((vxt[ii][t])^2)+((vyt[ii][t])^2)),{t,0,tMax},PlotPoints->300]

Analysis of the energy distribution

Definitions

Ekx[i_][t_]:=((m )/2)((vxt[i][t])^2)   (* Kinetic energy of the x-motion of particle i *)
Eky[i_][t_]:=((m )/2)((vyt[i][t])^2)   (* Kinetic energy of the y-motion of particle i *)
Ek[i_][t_]:=Ekx[i][t]+Eky[i][t]

EkxAvr[t_]:=(1/NPart)(∑_(i=1)^NPart)Ekx[i][t]      (* Average kinetic energies of the particles *)
EkyAvr[t_]:=(1/NPart)(∑_(i=1)^NPart)Eky[i][t]  
EkAvr[t_]:=(1/NPart)(∑_(i=1)^NPart)Ek[i][t] 
ΔEk[t_]:=(Sqrt[(1/NPart)(∑_(i=1)^NPart)(((Ek[i][t]-EkAvr[t]))^2)]);    (* Energy dispersion *)
ΔEk[t_]:=(Sqrt[(1/NPart)(∑_(i=1)^NPart)((Ek[i][t])^2)-((EkAvr[t])^2)]);    (* Energy dispersion rearranged *)

Without interaction

Direction of motion of one particle changes because of collision with the container walls and the energies of motion along x and y are converted into each other.

ii=1;
Plot[{Ekx[ii][t],Eky[ii][t]},{t,0,0.2tMax},PlotRange->All,PlotStyle->{{Blue,Thick},{Red,Thick}}]

However, the total kinetic energy is conserved, except for the dips due to collision with the container walls

ii=1;
Plot[Ek[ii][t],{t,0,0.1tMax},PlotRange->All,PlotStyle->{Black,Thick}]

After averaging over N particles one can see the chaotization that is not complete because the number of particles is still too low.

Plot[{EkxAvr[t],EkyAvr[t],EkAvr[t]},{t,0,tMax},PlotRange->All,PlotStyle->{{Blue,Thick},{Red,Thick},{Black,Thick}}]

Plot[ΔEk[t],{t,0,tMax},PlotRange->All,PlotStyle->{Black,Thick},PlotPoints->300]

With interaction

Direction of motion of one particle changes because of collision with the container walls and the energies of motion along x and y are converted into each other.

ii=1;
Plot[{Ekx[ii][t],Eky[ii][t]},{t,0,tMax},PlotRange->All,PlotStyle->{{Blue,Thick},{Red,Thick}},PlotPoints->300]

The total kinetic energy of a particle is not conserved because of collisions between the particles

ii=1;
Plot[Ek[ii][t],{t,0,tMax},PlotRange->All,PlotStyle->{Black,Thick}]

After averaging over N particles one can see the chaotization.

Plot[{EkxAvr[t],EkyAvr[t],EkAvr[t]},{t,0,tMax},PlotRange->All,PlotStyle->{{Blue,Thick},{Red,Thick},{Black,Thick}}]

Energy dispersion increases and approaches the asymptotic value (Subscript[k, B])T=1 in our case. This is relaxation due to collisions between particles.

Plot[ΔEk[t],{t,0,tMax},PlotRange->All,PlotStyle->{Black,Thick},PlotPoints->300]