[lammps-users] fix gcmc problem - water filling in slit pore

Dear LAMMPS community,

I’m using fix gcmc to insert flexible SPC water molecules (constrained by fix shake) into 3d periodic boxes at 300K to simulate sorption in a 4 nm-width silicate slit pores with 3 nm walls at saturation, with the 29 Sep 2021 version of LAMMPS. The “full_energy” option is enabled since the pair style is lj/cut/coul/long and PPPM solver used. The chemical potential I used was acquired by adjusting the molecule file and the value of mu in the example “in.gcmc.h2o” script and see when the liquid phase formed (about -10.5 kcal/mol, which is near to literature value). The insertion of water molecules in pore wall structures and on walls surfaces are working, but I expect full filling of the pore space, rather than a few layers adsorbed to the surface. I tried adopting -8 kcal/mol for chemical potential, and the full filling occurred, but this value seemed to be way too high. I have also tried original rigid SPC model, but it doesn’t work better. A large number of MC steps adopted because the initial configuration should be filled with water.

My questions are:

(1) Whether I had got the value of chemical potential wrong? Should I use a different chemical potential or adopt another way for deriving chemical potential? The exact saturation vapor state is important to me because I want to simulate sorption.

(2) I found that a previous thread in the mailing list stated that with “full_energy” option on, “the large and positive intramolecular energy which causes insertion events to have a very low probability of acceptance”, and the “intra_energy” option is thus provided. I am not sure that whether this case is suitable for SPC water models (rigid/flexible with fix shake). If so, how should I calculate and specify this “intra_energy” value?

(3) Another confusing issue is that some slightly different silicate pore wall (3 nm, slightly different ion distribution, file A and B) with unchanged geometry (4 nm void) may induce pore filling. In most cases the same scenario occurs (like in file A), but in a few cases (file B, rare, no regularity found), a lot of water molecules will quickly appears in the pore void (even with mu of -10.9 kcal/mol). Why is this happening? The further nvt relaxation or minimization of pore wall do not relieve the problem.

(4) Is there anything else I can do to solve this problem?

I attached two .lammpstrj and .log file to this e-mail. The file A and file B are same in wall composition, but little different in distribution of ions due to random selection. File A has a mu of -10.5 kcal/mol, but no full filling observed; while file B has a mu of -10.8891 kcal/mol, but full filling can be observed on the very first step.

Much thanks in advance!


fileA.lammpstrj (782 KB)

fileB.lammpstrj (950 KB)

fileA.log (7.07 KB)

fileB.log (11.2 KB)

Rather than solving your issue, I’ll answer the question no. 1.
I presume, you have a really dense system at some point of the insertion of molecules.
It has been proven that with standard Monte Carlo insertion/deletion moves it is hardly possible

to obtain the correct chemical potential since as the system gets denser, it is very unlikely for the insertion to occur.
This has been solved by the use of different biases, e.g.
Siepmann, Frenkel, MOLECULAR PHYSICS, 1992, VOL. 75, NO. 1, 59-70

or “cavity-biased” methods

Those are not implemented in the LAMMPS and therefore the use of really dense systems may lead to unphysical results.

Why don’t you somehow prepare your initial configuration and then simulate the sorption, instead of overcomplicating things?

Best wishes,

Lukasz Baran
Ph.D. Candidate
Department of Theoretical Chemistry
Institute of Chemical Sciences
Faculty of Chemistry
Maria Curie Skłodowska University in Lublin

ORCID: 0000-0003-1777-1998
ResearcherID: AAP-1152-2021


insert flexible SPC water molecules (constrained by fix shake)

This is very bad. The parameters for a flexible SPC water are chosen so they represent relevant properties well when the molecule is allowed to be flexible.
Its parameters are different from rigid SPC (or SPC/E or TIP3P etc.) for a reason. So you are effectively making up a new water model without any justification and one that is very likely to represent water worse than either the (published and validated) flexible or rigid SPC model.


I have also tried original rigid SPC model, but it doesn’t work better.

This is a broken argument: you “improve” just one specific property but at the same time, everything else become undefined and likely (much?) worse. Not good.


(4) Is there anything else I can do to solve this problem?

You are trying to use LAMMPS for something it was not designed to do well. You’ll be much better off using a real monte carlo simulation software, but also need to adopt a more rigorous approach to the physics involved.
Just picking different settings because they appear to give “better” results is not a good approach. You will always need to be able to provide a physical justification for that, and - like in the case of changing potential parameters/settings - will have to re-validate them to be suitable for the thermodynamic conditions that you study and to confirm that by improving the behavior in one aspect you are not making matters worse in (multiple?) others.

Thank you Dr. Kohlmeyer and Łukasz! I will try to use MC softwares and adopt a more physical approach.

Best wishes,