Hello, bonjour, sabah el kheir, pozdrav,
welcome to the third tutorial of Statistical Mechanics: Algorithms and Computations,
from the Physics Department of Ecole Normale Superieure.
This week's subject is a model of random clothes-pins
randomly positioned on a line between two poles, our boundaries.
This models is equivalent to one-dimensional hard spheres.
This model allows us to move forward
very far, along three directions
First: it is prime example for the Asakura-Oosawa depletion interaction
which is of great importance for macromolecules, polymers, blood cells and other systems.
This is what we studied in the lecture using the halo picture.
Second: the random clothes-pin model can be simulated very efficiently
we started out in the lecture with a naive algorithm
but Michael will show a providential method
that allows us to simulate hundred, thousand and millions of clothes-pin on a line
without rejection and without any effort.
Third: this fantastic algorithm will show us the way to solve the model analytically.
Alberto will show us how to compute the partition function,
how to sum up the virial expansion
and how to prove that this system in fact does not have a phase transition.
Then Vivien will investigate
the relationship between pins on a line with walls or poles
and pins with periodic boundary conditions on a ring.
He will characterize the correlations in the system
and show that they are short range and that this characterizes a liquid.
In what we will do, we will see that analytical results are strong
and the algorithms that Michael will present are superb.
But don't worry, the Python implementation will be about 5 to 15 lines long.
In this week's lecture, Werner presented the random clothes-pin model,
a model with hard-core interactions, just like the hard-disks we discussed earlier.
We consider N clothes-pins of width (2 sigma)
on a washing line, that is a line segment,
between boundaries at x=0 to x=L
as is shown here.
The pins are placed one after the other
but if an overlap is generated, we take them all off the line and start anew.
This is again our familiar tabula rasa principle
and once more it gives us a flat probability distribution.
It samples N-pins configurations
with a distribution pi(x_0,..,x_(N-1))
=1 if the configuration is legal
and =0 if the configuration is illegal.
This is in agreement with the equiprobability principle:
all legal configurations have the same energy
E=0.
Our pins are in one dimension what hard disks are in two dimensions and hard spheres are in three dimensions.
A configuration is illegal if two pins overlap.
This happens if their centers are less than (2 sigma) away from each other.
Likewise, pins are not allowed to overlap with the boundaries.
This happens if their centers are away less than sigma from x=0 or x=L.
Here is how the program can be implemented in Python
Of course this algorithm works and obeys the equiprobability principle
as it should.
but - just like in the case of two-dimensional hard-disks - it has a very high rejection rate.
Output of the program is shown here.
Is this the best algorithm for this problem?
Or can we be a bit smarter? Of course we can.
In fact, we are able to sample N-pins configurations without producing any illegal configuration on the way.
This gives us a rejection rate of exactly zero.
What you will see now is an example of a rejection free algorithm
and it is way more sophisticated of all the algorithms we have presented so far in this course.
But how can we possibly achieve a zero rejection rate?
This works because we already learned from our experience with hard-disks
that in the full configuration space
of legal and illegal configurations
the regions of illegal configurations are like the holes in the Swiss cheese.
For the clothes-pin problem, we can actually remove the holes from the cheese.
To illustrate the idea, let us consider this example
of 5 clothes-pins on a washing line.
The first pin could not be placed within the halo of the left wall.
Likewise, the second pin cannot be placed
within a distance of (2 sigma) from the first pin's center,
and so on, for the other pins.
Finally, no pin can be placed within the halo of the right wall.
We can therefore deflate the washing line by cutting out these locked regions.
If we focus on the pin centers, and glue the remaining line together,
we obtain a shorter washing line of length (L - 2 N sigma).
Let us call the pin centers on this shorter washing line y_0, y_1, ...
Then these y's are related to the original positions x_0, x_1, x_2, ...
on the longer washing line
by x_i = y_i + 2 (i+1) sigma.
Since every legal configuration can be deflated in the same way and gives a shorter washing line
of the same length,
we realize that the legal configurations
distribute the pin centers not on the long washing line we realize that the legal configurations
distribute the pin centers not on the long washing line
but on a shorter washing line of length L - 2 N sigma.
This becomes even clearer if we invert this procedure.
Any chosen five points on the shorter washing line, for example this five points here,
can be inflated back to a configuration of five pins on a longer washing line of length L.
Let us turn this into a Python program.
Start by randomly choosing N points between 0 and L - 2 N sigma
then sort them so that y_0 < y_1 < y_2 < ..
and then finally take care of the inflation in one step.
In Python, this gives the following program.
First we choose the y's randomly between 0 and L - 2 N sigma,
then we sort them, and this line here
takes care of the inflation,
that is: it related the positions y_0, y_1, y_2, ... to the positions x_0, x_1, x_2, ...
Notice that the program is called noreject
in honor of its remarkable rejection rate of exactly zero.
Finally you can see here the comparison of outputs of the programs direct_pins.py
and direct_pins_noreject.py, for 10 pins
at a density of 75%.
Visibly they are the same, but direct_pins.py
generated more than 500000 illegal and therefore useless configurations
on the way, while every single configuration generated by direct_pins_noreject.py is legal.
Before moving on, please take a moment or two to download, run and modify the two algorithms discussed in this section.
On the Coursera website there is the algorithm direct_pins.py
which mirrors last week's direct_disks_box.py
There's also the algorithm that Alberto will discuss in a few moments
and this algorithms is really beautiful: it allows you to sample configurations
with thousand, millions, billions clothe-pins without any rejection.
Unfortunately, it is very difficult to generalize this algorithm to 2 or even 3 dimensions
mainly because the crucial sorting steps then becomes impossible.
Michael just presented a rejection-free direct-sampling Monte Carlo algorithm
for a many-body interacting system.
This algorithm is rooted in the fact
that we can actually compute the exact partition function of the system
namely the number of legal configurations
given by the following integral,
where pi as usual is 0 if there is an overlap and 1 if there is no overlap.
In Michael's derivation, the sorting step played a crucial role.
We can do the sorting step in the above integral
and obtain the following result
where the theta function selects between the N factorial possible orderings of the pins' centers
the sorted one with x_0 < x_1 < x_2 < ..
We can now apply Michael's deflation step
and change the x variables to the y variables.
We obtain this integral, where there is no longer overlap condition
but only the ordering of the y variables imposed by the theta function.
Remember: the theta function picks out one of the N factorial permutations of the integrand
so we can remove both the theta function and the N factorial,
compute the integral and get the final result for the partition function.
This is a beautiful formula.
Which allows you to compute the acceptance probability
for the program direct_pin.py, that was introduced in this week's lecture.
The acceptance probability is the ratio of the volume of legal configurations to the volume of legal and illegal configurations.
To test your understanding on this concept, please let me ask you a few questions
Question 1: select from the following three possibilities
the correct expression for the acceptance probability of direct_pins.py.
The correct answer is option A, because in direct_pin.py
the centers are generated between sigma and L - sigma
and because the volume of legal configurations
is (L - 2 N sigma)^N
(the partition function that we have just computed).
Question 2: select from the following three possibilities
and estimate of the average number of rejected configurations
that should be generated before having a single legal configuration.
The correct answer is option B,
because the number of configurations that we have to wipe-out before generating a single legal configuration is on average the inverse of the acceptance probability
Question 3: compute this number for a density of eta = 0.75 and 15 pins
Is this number 5000000,
50000000 or 500000000?
As anticipated by Werner in the lecture, the correct answer is option C.
The acceptance probability is 2 times 10^(-9)
which means 500000000 rejected configurations.
Now, let us recall the virial expansion introduced by Vivien in Tutorial 2.
For clothes-pins, we can actually compute exactly the derivative of log(Z)
with respect to the volume, or rather with respect to the length L of the system.
And we obtain this exact expression
that can be written in an even simpler way 1- 1/eta, where eta is the pin density.
Note that this quantity can be expressed as an infinite virial expansion
1 + eta - eta^2 + eta^3 ...
This series converges for all densities smaller than one
this means that there is no phase transition.
However, we not the divergence at the point eta=1
a special point where we have long-range order.
Knowing the partition function, we can actually compute exactly other quantities
like the distribution pi(x) that shows the oscillatory behavior that you have seen in this figure.
In fact, pi(x) is the probability of having a pin center at x
which is also the probability, given the pin center at x, to put the other pins
on its left and its right.
Let's try to place k pins on the left, this is possible only if we have enough space.
The number of possibilities of doing so is the partition function of k pins in the interval x - sigma.
In the same way, we have to put N - 1 - k pins on the right.
Putting these pieces together gives the following formula,
where a combinatorial factor gives the number of choices of picking
k left pins out of the bucket of the N - 1 pins
that remain once we have placed the first one at x.
In Python, this gives the following program
If you have doubts on our calculation, you can compare this result with the result that I have just shown from the numerical simulation.
The function pi(x) increases steeply close to the poles
Our formula allows us to see that the distribution pi(x) for x=sigma
is a factor 1/(1 - eta) larger than in the center.
This becomes very large at high densities, where the pins are strongly attracted by the walls.
In a few moments, Vivien will show us that the attraction of the wall is not a simple boundary effect:
the same attraction acts also between two clothes-pins.
Before doing that, please take a moment
to download, run and modify the program we discussed in this section.
On the Coursera website you will find direct_pins_density.py
Its output should give you the same results as the numerical simulations performed with
direct_pins_noreject.py, discussed by Michael.
Several times throughout the lectures and the tutorials we have generated an annoying situation
by putting particles in a box rather than allowing for periodic boundary conditions.
For instance, in the adults' game in the heliport,
we could have allowed a particle which exits the system from the right
to enter back the system from the left.
Nevertheless, such strict but realistic rule of the wall has allowed us to learn about the Metropolis algorithm.
That is, how to treat rejections.
The same objection can be made for our system of interest: the pins on a line.
If the probability distribution pi(x) is oscillatory close to the boundaries, this must be due to the poles.
Now it is true that the probability depends on the position
because of the poles, because of the boundaries.
This raises the following question
with periodic boundary conditions, this probability is uniform?
So, does it mean that our discussion above was useless?
Not at all, as we now discuss.
by examining and comparing the two situations
where the pins are put on a line
to the situation where the pins are put on a circle
that is to say: with periodic boundary conditions.
Let us first write a program to simulate our periodic system
we start by putting a particle or a pin
at position L - sigma.
Its right side provides a pole at position x=0
due to the periodic boundary conditions
this leaves us with N - 1 pins
that we have to put between 0 and L - 2 sigma
where L - 2 sigma is the left side of the pin we have put initially.
All in all, this means that our system N pins on a circle of length L
is completely equivalent to a system of N - 1 pins that we put on a line of length L - 2 sigma.
Besides, to ensure the invariance by translations, we shift the centers of the pins
by the same number which is drawn uniformly between 0 and L.
In Python, this gives the following program
direct_pins_noreject_periodic.py
By running this program, we observe that the distribution pi(x) converges to a flat function
as the number of samples increases from 1000 to 10000 to 100000.
You notice a sharp difference compared to the previous boundary conditions:
the oscillations of the distribution close to the poles have vanished.
But beware, this does not mean that there are no correlations in the system between the pins.
To study such possible correlations, a good tool is provided by the pair correlation function
pi(x,x'), which represents the probability
of having a pins at position x and another pin at position x'.
This correlation function is a useful tool.
For instance it is trivially equal to a constant
if there are no correlations in the system.
Now let us analyze the pair correlation function using the following program
direct_pins_noreject_periodic_pair.py
Let us fix the density to eta=0.9
and increase the system size at fixed pin width 2 sigma = 1.
We observe that the pair correlation function quickly converges to its thermodynamic limit
in the example displayed here
for system size equal 100, 500, 1000 and 2000.
You notice that this pair correlation function presents oscillations in the same way as for the distribution
for the system on a line.
In fact, those two functions are the same
because the pin which is at position x acts as a pole for the pin at position x'.
Indeed, remember now how we devised our simulation
we first put a pin at position L - sigma
using the same reasoning you will convince yourself that the following formula holds
it relates the pair correlation function
for a system of N pins on a circle of length L
to the probability distribution pi for N - 1 pins
in a line of length L - 2 sigma.
Besides, Alberto explained to you previously
that this function can be analytically determined
Hence, to summarize, we have constructed a system which is exactly solvable, in which you can compute
in an analytical way the pair correlation function
in one dimension and for a system where there are non-trivial interactions between pins.
Let us now physically interpret the features of this correlation function.
As you see here, the oscillations are embedded in a decaying envelope
and the correlation goes to a constant at large distances.
To be more quantitative, let us plot this decay in a logarithmic scale.
As you now observe, the linear behavior indicated that the correlation function
decays exponentially as exp(- x / xi)
xi has a physical meaning
it represents the correlation lengths between the pins in the system.
At pin density eta=0.8, its value is approximately 2
Increasing to a density 0.9, we observe the same behavior
now with xi=6
And a dramatic increase to xi = 26 when eta=0.95.
In fact, one can show that the correlation length diverges as the density eta goes to one.
Alberto has shown that the partition function of our system of pins in one dimensions
can be computed and is analytic for all densities, which means that there are no phase transitions.
This is now supplemented by the fact
that the pair correlation function, which can also be computed,
decays exponentially on a length xi
This length xi, which is the length of the correlations between the pins,
remains finite even if the system size goes to infinite.
At distances larger than xi, the system is homogeneous
This is what defines a liquid state.
Long-range positional order is possible in our system only at density 1.
that is for L = 2 N sigma
that is for configurations where the pins are close-packed.
What we have shown conclude our study of this one-dimensional model
but in fact what we have found
(the exponential decay of correlation functions and the absence of phase transition)
is in fact generic in one-dimensional systems with local interactions.
Naturally many subjects remain to be studied, especially if one wants to carry out
the insights we arrived to from dimension 1 (realizing wires and tubes)
to dimension 2 (surfaces, interfaces) and even to dimension 3
(bulk materials).
In this tutorial, we have gotten quite far in our study of the random clothes-pin model
in other words: of one dimensional hard spheres.
We again obtain a beautiful relationship
between the algorithms and the analytical solutions.
The providential direct sampling Monte Carlo algorithm
showed us the way to the analytical solution of the model,
the calculation of the partition function, of the pair correlation function and of density distribution.
From a physical point of view, we saw that the system of random clothes-pin is a liquid
There are some correlations close to the position of one pin, but these correlations decay exponentially.
On large length-scales, the system is homogeneous
So again, at small distances we have structure, like in a crystal
but on large distances it is homogeneous, like in water.
So this system is a liquid.
This week and the previous week we studied hard sphere models
and maybe you wander why we never mentioned velocities or the temperature.
This was possible in fact because the systems we study are athermal
they are the same at all temperatures
Wait until next week, when we will discuss the Maxwell and Boltzmann distribution
and take up what was missing today and last week:
the discussion of velocities, of temperature and of finite interactions.
Until then, have fun with the homework session 3 of Statistical Mechanics: Algorithms and Computations.