wann_mmkb_vac.F 11.2 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
      MODULE m_wann_mmkb_vac
      use m_juDFT
c**************************************************************
c      Determines the overlap matrix Mmn(k) in the vacuum
c      in the film case for the wannier functions.
c      For more details see routine wannier.F 
c
c      Y. Mokrousov, F. Freimuth
c*************************************************************** 
      CONTAINS
      SUBROUTINE wann_mmkb_vac(
18
     >     vacchi,l_noco,nlotot,qss,
19 20
     >     nbnd,z1,nmzd,nv2d,k1d,k2d,k3d,n3d,nvac,
     >     ig,nmz,delz,ig2,area,bmat,
21 22
     >     bbmat,evac,evac_b,bkpt,bkpt_b,vz,vz_b,
     >     nslibd,nslibd_b,jspin,jspin_b,
23
     >     k1,k2,k3,k1_b,k2_b,k3_b,jspd,nvd,
24
     >     nbasfcn,neigd,zMat,zMat_b,nv,nv_b,omtil,gb,
25 26
     <     mmn)

27
      USE m_types
28
      use m_constants
29
      use m_intgr, only : intgz0
30 31 32
      USE m_vacuz
      USE m_vacudz

33
      implicit none
34 35 36

      TYPE(t_zmat), INTENT(IN) :: zMat, zMat_b

37 38
c     .. scalar Arguments..
      logical, intent (in) :: l_noco
39
      integer, intent (in) :: nlotot,jspin_b
40 41 42 43 44 45 46
      real,    intent (in) :: qss(3)
      integer, intent (in) :: nmzd,nv2d,k1d,k2d,k3d,n3d,nbnd
      integer, intent (in) :: nmz,nslibd,nslibd_b,nvac
      integer, intent (in) :: jspin,jspd,nvd
      integer, intent (in) :: nbasfcn,neigd
      integer, intent (in) :: gb(3)
      real,    intent (in) :: delz,z1,omtil,area
47
      complex, intent (in) :: vacchi
48 49

c     ..array arguments..
50
      real,    intent (in) :: bkpt(3),bkpt_b(3),evac(2),evac_b(2)
51 52
      integer, intent (in) :: ig(-k1d:k1d,-k2d:k2d,-k3d:k3d)
      integer, intent (in) :: ig2(n3d),nv(jspd),nv_b(jspd)
53 54
      real,    intent (in) :: vz(nmzd,2),vz_b(nmzd,2)
      real,    intent (in) :: bbmat(3,3),bmat(3,3)
55 56 57 58 59 60 61 62 63 64 65 66 67 68 69
      integer, intent (in) :: k1(nvd,jspd),k2(nvd,jspd),k3(nvd,jspd)
      integer,intent(in)::k1_b(nvd,jspd),k2_b(nvd,jspd),k3_b(nvd,jspd)
      complex, intent (inout) :: mmn(nbnd,nbnd)

c     ..basis wavefunctions in the vacuum
      complex, allocatable :: ac(:,:),bc(:,:),ac_b(:,:),bc_b(:,:)
      real,    allocatable :: dt(:),dte(:)
      real,    allocatable :: t(:),te(:),tei(:)
      real,    allocatable :: u(:,:),ue(:,:),v(:)
      real,    allocatable :: dt_b(:),dte_b(:)
      real,    allocatable :: t_b(:),te_b(:),tei_b(:)
      real,    allocatable :: u_b(:,:),ue_b(:,:)

c     ..local scalars..
      logical tail
