brzone2.f90 21.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
!--------------------------------------------------------------------------------
! 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.
!--------------------------------------------------------------------------------

MODULE m_brzone2
USE m_juDFT

CONTAINS
SUBROUTINE brzone2(rcmt,nsym,idrot,mface,nbsz,nv48,&
                   cpoint,xvec,ncorn,nedge,nface,fnorm,fdist)

   USE m_constants, ONLY : pimach

   IMPLICIT NONE

   ! This subroutine constructs the irreducible wedge of the Brillouin
   ! zone (IBZ). The subroutine signature is based on an earlier
   ! implementation of this functionality in brzone but the algorithm 
   ! is a new design as the old routine featured too many bugs. Some 
   ! code parts are taken from the old routine.
   !
   !                                  GM, 2016

   INTEGER, INTENT (IN) :: mface,nbsz,nv48
   INTEGER, INTENT (IN) :: nsym               ! number of symmetry elements
   REAL,    INTENT (IN) :: rcmt(3,3)          ! reciprocal lattice basis (2\pi/a.u.)
   REAL,    INTENT (IN) :: idrot(3,3,48)      ! rotation matrices in cartesian repr.

   INTEGER, INTENT (OUT) :: ncorn,nedge,nface ! number of corners, faces and edges of the IBZ
   REAL,    INTENT (OUT) :: fnorm(3,mface)    ! normal vector of the planes bordering the IBZ
   REAL,    INTENT (OUT) :: fdist(mface)      ! distance vector of the planes bordering the IBZ
   REAL,    INTENT (OUT) :: cpoint(3,mface)   ! cartesian coordinates of corner points of IBZ
   REAL,    INTENT (OUT) :: xvec(3)           ! arbitrary vector lying in the IBZ

   INTEGER               :: info, n1, n2, n3, n4, i, j, k, iPlane, iSym
   INTEGER               :: nPlanes, nOuterPlanes, nCorners, numIBZPlanes
   INTEGER               :: maxNumPlanes, nUniqueCorners, duplicateNum
Daniel Wortmann's avatar
Daniel Wortmann committed
40
   INTEGER               :: nbszLocal, planesDim
41
   INTEGER               :: ipiv(3), krecip(3)
Daniel Wortmann's avatar
Daniel Wortmann committed
42
   INTEGER               :: cornerPlanes(3,5000)
43 44
   INTEGER               :: numCornerPlanes(5000)
   LOGICAL               :: foundDuplicate, filterOut
Daniel Wortmann's avatar
Daniel Wortmann committed
45
   LOGICAL               :: isIBZCorner(5000)
46
   REAL, PARAMETER       :: eps09 = 1.0e-9
47
   REAL                  :: pi, maxRecDist, vecDist, recScale, scalarProduct
48
   REAL                  :: norm, normalScalarProduct
49 50
   REAL                  :: denominator, vec1Fac, vec2Fac, edgeDist
   REAL                  :: testVec(3), sk(3), vecA(3), vecB(3), vecC(3)
Daniel Wortmann's avatar
Daniel Wortmann committed
51
   REAL                  :: corners(3,5000)
52
   REAL                  :: planesMatrix(3,3), solutions(3,1)
53
   REAL                  :: equationSystem(3,3)
54 55 56
   REAL                  :: edgeDirec(3), edgeDistVec(3), distVec(3)

   INTEGER, ALLOCATABLE  :: cornerPlaneList(:,:)
57
   INTEGER, ALLOCATABLE  :: planeCorners(:,:)
