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, Ed ≥ 0. The probability P(Ed) of finding the demon with energy Ed is given by
P(Ed) ∝ exp(-βEd)
where β = 1/kT. Hence, the inverse slope of P(Ed) 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(Ed, Nd) ∝ exp[-β(Ed - μNd)].
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).
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.
μ = -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(Ed, Nd = 1) versus Ed should be approximately linear for small Ed. The inverse slope is related to the temperature. The plot of ln P(Ed = 1, Nd) versus Nd 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?
Java Classes Used
Updated 28 December 2009.