Python scipy.linalg.eigvalsh() Examples

def dos1d_ewindow(h,energies=np.linspace(-1.,1.,30),delta=None,info=False,
                    use_green=True,nk=300):
  """Calculate the density of states in certain energy window"""
  ys = [] # density of states
  if delta is None: # pick a good delta value
    delta = 0.1*(max(energies) - min(energies))/len(energies)
  if True: # do not use green function    
    import scipy.linalg as lg
    kxs = np.linspace(0.,1.,nk)
    hkgen= h.get_hk_gen() # get hamiltonian generator
    ys = energies*0.
    weight = 1./(nk)
    for ix in kxs:
      hk = hkgen(ix) # get hk hamiltonian
      evals = lg.eigvalsh(hk) # get eigenvalues
      ys += weight*calculate_dos(evals,energies,delta) # add this contribution
    if info: print("Done",ix)
    write_dos(energies,ys) # write in file
    return 

def eigvalsh(m):
    """Wrapper for linalg"""
    m = todense(m) # turn the matrix dense
    if np.max(np.abs(m.imag))<error: m = m.real # real matrix
    if not accelerate: return dlg.eigvalsh(m)
    # check if doing slices helps
    n = m.shape[0] # size of the matrix
    mo = m[0:n:2,1:n:2] # off diagonal is zero
#    if False: # assume block diagonal
    if np.max(np.abs(mo))<error: # assume block diagonal
        # detected block diagonal
        es0 = dlg.eigvalsh(m[0:n:2,0:n:2]) # recall
        es1 = dlg.eigvalsh(m[1:n:2,1:n:2]) # recall
        es = np.concatenate([es0,es1]) # concatenate array
        return es

    else:
      if np.max(np.abs(m.imag))<error: # assume real
          return dlg.eigvalsh(m.real) # diagonalize real matrix
      else: return dlg.eigvalsh(m) # diagonalize complex matrix 

def kmesh_density(h,nk=20,delta=0.01):
  """Function to compute a k-mesh density districution"""
  if h.dimensionality !=2: raise
  hkgen = h.get_hk_gen() # get generator
  out = [] # empty list
  xout = []
  yout = []
  for x in np.linspace(-0.1,1.1,nk,endpoint=True):
    for y in np.linspace(-0.1,1.1,nk,endpoint=True):
      hk = hkgen([x,y])
      es = lg.eigvalsh(hk) # eigenvalues
      e = np.min(es[es>0.]) - np.max(es[es<0.]) # gap
      out.append(e) # append gap
      xout.append(x)
      yout.append(y)
  # now create the function that interpolates
  import interpolation
  f = interpolation.interpolator2d(xout,yout,out) # interpolator function
  # now define a function that returns a probability density
  def fout(k):
    return delta/(delta**2 + f(k[0],k[1])) # output
  return fout # return function 

def diagonalize(h,nkpoints=100):
  """ Diagonalice a hamiltonian """
  import scipy.linalg as lg
  # for one dimensional systems
  if h.dimensionality==1:  # one simensional system
    klist = np.arange(0.0,1.0,1.0/nkpoints)  # create array with the klist
    if h.geometry.shift_kspace:
      klist = np.arange(-0.5,0.5,1.0/nkpoints)  # create array with the klist
    intra = h.intra  # assign intraterm
    inter = h.inter  # assign interterm
    energies = [] # list with the energies
    for k in klist: # loop over kpoints
      bf = np.exp(1j*np.pi*2.*k)  # bloch factor for the intercell terms
      inter_k = inter*bf  # bloch hopping
      hk = intra + inter_k + inter_k.H # k dependent hamiltonian
      energies += [lg.eigvalsh(hk)] # get eigenvalues of the current hamiltonian
    energies = np.array(energies).transpose() # each kpoint in a line
    return (klist,energies) # return the klist and the energies
# for zero dimensional systems system
  elif h.dimensionality==0:  
    intra = h.intra  # assign intraterm
    energies = lg.eigvalsh(intra) # get eigenvalues of the current hamiltonian
    return (range(len(intra)),energies) # return indexes and energies
  else: raise 

