I've been trying to use the Gibbs Ensemble Monte Carlo method for finding the liquid and vapor phases for a given particle interaction. This method sets up two isothermal subsystems that exchange particles and volume, but otherwise don't interact. In principle, these two systems should equilibrate so that one contains the liquid phase, and the other the vapor phase. This gives you the densities of each phase in a single simulation.
The problem is, when I test it with a Lennard-Jones potential (diameter $\sigma$), I can't reproduce the results given here (DOI 10.1007/BF01458815):
($L^*$ is the size and $\rho_g^*$ and $\rho_l^*$ are the densities $N\sigma^3/V$ for the two phases).
I find I get results pretty close to this for the liquid phase at $T^*=1$, but the gas phase is much denser than it should be -- $\rho^*_g=$.1 to .2. Worse, the density of the gas depends on the initial density of the whole system, which should not be the case.
I've built the simulation in two completely different languages and gotten the same results, so it isn't a simple bug. The systems correctly exchange volume and particles, and particles inside each arrange themselves as I would expect. The only thing I do differently than the algorithm is to rescale $\sigma$ to change volume (instead of changing the box size and rescaling the particle positions).
What might be going wrong here?