!{\src2tex{textfont=tt}}
!!****m* ABINIT/m_commutator_vkbr
!! NAME
!!  m_commutator_vkbr
!!
!! FUNCTION
!!  This module provides objects and methods used to calculate the matrix elements 
!!  of the commutator [H,r] needed for the correct treatment of the optical limit q-->0
!!  in the matrix elements <k-q,b1|e^{-iqr}|k,b2> when non-local pseudopotentials 
!!  are used or PAW with a nonlocal term in the all-electron Hamiltonian (eg. LDA+U).
!!
!! COPYRIGHT
!! Copyright (C) 2008-2014 ABINIT group (MG, FB)
!! This file is distributed under the terms of the
!! GNU General Public License, see ~abinit/COPYING
!! or http://www.gnu.org/copyleft/gpl.txt .
!!
!! SOURCE

#if defined HAVE_CONFIG_H
#include "config.h"
#endif

#include "abi_common.h"

MODULE m_commutator_vkbr

 use defs_basis
 use defs_datatypes
 use m_blas
 use m_errors
 use m_profiling

 use m_gwdefs,        only : czero_gw
 use m_sphharm,       only : ylmc, ylmcd
 use m_geometry,      only : normv
 use m_crystal,       only : crystal_t

 implicit none

 private

 public ::  ccgradvnl          ! Compute (grad_K+grad_K') Vnl(K,K') using Legendre polynomials.
 public ::  ccgradvnl_ylm      ! Compute (grad_K+grad_K') Vnl(K,K') using complex spherical harmonics.
 public ::  calc_vkb
 public ::  nc_ihr_comm
!!***

!!****t* m_commutator_vkbr/kb_form_factors
!! NAME
!! kb_form_factors
!!
!! FUNCTION
!!  The kb_form_factors data type stores basic dimensions and quantities 
!!  used in the GW part for the treatment of the non-analytic behavior of the 
!!  heads and wings of the irreducible polarizability in the long wave-length limit (i.e. q-->0).
!!
!! SOURCE

 type,public :: kb_form_factors

!scalars
  integer :: npw       
  ! Number of plane waves for this  k-point.

  integer :: lnmax     
  ! Max. number of (l,n) components over all type of pseudos.

  integer :: ntypat    
  ! Number of type of atoms.

!arrays
  integer,allocatable :: sign_dyad(:,:)
  ! sign_dyad(lnmax,ntypat). 
  ! sign of the KB dyadic product.

  real(dp),allocatable :: ff(:,:,:) 
  ! ff(npw,lnmax,ntypat) 
  ! KB form factor.

  real(dp),allocatable :: ffd(:,:,:)
  ! ffd(npw,lnmax,ntypat) 
  ! Derivative of ff wrt k+G.

 end type kb_form_factors
!!***

 !public :: correct_init_kb_ffactors
 public :: destroy_kb_ffactors

!----------------------------------------------------------------------

!!****t* m_commutator_vkbr/kb_potential
!! NAME
!!
!! FUNCTION
!!  Matrix elements in |k+G> space needed for the 
!!  evaluation of the matrix elements of the commutator [Vnl,r] for the 
!!  optical limit in <kb1|e^{-iqr}|kb2>.
!!
!! SOURCE                                                                                   

type,public :: kb_potential

 integer :: istwfk
 ! Storage mode of the G vectors for this k-point.

 integer :: ntypat
 ! Number of type of atoms

 integer :: natom
 ! Number of atoms

 integer :: mpsang
 ! Max l+1 over atoms

 integer :: npwwfn
 ! Number of G-vectors.

 integer :: inclvkb
 ! Option for calculating the matrix elements of [Vnl,r].

 !% integer :: has_fnlr
 ! 1 if fnlr are stored and allocated

 !% integer :: ngfft(18)
 ! Info on the FFT algorithm.

 real(dp) :: kpoint(3)
 ! The k-point in reduced coordinates.

 complex(gwpc),allocatable :: gradvnl(:,:,:)
 ! gradvnl(3,npwwfn,npwwfn)
 ! (grad_K + grad_K') Vnl(K,K') in reciprocal lattice units. used if inclvkb==0 (OBSOLETE)

 complex(gwpc),allocatable :: fnl(:,:,:)
 ! fnl(npwwfn,mpsang*mpsang,natom)

 !% complex(gwpc),allocatable :: fnlr(:,:,:) 
 ! fnlr(nfft,mpsang*mpsang,natom)
 ! Real space representation of fnl.

 complex(gwpc),allocatable :: fnld(:,:,:,:) 
 ! fnld(3,npwwfn,mpsang*mpsang,natom)

end type kb_potential

 public :: init_kb_potential
 public :: destroy_kb_potential
 public :: add_vnlr_commutator
!!**

 interface destroy_kb_potential
   module procedure destroy_kb_potential_0D
   module procedure destroy_kb_potential_1D
 end interface destroy_kb_potential

CONTAINS  !========================================================================================
!!***

!----------------------------------------------------------------------

!!****f* m_commutator_vkbr/correct_init_kb_ffactors
!! NAME
!!  correct_init_kb_ffactors
!!
!! FUNCTION
!!  Calculate KB form factors and derivatives required to evalute
!!  the matrix elements of the commutator [Vnl,r]-. 
!!  This term enters the expression for the oscillator strengths in 
!!  the optical limit q-->0. Pseudopotentials with more than one
!!  projector per angular channel are supported.
!!
!! INPUTS
!!  npw_k=Number of planewaves for this k-point.
!!  kpoint(3)=The kpoint in reduced coordinates.
!!  Psps<pseudopotential_type>Structure gathering info on the pseudopotentials.
!!  rprimd(3,3)=dimensional primitive translations for real space (bohr)
!!  kg_k(3,npw_k)=Set of G-vectors for this k-point.
!!
!! OUTPUT
!!  KBff_k<KB_form_factor>=Structure storing the Kleynmann-Bylander form factor and derivatives for a single k-point.
!!
!! TODO 
!!  Replace old implementation with this new routine. Matrix elements
!!  of the commutator should be calculated on-the-fly in screening only
!!  if really needed. This is the first step toward the elimination
!!  of the KSS file. Modifications in cchi0q0 are needed.
!!
!! PARENTS
!!
!! CHILDREN
!!      add_vnlr_commutator,metric,mkffnl,mkkin
!!
!! SOURCE

subroutine correct_init_kb_ffactors(KBff_k,kpoint,rprimd,npw_k,kg_k,Psps)

 use defs_datatypes

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'correct_init_kb_ffactors'
 use interfaces_41_geometry
 use interfaces_65_nonlocal
!End of the abilint section

 implicit none

!Arguments ------------------------------------
!scalars
 integer,intent(in) :: npw_k
 type(Pseudopotential_type),intent(in) :: Psps
 type(kb_form_factors),intent(inout) :: KBff_k
!arrays
 integer,intent(in) :: kg_k(3,npw_k)
 real(dp),intent(in) :: kpoint(3),rprimd(3,3)

!Local variables ------------------------------
!scalars
 integer :: dimffnl,ider,idir,itypat,nkpg
 integer :: il,ilmn,iln,iln0,nlmn,ig
 real(dp) :: fact,ucvol
 !real(dp) :: ,effmass,ecutsm
 character(len=500) :: msg
!arrays
 real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),kpg(3)
 real(dp),allocatable :: ffnl(:,:,:,:),kpg_dum(:,:),modkplusg(:)
 real(dp),allocatable :: ylm(:,:),ylm_gr(:,:,:),ylm_k(:,:)