Daniel Wortmann's avatar
Daniel Wortmann committed
58 59 60 61
   INTEGER, ALLOCATABLE  :: nPlaneCorners(:)
   LOGICAL, ALLOCATABLE  :: isIBZPlane(:)

   REAL, ALLOCATABLE     :: dvec(:,:), ddist(:)
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

   ! Subroutine plan:
   ! 0. Initializations
   ! 1. Construct all possible boundary planes. Remove duplicate planes.
   ! 1.1. Construct boundary planes according to bisections of the lines
   !      to reciprocal lattice vectors <> 0.
   ! 1.2. Construct all boundary planes due to symmetry relations.
   ! 2. Determine all corner points of the IBZ.
   ! 2.1. Determine all possibly relevant intersection points of each 
   !      subset of 3 boundary planes, respectively.
   ! 2.2. Select those points that are on the correct side of all of the
   !      possible boundary planes.
   ! 3. Remove all boundary planes that don't feature at least three 
   !    corner points. These planes are not part of the real boundary.
   !    Also remove the corner points associated with each removed plane.
   ! 4. Construct edges from the set of boundary points. Actually
   !    only the number of edges is needed. I might use the Euler 
   !    characteristic to calculate it.

   ! 0. Initializations

   pi = pimach()

   maxRecDist = 0.0
   DO n1 = -1,1,2
      DO n2 = -1,1,2
         DO n3 = -1,1,2
            testVec(:) = n1*rcmt(:,1)+n2*rcmt(:,2)+n3*rcmt(:,3)
            vecDist = SQRT(testVec(1)**2.0+testVec(2)**2.0+testVec(3)**2.0)
            IF(vecDist.GT.maxRecDist) maxRecDist = vecDist
         END DO
      END DO
   END DO
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
   maxRecDist = 1.01 * maxRecDist

   ! We ignore the nbsz (the number of neighboring cells to consider) passed to the routine
   ! and calculate an own nbszLocal on the basis of maxRecDist.

   nbszLocal = 1
   DO i = -1, 1, 2
      DO j = -1, 1, 2
         DO k = -1, 1, 2
            solutions(1,1) = i * maxRecDist
            solutions(2,1) = j * maxRecDist
            solutions(3,1) = k * maxRecDist
            equationSystem(:,:) = rcmt(:,:)
            ipiv = 0
            info = 0
            CALL DGETRF(3,3, equationSystem,3,ipiv,info)
            CALL DGETRS('N',3,1,equationSystem,3,ipiv,solutions,3,info)
            ! I assume that info == 0: The reciprocal lattice vectors should be linearly independent.
            if(nbszLocal.LT.ABS(solutions(1,1))) nbszLocal = CEILING(ABS(solutions(1,1)))
            if(nbszLocal.LT.ABS(solutions(2,1))) nbszLocal = CEILING(ABS(solutions(2,1)))
            if(nbszLocal.LT.ABS(solutions(3,1))) nbszLocal = CEILING(ABS(solutions(3,1)))
         END DO
      END DO
   END DO

   maxRecDist = 0.5 * maxRecDist
121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142

   ! 1. Construct all possible boundary planes. Remove duplicate planes.
   ! 1.1. Construct boundary planes according to bisections of the lines
   !      to reciprocal lattice vectors <> 0

   ! 1.1.1. Construct an arbitrary point xvec in the IBZ
   ! 1.1.1.1. Determine a "scale" for the point on the basis of the
   !          reciprocal lattice basis.
   DO i = 1,3
      sk(i) = 0.0
      DO j = 1,3
         sk(i)=sk(i)+rcmt(j,i)**2
      END DO
   END DO
   recScale = SQRT(MIN(sk(1),sk(2),sk(3)))*0.1

   ! 1.1.1.2. Use the scale to construct the arbitrary point xvec.
   xvec(1) = recScale
   xvec(2) = recScale/SQRT(pi)
   xvec(3) = recScale/pi

   ! loops over all neighboring lattice vectors. The size of the 
143
   ! neigborhood is defined by nbszLocal.
144

Daniel Wortmann's avatar
Daniel Wortmann committed
145 146 147 148 149
   planesDim = 2*nbszLocal+1
   planesDim = planesDim*planesDim*planesDim + nsym
   ALLOCATE(dvec(3,planesDim), ddist(planesDim))
   ALLOCATE(nPlaneCorners(planesDim),isIBZPlane(planesDim))

150
   nPlanes = 0
151
   DO n1 = -nbszLocal, nbszLocal
