Estimating the chemical potential using the Widom insertion method


For a system at a given temperature and volume the chemical potential is given by

μ = -kT ln(ZN+1/ZN),

where ZN is the partition function for N particles. For a classical system the partition function contains an integral over the momentum degrees of freedom and an integral over the position degrees of freedom. The integral over the momenta gives the ideal gas contribution to the chemical potential, and the integral over the coordinates of the particles yields the excess chemical potential μ*. That is, the total chemical potential is given by μ = μideal + μ*. It can be shown that the excess chemical potential can be expressed for sufficiently large N as

μ* = -kT ln<exp(-ΔU/kT)>,

where the average is computed by averaging over random potential additions of a particle. The potential energy ΔU is the change in potential energy of the system that would be found if we were to add a particle to a system of N particles. During the simulation no particle is actually added. The value of ΔU is computed by inserting a ghost particle of the same type at random positions. The average is computed by averaging over many random additions of the ghost particle.

    The program shows the positions of the particles at each Monte Carlo step per particle (mcs) and plots the mean excess chemical potential as a function of mcs. The ghost particle is shown in green. When you Stop the simulation, the program returns the particle density, the mean potential energy, the specific heat, the difference of the pressure <P> - Pideal, the mean excess chemical potential μ*, and the acceptance probability.


  1. Use the Metropolis algorithm to simulate a system of particles at fixed temperature, volume, and number of particles. In our program we simulate a system of Lennard-Jones particles in two dimensions.
  2. Periodically add a "ghost" particle at a random position and compute the change in energy ΔU that would result. In our simulation we do such a trial after each MC step.
  3. Compute the average of exp(-ΔU/kT) and estimate μ* = -kT<exp(-ΔU/kT)>.


  1. Run the simulation using the default parameters. After the average value for the mean excess chemical potential has reached an apparent equilibrium value, press Reset Data, and run until μ* has reached an approximately constant value with small fluctuations. Record your value for μ*. Is it positive or negative? Explain why it has the sign that it does. Remember that μ* is in addition to that of an ideal gas at the same temperature.
  2. Repeat your simulations for different temperatures. Does the excess chemical potential increase or decrease with temperature? Does it change sign? How does the dependence of μ* on temperature compare with the dependence of the potential energy on the temperature?
  3. Repeat your simulations for T = 1 but at different densities. For example try 30 × 30 and 20 × 20 cells. How does μ* depend on the density? Compare its dependence to that of the mean potential energy.


Java Classes

Updated 28 December 2009.