!
  !-----------------------------------------------------------------------
  subroutine solve_linter_dyn ( ir0, dvbare, dvscf, xxq, et, vrs, w, maxscf, maxbcgsolve, alpha_mix)
  !-----------------------------------------------------------------------
  !
  ! this works for one perturbation at a time
  !
  ! Fully converge dpsi for each one of the first "maxscf" dvscf iterations,    
  ! then perform only "maxbcgsolve" steps for dpsi at each dvscf iteration.
  ! This should be similar to the "direct minimization method" in Paratec
  ! setting maxscf = 10**10 we have the conventional minimization
  ! (full relaxation at every scf step)
  !
  !-----------------------------------------------------------------------
  !
  use parameters
  use constants
  use gspace
  use kspace
  implicit none
  !
  integer :: maxscf, maxbcgsolve
  integer :: maxbcgsolve0
  !
  complex(kind=DP) :: dvbare(nrs)
  ! the perturbation in real space
  complex(kind=DP) :: dvscf(nrs), vrs(nrs)
  real(DP) :: et(nbnd_occ, nks), w, alpha_mix
  !
  complex(kind=DP) :: vrs_dyn (nrs)
  ! local potential plus the dynamical part w + i * eta
  !
  integer :: ik, ikk, ikq, iter, ibnd, jbnd, ios, ig, ir, ir0
  integer :: i2,j2,k2,i1,j1,k1
  complex(kind=DP) :: evc (ngm, nbnd_occ), evq (ngm, nbnd_occ)
  real(kind=DP) :: g2kin(ngm), dr2, dr2_old, wgt, qg2, xxq(3), x
  !
  complex(kind=DP) :: dpsi(ngms,nbnd_occ), dvpsi(ngms,nbnd_occ), &
                      dvpsi0(ngms,nbnd_occ), dvscfin(nrs), hpsi(ngm), &
                      hpsi2(ngm,nbnd_occ), ps2(nbnd_occ,nbnd_occ), &
                      drhoscf(nrs), dvscfout(nrs)
  complex(kind=DP) :: dpsip(ngms,nbnd_occ), dpsim(ngms,nbnd_occ)
  complex(kind=DP) :: aux(nrs), aux1(nrs), aux2(nrs)
  complex(kind=DP) :: ps(nbnd_occ), auxg(ngms)
  complex(kind=DP) :: ZDOTC
  real(DP) :: eprec(nbnd_occ), h_diag(ngms, nbnd_occ), anorm, meandvb
  logical :: conv_root, convt, panic, neighbor
  integer :: lter
  external ch_psi_all, ch_psi_all_eta

  maxbcgsolve0 = maxbcgsolve
  panic = .false.
  dr2_old = 1.d10
  !
  if (ir0.gt.1) then
    k2 = (ir0-1)/nr1/nr2+1
    j2 = 1+(ir0-(k2-1)*nr1*nr2-1)/nr1
    i2 = ir0 - ( (j2-1) * nr1 + (k2-1) * nr1 * nr2  )
    !
    k1 = (ir0-1-1)/nr1/nr2+1
    j1 = 1+(ir0-1-(k1-1)*nr1*nr2-1)/nr1
    i1 = ir0-1 - ( (j1-1) * nr1 + (k1-1) * nr1 * nr2  )
    !
    if ( (abs(i1-i2).le.1) .and. (abs(j1-j2).le.1) .and. ((abs(k1-k2).le.1)) ) then
      neighbor = .true.
      write(6,*) 'neighbors: ',i1,j1,k1,i2,j2,k2
    endif
    !
  else
    neighbor = .false.
  endif
  !
  !  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) .and. .not.neighbor ) then
          !
          do ibnd = 1, nbnd_occ
            !
            !  dpsi and dvscfin are set to zero
            !
            dpsip = czero
            dpsim = 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, ngms
              aux ( nls ( ig ) ) = evc (ig, ibnd)
            enddo
            call cfft3s ( aux, nr1s, nr2s, nr3s,  1)
            do ir = 1, nrs
              aux (ir) = aux(ir) * dvbare (ir)
            enddo
            ! back to G-space (fft order of G-vectors)
            call cfft3s ( aux, nr1s, nr2s, nr3s, -1)
            ! switch to magnitude-order of G-vectors
            do ig = 1, ngms
              dvpsi(ig, ibnd) = aux( nls(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, ngms
              aux ( nls ( ig ) ) = evc (ig, ibnd)
            enddo
            call cfft3s ( aux, nr1s, nr2s, nr3s,  1)
            do ir = 1, nrs
               aux (ir) = aux (ir) * dvscfin (ir)
            enddo
            call cfft3s ( aux, nr1s, nr2s, nr3s, -1)
            do ig = 1, ngms
              dvpsi (ig, ibnd) = dvpsi (ig, ibnd) + aux ( nls(ig) )
            enddo
            !
          enddo
          !
          !  read dpsi for this k-point from iudwf
          !
          read ( iudwfp, rec = ik, iostat = ios) dpsip
          read ( iudwfm, rec = ik, iostat = ios) dpsim
          !
        endif
        !
        !  ( 1 - P_occ^{k+q} ) * dvpsi
        !
        do ibnd = 1, nbnd_occ
           auxg = czero
           do jbnd = 1, nbnd_occ
              ps(jbnd) = - ZDOTC(ngms, evq(:,jbnd), 1, dvpsi(:,ibnd), 1)
              call ZAXPY (ngms, ps (jbnd), evq (:, jbnd), 1, auxg, 1)
           enddo
           call DAXPY (2 * ngms, one, auxg, 1, dvpsi (:, ibnd), 1)
        enddo
        !
        !  change the sign of the known term
        !
        call DSCAL (2 * ngms * nbnd_occ, - 1.d0, dvpsi, 1)
        !
        ! iterative solution of the linear system 
        ! (H-et)*dpsi=   - ( 1 - P_occ^{k+q} ) * (dvbare+dvscf)*psi
        !              [                  dvpsi                     ]

        !
        ! I use the preconditioner of Teter, Payne, Allan [PRB 40, 12255 (1988)]
        !
        do ibnd = 1, nbnd_occ
           do ig = 1, ngms
               auxg (ig) = g2kin (ig) * evq (ig, ibnd)
           enddo
           eprec (ibnd) = 1.35d0 * ZDOTC (ngms, evq (1, ibnd), 1, auxg, 1)
        enddo
        do ibnd = 1, nbnd_occ
           do ig = 1, ngms
      !       x = g2kin(ig)/eprec(ibnd)
      !       h_diag(ig,ibnd) = (27.d0+18.d0*x+12.d0*x*x+8.d0*x**3.d0) &
      !                       / (27.d0+18.d0*x+12.d0*x*x+8.d0*x**3.d0+16.d0*x**4.d0)
              ! original preconditioning used in PH 
              !
      !       h_diag(ig,ibnd) = 1.d0/max(1.d0,x) 
              ! 
              ! NO preconditioning
              !
              h_diag(ig,ibnd) = 1.d0
              !
           enddo
        enddo
        !
        ! --- now obtain dpsi from cBiCG
        !
        !
        ! include the dynamnical part inside the local potential 
        ! (which is already complex)
        !
        vrs_dyn = vrs - w 
        !
!@@        if ( ( iter.le.maxscf ) .or. panic ) then 
          call bcgsolve_all   (ch_psi_all_eta, et(:,ikk), dvpsi, dpsip, h_diag, &
               ngms, ngms, tr_cgsolve, ik, lter, conv_root, anorm, nbnd_occ, &
               g2kin, vrs_dyn, evq )
 !@@       else
 !@@         call bcgsolve_all_fixed (ch_psi_all_eta, et(:,ikk), dvpsi, dpsip, h_diag, &
 !@@              ngms, ngms, tr_cgsolve, ik, lter, conv_root, anorm, nbnd_occ, &
 !@@              g2kin, vrs_dyn, evq, maxbcgsolve0 )
  !@@      endif
        !
        vrs_dyn = vrs + w
        !
        if ( ( iter.le.maxscf ) .or. panic ) then 
          call bcgsolve_all   (ch_psi_all_eta, et(:,ikk), dvpsi, dpsim, h_diag, &
               ngms, ngms, tr_cgsolve, ik, lter, conv_root, anorm, nbnd_occ, &
               g2kin, vrs_dyn, evq )
        else
          call bcgsolve_all_fixed (ch_psi_all_eta, et(:,ikk), dvpsi, dpsim, h_diag, &
               ngms, ngms, tr_cgsolve, ik, lter, conv_root, anorm, nbnd_occ, &
               g2kin, vrs_dyn, evq, maxbcgsolve0 )
        endif
        !
        dpsi = 0.5d0 * ( dpsip + dpsim )
        !


!       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 ( iudwfp, rec = ik, iostat = ios) dpsip
        write ( iudwfm, rec = ik, iostat = ios) dpsim
        !
        ! contribution to drhoscf from this kpoint
        !
        wgt = 2.d0 * wk(ikk) / omega
        !
        do ibnd = 1, nbnd_occ
          !
          aux1 = czero
          do ig = 1, ngms
            aux1 ( nls ( ig ) ) = evc (ig, ibnd)
          enddo
          call cfft3s ( aux1, nr1s, nr2s, nr3s,  1)
          !
          aux2 = czero
          do ig = 1, ngms
            aux2 ( nls ( ig ) ) = dpsi (ig, ibnd)
          enddo
          call cfft3s ( aux2, nr1s, nr2s, nr3s,  1)
          !
          do ir = 1, nrs
            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 cfft3s ( drhoscf, nr1s, nr2s, nr3s, -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(nrs)
     if (meandvb.lt.1.d-8) drhoscf ( nls(1) ) = 0.d0
     !
     do ig = 1, ngms
       qg2 = (g(1,ig)+xxq(1))**2 + (g(2,ig)+xxq(2))**2 + (g(3,ig)+xxq(3))**2
       if (qg2 > 1.d-8) &
         dvscfout ( nls(ig) ) =  e2 * fpi * drhoscf ( nls(ig) ) / (tpiba2 * qg2)

     enddo
     !
     call cfft3s ( dvscfout, nr1s, nr2s, nr3s,  1)
     !
     ! we mix with the old potential
     !
     ! modif broyden mixing, complex potential
     !
     call mix_potential_c ( nrs, 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
     !
     if (dr2.gt.2.d0*dr2_old) panic = .true.
     !
     dr2_old = 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_dyn
  !
  !-----------------------------------------------------------------------
  !