wann_mmnk_symm.f 11 KB
Newer Older
1 2 3 4 5 6
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------

7 8
      module m_wann_mmnk_symm
      use m_juDFT
9
      USE m_types
10 11 12 13 14 15 16 17
      private
      public:: wann_mmnk_symm
      contains
c******************************************************************
c     Find out minimal set of k-point-pairs that have to be
c     calculated; map symmetry-related k-point-pairs to this
c     minimal set.
c     Frank Freimuth
Daniel Wortmann's avatar
Daniel Wortmann committed
18
c******************************************************************
19
      subroutine wann_mmnk_symm(input,kpts,
20 21 22
     >               fullnkpts,nntot,bpt,gb,l_bzsym,
     >               irreduc,mapkoper,l_p0,film,nop,
     >               invtab,mrot,l_onedimens,tau,
Daniel Wortmann's avatar
Daniel Wortmann committed
23
     <               pair_to_do,maptopair,kdiff,l_q,param_file)
24

25 26
      USE m_constants

27
      implicit none
28 29 30

      TYPE(t_input),INTENT(IN)     :: input
      TYPE(t_kpts),INTENT(IN)      :: kpts
31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47
      integer,intent(in) :: nop
      integer,intent(in) :: mrot(3,3,nop)
      integer,intent(in) :: invtab(nop)
      integer,intent(in) :: fullnkpts
      integer,intent(in) :: nntot
      integer,intent(in) :: bpt(nntot,fullnkpts)
      integer,intent(in) :: gb(3,nntot,fullnkpts)
      logical,intent(in) :: l_bzsym
      integer,intent(in) :: irreduc(fullnkpts)
      integer,intent(in) :: mapkoper(fullnkpts)
      logical,intent(in) :: l_p0,film
      real,intent(in)    :: tau(3,nop)

      integer,intent(out):: pair_to_do(fullnkpts,nntot)
      integer,intent(out):: maptopair(3,fullnkpts,nntot)
      real,intent(out)   :: kdiff(3,nntot)

Daniel Wortmann's avatar
Daniel Wortmann committed
48 49 50
      logical, intent(in) :: l_q
      character(len=20), intent(in) :: param_file

51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87
      integer :: ikpt,ikpt_k,kptibz,ikpt_b
      integer :: num_pair,kptibz_b
      integer :: index(fullnkpts,nntot)
      integer :: oper,oper_b,num_rot,num_conj,repo,repo_b
      integer :: k,kx,ky,repkpt,repkpt_b,repkpt_bb
      integer :: sign,sign_b,ngis,ngis_b
      logical :: startloop
      logical :: l_file,l_onedimens
      real    :: kpoints(3,fullnkpts)
      real    :: kdiffvec(3)
      integer :: multtab(nop,nop)
      integer :: fullnkpts_tmp,kr
      real    :: scale
      real    :: brot(3)
      logical :: l_testnosymm,l_nosymm1,l_nosymm2

      index(:,:)=0
      num_pair=0
      num_rot=0
      num_conj=0
      pair_to_do(:,:)=0

      inquire(file='testnosymm',exist=l_testnosymm)

c-----Test for nonsymmorphic space groups
      if(l_bzsym)then
       do ikpt=1,fullnkpts
         oper=abs(mapkoper(ikpt))
         if( any( abs(tau(:,oper)).gt.1.e-6 ) ) l_testnosymm=.true.
       enddo
      endif

      if(l_bzsym)then
         call close_pt(nop,mrot,multtab)
      endif

      do 10 ikpt = 1,fullnkpts  ! loop by k-points starts
Daniel Wortmann's avatar
Daniel Wortmann committed
88 89
        l_nosymm1=.false.
        kptibz=ikpt
90 91 92 93 94 95 96 97 98
        if(l_bzsym) then
           kptibz=irreduc(ikpt)
           oper=mapkoper(ikpt)
           if(oper.lt.0)then
              oper=-oper
              sign=-1
           else
              sign=1
           endif
