-
Notifications
You must be signed in to change notification settings - Fork 676
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.
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
Add Chemfiles as a coordinate reader/writer #1862
Changes from 1 commit
240c3d3
1862eff
65c5c86
ad5d199
b450ee3
634e326
ca33ab2
5eadec7
b4410db
93ddfd7
6f195ae
c3d5321
eaf172c
359ea4e
b1aa7d2
55e97e7
8fcdc23
0338863
85b7ac1
570e191
7f6395f
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -47,15 +47,25 @@ | |
|
||
from . import base, core | ||
|
||
import chemfiles | ||
from chemfiles import Trajectory, Frame, Atom, UnitCell, Residue, Topology | ||
try: | ||
import chemfiles | ||
except ImportError: | ||
HAS_CHEMFILES = False | ||
else: | ||
HAS_CHEMFILES = True | ||
|
||
|
||
|
||
MIN_CHEMFILES_VERSION = LooseVersion("0.9") | ||
MAX_CHEMFILES_VERSION = LooseVersion("0.10") | ||
|
||
|
||
def check_chemfiles_version(): | ||
if not HAS_CHEMFILES: | ||
raise RuntimeError( | ||
"No Chemfiles package found. " | ||
"Try installing with 'pip install chemfiles'" | ||
) | ||
version = LooseVersion(chemfiles.__version__) | ||
if version < MIN_CHEMFILES_VERSION or version >= MAX_CHEMFILES_VERSION: | ||
raise RuntimeError( | ||
|
@@ -65,8 +75,7 @@ def check_chemfiles_version(): | |
|
||
|
||
class ChemfilesReader(base.ReaderBase): | ||
""" | ||
Coordinate reader using chemfiles. | ||
"""Coordinate reader using chemfiles. | ||
|
||
The file format to used is guessed based on the file extension. If no | ||
matching format is found, a ``ChemfilesError`` is raised. It is also | ||
|
@@ -79,12 +88,13 @@ def __init__(self, filename, chemfiles_format="", **kwargs): | |
""" | ||
Parameters | ||
---------- | ||
filename : str | ||
trajectory filename | ||
filename : chemfiles.Trajectory or str | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I am not sure I understand why passing a chemfiles.Trajectory should work. What would be the use case for this? Users creating chemfiles trajectories and then wanting to do analysis with MDA? There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Yeah, it's so you can open a chemfiles thing as you want then pass this as the trajectory MDA will use. Passing a filename and specifying |
||
the chemfiles object to read or filename to read | ||
chemfiles_format : str (optional) | ||
use the given format name instead of guessing from the extension. | ||
The `list of supported formats <formats>`_ and the associated names | ||
is available in chemfiles documentation. | ||
if *filename* was a string, use the given format name instead of | ||
guessing from the extension. The `list of supported formats | ||
<formats>`_ and the associated names is available in chemfiles | ||
documentation. | ||
**kwargs : dict | ||
General reader arguments. | ||
|
||
|
@@ -99,13 +109,21 @@ def __init__(self, filename, chemfiles_format="", **kwargs): | |
self.ts = self._Timestep(self.n_atoms, **self._ts_kwargs) | ||
self.next() | ||
|
||
@staticmethod | ||
def _format_hint(thing): | ||
"""Can this Reader read *thing*""" | ||
# nb, filename strings can still get passed through if | ||
# format='CHEMFILES' is used | ||
return HAS_CHEMFILES and isinstance(thing, chemfiles.Trajectory) | ||
|
||
def _open(self): | ||
self._closed = False | ||
self._step = 0 | ||
self._frame = -1 | ||
# Open file last to ensure that all other attributes are set | ||
# in case of exception | ||
self._file = Trajectory(self.filename, 'r', self._format) | ||
if isinstance(self.filename, chemfiles.Trajectory): | ||
self._file = self.filename | ||
else: | ||
self._file = chemfiles.Trajectory(self.filename, 'r', self._format) | ||
|
||
def close(self): | ||
"""close reader""" | ||
|
@@ -224,7 +242,7 @@ def __init__(self, filename, n_atoms=0, mode="w", chemfiles_format="", topology= | |
self.n_atoms = n_atoms | ||
if mode != "a" and mode != "w": | ||
raise IOError("Expected 'a' or 'w' as mode in ChemfilesWriter") | ||
self._file = Trajectory(self.filename, mode, chemfiles_format) | ||
self._file = chemfiles.Trajectory(self.filename, mode, chemfiles_format) | ||
self._closed = False | ||
if topology is not None: | ||
if isinstance(topology, string_types): | ||
|
@@ -293,21 +311,22 @@ def _timestep_to_chemfiles(self, ts): | |
""" | ||
Convert a Timestep to a chemfiles Frame | ||
""" | ||
frame = Frame() | ||
# TODO: CONVERTERS? | ||
frame = chemfiles.Frame() | ||
frame.resize(ts.n_atoms) | ||
if ts.has_positions: | ||
frame.positions[:] = ts.positions[:] | ||
if ts.has_velocities: | ||
frame.add_velocities() | ||
frame.velocities[:] = ts.velocities[:] | ||
Luthaf marked this conversation as resolved.
Show resolved
Hide resolved
|
||
frame.cell = UnitCell(*ts.dimensions) | ||
frame.cell = chemfiles.UnitCell(*ts.dimensions) | ||
return frame | ||
|
||
def _topology_to_chemfiles(self, obj, n_atoms): | ||
""" | ||
Convert an AtomGroup or an Universe to a chemfiles Topology | ||
""" | ||
topology = Topology() | ||
topology = chemfiles.Topology() | ||
if not hasattr(obj, "atoms"): | ||
# use an empty topology | ||
topology.resize(n_atoms) | ||
|
@@ -318,7 +337,7 @@ def _topology_to_chemfiles(self, obj, n_atoms): | |
for atom in obj.atoms: | ||
name = getattr(atom, 'name', "") | ||
type = getattr(atom, 'type', name) | ||
chemfiles_atom = Atom(name, type) | ||
chemfiles_atom = chemfiles.Atom(name, type) | ||
|
||
if hasattr(atom, 'altLoc'): | ||
chemfiles_atom["altloc"] = str(atom.altLoc) | ||
|
@@ -332,7 +351,7 @@ def _topology_to_chemfiles(self, obj, n_atoms): | |
if hasattr(atom, 'resid'): | ||
resname = getattr(atom, 'resname', "") | ||
if atom.resid not in residues.keys(): | ||
residues[atom.resid] = Residue(resname, atom.resid) | ||
residues[atom.resid] = chemfiles.Residue(resname, atom.resid) | ||
residue = residues[atom.resid] | ||
|
||
atom_idx = len(topology.atoms) | ||
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@richardjgowers you want a nicer error message than this? I notice the imports should be in the classes themselves to avoid an error when importing mda without chemfiles.