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.