!{\src2tex{textfont=tt}}
!!****f* ABINIT/uderiv
!! NAME
!! uderiv
!!
!! FUNCTION
!! This routine is called in scfcv.f to compute the derivative of
!! ground-state wavefunctions with respect to k (du/dk) by finite differencing
!! on neighbouring k points
!! Work for nsppol=1 or 2, but only accept nspinor=1,
!!
!! COPYRIGHT
!! Copyright (C) 2001-2014 ABINIT group (NSAI).
!!
!! INPUTS
!!  bdberry(4)=band limits for Berry phase contributions (or du/dk)
!!   spin up and spin down (bdberry(3:4) is irrelevant when nsppol=1)
!!  cg(2,mcg)=planewave coefficients of wavefunctions
!!  gprimd(3,3)=reciprocal space dimensional primitive translations
!!  hdr <type(hdr_type)>=the header of wf, den and pot files
!!  istwfk(nkpt_)=input option parameter that describes the storage of wfs
!!  kberry(3,20)= different delta k for Berry phases(or du/dk),
!!   in unit of kptrlatt only kberry(1:3,1:nberry) is relevant
!!  kg(3,mpw*mkmem)=reduced planewave coordinates
!!  kpt_(3,nkpt_)=reduced coordinates of k points generated by ABINIT,
!!               kpt_ samples half the BZ if time-reversal symetrie is used
!!  kptopt=2 when time-reversal symmetry is used
!!  kptrlatt(3,3)=k-point lattice specification
!!  mband=maximum number of bands
!!  mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
!!  mkmem=number of k points treated by this node.
!!  mpi_enreg=informations about MPI parallelization
!!  mpw=maximum dimensioned size of npw
!!  natom=number of atoms in cell
!!  nband(nkpt*nsppol)=number of bands at each k point, for each polarization
!!  nberry=number of Berry phases(or du/dk) to be computed
!!  nkpt=number of k points
!!  npwarr(nkpt)=number of planewaves in basis at this k point
!!  nspinor=number of spinorial components of the wavefunctions (on current proc)
!!  nsppol=1 for unpolarized, 2 for spin-polarized
!!  unddk=unit number for ddk file
!!
!! OUTPUT
!!  (the ddk wavefunctions are written on disk)
!!
!! SIDE EFFECTS
!!
!! TODO
!!  Cleaning, checking for rules
!!  Should allow for time-reversal symmetry (istwfk)
!!  WARNING : the use of nspinor is completely erroneous
!!
!! NOTES
!! Local Variables:
!!  cmatrix(:,:,:)= overlap matrix of size maxband*maxband
!!  cg_index(:,:,:)= unpacked cg index array for specific band,
!!   k point and polarization.
!!  det(2,2)= intermediate output of Lapack routine zgedi.f
!!  dk(3)= step taken to the next k mesh point along the kberry direction
!!  gpard(3)= dimensionalreciprocal lattice vector G along which the
!!          polarization is computed
!!  kg_kpt(:,:,:)= unpacked reduced planewave coordinates with subscript of
!!          planewave and k point
!!  kpt(3,nkpt)=reduced coordinates of k-point grid that samples the whole BZ
!!  kpt_flag(nkpt)=kpt_flag(ikpt)=0 when the wf was generated by the ABINIT
!!                 code
!!                 kpt_flag(ikpt) gives the indices of the k-point
!!                 related to ikpt by time revers
!!  maxband/minband= control the minimum and maximum band calculated in the
!!           overlap matrix
!!  npw_k= npwarr(ikpt), number of planewaves in basis at this k point
!!  shift_g_2(nkpt,nkpt)= .true. if the k point should be shifted by a G vector;
!!  .false. if not
!!  tr(2)=variable that changes k to -k
!!                              G to -G
!!                     $c_g$ to $c_g^*$ when time-reversal symetrie is used
!!
!! PARENTS
!!      elpolariz
!!
!! CHILDREN
!!      appdig,dzgedi,dzgefa,hdr_io,matr3inv,rwwf,waveformat,wffclose,wffopen
!!      wrtout,xdefineoff
!!
!! SOURCE

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

#include "abi_common.h"

subroutine uderiv(bdberry,cg,gprimd,hdr,istwfk,kberry,kg,kpt_,kptopt,kptrlatt,&
& mband,mcg,mkmem,mpi_enreg,mpw,natom,nband,nberry,npwarr,nspinor,nsppol,nkpt_,&
& unddk,fnameabo_1wf)

 use defs_basis
 use defs_abitypes
 use m_errors
 use m_wffile
 use m_errors
 use m_profiling

 use m_header,       only : hdr_skip, hdr_io

