!
! Copyright (C) 2003-2007 Quantum-Espresso group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
#include "f_defs.h"
!
!------------------------------------------------------------------------
program pw2wannier90
  ! This is the interface to the Wannier90 code: see http://www.wannier.org
  !------------------------------------------------------------------------
  !
  USE io_global,  ONLY : stdout, ionode, ionode_id
  USE mp_global,  ONLY : mpime, kunit
  USE mp,         ONLY : mp_bcast
  USE cell_base,  ONLY : at, bg
  use lsda_mod,   ONLY : nspin, isk
  use klist,      ONLY : nkstot
  use ktetra,     ONLY : k1, k2, k3, nk1, nk2, nk3
  use io_files,   ONLY : nd_nmbr, prefix, tmp_dir
  use noncollin_module, ONLY : noncolin
  use wvfct,            ONLY : gamma_only
  use wannier
  !
  implicit none
  integer :: ios
  CHARACTER(LEN=4) :: spin_component
  CHARACTER(len=256) :: outdir

  ! these are in wannier module.....-> integer :: ispinw, ikstart, ikstop, iknum
  namelist / inputpp / outdir, prefix, spin_component, wan_mode, &
       seedname, write_unk, write_amn, write_mmn, write_spn, &
       wvfn_formatted, reduce_unk
  !
  call start_postproc (nd_nmbr)
  !
  ! Read input on i/o node and broadcast to the rest
  !
  ios = 0
  if(ionode) then
     !
     ! Check to see if we are reading from a file
     !
     call input_from_file()
     !
     !   set default values for variables in namelist
     !
     outdir = './'
     prefix = ' '
     seedname = 'wannier'
     spin_component = 'none'
     wan_mode = 'standalone'
     wvfn_formatted = .false.
     write_unk = .false.
     write_amn = .true.
     write_mmn = .true.
     write_spn = .false.
     reduce_unk= .false.
     !
     !     reading the namelist inputpp
     !
     read (5, inputpp, iostat=ios)
     !
     !     Check of namelist variables
     !
     tmp_dir = TRIM(outdir) 
     ! back to all nodes
  end if
  !
  call mp_bcast(ios,ionode_id)    
  if (ios /= 0) call errore( 'pw2wannier90', 'reading inputpp namelist', abs(ios))
  !
  ! broadcast input variable to all nodes
  !
  call mp_bcast(outdir,ionode_id)    
  call mp_bcast(tmp_dir,ionode_id)
  call mp_bcast(prefix,ionode_id)
  call mp_bcast(seedname,ionode_id)
  call mp_bcast(spin_component,ionode_id)
  call mp_bcast(wan_mode,ionode_id)
  call mp_bcast(wvfn_formatted,ionode_id)
  call mp_bcast(write_unk,ionode_id)
  call mp_bcast(write_amn,ionode_id)
  call mp_bcast(write_mmn,ionode_id)
  call mp_bcast(write_spn,ionode_id)
  call mp_bcast(reduce_unk,ionode_id)
  !
  !   Now allocate space for pwscf variables, read and check them.
  !
  logwann = .true.
  write(stdout,*)
  write(stdout,*) ' Reading nscf_save data'
  call read_file  
  write(stdout,*)
  !
  if (noncolin.and.gamma_only) call errore('pw2wannier90',&
       'Non-collinear and gamma_only not implemented',1)

  ! Here we should trap restarts from a different number of nodes.
  ! or attempts at kpoint distribution
  !
  SELECT CASE ( TRIM( spin_component ) )
  CASE ( 'up' )
     write(stdout,*) ' Spin CASE ( up )'
     ispinw  = 1
     ikstart = 1
     ikstop  = nkstot/2
     iknum   = nkstot/2
  CASE ( 'down' )
     write(stdout,*) ' Spin CASE ( down )'
     ispinw = 2
     ikstart = nkstot/2 + 1
     ikstop  = nkstot
     iknum   = nkstot/2
  CASE DEFAULT
     if(noncolin) then
        write(stdout,*) ' Spin CASE ( non-collinear )'
     else
        write(stdout,*) ' Spin CASE ( default = unpolarized )'
     end if
     ispinw = 0
     ikstart = 1
     ikstop  = nkstot
     iknum   = nkstot
  END SELECT
  !
  write(stdout,*)
  write(stdout,*) ' Wannier mode is: ',wan_mode
  write(stdout,*)
  !
  if(wan_mode.eq.'standalone') then
     !
     write(stdout,*) ' -----------------'
     write(stdout,*) ' *** Reading nnkp '
     write(stdout,*) ' -----------------'
     write(stdout,*)
     call read_nnkp
     write(stdout,*) ' Opening pp-files '
     call openfil_pp
     call ylm_expansion
     write(stdout,*)
     write(stdout,*)
     if(write_amn) then
        write(stdout,*) ' ---------------'
        write(stdout,*) ' *** Compute  A '
        write(stdout,*) ' ---------------'
        write(stdout,*)
        call compute_amn
        write(stdout,*)
     else
        write(stdout,*) ' -----------------------------'
        write(stdout,*) ' *** A matrix is not computed '
        write(stdout,*) ' -----------------------------'
        write(stdout,*)
     endif
     if(write_mmn) then
        write(stdout,*) ' ---------------'
        write(stdout,*) ' *** Compute  M '
        write(stdout,*) ' ---------------'
        write(stdout,*) 
        call compute_mmn
        write(stdout,*)
     else
        write(stdout,*) ' -----------------------------'
        write(stdout,*) ' *** M matrix is not computed '
        write(stdout,*) ' -----------------------------'
        write(stdout,*)
     endif
     write(stdout,*) ' ----------------'
     write(stdout,*) ' *** Write bands '
     write(stdout,*) ' ----------------'
     write(stdout,*)
     call write_band
     write(stdout,*)
     if(write_unk) then
        write(stdout,*) ' --------------------'
        write(stdout,*) ' *** Write plot info '
        write(stdout,*) ' --------------------'
        write(stdout,*)
        call write_plot
        write(stdout,*)
     else
        write(stdout,*) ' -----------------------------'
        write(stdout,*) ' *** Plot info is not printed '
        write(stdout,*) ' -----------------------------'
        write(stdout,*)
     endif
     write(stdout,*) ' ------------'
     write(stdout,*) ' *** Stop pp '
     write(stdout,*) ' ------------' 
     write(stdout,*)
     call stop_pp
     !
  endif
  !
  if(wan_mode.eq.'library') then
     !
!     seedname='wannier'
     write(stdout,*) ' Setting up...'
     call setup_nnkp
     write(stdout,*)
     write(stdout,*) ' Opening pp-files '
     call openfil_pp
     write(stdout,*)
     write(stdout,*) ' Ylm expansion'
     call ylm_expansion
     write(stdout,*)
     call compute_amn
     call compute_mmn
     call write_band
     call run_wannier
     if(write_unk) call write_plot
     call lib_dealloc
     call stop_pp
     !
  endif
  !
  if(wan_mode.eq.'wannier2sic') then
     !
     call read_nnkp
     call wan2sic
     !
  endif
  !
  stop
end program pw2wannier90
!
!-----------------------------------------------------------------------
subroutine lib_dealloc
  !-----------------------------------------------------------------------
  !
  use wannier

  implicit none

  deallocate(m_mat,u_mat,u_mat_opt,a_mat,eigval)

  return
end subroutine lib_dealloc
!
!-----------------------------------------------------------------------
subroutine setup_nnkp
  !-----------------------------------------------------------------------
  !
  use io_global, only : stdout, ionode, ionode_id
  use kinds,     only : DP
  use constants, only : eps6, tpi, bohr => BOHR_RADIUS_ANGS
  use cell_base, only : at, bg, alat
  use gvect,     only : g, gg
  use ions_base, only : nat, tau, ityp, atm
  use klist,     only : xk
  USE mp,        ONLY : mp_bcast, mp_sum
  use mp_global, ONLY : intra_pool_comm
  use wvfct,     only : nbnd,npwx,gamma_only
  use noncollin_module, ONLY : noncolin
  use wannier

  implicit none
  real(DP) :: g_(3), gg_
  integer  :: ik, ib, ig, iw, ia, indexb, type
  integer, allocatable :: ig_check(:,:)
  real(DP) :: xnorm, znorm, coseno
  integer  :: exclude_bands(nbnd)

  ! aam: translations between PW2Wannier90 and Wannier90
  ! pw2wannier90   <==>   Wannier90
  !    nbnd                num_bands_tot
  !    n_wannier           num_wann
  !    num_bands           num_bands
  !    nat                 num_atoms
  !    iknum               num_kpts
  !    rlatt               transpose(real_lattice)
  !    glatt               transpose(recip_lattice)
  !    kpt_latt            kpt_latt
  !    nnb                 nntot
  !    kpb                 nnlist
  !    g_kpb               nncell
  !    mp_grid             mp_grid
  !    center_w            proj_site
  !    l_w,mr_w,r_w        proj_l,proj_m,proj_radial
  !    xaxis,zaxis         proj_x,proj_z
  !    alpha_w             proj_zona
  !    exclude_bands       exclude_bands
  !    atcart              atoms_cart
  !    atsym               atom_symbols

  allocate( kpt_latt(3,iknum) )
  allocate( atcart(3,nat), atsym(nat) )
  allocate( kpb(iknum,num_nnmax), g_kpb(3,iknum,num_nnmax) )
  allocate( center_w(3,nbnd), alpha_w(nbnd), l_w(nbnd), &
       mr_w(nbnd), r_w(nbnd), zaxis(3,nbnd), xaxis(3,nbnd) )
  allocate( excluded_band(nbnd) )

  ! real lattice (Cartesians, Angstrom)
  rlatt(:,:) = transpose(at(:,:))*alat*bohr
  ! reciprocal lattice (Cartesians, Angstrom)
  glatt(:,:) = transpose(bg(:,:))*tpi/(alat*bohr)
  ! convert Cartesian k-points to crystallographic co-ordinates
  kpt_latt(:,:)=xk(:,:)
  CALL cryst_to_cart(iknum,kpt_latt,at,-1)
  ! atom co-ordinates in Cartesian co-ords and Angstrom units
  atcart(:,:) = tau(:,:)*bohr*alat
  ! atom symbols
  do ia=1,nat
     type=ityp(ia)
     atsym(ia)=atm(type)
  enddo

  ! MP grid dimensions
  call find_mp_grid()

  write(stdout,'("  - Number of atoms is (",i3,")")') nat 

