Incorrect unwrapping for LAMMPS trajectories due to unused image flags for periodic boundary condition

[gitsirsha]

Is your feature request related to a problem?

MDAnalysis universe currently ignores LAMMPS image flags (ix, iy and iz) to properly unwrap wrapped coordinates according to periodic boundary condition when unwrapping coordinates from LAMMPSDUMP files. MDAnalysis solely does:

x_unwrapped = x_scaled * L + L/2 (where L is box dimension)

This works good only when atoms never crossed boundaries i.e., flag ix_i = 0 but this is incorrect unwrapping when the atoms have crossed boundaries and that leads to wrong results when computing MSD or diffusion coefficients or COM.

While comparing MSD results from the new method in PR#5066 (merged) vs our in-house code (that incorporates proper conversion accounting the image flags) there is a discrepancy.

x_unwrapped = x_scaled * L + L*ix (Correct convention with the in-house code)

We've further verified this when provided with an already pre-unwrapped LAMMPSSUMP file. The results are consistent between MDAnalysis and In-house code! Meaning that it's definitely the improper unwrapping only while creating the universe that creates the discrepancy.

Describe the solution you'd like

i. Properly take care of the periodic boundary condition image flags while unwrapping the universe
ii. Or at least for the time-being provide a warning to the user to provide unwrapped coordinates.

Describe alternatives you've considered

Pre-unwrapping the trajectory file (LAMMPSDUMP) following proper PBC flags

Additional context

Example below:

Notice the shoot up! Pre-unwrap to get correct results!

notebooks for this example can be found in the repository-link.

[mrshirts]

Let us know how we can help with any changes! I figure there are good periodic boundary tools in MD Analysis already, so it's hopefully just a matter of adding a couple of calls to those routines inside the MSD calculation.

[orbeckst]

The basic philosophy in MDA is to present the trajectory data as stored in the trajectory after conversion to standard units.

I'd say the standard approach for MSD is to "unwrap" the trajectory with the NoJump transformation on the fly (or pre-process).

Can you explain how the image flag in the dump file works?

• Is it recording the number of PBC crossings and is available for each step?
• Does the algorithm work for triclinic (non-orthogonal) unit cells?
• How would you want MDAnalysis to honor this flag, e.g., how would the user tell MDAnalysis that unwrapping ought to be happening?

[fxcoudert]

Doc on image flags is there: https://docs.lammps.org/read_data.html

LAMMPS always stores atom coordinates as values which are inside the simulation box. It also stores 3 flags which indicate which image of the simulation box (in each dimension) the atom would be in if its coordinates were unwrapped across periodic boundaries. An image flag of 0 means the atom is still inside the box when unwrapped. A value of 2 means add 2 box lengths to get the unwrapped coordinate. A value of -1 means subtract 1 box length to get the unwrapped coordinate. LAMMPS updates these flags as atoms cross periodic boundaries during the simulation. The dump command can output atom coordinates in wrapped or unwrapped form, as well as the 3 image flags.

LAMMPS stores coordinates as $(x, y, z, ix, iy, iz)$ where $(x,y,z)$ are Cartesian coordinates and are always inside the simulation box, and $(ix,iy,iz)$ is a triplet of integers, and tracks the dynamics of atoms across periodic boundaries. Taking a 1D example, if an atom starts at $(x=0.1, ix=0)$ is a box of size 10, and move 0.2 to the left, it will become $(x=9.9, ix=-1)$.

The algorithm works for all box types (in LAMMPS term, orthogonal or restricted triclinic or general triclinic), where the actual (unwrapped) position of one atom would be $(x, y, z) + ix \times \mathbf{a} + iy \times \mathbf{b} + iz \times \mathbf{c}$, and $(\mathbf{a}, \mathbf{b}, \mathbf{c})$ are the box vectors.

The only thing that needs to be checked is: LAMMPS dumps have so many options, it's possible to write with or without image flags, wrapped or unwrapped. It's unclear to me if there is a one-size-fits-all solution.

[gitsirsha]

@orbeckst Thanks for the reply. The nojump() transformation you'd suggested indeed seems to fix the problem. It just requires two lines of code to transform the universe before the other modules are used. Thank you for letting us know.

import MDAnalysis.transformations.nojump as nj
transformation = nj.NoJump()
u.trajectory.add_transformations(transformation)

Thank you @fxcoudert for the clarification on triclinic box!

[jnnsmpth]

For your last point @orbeckst, about how we want MDAnalysis to honor this flag - if a LAMMPS dump file has a ix iy iz column, then image flags need to be taken into account. If not, the analysis can proceed without unwrapping.

[fxcoudert]

Actually, doesn't the LAMMPS DumpReader (

mdanalysis/package/MDAnalysis/coordinates/LAMMPS.py

Line 534 in 5c7c480

class DumpReader(base.ReaderBase):

) already support this? In particular, the unwrap_images argument:

    lammps_coordinate_convention : str (optional) default="auto"
        Convention used in coordinates, can be one of the following according
        to the `LAMMPS documentation <https://docs.lammps.org/dump.html>`__:

         - "auto" - Detect coordinate type from file column header. If auto
           detection is used, the guessing checks whether the coordinates
           fit each convention in the order "unscaled", "scaled", "unwrapped",
           "scaled_unwrapped" and whichever set of coordinates is detected
           first will be used.
         - "scaled" - Coordinates wrapped in box and scaled by box length (see
            note below), i.e., xs, ys, zs
         - "scaled_unwrapped" - Coordinates unwrapped and scaled by box length,
           (see note below) i.e., xsu, ysu, zsu
         - "unscaled" - Coordinates wrapped in box, i.e., x, y, z
         - "unwrapped" - Coordinates unwrapped, i.e., xu, yu, zu

        If coordinates are given in the scaled coordinate convention (xs,ys,zs)
        or scaled unwrapped coordinate convention (xsu,ysu,zsu) they will
        automatically be converted from their scaled/fractional representation
        to their real values.
    unwrap_images : bool (optional) default=False
        If `True` and the dump file contains image flags, the coordinates
        will be unwrapped. See `read_data
        <https://docs.lammps.org/read_data.html>`__  in the lammps
        documentation for more information.

[orbeckst]

@gitsirsha did you try @fxcoudert ’s suggestion (I didn’t know..). You should be able to pass the unwarp_images=True argument in the Universe and the DumpReader should pick it up.