module  m_rsibl
  use m_ll,only:ll
  use m_hsibl,only: hsibl,hsibl1
  public rsibl,rsibl_ev
  private
contains
  subroutine rsibl(lfrce,isp,q,iq,ndimh,nspc,napw,igapw,nevec,evec,ewgt,k1,k2,k3,smpot,smrho,f) !Add smooth part of output density into smrho and forces.
    use m_struc_def  
    use m_lmfinit,only: lat_alat,nspec,nbas
    use m_lattic,only: lat_qlat, lat_vol,lat_plat
    use m_supot,only: lat_nabc, lat_ng, lat_gmax
    use m_uspecb,only:uspecb
    !i Inputs
    !i   lfrce :if nonzero, accumulate contribution to force
    !i   nbas  :size of basis
    !i   sspec :struct for species-specific information; see routine uspec
    !i     Elts read: ngcut
    !i     Stored:    *
    !i     Passed to: tbhsi rsibl1 uspecb
    !i   lfrce :1 calculate contribution to forces
    !i   nbas  :size of basis
    !i   q     :Bloch vector
    !i   iq    :index to current k-point
    !i   ndimh :dimension of hamiltonian
    !i   nspc  :2 for coupled spins; otherwise 1
    !i   napw  :number of augmented PWs in basis
    !i   igapw :vector of APWs, in units of reciprocal lattice vectors
    !i   nevec :number of eigenvectors with nonzero weights
    !i   evec  :eigenvectors
    !i   ewgt  :eigenvector weights
    !i   k1..3 :dimensions smpot,smrho
    !i   smpot :smooth potential on uniform mesh, needed for forces
    !o Outputs
    !o   smrho :smooth density accumulated for this qp
    !o   f     :force contribution accumulated for this qp
    !r Remarks
    !m MPI
    !m   Parallelise over the eigenvector loop. The vector block size is
    !m   chosen (in the range 6-16, by dstrbp.f) so as to distribute the
    !m   work optimally across processes. Two work arrays of the size of
    !m   smrho are allocated from the heap as buffers. Only one will be
    !m   needed under MPI-2. See comments in hsibl.
    !b Bugs
    !b    replace call to gvgvcomp and pass ipv as input
    !b    The non-F90 version should work, but it is no longer tested
    !u Updates
    !u   29 Dec 08 Unsuccessful attempt to make work with openmp
    !u   05 Jul 08 (T. Kotani) output density for new PW part
    !u   10 Sep 06 Added MPI parallelization in the spin-coupled case
    !u   23 Dec 04 Extended to spin-coupled case
    !u   25 Aug 04 Adapted to extended local orbitals
    !u   10 Apr 02 Redimensionsed eh,rsmh to accomodate larger lmax
    !u   15 Feb 02 (ATP) Added MPI parallelization
    !u   27 Aug 01 Extended to local orbitals.
    !u   12 Oct 00 Use q-dependent list of G vectors
    !u    6 Jul 00 attempt to vectorize by grouping eigenvectors in blocks
    !u   17 Jun 00 Spin polarized
    !u   23 May 00 Adapted from nfp rsif_q.f
    ! ----------------------------------------------------------------------
    implicit none
    integer :: procid, master, nproc, mpipid
    integer :: lfrce,isp,k1,k2,k3,ndimh,nevec,iq,nspc
    integer :: napw,igapw(3,napw)
    real(8):: q(3) , ewgt(nevec) , f(3,nbas)
    complex(8):: evec(ndimh,nspc,nevec),smrho(k1,k2,k3,isp), smpot(k1,k2,k3,isp)
    ! ... Local parameters
    integer :: n0,nkap0,nermx,npmx,nblk,nlmto
    parameter (n0=10,nkap0=3,nermx=100,npmx=128)
    integer:: ngabc(3) , n1 , n2 , n3 , nrt , net , ng , &
         nglob , ltop , nlmtop , ogq , og2 , ohe , ohr , oyl , oylw , oiv , iprint
    integer,allocatable :: iv_a_okv(:)
    real(8),allocatable :: rv_a_ogv(:,:)
    equivalence (n1,ngabc(1)),(n2,ngabc(2)),(n3,ngabc(3))
    integer :: iprt(n0,nkap0,nermx),ipet(n0,nkap0,nermx),i_copy_size
    double precision :: alat,qlat(3,3),plat(3,3),q0(3),gmax,xx
    complex(8):: w(1,1,1)!dummy 
    real(8):: wff(1)     !dummy
    double precision :: vol
    double precision :: etab(nermx),rtab(nermx)
    integer :: ivec,nvec
    integer,allocatable:: ivp(:)
    complex(8),allocatable::psi(:,:,:),psir(:,:,:),vpsi(:,:,:), &
         wk(:,:,:)
    real(8),allocatable:: cosi(:),sini(:),wk2(:)
    integer:: ivecini,ivecend,ig
    integer,allocatable:: w_oiv(:)
    real(8),allocatable:: w_ogq(:),w_oyl(:),w_oylw(:),w_og2(:),w_ohe(:),w_ohr(:)
    complex(8),allocatable:: w_osmbuf(:)
    real(8),allocatable:: w_ofrbuf(:)
    nproc  = mpipid(0)
    procid = mpipid(1)
    if (nevec <= 0) return
    call tcn('rsibl')
    nlmto = ndimh-napw
    alat=lat_alat
    plat = lat_plat
    qlat = lat_qlat
    gmax=lat_gmax
    ngabc=lat_nabc
    ng=lat_ng
    vol=lat_vol
    ! ... Setup for q-dependent gv ... also makes kv, gv+q and iv    !     NB: gv generated by gvlst2 has q already added to it!
