LCOV - code coverage report
Current view: top level - wannier - wann_mmnk_symm.f (source / functions) Coverage Total Hit
Test: FLEUR test coverage Lines: 44.6 % 168 75
Test Date: 2025-06-14 04:34:23 Functions: 50.0 % 2 1

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

Generated by: LCOV version 2.0-1