I am trying to use fix bond/react to model a system which has some reacting rigid molecules. I am having a hard time with fix bond/react because some of my molecules are linear.
To make the molecules rigid, I could use something like SHAKE or RATTLE or RIGID. But for linear (180 deg angle) molecules, SHAKE and RATTLE do not work, as described in the fix-shake doc page. So RIGID seems to be the only option for making molecules rigid.
When using fix bond/react, I first set up the reaction:
fix rxns all bond/react stabilization yes statted_grp .03 &
react reaction_stp1 all 1 0.0 3.0 mol.pre mol.post my.map
If I use RIGID/SMALL for time integration, with something like:
fix nvt_react statted_grp_REACT rigid/nvt/small molecule temp 300 300 100
I get an error because the group ‘statted_grp_REACT’ created by bond/react is a dynamic group (containing all the non-reaction-site particles), which it seems rigid/nvt/small cannot use (error: “Fix nvt_react does not allow use with a dynamic group”)
To recap my logic, I believe the only choice for rigid linear molecules is to use RIGID, but time integration with RIGID is not possible because the group created by bond/react is dynamic.
Is this train of logic correct? Or is there some way to simulate rigid, linear molecules using bond/react?
Thanks for your ideas!
I cannot talk much about fix bond/react, but I have some insight on the reasons for the issues you are seeing.
First off, using fix shake or fix rattle wouldn’t make a difference. The reason is that for both, shake/rattle constraints and rigid bodies, some information about the constrained/rigid objects is collected at the beginning of a run and then used throughout for efficiency reasons. That collides with applying them to dynamic groups. Instead, there are ways that fixes can tell fix shake or rigid/small that a rigid object has been added or removed (and those have to be defined as a molecule object). So that creates significant limitations.
So the only way that I can think of to implement a simulation like you are describing would be to use fix halt with some custom compute style that detects if a condition for a reaction is true. That would stop the simulation, you could manipulate your system as needed, and then continue while re-initializing the rigid bodies at the beginning of the run continuation.
But perhaps Jake has a better idea.
This issue has come up before, and you summarized it well. I briefly tried to add ‘dynamic group’ support to rigid bodies but ran into the complications Axel mentioned, and I don’t see a straightforward way to solve it from the bond/react side. Perhaps I’ll take another look eventually, but for now reactive rigid bodies and bond/react don’t mix well.
Hello again, thanks to both of you for your ideas! I had a feeling this would be the case. I’ll think through some work-arounds.