[lammps-users] How to set mole fraction of each component when using multiple gcmc commands?

Dear lammps-users:

I have a question about using multiple fix gcmc commands to simulate binary gas adsorption:
How to set specific mole fraction of each component when we use multiple gcmc commands?
For example,  I want to set the mole fraction of methane and butane as 0.9/0.1 in gas reservoir with fixed total pressure. 
May I implement this by setting pressure value of each component according to 'the molar fractions  of the components are equal to the fugacity ratios' ? Do I have to calculate chemical potential of each component from 'fix widom' to implement this?
Any suggestion could help a lot. Thanks.


I am not much of an expert in this kind of calculations, so please take the following advice with a sufficiently large grain of salt. This is mostly based on common sense and a little bit of technical knowledge of the code itself.

I am not aware of a ready made solution to your problem. Please keep in mind that when using multiple instances of fix gcmc that they don’t “know” about each other, thus I would try to set things up in a way that they are not called on the same timestep (whatever is the second fix would have systematically different conditions than the first) but on equal distances. That is, if both fixes would be called every 100 MD steps, then I would try to have the second one added/started 50 steps later, so that the first would execute on steps 0, 100, 200, 300, etc. and the second on steps 50, 150, 250, 350, etc.
however, I am not certain that using different settings to achieve different absorption rates is straightforward and would not know how to derive or justify the choices. One thought could be to use the same settings, but make it so that one fix executes more frequently than the other in order to represent that ratio.

Another thought would be to modify the source code so that there are two molecule IDs and two molecule specific groups to which atoms get added, so that you can specify the ratio in the reservoir and then first pick a random number to choose which molecule to pick based on your intended ratio.

My biggest concern would be how to correctly represent the ratio of the two molecules in your reservoir. If there is a significant difference in how easily one is inserted into the system compared to the other, then the ratio or 9:1 should only be reflected in how often the insertion of one molecule is tried over the other, but not by how many of each kind are finally inserted. All of that is assuming the reservoir is infinite, if you have a setup that would be more appropriately described by a finite reservoir, then also the change of the ratio in the reservoir would have to be considered somehow.

hope this helps,


Instead of starting with what you want, try starting with what the code actually does, as described in the documentation, or in any textbook description of the GCMC method. That way, you can avoid wasting time imagining calculations that can not be performed, and instead focus on which actual calculations might give you what you need. Here is a nice summary phrase from the first sentence of the doc page:

“exchanges atoms or molecules with an imaginary ideal gas reservoir at the specified T and chemical potential (mu) as discussed in (Frenkel).”

So, the only knob you really have is the chemical potential “mu”. If you invoke fix gcmc twice, then you have two knobs, one for each species. What those knobs actually do is highly dependent on things like temperature, interatomic potential, as well as more subtle things like system size and initial composition. However, in the low density limit, the effect of the knobs should be very predictiable because in that regime:

mu = kT log(N/V) + constant

This is true independently for each species, but the constant will be different for each. Knowing all of this, you are now in a good position to do you own computer experiments. Starting with a box of the right size to hold a few dozen molecules at some moderate density, not too big or small, you should quickly be able to figure out the mu-knobs do.


p.s. fix gcmc provides some modified knobs like “pressure” and “fugacity_coeff” that are really just alternative ways of defining mu.