LCOV - code coverage report
Current view: top level - wannier/uhu - wann_uHu_sph.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   Muffin tin contribution to uHu       c
       3             : c  < u_{k+b1} | H_{k}^{mt} | u_{k+b2} >  c
       4             : c****************************************c
       5             : c   Use T(lmp,lm) calculated before for  c
       6             : c   every pair (b1,b2) and all atoms n   c
       7             : c                                        c
       8             : c   acof  ,bcof  ,ccof  : coefficients   c
       9             : c                         at k+b1        c
      10             : c   acof_b,bcof_b,ccof_b: coefficients   c
      11             : c                         at k+b2        c
      12             : c                                        c
      13             : c   bkpt   : k-point                     c
      14             : c   bkpt_b : (k+b1)-point                c
      15             : c   bkpt_b2: (k+b2)-point                c
      16             : c****************************************c
      17             : c               J.-P. Hanke, Dec. 2015   c
      18             : c****************************************c
      19             :       MODULE m_wann_uHu_sph
      20             :       contains
      21           0 :       subroutine wann_uHu_sph(
      22             :      >        chi,nbnd,llod,nslibd,nslibd_b,nlod,natd,ntypd,lmd,jmtd,
      23           0 :      >        taual,nop,lmax,
      24           0 :      >        ntype,neq,nlo,llo,acof,bcof,ccof,bkpt_b2,
      25           0 :      >        acof_b,bcof_b,ccof_b,bkpt_b,bkpt,gb_b,gb_b2,
      26           0 :      >        tuu,tud,tdu,tdd,tuulo,tulou,tdulo,tulod,tuloulo,
      27           0 :      >        kdiff,kdiff2,nntot,nntot2,uHu)
      28             : #include "cpp_double.h"
      29             :       USE m_juDFT
      30             :       use m_constants, only : pimach
      31             :       use m_matmul   , only : matmul3,matmul3r
      32             : 
      33             :       implicit none
      34             : 
      35             : c     .. scalar arguments ..
      36             :       integer, intent (in) :: llod,nlod,natd,ntypd,lmd,nbnd
      37             :       integer, intent (in) :: nntot,nntot2
      38             :       integer, intent (in) :: ntype,nslibd,nslibd_b,nop,jmtd
      39             :       complex, intent (in) :: chi
      40             : 
      41             : c     .. array arguments ..
      42             :       integer, intent (in)  :: neq(:) !(ntypd)
      43             :       integer, intent (in)  :: lmax(:) !(ntypd)
      44             :       integer, intent (in)  :: nlo(:) !(ntypd)
      45             :       integer, intent (in)  :: llo(:,:) !(nlod,ntypd)
      46             :       real,    intent (in)  :: bkpt_b(:) !(3)
      47             :       real,    intent (in)  :: bkpt_b2(:) !(3)
      48             :       real,    intent (in)  :: taual(:,:) !(3,natd)
      49             :       real,    intent (in)  :: bkpt(:) !(3)
      50             :       integer, intent (in)  :: gb_b(:) !(3)
      51             :       integer, intent (in)  :: gb_b2(:) !(3)
      52             :       complex, intent (in)  :: ccof(-llod:,:,:,:) !(-llod:llod,nslibd,nlod,natd)
      53             :       complex, intent (in)  :: acof(:,0:,:) !(nslibd,0:lmd,natd)
      54             :       complex, intent (in)  :: bcof(:,0:,:) !(nslibd,0:lmd,natd)
      55             :       complex, intent (in)  :: ccof_b(-llod:,:,:,:) !(-llod:llod,nslibd_b,nlod,natd)
      56             :       complex, intent (in)  :: acof_b(:,0:,:) !(nslibd_b,0:lmd,natd)
      57             :       complex, intent (in)  :: bcof_b(:,0:,:) !(nslibd_b,0:lmd,natd)
      58             : 
      59             :       complex, intent (in)  :: tuu(0:,0:,:,:)
      60             :       complex, intent (in)  :: tud(0:,0:,:,:)
      61             :       complex, intent (in)  :: tdu(0:,0:,:,:)
      62             :       complex, intent (in)  :: tdd(0:,0:,:,:)
      63             :       complex, intent (in)  :: tuulo(0:,:,-llod:,:,:)
      64             :       complex, intent (in)  :: tulou(0:,:,-llod:,:,:)
      65             :       complex, intent (in)  :: tdulo(0:,:,-llod:,:,:)
      66             :       complex, intent (in)  :: tulod(0:,:,-llod:,:,:)
      67             :       complex, intent (in)  :: tuloulo(:,-llod:,:,-llod:,:,:)
      68             : 
      69             :       real, intent (in)     :: kdiff(:,:),kdiff2(:,:)
      70             :       complex,intent(inout) :: uHu(:,:)
      71             : 
      72             : 
      73             : c     .. local scalars ..
      74             :       integer i,lm,nn,n,na,j,lmp,l,lp,m,mp,lwn,lo,lop
      75             :       integer ll,llp,lmd2
      76             :       real rph,cph,tpi,sqpi16,th,t1nn,t2nn,t3nn
      77             :       integer nene,nene2,indexx
      78             :       complex :: fac1,fac2,fac3,fac4
      79           0 :       complex :: mat(0:lmd,nslibd_b)
      80             : C     ..
      81             : C     .. local arrays ..
      82             :       real bpt(3)
      83             :       real bpt2(3)
      84             : 
      85             : C     ..
      86             : C     .. intrinsic functions ..
      87             :       intrinsic conjg,cmplx,sqrt,cos,sin
      88             : 
      89           0 :       tpi = 2* pimach()
      90           0 :       sqpi16 = 4*tpi*tpi
      91           0 :       lmd2 = lmd+1
      92             :       
      93           0 :       na = 0
      94             : 
      95             :       ! find neighbor k+b1
      96           0 :       bpt(:) = bkpt_b(:) + gb_b(:) - bkpt(:)
      97           0 :       do nene=1,nntot
      98           0 :          if(all(abs(bpt(:)-kdiff(:,nene)).lt.1e-4)) exit
      99             :       enddo
     100           0 :       IF(nene==nntot+1) CALL juDFT_error
     101             :      +     ("cannot find matching nearest neighbor k+b1",calledby
     102           0 :      +     ="wann_uHu_sph")
     103             : 
     104             :       ! find neighbor k+b2
     105           0 :       bpt2(:) = bkpt_b2(:) + gb_b2(:) - bkpt(:)
     106           0 :       do nene2=1,nntot2
     107           0 :          if(all(abs(bpt2(:)-kdiff2(:,nene2)).lt.1e-4)) exit
     108             :       enddo  
     109           0 :       IF(nene2==nntot2+1) CALL juDFT_error
     110             :      +     ("cannot find matching nearest neighbor k+b2",calledby
     111           0 :      +     ="wann_uHu_sph")
     112             : 
     113           0 :       indexx=nene2+(nene-1)*nntot2
     114             : 
     115             : c      if(nene2.ne.1) stop 'nene2.ne.1'
     116             : c      if(indexx.ne.nene) stop 'nene.ne.indexx'
     117             : c      if(ANY(bpt2.ne.0.0)) stop 'bpt2.ne.0'
     118             : 
     119           0 :       do n=1,ntype
     120           0 :        lwn = lmax(n)
     121           0 :          do nn = 1,neq(n) ! cycle by the atoms within the atom type
     122           0 :          na = na + 1
     123             : c...set up phase factors ( 4pi e^{ib2\tau} 4pi e^{-ib1\tau} )
     124             : 
     125           0 :          t1nn =  tpi*taual(1,na)
     126           0 :          t2nn =  tpi*taual(2,na)
     127           0 :          t3nn =  tpi*taual(3,na)
     128             : 
     129             :          th = (bpt2(1)-bpt(1))*t1nn
     130             :      >       +(bpt2(2)-bpt(2))*t2nn
     131           0 :      >       +(bpt2(3)-bpt(3))*t3nn
     132           0 :          rph = sqpi16*cos(th)
     133           0 :          cph = sqpi16*sin(th)
     134             : 
     135             : 
     136             : c...apw-apw
     137             :          call CPP_BLAS_cgemm('T','C',lmd2,nslibd_b,lmd2,cmplx(rph,cph),
     138             :      >                    tuu(0,0,n,indexx),lmd2,
     139             :      >                    acof_b(1,0,na),nslibd_b,
     140           0 :      >                    cmplx(0.0),mat(0,1),lmd2)
     141             :          call CPP_BLAS_cgemm('N','N',nslibd,nslibd_b,lmd2,chi,
     142             :      >                    acof(1,0,na),nslibd,mat(0,1),lmd2,
     143           0 :      >                    cmplx(1.0),uHu,nbnd)
     144             : 
     145             :          call CPP_BLAS_cgemm('T','C',lmd2,nslibd_b,lmd2,cmplx(rph,cph),
     146             :      >                    tud(0,0,n,indexx),lmd2,
     147             :      >                    bcof_b(1,0,na),nslibd_b,
     148           0 :      >                    cmplx(0.0),mat(0,1),lmd2)
     149             :          call CPP_BLAS_cgemm('N','N',nslibd,nslibd_b,lmd2,chi,
     150             :      >                    acof(1,0,na),nslibd,mat(0,1),lmd2,
     151           0 :      >                    cmplx(1.0),uHu,nbnd)
     152             : 
     153             :          call CPP_BLAS_cgemm('T','C',lmd2,nslibd_b,lmd2,cmplx(rph,cph),
     154             :      >                    tdu(0,0,n,indexx),lmd2,
     155             :      >                    acof_b(1,0,na),nslibd_b,
     156           0 :      >                    cmplx(0.0),mat(0,1),lmd2)
     157             :          call CPP_BLAS_cgemm('N','N',nslibd,nslibd_b,lmd2,chi,
     158             :      >                    bcof(1,0,na),nslibd,mat(0,1),lmd2,
     159           0 :      >                    cmplx(1.0),uHu,nbnd)
     160             : 
     161             :          call CPP_BLAS_cgemm('T','C',lmd2,nslibd_b,lmd2,cmplx(rph,cph),
     162             :      >                    tdd(0,0,n,indexx),lmd2,
     163             :      >                    bcof_b(1,0,na),nslibd_b,
     164           0 :      >                    cmplx(0.0),mat(0,1),lmd2)
     165             :          call CPP_BLAS_cgemm('N','N',nslibd,nslibd_b,lmd2,chi,
     166             :      >                    bcof(1,0,na),nslibd,mat(0,1),lmd2,
     167           0 :      >                    cmplx(1.0),uHu,nbnd)
     168             : 
     169           0 :          if(nlo(n).ge.1) then
     170             : 
     171             : c...apw-lo
     172           0 :          do lo = 1,nlo(n)
     173           0 :           l = llo(lo,n)
     174           0 :           do m = -l, l
     175             :            
     176           0 :            do lp = 0, lwn
     177           0 :             llp = lp*(lp+1)
     178           0 :             do mp = -lp, lp
     179           0 :              lmp = llp + mp
     180             : 
     181           0 :              fac1=cmplx(rph,cph)*tulou(lmp,lo,m,n,indexx)*chi
     182           0 :              fac2=cmplx(rph,cph)*tulod(lmp,lo,m,n,indexx)*chi
     183           0 :              fac3=cmplx(rph,cph)*tuulo(lmp,lo,m,n,indexx)*chi
     184           0 :              fac4=cmplx(rph,cph)*tdulo(lmp,lo,m,n,indexx)*chi
     185             : 
     186           0 :              do i = 1,nslibd
     187           0 :                do j = 1,nslibd_b
     188             :                  uHu(i,j) = uHu(i,j)
     189             :      >             + ccof(m,i,lo,na)* fac1 *conjg(acof_b(j,lmp,na))
     190             :      >             + ccof(m,i,lo,na)* fac2 *conjg(bcof_b(j,lmp,na))
     191             :      >             + acof(i,lmp,na) * fac3 *conjg(ccof_b(m,j,lo,na))
     192           0 :      >             + bcof(i,lmp,na) * fac4 *conjg(ccof_b(m,j,lo,na))
     193             :                enddo 
     194             :              enddo  
     195             : 
     196             :             enddo !mp
     197             :            enddo  !lp
     198             : 
     199             :           enddo ! m
     200             :          enddo  ! lo
     201             : 
     202             : c...lo-lo
     203           0 :          do lo = 1,nlo(n)
     204           0 :           l = llo(lo,n)
     205           0 :           do m = -l, l
     206             :            
     207           0 :            do lop = 1,nlo(n)
     208           0 :             lp = llo(lop,n)
     209           0 :             do mp = -lp, lp
     210             : 
     211           0 :              fac1=cmplx(rph,cph)*tuloulo(lop,mp,lo,m,n,indexx)*chi
     212             : 
     213           0 :              do i = 1,nslibd
     214           0 :                do j = 1,nslibd_b
     215             :                  uHu(i,j) = uHu(i,j)
     216           0 :      >             + ccof(m,i,lo,na)*fac1*conjg(ccof_b(mp,j,lop,na))
     217             :                enddo 
     218             :              enddo  
     219             : 
     220             :             enddo !mp
     221             :            enddo  !lop
     222             : 
     223             :           enddo ! m
     224             :          enddo  ! lo
     225             : 
     226             :          endif !(nlo(n).ge.1)
     227             : 
     228             :       enddo  ! atoms in the type
     229             : 
     230             :       enddo  ! atom type 
     231             : 
     232           0 :       end subroutine wann_uHu_sph
     233             :       end module m_wann_uHu_sph

Generated by: LCOV version 1.13