!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 'uderiv'
 use interfaces_14_hidewrite
 use interfaces_28_numeric_noabirule
 use interfaces_32_util
 use interfaces_51_manage_mpi
 use interfaces_56_io_mpi
 use interfaces_67_common, except_this_one => uderiv
!End of the abilint section

 implicit none

!Arguments ------------------------------------
!scalars
 integer,intent(in) :: kptopt,mband,mcg,mkmem,mpw,natom,nberry,nkpt_,nspinor
 integer,intent(in) :: nsppol,unddk
 type(MPI_type),intent(inout) :: mpi_enreg
 type(hdr_type),intent(inout) :: hdr
!arrays
 integer,intent(in) :: bdberry(4),istwfk(nkpt_),kberry(3,20),kg(3,mpw*mkmem)
 integer,intent(in) :: kptrlatt(3,3),nband(nkpt_*nsppol),npwarr(nkpt_)
 real(dp),intent(in) :: cg(2,mcg),gprimd(1:3,1:3)
 real(dp),intent(in) :: kpt_(3,nkpt_)
 character(len=fnlen),intent(in) :: fnameabo_1wf

!Local variables -------------------------
!scalars
 integer :: accesswff,band_in,cg_index_iband,fform,flag1
 integer :: formeig,iband,iberry,icg,idir,ierr,ifor,ii,ikpt,ikpt2,ikpt_
 integer :: index,index1,info,ipert,ipw,isgn,isppol,jband,jj,jkpt,jkpt_
 integer :: master,maxband,mcg_disk,me,minband,nband_diff,nband_k
 integer :: nkpt,npw_k,pertcase,rdwr,read_k,spaceComm
 integer :: tim_rwwf
 real(dp) :: gmod,twodk
 character(len=500) :: message
 character(len=fnlen) :: fiwf1o
 type(wffile_type) :: wffddk
!arrays
 integer :: kg_jl(0,0,0),kpt_flag(2*nkpt_)
 integer,allocatable :: cg_index(:,:,:),ikpt_dk(:,:),ipvt(:)
 integer,allocatable :: kg_kpt(:,:,:)
 real(dp) :: det(2,2),diffk(3),diffk2(3),dk(3),gpard(3),klattice(3,3)
 real(dp) :: kptrlattr(3,3),tr(2)
 real(dp) :: cg_disk(0,0,0)
 real(dp),allocatable :: cmatrix(:,:,:),dudk(:,:)
 real(dp),allocatable :: eig_dum_2(:),kpt(:,:)
 real(dp),allocatable :: occ_dum_2(:),phi(:,:,:),u_tilde(:,:,:,:),zgwork(:,:)
 logical,allocatable :: shift_g_2(:,:)

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


 if(min(2,(1+mpi_enreg%paral_spinor)*nspinor)==2)then
   message = ' uderiv : does not yet work for nspinor=2'
   MSG_ERROR(message)
 end if

 if(maxval(istwfk(:))/=1)then
   write(message,'(3a)')&
&   '  Sorry, this routine does not work yet with istwfk/=1.',ch10,&
&   '  This should have been tested previously ...'
   MSG_BUG(message)
 end if

 if (kptopt==3) then
   nkpt = nkpt_
   ABI_ALLOCATE(kpt,(3,nkpt))
   kpt(:,:)=kpt_(:,:)
 else if (kptopt==2) then
   nkpt = nkpt_*2
   ABI_ALLOCATE(kpt,(3,nkpt))
   do ikpt = 1,nkpt/2
     kpt_flag(ikpt) = 0
     kpt(:,ikpt)=kpt_(:,ikpt)
   end do
   index = 0
   do ikpt = (nkpt/2+1),nkpt
     flag1 = 0
     do jkpt = 1, nkpt/2
       if (((abs(kpt_(1,ikpt-nkpt/2)+kpt_(1,jkpt))<1.0d-8).or.&
&       (abs(1-abs(kpt_(1,ikpt-nkpt/2)+kpt_(1,jkpt)))<1.0d-8))&
&       .and.((abs(kpt_(2,ikpt-nkpt/2)+kpt_(2,jkpt))<1.0d-8).or.&
&       (abs(1-abs(kpt_(2,ikpt-nkpt/2)+kpt_(2,jkpt)))<1.0d-8))&
&       .and.((abs(kpt_(3,ikpt-nkpt/2)+kpt_(3,jkpt))<1.0d-8).or.&
&       (abs(1-abs(kpt_(3,ikpt-nkpt/2)+kpt_(3,jkpt)))<1.0d-8))) then
         flag1 = 1
         index = index + 1
         exit
       end if
     end do
     if (flag1==0) then
       kpt_flag(ikpt-index)=ikpt-nkpt/2
       kpt(:,ikpt-index)=-kpt_(:,ikpt-nkpt/2)
     end if
   end do
   nkpt = nkpt - index
 end if

