ferhis.f90 11.6 KB
Newer Older
 Gregor Michalicek committed Jun 29, 2016 1 2 3 4 5 6 ``````!-------------------------------------------------------------------------------- ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany ! This file is part of FLEUR and available as free software under the conditions ! of the MIT license as expressed in the LICENSE file in more detail. !-------------------------------------------------------------------------------- `````` Markus Betzinger committed Apr 26, 2016 7 8 ``````MODULE m_ferhis CONTAINS `````` Daniel Wortmann committed Sep 10, 2018 9 10 `````` SUBROUTINE ferhis(input,kpts,mpi, index,idxeig,idxkpt,idxjsp,n,& nstef,ws,spindg,weight, e,ne,we, noco,cell,ef,seigv,w_iks,results) `````` Markus Betzinger committed Apr 26, 2016 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 `````` !*********************************************************************** ! ! This subroutine determines the fermi energy and the sum of the ! single particle eigenvalues by histogram method. ! ! ! Theory : zelec(nwd) = sum{ sum{ sum{ we * f(e) } } } ! sp e k ! ! ! seigv = sum{ sum{ sum{ w(k) * f(e) * e } ! sp e k ! ! the weights w(k) are normalized: sum{w(k)} = 1 ! k -6 ! a) 1 for kT < 10 ! we = { -1 -6 ! b){ 1+exp(e(k,nu) -ef)/kt) } for kt >=10 ! ! in case a) we choose the Fermi energy the highest ! valence state ! b) we choose as Fermi energy the arithmetric ! mean between the highest occupied and lowest ! unoccupied state if this energy difference ! Delta E <= kT, otherwise as a). ! ! stefan bl"ugel , kfa , oct 1991 ! ! free energy and extrapolation T -> 0 implemented ! (see M.J.Gillan, J.Phys.: Condens. Matter 1, ! (1989) 689-711 ) ! ! peter richard, kfa, jun. 1995 ! ! adapted to flapw7 ! ! philipp kurz, kfa, oct. 1995 ! entropy calculation changed ! ! r.pentcheva, kfa, may 1996 ! !*********************************************************************** USE m_efnewton USE m_types `````` Gregor Michalicek committed Jun 07, 2016 55 56 `````` USE m_xmlOutput USE m_constants `````` Markus Betzinger committed Apr 26, 2016 57 58 59 60 61 62 63 64 65 66 67 `````` IMPLICIT NONE TYPE(t_results),INTENT(INOUT) :: results TYPE(t_mpi),INTENT(IN) :: mpi TYPE(t_input),INTENT(IN) :: input TYPE(t_kpts),INTENT(IN) :: kpts TYPE(t_noco),INTENT(IN),OPTIONAL :: noco TYPE(t_cell),INTENT(IN),OPTIONAL :: cell ! .. ! .. Scalar Arguments .. INTEGER,INTENT(IN) :: n ,nstef REAL,INTENT(IN) :: spindg,ws,weight `````` Daniel Wortmann committed Sep 10, 2018 68 69 `````` REAL,INTENT(INOUT) :: ef,seigv REAL,INTENT(OUT) :: w_iks(:,:,:) `````` Markus Betzinger committed Apr 26, 2016 70 71 `````` ! .. ! .. Array Arguments .. `````` Daniel Wortmann committed Jan 17, 2017 72 73 74 75 76 77 78 `````` INTEGER, INTENT (IN) :: idxeig(:)!(dimension%neigd*kpts%nkpt*dimension%jspd) INTEGER, INTENT (IN) :: idxjsp(:)!(dimension%neigd*kpts%nkpt*dimension%jspd) INTEGER, INTENT (IN) :: idxkpt(:)!(dimension%neigd*kpts%nkpt*dimension%jspd) INTEGER, INTENT (IN) :: INDEX(:)!(dimension%neigd*kpts%nkpt*dimension%jspd) INTEGER, INTENT (IN) :: ne(:,:)!(kpts%nkpt,dimension%jspd) REAL, INTENT (IN) :: e(:)!(kpts%nkpt*dimension%neigd*dimension%jspd) REAL, INTENT (INOUT) :: we(:)!(kpts%nkpt*dimension%neigd*dimension%jspd) `````` Markus Betzinger committed Apr 26, 2016 79 80 81 82 83 84 85 86 87 `````` !--- J constants !--- J constants ! .. ! .. Local Scalars .. REAL,PARAMETER:: del=1.e-6 REAL :: efermi,emax,emin,entropy,fermikn,gap,& wfermi,wvals,w_below_emin,w_near_ef,tkb `````` Daniel Wortmann committed Apr 29, 2016 88 `````` INTEGER ink,inkem,j,js,k,kpt,nocc,nocst,i,nspins `````` Markus Betzinger committed Apr 26, 2016 89 90 91 `````` ! .. Local Arrays .. REAL :: qc(3) `````` Gregor Michalicek committed Jun 07, 2016 92 `````` CHARACTER(LEN=20) :: attributes(2) `````` Markus Betzinger committed Apr 26, 2016 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 `````` ! .. !*********************************************************************** !-------> ABBREVIATIONS ! ! eig : array of eigenvalues within all energy-windows ! wtkpt : list of the weights of each k-point (from inp-file) ! e : linear list of the eigenvalues within the highest ! energy-window ! we : list of weights of the eigenvalues in e ! w : array of weights (output, needed to calculate the ! new charge-density) ! zelec : number of electrons in a window ! 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 ! tkb : value of temperature (kt) broadening around fermi ! energy in htr units ! ef : fermi energy determined to obtain charge neutrality ! wfermi : weight of the occupation number for states at the ! fermi energy. ! fd : fermi dirac distribution ! fermikn : fermi dirac distribution for the k-th point ! and n-th state !********************************************************************** ! .. `````` Daniel Wortmann committed Apr 29, 2016 121 122 `````` nspins=input%jspins if (noco%l_noco) nspins=1 `````` Markus Betzinger committed Apr 26, 2016 123 124 125 126 127 128 `````` tkb=input%tkb !might be modified if we have an insulator IF ( mpi%irank == 0 ) THEN WRITE (6,FMT='(/)') WRITE (6,FMT='(''FERHIS: Fermi-Energy by histogram:'')') END IF `````` Daniel Wortmann committed Sep 10, 2018 129 `````` efermi = ef `````` Markus Betzinger committed Apr 26, 2016 130 `````` IF (nstef.LT.n) THEN `````` Daniel Wortmann committed Sep 10, 2018 131 `````` gap = e(INDEX(nstef+1)) - ef `````` Gregor Michalicek committed Mar 02, 2017 132 `````` results%bandgap = gap*hartree_to_ev_const `````` Gregor Michalicek committed Jun 07, 2016 133 134 `````` IF ( mpi%irank == 0 ) THEN attributes = '' `````` Gregor Michalicek committed Jun 09, 2016 135 136 `````` WRITE(attributes(1),'(f20.10)') gap*hartree_to_ev_const WRITE(attributes(2),'(a)') 'eV' `````` Gregor Michalicek committed Jun 07, 2016 137 138 139 `````` CALL writeXMLElement('bandgap',(/'value','units'/),attributes) WRITE (6,FMT=8050) gap END IF `````` Markus Betzinger committed Apr 26, 2016 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 `````` END IF IF ( mpi%irank == 0 ) THEN WRITE ( 6,FMT=8010) spindg* (ws-weight) END IF ! !---> DETERMINE OCCUPATION AT THE FERMI LEVEL ! wfermi = ws - weight ! -6 !======> DETERMINE FERMI ENERGY for kT >= 10 ! ! IF (tkb.GE.del) THEN ! !---> TEMPERATURE BROADENING ! IF (nstef.LT.n) THEN ! !---> STATES ABOVE EF AVAILABLE ! `````` Daniel Wortmann committed Sep 10, 2018 160 161 162 `````` ef = 0.5* (e(INDEX(nstef+1))+ef) emax = ef + 8.0*tkb emin = ef - 8.0*tkb `````` Markus Betzinger committed Apr 26, 2016 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 `````` w_near_ef = 0.0 w_below_emin = 0.0 inkem = 0 ink_loop: DO ink = 1,n IF (e(INDEX(ink)).LT.emin) THEN inkem = ink w_below_emin = w_below_emin + we(INDEX(ink)) ELSE IF (e(INDEX(ink)).GT.emax) THEN EXIT ink_loop END IF ENDDO ink_loop IF (ink>n) THEN IF ( mpi%irank == 0 ) THEN `````` Daniel Wortmann committed Apr 29, 2016 178 `````` WRITE (6,*) 'CAUTION!!! All calculated eigenvalues ', 'are below ef + 8kt.' `````` Markus Betzinger committed Apr 26, 2016 179 180 181 182 183 184 185 186 187 188 189 `````` END IF ENDIF w_near_ef = weight - w_below_emin IF (w_near_ef.GT.del) THEN ! !---> STATES BETWEEN EF-8kt AND EF+8kt AVAILABLE: !---> ADJUST FERMI-ENERGY BY NEWTON-METHOD ! nocst = ink - 1 `````` Daniel Wortmann committed Sep 10, 2018 190 `````` CALL ef_newton(n,mpi%irank, inkem,nocst,index,tkb,e, w_near_ef,ef,we) `````` Markus Betzinger committed Apr 26, 2016 191 192 `````` ! IF ( mpi%irank == 0 ) THEN `````` Daniel Wortmann committed Sep 10, 2018 193 `````` WRITE (6,FMT=8030) ef,spindg*weight, spindg*w_below_emin,spindg* (w_below_emin+w_near_ef) `````` Markus Betzinger committed Apr 26, 2016 194 195 196 197 198 199 200 201 202 `````` END IF ELSE ! !---> NO STATES BETWEEN EF-8kt AND EF+8kt AVAILABLE ! IF ( mpi%irank == 0 ) WRITE (6,FMT=8020) nocst = nstef we(INDEX(nocst)) = we(INDEX(nocst)) - wfermi `````` Daniel Wortmann committed Sep 10, 2018 203 `````` ef = efermi `````` Markus Betzinger committed Apr 26, 2016 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 `````` tkb = 0.0 END IF ELSE ! !---> NO STATES ABOVE EF AVAILABLE ! tkb = 0.0 nocst = nstef we(INDEX(nocst)) = we(INDEX(nocst)) - wfermi END IF ELSE ! !---> NO TEMPERATURE BROADENING IF tkb < del ! nocst = nstef we(INDEX(nocst)) = we(INDEX(nocst)) - wfermi END IF ! ! write(6,*) nocst,' nocst in ferhis' ! do ink = 1,nocst ! write(6,*) ink,index(ink),we(index(ink)), ! + ' ink,index(ink),we(index(ink)): weights for eigenvalues' ! end do ! ! !=======> DETERMINE OCCUPATION NUMBER AND WEIGHT OF EIGENVALUES ! FOR EACH K_POINT ! `````` Daniel Wortmann committed Sep 10, 2018 233 `````` w_iks(:,:,:) = 0.0 `````` Markus Betzinger committed Apr 26, 2016 234 235 236 `````` IF ( mpi%irank == 0 ) WRITE (6,FMT=8080) nocst DO i=1,nocst `````` Daniel Wortmann committed Sep 10, 2018 237 `````` w_iks(idxeig(INDEX(i)),idxkpt(INDEX(i)),idxjsp(INDEX(i))) = we(INDEX(i)) `````` Markus Betzinger committed Apr 26, 2016 238 239 240 241 242 243 `````` ENDDO ! !======> CHECK SUM OF VALENCE WEIGHTS ! wvals = 0.0 `````` Daniel Wortmann committed Apr 29, 2016 244 `````` DO js = 1,nspins `````` Markus Betzinger committed Apr 26, 2016 245 `````` DO k = 1,kpts%nkpt `````` Daniel Wortmann committed Sep 10, 2018 246 `````` wvals = wvals + SUM(w_iks(:ne(k,js),k,js)) `````` Markus Betzinger committed Apr 26, 2016 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 `````` ENDDO ENDDO IF ( mpi%irank == 0 ) WRITE (6,FMT=8070) wvals ! ! !=======> DETERMINE ENTROPY ! ! ! ---> formula for entropy: ! ! entropy = - two * sum wtkpt(kpt) * ! kpt ! { sum ( f(e(kpt,nu))*log(f(e(kpt,nu))) ! n ! +(1-f(e(kpt,nu)))*log(1-f(e(kpt,nu))) ) } ! ! here we have w(n,kpt,js)= spindg*wghtkp(kpt)*f(e(kpt,n)) ! entropy = 0.0 `````` Daniel Wortmann committed Apr 29, 2016 267 `````` DO js = 1,nspins `````` Markus Betzinger committed Apr 26, 2016 268 269 `````` DO kpt = 1 , kpts%nkpt DO nocc=1,ne(kpt,js) `````` Daniel Wortmann committed Sep 10, 2018 270 `````` fermikn = w_iks(nocc,kpt,js)/kpts%wtkpt(kpt) `````` Markus Betzinger committed Apr 26, 2016 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 `````` IF ( fermikn .GT. 0.0 .AND. fermikn .LT. 1.0 ) & entropy = entropy + kpts%wtkpt(kpt) * ( fermikn * LOG( fermikn) + ( 1.0 - fermikn) * LOG( 1.0 - fermikn) ) END DO END DO ENDDO entropy = -spindg*entropy results%ts = tkb*entropy IF ( mpi%irank == 0 ) WRITE (6,FMT=8060) entropy,entropy*3.0553e-6 !: boltzmann constant in htr/k ! !=======> DETERMINE SINGLE PARTICLE ENERGY ! ! `````` Daniel Wortmann committed Sep 10, 2018 287 `````` seigv = seigv+spindg*DOT_PRODUCT(e(INDEX(:nocst)),we(INDEX(:nocst))) `````` Gregor Michalicek committed Jun 16, 2016 288 289 `````` IF (mpi%irank == 0) THEN attributes = '' `````` Daniel Wortmann committed Sep 10, 2018 290 `````` WRITE(attributes(1),'(f20.10)') seigv `````` Gregor Michalicek committed Jun 16, 2016 291 292 `````` WRITE(attributes(2),'(a)') 'Htr' CALL writeXMLElement('sumValenceSingleParticleEnergies',(/'value','units'/),attributes) `````` Daniel Wortmann committed Sep 10, 2018 293 `````` WRITE (6,FMT=8040) seigv `````` Gregor Michalicek committed Jun 16, 2016 294 `````` END IF `````` Markus Betzinger committed Apr 26, 2016 295 296 297 298 299 300 301 302 303 304 305 306 307 `````` ! ! 7.12.95 r.pentcheva seigscv = seigsc + seigv will be ! calculated in fermie ! 8000 FORMAT (/,10x,'==>efrmhi: not enough wavefunctions.',i10,2e20.10) 8010 FORMAT (10x,'charge neutrality (T=0) :',f11.6,' (zero if ',& & 'the highest occ. eigenvalue is "entirely" occupied)') 8020 FORMAT (/,10x,'no eigenvalues within 8 tkb of ef',& & ' reverts to the t=0 k method.') 8030 FORMAT (/,5x,'--> new fermi energy :',f11.6,' htr',& & /,10x,'valence charge :',f11.6,' e ',/,10x,& `````` Daniel Wortmann committed Apr 29, 2016 308 309 `````` & 'actual charge blw ef-8kt :',f11.6,' e ',/,10x,& & 'actual charge blw ef+8kt :',f11.6,' e ') `````` Markus Betzinger committed Apr 26, 2016 310 311 312 313 314 315 316 317 318 319 ``````8040 FORMAT (/,10x,'sum of val. single particle energies: ',f20.10,& & ' htr',/) 8050 FORMAT (/,10x,'bandgap :',f11.6,' htr') 8060 FORMAT (10x,'entropy :',f11.6,' *kb htr/K =',& & f10.5,' htr/K') 8070 FORMAT (10x,'sum of the valence weights :',f11.6) 8080 FORMAT (10x,'number of occ. states :',i10) END SUBROUTINE ferhis END MODULE m_ferhis``````