module m_hsibl ! Interstitial matrix elements of smooth Bloch Hankels, smooth potential.
  use m_ll,only:ll
  public hsibl,hsibl1
  private
contains
  subroutine hsibl(k1,k2,k3,vsm,isp,q,ndimh,napw,igapw, h)
    use m_lmfinit,only: alat=>lat_alat,nspec,nbas,ispec
    use m_lattic,only: qlat=>lat_qlat,vol=>lat_vol,rv_a_opos
    use m_supot,only: ngabc=>lat_nabc, lat_ng, gmax=>lat_gmax
    use m_lattic,only:  plat=>lat_plat
    use m_uspecb,only:uspecb
    use m_orbl,only: Orblib1,Orblib2,ktab1,ltab1,offl1,norb1,ktab2,ltab2,offl2,norb2
    use m_ftox
    use m_sugcut,only: ngcut
    use m_lmfinit,only: nkaph,nl,ndimx
    !i Inputs
    !i   k1,k2,k3 dimensions of vsm
    !i   vsm   :smoothed potential, real-space mesh
    !i   isp   : spin (1 or 2)
    !i   q     :Bloch vector
    !i   ndimh :dimension of hamiltonian
    !i   napw  :number of augmented PWs in basis
    !i   igapw :vector of APWs, in units of reciprocal lattice vectors
    !o Outputs
    !o   h     :interstitial matrix elements of vsm added to h
    !r Remarks
    !r  *How orbital is extracted and employed.  See Remarks in smhsbl.f
    implicit none
    integer :: k1,k2,k3,isp,ndimh,napw,igapw(3,napw)
    integer:: i , ib1 , ib2 , ie , ofh1 , ofh2 , ip , ir , is1 , &
         is2 , j , nlm1 , nlm2 , ik1 , ik2 , l1 , iorb1 , l2 , l2t &
         , iorb2 , jorb2 , osin1 , osin2 , oc1 , ocf1 , oc2 , ocf2 , ocos1 &
         , ocos2 , of , owk , ndim1 , ndim2 , nkap1 , nkap2
    integer,parameter:: n0=10,nkap0=3,nkape=2,nermx=100
    real(8):: q(3)
    real(8) :: eh1(n0,nkap0),rsmh1(n0,nkap0)
    real(8) :: eh2(n0,nkap0),rsmh2(n0,nkap0)
    real(8) :: xx(n0),p1(3),p2(3)
    complex(8):: h(ndimh,ndimh),vsm(k1,k2,k3,isp)
    integer :: npmx,nlmto
    integer:: ltop , n1 , n2 , n3 , net , nglob , nlmtop , nrt , iprint
    real(8) ,allocatable :: gg(:), g2(:), gvv(:),he(:,:), hr(:,:),yl(:,:)
    real(8) :: q0(3), etab(nermx),rtab(nermx)
    integer ,allocatable :: kv(:),iv(:,:),ncuti(:)
    complex(8) ,allocatable :: wk2_zv(:)
    integer :: iprt(n0,nkap0,nermx),ipet(n0,nkap0,nermx), ncut(n0,nkap0)!,lh(nkap0)
    integer:: blks1(n0*nkap0),ntab1(n0*nkap0)
    integer:: blks2(n0*nkap0),ntab2(n0*nkap0)
    integer:: xxxx(nkap0), ig1,i1,ig2,i2,igx(3),igx1,igx2,igx3,oiv1, iloop,ng
    integer:: ibini,ibend,nnrl,lmri,li,nnrlx,nnrli,ik,ib,ndim,ilm,offi
    complex(8) ,allocatable :: h_zv(:),phase(:,:)
    complex(8),allocatable,target:: w_oc1(:,:),w_ocf1(:),w_oc2(:,:),w_ocf2(:),ff(:)
    real(8),allocatable:: w_ocos1(:,:), w_osin1(:,:),w_ocos2(:,:), w_osin2(:,:),w_owk(:,:)
    real(8):: gmin=0d0,fac1
    complex(8):: img=(0d0,1d0)
    real(8),parameter:: pi=4d0*atan(1d0),tpi=2d0*pi
    call tcn('hsibl')
    ng=lat_ng
    n1=ngabc(1)
    n2=ngabc(2)
    n3=ngabc(3)
    nlmto = ndimh - napw
    ! --- <MTO|V|MTO> and < MTO|V|PW> parts of h ---
    if (nlmto > 0) then
       ! ... 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,gmin,gmax,0,500,0,ng,xx,xx,xx)!,xx)
       allocate(gvv(ng*3))
       allocate(kv(ng*3))
       allocate(iv(ng,3))
       call gvlst2(alat, plat, q, n1, n2, n3, gmin, gmax, 0, 509, ng, ng, kv, gvv, iv)!, iv)
       call poppr
