Python scipy.sparse.bmat() Examples

def rotate_mag(m,angle=0.0,axis=[0,0,1]):
  """ Rotates the magnetization along a particular axis """
  nat = len(m)/2

  #####################
  ## rotate the dmat ##
  #####################
  from scipy.sparse import coo_matrix, bmat
  sx = np.matrix([[0.,1.],[1.,0.]])
  sy = np.matrix([[0.,1j],[-1j,0.]])
  sz = np.matrix([[1.,0.],[0.,-1.]])

  from scipy.linalg import expm
  rmat = expm(1j*np.pi*angle*(axis[0]*sx + axis[1]*sy + axis[2]*sz))
  rmat = coo_matrix(rmat) # to sparse
  rot = [[None for i in range(nat)] for j in range(nat)]
  for i in range(nat):
    rot[i][i] = rmat
  rot = bmat(rot).todense() # rotational matrix

  # rotate the matrix
  m = rot * m * rot.H

  return m 

def _compute_jacobian(self, J_eq, J_ineq, s):
        if self.n_ineq == 0:
            return J_eq
        else:
            if sps.issparse(J_eq) or sps.issparse(J_ineq):
                # It is expected that J_eq and J_ineq
                # are already `csr_matrix` because of
                # the way ``BoxConstraint``, ``NonlinearConstraint``
                # and ``LinearConstraint`` are defined.
                J_eq = sps.csr_matrix(J_eq)
                J_ineq = sps.csr_matrix(J_ineq)
                return self._assemble_sparse_jacobian(J_eq, J_ineq, s)
            else:
                S = np.diag(s)
                zeros = np.zeros((self.n_eq, self.n_ineq))
                # Convert to matrix
                if sps.issparse(J_ineq):
                    J_ineq = J_ineq.toarray()
                if sps.issparse(J_eq):
                    J_eq = J_eq.toarray()
                # Concatenate matrices
                return np.asarray(np.bmat([[J_eq, zeros],
                                           [J_ineq, S]])) 

def _compute_jacobian(self, J_eq, J_ineq, s):
        if self.n_ineq == 0:
            return J_eq
        else:
            if spc.issparse(J_eq) or spc.issparse(J_ineq):
                # It is expected that J_eq and J_ineq
                # are already `csr_matrix` because of
                # the way ``BoxConstraint``, ``NonlinearConstraint``
                # and ``LinearConstraint`` are defined.
                J_eq = spc.csr_matrix(J_eq)
                J_ineq = spc.csr_matrix(J_ineq)
                return self._assemble_sparse_jacobian(J_eq, J_ineq, s)
            else:
                S = np.diag(s)
                zeros = np.zeros((self.n_eq, self.n_ineq))
                # Convert to matrix
                if spc.issparse(J_ineq):
                    J_ineq = J_ineq.toarray()
                if spc.issparse(J_eq):
                    J_eq = J_eq.toarray()
                # Concatenate matrices
                return np.asarray(np.bmat([[J_eq, zeros],
                                           [J_ineq, S]])) 

def _merge_blocks(blocks):
        """
        Helper function that merges the _blocks attribute of a ds-array into
        a single ndarray / sparse matrix.
        """
        sparse = None
        b0 = blocks[0][0]
        if sparse is None:
            sparse = issparse(b0)

        if sparse:
            ret = sp.bmat(blocks, format=b0.getformat(), dtype=b0.dtype)
        else:
            ret = np.block(blocks)

        return ret 