Daniel Wortmann's avatar
Daniel Wortmann committed
99 100
           if(
     &          any( abs(tau(:,oper)).gt.1.e-6 )
101 102 103 104
     &         )l_nosymm1=.true.
        endif

        do 15 ikpt_b = 1,nntot
Daniel Wortmann's avatar
Daniel Wortmann committed
105
         l_nosymm2=.false.
106 107 108 109 110 111 112 113 114 115 116
         if(index(ikpt,ikpt_b).eq.1)cycle
         kptibz_b=bpt(ikpt_b,ikpt)
         if(l_bzsym) then
            oper_b=mapkoper(kptibz_b)
            kptibz_b=irreduc(kptibz_b)
            if(oper_b.lt.0)then
               oper_b=-oper_b
               sign_b=-1
            else
               sign_b=1
            endif
Daniel Wortmann's avatar
Daniel Wortmann committed
117 118
            if(
     &          any( abs(tau(:,oper_b)).gt.1.e-6 )
119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139
     &         )l_nosymm2=.true.
         endif

         if(l_testnosymm)goto 33
         if(l_nosymm1.or.l_nosymm2)goto 33

c***************************************************************
c..the conjugation selection rule, which speeds up significantly
c***************************************************************
        do ikpt_k = 1,nntot
         if((bpt(ikpt_k,bpt(ikpt_b,ikpt)).eq.ikpt).and.
     &       gb(1,ikpt_b,ikpt).eq.(-gb(1,ikpt_k,bpt(ikpt_b,ikpt))).and.
     &       gb(2,ikpt_b,ikpt).eq.(-gb(2,ikpt_k,bpt(ikpt_b,ikpt))).and.
     &       gb(3,ikpt_b,ikpt).eq.(-gb(3,ikpt_k,bpt(ikpt_b,ikpt))).and.
     &       (index(bpt(ikpt_b,ikpt),ikpt_k).eq.1))then
            index(ikpt,ikpt_b)=1
            maptopair(1,ikpt,ikpt_b)=bpt(ikpt_b,ikpt)
            maptopair(2,ikpt,ikpt_b)=ikpt_k
            maptopair(3,ikpt,ikpt_b)=1
c            print*,"conjugation"
            num_conj=num_conj+1
Daniel Wortmann's avatar
Daniel Wortmann committed
140
           goto 15
141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168
         endif
        enddo !ikpt_k
c****************************************************************
c     check whether k-point pairs can be mapped onto each other
c         by rotation
c****************************************************************
        if(l_bzsym)then
c         if(all(gb(:,ikpt_b,ikpt).eq.0))then
          do k=1,fullnkpts
           if(irreduc(k).eq.kptibz)then
             repkpt=k
             repo=mapkoper(k)
             if(repo.lt.0)then
                repo=-repo
                ngis=-1
             else
                ngis=1
             endif

             do kx=1,fullnkpts
              if(irreduc(kx).eq.kptibz_b)then
               repkpt_bb=kx
               repo_b=mapkoper(kx)
               if(repo_b.lt.0)then
                  repo_b=-repo_b
                  ngis_b=-1
               else
                  ngis_b=1
Daniel Wortmann's avatar
Daniel Wortmann committed
169
               endif
170 171 172 173 174 175 176 177 178 179 180 181 182 183
               do ky=1,nntot
                if(bpt(ky,repkpt).eq.repkpt_bb)then
                 repkpt_b=ky
                 if (index(repkpt,repkpt_b).eq.1)then
                  if(.not.all(gb(:,ikpt_b,ikpt).eq.0))then
                   brot(:)=0.0
                   do kr=1,3
                     brot(:)=brot(:)
     &                 -sign_b*mrot(kr,:,oper)*gb(kr,ikpt_b,ikpt)
     &                 +ngis_b*mrot(kr,:,oper_b)*gb(kr,repkpt_b,repkpt)
                   enddo
                   if( any(   abs(brot).gt.1e-6       )   )cycle
                  endif
                  if(sign*ngis*multtab(invtab(oper),repo).eq.
Daniel Wortmann's avatar
Daniel Wortmann committed
184
     &               sign_b*ngis_b*multtab(invtab(oper_b),repo_b))then
