kpttet.f 10.6 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 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298
      MODULE m_kpttet
      use m_juDFT
      CONTAINS
      SUBROUTINE kpttet(
     >                  iofile,ibfile,iokpt,
     >                  kpri,ktest,kmidtet,mkpt,ndiv3,
     >                  nreg,nfulst,rltv,voluni,
     >                  nsym,ccr,mdir,mface,
     >                  ncorn,nface,fdist,fnorm,cpoint,
     <                  voltet,ntetra,ntet,vktet,
     =                  nkpt,
     <                  divis,vkxyz,wghtkp)
c
c
c ---> This program generates k-points
c           in irreducible wedge of BZ  (for nreg=0)
c           in total BZ                 (for nreg=1)
c      (BZ = 1. Brillouin-zone) for all canonical Bravais lattices
c      in 3 dimensions,
c      using the basis vectors of the reciprocal lattice,
c      the corner points of the irreducible wedge of the BZ
c      and the bordering planes of the irreducible wedge.
c
c      The k-points are generated by the tetrahedron method.
c      by generating a set of k-points which are maximally far apart
c      for the rquired number of points.
c      The method and the subroutines were obtained via St.Bluegel.
c      The information about the irr wedge of the BZ
c      is taken from BRZONE.
c
c-----------------------------------------------------------------------
c    Meaning of variables:
c    INPUT:
c
c    Symmetry of lattice:
c    rltv     : cartesian coordinates of basis vectors for
c               reciprocal lattice: rltv(ix,jn), ix=1,3; jn=1,3
c    voluni   : volume of the Bravais lattice unit cell
c    nsym     : number of symmetry elements of points group
c    ccr     : rotation matrix for symmetry element
c                   in cartesian representation
c
c    representation of the irreducible part of the BZ:
c    fnorm    : normal vector of the planes bordering the irrBZ
c    fdist    : distance vector of the planes bordering the irrBZ
c    ncorn    : number of corners of the irrBZ
c    nface    : number of faces of the irrBZ
c    cpoint   : cartesian coordinates of corner points of irrBZ
c
c    characterization of the tetrahedron-method k-point set:
c    nreg     : 1 kpoints in full BZ; 0 kpoints in irrBZ
c    nfulst   : 1 kpoints ordered in full stars 
c                  (meaningful only for nreg =1; full BZ)
c    nkpt     : on input: required number of k-points inside irrBZ
c               to build the tetrahedrons
c    ntet     : number of tetrahedra generated
c    ntetra   : list of four points for each tetrahedron
c               containing the indices of the respective corner points
c    vktet    : corner points of tetrahedra
c
c    kmidtet  : key to generate mid-tetrahedron k-points
c               1 mid-points are generated; 0 not generated
c
c    OUTPUT: k-point set
c    nkpt     : number of k-points generated in set
c    vkxyz    : vector of kpoint generated; in cartesian representation
c    wghtkp   : weight associated with k-points for BZ integration
c    divis    : integer triple divis(i); i=1,4.
c               Used to find more accurate representation of k-points
c               vklmn(i,kpt)/divis(i) and weights as wght(kpt)/divis(4)
c
c-----------------------------------------------------------------------
      USE m_constants, ONLY : pimach
      USE m_tetcon
      USE m_kvecon
      USE m_fulstar
      IMPLICIT NONE
C
C-----> PARAMETER STATEMENTS
C
      INTEGER, INTENT (IN) :: mkpt,ndiv3,mface,mdir
c
c ---> file number for read and write
c
      INTEGER, INTENT (IN) :: iofile,iokpt,ibfile
c
c ---> running mode parameter
c
      INTEGER, INTENT (IN) :: kpri,ktest
C
C----->  Symmetry information
C
      INTEGER, INTENT (IN) :: nsym
      REAL,    INTENT (IN) :: ccr(3,3,48)
C
C----->  BRAVAIS LATTICE INFORMATION
C
      REAL,    INTENT (IN) ::  voluni
C
C----->  RECIPROCAL LATTICE INFORMATION
C
      INTEGER, INTENT (IN) :: ncorn,nface
      REAL,    INTENT (IN) :: rltv(3,3),fnorm(3,mface),fdist(mface)
      REAL,    INTENT (IN) :: cpoint(3,mface)
