Function Repository Resource:

HardSphereSimulation

Source Notebook

Simulate 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 hard spheres with simulation parameters specified by assoc.

Details

Input assoc must include the following:

"Positions"

initial particle positions

"Velocities"

initial particle velocities

"BoxSize"

side-length of simulation box

"Steps"

number of simulation time steps

"StepSize"

duration of a single time step

"BoundaryCondition"

interaction between particles and container walls

"Output"

simulation output

Input assoc may also include:

"ParticleRadius"

radius of spheres

"ParticleMass"

mass of spheres

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

"ParticleMass" can either be a single number for identical-mass spheres or a list for spheres of different masses.

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

"Temperature"

list of temperature at every time step

"CausalGraph"

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

Multiple outputs can be provided as a List

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

The option "RandomNonOverlapping" may be given for pos, in which case a random configuration of non-overlapping particle positions is chosen within the box.

Examples

Basic Examples (2)

Simulate the motion of two hard disks over five time steps, and extract the trajectory data in various ways. Get the positions over time:

In[1]:=

ResourceFunction[
"HardSphereSimulation"][<|"Positions" -> {{-2, -2}, {2, 2}}, "Velocities" -> {{1, 0}, {-1, 0}}, "BoxSize" -> 10, "ParticleRadius" -> 1, "Steps" -> 5, "StepSize" -> 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]=

Get the positions grouped by particle:

In[5]:=

Out[5]=

Get the velocities grouped by particle:

In[6]:=

Out[6]=

Get the speeds grouped by particle:

In[7]:=

Out[7]=

Get the trajectories grouped by particle:

In[8]:=

Out[8]=

Scope (7)

Use trajectory data to visualize frames from a hard sphere simulation:

In[9]:=

