module m_makusq !Accumulate coefficients (u,s,z) in all augmentation spheres for evec(:,iq,isp)
  use m_ll,only:ll
  public makusq
  private
  contains
subroutine makusq(nsites,isite,nev,isp,iq,q,evec, auszall)!Accumulate coefficients (u,s,z) in all augmentation spheres for evec(:,iq,isp)
  use m_lmfinit,only: ispec,sspec=>v_sspec,nbas,nlmax,nsp,nspc,nkapii,lhh,rsma
  use m_suham,only: ndham=>ham_ndham
  use m_igv2x,only: ndimh 
  use m_uspecb,only:uspecb
  use m_orbl,only: Orblib,ktab,ltab,offl,norb,blks
  use m_bstrux,only: bstrux_set,bstr
  implicit none
  intent(in)::    nsites,isite,nev,isp,iq,q,evec
  intent(out)::                                   auszall
  integer,parameter:: n0=10,nkap0=3
  integer:: isp,iq,nev,nsites,isite(nsites),nkape,ib,nkapi,is,nr,kmax,lmxa,lmxl,lmxh,i,nlma
  real(8):: q(3),eh(n0,nkap0),rsmh(n0,nkap0),a,rmt
  complex(8):: evec(ndimh,nspc,nev)
  complex(8),target:: auszall(nlmax,ndham*nspc,3,nsp,nsites,iq)
  auszall=0d0
  call tcn ('makusq')
  do  i = 1, nsites
     if (nsites == nbas) ib = i
     if (nsites /= nbas) ib = isite(i)
     is = ispec(ib)
     lmxa=sspec(is)%lmxa
     lmxl=sspec(is)%lmxl
     kmax=sspec(is)%kmxt
     nr=  sspec(is)%nr
     lmxh=sspec(is)%lmxb
     rmt =sspec(is)%rmt
     a =  sspec(is)%a
     if (lmxa == -1) cycle
     call uspecb(is,rsmh,eh)
     nkapi= nkapii(is)
     nkape= nkapii(is)
     nlma = (lmxa+1)**2
     radset: block !   --- Set up all radial head and tail functions, and their BC's
       real(8):: rofi_rv(nr), &
            fh_rv(nr*(lmxh+1)*nkapi),   xh_rv(nr*(lmxh+1)*nkapi),   vh(0:lmxh,nkapi), dh(0:lmxh,nkapi), &
            fp_rv(nr*(lmxa+1)*(kmax+1)),xp_rv(nr*(lmxa+1)*(kmax+1)),vp(0:lmxa,0:kmax),dp(0:lmxa,0:kmax)
       call radmsh(rmt,a,nr,rofi_rv )
       call fradhd(nkapi,eh,rsmh,lhh(:,is),lmxh,nr,rofi_rv,fh_rv,xh_rv, vh,dh) !head part
       call fradpk(kmax,rsma(is),lmxa,nr,rofi_rv,fp_rv,xp_rv, vp,dp)           !tail part
       puqs11:block
         integer:: ivec,ksp,ispc, io1,l1,ik1,nlm11,nlm12,ilm1,i1,ilma,k,l
         complex(8) ::cPkl(0:kmax,nlma)
         complex(8),pointer:: ausz(:,:,:,:)
         logical:: s12
         ausz => auszall(:,:,:,:,i,iq) !the coefficient for the projection onto (u,s,gz) for this site
         call bstrux_set(ib,q) !Get bstr
         ispcloop: do ispc = 1, nspc ! nspc=1 for so=0,2, nspc=2 for so=1. See m_lmfinit.f90)
            ksp = max(ispc,isp)
            do  ivec = 1, nev ! Loop over eigenfunctions
               do  k = 0, kmax
                  do  ilma = 1, nlma
                     cPkL(k,ilma) =  sum(evec(1:ndimh,ispc,ivec)*bstr(1:ndimh,ilma,k))
                  enddo
               enddo
               !call rlocb1(ndimh, nlma, kmax, evec(1,ispc,ivec), bstr,cPkl)
               call orblib(ib) !Return norb,ltab,ktab,offl
               do  io1 = 1, norb ! Contribution from head part
                  l1  = ltab(io1)
                  ik1 = ktab(io1)
                  nlm11 = l1**2+1
                  nlm12 = nlm11 + blks(io1)-1
                  i1 = offl(io1)-nlm11+1 !  i1 = hamiltonian offset for first orbital in block
                  if (ik1 <= nkape) s12=.true.
                  if (ik1 >  nkape) s12=.false.
                  do  ilm1 = nlm11, nlm12
                     l = ll(ilm1)
                     if(s12)  ausz(ilm1,ivec,1:2,ksp)= ausz(ilm1,ivec,1:2,ksp) + [vh(l,ik1),dh(l,ik1)] * evec(ilm1+i1,ispc,ivec)
                     if(.not.s12) ausz(ilm1,ivec,3,ksp) = ausz(ilm1,ivec,3,ksp) + evec(ilm1+i1,ispc,ivec)
                  enddo
               enddo
               do ilma = 1, nlma ! Contribution from tail part
                  l = ll(ilma)
                  ausz(ilma,ivec,1:2,ksp)=ausz(ilma,ivec,1:2,ksp)+[sum(vp(l,:)*cPkL(:,ilma)),sum(dp(l,:)*cPkL(:,ilma))]
               enddo
            enddo
         enddo ispcloop
       endblock puqs11
     endblock radset
  enddo
  call tcx('makusq')
