LCOV - code coverage report
Current view: top level - kpoints - julia.f90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 152 178 85.4 %
Date: 2019-09-08 04:53:50 Functions: 2 2 100.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             : MODULE m_julia
       7             : 
       8             : USE m_juDFT
       9             : 
      10             : CONTAINS
      11             : 
      12          27 : SUBROUTINE julia(sym,cell,input,noco,banddos,kpts,l_q,l_fillArrays)
      13             : !----------------------------------------------------------------------+
      14             : ! Generate a k-point file with approx. nkpt k-pts or a Monkhorst-Pack  |
      15             : ! set with nmod(i) divisions in i=x,y,z direction. Interface to kptmop |
      16             : ! and kpttet routines of the MD-programm.                              |
      17             : !                                                          G.B. 07/01  |
      18             : !----------------------------------------------------------------------+
      19             : 
      20             :    USE m_constants
      21             :    USE m_bravais
      22             :    USE m_divi
      23             :    USE m_brzone
      24             :    USE m_brzone2
      25             :    USE m_kptmop
      26             :    USE m_kpttet
      27             :    USE m_bandstr1
      28             :    USE m_types
      29             : 
      30             :    IMPLICIT NONE
      31             : 
      32             :    TYPE(t_sym),     INTENT(IN)    :: sym
      33             :    TYPE(t_cell),    INTENT(IN)    :: cell
      34             :    TYPE(t_input),   INTENT(IN)    :: input
      35             :    TYPE(t_noco),    INTENT(IN)    :: noco
      36             :    TYPE(t_banddos), INTENT(IN)    :: banddos
      37             :    TYPE(t_kpts),    INTENT(INOUT) :: kpts
      38             : 
      39             :    LOGICAL, INTENT (IN)           :: l_q, l_fillArrays
      40             :    
      41             :    INTEGER, PARAMETER :: nop48  = 48
      42             :    INTEGER, PARAMETER :: mface  = 51
      43             :    INTEGER, PARAMETER :: mdir   = 10
      44             :    INTEGER, PARAMETER :: nbsz   =  3
      45             :    INTEGER, PARAMETER :: ibfile =  6
      46             :    INTEGER, PARAMETER :: nv48   = (2*nbsz+1)**3+48
      47             : 
      48             :    INTEGER ndiv3              ! max. number of tetrahedrons (< 6*(kpts%nkpt+1)
      49             :    INTEGER ntet               ! actual number of tetrahedrons
      50             : 
      51          27 :    REAL, ALLOCATABLE    :: vkxyz(:,:)  ! vector of kpoint generated; in cartesian representation
      52          27 :    REAL, ALLOCATABLE    :: wghtkp(:)   !   associated with k-points for BZ integration
      53          27 :    INTEGER, ALLOCATABLE :: ntetra(:,:) ! corners of the tetrahedrons
      54          27 :    REAL, ALLOCATABLE    :: voltet(:)   ! voulmes of the tetrahedrons
      55          27 :    REAL, ALLOCATABLE    :: vktet(:,:)
      56             : 
      57             :    REAL    divis(4)           ! Used to find more accurate representation of k-points
      58             :                               ! vklmn(i,kpt)/divis(i) and weights as wght(kpt)/divis(4)
      59             :    INTEGER nkstar             ! number of stars for k-points generated in full stars
      60             :    REAL    bltv(3,3)          ! cartesian Bravais lattice basis (a.u.)
      61             :    REAL    rltv(3,3)          ! reciprocal lattice basis (2\pi/a.u.)
      62             :    REAL    ccr(3,3,nop48)     ! rotation matrices in cartesian repr.
      63             :    REAL    rlsymr(3,3,nop48)  ! rotation matrices in reciprocal lattice basis representation
      64             :    REAL    talfa(3,nop48)     ! translation vector associated with (non-symmorphic)
      65             :                               ! symmetry elements in Bravais lattice representation
      66             :    INTEGER ncorn,nedge,nface  ! number of corners, faces and edges of the IBZ
      67             :    REAL    fnorm(3,mface)     ! normal vector of the planes bordering the IBZ
      68             :    REAL    fdist(mface)       ! distance vector of the planes bordering t IBZ
      69             :    REAL    cpoint(3,mface)    ! cartesian coordinates of corner points of IBZ
      70             :    REAL    xvec(3)            ! arbitrary vector lying in the IBZ
      71             : 
      72             :    INTEGER idsyst   ! crystal system identification in MDDFT programs
      73             :    INTEGER idtype   ! lattice type identification in MDDFT programs
      74             : 
      75             :    INTEGER idimens  ! number of dimensions for k-point set (2 or 3)
      76             :    INTEGER nreg     ! 1 kpoints in full BZ; 0 kpoints in irrBZ
      77             :    INTEGER nfulst   ! 1 kpoints ordered in full stars
      78             :                     !    (meaningful only for nreg =1; full BZ)
      79             :    INTEGER nbound   ! 0 no primary points on BZ boundary;
      80             :                     ! 1 with boundary points (not for BZ integration!!!)
      81             :    INTEGER ikzero   ! 0 no shift of k-points;
      82             :                     ! 1 shift of k-points for better use of sym in irrBZ
      83             :    REAL    kzero(3) ! shifting vector to bring one k-point to or 
      84             :                     ! away from (0,0,0) (for even/odd nkpt3)
      85             : 
      86             :    INTEGER i,j,k,l,idiv,mkpt,addSym,nsym
      87             :    INTEGER iofile,iokpt,kpri,ktest,kmidtet
      88             :    INTEGER idivis(3)
      89             :    LOGICAL random,trias
      90             :    REAL help(3),binv(3,3),rlsymr1(3,3),ccr1(3,3)
      91             : 
      92          27 :    random = .false.  ! do not use random tetra-points
      93             : 
      94             :    !------------------------------------------------------------
      95             :    !
      96             :    !        idsyst         idtype 
      97             :    !
      98             :    !   1  cubic          primitive
      99             :    !   2  tetragonal     body centered
     100             :    !   3  orthorhombic   face centered
     101             :    !   4  hexagonal      A-face centered
     102             :    !   5  trigonal       B-face centered
     103             :    !   6  monoclinic     C-face centered
     104             :    !   7  triclinic 
     105             :    !
     106             :    ! --->   for 2 dimensions only the following Bravais lattices exist:
     107             :    !
     108             :    !    TYPE                    EQUIVALENT 3-DIM        idsyst/idtype
     109             :    !   square               = p-tetragonal ( 1+2 axis )      2/1
     110             :    !   rectangular          = p-orthorhomb ( 1+2 axis )      3/1
     111             :    !   centered rectangular = c-face-orthorhomb( 1+2 axis)   3/6
     112             :    !   hexagonal            = p-hexagonal  ( 1+2 axis )      4/1
     113             :    !   oblique              = p-monoclinic ( 1+2 axis )      6/1
     114             :    !
     115             :    !------------------------------------------------------------
     116             : 
     117          27 :    IF(l_q) THEN
     118           0 :       trias=input%tria
     119           0 :       if (input%tria) call judft_error("tria=T not implemented for q-point generator",calledby='julia')
     120             :       !input%tria=.false.
     121             :    ENDIF
     122             :        
     123          27 :    IF (cell%latnam.EQ.'squ') THEN
     124          15 :       idsyst = 2
     125          15 :       idtype = 1
     126          15 :          IF (.not.input%film) THEN
     127          13 :             IF (abs(cell%amat(1,1)-cell%amat(3,3)) < 0.0000001) THEN
     128           0 :                idsyst = 1
     129             :                idtype = 1
     130             :             END IF
     131             :          END IF
     132             :       END IF
     133          27 :    IF (cell%latnam.EQ.'p-r') THEN
     134           0 :       idsyst = 3
     135           0 :       idtype = 1
     136             :    END IF
     137          27 :    IF ((cell%latnam.EQ.'c-b').OR.(cell%latnam.EQ.'c-r')) THEN
     138           0 :       idsyst = 3
     139           0 :       idtype = 6
     140             :    END IF
     141          27 :    IF ((cell%latnam.EQ.'hex').OR.(cell%latnam.EQ.'hx3')) THEN
     142           6 :       idsyst = 4
     143           6 :       idtype = 1
     144             :    END IF
     145          27 :    IF (cell%latnam.EQ.'obl') THEN
     146           0 :       idsyst = 6
     147           0 :       idtype = 1
     148             :    END IF
     149          27 :    IF (cell%latnam.EQ.'any') THEN
     150           6 :       CALL bravais(cell%amat,idsyst,idtype) 
     151             :    END IF
     152          27 :    nsym = sym%nop
     153          27 :    IF (input%film) nsym = sym%nop2        
     154             : 
     155             :    ! Want to make a Bandstructure?
     156          27 :    IF (banddos%ndir == -4) THEN
     157           1 :       CALL bandstr1(idsyst,idtype,cell%bmat,kpts,input,l_fillArrays,banddos)
     158           1 :       RETURN
     159             :    END IF
     160             : 
     161             :    ! Some variables we do not use
     162             : 
     163          26 :    iofile = 6
     164          26 :    iokpt  = 6
     165          26 :    kpri   = 0 ! 3
     166          26 :    ktest  = 0 ! 5
     167          26 :    kmidtet = 0
     168          26 :    nreg    = 0
     169          26 :    nfulst  = 0
     170          26 :    ikzero  = 0
     171          26 :    kzero(1) = 0.0
     172          26 :    kzero(2) = 0.0
     173          26 :    kzero(3) = 0.0 
     174          26 :    nbound  = 0
     175          26 :    IF (input%tria) THEN
     176           1 :       IF (input%film) nbound  = 1
     177             :       ! IF ((idsyst==1).AND.(idtype==1)) nbound  = 1
     178             :       ! IF ((idsyst==2).AND.(idtype==1)) nbound  = 1
     179             :       ! IF ((idsyst==3).AND.(idtype==1)) nbound  = 1
     180             :       ! IF ((idsyst==3).AND.(idtype==6)) nbound  = 1
     181             :       ! IF ((idsyst==4).AND.(idtype==1)) nbound  = 1
     182           1 :       IF (nbound == 0) random = .true.
     183             :    END IF
     184          26 :    idimens = 3
     185          26 :    IF (input%film) idimens = 2
     186             : 
     187             :    ! Lattice information
     188             : 
     189         182 :    DO j = 1, 3
     190         572 :       DO k = 1, 3
     191         234 :          bltv(j,k) = cell%amat(k,j)
     192         234 :          binv(j,k) = cell%bmat(k,j) / tpi_const
     193         234 :          rltv(j,k) = cell%bmat(k,j)
     194        3552 :          DO i = 1,nsym
     195        3474 :             rlsymr(k,j,i) = real(sym%mrot(j,k,i))
     196             :          END DO
     197             :       END DO
     198             :    END DO
     199             : 
     200          26 :    ccr = 0.0
     201         386 :    DO i = 1, nsym
     202        2546 :       DO j = 1, 3
     203        1080 :          talfa(j,i) = 0.0
     204        4320 :          DO k = 1, 3
     205        3240 :             talfa(j,i) = bltv(j,k) * sym%tau(k,i)
     206        3240 :             help(k) = 0.0
     207       14040 :             DO l = 1, 3
     208       12960 :                help(k) =  help(k) + rlsymr(l,k,i) * binv(j,l)
     209             :             END DO
     210             :          END DO
     211        7920 :          DO k = 1, 3
     212        3240 :             ccr(j,k,i) = 0.0
     213       14040 :             DO l = 1, 3
     214       12960 :                ccr(j,k,i) = ccr(j,k,i) + bltv(l,k) * help(l)
     215             :             END DO
     216             :          END DO
     217             :       END DO
     218             :    END DO
     219         746 :    DO i = 1, nsym
     220         360 :       rlsymr1(:,:) = rlsymr(:,:,i)
     221        1440 :       ccr1(:,:)    = ccr(:,:,i)
     222        2546 :       DO j = 1, 3
     223        7920 :          DO k = 1, 3
     224        3240 :             rlsymr(k,j,i) = rlsymr1(j,k)
     225        4320 :             ccr(k,j,i)    = ccr1(j,k)
     226             :          END DO
     227             :       END DO
     228             :    END DO
     229             : 
     230          26 :    IF ((.not.noco%l_ss).AND.(.not.noco%l_soc).AND.(2*nsym<nop48)) THEN
     231          17 :       IF ((input%film.AND.(.not.sym%invs2)).OR.((.not.input%film).AND.(.not.sym%invs))) THEN
     232             :          addSym = 0
     233             :          ! Note: We have to add the negative of each symmetry operation
     234             :          !       to exploit time reversal symmetry. However, if the new
     235             :          !       symmetry operation is the identity matrix it is excluded.
     236             :          !       This is the case iff it is (-Id) + a translation vector.
     237         128 :          DO i = 1, nsym
     238             :             ! This test assumes that ccr(:,:,1) is the identity matrix.
     239          69 :             IF(.NOT.ALL(ABS(ccr(:,:,1)+ccr(:,:,i)).LT.10e-10) ) THEN
     240          59 :                ccr(:,:,nsym+addSym+1 ) = -ccr(:,:,i)
     241         236 :                rlsymr(:,:,nsym+addSym+1 ) = -rlsymr(:,:,i)
     242          59 :                addSym = addSym + 1
     243             :             END IF
     244             :          END DO
     245          10 :          nsym = nsym + addSym
     246             :       END IF
     247             :    END IF
     248             : 
     249             :    ! brzone and brzone2 find the corner-points, the edges, and the
     250             :    ! faces of the irreducible wedge of the brillouin zone (IBZ).
     251             :    ! In these subroutines many special cases can occur. Due to this the very 
     252             :    ! sophisticated old routine brzone had a few bugs. The new routine
     253             :    ! brzone2 was written with a different algorithm that is slightly slower
     254             :    ! but should be more stable. To make comparisons possible the old
     255             :    ! routine is only commented out. Both routines are directly 
     256             :    ! interchangable. GM, 2016.
     257             : 
     258             : !   CALL brzone(rltv,nsym,ccr,mface,nbsz,nv48,cpoint,xvec,ncorn,nedge,nface,fnorm,fdist)
     259             : 
     260          26 :    CALL brzone2(rltv,nsym,ccr,mface,nbsz,nv48,cpoint,xvec,ncorn,nedge,nface,fnorm,fdist)
     261             : 
     262          26 :    IF (input%tria.AND.random) THEN
     263             :       ! Calculate the points for tetrahedron method
     264           1 :       mkpt = kpts%nkpt
     265           1 :       ndiv3 = 6*(mkpt+1)
     266           1 :       ALLOCATE (vkxyz(3,mkpt),wghtkp(mkpt))
     267           1 :       ALLOCATE (voltet(ndiv3),vktet(3,mkpt),ntetra(4,ndiv3))
     268          21 :       vkxyz = 0.0
     269             :       CALL kpttet(iofile,ibfile,iokpt,kpri,ktest,kmidtet,mkpt,ndiv3,&
     270             :                   nreg,nfulst,rltv,cell%omtil,nsym,ccr,mdir,mface,&
     271             :                   ncorn,nface,fdist,fnorm,cpoint,voltet,ntetra,ntet,vktet,&
     272           1 :                   kpts%nkpt,divis,vkxyz,wghtkp)
     273             :    ELSE
     274             :       ! If just the total number of k-points is given, determine 
     275             :       ! the divisions in each direction (nkpt3):
     276             : 
     277             :       ! IF (tria) THEN
     278             :       !    nkpt = nkpt/4
     279             :       !    nkpt3(:) = nkpt3(:) / 2
     280             :       ! END IF
     281         100 :       IF (sum(kpts%nkpt3).EQ.0) THEN
     282          18 :          CALL divi(kpts%nkpt,cell%bmat,input%film,sym%nop,sym%nop2,kpts%nkpt3)
     283             :       END IF
     284             : 
     285             :       ! Now calculate Monkhorst-Pack k-points:
     286          25 :       IF (kpts%nkpt3(2).EQ.0) kpts%nkpt3(2) = kpts%nkpt3(1)
     287          25 :       IF ((.not.input%film).AND.(kpts%nkpt3(3).EQ.0)) kpts%nkpt3(3) = kpts%nkpt3(2)
     288          25 :       IF (nbound.EQ.1) THEN
     289           0 :          mkpt = (2*kpts%nkpt3(1)+1)*(2*kpts%nkpt3(2)+1)
     290           0 :          IF (.not.input%film) mkpt = mkpt*(2*kpts%nkpt3(3)+1)
     291             :       ELSE
     292          25 :          mkpt = kpts%nkpt3(1)*kpts%nkpt3(2)
     293          25 :          IF (.not.input%film) mkpt = mkpt*kpts%nkpt3(3)
     294             :       END IF
     295          25 :       ALLOCATE (vkxyz(3,mkpt),wghtkp(mkpt) )
     296        6061 :       vkxyz = 0.0
     297             : 
     298             :       CALL kptmop(iofile,iokpt,kpri,ktest,idsyst,idtype,kpts%nkpt3,ikzero,kzero,&
     299             :                   rltv,bltv,nreg,nfulst,nbound,idimens,xvec,fnorm,fdist,ncorn,nface,&
     300             :                   nedge,cpoint,nsym,ccr,rlsymr,talfa,mkpt,mface,mdir,&
     301          25 :                   kpts%nkpt,divis,vkxyz,nkstar,wghtkp)
     302             :    END IF
     303             : 
     304          26 :    idivis(1) = int(divis(1)) 
     305          26 :    idivis(2) = int(divis(2)) 
     306          26 :    idivis(3) = int(divis(3)) 
     307          26 :    idiv = lcm(3,idivis)
     308          26 :    IF (idiv.GE.200) idiv = 1
     309         368 :    DO j=1,kpts%nkpt
     310         342 :       wghtkp(j) = wghtkp(j) * divis(4)
     311        1368 :       DO k = 1,3
     312        1026 :          help(k) = 0.0
     313        4446 :          DO l = 1,3
     314        4104 :             help(k) = help(k) + cell%amat(l,k) * vkxyz(l,j)
     315             :          END DO
     316             :       END DO
     317        2420 :       DO i=1,3
     318        1368 :          vkxyz(i,j) = help(i) * idiv / tpi_const
     319             :       END DO
     320             :    END DO
     321             : 
     322             :    ! if (l_q) write qpts file:
     323          26 :    IF(l_q)THEN
     324           0 :       IF(input%film) THEN
     325           0 :          CALL juDFT_error("For the case of input%film q-points generator not implemented!", calledby = "julia")
     326             :       END IF
     327             :     
     328           0 :       OPEN(113,file='qpts',form='formatted',status='new')
     329           0 :       WRITE(113,'(i5)') kpts%nkpt+1
     330           0 :       WRITE(113,8050) 0.,0.,0.
     331           0 :       DO j = 1, kpts%nkpt
     332           0 :          WRITE (113,FMT=8050) (vkxyz(i,j)/real(idiv),i=1,3)
     333             :       END DO
     334           0 :       CLOSE(113)
     335             :          !input%tria=trias
     336           0 :          RETURN
     337             :    END IF
     338             :    8050 FORMAT (2(f14.10,1x),f14.10)
     339             : 
     340             :    ! write k-points file or write data into arrays
     341          26 :    IF (l_fillArrays) THEN
     342          24 :       IF (ALLOCATED(kpts%bk)) THEN
     343          16 :          DEALLOCATE(kpts%bk)
     344             :       END IF
     345          24 :       IF (ALLOCATED(kpts%wtkpt)) THEN
     346          16 :          DEALLOCATE(kpts%wtkpt)
     347             :       END IF
     348          24 :       ALLOCATE(kpts%bk(3,kpts%nkpt),kpts%wtkpt(kpts%nkpt))
     349          24 :       IF (idiv.NE.0) kpts%posScale = REAL(idiv)
     350         702 :       DO j = 1, kpts%nkpt
     351         339 :          kpts%bk(1,j) = vkxyz(1,j)
     352         339 :          kpts%bk(2,j) = vkxyz(2,j)
     353         339 :          kpts%bk(3,j) = vkxyz(3,j)
     354         363 :          kpts%wtkpt(j) = wghtkp(j)
     355             :       END DO
     356          24 :       IF (input%tria.AND.random) THEN
     357           1 :          kpts%ntet = ntet
     358           1 :          IF (ALLOCATED(kpts%ntetra)) THEN
     359           1 :             DEALLOCATE(kpts%ntetra)
     360             :          END IF
     361           1 :          IF (ALLOCATED(kpts%voltet)) THEN
     362           1 :             DEALLOCATE(kpts%voltet)
     363             :          END IF
     364           1 :          ALLOCATE(kpts%ntetra(4,kpts%ntet))
     365           1 :          ALLOCATE(kpts%voltet(kpts%ntet))
     366          45 :          DO j = 1, ntet
     367         396 :             DO i = 1, 4
     368         220 :                kpts%ntetra(i,j) = ntetra(i,j)
     369             :             END DO
     370          45 :             kpts%voltet(j) = ABS(voltet(j))
     371             :          END DO
     372             :       END IF
     373             :    ELSE
     374           2 :       OPEN (41,file='kpts',form='formatted',status='new')
     375           2 :       IF (input%film) THEN
     376           0 :          WRITE (41,FMT=8110) kpts%nkpt,real(idiv),.false.
     377           0 :          DO j = kpts%nkpt, 1, -1
     378           0 :             WRITE (41,FMT=8040) (vkxyz(i,j),i=1,2),wghtkp(j)
     379             :          END DO
     380             :       ELSE
     381           2 :          WRITE (41,FMT=8100) kpts%nkpt,real(idiv)
     382           5 :          DO j = 1, kpts%nkpt
     383           5 :             WRITE (41,FMT=8040) (vkxyz(i,j),i=1,3),wghtkp(j)
     384             :          END DO
     385           2 :          IF (input%tria.AND.random) THEN
     386           0 :             WRITE (41,'(i5)') ntet
     387           0 :             WRITE (41,'(4(4i6,4x))') ((ntetra(i,j),i=1,4),j=1,ntet)
     388           0 :             WRITE (41,'(4f20.13)') (ABS(voltet(j)),j=1,ntet)
     389             :          END IF
     390             :       END IF
     391             :       8100 FORMAT (i5,f20.10)
     392             :       8110 FORMAT (i5,f20.10,3x,l1)
     393             :       8040 FORMAT (4f10.5)
     394           2 :       CLOSE (41)
     395             :    END IF
     396             : 
     397          26 :    DEALLOCATE (vkxyz,wghtkp)
     398          26 :    IF (input%tria.AND..not.input%film) DEALLOCATE (voltet,vktet,ntetra)
     399          27 :    RETURN
     400             : 
     401             :    CONTAINS
     402             : 
     403          26 :    INTEGER FUNCTION lcm( n, ints )
     404             :       ! Compute least common multiple (lcm) of n positive integers.
     405             :       INTEGER, INTENT(IN) :: n
     406             :       INTEGER, INTENT(IN) :: ints(n)
     407             : 
     408             :       INTEGER :: i,j,m
     409             : 
     410         104 :       IF (any(ints(1:n) <= 0)) THEN
     411             :          m = 0
     412             :       ELSE
     413         104 :          m = maxval( ints(1:n) )
     414         182 :          DO i = 1, n
     415          78 :             DO j = 1, ints(i) / 2
     416          78 :                IF (mod(m*j,ints(i)) == 0) EXIT
     417             :             END DO
     418         104 :             m = m*j
     419             :          END DO
     420             :       END IF
     421          26 :       lcm = m
     422             :       RETURN
     423             :    END FUNCTION lcm
     424             : 
     425             : END SUBROUTINE julia
     426             : 
     427             : END MODULE m_julia

Generated by: LCOV version 1.13