!    call tcn('gvlst2')
    call pshpr(iprint()-30)
    call gvlst2(alat,plat,q,n1,n2,n3,0d0,gmax,0,000,0,ng,xx,xx,xx)
    allocate(rv_a_ogv(ng,3))
    allocate(w_oiv(ng*3))
    allocate(iv_a_okv(ng*3))
    call gvlst2(alat,plat, q, n1, n2, n3, 0d0,gmax,0,509, ng, ng, iv_a_okv, rv_a_ogv, w_oiv)
    call poppr
!    call tcx('gvlst2')
    if (napw > 0) then
       allocate(ivp(napw))
       call gvgvcomp(ng,w_oiv,napw,igapw,ivp)
    else
       allocate(ivp(1))
    endif
    ! --- Tables of energies, rsm, indices to them ---
    call tbhsi(nspec,nermx,net,etab,ipet,nrt,rtab,iprt,ltop)
    ! --- Allocate and occupy arrays for yl, energy factors, rsm factors ---
    nlmtop = (ltop+1)**2
    allocate(w_ogq(ng*3),w_oyl(ng*nlmtop),w_oylw(ng*nlmtop), w_og2(ng), w_ohe(ng*net), w_ohr(ng*nrt))
    ! ino H_L(G)= \frac{-4 pi}{e-G^2} {cal Y}_L(-iG) exp(gamma(e-G^2))
    ! ino hsibl1 calculaets he=1/(e-G^2) and hr=exp(-gamma G^2)
    ! ino the other parts are calculated in rsibl5.
    q0=0d0
    if (nlmto > 0) then
       call hsibl1 ( net , etab , nrt , rtab , ltop , alat , q0 , ng &
            , rv_a_ogv , w_ogq , w_og2 , w_oyl , w_ohe ,  w_ohr )
    endif
    deallocate(w_og2)
    nblk = 16
    !  --- Loop over eigenstates ---
    allocate(psi(ng,nspc,nblk),vpsi(ng,nspc,nblk),wk(ng,nspc,nblk))
    allocate(psir(k1,k2,k3),cosi(ng),sini(ng),wk2(ng))
    ivecini= 1
    ivecend= nevec
    do  ivec = ivecini,ivecend, nblk !blocked calculation for future
       nvec = min(nblk, nevec-ivec+1)
       call rsibl1(0,q,nbas,ng,w_ogq,w_oiv,n1,n2, &
            n3,qlat,cosi,sini,w_oyl,w_oylw,w_ohe,w_ohr,wk, &
            wk2,vol,iprt,ipet,etab,rtab,ndimh,nlmto,nspc, &
            ewgt,ivec,nvec,evec,w,psi,wff)
       ! ino    rsiblp adds PW(G) to psi
       !call rsiblp(ng,ndimh,nlmto,nspc,napw,ivp,nvec,dsqrt(vol), evec(1,1,ivec),psi)
       if(napw>0) psi(ivp(:),:,1:nvec) = psi(ivp(:),:,1:nvec) &
                  +      evec(nlmto+1:nlmto+napw,:,ivec:ivec+nvec-1)/vol**.5
       ! ino now psi= H(G) + PW(G)
       !   ... Add to real-space mesh, optionally make smpot*psi for forces
       ! ino rsibl2 executes FFT to get psi(r), which is F0
       ! ino and also calculates <psi|psi>(=F0F0) to get real space charge density.
       !call rsibl2 ( ng , nspc , nvec , psi , n1 , n2 , n3 , k1 , k2 &
       !     , k3 , iv_a_okv ,  ewgt ( ivec ) , lfrce , smpot ( &
       !     1 , 1 , 1 , isp ) , psir , smrho ( 1 , 1 , 1 , isp ) , vpsi  )
       rsibl2block: block
         integer:: i,ispc
         real(8):: wgt1
         !Confusing, but it works well probably. We need transparent logic. nspc is confusing.
         !   When nsp=nspc=2,isp=1 (=jsp=1 supplied from m_bandcal, coupled spin case),
         ! we have smrho(:,:,:,1:2).
         ! When nsp=2 nspc=1 and isp=1or2, we have smrho(:,:,:,isp).
         do  ispc = 1, nspc !nspc=2 for SO=1 and nsp=2 (isp=jsp=1 is supplied).
            do  i = 1, nvec
               call gvputf(ng,1,iv_a_okv,k1,k2,k3,psi(1,ispc,i),psir)
               call fftz3(psir,n1,n2,n3,k1,k2,k3,1,0,1)
               wgt1 = ewgt(ivec+i-1)
               smrho(:,:,:,isp+ispc-1)=smrho(:,:,:,isp+ispc-1)+ wgt1*dconjg(psir(:,:,:))*psir(:,:,:)
               if (lfrce /= 0) then
                  psir(:,:,:) = psir(:,:,:)*smpot(:,:,:,isp+ispc-1)
                  call fftz3(psir,n1,n2,n3,k1,k2,k3,1,0,-1)
                  call gvgetf(ng,1,iv_a_okv,k1,k2,k3,psir,vpsi(1,ispc,i))
               endif
            enddo
         enddo
       endblock rsibl2block
       !    --- Add to forces ---
       if (lfrce /= 0) then
          call rsibl1(1,q,nbas,ng,w_ogq,w_oiv,n1,n2, &
               n3,qlat,cosi,sini,w_oyl,w_oylw,w_ohe,w_ohr, &
               wk,wk2,vol,iprt,ipet,etab,rtab,ndimh,nlmto,nspc, &
               ewgt,ivec,nvec,evec,vpsi(:,:,1:nvec),psi(:,:,1:nvec),f)
       endif
    enddo
    deallocate(psi,vpsi,wk,psir,cosi,sini,wk2)
    if (allocated(rv_a_ogv)) deallocate(rv_a_ogv)
    if (allocated(iv_a_okv)) deallocate(iv_a_okv)
    deallocate(ivp)
    call tcx('rsibl')
  end subroutine rsibl

  subroutine rsibl1(mode,q,nbas,ng,gq,iv,n1,n2,n3, &
       qlat,cosgp,singp,yl,ylw,he,hr,psi0,wk2,vol,iprt,ipet,etab,rtab, &
       ndimh,nlmto,nspc,ewgt,ivec,nvec,evec,vpsi,psi,f)
    use m_uspecb,only:uspecb
    use m_struc_def 
    use m_orbl,only: Orblib,ktab,ltab,offl,norb
    use m_sugcut,only:ngcut
    use m_lattic,only:rv_a_opos
    use m_lmfinit,only: ispec
    !- Make wave function for a block of evecs, or add contr. to forces
    ! ----------------------------------------------------------------------
    !i Inputs
    !i   mode  :0 make wave function
    !i         :1 Add 2*Re( (v psi+) grad(psi) ) to f
    !i   q     :Bloch wave number
    !i   nbas  :size of basis
    !i   ng    :number of G-vectors
    !i   gq    :2*pi/alat * (q+G) for all G-vectors
    !i   iv    :g-vectors as integer multiples of qlat (suphs0)
    !i   n1..3 :size uniform mesh for smooth density and potential
    !i   qlat  :primitive reciprocal lattice vectors, in units of 2*pi/alat
    !i   cosgp :cos(phase) for each g-vector
    !i   singp :sin(phase) for each g-vector
    !i   yl    :spherical harmonics for ng vectors
    !i   ylw   :work array of same dimension as yl
    !i   he    :table of energy factors
    !i   hr    :table of smoothing radius factors
    !i   psi0  :work array (dim ng*2*nspc*nev): psi sans phase factors
    !i   wk2   :work array of dimension ng
    !i   vol   :cell volume
    !o   iprt  :index to which entry in rt a given orbital belongs
    !i   ipet  :index to which entry in etab a given orbital belongs
    !i   etab  :table of all inequivalent energies
    !i   rtab  :table of all inequivalent smoothing radii
    !i   ndimh :dimensions evec
    !i   nspc  :2 for coupled spins; otherwise 1
    !i   ewgt  :weights for each of the trial fermi levels
    !i   ivec  :first of current block of eigenvectors
    !i   nvec  :number of eigenstates to generate
    !i   evec  :eigenvectors
    !i   vspi  :potential * wave function, needed only for mode=1
    !o Outputs
    !o   psi   :wave function (mode=0); work area (mode=1)
    !o   f     :term added to forces (mode=1)
    !r Remarks
    !u Updates
    ! ----------------------------------------------------------------------
    implicit none
    real(8):: q(3)
    integer :: mode,nbas,ng,ndimh,nlmto,nspc,ivec,nvec,iv(ng,3), n1,n2,n3,n0,nkap0,i_copy_size
    parameter (n0=10,nkap0=3)
    integer :: iprt(n0,nkap0,*),ipet(n0,nkap0,*)
    double precision :: vol,yl(ng,*),ylw(ng,*),he(ng,*),hr(ng,*), &
         wk2(ng),cosgp(ng),singp(ng),etab(*),rtab(*),gq(ng,3),f(3,nbas),ewgt(nvec+ivec-1),qlat(3,3)
    complex(8):: psi0(:,:,:),psi(:,:,:), evec(ndimh,nspc,ivec),vpsi(:,:,:)  !psi(ng,nspc,nvec)
    integer:: blks(n0*nkap0),ntab(n0*nkap0),ncut(n0,nkap0),lh(nkap0),nkapi
    double precision :: e,rsm,eh(n0,nkap0),rsmh(n0,nkap0),f0(3)
    double precision :: xx(n0),wt,p(3)
    integer :: ib,is,io,jo,l2,kp,ie,ir,ioff,nlm1,nlm2,iq,kb,lt,i
    integer::ncutt
    complex(8):: phase(ng),img=(0d0,1d0)
    real(8),parameter:: pi = 4d0*datan(1d0),tpi = 2d0*pi
    psi=0d0 
    if(nlmto == 0) return
    do ib = 1, nbas
       is=ispec(ib) !ssite(ib)%spec
       p=rv_a_opos(:,ib) !ssite(ib)%pos
       ncut=ngcut(:,:,is)
