wann_uHu_int.F 5.29 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152
c*****************************************c
c   Interstitial contribution to uHu      c
c  < u_{k+b1} | H_{k}^{int} | u_{k+b2} >  c
c*****************************************c
c  k1_b,  k2_b,  k3_b : G-vectors at k+b1 c
c  k1_b2, k2_b2, k3_b2: G-vectors at k+b2 c
c                                         c
c  gb : reciprocal latt. vector taking    c
c       k+b1 back to first Brillouin zone c
c  gb2: reciprocal latt. vector taking    c
c       k+b2 back to first Brillouin zone c
c                                         c
c  z_b : eigenvectors at k+b1             c
c  z_b2: eigenvectors at k+b2             c
c*****************************************c
c                J.-P. Hanke, Dec. 2015   c
c*****************************************c
      module m_wann_uHu_int
      contains
      subroutine wann_uHu_int(
     >               chi,nvd,k1d,k2d,k3d,
     >               n3d,nv_b,nv_b2,nbnd,neigd,
     >               nslibd_b,nslibd_b2,nbasfcn,
     >               addnoco,addnoco2,
     >               k1_b ,k2_b ,k3_b,  gb,
     >               k1_b2,k2_b2,k3_b2,gb2,
     >               bkpt,bbmat,vpw,z_b,z_b2,
     >               rgphs,ustep,ig,l_kin,sign,uHu)
#include "cpp_double.h"
      use m_dotir, ONLY : dotirp
      implicit none

c     ..arguments..
      logical, intent(in) :: l_kin
      integer, intent(in) :: gb(3),gb2(3)
      integer, intent(in) :: k1d,k2d,k3d
      integer, intent(in) :: nvd,n3d,nv_b,nv_b2,nbnd,neigd
      integer, intent(in) :: nslibd_b,nslibd_b2,nbasfcn
      integer, intent(in) :: addnoco,addnoco2,sign
      integer, intent(in) :: k1_b(nvd) ,k2_b(nvd) ,k3_b(nvd)
      integer, intent(in) :: k1_b2(nvd),k2_b2(nvd),k3_b2(nvd)
      integer, intent(in) :: ig(-k1d:k1d,-k2d:k2d,-k3d:k3d)
      real,    intent(in) :: rgphs(-k1d:k1d,-k2d:k2d,-k3d:k3d)
      real,    intent(in) :: bkpt(3),bbmat(3,3)
      complex, intent(in) :: ustep(n3d),vpw(n3d)
      complex, intent(in) :: chi
#if ( !defined(CPP_INVERSION) || defined(CPP_SOC) )
      complex, intent(in) :: z_b(nbasfcn,neigd),z_b2(nbasfcn,neigd)
#else
      real,    intent(in) :: z_b(nbasfcn,neigd),z_b2(nbasfcn,neigd)
#endif 
      complex, intent(inout) :: uHu(nbnd,nbnd)

c     ..local variables..
#if (!defined(CPP_INVERSION)||defined(CPP_SOC))
      complex,allocatable :: vstep(:,:)
      complex,allocatable :: mat(:,:)
      complex :: th
#else
      real,allocatable :: vstep(:,:)
      real,allocatable :: mat(:,:)
      real,allocatable :: uHu_tmp(:,:)
      real :: th
#endif
      real,allocatable :: rk_b(:),rk_b2(:)
      real :: phase,phase2,ekk,s(3)
      integer :: i,j,j1,j2,j3,i1,i2,i3,in,ind

      allocate( vstep(nv_b2,nv_b) )
      allocate( mat(nv_b,nslibd_b2) )
      allocate( rk_b (nv_b)  )
      allocate( rk_b2(nv_b2) )
#if (!defined(CPP_INVERSION)||defined(CPP_SOC))
#else
      allocate( uHu_tmp(nslibd_b,nslibd_b2) )
#endif

      ! set up |k+G-G(k+b1)|^2
      do i=1,nv_b
       s(1) = bkpt(1) + k1_b(i) - gb(1)
       s(2) = bkpt(2) + k2_b(i) - gb(2)
       s(3) = bkpt(3) + k3_b(i) - gb(3)
       rk_b(i) = dotirp(s,s,bbmat)
      enddo

      ! set up |k+G'-G(k+b2)|^2
      do i=1,nv_b2
       s(1) = bkpt(1) + k1_b2(i) - gb2(1)
       s(2) = bkpt(2) + k2_b2(i) - gb2(2)
       s(3) = bkpt(3) + k3_b2(i) - gb2(3)
       rk_b2(i) = dotirp(s,s,bbmat)
      enddo

      ! construct vstep(g,g') ~ V(g-g') 
      !                       + Theta(g-g')*[rk_b+rk_b2]
      vstep(:,:) = 0.0
      do i=1,nv_b
       j1 = -k1_b(i) + gb(1) - gb2(1)
       j2 = -k2_b(i) + gb(2) - gb2(2)
       j3 = -k3_b(i) + gb(3) - gb2(3)
       do j=1,nv_b2
        i1 = j1 + k1_b2(j)
        i2 = j2 + k2_b2(j)
        i3 = j3 + k3_b2(j)
        in = ig(sign*i1,sign*i2,sign*i3)
        if(in.eq.0) cycle
        phase = rgphs(i1,i2,i3)    ! TODO: sign also here?
        phase2= rgphs(sign*i1,sign*i2,sign*i3)
        if(phase.ne.phase2) then
           stop 'rgphs in wann_uHu_int'
        endif

        ekk = rk_b(i) + rk_b2(j)
#if ( !defined(CPP_INVERSION) || defined(CPP_SOC) )        
        th = phase*conjg(vpw(in))
        if(l_kin) th = th + phase*0.25*ekk*conjg(ustep(in))
#else
        th = phase*real(vpw(in))
        if(l_kin) th = th + phase*0.25*ekk*real(ustep(in))
#endif
        vstep(j,i) = th
       enddo
      enddo

      ! complex conjugate of (z(k+b1,g))^* vstep(g,g') z(k+b2,g')
#if (!defined(CPP_INVERSION) || defined(CPP_SOC))
      call CPP_BLAS_cgemm('T','N',nv_b,nslibd_b2,nv_b2,cmplx(1.0),
     >            vstep,nv_b2,z_b2(1+addnoco2,1),nbasfcn,
     >            cmplx(0.0),mat,nv_b)
      mat = conjg(mat)
      call CPP_BLAS_cgemm('T','N',nslibd_b,nslibd_b2,nv_b,
     >           chi,z_b(1+addnoco,1),nbasfcn,
     >           mat,nv_b,cmplx(1.0),uHu,nbnd)
#else
      call CPP_BLAS_sgemm('T','N',nv_b,nslibd_b2,nv_b2,real(1.0),
     >            vstep,nv_b2,z_b2(1+addnoco2,1),nbasfcn,
     >            real(0.0),mat,nv_b)
      call CPP_BLAS_sgemm('T','N',nslibd_b,nslibd_b2,nv_b,
     >           real(1.0),z_b(1+addnoco,1),nbasfcn,
     >           mat,nv_b,real(0.0),uHu_tmp,nslibd_b)
      uHu(1:nslibd_b,1:nslibd_b2) = uHu(1:nslibd_b,1:nslibd_b2)
     >                            + uHu_tmp(1:nslibd_b,1:nslibd_b2)*chi
#endif

      deallocate( vstep, mat )
      deallocate( rk_b, rk_b2 )
#if (!defined(CPP_INVERSION) || defined(CPP_SOC))
#else
      deallocate( uHu_tmp )
#endif
      end subroutine
      end module m_wann_uHu_int