def diagonalize_kpath(h,kpath):
  """Diagonalice in a certain path"""
  energies = [] # empty list with energies
  import scipy.linalg as lg
  ik = 0.
  iks = [] # empty list
  for k in kpath:
    f = h.get_hk_gen() # get Hk generator
    hk = f(k)  # k dependent hamiltonian
    es = (lg.eigvalsh(hk)).tolist() # get eigenvalues for current hamiltonian
    energies += es # append energies 
    iks += [ik for i in es]
    ik += 1.
  iks = np.array(iks)
  iks = iks/max(iks) # normalize path
  return (iks,energies) 

def eigenvalues(h0,nk=10,notime=False):
    """Return all the eigenvalues of a Hamiltonian"""
    from . import klist
    h = h0.copy() # copy hamiltonian
    h.turn_dense()
    ks = klist.kmesh(h.dimensionality,nk=nk) # get grid
    hkgen = h.get_hk_gen() # get generator
    if parallel.cores==1:
      es = [] # empty list
      if not notime: est = timing.Testimator(maxite=len(ks))
      for k in ks: # loop
          if not notime: est.iterate()
          es += algebra.eigvalsh(hkgen(k)).tolist() # add
    else:
        f = lambda k: algebra.eigvalsh(hkgen(k)) # add
        es = parallel.pcall(f,ks) # call in parallel
        es = np.array(es)
        es = es.reshape(es.shape[0]*es.shape[1])
    return es # return all the eigenvalues 

def degenerate_points(h,n=0):
    """Return the points in the Brillouin zone that have a node
    in the bandstructure"""
    from scipy.optimize import differential_evolution
    bounds = [(0.,1.) for i in range(h.dimensionality)]
    hk_gen = h.get_hk_gen() # generator
    def get_point(x0):
      def f(k): # conduction band eigenvalues
        hk = hk_gen(k) # Hamiltonian
        es = lg.eigvalsh(hk) # get eigenvalues
        return abs(es[n]-es[n+1]) # gap
      res = differential_evolution(f,bounds=bounds) # minimize
      return res.x
    x0 = np.random.random(h.dimensionality) # inital vector
    return get_point(x0) # get the k-point 

def _check_precision_matrix(precision, covariance_type):
    """Check a precision matrix is symmetric and positive-definite."""
    if not (np.allclose(precision, precision.T) and
            np.all(linalg.eigvalsh(precision) > 0.)):
        raise ValueError("'%s precision' should be symmetric, "
                         "positive-definite" % covariance_type) 

def _validate_covars(covars, covariance_type, n_components):
    """Do basic checks on matrix covariance sizes and values."""
    from scipy import linalg
    if covariance_type == 'spherical':
        if len(covars) != n_components:
            raise ValueError("'spherical' covars have length n_components")
        elif np.any(covars <= 0):
            raise ValueError("'spherical' covars must be non-negative")
    elif covariance_type == 'tied':
        if covars.shape[0] != covars.shape[1]:
            raise ValueError("'tied' covars must have shape (n_dim, n_dim)")
        elif (not np.allclose(covars, covars.T)
              or np.any(linalg.eigvalsh(covars) <= 0)):
            raise ValueError("'tied' covars must be symmetric, "
                             "positive-definite")
    elif covariance_type == 'diag':
        if len(covars.shape) != 2:
            raise ValueError("'diag' covars must have shape "
                             "(n_components, n_dim)")
        elif np.any(covars <= 0):
            raise ValueError("'diag' covars must be non-negative")
    elif covariance_type == 'full':
        if len(covars.shape) != 3:
            raise ValueError("'full' covars must have shape "
                             "(n_components, n_dim, n_dim)")
        elif covars.shape[1] != covars.shape[2]:
            raise ValueError("'full' covars must have shape "
                             "(n_components, n_dim, n_dim)")
        for n, cv in enumerate(covars):
            if (not np.allclose(cv, cv.T)
                    or np.any(linalg.eigvalsh(cv) <= 0)):
                raise ValueError("component %d of 'full' covars must be "
                                 "symmetric, positive-definite" % n)
    else:
        raise ValueError("covariance_type must be one of " +
                         "'spherical', 'tied', 'diag', 'full'") 