!DEBUG
!write(101,*) 'beginning write kpt'
!do ikpt=1,nkpt
!write(101,*) kpt(:,ikpt)
!end do
!ENDDEBUG

 ABI_ALLOCATE(shift_g_2,(nkpt,nkpt))

!Compute primitive vectors of the k point lattice
!Copy to real(dp)
 kptrlattr(:,:)=kptrlatt(:,:)
!Go to reciprocal space (in reduced coordinates)
 call matr3inv(kptrlattr,klattice)

 do iberry=1,nberry

!  **************************************************************************
!  Determine the appended index for ddk 1WF files

   do idir=1,3
     if (kberry(idir,iberry) ==1) then
       ipert=natom+1
       pertcase=idir+(ipert-1)*3
     end if
   end do

!  open ddk 1WF file
   formeig=1

   call appdig(pertcase,fnameabo_1wf,fiwf1o)

   spaceComm=abinit_comm_serial ; me=0 ; master=0 ; accesswff=IO_MODE_FORTRAN
   call WffOpen(accesswff,spaceComm,fiwf1o,ierr,wffddk,master,me,unddk)

   rdwr=2 ; fform=2

   if (wffddk%accesswff /= IO_MODE_NETCDF) then
     call hdr_io(fform,hdr,rdwr,wffddk)
   else if (wffddk%accesswff==IO_MODE_ETSF) then
     MSG_ERROR("Missing ETSF I/O support in uderiv")
   else 
     MSG_ERROR("Unknown accesswff")
   end if

!  Define offsets, in case of MPI I/O
   call xdefineOff(formeig,wffddk,mpi_enreg,nband,npwarr,nspinor,nsppol,nkpt_)

!  *****************************************************************************
!  Calculate dimensional recip lattice vector along which P is calculated
!  dk =  step to the nearest k point along that direction
!  in reduced coordinates

   dk(:)=dble(kberry(1,iberry))*klattice(:,1)+&
&   dble(kberry(2,iberry))*klattice(:,2)+&
&   dble(kberry(3,iberry))*klattice(:,3)

   do idir=1,3
     if (dk(idir)/=0) then
       twodk=2*dk(idir)
     end if
   end do

   gpard(:)=dk(1)*gprimd(:,1)+dk(2)*gprimd(:,2)+dk(3)*gprimd(:,3)
   gmod=sqrt(dot_product(gpard,gpard))

!  ******************************************************************************
!  Select the k grid  points along the direction to compute dudk
!  dk =  step to the nearest k point along that direction

!  For each k point, find k_prim such that k_prim= k + dk mod(G)
!  where G is a vector of the reciprocal lattice
   ABI_ALLOCATE(ikpt_dk,(2,nkpt))
   ikpt_dk(1:2,1:nkpt)=0
   shift_g_2(:,:)= .false.

   do ikpt=1,nkpt
     do ikpt2=1,nkpt
       diffk(:)=abs(kpt(:,ikpt2)-kpt(:,ikpt)-dk(:))
       diffk2(:)=abs(kpt(:,ikpt2)-kpt(:,ikpt)+dk(:))
       if (sum(abs(diffk(:)-nint(diffk(:))))<3*tol8)then
         ikpt_dk(1,ikpt)=ikpt2
         if(sum(diffk(:))>=3*tol8) shift_g_2(ikpt,ikpt2) = .true.
       end if
       if (sum(abs(diffk2(:)-nint(diffk2(:))))<3*tol8)then
         ikpt_dk(2,ikpt)=ikpt2
         if(sum(diffk2(:))>=3*tol8) shift_g_2(ikpt,ikpt2) = .true.
       end if
     end do
   end do

   write(message,'(a,a,a,3f9.5,a,a,3f9.5,a)')ch10,&
&   ' Computing the derivative for reciprocal vector:',ch10,&
&   dk(:),' (in reduced coordinates)',ch10,&
&   gpard(1:3),' (in cartesian coordinates - atomic units)'
   call wrtout(ab_out,message,'COLL')
   call wrtout(std_out,message,'COLL')

   if(nsppol==1)then
     write(message, '(a,i5,a,i5)')&
&     ' From band number',bdberry(1),'  to band number',bdberry(2)
   else
     write(message, '(a,i5,a,i5,a,a,a,i5,a,i5,a)')&
