Monte Carlo simulation of a system of Lennard-Jones particles
Why use a Monte Carlo method to simulate a system of particles interacting by the Lennard-Jones potential when it is possible to use molecular dynamics methods? After all, molecular dynamics gives the same information about the static properties of the system such as the equation of state and gives dynamical properties in addition. The simple answer is that Monte Carlo methods are more flexible and sometimes more efficient. In particular, the initial configuration can be chosen at random if Monte Carlo methods are used, but we have to be very careful that no particles are too close to each other if we are going to solve Newton's equations of motion. Another reason is that Monte Carlo methods can be used to simulate any desired ensemble.
The program uses the Metropolis algorithm which simulates a system of particles in equilibrium with a heat bath at temperature T. The algorithm for simulating the configurations of the particles can
be summarized by the following steps:
- Establish an initial configuration. Because the velocities of the particles are irrelevant, we can ignore them and consider only the positions of the particles. The energy of the initial
configuration is not important.
- Choose a particle at random and make a trial change in its position.
- Compute ΔE, the
change in the potential energy of the system due to the trial move.
- If ΔE is less than or equal to zero, accept the new
configuration and go to step 8.
- If ΔE is positive, compute the quantity w = e-βΔE.
- Generate a uniform random number r in the unit interval [0,1].
- If r ≤ w, accept the trial move; otherwise retain
the previous microstate.
- Determine the value of the desired physical quantities.
- Repeat steps (2) through (8) to obtain a sufficient number of
- Accumulate data for various averages.
- Determine the ratio PA/NkT for several different densities starting with L = 24 and N = 64. What should this ratio be for an ideal gas? How does the ratio change as the density is increased?
- Determine the temperature dependence of the mean potential energy and the heat capacity.
- Describe the qualitative behavior of the radial distribution function g(r).
- Use the default parameters and run the simulation until the system is more or less in equilibrium. Then decrease the temperature gradually. Is there evidence of a phase transition between
a solid and a liquid? If you increase the number of particles and increase the linear dimension of the system and choose different values of the temperature , you will be able to see a configuration where
two phases coexist. Can you tell the difference between a liquid and a gas?
Updated 28 December 2009.