def laplacian_spectrum(G, weight='weight'):
    """Return eigenvalues of the Laplacian of G

    Parameters
    ----------
    G : graph
       A NetworkX graph

    weight : string or None, optional (default='weight')
       The edge data key used to compute each value in the matrix.
       If None, then each edge has weight 1.

    Returns
    -------
    evals : NumPy array
      Eigenvalues

    Notes
    -----
    For MultiGraph/MultiDiGraph, the edges weights are summed.
    See to_numpy_matrix for other options.

    See Also
    --------
    laplacian_matrix
    """
    from scipy.linalg import eigvalsh
    return eigvalsh(nx.laplacian_matrix(G, weight=weight).todense()) 

def normalized_laplacian_spectrum(G, weight='weight'):
    """Return eigenvalues of the normalized Laplacian of G

    Parameters
    ----------
    G : graph
       A NetworkX graph

    weight : string or None, optional (default='weight')
       The edge data key used to compute each value in the matrix.
       If None, then each edge has weight 1.

    Returns
    -------
    evals : NumPy array
      Eigenvalues

    Notes
    -----
    For MultiGraph/MultiDiGraph, the edges weights are summed.
    See to_numpy_matrix for other options.

    See Also
    --------
    normalized_laplacian_matrix
    """
    from scipy.linalg import eigvalsh
    return eigvalsh(nx.normalized_laplacian_matrix(G, weight=weight).todense()) 

def laplacian_spectrum(G, weight='weight'):
    """Returns eigenvalues of the Laplacian of G

    Parameters
    ----------
    G : graph
       A NetworkX graph

    weight : string or None, optional (default='weight')
       The edge data key used to compute each value in the matrix.
       If None, then each edge has weight 1.

    Returns
    -------
    evals : NumPy array
      Eigenvalues

    Notes
    -----
    For MultiGraph/MultiDiGraph, the edges weights are summed.
    See to_numpy_matrix for other options.

    See Also
    --------
    laplacian_matrix
    """
    from scipy.linalg import eigvalsh
    return eigvalsh(nx.laplacian_matrix(G, weight=weight).todense()) 

def eigenvalues(h,nk):
    """Return all the eigenvalues of a Hamiltonian"""
    import klist
    h.turn_dense()
    ks = klist.kmesh(h.dimensionality,nk=nk) # get grid
    hkgen = h.get_hk_gen() # get generator
    e0 = 0.0
    import timing
    est = timing.Testimator(maxite=len(ks))
    for k in ks: # loop
      est.iterate()
      es = eigvalsh(hkgen(k)).tolist() # add
      e0 += np.sum(es[es<0.]) # add contribution
    return e0/len(ks) # return the ground state energy

#
#def gsenergy(h,nk=10):
#  """Calculate the ground state energy"""
#  if h.dimensionality!=2: raise  # continue if two dimensional
#  hk_gen = h.get_hk_gen() # gets the function to generate h(k)
#  kxs = np.linspace(0.,1.,nk)  # generate kx
#  kys = np.linspace(0.,1.,nk)  # generate ky
#  e0 = 0.0 # initalize
#  for x in kxs:
#    for y in kys:
#      hk = hk_gen([x,y]) # get hamiltonian
#      es = eigvalsh(hk) # eigenvalues
#      e0 += np.sum(es[es<0.]) # add contribution
#  return e0/nk**2 

