Fixed bugs in XYZ caused by old syntax and updated nn searches to support atomic labels

Author: avmarchenko

.gitignore

@@ -25,6 +25,11 @@ var/*
 .installed.cfg
 *.egg
 *.pypirc
+win-32
+win-64
+linux-32
+linux-64
+osx-64
 
 # PyInstaller
 #  Usually these files are written by a python script from a template

exatomic/__init__.py

@@ -18,7 +18,7 @@
 
 .. _atomic: https://en.wikipedia.org/wiki/Atomic_units
 """
-__exatomic_version__ = (0, 3, 2)
+__exatomic_version__ = (0, 3, 3)
 __version__ = '.'.join((str(v) for v in __exatomic_version__))
 
 

exatomic/algorithms/neighbors.py

@@ -65,17 +65,21 @@ def _slice_atoms_molecules(universe, sources, restrictions, n):
     Initial check of the unvierse data and argument types and creation of atom
     and molecule table slices.
     """
-    if 'classification' not in universe.molecule.columns and any(len(source) > 3 for source in sources):
-        raise KeyErrror("Column 'classification' not in the molecule table, please classify molecules or select by symbols only.")
     if not isinstance(sources, list):
         sources = [sources]
     if not isinstance(restrictions, list) and restrictions is not None:
         restrictions = [restrictions]
     if isinstance(n, (int, np.int32, np.int64)):
         n = [n]
-    symbols = universe.atom['symbol'].unique()
-    classification = universe.molecule['classification'].unique()
-    if all(source in symbols for source in sources):
+    labels = universe.atom.get_atom_labels()
+    universe.atom['label'] = labels
+    symbols = universe.atom['symbol']
+    classification = universe.molecule['classification']
+    if all(source in labels for source in sources):
+        source_atoms = universe.atom[universe.atom['label'].isin(sources)]
+        mdx = source_atoms['molecule'].astype(np.int64)
+        source_molecules = universe.molecule[universe.molecule.index.isin(mdx)]
+    elif all(source in symbols for source in sources):
         source_atoms = universe.atom[universe.atom['symbol'].isin(sources)]
         mdx = source_atoms['molecule'].astype(np.int64)
         source_molecules = universe.molecule[universe.molecule.index.isin(mdx)]
@@ -85,13 +89,21 @@ def _slice_atoms_molecules(universe, sources, restrictions, n):
     else:
         classif = [source for source in sources if source in classification]
         syms = [source for source in sources if source in symbols]
+        lbls = [source for source in sources if source in labels]
         source_molecules = universe.molecule[universe.molecule['classification'].isin(classif)]
         source_atoms = universe.atom[universe.atom['molecule'].isin(source_molecules.index)]
-        source_atoms = source_atoms[source_atoms['symbol'].isin(syms)]
+        if len(syms) > 0:
+            source_atoms = source_atoms[source_atoms['symbol'].isin(syms)]
+        if len(lbls) > 0:
+            source_atoms = source_atoms[source_atoms['label'].isin(lbls)]
     other_molecules = universe.molecule[~universe.molecule.index.isin(source_molecules.index)]
     other_atoms = universe.atom[~universe.atom.index.isin(source_atoms.index)]
     if restrictions is not None:
-        if all(other in symbols for other in restrictions):
+        if all(other in labels for other in restrictions):
+            other_atoms = other_atoms[other_atoms['label'].isin(restrictions)]
+            mdx = other_atoms['molecule'].astype(np.int64)
+            other_molecules = other_molecules[other_molecules.index.isin(mdx)]
+        elif all(other in symbols for other in restrictions):
             other_atoms = other_atoms[other_atoms['symbol'].isin(restrictions)]
             mdx = other_atoms['molecule'].astype(np.int64)
             other_molecules = other_molecules[other_molecules.index.isin(mdx)]
@@ -101,9 +113,14 @@ def _slice_atoms_molecules(universe, sources, restrictions, n):
         else:
             classif = [other for other in restrictions if other in classification]
             syms = [other for other in restrictions if other in symbols]
+            lbls = [other for other in restrictions if other in labels]
             other_molecules = other_molecules[other_molecules['classification'].isin(classif)]
             other_atoms = other_atoms[other_atoms['molecule'].isin(other_molecules.index)]
-            other_atoms = other_atoms[other_atoms['symbol'].isin(syms)]
+            if len(syms) > 0:
+                other_atoms = other_atoms[other_atoms['symbol'].isin(syms)]
+            if len(lbls) > 0:
+                other_atoms = other_atoms[other_atoms['label'].isin(lbls)]
+    del universe.atom['label']
     return source_atoms, other_atoms, source_molecules, other_molecules, n
 
 

exatomic/algorithms/pcf.py

@@ -1,16 +1,18 @@
 # -*- coding: utf-8 -*-
-'''
+# Copyright (c) 2015-2016, Exa Analytics Development Team
+# Distributed under the terms of the Apache License 2.0
+"""
 Pair Correlation Functions
 ############################
-'''
+"""
 import numpy as np
 import pandas as pd
 from exatomic import Length
 
 
 def radial_pair_correlation(universe, a, b, dr=0.05, start=1.0, stop=13.0,
-                            length='A', window=1):
-    '''
+                            length="A", window=1):
+    """
     Compute the angularly independent pair correlation function.
 
     This function is sometimes called the pair radial distribution function. The
@@ -24,8 +26,8 @@ def radial_pair_correlation(universe, a, b, dr=0.05, start=1.0, stop=13.0,
 
     .. code-block:: Python
 
-        pcf = radial_pair_correlation(universe, 'O', 'O')
-        pcf.plot(secondary_y='Pair Count')
+        pcf = radial_pair_correlation(universe, "O", "O")
+        pcf.plot(secondary_y="Pair Count")
 
     .. math::
 
@@ -57,30 +59,30 @@ def radial_pair_correlation(universe, a, b, dr=0.05, start=1.0, stop=13.0,
         the volume sampled during computation of two body properties divided by
         the number of properties used in the histogram (the triple summation
         above, divided by the normalization for the radial distance outward).
-    '''
+    """
     bins = np.arange(start, stop, dr)                     # Discrete values of r for histogram
-    symbol = universe.atom['symbol'].astype(str)          # To select distances, map to symbols
-    symbol0 = universe.atom_two['atom0'].map(symbol)
-    symbol1 = universe.atom_two['atom1'].map(symbol)
+    symbol = universe.atom["symbol"].astype(str)          # To select distances, map to symbols
+    symbol0 = universe.atom_two["atom0"].map(symbol)
+    symbol1 = universe.atom_two["atom1"].map(symbol)
     symbols = symbol0 + symbol1
     indexes = symbols[symbols.isin([a + b, b + a])].index # Distances of interest or those that
-    distances = universe.atom_two.ix[indexes, 'distance'] # match symbol pairs
+    distances = universe.atom_two.ix[indexes, "distance"] # match symbol pairs
     hist, bins = np.histogram(distances, bins)            # Compute histogram
     nn = hist.sum()                                       # Number of observations
     bmax = bins.max()                                     # Note that bins is unchanged by np.hist..
-    rx, ry, rz = universe.frame[['rx', 'ry', 'rz']].mean().values
+    rx, ry, rz = universe.frame[["rx", "ry", "rz"]].mean().values
     ratio = (((bmax/rx + bmax/ry + bmax/rz)/3)**3).mean() # Variable actual vol and bin vol
     v_shell = bins[1:]**3 - bins[:-1]**3                  # Volume of each bin shell
-    v_cell = universe.frame['cell_volume'].mean()         # Actual volume
+    v_cell = universe.frame["cell_volume"].mean()         # Actual volume
     g = hist*v_cell*ratio/(v_shell*nn)                    # Compute pair correlation
-    na = universe.atom[universe.atom['symbol'] == a].groupby('frame').size().mean()
-    nb = universe.atom[universe.atom['symbol'] == b].groupby('frame').size().mean()
+    na = universe.atom[universe.atom["symbol"] == a].groupby("frame").size().mean()
+    nb = universe.atom[universe.atom["symbol"] == b].groupby("frame").size().mean()
     if a == b:
         nb -= 1
     n = hist.cumsum()/nn*na*nb                           # Compute pair count
-    r = (bins[1:] + bins[:-1])/2*Length['au', length]
-    unit = 'au'
-    if length in ['A', 'angstrom', 'ang']:
+    r = (bins[1:] + bins[:-1])/2*Length["au", length]
+    unit = "au"
+    if length in ["A", "angstrom", "ang"]:
         unit = r"\AA"
     rlabel = r"$r\ \mathrm{(" + unit + ")}$"
     glabel = r"$g_\mathrm{" + a + b + r"}(r)$"

exatomic/filetypes/xyz.py

@@ -75,7 +75,7 @@ def write(self, path, trajectory=True, float_format='%    .8f'):
             with open(path, 'w') as f:
                 f.write(str(self))
         else:
-            grps = self.atom.grouped()
+            grps = self.atom.cardinal_groupby()
             n = len(str(self.frame.index.max()))
             for frame, atom in grps:
                 filename = str(frame).zfill(n) + '.xyz'
@@ -90,16 +90,25 @@ def write(self, path, trajectory=True, float_format='%    .8f'):
                              quoting=csv.QUOTE_NONE, escapechar=' ')
 
     @classmethod
-    def from_universe(cls, universe, float_format='%    .8f'):
+    def from_universe(cls, universe, atom_table='atom', float_format='%    .8f'):
         """
         Create an xyz file editor from a given universe. If the universe has
         more than one frame, creates an xyz trajectory format editor.
+
+        Args:
+            universe: The universe
+            atom_table (str): One of 'atom', 'unit', or 'visual' corresponding to coordinates
+            float_format (str): Floating point format (for writing)
         """
         string = ''
-        grps = universe.atom.grouped()
+        grps = universe.atom.cardinal_groupby()
         for frame, atom in grps:
             string += cls._header.format(nat=len(atom), comment='frame: ' + str(frame))
             atom_copy = atom[cls._cols].copy()
+            if atom_table == 'unit':
+                atom_copy.update(universe.unit_atom)
+            elif atom_table == 'visual':
+                atom_copy.update(universe.visual_atom)
             atom_copy['x'] *= Length['au', 'A']
             atom_copy['y'] *= Length['au', 'A']
             atom_copy['z'] *= Length['au', 'A']

meta.yaml

@@ -23,6 +23,7 @@ requirements:
         - ipyparallel
         - sphinx
         - nose
+        - numba
         - exa
 
     run:
@@ -41,6 +42,7 @@ requirements:
         - ipyparallel
         - sphinx
         - nose
+        - numba
         - exa
 
 about: