Change integration in DOS.density() to make careful integration over nonlinear meshes

[aeantipov]

Hi all,

I was playing with DOS class and the HilbertTransform. To check that my dos is normalized I used DOS.density() method with mu as a big number. Right now the integration is done in a following way ( in pytriqs/DOS/DOS.py)

de = self.eps[1]-self.eps[0]
dens = (sum(self.rho[0:ind]) - self.rho[0]/2.0 - self.rho[ind-1]/2.0) * de

A better solution is to write something like

dens = numpy.trapz( x = self.eps[0:ind], y = self.rho[0:ind] )

which I attach as a pull request.

[aeantipov]

To check the results I used the following script ( A[0].eps is a logarithmic grid with several dense points ):

f=lambda x : exp(-x*x)/sqrt(pi)
G=DOS_from_function(lambda x : exp(-x*x)/sqrt(pi), xmin = -6, xmax = 6, Npts = len(A[0].eps))
val=[]
for eps in A[0].eps:
    val.append(f(eps))
G2 = DOS(eps = A[0].eps, rho = val)
print "Regular mesh     : ", G.density(10)
print "Logarithmic mesh : ", G2.density(10)

Before was:

Regular mesh    :  1.0
Logarithmic mesh :  68.0356596281

And with a patch I get

Regular mesh     :  1.0
Logarithmic mesh :  1.00002354885

[mferrero]

Before pulling this in the main branch, I need to make sure the modification is OK with the rest of the framework. Here's the issue: the DOS object has mainly been created to be used with the Hilbert transform. Your correction will indeed cure the density() method when you use a non-linear mesh. However, if you were to try to compute a Hilbert transform, it would not work because a linear mesh is still expected. So I think that one should do the things right and go all the way by modifying the code so that also the Hilbert transforms works with non-linear meshes.