!       call tcx('gvlst2')
       ! ... 
       call tbhsi(nspec,nermx,net,etab,ipet,nrt,rtab,iprt,ltop) !Tables of energies, rsm, indices to them
       !     ndimx = maximum hamiltonian dimension for any site (in m_lmfinit now)
       nlmtop = (ltop+1)**2
       allocate(gg(ng*3))
       allocate(yl(ng,nlmtop))
       allocate(g2(ng))
       allocate(he(ng,net))
       allocate(hr(ng,nrt),phase(ng,nbas))
       q0=0d0
       call hsibl1(net,etab,nrt,rtab,ltop,alat,q0,ng,gvv,  gg,g2,yl,he,hr)
       allocate( w_oc1( ng,ndimx), w_ocf1(ndimx))
       allocate( w_oc2( ng,ndimx), w_ocf2(ndimx), ff(k1*k2*k3))
       ibini=1
       ibend=nbas
       do  ib1 = ibini,ibend
          p1=rv_a_opos(:,ib1)
          phase(:,ib1)=exp(-img*tpi*sum(p1*q)) *exp(-img*tpi*matmul(p1, matmul(qlat, transpose(iv))))
       enddo
       ib1loop: do  iloop = ibini,ibend
          ib1=iloop
          ndim1 = 0
          is1=ispec(ib1)
          call orblib1(ib1) !norb1,ltab1,ktab1,offl1
          ofh1 = offl1(1)
          call uspecb(is1,rsmh1,eh1)
          call gtbsl1(1,norb1,ltab1,ktab1,rsmh1,eh1,ntab1,blks1) ![1,1,1,0,1]
          irob1loop: do  iorb1 = 1, norb1 !MTO orbital loop
             if (blks1(iorb1) == 0) cycle
             l1   = ltab1(iorb1)
             ik1  = ktab1(iorb1)
             nlm1 = l1**2+1
             nlm2 = nlm1 + blks1(iorb1)-1
             ie   = ipet(l1+1,ik1,is1) !energy index for etab(ie)
             ir   = iprt(l1+1,ik1,is1)
             fac1 = -4d0*pi*dexp(etab(ie)*rtab(ir)**2/4)/dsqrt(vol)
             offi = ndim1-nlm1+1
             do  ilm = nlm1, nlm2
                w_oc1(:,ilm+offi) = he(:,ie)*hr(:,ir)* yl(:,ilm)*phase(:,ib1) !PW coefficients for MTO(ilm+offi)
                w_ocf1(ilm+offi) = (-img)**ll(ilm) * fac1
             enddo
             ndim1 = ndim1 + max(blks1(iorb1),0)
          enddo irob1loop
          fvsm: block ! ... Multiply potential into wave functions for orbitals in ib1
            complex(8):: f(k1,k2,k3)
            do  i = 1, ndim1
               call gvputf(ng,1,kv,k1,k2,k3,w_oc1(1,i),f)
               call fftz3(f,n1,n2,n3,k1,k2,k3,1,0,1)
               f = f*vsm(:,:,:,isp)
               call fftz3(f,n1,n2,n3,k1,k2,k3,1,0,-1)
               call gvgetf(ng,1,kv,k1,k2,k3,f,w_oc1(1,i)) !w_oc1(ipw,i) = <MTO(i)|vsm for ith pw
            enddo
          endblock fvsm
          ib2loop: do 1010 ib2 = ib1, nbas ! Loop over second of (ib1,ib2) site pairs
             is2 =ispec(ib2)
             ncut=ngcut(:,:,is2)
             call orblib2(ib2) !norb2,ltab2,ktab2,offl2
             ofh2 = offl2(1)
             call uspecb(is2,rsmh2,eh2) 
             call gtbsl1(1,norb2,ltab2,ktab2,rsmh2,eh2,ntab2,blks2) ![1,1,1,0,1]
             ndim2 = 0
             do  iorb2 = 1, norb2
                if (blks2(iorb2) == 0) cycle
                nlm1 = l2**2+1
                nlm2 = nlm1 + blks2(iorb2)-1
                ndim2 = ndim2 + max(blks2(iorb2),0)
             enddo
             allocate(ncuti(ndim2))
             ndim2 = 0
             iorb2loop: do  iorb2 = 1, norb2
                if (blks2(iorb2) == 0) cycle
                jorb2 = ntab2(iorb2)
                l2t  = ltab2(jorb2)
                l2   = ltab2(iorb2)
                ik2  = ktab2(iorb2)
                nlm1 = l2**2+1
                nlm2 = nlm1 + blks2(iorb2)-1
                ie   = ipet(l2+1,ik2,is2)
                ir   = iprt(l2+1,ik2,is2)
                fac1 = -4d0*pi*dexp(etab(ie)*rtab(ir)**2/4)/dsqrt(vol)
                offi = ndim2-nlm1+1
                do  ilm = nlm1, nlm2
                   w_oc2(:,ilm+offi) = he(:,ie)*hr(:,ir)* yl(:,ilm)*phase(:,ib2) 
                   w_ocf2(ilm+offi) = (-img)**ll(ilm) * fac1
                enddo
                ncuti(ndim2+1:ndim2+nlm2-nlm1+1)=ncut(l2t+1,ik2)
                ndim2 = ndim2 + max(nlm2-nlm1+1,0)
             enddo iorb2loop
             !     ... Scalar products phi1*vsm*phi2 for all orbitals in (ib1,ib2)
             ! ncuti are only at Gamma point; thus symmetry can not be kept well for other k points.
             !call ncutcorrect ( ncuti , ndim2 , gvv , ng )
             hssblock: block
               integer::ncut,i2, io1,io2,ofw1,ofw2
               complex(8)::hss(ndim1,ndim2) 
               complex(8),pointer:: c1(:,:),c2(:),cf1(:),cf2(:)
               cf1=>w_ocf1(1:ndim1)
               cf2=>w_ocf2(1:ndim2)
               do  i2 = 1, ndim2
                  ncut = min(ng,ncuti(i2))
                  c1=>w_oc1(1:ncut,1:ndim1)
                  c2=>w_oc2(1:ncut,i2)
                  hss(:,i2)= dconjg(cf1)* matmul(dconjg(transpose(c1)),c2) *cf2(i2) ! = phi1*vsm*phi2
               enddo
               deallocate(ncuti)
               ofw1 = 0
               do  io1 = 1, norb1
                  if (blks1(io1) ==0) cycle
                  ofh1 = offl1(io1)
                  nlm1 = blks1(io1)
                  ofw2 = 0
                  do  io2 = 1, norb2
                     if (blks2(io2) == 0) cycle
                     ofh2 = offl2(io2)
                     nlm2 = blks2(io2)
                     do i1 = 1, nlm1
                        do  i2 = 1, nlm2
                           h(ofh1+i1,ofh2+i2) = h(ofh1+i1,ofh2+i2) + hss(ofw1+i1,ofw2+i2)
                        enddo
                     enddo
                     ofw2 = ofw2 + blks2(io2)
                  enddo
                  ofw1 = ofw1 + blks1(io1)
               enddo
             endblock hssblock
