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