185 186 187 188 189 190 191
                    maptopair(1,ikpt,ikpt_b)=repkpt
                    maptopair(2,ikpt,ikpt_b)=repkpt_b
                    maptopair(3,ikpt,ikpt_b)=2+(1-ngis*sign)/2
                    index(ikpt,ikpt_b)=1
                    num_rot=num_rot+1
                    goto 15
                  endif
Daniel Wortmann's avatar
Daniel Wortmann committed
192 193 194 195 196 197 198 199
                 endif
                endif
               enddo
              endif
             enddo
           endif
          enddo
c        endif !gb=0
200 201 202 203 204 205 206
        endif

 33     continue

        index(ikpt,ikpt_b)=1
        num_pair=num_pair+1
        pair_to_do(ikpt,ikpt_b)=num_pair
Daniel Wortmann's avatar
Daniel Wortmann committed
207

208
15    continue !loop over nearest neighbor k-points
Daniel Wortmann's avatar
Daniel Wortmann committed
209
10    continue ! end of cycle by the k-points
210 211

      if(l_p0)then
212 213 214 215 216
      write(oUnit,*)"pairs to calculate: ",num_pair
      write(oUnit,*)"maps by conjugation: ",num_conj
      write(oUnit,*)"maps by rotation:", num_rot
      write(oUnit,*)"num_pair+num_rot+num_conj:",
     +              num_pair+num_conj+num_rot
Daniel Wortmann's avatar
Daniel Wortmann committed
217
      if(.not.l_q) then
218
         write(oUnit,*)"fullnkpts*nntot:", fullnkpts*nntot
Daniel Wortmann's avatar
Daniel Wortmann committed
219
      else
220
          write(oUnit,*)"fullnqpts*nntot:", fullnkpts*nntot
Daniel Wortmann's avatar
Daniel Wortmann committed
221
      endif
222 223 224 225
      endif !l_p0

c*****************************************************************
c     determine difference vectors that occur on the k-mesh
Daniel Wortmann's avatar
Daniel Wortmann committed
226
c*****************************************************************
227 228
      if (l_bzsym) then
         l_file=.false.
Daniel Wortmann's avatar
Daniel Wortmann committed
229 230 231
         IF(.NOT.l_q)THEN
           inquire(file='w90kpts',exist=l_file)
           IF(.NOT.l_file)  CALL juDFT_error
232 233
     +        ("w90kpts not found, needed if bzsym",calledby
     +        ="wann_mmnk_symm")
Daniel Wortmann's avatar
Daniel Wortmann committed
234 235 236 237 238 239
           open(412,file='w90kpts',form='formatted')
         ELSE
           inquire(file='w90qpts',exist=l_file)
           IF(.NOT.l_file)  CALL juDFT_error
     +        ("w90qpts not found, needed if bzsym",calledby
     +        ="wann_mmnk_symm")
Daniel Wortmann's avatar
Daniel Wortmann committed
240
           open(412,file='w90qpts',form='formatted')
Daniel Wortmann's avatar
Daniel Wortmann committed
241
         ENDIF
242 243 244
         read(412,*)fullnkpts_tmp,scale
         do k=1,fullnkpts
               read(412,*)kpoints(:,k)
Daniel Wortmann's avatar
Daniel Wortmann committed
245
         enddo
246 247
         kpoints=kpoints/scale
         close(412)
Daniel Wortmann's avatar
Daniel Wortmann committed
248
      else
Daniel Wortmann's avatar
Daniel Wortmann committed
249
            IF(.not.l_q)THEN
250 251 252
                 fullnkpts_tmp = kpts%nkpt
                 do k=1,fullnkpts
                    kpoints(:,k) = kpts%bk(:,k)
