Source code for triqs.dos.dos

# Copyright (c) 2013-2016 Commissariat à l'énergie atomique et aux énergies alternatives (CEA)
# Copyright (c) 2013-2016 Centre national de la recherche scientifique (CNRS)
# Copyright (c) 2020 Simons Foundation
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You may obtain a copy of the License at
#     https:#www.gnu.org/licenses/gpl-3.0.txt
#
# Authors: Thomas Ayral, Michel Ferrero, Olivier Parcollet, Priyanka Seth, Nils Wentzell


import types,string,itertools
import numpy

[docs] class DOS : r""" * Stores a density of state of fermions .. math:: \rho (\epsilon) \equiv \sum'_k \delta( \epsilon - \epsilon_k) * The sum is normalized .. math:: \int_{-\infty}^{\infty} d\epsilon \rho (\epsilon) = 1 * Implement :ref:`Plot Protocol <plotting>`. """
[docs] def __init__(self, eps, rho, name = ''): """ Parameters ------------ eps : 1d array-type eps[i] is value of epsilon. rho : 1d array-type The corresponding value of the dos. name : string Name of the dos/orbital """ self.name = name try : self.eps = numpy.array( eps ) assert len(self.eps.shape) ==1 except : raise RuntimeError("Argument eps mismatch") try : self.rho = numpy.array( rho ) assert len(self.rho.shape) ==1 except : raise RuntimeError("Argument rho mismatch") assert self.eps.shape[0] == self.rho.shape[0], "Dimensions of eps and rho do not match"
#------------------------------------------------------------- def __reduce__(self) : return self.__class__, (self.eps,self.rho, self.name) def __reduce_to_dict__(self) : return {'epsilon' : self.eps, 'rho': self.rho} @classmethod def __factory_from_dict__(cls,name, D) : return cls(D['epsilon'],D['rho'], name) def __repr__(self) : return """ DOS object : """%self.__dict__
[docs] def copy(self): return DOS(self.eps,self.rho,self.name)
def _plot_(self, Options) : return [ {'label' : self.name, 'xlabel' :r'$\epsilon$', 'ylabel' : r'%s$(\epsilon)$'%self.name, 'xdata' : self.eps,'ydata' : self.rho } ]
[docs] def density(self,mu=0): """Calculates the density of free fermions for the given DOS for chemical potential mu.""" dens = 0.0 a = [ (e>mu) for e in self.eps ] try: ind = a.index(True) except: ind = self.eps.shape[0] de = self.eps[1]-self.eps[0] #for e,r in itertools.izip(self.eps[0:ind],self.rho[0:ind]): # dens += r dens = (sum(self.rho[0:ind]) - self.rho[0]/2.0 - self.rho[ind-1]/2.0) * de #dens2 = dens + (self.rho[ind-1]/2.0 + self.rho[ind]/2.0) * de if (ind<self.eps.shape[0]): dens += (mu-self.eps[ind-1]) * (self.rho[ind-1] + self.rho[ind])/2.0 return dens
##########################################################################
[docs] def dos_from_file(Filename, name = '', single_orbital = None): """ Read the DOS from a file :param Filename: a string : name of the file :param name: name of the DOS :param single_orbital: can be None or an integer. :rtype: * if single_orbital== None, returns a tuple of DOS (even if there is one dos !). * If single_orbital==i, return only ONE DOS corresponding to ith orbital (starting at 1). Format of the file : * N_orbitals +1 columns, * the first column is the value of epsilon * the N_orbitals other columns are the values of the dos for various orbitals """ f = open(Filename); s='' while not(s.strip()) : s= f.readline() assert s, "File is empty !" N_Orbitals = len (s.split()) - 1 assert N_Orbitals >0, "File : wrong format" # not very safe : fromfile routine can crashes if given non numerics r = numpy.fromfile(Filename,sep=' ') l,div = r.shape[0], N_Orbitals +1 assert l%(div)==0,"File does not contains N*%d numbers !"%(div) r.shape = l//(div) , div # reshape the array eps = r[:,0] if single_orbital : assert single_orbital>0 and single_orbital <= N_Orbitals, " single_orbital " return DOS (r[:,0] ,r[:,single_orbital], name) else : return [ DOS (r[:,0] ,r[:,i +1 ], name) for i in range (N_Orbitals)]
##########################################################################
[docs] class DOSFromFunction(DOS): """ * A DOS class, but constructed from a function. * The number of points can be variable and self-adjusted in the Hilbert transform to adapt precision. """
[docs] def __init__(self, function, x_min, x_max, n_pts=100, name=''): """ :param function: * a function :math:`\\epsilon \\rightarrow \\rho(\\epsilon)` * The result type can be anything from which a 1d-array can be constructed by numpy :param x_min,x_max: Bound of the mesh (domain of the function). :param n_pts: Number of points in the mesh. :param name: Name of the DOS. """ assert callable(function), "function is not callable" self.function,self.x_min,self.x_max = function,x_min,x_max try : e = function(0.001) len(numpy.array(e).shape) ==1 except : raise RuntimeError("Value of the function must be a 1d-array") self.__f(n_pts) # compute arrays DOS.__init__(self,self.eps,self.rho,name)
#------------------------------------------------------------- def __reduce__(self) : return self.__class__, (self.function,self.x_min, self.x_max, len(self.eps), self.name) #------------------------------------------------------------- def __f(self,N) : r = (self.x_max - self.x_min)/float(N-1) self.eps = numpy.array( [self.x_min + r* i for i in range(N) ] ) self.rho = numpy.array( [self.function(e) for e in self.eps])
#----------------------------------------------------- # Register the class for HDFArchive #----------------------------------------------------- from h5.formats import register_class register_class (DOS)