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

          Line data    Source code
       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3             : ! This file is part of FLEUR and available as free software under the conditions
       4             : ! of the MIT license as expressed in the LICENSE file in more detail.
       5             : !--------------------------------------------------------------------------------
       6             : 
       7             :       MODULE m_kptgen_hybrid
       8             : 
       9             :       USE m_juDFT
      10             : 
      11             :       CONTAINS
      12             :       
      13             : ! this programm generates an aequdistant kpoint set including the
      14             : ! Gamma point; it is reduced to IBZ and written in kpts (M.B.)
      15             : 
      16             :       !Modified for types D.W.
      17             : 
      18             : 
      19           0 :       SUBROUTINE kptgen_hybrid(input,cell,sym,kpts,l_soc)
      20             : 
      21             :       USE m_types
      22             :       USE m_divi
      23             : 
      24             :       IMPLICIT NONE
      25             : 
      26             :       TYPE(t_input), INTENT(IN)    :: input
      27             :       TYPE(t_cell),  INTENT(IN)    :: cell
      28             :       TYPE(t_sym),   INTENT(IN)    :: sym
      29             :       TYPE(t_kpts),  INTENT(INOUT) :: kpts
      30             :       ! - scalars -
      31             :       LOGICAL, INTENT(IN)   ::  l_soc
      32             :       ! - local scalars -
      33             :       INTEGER               ::  i,j,k,nkpt
      34             :       INTEGER               ::  ikpt,ikpt0,nkpti
      35             :       INTEGER               ::  nsym
      36             :       ! - local arrays -
      37           0 :       INTEGER,ALLOCATABLE   ::  rot(:,:,:),rrot(:,:,:)
      38           0 :       INTEGER,ALLOCATABLE   ::  invtab(:)
      39           0 :       INTEGER,ALLOCATABLE   ::  neqkpt(:)
      40           0 :       INTEGER,ALLOCATABLE   ::  pkpt(:,:,:),kptp(:),symkpt(:),iarr(:),
      41           0 :      &                          iarr2(:)
      42           0 :       REAL,ALLOCATABLE      ::  rtau(:,:)
      43           0 :       REAL,ALLOCATABLE      ::  bk(:,:),bkhlp(:,:)
      44           0 :       REAL,ALLOCATABLE      ::  rarr(:)
      45             :       LOGICAL               ::  ldum
      46             : 
      47           0 :       IF (sum(kpts%nkpt3).EQ.0) THEN
      48             :          CALL divi(kpts%nkpt,cell%bmat,input%film,sym%nop,
      49           0 :      &             sym%nop2,kpts%nkpt3)
      50             :       END IF
      51             : 
      52           0 :       nkpt=kpts%nkpt3(1)*kpts%nkpt3(2)*kpts%nkpt3(3)
      53           0 :       ALLOCATE( bk(3,nkpt),bkhlp(3,nkpt) )
      54             : 
      55             :       ikpt = 0
      56           0 :       DO i=0,kpts%nkpt3(1)-1
      57           0 :         DO j=0,kpts%nkpt3(2)-1
      58           0 :           DO k=0,kpts%nkpt3(3)-1
      59           0 :             ikpt       = ikpt + 1
      60             :             bk(:,ikpt) = (/ 1.0*i/kpts%nkpt3(1),1.0*j/kpts%nkpt3(2),
      61           0 :      &                                     1.0*k/kpts%nkpt3(3) /)
      62             :           END DO
      63             :         END DO
      64             :       END DO
      65             :       
      66           0 :       IF( ikpt .ne. nkpt) STOP 'failure: number of k-points'
      67             : 
      68           0 :       IF( sym%invs .or. l_soc ) THEN
      69           0 :         nsym = sym%nop
      70             :       ELSE
      71           0 :         nsym = 2*sym%nop
      72             :       END IF
      73             :       
      74           0 :       ALLOCATE( rot(3,3,nsym),rtau(3,nsym) )
      75             : 
      76           0 :       DO i=1,sym%nop
      77           0 :          rot(:,:,i) = sym%mrot(:,:,i)
      78           0 :         rtau(  :,i) = sym%tau(:,i)
      79             :       END DO
      80             : 
      81           0 :       DO i = sym%nop+1,nsym
      82           0 :          rot(:,:,i) =  rot(:,:,i-sym%nop)
      83           0 :         rtau(  :,i)   = rtau(  :,i-sym%nop)
      84             :       END DO
      85             : 
      86           0 :       IF(any(rot(:,:,1)-reshape((/1,0,0,0,1,0,0,0,1/),(/3,3/)).ne.0))
      87           0 :      &  STOP 'kptgen: First symmetry operation is not the identity.'
      88             : 
      89           0 :       ALLOCATE( rrot(3,3,nsym),invtab(nsym) )    
      90             : 
      91           0 :       invtab = 0
      92             : 
      93           0 :       DO i = 1,sym%nop
      94           0 :         DO j = 1,sym%nop
      95             : 
      96             :           IF(    all(    matmul(rot(:,:,i),rot(:,:,j))
      97             :      &               .eq.reshape((/1,0,0,0,1,0,0,0,1/),(/3,3/)))
      98           0 :      &     .and.all(modulo(matmul(rot(:,:,i),rtau(:,j))+rtau(:,i),1.0)
      99           0 :      &               .lt.1e-10)                                  )THEN
     100           0 :             IF(invtab(i).ne.0) STOP 'kptgen: inverse operation
     101           0 :      &                              & already defined.'
     102           0 :             invtab(i)   = j
     103           0 :             rrot(:,:,i) = transpose_int ( rot(:,:,j) ) ! temporary fix for ifc
     104             :           END IF
     105             :         END DO
     106           0 :         IF(invtab(i).eq.0) STOP 'kptgen: inverse operation not found.'
     107             :         
     108             :       END DO
     109             : 
     110             : 
     111           0 :       DO i = sym%nop+1,nsym
     112           0 :          rrot(:,:,i) = - rrot(:,:,i-sym%nop)
     113             :       END DO
     114             : 
     115           0 :       ALLOCATE ( kptp(nkpt),symkpt(nkpt),rarr(3),iarr2(3),iarr(nkpt) )
     116           0 :       ALLOCATE ( pkpt(kpts%nkpt3(1)+1,kpts%nkpt3(2)+1,kpts%nkpt3(3)+1) )
     117           0 :       pkpt = 0
     118           0 :       DO ikpt = 1,nkpt
     119           0 :         iarr2 = nint ( bk(:,ikpt) * kpts%nkpt3 ) + 1
     120           0 :         pkpt(iarr2(1),iarr2(2),iarr2(3)) = ikpt
     121             :       END DO
     122             : 
     123           0 :       pkpt(kpts%nkpt3(1)+1,    :     ,    :     ) = pkpt(1,:,:)
     124           0 :       pkpt(    :     ,kpts%nkpt3(2)+1,    :     ) = pkpt(:,1,:)
     125           0 :       pkpt(    :     ,    :     ,kpts%nkpt3(3)+1) = pkpt(:,:,1)
     126             :       
     127           0 :       IF(any(pkpt.eq.0)) THEN
     128             :          CALL juDFT_error('kptgen: Definition of pkpt-pointer failed.',
     129           0 :      &                    calledby='kptgen_hybrid')
     130             :       END IF
     131           0 :       iarr = 1
     132             :       ldum = .false.
     133           0 :       DO i = 1,nkpt
     134           0 :         IF(iarr(i).eq.0) CYCLE
     135           0 :         kptp(i)   = i
     136           0 :         symkpt(i) = 1
     137           0 :         DO k = 2,nsym
     138           0 :           rarr  = matmul(rrot(:,:,k),bk(:,i)) * kpts%nkpt3
     139           0 :           iarr2 = nint(rarr)
     140           0 :           IF(any(abs(iarr2-rarr).gt.1e-10)) THEN
     141           0 :             WRITE(6,'(A,I3,A)') 'kptgen: Symmetry operation',k,
     142           0 :      &                        ' incompatible with k-point set.'
     143           0 :             ldum = .true.
     144             :           END IF
     145           0 :           iarr2 = modulo(iarr2,kpts%nkpt3) + 1
     146           0 :           IF(any(iarr2.gt.kpts%nkpt3)) 
     147           0 :      &    STOP 'kptgen: pointer indices exceed pointer dimensions.'
     148           0 :           j     = pkpt(iarr2(1),iarr2(2),iarr2(3))
     149           0 :           IF(j.eq.0) STOP 'kptgen: k-point index is zero (bug?)'
     150           0 :           IF(iarr(j).eq.0.or.j.eq.i) CYCLE
     151           0 :           iarr(j)   = 0
     152           0 :           kptp(j)   = i
     153           0 :           symkpt(j) = k
     154             :         END DO
     155             :       END DO
     156           0 :       IF(ldum) 
     157             :      &STOP 'kptgen: Some symmetry operations are incompatible
     158           0 :      & with k-point set.'
     159             :       i = 0
     160           0 :       DO ikpt = 1,nkpt
     161           0 :         IF(iarr(ikpt).eq.1) THEN
     162           0 :           i          = i + 1
     163           0 :           iarr(ikpt) = i
     164             :         END IF
     165             :       END DO
     166             :       nkpti = i
     167           0 :       DO ikpt = 1,nkpt
     168           0 :         IF(iarr(ikpt).eq.0) THEN
     169           0 :           i          = i + 1
     170           0 :           iarr(ikpt) = i
     171             :         END IF
     172             :       END DO
     173           0 :       bk(:,iarr)   = bk
     174           0 :       kptp         = iarr(kptp)
     175           0 :       kptp(iarr)   = kptp
     176           0 :       symkpt(iarr) = symkpt
     177           0 :       DO i=1,kpts%nkpt3(1)+1 
     178           0 :         DO j=1,kpts%nkpt3(2)+1
     179           0 :           DO k=1,kpts%nkpt3(3)+1 
     180           0 :             pkpt(i,j,k) = iarr(pkpt(i,j,k))
     181             :           END DO
     182             :         END DO
     183             :       END DO
     184             :    
     185           0 :       ALLOCATE( neqkpt(nkpti) )
     186           0 :       neqkpt = 0
     187           0 :       DO ikpt0 = 1,nkpti
     188           0 :         DO ikpt = 1,nkpt
     189           0 :           IF( kptp(ikpt) .eq. ikpt0 ) neqkpt(ikpt0) = neqkpt(ikpt0) + 1
     190             :         END DO
     191             :       END DO
     192             : 
     193             : !     Do not do any IO, but store in kpts
     194           0 :       kpts%nkpt=nkpti
     195           0 :       if (allocated(kpts%bk)) deallocate(kpts%bk)
     196           0 :       if (allocated(kpts%wtkpt)) deallocate(kpts%wtkpt)
     197           0 :       ALLOCATE(kpts%bk(3,kpts%nkpt),kpts%wtkpt(kpts%nkpt))
     198             : 
     199             :       
     200           0 :       DO ikpt=1,nkpti
     201           0 :          kpts%bk(:,ikpt)=bk(:,ikpt)
     202           0 :          kpts%wtkpt(ikpt)=neqkpt(ikpt)
     203             :       END DO
     204           0 :       kpts%posScale=1.0
     205             :       
     206             :       CONTAINS
     207             : 
     208             :       ! Returns least common multiple of the integers iarr(1:n).
     209             :       FUNCTION kgv(iarr,n)
     210             :       IMPLICIT NONE
     211             :       INTEGER              :: kgv
     212             :       INTEGER, INTENT(IN)  :: n,iarr(n)
     213             :       LOGICAL              :: lprim(2:maxval(iarr))
     214             :       INTEGER, ALLOCATABLE :: prim(:),expo(:)
     215             :       INTEGER              :: nprim,marr
     216             :       INTEGER              :: i,j,ia,k
     217             :       ! Determine prime numbers
     218             :       marr  = maxval(iarr)
     219             :       lprim = .true.
     220             :       DO i = 2,marr
     221             :         j = 2
     222             :         DO WHILE (i*j.le.marr)
     223             :           lprim(i*j) = .false.
     224             :           j          = j + 1
     225             :         END DO
     226             :       END DO
     227             :       nprim = count(lprim)
     228             :       ALLOCATE ( prim(nprim),expo(nprim) )
     229             :       j = 0
     230             :       DO i = 2,marr
     231             :         IF(lprim(i)) THEN
     232             :           j       = j + 1
     233             :           prim(j) = i
     234             :         END IF
     235             :       END DO
     236             :       ! Determine least common multiple
     237             :       expo = 0
     238             :       DO i = 1,n
     239             :         ia = iarr(i)
     240             :         IF(ia.eq.0) CYCLE
     241             :         DO j = 1,nprim
     242             :           k = 0
     243             :           DO WHILE(ia/prim(j)*prim(j).eq.ia)
     244             :             k  = k + 1
     245             :             ia = ia / prim(j)
     246             :           END DO
     247             :           expo(j) = max(expo(j),k)
     248             :         END DO
     249             :       END DO
     250             :       kgv = 1
     251             :       DO j = 1,nprim
     252             :         kgv = kgv * prim(j)**expo(j)
     253             :       END DO
     254             :       DEALLOCATE ( prim,expo )
     255             :       END FUNCTION kgv
     256             :       
     257             : 
     258             : c     ifc seems to have problems transposing integer arrays. this is a fix.
     259           0 :       FUNCTION transpose_int ( a )
     260             :       IMPLICIT NONE
     261             :       integer transpose_int(3,3),a(3,3)
     262             :       integer i,j
     263           0 :       DO i = 1,3
     264           0 :         DO j = 1,3
     265           0 :           transpose_int(i,j) = a(j,i)
     266             :         END DO
     267             :       END DO
     268             :       END FUNCTION transpose_int
     269             : 
     270             : c     function modulo1 maps kpoint into first BZ
     271             :       FUNCTION modulo1(kpoint,nkpt,a,b,c)
     272             : 
     273             :       IMPLICIT NONE
     274             : 
     275             :       INTEGER,INTENT(IN)  :: nkpt,a,b,c
     276             :       REAL, INTENT(IN)    :: kpoint(3)
     277             :       REAL                :: modulo1(3)
     278             :       INTEGER             :: help(3),nkpt3(3)
     279             : 
     280             :       nkpt3 = (/a,b,c/)
     281             :       modulo1 = kpoint*nkpt3
     282             :       help    = nint(modulo1)
     283             :       IF(any(abs(help-modulo1).gt.1e-8)) THEN
     284             :         modulo1 = kpoint*nkpt3
     285             :         WRITE(*,*) modulo1
     286             :         help    = nint(modulo1)
     287             :         WRITE(*,*) help
     288             :         WRITE(6,'(A,F5.3,2('','',F5.3),A)') 'modulo1: argument (',
     289             :      &           kpoint,') is not an element of the k-point set.'
     290             :         CALL juDFT_error(
     291             :      +     'modulo1: argument not an element of k-point set.', 
     292             :      +     calledby = 'kptgen_hybrid:modulo1')
     293             :       END IF
     294             :       modulo1 = modulo(help,nkpt3)*1.0/nkpt3
     295             : 
     296             :       END FUNCTION modulo1
     297             : 
     298             :       END SUBROUTINE kptgen_hybrid
     299             :       
     300             :       END MODULE m_kptgen_hybrid

Generated by: LCOV version 1.13