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