152
      krecip(1) = n1
153
      DO n2 = -nbszLocal, nbszLocal
154
         krecip(2) = n2
155
         DO n3 = -nbszLocal, nbszLocal
156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178
            IF ( .NOT.(n1.EQ.0.AND.n2.EQ.0.AND.n3.EQ.0) ) THEN
               krecip(3) = n3
               nPlanes = nPlanes + 1

               ! determine distance vector dvec to the selected reciprocal
               ! lattice vector
               DO i = 1,3
                  dvec(i,nPlanes) = 0.0
                  DO j = 1,3
                     dvec(i,nPlanes) = dvec(i,nPlanes) + rcmt(i,j)*krecip(j)
                  END DO
               END DO

               ! Determine the norm of the distance vector
               norm = 0.0
               DO i = 1,3
                  norm = norm + dvec(i,nPlanes)**2
               END DO
               norm = SQRT(norm)

               ! Deterime the distance for the intersection
               ddist(nPlanes) = 0.5*norm

179 180 181 182
               IF(ddist(nPlanes).GT.maxRecDist) THEN
                  nPlanes = nPlanes - 1
                  CYCLE
               END IF
183 184 185 186 187 188 189 190

               ! Normalize distance vector. Note that this is now the
               ! normal to the intersecting plane. The intersecting plane 
               ! is defined by the normal and the distance.
               ! Note: With the current sign the vector points outwards.
               dvec(:,nPlanes) = dvec(:,nPlanes) / norm

               ! Remove plane again if it is a duplicate of another plane
191 192
               ! or only use the plane nearer to the origin if it is parallel 
               ! to another plane and has the same direction of the normal.
193 194
               foundDuplicate = .FALSE.
               DO iPlane = 1, nPlanes-1
195 196 197
                  IF (ABS(dvec(1,iPlane)-dvec(1,nPlanes)).GT.eps09) CYCLE
                  IF (ABS(dvec(2,iPlane)-dvec(2,nPlanes)).GT.eps09) CYCLE
                  IF (ABS(dvec(3,iPlane)-dvec(3,nPlanes)).GT.eps09) CYCLE
198
                  IF (ddist(nPlanes).LT.ddist(iPlane)) ddist(iPlane) = ddist(nPlanes)
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
                  foundDuplicate = .TRUE.
                  EXIT
               END DO
               IF(foundDuplicate) nPlanes = nPlanes - 1
            END IF
         END DO
      END DO
   END DO

   nOuterPlanes = nPlanes

   ! 1.2. Construct all boundary planes due to symmetry relations.
   !      These are the planes bisecting the line connecting xvec
   !      with an element of the star of xvec ( <> xvec )
   !      Note that the loop over the symmetry operations actually
   !      is a loop over all group element of all the relevant 
   !      symmetry groups. These are not only the generating elements
   !      and therefore it is not neccessary to invert the matrices to
   !      also obtain boundary planes on the other side.

   DO iSym = 2, nsym ! leave out identity symmetry operation
      nPlanes = nPlanes + 1

      ! The origin is part of all the planes determined in this way.
      ! -> distance = 0.
      ddist(nPlanes) = 0.0

      ! Determine the vector connecting xvec and its image after performing
      ! the symmetry operation.
      DO i = 1,3
         dvec(i,nPlanes)=-xvec(i)
         DO j=1,3
            dvec(i,nPlanes) = dvec(i,nPlanes) + idrot(i,j,iSym)*xvec(j)
         END DO
      END DO

      ! Normalize the vector (it is normal to the possible boundary plane).
      ! Note that the vector points away from xvec.
      norm = 0.0
      DO i = 1,3
         norm = norm + dvec(i,nPlanes)**2
      END DO
      norm = SQRT(norm)
      dvec(:,nPlanes)=dvec(:,nPlanes) / norm

      ! Remove plane again if it is a duplicate of another plane
      foundDuplicate = .FALSE.
      DO iPlane = nOuterPlanes+1, nPlanes-1
