Line data Source code
1 : !!****m* ABINIT/m_occ
2 : !! NAME
3 : !! m_occ
4 : !!
5 : !! FUNCTION
6 : !! Low-level functions for occupation factors.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2008-2026 ABINIT group (XG, AF, MG)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_occ
23 :
24 : use defs_basis
25 : use m_errors
26 : use m_abicore
27 : use m_splines
28 : use m_xmpi
29 : use m_extfpmd
30 :
31 : use m_time, only : timab, cwtime, cwtime_report
32 : use m_fstrings, only : sjoin, itoa
33 : use m_rcpaw, only : rcpaw_type
34 :
35 : implicit none
36 :
37 : private
38 : !!***
39 :
40 : public :: getnel ! Compute total number of electrons from efermi or DOS
41 : public :: newocc ! Compute new occupation numbers at each k point,
42 : public :: occeig ! (occ_{k,q}(m)-occ_k(n))/(eig0_{k,q}(m)-eig0_k(n))$,
43 : public :: occ_fd ! Fermi-Dirac statistics 1 / [(exp((e - mu)/ KT) + 1]
44 : public :: occ_dfde ! Derivative of Fermi-Dirac statistics wrt e: (exp((e - mu)/ KT) / KT[(exp((e - mu)/ KT) + 1]^2
45 : public :: occ_be ! Bose-Einstein statistics 1 / [(exp((e - mu)/ KT) - 1]
46 : public :: occ_dbe ! Derivative of Bose-Einstein statistics (exp((e - mu)/ KT) / KT[(exp((e - mu)/ KT) - 1]^2
47 : public :: dos_hdr_write
48 : public :: get_fact_spin_tol_empty
49 :
50 :
51 : integer,parameter :: nptsdiv2_def=6000
52 : ! This parameter is used in init_occ_ent and getnel
53 : ! nptsdiv2 is the number of integration points, divided by 2.
54 :
55 : real(dp),parameter :: huge_tsmearinv = 1e50_dp
56 : real(dp),parameter :: maxFDarg = 500.0_dp
57 : real(dp),parameter :: maxDFDarg = 200.0_dp
58 : real(dp),parameter :: maxBEarg = 600.0_dp
59 : real(dp),parameter :: maxDBEarg = 200.0_dp
60 :
61 :
62 : contains
63 : !!***
64 :
65 : !!****f* m_abinit/getnel
66 : !! NAME
67 : !! getnel
68 : !!
69 : !! FUNCTION
70 : !! Option = 1:
71 : !! Get the total number of electrons nelect, given a trial fermienergy fermie.
72 : !! For this, compute new occupation numbers at each k point,
73 : !! from eigenenergies eigen, according to the
74 : !! smearing scheme defined by occopt (and smearing width tsmear or tphysel).
75 : !!
76 : !! Option = 2:
77 : !! Compute and output the smeared density of states, and the integrated density
78 : !! of states, then write these data
79 : !!
80 : !! Warning: this routine assumes checks have been done in the calling
81 : !! routine, and that the values of the arguments are sensible
82 : !!
83 : !! NOTE
84 : !! In order to speed the calculation, it would be easy to
85 : !! compute the entropy only when the fermi energy is well converged
86 : !!
87 : !! INPUTS
88 : !! dosdeltae= DOS delta of Energy (needed if Option=2)
89 : !! eigen(mband*nkpt*nsppol)=eigenvalues (input or init to large number), hartree
90 : !! fermie= fermi energy/ fermi energy for excited electrons if occopt = 9 (Hartree) ! CP description modified
91 : !! fermih= fermi energy for excited holes (Hartree)
92 : !! maxocc=asymptotic maximum occupation number per band
93 : !! mband=maximum number of bands
94 : !! nband(nkpt*nsppol)=number of bands at each k point
95 : !! nkpt=number of k points
96 : !! nsppol=1 for unpolarized, 2 for spin-polarized
97 : !! occopt=option for occupancies, or re-smearing scheme if dblsmr /= 0
98 : !! option=see above
99 : !! tphysel="physical" electronic temperature with FD occupations
100 : !! tsmear=smearing width (or temperature)
101 : !! unitdos=unit number of output of the DOS. Not needed if option==1
102 : !! wtk(nkpt)=k point weights
103 : !! iB1, iB2 = band min and max between which to calculate the number of electrons
104 : !! extfpmd_nbdbuf=--optional--number of bands forced to be unoccupied for extfpmd calculations
105 : !!
106 : !! OUTPUT
107 : !! doccde(mband*nkpt*nsppol)=derivative of occupancies wrt the energy for each band and k point.
108 : !! entropy= entropy associated with the smearing (adimensional)
109 : !! nelect=number of electrons per unit cell
110 : !! occ(mband*nkpt*nsppol)=occupancies for each band and k point.
111 : !!
112 : !! NOTES
113 : !! Modified beginning 23/11/2000 by MV
114 : !! Add an additional smearing on top of a FD type, in order to improve k-point
115 : !! convergence: tsmear = 0 and tphysel ~= 2.e-3 corresponds to a small (300K)
116 : !! temperature on the electrons insufficient for convergence purposes.
117 : !! Feed re-smeared "Dirac delta" to the rest of ABINIT with only one parameter,
118 : !! tphysel, which is the physical temperature.
119 : !! encorr = correction to energy for terms of order tsmear^2:
120 : !!
121 : !! $ E_{phys} = E_{free} - encorr*(E_{int}-E_{free}) + O(tsmear^3) $
122 : !!
123 : !! SOURCE
124 :
125 549521 : subroutine getnel(doccde,dosdeltae,eigen,entropy,fermie,fermih,maxocc,mband,nband,&
126 549521 : nelect,nkpt,nsppol,occ,occopt,option,tphysel,tsmear,unitdos,wtk,&
127 : iB1,iB2,extfpmd_nbdbuf,rcpaw) ! optional parameters
128 :
129 : !Arguments ------------------------------------
130 : !scalars
131 : integer,intent(in) :: mband,nkpt,nsppol,occopt,option,unitdos
132 : real(dp),intent(in) :: dosdeltae,fermie,fermih,maxocc,tphysel,tsmear
133 : real(dp),intent(out) :: entropy,nelect
134 : !arrays
135 : integer,intent(in) :: nband(nkpt*nsppol)
136 : real(dp),intent(in) :: eigen(mband*nkpt*nsppol),wtk(nkpt)
137 : real(dp),intent(out) :: doccde(mband*nkpt*nsppol)
138 : real(dp),intent(inout) :: occ(mband*nkpt*nsppol)
139 : integer, intent(in), optional:: iB1, iB2 !! CP: added optional arguments to get number of electrons between bands iB1 and iB2
140 : integer, intent(in), optional :: extfpmd_nbdbuf
141 : !! Used only when occopt = 9
142 : type(rcpaw_type),pointer,intent(inout),optional :: rcpaw
143 :
144 : !Local variables-------------------------------
145 : ! nptsdiv2 is the number of integration points, divided by 2.
146 : ! tratio = ratio tsmear/tphysel for convoluted smearing function
147 : ! save values so we can impose recalculation of smdfun when
148 : ! the smearing or electronic temperature change between datasets
149 : ! corresponds roughly to delta_FD (maxFDarg) = 1.0d-100
150 : !
151 : ! return fermi-dirac smearing function analytically
152 : ! real(dp) :: smdFD
153 : ! smdFD (tt) = 1.0_dp / (exp(-tt/2.0_dp) + exp(tt/2.0_dp))**2
154 : !scalars
155 : integer,parameter :: prtdos1=1
156 : integer :: iband,iene,ikpt,index,index_tot,index_start,isppol,nene,nptsdiv2
157 : integer :: low_band_index, high_band_index,number_of_bands,itypat,iln,isppol_cor
158 : real(dp) :: buffer,deltaene,dosdbletot,doshalftot,dostot, wk
159 : real(dp) :: enemax,enemin,enex,intdostot,limit,tsmearinv,tsmear_eff
160 : !real(dp) :: cpu, wall, gflops
161 : character(len=500) :: msg
162 : !arrays
163 549521 : real(dp),allocatable :: entfun(:,:),occfun(:,:), smdfun(:,:),xgrid(:)
164 549521 : real(dp),allocatable :: arg(:),derfun(:),dos(:),dosdble(:),doshalf(:),ent(:), intdos(:)
165 549521 : real(dp),allocatable :: occ_tmp(:),ent_tmp(:),doccde_tmp(:)
166 549521 : real(dp),allocatable :: dos_core(:),dosdble_core(:),doshalf_core(:),intdos_core(:)
167 549521 : real(dp),allocatable :: occ_tmp_core(:),doccde_tmp_core(:),arg_core(:)
168 549521 : real(dp),allocatable :: ent_core(:), derfun_core(:)
169 : ! *************************************************************************
170 :
171 : !call cwtime(cpu, wall, gflops, "start")
172 :
173 549521 : if (option/=1 .and. option/=2)then
174 0 : ABI_BUG(sjoin('Option must be either 1 or 2. It is:', itoa(option)))
175 : end if
176 :
177 : ! Initialize the occupation function and generalized entropy function,
178 : ! at the beginning, or if occopt changed
179 :
180 549521 : if (occopt==9) then
181 1000 : low_band_index = iB1
182 1000 : high_band_index = iB2
183 1000 : number_of_bands = (iB2-iB1+1)*nkpt*nsppol
184 : else
185 548521 : low_band_index = 1
186 548521 : high_band_index = nband(1)
187 8494557 : number_of_bands = sum(nband(:))
188 : end if
189 1648563 : ABI_MALLOC(occ_tmp,(number_of_bands))
190 1099042 : ABI_MALLOC(ent_tmp,(number_of_bands))
191 1099042 : ABI_MALLOC(doccde_tmp,(number_of_bands))
192 :
193 : ! Just get the number nptsdiv2 and allocate entfun, occfun, smdfun and xgrid accordingly
194 549521 : nptsdiv2 = nptsdiv2_def
195 :
196 549521 : ABI_MALLOC(entfun,(-nptsdiv2:nptsdiv2,2))
197 549521 : ABI_MALLOC(occfun,(-nptsdiv2:nptsdiv2,2))
198 549521 : ABI_MALLOC(smdfun,(-nptsdiv2:nptsdiv2,2))
199 549521 : ABI_MALLOC(xgrid,(-nptsdiv2:nptsdiv2))
200 :
201 549521 : call init_occ_ent(entfun, limit, nptsdiv2, occfun, occopt, option, smdfun, tphysel, tsmear, tsmearinv, xgrid)
202 : ! The initialisation of occfun and entfun is done
203 :
204 : !---------------------------------------------------------------------
205 :
206 : ! write(std_out,*)' getnel : debug tphysel, tsmear = ', tphysel, tsmear
207 1099042 : ABI_MALLOC(arg,(number_of_bands))
208 1099042 : ABI_MALLOC(derfun,(number_of_bands))
209 1099042 : ABI_MALLOC(ent,(number_of_bands))
210 549521 : if (option==1) then
211 : ! normal evaluation of occupations and entropy
212 :
213 : index = 0
214 : index_tot = 0
215 1181608 : do isppol=1,nsppol
216 9137468 : do ikpt=1,nkpt
217 7955860 : if (occopt == 2) high_band_index=nband(ikpt+nkpt*(isppol-1))
218 91564139 : do iband=low_band_index,high_band_index
219 83608279 : index = index + 1
220 91564139 : if (tsmear==0) then
221 262533 : arg(index) = sign(huge_tsmearinv,fermie-eigen(index_tot + iband))
222 : else
223 83345746 : arg(index)=(fermie-eigen(index_tot + iband))*tsmearinv
224 : end if
225 : end do
226 8587956 : index_tot = index_tot + nband(ikpt+nkpt*(isppol-1))
227 : end do
228 : end do
229 :
230 : ! MG TODO: This part is expensive for dense k-meshes
231 : ! Compute the values of the occupation function, and the entropy function
232 : ! Note: splfit also takes care of the points outside of the interval,
233 : ! and assign to them the value of the closest extremal point,
234 : ! which is what is needed here.
235 :
236 549512 : call splfit(xgrid, doccde_tmp, occfun, 1, arg, occ_tmp, (2*nptsdiv2+1), number_of_bands)
237 549512 : call splfit(xgrid, derfun, entfun, 0, arg, ent, (2*nptsdiv2+1), number_of_bands)
238 :
239 : ! Normalize occ and ent, and sum number of electrons and entropy
240 : ! Use different loops for nelect and entropy because bantot may be quite large in the EPH code
241 : ! when we use very dense k-meshes.
242 :
243 : ! Manage number of bands in buffer for extfpmd calculation, when extfpmd_nbdbuf not 0.
244 : ! Set occupation and entropy of buffered bands to zero.
245 549512 : if(present(extfpmd_nbdbuf)) then
246 : index=0
247 : index_tot=0
248 1178591 : do isppol=1,nsppol
249 9077567 : do ikpt=1,nkpt
250 7898976 : if (occopt == 2) high_band_index=nband(ikpt+nkpt*(isppol-1))
251 91164791 : do iband=low_band_index,high_band_index
252 83265815 : index=index+1
253 91164791 : if (iband>high_band_index-extfpmd_nbdbuf) then
254 56680 : ent(index) = zero
255 56680 : occ_tmp(index) = zero
256 : end if
257 : end do
258 8529549 : index_tot=index_tot+nband(ikpt+nkpt*(isppol-1))
259 : end do
260 : end do
261 : end if
262 :
263 549512 : nelect=zero; entropy=zero
264 549512 : index=0
265 549512 : index_tot = 0
266 1181608 : do isppol=1,nsppol
267 9137468 : do ikpt=1,nkpt
268 7955860 : wk = wtk(ikpt)
269 7955860 : if (occopt == 2) high_band_index=nband(ikpt+nkpt*(isppol-1))
270 91564139 : do iband=low_band_index,high_band_index
271 83608279 : index = index + 1
272 83608279 : ent(index) = ent(index)*maxocc
273 83608279 : occ(iband + index_tot) = occ_tmp(index)*maxocc
274 83608279 : doccde(iband + index_tot) = -doccde_tmp(index)*maxocc*tsmearinv
275 83608279 : entropy = entropy + wk*ent(index)
276 91564139 : nelect = nelect + wk*occ(iband + index_tot)
277 : end do
278 8587956 : index_tot = index_tot + nband(ikpt+nkpt*(isppol-1))
279 : end do
280 : end do
281 :
282 549512 : if(present(rcpaw)) then
283 538896 : if(associated(rcpaw)) then
284 700 : if(.not.rcpaw%frocc.or.(rcpaw%frocc.and.rcpaw%istep<=rcpaw%updateocc)) then
285 700 : rcpaw%entropy=zero
286 700 : rcpaw%nelect_core=zero
287 1400 : do itypat=1,rcpaw%ntypat
288 1400 : if(rcpaw%atm(itypat)%zcore_orig>zero) then
289 700 : rcpaw%atm(itypat)%zcore=zero
290 1400 : do isppol=1,rcpaw%atm(itypat)%nsppol
291 2100 : ABI_MALLOC(occ_tmp_core,(rcpaw%atm(itypat)%ln_size))
292 1400 : ABI_MALLOC(doccde_tmp_core,(rcpaw%atm(itypat)%ln_size))
293 1400 : ABI_MALLOC(arg_core,(rcpaw%atm(itypat)%ln_size))
294 1400 : ABI_MALLOC(derfun_core,(rcpaw%atm(itypat)%ln_size))
295 1400 : ABI_MALLOC(ent_core,(rcpaw%atm(itypat)%ln_size))
296 2800 : do iln=1,rcpaw%atm(itypat)%ln_size
297 2800 : if (tsmear==0) then
298 0 : arg_core(iln)=sign(huge_tsmearinv,fermie-rcpaw%atm(itypat)%eig(iln,isppol))
299 : else
300 2100 : arg_core(iln)=(fermie-rcpaw%atm(itypat)%eig(iln,isppol))*tsmearinv
301 : end if
302 : enddo
303 700 : call splfit(xgrid, doccde_tmp_core, occfun, 1,arg_core,occ_tmp_core,(2*nptsdiv2+1),rcpaw%atm(itypat)%ln_size)
304 700 : call splfit(xgrid, derfun_core, entfun, 0, arg_core, ent_core,(2*nptsdiv2+1),rcpaw%atm(itypat)%ln_size)
305 2800 : do iln=1,rcpaw%atm(itypat)%ln_size
306 2100 : rcpaw%atm(itypat)%occ(iln,isppol)=rcpaw%atm(itypat)%max_occ(iln,isppol)*occ_tmp_core(iln)
307 2100 : rcpaw%atm(itypat)%zcore=rcpaw%atm(itypat)%zcore+rcpaw%atm(itypat)%occ(iln,isppol)
308 2800 : rcpaw%entropy=rcpaw%entropy+ ent_core(iln)*rcpaw%atm(itypat)%max_occ(iln,isppol)*rcpaw%atm(itypat)%mult
309 : enddo
310 700 : ABI_FREE(occ_tmp_core)
311 700 : ABI_FREE(doccde_tmp_core)
312 700 : ABI_FREE(arg_core)
313 700 : ABI_FREE(derfun_core)
314 1400 : ABI_FREE(ent_core)
315 : enddo
316 700 : rcpaw%nelect_core=rcpaw%nelect_core+rcpaw%atm(itypat)%zcore*rcpaw%atm(itypat)%mult
317 : endif
318 : end do
319 700 : nelect=nelect+rcpaw%nelect_core
320 : endif
321 700 : entropy=entropy+rcpaw%entropy
322 : endif
323 : endif
324 :
325 : !write(std_out,*) ' getnel : debug wtk, occ, eigen = ', wtk, occ, eigen
326 : !write(std_out,*)xgrid(-nptsdiv2),xgrid(nptsdiv2)
327 : !write(std_out,*)'fermie',fermie
328 : !do ii=1,bantot
329 : !write(std_out,*)ii,arg(ii),doccde(ii)
330 : !end do
331 : !write(std_out,*)'eigen',eigen(:)
332 : !write(std_out,*)'arg',arg(:)
333 : !write(std_out,*)'occ',occ(:)
334 : !write(std_out,*)'nelect',nelect
335 :
336 9 : else if (option==2) then
337 : ! evaluate DOS for smearing, half smearing, and double.
338 9 : if(tsmearinv>two*Ha_eV) then
339 8 : tsmear_eff=tsmear
340 : else ! Dirty fix for high temperatures (T>0.5eV)
341 1 : tsmearinv=two*Ha_eV
342 1 : tsmear_eff=half*ev_Ha
343 : endif
344 9 : buffer=limit/tsmearinv*.5_dp
345 : ! A Similar section is present is dos_calcnwrite. Should move all DOS stuff to m_ebands
346 : ! Choose the lower and upper energies
347 1744 : enemax=maxval(eigen(1:number_of_bands))+buffer
348 1744 : enemin=minval(eigen(1:number_of_bands))-buffer
349 9 : if(present(rcpaw)) then
350 9 : if(associated(rcpaw)) then
351 0 : do itypat=1,rcpaw%ntypat
352 0 : enemin=min(enemin,rcpaw%atm(itypat)%eig(1,1)-buffer)
353 : enddo
354 : endif
355 : endif
356 :
357 : ! Extend the range to a nicer value
358 9 : enemax=0.1_dp*ceiling(enemax*10._dp)
359 9 : enemin=0.1_dp*floor(enemin*10._dp)
360 :
361 : ! Choose the energy increment
362 9 : if(abs(dosdeltae)<tol10)then
363 9 : deltaene=0.001_dp
364 : if(prtdos1>=2)deltaene=0.0005_dp ! Higher resolution possible (and wanted) for tetrahedron
365 : else
366 0 : deltaene=dosdeltae
367 : end if
368 9 : nene=nint((enemax-enemin)/deltaene)+1
369 :
370 :
371 : ! Write the header of the DOS file, and also decides the energy range and increment
372 : call dos_hdr_write(deltaene,eigen,enemax,enemin,fermie,fermih,mband,nband,nene,&
373 9 : nkpt,nsppol,occopt,prtdos1,tphysel,tsmear_eff,unitdos)
374 :
375 18 : ABI_MALLOC(dos,(number_of_bands))
376 18 : ABI_MALLOC(dosdble,(number_of_bands))
377 18 : ABI_MALLOC(doshalf,(number_of_bands))
378 18 : ABI_MALLOC(intdos,(number_of_bands))
379 :
380 21 : do isppol=1,nsppol
381 :
382 12 : if (nsppol==2) then
383 6 : if(isppol==1) write(msg,'(a,16x,a)') '#','Spin-up DOS'
384 6 : if(isppol==2) write(msg,'(2a,16x,a)') ch10,'#','Spin-dn DOS '
385 6 : call wrtout(unitdos,msg)
386 : end if
387 12 : index_start=0
388 12 : if(isppol==2)then
389 23 : do ikpt=1,nkpt
390 20 : if (occopt == 2) high_band_index = nband(ikpt + nkpt*(isppol - 1))
391 23 : index_start=index_start + high_band_index - low_band_index + 1
392 : end do
393 : end if
394 :
395 12 : enex=enemin
396 39933 : do iene=1,nene
397 :
398 : ! Compute the arguments of the dos and occupation function
399 11342918 : arg(:)=(enex-eigen(1:number_of_bands))*tsmearinv
400 :
401 39912 : call splfit(xgrid,derfun,smdfun,0,arg,dos,(2*nptsdiv2+1),number_of_bands)
402 39912 : call splfit(xgrid,derfun,occfun,0,arg,intdos,(2*nptsdiv2+1),number_of_bands)
403 :
404 : ! Also compute the dos with tsmear halved and doubled
405 11342918 : arg(:)=arg(:)*2.0_dp
406 : !call splfit(xgrid,derfun,smdfun,0,arg,doshalf,(2*nptsdiv2+1),bantot)
407 39912 : call splfit(xgrid,derfun,smdfun,0,arg,doshalf,(2*nptsdiv2+1),number_of_bands)
408 :
409 : ! Since arg was already doubled, must divide by four
410 11342918 : arg(:)=arg(:)*0.25_dp
411 : !call splfit(xgrid,derfun,smdfun,0,arg,dosdble,(2*nptsdiv2+1),bantot)
412 39912 : call splfit(xgrid,derfun,smdfun,0,arg,dosdble,(2*nptsdiv2+1),number_of_bands)
413 :
414 : ! Now, accumulate the contribution from each eigenenergy
415 39912 : dostot=zero
416 39912 : intdostot=zero
417 39912 : doshalftot=zero
418 39912 : dosdbletot=zero
419 39912 : index=index_start
420 :
421 455288 : do ikpt=1,nkpt
422 415376 : if (occopt == 2) high_band_index=nband(ikpt+nkpt*(isppol-1))
423 6869414 : do iband=low_band_index,high_band_index
424 6414126 : index=index+1
425 6414126 : dostot=dostot+wtk(ikpt)*maxocc*dos(index)*tsmearinv
426 6414126 : intdostot=intdostot+wtk(ikpt)*maxocc*intdos(index)
427 6414126 : doshalftot=doshalftot+wtk(ikpt)*maxocc*doshalf(index)*tsmearinv*2.0_dp
428 6829502 : dosdbletot=dosdbletot+wtk(ikpt)*maxocc*dosdble(index)*tsmearinv*0.5_dp
429 : end do
430 : end do
431 :
432 39912 : if(present(rcpaw)) then
433 39912 : if(associated(rcpaw)) then
434 0 : do itypat=1,rcpaw%ntypat
435 0 : if(rcpaw%atm(itypat)%zcore_orig>zero) then
436 0 : do isppol_cor=1,rcpaw%atm(itypat)%nsppol
437 0 : ABI_MALLOC(arg_core,(rcpaw%atm(itypat)%ln_size))
438 0 : ABI_MALLOC(dos_core,(rcpaw%atm(itypat)%ln_size))
439 0 : ABI_MALLOC(intdos_core,(rcpaw%atm(itypat)%ln_size))
440 0 : ABI_MALLOC(doshalf_core,(rcpaw%atm(itypat)%ln_size))
441 0 : ABI_MALLOC(dosdble_core,(rcpaw%atm(itypat)%ln_size))
442 0 : do iln=1,rcpaw%atm(itypat)%ln_size
443 0 : arg_core(iln)=(enex-rcpaw%atm(itypat)%eig(iln,isppol_cor))*tsmearinv
444 : enddo
445 0 : call splfit(xgrid, derfun_core, smdfun, 0, arg_core,dos_core,(2*nptsdiv2+1),rcpaw%atm(itypat)%ln_size)
446 0 : call splfit(xgrid, derfun_core, occfun, 0, arg_core,intdos_core,(2*nptsdiv2+1),rcpaw%atm(itypat)%ln_size)
447 0 : arg_core=arg_core*2.0_dp
448 0 : call splfit(xgrid, derfun_core, smdfun,0, arg_core,doshalf_core,(2*nptsdiv2+1),rcpaw%atm(itypat)%ln_size)
449 0 : arg_core=arg_core*0.25_dp
450 0 : call splfit(xgrid, derfun_core, smdfun,0,arg_core,dosdble_core,(2*nptsdiv2+1),rcpaw%atm(itypat)%ln_size)
451 0 : do iln=1,rcpaw%atm(itypat)%ln_size
452 : dostot=dostot+rcpaw%atm(itypat)%max_occ(iln,isppol_cor)*dos_core(iln)*tsmearinv*rcpaw%atm(itypat)%mult/&
453 0 : nsppol
454 : intdostot=intdostot+rcpaw%atm(itypat)%max_occ(iln,isppol_cor)*intdos_core(iln)*rcpaw%atm(itypat)%mult/&
455 0 : nsppol
456 : doshalftot=doshalftot+rcpaw%atm(itypat)%max_occ(iln,isppol_cor)*doshalf_core(iln)*tsmearinv*2.0_dp*&
457 0 : rcpaw%atm(itypat)%mult/nsppol
458 : dosdbletot=dosdbletot+rcpaw%atm(itypat)%max_occ(iln,isppol_cor)*dosdble_core(iln)*tsmearinv*0.5_dp*&
459 0 : rcpaw%atm(itypat)%mult/nsppol
460 : enddo
461 0 : ABI_FREE(arg_core)
462 0 : ABI_FREE(dos_core)
463 0 : ABI_FREE(intdos_core)
464 0 : ABI_FREE(doshalf_core)
465 0 : ABI_FREE(dosdble_core)
466 : enddo
467 : endif
468 : enddo
469 : endif
470 : endif
471 :
472 : ! Print the data for this energy
473 39912 : write(unitdos, '(f8.3,2f14.6,2f14.3)' )enex,dostot,intdostot,doshalftot,dosdbletot
474 :
475 39924 : enex=enex+deltaene
476 : end do ! iene
477 : end do ! isppol
478 :
479 9 : ABI_FREE(dos)
480 9 : ABI_FREE(dosdble)
481 9 : ABI_FREE(doshalf)
482 9 : ABI_FREE(intdos)
483 :
484 : ! MG: It does not make sense to close the unit here since the routines
485 : ! did not open the file here!
486 : ! Close the DOS file
487 9 : close(unitdos)
488 : end if
489 :
490 549521 : ABI_FREE(arg)
491 549521 : ABI_FREE(derfun)
492 549521 : ABI_FREE(ent)
493 549521 : ABI_FREE(entfun)
494 549521 : ABI_FREE(occfun)
495 549521 : ABI_FREE(smdfun)
496 549521 : ABI_FREE(xgrid)
497 549521 : ABI_FREE(occ_tmp)
498 549521 : ABI_FREE(doccde_tmp)
499 549521 : ABI_FREE(ent_tmp)
500 :
501 : !call cwtime_report(" getnel", cpu, wall, gflops, end_str=ch10)
502 :
503 549521 : end subroutine getnel
504 : !!***
505 :
506 : !!****f* m_occ/newocc
507 : !! NAME
508 : !! newocc
509 : !!
510 : !! FUNCTION
511 : !! Compute new occupation numbers at each k point,
512 : !! from eigenenergies eigen, according to the
513 : !! smearing scheme defined by occopt (smearing width tsmear and
514 : !! physical temperature tphysel),
515 : !! with the constraint of number of valence electrons per unit cell nelect.
516 : !!
517 : !! INPUTS
518 : !! eigen(mband*nkpt*nsppol)=eigenvalues (input or init to large number), hartree
519 : !! spinmagntarget=if differ from -99.99_dp, fix the magnetic moment (in Bohr magneton)
520 : !! mband=maximum number of bands
521 : !! nband(nkpt)=number of bands at each k point
522 : !! nelect=number of electrons per unit cell
523 : !! ne_qFD, nh_qFD=number of thermalized excited electrons (resp. holes) in bands > ivalence (resp. <= ivalence)
524 : !! ivalence= band index of the last valence band
525 : !! nkpt=number of k points
526 : !! nspinor=number of spinorial components of the wavefunctions
527 : !! nsppol=1 for unpolarized, 2 for spin-polarized
528 : !! occopt=option for occupancies
529 : !! prtstm=optional, might govern the band-by-band decomposition of the stm charge density
530 : !! prtvol=control print volume and debugging output
531 : !! stmbias= optional, if non-zero, compute occupation numbers for STM (non-zero around the Fermi energy)
532 : !! NOTE: in this case, only fermie and occ are meaningful outputs.
533 : !! extfpmd <type(extfpmd_type)>=--optional--extended first-principles molecular dynamics type
534 : !! tphysel="physical" electronic temperature with FD occupations
535 : !! tsmear=smearing width (or temperature)
536 : !! wtk(nkpt)=k point weights
537 : !!
538 : !! OUTPUT
539 : !! doccde(maxval(nband(:))*nkpt*nsppol)=derivative of occupancies wrt
540 : !! the energy for each band and k point
541 : !! entropy= entropy associated with the smearing (adimensional)
542 : !! fermie= fermi energy (Hartree)/fermi level for thermalized excited electrons in bands > ivalence when occopt=9
543 : !! fermih= fermi level for thermalized excited holes in bands <= ivalence
544 : !! occ(maxval(nband(:))*nkpt*nsppol)=occupancies for each band and k point
545 : !!
546 : !! SOURCE
547 :
548 14188 : subroutine newocc(doccde,eigen,entropy,fermie,fermih,ivalence,spinmagntarget,mband,nband,&
549 14188 : nelect,ne_qFD,nh_qFD,nkpt,nspinor,nsppol,occ,occopt,prtvol,tphysel,tsmear,wtk,&
550 : prtstm,stmbias,extfpmd,rcpaw) ! Optional argument
551 :
552 : !Arguments ------------------------------------
553 : !scalars
554 : integer,intent(in) :: ivalence,mband,nkpt,nspinor,nsppol,occopt,prtvol
555 : integer,intent(in),optional :: prtstm
556 : real(dp),intent(in) :: spinmagntarget,nelect,tphysel,tsmear,ne_qFD, nh_qFD
557 : real(dp),intent(in),optional :: stmbias
558 : real(dp),intent(out) :: entropy,fermie,fermih
559 : type(extfpmd_type),pointer,intent(inout),optional :: extfpmd
560 : type(rcpaw_type),pointer,intent(inout),optional :: rcpaw
561 : !arrays
562 : integer,intent(in) :: nband(nkpt*nsppol)
563 : real(dp),intent(in) :: eigen(mband*nkpt*nsppol),wtk(nkpt)
564 : real(dp),intent(out) :: doccde(mband*nkpt*nsppol)
565 : real(dp),intent(inout) :: occ(mband*nkpt*nsppol)
566 :
567 : !Local variables-------------------------------
568 : integer,parameter :: niter_max=120,nkpt_max=2,fake_unit=-666,option1=1
569 : integer :: cnt,cnt2,cnt3,ib,iban,ibantot,ii,ik,ikpt,is,isppol,nban,nkpt_eff,sign
570 : integer :: extfpmd_nbdbuf=0,itypat
571 14188 : integer,allocatable :: nbandt(:)
572 : real(dp),parameter :: tol = tol14
573 : type(rcpaw_type),pointer :: rcpaw_getnel => null()
574 : !real(dp),parameter :: tol = tol10
575 : real(dp) :: dosdeltae,entropy_tmp,fermie_hi,fermie_lo,fermie_mid,fermie_mid_tmp
576 : real(dp) :: fermih_lo,fermih_mid,fermih_hi
577 : real(dp) :: fermie_biased,maxocc,rcpaw_nelect
578 : real(dp) :: nelect_tmp,nelecthi,nelectlo,nelectmid,nelect_biased
579 : real(dp) :: nholeshi,nholeslo,nholesmid
580 : real(dp) :: entropyet(2),fermie_hit(2),fermie_lot(2),fermie_midt(2),nelecthit(2)
581 : real(dp) :: nelectlot(2),nelectt(2),tsec(2)
582 : real(dp) :: entropye, entropyh
583 14188 : real(dp),allocatable :: doccdet(:),eigent(:),occt(:)
584 : character(len=500) :: msg
585 : logical:: not_enough_bands=.false.
586 : ! *************************************************************************
587 :
588 : DBG_ENTER("COLL")
589 :
590 14188 : call timab(74,1,tsec)
591 :
592 14188 : rcpaw_nelect=zero
593 14188 : if(present(rcpaw)) then
594 13651 : if(associated(rcpaw)) then
595 21 : rcpaw_getnel=>rcpaw
596 21 : rcpaw_nelect=rcpaw%nelect_core
597 42 : do itypat=1,rcpaw%ntypat
598 42 : if(rcpaw%atm(itypat)%zcore_orig>zero) then
599 126 : rcpaw%atm(itypat)%occ_res=-rcpaw%atm(itypat)%occ
600 : endif
601 : enddo
602 : endif
603 : endif
604 :
605 : ! Here treat the case where occopt does not correspond to a metallic occupation scheme
606 14188 : if (occopt < 3 .or. occopt > 9) then
607 0 : ABI_BUG(sjoin(' occopt= ',itoa(occopt),', a value not allowed in newocc.'))
608 : end if
609 :
610 : ! Check whether nband is a constant for all k point and spin-pol
611 30187 : do isppol=1,nsppol
612 220744 : do ikpt=1,nkpt
613 206556 : if(nband(ikpt+(isppol-1)*nkpt)/=nband(1)) then
614 : write(msg,'(3a,i0,a,i0,a,i0,a)')&
615 0 : 'The number of bands must be the same for all k-points ',ch10,&
616 0 : 'but nband(1)= ',nband(1),' is different of nband(',ikpt+(isppol-1)*nkpt,') = ',nband(ikpt+(isppol-1)*nkpt),'.'
617 0 : ABI_BUG(msg)
618 : end if
619 : end do
620 : end do
621 :
622 : ! Check whether nelect is strictly positive
623 14188 : if (nelect <= zero) then
624 : write(msg,'(3a,es16.8,a)')&
625 0 : 'nelect must be a positive number, while ',ch10, 'the calling routine asks nelect= ',nelect,'.'
626 0 : ABI_BUG(msg)
627 : end if
628 :
629 : ! Check whether the number of holes and electrons if positive
630 14188 : if (occopt == 9) then
631 10 : if ( (ne_qFD < zero) .or. (nh_qFD < zero) ) then
632 : write(msg,'(3a,es16.8,a,es16.8,a)')&
633 0 : & 'ne_qFD or nh_qFD must be positive numbers, while ',ch10,&
634 0 : & 'the calling routine asks ne_qFD= ',ne_qFD,' and nh_qFD= ',nh_qFD, '.'
635 0 : ABI_BUG(msg)
636 : end if
637 : end if
638 :
639 14188 : maxocc = two / (nsppol * nspinor)
640 :
641 : ! Check whether nelect is coherent with nband (nband(1) is enough,
642 : ! since it was checked that nband is independent of k-point and spin-pol
643 14188 : if (nelect > nband(1) * nsppol * maxocc+rcpaw_nelect) then
644 : write(msg,'(3a,es16.8,a,i0,a,es16.8,a)' )&
645 0 : 'nelect must be smaller than nband*maxocc, while ',ch10,&
646 0 : 'the calling routine gives nelect= ',nelect,', nband= ',nband(1),' and maxocc= ',maxocc,'.'
647 0 : ABI_BUG(msg)
648 : end if
649 :
650 : ! Providing additional checks to ensure that there are enough valence and conduction bands
651 : ! to accomodate ne_qFD and nh_qFD
652 14188 : if( occopt==9 .and. ne_qFD > (nband(1)-ivalence)*nsppol*maxocc )then
653 0 : write(msg,'(a,es16.8,2a,es16.8,a)') 'ne_qFD = ', ne_qFD ,ch10, &
654 0 : & 'must be smaller than (nband-ivalence)*maxocc*nsppol = ', &
655 0 : & (nband(1)-ivalence)*nsppol*maxocc,'.'
656 0 : ABI_BUG(msg)
657 14188 : else if( occopt==9 .and. (nh_qFD > ivalence*nsppol*maxocc .or. &
658 : & nelect - nh_qFD > ivalence*nsppol*maxocc ) )then
659 0 : write(msg,'(a,es16.8,2a,es16.8,2a,es16.8,a)') 'nh_qFD = ', nh_qFD ,ch10, &
660 0 : & 'and nelect-nh_qFD = ', nelect - nh_qFD,ch10, ' must be smaller than ivalence*maxocc*nsppol = ', &
661 0 : & ivalence*nsppol*maxocc,'.'
662 0 : ABI_BUG(msg)
663 : end if
664 :
665 : ! Set extfpmd band buffer if needed
666 14188 : if(present(extfpmd)) then
667 13675 : if(associated(extfpmd)) then
668 90 : extfpmd_nbdbuf=extfpmd%nbdbuf
669 : end if
670 : end if
671 :
672 : ! Use bisection algorithm to find fermi energy
673 : ! This choice is due to the fact that it will always give sensible
674 : ! result (because the answer is bounded, even if the smearing function
675 : ! is non-monotonic (which is the case for occopt=4 or 6)
676 : ! Might speed up it, if needed !
677 :
678 : ! Lowest and largest trial fermi energies, and corresponding number of electrons
679 : ! They are obtained from the smallest or largest eigenenergy, plus a range of
680 : ! energy that allows for complete occupation of all bands, or, on the opposite,
681 : ! for zero occupation of all bands (see getnel.f)
682 :
683 14188 : dosdeltae = zero ! the DOS is not computed, with option=1
684 2066602 : fermie_lo = minval(eigen(1:nband(1)*nkpt*nsppol)) - 6.001_dp * tsmear ! fermi_lo ->fermie_lo
685 14188 : if (occopt == 3 .or. occopt==9) fermie_lo = fermie_lo - 24.0_dp * tsmear
686 14188 : if(occopt==9) fermih_lo = fermie_lo ! Take into account holes
687 :
688 14188 : if(occopt >= 3 .and. occopt <= 8) then
689 : call getnel(doccde,dosdeltae,eigen,entropye,fermie_lo,fermie_lo,maxocc,mband,nband,&
690 : & nelectlo,nkpt,nsppol,occ,occopt,option1,tphysel,tsmear,fake_unit,wtk,1,nband(1),&
691 14178 : & extfpmd_nbdbuf=extfpmd_nbdbuf,rcpaw=rcpaw_getnel)
692 10 : else if (occopt == 9) then
693 : call getnel(doccde,dosdeltae,eigen,entropye,fermie_lo,fermie_lo,maxocc,mband,nband,&
694 10 : & nelectlo,nkpt,nsppol,occ,occopt,option1,tphysel,tsmear,fake_unit,wtk, ivalence+1, nband(1)) ! Excited electrons
695 : call getnel(doccde,dosdeltae,eigen,entropyh,fermih_lo,fermih_lo,maxocc,mband,nband,&
696 10 : & nholeslo,nkpt,nsppol,occ,occopt,option1,tphysel,tsmear,fake_unit,wtk,1, ivalence)
697 : end if
698 :
699 2066602 : fermie_hi = maxval(eigen(1:nband(1)*nkpt*nsppol)) + 6.001_dp * tsmear
700 : ! Safety value
701 14188 : fermie_hi = min(fermie_hi, 1.e6_dp)
702 14188 : if(occopt == 3 .or. occopt == 9) fermie_hi = fermie_hi + 24.0_dp * tsmear
703 14188 : if(occopt == 9) fermih_hi=fermie_hi
704 :
705 14188 : if (occopt >= 3 .and. occopt <= 8) then
706 : call getnel(doccde,dosdeltae,eigen,entropye,fermie_hi,fermie_hi,maxocc,mband,nband,&
707 : & nelecthi,nkpt,nsppol,occ,occopt,option1,tphysel,tsmear,fake_unit,wtk,1,nband(1),&
708 14178 : & extfpmd_nbdbuf=extfpmd_nbdbuf,rcpaw=rcpaw_getnel)
709 10 : else if (occopt == 9) then
710 : call getnel(doccde,dosdeltae,eigen,entropye,fermie_hi,fermie_hi,maxocc,mband,nband,&
711 10 : & nelecthi,nkpt,nsppol,occ,occopt,option1,tphysel,tsmear,fake_unit,wtk, ivalence+1, nband(1)) ! Excited electrons
712 10 : fermih_hi=fermie_hi
713 : call getnel(doccde,dosdeltae,eigen,entropyh,fermih_hi,fermih_hi,maxocc,mband,nband,&
714 10 : & nholeshi,nkpt,nsppol,occ,occopt,option1,tphysel,tsmear,fake_unit,wtk,1, ivalence)
715 : end if
716 :
717 : ! Compute the number of free electrons with corresponding chemical
718 : ! potential and add to nelect bounds.
719 14188 : if(present(extfpmd)) then
720 13675 : if(associated(extfpmd)) then
721 90 : call extfpmd%compute_nelect(fermie_lo,nband,nelectlo,nkpt,nspinor,nsppol,wtk)
722 90 : call extfpmd%compute_nelect(fermie_hi,nband,nelecthi,nkpt,nspinor,nsppol,wtk)
723 : end if
724 : end if
725 :
726 : !Prepare fixed moment calculation
727 14188 : if(abs(spinmagntarget+99.99_dp)>1.0d-10)then
728 132 : if (occopt==9)then
729 0 : ABI_ERROR('occopt=9 and spinmagntarget not implemented.')
730 : end if
731 132 : sign = 1
732 396 : do is = 1, nsppol
733 264 : fermie_hit(is) = fermie_hi
734 264 : fermie_lot(is) = fermie_lo
735 264 : nelectt(is) = half*(nelect+sign*spinmagntarget)
736 264 : sign = -sign
737 264 : nelecthit(is) = nelecthi
738 396 : nelectlot(is) = nelectlo
739 : end do
740 : end if
741 :
742 :
743 : ! If the target nelect is not between nelectlo and nelecthi, exit
744 14188 : if ((nelect < nelectlo - tol8 .or. nelect > nelecthi + tol8) .and. (occopt <= 8)) then
745 0 : not_enough_bands = .true.
746 0 : write(msg, '(a,a,a,a,d16.8,a,a,d16.8,a,d16.8,a,a,d16.8,a,d16.8)') ch10,&
747 0 : ' newocc: ',ch10,&
748 0 : ' The calling routine gives nelect= ',nelect,ch10,&
749 0 : ' The lowest bound is ',fermie_lo,', with nelect=',nelectlo,ch10,&
750 0 : ' The highest bound is ',fermie_hi,', with nelect=',nelecthi
751 0 : call wrtout(std_out, msg)
752 0 : ABI_BUG(msg)
753 : end if
754 :
755 14188 : if( occopt==9 ) then
756 10 : if ((nelect-nh_qFD)<nholeslo .or. (nelect-nh_qFD)>nholeshi) then
757 0 : not_enough_bands = .true.
758 0 : write(msg,'(a,a,a,d16.8,a,a,d16.8,a,d16.8,a)') 'newocc : ',ch10, &
759 0 : 'The calling routine gives nelect-nh_qFD = ', nelect-nh_qFD, ch10, &
760 0 : 'The lowest (highest resp.) bound for nelect-nh_qFD is ', nholeslo, ' ( ', nholeshi, ' ).'
761 0 : ABI_BUG(msg)
762 : endif
763 10 : if ((ne_qFD < nelectlo) .or. (ne_qFD > nelecthi) ) then
764 0 : not_enough_bands = .true.
765 0 : write(msg,'(a,a,a,d16.8,a,a,d16.8,a,d16.8,a)') 'newocc : ',ch10, &
766 0 : 'The calling routine gives ne_qFD = ', ne_qFD, ch10, 'The lowest (highest resp.) bound for ne_qFD are ',&
767 0 : nelectlo, ' ( ', nelecthi, ' ) .'
768 0 : ABI_BUG(msg)
769 : endif
770 :
771 10 : if (not_enough_bands) then
772 : write(msg, '(11a)' )&
773 0 : 'In order to get the right number of carriers,',ch10,&
774 0 : 'it seems that the Fermi energies must be outside the range',ch10,&
775 0 : 'of eigenenergies, plus 6 or 30 times the smearing, which is strange.',ch10,&
776 0 : 'It might be that your number of bands (nband) corresponds to the strictly',ch10,&
777 0 : 'minimum number of bands to accomodate your electrons (so, OK for an insulator),',ch10,&
778 0 : 'while you are trying to describe a metal. In this case, increase nband, otherwise ...'
779 0 : ABI_BUG(msg)
780 : end if
781 : end if
782 :
783 14188 : if( abs(spinmagntarget+99.99_dp) < tol10) then
784 :
785 : ! Usual bisection loop
786 511020 : do ii=1,niter_max
787 511020 : fermie_mid = (fermie_hi + fermie_lo) * half
788 511020 : if (occopt == 9) fermih_mid=(fermih_hi+fermih_lo)*half
789 : ! Produce nelectmid from fermimid
790 511020 : if (occopt /= 9) then
791 :
792 : call getnel(doccde,dosdeltae,eigen,entropye,fermie_mid,fermie_mid,maxocc,mband,nband,&
793 : nelectmid,nkpt,nsppol,occ,occopt,option1,tphysel,tsmear,fake_unit,wtk, 1, nband(1),&
794 510540 : extfpmd_nbdbuf=extfpmd_nbdbuf,rcpaw=rcpaw_getnel)
795 :
796 : ! Compute the number of free electrons of the extfpmd model
797 : ! with corresponding chemical potential and add to nelect bounds.
798 510540 : if(present(extfpmd)) then
799 495002 : if(associated(extfpmd)) then
800 4461 : call extfpmd%compute_nelect(fermie_mid,nband,nelectmid,nkpt,nspinor,nsppol,wtk)
801 : end if
802 : end if
803 :
804 : !write(std_out,'(a,i0,1x, 3(a,es13.5))' ) " iter: ", ii, &
805 : ! ' fermi_mid: ',fermimid * Ha_eV, ', n_mid: ',nelectmid, &
806 : ! ", (n_mid-nelect)/nelect: ", (nelectmid - nelect) / nelect
807 :
808 : !if (nelectmid > nelect * (one - tol)) then
809 : ! fermihi = fermimid
810 : ! nelecthi = nelectmid
811 : !end if
812 : !if (nelectmid < nelect * (one + tol)) then
813 : ! fermilo = fermimid
814 : ! nelectlo = nelectmid
815 : !end if
816 510540 : if(nelectmid>nelect*(one-tol14))then
817 254597 : fermie_hi=fermie_mid
818 254597 : nelecthi=nelectmid
819 : end if
820 510540 : if(nelectmid<nelect*(one+tol14))then
821 267926 : fermie_lo=fermie_mid
822 267926 : nelectlo=nelectmid
823 : end if
824 :
825 : else
826 :
827 : call getnel(doccde,dosdeltae,eigen,entropye,fermie_mid,fermie_mid,maxocc,mband,nband,&
828 480 : nelectmid,nkpt,nsppol,occ,occopt,option1,tphysel,tsmear,fake_unit,wtk, ivalence+1, nband(1))
829 : call getnel(doccde,dosdeltae,eigen,entropyh,fermih_mid,fermih_mid,maxocc,mband,nband,&
830 480 : nholesmid,nkpt,nsppol,occ,occopt,option1,tphysel,tsmear,fake_unit,wtk,1,ivalence)
831 :
832 480 : if(nelectmid>ne_qFD*(one-tol14))then
833 260 : fermie_hi = fermie_mid
834 260 : nelecthi = nelectmid
835 220 : else if (nelectmid<ne_qFD*(one-tol14))then
836 220 : fermie_lo = fermie_mid
837 220 : nelectlo = nelectmid
838 : end if
839 480 : if(nholesmid>(nelect-nh_qFD)*(one-tol14))then
840 266 : fermih_hi = fermih_mid
841 266 : nholeshi = nholesmid
842 214 : else if(nholesmid<(nelect-nh_qFD)*(one+tol14))then
843 214 : fermih_lo = fermih_mid
844 214 : nholeslo = nholesmid
845 : end if
846 :
847 : end if
848 :
849 : !if (abs(nelectmid - nelect) <= nelect*two*tol) exit
850 : !write(std_out,'(2(a,es13.5))' )' bisection move: fermi_lo: ',fermilo * Ha_eV,", fermi_hi: ", fermihi * Ha_eV
851 :
852 : ! if (abs(nelecthi - nelectlo) <= nelect*two*tol .or. &
853 : ! abs(fermihi - fermilo) <= tol * abs(fermihi + fermilo) ) exit
854 511020 : if (occopt /= 9) then
855 510540 : if( abs(nelecthi-nelectlo) <= nelect*two*tol14 .or. abs(fermie_hi-fermie_lo) <= tol14*abs(fermie_hi+fermie_lo) ) exit
856 : else
857 : if( ( abs(nelecthi-nelectlo) <= ne_qFD*two*tol14 .or. &
858 480 : & abs(fermie_hi-fermie_lo) <= tol14*abs(fermie_hi+fermie_lo) ) .and. &
859 : ( abs(nholeshi-nholeslo) <= (nelect-nh_qFD)*two*tol14 .or. &
860 : & abs(fermih_hi-fermih_lo) <= tol14*abs(fermih_hi+fermih_lo) ) ) exit
861 : end if
862 :
863 511020 : if (ii == niter_max) then
864 : write(msg,'(a,i0,3a,es22.14,a,es22.14,a)')&
865 0 : 'It was not possible to find Fermi energy in ',niter_max,' max bisections.',ch10,&
866 0 : 'nelecthi: ',nelecthi,', and nelectlo: ',nelectlo,'.'
867 0 : ABI_BUG(msg)
868 0 : if (occopt == 9) then
869 : write(msg,'(a,es22.14,a,es22.14,a)')&
870 0 : 'nholesi = ',nholeshi,', and holeslo = ',nholeslo,'.'
871 : end if
872 : end if
873 : end do ! End of bisection loop
874 :
875 14056 : fermie = fermie_mid
876 14056 : entropy= entropye
877 :
878 14056 : if (occopt /= 9) then
879 : write(msg, '(2(a,f14.6),a,i0)' ) &
880 14046 : & ' newocc: new Fermi energy is ',fermie,' , with nelect=',nelectmid,', Number of bisection calls: ',ii
881 : else
882 10 : fermih=fermih_mid
883 10 : entropy = entropy + entropyh ! CP: adding entropy of the holes subsystem
884 : write(msg, '(2(a,f14.6),a,i0)' ) &
885 10 : & ' newocc: new Fermi energy for excited electrons is ',fermie,' , with ne_qFD=',nelectmid,', Number of bisection calls: ',ii
886 10 : call wrtout(std_out,msg,'COLL')
887 : write(msg, '(2(a,f14.6),a,i0)' ) &
888 10 : & ' newocc: new Fermi energy for excited holes is ',fermih,' , with nh_qFD=',nelect-nholesmid,&
889 20 : & ', Number of bisection calls: ',ii
890 : end if
891 14056 : call wrtout(std_out,msg)
892 :
893 : ! Compute occupation numbers for prtstm/=0, close to the Fermi energy
894 14056 : if (present(stmbias)) then
895 :
896 13575 : if (abs(stmbias) > tol10) then
897 :
898 : ! Prevent use with occopt = 9 so far
899 : ! XG220804: This test is not needed, as prtstm/=0 must be used with occopt==7, as tested in chkinp.F90
900 3 : if (occopt == 9) then
901 0 : ABI_ERROR('Occopt 9 and prtstm /=0 not implemented together. Change occopt or prtstm.')
902 : end if
903 :
904 3 : fermie_biased = fermie - stmbias
905 9 : ABI_MALLOC(occt,(mband*nkpt*nsppol))
906 :
907 : call getnel(doccde,dosdeltae,eigen,entropy,fermie_biased,fermie_biased,maxocc,mband,nband,&
908 : & nelect_biased,nkpt,nsppol,occt,occopt,option1,tphysel,tsmear,fake_unit,wtk,1,nband(1),&
909 3 : & extfpmd_nbdbuf=extfpmd_nbdbuf)
910 273 : occ(:)=occ(:)-occt(:)
911 :
912 : ! Possibly filter a specific band contribution
913 3 : if (present(prtstm)) then
914 3 : if (prtstm < 0)then
915 : ibantot=1
916 2 : do isppol=1,nsppol
917 8 : do ikpt=1,nkpt
918 6 : nban=nband(ikpt+(isppol-1)*nkpt)
919 97 : do iban=1,nban
920 90 : if(iban/=abs(prtstm)) occ(ibantot)=zero
921 96 : ibantot=ibantot+1
922 : end do ! iban
923 : end do ! ikpt
924 : end do ! isppol
925 : end if ! prtstm < 0
926 : end if ! present(prtstm)
927 :
928 3 : nelect_biased = abs(nelectmid - nelect_biased)
929 : ! Here, arrange to have globally positive occupation numbers, irrespective of the stmbias sign
930 93 : if (-stmbias > tol10) occ(:) = -occ(:)
931 3 : ABI_FREE(occt)
932 :
933 3 : write(msg,'(a,f14.6)')' newocc: the number of electrons in the STM range is nelect_biased=',nelect_biased
934 3 : call wrtout(std_out,msg)
935 : end if
936 : endif ! present(stmbias)
937 :
938 : else
939 : ! Calculations with a specified moment
940 : ! Bisection loop
941 132 : cnt2=0
942 132 : cnt3=0
943 132 : entropy=zero
944 132 : maxocc=one
945 396 : ABI_MALLOC(doccdet,(nkpt*mband))
946 264 : ABI_MALLOC(eigent,(nkpt*mband))
947 264 : ABI_MALLOC(occt,(nkpt*mband))
948 396 : ABI_MALLOC(nbandt,(nkpt))
949 :
950 396 : do is = 1, nsppol
951 264 : nelect_tmp = nelectt(is)
952 264 : fermie_hi = fermie_hit(is) ! CP modify name
953 264 : fermie_lo = fermie_lot(is) ! CP modify name
954 264 : nelecthi = nelecthit(is)
955 264 : nelectlo = nelectlot(is)
956 : ! write(std_out,'(a,i1,3(f8.4,1x))') "Spin, N(spin):", is, nelect, fermihi, fermilo
957 : ! write(std_out,'(a,2(f8.4,1x))') "Hi, lo:", nelecthi, nelectlo
958 :
959 9119 : do ii=1,niter_max
960 9119 : fermie_mid_tmp=(fermie_hi+fermie_lo)/2.0_dp ! CP modify name
961 : ! temporary arrays
962 9119 : cnt = 0
963 25965 : do ik = 1, nkpt
964 16846 : nbandt(ik) = mband
965 222689 : do ib = 1, mband
966 196724 : cnt = cnt + 1
967 196724 : eigent(cnt) = eigen(cnt+cnt2)
968 196724 : occt(cnt) = occ(cnt+cnt2)
969 213570 : doccdet(cnt) = doccde(cnt+cnt2)
970 : end do
971 : end do
972 :
973 : ! Produce nelectmid from fermimid
974 : call getnel(doccdet,dosdeltae,eigent,entropy_tmp,fermie_mid_tmp,fermie_mid_tmp,maxocc,mband,nbandt,&
975 : nelectmid,nkpt,1,occt,occopt,option1,tphysel,tsmear,fake_unit,wtk,1,nband(1),&
976 9119 : extfpmd_nbdbuf=extfpmd_nbdbuf)
977 :
978 9119 : entropyet(is) = entropy_tmp
979 : fermie_midt(is) = fermie_mid_tmp
980 9119 : fermie_mid = fermie_midt(is)
981 :
982 : ! temporary arrays
983 9119 : cnt = 0
984 25965 : do ik = 1, nkpt
985 222689 : do ib = 1, mband
986 196724 : cnt = cnt + 1
987 196724 : occ(cnt+cnt2) = occt(cnt)
988 213570 : doccde(cnt+cnt2) = doccdet(cnt)
989 : end do
990 : end do
991 : ! write(std_out,'(a,es24.16,a,es24.16)' )' newocc: from fermi=',fermimid,', getnel gives nelect=',nelectmid
992 :
993 9119 : if(nelectmid>=nelect_tmp)then
994 4474 : fermie_hi=fermie_mid_tmp
995 4474 : nelecthi=nelectmid
996 : else
997 4645 : fermie_lo=fermie_mid_tmp
998 4645 : nelectlo=nelectmid
999 : end if
1000 9119 : if( abs(nelecthi-nelectlo) <= 1.0d-13 .or. abs(fermie_hi-fermie_lo) <= 0.5d-14*abs(fermie_hi+fermie_lo) ) exit
1001 :
1002 9119 : if(ii==niter_max)then
1003 : write(msg,'(a,i3,3a,es22.14,a,es22.14,a)')&
1004 0 : 'It was not possible to find Fermi energy in ',niter_max,' bisections.',ch10,&
1005 0 : 'nelecthi: ',nelecthi,', and nelectlo: ',nelectlo,'.'
1006 0 : ABI_BUG(msg)
1007 : end if
1008 : end do ! End of bisection loop
1009 :
1010 264 : cnt2 = cnt2 + nkpt*mband
1011 264 : entropy = entropy + entropyet(is)
1012 264 : fermie=fermie_mid
1013 : write(msg, '(a,i2,a,f14.6,a,f14.6,a,a,i4)' ) &
1014 264 : ' newocc: new Fermi energy for spin ', is, ' is ',fermie,' , with nelect: ',nelectmid,ch10,&
1015 528 : ' Number of bisection calls =',ii
1016 396 : call wrtout(std_out,msg)
1017 :
1018 : end do ! spin
1019 :
1020 132 : ABI_FREE(doccdet)
1021 132 : ABI_FREE(eigent)
1022 132 : ABI_FREE(nbandt)
1023 132 : ABI_FREE(occt)
1024 :
1025 : end if ! End of logical on fixed moment calculations
1026 :
1027 14188 : if(present(rcpaw)) then
1028 13651 : if(associated(rcpaw)) then
1029 42 : do itypat=1,rcpaw%ntypat
1030 42 : if(rcpaw%atm(itypat)%zcore_orig>zero) then
1031 126 : rcpaw%atm(itypat)%occ_res=rcpaw%atm(itypat)%occ_res+rcpaw%atm(itypat)%occ
1032 : endif
1033 : enddo
1034 : endif
1035 : endif
1036 :
1037 : !write(std_out,*) "kT*Entropy:", entropy*tsmear
1038 :
1039 : ! MG: If you are wondering why this part is now disabled by default consider that this output
1040 : ! is produced many times in the SCF cycle and in EPH we have to call this routine for
1041 : ! several temperature and the log becomes unreadable.
1042 : ! If you really need to look at the occupation factors use prtvol > 0.
1043 14188 : nkpt_eff = nkpt
1044 14188 : if (prtvol == 0) nkpt_eff = 0
1045 14188 : if (prtvol == 1) nkpt_eff = min(nkpt_max, nkpt)
1046 :
1047 14188 : if (nsppol == 1)then
1048 :
1049 12377 : if (nkpt_eff /= 0) then
1050 2489 : write(msg, '(a,i0,a)' )' newocc: computed new occ. numbers for occopt= ',occopt,' , spin-unpolarized case. '
1051 2489 : call wrtout(std_out,msg)
1052 20346 : do ikpt=1,nkpt_eff
1053 17857 : write(msg,'(a,i4,a)' ) ' k-point number ',ikpt,' :'
1054 45481 : do ii=0,(nband(1)-1)/12
1055 25299 : if (ii == 3 .and. prtvol /= 0) exit
1056 25135 : write(msg,'(12f6.3)') occ(1+ii*12+(ikpt-1)*nband(1):min(12+ii*12,nband(1))+(ikpt-1)*nband(1))
1057 42992 : call wrtout(std_out,msg)
1058 : end do
1059 : end do
1060 2489 : if (nkpt /= nkpt_eff) call wrtout(std_out,' newocc: prtvol=0, stop printing more k-point information')
1061 :
1062 : !call wrtout(std_out,' newocc: corresponding derivatives are ')
1063 : !do ikpt=1,nkpt_eff
1064 : !write(msg,'(a,i4,a)' ) ' k-point number ',ikpt,' :'
1065 : !do ii=0,(nband(1)-1)/12
1066 : !write(msg,'(12f6.1)') doccde(1+ii*12+(ikpt-1)*nband(1):min(12+ii*12,nband(1))+(ikpt-1)*nband(1))
1067 : !call wrtout(std_out,msg)
1068 : !end do
1069 : !end do
1070 : !if(nkpt/=nkpt_eff)then
1071 : ! call wrtout(std_out,'newocc: prtvol=0, stop printing more k-point information')
1072 : !end if
1073 : end if
1074 :
1075 : else
1076 :
1077 1811 : if (nkpt_eff /= 0) then
1078 721 : write(msg, '(a,i0,2a)' )' newocc: computed new occupation numbers for occopt= ',occopt,ch10,' (1) spin up values '
1079 721 : call wrtout(std_out, msg)
1080 6291 : do ikpt=1,nkpt_eff
1081 5570 : write(msg,'(a,i0,a)' ) ' k-point number ',ikpt,':'
1082 16282 : do ii=0,(nband(1)-1)/12
1083 9991 : if (ii == 3 .and. prtvol /= 0) exit
1084 9991 : write(msg,'(12f6.3)') occ(1+ii*12+(ikpt-1)*nband(1):min(12+ii*12,nband(1))+(ikpt-1)*nband(1))
1085 15561 : call wrtout(std_out,msg)
1086 : end do
1087 : end do
1088 721 : if (nkpt/=nkpt_eff) call wrtout(std_out,' newocc: prtvol=0, stop printing more k-point information')
1089 :
1090 721 : call wrtout(std_out,' (2) spin down values ')
1091 6291 : do ikpt=1,nkpt_eff
1092 16282 : do ii=0,(nband(1)-1)/12
1093 9991 : if (ii == 3 .and. prtvol /= 0) exit
1094 9991 : write(msg,'(12f6.3)') occ( 1+ii*12+(ikpt-1+nkpt)*nband(1):min(12+ii*12,nband(1))+(ikpt-1+nkpt)*nband(1) )
1095 15561 : call wrtout(std_out,msg)
1096 : end do
1097 : end do
1098 721 : if(nkpt/=nkpt_eff) call wrtout(std_out,' newocc: prtvol=0, stop printing more k-point information')
1099 : end if
1100 :
1101 : end if ! End choice based on spin
1102 :
1103 14188 : rcpaw_getnel=>null()
1104 :
1105 14188 : call timab(74,2,tsec)
1106 :
1107 : DBG_EXIT("COLL")
1108 :
1109 14188 : end subroutine newocc
1110 : !!***
1111 :
1112 : !!****f* m_occ/init_occ_ent
1113 : !! NAME
1114 : !! init_occ_ent
1115 : !!
1116 : !! FUNCTION
1117 : !!
1118 : !! INPUTS
1119 : !!
1120 : !! OUTPUT
1121 : !!
1122 : !! SOURCE
1123 :
1124 549521 : subroutine init_occ_ent(entfun,limit,nptsdiv2,occfun,occopt,option,smdfun,tphysel,tsmear,tsmearinv,xgrid)
1125 :
1126 : !Arguments ------------------------------------
1127 : !scalars
1128 : integer,intent(in) :: occopt,option
1129 : real(dp),intent(in) :: tphysel,tsmear
1130 : integer,intent(inout) :: nptsdiv2
1131 : real(dp),intent(out) :: limit,tsmearinv
1132 : real(dp),intent(inout) :: entfun(-nptsdiv2:nptsdiv2,2),occfun(-nptsdiv2:nptsdiv2,2)
1133 : real(dp),intent(inout) :: smdfun(-nptsdiv2:nptsdiv2,2),xgrid(-nptsdiv2:nptsdiv2)
1134 :
1135 : !Local variables-------------------------------
1136 : !scalars
1137 : integer :: algo,ii,jj,nconvd2,nmaxFD,nminFD
1138 : integer,save :: dblsmr,occopt_prev=-9999
1139 : real(dp),save :: convlim,incconv,limit_occ,tphysel_prev=-9999,tsmear_prev=-9999
1140 : real(dp) :: aa,dsqrpi,encorr,factor
1141 : real(dp) :: expinc,expx22,expxo2,gauss,increm
1142 : real(dp) :: resFD1,resFD2,resFD3,resFD4,resmom,resmom1,resmom2
1143 : real(dp) :: resmom3,resmom4,secmom,smom1,smom2,thdmom,tmom1,tmom2,tmpexpsum
1144 : real(dp) :: tmpsmdfun,tratio,tt,xx,yp1,ypn
1145 : character(len=500) :: msg
1146 : !arrays
1147 : real(dp),save :: entfun_prev(-nptsdiv2_def:nptsdiv2_def,2),occfun_prev(-nptsdiv2_def:nptsdiv2_def,2)
1148 : real(dp),save :: smdfun_prev(-nptsdiv2_def:nptsdiv2_def,2),xgrid_prev(-nptsdiv2_def:nptsdiv2_def)
1149 549521 : real(dp),allocatable :: entder(:),occder(:),smd1(:),smd2(:)
1150 549521 : real(dp),allocatable :: smdder(:),tgrid(:),work(:),workfun(:)
1151 : ! *************************************************************************
1152 :
1153 : ! Initialize the occupation function and generalized entropy function,
1154 : ! at the beginning, or if occopt changed
1155 :
1156 549521 : if(option==-1)then
1157 0 : nptsdiv2 = nptsdiv2_def
1158 0 : return
1159 : end if
1160 :
1161 549521 : if (occopt_prev/=occopt .or. abs(tsmear_prev-tsmear) >tol12 .or. abs(tphysel_prev-tphysel)>tol12) then
1162 753 : occopt_prev=occopt
1163 753 : tsmear_prev=tsmear
1164 753 : tphysel_prev=tphysel
1165 :
1166 : ! Check whether input values of tphysel tsmear and occopt are consistent
1167 753 : dblsmr = 0
1168 753 : if (abs(tphysel)>tol12) then
1169 : ! Use re-smearing scheme
1170 10 : if (abs(tsmear)>tol12) then
1171 10 : dblsmr = 1
1172 : ! Use FD occupations (one smearing) only with "physical" temperature tphysel
1173 : ! CP modify
1174 : !else if (occopt /= 3) then
1175 : ! write(msg, '(a,i6,a)' )' tphysel /= 0, tsmear == 0, but occopt is not = 3, but ',occopt,'.'
1176 0 : else if (occopt /= 3 .and. occopt/=9) then
1177 0 : write(msg, '(a,i6,a)' )' tphysel /= 0, tsmear == 0, but occopt is not = 3 or 9, but ',occopt,'.'
1178 : ! End CP modify
1179 0 : ABI_ERROR(msg)
1180 : end if
1181 : end if
1182 :
1183 753 : ABI_MALLOC(entder,(-nptsdiv2_def:nptsdiv2_def))
1184 753 : ABI_MALLOC(occder,(-nptsdiv2_def:nptsdiv2_def))
1185 753 : ABI_MALLOC(smdder,(-nptsdiv2_def:nptsdiv2_def))
1186 753 : ABI_MALLOC(workfun,(-nptsdiv2_def:nptsdiv2_def))
1187 753 : ABI_MALLOC(work,(-nptsdiv2_def:nptsdiv2_def))
1188 :
1189 : ! Prepare the points on the grid
1190 : ! limit is the value of the argument that will give 0.0 or 1.0 , with
1191 : ! less than about 1.0d-15 error for 4<=occopt<=8, and less than about 1.0d-12
1192 : ! error for occopt==3. It is not worth to compute the function beyond
1193 : ! that point. Even with a less severe requirement, it is significantly
1194 : ! larger for occopt==3, with an exponential
1195 : ! tail, than for the other occupation functions, with a Gaussian tail.
1196 : ! Note that these values are useful in newocc.f also.
1197 753 : limit_occ=6.0_dp
1198 : ! CP modify
1199 : !if(occopt==3)limit_occ=30.0_dp
1200 753 : if(occopt==3 .or. occopt==9)limit_occ=30.0_dp
1201 : ! End CP modify
1202 753 : if(dblsmr /= 0) then
1203 10 : tratio = tsmear / tphysel
1204 10 : limit_occ=30.0_dp + 6.0_dp*tratio
1205 : end if
1206 :
1207 : ! With nptsdiv2_def=6000 (thus increm=0.001 for 4<=occopt<=8,
1208 : ! and increm=0.005 for occopt==3, the O(1/N4) algorithm gives 1.0d-12
1209 : ! accuracy on the stored values occfun and entfun. These, together
1210 : ! with smdfun and xgrid_prev, need permanently about 0.67 MB, which is affordable.
1211 753 : increm=limit_occ/nptsdiv2_def
1212 9037506 : do ii=-nptsdiv2_def,nptsdiv2_def
1213 9037506 : xgrid_prev(ii)=ii*increm
1214 : end do
1215 :
1216 : ! ---------------------------------------------------------
1217 : ! Ordinary (unique) smearing function
1218 : ! ---------------------------------------------------------
1219 753 : if (dblsmr == 0) then
1220 :
1221 : ! Compute the unnormalized smeared delta function between -limit_occ and +limit_occ
1222 : ! (well, they are actually normalized ...)
1223 :
1224 : ! CP modify
1225 : !if(occopt==3)then
1226 : if(occopt==3 .or. occopt==9)then
1227 : ! End CP modify
1228 : ! Fermi-Dirac
1229 2448816 : do ii=0,nptsdiv2_def
1230 2448408 : xx=xgrid_prev(ii)
1231 2448408 : smdfun_prev( ii,1)=0.25_dp/(cosh(xx/2.0_dp)**2)
1232 2448816 : smdfun_prev(-ii,1)=smdfun_prev(ii,1)
1233 : end do
1234 :
1235 : else if(occopt==4 .or. occopt==5)then
1236 : ! Cold smearing of Marzari, two values of the "a" parameter being possible
1237 : ! first value gives minimization of the bump
1238 : if(occopt==4)aa=-.5634
1239 : ! second value gives monotonic occupation function
1240 122 : if(occopt==5)aa=-.8165
1241 :
1242 122 : dsqrpi=1.0_dp/sqrt(pi)
1243 732244 : do ii=0,nptsdiv2_def
1244 732122 : xx=xgrid_prev(ii)
1245 732122 : gauss=dsqrpi*exp(-xx**2)
1246 732122 : smdfun_prev( ii,1)=gauss*(1.5_dp+xx*(-aa*1.5_dp+xx*(-1.0_dp+aa*xx)))
1247 732244 : smdfun_prev(-ii,1)=gauss*(1.5_dp+xx*( aa*1.5_dp+xx*(-1.0_dp-aa*xx)))
1248 : end do
1249 :
1250 : else if(occopt==6)then
1251 :
1252 : ! First order Hermite-Gaussian of Paxton and Methfessel
1253 30010 : dsqrpi=1.0_dp/sqrt(pi)
1254 30010 : do ii=0,nptsdiv2_def
1255 30005 : xx=xgrid_prev(ii)
1256 30005 : smdfun_prev( ii,1)=dsqrpi*(1.5_dp-xx**2)*exp(-xx**2)
1257 30010 : smdfun_prev(-ii,1)=smdfun_prev(ii,1)
1258 : end do
1259 :
1260 : else if(occopt==7)then
1261 :
1262 : ! Gaussian smearing
1263 1248416 : dsqrpi=1.0_dp/sqrt(pi)
1264 1248416 : do ii=0,nptsdiv2_def
1265 1248208 : xx=xgrid_prev(ii)
1266 1248208 : smdfun_prev( ii,1)=dsqrpi*exp(-xx**2)
1267 1248416 : smdfun_prev(-ii,1)=smdfun_prev(ii,1)
1268 : end do
1269 :
1270 : else if(occopt==8)then
1271 :
1272 : ! Constant value of the delta function over the smearing interval, for testing purposes only.
1273 0 : do ii=0,nptsdiv2_def
1274 0 : xx=xgrid_prev(ii)
1275 0 : if(xx>half+tol8)then
1276 0 : smdfun_prev( ii,1)=zero
1277 0 : else if(xx<half-tol8)then
1278 0 : smdfun_prev( ii,1)=one
1279 : else
1280 0 : smdfun_prev( ii,1)=half
1281 : end if
1282 0 : smdfun_prev(-ii,1)=smdfun_prev(ii,1)
1283 : end do
1284 :
1285 : else
1286 0 : ABI_BUG(sjoin('Occopt: ', itoa(occopt),' is not allowed in getnel.'))
1287 : end if
1288 :
1289 : else if (dblsmr /= 0) then
1290 : ! ---------------------------------------------------------
1291 : ! smear FD delta with occopt delta calculated in smdfun_prev
1292 : ! ---------------------------------------------------------
1293 :
1294 10 : nconvd2 = 6000
1295 10 : convlim = 10.0_dp
1296 10 : incconv = convlim / nconvd2
1297 :
1298 : ! store smearing functions in smd1 and smd2
1299 10 : ABI_MALLOC(smd1,(-nconvd2:nconvd2))
1300 10 : ABI_MALLOC(smd2,(-nconvd2:nconvd2))
1301 10 : ABI_MALLOC(tgrid,(-nconvd2:nconvd2))
1302 :
1303 : ! FD function in smd1( ii) and second smearing delta in smd2( ii)
1304 : !
1305 : ! smd1(:) contains delta_FD ( x )
1306 60020 : do ii=0,nconvd2
1307 60010 : tgrid(ii)=ii*incconv
1308 60010 : tgrid(-ii)=-tgrid(ii)
1309 60010 : tt=tgrid(ii)
1310 60010 : smd1( ii)=0.25_dp/(cosh(tt/2.0_dp)**2)
1311 60020 : smd1(-ii)=smd1(ii)
1312 : end do
1313 :
1314 : ! check input values of occopt and fill smd2(:) with appropriate data:
1315 : ! smd2(:) contains delta_resmear ( x )
1316 : ! CP modify
1317 : !if(occopt == 3) then
1318 : if(occopt == 3 .or. occopt==9) then
1319 : ! End CP modify
1320 : write(msg, '(a,a)' )&
1321 0 : 'Occopt=3 is not allowed as a re-smearing.', &
1322 0 : 'Use a single FD, or re-smear with a different delta type (faster cutoff). '
1323 0 : ABI_ERROR(msg)
1324 : else if(occopt==4 .or. occopt==5)then
1325 : ! Cold smearing of Marzari, two values of the "a" parameter being possible
1326 : ! first value gives minimization of the bump
1327 : if(occopt==4)aa=-.5634
1328 : ! second value gives monotonic occupation function
1329 10 : if(occopt==5)aa=-.8165
1330 :
1331 10 : dsqrpi=1.0_dp/sqrt(pi)
1332 60020 : do ii=0,nconvd2
1333 60010 : tt=tgrid(ii)
1334 60010 : gauss=dsqrpi*exp(-tt**2)
1335 60010 : smd2( ii)=gauss*(1.5_dp+tt*(-aa*1.5_dp+tt*(-1.0_dp+aa*tt)))
1336 60020 : smd2(-ii)=gauss*(1.5_dp+tt*( aa*1.5_dp+tt*(-1.0_dp-aa*tt)))
1337 : end do
1338 : else if(occopt==6)then
1339 0 : dsqrpi=1.0_dp/sqrt(pi)
1340 0 : do ii=0,nconvd2
1341 0 : tt=tgrid(ii)
1342 0 : smd2( ii)=dsqrpi*(1.5_dp-tt**2)*exp(-tt**2)
1343 0 : smd2(-ii)=smd2(ii)
1344 : end do
1345 : else if(occopt==7)then
1346 0 : dsqrpi=1.0_dp/sqrt(pi)
1347 0 : do ii=0,nconvd2
1348 0 : tt=tgrid(ii)
1349 0 : smd2( ii)=dsqrpi*exp(-tt**2)
1350 0 : smd2(-ii)=smd2(ii)
1351 : end do
1352 : else if(occopt==8)then
1353 0 : do ii=0,nconvd2
1354 0 : tt=tgrid(ii)
1355 0 : if(tt>half+tol8)then
1356 0 : smd2( ii)=zero
1357 0 : else if(tt<half-tol8)then
1358 0 : smd2( ii)=one
1359 : else
1360 0 : smd2( ii)=half
1361 : end if
1362 0 : smd2(-ii)=smd2(ii)
1363 : end do
1364 : else
1365 0 : ABI_BUG(sjoin('Occopt: ', itoa(occopt),' is not allowed in getnel.'))
1366 : end if
1367 :
1368 : ! Use O(1/N4) algorithm from Num Rec (see below)
1369 : !
1370 : ! The grid for the convoluted delta is taken (conservatively)
1371 : ! to be that for the FD delta ie 6000 pts in [-limit_occ;limit_occ]
1372 : ! Smearing functions are given on [-dbllim;dbllim] and the grid must
1373 : ! superpose the normal grid on [-limit_occ:limit_occ]
1374 : ! The maximal interval for integration of the convolution is
1375 : ! [-dbllim+limit_occ+lim(delta2);dbllim-limit_occ-lim(delta2)] =
1376 : ! [-dbllim+36;dbllim-36]
1377 :
1378 : ! test the smdFD function for extreme values:
1379 : ! do jj=-nptsdiv2_def,-nptsdiv2_def
1380 : ! do ii=-nconvd2+4,nconvd2
1381 : ! call smdFD(xgrid_prev(jj) - tgrid(ii)*tratio, resFD)
1382 : ! write(std_out,*) 'ii jj = ', ii,jj, ' smdFD (', xgrid_prev(jj) - tgrid(ii)*tratio, ') ', resFD
1383 : ! end do
1384 : ! end do
1385 :
1386 10 : expinc = exp(half*incconv*tratio)
1387 :
1388 : ! jj = position of point at which we are calculating smdfun_prev
1389 120020 : do jj=-nptsdiv2_def,nptsdiv2_def
1390 : ! Do not care about the 8 boundary points,
1391 : ! where the values should be extremely small anyway
1392 120010 : smdfun_prev(jj,1)=0.0_dp
1393 : ! only add contribution with delta_FD > 1.0d-100
1394 120010 : nmaxFD = floor (( maxFDarg+xgrid_prev(jj)) / tratio / incconv )
1395 120010 : nmaxFD = min (nmaxFD, nconvd2)
1396 120010 : nminFD = ceiling((-maxFDarg+xgrid_prev(jj)) / tratio / incconv )
1397 120010 : nminFD = max (nminFD, -nconvd2)
1398 :
1399 : ! Calculate the Fermi-Dirac distrib at point xgrid_prev(jj)-tgrid(ii)*tratio
1400 120010 : expxo2 = exp (-half*(xgrid_prev(jj) - (nminFD)*incconv*tratio))
1401 120010 : expx22 = expxo2*expxo2
1402 120010 : tmpexpsum = expxo2 / (expx22 + 1.0_dp)
1403 120010 : resFD4 = tmpexpsum * tmpexpsum
1404 120010 : expxo2 = expxo2*expinc
1405 120010 : expx22 = expxo2*expxo2
1406 120010 : tmpexpsum = expxo2 / (expx22 + 1.0_dp)
1407 120010 : resFD3 = tmpexpsum * tmpexpsum
1408 120010 : expxo2 = expxo2*expinc
1409 120010 : expx22 = expxo2*expxo2
1410 120010 : tmpexpsum = expxo2 / (expx22 + 1.0_dp)
1411 120010 : resFD2 = tmpexpsum * tmpexpsum
1412 120010 : expxo2 = expxo2*expinc
1413 120010 : expx22 = expxo2*expxo2
1414 120010 : tmpexpsum = expxo2 / (expx22 + 1.0_dp)
1415 120010 : resFD1 = tmpexpsum * tmpexpsum
1416 :
1417 : ! core contribution to the integral with constant weight (48)
1418 120010 : tmpsmdfun = 0.0_dp
1419 1339688356 : do ii=nminFD+4,nmaxFD-4
1420 1339568346 : expxo2 = expxo2*expinc
1421 : ! tmpexpsum = 1.0_dp / (expxo2 + 1.0_dp / expxo2 )
1422 1339568346 : expx22 = expxo2*expxo2
1423 1339568346 : tmpexpsum = expxo2 / (expx22 + 1.0_dp)
1424 1339688356 : tmpsmdfun = tmpsmdfun + smd2(ii) * tmpexpsum * tmpexpsum
1425 : end do
1426 :
1427 : ! Add on end contributions for show (both functions smd and smdFD are very small
1428 : smdfun_prev(jj,1)=smdfun_prev(jj,1) +48.0_dp*tmpsmdfun &
1429 : + 31.0_dp*smd2(nminFD+3)*resFD1 -11.0_dp*smd2(nminFD+2)*resFD2 &
1430 120010 : + 5.0_dp*smd2(nminFD+1)*resFD3 - smd2(nminFD)*resFD4
1431 :
1432 120010 : expxo2 = expxo2*expinc
1433 120010 : expx22 = expxo2*expxo2
1434 120010 : tmpexpsum = expxo2 / (expx22 + 1.0_dp)
1435 120010 : resFD1 = tmpexpsum * tmpexpsum
1436 120010 : expxo2 = expxo2*expinc
1437 120010 : expx22 = expxo2*expxo2
1438 120010 : tmpexpsum = expxo2 / (expx22 + 1.0_dp)
1439 120010 : resFD2 = tmpexpsum * tmpexpsum
1440 120010 : expxo2 = expxo2*expinc
1441 120010 : expx22 = expxo2*expxo2
1442 120010 : tmpexpsum = expxo2 / (expx22 + 1.0_dp)
1443 120010 : resFD3 = tmpexpsum * tmpexpsum
1444 120010 : expxo2 = expxo2*expinc
1445 120010 : expx22 = expxo2*expxo2
1446 120010 : tmpexpsum = expxo2 / (expx22 + 1.0_dp)
1447 120010 : resFD4 = tmpexpsum * tmpexpsum
1448 :
1449 : ! Contribution above
1450 : smdfun_prev(jj,1)=smdfun_prev(jj,1) &
1451 : + 31.0_dp*smd2(nmaxFD-3)*resFD1 -11.0_dp*smd2(nmaxFD-2)*resFD2 &
1452 120010 : + 5.0_dp*smd2(nmaxFD-1)*resFD3 - smd2(nmaxFD)*resFD4
1453 120020 : smdfun_prev(jj,1)=incconv*smdfun_prev(jj,1)/48.0_dp
1454 : end do
1455 :
1456 10 : secmom = 0.0_dp
1457 10 : thdmom = 0.0_dp
1458 10 : resmom4 = xgrid_prev(-nptsdiv2_def )*xgrid_prev(-nptsdiv2_def )*smdfun_prev(-nptsdiv2_def , 1)
1459 10 : resmom3 = xgrid_prev(-nptsdiv2_def+1)*xgrid_prev(-nptsdiv2_def+1)*smdfun_prev(-nptsdiv2_def+1, 1)
1460 10 : resmom2 = xgrid_prev(-nptsdiv2_def+2)*xgrid_prev(-nptsdiv2_def+2)*smdfun_prev(-nptsdiv2_def+2, 1)
1461 10 : resmom1 = xgrid_prev(-nptsdiv2_def+3)*xgrid_prev(-nptsdiv2_def+3)*smdfun_prev(-nptsdiv2_def+3, 1)
1462 10 : resmom = xgrid_prev(-nptsdiv2_def+4)*xgrid_prev(-nptsdiv2_def+4)*smdfun_prev(-nptsdiv2_def+4, 1)
1463 10 : do ii=-nptsdiv2_def+4,nptsdiv2_def-1
1464 : secmom = secmom + &
1465 : & ( 17.0_dp*xgrid_prev(ii) *xgrid_prev(ii) *smdfun_prev(ii, 1) &
1466 : & +42.0_dp*xgrid_prev(ii-1)*xgrid_prev(ii-1)*smdfun_prev(ii-1,1) &
1467 : & -16.0_dp*xgrid_prev(ii-2)*xgrid_prev(ii-2)*smdfun_prev(ii-2,1) &
1468 : & + 6.0_dp*xgrid_prev(ii-3)*xgrid_prev(ii-3)*smdfun_prev(ii-3,1) &
1469 : & - xgrid_prev(ii-4)*xgrid_prev(ii-4)*smdfun_prev(ii-4,1) )
1470 : resmom4 = resmom3
1471 : resmom3 = resmom2
1472 : resmom2 = resmom1
1473 : resmom1 = resmom
1474 : resmom = xgrid_prev(ii+1) *xgrid_prev(ii+1) *smdfun_prev(ii+1, 1)
1475 : end do
1476 10 : secmom=increm * secmom / 48.0_dp
1477 : ! thdmom=increm * thdmom / 48.0_dp
1478 : !
1479 : ! smom1 = second moment of delta in smd1(:)
1480 : ! smom2 = second moment of delta in smd2(:)
1481 : !
1482 : smom1 = 0.0_dp
1483 : smom2 = 0.0_dp
1484 : tmom1 = 0.0_dp
1485 : tmom2 = 0.0_dp
1486 : do ii=-nconvd2+4,nconvd2
1487 : smom1 = smom1+ &
1488 : & ( 17.0_dp*tgrid(ii) *tgrid(ii) *smd1(ii) &
1489 : & +42.0_dp*tgrid(ii-1)*tgrid(ii-1)*smd1(ii-1) &
1490 : & -16.0_dp*tgrid(ii-2)*tgrid(ii-2)*smd1(ii-2) &
1491 : & + 6.0_dp*tgrid(ii-3)*tgrid(ii-3)*smd1(ii-3) &
1492 : & - tgrid(ii-4)*tgrid(ii-4)*smd1(ii-4) )
1493 : smom2 = smom2+ &
1494 : & ( 17.0_dp*tgrid(ii) *tgrid(ii) *smd2(ii ) &
1495 : & +42.0_dp*tgrid(ii-1)*tgrid(ii-1)*smd2(ii-1) &
1496 : & -16.0_dp*tgrid(ii-2)*tgrid(ii-2)*smd2(ii-2) &
1497 : & + 6.0_dp*tgrid(ii-3)*tgrid(ii-3)*smd2(ii-3) &
1498 : & - tgrid(ii-4)*tgrid(ii-4)*smd2(ii-4) )
1499 : end do
1500 10 : smom1 =incconv * smom1 / 48.0_dp
1501 10 : smom2 =incconv * smom2 / 48.0_dp
1502 : ! tmom1 =incconv * tmom1 / 48.0_dp
1503 : ! tmom2 =incconv * tmom2 / 48.0_dp
1504 :
1505 10 : encorr = smom2*tratio*tratio/secmom
1506 :
1507 10 : ABI_FREE(tgrid)
1508 10 : ABI_FREE(smd1)
1509 10 : ABI_FREE(smd2)
1510 :
1511 : end if
1512 :
1513 : ! --------------------------------------------------------
1514 : ! end of smearing function initialisation, dblsmr case
1515 : ! --------------------------------------------------------
1516 :
1517 :
1518 : ! Now that the smeared delta function has been initialized, compute the
1519 : ! occupation function
1520 753 : occfun_prev(-nptsdiv2_def,1)=zero
1521 753 : entfun_prev(-nptsdiv2_def,1)=zero
1522 :
1523 : ! Different algorithms are possible, corresponding to the formulas
1524 : ! (4.1.11), (4.1.12) and (4.1.14) in Numerical recipes (pp 107 and 108),
1525 : ! with respective O(1/N2), O(1/N3), O(1/N4) convergence, where N is the
1526 : ! number of points in the interval.
1527 753 : algo=4
1528 :
1529 : if(algo==2)then
1530 :
1531 : ! Extended trapezoidal rule (4.1.11), taken in a cumulative way
1532 : do ii=-nptsdiv2_def+1,nptsdiv2_def
1533 : occfun_prev(ii,1)=occfun_prev(ii-1,1)+increm*(smdfun_prev(ii,1)+smdfun_prev(ii-1,1))/2.0_dp
1534 : entfun_prev(ii,1)=entfun_prev(ii-1,1)+increm*&
1535 : & ( -xgrid_prev(ii)*smdfun_prev(ii,1) -xgrid_prev(ii-1)*smdfun_prev(ii-1,1) )/2.0_dp
1536 : end do
1537 :
1538 : else if(algo==3)then
1539 :
1540 : ! Derived from (4.1.12). Converges as O(1/N3).
1541 : ! Do not care about the following points,
1542 : ! where the values are extremely small anyway
1543 : occfun_prev(-nptsdiv2_def+1,1)=0.0_dp ; entfun_prev(-nptsdiv2_def+1,1)=0.0_dp
1544 : do ii=-nptsdiv2_def+2,nptsdiv2_def
1545 : occfun_prev(ii,1)=occfun_prev(ii-1,1)+increm*&
1546 : & ( 5.0_dp*smdfun_prev(ii,1) + 8.0_dp*smdfun_prev(ii-1,1) - smdfun_prev(ii-2,1) )/12.0_dp
1547 : entfun_prev(ii,1)=entfun_prev(ii-1,1)+increm*&
1548 : & ( 5.0_dp*(-xgrid_prev(ii) )*smdfun_prev(ii,1) &
1549 : & +8.0_dp*(-xgrid_prev(ii-1))*smdfun_prev(ii-1,1)&
1550 : & - (-xgrid_prev(ii-2))*smdfun_prev(ii-2,1) )/12.0_dp
1551 : end do
1552 :
1553 : else if(algo==4)then
1554 :
1555 : ! Derived from (4.1.14)- alternative extended Simpsons rule. Converges as O(1/N4).
1556 : ! Do not care about the following points,
1557 : ! where the values are extremely small anyway
1558 753 : occfun_prev(-nptsdiv2_def+1,1)=0.0_dp ; entfun_prev(-nptsdiv2_def+1,1)=0.0_dp
1559 753 : occfun_prev(-nptsdiv2_def+2,1)=0.0_dp ; entfun_prev(-nptsdiv2_def+2,1)=0.0_dp
1560 753 : occfun_prev(-nptsdiv2_def+3,1)=0.0_dp ; entfun_prev(-nptsdiv2_def+3,1)=0.0_dp
1561 9034494 : do ii=-nptsdiv2_def+4,nptsdiv2_def
1562 : occfun_prev(ii,1)=occfun_prev(ii-1,1)+increm*&
1563 : & ( 17.0_dp*smdfun_prev(ii,1) &
1564 : & +42.0_dp*smdfun_prev(ii-1,1)&
1565 : & -16.0_dp*smdfun_prev(ii-2,1)&
1566 : & + 6.0_dp*smdfun_prev(ii-3,1)&
1567 9033741 : & - smdfun_prev(ii-4,1) )/48.0_dp
1568 : entfun_prev(ii,1)=entfun_prev(ii-1,1)+increm*&
1569 : & ( 17.0_dp*(-xgrid_prev(ii) )*smdfun_prev(ii,1) &
1570 : & +42.0_dp*(-xgrid_prev(ii-1))*smdfun_prev(ii-1,1)&
1571 : & -16.0_dp*(-xgrid_prev(ii-2))*smdfun_prev(ii-2,1)&
1572 : & + 6.0_dp*(-xgrid_prev(ii-3))*smdfun_prev(ii-3,1)&
1573 9034494 : & - (-xgrid_prev(ii-4))*smdfun_prev(ii-4,1) )/48.0_dp
1574 : end do
1575 :
1576 : end if ! End of choice between different algorithms for integration
1577 :
1578 : ! Normalize the functions (actually not needed for occopt=3..7)
1579 753 : factor=1.0_dp/occfun_prev(nptsdiv2_def,1)
1580 9037506 : smdfun_prev(:,1)=smdfun_prev(:,1)*factor
1581 9037506 : occfun_prev(:,1)=occfun_prev(:,1)*factor
1582 9037506 : entfun_prev(:,1)=entfun_prev(:,1)*factor
1583 :
1584 : ! Compute the cubic spline fitting of the smeared delta function
1585 753 : yp1=0.0_dp ; ypn=0.0_dp
1586 9037506 : workfun(:)=smdfun_prev(:,1)
1587 753 : call spline(xgrid_prev, workfun, (2*nptsdiv2_def+1), yp1, ypn, smdder)
1588 9037506 : smdfun_prev(:,2)=smdder(:)
1589 :
1590 : ! Compute the cubic spline fitting of the occupation function
1591 : yp1=0.0_dp ; ypn=0.0_dp
1592 9037506 : workfun(:)=occfun_prev(:,1)
1593 753 : call spline(xgrid_prev, workfun, (2*nptsdiv2_def+1), yp1, ypn, occder)
1594 9037506 : occfun_prev(:,2)=occder(:)
1595 :
1596 : ! Compute the cubic spline fitting of the entropy function
1597 : yp1=0.0_dp ; ypn=0.0_dp
1598 9037506 : workfun(:)=entfun_prev(:,1)
1599 753 : call spline(xgrid_prev, workfun, (2*nptsdiv2_def+1), yp1, ypn, entder)
1600 9037506 : entfun_prev(:,2)=entder(:)
1601 :
1602 753 : ABI_FREE(entder)
1603 753 : ABI_FREE(occder)
1604 753 : ABI_FREE(smdder)
1605 753 : ABI_FREE(work)
1606 753 : ABI_FREE(workfun)
1607 :
1608 : end if
1609 :
1610 549521 : if (abs(tphysel)<tol12) then
1611 547224 : if (tsmear == zero) then
1612 325 : tsmearinv = huge_tsmearinv
1613 : else
1614 546899 : tsmearinv=one/tsmear
1615 : end if
1616 : else
1617 2297 : tsmearinv=one/tphysel
1618 : end if
1619 :
1620 13191251605 : entfun(:,:) = entfun_prev(:,:)
1621 13191251605 : occfun(:,:) = occfun_prev(:,:)
1622 13191251605 : smdfun(:,:) = smdfun_prev(:,:)
1623 6595351042 : xgrid(:) = xgrid_prev(:)
1624 549521 : limit = limit_occ
1625 549521 : nptsdiv2 = nptsdiv2_def
1626 :
1627 : end subroutine init_occ_ent
1628 : !!***
1629 :
1630 : !!****f* m_occ/occeig
1631 : !! NAME
1632 : !! occeig
1633 : !!
1634 : !! FUNCTION
1635 : !! For each pair of active bands (m,n), generates ratios
1636 : !! that depend on the difference between occupation numbers and eigenvalues.
1637 : !!
1638 : !! INPUTS
1639 : !! doccde_k(nband_k)=derivative of occ_k wrt the energy
1640 : !! doccde_kq(nband_k)=derivative of occ_kq wrt the energy
1641 : !! eig0_k(nband_k)=GS eigenvalues at k
1642 : !! eig0_kq(nband_k)=GS eigenvalues at k+q
1643 : !! nband_k=number of bands
1644 : !! occopt=option for occupancies
1645 : !! occ_k(nband_k)=occupation number for each band at k
1646 : !! occ_kq(nband_k)=occupation number for each band at k+q
1647 : !!
1648 : !! OUTPUT
1649 : !! rocceig(nband_k,nband_k)$= (occ_{k,q}(m)-occ_k(n))/(eig0_{k,q}(m)-eig0_k(n))$,
1650 : !! if this ratio has been attributed to the band n, 0.0_dp otherwise
1651 : !!
1652 : !! NOTES
1653 : !! Supposing the occupations numbers differ:
1654 : !! if $abs(occ_{k,q}(m)) < abs(occ_k(n))$
1655 : !! $rocceig(m,n)=(occ_{k,q}(m)-occ_k(n))/(eig0_{k,q}(m)-eig0_k(n)) $
1656 : !! if $abs(occ_{k,q}(m))>abs(occ_k(n))$
1657 : !! rocceig(m,n)=0.0_dp
1658 : !!
1659 : !! If the occupation numbers are close enough, then
1660 : !! if the eigenvalues are also close, take the derivative
1661 : !! $ rocceig(m,n)=\frac{1}{2}*docc/deig0 $
1662 : !! otherwise,
1663 : !! $ rocceig(m,n)=\frac{1}{2}*(occ_{k,q}(m)-occ_k(n))/(eig0_{k,q}(m)-eig0_k(n))$
1664 : !!
1665 : !! SOURCE
1666 :
1667 855197 : subroutine occeig(doccde_k,doccde_kq,eig0_k,eig0_kq,nband_k,occopt,occ_k,occ_kq,rocceig)
1668 :
1669 : !Arguments ------------------------------------
1670 : !scalars
1671 : integer,intent(in) :: nband_k,occopt
1672 : !arrays
1673 : real(dp),intent(in) :: doccde_k(nband_k),doccde_kq(nband_k),eig0_k(nband_k)
1674 : real(dp),intent(in) :: eig0_kq(nband_k),occ_k(nband_k),occ_kq(nband_k)
1675 : real(dp),intent(out) :: rocceig(nband_k,nband_k)
1676 :
1677 : !Local variables-------------------------------
1678 : !scalars
1679 : integer :: ibandk,ibandkq
1680 : real(dp) :: diffabsocc,diffeig,diffocc,ratio,sumabsocc
1681 : character(len=500) :: msg
1682 : ! *************************************************************************
1683 :
1684 : ! The parameter tol5 defines the treshhold for degeneracy, and the width of the step function
1685 :
1686 40873801 : rocceig(:,:) = zero
1687 :
1688 5388783 : do ibandk=1,nband_k
1689 40873801 : do ibandkq=1,nband_k
1690 35485018 : diffeig=eig0_kq(ibandkq)-eig0_k(ibandk)
1691 35485018 : diffocc=occ_kq(ibandkq)-occ_k(ibandk)
1692 :
1693 35485018 : if( abs(diffeig) > tol5 ) then
1694 32384630 : ratio=diffocc/diffeig
1695 : else
1696 3100388 : if(occopt<3)then
1697 : ! In a non-metallic case, if the eigenvalues are degenerate,
1698 : ! the occupation numbers must also be degenerate, in which
1699 : ! case there is no contribution from this pair of bands
1700 2537164 : if( abs(diffocc) > tol5 ) then
1701 : write(msg,'(a,a,a,a,a,a,a,2(a,i4,a,es16.6,a,es16.6,a,a),a)' ) &
1702 0 : 'In a non-metallic case (occopt<3), for a RF calculation,',ch10,&
1703 0 : 'if the eigenvalues are degenerate,',' the occupation numbers must also be degenerate.',ch10,&
1704 0 : 'However, the following pair of states gave :',ch10,&
1705 0 : 'k -state, band number',ibandk,', occ=',occ_k(ibandk),'eigenvalue=',eig0_k(ibandk),',',ch10,&
1706 0 : ' kq-state, band number',ibandkq,', occ=',occ_kq(ibandkq),', eigenvalue=',eig0_kq(ibandkq),'.',ch10,&
1707 0 : 'Action: change occopt, consistently, in GS and RF calculations.'
1708 0 : ABI_ERROR(msg)
1709 : end if
1710 : ratio=0.0_dp
1711 : else
1712 : ! In the metallic case, one can compute a better approximation of the
1713 : ! ratio by using derivatives doccde
1714 563224 : ratio=0.5_dp*(doccde_kq(ibandkq)+doccde_k(ibandk))
1715 : ! write(std_out,*)' occeig : ibandkq,doccde_kq(ibandkq)',ibandkq,doccde_kq(ibandkq)
1716 : ! write(std_out,*)' ibandk ,doccde_k (ibandk )',ibandk,doccde_k(ibandk)
1717 : end if
1718 : end if
1719 :
1720 : ! Here, must pay attention to the smallness of some coefficient
1721 35485018 : diffabsocc=abs(occ_k(ibandk))-abs(occ_kq(ibandkq))
1722 35485018 : sumabsocc=abs(occ_k(ibandk))+abs(occ_kq(ibandkq))
1723 40018604 : if(sumabsocc>tol8)then
1724 25770634 : if( diffabsocc > sumabsocc*tol5 ) then
1725 4057275 : rocceig(ibandkq,ibandk)=ratio
1726 21713359 : else if ( diffabsocc >= -sumabsocc*tol5 ) then
1727 17655526 : rocceig(ibandkq,ibandk)=0.5_dp*ratio
1728 : else
1729 4057833 : rocceig(ibandkq,ibandk)=0.0_dp
1730 : end if
1731 : end if
1732 :
1733 : end do ! ibandkq
1734 : end do ! ibandk
1735 :
1736 855197 : end subroutine occeig
1737 : !!***
1738 :
1739 : !----------------------------------------------------------------------
1740 :
1741 : !!****f* m_occ/occ_fd
1742 : !! NAME
1743 : !! occ_fd
1744 : !!
1745 : !! FUNCTION
1746 : !! Fermi-Dirac statistics: 1 / [(exp((e - mu)/ KT) + 1]
1747 : !! Note that occ_fs in [0, 1] so the spin factor is not included, unlike the
1748 : !! occupations stored in ebands%occ.
1749 : !!
1750 : !! INPUTS
1751 : !! ee=Single particle energy in Ha
1752 : !! kT=Value of K_Boltzmann x T in Ha.
1753 : !! mu=Chemical potential in Ha.
1754 : !!
1755 : !! SOURCE
1756 :
1757 41421098 : elemental real(dp) function occ_fd(ee, kT, mu)
1758 :
1759 : !Arguments ------------------------------------
1760 : real(dp),intent(in) :: ee, kT, mu
1761 :
1762 : !Local variables ------------------------------
1763 : real(dp) :: ee_mu,arg
1764 : ! *************************************************************************
1765 :
1766 41421098 : ee_mu = ee - mu
1767 :
1768 : ! 1 kelvin [K] = 3.16680853419133E-06 Hartree
1769 41421098 : if (kT > tol6) then
1770 23910079 : arg = ee_mu / kT
1771 23910079 : if (arg > maxFDarg) then
1772 : occ_fd = zero
1773 15098260 : else if (arg < -maxFDarg) then
1774 : occ_fd = one
1775 : else
1776 6984950 : occ_fd = one / (exp(arg) + one)
1777 : end if
1778 : else
1779 : ! Heaviside
1780 17511019 : if (ee_mu > zero) then
1781 : occ_fd = zero
1782 8791513 : else if (ee_mu < zero) then
1783 : occ_fd = one
1784 : else
1785 41421098 : occ_fd = half
1786 : end if
1787 : end if
1788 :
1789 41421098 : end function occ_fd
1790 : !!***
1791 :
1792 : !----------------------------------------------------------------------
1793 :
1794 : !!****f* m_occ/occ_dfde
1795 : !! NAME
1796 : !! occ_dfde
1797 : !!
1798 : !! FUNCTION
1799 : !! Derivative of Fermi-Dirac statistics: - (exp((e - mu)/ KT) / KT[(exp((e - mu)/ KT) + 1]^2
1800 : !! Note that kT is given in Hartree so the derivative as well
1801 : !!
1802 : !! INPUTS
1803 : !! ee=Single particle energy in Ha
1804 : !! kT=Value of K_Boltzmann x T in Ha.
1805 : !! mu=Chemical potential in Ha.
1806 : !!
1807 : !! SOURCE
1808 :
1809 204693600 : elemental real(dp) function occ_dfde(ee, kT, mu)
1810 :
1811 : !Arguments ------------------------------------
1812 : real(dp),intent(in) :: ee, kT, mu
1813 :
1814 : !Local variables ------------------------------
1815 : real(dp) :: ee_mu,arg
1816 : ! *************************************************************************
1817 :
1818 204693600 : ee_mu = ee - mu
1819 :
1820 : ! 1 kelvin [K] = 3.16680853419133E-06 Hartree
1821 204693600 : if (kT > tol6) then
1822 102736146 : arg = ee_mu / kT
1823 102736146 : if (arg > maxDFDarg) then
1824 : occ_dfde = zero
1825 57105075 : else if (arg < -maxDFDarg) then
1826 : occ_dfde = zero
1827 : else
1828 11474029 : occ_dfde = - exp(arg) / (exp(arg) + one)**2 / kT
1829 : end if
1830 : else
1831 : occ_dfde = zero
1832 : end if
1833 :
1834 204693600 : end function occ_dfde
1835 : !!***
1836 :
1837 : !----------------------------------------------------------------------
1838 :
1839 : !!****f* m_occ/occ_be
1840 : !! NAME
1841 : !! occ_be
1842 : !!
1843 : !! FUNCTION
1844 : !! Bose-Einstein statistics 1 / [(exp((e - mu)/ KT) - 1]
1845 : !!
1846 : !! INPUTS
1847 : !! ee=Single particle energy in Ha
1848 : !! kT=Value of K_Boltzmann x T in Ha.
1849 : !! mu=Chemical potential in Ha (usually zero)
1850 : !!
1851 : !! SOURCE
1852 :
1853 4470992 : elemental real(dp) function occ_be(ee, kT, mu)
1854 :
1855 : !Arguments ------------------------------------
1856 : real(dp),intent(in) :: ee, kT, mu
1857 :
1858 : !Local variables ------------------------------
1859 : real(dp) :: ee_mu, arg
1860 : ! *************************************************************************
1861 :
1862 4470992 : ee_mu = ee - mu
1863 :
1864 : ! 1 kelvin [K] = 3.16680853419133E-06 Hartree
1865 4470992 : if (kT > tol12) then
1866 4193278 : arg = ee_mu / kT
1867 4193278 : if (arg > tol12 .and. arg < maxBEarg) then
1868 4193278 : occ_be = one / (exp(arg) - one)
1869 : else
1870 : occ_be = zero
1871 : end if
1872 : else
1873 : ! No condensate for T --> 0
1874 : occ_be = zero
1875 : end if
1876 :
1877 4470992 : end function occ_be
1878 : !!***
1879 :
1880 : !----------------------------------------------------------------------
1881 :
1882 : !!****f* m_occ/occ_dbe
1883 : !! NAME
1884 : !! occ_dbe
1885 : !!
1886 : !! FUNCTION
1887 : !! Derivative of Bose-Einstein statistics (exp((e - mu)/ KT) / KT[(exp((e - mu)/ KT) - 1]^2
1888 : !! Note that kT is given in Hartree so the derivative as well
1889 : !!
1890 : !! INPUTS
1891 : !! ee=Single particle energy in Ha
1892 : !! kT=Value of K_Boltzmann x T in Ha.
1893 : !! mu=Chemical potential in Ha (usually zero)
1894 : !!
1895 : !! SOURCE
1896 :
1897 0 : elemental real(dp) function occ_dbe(ee, kT, mu)
1898 :
1899 : !Arguments ------------------------------------
1900 : real(dp),intent(in) :: ee, kT, mu
1901 :
1902 : !Local variables ------------------------------
1903 : real(dp) :: ee_mu, arg
1904 : ! *************************************************************************
1905 :
1906 0 : ee_mu = ee - mu
1907 :
1908 : ! 1 kelvin [K] = 3.16680853419133E-06 Hartree
1909 0 : if (kT > tol12) then
1910 0 : arg = ee_mu / kT
1911 0 : if (arg > tol12 .and. arg < maxDBEarg) then
1912 0 : occ_dbe = exp(arg) / (kT * (exp(arg) - one)**2)
1913 : else
1914 : occ_dbe = zero
1915 : end if
1916 : else
1917 : ! No condensate for T --> 0
1918 : occ_dbe = zero
1919 : end if
1920 :
1921 0 : end function occ_dbe
1922 : !!***
1923 :
1924 : !----------------------------------------------------------------------
1925 :
1926 : !!****f* m_occ/dos_hdr_write
1927 : !!
1928 : !! NAME
1929 : !! dos_hdr_write
1930 : !!
1931 : !! FUNCTION
1932 : !! Write the header of the DOS files, for both smearing and tetrahedron methods.
1933 : !!
1934 : !! INPUTS
1935 : !! deltaene=increment of DOS energy arguments
1936 : !! enemax=maximal value of the DOS energy argument
1937 : !! enemin=minimal value of the DOS energy argument
1938 : !! nene=number of DOS energy argument
1939 : !! eigen(mband*nkpt*nsppol)=eigenvalues (input or init to large number), hartree
1940 : !! fermie=fermi energy useful for band alignment...
1941 : !! fermih= fermi energy of thermalized excited holes when occopt = 9
1942 : !! mband=maximum number of bands
1943 : !! nband(nkpt*nsppol)=number of bands at each k point
1944 : !! nkpt=number of k points
1945 : !! nsppol=1 for unpolarized, 2 for spin-polarized
1946 : !! occopt=option for occupancies, or re-smearing scheme if dblsmr /= 0
1947 : !! prtdos=1 for smearing technique, 2 or 3 for tetrahedron technique
1948 : !! tphysel="physical" electronic temperature with FD occupations
1949 : !! tsmear=smearing width (or temperature)
1950 : !! unitdos=unit number of output of the DOS.
1951 : !!
1952 : !! OUTPUT
1953 : !! Only writing.
1954 : !!
1955 : !! SOURCE
1956 :
1957 76 : subroutine dos_hdr_write(deltaene,eigen,enemax,enemin,fermie,fermih,mband,nband,nene,&
1958 : nkpt,nsppol,occopt,prtdos,tphysel,tsmear,unitdos)
1959 :
1960 : !Arguments ------------------------------------
1961 : !scalars
1962 : integer,intent(in) :: mband,nkpt,nsppol,occopt,prtdos,unitdos,nene
1963 : real(dp),intent(in) :: fermie,fermih,tphysel,tsmear
1964 : real(dp),intent(in) :: deltaene,enemax,enemin
1965 : !arrays
1966 : integer,intent(in) :: nband(nkpt*nsppol)
1967 : real(dp),intent(in) :: eigen(mband*nkpt*nsppol)
1968 :
1969 : !Local variables-------------------------------
1970 : character(len=500) :: msg
1971 : ! *************************************************************************
1972 :
1973 : ! Write the DOS file
1974 76 : write(msg, '(7a,i2,a,i5,a,i4)' ) "#",ch10, &
1975 76 : '# ABINIT package : DOS file ',ch10,"#",ch10,&
1976 152 : '# nsppol =',nsppol,', nkpt =',nkpt,', nband(1)=',nband(1)
1977 76 : call wrtout(unitdos, msg)
1978 :
1979 76 : if (any(prtdos== [1,4])) then
1980 : write(msg, '(a,i2,a,f6.3,a,f6.3,a)' ) &
1981 9 : '# Smearing technique, occopt =',occopt,', tsmear=',tsmear,' Hartree, tphysel=',tphysel,' Hartree'
1982 : else
1983 67 : write(msg, '(a)' ) '# Tetrahedron method '
1984 : end if
1985 76 : call wrtout(unitdos, msg)
1986 :
1987 76 : if (mband*nkpt*nsppol>=3) then
1988 76 : write(msg, '(a,3f8.3,2a)' )'# For identification : eigen(1:3)=',eigen(1:3),ch10,"#"
1989 : else
1990 0 : write(msg, '(a,3f8.3)' ) '# For identification : eigen=',eigen
1991 0 : write(msg, '(3a)')trim(msg),ch10,"#"
1992 : end if
1993 76 : call wrtout(unitdos, msg)
1994 :
1995 76 : if (occopt == 9) then
1996 0 : write(msg, '(a,f16.8, f16.8)' ) '# Fermi energy for electrons and holes ', fermie, fermih
1997 : else
1998 76 : write(msg, '(a,f16.8)' ) '# Fermi energy : ', fermie
1999 : end if
2000 76 : call wrtout(unitdos, msg)
2001 :
2002 76 : if (prtdos==1) then
2003 9 : write(msg, '(5a)' ) "#",ch10,&
2004 9 : '# The DOS (in electrons/Hartree/cell) and integrated DOS (in electrons/cell),',&
2005 18 : ch10,'# as well as the DOS with tsmear halved and doubled, are computed,'
2006 :
2007 67 : else if (prtdos==2)then
2008 4 : write(msg, '(3a)' ) "#",ch10,&
2009 8 : '# The DOS (in electrons/Hartree/cell) and integrated DOS (in electrons/cell) are computed,'
2010 :
2011 63 : else if (any(prtdos == [3, 4])) then
2012 62 : write(msg, '(5a)' ) "#",ch10,&
2013 62 : '# The local DOS (in electrons/Hartree for one atomic sphere)',ch10,&
2014 124 : '# and integrated local DOS (in electrons for one atomic sphere) are computed.'
2015 :
2016 1 : else if (prtdos==5)then
2017 1 : write(msg, '(9a)' ) "#",ch10,&
2018 1 : '# The spin component DOS (in electrons/Hartree/cell)',ch10,&
2019 1 : '# and integrated spin component DOS (in electrons/cell) are computed.',ch10,&
2020 1 : '# Remember that the wf are eigenstates of S_z and S^2, not S_x and S_y',ch10,&
2021 2 : '# so the latter will not always sum to 0 for paired electronic states.'
2022 : end if
2023 76 : call wrtout(unitdos, msg)
2024 :
2025 : write(msg, '(a,i5,a,a,a,f9.4,a,f9.4,a,f8.5,a,a,a)' )&
2026 76 : '# at ',nene,' energies (in Hartree) covering the interval ',ch10,&
2027 152 : '# between ',enemin,' and ',enemax,' Hartree by steps of ',deltaene,' Hartree.',ch10,"#"
2028 76 : call wrtout(unitdos, msg)
2029 :
2030 76 : if (prtdos==1) then
2031 : write(msg, '(a,a)' )&
2032 9 : '# energy DOS Integr. DOS ',' DOS DOS '
2033 9 : call wrtout(unitdos,msg)
2034 :
2035 : write(msg, '(a)' )&
2036 9 : '# (tsmear/2) (tsmear*2) '
2037 9 : call wrtout(unitdos,msg)
2038 : else
2039 67 : write(msg, '(a)' ) '# energy DOS '
2040 : end if
2041 :
2042 76 : end subroutine dos_hdr_write
2043 : !!***
2044 :
2045 : !!****f* m_occ/get_fact_spin_tol_empty
2046 : !! NAME
2047 : !! get_fact_spin_tol_empty
2048 : !!
2049 : !! FUNCTION
2050 : !!
2051 : !! INPUTS
2052 : !!
2053 : !! OUTPUT
2054 : !!
2055 : !! SOURCE
2056 :
2057 1118 : subroutine get_fact_spin_tol_empty(nsppol, nspinor, tol_empty_in, fact_spin, tol_empty)
2058 :
2059 : integer,intent(in) :: nsppol, nspinor
2060 : real(dp),intent(in) :: tol_empty_in
2061 : real(dp),intent(out) :: fact_spin, tol_empty
2062 :
2063 2222 : select case (nsppol)
2064 : case (1)
2065 1104 : fact_spin = half; tol_empty = tol_empty_in ! below this value the state is assumed empty
2066 1104 : if (nspinor == 2) then
2067 17 : fact_spin = one; tol_empty = half * tol_empty_in ! below this value the state is assumed empty
2068 : end if
2069 : case (2)
2070 14 : fact_spin = one; tol_empty = half * tol_empty_in ! to be consistent and obtain similar results if a metallic
2071 : case default ! spin unpolarized system is treated using nsppol==2
2072 1118 : ABI_BUG(sjoin('Wrong nsppol:', itoa(nsppol)))
2073 : end select
2074 :
2075 1118 : end subroutine get_fact_spin_tol_empty
2076 : !!***
2077 :
2078 : end module m_occ
2079 : !!***
|
