LCOV - code coverage report
Current view: top level - fermi - tetra_ef.f (source / functions) Hit Total Coverage
Test: combined.info Lines: 101 109 92.7 %
Date: 2019-09-08 04:53:50 Functions: 1 1 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             : 
       7             :       MODULE m_tetraef
       8             :       use m_juDFT
       9             : ! -----------------------------------------------------------------------
      10             : ! This subroutine evaluates the density of states with the tetrahedron method
      11             : ! and sets the weight factors needed for the charge density for bulk systems.
      12             : ! Adapted for the FLEUR code                                          GB 2000
      13             : ! -----------------------------------------------------------------------
      14             :       CONTAINS
      15           3 :       SUBROUTINE tetra_ef(
      16             :      >                    jspins,nkpt,
      17           3 :      >                    lb,ub,eig,zc,xfac,
      18           3 :      >                    ntetra,itetra,voltet,
      19           3 :      <                    efermi,w)
      20             : c
      21             :       IMPLICIT NONE
      22             : c
      23             : C     ..Scalar Arguments ..
      24             :       INTEGER, INTENT (IN)  :: jspins,nkpt,ntetra
      25             :       REAL,    INTENT (IN)  :: lb,ub,zc,xfac
      26             :       REAL,    INTENT (OUT) :: efermi
      27             : C     ..
      28             : C     .. Array Arguments ..
      29             :       INTEGER, INTENT (IN)    :: itetra(:,:) !(4,6*nkptd)
      30             :       REAL,    INTENT (IN)    :: voltet(:) !(6*nkpt)
      31             :       REAL,    INTENT (OUT)   ::   w(:,:,:) !(neigd,nkptd,jspd)
      32             :       REAL,    INTENT (INOUT) :: eig(:,:,:) !(neigd,nkptd,jspd)
      33             : C     ..
      34             : C     .. Local Variables ..
      35             :       INTEGER i,j,jspin,neig,nk,nelec,ncr,ntet,it
      36             :       REAL    elow,dlow,eup,dup,ttt,dfermi,wgs
      37             : C     ..
      38             : C     .. Local Arrays ..
      39           6 :       REAL weight(4),ecmax(2,size(w,1))
      40           6 :       REAL wght(2,size(w,2),size(w,1)),eval(4)
      41             : C     ..
      42             : 
      43             : 
      44          57 :       DO neig = 1,size(w,1)
      45         111 :         DO jspin = 1,jspins
      46          54 :           ecmax(jspin,neig) = -1.0e25
      47        1188 :           DO nk = 1,nkpt
      48        1080 :             wght(jspin,nk,neig) = 0.0e0
      49        1080 :                w(neig,nk,jspin) = 0.0e0
      50        1080 :             IF ( eig(neig,nk,jspin).GT.ecmax(jspin,neig) )
      51         172 :      $              ecmax(jspin,neig) = eig(neig,nk,jspin)
      52             :           ENDDO
      53             :         ENDDO
      54             :       ENDDO
      55             : c
      56             : c  check for energy degeneracies in tetrahedrons
      57             : c
      58           6 :       DO jspin = 1,jspins
      59         138 :         DO ntet = 1,ntetra
      60        4887 :           DO neig = 1,size(w,1)
      61        9636 :             DO i = 1,3
      62       21384 :               DO j = i+1,4
      63       14256 :                 IF (abs(eig(neig,itetra(i,ntet),jspin) -
      64        7128 :      +                  eig(neig,itetra(j,ntet),jspin)).LT.1.0e-7) THEN
      65             :                   eig(neig,itetra(i,ntet),jspin) =
      66          24 :      +            eig(neig,itetra(i,ntet),jspin) + 1.0e-7
      67             :                   eig(neig,itetra(j,ntet),jspin) =
      68          24 :      +            eig(neig,itetra(j,ntet),jspin) - 1.0e-7
      69             :                 ENDIF
      70             :               ENDDO
      71             :             ENDDO
      72             :           ENDDO
      73             :         ENDDO
      74             :       ENDDO
      75             : c
      76             : c calculate weight factors
      77             : c
      78         135 :       DO ntet=1,ntetra
      79        4887 :         DO neig=1,size(w,1)
      80        7260 :           DO jspin=1,jspins
      81       21384 :             DO i=1,4
      82       11880 :               eval(i) = eig(neig,itetra(i,ntet),jspin)   
      83             :             ENDDO
      84        4752 :             IF (max(eval(1),eval(2),eval(3),eval(4)).LT.9999.9) THEN
      85       21384 :               DO i=1,4
      86        9504 :                 weight(i) = 1.0
      87       49896 :                 DO j=1,4
      88       47520 :                   IF (i.NE.j) weight(i) = weight(i)*(eval(j)-eval(i))
      89             :                 ENDDO
      90             :               ENDDO
      91       21384 :               DO i=1,4
      92        9504 :                 weight(i) = 6.0*voltet(ntet)/weight(i)
      93             :                 wght(jspin,itetra(i,ntet),neig) =
      94       11880 :      +          wght(jspin,itetra(i,ntet),neig) + weight(i)
      95             :               ENDDO
      96             :             ENDIF
      97             :           ENDDO
      98             :         ENDDO
      99             :       ENDDO
     100             : c
     101             : !      xfac = 2.0/jspins
     102         111 :       DO neig = 1,size(w,1)
     103        1137 :         DO nk = 1,nkpt
     104        3294 :           DO jspin = 1,jspins
     105        2160 :             wght(jspin,nk,neig)=xfac*wght(jspin,nk,neig)
     106             :           ENDDO
     107             :         ENDDO
     108             :       ENDDO
     109             : c
     110             : c---------------------------------------------------
     111             : c determine fermi energy
     112             : c---------------------------------------------------
     113             : c
     114           3 :       nelec = zc                                     ! determine # of electrons
     115           3 :       ncr   = 0                                      ! (modified gb2000)
     116             : 
     117           3 :       elow = lb                                      ! determine lower bound 
     118           3 :       dlow = ncr
     119          63 :       DO nk = 1,nkpt
     120        2223 :         DO neig = 1,size(w,1)
     121        3300 :           DO jspin = 1,jspins
     122        1080 :             ttt = elow - eig(neig,nk,jspin)
     123        1080 :               IF ( elow.GT.ecmax(jspin,neig) )
     124           0 :      +          ttt = ecmax(jspin,neig) - eig(neig,nk,jspin)
     125        1080 :               IF (ttt.LT.0.0e0) ttt = 0.0e0
     126        2160 :               dlow = dlow + wght(jspin,nk,neig)*ttt*ttt*ttt/6
     127             :           ENDDO
     128             :         ENDDO
     129             :       ENDDO
     130           3 :       IF (dlow.GT.nelec) THEN    
     131           0 :         WRITE (6,180) elow,dlow,nelec
     132             :         CALL juDFT_error("dos: valence band too high ",calledby
     133           0 :      +       ="tetra_ef")
     134             :       ENDIF
     135             :   180 FORMAT (' valence band too high ',/,
     136             :      + '  elow ',f10.5,' dlow ',f10.5,' nelec ',i5)
     137             : 
     138           3 :       it  = 0
     139           3 :       eup = ub                                      ! determine upper bound
     140           3 :  10   dup = ncr
     141          63 :       DO nk = 1,nkpt
     142        2223 :         DO neig = 1,size(w,1)
     143        3300 :           DO jspin = 1,jspins
     144        1080 :             ttt = eup - eig(neig,nk,jspin)  
     145        1080 :             IF ( eup.GT.ecmax(jspin,neig) )
     146         600 :      +         ttt = ecmax(jspin,neig) - eig(neig,nk,jspin)  
     147        1080 :             IF (ttt.LT.0.0e0) ttt = 0.0e0
     148        2160 :             dup = dup + wght(jspin,nk,neig)*ttt*ttt*ttt/6
     149             :           ENDDO
     150             :         ENDDO
     151             :       ENDDO
     152           3 :       IF ( (dup-nelec).LT.0.00001 ) THEN 
     153           0 :         eup = eup + 0.2
     154           0 :         it  = it + 1
     155           0 :         IF( it .gt. 10 ) THEN
     156           0 :           WRITE (6,200) eup,dup,nelec
     157             :           CALL juDFT_error("dos: valence band too low ",
     158           0 :      +               calledby ="tetra_ef")
     159             :         END IF
     160             :         GOTO 10
     161             :       ENDIF
     162             :   200 FORMAT (' valence band too low ',/,
     163             :      + '  eup  ',f10.5,' dup  ',f10.5,' nelec ',i5)
     164             : 
     165         102 :       DO WHILE ( (eup-elow).GT.1.0e-10 )          ! iterate for fermi-energy
     166          99 :         efermi = 0.5*(elow+eup)
     167          99 :         dfermi = real(ncr)
     168        2079 :         DO nk = 1,nkpt
     169       73359 :           DO neig = 1,size(w,1)
     170      108900 :             DO jspin = 1,jspins
     171       35640 :               ttt  =efermi-eig(neig,nk,jspin)  
     172       35640 :               IF ( efermi.GT.ecmax(jspin,neig) )
     173       18960 :      +              ttt = ecmax(jspin,neig) - eig(neig,nk,jspin)  
     174       35640 :               IF (ttt.LT.0.0e0) ttt = 0.0e0
     175       71280 :               dfermi = dfermi + wght(jspin,nk,neig)*ttt*ttt*ttt/6
     176             :             ENDDO
     177             :           ENDDO
     178             :         ENDDO
     179         102 :         IF (dfermi.GT.nelec) THEN
     180          44 :           eup = efermi
     181             :         ELSE
     182          55 :           elow = efermi
     183             :         ENDIF
     184             :       ENDDO
     185           3 :       WRITE (6,220) efermi,dfermi,nelec
     186             :   220 FORMAT (//,'>>> D O S <<<',//,'   fermi energy ',f10.5,
     187             :      +       ' dtot ',f10.5,' nelec ',i5)
     188             : c
     189             : c---------------------------------------------------
     190             : c calculate weight factors for charge density integration
     191             : c---------------------------------------------------
     192             : c
     193         135 :       DO ntet = 1,ntetra
     194        4887 :         DO neig = 1,size(w,1)
     195        7260 :           DO jspin = 1,jspins
     196       21384 :             DO i=1,4
     197       11880 :               eval(i)=eig(neig,itetra(i,ntet),jspin)
     198             :             ENDDO
     199        4752 :             IF (max(eval(1),eval(2),eval(3),eval(4)).LT.9999.9) THEN
     200             : 
     201       21384 :               DO i = 1,4
     202        9504 :                 weight(i) = 1.0
     203       47520 :                 DO j = 1,4
     204       47520 :                   IF (i.NE.j) THEN
     205       28512 :                     weight(i) = weight(i) * (eval(j) - eval(i))
     206             :                   ENDIF
     207             :                 ENDDO
     208       11880 :                 weight(i) = 6.0 * voltet(ntet) / weight(i)
     209             :               ENDDO
     210             : 
     211             :               wgs = 0.0e0
     212       21384 :               DO i = 1,4
     213        9504 :                 ttt = efermi - eval(i)
     214        9504 :                 IF (efermi.GT.ecmax(jspin,neig)) 
     215        5280 :      +                     ttt = ecmax(jspin,neig) - eval(i)
     216        9504 :                 IF ( ttt.LT.0.0e0 ) ttt = 0.0e0
     217       11880 :                 wgs = wgs + ttt**3*weight(i)
     218             :               ENDDO
     219        2376 :               wgs = wgs / 24.0
     220             : 
     221       11880 :               DO i = 1,4
     222             :                 w(neig,itetra(i,ntet),jspin) =
     223       11880 :      $          w(neig,itetra(i,ntet),jspin) + wgs
     224             :               ENDDO
     225             : 
     226             :             ENDIF
     227             :           ENDDO
     228             :         ENDDO
     229             :       ENDDO
     230             : 
     231             : c      DO jspin = 1,jspins
     232             : c        DO nk = 1,nkpt
     233             : c          DO neig = 1,size(w,1)
     234             : c            w(neig,nk,jspin) = xfac * w(neig,nk,jspin)
     235             : c          ENDDO
     236             : c        ENDDO
     237             : c      ENDDO
     238             : 
     239           3 :       RETURN
     240             :       END SUBROUTINE tetra_ef
     241             :       END MODULE m_tetraef

Generated by: LCOV version 1.13