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

Generated by: LCOV version 2.0-1