! !----------------------------------------------------------------------- subroutine green ( ig1, ig2, g2kin, vr, gr) !----------------------------------------------------------------------- ! ! Green's function from Haydock's recursion: G_k (G,G',w) for a given k ! the k-dependence is inside the kinetic energy g2kin ! ! for some reasons the dot_product(conjg(),()) does not work ! need to check this out ! !----------------------------------------------------------------------- ! use gspace use constants use parameters implicit none ! complex(kind=DP) :: vr (nr) real(kind=DP) :: g2kin (ngm) integer :: ig1, ig2, i, j, k, ir ! ! notation as in Haydock - review (pag 227) ! e.g.: bnpunp = b_{n+1} * u_{n+1}, hun = H * u_n ! complex(DP) :: d1 (ngm), d2 (ngm), u (ngm) ! complex(DP) :: ZDOTC real(DP) :: freq(nw), norm integer :: ig, iw, ibnd complex(DP) :: gr(nw), gr_xx(nw) ! ! variables for debugging ! real(DP) :: eval(nbnd) complex(DP) :: ctmp1, ctmp2, gr_exp(nw), psi(ngm,nbnd) ! ! frequency bins for the Green's function ! do iw = 1, nw freq(iw) = wmin + ( wmax - wmin ) / float (nw-1) * float(iw-1) enddo ! ! ! DEFINE PLANEWAVES ! ! d1 = exp(i*G1*r), d2 = exp(i*G2*r) ! d1 = czero d1 ( ig1 ) = cone d2 = czero d2 ( ig2 ) = cone ! ! ! DEBUG: check normalization, should be 1 (OK) ! ! ! norm = DREAL ( ZDOTC (ngm, d1, 1, d1, 1) ) ! norm = DREAL ( ZDOTC (ngm, d2, 1, d2, 1) ) ! call DSCAL (2 * ngm, one / norm, d2, 1) ! ! ! ! ! ! DEBUG - calculate the Green's function using the sum over states ! ! ! ! ! call eigenstates_all ( vr, g2kin, psi, eval ) ! gr_exp = czero ! do ibnd = 1, nbnd ! ctmp1 = ZDOTC (ngm, psi(:,ibnd), 1, d1, 1) ! ctmp2 = ZDOTC (ngm, psi(:,ibnd), 1, d2, 1) ! gr_exp = gr_exp + conjg(ctmp1)*ctmp2 / ( freq - eval(ibnd)*ryd2ev + ci * delta) ! enddo ! do iw = 1, nw ! write(200,'(3f15.10)') freq(iw), gr_exp(iw) ! enddo ! ! ! ! G_12 = G_uu - G_vv - i * (G_ww - G_zz) ! ! u = 0.5 * (d1+ d2) ! v = 0.5 * (d1- d2) ! w = 0.5 * (d1+id2) ! z = 0.5 * (d1-id2) ! ! NOTE: when ig1=ig2 the procedure at (d1-d2)/2 ! finds a nihil eigenvalue and the ! Green's function is peaked at 0 - spurious ! ! IF I DON'T NORMALIZE IT DOESN NOT WORK - this is due to the fact that I am ! using wrong psinc functions - these are not orthogonal! ! gr = czero ! u = 0.5d0 * ( d1 + d2 ) norm = DREAL ( ZDOTC (ngm, u, 1, u, 1) ) norm = sqrt(norm) call DSCAL (2 * ngm, one / norm, u, 1) call haydock ( g2kin, vr, freq, u, gr_xx) gr_xx = norm * norm * gr_xx gr = gr + gr_xx ! if (ig1.ne.ig2) then ! ! trivial solution - spurios peak at 0 in the GF ! u = 0.5d0 * ( d1 - d2 ) norm = DREAL ( ZDOTC (ngm, u, 1, u, 1) ) norm = sqrt(norm) call DSCAL (2 * ngm, one / norm, u, 1) call haydock ( g2kin, vr, freq, u, gr_xx) gr_xx = norm * norm * gr_xx gr = gr - gr_xx ! endif ! u = 0.5d0 * ( d1 + ci * d2 ) norm = DREAL ( ZDOTC (ngm, u, 1, u, 1) ) norm = sqrt(norm) call DSCAL (2 * ngm, one / norm, u, 1) call haydock ( g2kin, vr, freq, u, gr_xx) gr_xx = norm * norm * gr_xx gr = gr - ci * gr_xx ! u = 0.5d0 * ( d1 - ci * d2 ) norm = DREAL ( ZDOTC (ngm, u, 1, u, 1) ) norm = sqrt(norm) call DSCAL (2 * ngm, one / norm, u, 1) call haydock ( g2kin, vr, freq, u, gr_xx) gr_xx = norm * norm * gr_xx gr = gr + ci * gr_xx ! do iw = 1, nw write(201,'(3f15.10)') freq(iw), gr(iw) enddo ! return end subroutine green ! !----------------------------------------------------------------------- subroutine haydock ( g2kin, vr, w, un, gr) !----------------------------------------------------------------------- ! use gspace use constants use parameters implicit none ! complex(kind=DP) :: vr (nr) real(kind=DP) :: g2kin (ngm) integer :: i, j, k, ir ! ! notation as in Haydock - review (pag 227) ! e.g.: bnpunp = b_{n+1} * u_{n+1}, hun = H * u_n complex(DP) :: unp (ngm), un (ngm), unm (ngm), hun (ngm), bnpunp (ngm) ! real(DP) :: w(nw), an, bn, bnp, a(nstep+2), b(nstep+2), bnp2 complex(DP) :: gr(nw), ZDOTC, DZNRM2 integer :: ig, istep, iw, nterm complex(DP) :: Z, aux real(DP) :: D(nstep+1), E(nstep), WORK, tmp integer :: INFO, N ! ! variables for the terminator ! real(DP) :: a_term, b2_term, term(nw) ! ! ! initialize the recursion ! a = zero b = zero bn = zero unm = czero ! do istep = 0, nstep ! ! basic Haydock scheme [Eqs. (5.2)-(5.13) of SSP] ! call h_psi_c ( un, hun, g2kin, vr) an = DREAL ( ZDOTC (ngm, un, 1, hun, 1) ) ! ! bnpunp = hun - an * un - bn * unm ! unp = bunpunp / norm(bnpunp) ! ! NOTE it is crytical to enter ZAXPY with complex coefficients ! otehrwise unbelievable things may happen... ! aux = dcmplx(-an,0) call ZAXPY (ngm, aux, un , 1, hun, 1) aux = dcmplx(-bn,0) call ZAXPY (ngm, aux, unm, 1, hun, 1) bnp = DREAL ( ZDOTC (ngm, hun, 1, hun, 1) ) bnp = sqrt (bnp) call DSCAL (2 * ngm, one / bnp, hun, 1) ! call ZCOPY (ngm, un , 1, unm, 1) ! unm <- un call ZCOPY (ngm, hun, 1, un , 1) ! un <- unp ! a (istep+1) = an b (istep+2) = bnp bn = bnp ! ! if (istep.gt.0) then ! ! ! ! DEBUG: eigenvalues of the Hamiltonian at this iteration ! ! ! write(6,'(/4x,a)') repeat('-',67) ! write(6,'(4x,"Eigenvalues (eV) of the tridiagonal matrix at iteration",i5)') & ! istep + 1 ! write(6,'(4x,a/)') repeat('-',67) ! ! ! N = istep + 1 ! D = zero ! E = zero ! D(1:N) = a(1:N) ! E(1:N-1) = b(2:N) ! call DSTEV( 'N', N, D, E, Z, 1, WORK, INFO ) ! write(6,'(4x,7f9.3)') D(1:N) * ryd2ev ! ! ! ! define the terminator using the extrema of the eigenvalue spectrum ! ! (this should be much more stable than trying to fit the fluctuating ! ! coefficients to constant values) - Giannozzi et al, Appl. Num. Math. ! ! 4, 273 (1988) - Eq. (7.1) ! ! ! a_term = 0.5d0 * ( D(1) + D(N) ) * ryd2ev ! b2_term = ( 0.25d0 * ( D(N) - D(1) ) * ryd2ev ) ** 2.d0 ! ! ! endif ! ! ! ! DEBUG: green's function at each step ! ! ! gr = czero ! ! ! do i = istep,0,-1 ! gr = w + cmplx(zero, delta) - a(i+1) * ryd2ev - ( ryd2ev * b(i+2) ) ** two * gr ! gr = cone / gr ! enddo ! do iw = 1, nw ! write(100,'(i5,2f15.10)') istep, w(iw), - one/pi * aimag ( gr(iw) ) ! enddo ! write(100,*) ! ! enddo ! ! Green's function projected on the initial state ! gr = czero ! ! ! ! start with a costant chain defined by a_term and b2_term (already in eV) ! ! ! ! THIS IS BAD ! IF I LOOK AT gr GIVEN ONLY BY THE TERMINATOR IT'S A MESS... ! WHAT THE HELL !!! ! ! nterm = 500 ! do i = 1, nterm ! gr = w + cmplx(zero, delta) - a_term - b2_term * gr ! gr = cone / gr ! enddo ! do i = nstep,0,-1 gr = w + cmplx(zero, delta) - a(i+1) * ryd2ev - ( ryd2ev * b(i+2) ) ** two * gr gr = cone / gr enddo ! ! ! ! write down the recursion coefficients ! ! ! write(6,'(/4x,a)') repeat('-',67) ! write(6,'(4x,"Recursion coefficients")') ! write(6,'(4x,a/)') repeat('-',67) ! do istep = 0, nstep ! write(6,'(4x,2f15.5)') a(istep+1), b(istep+1) ! enddo ! return end subroutine haydock !