I have a question related to the use and sizing of SHAKE clusters in LAMMPS version 23 June 2022.
I am attempting to replicate a model for the polyatomic ion SO4, that uses the SHAKE algorithm to facilitate the ion’s shape. My understanding after reading the documentation with both the “fix shake” and “molecule” commands is that a shake cluster can only be up to 4 atoms, meaing that SO4, with 4 Oxygen atoms around the central Sulfur atom cannot be used. Could someone please clarify the limitations associated with LAMMPS’ SHAKE command and if there is a method that I am missing that would allow me to use SHAKE and for a bonded set with more than 4 atoms?
The problem with having connected constraints is that the system of linear equations become quickly overdefined and then has no unique solution. This already happens if you have, e.g. two bonds and an 180 degree angle.
The options you have is 1 central atom with 1, 2, 3, or 4 bonds connected to it that are constrained, or two bonds and an angle where the middle atom of the angle is the central atom of the cluster (internally, there is no angle constraint, but the system becomes a triangle of 3 bonds)
None of the bond atoms may be the central atom for another SHAKE cluster.
So, you could have a shake constraint for all for S-O bonds, but not on the O-O distances or O-S-O angles.
To keep such a construct rigid in LAMMPS, you must use fix rigid or fix rigid/small, where you would assign a unique molecule ID to each \textsf{SO}_4^{2-} anion and then use the “molecule” keyword to define the rigid clusters.
Thank you for the response. If I understand correctly, I can use I can use the SHAKE algorithm to constrain the SO4 anion with 4 bonds all connected to a central S atom if I ignore the angle restraints and also include a “fix rigid” command with a “molecule” template to further constrain the anion. Is that correct?
If I have understood, I have a follow up question pertaining the “molecule” command template formatting. The documentation for the shake flags, atoms, and bond types seems to focus on a maximum of 3 bonds when the maximum flag (4) is used.
No. Either you only constrain the S-O bonds or you use fix rigid.
You don’t need to use a molecule template unless you want to deposit sulfate anions, just use the “molecule” keyword to define which are the rigid bodies.
You haven’t and the followup question is irrelevant