0:11

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.

7:02

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.

9:11

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.

10:42

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.

14:25

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