247 248 249 250
         IF (ABS(ddist(iPlane)-ddist(nPlanes)).GT.eps09) CYCLE
         IF (ABS(dvec(1,iPlane)-dvec(1,nPlanes)).GT.eps09) CYCLE
         IF (ABS(dvec(2,iPlane)-dvec(2,nPlanes)).GT.eps09) CYCLE
         IF (ABS(dvec(3,iPlane)-dvec(3,nPlanes)).GT.eps09) CYCLE
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
         foundDuplicate = .TRUE.
         EXIT
      END DO
      IF(foundDuplicate) nPlanes = nPlanes - 1
   END DO

   ! 2. Determine all corner points of the IBZ.
   ! 2.1. Determine all possibly relevant intersection points of each 
   !      subset of 3 boundary planes, respectively.


   nCorners = 0
   nPlaneCorners = 0
   DO n1 = 1, nOuterPlanes
      DO n2 = n1+1, nPlanes
         ! 2.1.1. Calculate intersection edge between planes n1 and n2
         !        and only consider those n1 and n2 where the cuttig 
         !        edge is not too far from the origin. (This would mean
         !        that the crossing point of these two planes and a third 
         !        plane is never relevant for the construction of the
         !        IBZ.)

         ! The direction of the intersection edge is the cross product
         ! of the normals on planes n1 and n2:
         edgeDirec(1) = dvec(2,n1)*dvec(3,n2) - dvec(3,n1)*dvec(2,n2)
         edgeDirec(2) = dvec(3,n1)*dvec(1,n2) - dvec(1,n1)*dvec(3,n2)
         edgeDirec(3) = dvec(1,n1)*dvec(2,n2) - dvec(2,n1)*dvec(1,n2)

         ! Ignore parallel planes
280 281 282
         IF ((ABS(edgeDirec(1)).LT.eps09).AND.&
             (ABS(edgeDirec(2)).LT.eps09).AND.&
             (ABS(edgeDirec(3)).LT.eps09)) CYCLE
283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331

         ! The distance vector of the intersection edge to the origin is given
         ! by (since dvec is normalized):
         !
         !       ddist(n1) - ddist(n2) * <dvec(:,n1)|dvec(:,n2)>              ddist(n2) - ddist(n1) * <dvec(:,n1)|dvec(:,n2)>
         ! d_e = ----------------------------------------------- dvec(:,n1) + ----------------------------------------------- dvec(:,n2)
         !              1.0 - <dvec(:,n1)|dvec(:,n2)>^2                               1.0 - <dvec(:,n1)|dvec(:,n2)>^2

         normalScalarProduct = dvec(1,n1)*dvec(1,n2)+&
                               dvec(2,n1)*dvec(2,n2)+&
                               dvec(3,n1)*dvec(3,n2)
         denominator = 1.0 - (normalScalarProduct**2.0)
         vec1Fac = ddist(n1) - ddist(n2) * normalScalarProduct
         vec2Fac = ddist(n2) - ddist(n1) * normalScalarProduct
         edgeDistVec(:) = vec1Fac*dvec(:,n1) + vec2Fac*dvec(:,n2)
         edgeDistVec(:) = edgeDistVec(:) / denominator
         edgeDist = SQRT(edgeDistVec(1)**2.0+edgeDistVec(2)**2.0+edgeDistVec(3)**2.0)

         ! Ignore planes if intersection edge is too far from origin
         ! (is this criterion ok?)
         IF (edgeDist.GT.maxRecDist) CYCLE

         ! Now calculate the possibly relevant crossing points
         innerPlaneLoop: DO n3 = n2+1, nPlanes
            ! Set up system of linear equations to determine the crossing point
            DO i = 1, 3
               planesMatrix(1,i) = dvec(i,n1)
               planesMatrix(2,i) = dvec(i,n2)
               planesMatrix(3,i) = dvec(i,n3)
            END DO
            solutions(1,1) = ddist(n1)
            solutions(2,1) = ddist(n2)
            solutions(3,1) = ddist(n3)
            ! Solve system of linear equations and cycle if no solution is found
            ! or solution is not in relevant region.
            ipiv = 0
            info = 0
            CALL DGETRF(3,3, planesMatrix,3,ipiv,info)
            IF(info.NE.0) CYCLE
            CALL DGETRS('N',3,1,planesMatrix,3,ipiv,solutions,3,info)
            IF(info.NE.0) CYCLE
            vecDist = SQRT(solutions(1,1)**2.0+solutions(2,1)**2.0+solutions(3,1)**2.0)
            IF(vecDist.GT.maxRecDist) CYCLE

   ! 2.2. Select those points that are on the correct side of all of the
   !      possible boundary planes.

            DO n4 = 1, nPlanes
               vecDist = dvec(1,n4)*solutions(1,1)+dvec(2,n4)*solutions(2,1)+dvec(3,n4)*solutions(3,1) - ddist(n4)
