! !----------------------------------------------------------------------- subroutine green ( ik, g2kin, vr) !----------------------------------------------------------------------- ! ! 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 ! ! ik is the index of this k-point in the {k0-q} list ! !----------------------------------------------------------------------- ! 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, ih, ik ! ! 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), phi (ngm) ! real(DP) :: a(nstep+2,4), b(nstep+2,4), norm(4), a_term(4), b_term(4) ! a(:,1) are the coefficients for ! a(:,2) are the coefficients for ! ... ! norm(1) is the norm ^1/2 ! ... ! a_term(1) is the constant coefficient for the terminator for ! ... ! complex(DP) :: prefac(4) ! complex(DP) :: ZDOTC real(DP) :: norm1, w integer :: ig, iw, ibnd complex(DP) :: gr, gr_tmp ! ! variables for debugging ! real(DP) :: eval(nbnd) complex(DP) :: ctmp1, ctmp2, gr_exp(nw), psi(ngm,nbnd) ! prefac = (/ cone, -cone, ci, -ci/) ! a = zero b = zero a_term = zero b_term = zero norm = zero ! do ig1 = 1, 2 !@ ngms do ig2 = 1, 2 !@ ngms ! ! define planewaves for projecting the Green's function ! ! d1 = exp(i*G1*r), d2 = exp(i*G2*r) ! d1 = czero d1 ( ig1 ) = cone d2 = czero d2 ( ig2 ) = cone ! ! check normalization, should be 1 (OK) ! ! norm1 = DREAL ( ZDOTC (ngm, d1, 1, d1, 1) ) ! norm1 = DREAL ( ZDOTC (ngm, d2, 1, 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) do iw = 1, nw w = wmin + ( wmax - wmin ) / float (nw-1) * float(iw-1) gr_exp(iw) = gr_exp(iw) + conjg(ctmp1)*ctmp2 / ( w - eval(ibnd)*ryd2ev + ci * delta) enddo enddo do iw = 1, nw w = wmin + ( wmax - wmin ) / float (nw-1) * float(iw-1) write(200,'(3f15.10)') w, gr_exp(iw) enddo ! ! Haydock's formula for off-diagonal elements [insert ref. to book] ! ! G_12 = G_uu - G_vv - i * (G_ww - G_zz) ! G_21 = 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 ! ! NORMALIZATION: If I use a non-normalized function ! Haydock does not work. Keep this in mind for the ! local-orbital implementation ! ! the 4 Haydock sequences ! do ih = 1, 4 if (ih.ne.2 .or. ig1.ne.ig2) then ! ! in the case ih=2 we need to remove the spurious ! solution corresponding to ig1=ig2 (phi = d1-d2 = 0) ! phi = 0.5d0 * ( d1 + prefac(ih) * d2 ) call normalize ( phi, norm(ih) ) call haydock ( g2kin, vr, phi, a(:,ih), b(:,ih), a_term(ih), b_term(ih)) endif enddo ! call rw_haydock ( ik, ig1, ig2, a, b, a_term, b_term, norm, +1) ! ! enddo enddo ! ! do ig1 = 1, 2 !@ ngms do ig2 = 1, 2 !@ ngms ! ! call rw_haydock ( ik, ig1, ig2, a, b, a_term, b_term, norm, -1) ! do iw = 1, nw ! w = wmin + ( wmax - wmin ) / float (nw-1) * float(iw-1) gr = czero ! do ih = 1, 4 ! call recfrac ( w, a(:,ih), b(:,ih), a_term(ih), b_term(ih), gr_tmp) ! ! in the case ih=2 we need to remove the spurious ! solution corresponding to ig1=ig2 (phi = d1-d2 = 0) ! if ( ih.ne.2 .or. ig1.ne.ig2 ) & gr = gr + prefac(ih) * gr_tmp * norm(ih) * norm(ih) ! enddo ! write(201,'(3f15.10)') w, gr ! enddo ! enddo enddo ! return end subroutine green ! !----------------------------------------------------------------------- subroutine haydock ( g2kin, vr, un, a, b, a_term, b_term) !----------------------------------------------------------------------- ! 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, b_term ! do iw = 1, nw w(iw) = wmin + ( wmax - wmin ) / float (nw-1) * float(iw-1) enddo ! ! 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 ! ! 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) ! ! eigenvalues of the tridiagonal 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 ! a_term = 0.5d0 * ( D(1) + D(N) ) b_term = 0.25d0 * ( D(N) - D(1) ) ! ! write(44,*) a_term, b_term ! endif ! 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,4f15.5)') a(istep+1), b(istep+1), a_term, b_term ! enddo ! return end subroutine haydock ! !----------------------------------------------------------------------- subroutine recfrac ( w, a, b, a_term, b_term, gr) !----------------------------------------------------------------------- ! ! the recursive fraction ! ! Note that the use of the terminator does remove some wild oscillations ! at high energy, but it cannot fix in any way the location of the poles. ! I have the impression that in order to converge the Green's function ! to the energy E, we need to apply the Hamiltonian (# steps in Haydock) ! the number of times required to reach the eigenvalue close to E. ! use constants use parameters implicit none ! complex(DP) :: gr real(DP) :: a(nstep+2), b(nstep+2), w, a_term, b_term integer :: i integer, parameter :: nterm = 300 ! nterm is the number of fractions used to generate the terminator real(DP) :: a1, b2 ! gr = czero ! ! the terminator ! do i = nterm+nstep,nstep+1,-1 gr = w + cmplx(zero, delta) - a_term * ryd2ev - ( ryd2ev * b_term ) ** two * gr gr = cone / gr enddo ! ! the actual fraction ! do i = nstep,0,-1 gr = w + cmplx(zero, delta) - a(i+1) * ryd2ev - ( ryd2ev * b(i+2) ) ** two * gr gr = cone / gr enddo ! return end subroutine recfrac ! !----------------------------------------------------------------------- subroutine normalize ( phi, norm ) !----------------------------------------------------------------------- ! ! simply normalize a vector in G-space ! use constants use parameters use gspace implicit none ! complex(DP) :: phi(ngm) real(DP) :: norm complex(DP) :: ZDOTC ! norm = DREAL ( ZDOTC (ngm, phi, 1, phi, 1) ) norm = sqrt(norm) call DSCAL (2 * ngm, one / norm, phi, 1) ! return end subroutine normalize !