Merge pull request #118 from tjduigna/master

Update to orbital viz., tidied up gaussian parsers

Author: avmarchenko

docs/source/notebooks/04_orbitals.ipynb

@@ -297,13 +297,12 @@ def evaluate(self, xs=None, ys=None, zs=None, irrep=None, verbose=False):
         if self._meta['gaussian']:
             if self._meta.get('symmetrized', False):
                 func = self._evaluate_gau_bso_sym
-            elif self._meta['program'] in ['nwchem']:
-                func = self._evaluate_gau_bso
             elif self._meta['program'] in ['molcas']:
                 func = self._evaluate_gau_mag
+            else:
+                func = self._evaluate_gau_bso
         else:
             func = self._evaluate_sto
-        # if verbose: print('calling {}'.format(func.__name__))
         return func(xs=xs, ys=ys, zs=zs, irrep=irrep)
 
 

exatomic/algorithms/basis.py

@@ -8,6 +8,7 @@
 import os
 from exa.util import isotopes
 from platform import system
+from IPython.display import display_html
 
 # For numba compiled functions
 sysname= system().lower()
@@ -69,3 +70,10 @@ def list_resources():
     for path, _, files_ in os.walk(staticdir()):
         files.extend(files_)
     return files
+
+
+def display_side_by_side(*args):
+    """Simple function to display 2 dataframes side by side in a notebook."""
+    html_str = ''.join([df.to_html() for df in args])
+    display_html(html_str.replace('table','table style=\"display:inline\"'),
+                 raw=True)

exatomic/base.py

@@ -162,6 +162,59 @@ def primitives(self, spherical):
         #self.gaussian = gaussian
 
 
+def deduplicate_basis_sets(sets, sp=False):
+    """Deduplicate identical basis sets on different centers.
+
+    Args:
+        sets (pd.DataFrame): non-unique basis sets
+        sp (bool): Whether or not to call _expand_sp (gaussian program only)
+
+    Returns:
+        tup (tuple): deduplicated basis sets and basis set map for atom table
+    """
+    unique, setmap, cnt = [], {}, 0
+    sets = sets.groupby('center')
+    chk = ['alpha', 'd']
+    for center, seht in sets:
+        for i, other in enumerate(unique):
+            if other.shape != seht.shape: continue
+            if np.allclose(other[chk], seht[chk]):
+                setmap[center] = i
+                break
+        else:
+            unique.append(seht)
+            setmap[center] = cnt
+            cnt += 1
+    if sp: unique = _expand_sp(unique)
+    sets = pd.concat(unique).reset_index(drop=True)
+    try: sets.drop([2, 3], axis=1, inplace=True)
+    except (KeyError, ValueError): pass
+    sets.rename(columns={'center': 'set'}, inplace=True)
+    sets['set'] = sets['set'].map(setmap)
+    sets['frame'] = 0
+    return sets, setmap
+
+def _expand_sp(unique):
+    """Currently only used when 'program' == 'gaussian'."""
+    expand = []
+    for seht in unique:
+        if np.isnan(seht[2]).sum() == seht.shape[0]:
+            expand.append(seht)
+            continue
+        sps = seht[2][~np.isnan(seht[2])].index
+        shls = len(seht.ix[sps]['shell'].unique())
+        dupl = seht.ix[sps[0]:sps[-1]].copy()
+        dupl[1] = dupl[2]
+        dupl['L'] = 1
+        dupl['shell'] += shls
+        last = seht.ix[sps[-1] + 1:].copy()
+        last['shell'] += shls
+        expand.append(pd.concat([seht.ix[:sps[0] - 1],
+                                 seht.ix[sps[0]:sps[-1]],
+                                 dupl, last]))
+    return expand
+
+
 class BasisSetOrder(DataFrame):
     """
     BasisSetOrder uniquely determines the basis function ordering scheme for
@@ -223,22 +276,38 @@ class Overlap(DataFrame):
     _columns = ['chi0', 'chi1', 'coef', 'frame']
     _index = 'index'
 
-    #@property
-    #def _constructor(self):
-    #    return Overlap
 
-    def square(self, frame=0, column='coef', mocoefs=None):
+    def square(self, frame=0, column='coef', mocoefs=None, irrep=None):
         """Return a 'square' matrix DataFrame of the Overlap.
 
         Args:
             column (str): column of coefficients to reshape
             mocoefs (str): alias for `column`
             frame (int): default 0
+            irrep (int): irreducible representation if symmetrized
         """
         if mocoefs is not None: column = mocoefs
+        if 'irrep' in self.columns:
+            if irrep is None:
+                irreps, i, j = self.groupby('irrep'), 0, 0
+                norb = (irreps.chi0.max() + 1).sum()
+                nchi = (irreps.chi1.max() + 1).sum()
+                cmat = np.zeros((nchi, norb))
+                for irrep, grp in irreps:
+                    piv = grp.pivot('chi0', 'chi1', column)
+                    ii, jj = piv.shape
+                    cmat[i : i + ii, j : j + jj] = piv.values
+                    i += ii
+                    j += jj
+                idx = pd.Index(range(nchi), name='chi0')
+                orb = pd.Index(range(norb), name='chi1')
+                return pd.DataFrame(cmat, index=idx, columns=orb)
+            return self.groupby('irrep').get_group(irrep
+                        ).pivot('chi', 'orbital', column)
         sq = _square(self[column].values)
-        return pd.DataFrame(sq, index=pd.Index(range(sq.shape[0]), name='chi0'),
-                            columns=pd.Index(range(sq.shape[1]), name='chi1'))
+        idx = pd.Index(range(sq.shape[0]), name='chi0')
+        orb = pd.Index(range(sq.shape[1]), name='chi1')
+        return pd.DataFrame(sq, index=idx, columns=orb)
 
     @classmethod
     def from_column(cls, source):

exatomic/core/basis.py

@@ -204,7 +204,7 @@ def enumerate_shells(self, frame=0):
             frame (int): state of the universe (default 0)
         """
         atom = self.atom.groupby('frame').get_group(frame)
-        if self.meta['program'] not in ['molcas', 'adf', 'nwchem']:
+        if self.meta['program'] not in ['molcas', 'adf', 'nwchem', 'gaussian']:
             print('Warning: Check spherical shell parameter for {} '
                   'molecular orbital generation'.format(self.meta['program']))
         shls = self.basis_set.shells(self.meta['program'],