70
      real wronk,arg,zks,tpi,vz0(2),vz0_b(2),scale,evacp,ev,const
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
      real uu,ud,du,dd,xx(nmz),xximag(nmz)
      real :: uuimag,udimag,duimag,ddimag,qss1,qss2
      integer i,m,l,j,k,n,nv2,nv2_b,ivac,n2,n2_b,sign,ik
      integer :: lprime,np1,addnoco,addnoco2
      complex :: av,bv,ic,c_1
      integer, allocatable :: kvac1(:),kvac2(:),map2(:)
      integer, allocatable :: kvac1_b(:),kvac2_b(:),map2_b(:)
      integer symvac,symvacvac

      allocate ( ac(nv2d,nslibd),bc(nv2d,nslibd),
     +           ac_b(nv2d,nslibd_b),bc_b(nv2d,nslibd_b),
     +           dt(nv2d),dte(nv2d),t(nv2d),te(nv2d),
     +           tei(nv2d),u(nmzd,nv2d),ue(nmzd,nv2d),
     +           dt_b(nv2d),dte_b(nv2d),t_b(nv2d),te_b(nv2d),
     +           tei_b(nv2d),u_b(nmzd,nv2d),ue_b(nmzd,nv2d),
     +           v(3),kvac1(nv2d),kvac2(nv2d),map2(nvd),
     +           kvac1_b(nv2d),kvac2_b(nv2d),map2_b(nvd) )

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

      tail = .true.
      np1 = nmz + 1

c.. determining the indexing array (in-plane stars)
c.. for the k-point

      wronk = 2.0
      const = 1.0 / ( sqrt(omtil)*wronk )

      do ivac = 1,2
         vz0(ivac) = vz(nmz,ivac)
102
         vz0_b(ivac) = vz_b(nmz,ivac)
103 104 105 106 107 108 109 110
      enddo

      n2 = 0 ; n2_b = 0

      addnoco=0
      addnoco2=0
      if(l_noco.and.jspin.eq.2)then
        addnoco=nv(1)+nlotot
111 112
      endif
      if(l_noco.and.jspin_b.eq.2)then
113 114 115 116 117 118 119 120 121 122 123 124
        addnoco2=nv_b(1)+nlotot
      endif

      do 40 k = 1,nv(jspin)
         do 30 j = 1,n2
            if ( k1(k,jspin).eq.kvac1(j) .and.
     +          k2(k,jspin).eq.kvac2(j) ) then
                map2(k) = j
                goto 40
             endif 
 30      continue
         n2 = n2 + 1
125 126 127 128 129
         
         IF(n2>nv2d) then
            write(*,*)n2,nv2d,'jspin',jspin
         endif
 
130 131 132 133 134 135 136
         IF (n2>nv2d)  CALL juDFT_error("wannier Mmn vac",calledby
     +        ="wann_mmkb_vac")

         kvac1(n2) = k1(k,jspin)
         kvac2(n2) = k2(k,jspin)
         map2(k) = n2
 40   continue
137
         !write(*,*)'ok',n2,nv2d,'jspin',jspin
138 139 140

c.. and for the b-point
 
141
      do 41 k = 1,nv_b(jspin_b)
142
         do 31 j = 1,n2_b
143 144
            if ( k1_b(k,jspin_b).eq.kvac1_b(j) .and.
     +          k2_b(k,jspin_b).eq.kvac2_b(j) ) then
145 146 147 148 149 150
                map2_b(k) = j
                goto 41
             endif
 31      continue
         n2_b = n2_b + 1
     
151 152 153 154
         IF(n2_b>nv2d) then
            write(*,*)n2_b,nv2d,'jspin_b',jspin_b
         endif

155 156 157
         IF (n2_b>nv2d)  CALL juDFT_error("wannier Mmn vac",calledby
     +        ="wann_mmkb_vac")

158 159
         kvac1_b(n2_b) = k1_b(k,jspin_b)
         kvac2_b(n2_b) = k2_b(k,jspin_b)
160 161
         map2_b(k) = n2_b
 41   continue
162
         !write(*,*)'ok',n2_b,nv2d,'jspin_b',jspin_b
163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189

c...cycle by the vacua
      do 140 ivac = 1,nvac



      sign = 3. - 2.*ivac
      evacp = evac(ivac)

      nv2 = n2 ; nv2_b = n2_b

c.. the body of the routine

      qss1=0.0
      qss2=0.0
      if(l_noco.and.jspin.eq.1)then
        qss1=-qss(1)/2.0
        qss2=-qss(2)/2.0
      elseif(l_noco.and.jspin.eq.2)then
        qss1=qss(1)/2.0
        qss2=qss(2)/2.0
      endif

      do ik = 1,nv2
         v(1) = bkpt(1) + kvac1(ik) + qss1
         v(2) = bkpt(2) + kvac2(ik) + qss2
         v(3) = 0.
