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

          Line data    Source code
       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3             : ! This file is part of FLEUR and available as free software under the conditions
       4             : ! of the MIT license as expressed in the LICENSE file in more detail.
       5             : !--------------------------------------------------------------------------------
       6             : 
       7             :       MODULE m_wann_mmkb_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) :: ig2(n3d),nv(jspd),nv_b(jspd)
      53             :       real,    intent (in) :: vz(nmzd,2),vz_b(nmzd,2)
      54             :       real,    intent (in) :: bbmat(3,3),bmat(3,3)
      55             :       integer, intent (in) :: k1(nvd,jspd),k2(nvd,jspd),k3(nvd,jspd)
      56             :       integer,intent(in)::k1_b(nvd,jspd),k2_b(nvd,jspd),k3_b(nvd,jspd)
      57             :       complex, intent (inout) :: mmn(nbnd,nbnd)
      58             : 
      59             : c     ..basis wavefunctions in the vacuum
      60           0 :       complex, allocatable :: ac(:,:),bc(:,:),ac_b(:,:),bc_b(:,:)
      61           0 :       real,    allocatable :: dt(:),dte(:)
      62           0 :       real,    allocatable :: t(:),te(:),tei(:)
      63           0 :       real,    allocatable :: u(:,:),ue(:,:),v(:)
      64           0 :       real,    allocatable :: dt_b(:),dte_b(:)
      65           0 :       real,    allocatable :: t_b(:),te_b(:),tei_b(:)
      66           0 :       real,    allocatable :: u_b(:,:),ue_b(:,:)
      67             : 
      68             : c     ..local scalars..
      69             :       logical tail
      70             :       real wronk,arg,zks,tpi,vz0(2),vz0_b(2),scale,evacp,ev,const
      71           0 :       real uu,ud,du,dd,xx(nmz),xximag(nmz)
      72             :       real :: uuimag,udimag,duimag,ddimag,qss1,qss2
      73             :       integer i,m,l,j,k,n,nv2,nv2_b,ivac,n2,n2_b,sign,ik
      74             :       integer :: lprime,np1,addnoco,addnoco2
      75             :       complex :: av,bv,ic,c_1
      76           0 :       integer, allocatable :: kvac1(:),kvac2(:),map2(:)
      77           0 :       integer, allocatable :: kvac1_b(:),kvac2_b(:),map2_b(:)
      78             :       integer symvac,symvacvac
      79             : 
      80             :       allocate ( ac(nv2d,nslibd),bc(nv2d,nslibd),
      81             :      +           ac_b(nv2d,nslibd_b),bc_b(nv2d,nslibd_b),
      82             :      +           dt(nv2d),dte(nv2d),t(nv2d),te(nv2d),
      83             :      +           tei(nv2d),u(nmzd,nv2d),ue(nmzd,nv2d),
      84             :      +           dt_b(nv2d),dte_b(nv2d),t_b(nv2d),te_b(nv2d),
      85             :      +           tei_b(nv2d),u_b(nmzd,nv2d),ue_b(nmzd,nv2d),
      86             :      +           v(3),kvac1(nv2d),kvac2(nv2d),map2(nvd),
      87           0 :      +           kvac1_b(nv2d),kvac2_b(nv2d),map2_b(nvd) )
      88             : 
      89           0 :       tpi = 2 * pimach() ; ic = cmplx(0.,1.)
      90             : 
      91           0 :       tail = .true.
      92           0 :       np1 = nmz + 1
      93             : 
      94             : c.. determining the indexing array (in-plane stars)
      95             : c.. for the k-point
      96             : 
      97           0 :       wronk = 2.0
      98           0 :       const = 1.0 / ( sqrt(omtil)*wronk )
      99             : 
     100           0 :       do ivac = 1,2
     101           0 :          vz0(ivac) = vz(nmz,ivac)
     102           0 :          vz0_b(ivac) = vz_b(nmz,ivac)
     103             :       enddo
     104             : 
     105           0 :       n2 = 0 ; n2_b = 0
     106             : 
     107           0 :       addnoco=0
     108           0 :       addnoco2=0
     109           0 :       if(l_noco.and.jspin.eq.2)then
     110           0 :         addnoco=nv(1)+nlotot
     111             :       endif
     112             :       if(l_noco.and.jspin_b.eq.2)then
     113             :         addnoco2=nv_b(1)+nlotot
     114             :       endif
     115             : 
     116           0 :       do 40 k = 1,nv(jspin)
     117           0 :          do 30 j = 1,n2
     118           0 :             if ( k1(k,jspin).eq.kvac1(j) .and.
     119             :      +          k2(k,jspin).eq.kvac2(j) ) then
     120           0 :                 map2(k) = j
     121           0 :                 goto 40
     122             :              endif 
     123           0 :  30      continue
     124           0 :          n2 = n2 + 1
     125             :          
     126           0 :          IF(n2>nv2d) then
     127           0 :             write(*,*)n2,nv2d,'jspin',jspin
     128             :          endif
     129             :  
     130           0 :          IF (n2>nv2d)  CALL juDFT_error("wannier Mmn vac",calledby
     131           0 :      +        ="wann_mmkb_vac")
     132             : 
     133           0 :          kvac1(n2) = k1(k,jspin)
     134           0 :          kvac2(n2) = k2(k,jspin)
     135           0 :          map2(k) = n2
     136           0 :  40   continue
     137             :          !write(*,*)'ok',n2,nv2d,'jspin',jspin
     138             : 
     139             : c.. and for the b-point
     140             :  
     141           0 :       do 41 k = 1,nv_b(jspin_b)
     142           0 :          do 31 j = 1,n2_b
     143           0 :             if ( k1_b(k,jspin_b).eq.kvac1_b(j) .and.
     144             :      +          k2_b(k,jspin_b).eq.kvac2_b(j) ) then
     145           0 :                 map2_b(k) = j
     146           0 :                 goto 41
     147             :              endif
     148           0 :  31      continue
     149           0 :          n2_b = n2_b + 1
     150             :      
     151           0 :          IF(n2_b>nv2d) then
     152           0 :             write(*,*)n2_b,nv2d,'jspin_b',jspin_b
     153             :          endif
     154             : 
     155           0 :          IF (n2_b>nv2d)  CALL juDFT_error("wannier Mmn vac",calledby
     156           0 :      +        ="wann_mmkb_vac")
     157             : 
     158           0 :          kvac1_b(n2_b) = k1_b(k,jspin_b)
     159           0 :          kvac2_b(n2_b) = k2_b(k,jspin_b)
     160           0 :          map2_b(k) = n2_b
     161           0 :  41   continue
     162             :          !write(*,*)'ok',n2_b,nv2d,'jspin_b',jspin_b
     163             : 
     164             : c...cycle by the vacua
     165           0 :       do 140 ivac = 1,nvac
     166             : 
     167             : 
     168             : 
     169           0 :       sign = 3. - 2.*ivac
     170           0 :       evacp = evac(ivac)
     171             : 
     172           0 :       nv2 = n2 ; nv2_b = n2_b
     173             : 
     174             : c.. the body of the routine
     175             : 
     176           0 :       qss1=0.0
     177           0 :       qss2=0.0
     178           0 :       if(l_noco.and.jspin.eq.1)then
     179           0 :         qss1=-qss(1)/2.0
     180           0 :         qss2=-qss(2)/2.0
     181           0 :       elseif(l_noco.and.jspin.eq.2)then
     182           0 :         qss1=qss(1)/2.0
     183           0 :         qss2=qss(2)/2.0
     184             :       endif
     185             : 
     186           0 :       do ik = 1,nv2
     187           0 :          v(1) = bkpt(1) + kvac1(ik) + qss1
     188           0 :          v(2) = bkpt(2) + kvac2(ik) + qss2
     189           0 :          v(3) = 0.
     190           0 :          ev = evacp - 0.5*dot_product(v,matmul(bbmat,v))
     191             :          call vacuz(ev,vz(1,ivac),vz0(ivac),nmz,delz,t(ik),dt(ik),
     192           0 :      +        u(1,ik))
     193             :          call vacudz(ev,vz(1,ivac),vz0(ivac),nmz,delz,te(ik),
     194             :      +        dte(ik),tei(ik),ue(1,ik),dt(ik),
     195           0 :      +        u(1,ik))
     196           0 :          scale = wronk/ (te(ik)*dt(ik)-dte(ik)*t(ik))
     197           0 :          te(ik) = scale*te(ik)
     198           0 :          dte(ik) = scale*dte(ik)
     199           0 :          tei(ik) = scale*tei(ik)
     200           0 :          do j = 1,nmz
     201           0 :             ue(j,ik) = scale*ue(j,ik)
     202             :          enddo
     203             :       enddo
     204             : c-----> construct a and b coefficients for the k-point
     205           0 :        symvacvac=1
     206           0 :        if (nvac==1) symvacvac=2
     207           0 :        do symvac=1,symvacvac
     208           0 :          do i = 1,nv2d
     209           0 :             do n = 1,nslibd
     210           0 :                ac(i,n) = cmplx(0.0,0.0)
     211           0 :                bc(i,n) = cmplx(0.0,0.0)
     212             :             enddo   
     213           0 :             do n = 1,nslibd_b
     214           0 :                ac_b(i,n) = cmplx(0.0,0.0)
     215           0 :                bc_b(i,n) = cmplx(0.0,0.0)
     216             :             enddo   
     217             :          enddo   
     218             : 
     219           0 :         if (symvac==2) sign=-1.0   
     220             : 
     221           0 :       do k = 1,nv(jspin)
     222           0 :          l = map2(k)
     223           0 :          zks = k3(k,jspin)*bmat(3,3)*sign
     224           0 :          arg = zks*z1
     225           0 :          c_1 = cmplx(cos(arg),sin(arg)) * const
     226           0 :          av = -c_1 * cmplx( dte(l),zks*te(l) )
     227           0 :          bv =  c_1 * cmplx(  dt(l),zks* t(l) )
     228             : c-----> loop over basis functions
     229           0 :          IF(zMat%l_real) THEN
     230           0 :             do n = 1,nslibd
     231           0 :                ac(l,n) = ac(l,n) + zMat%data_r(k+addnoco,n)*av
     232           0 :                bc(l,n) = bc(l,n) + zMat%data_r(k+addnoco,n)*bv
     233             :             enddo
     234             :          ELSE
     235           0 :             do n = 1,nslibd
     236           0 :                ac(l,n) = ac(l,n) + zMat%data_c(k+addnoco,n)*av
     237           0 :                bc(l,n) = bc(l,n) + zMat%data_c(k+addnoco,n)*bv
     238             :             enddo
     239             :          END IF
     240             :       enddo
     241             : 
     242             : c...now for the k+b point
     243             : 
     244           0 :       evacp = evac_b(ivac)
     245           0 :       do ik = 1,nv2_b
     246           0 :          v(1) = bkpt_b(1) + kvac1_b(ik) + qss1
     247           0 :          v(2) = bkpt_b(2) + kvac2_b(ik) + qss2
     248           0 :          v(3) = 0.
     249           0 :          ev = evacp - 0.5*dot_product(v,matmul(bbmat,v))
     250             :          call vacuz(ev,vz_b(1,ivac),vz0_b(ivac),nmz,delz,t_b(ik),
     251           0 :      +        dt_b(ik),u_b(1,ik))
     252             :          call vacudz(ev,vz_b(1,ivac),vz0_b(ivac),nmz,delz,te_b(ik),
     253             :      +        dte_b(ik),tei_b(ik),ue_b(1,ik),dt_b(ik),
     254           0 :      +        u_b(1,ik))
     255           0 :          scale = wronk/ (te_b(ik)*dt_b(ik)-dte_b(ik)*t_b(ik))
     256           0 :          te_b(ik) = scale*te_b(ik)
     257           0 :          dte_b(ik) = scale*dte_b(ik)
     258           0 :          tei_b(ik) = scale*tei_b(ik)
     259           0 :          do j = 1,nmz
     260           0 :             ue_b(j,ik) = scale*ue_b(j,ik)
     261             :          enddo
     262             :       enddo
     263             : c-----> construct a and b coefficients for the k+b point
     264             : 
     265           0 :       do k = 1,nv_b(jspin_b)
     266           0 :          l = map2_b(k)
     267           0 :          zks = k3_b(k,jspin_b)*bmat(3,3)*sign
     268           0 :          arg = zks*z1
     269           0 :          c_1 = cmplx(cos(arg),sin(arg)) * const
     270           0 :          av = -c_1 * cmplx( dte_b(l),zks*te_b(l) )
     271           0 :          bv =  c_1 * cmplx( dt_b(l),zks*t_b(l) )
     272             : c-----> loop over basis functions
     273           0 :          IF(zMat_b%l_real) THEN
     274           0 :             do n = 1,nslibd_b
     275           0 :                ac_b(l,n) = ac_b(l,n) + zMat_b%data_r(k+addnoco,n)*av
     276           0 :                bc_b(l,n) = bc_b(l,n) + zMat_b%data_r(k+addnoco,n)*bv
     277             :             enddo
     278             :          ELSE
     279           0 :             do n = 1,nslibd_b
     280           0 :                ac_b(l,n) = ac_b(l,n) + zMat_b%data_c(k+addnoco,n)*av
     281           0 :                bc_b(l,n) = bc_b(l,n) + zMat_b%data_c(k+addnoco,n)*bv
     282             :             enddo
     283             :          END IF
     284             :       enddo
     285             : 
     286             : 
     287           0 :       do l = 1,nv2
     288           0 :       do lprime = 1,nv2_b
     289             :       if (kvac1(l).eq.(kvac1_b(lprime)-gb(1))
     290           0 :      & .and. kvac2(l).eq.(kvac2_b(lprime)-gb(2)))then
     291           0 :          zks = gb(3)*bmat(3,3)*sign
     292             : 
     293           0 :          do i = 1,nmz
     294             :              xx(np1-i) = u(i,l)*u_b(i,lprime)*
     295           0 :      *          cos((z1+(i-1)*delz)*zks)
     296             :             xximag(np1-i) = u(i,l)*u_b(i,lprime)*
     297           0 :      *          sin((z1+(i-1)*delz)*zks)
     298             :          enddo   
     299           0 :          call intgz0(xx,delz,nmz,uu,tail)
     300           0 :          call intgz0(xximag,delz,nmz,uuimag,tail)
     301             : 
     302           0 :          do i = 1,nmz
     303             :             xx(np1-i) = u(i,l)*ue_b(i,lprime)*
     304           0 :      *   cos((z1+(i-1)*delz)*zks)
     305             :             xximag(np1-i) = u(i,l)*ue_b(i,lprime)*
     306           0 :      *          sin((z1+(i-1)*delz)*zks)
     307             :          enddo   
     308           0 :          call intgz0(xx,delz,nmz,ud,tail)
     309           0 :          call intgz0(xximag,delz,nmz,udimag,tail)
     310             : 
     311           0 :          do i = 1,nmz
     312             :             xx(np1-i) = ue(i,l)*u_b(i,lprime)*
     313           0 :      *   cos((z1+(i-1)*delz)*zks)
     314             :             xximag(np1-i) = ue(i,l)*u_b(i,lprime)*
     315           0 :      *          sin((z1+(i-1)*delz)*zks)
     316             :          enddo   
     317           0 :          call intgz0(xx,delz,nmz,du,tail)
     318           0 :          call intgz0(xximag,delz,nmz,duimag,tail)
     319           0 :          do i = 1,nmz
     320             :             xx(np1-i) = ue(i,l)*ue_b(i,lprime)*
     321           0 :      *   cos((z1+(i-1)*delz)*zks)
     322             :             xximag(np1-i) = ue(i,l)*ue_b(i,lprime)*
     323           0 :      *          sin((z1+(i-1)*delz)*zks)
     324             :          enddo   
     325           0 :          call intgz0(xx,delz,nmz,dd,tail)
     326           0 :          call intgz0(xximag,delz,nmz,ddimag,tail)
     327             : 
     328           0 :          do i = 1,nslibd
     329           0 :             do j = 1,nslibd_b
     330             :                mmn(i,j) = mmn(i,j) + area*(
     331             :      *  ac(l,i)*conjg(ac_b(lprime,j))*cmplx(uu,uuimag) +
     332             :      +  ac(l,i)*conjg(bc_b(lprime,j))*cmplx(ud,udimag) +
     333             :      *  bc(l,i)*conjg(ac_b(lprime,j))*cmplx(du,duimag) +
     334           0 :      +  bc(l,i)*conjg(bc_b(lprime,j))*cmplx(dd,ddimag) )*vacchi
     335             :             enddo
     336             :          enddo
     337             : 
     338             : 
     339             : 
     340             : 
     341             : 
     342             :       endif ! kvac1=kvac1_b and kvac2=kvac2_b
     343             :       enddo ! lprime
     344             :       enddo ! l
     345             : 
     346             :       enddo !symvac
     347             : c... cycle by the vacua finishes
     348           0 :  140  enddo      
     349             : 
     350           0 :       deallocate ( ac,bc,dt,dte,t,te,tei,u,ue,
     351           0 :      +             v,kvac1,kvac2,map2 )
     352           0 :       deallocate ( ac_b,bc_b,dt_b,dte_b,t_b,te_b,tei_b,u_b,ue_b,
     353           0 :      +             kvac1_b,kvac2_b,map2_b )
     354             : 
     355           0 :       END SUBROUTINE wann_mmkb_vac
     356             :       END MODULE m_wann_mmkb_vac

Generated by: LCOV version 1.13