! *************************************************************************

 DBG_ENTER("COLL")

 ABI_CHECK(Psps%usepaw==0,"You should not be here!")

 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
 !
 ! === Save KB dyadic sign (integer-valued) ===
 ! * Notice how the ordering is chosen correctly unlike in outkss.
 ! * More than one projector per angular channel is allowed but changes in cchi0q0 are needed
 !  allocate(KBff_ksign(Psps%mpsang,ntypat))  THIS THE OLD IMPLEMENTATION
 ABI_MALLOC(KBff_k%sign_dyad,(Psps%lnmax,Psps%ntypat))
 KBff_k%sign_dyad(:,:)=0

 do itypat=1,Psps%ntypat
   nlmn=COUNT(Psps%indlmn(3,:,itypat)>0)
   iln0=0 
   do ilmn=1,nlmn
     iln=Psps%indlmn(5,ilmn,itypat)
     if (iln>iln0) then
       iln0=iln
       KBff_k%sign_dyad(iln,itypat)=NINT(DSIGN(one,Psps%ekb(ilmn,itypat)))
     end if
   end do
 end do

 KBff_k%npw    = npw_k
 KBff_k%lnmax  = Psps%lnmax
 KBff_k%ntypat = Psps%ntypat

 ! === Allocate KB form factor and derivative wrt k+G ===
 ! * Also here we use correct ordering for dimensions
 ABI_MALLOC(KBff_k%ff ,(npw_k,Psps%lnmax,Psps%ntypat))
 ABI_MALLOC(KBff_k%ffd,(npw_k,Psps%lnmax,Psps%ntypat))
 KBff_k%ff(:,:,:)=zero ; KBff_k%ffd(:,:,:)=zero
 
 ider=1 ; dimffnl=2 ! To retrieve the first derivative.
 idir=0 ; nkpg=0
 !
 ! Quantities used only if useylm==1
 ABI_MALLOC(ylm,(npw_k,Psps%mpsang**2*Psps%useylm))
 ABI_MALLOC(ylm_gr,(npw_k,3+6*(ider/2),Psps%mpsang**2*Psps%useylm))
 ABI_MALLOC(ylm_k,(npw_k,Psps%mpsang**2*Psps%useylm))
 ABI_MALLOC(kpg_dum,(npw_k,nkpg))

 ABI_MALLOC(ffnl,(npw_k,dimffnl,Psps%lmnmax,Psps%ntypat))

 call mkffnl(Psps%dimekb,dimffnl,Psps%ekb,ffnl,Psps%ffspl,gmet,gprimd,ider,idir,Psps%indlmn,&
&  kg_k,kpg_dum,kpoint,Psps%lmnmax,Psps%lnmax,Psps%mpsang,Psps%mqgrid_ff,nkpg,npw_k,& 
&  Psps%ntypat,Psps%pspso,Psps%qgrid_ff,rmet,Psps%usepaw,Psps%useylm,ylm_k,ylm_gr)

 ABI_FREE(kpg_dum)
 ABI_FREE(ylm)
 ABI_FREE(ylm_gr)
 ABI_FREE(ylm_k)

 ABI_MALLOC(modkplusg,(npw_k))

 !effmass=one; ecutsm=zero
 !call mkkin(ecut,ecutsm,effmass,gmet,kg_k,modkplusg,kpoint,npw_k)
 !modkplusg(:)=SQRT(half/pi**2*modkplusg(:))
 !modkplusg(:)=MAX(modkplusg(:),tol10)

 do ig=1,npw_k
   kpg(:)= kpoint(:)+kg_k(:,ig)
   modkplusg(ig) = normv(kpg,gmet,"G")
 end do

 do itypat=1,Psps%ntypat
   nlmn=COUNT(Psps%indlmn(3,:,itypat)>0)
   iln0=0 
   do ilmn=1,nlmn
     il= Psps%indlmn(1,ilmn,itypat)+1
     iln=Psps%indlmn(5,ilmn,itypat)

     if (iln>iln0) then
       iln0=iln

       if (ABS(Psps%ekb(ilmn,itypat))>1.0d-10) then
         SELECT CASE (il)
         CASE (1)
           KBff_k%ff (1:npw_k,iln,itypat) = ffnl(:,1,ilmn,itypat)
           KBff_k%ffd(1:npw_k,iln,itypat) = ffnl(:,2,ilmn,itypat)*modkplusg(:)/two_pi

         CASE (2)
           KBff_k%ff (1:npw_k,iln,itypat) =   ffnl(:,1,ilmn,itypat)*modkplusg(:)
           KBff_k%ffd(1:npw_k,iln,itypat) = ((ffnl(:,2,ilmn,itypat)*modkplusg(:)**2)+&
&                                             ffnl(:,1,ilmn,itypat))/two_pi
         CASE (3)
           KBff_k%ff (1:npw_k,iln,itypat) =  ffnl(:,1,ilmn,itypat)*modkplusg(:)**2
           KBff_k%ffd(1:npw_k,iln,itypat) = (ffnl(:,2,ilmn,itypat)*modkplusg(:)**3+&
&                                          2*ffnl(:,1,ilmn,itypat)*modkplusg(:))/two_pi
         CASE (4)
           KBff_k%ff (1:npw_k,iln,itypat) =  ffnl(:,1,ilmn,itypat)*modkplusg(:)**3
           KBff_k%ffd(1:npw_k,iln,itypat) = (ffnl(:,2,ilmn,itypat)*modkplusg(:)**4+&
&                                          3*ffnl(:,1,ilmn,itypat)*modkplusg(:)**2)/two_pi
         CASE DEFAULT
           msg=' l greater than g not implemented. '
           MSG_ERROR(msg)
         END SELECT

         fact = SQRT(four_pi/ucvol*(2*il-1)*ABS(Psps%ekb(ilmn,itypat)))
         KBff_k%ff (:,iln,itypat) = fact * KBff_k%ff (:,iln,itypat)
         KBff_k%ffd(:,iln,itypat) = fact * KBff_k%ffd(:,iln,itypat)

       else ! ekb==0
         KBff_k%ff (:,iln,itypat)=zero
         KBff_k%ffd(:,iln,itypat)=zero
       end if

     end if
   end do
 end do

 ABI_FREE(ffnl)
 ABI_FREE(modkplusg)

 DBG_EXIT("COLL")

end subroutine correct_init_kb_ffactors
!!***

!----------------------------------------------------------------------

!!****f* m_commutator_vkbr/destroy_kb_ffactors
!! NAME
!!  destroy_kb_ffactors
!!
!! FUNCTION
!!  Free the memory allocated in a structure of type kb_form_factors
!!
!! INPUTS
!!
!! OUTPUT
!!
!! PARENTS
!!
!! CHILDREN
!!      add_vnlr_commutator,metric,mkffnl,mkkin
!!
!! SOURCE

subroutine destroy_kb_ffactors(KBff_k)


!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'destroy_kb_ffactors'
!End of the abilint section

 implicit none

!Arguments ------------------------------------
!scalars
 type(kb_form_factors),intent(inout) :: KBff_k

!************************************************************************

 !@kb_form_factors

!integer 
 if (allocated(KBff_k%sign_dyad)) then
   ABI_FREE(KBff_k%sign_dyad)
 end if

!real
 if (allocated(KBff_k%ff )) then
   ABI_FREE(KBff_k%ff)
 end if
 if (allocated(KBff_k%ffd)) then
   ABI_FREE(KBff_k%ffd)
 end if

end subroutine destroy_kb_ffactors
!!***

!----------------------------------------------------------------------

