LCOV - code coverage report
Current view: top level - fermi - fertri.f (source / functions) Hit Total Coverage
Test: combined.info Lines: 58 107 54.2 %
Date: 2019-09-08 04:53:50 Functions: 1 1 100.0 %

          Line data    Source code
       1             :       MODULE m_fertri
       2             :       use m_juDFT
       3             : !
       4             : !     calculates fermi energy and weights using triangular method
       5             : !
       6             :       CONTAINS
       7           3 :       SUBROUTINE fertri(
       8             :      >                  input,kpts,irank,
       9           3 :      >                  ne,nkpt,jspins,zc,eig,bk,sfac,
      10             :      X                  ef,
      11           3 :      <                  seigv,w)
      12             : 
      13             :       USE m_triang
      14             :       USE m_maketetra
      15             :       USE m_tetraef
      16             :       USE m_dosef
      17             :       USE m_dosint
      18             :       USE m_doswt
      19             :       USE m_types
      20             : !     USE m_bzints
      21             : 
      22             :       IMPLICIT NONE
      23             : 
      24             :       TYPE(t_input),INTENT(IN):: input
      25             :       TYPE(t_kpts), INTENT(IN):: kpts
      26             : !     ..
      27             : !     .. Scalar Arguments ..
      28             :       INTEGER, INTENT (IN)    :: nkpt,jspins,irank
      29             :       REAL,    INTENT (IN)    :: zc,sfac
      30             :       REAL,    INTENT (OUT)   :: seigv
      31             :       REAL,    INTENT (INOUT) :: ef
      32             : !     ..
      33             : !     .. Array Arguments ..
      34             :       INTEGER, INTENT (IN)    :: ne(:,:)!(nkptd,jspd)
      35             :       REAL,    INTENT (IN)    :: bk(:,:) !(3,nkptd)
      36             :       REAL,    INTENT (OUT)   :: w(:,:,:) !(neigd,nkptd,jspd)
      37             :       REAL,    INTENT (INOUT) :: eig(:,:,:)!(neigd,nkptd,jspd)
      38             : !     ..
      39             : !     .. Local Scalars ..
      40             :       REAL chmom,ct,de,del,dez,ei,emax,emin,s,s1,workf
      41             :       REAL lb,ub,e_set
      42             :       LOGICAL film
      43             :       INTEGER i,ic,j,jsp,k,neig
      44             :       INTEGER ntria,ntetra      ! number of triangles & tetrahedrons
      45             :       REAL as                   ! total area covered by triangles
      46             : !     ..
      47             : !     .. Local Arrays ..
      48           6 :       INTEGER itria(3,2*size(w,2))  ! index of k-points that are corner points of a triangle
      49           6 :       REAL    atr(2*size(w,2))      ! area of a triangle
      50           6 :       INTEGER itetra(4,6*size(w,2)) ! ditto for tetrahedrons
      51           6 :       REAL    voltet(6*nkpt)
      52             :       INTEGER nemax(2)
      53             : !     ..
      54             : !     .. Data statements ..
      55             :       DATA de/5.0e-3/
      56             : 
      57           3 :       IF ( irank == 0 ) THEN
      58           3 :         WRITE (6,FMT=8000)
      59             :       END IF
      60             :  8000 FORMAT (/,/,10x,'linear triangular method')
      61             : c
      62           3 :       film = .true.
      63             :       CALL triang(
      64             :      >            bk,nkpt,
      65           3 :      <            itria,ntria,atr,as,film)!keep
      66             : c
      67             : c--->   clear w and set eig=-9999.9
      68           3 :       e_set = -9999.9
      69           3 :       IF (.NOT.film) e_set = 1.0e10
      70           6 :       DO jsp = 1,jspins
      71           3 :          nemax(jsp) = 0.0
      72          66 :          DO k = 1,nkpt
      73          60 :             nemax(jsp) = max0(nemax(jsp),ne(k,jsp))
      74        1140 :             DO i = 1,ne(k,jsp)
      75        1140 :                w(i,k,jsp) = 0.
      76             :             ENDDO
      77          63 :             DO i = ne(k,jsp)+1,size(w,1)
      78           0 :                w(i,k,jsp) = 0.
      79          60 :                eig(i,k,jsp) = e_set
      80             :             ENDDO
      81             :          ENDDO
      82             :       ENDDO
      83             : c
      84             : !      sfac = 2.0/real(jspins)
      85             : c
      86             : c--->   write results of triang
      87             : 
      88           3 :       IF (.not.film) THEN
      89           3 :          IF(input%l_inpXML) THEN
      90           1 :             ntetra = kpts%ntet
      91          45 :             DO j = 1, ntetra
      92          44 :                itetra(1:4,j) = kpts%ntetra(1:4,j)
      93          45 :                voltet(j) = kpts%voltet(j) / ntetra
      94             :             END DO
      95             :          ELSE
      96           2 :             IF ( irank == 0 ) THEN
      97           2 :               WRITE (6,*)  'reading tetrahedrons from file kpts'
      98             :             END IF
      99           2 :             OPEN (41,file='kpts',FORM='formatted',STATUS='old')
     100          44 :             DO i = 1, nkpt+1
     101          44 :               READ (41,*)
     102             :             ENDDO
     103           2 :             READ (41,'(i5)',ERR=66,END=66) ntetra
     104           2 :             IF (ntetra>6*nkpt)  CALL juDFT_error("ntetra > 6 nkpt"
     105           0 :      +           ,calledby ="fertri")
     106           2 :             READ (41,'(4(4i6,4x))') ((itetra(i,j),i=1,4),j=1,ntetra)
     107           2 :             READ (41,'(4f20.13)') (voltet(j),j=1,ntetra)
     108           2 :             voltet(1:ntetra) = voltet(1:ntetra) / ntetra
     109             :             GOTO 67
     110             :  66         CONTINUE                       ! no tetrahedron-information of file
     111             :             CALL make_tetra(
     112             :      >                      nkpt,bk,ntria,itria,atr,
     113           0 :      <                      ntetra,itetra,voltet)!keep
     114             :  
     115             :  67         CONTINUE                       ! tetrahedron-information read or created
     116           2 :             CLOSE(41)
     117             :          END IF
     118           6 :          lb = MINVAL(eig(:,:,:)) - 0.01
     119           3 :          ub = ef + 0.2
     120             :          CALL tetra_ef(
     121             :      >                 jspins,nkpt,
     122             :      >                 lb,ub,eig,zc,sfac,
     123             :      >                 ntetra,itetra,voltet,
     124           3 :      <                 ef,w)!keep
     125             :       ELSE
     126             : 
     127           0 :         DO i = 1,ntria
     128           0 :            atr(i) = atr(i)/as
     129             :         ENDDO
     130           0 :         IF ( irank == 0 ) THEN
     131           0 :           WRITE (6,FMT=8010) ntria,as
     132           0 :           DO i = 1,ntria
     133           0 :             WRITE (6,FMT=8020) i, (itria(j,i),j=1,3),atr(i)
     134             :           ENDDO
     135             :         END IF
     136             :  8010   FORMAT (/,10x,'triangular decomposition of brillouin zone:',/,
     137             :      +          10x,'number of triangles=',i3,/,10x,
     138             :      +          'total area of triangles=',f12.6,/,10x,
     139             :      +          'no.,corners and (normalized) area of each triangle:',/)
     140             :  8020   FORMAT (10x,i3,3x,3i3,f14.6)
     141           0 :         IF ( irank == 0 ) THEN
     142           0 :           WRITE (6,FMT=*) 'ef_hist=',ef
     143             :         END IF
     144           0 :         ei = ef
     145             : cjr     emin = -9999.9
     146           0 :         emin = +9999.9
     147           0 :         emax = -emin
     148           0 :         ic = 1
     149           0 :    90   IF (ic.GT.100) GO TO 230
     150           0 :         ic = ic + 1
     151             : c
     152             : c     results from triang are included here
     153             : c
     154             :         CALL dosint(
     155             :      >              ei,nemax,jspins,sfac,ntria,itria,atr,eig,
     156           0 :      <              ct)
     157             : c
     158           0 :         IF ( irank == 0 ) WRITE (6,FMT=*) 'ct=',ct
     159             : 
     160           0 :         IF (ct.LT.zc) THEN            ! ei < ef
     161           0 :           emin = ei
     162           0 :           ei = ei + de
     163           0 :           IF (emin.GT.emax) GO TO 90
     164           0 :         ELSEIF (ct.GT.zc) THEN        ! ei > ef
     165           0 :           emax = ei
     166           0 :           ei = ei - de
     167           0 :           IF (emin.GT.emax) GO TO 90
     168             :         ENDIF
     169           0 :         IF (ct.NE.zc) THEN
     170           0 :           IF ( irank == 0 ) WRITE (6,FMT=*) '2nd dosint'
     171             : c--->     refine ef to a value of 5 mry * (2**-20)
     172           0 :           iterate : DO i = 1, 40
     173           0 :             ei = 0.5* (emin+emax)
     174             : c
     175             :             CALL dosint(
     176             :      >               ei,nemax,jspins,sfac,ntria,itria,atr,eig,
     177           0 :      <               ct)
     178             : c
     179           0 :             IF ( irank == 0 ) WRITE (6,FMT=*) 'i=',i,', ct=',ct
     180           0 :             IF ( ct == zc ) THEN
     181             :               EXIT iterate
     182           0 :             ELSEIF ( ct > zc ) THEN
     183           0 :               emax = ei
     184             :             ELSE
     185           0 :               emin = ei
     186             :             ENDIF
     187             :           ENDDO iterate
     188             :         ENDIF
     189           0 :         ef = ei
     190           0 :         del = emax - emin
     191           0 :         dez = zc - ct
     192           0 :         workf = -13.6058*2*ef
     193           0 :         IF ( irank == 0 ) THEN
     194           0 :           WRITE (6,FMT=8030) ef,workf,del,dez
     195             :         END IF
     196             :  8030   FORMAT(/,10x,'fermi energy=',f10.5,' har',/,10x,'work function='
     197             :      +         ,f10.5,' ev',/,10x,'uncertainity in energy and weights=',
     198             :      +         2e16.6)
     199             : c
     200             : c--->   obtain dos at ef
     201             : c
     202             :         CALL dosef(
     203           0 :      >             ei,nemax,jspins,sfac,ntria,itria,atr,eig)
     204             : c
     205             : c--->   obtain weights needed for integration
     206             : c
     207             :         CALL doswt(
     208             :      >             ei,nemax,jspins,ntria,itria,atr,eig,
     209           0 :      <             w)
     210             : 
     211             :       ENDIF ! .NOT.film
     212             : c
     213             : c--->   write weights
     214             : c
     215             : c      DO 190 jsp = 1,jspins
     216             : c         neig = nemax(jsp)
     217             : c         DO 180 i = 1,neig
     218             : c            DO 170 k = 1,nkpt
     219             : c               WRITE (6,FMT=*) 'w(',i,',',k,',',jsp,')=',w(i,k,jsp)
     220             : c  170       CONTINUE
     221             : c  180    CONTINUE
     222             : c  190 CONTINUE
     223             : c
     224             : c--->   obtain sum of weights and valence eigenvalues
     225             : c
     226           3 :       s1 = 0.
     227           3 :       seigv = 0.
     228           6 :       DO 220 jsp = 1,jspins
     229           3 :          s = 0.
     230           3 :          neig = nemax(jsp)
     231          57 :          DO 210 i = 1,neig
     232        1134 :             DO 200 k = 1,nkpt
     233        1080 :                s = s + w(i,k,jsp)
     234        1080 :                seigv = seigv + w(i,k,jsp)*eig(i,k,jsp)
     235          54 :   200       CONTINUE
     236           3 :   210    CONTINUE
     237           3 :          s1 = s1 + s
     238           3 :   220 CONTINUE
     239           3 :       seigv = sfac*seigv
     240           3 :       chmom = s1 - jspins*s
     241           3 :       IF ( irank == 0 ) THEN
     242           3 :         WRITE (6,FMT=8040) seigv,s1,chmom
     243             :       END IF
     244             :  8040 FORMAT (/,10x,'sum of valence eigenvalues=',f20.6,5x,
     245             :      +       'sum of weights=',f10.6,/,10x,'moment=',f12.6)
     246           0 :       RETURN
     247             : c
     248           0 :   230 IF ( irank == 0 ) THEN
     249           0 :         WRITE (6,FMT=8050) ei,ef,emin,emax,ct,zc
     250             :       END IF
     251             :  8050 FORMAT (/,/,10x,'error fertri: initial guess of ef off by 25 mry',
     252             :      +       ' ei,ef,emin,emax,ct,zc',/,10x,6e16.7,/,10x,
     253             :      +       'check number of bands')
     254             :       CALL juDFT_error("initial guess of ef off by 25 mry",calledby
     255           0 :      +     ="fertri")
     256             :       END SUBROUTINE fertri
     257             :       END MODULE m_fertri

Generated by: LCOV version 1.13