LCOV - code coverage report
Current view: top level - wannier/uhu - wann_uHu_commat.f (source / functions) Hit Total Coverage
Test: FLEUR test coverage Lines: 0 207 0.0 %
Date: 2024-03-29 04:21:46 Functions: 0 1 0.0 %

          Line data    Source code
       1             :       module m_wann_uHu_commat
       2             :       USE m_juDFT
       3             :       implicit none
       4             :       contains
       5             : 
       6           0 :       subroutine wann_uHu_commat(nkpts,nntot,bpt,nqpts,nntot_q,
       7           0 :      >                           bpt_q,gb,gb_q,latt_const_q,
       8             :      >                           l_unformatted,l_dim,nparampts,
       9             :      >                           param_vec)
      10             :       use m_wann_gwf_tools
      11             :       use m_wann_gwf_auxovlp
      12             :       use m_constants, only : pimach
      13             : 
      14             :       implicit none
      15             :       
      16             :       ! input parameters
      17             :       logical,intent(in) :: l_unformatted,l_dim(3)
      18             :       integer,intent(in) :: nkpts,nntot,nqpts,nntot_q,nparampts
      19             :       integer,intent(in) :: bpt(nntot,nkpts),bpt_q(nntot_q,nqpts)
      20             :       integer,intent(in) :: gb(3,nntot,nkpts),gb_q(3,nntot_q,nqpts)
      21             :       real,   intent(in) :: latt_const_q,param_vec(3,nparampts)
      22             : 
      23             :       ! allocatable arrays
      24           0 :       logical,allocatable :: l_fm(:),l_fm2(:)
      25           0 :       integer,allocatable :: gb_kq(:,:,:)
      26           0 :       integer,allocatable :: bpt_kq(:,:)
      27           0 :       integer,allocatable :: g1(:,:,:),g2(:,:)
      28           0 :       integer,allocatable :: b1(:,:),b2(:)
      29             :       complex,allocatable :: tmp_a(:,:,:),tmp_m(:,:,:,:)
      30           0 :       complex,allocatable :: uHukk(:,:,:,:,:)
      31           0 :       complex,allocatable :: uHukq(:,:,:,:,:)
      32           0 :       complex,allocatable :: mkk(:,:,:,:)
      33           0 :       complex,allocatable :: mkq(:,:,:,:)
      34           0 :       complex,allocatable :: uHu(:,:,:)
      35           0 :       character(len=20),allocatable :: fuHu(:),fuHu2(:)
      36             : 
      37             :       ! further variables
      38             :       integer :: nbnd,nwfs,kb,qb,kb2,qb2,nn_kq,kqb,kqb2,arr_len,shift(3)
      39             :       integer :: i,j,ikqpt,ikqpt_b,nk,nq,nk2,nq2,g(3)
      40             :       integer :: nbnd_arti,nqpts_arti,nntot_arti,nnkq_len,nnkq_ind
      41             :       integer :: nkqpts,nntot_kq,nn,nn2,ikqpt_help,testnn
      42             :       integer :: dump,nn_arti,nkp,k,b,q,nwf,n,is,ie,countkq,countqk
      43             :       real :: mmn_r,mmn_i,tau,tpi,a_arti,b_arti,eig
      44             :       complex :: ovaux
      45             :       logical :: l_exist,l_miss,l_fe,l_fa,l_proj
      46             :       character(len=20) :: fq
      47             :       character(len=32) :: fmt,fmt2
      48             : 
      49           0 :       write(*,*)'create uHu for HDWF case...'
      50             : 
      51           0 :       call get_dimension(l_dim,arr_len)
      52           0 :       call get_shift(l_dim,shift)
      53             :       
      54           0 :       write(*,*)'dimension:',arr_len
      55           0 :       write(*,*)'x?,y?,z? :',l_dim(1:3)
      56           0 :       if(arr_len.le.3) call juDFT_error("dimension<4",
      57           0 :      >                                  calledby='wann_gwf_commat')
      58             : 
      59           0 :       nkqpts=nkpts*nqpts
      60           0 :       tpi = 2.0*pimach()
      61           0 :       a_arti = latt_const_q
      62           0 :       b_arti = 0.98*latt_const_q
      63           0 :       tau = tpi/real(nqpts)/a_arti
      64           0 :       ovaux = 2.0*tpi*tpi*sin(tau*b_arti/2.0)
      65           0 :       ovaux = ovaux/(tpi*tpi-tau*tau*b_arti*b_arti)
      66           0 :       ovaux = ovaux/(tau*b_arti)
      67           0 :       write(*,*)'ov(0)=',ovaux
      68             : 
      69           0 :       allocate(fuHu(nqpts),fuHu2(nqpts),l_fm(nqpts),l_fm2(nqpts))
      70             : 
      71             : ! are all necessary files present?
      72           0 :       l_miss=.false.
      73           0 :       do q=1,nqpts
      74           0 :          write(fq,'(i4.4)')q
      75           0 :          fuHu(q) = 'WF1_'//trim(fq)//'.uHu'    ! kk-overlaps
      76           0 :          fuHu2(q)= 'WF1_'//trim(fq)//'.uHu_kq'  ! kq-overlaps
      77             : 
      78           0 :          inquire(file=fuHu(q) ,exist=l_fm(q) )
      79           0 :          inquire(file=fuHu2(q),exist=l_fm2(q))
      80             : 
      81           0 :          if(.not.l_fm(q) ) write(*,*)'missing: ',fuHu(q)
      82           0 :          if(.not.l_fm2(q)) write(*,*)'missing: ',fuHu2(q)
      83           0 :          if(.not.(l_fm(q).or.l_fm2(q))) l_miss=.true.
      84             :       enddo
      85           0 :       inquire(file='bkqpts',exist=l_exist)
      86           0 :       if(.not.l_exist) write(*,*)'missing: bkqpts'
      87           0 :       inquire(file='proj',exist=l_proj)
      88           0 :       if(.not.l_proj) write(*,*)'missing: proj'
      89           0 :       if(l_miss.or.(.not.l_exist).or.(.not.l_proj))
      90           0 :      >   call juDFT_error("missing file(s) for HDWFs")
      91             : 
      92             : ! get number of bands and wfs from proj
      93           0 :       open(405,file='proj',status='old')
      94           0 :       read(405,*)nwfs,nbnd
      95           0 :       close(405)
      96           0 :       write(*,*)'nbnd=',nbnd
      97           0 :       write(*,*)'nwfs=',nwfs
      98             : 
      99             : 
     100             : c*****************************c
     101             : c     COMPOSITE .MMN FILE     c
     102             : c*****************************c
     103             : 
     104           0 :       write(fmt,'(a,i1,a)')'(2i6,3x,',arr_len,'i4)'
     105           0 :       write(fmt2,'(a,i1,a)')'(i7,i7,3x,',arr_len,'i4)'
     106             : 
     107           0 :       open (202,file='bkqpts',form='formatted',status='old')
     108           0 :       rewind (202)
     109           0 :       read (202,'(i4)') nntot_kq
     110           0 :       write (*,*) 'nntot_kq=',nntot_kq
     111           0 :       allocate ( gb_kq(arr_len,nntot_kq,nkqpts),
     112           0 :      &           bpt_kq(nntot_kq,nkqpts))
     113           0 :       do ikqpt=1,nkqpts
     114           0 :        do nn=1,nntot_kq
     115             :           read (202,fmt)
     116           0 :      &      ikqpt_help,bpt_kq(nn,ikqpt),(gb_kq(i,nn,ikqpt),i=1,arr_len)
     117             :        enddo
     118             :       enddo
     119           0 :       close (202)
     120             : 
     121           0 :       if(l_unformatted) then
     122           0 :          nnkq_len = (nntot_kq+1)*nntot_kq/2
     123           0 :          allocate(uHukk(nbnd,nbnd,nntot,nntot,nkpts))
     124           0 :          allocate(uHukq(nbnd,nbnd,nntot_q,nntot,nkpts))
     125           0 :          allocate(g1(3,nntot,nkpts))
     126           0 :          allocate(b1(nntot,nkpts))
     127             : c         allocate(uHu(nbnd,nbnd,nntot_kq,nntot_kq))!,nkqpts))
     128           0 :          allocate(uHu(nbnd,nbnd,nnkq_len))
     129             : 
     130           0 :          open(307,file='WF1_gwf.uHu',form='unformatted')
     131           0 :          write(307)nbnd,nkqpts,nntot_kq,nnkq_len
     132           0 :          write(307)bpt_kq,gb_kq
     133             : 
     134             :          !open(853,file='debug_commat')
     135             : 
     136           0 :          countkq=0
     137           0 :          countqk=0
     138             : 
     139           0 :          do q=1,nqpts
     140             :            is = get_index_kq(1,q,nkpts)
     141             :            ie = get_index_kq(nkpts,q,nkpts)
     142             :          
     143           0 :            if(l_fm(q)) then
     144           0 :             open(405,file=fuHu(q),form='unformatted')
     145           0 :             read(405)b,k,nn
     146           0 :             read(405)b1,g1
     147           0 :             read(405)uHukk
     148           0 :             close(405,status='delete')
     149             :            else
     150           0 :             uHukk=cmplx(0.,0.)
     151           0 :             write(*,*)'k-k overlaps set to zero for q=',q
     152             :            endif
     153             : 
     154           0 :            if(l_fm2(q)) then
     155           0 :             open(405,file=fuHu2(q),form='unformatted')
     156           0 :             read(405)b,k,nn
     157           0 :             read(405)b1,g1
     158           0 :             read(405)uHukq
     159           0 :             close(405,status='delete')
     160             :            else
     161           0 :             uHukq=cmplx(0.,0.)
     162           0 :             write(*,*)'k-q overlaps set to zero for q=',q
     163             :            endif
     164             : 
     165             : c           mmnq = mmnq*conjg(mmn_arti)
     166             : 
     167           0 :            do k=1,nkpts
     168           0 :              uHu = cmplx(0.,0.)
     169           0 :              ikqpt=get_index_kq(k,q,nkpts)
     170           0 :              nnkq_ind=0 
     171           0 :              do nn=1,nntot_kq
     172           0 :                kqb = bpt_kq(nn,ikqpt)
     173           0 :                kb=get_index_k(kqb,nkpts)
     174           0 :                qb=get_index_q(kqb,nkpts)
     175             : 
     176             :                nk=get_index_nn_k(bpt(:,k),nntot,kb,
     177             :      >                         gb_kq(:,nn,ikqpt),
     178           0 :      >                         gb(1:3,1:nntot,k),arr_len)
     179             :                nq=get_index_nn_q(bpt_q(:,q),nntot_q,qb,
     180             :      >                         gb_kq(:,nn,ikqpt),
     181             :      >                         gb_q(1:3,1:nntot_q,q),
     182           0 :      >                         arr_len,shift,l_dim)
     183             : 
     184           0 :                do nn2=1,nn!nntot_kq
     185           0 :                nnkq_ind = nnkq_ind+1
     186           0 :                kqb2= bpt_kq(nn2,ikqpt)
     187           0 :                kb2=get_index_k(kqb2,nkpts)
     188           0 :                qb2=get_index_q(kqb2,nkpts)
     189             : 
     190             :                nk2=get_index_nn_k(bpt(:,k),nntot,kb2,
     191             :      >                         gb_kq(:,nn2,ikqpt),
     192           0 :      >                         gb(1:3,1:nntot,k),arr_len)
     193             :                nq2=get_index_nn_q(bpt_q(:,q),nntot_q,qb2,
     194             :      >                         gb_kq(:,nn2,ikqpt),
     195             :      >                         gb_q(1:3,1:nntot_q,q),
     196           0 :      >                         arr_len,shift,l_dim)
     197           0 :                testnn = nk/abs(nk)+nk2/abs(nk2)+nq/abs(nq)+nq2/abs(nq2)
     198           0 :                if(testnn.ne.0) stop 'testnn.ne.0'
     199             : c               write(*,*)'kq',ikqpt,nn,nn2
     200             : c               write(*,*)'k ',k,nk,nk2
     201             : c               write(*,*)'q ',q,nq,nq2
     202             : 
     203             :                !write(853,'(a,1x,5(5i))')'kq,nkq',ikqpt,kqb,kqb2,nn,nn2
     204             :                !write(853,'(a,1x,5(5i))')'k,nk',k,kb,kb2,nk,nk2
     205             :                !write(853,'(a,1x,5(5i))')'q,nq',q,qb,qb2,nq,nq2
     206             :                !write(853,*)'nnkq_ind',nnkq_ind,nnkq_len
     207             : 
     208             : c               CALL wann_gwf_auxovlp(param_vec(:,qb),param_vec(:,qb2),
     209             : c     >                               latt_const_q,ovaux)
     210             : c               if(ovaux.ne.1.0) write(*,*)'ovaux',ovaux
     211             : 
     212           0 :                if((nq.lt.0) .and. (nq2.lt.0)) then ! kk
     213           0 :                   uHu(:,:,nnkq_ind)=uHukk(:,:,nk2,nk,k)!*ovaux
     214             : c                  write(*,'(a,5i8)')'<k|k>'
     215             :                   !write(853,*)'<k|k>'
     216           0 :                elseif((nk.lt.0) .and. (nk2.lt.0)) then ! qq
     217           0 :                   uHu(:,:,nnkq_ind)=cmplx(0.,0.)
     218             : c                  write(*,'(a,5i8)')'<q|q>'
     219             :                   !write(853,*)'<q|q>'
     220             :                else
     221             :                   if((nk.lt.0) .and. (nk2.ge.0) .and.
     222           0 :      >               (nq.ge.0) .and. (nq2.lt.0) ) then
     223             :                      uHu(:,:,nnkq_ind)
     224             :      >                  = transpose(conjg(uHukq(:,:,nq,nk2,k)))
     225           0 :      >                    *conjg(ovaux)
     226             :                   !write(853,*)'<q|k>'
     227             : c                     write(*,'(a,8i8)')'<q|k>'
     228             : c                  write(*,*)nn,nn2
     229             : c                  write(*,*)'<q|k>',q,bpt_q(nq,q),k,bpt(nk2,k)
     230           0 :                   if(q.eq.bpt_q(nq,q)) stop 'q.eq.bpt_q'
     231           0 :                   countqk=countqk+1
     232             :                   elseif((nk.ge.0) .and. (nk2.lt.0) .and. 
     233           0 :      >                   (nq.lt.0) .and. (nq2.ge.0) ) then
     234             :                      uHu(:,:,nnkq_ind)
     235           0 :      >                  = uHukq(:,:,nq2,nk,k) * ovaux
     236             :                   !write(853,*)'<k|q>'
     237             : c                     write(*,'(a,8i8)')'<k|q>'
     238             : c                  write(*,*)nn,nn2
     239             : c                  write(*,*)'<k|q>',k,bpt(nk,k),q,bpt_q(nq2,q)
     240           0 :                   if(q.eq.bpt_q(nq2,q)) stop 'q.eq.bpt_q'
     241           0 :                   countkq=countkq+1
     242             :                   else
     243             :                      call juDFT_error("problem k,q neighbors",
     244           0 :      >                                calledby="wann_uHu_commat")
     245             :                   endif
     246             :                endif
     247             :                enddo
     248             : 
     249             :              enddo!nn
     250           0 :              write(307)ikqpt,uHu
     251             :            enddo!k
     252             :          enddo!q
     253           0 :          deallocate(b1,g1,uHukk,uHukq)
     254             : 
     255             :          !close(853)
     256           0 :          close(307)
     257           0 :          deallocate(uHu)
     258             : 
     259           0 :          write(*,*)'countkq',countkq
     260           0 :          write(*,*)'countqk',countqk
     261             : 
     262             :       else
     263           0 :          allocate(mkk(nbnd,nbnd,nntot,nntot))
     264           0 :          allocate(mkq(nbnd,nbnd,nntot_q,nntot))
     265           0 :          open(307,file='WF1_gwf.uHu')
     266           0 :          write(307,*)'uHu for HDWFs'
     267           0 :          write(307,'(i5,i7,i5)')nbnd,nkqpts,nntot_kq
     268             : 
     269           0 :          do q=1,nqpts      
     270             :  
     271           0 :            if(l_fm(q)) then
     272           0 :             open(405,file=fuHu(q))
     273           0 :             read(405,*)!title
     274           0 :             read(405,*)b,k,nn
     275             :            endif
     276           0 :            if(l_fm2(q)) then
     277           0 :             open(505,file=fuHu2(q))
     278           0 :             read(505,*)!title
     279           0 :             read(505,*)b,k,nn
     280             :            endif
     281             : 
     282           0 :            do k=1,nkpts
     283           0 :              ikqpt=get_index_kq(k,q,nkpts)           
     284             : 
     285           0 :              mkk = cmplx(0.,0.)
     286           0 :              if(l_fm(q)) then
     287           0 :               do nn=1,nntot
     288           0 :                do nn2=1,nntot
     289           0 :                 read(405,*)nkp,b,i!g(1:3)
     290           0 :                 do i=1,nbnd
     291           0 :                    do j=1,nbnd
     292           0 :                       read(405,*)mmn_r,mmn_i
     293           0 :                       mkk(j,i,nn2,nn)=cmplx(mmn_r,mmn_i)
     294             :                    enddo
     295             :                 enddo              
     296             :                enddo 
     297             :               enddo
     298             :              endif
     299             : 
     300           0 :              mkq = cmplx(0.,0.)
     301           0 :              if(l_fm2(q)) then
     302           0 :               do nn=1,nntot
     303           0 :                do nn2=1,nntot_q
     304           0 :                 read(505,*)nkp,b,i!g(1:3)
     305           0 :                 do i=1,nbnd
     306           0 :                    do j=1,nbnd
     307           0 :                       read(505,*)mmn_r,mmn_i
     308           0 :                       mkq(j,i,nn2,nn)=cmplx(mmn_r,mmn_i)
     309             :                    enddo
     310             :                 enddo              
     311             :                enddo 
     312             :               enddo
     313           0 :               mkq = mkq*ovaux
     314             :              endif
     315             : 
     316           0 :              do nn=1,nntot_kq ! write composite overlaps
     317           0 :                kqb = bpt_kq(nn,ikqpt)
     318           0 :                kb=get_index_k(kqb,nkpts)
     319           0 :                qb=get_index_q(kqb,nkpts)
     320             : 
     321           0 :                do nn2=1,nntot_kq
     322           0 :                kqb2 = bpt_kq(nn2,ikqpt)
     323           0 :                kb2=get_index_k(kqb2,nkpts)
     324           0 :                qb2=get_index_q(kqb2,nkpts)
     325             : 
     326           0 :                write(307,'(3i7)')ikqpt,kqb,kqb2!gb_kq(1:arr_len,nn,ikqpt)
     327           0 :                if((q.eq.qb).and.(q.eq.qb2)) then ! k-overlap
     328             :                   !write(307,*)'<k|k>'
     329             :                   nk=get_index_nn_k(bpt(:,k),nntot,kb,
     330             :      >                         gb_kq(:,nn,ikqpt),
     331           0 :      >                         gb(1:3,1:nntot,k),arr_len)
     332             :                   nk2=get_index_nn_k(bpt(:,k),nntot,kb2,
     333             :      >                         gb_kq(:,nn2,ikqpt),
     334           0 :      >                         gb(1:3,1:nntot,k),arr_len)
     335           0 :                   do i=1,nbnd
     336           0 :                      do j=1,nbnd
     337             :                         write(307,'(2f24.18)')
     338           0 :      >                      real(mkk(j,i,nk2,nk)),aimag(mkk(j,i,nk2,nk))
     339             :                      enddo
     340             :                   enddo
     341           0 :                elseif((k.eq.kb).and.(k.eq.kb2)) then ! q-overlap
     342             : 
     343             : c                  write(*,*)'<q|q>',qb,qb2
     344             : 
     345             :                   nq=get_index_nn_q(bpt_q(:,q),nntot_q,qb,
     346             :      >                         gb_kq(:,nn,ikqpt),
     347             :      >                         gb_q(1:3,1:nntot_q,q),
     348           0 :      >                         arr_len,shift,l_dim)
     349             :                   nq2=get_index_nn_q(bpt_q(:,q),nntot_q,qb2,
     350             :      >                         gb_kq(:,nn2,ikqpt),
     351             :      >                         gb_q(1:3,1:nntot_q,q),
     352           0 :      >                         arr_len,shift,l_dim)
     353             : 
     354             : c                  CALL wann_gwf_auxovlp( param_vec(:,qb)+gb_q(:,nq,q),
     355             : c     >              param_vec(:,qb2)+gb_q(:,nq2,q),latt_const_q,ovaux)
     356             : 
     357           0 :                   do i=1,nbnd
     358           0 :                      do j=1,nbnd
     359           0 :                         write(307,'(2f24.18)') 0.0,0.0
     360             : c     >                      real(mq(j,i,nq)),aimag(mq(j,i,nq))
     361             :                      enddo
     362             :                   enddo
     363             : 
     364             :                else 
     365             : 
     366             :                   if((k.eq.kb).and.(k.ne.kb2).and.
     367           0 :      >               (q.ne.qb).and.(q.eq.qb2)) then
     368             : c                  write(*,*)'<q|k>',qb,qb2
     369             :                   nk2=get_index_nn_k(bpt(:,k),nntot,kb2,
     370             :      >                         gb_kq(:,nn2,ikqpt),
     371           0 :      >                         gb(1:3,1:nntot,k),arr_len)
     372             :                   nq=get_index_nn_q(bpt_q(:,q),nntot_q,qb,
     373             :      >                         gb_kq(:,nn,ikqpt),
     374             :      >                         gb_q(1:3,1:nntot_q,q),
     375           0 :      >                         arr_len,shift,l_dim)
     376             : 
     377             : c                  CALL wann_gwf_auxovlp( param_vec(:,qb)+gb_q(:,nq,q),
     378             : c     >              param_vec(:,qb2),latt_const_q,ovaux)
     379             : 
     380           0 :                   do i=1,nbnd
     381           0 :                      do j=1,nbnd
     382             :                         write(307,'(2f24.18)')
     383           0 :      >                    real(mkq(i,j,nq,nk2)),-aimag(mkq(i,j,nq,nk2))
     384             :                      enddo
     385             :                   enddo
     386             :                   
     387             :                   elseif((k.ne.kb).and.(k.eq.kb2).and.
     388           0 :      >                   (q.eq.qb).and.(q.ne.qb2)) then
     389             : c                  write(*,*)'<k|q>',qb,qb2
     390             :                   nk=get_index_nn_k(bpt(:,k),nntot,kb,
     391             :      >                         gb_kq(:,nn,ikqpt),
     392           0 :      >                         gb(1:3,1:nntot,k),arr_len)
     393             :                   nq2=get_index_nn_q(bpt_q(:,q),nntot_q,qb2,
     394             :      >                         gb_kq(:,nn2,ikqpt),
     395             :      >                         gb_q(1:3,1:nntot_q,q),
     396           0 :      >                         arr_len,shift,l_dim)
     397             : 
     398             : c                  CALL wann_gwf_auxovlp( param_vec(:,qb),
     399             : c     >              param_vec(:,qb2)+gb_q(:,nq2,q),latt_const_q,ovaux)
     400             : 
     401           0 :                   do i=1,nbnd
     402           0 :                      do j=1,nbnd
     403             :                         write(307,'(2f24.18)')
     404           0 :      >                    real(mkq(j,i,nq2,nk)),aimag(mkq(j,i,nq2,nk))
     405             :                      enddo
     406             :                   enddo
     407             : 
     408             :                   else
     409             :                      call juDFT_error("problem k,q neighbors",
     410           0 :      >                                calledby="wann_uHu_commat")
     411             :                   endif
     412             :                endif
     413             :                enddo
     414             : 
     415             :              enddo!nn
     416             :            enddo!k
     417           0 :            close(405)
     418           0 :            close(505)
     419             :          enddo!q           
     420           0 :          close(307)
     421           0 :          deallocate(mkk,mkq)
     422             :       endif!l_unformatted
     423             : 
     424           0 :       deallocate(gb_kq,bpt_kq)
     425           0 :       deallocate(fuHu,l_fm)
     426           0 :       deallocate(fuHu2,l_fm2)
     427             : 
     428           0 :       end subroutine wann_uHu_commat    
     429             :       end module m_wann_uHu_commat

Generated by: LCOV version 1.14