def boolean_fermi_surface(h,write=True,output_file="BOOL_FERMI_MAP.OUT",
                    e=0.0,nk=50,nsuper=1,reciprocal=False,
                    delta=None):
  """Calculates the Fermi surface of a 2d system"""
  if h.dimensionality!=2: raise  # continue if two dimensional
  hk_gen = h.get_hk_gen() # gets the function to generate h(k)
  kxs = np.linspace(-nsuper,nsuper,nk)  # generate kx
  kys = np.linspace(-nsuper,nsuper,nk)  # generate ky
  kdos = [] # empty list
  kxout = []
  kyout = []
  if reciprocal: R = h.geometry.get_k2K() # get matrix
  # setup a reasonable value for delta
  if delta is None:
    delta = 8./np.max(np.abs(h.intra))/nk
  for x in kxs:
    for y in kxs:
      r = np.matrix([x,y,0.]).T # real space vectors
      k = np.array((R*r).T)[0] # change of basis
      hk = hk_gen(k) # get hamiltonian
      evals = lg.eigvalsh(hk) # diagonalize
      de = np.abs(evals - e) # difference with respect to fermi
      de = de[de<delta] # energies close to fermi
      if len(de)>0: kdos.append(1.0) # add to the list
      else: kdos.append(0.0) # add to the list
      kxout.append(x)
      kyout.append(y)
  if write:  # optionally, write in file
    f = open(output_file,"w") 
    for (x,y,d) in zip(kxout,kyout,kdos):
      f.write(str(x)+ "   "+str(y)+"   "+str(d)+"\n")
    f.close() # close the file
  return (kxout,kyout,d) # return result 

def _validate_covars(covars, covariance_type, n_components):
    """Do basic checks on matrix covariance sizes and values
    """
    from scipy import linalg
    if covariance_type == 'spherical':
        if len(covars) != n_components:
            raise ValueError("'spherical' covars have length n_components")
        elif np.any(covars <= 0):
            raise ValueError("'spherical' covars must be non-negative")
    elif covariance_type == 'tied':
        if covars.shape[0] != covars.shape[1]:
            raise ValueError("'tied' covars must have shape (n_dim, n_dim)")
        elif (not np.allclose(covars, covars.T)
              or np.any(linalg.eigvalsh(covars) <= 0)):
            raise ValueError("'tied' covars must be symmetric, "
                             "positive-definite")
    elif covariance_type == 'diag':
        if len(covars.shape) != 2:
            raise ValueError("'diag' covars must have shape "
                             "(n_components, n_dim)")
        elif np.any(covars <= 0):
            raise ValueError("'diag' covars must be non-negative")
    elif covariance_type == 'full':
        if len(covars.shape) != 3:
            raise ValueError("'full' covars must have shape "
                             "(n_components, n_dim, n_dim)")
        elif covars.shape[1] != covars.shape[2]:
            raise ValueError("'full' covars must have shape "
                             "(n_components, n_dim, n_dim)")
        for n, cv in enumerate(covars):
            if (not np.allclose(cv, cv.T)
                    or np.any(linalg.eigvalsh(cv) <= 0)):
                raise ValueError("component %d of 'full' covars must be "
                                 "symmetric, positive-definite" % n)
    else:
        raise ValueError("covariance_type must be one of " +
                         "'spherical', 'tied', 'diag', 'full'") 

def distance_riemann(A, B):
    """Riemannian distance between two covariance matrices A and B.

    .. math::
            d = {\left( \sum_i \log(\lambda_i)^2 \\right)}^{-1/2}

    where :math:`\lambda_i` are the joint eigenvalues of A and B

    :param A: First covariance matrix
    :param B: Second covariance matrix
    :returns: Riemannian distance between A and B

    """
    return numpy.sqrt((numpy.log(eigvalsh(A, B))**2).sum()) 

def laplacian_spectrum(G, weight='weight'):
    """Return eigenvalues of the Laplacian of G

    Parameters
    ----------
    G : graph
       A NetworkX graph

    weight : string or None, optional (default='weight')
       The edge data key used to compute each value in the matrix.
       If None, then each edge has weight 1.

    Returns
    -------
    evals : NumPy array
      Eigenvalues

    Notes
    -----
    For MultiGraph/MultiDiGraph, the edges weights are summed.
    See to_numpy_matrix for other options.

    See Also
    --------
    laplacian_matrix
    """
    from scipy.linalg import eigvalsh
    return eigvalsh(nx.laplacian_matrix(G,weight=weight).todense()) 

