Merge pull request #78 from avmarchenko/master

tjduigna · web-flow · commit 735b901d1b6e · 2016-08-12T10:18:55.000-04:00
Bugfixes for basis sets and version bump to 0.3.0
diff --git a/README.rst b/README.rst
@@ -21,6 +21,7 @@ Getting Started
 ###################
 | |docs|
 | |gitter|
+| |doi|
 Documentation can be built using `sphinx`_:
 
 .. code-block:: bash
@@ -81,3 +82,7 @@ Legal
 .. |lic| image:: http://img.shields.io/:license-apache-blue.svg?style=flat-square
     :target: http://www.apache.org/licenses/LICENSE-2.0
     :alt: License
+
+.. |doi| image:: https://zenodo.org/badge/23807/exa-analytics/exatomic.svg
+    :target: https://zenodo.org/badge/latestdoi/23807/exa-analytics/exatomic
+    :alt: DOI
diff --git a/exatomic/__init__.py b/exatomic/__init__.py
@@ -18,7 +18,7 @@
 
 .. _atomic: https://en.wikipedia.org/wiki/Atomic_units
 """
-__exatomic_version__ = (0, 2, 14)
+__exatomic_version__ = (0, 3, 0)
 __version__ = '.'.join((str(v) for v in __exatomic_version__))
 
 
diff --git a/exatomic/algorithms/basis.py b/exatomic/algorithms/basis.py
@@ -241,14 +241,14 @@ def _overlap(x1, x2, y1, y2, z1, z2, l1, l2, m1, m2, n1, n2, N1, N2, alpha1, alp
     py2 = yp - y2
     pz2 = zp - z2
     pg12 = np.sqrt(np.pi / gamma)
-    xIx = 0
-    yIy = 0
-    zIz = 0
+    xix = 0
+    yiy = 0
+    ziz = 0
     ltot = l1 + l2
     mtot = m1 + m2
     ntot = n1 + n2
     if ltot == 0:
-        xIx = pg12
+        xix = pg12
     else:
         iii = (ltot - 1) // 2 if ltot % 2 != 0 else ltot // 2
         for i in range(iii):
@@ -263,9 +263,9 @@ def _overlap(x1, x2, y1, y2, z1, z2, l1, l2, m1, m2, n1, n2, N1, N2, alpha1, alp
                 newt1 = fac(l1) / fac(j) / fac(l1 - j)
                 newt2 = fac(l2) / fac(k) / fac(l2 - k)
                 fk += newt1 * newt2 * (px1 ** (l1 - j)) * (px2 ** (l2 - k))
-            xIx += prod * fk
+            xix += prod * fk
     if mtot == 0:
-        yIy = pg12
+        yiy = pg12
     else:
         iii = (mtot - 1) // 2 if mtot % 2 != 0 else mtot // 2
         for i in range(iii):
@@ -280,9 +280,9 @@ def _overlap(x1, x2, y1, y2, z1, z2, l1, l2, m1, m2, n1, n2, N1, N2, alpha1, alp
                 newt1 = fac(m1) / fac(j) / fac(m1 - j)
                 newt2 = fac(m2) / fac(k) / fac(m2 - k)
                 fk += newt1 * newt2 * (py1 ** (m1 - j)) * (py2 ** (m2 - k))
-            yIy += prod * fk
+            yiy += prod * fk
     if ntot == 0:
-        zIz = pg12
+        ziz = pg12
     else:
         iii = (ntot - 1) // 2 if ntot % 2 != 0 else ntot // 2
         for i in range(iii):
@@ -297,9 +297,9 @@ def _overlap(x1, x2, y1, y2, z1, z2, l1, l2, m1, m2, n1, n2, N1, N2, alpha1, alp
                 newt1 = fac(n1) / fac(j) / fac(n1 - j)
                 newt2 = fac(n2) / fac(k) / fac(n2 - k)
                 fk += newt1 * newt2 * (pz1 ** (n1 - j)) * (pz2 ** (n2 - k))
-            zIz += prod * fk
+            ziz += prod * fk
     exponent = alpha1 * alpha2 * ab2 / gamma
-    s12 = N1 * N2 * np.exp(-exponent) * xIx * yIy * zIz
+    s12 = N1 * N2 * np.exp(-exponent) * xix * yiy * ziz
     return s12
 
 
diff --git a/exatomic/basis.py b/exatomic/basis.py
@@ -70,14 +70,7 @@ def shells(self):
         Returns a list of all the shells in the basis set
         """
         cols = [col.replace('bas_', '') for col in self if 'bas_' in col]
