kptmop.f 17.4 KB
Newer Older
1 2 3
      MODULE m_kptmop
      use m_juDFT
!-----------------------------------------------------------------------
Daniel Wortmann's avatar
Daniel Wortmann committed
4 5 6 7 8 9 10 11 12 13 14 15 16
!     ---> This program generates k-points
!     in irreducible wedge of BZ  
!     (BZ = 1. Brillouin-zone) for all canonical Bravais lattices
!     in 2 and 3 dimensions,
!     using the basis vectors of the reciprocal lattice
!     and the bordering planes of the irreducible wedge.
!     
!     The k-points are generated by the Monkhorst--Pack method.
!     The information on the bordering planes of the irr wedge
!     is taken from BRZONE.
!     
!     The program checks the compatibility of the dimension and
!     symmetry and of the provided Monkhorst-Pack-parameters.
17 18 19
!-----------------------------------------------------------------------
      CONTAINS
      SUBROUTINE kptmop(
Daniel Wortmann's avatar
Daniel Wortmann committed
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
     >     idsyst,idtype,nmop,
     >     rltv,bltv,nbound,idimens,
     >     xvec,fnorm,fdist,ncorn,nface,nedge,cpoint,
     >     nsym,ccr,rlsymr,talfa,mkpt,mface,mdir,
     <     nkpt,vkxyz,wghtkp)
c     
c     Meaning of variables:
c     INPUT:
c     
c     Symmetry of lattice:
c     idsyst   : crystal system identification in MDDFT programs
c     idtype   : lattice type identification in MDDFT programs
c     bltv     : cartesian coordinates of basis vectors for
c     Bravais lattice: bltv(ix,jn), ix=1,3; jn=1,3
c     rltv     : cartesian coordinates of basis vectors for
c     reciprocal lattice: rltv(ix,jn), ix=1,3; jn=1,3
c     nsym     : number of symmetry elements of points group
c     ccr     : rotation matrix for symmetry element
c     in cartesian representation
c     rlsymr   : rotation matrix for symmetry element
c     in reciprocal lattice basis representation
c     talfa    : translation vector associated with (non-symmorphic)
c     symmetry elements in Bravais lattice 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     nedge    : number of edges of the irrBZ
c     xvec     : arbitrary vector lying in the irrBZ (FOR SURE!!)
c     components are:
c     
c     characterization of the Monkhorst-Pack k-point set:
c     idimens  : number of dimensions for k-point set (2 or 3)
c     nbound   : 0 no primary points on BZ boundary; 
c     1 with boundary points (not for BZ integration!!!)
c     nmop     : integer number triple: nmop(i), i=1,3; nmop(i)
c     determines number of k-points in direction of rltv(ix,i)
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     nkstar   : number of stars for k-points generated in full stars
c     
69 70 71 72 73
c-----------------------------------------------------------------------
      USE m_constants, ONLY : pimach
      USE m_ordstar
      USE m_fulstar
      IMPLICIT NONE
Daniel Wortmann's avatar
Daniel Wortmann committed
74
C     
75
C-----> PARAMETER STATEMENTS
Daniel Wortmann's avatar
Daniel Wortmann committed
76
C     
77
      INTEGER, INTENT (IN) :: mkpt,mface,mdir
Daniel Wortmann's avatar
Daniel Wortmann committed
78
c     
79
C----->  Symmetry information
Daniel Wortmann's avatar
Daniel Wortmann committed
80
C     
81 82 83
      INTEGER, INTENT (IN) :: nsym,idsyst,idtype
      REAL,    INTENT (IN) :: ccr(3,3,48)
      REAL,    INTENT (IN) :: rlsymr(3,3,48),talfa(3,48)
Daniel Wortmann's avatar
Daniel Wortmann committed
84
C     
85
C----->  BRAVAIS LATTICE INFORMATION
Daniel Wortmann's avatar
Daniel Wortmann committed
86
C     
87
      REAL,    INTENT (IN) :: bltv(3,3),cpoint(3,mface)
Daniel Wortmann's avatar
Daniel Wortmann committed
88
C     
89
C----->  RECIPROCAL LATTICE INFORMATION
Daniel Wortmann's avatar
Daniel Wortmann committed
90
C     
91 92 93
      INTEGER, INTENT (IN) :: ncorn,nface,nedge
      REAL,    INTENT (IN) :: xvec(3),rltv(3,3)
      REAL,    INTENT (IN) :: fnorm(3,mface),fdist(mface)
