Tracking periodic boundary crossings

[bbye98]

Hello!

I'd like to get some clarification on how OpenMM handles boundary crossings with periodic boundary conditions.

For reference, when dumping a trajectory in LAMMPS, you can request the unwrapped particle positions xu, yu, and zu using dump_modify unwrap yes. The logic behind this is simple: LAMMPS keeps track of both wrapped particle positions (x, y, z) and the number of times the particles have crossed the boundaries (ix, iy, iz). For example, in the x-direction, if the particle crosses the left x = 0 (resp. x = L_x) boundary to jump from the left (resp. right) side of the simulation box to the right (resp. left) side, then 1 is subtracted from (resp. added to) ix. The unwrapped positions are then assumedly calculated using just xu = x + ix * dimx, yu = y + iy * dimy, etc., where dimx is the box length in the x-direction and so on.

I need the unwrapped positions because I am trying to calculate diffusivity/transport coefficients. My trajectories currently have only the wrapped positions (easy to tell from visualizing in OVITO since the particles are contained within the simulation box). I am manually unwrapping them in post-processing by checking if a particle has travelled more than half of the box size, but this is not a fool-proof approach since it may be possible for my particles to actually travel more than half of the box size with a large enough timestep and a small enough output frequency.

I was hoping OpenMM would report the number of periodic boundary crossings, but I don't believe this is a feature yet. The enforcePeriodicBox argument for the reporters and Simulation.Context.getState() only "unwraps" molecules/bonded particles on the boundary, if my understanding is correct (I saw this issue).

Is there an easy way to get the real unwrapped particle positions the way that LAMMPS currently has it implemented?

Thank you in advance for any feedback or suggestions!

[peastman]

See the FAQ. I think it will clarify the matter for you. Positions never get "unwrapped". They are what they are. Every particle has a position, and that doesn't need to have any connection to the box dimensions. When retrieving positions with getState(), you can optionally tell it to wrap the positions and bring all the molecules into a single periodic box. If you don't do that, you get the positions with no wrapping applied, whatever they happen to be.

[bbye98]

I read the section on how periodic boundary conditions work in the FAQ, and I concur--particle positions are what they are.

However, let's consider a simple toy example. Suppose we have a box that is 10x10x10 large (irrelevant units). I have two (bonded) particles at (0, 0, 0) and (1, 0, 0).

Suppose they both move to the "left" by 0.5. This means their new "real" positions are (-0.5, 0, 0) and (0.5, 0, 0).

Presumably, if I specify enforcePeriodicBox=True for the reporter, the outputted coordinates will be (-0.5, 0, 0) and (0.5, 0, 0) because the particles are bonded and the whole molecule is shifted so that they are in the same "image". This is fine and what I want.

If I specify enforcePeriodicBox=False, then I will probably get (9.5, 0, 0) and (0.5, 0, 0). Still fine since the first particle moved by well over half the box x-length, which means that I can tell that its position needs to be unwrapped from (9.5, 0, 0) to (-0.5, 0, 0). The post-processing code I've written takes care of this.

However, that's for a single timestep. Let's say I output to a trajectory file infrequently, say every 12 timesteps, with the particle travelling directly left in the x-direction at the same velocity v=-0.5. Then, the final "real" unwrapped positions should be (-6, 0, 0) and (-5, 0, 0).

In this case, regardless of what enforcePeriodicBox is, the outputted coordinates will be (4, 0, 0) and (5, 0, 0). Taking these coordinates at face value would suggest that the particles have moved "right" by 4 from the initial positions, instead of moving "left" by 6. This will give me the wrong MSD/diffusion coefficient if I calculate it using the coordinates in my trajectory (4^2/dt vs. (-6)^2/dt).

I am looking for a way to somehow recover (-6, 0, 0) and (-5, 0, 0) so that I get the correct displacement, but I'm not sure if it's possible or I'm missing something simple. I believe this is what the LAMMPS dump_modify unwrap yes does. Effectively, I'm interested in the "correct" displacement of a single particle, not between a pair of particles. I am aware that GROMACS suffers from the same issue, even with trjconv -pbc nojump, but I'm hoping to find a quick fix for OpenMM.

Thanks as always for the fast response!

[peastman]

If I specify enforcePeriodicBox=False, then I will probably get (9.5, 0, 0) and (0.5, 0, 0).

No, you will get (-0.5, 0, 0) and (0.5, 0, 0). The first particle started at (0, 0, 0) and moved by -0.5 along the x axis. That puts it at -0.5, not 9.5. After 12 steps of that size it will be at (-6, 0, 0). That's the position it will return. Since you specify enforcePeriodicBox=False it won't do anything to move the particle into a particular periodic box. It will just return the position.

Remember this important point from the FAQ: periodic boundary conditions are applied to displacements, not positions. Unless you specifically tell it to wrap the positions (with enforcePeriodicBox=True), the presence of periodic boundary conditions has no impact on the position of any particle. It can still have any position.

[bbye98]

Interesting, looks like my understanding of how enforcePeriodicBox works was incorrect, and enforcePeriodicBox=False indeed is the behavior that I was looking for (giving me "raw" "unwrapped" positions always, regardless of connectivity).

Super appreciate your help as always. Keep up the awesome work--looking forward to the release of OpenMM 8.0!