-        shells = []
-        for l in lorder:
-            for ll in cols:
-                if l == ll:
-                    shells.append(l)
-        return shells
-
-
+        return lorder[:len(cols)]
 
 
 class BasisSet(DataFrame):
@@ -151,31 +144,6 @@ class GaussianBasisSet(BasisSet):
     _precision = {'alpha': 8, 'd': 8}
     _categories = {'set': np.int64, 'L': np.int64}#, 'shell_function': np.int64}
 
-#    def _custom_traits(self):
-#        g = self.grpd
-#        alphas = g.apply(
-#                 lambda x: x.groupby('set').apply(
-#                 lambda y: y.groupby('shell_function').apply(
-#                 lambda z: z['alpha'].values).values)).to_json(orient='values')
-#        #alphas = Unicode(''.join(['[', alphas, ']'])).tag(sync=True)
-#        alphas = Unicode(alphas).tag(sync=True)
-#
-#        ds = g.apply(
-#             lambda x: x.groupby('set').apply(
-#             lambda y: y.groupby('shell_function').apply(
-#             lambda z: z['d'].values).values)).to_json(orient='values')
-#        #ds = Unicode(''.join(['[', ds, ']'])).tag(sync=True)
-#        ds = Unicode(ds).tag(sync=True)
-#
-#        ls = g.apply(
-#             lambda x: x.groupby('set').apply(
-#             lambda y: y.groupby('shell_function').apply(
-#             lambda z: z['L'].astype(np.int64).values).values)).to_json(orient='values')
-#        #ls = Unicode(''.join(['[', ls, ']'])).tag(sync=True)
-#        ls = Unicode(ls).tag(sync=True)
-#
-#        return {'gaussianbasisset_d': ds, 'gaussianbasisset_l': ls,
-#                'gaussianbasisset_alpha': alphas}
 
     def basis_count(self):
         """
@@ -191,39 +159,40 @@ def _check(self):
         if 'N' not in self.columns:
             self._normalize()
 
-#    @classmethod
-#    def expand(cls, universe, inplace=False):
-#        '''
-#        The minimum information specified by a basis set does not include
-#        expansion due to degeneracy from m_l. This will expand the basis in a
-#        systematic cartesian ordering convention to generate the full cartesian
-#        basis. The universe argument must already have a universe with atom,
-#        basis_set_summary, and gaussian_basis_set attributes.
-#        '''
-#        bases = universe.gaussian_basis_set[abs(universe.gaussian_basis_set['d']) > 0].groupby('set')
-#        primdf = []
-#        shfunc, func = -1, -1
-#        for seht, x, y, z in zip(universe.atom['set'], universe.atom['x'],
-#                                 universe.atom['y'], universe.atom['z']):
-#            summ = universe.basis_set_summary.ix[seht]
-#            b = bases.get_group(seht).groupby('shell_function')
-#            for sh in range(len(b)):
-#                prims = b.get_group(sh)
-#                l = prims['L'].cat.as_ordered().max()
-#                shfunc += 1
-#                for l, m, n in enum_cartesian[l]:
-#                    func += 1
-#                    for alpha, d in zip(prims['alpha'], prims['d']):
-#                        primdf.append([x, y, z, alpha, d, l, m, n, l + m + n, sh, shfunc, func, seht])
-#        primdf = pd.DataFrame(primdf)
-#        primdf.columns = ['xa', 'ya', 'za', 'alpha', 'd', 'l', 'm', 'n', 'L', 'shell_function', 'shell', 'func', 'set']
-#        if inplace:
-#            universe.gaussian_basis_set = primdf
-#        else:
-#            return cls(primdf)
-
 
 class Primitive(DataFrame):
+    """
+    Contains the required information to perform molecular integrals. Some
+    repetition of data with GaussianBasisSet but for convenience also stored
+    here.
+
+    Currently has the capability to compute the primitive overlap matrix
+    and reduce the dimensionality to the contracted cartesian overlap
+    matrix. Does not have the functionality to convert to the contracted
+    spherical overlap matrix (the fully contracted basis set of routine
+    gaussian type calculations).
+    +-------------------+----------+-------------------------------------------+
+    | Column            | Type     | Description                               |
+    +===================+==========+===========================================+
+    | xa                | float    | center in x direction of primitive        |
+    +-------------------+----------+-------------------------------------------+
+    | ya                | float    | center in y direction of primitive        |
+    +-------------------+----------+-------------------------------------------+
+    | za                | float    | center in z direction of primitive        |
+    +-------------------+----------+-------------------------------------------+
+    | alpha             | float    | value of :math:`\\alpha`, the exponent    |
+    +-------------------+----------+-------------------------------------------+
+    | d                 | float    | value of the contraction coefficient      |
+    +-------------------+----------+-------------------------------------------+
+    | l                 | int      | pre-exponential power of x                |
+    +-------------------+----------+-------------------------------------------+
+    | m                 | int      | pre-exponential power of y                |
+    +-------------------+----------+-------------------------------------------+
+    | n                 | int      | pre-exponential power of z                |
+    +-------------------+----------+-------------------------------------------+
+    | L                 | int/cat  | sum of l + m + n                          |
+    +-------------------+----------+-------------------------------------------+
+    """
     _columns = ['xa', 'ya', 'za', 'alpha', 'd', 'l', 'm', 'n', 'L']
     _indices = ['primitive']
     _categories = {'l': np.int64, 'm': np.int64, 'n': np.int64, 'L': np.int64}
@@ -277,9 +246,9 @@ def _spherical_from_cartesian(self):
         '''
         print('warning: this is not correct')
         lmax = self['L'].cat.as_ordered().max()
-        primS = self.primitive_overlap().square()
+        prim_ovl = self.primitive_overlap().square()
         cartprim, ls = self._cartesian_contraction_matrix(l=True)
-        contracted = pd.DataFrame(np.dot(np.dot(cartprim.T, primS), cartprim))
+        contracted = pd.DataFrame(np.dot(np.dot(cartprim.T, prim_ovl), cartprim))
         sh = solid_harmonics(lmax)
         sphtrans = car2sph_transform_matrices(sh, lmax)
         bfns = self.groupby('func')
@@ -308,9 +277,7 @@ def _spherical_from_cartesian(self):
 
 
     def primitive_overlap(self):
-        '''
-        Computes the complete primitive cartesian overlap matrix.
-        '''
+        """Computes the complete primitive cartesian overlap matrix."""
         if 'N' not in self.columns:
             self._normalize()
         chi1, chi2, overlap =  _wrap_overlap(self['xa'].values,
@@ -326,39 +293,16 @@ def primitive_overlap(self):
 
 
     def contracted_cartesian_overlap(self):
-        primS = self.primitive_overlap().square()
+        """Returns the contracted cartesian overlap matrix."""
+        prim_ovl = self.primitive_overlap().square()
         contprim = self._cartesian_contraction_matrix()
-        square = pd.DataFrame(np.dot(np.dot(contprim.T, primS), contprim))
+        square = pd.DataFrame(np.dot(np.dot(contprim.T, prim_ovl), contprim))
         return Overlap.from_square(square)
 
     def contracted_spherical_overlap(self):
         return self._spherical_from_cartesian()
 
 
-#    @classmethod
-#    def from_universe(cls, universe, inplace=False):
-#        bases = universe.gaussian_basis_set[abs(universe.gaussian_basis_set['d']) > 0].groupby('set')
-#        primdf = []
-#        shfunc, func = -1, -1
-#        for seht, x, y, z in zip(universe.atom['set'], universe.atom['x'],
-#                                 universe.atom['y'], universe.atom['z']):
-#            summ = universe.basis_set_summary.ix[seht]
-#            b = bases.get_group(seht).groupby('shell_function')
-#            for sh in range(len(b)):
-#                prims = b.get_group(sh)
-#                l = prims['l'].cat.as_ordered().max()
-#                shfunc += 1
-#                for l, m, n in enum_cartesian[l]:
-#                    func += 1
-#                    for alpha, d in zip(prims['alpha'].values, prims['d'].values):
-#                        primdf.append([x, y, z, alpha, d, l, m, n, l + m + n, shfunc, func])
-#        primdf = pd.DataFrame(primdf)
-#        primdf.columns = ['xa', 'ya', 'za', 'alpha', 'd', 'l', 'm', 'n', 'L', 'shell', 'func']
-#        if inplace:
-#            universe.primitive = primdf
-#        else:
-#            return cls(primdf)
-
     @classmethod
     def from_universe(cls, universe, inplace=False):
         '''
@@ -375,8 +319,8 @@ def from_universe(cls, universe, inplace=False):
                                  universe.atom['y'], universe.atom['z']):
             summ = universe.basis_set_summary.ix[seht]
             b = bases.get_group(seht).groupby('shell_function')
-            for sh in range(len(b)):
-                prims = b.get_group(sh)
+            for sh, prims in b:
+                if len(prims) == 0: continue
                 l = prims['L'].cat.as_ordered().max()
                 shfunc += 1
                 for l, m, n in enum_cartesian[l]:
@@ -391,10 +335,6 @@ def from_universe(cls, universe, inplace=False):
             return cls(primdf)
 
 
-
-
-
-
 class BasisSetOrder(BasisSet):
     """
     BasisSetOrder uniquely determines the basis function ordering scheme for
