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

Generated by: LCOV version 1.13