Daniel Wortmann's avatar
Daniel Wortmann committed
94
C     
95
C----->  BRILLOUINE ZONE INTEGRATION
Daniel Wortmann's avatar
Daniel Wortmann committed
96 97
C     
      INTEGER, INTENT (IN) :: nbound,idimens
98
      INTEGER, INTENT (INOUT) :: nmop(3)
Daniel Wortmann's avatar
Daniel Wortmann committed
99 100 101 102 103
      INTEGER, INTENT (OUT):: nkpt
      REAL,    INTENT (OUT):: vkxyz(3,mkpt),wghtkp(mkpt)
C     
C     --->  local variables
c     
104
      CHARACTER*80 blank
Daniel Wortmann's avatar
Daniel Wortmann committed
105
      INTEGER  nkstar,divis(4)
106 107 108 109 110 111 112 113 114 115 116 117 118
      INTEGER  i,idim,i1,i2,i3,ii,ij,ik,is,isym,ifac, iik,iiik
      INTEGER  ikc, i1red,nred,isumkpt,nleft,nirrbz
      INTEGER  dirmin,dirmax,ndir1,ndir2,idir
      INTEGER  kpl,kpm,kpn,nstar(mdir),nstnew
      INTEGER  iplus,iminus,nc2d,n
      REAL     invtpi,zero,one,half,eps,eps1,orient,sum,denom,aivnkpt

      INTEGER  nfract(3),lim(3),isi(3)
      REAL     cp2d(3,mface)
      REAL     vktra(3),vkstar(3,48),ainvnmop(3),fsig(2),vktes(3)
      INTEGER, ALLOCATABLE :: ikpn(:,:),irrkpn(:),nkrep(:)
      INTEGER, ALLOCATABLE :: iside(:),iostar(:)
      REAL,    ALLOCATABLE :: fract(:,:), vkrep(:,:)
Daniel Wortmann's avatar
Daniel Wortmann committed
119 120 121
C     
C     --->  intrinsic functions
c     
122
      INTRINSIC   real,abs
Daniel Wortmann's avatar
Daniel Wortmann committed
123 124 125
C     
C     --->  save and data statements
c     
126 127
      SAVE     one,zero,half,eps,eps1,iplus,iminus
      DATA     zero/0.00/,one/1.00/,half/0.50/,
Daniel Wortmann's avatar
Daniel Wortmann committed
128 129
     +     eps/1.0e-8/,eps1/1.0e-6/,iplus/1/,iminus/-1/
c     
130
c-----------------------------------------------------------------------
Daniel Wortmann's avatar
Daniel Wortmann committed
131
c     
132 133
      ALLOCATE (fract(mkpt,3),vkrep(3,mkpt),ikpn(48,mkpt),irrkpn(mkpt))
      ALLOCATE (nkrep(mkpt),iostar(mkpt),iside(mface))
Daniel Wortmann's avatar
Daniel Wortmann committed
134 135 136 137 138 139 140 141 142 143
      
