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

          Line data    Source code
       1             :       MODULE m_kvecon
       2             :       use m_juDFT
       3             : !
       4             : ! This subroutine determines the k-points with which we
       5             : ! will calculate the band-structure. The first ncorn
       6             : ! points are the corners of the irreducible wedge, as
       7             : ! determined in subroutine bzone. We then generate
       8             : ! nkpt-ncorn points inside the wedge, distributed in
       9             : ! such a way that the minimal distance between the points
      10             : ! is maximal. This subroutine does not find the optimal
      11             : ! distribution of k-points, it only finds a reasonable one.
      12             : !
      13             :       CONTAINS
      14           1 :       SUBROUTINE kvecon(
      15             :      >                  iofile,ibfile,mkpt,mface,
      16           1 :      >                  nkpt,ncorn,nsym,nface,rcmt,fdist,fnorm,cpoint,
      17           1 :      <                  kvc )
      18             : 
      19             :       IMPLICIT NONE
      20             : 
      21             : ! ... Arguments ...
      22             :       INTEGER, INTENT (IN) :: mkpt,mface
      23             :       INTEGER, INTENT (IN) :: iofile,ibfile
      24             :       INTEGER, INTENT (IN) :: nkpt,ncorn,nsym,nface
      25             :       REAL,    INTENT (IN) :: fdist(mface),fnorm(3,mface)
      26             :       REAL,    INTENT (IN) :: rcmt(3,3),cpoint(3,mface)
      27             :       REAL,    INTENT (OUT) :: kvc(3,mkpt)
      28             : 
      29             : ! ... Locals ...
      30             :       INTEGER nk,i,lmin,l,j,n1,n2,n3,i1,i2,i3,ncd
      31             :       REAL d,dm,dist,dmax,xmin,cn,alpha,thrd
      32           2 :       REAL knew(3),kc(3),cnorm(3),xnorm(3),cand(3,48*mkpt)
      33             : !
      34             : ! ... Intrinsic Functions ...
      35             :       INTRINSIC abs,sqrt
      36             : 
      37           1 :       IF ( nkpt.LT.ncorn ) THEN
      38           0 :         WRITE (iofile,'(1x,''nkpt='',i4,'' ,ncorn='',i4)') nkpt,ncorn
      39           0 :          CALL juDFT_error("nkpt<ncorn ",calledby="kvecon")
      40             :       ENDIF
      41           1 :       IF ( nkpt.GT.mkpt ) THEN
      42           0 :         WRITE (iofile,'(1x,''nkpt='',i4,'' , mkpt='',i4)') nkpt, mkpt
      43           0 :          CALL juDFT_error("nkpt>mkpt ",calledby="kvecon")
      44             :       ENDIF
      45             :       thrd=1.00/3.00
      46             : !
      47             : ! Construct list of candidate points
      48             : !
      49             :       d = 1.00
      50           7 :       DO  i = 1, 3
      51             :         dist = 0.00
      52          21 :         DO  j=1,3
      53          12 :           dist = dist + rcmt(j,i)**2
      54             :         ENDDO
      55           3 :         dist = sqrt(dist)
      56           3 :         xnorm(i) = dist
      57           4 :         d = d * dist
      58             :       ENDDO
      59             : 
      60           1 :       d = thrd*(d/(nkpt*nsym))**thrd
      61           1 :       n1 = int(xnorm(1)/d)+1
      62           1 :       n2 = int(xnorm(2)/d)+1
      63           1 :       n3 = int(xnorm(3)/d)+1
      64             :       WRITE( iofile,'('' n1 = '',i5,'' n2 = '',i5,'' n3 = '',i5)') 
      65           1 :      +                                                    n1,n2,n3
      66           1 :       ncd=0
      67             : 
      68           1 :       WRITE (iofile,'('' $$$$$ check $$$$$ '')')
      69           1 :       WRITE (iofile,'(''   ncorn = '',i4)') ncorn
      70           1 :       WRITE (ibfile,'(''    cpoints '')')
      71           7 :       DO i = 1, ncorn
      72           7 :         WRITE (ibfile,8901) cpoint(1,i),cpoint(2,i),cpoint(3,i)
      73             :       ENDDO
      74             :  8901 FORMAT ('    ( ',2(f13.6,','),f13.6,' )',/)
      75           1 :       WRITE (ibfile,'(''   nface = '',i4)') nface
      76           6 :       DO i = 1, nface
      77           5 :         WRITE (ibfile,8902) fdist(i)
      78           6 :         WRITE (ibfile,8901) fnorm(1,i),fnorm(2,i),fnorm(3,i)
      79             :       ENDDO
      80             :  8902 FORMAT ('  fdist = ',f13.6)
      81             : C
      82          50 :       DO i1=-n1,n1
      83        2451 :         DO i2=-n2,n2
      84       86485 :           DO i3=-n3,n3
      85      588245 :             DO i = 1, 3
      86             :               knew(i) = i1*rcmt(i,1)/n1 +
      87             :      +                  i2*rcmt(i,2)/n2 + 
      88      336140 :      +                  i3*rcmt(i,3)/n3
      89             :             ENDDO
      90      357633 :             DO l = 1, nface
      91             :               alpha=0.0e0
      92     1540105 :               DO i = 1, 3
      93      880060 :                 alpha = alpha + knew(i)*fnorm(i,l)
      94             :               ENDDO
      95      220834 :               IF ( alpha.GT.( fdist(l)+1.0e-6*d ) ) GOTO 200
      96             :             ENDDO
      97         819 :             ncd = ncd + 1
      98         819 :             IF (ncd>48*mkpt)  CALL juDFT_error("ncd>ncmax",calledby
      99           0 :      +           ="kvecon")
     100        5733 :             DO i = 1, 3
     101       84035 :               cand(i,ncd)=knew(i)
     102             :             ENDDO 
     103        2401 :   200       CONTINUE
     104             :           ENDDO
     105             :         ENDDO
     106             :       ENDDO
     107           1 :       IF (ncd<nkpt)  CALL juDFT_error("ncd<nkpt",calledby ="kvecon")
     108             : 
     109             : ! Initialize the kpoints
     110          13 :       DO nk=1,ncorn
     111          43 :         DO i=1,3
     112          24 :           kvc(i,nk) = cpoint(i,nk)
     113             :         ENDDO
     114             :       ENDDO
     115          15 :       DO nk=ncorn+1,nkpt
     116          99 :         DO i=1,3
     117          56 :           kvc(i,nk) = cpoint(i,1)
     118             :         ENDDO
     119             :       ENDDO
     120           1 :       IF ( nkpt.eq.ncorn ) GOTO 2
     121             : !
     122             : ! Enter dynamical cycle. We find a k-point with the shortest
     123             : ! distance to the other k-points, take it out, and put it in
     124             : ! a position where it has a maximal minimal distance to all
     125             : ! the other kpoints
     126             : !
     127             :     1 continue
     128             : C
     129             : C FIND MINIMAL DISTANCE K-POINT , EXCLUDE CORNER POINTS
     130             : C
     131          15 :       lmin=0
     132          15 :       d=1.0e6
     133         225 :       do 2000 l=ncorn+1,nkpt
     134             :       dm=1.0e6
     135        8610 :       do 2100 j=1,nkpt
     136        4200 :       if ( j.eq.l ) goto 2100
     137             :       dist=0.0e0
     138       27930 :       do 2200 i=1,3
     139       11970 :       dist=dist+(kvc(i,j)-kvc(i,l))**2
     140        3990 :  2200 continue
     141        3990 :       if ( dist.lt.(dm*(1.0-1e-6))  ) dm=dist
     142         210 :  2100 continue
     143         210 :       if ( d.lt.(dm*(1.0+1e-6))  ) goto 2000
     144         106 :       d=dm
     145         210 :       lmin=l
     146          15 :  2000 continue
     147          15 :       if ( lmin.eq.0 )  CALL juDFT_error(" lmin=0 ",calledby="kvecon")
     148             : C
     149             : C Find the point which is furthest from all k-points. We
     150             : C restrict our points to be on a finite grid determined by
     151             : C the corner points.
     152             : C
     153          15 :       dmax=0.0e0
     154             :       write(iofile,'('' kpoint ('',i3,'')  ncd = '',i8,
     155          15 :      +               '' max. allowed is '',i8)') lmin,ncd,48*mkpt
     156       12300 :       do 3000 n1=1,ncd
     157       85995 :       do 3100 i=1,3
     158       36855 :       knew(i)=cand(i,n1)
     159       12285 :  3100 continue
     160             :       dm=1.0e6
     161      503685 :       do 3200 j=1,nkpt
     162      245700 :       if ( j.eq.lmin ) goto 3200
     163             :       dist=0.0e0
     164     1633905 :       do 3300 i=1,3
     165      700245 :       dist=dist+(kvc(i,j)-knew(i))**2
     166      233415 :  3300 continue
     167      233415 :       if ( dist.lt.(dm*(1.0-1e-6))  ) dm=dist
     168       12285 :  3200 continue
     169       12285 :       if ( dm.lt.(dmax*(1.0+1e-6)) ) goto 3000
     170             :       dmax=dm
     171        1512 :       do 3400 i=1,3
     172         648 :       kc(i)=knew(i)
     173       12285 :  3400 continue
     174          15 :  3000 continue
     175             : C
     176             : C WE HAVE FOUND A NEW POINT KC . IF THE MINIMAL DISTANCE OF
     177             : C THIS POINT IS LESS THAN THE DISTANCE D WE ALREADY HAVE ,
     178             : C WE STOP THE DYNAMIC CYCLE .
     179             : C
     180          15 :       if ( dmax.lt.1.0001*d ) goto 4000
     181          98 :       do 3500 i=1,3
     182          42 :       kvc(i,lmin)=kc(i)
     183          14 :  3500 continue
     184           1 :       goto 1
     185             :  4000 continue
     186             : C
     187             : C WE HAVE FOUND A REASONABLE DISTRIBUTION . TO AVOID TETRAHEDRA
     188             : C WITH A FUNNY SHAPE , WE PUT ALL THE POINTS WHICH ARE NEAR A
     189             : C SIDE ON THIS SIDE .
     190             : C
     191           1 :       d=0.49*sqrt(d)
     192          15 :       do 4100 l=ncorn+1,nkpt
     193          14 :       xmin=1.0e6
     194          70 :       do n1=1,ncorn-2
     195         210 :        do n2=n1+1,ncorn-1
     196         476 :         do 4200 n3=n2+1,ncorn
     197             :         cnorm(1)=(cpoint(2,n1)-cpoint(2,n2))*(cpoint(3,n2)-cpoint(3,n3))
     198         280 :      &          -(cpoint(3,n1)-cpoint(3,n2))*(cpoint(2,n2)-cpoint(2,n3))
     199             :         cnorm(2)=(cpoint(3,n1)-cpoint(3,n2))*(cpoint(1,n2)-cpoint(1,n3))
     200         280 :      &          -(cpoint(1,n1)-cpoint(1,n2))*(cpoint(3,n2)-cpoint(3,n3))
     201             :         cnorm(3)=(cpoint(1,n1)-cpoint(1,n2))*(cpoint(2,n2)-cpoint(2,n3))
     202         280 :      &          -(cpoint(2,n1)-cpoint(2,n2))*(cpoint(1,n2)-cpoint(1,n3))
     203         280 :         cn=0.0e0
     204        1120 :         do i=1,3
     205        1120 :           cn=cn+cnorm(i)**2
     206             :         end do
     207         280 :         cn=1.0e0/sqrt(cn)
     208        1120 :         do i=1,3
     209        1120 :           cnorm(i)=cnorm(i)*cn
     210             :         end do
     211             :         alpha=0.0e0
     212        1960 :         do i=1,3
     213        1120 :           alpha=alpha+cnorm(i)*(kvc(i,l)-cpoint(i,n1))
     214             :         end do
     215         280 :         if ( abs(alpha).gt.abs(xmin) ) goto 4200
     216             :         xmin=alpha
     217         518 :         do i=1,3
     218         280 :           xnorm(i)=cnorm(i)
     219             :         end do
     220         140 :  4200   continue
     221             :        end do
     222             :       end do
     223          14 :       if ( abs(xmin).gt.d ) goto 4100
     224          98 :       do 4700 i=1,3
     225          42 :       kvc(i,l)=kvc(i,l)-xmin*xnorm(i)
     226          14 :  4700 continue
     227           1 :  4100 continue
     228             : C
     229             : C IF A POINT IS TOO CLOSE TO AN EDGE , WE PUT IT ON THAT EDGE .
     230             : C
     231          29 :       do 5000 l=ncorn+1,nkpt
     232          14 :       xmin=1.0e6
     233          84 :       do n1=1,ncorn-1
     234         294 :         do 5100 n2=n1+1,ncorn
     235             :           cnorm(1)=(cpoint(2,n1)-cpoint(2,n2))*(cpoint(3,n2)-kvc(3,l))
     236         210 :      &            -(cpoint(3,n1)-cpoint(3,n2))*(cpoint(2,n2)-kvc(2,l))
     237             :           cnorm(2)=(cpoint(3,n1)-cpoint(3,n2))*(cpoint(1,n2)-kvc(1,l))
     238         210 :      &            -(cpoint(1,n1)-cpoint(1,n2))*(cpoint(3,n2)-kvc(3,l))
     239             :           cnorm(3)=(cpoint(1,n1)-cpoint(1,n2))*(cpoint(2,n2)-kvc(2,l))
     240         210 :      &            -(cpoint(2,n1)-cpoint(2,n2))*(cpoint(1,n2)-kvc(1,l))
     241         210 :           cn=0.0e0
     242         840 :           do i=1,3
     243         210 :             cn=cn+cnorm(i)**2
     244             :           end do
     245             :           dist=0.0e0
     246        1470 :           do i=1,3
     247         840 :             dist=dist+(cpoint(i,n1)-cpoint(i,n2))**2
     248             :           end do
     249         210 :           cn=sqrt(cn/dist)
     250         210 :           if ( cn.gt.xmin ) goto 5100
     251             :           xmin=cn
     252             :           cn=0.0e0
     253         350 :           do i=1,3
     254         200 :             cn=cn+(cpoint(i,n2)-cpoint(i,n1))*(kvc(i,l)-cpoint(i,n1))
     255             :           end do
     256          50 :           cn=cn/dist
     257         200 :           do i=1,3
     258         210 :             knew(i)=cpoint(i,n1)+cn*(cpoint(i,n2)-cpoint(i,n1))
     259             :           end do
     260          70 :  5100   continue
     261             :       end do
     262          14 :       if ( xmin.gt.d ) goto 5000
     263          49 :       do 5600 i=1,3
     264          21 :       kvc(i,l)=knew(i)
     265          14 :  5600 continue
     266           1 :  5000 continue
     267             :     2 CONTINUE
     268             : C
     269             : C OUTPUT OF KVECTORS
     270             : C
     271           1 :       write(iofile,6000) 2.0e0*d,nkpt
     272             :  6000 format(//,'   minimal distance  ',f10.5,
     273             :      & //,'    coordinates of',i5,'  kpoints in cart. units ',//)
     274           1 :       write(iofile,6100) ((kvc(i,l),i=1,3),l=1,nkpt)
     275             :  6100 format(2(3f13.6,3x),3f13.6)
     276           1 :       RETURN
     277             :       END SUBROUTINE kvecon
     278             :       END MODULE m_kvecon

Generated by: LCOV version 1.13