C
C----->  BRILLOUINE ZONE INTEGRATION
C
      INTEGER, INTENT (IN) :: nreg,nfulst,kmidtet
      INTEGER, INTENT (INOUT) :: nkpt
      INTEGER, INTENT (OUT) :: ntetra(4,ndiv3),ntet
      REAL,    INTENT (OUT) :: voltet(ndiv3),vktet(3,mkpt)
      REAL,    INTENT (OUT) :: vkxyz(3,mkpt),wghtkp(mkpt),divis(4)

C
C --->  local variables
c
      INTEGER   i,j,ii, nkstar
      REAL      sumwght,eps,one,tpi,sumvol,volirbz
      REAL      vkmid(3,mkpt)
C
C --->  set local constants
c
      SAVE      eps,one
      DATA      eps/1.0e-9/,one/1.0/
c
c======================================================================
c
      tpi = 2.0 * pimach()
c
      IF (kpri.GE.1) THEN
        WRITE (iofile,'(3x,'' *<* kpttet *>* '')')
        WRITE (iofile,'(3x,'' generate k-vectors'')')
        WRITE (iofile,'(3x,'' ~~~~~~~~~~~~~~~~~~'')')
        WRITE (iofile,'(3x,'' by tetrahedron-method'')')
        WRITE (iofile,'(3x,'' ~~~~~~~~~~~~~~~~~~~~~'')')
      ENDIF
 
      WRITE (iofile,'('' k-points generated with tetrahedron '',
     >                                              ''method'')')
      WRITE (iokpt,'(''# k-points generated with tetrahedron '',
     >                                              ''method'')')
      IF (nreg .EQ. 0) THEN
        WRITE (iofile,'(3x,'' in irred wedge of 1. Brillouin zone'')')
        WRITE (iofile,'(3x,'' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'')')
        WRITE (iokpt,'(''#'',i4,21x,''nreg: k-points in irrBZ'')') nreg
      ELSEIF (nreg .eq. 1 .and. nfulst .eq. 1) then
        WRITE (iofile,'(3x,'' in 1. Brillouin zone'')')
        WRITE (iofile,'(3x,'' ~~~~~~~~~~~~~~~~~~~~'')')
        WRITE (iofile,'(3x,'' full stars generated'')')
        WRITE (iofile,'(3x,'' ~~~~~~~~~~~~~~~~~~~~'')')
        WRITE (iokpt,'(''#'',2(i4,1x),14x,'' nreg,nfulst: '',
     >      ''k-points in totBZ, ordered in full stars'')') nreg,nfulst
      ELSE
        WRITE (iofile,'(2(1x,i4),4x,'' nreg,nfulst: wrong choice;'',
     >       /,27x,'' allowed combinations: (1,1); (0,0),(0,1)'')' )
     >                                                   nreg,nfulst
         CALL juDFT_error("nreg,nfulst: wrong choice",calledby="kpttet")
      ENDIF

      CALL kvecon(
     >            iofile,ibfile,mkpt,mface,
     >            nkpt,ncorn,nsym,nface,rltv,fdist,fnorm,cpoint,
     <            vktet )
!
! --->  generate tetrahedra and mid-tetrahedron k-points
!
! --->  (a) Determine the corner K-POINTs for X number of Tetrahedra for
!           doing a very pretty Brillouine zone Integration;
! --->      determine the volume of each tetrahedron
!
      CALL tetcon(
     >            iofile,ibfile,mkpt,ndiv3,
     >            nkpt,voluni,vktet,
     =            nsym,
     <            ntet,voltet,ntetra)
c
      WRITE (iofile,'('' the number of tetrahedra '')')
      WRITE (iofile,*) ntet
      WRITE (iofile,'('' volumes of the tetrahedra '')')
      WRITE (iofile,'(e19.12,1x,i5,5x,''voltet(i),i'')')
     >                               (voltet(i),i,i=1,ntet)
      WRITE (iofile,'('' corners of the tetrahedra '')')
      WRITE (iofile, 999) ((ntetra(j,i),j=1,4),i=1,ntet)
      WRITE (iofile,'('' the # of different k-points '')')
      WRITE (iofile,*) nkpt
      WRITE (iofile,'('' k-points used to construct tetrahedra'')')
      WRITE (iofile,'(3(4x,f10.6))') ((vktet(i,j),i=1,3),j=1,nkpt)
  999 FORMAT (4(3x,4i4))
