wann_wigner_seitz.f 2.53 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 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
      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)
      real                :: dist(125),summa,dist_min
      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  
      do k1=-num(1),num(1)  
       do k2=-num(2),num(2)  
        do k3=-num(3),num(3)  
         count=0  
         do i1=-2,2  
          do i2=-2,2  
           do i3=-2,2  
            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)
         if (abs(dist(63) - dist_min ) .lt. eps7 ) then
                rvecnum = rvecnum + 1  
                if(.not. l_get_rvecnum) then
c                   if(.not.allocated(ndegen))
c     &                 allocate(ndegen(rvecnum_in))
                   ndegen(rvecnum)=0
                   do i=1,125
                      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
         CALL juDFT_error("problem finding Wigner-Seitz points",calledby
     +        ="wann_wigner_seitz")
      endif

      end subroutine wann_wigner_seitz

      end module m_wann_wigner_seitz