332
               IF(vecDist.GT.eps09) THEN
333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359
                  CYCLE innerPlaneLoop
               END IF
            END DO

            !Add new crossing point to list and associate it with the 3 planes

            nCorners = nCorners + 1
            corners(:,nCorners) = solutions(:,1)
            cornerPlanes(1,nCorners) = n1
            cornerPlanes(2,nCorners) = n2
            cornerPlanes(3,nCorners) = n3
            nPlaneCorners(n1) = nPlaneCorners(n1) + 1
            nPlaneCorners(n2) = nPlaneCorners(n2) + 1
            nPlaneCorners(n3) = nPlaneCorners(n3) + 1
         END DO innerPlaneLoop
      END DO
   END DO

   ! 3. Remove all boundary planes that don't feature at least three 
   !    corner points. These planes are not part of the real boundary.
   !    Also remove the corner points associated with each removed plane.
   !    ...and also those planes and corners that are irrelevant due to 
   !    other reasons.

   ! Corners might be present multiple times (for different triples of planes). 
   ! Remove duplicates. Collect planes intersecting at the respective point.
   maxNumPlanes = MAX(nPlanes-nOuterPlanes,3*nCorners) ! Just a save guess for the array dimension in the next line
360
   ALLOCATE(cornerPlaneList(maxNumPlanes,nCorners+1)) ! The +1 is because of the origin that is added later.
361 362 363 364 365 366 367 368 369
   cornerPlaneList = 0
   numCornerPlanes = 0

   nUniqueCorners = 0
   DO i = 1, nCorners
      duplicateNum = -1
      compareLoop: DO j = 1, i-1
         distVec(:) = corners(:,i) - corners(:,j)
         vecDist = SQRT(distVec(1)**2+distVec(2)**2+distVec(3)**2)
370
         IF (vecDist.LT.eps09) THEN
371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399
            duplicateNum = j
            EXIT compareLoop
         END IF 
      END DO compareLoop
      IF (duplicateNum.EQ.-1) THEN
         nUniqueCorners = nUniqueCorners + 1
         numCornerPlanes(nUniqueCorners) = 3
         corners(:,nUniqueCorners) = corners(:,i)
         cornerPlaneList(1,nUniqueCorners) = cornerPlanes(1,i)
         cornerPlaneList(2,nUniqueCorners) = cornerPlanes(2,i)
         cornerPlaneList(3,nUniqueCorners) = cornerPlanes(3,i)
      ELSE
         DO j = 1, 3
            foundDuplicate = .FALSE.
            DO k = 1, numCornerPlanes(duplicateNum)
               IF (cornerPlaneList(k,duplicateNum).EQ.cornerPlanes(j,i)) THEN
                  foundDuplicate = .TRUE.
               END IF
            END DO
            IF (.NOT.foundDuplicate) THEN
               numCornerPlanes(duplicateNum) = numCornerPlanes(duplicateNum) + 1
               cornerPlaneList(numCornerPlanes(duplicateNum),duplicateNum) = cornerPlanes(j,i)
            END IF
         END DO
      END IF
   END DO
   nCorners = nUniqueCorners

   ! Add origin to corner points