190
         ev = evacp - 0.5*dot_product(v,matmul(bbmat,v))
191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228
         call vacuz(ev,vz(1,ivac),vz0(ivac),nmz,delz,t(ik),dt(ik),
     +        u(1,ik))
         call vacudz(ev,vz(1,ivac),vz0(ivac),nmz,delz,te(ik),
     +        dte(ik),tei(ik),ue(1,ik),dt(ik),
     +        u(1,ik))
         scale = wronk/ (te(ik)*dt(ik)-dte(ik)*t(ik))
         te(ik) = scale*te(ik)
         dte(ik) = scale*dte(ik)
         tei(ik) = scale*tei(ik)
         do j = 1,nmz
            ue(j,ik) = scale*ue(j,ik)
         enddo
      enddo
c-----> construct a and b coefficients for the k-point
       symvacvac=1
       if (nvac==1) symvacvac=2
       do symvac=1,symvacvac
         do i = 1,nv2d
            do n = 1,nslibd
               ac(i,n) = cmplx(0.0,0.0)
               bc(i,n) = cmplx(0.0,0.0)
            enddo   
            do n = 1,nslibd_b
               ac_b(i,n) = cmplx(0.0,0.0)
               bc_b(i,n) = cmplx(0.0,0.0)
            enddo   
         enddo   

        if (symvac==2) sign=-1.0   

      do k = 1,nv(jspin)
         l = map2(k)
         zks = k3(k,jspin)*bmat(3,3)*sign
         arg = zks*z1
         c_1 = cmplx(cos(arg),sin(arg)) * const
         av = -c_1 * cmplx( dte(l),zks*te(l) )
         bv =  c_1 * cmplx(  dt(l),zks* t(l) )
c-----> loop over basis functions
229 230 231 232 233 234 235 236 237 238 239
         IF(zMat%l_real) THEN
            do n = 1,nslibd
               ac(l,n) = ac(l,n) + zMat%z_r(k+addnoco,n)*av
               bc(l,n) = bc(l,n) + zMat%z_r(k+addnoco,n)*bv
            enddo
         ELSE
            do n = 1,nslibd
               ac(l,n) = ac(l,n) + zMat%z_c(k+addnoco,n)*av
               bc(l,n) = bc(l,n) + zMat%z_c(k+addnoco,n)*bv
            enddo
         END IF
240 241 242 243
      enddo

c...now for the k+b point

244
      evacp = evac_b(ivac)
245 246 247 248
      do ik = 1,nv2_b
         v(1) = bkpt_b(1) + kvac1_b(ik) + qss1
         v(2) = bkpt_b(2) + kvac2_b(ik) + qss2
         v(3) = 0.
249
         ev = evacp - 0.5*dot_product(v,matmul(bbmat,v))
250 251 252
         call vacuz(ev,vz_b(1,ivac),vz0_b(ivac),nmz,delz,t_b(ik),
     +        dt_b(ik),u_b(1,ik))
         call vacudz(ev,vz_b(1,ivac),vz0_b(ivac),nmz,delz,te_b(ik),
253 254 255 256 257 258 259 260 261 262 263 264
     +        dte_b(ik),tei_b(ik),ue_b(1,ik),dt_b(ik),
     +        u_b(1,ik))
         scale = wronk/ (te_b(ik)*dt_b(ik)-dte_b(ik)*t_b(ik))
         te_b(ik) = scale*te_b(ik)
         dte_b(ik) = scale*dte_b(ik)
         tei_b(ik) = scale*tei_b(ik)
         do j = 1,nmz
            ue_b(j,ik) = scale*ue_b(j,ik)
         enddo
      enddo
c-----> construct a and b coefficients for the k+b point

265
      do k = 1,nv_b(jspin_b)
266
         l = map2_b(k)
267
         zks = k3_b(k,jspin_b)*bmat(3,3)*sign
