[MUSIC] Hi, and welcome again to our class, Simulation and Modeling of Natural Processes. In the present module, we go back to one of the potentials which we have introduced, the Lennard-Jones potential. And we'll talk about about how to solve interaction with this potential in an efficient way by introducing a cut off distance. You remember from the previous module that I showed you a way to integrate Newton's equations of motion, in the case of a system of many bodies or many particles. We introduced a Verlet scheme in which we took the system from non-discrete time to the next discrete time. Now we saw that this discrete scheme was very efficient computationally, but there was one expensive operation hidden in this scheme. This expensive operation is the calculation of the forces acting on the bodies around the particles which we use to evaluate the acceleration of the particles then integrated. Let's go back to the case of N particles interacting with each other through a potential, in this case the Lennard-Jones potential. As a matter of fact, it could be any other potential. In principle if you want to calculate the force acting on a given body, in this case on body one, you need to take into account the force which comes from all the other bodies which will perpetually act on the given body. I show you now an example with a total of ten bodies, which means if you want to know the force acting on body one you need to calculate nine interaction terms. Of course, in a bigger system if you have a galaxy with 400 billion bodies, you will have to calculate 400 billion forces acting on the given body if you do it in a brute force approach. We say that, in this case, the algorithm has a complexity ends clear which means the algorithm are just computing all the forces on all bodies acting on the system at a given discrete time, and takes a number of operations which is proportionate to n squared. Because, first of all, you have to calculate the force for every one of these n bodies. And then for every one of these n bodies, you have to calculate the interaction with all the n minus one other bodies. Of course in the case of gravity, and as a matter of fact in the case of all the other forces which we introduced here, because of Newton's Third Law, the forces are symmetric. The force of body A acting on Body B Is equal in norm to the force of body B acting on body A. So we really only have to evaluate half of these forces. So, the total of amount of operations is n(n-1)/2 but still the leading term of this expression is n square so this is a n square algorithm which can be extremely expensive, except if you apply an intelligence scheme to shorten the computation. Let's see what a complexity of n square means. On this graph, I have plotted on the x-axis a value of n which ranges from zero to ten to the power 12. Along the y-axis I have simply plotted the square of this value. This graph is plotted on a log-log scale, which means that both the x and the y axis are logarithmic. On the log-log scale, the n square curve is a linear curve. It shows us how many calculations we need to do to perform to compute all forces at one iteration of our time stepping scheme. And I have showed you to give you an order of magnitude which type of computers we have available nowadays can compute which amount of interactions. Of course, the number of interactions a computer can evaluate depends on many factors. It depends of the type of potential which is being used. It also depends on the kind of algorithm you use to implement interaction between the bodies on the data structure which is being implemented. And you will see shortly different type of data structures which we will use in different types of contexts. So the computational power, which I am estimating here, is really just a rough estimation, and you should take it as an order of magnitude and not an exact value. You can see here that if you are calculating these interactions with a desktop computer, you can just compute the interaction between a few thousand bodies or a few thousand particles, and that's all. The worlds fastest super computer, if it uses a brute force n square algorithm to evaluate all interaction, can just handle a system of around 1million particles and not more. This is not allowed. The number of stars you need to simulate fi you want to simulate our galaxy is approximately 400 billion. In this case, if you take the square of 400 billion, you will see that you have so many interactions to take into account that you will need 100 billion of the world's fastest supercomputer just to evaluate these interactions at one time step in a reasonable time. This is completely unrealistic. It's unrealistic nowadays, and also in the near future, you will not be able to reach this computational power. This shows you that it is not possible for interesting problems like the stars in the galaxy, or even where the molecules which you have in a gas, to use a computer to evaluate all the n square interactions. We need to find shortcuts, ways to simplify the problem to make the computations faster. Again, to summarize the problem we are treating, a brute force approach for treating a system of n objects is ten square. We cannot solve that nowadays with the available power. Strategies are required to reduce the valid calculations. But let's go back to the Lennard-Jones potential. Because this is the first potential which you are going to treat you remember that I told you that, most of the time, molecules which share a common space with other molecules don't see each other. As opposed to stars, which constantly are influenced by all other stars in the system. Molecules, most of the time, just travel along straight lines. They are not affected by the other molecules. That's because the interaction potential of the Lennard Jones model is a short range interaction model. It decreases to the power six of the distance between two molecules. As you can see on this graph, if you go just at a distance of two to five times the radius of the molecule, the potential is almost zero. Actually, it amounts to less than 1% of the potential of epsilon. It becomes negligible. Therefore, it is common practice when you simulate such a system to introduce a cutoff length. A cutoff length after which you will neglect all the other molecules, because you consider their influence on the present molecule to be zero. Common value of this cutoff length is a length of d equals 2.5 times sigma. At which, as I said, the value of the interaction is only 1% of the depth of the potential well and becomes negligible. It can, of course, be another value, if you want to improve the numerical accuracy of this approximation. But the cut off value, does not give you an an answer for the most fundamental problem. How do you actually find nearby particles? How do you find particles which fall inside the radius of the cutoff value of dc? If you were to check all the particles to select the ones which are within a reasonable range, you will need to do n operations, loop all the other particles to do so. You would not reduce the overall algorithmic complexity of n square because for every particle, which means an operation of n, you would need to check all the other particles. Again an operation of order n. You will again have an n square algorithm. You will gain nothing. You need a smart data structure, which makes it possible to quickly find particles which are in the neighborhood of a given particle. The solution here, is to use a grid based method, which means we take the space occupied by the particles, and overlap it with a grid, which has equal-sized cells. In a computer simulation, we represent this grid by a matrix. Every element of this matrix represents a cell of the grid which is overlapped with the physical space. Let's take for example the grid which has the coordinates one, two in this system, and let's assume that there are five particles currently present in this cell. These particles have a given ID, let's say for example ID [2, 3, 7, 11 and 19]. Which means that when you are doing a computer simulation here, every grid cell will contain a list of IDs, of particles, which are present on this grid. Grid cell one, two will have a list with the values [2, 3, 7, 11, and 19]. We're not going to store the position and the velocity of the particles inside this matrix. This matrix is only here to identify particles at a given position in space. We use these IDs to look up the positions and velocities in some other kind of data structure Now, as before, when we talk about data structures, we need to create this data structure. At some point, these particles move, and whenever they move, we need to reassign them to the right cell inside this grid. Unfortunately, this is quite a fast preparation to perform because given that a particle which has changed its position, you can find out to which grid cell it is being assigned in constant-time, because the space occupied by grid cells is a constant value which we know in advance. We can take the particle position, and in constant time, assign it to the range of a given grid cell, and then add it to the list contained in the matrix an element of this grid cell. Assigning all particles of a given time stat to their corresponding grid position inside the matrix is a linear operation, an order n operation, which is much better than an order n squared operation. Once they have been assigned to the grid, we can exploit this matrix to find nearby particles and evaluate the forces acting on a given particle. Let's take, for example, the particle which in this image is red. It's the particle with the ID 11. The green circle shows the distance of the cutoff length dc. We will only consider the interaction of the red particle with particles which are contained inside the circle. To find these particles, instead of looping over all other possible particles, we will use the grid to find possible candidates. In this case, we just need to access the list of particles of grid cell 1,1, 1,2, 2,1, and 2,2. We access four grid cells, access their lists, and take all the particles contained in these lists as possible candidates for particles contained within the circle of cutoff length dc.. Once we have selected them, we need to do an additional operation to check if they actually are contained within the circle, and then select them to calculate the forces acting on the red particle, particle number 11. What's the complexity of doing so? Let's assume that the single cell of our grid is larger in width than the cutoff distance. This means that, in the worst case, when we need to check for possible dates of neighboring particles, we need to check in a 2D simulation at most, nine cells. The cell in which the particle is located on to which we are calculating the forces, plus, it's eight neighboring cells. In 3D, this would be 27 cells. If we have Nc cells in total, the average number of check is n for the amount of particles times nine for the neighborhood, 9 N divided by Nc where the total amount of particles is divided by the amount of cells, to get the average number of particles per cell. If we have many cells in the system, well if the number of cells in the system is of the order of the number of total particles, and this is the case if the cutoff lengths is much smaller than the size of the full domain, then this algorithm has a complexity of order N, which means it is a linear algorithm. It is an algorithm which is much more efficient than the root force and square algorithm This ends the module in which we talked about the Lennard-Jones potential, and the technique of the rate based method and cutoff distance to reduce the calculations from order N square to order N. Stay tuned for next module.