LCOV - code coverage report
Current view: top level - wannier - wann_mmnk_symm.f (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 169 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 2 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_mmnk_symm
       8             :       use m_juDFT
       9             :       USE m_types
      10             :       private
      11             :       public:: wann_mmnk_symm
      12             :       contains
      13             : c******************************************************************
      14             : c     Find out minimal set of k-point-pairs that have to be
      15             : c     calculated; map symmetry-related k-point-pairs to this
      16             : c     minimal set.
      17             : c     Frank Freimuth
      18             : c******************************************************************      
      19           0 :       subroutine wann_mmnk_symm(input,kpts,
      20           0 :      >               fullnkpts,nntot,bpt,gb,l_bzsym,
      21           0 :      >               irreduc,mapkoper,l_p0,film,nop,
      22           0 :      >               invtab,mrot,l_onedimens,tau,
      23           0 :      <               pair_to_do,maptopair,kdiff,l_q,param_file)
      24             : 
      25             :       implicit none
      26             : 
      27             :       TYPE(t_input),INTENT(IN)     :: input
      28             :       TYPE(t_kpts),INTENT(IN)      :: kpts
      29             :       integer,intent(in) :: nop
      30             :       integer,intent(in) :: mrot(3,3,nop)
      31             :       integer,intent(in) :: invtab(nop)
      32             :       integer,intent(in) :: fullnkpts
      33             :       integer,intent(in) :: nntot
      34             :       integer,intent(in) :: bpt(nntot,fullnkpts)
      35             :       integer,intent(in) :: gb(3,nntot,fullnkpts)
      36             :       logical,intent(in) :: l_bzsym
      37             :       integer,intent(in) :: irreduc(fullnkpts)
      38             :       integer,intent(in) :: mapkoper(fullnkpts)
      39             :       logical,intent(in) :: l_p0,film
      40             :       real,intent(in)    :: tau(3,nop)
      41             : 
      42             :       integer,intent(out):: pair_to_do(fullnkpts,nntot)
      43             :       integer,intent(out):: maptopair(3,fullnkpts,nntot)
      44             :       real,intent(out)   :: kdiff(3,nntot)
      45             : 
      46             :       logical, intent(in) :: l_q
      47             :       character(len=20), intent(in) :: param_file
      48             : 
      49             :       integer :: ikpt,ikpt_k,kptibz,ikpt_b
      50             :       integer :: num_pair,kptibz_b
      51           0 :       integer :: index(fullnkpts,nntot)
      52             :       integer :: oper,oper_b,num_rot,num_conj,repo,repo_b
      53             :       integer :: k,kx,ky,repkpt,repkpt_b,repkpt_bb
      54             :       integer :: sign,sign_b,ngis,ngis_b
      55             :       logical :: startloop
      56             :       logical :: l_file,l_onedimens
      57           0 :       real    :: kpoints(3,fullnkpts)
      58             :       real    :: kdiffvec(3)
      59           0 :       integer :: multtab(nop,nop)
      60             :       integer :: fullnkpts_tmp,kr
      61             :       real    :: scale
      62             :       real    :: brot(3)
      63             :       logical :: l_testnosymm,l_nosymm1,l_nosymm2
      64             : 
      65           0 :       index(:,:)=0
      66           0 :       num_pair=0
      67           0 :       num_rot=0
      68           0 :       num_conj=0
      69           0 :       pair_to_do(:,:)=0
      70             : 
      71           0 :       inquire(file='testnosymm',exist=l_testnosymm)
      72             : 
      73             : c-----Test for nonsymmorphic space groups
      74           0 :       if(l_bzsym)then
      75           0 :        do ikpt=1,fullnkpts
      76           0 :          oper=abs(mapkoper(ikpt))
      77           0 :          if( any( abs(tau(:,oper)).gt.1.e-6 ) ) l_testnosymm=.true.
      78             :        enddo
      79             :       endif
      80             : 
      81           0 :       if(l_bzsym)then
      82           0 :          call close_pt(nop,mrot,multtab)
      83             :       endif
      84             : 
      85           0 :       do 10 ikpt = 1,fullnkpts  ! loop by k-points starts
      86           0 :         l_nosymm1=.false. 
      87           0 :         kptibz=ikpt 
      88           0 :         if(l_bzsym) then
      89           0 :            kptibz=irreduc(ikpt)
      90           0 :            oper=mapkoper(ikpt)
      91           0 :            if(oper.lt.0)then
      92           0 :               oper=-oper
      93           0 :               sign=-1
      94             :            else
      95             :               sign=1
      96             :            endif
      97             :            if( 
      98           0 :      &          any( abs(tau(:,oper)).gt.1.e-6 ) 
      99           0 :      &         )l_nosymm1=.true.
     100             :         endif
     101             : 
     102           0 :         do 15 ikpt_b = 1,nntot
     103           0 :          l_nosymm2=.false.  
     104           0 :          if(index(ikpt,ikpt_b).eq.1)cycle
     105           0 :          kptibz_b=bpt(ikpt_b,ikpt)
     106           0 :          if(l_bzsym) then
     107           0 :             oper_b=mapkoper(kptibz_b)
     108           0 :             kptibz_b=irreduc(kptibz_b)
     109           0 :             if(oper_b.lt.0)then
     110           0 :                oper_b=-oper_b
     111           0 :                sign_b=-1
     112             :             else
     113             :                sign_b=1
     114             :             endif
     115             :             if( 
     116           0 :      &          any( abs(tau(:,oper_b)).gt.1.e-6 ) 
     117           0 :      &         )l_nosymm2=.true.
     118             :          endif
     119             : 
     120           0 :          if(l_testnosymm)goto 33
     121           0 :          if(l_nosymm1.or.l_nosymm2)goto 33
     122             : 
     123             : c***************************************************************
     124             : c..the conjugation selection rule, which speeds up significantly
     125             : c***************************************************************
     126           0 :         do ikpt_k = 1,nntot
     127             :          if((bpt(ikpt_k,bpt(ikpt_b,ikpt)).eq.ikpt).and.
     128             :      &       gb(1,ikpt_b,ikpt).eq.(-gb(1,ikpt_k,bpt(ikpt_b,ikpt))).and.
     129             :      &       gb(2,ikpt_b,ikpt).eq.(-gb(2,ikpt_k,bpt(ikpt_b,ikpt))).and.
     130           0 :      &       gb(3,ikpt_b,ikpt).eq.(-gb(3,ikpt_k,bpt(ikpt_b,ikpt))).and.
     131           0 :      &       (index(bpt(ikpt_b,ikpt),ikpt_k).eq.1))then
     132           0 :             index(ikpt,ikpt_b)=1
     133           0 :             maptopair(1,ikpt,ikpt_b)=bpt(ikpt_b,ikpt)
     134           0 :             maptopair(2,ikpt,ikpt_b)=ikpt_k
     135           0 :             maptopair(3,ikpt,ikpt_b)=1
     136             : c            print*,"conjugation"
     137           0 :             num_conj=num_conj+1
     138           0 :            goto 15  
     139             :          endif
     140             :         enddo !ikpt_k
     141             : c****************************************************************
     142             : c     check whether k-point pairs can be mapped onto each other
     143             : c         by rotation
     144             : c****************************************************************
     145           0 :         if(l_bzsym)then
     146             : c         if(all(gb(:,ikpt_b,ikpt).eq.0))then
     147           0 :           do k=1,fullnkpts
     148           0 :            if(irreduc(k).eq.kptibz)then
     149           0 :              repkpt=k
     150           0 :              repo=mapkoper(k)
     151           0 :              if(repo.lt.0)then
     152           0 :                 repo=-repo
     153           0 :                 ngis=-1
     154             :              else
     155             :                 ngis=1
     156             :              endif
     157             : 
     158           0 :              do kx=1,fullnkpts
     159           0 :               if(irreduc(kx).eq.kptibz_b)then
     160           0 :                repkpt_bb=kx
     161           0 :                repo_b=mapkoper(kx)
     162           0 :                if(repo_b.lt.0)then
     163           0 :                   repo_b=-repo_b
     164           0 :                   ngis_b=-1
     165             :                else
     166             :                   ngis_b=1
     167             :                endif 
     168           0 :                do ky=1,nntot
     169           0 :                 if(bpt(ky,repkpt).eq.repkpt_bb)then
     170           0 :                  repkpt_b=ky
     171           0 :                  if (index(repkpt,repkpt_b).eq.1)then
     172           0 :                   if(.not.all(gb(:,ikpt_b,ikpt).eq.0))then
     173           0 :                    brot(:)=0.0
     174           0 :                    do kr=1,3
     175             :                      brot(:)=brot(:)
     176             :      &                 -sign_b*mrot(kr,:,oper)*gb(kr,ikpt_b,ikpt)
     177           0 :      &                 +ngis_b*mrot(kr,:,oper_b)*gb(kr,repkpt_b,repkpt)
     178             :                    enddo
     179           0 :                    if( any(   abs(brot).gt.1e-6       )   )cycle
     180             :                   endif
     181           0 :                   if(sign*ngis*multtab(invtab(oper),repo).eq.
     182             :      &               sign_b*ngis_b*multtab(invtab(oper_b),repo_b))then  
     183           0 :                     maptopair(1,ikpt,ikpt_b)=repkpt
     184           0 :                     maptopair(2,ikpt,ikpt_b)=repkpt_b
     185           0 :                     maptopair(3,ikpt,ikpt_b)=2+(1-ngis*sign)/2
     186           0 :                     index(ikpt,ikpt_b)=1
     187           0 :                     num_rot=num_rot+1
     188           0 :                     goto 15
     189             :                   endif
     190             :                  endif                   
     191             :                 endif   
     192             :                enddo    
     193             :               endif   
     194             :              enddo 
     195             :            endif   
     196             :           enddo   
     197             : c        endif !gb=0   
     198             :         endif
     199             : 
     200             :  33     continue
     201             : 
     202           0 :         index(ikpt,ikpt_b)=1
     203           0 :         num_pair=num_pair+1
     204           0 :         pair_to_do(ikpt,ikpt_b)=num_pair
     205             :         
     206           0 : 15    continue !loop over nearest neighbor k-points
     207           0 : 10    continue ! end of cycle by the k-points  
     208             : 
     209           0 :       if(l_p0)then
     210           0 :       write(6,*)"pairs to calculate: ",num_pair
     211           0 :       write(6,*)"maps by conjugation: ",num_conj
     212           0 :       write(6,*)"maps by rotation:", num_rot
     213           0 :       write(6,*)"num_pair+num_rot+num_conj:",num_pair+num_conj+num_rot
     214           0 :       if(.not.l_q) then
     215           0 :          write(6,*)"fullnkpts*nntot:", fullnkpts*nntot
     216             :       else
     217           0 :           write(6,*)"fullnqpts*nntot:", fullnkpts*nntot
     218             :       endif
     219             :       endif !l_p0
     220             : 
     221             : c*****************************************************************
     222             : c     determine difference vectors that occur on the k-mesh
     223             : c*****************************************************************      
     224           0 :       if (l_bzsym) then
     225           0 :          l_file=.false.
     226           0 :          IF(.NOT.l_q)THEN
     227           0 :            inquire(file='w90kpts',exist=l_file)
     228           0 :            IF(.NOT.l_file)  CALL juDFT_error
     229             :      +        ("w90kpts not found, needed if bzsym",calledby
     230           0 :      +        ="wann_mmnk_symm")
     231           0 :            open(412,file='w90kpts',form='formatted')
     232             :          ELSE
     233           0 :            inquire(file='w90qpts',exist=l_file)
     234           0 :            IF(.NOT.l_file)  CALL juDFT_error
     235             :      +        ("w90qpts not found, needed if bzsym",calledby
     236           0 :      +        ="wann_mmnk_symm")
     237           0 :            open(412,file='w90qpts',form='formatted')            
     238             :          ENDIF
     239           0 :          read(412,*)fullnkpts_tmp,scale
     240           0 :          do k=1,fullnkpts
     241           0 :                read(412,*)kpoints(:,k)
     242             :          enddo   
     243           0 :          kpoints=kpoints/scale
     244           0 :          close(412)
     245             :       else   
     246           0 :             IF(.not.l_q)THEN
     247           0 :               IF(.NOT.input%l_inpXML) THEN
     248           0 :                  open(412,file='kpts',form='formatted')
     249           0 :                  read(412,*)fullnkpts_tmp,scale
     250           0 :                  do k=1,fullnkpts
     251           0 :                     read(412,*)kpoints(:,k)
     252             :                  enddo   
     253           0 :                  kpoints(:,:)=kpoints/scale
     254             :               ELSE
     255           0 :                  fullnkpts_tmp = kpts%nkpt
     256           0 :                  do k=1,fullnkpts
     257           0 :                     kpoints(:,k) = kpts%bk(:,k)
     258             :                  enddo   
     259             :               END IF
     260             :             ELSE
     261           0 :               open(412,file=param_file,form='formatted')
     262           0 :               read(412,*)fullnkpts_tmp,scale
     263           0 :               do k=1,fullnkpts
     264           0 :                  read(412,*)kpoints(:,k)
     265             :               enddo   
     266           0 :               kpoints(:,:)=kpoints/scale
     267             :             ENDIF
     268             : 
     269           0 :             if (film.and..not.l_onedimens) kpoints(3,:)=0.0
     270           0 :             close(412)
     271             :       endif
     272             :       
     273           0 :       if(l_p0)then
     274           0 :          IF(.not.l_q) THEN
     275           0 :            print*,"vectors combining nearest neighbor k-points:"
     276             :          ELSE
     277           0 :            print*,"vectors combining nearest neighbor q-points:"
     278             :          ENDIF
     279             :       endif   
     280           0 :       ky=1
     281           0 :       do k=1,fullnkpts
     282           0 :          do kx=1,nntot
     283           0 :             kdiffvec=kpoints(:,bpt(kx,k))+gb(:,kx,k)-kpoints(:,k)
     284           0 :             do ikpt=1,ky-1
     285           0 :                if(all(abs(kdiff(:,ikpt)-kdiffvec).le.0.0001))goto 200
     286             :             enddo 
     287           0 :             IF(ky>nntot)  CALL juDFT_error("problem in wann_mmnk_symm"
     288           0 :      +           ,calledby ="wann_mmnk_symm")
     289           0 :             kdiff(:,ky)=kdiffvec(:)
     290           0 :             if(l_p0)then
     291           0 :                print*,ky,k,kx,kdiff(:,ky)
     292             :             endif
     293             : 
     294           0 :             ky=ky+1
     295           0 :  200        continue
     296             : 
     297             :          enddo
     298             :       enddo   
     299           0 :       end subroutine
     300           0 :       SUBROUTINE close_pt(
     301           0 :      >                    nops,mrot,
     302           0 :      <                    mtable)
     303             : 
     304             :       IMPLICIT NONE
     305             : 
     306             :       INTEGER, INTENT (IN)  :: nops,mrot(3,3,nops)
     307             :       INTEGER, INTENT (OUT) :: mtable(nops,nops)   ! table(i,j) = {R_i|0}{R_j|0}
     308             : 
     309           0 :       INTEGER              :: i,j,k,mp(3,3),map(nops)
     310             : 
     311             : !---> loop over all operations
     312           0 :       DO j=1,nops
     313             : 
     314           0 :          map(1:nops) = 0
     315             : 
     316             : !--->    multiply {R_j|0}{R_i|0}
     317           0 :          DO i=1,nops
     318           0 :             mp = matmul( mrot(:,:,j) , mrot(:,:,i) )
     319             : 
     320             : !--->       determine which operation this is
     321           0 :             DO k = 1, nops
     322           0 :               IF ( all( mp(:,:)==mrot(:,:,k) ) ) THEN
     323           0 :                  IF ( map(i) .eq. 0 ) THEN
     324           0 :                     map(i) = k
     325             :                  ELSE
     326           0 :                     WRITE (6,'(" Symmetry error : multiple ops")')
     327             :                     CALL juDFT_error("close_pt: Multiple ops (Bravais)"
     328           0 :      +                   ,calledby ="wann_mmnk_symm")
     329             :                  ENDIF
     330             :               ENDIF
     331             :             ENDDO
     332             : 
     333           0 :             IF (map(i).eq.0) THEN
     334           0 :                WRITE (6,'(" Group not closed (Bravais lattice)")')
     335             :                WRITE (6,'(" operation j=",i2,"  map=",12i4,:/,
     336           0 :      &                  (21x,12i4))')  j, map(1:nops)
     337             :                CALL juDFT_error("close_pt: Not closed",calledby
     338           0 :      +              ="wann_mmnk_symm")
     339             :             ENDIF
     340             :          ENDDo
     341           0 :          mtable(j,1:nops) = map(1:nops)
     342             :       ENDDO
     343             : 
     344           0 :       END SUBROUTINE close_pt
     345             : 
     346             :       end module m_wann_mmnk_symm

Generated by: LCOV version 1.13