I have the following question regarding a GCMC simulation of SPC/E water using the fix gcmc features. In the ~/examples/gcmc directory, the in.gcmc.h2o script certainly runs perfect, however I do notice that when I set the chemical potential (mu) to very small (negative) value (e.g. -200 kcal/mol), the output seems suggest that there is always 1 water molecule left in the box without being deleted.
My understanding of the thermodynamics is that a negative “infinite” chemical potential should correspond to an empty box (a system having no ability to give away its particle). Can someone please let me know what I miss here?
the issue seems related to if the exchanged particle is an atom or a molecule. Instead of SPC/E water, I tried a fix gcmc run for LJ particles. If the LJ particle is presented to LAMMPS as a 1-atom molecule (i.e. atom_style full) and the fix gcmc invokes the “mol” keyword with a template defined using “molecule” command, 1 molecule (i.e. the 1-atom LJ molecule) always remains without being deleted in the box at negative chemical potentials (I tried -500 kcal/mol).
However, if the LJ particle is presented to LAMMPS as an atom (i.e. atom_style atomic) and the fix gcmc simply choose type “1” (the LJ particle is defined as type “1” in the simulation) as the species to be exchanged with reservoir without using the “mol” keyword, the box can be emptied at low chemical potential, which in my own view is the correct behavior.
The LJ particle being an atom or a 1-atom molecule should not matter, and LAMMPS (at least based on my test) seems to have different behaviors. I don’t want to believe this is a bug in the code, but could someone please comment on this?
At first I thought this was not a very well-defined question. But your test of LJ being an atom or a molecule is a good one and indicates that there really is an issue, probably a bug/feature in the implementation of molecule support (which involves a lot of extra code). I confirmed that this is happening. It is an undocumented feature of using the mol keyword. The molecular count is never allowed to drop to zero, as this creates problems for fix rigid/small. This restriction should be limited to rigidflag==TRUE and should be documented.
Thanks so much for your clarification and confirmation of the problem!
Based on your comment, I get the impression that fix gcmc is not suitable (or at least should be handled with care) for simulations of vapor phase and calculation of vapor-liquid coexistence under the grand canonical ensemble, as the vapor is frequently depleted unless the box is made so large or the vapor pressure is so high that there are always molecules in box. For water at not so high temperature, this seems not very practical. In other words, there should be a lower limit of the chemical potential one can use in the input file, and that limit is dependent on the box volume.
I do not view this as a practical limitation on the implementation of GCMC in LAMMPS. Any molecular simulation that is attempting to resolve the difference between <N> ~= 1 and <N> << 1 is of questionable value. Moreover, any system for which <N> << 1 can be better analyzed using non-interacting gas molecules e.g. Widom insertion (see fix widom).
I agree with you that for many simulations << 1 belongs to part of the phase space that is of limited interest, however, the suboptimal behavior of LAMMPS to correctly sample small chemical potential values that resulting << 1 limits its ability to perform histogram reweighting to construct phase coexistence curve.