&     ' From band number',bdberry(1),'  to band number',bdberry(2),' for spin up,',&
&     ch10,&
&     ' from band number',bdberry(3),'  to band number',bdberry(4),' for spin down.'
   end if
   call wrtout(ab_out,message,'COLL')
   call wrtout(std_out,message,'COLL')

!  *****************************************************************************
   ABI_ALLOCATE(dudk,(2,mpw*nspinor*mband*nsppol))
   ABI_ALLOCATE(eig_dum_2,((2*mband)**formeig*mband))
   ABI_ALLOCATE(occ_dum_2,((2*mband)**formeig*mband))
   dudk(1:2,:)=0.0_dp
   eig_dum_2=0.0_dp
   occ_dum_2=0.0_dp

   if (mkmem/=0) then

!    Find the location of each wavefunction

     ABI_ALLOCATE(cg_index,(mband,nkpt_,nsppol))
     icg = 0
     do isppol=1,nsppol
       do ikpt=1,nkpt_
         nband_k=nband(ikpt+(isppol-1)*nkpt_)
         npw_k=npwarr(ikpt)
         do iband=1,nband_k
           cg_index(iband,ikpt,isppol)=(iband-1)*npw_k*nspinor+icg
         end do
         icg=icg+npw_k*nspinor*nband(ikpt)
       end do
     end do

!    Find the planewave vectors for each k point
!    SHOULD BE REMOVED WHEN ANOTHER INDEXING TECHNIQUE WILL BE USED FOR kg
     ABI_ALLOCATE(kg_kpt,(3,mpw*nspinor,nkpt_))
     kg_kpt(:,:,:) = 0
     index1 = 0
     do ikpt=1,nkpt_
       npw_k=npwarr(ikpt)
       do ipw=1,npw_k*nspinor
         kg_kpt(1:3,ipw,ikpt)=kg(1:3,ipw+index1)
       end do
       index1=index1+npw_k*nspinor
     end do
   end if

!  *************************************************************************
!  Loop over spins
   do isppol=1,nsppol

     minband=bdberry(2*isppol-1)
     maxband=bdberry(2*isppol)

     if(minband<1)then
       write(message,'(a,i0,a)')'  The band limit minband=',minband,', is lower than 0.'
       MSG_BUG(message)
     end if

     if(maxband<1)then
       write(message,'(a,i0,a)')' The band limit maxband=',maxband,', is lower than 0.'
       MSG_BUG(message)
     end if

     if(maxband<minband)then
       write(message,'(a,i0,a,i0)')&
&       '  maxband=',maxband,', is lower than minband=',minband
       MSG_BUG(message)
     end if

!    Loop over k points
     do ikpt_=1,nkpt_

       read_k = 0

       ikpt=ikpt_
       tr(1) = 1.0_dp

       if (kptopt==2) then
         if (read_k == 0) then
           if (kpt_flag(ikpt_)/=0) then
             tr(1) = -1.0_dp
             ikpt= kpt_flag(ikpt_)
           end if
         else           !read_k
           if (kpt_flag(ikpt_)/=0) then
             tr(-1*read_k+3) = -1.0_dp
             ikpt= kpt_flag(ikpt_)
           end if
         end if       !read_k
       end if           !kptopt

       nband_k=nband(ikpt+(isppol-1)*nkpt_)

       if(nband_k<maxband)then
         write(message,'(a,i0,a,i0)')'  maxband=',maxband,', is larger than nband(i,isppol)=',nband_k
         MSG_BUG(message)
       end if

       npw_k=npwarr(ikpt)

       ABI_ALLOCATE(u_tilde,(2,npw_k,maxband,2))
       u_tilde(1:2,1:npw_k,1:maxband,1:2)=0.0_dp

!      ifor = 1,2 represents forward and backward neighbouring k points of ikpt
!      respectively along dk direction

       do ifor=1,2

         ABI_ALLOCATE(phi,(2,mpw,mband))
         ABI_ALLOCATE(cmatrix,(2,maxband,maxband))
         phi(1:2,1:mpw,1:mband)=0.0_dp; cmatrix(1:2,1:maxband,1:maxband)=0.0_dp

         isgn=(-1)**ifor
         jkpt_= ikpt_dk(ifor,ikpt_)

         tr(2) = 1.0_dp

         jkpt=jkpt_

         if (kptopt==2) then
           if (read_k == 0) then
             if (kpt_flag(jkpt_)/=0) then
               tr(2) = -1.0_dp
               jkpt= kpt_flag(jkpt_)
             end if
           else           !read_k
             if (kpt_flag(jkpt_)/=0) then
               tr(read_k) = -1.0_dp
               jkpt= kpt_flag(jkpt_)
             end if
           end if       !read_k
         end if           !kptopt

         if (ifor==1) read_k = 2

         jj = read_k
         ii = -1*read_k+3

         call waveformat(cg,cg_disk,cg_index,phi,dk,ii,ikpt,&
&         ikpt_,isgn,isppol,jj,jkpt,jkpt_,kg_kpt,kpt,kg_jl,maxband,mband,mcg,mcg_disk,&
&         minband,mkmem,mpw,nkpt,nkpt_,npwarr,nsppol,nspinor,shift_g_2,tr)

