fermie.F90 9.79 KB
Newer Older
1 2 3 4 5 6 7 8 9 10
MODULE m_fermie
  USE m_juDFT
  !-----------------------------------------------------------------------
  !     determines the fermi energy by
  !            gaussian-integration method                          c.l.fu
  !            triangular method (or tetrahedrons)
  !            or fermi-function                                    p.kurz
  !----------------------------------------------------------------------
CONTAINS
  SUBROUTINE fermie(eig_id, mpi,kpts,obsolete,&
11
       input, noco,e_min,cell,results)
12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30

    !---------------------------------------------------f--------------------
    !
    !     a fist (T=0) approximation to the fermi-energy is determined
    !     by:
    !           zelec = sum { spindg * we }
    !                       e=<ef
    !
    !     TREE STRUCTURE: fermie----sort
    !                           ----fergwt
    !                           ----fertri---triang
    !                                     |--dosint
    !                                     |--dosef -- trisrt
    !                                     +--doswt
    !                           ----ferhis---ef_newton
    !
    !-----------------------------------------------------------------------


Daniel Wortmann's avatar
Daniel Wortmann committed
31
    USE m_eig66_io, ONLY : read_eig,write_eig
32
#if defined(CPP_MPI)&&defined(CPP_NEVER)
33 34 35 36 37 38 39
    USE m_mpi_col_eigJ
#endif
    USE m_sort
    USE m_fertri
    USE m_ferhis
    USE m_fergwt
    USE m_types
40
    USE m_xmlOutput
41 42 43 44 45 46 47 48 49 50 51 52 53 54
    IMPLICIT NONE
    TYPE(t_results),INTENT(INOUT)   :: results
    TYPE(t_mpi),INTENT(IN)   :: mpi
    TYPE(t_obsolete),INTENT(IN)   :: obsolete
    TYPE(t_input),INTENT(IN)   :: input
    TYPE(t_noco),INTENT(IN)   :: noco
    TYPE(t_cell),INTENT(IN)   :: cell
    TYPE(t_kpts),INTENT(IN)   :: kpts
    !     ..
    !     .. Scalar Arguments ..
    INTEGER, INTENT (IN) :: eig_id
    REAL,INTENT(IN)      :: e_min
    !     ..
    !     .. Array Arguments ..
Daniel Wortmann's avatar
Daniel Wortmann committed
55
    !REAL,    INTENT (OUT):: w(:,:,:) !(dimension%neigd,kpts%nkpt,dimension%jspd)
56 57
    !     ..
    !     .. Local Scalars ..
58
    REAL del  ,spindg,ssc ,ws,zc,weight,efermi
59
    INTEGER i,idummy,j,jsp,k,l,n,nbands,nstef,nv,nmat,nspins
60
    INTEGER n_help,m_spins,mspin,sslice(2)
61 62 63 64 65
    !     ..
    !     .. Local Arrays ..
    !
    INTEGER, ALLOCATABLE :: idxeig(:),idxjsp(:),idxkpt(:),INDEX(:)
    REAL,    ALLOCATABLE :: e(:),eig(:,:,:),we(:)
Daniel Wortmann's avatar
Daniel Wortmann committed
66
    INTEGER ne(kpts%nkpt,SIZE(results%w_iks,3))
67
    CHARACTER(LEN=20)    :: attributes(5)
68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124

    !--- J constants
    !--- J constants

#ifdef CPP_MPI
    INCLUDE 'mpif.h'
    INTEGER, PARAMETER :: comm = MPI_COMM_SELF
    INTEGER*4 :: nv_mpi(2),idum1d(0),idum2d(0,0)
#endif

    !     ..
    !     ..
    !***********************************************************************
    !                  ABBREVIATIONS
    !
    !     eig        : array of eigenvalues 
    !     wtkpt      : list of the weights of each k-point (from inp-file)
    !     e          : linear list of the eigenvalues 
    !     we         : list of weights of the eigenvalues in e
    !     zelec      : number of electrons 
    !     spindg     : spindegeneracy (2 in nonmagnetic calculations)
    !     seigv      : weighted sum of the occupied valence eigenvalues
    !     seigsc     : weighted sum of the semi-core eigenvalues
    !     seigscv    : sum of seigv and seigsc
    !     ts         : entropy contribution to the free energy
    !
    !***********************************************************************
    !     .. Data statements ..
    DATA del/1.0e-6/
    !     ..
    n=SIZE(results%w_iks) !size of list of all eigenvalues
    ALLOCATE (idxeig(n),idxjsp(n),idxkpt(n),INDEX(n),e(n),we(n) )
    ALLOCATE (eig(SIZE(results%w_iks,1),SIZE(results%w_iks,2),SIZE(results%w_iks,3)))

    ! initiliaze e
    e = 0


    IF ( mpi%irank == 0 ) WRITE (6,FMT=8000)