!       call suphas(q,p,ng,iv,n1,n2,n3,qlat,cosgp,singp)
       phase = exp(-img*tpi*sum(p*q)) * exp(-img*tpi*matmul(p, matmul(qlat, transpose(iv))))
!       cosgp=dreal(cc)
!       singp=dimag(cc)
       call orblib(ib) !Return norb,ltab,ktab,offl
       call uspecb(is,rsmh,eh)
       call gtbsl1(1,norb,ltab,ktab,rsmh,eh,ntab,blks)
       psi0=0d0
       if(mode == 1) psi=0d0
       do io = 1, norb
          if (blks(io) == 0) cycle
          jo = ntab(io)
          l2 = ltab(io)
          lt = ltab(jo)
          kp = ktab(io)
          ie = ipet(l2+1,kp,is)
          ir = iprt(l2+1,kp,is)
          ioff = offl(io)
          nlm1 = l2**2+1
          nlm2 = nlm1 + blks(io)-1
          rsm = rtab(ir)
          e   = etab(ie)
          ncutt=ncut(lt+1,kp) ! takao Apr2009
!          call ncutcorrect(ncutt,1,gq,ng)
          ! kino rsibl5 calculates 4 pi exp(e gamma) and {cal Y}_L(-iG)
          !      and make psi0(G), which is H(-iG), with hr and he.
          ! ... Phase and other factors
          rsibl5block: block !- Add contribution to wave function from one block of orbitals
            real(8),parameter:: pi=3.1415926535897931d0
            real(8):: fac
            integer:: ilm,l,i
            fac = 4d0*pi*dexp(e*rsm*rsm*0.25d0)/vol
            ! ... Combine G-dependent energy, rsm and YL factors
            do  ilm = nlm1, nlm2 ! ... Make vector evec*phase
               l=ll(ilm)
               do i = 1, min(ng,ncutt)
                  psi0(i,1:nspc,1:nvec) = psi0(i,1:nspc,1:nvec)  + he(i,ie)*hr(i,ir)*yl(i,ilm) &
                  *(0d0,-1d0)**(l+2)*fac*evec(ilm-nlm1+ioff+1,1:nspc,ivec:ivec+nvec-1)
               enddo
            enddo
          endblock rsibl5block 
       enddo
       do  i = 1,ng  ! psi= exp(i G R_i) * psi0
          psi(i,:,:) = psi(i,:,:)+ psi0(i,:,:)*phase(i) !dcmplx(cosgp(i),singp(i)) 
       enddo
       if (mode == 1) then
          rsibl4block: block
          real(8):: xx(ng)
          do i = 1, nvec
             xx = sum(dimag(vpsi(:,:,i))*dreal(psi(:,:,i)) - dreal(vpsi(:,:,i))*dimag(psi(:,:,i)),dim=2)
             f0(:) = 2d0*vol*matmul(xx,gq(:,:))
             wt = ewgt(i+ivec-1)
             f(:,ib) = f(:,ib) + wt*f0(:)
             do  kb = 1, nbas !!             This shouldn't be necessary
                f(:,kb) = f(:,kb) - wt*f0(:)/nbas
             enddo
          enddo
          endblock rsibl4block
       endif
    enddo
  end subroutine rsibl1

  subroutine gvgvcomp(ng,igv,napw,igapw,ivp)
    !- Find pointer ivp that maps igv to igapw.
    ! ----------------------------------------------------------------------
    !i Inputs
    !i   ng    :number of G-vectors
    !o   igv   :list of reciprocal lattice vectors G (gvlist.f)
    !i   napw   :number of R.L.V for PW basis
    !i   igapw :reciprocal lattice vectors for PW basis.
    !o Outputs
    !o   ivp   :if ig = ivp(jg), igv(jg) and nvec(ig) are same vector
    !r Remarks
    !r  This routine should be be cleaned up and ivp
    !r  used by rest of program in place of igapw
    !u Updates
    !u   05 Jul 08 (T. Kotani) first created
    ! ----------------------------------------------------------------------
    implicit none
    integer :: ng,igv(ng,3),napw,igapw(3,napw),ivp(napw), jg,ig
    do  ig = 1, napw
       do  jg = 1, ng
          if(sum(abs(igv(jg,:)-igapw(:,ig))) == 0) then
             ivp(ig) = jg
             goto 333
          endif
       enddo
       call rx('gvgvcomp wrong 111! maybe enlarge GMAX or so')
