LCOV - code coverage report
Current view: top level - wannier/uhu - wann_uHu_od_vac.F (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 267 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 1 0.0 %

          Line data    Source code
       1             : c***************************************c
       2             : c   Vacuum contribution to uHu matrix   c
       3             : c   in one-dimensional FLAPW mode       c
       4             : c***************************************c
       5             : c   < u_{k+b1} | H_{k} | u_{k+b2} >     c
       6             : c***************************************c
       7             : c              J.-P. Hanke, Dec. 2015   c
       8             : c***************************************c
       9             :       MODULE m_wann_uHu_od_vac
      10             :       CONTAINS
      11           0 :       SUBROUTINE wann_uHu_od_vac(
      12             :      >     DIMENSION,oneD,vacuum,stars,cell,
      13             :      >     chi,l_noco,l_soc,jspins,nlotot,
      14             :      >     nbnd,z1,nmzxyd,nmzd,nv2d,k1d,k2d,k3d,n2d,n3d,
      15           0 :      >     ig,nmzxy,nmz,delz,ig2,
      16           0 :      >     bbmat,evac,evac_b,bkpt,bkpt_b,odi,vxy,vz,
      17             :      >     nslibd,nslibd_b,jspin,jspin_b,ico,
      18           0 :      >     k1,k2,k3,k1_b,k2_b,k3_b,
      19             :      >     jspd,nvd,area,nbasfcn,neigd,
      20           0 :      >     zMat,zMat_b,nv,nv_b,sk2,phi2,omtil,gb,gb2,qss,sign2,
      21           0 :      <     uHu)
      22             : 
      23             :       use m_constants
      24             :       use m_types
      25             :       use m_od_abvac
      26             :       use m_cylbes
      27             :       use m_dcylbs
      28             :       use m_intgr, only : intgz0
      29             :       use m_d2fdz2
      30             : 
      31             :       implicit none
      32             : 
      33             :       TYPE(t_dimension),INTENT(IN)   :: DIMENSION
      34             :       TYPE(t_oneD),INTENT(IN)        :: oneD
      35             :       TYPE(t_vacuum),INTENT(IN)      :: vacuum
      36             :       TYPE(t_stars),INTENT(IN)       :: stars
      37             :       TYPE(t_cell),INTENT(IN)        :: cell
      38             : 
      39             :       TYPE(t_zmat), INTENT(IN) :: zMat, zMat_b
      40             : 
      41             : c     .. scalar Arguments..
      42             :       logical, intent (in) :: l_noco,l_soc
      43             :       integer, intent (in) :: jspins,nlotot,ico
      44             :       integer, intent (in) :: nmzxyd,nmzd,nv2d,k1d,k2d,k3d,n3d
      45             :       integer, intent (in) :: nmzxy,nmz,n2d,nbnd
      46             :       integer, intent (in) :: nslibd,nslibd_b
      47             :       integer, intent (in) :: jspin,jspin_b,jspd,nvd
      48             :       integer, intent (in) :: nbasfcn,neigd
      49             :       real,    intent (in) :: delz,z1,evac,area,omtil,evac_b
      50             :       complex, intent (in) :: chi
      51             :       type (od_inp), intent (in) :: odi
      52             : 
      53             : 
      54             : c     ..array arguments..
      55             :       real,    intent (in) :: bkpt(:),bkpt_b(:),qss(:) !bkpt(3),bkpt_b(3),qss(3)
      56             :       real,    intent (in) :: sk2(:),phi2(:) !sk2(n2d),phi2(n2d)
      57             :       integer, intent (in) :: ig(-k1d:,-k2d:,-k3d:) !ig(-k1d:k1d,-k2d:k2d,-k3d:k3d)
      58             :       integer, intent (in) :: ig2(:),nv(:),nv_b(:) !ig2(n3d),nv(jspd),nv_b(jspd)
      59             :       integer, intent (in) :: gb(3),gb2(3)
      60             :       real,    intent (in) :: vz(nmzd,2,4),bbmat(3,3)
      61             :       integer, intent (in) :: k1(:,:),k2(:,:),k3(:,:) !k1(nvd,jspd),k2(nvd,jspd),k3(nvd,jspd)
      62             :       integer, intent (in) :: k1_b(:,:),k2_b(:,:),k3_b(:,:) !k1_b(nvd,jspd),k2_b(nvd,jspd),k3_b(nvd,jspd)
      63             :       complex, intent (inout) :: uHu(:,:) !uHu(nbnd,nbnd)
      64             : 
      65             :       complex, intent (in):: vxy(nmzxyd,odi%n2d-1,2)
      66             : 
      67             : 
      68             : c     ..basis wavefunctions in the vacuum
      69           0 :       real,    allocatable :: udz(:,:)
      70           0 :       real,    allocatable :: uz(:,:)
      71           0 :       real,    allocatable :: dudz(:,:)
      72           0 :       real,    allocatable :: duz(:,:)
      73           0 :       real,    allocatable :: u(:,:,:)
      74           0 :       real,    allocatable :: ud(:,:,:)
      75           0 :       real,    allocatable :: ddnv(:,:)
      76             : 
      77           0 :       real,    allocatable :: udz_b(:,:)
      78           0 :       real,    allocatable :: uz_b(:,:)
      79           0 :       real,    allocatable :: dudz_b(:,:)
      80           0 :       real,    allocatable :: duz_b(:,:)
      81           0 :       real,    allocatable :: u_b(:,:,:)
      82           0 :       real,    allocatable :: ud_b(:,:,:)
      83           0 :       real,    allocatable :: ddnv_b(:,:)
      84             : 
      85             : c     ..local scalars..
      86             :       logical tail
      87             :       real wronk,wronk1,arg,zks,tpi,zz,rk,rk_b
      88             :       real zks0,zks0_b,arg_b,fac1,xv,yv,s(3)
      89             :       integer i,m,l,j,k,irec3,irec2,n,nv2,mp,ispin,np0,ind1,ind3
      90             :       integer nv2_b,np1,lp,addnoco,addnoco2,i3,j3
      91             :       integer sign,ivac,m1,m2,jspin2,jspin2_b,jspin2H,sign2
      92             :       complex tuu,tud,tdu,tdd
      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(:,:),gbess_b(:,:),besss(:)
     100           0 :       real, allocatable :: zmsh(:),xx(:),xximag(:),v1(:)
     101           0 :       real, allocatable :: gdu(:),gdud(:),gdu_b(:),gdud_b(:)
     102             :       REAL qssbti(3,2)
     103             : 
     104           0 :       addnoco=0
     105           0 :       addnoco2=0
     106           0 :       if(l_noco.and.(jspin.eq.2))then
     107           0 :          addnoco  = nv(1)   + nlotot
     108             :       endif
     109           0 :       if(l_noco.and.(jspin_b.eq.2))then
     110           0 :          addnoco2 = nv_b(1) + nlotot
     111             :       endif
     112             : 
     113             :       allocate ( udz(nv2d,-odi%mb:odi%mb),uz(nv2d,-odi%mb:odi%mb),
     114             :      +           dudz(nv2d,-odi%mb:odi%mb),
     115             :      +           duz(nv2d,-odi%mb:odi%mb),u(nmzd,nv2d,-odi%mb:odi%mb),
     116             :      +           ud(nmzd,nv2d,-odi%mb:odi%mb),ddnv(nv2d,-odi%mb:odi%mb),
     117             :      +           udz_b(nv2d,-odi%mb:odi%mb),uz_b(nv2d,-odi%mb:odi%mb),
     118             :      +           dudz_b(nv2d,-odi%mb:odi%mb),
     119             :      +           duz_b(nv2d,-odi%mb:odi%mb),
     120             :      +           u_b(nmzd,nv2d,-odi%mb:odi%mb),
     121             :      +           ud_b(nmzd,nv2d,-odi%mb:odi%mb),
     122             :      +           ddnv_b(nv2d,-odi%mb:odi%mb),
     123             :      +           bess(-odi%mb:odi%mb),dbss(-odi%mb:odi%mb),
     124             :      +           besss(-odi%M:odi%M),gbess(-odi%M:odi%M,nmzd),
     125             :      +           gbess_b(-odi%M:odi%M,nmzd),
     126             :      +           acof(nv2d,-odi%mb:odi%mb,nslibd),
     127             :      +           bcof(nv2d,-odi%mb:odi%mb,nslibd),
     128             :      +           acof_b(nv2d,-odi%mb:odi%mb,nslibd_b),
     129             :      +           bcof_b(nv2d,-odi%mb:odi%mb,nslibd_b),
     130             :      +           kvac3(nv2d),map1(nvd), 
     131             :      +           kvac3_b(nv2d),map1_b(nvd),
     132             :      +           zmsh(nmz),xx(nmz),xximag(nmz),v1(nmzd),
     133           0 :      +           gdu(nmzd),gdud(nmzd),gdu_b(nmzd),gdud_b(nmzd) )
     134             : 
     135           0 :       tpi = 2 * pimach() ; ic = cmplx(0.,1.)
     136             : 
     137           0 :       tail = .true.
     138           0 :       np0 = nmzxy + 1
     139           0 :       np1 = nmz + 1
     140           0 :       ivac = 1
     141             : 
     142           0 :       jspin2 = jspin
     143           0 :       jspin2_b = jspin_b
     144           0 :       jspin2H = ico
     145           0 :       if(l_soc.and.jspins.eq.1) then
     146             :          jspin2 = 1
     147             :          jspin2_b = 1
     148             :          jspin2H = 1
     149             :       endif
     150             : 
     151           0 :       sign = 1.0
     152           0 :       IF (ico.EQ.4) sign=-1.0  ! for c.c. of vz
     153             : 
     154             : 
     155           0 :       acof(:,:,:) = cmplx(0.,0.) ; bcof(:,:,:) = cmplx(0.,0.)
     156           0 :       acof_b(:,:,:) = cmplx(0.,0.) ; bcof_b(:,:,:) = cmplx(0.,0.)
     157             : 
     158           0 :       nv2 = 0 ; nv2_b = 0
     159             : 
     160           0 :       do 20 k = 1,nv(jspin)
     161           0 :          do 10 j = 1,nv2
     162           0 :             if (k3(k,jspin).eq.kvac3(j)) then
     163           0 :                map1(k) = j
     164           0 :                goto 20
     165             :             endif
     166           0 :  10      continue
     167           0 :          nv2 = nv2 + 1
     168           0 :          if (nv2.gt.nv2d) stop 'nv2d'
     169           0 :          kvac3(nv2) = k3(k,jspin)
     170           0 :          map1(k) = nv2
     171           0 :  20   continue
     172             : 
     173           0 :       do 21 k = 1,nv_b(jspin_b)
     174           0 :          do 11 j = 1,nv2_b
     175           0 :             if (k3_b(k,jspin_b).eq.kvac3_b(j)) then
     176           0 :                map1_b(k) = j
     177           0 :                goto 21
     178             :             endif
     179           0 :  11      continue
     180           0 :          nv2_b = nv2_b + 1
     181           0 :          if (nv2_b.gt.nv2d) stop 'nv2d'
     182           0 :          kvac3_b(nv2_b) = k3_b(k,jspin_b)
     183           0 :          map1_b(k) = nv2_b
     184           0 :  21   continue
     185             : 
     186           0 :       wronk = 2.0
     187             : 
     188             : c...for the k-point
     189             : 
     190           0 :       qssbti(1,1) = - qss(1)/2.     
     191           0 :       qssbti(2,1) = - qss(2)/2.     
     192           0 :       qssbti(1,2) = + qss(1)/2.     
     193           0 :       qssbti(2,2) = + qss(2)/2.
     194           0 :       qssbti(3,1) = - qss(3)/2.
     195           0 :       qssbti(3,2) = + qss(3)/2.
     196           0 :       DO ispin = 1,1 ! jspins
     197             :       CALL od_abvac(
     198             :      >      cell,vacuum,DIMENSION,stars,oneD,
     199             :      >      qssbti(3,jspin),odi%n2d,
     200             :      >      wronk,evac,bkpt,odi%M,odi%mb,
     201             :      >      vz(1,ivac,jspin2),kvac3,nv2,
     202             :      <      uz(1,-odi%mb),duz(1,-odi%mb),u(1,1,-odi%mb),udz(1,-odi%mb),
     203           0 :      <      dudz(1,-odi%mb),ddnv(1,-odi%mb),ud(1,1,-odi%mb))
     204             :       ENDDO
     205             : 
     206           0 :       do k = 1,nv(jspin)
     207           0 :          l = map1(k)
     208           0 :          irec3 = ig(k1(k,jspin),k2(k,jspin),k3(k,jspin))
     209           0 :          if (irec3.ne.0) then
     210           0 :             irec2 = ig2(irec3)
     211           0 :             zks = sk2(irec2)*z1
     212           0 :             arg = phi2(irec2)
     213           0 :             call cylbes(odi%mb,zks,bess)
     214           0 :             call dcylbs(odi%mb,zks,bess,dbss)
     215           0 :             do m = -odi%mb,odi%mb
     216             :                wronk1 = uz(l,m)*dudz(l,m) -
     217           0 :      -              udz(l,m)*duz(l,m)
     218             :                avac = exp(-cmplx(0.0,m*arg))*(ic**m)*
     219             :      *              cmplx(dudz(l,m)*bess(m) -
     220             :      +              udz(l,m)*sk2(irec2)*dbss(m),0.0)/
     221           0 :      /              ((wronk1)*sqrt(omtil))
     222             :                bvac = exp(-cmplx(0.0,m*arg))*(ic**m)*
     223             :      *              cmplx(-duz(l,m)*bess(m) +
     224             :      -              uz(l,m)*sk2(irec2)*dbss(m),0.0)/
     225           0 :      /              ((wronk1)*sqrt(omtil))
     226           0 :                IF(zMat%l_real) THEN
     227           0 :                   do n = 1,nslibd
     228             :                       acof(l,m,n) = acof(l,m,n) +
     229           0 :      +                   zMat%z_r(k+addnoco,n)*avac
     230             : c     +                    conjg(zMat%z_r(k,n))*avac
     231             :                       bcof(l,m,n) = bcof(l,m,n) +
     232           0 :      +                   zMat%z_r(k+addnoco,n)*bvac
     233             : c     +                    conjg(zMat%z_r(k,n))*bvac
     234             :                   enddo
     235             :                ELSE
     236           0 :                   do n = 1,nslibd
     237             :                       acof(l,m,n) = acof(l,m,n) +
     238           0 :      +                   zMat%z_c(k+addnoco,n)*avac
     239             : c     +                    conjg(zMat%z_c(k,n))*avac
     240             :                       bcof(l,m,n) = bcof(l,m,n) +
     241           0 :      +                   zMat%z_c(k+addnoco,n)*bvac
     242             : c     +                    conjg(zMat%z_c(k,n))*bvac
     243             :                   enddo
     244             :                END IF
     245             :             enddo      ! -mb:mb
     246             :          endif
     247             :       enddo
     248             : 
     249             : c...for the b-point
     250             : 
     251           0 :       DO ispin = 1,1 ! jspins
     252             :       call od_abvac(
     253             :      >      cell,vacuum,DIMENSION,stars,oneD,
     254             :      >      qssbti(3,jspin_b),odi%n2d,
     255             :      >      wronk,evac_b,bkpt_b,odi%M,odi%mb,
     256             :      >      vz(1,ivac,jspin2_b),kvac3_b,nv2_b,
     257             :      <      uz_b(1,-odi%mb),duz_b(1,-odi%mb),u_b(1,1,-odi%mb),
     258             :      <      udz_b(1,-odi%mb),
     259           0 :      <      dudz_b(1,-odi%mb),ddnv_b(1,-odi%mb),ud_b(1,1,-odi%mb))
     260             :       ENDDO
     261             : 
     262           0 :       do k = 1,nv_b(jspin_b)
     263           0 :          l = map1_b(k)
     264           0 :          irec3 = ig(k1_b(k,jspin_b),k2_b(k,jspin_b),k3_b(k,jspin_b))
     265           0 :          if (irec3.ne.0) then
     266           0 :             irec2 = ig2(irec3)
     267           0 :             zks = sk2(irec2)*z1
     268           0 :             arg = phi2(irec2)
     269           0 :             call cylbes(odi%mb,zks,bess)
     270           0 :             call dcylbs(odi%mb,zks,bess,dbss)
     271           0 :             do m = -odi%mb,odi%mb
     272             :                wronk1 = uz_b(l,m)*dudz_b(l,m) -
     273           0 :      -              udz_b(l,m)*duz_b(l,m)
     274             :                avac = exp(-cmplx(0.0,m*arg))*(ic**m)*
     275             :      *              cmplx(dudz_b(l,m)*bess(m) -
     276             :      +              udz_b(l,m)*sk2(irec2)*dbss(m),0.0)/
     277           0 :      /              ((wronk1)*sqrt(omtil))
     278             :                bvac = exp(-cmplx(0.0,m*arg))*(ic**m)*
     279             :      &              cmplx(-duz_b(l,m)*bess(m) +
     280             :      -              uz_b(l,m)*sk2(irec2)*dbss(m),0.0)/
     281           0 :      /              ((wronk1)*sqrt(omtil))
     282             : 
     283           0 :                IF(zMat_b%l_real) THEN
     284           0 :                   do n = 1,nslibd_b
     285             :                       acof_b(l,m,n) = acof_b(l,m,n) +
     286           0 :      +                   zMat_b%z_r(k+addnoco2,n)*avac
     287             : c     +                    conjg(zMat_b%z_r(k,n))*avac
     288             :                       bcof_b(l,m,n) = bcof_b(l,m,n) +
     289           0 :      +                   zMat_b%z_r(k+addnoco2,n)*bvac
     290             : c     +                    conjg(zMat_b%z_r(k,n))*bvac
     291             :                   enddo
     292             :                ELSE
     293           0 :                   do n = 1,nslibd_b
     294             :                       acof_b(l,m,n) = acof_b(l,m,n) +
     295           0 :      +                   zMat_b%z_c(k+addnoco2,n)*avac
     296             : c     +                    conjg(zMat_b%z_c(k,n))*avac
     297             :                       bcof_b(l,m,n) = bcof_b(l,m,n) +
     298           0 :      +                   zMat_b%z_c(k+addnoco2,n)*bvac
     299             : c     +                    conjg(zMat_b%z_c(k,n))*bvac
     300             :                   enddo
     301             :                END IF
     302             :             enddo      ! -mb:mb
     303             :          endif 
     304             :       enddo          ! k = 1,nv
     305             : 
     306             : c  now actually computing the uHu matrix
     307             : 
     308           0 :       irec3 = ig(gb(1),gb(2),gb(3))
     309           0 :       if (irec3.eq.0) stop 'Gb is not in the list of Gs'
     310           0 :       irec2 = ig2(irec3)
     311           0 :       zks0 = sk2(irec2)
     312           0 :       arg = phi2(irec2)
     313             : 
     314           0 :       irec3 = ig(gb2(1),gb2(2),gb2(3))
     315           0 :       if (irec3.eq.0) stop 'Gb2 is not in the list of Gs'
     316           0 :       irec2 = ig2(irec3)
     317           0 :       zks0_b = sk2(irec2)
     318           0 :       arg_b = phi2(irec2)
     319             : 
     320           0 :       gbess(:,:) = 0.
     321           0 :       gbess_b(:,:) = 0.
     322           0 :       do i = 1,nmz
     323             :          ! set up mesh
     324           0 :          zmsh(i) = z1 + (i-1)*delz
     325           0 :          zz = sqrt(zmsh(i))
     326             : 
     327             :          ! transformation u --> v = sqrt(z)*u
     328             :          ! thus we can use simplified 'pseudopotential'
     329             :          ! and skip r in integration dx dy = r dr dphi
     330           0 :          u(i,:,:) = zz*u(i,:,:)              
     331           0 :          ud(i,:,:) = zz*ud(i,:,:)
     332           0 :          u_b(i,:,:) = zz*u_b(i,:,:)
     333           0 :          ud_b(i,:,:) = zz*ud_b(i,:,:)
     334             : 
     335             :          ! cyl. Bessel at G(k+b1)
     336           0 :          zks = zks0*zmsh(i)
     337           0 :          besss(:) = 0.
     338           0 :          call cylbes(odi%M,zks,besss)
     339           0 :          do m = -odi%M,odi%M
     340           0 :             gbess(m,i) = besss(m)
     341             :          enddo
     342             : 
     343             :          ! cyl. Bessel at G(k+b2)
     344           0 :          zks = zks0_b*zmsh(i)
     345           0 :          besss(:) = 0.
     346           0 :          call cylbes(odi%M,zks,besss)
     347           0 :          do m = -odi%M,odi%M
     348           0 :             gbess_b(m,i) = besss(m)
     349             :          enddo
     350             :       enddo
     351             : 
     352             :      
     353             :       ! calculate uHu matrix elements
     354           0 :       do l = 1,nv2
     355             : 
     356           0 :          j3 = kvac3(l) - gb(3)
     357             : 
     358           0 :       do lp = 1,nv2_b
     359             : 
     360           0 :          i3 = j3 - kvac3_b(lp) + gb2(3)
     361             : 
     362           0 :          ind3 = ig(0,0,sign2*i3)
     363           0 :          IF (ind3.EQ.0) CYCLE
     364             : 
     365           0 :          DO m = -odi%mb, odi%mb
     366           0 :           DO mp = -odi%mb, odi%mb
     367           0 :            IF ((mp.EQ.m) .OR. ((iabs(m ).LE.odi%m_cyl) .AND.
     368           0 :      >                         (iabs(mp).LE.odi%m_cyl))) THEN
     369             : 
     370           0 :             DO m1 = -odi%mb, odi%mb
     371           0 :              DO m2 = -odi%mb, odi%mb
     372             : 
     373           0 :               ind1 = odi%ig(sign2*i3,m-m1-mp+m2)  ! TODO: sign also for m-part?
     374           0 :               if(sign2.ne.1) stop 'sign2, check m-part, wann_uHu_od_vac'
     375             : 
     376           0 :               IF(ind1.EQ.0) CYCLE
     377             : 
     378           0 :               ind1 = ind1 - 1
     379           0 :               IF(ind1.NE.0) THEN ! warping components Gz=/=0, m=/=0
     380             :                                  ! need integral with vxy and cyl. Bessel
     381             :                  ! tuu
     382           0 :                  DO i=1,nmzxy
     383             :                     fac1 = u  (i,l ,m )*gbess  (m1,i)
     384           0 :      >                    *u_b(i,lp,mp)*gbess_b(m2,i)
     385           0 :                     xx    (np0-i) = fac1*real (vxy(i,ind1,ivac))
     386           0 :                     xximag(np0-i) = fac1*aimag(vxy(i,ind1,ivac))
     387             :                  ENDDO
     388           0 :                  call intgz0(xx,delz,nmzxy,xv,tail)
     389           0 :                  call intgz0(xximag,delz,nmzxy,yv,tail)
     390           0 :                  tuu = cmplx(xv,-yv)                 
     391             : 
     392             :                  ! tud
     393           0 :                  DO i=1,nmzxy
     394             :                     fac1 = u   (i,l ,m )*gbess  (m1,i)
     395           0 :      >                    *ud_b(i,lp,mp)*gbess_b(m2,i)
     396           0 :                     xx    (np0-i) = fac1*real (vxy(i,ind1,ivac))
     397           0 :                     xximag(np0-i) = fac1*aimag(vxy(i,ind1,ivac))
     398             :                  ENDDO
     399           0 :                  call intgz0(xx,delz,nmzxy,xv,tail)
     400           0 :                  call intgz0(xximag,delz,nmzxy,yv,tail)
     401           0 :                  tud = cmplx(xv,-yv)                 
     402             : 
     403             :                  ! tdu
     404           0 :                  DO i=1,nmzxy
     405             :                     fac1 = ud (i,l ,m )*gbess  (m1,i)
     406           0 :      >                    *u_b(i,lp,mp)*gbess_b(m2,i)
     407           0 :                     xx    (np0-i) = fac1*real (vxy(i,ind1,ivac))
     408           0 :                     xximag(np0-i) = fac1*aimag(vxy(i,ind1,ivac))
     409             :                  ENDDO
     410           0 :                  call intgz0(xx,delz,nmzxy,xv,tail)
     411           0 :                  call intgz0(xximag,delz,nmzxy,yv,tail)
     412           0 :                  tdu = cmplx(xv,-yv)                 
     413             : 
     414             :                  ! tdd
     415           0 :                  DO i=1,nmzxy
     416             :                     fac1 = ud  (i,l ,m )*gbess  (m1,i)
     417           0 :      >                    *ud_b(i,lp,mp)*gbess_b(m2,i)
     418           0 :                     xx    (np0-i) = fac1*real (vxy(i,ind1,ivac))
     419           0 :                     xximag(np0-i) = fac1*aimag(vxy(i,ind1,ivac))
     420             :                  ENDDO
     421           0 :                  call intgz0(xx,delz,nmzxy,xv,tail)
     422           0 :                  call intgz0(xximag,delz,nmzxy,yv,tail)
     423           0 :                  tdd = cmplx(xv,-yv)                
     424             :  
     425             :               ELSE ! non-warping components Gz==0, m==0
     426             :                    ! need integral with XXXXXX
     427             : 
     428           0 :                  IF ((ico.EQ.1) .OR. (ico.EQ.2)) THEN ! spin-diagonal
     429             : 
     430             :                     ! determine second derivative of (u*gbess) etc.
     431             :                     CALL d2fdz2(nmzd,nmz,zmsh,delz,
     432           0 :      >                          u(:,l,m),gbess(m1,:),gdu)
     433             :                     CALL d2fdz2(nmzd,nmz,zmsh,delz,
     434           0 :      >                          ud(:,l,m),gbess(m1,:),gdud)
     435             :                     CALL d2fdz2(nmzd,nmz,zmsh,delz,
     436           0 :      >                          u_b(:,lp,mp),gbess_b(m2,:),gdu_b)
     437             :                     CALL d2fdz2(nmzd,nmz,zmsh,delz,
     438           0 :      >                          ud_b(:,lp,mp),gbess_b(m2,:),gdud_b)
     439             : 
     440             :                     ! determine |G+k+b1|^2 and |G'+k+b2|^2
     441           0 :                     s(1) = 0.0
     442           0 :                     s(2) = 0.0
     443           0 :                     s(3) = bkpt(3) + kvac3(l) + qssbti(3,jspin)!-gb(3)
     444           0 :                     rk = dot_product(s,matmul(bbmat,s))
     445             : !                    rk   = dotirp(s,s,bbmat)
     446             : 
     447             :                     s(1) = 0.0
     448             :                     s(2) = 0.0
     449           0 :                     s(3) = bkpt_b(3) + kvac3_b(lp) + qssbti(3,jspin_b)!-gb2(3)
     450           0 :                     rk_b = dot_product(s,matmul(bbmat,s))
     451             : !                    rk_b = dotirp(s,s,bbmat)
     452             : 
     453             :                     ! construct symmetrized 'pseudopotential'
     454             :                     ! TODO: valid to simply symmetrize?
     455           0 :                     DO i=1,nmzd
     456             :                        v1(i) = vz(i,ivac,jspin2H) 
     457             :      >                       + (m*m+mp*mp)/(4.*zmsh(i)*zmsh(i))
     458           0 :      >                       - 1./(8.*zmsh(i)*zmsh(i))
     459             :                     ENDDO
     460             : 
     461             :                     ! tuu
     462           0 :                     DO i=1,nmz
     463             :                        fac1 = u  (i,l ,m )*gbess  (m1,i)
     464           0 :      >                       *u_b(i,lp,mp)*gbess_b(m2,i)
     465             :                        xx(np1-i) = fac1*( 0.25*(rk+rk_b) + v1(i) )
     466             :      >                   -0.25*( gdu  (i)*u_b(i,lp,mp)*gbess_b(m2,i)
     467           0 :      >                          +gdu_b(i)*u  (i,l ,m )*gbess  (m1,i) )
     468             :                     ENDDO
     469           0 :                     call intgz0(xx,delz,nmz,xv,tail)
     470           0 :                     tuu = cmplx(xv,0.0)
     471             : 
     472             :                     ! tud
     473           0 :                     DO i=1,nmz
     474             :                        fac1 = u   (i,l ,m )*gbess  (m1,i)
     475           0 :      >                       *ud_b(i,lp,mp)*gbess_b(m2,i)
     476             :                        xx(np1-i) = fac1*( 0.25*(rk+rk_b) + v1(i) )
     477             :      >                   -0.25*( gdu   (i)*ud_b(i,lp,mp)*gbess_b(m2,i)
     478           0 :      >                          +gdud_b(i)*u   (i,l ,m )*gbess  (m1,i) )
     479             :                     ENDDO
     480           0 :                     call intgz0(xx,delz,nmz,xv,tail)
     481           0 :                     tud = cmplx(xv,0.0)
     482             : 
     483             :                     ! tdu
     484           0 :                     DO i=1,nmz
     485             :                        fac1 = ud (i,l ,m )*gbess  (m1,i)
     486           0 :      >                       *u_b(i,lp,mp)*gbess_b(m2,i)
     487             :                        xx(np1-i) = fac1*( 0.25*(rk+rk_b) + v1(i) )
     488             :      >                   -0.25*( gdud (i)*u_b(i,lp,mp)*gbess_b(m2,i)
     489           0 :      >                          +gdu_b(i)*ud (i,l ,m )*gbess  (m1,i) )
     490             :                     ENDDO
     491           0 :                     call intgz0(xx,delz,nmz,xv,tail)
     492           0 :                     tdu = cmplx(xv,0.0)
     493             : 
     494             :                     ! tdd
     495           0 :                     DO i=1,nmz
     496             :                        fac1 = ud  (i,l ,m )*gbess  (m1,i)
     497           0 :      >                       *ud_b(i,lp,mp)*gbess_b(m2,i)
     498             :                        xx(np1-i) = fac1*( 0.25*(rk+rk_b) + v1(i) )
     499             :      >                   -0.25*( gdud  (i)*ud_b(i,lp,mp)*gbess_b(m2,i)
     500           0 :      >                          +gdud_b(i)*ud  (i,l ,m )*gbess  (m1,i) )
     501             :                     ENDDO
     502           0 :                     call intgz0(xx,delz,nmz,xv,tail)
     503           0 :                     tdd = cmplx(xv,0.0)
     504             : 
     505             :                  ELSE ! spin-off-diagonal
     506             : 
     507             :                     ! tuu
     508           0 :                     DO i=1,nmz
     509             :                        fac1 = u  (i,l ,m )*gbess  (m1,i)
     510           0 :      >                       *u_b(i,lp,mp)*gbess_b(m2,i)
     511           0 :                        xx    (np1-i) = fac1*vz(i,ivac,3)
     512           0 :                        xximag(np1-i) = fac1*vz(i,ivac,4)*sign
     513             :                     ENDDO
     514           0 :                     call intgz0(xx,delz,nmz,xv,tail)
     515           0 :                     call intgz0(xximag,delz,nmz,yv,tail)
     516           0 :                     tuu = cmplx(xv,-yv)
     517             : 
     518             :                     ! tud
     519           0 :                     DO i=1,nmz
     520             :                        fac1 = u   (i,l ,m )*gbess  (m1,i)
     521           0 :      >                       *ud_b(i,lp,mp)*gbess_b(m2,i)
     522           0 :                        xx    (np1-i) = fac1*vz(i,ivac,3)
     523           0 :                        xximag(np1-i) = fac1*vz(i,ivac,4)*sign
     524             :                     ENDDO
     525           0 :                     call intgz0(xx,delz,nmz,xv,tail)
     526           0 :                     call intgz0(xximag,delz,nmz,yv,tail)
     527           0 :                     tud = cmplx(xv,-yv)
     528             : 
     529             :                     ! tdu
     530           0 :                     DO i=1,nmz
     531             :                        fac1 = ud (i,l ,m )*gbess  (m1,i)
     532           0 :      >                       *u_b(i,lp,mp)*gbess_b(m2,i)
     533           0 :                        xx    (np1-i) = fac1*vz(i,ivac,3)
     534           0 :                        xximag(np1-i) = fac1*vz(i,ivac,4)*sign
     535             :                     ENDDO
     536           0 :                     call intgz0(xx,delz,nmz,xv,tail)
     537           0 :                     call intgz0(xximag,delz,nmz,yv,tail)
     538           0 :                     tdu = cmplx(xv,-yv)
     539             : 
     540             :                     ! tdd
     541           0 :                     DO i=1,nmz
     542             :                        fac1 = ud  (i,l ,m )*gbess  (m1,i)
     543           0 :      >                       *ud_b(i,lp,mp)*gbess_b(m2,i)
     544           0 :                        xx    (np1-i) = fac1*vz(i,ivac,3)
     545           0 :                        xximag(np1-i) = fac1*vz(i,ivac,4)*sign
     546             :                     ENDDO
     547           0 :                     call intgz0(xx,delz,nmz,xv,tail)
     548           0 :                     call intgz0(xximag,delz,nmz,yv,tail)
     549           0 :                     tdd = cmplx(xv,-yv)
     550             : 
     551             :                  ENDIF ! ((ico.EQ.1) .OR. (ico.EQ.2))
     552             :             
     553             :               ENDIF ! (ind1.NE.0)
     554             :         
     555             :               ! determine phase factor
     556           0 :               phasfc = chi*exp(-cmplx(0.0,m2*arg_b-m1*arg)) 
     557             : 
     558             :               ! contraction of integrals with a,b coefficients
     559             :               ! yields contribution to uHu matrix
     560           0 :               do i = 1,nslibd
     561           0 :                 do j = 1,nslibd_b
     562             :                    uHu(i,j) = uHu(i,j) + phasfc*area*(
     563             :      *             acof(l,m,i)*conjg(acof_b(lp,mp,j))*tuu +
     564             :      +             acof(l,m,i)*conjg(bcof_b(lp,mp,j))*tud +
     565             :      +             bcof(l,m,i)*conjg(acof_b(lp,mp,j))*tdu +
     566           0 :      +             bcof(l,m,i)*conjg(bcof_b(lp,mp,j))*tdd ) 
     567             :                 enddo 
     568             :               enddo
     569             : 
     570             :              ENDDO ! m2
     571             :             ENDDO  ! m1
     572             :            ENDIF ! noncyl. contributions
     573             :           ENDDO ! mp
     574             :          ENDDO  ! m
     575             : 
     576             :       enddo  ! lp
     577             :       enddo  ! l
     578             : 
     579           0 :       deallocate ( udz,uz,dudz,duz,u,ud,ddnv,bess,dbss,acof,bcof )
     580           0 :       deallocate ( udz_b,uz_b,dudz_b,duz_b,u_b,ud_b,ddnv_b,gbess,besss )
     581           0 :       deallocate ( acof_b,bcof_b,gbess_b )
     582           0 :       deallocate ( zmsh,xx,xximag,v1 )
     583           0 :       deallocate ( gdu, gdud, gdu_b, gdud_b )
     584             : 
     585           0 :       END SUBROUTINE wann_uHu_od_vac
     586             :       END MODULE m_wann_uHu_od_vac

Generated by: LCOV version 1.13