!!****f* ABINIT/ccgradvnl_ylm
!! NAME
!! ccgradvnl_ylm
!!
!! FUNCTION
!!  Compute Vnl(K) and grad_K Vnl(K) three reciprocal lattice units components
!!  using spherical harmonics instead of Legendre polynomials
!!  Needed for chi0(q=0)
!!
!! INPUTS
!!  Cryst<crystal_t>=Unit cell and symmetries
!!    %rprimd(3,3)=the three primitive vectors
!!    %ntypat=number of types of atoms
!!    %natom=number of atoms
!!    %xcart(3,natom)=cartesian coordinates of nuclei
!!    %typat(natom)=type of each atom
!!  gprimd(3,3)=dimensional primitive translations for reciprocal space ($\textrm{bohr}^{-1}$)
!!  gvec(3,npwwfn)=integer coordinates of each plane wave in reciprocal space
!!  kibz(3,nkibz)=coordinates of all k points in the irreducible Brillouin zone
!!  mpsang=1+maximum angular momentum for nonlocal pseudopotentials
!!  nkibz=number of k points in the irreducible Brillouin zone
!!  npwwfn=number of planewaves for wavefunctions
!!  vkb(npwwfn,ntypat,mpsang,nkibz)=KB projector function
!!  vkbd(npwwfn,ntypat,mpsang,nkibz)=derivative of the KB projector function in reciprocal space
!!  vkbsign(mpsang,ntypat)=sign of each KB dyadic product
!!
!! OUTPUT
!!  l_fnl(npwwfn,mpsang*mpsang,natom,nkibz),
!!  l_fnld(3,npwwfn,mpsang*mpsang,natom,nkibz)
!!
!! NOTES
!!  Subroutine taken from the EXC code  
!!  All the calculations are done in double precision, but the output arrays l_fnl and l_fnld 
!!  are in single precision, should use double precision after modification of the
!!  other subroutines 
!!  the subroutine does not work wity pseudo with more that one projector per angular state TODO
!!
!! PARENTS
!!      m_commutator_vkbr
!!
!! CHILDREN
!!      add_vnlr_commutator,metric,mkffnl,mkkin
!!
!! SOURCE

subroutine ccgradvnl_ylm(npwwfn,Cryst,gvec,kpoint,mpsang,vkbsign,vkb,vkbd,l_fnl,l_fnld)

 use defs_basis

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'ccgradvnl_ylm'
!End of the abilint section

 implicit none

!Arguments ------------------------------------
!scalars
 integer,intent(in) :: mpsang,npwwfn
 type(crystal_t),intent(in) :: Cryst
!arrays
 integer,intent(in) :: gvec(3,npwwfn)
 real(dp),intent(in) :: kpoint(3)
 real(dp),intent(in) :: vkb(npwwfn,Cryst%ntypat,mpsang)
 real(dp),intent(in) :: vkbd(npwwfn,Cryst%ntypat,mpsang),vkbsign(mpsang,Cryst%ntypat)
 complex(gwpc),intent(out) :: l_fnl(npwwfn,mpsang**2,Cryst%natom)
 complex(gwpc),intent(out) :: l_fnld(3,npwwfn,mpsang**2,Cryst%natom)

!Local variables-------------------------------
!scalars
 integer,parameter :: nlx=4
 integer :: i,ia,ig,il,im,iml,is,lmax
 real(dp),parameter :: ppad=tol8
 real(dp) :: cosphi,costh,factor,mkg,mkg2,sinphi,sinth,sq,xdotg
 complex(dpc) :: dphi,dth,sfac
 character(len=500) :: msg
!arrays
 real(dp) :: gcart(3),kcart(3),kg(3)
 real(dp) :: b1(3),b2(3),b3(3),a1(3),a2(3),a3(3)
 complex(dpc) :: dylmcart(3),dylmcrys(3),gradphi(3),gradth(3)
!************************************************************************

 DBG_ENTER("COLL")

 lmax=mpsang
 if (mpsang>nlx) then
   write(msg,'(3a)')&
&   'Number of angular momentum components bigger than programmed.',ch10,&
&   'Taking into account only s p d f ' 
   MSG_WARNING(msg)
   lmax=nlx
 end if

 a1=Cryst%rprimd(:,1) ; b1=two_pi*Cryst%gprimd(:,1)
 a2=Cryst%rprimd(:,2) ; b2=two_pi*Cryst%gprimd(:,2)
 a3=Cryst%rprimd(:,3) ; b3=two_pi*Cryst%gprimd(:,3)
 !
 !
 ! === Calculate Kleiman-Bylander factor and first derivative ===
 l_fnl=czero_gw; l_fnld=czero_gw
 do ig=1,npwwfn
   !  
   ! === Get kcart =k+G in Cartesian coordinates ===
   kg(:)= kpoint(:)+REAL(gvec(:,ig))
   kcart(:)=kg(1)*b1(:)+kg(2)*b2(:)+kg(3)*b3(:)
   ! * Solve the problem with sinth=0. or sinphi=0
   if (ABS(kcart(2))<ppad) kcart(2)=kcart(2)+ppad

   mkg2=kcart(1)**2+kcart(2)**2+kcart(3)**2
   mkg=SQRT(mkg2)
   ! The next to solve the problem with k=Gamma.
   !if (mkg < 0.0001) cycle

   sq=SQRT(kcart(1)**2+kcart(2)**2)

   gcart(:)=  REAL(gvec(1,ig))*b1(:)&
&            +REAL(gvec(2,ig))*b2(:)&
&            +REAL(gvec(3,ig))*b3(:)

   ! === Calculate spherical coordinates (th,phi) ===
   costh = kcart(3)/mkg
   sinth = sq/mkg
   cosphi= kcart(1)/sq
   sinphi= kcart(2)/sq
   
   gradth(1)  = kcart(1)*kcart(3)/mkg**3/sinth
   gradth(2)  = kcart(2)*kcart(3)/mkg**3/sinth
   gradth(3)  = -(one/mkg-kcart(3)**2/mkg**3)/sinth
   gradphi(1) = -(one/sq - kcart(1)**2/sq**3)/sinphi
   gradphi(2) = kcart(2)*kcart(1)/sq**3/sinphi
   gradphi(3) = czero
   
   do ia=1,Cryst%natom
     is=Cryst%typat(ia)
     xdotg=gcart(1)*Cryst%xcart(1,ia)+gcart(2)*Cryst%xcart(2,ia)+gcart(3)*Cryst%xcart(3,ia)
     ! Remember that in the GW code the reciprocal vectors 
     ! are defined such as a_i*b_j = 2pi delta_ij, no need to introduce 2pi
     sfac=CMPLX(COS(xdotg),SIN(xdotg)) 

     do il=1,lmax
       factor=SQRT(four_pi/REAL(2*(il-1)+1))
       do im=1,2*(il-1)+1
         ! Index of im and il
         iml=im+(il-1)*(il-1)
         !     
         ! Calculate the first KB factor, note that l_fnl is simple precision complex
         l_fnl(ig,iml,ia)=factor*sfac*ylmc(il-1,im-il,kcart)*vkb(ig,is,il)*vkbsign(il,is)
         !     
         ! Calculate the second KB factor (involving first derivatives)
         ! dYlm/dK = dYlm/dth * grad_K th + dYlm/dphi + grad_K phi
         call ylmcd(il-1,im-il,kcart,dth,dphi)
         dylmcart(:) = dth*gradth(:) + dphi*gradphi(:)
         !     
         ! Cartesian to crystallographic axis
         ! Notice: a bug was discovered by Marco Cazzaniga, december 2009
         ! the transformation matrix A=(a1,a2,a3) must act *on its left* on the
         ! covariant vector dylmcart (a *row* vector), hence:
         dylmcrys(1)=(a1(1)*dylmcart(1)+a1(2)*dylmcart(2)+a1(3)*dylmcart(3))/(two_pi)
         dylmcrys(2)=(a2(1)*dylmcart(1)+a2(2)*dylmcart(2)+a2(3)*dylmcart(3))/(two_pi)
         dylmcrys(3)=(a3(1)*dylmcart(1)+a3(2)*dylmcart(2)+a3(3)*dylmcart(3))/(two_pi)
         ! the previous implementation assumed A acting on its right
         ! on a column vector, yielding wrong results for the (small)
         ! non local contributions to the spectra,
         ! such as a spurious anisotropy in isotropic systems.
         !
         ! Note that l_fnld is simple precision complex, it could be possible use double precision
         do i=1,3
           l_fnld(i,ig,iml,ia) = factor*sfac*&
