# 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
#
# Authors: Michel Ferrero, Olivier Parcollet, Priyanka Seth, Nils Wentzell, tayral

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()) :
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)