LCOV - code coverage report
Current view: top level - wannier - wann_kptsreduc2.f (source / functions) Coverage Total Hit
Test: FLEUR test coverage Lines: 0.0 % 111 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              :       module m_wann_kptsreduc2
       8              :       use m_juDFT
       9              :       contains
      10            0 :       subroutine wann_kptsreduc2(
      11              :      >               mhp,
      12            0 :      >               nop,mrot,bmat,tau,film,
      13              :      >               l_nocosoc)
      14              : c*****************************************************************
      15              : c     Apply the symmetries to reduce the number of k-points.
      16              : c     Frank Freimuth
      17              : c*****************************************************************
      18              : 
      19              :       USE m_constants
      20              : 
      21              :       implicit none
      22              : 
      23              :       integer, intent(in) :: mhp(3)
      24              :       integer, intent(in) :: nop
      25              :       logical,intent(in)  :: film
      26              :       real, intent(in)    :: bmat(3,3)
      27              :       real,intent(in)     :: tau(3,nop)
      28              :       integer, intent(in) :: mrot(3,3,nop)
      29              :       logical,intent(in)  :: l_nocosoc
      30              : 
      31            0 :       real,allocatable    :: weight(:)
      32              :       integer             :: reduznumk,ikpt2
      33            0 :       integer,allocatable :: mapk(:),mapkoper(:)
      34              :       integer             :: ikpt,nmat,k,j
      35              :       integer             :: k1,k2,k3
      36              :       real                :: bbkpt(3)
      37              :       real                :: wk
      38              :       integer             :: brot(3)
      39              :       integer             :: bkpt(3),bkpt2(3)
      40              :       real                :: kpoint(3)
      41              :       integer             :: oper,iter
      42              :       logical             :: l_file
      43              :       real                :: scale
      44              :       integer             :: nkpts
      45              :       integer             :: sumweights
      46            0 :       integer,allocatable :: irreduc(:)
      47            0 :       integer,allocatable :: shiftkpt(:,:)
      48              :       real                :: bbmat(3,3)
      49            0 :       real,allocatable    :: kptslen(:)
      50              :       integer             :: i,kk1,kk2,kk3
      51              : 
      52            0 :       call timestart("wann_kptsreduc2")
      53            0 :       bbmat=matmul(bmat,transpose(bmat))
      54            0 :       write(oUnit,*) "Apply the symmetries to w90kpts"
      55              : c**************************************************************
      56              : c     The array 'mhp' specifies the Monkhorst-Pack mesh.
      57              : c**************************************************************
      58            0 :       nkpts = mhp(1) * mhp(2) * mhp(3)
      59              : 
      60            0 :       allocate(mapk(nkpts),mapkoper(nkpts))
      61            0 :       allocate(weight(nkpts),irreduc(nkpts))
      62              : 
      63              : c*********************************************************
      64              : c     determine lengths of kpoints
      65              : c*********************************************************
      66            0 :       allocate( kptslen(nkpts) )
      67            0 :       iter=0
      68            0 :       do k3=0,mhp(3)-1
      69            0 :        do k2=0,mhp(2)-1
      70            0 :         do k1=0,mhp(1)-1  
      71            0 :          iter=iter+1
      72            0 :          kpoint(1)=real(k1)/real(mhp(1))
      73            0 :          kpoint(2)=real(k2)/real(mhp(2))
      74            0 :          kpoint(3)=real(k3)/real(mhp(3))
      75              :          kptslen(iter)=dot_product( axmintx(kpoint(:)),
      76            0 :      &                matmul(bbmat,axmintx(kpoint(:))) )
      77              :         enddo !k1 
      78              :        enddo !k2
      79              :       enddo !k3
      80              : 
      81              : c***********************************************
      82              : c     determine irreducible part
      83              : c***********************************************
      84            0 :       weight(:)=1.0
      85            0 :       mapk(:)=0
      86            0 :       reduznumk=0
      87            0 :       mapkoper(:)=0
      88            0 :       irreduc(:)=0
      89            0 :       shiftkpt(:,:)=0
      90            0 :       ikpt=0
      91            0 :       do k3=0,mhp(3)-1
      92            0 :        do k2=0,mhp(2)-1
      93            0 :         do k1=0,mhp(1)-1
      94            0 :          ikpt=ikpt+1
      95            0 :          if(mapk(ikpt).ne.0) cycle
      96            0 :          reduznumk=reduznumk+1
      97            0 :          irreduc(ikpt)=reduznumk
      98              : c         bkpt(:)=kpoints(:,ikpt)         
      99            0 :          bkpt(1)=k1*mhp(2)*mhp(3)
     100            0 :          bkpt(2)=k2*mhp(1)*mhp(3)
     101            0 :          bkpt(3)=k3*mhp(1)*mhp(2)
     102            0 :          do oper=1,nop
     103              :           do i=1,3
     104              :            brot(i)=0.0
     105            0 :            do k=1,3
     106              :             brot(i)=brot(i)+mrot(k,i,oper)*bkpt(k)
     107              :            enddo
     108              :           enddo
     109              :           ikpt2=0
     110            0 :           do kk3=0,mhp(3)-1
     111            0 :            do kk2=0,mhp(2)-1
     112            0 :             do kk1=0,mhp(1)-1
     113            0 :              ikpt2=ikpt2+1  
     114            0 :              if(ikpt2.le.ikpt)cycle  
     115            0 :              bkpt2(1)=kk1*mhp(2)*mhp(3)
     116            0 :              bkpt2(2)=kk2*mhp(1)*mhp(3)
     117            0 :              bkpt2(3)=kk3*mhp(1)*mhp(2)
     118            0 :              if(mapk(ikpt2).ne.0)cycle
     119            0 :              if(abs(kptslen(ikpt2)-kptslen(ikpt)).gt.1.e-8)cycle
     120            0 :              kpoint(1)=real(bkpt(1)-bkpt2(1))/real(mhp(1)*mhp(2)*mhp(3))
     121            0 :              kpoint(2)=real(bkpt(2)-bkpt2(2))/real(mhp(1)*mhp(2)*mhp(3))
     122            0 :              kpoint(3)=real(bkpt(3)-bkpt2(3))/real(mhp(1)*mhp(2)*mhp(3))
     123            0 :              if( all(   axmintx(kpoint(:)) .lt.1e-6     ) ) then
     124            0 :                     weight(ikpt)=weight(ikpt)+1.0
     125            0 :                     mapk(ikpt2)=ikpt
     126            0 :                     mapkoper(ikpt2)=oper
     127              :              endif
     128              :             enddo !kk1
     129              :            enddo !kk2
     130              :           enddo !kk3
     131              :          enddo !oper
     132            0 :          if(.not.l_nocosoc)then
     133            0 :           do oper=1,nop
     134              :            do i=1,3
     135              :             brot(i)=0.0
     136            0 :             do k=1,3
     137              :              brot(i)=brot(i)+mrot(k,i,oper)*bkpt(k)
     138              :             enddo
     139              :            enddo
     140              :            ikpt2=0
     141            0 :            do kk3=0,mhp(3)-1
     142            0 :             do kk2=0,mhp(2)-1
     143            0 :              do kk1=0,mhp(1)-1
     144            0 :               ikpt2=ikpt2+1  
     145            0 :               if(ikpt2.le.ikpt)cycle  
     146            0 :               bkpt2(1)=kk1*mhp(2)*mhp(3)
     147            0 :               bkpt2(2)=kk2*mhp(1)*mhp(3)
     148            0 :               bkpt2(3)=kk3*mhp(1)*mhp(2)
     149            0 :               if(mapk(ikpt2).ne.0)cycle
     150            0 :               if(abs(kptslen(ikpt2)-kptslen(ikpt)).gt.1.e-8)cycle
     151              :               kpoint(1)=
     152            0 :      &           real(bkpt(1)-bkpt2(1))/real(mhp(1)*mhp(2)*mhp(3))
     153              :               kpoint(2)=
     154            0 :      &           real(bkpt(2)-bkpt2(2))/real(mhp(1)*mhp(2)*mhp(3))
     155              :               kpoint(3)=
     156            0 :      &           real(bkpt(3)-bkpt2(3))/real(mhp(1)*mhp(2)*mhp(3))
     157            0 :               if( all(  axmintx(kpoint(:)).lt.1e-6  ) ) then
     158            0 :                     weight(ikpt)=weight(ikpt)+1.0
     159            0 :                     mapk(ikpt2)=ikpt
     160            0 :                     mapkoper(ikpt2)=oper
     161              :               endif
     162              :              enddo !kk1
     163              :             enddo !kk2
     164              :            enddo !kk3
     165              :           enddo !oper
     166              :          endif !l_nocosoc 
     167              :         enddo !k1
     168              :        enddo !k2 
     169              :       enddo !k3  
     170              : 
     171              : c****************************************************
     172              : c     write results to files
     173              : c     w90kpts: whole Brillouin zone (for w90)
     174              : c     kpts: irreducible part (for fleur)
     175              : c     kptsmap: mapping from w90kpts to kpts
     176              : c****************************************************
     177              : c      open(117,file='kptsmap',form='formatted',recl=1000)
     178              : c      do ikpt=1,nkpts
     179              : c         if(mapk(ikpt)==0)then
     180              : c            write(117,*)ikpt,irreduc(ikpt),1,0,0,0
     181              : c         else
     182              : c            write(117,*)ikpt,irreduc(mapk(ikpt)),mapkoper(ikpt),
     183              : c     &                  shiftkpt(:,ikpt)
     184              : c         endif  
     185              : c      enddo   
     186              : c      close(117)
     187              : 
     188            0 :       scale=1.000
     189              : 
     190            0 :       open(119,file='kpts',form='formatted')
     191            0 :       if (film)then
     192            0 :         write(119,'(i5,f20.10,3x,l1)')reduznumk,scale,.false.         
     193              :       else
     194            0 :         write(119,'(i5,f20.10)')reduznumk,scale
     195              :       endif   
     196            0 :       sumweights=0
     197              : 
     198            0 :       ikpt=0
     199            0 :       do k3=0,mhp(3)-1
     200            0 :        do k2=0,mhp(2)-1
     201            0 :         do k1=0,mhp(1)-1  
     202            0 :          ikpt=ikpt+1
     203            0 :          kpoint(1)=real(k1)/real(mhp(1))
     204            0 :          kpoint(2)=real(k2)/real(mhp(2))
     205            0 :          kpoint(3)=real(k3)/real(mhp(3))
     206              : 
     207              : 
     208            0 :          if (mapk(ikpt)==0)then
     209            0 :             sumweights=sumweights+weight(ikpt)
     210            0 :             write(oUnit,*)"ikpt=",ikpt
     211            0 :             write(oUnit,*)"irreducible"
     212            0 :             write(oUnit,fmt='(a10,3f9.6)')"internal: ",kpoint(:)/scale
     213              : 
     214            0 :             if (film)then
     215            0 :                write(119,'(3f10.5)')kpoint(1:2),weight(ikpt)
     216              :             else
     217            0 :                write(119,'(4f10.5)')kpoint(:),weight(ikpt)
     218              :             endif   
     219              :       
     220              :          elseif(mapkoper(ikpt).gt.0)then
     221              : c            write(oUnit,*)"ikpt=",ikpt
     222              : c            write(oUnit,*)"reducible"
     223              : c            write(oUnit,*)"map=",mapk(ikpt)
     224              : c            brot(:)=0.0
     225              : c            do k=1,3
     226              : c               brot(:)=brot(:)+
     227              : c     + mrot(k,:,mapkoper(ikpt))*kpoints(k,mapk(ikpt))
     228              : c            enddo   
     229              : c            write(oUnit,'(a19,3f9.6)')"rotated internal: ",brot(:)/scale
     230              :          elseif(mapkoper(ikpt).lt.0)then
     231              : c            write(oUnit,*)"ikpt=",ikpt
     232              : c            write(oUnit,*)"reducible"
     233              : c            write(oUnit,*)"map=",mapk(ikpt)
     234              : c            brot(:)=0.0
     235              : c            do k=1,3
     236              : c               brot(:)=brot(:)-
     237              : c     + mrot(k,:,-mapkoper(ikpt))*kpoints(k,mapk(ikpt))
     238              : c            enddo   
     239              : c            write(oUnit,'(a19,3f9.6)')"rotated internal: ",brot(:)/scale
     240              :          endif   
     241              :         enddo !k1 
     242              :        enddo !k2 
     243              :       enddo !k3  
     244            0 :       close(119)
     245              : 
     246            0 :       write(oUnit,*)"reduznumk=",reduznumk     
     247            0 :       write(oUnit,*)"nkpts=",nkpts
     248            0 :       write(oUnit,*)"sumweights=",sumweights
     249              :       
     250            0 :       IF(sumweights/=nkpts) CALL juDFT_error
     251              :      +     ("sum of weights differs from nkpts",calledby
     252            0 :      +     ="wann_kptsreduc2")
     253            0 :       deallocate(mapk,mapkoper,weight,irreduc)
     254              : 
     255            0 :       call timestop("wann_kptsreduc2")
     256              :       contains
     257            0 :       real elemental function axmintx(x)
     258              :       implicit none
     259              :       real,intent(in) :: x
     260            0 :       axmintx = abs(x-nint(x))
     261              :       end function     
     262              : 
     263              :       end subroutine wann_kptsreduc2
     264              :       end module m_wann_kptsreduc2
        

Generated by: LCOV version 2.0-1