# Finding liquid/vapor phases of LJ system with Gibbs Ensemble

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?

• Your question is a little vague in that you don't provide a lot of specific details about your implementation. What sort of debugging have you done? Have you tried simplifying the simulation as much as possible (e.g. with a trivial test case) and validated it? You mention you don't get the expected densities, is mass conserved during the simulation? Hopefully this can lead you to a working simulation – nluigi Oct 15 '15 at 22:24

I finally discovered the problem, and I'll put it here because it might trip up others in the future. I was iterating by particle: that is, I would choose a particle, decide the kind of move to try, and attempt it. The paper I cited assumes you decide the kind of move to try, then pick a particle to try that move with. For translation moves this doesn't matter. But for exchanging particles, this matters a lot: if you have twice as many particles in one container as the other, then iterating by particles means you will try to move from the more populous container twice as often as you should.

I solved this by changing the entropic part of the move acceptance energy for exchange moves from system 1->2 from $\ln \frac{(N_2+1)V_1}{N_1V_2}$ to $\ln \frac{V_1}{V_2}$ and was able to reproduce the results in the table.