LCOV - code coverage report
Current view: top level - kpoints - brzone2.f90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 226 246 91.9 %
Date: 2019-09-08 04:53:50 Functions: 1 1 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3             : ! This file is part of FLEUR and available as free software under the conditions
       4             : ! of the MIT license as expressed in the LICENSE file in more detail.
       5             : !--------------------------------------------------------------------------------
       6             : 
       7             : MODULE m_brzone2
       8             : USE m_juDFT
       9             : 
      10             : CONTAINS
      11          26 : SUBROUTINE brzone2(rcmt,nsym,idrot,mface,nbsz,nv48,&
      12          26 :                    cpoint,xvec,ncorn,nedge,nface,fnorm,fdist)
      13             : 
      14             :    USE m_constants, ONLY : pimach
      15             : 
      16             :    IMPLICIT NONE
      17             : 
      18             :    ! This subroutine constructs the irreducible wedge of the Brillouin
      19             :    ! zone (IBZ). The subroutine signature is based on an earlier
      20             :    ! implementation of this functionality in brzone but the algorithm 
      21             :    ! is a new design as the old routine featured too many bugs. Some 
      22             :    ! code parts are taken from the old routine.
      23             :    !
      24             :    !                                  GM, 2016
      25             : 
      26             :    INTEGER, INTENT (IN) :: mface,nbsz,nv48
      27             :    INTEGER, INTENT (IN) :: nsym               ! number of symmetry elements
      28             :    REAL,    INTENT (IN) :: rcmt(3,3)          ! reciprocal lattice basis (2\pi/a.u.)
      29             :    REAL,    INTENT (IN) :: idrot(3,3,48)      ! rotation matrices in cartesian repr.
      30             : 
      31             :    INTEGER, INTENT (OUT) :: ncorn,nedge,nface ! number of corners, faces and edges of the IBZ
      32             :    REAL,    INTENT (OUT) :: fnorm(3,mface)    ! normal vector of the planes bordering the IBZ
      33             :    REAL,    INTENT (OUT) :: fdist(mface)      ! distance vector of the planes bordering the IBZ
      34             :    REAL,    INTENT (OUT) :: cpoint(3,mface)   ! cartesian coordinates of corner points of IBZ
      35             :    REAL,    INTENT (OUT) :: xvec(3)           ! arbitrary vector lying in the IBZ
      36             : 
      37             :    INTEGER               :: info, n1, n2, n3, n4, i, j, k, iPlane, iSym
      38             :    INTEGER               :: nPlanes, nOuterPlanes, nCorners, numIBZPlanes
      39             :    INTEGER               :: maxNumPlanes, nUniqueCorners, duplicateNum
      40             :    INTEGER               :: nbszLocal, planesDim
      41             :    INTEGER               :: ipiv(3), krecip(3)
      42             :    INTEGER               :: cornerPlanes(3,5000)
      43             :    INTEGER               :: numCornerPlanes(5000)
      44             :    LOGICAL               :: foundDuplicate, filterOut
      45             :    LOGICAL               :: isIBZCorner(5000)
      46             :    REAL, PARAMETER       :: eps09 = 1.0e-9
      47             :    REAL                  :: pi, maxRecDist, vecDist, recScale, scalarProduct
      48             :    REAL                  :: norm, normalScalarProduct
      49             :    REAL                  :: denominator, vec1Fac, vec2Fac, edgeDist
      50             :    REAL                  :: testVec(3), sk(3), vecA(3), vecB(3), vecC(3)
      51             :    REAL                  :: corners(3,5000)
      52             :    REAL                  :: planesMatrix(3,3), solutions(3,1)
      53             :    REAL                  :: equationSystem(3,3)
      54             :    REAL                  :: edgeDirec(3), edgeDistVec(3), distVec(3)
      55             : 
      56          26 :    INTEGER, ALLOCATABLE  :: cornerPlaneList(:,:)
      57          26 :    INTEGER, ALLOCATABLE  :: planeCorners(:,:)
      58          26 :    INTEGER, ALLOCATABLE  :: nPlaneCorners(:)
      59          26 :    LOGICAL, ALLOCATABLE  :: isIBZPlane(:)
      60             : 
      61          26 :    REAL, ALLOCATABLE     :: dvec(:,:), ddist(:)
      62             : 
      63             :    ! Subroutine plan:
      64             :    ! 0. Initializations
      65             :    ! 1. Construct all possible boundary planes. Remove duplicate planes.
      66             :    ! 1.1. Construct boundary planes according to bisections of the lines
      67             :    !      to reciprocal lattice vectors <> 0.
      68             :    ! 1.2. Construct all boundary planes due to symmetry relations.
      69             :    ! 2. Determine all corner points of the IBZ.
      70             :    ! 2.1. Determine all possibly relevant intersection points of each 
      71             :    !      subset of 3 boundary planes, respectively.
      72             :    ! 2.2. Select those points that are on the correct side of all of the
      73             :    !      possible boundary planes.
      74             :    ! 3. Remove all boundary planes that don't feature at least three 
      75             :    !    corner points. These planes are not part of the real boundary.
      76             :    !    Also remove the corner points associated with each removed plane.
      77             :    ! 4. Construct edges from the set of boundary points. Actually
      78             :    !    only the number of edges is needed. I might use the Euler 
      79             :    !    characteristic to calculate it.
      80             : 
      81             :    ! 0. Initializations
      82             : 
      83          26 :    pi = pimach()
      84             : 
      85          26 :    maxRecDist = 0.0
      86          78 :    DO n1 = -1,1,2
      87         182 :       DO n2 = -1,1,2
      88         364 :          DO n3 = -1,1,2
      89         208 :             testVec(:) = n1*rcmt(:,1)+n2*rcmt(:,2)+n3*rcmt(:,3)
      90         208 :             vecDist = SQRT(testVec(1)**2.0+testVec(2)**2.0+testVec(3)**2.0)
      91         312 :             IF(vecDist.GT.maxRecDist) maxRecDist = vecDist
      92             :          END DO
      93             :       END DO
      94             :    END DO
      95          26 :    maxRecDist = 1.01 * maxRecDist
      96             : 
      97             :    ! We ignore the nbsz (the number of neighboring cells to consider) passed to the routine
      98             :    ! and calculate an own nbszLocal on the basis of maxRecDist.
      99             : 
     100          26 :    nbszLocal = 1
     101          78 :    DO i = -1, 1, 2
     102         182 :       DO j = -1, 1, 2
     103         364 :          DO k = -1, 1, 2
     104         208 :             solutions(1,1) = i * maxRecDist
     105         208 :             solutions(2,1) = j * maxRecDist
     106         208 :             solutions(3,1) = k * maxRecDist
     107         208 :             equationSystem(:,:) = rcmt(:,:)
     108         208 :             ipiv = 0
     109         208 :             info = 0
     110         208 :             CALL DGETRF(3,3, equationSystem,3,ipiv,info)
     111         208 :             CALL DGETRS('N',3,1,equationSystem,3,ipiv,solutions,3,info)
     112             :             ! I assume that info == 0: The reciprocal lattice vectors should be linearly independent.
     113         208 :             if(nbszLocal.LT.ABS(solutions(1,1))) nbszLocal = CEILING(ABS(solutions(1,1)))
     114         208 :             if(nbszLocal.LT.ABS(solutions(2,1))) nbszLocal = CEILING(ABS(solutions(2,1)))
     115         312 :             if(nbszLocal.LT.ABS(solutions(3,1))) nbszLocal = CEILING(ABS(solutions(3,1)))
     116             :          END DO
     117             :       END DO
     118             :    END DO
     119             : 
     120          26 :    maxRecDist = 0.5 * maxRecDist
     121             : 
     122             :    ! 1. Construct all possible boundary planes. Remove duplicate planes.
     123             :    ! 1.1. Construct boundary planes according to bisections of the lines
     124             :    !      to reciprocal lattice vectors <> 0
     125             : 
     126             :    ! 1.1.1. Construct an arbitrary point xvec in the IBZ
     127             :    ! 1.1.1.1. Determine a "scale" for the point on the basis of the
     128             :    !          reciprocal lattice basis.
     129         104 :    DO i = 1,3
     130          78 :       sk(i) = 0.0
     131         338 :       DO j = 1,3
     132         312 :          sk(i)=sk(i)+rcmt(j,i)**2
     133             :       END DO
     134             :    END DO
     135          26 :    recScale = SQRT(MIN(sk(1),sk(2),sk(3)))*0.1
     136             : 
     137             :    ! 1.1.1.2. Use the scale to construct the arbitrary point xvec.
     138          26 :    xvec(1) = recScale
     139          26 :    xvec(2) = recScale/SQRT(pi)
     140          26 :    xvec(3) = recScale/pi
     141             : 
     142             :    ! loops over all neighboring lattice vectors. The size of the 
     143             :    ! neigborhood is defined by nbszLocal.
     144             : 
     145          26 :    planesDim = 2*nbszLocal+1
     146          26 :    planesDim = planesDim*planesDim*planesDim + nsym
     147          26 :    ALLOCATE(dvec(3,planesDim), ddist(planesDim))
     148          26 :    ALLOCATE(nPlaneCorners(planesDim),isIBZPlane(planesDim))
     149             : 
     150          26 :    nPlanes = 0
     151         258 :    DO n1 = -nbszLocal, nbszLocal
     152         232 :       krecip(1) = n1
     153        2636 :       DO n2 = -nbszLocal, nbszLocal
     154        2378 :          krecip(2) = n2
     155       30442 :          DO n3 = -nbszLocal, nbszLocal
     156       30210 :             IF ( .NOT.(n1.EQ.0.AND.n2.EQ.0.AND.n3.EQ.0) ) THEN
     157       27806 :                krecip(3) = n3
     158       27806 :                nPlanes = nPlanes + 1
     159             : 
     160             :                ! determine distance vector dvec to the selected reciprocal
     161             :                ! lattice vector
     162      111224 :                DO i = 1,3
     163       83418 :                   dvec(i,nPlanes) = 0.0
     164      361478 :                   DO j = 1,3
     165      333672 :                      dvec(i,nPlanes) = dvec(i,nPlanes) + rcmt(i,j)*krecip(j)
     166             :                   END DO
     167             :                END DO
     168             : 
     169             :                ! Determine the norm of the distance vector
     170             :                norm = 0.0
     171      194642 :                DO i = 1,3
     172      111224 :                   norm = norm + dvec(i,nPlanes)**2
     173             :                END DO
     174       27806 :                norm = SQRT(norm)
     175             : 
     176             :                ! Deterime the distance for the intersection
     177       27806 :                ddist(nPlanes) = 0.5*norm
     178             : 
     179       27806 :                IF(ddist(nPlanes).GT.maxRecDist) THEN
     180             :                   nPlanes = nPlanes - 1
     181             :                   CYCLE
     182             :                END IF
     183             : 
     184             :                ! Normalize distance vector. Note that this is now the
     185             :                ! normal to the intersecting plane. The intersecting plane 
     186             :                ! is defined by the normal and the distance.
     187             :                ! Note: With the current sign the vector points outwards.
     188        1186 :                dvec(:,nPlanes) = dvec(:,nPlanes) / norm
     189             : 
     190             :                ! Remove plane again if it is a duplicate of another plane
     191             :                ! or only use the plane nearer to the origin if it is parallel 
     192             :                ! to another plane and has the same direction of the normal.
     193             :                foundDuplicate = .FALSE.
     194       59572 :                DO iPlane = 1, nPlanes-1
     195       29271 :                   IF (ABS(dvec(1,iPlane)-dvec(1,nPlanes)).GT.eps09) CYCLE
     196        3535 :                   IF (ABS(dvec(2,iPlane)-dvec(2,nPlanes)).GT.eps09) CYCLE
     197         555 :                   IF (ABS(dvec(3,iPlane)-dvec(3,nPlanes)).GT.eps09) CYCLE
     198          78 :                   IF (ddist(nPlanes).LT.ddist(iPlane)) ddist(iPlane) = ddist(nPlanes)
     199             :                   foundDuplicate = .TRUE.
     200       30301 :                   EXIT
     201             :                END DO
     202             :                IF(foundDuplicate) nPlanes = nPlanes - 1
     203             :             END IF
     204             :          END DO
     205             :       END DO
     206             :    END DO
     207             : 
     208          26 :    nOuterPlanes = nPlanes
     209             : 
     210             :    ! 1.2. Construct all boundary planes due to symmetry relations.
     211             :    !      These are the planes bisecting the line connecting xvec
     212             :    !      with an element of the star of xvec ( <> xvec )
     213             :    !      Note that the loop over the symmetry operations actually
     214             :    !      is a loop over all group element of all the relevant 
     215             :    !      symmetry groups. These are not only the generating elements
     216             :    !      and therefore it is not neccessary to invert the matrices to
     217             :    !      also obtain boundary planes on the other side.
     218             : 
     219         419 :    DO iSym = 2, nsym ! leave out identity symmetry operation
     220         393 :       nPlanes = nPlanes + 1
     221             : 
     222             :       ! The origin is part of all the planes determined in this way.
     223             :       ! -> distance = 0.
     224         393 :       ddist(nPlanes) = 0.0
     225             : 
     226             :       ! Determine the vector connecting xvec and its image after performing
     227             :       ! the symmetry operation.
     228        1572 :       DO i = 1,3
     229        1179 :          dvec(i,nPlanes)=-xvec(i)
     230        5109 :          DO j=1,3
     231        4716 :             dvec(i,nPlanes) = dvec(i,nPlanes) + idrot(i,j,iSym)*xvec(j)
     232             :          END DO
     233             :       END DO
     234             : 
     235             :       ! Normalize the vector (it is normal to the possible boundary plane).
     236             :       ! Note that the vector points away from xvec.
     237             :       norm = 0.0
     238        2751 :       DO i = 1,3
     239        1572 :          norm = norm + dvec(i,nPlanes)**2
     240             :       END DO
     241         393 :       norm = SQRT(norm)
     242         393 :       dvec(:,nPlanes)=dvec(:,nPlanes) / norm
     243             : 
     244             :       ! Remove plane again if it is a duplicate of another plane
     245         393 :       foundDuplicate = .FALSE.
     246        5133 :       DO iPlane = nOuterPlanes+1, nPlanes-1
     247        4740 :          IF (ABS(ddist(iPlane)-ddist(nPlanes)).GT.eps09) CYCLE
     248        4740 :          IF (ABS(dvec(1,iPlane)-dvec(1,nPlanes)).GT.eps09) CYCLE
     249         148 :          IF (ABS(dvec(2,iPlane)-dvec(2,nPlanes)).GT.eps09) CYCLE
     250           6 :          IF (ABS(dvec(3,iPlane)-dvec(3,nPlanes)).GT.eps09) CYCLE
     251             :          foundDuplicate = .TRUE.
     252        5133 :          EXIT
     253             :       END DO
     254         419 :       IF(foundDuplicate) nPlanes = nPlanes - 1
     255             :    END DO
     256             : 
     257             :    ! 2. Determine all corner points of the IBZ.
     258             :    ! 2.1. Determine all possibly relevant intersection points of each 
     259             :    !      subset of 3 boundary planes, respectively.
     260             : 
     261             : 
     262             :    nCorners = 0
     263       28277 :    nPlaneCorners = 0
     264        1134 :    DO n1 = 1, nOuterPlanes
     265       45296 :       DO n2 = n1+1, nPlanes
     266             :          ! 2.1.1. Calculate intersection edge between planes n1 and n2
     267             :          !        and only consider those n1 and n2 where the cuttig 
     268             :          !        edge is not too far from the origin. (This would mean
     269             :          !        that the crossing point of these two planes and a third 
     270             :          !        plane is never relevant for the construction of the
     271             :          !        IBZ.)
     272             : 
     273             :          ! The direction of the intersection edge is the cross product
     274             :          ! of the normals on planes n1 and n2:
     275       44188 :          edgeDirec(1) = dvec(2,n1)*dvec(3,n2) - dvec(3,n1)*dvec(2,n2)
     276       44188 :          edgeDirec(2) = dvec(3,n1)*dvec(1,n2) - dvec(1,n1)*dvec(3,n2)
     277       44188 :          edgeDirec(3) = dvec(1,n1)*dvec(2,n2) - dvec(2,n1)*dvec(1,n2)
     278             : 
     279             :          ! Ignore parallel planes
     280             :          IF ((ABS(edgeDirec(1)).LT.eps09).AND.&
     281       44188 :              (ABS(edgeDirec(2)).LT.eps09).AND.&
     282             :              (ABS(edgeDirec(3)).LT.eps09)) CYCLE
     283             : 
     284             :          ! The distance vector of the intersection edge to the origin is given
     285             :          ! by (since dvec is normalized):
     286             :          !
     287             :          !       ddist(n1) - ddist(n2) * <dvec(:,n1)|dvec(:,n2)>              ddist(n2) - ddist(n1) * <dvec(:,n1)|dvec(:,n2)>
     288             :          ! d_e = ----------------------------------------------- dvec(:,n1) + ----------------------------------------------- dvec(:,n2)
     289             :          !              1.0 - <dvec(:,n1)|dvec(:,n2)>^2                               1.0 - <dvec(:,n1)|dvec(:,n2)>^2
     290             : 
     291             :          normalScalarProduct = dvec(1,n1)*dvec(1,n2)+&
     292             :                                dvec(2,n1)*dvec(2,n2)+&
     293       43418 :                                dvec(3,n1)*dvec(3,n2)
     294       43418 :          denominator = 1.0 - (normalScalarProduct**2.0)
     295       43418 :          vec1Fac = ddist(n1) - ddist(n2) * normalScalarProduct
     296       43418 :          vec2Fac = ddist(n2) - ddist(n1) * normalScalarProduct
     297      173672 :          edgeDistVec(:) = vec1Fac*dvec(:,n1) + vec2Fac*dvec(:,n2)
     298      173672 :          edgeDistVec(:) = edgeDistVec(:) / denominator
     299       43418 :          edgeDist = SQRT(edgeDistVec(1)**2.0+edgeDistVec(2)**2.0+edgeDistVec(3)**2.0)
     300             : 
     301             :          ! Ignore planes if intersection edge is too far from origin
     302             :          ! (is this criterion ok?)
     303       43418 :          IF (edgeDist.GT.maxRecDist) CYCLE
     304             : 
     305             :          ! Now calculate the possibly relevant crossing points
     306      389134 :          innerPlaneLoop: DO n3 = n2+1, nPlanes
     307             :             ! Set up system of linear equations to determine the crossing point
     308     2605246 :             DO i = 1, 3
     309     1116534 :                planesMatrix(1,i) = dvec(i,n1)
     310     1116534 :                planesMatrix(2,i) = dvec(i,n2)
     311     1488712 :                planesMatrix(3,i) = dvec(i,n3)
     312             :             END DO
     313      372178 :             solutions(1,1) = ddist(n1)
     314      372178 :             solutions(2,1) = ddist(n2)
     315      372178 :             solutions(3,1) = ddist(n3)
     316             :             ! Solve system of linear equations and cycle if no solution is found
     317             :             ! or solution is not in relevant region.
     318      372178 :             ipiv = 0
     319      372178 :             info = 0
     320      372178 :             CALL DGETRF(3,3, planesMatrix,3,ipiv,info)
     321      372178 :             IF(info.NE.0) CYCLE
     322      346638 :             CALL DGETRS('N',3,1,planesMatrix,3,ipiv,solutions,3,info)
     323      346638 :             IF(info.NE.0) CYCLE
     324      346638 :             vecDist = SQRT(solutions(1,1)**2.0+solutions(2,1)**2.0+solutions(3,1)**2.0)
     325      346638 :             IF(vecDist.GT.maxRecDist) CYCLE
     326             : 
     327             :    ! 2.2. Select those points that are on the correct side of all of the
     328             :    !      possible boundary planes.
     329             : 
     330     4595287 :             DO n4 = 1, nPlanes
     331     2335210 :                vecDist = dvec(1,n4)*solutions(1,1)+dvec(2,n4)*solutions(2,1)+dvec(3,n4)*solutions(3,1) - ddist(n4)
     332     2337135 :                IF(vecDist.GT.eps09) THEN
     333             :                   CYCLE innerPlaneLoop
     334             :                END IF
     335             :             END DO
     336             : 
     337             :             !Add new crossing point to list and associate it with the 3 planes
     338             : 
     339        1925 :             nCorners = nCorners + 1
     340        1925 :             corners(:,nCorners) = solutions(:,1)
     341        1925 :             cornerPlanes(1,nCorners) = n1
     342        1925 :             cornerPlanes(2,nCorners) = n2
     343        1925 :             cornerPlanes(3,nCorners) = n3
     344        1925 :             nPlaneCorners(n1) = nPlaneCorners(n1) + 1
     345        1925 :             nPlaneCorners(n2) = nPlaneCorners(n2) + 1
     346      415596 :             nPlaneCorners(n3) = nPlaneCorners(n3) + 1
     347             :          END DO innerPlaneLoop
     348             :       END DO
     349             :    END DO
     350             : 
     351             :    ! 3. Remove all boundary planes that don't feature at least three 
     352             :    !    corner points. These planes are not part of the real boundary.
     353             :    !    Also remove the corner points associated with each removed plane.
     354             :    !    ...and also those planes and corners that are irrelevant due to 
     355             :    !    other reasons.
     356             : 
     357             :    ! Corners might be present multiple times (for different triples of planes). 
     358             :    ! Remove duplicates. Collect planes intersecting at the respective point.
     359          26 :    maxNumPlanes = MAX(nPlanes-nOuterPlanes,3*nCorners) ! Just a save guess for the array dimension in the next line
     360          26 :    ALLOCATE(cornerPlaneList(maxNumPlanes,nCorners+1)) ! The +1 is because of the origin that is added later.
     361        1977 :    cornerPlaneList = 0
     362          26 :    numCornerPlanes = 0
     363             : 
     364          26 :    nUniqueCorners = 0
     365        1951 :    DO i = 1, nCorners
     366             :       duplicateNum = -1
     367       17159 :       compareLoop: DO j = 1, i-1
     368        9368 :          distVec(:) = corners(:,i) - corners(:,j)
     369        9368 :          vecDist = SQRT(distVec(1)**2+distVec(2)**2+distVec(3)**2)
     370        9542 :          IF (vecDist.LT.eps09) THEN
     371             :             duplicateNum = j
     372             :             EXIT compareLoop
     373             :          END IF 
     374             :       END DO compareLoop
     375        1951 :       IF (duplicateNum.EQ.-1) THEN
     376         174 :          nUniqueCorners = nUniqueCorners + 1
     377         174 :          numCornerPlanes(nUniqueCorners) = 3
     378         174 :          corners(:,nUniqueCorners) = corners(:,i)
     379         174 :          cornerPlaneList(1,nUniqueCorners) = cornerPlanes(1,i)
     380         174 :          cornerPlaneList(2,nUniqueCorners) = cornerPlanes(2,i)
     381         174 :          cornerPlaneList(3,nUniqueCorners) = cornerPlanes(3,i)
     382             :       ELSE
     383       12257 :          DO j = 1, 3
     384        5253 :             foundDuplicate = .FALSE.
     385       40150 :             DO k = 1, numCornerPlanes(duplicateNum)
     386       40150 :                IF (cornerPlaneList(k,duplicateNum).EQ.cornerPlanes(j,i)) THEN
     387        4899 :                   foundDuplicate = .TRUE.
     388             :                END IF
     389             :             END DO
     390        7004 :             IF (.NOT.foundDuplicate) THEN
     391         354 :                numCornerPlanes(duplicateNum) = numCornerPlanes(duplicateNum) + 1
     392         354 :                cornerPlaneList(numCornerPlanes(duplicateNum),duplicateNum) = cornerPlanes(j,i)
     393             :             END IF
     394             :          END DO
     395             :       END IF
     396             :    END DO
     397          26 :    nCorners = nUniqueCorners
     398             : 
     399             :    ! Add origin to corner points
     400          26 :    IF ((nPlanes-nOuterPlanes).GE.3) THEN ! The origin is only a corner if at least 3 planes feature this point
     401          24 :       nCorners = nCorners + 1
     402          24 :       corners(:,nCorners) = 0.0
     403         416 :       DO i = nOuterPlanes + 1, nPlanes
     404         416 :          cornerPlaneList(i-nOuterPlanes,nCorners) = i
     405             :       END DO
     406          24 :       numCornerPlanes(nCorners) = nPlanes - nOuterPlanes
     407             :    END IF
     408             : 
     409             :    ! Filter out "corners" found for sets of planes that do not meet in a single
     410             :    ! point but have a common intersection edge.
     411             :    ! (This absurd case actually occured so it has to be treated. LAPACK might
     412             :    ! find a "corner point" in this situation. Also the origin can be such a point.)
     413          26 :    nUniqueCorners = 0
     414         224 :    DO i = 1, nCorners
     415         198 :       filterOut = .FALSE.
     416         198 :       IF(numCornerPlanes(i).GE.3) THEN
     417         198 :          filterOut = .TRUE.
     418         198 :          vecA(:) = dvec(:,cornerPlaneList(1,i))
     419         198 :          vecB(:) = dvec(:,cornerPlaneList(2,i))
     420         366 :          cornerPlaneLoop: DO n3 = 3, numCornerPlanes(i)
     421         280 :             vecC(:) = dvec(:,cornerPlaneList(n3,i))
     422         280 :             testVec(1) = vecA(2)*vecB(3)-vecA(3)*vecB(2)
     423         280 :             testVec(2) = vecA(3)*vecB(1)-vecA(1)*vecB(3)
     424         280 :             testVec(3) = vecA(1)*vecB(2)-vecA(2)*vecB(1)
     425         280 :             scalarProduct = testVec(1)*vecC(1) + testVec(2)*vecC(2) + testVec(3)*vecC(3)
     426         282 :             IF (ABS(scalarProduct).GT.eps09) THEN
     427             :                filterOut = .FALSE.
     428             :                EXIT cornerPlaneLoop
     429             :             END IF
     430             :          END DO cornerPlaneLoop
     431             :       END IF
     432         198 :       IF(filterOut) THEN
     433             :          CYCLE
     434             :       END IF
     435         196 :       nUniqueCorners = nUniqueCorners + 1
     436         222 :       IF(nUniqueCorners.NE.i) THEN
     437           0 :          numCornerPlanes(nUniqueCorners) = numCornerPlanes(i)
     438           0 :          DO j = 1, numCornerPlanes(i)
     439           0 :             cornerPlaneList(j,nUniqueCorners) = cornerPlaneList(j,i)
     440             :          END DO
     441           0 :          corners(:,nUniqueCorners) = corners(:,i)
     442             :       END IF
     443             :    END DO
     444             :    nCorners = nUniqueCorners
     445             : 
     446             :    ! Count the number of corners for each plane
     447       28277 :    nPlaneCorners = 0
     448          26 :    ALLOCATE(planeCorners(nPlanes,nCorners))
     449         222 :    planeCorners = 0
     450         418 :    DO i = 1, nCorners
     451        1476 :       DO j = 1, numCornerPlanes(i)
     452        1254 :          nPlaneCorners(cornerPlaneList(j,i)) = nPlaneCorners(cornerPlaneList(j,i)) + 1
     453        1450 :          planeCorners(cornerPlaneList(j,i),nPlaneCorners(cornerPlaneList(j,i))) = i
     454             :       END DO
     455             :    END DO
     456             : 
     457             :    ! Remove irrelevant planes:
     458          26 :    nface = 0
     459       28277 :    isIBZPlane(:) = .TRUE.
     460        3028 :    DO n1 = 1, nPlanes
     461        1501 :       IF (nPlaneCorners(n1).LE.2) THEN
     462        1351 :          isIBZPlane(n1) = .FALSE.
     463        1351 :          CYCLE
     464             :       END IF
     465             : !      WRITE(*,*) 'plane ', n1
     466             : !      WRITE(*,'(4f20.13)') dvec(:,n1), ddist(n1)
     467             : !      WRITE(*,*) 'corners:'
     468             : !      DO i = 1, nPlaneCorners(n1)
     469             : !         WRITE(*,'(i5,3f20.13)') planeCorners(n1,i), corners(:,planeCorners(n1,i))
     470             : !      END DO
     471         176 :       nface = nface + 1
     472             :    END DO
     473             :    
     474             :    ! Remove irrelevant corners:
     475          26 :    ncorn = 0
     476      130026 :    isIBZCorner(:) = .TRUE.
     477             : !   WRITE(*,*) 'IBZ corners:'
     478         418 :    DO i = 1, nCorners
     479         196 :       numIBZPlanes = 0
     480        1450 :       DO j = 1, numCornerPlanes(i)
     481        1450 :          IF(isIBZPlane(cornerPlaneList(j,i))) THEN
     482         588 :             numIBZPlanes = numIBZPlanes + 1
     483             :          END IF
     484             :       END DO
     485         196 :       IF(numIBZPlanes.LE.2) isIBZCorner(i) = .FALSE.
     486         196 :       IF(.NOT.isIBZCorner(i)) CYCLE
     487          26 :       ncorn = ncorn + 1
     488             : !      WRITE(*,'(i5,3f20.13)') i, corners(:,i)
     489             :    END DO
     490             : 
     491          26 :    DEALLOCATE(cornerPlaneList)
     492             : 
     493             :    ! 4. Construct edges from the set of boundary points. Actually
     494             :    !    only the number of edges is needed.
     495             : 
     496             :    ! Count number of edges based on number of corners for each face
     497             :    ! Each face has the same number of edges as it has corners.
     498             :    ! Counting the number of edges in this way counts each edge twice.
     499          26 :    nedge = 0
     500        1527 :    DO n1 = 1, nPlanes
     501        1501 :       IF(.NOT.isIBZPlane(n1)) CYCLE
     502        1377 :       nedge = nedge + nPlaneCorners(n1)
     503             :    END DO
     504          26 :    nedge = nedge / 2
     505             : 
     506             :    ! 5. Fill the output arrays.
     507             : 
     508          26 :    IF((mface.LT.ncorn).OR.(mface.LT.nface)) THEN
     509           0 :       WRITE(*,*) "mface: ", mface
     510           0 :       WRITE(*,*) "ncorn: ", ncorn
     511           0 :       WRITE(*,*) "nface: ", nface
     512           0 :       WRITE(*,*) "mface has to be larger or equal to nface and ncorn."
     513           0 :       WRITE(*,*) "Its value is hardcoded. Adapt it."
     514           0 :       CALL juDFT_error("mface was chosen too small.",calledby ="brzone2")
     515             :    END IF
     516             : 
     517        1352 :    cpoint = 0.0
     518             :    j = 0
     519         418 :    DO i = 1, nCorners
     520         196 :       IF (.NOT.isIBZCorner(i)) CYCLE
     521         196 :       j = j + 1
     522         222 :       cpoint(:,j) = corners(:,i)
     523             :    END DO
     524             : 
     525        1352 :    fdist = 0.0
     526        1352 :    fnorm = 0.0
     527             :    j = 0
     528        3028 :    DO i = 1, nPlanes
     529        1501 :       IF(.NOT.isIBZPlane(i)) CYCLE
     530         150 :       j = j + 1
     531         150 :       fdist(j) = ddist(i)
     532         176 :       fnorm(:,j) = dvec(:,i)
     533             :    END DO
     534             : 
     535             : !   WRITE(*,*) 'ncorn', ncorn
     536             : !   WRITE(*,*) 'nedge', nedge
     537             : !   WRITE(*,*) 'nface', nface
     538             : !   WRITE(*,*) 'faces:'
     539             : !   DO i  =1,nface
     540             : !      WRITE(*,'(4f20.13)') fnorm(:,i), fdist(i)
     541             : !   END DO
     542             : !   WRITE(*,*) 'coners:'
     543             : !   DO i = 1,ncorn
     544             : !      WRITE(*,'(3f20.13)') cpoint(:,i)
     545             : !   END DO
     546             : 
     547             :    ! 5.1. Check Euler characteristic
     548             : 
     549          26 :    IF ((ncorn + nface - nedge).NE.2) THEN
     550           0 :       WRITE(*,*) "ncorn: ", ncorn
     551           0 :       WRITE(*,*) "nface: ", nface
     552           0 :       WRITE(*,*) "nedge: ", nedge
     553           0 :       WRITE(*,*) "corners: "
     554           0 :       DO i = 1, ncorn
     555           0 :          WRITE(*,'(3f20.13)') cpoint(:,i)
     556             :       END DO
     557           0 :       WRITE(*,*) "faces: "
     558           0 :       DO i = 1, nface
     559           0 :          WRITE(*,'(4f20.13)') fnorm(:,i), fdist(i)
     560             :       END DO
     561           0 :       CALL juDFT_error("Brillouin zone does not fulfill Euler characterisic.",calledby ="brzone2")
     562             :    END IF
     563             : 
     564          26 :    DEALLOCATE (planeCorners,dvec,ddist,nPlaneCorners,isIBZPlane)
     565             : 
     566          26 : END SUBROUTINE brzone2
     567             : 
     568             : END MODULE m_brzone2

Generated by: LCOV version 1.13