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