def _validate_covars(covars, covariance_type, n_components):
    """Do basic checks on matrix covariance sizes and values."""
    from scipy import linalg
    if covariance_type == 'spherical':
        if len(covars) != n_components:
            raise ValueError("'spherical' covars have length n_components")
        elif np.any(covars <= 0):
            raise ValueError("'spherical' covars must be positive")
    elif covariance_type == 'tied':
        if covars.shape[0] != covars.shape[1]:
            raise ValueError("'tied' covars must have shape (n_dim, n_dim)")
        elif (not np.allclose(covars, covars.T)
              or np.any(linalg.eigvalsh(covars) <= 0)):
            raise ValueError("'tied' covars must be symmetric, "
                             "positive-definite")
    elif covariance_type == 'diag':
        if len(covars.shape) != 2:
            raise ValueError("'diag' covars must have shape "
                             "(n_components, n_dim)")
        elif np.any(covars <= 0):
            raise ValueError("'diag' covars must be positive")
    elif covariance_type == 'full':
        if len(covars.shape) != 3:
            raise ValueError("'full' covars must have shape "
                             "(n_components, n_dim, n_dim)")
        elif covars.shape[1] != covars.shape[2]:
            raise ValueError("'full' covars must have shape "
                             "(n_components, n_dim, n_dim)")
        for n, cv in enumerate(covars):
            if (not np.allclose(cv, cv.T)
                    or np.any(linalg.eigvalsh(cv) <= 0)):
                raise ValueError("component %d of 'full' covars must be "
                                 "symmetric, positive-definite" % n)
    else:
        raise ValueError("covariance_type must be one of " +
                         "'spherical', 'tied', 'diag', 'full'")


# Copied from scikit-learn 0.19. 

def _check_precision_matrix(precision, covariance_type):
    """Check a precision matrix is symmetric and positive-definite."""
    if not (np.allclose(precision, precision.T) and
            np.all(linalg.eigvalsh(precision) > 0.)):
        raise ValueError("'%s precision' should be symmetric, "
                         "positive-definite" % covariance_type) 

def _validate_covars(covars, covariance_type, n_components):
    """Do basic checks on matrix covariance sizes and values."""
    from scipy import linalg
    if covariance_type == 'spherical':
        if len(covars) != n_components:
            raise ValueError("'spherical' covars have length n_components")
        elif np.any(covars <= 0):
            raise ValueError("'spherical' covars must be non-negative")
    elif covariance_type == 'tied':
        if covars.shape[0] != covars.shape[1]:
            raise ValueError("'tied' covars must have shape (n_dim, n_dim)")
        elif (not np.allclose(covars, covars.T)
              or np.any(linalg.eigvalsh(covars) <= 0)):
            raise ValueError("'tied' covars must be symmetric, "
                             "positive-definite")
    elif covariance_type == 'diag':
        if len(covars.shape) != 2:
            raise ValueError("'diag' covars must have shape "
                             "(n_components, n_dim)")
        elif np.any(covars <= 0):
            raise ValueError("'diag' covars must be non-negative")
    elif covariance_type == 'full':
        if len(covars.shape) != 3:
            raise ValueError("'full' covars must have shape "
                             "(n_components, n_dim, n_dim)")
        elif covars.shape[1] != covars.shape[2]:
            raise ValueError("'full' covars must have shape "
                             "(n_components, n_dim, n_dim)")
        for n, cv in enumerate(covars):
            if (not np.allclose(cv, cv.T)
                    or np.any(linalg.eigvalsh(cv) <= 0)):
                raise ValueError("component %d of 'full' covars must be "
                                 "symmetric, positive-definite" % n)
    else:
        raise ValueError("covariance_type must be one of " +
                         "'spherical', 'tied', 'diag', 'full'") 