Module[{trajectory}, (*Compute trajectories*)
trajectory = ResourceFunction[
"HardSphereSimulation"][<|
"Positions" -> Catenate@Table[{-14 + 2.5 i, 14 - 2.5 j}, {i, 4}, {j, 4}], "Velocities" -> Table[.5 {1, -1}, 16], "BoxSize" -> 30, "StepSize" -> 0.1, "ParticleRadius" -> 1, "Steps" -> 2000, "BoundaryCondition" -> "Reflecting", "Output" -> "TrajectoriesByParticle"|>]; (*Display snapshots of disks with velocity vectors*)
GraphicsGrid[Partition[Table[
Graphics[{Style[Disk[#[[1]], 1], Orange], Arrow[{#[[1]], #[[1]] + 4 #[[2]]}]} & /@ trajectory[[All, t]],
PlotRange -> 15, Frame -> True, FrameTicks -> None, ImageSize -> 100], {t, 1, 2000, 100}], 5], Spacings -> {5, 5}]]

Out[9]=

Use "Visualize" to generate an interactive visualization of many hard spheres moving in a box:

In[10]:=

ResourceFunction[
"HardSphereSimulation"][<|"Positions" -> 3.2 N[CirclePoints[10]], "Velocities" -> N[CirclePoints[10]], "BoxSize" -> 20, "ParticleRadius" -> 1, "Steps" -> 200, "StepSize" -> 0.1, "BoundaryCondition" -> "Reflecting", "Output" -> "Visualize"|>]

Out[10]=

Use "RandomNonOverlapping" to generate random non-overlapping initial positions in the box:

In[11]:=

ResourceFunction[
"HardSphereSimulation"][<|"Positions" -> "RandomNonOverlapping", "Velocities" -> RandomPoint[Ball[{0, 0}, 1], 15], "BoxSize" -> 10, "ParticleRadius" -> 1, "Steps" -> 200, "StepSize" -> 0.1, "BoundaryCondition" -> "Reflecting", "Output" -> "Visualize"|>]

Out[11]=

Visualization works in 1D:

In[12]:=

ResourceFunction[
"HardSphereSimulation"][<|"Positions" -> "RandomNonOverlapping", "Velocities" -> RandomPoint[Ball[{0}, 1], 15], "BoxSize" -> 60, "ParticleRadius" -> 1, "Steps" -> 500, "StepSize" -> 0.1, "BoundaryCondition" -> "Reflecting", "Output" -> "Visualize"|>]

Out[12]=

Visualization works in 3D as well:

In[13]:=

ResourceFunction[
"HardSphereSimulation"][<|"Positions" -> "RandomNonOverlapping", "Velocities" -> RandomPoint[Ball[{0, 0, 0}, 1], 10], "BoxSize" -> 6,
"ParticleRadius" -> 1, "Steps" -> 500, "StepSize" -> 0.1, "BoundaryCondition" -> "Reflecting", "Output" -> "Visualize"|>]

Out[13]=

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

In[14]:=

Out[14]=

Infinite boundary conditions eliminate particle-wall interactions altogether:

In[15]:=

Out[15]=

Multiple outputs can be returned at once:

In[16]:=

Out[16]=

Applications (3)

Particle speeds should obey Maxwell-Boltzmann distribution:

In[17]:=

$Module[{speeds, Temp}, {speeds, Temp} = {Flatten[#[[1]]], First[#[[2]]]} &@ ResourceFunction["HardSphereSimulation"][ <|"Positions" -> "RandomNonOverlapping", "Velocities" -> RandomPoint[Ball[{0, 0}, 1], 15], "BoxSize" -> 10, "ParticleRadius" -> 1, "Steps" -> 5000, "StepSize" -> 0.1, "BoundaryCondition" -> "Reflecting", "Output" -> {"SpeedsByTime", "Temperature"}|>]; Show[ (*Empirical histogram of speeds*) Histogram[speeds, {0.1`}, "PDF", Frame -> True, FrameLabel -> (Style[#1, Black, 14] &) /@ {"v", "\[Rho](v)"}, PlotLabel -> Style["Speed Distribution with Maxwell-Boltzmann Fit", Black, 15], ImageSize -> 400], (*Plot of 2D Maxwell-Boltzmann distribution with best- fit temperature.*) Plot[1/Temp v E^(-(v^2/(2 Temp))), {v, 0, 2}] ] ]$

Out[17]=

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

In[18]:=

ResourceFunction["HardSphereSimulation"][
<|"Positions" -> {{-8, 0}, {-4, 0}, {0, 0}, {4, 0}, {8, 0}},
"Velocities" -> {{1, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}},
"BoxSize" -> 20, "ParticleRadius" -> 1, "Steps" -> 200, "StepSize" -> 0.1, "BoundaryCondition" -> "Reflecting", "Output" -> "Collisions"|>]

Out[18]=

From particle-particle collisions, construct the causal graph:

In[19]:=

Out[19]=

Construct a causal graph for a disordered system:

In[20]:=

Out[20]=

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[21]:=

Out[21]=

Properties and Relations (2)

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

In[22]:=

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

Out[22]=

Momentum should be conserved during collisions between particles. Collision of equal mass particles:

In[23]:=

ResourceFunction[
"HardSphereSimulation"][<|
"Positions" -> {{-8, 0}, {-4, 0}, {0, 0}, {4, 0}, {8, 0}},
"Velocities" -> {{1, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}},
"BoxSize" -> 20, "ParticleRadius" -> 1, "ParticleMass" -> 2.0, "Steps" -> 200, "StepSize" -> 0.1, "BoundaryCondition" -> "Reflecting", "Output" -> "Visualize"|>]

Out[23]=

Collision of an incoming heavy particle with a much lighter stationary particle causes the light particle to go out with double the speed of the incoming one:

In[24]:=

$ResourceFunction[ "HardSphereSimulation"][<| "Positions" -> {{-2.1, 0}, {0, 0}, {2.1, 0}}, "Velocities" -> {{1, 0}, {0, 0}, {0, 0}}, "BoxSize" -> 50, "ParticleRadius" -> 1, "Steps" -> 4, "StepSize" -> 0.1, "BoundaryCondition" -> "Reflecting", "Output" -> "SpeedsByTime", "ParticleMass" -> {10^10, 10^5, 1}|>]$

Out[24]=

Collision of a very light particle with a much heavier stationary particle leads to reflection of the incoming particle:

In[25]:=

$ResourceFunction[ "HardSphereSimulation"][<|"Positions" -> {{0, 0}, {2.1, 0}}, "Velocities" -> {{1, 0}, {0, 0}}, "BoxSize" -> 20, "ParticleRadius" -> 1, "ParticleMass" -> {1, 10^5}, "Steps" -> 3, "StepSize" -> 0.1, "BoundaryCondition" -> "Reflecting", "Output" -> "VelocitiesByTime"|>]$

Out[25]=

Possible Issues (2)

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

In[26]:=