Does OpenMM automatically image the molecule's positions?

[She Zhang]

Hi All,

I am trying to run a simulation with periodic boundary conditions. I have nonbonded forces in the system with a cutoff of 10 angstroms. I would like to obtain unwrapped positions of the atoms for some downstream analyses, but I found that the positions obtained from my simulation seem to be automatically wrapped in both trajectory and state files. I turned off enforcePeriodicBox in the reporter (which didn't work), and saveState has the option, enforcePeriodicBox, default to False. If I understand it correctly, getState/getPosition should return the "raw" positions without wrapping them when enforcePeriodicBox=False.

At this point, I am suspecting the simulation itself is wrapping the positions. Does anyone know if this is the case? If so, is there a way to turn off the wrapping while still apply the periodic boundary conditions when applying the forces?

I am using OpenMM 7.5.1 on Linux.

Thanks a lot in advance,
She

[Peter Eastman]

If you're running a constant pressure simulation, the barostat will wrap the molecules into a single periodic box every time it rescales the box.

[She Zhang]

If you're running a constant pressure simulation, the barostat will wrap the molecules into a single periodic box every time it rescales the box.

Thanks a lot for the reply. My simulation is indeed NPT. Is there an option to turn off the wrapping without turning off PBC altogether in barostat, by any chance?

[Peter Eastman]

No, it's built into the barostat algorithm. Sorry!

[Emilio Galllicchio]

Can you explain why the barostat algorithm requires the wrapping? In principle, the shape and volume of the box are inputs of the potential energy function, not a property of the configuration of the molecules.
Thanks
Emilio

[Peter Eastman]

That's just how the code is written. It might be possible to restructure the algorithm in a way that wouldn't require it.

[She Zhang]

No, it's built into the barostat algorithm. Sorry!

No worries. Thanks for confirming it. I can try to unwrap the trajectories in post-analyses, but it would be nice if I can obtain the unwrapped trajectories from OpenMM directly.

Thanks a lot for the prompt response too!

[Emilio Galllicchio]

That's just how the code is written. It might be possible to restructure the algorithm in a way that wouldn't require it.

We have restraint potentials implemented using CustomCentroidBondForce to keep two molecules near a certain vector displacement (offx, offy, offz) based on for example:

d12 = sqrt((x1 - offx - x2)^2 + (y1 - offy - y2)^2 + (z1 - offz - z2)^2 )

which I think would break if the two molecules are wrapped at two positions more than the minimum image distance away.

Is there you think a workaround based on periodicdistance() in this case?

(It seems to me that one would need a function that returns the minimum image vector displacement.)

[Peter Eastman]

If there's a CustomCentroidBondForce (or any other bonded force) connecting them, it will consider them to be a single molecule and always wrap them as a unit.

[Emilio Galllicchio]

If there's a CustomCentroidBondForce (or any other bonded force) connecting them, it will consider them to be a single molecule and always wrap them as a unit.

Great. It makes sense. Thanks!