fix GCMC for water in slit shaped pores

LAMMPS LAMMPS Mailing List Mirror
1

Greetings LAMMPS users,

I have a few questions about the MC package “fix gcmc” command. I have both conceptual and implementation related questions. My geometry consists of a graphene slit pore (periodic in x and y direction, confined in z-direction) inside which I want water molecules to be filled for a given thermodynamic state (ideal variables that I would like to specify are Temperature and Pressure). I would like to compute the adsorption isotherm. After reading the description of the “fix gcmc” in the user manual, here is a list of questions that I have come up with:

  1. Most GCMC bare-bones algorithms that I have read either use chemical potential “mu” or fugacity “f”. Fugacity can be obtained if fugacity coefficient and Pressure (of the bath or reservoir) are known. Fix gcmc uses both “mu” (as a mandatory input) and pressure and fugacity coefficient as keyword/value (optional) pair. According to the manual, the usage of one parameter suffices for the simulation. From the isotherm data, I wish to compute “partial (ln f) / partial (ln pore_density)”. Should I just vary the fugacity coefficient setting keyword pressure in “fix gcmc” to 1 bar to obtain the “fugacity vs pore density” data? Would this be physically accurate?
  2. How do I write 1 molecule topology (model molecule in data file)? This is required since I would be simulating SPC/E water. I have performed MD simulations with LAMMPS before, and I always specified the coordinate file with “N_frozen_solid” carbon atoms and “3*N_water_molecules” water constituents (oxygen and 2 hydrogens), with angles and bonds specified. Then in the input file, I group together all carbons and all water constituents as separate groups. Should I keep the solid atoms as is, and just specify 1 water molecule, with 2 bonds and 1 angle. Would it be all right? Is there any restriction on where it should be placed relative to the carbon atoms?
  3. In order to avoid the interaction of periodic images, boundary is set to be “p p p” (where z-direction is the confining direction). One cannot set it as “p p f” since there are PPPM calculations involved, hence it would return an error asking for periodic boundaries. As a workaround, one can set some distance between the periodic images in z-direction (by increasing box length beyond wall), so that there is no interaction between the periodic images in z-direction. My question is, if a buffer region exists between two periodic images, how would one ensure that the water molecules are only manipulated inside the physical domain in z-direction, and not in the buffer region. Would the keyword “region” in fix gcmc provide any resolve? If so, how exactly? Or am I completely missing something here?

I understand that some of these questions could be very trivial, if that is indeed the case I would be really grateful if you would point me towards proper references to look into. Your help/comments are highly appreciated. Thanks.

Ravi Bhadauria
University of Illinois
Urbana

2

Paul can answer these.

Steve

3

  1. The approach you suggest sounds fine to me. Note that the fix GCMC code uses the fugacity_coefficient and pressure only to compute the fugacity as the product of the pressure and fugacity_coefficient. The default fugacity_coefficient value is unity. So, if you wish to set the fugacity, just make sure that fugacity = pressure * fugacity coefficient. Here’s the line from the code where pressure and fugacity_coeff are used: “if (pressure_flag) zz = pressurefugacity_coeffbeta/force->nktv2p;”

  2. Again your approach sounds OK. It doesn’t matter where you place that first water molecule (any sensible location), and you can even delete it after initializing fix GCMC.

  3. Note that fix GCMC does not include long-range interactions, so I wouldn’t recommend running fix GCMC with a kspace style on. Also note that your geometry is what I’d call a slab geometry and there’s a special flag for doing slab geometries with kspace on. See the slab option under the kspace_modify command: http://lammps.sandia.gov/doc/kspace_modify.html So you can use “boundary p p f” with PPPM. To insert/move/rotate/delete water molecules only inside your slit pore, you can either use the slab command, or the fix GCMC region command, or both. The slab command should be sufficient as long as your slit pore encompasses the entire physical domain.

  4. One final note: the fix GCMC molecule option is still fairly new, and has not been thoroughly verified, so users should check their results for correctness and report back any issues they discover.

Paul

4

Thank you Paul and Steve. I would run some benchmarks and post again if I run into any trouble.

5

Paul and others,

I have few issues with the GCMC routine again. Upon inspection of snapshots after some intervals, I see that a couple of water molecules have moved out of physical domain, beyond the graphene wall. I have kept the boundary “p p p”, with some buffer region in z-direction so that periodic images in z do not interact. I believe that this is happening in the translation of gcmc group as written in the manual:

“This command may optionally use the region keyword to define an exchange and move volume. The specified region must have been previously defined with a region command. It must be defined with side = in. Insertion attempts occur only within the specified region. Move and deletion attempt candidates are selected from gas atoms or molecules within the region. If no candidate can be found within the specified region after randomly selecting candidates 1000 times, the move or deletion attempt is considered a failure. Moves must start within the specified region, but may move the atom or molecule slightly outside of the region.”

I cannot avoid this, usage of a “p p f” boundary style, with no buffer region results in “ERROR: Lost atoms”, possibly happening because the translation moves them out of simulation box. How can this be circumvented?

Other issue is with the usage of fix/wall command. Instead of having an atomistic description of graphene, I would also like it to compare the results with 10-4 type potential (taken only first two terms from this expression: http://www.sklogwiki.org/SklogWiki/index.php/10-4-3_Lennard-Jones_potential), which I have already hard coded as a fix (only with real units). Again, “p p f” has to be used (mandated by fix/wall), but I can choose whether I want a buffer region or not. In either case, the simulation is susceptible to a crash. When no buffer is specified, it crashes due to lost atoms. If a buffer is specified, as soon as the molecule translates to buffer region, the simulation crashes with error “Particle on or inside fix wall surface.”

In my opinion, the region handling in GCMC should be robust in not letting a particle translate outside the bounds. If there is a workaround, I would be more than happy to know about it. Please let me know, if you want to test the repeatability of the errors so I can provide you with appropriate setup files.

Sincerely,
Ravi Bhadauria
University of Illinois

6

Ravi,

The issues you describe are not altogether unexpected, but I can see how the behavior you’re observing is not what you want. It may boil down to the fact that the fix GCMC moves/insertions/deletions/rotations only consider intermolecular interactions and ignore all other imposed energies like the wall you mention. In determining whether or not to accept a move, those fix wall interactions are not considered.

You’re right that use of the region option in fix GCMC does allow moves of molecules or atoms that start inside the region to move out of the region. This is not an issue with robustness, but rather a choice in the design of how the region option works. The code could be hacked to do otherwise.

So in summary, I don’t see these issues you describe as errors, but rather that the code isn’t doing what you wanted. Perhaps the best path forward would be for you to make code modifications as you deem appropriate, and send it back to us if you’d like. If we agree that they are improvements or useful options, we can add it to the distribution.

Thanks,

Paul

7

Thanks Paul. I would try to read the fix_gcmc source code and would ask you if I have any questions.

Ravi