!        Compute the overlap matrix <u_k|u_k+b>

         do iband=minband,maxband
           cg_index_iband=cg_index(iband,ikpt,isppol)
           do jband=minband,maxband
             do ipw=1,npwarr(ikpt)
               cmatrix(1,iband,jband)=cmatrix(1,iband,jband)+&
&               cg(1,ipw+cg_index_iband)*phi(1,ipw,jband)+&
&               tr(ii)*cg(2,ipw+cg_index_iband)*tr(jj)*phi(2,ipw,jband)

               cmatrix(2,iband,jband)=cmatrix(2,iband,jband)+&
&               cg(1,ipw+cg_index_iband)*tr(jj)*phi(2,ipw,jband)-&
&               tr(ii)*cg(2,ipw+cg_index_iband)*phi(1,ipw,jband)
             end do
           end do
         end do

!        Compute the inverse of cmatrix(1:2,minband:maxband, minband:maxband)

         band_in = maxband - minband + 1
         ABI_ALLOCATE(ipvt,(maxband))
         ABI_ALLOCATE(zgwork,(2,1:maxband))

!        Last argument of zgedi means calculate inverse only
         call dzgefa(cmatrix(1,minband,minband),maxband, band_in,ipvt,info)
         call dzgedi(cmatrix(1,minband,minband),maxband, band_in,ipvt,det,zgwork,01)

         ABI_DEALLOCATE(zgwork)
         ABI_DEALLOCATE(ipvt)

!        Compute the product of Inverse overlap matrix with the wavefunction

         do iband=minband,maxband
           do ipw=1,npwarr(ikpt)
             u_tilde(1,ipw,iband,ifor)= &
&             dot_product(cmatrix(1,minband:maxband,iband),&
&             phi(1,ipw,minband:maxband))-&
&             dot_product(cmatrix(2,minband:maxband,iband),&
&             tr(jj)*phi(2,ipw,minband:maxband))
             u_tilde(2,ipw,iband,ifor)= &
&             dot_product(cmatrix(1,minband:maxband,iband),&
&             tr(jj)*phi(2,ipw,minband:maxband))+&
&             dot_product(cmatrix(2,minband:maxband,iband),&
&             phi(1,ipw,minband:maxband))
           end do
         end do
         ABI_DEALLOCATE(cmatrix)
         ABI_DEALLOCATE(phi)

       end do !ifor

!      Compute dudk for ikpt

       npw_k=npwarr(ikpt)

       do iband=minband,maxband

         icg=(iband-minband)*npw_k

         dudk(1,1+icg:npw_k+icg)=(u_tilde(1,1:npw_k,iband,1)-&
&         u_tilde(1,1:npw_k,iband,2))/twodk

         dudk(2,1+icg:npw_k+icg)=(u_tilde(2,1:npw_k,iband,1)-&
&         u_tilde(2,1:npw_k,iband,2))/twodk

       end do

       tim_rwwf=0
       mcg_disk=mpw*nspinor*mband
       nband_diff=maxband-minband+1
       call rwwf(dudk,eig_dum_2,formeig,0,0,ikpt,isppol,kg_kpt(:,:,ikpt),&
&       mband,mcg_disk,mpi_enreg,nband_diff,nband_diff,&
&       npw_k,nspinor,occ_dum_2,2,1,tim_rwwf,wffddk)

       ABI_DEALLOCATE(u_tilde)

     end do !ikpt
   end do  !isppol

   ABI_DEALLOCATE(eig_dum_2)
   ABI_DEALLOCATE(occ_dum_2)
   ABI_DEALLOCATE(dudk)

   call WffClose(wffddk,ierr)

   ABI_DEALLOCATE(kg_kpt)
   ABI_DEALLOCATE(cg_index)
   ABI_DEALLOCATE(ikpt_dk)

 end do ! iberry

 ABI_DEALLOCATE(shift_g_2)
 ABI_DEALLOCATE(kpt)

 write(std_out,*) 'uderiv:  exit '

end subroutine uderiv
!!***
