Using a demon to determine the temperature and the chemical potential

Introduction

It is suggested that you first understand the usual demon algorithm before considering its generalization here and why the demon is an ideal thermometer. In this program we generalize the demon so that it can give us information about the temperature and the chemical potential.

    We can picture the usual demon as carrying a sack of energy. The demon makes a trial change in the system and computes the resultant change in energy ΔE of the system. If ΔE ≤ 0, the change is accepted and the demon takes the extra energy. If ΔE > 0, the change is accepted if the demon has sufficient energy to give to the system. The only constraint is that the demon cannot go into debt, that is, E_d ≥ 0. The probability P(E_d) of finding the demon with energy E_d is given by

P(E_d) ∝ exp(-βE_d)

where β = 1/kT. Hence, the inverse slope of P(E_d) is a measure of the temperature T. And because the demon is only one degree of freedom, it does not significantly disturb the system of interest.

    We now generalize the demon algorithm so that it carries two sacks, one for energy and one for particles. The idea is to use this generalized demon algorithm to determine the chemical potential μ. For simplicity, we consider a one-dimensional system of particles, and assume that the positions of the particles are restricted to a grid. We also take the momentum degrees of freedom to be discrete. Thus, the phase space of this system is two-dimensional, with position in one direction and momentum in the other. We can extract the temperature and chemical potential from the probability distribution of the demon energy and particle number, which is given by the Gibbs distribution:

P(E_d, N_d) ∝ exp[-β(E_d - μN_d)].

    The interaction between particles in the system may include a hard core, which implies that no more than one particle can have the same position, and an attractive interaction between two particles that are in nearest neighbor sites (in position).

Algorithm

1. Set up an initial microstate of the system with the desired energy E and number of particles N. The N particles are randomly placed on the phase space lattice with no more than one particle on a lattice site (in phase space). The initial demon energy and particle number is set to zero for convenience.
2. Choose a random site in the phase space lattice. If there is no particle at that site, take a particle from the demon and add it to the site if the demon has sufficient energy to do so and at least one particle.
3. If there is a particle at the random site, remove the particle and give it to the demon if the system energy is lowered or if the demon has sufficient energy to give to the system. Removing a particle can lower the energy of the system because particles have positive kinetic energy or can raise the energy if there is an attractive interaction with a neighboring particle.
4. Repeat steps 2 and 3 many times.
5. Compute the probability of the demon having energy E_d and particle number N_d once the system and the demon have reached equilibrium.

    Note that the total energy of the demon plus the system is fixed, and the total number of particles is fixed. The demon is a facilitator and allows the system's energy and particle number to change by exchanging particles. We could consider a separate Monte Carlo move to move a particle from one phase space lattice site to another, but this move is equivalent to removing and adding a particle in two steps and thus is not necessary.

    We expect the demon to act like a system in equilibrium with a heat bath and a particle reservoir and the demon energy and particle number to be distributed according to the Gibbs distribution.

    After stopping the simulation press Plot Demon Distributions to compute the demon energy and particle number distributions.

Problems

1. Run the simulation with the default parameters, number of particles N = 100, energy E = 200, position space dimension L = 100, and maximum momentum p_max = 10. Begin with no hard core or attractive nearest neighbor potential so that the the system is an ideal gas. We have chosen the mass m of a particle to be m = 1/2 and units such that Boltzmann's constant k = 1, and lattice spacings so that ΔxΔp = 1. In this way the momentum and energy are integers. Discuss the phase space plot. Is p_max sufficiently large so that particles do not reach the edge of the phase space lattice? Are there particles with the same position? Where are most of the particles located? Rerun the simulation for different values of E but the same values for the other parameters. Discuss the nature of the distribution of particles in phase space.
2. Explain why the simulation in Problem 1 is identical to a simulation of an ideal gas in the semiclassical limit. Show that in this limit the chemical potential is given by

   μ = -T ln[(L/N)(πT)^1/2].

   (Note that ΔxΔp = 1 in our units.) Use the default parameters in Problem 1, and run for about 200 Monte Carlo steps per particle (mcs) for equilibration. Then click Zero Averages and average over at least 1000 mcs. (You can speed up the simulation by increasing the steps per display to 10 or 100.) After stopping the simulation press Plot Demon Distributions to compute the demon energy and particle number distributions. The plot of ln P(E_d, N_d = 1) versus E_d should be approximately linear for small E_d. The inverse slope is related to the temperature. The plot of ln P(E_d = 1, N_d) versus N_d should give a slope equal to μ/T. Compare your results with the exact result. What does a negative chemical potential mean? Are there any aspects of the simulation that are not assumed in the textbook derivations?
3. Repeat your simulations in Problem 2 for L = 1000. For this value of L the phase space visualization will not be useful because the particle size is too small. Calculate the distributions. Are the results closer to the exact result? The default setting computes the distributions P(E_d,1) and P(1,N_d). Use other values for the energy for P(E_d, N_d) versus N_d and the particle number for P(E_d, N_d) versus E_d. Do you obtain approximately the same values for temperature and chemical potential? Should you? Vary the values of E and compare with exact results for other sets of parameters.
4. Repeat Problem 1 with a hard core. How can you tell this condition is satisfied by looking at the phase space plot? Repeat the hard core simulation with the other parameters as given in Problem 2. How does the chemical potential change? Vary the total energy until the temperature is about the same as in Problem 2. Is the chemical potential greater or less than that found for no hard core? Explain your results.
5. Repeat Problem 4 but with an attractive well in addition to a hard core. (Enter a positive number for the strength of the attractive well.)

References

• Jan Tobochnik, Harvey Gould, and Jonathan Machta, "Understanding temperature and chemical potential using computer simulations," Am. J. Phys. 73, 708–716 (2005).

Java Classes Used

• DemonLatticeApp
• DemonLattice

Updated 28 December 2009.