Entry 4705

subroutine electronic_kdotp_matrices(wvfn,dens,eigenvalues,nk,Gfrac,Bfrac)
    use constants, only: dp,cmplx_0,cmplx_i,cmplx_1,two_pi
    use comms, only: on_root_node
    use basis, only: pw_gk_vector,num_plane_waves_kp
    use cell, only: current_cell,kpoints_on_node,kpoint_weights_on_node
    use density, only: electron_density
    use ion, only: max_ps_projectors
    use pot, only : potential,pot_allocate,pot_deallocate
    use locpot, only : locpot_calculate
    use nlpot, only : nlpot_apply_velocity,nlpot_calculate_d
    use parameters, only : seedname
    use wave, only: wavefunction,band,wave_allocate,wave_initialise,wave_deallocate
    use trace, only: trace_entry,trace_exit
    use io, only: stdout,io_abort,io_atomic_to_unit,io_allocate_abort
    
    implicit none

    ! arguments
    type(wavefunction), intent(inout) :: wvfn !inout for nlpot_apply_velocity
    type(electron_density),intent(in) :: dens
    real(kind=dp),dimension(:,:,:), intent(in) :: eigenvalues
    integer, intent(in) :: nk
    complex(kind=dp),dimension(:,:,:,:),intent(out) :: Gfrac   ! Fractional Gradient vector
    complex(kind=dp),dimension(:,:,:,:,:),intent(out) :: Bfrac ! Fractional Curvature tensor


    ! local variables
    character(len=1), dimension(3) :: xyz=(/'X','Y','Z'/)
    integer :: ns,nb1,nb2,nw1,nc1,nc2  ! loop indices
    integer :: status
    type(wavefunction) :: tmp_wvfn
    complex(kind=dp),allocatable,dimension(:,:,:,:,:) :: momentum
    complex(kind=dp),allocatable,dimension(:,:,:,:) :: G       ! Gradient vector
    complex(kind=dp),allocatable,dimension(:,:,:,:,:) :: B     ! Curvature tensor
    complex(kind=dp) :: Btmp(3,3)                              ! Temporary Curvature
    complex(kind=dp),allocatable,dimension(:,:,:,:) :: kB      ! Curvature tensor premultiplied by k direction

    real(kind=dp) :: jacdet ! Fractional -> Cartesian Jacobian Determinant
    real(kind=dp) :: prodxyz ! Volume element

    real(kind=dp) :: diff_fdk(3)
    real(kind=dp) :: prod_dfdk(3)
    real(kind=dp) :: diff_fdk_sq(3)
    real(kind=dp) :: prod_dfdk_sq(3)
    real(kind=dp) :: diff_fdk_cu(3)


    ! nl x
    real(kind=dp), allocatable, dimension(:,:,:,:,:)::nl_d
    ! Potentials
    type(potential)::pot_loc     ! Up potential (or total potential for nspins=1 calculations)
    real(kind=dp) :: pot_energy

    write(*,*) "wvfn,nbands :", size(wvfn%nbands,1),size(wvfn%nbands,2), "<-- ja531"

    call trace_entry('electronic_kdotp_matrices',status)

    allocate(momentum(wvfn%nbands_max,wvfn%nbands_max,wvfn%nkpts,wvfn%nspins,3))
    allocate(G(wvfn%nbands_max,wvfn%nkpts,wvfn%nspins,3))
    allocate(B(wvfn%nbands_max,wvfn%nkpts,wvfn%nspins,3,3))
    allocate(kB(wvfn%nbands_max,wvfn%nkpts,wvfn%nspins,3))
 
    call wave_allocate(tmp_wvfn,trim(seedname)//'.kpwfn',nb=wvfn%nbands_max,nk=wvfn%nkpts,ns=wvfn%nspins)

    call wave_initialise(tmp_wvfn,'Z')

    ! ** Allocate non-local weights

    allocate(nl_d(max_ps_projectors,max_ps_projectors,current_cell%max_ions_in_species, &
         current_cell%num_species,wvfn%nspins),stat=status)
    if (status/=0) call io_allocate_abort('nl_d','electronic_kdotp_matrices')

    ! ** Allocate potentials

    call pot_allocate(pot_loc,dens%nspins,'R')

    ! ** Calculate local potential and non-local weights

    call locpot_calculate(dens,pot_loc,pot_energy)
    call nlpot_calculate_d(pot_loc,nl_d)


    momentum=cmplx_0
    G=cmplx_0
    B=cmplx_0
    kB=cmplx_0

!write(*,*) "--> size", size(Gfrac,1),size(Gfrac,2),size(Gfrac,3),size(Gfrac,4)
    ! Calculate momentum matrix
    do ns=1,wvfn%nspins
!       do nk=1,wvfn%nkpts
          do nc1=1,3
             call nlpot_apply_velocity(xyz(nc1),wvfn,nk,ns,nl_d,eigenvalue=eigenvalues(:,nk,ns),Vel_wvfn=tmp_wvfn)
             do nb1=1,wvfn%nbands_max
                do nb2=1,wvfn%nbands_max
                   do nw1=1,num_plane_waves_kp(nk)
                      momentum(nb1,nb2,nk,ns,nc1)=momentum(nb1,nb2,nk,ns,nc1)+(tmp_wvfn%coeffs(nw1,nb1,nk,ns)*&
                           &conjg(wvfn%coeffs(nw1,nb2,nk,ns)))
                   end do
                end do
             end do
          end do
!       end do
    end do

if(any(real(momentum,dp)<0.0_dp)) write(*,*) "count:", count(real(momentum,dp)<0.0_dp),&
     & count(real(momentum,dp)>0.0_dp), count(real(momentum,dp).eq.0.0_dp)
!if(any(aimag(momentum)>epsilon(1.0_dp))) write(*,*) "imagy:", maxval(imagpart(momentum))


    ! Calculate G (Gradient) vector
    G=0.0_dp
    do ns=1,wvfn%nspins
 !      do nk=1,wvfn%nkpts
          do nb1=1,wvfn%nbands_max
             do nc1=1,3
                G(nb1,nk,ns,nc1)= momentum(nb1,nb1,nk,ns,nc1)
             end do
          end do
 !      end do
    end do

    ! Transform G into fractional basis
    Gfrac=0.0_dp
    do ns=1,wvfn%nspins
 !      do nk=1,wvfn%nkpts
          do nb1=1,wvfn%nbands_max
             !do nc1=1,3
             !   Gfrac(nb1,nk,ns,nc1)= dot_product(current_cell%recip_lattice(nc1,:),G(nb1,nk,ns,:))
             !end do
             Gfrac(nb1,nk,ns,:)=matmul(current_cell%recip_lattice,G(nb1,nk,ns,:))
          end do
 !      end do
    end do

Gfrac=Gfrac/two_pi

    !B=1.0_dp
    ! set diagonal elements of B to 1.0 to begin with.
    do ns=1,wvfn%nspins
 !      do nk=1,wvfn%nkpts
          do nb1=1,wvfn%nbands_max
!             do nb2=1,wvfn%nbands_max
                do nc1=1,3
                   B(nb1,nk,ns,nc1,nc1)=+cmplx_1
                end do
 !            end do
          end do
 !      end do
    end do

    ! Calculate B (Curvature) Matrix
    do ns=1,wvfn%nspins
 !      do nk=1,wvfn%nkpts
          do nb1=1,wvfn%nbands_max
             do nb2=1,wvfn%nbands_max
                if(nb1.eq.nb2) cycle ! catch 1st divide by zero
                if(abs(eigenvalues(nb1,nk,ns)-eigenvalues(nb2,nk,ns)).le.1e-7) cycle !epsilon(1.0_dp)) cycle ! catch 2nd divide by zero
                do nc1=1,3
                   do nc2=1,3
                      B(nb1,nk,ns,nc1,nc2)=B(nb1,nk,ns,nc1,nc2)+(((momentum(nb1,nb2,nk,ns,nc1)*&
                           &conjg(momentum(nb2,nb1,nk,ns,nc2)))+(momentum(nb1,nb2,nk,ns,nc2)*&
                           &conjg(momentum(nb2,nb1,nk,ns,nc1))))/(eigenvalues(nb1,nk,ns)-&
                           &eigenvalues(nb2,nk,ns)))
                   end do
                end do