!     
!     --->   for 2 dimensions only the following Bravais lattices exist:
!     TYPE                    EQUIVALENT 3-DIM        idsyst/idtype
!     square               = p-tetragonal ( 1+2 axis )      2/1
!     rectangular          = p-orthorhomb ( 1+2 axis )      3/1
!     centered rectangular = c-face-orthorhomb( 1+2 axis)   3/6
!     hexagonal            = p-hexagonal  ( 1+2 axis )      4/1
!     oblique              = p-monoclinic ( 1+2 axis )      6/1
!     
144
      IF (idimens .EQ. 2) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
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
!     
!     --->   identify the allowed symmetries
!     and check the consistency of the Monkhorst-Pack-parameters
!     
         IF (idsyst.EQ.2 .OR. idsyst.EQ.4) THEN
            IF (idtype.EQ.1) THEN
               IF (nmop(1).NE.nmop(2) .OR. nmop(3).NE.0) THEN
                  nmop(2) = nmop(1)
                  nmop(3) = 0
                  WRITE (6,'(1x,''WARNING!!!!!!!'',/,
     +''nmop-Parameters not in accordance with symmetry'',/,
     +2(1x,i4),/,
     +'' we have set nmop(2) = nmop(1)'',/,
     +'' and/or nmop(3) = 0'')') idsyst, idtype
                  WRITE (6,'(3(1x,i4),'' new val for nmop: '')')
     +                 (nmop(i),i=1,3)
               ELSE
                  WRITE (6,'('' values accepted unchanged'')')
                  WRITE (6,'(3(1x,i4),14x,''nmop(i),i=1,3'')')
     +                 (nmop(i),i=1,3)
               ENDIF
            ENDIF
         ELSEIF (idsyst.EQ.3) THEN
            IF (idtype.EQ.1 .OR. idtype.EQ.6) THEN
               IF (nmop(3).NE.0) THEN
                  nmop(3) = 0
                  WRITE (6,'(1x,''WARNING!!!!!!!'',/,
     +''nmop-Parameters not in accordance with symmetry'',/,
     +2(1x,i4),/,
     +'' we have set nmop(3) = 0'')') idsyst, idtype
                  WRITE (6,'(3(1x,i4),'' new val for nmop: '')')
     +                 (nmop(i),i=1,3)
               ELSE
                  WRITE (6,'('' values accepted unchanged'')')
                  WRITE (6,'(3(1x,i4),14x,''nmop(i),i=1,3'')')
     +                 (nmop(i),i=1,3)
               ENDIF
            ENDIF
         ELSEIF (idsyst.EQ.6) THEN
            IF (idtype.EQ.1) THEN
               IF (nmop(3).NE.0) THEN
                  nmop(3) = 0
                  WRITE (6,'(1x,''WARNING!!!!!!!'',/,
     +''nmop-Parameters not in accordance with symmetry'',/,
     +2(1x,i4),/,
     +'' we have set nmop(3) = 0'')') idsyst, idtype
                  WRITE (6,'(3(1x,i4),'' new val for nmop: '')')
     +                 (nmop(i),i=1,3)
               ELSE
                  WRITE (6,'('' values accepted unchanged'')')
                  WRITE (6,'(3(1x,i4),14x,''nmop(i),i=1,3'')')
     +                 (nmop(i),i=1,3)
               ENDIF
            ENDIF
         ELSE
!     
!     --->   in all other cases:
!     
            WRITE (6,'(3(1x,i4),20x,'' idimens,idsyst,idtype: '',
     >''wrong choice for 2-dimensional crystal structure'')')
     >           idimens,idsyst,idtype
            CALL juDFT_error("2-dim crystal",calledby="kptmop")
         ENDIF 
!     
!     --->   check consistency of nmop-parameters with crystal symmetry
!     
211 212
      ELSEIF (idimens.EQ.3) THEN
         IF (idsyst.EQ.1 .OR. idsyst.EQ.5) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
213 214 215 216 217 218 219 220 221 222 223 224 225 226 227
            IF (nmop(1).NE.nmop(2) .OR. nmop(1).NE.nmop(3)
     +           .OR. nmop(2).NE.nmop(3)) THEN
               nmop(3) = nmop(1)
               nmop(2) = nmop(1)
               WRITE (6,'(1x,''WARNING!!!!!!!'',/,
     +''nmop-Parameters not in accordance with symmetry'',/,
     +2(1x,i4),/,
     +'' we have set all nmop(i) = nmop(1)'')') idsyst, idtype
               WRITE (6,'(3(1x,i4),'' new val for nmop(i): '')')
     +              (nmop(i),i=1,3)
            ELSE
               WRITE (6,'('' values accepted unchanged'')')
               WRITE (6,'(3(1x,i4),14x,''nmop(i),i=1,3'')')
     +              (nmop(i),i=1,3)
            ENDIF
228
         ELSEIF (idsyst.EQ.2 .OR. idsyst.eq.4) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245
            if((nmop(3).eq.nmop(2)).and.idsyst.eq.2)then
               WRITE (6,'('' values accepted unchanged'')')
               WRITE (6,'(3(1x,i4),14x,''nmop(i),i=1,3'')')
     +              (nmop(i),i=1,3)
            elseif(nmop(1).NE.nmop(2)) THEN
               nmop(2) = nmop(1)
               WRITE (6,'(1x,''WARNING!!!!!!!'',/,
     +''nmop-Parameters not in accordance with symmetry'',/,
     +2(1x,i4),/,
     +'' we have set nmop(2) = nmop(1)'')') idsyst, idtype
               WRITE (6,'(3(1x,i4),'' new val for nmop: '')')
     +              (nmop(i),i=1,3)
            ELSE
               WRITE (6,'('' values accepted unchanged'')')
               WRITE (6,'(3(1x,i4),14x,''nmop(i),i=1,3'')')
     +              (nmop(i),i=1,3)
            ENDIF
