Skip to content

LAMMPS DumpReader does not return the correct coordinates for unwrapped trajectories. #3358

@hmacdope

Description

@hmacdope

Expected behavior

Unscaled coordinates written with LAMMPS in a lammpstrj dump format should be parsed correctly.

Actual behavior

The DumpReader Reader does not check what convention coordinates are written in. (xs,ys,zs) vs (xu,yu,zu). The transform distances.transform_StoR(ts.positions, ts.dimensions) is also called in _read_next_timestep regardless of whether the coordinates are scaled or not, resulting in insanely large coordinates if unscaled (xu,yu,zu) coordinates are input.

Code to reproduce the behavior

import MDAnalysis as mda
import numpy as np
from numpy import testing

n_atoms = 1536
# do it by reading the file by hand
def parse_lammpstrj(file, n_atoms):
    with open(file, "r") as f:
        lines = f.readlines()
        atomsegs = []
        for i, line in enumerate(lines):
            if "ITEM: ATOMS" in line:
                atomseg = lines[i+1:i+n_atoms+1]
                atomsegs.append(atomseg)
        coords = np.empty([len(atomsegs),n_atoms,3])
        for frame in range(len(atomsegs)):
            atom_frame = atomsegs[frame]
            for atom in range(len(atom_frame)):
                crd_line = atom_frame[atom]
                idx, _, xu, yu, zu = crd_line.split()
                xu = float(xu)
                yu = float(yu)
                zu = float(zu)
                coords[frame, atom,:] = [xu,yu,zu]
        return coords

coords = parse_lammpstrj("wat_unwraped.lammpsdump",n_atoms)

# do it with MDA
u = mda.Universe("wat_unwraped.lammpsdump")
coords_mda = np.empty([len(u.trajectory),n_atoms,3])
for i,ts in enumerate(u.trajectory):
    coords_mda[i,:,:] = ts.positions

np.testing.assert_allclose(coords_mda,coords)
....

Current version of MDAnalysis

  • Which version are you using: 2.0.0-dev0
  • Which version of Python: 3.7.5
  • Which operating system: Ubuntu Eoan Ermine

Metadata

Metadata

Assignees

Type

No type

Projects

No projects

Milestone

No milestone

Relationships

None yet

Development

No branches or pull requests

Issue actions