&            ( kg(i)/mkg*ylmc(il-1,im-il,kcart)*vkbd(ig,is,il) + dylmcrys(i)*vkb(ig,is,il) )
         end do 
       end do !im
     end do !il
   end do !ia
 end do !ig

 DBG_EXIT("COLL")

end subroutine ccgradvnl_ylm
!!***

!----------------------------------------------------------------------

!!****f* ABINIT/ccgradvnl
!! NAME
!! ccgradvnl
!!
!! FUNCTION
!! Compute the (grad_K+grad_K') Vnl(K,K') (three reciprocal lattice units components)
!! Needed for chi0(q=0)
!!
!! INPUTS
!!  gvec(3,npwwfn)=integer coordinates of each plane wave in reciprocal space
!!  kibz(3,nkibz)=coordinates of all k points in the irreducible Brillouin zone
!!  mpsang=1+maximum angular momentum for nonlocal pseudopotentials
!!  Cryst<crystal_t>=data type gathering information on unit cell and symmetries
!!    %natom=number of atoms
!!    %typat(natom)=type of each atom
!!    %ntypat=number of types of atoms
!!    %xcart(3,natom)=cartesian coordinates of nuclei
!!  nkibz=number of k points in the irreducible Brillouin zone
!!  npwwfn=number of planewaves for wavefunctions
!!  vkb(npwwfn,ntypat,mpsang,nkibz)=KB projector function
!!  vkbd(npwwfn,ntypat,mpsang,nkibz)=derivative of the KB projector function in reciprocal space
!!  vkbsign(mpsang,ntypat)=sign of each KB dyadic product
!!
!! OUTPUT
!!  gradvnl =(grad_K + grad_K') Vnl(K,K') in reciprocal lattice units
!!  
!! PARENTS
!!      m_commutator_vkbr
!!
!! CHILDREN
!!      add_vnlr_commutator,metric,mkffnl,mkkin
!!
!! SOURCE

subroutine ccgradvnl(npwwfn,gvec,kpoint,Cryst,mpsang,vkbsign,vkb,vkbd,gradvnl)

 use defs_basis

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'ccgradvnl'
 use interfaces_14_hidewrite
!End of the abilint section

 implicit none

!Arguments ------------------------------------
!scalars
 integer,intent(in) :: mpsang,npwwfn
 type(crystal_t),intent(in) :: Cryst
!arrays
 integer,intent(in) :: gvec(3,npwwfn)
 real(dp),intent(in) :: kpoint(3),vkb(npwwfn,Cryst%ntypat,mpsang)
 real(dp),intent(in) :: vkbd(npwwfn,Cryst%ntypat,mpsang) 
 real(dp),intent(in) :: vkbsign(mpsang,Cryst%ntypat)
 complex(gwpc),intent(out) :: gradvnl(3,npwwfn,npwwfn)

!Local variables ------------------------------
!scalars
 integer :: ii,ia,ig,igd1,igd2,igd3,igp,il,is,lmax
 real(dp) :: mkg,mkg2,mkgkgp,mkgp,mkgp2,rgdx,rgdy,rgdz,taugd,x,x2,x3,x4,x5,x6,x7,xcostheta
 complex(dpc) :: cs,ct
 character(len=500) :: msg
!arrays
 complex(dpc) :: sfac(Cryst%ntypat)
 integer,parameter :: nlx=8
 real(dp) :: pl(nlx),pld(nlx)
 real(dp) :: kg(3),kgp(3)
 real(dp) :: b1(3),b2(3),b3(3)
!************************************************************************

 call wrtout(std_out,' limit q->0, including term <n,k|[Vnl,iqr]|n"k>','COLL')

 b1=two_pi*Cryst%gprimd(:,1)
 b2=two_pi*Cryst%gprimd(:,2)
 b3=two_pi*Cryst%gprimd(:,3)

 lmax=mpsang
 if(mpsang>nlx) then
   write(msg,'(3a)')&
&    'Number of angular momentum components bigger than programmed! ',ch10,&
&    'Taking into account only s p d f g h i j '
   MSG_WARNING(msg)
   lmax=nlx
 end if
 !
 ! Legendre polynomials and their first derivatives
 ! s p d f g h i j  so up to PL_7 = pl(8)
 !
 pl(1)  = one
 pld(1) = zero
 !pl(2) = costheta
 pld(2) = 1.0
 !pl(3) = 1/2 ( 3 costheta**2 - 1 )
 !pld(3)= 3 costheta

 gradvnl = czero_gw

 do ig=1,npwwfn
   kg(:)=kpoint(:) + real(gvec(:,ig))
   mkg2=scpdt(kg,kg,b1,b2,b3)
   mkg=sqrt(mkg2)
   ! The next to solve the problem with k=Gamma.
   if (mkg < 0.0001) cycle

   do igp=1,npwwfn
     kgp(:)=kpoint(:) + real(gvec(:,igp))
     mkgp2=scpdt(kgp,kgp,b1,b2,b3)
     mkgp=sqrt(mkgp2)
     ! The next to solve the problem with k=Gamma.
     if (mkgp < 0.0001) cycle

     mkgkgp = mkg*mkgp
     xcostheta = scpdt(kg,kgp,b1,b2,b3) / mkgkgp
     x = xcostheta
     x2 = x * x
     x3 = x2 * x
     x4 = x3 * x
     x5 = x4 * x
     x6 = x5 * x
     x7 = x6 * x
     !
     ! Calculate legendre polynomial PL_0 = pl(1)
     pl(2) = x
     pl(3) = (3.0/2.0) * x2 - (1.0/2.0)
     pl(4) = (5.0/2.0) * x3 - (3.0/2.0) * x
     pl(5) = (35.0/8.0) * x4 - (30.0/8.0) * x2 + (3.0/8.0)
     pl(6) = (63.0/8.0) * x5 - (70.0/8.0) * x3 + (15.0/8.0) * x
     pl(7) = (231.0/16.0) * x6 - (315.0/16.0) * x4 + (105.0/16.0) * x2 - (5.0/16.0)
     pl(8) = (429.0/16.0) * x7 - (693.0/16.0) * x5 + (315.0/16.0) * x3 - (35.0/16.0) * x
     !
     ! Calculate legendre polynomial derivative
     pld(3) = 3.0 * x
     pld(4) = (15.0/2.0) * x2 - (3.0/2.0)
     pld(5) = (35.0/2.0) * x3 - (15.0/2.0) * x
     pld(6) = (315.0/8.0) * x4 - (210.0/8.0) * x2 + (15.0/8.0)
     pld(7) = (693.0/8.0) * x5 - (315.0/4.0) * x3 + (105.0/8.0) * x
     pld(8) = (3003.0/16.0) * x6 - (3465.0/16.0) * x4 + (945.0/16.0) * x2 - (35.0/16.0)

     igd1 = gvec(1,ig)-gvec(1,igp)
     igd2 = gvec(2,ig)-gvec(2,igp)
     igd3 = gvec(3,ig)-gvec(3,igp)
     rgdx = igd1*b1(1)+igd2*b2(1)+igd3*b3(1)
     rgdy = igd1*b1(2)+igd2*b2(2)+igd3*b3(2)
     rgdz = igd1*b1(3)+igd2*b2(3)+igd3*b3(3)

     do is=1,Cryst%ntypat
       sfac(is) = czero
       do ia=1,Cryst%natom
         if (Cryst%typat(ia)/=is) CYCLE
         taugd = rgdx*Cryst%xcart(1,ia)+rgdy*Cryst%xcart(2,ia)+ rgdz*Cryst%xcart(3,ia)
         sfac(is) = sfac(is) + cmplx(cos(taugd),-sin(taugd))
       end do
     end do

     do ii = 1, 3
       gradvnl(ii,ig,igp) = 0.0
       do is=1,Cryst%ntypat
         do il = 1, lmax
           ct =( kg(ii)*(1/mkgkgp - xcostheta/mkg2 ) + &
&           kgp(ii)*(1/mkgkgp - xcostheta/mkgp2 ) ) * &
&           pld(il) * vkbsign(il,is) * vkb(ig,is,il) * vkb(igp,is,il)
           
           cs = pl(il) * vkbsign(il,is) * &
