Solving density divergence issues with MD + Fix GCMC + Region + Molecules ?

Dear LAMMPS users,

I would like to share with you what can possibly fix some issues one gets while running hybrid GCMD simulations using the fix gcmc with both the '’mol’’ and ‘’region’’ keywords.

I am using the 17 November 2016 version. I have performed some test simulations in which the fix gcmc is used to insert/delete molecules in a given region of the simulation box while atoms positions are integrated by means of NVT time integration (fix nvt). Only one type of molecule (CH4) is present in these simulations. Initially, the system is composed of a given number of molecules, which were equilibrated in a previous classical NVT MD simulation run.

During the GCMD simulations, I observed that the number of molecules in the box increased almost monotonically and could not reach a reasonable stationary value (even lowering the chemical potential value in some unphysical way could not help). After some other tests I found out that molecules mapped to periodic images of the simulation box were never deleted by the GCMC algorithm, resulting in the number of accepted insertion moves always superior to the number of accepted deletion moves.

Then I had a look at the source code, and more specifically to the ‘’update_gas_atom_list()’’ routine coded in the ‘’fix_gcmc.cpp’’, which creates and updates the list of molecules that fix gcmc is allowed to move or delete. It turns out that when the keyword ‘’region’’ is on, only those molecules with a center of mass (com) inside the specified region are added to this list. However the com position is computed in its unwrapped form, which explains why molecules that have moved outside the periodic cell are not added to the list if the region does not extend far enough outside the simulation box.

I might have fixed the issue by inserting the command ‘’domain->remap(com);’’ on the next line after line 2250 in the source file ‘’fix_gcmc.cpp’’. This remaps the computed com position inside the simulation box and yields a correct list of movable/deletable molecules, which seems to solve the issue of the ever increasing density during the GCMD simulations.

Is there any specific reason why the com position is left in its unwrapped form when update_gas_atom_list() performs the region->match test (in other words, will my solution trigger other problems that I did not anticipate) ? Or is this just a mistake in the source code ?

Best,
Romain

Hi Romain,

Yes, I think this is a bug.

It arises for molecules that exit the gcmc region, and subsequently enter a different periodic image of the region. I suspect that this scenario was overlooked in the original implementation of update_gas_atoms_list() for the region/molecule combo, and never showed up in testing.

For a sufficiently long simulation, the proportion of molecules in the gcmc region that are actually in a different periodic image of the region will grow larger and larger, and so the number of molecules in the region will appear to get smaller and smaller. This will correspond to a GCMC simulation with very large or even infinite chemical potential. This problem does not occur for atomic gcmc, because the group->xcm() calculation is not required.

I think your solution is exactly right. I will add it.

Aidan