wann_mmk0_od_vac.F 6.88 KB
Newer Older
1 2 3 4 5 6
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------

7 8 9 10 11 12 13 14 15 16 17 18 19
      MODULE m_wann_mmk0_od_vac
      use m_juDFT
c**************************************************************
c      Determines the overlap matrix Mmn(k) in the vacuum
c      for the wannier functions.
c      For more details see routine wannier.F 
c
c      Y. Mokrousov, F. Freimuth 
c*************************************************************** 

      CONTAINS

      SUBROUTINE wann_mmk0_od_vac(
20
     >     DIMENSION, oneD, vacuum, stars, cell,
21 22 23 24 25
     >     l_noco,nlotot,
     >     z1,nmzxyd,nmzd,nv2d,k1d,k2d,k3d,n2d,n3d,
     >     ig,nmzxy,nmz,delz,ig2,n2d_1,
     >     bbmat,evac,bkpt,MM,vM,vz,odi,
     >     nslibd,jspin,k1,k2,k3,jspd,nvd,area,
26
     >     nbasfcn,neigd,zMat,nv,sk2,phi2,omtil,qss,
27 28 29
     <     mmn)

      use m_constants, only : pimach
30
      use m_types
31 32 33
      use m_od_abvac
      use m_cylbes
      use m_dcylbs
34
      USE m_types
35 36

      implicit none
37

38
      TYPE(t_mat), INTENT(IN) :: zMat
39

40 41 42 43 44 45
      TYPE(t_dimension),INTENT(IN)   :: DIMENSION
      TYPE(t_oneD),INTENT(IN)        :: oneD
      TYPE(t_vacuum),INTENT(IN)      :: vacuum
      TYPE(t_stars),INTENT(IN)       :: stars
      TYPE(t_cell),INTENT(IN)        :: cell

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
c     .. scalar Arguments..
      logical, intent (in) :: l_noco
      integer, intent (in) :: nlotot
      integer, intent (in) :: nmzxyd,nmzd,nv2d,k1d,k2d,k3d,n3d
      integer, intent (in) :: nmzxy,nmz,MM,n2d,vM,nslibd
      integer, intent (in) :: n2d_1,jspin,jspd,nvd
      integer, intent (in) :: nbasfcn,neigd
      real,    intent (in) :: delz,z1,evac,area,omtil
      type (od_inp), intent (in) :: odi

c     ..array arguments..
      real,    intent (in) :: bkpt(3),qss(3)
      real,    intent (in) :: sk2(:),phi2(:) !sk2(n2d),phi2(n2d)
      integer, intent (in) :: ig(-k1d:,-k2d:,-k3d:) !ig(-k1d:k1d,-k2d:k2d,-k3d:k3d)
      integer, intent (in) :: ig2(:),nv(:) !ig2(n3d),nv(jspd)
      real,    intent (in) :: vz(:),bbmat(3,3)!vz(nmzd),bbmat(3,3)
      integer, intent (in) :: k1(:,:),k2(:,:),k3(:,:)!k1(nvd,jspd),k2(nvd,jspd),k3(nvd,jspd)
      complex, intent (inout) :: mmn(:,:)!mmn(nslibd,nslibd)

c     ..basis wavefunctions in the vacuum
      real,    allocatable :: udz(:,:)
      real,    allocatable :: uz(:,:)
      real,    allocatable :: dudz(:,:)
      real,    allocatable :: duz(:,:)
      real,    allocatable :: u(:,:,:)
      real,    allocatable :: ud(:,:,:)
      real,    allocatable :: ddnv(:,:)

c     ..local scalars..
      real :: wronk,wronk1,arg,zks,tpi
      integer :: i,m,l,j,k,irec3,irec2,n,nv2,ispin,addnoco
      complex :: avac,bvac,ic
      complex, allocatable :: acof(:,:,:),bcof(:,:,:)
      integer, allocatable :: kvac3(:),map1(:)
      real, allocatable :: bess(:),dbss(:)
      real :: qssbti(3,2)

