Molecular dynamics simulation of a system of Lennard-Jones particles


This program implements a molecular dynamics simulation of molecules in two dimensions interacting via the Lennard-Jones potential. The output includes the velocity distribution, mean pressure, temperature, and the heat capacity.


The algorithm for simulating the evolution of the particles can be summarized by the following steps:

  1. Compute the forces on all the particles due to the other particles. By Newton's second law the acceleration of each particle is given by the total force on the particle divided by its mass, which we set to unity.
  2. Use a numerical algorithm to update the positions and velocities of all the particles at a time t + Δt later. The program uses the Verlet algorithm.
  3. Accumulate data for various averages.


  1. Choose the initial configuration to be random (anything but crystalline). Run the simulation and observe the velocity distribution. Does it appear to be a Gaussian? Estimate the width of the distribution and record the average temperature. Repeat for several different densities and initial kinetic energies per particle. Plot the temperature versus the width of the velocity distribution. You should obtain a linear dependence.
  2. Determine the ratio PA/NkT for several different densities starting with L = 24. What should this ratio be for an ideal gas? How does the ratio change as the density is increased? For N = 64 and a random initial configuration the linear dimension cannot be chosen to be much less than L = 18 with Δt = 0.01. Why?
  3. Use the default parameters and run the simulation until the system is more or less in equilibrium. Then change the quench rate λ to 0.999. The effect of the quench rate is to multiply each velocity component by λ after each time step. What happens to the temperature of the system when λ ≠ 1? 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 quench rate, you will be able to see a configuration where two phases coexist. Can you tell the difference between a liquid and a gas?
  4. Estimate the magnitude of the temperature fluctuations. Are they approximately proportional to the heat capacity? The latter can be calculated approximately from the numerical derivative of the temperature versus energy curve.

Java Classes

Updated 82 December 2009.