Function Repository Resource:

HardSphereSimulation (2.0.0) current version: 2.1.0 »

Simulate identical unit-mass hard spheres moving in an N-dimensional box

Contributed by: Matt Kafker and Christopher Wolfram

ResourceFunction["HardSphereSimulation"][assoc]

generates a simulation of the motion of a system of identical unit-mass hard spheres with simulation parameters specified by assoc.

Details

Input assoc must include the following:

"Positions"

initial particle positions

"Velocities"

initial particle velocities

"BoxRadius"

half-side-length of simulation box

"ParticleRadius"

radius of spheres

"Steps"

number of simulation time steps

"dt"

duration of a single time step

"BoundaryCondition"

interaction between particles and container walls

"Output"

simulation output

The dimensions of vectors in "Positions" and "Velocities" determine the spatial dimension of the system (they can be any positive integer dimension).

The box is centered on the origin.

Possible named boundary conditions include:

"Reflecting"

particle velocities reflect upon contact with container walls

"Periodic"

particles are translated to opposite side of box upon reaching container walls

"Infinite"

particles do not alter motion upon reaching container walls

Possible named outputs include:

"PositionsByTime"

list of all particle positions at every time step

"VelocitiesByTime"

list of all particle velocities at every time step

"SpeedsByTime"

list of all particle speeds at every time step

"TrajectoriesByTime"

list of all particle positions and velocities at every time step

"PositionsByParticle"

list of positions over time for each particle

"VelocitiesByParticle"

list of velocities over time for each particle

"SpeedsByParticle"

list of speeds over time for each particle

"TrajectoriesByParticle"

list of positions and velocities over time for each particle

"Collisions"

list of particle-particle collisions

"All"

list containing “TrajectoriesByTime” output as the first entry, and “Collisions” as the second

"Visualize"

Graphics object with a Manipulate slider for visualizing the time-evolution of the system in one, two, or three spatial dimensions

"KineticEnergy"

list of total kinetic energy at every time step

"FinalFrame"

Graphics object depicting the state of the system at the final time step

"CausalGraph"

Graph object with nodes representing collisions and directed edges indicating that one particle was contiguously involved in two collisions

Each collision takes the form {particle_i,particle_j,time of collision}.

The option "Random" may be given for pos, in which case particle positions are chosen at random within the box.

Examples

Basic Examples (2)

Find the position of two hard disks over five time steps:

In[1]:=

ResourceFunction["HardSphereSimulation", ResourceVersion->"2.0.0"][<|"Positions" -> {{-2, -2}, {2, 2}}, "Velocities" -> {{1, 0}, {-1, 0}}, "BoxRadius" -> 5, "ParticleRadius" -> 1, "Steps" -> 5, "dt" -> 1, "BoundaryCondition" -> "Reflecting",
"Output" -> "PositionsByTime"|>]

Out[1]=

Get the velocities over time:

In[2]:=

Out[2]=

Get the speeds over time:

In[3]:=

Out[3]=

Get the trajectories over time:

In[4]:=

Out[4]=

Return the positions grouped by particle:

In[5]:=

Out[5]=

Get the velocities by particle:

In[6]:=

Out[6]=

Get the speeds by particle:

In[7]:=

Out[7]=

Get the trajectories by particle:

In[8]:=

Out[8]=

Scope (5)

Use HardSphereSimulation to visualize many hard spheres moving in a box:

In[9]:=

ResourceFunction["HardSphereSimulation", ResourceVersion->"2.0.0"][<|"Positions" -> 3.2 N[CirclePoints[10]], "Velocities" -> N[CirclePoints[10]], "BoxRadius" -> 10, "ParticleRadius" -> 1, "Steps" -> 200, "dt" -> 0.1, "BoundaryCondition" -> "Reflecting", "Output" -> "Visualize"|>]

Out[9]=

Use "Random" to generate random initial positions in the box:

In[10]:=

ResourceFunction["HardSphereSimulation", ResourceVersion->"2.0.0"][<|"Positions" -> "Random", "Velocities" -> RandomPoint[Ball[{0, 0}, 1], 15], "BoxRadius" -> 5, "ParticleRadius" -> 1, "Steps" -> 500, "dt" -> 0.1, "BoundaryCondition" -> "Reflecting", "Output" -> "Visualize"|>]

Out[10]=

Visualization works in 1D and 3D as well:

In[11]:=

ResourceFunction["HardSphereSimulation", ResourceVersion->"2.0.0"][<|"Positions" -> "Random", "Velocities" -> RandomPoint[Ball[{0}, 1], 15], "BoxRadius" -> 30, "ParticleRadius" -> 1, "Steps" -> 500, "dt" -> 0.1, "BoundaryCondition" -> "Reflecting", "Output" -> "Visualize"|>]

Out[11]=

In[12]:=

ResourceFunction["HardSphereSimulation", ResourceVersion->"2.0.0"][<|"Positions" -> "Random", "Velocities" -> RandomPoint[Ball[{0, 0, 0}, 1], 10], "BoxRadius" -> 3, "ParticleRadius" -> 1, "Steps" -> 500, "dt" -> 0.1, "BoundaryCondition" -> "Reflecting", "Output" -> "Visualize"|>]

Out[12]=

Periodic boundary conditions identify opposite walls of the box with one another:

In[13]:=

Out[13]=

Infinite boundary conditions eliminate particle-wall interactions altogether:

In[14]:=

Out[14]=