&           ( kg(ii)/mkg * vkbd(ig,is,il) * vkb(igp,is,il) + &
&           kgp(ii)/mkgp * vkb(ig,is,il) * vkbd(igp,is,il) )
           
           gradvnl(ii,ig,igp) = gradvnl(ii,ig,igp) + sfac(is) * (ct + cs)
         end do !il
       end do !is
     end do !ii

   end do !igp
 end do !ig

 contains
!!***

!!****f* ccgradvnl/scpdt
!! NAME
!! scpdt
!!
!! FUNCTION
!! Compute scalar product of two vectors u and v in reciprocal space
!!
!! INPUTS
!!  b1(3),b2(3),b3(3)=the three primitive vectors in reciprocal space
!!  u(3),v(3)=the two vectors
!!
!! OUTPUT
!!  function scpdt=scalar product of u and v in reciprocal space
!!
!! PARENTS
!!
!! CHILDREN
!!
!! SOURCE

function scpdt(u,v,b1,b2,b3)

 use defs_basis

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'scpdt'
!End of the abilint section

 implicit none

!Arguments ------------------------------------
!scalars
 real(dp) :: scpdt
!arrays
 real(dp),intent(in) :: b1(3),b2(3),b3(3),u(3),v(3)
! *************************************************************************
 scpdt=&
& (u(1)*b1(1)+u(2)*b2(1)+u(3)*b3(1))*(v(1)*b1(1)+v(2)*b2(1)+v(3)*b3(1))+&
& (u(1)*b1(2)+u(2)*b2(2)+u(3)*b3(2))*(v(1)*b1(2)+v(2)*b2(2)+v(3)*b3(2))+&
& (u(1)*b1(3)+u(2)*b2(3)+u(3)*b3(3))*(v(1)*b1(3)+v(2)*b2(3)+v(3)*b3(3))

 end function scpdt

end subroutine ccgradvnl
!!***

!----------------------------------------------------------------------

!!****f* m_commutator_vkbr/init_kb_potential
!! NAME
!!  init_kb_potential
!!
!! FUNCTION
!!  Creation method the the kb_potential structures datatype.
!!
!! INPUTS
!!  Cryst<crystal_t>=Datatype gathering info on the crystal structure.
!!  Psps<pseudopotential_type>Structure gathering info on the pseudopotentials.
!!  inclvkb=Option defining the algorithm used for the application of [Vnl,r].
!!   == 1 for Legendre polynomials (obsolete, very memory demanding and slow)
!!   == 2 for Spherical harmonics (much faster and less memort demanding).
!!  npwwfn=Number of planewaves in <k+G1|[Vnl,r]|k+G2>
!!  inclvkb=Storage mode for the wavefunctions at this k-point.
!!  kpoint(3)=K-point of interest in reduced coordinates.
!!  gvec(3,npwwfn)=Reduced coordinates of the G-vectors.
!!
!! OUTPUT
!!  KBgrad_k<kb_potential>=Structure containing arrays needed for calculating <\psi_1|[Vnl,r]\psi_2>.
!!    Completely initialized in output.
!!
!! PARENTS
!!      calc_optical_mels,cchi0q0,cchi0q0_intraband,m_eet
!!
!! CHILDREN
!!      add_vnlr_commutator,metric,mkffnl,mkkin
!!
!! SOURCE

subroutine init_kb_potential(KBgrad_k,Cryst,Psps,inclvkb,istwfk,npwwfn,kpoint,gvec)
 
 use defs_basis

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'init_kb_potential'
!End of the abilint section

 implicit none

!Arguments ------------------------------------
!scalars
 integer,intent(in) :: npwwfn,inclvkb,istwfk
 type(crystal_t),intent(in) :: Cryst
 type(kb_potential),intent(inout) :: KBgrad_k
 type(pseudopotential_type),intent(in) :: Psps
!arrays
 integer,intent(in) :: gvec(3,npwwfn)
 real(dp),intent(in) :: kpoint(3)

!Local variables-------------------------------
!scalars
 character(len=500) :: msg
!arrays
 real(dp),allocatable :: vkb(:,:,:),vkbd(:,:,:),vkbsign(:,:)

!************************************************************************

 !@kb_potential
 KBgrad_k%istwfk  = istwfk
 KBgrad_k%ntypat  = Cryst%ntypat
 KBgrad_k%natom   = Cryst%natom
 KBgrad_k%mpsang  = Psps%mpsang
 KBgrad_k%npwwfn  = npwwfn
 KBgrad_k%inclvkb = inclvkb
 KBgrad_k%kpoint  = kpoint
 !
 ! Calculate KB form factors and derivatives.
 !
 ! The following arrays should be allocated using lnmax instead of mpsang
 ! in order to support pseudos with more than projector.
 ! Moreover they should be calculated on-the-fly using calc_vkb
 ! For the moment, we opt for a quick an dirty implementation.
 ABI_MALLOC(vkbsign,(Psps%mpsang,Cryst%ntypat))
 ABI_MALLOC(vkb ,(npwwfn,Cryst%ntypat,Psps%mpsang))
 ABI_MALLOC(vkbd,(npwwfn,Cryst%ntypat,Psps%mpsang))

 call calc_vkb(Psps,kpoint,npwwfn,gvec,Cryst%rprimd,vkbsign,vkb,vkbd)
 
 SELECT CASE (inclvkb)

 CASE (1) ! * Legendre polynomials (CPU and mem ~npwwfn^2) ===

  ! gradvnl = (grad_K+grad_Kp) Vnl(K,Kp) in reciprocal lattice units.
  write(msg,'(a,f12.1)')'out-of-memory gradvnl; Mb= ',3*npwwfn**2*2*gwpc*b2Mb
  ABI_MALLOC(KBgrad_k%gradvnl,(3,npwwfn,npwwfn))
  ABI_CHECK_ALLOC(msg)

  call ccgradvnl(npwwfn,gvec,kpoint,Cryst,Psps%mpsang,vkbsign,vkb,vkbd,Kbgrad_k%gradvnl)

 CASE (2) ! * Complex spherical harmonics (CPU and mem \propto npwwfn).

  write(msg,'(a,f12.1)')'out-of-memory fnl; Mb= ',npwwfn*Psps%mpsang**2*Cryst%natom*2*gwpc*b2Mb
  ABI_MALLOC(KBgrad_k%fnl,(npwwfn,Psps%mpsang**2,Cryst%natom))
  ABI_CHECK_ALLOC(msg)

  write(msg,'(a,f12.1)')'out-of-memory fnld; Mb= ',3*npwwfn*Psps%mpsang**2*Cryst%natom*2*gwpc*b2Mb
  ABI_MALLOC(KBgrad_k%fnld,(3,npwwfn,Psps%mpsang**2,Cryst%natom))
  ABI_CHECK_ALLOC(msg)

  call ccgradvnl_ylm(npwwfn,Cryst,gvec,kpoint,Psps%mpsang,vkbsign,vkb,vkbd,&
&   KBgrad_k%fnl,KBgrad_k%fnld)


 CASE DEFAULT
   write(msg,'(a,i0)')" Wrong value for inclvkb= ",inclvkb
   MSG_ERROR(msg)
 END SELECT 

 ABI_FREE(vkbsign)
 ABI_FREE(vkb)
 ABI_FREE(vkbd)

