gkptwgt.f90 3.95 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
      MODULE m_gkptwgt
      CONTAINS
        SUBROUTINE gkptwgt(&
             &                   kpts,cell)
          !                                                             
          !     this subroutine generates the weight factor of a k-point 
          !     in the irreducible wedge of a 2d-brillouin zone         
          !     ver are vertex coordinates in counter clockwise order and 
          !     in units of pi/a(1) and pi/a(2)                              
          !     it is assumed that each k-point has the same 'area'. D.S.Wang     
          !     
          !     changed by                        Stefan Bl"ugel, IFF, Jan.96
          !                                                           
          USE m_constants
          USE m_types
          USE m_juDFT
          IMPLICIT NONE
          TYPE(t_cell),INTENT(IN)   :: cell
          TYPE(t_kpts),INTENT(INOUT):: kpts

          !                                                        
          !     .. was an Argument
          REAL    :: wt
          !     ..
          !     .. Array Arguments ..
          !     ..
          !     .. Local Scalars ..
          REAL cross,dot,eps,x1,x2,y1,y2
          REAL s1,s2
          INTEGER i,j,ikpt,nver
          !     ..
          !     .. Local Arrays ..
33
          REAL ver(2,4),dummy(2,kpts%nkpt)
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
          !     ..
          !     .. Intrinsic Functions ..
          INTRINSIC abs,atan2
          !     ..
          !     .. Data Statements ..
          DATA ver/0.e0,0.e0,1.e0,0.e0,1.e0,1.e0,0.e0,1.e0/
          DATA eps/1.e-6/
          !     ..

          !      nver = 3
          IF ( cell%latnam.EQ.'squ' ) THEN
             nver = 3
          ELSEIF ( cell%latnam.EQ.'p-r' .OR. cell%latnam.EQ.'c-r' ) THEN
             nver = 4
          ELSEIF ( cell%latnam.EQ.'hex' ) THEN
             nver = 3
             ver(2,3) = 1./3.
          ELSEIF ( cell%latnam.EQ.'hx3' .OR. cell%latnam.EQ. 'obl' ) THEN
             CALL juDFT_error("weights for hx3 or obl not defined" ,calledby&
                  &        ="gkptwgt")
          ENDIF
          !                                                          
          !     transform from internal coordinates to xy-coordinates
          !
          !                            changed by shz Feb.96
          DO 10 ikpt = 1 , kpts%nkpt
60
             kpts%wtkpt(ikpt) = 0
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 87 88
             dummy(1,ikpt)=kpts%bk(1,ikpt)
             dummy(2,ikpt)=kpts%bk(2,ikpt)
             s1 = 0.0
             s2 = 0.0
             DO i = 1,2 
                s1 = s1+cell%bmat(i,1)*kpts%bk(i,ikpt)
                s2 = s2+cell%bmat(i,2)*kpts%bk(i,ikpt)
             ENDDO
             IF (cell%latnam.EQ.'hex') THEN
                kpts%bk(1,ikpt) = s1*cell%amat(2,2)/tpi_const
                kpts%bk(2,ikpt) = s2*cell%amat(1,1)/pi_const
             ELSE
                kpts%bk(1,ikpt) = s1*cell%amat(1,1)/pi_const
                kpts%bk(2,ikpt) = s2*cell%amat(2,2)/pi_const
             ENDIF
10        ENDDO

          DO 20 ikpt = 1 , kpts%nkpt
             DO 30 i = 1,nver
                x1 = ( ver(1,i)-kpts%bk(1,ikpt) ) / cell%amat(1,1)
                y1 = ( ver(2,i)-kpts%bk(2,ikpt) ) / cell%amat(2,2)
                j  = i + 1
                IF ( j.GT.nver ) j = 1
                x2 = ( ver(1,j)-kpts%bk(1,ikpt) ) / cell%amat(1,1)
                y2 = ( ver(2,j)-kpts%bk(2,ikpt) ) / cell%amat(2,2)
                dot = x1*x2 + y1*y2
                cross = x1*y2 - y1*x2
                IF ( ABS(cross).GE.eps ) THEN
89
                   kpts%wtkpt(ikpt) = kpts%wtkpt(ikpt) + ATAN2(cross,dot)
90 91 92 93 94
                ENDIF
30           ENDDO
20        ENDDO
          !
          DO ikpt = 1 , kpts%nkpt
95
             kpts%wtkpt(ikpt) = kpts%wtkpt(ikpt) /tpi_const
96 97 98 99
          ENDDO
          !   
          wt = 0.0
          DO ikpt = 1,kpts%nkpt
100
             wt = wt + kpts%wtkpt(ikpt)
101 102 103 104 105 106 107
             kpts%bk(1,ikpt)=dummy(1,ikpt)
             kpts%bk(2,ikpt)=dummy(2,ikpt)
          ENDDO

          RETURN
        END SUBROUTINE gkptwgt
      END