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

          Line data    Source code
       1             :       MODULE m_kpttet
       2             :       use m_juDFT
       3             :       CONTAINS
       4           1 :       SUBROUTINE kpttet(
       5             :      >                  iofile,ibfile,iokpt,
       6             :      >                  kpri,ktest,kmidtet,mkpt,ndiv3,
       7             :      >                  nreg,nfulst,rltv,voluni,
       8             :      >                  nsym,ccr,mdir,mface,
       9           1 :      >                  ncorn,nface,fdist,fnorm,cpoint,
      10           1 :      <                  voltet,ntetra,ntet,vktet,
      11             :      =                  nkpt,
      12           1 :      <                  divis,vkxyz,wghtkp)
      13             : c
      14             : c
      15             : c ---> This program generates k-points
      16             : c           in irreducible wedge of BZ  (for nreg=0)
      17             : c           in total BZ                 (for nreg=1)
      18             : c      (BZ = 1. Brillouin-zone) for all canonical Bravais lattices
      19             : c      in 3 dimensions,
      20             : c      using the basis vectors of the reciprocal lattice,
      21             : c      the corner points of the irreducible wedge of the BZ
      22             : c      and the bordering planes of the irreducible wedge.
      23             : c
      24             : c      The k-points are generated by the tetrahedron method.
      25             : c      by generating a set of k-points which are maximally far apart
      26             : c      for the rquired number of points.
      27             : c      The method and the subroutines were obtained via St.Bluegel.
      28             : c      The information about the irr wedge of the BZ
      29             : c      is taken from BRZONE.
      30             : c
      31             : c-----------------------------------------------------------------------
      32             : c    Meaning of variables:
      33             : c    INPUT:
      34             : c
      35             : c    Symmetry of lattice:
      36             : c    rltv     : cartesian coordinates of basis vectors for
      37             : c               reciprocal lattice: rltv(ix,jn), ix=1,3; jn=1,3
      38             : c    voluni   : volume of the Bravais lattice unit cell
      39             : c    nsym     : number of symmetry elements of points group
      40             : c    ccr     : rotation matrix for symmetry element
      41             : c                   in cartesian representation
      42             : c
      43             : c    representation of the irreducible part of the BZ:
      44             : c    fnorm    : normal vector of the planes bordering the irrBZ
      45             : c    fdist    : distance vector of the planes bordering the irrBZ
      46             : c    ncorn    : number of corners of the irrBZ
      47             : c    nface    : number of faces of the irrBZ
      48             : c    cpoint   : cartesian coordinates of corner points of irrBZ
      49             : c
      50             : c    characterization of the tetrahedron-method k-point set:
      51             : c    nreg     : 1 kpoints in full BZ; 0 kpoints in irrBZ
      52             : c    nfulst   : 1 kpoints ordered in full stars 
      53             : c                  (meaningful only for nreg =1; full BZ)
      54             : c    nkpt     : on input: required number of k-points inside irrBZ
      55             : c               to build the tetrahedrons
      56             : c    ntet     : number of tetrahedra generated
      57             : c    ntetra   : list of four points for each tetrahedron
      58             : c               containing the indices of the respective corner points
      59             : c    vktet    : corner points of tetrahedra
      60             : c
      61             : c    kmidtet  : key to generate mid-tetrahedron k-points
      62             : c               1 mid-points are generated; 0 not generated
      63             : c
      64             : c    OUTPUT: k-point set
      65             : c    nkpt     : number of k-points generated in set
      66             : c    vkxyz    : vector of kpoint generated; in cartesian representation
      67             : c    wghtkp   : weight associated with k-points for BZ integration
      68             : c    divis    : integer triple divis(i); i=1,4.
      69             : c               Used to find more accurate representation of k-points
      70             : c               vklmn(i,kpt)/divis(i) and weights as wght(kpt)/divis(4)
      71             : c
      72             : c-----------------------------------------------------------------------
      73             :       USE m_constants, ONLY : pimach
      74             :       USE m_tetcon
      75             :       USE m_kvecon
      76             :       USE m_fulstar
      77             :       IMPLICIT NONE
      78             : C
      79             : C-----> PARAMETER STATEMENTS
      80             : C
      81             :       INTEGER, INTENT (IN) :: mkpt,ndiv3,mface,mdir
      82             : c
      83             : c ---> file number for read and write
      84             : c
      85             :       INTEGER, INTENT (IN) :: iofile,iokpt,ibfile
      86             : c
      87             : c ---> running mode parameter
      88             : c
      89             :       INTEGER, INTENT (IN) :: kpri,ktest
      90             : C
      91             : C----->  Symmetry information
      92             : C
      93             :       INTEGER, INTENT (IN) :: nsym
      94             :       REAL,    INTENT (IN) :: ccr(3,3,48)
      95             : C
      96             : C----->  BRAVAIS LATTICE INFORMATION
      97             : C
      98             :       REAL,    INTENT (IN) ::  voluni
      99             : C
     100             : C----->  RECIPROCAL LATTICE INFORMATION
     101             : C
     102             :       INTEGER, INTENT (IN) :: ncorn,nface
     103             :       REAL,    INTENT (IN) :: rltv(3,3),fnorm(3,mface),fdist(mface)
     104             :       REAL,    INTENT (IN) :: cpoint(3,mface)
     105             : C
     106             : C----->  BRILLOUINE ZONE INTEGRATION
     107             : C
     108             :       INTEGER, INTENT (IN) :: nreg,nfulst,kmidtet
     109             :       INTEGER, INTENT (INOUT) :: nkpt
     110             :       INTEGER, INTENT (OUT) :: ntetra(4,ndiv3),ntet
     111             :       REAL,    INTENT (OUT) :: voltet(ndiv3),vktet(3,mkpt)
     112             :       REAL,    INTENT (OUT) :: vkxyz(3,mkpt),wghtkp(mkpt),divis(4)
     113             : 
     114             : C
     115             : C --->  local variables
     116             : c
     117             :       INTEGER   i,j,ii, nkstar
     118             :       REAL      sumwght,eps,one,tpi,sumvol,volirbz
     119           2 :       REAL      vkmid(3,mkpt)
     120             : C
     121             : C --->  set local constants
     122             : c
     123             :       SAVE      eps,one
     124             :       DATA      eps/1.0e-9/,one/1.0/
     125             : c
     126             : c======================================================================
     127             : c
     128           1 :       tpi = 2.0 * pimach()
     129             : c
     130           1 :       IF (kpri.GE.1) THEN
     131           0 :         WRITE (iofile,'(3x,'' *<* kpttet *>* '')')
     132           0 :         WRITE (iofile,'(3x,'' generate k-vectors'')')
     133           0 :         WRITE (iofile,'(3x,'' ~~~~~~~~~~~~~~~~~~'')')
     134           0 :         WRITE (iofile,'(3x,'' by tetrahedron-method'')')
     135           0 :         WRITE (iofile,'(3x,'' ~~~~~~~~~~~~~~~~~~~~~'')')
     136             :       ENDIF
     137             :  
     138             :       WRITE (iofile,'('' k-points generated with tetrahedron '',
     139           1 :      >                                              ''method'')')
     140             :       WRITE (iokpt,'(''# k-points generated with tetrahedron '',
     141           1 :      >                                              ''method'')')
     142           1 :       IF (nreg .EQ. 0) THEN
     143           1 :         WRITE (iofile,'(3x,'' in irred wedge of 1. Brillouin zone'')')
     144           1 :         WRITE (iofile,'(3x,'' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'')')
     145           1 :         WRITE (iokpt,'(''#'',i4,21x,''nreg: k-points in irrBZ'')') nreg
     146           0 :       ELSEIF (nreg .eq. 1 .and. nfulst .eq. 1) then
     147           0 :         WRITE (iofile,'(3x,'' in 1. Brillouin zone'')')
     148           0 :         WRITE (iofile,'(3x,'' ~~~~~~~~~~~~~~~~~~~~'')')
     149           0 :         WRITE (iofile,'(3x,'' full stars generated'')')
     150           0 :         WRITE (iofile,'(3x,'' ~~~~~~~~~~~~~~~~~~~~'')')
     151             :         WRITE (iokpt,'(''#'',2(i4,1x),14x,'' nreg,nfulst: '',
     152           0 :      >      ''k-points in totBZ, ordered in full stars'')') nreg,nfulst
     153             :       ELSE
     154             :         WRITE (iofile,'(2(1x,i4),4x,'' nreg,nfulst: wrong choice;'',
     155             :      >       /,27x,'' allowed combinations: (1,1); (0,0),(0,1)'')' )
     156           0 :      >                                                   nreg,nfulst
     157           0 :          CALL juDFT_error("nreg,nfulst: wrong choice",calledby="kpttet")
     158             :       ENDIF
     159             : 
     160             :       CALL kvecon(
     161             :      >            iofile,ibfile,mkpt,mface,
     162             :      >            nkpt,ncorn,nsym,nface,rltv,fdist,fnorm,cpoint,
     163           1 :      <            vktet )
     164             : !
     165             : ! --->  generate tetrahedra and mid-tetrahedron k-points
     166             : !
     167             : ! --->  (a) Determine the corner K-POINTs for X number of Tetrahedra for
     168             : !           doing a very pretty Brillouine zone Integration;
     169             : ! --->      determine the volume of each tetrahedron
     170             : !
     171             :       CALL tetcon(
     172             :      >            iofile,ibfile,mkpt,ndiv3,
     173             :      >            nkpt,voluni,vktet,
     174             :      =            nsym,
     175           1 :      <            ntet,voltet,ntetra)
     176             : c
     177           1 :       WRITE (iofile,'('' the number of tetrahedra '')')
     178           1 :       WRITE (iofile,*) ntet
     179           1 :       WRITE (iofile,'('' volumes of the tetrahedra '')')
     180             :       WRITE (iofile,'(e19.12,1x,i5,5x,''voltet(i),i'')')
     181           1 :      >                               (voltet(i),i,i=1,ntet)
     182           1 :       WRITE (iofile,'('' corners of the tetrahedra '')')
     183           1 :       WRITE (iofile, 999) ((ntetra(j,i),j=1,4),i=1,ntet)
     184           1 :       WRITE (iofile,'('' the # of different k-points '')')
     185           1 :       WRITE (iofile,*) nkpt
     186           1 :       WRITE (iofile,'('' k-points used to construct tetrahedra'')')
     187           1 :       WRITE (iofile,'(3(4x,f10.6))') ((vktet(i,j),i=1,3),j=1,nkpt)
     188             :   999 FORMAT (4(3x,4i4))
     189             : c
     190             : c --->   calculate weights from volume of tetrahedra
     191             : c
     192           1 :       volirbz =  tpi**3 /(real(nsym)*voluni)
     193           1 :       sumvol = 0.0
     194          45 :       DO i = 1, ntet
     195          44 :          sumvol = sumvol + voltet(i)
     196          45 :          voltet(i) = ntet * voltet(i) / volirbz 
     197             :       ENDDO
     198             : c
     199           1 :       IF ((sumvol-volirbz)/volirbz .LE. eps) THEN
     200           1 :         IF (kmidtet.EQ.1) THEN
     201           0 :           DO i = 1, ntet
     202           0 :             wghtkp(i) = voltet(i)/sumvol
     203             :           ENDDO
     204             :         ELSE
     205          21 :           DO i = 1, nkpt
     206           1 :             wghtkp(i) = 1./nkpt
     207             :           ENDDO
     208             :         ENDIF
     209             :       ELSE
     210             :         WRITE (iofile, '(2(e19.12,1x),5x,''summvol.ne.volirbz'')')
     211           0 :      >                                     sumvol,volirbz
     212           0 :          CALL juDFT_error("sumvol =/= volirbz",calledby="kpttet")
     213             :       ENDIF
     214             : c
     215             : c --->  prepare the final set of kpoints in irrBZ (depending on kmidtet)
     216             : c
     217           1 :       IF ( kmidtet.EQ.0) THEN
     218             : c
     219          21 :         DO i = 1, nkpt
     220          20 :            vkxyz(:,i) = vktet(:,i)
     221             :            WRITE (iofile,'(3(f10.7,1x),f12.10,1x,i4,3x,
     222          21 :      +           ''vkxyz, wghtkp'')') (vkxyz(ii,i),ii=1,3),wghtkp(i),i
     223             :         ENDDO  
     224           1 :         nkstar = nkpt
     225             : 
     226           0 :       ELSEIF ( kmidtet.EQ.1) THEN
     227             : c
     228             : c --->  (b) calculate mid-tetrahedron k-points
     229             : c
     230           0 :          DO i=1,ntet
     231           0 :             vkmid(:,i) = 0.0
     232           0 :             DO j=1,4
     233           0 :                vkmid(:,i) = vkmid(:,i) + vktet(:,ntetra(j,i))
     234             :             ENDDO
     235           0 :             vkmid(:,i) = vkmid(:,i) * 0.25
     236             :          ENDDO
     237             : 
     238           0 :          nkpt = ntet
     239           0 :          WRITE (iofile,'('' the new number of k-points is '',i4)') nkpt
     240             :          WRITE (iofile,'('' the new k-points are the '',
     241           0 :      +                        ''mid-tetrahedron-points '')')
     242             :          WRITE (iokpt,'(''# the new k-points are the '',
     243           0 :      +                        ''mid-tetrahedron-points '')')
     244           0 :          sumwght = 0.00
     245           0 :          DO i=1,ntet
     246           0 :            vkxyz(:,i) = vkmid(:,i)
     247           0 :            sumwght = sumwght + wghtkp(i)
     248             :          ENDDO
     249             : !
     250             : ! ---> check sumwght; if abs(sumwght-1).lt.eps print kpoints and weights
     251             : !
     252           0 :          IF ( abs(sumwght - one).LT.eps) THEN
     253             :             WRITE (iofile,'(1x,f12.10,1x,'' sumwght .eq. one'')')
     254           0 :      +                                                   sumwght
     255           0 :             DO i=1,nkpt
     256             :                WRITE (iofile,'(3(f10.7,1x),f12.10,1x,i4,3x,
     257           0 :      +            ''vkxyz, wghtkp'')') (vkxyz(ii,i),ii=1,3),wghtkp(i), i
     258             :             ENDDO
     259           0 :             nkstar = ntet
     260             : 
     261             :          ELSE
     262             :             WRITE (iofile,'(1x,f12.10,1x,'' sumwght .ne. one'')')
     263           0 :      +                                                   sumwght
     264           0 :              CALL juDFT_error("sumwght",calledby="kpttet")
     265             :          ENDIF
     266             : 
     267             :       END IF ! end of generation of mid-tetrahedron k-points
     268             : 
     269             : !
     270             : ! --->   set denominators for more accurate k-point representation
     271             : !
     272           5 :       DO i=1,4
     273           5 :          divis(i) = real(nkpt)
     274             :       ENDDO
     275             : 
     276           1 :       IF ( nreg.EQ.1 .AND. nfulst.EQ.1 ) THEN
     277             : 
     278             : ! --->   generate full stars for all representative k-points
     279             : !        - for nreg=1 and nfulst=1:
     280             : !              - determine order of full star ifstar(kpn).le.nsym
     281             : !              - assign nkpt= sum {ifstar(ik)} (ik=1,ntet)
     282             : !              - assign vkxyz(ix,kpn) = vkstar(ix,ikpn(is,ik));
     283             : !                        ix=1,3; kpn=1,nkpt; ik=1,ntet; is=1,ifstar(ik)
     284             : !              - calculate wghtkp(kpn)=wghtkp_old(ik)/ifstar(ik)
     285             : !                                kpn=1,nkpt; ik=1,ntet
     286             : 
     287             :          CALL fulstar(
     288             :      >                iofile,iokpt,kpri,ktest,
     289             :      >                ccr,nsym,
     290             :      >                vkxyz,nkstar,mkpt,mface,mdir,
     291           0 :      =                nkpt,vkxyz,wghtkp)
     292             : 
     293           0 :          divis(4) = divis(4) * nsym
     294             :       ENDIF
     295             : 
     296           1 :       RETURN
     297             :       END SUBROUTINE kpttet
     298             :       END MODULE m_kpttet

Generated by: LCOV version 1.13