end subroutine init_kb_potential 
!!***

!----------------------------------------------------------------------

!!****f* m_commutator_vkbr/destroy_kb_potential_0D
!! NAME
!!  destroy_kb_potential_0D
!!
!! FUNCTION
!!  Free all memory allocated in a structure of type kb_potential
!!
!! PARENTS
!!      m_commutator_vkbr
!!
!! CHILDREN
!!      add_vnlr_commutator,metric,mkffnl,mkkin
!!
!! SOURCE

subroutine destroy_kb_potential_0D(KBgrad_k)
 
 use defs_basis

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'destroy_kb_potential_0D'
!End of the abilint section

 implicit none

!Arguments ------------------------------------
!scalars
 type(kb_potential),intent(inout) :: KBgrad_k

!************************************************************************

 !@kb_potential

!complex
 if (allocated(KBgrad_k%gradvnl)) then
   ABI_FREE(KBgrad_k%gradvnl)
 end if
 if (allocated(KBgrad_k%fnl)) then
   ABI_FREE(KBgrad_k%fnl)
 end if
 if (allocated(KBgrad_k%fnld)) then
   ABI_FREE(KBgrad_k%fnld)
 end if

end subroutine destroy_kb_potential_0D
!!***

!----------------------------------------------------------------------

!!****f* m_commutator_vkbr/destroy_kb_potential_1D
!! NAME
!!  destroy_kb_potential_1D
!!
!! FUNCTION
!!  Free all memory allocated in a structure of type kb_potential
!!
!! PARENTS
!!
!! CHILDREN
!!      add_vnlr_commutator,metric,mkffnl,mkkin
!!
!! SOURCE

subroutine destroy_kb_potential_1D(KBgrad_k)
 
 use defs_basis

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'destroy_kb_potential_1D'
!End of the abilint section

 implicit none

!Arguments ------------------------------------
!arrays
 type(kb_potential),intent(inout) :: KBgrad_k(:)

!Local variables ------------------------------
!scalars
 integer :: ii

!************************************************************************

 do ii=1,SIZE(KBgrad_k)
   call destroy_kb_potential_0D(KBgrad_k(ii))
 end do

end subroutine destroy_kb_potential_1D
!!***

!----------------------------------------------------------------------

!!****f* m_commutator_vkbr/add_vnlr_commutator
!! NAME
!!  add_vnlr_commutator
!!
!! FUNCTION
!!  Calculate the matrix elements of the dipole operator <phi1|r|phi2>.
!!  For norm conserving potentials the commutator [Vnl,r] is included according to inclvkb. 
!!
!! INPUTS
!!  KBgrad_k<kb_potential>
!!  npwwfn=Number of G for wavefunctions.
!!  nspinor=Number of spinorial components.
!!  wfg1(nspinor*npwwfn)=Left wavefunction.
!!  wfg2(nspinor*npwwfn)=Right wavefunction
!!
!! SIDE EFFECTS
!!  rhotwx(3,nspinor**2)= Updated. Matrix elements in reduced coordinates, see NOTES below.
!!
!! NOTES
!!   1) <k b1|e^{-iq.r}|k b2> = \delta_{b1 b2} -iq <k b1|r|k b2> =  \delta_{b1 b2} -iq ( <k b1| [H,r] |k b2> / (e1-e2) ).
!!
!!      This routine calculates the matrix elements of ir*(e1-e2) 
!!      Remember that [H,r] = -\nabla + [V_nl,r]
!!
!!  2) The Fourier transform of a two-point real function f(r1,r2) satisfies:
!!      a) f_{\Gamma}(G1,G2) = f_{\Gamma}(-G1,-G2)^*
!!      b) f_{G0/2}  (G1,G2) = f_{G0/2}(-G1-G0,-G2-G0)^*
!!
!! TODO
!!  *) Spinorial case is not implemented.
!!
!! PARENTS
!!      m_commutator_vkbr
!!
!! CHILDREN
!!      add_vnlr_commutator,metric,mkffnl,mkkin
!!
!! SOURCE

subroutine add_vnlr_commutator(KBgrad_k,npwwfn,nspinor,wfg1,wfg2,rhotwx)

 use defs_basis

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'add_vnlr_commutator'
!End of the abilint section

 implicit none

!Arguments ------------------------------------
!scalars
 integer,intent(in) :: npwwfn,nspinor
 type(kb_potential),intent(in) :: KBgrad_k
!arrays
 complex(gwpc),target,intent(in) :: wfg1(npwwfn*nspinor),wfg2(npwwfn*nspinor)
 complex(gwpc),intent(inout) :: rhotwx(3,nspinor**2)

!Local variables ------------------------------
!scalars
 integer :: ig1,ig2,iat,ig,ilm
 complex(gwpc) :: cta1,cta4,ct
 character(len=500) :: msg
!arrays
 complex(gwpc) :: dum(3),cta2(3),cta3(3),gamma_term(3)

!************************************************************************

 ! === Adding term i <c,k|[Vnl,r]|v,k> === 
 ! * Two different algorithms are coded, the second one is much faster
 if (KBgrad_k%istwfk/=1) then 
   MSG_ERROR("istwfk/=1 not tested")
 end if

 SELECT CASE (KBgrad_k%inclvkb)

 CASE (1) ! Legendre Polynomials. CPU and MEM ~ npwwfn**2 ===
          ! res(1:3)= sum_{G1,G2} c(G1)^\dagger c(G2) V^{nl}_{1:3,G1,G2}
  dum(:)=czero_gw
  do ig1=1,npwwfn
    do ig2=1,npwwfn
      ct=CONJG(wfg1(ig1))*wfg2(ig2)
      dum(:)=dum(:)+ct*KBgrad_k%gradvnl(:,ig1,ig2)
    end do
  end do

  if (KBgrad_k%istwfk>1) then
    dum=two*j_dpc*AIMAG(dum)
    if (KBgrad_k%istwfk==2) dum=dum - CONJG(wfg1(1))*wfg2(1) * KBgrad_k%gradvnl(:,1,1)
  end if

  rhotwx(:,1)=rhotwx(:,1)+dum(:)

 CASE (2)  ! Complex spherical harmonics (much faster!).

  dum=czero_gw; gamma_term=czero
