Skip to main content
Loading Explore guides…
Monte Carlo Sampling · Acedra
Home
Explore
Physics
Monte Carlo Sampling
Sign in to save
Chat
If YouTube asks you to sign in, open the video on YouTube.
Watch on YouTube
Reload player
Transcript
83 segments
Copy
Hide video
0:24
Hello all of you, so, in the last lecture we have outlined the basic idea behind the Monte Carlo simulation but I have not really said how exactly we do that this will be the part that we discuss now. So, the essential idea and I am again redrawing the phase space
0:48
that the essential idea here is that we want to sample the points in the phase space
1:00
unlike the molecular dynamic simulation where we follow the trajectory starting from a phase point.
1:09
We take for our advantage the fact that the probability density actually goes to a stationary solution in the phase space and therefore instead of like following the trajectory, we can simply sample the phase space and as long as I am sampling according to this particular probability density, I am essentially getting the properties that I am interested in at equilibrium condition.
1:49
So, we need to find a mechanism in which we can sample this particular phase space. So, this is like I would say the physical idea behind it but, let us now try to understand what is the mathematical problem that we are trying to solve? So, I have said that, let us say if I am working in the canonical ensemble, we
2:19
will develop the ideas in the canonical ensemble. But we can apply Monte Carlo in other ensembles as well but I am developing in the NVT ensemble there are some details that need to be considered when we are going from 1 ensemble to the other. But I will not get into that for a moment, so I am working in the canonical ensemble for which the partition function
2:46
is given as, sum over the states exponential of minus βEj or minus Ej by kB T.
3:00
Now I have said in classical mechanics we are pretty much covering the entire phase space so as to speak. So, instead of looking at discrete energy states we can talk in terms of a continuous space because p q phase space is a continuous space so therefore energy may also be assumed to vary in the continuous space and therefore, we started looking at the Hamiltonian that is a continuous function of q and p. And then I can also replace this summation
3:39
by an integration and what we can write is my partition function is something like integration over the q integration over the p exponential of minus Hamiltonian q, p by kB T.
4:08
It is pretty much same as the previous expression except we have replaced the summation by the integration. Since there are 6N variables in the problem grouped into 2 p and q each containing 3 n variables therefore, I am integrating with respect to q with respect to p and this is essentially the Boltzmann factor almost the same as there except that I have replaced the discrete energies Ej with
4:39
a continuous energy function as a function of q and p, so this is what we are trying to evaluate. Now if we are somewhat naive we may imagine that we simply have to integrate this what is a big deal and indeed that can be easily thought of let us first try to simply integrate this and the method that comes to our mind is let us say I want to integrate using some trapezoid rule or Simpson rule or any kind of a quadrature rule let say Gaussian cut quadrature or whatever
5:14
so if I want to do this kind of an integral what essentially I have to do is I will discretize
5:22
each of my variables into some discrete points. Let us say I am looking along some qi variable then of course that qi variable that is 1 of the coordinate variable will be going to be in a range that may be a very large range may be 0 to infinity or whatever but it is going to be in a range. So, let us first take of argument if we want to do a trapezoid rule on that what we will do is we will divide
5:57
that range or domain into discrete points.
6:04
So, let us say we have only 1 variable in the problem then this is what the grid will look like. If on the other hand we have 2 variables let us say we are looking at q1 and q2
6:21
and I am interested in integrating with respect to q1 and q2 then I will draw a 2D grid
6:36
and I will evaluate the function
6:41
at all of these grid locations and the key idea is let us say if I am doing a 1 dimensional integral. Let us say if I have to integrate some function of qi with respect to dqi let us say from
7:00
qa to qb this can be represented as sum of the values at these locations something like
7:15
for trapezoid rule this can be half of f qa f qb plus sum over
7:33
all the k value f qk where k represents the grid locations and these are my grid locations. If I want it to be more accurate, I will just make the grid finer and I will have better accuracy in the integration. I can also use Simpsons rule that is slightly more complicated version of that but nonetheless we are replacing the integration with summation over the grid point that is the key idea that we have, we may not need to have use this particular homogeneous grid we can assume
8:12
a non-uniform grid even then we are replacing integration by a summation that is a key point. Now as soon as we are integrating with respect to two variables then
8:33
this has to be done again within some range let us say q1a, q1b, q2a, q2b. So, these 2 integrations will therefore be replaced by 2 summations and the function values at
9:02
those locations.
9:07
So, long this is sort let us say for example, if you have N grid points
9:17
in 1 dimension for a 1 variable you will have something like N square grid points
9:25
when you are doing 2 variables. Now let us see I want to apply it for this integration so how many variables we have we have 3N coordinates and 3N momentum. So, therefore the number of grid points that we need will be something like let me call M as the number of grid points
9:55
because I am using here N as the number of particles so this will be M to the power
10:03
6N. In other words, we have to imagine a 6N dimensional grid and as you may imagine going from 1 to 2 dimensions, already we went from something like N to N square points. If N was the number of grid points now, we are going from 2 dimension to 6N dimension and keep in mind that this N that is the number of particles is something like 1023.
10:32
So, this is an impossible problem if we are following the standard rules of integration that is one part of the story, the second part of the story is let us say by some imaginary computer that is not there we can do that if we are doing that it is going to be highly inefficient and the reason is because I will be generating points uniformly in this particular phase space as something like this.
11:09
I can also think of some intelligent non uniform grid but nonetheless we are pretty much drawing a grid in the phase space and therefore I am pretty much looking at all
11:26
the locations of grid points in my phase space and my objective should be to cover the phase space as much as I can. But the key point here is this particular integrand that we have this integrand is probably having very small values at many of these phase points because we have an exponential function and we have a negative before Hamiltonian there could be and in most cases, there would be many points actually a majority of points on the
12:00
grid which will have large values of Hamiltonian giving rise to a small value in the integrand and therefore they will contribute very less to the integration. So I will be wasting my effort by computing this quantity at many points in the phase space trying to cover the entire phase space but there are a majority of places in the phase space where the integrand happens to be very small and therefore, we are simply wasting our computational effort.
12:36
So we cannot do this and we do not want to do it even if it is hypothetically possible with current computers we cannot do that, but this is even if it is possible it is not a recommended strategy as soon as the number of particles becomes very large and that is the case for Monte Carlo simulation or any molecular simulation.
13:01
So, with this kind of a glitch now I have to think of how can I integrate apart from the quadrature methods like trapezoidal rule, Simpsons rule and all that which I am more used to doing. So, essentially speaking Monte Carlo simulation is not a method of simulation it is a method of integration because in the end even when I am interested in finding some property it is derived from the partition function and therefore, we will be having some integration in whatever
13:36
property we are interested in. So it is those integrations that Monte Carlo is actually doing at the core of it. Unlike say molecular dynamics where we are solving the Newton laws of motion
13:50
and therefore, that is more sort of a physically driven kind of a scheme Monte Carlo simulation is really a mathematical scheme and it does not have to be a thermodynamic system to use Monte Carlo simulation. In fact in most applications Monte Carlo is applied is beyond thermodynamics in purely mathematical terms. Let me start with a very simple example and then we return back to the problem that we
14:17
are trying to solve. So, one of the mathematical problems that we can imagine is I want to compute
14:31
the value of π. So, we know that π is something like 3.1415 and so on but how do I know that it must has to come from somewhere and I can tell that Monte Carlo is one way we can get a very good estimate of π. So, how will I do that? So first of all I know from where πappears in our understanding. So, π r2 is the area of a circle. So if I draw a unit circle
15:09
that is of radius 1 it is going to have an area of
15:17
π. So, therefore if I can find the area of circle, I can find π.
15:24
So, let us try to use Monte Carlo simulation to find the area of a unit circle. Now how do we define the area so area is area under the any curve is let us say for example, this is my y and this is my x for the circle the equation is x2 + y2 is equal to r2 in this case r is equal to 1. So, x2 + y2 is equal to 1.
15:53
So, essentially my y is therefore plus minus square root of 1 minus x2 and it is that integration of y with respect to x that will give me the area under the graph. Let us say for example, if I look in the top semicircle that is represented by y is equal to under root 1 minus x square and the bottom semicircle is represented by y is equal to under root 1 minus x2 with a minus sign. So, the area under the curve for the top is given by something like
16:44
minus 1 to 1 under root 1 minus x2dx and for the bottom it is given as minus 1 to 1 minus under root 1 minus x2dx and this is the area of the circle. So, if I can evaluate these two integrals then I can get the area of the circle but since for this particular case there is an easier way of doing it and that is, we simply compute the area of a quadrant because ultimately if for example, I compute the area of this quadrant
17:22
the area of the circle is going to be 4 times that. So, instead of doing 2 integrations we can do just 1 integration for this particular quadrant. So we can integrate from 0 to 1
17:40
under root 1 minus x square dx and if I multiply it 4 times I will get the area of the circle. So I get the value of; so A by 4 essentially pi by 4 is 0 to 1 under root 1 minus x square dx.
18:05
So, let us say I want to perform this integration using Monte Carlo simulation so of course we could have used trapezoidal rule in this case and that works fine here because there is only 1 dimension but we want to demonstrate the application of a Monte Carlo scheme. Even for this case we can show that the Monte Carlo is much more efficient for doing the computation of π.
18:31
So, let us get started.
18:34
So, now I am computing
18:39
the area of a quadrant that is from 0 to 1. So, I want to compute the area under this.
18:53
What I can do is, I can imagine a unit square something like this
19:08
that will be a 1 1 point and we know that the area of
19:16
unit square is simply base by height that is equal to 1 this is the construction that I am taking Now instead of trying to find the area under the curve using trapezoidal rule or whatever what I do is I do some kind of a numerical experiment and the experiment goes like following what we say I will generate
19:43
random values of x y in the range 0 to 1 that is within the unit square. So, I randomly generate points between 0 to 1 of x and 0 to 1 of y. So, let us say if I generate 5 points and these are,
20:22
and I am not really when I generate the points, I am not really bothered whether we have a circle passing through that or not. I simply generate points randomly in the domain 0 to 1.
20:34
But now I ask how many of these points will lie within the circle in this particular case in the drawing you have 4 points in the serial regime and 1 point outside the circle. So, if I do this and I ask what is the ratio of the number of points in the circle by total number of points look like?
21:05
So, if you are having difficulty imagining this, think of that there is a dart board and you are simply throwing arrows on the dart board so every time you are hitting the bulls eye you are inside a particular regime every time it is outside that then you are not in the regime. So you simply have to count how many darts you throw go and hit the bull’s eye. And that will tell you essentially the probability of hitting the bull’s eye and it is that probability
21:41
that should depend on the area because if for example, the area of the circle would have been less in comparison to the square if I compare the circle with for example a curve like this
21:57
and let us say I am now interested in area under this curve. So clearly since the area of the circle is higher than the area under this line that is the red region is higher than the green region. It is more likely that I would hit within the red region as opposed to within the green region you can also imagine a smaller area. So, let us say for example we look at area within this that means I am interested in
22:34
the area in the blue regime. So clearly now still if I hit the same number of points it is even lesser likely to be in the within the blue regime. So as the area of the regime of interest or reason of interest is decreasing the likelihood of hitting that will decrease. So, if I really want to have an intelligent game then I should make the bull’s eye on the dart board smaller when it is large it is very likely that you will hit it only when it is smaller the likelihood decreases.
23:14
So therefore what we can say is that the probability of hitting
23:21
a point within the circle should be proportional to the area of that and in fact we can find the probability as the area divided by the total area in which it is thrown because in the limit when I take the square itself of course the probability is equal to 1 for hitting within the square. So, the probability is normalized for the entire square and as we look in smaller regions the probability decreases according to the area the probability is essentially a ratio of the area
24:03
and the area of the unit square. So, in this case this would be the π by 4 that is the area that I am interested in divided by the area of the square that is equal to 1 and therefore my π is going to be 4 multiplied with number of points in circle divided by total number of points.
24:40
And what we have been able to do is in fact integration without even dividing the domain all we did is we simply threw points in that regime. Now clearly if I just throw 5 points you may imagine that the answer is not going to be accurate because it may very well be that if I repeat the experiment with 5, I may have 2 points outside the circle and 3 within circle, it is not that every time you will have 4 within the circle in this example that I am discussing.
25:17
So, therefore we really have to do it for large number of points and I will actually show the code of this in the in the in the coming week where you will see that as I increase the number of points the accuracy with which I can determine π also increases with 10 points I may be around 3. With say 100 points I may get close to 3.1 that is the first decimal, with thousand points I may get
25:47
close to 3.14 with say 10 to the power 6 points I may get close to 3 decimals correct or 4 decimals correct and so on. I simply have to increase the number of points to improve the accuracy with which I can do that so this is the idea on which the Monte Carlo simulation is built on. It is simply how we are integrating the function that we are interested in, it turns out there are details by which we are sampling this. So in this case we are doing what we call a random sampling
26:27
and it turns out not very intelligent for the thermodynamic integration that we are interested in but nonetheless the basic scheme remains the same we want to perform a numerical experiment and do an integration that we are interested in. So, very quickly we can at least develop the ingredients of doing a different kind of sampling that is called the importance sampling.
27:07
So, let us say I want to integrate
27:12
a function with respect to x. So, right now what we have been trying to do is we are simply sampling the regime between ‘a’ to ‘b’. So, I randomly generate points
27:41
in ‘a’ to ‘b’ and using that I am calculating the integral, I mean like I am not really doing the complete math of it but that is the idea that we are following in the in the Monte Carlo simulation. So, now instead of doing that what we could have done is we could have divided this function f by
28:07
some function ρ of x and multiplied with ρ of x dx that is almost the same function. So, in the first instance essentially, we are computing the fi values at for different values of the points that I pick but now what I will do is I will compute
28:38
the value of f by ρ for the points that I am picking.
28:45
And now instead of doing the so called uniform
28:52
probability of picking any point in this regime ‘a’ to ‘b’ so in the previous scheme I can pick any point within ‘a’ to ‘b’ with an equal probability so there is no distinction between different range of different values of the x. But now we can say that we will
29:14
generate or pick points according to the probability density, probability density ρ(x)
29:49
So, now I am essentially doing some kind of a waiting. So, for example,
29:55
let us say if I look in the range of
29:58
‘a’ to ‘b’. Let us say if this is my x axis and let us say this is the probability of picking a point in the first instance it was something like that for random sampling.
30:18
So, I can pick points uniformly in the entire range ‘a’ to ‘b’. Now I can imagine any kind of a waiting function so this becomes my ρ (x) and I will generate points according to that. So, that would mean that I will generate fewer point here fewer points here and more points here
30:46
and why would I do that, I would do that because the value of the integrand is higher in these regimes. So, this is like when the integrand is high or we can say it is a more important regime
31:10
because in this part and this part if the integrand is small then that will not significantly contribute to the integration anyway. So, why I keep sampling points in the regime for which the function value is very small it is not worth doing this was the same motivation I was giving you when I was saying why we should not use trapezoid rule even if it is possible because we will end up sampling many points for which the Hamiltonian value is very,
31:43
very small or very, very large and with so that it does not contribute to the integrand. So, therefore we should focus more in the regime where the integrand is large and therefore, we should generate more points in the reason of interest as opposed to just generating points randomly everywhere and that is the key idea of important sampling that is, I would say what makes Monte Carlo an intelligent method to explore the phase space in thermodynamic systems. So, we will
32:20
never use random sampling for the case of interest in thermodynamics that we are interested in, it can be useful for simple integrations but it is the important sampling that makes the Monte Carlo simulations highly powerful in exploring the phase space in a much better way often than compared to the molecular dynamic simulations. So, with this idea I want to stop here, thank you.
Here to help you learn
Sign in to keep a private conversation with your course tutor.
Ask me anything about the material…
Sign in to chat