LCOV - code coverage report
Current view: top level - wannier/uhu - wann_uHu_int.F (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 66 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 1 0.0 %

          Line data    Source code
       1             : c*****************************************c
       2             : c   Interstitial contribution to uHu      c
       3             : c  < u_{k+b1} | H_{k}^{int} | u_{k+b2} >  c
       4             : c*****************************************c
       5             : c  k1_b,  k2_b,  k3_b : G-vectors at k+b1 c
       6             : c  k1_b2, k2_b2, k3_b2: G-vectors at k+b2 c
       7             : c                                         c
       8             : c  gb : reciprocal latt. vector taking    c
       9             : c       k+b1 back to first Brillouin zone c
      10             : c  gb2: reciprocal latt. vector taking    c
      11             : c       k+b2 back to first Brillouin zone c
      12             : c                                         c
      13             : c  z_b : eigenvectors at k+b1             c
      14             : c  z_b2: eigenvectors at k+b2             c
      15             : c*****************************************c
      16             : c                J.-P. Hanke, Dec. 2015   c
      17             : c*****************************************c
      18             :       module m_wann_uHu_int
      19             :       contains
      20           0 :       subroutine wann_uHu_int(
      21             :      >               chi,nvd,k1d,k2d,k3d,
      22             :      >               n3d,nv_b,nv_b2,nbnd,neigd,
      23             :      >               nslibd_b,nslibd_b2,nbasfcn,
      24             :      >               addnoco,addnoco2,
      25           0 :      >               k1_b ,k2_b ,k3_b, gb,
      26           0 :      >               k1_b2,k2_b2,k3_b2,gb2,
      27           0 :      >               bkpt,bbmat,vpw,zMat_b,zMat_b2,
      28           0 :      >               rgphs,ustep,ig,l_kin,sign,uHu)
      29             : #include "cpp_double.h"
      30             : 
      31             :       USE m_types
      32             : 
      33             :       implicit none
      34             : 
      35             :       TYPE(t_zmat), INTENT(IN) :: zMat_b, zMat_b2
      36             : 
      37             : c     ..arguments..
      38             :       logical, intent(in) :: l_kin
      39             :       integer, intent(in) :: gb(3),gb2(3)
      40             :       integer, intent(in) :: k1d,k2d,k3d
      41             :       integer, intent(in) :: nvd,n3d,nv_b,nv_b2,nbnd,neigd
      42             :       integer, intent(in) :: nslibd_b,nslibd_b2,nbasfcn
      43             :       integer, intent(in) :: addnoco,addnoco2,sign
      44             :       integer, intent(in) :: k1_b(nvd) ,k2_b(nvd) ,k3_b(nvd)
      45             :       integer, intent(in) :: k1_b2(nvd),k2_b2(nvd),k3_b2(nvd)
      46             :       integer, intent(in) :: ig(-k1d:k1d,-k2d:k2d,-k3d:k3d)
      47             :       complex, intent(in) :: rgphs(-k1d:k1d,-k2d:k2d,-k3d:k3d)
      48             :       real,    intent(in) :: bkpt(3),bbmat(3,3)
      49             :       complex, intent(in) :: ustep(n3d),vpw(n3d)
      50             :       complex, intent(in) :: chi
      51             : 
      52             :       complex, intent(inout) :: uHu(nbnd,nbnd)
      53             : 
      54             : c     ..local variables..
      55           0 :       complex,allocatable :: vstep_c(:,:)
      56           0 :       complex,allocatable :: mat_c(:,:)
      57             :       complex :: th_c
      58           0 :       real,allocatable :: vstep_r(:,:)
      59           0 :       real,allocatable :: mat_r(:,:)
      60           0 :       real,allocatable :: uHu_tmp(:,:)
      61             :       real :: th_r
      62             : 
      63           0 :       real,allocatable :: rk_b(:),rk_b2(:)
      64             :       real :: phase,phase2,ekk,s(3)
      65             :       integer :: i,j,j1,j2,j3,i1,i2,i3,in,ind
      66             : 
      67           0 :       allocate( rk_b (nv_b)  )
      68           0 :       allocate( rk_b2(nv_b2) )
      69             : 
      70           0 :       IF(zMat_b%l_real) THEN
      71           0 :          allocate( uHu_tmp(nslibd_b,nslibd_b2) )
      72           0 :          allocate( vstep_r(nv_b2,nv_b) )
      73           0 :          allocate( mat_r(nv_b,nslibd_b2) )
      74           0 :          vstep_r(:,:) = 0.0
      75             :       ELSE
      76           0 :          allocate( vstep_c(nv_b2,nv_b) )
      77           0 :          allocate( mat_c(nv_b,nslibd_b2) )
      78           0 :          vstep_c(:,:) = CMPLX(0.0,0.0)
      79             :       END IF
      80             : 
      81             :       ! set up |k+G-G(k+b1)|^2
      82           0 :       do i=1,nv_b
      83           0 :        s(1) = bkpt(1) + k1_b(i) - gb(1)
      84           0 :        s(2) = bkpt(2) + k2_b(i) - gb(2)
      85           0 :        s(3) = bkpt(3) + k3_b(i) - gb(3)
      86           0 :        rk_b(i) = dot_product(s,matmul(bbmat,s))
      87             : !       rk_b(i) = dotirp(s,s,bbmat)
      88             :       enddo
      89             : 
      90             :       ! set up |k+G'-G(k+b2)|^2
      91           0 :       do i=1,nv_b2
      92           0 :        s(1) = bkpt(1) + k1_b2(i) - gb2(1)
      93           0 :        s(2) = bkpt(2) + k2_b2(i) - gb2(2)
      94           0 :        s(3) = bkpt(3) + k3_b2(i) - gb2(3)
      95           0 :        rk_b2(i) = dot_product(s,matmul(bbmat,s))
      96             : !       rk_b2(i) = dotirp(s,s,bbmat)
      97             :       enddo
      98             : 
      99             :       ! construct vstep(g,g') ~ V(g-g') 
     100             :       !                       + Theta(g-g')*[rk_b+rk_b2]
     101           0 :       do i=1,nv_b
     102           0 :        j1 = -k1_b(i) + gb(1) - gb2(1)
     103           0 :        j2 = -k2_b(i) + gb(2) - gb2(2)
     104           0 :        j3 = -k3_b(i) + gb(3) - gb2(3)
     105           0 :        do j=1,nv_b2
     106           0 :         i1 = j1 + k1_b2(j)
     107           0 :         i2 = j2 + k2_b2(j)
     108           0 :         i3 = j3 + k3_b2(j)
     109           0 :         in = ig(sign*i1,sign*i2,sign*i3)
     110           0 :         if(in.eq.0) cycle
     111           0 :         phase = rgphs(i1,i2,i3)    ! TODO: sign also here?
     112           0 :         phase2= rgphs(sign*i1,sign*i2,sign*i3)
     113           0 :         if(phase.ne.phase2) then
     114           0 :            stop 'rgphs in wann_uHu_int'
     115             :         endif
     116             : 
     117           0 :         ekk = rk_b(i) + rk_b2(j)
     118             : 
     119           0 :         IF(zMat_b%l_real) THEN
     120           0 :            th_r = phase*real(vpw(in))
     121           0 :            if(l_kin) th_r = th_r + phase*0.25*ekk*real(ustep(in))
     122           0 :            vstep_r(j,i) = th_r
     123             :         ELSE
     124           0 :            th_c = phase*conjg(vpw(in))
     125           0 :            if(l_kin) th_c = th_c + phase*0.25*ekk*conjg(ustep(in))
     126           0 :            vstep_c(j,i) = th_c
     127             :         END IF
     128             :        enddo
     129             :       enddo
     130             : 
     131             :       ! complex conjugate of (z(k+b1,g))^* vstep(g,g') z(k+b2,g')
     132           0 :       IF(zMat_b%l_real) THEN
     133             :          call CPP_BLAS_sgemm('T','N',nv_b,nslibd_b2,nv_b2,real(1.0),
     134             :      >               vstep_r,nv_b2,zMat_b2%z_r(1+addnoco2,1),nbasfcn,
     135           0 :      >               real(0.0),mat_r,nv_b)
     136             :          call CPP_BLAS_sgemm('T','N',nslibd_b,nslibd_b2,nv_b,
     137             :      >              real(1.0),zMat_b%z_r(1+addnoco,1),nbasfcn,
     138           0 :      >              mat_r,nv_b,real(0.0),uHu_tmp,nslibd_b)
     139             :          uHu(1:nslibd_b,1:nslibd_b2) = uHu(1:nslibd_b,1:nslibd_b2)
     140           0 :      >                            + uHu_tmp(1:nslibd_b,1:nslibd_b2)*chi
     141             :       ELSE
     142             :          call CPP_BLAS_cgemm('T','N',nv_b,nslibd_b2,nv_b2,cmplx(1.0),
     143             :      >               vstep_c,nv_b2,zMat_b2%z_c(1+addnoco2,1),nbasfcn,
     144           0 :      >               cmplx(0.0),mat_c,nv_b)
     145           0 :          mat_c = conjg(mat_c)
     146             :          call CPP_BLAS_cgemm('T','N',nslibd_b,nslibd_b2,nv_b,
     147             :      >              chi,zMat_b%z_c(1+addnoco,1),nbasfcn,
     148           0 :      >              mat_c,nv_b,cmplx(1.0),uHu,nbnd)
     149             :       END IF
     150             : 
     151           0 :       deallocate( rk_b, rk_b2 )
     152             : 
     153           0 :       IF(zMat_b%l_real) THEN
     154           0 :          deallocate( uHu_tmp )
     155           0 :          deallocate( vstep_r, mat_r )
     156             :       ELSE
     157           0 :          deallocate( vstep_c, mat_c )
     158             :       END IF
     159             : 
     160           0 :       end subroutine
     161             :       end module m_wann_uHu_int

Generated by: LCOV version 1.13