LCOV - code coverage report
Current view: top level - wannier - wann_mmkb_vac.F (source / functions) Hit Total Coverage
Test: FLEUR test coverage Lines: 0 172 0.0 %
Date: 2024-03-28 04:22:06 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_vac
       8             :       use m_juDFT
       9             : c**************************************************************
      10             : c      Determines the overlap matrix Mmn(k) in the vacuum
      11             : c      in the film case for the wannier functions.
      12             : c      For more details see routine wannier.F 
      13             : c
      14             : c      Y. Mokrousov, F. Freimuth
      15             : c*************************************************************** 
      16             :       CONTAINS
      17           0 :       SUBROUTINE wann_mmkb_vac(
      18             :      >     vacchi,l_noco,nlotot,qss,
      19             :      >     nbnd,z1,nmzd,nv2d,k1d,k2d,k3d,n3d,nvac,
      20             :      >     ig,nmz,delz,ig2,area,bmat,
      21           0 :      >     bbmat,evac,evac_b,bkpt,bkpt_b,vz,vz_b,
      22             :      >     nslibd,nslibd_b,jspin,jspin_b,
      23           0 :      >     k1,k2,k3,k1_b,k2_b,k3_b,jspd,nvd,
      24           0 :      >     nbasfcn,neigd,zMat,zMat_b,nv,nv_b,omtil,gb,
      25           0 :      <     mmn)
      26             : 
      27             :       USE m_types
      28             :       use m_constants
      29             :       use m_intgr, only : intgz0
      30             :       USE m_vacuz
      31             :       USE m_vacudz
      32             : 
      33             :       implicit none
      34             : 
      35             :       TYPE(t_mat), INTENT(IN) :: zMat, zMat_b
      36             : 
      37             : c     .. scalar Arguments..
      38             :       logical, intent (in) :: l_noco
      39             :       integer, intent (in) :: nlotot,jspin_b
      40             :       real,    intent (in) :: qss(3)
      41             :       integer, intent (in) :: nmzd,nv2d,k1d,k2d,k3d,n3d,nbnd
      42             :       integer, intent (in) :: nmz,nslibd,nslibd_b,nvac
      43             :       integer, intent (in) :: jspin,jspd,nvd
      44             :       integer, intent (in) :: nbasfcn,neigd
      45             :       integer, intent (in) :: gb(3)
      46             :       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             : !      integer, intent (in) :: ig(-k1d:k1d,-k2d:k2d,-k3d:k3d)
      52             :       integer, intent (in) :: ig(-k1d:,-k2d:,-k3d:)
      53             : !      integer, intent (in) :: ig2(n3d)
      54             :       integer, intent (in) :: ig2(:)!ig2(n3d)      
      55             : !      integer, intent (in) :: nv(jspd),nv_b(jspd)
      56             :       integer, intent (in) :: nv(:),nv_b(:)      
      57             :       real,    intent (in) :: vz(nmzd,2),vz_b(nmzd,2)
      58             :       real,    intent (in) :: bbmat(3,3),bmat(3,3)
      59             : !      integer, intent (in) :: k1(nvd,jspd),k2(nvd,jspd),k3(nvd,jspd)
      60             :       integer, intent (in) :: k1(:,:),k2(:,:),k3(:,:)      
      61             : !      integer,intent(in)::k1_b(nvd,jspd),k2_b(nvd,jspd),k3_b(nvd,jspd)
      62             :       integer,intent(in)::k1_b(:,:),k2_b(:,:),k3_b(:,:)      
      63             : !      complex, intent (inout) :: mmn(nbnd,nbnd)
      64             :       complex, intent (inout) :: mmn(:,:)
      65             : 
      66             : c     ..basis wavefunctions in the vacuum
      67           0 :       complex, allocatable :: ac(:,:),bc(:,:),ac_b(:,:),bc_b(:,:)
      68           0 :       real,    allocatable :: dt(:),dte(:)
      69           0 :       real,    allocatable :: t(:),te(:),tei(:)
      70           0 :       real,    allocatable :: u(:,:),ue(:,:),v(:)
      71           0 :       real,    allocatable :: dt_b(:),dte_b(:)
      72           0 :       real,    allocatable :: t_b(:),te_b(:),tei_b(:)
      73           0 :       real,    allocatable :: u_b(:,:),ue_b(:,:)
      74             : 
      75             : c     ..local scalars..
      76             :       logical tail
      77             :       real wronk,arg,zks,tpi,vz0(2),vz0_b(2),scale,evacp,ev,const
      78           0 :       real uu,ud,du,dd,xx(nmz),xximag(nmz)
      79             :       real :: uuimag,udimag,duimag,ddimag,qss1,qss2
      80             :       integer i,m,l,j,k,n,nv2,nv2_b,ivac,n2,n2_b,sign,ik
      81             :       integer :: lprime,np1,addnoco,addnoco2
      82             :       complex :: av,bv,ic,c_1
      83           0 :       integer, allocatable :: kvac1(:),kvac2(:),map2(:)
      84           0 :       integer, allocatable :: kvac1_b(:),kvac2_b(:),map2_b(:)
      85             :       integer symvac,symvacvac
      86           0 :       call timestart("wann_mmkb_vac")
      87             :       allocate ( ac(nv2d,nslibd),bc(nv2d,nslibd),
      88             :      +           ac_b(nv2d,nslibd_b),bc_b(nv2d,nslibd_b),
      89             :      +           dt(nv2d),dte(nv2d),t(nv2d),te(nv2d),
      90             :      +           tei(nv2d),u(nmzd,nv2d),ue(nmzd,nv2d),
      91             :      +           dt_b(nv2d),dte_b(nv2d),t_b(nv2d),te_b(nv2d),
      92             :      +           tei_b(nv2d),u_b(nmzd,nv2d),ue_b(nmzd,nv2d),
      93             :      +           v(3),kvac1(nv2d),kvac2(nv2d),map2(nvd),
      94           0 :      +           kvac1_b(nv2d),kvac2_b(nv2d),map2_b(nvd) )
      95             : 
      96           0 :       tpi = 2 * pimach() ; ic = cmplx(0.,1.)
      97             : 
      98           0 :       tail = .true.
      99           0 :       np1 = nmz + 1
     100             : 
     101             : c.. determining the indexing array (in-plane stars)
     102             : c.. for the k-point
     103             : 
     104           0 :       wronk = 2.0
     105           0 :       const = 1.0 / ( sqrt(omtil)*wronk )
     106             : 
     107           0 :       do ivac = 1,2
     108           0 :          vz0(ivac) = vz(nmz,ivac)
     109           0 :          vz0_b(ivac) = vz_b(nmz,ivac)
     110             :       enddo
     111             : 
     112           0 :       n2 = 0 ; n2_b = 0
     113             : 
     114           0 :       addnoco=0
     115           0 :       addnoco2=0
     116           0 :       if(l_noco.and.jspin.eq.2)then
     117           0 :         addnoco=nv(1)+nlotot
     118             :       endif
     119             :       if(l_noco.and.jspin_b.eq.2)then
     120             :         addnoco2=nv_b(1)+nlotot
     121             :       endif
     122             : 
     123           0 :       do 40 k = 1,nv(jspin)
     124           0 :          do 30 j = 1,n2
     125           0 :             if ( k1(k,jspin).eq.kvac1(j) .and.
     126             :      +          k2(k,jspin).eq.kvac2(j) ) then
     127           0 :                 map2(k) = j
     128           0 :                 goto 40
     129             :              endif 
     130           0 :  30      continue
     131           0 :          n2 = n2 + 1
     132             :          
     133           0 :          IF(n2>nv2d) then
     134           0 :             write(*,*)n2,nv2d,'jspin',jspin
     135             :          endif
     136             :  
     137           0 :          IF (n2>nv2d)  CALL juDFT_error("wannier Mmn vac",calledby
     138           0 :      +        ="wann_mmkb_vac")
     139             : 
     140           0 :          kvac1(n2) = k1(k,jspin)
     141           0 :          kvac2(n2) = k2(k,jspin)
     142           0 :          map2(k) = n2
     143           0 :  40   continue
     144             :          !write(*,*)'ok',n2,nv2d,'jspin',jspin
     145             : 
     146             : c.. and for the b-point
     147             :  
     148           0 :       do 41 k = 1,nv_b(jspin_b)
     149           0 :          do 31 j = 1,n2_b
     150           0 :             if ( k1_b(k,jspin_b).eq.kvac1_b(j) .and.
     151             :      +          k2_b(k,jspin_b).eq.kvac2_b(j) ) then
     152           0 :                 map2_b(k) = j
     153           0 :                 goto 41
     154             :              endif
     155           0 :  31      continue
     156           0 :          n2_b = n2_b + 1
     157             :      
     158           0 :          IF(n2_b>nv2d) then
     159           0 :             write(*,*)n2_b,nv2d,'jspin_b',jspin_b
     160             :          endif
     161             : 
     162           0 :          IF (n2_b>nv2d)  CALL juDFT_error("wannier Mmn vac",calledby
     163           0 :      +        ="wann_mmkb_vac")
     164             : 
     165           0 :          kvac1_b(n2_b) = k1_b(k,jspin_b)
     166           0 :          kvac2_b(n2_b) = k2_b(k,jspin_b)
     167           0 :          map2_b(k) = n2_b
     168           0 :  41   continue
     169             :          !write(*,*)'ok',n2_b,nv2d,'jspin_b',jspin_b
     170             : 
     171             : c...cycle by the vacua
     172           0 :       do 140 ivac = 1,nvac
     173             : 
     174             : 
     175             : 
     176           0 :       sign = 3. - 2.*ivac
     177           0 :       evacp = evac(ivac)
     178             : 
     179           0 :       nv2 = n2 ; nv2_b = n2_b
     180             : 
     181             : c.. the body of the routine
     182             : 
     183           0 :       qss1=0.0
     184           0 :       qss2=0.0
     185           0 :       if(l_noco.and.jspin.eq.1)then
     186           0 :         qss1=-qss(1)/2.0
     187           0 :         qss2=-qss(2)/2.0
     188           0 :       elseif(l_noco.and.jspin.eq.2)then
     189           0 :         qss1=qss(1)/2.0
     190           0 :         qss2=qss(2)/2.0
     191             :       endif
     192             : 
     193           0 :       do ik = 1,nv2
     194           0 :          v(1) = bkpt(1) + kvac1(ik) + qss1
     195           0 :          v(2) = bkpt(2) + kvac2(ik) + qss2
     196           0 :          v(3) = 0.
     197           0 :          ev = evacp - 0.5*dot_product(v,matmul(bbmat,v))
     198             :          call vacuz(ev,vz(1,ivac),vz0(ivac),nmz,delz,t(ik),dt(ik),
     199           0 :      +        u(1,ik))
     200             :          call vacudz(ev,vz(1,ivac),vz0(ivac),nmz,delz,te(ik),
     201             :      +        dte(ik),tei(ik),ue(1,ik),dt(ik),
     202           0 :      +        u(1,ik))
     203           0 :          scale = wronk/ (te(ik)*dt(ik)-dte(ik)*t(ik))
     204           0 :          te(ik) = scale*te(ik)
     205           0 :          dte(ik) = scale*dte(ik)
     206           0 :          tei(ik) = scale*tei(ik)
     207           0 :          do j = 1,nmz
     208           0 :             ue(j,ik) = scale*ue(j,ik)
     209             :          enddo
     210             :       enddo
     211             : c-----> construct a and b coefficients for the k-point
     212           0 :        symvacvac=1
     213           0 :        if (nvac==1) symvacvac=2
     214           0 :        do symvac=1,symvacvac
     215           0 :          do i = 1,nv2d
     216           0 :             do n = 1,nslibd
     217           0 :                ac(i,n) = cmplx(0.0,0.0)
     218           0 :                bc(i,n) = cmplx(0.0,0.0)
     219             :             enddo   
     220           0 :             do n = 1,nslibd_b
     221           0 :                ac_b(i,n) = cmplx(0.0,0.0)
     222           0 :                bc_b(i,n) = cmplx(0.0,0.0)
     223             :             enddo   
     224             :          enddo   
     225             : 
     226           0 :         if (symvac==2) sign=-1.0   
     227             : 
     228           0 :       do k = 1,nv(jspin)
     229           0 :          l = map2(k)
     230           0 :          zks = k3(k,jspin)*bmat(3,3)*sign
     231           0 :          arg = zks*z1
     232           0 :          c_1 = cmplx(cos(arg),sin(arg)) * const
     233           0 :          av = -c_1 * cmplx( dte(l),zks*te(l) )
     234           0 :          bv =  c_1 * cmplx(  dt(l),zks* t(l) )
     235             : c-----> loop over basis functions
     236           0 :          IF(zMat%l_real) THEN
     237           0 :             do n = 1,nslibd
     238           0 :                ac(l,n) = ac(l,n) + zMat%data_r(k+addnoco,n)*av
     239           0 :                bc(l,n) = bc(l,n) + zMat%data_r(k+addnoco,n)*bv
     240             :             enddo
     241             :          ELSE
     242           0 :             do n = 1,nslibd
     243           0 :                ac(l,n) = ac(l,n) + zMat%data_c(k+addnoco,n)*av
     244           0 :                bc(l,n) = bc(l,n) + zMat%data_c(k+addnoco,n)*bv
     245             :             enddo
     246             :          END IF
     247             :       enddo
     248             : 
     249             : c...now for the k+b point
     250             : 
     251           0 :       evacp = evac_b(ivac)
     252           0 :       do ik = 1,nv2_b
     253           0 :          v(1) = bkpt_b(1) + kvac1_b(ik) + qss1
     254           0 :          v(2) = bkpt_b(2) + kvac2_b(ik) + qss2
     255           0 :          v(3) = 0.
     256           0 :          ev = evacp - 0.5*dot_product(v,matmul(bbmat,v))
     257             :          call vacuz(ev,vz_b(1,ivac),vz0_b(ivac),nmz,delz,t_b(ik),
     258           0 :      +        dt_b(ik),u_b(1,ik))
     259             :          call vacudz(ev,vz_b(1,ivac),vz0_b(ivac),nmz,delz,te_b(ik),
     260             :      +        dte_b(ik),tei_b(ik),ue_b(1,ik),dt_b(ik),
     261           0 :      +        u_b(1,ik))
     262           0 :          scale = wronk/ (te_b(ik)*dt_b(ik)-dte_b(ik)*t_b(ik))
     263           0 :          te_b(ik) = scale*te_b(ik)
     264           0 :          dte_b(ik) = scale*dte_b(ik)
     265           0 :          tei_b(ik) = scale*tei_b(ik)
     266           0 :          do j = 1,nmz
     267           0 :             ue_b(j,ik) = scale*ue_b(j,ik)
     268             :          enddo
     269             :       enddo
     270             : c-----> construct a and b coefficients for the k+b point
     271             : 
     272           0 :       do k = 1,nv_b(jspin_b)
     273           0 :          l = map2_b(k)
     274           0 :          zks = k3_b(k,jspin_b)*bmat(3,3)*sign
     275           0 :          arg = zks*z1
     276           0 :          c_1 = cmplx(cos(arg),sin(arg)) * const
     277           0 :          av = -c_1 * cmplx( dte_b(l),zks*te_b(l) )
     278           0 :          bv =  c_1 * cmplx( dt_b(l),zks*t_b(l) )
     279             : c-----> loop over basis functions
     280           0 :          IF(zMat_b%l_real) THEN
     281           0 :             do n = 1,nslibd_b
     282           0 :                ac_b(l,n) = ac_b(l,n) + zMat_b%data_r(k+addnoco,n)*av
     283           0 :                bc_b(l,n) = bc_b(l,n) + zMat_b%data_r(k+addnoco,n)*bv
     284             :             enddo
     285             :          ELSE
     286           0 :             do n = 1,nslibd_b
     287           0 :                ac_b(l,n) = ac_b(l,n) + zMat_b%data_c(k+addnoco,n)*av
     288           0 :                bc_b(l,n) = bc_b(l,n) + zMat_b%data_c(k+addnoco,n)*bv
     289             :             enddo
     290             :          END IF
     291             :       enddo
     292             : 
     293             : 
     294           0 :       do l = 1,nv2
     295           0 :       do lprime = 1,nv2_b
     296             :       if (kvac1(l).eq.(kvac1_b(lprime)-gb(1))
     297           0 :      & .and. kvac2(l).eq.(kvac2_b(lprime)-gb(2)))then
     298           0 :          zks = gb(3)*bmat(3,3)*sign
     299             : 
     300           0 :          do i = 1,nmz
     301             :              xx(np1-i) = u(i,l)*u_b(i,lprime)*
     302           0 :      *          cos((z1+(i-1)*delz)*zks)
     303             :             xximag(np1-i) = u(i,l)*u_b(i,lprime)*
     304           0 :      *          sin((z1+(i-1)*delz)*zks)
     305             :          enddo   
     306           0 :          call intgz0(xx,delz,nmz,uu,tail)
     307           0 :          call intgz0(xximag,delz,nmz,uuimag,tail)
     308             : 
     309           0 :          do i = 1,nmz
     310             :             xx(np1-i) = u(i,l)*ue_b(i,lprime)*
     311           0 :      *   cos((z1+(i-1)*delz)*zks)
     312             :             xximag(np1-i) = u(i,l)*ue_b(i,lprime)*
     313           0 :      *          sin((z1+(i-1)*delz)*zks)
     314             :          enddo   
     315           0 :          call intgz0(xx,delz,nmz,ud,tail)
     316           0 :          call intgz0(xximag,delz,nmz,udimag,tail)
     317             : 
     318           0 :          do i = 1,nmz
     319             :             xx(np1-i) = ue(i,l)*u_b(i,lprime)*
     320           0 :      *   cos((z1+(i-1)*delz)*zks)
     321             :             xximag(np1-i) = ue(i,l)*u_b(i,lprime)*
     322           0 :      *          sin((z1+(i-1)*delz)*zks)
     323             :          enddo   
     324           0 :          call intgz0(xx,delz,nmz,du,tail)
     325           0 :          call intgz0(xximag,delz,nmz,duimag,tail)
     326           0 :          do i = 1,nmz
     327             :             xx(np1-i) = ue(i,l)*ue_b(i,lprime)*
     328           0 :      *   cos((z1+(i-1)*delz)*zks)
     329             :             xximag(np1-i) = ue(i,l)*ue_b(i,lprime)*
     330           0 :      *          sin((z1+(i-1)*delz)*zks)
     331             :          enddo   
     332           0 :          call intgz0(xx,delz,nmz,dd,tail)
     333           0 :          call intgz0(xximag,delz,nmz,ddimag,tail)
     334             : 
     335           0 :          do i = 1,nslibd
     336           0 :             do j = 1,nslibd_b
     337             :                mmn(i,j) = mmn(i,j) + area*(
     338             :      *  ac(l,i)*conjg(ac_b(lprime,j))*cmplx(uu,uuimag) +
     339             :      +  ac(l,i)*conjg(bc_b(lprime,j))*cmplx(ud,udimag) +
     340             :      *  bc(l,i)*conjg(ac_b(lprime,j))*cmplx(du,duimag) +
     341           0 :      +  bc(l,i)*conjg(bc_b(lprime,j))*cmplx(dd,ddimag) )*vacchi
     342             :             enddo
     343             :          enddo
     344             : 
     345             : 
     346             : 
     347             : 
     348             : 
     349             :       endif ! kvac1=kvac1_b and kvac2=kvac2_b
     350             :       enddo ! lprime
     351             :       enddo ! l
     352             : 
     353             :       enddo !symvac
     354             : c... cycle by the vacua finishes
     355           0 :  140  enddo      
     356             : 
     357           0 :       deallocate ( ac,bc,dt,dte,t,te,tei,u,ue,
     358           0 :      +             v,kvac1,kvac2,map2 )
     359           0 :       deallocate ( ac_b,bc_b,dt_b,dte_b,t_b,te_b,tei_b,u_b,ue_b,
     360           0 :      +             kvac1_b,kvac2_b,map2_b )
     361             : 
     362           0 :       call timestop("wann_mmkb_vac")
     363           0 :       END SUBROUTINE wann_mmkb_vac
     364             :       END MODULE m_wann_mmkb_vac

Generated by: LCOV version 1.14