def edge_dos(intra0,inter0,scale=4.,w=20,npol=300,ne=500,bulk=False,
                use_random=True,nrand=20):
  """Calculated the edge DOS using the KPM"""
  h = [[None for j in range(w)] for i in range(w)]
  intra = csc_matrix(intra0)
  inter = csc_matrix(inter0)
  for i in range(w): h[i][i] = intra
  for i in range(w-1): 
    h[i+1][i] = inter.H
    h[i][i+1] = inter
  h = bmat(h) # sparse hamiltonian
  ds = np.zeros(ne)
  dsb = np.zeros(ne)
  norb = intra0.shape[0] # orbitals ina cell
  for i in range(norb):
    (xs,ys) = ldos0d(h,i=i,scale=scale,npol=npol,ne=ne) 
    ds += ys # store
    if bulk:
      (xs,zs) = ldos0d(h,i=w*norb//2 + i,scale=scale,npol=npol,ne=ne) 
      dsb += zs # store
  if not bulk: return (xs,ds/w)
  else: return (xs,ds/w,dsb/w) 

def add_interlayer(t,ri,rj,has_spin=True,is_sparse=True):
  """Calculate interlayer coupling"""
  m = np.matrix([[0. for i in ri] for j in rj])
  if has_spin: m = bmat([[csc(m),None],[None,csc(m)]]).todense()
  zi = [r[2] for r in ri]
  zj = [r[2] for r in rj]
  for i in range(len(ri)): # add the interlayer
    for j in range(len(ri)): # add the interlayer
      rij = ri[i] - rj[j] # distance between atoms
      dz = zi[i] - zj[j] # vertical distance
      if (3.99<rij.dot(rij)<4.01) and (3.99<(dz*dz)<4.01): # check if connect
        if has_spin: # spin polarized
          m[2*i,2*j] = t
          m[2*i+1,2*j+1] = t
        else:  # spin unpolarized
          m[i,j] = t
  return m 

def bulk2ribbon_zz(h,n=10):
  """Converts a hexagonal bulk hamiltonian into a ribbon hamiltonian"""
  h = honeycomb2square(h) # create a square 2d geometry
  go = h.geometry.copy() # copy geometry
  ho = h.copy() # copy hamiltonian
  ho.dimensionality = 1 # reduce dimensionality
  go.dimensionality = 1 # reduce dimensionality
  intra = [[None for i in range(n)] for j in range(n)]
  inter = [[None for i in range(n)] for j in range(n)]
  for i in range(n): # to the the sam index
    intra[i][i] = csc(h.intra) 
    inter[i][i] = csc(h.tx) 
  for i in range(n-1): # one more or less
    intra[i][i+1] = csc(h.ty)  
    intra[i+1][i] = csc(h.ty.H)  
    inter[i+1][i] = csc(h.txmy) 
    inter[i][i+1] = csc(h.txy) 
  ho.intra = bmat(intra).todense()
  ho.inter = bmat(inter).todense()
  return ho 

def ldos_finite(h,e=0.0,n=10,nwf=4,delta=0.0001):
  """Calculate the density of states for a finite system"""
  if h.dimensionality!=1: raise # if it is not one dimensional
  intra = csc(h.intra) # convert to sparse
  inter = csc(h.inter) # convert to sparse
  interH = inter.H # hermitian
  m = [[None for i in range(n)] for j in range(n)] # full matrix
  for i in range(n): # add intracell
    m[i][i] = intra
  for i in range(n-1): # add intercell
    m[i][i+1] = inter
    m[i+1][i] = interH
  m = bmat(m) # convert to matrix
  (ene,wfs) = slg.eigsh(m,k=nwf,which="LM",sigma=0.0) # diagonalize
  wfs = wfs.transpose() # transpose wavefunctions
  dos = (wfs[0].real)*0.0 # calculate dos
  for (ie,f) in zip(ene,wfs): # loop over waves
    c = 1./(1.+((ie-e)/delta)**2) # calculate coefficient
    dos += np.abs(f)*c # add contribution
  odos = spatial_dos(h,dos) # get the spatial distribution
  go = h.geometry.supercell(n) # get the supercell
  write_ldos(go.x,go.y,odos) # write in a file
  return dos # return the dos 

def block_inverse(m,i=0,j=0):
  """ Calculate a certain element of the inverse of a block matrix"""
  from scipy.sparse import csc_matrix,bmat
  nb = len(m) # number of blocks
  if i<0: i += nb 
  if j<0: j += nb 
  mt = [[None for ii in range(nb)] for jj in range(nb)]
  for ii in range(nb): # diagonal part
    mt[ii][ii] = csc_matrix(m[ii][ii])
  for ii in range(nb-1):
    mt[ii][ii+1] = csc_matrix(m[ii][ii+1])
    mt[ii+1][ii] = csc_matrix(m[ii+1][ii])
  mt = bmat(mt).todense() # create dense matrix
  # select which elements you need
  ilist = [m[ii][ii].shape[0] for ii in range(i)] 
  jlist = [m[jj][jj].shape[1] for jj in range(j)] 
  imin = sum(ilist)
  jmin = sum(jlist)
  mt = mt.I # calculate inverse
  imax = imin + m[i][i].shape[0]
  jmax = jmin + m[j][j].shape[1]
  mo = [ [mt[ii,jj] for jj in range(jmin,jmax)] for ii in range(imin,imax) ] 
  mo = np.matrix(mo)
  return mo 

def interface(h1,h2,k=[0.0,0.,0.],energy=0.0,delta=0.01):
  """Get the Green function of an interface"""
  from scipy.sparse import csc_matrix as csc
  from scipy.sparse import bmat
  gs1,sf1 = green_kchain(h1,k=k,energy=energy,delta=delta,
                   only_bulk=False,reverse=True) # surface green function 
  gs2,sf2 = green_kchain(h2,k=k,energy=energy,delta=delta,
                   only_bulk=False,reverse=False) # surface green function 
  #############
  ## 1  C  2 ##
  #############
  # Now apply the Dyson equation
  (ons1,hop1) = get1dhamiltonian(h1,k,reverse=True) # get 1D Hamiltonian
  (ons2,hop2) = get1dhamiltonian(h2,k,reverse=False) # get 1D Hamiltonian
  havg = (hop1.H + hop2)/2. # average hopping
  ons = bmat([[csc(ons1),csc(havg)],[csc(havg.H),csc(ons2)]]) # onsite
  self2 = bmat([[csc(ons1)*0.0,None],[None,csc(hop2@sf2@hop2.H)]])
  self1 = bmat([[csc(hop1@sf1@hop1.H),None],[None,csc(ons2)*0.0]])
  # Dyson equation
  ez = (energy+1j*delta)*np.identity(ons1.shape[0]+ons2.shape[0]) # energy
  ginter = (ez - ons - self1 - self2).I # Green function
  # now return everything, first, second and hybrid
  return (gs1,sf1,gs2,sf2,ginter) 

def bulk_and_surface(h1,nk=100,energies=np.linspace(-1.,1.,100),**kwargs):
  """Get the surface DOS of an interface"""
  from scipy.sparse import csc_matrix,bmat
  if h1.dimensionality==2:
      kpath = [[k,0.,0.] for k in np.linspace(0.,1.,nk)]
  else: raise
  ik = 0
  h1 = h1.get_multicell() # multicell Hamiltonian
  tr = timing.Testimator("DOS") # generate object
  dos_bulk = energies*0.0
  dos_sur = energies*0.0
  for k in kpath:
    tr.remaining(ik,len(kpath)) # generate object
    ik += 1
    outs = green.surface_multienergy(h1,k=k,energies=energies,**kwargs)
    dos_bulk += np.array([-algebra.trace(g[1]).imag for g in outs])
    dos_sur += np.array([-algebra.trace(g[0]).imag for g in outs])
  dos_bulk /= len(kpath)
  dos_sur /= len(kpath)
  np.savetxt("DOS.OUT",np.array([energies,dos_bulk,dos_sur]).T)
  return energies,dos_bulk,dos_sur 

def get_interface(h,fun=None):
  """Return an operator that projects onte the interface"""
  dind = 1 # index to which divide the positions
  if h.has_spin:  dind *= 2 # duplicate for spin
  if h.has_eh:  dind *= 2  # duplicate for eh
  iden = csc(np.matrix(np.identity(dind,dtype=np.complex))) # identity matrix
  r = h.geometry.r # positions
  out = [[None for ri in r] for rj in r] # initialize
  if fun is None: # no input function
    cut = 2.0 # cutoff
    if h.dimensionality==1: index = 1
    elif h.dimensionality==2: index = 2
    else: raise
    def fun(ri): # define the function
      if np.abs(ri[index])<cut: return 1.0
      else: return 0.0
  for i in range(len(r)): # loop over positions
    out[i][i] = fun(r[i])*iden 
  return bmat(out) # return matrix 

def get_si(h,i=1):
  """Return a certain Pauli matrix for the full Hamiltonian"""
  if not h.has_spin: return None # no spin
  if i==1: si = sx # sx matrix
  elif i==2: si = sy # sy matrix
  elif i==3: si = sz # sz matrix
  else: raise # unknown pauli matrix
  if h.has_eh: ndim = h.intra.shape[0]//4 # half the dimension
  else: ndim = h.intra.shape[0]//2 # dimension
  if h.has_spin: # spinful system
    op = [[None for i in range(ndim)] for j in range(ndim)] # initialize
    for i in range(ndim): op[i][i] = si # store matrix
    op = bmat(op) # create matrix
  if h.has_eh: op = build_eh(op,is_sparse=True) # add electron and hole parts 
  return op

# define the functions for the three spin components 

def eh_operator(m):
  """Return the electron hole symmetry operator, as a function"""
  n = m.shape[0]//4 # number of sites 
  from .hamiltonians import sy
  msy = [[None for ri in range(n)] for j in range(n)]
  for i in range(n): msy[i][i] = sy # add sy
  msy = bmat(msy) # sy matrix in the electron subspace
  out = [[None,1j*msy],[-1j*msy,None]]
#  out = [[None for i in range(n)] for j in range(n)] # output matrix
#  tau = csc_matrix([[0.,0.,0.,-1.],[0,0,1,0],[0,1,0,0],[-1,0,0,0]])
#  for i in range(n):
#    out[i][i] = tau
  out = bmat(out) # sparse matrix
  out = reorder(out) # reshufle the matrix
  def ehop(inm):
    """Function that applies electron-hole symmetry to a matrix"""
    return out*np.conjugate(inm)*out # return matrix
  return ehop # return operator 

def floquet_selfenergy(selfgen,e,omega,n=20,less=False,delta=0.01):
  """Generate a floquet selfenergy starting from
  a function that returns selfenergies"""
  out = [[None for i in range(2*n+1)] for i in range(2*n+1)] # output
  for i in range(-n,n+1): # loop
    ebar = e+i*omega # floquet energy
    term = selfgen(ebar) # call and store
    if less: # special propagator with the < symbol
#      bfac = 1./(np.exp(ebar/delta)+1) # boltzman factor
#      term *= bfac
      if ebar<0.: term *= 0.
      else: term = -(term-term.H)
    out[i+n][i+n] = term # store
  return bmat(out) # return dense matrix 

def spiralhopping(m,ri,rj,svector = np.array([0.,0.,1.]),
        qvector=[1.,0.,0.]): 
  """ Rotates a hopping matrix to create a spin spiral
  antsaz
      - ri and rj must be coordinates in lattice constants
      - svector is the axis of the rotation
      - qvector is the vector of the spin spiral
  """
  from scipy.sparse import csc_matrix,bmat
  iden = csc_matrix([[1.,0.],[0.,1.]]) # identity matrix
  def getR(r):
      """Return a rotation matrix"""
      n = len(r) # number of sites 
      R = [[None for i in range(n)] for j in range(n)] # rotation matrix
      u = np.array(svector) # rotation direction
      u = u/np.sqrt(u.dot(u)) # normalize rotation direction
      for i in range(n): # loop over sites
         rot = u[0]*sx + u[1]*sy + u[2]*sz 
         angle = np.array(qvector).dot(np.array(r[i])) # angle of rotation
         # a factor 2 is taken out due to 1/2 of S
         # a factor 2 is added to have BZ in the interval 0,1
         R[i][i] = lg.expm(2.*np.pi*1j*rot*angle/2.0)
      return bmat(R)  # convert to full sparse matrix
  Roti = getR(ri) # get the first rotation matrix
  Rotj = getR(rj) # get the second rotation matrix
#  print(Roti@Rotj.H)
#  print(ri,rj)
  return Rotj @ m @ Roti.H # return the rotated matrix 

def pseudo_hashimoto(graph):
    """Return the pseudo-Hashimoto matrix.

    The pseudo Hashimoto matrix of a graph is the block matrix defined as
    B' = [0  D-I]
         [-I  A ]

    Where D is the degree-diagonal matrix, I is the identity matrix and A
    is the adjacency matrix.  The eigenvalues of B' are always eigenvalues
    of B, the non-backtracking or Hashimoto matrix.

    Parameters
    ----------

    graph (nx.Graph): A NetworkX graph object.

    Returns
    -------

    A sparse matrix in csr format.

    """
    # Note: the rows of nx.adjacency_matrix(graph) are in the same order as
    # the list returned by graph.nodes().
    degrees = graph.degree()
    degrees = sparse.diags([degrees[n] for n in graph.nodes()])
    adj = nx.adjacency_matrix(graph)
    ident = sparse.eye(graph.order())
    pseudo = sparse.bmat([[None, degrees - ident], [-ident, adj]])
    return pseudo.asformat('csr') 

def intra_super2d(h,n=1,central=None):
  """ Calculates the intra of a 2d system"""
  from scipy.sparse import bmat
  from scipy.sparse import csc_matrix as csc
  tx = csc(h.tx)
  ty = csc(h.ty)
  txy = csc(h.txy)
  txmy = csc(h.txmy)
  intra = csc(h.intra)
  for i in range(n):
    intrasuper[i][i] = intra # intracell
    (x1,y1) = inds[i]
    for j in range(n):
      (x2,y2) = inds[j]
      dx = x2-x1
      dy = y2-y1
      if dx==1 and  dy==0: intrasuper[i][j] = tx
      if dx==-1 and dy==0: intrasuper[i][j] = tx.H
      if dx==0 and  dy==1: intrasuper[i][j] = ty
      if dx==0 and  dy==-1: intrasuper[i][j] = ty.H
      if dx==1 and  dy==1: intrasuper[i][j] = txy
      if dx==-1 and dy==-1: intrasuper[i][j] = txy.H
      if dx==1 and  dy==-1: intrasuper[i][j] = txmy
      if dx==-1 and dy==1: intrasuper[i][j] = txmy.H
  # substitute the central cell if it is the case
  if central!=None:
    ic = (n-1)/2
    intrasuper[ic][ic] = central
  # now convert to matrix
  intrasuper = bmat(intrasuper).todense() # supercell
  return intrasuper 

def rotation_matrix(m,vectors):
  """ Rotates a matrix, to align its components with the direction
  of the magnetism """
  if not len(m)==2*len(vectors): # stop if they don't have
                                  # compatible dimensions
     raise
  # pauli matrices
  n = len(m)/2 # number of sites
  R = [[None for i in range(n)] for j in range(n)] # rotation matrix
  from scipy.linalg import expm  # exponenciate matrix
  for (i,v) in zip(range(n),vectors): # loop over sites
    vv = np.sqrt(v.dot(v)) # norm of v
    if vv>0.000001: # if nonzero scale
      u = v/vv
    else: # if zero put to zero
      u = np.array([0.,0.,0.,])
#    rot = u[0]*sx + u[1]*sy + u[2]*sz 
    uxy = np.sqrt(u[0]**2 + u[1]**2) # component in xy plane
    phi = np.arctan2(u[1],u[0])
    theta = np.arctan2(uxy,u[2])
    r1 =  phi*sz/2.0 # rotate along z
    r2 =  theta*sy/2.0 # rotate along y
    # a factor 2 is taken out due to 1/2 of S
    rot = expm(1j*r2) * expm(1j*r1)   
    R[i][i] = rot  # save term
  R = bmat(R)  # convert to full sparse matrix
  return R.todense() 

def global_spin_rotation(m,vector = np.array([0.,0.,1.]),angle = 0.0,
                             spiral = False,atoms = None):
  """ Rotates a matrix along a certain qvector """
  # pauli matrices
  from scipy.sparse import csc_matrix,bmat
  iden = csc_matrix([[1.,0.],[0.,1.]])
  n = m.shape[0]//2 # number of sites
  if atoms==None: 
    atoms = range(n) # all the atoms
  else: 
    raise
  R = [[None for i in range(n)] for j in range(n)] # rotation matrix
  for i in range(n): # loop over sites
    u = np.array(vector) # rotation direction
    u = u/np.sqrt(u.dot(u)) # normalize rotation direction
    rot = (u[0]*sx + u[1]*sy + u[2]*sz)/2. # rotation
    # a factor 2 is taken out due to 1/2 of S
    # a factor 2 is added to have BZ in the interval 0,1
    rot = lg.expm(2.*np.pi*1j*rot*angle/2.0)
#    if i in atoms:
    R[i][i] = rot  # save term
#    else:
#      R[i][i] = iden  # save term
  R = bmat(R)  # convert to full sparse matrix
  if spiral:  # for spin spiral
    mout = R @ m  # rotate matrix
  else:  # normal global rotation
    mout = R @ m @ R.H  # rotate matrix
  return mout # return dense matrix 

def floquet_tauz(m,n=20):
  """Calculate the floquet Hamiltonian"""
  # hamgen is a function that return a time dependent hamiltonian
  out = [[None for i in range(2*n+1)] for i in range(2*n+1)] # output
  for i in range(-n,n+1): # loop
    out[i+n][i+n] = m # store
  return bmat(out) # return dense matrix 

def element(i,n,p,d=2,j=None):
  if j is None: j=i
  o = csc_matrix(([1.0],[[d*i+p[0]],[d*j+p[1]]]),
                   shape=(n*d,n*d),dtype=np.complex)  
  return o
  o = csc_matrix(([1.0],[[p[0]],[p[1]]]),shape=(d,d),dtype=np.complex)
  m = [[None for i1 in range(n)] for i2 in range(n)]
  for i1 in range(n): m[i1][i1] = zero(d)
  m[i][j] = o.copy()
  return bmat(m) # return matrix 

def _assemble_sparse_jacobian(self, J_eq, J_ineq, s):
        """Assemble sparse jacobian given its components.

        Given ``J_eq``, ``J_ineq`` and ``s`` returns:
            jacobian = [ J_eq,     0     ]
                       [ J_ineq, diag(s) ]

        It is equivalent to:
            sps.bmat([[ J_eq,   None    ],
                      [ J_ineq, diag(s) ]], "csr")
        but significantly more efficient for this
        given structure.
        """
        n_vars, n_ineq, n_eq = self.n_vars, self.n_ineq, self.n_eq
        J_aux = sps.vstack([J_eq, J_ineq], "csr")
        indptr, indices, data = J_aux.indptr, J_aux.indices, J_aux.data
        new_indptr = indptr + np.hstack((np.zeros(n_eq, dtype=int),
                                         np.arange(n_ineq+1, dtype=int)))
        size = indices.size+n_ineq
        new_indices = np.empty(size)
        new_data = np.empty(size)
        mask = np.full(size, False, bool)
        mask[new_indptr[-n_ineq:]-1] = True
        new_indices[mask] = n_vars+np.arange(n_ineq)
        new_indices[~mask] = indices
        new_data[mask] = s
        new_data[~mask] = data
        J = sps.csr_matrix((new_data, new_indices, new_indptr),
                           (n_eq + n_ineq, n_vars + n_ineq))
        return J 

def generate_soc(specie,value,input_file="wannier.win",nsuper=1,path=None):
  """Add SOC to a hamiltonian based on wannier.win"""
  if path is not None: 
      inipath = os.getcwd() # current path
      os.chdir(path) # go there
  o = open(".soc.status","w")
  iat = 1 # atom counter
  orbnames = names_soc_orbitals(specie) # get which are the orbitals
  ls = soc_l((len(orbnames)-1)/2) # get the SOC matrix
  norb = get_num_wannier(input_file) # number of wannier orbitals
  m = np.matrix([[0.0j for i in range(norb*2)] for j in range(norb*2)]) # matrix
  nat = get_num_atoms(specie,input_file) # number of atoms of this specie
  for iat in range(nat):
   # try:
#      fo.write("Attempting "+specie+"  "+str(iat+1)+"\n")
      for i in range(len(orbnames)): # loop over one index
        orbi = orbnames[i] 
        for j in range(len(orbnames)): # loop over other index
          orbj = orbnames[j]
          ii = get_index_orbital(specie,iat+1,orbi)  # index in wannier
          jj = get_index_orbital(specie,iat+1,orbj) # index in wannier
#          fo.write(str(ii)+"   "+str(jj)+"\n")
          ii += -1 # python starts in 0
          jj += -1 # python starts in 0
          m[2*ii,2*jj] = ls[2*i,2*j] # store the soc coupling
          m[2*ii+1,2*jj] = ls[2*i+1,2*j] # store the soc coupling
          m[2*ii,2*jj+1] = ls[2*i,2*j+1] # store the soc coupling
          m[2*ii+1,2*jj+1] = ls[2*i+1,2*j+1] # store the soc coupling
   #   return
   # except: break
#  fo.close()
  n = nsuper**2 # supercell
  mo = [[None for i in range(n)] for j in range(n)]
  for i in range(n): mo[i][i] = csc(m) # diagonal elements
  mo = bmat(mo).todense() # dense matrix
  if path is not None: 
      os.chdir(inipath) # go there
      print(np.max(np.abs(mo)))
  return np.matrix(mo*value) # return matrix
#  for name in atoms: # loop over atoms 

def surface(h1,energies=np.linspace(-1.,1.,100),operator=None,
                    delta=0.01,kpath=None,hs=None):
  """Get the surface DOS of an interface"""
  from scipy.sparse import csc_matrix,bmat
  if kpath is None: 
    if h1.dimensionality==3:
      g2d = h1.geometry.copy() # copy Hamiltonian
      g2d = sculpt.set_xy_plane(g2d)
      kpath = klist.default(g2d,nk=100)
    elif h1.dimensionality==2:
      kpath = [[k,0.,0.] for k in np.linspace(0.,1.,40)]
    else: raise
  fo = open("KDOS.OUT","w")
  fo.write("# k, E, Surface, Bulk\n")
  tr = timing.Testimator("KDOS") # generate object
  ik = 0
  h1 = h1.get_multicell() # multicell Hamiltonian
  for k in kpath:
    tr.remaining(ik,len(kpath)) # generate object
    ik += 1
    outs = green.surface_multienergy(h1,k=k,energies=energies,delta=delta,hs=hs)
    for (energy,out) in zip(energies,outs):
      # write everything
      fo.write(str(ik)+"   "+str(energy)+"   ")
      for g in out: # loop
        if operator is None: d = -g.trace()[0,0].imag # only the trace 
        elif callable(operator): d = operator(g,k=k) # call the operator
        else:  d = -(g@operator).trace()[0,0].imag # assume it is a matrix
        fo.write(str(d)+"   ") # write in a file
      fo.write("\n") # next line
      fo.flush() # flush
  fo.close() 

def densebmat(ms):
    """Turn a block matrix dense"""
    ms = [[todense(mi) for mi in mij] for mij in m]
    return todense(bmat(ms)) # return block matrix 

def parametric_hopping_spinful(r1,r2,fc,is_sparse=False):
  """ Generates a parametric hopping based on a function, that returns
  a 2x2 matrix"""
  m = [[None for i in range(len(r2))] for j in range(len(r1))]
  for i in range(len(r1)):
    for j in range(len(r2)):
      val = fc(r1[i],r2[j]) # add hopping based on function
      m[i][j] = val # store this result
  m = bmat(m) # convert to matrix
  if not is_sparse: m = m.todense() # dense matrix
  return m 

def get_pairing(h,ptype="s"):
  """Return an operator that calculates the expectation value of the
  s-wave pairing"""
  if not h.has_eh: raise # only for e-h systems
  if ptype=="s": op = superconductivity.spair
  elif ptype=="deltax": op = superconductivity.deltax
  elif ptype=="deltay": op = superconductivity.deltay
  elif ptype=="deltaz": op = superconductivity.deltaz
  else: raise
  r = h.geometry.r
  out = [[None for ri in r] for rj in r]
  for i in range(len(r)): # loop over positions
    out[i][i] = op
  return bmat(out) # return matrix 

def get_electron(h):
  """Operator to project on the electron sector"""
  if not h.has_eh:
      return np.identity(h.intra.shape[0])
  elif h.check_mode("spinful_nambu"): # only for e-h systems
      op = superconductivity.proje
      r = h.geometry.r
      out = [[None for ri in r] for rj in r]
      for i in range(len(r)): # loop over positions
        out[i][i] = op
      return bmat(out)
  elif h.check_mode("spinless_nambu"):
      from .sctk import spinless
      return spinless.proje(h.intra.shape[0])
  else: raise 

def pseudo_hashimoto(graph):
    """
    Return the pseudo-Hashimoto matrix.

    The pseudo Hashimoto matrix of a graph is the block matrix defined as
    B' = [0  D-I]
         [-I  A ]

    Where D is the degree-diagonal matrix, I is the identity matrix and A
    is the adjacency matrix.  The eigenvalues of B' are always eigenvalues
    of B, the non-backtracking or Hashimoto matrix.

    Parameters
    ----------

    graph (nx.Graph): A NetworkX graph object.

    Returns
    -------

    A sparse matrix in csr format.

    NOTE: duplicated from "nbd.py" to avoid excessive imports.

    """
    # Note: the rows of nx.adjacency_matrix(graph) are in the same order as
    # the list returned by graph.nodes().
    degrees = graph.degree()
    degrees = sp.diags([degrees[n] for n in graph.nodes()])
    adj = nx.adjacency_matrix(graph)
    ident = sp.eye(graph.order())
    pseudo = sp.bmat([[None, degrees - ident], [-ident, adj]])
    return pseudo.asformat('csr')