def _validate_covars(covars, covariance_type, n_components):
    '''Do basic checks on matrix covariance sizes and values
    '''
    if covariance_type == 'spherical':
        if len(covars) != n_components:
            raise ValueError("'spherical' covars have length n_components")
        elif np.any(covars <= 0):
            raise ValueError("'spherical' covars must be non-negative")
    elif covariance_type == 'tied':
        if covars.shape[0] != covars.shape[1]:
            raise ValueError(
                "'tied' covars must have shape (n_dimensions, n_dimensions)")
        elif (not np.allclose(covars, covars.T)
              or np.any(linalg.eigvalsh(covars) <= 0)):
            raise ValueError("'tied' covars must be symmetric, "
                             "positive-definite")
    elif covariance_type == 'diag':
        if len(covars.shape) != 2:
            raise ValueError("'diag' covars must have shape "
                             "(n_components, n_dimensions)")
        elif np.any(covars <= 0):
            raise ValueError("'diag' covars must be non-negative")
    elif covariance_type == 'full':
        if len(covars.shape) != 3:
            raise ValueError("'full' covars must have shape "
                             "(n_components, n_dimensions, n_dimensions)")
        elif covars.shape[1] != covars.shape[2]:
            raise ValueError("'full' covars must have shape "
                             "(n_components, n_dimensions, n_dimensions)")
        for n, cv in enumerate(covars):
            if (not np.allclose(cv, cv.T)
                    or np.any(linalg.eigvalsh(cv) <= 0)):
                raise ValueError("component %d of 'full' covars must be "
                                 "symmetric, positive-definite" % n)
    else:
        raise ValueError("covariance_type must be one of " +
                         "'spherical', 'tied', 'diag', 'full'") 

def _check_precision_matrix(precision, covariance_type):
    """Check a precision matrix is symmetric and positive-definite."""
    if not (np.allclose(precision, precision.T) and
            np.all(linalg.eigvalsh(precision) > 0.)):
        raise ValueError("'%s precision' should be symmetric, "
                         "positive-definite" % covariance_type) 

def eigenvalues(hetero,numeig=10,effective=False,gf=None,full=False):
  """ Gets the lowest eigenvalues of the central part of the hamiltonian"""
  if not hetero.block_diagonal:
    print(""" heterounction in eigenvalues must be block diagonal""")
    raise
  # if effective hamiltonian, just calculate the eigenvalues
  if effective: # effective hamiltonian
    print("Calculating eigenvalues of effective hamiltonian...")
    if hetero.heff is None:  #if hasn't been calculated so far
      effective_central_hamiltonian(hetero,write=False)
    heff = hetero.heff # store the list of list with central EFF ham
    import scipy.sparse.linalg as lg
    eig,eigvec = lg.eigs(heff,k=numeig,which="LM",sigma=0.0)
    return eig

  from scipy.sparse import csc_matrix,bmat
  import scipy.sparse.linalg as lg
  if not effective: # do not use the effective hamiltonian, only the central
    intra = hetero.central_intra # store the list of list with central ham
  numb = len(intra) # number of central blocks
  intrasp = [[None for i in range(numb)] for j in range(numb)]
  # assign onsite and couplings in sparse form
  if not hetero.is_sparse:
    for i in range(numb):
      intrasp[i][i] = csc_matrix(intra[i][i])
    for i in range(numb-1):
      intrasp[i][i+1] = csc_matrix(intra[i][i+1])
      intrasp[i+1][i] = csc_matrix(intra[i+1][i])
  if hetero.is_sparse:
    for i in range(numb):
      intrasp[i][i] = intra[i][i]
    for i in range(numb-1):
      intrasp[i][i+1] = intra[i][i+1]
      intrasp[i+1][i] = intra[i+1][i]
  intrasp = bmat(intrasp) # create sparse matrix
  if effective:
    evals,evecs = lg.eigs(intrasp,k=numeig,which="LM",sigma=0.0)
  if not effective:
    if full:  # full diagonalization
      from scipy.linalg import eigvalsh
      evals = eigvalsh(intrasp.todense())
    else:
      evals,evecs = lg.eigsh(intrasp,k=numeig,which="LM",sigma=0.0)
  return evals 

