LCOV - code coverage report
Current view: top level - wannier/uhu - wann_uHu_vac.F (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 290 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 FLAPW film 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_vac
      10             :       USE m_juDFT
      11             : 
      12             :       CONTAINS
      13           0 :       SUBROUTINE wann_uHu_vac(
      14             :      >     chi,l_noco,l_soc,zrfs,jspins,nlotot,qss,nbnd,z1,nmzxyd,
      15           0 :      >     nmzd,n2d,nv2d,k1d,k2d,k3d,n3d,nvac,ig,rgphs,nmzxy,
      16           0 :      >     nmz,delz,ig2,nq2,kv2,area,bmat,bbmat,evac,evac_b,
      17           0 :      >     bkpt,bkpt_b,vzxy,vz,nslibd,nslibd_b,jspin,jspin_b,
      18           0 :      >     ico,k1,k2,k3,k1_b,k2_b,k3_b,jspd,nvd,nbasfcn,neigd,
      19           0 :      >     zMat,zMat_b,nv,nv_b,omtil,gb,gb2,sign2,uHu)
      20             : #include "cpp_double.h"
      21             : 
      22             :       USE m_types
      23             :       use m_constants, only : pimach
      24             :       use m_intgr, only : intgz0
      25             :       use m_d2fdz2cmplx
      26             :       USE m_vacuz
      27             :       USE m_vacudz
      28             : 
      29             :       implicit none
      30             : 
      31             :       TYPE(t_zmat), INTENT(IN) :: zMat, zMat_b
      32             : 
      33             : c     .. scalar Arguments..
      34             :       logical, intent (in) :: l_noco,l_soc,zrfs
      35             :       integer, intent (in) :: nlotot,jspin_b,n2d,jspins,ico
      36             :       integer, intent (in) :: nmzxyd,nmzd,nv2d,k1d,k2d,k3d,n3d,nbnd
      37             :       integer, intent (in) :: nmzxy,nmz,nslibd,nslibd_b,nvac
      38             :       integer, intent (in) :: jspin,jspd,nvd,nq2
      39             :       integer, intent (in) :: nbasfcn,neigd
      40             :       real,    intent (in) :: delz,z1,omtil,area
      41             :       complex, intent (in) :: chi
      42             : 
      43             : c     ..array arguments..
      44             :       real,    intent (in) :: bkpt(3),bkpt_b(3),evac(2),evac_b(2)
      45             :       real,    intent (in) :: qss(3)
      46             :       real,    intent (in) :: vz(nmzd,2,4)
      47             :       real,    intent (in) :: bbmat(3,3),bmat(3,3)
      48             :       complex, intent (in) :: rgphs(-k1d:k1d,-k2d:k2d,-k3d:k3d)
      49             :       integer, intent (in) :: ig(-k1d:k1d,-k2d:k2d,-k3d:k3d)
      50             :       integer, intent (in) :: ig2(n3d),nv(jspd),nv_b(jspd)
      51             :       integer, intent (in) :: kv2(2,n2d)
      52             :       integer, intent (in) :: gb(3),gb2(3)
      53             :       integer, intent (in) :: k1(nvd,jspd),k2(nvd,jspd),k3(nvd,jspd)
      54             :       integer, intent (in) :: k1_b(nvd,jspd),k2_b(nvd,jspd)
      55             :       integer, intent (in) :: k3_b(nvd,jspd)
      56             :       complex, intent (in) :: vzxy(nmzxyd,n2d-1,2)
      57             :       complex, intent (inout) :: uHu(nbnd,nbnd)
      58             : 
      59             : c     ..local arrays..
      60           0 :       complex, allocatable :: ac(:,:),bc(:,:),ac_b(:,:),bc_b(:,:)
      61           0 :       complex, allocatable :: vv(:),fac(:),fac_b(:),fac1(:)
      62           0 :       complex, allocatable :: du(:),due(:),du_b(:),due_b(:)
      63           0 :       complex, allocatable :: vxy(:,:),vxy_help(:)
      64           0 :       real,    allocatable :: dt(:),dte(:)
      65           0 :       real,    allocatable :: t(:),te(:),tei(:)
      66           0 :       real,    allocatable :: u(:,:),ue(:,:),v(:)
      67           0 :       real,    allocatable :: dt_b(:),dte_b(:)
      68           0 :       real,    allocatable :: t_b(:),te_b(:),tei_b(:)
      69           0 :       real,    allocatable :: u_b(:,:),ue_b(:,:)
      70           0 :       real,    allocatable :: zmsh(:)
      71             : 
      72             : c     ..local scalars..
      73             :       logical tail
      74             :       real wronk,arg,zks,zks_b,tpi,vz0(2),vz0_b(2),scale,evacp,ev,const
      75           0 :       real xx(nmz),xximag(nmz),phase,phase2,rk,rk_b
      76             :       real :: qss1,qss2,xv,yv,s(3)
      77             :       integer i,m,l,j,k,n,nv2,nv2_b,ivac,n2,n2_b,sign,ik,sign2
      78             :       integer :: lp,np0,np1,addnoco,addnoco2,jspin2,jspin2_b,jspin2H
      79             :       complex :: av,bv,ic,c_1
      80             : c      complex :: tuu,tud,tdu,tdd
      81           0 :       complex, allocatable :: tuu(:,:),tud(:,:),tdu(:,:),tdd(:,:)
      82           0 :       complex, allocatable :: mat(:,:)
      83           0 :       integer, allocatable :: kvac1(:),kvac2(:),map2(:)
      84           0 :       integer, allocatable :: kvac1_b(:),kvac2_b(:),map2_b(:)
      85             :       integer symvac,symvacvac,igvm2,igvm2i
      86             :       integer :: j1,j2,i1,i2,ind3,ind2
      87             :       character(len=20) :: fname
      88             : c     ..intrinsic functions..
      89             :       intrinsic aimag,cmplx,conjg,real,sqrt
      90             : 
      91             :       allocate ( ac(nv2d,nslibd),bc(nv2d,nslibd),
      92             :      +           ac_b(nv2d,nslibd_b),bc_b(nv2d,nslibd_b),
      93             :      +           dt(nv2d),dte(nv2d),t(nv2d),te(nv2d),
      94             :      +           tei(nv2d),u(nmzd,nv2d),ue(nmzd,nv2d),
      95             :      +           dt_b(nv2d),dte_b(nv2d),t_b(nv2d),te_b(nv2d),
      96             :      +           tei_b(nv2d),u_b(nmzd,nv2d),ue_b(nmzd,nv2d),
      97             :      +           v(3),kvac1(nv2d),kvac2(nv2d),map2(nvd),
      98             :      +           kvac1_b(nv2d),kvac2_b(nv2d),map2_b(nvd),
      99             :      +           vv(nmzd),zmsh(nmzd),vxy(nmzxyd,n2d-1),vxy_help(n2d-1),
     100             :      +           fac1(nmzd), fac(nmzd), fac_b(nmzd),
     101             :      +           du(nmzd), due(nmzd), du_b(nmzd), due_b(nmzd),
     102             :      +           tuu(nv2d,nv2d),tud(nv2d,nv2d),tdu(nv2d,nv2d),
     103           0 :      +           tdd(nv2d,nv2d),mat(nslibd_b,nv2d) )
     104             : 
     105           0 :       tpi = 2 * pimach() ; ic = cmplx(0.,1.)
     106           0 :       tail = .true.
     107           0 :       np0 = nmzxy + 1
     108           0 :       np1 = nmz + 1
     109             : 
     110           0 :       jspin2 = jspin
     111           0 :       jspin2_b = jspin_b
     112           0 :       jspin2H = ico
     113           0 :       if(l_soc.and.jspins.eq.1) then
     114             :          jspin2 = 1
     115             :          jspin2_b = 1
     116             :          jspin2H = 1
     117             :       endif
     118             : 
     119             : c.. determining the indexing array (in-plane stars)
     120             : c.. for the k-point
     121             : 
     122           0 :       wronk = 2.0
     123           0 :       const = 1.0 / ( sqrt(omtil)*wronk )
     124             : 
     125           0 :       do ivac = 1,2
     126           0 :          vz0(ivac) = vz(nmz,ivac,jspin2)
     127           0 :          vz0_b(ivac) = vz(nmz,ivac,jspin2_b)
     128             : c         write(*,*)'vz0',ivac,vz0(ivac),vz0_b(ivac)
     129             :       enddo
     130             : 
     131           0 :       n2 = 0 ; n2_b = 0
     132             : 
     133           0 :       addnoco=0
     134           0 :       addnoco2=0
     135           0 :       if(l_noco.and.jspin.eq.2)then
     136           0 :         addnoco=nv(1)+nlotot
     137             :       endif
     138           0 :       if(l_noco.and.jspin_b.eq.2)then
     139           0 :         addnoco2=nv_b(1)+nlotot
     140             :       endif
     141             : 
     142           0 :       do 40 k = 1,nv(jspin)
     143           0 :          do 30 j = 1,n2
     144           0 :             if ( k1(k,jspin).eq.kvac1(j) .and.
     145             :      +          k2(k,jspin).eq.kvac2(j) ) then
     146           0 :                 map2(k) = j
     147           0 :                 goto 40
     148             :              endif 
     149           0 :  30      continue
     150           0 :          n2 = n2 + 1
     151             :          
     152           0 :          IF(n2>nv2d) then
     153           0 :             write(*,*)n2,nv2d,'jspin',jspin
     154             :          endif
     155             :  
     156           0 :          IF (n2>nv2d)  CALL juDFT_error("wannier uHu vac",calledby
     157           0 :      +        ="wann_uHu_vac")
     158             : 
     159           0 :          kvac1(n2) = k1(k,jspin)
     160           0 :          kvac2(n2) = k2(k,jspin)
     161           0 :          map2(k) = n2
     162           0 :  40   continue
     163             : 
     164             : 
     165             : c.. and for the b-point
     166             :  
     167           0 :       do 41 k = 1,nv_b(jspin_b)
     168           0 :          do 31 j = 1,n2_b
     169           0 :             if ( k1_b(k,jspin_b).eq.kvac1_b(j) .and.
     170             :      +          k2_b(k,jspin_b).eq.kvac2_b(j) ) then
     171           0 :                 map2_b(k) = j
     172           0 :                 goto 41
     173             :              endif
     174           0 :  31      continue
     175           0 :          n2_b = n2_b + 1
     176             :      
     177           0 :          IF(n2_b>nv2d) then
     178           0 :             write(*,*)n2_b,nv2d,'jspin_b',jspin_b
     179             :          endif
     180             : 
     181           0 :          IF (n2_b>nv2d)  CALL juDFT_error("wannier uHu vac",calledby
     182           0 :      +        ="wann_uHu_vac")
     183             : 
     184           0 :          kvac1_b(n2_b) = k1_b(k,jspin_b)
     185           0 :          kvac2_b(n2_b) = k2_b(k,jspin_b)
     186           0 :          map2_b(k) = n2_b
     187           0 :  41   continue
     188             : 
     189             : 
     190             : 
     191             : 
     192             : c...cycle by the vacua
     193           0 :       do 140 ivac = 1,2!nvac
     194             : 
     195           0 :       sign = 3. - 2.*ivac
     196           0 :       evacp = evac(ivac)
     197             : 
     198           0 :       vxy(:,:) = vzxy(:,:,ivac)
     199           0 :       IF(l_noco) THEN ! symmetrization
     200           0 :        IF (nvac.EQ.1 .and. ivac.EQ.2 .AND.(.NOT.zrfs) ) THEN
     201           0 :         DO i=1,nmzxy
     202           0 :          DO igvm2 = 2,nq2
     203           0 :           igvm2i = ig2(ig(-kv2(1,igvm2),-kv2(2,igvm2),0))
     204           0 :           vxy_help(igvm2-1) = vxy(i,igvm2i-1)
     205             :          ENDDO
     206           0 :          DO igvm2 = 2,nq2
     207           0 :           vxy(i,igvm2-1) = vxy_help(igvm2-1)
     208             :          ENDDO
     209             :         ENDDO
     210             :        ENDIF
     211             :       ENDIF
     212             : 
     213           0 :       tuu = cmplx(0.,0.)
     214           0 :       tud = cmplx(0.,0.)
     215           0 :       tdu = cmplx(0.,0.)
     216           0 :       tdd = cmplx(0.,0.)
     217             : 
     218           0 :       nv2 = n2 ; nv2_b = n2_b
     219             : 
     220             : c.. the body of the routine
     221             : 
     222           0 :       qss1=0.0
     223           0 :       qss2=0.0
     224           0 :       if(l_noco.and.jspin.eq.1)then
     225           0 :         qss1=-qss(1)/2.0
     226           0 :         qss2=-qss(2)/2.0
     227           0 :       elseif(l_noco.and.jspin.eq.2)then
     228           0 :         qss1=qss(1)/2.0
     229           0 :         qss2=qss(2)/2.0
     230             :       endif
     231             : 
     232           0 :       do ik = 1,nv2
     233           0 :          v(1) = bkpt(1) + kvac1(ik) + qss1
     234           0 :          v(2) = bkpt(2) + kvac2(ik) + qss2
     235           0 :          v(3) = 0.
     236           0 :          ev = evacp - 0.5*dot_product(v,matmul(bbmat,v))
     237             : !         ev = evacp - 0.5*dotirp(v,v,bbmat)
     238             :          call vacuz(ev,vz(1,ivac,jspin2),vz0(ivac),nmz,delz,t(ik),
     239           0 :      +        dt(ik),u(1,ik))
     240             :          call vacudz(ev,vz(1,ivac,jspin2),vz0(ivac),nmz,delz,te(ik),
     241             :      +        dte(ik),tei(ik),ue(1,ik),dt(ik),
     242           0 :      +        u(1,ik))
     243           0 :          scale = wronk/ (te(ik)*dt(ik)-dte(ik)*t(ik))
     244           0 :          te(ik) = scale*te(ik)
     245           0 :          dte(ik) = scale*dte(ik)
     246           0 :          tei(ik) = scale*tei(ik)
     247           0 :          do j = 1,nmz
     248           0 :             ue(j,ik) = scale*ue(j,ik)
     249             :          enddo
     250             :       enddo
     251             : c-----> construct a and b coefficients for the k-point
     252             :        symvacvac=1
     253             : c       if (nvac==1) symvacvac=2
     254           0 :        do symvac=1,symvacvac
     255           0 :          do i = 1,nv2d
     256           0 :             do n = 1,nslibd
     257           0 :                ac(i,n) = cmplx(0.0,0.0)
     258           0 :                bc(i,n) = cmplx(0.0,0.0)
     259             :             enddo   
     260           0 :             do n = 1,nslibd_b
     261           0 :                ac_b(i,n) = cmplx(0.0,0.0)
     262           0 :                bc_b(i,n) = cmplx(0.0,0.0)
     263             :             enddo   
     264             :          enddo   
     265             : 
     266             :         if (symvac==2) sign=-1.0   
     267             : 
     268           0 :       do k = 1,nv(jspin)
     269           0 :          l = map2(k)
     270           0 :          zks = k3(k,jspin)*bmat(3,3)*sign
     271           0 :          arg = zks*z1
     272           0 :          c_1 = cmplx(cos(arg),sin(arg)) * const
     273           0 :          av = -c_1 * cmplx( dte(l),zks*te(l) )
     274           0 :          bv =  c_1 * cmplx(  dt(l),zks* t(l) )
     275             : c-----> loop over basis functions
     276             : 
     277           0 :          IF(zMat%l_real) THEN
     278           0 :             do n = 1,nslibd
     279           0 :                ac(l,n) = ac(l,n) + zMat%z_r(k+addnoco,n)*av
     280           0 :                bc(l,n) = bc(l,n) + zMat%z_r(k+addnoco,n)*bv
     281             :             enddo
     282             :          ELSE
     283           0 :             do n = 1,nslibd
     284           0 :                ac(l,n) = ac(l,n) + zMat%z_c(k+addnoco,n)*av
     285           0 :                bc(l,n) = bc(l,n) + zMat%z_c(k+addnoco,n)*bv
     286             :             enddo
     287             :          END IF
     288             :       enddo
     289             : 
     290             : 
     291             : c...now for the k+b point
     292             : 
     293           0 :       evacp = evac_b(ivac)
     294           0 :       do ik = 1,nv2_b
     295           0 :          v(1) = bkpt_b(1) + kvac1_b(ik) + qss1
     296           0 :          v(2) = bkpt_b(2) + kvac2_b(ik) + qss2
     297           0 :          v(3) = 0.
     298           0 :          ev = evacp - 0.5*dot_product(v,matmul(bbmat,v))
     299             : !         ev = evacp - 0.5*dotirp(v,v,bbmat)
     300             :          call vacuz(ev,vz(1,ivac,jspin2_b),vz0_b(ivac),nmz,delz,
     301           0 :      +        t_b(ik),dt_b(ik),u_b(1,ik))
     302             :          call vacudz(ev,vz(1,ivac,jspin2_b),vz0_b(ivac),nmz,delz,
     303             :      +        te_b(ik),dte_b(ik),tei_b(ik),ue_b(1,ik),dt_b(ik),
     304           0 :      +        u_b(1,ik))
     305           0 :          scale = wronk/ (te_b(ik)*dt_b(ik)-dte_b(ik)*t_b(ik))
     306           0 :          te_b(ik) = scale*te_b(ik)
     307           0 :          dte_b(ik) = scale*dte_b(ik)
     308           0 :          tei_b(ik) = scale*tei_b(ik)
     309           0 :          do j = 1,nmz
     310           0 :             ue_b(j,ik) = scale*ue_b(j,ik)
     311             :          enddo
     312             :       enddo
     313             : c-----> construct a and b coefficients for the k+b point
     314             : 
     315           0 :       do k = 1,nv_b(jspin_b)
     316           0 :          l = map2_b(k)
     317           0 :          zks = k3_b(k,jspin_b)*bmat(3,3)*sign
     318           0 :          arg = zks*z1
     319           0 :          c_1 = cmplx(cos(arg),sin(arg)) * const
     320           0 :          av = -c_1 * cmplx( dte_b(l),zks*te_b(l) )
     321           0 :          bv =  c_1 * cmplx( dt_b(l),zks*t_b(l) )
     322             : c-----> loop over basis functions
     323             : 
     324           0 :          IF(zMat_b%l_real) THEN
     325           0 :             do n = 1,nslibd_b
     326           0 :                ac_b(l,n) = ac_b(l,n) + zMat_b%z_r(k+addnoco2,n)*av
     327           0 :                bc_b(l,n) = bc_b(l,n) + zMat_b%z_r(k+addnoco2,n)*bv
     328             :             enddo
     329             :          ELSE
     330           0 :             do n = 1,nslibd_b
     331           0 :                ac_b(l,n) = ac_b(l,n) + zMat_b%z_c(k+addnoco2,n)*av
     332           0 :                bc_b(l,n) = bc_b(l,n) + zMat_b%z_c(k+addnoco2,n)*bv
     333             :             enddo
     334             :          END IF
     335             :       enddo
     336             : 
     337             : 
     338             :       ! set up z-mesh and plane-wave factors
     339           0 :       zks  =  gb(3) *bmat(3,3)   ! TODO: different sign before
     340           0 :       zks_b= -gb2(3)*bmat(3,3)   ! TODO: different sign before
     341           0 :       do i=1,nmz
     342           0 :        zmsh(i)= ( z1 + (i-1)*delz )*sign
     343           0 :        fac(i)  = cmplx( cos(zmsh(i)*zks), sin(zmsh(i)*zks) )
     344           0 :        fac_b(i)= cmplx( cos(zmsh(i)*zks_b),sin(zmsh(i)*zks_b))
     345           0 :        fac1(i) = fac(i) * fac_b(i)
     346             :       enddo
     347             : 
     348             : 
     349             :       ! calculate uHu matrix elements
     350           0 :       do l = 1,nv2
     351             : 
     352           0 :        j1 = kvac1(l) - gb(1)
     353           0 :        j2 = kvac2(l) - gb(2)
     354             : 
     355           0 :        do lp = 1,nv2_b
     356           0 :         i1 = j1 - kvac1_b(lp) + gb2(1)
     357           0 :         i2 = j2 - kvac2_b(lp) + gb2(2)
     358             : 
     359             : c        i1 = -i1
     360             : c        i2 = -i2
     361             : 
     362           0 :         ind3 = ig(sign2*i1,sign2*i2,0)
     363           0 :         IF (ind3.EQ.0) CYCLE
     364             : 
     365           0 :         phase = rgphs(i1,i2,0)   ! TODO: sign also here?
     366           0 :         phase2= rgphs(sign2*i1,sign2*i2,0)
     367           0 :         if(phase.ne.phase2) then
     368           0 :            stop 'rgphs in wann_uHu_vac'
     369             :         endif
     370             : 
     371           0 :         ind2 = ig2(ind3)
     372           0 :         IF (ind2.EQ.0) CALL juDFT_error("error in map2 for 2-d stars",
     373           0 :      >                                calledby="wann_uHu_vac")
     374           0 :         ind2 = ind2 - 1
     375           0 :         IF (ind2.NE.0) THEN ! warping components G=/=0
     376             :                             ! need integral with vxy and plane-waves
     377             :            ! tuu
     378           0 :            DO i=1,nmzxy
     379             :               xx(np0-i) = u(i,l)*u_b(i,lp)
     380           0 :      >                  * real( fac1(i)*vxy(i,ind2) )
     381             :               xximag(np0-i) = u(i,l)*u_b(i,lp)
     382           0 :      >                      * aimag( fac1(i)*vxy(i,ind2) )
     383             :            ENDDO
     384           0 :            CALL intgz0(xx,delz,nmzxy,xv,tail)
     385           0 :            CALL intgz0(xximag,delz,nmzxy,yv,tail)
     386           0 :            tuu(lp,l) = phase*cmplx(xv,yv)
     387             : 
     388             :            ! tud
     389           0 :            DO i=1,nmzxy
     390             :               xx(np0-i) = u(i,l)*ue_b(i,lp)
     391           0 :      >                  * real( fac1(i)*vxy(i,ind2) )
     392             :               xximag(np0-i) = u(i,l)*ue_b(i,lp)
     393           0 :      >                      * aimag( fac1(i)*vxy(i,ind2) )
     394             :            ENDDO
     395           0 :            CALL intgz0(xx,delz,nmzxy,xv,tail)
     396           0 :            CALL intgz0(xximag,delz,nmzxy,yv,tail)
     397           0 :            tud(lp,l) = phase*cmplx(xv,yv)
     398             : 
     399             :            ! tdu
     400           0 :            DO i=1,nmzxy
     401             :               xx(np0-i) = ue(i,l)*u_b(i,lp) 
     402           0 :      >                  * real( fac1(i)*vxy(i,ind2) )
     403             :               xximag(np0-i) = ue(i,l)*u_b(i,lp)
     404           0 :      >                      * aimag( fac1(i)*vxy(i,ind2) )
     405             :            ENDDO
     406           0 :            CALL intgz0(xx,delz,nmzxy,xv,tail)
     407           0 :            CALL intgz0(xximag,delz,nmzxy,yv,tail)
     408           0 :            tdu(lp,l) = phase*cmplx(xv,yv)
     409             : 
     410             :            ! tdd
     411           0 :            DO i=1,nmzxy
     412             :               xx(np0-i) = ue(i,l)*ue_b(i,lp)
     413           0 :      >                  * real( fac1(i)*vxy(i,ind2) )
     414             :               xximag(np0-i) = ue(i,l)*ue_b(i,lp)
     415           0 :      >                      * aimag( fac1(i)*vxy(i,ind2) )
     416             :            ENDDO
     417           0 :            CALL intgz0(xx,delz,nmzxy,xv,tail)
     418           0 :            CALL intgz0(xximag,delz,nmzxy,yv,tail)
     419           0 :            tdd(lp,l) = phase*cmplx(xv,yv)
     420             : 
     421             :         ELSE ! non-warping components G==0
     422             :              ! need integral with H(z) = -1/2 d2/dz2 + vz + (G+k)^2
     423             : 
     424           0 :            IF ( (ico.EQ.1) .OR. (ico.EQ.2) ) THEN ! spin-diagonal
     425             : 
     426             :               ! determine second derivative of (u*fac) etc.
     427             :               CALL d2fdz2cmplx(nmzd,nmz,zmsh,delz,
     428           0 :      >                         u(:,l),fac,du)
     429             :               CALL d2fdz2cmplx(nmzd,nmz,zmsh,delz,
     430           0 :      >                         ue(:,l),fac,due)
     431             :               CALL d2fdz2cmplx(nmzd,nmz,zmsh,delz,
     432           0 :      >                         u_b(:,lp),fac_b,du_b)
     433             :               CALL d2fdz2cmplx(nmzd,nmz,zmsh,delz,
     434           0 :      >                         ue_b(:,lp),fac_b,due_b)
     435             : 
     436             :               ! determine |G+k+b1|^2 and |G'+k+b2|^2
     437           0 :               s(1) = bkpt(1) + kvac1(l)
     438           0 :               s(2) = bkpt(2) + kvac2(l)
     439           0 :               s(3) = 0.0
     440           0 :               rk = dot_product(s,matmul(bbmat,s))
     441             : !              rk = dotirp(s,s,bbmat)
     442             : 
     443           0 :               s(1) = bkpt_b(1) + kvac1_b(lp)
     444           0 :               s(2) = bkpt_b(2) + kvac2_b(lp)
     445             :               s(3) = 0.0 
     446           0 :               rk_b = dot_product(s,matmul(bbmat,s))
     447             : !              rk_b = dotirp(s,s,bbmat)
     448             : 
     449             : c              ! no kinetic energy if jspin =/= jspin_b
     450             : c              if(jspin.ne.jspin_b) then
     451             : c               rk = 0.0
     452             : c               rk_b = 0.0
     453             : c               du=cmplx(0.0,0.0)
     454             : c               due=cmplx(0.0,0.0)
     455             : c               du_b=cmplx(0.0,0.0)
     456             : c               due_b=cmplx(0.0,0.0)
     457             : c              endif
     458             : 
     459             : c              rk_b=0.0
     460             : c              du_b=cmplx(0.,0.)
     461             : c              due_b=cmplx(0.,0.)
     462             : 
     463             :               ! tuu
     464           0 :               DO i=1,nmz
     465             :                xx(np1-i) = u(i,l)*u_b(i,lp) 
     466             :      >           * real(fac1(i)) * ( vz(i,ivac,jspin2H)+0.25*(rk+rk_b) )
     467             :      >           - 0.25*real(  fac(i)*u(i,l)*du_b(i) 
     468           0 :      >                       + du(i)*u_b(i,lp)*fac_b(i) )
     469             : 
     470             :                xximag(np1-i) = u(i,l)*u_b(i,lp) 
     471             :      >           * aimag(fac1(i)) *( vz(i,ivac,jspin2H)+0.25*(rk+rk_b) )
     472             :      >           - 0.25*aimag(  fac(i)*u(i,l)*du_b(i) 
     473           0 :      >                        + du(i)*u_b(i,lp)*fac_b(i) )
     474             :               ENDDO
     475           0 :               CALL intgz0(xx,delz,nmz,xv,tail)
     476           0 :               CALL intgz0(xximag,delz,nmz,yv,tail)
     477           0 :               tuu(lp,l) = cmplx(xv,yv)
     478             : 
     479             :               ! tud
     480           0 :               DO i=1,nmz
     481             :                xx(np1-i) = u(i,l)*ue_b(i,lp) 
     482             :      >           * real(fac1(i)) * ( vz(i,ivac,jspin2H)+0.25*(rk+rk_b) )
     483             :      >           - 0.25*real(  fac(i)*u(i,l)*due_b(i) 
     484           0 :      >                       + du(i)*ue_b(i,lp)*fac_b(i) )
     485             : 
     486             :                xximag(np1-i) = u(i,l)*ue_b(i,lp) 
     487             :      >           * aimag(fac1(i)) *( vz(i,ivac,jspin2H)+0.25*(rk+rk_b) )
     488             :      >           - 0.25*aimag(  fac(i)*u(i,l)*due_b(i) 
     489           0 :      >                        + du(i)*ue_b(i,lp)*fac_b(i) )
     490             :               ENDDO
     491           0 :               CALL intgz0(xx,delz,nmz,xv,tail)
     492           0 :               CALL intgz0(xximag,delz,nmz,yv,tail)
     493           0 :               tud(lp,l) = cmplx(xv,yv)
     494             : 
     495             :               ! tdu
     496           0 :               DO i=1,nmz
     497             :                xx(np1-i) = ue(i,l)*u_b(i,lp) 
     498             :      >           * real(fac1(i)) * ( vz(i,ivac,jspin2H)+0.25*(rk+rk_b) )
     499             :      >           - 0.25*real(  fac(i)*ue(i,l)*du_b(i) 
     500           0 :      >                       + due(i)*u_b(i,lp)*fac_b(i) )
     501             : 
     502             :                xximag(np1-i) = ue(i,l)*u_b(i,lp) 
     503             :      >           * aimag(fac1(i)) *( vz(i,ivac,jspin2H)+0.25*(rk+rk_b) )
     504             :      >           - 0.25*aimag(  fac(i)*ue(i,l)*du_b(i) 
     505           0 :      >                        + due(i)*u_b(i,lp)*fac_b(i) )
     506             :               ENDDO
     507           0 :               CALL intgz0(xx,delz,nmz,xv,tail)
     508           0 :               CALL intgz0(xximag,delz,nmz,yv,tail)
     509           0 :               tdu(lp,l) = cmplx(xv,yv)
     510             : 
     511             :               ! tdd
     512           0 :               DO i=1,nmz
     513             :                xx(np1-i) = ue(i,l)*ue_b(i,lp) 
     514             :      >           * real(fac1(i)) * ( vz(i,ivac,jspin2H)+0.25*(rk+rk_b) )
     515             :      >           - 0.25*real(  fac(i)*ue(i,l)*due_b(i) 
     516           0 :      >                       + due(i)*ue_b(i,lp)*fac_b(i) )
     517             : 
     518             :                xximag(np1-i) = ue(i,l)*ue_b(i,lp) 
     519             :      >           * aimag(fac1(i)) *( vz(i,ivac,jspin2H)+0.25*(rk+rk_b) )
     520             :      >           - 0.25*aimag(  fac(i)*ue(i,l)*due_b(i) 
     521           0 :      >                        + due(i)*ue_b(i,lp)*fac_b(i) )
     522             :               ENDDO
     523           0 :               CALL intgz0(xx,delz,nmz,xv,tail)
     524           0 :               CALL intgz0(xximag,delz,nmz,yv,tail)
     525           0 :               tdd(lp,l) = cmplx(xv,yv)
     526             : 
     527             :            ELSE ! spin-off-diagonal
     528             : 
     529           0 :               DO i=1,nmz
     530           0 :                  vv(i) = cmplx( vz(i,ivac,3), vz(i,ivac,4) )
     531             :               ENDDO
     532           0 :               IF ( ico.EQ.4 ) vv = conjg(vv)
     533             : 
     534             :               ! tuu
     535           0 :               DO i=1,nmz
     536             :                  xx(np1-i) = u(i,l)*u_b(i,lp)
     537           0 :      >                     * real( fac1(i)*vv(i) )
     538             :                  xximag(np1-i) = u(i,l)*u_b(i,lp)
     539           0 :      >                         * aimag( fac1(i)*vv(i) )
     540             :               ENDDO
     541           0 :               CALL intgz0(xx,delz,nmz,xv,tail)
     542           0 :               CALL intgz0(xximag,delz,nmz,yv,tail)
     543           0 :               tuu(lp,l) = cmplx(xv,yv)
     544             : 
     545             :               ! tud
     546           0 :               DO i=1,nmz
     547             :                  xx(np1-i) = u(i,l)*ue_b(i,lp)
     548           0 :      >                     * real( fac1(i)*vv(i) )
     549             :                  xximag(np1-i) = u(i,l)*ue_b(i,lp)
     550           0 :      >                         * aimag( fac1(i)*vv(i) )
     551             :               ENDDO
     552           0 :               CALL intgz0(xx,delz,nmz,xv,tail)
     553           0 :               CALL intgz0(xximag,delz,nmz,yv,tail)
     554           0 :               tud(lp,l) = cmplx(xv,yv)
     555             : 
     556             :               ! tdu
     557           0 :               DO i=1,nmz
     558             :                  xx(np1-i) = ue(i,l)*u_b(i,lp)
     559           0 :      >                     * real( fac1(i)*vv(i) )
     560             :                  xximag(np1-i) = ue(i,l)*u_b(i,lp)
     561           0 :      >                         * aimag( fac1(i)*vv(i) )
     562             :               ENDDO
     563           0 :               CALL intgz0(xx,delz,nmz,xv,tail)
     564           0 :               CALL intgz0(xximag,delz,nmz,yv,tail)
     565           0 :               tdu(lp,l) = cmplx(xv,yv)
     566             : 
     567             :               ! tdd
     568           0 :               DO i=1,nmz
     569             :                  xx(np1-i) = ue(i,l)*ue_b(i,lp)
     570           0 :      >                     * real( fac1(i)*vv(i) )
     571             :                  xximag(np1-i) = ue(i,l)*ue_b(i,lp)
     572           0 :      >                         * aimag( fac1(i)*vv(i) )
     573             :               ENDDO
     574           0 :               CALL intgz0(xx,delz,nmz,xv,tail)
     575           0 :               CALL intgz0(xximag,delz,nmz,yv,tail)
     576           0 :               tdd(lp,l) = cmplx(xv,yv)
     577             : 
     578             :            ENDIF !(ico.eq.1).or.(ico.eq.2)
     579             : 
     580             :         ENDIF !(ind2.ne.0)
     581             : 
     582             : 
     583             :         ! contraction of integrals with a,b coefficients
     584             :         ! yields contribution to uHu matrix        
     585             : c        do i = 1,nslibd
     586             : c           do j = 1,nslibd_b
     587             : c              uHu(i,j) = uHu(i,j) + area*
     588             : c     *            ( ac(l,i)*conjg(ac_b(lp,j))*tuu +
     589             : c     +              ac(l,i)*conjg(bc_b(lp,j))*tud +
     590             : c     *              bc(l,i)*conjg(ac_b(lp,j))*tdu +
     591             : c     +              bc(l,i)*conjg(bc_b(lp,j))*tdd )
     592             : c           enddo
     593             : c        enddo
     594             : 
     595             : 
     596             :        enddo ! lp
     597             :       enddo ! l
     598             : 
     599             :       if(.false.) then
     600             :       call CPP_BLAS_cgemm('C','N',nslibd_b,nv2d,nv2d,cmplx(area),
     601             :      >                    ac_b(1,1),nv2d,tuu(1,1),nv2d,cmplx(0.0),
     602             :      >                    mat(1,1),nslibd_b)
     603             :       call CPP_BLAS_cgemm('T','T',nslibd,nslibd_b,nv2d,chi,
     604             :      >                    ac(1,1),nv2d,mat(1,1),nslibd_b,cmplx(1.0),
     605             :      >                    uHu,nbnd)
     606             : 
     607             :       call CPP_BLAS_cgemm('C','N',nslibd_b,nv2d,nv2d,cmplx(area),
     608             :      >                    bc_b(1,1),nv2d,tud(1,1),nv2d,cmplx(0.0),
     609             :      >                    mat(1,1),nslibd_b)
     610             :       call CPP_BLAS_cgemm('T','T',nslibd,nslibd_b,nv2d,chi,
     611             :      >                    ac(1,1),nv2d,mat(1,1),nslibd_b,cmplx(1.0),
     612             :      >                    uHu,nbnd)
     613             : 
     614             :       call CPP_BLAS_cgemm('C','N',nslibd_b,nv2d,nv2d,cmplx(area),
     615             :      >                    ac_b(1,1),nv2d,tdu(1,1),nv2d,cmplx(0.0),
     616             :      >                    mat(1,1),nslibd_b)
     617             :       call CPP_BLAS_cgemm('T','T',nslibd,nslibd_b,nv2d,chi,
     618             :      >                    bc(1,1),nv2d,mat(1,1),nslibd_b,cmplx(1.0),
     619             :      >                    uHu,nbnd)
     620             : 
     621             :       call CPP_BLAS_cgemm('C','N',nslibd_b,nv2d,nv2d,cmplx(area),
     622             :      >                    bc_b(1,1),nv2d,tdd(1,1),nv2d,cmplx(0.0),
     623             :      >                    mat(1,1),nslibd_b)
     624             :       call CPP_BLAS_cgemm('T','T',nslibd,nslibd_b,nv2d,chi,
     625             :      >                    bc(1,1),nv2d,mat(1,1),nslibd_b,cmplx(1.0),
     626             :      >                    uHu,nbnd)
     627             :       endif
     628             : 
     629           0 :       if(.true.)then ! standard so far
     630             :       call CPP_BLAS_cgemm('T','N',nslibd_b,nv2d,nv2d,cmplx(area),
     631             :      >                    ac_b(1,1),nv2d,tuu(1,1),nv2d,cmplx(0.0),
     632           0 :      >                    mat(1,1),nslibd_b)
     633             :       call CPP_BLAS_cgemm('T','C',nslibd,nslibd_b,nv2d,chi,
     634             :      >                    ac(1,1),nv2d,mat(1,1),nslibd_b,cmplx(1.0),
     635           0 :      >                    uHu,nbnd)
     636             : 
     637             :       call CPP_BLAS_cgemm('T','N',nslibd_b,nv2d,nv2d,cmplx(area),
     638             :      >                    bc_b(1,1),nv2d,tud(1,1),nv2d,cmplx(0.0),
     639           0 :      >                    mat(1,1),nslibd_b)
     640             :       call CPP_BLAS_cgemm('T','C',nslibd,nslibd_b,nv2d,chi,
     641             :      >                    ac(1,1),nv2d,mat(1,1),nslibd_b,cmplx(1.0),
     642           0 :      >                    uHu,nbnd)
     643             : 
     644             :       call CPP_BLAS_cgemm('T','N',nslibd_b,nv2d,nv2d,cmplx(area),
     645             :      >                    ac_b(1,1),nv2d,tdu(1,1),nv2d,cmplx(0.0),
     646           0 :      >                    mat(1,1),nslibd_b)
     647             :       call CPP_BLAS_cgemm('T','C',nslibd,nslibd_b,nv2d,chi,
     648             :      >                    bc(1,1),nv2d,mat(1,1),nslibd_b,cmplx(1.0),
     649           0 :      >                    uHu,nbnd)
     650             : 
     651             :       call CPP_BLAS_cgemm('T','N',nslibd_b,nv2d,nv2d,cmplx(area),
     652             :      >                    bc_b(1,1),nv2d,tdd(1,1),nv2d,cmplx(0.0),
     653           0 :      >                    mat(1,1),nslibd_b)
     654             :       call CPP_BLAS_cgemm('T','C',nslibd,nslibd_b,nv2d,chi,
     655             :      >                    bc(1,1),nv2d,mat(1,1),nslibd_b,cmplx(1.0),
     656           0 :      >                    uHu,nbnd)
     657             :       endif
     658             : 
     659             :       enddo !symvac
     660             : c... cycle by the vacua finishes
     661           0 :  140  enddo      
     662             : 
     663           0 :       deallocate ( ac,bc,dt,dte,t,te,tei,u,ue,
     664           0 :      +             v,kvac1,kvac2,map2 )
     665           0 :       deallocate ( ac_b,bc_b,dt_b,dte_b,t_b,te_b,tei_b,u_b,ue_b,
     666           0 :      +             kvac1_b,kvac2_b,map2_b )
     667           0 :       deallocate ( vxy, vxy_help, vv, zmsh, fac1, fac, fac_b )
     668           0 :       deallocate ( du,due,du_b,due_b )
     669           0 :       deallocate ( tuu,tud,tdu,tdd,mat )
     670             : 
     671             : c      call fleur_end("stop")
     672             : 
     673           0 :       END SUBROUTINE wann_uHu_vac
     674             :       END MODULE m_wann_uHu_vac

Generated by: LCOV version 1.13