!
!-----------------------------------------------------------------------
subroutine solve_linter ( dvbare, dvscf, xq, et, vr, w)
!-----------------------------------------------------------------------
!
! this works for one perturbation at a time
!
!-----------------------------------------------------------------------
!
use parameters
use constants
use gspace
use kspace
implicit none
!
complex(kind=DP) :: dvbare(nr)
! the perturbation in real space
complex(kind=DP) :: dvscf(nr), vr(nr)
real(DP) :: et(nbnd_occ, nks), w
!
complex(kind=DP) :: vr_dyn (nr)
! local potential plus teh dynamicla part w + i * eta
!
integer :: ik, ikk, ikq, iter, ibnd, jbnd, ios, ig, ir
complex(kind=DP) :: evc (ngm, nbnd_occ), evq (ngm, nbnd_occ)
real(kind=DP) :: g2kin(ngm), dr2, wgt, qg2, xq(3)
!
complex(kind=DP) :: dpsi(ngm,nbnd_occ), dvpsi(ngm,nbnd_occ), &
dvpsi0(ngm,nbnd_occ), dvscfin(nr), hpsi(ngm), &
hpsi2(ngm,nbnd_occ), ps2(nbnd_occ,nbnd_occ), &
dpsi0(ngm,nbnd_occ), drhoscf(nr), dvscfout(nr)
complex(kind=DP) :: aux(nr), aux1(nr), aux2(nr)
complex(kind=DP) :: ps(nbnd_occ), auxg(ngm)
complex(kind=DP) :: ZDOTC
real(DP) :: eprec(nbnd_occ), h_diag(ngm, nbnd_occ), anorm, meandvb
logical :: conv_root, convt
integer :: lter
external ch_psi_all, ch_psi_all_eta
!
! include the dynamnical part inside the local potential
! (which is already complex)
!
vr_dyn = vr - w
!
!
! loop over the iterations
!
iter = 0
convt = .false.
do while (iter.lt.nmax_iter .and. .not.convt)
!
iter = iter + 1
drhoscf = czero
!
do ik = 1, nksq
!
ikk = 2 * ik - 1
ikq = ikk + 1
!
! reads unperturbed wavefuctions u(k) and u(k+q)
!
read ( iunwfc, rec = ikk, iostat = ios) evc
read ( iunwfc, rec = ikq, iostat = ios) evq
!
! compute the kinetic energy for k+q
!
do ig = 1, ngm
g2kin(ig) = ( (xk(1,ikq) + g(1,ig))**2.d0 + &
(xk(2,ikq) + g(2,ig))**2.d0 + &
(xk(3,ikq) + g(3,ig))**2.d0 ) * tpiba2
enddo
!
if (iter.eq.1) then
!
do ibnd = 1, nbnd_occ
!
! dpsi and dvscfin are set to zero
!
dpsi = czero
dvscfin = czero
!
! dvbare*psi is calculated for this k point and all bands...
! (we compute the product in real space)
!
aux = czero
do ig = 1, ngm
aux ( nl ( ig ) ) = evc (ig, ibnd)
enddo
call cfft3 ( aux, nr1, nr2, nr3, 1)
do ir = 1, nr
aux (ir) = aux(ir) * dvbare (ir)
enddo
! back to G-space (fft order of G-vectors)
call cfft3 ( aux, nr1, nr2, nr3, -1)
! switch to magnitude-order of G-vectors
do ig = 1, ngm
dvpsi(ig, ibnd) = aux( nl(ig) )
enddo
!
enddo
!
! writes dvpsi for this k-point on iunit iubar
!
write ( iubar, rec = ik, iostat = ios) dvpsi
!
else
!
! read dvbare*psi for this k-point on iunit iubar
!
read ( iubar, rec = ik, iostat = ios) dvpsi
!
! dvpsi = dvbare*psi + dvscfin*psi
!
do ibnd = 1, nbnd_occ
!
aux = czero
do ig = 1, ngm
aux ( nl ( ig ) ) = evc (ig, ibnd)
enddo
call cfft3 ( aux, nr1, nr2, nr3, 1)
do ir = 1, nr
aux (ir) = aux (ir) * dvscfin (ir)
enddo
call cfft3 ( aux, nr1, nr2, nr3, -1)
do ig = 1, ngm
dvpsi (ig, ibnd) = dvpsi (ig, ibnd) + aux ( nl(ig) )
enddo
!
enddo
!
! read dpsi for this k-point from iudwf
!
read ( iudwf, rec = ik, iostat = ios) dpsi
!
endif
!
! ( 1 - P_occ^{k+q} ) * dvpsi
!
do ibnd = 1, nbnd_occ
auxg = czero
do jbnd = 1, nbnd_occ
ps(jbnd) = - ZDOTC(ngm, evq(:,jbnd), 1, dvpsi(:,ibnd), 1)
call ZAXPY (ngm, ps (jbnd), evq (:, jbnd), 1, auxg, 1)
enddo
call DAXPY (2 * ngm, one, auxg, 1, dvpsi (:, ibnd), 1)
enddo
!
! change the sign of the known term
!
call DSCAL (2 * ngm * nbnd_occ, - 1.d0, dvpsi, 1)
!
! iterative solution of the linear system
! (H-et)*dpsi= - ( 1 - P_occ^{k+q} ) * (dvbare+dvscf)*psi
! [ dvpsi ]
!
! !
! ! this preconditioner sounds complicated...
! !
! do ibnd = 1, nbnd_occ
! do ig = 1, ngm
! auxg (ig) = g2kin (ig) * evq (ig, ibnd)
! enddo
! eprec (ibnd) = 1.35d0 * ZDOTC (ngm, evq (1, ibnd), 1, auxg, 1)
! enddo
! do ibnd = 1, nbnd_occ
! do ig = 1, ngm
! h_diag(ig,ibnd)=1.d0/max(1.0d0,g2kin(ig)/eprec(ibnd))
! enddo
! enddo
!
! this seems to work even better...
!
h_diag=1.d0
!
!
! now obtain dpsi from cgsolve_all
!
!@ call cgsolve_all (ch_psi_all, et(:,ikk), dvpsi, dpsi, h_diag, &
!@ ngm, ngm, tr_cgsolve, ik, lter, conv_root, anorm, nbnd_occ, &
!@ g2kin, vr_dyn, evq )
call bcgsolve_all (ch_psi_all_eta, et(:,ikk), dvpsi, dpsi, h_diag, &
ngm, ngm, tr_cgsolve, ik, lter, conv_root, anorm, nbnd_occ, &
g2kin, vr_dyn, evq )
! if (.not.conv_root) &
! write( 6, '(4x,"ik",i4," linter: one or more roots not converged ",e10.3)') &
! ik , anorm
! write(6,'("cgsolve_all:",2x,3i5,3x,e9.3)') iter, ik, lter, anorm
!
! !
! ! DEBUG: calculate (H-et+alpha*Pv)*dpsi-dvpsi, this should be zero
! ! if dpsi is the correct solution - O.K. this is checked
! !
! call ZGEMM ('C', 'N', nbnd_occ , nbnd_occ, ngm, (1.d0, 0.d0) , evq, &
! ngm, dpsi, ngm, (0.d0, 0.d0) , ps2, nbnd_occ)
! call ZGEMM ('N', 'N', ngm, nbnd_occ, nbnd_occ, dcmplx(alpha_pv,0.d0), evq, &
! ngm, ps2, nbnd_occ, czero, hpsi2, ngm)
! do ibnd = 1, nbnd_occ
! call h_psi_c ( dpsi(:,ibnd), hpsi, g2kin, vr_dyn)
! hpsi = hpsi - et(ibnd,ikk) * dpsi(:,ibnd) - dvpsi(:,ibnd)
! hpsi = hpsi + hpsi2 (:, ibnd)
! write(6,*) '--------------------'
! do ig = 1, ngm
! write(6,'(3i5,3(2x,2f15.5))') &
! ik, ibnd, ig, 100*dpsi(ig,ibnd), 100*dvpsi(ig,ibnd), 100*hpsi(ig)
! enddo
! enddo
!
! writes dpsi for this k point on iunit iudwf
!
write ( iudwf, rec = ik, iostat = ios) dpsi
!
! contribution to drhoscf from this kpoint
!
wgt = 2.d0 * wk(ikk) / omega
!
do ibnd = 1, nbnd_occ
!
aux1 = czero
do ig = 1, ngm
aux1 ( nl ( ig ) ) = evc (ig, ibnd)
enddo
call cfft3 ( aux1, nr1, nr2, nr3, 1)
!
aux2 = czero
do ig = 1, ngm
aux2 ( nl ( ig ) ) = dpsi (ig, ibnd)
enddo
call cfft3 ( aux2, nr1, nr2, nr3, 1)
!
do ir = 1, nr
drhoscf (ir) = drhoscf (ir) + wgt * conjg (aux1 (ir) ) * aux2 (ir)
enddo
!
enddo
!
enddo
!
! here we have drhoscf for this iteration
! compute the corresponding Hartree potential (RPA)
!
dvscfout = czero
call cfft3 ( drhoscf, nr1, nr2, nr3, -1)
!
! here we enforce zero average variation of the charge density
! if the bare perturbation does not have a constant term
! (otherwise the numerical error, coupled with a small denominator
! in the coulomb term, gives rise to a spurious dvscf response)
!
meandvb = sqrt ( (sum(dreal(dvbare)))**2.d0 + (sum(aimag(dvbare)))**2.d0 ) / float(nr)
if (meandvb.lt.1.d-8) drhoscf ( nl(1) ) = 0.d0
!
do ig = 1, ngm
qg2 = (g(1,ig)+xq(1))**2 + (g(2,ig)+xq(2))**2 + (g(3,ig)+xq(3))**2
if (qg2 > 1.d-8) &
dvscfout ( nl(ig) ) = e2 * fpi * drhoscf ( nl(ig) ) / (tpiba2 * qg2)
enddo
!
call cfft3 ( dvscfout, nr1, nr2, nr3, 1)
!
! we mix with the old potential
!
! modif broyden mixing, complex potential
!
call mix_potential_c ( nr, dvscfout, dvscfin, alpha_mix, dr2, tr2_ph, iter, nmix_ph, convt)
!
! convergence criterion
!
convt = dr2.lt.tr2_ph
!
write(6,'(4x, "scf iteration ",i3,": dr2 = ",e8.2)') iter, dr2
!
enddo
!
! at this point dvscfin is the converged Hartree screening.
! the screened coulomb interaction corresponds to dv_bare + dv_hartree (RPA)
!
dvscf = dvscfin + dvbare
!
return
end subroutine solve_linter
!
!-----------------------------------------------------------------------
!