400 401 402 403 404 405 406 407
   IF ((nPlanes-nOuterPlanes).GE.3) THEN ! The origin is only a corner if at least 3 planes feature this point
      nCorners = nCorners + 1
      corners(:,nCorners) = 0.0
      DO i = nOuterPlanes + 1, nPlanes
         cornerPlaneList(i-nOuterPlanes,nCorners) = i
      END DO
      numCornerPlanes(nCorners) = nPlanes - nOuterPlanes
   END IF
408 409 410 411 412 413 414 415 416 417

   ! Filter out "corners" found for sets of planes that do not meet in a single
   ! point but have a common intersection edge.
   ! (This absurd case actually occured so it has to be treated. LAPACK might
   ! find a "corner point" in this situation. Also the origin can be such a point.)
   nUniqueCorners = 0
   DO i = 1, nCorners
      filterOut = .FALSE.
      IF(numCornerPlanes(i).GE.3) THEN
         filterOut = .TRUE.
418 419 420 421 422 423 424 425 426 427 428 429
         vecA(:) = dvec(:,cornerPlaneList(1,i))
         vecB(:) = dvec(:,cornerPlaneList(2,i))
         cornerPlaneLoop: DO n3 = 3, numCornerPlanes(i)
            vecC(:) = dvec(:,cornerPlaneList(n3,i))
            testVec(1) = vecA(2)*vecB(3)-vecA(3)*vecB(2)
            testVec(2) = vecA(3)*vecB(1)-vecA(1)*vecB(3)
            testVec(3) = vecA(1)*vecB(2)-vecA(2)*vecB(1)
            scalarProduct = testVec(1)*vecC(1) + testVec(2)*vecC(2) + testVec(3)*vecC(3)
            IF (ABS(scalarProduct).GT.eps09) THEN
               filterOut = .FALSE.
               EXIT cornerPlaneLoop
            END IF
430 431 432 433
         END DO cornerPlaneLoop
      END IF
      IF(filterOut) THEN
         CYCLE
434 435 436 437 438 439 440
      END IF
      nUniqueCorners = nUniqueCorners + 1
      IF(nUniqueCorners.NE.i) THEN
         numCornerPlanes(nUniqueCorners) = numCornerPlanes(i)
         DO j = 1, numCornerPlanes(i)
            cornerPlaneList(j,nUniqueCorners) = cornerPlaneList(j,i)
         END DO
441
         corners(:,nUniqueCorners) = corners(:,i)
442 443 444 445 446 447
      END IF
   END DO
   nCorners = nUniqueCorners

   ! Count the number of corners for each plane
   nPlaneCorners = 0
448 449
   ALLOCATE(planeCorners(nPlanes,nCorners))
   planeCorners = 0
450 451 452
   DO i = 1, nCorners
      DO j = 1, numCornerPlanes(i)
         nPlaneCorners(cornerPlaneList(j,i)) = nPlaneCorners(cornerPlaneList(j,i)) + 1
453
         planeCorners(cornerPlaneList(j,i),nPlaneCorners(cornerPlaneList(j,i))) = i
454 455 456 457 458 459 460 461 462 463 464
      END DO
   END DO

   ! Remove irrelevant planes:
   nface = 0
   isIBZPlane(:) = .TRUE.
   DO n1 = 1, nPlanes
      IF (nPlaneCorners(n1).LE.2) THEN
         isIBZPlane(n1) = .FALSE.
         CYCLE
      END IF
465 466 467 468 469 470
!      WRITE(*,*) 'plane ', n1
!      WRITE(*,'(4f20.13)') dvec(:,n1), ddist(n1)
!      WRITE(*,*) 'corners:'
!      DO i = 1, nPlaneCorners(n1)
!         WRITE(*,'(i5,3f20.13)') planeCorners(n1,i), corners(:,planeCorners(n1,i))
!      END DO
471 472 473 474 475 476
      nface = nface + 1
   END DO
   
   ! Remove irrelevant corners:
   ncorn = 0
   isIBZCorner(:) = .TRUE.
