enh: many features enabled for the tbtrans files, and many fixes

Changed the reciprocal cell to include 2 * pi. For newcomers this
seems like it is much more intuitive.

Updated documentation a couple of places (we need to go through all
of it and abide to:

        Parameters
        ----------
        ia : `list` of `int`
           Atomic indices
        all : `bool = False`
           `False`, return only the first orbital ...

Removed the need for opening/closing NetCDF files in the SileCDF.
This was redundant as the NetCDF file is not a linear file.
The SileCDF also implements a few of shorthands for retrieving values
vs. variables.

The sisl/io/tbtrans.py has been overhauled. Many new routines
for mingling with the orbital currents.
Now one may request direct unit-cell -> supercell bond/orbital currents
for increased flexibility.

Signed-off-by: Nick Papior <nickpapior@gmail.com>

sisl/geometry.py

@@ -401,8 +401,8 @@ def cut(self, seps, axis, seg=0, rtol=1e-4, atol=1e-4):
         It will then _only_ return the first cut.
 
         This will effectively change the unit-cell in the ``axis`` as-well
-        as removing ``self.na_u/seps`` atoms.
-        It requires that ``self.na_u % seps == 0``.
+        as removing ``self.na/seps`` atoms.
+        It requires that ``self.na % seps == 0``.
 
         REMARK: You need to ensure that all atoms within the first
         cut out region are within the primary unit-cell.
@@ -1252,41 +1252,46 @@ def a2o(self, ia, all=False):
         This is particularly handy if you want to create
         TB models with more than one orbital per atom.
 
+        Note that this will preserve the super-cell offsets.
+
         Parameters
         ----------
-        ia : list, int
+        ia : `list` of `int`
              Atomic indices
-        all: False, bool
-             ``False``, return only the first orbital corresponding to the atom,
-             ``True``, returns list of the full atom
+        all : `bool = False`
+             `False`, return only the first orbital corresponding to the atom,
+             `True`, returns list of the full atom
         """
         if not all:
             return self.lasto[ia % self.na] + (ia // self.na) * self.no
         ia = np.asarray(ia, np.int32)
         ob = self.a2o(ia)
         oe = self.a2o(ia + 1)
-        arange = np.arange
+
         # Create ranges
         if isinstance(ob, Integral):
-            return arange(ob, oe, dtype=np.int32)
+            return np.arange(ob, oe, dtype=np.int32)
 
         # Several ranges
         o = np.empty([np.sum(oe - ob)], np.int32)
         n = 0
+        narange = np.arange
         for i in range(len(ob)):
-            o[n:n + oe[i] - ob[i]] = arange(ob[i], oe[i], dtype=np.int32)
+            o[n:n + oe[i] - ob[i]] = narange(ob[i], oe[i], dtype=np.int32)
             n += oe[i] - ob[i]
         return o
 
     def o2a(self, io):
         """
         Returns an atomic index corresponding to the orbital indicies.
 
-        This is a particurlaly slow algorithm.
+        This is a particurlaly slow algorithm due to for-loops.
+
+        Note that this will preserve the super-cell offsets.
 
         Parameters
         ----------
-        io: list, int
+        io: `list` of `int`
              List of indices to return the atoms for
         """
         rlasto = self.lasto[::-1]
@@ -1456,7 +1461,7 @@ def __call__(self, parser, ns, value, option_string=None):
                     ns._geometry = ns._geometry.translate(-tmp + np.array([d,d,d]))
                 elif args.unit_cell in ['mod']:
                     # Change all coordinates using the reciprocal cell
-                    rcell = ns._geometry.rcell
+                    rcell = ns._geometry.rcell / ( 2. * np.pi )
                     idx = np.abs(np.array(np.dot(ns._geometry.xyz, rcell),np.int32))
                     # change supercell
                     nsc = np.amax(idx * 2 + 1,axis=0)

sisl/grid.py

@@ -391,11 +391,12 @@ def index(self, coord, axis=None):
 
         # if the axis is none, we do this for all axes
         if axis is None:
+            rcell = self.rcell / ( 2. * np.pi )
             # Loop over each direction
             idx = np.empty([3], np.int32)
             for i in [0, 1, 2]:
                 # get the coordinate along the direction of the cell vector
-                c = np.dot(self.rcell[i, :], coord) * self.cell[i,:]
+                c = np.dot(rcell[i, :], coord) * self.cell[i,:]
                 # Calculate the index corresponding to this placement
                 idx[i] = self.index(c, i)
             return idx

sisl/io/siesta/siesta.py

@@ -23,21 +23,19 @@ class ncSileSiesta(SileCDFSIESTA):
     """ SIESTA file object """
 
 
-    @Sile_fh_open
     def read_sc(self):
         """ Returns a SuperCell object from a SIESTA.nc file
         """
-        cell = np.array(self.variables['cell'][:], np.float64)
+        cell = np.array(self._value('cell'), np.float64)
         # Yes, this is ugly, I really should implement my unit-conversion tool
         cell *= Bohr2Ang
         cell.shape = (3, 3)
 
-        nsc = np.array(self.variables['nsc'][:], np.int32)
+        nsc = np.array(self._value('nsc'), np.int32)
 
         return SuperCell(cell, nsc=nsc)
 
 
-    @Sile_fh_open
     def read_geom(self):
         """ Returns Geometry object from a SIESTA.nc file
 
@@ -47,13 +45,13 @@ def read_geom(self):
         # Read supercell
         sc = self.read_sc()
 
-        xyz = np.array(self.variables['xa'][:], np.float64)
+        xyz = np.array(self._value('xa'), np.float64)
         xyz.shape = (-1, 3)
 
         if 'BASIS' in self.groups:
             bg = self.groups['BASIS']
             # We can actually read the exact basis-information
-            b_idx = np.array(bg.variables['basis'][:], np.int32)
+            b_idx = np.array(bg._value('basis'), np.int32)
 
             # Get number of different species
             n_b = len(bg.groups)
@@ -83,15 +81,14 @@ def read_geom(self):
         return geom
 
 
-    @Sile_fh_open
     def read_es(self, **kwargs):
         """ Returns a tight-binding model from the underlying NetCDF file """
 
         # Get the default spin channel
         ispin = kwargs.get('ispin', -1)
         spin = 1
         if ispin == -1:
-            spin = len(self.dimensions['spin'])
+            spin = len(self._dimension('spin'))
 
         # First read the geometry
         geom = self.read_geom()
@@ -108,7 +105,7 @@ def read_es(self, **kwargs):
         ham = Hamiltonian(geom, nnzpr=1, ortho=False, spin=spin)
 
         # Use Ef to move H to Ef = 0
-        Ef = float(self.variables['Ef'][0]) * Ry2eV ** ham._E_order
+        Ef = float(self._value('Ef')[0]) * Ry2eV ** ham._E_order
         S = np.array(sp.variables['S'][:], np.float64)
 
         ncol = np.array(sp.variables['n_col'][:], np.int32)
@@ -149,7 +146,6 @@ def read_es(self, **kwargs):
         return ham
 
 
-    @Sile_fh_open
     def grids(self):
         """ Return a list of available grids in this file. """
 
@@ -159,7 +155,6 @@ def grids(self):
 
         return grids
 
-    @Sile_fh_open
     def read_grid(self, name, idx=0):
         """ Reads a grid in the current SIESTA.nc file
 
@@ -203,7 +198,6 @@ def read_grid(self, name, idx=0):
 
         return grid
 
-    @Sile_fh_open
     def write_geom(self, geom):
         """
         Creates the NetCDF file and writes the geometry information
@@ -271,7 +265,6 @@ def write_geom(self, geom):
         self.variables['lasto'][:] = np.cumsum(orbs)
 
 
-    @Sile_fh_open
     def write_es(self, ham, **kwargs):
         """ Writes Hamiltonian model to file
 

sisl/io/siesta/siesta_grid.py

@@ -21,18 +21,16 @@
 class gridncSileSiesta(SileCDFSIESTA):
     """ SIESTA Grid file object """
 
-    @Sile_fh_open
     def read_sc(self):
         """ Returns a SuperCell object from a SIESTA.grid.nc file
         """
-        cell = np.array(self.variables['cell'][:], np.float64)
+        cell = np.array(self._value('cell'), np.float64)
         # Yes, this is ugly, I really should implement my unit-conversion tool
         cell *= Bohr2Ang
         cell.shape = (3, 3)
 
         return SuperCell(cell)
 
-    @Sile_fh_open
     def read_grid(self, name='gridfunc', idx=0, *args, **kwargs):
         """ Reads a grid in the current SIESTA.grid.nc file
 
@@ -42,14 +40,14 @@ def read_grid(self, name='gridfunc', idx=0, *args, **kwargs):
         sc = self.read_sc().swapaxes(0, 2)
 
         # Create the grid
-        nx = len(self.dimensions['n1'])
-        ny = len(self.dimensions['n2'])
-        nz = len(self.dimensions['n3'])
+        nx = len(self._dimension('n1'))
+        ny = len(self._dimension('n2'))
+        nz = len(self._dimension('n3'))
 
         if name is None:
-            v = self.variables['gridfunc']
+            v = self._variable('gridfunc')
         else:
-            v = self.variables[name]
+            v = self._variable(name)
 
         # Create the grid, SIESTA uses periodic, always
         grid = Grid([nz, ny, nx], bc=Grid.Periodic, sc=sc,