8000 FORMAT (/,/,1x,'fermi energy and band-weighting factors:')
    !
    !---> READ IN EIGENVALUES
    !
    spindg = 2.0/REAL(input%jspins)
    n = 0
    results%seigsc = 0.0
    ssc = 0.0
    n_help = 0
    !
    !---> pk non-collinear
    IF (noco%l_noco) THEN
       nspins = 1
    ELSE
       nspins = input%jspins
    ENDIF
    !---> pk non-collinear
    !
125 126
    IF (mpi%irank == 0) CALL openXMLElementNoAttributes('eigenvalues')
    DO jsp = 1,nspins
127 128 129 130 131 132
       DO  k = 1,kpts%nkpt
          CALL read_eig(eig_id,k,jsp,neig=ne(k,jsp),eig=eig(:,k,jsp))
          IF ( mpi%irank == 0 ) THEN
             WRITE (6,'(a2,3f10.5,f12.6)') 'at',kpts%bk(:,k),kpts%wtkpt(k)
             WRITE (6,'(i5,a14)') ne(k,jsp),' eigenvalues :' 
             WRITE (6,'(8f12.6)') (eig(i,k,jsp),i=1,ne(k,jsp))
133
             attributes = ''
134 135 136 137 138 139
             WRITE(attributes(1),'(i0)') jsp
             WRITE(attributes(2),'(i0)') k
             WRITE(attributes(3),'(f15.8)') kpts%bk(1,k)
             WRITE(attributes(4),'(f15.8)') kpts%bk(2,k)
             WRITE(attributes(5),'(f15.8)') kpts%bk(3,k)
             CALL writeXMLElementPoly('eigenvaluesAt',(/'spin','ikpt','k_x ','k_y ','k_z '/),attributes,eig(1:ne(k,jsp),k,jsp))
140
          END IF
141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164
       END DO
    ENDDO
    !finished reading of eigenvalues
    IF (mpi%irank == 0) CALL closeXMLElement('eigenvalues')

    if (abs(input%fixed_moment)<1E-6) THEN
       !this is a standard calculation
       m_spins=1
    else
       !total moment is fixed
       m_spins=2
    end if

    do mspin=1,m_spins
       IF (m_spins    == 1) THEN
          sslice = (/1,nspins/)
       ELSE
          sslice = (/mspin,mspin/)
          nspins = 1
       ENDIF
       n = 0
       DO jsp = sslice(1),sslice(2)
          !Generate a list of energies
          DO  k = 1,kpts%nkpt
165 166 167 168 169 170 171 172 173 174
          !
          !--->          STORE EIGENVALUES AND WEIGHTS IN A LINEAR LIST. AND MEMORIZE 
          !--->          CONECTION TO THE ORIGINAL ARRAYS
          !
          DO  j = 1,ne(k,jsp)
             e(n+j) = eig(j,k,jsp)
             we(n+j) = kpts%wtkpt(k)
             idxeig(n+j) = j+n_help
             idxkpt(n+j) = k
             idxjsp(n+j) = jsp
175
          END DO
176 177
          !--->          COUNT THE NUMBER OF EIGENVALUES
          n = n + ne(k,jsp)
178 179 180
       END DO
    END DO

181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197
    CALL sort(n,e,index)

    !     Check if no deep eigenvalue is found
    IF (e_min-MINVAL(e(1:n))>1.0) THEN
       WRITE(6,*) 'WARNING: Too low eigenvalue detected:'
       WRITE(6,*) 'min E=', MINVAL(e(1:n)),' min(enpara)=',&
            &             e_min
       CALL juDFT_warn("Too low eigenvalue detected",calledby="fermi" &
            &     ,hint ="If the lowest eigenvalue is more than 1Htr below "//&
            &     "the lowest energy parameter, you probably have picked up"//&
            &     " a ghoststate")
    END IF
    !
    !---> DETERMINE EF BY SUMMING WEIGHTS
    !
    weight = input%zelec/spindg
    results%seigv = 0.0e0
198
    IF(m_spins /= 1) weight = weight/2.0  -(mspin-1.5)*input%fixed_moment
199 200 201 202 203 204 205 206
    ws = 0.0e0
    l = 0
    DO WHILE ((ws+del).LT.weight)
       l = l + 1
       IF (l.GT.n) THEN
          IF ( mpi%irank == 0 ) THEN
             WRITE (6,FMT=8010) n,ws,weight
          END IF