!TODO this section causes segmentation faults
  !$OMP PARALLEL DO PRIVATE(cta1,cta2,cta3,cta4) COLLAPSE(2) REDUCTION(+:dum,gamma_term)
  do iat=1,KBgrad_k%natom
    do ilm=1,KBgrad_k%mpsang**2
      cta1=czero_gw; cta2(:)=czero_gw
      cta4=czero_gw; cta3(:)=czero_gw
      do ig=1,npwwfn ! Here we take advantage of the property Y_(l-m)= (-i)^m Y_lm^*.
        cta1   = cta1    + wfg1(ig) * KBgrad_k%fnl (ig,ilm,iat)
        cta2(:)= cta2(:) + wfg2(ig) * KBgrad_k%fnld(:,ig,ilm,iat)
        cta3(:)= cta3(:) + wfg1(ig) * KBgrad_k%fnld(:,ig,ilm,iat)
        cta4   = cta4    + wfg2(ig) * KBgrad_k%fnl (ig,ilm,iat)
        if (ig==1) gamma_term = gamma_term + CONJG(cta1)*cta2(:) +CONJG(cta3(:))*cta4
      end do
      dum(:)= dum(:) +CONJG(cta1)*cta2(:) +CONJG(cta3(:))*cta4
    end do
  end do

  if (KBgrad_k%istwfk>1) then
    dum=two*j_dpc*AIMAG(dum); if (KBgrad_k%istwfk==2) dum=dum-j_dpc*AIMAG(gamma_term)
  end if

  rhotwx(:,1)=rhotwx(:,1)+dum(:)

 CASE DEFAULT
   write(msg,'(a,i0)')" Wrong value for inclvkb= ",KBgrad_k%inclvkb
   MSG_ERROR(msg)
 END SELECT

end subroutine add_vnlr_commutator
!!***

!----------------------------------------------------------------------

!!****f* m_commutator_vkbr/calc_vkb
!! NAME
!!  calc_vkb
!!
!! FUNCTION
!!  This routine calculates the Kleynman-Bylander form factors and its derivatives 
!!  needed for the evaluation of the matrix elements of the dipole operator <phi1|r|phi2>.
!!
!! INPUTS
!!  npw_k=Number of plane waves for this k-point.
!!  Psps<pseudopotential_type>=Structured datatype gathering information on the pseudopotentials.
!!  kg_k(3,npw_k)=Reduced coordinates of the G-vectors.
!!  kpoint(3)=The k-point in reduced coordinates.
!!  rprimd(3,3)=dimensional primitive translations for real space (bohr)
!!
!! OUTPUT
!!  vkb (npw_k,Psps%ntypat,Psps%mpsang)=KB form factors.
!!  vkbd(npw_k,Psps%ntypat,Psps%mpsang)=KB form factor derivatives.
!!  vkbsign(Psps%mpsang,Psps%ntypat)   =KS dyadic sign.
!!
!! NOTES
!!  This piece of code has been extracted from outkss.F90. The implementation is consistent
!!  with the KSS file formata (Fortran version) but it presents two design flaws.
!!
!!   1) Pseudo with more that one projector per l-channel are not supported.
!!   2) Ordering of dimensions in vkb and vkbd is not optimal. We are not programming C!!!
!!
!! TODO
!!  *) Spinorial case is not implemented.
!!  *) Fix the above mentioned programming sins (KSS FORTRAN fileformat has to be modified though)
!!
!! PARENTS
!!      m_commutator_vkbr,m_io_kss
!!
!! CHILDREN
!!      add_vnlr_commutator,metric,mkffnl,mkkin
!!
!! SOURCE

subroutine calc_vkb(Psps,kpoint,npw_k,kg_k,rprimd,vkbsign,vkb,vkbd)


!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'calc_vkb'
 use interfaces_41_geometry
 use interfaces_56_recipspace
 use interfaces_65_nonlocal
!End of the abilint section

 implicit none

!Arguments ------------------------------------
!scalars
 integer,intent(in) :: npw_k
 type(Pseudopotential_type),intent(in) :: Psps
!arrays
 integer,intent(in) :: kg_k(3,npw_k)
 real(dp),intent(in) :: kpoint(3),rprimd(3,3)
 real(dp),intent(out) :: vkb (npw_k,Psps%ntypat,Psps%mpsang)
 real(dp),intent(out) :: vkbd(npw_k,Psps%ntypat,Psps%mpsang)
 real(dp),intent(out) :: vkbsign(Psps%mpsang,Psps%ntypat)

!Local variables ------------------------------
!scalars
 integer :: dimffnl,ider,idir,itypat,nkpg,il0,in
 integer :: il,ilmn,ig,is
 real(dp) :: ucvol,effmass,ecutsm,ecut
!arrays
 real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
 real(dp),allocatable :: ffnl(:,:,:,:),kpg_dum(:,:),modkplusg(:)
 real(dp),allocatable :: ylm(:,:),ylm_gr(:,:,:),ylm_k(:,:)