end subroutine makusq
end module m_makusq
! Roughly examined at 2023jan ---------------------------------------------
!i Inputs
!i   nbas  :number of basis atoms
!i   isite :sites at which to calculate coefficients; see nsites
!i   nlmax :1st dimension of ausz (maximum nlma over all sites)
!i   ndham :dimensions ausz
!i   ndimh :dimensions evec
!i   napw  :number of G vectors in PW basis (gvlst2.f)
!i   nev   :number of eigenvectors for which to accumulate ausz
!i   nsp   :2 for spin-polarized case, otherwise 1
!i   nspc  :2 for so=1; otherwise 1
!i   isp   :spin channel
!i   iq    :qp index
!i   q     :Bloch vector
!i   evec  :eigenvectors for this q
!o Outputs
!o   ausz   :val,slo of w.f. at MT sphere surface added to ausz; see Remarks
!
!l Local variables
!l   ispc  :the current spin index in the coupled spins case.
!l         :Some quantities have no separate address space for each
!l         :spin in the indepedent-spins case (evec,evl,ewgt) but do
!l         :in the coupled-spins case.  A separate loop ispc=1..nspc
!l         :must be added for the latter case
!l         :ispc is the appropriate index for objects which distinguish
!l         :spins in the spin-coupled case only
!l   isp   :isp  is the appropriate index for objects which distinguish spins in the spin-uncoupled case only
!l   ksp   :the current spin index in both independent and coupled cases.
!l         :ksp is appropriate spin index for quantities that have
!l         :separate address space for each spin in every case
!l         :(potential- and density-like objects).
!r Remarks
!r   Makes coefficients for projection of wave function onto
!r   augmented functions (u,s,gz) which is valid inside the MT spheres.
!r   u and s are linear combinations of and phi,phidot defined as:
!r   u has val=1, slo=0 at rmax, s has val=0, slo=1
!r
!r   For example, for EELS matrix elements <nk|r|core> we will need
!r    |nk> = \sum_L(au_nkL*u_l*Y_L + as_nkL*s_l*Y_L + az_nkl*gz_l*Y_L)
!r
!r   These are generated from the potential later (see vcdmel)
!r   makusq returns the au_nkL and as_nkL az_nkl at one spin and k-pt for
!r   each of the sites in the unit cell.
!r   If nsites/=nbas, coeffs are made just for the nsites sites listed in isite.
!   ! ----------------------------------------------------------------------
!   !     In noncollinear case, isp=1 always => need internal ispc=1..2
!   !     ksp is the current spin index in both cases:
!   !     ksp = isp  in the collinear case
!   !         = ispc in the noncollinear case
!   !     whereas ispc=1 for independent spins, and spin index when nspc=2
