Hello, bonjour, ni hao, hola,
welcome to the second week of Statistical Mechanics: Algorithms and Computations
from the Physics department of Ecole Normale Superieure.
This week's subject is the physics of hard disks.
Remember, hard disks: two-dimensional idealization of billiard balls,
not the hard disk in your computer.
These hard disks - or at least their idealization -
were studied with two approaches.
Molecular dynamics, the solution of Newtonian deterministic mechanics,
and Monte Carlo, the solution of Boltzmann statistical mechanics.
As discussed in the lecture,
for hard disks, Boltzmann statistical mechanics means equal weights
for all legal configurations.
We also discussed that Boltzmann and Newton
are equivalent for all thermodynamic properties,
that means for all properties that do not explicitly depend on time.
In this tutorial, we will study direct sampling for hard disks.
We will find out that
although the program direct_disks.py is not a very powerful program
it is a VIP, a very important program,
as it contains important connections between physics and computation.
In what follows,
Michael will scrutinize the equal probability principle
and the tabula rasa rule.
Then Alberto will explain the profound link
between computation and physics,
or rather, between the acceptance probability
of the direct sampling algorithm
and the partition function of the physical system.
Finally, Vivien will explore an analytical method
for studying the hard disks system,
the virial expansion.
The virial expansion has existed a long time
before numerical calculations became available.
Ludwig Boltzmann himself, the founder of statistical mechanics
wrote a famous paper in 1874
on the virial expansion of the hard sphere system.
So you see, we also keep to the center of the road of statistical mechanics.
In this week's lecture, Lecture 2, Werner explained
that the equiprobability principle is one of the pillars of statistical physics.
This principle states that the probability pi(a)
of a configuration is a function
of the energy of the configuration, that is
pi(a) equals pi of the energy of a.
For hard disks, all allowed configurations
have the same energy E=0.
This gives the same statistical weight to each configuration
pi(a) = pi(b) = pi(c).
The Monte Carlo algorithm must sample these configurations with equal probability,
this was achieved by the direct sampling Monte Carlo algorithm.
Remember that the algorithm direct_disks_box.py
first places one disk at a randomly chosen position
inside the box, then a second one, then a third one,
but in the case of an overlap applies the tabula rasa rule
and wipes out the entire configuration,
and starts afresh.
But why don't we simply take the last disk
(the one that caused the overlap)
and try to put it a second time or a third time
until it finally finds its place?
This algorithm is called random sequential deposition,
and it is important for adhesion and catalysis,
but not for equilibrium statistical physics.
Its configurations are not equally probable
as we will now see, using two arguments.
To do so, let us first consider a discrete hard rods model,
a one dimensional discrete version of hard disks.
On the line shown here,
suppose that the rod centers
can be placed on the sites 0, 1, 2, 3 and 4.
Now assume that there also is an overlap condition,
just like for the two-dimensional hard disks.
We say that two adjacent rod centers
must differ by at least three,
this means that there must be at least two free sites
in between two adjacent rods.
For example we cannot place rods on the sites 0 and 2,
but there's no problem to place them on the sites 1 and 4.
Let us do a simple simulation of random sequential deposition,
where we impose that we end up with two rods:
the red rod and the blue rod.
We place the red rod, then we place the blue rod,
but we have to pull it off again because there's an overlap
and we put it again, we pull it off again, we put it again
until we succeed.
Here you can see a second simulation.
We now analyze the algorithm:
we first place the red rod
and it'll fall with the probability of 1/5
on any of the sites 0, 1, 2, 3, 4.
Rather, with a probability of 1/4
because there is no two-disks configuration
with one of the disks on site 2.
Now, if the red rod is on site 0,
the blue rod can be placed with a probability of 1/2
on the sites 3 and 4.
This gives the allowed configurations a and b.
The same reasoning is valid for the configurations e and f.
These four allowed configurations
have the probability 1/4 times 1/2,
equal 1/8.
In contrast, if the red rod is placed
either on site 1 or on site 3,
there's only one site left
for the blue rod to go.
This leads to the configurations c and d,
which have the probabilities 1/4 times 1, equals to 1/4.
Clearly, we have violated the equiprobability principle.
What would happen if we apply - instead of random sequential deposition -
the tabula rasa strategy of the Monte Carlo algorithm?
In Python, this would give the program shown here
as you will discover in the following little quiz.
Question 1: how many configurations
(legal and illegal ones)
of red and blue rods can be generated by the program direct_discrete.py?
The answer is 25,
as there are five positions for the red rod and five positions for the blue rod to be placed.
Question 2:
how many legal configurations are among these 25 configurations?
Remember that there must be at least two free sites
in between the centers of the red and the blue rods.
The answer is 6, as only the 6 configurations
a, b, c, d, e and f (that we discussed earlier)
survive the tabula rasa wipe out.
Question 3:
what is the probability of the legal configurations a,...,f?
The legal configurations a,...,f have the probabilities
pi(a) = pi(b) = .. = 1/6.
The direct sampling Monte Carlo algorithm
generates the 25 legal and illegal configurations
with equal probabilities 1/25.
The legal configurations a,..,f have
equal probabilities before the tabula rasa wipe out
and consequently have equal probability after.
The argument we just used for the discrete model
applies also for the four hard disks
in the two dimensional box, if we reprogram
direct_disks_box.py
so that it places the four disks at randomly chosen positions
without checking for an overlap.
Now, there's no doubt that these four disks configurations
are uniformly distributed, in an eight-dimensional space
as we have two coordinates (namely x and y)
for each of the four hard disks.
In fact, as we do not check for overlaps,
each point in this eight-dimensional space
is sampled with equal probability.
It is only after we have generated all these configurations
that we cut out pieces.
In this program, the vector L
is a uniform random vector in an eight-dimensional hypercube
from sigma to (1-sigma) in all dimensions.
We then cut out pieces from this hypercube
with the overlap condition.
This is like in the discrete model, where we threw away
19 of the 25 configurations.
In the next section, Alberto will deepen our understanding
of the eight-dimensional space picture
for hard disks, but no longer with walls,
but with periodic boundary condition in the x and y directions.
Before following these exciting developments
you should download, run and modify the relevant programs.
On the Coursera website, you will find programs
for the discrete rods model.
First, the program for sequential random deposition
which also outputs a nice histogram, and as usual
also a random sequential deposition program with nice graphics.
You can also find the program for direct sampling of the four hard disks in the box,
if you have not already downloaded it in Lecture 2.
And even the slower program
which first samples all the positions, and checks for overlaps only in the end.
Make sure that you understand
that direct_disks_box and direct_disks_box_slow
are statistically speaking strictly identical in their output.
Let's move to direct sampling of hard disks
with periodic boundary conditions.
As you can see here,
when a disk is close to the right hand side
of the box, it will appear also on the left hand side
of the box. If you have two disks
they overlap not only if they're close together
but also if for example one disk is on the lower right corner
and the other disk is in the upper right corner
In a system with periodic boundary conditions, walls play no role
Here we show the multirun Python version with periodic boundary conditions.
Many legal configurations are generated
Check out for yourself that the function dist()
correctly implements the periodic boundary conditions.
Now we consider more than four disks.
Let's say N=16 disks, and we fix also the density eta,
which is the fraction of space occupied by disks.
Here, we show the output
for a density eta=0.3,
only 30% of the space is occupied by disks.
For our box of side 1,
the radius of the hard disks can be obtained by the relation
16 x pi x sigma^2
equal to 0.3.
Something chilling is happening:
we typically have to wipe off our table 100000 times
before having a single legal configuration.
This means that our algorithm is exact
but is very bad, because the acceptance ratio is very small.
For a system of 16 disks with a density of 0.3
the acceptance ratio is of the order of 5 x 10^(-6).
Michael spoke a few minutes ago
about the eight-dimensional space of configurations of four disks
These configurations are described by the eight coordinates of the disks
(x_1, y_1), (x_2, y_2), (x_3, y_3), (x_4, y_4).
Here we have 16 disks, which means that the dimension of the space is 32
Furthermore, if our box has a volume (or an area) equal to V
the 32-dimensional space has a volume V^16
or V^N, where N is the number of disks.
If the system has a density equal to zero, which means that the radius of the hard disks is equal to zero,
all the points of the configurational space are legal configurations.
So the number of legal configurations, which is the partition function of the system,
is V^N. At density 0, the acceptance ratio is 1.
We found that at density 0.3
the acceptance ratio is 5 x 10^(-6).
This means that the number of legal configurations
is 200000 smaller than the volume
of the original 32-dimensional space (V^N).
The number of legal configurations at density eta
- which is the partition function of the hard disks -
is equal to the product of the partition function at density 0
times the acceptance ratio at density eta.
The partition function is a central object in equilibrium statistical mechanics
and is a central object in Statistical Mechanics: Algorithms and Computations.
In a few moments, Vivien will discuss this in more details.
and will provide an approach (an analytical approach)
to compute the acceptance ratio and the partition function
for hard disks.
Before doing this, please take a moment to download, run and modify
the programs we discussed in this section
On the Coursera website, you will find the algorithm
for direct sampling with periodic boundary conditions.
This algorithm allows you to pick out legal configurations
from the hypercube of dimension V^N.
These configurations are as rare as diamonds,
and if you wish they are almost as beautiful.
Note that all coordinates are now random numbers between 0 and 1
because no wall has to be considered.
You can inspect these legal configurations, the diamonds of statistical physics
using, as usual, the graphic version of the program.
Let us go one step further and determine the acceptance ratio
of the direct sampling algorithm at any density eta.
You may be tempted to run the program direct_disks_multirun
for several runs and maybe on many computers for different values of eta.
But I'd rather look at the following configuration of N=16 points
in a box with periodic boundary conditions.
At zero density, that is to say for zero radius
it clearly leads to a valid configuration of hard disks
and this is also the case for a density equal to 0.0001
In contrast, at higher density the radius is larger,
and the points lead to a configuration of hard disks which is illegal.
In fact, from these 16 points in the box
you can obtain a hard disks configuration
for all values of the radius which are smaller
than the maximal radius, equal to half of the minimal distance between two of these points.
In this present configuration, the value sigma_max of this maximal radius is equal to 0.025
and it corresponds to a maximal density eta_max = 0.03143
You may consider the probability distribution p(eta_max)
associated to this maximal density
and from it - as displayed on this formula -
you can reconstitute the probability distribution pi_accept(eta)
of the acceptance ratio of the direct sampling algorithm.
You can find this formula on the fact-sheet associated to this class session.
For 16 disks you obtain the following graph
which has been obtained from a few millions realizations of the program.
It describes the acceptance rate as a function of the density eta
and you observe that it decays faster than exponentially with eta.
We will now obtain an analytical approximation
of the acceptance rate and of the partition function of hard disks.
It will be excellent in the limit of small densities, and valid for any value of N
To do so, let's consider the following N=16
points in the box with periodic boundary conditions
This configuration is described by the set of vectors
x_0, x_1, .., x_(N-1)
each of these vectors contains the two spatial coordinates of the point, x and y.
For relatively small radius, this configuration will be allowed
In general this configuration will be valid provided that the distances
between disk 0 and disk 1, disk 0 and disk 2 up to N-1 and N-2
are all larger than 2 sigma.
If one of these conditions is not valid,
then the configuration is illegal.
Following Alberto a few minutes ago, we will rewrite this condition in integral form, as follows.
This integral runs over the 2 N components of the N vectors of the configuration
and Upsilon is equal to 1 if there is an overlap, and 0 otherwise.
It is important to note that if one of these brackets is equal to zero
then the contribution of the configuration vanishes in the integral
This expression of the partition function is exact, but no one has been able to compute it analytically.
Let's consider all the brackets, each of them is equal to
1 if there is no overlap, and 0 otherwise.
We can expand and multiply out all these brackets
For instance if one picks the ones in each of these brackets
and then one computes the integral, one obtains V^N.
This is in fact the partition function - as Alberto showed -
of the ideal gas of hard disks of radius equal to zero.
To go to the next order, let's consider all the brackets
there are (N-1)N/2 of them,
so what you can do is to pick one of the Upsilon
in a single term, and to pick one in all the other brackets.
Doing so, you obtain a four dimensional integral, which is displayed here.
Thanks to the invariance by translation and the periodic boundary condition, it reduces to V times a two-dimensional integral.
We observe - as displayed in this figure - that the excluded region corresponds to a disk of radius 2sigma.
This allows to compute the integral as follows.
Now, we can come back to the full expression of the partition function,
and if we stop the evaluation at this term
we obtain the following expression.
For large N, it behaves as V^N
times the exponential of - (N-1) times eta.
This means that the acceptance probability decays exponentially with both the density and the number N of particles.
Remember that we have a computer program called direct_disks_any.py,
which allows to evaluate numerically the acceptance rate.
We can thus compare it to our prediction, and this is displayed on the following picture
As you observe the analytical prediction matches very well
the numerical results in the low density regime.
Let's now remark that the partition function is proportional to V^N
To obtain an extensive quantity, in physics we take the logarithm of this quantity
and you obtain the following result.
Let us now explain why it corresponds to a virial expansion.
To do so, we can first differentiate log(Z(eta)) with respect to V
and then multiply by V/N.
You obtain the following expression
the first term is a constant equal to 1
and it corresponds to the partition function V^N
of the ideal hard disks gas of radius zero, that Alberto explained.
The next term is proportional to 1/V
In fact, we obtain this term in 1/V because we consider only the term with one Upsilon in the expansion.
We could go on in this expansion of the brackets and consider pairs of Upsilon, triplets and so on
This would yield an expansion - as previously - in powers of 1/V
However, this is not the state of the art of statistical mechanics for the 21st century.
Indeed, Ludwig Boltzmann in 1874 computed the 4th term of the virial expansion,
that is to say the term in 1/V^3.
It was for three-dimensional hard spheres.
This computation was really complicated,
it was subject to years of debate, but in the end he was right and others were not.
The virial expansion was once believed to provide a systematic access to the thermodynamics of liquids and gases, in our case hard disks.
In the same as the expansion of the exponential function
exp(x) = 1 + x + x^2/2 + x^3/6 and so on.
This expansion is valid for all values of x, real or complex.
However, for the virial expansion this is not the case any more.
What is proven today is that this expansion is valid
only up to some finite values of the density eta.
However it cannot predict the existence of a phase transition.
Those phase transitions will be the subject of our next week.
Now, please, download run and modify the programs we have seen during this session.
You may also have a look to the fact-sheet that corresponds to the computation
of the virial expansion we have just seen.
On the Coursera website you may also have a look on the following program.
It is sweet and short,
it allows to determine the acceptance ratio for any density eta.
Have fun going through this
those are 15 lines, but not easy to understand.
In conclusion, we have moved in this tutorial
from the tabula rasa rule to the discussion of the high-dimensional abstract spaces of our hard disk configurations.
We moved towards the center of statistical physics: the partition function and the virial expansion.
Our discussion revealed a beautiful connection
between algorithms and the physical description of our systems.
But our analytical approach remained perturbative.
This is a true limitation
We will see next week that our algorithms have the potential to carry us beyond this horizon and this limitation of the perturbative approach
This will allow us to study phase transitions,
when the liquid of hard disks becomes solid,
and this is a major revolution in physics.
In the meantime, have fun with the homework session of this week
and see you again next week at Statistical Mechanics: Algorithms and Computations.