Daniel Wortmann's avatar
Daniel Wortmann committed
253
                 enddo
Daniel Wortmann's avatar
Daniel Wortmann committed
254 255 256
            ELSE
              open(412,file=param_file,form='formatted')
              read(412,*)fullnkpts_tmp,scale
257 258
              do k=1,fullnkpts
                 read(412,*)kpoints(:,k)
Daniel Wortmann's avatar
Daniel Wortmann committed
259
              enddo
260
              kpoints(:,:)=kpoints/scale
Daniel Wortmann's avatar
Daniel Wortmann committed
261
            ENDIF
262

263 264 265
            if (film.and..not.l_onedimens) kpoints(3,:)=0.0
            close(412)
      endif
Daniel Wortmann's avatar
Daniel Wortmann committed
266

267
      if(l_p0)then
Daniel Wortmann's avatar
Daniel Wortmann committed
268 269 270 271 272
         IF(.not.l_q) THEN
           print*,"vectors combining nearest neighbor k-points:"
         ELSE
           print*,"vectors combining nearest neighbor q-points:"
         ENDIF
Daniel Wortmann's avatar
Daniel Wortmann committed
273
      endif
274 275 276 277 278 279
      ky=1
      do k=1,fullnkpts
         do kx=1,nntot
            kdiffvec=kpoints(:,bpt(kx,k))+gb(:,kx,k)-kpoints(:,k)
            do ikpt=1,ky-1
               if(all(abs(kdiff(:,ikpt)-kdiffvec).le.0.0001))goto 200
Daniel Wortmann's avatar
Daniel Wortmann committed
280
            enddo
281 282 283 284 285 286 287 288 289 290 291
            IF(ky>nntot)  CALL juDFT_error("problem in wann_mmnk_symm"
     +           ,calledby ="wann_mmnk_symm")
            kdiff(:,ky)=kdiffvec(:)
            if(l_p0)then
               print*,ky,k,kx,kdiff(:,ky)
            endif

            ky=ky+1
 200        continue

         enddo
Daniel Wortmann's avatar
Daniel Wortmann committed
292
      enddo
293
      end subroutine
294

295 296 297 298
      SUBROUTINE close_pt(
     >                    nops,mrot,
     <                    mtable)

299 300
      USE m_constants

301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322
      IMPLICIT NONE

      INTEGER, INTENT (IN)  :: nops,mrot(3,3,nops)
      INTEGER, INTENT (OUT) :: mtable(nops,nops)   ! table(i,j) = {R_i|0}{R_j|0}

      INTEGER              :: i,j,k,mp(3,3),map(nops)

!---> loop over all operations
      DO j=1,nops

         map(1:nops) = 0

!--->    multiply {R_j|0}{R_i|0}
         DO i=1,nops
            mp = matmul( mrot(:,:,j) , mrot(:,:,i) )

!--->       determine which operation this is
            DO k = 1, nops
              IF ( all( mp(:,:)==mrot(:,:,k) ) ) THEN
                 IF ( map(i) .eq. 0 ) THEN
                    map(i) = k
                 ELSE
323
                    WRITE (oUnit,'(" Symmetry error : multiple ops")')
324 325 326 327 328 329 330
                    CALL juDFT_error("close_pt: Multiple ops (Bravais)"
     +                   ,calledby ="wann_mmnk_symm")
                 ENDIF
              ENDIF
            ENDDO

            IF (map(i).eq.0) THEN
331 332
               WRITE (oUnit,'(" Group not closed (Bravais lattice)")')
               WRITE (oUnit,'(" operation j=",i2,"  map=",12i4,:/,
333 334 335 336 337 338 339 340 341 342 343
     &                  (21x,12i4))')  j, map(1:nops)
               CALL juDFT_error("close_pt: Not closed",calledby
     +              ="wann_mmnk_symm")
            ENDIF
         ENDDo
         mtable(j,1:nops) = map(1:nops)
      ENDDO

      END SUBROUTINE close_pt

      end module m_wann_mmnk_symm