#ifdef __WANLIB
  if (ionode) then
     call wannier_setup(seedname,mp_grid,iknum,rlatt, &               ! input
          glatt,kpt_latt,nbnd,nat,atsym,atcart,gamma_only,noncolin, & ! input
          nnb,kpb,g_kpb,num_bands,n_wannier,center_w, &               ! output
          l_w,mr_w,r_w,zaxis,xaxis,alpha_w,exclude_bands)             ! output
  endif
#endif

  call mp_bcast(nnb,ionode_id)
  call mp_bcast(kpb,ionode_id)
  call mp_bcast(g_kpb,ionode_id)
  call mp_bcast(num_bands,ionode_id)
  call mp_bcast(n_wannier,ionode_id)
  call mp_bcast(center_w,ionode_id)
  call mp_bcast(l_w,ionode_id)
  call mp_bcast(mr_w,ionode_id)
  call mp_bcast(r_w,ionode_id)
  call mp_bcast(zaxis,ionode_id)
  call mp_bcast(xaxis,ionode_id)
  call mp_bcast(alpha_w,ionode_id)
  call mp_bcast(exclude_bands,ionode_id)

  if(noncolin) then
     n_proj=n_wannier/2
  else
     n_proj=n_wannier
  end if

  allocate( gf(npwx,n_proj), csph(16,n_proj) ) 

  write(stdout,'("  - Number of wannier functions is (",i3,")")') n_wannier 

  excluded_band(1:nbnd)=.false.
  nexband=0
  band_loop: do ib=1,nbnd
     indexb=exclude_bands(ib)
     if (indexb>nbnd .or. indexb<0) then
        call errore('setup_nnkp',' wrong excluded band index ', 1)
     elseif (indexb.eq.0) then 
        exit band_loop
     else
        nexband=nexband+1
        excluded_band(indexb)=.true.
     endif
  enddo band_loop

  if ( (nbnd-nexband).ne.num_bands ) &
       call errore('setup_nnkp',' something wrong with num_bands',1)

  do iw=1,n_proj
     xnorm = sqrt(xaxis(1,iw)*xaxis(1,iw) + xaxis(2,iw)*xaxis(2,iw) + &
          xaxis(3,iw)*xaxis(3,iw))
     if (xnorm < eps6) call errore ('setup_nnkp',' |xaxis| < eps ',1)
     znorm = sqrt(zaxis(1,iw)*zaxis(1,iw) + zaxis(2,iw)*zaxis(2,iw) + &
          zaxis(3,iw)*zaxis(3,iw))
     if (znorm < eps6) call errore ('setup_nnkp',' |zaxis| < eps ',1)
     coseno = (xaxis(1,iw)*zaxis(1,iw) + xaxis(2,iw)*zaxis(2,iw) + &
          xaxis(3,iw)*zaxis(3,iw))/xnorm/znorm
     if (abs(coseno) > eps6) &
          call errore('setup_nnkp',' xaxis and zaxis are not orthogonal !',1)
     if (alpha_w(iw) < eps6) &
          call errore('setup_nnkp',' zona value must be positive', 1)
     ! convert wannier center in cartesian coordinates (in unit of alat)
     CALL cryst_to_cart( 1, center_w(:,iw), at, 1 )
  enddo
  write(stdout,*) ' - All guiding functions are given '

  nnbx=0
  nnb=max(nnbx,nnb)

  allocate( ig_(iknum,nnb), ig_check(iknum,nnb) )

  do ik=1, iknum
     do ib = 1, nnb
        g_(:) = REAL( g_kpb(:,ik,ib) )
        call trnvect (g_, at, bg, 1)
        gg_ = g_(1)*g_(1) + g_(2)*g_(2) + g_(3)*g_(3)
        ig_(ik,ib) = 0
        ig = 1
        do while  (gg(ig) <= gg_ + eps6) 
           if ( (abs(g(1,ig)-g_(1)) < eps6) .and.  &
                (abs(g(2,ig)-g_(2)) < eps6) .and.  &
                (abs(g(3,ig)-g_(3)) < eps6)  ) ig_(ik,ib) = ig
           ig= ig +1
        end do
     end do
  end do

  ig_check(:,:) = ig_(:,:)
  CALL mp_sum( ig_check, intra_pool_comm )
  do ik=1, iknum
     do ib = 1, nnb
        if (ig_check(ik,ib) ==0) &
          call errore('setup_nnkp', &
                      ' g_kpb vector is not in the list of Gs', 100*ik+ib )
     end do 
  end do 
  deallocate (ig_check)

  write(stdout,*) ' - All neighbours are found '
  write(stdout,*)
  
  return
end subroutine setup_nnkp
 !
 !-----------------------------------------------------------------------
subroutine run_wannier
  !-----------------------------------------------------------------------
  !
  use io_global, only : ionode, ionode_id
  use ions_base, only : nat
  use mp,        only : mp_bcast
  use wvfct,     only : gamma_only
  use wannier

  implicit none

  allocate(u_mat(n_wannier,n_wannier,iknum))
  allocate(u_mat_opt(num_bands,n_wannier,iknum))
  allocate(lwindow(num_bands,iknum))
  allocate(wann_centers(3,n_wannier))
  allocate(wann_spreads(n_wannier))

#ifdef __WANLIB
  if (ionode) then
     call wannier_run(seedname,mp_grid,iknum,rlatt, &                ! input
          glatt,kpt_latt,num_bands,n_wannier,nnb,nat, &              ! input
          atsym,atcart,gamma_only,m_mat,a_mat,eigval, &              ! input
          u_mat,u_mat_opt,lwindow,wann_centers,wann_spreads,spreads) ! output
  endif
#endif

  call mp_bcast(u_mat,ionode_id)
  call mp_bcast(u_mat_opt,ionode_id)
  call mp_bcast(lwindow,ionode_id)
  call mp_bcast(wann_centers,ionode_id)
  call mp_bcast(wann_spreads,ionode_id)
  call mp_bcast(spreads,ionode_id)

  return
end subroutine run_wannier
!-----------------------------------------------------------------------
!
subroutine find_mp_grid()
  !-----------------------------------------------------------------------
  !
  use io_global, only : stdout
  use kinds,     only: DP
  use wannier

  implicit none

  ! <<<local variables>>>
  integer  :: ik,ntemp,ii
  real(DP) :: min_k,temp(3,iknum),mpg1

  min_k=minval(kpt_latt(1,:))
  ii=0
  do ik=1,iknum
     if (kpt_latt(1,ik).eq.min_k) then
        ii=ii+1
        temp(:,ii)=kpt_latt(:,ik)
     endif
  enddo
  ntemp=ii

  min_k=minval(temp(2,1:ntemp))
  ii=0
  do ik=1,ntemp
     if (temp(2,ik).eq.min_k) then
        ii=ii+1
     endif
  enddo
  mp_grid(3)=ii

  min_k=minval(temp(3,1:ntemp))
  ii=0
  do ik=1,ntemp
     if (temp(3,ik).eq.min_k) then
        ii=ii+1
     endif
  enddo
  mp_grid(2)=ii

  if ( (mp_grid(2).eq.0) .or. (mp_grid(3).eq.0) ) &
       call errore('find_mp_grid',' one or more mp_grid dimensions is zero', 1)  

  mpg1=iknum/(mp_grid(2)*mp_grid(3))

  mp_grid(1) = nint(mpg1)

  write(stdout,*)
  write(stdout,'(3(a,i3))') '  MP grid is ',mp_grid(1),' x',mp_grid(2),' x',mp_grid(3)

  if (real(mp_grid(1),kind=DP).ne.mpg1) &
       call errore('find_mp_grid',' determining mp_grid failed', 1)

  return
