LCOV - code coverage report
Current view: top level - wannier - wann_mmkb_od_vac.F (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 156 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_od_vac
       8             : c**************************************************************
       9             : c      Determines the overlap matrix Mmn(k,k+b) in the vacuum
      10             : c      for the wannier functions. 
      11             : c      For more details see routine wannier.F and wann_mmk0_od_vac.F
      12             : c
      13             : c      Y. Mokrousov, F. Freimuth
      14             : c*************************************************************** 
      15             :       CONTAINS
      16           0 :       SUBROUTINE wann_mmkb_od_vac(
      17             :      >     DIMENSION,oneD,vacuum,stars,cell,
      18             :      >     vacchi,l_noco,nlotot,
      19             :      >     nbnd,z1,nmzxyd,nmzd,nv2d,k1d,k2d,k3d,n2d,n3d,
      20           0 :      >     ig,nmzxy,nmz,delz,ig2,n2d_1,
      21           0 :      >     bbmat,evac,evac_b,bkpt,bkpt_b,MM,vM,vz,vz_b,odi,
      22           0 :      >     nslibd,nslibd_b,jspin,jspin_b,k1,k2,k3,k1_b,k2_b,k3_b,
      23             :      >     jspd,nvd,area,nbasfcn,neigd,
      24           0 :      >     zMat,zMat_b,nv,nv_b,sk2,phi2,omtil,gb,qss,
      25             :      >     l_q, sign_q,
      26           0 :      <     mmn)
      27             :       use m_constants
      28             :       use m_types
      29             :       use m_od_abvac
      30             :       use m_cylbes
      31             :       use m_dcylbs
      32             :       use m_intgr, only : intgz0
      33             : 
      34             :       implicit none
      35             : 
      36             :       TYPE(t_dimension),INTENT(IN)   :: DIMENSION
      37             :       TYPE(t_oneD),INTENT(IN)        :: oneD
      38             :       TYPE(t_vacuum),INTENT(IN)      :: vacuum
      39             :       TYPE(t_stars),INTENT(IN)       :: stars
      40             :       TYPE(t_cell),INTENT(IN)        :: cell
      41             :       TYPE(t_mat), INTENT(IN)        :: zMat, zMat_b
      42             : 
      43             : c     .. scalar Arguments..
      44             :       logical, intent (in) :: l_noco
      45             :       integer, intent (in) :: nlotot
      46             :       integer, intent (in) :: nmzxyd,nmzd,nv2d,k1d,k2d,k3d,n3d
      47             :       integer, intent (in) :: nmzxy,nmz,MM,n2d,vM,nbnd
      48             :       integer, intent (in) :: nslibd,nslibd_b
      49             :       integer, intent (in) :: n2d_1,jspin,jspin_b,jspd,nvd
      50             :       integer, intent (in) :: nbasfcn,neigd
      51             :       real,    intent (in) :: delz,z1,evac,area,omtil,evac_b
      52             :       complex, intent (in) :: vacchi
      53             :       type (od_inp), intent (in) :: odi
      54             : 
      55             :       logical, intent (in) :: l_q
      56             :       integer, intent (in) :: sign_q
      57             : 
      58             : c     ..array arguments..
      59             :       real,    intent (in) :: bkpt(:),bkpt_b(:),qss(:) !bkpt(3),bkpt_b(3),qss(3)
      60             :       real,    intent (in) :: sk2(:),phi2(:) !sk2(n2d),phi2(n2d)
      61             :       integer, intent (in) :: ig(-k1d:,-k2d:,-k3d:) !ig(-k1d:k1d,-k2d:k2d,-k3d:k3d)
      62             :       integer, intent (in) :: ig2(:),nv(:),nv_b(:) !ig2(n3d),nv(jspd),nv_b(jspd)
      63             :       integer, intent (in) :: gb(:) !gb(3)
      64             :       real,    intent (in) :: vz(:),vz_b(:),bbmat(:,:) !vz(nmzd),bbmat(3,3)
      65             :       integer, intent (in) :: k1(:,:),k2(:,:),k3(:,:) !k1(nvd,jspd),k2(nvd,jspd),k3(nvd,jspd)
      66             :       integer, intent (in) :: k1_b(:,:),k2_b(:,:),k3_b(:,:) !k1_b(nvd,jspd),k2_b(nvd,jspd),k3_b(nvd,jspd)
      67             :       complex, intent (inout) :: mmn(:,:) !mmn(nbnd,nbnd)
      68             : 
      69             :       !logical, intent (in) :: l_q   ! dealing with q points?
      70             : 
      71             : c     ..basis wavefunctions in the vacuum
      72           0 :       real,    allocatable :: udz(:,:)
      73           0 :       real,    allocatable :: uz(:,:)
      74           0 :       real,    allocatable :: dudz(:,:)
      75           0 :       real,    allocatable :: duz(:,:)
      76           0 :       real,    allocatable :: u(:,:,:)
      77           0 :       real,    allocatable :: ud(:,:,:)
      78           0 :       real,    allocatable :: ddnv(:,:)
      79             : 
      80           0 :       real,    allocatable :: udz_b(:,:)
      81           0 :       real,    allocatable :: uz_b(:,:)
      82           0 :       real,    allocatable :: dudz_b(:,:)
      83           0 :       real,    allocatable :: duz_b(:,:)
      84           0 :       real,    allocatable :: u_b(:,:,:)
      85           0 :       real,    allocatable :: ud_b(:,:,:)
      86           0 :       real,    allocatable :: ddnv_b(:,:)
      87             : 
      88             : c     ..local scalars..
      89             :       logical tail
      90           0 :       real wronk,wronk1,arg,zks,tpi,uuo,udo,duo,ddo,zz,xx(nmz),zks0
      91             :       integer i,m,l,j,k,irec3,irec2,n,nv2,mp,ispin
      92             :       integer nv2_b,np1,lprime,addnoco,addnoco2
      93             :       complex avac,bvac,ic,phasfc
      94           0 :       complex, allocatable :: acof(:,:,:),bcof(:,:,:)
      95           0 :       integer, allocatable :: kvac3(:),map1(:)
      96           0 :       complex, allocatable :: acof_b(:,:,:),bcof_b(:,:,:)
      97           0 :       integer, allocatable :: kvac3_b(:),map1_b(:)
      98           0 :       real, allocatable :: bess(:),dbss(:)
      99           0 :       real, allocatable :: gbess(:,:),besss(:)
     100             :       REAL qssbti(3,2)
     101             : 
     102           0 :       addnoco=0
     103           0 :       addnoco2=0
     104           0 :       if(l_noco.and.(jspin.eq.2))then
     105           0 :          addnoco  = nv(1)   + nlotot
     106             :       endif
     107           0 :       if(l_noco.and.(jspin_b.eq.2))then
     108           0 :          addnoco2 = nv_b(1) + nlotot
     109             :       endif
     110             : 
     111             :       allocate ( udz(nv2d,-vM:vM),uz(nv2d,-vM:vM),dudz(nv2d,-vM:vM),
     112             :      +           duz(nv2d,-vM:vM),u(nmzd,nv2d,-vM:vM),
     113             :      +           ud(nmzd,nv2d,-vM:vM),ddnv(nv2d,-vM:vM),
     114             :      +           udz_b(nv2d,-vM:vM),uz_b(nv2d,-vM:vM),
     115             :      +           dudz_b(nv2d,-vM:vM),
     116             :      +           duz_b(nv2d,-vM:vM),u_b(nmzd,nv2d,-vM:vM),
     117             :      +           ud_b(nmzd,nv2d,-vM:vM),ddnv_b(nv2d,-vM:vM),
     118             :      +           bess(-odi%mb:odi%mb),dbss(-odi%mb:odi%mb),
     119             :      +           besss(-odi%M:odi%M),gbess(-odi%M:odi%M,nmzd),
     120             :      +           acof(nv2d,-odi%mb:odi%mb,nslibd),
     121             :      +           bcof(nv2d,-odi%mb:odi%mb,nslibd),
     122             :      +           acof_b(nv2d,-odi%mb:odi%mb,nslibd_b),
     123             :      +           bcof_b(nv2d,-odi%mb:odi%mb,nslibd_b),
     124             :      +           kvac3(nv2d),map1(nvd), 
     125           0 :      +           kvac3_b(nv2d),map1_b(nvd) )
     126             : 
     127           0 :       tpi = 2 * pimach() ; ic = cmplx(0.,1.)
     128             : 
     129           0 :       tail = .true.
     130           0 :       np1 = nmz + 1
     131             : 
     132           0 :       acof(:,:,:) = cmplx(0.,0.) ; bcof(:,:,:) = cmplx(0.,0.)
     133           0 :       acof_b(:,:,:) = cmplx(0.,0.) ; bcof_b(:,:,:) = cmplx(0.,0.)
     134             : 
     135           0 :       nv2 = 0 ; nv2_b = 0
     136             : 
     137           0 :       do 20 k = 1,nv(jspin)
     138           0 :          do 10 j = 1,nv2
     139           0 :             if (k3(k,jspin).eq.kvac3(j)) then
     140           0 :                map1(k) = j
     141           0 :                goto 20
     142             :             endif
     143           0 :  10      continue
     144           0 :          nv2 = nv2 + 1
     145           0 :          if (nv2.gt.nv2d) stop 'nv2d'
     146           0 :          kvac3(nv2) = k3(k,jspin)
     147           0 :          map1(k) = nv2
     148           0 :  20   continue
     149             : 
     150           0 :       do 21 k = 1,nv_b(jspin_b)
     151           0 :          do 11 j = 1,nv2_b
     152           0 :             if (k3_b(k,jspin_b).eq.kvac3_b(j)) then
     153           0 :                map1_b(k) = j
     154           0 :                goto 21
     155             :             endif
     156           0 :  11      continue
     157           0 :          nv2_b = nv2_b + 1
     158           0 :          if (nv2_b.gt.nv2d) stop 'nv2d'
     159           0 :          kvac3_b(nv2_b) = k3_b(k,jspin_b)
     160           0 :          map1_b(k) = nv2_b
     161           0 :  21   continue
     162             : 
     163           0 :       wronk = 2.0
     164             : 
     165             : c...for the k-point
     166             : 
     167           0 :       qssbti(1,1) = - qss(1)/2.     
     168           0 :       qssbti(2,1) = - qss(2)/2.     
     169           0 :       qssbti(1,2) = + qss(1)/2.     
     170           0 :       qssbti(2,2) = + qss(2)/2.
     171           0 :       qssbti(3,1) = - qss(3)/2.
     172           0 :       qssbti(3,2) = + qss(3)/2.
     173           0 :       DO ispin = 1,1 ! jspins
     174             :       CALL od_abvac(
     175             :      >      cell,vacuum,DIMENSION,stars,oneD,
     176             :      >      qssbti(3,jspin),odi%n2d,
     177             :      >      wronk,evac,bkpt,odi%M,odi%mb,
     178             :      >      vz,kvac3,nv2,
     179             :      <      uz(1,-vM),duz(1,-vM),u(1,1,-vM),udz(1,-vM),
     180           0 :      <      dudz(1,-vM),ddnv(1,-vM),ud(1,1,-vM))
     181             :       ENDDO
     182             : 
     183           0 :       do k = 1,nv(jspin)
     184           0 :          l = map1(k)
     185           0 :          irec3 = ig(k1(k,jspin),k2(k,jspin),k3(k,jspin))
     186           0 :          if (irec3.ne.0) then
     187           0 :             irec2 = ig2(irec3)
     188           0 :             zks = sk2(irec2)*z1
     189           0 :             arg = phi2(irec2)
     190           0 :             call cylbes(odi%mb,zks,bess)
     191           0 :             call dcylbs(odi%mb,zks,bess,dbss)
     192           0 :             do m = -odi%mb,odi%mb
     193             :                wronk1 = uz(l,m)*dudz(l,m) -
     194           0 :      -              udz(l,m)*duz(l,m)
     195             :                avac = exp(-cmplx(0.0,m*arg))*(ic**m)*
     196             :      *              cmplx(dudz(l,m)*bess(m) -
     197             :      +              udz(l,m)*sk2(irec2)*dbss(m),0.0)/
     198           0 :      /              ((wronk1)*sqrt(omtil))
     199             :                bvac = exp(-cmplx(0.0,m*arg))*(ic**m)*
     200             :      *              cmplx(-duz(l,m)*bess(m) +
     201             :      -              uz(l,m)*sk2(irec2)*dbss(m),0.0)/
     202           0 :      /              ((wronk1)*sqrt(omtil))
     203           0 :                IF (zMat%l_real) THEN
     204           0 :                   do n = 1,nslibd
     205             :                       acof(l,m,n) = acof(l,m,n) +
     206           0 :      +                   zMat%data_r(k+addnoco,n)*avac
     207             : c     +                    conjg(zMat%data_r(k,n))*avac
     208             :                       bcof(l,m,n) = bcof(l,m,n) +
     209           0 :      +                   zMat%data_r(k+addnoco,n)*bvac
     210             : c     +                    conjg(zMat%data_r(k,n))*bvac
     211             :                   enddo
     212             :                ELSE
     213           0 :                   do n = 1,nslibd
     214             :                       acof(l,m,n) = acof(l,m,n) +
     215           0 :      +                   zMat%data_c(k+addnoco,n)*avac
     216             : c     +                    conjg(zMat%data_c(k,n))*avac
     217             :                       bcof(l,m,n) = bcof(l,m,n) +
     218           0 :      +                   zMat%data_c(k+addnoco,n)*bvac
     219             : c     +                    conjg(zMat%data_c(k,n))*bvac
     220             :                   enddo
     221             :                END IF
     222             :             enddo      ! -mb:mb
     223             :          endif
     224             :       enddo
     225             : 
     226             : c...for the b-point
     227             : 
     228           0 :       DO ispin = 1,1 ! jspins
     229             :       call od_abvac(
     230             :      >      cell,vacuum,DIMENSION,stars,oneD,
     231             :      >      qssbti(3,jspin_b),odi%n2d,
     232             :      >      wronk,evac_b,bkpt_b,odi%M,odi%mb,
     233             :      >      vz_b,kvac3_b,nv2_b,
     234             :      <      uz_b(1,-vM),duz_b(1,-vM),u_b(1,1,-vM),udz_b(1,-vM),
     235           0 :      <      dudz_b(1,-vM),ddnv_b(1,-vM),ud_b(1,1,-vM))
     236             :       ENDDO
     237             : 
     238           0 :       do k = 1,nv_b(jspin_b)
     239           0 :          l = map1_b(k)
     240           0 :          irec3 = ig(k1_b(k,jspin_b),k2_b(k,jspin_b),k3_b(k,jspin_b))
     241           0 :          if (irec3.ne.0) then
     242           0 :             irec2 = ig2(irec3)
     243           0 :             zks = sk2(irec2)*z1
     244           0 :             arg = phi2(irec2)
     245           0 :             call cylbes(odi%mb,zks,bess)
     246           0 :             call dcylbs(odi%mb,zks,bess,dbss)
     247           0 :             do m = -odi%mb,odi%mb
     248             :                wronk1 = uz_b(l,m)*dudz_b(l,m) -
     249           0 :      -              udz_b(l,m)*duz_b(l,m)
     250             :                avac = exp(-cmplx(0.0,m*arg))*(ic**m)*
     251             :      *              cmplx(dudz_b(l,m)*bess(m) -
     252             :      +              udz_b(l,m)*sk2(irec2)*dbss(m),0.0)/
     253           0 :      /              ((wronk1)*sqrt(omtil))
     254             :                bvac = exp(-cmplx(0.0,m*arg))*(ic**m)*
     255             :      &              cmplx(-duz_b(l,m)*bess(m) +
     256             :      -              uz_b(l,m)*sk2(irec2)*dbss(m),0.0)/
     257           0 :      /              ((wronk1)*sqrt(omtil))
     258           0 :                IF (zMat_b%l_real) THEN
     259           0 :                   do n = 1,nslibd_b
     260             :                       acof_b(l,m,n) = acof_b(l,m,n) +
     261           0 :      +                   zMat_b%data_r(k+addnoco2,n)*avac
     262             : c     +                    conjg(zMat_b%data_r(k,n))*avac
     263             :                       bcof_b(l,m,n) = bcof_b(l,m,n) +
     264           0 :      +                   zMat_b%data_r(k+addnoco2,n)*bvac
     265             : c     +                    conjg(zMat_b%data_r(k,n))*bvac
     266             :                   enddo
     267             :                ELSE
     268           0 :                   do n = 1,nslibd_b
     269             :                       acof_b(l,m,n) = acof_b(l,m,n) +
     270           0 :      +                   zMat_b%data_c(k+addnoco2,n)*avac
     271             : c     +                    conjg(zMat_b%data_c(k,n))*avac
     272             :                       bcof_b(l,m,n) = bcof_b(l,m,n) +
     273           0 :      +                   zMat_b%data_c(k+addnoco2,n)*bvac
     274             : c     +                    conjg(zMat_b%data_c(k,n))*bvac
     275             :                   enddo
     276             :                END IF
     277             :             enddo      ! -mb:mb
     278             :          endif 
     279             :       enddo          ! k = 1,nv
     280             : 
     281             : c  now actually computing the Mmn matrix
     282             : 
     283           0 :       irec3 = ig(gb(1),gb(2),gb(3))
     284           0 :       if (irec3.eq.0) stop 'Gb is not in the list of Gs'
     285           0 :       irec2 = ig2(irec3)
     286           0 :       zks0 = sk2(irec2)
     287           0 :       arg = phi2(irec2)
     288             : 
     289             :       !write(*,*)'mmkb_od_vac zks0',zks0
     290             :       
     291           0 :       gbess(:,:) = 0.
     292           0 :       do i = 1,nmz
     293           0 :          zz = z1+(i-1)*delz
     294           0 :          zks = zks0*zz
     295           0 :          besss(:) = 0.
     296           0 :          call cylbes(odi%M,zks,besss)
     297           0 :          do m = -odi%M,odi%M
     298           0 :             gbess(m,i) = besss(m)
     299             : c           gbess(i) = 1.
     300             :          enddo
     301             :       enddo
     302             : 
     303           0 :       do l = 1,nv2
     304           0 :       do lprime = 1,nv2_b
     305             : 
     306           0 :       if (kvac3(l).eq.(kvac3_b(lprime)-gb(3))) then      
     307             : 
     308           0 :        do m = -odi%mb,odi%mb
     309           0 :         do mp = -odi%mb,odi%mb
     310             : 
     311           0 :          phasfc = exp(-cmplx(0.0,(mp-m)*arg))
     312           0 :          if(l_q) phasfc =  (sign_q)**(mp-m) * phasfc
     313             : 
     314           0 :          do i = 1,nmz
     315           0 :            zz = z1+(i-1)*delz
     316           0 :            xx(np1-i) = zz*u(i,l,m)*u_b(i,lprime,mp)*gbess(mp-m,i)
     317             :          enddo    
     318           0 :          call intgz0(xx,delz,nmz,uuo,tail)
     319             : 
     320           0 :          do i = 1,nmz
     321           0 :            zz = z1+(i-1)*delz
     322           0 :            xx(np1-i) = zz*u(i,l,m)*ud_b(i,lprime,mp)*gbess(mp-m,i)
     323             :          enddo    
     324           0 :          call intgz0(xx,delz,nmz,udo,tail)
     325             : 
     326           0 :          do i = 1,nmz
     327           0 :            zz = z1+(i-1)*delz
     328           0 :            xx(np1-i) = zz*ud(i,l,m)*u_b(i,lprime,mp)*gbess(mp-m,i)
     329             :          enddo     
     330           0 :          call intgz0(xx,delz,nmz,duo,tail)
     331             : 
     332           0 :          do i = 1,nmz
     333           0 :            zz = z1+(i-1)*delz
     334           0 :            xx(np1-i) = zz*ud(i,l,m)*ud_b(i,lprime,mp)*gbess(mp-m,i)
     335             :          enddo   
     336           0 :          call intgz0(xx,delz,nmz,ddo,tail)
     337             : 
     338           0 :          do i = 1,nslibd
     339           0 :            do j = 1,nslibd_b
     340             :             mmn(i,j) = mmn(i,j) + phasfc*area*(
     341             :      *          acof(l,m,i)*conjg(acof_b(lprime,mp,j))*uuo +
     342             :      +          acof(l,m,i)*conjg(bcof_b(lprime,mp,j))*udo +
     343             :      +          bcof(l,m,i)*conjg(acof_b(lprime,mp,j))*duo +
     344           0 :      +          bcof(l,m,i)*conjg(bcof_b(lprime,mp,j))*ddo )*vacchi  
     345             :            enddo 
     346             :          enddo
     347             : 
     348             :         enddo !mp
     349             :        enddo  !m
     350             : 
     351             :       endif  ! kvac(k)=kvac(k+b)
     352             : 
     353             :       enddo  ! lprime
     354             :       enddo  ! l
     355             : 
     356           0 :       deallocate ( udz,uz,dudz,duz,u,ud,ddnv,bess,dbss,acof,bcof )
     357           0 :       deallocate ( udz_b,uz_b,dudz_b,duz_b,u_b,ud_b,ddnv_b,gbess,besss )
     358           0 :       deallocate ( acof_b,bcof_b )
     359             : 
     360           0 :       END SUBROUTINE wann_mmkb_od_vac
     361             :       END MODULE m_wann_mmkb_od_vac

Generated by: LCOV version 1.13