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

          Line data    Source code
       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3             : ! This file is part of FLEUR and available as free software under the conditions
       4             : ! of the MIT license as expressed in the LICENSE file in more detail.
       5             : !--------------------------------------------------------------------------------
       6             : 
       7             :       module m_wann_ujugaunt
       8             :       contains     
       9           0 :       SUBROUTINE wann_ujugaunt(
      10             : c***********************************************************************
      11             : c    wann_ujugaunt calculates integrals of radial wave functions with
      12             : c    bessel functions and multiplies them with an angular factor.
      13             : c    Calculating them only once gives some speed-up of wann_mmkb_sph.
      14             : c    Frank Freimuth, October 2006
      15             : c*********************************************************************** 
      16           0 :      >                          llod,nntot,kdiff,lmax,
      17           0 :      >                          ntype,ntypd,bbmat,bmat,
      18           0 :      >                          nlod,nlo,llo,flo,flo_b,f,f_b,g,g_b,
      19           0 :      >                          jri,rmsh,dx,jmtd,
      20             :      >                          lmaxd,lmd,
      21           0 :      <                          ujug,ujdg,djug,djdg,ujulog,djulog,
      22           0 :      <                          ulojug,ulojdg,ulojulog,l_q,sign_q)
      23             :       use m_constants, only : pimach
      24             :       use m_matmul   , only : matmul3,matmul3r
      25             :       use m_sphbes
      26             :       use m_ylm
      27             :       use m_types, only : od_inp, od_sym
      28             :       use m_intgr, only : intgr3
      29             :       use m_gaunt, only: gaunt1
      30             : 
      31             :       IMPLICIT NONE
      32             :       integer, intent (in)  :: llod
      33             :       INTEGER, INTENT (IN)  :: nntot,ntype,ntypd
      34             :       INTEGER, INTENT (IN)  :: lmaxd,jmtd,lmd
      35             :       real,    intent (in)  :: kdiff(:,:) !(3,nntot)
      36             :       real,    intent (in)  :: bbmat(:,:) !(3,3)
      37             :       real,    intent (in)  :: bmat(:,:) !(3,3)
      38             :       integer, intent (in)  :: lmax(:) !(ntypd)
      39             :       integer, intent (in)  :: nlod
      40             :       integer, intent (in)  :: jri(:) !(ntypd)
      41             :       integer, intent (in)  :: nlo(:) !(ntypd)
      42             :       integer, intent (in)  :: llo(:,:) !(nlod,ntypd)    
      43             :       real,    intent (in)  :: f(:,:,:,0:) !(ntypd,jmtd,2,0:lmaxd)
      44             :       real,    intent (in)  :: f_b(:,:,:,0:) !(ntypd,jmtd,2,0:lmaxd)
      45             :       real,    intent (in)  :: g(:,:,:,0:) !(ntypd,jmtd,2,0:lmaxd)
      46             :       real,    intent (in)  :: g_b(:,:,:,0:) !(ntypd,jmtd,2,0:lmaxd)
      47             :       real,    intent (in)  :: flo(:,:,:,:) !(ntypd,jmtd,2,nlod)
      48             :       real,    intent (in)  :: flo_b(:,:,:,:) !(ntypd,jmtd,2,nlod)
      49             :       real,    intent (in)  :: rmsh(:,:) !(jmtd,ntypd)
      50             :       real,    intent (in)  :: dx(:) !(ntypd)
      51             : 
      52             :       logical,    intent (in)  :: l_q    ! if true, we deal with q points
      53             :       integer,    intent (in)  :: sign_q  ! if we deal with q points, we might pick up an additional sign
      54             : 
      55             :       complex, intent (out) :: ujug(0:,0:,:,:) !(0:lmd,0:lmd,1:ntype,1:nntot)
      56             :       complex, intent (out) :: ujdg(0:,0:,:,:) !(0:lmd,0:lmd,1:ntype,1:nntot)
      57             :       complex, intent (out) :: djug(0:,0:,:,:) !(0:lmd,0:lmd,1:ntype,1:nntot)
      58             :       complex, intent (out) :: djdg(0:,0:,:,:) !(0:lmd,0:lmd,1:ntype,1:nntot)
      59             :       complex, intent (out) :: ujulog(0:,:,-llod:,:,:) !(0:lmd,nlod,-llod:llod,1:ntype,1:nntot)
      60             :       complex, intent (out) :: djulog(0:,:,-llod:,:,:) !(0:lmd,nlod,-llod:llod,1:ntype,1:nntot)
      61             :       complex, intent (out) :: ulojug(0:,:,-llod:,:,:) !(0:lmd,nlod,-llod:llod,1:ntype,1:nntot)
      62             :       complex, intent (out) :: ulojdg(0:,:,-llod:,:,:) !(0:lmd,nlod,-llod:llod,1:ntype,1:nntot)
      63             :       complex, intent (out) :: ulojulog(:,-llod:,:,-llod:,:,:) !(1:nlod,-llod:llod,1:nlod,-llod:llod,1:ntype,1:nntot)
      64             : 
      65           0 :       real, allocatable :: djd(:,:,:),ujd(:,:,:),uju(:,:,:)
      66           0 :       real, allocatable :: dju(:,:,:)
      67           0 :       real, allocatable :: ujulo(:,:,:),djulo(:,:,:),ulojulo(:,:,:)
      68           0 :       real, allocatable :: uloju(:,:,:),ulojd(:,:,:)
      69             :       integer           :: ikpt_b,i,lwn,n,lpp,lop,lo,l,lp
      70             :       integer           :: lmini,lmaxi,m,mp,llpp,mpp
      71             :       integer           :: lmpp,lminp,lmaxp,lm,lpmp
      72           0 :       real              :: rk,bpt(3),gs,jlpp(0:lmaxd)
      73           0 :       real              :: jj(0:lmaxd,jmtd),x(jmtd)
      74             :       real              :: bkrot(3)
      75           0 :       complex           :: ylmpp((lmaxd+1)**2),factor,ic
      76             :       complex           :: factor2
      77             : 
      78           0 :       ic = cmplx(0.,1.)
      79             : 
      80           0 :       allocate( djd(0:lmaxd,0:lmaxd,0:lmaxd) )
      81           0 :       allocate( ujd(0:lmaxd,0:lmaxd,0:lmaxd) )
      82           0 :       allocate( dju(0:lmaxd,0:lmaxd,0:lmaxd) )
      83           0 :       allocate( uju(0:lmaxd,0:lmaxd,0:lmaxd) ) 
      84             : 
      85           0 :       allocate( ujulo(nlod,0:lmaxd,0:lmaxd) )
      86           0 :       allocate( djulo(nlod,0:lmaxd,0:lmaxd) )
      87           0 :       allocate( uloju(nlod,0:lmaxd,0:lmaxd) )
      88           0 :       allocate( ulojd(nlod,0:lmaxd,0:lmaxd) )
      89             : 
      90           0 :       allocate( ulojulo(nlod,nlod,0:lmaxd) )
      91             : 
      92           0 :       ujug = cmplx(0.0,0.0)
      93           0 :       ujdg = cmplx(0.0,0.0)
      94           0 :       djug = cmplx(0.0,0.0)
      95           0 :       djdg = cmplx(0.0,0.0)
      96             : 
      97           0 :       ujulog = cmplx(0.0,0.0)
      98           0 :       djulog = cmplx(0.0,0.0)
      99           0 :       ulojug = cmplx(0.0,0.0)
     100           0 :       ulojdg = cmplx(0.0,0.0)
     101             : 
     102           0 :       ulojulog = cmplx(0.0,0.0)
     103             :       
     104           0 :       do ikpt_b=1,nntot
     105           0 :         bpt(:)=kdiff(:,ikpt_b)
     106           0 :         rk = sqrt(dot_product(bpt,matmul(bbmat,bpt)))
     107             :         !write(*,*)'ujugaunt rk',rk
     108             : 
     109           0 :         do n=1,ntype
     110           0 :          lwn = lmax(n)
     111             : c...generate the j_lpp(br) on the radial grid   
     112           0 :          do i = 1,jri(n)
     113           0 :            gs = rk*rmsh(i,n)
     114           0 :            call sphbes(lwn,gs,jlpp)
     115           0 :            jj(:,i) = jlpp(:)
     116             :          enddo
     117           0 :          do lpp = 0,lwn   ! lpp is the ang. momentum of the bessel function
     118             : c***************************************************************************
     119             : c...the local orbitals overlaps
     120             : c***************************************************************************
     121           0 :           if (nlo(n).GE.1) then
     122           0 :            do lop = 1,nlo(n)
     123           0 :             do lo = 1,nlo(n)
     124           0 :              l = llo(lo,n)
     125           0 :              lp = llo(lop,n)
     126           0 :              lmini = abs(lp - l)
     127           0 :              lmaxi = lp + l
     128             : c..the gaunt conditions
     129           0 :              if ((mod(l+lp+lpp,2).eq.1) .or. (lpp.LT.lmini) .or.
     130           0 :      +             (lpp.gt.lmaxi)) then
     131           0 :                ulojulo(lo,lop,lpp) = 0. 
     132             :              else 
     133           0 :               do i = 1,jri(n)
     134             :                 x(i) = ( flo(n,i,1,lo)*flo_b(n,i,1,lop)+
     135           0 :      +                   flo(n,i,2,lo)*flo_b(n,i,2,lop) )*jj(lpp,i)
     136             :               enddo 
     137           0 :               call intgr3(x,rmsh(1:,n),dx(n),jri(n),ulojulo(lo,lop,lpp))
     138             :              endif
     139             :             enddo
     140             :            enddo
     141             :           endif ! local orbitals 
     142             : c**************************************************************************
     143             : c...overlaps of the apws only
     144             : c**************************************************************************
     145           0 :           do lp = 0,lwn
     146           0 :            do l = 0,lwn
     147           0 :             lmini = abs(lp-l)
     148           0 :             lmaxi = lp + l
     149             : c..gaunt conditions
     150           0 :             if ((mod(l+lp+lpp,2).eq.1) .or. (lpp.LT.lmini) .or.
     151           0 :      +             (lpp.gt.lmaxi)) then
     152           0 :              uju(l,lp,lpp) = 0.
     153           0 :              ujd(l,lp,lpp) = 0.
     154           0 :              dju(l,lp,lpp) = 0.
     155           0 :              djd(l,lp,lpp) = 0.
     156             :             else
     157           0 :              do i = 1,jri(n)
     158             :                 x(i) = ( f(n,i,1,l)*f_b(n,i,1,lp)+
     159           0 :      +                   f(n,i,2,l)*f_b(n,i,2,lp) )*jj(lpp,i)
     160             :              enddo      
     161           0 :              call intgr3(x,rmsh(1:,n),dx(n),jri(n),uju(l,lp,lpp))
     162             : 
     163           0 :              do i = 1,jri(n)
     164             :                 x(i) = ( f(n,i,1,l)*g_b(n,i,1,lp)+
     165           0 :      +                   f(n,i,2,l)*g_b(n,i,2,lp) )*jj(lpp,i)
     166             :              enddo      
     167           0 :              call intgr3(x,rmsh(1:,n),dx(n),jri(n),ujd(l,lp,lpp))
     168             : 
     169           0 :              do i = 1,jri(n)
     170             :                 x(i) = ( g(n,i,1,l)*f_b(n,i,1,lp)+
     171           0 :      +                   g(n,i,2,l)*f_b(n,i,2,lp) )*jj(lpp,i)
     172             :              enddo      
     173           0 :              call intgr3(x,rmsh(1:,n),dx(n),jri(n),dju(l,lp,lpp))
     174             : 
     175           0 :              do i = 1,jri(n)
     176             :                 x(i) = ( g(n,i,1,l)*g_b(n,i,1,lp)+
     177           0 :      +                   g(n,i,2,l)*g_b(n,i,2,lp) )*jj(lpp,i)
     178             :              enddo     
     179           0 :              call intgr3(x,rmsh(1:,n),dx(n),jri(n),djd(l,lp,lpp))
     180             :             endif
     181             :            enddo ! l
     182             : 
     183             : c********************************************************************
     184             : c...overlaps of the lo's with the apws 
     185             : c********************************************************************
     186           0 :            if (nlo(n).GE.1) then
     187           0 :             do lo = 1,nlo(n)
     188           0 :              l = llo(lo,n)
     189           0 :              lmini = abs(lp-l)
     190           0 :              lmaxi = lp + l
     191             : c..gaunt conditions
     192           0 :              if ((mod(l+lp+lpp,2).eq.1) .OR. (lpp.lt.lmini) .or.
     193           0 :      +             (lpp.gt.lmaxi)) then
     194           0 :                ujulo(lo,lp,lpp) = 0.
     195           0 :                djulo(lo,lp,lpp) = 0.
     196           0 :                uloju(lo,lp,lpp) = 0.
     197           0 :                ulojd(lo,lp,lpp) = 0.
     198             :              else
     199           0 :               do i = 1,jri(n)
     200             :                x(i) = ( flo(n,i,1,lo)*f_b(n,i,1,lp)+
     201           0 :      +                  flo(n,i,2,lo)*f_b(n,i,2,lp) )*jj(lpp,i)
     202             :               enddo 
     203           0 :               call intgr3(x,rmsh(1:,n),dx(n),jri(n),ujulo(lo,lp,lpp))
     204           0 :               do i = 1,jri(n)
     205             :                x(i) = ( flo(n,i,1,lo)*g_b(n,i,1,lp)+
     206           0 :      +                  flo(n,i,2,lo)*g_b(n,i,2,lp) )*jj(lpp,i)
     207             :               enddo 
     208           0 :               call intgr3(x,rmsh(1:,n),dx(n),jri(n),djulo(lo,lp,lpp))
     209           0 :               do i = 1,jri(n)
     210             :                x(i) = ( flo_b(n,i,1,lo)*f(n,i,1,lp)+
     211           0 :      +                  flo_b(n,i,2,lo)*f(n,i,2,lp) )*jj(lpp,i)
     212             :               enddo 
     213           0 :               call intgr3(x,rmsh(1:,n),dx(n),jri(n),uloju(lo,lp,lpp))
     214           0 :               do i = 1,jri(n)
     215             :                x(i) = ( flo_b(n,i,1,lo)*g(n,i,1,lp)+
     216           0 :      +                  flo_b(n,i,2,lo)*g(n,i,2,lp) )*jj(lpp,i)
     217             :               enddo 
     218           0 :               call intgr3(x,rmsh(1:,n),dx(n),jri(n),ulojd(lo,lp,lpp))
     219             :              endif
     220             :             enddo !lo  
     221             :            endif  ! local orbitals  
     222             :           enddo !lp
     223             :          enddo !lpp
     224             : c********************************************************************
     225             : c       multiply with gaunt-coefficient (apw-apw)
     226             : c********************************************************************
     227           0 :          bkrot=matmul(bpt,bmat)
     228           0 :          call ylm4(lwn,bkrot,ylmpp)
     229           0 :          do l = 0,lwn
     230           0 :           do m = -l,l
     231           0 :            lm=l*(l+1)+m  
     232           0 :            do lp = 0,lwn
     233           0 :             do mp = -lp,lp
     234           0 :              lpmp=lp*(lp+1)+mp  
     235           0 :              do lpp = 0,lwn
     236           0 :                llpp = lpp*(lpp+1)
     237             : 
     238           0 :                mpp = mp - m
     239             : 
     240           0 :                lmpp = llpp + mpp 
     241           0 :                lmini = abs(l-lpp)
     242           0 :                lmaxi = l+lpp
     243             :                if ((lmini.le.lp).and.(lp.le.lmaxi).and.
     244           0 :      &            (mod(l+lp+lpp,2).eq.0).and.(abs(mpp).LE.lpp))then  
     245             :                   
     246             :                   factor=conjg(ylmpp(lmpp+1))*(ic**(l+lpp-lp))*
     247           0 :      *                     gaunt1(lp,lpp,l,mp,mpp,m,lmaxd)
     248             : 
     249             : c                  if(factor.ne.0 .and. uju(l,lp,lpp).ne.0)
     250             : c     >               write(*,*)lpp,lp,l,mp,m
     251             : 
     252           0 :                   if(l_q) then
     253           0 :                     factor=(sign_q**lpp)*factor            ! additional sign for q points
     254             :                   endif
     255             :      
     256             :                 ujug(lpmp,lm,n,ikpt_b)=ujug(lpmp,lm,n,ikpt_b)+
     257           0 :      +               factor*uju(l,lp,lpp)
     258             :                 ujdg(lpmp,lm,n,ikpt_b)=ujdg(lpmp,lm,n,ikpt_b)+
     259           0 :      +               factor*ujd(l,lp,lpp)
     260             :                 djug(lpmp,lm,n,ikpt_b)=djug(lpmp,lm,n,ikpt_b)+
     261           0 :      +               factor*dju(l,lp,lpp)
     262             :                 djdg(lpmp,lm,n,ikpt_b)=djdg(lpmp,lm,n,ikpt_b)+
     263           0 :      +               factor*djd(l,lp,lpp)
     264             :               
     265             :                endif
     266             :               enddo  ! lpp
     267             :             enddo ! mp
     268             :            enddo  ! lp
     269             :           enddo  ! m
     270             :          enddo   ! l
     271             : c******************************************************************
     272             : c       multiply with the gaunt-coefficient (apw-lo)
     273             : c******************************************************************
     274           0 :          if (nlo(n).ge.1) then 
     275           0 :          do lo = 1,nlo(n) 
     276           0 :           l = llo(lo,n)
     277           0 :           do m = -l,l
     278           0 :            lm=l*(l+1)+m  
     279           0 :            do lp = 0,lwn
     280           0 :             do mp = -lp,lp
     281           0 :               lpmp=lp*(lp+1)+mp 
     282           0 :               do lpp = 0,lwn
     283           0 :                llpp = lpp*(lpp+1)
     284           0 :                lmini = abs(l-lpp)
     285           0 :                lmaxi = l+lpp
     286           0 :                lminp = abs(lp-lpp)
     287           0 :                lmaxp = lp+lpp
     288             :                if ((lmini.le.lp).and.(lp.le.lmaxi).and.
     289           0 :      &            (mod(l+lp+lpp,2).eq.0).and.(abs(mp-m).LE.lpp)) then
     290           0 :                 mpp = mp - m
     291           0 :                 lmpp = llpp + mpp
     292             :                 factor=conjg(ylmpp(lmpp+1))*(ic**(l+lpp-lp))*
     293           0 :      *                   gaunt1(lp,lpp,l,mp,mpp,m,lmaxd)
     294           0 :                 if(l_q) then
     295           0 :                   factor=(sign_q**lpp)*factor            ! additional sign for q points
     296             :                 endif
     297             : 
     298             :                 ujulog(lpmp,lo,m,n,ikpt_b)=ujulog(lpmp,lo,m,n,ikpt_b)+
     299           0 :      +               factor*ujulo(lo,lp,lpp)
     300             :                 djulog(lpmp,lo,m,n,ikpt_b)=djulog(lpmp,lo,m,n,ikpt_b)+
     301           0 :      +               factor*djulo(lo,lp,lpp)
     302             :                endif
     303             : 
     304             :                if ((lminp.le.l).and.(l.le.lmaxp).and.
     305           0 :      &            (mod(l+lp+lpp,2).eq.0).and.(abs(m-mp).LE.lpp)) then
     306           0 :                 mpp = m - mp
     307           0 :                 lmpp = llpp + mpp
     308             :                 factor=conjg(ylmpp(lmpp+1))*(ic**(lp+lpp-l))*
     309           0 :      *                   gaunt1(l,lpp,lp,m,mpp,mp,lmaxd)
     310           0 :                 if(l_q) then
     311           0 :                    factor=(sign_q**lpp)*factor
     312             :                 endif
     313             :              
     314             :                 ulojug(lpmp,lo,m,n,ikpt_b)=ulojug(lpmp,lo,m,n,ikpt_b)+
     315           0 :      +               factor*uloju(lo,lp,lpp)         
     316             :                 ulojdg(lpmp,lo,m,n,ikpt_b)=ulojdg(lpmp,lo,m,n,ikpt_b)+
     317           0 :      +               factor*ulojd(lo,lp,lpp)        
     318             :                endif
     319             :               enddo  ! lpp
     320             :             enddo ! mp
     321             :            enddo  ! lp
     322             :           enddo ! m lo
     323             :          enddo  ! lo
     324             : c*************************************************************
     325             : c         multiply with the gaunt-coefficient (lo-lo)
     326             : c*************************************************************
     327           0 :          do lo = 1,nlo(n)
     328           0 :           l = llo(lo,n)
     329           0 :           do m = -l,l
     330           0 :            lm=l*(l+1)+m  
     331           0 :            do lop = 1,nlo(n)
     332           0 :             lp = llo(lop,n)
     333           0 :             do mp = -lp,lp
     334           0 :               lpmp=lp*(lp+1)+mp 
     335           0 :               do lpp = 0,lwn
     336           0 :                llpp = lpp*(lpp+1)
     337           0 :                mpp = mp - m 
     338           0 :                lmpp = llpp + mpp
     339           0 :                lmini = abs(l-lpp)
     340           0 :                lmaxi = l+lpp
     341             :                if ((lmini.le.lp).and.(lp.le.lmaxi).and.
     342           0 :      &            (mod(l+lp+lpp,2).eq.0).and.(abs(mpp).LE.lpp))then  
     343             :                   
     344             :                   factor= conjg(ylmpp(lmpp+1))*(ic**(l+lpp-lp))*
     345           0 :      *                 gaunt1(lp,lpp,l,mp,mpp,m,lmaxd)
     346           0 :                   if(l_q) then
     347           0 :                     factor=(sign_q**lpp)*factor            ! additional sign for q points
     348             :                   endif
     349             : 
     350             :                   ulojulog(lop,mp,lo,m,n,ikpt_b)=
     351             :      =                 ulojulog(lop,mp,lo,m,n,ikpt_b)+
     352           0 :      +                    ulojulo(lo,lop,lpp)*factor
     353             : !     +                 conjg(ylmpp(lmpp+1))*(ic**(l+lpp-lp))*
     354             : !     *                 gaunt1(lp,lpp,l,mp,mpp,m,lmaxd)*
     355             : !     *                      ulojulo(lo,lop,lpp)
     356             :                endif
     357             :               enddo  ! lpp
     358             :             enddo ! mp lop
     359             :            enddo  ! lop
     360             :           enddo ! m lo
     361             :          enddo  ! lo           
     362             :          endif ! local orbitals on this atom
     363             :         enddo !ntype
     364             :       enddo !ikpt_b   
     365           0 :       deallocate(djd,ujd,uju,ujulo,djulo,ulojulo)
     366           0 :       deallocate(dju,uloju,ulojd)
     367             : 
     368           0 :       end subroutine wann_ujugaunt
     369             :       end module m_wann_ujugaunt

Generated by: LCOV version 1.13