333    continue
       if(sum(abs(igapw(:,ig)-igv(ivp(ig),:)))/=0) call rx('bug in gvgvcomp.  Cannot find ivp')
    enddo
  end subroutine gvgvcomp

  subroutine rsibl_ev(isp,q,iq,ndimh,nspc,napw,igapw,nevec,evec,k1,k2,k3, n_eiglist,eiglist)
    use m_struc_def
    use m_w_psir
    use m_lmfinit,only:  lat_alat,nspec,ispec,sspec=>v_sspec,nbas
    use m_lattic,only: lat_qlat
    use m_lattic,only: lat_vol
    use m_supot,only: lat_nabc
    use m_supot,only: lat_ng
    use m_supot,only: lat_gmax
    use m_lattic,only: lat_plat,rv_a_opos
    !- Add smooth part of output density into smrho and forces.
    ! ----------------------------------------------------------------------
    !i Inputs
    !i  x lfrce :if nonzero, accumulate contribution to force
    !i   nbas  :size of basis
    !i   ssite :struct for site-specific information; see routine usite
    !i     Elts read: spec pos
    !i     Stored:    *
    !i     Passed to: rsibl1
    !i   sspec :struct for species-specific information; see routine uspec
    !i     Elts read: ngcut
    !i     Stored:    *
    !i     Passed to: tbhsi rsibl1 uspecb
    !i   slat  :struct for lattice information; see routine ulat
    !i     Elts read: alat plat qlat gmax nabc ng ogv okv vol
    !i     Stored:    *
    !i     Passed to: *
    !i   lfrce :1 calculate contribution to forces
    !i   nbas  :size of basis
    !i   q     :Bloch vector
    !i   iq    :index to current k-point
    !i   ndimh :dimension of hamiltonian
    !i   nspc  :2 for coupled spins; otherwise 1
    !i   napw  :number of augmented PWs in basis
    !i   igapw :vector of APWs, in units of reciprocal lattice vectors
    !i   nevec :number of eigenvectors with nonzero weights
    !i   evec  :eigenvectors
    !i   ewgt  :eigenvector weights
    !i   k1..3 :dimensions smpot,smrho
    !i   smpot :smooth potential on uniform mesh, needed for forces
    !o Outputs
    !o   smrho :smooth density accumulated for this qp
    !o   f     :force contribution accumulated for this qp
    !r Remarks
    !m MPI
    !m   Parallelise over the eigenvector loop. The vector block size is
    !m   chosen (in the range 6-16, by dstrbp.f) so as to distribute the
    !m   work optimally across processes. Two work arrays of the size of
    !m   smrho are allocated from the heap as buffers. Only one will be
    !m   needed under MPI-2. See comments in hsibl.
    !b Bugs
    !b    replace call to gvgvcomp and pass ipv as input
    !b    The non-F90 version should work, but it is no longer tested
    !u Updates
    !u   29 Dec 08 Unsuccessful attempt to make work with openmp
    !u   05 Jul 08 (T. Kotani) output density for new PW part
    !u   10 Sep 06 Added MPI parallelization in the spin-coupled case
    !u   23 Dec 04 Extended to spin-coupled case
    !u   25 Aug 04 Adapted to extended local orbitals
    !u   10 Apr 02 Redimensionsed eh,rsmh to accomodate larger lmax
    !u   15 Feb 02 (ATP) Added MPI parallelization
    !u   27 Aug 01 Extended to local orbitals.
    !u   12 Oct 00 Use q-dependent list of G vectors
    !u    6 Jul 00 attempt to vectorize by grouping eigenvectors in blocks
    !u   17 Jun 00 Spin polarized
    !u   23 May 00 Adapted from nfp rsif_q.f
    ! ----------------------------------------------------------------------
    implicit none
    integer :: procid, master, nproc, mpipid
    integer :: isp,k1,k2,k3,ndimh,nevec,iq,nspc
    integer :: napw,igapw(3,napw)
    real(8):: q(3)
    double complex evec(ndimh,nspc,nevec)
    integer:: n_eiglist
    integer:: eiglist(n_eiglist)
    integer :: n0,nkap0,nermx,npmx,nblk,nlmto
    parameter (n0=10,nkap0=3,nermx=100,npmx=128)
    integer:: ngabc(3) , n1 , n2 , n3 , nrt , net , ng , &
         nglob , ltop , nlmtop , ogq , og2 , ohe , ohr , oyl , oylw , &
         oiv , iprint
    integer,allocatable :: iv_a_okv(:)
    real(8),allocatable :: rv_a_ogv(:)
    equivalence (n1,ngabc(1)),(n2,ngabc(2)),(n3,ngabc(3))
    integer :: iprt(n0,nkap0,nermx),ipet(n0,nkap0,nermx)
    double precision :: alat,qlat(3,3),plat(3,3),q0(3),gmax,xx
    double precision :: vol
    double precision :: etab(nermx),rtab(nermx)
    integer :: ivec,nvec
    integer,allocatable:: ivp(:)
    complex(8),allocatable::psi(:,:,:),psir(:,:,:),vpsi(:,:,:), wk(:,:,:)
    real(8),allocatable:: cosi(:),sini(:),wk2(:)
    integer:: ivecini,ivecend
    integer,allocatable:: w_oiv(:)
    real(8),allocatable:: w_ogq(:),w_oyl(:),w_oylw(:),w_og2(:),w_ohe(:),w_ohr(:)
    complex(8),allocatable:: w_osmbuf(:)
    real(8),allocatable:: w_ofrbuf(:)
    real(8),allocatable:: ewgt(:)
    complex(8),allocatable:: smrho(:,:,:,:), smpot(:,:,:,:)
    logical:: cmdopt
    real(8):: w(1)!dummy
    complex(8):: wpsidummy(1,1,1)
    nproc  = mpipid(0)
    procid = mpipid(1)
    if (nevec <= 0) return
    nlmto = ndimh-napw
    alat = lat_alat
    plat = lat_plat
    qlat = lat_qlat
    gmax = lat_gmax
    ngabc= lat_nabc
    ng = lat_ng
    vol= lat_vol
    ! ... Setup for q-dependent gv ... also makes kv, gv+q and iv
    !     NB: gv generated by gvlst2 has q already added to it!
    call pshpr(iprint()-30)
    call gvlst2(alat,plat,q,n1,n2,n3,0d0,gmax,0,500,0,ng,xx,xx,xx)
    allocate(rv_a_ogv(abs(ng*3)))
    rv_a_ogv(:)=0.0d0
    allocate(iv_a_okv(abs(ng*3)))
    iv_a_okv(:)=0
    allocate(smrho(k1,k2,k3,isp))
    allocate(smpot(k1,k2,3,isp))
    allocate(w_oiv(ng*3))
    call gvlst2(alat, plat, q, n1, n2, n3, 0d0,gmax,0,509, ng, ng, iv_a_okv, rv_a_ogv, w_oiv)
    call poppr
    !     For PW basis ... for now.
    if (napw > 0) then
       allocate(ivp(napw))
       call gvgvcomp(ng,w_oiv,napw,igapw,ivp)
    else
       allocate(ivp(1))
    endif
    ! --- Tables of energies, rsm, indices to them ---
    call tbhsi(nspec,nermx,net,etab,ipet,nrt,rtab,iprt,ltop)
    ! --- Allocate and occupy arrays for yl, energy factors, rsm factors ---
    nlmtop = (ltop+1)**2
    allocate(w_ogq(ng*3),w_oyl(ng*nlmtop),w_oylw(ng*nlmtop), w_og2(ng), w_ohe(ng*net), w_ohr(ng*nrt))
    ! ino H_L(G)= \frac{-4 pi}{e-G^2} {cal Y}_L(-iG) exp(gamma(e-G^2))
    ! ino hsibl1 calculaets he=1/(e-G^2) and hr=exp(-gamma G^2)
    ! ino the other parts are calculated in rsibl5.
    q0=0d0
    if(nlmto>0) then
       call hsibl1 ( net , etab , nrt , rtab , ltop , alat , q0 , ng &
            , rv_a_ogv , w_ogq , w_og2 , w_oyl , w_ohe ,  w_ohr )
    endif
    deallocate(w_og2)
    nblk = 16 !nblk is for futre parallelization
    !  --- Loop over eigenstates ---
    allocate(psi(ng,nspc,nblk),vpsi(ng,nspc,nblk),wk(ng,nspc,nblk))
    allocate(psir(k1,k2,k3),cosi(ng),sini(ng),wk2(ng))
    ivecini= 1
    ivecend= nevec
    allocate(ewgt(nevec))
    ewgt=1.0d0
    write(*,'(a,9i5)') 'ivecinic,ivecend,nblk,nspc,nspec=',ivecini,ivecend,nblk,nspc,nspec
    ivecloop: do  ivec = ivecini,ivecend, nblk
       nvec = min(nblk, nevec-ivec+1)
       !   ... Add together total smooth wavefunction
       ! ino   rsibl1 calculates fourier transformed smoothed Hankel, H(G), to psi
       call rsibl1(0,q,nbas,ng,w_ogq,w_oiv,n1,n2, &
            n3,qlat,cosi,sini,w_oyl,w_oylw,w_ohe,w_ohr,wk, &
            wk2,vol,iprt,ipet,etab,rtab,ndimh,nlmto,nspc, &
            ewgt,ivec,nvec,evec,wpsidummy,psi,w)
       ! ino    rsiblp adds PW(G) to psi
       !call rsiblp(ng,ndimh,nlmto,nspc,napw,ivp,nvec,dsqrt(vol), evec(1,1,ivec),psi)
       if(napw>0) psi(ivp(:),:,1:nvec) = psi(ivp(:),:,1:nvec) &
            +     evec(nlmto+1:nlmto+napw,:,ivec:ivec+nvec-1)/vol**.5 
       ! ino now psi= H(G) + PW(G)
       ! ino write psi=F0 part
       if ( .TRUE. ) then
          call w_psir(ng , nspc , nvec , psi , n1 , n2 , n3 , k1 , k2 &
               , k3 , iv_a_okv , isp, q ,iq,  n_eiglist,eiglist &
               , plat,alat, nbas, rv_a_opos, sspec(ispec(1:nbas))%z, psir )
       endif
    enddo ivecloop
    deallocate(psi,vpsi,wk,psir,cosi,sini,wk2)
    if (allocated(ewgt)) deallocate(ewgt)
    if (allocated(rv_a_ogv)) deallocate(rv_a_ogv)
    if (allocated(iv_a_okv)) deallocate(iv_a_okv)
    deallocate(ivp)
  end subroutine rsibl_ev
end module m_rsibl
