LCOV - code coverage report
Current view: top level - fermi - fermie.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 88 110 80.0 %
Date: 2019-09-08 04:53:50 Functions: 1 1 100.0 %

          Line data    Source code
       1             : MODULE m_fermie
       2             :   USE m_juDFT
       3             :   !-----------------------------------------------------------------------
       4             :   !     determines the fermi energy by
       5             :   !            gaussian-integration method                          c.l.fu
       6             :   !            triangular method (or tetrahedrons)
       7             :   !            or fermi-function                                    p.kurz
       8             :   !----------------------------------------------------------------------
       9             : CONTAINS
      10         340 :   SUBROUTINE fermie(eig_id,mpi,kpts,input,noco,e_min,cell,results)
      11             : 
      12             :     !---------------------------------------------------f--------------------
      13             :     !
      14             :     !     a fist (T=0) approximation to the fermi-energy is determined
      15             :     !     by:
      16             :     !           zelec = sum { spindg * we }
      17             :     !                       e=<ef
      18             :     !
      19             :     !     TREE STRUCTURE: fermie----sort
      20             :     !                           ----fergwt
      21             :     !                           ----fertri---triang
      22             :     !                                     |--dosint
      23             :     !                                     |--dosef -- trisrt
      24             :     !                                     +--doswt
      25             :     !                           ----ferhis---ef_newton
      26             :     !
      27             :     !-----------------------------------------------------------------------
      28             : 
      29             : 
      30             :     USE m_eig66_io, ONLY : read_eig,write_eig
      31             : #if defined(CPP_MPI)&&defined(CPP_NEVER)
      32             :     USE m_mpi_col_eigJ
      33             : #endif
      34             :     USE m_sort
      35             :     USE m_fertri
      36             :     USE m_ferhis
      37             :     USE m_fergwt
      38             :     USE m_types
      39             :     USE m_xmlOutput
      40             :     IMPLICIT NONE
      41             :     TYPE(t_results), INTENT(INOUT) :: results
      42             :     TYPE(t_mpi),     INTENT(IN)    :: mpi
      43             :     TYPE(t_input),   INTENT(IN)    :: input
      44             :     TYPE(t_noco),    INTENT(IN)    :: noco
      45             :     TYPE(t_cell),    INTENT(IN)    :: cell
      46             :     TYPE(t_kpts),    INTENT(IN)    :: kpts
      47             :     !     ..
      48             :     !     .. Scalar Arguments ..
      49             :     INTEGER, INTENT (IN) :: eig_id
      50             :     REAL,INTENT(IN)      :: e_min
      51             :     !     ..
      52             :     !     .. Array Arguments ..
      53             :     !REAL,    INTENT (OUT):: w(:,:,:) !(dimension%neigd,kpts%nkpt,dimension%jspd)
      54             :     !     ..
      55             :     !     .. Local Scalars ..
      56             :     REAL del  ,spindg,ssc ,ws,zc,weight,efermi,seigv
      57             :     INTEGER i,idummy,j,jsp,k,l,n,nbands,nstef,nv,nmat,nspins
      58             :     INTEGER n_help,m_spins,mspin,sslice(2)
      59             :     !     ..
      60             :     !     .. Local Arrays ..
      61             :     !
      62         680 :     INTEGER :: idxeig(SIZE(results%w_iks)),idxjsp(SIZE(results%w_iks)),idxkpt(SIZE(results%w_iks)),INDEX(SIZE(results%w_iks))
      63        1360 :     REAL    :: e(SIZE(results%w_iks)),we(SIZE(results%w_iks))
      64         340 :     REAL,    ALLOCATABLE :: eig(:,:,:)
      65         680 :     INTEGER ne(kpts%nkpt,SIZE(results%w_iks,3))
      66             :     CHARACTER(LEN=20)    :: attributes(5)
      67             : 
      68             :     !--- J constants
      69             :     !--- J constants
      70             : 
      71             : #ifdef CPP_MPI
      72             :     INCLUDE 'mpif.h'
      73             :     INTEGER, PARAMETER :: comm = MPI_COMM_SELF
      74             :     INTEGER*4 :: nv_mpi(2),idum1d(0),idum2d(0,0)
      75             :     INTEGER ierr
      76             : #endif
      77             : 
      78             :     !     ..
      79             :     !     ..
      80             :     !***********************************************************************
      81             :     !                  ABBREVIATIONS
      82             :     !
      83             :     !     eig        : array of eigenvalues 
      84             :     !     wtkpt      : list of the weights of each k-point (from inp-file)
      85             :     !     e          : linear list of the eigenvalues 
      86             :     !     we         : list of weights of the eigenvalues in e
      87             :     !     zelec      : number of electrons 
      88             :     !     spindg     : spindegeneracy (2 in nonmagnetic calculations)
      89             :     !     seigv      : weighted sum of the occupied valence eigenvalues
      90             :     !     seigsc     : weighted sum of the semi-core eigenvalues
      91             :     !     seigscv    : sum of seigv and seigsc
      92             :     !     ts         : entropy contribution to the free energy
      93             :     !
      94             :     !***********************************************************************
      95             :     !     .. Data statements ..
      96             :     DATA del/1.0e-6/
      97             : 
      98         340 :     ALLOCATE (eig(SIZE(results%w_iks,1),SIZE(results%w_iks,2),SIZE(results%w_iks,3)))
      99             : 
     100             :     ! initiliaze e
     101       79108 :     e = 0
     102             : 
     103         340 :     IF ( mpi%irank == 0 ) WRITE (6,FMT=8000)
     104             : 8000 FORMAT (/,/,1x,'fermi energy and band-weighting factors:')
     105             :     !
     106             :     !---> READ IN EIGENVALUES
     107             :     !
     108         340 :     spindg = 2.0/REAL(input%jspins)
     109         340 :     n = 0
     110         340 :     results%seigsc = 0.0
     111         340 :     ssc = 0.0
     112         340 :     n_help = 0
     113             :     !
     114             :     !---> pk non-collinear
     115         340 :     IF (noco%l_noco) THEN
     116         248 :        nspins = 1
     117             :     ELSE
     118          92 :        nspins = input%jspins
     119             :     ENDIF
     120             :     !---> pk non-collinear
     121             :     !
     122         340 :     IF (mpi%irank == 0) CALL openXMLElementNoAttributes('eigenvalues')
     123         700 :     DO jsp = 1,nspins
     124        2192 :        DO  k = 1,kpts%nkpt
     125        1492 :           IF (mpi%irank == 0) THEN
     126         746 :              CALL read_eig(eig_id,k,jsp,neig=ne(k,jsp),eig=eig(:,k,jsp))
     127         746 :              WRITE (6,'(a2,3f10.5,f12.6)') 'at',kpts%bk(:,k),kpts%wtkpt(k)
     128         746 :              WRITE (6,'(i5,a14)') ne(k,jsp),' eigenvalues :' 
     129         746 :              WRITE (6,'(8f12.6)') (eig(i,k,jsp),i=1,ne(k,jsp))
     130         746 :              IF(.NOT.judft_was_argument("-minimalOutput")) THEN
     131        4476 :                 attributes = ''
     132         746 :                 WRITE(attributes(1),'(i0)') jsp
     133         746 :                 WRITE(attributes(2),'(i0)') k
     134         746 :                 WRITE(attributes(3),'(f15.8)') kpts%bk(1,k)
     135         746 :                 WRITE(attributes(4),'(f15.8)') kpts%bk(2,k)
     136         746 :                 WRITE(attributes(5),'(f15.8)') kpts%bk(3,k)
     137         746 :                 CALL writeXMLElementPoly('eigenvaluesAt',(/'spin','ikpt','k_x ','k_y ','k_z '/),attributes,eig(1:ne(k,jsp),k,jsp))
     138             :              END IF
     139             :           END IF
     140             : #ifdef CPP_MPI
     141        1852 :           CALL MPI_BARRIER(mpi%mpi_comm,ierr)
     142             : #endif
     143             :        END DO
     144             :     ENDDO
     145             :     !finished reading of eigenvalues
     146         340 :     IF (mpi%irank == 0) CALL closeXMLElement('eigenvalues')
     147         340 :   IF (mpi%irank == 0) THEN 
     148             : 
     149         170 :     IF (ABS(input%fixed_moment)<1E-6) THEN
     150             :        !this is a standard calculation
     151             :        m_spins=1
     152             :     else
     153             :        !total moment is fixed
     154           0 :        m_spins=2
     155             :     END IF
     156             : 
     157         170 :     results%seigv = 0.0e0
     158         340 :     do mspin=1,m_spins
     159         170 :        IF (m_spins    == 1) THEN
     160         170 :           sslice = (/1,nspins/)
     161             :        ELSE
     162           0 :           sslice = (/mspin,mspin/)
     163           0 :           nspins = 1
     164             :        ENDIF
     165         170 :        n = 0
     166         350 :        DO jsp = sslice(1),sslice(2)
     167             :           !Generate a list of energies
     168        1096 :           DO  k = 1,kpts%nkpt
     169             :              !
     170             :              !--->          STORE EIGENVALUES AND WEIGHTS IN A LINEAR LIST. AND MEMORIZE 
     171             :              !--->          CONECTION TO THE ORIGINAL ARRAYS
     172             :              !
     173       30322 :              DO  j = 1,ne(k,jsp)
     174       29576 :                 e(n+j) = eig(j,k,jsp)
     175       29576 :                 we(n+j) = kpts%wtkpt(k)
     176       29576 :                 idxeig(n+j) = j+n_help
     177       29576 :                 idxkpt(n+j) = k
     178       30322 :                 idxjsp(n+j) = jsp
     179             :              END DO
     180             :              !--->          COUNT THE NUMBER OF EIGENVALUES
     181         926 :              n = n + ne(k,jsp)
     182             :           END DO
     183             :        END DO
     184             : 
     185         170 :        CALL sort(index(:n),e)
     186             : 
     187             :        !     Check if no deep eigenvalue is found
     188         170 :        IF (e_min-MINVAL(e(1:n))>1.0) THEN
     189           0 :           WRITE(6,*) 'WARNING: Too low eigenvalue detected:'
     190           0 :           WRITE(6,*) 'min E=', MINVAL(e(1:n)),' min(enpara)=',&
     191           0 :                &             e_min
     192             :           CALL juDFT_warn("Too low eigenvalue detected",calledby="fermi" &
     193             :                &     ,hint ="If the lowest eigenvalue is more than 1Htr below "//&
     194             :                &     "the lowest energy parameter, you probably have picked up"//&
     195           0 :                &     " a ghoststate")
     196             :        END IF
     197             :        !
     198             :        !---> DETERMINE EF BY SUMMING WEIGHTS
     199             :        !
     200         170 :        weight = input%zelec/spindg
     201         170 :        seigv=0.0
     202         170 :        IF(m_spins /= 1) weight = weight/2.0  -(mspin-1.5)*input%fixed_moment
     203         170 :        ws = 0.0e0
     204         170 :        l = 0
     205       32286 :        DO WHILE ((ws+del).LT.weight)
     206       16058 :           l = l + 1
     207       16058 :           IF (l.GT.n) THEN
     208           0 :              IF ( mpi%irank == 0 ) THEN
     209           0 :                 WRITE (6,FMT=8010) n,ws,weight
     210             :              END IF
     211           0 :              CALL juDFT_error("Not enough eavefunctions",calledby="fermie")
     212             : 8010         FORMAT (/,10x,'error: not enough wavefunctions.',i10,2d20.10)
     213             :           END IF
     214       16058 :           ws = ws + we(INDEX(l))
     215       16228 :           seigv =seigv + e(INDEX(l))*we(INDEX(l))*spindg
     216             :           !         WRITE (6,FMT='(2f10.7)') e(index(l)),we(index(l))
     217             :        END DO
     218         170 :        results%ef = e(INDEX(l))
     219         170 :        nstef = l
     220         170 :        zc = input%zelec
     221         170 :        IF(m_spins /= 1) THEN
     222           0 :           zc = zc/2.0-(mspin-1.5)*input%fixed_moment
     223           0 :           idxjsp = 1 !assume single spin in following calculations
     224           0 :           IF (mspin == 1) THEN
     225           0 :              WRITE(6,*) "Fixed total moment calculation"
     226           0 :              WRITE(6,*) "Moment:",input%fixed_moment
     227           0 :              write(6,*) "First Spin:"
     228             :           ELSE
     229           0 :              WRITE(6,*) "Second Spin:"
     230             :           ENDIF
     231             :        ENDIF
     232             : 
     233         170 :        IF ( mpi%irank == 0 ) WRITE (6,FMT=8020) results%ef,nstef,seigv,ws,results%seigsc,ssc
     234             : 
     235             :        !+po
     236         170 :        results%ts = 0.0
     237             :        !-po
     238         170 :        results%w_iks(:,:,sslice(1):sslice(2)) = 0.0
     239         170 :        results%bandgap = 0.0
     240         170 :        IF (input%gauss) THEN
     241           0 :           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)
     242         170 :        ELSE IF (input%tria) THEN
     243             :           CALL fertri(input,kpts,mpi%irank, ne(:,sslice(1):sslice(2)),kpts%nkpt,nspins,zc,eig(:,:,sslice(1):sslice(2)),kpts%bk,spindg,&
     244           3 :                results%ef,results%seigv,results%w_iks(:,:,sslice(1):sslice(2)))
     245             :        ELSE
     246             :           CALL ferhis(input,kpts,mpi,index,idxeig,idxkpt,idxjsp, n,&
     247         167 :                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)
     248             :        END IF
     249         170 :        results%seigscv = results%seigsc + results%seigv
     250             : 
     251         170 :        IF (mspin == 2) THEN
     252           0 :           WRITE(6,*) "Different Fermi-energies for both spins:"
     253           0 :           WRITE(6,"(a,f0.3,a,f0.4,a,f0.4,a,f0.4)") "Fixed Moment:" &
     254           0 :                ,input%fixed_moment,"   Difference(EF):",efermi," - ",results%ef,"="&
     255           0 :                ,efermi-results%ef
     256             :        ENDIF
     257         340 :        efermi = results%ef
     258             :     enddo
     259         170 :     DEALLOCATE (eig)
     260             : 
     261        1020 :     attributes = ''
     262         170 :     WRITE(attributes(1),'(f20.10)') results%ef
     263         170 :     WRITE(attributes(2),'(a)') 'Htr'
     264         170 :     IF (mpi%irank.EQ.0) CALL writeXMLElement('FermiEnergy',(/'value','units'/),attributes(1:2))
     265             :  ENDIF
     266             : 
     267             :     !Put w_iks into eig-file
     268         700 :     DO jsp = 1,nspins
     269        2192 :        DO  k = 1,kpts%nkpt
     270        1492 :           IF (mpi%irank == 0) CALL write_eig(eig_id,k,jsp,w_iks=results%w_iks(:,k,jsp))
     271             : #ifdef CPP_MPI
     272        1852 :           CALL MPI_BARRIER(mpi%mpi_comm,ierr)
     273             : #endif
     274             :        ENDDO
     275             :     ENDDO
     276             : 
     277         340 :     RETURN
     278             : 8020 FORMAT (/,'FERMIE:',/,&
     279             :          &       10x,'first approx. to ef    (T=0)  :',f10.6,' htr',&
     280             :          &       '   (energy of the highest occ. eigenvalue)',/,&
     281             :          &       10x,'number of occ. states  (T=0)  :',i10,/,&
     282             :          &       10x,'first approx. to seigv (T=0)  :',f10.6,' htr',/,&
     283             :          &       10x,'sum of weights of occ. states :',f10.6,/,&
     284             :          &       10x,'sum of semicore eigenvalues   :',f10.6,' htr',/,&
     285             :          &       10x,'sum of semicore charge        :',f10.6,' e',/)
     286         340 :   END SUBROUTINE fermie
     287             : END MODULE m_fermie

Generated by: LCOV version 1.13