SHAKE clusters - one central atom with four bonds

Hi all-

I am simulating the diffusion of gas molecules in microporous polymers using MD in LAMMPS. I would like to use the SHAKE algorithm to constrain the bond lengths in my gas molecules so that I can transition between MD in LAMMPS and GCMC in Cassandra, which requires fixed bond lengths.

After attempting use SHAKE to constrain only the bonds on my 5-site methane model with the following command,

fix gasbonds gas shake 0.0001 20 10000 b 15

where 15 is the bond type I care about constraining, get an error stating that clusters of only 4 atoms max are allowed.

In this post from last year, referencing a similar 5-site molecule, Axel helpfully noted that SHAKE can be used for the exact kind of system I have:

1 central atom with 1, 2, 3, or 4 bonds connected to it that are constrained

However, I cannot figure out how to actually do this, and it seems to contradict the documentation.

I figured I would ask, before attempting to extend fix_shake.cpp to handle 5-site clusters, whether anyone has done this using the most recent (February 4 2025) version of LAMMPS.

Yes, so I was wrong. Happens sometimes.

Note the alternative fix rigid/* that handles clusters of any size.

@akohlmey Axel, no worries at all - I would wager a guess that I am wrong about LAMMPS far more often that you, which was why I wanted to check.

@srtee Shern, I had tried using fix/rigid/small/nvt as my first line of defense, but have had issues trying to validate my dynamics against experiments - observing energy leakage when testing with fix nve and needing to use a smaller timestep being two of those issues. I suspect that using fix momentum in my simulations, even every 1000 timesteps, may also be related. Of course this is all part of the debugging process.

Because the model I am using was parametrized as a flexible one, I am going to try seeing if just doing it that way will work. Perhaps I should have started there.

Appreciate the feedback from both of you.