def calculate_dos_hkgen(hkgen,ks,ndos=100,delta=None,
         is_sparse=False,numw=10,window=None,energies=None):
  """Calculate density of states using the ks given on input"""
  if not is_sparse: # if not is sparse
      m = hkgen([0,0,0]) # get the matrix
      if algebra.issparse(m): 
          print("Hamiltonian is sparse, selecting sparse mode in DOS")
          is_sparse = True
  es = np.zeros((len(ks),hkgen(ks[0]).shape[0])) # empty list
  tr = timing.Testimator("DOS",maxite=len(ks))
  if delta is None: delta = 5./len(ks) # automatic delta
  from . import parallel
  def fun(k): # function to execute
    if parallel.cores==1: tr.iterate() # print the info
    hk = hkgen(k) # Hamiltonian
    t0 = time.perf_counter() # time
    if is_sparse: # sparse Hamiltonian 
      es = algebra.smalleig(hk,numw=numw,tol=delta/1e3) # eigenvalues
      ws = np.zeros(es.shape[0])+1.0 # weight
    else: # dense Hamiltonian
      es = algebra.eigvalsh(hk) # get eigenvalues
      ws = np.zeros(es.shape[0])+1.0 # weight
    return es # return energies
#  for ik in range(len(ks)):  
  out = parallel.pcall(fun,ks) # launch all the processes
  es = [] # empty list
  for o in out: 
      es = np.concatenate([es,o]) # concatenate
#    tr.remaining(ik,len(ks))
#  es = es.reshape(len(es)*len(es[0])) # 1d array
  es = np.array(es) # convert to array
  nk = len(ks) # number of kpoints
  if energies is not None: # energies given on input
    xs = energies
  else:
    if window is None:
      xs = np.linspace(np.min(es)-.5,np.max(es)+.5,ndos) # create x
    else:
      xs = np.linspace(-window,window,ndos) # create x
  ys = calculate_dos(es,xs,delta) # use the Fortran routine
  ys /= nk # normalize by the number of k-points
  ys *= 1./np.pi # normalization of the Lorentzian
  write_dos(xs,ys) # write in file
  print("\nDOS finished")
  return (xs,ys) # return result 

def dist(self, G1, G2, normalization=None, timescales=None):
        """NetLSD: Hearing the Shape of a Graph.

        A network similarity measure based on spectral node signature
        distributions.

        The results dictionary includes the underlying signature vectors in
        `'signatures'`.

        Parameters
        ----------

        G1, G2 (nx.Graph)
            two undirected networkx graphs to be compared.

        normalization (str)
            type of normalization of the heat kernel vectors. either
            `'complete'`, `'empty'` or `'none'`

        timescales (np.ndarray)
            timescales for the comparison. None yields default.

        Returns
        -------

        dist (float)
            the distance between `G1` and `G2`.

        References
        ----------

        .. [1] A. Tsitsulin, D. Mottin, P. Karras, A. Bronstein &
               E. Müller. NetLSD: Hearing the Shape of a Graph. KDD 2018

        """
        if normalization is None:
            normalization = 'none'
        if timescales is None:
            timescales = np.logspace(-2, 2, 256)
        assert isinstance(
            normalization, str
        ), 'Normalization parameter must be of string type'

        lap1 = nx.normalized_laplacian_matrix(G1)
        lap2 = nx.normalized_laplacian_matrix(G2)

        # Note: this is O(n^3) worst-case.
        eigs1 = spl.eigvalsh(lap1.todense())
        eigs2 = spl.eigvalsh(lap2.todense())

        hkt1 = _lsd_signature(eigs1, timescales, normalization)
        hkt2 = _lsd_signature(eigs2, timescales, normalization)

        self.results['signatures'] = (hkt1, hkt2)
        self.results['dist'] = np.linalg.norm(hkt1 - hkt2)

        return self.results['dist']