LCOV - code coverage report
Current view: top level - wannier - wann_rmat.f (source / functions) Coverage Total Hit
Test: FLEUR test coverage Lines: 0.0 % 174 0
Test Date: 2025-06-14 04:34:23 Functions: 0.0 % 1 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_rmat
       8              :       USE m_juDFT
       9              :       contains
      10            0 :       subroutine wann_rmat(
      11            0 :      >          bmat,amat, 
      12            0 :      >          rvecnum,rvec,kpoints,
      13              :      >          jspins_in,nkpts,l_bzsym,film,
      14              :      >          l_nocosoc,band_min,band_max,neigd,
      15              :      >          l_socmmn0,wan90version)
      16              : c***********************************************************
      17              : c     Calculate the matrix elements of the position operator
      18              : c     between wannier functions.
      19              : c
      20              : c     FF, 2009
      21              : c***********************************************************
      22              :       use m_constants
      23              :       use m_wann_read_umatrix
      24              : 
      25              :       implicit none
      26              : 
      27              :       real, intent(in) :: bmat(:,:)
      28              :       real, intent(in) :: amat(:,:)
      29              : 
      30              :       integer, intent(in) :: rvecnum
      31              :       integer, intent(in) :: rvec(:,:)
      32              :       real,    intent(in) :: kpoints(:,:)
      33              :       integer, intent(in) :: jspins_in
      34              :       integer, intent(in) :: nkpts
      35              :       logical, intent(in) :: l_bzsym
      36              :       logical, intent(in) :: film
      37              :    
      38              :       logical, intent(in) :: l_nocosoc
      39              :       integer, intent(in) :: band_min(2),band_max(2),neigd
      40              : 
      41              :       logical, intent(in) :: l_socmmn0
      42              :       integer, intent(in) :: wan90version
      43              : 
      44            0 :       real,allocatable    :: bpunkt(:,:)
      45              :       real                :: kdiff(3)
      46            0 :       real,allocatable    :: wb(:)
      47            0 :       integer,allocatable :: gb(:,:,:),bpt(:,:)
      48              :       integer             :: ikpt,jspins,ikpt_help
      49              :       integer             :: kpts,nntot,nn
      50              :       logical             :: l_file
      51              :       integer             :: num_wann,num_kpts,num_nnmax,jspin
      52              :       integer             :: kspin,kkspin
      53              :       integer             :: num_wann2
      54              :       integer             :: i,j,k,m,info,r1,r2,r3,dummy1
      55              :       integer             :: dummy2,dummy3,dummy4
      56              :       integer             :: hopmin,hopmax,counter,m1,m2
      57              :       integer,allocatable :: iwork(:)
      58            0 :       real,allocatable    :: energy(:,:),ei(:)
      59              :       real,allocatable    :: eigw(:,:),rwork(:)
      60              :       complex,allocatable :: work(:),vec(:,:)
      61            0 :       complex,allocatable :: u_matrix(:,:,:),m_matrix(:,:,:,:)
      62            0 :       complex,allocatable :: posop(:,:,:,:),posop2(:,:,:,:)
      63              :       complex             :: fac,eulav,eulav1
      64              :       real                :: tmp_omi,rdotk,tpi,minenerg,maxenerg
      65              :       real, allocatable   :: minieni(:),maxieni(:)
      66              :       character           :: jobz,uplo
      67              :       integer             :: kpt,band,lee,lwork,lrwork,liwork,n,lda
      68              :       real                :: rvalue
      69              :       logical             :: um_format
      70              :       logical             :: repro_eig
      71              :       logical             :: l_chk,l_proj
      72              :       logical             :: have_disentangled
      73            0 :       integer,allocatable :: ndimwin(:)
      74            0 :       logical,allocatable :: lwindow(:,:)
      75              :       integer             :: chk_unit,nkp,ntmp,ierr
      76              :       character(len=33)   :: header
      77              :       character(len=20)   :: checkpoint
      78              :       real                :: tmp_latt(3,3), tmp_kpt_latt(3,nkpts)
      79              :       real                :: omega_invariant
      80            0 :       complex,allocatable :: u_matrix_opt(:,:,:)
      81              :       integer             :: num_bands,nn2
      82              :       logical             :: l_umdat
      83            0 :       real,allocatable    :: eigval2(:,:)
      84            0 :       real,allocatable    :: eigval_opt(:,:)
      85              :       real                :: scale,a,b
      86              :       character(len=2)    :: spinspin12(0:2)
      87              :       character(len=3)    :: spin12(2)
      88              :       character(len=6)    :: filename
      89              :       integer             :: jp,mp,kk
      90              :       integer             :: rvecind,rvecind_0
      91              :       real                :: realvec(3)
      92              :       complex,allocatable :: mmnk(:,:,:,:)
      93              :       logical             :: l_worksout,l_she
      94              : 
      95              :       data spinspin12/'  ','.1' , '.2'/
      96              :       data spin12/'WF1','WF2'/
      97              : 
      98            0 :       call timestart("wann_rmat")
      99              : 
     100            0 :       tpi=2*pimach()
     101              : 
     102            0 :       jspins=jspins_in
     103            0 :       if(l_nocosoc)jspins=1
     104              : 
     105            0 :       write(oUnit,*)"nkpts=",nkpts
     106              : 
     107            0 :       do jspin=1,jspins  !spin loop
     108              : c*****************************************************
     109              : c        get num_bands and num_wann from the proj file
     110              : c*****************************************************
     111            0 :          do j=jspin,0,-1
     112            0 :           inquire(file=trim('proj'//spinspin12(j)),exist=l_file)
     113            0 :           if(l_file)then
     114            0 :             filename='proj'//spinspin12(j)
     115            0 :             exit
     116              :           endif
     117              :          enddo
     118            0 :          if(l_file)then
     119            0 :           open (203,file=trim(filename),status='old')
     120            0 :           rewind (203)
     121              :          else
     122              :             CALL judft_error("no proj/proj.1/proj.2",calledby
     123            0 :      +           ="wann_rmat")
     124              :          endif
     125            0 :          read (203,*) num_wann,num_bands
     126            0 :          close (203)
     127            0 :          write(oUnit,*)'According to proj there are ',
     128            0 :      +                 num_bands,' bands'
     129            0 :          write(oUnit,*)"and ",num_wann," wannier functions."
     130              : 
     131              : c****************************************************************
     132              : c        get nntot and bk and wb from bkpts file
     133              : c****************************************************************
     134            0 :          inquire (file='bkpts',exist=l_file)
     135            0 :          if (.not.l_file)  CALL judft_error("need bkpts"
     136            0 :      +      ,calledby ="wann_rmat")
     137            0 :          open (202,file='bkpts',form='formatted',status='old')
     138            0 :          read (202,'(i4)') nntot
     139            0 :          allocate ( gb(3,nntot,nkpts) )
     140            0 :          allocate ( bpt(nntot,nkpts) )
     141            0 :          do ikpt=1,nkpts
     142            0 :           do nn=1,nntot
     143              :            read (202,'(2i6,3x,3i4)')
     144            0 :      &      ikpt_help,bpt(nn,ikpt),(gb(i,nn,ikpt),i=1,3)
     145              :           enddo !nn
     146              :          enddo !ikpt
     147            0 :          allocate( wb(nntot) )
     148            0 :          allocate( bpunkt(3,nntot) )
     149            0 :          do nn=1,nntot
     150            0 :             read(202,*)bpunkt(:,nn),wb(nn)
     151              :          enddo
     152            0 :          close (202)
     153              : 
     154            0 :          do ikpt=1,nkpts
     155            0 :           do nn=1,nntot
     156              :            kdiff(1)=kpoints(1,bpt(nn,ikpt))+
     157            0 :      +            gb(1,nn,ikpt) - kpoints(1,ikpt)
     158              :            kdiff(2)=kpoints(2,bpt(nn,ikpt))+
     159            0 :      +            gb(2,nn,ikpt) - kpoints(2,ikpt)
     160              :            kdiff(3)=kpoints(3,bpt(nn,ikpt))+
     161            0 :      +            gb(3,nn,ikpt) - kpoints(3,ikpt)
     162            0 :            kdiff=matmul(transpose(bmat),kdiff)
     163            0 :            l_worksout=.false.
     164            0 :            do nn2=1,nntot
     165            0 :              if (abs(kdiff(1)-bpunkt(1,nn2)).lt.1.e-5)then
     166            0 :                if (abs(kdiff(2)-bpunkt(2,nn2)).lt.1.e-5)then
     167            0 :                  if (abs(kdiff(3)-bpunkt(3,nn2)).lt.1.e-5)then
     168              :                      l_worksout=.true. 
     169              :                      exit
     170              :                  endif          
     171              :                endif
     172              :              endif 
     173              :            enddo !nn2
     174            0 :            if(.not.l_worksout)then
     175            0 :               write(*,*)"ikpt,nn=",ikpt,nn
     176            0 :               write(*,*)"kdiff(1)=",kdiff(1)
     177            0 :               write(*,*)"kdiff(2)=",kdiff(2)
     178            0 :               write(*,*)"kdiff(3)=",kdiff(3)
     179            0 :               stop 'worksout'
     180              :            endif   
     181              :           enddo !nn
     182              :          enddo !ikpt 
     183              : 
     184              : c****************************************************************
     185              : c        read in chk
     186              : c****************************************************************
     187            0 :          num_kpts=nkpts
     188            0 :          allocate( u_matrix_opt(num_bands,num_wann,nkpts) )
     189            0 :          allocate( u_matrix(num_wann,num_wann,nkpts) )
     190            0 :          allocate( lwindow(num_bands,nkpts) )
     191            0 :          allocate( ndimwin(nkpts) )
     192            0 :          allocate( m_matrix(num_wann,num_wann,nntot,num_kpts) )
     193              :          call wann_read_umatrix2(
     194              :      >            nkpts,num_wann,num_bands,
     195              :      >            um_format,jspin,wan90version,
     196              :      <            have_disentangled,
     197              :      <            lwindow,ndimwin,
     198            0 :      <            u_matrix_opt,u_matrix,m_matrix)
     199              : c****************************************************************
     200              : c        read in eig-file
     201              : c****************************************************************
     202            0 :          write(oUnit,*)"read in eig-file"
     203            0 :          allocate(energy(num_bands,num_kpts))
     204            0 :          inquire(file=spin12(jspin)//'.eig',exist=l_umdat)
     205            0 :          IF(.NOT.l_umdat)  CALL judft_error
     206              :      +        ("Thou shall not hide your eig file",calledby
     207            0 :      +        ="wann_hopping")
     208            0 :          open(300,file=spin12(jspin)//'.eig',form='formatted')
     209            0 :          do i=1,num_kpts
     210            0 :            do j=1,num_bands
     211            0 :               read(300,*)band,kpt,energy(j,i)
     212              :            enddo
     213              :          enddo
     214            0 :          close(300)
     215              : 
     216            0 :          minenerg=minval(energy(:,:))
     217            0 :          maxenerg=maxval(energy(:,:))
     218            0 :          write(oUnit,*)"minenerg=",minenerg
     219            0 :          write(oUnit,*)"maxenerg=",maxenerg
     220              : 
     221              : 
     222            0 :          allocate(eigval_opt(num_bands,nkpts))
     223            0 :          allocate(eigval2(num_wann,nkpts))
     224            0 :          eigval_opt=0.0
     225            0 :          eigval2=0.0
     226              : 
     227            0 :          if(have_disentangled) then
     228              : 
     229            0 :            do nkp=1,num_kpts
     230            0 :             counter=0
     231            0 :             do j=1,num_bands
     232            0 :               if(lwindow(j,nkp)) then
     233            0 :                 counter=counter+1
     234            0 :                 eigval_opt(counter,nkp)=energy(j,nkp)
     235              :               end if
     236              :             end do
     237              :            end do
     238              :        
     239            0 :            do nkp=1,num_kpts
     240            0 :             do j=1,num_wann
     241            0 :              do m=1,ndimwin(nkp)
     242              :                 eigval2(j,nkp)=eigval2(j,nkp)+eigval_opt(m,nkp)* 
     243            0 :      &    real(conjg(u_matrix_opt(m,j,nkp))*u_matrix_opt(m,j,nkp))
     244              :              enddo
     245              :             enddo
     246              :            enddo
     247              : 
     248              :          else
     249            0 :            eigval2 = energy
     250              :          end if                    !have_disentangled
     251              : 
     252            0 :          deallocate(eigval_opt)
     253            0 :          deallocate(energy)
     254              : 
     255              : c****************************************************************
     256              : c        Set up posop.
     257              : c****************************************************************
     258            0 :          write(oUnit,*)"Set up posop."
     259            0 :          allocate( posop2(3,num_wann,num_wann,num_kpts) )
     260            0 :          posop2=cmplx(0.0,0.0)
     261              :          
     262            0 :          do ikpt=1,nkpts
     263            0 :           do nn=1,nntot
     264              :            kdiff(1)=kpoints(1,bpt(nn,ikpt))+
     265            0 :      +            gb(1,nn,ikpt) - kpoints(1,ikpt)
     266              :            kdiff(2)=kpoints(2,bpt(nn,ikpt))+
     267            0 :      +            gb(2,nn,ikpt) - kpoints(2,ikpt)
     268              :            kdiff(3)=kpoints(3,bpt(nn,ikpt))+
     269            0 :      +            gb(3,nn,ikpt) - kpoints(3,ikpt)
     270            0 :            kdiff=matmul(transpose(bmat),kdiff)          
     271            0 :            do i=1,num_wann
     272            0 :             do j=1,num_wann
     273            0 :              if(j.eq.i)then
     274            0 :               do kk=1,3
     275              :                posop2(kk,j,i,ikpt)=posop2(kk,j,i,ikpt)+ImagUnit*
     276            0 :      &             wb(nn)*kdiff(kk)*(m_matrix(j,i,nn,ikpt)-1.0)  
     277              :               enddo !kk
     278              :              else
     279            0 :               do kk=1,3
     280              :                posop2(kk,j,i,ikpt)=posop2(kk,j,i,ikpt)+ImagUnit*
     281            0 :      &             wb(nn)*kdiff(kk)*m_matrix(j,i,nn,ikpt)  
     282              :               enddo !kk
     283              :              endif
     284              :             enddo !j
     285              :            enddo !i
     286              :           enddo !nn   
     287              :          enddo !ikpt
     288              : 
     289            0 :          do ikpt=1,nkpts
     290            0 :           do i=1,num_wann
     291            0 :            do j=1,i
     292            0 :             do kk=1,3
     293              :               posop2(kk,j,i,ikpt)=(posop2(kk,j,i,ikpt)+
     294            0 :      &                        conjg(posop2(kk,i,j,ikpt)))/2.0
     295            0 :               posop2(kk,i,j,ikpt)=posop2(kk,j,i,ikpt)
     296              :             enddo !kk
     297              :            enddo !j
     298              :           enddo !i
     299              :          enddo !ikpt
     300              : 
     301            0 :          allocate( posop(3,num_wann,num_wann,rvecnum) )
     302            0 :          posop=cmplx(0.0,0.0)
     303            0 :          do rvecind=1,rvecnum
     304            0 :             do k=1,nkpts  
     305              :               rdotk=tpi*( kpoints(1,k)*rvec(1,rvecind)+
     306              :      &                    kpoints(2,k)*rvec(2,rvecind)+
     307            0 :      &                    kpoints(3,k)*rvec(3,rvecind) )
     308            0 :               fac=cmplx(cos(rdotk),-sin(rdotk))
     309            0 :               do i=1,num_wann
     310            0 :                do j=1,num_wann
     311            0 :                 do kk=1,3
     312              :                  posop(kk,j,i,rvecind)=
     313              :      &                   posop(kk,j,i,rvecind)+
     314            0 :      &                   fac*posop2(kk,j,i,k)
     315              :                 enddo !kk
     316              :                enddo !j
     317              :               enddo !i
     318              :             enddo !k
     319              :          enddo !rvecind
     320            0 :          posop=posop/cmplx(real(nkpts),0.0)
     321            0 :          deallocate( posop2 )
     322              : 
     323            0 :          do rvecind=1,rvecnum
     324            0 :           if(rvec(1,rvecind).eq.0)then
     325            0 :            if(rvec(2,rvecind).eq.0)then
     326            0 :             if(rvec(3,rvecind).eq.0)then
     327            0 :                rvecind_0=rvecind
     328              :                goto 123
     329              :             endif
     330              :            endif
     331              :           endif
     332              :          enddo !rvecind
     333            0 :          stop 'Ou est ce point-la ?'
     334              :  123     continue
     335              : 
     336            0 :          do i=1,num_wann
     337            0 :           do kk=1,3 
     338            0 :              posop(kk,i,i,rvecind_0)=0.0
     339              :           enddo !kk
     340              :          enddo !i
     341              : 
     342            0 :          do ikpt=1,nkpts
     343            0 :           do nn=1,nntot  
     344              :            kdiff(1)=kpoints(1,bpt(nn,ikpt))+
     345            0 :      +            gb(1,nn,ikpt) - kpoints(1,ikpt)
     346              :            kdiff(2)=kpoints(2,bpt(nn,ikpt))+
     347            0 :      +            gb(2,nn,ikpt) - kpoints(2,ikpt)
     348              :            kdiff(3)=kpoints(3,bpt(nn,ikpt))+
     349            0 :      +            gb(3,nn,ikpt) - kpoints(3,ikpt)
     350            0 :            kdiff=matmul(transpose(bmat),kdiff)/real(nkpts)             
     351            0 :            do i=1,num_wann
     352            0 :             do kk=1,3 
     353              :              posop(kk,i,i,rvecind_0)=posop(kk,i,i,rvecind_0)-
     354            0 :      &        wb(nn)*kdiff(kk)*aimag(log(m_matrix(i,i,nn,ikpt)))
     355              :             enddo !kk
     356              :            enddo !i
     357              :           enddo !nn
     358              :          enddo !ikpt
     359              : 
     360              : c********************************************
     361              : c        Print posop.
     362              : c********************************************
     363            0 :          open(321,file='posop'//spinspin12(jspin),form='formatted')
     364            0 :          do rvecind=1,rvecnum
     365            0 :             r3=rvec(3,rvecind)
     366            0 :             r2=rvec(2,rvecind)
     367            0 :             r1=rvec(1,rvecind)             
     368            0 :             do j=1,num_wann
     369            0 :              do i=1,num_wann
     370            0 :               do kk=1,3  
     371              :                write(321,'(i3,1x,i3,1x,i3,1x,i3,
     372              :      &                 1x,i3,1x,i3,1x,f20.8,1x,f20.8)')
     373            0 :      &                 r1,r2,r3,i,j,kk,posop(kk,i,j,rvecind) 
     374              :               enddo !kk  
     375              :              enddo !j
     376              :             enddo !i
     377              :          enddo !rvecind  
     378            0 :          close(321)
     379              : 
     380            0 :          deallocate( lwindow,u_matrix_opt,ndimwin )
     381            0 :          deallocate( u_matrix,m_matrix )
     382            0 :          deallocate( posop,eigval2 )
     383            0 :          deallocate( gb,bpt,wb,bpunkt )
     384              :       enddo !jspin
     385              : 
     386            0 :       call timestop("wann_rmat")
     387            0 :       end subroutine wann_rmat
     388            0 :       end module m_wann_rmat
        

Generated by: LCOV version 2.0-1