1 00:00:03,220 --> 00:00:08,960 Hello, Bonjour, kumusta ka, sawa dee-krap! 2 00:00:09,040 --> 00:00:18,900 Welcome to the eighth week of Statistical Mechanics: Algorithms and Computations from the Physics Department of Ecole normale supérieure. 3 00:00:18,900 --> 00:00:28,660 This week and the next one, our subject will be the statistical mechanics and computational physics of the Ising model, 4 00:00:28,660 --> 00:00:33,140 That has inspired generations of physicists 5 00:00:33,140 --> 00:00:41,700 This very simple model of spins on a lattice with simple pair interactions 6 00:00:42,160 --> 00:00:48,640 Undergoes an order/disorder transition in dimensions 2 and higher 7 00:00:48,640 --> 00:00:57,660 and, in two dimensions, it can be solved exactly, as first shown by Onsager in 1944. 8 00:00:57,660 --> 00:01:05,520 In this week's lecture, after a short introduction, we will treat Monte-Carlo algorithms for the Ising Monte-Carlo. 9 00:01:05,520 --> 00:01:13,320 First, a local metropolis algorithm, and then, a global cluster algorithm. 10 00:01:13,320 --> 00:01:22,960 Cluster Monte-Carlo algorithms originated here in the Ising Model, and they have since revolutionized computations 11 00:01:22,960 --> 00:01:27,300 in many fields of classical and quantum physics. 12 00:01:28,000 --> 00:01:36,640 The algorithm we present here, fortunately for us, can be explained in a few minutes, 13 00:01:36,640 --> 00:01:42,060 and implemented in just over a dozen lines of Python code. 14 00:01:42,060 --> 00:01:51,580 In the tutorial, we will consider the Heatbath algorithm, and illustrate one of its surprising features, 15 00:01:51,580 --> 00:01:59,080 namely that it couples and that it allows for perfect simulations. 16 00:01:59,080 --> 00:02:08,300 The Ising model is a unique meeting point for Mathematics, Computer Science, and many branches of Physics 17 00:02:08,300 --> 00:02:15,900 It as a less immediate connection with classical mechanics than the Hard disks model 18 00:02:16,120 --> 00:02:22,200 because there is no molecular dynamics, and no kinetic energy, 19 00:02:22,580 --> 00:02:28,960 but the Ising model phase transition is much better understood, 20 00:02:28,960 --> 00:02:36,680 and in two dimensions, there are exact solutions, even for finite systems, 21 00:02:36,680 --> 00:02:44,400 so this allows for very fine-tests of the algorithms on lattices of any size. 22 00:02:44,400 --> 00:02:52,400 So, let's get started, with the Statistical mechanics, and computational physics of the Ising model, 23 00:02:52,400 --> 00:03:02,400 and with week 8 of Statistical Mechanics: algorithms and computations. 24 00:03:02,860 --> 00:03:12,400 The Ising model describes classical spins sigma that can be +1 or -1 on a lattice 25 00:03:12,400 --> 00:03:16,980 like the two-dimentionnal square lattice shown here. 26 00:03:16,980 --> 00:03:26,980 We consider the ferromagnetic Ising model, where neighbooring spins prefer to be aligned. 27 00:03:27,640 --> 00:03:36,980 Two neighbooring up spins or two neighbooring down spins have an energy -1, 28 00:03:36,980 --> 00:03:44,700 whereas two spins like this or like that have energy +1. 29 00:03:44,700 --> 00:03:54,700 This gives an energy of each configuration like this, where the sum is over all pairs. 30 00:03:54,700 --> 00:04:02,020 But attention: we count each pair only once. 31 00:04:02,020 --> 00:04:12,180 For concreteness, here is a list of the 16 configurations of the 2 by 2 Ising model in two dimensions. 32 00:04:12,920 --> 00:04:24,120 Without periodic boundary conditions, these two configurations have energy +4, these two have an energy -4, 33 00:04:24,540 --> 00:04:29,740 and all other configurations have energy 0. 34 00:04:30,040 --> 00:04:38,280 For small systems, we can enumerate all configurations explicitely, 35 00:04:38,280 --> 00:04:46,980 and a nice little algorithm implementing Gray code enumeration is shown here. 36 00:04:46,980 --> 00:04:53,320 Here is the output of the movie version of the movie version of enumerate_ising.py, 37 00:04:53,320 --> 00:04:58,100 for 4*4 Ising Model for periodic boundary conditions. 38 00:04:58,100 --> 00:05:05,360 Each configuration differs from the previous one on a single site only, 39 00:05:05,360 --> 00:05:13,800 and we don't have to recalculate the energy from scratch at each step. 40 00:05:13,800 --> 00:05:23,800 It is sufficient to compute the field H on the spin that is going to flip 41 00:05:23,800 --> 00:05:33,800 and the new energy is then equal to the old energy plus 2*H*sigma. 42 00:05:33,800 --> 00:05:43,800 Keeping track of the energy like this, we can compute the partition function as a function of the temperature, 43 00:05:43,800 --> 00:05:51,580 the mean energy, the mean square energy, and the specific heat: 44 00:05:51,580 --> 00:05:56,320 the change of the mean energy with temperature. 45 00:05:56,360 --> 00:06:08,580 There are 2 things to notice: first of all, the specific heat, altough it is a derivative with temperature of the energy, 46 00:06:08,580 --> 00:06:13,320 does not have to be computed by numerical derivation. 47 00:06:13,320 --> 00:06:23,320 It is simply given by a quantity that is proportional to the variance of the energy. 48 00:06:24,680 --> 00:06:35,660 Secondly, to compute the thermodynamic quantities, it is not necessary to do the enumeration 49 00:06:35,660 --> 00:06:46,000 at all temperatures. It's a lot better to compute the number of configurations per energy, 50 00:06:46,020 --> 00:06:48,420 the density of states. 51 00:06:49,300 --> 00:06:58,860 So, compute the density of states once and for all, like here in enumerate_ising.py, 52 00:06:59,580 --> 00:07:10,680 and then use this density of states for all temperatures to compute the thermodynamics, as in in thermo_ising.py. 53 00:07:11,640 --> 00:07:18,580 So, here, you see the specific heat for small systems, 54 00:07:19,580 --> 00:07:28,580 and the curve already indicates the phase transition that is known to take place 55 00:07:28,580 --> 00:07:36,080 at a temperature of 2/log(1+sqrt(2)). 56 00:07:37,120 --> 00:07:47,760 We can do our enumerations up to system sizes of 4*4, or, with a little bit of effort, to 6*6. 57 00:07:48,620 --> 00:07:57,760 But our pure Python program is a bit to slow for real number crunshing 58 00:07:57,760 --> 00:08:04,580 so, let's go ahead, and consider the sampling approach. 59 00:08:04,580 --> 00:08:15,240 Before doing so, however, please take a moment to download, run and modify the programs that we discussed so far. 60 00:08:15,240 --> 00:08:25,240 There is this nice little program enumerate_ising.py, that uses the Grey code enumeration 61 00:08:25,240 --> 00:08:35,240 for all the configurations. Try to understand what the Grey code does, you may need it later. 62 00:08:39,460 --> 00:08:45,240 In the Ising model, we can get very far by counting configurations 63 00:08:45,240 --> 00:08:50,140 even tough the listing of them becomes very difficult. 64 00:08:50,740 --> 00:09:01,940 But in general, for large systems and systems that don't exactly correspond to the Ising energy - sigma_i sigma_j, 65 00:09:01,940 --> 00:09:05,740 the sampling approach is more reasonable. 66 00:09:06,600 --> 00:09:19,020 Analogously to what we did for hard disks, we select one site and attempt to flip that spin on that site. 67 00:09:19,980 --> 00:09:31,060 The move from configuration a to the flipped configuration b must be accepted with the Metropolis acceptance probability... 68 00:09:38,860 --> 00:09:48,460 This algorithm in Python is the following: markov_ising.py. 69 00:09:48,460 --> 00:09:58,460 It is even simpler than the corresponding version for hard disks: markov_disks.py 70 00:09:58,460 --> 00:10:07,200 Output of the graphics version of markov_ising.py is shown here. 71 00:10:07,200 --> 00:10:17,840 This program can be checked against the exact enumeration for small lattices, but I can assure you that everything comes out all right. 72 00:10:18,580 --> 00:10:30,020 This program also very easily recovers the phase transition that takes place between the ferromagnet and the paramagnetic phase at a temperature just above 2. 73 00:10:30,440 --> 00:10:38,860 You can trace this phase transition by plotting the absolute mean magnetization as a function of the temperature, 74 00:10:38,860 --> 00:10:48,860 even though there are more expert ways of locating Tc that you will discover in this week's homework session. 75 00:10:48,860 --> 00:11:00,840 You can also simply look at the configurations as a function of the temperature, and you see a very clear qualitative difference between 76 00:11:00,840 --> 00:11:08,720 the paramagnetic regime at high temperature, and the ferromagnetic regime a lower temperatures. 77 00:11:08,720 --> 00:11:16,000 Close to the critical temperature, the local Monte-Carlo algorithm slows down increasingly. 78 00:11:16,000 --> 00:11:23,260 In a nutshell, this is due to the increase of the correlation length in the system, 79 00:11:23,260 --> 00:11:33,260 and to the fact that turning around a spin in a larger correlated region becomes increasingly difficult. 80 00:11:33,260 --> 00:11:45,860 The large-scale correlations of the local Monte-Carlo algorithm close to the phase transition is also known in experiments 81 00:11:46,060 --> 00:11:58,700 and it is called critical slowing down. Critical slowing down is overcomed by the cluster algorithms that we next discuss. 82 00:11:59,400 --> 00:12:11,080 But vefore going to the next session, please take a moment to download, run ad to modify the programs that we just discussed: 83 00:12:11,080 --> 00:12:18,980 this is the algorithm markov-ising.py, with a local Markov-Chain Monte-Carlo algorithm, 84 00:12:18,980 --> 00:12:25,980 and its graphics version, markov_ising_movie.py. 85 00:12:30,400 --> 00:12:42,620 The traditional simulations methods for the Ising model, markov_ising.py and heatbath_ising.py that we will discuss in this week's tutorial, 86 00:12:42,720 --> 00:12:50,140 have gradually given way to cluster algorithms that converge much faster. 87 00:12:50,140 --> 00:12:56,800 Rather than flipping a spin indiviually at each timesteps, 88 00:12:56,800 --> 00:13:07,660 these algorithms construct larger ensembles, clusters, of the length scale of the system, and flip them in one step. 89 00:13:08,620 --> 00:13:20,560 Cluster algorithms are truly elegant, very general, and you will see, they turn out to be very easily programmed. 90 00:13:20,780 --> 00:13:30,560 Suppose that, starting from a random initial site, rather than flipping the site on that site, 91 00:13:31,460 --> 00:13:42,060 we attempted to construct a cluster of connected up spins and flip them later. 92 00:13:43,060 --> 00:13:50,340 Clearly, we cannot construct a cluster from all the connected up-spins, 93 00:13:50,340 --> 00:14:00,340 flip then, then construct a cluster of all the connected down spins, and flip them again. 94 00:14:00,340 --> 00:14:10,340 Very soon, we would be flipping between the all-up and the all-down configuration 95 00:14:10,340 --> 00:14:17,020 and the Boltzmann probabilities would be violated. 96 00:14:19,000 --> 00:14:29,860 So, from a spin already in a cluster, we should accept a spin outside of the cluster with probability p: 97 00:14:29,860 --> 00:14:35,080 sometimes we grow the cluster and sometimes we dont. 98 00:14:35,080 --> 00:14:45,080 So here, in this example, we added all three new spins to the site already present 99 00:14:45,080 --> 00:14:48,920 all of them with probability p. 100 00:14:49,980 --> 00:15:00,880 We added them to the cluster and also to a pocket, just to remember that we haven't checked their neighboors yet. 101 00:15:00,880 --> 00:15:10,880 At some moment, the pocket is empty, for example for the cluster shown here. 102 00:15:10,880 --> 00:15:23,840 You see, and can count for yourself, that accross the boundaries of the cluster, there are 18 links plus-minus 103 00:15:23,840 --> 00:15:28,560 and 14 links plus-plus. 104 00:15:30,020 --> 00:15:39,400 For each of these 14 links, the extension of the cluster was rejected 105 00:15:40,180 --> 00:15:49,400 because at some time, a random number was drawn and was larger than p. 106 00:15:50,820 --> 00:16:03,360 So then, we can flip this cluster and arrive from the configuration a to the configuration b. 107 00:16:03,360 --> 00:16:16,060 Notice that in the configuration b, there are now 18 links minus-minus, 108 00:16:16,060 --> 00:16:26,360 and 14 links minus-plus. Now you must take into account the detailed balance condition 109 00:16:27,100 --> 00:16:34,460 The statistical weight of configuration b times the probability to move from a to b 110 00:16:34,460 --> 00:16:43,420 must be equal to the statistical weight of configuration b times the probability to move from b to a. 111 00:16:44,420 --> 00:16:56,440 The probability to move from a to b is given by the probability to construct this cluster 112 00:16:57,060 --> 00:17:09,440 and then to stop the construction of the cluster at this stage, times the probability to accept the move from a to b 113 00:17:09,440 --> 00:17:22,640 Likewise, the probability to move from b to a is given by the probability of constructing this cluster 114 00:17:22,640 --> 00:17:33,080 the flipped cluster rather than another cluster, to stop the construction of this cluster at that stage, 115 00:17:33,080 --> 00:17:38,580 times the probability of accepting the move from b to a. 116 00:17:38,580 --> 00:17:48,580 Now, for the a priori probability A(a->b) to construct this cluster rather than another cluster 117 00:17:48,580 --> 00:17:52,340 and to stop the construction of the cluster right here: 118 00:17:52,340 --> 00:18:02,700 In these a-priori probabilities, there are terms that are the same in a and b, and they drop out of the picture, for example 119 00:18:02,700 --> 00:18:12,700 the probability to add a plus spin to a plus spin already present in a is the same as the probability 120 00:18:12,700 --> 00:18:18,720 of adding a minus spin in b to a spin already there. 121 00:18:18,720 --> 00:18:29,800 The game is really played at the boundary: the probability to stop the cluster construction here in a 122 00:18:29,800 --> 00:18:35,900 differs from the probability to stop cluster construction in b. 123 00:18:35,900 --> 00:18:44,520 So, each site on the boundary of the cluster in a was once a pocket site, 124 00:18:44,520 --> 00:18:57,500 and the construction of this cluster came to a halt because none of the edges accross the cluster was accepted, 125 00:18:57,500 --> 00:19:06,740 precisely, in this cluster a, they are 14 edges plus-plus 126 00:19:06,740 --> 00:19:13,900 and the extension was rejected with a probability 1-p for each of them. 127 00:19:13,900 --> 00:19:26,020 So the a-priori probability from a to b is proportional to (1-p)^14. 128 00:19:27,360 --> 00:19:36,660 Analogously, for the cluster in b, also the construction came to a halt 129 00:19:36,660 --> 00:19:44,340 because none of the possible minus-minus edges accross the boundaries were accepted. 130 00:19:44,340 --> 00:19:55,620 There are 18 minus-minus edges in b, so the a-priori probability from b to a 131 00:19:55,620 --> 00:20:02,480 is proportional to (1-p)^18. 132 00:20:02,480 --> 00:20:10,040 Now, let us consider the statistical weights of configurations a and b: 133 00:20:10,040 --> 00:20:17,280 again, there are terms outside the cluster, but they are the same in a and b. 134 00:20:17,280 --> 00:20:27,660 There are terms inside the cluster, also they are the same in a and b. The game is played at the boundary. 135 00:20:27,660 --> 00:20:40,560 Now, in configuration a, threy are 18 terms plus-minus, hence 18 terms +1 136 00:20:40,560 --> 00:20:45,440 and 14 terms -1 in the energy, 137 00:20:45,440 --> 00:20:55,440 or if we write N1 instead of +1 and N2 instead of the terms plus-plus, 138 00:20:55,440 --> 00:21:04,320 the energy accross the boundary is N1-N2 in a, 139 00:21:04,320 --> 00:21:09,100 and it is N2-N1 in b. 140 00:21:09,100 --> 00:21:16,000 Now, we have all that it takes to write down the detailed balance condition, 141 00:21:16,000 --> 00:21:26,200 like here, we have the statistical weight of configuration a, times the a-priori probability to construct this cluster 142 00:21:26,200 --> 00:21:30,600 times the acceptance probability to move from a to be, 143 00:21:30,600 --> 00:21:40,600 equals the statistical weight times the construction probability times the acceptance probability from be to a. 144 00:21:40,600 --> 00:21:46,160 This gives the Metropolis acceptance rate shown here, 145 00:21:47,720 --> 00:21:58,780 and you see we can simplify this acceptance probability and write it as... 146 00:22:14,760 --> 00:22:21,780 You see that this equation again simplifies for the magic value of... 147 00:22:28,020 --> 00:22:33,620 For this magic value, we have no rejections. 148 00:22:34,220 --> 00:22:46,380 So with this value of p, we may simply construct the cluster, then start another one, construct the cluster and flip it, 149 00:22:47,300 --> 00:22:56,380 another construction flip, and so on. This is the famous cluster algorithm 150 00:22:56,380 --> 00:23:06,380 by Wolff from 1989. In Python this gives the following 19-lines program: 151 00:23:06,380 --> 00:23:15,700 so in this algorithm, you take a random spin and you put it into the custer and into the pocket, 152 00:23:15,820 --> 00:23:25,080 then, a complete cycle of this algorithm consists in taking a random element out of the pocket, 153 00:23:25,080 --> 00:23:35,880 and checking all of its neighboors. For each of the neighboors, you check wether it is not already in the cluster, 154 00:23:36,540 --> 00:23:46,680 wether it has the same sign, and wether a random number between 0 and 1 is smaller than p. 155 00:23:46,680 --> 00:23:59,000 If these 3 conditions are satisfied, you add the new spin to the cluster and to the pocket. 156 00:23:59,000 --> 00:24:08,440 Once you have checked all the neighboors of this spin, you take it out of the pocket. 157 00:24:08,440 --> 00:24:16,560 Once the pocket is empty, you flip the entire cluster 158 00:24:17,580 --> 00:24:26,560 This is all there is to one of the most influencial algorithms in all statistical physics. 159 00:24:26,560 --> 00:24:32,020 An example of how this cluster algorithm works is shown here. 160 00:24:32,020 --> 00:24:40,960 Close to the transition temperature of the Ising model, enormous clusters are constructed, 161 00:24:40,960 --> 00:24:45,400 and they are flipped without thinking twice. 162 00:24:46,740 --> 00:24:57,420 It moves though configuration space with breathtaking speed, and far outpaces the local Monte-Carlo algorithm. 163 00:24:58,660 --> 00:25:09,060 Very important it is to realize that the speedup it realizes with respect to markov_ising.py 164 00:25:09,060 --> 00:25:21,280 increases with system size, so for a small system it may be 10 times faster, for a larger system 100, 1000, 1000000... times faster 165 00:25:21,280 --> 00:25:32,200 than the local algorithm. So before continuing, please take a moment to go over a derivation of this cluster algorithm 166 00:25:32,200 --> 00:25:43,940 and then to download, run and modify cluster_ising.py, to see how these simple ideas are implemented in Python. 167 00:25:46,800 --> 00:25:56,620 In conclusion, we have introduced in this lecture to the Ising model and its associated Monte-Carlo algorithm. 168 00:25:57,100 --> 00:26:09,060 What we have presented here, and what we can program in just over a dozen lines today, has taken several decades to develop. 169 00:26:09,900 --> 00:26:18,240 Roughly speaking, in Ising model computations, there have been 2 crucial periods: 170 00:26:18,240 --> 00:26:29,880 The first crucial period was the following: once it was realized that Monte-Carlo calculations were rather imprecise 171 00:26:29,880 --> 00:26:40,140 in the viscinity of the critical point, it has taken a long time to understand that there were two reasons: 172 00:26:40,140 --> 00:26:52,280 one was that it is quite complicated to compute observables close to the critical point, because of the critical slowing down as we discussed. 173 00:26:53,260 --> 00:27:01,400 But there is a second reason, and the second reason is that the observables of the infinite systems 174 00:27:01,400 --> 00:27:07,640 are quite different from the observables in the finite system. 175 00:27:07,640 --> 00:27:18,140 The extrapolation from the finite system where we can do our calculations to the infinite system which we are really interested in, is halas non trivial. 176 00:27:19,660 --> 00:27:29,800 The second crucial period of Monte-Carlo calculations came in the 1980, 177 00:27:30,600 --> 00:27:42,960 when it was found out that critical slowing down was not an inescapable faith for algorithms as it is for experiments, 178 00:27:44,280 --> 00:27:55,080 and this was realized when the first cluster algorthms came out. These cluster algorithms realized non-physical dynamics 179 00:27:55,080 --> 00:28:03,620 but they satisfied detailed balance, so there final state is exactly the same as the one of the local algorithm. 180 00:28:03,620 --> 00:28:14,880 As we have seen, they do through configurations at breathtaking speed without suffering from the critical slowing down. 181 00:28:14,880 --> 00:28:23,620 We will immerse ourselves farther in this fascinating subject in this week's tutorial and homework session 182 00:28:23,620 --> 00:28:29,140 but for the time being, let me thank you once again for your attention. 183 00:28:29,140 --> 00:28:41,100 See you again next week, in week 9 in Statistical Mechanics: Algorithms and Computations, for faster-than-the-clock algorithms. 184 00:28:41,100 --> 00:28:51,100 So good luck with this week's homework session, homework session 8.