Hello, bonjour, nomoshkar, ola, estou?
Welcome to the second week of Statistical Mechanics: Algorithms and Computations
from the Physics Department of Ecole Normale Superieure.
Last week, we studied simple problems in statistics
children on the beach, discrete pebble games,
people walking around on the Monte Carlo heliport.
In this second week, we move on to study and to simulate
model systems in physics
particle with positions, velocities and interactions
using a statistical approach
that this is at all possible and exact
is one of the great legacies of physics.
The entire program of this week
is to motivate, define and illustrate
the statistical approach in the very famous model
of two dimensional hard spheres
so called hard disks.
This model is at the center of physics
in particular of statistical mechanics
Note, when we speak about hard disks
we mean idealizations of this object
not of computer hard disks
the memories in our laptop
In this lecture - Lecture 2 -
we will introduce the model of hard disks
and also will start simulating
hard disks, which are idealizations of billiard balls
move as free particles whenever they are not in contact with each other
or with a wall
and obey simple reflection rules
whenever they are in contact with each other, as you see here
in an actual simulation of Newtonian dynamics
Building on what we did last week
We will also introduce
direct sampling and Markov chain sampling algorithms
we thus scrutinize the interface
between classical Newtonian mechanics
and statistical mechanics.
We will see that in the statistical approach
any two configurations of hard disks
must be sampled with equal probability
just as we did last week
for pebbles on the heliport
or on the Monte Carlo beach.
In Tutorial 2,
we will connect the equiprobability concept
- a basic version of the celebrated Boltzmann distribution -
with the central notions of statistical mechanics
namely the partition function, the free energy,
the phase space volume
and the virial expansion.
In Homework session 2
our main concern will be - just like last week - on practical computation.
We will actually check that the Newtonian dynamics
of four hard disks perfectly agrees
with the statistical approach.
We will then analyze and interpret
the output for this very simple problem
The material of this week is roughly as difficult as the one of last week.
We'll move to the central notion of statistical mechanics
but we will discover and develop all the concepts using our algorithmic approach.
We will remain perfectly self-contained
and you don't have to have already taken a course
in statistical physics or thermal physics in order to follow.
So, let's get started with week 2 of Statistical Mechanics: Algorithms and Computations.
What you see here
is a configuration of four disks
red, blue, yellow and green
with given positions and velocities
at time t=0
in a box going in x from 0 to 1 and in y from 0 to 1.
Here is the time evolution of the disks
moving about the box, like billiard balls.
The dynamics of this system
- for example from a configuration at t=0 -
consists in straight line evolution
up to the next event.
This event may be a wall collision,
with incoming angle equal to the outgoing angle.
Like with real disks, there are also pair collisions
where the velocities change
according to the rules of an idealized billiard
without the complications of real systems
such as torque or friction.
The collision rules are quite easy to write down.
Let us now describe in detail
the event-driven molecular dynamics algorithm
invented by Alder and Wainwright in 1957.
It solves the Newtonian dynamics of the system without any approximation.
Starting from an initial configuration
- at t=0 -
we must compute the next event.
This event may be a wall collision
or a pair collision.
If it is a wall collision
its time is given by the minimum of the wall collision times
of the red, the blue, the yellow and the green disks,
taken individually.
Look here, the red disk
is moving down and to the left.
Its wall-collision time - disregarding the three other particles -
is the minimum of the times at which it would hit the line x=0
and y=0.
The second type of events is a pair collision,
for example here of the red and the blue disks.
A pair collision takes place when at some future time
the distance between the two disks is equal to twice their radius.
We can compute this time using a quadratic equation.
When the two disks are approaching each other,
the pair collision time
is the minimum of the two positive solutions of this quadratic equation,
otherwise the pair collision time is infinite.
The minimum of the six pair collision times
and the four wall collision times
gives us the next event in the system,
here, at time t=0.45813
a number that we can compute to as high precision as we want.
We can then click through the events.
Here the first event, a pair collision between the red and the yellow disks,
closely followed by the second event,
a wall collision of the blue disk,
then the third event, the fourth event, the fifth event, the sixth event, and so on..
To conclude our description of the event-driven molecular dynamics algorithm,
we realize that we must not collect events,
but rather make a movie
going from time t to t+delta_t
t+2 delta_t, .. or to make things simple, from t=0 to t=1, 2, 3, 4, 5, 6 and so on..
Let us now ask a few questions
Question 1
why is Alder and Wainwright's algorithm for molecular dynamics called event-driven?
Is it because
Answer 1: its invention was a major event in the history of science?
or is it (Answer 2) because the algorithm moves forth in a sequence of collisions called events?
Of course the correct answer is Answer 2
although it is true that the invention of molecular dynamics turned out to be a real big thing
Question 2
consider the sequence of configurations at regular time intervals,
that we just saw rushing through from time t=0 to t=1, 2, 3 and so on, up to infinity/
On which fraction of these pictures do we see two disks in contact
or one disk in contact with the wall?
Answer 1:
the probability is equal to zero, exactly.
Answer 2: the probability is very small,
between 1% and 10%,
depending on the time interval between frames.
Answer 3: the probability of contact is equal to 100%,
because the algorithm is driven forward by collisions.
Answer 1 is correct.
Common sense and powerful mathematical theorems
tell us that there is a finite number of events per unit time interval
so the probability to have an event at time equal to 1.00000 and so on
or at time 16.000000 and so on is exactly equal to zero.
You will not see any collision event
in these frames taken at integer times.
Question 3:
in the event-driven molecular dynamics algorithm
we suppose that a pair collision
or a collision of a particle with a wall
constitute the next event.
Isn't it possible that the next event
is rather a collision of three disks
or a pair collision
taking place at exactly the same time as another wall collision?
No, this is not possible
Common sense and powerful mathematical theorems
again - just like in Question 2 - tell us
that the probability of one event taking place exactly at the same time as another event is equal to zero.
With finite precision arithmetic on our computers
these exceptional events may have to be taken into account in the treatment of errors.
The event-driven molecular dynamics algorithm
exactly solves the hard disks equations of motion
under the assumption
that the calculation of collision times
velocity changes, positions and so on.. is done with infinite precision.
This cannot really be achieved on a computer
but the question arises whether it matters
that in Python - as in other languages - real numbers are truncated to finite precision.
It is easiest to answer this question
by pitting different versions of the same event-driven molecular dynamics algorithm
against each other.
The two versions of the four disks molecular dynamics calculations
that are identical at t=0
remain very similar at initial collisions
but then completely get out of step
after a few dozens collisions.
This is an extremely small number
compared to the millions or so
collisions that we can compute on our computers
per second.
This is quite strange.
Our results for computing times as small as a few microseconds
are exact yet uncontrolled.
We may drive up the precision levels
using special number formats,
like 100 significant digits, 1000 significant digits,
1000000 significant digits,
that are available in Python as in other computer languages.
But this strategy cannot defeat the onset of chaos,
that is: cure the extreme sensitivity of our calculations
to the numerical precision.
It is impossible
to control an event-driven molecular dynamics calculation
for a few billions iterations
Chaos, in a hard disks system,
has its origin in the negative curvature of the disks' surface.
This magnifies
small differences of the trajectory
at each pair collision
and is very well known to whoever has played billiards.
On the other hand, chaos
is the main reason why even for finite system of hard disks
we can describe it by statistical mechanics, as we will discuss in a few minutes.
For the moment, we considered Newtonian dynamics, which is completely deterministic.
Before moving on to the statistical approach,
please take a moment to download, run and modify the programs.
On the Coursera website,
you'll find the program event_disks_box.py
which computes the time evolution
and the program event_disks_box_movie.py
which produced the nice graphic output that you saw all along this section.
Check them out for yourself.
Both the wall collisions and the pair collisions
are programmed in a few lines,
and there are absolutely no approximations
besides the finite precision of our calculations.
The event-driven molecular dynamics algorithm
of the last section
can be written down in a few dozens lines
- as we saw -
and it allows us to follow the ballet
of disks approaching and flying away from each other,
along intricate, even unpredictable trajectories.
In doing so however,
we engage on a computational project
which is much more complicated than what we need.
To continue, we again consider the configurations
that we generated just a few moments ago.
From a physical point of view, these configurations
have kinetic energy
- as we will discuss two weeks from now -
but their potential energy is zero if they have no overlap,
and is infinite if they have an overlap.
Configurations with overlaps
are simply forbidden.
Ludwig Boltzmann - at the end of the 19th century -
boldly conjectured that any configurations with the same energy
also have the same probability to appear
For hard disks, it is the quintessence
of Boltzmann's statistical mechanics
that all configurations a and b of hard disks
have the same statistical weight pi(a)=pi(b).
We already encountered this equation last week
in our Monaco games
and a and b were pebble positions.
We sampled this pebble positions with a random number in x and a random number in y.
Now, we generalize from the children's algorithm on the beach
to a direct sampling Monte Carlo algorithm for hard disks.
We first place the red disk
as a random number in x between sigma (the radius) and 1-sigma
and y: another random number between sigma and 1-sigma.
The red disk has been placed at a random position.
We do the same for the blue disk, the green disk and the yellow disk.
Fantastic! Congratulations!
We have placed the four disks at random positions.
In total, we have produced a random configuration of the four disks,
as called for by Ludwig Boltzmann.
Let's do it again.
Let's place the red disk at a new random position inside the box.
Then the blue disk, then the green disk, ouch!
We have generated an overlap. What should we do now?
You might be tempted by simply taking away the green disk
and trying to place it again
but this is wrong.
The correct solution is called 'tabula rasa'
It consists in wiping out the whole configuration
and in starting again from an empty box
with the red, the blue, the green, the yellow disks
until you succeed.
We will spend a good part of this week's tutorial
to make sure we completely understand
why we have to implement the tabula rasa rule
and why this allows us
to sample a random configuration of hard disks.
Before moving on, please take a moment
to download, run and modify the programs.
On the Coursera website, you'll find the program direct_disks_box.py
that implements the direct sampling Monte Carlo algorithm.
It samples one configuration of four disks
its key element are the random positions in x
and in y, between sigma and 1-sigma.
Four times we generate such a random position
and check for the minimum distance
with the disks we placed earlier.
If this minimum distance is smaller than (2 x sigma)
- twice the radius - we know that we have an overlap.
We then break out of the construction,
and start again, from scratch.
On the website, you also find the program
direct_disks_box_multirun.py
that generates many such configurations,
and as usual the program direct_disks_box_movie.py,
that generated nice graphics output that you see here.
I can ensure you that the sequence of configurations that you see here
is equivalent to the output of the molecular dynamics calculation
if you erase time information, scramble the sequence
and also erase the velocities.
Then, these two sequences are equivalent
and we have arrived at the center of statistical mechanics.
Direct sampling for hard disks works only at low densities and for small particle numbers
as we will discuss in detail in this week's tutorial
To make progress, we move to the more general Markov chain Monte Carlo algorithm
Let me first remind you of the algorithm markov_pi that we discussed last week
From a configuration with position at (x,y)
we moved to a configuration b
by little moves delta_x and delta_y that could be positive or negative.
Detailed balance imposed
that the transition probability from a to b
had to be equal to the transition probability from b to a
in order to satisfy equiprobability
pi(a) = pi(b).
Let us now apply the same reasoning to hard disks.
Consider here configuration a of four hard disks
the move consists in picking one of the four disks - for example the blue disk -
and making a little displacement delta_x and delta_y - that can be again positive or negative - to a configuration b.
The detailed balance condition is satisfied
- with pi(a)=pi(b) -
if we move with the same probability from a to b
as we move from b to a.
It can also happen that the move from configuration a is rejected, as it creates an overlap.
In that case, we remain at the configuration a for the next iteration.
In Python this gives the following program markov_disks_box.py
that you should download from the Coursera website, run and modify.
The key element is the random choice of one disk
whose coordinates we modify a little bit.
The new configuration is tested for overlap
and if two disks are getting too close
or if we have an overlap with the wall
we simply remain with the configuration a
otherwise let's add the new configuration
This precise algorithm for this model of hard disks
was introduced by Metropolis et al.
in a famous article of 1953.
On the website you'll also find the program markov_disk_box_movie.py
that produces nice graphics output shown here.
You'll see that there are correlations between subsequent configurations
and sometimes two subsequent configurations are identical
This corresponds to piles of pebbles that we created the last week.
This algorithm samples the same distribution as the direct_disks algorithm.
As discussed in last week's tutorial by Vivien,
we must check - in addition to the detailed balance condition -
that we satisfy irreducibility
and aperiodicity.
Aperiodicity is trivial to check, for our algorithm
but in order to satisfy irreducibility
we must avoid situations as the one shown here
where the blue disk remains in the upper right corner up to infinite times.
For small systems, with a small number of particles,
we must go to small disk radii
as shown here where the disks move around the whole system
But for larger system this is not a problem
and the markov_disks algorithm remains irreducible up to very high densities.
In this week's homework session, you will show - using numerical simulations -
and up to the numerical precision
that the event-driven molecular dynamics calculations
and the Monte Carlo simulation give equivalent result for thermodynamic quantities.
This means that the equiprobability hypothesis formulated by Boltzmann is satisfied.
This equivalence between Newtonian deterministic mechanics
and statistical mechanics is called the ergodicity hypothesis.
It certainly does not follow from simple principles
like time-reversal invariance or detailed balance.
The proof of equiprobability
has presented a formidable mathematical challenge.
It was once believed to hold only in the limit of infinite number of particles.
Mathematical research has gone very far in actually proving the ergodicity hypothesis
for the special case of hard disks and hard spheres and already for finite number of particles.
The first milestone was obtained by Sinai in 1970.
He proved that for two disks with periodic boundary conditions
the hypothesis was actually satisfied.
Very recently, mathematicians like Simanyi have proven
that a finite number number of hard disks have to satisfy the ergodicity hypothesis.
So you can be sure that what you will do for your homework this week stands on very solid mathematical ground.
In conclusion we have plunged, in Lecture 2 of Statistical Mechanics: Algorithms and Computations,
into the foundations of statistical mechanics.
The hard disk model that we have studied
is special in that all legal configurations
have the same potential energy and therefore the same statistical weight.
This model has been instrumental for the development of statistical mechanics
and also for the development of molecular dynamics
and of the Monte Carlo method.
Next week, we will study the very rich phenomenology of this model
which is at the basis of the physics of liquids, among others.
In the meantime, have fun with Tutorial 2
and with this week's homework.
It remains for me to thank you for your attention for this session
See you again in other sessions of Statistical Mechanics: Algorithms and Computations.