! *************************************************************************

 DBG_ENTER("COLL")

 ABI_CHECK(Psps%usepaw==0,"You should not be here!")
 ABI_CHECK(Psps%useylm==0,"useylm/=0 not considered!")

 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
 !
 ! === Save KB dyadic sign (integer-valued) ===
 vkbsign=zero
 do itypat=1,Psps%ntypat
   il0=0 
   do ilmn=1,Psps%lmnmax
     il=1+Psps%indlmn(1,ilmn,itypat)
     in=Psps%indlmn(3,ilmn,itypat)
     if (il/=il0 .and. in==1) then
       il0=il
       vkbsign(il,itypat)=DSIGN(one,Psps%ekb(ilmn,itypat))
     end if
   end do
 end do

 ! === Allocate KB form factor and derivative wrt k+G ===
 ! * Here we do not use correct ordering for dimensions
 
 ider=1; dimffnl=2 ! To retrieve the first derivative.
 idir=0; nkpg=0
 !
 ! Quantities used only if useylm==1
 ABI_MALLOC(ylm,(npw_k,Psps%mpsang**2*Psps%useylm))
 ABI_MALLOC(ylm_gr,(npw_k,3+6*(ider/2),Psps%mpsang**2*Psps%useylm))
 ABI_MALLOC(ylm_k,(npw_k,Psps%mpsang**2*Psps%useylm))
 ABI_MALLOC(kpg_dum,(npw_k,nkpg))

 ABI_MALLOC(ffnl,(npw_k,dimffnl,Psps%lmnmax,Psps%ntypat))

 call mkffnl(Psps%dimekb,dimffnl,Psps%ekb,ffnl,Psps%ffspl,gmet,gprimd,ider,idir,Psps%indlmn,&
&  kg_k,kpg_dum,kpoint,Psps%lmnmax,Psps%lnmax,Psps%mpsang,Psps%mqgrid_ff,nkpg,npw_k,& 
&  Psps%ntypat,Psps%pspso,Psps%qgrid_ff,rmet,Psps%usepaw,Psps%useylm,ylm_k,ylm_gr)

 ABI_FREE(kpg_dum)
 ABI_FREE(ylm)
 ABI_FREE(ylm_gr)
 ABI_FREE(ylm_k)

 ABI_MALLOC(modkplusg,(npw_k))

 effmass=one; ecutsm=zero; ecut=HUGE(one)
 call mkkin(ecut,ecutsm,effmass,gmet,kg_k,modkplusg,kpoint,npw_k)
 modkplusg(:)=SQRT(half/pi**2*modkplusg(:))
 modkplusg(:)=MAX(modkplusg(:),tol10)

 !do ig=1,npw_k
 ! kpg(:)= kpoint(:)+kg_k(:,ig)
 ! modkplusg(ig) = normv(kpg,gmet,"G")
 !end do

 ! Calculate matrix elements.
 vkb=zero; vkbd=zero

 do is=1,Psps%ntypat
   il0=0
   do ilmn=1,Psps%lmnmax
     il=1+Psps%indlmn(1,ilmn,is)
     in=Psps%indlmn(3,ilmn,is)
     if ((il/=il0).and.(in==1)) then
       il0=il
       if (ABS(Psps%ekb(ilmn,is))>1.0d-10) then
         if (il==1) then
           vkb (1:npw_k,is,il) = ffnl(:,1,ilmn,is)
           vkbd(1:npw_k,is,il) = ffnl(:,2,ilmn,is)*modkplusg(:)/two_pi
         else if (il==2) then
           vkb(1:npw_k,is,il)  = ffnl(:,1,ilmn,is)*modkplusg(:)
           do ig=1,npw_k
             vkbd(ig,is,il) = ((ffnl(ig,2,ilmn,is)*modkplusg(ig)*modkplusg(ig))+&
&             ffnl(ig,1,ilmn,is) )/two_pi
           end do
         else if (il==3) then
           vkb (1:npw_k,is,il) =  ffnl(:,1,ilmn,is)*modkplusg(:)**2
           vkbd(1:npw_k,is,il) = (ffnl(:,2,ilmn,is)*modkplusg(:)**3+&
&           2*ffnl(:,1,ilmn,is)*modkplusg(:) )/two_pi
         else if (il==4) then
           vkb (1:npw_k,is,il) =  ffnl(:,1,ilmn,is)*modkplusg(:)**3
           vkbd(1:npw_k,is,il) = (ffnl(:,2,ilmn,is)*modkplusg(:)**4+&
&           3*ffnl(:,1,ilmn,is)*modkplusg(:)**2 )/two_pi
         end if
         vkb (:,is,il) = SQRT(4*pi/ucvol*(2*il-1)*ABS(Psps%ekb(ilmn,is)))*vkb (:,is,il)
         vkbd(:,is,il) = SQRT(4*pi/ucvol*(2*il-1)*ABS(Psps%ekb(ilmn,is)))*vkbd(:,is,il)
       else
         vkb (:,is,il)=zero
         vkbd(:,is,il)=zero
       end if
     end if
   end do
 end do

 ABI_FREE(ffnl)
 ABI_FREE(modkplusg)

 DBG_EXIT("COLL")

end subroutine calc_vkb
!!***

!----------------------------------------------------------------------

!!****f* m_commutator_vkbr/nc_ihr_comm
!! NAME
!!  nc_pwihr_comm
!!
!! FUNCTION
!!  Calculate the matrix elements of the commutator i[H,r]  
!!  For norm conserving potentials the commutator i[Vnl,r] is included depending on inclvkb. 
!!
!! INPUTS
!!  nspinor=Number of spinorial components.
!!  npw=Number of G for wavefunctions.
!!  istwfk=Storage mode for wavefunctions.
!!  inclvkb=Option defining whether [Vnl,r] is added or not.
!!  kpoint(3)=k-point in reduced coordinates.
!!  KBgrad_k<kb_potential>
!!  ug1(npw*nspinor)=Left wavefunction.
!!  ug2(npw*nspinor)=Right wavefunction
!!  gvec(3,npw)=Planes waves for wavefunctions.
!!
!! OUTPUT
!!  ihr_comm(3,nspinor**2)= Matrix elements of the commutator i[H,r] between the input states.
!!   Result is in reduced coordinates. ug1 and ug2 are supposed to be orthogonal. 
!!
!! NOTES
!!  <k b1|e^{-iq.r}|k b2> = \delta_{b1 b2} -iq <k b1|r|k b2> =  \delta_{b1 b2} -iq ( <k b1| [H,r] |k b2> / (e1-e2) ).
!!  Remember that [H,r] = -\nabla + [V_nl,r]
!!
!! TODO
!!  *) Spinorial case is not implemented.
!!
!! PARENTS
!!
!! CHILDREN
!!
!! SOURCE

function nc_ihr_comm(nspinor,npw,istwfk,inclvkb,kpoint,KBgrad_k,ug1,ug2,gvec) result(ihr_comm)

 use defs_basis

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'nc_ihr_comm'
!End of the abilint section

 implicit none

!Arguments ------------------------------------
!scalars
 integer,intent(in) :: npw,nspinor,inclvkb,istwfk
 type(kb_potential),intent(in) :: KBgrad_k
!arrays
 integer,intent(in) :: gvec(3,npw)
 real(dp),intent(in) :: kpoint(3)
 complex(gwpc),intent(in) :: ug1(npw*nspinor),ug2(npw*nspinor)
 complex(gwpc) :: ihr_comm(3,nspinor**2)

!Local variables ------------------------------
!scalars
 integer :: ig,iab,spad1,spad2
 character(len=500) :: msg
 complex(dpc) :: c_tmp
!arrays
 integer :: spinorwf_pad(2,4)

!************************************************************************

 ! [H, r] = -\nabla + [V_{nl}, r]
 ! V_nl is present only in the case of NC pseudos but 
 ! not in PAW unless even the AE Hamiltonian in non-local e.g. LDA+U or LEXX.

 ! -i <c,k|\nabla_r|v,k> in reduced coordinates is always included.
 ! -i <c,k|\nabla_r|v,k> = \sum_G u_{ck}^*(G) [k+G] u_{vk}(G)
 ! Note that here we assume c/=v, moreover the ug are supposed to be orthonormal and 
 ! hence k+G can be replaced by G.

 spinorwf_pad=RESHAPE((/0,0,npw,npw,0,npw,npw,0/),(/2,4/))

 ihr_comm = czero

 if (istwfk==1) then
   do iab=1,nspinor**2
     spad1 = spinorwf_pad(1,iab)
     spad2 = spinorwf_pad(2,iab)
     do ig=1,npw 
       c_tmp = CONJG(ug1(ig+spad1)) * ug2(ig+spad2)
       ihr_comm(:,iab) = ihr_comm(:,iab) + c_tmp*gvec(:,ig)
     end do
   end do
 else
   ! Symmetrized expression: \sum_G  (k+G) 2i Im [ u_a^*(G) u_b(G) ]. (k0,G0) term is null.
   do ig=1,npw 
     c_tmp = CONJG(ug1(ig)) * ug2(ig)
     ihr_comm(:,1) = ihr_comm(:,1) + two*j_dpc * AIMAG(c_tmp) * (kpoint + gvec(:,ig))
   end do  
 end if

 if (inclvkb/=0) then ! Add second term $i <c,k|[Vnl,r]|v,k> in$ reduced cordinates.
   if (nspinor/=1) then
     MSG_ERROR("nspinor/=1 not coded")
   end if
   !if (istwfk/=1) then
   !  MSG_ERROR("Not Implemented Error")
   !end if
   if (istwfk/=KBgrad_k%istwfk) then
     write(msg,"(a,2(i0,1x))")" input istwfk /= KBgrad_k%istwfk ",istwfk,KBgrad_k%istwfk
     MSG_ERROR(msg)
   end if
   call add_vnlr_commutator(KBgrad_k,npw,nspinor,ug1,ug2,ihr_comm)
 end if

end function nc_ihr_comm
!!***

!----------------------------------------------------------------------

END MODULE m_commutator_vkbr
!!***
