wann_wigner_seitz.f 2.95 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 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
      module m_wann_wigner_seitz
      use m_juDFT
      contains
      subroutine wann_wigner_seitz(
     >          l_get_rvecnum,num,amat,
     >          rvecnum_in,
     <          rvecnum,rvec,ndegen)

      implicit none
      logical, intent(in) :: l_get_rvecnum
      integer, intent(in) :: num(:)
      real, intent(in)    :: amat(:,:)
      integer,intent(in)  :: rvecnum_in
      integer, intent(out):: rvecnum
      integer, intent(out):: rvec(:,:)
      integer, intent(out):: ndegen(:)	

      integer             :: idist (3)
Daniel Wortmann's avatar
Daniel Wortmann committed
25
      real                :: dist(1331),summa,dist_min
26 27 28 29 30 31 32 33 34 35 36
      integer             :: k1,k2,k3,i1,i2,i3,count,i,j
      real                :: eps7,eps8
      real                :: metric(3,3)

      eps7=1.e-7
      eps8=1.e-8
      rvecnum=0

      metric=matmul(transpose(amat),amat)

      rvecnum = 0  
Daniel Wortmann's avatar
Daniel Wortmann committed
37 38 39
      do k1=-2*num(1),2*num(1)  
       do k2=-2*num(2),2*num(2)  
        do k3=-2*num(3),2*num(3)  
40
         count=0  
Daniel Wortmann's avatar
Daniel Wortmann committed
41 42 43
         do i1=-5,5  
          do i2=-5,5  
           do i3=-5,5  
44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59
            count=count+1  
            ! Get |r-R|^2
            idist(1)=k1-i1*num(1)  
            idist(2)=k2-i2*num(2)  
            idist(3)=k3-i3*num(3)  
            dist(count)=0.0 
            do i=1,3  
             do j=1,3  
              dist(count)=dist(count)+ 
     +         real(idist(i))*metric(i,j)*real(idist(j))
             enddo !i
            enddo !j
           enddo !i3
          enddo !i2
         enddo !i1
         dist_min=minval(dist)
Daniel Wortmann's avatar
Daniel Wortmann committed
60
         if (abs(dist(666) - dist_min ) .lt. eps7 ) then
61 62 63 64 65
                rvecnum = rvecnum + 1  
                if(.not. l_get_rvecnum) then
c                   if(.not.allocated(ndegen))
c     &                 allocate(ndegen(rvecnum_in))
                   ndegen(rvecnum)=0
Daniel Wortmann's avatar
Daniel Wortmann committed
66
                   do i=1,1331
67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85
                      if (abs (dist (i) - dist_min) .lt. eps7 ) 
     &                   ndegen(rvecnum)=ndegen(rvecnum)+1
                   end do
                   rvec(1,rvecnum) = k1  
                   rvec(2,rvecnum) = k2   
                   rvec(3,rvecnum) = k3   
                endif
         endif
        enddo !k3
       enddo !k2
      enddo !k1
      if(l_get_rvecnum) return

c------ Consistency Check.
      summa = 0.0
      do i = 1, rvecnum  
       summa = summa + 1.0/real(ndegen(i))  
      enddo
      if (abs (summa - real(num(1)*num(2)*num(3)) ) > eps8) then
Daniel Wortmann's avatar
Daniel Wortmann committed
86 87
         CALL juDFT_error("problem finding Wigner-Seitz points",
     +                    calledby ="wann_wigner_seitz")
88 89 90 91 92
      endif

      end subroutine wann_wigner_seitz

      end module m_wann_wigner_seitz