Branch data 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 : 552699 : subroutine getnel(doccde,dosdeltae,eigen,entropy,fermie,fermih,maxocc,mband,nband,&
126 : 552699 : 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 : 552699 : real(dp),allocatable :: entfun(:,:),occfun(:,:), smdfun(:,:),xgrid(:)
164 : 552699 : real(dp),allocatable :: arg(:),derfun(:),dos(:),dosdble(:),doshalf(:),ent(:), intdos(:)
165 : 552699 : real(dp),allocatable :: occ_tmp(:),ent_tmp(:),doccde_tmp(:)
166 : 552699 : real(dp),allocatable :: dos_core(:),dosdble_core(:),doshalf_core(:),intdos_core(:)
167 : 552699 : real(dp),allocatable :: occ_tmp_core(:),doccde_tmp_core(:),arg_core(:)
168 : 552699 : real(dp),allocatable :: ent_core(:), derfun_core(:)
169 : : ! *************************************************************************
170 : :
171 : : !call cwtime(cpu, wall, gflops, "start")
172 : :
173 [ - + ]: 552699 : 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 [ + + ]: 552699 : 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 : 551699 : low_band_index = 1
186 : 551699 : high_band_index = nband(1)
187 [ + + ]: 8523024 : number_of_bands = sum(nband(:))
188 : : end if
189 [ - + + - : 1658097 : ABI_MALLOC(occ_tmp,(number_of_bands))
+ - + - -
+ - + ]
190 [ - + + - : 1105398 : ABI_MALLOC(ent_tmp,(number_of_bands))
- + - + ]
191 [ - + + - : 1105398 : ABI_MALLOC(doccde_tmp,(number_of_bands))
- + - + ]
192 : :
193 : : ! Just get the number nptsdiv2 and allocate entfun, occfun, smdfun and xgrid accordingly
194 : 552699 : nptsdiv2 = nptsdiv2_def
195 : :
196 [ - + ]: 552699 : ABI_MALLOC(entfun,(-nptsdiv2:nptsdiv2,2))
197 [ - + ]: 552699 : ABI_MALLOC(occfun,(-nptsdiv2:nptsdiv2,2))
198 [ - + ]: 552699 : ABI_MALLOC(smdfun,(-nptsdiv2:nptsdiv2,2))
199 [ - + ]: 552699 : ABI_MALLOC(xgrid,(-nptsdiv2:nptsdiv2))
200 : :
201 : 552699 : 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 [ - + + - : 1105398 : ABI_MALLOC(arg,(number_of_bands))
- + - + ]
208 [ - + + - : 1105398 : ABI_MALLOC(derfun,(number_of_bands))
- + - + ]
209 [ - + + - : 1105398 : ABI_MALLOC(ent,(number_of_bands))
- + - + ]
210 [ + + ]: 552699 : if (option==1) then
211 : : ! normal evaluation of occupations and entropy
212 : :
213 : : index = 0
214 : : index_tot = 0
215 [ + + ]: 1187964 : do isppol=1,nsppol
216 [ + + ]: 9169113 : do ikpt=1,nkpt
217 [ - + ]: 7981149 : if (occopt == 2) high_band_index=nband(ikpt+nkpt*(isppol-1))
218 [ + + ]: 92195230 : do iband=low_band_index,high_band_index
219 : 84214081 : index = index + 1
220 [ + + ]: 92195230 : if (tsmear==0) then
221 : 262533 : arg(index) = sign(huge_tsmearinv,fermie-eigen(index_tot + iband))
222 : : else
223 : 83951548 : arg(index)=(fermie-eigen(index_tot + iband))*tsmearinv
224 : : end if
225 : : end do
226 : 8616423 : 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 : 552690 : call splfit(xgrid, doccde_tmp, occfun, 1, arg, occ_tmp, (2*nptsdiv2+1), number_of_bands)
237 : 552690 : 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 [ + + ]: 552690 : if(present(extfpmd_nbdbuf)) then
246 : : index=0
247 : : index_tot=0
248 [ + + ]: 1184947 : do isppol=1,nsppol
249 [ + + ]: 9109212 : do ikpt=1,nkpt
250 [ - + ]: 7924265 : if (occopt == 2) high_band_index=nband(ikpt+nkpt*(isppol-1))
251 [ + + ]: 91795882 : do iband=low_band_index,high_band_index
252 : 83871617 : index=index+1
253 [ + + ]: 91795882 : 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 : 8558016 : index_tot=index_tot+nband(ikpt+nkpt*(isppol-1))
259 : : end do
260 : : end do
261 : : end if
262 : :
263 : 552690 : nelect=zero; entropy=zero
264 : 552690 : index=0
265 : 552690 : index_tot = 0
266 [ + + ]: 1187964 : do isppol=1,nsppol
267 [ + + ]: 9169113 : do ikpt=1,nkpt
268 : 7981149 : wk = wtk(ikpt)
269 [ - + ]: 7981149 : if (occopt == 2) high_band_index=nband(ikpt+nkpt*(isppol-1))
270 [ + + ]: 92195230 : do iband=low_band_index,high_band_index
271 : 84214081 : index = index + 1
272 : 84214081 : ent(index) = ent(index)*maxocc
273 : 84214081 : occ(iband + index_tot) = occ_tmp(index)*maxocc
274 : 84214081 : doccde(iband + index_tot) = -doccde_tmp(index)*maxocc*tsmearinv
275 : 84214081 : entropy = entropy + wk*ent(index)
276 : 92195230 : nelect = nelect + wk*occ(iband + index_tot)
277 : : end do
278 : 8616423 : index_tot = index_tot + nband(ikpt+nkpt*(isppol-1))
279 : : end do
280 : : end do
281 : :
282 [ + + ]: 552690 : if(present(rcpaw)) then
283 [ + + ]: 542074 : 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 [ - + ]: 552699 : ABI_FREE(arg)
491 [ - + ]: 552699 : ABI_FREE(derfun)
492 [ - + ]: 552699 : ABI_FREE(ent)
493 [ - + ]: 552699 : ABI_FREE(entfun)
494 [ - + ]: 552699 : ABI_FREE(occfun)
495 [ - + ]: 552699 : ABI_FREE(smdfun)
496 [ - + ]: 552699 : ABI_FREE(xgrid)
497 [ - + ]: 552699 : ABI_FREE(occ_tmp)
498 [ - + ]: 552699 : ABI_FREE(doccde_tmp)
499 [ - + ]: 552699 : ABI_FREE(ent_tmp)
500 : :
501 : : !call cwtime_report(" getnel", cpu, wall, gflops, end_str=ch10)
502 : :
503 : 552699 : 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 : 14255 : subroutine newocc(doccde,eigen,entropy,fermie,fermih,ivalence,spinmagntarget,mband,nband,&
549 : 14255 : 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 : 14255 : 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,tol_nel
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 : 14255 : 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 : 14255 : call timab(74,1,tsec)
591 : :
592 : 14255 : rcpaw_nelect=zero
593 [ + + ]: 14255 : if(present(rcpaw)) then
594 [ + + ]: 13715 : 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 [ - + ]: 14255 : 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 [ + + ]: 30321 : do isppol=1,nsppol
612 [ + + ]: 221399 : do ikpt=1,nkpt
613 [ - + ]: 207144 : 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 [ - + ]: 14255 : 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 [ + + ]: 14255 : 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 : 14255 : 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 [ - + ]: 14255 : 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 [ + + - + ]: 14255 : 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 [ + + + - : 14255 : 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 [ + + ]: 14255 : if(present(extfpmd)) then
667 [ + + ]: 13739 : 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 : 14255 : dosdeltae = zero ! the DOS is not computed, with option=1
684 [ + - - + : 2079114 : fermie_lo = minval(eigen(1:nband(1)*nkpt*nsppol)) - 6.001_dp * tsmear ! fermi_lo ->fermie_lo
- - + + ]
685 [ + + ]: 14255 : if (occopt == 3 .or. occopt==9) fermie_lo = fermie_lo - 24.0_dp * tsmear
686 [ + + ]: 14255 : if(occopt==9) fermih_lo = fermie_lo ! Take into account holes
687 : :
688 [ + + ]: 14255 : 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 : 14245 : & 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 [ + - - + : 2079114 : fermie_hi = maxval(eigen(1:nband(1)*nkpt*nsppol)) + 6.001_dp * tsmear
- - + + ]
700 : : ! Safety value
701 : 14255 : fermie_hi = min(fermie_hi, 1.e6_dp)
702 [ + + ]: 14255 : if(occopt == 3 .or. occopt == 9) fermie_hi = fermie_hi + 24.0_dp * tsmear
703 [ + + ]: 14255 : if(occopt == 9) fermih_hi=fermie_hi
704 : :
705 [ + + ]: 14255 : 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 : 14245 : & 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 [ + + ]: 14255 : if(present(extfpmd)) then
720 [ + + ]: 13739 : 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 [ + + ]: 14255 : 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 [ + - - + : 14255 : 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 [ + + ]: 14255 : 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 [ + + ]: 14255 : if( abs(spinmagntarget+99.99_dp) < tol10) then
784 : :
785 : : ! The bisection below stops at the first mu with |N(mu) - nelect| < tol_nel * nelect.
786 : : ! Only the e-ph codes call newocc with tsmear = 0, i.e. at T = 0. Then N(mu) = nelect for every mu in a gap,
787 : : ! but getnel sums nkpt*nelect/maxocc nonzero terms, each with a round-off error of up to epsilon * nelect.
788 : : ! For k-point weights that are not exact in binary, such as those of a 24^3 mesh, this error exceeds tol14 * nelect,
789 : : ! so with tol_nel = tol14 no mu in the gap meets the stopping condition above, and the bisection wrongly ends at a
790 : : ! band edge (e.g. mu = VBM - 2e-15 Ha), where round-off decides whether the band-edge state is occupied.
791 : : ! At T = 0 we therefore take tol_nel = epsilon * (number of nonzero terms), but not less than tol14.
792 : 14123 : tol_nel = tol14
793 [ + + ]: 14123 : if (tsmear == zero) tol_nel = max(tol14, epsilon(one) * nkpt * nelect / maxocc)
794 : :
795 : : ! Usual bisection loop
796 [ + - ]: 514064 : do ii=1,niter_max
797 : 514064 : fermie_mid = (fermie_hi + fermie_lo) * half
798 [ + + ]: 514064 : if (occopt == 9) fermih_mid=(fermih_hi+fermih_lo)*half
799 : : ! Produce nelectmid from fermimid
800 [ + + ]: 514064 : if (occopt /= 9) then
801 : :
802 : : call getnel(doccde,dosdeltae,eigen,entropye,fermie_mid,fermie_mid,maxocc,mband,nband,&
803 : : nelectmid,nkpt,nsppol,occ,occopt,option1,tphysel,tsmear,fake_unit,wtk, 1, nband(1),&
804 : 513584 : extfpmd_nbdbuf=extfpmd_nbdbuf,rcpaw=rcpaw_getnel)
805 : :
806 : : ! Compute the number of free electrons of the extfpmd model
807 : : ! with corresponding chemical potential and add to nelect bounds.
808 [ + + ]: 513584 : if(present(extfpmd)) then
809 [ + + ]: 498025 : if(associated(extfpmd)) then
810 : 4461 : call extfpmd%compute_nelect(fermie_mid,nband,nelectmid,nkpt,nspinor,nsppol,wtk)
811 : : end if
812 : : end if
813 : :
814 : : !write(std_out,'(a,i0,1x, 3(a,es13.5))' ) " iter: ", ii, &
815 : : ! ' fermi_mid: ',fermimid * Ha_eV, ', n_mid: ',nelectmid, &
816 : : ! ", (n_mid-nelect)/nelect: ", (nelectmid - nelect) / nelect
817 : :
818 : : !if (nelectmid > nelect * (one - tol)) then
819 : : ! fermihi = fermimid
820 : : ! nelecthi = nelectmid
821 : : !end if
822 : : !if (nelectmid < nelect * (one + tol)) then
823 : : ! fermilo = fermimid
824 : : ! nelectlo = nelectmid
825 : : !end if
826 [ + + ]: 513584 : if(nelectmid>nelect*(one-tol_nel))then
827 : 256184 : fermie_hi=fermie_mid
828 : 256184 : nelecthi=nelectmid
829 : : end if
830 [ + + ]: 513584 : if(nelectmid<nelect*(one+tol_nel))then
831 : 269440 : fermie_lo=fermie_mid
832 : 269440 : nelectlo=nelectmid
833 : : end if
834 : :
835 : : else
836 : :
837 : : call getnel(doccde,dosdeltae,eigen,entropye,fermie_mid,fermie_mid,maxocc,mband,nband,&
838 : 480 : nelectmid,nkpt,nsppol,occ,occopt,option1,tphysel,tsmear,fake_unit,wtk, ivalence+1, nband(1))
839 : : call getnel(doccde,dosdeltae,eigen,entropyh,fermih_mid,fermih_mid,maxocc,mband,nband,&
840 : 480 : nholesmid,nkpt,nsppol,occ,occopt,option1,tphysel,tsmear,fake_unit,wtk,1,ivalence)
841 : :
842 [ + + ]: 480 : if(nelectmid>ne_qFD*(one-tol14))then
843 : 260 : fermie_hi = fermie_mid
844 : 260 : nelecthi = nelectmid
845 [ + - ]: 220 : else if (nelectmid<ne_qFD*(one-tol14))then
846 : 220 : fermie_lo = fermie_mid
847 : 220 : nelectlo = nelectmid
848 : : end if
849 [ + + ]: 480 : if(nholesmid>(nelect-nh_qFD)*(one-tol14))then
850 : 266 : fermih_hi = fermih_mid
851 : 266 : nholeshi = nholesmid
852 [ + - ]: 214 : else if(nholesmid<(nelect-nh_qFD)*(one+tol14))then
853 : 214 : fermih_lo = fermih_mid
854 : 214 : nholeslo = nholesmid
855 : : end if
856 : :
857 : : end if
858 : :
859 : : !if (abs(nelectmid - nelect) <= nelect*two*tol) exit
860 : : !write(std_out,'(2(a,es13.5))' )' bisection move: fermi_lo: ',fermilo * Ha_eV,", fermi_hi: ", fermihi * Ha_eV
861 : :
862 : : ! if (abs(nelecthi - nelectlo) <= nelect*two*tol .or. &
863 : : ! abs(fermihi - fermilo) <= tol * abs(fermihi + fermilo) ) exit
864 [ + + ]: 514064 : if (occopt /= 9) then
865 [ + + + + ]: 513584 : if( abs(nelecthi-nelectlo) <= nelect*two*tol14 .or. abs(fermie_hi-fermie_lo) <= tol14*abs(fermie_hi+fermie_lo) ) exit
866 : : else
867 : : if( ( abs(nelecthi-nelectlo) <= ne_qFD*two*tol14 .or. &
868 [ + - + + : 480 : & abs(fermie_hi-fermie_lo) <= tol14*abs(fermie_hi+fermie_lo) ) .and. &
+ - - - ]
869 : : ( abs(nholeshi-nholeslo) <= (nelect-nh_qFD)*two*tol14 .or. &
870 : : & abs(fermih_hi-fermih_lo) <= tol14*abs(fermih_hi+fermih_lo) ) ) exit
871 : : end if
872 : :
873 [ - + ]: 514064 : if (ii == niter_max) then
874 : : write(msg,'(a,i0,3a,es22.14,a,es22.14,a)')&
875 : 0 : 'It was not possible to find Fermi energy in ',niter_max,' max bisections.',ch10,&
876 : 0 : 'nelecthi: ',nelecthi,', and nelectlo: ',nelectlo,'.'
877 : 0 : ABI_BUG(msg)
878 [ # # ]: 0 : if (occopt == 9) then
879 : : write(msg,'(a,es22.14,a,es22.14,a)')&
880 : 0 : 'nholesi = ',nholeshi,', and holeslo = ',nholeslo,'.'
881 : : end if
882 : : end if
883 : : end do ! End of bisection loop
884 : :
885 : 14123 : fermie = fermie_mid
886 : 14123 : entropy= entropye
887 : :
888 [ + + ]: 14123 : if (occopt /= 9) then
889 : : write(msg, '(2(a,f14.6),a,i0)' ) &
890 : 14113 : & ' newocc: new Fermi energy is ',fermie,' , with nelect=',nelectmid,', Number of bisection calls: ',ii
891 : : else
892 : 10 : fermih=fermih_mid
893 : 10 : entropy = entropy + entropyh ! CP: adding entropy of the holes subsystem
894 : : write(msg, '(2(a,f14.6),a,i0)' ) &
895 : 10 : & ' newocc: new Fermi energy for excited electrons is ',fermie,' , with ne_qFD=',nelectmid,', Number of bisection calls: ',ii
896 : 10 : call wrtout(std_out,msg,'COLL')
897 : : write(msg, '(2(a,f14.6),a,i0)' ) &
898 : 10 : & ' newocc: new Fermi energy for excited holes is ',fermih,' , with nh_qFD=',nelect-nholesmid,&
899 : 20 : & ', Number of bisection calls: ',ii
900 : : end if
901 : 14123 : call wrtout(std_out,msg)
902 : :
903 : : ! Compute occupation numbers for prtstm/=0, close to the Fermi energy
904 [ + + ]: 14123 : if (present(stmbias)) then
905 : :
906 [ + + ]: 13640 : if (abs(stmbias) > tol10) then
907 : :
908 : : ! Prevent use with occopt = 9 so far
909 : : ! XG220804: This test is not needed, as prtstm/=0 must be used with occopt==7, as tested in chkinp.F90
910 [ - + ]: 3 : if (occopt == 9) then
911 : 0 : ABI_ERROR('Occopt 9 and prtstm /=0 not implemented together. Change occopt or prtstm.')
912 : : end if
913 : :
914 : 3 : fermie_biased = fermie - stmbias
915 [ - + + - : 9 : ABI_MALLOC(occt,(mband*nkpt*nsppol))
+ - + - -
+ - + ]
916 : :
917 : : call getnel(doccde,dosdeltae,eigen,entropy,fermie_biased,fermie_biased,maxocc,mband,nband,&
918 : : & nelect_biased,nkpt,nsppol,occt,occopt,option1,tphysel,tsmear,fake_unit,wtk,1,nband(1),&
919 : 3 : & extfpmd_nbdbuf=extfpmd_nbdbuf)
920 [ + + ]: 273 : occ(:)=occ(:)-occt(:)
921 : :
922 : : ! Possibly filter a specific band contribution
923 [ + - ]: 3 : if (present(prtstm)) then
924 [ + + ]: 3 : if (prtstm < 0)then
925 : : ibantot=1
926 [ + + ]: 2 : do isppol=1,nsppol
927 [ + + ]: 8 : do ikpt=1,nkpt
928 : 6 : nban=nband(ikpt+(isppol-1)*nkpt)
929 [ + + ]: 97 : do iban=1,nban
930 [ + + ]: 90 : if(iban/=abs(prtstm)) occ(ibantot)=zero
931 : 96 : ibantot=ibantot+1
932 : : end do ! iban
933 : : end do ! ikpt
934 : : end do ! isppol
935 : : end if ! prtstm < 0
936 : : end if ! present(prtstm)
937 : :
938 : 3 : nelect_biased = abs(nelectmid - nelect_biased)
939 : : ! Here, arrange to have globally positive occupation numbers, irrespective of the stmbias sign
940 [ + + + + ]: 93 : if (-stmbias > tol10) occ(:) = -occ(:)
941 [ - + ]: 3 : ABI_FREE(occt)
942 : :
943 : 3 : write(msg,'(a,f14.6)')' newocc: the number of electrons in the STM range is nelect_biased=',nelect_biased
944 : 3 : call wrtout(std_out,msg)
945 : : end if
946 : : endif ! present(stmbias)
947 : :
948 : : else
949 : : ! Calculations with a specified moment
950 : : ! Bisection loop
951 : 132 : cnt2=0
952 : 132 : cnt3=0
953 : 132 : entropy=zero
954 : 132 : maxocc=one
955 [ - + + - : 396 : ABI_MALLOC(doccdet,(nkpt*mband))
+ - + - -
+ - + ]
956 [ - + + - : 264 : ABI_MALLOC(eigent,(nkpt*mband))
- + - + ]
957 [ - + + - : 264 : ABI_MALLOC(occt,(nkpt*mband))
- + - + ]
958 [ - + + - : 396 : ABI_MALLOC(nbandt,(nkpt))
+ - + - -
+ - + ]
959 : :
960 [ + + ]: 396 : do is = 1, nsppol
961 : 264 : nelect_tmp = nelectt(is)
962 : 264 : fermie_hi = fermie_hit(is) ! CP modify name
963 : 264 : fermie_lo = fermie_lot(is) ! CP modify name
964 : 264 : nelecthi = nelecthit(is)
965 : 264 : nelectlo = nelectlot(is)
966 : : ! write(std_out,'(a,i1,3(f8.4,1x))') "Spin, N(spin):", is, nelect, fermihi, fermilo
967 : : ! write(std_out,'(a,2(f8.4,1x))') "Hi, lo:", nelecthi, nelectlo
968 : :
969 [ + - ]: 9119 : do ii=1,niter_max
970 : 9119 : fermie_mid_tmp=(fermie_hi+fermie_lo)/2.0_dp ! CP modify name
971 : : ! temporary arrays
972 : 9119 : cnt = 0
973 [ + + ]: 25965 : do ik = 1, nkpt
974 : 16846 : nbandt(ik) = mband
975 [ + + ]: 222689 : do ib = 1, mband
976 : 196724 : cnt = cnt + 1
977 : 196724 : eigent(cnt) = eigen(cnt+cnt2)
978 : 196724 : occt(cnt) = occ(cnt+cnt2)
979 : 213570 : doccdet(cnt) = doccde(cnt+cnt2)
980 : : end do
981 : : end do
982 : :
983 : : ! Produce nelectmid from fermimid
984 : : call getnel(doccdet,dosdeltae,eigent,entropy_tmp,fermie_mid_tmp,fermie_mid_tmp,maxocc,mband,nbandt,&
985 : : nelectmid,nkpt,1,occt,occopt,option1,tphysel,tsmear,fake_unit,wtk,1,nband(1),&
986 : 9119 : extfpmd_nbdbuf=extfpmd_nbdbuf)
987 : :
988 : 9119 : entropyet(is) = entropy_tmp
989 : : fermie_midt(is) = fermie_mid_tmp
990 : 9119 : fermie_mid = fermie_midt(is)
991 : :
992 : : ! temporary arrays
993 : 9119 : cnt = 0
994 [ + + ]: 25965 : do ik = 1, nkpt
995 [ + + ]: 222689 : do ib = 1, mband
996 : 196724 : cnt = cnt + 1
997 : 196724 : occ(cnt+cnt2) = occt(cnt)
998 : 213570 : doccde(cnt+cnt2) = doccdet(cnt)
999 : : end do
1000 : : end do
1001 : : ! write(std_out,'(a,es24.16,a,es24.16)' )' newocc: from fermi=',fermimid,', getnel gives nelect=',nelectmid
1002 : :
1003 [ + + ]: 9119 : if(nelectmid>=nelect_tmp)then
1004 : 4474 : fermie_hi=fermie_mid_tmp
1005 : 4474 : nelecthi=nelectmid
1006 : : else
1007 : 4645 : fermie_lo=fermie_mid_tmp
1008 : 4645 : nelectlo=nelectmid
1009 : : end if
1010 [ + + + + ]: 9119 : if( abs(nelecthi-nelectlo) <= 1.0d-13 .or. abs(fermie_hi-fermie_lo) <= 0.5d-14*abs(fermie_hi+fermie_lo) ) exit
1011 : :
1012 [ - + ]: 9119 : if(ii==niter_max)then
1013 : : write(msg,'(a,i3,3a,es22.14,a,es22.14,a)')&
1014 : 0 : 'It was not possible to find Fermi energy in ',niter_max,' bisections.',ch10,&
1015 : 0 : 'nelecthi: ',nelecthi,', and nelectlo: ',nelectlo,'.'
1016 : 0 : ABI_BUG(msg)
1017 : : end if
1018 : : end do ! End of bisection loop
1019 : :
1020 : 264 : cnt2 = cnt2 + nkpt*mband
1021 : 264 : entropy = entropy + entropyet(is)
1022 : 264 : fermie=fermie_mid
1023 : : write(msg, '(a,i2,a,f14.6,a,f14.6,a,a,i4)' ) &
1024 : 264 : ' newocc: new Fermi energy for spin ', is, ' is ',fermie,' , with nelect: ',nelectmid,ch10,&
1025 : 528 : ' Number of bisection calls =',ii
1026 : 396 : call wrtout(std_out,msg)
1027 : :
1028 : : end do ! spin
1029 : :
1030 [ - + ]: 132 : ABI_FREE(doccdet)
1031 [ - + ]: 132 : ABI_FREE(eigent)
1032 [ - + ]: 132 : ABI_FREE(nbandt)
1033 [ - + ]: 132 : ABI_FREE(occt)
1034 : :
1035 : : end if ! End of logical on fixed moment calculations
1036 : :
1037 [ + + ]: 14255 : if(present(rcpaw)) then
1038 [ + + ]: 13715 : if(associated(rcpaw)) then
1039 [ + + ]: 42 : do itypat=1,rcpaw%ntypat
1040 [ + - ]: 42 : if(rcpaw%atm(itypat)%zcore_orig>zero) then
1041 [ - + + + : 126 : rcpaw%atm(itypat)%occ_res=rcpaw%atm(itypat)%occ_res+rcpaw%atm(itypat)%occ
+ + ]
1042 : : endif
1043 : : enddo
1044 : : endif
1045 : : endif
1046 : :
1047 : : !write(std_out,*) "kT*Entropy:", entropy*tsmear
1048 : :
1049 : : ! MG: If you are wondering why this part is now disabled by default consider that this output
1050 : : ! is produced many times in the SCF cycle and in EPH we have to call this routine for
1051 : : ! several temperature and the log becomes unreadable.
1052 : : ! If you really need to look at the occupation factors use prtvol > 0.
1053 : 14255 : nkpt_eff = nkpt
1054 [ + + ]: 14255 : if (prtvol == 0) nkpt_eff = 0
1055 [ + + ]: 14255 : if (prtvol == 1) nkpt_eff = min(nkpt_max, nkpt)
1056 : :
1057 [ + + ]: 14255 : if (nsppol == 1)then
1058 : :
1059 [ + + ]: 12444 : if (nkpt_eff /= 0) then
1060 : 2489 : write(msg, '(a,i0,a)' )' newocc: computed new occ. numbers for occopt= ',occopt,' , spin-unpolarized case. '
1061 : 2489 : call wrtout(std_out,msg)
1062 [ + + ]: 20346 : do ikpt=1,nkpt_eff
1063 : 17857 : write(msg,'(a,i4,a)' ) ' k-point number ',ikpt,' :'
1064 [ + + ]: 45481 : do ii=0,(nband(1)-1)/12
1065 [ + + + - ]: 25299 : if (ii == 3 .and. prtvol /= 0) exit
1066 : 25135 : write(msg,'(12f6.3)') occ(1+ii*12+(ikpt-1)*nband(1):min(12+ii*12,nband(1))+(ikpt-1)*nband(1))
1067 : 42992 : call wrtout(std_out,msg)
1068 : : end do
1069 : : end do
1070 [ + + ]: 2489 : if (nkpt /= nkpt_eff) call wrtout(std_out,' newocc: prtvol=0, stop printing more k-point information')
1071 : :
1072 : : !call wrtout(std_out,' newocc: corresponding derivatives are ')
1073 : : !do ikpt=1,nkpt_eff
1074 : : !write(msg,'(a,i4,a)' ) ' k-point number ',ikpt,' :'
1075 : : !do ii=0,(nband(1)-1)/12
1076 : : !write(msg,'(12f6.1)') doccde(1+ii*12+(ikpt-1)*nband(1):min(12+ii*12,nband(1))+(ikpt-1)*nband(1))
1077 : : !call wrtout(std_out,msg)
1078 : : !end do
1079 : : !end do
1080 : : !if(nkpt/=nkpt_eff)then
1081 : : ! call wrtout(std_out,'newocc: prtvol=0, stop printing more k-point information')
1082 : : !end if
1083 : : end if
1084 : :
1085 : : else
1086 : :
1087 [ + + ]: 1811 : if (nkpt_eff /= 0) then
1088 : 721 : write(msg, '(a,i0,2a)' )' newocc: computed new occupation numbers for occopt= ',occopt,ch10,' (1) spin up values '
1089 : 721 : call wrtout(std_out, msg)
1090 [ + + ]: 6291 : do ikpt=1,nkpt_eff
1091 : 5570 : write(msg,'(a,i0,a)' ) ' k-point number ',ikpt,':'
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)*nband(1):min(12+ii*12,nband(1))+(ikpt-1)*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 : :
1100 : 721 : call wrtout(std_out,' (2) spin down values ')
1101 [ + + ]: 6291 : do ikpt=1,nkpt_eff
1102 [ + + ]: 16282 : do ii=0,(nband(1)-1)/12
1103 [ - + - - ]: 9991 : if (ii == 3 .and. prtvol /= 0) exit
1104 : 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) )
1105 : 15561 : call wrtout(std_out,msg)
1106 : : end do
1107 : : end do
1108 [ + + ]: 721 : if(nkpt/=nkpt_eff) call wrtout(std_out,' newocc: prtvol=0, stop printing more k-point information')
1109 : : end if
1110 : :
1111 : : end if ! End choice based on spin
1112 : :
1113 : 14255 : rcpaw_getnel=>null()
1114 : :
1115 : 14255 : call timab(74,2,tsec)
1116 : :
1117 : : DBG_EXIT("COLL")
1118 : :
1119 : 14255 : end subroutine newocc
1120 : : !!***
1121 : :
1122 : : !!****f* m_occ/init_occ_ent
1123 : : !! NAME
1124 : : !! init_occ_ent
1125 : : !!
1126 : : !! FUNCTION
1127 : : !!
1128 : : !! INPUTS
1129 : : !!
1130 : : !! OUTPUT
1131 : : !!
1132 : : !! SOURCE
1133 : :
1134 : 552699 : subroutine init_occ_ent(entfun,limit,nptsdiv2,occfun,occopt,option,smdfun,tphysel,tsmear,tsmearinv,xgrid)
1135 : :
1136 : : !Arguments ------------------------------------
1137 : : !scalars
1138 : : integer,intent(in) :: occopt,option
1139 : : real(dp),intent(in) :: tphysel,tsmear
1140 : : integer,intent(inout) :: nptsdiv2
1141 : : real(dp),intent(out) :: limit,tsmearinv
1142 : : real(dp),intent(inout) :: entfun(-nptsdiv2:nptsdiv2,2),occfun(-nptsdiv2:nptsdiv2,2)
1143 : : real(dp),intent(inout) :: smdfun(-nptsdiv2:nptsdiv2,2),xgrid(-nptsdiv2:nptsdiv2)
1144 : :
1145 : : !Local variables-------------------------------
1146 : : !scalars
1147 : : integer :: algo,ii,jj,nconvd2,nmaxFD,nminFD
1148 : : integer,save :: dblsmr,occopt_prev=-9999
1149 : : real(dp),save :: convlim,incconv,limit_occ,tphysel_prev=-9999,tsmear_prev=-9999
1150 : : real(dp) :: aa,dsqrpi,encorr,factor
1151 : : real(dp) :: expinc,expx22,expxo2,gauss,increm
1152 : : real(dp) :: resFD1,resFD2,resFD3,resFD4,resmom,resmom1,resmom2
1153 : : real(dp) :: resmom3,resmom4,secmom,smom1,smom2,thdmom,tmom1,tmom2,tmpexpsum
1154 : : real(dp) :: tmpsmdfun,tratio,tt,xx,yp1,ypn
1155 : : character(len=500) :: msg
1156 : : !arrays
1157 : : real(dp),save :: entfun_prev(-nptsdiv2_def:nptsdiv2_def,2),occfun_prev(-nptsdiv2_def:nptsdiv2_def,2)
1158 : : real(dp),save :: smdfun_prev(-nptsdiv2_def:nptsdiv2_def,2),xgrid_prev(-nptsdiv2_def:nptsdiv2_def)
1159 : 552699 : real(dp),allocatable :: entder(:),occder(:),smd1(:),smd2(:)
1160 : 552699 : real(dp),allocatable :: smdder(:),tgrid(:),work(:),workfun(:)
1161 : : ! *************************************************************************
1162 : :
1163 : : ! Initialize the occupation function and generalized entropy function,
1164 : : ! at the beginning, or if occopt changed
1165 : :
1166 [ - + ]: 552699 : if(option==-1)then
1167 : 0 : nptsdiv2 = nptsdiv2_def
1168 : 0 : return
1169 : : end if
1170 : :
1171 [ + + + + : 552699 : if (occopt_prev/=occopt .or. abs(tsmear_prev-tsmear) >tol12 .or. abs(tphysel_prev-tphysel)>tol12) then
+ + ]
1172 : 757 : occopt_prev=occopt
1173 : 757 : tsmear_prev=tsmear
1174 : 757 : tphysel_prev=tphysel
1175 : :
1176 : : ! Check whether input values of tphysel tsmear and occopt are consistent
1177 : 757 : dblsmr = 0
1178 [ + + ]: 757 : if (abs(tphysel)>tol12) then
1179 : : ! Use re-smearing scheme
1180 [ + - ]: 10 : if (abs(tsmear)>tol12) then
1181 : 10 : dblsmr = 1
1182 : : ! Use FD occupations (one smearing) only with "physical" temperature tphysel
1183 : : ! CP modify
1184 : : !else if (occopt /= 3) then
1185 : : ! write(msg, '(a,i6,a)' )' tphysel /= 0, tsmear == 0, but occopt is not = 3, but ',occopt,'.'
1186 [ # # ]: 0 : else if (occopt /= 3 .and. occopt/=9) then
1187 : 0 : write(msg, '(a,i6,a)' )' tphysel /= 0, tsmear == 0, but occopt is not = 3 or 9, but ',occopt,'.'
1188 : : ! End CP modify
1189 : 0 : ABI_ERROR(msg)
1190 : : end if
1191 : : end if
1192 : :
1193 [ - + ]: 757 : ABI_MALLOC(entder,(-nptsdiv2_def:nptsdiv2_def))
1194 [ - + ]: 757 : ABI_MALLOC(occder,(-nptsdiv2_def:nptsdiv2_def))
1195 [ - + ]: 757 : ABI_MALLOC(smdder,(-nptsdiv2_def:nptsdiv2_def))
1196 [ - + ]: 757 : ABI_MALLOC(workfun,(-nptsdiv2_def:nptsdiv2_def))
1197 [ - + ]: 757 : ABI_MALLOC(work,(-nptsdiv2_def:nptsdiv2_def))
1198 : :
1199 : : ! Prepare the points on the grid
1200 : : ! limit is the value of the argument that will give 0.0 or 1.0 , with
1201 : : ! less than about 1.0d-15 error for 4<=occopt<=8, and less than about 1.0d-12
1202 : : ! error for occopt==3. It is not worth to compute the function beyond
1203 : : ! that point. Even with a less severe requirement, it is significantly
1204 : : ! larger for occopt==3, with an exponential
1205 : : ! tail, than for the other occupation functions, with a Gaussian tail.
1206 : : ! Note that these values are useful in newocc.f also.
1207 : 757 : limit_occ=6.0_dp
1208 : : ! CP modify
1209 : : !if(occopt==3)limit_occ=30.0_dp
1210 [ + + ]: 757 : if(occopt==3 .or. occopt==9)limit_occ=30.0_dp
1211 : : ! End CP modify
1212 [ + + ]: 757 : if(dblsmr /= 0) then
1213 : 10 : tratio = tsmear / tphysel
1214 : 10 : limit_occ=30.0_dp + 6.0_dp*tratio
1215 : : end if
1216 : :
1217 : : ! With nptsdiv2_def=6000 (thus increm=0.001 for 4<=occopt<=8,
1218 : : ! and increm=0.005 for occopt==3, the O(1/N4) algorithm gives 1.0d-12
1219 : : ! accuracy on the stored values occfun and entfun. These, together
1220 : : ! with smdfun and xgrid_prev, need permanently about 0.67 MB, which is affordable.
1221 : 757 : increm=limit_occ/nptsdiv2_def
1222 [ + + ]: 9085514 : do ii=-nptsdiv2_def,nptsdiv2_def
1223 : 9085514 : xgrid_prev(ii)=ii*increm
1224 : : end do
1225 : :
1226 : : ! ---------------------------------------------------------
1227 : : ! Ordinary (unique) smearing function
1228 : : ! ---------------------------------------------------------
1229 [ + + ]: 757 : if (dblsmr == 0) then
1230 : :
1231 : : ! Compute the unnormalized smeared delta function between -limit_occ and +limit_occ
1232 : : ! (well, they are actually normalized ...)
1233 : :
1234 : : ! CP modify
1235 : : !if(occopt==3)then
1236 : : if(occopt==3 .or. occopt==9)then
1237 : : ! End CP modify
1238 : : ! Fermi-Dirac
1239 [ + + ]: 2454818 : do ii=0,nptsdiv2_def
1240 : 2454409 : xx=xgrid_prev(ii)
1241 : 2454409 : smdfun_prev( ii,1)=0.25_dp/(cosh(xx/2.0_dp)**2)
1242 : 2454818 : smdfun_prev(-ii,1)=smdfun_prev(ii,1)
1243 : : end do
1244 : :
1245 : : else if(occopt==4 .or. occopt==5)then
1246 : : ! Cold smearing of Marzari, two values of the "a" parameter being possible
1247 : : ! first value gives minimization of the bump
1248 : : if(occopt==4)aa=-.5634
1249 : : ! second value gives monotonic occupation function
1250 [ + + ]: 122 : if(occopt==5)aa=-.8165
1251 : :
1252 : 122 : dsqrpi=1.0_dp/sqrt(pi)
1253 [ + + ]: 732244 : do ii=0,nptsdiv2_def
1254 : 732122 : xx=xgrid_prev(ii)
1255 : 732122 : gauss=dsqrpi*exp(-xx**2)
1256 : 732122 : smdfun_prev( ii,1)=gauss*(1.5_dp+xx*(-aa*1.5_dp+xx*(-1.0_dp+aa*xx)))
1257 : 732244 : smdfun_prev(-ii,1)=gauss*(1.5_dp+xx*( aa*1.5_dp+xx*(-1.0_dp-aa*xx)))
1258 : : end do
1259 : :
1260 : : else if(occopt==6)then
1261 : :
1262 : : ! First order Hermite-Gaussian of Paxton and Methfessel
1263 : 30010 : dsqrpi=1.0_dp/sqrt(pi)
1264 [ + + ]: 30010 : do ii=0,nptsdiv2_def
1265 : 30005 : xx=xgrid_prev(ii)
1266 : 30005 : smdfun_prev( ii,1)=dsqrpi*(1.5_dp-xx**2)*exp(-xx**2)
1267 : 30010 : smdfun_prev(-ii,1)=smdfun_prev(ii,1)
1268 : : end do
1269 : :
1270 : : else if(occopt==7)then
1271 : :
1272 : : ! Gaussian smearing
1273 : 1266422 : dsqrpi=1.0_dp/sqrt(pi)
1274 [ + + ]: 1266422 : do ii=0,nptsdiv2_def
1275 : 1266211 : xx=xgrid_prev(ii)
1276 : 1266211 : smdfun_prev( ii,1)=dsqrpi*exp(-xx**2)
1277 : 1266422 : smdfun_prev(-ii,1)=smdfun_prev(ii,1)
1278 : : end do
1279 : :
1280 : : else if(occopt==8)then
1281 : :
1282 : : ! Constant value of the delta function over the smearing interval, for testing purposes only.
1283 [ # # ]: 0 : do ii=0,nptsdiv2_def
1284 : 0 : xx=xgrid_prev(ii)
1285 [ # # ]: 0 : if(xx>half+tol8)then
1286 : 0 : smdfun_prev( ii,1)=zero
1287 [ # # ]: 0 : else if(xx<half-tol8)then
1288 : 0 : smdfun_prev( ii,1)=one
1289 : : else
1290 : 0 : smdfun_prev( ii,1)=half
1291 : : end if
1292 : 0 : smdfun_prev(-ii,1)=smdfun_prev(ii,1)
1293 : : end do
1294 : :
1295 : : else
1296 : 0 : ABI_BUG(sjoin('Occopt: ', itoa(occopt),' is not allowed in getnel.'))
1297 : : end if
1298 : :
1299 : : else if (dblsmr /= 0) then
1300 : : ! ---------------------------------------------------------
1301 : : ! smear FD delta with occopt delta calculated in smdfun_prev
1302 : : ! ---------------------------------------------------------
1303 : :
1304 : 10 : nconvd2 = 6000
1305 : 10 : convlim = 10.0_dp
1306 : 10 : incconv = convlim / nconvd2
1307 : :
1308 : : ! store smearing functions in smd1 and smd2
1309 [ - + ]: 10 : ABI_MALLOC(smd1,(-nconvd2:nconvd2))
1310 [ - + ]: 10 : ABI_MALLOC(smd2,(-nconvd2:nconvd2))
1311 [ - + ]: 10 : ABI_MALLOC(tgrid,(-nconvd2:nconvd2))
1312 : :
1313 : : ! FD function in smd1( ii) and second smearing delta in smd2( ii)
1314 : : !
1315 : : ! smd1(:) contains delta_FD ( x )
1316 [ + + ]: 60020 : do ii=0,nconvd2
1317 : 60010 : tgrid(ii)=ii*incconv
1318 : 60010 : tgrid(-ii)=-tgrid(ii)
1319 : 60010 : tt=tgrid(ii)
1320 : 60010 : smd1( ii)=0.25_dp/(cosh(tt/2.0_dp)**2)
1321 : 60020 : smd1(-ii)=smd1(ii)
1322 : : end do
1323 : :
1324 : : ! check input values of occopt and fill smd2(:) with appropriate data:
1325 : : ! smd2(:) contains delta_resmear ( x )
1326 : : ! CP modify
1327 : : !if(occopt == 3) then
1328 : : if(occopt == 3 .or. occopt==9) then
1329 : : ! End CP modify
1330 : : write(msg, '(a,a)' )&
1331 : 0 : 'Occopt=3 is not allowed as a re-smearing.', &
1332 : 0 : 'Use a single FD, or re-smear with a different delta type (faster cutoff). '
1333 : 0 : ABI_ERROR(msg)
1334 : : else if(occopt==4 .or. occopt==5)then
1335 : : ! Cold smearing of Marzari, two values of the "a" parameter being possible
1336 : : ! first value gives minimization of the bump
1337 : : if(occopt==4)aa=-.5634
1338 : : ! second value gives monotonic occupation function
1339 [ - + ]: 10 : if(occopt==5)aa=-.8165
1340 : :
1341 : 10 : dsqrpi=1.0_dp/sqrt(pi)
1342 [ + + ]: 60020 : do ii=0,nconvd2
1343 : 60010 : tt=tgrid(ii)
1344 : 60010 : gauss=dsqrpi*exp(-tt**2)
1345 : 60010 : smd2( ii)=gauss*(1.5_dp+tt*(-aa*1.5_dp+tt*(-1.0_dp+aa*tt)))
1346 : 60020 : smd2(-ii)=gauss*(1.5_dp+tt*( aa*1.5_dp+tt*(-1.0_dp-aa*tt)))
1347 : : end do
1348 : : else if(occopt==6)then
1349 : 0 : dsqrpi=1.0_dp/sqrt(pi)
1350 [ # # ]: 0 : do ii=0,nconvd2
1351 : 0 : tt=tgrid(ii)
1352 : 0 : smd2( ii)=dsqrpi*(1.5_dp-tt**2)*exp(-tt**2)
1353 : 0 : smd2(-ii)=smd2(ii)
1354 : : end do
1355 : : else if(occopt==7)then
1356 : 0 : dsqrpi=1.0_dp/sqrt(pi)
1357 [ # # ]: 0 : do ii=0,nconvd2
1358 : 0 : tt=tgrid(ii)
1359 : 0 : smd2( ii)=dsqrpi*exp(-tt**2)
1360 : 0 : smd2(-ii)=smd2(ii)
1361 : : end do
1362 : : else if(occopt==8)then
1363 [ # # ]: 0 : do ii=0,nconvd2
1364 : 0 : tt=tgrid(ii)
1365 [ # # ]: 0 : if(tt>half+tol8)then
1366 : 0 : smd2( ii)=zero
1367 [ # # ]: 0 : else if(tt<half-tol8)then
1368 : 0 : smd2( ii)=one
1369 : : else
1370 : 0 : smd2( ii)=half
1371 : : end if
1372 : 0 : smd2(-ii)=smd2(ii)
1373 : : end do
1374 : : else
1375 : 0 : ABI_BUG(sjoin('Occopt: ', itoa(occopt),' is not allowed in getnel.'))
1376 : : end if
1377 : :
1378 : : ! Use O(1/N4) algorithm from Num Rec (see below)
1379 : : !
1380 : : ! The grid for the convoluted delta is taken (conservatively)
1381 : : ! to be that for the FD delta ie 6000 pts in [-limit_occ;limit_occ]
1382 : : ! Smearing functions are given on [-dbllim;dbllim] and the grid must
1383 : : ! superpose the normal grid on [-limit_occ:limit_occ]
1384 : : ! The maximal interval for integration of the convolution is
1385 : : ! [-dbllim+limit_occ+lim(delta2);dbllim-limit_occ-lim(delta2)] =
1386 : : ! [-dbllim+36;dbllim-36]
1387 : :
1388 : : ! test the smdFD function for extreme values:
1389 : : ! do jj=-nptsdiv2_def,-nptsdiv2_def
1390 : : ! do ii=-nconvd2+4,nconvd2
1391 : : ! call smdFD(xgrid_prev(jj) - tgrid(ii)*tratio, resFD)
1392 : : ! write(std_out,*) 'ii jj = ', ii,jj, ' smdFD (', xgrid_prev(jj) - tgrid(ii)*tratio, ') ', resFD
1393 : : ! end do
1394 : : ! end do
1395 : :
1396 : 10 : expinc = exp(half*incconv*tratio)
1397 : :
1398 : : ! jj = position of point at which we are calculating smdfun_prev
1399 [ + + ]: 120020 : do jj=-nptsdiv2_def,nptsdiv2_def
1400 : : ! Do not care about the 8 boundary points,
1401 : : ! where the values should be extremely small anyway
1402 : 120010 : smdfun_prev(jj,1)=0.0_dp
1403 : : ! only add contribution with delta_FD > 1.0d-100
1404 [ + + ]: 120010 : nmaxFD = floor (( maxFDarg+xgrid_prev(jj)) / tratio / incconv )
1405 : 120010 : nmaxFD = min (nmaxFD, nconvd2)
1406 [ + + ]: 120010 : nminFD = ceiling((-maxFDarg+xgrid_prev(jj)) / tratio / incconv )
1407 : 120010 : nminFD = max (nminFD, -nconvd2)
1408 : :
1409 : : ! Calculate the Fermi-Dirac distrib at point xgrid_prev(jj)-tgrid(ii)*tratio
1410 : 120010 : expxo2 = exp (-half*(xgrid_prev(jj) - (nminFD)*incconv*tratio))
1411 : 120010 : expx22 = expxo2*expxo2
1412 : 120010 : tmpexpsum = expxo2 / (expx22 + 1.0_dp)
1413 : 120010 : resFD4 = tmpexpsum * tmpexpsum
1414 : 120010 : expxo2 = expxo2*expinc
1415 : 120010 : expx22 = expxo2*expxo2
1416 : 120010 : tmpexpsum = expxo2 / (expx22 + 1.0_dp)
1417 : 120010 : resFD3 = tmpexpsum * tmpexpsum
1418 : 120010 : expxo2 = expxo2*expinc
1419 : 120010 : expx22 = expxo2*expxo2
1420 : 120010 : tmpexpsum = expxo2 / (expx22 + 1.0_dp)
1421 : 120010 : resFD2 = tmpexpsum * tmpexpsum
1422 : 120010 : expxo2 = expxo2*expinc
1423 : 120010 : expx22 = expxo2*expxo2
1424 : 120010 : tmpexpsum = expxo2 / (expx22 + 1.0_dp)
1425 : 120010 : resFD1 = tmpexpsum * tmpexpsum
1426 : :
1427 : : ! core contribution to the integral with constant weight (48)
1428 : 120010 : tmpsmdfun = 0.0_dp
1429 [ + + ]: 1339688356 : do ii=nminFD+4,nmaxFD-4
1430 : 1339568346 : expxo2 = expxo2*expinc
1431 : : ! tmpexpsum = 1.0_dp / (expxo2 + 1.0_dp / expxo2 )
1432 : 1339568346 : expx22 = expxo2*expxo2
1433 : 1339568346 : tmpexpsum = expxo2 / (expx22 + 1.0_dp)
1434 : 1339688356 : tmpsmdfun = tmpsmdfun + smd2(ii) * tmpexpsum * tmpexpsum
1435 : : end do
1436 : :
1437 : : ! Add on end contributions for show (both functions smd and smdFD are very small
1438 : : smdfun_prev(jj,1)=smdfun_prev(jj,1) +48.0_dp*tmpsmdfun &
1439 : : + 31.0_dp*smd2(nminFD+3)*resFD1 -11.0_dp*smd2(nminFD+2)*resFD2 &
1440 : 120010 : + 5.0_dp*smd2(nminFD+1)*resFD3 - smd2(nminFD)*resFD4
1441 : :
1442 : 120010 : expxo2 = expxo2*expinc
1443 : 120010 : expx22 = expxo2*expxo2
1444 : 120010 : tmpexpsum = expxo2 / (expx22 + 1.0_dp)
1445 : 120010 : resFD1 = tmpexpsum * tmpexpsum
1446 : 120010 : expxo2 = expxo2*expinc
1447 : 120010 : expx22 = expxo2*expxo2
1448 : 120010 : tmpexpsum = expxo2 / (expx22 + 1.0_dp)
1449 : 120010 : resFD2 = tmpexpsum * tmpexpsum
1450 : 120010 : expxo2 = expxo2*expinc
1451 : 120010 : expx22 = expxo2*expxo2
1452 : 120010 : tmpexpsum = expxo2 / (expx22 + 1.0_dp)
1453 : 120010 : resFD3 = tmpexpsum * tmpexpsum
1454 : 120010 : expxo2 = expxo2*expinc
1455 : 120010 : expx22 = expxo2*expxo2
1456 : 120010 : tmpexpsum = expxo2 / (expx22 + 1.0_dp)
1457 : 120010 : resFD4 = tmpexpsum * tmpexpsum
1458 : :
1459 : : ! Contribution above
1460 : : smdfun_prev(jj,1)=smdfun_prev(jj,1) &
1461 : : + 31.0_dp*smd2(nmaxFD-3)*resFD1 -11.0_dp*smd2(nmaxFD-2)*resFD2 &
1462 : 120010 : + 5.0_dp*smd2(nmaxFD-1)*resFD3 - smd2(nmaxFD)*resFD4
1463 : 120020 : smdfun_prev(jj,1)=incconv*smdfun_prev(jj,1)/48.0_dp
1464 : : end do
1465 : :
1466 : 10 : secmom = 0.0_dp
1467 : 10 : thdmom = 0.0_dp
1468 : 10 : resmom4 = xgrid_prev(-nptsdiv2_def )*xgrid_prev(-nptsdiv2_def )*smdfun_prev(-nptsdiv2_def , 1)
1469 : 10 : resmom3 = xgrid_prev(-nptsdiv2_def+1)*xgrid_prev(-nptsdiv2_def+1)*smdfun_prev(-nptsdiv2_def+1, 1)
1470 : 10 : resmom2 = xgrid_prev(-nptsdiv2_def+2)*xgrid_prev(-nptsdiv2_def+2)*smdfun_prev(-nptsdiv2_def+2, 1)
1471 : 10 : resmom1 = xgrid_prev(-nptsdiv2_def+3)*xgrid_prev(-nptsdiv2_def+3)*smdfun_prev(-nptsdiv2_def+3, 1)
1472 : 10 : resmom = xgrid_prev(-nptsdiv2_def+4)*xgrid_prev(-nptsdiv2_def+4)*smdfun_prev(-nptsdiv2_def+4, 1)
1473 : 10 : do ii=-nptsdiv2_def+4,nptsdiv2_def-1
1474 : : secmom = secmom + &
1475 : : & ( 17.0_dp*xgrid_prev(ii) *xgrid_prev(ii) *smdfun_prev(ii, 1) &
1476 : : & +42.0_dp*xgrid_prev(ii-1)*xgrid_prev(ii-1)*smdfun_prev(ii-1,1) &
1477 : : & -16.0_dp*xgrid_prev(ii-2)*xgrid_prev(ii-2)*smdfun_prev(ii-2,1) &
1478 : : & + 6.0_dp*xgrid_prev(ii-3)*xgrid_prev(ii-3)*smdfun_prev(ii-3,1) &
1479 : : & - xgrid_prev(ii-4)*xgrid_prev(ii-4)*smdfun_prev(ii-4,1) )
1480 : : resmom4 = resmom3
1481 : : resmom3 = resmom2
1482 : : resmom2 = resmom1
1483 : : resmom1 = resmom
1484 : : resmom = xgrid_prev(ii+1) *xgrid_prev(ii+1) *smdfun_prev(ii+1, 1)
1485 : : end do
1486 : 10 : secmom=increm * secmom / 48.0_dp
1487 : : ! thdmom=increm * thdmom / 48.0_dp
1488 : : !
1489 : : ! smom1 = second moment of delta in smd1(:)
1490 : : ! smom2 = second moment of delta in smd2(:)
1491 : : !
1492 : : smom1 = 0.0_dp
1493 : : smom2 = 0.0_dp
1494 : : tmom1 = 0.0_dp
1495 : : tmom2 = 0.0_dp
1496 : : do ii=-nconvd2+4,nconvd2
1497 : : smom1 = smom1+ &
1498 : : & ( 17.0_dp*tgrid(ii) *tgrid(ii) *smd1(ii) &
1499 : : & +42.0_dp*tgrid(ii-1)*tgrid(ii-1)*smd1(ii-1) &
1500 : : & -16.0_dp*tgrid(ii-2)*tgrid(ii-2)*smd1(ii-2) &
1501 : : & + 6.0_dp*tgrid(ii-3)*tgrid(ii-3)*smd1(ii-3) &
1502 : : & - tgrid(ii-4)*tgrid(ii-4)*smd1(ii-4) )
1503 : : smom2 = smom2+ &
1504 : : & ( 17.0_dp*tgrid(ii) *tgrid(ii) *smd2(ii ) &
1505 : : & +42.0_dp*tgrid(ii-1)*tgrid(ii-1)*smd2(ii-1) &
1506 : : & -16.0_dp*tgrid(ii-2)*tgrid(ii-2)*smd2(ii-2) &
1507 : : & + 6.0_dp*tgrid(ii-3)*tgrid(ii-3)*smd2(ii-3) &
1508 : : & - tgrid(ii-4)*tgrid(ii-4)*smd2(ii-4) )
1509 : : end do
1510 : 10 : smom1 =incconv * smom1 / 48.0_dp
1511 : 10 : smom2 =incconv * smom2 / 48.0_dp
1512 : : ! tmom1 =incconv * tmom1 / 48.0_dp
1513 : : ! tmom2 =incconv * tmom2 / 48.0_dp
1514 : :
1515 : 10 : encorr = smom2*tratio*tratio/secmom
1516 : :
1517 [ - + ]: 10 : ABI_FREE(tgrid)
1518 [ - + ]: 10 : ABI_FREE(smd1)
1519 [ - + ]: 10 : ABI_FREE(smd2)
1520 : :
1521 : : end if
1522 : :
1523 : : ! --------------------------------------------------------
1524 : : ! end of smearing function initialisation, dblsmr case
1525 : : ! --------------------------------------------------------
1526 : :
1527 : :
1528 : : ! Now that the smeared delta function has been initialized, compute the
1529 : : ! occupation function
1530 : 757 : occfun_prev(-nptsdiv2_def,1)=zero
1531 : 757 : entfun_prev(-nptsdiv2_def,1)=zero
1532 : :
1533 : : ! Different algorithms are possible, corresponding to the formulas
1534 : : ! (4.1.11), (4.1.12) and (4.1.14) in Numerical recipes (pp 107 and 108),
1535 : : ! with respective O(1/N2), O(1/N3), O(1/N4) convergence, where N is the
1536 : : ! number of points in the interval.
1537 : 757 : algo=4
1538 : :
1539 : : if(algo==2)then
1540 : :
1541 : : ! Extended trapezoidal rule (4.1.11), taken in a cumulative way
1542 : : do ii=-nptsdiv2_def+1,nptsdiv2_def
1543 : : occfun_prev(ii,1)=occfun_prev(ii-1,1)+increm*(smdfun_prev(ii,1)+smdfun_prev(ii-1,1))/2.0_dp
1544 : : entfun_prev(ii,1)=entfun_prev(ii-1,1)+increm*&
1545 : : & ( -xgrid_prev(ii)*smdfun_prev(ii,1) -xgrid_prev(ii-1)*smdfun_prev(ii-1,1) )/2.0_dp
1546 : : end do
1547 : :
1548 : : else if(algo==3)then
1549 : :
1550 : : ! Derived from (4.1.12). Converges as O(1/N3).
1551 : : ! Do not care about the following points,
1552 : : ! where the values are extremely small anyway
1553 : : occfun_prev(-nptsdiv2_def+1,1)=0.0_dp ; entfun_prev(-nptsdiv2_def+1,1)=0.0_dp
1554 : : do ii=-nptsdiv2_def+2,nptsdiv2_def
1555 : : occfun_prev(ii,1)=occfun_prev(ii-1,1)+increm*&
1556 : : & ( 5.0_dp*smdfun_prev(ii,1) + 8.0_dp*smdfun_prev(ii-1,1) - smdfun_prev(ii-2,1) )/12.0_dp
1557 : : entfun_prev(ii,1)=entfun_prev(ii-1,1)+increm*&
1558 : : & ( 5.0_dp*(-xgrid_prev(ii) )*smdfun_prev(ii,1) &
1559 : : & +8.0_dp*(-xgrid_prev(ii-1))*smdfun_prev(ii-1,1)&
1560 : : & - (-xgrid_prev(ii-2))*smdfun_prev(ii-2,1) )/12.0_dp
1561 : : end do
1562 : :
1563 : : else if(algo==4)then
1564 : :
1565 : : ! Derived from (4.1.14)- alternative extended Simpsons rule. Converges as O(1/N4).
1566 : : ! Do not care about the following points,
1567 : : ! where the values are extremely small anyway
1568 : 757 : occfun_prev(-nptsdiv2_def+1,1)=0.0_dp ; entfun_prev(-nptsdiv2_def+1,1)=0.0_dp
1569 : 757 : occfun_prev(-nptsdiv2_def+2,1)=0.0_dp ; entfun_prev(-nptsdiv2_def+2,1)=0.0_dp
1570 : 757 : occfun_prev(-nptsdiv2_def+3,1)=0.0_dp ; entfun_prev(-nptsdiv2_def+3,1)=0.0_dp
1571 [ + + ]: 9082486 : do ii=-nptsdiv2_def+4,nptsdiv2_def
1572 : : occfun_prev(ii,1)=occfun_prev(ii-1,1)+increm*&
1573 : : & ( 17.0_dp*smdfun_prev(ii,1) &
1574 : : & +42.0_dp*smdfun_prev(ii-1,1)&
1575 : : & -16.0_dp*smdfun_prev(ii-2,1)&
1576 : : & + 6.0_dp*smdfun_prev(ii-3,1)&
1577 : 9081729 : & - smdfun_prev(ii-4,1) )/48.0_dp
1578 : : entfun_prev(ii,1)=entfun_prev(ii-1,1)+increm*&
1579 : : & ( 17.0_dp*(-xgrid_prev(ii) )*smdfun_prev(ii,1) &
1580 : : & +42.0_dp*(-xgrid_prev(ii-1))*smdfun_prev(ii-1,1)&
1581 : : & -16.0_dp*(-xgrid_prev(ii-2))*smdfun_prev(ii-2,1)&
1582 : : & + 6.0_dp*(-xgrid_prev(ii-3))*smdfun_prev(ii-3,1)&
1583 : 9082486 : & - (-xgrid_prev(ii-4))*smdfun_prev(ii-4,1) )/48.0_dp
1584 : : end do
1585 : :
1586 : : end if ! End of choice between different algorithms for integration
1587 : :
1588 : : ! Normalize the functions (actually not needed for occopt=3..7)
1589 : 757 : factor=1.0_dp/occfun_prev(nptsdiv2_def,1)
1590 [ + + ]: 9085514 : smdfun_prev(:,1)=smdfun_prev(:,1)*factor
1591 [ + + ]: 9085514 : occfun_prev(:,1)=occfun_prev(:,1)*factor
1592 [ + + ]: 9085514 : entfun_prev(:,1)=entfun_prev(:,1)*factor
1593 : :
1594 : : ! Compute the cubic spline fitting of the smeared delta function
1595 : 757 : yp1=0.0_dp ; ypn=0.0_dp
1596 [ + + ]: 9085514 : workfun(:)=smdfun_prev(:,1)
1597 : 757 : call spline(xgrid_prev, workfun, (2*nptsdiv2_def+1), yp1, ypn, smdder)
1598 [ + + ]: 9085514 : smdfun_prev(:,2)=smdder(:)
1599 : :
1600 : : ! Compute the cubic spline fitting of the occupation function
1601 : : yp1=0.0_dp ; ypn=0.0_dp
1602 [ + + ]: 9085514 : workfun(:)=occfun_prev(:,1)
1603 : 757 : call spline(xgrid_prev, workfun, (2*nptsdiv2_def+1), yp1, ypn, occder)
1604 [ + + ]: 9085514 : occfun_prev(:,2)=occder(:)
1605 : :
1606 : : ! Compute the cubic spline fitting of the entropy function
1607 : : yp1=0.0_dp ; ypn=0.0_dp
1608 [ + + ]: 9085514 : workfun(:)=entfun_prev(:,1)
1609 : 757 : call spline(xgrid_prev, workfun, (2*nptsdiv2_def+1), yp1, ypn, entder)
1610 [ + + ]: 9085514 : entfun_prev(:,2)=entder(:)
1611 : :
1612 [ - + ]: 757 : ABI_FREE(entder)
1613 [ - + ]: 757 : ABI_FREE(occder)
1614 [ - + ]: 757 : ABI_FREE(smdder)
1615 [ - + ]: 757 : ABI_FREE(work)
1616 [ - + ]: 757 : ABI_FREE(workfun)
1617 : :
1618 : : end if
1619 : :
1620 [ + + ]: 552699 : if (abs(tphysel)<tol12) then
1621 [ + + ]: 550402 : if (tsmear == zero) then
1622 : 325 : tsmearinv = huge_tsmearinv
1623 : : else
1624 : 550077 : tsmearinv=one/tsmear
1625 : : end if
1626 : : else
1627 : 2297 : tsmearinv=one/tphysel
1628 : : end if
1629 : :
1630 [ + + + + ]: 13267539495 : entfun(:,:) = entfun_prev(:,:)
1631 [ + + + + ]: 13267539495 : occfun(:,:) = occfun_prev(:,:)
1632 [ + + + + ]: 13267539495 : smdfun(:,:) = smdfun_prev(:,:)
1633 [ + + ]: 6633493398 : xgrid(:) = xgrid_prev(:)
1634 : 552699 : limit = limit_occ
1635 : 552699 : nptsdiv2 = nptsdiv2_def
1636 : :
1637 : : end subroutine init_occ_ent
1638 : : !!***
1639 : :
1640 : : !!****f* m_occ/occeig
1641 : : !! NAME
1642 : : !! occeig
1643 : : !!
1644 : : !! FUNCTION
1645 : : !! For each pair of active bands (m,n), generates ratios
1646 : : !! that depend on the difference between occupation numbers and eigenvalues.
1647 : : !!
1648 : : !! INPUTS
1649 : : !! doccde_k(nband_k)=derivative of occ_k wrt the energy
1650 : : !! doccde_kq(nband_k)=derivative of occ_kq wrt the energy
1651 : : !! eig0_k(nband_k)=GS eigenvalues at k
1652 : : !! eig0_kq(nband_k)=GS eigenvalues at k+q
1653 : : !! nband_k=number of bands
1654 : : !! occopt=option for occupancies
1655 : : !! occ_k(nband_k)=occupation number for each band at k
1656 : : !! occ_kq(nband_k)=occupation number for each band at k+q
1657 : : !!
1658 : : !! OUTPUT
1659 : : !! rocceig(nband_k,nband_k)$= (occ_{k,q}(m)-occ_k(n))/(eig0_{k,q}(m)-eig0_k(n))$,
1660 : : !! if this ratio has been attributed to the band n, 0.0_dp otherwise
1661 : : !!
1662 : : !! NOTES
1663 : : !! Supposing the occupations numbers differ:
1664 : : !! if $abs(occ_{k,q}(m)) < abs(occ_k(n))$
1665 : : !! $rocceig(m,n)=(occ_{k,q}(m)-occ_k(n))/(eig0_{k,q}(m)-eig0_k(n)) $
1666 : : !! if $abs(occ_{k,q}(m))>abs(occ_k(n))$
1667 : : !! rocceig(m,n)=0.0_dp
1668 : : !!
1669 : : !! If the occupation numbers are close enough, then
1670 : : !! if the eigenvalues are also close, take the derivative
1671 : : !! $ rocceig(m,n)=\frac{1}{2}*docc/deig0 $
1672 : : !! otherwise,
1673 : : !! $ rocceig(m,n)=\frac{1}{2}*(occ_{k,q}(m)-occ_k(n))/(eig0_{k,q}(m)-eig0_k(n))$
1674 : : !!
1675 : : !! SOURCE
1676 : :
1677 : 855845 : subroutine occeig(doccde_k,doccde_kq,eig0_k,eig0_kq,nband_k,occopt,occ_k,occ_kq,rocceig)
1678 : :
1679 : : !Arguments ------------------------------------
1680 : : !scalars
1681 : : integer,intent(in) :: nband_k,occopt
1682 : : !arrays
1683 : : real(dp),intent(in) :: doccde_k(nband_k),doccde_kq(nband_k),eig0_k(nband_k)
1684 : : real(dp),intent(in) :: eig0_kq(nband_k),occ_k(nband_k),occ_kq(nband_k)
1685 : : real(dp),intent(out) :: rocceig(nband_k,nband_k)
1686 : :
1687 : : !Local variables-------------------------------
1688 : : !scalars
1689 : : integer :: ibandk,ibandkq
1690 : : real(dp) :: diffabsocc,diffeig,diffocc,ratio,sumabsocc
1691 : : character(len=500) :: msg
1692 : : ! *************************************************************************
1693 : :
1694 : : ! The parameter tol5 defines the treshhold for degeneracy, and the width of the step function
1695 : :
1696 [ + + + + ]: 41411569 : rocceig(:,:) = zero
1697 : :
1698 [ + + ]: 5406999 : do ibandk=1,nband_k
1699 [ + + ]: 41411569 : do ibandkq=1,nband_k
1700 : 36004570 : diffeig=eig0_kq(ibandkq)-eig0_k(ibandk)
1701 : 36004570 : diffocc=occ_kq(ibandkq)-occ_k(ibandk)
1702 : :
1703 [ + + ]: 36004570 : if( abs(diffeig) > tol5 ) then
1704 : 32883542 : ratio=diffocc/diffeig
1705 : : else
1706 [ + + ]: 3121028 : if(occopt<3)then
1707 : : ! In a non-metallic case, if the eigenvalues are degenerate,
1708 : : ! the occupation numbers must also be degenerate, in which
1709 : : ! case there is no contribution from this pair of bands
1710 [ - + ]: 2557804 : if( abs(diffocc) > tol5 ) then
1711 : : write(msg,'(a,a,a,a,a,a,a,2(a,i4,a,es16.6,a,es16.6,a,a),a)' ) &
1712 : 0 : 'In a non-metallic case (occopt<3), for a RF calculation,',ch10,&
1713 : 0 : 'if the eigenvalues are degenerate,',' the occupation numbers must also be degenerate.',ch10,&
1714 : 0 : 'However, the following pair of states gave :',ch10,&
1715 : 0 : 'k -state, band number',ibandk,', occ=',occ_k(ibandk),'eigenvalue=',eig0_k(ibandk),',',ch10,&
1716 : 0 : ' kq-state, band number',ibandkq,', occ=',occ_kq(ibandkq),', eigenvalue=',eig0_kq(ibandkq),'.',ch10,&
1717 : 0 : 'Action: change occopt, consistently, in GS and RF calculations.'
1718 : 0 : ABI_ERROR(msg)
1719 : : end if
1720 : : ratio=0.0_dp
1721 : : else
1722 : : ! In the metallic case, one can compute a better approximation of the
1723 : : ! ratio by using derivatives doccde
1724 : 563224 : ratio=0.5_dp*(doccde_kq(ibandkq)+doccde_k(ibandk))
1725 : : ! write(std_out,*)' occeig : ibandkq,doccde_kq(ibandkq)',ibandkq,doccde_kq(ibandkq)
1726 : : ! write(std_out,*)' ibandk ,doccde_k (ibandk )',ibandk,doccde_k(ibandk)
1727 : : end if
1728 : : end if
1729 : :
1730 : : ! Here, must pay attention to the smallness of some coefficient
1731 : 36004570 : diffabsocc=abs(occ_k(ibandk))-abs(occ_kq(ibandkq))
1732 : 36004570 : sumabsocc=abs(occ_k(ibandk))+abs(occ_kq(ibandkq))
1733 [ + + ]: 40555724 : if(sumabsocc>tol8)then
1734 [ + + ]: 26269450 : if( diffabsocc > sumabsocc*tol5 ) then
1735 : 4140219 : rocceig(ibandkq,ibandk)=ratio
1736 [ + + ]: 22129231 : else if ( diffabsocc >= -sumabsocc*tol5 ) then
1737 : 17988454 : rocceig(ibandkq,ibandk)=0.5_dp*ratio
1738 : : else
1739 : 4140777 : rocceig(ibandkq,ibandk)=0.0_dp
1740 : : end if
1741 : : end if
1742 : :
1743 : : end do ! ibandkq
1744 : : end do ! ibandk
1745 : :
1746 : 855845 : end subroutine occeig
1747 : : !!***
1748 : :
1749 : : !----------------------------------------------------------------------
1750 : :
1751 : : !!****f* m_occ/occ_fd
1752 : : !! NAME
1753 : : !! occ_fd
1754 : : !!
1755 : : !! FUNCTION
1756 : : !! Fermi-Dirac statistics: 1 / [(exp((e - mu)/ KT) + 1]
1757 : : !! Note that occ_fs in [0, 1] so the spin factor is not included, unlike the
1758 : : !! occupations stored in ebands%occ.
1759 : : !!
1760 : : !! INPUTS
1761 : : !! ee=Single particle energy in Ha
1762 : : !! kT=Value of K_Boltzmann x T in Ha.
1763 : : !! mu=Chemical potential in Ha.
1764 : : !!
1765 : : !! SOURCE
1766 : :
1767 : 41421098 : elemental real(dp) function occ_fd(ee, kT, mu)
1768 : :
1769 : : !Arguments ------------------------------------
1770 : : real(dp),intent(in) :: ee, kT, mu
1771 : :
1772 : : !Local variables ------------------------------
1773 : : real(dp) :: ee_mu,arg
1774 : : ! *************************************************************************
1775 : :
1776 : 41421098 : ee_mu = ee - mu
1777 : :
1778 : : ! 1 kelvin [K] = 3.16680853419133E-06 Hartree
1779 [ + + ]: 41421098 : if (kT > tol6) then
1780 : 23910079 : arg = ee_mu / kT
1781 [ + + ]: 23910079 : if (arg > maxFDarg) then
1782 : : occ_fd = zero
1783 [ + + ]: 15098260 : else if (arg < -maxFDarg) then
1784 : : occ_fd = one
1785 : : else
1786 : 6984950 : occ_fd = one / (exp(arg) + one)
1787 : : end if
1788 : : else
1789 : : ! Heaviside
1790 [ + + ]: 17511019 : if (ee_mu > zero) then
1791 : : occ_fd = zero
1792 [ + + ]: 8791513 : else if (ee_mu < zero) then
1793 : : occ_fd = one
1794 : : else
1795 : 41421098 : occ_fd = half
1796 : : end if
1797 : : end if
1798 : :
1799 : 41421098 : end function occ_fd
1800 : : !!***
1801 : :
1802 : : !----------------------------------------------------------------------
1803 : :
1804 : : !!****f* m_occ/occ_dfde
1805 : : !! NAME
1806 : : !! occ_dfde
1807 : : !!
1808 : : !! FUNCTION
1809 : : !! Derivative of Fermi-Dirac statistics: - (exp((e - mu)/ KT) / KT[(exp((e - mu)/ KT) + 1]^2
1810 : : !! Note that kT is given in Hartree so the derivative as well
1811 : : !!
1812 : : !! INPUTS
1813 : : !! ee=Single particle energy in Ha
1814 : : !! kT=Value of K_Boltzmann x T in Ha.
1815 : : !! mu=Chemical potential in Ha.
1816 : : !!
1817 : : !! SOURCE
1818 : :
1819 : 204693600 : elemental real(dp) function occ_dfde(ee, kT, mu)
1820 : :
1821 : : !Arguments ------------------------------------
1822 : : real(dp),intent(in) :: ee, kT, mu
1823 : :
1824 : : !Local variables ------------------------------
1825 : : real(dp) :: ee_mu,arg
1826 : : ! *************************************************************************
1827 : :
1828 : 204693600 : ee_mu = ee - mu
1829 : :
1830 : : ! 1 kelvin [K] = 3.16680853419133E-06 Hartree
1831 [ + + ]: 204693600 : if (kT > tol6) then
1832 : 102736146 : arg = ee_mu / kT
1833 [ + + ]: 102736146 : if (arg > maxDFDarg) then
1834 : : occ_dfde = zero
1835 [ + + ]: 57105075 : else if (arg < -maxDFDarg) then
1836 : : occ_dfde = zero
1837 : : else
1838 : 11474029 : occ_dfde = - exp(arg) / (exp(arg) + one)**2 / kT
1839 : : end if
1840 : : else
1841 : : occ_dfde = zero
1842 : : end if
1843 : :
1844 : 204693600 : end function occ_dfde
1845 : : !!***
1846 : :
1847 : : !----------------------------------------------------------------------
1848 : :
1849 : : !!****f* m_occ/occ_be
1850 : : !! NAME
1851 : : !! occ_be
1852 : : !!
1853 : : !! FUNCTION
1854 : : !! Bose-Einstein statistics 1 / [(exp((e - mu)/ KT) - 1]
1855 : : !!
1856 : : !! INPUTS
1857 : : !! ee=Single particle energy in Ha
1858 : : !! kT=Value of K_Boltzmann x T in Ha.
1859 : : !! mu=Chemical potential in Ha (usually zero)
1860 : : !!
1861 : : !! SOURCE
1862 : :
1863 : 4470992 : elemental real(dp) function occ_be(ee, kT, mu)
1864 : :
1865 : : !Arguments ------------------------------------
1866 : : real(dp),intent(in) :: ee, kT, mu
1867 : :
1868 : : !Local variables ------------------------------
1869 : : real(dp) :: ee_mu, arg
1870 : : ! *************************************************************************
1871 : :
1872 : 4470992 : ee_mu = ee - mu
1873 : :
1874 : : ! 1 kelvin [K] = 3.16680853419133E-06 Hartree
1875 [ + + ]: 4470992 : if (kT > tol12) then
1876 : 4193278 : arg = ee_mu / kT
1877 [ + - + - ]: 4193278 : if (arg > tol12 .and. arg < maxBEarg) then
1878 : 4193278 : occ_be = one / (exp(arg) - one)
1879 : : else
1880 : : occ_be = zero
1881 : : end if
1882 : : else
1883 : : ! No condensate for T --> 0
1884 : : occ_be = zero
1885 : : end if
1886 : :
1887 : 4470992 : end function occ_be
1888 : : !!***
1889 : :
1890 : : !----------------------------------------------------------------------
1891 : :
1892 : : !!****f* m_occ/occ_dbe
1893 : : !! NAME
1894 : : !! occ_dbe
1895 : : !!
1896 : : !! FUNCTION
1897 : : !! Derivative of Bose-Einstein statistics (exp((e - mu)/ KT) / KT[(exp((e - mu)/ KT) - 1]^2
1898 : : !! Note that kT is given in Hartree so the derivative as well
1899 : : !!
1900 : : !! INPUTS
1901 : : !! ee=Single particle energy in Ha
1902 : : !! kT=Value of K_Boltzmann x T in Ha.
1903 : : !! mu=Chemical potential in Ha (usually zero)
1904 : : !!
1905 : : !! SOURCE
1906 : :
1907 : 0 : elemental real(dp) function occ_dbe(ee, kT, mu)
1908 : :
1909 : : !Arguments ------------------------------------
1910 : : real(dp),intent(in) :: ee, kT, mu
1911 : :
1912 : : !Local variables ------------------------------
1913 : : real(dp) :: ee_mu, arg
1914 : : ! *************************************************************************
1915 : :
1916 : 0 : ee_mu = ee - mu
1917 : :
1918 : : ! 1 kelvin [K] = 3.16680853419133E-06 Hartree
1919 [ # # ]: 0 : if (kT > tol12) then
1920 : 0 : arg = ee_mu / kT
1921 [ # # # # ]: 0 : if (arg > tol12 .and. arg < maxDBEarg) then
1922 : 0 : occ_dbe = exp(arg) / (kT * (exp(arg) - one)**2)
1923 : : else
1924 : : occ_dbe = zero
1925 : : end if
1926 : : else
1927 : : ! No condensate for T --> 0
1928 : : occ_dbe = zero
1929 : : end if
1930 : :
1931 : 0 : end function occ_dbe
1932 : : !!***
1933 : :
1934 : : !----------------------------------------------------------------------
1935 : :
1936 : : !!****f* m_occ/dos_hdr_write
1937 : : !!
1938 : : !! NAME
1939 : : !! dos_hdr_write
1940 : : !!
1941 : : !! FUNCTION
1942 : : !! Write the header of the DOS files, for both smearing and tetrahedron methods.
1943 : : !!
1944 : : !! INPUTS
1945 : : !! deltaene=increment of DOS energy arguments
1946 : : !! enemax=maximal value of the DOS energy argument
1947 : : !! enemin=minimal value of the DOS energy argument
1948 : : !! nene=number of DOS energy argument
1949 : : !! eigen(mband*nkpt*nsppol)=eigenvalues (input or init to large number), hartree
1950 : : !! fermie=fermi energy useful for band alignment...
1951 : : !! fermih= fermi energy of thermalized excited holes when occopt = 9
1952 : : !! mband=maximum number of bands
1953 : : !! nband(nkpt*nsppol)=number of bands at each k point
1954 : : !! nkpt=number of k points
1955 : : !! nsppol=1 for unpolarized, 2 for spin-polarized
1956 : : !! occopt=option for occupancies, or re-smearing scheme if dblsmr /= 0
1957 : : !! prtdos=1 for smearing technique, 2 or 3 for tetrahedron technique
1958 : : !! tphysel="physical" electronic temperature with FD occupations
1959 : : !! tsmear=smearing width (or temperature)
1960 : : !! unitdos=unit number of output of the DOS.
1961 : : !!
1962 : : !! OUTPUT
1963 : : !! Only writing.
1964 : : !!
1965 : : !! SOURCE
1966 : :
1967 : 76 : subroutine dos_hdr_write(deltaene,eigen,enemax,enemin,fermie,fermih,mband,nband,nene,&
1968 : : nkpt,nsppol,occopt,prtdos,tphysel,tsmear,unitdos)
1969 : :
1970 : : !Arguments ------------------------------------
1971 : : !scalars
1972 : : integer,intent(in) :: mband,nkpt,nsppol,occopt,prtdos,unitdos,nene
1973 : : real(dp),intent(in) :: fermie,fermih,tphysel,tsmear
1974 : : real(dp),intent(in) :: deltaene,enemax,enemin
1975 : : !arrays
1976 : : integer,intent(in) :: nband(nkpt*nsppol)
1977 : : real(dp),intent(in) :: eigen(mband*nkpt*nsppol)
1978 : :
1979 : : !Local variables-------------------------------
1980 : : character(len=500) :: msg
1981 : : ! *************************************************************************
1982 : :
1983 : : ! Write the DOS file
1984 : 76 : write(msg, '(7a,i2,a,i5,a,i4)' ) "#",ch10, &
1985 : 76 : '# ABINIT package : DOS file ',ch10,"#",ch10,&
1986 : 152 : '# nsppol =',nsppol,', nkpt =',nkpt,', nband(1)=',nband(1)
1987 : 76 : call wrtout(unitdos, msg)
1988 : :
1989 [ + + ]: 76 : if (any(prtdos== [1,4])) then
1990 : : write(msg, '(a,i2,a,f6.3,a,f6.3,a)' ) &
1991 : 9 : '# Smearing technique, occopt =',occopt,', tsmear=',tsmear,' Hartree, tphysel=',tphysel,' Hartree'
1992 : : else
1993 : 67 : write(msg, '(a)' ) '# Tetrahedron method '
1994 : : end if
1995 : 76 : call wrtout(unitdos, msg)
1996 : :
1997 [ + - ]: 76 : if (mband*nkpt*nsppol>=3) then
1998 : 76 : write(msg, '(a,3f8.3,2a)' )'# For identification : eigen(1:3)=',eigen(1:3),ch10,"#"
1999 : : else
2000 : 0 : write(msg, '(a,3f8.3)' ) '# For identification : eigen=',eigen
2001 : 0 : write(msg, '(3a)')trim(msg),ch10,"#"
2002 : : end if
2003 : 76 : call wrtout(unitdos, msg)
2004 : :
2005 [ - + ]: 76 : if (occopt == 9) then
2006 : 0 : write(msg, '(a,f16.8, f16.8)' ) '# Fermi energy for electrons and holes ', fermie, fermih
2007 : : else
2008 : 76 : write(msg, '(a,f16.8)' ) '# Fermi energy : ', fermie
2009 : : end if
2010 : 76 : call wrtout(unitdos, msg)
2011 : :
2012 [ + + ]: 76 : if (prtdos==1) then
2013 : 9 : write(msg, '(5a)' ) "#",ch10,&
2014 : 9 : '# The DOS (in electrons/Hartree/cell) and integrated DOS (in electrons/cell),',&
2015 : 18 : ch10,'# as well as the DOS with tsmear halved and doubled, are computed,'
2016 : :
2017 [ + + ]: 67 : else if (prtdos==2)then
2018 : 4 : write(msg, '(3a)' ) "#",ch10,&
2019 : 8 : '# The DOS (in electrons/Hartree/cell) and integrated DOS (in electrons/cell) are computed,'
2020 : :
2021 [ + + ]: 63 : else if (any(prtdos == [3, 4])) then
2022 : 62 : write(msg, '(5a)' ) "#",ch10,&
2023 : 62 : '# The local DOS (in electrons/Hartree for one atomic sphere)',ch10,&
2024 : 124 : '# and integrated local DOS (in electrons for one atomic sphere) are computed.'
2025 : :
2026 [ + - ]: 1 : else if (prtdos==5)then
2027 : 1 : write(msg, '(9a)' ) "#",ch10,&
2028 : 1 : '# The spin component DOS (in electrons/Hartree/cell)',ch10,&
2029 : 1 : '# and integrated spin component DOS (in electrons/cell) are computed.',ch10,&
2030 : 1 : '# Remember that the wf are eigenstates of S_z and S^2, not S_x and S_y',ch10,&
2031 : 2 : '# so the latter will not always sum to 0 for paired electronic states.'
2032 : : end if
2033 : 76 : call wrtout(unitdos, msg)
2034 : :
2035 : : write(msg, '(a,i5,a,a,a,f9.4,a,f9.4,a,f8.5,a,a,a)' )&
2036 : 76 : '# at ',nene,' energies (in Hartree) covering the interval ',ch10,&
2037 : 152 : '# between ',enemin,' and ',enemax,' Hartree by steps of ',deltaene,' Hartree.',ch10,"#"
2038 : 76 : call wrtout(unitdos, msg)
2039 : :
2040 [ + + ]: 76 : if (prtdos==1) then
2041 : : write(msg, '(a,a)' )&
2042 : 9 : '# energy DOS Integr. DOS ',' DOS DOS '
2043 : 9 : call wrtout(unitdos,msg)
2044 : :
2045 : : write(msg, '(a)' )&
2046 : 9 : '# (tsmear/2) (tsmear*2) '
2047 : 9 : call wrtout(unitdos,msg)
2048 : : else
2049 : 67 : write(msg, '(a)' ) '# energy DOS '
2050 : : end if
2051 : :
2052 : 76 : end subroutine dos_hdr_write
2053 : : !!***
2054 : :
2055 : : !!****f* m_occ/get_fact_spin_tol_empty
2056 : : !! NAME
2057 : : !! get_fact_spin_tol_empty
2058 : : !!
2059 : : !! FUNCTION
2060 : : !!
2061 : : !! INPUTS
2062 : : !!
2063 : : !! OUTPUT
2064 : : !!
2065 : : !! SOURCE
2066 : :
2067 : 1128 : subroutine get_fact_spin_tol_empty(nsppol, nspinor, tol_empty_in, fact_spin, tol_empty)
2068 : :
2069 : : integer,intent(in) :: nsppol, nspinor
2070 : : real(dp),intent(in) :: tol_empty_in
2071 : : real(dp),intent(out) :: fact_spin, tol_empty
2072 : :
2073 : 2242 : select case (nsppol)
2074 : : case (1)
2075 : 1114 : fact_spin = half; tol_empty = tol_empty_in ! below this value the state is assumed empty
2076 [ + + ]: 1114 : if (nspinor == 2) then
2077 : 17 : fact_spin = one; tol_empty = half * tol_empty_in ! below this value the state is assumed empty
2078 : : end if
2079 : : case (2)
2080 : 14 : fact_spin = one; tol_empty = half * tol_empty_in ! to be consistent and obtain similar results if a metallic
2081 : : case default ! spin unpolarized system is treated using nsppol==2
2082 [ + + - ]: 1128 : ABI_BUG(sjoin('Wrong nsppol:', itoa(nsppol)))
2083 : : end select
2084 : :
2085 : 1128 : end subroutine get_fact_spin_tol_empty
2086 : : !!***
2087 : :
2088 : : end module m_occ
2089 : : !!***
|