268 269 270 271 272
         arg = zks*z1
         c_1 = cmplx(cos(arg),sin(arg)) * const
         av = -c_1 * cmplx( dte_b(l),zks*te_b(l) )
         bv =  c_1 * cmplx( dt_b(l),zks*t_b(l) )
c-----> loop over basis functions
273 274 275 276 277 278 279 280 281 282 283
         IF(zMat_b%l_real) THEN
            do n = 1,nslibd_b
               ac_b(l,n) = ac_b(l,n) + zMat_b%z_r(k+addnoco,n)*av
               bc_b(l,n) = bc_b(l,n) + zMat_b%z_r(k+addnoco,n)*bv
            enddo
         ELSE
            do n = 1,nslibd_b
               ac_b(l,n) = ac_b(l,n) + zMat_b%z_c(k+addnoco,n)*av
               bc_b(l,n) = bc_b(l,n) + zMat_b%z_c(k+addnoco,n)*bv
            enddo
         END IF
284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333
      enddo


      do l = 1,nv2
      do lprime = 1,nv2_b
      if (kvac1(l).eq.(kvac1_b(lprime)-gb(1))
     & .and. kvac2(l).eq.(kvac2_b(lprime)-gb(2)))then
         zks = gb(3)*bmat(3,3)*sign

         do i = 1,nmz
             xx(np1-i) = u(i,l)*u_b(i,lprime)*
     *          cos((z1+(i-1)*delz)*zks)
            xximag(np1-i) = u(i,l)*u_b(i,lprime)*
     *          sin((z1+(i-1)*delz)*zks)
         enddo   
         call intgz0(xx,delz,nmz,uu,tail)
         call intgz0(xximag,delz,nmz,uuimag,tail)

         do i = 1,nmz
            xx(np1-i) = u(i,l)*ue_b(i,lprime)*
     *   cos((z1+(i-1)*delz)*zks)
            xximag(np1-i) = u(i,l)*ue_b(i,lprime)*
     *          sin((z1+(i-1)*delz)*zks)
         enddo   
         call intgz0(xx,delz,nmz,ud,tail)
         call intgz0(xximag,delz,nmz,udimag,tail)

         do i = 1,nmz
            xx(np1-i) = ue(i,l)*u_b(i,lprime)*
     *   cos((z1+(i-1)*delz)*zks)
            xximag(np1-i) = ue(i,l)*u_b(i,lprime)*
     *          sin((z1+(i-1)*delz)*zks)
         enddo   
         call intgz0(xx,delz,nmz,du,tail)
         call intgz0(xximag,delz,nmz,duimag,tail)
         do i = 1,nmz
            xx(np1-i) = ue(i,l)*ue_b(i,lprime)*
     *   cos((z1+(i-1)*delz)*zks)
            xximag(np1-i) = ue(i,l)*ue_b(i,lprime)*
     *          sin((z1+(i-1)*delz)*zks)
         enddo   
         call intgz0(xx,delz,nmz,dd,tail)
         call intgz0(xximag,delz,nmz,ddimag,tail)

         do i = 1,nslibd
            do j = 1,nslibd_b
               mmn(i,j) = mmn(i,j) + area*(
     *  ac(l,i)*conjg(ac_b(lprime,j))*cmplx(uu,uuimag) +
     +  ac(l,i)*conjg(bc_b(lprime,j))*cmplx(ud,udimag) +
     *  bc(l,i)*conjg(ac_b(lprime,j))*cmplx(du,duimag) +
334
     +  bc(l,i)*conjg(bc_b(lprime,j))*cmplx(dd,ddimag) )*vacchi
335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356
            enddo
         enddo





      endif ! kvac1=kvac1_b and kvac2=kvac2_b
      enddo ! lprime
      enddo ! l

      enddo !symvac
c... cycle by the vacua finishes
 140  enddo      

      deallocate ( ac,bc,dt,dte,t,te,tei,u,ue,
     +             v,kvac1,kvac2,map2 )
      deallocate ( ac_b,bc_b,dt_b,dte_b,t_b,te_b,tei_b,u_b,ue_b,
     +             kvac1_b,kvac2_b,map2_b )

      END SUBROUTINE wann_mmkb_vac
      END MODULE m_wann_mmkb_vac