Skip to content

Commit

Permalink
Add simple xyz reader (#423)
Browse files Browse the repository at this point in the history
* Add simple xyz reader

* Remove old code

* Small pep8 fix

* Add io test

* Allow python2 compatability

* Various fixes

* Rework test with topology

* Clean up test

* Add a test for the case of an incorrect number of atoms

* Attempt to handle the case of an incorrect number of atoms

* Handle and test two cases of incorrect n_atoms

* Switch to raising MBuildError

* Add xyz to manifest
  • Loading branch information
mattwthompson authored and summeraz committed May 29, 2018
1 parent 5ebdde0 commit 1468900
Show file tree
Hide file tree
Showing 6 changed files with 120 additions and 1 deletion.
2 changes: 1 addition & 1 deletion MANIFEST.in
@@ -1 +1 @@
global-include *.pdb *.mol2 *.smi *.xml
global-include *.pdb *.mol2 *.smi *.xml *.xyz
7 changes: 7 additions & 0 deletions mbuild/compound.py
Expand Up @@ -22,6 +22,7 @@
from mbuild.bond_graph import BondGraph
from mbuild.box import Box
from mbuild.exceptions import MBuildError
from mbuild.formats.xyz import read_xyz
from mbuild.formats.hoomdxml import write_hoomdxml
from mbuild.formats.lammpsdata import write_lammpsdata
from mbuild.formats.gsdwriter import write_gsd
Expand Down Expand Up @@ -78,6 +79,12 @@ def load(filename, relative_to_module=None, compound=None, coords_only=False,
if compound is None:
compound = Compound()

# Handle the case of a xyz file, which must use an internal reader
extension = os.path.splitext(filename)[-1]
if extension == '.xyz' and not 'top' in kwargs:
compound = read_xyz(filename)
return compound

if use_parmed:
warn("use_parmed set to True. Bonds may be inferred from inter-particle "
"distances and standard residue templates!")
Expand Down
62 changes: 62 additions & 0 deletions mbuild/formats/xyz.py
@@ -0,0 +1,62 @@
import numpy as np

import mbuild as mb
from mbuild.exceptions import MBuildError

__all__ = ['read_xyz']


def read_xyz(filename):
"""Read an XYZ file. The expected format is as follows:
The first line contains the number of atoms in the file The second line
contains a comment, which is not read. Remaining lines, one for each
atom in the file, include an elemental symbol followed by X, Y, and Z
coordinates in Angstroms. Columns are expected tbe separated by
whitespace. See https://openbabel.org/wiki/XYZ_(format).
Parameters
----------
filename : str
Path of the input file
Returns
-------
compound : mb.Compound
Notes
-----
The XYZ file format neglects many important details, notably as bonds,
residues, and box information.
There are some other flavors of the XYZ file format and not all are
guaranteed to be compatible with this reader. For example, the TINKER
XYZ format is not expected to be properly read.
"""

compound = mb.Compound()

with open(filename, 'r') as xyz_file:
n_atoms = int(xyz_file.readline())
xyz_file.readline()
coords = np.zeros(shape=(n_atoms, 3), dtype=np.float64)
for row, _ in enumerate(coords):
line = xyz_file.readline().split()
if not line:
msg = ('Incorrect number of lines in input file. Based on the '
'number in the first line of the file, {} rows of atoms '
'were expected, but at least one fewer was found.')
raise MBuildError(msg.format(n_atoms))
coords[row] = line[1:4]
coords[row] *= 0.1
particle = mb.Compound(pos=coords[row], name=line[0])
compound.add(particle)

# Verify we have read the last line by ensuring the next line in blank
line = xyz_file.readline().split()
if line:
msg = ('Incorrect number of lines in input file. Based on the '
'number in the first line of the file, {} rows of atoms '
'were expected, but at least one more was found.')
raise MBuildError(msg.format(n_atoms))

return compound
30 changes: 30 additions & 0 deletions mbuild/tests/test_xyz.py
@@ -0,0 +1,30 @@
import numpy as np
import pytest

import mbuild as mb
from mbuild.utils.io import get_fn
from mbuild.tests.base_test import BaseTest
from mbuild.exceptions import MBuildError


class TestXYZ(BaseTest):
def test_load_no_top(self, ethane):
ethane.save(filename='ethane.xyz')
ethane_in = mb.load('ethane.xyz')
assert len(ethane_in.children) == 8
assert ethane_in.n_bonds == 0
assert set([child.name for child in ethane_in.children]) == {'C', 'H'}

def test_load_with_top(self, ethane):
ethane.save(filename='ethane.xyz')
ethane.save(filename='ethane.mol2')
ethane_in = mb.load('ethane.xyz', top='ethane.mol2')
assert len(ethane_in.children) == 8
assert ethane_in.n_bonds == 7
assert set([child.name for child in ethane_in.children]) == {'C', 'H'}

def test_wrong_n_atoms(self):
with pytest.raises(MBuildError):
mb.load(get_fn('too_few_atoms.xyz'))
with pytest.raises(MBuildError):
mb.load(get_fn('too_many_atoms.xyz'))
10 changes: 10 additions & 0 deletions mbuild/utils/reference/too_few_atoms.xyz
@@ -0,0 +1,10 @@
10
Created with MDTraj 1.8.0, 2018-05-09
C 0.000 -1.400 -0.000
H -1.070 -1.400 0.000
H 0.357 -2.169 0.653
H 0.357 -1.581 -0.993
C 0.000 0.000 0.000
H 1.070 0.000 0.000
H -0.357 0.769 0.653
H -0.357 0.181 -0.993
10 changes: 10 additions & 0 deletions mbuild/utils/reference/too_many_atoms.xyz
@@ -0,0 +1,10 @@
6
Created with MDTraj 1.8.0, 2018-05-09
C 0.000 -1.400 -0.000
H -1.070 -1.400 0.000
H 0.357 -2.169 0.653
H 0.357 -1.581 -0.993
C 0.000 0.000 0.000
H 1.070 0.000 0.000
H -0.357 0.769 0.653
H -0.357 0.181 -0.993

0 comments on commit 1468900

Please sign in to comment.