LCOV - code coverage report
Current view: top level - wannier - wann_mmkb_sph.F (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 64 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_mmkb_sph
       8             :       contains
       9           0 :       subroutine wann_mmkb_sph(
      10             : c***********************************************************************
      11             : c   computes the Mmn(k,b) matrix elements which are the overlaps
      12             : c   between the Bloch wavefunctions at the state m and n
      13             : c   at the k-point and its nearest neighbor b (bpt), in the spheres
      14             : c                                Y.Mokrousov 15.6.06
      15             : c***********************************************************************
      16             : c   the gaunt coefficients impose certain conditions (which can be
      17             : c   seen in the gaunt.F) on the l,lp,lpp angular momenta, which can be
      18             : c   used in order to speedup the procedure. 
      19             : c   The conditions are the following:
      20             : c   1. |l - l''| <= l' <= l + l''
      21             : c   2. m' = m + m''
      22             : c   3. l + l' + l''  - is even
      23             : c*********************************************************************** 
      24             : c   modified for use in conjunction with wann_ujugaunt, which
      25             : c   calculates certain matrix elements that have to be treated only
      26             : c   once per calculation
      27             : c   Frank Freimuth, October 2006
      28             : c***********************************************************************
      29             : c   use BLAS matrix-matrix multiplication in apw-apw part to be faster
      30             : c   J.-P. Hanke, Dec. 2015
      31             : c***********************************************************************
      32             :      >        nbnd,llod,nslibd,nslibd_b,nlod,natd,ntypd,lmd,jmtd,
      33           0 :      >        taual,nop,lmax,
      34           0 :      >        ntype,neq,nlo,llo,acof,bcof,ccof,bbpt,
      35           0 :      >        acof_b,bcof_b,ccof_b,gb,bkpt,ujug,ujdg,djug,djdg,ujulog,
      36           0 :      >        djulog,ulojug,ulojdg,ulojulog,kdiff,nntot,chi,
      37           0 :      =        mmn)
      38             : #include "cpp_double.h"
      39             :       use m_juDFT
      40             :       use m_constants, only : pimach
      41             :       use m_matmul   , only : matmul3,matmul3r
      42             :       use m_sphbes
      43             :       use m_ylm
      44             :       use m_intgr, only : intgr3
      45             :       use m_gaunt, only: gaunt1
      46             :     
      47             :       implicit none
      48             : 
      49             : c     .. scalar arguments ..
      50             :       integer, intent (in) :: llod,nlod,natd,ntypd,lmd,nbnd
      51             :       integer, intent (in) :: nntot
      52             :       integer, intent (in) :: ntype,nslibd,nslibd_b,nop,jmtd
      53             : 
      54             : c     .. array arguments ..
      55             :       integer, intent (in)  :: neq(:) !(ntypd)
      56             :       integer, intent (in)  :: lmax(:) !(ntypd)
      57             :       integer, intent (in)  :: nlo(:) !(ntypd)
      58             :       integer, intent (in)  :: llo(:,:) !(nlod,ntypd)
      59             :       complex, intent (in)  :: chi(ntypd)
      60             :       real,    intent (in)  :: bbpt(:) !(3)
      61             :       real,    intent (in)  :: taual(:,:) !(3,natd)
      62             :       real,    intent (in)  :: bkpt(:) !(3)
      63             :       integer, intent (in)  :: gb(:) !(3)
      64             :       complex, intent (in)  :: ccof(-llod:,:,:,:) !(-llod:llod,nslibd,nlod,natd)
      65             :       complex, intent (in)  :: acof(:,0:,:) !(nslibd,0:lmd,natd)
      66             :       complex, intent (in)  :: bcof(:,0:,:) !(nslibd,0:lmd,natd)
      67             :       complex, intent (in)  :: ccof_b(-llod:,:,:,:) !(-llod:llod,nslibd_b,nlod,natd)
      68             :       complex, intent (in)  :: acof_b(:,0:,:) !(nslibd_b,0:lmd,natd)
      69             :       complex, intent (in)  :: bcof_b(:,0:,:) !(nslibd_b,0:lmd,natd)
      70             : 
      71             :       complex, intent (in)  :: ujug(0:,0:,:,:) !(0:lmd,0:lmd,1:ntype,1:nntot)
      72             :       complex, intent (in)  :: ujdg(0:,0:,:,:) !(0:lmd,0:lmd,1:ntype,1:nntot)
      73             :       complex, intent (in)  :: djug(0:,0:,:,:) !(0:lmd,0:lmd,1:ntype,1:nntot)
      74             :       complex, intent (in)  :: djdg(0:,0:,:,:) !(0:lmd,0:lmd,1:ntype,1:nntot)
      75             :       complex, intent (in)  :: ujulog(0:,:,-llod:,:,:) !(0:lmd,nlod,-llod:llod,1:ntype,1:nntot)
      76             :       complex, intent (in)  :: djulog(0:,:,-llod:,:,:) !(0:lmd,nlod,-llod:llod,1:ntype,1:nntot)
      77             :       complex, intent (in)  :: ulojug(0:,:,-llod:,:,:) !(0:lmd,nlod,-llod:llod,1:ntype,1:nntot)
      78             :       complex, intent (in)  :: ulojdg(0:,:,-llod:,:,:) !(0:lmd,nlod,-llod:llod,1:ntype,1:nntot)
      79             :       complex, intent (in)  :: ulojulog(:,-llod:,:,-llod:,:,:) !(1:nlod,-llod:llod,1:nlod,-llod:llod,1:ntype,1:nntot)
      80             : 
      81             :       real, intent (in)     :: kdiff(:,:) !(3,nntot)
      82             :       complex,intent(inout) :: mmn(:,:)!nbnd,nbnd) !(nbnd,nbnd)
      83             : 
      84             : 
      85             : c     .. local scalars ..
      86             :       integer i,lm,nn,n,na,j,lmp,l,lp,m,mp,lwn,lo,lop
      87             :       integer ll,llp
      88             :       real rph,cph,tpi,th,t1nn,t2nn,t3nn,dummy
      89             :       complex ic
      90             :       integer nene,lmd2
      91             :       complex :: fac1,fac2,fac3,fac4
      92             :       complex :: fac5,fac6,fac7,fac8
      93           0 :       complex :: mat(0:lmd,nslibd_b)
      94             : C     ..
      95             : C     .. local arrays ..
      96             :       real bpt(3)
      97             : 
      98             : C     ..
      99             : C     .. intrinsic functions ..
     100             :       intrinsic conjg,cmplx,sqrt,cos,sin
     101             : 
     102           0 :       ic = cmplx(0.,1.)
     103           0 :       tpi = 2* pimach()
     104           0 :       lmd2 = lmd+1
     105             :       
     106           0 :       na = 0
     107             : 
     108           0 :       bpt(:) = bbpt(:) + gb(:) - bkpt(:)
     109           0 :       do nene=1,nntot
     110           0 :          if(all(abs(bpt(:)-kdiff(:,nene)).lt.1e-4)) exit
     111             :       enddo  
     112           0 :       IF(nene==nntot+1) CALL juDFT_error
     113             :      +     ("cannot find matching nearest neighbor k",calledby
     114           0 :      +     ="wann_mmkb_sph")
     115             : 
     116           0 :       do n=1,ntype
     117           0 :        lwn = lmax(n)
     118           0 :          do nn = 1,neq(n) ! cycle by the atoms within the atom type
     119           0 :          na = na + 1
     120             : c...set up phase factors ( e^{ib\tau} )
     121             : 
     122           0 :          t1nn =  tpi*taual(1,na)
     123           0 :          t2nn =  tpi*taual(2,na)
     124           0 :          t3nn =  tpi*taual(3,na)
     125             : 
     126           0 :          th = bpt(1)*t1nn + bpt(2)*t2nn + bpt(3)*t3nn
     127           0 :          rph = 2*tpi*cos(th)
     128           0 :          cph = 2*tpi*sin(th)
     129             : 
     130             : c...contributions from the apws
     131             :          call CPP_BLAS_cgemm(
     132             :      >          'T','C',lmd2,nslibd_b,lmd2,cmplx(rph,cph),
     133             :      >          ujug(0,0,n,nene),lmd2,acof_b(1,0,na),nslibd_b,
     134           0 :      >          cmplx(0.0),mat(0,1),lmd2)
     135             :          call CPP_BLAS_cgemm(
     136             :      >          'N','N',nslibd,nslibd_b,lmd2,chi(n),
     137             :      >          acof(1,0,na),nslibd,mat(0,1),lmd2,
     138           0 :      >          cmplx(1.0),mmn,nbnd)
     139             : 
     140             :          call CPP_BLAS_cgemm(
     141             :      >          'T','C',lmd2,nslibd_b,lmd2,cmplx(rph,cph),
     142             :      >          ujdg(0,0,n,nene),lmd2,bcof_b(1,0,na),nslibd_b,
     143           0 :      >          cmplx(0.0),mat(0,1),lmd2)
     144             :          call CPP_BLAS_cgemm(
     145             :      >          'N','N',nslibd,nslibd_b,lmd2,chi(n),
     146             :      >          acof(1,0,na),nslibd,mat(0,1),lmd2,
     147           0 :      >          cmplx(1.0),mmn,nbnd)
     148             : 
     149             :          call CPP_BLAS_cgemm(
     150             :      >          'T','C',lmd2,nslibd_b,lmd2,cmplx(rph,cph),
     151             :      >          djug(0,0,n,nene),lmd2,acof_b(1,0,na),nslibd_b,
     152           0 :      >          cmplx(0.0),mat(0,1),lmd2)
     153             :          call CPP_BLAS_cgemm(
     154             :      >          'N','N',nslibd,nslibd_b,lmd2,chi(n),
     155             :      >          bcof(1,0,na),nslibd,mat(0,1),lmd2,
     156           0 :      >          cmplx(1.0),mmn,nbnd)
     157             : 
     158             :          call CPP_BLAS_cgemm(
     159             :      >          'T','C',lmd2,nslibd_b,lmd2,cmplx(rph,cph),
     160             :      >          djdg(0,0,n,nene),lmd2,bcof_b(1,0,na),nslibd_b,
     161           0 :      >          cmplx(0.0),mat(0,1),lmd2)
     162             :          call CPP_BLAS_cgemm(
     163             :      >          'N','N',nslibd,nslibd_b,lmd2,chi(n),
     164             :      >          bcof(1,0,na),nslibd,mat(0,1),lmd2,
     165           0 :      >          cmplx(1.0),mmn,nbnd)
     166             : 
     167             : 
     168             : 
     169           0 :          if (nlo(n).ge.1) then 
     170             :          
     171           0 :          do lo = 1,nlo(n) 
     172           0 :           l = llo(lo,n)
     173           0 :           ll = l*(l+1)
     174           0 :           do m = -l,l
     175             :            lm = ll + m
     176             : 
     177           0 :            do lp = 0,lwn
     178           0 :             llp = lp*(lp+1)
     179           0 :             do mp = -lp,lp
     180           0 :              lmp = llp + mp
     181             : 
     182           0 :              fac1=cmplx(rph,cph)*ujulog(lmp,lo,m,n,nene)*chi(n)
     183           0 :              fac2=cmplx(rph,cph)*djulog(lmp,lo,m,n,nene)*chi(n)
     184           0 :              fac3=cmplx(rph,cph)*ulojdg(lmp,lo,m,n,nene)*chi(n)
     185           0 :              fac4=cmplx(rph,cph)*ulojug(lmp,lo,m,n,nene)*chi(n)
     186             :              
     187             : c             if(fac1.ne.0)write(*,*)'fac1',lmp,lm
     188             : c             if(fac2.ne.0)write(*,*)'fac2',lmp,lm
     189             : c             if(fac3.ne.0)write(*,*)'fac3',lmp,lm
     190             : c             if(fac4.ne.0)write(*,*)'fac4',lmp,lm
     191             : 
     192           0 :                 do j = 1,nslibd_b
     193           0 :                  do i = 1,nslibd
     194             :                   mmn(i,j) = mmn(i,j) +
     195             :      +                 ccof(m,i,lo,na)*
     196             :      *              (conjg(acof_b(j,lmp,na))*fac1+
     197             :      +               conjg(bcof_b(j,lmp,na))*fac2 )+
     198             : 
     199             :      +                 conjg(ccof_b(m,j,lo,na))*
     200             :      *              ( bcof(i,lmp,na)*fac3 +
     201           0 :      +                acof(i,lmp,na)*fac4 )
     202             :                  enddo 
     203             :                 enddo  
     204             : 
     205             : 
     206             :             enddo ! mp
     207             :            enddo  ! lp
     208             : 
     209             :           enddo ! lo
     210             :          enddo  ! m lo
     211             : 
     212             : c...contributions from lo*lo
     213             :           
     214           0 :          do lo = 1,nlo(n)
     215           0 :           l = llo(lo,n)
     216           0 :           ll = l*(l+1)
     217           0 :           do m = -l,l
     218             :            lm = ll + m
     219             : 
     220           0 :            do lop = 1,nlo(n)
     221           0 :             lp = llo(lop,n)
     222           0 :             llp = lp*(lp+1)
     223           0 :             do mp = -lp,lp
     224           0 :              lmp = llp + mp
     225             : 
     226           0 :              fac1=cmplx(rph,cph)*ulojulog(lop,mp,lo,m,n,nene)*chi(n)
     227             : 
     228             : c             if(fac1.ne.0)write(*,*)'fac1',lop,lo,mp,m
     229             : 
     230           0 :                 do j = 1,nslibd_b
     231           0 :                  do i = 1,nslibd
     232             : 
     233             :                   mmn(i,j) = mmn(i,j) +
     234             :      +               ccof(m,i,lo,na)*conjg(ccof_b(mp,j,lop,na))*
     235           0 :      *                        fac1
     236             : 
     237             :                  enddo 
     238             :                 enddo  
     239             :                 
     240             :             enddo ! mp lop
     241             :            enddo  ! lop
     242             : 
     243             :           enddo ! m lo
     244             :          enddo  ! lo           
     245             :          
     246             :          endif ! local orbitals on this atom
     247             : 
     248             :       enddo  ! atoms in the type
     249             : 
     250             :       enddo  ! atom type   
     251             : 
     252             : 
     253           0 :       end subroutine wann_mmkb_sph
     254             :       end module m_wann_mmkb_sph

Generated by: LCOV version 1.13