GCMD Simulation Dilemma at High Pressure (10-30 atm): MAXENERGYTEST Stalls (PPPM/DSF) vs. Severe Jamming (Cutoff)

LAMMPS LAMMPS Beginners
1

Hi, all.

I am performing hybrid GCMD (Grand Canonical Molecular Dynamics) simulations to study the adsorption and permeation of CO2 (TraPPE model) through a MOF (with open metal sites) along the Z-axis at a high pressure of 10-30 atm.

My LAMMPS version is 22 Jul 2025 - Update 3.

I am facing a fundamental dilemma regarding the handling of electrostatics and the fix gcmc command. My system consists of the MOF structure and an empty Feed reservoir (buffer space) outside the MOF along the Z-axis.

Here is what happens depending on the electrostatic settings:

1. Approach A: Using kspace pppm (with full_energy) The simulation runs normally at first. However, when CO2 molecules adsorb onto specific sites within the MOF (likely the open metal sites), it triggers a massive energy spike. This results in the following warning: WARNING: Energy of old configuration in fix gcmc is > MAXENERGYTEST. (src/MC/fix_gcmc.cpp:798) Once this happens, the full_energy evaluation causes the GCMC engine to completely stall (rejecting all insertions in the feed region).

2. Approach B: Using simple coul/cut 12.5 To avoid the full_energy global rejection, I switched to a simple cutoff model. Interestingly, this approach does not throw the MAXENERGYTEST warning. Gas insertion works continuously, but due to the lack of long-range repulsion, GCMC severely over-inserts CO2. The molecules jam completely at the feed/surface interface and cannot permeate into the pore channel at all. Even when I dropped the pressure to 10 atm, this severe physical jamming persisted.

3. Approach C: Using coul/dsf 0.05 12.5 I then attempted to use Damped Shifted Force, hoping it would prevent the over-packing of coul/cut without the global energy issues of PPPM. However, the result is exactly the same as Approach A. No physical jamming is observed, but when CO2 attaches to specific MOF sites, it throws the exact same MAXENERGYTEST warning, and the GCMC engine completely stalls.

Current Status & Workaround: While researching, I found a related discussion (GitHub issue #409) suggesting that normal MD fluctuations might trigger this warning if overlap_cutoff is too large. I have just reduced the overlap_cutoff in fix gcmc from 2.5 Ă… to 1.5 Ă… and am currently running a test to see if this prevents the false MAXENERGYTEST stalls for Approaches 1 and 3.

My Questions:

  1. Have others experienced this sudden MAXENERGYTEST stall when gas molecules attach to specific highly-charged sites in MOFs (Approaches 1 & 3)? Does tuning overlap_cutoff down to 1.5 usually resolve this, or is there another fundamental algorithmic fix I am missing?

  2. To stabilize these localized energy spikes, is “Site Passivation” (artificially setting the charges of the specific outermost metal sites to 0.0) an acceptable and valid practice in the MD community for such GCMD setups?

  3. If I am forced to use coul/cut (Approach 2) to prevent stalling, how do researchers typically calibrate the GCMC pressure or fugacity_coeff to avoid the severe over-insertion and jamming caused by the missing long-range repulsion?

Any insights, feedback on my approaches, or references would be greatly appreciated. Thank you!

2

Since nobody else seems to be willing to give some advice, here are some general suggestions:

  • first off, what you are asking about is mostly related to your research and not so much about LAMMPS. Thus you should rather consult with experts on that kind of research topic and less so post in a forum focused about using LAMMPS.
  • trying to insert molecules without issues into a dense system is notoriously difficult and requires special functionality to identify suitable insertion locations and avoid extreme forces from insertions. Fix gcmc was conceived to study diffusion of dilute gasses through porous system, so it has none of that
  • you are looking at many options to “fix” issues after the fact, but I don’t see any mention that you have validated your force field. Is the TraPPE CO2 molecule parameterization compatible with your MOF force field? Can you run simulations of systems with just some statically inserted CO2 molecules reliably?
  • reducing the overlap_cutoff parameter seems like the wrong choice to me. If anything, you might need to increase" it to make sure there is sufficient space for a CO2 molecule,
  • to gain experience with these kinds of calculations, it might be beneficial to start with a simulation that uses a noble gas instead of CO2 and learn from what you need to do.

Furthermore,

I think this is an incorrect interpretation. How much CO2 is inserted is predominantly controlled by the chemical potential and even small changes there can have large effects.

You can only find that out by studying the published literature or talking to experts. It sounds like an ugly hack to me.

3

I think there is a possibility that the overlap_cutoff is too large, that even in regular MD the minimal distance between molecules is lower than that. Since MC trials are generated from the MD trajectory, the MC trials would also include molecules closer than overlap_cutoff (inherited from MD), leading to immediate rejections.

Of course, it’s also possible that the overlap_cutoff is too small, as atoms too close to each other can lead to very large energies. You need to find this out by yourself.

1 Like