Fix rigid does not keep bond lengths constant


I am testing to simulate SWM4-ndp water molecules using fix rigid, as in the examples/PACKAGES/drude/swm4-ndp/in.swm4-ndp.nh file.

I have used fix rigid/npt instead of fix rigid/npt/small in the example because the simulation would not run with fix rigid/npt/small. I think this problem was reported as #3313 on Github.

I have simply added a compute to calculate the bond length between O-M and O-H in the water molecules. It seems that the bond lengths are drifting slightly. Over longer time (millions of time steps) they can drift quite significantly, up to 0.3 Å (I have found). I suspect there is something wrong here? At least I have found the properties of the SWM4-ndp water model are a bit strange, and not what they are supposed to be. Here are O-H bond lengths during this short simulation: OH_bond.dump (18.9 KB)

Here is the input file I have used:
in.example (1.2 KB)

Should I create a bug report on Github for this?

If I could use fix shake instead of fix rigid I would do this, but it seems the simulation is unstable when I try to make the SWM4-ndp water molecules rigid using fix shake.

It seems that the water molecules remain rigid (constant bond lengths and angle) when using fix langevin/drude instead of the Nosé-Hoover termostat. So that is a way of avoiding this issue.

Since you have a working example with fix langevin, I would consider this a potential bug in fix rigid/npt worth reporting. If there is a page in the documentation which either says fix rigid/npt should do something different or does not say that it does what it does here, we could at the very least update documentation while having a look at the code.

1 Like