477
!   WRITE(*,*) 'IBZ corners:'
478 479 480 481 482 483 484 485 486 487
   DO i = 1, nCorners
      numIBZPlanes = 0
      DO j = 1, numCornerPlanes(i)
         IF(isIBZPlane(cornerPlaneList(j,i))) THEN
            numIBZPlanes = numIBZPlanes + 1
         END IF
      END DO
      IF(numIBZPlanes.LE.2) isIBZCorner(i) = .FALSE.
      IF(.NOT.isIBZCorner(i)) CYCLE
      ncorn = ncorn + 1
488
!      WRITE(*,'(i5,3f20.13)') i, corners(:,i)
489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534
   END DO

   DEALLOCATE(cornerPlaneList)

   ! 4. Construct edges from the set of boundary points. Actually
   !    only the number of edges is needed.

   ! Count number of edges based on number of corners for each face
   ! Each face has the same number of edges as it has corners.
   ! Counting the number of edges in this way counts each edge twice.
   nedge = 0
   DO n1 = 1, nPlanes
      IF(.NOT.isIBZPlane(n1)) CYCLE
      nedge = nedge + nPlaneCorners(n1)
   END DO
   nedge = nedge / 2

   ! 5. Fill the output arrays.

   IF((mface.LT.ncorn).OR.(mface.LT.nface)) THEN
      WRITE(*,*) "mface: ", mface
      WRITE(*,*) "ncorn: ", ncorn
      WRITE(*,*) "nface: ", nface
      WRITE(*,*) "mface has to be larger or equal to nface and ncorn."
      WRITE(*,*) "Its value is hardcoded. Adapt it."
      CALL juDFT_error("mface was chosen too small.",calledby ="brzone2")
   END IF

   cpoint = 0.0
   j = 0
   DO i = 1, nCorners
      IF (.NOT.isIBZCorner(i)) CYCLE
      j = j + 1
      cpoint(:,j) = corners(:,i)
   END DO

   fdist = 0.0
   fnorm = 0.0
   j = 0
   DO i = 1, nPlanes
      IF(.NOT.isIBZPlane(i)) CYCLE
      j = j + 1
      fdist(j) = ddist(i)
      fnorm(:,j) = dvec(:,i)
   END DO

535 536 537 538 539 540 541 542 543 544 545 546
!   WRITE(*,*) 'ncorn', ncorn
!   WRITE(*,*) 'nedge', nedge
!   WRITE(*,*) 'nface', nface
!   WRITE(*,*) 'faces:'
!   DO i  =1,nface
!      WRITE(*,'(4f20.13)') fnorm(:,i), fdist(i)
!   END DO
!   WRITE(*,*) 'coners:'
!   DO i = 1,ncorn
!      WRITE(*,'(3f20.13)') cpoint(:,i)
!   END DO

547 548 549 550 551 552
   ! 5.1. Check Euler characteristic

   IF ((ncorn + nface - nedge).NE.2) THEN
      WRITE(*,*) "ncorn: ", ncorn
      WRITE(*,*) "nface: ", nface
      WRITE(*,*) "nedge: ", nedge
553 554 555 556 557 558 559 560
      WRITE(*,*) "corners: "
      DO i = 1, ncorn
         WRITE(*,'(3f20.13)') cpoint(:,i)
      END DO
      WRITE(*,*) "faces: "
      DO i = 1, nface
         WRITE(*,'(4f20.13)') fnorm(:,i), fdist(i)
      END DO
561 562 563
      CALL juDFT_error("Brillouin zone does not fulfill Euler characterisic.",calledby ="brzone2")
   END IF

Daniel Wortmann's avatar
Daniel Wortmann committed
564
   DEALLOCATE (planeCorners,dvec,ddist,nPlaneCorners,isIBZPlane)
565

566 567 568
END SUBROUTINE brzone2

END MODULE m_brzone2