!                write(*,*) "eigthing-->",eigenvalues(nb1,nk,ns)-eigenvalues(nb2,nk,ns)
             end do
          end do
!       end do
    end do
if(any(real(B,dp)<0.0_dp)) write(*,*) "count:", count(real(B,dp)<0.0_dp), count(real(B,dp)>0.0_dp), count(real(B,dp).eq.0.0_dp)
    B=B*0.5_dp
 
    ! Transform B into fractional basis
    do ns=1,wvfn%nspins
!       do nk=1,wvfn%nkpts
          do nb1=1,wvfn%nbands_max
             do nc1=1,3
                do nc2=1,3
                   Btmp=matmul((current_cell%recip_lattice),B(nb1,nk,ns,:,:))
                   Bfrac(nb1,nk,ns,:,:)=matmul(Btmp,transpose(current_cell%recip_lattice))
!                   call dgemm("T","N",3,3,3,1.0_dp,current_cell%recip_lattice,size(current_cell%recip_lattice,1),B(nb1,nk,ns,:,:)...)
                end do
             end do
!             write(*,*) "maxB-->",maxval(abs(Bfrac(nb1,nk,ns,:,:))),maxloc(abs(Bfrac(nb1,nk,ns,:,:)))
          end do
!       end do
    end do
    
if(any(real(Bfrac,dp)<0.0_dp)) write(*,*) "countfrac:", count(real(Bfrac,dp)<0.0_dp), &
     &count(real(Bfrac,dp)>0.0_dp), count(real(Bfrac,dp).eq.0.0_dp)
!Bfrac=Bfrac/two_pi
Bfrac=Bfrac/(two_pi**2)

    deallocate(momentum,stat=status)
    if (status/=0) call io_abort('error deallocating momentum in electronic_kdotp_matrices')
    deallocate(G,stat=status)
    if (status/=0) call io_abort('error deallocating G in electronic_kdotp_matrices')
    deallocate(B,stat=status)
    if (status/=0) call io_abort('error deallocating B in electronic_kdotp_matrices')
    deallocate(kB,stat=status)
    if (status/=0) call io_abort('error deallocating kB in electronic_kdotp_matrices')

    call wave_deallocate(tmp_wvfn)

    deallocate(nl_d,stat=status)
    if (status/=0) call io_abort('error deallocating nl_d in electronic_kdotp_matrices')

    call pot_deallocate(pot_loc)

    call trace_exit('electronic_kdotp_matrices',status)

end subroutine electronic_kdotp_matrices