207 208
          CALL juDFT_error("Not enough eavefunctions",calledby="fermie")
8010      FORMAT (/,10x,'error: not enough wavefunctions.',i10,2d20.10)
209 210 211 212 213 214 215 216
       END IF
       ws = ws + we(INDEX(l))
       results%seigv = results%seigv + e(INDEX(l))*we(INDEX(l))*spindg
       !         WRITE (6,FMT='(2f10.7)') e(index(l)),we(index(l))
    END DO
    results%ef = e(INDEX(l))
    nstef = l
    zc = input%zelec
217 218 219 220 221 222 223 224 225 226 227 228
    IF(m_spins /= 1) THEN
       zc = zc/2.0-(mspin-1.5)*input%fixed_moment
       idxjsp = 1 !assume single spin in following calculations
       IF (mspin == 1) THEN
          WRITE(6,*) "Fixed total moment calculation"
          WRITE(6,*) "Moment:",input%fixed_moment
          write(6,*) "First Spin:"
       ELSE
          WRITE(6,*) "Second Spin:"
       ENDIF
    ENDIF
       
229 230 231 232 233
    IF ( mpi%irank == 0 ) WRITE (6,FMT=8020) results%ef,nstef,results%seigv,ws,results%seigsc,ssc

    !+po
    results%ts = 0.0
    !-po
234
    results%w_iks(:,:,sslice(1):sslice(2)) = 0.0
235
    results%bandgap = 0.0
236
    IF (input%gauss) THEN
237
       CALL fergwt(kpts,input,mpi,ne(:,sslice(1):sslice(2)), eig(:,:,sslice(1):sslice(2)),results%ef,results%w_iks(:,:,sslice(1):sslice(2)),results%seigv)
238
    ELSE IF (input%tria) THEN
239 240
       CALL fertri(input,kpts,mpi%irank, ne(:,sslice(1):sslice(2)),kpts%nkpt,nspins,zc,eig(:,:,sslice(1):sslice(2)),kpts%bk,spindg,&
            results%ef,results%seigv,results%w_iks(:,:,sslice(1):sslice(2)))
241
    ELSE
242 243
       CALL ferhis(input,kpts,mpi,index,idxeig,idxkpt,idxjsp, n,&
            nstef,ws,spindg,weight,e,ne(:,sslice(1):sslice(2)),we, noco,cell,results%ef,results%seigv,results%w_iks(:,:,sslice(1):sslice(2)),results)
244 245
    END IF
    results%seigscv = results%seigsc + results%seigv
246 247 248 249 250 251 252 253 254

    IF (mspin == 2) THEN
       WRITE(6,*) "Different Fermi-energies for both spins:"
       WRITE(6,"(a,f0.3,a,f0.4,a,f0.4,a,f0.4)") "Fixed Moment:" &
            ,input%fixed_moment,"   Difference(EF):",efermi," - ",results%ef,"="&
            ,efermi-results%ef
    ENDIF
    efermi = results%ef
 enddo
255
    DEALLOCATE ( idxeig,idxjsp,idxkpt,index,e,eig,we )
256 257

    attributes = ''
258
    WRITE(attributes(1),'(f20.10)') results%ef
259
    WRITE(attributes(2),'(a)') 'Htr'
260 261
    IF (mpi%irank.EQ.0) CALL writeXMLElement('FermiEnergy',(/'value','units'/),attributes(1:2))

Daniel Wortmann's avatar
Daniel Wortmann committed
262 263 264 265 266 267 268
    !Put w_iks into eig-file
    DO jsp = 1,nspins
       DO  k = 1,kpts%nkpt
          CALL write_eig(eig_id,k,jsp,w_iks=results%w_iks(:,k,jsp))
       ENDDO
    ENDDO
    
269 270
    RETURN
8020 FORMAT (/,'FERMIE:',/,&
Daniel Wortmann's avatar
Daniel Wortmann committed
271
         &       10x,'first approx. to ef    (T=0)  :',f10.6,' htr',&
272 273
         &       '   (energy of the highest occ. eigenvalue)',/,&
         &       10x,'number of occ. states  (T=0)  :',i10,/,&
Daniel Wortmann's avatar
Daniel Wortmann committed
274
         &       10x,'first approx. to seigv (T=0)  :',f10.6,' htr',/,&
275 276 277 278 279
         &       10x,'sum of weights of occ. states :',f10.6,/,&
         &       10x,'sum of semicore eigenvalues   :',f10.6,' htr',/,&
         &       10x,'sum of semicore charge        :',f10.6,' e',/)
  END SUBROUTINE fermie
     END MODULE m_fermie