LCOV - code coverage report
Current view: top level - wannier - wann_perpmag_rs.f (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 198 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_perpmag_rs
       8             :       use m_juDFT
       9             :       contains 
      10           0 :       subroutine wann_perpmag_rs(
      11           0 :      >               rvecnum,rvec,kpoints,
      12             :      >               jspins_in,nkpts,l_bzsym,film,l_onedimens,
      13             :      >               l_soc,band_min,band_max,neigd,
      14           0 :      >               l_socmmn0,l_ndegen,ndegen,wan90version,
      15             :      >               l_unformatted)
      16             : c*************************************************
      17             : c     Calculate the matrix elements of  
      18             : c     <n|B_{eff}|m> in real space from the
      19             : c     files WF1.chk (and WF1_um.dat) produced
      20             : c     by wannier90.
      21             : c
      22             : c     Frank Freimuth
      23             : c*************************************************
      24             :       use m_constants, only:pimach, ImagUnit
      25             :       use m_wann_read_umatrix
      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             :       logical, intent(in) :: l_ndegen
      38             :       integer, intent(in) :: ndegen(:)
      39             :       integer, intent(in) :: wan90version
      40             :       logical, intent(in) :: l_unformatted      
      41             : 
      42             :       integer             :: ikpt,jspins
      43             :       integer             :: kpts
      44             :       logical             :: l_file
      45             :       integer             :: num_wann,num_kpts,num_nnmax,jspin
      46             :       integer             :: kspin,kkspin
      47             :       integer             :: wann_shift,num_wann2
      48             :       integer             :: i,j,k,m,info,r1,r2,r3,dummy1
      49             :       integer             :: dummy2,dummy3
      50             :       integer             :: hopmin,hopmax,counter,m1,m2
      51             :       integer             :: num_bands2
      52             :       integer,allocatable :: iwork(:)
      53             :       real,allocatable    :: energy(:,:),ei(:)
      54             :       real,allocatable    :: eigw(:,:),rwork(:)
      55             :       complex,allocatable :: work(:),vec(:,:)
      56           0 :       complex,allocatable :: u_matrix(:,:,:,:),hwann(:,:,:,:)
      57           0 :       complex,allocatable :: hwannmix(:,:,:,:)
      58           0 :       complex,allocatable :: hreal(:,:,:,:)
      59           0 :       complex,allocatable :: paulimat(:,:,:,:)
      60           0 :       complex,allocatable :: paulimatmix(:,:,:,:)
      61           0 :       complex,allocatable :: paulimat_opt(:,:,:,:)
      62           0 :       complex,allocatable :: paulimat2(:,:,:,:)
      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             :       complex             :: value(4)
      69             :       logical             :: um_format
      70             :       logical             :: repro_eig
      71             :       logical             :: l_chk,l_proj,l_onedimens
      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,counter1,counter2
      82             :       logical             :: l_umdat
      83             :       real,allocatable    :: eigval2(:,:)
      84             :       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             :: hopmin_x,hopmin_y,hopmin_z
      91             :       integer             :: hopmax_x,hopmax_y,hopmax_z
      92             :       integer             :: rvecind
      93             :       integer             :: ii,int_dummy
      94             :       integer             :: dir,num(3)
      95             :       integer             :: spin1,spin2
      96             :       logical             :: l_addint
      97             :       integer             :: nbnd,nwfs,fullnkpts
      98           0 :       complex,allocatable :: amn(:,:,:)
      99           0 :       complex,allocatable :: perpmag(:,:,:,:)   
     100             : 
     101             :       data spinspin12/'  ','.1' , '.2'/
     102             :       data spin12/'WF1','WF2'/
     103             : 
     104           0 :       tpi=2*pimach()
     105             : 
     106           0 :       jspins=jspins_in
     107           0 :       if(l_soc)jspins=1
     108             : 
     109           0 :       write(6,*)"nkpts=",nkpts
     110             : 
     111             : 
     112             : c*****************************************************
     113             : c     get num_bands and num_wann from the proj file
     114             : c*****************************************************
     115           0 :       do j=jspins,0,-1
     116           0 :           inquire(file=trim('proj'//spinspin12(j)),exist=l_file)
     117           0 :           if(l_file)then
     118           0 :             filename='proj'//spinspin12(j)
     119           0 :             exit
     120             :           endif
     121             :       enddo
     122           0 :       if(l_file)then
     123           0 :           open (203,file=trim(filename),status='old')
     124           0 :           rewind (203)
     125             :       else
     126             :             CALL juDFT_error("no proj/proj.1/proj.2",calledby
     127           0 :      +           ="wann_perpmag_rs")
     128             :       endif
     129           0 :       read (203,*) num_wann,num_bands
     130           0 :       close (203)
     131           0 :       write(6,*)'According to proj there are ',num_bands,' bands'
     132           0 :       write(6,*)"and ",num_wann," wannier functions."
     133             : 
     134           0 :       allocate( u_matrix_opt(num_bands,num_wann,nkpts,jspins) )
     135           0 :       allocate( u_matrix(num_wann,num_wann,nkpts,jspins) )
     136           0 :       allocate( lwindow(num_bands,nkpts,jspins) )
     137           0 :       allocate( ndimwin(nkpts,jspins) )
     138             : 
     139           0 :       do jspin=1,jspins  !spin loop
     140             : c****************************************************************
     141             : c        read in chk
     142             : c****************************************************************
     143           0 :          num_kpts=nkpts
     144             : 
     145             :          call wann_read_umatrix2(
     146             :      >            nkpts,num_wann,num_bands,
     147             :      >            um_format,jspin,wan90version,
     148             :      <            have_disentangled,
     149             :      <            lwindow(:,:,jspin),
     150             :      <            ndimwin(:,jspin),
     151             :      <            u_matrix_opt(:,:,:,jspin),
     152           0 :      <            u_matrix(:,:,:,jspin) )
     153             : 
     154             : 
     155             : 
     156             :       enddo !jspin   
     157             :  
     158           0 :       jspin=1
     159             : 
     160           0 :       num_bands2=jspins*num_bands
     161           0 :       num_wann2=jspins*num_wann
     162             : 
     163           0 :       allocate( paulimat(3,num_bands2,num_bands2,nkpts) )
     164           0 :       paulimat=0.0
     165             : 
     166             : c****************************************************
     167             : c        Read the file "updown.perpmag".
     168             : c****************************************************
     169           0 :       if(l_unformatted)then
     170           0 :        inquire(file='updown.perpmagint_unf',exist=l_addint)
     171           0 :        open(304,file='updown.perpmag_unf',form='unformatted')
     172           0 :         if(l_addint)then
     173           0 :          open(307,file='updown.perpmagint_unf',form='unformatted')
     174             :         endif
     175           0 :         read(304)nbnd,fullnkpts,nwfs
     176           0 :         allocate(amn(nbnd,nwfs,fullnkpts))
     177           0 :         read(304)amn(1:nbnd,1:nwfs,1:fullnkpts)
     178           0 :         do nkp=1,num_kpts
     179           0 :          do i=1,num_bands
     180           0 :           do j=1,num_bands
     181           0 :            paulimat(1,j,i+num_bands*(jspins-1),nkp)=amn(j,i,nkp)
     182             :           enddo !j
     183             :          enddo !i
     184             :         enddo !nkp
     185           0 :         if(l_addint)then
     186           0 :           read(307)nbnd,fullnkpts,nwfs
     187           0 :           read(307)amn(1:nbnd,1:nwfs,1:fullnkpts)
     188           0 :           do nkp=1,num_kpts
     189           0 :            do i=1,num_bands
     190           0 :             do j=1,num_bands
     191             :              paulimat(1,j,i+num_bands*(jspins-1),nkp)=
     192           0 :      &       paulimat(1,j,i+num_bands*(jspins-1),nkp)+amn(j,i,nkp)
     193             :             enddo !j
     194             :            enddo !i
     195             :           enddo !nkp
     196             :         endif
     197           0 :         deallocate(amn)
     198             :       else              
     199           0 :        inquire(file='updown.perpmagint',exist=l_addint)
     200           0 :        open(304,file='updown.perpmag',form='formatted')
     201           0 :        read(304,*)
     202           0 :        read(304,*)
     203           0 :        if(l_addint)then
     204           0 :            open(307,file='updown.perpmagint',form='formatted')
     205           0 :            read(307,*)
     206           0 :            read(307,*)
     207             :        endif
     208           0 :        do nkp=1,num_kpts
     209           0 :          do i=1,num_bands
     210           0 :            do j=1,num_bands
     211           0 :              read(304,*)dummy1,dummy2,dummy3,a,b
     212           0 :              paulimat(1,j,i+num_bands*(jspins-1),nkp)=cmplx(a,b)
     213             :            enddo !j
     214             :          enddo !i
     215           0 :          if(l_addint)then
     216           0 :           do i=1,num_bands
     217           0 :            do j=1,num_bands
     218           0 :              read(307,*)dummy1,dummy2,dummy3,a,b
     219             :              paulimat(1,j,i+num_bands*(jspins-1),nkp)=
     220             :      =       paulimat(1,j,i+num_bands*(jspins-1),nkp)+
     221           0 :      +                 cmplx(a,b)
     222             :            enddo !j
     223             :           enddo !i
     224             :          endif
     225             :        enddo !nkp
     226             :       endif
     227           0 :       close(304)
     228           0 :       if(l_addint)then
     229           0 :         close(307)
     230             :       endif
     231             : 
     232           0 :       do nkp=1,num_kpts
     233           0 :           paulimat(2,:,:,nkp)=ImagUnit*paulimat(1,:,:,nkp)
     234             :           paulimat(1,:,:,nkp)=paulimat(1,:,:,nkp)+
     235           0 :      &        transpose(conjg( paulimat(1,:,:,nkp) ))
     236             :           paulimat(2,:,:,nkp)=paulimat(2,:,:,nkp)+
     237           0 :      &        transpose(conjg( paulimat(2,:,:,nkp) ))
     238             :       enddo     
     239             : 
     240             : c****************************************************************
     241             : c        Calculate matrix elements of B_{eff} in the basis of
     242             : c        rotated Bloch functions.
     243             : c****************************************************************
     244           0 :       allocate( paulimat2(3,num_wann2,num_wann2,nkpts) )
     245             :       write(6,*)"calculate matrix elements of exchange field
     246           0 :      & between wannier orbitals"
     247             : 
     248           0 :       if(have_disentangled) then       
     249             : 
     250           0 :           allocate( paulimat_opt(3,num_bands2,num_bands2,nkpts) )
     251           0 :           allocate( paulimatmix(3,num_wann2,num_bands2,nkpts))
     252             : 
     253           0 :           do nkp=1,num_kpts
     254             :            counter1=0
     255           0 :            do m=1,num_bands2
     256           0 :             spin1=(m-1)/num_bands 
     257           0 :             if(lwindow(m-spin1*num_bands,nkp,spin1+1))then
     258           0 :              counter1=counter1+1  
     259           0 :              counter2=0
     260           0 :              do mp=1,num_bands2
     261           0 :               spin2=(mp-1)/num_bands 
     262           0 :               if(lwindow(mp-spin2*num_bands,nkp,spin2+1))then
     263           0 :                counter2=counter2+1
     264           0 :                do k=1,3
     265             :                 paulimat_opt(k,counter2,counter1,nkp)=
     266           0 :      &          paulimat(k,mp,m,nkp)  
     267             :                enddo
     268             :               endif
     269             :              enddo !mp
     270             :             endif
     271             :            enddo !m 
     272             :           enddo
     273             : 
     274           0 :           paulimatmix=0.0  
     275           0 :           do spin1=0,jspins-1
     276           0 :             do spin2=0,jspins-1
     277           0 :              do nkp=1,num_kpts
     278           0 :               do jp=1,num_wann  
     279           0 :                do m=1,ndimwin(nkp,spin2+1)
     280           0 :                 do mp=1,ndimwin(nkp,spin1+1)
     281           0 :                  do k=1,3  
     282             : 
     283             : 
     284             :           paulimatmix(k,jp+spin1*num_wann,m+spin2*ndimwin(nkp,1),nkp)=
     285             :      =    paulimatmix(k,jp+spin1*num_wann,m+spin2*ndimwin(nkp,1),nkp)+
     286             :      &                 conjg(u_matrix_opt(mp,jp,nkp,spin1+1))*
     287             :      &                        paulimat_opt(k,mp+spin1*ndimwin(nkp,1),
     288           0 :      &                                   m+spin2*ndimwin(nkp,1),nkp)
     289             :                  enddo !k 
     290             :                 enddo !mp   
     291             :                enddo !m
     292             :               enddo !jp 
     293             :              enddo !nkp
     294             :             enddo !spin2
     295             :           enddo !spin1
     296             : 
     297           0 :           paulimat2=0.0  
     298           0 :           do spin1=0,jspins-1
     299           0 :            do spin2=0,jspins-1
     300           0 :             do nkp=1,num_kpts
     301           0 :              do j=1,num_wann
     302           0 :               do jp=1,num_wann  
     303           0 :                do m=1,ndimwin(nkp,spin2+1)
     304             : 
     305           0 :                  do k=1,3  
     306             : 
     307             : 
     308             :                 paulimat2(k,jp+spin1*num_wann,j+spin2*num_wann,nkp)=
     309             :      =          paulimat2(k,jp+spin1*num_wann,j+spin2*num_wann,nkp)+ 
     310             :      &                        paulimatmix(k,jp+spin1*num_wann,
     311             : 
     312             :      &                                   m+spin2*ndimwin(nkp,1),nkp)*
     313           0 :      &                       u_matrix_opt(m,j,nkp,spin2+1)
     314             :                  enddo !k 
     315             :                enddo !m
     316             :               enddo !jp 
     317             :              enddo !j
     318             :             enddo !nkp
     319             :            enddo !spin2
     320             :           enddo !spin1   
     321             : 
     322             : 
     323           0 :          deallocate(paulimat_opt)
     324           0 :         deallocate(paulimatmix)
     325             : 
     326             :       else
     327           0 :           paulimat2=paulimat
     328             :       end if !have_disentangled
     329             : 
     330           0 :       allocate(hwann(3,num_wann2,num_wann2,num_kpts))
     331           0 :       hwann=cmplx(0.0,0.0)
     332           0 :       wann_shift=0
     333             : 
     334             : 
     335           0 :       allocate(hwannmix(3,num_wann2,num_wann2,num_kpts))
     336           0 :       hwannmix=cmplx(0.0,0.0)
     337             : 
     338           0 :       do spin1=0,jspins-1
     339           0 :        do spin2=0,jspins-1
     340           0 :         do k=1,num_kpts
     341           0 :          do i=1,num_wann
     342           0 :           do mp=1,num_wann
     343           0 :            do j=1,num_wann
     344             : 
     345           0 :              do kk=1,3  
     346             :               hwannmix(kk,mp+spin1*num_wann,i+spin2*num_wann,k)=
     347             :      =         hwannmix(kk,mp+spin1*num_wann,i+spin2*num_wann,k)+
     348             :      *           conjg(u_matrix(j,mp,k,spin1+1))*
     349             :      *               paulimat2(kk,j+spin1*num_wann,
     350           0 :      *                            i+spin2*num_wann,k)
     351             : 
     352             :              enddo !kk
     353             :            enddo !j
     354             :           enddo !mp
     355             : 
     356             :          enddo !i     
     357             :         enddo !k
     358             :        enddo !spin2
     359             :       enddo !spin1 
     360             : 
     361           0 :       do spin1=0,jspins-1
     362           0 :        do spin2=0,jspins-1
     363           0 :         do k=1,num_kpts
     364           0 :          do m=1,num_wann
     365           0 :           do i=1,num_wann
     366           0 :            do mp=1,num_wann
     367           0 :             do kk=1,3  
     368             :               hwann(kk,mp+spin1*num_wann,m+spin2*num_wann,k)=
     369             :      =         hwann(kk,mp+spin1*num_wann,m+spin2*num_wann,k)+
     370             :      *               hwannmix(kk,mp+spin1*num_wann,
     371             :      *                            i+spin2*num_wann,k)*
     372           0 :      *                 u_matrix(i,m,k,spin2+1)
     373             :             enddo !kk
     374             :            enddo !mp
     375             :           enddo !i     
     376             :          enddo !m
     377             :         enddo !k
     378             :        enddo !spin2
     379             :       enddo !spin1
     380             : 
     381           0 :       deallocate(hwannmix)
     382             : 
     383             : 
     384             : 
     385             : c************************************************************
     386             : c        Calculate matrix elements in real space.
     387             : c***********************************************************      
     388           0 :        write(6,*)"calculate pauli-mat in rs"
     389             : 
     390             : 
     391           0 :        allocate( hreal(3,num_wann2,num_wann2,rvecnum) )
     392           0 :        hreal=cmplx(0.0,0.0)
     393           0 :        if(l_ndegen)then
     394           0 :         do rvecind=1,rvecnum
     395           0 :           do k=1,nkpts  
     396             :            rdotk=tpi*( kpoints(1,k)*rvec(1,rvecind)+
     397             :      &                 kpoints(2,k)*rvec(2,rvecind)+
     398           0 :      &                 kpoints(3,k)*rvec(3,rvecind) )
     399           0 :            fac=cmplx(cos(rdotk),-sin(rdotk))/real(ndegen(rvecind))
     400           0 :            do m2=1,num_wann2
     401           0 :             do m1=1,num_wann2
     402           0 :              do dir=1,3  
     403             :               hreal(dir,m1,m2,rvecind)=hreal(dir,m1,m2,rvecind)+
     404           0 :      &                   fac*hwann(dir,m1,m2,k)
     405             :              enddo !dir 
     406             :             enddo !m1  
     407             :            enddo !m2  
     408             :           enddo !k
     409             :         enddo !rvecind
     410             :        else 
     411           0 :         do rvecind=1,rvecnum
     412           0 :           do k=1,nkpts  
     413             :            rdotk=tpi*( kpoints(1,k)*rvec(1,rvecind)+
     414             :      &                 kpoints(2,k)*rvec(2,rvecind)+
     415           0 :      &                 kpoints(3,k)*rvec(3,rvecind) )
     416           0 :            fac=cmplx(cos(rdotk),-sin(rdotk))
     417           0 :            do m2=1,num_wann2
     418           0 :             do m1=1,num_wann2
     419           0 :              do dir=1,3  
     420             :               hreal(dir,m1,m2,rvecind)=hreal(dir,m1,m2,rvecind)+
     421           0 :      &                   fac*hwann(dir,m1,m2,k)
     422             :              enddo !dir 
     423             :             enddo !m1  
     424             :            enddo !m2  
     425             :           enddo !k
     426             :         enddo !rvecind
     427             :        endif   
     428           0 :        hreal=hreal/cmplx(real(nkpts),0.0)
     429             : 
     430             : c       open(321,file='rsperpmag'//spinspin12(1),form='formatted')
     431             : c       do rvecind=1,rvecnum
     432             : c         r3=rvec(3,rvecind)
     433             : c         r2=rvec(2,rvecind)
     434             : c         r1=rvec(1,rvecind)
     435             : c         do j=1,num_wann2
     436             : c          do i=1,num_wann2
     437             : c           do kk=1,3   
     438             : c            write(321,'(i3,1x,i3,1x,i3,1x,i3,
     439             : c     &           1x,i3,1x,i3,1x,f20.8,1x,f20.8)')
     440             : c     &          r1,r2,r3,i,j,kk,hreal(kk,i,j,rvecind) 
     441             : c           enddo !kk 
     442             : c          enddo !i
     443             : c         enddo !j   
     444             : c       enddo !rvecnum 
     445             : c       close(321)
     446           0 :        if(l_unformatted)then
     447           0 :         allocate(perpmag(num_wann2,num_wann2,3,rvecnum))
     448           0 :         if(l_ndegen)then
     449             :            open(321,file='tunfndegen',form='unformatted',
     450           0 :      &                     convert='BIG_ENDIAN')
     451             :         else
     452             :            open(321,file='tunf',form='unformatted',
     453           0 :      &                     convert='BIG_ENDIAN')
     454             :         endif
     455           0 :         do rvecind=1,rvecnum
     456           0 :          do j=1,num_wann2
     457           0 :           do i=1,num_wann2
     458           0 :            do kk=1,3
     459           0 :              perpmag(i,j,kk,rvecind)=hreal(kk,i,j,rvecind)
     460             :            enddo !kk
     461             :           enddo !i
     462             :          enddo !j
     463             :         enddo !rvecind
     464           0 :         write(321)perpmag
     465           0 :         deallocate(perpmag)
     466             :        else             
     467           0 :         if(l_ndegen)then
     468             : 
     469             :          open(321,file='rsperpmag_ndegen'//spinspin12(1),
     470             :      &           form='formatted',
     471           0 :      &        recl=1000)
     472             :         else
     473             :          open(321,file='rsperpmag'//spinspin12(1),form='formatted',
     474           0 :      &        recl=1000)
     475             :         endif
     476             : 
     477           0 :         do rvecind=1,rvecnum
     478           0 :          r3=rvec(3,rvecind)
     479           0 :          r2=rvec(2,rvecind)
     480           0 :          r1=rvec(1,rvecind)
     481           0 :          do j=1,num_wann2
     482           0 :           do i=1,num_wann2
     483           0 :            do kk=1,3
     484           0 :             write(321,*) r1,r2,r3,i,j,kk,real(hreal(kk,i,j,rvecind)),
     485           0 :      &                                   aimag(hreal(kk,i,j,rvecind))
     486             :            enddo !kk 
     487             :           enddo !i
     488             :          enddo !j   
     489             :         enddo !rvecnum
     490             :        endif !l_unformatted      
     491           0 :        close(321)
     492             : 
     493             : 
     494           0 :        deallocate(lwindow,u_matrix_opt,ndimwin)
     495           0 :        deallocate(u_matrix,hwann,hreal)
     496             : 
     497           0 :       end subroutine wann_perpmag_rs
     498             :       end module m_wann_perpmag_rs

Generated by: LCOV version 1.13