1010      enddo ib2loop
          hsmvsmpw: block !   ... Matrix elements <Hsm| Vsm |PW>
            integer:: i2x,ig,io1,ofw1
            do ig = 1, napw
               i2  = nlmto+ig
               i2x = ifindiv(igapw(1,ig),iv,ng) ! index matching igv,igapw
               ofw1 = 0
               do  io1 = 1, norb1
                  if (blks1(io1) == 0) cycle
                  ofh1 = offl1(io1)
                  nlm1 = blks1(io1)
                  do  i1 = 1, nlm1
                     h(ofh1+i1,i2)=h(ofh1+i1,i2) + dconjg( w_ocf1(ofw1+i1)*w_oc1(i2x, ofw1+i1) )
                  enddo
                  ofw1 = ofw1 + blks1(io1)
               enddo
            enddo
          endblock hsmvsmpw
       enddo ib1loop
       deallocate(hr, he, g2, yl, gg, iv, kv, gvv, w_oc1,w_ocf1, w_oc2, w_ocf2,ff) 
    endif
    ! --- <e^i qpG | V |e^i qpG'>/vol = V(G'-G) ---
    if (napw > 0) then
       call fftz3(vsm(1,1,1,isp),n1,n2,n3,k1,k2,k3,1,0,-1)
       do  ig1 = 1, napw
          i1 = ig1+nlmto
          do  ig2 = ig1, napw
             i2 = ig2+nlmto
             igx = igapw(:,ig1) - igapw(:,ig2)
             igx1 = mod(igx(1)+10*n1,n1)
             igx2 = mod(igx(2)+10*n2,n2)
             igx3 = mod(igx(3)+10*n3,n3)
             if (igx1<0 .OR. igx2<0 .OR. igx3<0) call rx('igx<0')
             h(i1,i2) = h(i1,i2) + vsm(igx1+1,igx2+1,igx3+1,isp)
          enddo
       enddo
       call fftz3(vsm(1,1,1,isp),n1,n2,n3,k1,k2,k3,1,0,1)
    endif
    do  i = 1, ndimh
       h(i:ndimh,i) = dconjg(h(i,i:ndimh)) ! ... Occupy second half of matrix
    enddo
    call tcx('hsibl')
  end subroutine hsibl
  subroutine hsibl1(net,etab,nrt,rtab,ltop,alat,q0,ng,gvv, gg,g2,yl,he,hr) !Make yl's, energy and rsm factors for list of G vectors
    use m_ropyln,only: ropyln
    !i Inputs
    !i   net   :size of table et
    !i   etab    :table of all inequivalent energies
    !i   nrt   :size of table rt
    !i   rtab    :table of all inequivalent smoothing radii
    !i   ltop  :largest l at any site
    !i   alat  :length scale of lattice and basis vectors, a.u.
    !i   q     :Bloch wave number
    !i   ng    :number of G-vectors
    !i   gvv    :list of reciprocal lattice vectors G (gvlist.f)
    !o Outputs
    !o   gg     :2*pi/alat * (q+gv) for all g-vectors
    !o   g2    :g**2
    !o   yl    :Y_L
    !o   he    :1/(etab-g2) for all inequivalent e's and g-vectors
    !o   hr    :dexp(-(rsm/2)**2*g2(i)) for all inequivalent rsm and g-vecs
    ! ----------------------------------------------------------------------
    implicit none
    integer:: ltop,net,ng,nrt,i,ie,ir
    real(8):: alat,q0(3),gvv(ng,3),gg(ng,3),yl(ng,1),he(ng,net),hr(ng,nrt),g2(ng),etab(net),rtab(nrt)
    real(8):: tpiba,gam
    real(8),parameter:: pi = 4d0*datan(1d0)
    do  i = 1, ng
       gg(i,:) = 2d0*pi/alat*(gvv(i,:)+q0(:)) 
    enddo
    call ropyln(ng,gg(1,1),gg(1,2),gg(1,3),ltop,ng, yl,g2)! ... Get yl and g2
    do  ie = 1, net ! ... Make the energy factors
       he(:,ie) = 1d0/(etab(ie)-g2(:))
    enddo
    do  ir = 1, nrt ! ... Make the rsm factors
       hr(:,ir) = dexp(-0.25d0*rtab(ir)**2 *g2(:))
    enddo
  end subroutine hsibl1
  integer function ifindiv(igapw,igv,ng) ! Find index in igv that corresponds to igapw
    !i   igapw :vector of APWs, in units of reciprocal lattice vectors
    !i   igv   :List of G vectors
    !i   ng    :number of group operations
    !o Outputs
    !o   ifindiv:index to igv that matches igapw
    implicit none
    integer :: ng,igapw(3),igv(ng,3)
    integer :: ig
    ifindiv=-999999
    do  ig = 1, ng
       if (igapw(1) /= igv(ig,1)) cycle
       if (igapw(2) /= igv(ig,2)) cycle
       if (igapw(3) /= igv(ig,3)) cycle
       ifindiv = ig
       return
    enddo
    call rx('ifindiv: igapw not found in igv')
  end function ifindiv
end module m_hsibl
