Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

WIP: Lammpsdump velocities (and other fields) #3448

Closed
wants to merge 5 commits into from
Closed

WIP: Lammpsdump velocities (and other fields) #3448

Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
387f0ad
Parse velocities from LAMMPSDUMP
schlaicha Oct 26, 2021
6f86e50
Add unit conversion for LAMMPS real units
schlaicha Oct 29, 2021
10723cf
Correct frame numbering
schlaicha Oct 29, 2021
c130c18
Adjust documentation
schlaicha Oct 29, 2021
a985b43
PEP8
schlaicha Oct 29, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
25 changes: 22 additions & 3 deletions package/MDAnalysis/coordinates/LAMMPS.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -55,8 +55,9 @@
Dump files
----------

The DumpReader expects ascii dump files written with the default
`LAMMPS dump format`_ of 'atom'
The DumpReader expects ascii dump files. Currently supported are
coordinates (dump style atom) and velocities for custom dump styles.
Other columns are ignored.


Example: Loading a LAMMPS simulation
Expand Down Expand Up @@ -462,6 +463,7 @@ class DumpReader(base.ReaderBase):
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.
Velocities (vx, vy, vz) can be parsed, too.


.. versionchanged:: 2.0.0
Expand Down Expand Up @@ -502,7 +504,7 @@ def _reopen(self):
self.close()
self._file = util.anyopen(self.filename)
self.ts = self._Timestep(self.n_atoms, **self._ts_kwargs)
self.ts.frame = -1
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this is deliberately -1 as it gets incremented in the loading process

self.ts.frame = 0

@property
@cached('n_atoms')
Expand Down Expand Up @@ -606,15 +608,32 @@ def _read_next_timestep(self):

coord_cols = convention_to_col_ix[self.lammps_coordinate_convention]

# pass possible additional fields from LAMMPS dump
# pass velocities
velocity_col_names = ["vx", "vy", "vz"]
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would make this a class attribute, much like self._conventions is, as the data is not mutable or instance dependent. :)

if any(item in attrs for item in velocity_col_names):
for cv_name, cv_col_names in enumerate(velocity_col_names):
try:
velocity_cols = [attr_to_col_ix[x] for x in velocity_col_names]
ts.has_velocities = True
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This works, but I think we only want to set ts.has_velocities = True on exit from this loop, slightly cleaner IMO.

except KeyError:
raise ValueError("Velocities detected in LAMMPS Dump file but MDAnalysis can only deal with all entries vx, vy, vz present")

ids = "id" in attr_to_col_ix
for i in range(self.n_atoms):
fields = f.readline().split()
if ids:
indices[i] = fields[attr_to_col_ix["id"]]
ts.positions[i] = [fields[dim] for dim in coord_cols]
if ts.has_velocities:
ts.velocities[i] = [fields[dim] for dim in velocity_cols]

order = np.argsort(indices)
ts.positions = ts.positions[order]
if ts.has_velocities:
ts.velocities = ts.velocities[order]
# LAMMPS reader only supports real units
ts.velocities *= units.speedUnit_factor['Angstrom/femtosecond']
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Query @richardjgowers on units? I'm no lammps expert, what are the range of possible units?

if (self.lammps_coordinate_convention.startswith("scaled")):
# if coordinates are given in scaled format, undo that
ts.positions = distances.transform_StoR(ts.positions,
Expand Down