c     ..intrinsic functions..
      intrinsic aimag,cmplx,conjg,real,sqrt

      allocate ( udz(nv2d,-vM:vM),uz(nv2d,-vM:vM),
     +           dudz(nv2d,-vM:vM),
     +           duz(nv2d,-vM:vM),u(nmzd,nv2d,-vM:vM),
     +           ud(nmzd,nv2d,-vM:vM),
     +           ddnv(nv2d,-vM:vM),
     +           bess(-odi%mb:odi%mb),dbss(-odi%mb:odi%mb),
     +           acof(nv2d,-odi%mb:odi%mb,nslibd),
     +           bcof(nv2d,-odi%mb:odi%mb,nslibd),
     +           kvac3(nv2d),map1(nvd) )

      acof=cmplx(0.0,0.0)
      bcof=cmplx(0.0,0.0)

      tpi = 2 * pimach() ; ic = cmplx(0.,1.)

      nv2 = 0

      do 20 k = 1,nv(jspin)
         do 10 j = 1,nv2
            if (k3(k,jspin).eq.kvac3(j)) then
               map1(k) = j
               goto 20
            endif
 10      continuE
         nv2 = nv2 + 1
         IF (nv2>nv2d)  CALL juDFT_error("nv2d",calledby
     +        ="wann_mmk0_od_vac")
         kvac3(nv2) = k3(k,jspin)
         map1(k) = nv2
 20   continue

      wronk = 2.0

      qssbti(1,1) = - qss(1)/2.     ! noco (ss) case not
      qssbti(2,1) = - qss(2)/2.     ! implemented, just
      qssbti(1,2) = + qss(1)/2.     ! for compatibility
      qssbti(2,2) = + qss(2)/2.
      qssbti(3,1) = - qss(3)/2.
      qssbti(3,2) = + qss(3)/2.
      DO ispin = 1,1 ! jspins
      CALL od_abvac(
127 128 129
     >      cell,vacuum,DIMENSION,stars,oneD,
     >      qssbti(3,jspin),odi%n2d,
     >      wronk,evac,bkpt,odi%M,odi%mb,
130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158
     >      vz,kvac3(1),nv2,
     <      uz(1,-vM),duz(1,-vM),u(1,1,-vM),udz(1,-vM),
     <      dudz(1,-vM),ddnv(1,-vM),ud(1,1,-vM))
      ENDDO

      addnoco=0
      if(l_noco.and.(jspin.eq.2))then
         addnoco= nv(1)+nlotot
      endif
      do k = 1,nv(jspin)
         l = map1(k)
         irec3 = ig(k1(k,jspin),k2(k,jspin),k3(k,jspin))
         if (irec3.ne.0) then
            irec2 = ig2(irec3)
            zks = sk2(irec2)*z1
            arg = phi2(irec2)
            call cylbes(odi%mb,zks,bess)
            call dcylbs(odi%mb,zks,bess,dbss)
            do m = -odi%mb,odi%mb
               wronk1 = uz(l,m)*dudz(l,m) -
     -              udz(l,m)*duz(l,m)
               avac = exp(-cmplx(0.0,m*arg))*(ic**m)*
     *              cmplx(dudz(l,m)*bess(m) -
     +              udz(l,m)*sk2(irec2)*dbss(m),0.0)/
     /              ((wronk1)*sqrt(omtil))
               bvac = exp(-cmplx(0.0,m*arg))*(ic**m)*
     *              cmplx(-duz(l,m)*bess(m) +
     -              uz(l,m)*sk2(irec2)*dbss(m),0.0)/
     /              ((wronk1)*sqrt(omtil))
159 160 161
               IF(zMat%l_real) THEN
                  do n = 1,nslibd
                      acof(l,m,n) = acof(l,m,n) +
162 163
     +                   zMat%data_r(k+addnoco,n)*avac
c     +                    conjg(zMat%data_r(k,n))*avac
164
                      bcof(l,m,n) = bcof(l,m,n) +
165 166
     +                   zMat%data_r(k+addnoco,n)*bvac
c     +                    conjg(zMat%data_r(k,n))*bvac
167 168 169 170
                  enddo
               ELSE
                  do n = 1,nslibd
                      acof(l,m,n) = acof(l,m,n) +
171 172
     +                   zMat%data_c(k+addnoco,n)*avac
c     +                    conjg(zMat%data_c(k,n))*avac
173
                      bcof(l,m,n) = bcof(l,m,n) +
174 175
     +                   zMat%data_c(k+addnoco,n)*bvac
c     +                    conjg(zMat%data_c(k,n))*bvac
176 177
                  enddo
               END IF
178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199
            enddo      ! -mb:mb
         endif 
      enddo          ! k = 1,nv

c  now actually computing the Mmn matrix

      do l = 1,nv2
       do m = -odi%mb,odi%mb
         do i = 1,nslibd
           do j = 1,nslibd 
            mmn(i,j) = mmn(i,j) + 
     +                area*(acof(l,m,i)*conjg(acof(l,m,j))
     +              + ddnv(l,m)*bcof(l,m,i)*conjg(bcof(l,m,j)))  
           enddo 
         enddo
       enddo
      enddo

      deallocate ( udz,uz,dudz,duz,u,ud,ddnv,bess,dbss,acof,bcof )

      END SUBROUTINE wann_mmk0_od_vac
      END MODULE m_wann_mmk0_od_vac