fix gcmc does not work correctly with MACE/MLIAP

[egor-demidov]

Summary

When attempting to run a GCMC simulation with MACE potential implemented via pair_style mliap, deletion of molecules is strongly favored over insertion. Upon investigation, it was revealed that fix gcmc does not delete atoms from the system when attempting a deletion, but rather it adds them to an "exclusion group" to calculate the potential energy of the system without them. Only if deletion is accepted do the atoms actually get deleted. While this approach works for regular pair_styles, the exclusion group is ignored when running gcmc with pair_style mliap. The apparent potential energy of the system without with or without the proposed deletion is the same. Since a molecule also has it internal potential energy $E_{\rm intra}$, the potential energy change of deletion $\Delta E$ becomes:

$$E_{\rm before} = E_{\rm after}$$ $$\Delta E = E_{\rm after} - E_{\rm before} + E_{\rm intra}$$ $$\Delta E = E_{\rm intra}$$

Because $E_{\rm intra}$ is typically a highly negative value, deletion of molecules becomes unphysically favorable over insertion.

Expected Behavior

If molecules in the box are interacting with each other, the potential energy change of deletion $\Delta E$ should be equal to the potential energy cost of breaking the interactions. If the interactions are attractive, $\Delta E>0$. If they are repulsive, $\Delta E<0$. If the molecules do not interact and the substance is an ideal gas, $\Delta E=0$.

LAMMPS Version and Platform

22 Jul 2025 - Update 4
Linux "Ubuntu 24.04.4 LTS" 6.17.0-35-generic x86_64

Full output from lmp -h is attached in lammps_version.txt

Steps to Reproduce

Path source file "src/MC/fix_gcmc.cpp" to modify the method FixGCMC::attempt_molecule_deletion_full() to output energy before and after each attempted deletion to a text file "gcmc_delete_debug.dat" (fix_gcmc.patch / fix_gcmc.cpp). Run a very short GCMC simulation trying to insert water molecules into a box with ~50 exchange attempts. Check "gcmc_delete_debug.dat" to see that Eafter is equal to Ebefore for each attempted deletion. An example input script and its corresponding output are attached in bug-preproduction.tar.gz. The potential used in this example can be downloaded from MACE foundation model repository.

[Luthaf]

Writing down some findings since this issue prompted me to look deeper into how fix_gcmc works and how to handle it with MLIPs.

As noted above, fix_gcmc does not actually remove the atoms but instead move them out of the "all" group. This is then handled by the neighbor list builder in LAMMPS, which ignores any pair between an atom in "all", so pure pair potential do the right thing and compute the energy without the contributions from the molecule. But most/all MLIPs also have some kind of per-atom baseline, which is not removed through this and likely the root cause of this deletion bias.

In our case (not pair_style mliap), this was fixed by checking atom->mask[i] & 1 to only include atoms that are in the "all" group in the calculation metatensor/lammps@43a023c. We do have a MACE interface as well, but it is still very experimental. Feel free to contact me if you would like to give it a try.

Regardless, this looks like something that should be fixed directly in pair_style MLIAP, and not inside the MACE code.

[bernstei]

@Luthaf I can't quite decide if what you're saying is actually compatible with the original description: " The apparent potential energy of the system without with or without the proposed deletion is the same.", assuming @egor-demidov means this precisely. I agree that turning off interactions without removing the isolated atom energy would make the current GCMC implementation fail for MLIPs (because of not getting rid of the isolated atom energy), but it wouldn't in general lead to the same energy (unless you're in the dilute limit). @egor-demidov could this be what's happening?

Also, does anyone happen to know how widely tests GCMC is for EAMs? I guess the statement about full_energy and EAMs in the docs suggests that the developers think that it does work. However, with EAMs the energy for an isolated atom (E_density(rho = 0) + (E_pair = 0) doesn't have to be 0, although perhaps by convention it usually, or always, is.

[Luthaf]

I did not try the patch from the initial post, but unless the MLIP re-computes its own NL and ignore the LAMMPS one, I don't see how the energy could contain anything else than the per-atom baseline?

[egor-demidov]

Here is the energy log from fix gcmc for the first deletion attempt with MACE-MLIAP:
DEL step 1 natoms_trial 15 ... Ebefore_stored -71.188074411111359 Ebefore_fresh -71.188074411111359 Eafter -71.188074411111359 Eintra -14.237614882222299 ... accepted 1

For comparison, here is the energy log for the fist deletion of a fictional molecule with lj/cut/coul/cut interactions with special_bonds lj/coul 1.0 1.0 1.0 and full_energy fix gcmc:
DEL step 1 natoms_trial 15 ngas_before 15 ... Ebefore_stored -44.048265994452862 Ebefore_fresh -44.048265994452862 Eafter -35.238612795562297 Eintra -8.8096531988905387 ... accepted 0

In the first case, the deletion attempt results in the same Ebefore and Eafter. In the second case, everything works as expected.

test_lj.tar.gz

[egor-demidov]

@Luthaf @bernstei This is entirely speculative, but is it possible that the failure point is conversion from LAMMPS neighbor list to MLIAP neighbor data? LAMMPS neighbor lists enumerate pairs of atoms with separation less than the cutoff distance. What if excluded pairs do get removed from neighbor lists, but excluded atoms remain as central atoms? And maybe MLIAP potentials do not rely on pairs in LAMMPS neighbor lists, but instead use central atoms to build atomic environments and evaluate energies?