246
         ELSEIF (idsyst.LT.1 .OR. idsyst.GT.7) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
247 248 249 250
            WRITE (6,'(1x,''wrong choice of symmetry'',/,
     +2(1x,i4))') idsyst, idtype
            WRITE (6,'(''only values 1.le.idsyst.le.7 allowed'')')
            CALL juDFT_error("wrong idsyst",calledby="kptmop")
251
         ELSE
Daniel Wortmann's avatar
Daniel Wortmann committed
252 253 254
            WRITE (6,'('' values accepted unchanged'')')
            WRITE (6,'(3(1x,i4),11x,''nmop(i),i=1,3'')')
     +           (nmop(i),i=1,3)
255 256
         ENDIF
      ELSE
Daniel Wortmann's avatar
Daniel Wortmann committed
257
         CALL juDFT_error("idimens =/= 2,3 ",calledby="kptmop")
258
      ENDIF
Daniel Wortmann's avatar
Daniel Wortmann committed
259 260 261 262 263 264 265 266 267 268 269 270 271 272
!     
!     --->   start calculation
!     =====================================================================
!     
!     ---> set sign constants
      isi(1) = 0
      isi(2) = iminus
      isi(3) = iplus
!     
!     ---> calc orientation of boundary faces of irr wedge of BZ
!     characterized by
!     iside(i)= sign( (xvec,fnorm(i))-fdist(i) ) ;(i=1,nface )
!     
      WRITE (6,'(1x,''orientation of boundary faces'')')
273 274 275 276 277 278 279 280
      DO ifac = 1, nface
         orient = zero
         iside(ifac) = iplus
         DO ii = 1, 3
            orient = orient + xvec(ii)*fnorm(ii,ifac)
         ENDDO
         orient = orient - fdist(ifac)
         IF (orient .LT. 0) iside(ifac) = iminus
Daniel Wortmann's avatar
Daniel Wortmann committed
281 282
         WRITE (6,'(1x,2(i4,2x),f10.7,10x,''ifac,iside,orient'',
     +'' for xvec'')') ifac,iside(ifac),orient
283 284 285 286
      ENDDO

      invtpi = one / ( 2.0 * pimach() )

Daniel Wortmann's avatar
Daniel Wortmann committed
287
      WRITE (6,'(''Bravais lattice vectors'')' )
288
      DO ii = 1, 3
Daniel Wortmann's avatar
Daniel Wortmann committed
289
         WRITE (6,'(43x,3(1x,f11.6))') (bltv(ii,ikc), ikc=1,3)
290
      ENDDO
Daniel Wortmann's avatar
Daniel Wortmann committed
291
      WRITE (6,'(''reciprocal lattice vectors'')' )
292
      DO ii = 1, 3
Daniel Wortmann's avatar
Daniel Wortmann committed
293
         WRITE (6,'(43x,3(1x,f11.6))' ) (rltv(ii,ikc), ikc=1,3)
294
      ENDDO
Daniel Wortmann's avatar
Daniel Wortmann committed
295 296 297 298 299 300 301 302
!     
!     ---> nmop(i) are Monkhorst-Pack parameters; they determine the
!     number of k-points in i-direction
!     if basis vector lengths are not related by symmetry,
!     we can use independent fractions for each direction
!     
      WRITE (6,'(3(1x,i4),10x,'' Monkhorst-Pack-parameters'')')
     +     (nmop(i1),i1=1,3)
303 304 305 306 307

      DO idim = 1, idimens
         IF (nmop(idim).GT.0) THEN
            ainvnmop(idim) = one/ real(nmop(idim))
         ELSE
Daniel Wortmann's avatar
Daniel Wortmann committed
308 309 310
            WRITE (6,'('' nmop('',i4,'') ='',i4,
     +'' not allowed'')') idim, nmop(idim)
            CALL juDFT_error("nmop wrong",calledby="kptmop")
311 312 313
         ENDIF
      ENDDO

Daniel Wortmann's avatar
Daniel Wortmann committed
314 315 316 317 318 319
      WRITE (6,'(1x,''Monkhorst-Pack-fractions'')' )
!     
!     ---> nbound=1: k-points are generated on boundary of BZ
!     include  fract(1) =       -1/2
!     and  fract(2*nmop+1) = 1/2     for surface points of BZ
!     
320
      IF ( nbound .EQ. 1) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