Applications (3)

Particle speeds should obey Maxwell-Boltzmann distribution:

In[15]:=

$Module[{speeds, T}, (*Run hard sphere simulation*) speeds = Flatten@ResourceFunction["HardSphereSimulation"][ <|"Positions" -> "Random", "Velocities" -> RandomPoint[Ball[{0, 0}, 1], 15], "BoxRadius" -> 5, "ParticleRadius" -> 1, "Steps" -> 5000, "dt" -> 0.1, "BoundaryCondition" -> "Reflecting", "Output" -> "SpeedsByTime"|>]; (*Equipartition/kinetic temperature.*) T = 1/2 Mean[speeds^2]; Show[ (*Empirical histogram of speeds*) Histogram[speeds, Sequence[{0.1}, "PDF", Frame -> True, FrameLabel -> Map[Style[#, Black, 14]& , {"v", "\[Rho](v)"}], PlotLabel -> Style[ "Speed Distribution with Maxwell-Boltzmann Fit", Black, 15], ImageSize -> Large]], (*Plot of 2D Maxwell-Boltzmann distribution with best- fit temperature.*) Plot[1/T v E^(-(v^2/(2 T))), {v, 0, 2}] ] ]$

Out[15]=

Extract a collision list from the simulation, where collisions take the form {particle₁,particle₂,time of collision}:

In[16]:=

With[{boxSize = 10, rad = 1},
ResourceFunction["HardSphereSimulation"][
<|"Positions" -> Map[{#, 0} &, Range[-(boxSize - 2 rad), (boxSize - 2 rad), 4*rad]],
"Velocities" -> Evaluate[
Table[{KroneckerDelta[n, 1], 0}, {n, Length[Range[-(boxSize - 2 rad), (boxSize - 2 rad), 4*rad]]}]],
"BoxRadius" -> boxSize, "ParticleRadius" -> rad, "Steps" -> 200,
"dt" -> 0.1, "BoundaryCondition" -> "Reflecting", "Output" -> #|>] & /@ {"Visualize", "Collisions"}
] // Column

Out[16]=

From particle-particle collisions, construct the causal graph:

In[17]:=

Out[17]=

Construct a causal graph for a disordered system (ensuring initial positions do not overlap to avoid spurious collisions):

In[18]:=

Module[{num = 15, boxSize = 5, rad = 1, InitPos}, (*Ensure initial positions are not overlapping by running a short simulation first.*)
InitPos = Last@ResourceFunction["HardSphereSimulation", ResourceVersion->"2.0.0"][<|"Positions" -> "Random", "Velocities" -> RandomPoint[Ball[{0, 0}, 1], num],
"BoxRadius" -> boxSize, "ParticleRadius" -> rad, "Steps" -> 50, "dt" -> 0.1, "BoundaryCondition" -> "Reflecting", "Output" -> "PositionsByTime"|>]; (*Run a simulation and extract the causal graph*)
ResourceFunction["HardSphereSimulation"][
<|"Positions" -> InitPos, "Velocities" -> RandomPoint[Ball[{0, 0}, 1], num],
"BoxRadius" -> boxSize, "ParticleRadius" -> rad, "Steps" -> 100, "dt" -> 0.1, "BoundaryCondition" -> "Reflecting", "Output" -> "CausalGraph"|>] // Graph[#, VertexLabels -> None] &
]

Out[18]=

As the nodes of the causal graph represent collisions between two particles and the edges represent individual particles, we expect all nodes to have indegree and outdegree 2, besides from the graph boundary (i.e., the beginning and end of the simulation):

In[19]:=

Out[19]=

Properties and Relations (2)

Total kinetic energy of spheres should be conserved in time. Show that the error stays within machine precision:

In[20]:=

$ListLinePlot[# - #[[1]] &[ ResourceFunction["HardSphereSimulation", ResourceVersion->"2.0.0"][<|"Positions" -> "Random", "Velocities" -> RandomPoint[Ball[{0, 0}, 1], 15], "BoxRadius" -> 5, "ParticleRadius" -> 1, "Steps" -> 500, "dt" -> 0.1, "BoundaryCondition" -> "Reflecting", "Output" -> "KineticEnergy"|>] ], Sequence[ Frame -> True, FrameLabel -> {"Time Step", "E-\!$\*SubscriptBox[\(E$, $initial$]\)"}, RotateLabel -> False, PlotLabel -> "Energy Conservation vs. Time", PlotRange -> {{0, 500}, Automatic}]]$

Out[20]=

Momentum should be conserved during collisions between particles:

In[21]:=

With[{boxSize = 10, rad = 1}, ResourceFunction["HardSphereSimulation", ResourceVersion->"2.0.0"][<|
"Positions" -> Map[{#, 0} &, Range[-(boxSize - 2 rad), (boxSize - 2 rad), 4*rad]], "Velocities" -> Table[{KroneckerDelta[n, 1], 0}, {n, Length[Range[-(boxSize - 2 rad), (boxSize - 2 rad), 4*rad]]}],
"BoxRadius" -> boxSize, "ParticleRadius" -> rad, "Steps" -> 400, "dt" -> 0.1, "BoundaryCondition" -> "Reflecting", "Output" -> "Visualize"|>]]

Out[21]=

Possible Issues (1)

Particle overlaps are resolved pairwise, so situations involving multiple overlapping spheres may lead to unexpected results (such as a disruption of symmetry):

In[22]:=