Using pressure in GCMC

Hi all,

I am checking the fix gcmc feature on a very simple problem, before I try it in a more complicated setup.

My simulation is rather simple: a cubic box of lj atoms, with pair coefficients eps=1.0 and sigma=1.0 (all in lj units), and cut-off 2.5. If I run an NPT simulation of that system at pressure 1.0 and temperature 1.0 (lj units), the resulting equilibrium density is 0.8.

Now, if I run an NVT simulation of the same system at T=1.0, starting at a lower density (let's say 0.5) and coupling my system with a reservoir at T=1.0 and P=1.0 through the command:

fix 3 all gcmc 100 100 500 1 4839843 1.0 1.0 0.01 pressure 1.0

The density of my system can go up to 0.84 (and even higher). From the top of my head, I can't remember what is the jamming density for a LJ fluid, but I must be getting closer. How is that possible? Am I missing something? Is there a bug in the insertion/deletion algorithm? Have I used the wrong P in the GCMC reservoir?

Any help will be highly appreciated.

Kind Regards,


Hello Jorge. Isn’t clear to me that there’s any problem with your input or with the fix GCMC code. I’d recommend checking your simulation result against an LJ equation of state. If there’s a clear disagreement, then I’d recommend double checking your input. If that still doesn’t solve it, there may be an issue with the code that we can look into. If you think this is the case, please send me your input and what you think the right answer should be, and we’ll go from there.

The molecule and pressure features of fix GCMC are fairly new, but the atomic and chemical potential specification features are not and have been pretty well tested. So I’d recommend trying without specifying the pressure on your atomic system to start out. It would also be simpler to start with pure GCMC and no MD steps by not specifying an MD integration fix.