Thursday, May 3, 2012

Programming the Ideal Gas

You've heard about it before. Many times. It's a model system used a lot in chemistry. The Ideal Gas law: A system where the particles can be treated as non-interacting. How do you program such a system? What do you need to describe it? Given 1) the basic introduction covered by Prof. Jan Jensen in another lecture, and 2) the following assignment (notice it is part 2), here is what the students are supposed to do once this exercise is complete:

  1. Take a box with dimensions $(-1<x<1$ and $-1<y<1)$.
  2. Make npart particles placed randomly in the right side of a box ($0<x<1$) but evenly from the top through the bottom
  3. Give the npart particles random velocities in both directions and make it so they can move in both the positive and negative direction.
  4. Move the particles around via the equation: $X_i^{n+1} = X_i^{n}+dt*dX_i^{n}$
  5. Keep them inside the box at all times.
  6. Plot the result.
The students start out with 1, 2 and 6 presented to them, both for inspiration and because they contain some subtle hints. The code is

One of the tricky parts is knowing "what is the i'th particle" because we need this information to update the system according to the equation above in point 4. Once you realize that you can access this information in the X, Y, dX and dY arrays as X[i] under a loop, then there is really no problem with that point either, for instance, to update all particles according to the equation in 4 in one go, we could write

Assigning a loop to iterate over the particles is not really troublesome, but understanding why is hard - even to explain it to yourself. A simple for-loop might not be the smartest thing to do, but I think it is the easiest one to go about and more importantly I often feel that many things can get lost in the whole numpy magic, especially for newcomers.

The next part is, that we should keep the particles inside the box at all times. The exercise has a little discussion about which version is the better, but in my opinion the particles should always be inside at all times. To keep them in, we ask whether a step (equation in point 4) in X and Y will result in the particle being outside the box and if that is the case, we reverse velocity that potentially brought the particle outside. We could use the following code to do that

Inserting the above two snippets of code, under a loop over the number of steps we want to take, we finally end up with the code that simulates non-interacting particles in a container.

Try and code it for yourself or use it in a course you teach. The students do find it interesting, both in the form of "what do I need to tell the computer to make it do stuff?" and "oh - the particles do move!".

Disclaimer: This code is part of a course I help teach at the University of Copenhagen - just imagine you never saw python before.


  1. Really funny :)

    In your final script you have renamed stepsize to dt, but not on the lines where you check if the atoms have moved outside the box. :)

  2. Søren,

    thanks for spotting that! :) it has been corrected now.