From 37b60e49d943bf8b1405a15ed5bfad8d9768bd54 Mon Sep 17 00:00:00 2001
From: Gregor Michalicek <g.michalicek@fz-juelich.de>
Date: Thu, 9 Feb 2017 11:27:42 +0100
Subject: [PATCH] Fixed bug in init/brzone2.f90

The problem was that in certain situations the origin was identified
as a corner point of the IBZ while it actually was not. It is not a
corner point if less than 3 faces of the IBZ feature this point.
---
 init/brzone2.f90 | 36 ++++++++++++++++++++++++++++++------
 1 file changed, 30 insertions(+), 6 deletions(-)

diff --git a/init/brzone2.f90 b/init/brzone2.f90
index d53abbc0..78902249 100644
--- a/init/brzone2.f90
+++ b/init/brzone2.f90
@@ -54,6 +54,7 @@ SUBROUTINE brzone2(rcmt,nsym,idrot,mface,nbsz,nv48,&
    REAL                  :: edgeDirec(3), edgeDistVec(3), distVec(3)
 
    INTEGER, ALLOCATABLE  :: cornerPlaneList(:,:)
+   INTEGER, ALLOCATABLE  :: planeCorners(:,:)
 
    ! Subroutine plan:
    ! 0. Initializations
@@ -387,12 +388,14 @@ SUBROUTINE brzone2(rcmt,nsym,idrot,mface,nbsz,nv48,&
    nCorners = nUniqueCorners
 
    ! Add origin to corner points
-   nCorners = nCorners + 1
-   corners(:,nCorners) = 0.0
-   DO i = nOuterPlanes + 1, nPlanes
-      cornerPlaneList(i-nOuterPlanes,nCorners) = i
-   END DO
-   numCornerPlanes(nCorners) = nPlanes - nOuterPlanes
+   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
 
    ! Filter out "corners" found for sets of planes that do not meet in a single
    ! point but have a common intersection edge.
@@ -433,9 +436,12 @@ SUBROUTINE brzone2(rcmt,nsym,idrot,mface,nbsz,nv48,&
 
    ! Count the number of corners for each plane
    nPlaneCorners = 0
+   ALLOCATE(planeCorners(nPlanes,nCorners))
+   planeCorners = 0
    DO i = 1, nCorners
       DO j = 1, numCornerPlanes(i)
          nPlaneCorners(cornerPlaneList(j,i)) = nPlaneCorners(cornerPlaneList(j,i)) + 1
+         planeCorners(cornerPlaneList(j,i),nPlaneCorners(cornerPlaneList(j,i))) = i
       END DO
    END DO
 
@@ -447,12 +453,19 @@ SUBROUTINE brzone2(rcmt,nsym,idrot,mface,nbsz,nv48,&
          isIBZPlane(n1) = .FALSE.
          CYCLE
       END IF
+!      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
       nface = nface + 1
    END DO
    
    ! Remove irrelevant corners:
    ncorn = 0
    isIBZCorner(:) = .TRUE.
+!   WRITE(*,*) 'IBZ corners:'
    DO i = 1, nCorners
       numIBZPlanes = 0
       DO j = 1, numCornerPlanes(i)
@@ -463,6 +476,7 @@ SUBROUTINE brzone2(rcmt,nsym,idrot,mface,nbsz,nv48,&
       IF(numIBZPlanes.LE.2) isIBZCorner(i) = .FALSE.
       IF(.NOT.isIBZCorner(i)) CYCLE
       ncorn = ncorn + 1
+!      WRITE(*,'(i5,3f20.13)') i, corners(:,i)
    END DO
 
    DEALLOCATE(cornerPlaneList)
@@ -527,9 +541,19 @@ SUBROUTINE brzone2(rcmt,nsym,idrot,mface,nbsz,nv48,&
       WRITE(*,*) "ncorn: ", ncorn
       WRITE(*,*) "nface: ", nface
       WRITE(*,*) "nedge: ", nedge
+      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
       CALL juDFT_error("Brillouin zone does not fulfill Euler characterisic.",calledby ="brzone2")
    END IF
 
+   DEALLOCATE (planeCorners)
+
 END SUBROUTINE brzone2
 
 END MODULE m_brzone2
-- 
GitLab