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

          Line data    Source code
       1             :       MODULE m_tetcon
       2             :       use m_juDFT
       3             : !
       4             : ! This subroutine constructs the tetrahedra for the
       5             : ! Brillouin zone integration
       6             : !
       7             :       CONTAINS
       8           1 :       SUBROUTINE tetcon(
       9             :      >                  iofile,ibfile,mkpt,ndiv3,
      10           1 :      >                  nkpt,omega,kvc,nsym,
      11           2 :      <                  nt,voltet,ntetra)
      12             : 
      13             :       IMPLICIT NONE
      14             : c
      15             :       INTEGER, INTENT (IN) :: iofile,ibfile
      16             :       INTEGER, INTENT (IN) :: mkpt,ndiv3,nkpt,nsym
      17             :       REAL,    INTENT (IN) :: omega
      18             :       REAL,    INTENT (IN) :: kvc(3,mkpt)
      19             :       INTEGER, INTENT (OUT) :: nt
      20             :       INTEGER, INTENT (OUT) :: ntetra(4,ndiv3)
      21             :       REAL,    INTENT (OUT) :: voltet(ndiv3)
      22             : 
      23             :       INTEGER i,it,j,ic,icom,ik,nkp,jj,nnt,nsid,nkq
      24             :       INTEGER i2,i3,i4,nk2,k,l,kk,is1,is2,is3,js1,js2,js3
      25           2 :       INTEGER nside(4),nside2(4),nkadd(mkpt)
      26             :       REAL dnorm,pi,vav,vsq,sav,ssq,dm,dist,dl,vect
      27             :       REAL vol,vt,length,minlen,eps,eps1
      28             :       REAL vmin,vmax,smin,smax,volnew,sum,sss
      29           2 :       REAL ddist(mkpt),kcorn(3,4),norm(3),d(3,16),dnm(3),cn(3)
      30             : C
      31             : C----->  Intrinsic Functions
      32             : C
      33             :        INTRINSIC abs,max,min,sqrt
      34             : C
      35             : C*************************************************************************
      36             : C
      37             : C     All commons have been removed from the historical program and re-
      38             : C     placed by direct subroutine calls.  Those statements immediately
      39             : C     below preceded by stars your removed prior to my changes
      40             : C
      41             : C                                          June 9, 1992
      42             : C                                                   Fred Hutson
      43             : C
      44             : C*************************************************************************
      45             : C
      46             :       data eps/1e-10/
      47             :       data eps1/1e-5/
      48           1 :       pi = 4.0*atan(1.0)
      49             : C
      50             : C CONSTRUCT THE FIRST TETRAHEDRON
      51             : C
      52           1 :       ntetra(1,1)=1
      53             : C
      54             : C FIND K-POINT NEAREST TO 1
      55             : C
      56           1 :       dm=10.0**9
      57           1 :       i2=0
      58          20 :       do 100 i=2,nkpt
      59             :       dist=0.0
      60         133 :       do 110 icom=1,3
      61          57 :       dist=dist+(kvc(icom,1)-kvc(icom,i))**2
      62          19 :   110 continue
      63          19 :       if ( dist.gt.dm ) goto 100
      64           4 :       i2=i
      65          19 :       dm=dist
      66           1 :   100 continue
      67           1 :       IF ( i2==0 )  CALL juDFT_error(" tetcon1 ",calledby ="tetcon")
      68           1 :       dnorm=sqrt(dm)
      69           1 :       ntetra(2,1)=i2
      70             : C
      71             : C FIND POINT NEAREST TO (1+I2)/2 , NOT ON THE LINE
      72             : C CONNECTING 1 AND I2
      73             : C
      74           1 :       dm=10.0**9
      75           1 :       i3=0
      76          20 :       do 200 i=2,nkpt
      77          19 :       if ( i.eq.i2 ) goto 200
      78             :       dl=0.0
      79             :       dist=0.0
      80         126 :       do 210 icom=1,3
      81          54 :       vect=kvc(icom,i)-0.5*(kvc(icom,1)+kvc(icom,i2))
      82          54 :       dist=dist+vect*vect
      83          54 :       dl=dl+vect*(kvc(icom,1)-kvc(icom,i2))
      84          18 :   210 continue
      85          18 :       dl=dl/dnorm
      86          18 :       if ( abs(dist-dl*dl).lt.0.01*dist ) goto 200
      87          17 :       if ( dist.gt.dm ) goto 200
      88           3 :       dm=dist
      89          19 :       i3=i
      90           1 :   200 continue
      91           1 :       IF ( i3==0 )  CALL juDFT_error(" tetcon2 ",calledby ="tetcon")
      92           1 :       ntetra(3,1)=i3
      93             : C
      94             : C FIND POINT NEAREST TO (1+I2+I3)/3 , NOT IN THE
      95             : C PLANE (1,I2,I3)
      96             : C
      97             :       cn(1)=(kvc(2,1)-kvc(2,i2))*(kvc(3,i2)-kvc(3,i3))-
      98           1 :      $      (kvc(3,1)-kvc(3,i2))*(kvc(2,i2)-kvc(2,i3))
      99             :       cn(2)=(kvc(3,1)-kvc(3,i2))*(kvc(1,i2)-kvc(1,i3))-
     100           1 :      $      (kvc(1,1)-kvc(1,i2))*(kvc(3,i2)-kvc(3,i3))
     101             :       cn(3)=(kvc(1,1)-kvc(1,i2))*(kvc(2,i2)-kvc(2,i3))-
     102           1 :      $      (kvc(2,1)-kvc(2,i2))*(kvc(1,i2)-kvc(1,i3))
     103           1 :       dnorm=0.0
     104           4 :       do 300 icom=1,3
     105           3 :       dnorm=dnorm+cn(icom)**2
     106           1 :   300 continue
     107           1 :       dnorm=1.0/sqrt(dnorm)
     108           4 :       do 310 icom=1,3
     109           3 :       cn(icom)=cn(icom)*dnorm
     110           1 :   310 continue
     111             :       i4=0
     112             :       dm=10.0**9
     113          39 :       do 400 i=2,nkpt
     114          19 :       if ( (i.eq.i2).or.(i.eq.i3) ) goto 400
     115             :       dist=0.0
     116             :       dl=0.0
     117         119 :       do 410 icom=1,3
     118             :       vect=kvc(icom,i)-(kvc(icom,1)+kvc(icom,i2)+
     119          51 :      $     kvc(icom,i3))/3.0
     120          51 :       dist=dist+vect**2
     121          51 :       dl=dl+vect*cn(icom)
     122          17 :   410 continue
     123          17 :       if ( dl*dl.lt.0.01*dist ) goto 400
     124           9 :       if ( dist.gt.dm ) goto 400
     125           3 :       i4=i
     126           3 :       dm=dist
     127          19 :       vt=dl/(dnorm*6.0)
     128           1 :   400 continue
     129           1 :       IF ( i4==0 )  CALL juDFT_error(" tetcon3 ",calledby ="tetcon")
     130           1 :       ntetra(4,1)=i4
     131           1 :       voltet(1)=abs(vt)
     132             : C
     133             : C ENTER LOOP FOR CONSTRUCTION OF TETRAHEDRA:
     134             : C
     135           1 :       nt=1
     136          44 :       it=0
     137             :  1000 continue
     138          44 :       it=it+1
     139             : C
     140             : C CHECK THE SIDES OF TETRAHEDRON IT
     141             : C
     142         220 :       do 1100 j=1,4
     143             : C
     144             : C CHECK SIDE OPPOSITE TO CORNER J:
     145             : C
     146             :       ic=0
     147        1584 :       do 1200 i=1,4
     148         704 :       if ( i.eq.j ) goto 1200
     149         528 :       ic=ic+1
     150         528 :       nside(ic)=ntetra(i,it)
     151        2112 :       do 1300 icom=1,3
     152        1584 :       kcorn(icom,ic)=kvc(icom,ntetra(i,it))
     153         704 :  1300 continue
     154         176 :  1200 continue
     155         176 :       nside(4)=ntetra(j,it)
     156         176 :       is1=min(nside(1),nside(2),nside(3))
     157         176 :       is3=max(nside(1),nside(2),nside(3))
     158         176 :       is2=nside(1)+nside(2)+nside(3)-is1-is3
     159             : C
     160             : C CHECK IF THERE IS ALREADY A TETRAHEDRON CONNECTED
     161             : C TO THIS SIDE:
     162             : C
     163        4157 :       do nnt=1,nt
     164       36315 :         do 1310 nsid=1,4
     165       16178 :           if ( nnt.eq.it ) goto 1310
     166             :           ic=0
     167      141822 :           do 1320 i=1,4
     168       63032 :             if ( i.eq.nsid ) goto 1320
     169       47274 :             ic=ic+1
     170       63032 :             nside2(ic)=ntetra(i,nnt)
     171       15758 :  1320     continue
     172       15758 :           js1=min(nside2(1),nside2(2),nside2(3))
     173       15758 :           js3=max(nside2(1),nside2(2),nside2(3))
     174       15758 :           js2=nside2(1)+nside2(2)+nside2(3)-js1-js3
     175       15758 :           if ( (is1.eq.js1).and.(is2.eq.js2).and.(is3.eq.js3) )
     176             :      $      goto 1100
     177        3981 :  1310   continue
     178             :       end do
     179             : C
     180             : C CONSTRUCT THE OUTWARD NORMAL ON THIS SIDE
     181             : C
     182             :       norm(1)=(kcorn(2,2)-kcorn(2,1))*(kcorn(3,3)-kcorn(3,1))-
     183          77 :      $        (kcorn(3,2)-kcorn(3,1))*(kcorn(2,3)-kcorn(2,1))
     184             :       norm(2)=(kcorn(3,2)-kcorn(3,1))*(kcorn(1,3)-kcorn(1,1))-
     185          77 :      $        (kcorn(1,2)-kcorn(1,1))*(kcorn(3,3)-kcorn(3,1))
     186             :       norm(3)=(kcorn(1,2)-kcorn(1,1))*(kcorn(2,3)-kcorn(2,1))-
     187          77 :      $        (kcorn(2,2)-kcorn(2,1))*(kcorn(1,3)-kcorn(1,1))
     188          77 :       vol=0.0
     189         308 :       do 1400 i=1,3
     190         231 :       vol=vol+norm(i)*(kvc(i,ntetra(j,it))-kcorn(i,1))
     191          77 :  1400 continue
     192          77 :       vol=vol/6.0
     193          77 :       if ( vol.lt.0.0 ) goto 1500
     194         301 :       do 1600 i=1,3
     195         129 :       norm(i)=-norm(i)
     196          77 :  1600 continue
     197             :  1500 continue
     198          77 :       vol=abs(vol)
     199             : C
     200             : C STORE THE K-POINT ADDRESS IN NKADD ARRAY ACCORDING TO THE
     201             : C ORDER OF DISTANCE BETWEEN THE K-POINT AND THIS FACE.
     202             : C
     203        1617 :       do 1800 nkp=1,nkpt
     204             :       length=0.0
     205       10780 :       do 1850 i=1,3
     206             :       length=length+(kvc(i,nkp)
     207        4620 :      &     -(kcorn(i,1)+kcorn(i,2)+kcorn(i,3))/3.0)**2
     208        1540 :  1850 continue
     209        1540 :       ddist(nkp)=length
     210        1540 :       nkadd(nkp)=nkp
     211          77 :  1800 continue
     212        3003 :       do 1900 nkp=1,nkpt-1
     213        1463 :       minlen=ddist(nkp)
     214        1463 :       ik=nkp
     215       16093 :       do 1950 nkq=nkp+1,nkpt
     216       14630 :       if(ddist(nkq).gt.minlen-eps) goto 1950
     217        2901 :       minlen=ddist(nkq)
     218       14630 :       ik=nkq
     219        1463 :  1950 continue
     220        1463 :       ddist(ik)=ddist(nkp)
     221        1463 :       k=nkadd(nkp)
     222        1463 :       nkadd(nkp)=nkadd(ik)
     223        1463 :       nkadd(ik)=k
     224          77 :  1900 continue
     225             : C
     226             : C CONSTRUCT A TETRAHEDRON WHICH CONNECT THIS FACE TO
     227             : C A K-POINT.
     228             : C
     229             :       ik=0
     230        1805 :       do 2000 nkp=1,nkpt
     231        5455 :       do 2010 i=1,3
     232        2497 :       if ( nkadd(nkp).eq.nside(i) ) goto 2000
     233         684 :  2010 continue
     234        4788 :       do 2050 i=1,3
     235        2052 :       kcorn(i,4)=kvc(i,nkadd(nkp))
     236         684 :  2050 continue
     237         684 :       nside(4)=nkadd(nkp)
     238         684 :       vt=0.0
     239        2736 :       do 2100 i=1,3
     240        2052 :       vt=vt+norm(i)*(kcorn(i,4)-kcorn(i,1))
     241         684 :  2100 continue
     242         684 :       vt=vt/6.0
     243             : C
     244             : C REJECT POINT NKP IF IT IS ON THE WRONG SIDE .
     245             : C
     246         684 :       if ( vt.lt.eps*vol ) goto 2000
     247             : C
     248             : C CHECK IF THIS TETRAHEDRON INTERSECTS AN EXISTING ONE
     249             : C
     250        1093 :       do 3100 nk2=1,nt
     251        1050 :       if(nk2.eq.it) goto 3100
     252             : C
     253             : C FIRST CHECK IF TETRAHEDRON NK2 IS ON THE DANGEROUS SIDE
     254             : C
     255        7066 :       do 3150 i=1,4
     256             :       sum=0.0
     257       23506 :       do 3160 jj=1,3
     258       10074 :       sum=sum+norm(jj)*(kvc(jj,ntetra(i,nk2))-kcorn(jj,1))
     259        3358 :  3160 continue
     260        3358 :       if ( sum.gt.eps ) goto 3170
     261         677 :  3150 continue
     262         327 :       goto 3100
     263             :  3170 continue
     264             : C
     265             : C TETRAHEDRON NK2 IS POTENTIALLY SUSPECT
     266             : C
     267             :       l=0
     268        2943 :       do i=1,4
     269       12099 :         do 3200 jj=1,4
     270        5232 :           if(ntetra(i,nk2).eq.nside(jj)) goto 3200
     271             :           sum=0.0
     272       33859 :           do icom=1,3
     273       14511 :             cn(icom)=kvc(icom,ntetra(i,nk2))-kcorn(icom,jj)
     274       19348 :             sum=sum+cn(icom)*cn(icom)
     275             :           end do
     276       33859 :           do icom=1,3
     277       19348 :             cn(icom)=cn(icom)/sqrt(sum)
     278             :           end do
     279       69399 :           do k=1,l
     280             :             sum=0.0
     281      226772 :             do icom=1,3
     282      129584 :               sum=sum+cn(icom)*d(icom,k)
     283             :             end do
     284       37118 :             if((1.0-sum).lt.eps) goto 3200
     285             : C
     286             : C EXCLUDE THE VECTOR WHICH HAS THE SAME DIRECTION AS AN EXISTING ONE.
     287             : C
     288             :           end do
     289        4722 :           l=l+1
     290       18888 :           do icom=1,3
     291        5232 :             d(icom,l)=cn(icom)
     292             :           end do
     293        1308 :  3200   continue
     294             :       end do
     295         327 :       IF(l<4)  CALL juDFT_error(" tetcon9 ",calledby ="tetcon")
     296             : C
     297             : C HERE, WE HAVE A SET OF D-VECTORS WHICH CONNECT THE CORNER POINTS
     298             : C OF ONE TETRAHEDRON NK2 WITH THOSE OF THE CURRENT TETRAHEDRON TO
     299             : C BE CHECKED.
     300             : C FIND A PLANE OF WHICH ALL D-VECTORS EXIST IN ONE SIDE.
     301             : C
     302         793 :       do 3510 i=1,l-1
     303        6579 :         do 3500 jj=i+1,l
     304        6113 :           dnm(1)=d(2,i)*d(3,jj)-d(3,i)*d(2,jj)
     305        6113 :           dnm(2)=d(3,i)*d(1,jj)-d(1,i)*d(3,jj)
     306        6113 :           dnm(3)=d(1,i)*d(2,jj)-d(2,i)*d(1,jj)
     307             : C
     308             : C DNM IS THE NORMAL VECTOR TO THE PLANE GIVEN BY
     309             : C I-TH AND JJ-TH D-VECTORS.
     310             : C
     311             :           if((abs(dnm(1)).lt.eps).and.(abs(dnm(2)).lt.eps)
     312        6113 :      &       .and.(abs(dnm(3)).lt.eps)) goto 3500
     313       13514 :           do 3600 k=1,l
     314        9804 :             if((k.eq.i).or.(k.eq.jj)) goto 3600
     315             :             sum=0.0
     316       45682 :             do icom=1,3
     317       26104 :               sum=sum+d(icom,k)*dnm(icom)
     318             :             end do
     319        6526 :             if(abs(sum).lt.eps) goto 3600
     320       64778 :             do 3700 kk=1,l
     321       35113 :               if((kk.eq.i).or.(kk.eq.jj).or.(kk.eq.k)) goto 3700
     322             :               sss=0.0
     323      150962 :               do icom=1,3
     324       86264 :                 sss=sss+d(icom,kk)*dnm(icom)
     325             :               end do
     326       21566 :               if(abs(sss).lt.eps) goto 3700
     327       19082 :               if(sss*sum.lt.0.0) goto 3500
     328             : C
     329             : C IF K-TH AND KK-TH D-VECTORS EXIST IN OPPOSITE SIDE WITH
     330             : C RESPECT TO THE PLANE GIVEN BY I-TH AND JJ-TH D-VECTORS,
     331             : C WE WILL ATTEMPT THE NEXT PLANE. (GOTO 3500)
     332             : C
     333         323 :  3700       continue
     334        3710 :             goto 3100
     335             : C
     336             : C WE SUCCEED WHEN THE PLANE SATISFIES THE CONDITION,
     337             : C I.E., THE CURRENT TETRAHEDRON DOES NOT INTERSECTS NK2-TH ONE:
     338             : C WE TRY TO CHECK THE NEXT TETRAHEDRON. (GOTO 3100)
     339             : C
     340        5771 :  3600     continue
     341         466 :  3500   continue
     342           4 :  3510 continue
     343          46 :       goto 2000
     344             : C
     345             : C HERE, A TETRAHEDRON MADE OF J-TH FACE AND NKADD(NKP) K-POINT
     346             : C INTERSECTS AT LEAST ONE OF THE EXISTING TETRAHEDRONS:
     347             : C WE TRY THE NEXT K-POINT. (GOTO 2000)
     348             : C
     349          43 :  3100 continue
     350          43 :       ik=nkadd(nkp)
     351          43 :       volnew=vt
     352         860 :       goto 2500
     353             : C
     354             : C HERE, WE GET THE NEW TETRAHEDRON WHICH IS MADE OF J-TH FACE
     355             : C AND IK-TH K-POINT.
     356             : C
     357          34 :  2000 continue
     358             :       goto 1100
     359             : C
     360             : C HERE, WE COULD NOT FIND ANY NEW TETRAHEDRON USING J-TH FACE.
     361             : C WE TRY THE NEXT FACE. (GOTO 1100)
     362             : C
     363             :  2500 continue
     364          43 :       nt=nt+1
     365         430 :       do 2300 i=1,4
     366         172 :       ntetra(i,nt)=ntetra(i,it)
     367          43 :  2300 continue
     368          43 :       ntetra(j,nt)=ik
     369          43 :       voltet(nt)=abs(volnew)
     370          44 :  1100 continue
     371          44 :       if ( it.lt.nt ) goto 1000
     372           1 :       IF (nt>ndiv3 )  CALL juDFT_error(" nt>ndiv3",calledby ="tetcon")
     373             : C
     374             : C DETERMINE THE CHARACTERISTICS OF THIS DIVISION INTO
     375             : C TETRAHEDRA .
     376             : C
     377           1 :       vav=0.0
     378           1 :       vsq=0.0
     379           1 :       sav=0.0
     380           1 :       ssq=0.0
     381           1 :       vmin=10.0**9
     382           1 :       vmax=-vmin
     383           1 :       smin=vmin
     384           1 :       smax=vmax
     385          45 :       do it=1,nt
     386          44 :         vt=voltet(it)
     387          44 :         vav=vav+vt
     388          44 :         vsq=vsq+vt**2
     389          44 :         if ( vt.gt.vmax ) vmax=vt
     390          44 :         if ( vt.lt.vmin ) vmin=vt
     391         177 :         do i=1,3
     392         396 :           do j=i+1,4
     393             :             dist=0.0
     394        1848 :             do icom=1,3
     395             :               dist=dist+(kvc(icom,ntetra(i,it))-
     396        1056 :      $             kvc(icom,ntetra(j,it)))**2
     397             :             end do
     398         264 :             dist=sqrt(dist)
     399         264 :             sav=sav+dist
     400         264 :             ssq=ssq+dist*dist
     401         264 :             if ( dist.gt.smax ) smax=dist
     402         396 :             if ( dist.lt.smin ) smin=dist
     403             :           end do
     404             :         end do
     405             :       end do
     406           1 :       vav=vav/nt
     407           1 :       sav=sav/(6*nt)
     408           1 :       vsq=vsq/nt
     409           1 :       ssq=ssq/(6*nt)
     410           1 :       vsq=sqrt(vsq-vav*vav)
     411           1 :       ssq=sqrt(ssq-sav*sav)
     412           1 :       write(iofile,5000) nt,vav,vsq,vmin,vmax,sav,ssq,smin,smax
     413             :  5000 format(/,'   division into tetrahedra  ',/,
     414             :      $  '  there are      ',i5,'  tetrahedra ',/,
     415             :      $  '  volume         ',f15.10,'  +/-  ',f10.5,3x,2f10.5,/,
     416             :      $  '  side           ',f15.10,'  +/-  ',f10.5,3x,2f10.5,/)
     417             : c     write(ibfile,5000) nt,vav,vsq,vmin,vmax,sav,ssq,smin,smax
     418           1 :       write(iofile,5100) ((ntetra(j,i),j=1,4),i=1,nt)
     419             :  5100 format(4(4x,4i4))
     420             : c     write(ibfile,5100) ((ntetra(j,i),j=1,4),i=1,nt)
     421             : C CHECK IF WE HAVE THE CORRECT TOTAL VOLUME
     422           1 :       vt=omega*vav*nt*nsym/(2*pi)**3-1.0
     423           1 :       write(iofile,5200) vt
     424             : c     write(ibfile,5200) vt
     425           1 :       do 5300 i =1,nt
     426             : c     write(ibfile,'('' tetrahedra # '',i5,'' is '',d12.4)') i,voltet(i)
     427             :  5300 continue
     428             :  5200 format(/,'  voltetsum/volBZ - 1  ',d12.4)
     429             : C
     430             : C     The following statement used to have a stop in it.
     431             : C     If the word TETCON5 appears you have failed the < 1.0D-5 test.
     432             : C
     433           1 :       if ( abs(vt).gt.eps1 ) write(iofile,'(''  tetcon5  '')')
     434           1 :       RETURN
     435             :       END SUBROUTINE tetcon
     436             :       END MODULE m_tetcon

Generated by: LCOV version 1.13