end subroutine find_mp_grid
!-----------------------------------------------------------------------
!
subroutine read_nnkp
  !-----------------------------------------------------------------------
  !
  USE io_global, ONLY : stdout, ionode, ionode_id
  use kinds,     only: DP
  use constants, only : eps6, tpi, bohr => BOHR_RADIUS_ANGS
  use cell_base, only : at, bg, alat
  use gvect,     only : g, gg
  use io_files,  only : find_free_unit
  use klist,     only : nkstot, xk
  USE mp,        ONLY : mp_bcast, mp_sum
  use mp_global, ONLY : intra_pool_comm
  use wvfct,     only : npwx, nbnd
  use noncollin_module, ONLY : noncolin
  use wannier

  implicit none

  real(DP) :: g_(3), gg_
  integer :: ik, ib, ig, ipol, iw, idum, indexb
  integer numk, i, j
  integer, allocatable :: ig_check(:,:)
  real(DP) :: xx(3), xnorm, znorm, coseno
  CHARACTER(LEN=80) :: line1, line2
  logical :: have_nnkp

  if (ionode) then  ! Read nnkp file on ionode only

     inquire(file=TRIM(seedname)//".nnkp",exist=have_nnkp)
     if(.not. have_nnkp) then
        write(stdout,*) ' Could not find the file '//TRIM(seedname)//'.nnkp'
        stop
     end if

     iun_nnkp = find_free_unit()
     open (unit=iun_nnkp, file=TRIM(seedname)//".nnkp",form='formatted')

  endif

  nnbx=0
  
  !   check the information from *.nnkp with the nscf_save data
  write(stdout,*) ' Checking info from wannier.nnkp file' 
  write(stdout,*)
  
  if (ionode) then   ! read from ionode only

     call scan_file_to('real_lattice')
     do j=1,3
        read(iun_nnkp,*) (rlatt(i,j),i=1,3)
        do i = 1,3
           rlatt(i,j) = rlatt(i,j)/(alat*bohr)
        enddo
     enddo
     do j=1,3
        do i=1,3
           if(abs(rlatt(i,j)-at(i,j)).gt.eps6) then
              write(stdout,*)  ' Something wrong! '
              write(stdout,*)  ' rlatt(i,j) =',rlatt(i,j),  ' at(i,j)=',at(i,j)
              stop
           endif
        enddo
     enddo
     write(stdout,*) ' - Real lattice is ok'

     call scan_file_to('recip_lattice')
     do j=1,3
        read(iun_nnkp,*) (glatt(i,j),i=1,3)
        do i = 1,3
           glatt(i,j) = (alat*bohr)*glatt(i,j)/tpi
        enddo
     enddo
     do j=1,3
        do i=1,3
           if(abs(glatt(i,j)-bg(i,j)).gt.eps6) then
              write(stdout,*)  ' Something wrong! '
              write(stdout,*)  ' glatt(i,j)=',glatt(i,j), ' bg(i,j)=',bg(i,j)
              stop
           endif
        enddo
     enddo
     write(stdout,*) ' - Reciprocal lattice is ok'

     call scan_file_to('kpoints')
     read(iun_nnkp,*) numk
     if(numk.ne.iknum) then
        write(stdout,*)  ' Something wrong! '
        write(stdout,*)  ' numk=',numk, ' iknum=',iknum
        stop
     endif
     do i=1,numk
        read(iun_nnkp,*) xx(1), xx(2), xx(3)
        CALL cryst_to_cart( 1, xx, bg, 1 )
        if(abs(xx(1)-xk(1,i)).gt.eps6.or. &
             abs(xx(2)-xk(2,i)).gt.eps6.or. &
             abs(xx(3)-xk(3,i)).gt.eps6) then
           write(stdout,*)  ' Something wrong! '
           write(stdout,*) ' k-point ',i,' is wrong'
           write(stdout,*) xx(1), xx(2), xx(3) 
           write(stdout,*) xk(1,i), xk(2,i), xk(3,i)
           stop
        endif
     enddo
     write(stdout,*) ' - K-points are ok'

  endif ! ionode

  ! Broadcast
  call mp_bcast(rlatt,ionode_id)
  call mp_bcast(glatt,ionode_id)
  
  if (ionode) then   ! read from ionode only
     call scan_file_to('projections')
     read(iun_nnkp,*) n_proj
  endif

  ! Broadcast
  call mp_bcast(n_proj,ionode_id)

  if(noncolin) then
     n_wannier=n_proj*2
  else
     n_wannier=n_proj
  end if



  allocate( center_w(3,n_proj), alpha_w(n_proj), gf(npwx,n_proj), &
       l_w(n_proj), mr_w(n_proj), r_w(n_proj), &
       zaxis(3,n_proj), xaxis(3,n_proj), csph(16,n_proj) )

  write(stdout,'("  - Number of wannier functions is ok (",i3,")")') n_wannier 

  if (ionode) then   ! read from ionode only
     do iw=1,n_proj
        read(iun_nnkp,*) (center_w(i,iw), i=1,3), l_w(iw), mr_w(iw), r_w(iw)
        read(iun_nnkp,*) (zaxis(i,iw),i=1,3),(xaxis(i,iw),i=1,3),alpha_w(iw)
        xnorm = sqrt(xaxis(1,iw)*xaxis(1,iw) + xaxis(2,iw)*xaxis(2,iw) + &
             xaxis(3,iw)*xaxis(3,iw))
        if (xnorm < eps6) call errore ('read_nnkp',' |xaxis| < eps ',1)
        znorm = sqrt(zaxis(1,iw)*zaxis(1,iw) + zaxis(2,iw)*zaxis(2,iw) + &
             zaxis(3,iw)*zaxis(3,iw))
        if (znorm < eps6) call errore ('read_nnkp',' |zaxis| < eps ',1)
        coseno = (xaxis(1,iw)*zaxis(1,iw) + xaxis(2,iw)*zaxis(2,iw) + &
             xaxis(3,iw)*zaxis(3,iw))/xnorm/znorm
        if (abs(coseno) > eps6) &
             call errore('read_nnkp',' xaxis and zaxis are not orthogonal !',1)
        if (alpha_w(iw) < eps6) &
             call errore('read_nnkp',' zona value must be positive', 1)
        ! convert wannier center in cartesian coordinates (in unit of alat)
        CALL cryst_to_cart( 1, center_w(:,iw), at, 1 )
     enddo
  endif

  write(stdout,*) ' - All guiding functions are given '

  ! Broadcast
  call mp_bcast(center_w,ionode_id)
  call mp_bcast(l_w,ionode_id)
  call mp_bcast(mr_w,ionode_id)
  call mp_bcast(r_w,ionode_id)
  call mp_bcast(zaxis,ionode_id)
  call mp_bcast(xaxis,ionode_id)
  call mp_bcast(alpha_w,ionode_id)

  !
  write(stdout,*)
  write(stdout,*) 'Projections:'
  do iw=1,n_proj
     write(stdout,'(3f12.6,3i3,f12.6)') &
          center_w(1:3,iw),l_w(iw),mr_w(iw),r_w(iw),alpha_w(iw)
  enddo

  if (ionode) then   ! read from ionode only
     call scan_file_to('nnkpts')
     read (iun_nnkp,*) nnb
  endif

  ! Broadcast
  call mp_bcast(nnb,ionode_id)
  !
  nnbx = max (nnbx, nnb )
  !
  allocate ( kpb(iknum,nnbx), g_kpb(3,iknum,nnbx),&
             ig_(iknum,nnbx), ig_check(iknum,nnbx) )

  !  read data about neighbours
  write(stdout,*)
  write(stdout,*) ' Reading data about k-point neighbours '
  write(stdout,*)
  
  if (ionode) then
     do ik=1, iknum
        do ib = 1, nnb
           read(iun_nnkp,*) idum, kpb(ik,ib), (g_kpb(ipol,ik,ib), ipol =1,3)
        enddo
     enddo
  endif

  ! Broadcast
  call mp_bcast(kpb,ionode_id)
  call mp_bcast(g_kpb,ionode_id)

  do ik=1, iknum
     do ib = 1, nnb
        g_(:) = REAL( g_kpb(:,ik,ib) )
        call trnvect (g_, at, bg, 1)
        gg_ = g_(1)*g_(1) + g_(2)*g_(2) + g_(3)*g_(3)
        ig_(ik,ib) = 0
        ig = 1
        do while  (gg(ig) <= gg_ + eps6) 
           if ( (abs(g(1,ig)-g_(1)) < eps6) .and.  &
                (abs(g(2,ig)-g_(2)) < eps6) .and.  &
                (abs(g(3,ig)-g_(3)) < eps6)  ) ig_(ik,ib) = ig
           ig= ig +1
        end do
     end do
  end do
  ig_check(:,:) = ig_(:,:)
  CALL mp_sum( ig_check, intra_pool_comm )
  do ik=1, iknum
     do ib = 1, nnb
        if (ig_check(ik,ib) ==0) &
          call errore('read_nnkp', &
                      ' g_kpb vector is not in the list of Gs', 100*ik+ib )
     end do
  end do
  deallocate (ig_check)

  write(stdout,*) ' All neighbours are found '
  write(stdout,*)

  allocate( excluded_band(nbnd) )

  if (ionode) then     ! read from ionode only
     call scan_file_to('exclude_bands')
     read (iun_nnkp,*) nexband
     excluded_band(1:nbnd)=.false.
     do i=1,nexband
        read(iun_nnkp,*) indexb
        if (indexb<1 .or. indexb>nbnd) &
             call errore('read_nnkp',' wrong excluded band index ', 1)
        excluded_band(indexb)=.true.
     enddo
  endif

  ! Broadcast
  call mp_bcast(nexband,ionode_id)
  call mp_bcast(excluded_band,ionode_id)

  if (ionode) close (iun_nnkp)   ! ionode only

  return
end subroutine read_nnkp
!
!-----------------------------------------------------------------------
subroutine scan_file_to (keyword)
   !-----------------------------------------------------------------------
   !
   use wannier, only :iun_nnkp
   USE io_global,  ONLY : stdout
   implicit none
   character(len=*) :: keyword
   character(len=80) :: line1, line2
!
! by uncommenting the following line the file scan restarts every time 
! from the beginning thus making the reading independent on the order 
! of data-blocks
!   rewind (iun_nnkp)
!
10 continue
   read(iun_nnkp,*,end=20) line1, line2
   if(line1.ne.'begin')  goto 10
   if(line2.ne.keyword) goto 10
   return
20 write (stdout,*) keyword," data-block missing "
   stop
end subroutine scan_file_to
!
!-----------------------------------------------------------------------
subroutine compute_mmn
   !-----------------------------------------------------------------------
   !
   USE io_global,  ONLY : stdout, ionode
   use kinds,           only: DP
   use wvfct,           only : nbnd, npw, npwx, igk, g2kin, gamma_only
   use wavefunctions_module, only : evc, psic, evc_nc, psic_nc
   use gsmooth,         only: nls, nlsm, nrxxs, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s
   use klist,           only : nkstot, xk
   use io_files,        only : nwordwfc, iunwfc
   use io_files,        only : find_free_unit
   use gvect,           only : g, ngm, ecutwfc, gstart
   use cell_base,       only : tpiba2, omega, alat, tpiba, at, bg
   USE ions_base,       only : nat, ntyp => nsp, ityp, tau
   use constants,       only : tpi
   use uspp,            only : nkb, vkb
   USE uspp_param,      ONLY : nh, tvanp, lmaxq
   use becmod,          only : becp, rbecp
   USE noncollin_module,ONLY : noncolin, npol
   use wannier

   implicit none

   integer :: mmn_tot, ik, ikp, ipol, ib, npwq, i, m, n
   integer :: ikb, jkb, ih, jh, na, nt, ijkb0, ind, nbt
   integer :: ikevc, ikpevcq
   complex(DP), allocatable :: phase(:), aux(:), aux2(:), evcq(:,:), &
                               becp2(:,:), Mkb(:,:),evcq_nc(:,:,:),aux_nc(:,:)
   real(DP), allocatable    :: rbecp2(:,:)
   complex(DP), allocatable :: qb(:,:,:,:), qgm(:)
   real(DP), allocatable    :: qg(:), ylm(:,:), dxk(:,:)
   integer, allocatable     :: igkq(:)
   complex(DP)              :: mmn, ZDOTC, phase1
   real(DP)                 :: aa, arg, g_(3)
   character (len=9)        :: cdate,ctime
   character (len=60)       :: header
   logical                  :: any_uspp
   integer                  :: nn,inn,loop,loop2
   logical                  :: nn_found
   complex(DP)              :: spin

   any_uspp = ANY(tvanp(1:ntyp))

   if(any_uspp .and. noncolin) call errore('pw2wannier90',&
       'NCLS calculation not implimented with USP',1)

   allocate( phase(nrxxs), igkq(npwx) )
   if(noncolin) then
      allocate( aux_nc(npwx,npol), evcq_nc(npwx,npol,nbnd) )
   else
      allocate( aux(npwx), evcq(npwx,nbnd) )
   end if

   if (gamma_only) allocate(aux2(npwx))

   if (wan_mode.eq.'library') allocate(m_mat(num_bands,num_bands,nnb,iknum))
   
   if (wan_mode.eq.'standalone') then
      iun_mmn = find_free_unit()
      if (ionode) open (unit=iun_mmn, file=TRIM(seedname)//".mmn",form='formatted')
   endif

   if(write_spn.and.noncolin) then
      iun_spn = find_free_unit()
       if (ionode) open (unit=iun_spn, file=TRIM(seedname)//".spn",form='formatted')
    end if

   mmn_tot = 0
   do ik=1,iknum
      mmn_tot = mmn_tot + nnb * nbnd * nbnd
   end do
   !
   !   USPP
   !
   !
   if(any_uspp) then
      CALL init_us_1
      if (gamma_only) then
         allocate ( rbecp(nkb,nbnd),rbecp2(nkb,nbnd))
      else
         allocate ( becp(nkb,nbnd),becp2(nkb,nbnd))
      end if
   end if
   !
   !     qb is  FT of Q(r) 
   !
!   nbt = 0
!   do ik=1, iknum
!      nbt = nbt + nnb
!   enddo
   nbt = nnb * iknum
   !
   allocate( qg(nbt) )
   allocate (dxk(3,nbt))
   !
   ind = 0
   do ik=1,iknum
      do ib=1,nnb
         ind = ind + 1
         ikp = kpb(ik,ib) 
         !
         g_(:) = REAL( g_kpb(:,ik,ib) )
         call trnvect (g_, at, bg, 1)
         dxk(:,ind) = xk(:,ikp) +g_(:) - xk(:,ik) 
         qg(ind) = dxk(1,ind)*dxk(1,ind)+dxk(2,ind)*dxk(2,ind)+dxk(3,ind)*dxk(3,ind)
      enddo
!      write (stdout,'(i3,12f8.4)')  ik, qg((ik-1)*nnb+1:ik*nnb)
   enddo 
   !
   !  USPP
   !
   if(any_uspp) then

      allocate( ylm(nbt,lmaxq*lmaxq), qgm(nbt) )
      allocate( qb (nkb, nkb, ntyp, nbt) )
      !
      call ylmr2 (lmaxq*lmaxq, nbt, dxk, qg, ylm)
      qg(:) = sqrt(qg(:)) * tpiba
      !
      do nt = 1, ntyp
         if (tvanp (nt) ) then 
            do ih = 1, nh (nt)
               do jh = 1, nh (nt)
                  CALL qvan2 (nbt, ih, jh, nt, qg, qgm, ylm)
                  qb (ih, jh, nt, 1:nbt) = omega * qgm(1:nbt)
               enddo 
            enddo 
         endif 
      enddo 
      !
      deallocate (qg, qgm, ylm )
      !
   end if
   write (stdout,*) "MMN"

   if (wan_mode.eq.'standalone') then
      CALL date_and_tim( cdate, ctime )
      header='Created on '//cdate//' at '//ctime
      if (ionode) then
         write (iun_mmn,*) header
         write (iun_mmn,*) nbnd-nexband, iknum, nnb
      endif
   endif
   if(write_spn.and.noncolin) then
      CALL date_and_tim( cdate, ctime )
      header='Created on '//cdate//' at '//ctime
      if (ionode) then
         write (iun_spn,*) header
         write (iun_spn,*) nbnd-nexband, iknum, nnb
      endif
   end if
   !
   allocate( Mkb(nbnd,nbnd) )
   !
   write(stdout,'(a,i8)') ' iknum = ',iknum

   ind = 0
   do ik=1,iknum
      write (stdout,'(i8)') ik
      ikevc = ik + ikstart - 1 
      if(noncolin) then
         call davcio (evc_nc, nwordwfc, iunwfc, ikevc, -1 )
      else
         call davcio (evc, nwordwfc, iunwfc, ikevc, -1 )
      end if
      call gk_sort (xk(1,ik), ngm, g, ecutwfc / tpiba2, npw, igk, g2kin)
      if(write_spn.and.noncolin) then
         do loop=1,nbnd
            do loop2=1,nbnd
               spin=ZDOTC (npw, evc_nc(1,1,loop2),1,evc_nc(1,1,loop),1) - &
                    ZDOTC (npw, evc_nc(1,2,loop2),1,evc_nc(1,2,loop),1) 
               call reduce(2,spin)
               if (ionode) write(iun_spn,'(3i7,2es22.12)') loop2,loop, ik,spin
            end do
         end do
      endif

      !
      !  USPP
      !
      if(any_uspp) then
         call init_us_2 (npw, igk, xk(1,ik), vkb)
         ! below we compute the product of beta functions with |psi> 
         if (gamma_only) then
            call ccalbec (nkb, npwx, npw, nbnd, rbecp, vkb, evc)
         else
            call ccalbec (nkb, npwx, npw, nbnd, becp, vkb, evc)
         end if 
      end if
      !
      !
      !do ib=1,nnb(ik)
      do ib=1,nnb
         ind = ind + 1
         ikp = kpb(ik,ib)
! read wfc at k+b
         ikpevcq = ikp + ikstart - 1
         if(noncolin) then
            call davcio (evcq_nc, nwordwfc, iunwfc, ikpevcq, -1 )
         else
            call davcio (evcq, nwordwfc, iunwfc, ikpevcq, -1 )
         end if
         call gk_sort (xk(1,ikp), ngm, g, ecutwfc / tpiba2, npwq, igkq, g2kin)
! compute the phase
         phase(:) = (0.d0,0.d0)
         if ( ig_(ik,ib)>0) phase( nls(ig_(ik,ib)) ) = (1.d0,0.d0)
         call cft3s (phase, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, +2)
         !
         !  USPP
         !
         if(any_uspp) then
            call init_us_2 (npwq, igkq, xk(1,ikp), vkb)
            ! below we compute the product of beta functions with |psi> 
            if (gamma_only) then
               call ccalbec (nkb, npwx, npwq, nbnd, rbecp2, vkb, evcq)
            else
               call ccalbec (nkb, npwx, npwq, nbnd, becp2, vkb, evcq)
            end if
         end if
         !
         !
         Mkb(:,:) = (0.0d0,0.0d0) 
         !
         if (any_uspp) then
            ijkb0 = 0
            do nt = 1, ntyp
               if ( tvanp(nt) ) then
                  do na = 1, nat
                     !
                     arg = DOT_PRODUCT( dxk(:,ind), tau(:,na) ) * tpi 
                     phase1 = CMPLX ( COS(arg), -SIN(arg) )
                     !
                     if ( ityp(na) == nt ) then
                        do jh = 1, nh(nt)
                           jkb = ijkb0 + jh
                           do ih = 1, nh(nt)
                              ikb = ijkb0 + ih
                              !
                              do m = 1,nbnd
                                 if (excluded_band(m)) cycle
                                 if (gamma_only) then
                                    do n=1,m ! Mkb(m,n) is symmetric in m and n for gamma_only case
                                       if (excluded_band(n)) cycle
                                       Mkb(m,n) = Mkb(m,n) + &
                                            phase1 * qb(ih,jh,nt,ind) * &
                                            rbecp(ikb,m)  * rbecp2(jkb,n)                                     
                                    enddo
                                 else
                                    do n=1,nbnd
                                       if (excluded_band(n)) cycle
                                       Mkb(m,n) = Mkb(m,n) + &
                                            phase1 * qb(ih,jh,nt,ind) * &
                                            CONJG( becp(ikb,m) ) * becp2(jkb,n) 
                                    enddo
                                 endif
                              enddo ! m
                           enddo !ih
                        enddo !jh
                        ijkb0 = ijkb0 + nh(nt)
                     endif  !ityp
                  enddo  !nat 
               else  !tvanp
                  do na = 1, nat
                     if ( ityp(na) == nt ) ijkb0 = ijkb0 + nh(nt)
                  enddo
               endif !tvanp
            enddo !ntyp
         end if ! any_uspp
         !
         !
! loops on bands
         !
         if (wan_mode.eq.'standalone') then
            if (ionode) write (iun_mmn,'(7i5)') ik, ikp, (g_kpb(ipol,ik,ib), ipol=1,3)
         endif
         !
         do m=1,nbnd
            if (excluded_band(m)) cycle
            !
            if(noncolin) then
               psic_nc(:,:) = (0.d0, 0.d0)
               do ipol=1,2!npol    
                  psic_nc(nls (igk (1:npw) ),ipol ) = evc_nc(1:npw,ipol, m)
                  call cft3s (psic_nc(1,ipol), nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, +2)
                  psic_nc(1:nrxxs,ipol) = psic_nc(1:nrxxs,ipol) * phase(1:nrxxs)
                  call cft3s (psic_nc(1,ipol), nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, -2)
                  aux_nc(1:npwq,ipol) = psic_nc(nls (igkq(1:npwq) ),ipol )
               end do
            else
               psic(:) = (0.d0, 0.d0)
               psic(nls (igk (1:npw) ) ) = evc (1:npw, m)
               if(gamma_only) psic(nlsm(igk (1:npw) ) ) = conjg(evc (1:npw, m))
               call cft3s (psic, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, +2)
               psic(1:nrxxs) = psic(1:nrxxs) * phase(1:nrxxs)
               call cft3s (psic, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, -2)
               aux(1:npwq)  = psic(nls (igkq(1:npwq) ) )
            end if
            if(gamma_only) then
               if (gstart==2) psic(nlsm(1)) = (0.d0,0.d0)
               aux2(1:npwq) = conjg(psic(nlsm(igkq(1:npwq) ) ) )
            end if
            aa = 0.d0
            !
            !  Mkb(m,n) = Mkb(m,n) + \sum_{ijI} qb_{ij}^I * e^-i(b*tau_I)
            !             <psi_m,k1| beta_i,k1 > < beta_j,k2 | psi_n,k2 > 
            !
            if (gamma_only) then
               do n=1,m ! Mkb(m,n) is symmetric in m and n for gamma_only case
                  if (excluded_band(n)) cycle
                  mmn = ZDOTC (npwq, aux,1,evcq(1,n),1) &
                       + CONJG(ZDOTC(npwq,aux2,1,evcq(1,n),1))
                  call reduce(2,mmn)
                  Mkb(m,n) = mmn + Mkb(m,n)
                  if (m.ne.n) Mkb(n,m) = Mkb(m,n) ! fill other half of matrix by symmetry
                  aa = aa + abs(mmn)**2
               enddo
            elseif(noncolin) then
               do n=1,nbnd
                  if (excluded_band(n)) cycle
                  mmn=(0.d0, 0.d0)
                  do ipol=1,2
                     mmn = mmn+ZDOTC (npwq, aux_nc(1,ipol),1,evcq_nc(1,ipol,n),1)
                  end do
                  call reduce(2,mmn)
                  Mkb(m,n) = mmn + Mkb(m,n)
                  aa = aa + abs(mmn)**2
               enddo
            else
               do n=1,nbnd
                  if (excluded_band(n)) cycle
                  mmn = ZDOTC (npwq, aux,1,evcq(1,n),1)
                  call reduce(2,mmn)
                  Mkb(m,n) = mmn + Mkb(m,n)
                  aa = aa + abs(mmn)**2
               enddo
            endif
         end do   ! m


         do n=1,nbnd
            if (excluded_band(n)) cycle
            do m=1,nbnd
               if (excluded_band(m)) cycle
               if (wan_mode.eq.'standalone') then
                  if (ionode) write (iun_mmn,'(2f18.12)') Mkb(m,n)
               elseif (wan_mode.eq.'library') then
                  m_mat(m,n,ib,ik)=Mkb(m,n)
               else
                  call errore('compute_mmn',' value of wan_mode not recognised',1)
               endif
            enddo
         enddo
            
      end do !ib
   end do  !ik

   if (ionode .and. wan_mode.eq.'standalone') close (iun_mmn)
   if (ionode .and. write_spn .and. noncolin) close (iun_spn)
! 
   if (gamma_only) deallocate(aux2)
   deallocate (Mkb, dxk, phase, igkq)
   if(noncolin) then
      deallocate(aux_nc,evcq_nc)
   else
      deallocate(aux,evcq)
   end if

   if(any_uspp) then
      deallocate (  qb)
      if (gamma_only) then
          deallocate (rbecp,rbecp2)
       else
          deallocate(becp,becp2)
       end if
    end if
!
   write(stdout,*)
   write(stdout,*) ' MMN calculated'

   return
end subroutine compute_mmn
!
!-----------------------------------------------------------------------
subroutine compute_amn
   !-----------------------------------------------------------------------
   !
   USE io_global,  ONLY : stdout, ionode
   use kinds,           only : DP
   use klist,           only : nkstot, xk
   use wvfct,           only : nbnd, npw, npwx, igk, g2kin, gamma_only
   use wavefunctions_module, only : evc,evc_nc
   use io_files,        only : nwordwfc, iunwfc
   use io_files,        only : find_free_unit
   use gvect,           only : g, ngm, ecutwfc, gstart
   use cell_base,       only : tpiba2
   use uspp,            only : nkb, vkb
   use becmod,          only : becp, rbecp
   use wannier
   USE ions_base,       only : nat, ntyp => nsp, ityp, tau
   USE uspp_param,      ONLY : tvanp
   USE noncollin_module,ONLY : noncolin,npol

   implicit none

   complex(DP) :: amn, ZDOTC
   real(DP):: DDOT
   complex(DP), allocatable :: sgf(:,:)
   integer :: ik, ibnd, ibnd1, iw,i, ikevc, nt,ipol
   character (len=9)  :: cdate,ctime
   character (len=60) :: header
   logical            :: any_uspp

   !nocolin: we have half as many projections g(r) defined as wannier
   !         functions. We project onto (1,0) (ie up spin) and then onto
   !         (0,1) to obtain num_wann projections. jry


   !call read_gf_definition.....>   this is done at the beging

   any_uspp =ANY (tvanp(1:ntyp)) 

   if(any_uspp .and. noncolin) call errore('pw2wannier90',&
       'NCLS calculation not implimented with USP',1)

   if (wan_mode.eq.'library') allocate(a_mat(num_bands,n_wannier,iknum))

   if (wan_mode.eq.'standalone') then
      iun_amn = find_free_unit()
      if (ionode) open (unit=iun_amn, file=TRIM(seedname)//".amn",form='formatted')
   endif

   write (stdout,*) "AMN"

   if (wan_mode.eq.'standalone') then
      CALL date_and_tim( cdate, ctime ) 
      header='Created on '//cdate//' at '//ctime
      if (ionode) then
         write (iun_amn,*) header 
         write (iun_amn,*) nbnd-nexband,  iknum, n_wannier
      endif
   endif
   !
   allocate( sgf(npwx,n_proj))
   !
   if (any_uspp) then
      if(gamma_only) then
          allocate ( rbecp(nkb,n_wannier))
      else
          allocate ( becp(nkb,n_wannier))
      end if
      CALL init_us_1
   end if
   !
   write(stdout,'(a,i8)') ' iknum = ',iknum
   do ik=1,iknum
      write (stdout,'(i8)') ik
      ikevc = ik + ikstart - 1
      if(noncolin) then
         call davcio (evc_nc, nwordwfc, iunwfc, ikevc, -1 )
      else
         call davcio (evc, nwordwfc, iunwfc, ikevc, -1 )
      end if
      call gk_sort (xk(1,ik), ngm, g, ecutwfc / tpiba2, npw, igk, g2kin)
      call generate_guiding_functions(ik)   ! they are called gf(npw,n_proj)
      !
      !  USPP
      !
      if(any_uspp) then
         call init_us_2 (npw, igk, xk (1, ik), vkb)
         ! below we compute the product of beta functions with trial func.
         if (gamma_only) then
!!$            CALL pw_gemm( 'Y', nkb, n_wannier, npw, vkb, npwx,gf, npwx, rbecp, nkb )
            call ccalbec (nkb, npwx, npw, n_proj, rbecp, vkb, gf)
         else
            call ccalbec (nkb, npwx, npw, n_proj, becp, vkb, gf)
         end if
         ! and we use it for the product S|trial_func>
         call s_psi (npwx, npw, n_proj, gf, sgf)  
      else
         sgf(:,:) = gf(:,:)
      endif
      !
      if(noncolin) then
         ! we do the projection as g(r)*a(r) and g(r)*b(r)
         do ipol=1,npol
            do iw = 1,n_proj
               ibnd1 = 0 
               do ibnd = 1,nbnd
                  if (excluded_band(ibnd)) cycle
                  amn=(0.0_dp,0.0_dp)
                  amn = ZDOTC(npw,evc_nc(1,ipol,ibnd),1,sgf(1,iw),1)
                  call reduce(2,amn)
                  ibnd1=ibnd1+1
                  if (wan_mode.eq.'standalone') then
                     if (ionode) write(iun_amn,'(3i5,2f18.12)') ibnd1, iw+n_proj*(ipol-1), ik, amn
                  elseif (wan_mode.eq.'library') then
                     a_mat(ibnd1,iw+n_proj*(ipol-1),ik) = amn
                  else
                     call errore('compute_amn',' value of wan_mode not recognised',1)
                 endif
               end do
            end do
         end do
      else
         do iw = 1,n_proj
            ibnd1 = 0 
            do ibnd = 1,nbnd
               if (excluded_band(ibnd)) cycle
               if (gamma_only) then
                  amn = 2.0_dp*DDOT(2*npw,evc(1,ibnd),1,sgf(1,iw),1)
                  if (gstart==2) amn = amn - real(conjg(evc(1,ibnd))*sgf(1,iw))
               else
                  amn = ZDOTC(npw,evc(1,ibnd),1,sgf(1,iw),1) 
               end if
               call reduce(2,amn)
               ibnd1=ibnd1+1
               if (wan_mode.eq.'standalone') then
                  if (ionode) write(iun_amn,'(3i5,2f18.12)') ibnd1, iw, ik, amn
               elseif (wan_mode.eq.'library') then
                  a_mat(ibnd1,iw,ik) = amn
               else
                  call errore('compute_amn',' value of wan_mode not recognised',1)
               endif
            end do
         end do
      end if
   end do  ! k-points
   deallocate (sgf,csph)
   if(any_uspp) then 
     if (gamma_only) then 
        deallocate (rbecp)
     else
       deallocate (becp)
     end if
   end if
   !
   if (ionode .and. wan_mode.eq.'standalone') close (iun_amn)
   
   write(stdout,*)
   write(stdout,*) ' AMN calculated'

   return
end subroutine compute_amn
!
!
subroutine generate_guiding_functions(ik)
   !
   USE io_global,  ONLY : stdout
   use constants, only : pi, tpi, fpi, eps8
   use wvfct, only : npw, g2kin, igk, gamma_only
   use gvect, only : ig1, ig2, ig3, g, gstart
   use cell_base,  ONLY : tpiba2, omega, tpiba
   use wannier
   use klist,      only : xk 
   USE cell_base, ONLY : bg

   implicit none

   integer, parameter :: lmax=3, lmax2=(lmax+1)**2
   integer :: iw, ig, ik, bgtau(3), isph, l, mesh_r
   integer :: lmax_iw, lm, ipol, n1, n2, n3, nr1, nr2, nr3, iig
   real(DP) :: arg, anorm, fac, alpha_w2, yy, alfa, DDOT
   complex(DP) :: ZDOTC, kphase, lphase, gff, lph
   real(DP), allocatable :: gk(:,:), qg(:), ylm(:,:), radial(:,:)
   complex(DP), allocatable :: sk(:) 
   !
   allocate( gk(3,npw), qg(npw), ylm(npw,lmax2), sk(npw), radial(npw,0:lmax) )
   !
   do ig = 1, npw
      gk (1,ig) = xk(1, ik) + g(1, igk(ig) )
      gk (2,ig) = xk(2, ik) + g(2, igk(ig) )
      gk (3,ig) = xk(3, ik) + g(3, igk(ig) )
      qg(ig) = gk(1, ig)**2 +  gk(2, ig)**2 + gk(3, ig)**2
   enddo

   call ylmr2 (lmax2, npw, gk, qg, ylm)
   ! define qg as the norm of (k+g) in a.u.
   qg(:) = sqrt(qg(:)) * tpiba

   do iw = 1, n_proj
      !
      gf(:,iw) = (0.d0,0.d0)

      call radialpart(npw, qg, alpha_w(iw), r_w(iw), lmax, radial) 

      do lm = 1, lmax2
         if ( abs(csph(lm,iw)) < eps8 ) cycle
         l = int (sqrt( lm-1.d0))
         lphase = (0.d0,-1.d0)**l
         !
         do ig=1,npw
            gf(ig,iw) = gf(ig,iw) + csph(lm,iw) * ylm(ig,lm) * radial(ig,l) * lphase
         end do !ig
      end do ! lm
      do ig=1,npw
         iig = igk(ig)
         arg = ( gk(1,ig)*center_w(1,iw) + gk(2,ig)*center_w(2,iw) + &
                                           gk(3,ig)*center_w(3,iw) ) * tpi
         ! center_w are cartesian coordinates in units of alat 
         sk(ig) = CMPLX(cos(arg), -sin(arg) )
         gf(ig,iw) = gf(ig,iw) * sk(ig) 
      end do
      if (gamma_only) then
          anorm = 2.0_dp*DDOT(2*npw,gf(1,iw),1,gf(1,iw),1)
          if (gstart==2) anorm = anorm - abs(gf(1,iw))**2
      else
          anorm = REAL(ZDOTC(npw,gf(1,iw),1,gf(1,iw),1))
      end if
      call reduce(1,anorm)
!      write (stdout,*) ik, iw, anorm
      gf(:,iw) = gf(:,iw) / dsqrt(anorm)
   end do
   !
   deallocate ( gk, qg, ylm, sk, radial)
   return
end subroutine generate_guiding_functions

subroutine write_band
   USE io_global,  ONLY : stdout, ionode
   use wvfct, only : nbnd, et
   use klist, only : nkstot
   use constants, only: rytoev
   use io_files, only : find_free_unit
   use wannier

   implicit none

   integer ik, ibnd, ibnd1, ikevc

   if (wan_mode.eq.'standalone') then
      iun_band = find_free_unit()
      if (ionode) open (unit=iun_band, file=TRIM(seedname)//".eig",form='formatted')
   endif

   if (wan_mode.eq.'library') allocate(eigval(num_bands,iknum))

   do ik=ikstart,ikstop
      ikevc = ik - ikstart + 1
      ibnd1=0
      do ibnd=1,nbnd
         if (excluded_band(ibnd)) cycle
         ibnd1=ibnd1 + 1
         if (wan_mode.eq.'standalone') then
            if (ionode) write (iun_band,'(2i5,f18.12)') ibnd1, ikevc, et(ibnd,ik)*rytoev
         elseif (wan_mode.eq.'library') then
            eigval(ibnd1,ikevc) = et(ibnd,ik)*rytoev
         else
            call errore('write_band',' value of wan_mode not recognised',1)
         endif
      end do
   end do
   return
end subroutine write_band

subroutine write_plot
   USE io_global,  ONLY : stdout, ionode
   use wvfct, only : nbnd, npw, igk, g2kin, gamma_only
   use wavefunctions_module, only : evc, psic
   use io_files, only : find_free_unit, nwordwfc, iunwfc
   use wannier
   use gsmooth,         only : nls, nlsm, nrxxs, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s
   use klist,           only : nkstot, xk
   use gvect,           only : g, ngm, ecutwfc
   use cell_base,       only : tpiba2
   USE noncollin_module,ONLY : noncolin

   implicit none
   integer ik, ibnd, ibnd1, ikevc, i1, j, spin
   character*20 wfnname

   ! aam: 1/5/06: for writing smaller unk files 
   integer :: n1by2,n2by2,n3by2,i,k,idx,pos
   COMPLEX(DP),allocatable :: psic_small(:)   
   !-------------------------------------------!

#ifdef __PARA
   integer nxxs
   COMPLEX(DP),allocatable :: psic_all(:)
   nxxs = nrx1s * nrx2s * nrx3s
   allocate(psic_all(nxxs) )
#endif

   if(noncolin) call errore('pw2wannier90',&
       'write_unk not implemented with ncls',1)

   if (reduce_unk) then
      write(stdout,'(3(a,i5))') 'nr1s =',nr1s,'nr2s=',nr2s,'nr3s=',nr3s
      n1by2=(nr1s+1)/2;n2by2=(nr2s+1)/2;n3by2=(nr3s+1)/2
      write(stdout,'(3(a,i5))') 'n1by2=',n1by2,'n2by2=',n2by2,'n3by2=',n3by2
      allocate(psic_small(n1by2*n2by2*n3by2))   
   endif

   do ik=ikstart,ikstop

      ikevc = ik - ikstart + 1

      iun_plot = find_free_unit()
      !write(wfnname,200) p,spin
      spin=ispinw
      if(ispinw.eq.0) spin=1
      write(wfnname,200) ikevc, spin
200   format ('UNK',i5.5,'.',i1)

   if (ionode) then
      if(wvfn_formatted) then
         open (unit=iun_plot, file=wfnname,form='formatted')
         if (reduce_unk) then
            write(iun_plot,*)  n1by2,n2by2,n3by2, ikevc, nbnd-nexband
         else
            write(iun_plot,*)  nr1s,nr2s,nr3s, ikevc, nbnd-nexband
         endif
      else
         open (unit=iun_plot, file=wfnname,form='unformatted')
         if (reduce_unk) then
            write(iun_plot)  n1by2,n2by2,n3by2, ikevc, nbnd-nexband
         else
            write(iun_plot)  nr1s,nr2s,nr3s, ikevc, nbnd-nexband
         endif
      endif
   end if

      call davcio (evc, nwordwfc, iunwfc, ikevc, -1 )
      call gk_sort (xk(1,ik), ngm, g, ecutwfc / tpiba2, npw, igk, g2kin)

      ibnd1 = 0
      do ibnd=1,nbnd
         if (excluded_band(ibnd)) cycle
         ibnd1=ibnd1 + 1
         psic(:) = (0.d0, 0.d0)
         psic(nls (igk (1:npw) ) ) = evc (1:npw, ibnd)
         if (gamma_only)  psic(nlsm(igk (1:npw) ) ) = conjg(evc (1:npw, ibnd))
         call cft3s (psic, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, +2)
         if (reduce_unk) pos=0
#ifdef __PARA
         call cgather_smooth(psic,psic_all)
         if (reduce_unk) then
            do k=1,nr3s,2
               do j=1,nr2s,2
                  do i=1,nr1s,2
                     idx = (k-1)*nr2s*nr1s + (j-1)*nr1s + i
                     pos=pos+1
                     psic_small(pos) = psic_all(idx) 
                  enddo
               enddo
            enddo
         endif
      if (ionode) then
         if(wvfn_formatted) then
            if (reduce_unk) then
               write (iun_plot,'(2ES20.10)') (psic_small(j),j=1,n1by2*n2by2*n3by2)
            else
               write (iun_plot,*) (psic_all(j),j=1,nr1s*nr2s*nr3s)
            endif
         else
            if (reduce_unk) then
               write (iun_plot) (psic_small(j),j=1,n1by2*n2by2*n3by2)
            else
               write (iun_plot) (psic_all(j),j=1,nr1s*nr2s*nr3s)
            endif
         endif
      end if
#else
         if (reduce_unk) then
            do k=1,nr3s,2
               do j=1,nr2s,2
                  do i=1,nr1s,2
                     idx = (k-1)*nr2s*nr1s + (j-1)*nr1s + i
                     pos=pos+1
                     psic_small(pos) = psic(idx) 
                  enddo
               enddo
            enddo
         endif
         if(wvfn_formatted) then 
            if (reduce_unk) then
               write (iun_plot,'(2ES20.10)') (psic_small(j),j=1,n1by2*n2by2*n3by2)
            else
               write (iun_plot,*) (psic(j),j=1,nr1s*nr2s*nr3s)
            endif
         else
            if (reduce_unk) then
               write (iun_plot) (psic_small(j),j=1,n1by2*n2by2*n3by2)
            else
               write (iun_plot) (psic(j),j=1,nr1s*nr2s*nr3s)
            endif
         endif
#endif
      end do !ibnd

      if(ionode) close (unit=iun_plot)

   end do  !ik
   
   if (reduce_unk) deallocate(psic_small)   

#ifdef __PARA
   deallocate( psic_all )
#endif
   return
end subroutine write_plot

subroutine wan2sic 

  USE io_global,  ONLY : stdout
  USE kinds, only : DP
  use io_files, only : iunwfc, iunatsicwfc, nwordwfc, nwordwann
  USE cell_base, only : omega, tpiba2
  use gvect, only : g, ngm, ecutwfc
  use gsmooth, only: nls, nrxxs, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s
  use wavefunctions_module, only : evc, psic
  use wvfct, only : nbnd, npwx, npw, igk, g2kin
  use klist, only : nkstot, xk, wk
  use wannier

  integer :: i, j, nn, ik, ibnd, iw, ikevc 
  complex(DP), allocatable :: orbital(:,:), orb(:,:), u_matrix(:,:,:) 

  open (20, file = TRIM(seedname)//".dat" , form = 'formatted', status = 'unknown')
  write(stdout,*) ' wannier plot '

  allocate ( u_matrix( n_wannier, n_wannier, nkstot) )
  allocate ( orbital( npwx, n_wannier), orb( nrxxs, n_wannier))

  !
  do i = 1, n_wannier
     do j = 1, n_wannier
        do ik = 1, nkstot
           read (20, * ) u_matrix(i,j,ik)
           !do nn = 1, nnb(ik)
           do nn = 1, nnb
              read (20, * ) ! m_matrix (i,j,nkp,nn)
           enddo
        enddo  !nkp
     enddo !j
  enddo !i
  !
  orb(:,:) = (0.0d0,0.0d0)
  do ik=1,iknum
     ikevc = ik + ikstart - 1
     call davcio (evc, nwordwfc, iunwfc, ikevc, -1)
     call gk_sort (xk(1,ik), ngm, g, ecutwfc/tpiba2, npw, igk, g2kin)
     write(stdout,*) 'npw ',npw
     do iw=1,n_wannier
        do j=1,npw
           orbital(j,iw) = (0.0d0,0.0d0)
           do ibnd=1,n_wannier
              orbital(j,iw) = orbital(j,iw) + u_matrix(iw,ibnd,ik)*evc(j,ibnd)
              write(stdout,*) j, iw, ibnd, ik, orbital(j,iw), &
                              u_matrix(iw,ibnd,ik), evc(j,ibnd)
           enddo !ibnd
        end do  !j
     end do !wannier
     call davcio (orbital, nwordwann, iunatsicwfc, ikevc, +1)
  end do ! k-points

  deallocate ( u_matrix) 
  write(stdout,*) ' dealloc u '
  deallocate (  orbital)
  write(stdout,*) ' dealloc orbital '
  deallocate ( orb )
  write(stdout,*) ' dealloc orb '
  !
end subroutine wan2sic 

subroutine ylm_expansion 
   USE io_global,  ONLY : stdout
   use kinds, ONLY :  DP
   USE random_numbers,       ONLY : rndm
   use wannier
   implicit none
   ! local variables
   integer, parameter :: lmax2=16
   integer ::  lm, i, ir, iw, m
   real(DP) :: capel
   real(DP), allocatable :: r(:,:), rr(:), rp(:,:), ylm_w(:), ylm(:,:), mly(:,:)
   real(DP) :: u(3,3)

   allocate (r(3,lmax2), rp(3,lmax2), rr(lmax2), ylm_w(lmax2))
   allocate (ylm(lmax2,lmax2), mly(lmax2,lmax2) )

   ! generate a set of nr=lmax2 random vectors
   do ir=1,lmax2
      do i=1,3
         r(i,ir) = rndm() -0.5d0
      end do
   end do
   rr(:) = r(1,:)*r(1,:) + r(2,:)*r(2,:) + r(3,:)*r(3,:)
   !- compute ylm(ir,lm)
   call ylmr2(lmax2, lmax2, r, rr, ylm)
   !- store the inverse of ylm(ir,lm) in mly(lm,ir)
   call invmat(lmax2, ylm, mly, capel)
   !- check that r points are independent
   call check_inverse(lmax2, ylm, mly)

   do iw=1, n_proj

      !- define the u matrix that rotate the reference frame
      call set_u_matrix (xaxis(:,iw),zaxis(:,iw),u)
      !- find rotated r-vectors 
      rp(:,:) = matmul ( u(:,:) , r(:,:) )
      !- set ylm funtion according to wannier90 (l,mr) indexing in the rotaterd points
      call ylm_wannier(ylm_w,l_w(iw),mr_w(iw),rp,lmax2) 

      csph(:,iw) = matmul (mly(:,:), ylm_w(:))

!      write (stdout,*) 
!      write (stdout,'(2i4,2(2x,3f6.3))') l_w(iw), mr_w(iw), xaxis(:,iw), zaxis(:,iw)
!      write (stdout,'(16i6)')   (lm, lm=1,lmax2)
!      write (stdout,'(16f6.3)') (csph(lm,iw), lm=1,lmax2)

   end do
   deallocate (r, rp, rr, ylm_w, ylm, mly )

   return
end subroutine ylm_expansion

subroutine check_inverse(lmax2, ylm, mly)
   use kinds, ONLY :  DP
   use constants, ONLY :  eps8
   implicit none
   ! I/O variables
   integer :: lmax2
   real(DP) :: ylm(lmax2,lmax2), mly(lmax2,lmax2)
   ! local variables
   real(DP), allocatable :: uno(:,:)
   real(DP) :: capel
   integer :: lm
   !
   allocate (uno(lmax2,lmax2) )
   uno = matmul(mly, ylm)
   capel = 0.d0
   do lm = 1, lmax2
      uno(lm,lm) = uno(lm,lm) - 1.d0
   end do
   capel = capel + SUM ( abs(uno(1:lmax2,1:lmax2) ) )
!   write (stdout,*) "capel = ", capel
   if (capel > eps8) call errore('ylm_expansion', &
                    ' inversion failed: r(*,1:nr) are not all independent !!',1)
   deallocate (uno)
   return
end subroutine check_inverse
   
subroutine set_u_matrix(x,z,u)
   use kinds, ONLY :  DP
   use constants, ONLY : eps6
   implicit none
   ! I/O variables
   real(DP) :: x(3),z(3),u(3,3)
   ! local variables
   real(DP) :: xx, zz, y(3), coseno

   xx = sqrt(x(1)*x(1) + x(2)*x(2) + x(3)*x(3))
   if (xx < eps6) call errore ('set_u_matrix',' |xaxis| < eps ',1)
!   x(:) = x(:)/xx
   zz = sqrt(z(1)*z(1) + z(2)*z(2) + z(3)*z(3))
   if (zz < eps6) call errore ('set_u_matrix',' |zaxis| < eps ',1)
!   z(:) = z(:)/zz

   coseno = (x(1)*z(1) + x(2)*z(2) + x(3)*z(3))/xx/zz
   if (abs(coseno) > eps6) call errore('set_u_matrix',' xaxis and zaxis are not orthogonal !',1)

   y(1) = (z(2)*x(3) - x(2)*z(3))/xx/zz
   y(2) = (z(3)*x(1) - x(3)*z(1))/xx/zz
   y(3) = (z(1)*x(2) - x(1)*z(2))/xx/zz

   u(1,:) = x(:)/xx
   u(2,:) = y(:)
   u(3,:) = z(:)/zz

!   write (stdout,'(3f10.7)') u(:,:)

   return

end subroutine set_u_matrix

subroutine ylm_wannier(ylm,l,mr,r,nr) 
!
! this routine returns in ylm(r) the values at the nr points r(1:3,1:nr) 
! of the spherical harmonic identified  by indices (l,mr) 
! in table 3.1 of the wannierf90 specification.
! 
! No reference to the particular ylm ordering internal to quantum-espresso
! is assumed. 
!
! If ordering in wannier90 code is changed or extended this should be the 
! only place to be modified accordingly
!
   use kinds, ONLY :  DP
   use constants, ONLY : pi, fpi, eps8
   implicit none
! I/O variables
!
   integer :: l, mr, nr
   real(DP) :: ylm(nr), r(3,nr)
!
! local variables
!
   real(DP), external :: s, p_z,px,py, dz2, dxz, dyz, dx2my2, dxy, &
                        fz3, fxz2, fyz2, fzx2my2, fxyz, fxx2m3y2, fy3x2my2
   real(DP) :: rr, cost, phi
   integer :: ir
   real(DP) :: bs2, bs3, bs6, bs12
   bs2 = 1.d0/sqrt(2.d0)
   bs3=1.d0/sqrt(3.d0)
   bs6 = 1.d0/sqrt(6.d0)
   bs12 = 1.d0/sqrt(12.d0)
!
   if (l > 3 .OR. l < -5 ) call errore('ylm_wannier',' l out of range ', 1)
   if (l>=0) then
      if (mr < 1 .OR. mr > 2*l+1) call errore('ylm_wannier','mr out of range' ,1)
   else
      if (mr < 1 .OR. mr > abs(l)+1 ) call errore('ylm_wannier','mr out of range',1)
   end if

   do ir=1, nr
      rr = sqrt( r(1,ir)*r(1,ir) +  r(2,ir)*r(2,ir) + r(3,ir)*r(3,ir) )
      if (rr < eps8) call errore('ylm_wannier',' rr too small ',1)

      cost =  r(3,ir) / rr
      !
      !  beware the arc tan, it is defined modulo pi
      !
      if (r(1,ir) > eps8) then
         phi = atan( r(2,ir)/r(1,ir) )
      else if (r(1,ir) < -eps8 ) then
         phi = atan( r(2,ir)/r(1,ir) ) + pi
      else
         phi = sign( pi/2.d0,r(2,ir) )
      end if

    
      if (l==0) then   ! s orbital
                    ylm(ir) = s(cost,phi)  
      end if
      if (l==1) then   ! p orbitals
         if (mr==1) ylm(ir) = p_z(cost,phi) 
         if (mr==2) ylm(ir) = px(cost,phi)
         if (mr==3) ylm(ir) = py(cost,phi)
      end if
      if (l==2) then   ! d orbitals
         if (mr==1) ylm(ir) = dz2(cost,phi)
         if (mr==2) ylm(ir) = dxz(cost,phi)
         if (mr==3) ylm(ir) = dyz(cost,phi)
         if (mr==4) ylm(ir) = dx2my2(cost,phi)
         if (mr==5) ylm(ir) = dxy(cost,phi)
      endif
      if (l==3) then   ! f orbitals
         if (mr==1) ylm(ir) = fz3(cost,phi)
         if (mr==2) ylm(ir) = fxz2(cost,phi)
         if (mr==3) ylm(ir) = fyz2(cost,phi)
         if (mr==4) ylm(ir) = fzx2my2(cost,phi)
         if (mr==5) ylm(ir) = fxyz(cost,phi)
         if (mr==6) ylm(ir) = fxx2m3y2(cost,phi)
         if (mr==7) ylm(ir) = fy3x2my2(cost,phi)
      endif
      if (l==-1) then  !  sp hybrids
         if (mr==1) ylm(ir) = bs2 * ( s(cost,phi) + px(cost,phi) ) 
         if (mr==2) ylm(ir) = bs2 * ( s(cost,phi) - px(cost,phi) ) 
      end if
      if (l==-2) then  !  sp2 hybrids 
         if (mr==1) ylm(ir) = bs3*s(cost,phi)-bs6*px(cost,phi)+bs2*py(cost,phi)
         if (mr==2) ylm(ir) = bs3*s(cost,phi)-bs6*px(cost,phi)-bs2*py(cost,phi)
         if (mr==3) ylm(ir) = bs3*s(cost,phi) +2.d0*bs6*px(cost,phi) 
      end if
      if (l==-3) then  !  sp3 hybrids
         if (mr==1) ylm(ir) = 0.5d0*(s(cost,phi)+px(cost,phi)+py(cost,phi)+p_z(cost,phi))
         if (mr==2) ylm(ir) = 0.5d0*(s(cost,phi)+px(cost,phi)-py(cost,phi)-p_z(cost,phi))
         if (mr==3) ylm(ir) = 0.5d0*(s(cost,phi)-px(cost,phi)+py(cost,phi)-p_z(cost,phi))
         if (mr==4) ylm(ir) = 0.5d0*(s(cost,phi)-px(cost,phi)-py(cost,phi)+p_z(cost,phi))
      end if
      if (l==-4) then  !  sp3d hybrids
         if (mr==1) ylm(ir) = bs3*s(cost,phi)-bs6*px(cost,phi)+bs2*py(cost,phi)
         if (mr==2) ylm(ir) = bs3*s(cost,phi)-bs6*px(cost,phi)-bs2*py(cost,phi)
         if (mr==3) ylm(ir) = bs3*s(cost,phi) +2.d0*bs6*px(cost,phi) 
         if (mr==4) ylm(ir) = bs2*p_z(cost,phi)+bs2*dz2(cost,phi)
         if (mr==5) ylm(ir) =-bs2*p_z(cost,phi)+bs2*dz2(cost,phi)
      end if
      if (l==-5) then  ! sp3d2 hybrids
         if (mr==1) ylm(ir) = bs6*s(cost,phi)-bs2*px(cost,phi)-bs12*dz2(cost,phi)+.5d0*dx2my2(cost,phi)
         if (mr==2) ylm(ir) = bs6*s(cost,phi)+bs2*px(cost,phi)-bs12*dz2(cost,phi)+.5d0*dx2my2(cost,phi)
         if (mr==3) ylm(ir) = bs6*s(cost,phi)-bs2*py(cost,phi)-bs12*dz2(cost,phi)-.5d0*dx2my2(cost,phi)
         if (mr==4) ylm(ir) = bs6*s(cost,phi)+bs2*py(cost,phi)-bs12*dz2(cost,phi)-.5d0*dx2my2(cost,phi)
         if (mr==5) ylm(ir) = bs6*s(cost,phi)-bs2*p_z(cost,phi)+bs3*dz2(cost,phi)
         if (mr==6) ylm(ir) = bs6*s(cost,phi)+bs2*p_z(cost,phi)+bs3*dz2(cost,phi)
      end if

   end do

   return

end subroutine ylm_wannier

!======== l = 0 =====================================================================
function s(cost,phi)
   use kinds, ONLY :  DP
   use constants, ONLY : fpi
   implicit none
   real(DP) :: s, cost,phi
   s = 1.d0/ sqrt(fpi)
   return
end function s
!======== l = 1 =====================================================================
function p_z(cost,phi)
   use kinds, ONLY :  DP
   use constants, ONLY : fpi
   implicit none
   real(DP) ::p_z, cost,phi
   p_z =  sqrt(3.d0/fpi) * cost
   return
end function p_z
function px(cost,phi)
   use kinds, ONLY :  DP
   use constants, ONLY : fpi
   implicit none
   real(DP) ::px, cost, phi, sint
   sint = sqrt(abs(1.d0 - cost*cost))
   px =  sqrt(3.d0/fpi) * sint * cos(phi)
   return
end function px
function py(cost,phi)
   use kinds, ONLY :  DP
   use constants, ONLY : fpi
   implicit none
   real(DP) ::py, cost, phi, sint
   sint = sqrt(abs(1.d0 - cost*cost))
   py =  sqrt(3.d0/fpi) * sint * sin(phi)
   return
end function py
!======== l = 2 =====================================================================
function dz2(cost,phi)
   use kinds, ONLY :  DP
   use constants, ONLY : fpi
   implicit none
   real(DP) ::dz2, cost, phi
   dz2 =  sqrt(1.25d0/fpi) * (3.d0* cost*cost-1.d0)
   return
end function dz2
function dxz(cost,phi)
   use kinds, ONLY :  DP
   use constants, ONLY : fpi
   implicit none
   real(DP) ::dxz, cost, phi, sint
   sint = sqrt(abs(1.d0 - cost*cost))
   dxz =  sqrt(15.d0/fpi) * sint*cost * cos(phi)
   return
end function dxz
function dyz(cost,phi)
   use kinds, ONLY :  DP
   use constants, ONLY : fpi
   implicit none
   real(DP) ::dyz, cost, phi, sint
   sint = sqrt(abs(1.d0 - cost*cost))
   dyz =  sqrt(15.d0/fpi) * sint*cost * sin(phi)
   return
end function dyz
function dx2my2(cost,phi)
   use kinds, ONLY :  DP
   use constants, ONLY : fpi
   implicit none
   real(DP) ::dx2my2, cost, phi, sint
   sint = sqrt(abs(1.d0 - cost*cost))
   dx2my2 =  sqrt(3.75d0/fpi) * sint*sint * cos(2.d0*phi)
   return
end function dx2my2
function dxy(cost,phi)
   use kinds, ONLY :  DP
   use constants, ONLY : fpi
   implicit none
   real(DP) ::dxy, cost, phi, sint
   sint = sqrt(abs(1.d0 - cost*cost))
   dxy =  sqrt(3.75d0/fpi) * sint*sint * sin(2.d0*phi)
   return
end function dxy
!======== l = 3 =====================================================================
function fz3(cost,phi)
   use kinds, ONLY :  DP
   use constants, ONLY : pi
   implicit none
   real(DP) ::fz3, cost, phi
   fz3 =  0.25d0*sqrt(7.d0/pi) * ( 5.d0 * cost * cost - 3.d0 ) * cost
   return
end function fz3
function fxz2(cost,phi)
   use kinds, ONLY :  DP
   use constants, ONLY : pi
   implicit none
   real(DP) ::fxz2, cost, phi, sint
   sint = sqrt(abs(1.d0 - cost*cost))
   fxz2 =  0.25d0*sqrt(10.5d0/pi) * ( 5.d0 * cost * cost - 1.d0 ) * sint * cos(phi)
   return
end function fxz2
function fyz2(cost,phi)
   use kinds, ONLY :  DP
   use constants, ONLY : pi
   implicit none
   real(DP) ::fyz2, cost, phi, sint
   sint = sqrt(abs(1.d0 - cost*cost))
   fyz2 =  0.25d0*sqrt(10.5d0/pi) * ( 5.d0 * cost * cost - 1.d0 ) * sint * sin(phi)
   return
end function fyz2
function fzx2my2(cost,phi)
   use kinds, ONLY :  DP
   use constants, ONLY : pi
   implicit none
   real(DP) ::fzx2my2, cost, phi, sint
   sint = sqrt(abs(1.d0 - cost*cost))
   fzx2my2 =  0.25d0*sqrt(105d0/pi) * sint * sint * cost * cos(2.d0*phi)
   return
end function fzx2my2
function fxyz(cost,phi)
   use kinds, ONLY :  DP
   use constants, ONLY : pi
   implicit none
   real(DP) ::fxyz, cost, phi, sint
   sint = sqrt(abs(1.d0 - cost*cost))
   fxyz =  0.25d0*sqrt(105d0/pi) * sint * sint * cost * sin(2.d0*phi)
   return
end function fxyz
function fxx2m3y2(cost,phi)
   use kinds, ONLY :  DP
   use constants, ONLY : pi
   implicit none
   real(DP) ::fxx2m3y2, cost, phi, sint
   sint = sqrt(abs(1.d0 - cost*cost))
   fxx2m3y2 =  0.25d0*sqrt(17.5d0/pi) * sint * sint * sint * cos(3.d0*phi)
   return
end function fxx2m3y2
function fy3x2my2(cost,phi)
   use kinds, ONLY :  DP
   use constants, ONLY : pi
   implicit none
   real(DP) ::fy3x2my2, cost, phi, sint
   sint = sqrt(abs(1.d0 - cost*cost))
   fy3x2my2 =  0.25d0*sqrt(17.5d0/pi) * sint * sint * sint * sin(3.d0*phi)
   return
end function fy3x2my2
!
!
!-----------------------------------------------------------------------
subroutine radialpart(ng, q, alfa, rvalue, lmax, radial)
  !-----------------------------------------------------------------------
  !
  ! This routine computes a table with the radial Fourier transform 
  ! of the radial functions.
  !
  USE kinds,      ONLY : dp
  USE constants,  ONLY : fpi
  USE cell_base,  ONLY : omega
  !
  implicit none
  ! I/O
  integer :: ng, rvalue, lmax
  real(DP) :: q(ng), alfa, radial(ng,0:lmax)
  ! local variables
  real(DP), parameter :: xmin=-6.d0, dx=0.025d0, rmax=10.d0

  real(DP) :: rad_int, pref, x
  integer :: l, lp1, ir, ig, mesh_r
  real(DP), allocatable :: bes(:), func_r(:), r(:), rij(:), aux(:)

  mesh_r = nint ( ( log ( rmax ) - xmin ) / dx + 1 )
  allocate ( bes(mesh_r), func_r(mesh_r), r(mesh_r), rij(mesh_r) )
  allocate ( aux(mesh_r))
  !
  !    compute the radial mesh
  !
  do ir = 1, mesh_r
     x = xmin  + DBLE (ir - 1) * dx 
     r (ir) = exp (x) / alfa
     rij (ir) = dx  * r (ir)
  enddo
  !
  if (rvalue==1) func_r(:) = 2.d0 * alfa**(3.d0/2.d0) * exp(-alfa*r(:))
  if (rvalue==2) func_r(:) = 1.d0/sqrt(8.d0) * alfa**(3.d0/2.d0) * & 
                     (2.0d0 - alfa*r(:)) * exp(-alfa*r(:)*0.5d0)
  if (rvalue==3) func_r(:) = sqrt(4.d0/27.d0) * alfa**(2.0d0/3.0d0) * &
                     (1.d0 - 1.5d0*alfa*r(:) + 2.d0*(alfa*r(:))**2/27.d0) * &
                                           exp(-alfa*r(:)/3.0d0)
  pref = fpi/sqrt(omega)
  !
  do l = 0, lmax
     do ig=1,ng
       call sph_bes (mesh_r, r(1), q(ig), l, bes)
       aux(:) = bes(:) * func_r(:) * r(:)
       call simpson (mesh_r, aux, rij, rad_int)
       radial(ig,l) = rad_int * pref
     enddo
  enddo

  deallocate (bes, func_r, r, rij, aux )
  return
end subroutine radialpart


