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_wigner_seitz
       8             :       use m_juDFT
       9             :       contains
      10           0 :       subroutine wann_wigner_seitz(
      11           0 :      >          l_get_rvecnum,num,amat,
      12             :      >          rvecnum_in,
      13           0 :      <          rvecnum,rvec,ndegen)
      14             : 
      15             :       implicit none
      16             :       logical, intent(in) :: l_get_rvecnum
      17             :       integer, intent(in) :: num(:)
      18             :       real, intent(in)    :: amat(:,:)
      19             :       integer,intent(in)  :: rvecnum_in
      20             :       integer, intent(out):: rvecnum
      21             :       integer, intent(out):: rvec(:,:)
      22             :       integer, intent(out):: ndegen(:)  
      23             : 
      24             :       integer             :: idist (3)
      25             :       real                :: dist(1331),summa,dist_min
      26             :       integer             :: k1,k2,k3,i1,i2,i3,count,i,j
      27             :       real                :: eps7,eps8
      28             :       real                :: metric(3,3)
      29             : 
      30           0 :       call timestart("wann_wigner_seitz")
      31           0 :       eps7=1.e-7
      32           0 :       eps8=1.e-8
      33           0 :       rvecnum=0
      34             : 
      35           0 :       metric=matmul(transpose(amat),amat)
      36             : 
      37           0 :       rvecnum = 0  
      38           0 :       do k1=-2*num(1),2*num(1)  
      39           0 :        do k2=-2*num(2),2*num(2)  
      40           0 :         do k3=-2*num(3),2*num(3)  
      41             :          count=0  
      42           0 :          do i1=-5,5  
      43           0 :           do i2=-5,5  
      44           0 :            do i3=-5,5  
      45           0 :             count=count+1  
      46             :             ! Get |r-R|^2
      47           0 :             idist(1)=k1-i1*num(1)  
      48           0 :             idist(2)=k2-i2*num(2)  
      49           0 :             idist(3)=k3-i3*num(3)  
      50           0 :             dist(count)=0.0 
      51           0 :             do i=1,3  
      52           0 :              do j=1,3  
      53             :               dist(count)=dist(count)+ 
      54           0 :      +         real(idist(i))*metric(i,j)*real(idist(j))
      55             :              enddo !i
      56             :             enddo !j
      57             :            enddo !i3
      58             :           enddo !i2
      59             :          enddo !i1
      60           0 :          dist_min=minval(dist)
      61           0 :          if (abs(dist(666) - dist_min ) .lt. eps7 ) then
      62           0 :                 rvecnum = rvecnum + 1  
      63           0 :                 if(.not. l_get_rvecnum) then
      64             : c                   if(.not.allocated(ndegen))
      65             : c     &                 allocate(ndegen(rvecnum_in))
      66           0 :                    ndegen(rvecnum)=0
      67           0 :                    do i=1,1331
      68           0 :                       if (abs (dist (i) - dist_min) .lt. eps7 ) 
      69           0 :      &                   ndegen(rvecnum)=ndegen(rvecnum)+1
      70             :                    end do
      71           0 :                    rvec(1,rvecnum) = k1  
      72           0 :                    rvec(2,rvecnum) = k2   
      73           0 :                    rvec(3,rvecnum) = k3   
      74             :                 endif
      75             :          endif
      76             :         enddo !k3
      77             :        enddo !k2
      78             :       enddo !k1
      79           0 :       if(l_get_rvecnum) then
      80           0 :             call timestop("wann_wigner_seitz")
      81           0 :             return
      82             :       endif
      83             : 
      84             : c------ Consistency Check.
      85           0 :       summa = 0.0
      86           0 :       do i = 1, rvecnum  
      87           0 :        summa = summa + 1.0/real(ndegen(i))  
      88             :       enddo
      89           0 :       if (abs (summa - real(num(1)*num(2)*num(3)) ) > eps8) then
      90             :          CALL juDFT_error("problem finding Wigner-Seitz points",
      91           0 :      +                    calledby ="wann_wigner_seitz")
      92             :       endif
      93             : 
      94           0 :       call timestop("wann_wigner_seitz")
      95             :       end subroutine wann_wigner_seitz
      96             : 
      97             :       end module m_wann_wigner_seitz