321 322 323 324 325 326 327 328 329
         WRITE (6,'(1x,i4,10x,''nbound; k-points on boundary'',
     +'' of BZ included'')' ) nbound
!     
!     ---> irregular Monkhorst--Pack--fractions
!     fract(r) = r / (2*nmop)
!     
         DO idim = 1,idimens
            denom = half*ainvnmop(idim)
            divis(idim) = one / denom
330

Daniel Wortmann's avatar
Daniel Wortmann committed
331 332 333 334 335 336 337 338 339 340 341 342 343 344 345
            DO kpn = -nmop(idim),nmop(idim)
               fract(kpn+nmop(idim)+1,idim) = denom * real (kpn)
               WRITE (6,'(10x,f10.7)' ) fract(kpn+nmop(idim)+1,idim)
            ENDDO
            nfract(idim) = 2*nmop(idim) + 1
         ENDDO
         IF (idimens .eq. 2) THEN
            nfract(3) = 1
            fract(1,3) = 0
            divis(3) = one
         END IF
!     
!     ---> nbound=0: k-points are NOT generated on boundary of BZ
!     This is the regular Monkhorst-Pack-method
!     
346
      ELSEIF ( nbound .eq. 0) then
Daniel Wortmann's avatar
Daniel Wortmann committed
347 348 349 350 351 352
         WRITE (6,'(1x,i4,10x,''nbound; no k-points '',
     +'' on boundary of BZ'')' ) nbound
!     
!     --->   regular Monkhorst--Pack--fractions
!     fract(r) =(2*r-nmop-1) / (2*nmop)
!     
353 354 355
         DO idim = 1,idimens
            denom = half*ainvnmop(idim)
            divis(idim) = one / denom
Daniel Wortmann's avatar
Daniel Wortmann committed
356
            WRITE(6,'(5x,i4,5x,''idim'')' ) idim
357
            DO  kpn = 1,nmop(idim)
Daniel Wortmann's avatar
Daniel Wortmann committed
358 359
               fract(kpn,idim) = denom * real (2*kpn -nmop(idim)-1)
               write(6,'(10x,f10.7)' ) fract(kpn,idim)
360 361 362 363 364 365 366 367 368 369 370
            ENDDO
            nfract(idim) = nmop(idim)
         ENDDO

         IF (idimens .EQ. 2) THEN
            nfract(3) = 1
            fract(1,3) = 0
            divis(3) = one
         ENDIF

      ELSE
Daniel Wortmann's avatar
Daniel Wortmann committed
371 372 373
         WRITE (6,'(3x,'' wrong choice of nbound:'', i4)') nbound
         WRITE (6,'(3x,'' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'')')
         CALL juDFT_error("nbound",calledby="kptmop")
374 375
      ENDIF

Daniel Wortmann's avatar
Daniel Wortmann committed
376 377 378 379
!     
!     
!     --->   initialize k-points = zero and weights = 1.0
!     
380 381 382 383 384 385
      DO  kpn = 1,mkpt
         vkxyz(1,kpn) = zero
         vkxyz(2,kpn) = zero
         vkxyz(3,kpn) = zero
         wghtkp(kpn) = one
      ENDDO
Daniel Wortmann's avatar
Daniel Wortmann committed
386 387 388
!     
!     ---> generate equidistant k-vectors in cartesian coordinates
!     
389 390 391 392
      nkpt = 0
      DO i3 = 1,nfract(3)
         DO i2 = 1,nfract(2)
            DO i1 = 1,nfract(1)
Daniel Wortmann's avatar
Daniel Wortmann committed
393 394 395 396 397 398 399 400 401 402 403 404
               nkpt = nkpt + 1
               IF (nkpt>mkpt)  CALL juDFT_error("nkpt > mkpt",calledby
     +              ="kptmop")
               vkxyz(1,nkpt) = rltv(1,1)*fract(i1,1)
     +              + rltv(1,2)*fract(i2,2)
     +              + rltv(1,3)*fract(i3,3)
               vkxyz(2,nkpt) = rltv(2,1)*fract(i1,1)
     +              + rltv(2,2)*fract(i2,2)
     +              + rltv(2,3)*fract(i3,3)
               vkxyz(3,nkpt) = rltv(3,1)*fract(i1,1)
     +              + rltv(3,2)*fract(i2,2)
     +              + rltv(3,3)*fract(i3,3)
405 406 407
            ENDDO
         ENDDO
      ENDDO
