::::::::::::::
gwdfpt.note.6.11.08
::::::::::::::

 I don't think we can achieve linear scaling as suggested
 in Baroni&Giannozzi92. In their case the kinetic energy
 was calculated by finite differences, so it's strictly
 localized in R-space.
 Using psinc functions instead, the kinetic energy is still
 localized, but it non-vanishing on further neighbors.
 Mostofi et al. mention that the application of the kinetic
 energy is still a global operation to be performed with
 FFTs (http://www.tcm.phy.cam.ac.uk/onetep/intro/node6.html).
 In other words, if we want to keep the G-space accuracy, 
 we NEED to use the full representation of the kinetic
 energy WITHOUT truncation.
 This means that it is useless to compute the matrix 
 elements of the kinetic energy in the psinc basis, since
 they don't decay so fast.
 For Haydock, we still need a localized object to start
 with, and the psinc choice seems still to be the most
 convenient to switch to/from the G-space.

::::::::::::::
gwdfpt.note.6.12.08
::::::::::::::

 In order to have a complete mapping of the psinc basis into the G-space,
 we need to keep the entire G-bx, and not to trunctate it within the
 kinetic energy cutoff !!

 All the nice properties of the psinc's are based on this full mapping,
 if we truncate the orthogonality and all that vanish.
 
 I think we need to define a grid spacing such that the corresponding
 G-box is CONTAINED within the sphere of the G-vectors smaller than
 the cutoff. In practice the calculation HAS TO BE RUN with a cutoff
 smaller than for the wavefunctions. I think this is still fine, since
 calculations of epsilon are usually done with 5 or 10 Ry.

::::::::::::::
gwdfpt.note.6.17.08
::::::::::::::

 Ok I think I am able to calculate the eigenvalues of H with the recursion method 
 for the special case k=0 and psinc(1) (diagonal element, hermitian Hamiltonian).
 The eigenvalues do coincide with those of the CG method to within 1 meV after
 about 150 iterations.

 I obtained the eigenvalues by diagonalizing the tridiagonal matrix obtained from
 the recursion coefficients. Some non0degenerate eigenvalues appear multiple times,
 I don't know at this stage whether this is a problem for the computation of the
 Green's function or not. i would say no, since in the recursion methods there is
 no assumption about the multiplicity of the calculated spectrum.

 The Green's function has poles at SOME OF the eigenvalues, the strenght is nonzero
 when the initial state has components on that eigenstate.

 Now I need to try the nondiagonal elements, and the nonzero k vectors.

::::::::::::::
gwdfpt.note.6.18.08
::::::::::::::

 I tested the Green's function G(r,r',k,w) against the sum-over-state expression.
  
 The diagonal element on psinc(1) and k=0 gives EXACTLY the same result as for the
 sum over states (I compared the imaginary parts to avoid calculating many empty states).
 This is great, now I need to do the off-diagonal matrix elements and the case k/=0.

::::::::::::::
gwdfpt.note.6.19.08
::::::::::::::

 I tested the Green's function G(r,r',k,w) against the sum-over-state expression.

 OK, the off-diagonal elements work perfectly.
 Also the case k/=0 works perfectly.

 I need about 100 iterations to have the Green's function converged with delta = 0.1 eV.

 It looks like the Green's function converges from the bottom, i.e. the more iterations
 we do, the larger the energy range where it matches the full calculation based on
 the sum-over-states.