@@ -415,30 +355,6 @@ class BasisSetOrder(BasisSet):
     _index = 'chi'
     _categories = {'center': np.int64, 'symbol': str}
 
-#class BasisSetMap(BasisSet):
-#    """
-#    BasisSetMap provides the auxiliary information about relational mapping
-#    between the complete uncontracted primitive basis set and the resultant
-#    contracted basis set within an :class:`~exatomic.universe.Universe`.
-#
-#    +-------------------+----------+-------------------------------------------+
-#    | Column            | Type     | Description                               |
-#    +===================+==========+===========================================+
-#    | tag               | str      | basis set identifier                      |
-#    +-------------------+----------+-------------------------------------------+
-#    | l                 | int      | oribtal angular momentum quantum number   |
-#    +-------------------+----------+-------------------------------------------+
-#    | nprim             | int      | number of primitives within shell         |
-#    +-------------------+----------+-------------------------------------------+
-#    | nbasis            | int      | number of basis functions within shell    |
-#    +-------------------+----------+-------------------------------------------+
-#    | degen             | bool     | False if cartesian, True if spherical     |
-#    +-------------------+----------+-------------------------------------------+
-#    """
-#    _columns = ['tag', 'nprim', 'nbasis', 'degen']
-#    _indices = ['index']
-#    #_categories = {'tag': str, 'shell': str, 'nbasis': np.int64, 'degen': bool}
-#
 
 class Overlap(DataFrame):
     """
@@ -522,18 +438,6 @@ class CartesianGTFOrder(DataFrame):
     _traits = ['l']
     _categories = {'l': np.int64, 'x': np.int64, 'y': np.int64, 'z': np.int64}
 
-    def _custom_traits(self):
-        #print(self.groupby('l').apply(lambda y: y['x'].values))
-        #print(self.groupby('l').apply(lambda y: y['x'].values).to_json(orient='values'))
-        #cgto_x = self.groupby('l').apply(lambda x: x['x'].values).to_json(orient='values')
-        #cgto_x = Unicode(''.join(['[', cgto_x, ']'])).tag(sync=True)
-        #cgto_y = self.groupby('l').apply(lambda x: x['y'].values).to_json(orient='values')
-        #cgto_y = Unicode(''.join(['[', cgto_y, ']'])).tag(sync=True)
-        #cgto_z = self.groupby('l').apply(lambda x: x['z'].values).to_json(orient='values')
-        #cgto_z = Unicode(''.join(['[', cgto_z, ']'])).tag(sync=True)
-        #return {'cartesiangtforder_x': cgto_x, 'cartesiangtforder_y': cgto_y,
-        #        'cartesiangtforder_z': cgto_z}
-        return {}
 
     @classmethod
     def from_lmax_order(cls, lmax, ordering_function):
@@ -579,13 +483,6 @@ class SphericalGTFOrder(DataFrame):
     _traits = ['l']
     _index = 'spherical_order'
 
-    def _custom_traits(self):
-        sgto = self.groupby('frame').apply(lambda x: x.groupby('l').apply( lambda y: y['ml'].values))
-        sgto = Unicode(sgto.to_json(orient='values')).tag(sync=True)
-        return {'sphericalgtforder_ml': sgto}
-        #Unicode('[' + self.groupby('l').apply(
-        #        lambda x: x['ml'].values).to_json(orient='values') + ']').tag(sync=True)}
-
     @classmethod
     def from_lmax_order(cls, lmax, ordering_function):
         """
diff --git a/exatomic/editor.py b/exatomic/editor.py
@@ -35,6 +35,11 @@ def to_universe(self, name=None, description=None, meta=None):
         """
         Convert the editor to a :class:`~exatomic.container.Universe` object.
         """
+        if self.meta is not None:
+            if meta is not None:
+                meta.update(self.meta)
+            else:
+                meta = self.meta
         kwargs = {'name': name, 'description': description, 'meta': meta}
         attrs = [attr.replace('parse_', '') for attr in vars(self.__class__).keys() if attr.startswith('parse_')]
         for attr in attrs:
diff --git a/exatomic/orbital.py b/exatomic/orbital.py
@@ -160,6 +160,7 @@ class DensityMatrix(DataFrame):
     _index = 'index'
 
     def square(self, frame=0):
+        """Returns a square dataframe of the density matrix."""
         denvec = self[self['frame'] == frame]['coefficient'].values
         square = pd.DataFrame(density_as_square(denvec))
         square.index.name = 'chi1'