Daniel Wortmann's avatar
Daniel Wortmann committed
408 409 410 411 412 413 414
!     
!     --->   calculate weights of k-points and print out k-points
!     wghtkp = 1/nkpt
!     ( = 1/(nmop(1)*nmop(2)*nmop(3)) for reg Monk-Pack-method)
!     
!      divis(4) = real(nkpt)
!      aivnkpt  = one/real(nkpt)
415

Daniel Wortmann's avatar
Daniel Wortmann committed
416 417 418
!      DO  kpn= 1,nkpt
!         wghtkp(kpn) = wghtkp(kpn)*aivnkpt
!      ENDDO
419

Daniel Wortmann's avatar
Daniel Wortmann committed
420 421 422 423 424 425 426 427 428 429 430
!     
!     ====================================================================
!     
!     --->   order generated k-points in stars by applying symmetry:
!     - determine number of different stars nkstar .le. nkpt
!     - determine order of star iostar(kpn) .le. nsym
!     - assign pointer ikpn(i,ik); i=1,iostar(ik); ik=1,nkstar
!     - determine representative vector in irrBZ for each star:
!     vkrep(ix,ik); ix=1,3; ik=1,nkstar
!     
      CALL ordstar(
431
     >     6,0,0,
Daniel Wortmann's avatar
Daniel Wortmann committed
432 433 434 435 436 437 438 439 440 441 442
     >     fnorm,fdist,nface,iside,
     >     nsym,ccr,rltv,mkpt,mface,mdir,
     =     nkpt,vkxyz,
     <     nkstar,iostar,ikpn,vkrep,nkrep)
!     
!     
!     (a) calculate weights for k-points in irrBZ
!     - wghtkp(ik)=iostar(ik)/nkpt_old ; ik=1,nkstar
!     
      DO ik = 1, nkstar
         wghtkp(ik) = wghtkp(ik)*iostar(ik)
443
      ENDDO
Daniel Wortmann's avatar
Daniel Wortmann committed
444 445 446 447 448
!     
!     (b) final preparation of k-points for transfer to file
!     - assign nkpt= nkstar
!     - assign vkxyz(ix,ik) = vkrep(ix,ik); ix=1,3; ik=1,nkstar
!     
449

Daniel Wortmann's avatar
Daniel Wortmann committed
450
         DO i1 = 1,3
451
            DO ik = 1,nkstar
Daniel Wortmann's avatar
Daniel Wortmann committed
452
               vkxyz(i1,ik) = vkrep(i1,ik)
453
            ENDDO
Daniel Wortmann's avatar
Daniel Wortmann committed
454 455 456 457 458
         ENDDO
         nkpt = nkstar
!     
!     --> check for corner points, include them into k-point set:
!     
459
      IF (nbound.EQ.1) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
460 461 462 463 464
         n = 1
         nc2d = 1               ! determine 2D corner points
         cp2d(:,nc2d) = cpoint(:,n)
         corn: DO n = 2, ncorn
         DO i = 1, n-1
465
            IF ((abs(cpoint(1,n)-cpoint(1,i)).LT.0.0001).AND.
Daniel Wortmann's avatar
Daniel Wortmann committed
466 467 468 469 470 471 472 473 474 475 476 477 478 479 480
     +           (abs(cpoint(2,n)-cpoint(2,i)).LT.0.0001)) CYCLE corn
         ENDDO
         nc2d = nc2d + 1
         cp2d(:,nc2d) = cpoint(:,n)
      ENDDO corn
      WRITE (6,'(''2D corner points in internal units'')')
      corn2d: DO n = 1, nc2d 
      WRITE (6,'(i3,3x,2(f10.7,1x))') n,cp2d(1,n),cp2d(2,n)
      DO i = 1, nkpt
         IF ((abs(cp2d(1,n)-vkxyz(1,i)).LT.0.0001).AND.
     +        (abs(cp2d(2,n)-vkxyz(2,i)).LT.0.0001)) CYCLE corn2d
      ENDDO
      nkpt = nkpt + 1
      vkxyz(:,nkpt) = cp2d(:,n)
      ENDDO corn2d
481
      ENDIF 
Daniel Wortmann's avatar
Daniel Wortmann committed
482 483 484
!     
!     --->   print out k-points and weights
!     
485 486 487 488 489
      DEALLOCATE (fract,vkrep,ikpn,irrkpn,nkrep,iostar,iside)

      RETURN
      END SUBROUTINE kptmop
      END MODULE m_kptmop