LCOV - code coverage report
Current view: top level - kpoints - brzone.f (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 318 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 1 0.0 %

          Line data    Source code
       1             :       MODULE m_brzone
       2             :       use m_juDFT
       3             : !
       4             : ! This subroutine finds the corner-points, the edges, and the
       5             : ! faces of the irreducible wedge of the brillouin zone (IBZ).
       6             : !
       7             :       CONTAINS
       8           0 :       SUBROUTINE brzone(
       9             :      >                   rcmt,nsym,idrot,mface,nbsz,nv48,
      10           0 :      =                   cpoint,
      11           0 :      <                   xvec,ncorn,nedge,nface,fnorm,fdist)
      12             : 
      13             :       USE m_constants, ONLY : pimach
      14             :       IMPLICIT NONE
      15             : 
      16             :       INTEGER, PARAMETER :: ibfile = 42
      17             : 
      18             :       INTEGER, INTENT (IN) :: mface,nbsz,nv48
      19             :       INTEGER, INTENT (IN) :: nsym               ! number of symmetry elements
      20             :       REAL,    INTENT (IN) :: rcmt(3,3)          ! reciprocal lattice basis (2\pi/a.u.)
      21             :       REAL,    INTENT (IN) :: idrot(3,3,48)      ! rotation matrices in cartesian repr.
      22             : 
      23             :       INTEGER, INTENT (OUT) :: ncorn,nedge,nface ! number of corners, faces and edges of the IBZ
      24             :       REAL,    INTENT (OUT) :: fnorm(3,mface)    ! normal vector of the planes bordering the IBZ
      25             :       REAL,    INTENT (OUT) :: fdist(mface)      ! distance vector of the planes bordering the IBZ
      26             :       REAL,    INTENT (OUT) :: cpoint(3,mface)   ! cartesian coordinates of corner points of IBZ
      27             :       REAL,    INTENT (OUT) :: xvec(3)           ! arbitrary vector lying in the IBZ
      28             : C
      29             : C   LOCAL variables
      30             : C
      31             :       REAL pi
      32             :       REAL scale,sum,amin,alpha
      33             :       REAL bmin,beta,cmin,cmax,gamma
      34             :       REAL sx,xmin
      35             :       INTEGER ntl,krecip(3),ntot,ip,i,j,l,n,m,nfp
      36             :       INTEGER nmin,mmin,lmin,lmax,nf,ncf
      37             :       INTEGER n1,n2,n3,nn,k,ii
      38             : C
      39             : C     Local working arrays and pointers
      40             : C
      41           0 :       REAL epoint(3,mface),fpoint(3,mface),cstart(3,2,mface)
      42             :       REAL fvec(3),evec(3),dir(3),c0(3),c1(3),c2(3),csum(3)
      43           0 :       REAL sk(3),yvec(3),ddist(nv48),dvec(3,nv48)
      44           0 :       INTEGER nplane(mface)
      45             : C
      46             : C----->  Intrinsic Functions
      47             : C
      48             :       INTRINSIC min,sqrt
      49             : C
      50           0 :       OPEN (ibfile,form='formatted',status='scratch')
      51             : c
      52           0 :       WRITE (ibfile,'('' brzone '')')
      53           0 :       ntot = (2*nbsz + 1)**3
      54           0 :       WRITE (ibfile,'('' ntot = '',i4,'' nsym = '',i4)') ntot,nsym
      55           0 :       ntl = ntot + nsym - 2
      56           0 :       WRITE (ibfile,'('' ntl = '',i4)') ntl
      57           0 :       WRITE (ibfile,'('' rcmt '',/)')
      58             : c     WRITE (ibfile,*) rcmt
      59           0 :       WRITE (ibfile,101) ((rcmt(i,j),j=1,3),i=1,3)
      60             :  101  FORMAT(/5x,3(f10.6,3x),2(/5x,3(f10.6,3x)))
      61             : C
      62             : C construct all boundary-planes
      63             : C first the planes that determine the first brillouin zone
      64             : C that is, the planes bisecting the line connecting the
      65             : C origin with a reciprocal lattice vector ( <> 0 )
      66             : C
      67           0 :       pi = pimach()
      68           0 :       DO i = 1,3
      69           0 :        sk(i) = 0.0
      70           0 :         DO j = 1,3
      71           0 :          sk(i)=sk(i)+rcmt(j,i)*rcmt(j,i)
      72             :         ENDDO
      73             :       ENDDO
      74           0 :       WRITE (ibfile,'('' sk(1) = j=1,3 of rcmt(j,1)**2 '')')
      75           0 :       WRITE (ibfile,97) (sk (ii),ii=1,3)
      76             :  97   FORMAT (/5x,'  sk(i) ',3(f13.6,2x))
      77             : 
      78           0 :       scale = sqrt(min(sk(1),sk(2),sk(3)))*0.1
      79           0 :       xvec(1) = scale
      80           0 :       xvec(2) = scale/sqrt(pi)
      81           0 :       xvec(3) = scale/pi
      82           0 :       WRITE (ibfile,98) (xvec(ii),ii=1,3)
      83             :  98   FORMAT (/5x,' xvec(i) ',3(f13.6,2x))
      84             : 
      85           0 :       n = 0
      86           0 :       DO n1 = -nbsz,nbsz
      87           0 :         krecip(1) = n1
      88           0 :         DO n2 = -nbsz,nbsz
      89           0 :         krecip(2) = n2
      90           0 :           DO n3 = -nbsz,nbsz
      91           0 :             IF ( .NOT.(n1.EQ.0.AND.n2.EQ.0.AND.n3.EQ.0) ) THEN
      92           0 :             krecip(3) = n3
      93             : 
      94           0 :             n = n + 1
      95           0 :             DO i = 1,3
      96           0 :               dvec(i,n) = 0.0
      97           0 :               DO j = 1,3
      98           0 :                 dvec(i,n) = dvec(i,n) + rcmt(i,j)*krecip(j)
      99             :               ENDDO
     100             :             ENDDO
     101           0 :             WRITE (ibfile,99) n,(dvec(k,n),k=1,3)
     102             :  99         FORMAT(/5x,'  dvec(k,',i4,') ',3(f13.6,2x))
     103             : 
     104           0 :             sum = 0.0
     105           0 :             DO i = 1,3
     106           0 :               sum = sum + dvec(i,n)**2
     107           0 :               WRITE (ibfile,'('' sum = dvec**2 = '',f13.6)') sum
     108             :             ENDDO
     109           0 :             sum = sqrt(sum)
     110           0 :             ddist(n) = 0.5*sum
     111             :             WRITE (ibfile,'(/'' ddist('',i3,'')=(.5*sum**.5) '',f13.6)')
     112           0 :      >                                                        n,ddist(n)
     113           0 :             sum = 1.0/sum
     114           0 :             WRITE (ibfile,'(/'' sum = ( 1/(.5*sum**.5) )'',f13.6)') sum
     115           0 :             DO i = 1,3
     116           0 :               dvec(i,n) = dvec(i,n)*sum
     117             :             ENDDO
     118           0 :             WRITE (ibfile,'('' dvec(i,n) * latest sum '')')
     119           0 :             WRITE (ibfile,99) n,(dvec(k,n),k=1,3)
     120             : 
     121             :             ENDIF
     122             :           ENDDO
     123             :         ENDDO
     124             :       ENDDO
     125             : 
     126             : C
     127             : C construct the planes that determine the irreducible wedge
     128             : C that is, the planes bisecting the line connecting xvec
     129             : C with an element of the star of xvec ( <> xvec )
     130             : C
     131           0 :       WRITE (ibfile,'('' working on star of xvec '')')
     132           0 :       WRITE (ibfile,'('' ntot = '',i4,'' ntl = '',i4,/)') ntot,ntl
     133           0 :       DO n = ntot,ntl
     134           0 :         ddist(n) = 0.0
     135           0 :         WRITE (ibfile,'(/)')
     136             : 
     137           0 :         DO i = 1,3
     138           0 :           dvec(i,n)=-xvec(i)
     139             :           WRITE (ibfile,'('' dvec('',i3,i4,'')=(here-xvec('',i2,'') '',
     140           0 :      +                                         f13.6)') i,n,i,dvec(i,n)
     141           0 :           DO j=1,3
     142           0 :             dvec(i,n) = dvec(i,n) + idrot(i,j,n+2-ntot)*xvec(j)
     143             :             WRITE (ibfile,'('' idrot('',i3,i3,i4,'') = '',f10.6)') 
     144           0 :      +                            i,j,n+2-ntot,idrot(i,j,n+2-ntot)
     145           0 :             WRITE (ibfile,'('' xvec('',i3,'') = '',f6.4)') j,xvec(j)
     146             :             WRITE (ibfile,'('' dvec('',i3,i4,'') = '',f13.6,/)') 
     147           0 :      +                                             i,n,dvec(i,n)
     148             :           ENDDO
     149             :         ENDDO
     150             : 
     151           0 :         sum = 0.0
     152           0 :         DO i = 1,3
     153           0 :           sum =sum + dvec(i,n)**2
     154           0 :           WRITE (ibfile,'('' sum = dvec**2 = '',f13.6)') sum
     155             :         ENDDO
     156           0 :         sum = 1.0/sqrt(sum)
     157           0 :         WRITE (ibfile,'(/'' sum = ( 1/(sum**.5) )'',f13.6)') sum
     158           0 :         DO i = 1,3
     159           0 :             dvec(i,n)=dvec(i,n)*sum
     160             :         ENDDO
     161           0 :         WRITE (ibfile,'('' dvec(i,n) * latest sum '')')
     162           0 :         WRITE (ibfile,99) n,(dvec(k,n),k=1,3)
     163             :       ENDDO
     164           0 :       nn = ntl - ntot + 1
     165             : C
     166             : C find the point on the line determined by the origin and xvec
     167             : C which is on the nearest boundary plane
     168             : C
     169           0 :       WRITE (ibfile,'(/,'' find points on nearest boundary plane '')')
     170           0 :       amin = scale*99999.9
     171           0 :       nmin = 0
     172           0 :       DO n = 1,ntl
     173           0 :         sum = 0.0
     174           0 :         DO i = 1,3
     175           0 :           sum = sum + xvec(i)*dvec(i,n)
     176             :         ENDDO
     177           0 :         WRITE (ibfile,'('' sum('',i4,'') = '',f13.6)') n,sum
     178           0 :         IF ( abs(sum).GT.1.0e-10 ) THEN
     179           0 :           alpha=ddist(n)/sum
     180           0 :           WRITE (ibfile,'('' alpha('',i4,'') = '',f13.6)') n,alpha
     181           0 :           IF ( .NOT.((alpha.LE.0.0).OR.(alpha.GT.amin)) ) THEN 
     182           0 :             amin = alpha
     183           0 :             nmin = n
     184           0 :             WRITE (ibfile,'('' nmin = '',i4)') n
     185             :           ENDIF
     186             :         ENDIF
     187             :       ENDDO
     188           0 :       IF ( nmin==0 )  CALL juDFT_error("bzone1",calledby ="brzone")
     189           0 :       WRITE (ibfile,'('' amin = '',f13.6)') amin
     190           0 :       DO i = 1,3
     191           0 :         fvec(i) = amin*xvec(i)
     192             :       ENDDO
     193           0 :       WRITE (ibfile,'('' fvec('',i3,'') = '',f13.6)') (i,fvec(i),i=1,3)
     194           0 :       nplane(1) = nmin
     195             : C
     196             : C find the nearest edge in this plane, along the line connecting
     197             : C fvec and the center of the plane, given by dvec*ddist
     198             : C
     199           0 :       WRITE (ibfile,'(/,'' find nearest edge in this plane '')')
     200           0 :       bmin = scale*99999.9
     201           0 :       mmin = 0
     202           0 :       DO m = 1,ntl
     203           0 :         IF ( m.NE.nmin ) THEN
     204           0 :           sum=0.0
     205           0 :           DO i = 1,3
     206           0 :             sum = sum + dvec(i,m)*(fvec(i)-dvec(i,nmin)*ddist(nmin))
     207             :           ENDDO
     208           0 :           WRITE (ibfile,'('' sum('',i4,'') = '',f13.6)') m,sum
     209           0 :           IF ( abs(sum).GT.1.0e-10 ) THEN
     210           0 :             beta = ddist(m)
     211           0 :             WRITE (ibfile,'('' beta('',i4,'') = '',f13.6)') m,ddist(m)
     212           0 :             DO i = 1,3
     213           0 :               beta=beta-fvec(i)*dvec(i,m)
     214             :             ENDDO
     215           0 :             WRITE (ibfile,'('' beta-fvec(i)*dvec(i,m) = '',f13.6)') beta
     216           0 :             beta = beta/sum
     217           0 :             IF ( .NOT.((beta.LT.0.0).OR.(beta.GT.bmin)) ) THEN
     218           0 :               bmin=beta
     219           0 :               mmin=m
     220             :             ENDIF
     221             :           ENDIF
     222             :         ENDIF
     223             :       ENDDO
     224           0 :       IF ( mmin==0 )  CALL juDFT_error("bzone2",calledby ="brzone")
     225           0 :       DO i = 1,3
     226           0 :         evec(i) = fvec(i) + bmin*(fvec(i)-dvec(i,nmin)*ddist(nmin))
     227             :       ENDDO
     228             :       WRITE (ibfile,'(/,'' evec('',i3,'') = '',f13.6)') 
     229           0 :      +                                (i,evec(i),i=1,3)
     230             : C
     231             : C find innermost boundary plane for this edge
     232             : C
     233           0 :       WRITE (ibfile,'(/,'' find innermost boundary plane '')')
     234           0 :       xmin = scale*99999.9
     235           0 :       mmin = 0
     236           0 :       DO  m = 1,ntl
     237           0 :         IF ( m.NE.nmin ) THEN
     238           0 :           sum = ddist(m)
     239           0 :           sx  = 0.0
     240           0 :           DO i = 1,3
     241           0 :             sum = sum - dvec(i,m)*evec(i)
     242           0 :             sx  = sx  + dvec(i,m)*dvec(i,nmin)
     243             :           ENDDO
     244           0 :           IF ( .NOT.((abs(sum).GT.1.e-10).OR.(sx.GT.xmin)) ) THEN
     245           0 :             xmin = sx
     246           0 :             mmin = m
     247             :           ENDIF
     248             :         ENDIF
     249             :       ENDDO
     250           0 :       IF ( mmin.EQ.0 )  CALL juDFT_error("bzone25",calledby="brzone")
     251           0 :       nplane(2) = mmin
     252             : C
     253             : C find direction of the edge
     254             : C
     255           0 :       dir(1) = dvec(2,nmin)*dvec(3,mmin) - dvec(3,nmin)*dvec(2,mmin)
     256           0 :       dir(2) = dvec(3,nmin)*dvec(1,mmin) - dvec(1,nmin)*dvec(3,mmin)
     257           0 :       dir(3) = dvec(1,nmin)*dvec(2,mmin) - dvec(2,nmin)*dvec(1,mmin)
     258           0 :       WRITE (ibfile,'('' dir('',i3,'') = '',f13.6)') (i,dir(i),i=1,3)
     259             : C
     260             : C find the corner points on this edge
     261             : C
     262           0 :       WRITE (ibfile,'('' find corner points on this edge '')')
     263           0 :       cmin = scale*99999.9
     264           0 :       cmax = -cmin
     265           0 :       lmin = 0
     266           0 :       lmax = 0
     267           0 :       DO l=1,ntl
     268           0 :         IF ( (l.EQ.nmin).OR.(l.EQ.mmin) ) GOTO 2700
     269           0 :         sum = 0.0
     270           0 :         DO i=1,3
     271           0 :           sum = sum + dir(i)*dvec(i,l)
     272             :         ENDDO
     273           0 :         IF ( abs(sum).LT.1.0e-10 ) GOTO 2700
     274           0 :         gamma=ddist(l)
     275           0 :         DO i = 1,3
     276           0 :           gamma = gamma - evec(i)*dvec(i,l)
     277             :         ENDDO
     278           0 :         gamma = gamma/sum
     279           0 :         IF ( gamma.GE.0.0 ) THEN
     280           0 :           IF (gamma.gt.cmin) GOTO 2700
     281           0 :           cmin=gamma
     282           0 :           lmin=l
     283           0 :           GOTO 2700
     284             :         ENDIF
     285           0 :         IF ( gamma.GE.cmax ) THEN
     286           0 :           cmax=gamma
     287           0 :           lmax=l
     288             :         ENDIF
     289           0 :  2700   CONTINUE
     290             :       ENDDO
     291           0 :       IF ( lmax*lmin.EQ.0 ) CALL juDFT_error("bzone3",calledby="brzone")
     292           0 :       WRITE (ibfile,'('' cmax = '',f13.6)') cmax
     293           0 :       WRITE (ibfile,'('' cmin = '',f13.6)') cmin
     294           0 :       DO i=1,3
     295           0 :          c0(i) = evec(i)+cmax*dir(i)
     296           0 :          WRITE (ibfile,'('' dir('',i3,'') = '',f13.6)') i,dir(i)
     297           0 :          WRITE (ibfile,'('' evec('',i3,'') = '',f13.6)') i,evec(i)
     298           0 :          WRITE (ibfile,'('' c0('',i3,'') = '',f13.6)') i,c0(i)
     299           0 :          c1(i)=evec(i)+cmin*dir(i)
     300           0 :          WRITE (ibfile,'('' c1('',i3,'') = '',f13.6)') i,c1(i)
     301             :       ENDDO
     302             : C
     303             : C prepare the list of corner points, etc, for the
     304             : C general scheme of finding the boundaries of the
     305             : C irreducible wedge of the first brillouin zone
     306             : C
     307           0 :       WRITE (ibfile,'(/,'' prepare list of corner points '')')
     308           0 :       DO i = 1,3
     309           0 :         cstart(i,1,1) = c0(i)
     310           0 :         cstart(i,2,1) = c1(i)
     311           0 :         cstart(i,1,2) = c1(i)
     312           0 :         cstart(i,2,2) = c0(i)
     313           0 :         cpoint(i,1)   = c0(i)
     314           0 :         cpoint(i,2)   = c1(i)
     315           0 :         epoint(i,1)   = 0.5*(c0(i)+c1(i))
     316           0 :         WRITE (ibfile,'('' cstart('',i2,'',1,1) = '',f13.6)') i,c0(i)
     317           0 :         WRITE (ibfile,'('' cstart('',i2,'',2,1) = '',f13.6)') i,c1(i)
     318           0 :         WRITE (ibfile,'('' cstart('',i2,'',1,2) = '',f13.6)') i,c1(i)
     319           0 :         WRITE (ibfile,'('' cstart('',i2,'',2,2) = '',f13.6)') i,c0(i)
     320           0 :         WRITE (ibfile,'('' cpoint('',i2,'',1) = '',f13.6)') i,c0(i)
     321           0 :         WRITE (ibfile,'('' cpoint('',i2,'',2) = '',f13.6)') i,c1(i)
     322           0 :         WRITE (ibfile,'('' epoint('',i2,'',1) = '',f13.6)')i,epoint(i,1)
     323             :       ENDDO
     324           0 :       ncorn = 2
     325           0 :       nedge = 1
     326           0 :       nface = 2
     327           0 :       nf = 0
     328             : C
     329             : C enter general loop which determines all corners and all edges
     330             : C of all faces , new faces are added to the list nplane
     331             : C
     332             :  4000 CONTINUE
     333           0 :       nf  = nf + 1
     334           0 :       nfp = nplane(nf)
     335             : C
     336             : C we consider face number nf
     337             : C start with the corner points of cstart , notice that the order
     338             : C of the corner points is important and is determined by the
     339             : C order in the outer product of the vectors dvec
     340             : C
     341           0 :       DO i=1,3
     342           0 :         c0(i)   = cstart(i,1,nf)
     343           0 :         c1(i)   = cstart(i,1,nf)
     344           0 :         c2(i)   = cstart(i,2,nf)
     345           0 :         csum(i) = cstart(i,1,nf) + cstart(i,2,nf)
     346             :       ENDDO
     347             :       ncf = 2
     348             :  4200 CONTINUE
     349             : C
     350             : C determine the point fvec
     351             : C
     352           0 :       fvec(1) = dvec(2,nfp)*(c2(3)-c1(3))-dvec(3,nfp)*(c2(2)-c1(2))
     353           0 :       fvec(2) = dvec(3,nfp)*(c2(1)-c1(1))-dvec(1,nfp)*(c2(3)-c1(3))
     354           0 :       fvec(3) = dvec(1,nfp)*(c2(2)-c1(2))-dvec(2,nfp)*(c2(1)-c1(1))
     355             : c     WRITE (ibfile,'('' pt fvec('',i3,'') = '',f13.6)') (i,fvec(i),i=1,3)
     356           0 :       DO i = 1,3
     357           0 :         fvec(i) = 0.5*(c2(i)+c1(i)) + 0.001*fvec(i)
     358             :       ENDDO
     359             : C
     360             : C determine the edge connected to c2 by moving outwards on c2-c1
     361             : C and finding the nearest intersection with a boundary plane
     362             : C on the line connecting this point and fvec , which is
     363             : C on the correct side of the line c2-c1 , by construction ,
     364             : C because of the way we order the corner points
     365             : C
     366           0 :       DO i = 1,3
     367           0 :         yvec(i) = c2(i) + 1.0e-5*(c2(i)-c1(i))
     368             :       ENDDO
     369             : C
     370             : C find nearest boundary plane
     371             : C
     372           0 :       bmin = scale*99999.9
     373           0 :       mmin = 0
     374           0 :       DO m = 1,ntl
     375           0 :        IF ( m.NE.nplane(nf) ) THEN
     376           0 :          sum = 0.0
     377           0 :          DO i = 1,3
     378           0 :            sum = sum + dvec(i,m)*(yvec(i)-fvec(i))
     379             :          ENDDO
     380           0 :          IF ( abs(sum).GE.1.0e-10 ) THEN
     381           0 :            beta=ddist(m)
     382           0 :            DO i = 1,3
     383           0 :              beta = beta - fvec(i)*dvec(i,m)
     384             :            ENDDO
     385           0 :            beta = beta/sum
     386           0 :            IF ( .NOT.((beta.LT.0.0).OR.(beta.GT.bmin)) ) THEN
     387           0 :              bmin = beta
     388           0 :              mmin = m
     389             :             ENDIF
     390             :           ENDIF
     391             :         ENDIF
     392             :       ENDDO
     393           0 :       IF ( mmin.EQ.0 )  CALL juDFT_error("bzone4",calledby="brzone")
     394             : C
     395             : C construct direction of this edge
     396             : C
     397           0 :       dir(1) = dvec(2,nfp)*dvec(3,mmin) - dvec(3,nfp)*dvec(2,mmin)
     398           0 :       dir(2) = dvec(3,nfp)*dvec(1,mmin) - dvec(1,nfp)*dvec(3,mmin)
     399           0 :       dir(3) = dvec(1,nfp)*dvec(2,mmin) - dvec(2,nfp)*dvec(1,mmin)
     400           0 :       WRITE (ibfile,'(''2 dir('',i3,'') = '',f13.6)') (i,dir(i),i=1,3)
     401             : C
     402             : C find other corner point on this edge
     403             : C
     404           0 :       cmin = scale*99999.9
     405           0 :       lmin = 0
     406           0 :       DO l = 1,ntl
     407           0 :         IF ( .NOT.((l.EQ.nplane(nf)).OR.(l.EQ.mmin)) ) THEN 
     408           0 :           sum = 0.0
     409           0 :           DO i = 1,3
     410           0 :             sum = sum + dir(i)*dvec(i,l)
     411             :           ENDDO
     412           0 :           IF ( abs(sum).GE.1.0e-10 )  THEN
     413           0 :             gamma=ddist(l)
     414           0 :             DO i = 1,3
     415           0 :               gamma = gamma - dvec(i,l)*c2(i)
     416             :             ENDDO
     417           0 :             gamma = gamma/sum
     418           0 :             IF ( .NOT.((gamma.LT.1.0e-9).OR.(gamma.GT.cmin)) ) THEN
     419           0 :               cmin = gamma
     420           0 :               lmin = l
     421             :             ENDIF
     422             :           ENDIF
     423             :         ENDIF
     424             :       ENDDO
     425           0 :       IF ( lmin.EQ.0 )  CALL juDFT_error("bzone5",calledby="brzone")
     426             : C
     427             : C move c2 and c1
     428             : C
     429           0 :       DO i = 1,3
     430           0 :         c1(i)   = c2(i)
     431           0 :         c2(i)   = c1(i) + cmin*dir(i)
     432           0 :         evec(i) = 0.5*( c1(i)+c2(i) )
     433             :       ENDDO
     434           0 :       WRITE (ibfile,'(''corner c1('',i3,'')='',f13.6)') (i,c1(i),i=1,3)
     435           0 :       WRITE (ibfile,'(''corner c2('',i3,'')='',f13.6)') (i,c2(i),i=1,3)
     436           0 :       WRITE (ibfile,'(''evec('',i3,'') = '',f13.6)') (i,evec(i),i=1,3)
     437             : C
     438             : C find innermost boundary plane for this edge
     439             : C
     440           0 :       xmin = scale*99999.9
     441           0 :       mmin = 0
     442           0 :       WRITE (ibfile,'(/,''bzone55 loop ntl='',i4,'' nfp='',i4)') ntl,nfp
     443           0 :       DO m = 1,ntl
     444           0 :         IF ( m.NE.nfp ) THEN
     445           0 :           sum = ddist(m)
     446           0 :           sx  = 0.0
     447           0 :           DO i=1,3
     448           0 :            sum = sum - dvec(i,m)*evec(i)
     449           0 :            sx  = sx  + dvec(i,m)*dvec(i,nfp)
     450             :           ENDDO
     451           0 :           IF ( .NOT.((abs(sum).GT.1.0e-6).OR.(sx.GT.xmin)) ) THEN
     452           0 :             xmin = sx
     453           0 :             mmin = m
     454             :             WRITE (ibfile,'('' m = '',i4,'' xmin = '',f16.12,'' nfp = ''
     455           0 :      +                                                ,i4)')  m,xmin,nfp
     456             :           ENDIF
     457             :         ENDIF
     458             :       ENDDO
     459             :       WRITE (ibfile,'('' m = '',i4,'' xmin = '',f16.12,'' nfp = '',i4)')
     460           0 :      +                                                        m,xmin,nfp
     461           0 :       IF ( mmin.EQ.0 )  CALL juDFT_error("bzone55",calledby="brzone")
     462             : C
     463             : C check if we have found a new face or not
     464             : C
     465           0 :       DO ip = 1,nface
     466           0 :          IF (nplane(ip).EQ.mmin) GOTO 5400
     467             :       ENDDO
     468           0 :       nface = nface + 1
     469           0 :       WRITE (ibfile,'('' nface = '',i4)') nface
     470           0 :       nplane(nface) = mmin
     471           0 :       DO i = 1,3
     472           0 :         cstart(i,1,nface) = c2(i)
     473           0 :         cstart(i,2,nface) = c1(i)
     474             :         WRITE (ibfile,'('' cstart('',i3,'', 1,'',i3,'') = '',f13.6)')
     475           0 :      +   i,nface,c2(i)
     476             :         WRITE (ibfile,'('' cstart('',i3,'', 2,'',i3,'') = '',f13.6)')
     477           0 :      +   i,nface,c1(i)
     478             :       ENDDO 
     479             :  5400 CONTINUE
     480             : C
     481             : C check if the new corner and edge points are contained
     482             : C in the list of existing points
     483             : C
     484           0 :       DO ip = 1,ncorn
     485           0 :         sum = 0.00
     486           0 :         DO i = 1,3
     487           0 :           sum = sum + (c2(i) - cpoint(i,ip))**2
     488             :         ENDDO
     489           0 :         IF ( abs(sum).LT.1.0e-10 ) GOTO 6300
     490             :       ENDDO 
     491           0 :       ncorn = ncorn + 1
     492           0 :       WRITE (ibfile,'('' ncorn = '',i5)') ncorn
     493           0 :       DO i = 1,3
     494           0 :         cpoint(i,ncorn) = c2(i)
     495             :       ENDDO
     496             :  6300 CONTINUE
     497             : c
     498           0 :       DO ip = 1,nedge
     499           0 :         sum = 0.0
     500           0 :         DO i = 1,3
     501           0 :           sum = sum + (evec(i) - epoint(i,ip))**2
     502             :         ENDDO 
     503           0 :       IF ( abs(sum).LT.1.0e-10 ) GOTO 6700
     504             :       ENDDO
     505           0 :       nedge = nedge + 1
     506           0 :       WRITE (ibfile,'('' nedge = '',i5)') nedge
     507           0 :       DO i = 1,3
     508           0 :        epoint(i,nedge) = evec(i)
     509             :       ENDDO
     510             :  6700 CONTINUE
     511             : C
     512             : C check if we have all points on this face
     513             : C
     514           0 :       sum = 0.0
     515           0 :       DO i = 1,3
     516           0 :        sum = sum + ( c2(i) - c0(i) )**2
     517             :       ENDDO
     518           0 :       IF ( abs(sum).GT.1.0e-10 ) THEN
     519           0 :         ncf = ncf + 1
     520           0 :         WRITE (ibfile,'('' nface = '',i4)') nface
     521           0 :         GOTO 4200
     522             :       ENDIF
     523             : C
     524             : C we have found all corner points on this face
     525             : C determine the center of gravity of this face
     526             : C
     527           0 :       DO i = 1,3
     528           0 :         fpoint(i,nf) = csum(i)/ncf
     529             :       ENDDO
     530           0 :       IF ( nf.LT.nface ) GOTO 4000
     531             : c
     532           0 :       DO ip  =1,nface
     533           0 :         nf = nplane(ip)
     534           0 :         fdist(ip) = ddist(nf)
     535           0 :         DO i=1,3
     536           0 :            fnorm(i,ip) = dvec(i,nf)
     537             :         ENDDO
     538             :       ENDDO
     539             : c
     540             : 
     541             : !      WRITE(*,*) 'ncorn', ncorn
     542             : !      WRITE(*,*) 'nedge', nedge
     543             : !      WRITE(*,*) 'nface', nface
     544             : !      WRITE(*,*) 'faces:'
     545             : !      DO ip  =1,nface
     546             : !         WRITE(*,'(4f20.13)') fnorm(:,ip), fdist(ip)
     547             : !      END DO
     548             : !      WRITE(*,*) 'coners:'
     549             : !      DO ip = 1,ncorn
     550             : !         WRITE(*,'(3f20.13)') cpoint(:,ip)
     551             : !      END DO
     552             : 
     553           0 :       WRITE (6,7100) ncorn,nedge,nface
     554           0 :       WRITE (ibfile,7100) ncorn,nedge,nface
     555             :  7100 FORMAT (///,'  the irreducible wedge of the first brillouin'
     556             :      $,' zone has :  ',/,
     557             :      $     i10,'     corner points   ',/,
     558             :      $     i10,'     edges           ',/,
     559             :      $     i10,'     faces           ')
     560           0 :       IF ( (ncorn + nface - nedge)/=2 )  CALL juDFT_error("bzone6"
     561           0 :      +     ,calledby ="brzone")
     562           0 :       WRITE (6,7200) ((cpoint(i,ip),i=1,3),ip=1,ncorn)
     563           0 :       WRITE (ibfile,7200) ((cpoint(i,ip),i=1,3),ip=1,ncorn)
     564             :  7200 FORMAT(//,'    corner points in carthesian units ',
     565             :      $     99(/,3f10.5))
     566             : 
     567           0 :       CLOSE (ibfile)
     568           0 :       RETURN
     569             :       END SUBROUTINE brzone
     570             :       END MODULE m_brzone

Generated by: LCOV version 1.13