c
c --->   calculate weights from volume of tetrahedra
c
      volirbz =  tpi**3 /(real(nsym)*voluni)
      sumvol = 0.0
      DO i = 1, ntet
         sumvol = sumvol + voltet(i)
         voltet(i) = ntet * voltet(i) / volirbz 
      ENDDO
c
      IF ((sumvol-volirbz)/volirbz .LE. eps) THEN
        IF (kmidtet.EQ.1) THEN
          DO i = 1, ntet
            wghtkp(i) = voltet(i)/sumvol
          ENDDO
        ELSE
          DO i = 1, nkpt
            wghtkp(i) = 1./nkpt
          ENDDO
        ENDIF
      ELSE
        WRITE (iofile, '(2(e19.12,1x),5x,''summvol.ne.volirbz'')')
     >                                     sumvol,volirbz
         CALL juDFT_error("sumvol =/= volirbz",calledby="kpttet")
      ENDIF
c
c --->  prepare the final set of kpoints in irrBZ (depending on kmidtet)
c
      IF ( kmidtet.EQ.0) THEN
c
        DO i = 1, nkpt
           vkxyz(:,i) = vktet(:,i)
           WRITE (iofile,'(3(f10.7,1x),f12.10,1x,i4,3x,
     +           ''vkxyz, wghtkp'')') (vkxyz(ii,i),ii=1,3),wghtkp(i),i
        ENDDO  
        nkstar = nkpt

      ELSEIF ( kmidtet.EQ.1) THEN
c
c --->  (b) calculate mid-tetrahedron k-points
c
         DO i=1,ntet
            vkmid(:,i) = 0.0
            DO j=1,4
               vkmid(:,i) = vkmid(:,i) + vktet(:,ntetra(j,i))
            ENDDO
            vkmid(:,i) = vkmid(:,i) * 0.25
         ENDDO

         nkpt = ntet
         WRITE (iofile,'('' the new number of k-points is '',i4)') nkpt
         WRITE (iofile,'('' the new k-points are the '',
     +                        ''mid-tetrahedron-points '')')
         WRITE (iokpt,'(''# the new k-points are the '',
     +                        ''mid-tetrahedron-points '')')
         sumwght = 0.00
         DO i=1,ntet
           vkxyz(:,i) = vkmid(:,i)
           sumwght = sumwght + wghtkp(i)
         ENDDO
!
! ---> check sumwght; if abs(sumwght-1).lt.eps print kpoints and weights
!
         IF ( abs(sumwght - one).LT.eps) THEN
            WRITE (iofile,'(1x,f12.10,1x,'' sumwght .eq. one'')')
     +                                                   sumwght
            DO i=1,nkpt
               WRITE (iofile,'(3(f10.7,1x),f12.10,1x,i4,3x,
     +            ''vkxyz, wghtkp'')') (vkxyz(ii,i),ii=1,3),wghtkp(i), i
            ENDDO
            nkstar = ntet

         ELSE
            WRITE (iofile,'(1x,f12.10,1x,'' sumwght .ne. one'')')
     +                                                   sumwght
             CALL juDFT_error("sumwght",calledby="kpttet")
         ENDIF

      END IF ! end of generation of mid-tetrahedron k-points

!
! --->   set denominators for more accurate k-point representation
!
      DO i=1,4
         divis(i) = real(nkpt)
      ENDDO

      IF ( nreg.EQ.1 .AND. nfulst.EQ.1 ) THEN

! --->   generate full stars for all representative k-points
!        - for nreg=1 and nfulst=1:
!              - determine order of full star ifstar(kpn).le.nsym
!              - assign nkpt= sum {ifstar(ik)} (ik=1,ntet)
!              - assign vkxyz(ix,kpn) = vkstar(ix,ikpn(is,ik));
!                        ix=1,3; kpn=1,nkpt; ik=1,ntet; is=1,ifstar(ik)
!              - calculate wghtkp(kpn)=wghtkp_old(ik)/ifstar(ik)
!                                kpn=1,nkpt; ik=1,ntet

         CALL fulstar(
     >                iofile,iokpt,kpri,ktest,
     >                ccr,nsym,
     >                vkxyz,nkstar,mkpt,mface,mdir,
     =                nkpt,vkxyz,wghtkp)

         divis(4) = divis(4) * nsym
      ENDIF

      RETURN
      END SUBROUTINE kpttet
      END MODULE m_kpttet