simtk.org

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

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.

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?

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

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

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

peastman wrote: ↑

Tue May 11, 2021 12:33 pm

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!

peastman wrote: ↑

Tue May 11, 2021 1:06 pm

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.)

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.

peastman wrote: ↑

Thu May 13, 2021 9:00 am

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!