TABLE OF CONTENTS
ABINIT/newocc [ Functions ]
NAME
newocc
FUNCTION
Compute new occupation numbers at each k point, from eigenenergies eigen, according to the smearing scheme defined by occopt (smearing width tsmear and physical temperature tphysel), with the constraint of number of valence electrons per unit cell nelect.
COPYRIGHT
Copyright (C) 1998-2017 ABINIT group (XG) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
INPUTS
eigen(mband*nkpt*nsppol)=eigenvalues (input or init to large number), hartree spinmagntarget=if differ from -99.99_dp, fix the magnetic moment (in Bohr magneton) mband=maximum number of bands nband(nkpt)=number of bands at each k point nelect=number of electrons per unit cell nkpt=number of k points nspinor=number of spinorial components of the wavefunctions nsppol=1 for unpolarized, 2 for spin-polarized occopt=option for occupancies prtvol=control print volume and debugging output stmbias=if non-zero, compute occupation numbers for STM (non-zero around the Fermi energy) NOTE : in this case, only fermie and occ are meaningful outputs. tphysel="physical" electronic temperature with FD occupations tsmear=smearing width (or temperature) wtk(nkpt)=k point weights
OUTPUT
doccde(maxval(nband(:))*nkpt*nsppol)=derivative of occupancies wrt
the energy for each band and k point
entropy= entropy associated with the smearing (adimensional)
fermie= fermi energy (Hartree)
occ(maxval(nband(:))*nkpt*nsppol)=occupancies for each band and k point
SIDE EFFECTS
PARENTS
gstate,m_ebands,respfn,vtorho
CHILDREN
getnel,timab,wrtout
SOURCE
53 #if defined HAVE_CONFIG_H 54 #include "config.h" 55 #endif 56 57 #include "abi_common.h" 58 59 subroutine newocc(doccde,eigen,entropy,fermie,spinmagntarget,mband,nband,& 60 & nelect,nkpt,nspinor,nsppol,occ,occopt,prtvol,stmbias,tphysel,tsmear,wtk) 61 62 use defs_basis 63 use m_profiling_abi 64 use m_errors 65 66 !This section has been created automatically by the script Abilint (TD). 67 !Do not modify the following lines by hand. 68 #undef ABI_FUNC 69 #define ABI_FUNC 'newocc' 70 use interfaces_14_hidewrite 71 use interfaces_18_timing 72 use interfaces_61_occeig, except_this_one => newocc 73 !End of the abilint section 74 75 implicit none 76 77 !Arguments ------------------------------------ 78 !scalars 79 integer,intent(in) :: mband,nkpt,nspinor,nsppol,occopt,prtvol 80 real(dp),intent(in) :: spinmagntarget,nelect,stmbias,tphysel,tsmear 81 real(dp),intent(out) :: entropy,fermie 82 !arrays 83 integer,intent(in) :: nband(nkpt*nsppol) 84 real(dp),intent(in) :: eigen(mband*nkpt*nsppol),wtk(nkpt) 85 real(dp),intent(out) :: doccde(mband*nkpt*nsppol) 86 real(dp),intent(inout) :: occ(mband*nkpt*nsppol) !vz_i 87 88 !Local variables------------------------------- 89 integer,parameter :: niter_max=120,nkpt_max=50,fake_unit=-666,option1=1 90 integer :: cnt,cnt2,cnt3,ib,ii,ik,ikpt,is,isppol,nkpt_eff 91 integer :: sign 92 integer,allocatable :: nbandt(:) 93 real(dp) :: dosdeltae,entropy_tmp,fermihi,fermilo,fermimid,fermimid_tmp 94 real(dp) :: fermi_biased,maxocc 95 real(dp) :: nelect_tmp,nelecthi,nelectlo,nelectmid,nelect_biased 96 real(dp) :: entropyt(2),fermihit(2),fermilot(2),fermimidt(2),nelecthit(2) 97 real(dp) :: nelectlot(2),nelectt(2),tsec(2) 98 real(dp),allocatable :: doccdet(:),eigent(:),occt(:) 99 character(len=500) :: message 100 101 ! ************************************************************************* 102 103 DBG_ENTER("COLL") 104 105 call timab(74,1,tsec) 106 107 !Here treat the case where occopt does not correspond to a metallic occupation scheme 108 if (occopt<3 .or. occopt>8) then 109 write(message,'(a,i0,a)')' occopt= ',occopt,', a value not allowed in newocc.' 110 MSG_BUG(message) 111 end if ! test of metallic occopt 112 113 !Check whether nband is a constant for all k point and spin-pol 114 do isppol=1,nsppol 115 do ikpt=1,nkpt 116 if(nband(ikpt+(isppol-1)*nkpt)/=nband(1))then 117 write(message,'(3a,i0,a,i0,a,i0,a)')& 118 & 'The number of bands must be the same for all k-points ',ch10,& 119 & 'but nband(1)= ',nband(1),' is different of nband(',& 120 & ikpt+(isppol-1)*nkpt,') = ',nband(ikpt+(isppol-1)*nkpt),'.' 121 MSG_BUG(message) 122 end if 123 end do 124 end do 125 126 !Check whether nelect is strictly positive 127 if(nelect<=zero)then 128 write(message,'(3a,es16.8,a)')& 129 & 'nelect must be a positive number, while ',ch10,& 130 & 'the calling routine ask nelect=',nelect,'.' 131 MSG_BUG(message) 132 end if 133 134 maxocc=two/(nsppol*nspinor) 135 !Check whether nelect is coherent with nband (nband(1) is enough, 136 !since it was checked that nband is independent of k-point and spin-pol 137 if( nelect > nband(1)*nsppol*maxocc )then 138 write(message,'(3a,es16.8,a,i0,a,es16.8,a)' )& 139 & 'nelect must be smaller than nband*maxocc, while ',ch10,& 140 & 'the calling routine gives nelect= ',nelect,', nband= ',nband(1),' and maxocc= ',maxocc,'.' 141 MSG_BUG(message) 142 end if 143 144 !Use bissection algorithm to find fermi energy 145 !This choice is due to the fact that it will always give sensible 146 !result (because the answer is bounded, even if the smearing function 147 !is non-monotonic (which is the case for occopt=4 or 6) 148 !Might speed up it, if needed ! 149 150 !Lowest and largest trial fermi energies, and corresponding number of electrons 151 !They are obtained from the smallest or largest eigenenergy, plus a range of 152 !energy that allows for complete occupation of all bands, or, on the opposite, 153 !for zero occupation of all bands (see getnel.f) 154 dosdeltae=zero ! the DOS is not computed, with option=1 155 fermilo=minval(eigen(1:nband(1)*nkpt*nsppol))-6.001_dp*tsmear 156 if(occopt==3)fermilo=fermilo-24.0_dp*tsmear 157 158 call getnel(doccde,dosdeltae,eigen,entropy,fermilo,maxocc,mband,nband,& 159 & nelectlo,nkpt,nsppol,occ,occopt,option1,tphysel,tsmear,fake_unit,wtk) 160 161 fermihi=maxval(eigen(1:nband(1)*nkpt*nsppol))+6.001_dp*tsmear 162 !safety value 163 fermihi = min(fermihi, 1.e6_dp) 164 if(occopt==3)fermihi=fermihi+24.0_dp*tsmear 165 166 call getnel(doccde,dosdeltae,eigen,entropy,fermihi,maxocc,mband,nband,& 167 & nelecthi,nkpt,nsppol,occ,occopt,option1,tphysel,tsmear,fake_unit,wtk) 168 169 !Prepare fixed moment calculation 170 if(abs(spinmagntarget+99.99_dp)>1.0d-10)then 171 sign = 1 172 do is = 1, nsppol 173 fermihit(is) = fermihi 174 fermilot(is) = fermilo 175 nelectt(is) = half*(nelect+sign*spinmagntarget) 176 sign = -sign 177 nelecthit(is) = nelecthi 178 nelectlot(is) = nelectlo 179 end do 180 end if 181 182 !If the target nelect is not between nelectlo and nelecthi, exit 183 if(nelect<nelectlo .or. nelect>nelecthi)then 184 write(message, '(a,a,a,a,d16.8,a,a,d16.8,a,d16.8,a,a,d16.8,a,d16.8)') ch10,& 185 & ' newocc : ',ch10,& 186 & ' The calling routine gives nelect=',nelect,ch10,& 187 & ' The lowest bound is ',fermilo,', with nelect=',nelectlo,ch10,& 188 & ' The highest bound is ',fermihi,', with nelect=',nelecthi 189 call wrtout(std_out,message,'COLL') 190 191 write(message, '(11a)' )& 192 & 'In order to get the right number of electrons,',ch10,& 193 & 'it seems that the Fermi energy must be outside the range',ch10,& 194 & 'of eigenenergies, plus 6 or 30 times the smearing, which is strange.',ch10,& 195 & 'It might be that your number of bands (nband) corresponds to the strictly',ch10,& 196 & 'minimum number of bands to accomodate your electrons (so, OK for an insulator),',ch10,& 197 & 'while you are trying to describe a metal. In this case, increase nband, otherwise ...' 198 MSG_BUG(message) 199 end if 200 201 if( abs(spinmagntarget+99.99_dp) < tol10 ) then 202 203 ! Usual bissection loop 204 do ii=1,niter_max 205 fermimid=(fermihi+fermilo)*half 206 ! Produce nelectmid from fermimid 207 call getnel(doccde,dosdeltae,eigen,entropy,fermimid,maxocc,mband,nband,& 208 & nelectmid,nkpt,nsppol,occ,occopt,option1,tphysel,tsmear,fake_unit,wtk) 209 ! write(std_out,'(a,es24.16,a,es24.16)' )' newocc : from fermi=',fermimid,', getnel gives nelect=',nelectmid 210 if(nelectmid>nelect*(one-tol14))then 211 fermihi=fermimid 212 nelecthi=nelectmid 213 end if 214 if(nelectmid<nelect*(one+tol14))then 215 fermilo=fermimid 216 nelectlo=nelectmid 217 end if 218 if( abs(nelecthi-nelectlo) <= nelect*two*tol14 .or. & 219 & abs(fermihi-fermilo) <= tol14*abs(fermihi+fermilo) ) exit 220 if(ii==niter_max)then 221 write(message,'(a,i0,3a,es22.14,a,es22.14,a)')& 222 & 'It was not possible to find Fermi energy in ',niter_max,' bissections.',ch10,& 223 & 'nelecthi = ',nelecthi,', and nelectlo = ',nelectlo,'.' 224 MSG_BUG(message) 225 end if 226 end do ! End of bissection loop 227 228 fermie=fermimid 229 write(message, '(a,f14.6,a,f14.6,a,a,i4)' ) & 230 & ' newocc: new Fermi energy is ',fermie,' , with nelect=',nelectmid,ch10,& 231 & ' Number of bissection calls =',ii 232 call wrtout(std_out,message,'COLL') 233 234 ! Compute occupation numbers for prtstm/=0, close to the Fermi energy 235 if(abs(stmbias)>tol10)then 236 fermi_biased=fermie-stmbias 237 ABI_ALLOCATE(occt,(mband*nkpt*nsppol)) 238 call getnel(doccde,dosdeltae,eigen,entropy,fermi_biased,maxocc,mband,nband,& 239 & nelect_biased,nkpt,nsppol,occt,occopt,option1,tphysel,tsmear,fake_unit,wtk) 240 occ(:)=occ(:)-occt(:) 241 nelect_biased=abs(nelectmid-nelect_biased) 242 ! Here, arrange to have globally positive occupation numbers, 243 ! irrespective of the stmbias sign 244 if(-stmbias>tol10)occ(:)=-occ(:) 245 ABI_DEALLOCATE(occt) 246 247 write(message,'(a,f14.6)')' newocc : the number of electrons in the STM range is nelect_biased=',nelect_biased 248 call wrtout(std_out,message,'COLL') 249 250 end if 251 252 else ! Calculations with a specified moment 253 254 ! Bissection loop 255 cnt2=0 256 cnt3=0 257 entropy=zero 258 maxocc=one 259 ABI_ALLOCATE(doccdet,(nkpt*mband)) 260 ABI_ALLOCATE(eigent,(nkpt*mband)) 261 ABI_ALLOCATE(occt,(nkpt*mband)) 262 ABI_ALLOCATE(nbandt,(nkpt)) 263 264 do is = 1, nsppol 265 nelect_tmp = nelectt(is) 266 fermihi = fermihit(is) 267 fermilo = fermilot(is) 268 nelecthi = nelecthit(is) 269 nelectlo = nelectlot(is) 270 ! DEBUG 271 ! write(std_out,'(a,i1,3(f8.4,1x))') "Spin, N(spin):", is, nelect, fermihi, fermilo 272 ! write(std_out,'(a,2(f8.4,1x))') "Hi, lo:", nelecthi, nelectlo 273 ! ENDDEBUG 274 275 do ii=1,niter_max 276 fermimid_tmp=(fermihi+fermilo)/2.0_dp 277 ! temporary arrays 278 cnt = 0 279 do ik = 1, nkpt 280 nbandt(ik) = mband 281 do ib = 1, mband 282 cnt = cnt + 1 283 eigent(cnt) = eigen(cnt+cnt2) 284 occt(cnt) = occ(cnt+cnt2) 285 doccdet(cnt) = doccde(cnt+cnt2) 286 end do 287 end do 288 289 ! Produce nelectmid from fermimid 290 call getnel(doccdet,dosdeltae,eigent,entropy_tmp,fermimid_tmp,maxocc,mband,nbandt,& 291 & nelectmid,nkpt,1,occt,occopt,option1,tphysel,tsmear,fake_unit,wtk) 292 entropyt(is) = entropy_tmp 293 fermimidt(is) = fermimid_tmp 294 fermimid = fermimidt(is) 295 ! temporary arrays 296 cnt = 0 297 do ik = 1, nkpt 298 do ib = 1, mband 299 cnt = cnt + 1 300 occ(cnt+cnt2) = occt(cnt) 301 doccde(cnt+cnt2) = doccdet(cnt) 302 end do 303 end do 304 305 ! DEBUG 306 ! write(std_out,'(a,es24.16,a,es24.16)' )& 307 ! & ' newocc : from fermi=',fermimid,', getnel gives nelect=',nelectmid 308 ! ENDDEBUG 309 310 if(nelectmid>=nelect_tmp)then 311 fermihi=fermimid_tmp 312 nelecthi=nelectmid 313 else 314 fermilo=fermimid_tmp 315 nelectlo=nelectmid 316 end if 317 if( abs(nelecthi-nelectlo) <= 1.0d-13 .or. abs(fermihi-fermilo) <= 0.5d-14*abs(fermihi+fermilo) ) exit 318 319 if(ii==niter_max)then 320 write(message,'(a,i3,3a,es22.14,a,es22.14,a)')& 321 & ' It was not possible to find Fermi energy in ',niter_max,' bissections.',ch10,& 322 & ' nelecthi=',nelecthi,', and nelectlo=',nelectlo,'.' 323 MSG_BUG(message) 324 end if 325 end do ! End of bissection loop 326 327 cnt2 = cnt2 + nkpt*mband 328 entropy = entropy + entropyt(is) 329 fermie=fermimid 330 write(message, '(a,i2,a,f14.6,a,f14.6,a,a,i4)' ) & 331 & ' newocc : new Fermi energy for spin ', is, ' is ',fermie,' , with nelect=',nelectmid,ch10,& 332 & ' Number of bissection calls =',ii 333 call wrtout(std_out,message,'COLL') 334 335 end do ! spin 336 337 ABI_DEALLOCATE(doccdet) 338 ABI_DEALLOCATE(eigent) 339 ABI_DEALLOCATE(nbandt) 340 ABI_DEALLOCATE(occt) 341 342 end if ! End of logical on fixed moment calculations 343 344 !write(std_out,*) "kT*Entropy:", entropy*tsmear 345 346 nkpt_eff=nkpt 347 if(prtvol==0)nkpt_eff=min(nkpt_max,nkpt) 348 349 if(nsppol==1)then 350 write(message, '(a,i0,a)' ) & 351 & ' newocc : computed new occ. numbers for occopt= ',occopt,& 352 & ' , spin-unpolarized case. ' 353 call wrtout(std_out,message,'COLL') 354 do ikpt=1,nkpt_eff 355 write(message,'(a,i4,a)' ) ' k-point number ',ikpt,' :' 356 do ii=0,(nband(1)-1)/12 357 write(message,'(12f6.3)') & 358 & occ(1+ii*12+(ikpt-1)*nband(1):min(12+ii*12,nband(1))+(ikpt-1)*nband(1)) 359 call wrtout(std_out,message,'COLL') 360 end do 361 end do 362 if(nkpt/=nkpt_eff)then 363 call wrtout(std_out,' newocc: prtvol=0, stop printing more k-point information','COLL') 364 end if 365 366 ! DEBUG 367 ! write(message, '(a)' ) & 368 ! & ' newocc : corresponding derivatives are ' 369 ! call wrtout(std_out,message,'COLL') 370 ! do ikpt=1,nkpt_eff 371 ! write(message,'(a,i4,a)' ) ' k-point number ',ikpt,' :' 372 ! do ii=0,(nband(1)-1)/12 373 ! write(message,'(12f6.1)') & 374 ! & doccde(1+ii*12+(ikpt-1)*nband(1):min(12+ii*12,nband(1))+(ikpt-1)*nband(1)) 375 ! call wrtout(std_out,message,'COLL') 376 ! end do 377 ! end do 378 ! if(nkpt/=nkpt_eff)then 379 ! write(message,'(a)') & 380 ! & ' newocc : prtvol=0, stop printing more k-point informations' 381 ! call wrtout(std_out,message,'COLL') 382 ! end if 383 ! ENDDEBUG 384 else 385 write(message, '(a,i0,a,a)' ) & 386 & ' newocc : computed new occupation numbers for occopt= ',occopt,& 387 & ch10,' (1) spin up values ' 388 call wrtout(std_out,message,'COLL') 389 do ikpt=1,nkpt_eff 390 write(message,'(a,i0,a)' ) ' k-point number ',ikpt,':' 391 do ii=0,(nband(1)-1)/12 392 write(message,'(12f6.3)') & 393 & occ(1+ii*12+(ikpt-1)*nband(1):min(12+ii*12,nband(1))+(ikpt-1)*nband(1)) 394 call wrtout(std_out,message,'COLL') 395 end do 396 end do 397 if(nkpt/=nkpt_eff)then 398 call wrtout(std_out,'newocc: prtvol=0, stop printing more k-point information','COLL') 399 end if 400 401 call wrtout(std_out,' (2) spin down values ','COLL') 402 do ikpt=1,nkpt_eff 403 do ii=0,(nband(1)-1)/12 404 write(message,'(12f6.3)') & 405 & occ( 1+ii*12+(ikpt-1+nkpt)*nband(1) : & 406 & min(12+ii*12,nband(1))+(ikpt-1+nkpt)*nband(1) ) 407 call wrtout(std_out,message,'COLL') 408 end do 409 end do 410 if(nkpt/=nkpt_eff)then 411 call wrtout(std_out,' newocc: prtvol=0, stop printing more k-point information','COLL') 412 end if 413 414 end if ! End choice based on spin 415 416 call timab(74,2,tsec) 417 418 DBG_EXIT("COLL") 419 420 end subroutine newocc