TABLE OF CONTENTS
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
COPYRIGHT
Copyright (C) 1998-2017 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 . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
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 (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
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(tseamr^3) $
PARENTS
cchi0q0_intraband,clnup1,conducti_nc,dfpt_looppert,m_ebands,newocc
CHILDREN
dos_hdr_write,init_occ_ent,splfit,wrtout
SOURCE
71 #if defined HAVE_CONFIG_H 72 #include "config.h" 73 #endif 74 75 #include "abi_common.h" 76 77 subroutine getnel(doccde,dosdeltae,eigen,entropy,fermie,maxocc,mband,nband,& 78 & nelect,nkpt,nsppol,occ,occopt,option,tphysel,tsmear,unitdos,wtk) 79 80 use defs_basis 81 use m_profiling_abi 82 use m_errors 83 use m_splines 84 85 use m_fstrings, only : sjoin, itoa 86 87 !This section has been created automatically by the script Abilint (TD). 88 !Do not modify the following lines by hand. 89 #undef ABI_FUNC 90 #define ABI_FUNC 'getnel' 91 use interfaces_14_hidewrite 92 use interfaces_61_occeig, except_this_one => getnel 93 !End of the abilint section 94 95 implicit none 96 97 !Arguments ------------------------------------ 98 !scalars 99 integer,intent(in) :: mband,nkpt,nsppol,occopt,option,unitdos 100 real(dp),intent(in) :: dosdeltae,fermie,maxocc,tphysel,tsmear 101 real(dp),intent(out) :: entropy,nelect 102 !arrays 103 integer,intent(in) :: nband(nkpt*nsppol) 104 real(dp),intent(in) :: eigen(mband*nkpt*nsppol),wtk(nkpt) 105 real(dp),intent(out) :: doccde(mband*nkpt*nsppol) !vz_i 106 real(dp),intent(inout) :: occ(mband*nkpt*nsppol) !vz_i 107 108 !Local variables------------------------------- 109 ! nptsdiv2 is the number of integration points, divided by 2. 110 ! tratio = ratio tsmear/tphysel for convoluted smearing function 111 ! save values so we can impose recalculation of smdfun when 112 ! the smearing or electronic temperature change between 113 ! datasets 114 ! corresponds roughly to delta_FD (maxFDarg) = 1.0d-100 115 ! 116 ! return fermi-dirac smearing function analytically 117 ! 118 ! real(dp) :: smdFD 119 ! smdFD (tt) = 1.0_dp / (exp(-tt/2.0_dp) + exp(tt/2.0_dp))**2 120 !scalars 121 ! TODO: This parameter is defined in init_occ_ent but we cannot call the 122 ! routine to get this value since the same variable is used to dimension the 123 ! arrays! This Constants should be stored somewhere in a module. 124 integer,parameter :: nptsdiv2_def=6000 125 integer,parameter :: prtdos1=1 126 integer :: bantot,iband,iene,ikpt,index,index_start,isppol 127 integer :: nene,nptsdiv2 128 real(dp) :: buffer,deltaene,dosdbletot,doshalftot,dostot 129 real(dp) :: enemax,enemin,enex,intdostot,limit,tsmearinv 130 character(len=500) :: message 131 !arrays 132 real(dp),allocatable :: entfun(:,:),occfun(:,:) 133 real(dp),allocatable :: smdfun(:,:),xgrid(:) 134 real(dp),allocatable :: arg(:),derfun(:),dos(:),dosdble(:),doshalf(:),ent(:) 135 real(dp),allocatable :: intdos(:) 136 137 ! ************************************************************************* 138 139 DBG_ENTER("COLL") 140 141 if(option/=1 .and. option/=2)then 142 MSG_BUG(sjoin('Option must be either 1 or 2. It is:', itoa(option))) 143 end if 144 145 !Initialize the occupation function and generalized entropy function, 146 !at the beginning, or if occopt changed 147 148 !Just get the number nptsdiv2 and allocate entfun, occfun, 149 !smdfun and xgrid accordingly 150 nptsdiv2 = nptsdiv2_def 151 152 ! call init_occ_ent(entfun, limit, & 153 !& nptsdiv2, occfun, occopt, -1, smdfun, tphysel, & 154 !& tsmear, tsmearinv, xgrid) 155 156 ABI_ALLOCATE(entfun,(-nptsdiv2:nptsdiv2,2)) 157 ABI_ALLOCATE(occfun,(-nptsdiv2:nptsdiv2,2)) 158 ABI_ALLOCATE(smdfun,(-nptsdiv2:nptsdiv2,2)) 159 ABI_ALLOCATE(xgrid,(-nptsdiv2:nptsdiv2)) 160 161 !Call to init_occ_ent 162 call init_occ_ent(entfun, limit, & 163 & nptsdiv2, occfun, occopt, option, smdfun, tphysel, & 164 & tsmear, tsmearinv, xgrid) 165 166 !The initialisation of occfun and entfun is done 167 168 !--------------------------------------------------------------------- 169 170 !write(std_out,*)' getnel : debug tphysel, tsmear = ', tphysel, tsmear 171 bantot=sum(nband(:)) 172 173 ABI_ALLOCATE(arg,(bantot)) 174 ABI_ALLOCATE(derfun,(bantot)) 175 ABI_ALLOCATE(ent,(bantot)) 176 177 if(option==1)then 178 !normal evaluation of occupations and entropy 179 180 ! Compute the arguments of the occupation and entropy functions 181 arg(:)=(fermie-eigen(1:bantot))*tsmearinv 182 183 ! Compute the values of the occupation function, and the entropy function 184 ! Note : splfit also takes care of the points outside of the interval, 185 ! and assign to them the value of the closest extremal point, 186 ! which is what is needed here. 187 188 call splfit(xgrid,doccde,occfun,1,arg,occ,(2*nptsdiv2+1),bantot) 189 call splfit(xgrid,derfun,entfun,0,arg,ent,(2*nptsdiv2+1),bantot) 190 191 ! Normalize occ and ent, and sum number of electrons and entropy 192 nelect=zero; entropy=zero 193 index=0 194 do isppol=1,nsppol 195 do ikpt=1,nkpt 196 do iband=1,nband(ikpt+nkpt*(isppol-1)) 197 index=index+1 198 ent(index)=ent(index)*maxocc 199 occ(index)=occ(index)*maxocc 200 doccde(index)=-doccde(index)*maxocc*tsmearinv 201 entropy=entropy+wtk(ikpt)*ent(index) 202 nelect=nelect+wtk(ikpt)*occ(index) 203 end do 204 end do 205 end do 206 207 ! write(std_out,*) ' getnel : debug wtk, occ, eigen = ', wtk, occ, eigen 208 ! write(std_out,*)xgrid(-nptsdiv2),xgrid(nptsdiv2) 209 ! write(std_out,*)'fermie',fermie 210 ! do ii=1,bantot 211 ! write(std_out,*)ii,arg(ii),doccde(ii) 212 ! end do 213 ! write(std_out,*)'eigen',eigen(:) 214 ! write(std_out,*)'arg',arg(:) 215 ! write(std_out,*)'occ',occ(:) 216 ! write(std_out,*)'nelect',nelect 217 218 else if(option==2)then 219 ! evaluate DOS for smearing, half smearing, and double. 220 221 buffer=limit/tsmearinv*.5_dp 222 223 ! A Similar section is present is dos_calcnwrite. Should move all DOS stuff to m_ebands 224 ! Choose the lower and upper energies 225 enemax=maxval(eigen(1:bantot))+buffer 226 enemin=minval(eigen(1:bantot))-buffer 227 228 ! Extend the range to a nicer value 229 enemax=0.1_dp*ceiling(enemax*10._dp) 230 enemin=0.1_dp*floor(enemin*10._dp) 231 232 ! Choose the energy increment 233 if(abs(dosdeltae)<tol10)then 234 deltaene=0.001_dp 235 if(prtdos1>=2)deltaene=0.0005_dp ! Higher resolution possible (and wanted) for tetrahedron 236 else 237 deltaene=dosdeltae 238 end if 239 nene=nint((enemax-enemin)/deltaene)+1 240 241 ! Write the header of the DOS file, and also decides the energy range and increment 242 call dos_hdr_write(buffer,deltaene,dosdeltae,eigen,enemax,enemin,fermie,mband,nband,nene,& 243 & nkpt,nsppol,occopt,prtdos1,tphysel,tsmear,unitdos) 244 245 ABI_ALLOCATE(dos,(bantot)) 246 ABI_ALLOCATE(dosdble,(bantot)) 247 ABI_ALLOCATE(doshalf,(bantot)) 248 ABI_ALLOCATE(intdos,(bantot)) 249 250 do isppol=1,nsppol 251 252 if (nsppol==2) then 253 if(isppol==1) write(message,'(a,16x,a)') '#','Spin-up DOS' 254 if(isppol==2) write(message,'(2a,16x,a)') ch10,'#','Spin-dn DOS ' 255 call wrtout(unitdos,message,'COLL') 256 end if 257 index_start=0 258 if(isppol==2)then 259 do ikpt=1,nkpt 260 index_start=index_start+nband(ikpt) 261 end do 262 end if 263 264 enex=enemin 265 do iene=1,nene 266 267 ! Compute the arguments of the dos and occupation function 268 arg(:)=(enex-eigen(1:bantot))*tsmearinv 269 270 call splfit(xgrid,derfun,smdfun,0,arg,dos,(2*nptsdiv2+1),bantot) 271 call splfit(xgrid,derfun,occfun,0,arg,intdos,(2*nptsdiv2+1),bantot) 272 ! Also compute the dos with tsmear halved and doubled 273 arg(:)=arg(:)*2.0_dp 274 call splfit(xgrid,derfun,smdfun,0,arg,doshalf,(2*nptsdiv2+1),bantot) 275 ! Since arg was already doubled, must divide by four 276 arg(:)=arg(:)*0.25_dp 277 call splfit(xgrid,derfun,smdfun,0,arg,dosdble,(2*nptsdiv2+1),bantot) 278 279 ! Now, accumulate the contribution from each eigenenergy 280 dostot=zero 281 intdostot=zero 282 doshalftot=zero 283 dosdbletot=zero 284 index=index_start 285 286 ! write(std_out,*)' eigen, arg, dos, intdos, doshalf, dosdble' 287 do ikpt=1,nkpt 288 do iband=1,nband(ikpt+nkpt*(isppol-1)) 289 index=index+1 290 dostot=dostot+wtk(ikpt)*maxocc*dos(index)*tsmearinv 291 intdostot=intdostot+wtk(ikpt)*maxocc*intdos(index) 292 doshalftot=doshalftot+wtk(ikpt)*maxocc*doshalf(index)*tsmearinv*2.0_dp 293 dosdbletot=dosdbletot+wtk(ikpt)*maxocc*dosdble(index)*tsmearinv*0.5_dp 294 end do 295 end do 296 297 ! Print the data for this energy 298 write(unitdos, '(f8.3,2f14.6,2f14.3)' )enex,dostot,intdostot,doshalftot,dosdbletot 299 300 enex=enex+deltaene 301 end do ! iene 302 end do ! isppol 303 304 ABI_DEALLOCATE(dos) 305 ABI_DEALLOCATE(dosdble) 306 ABI_DEALLOCATE(doshalf) 307 ABI_DEALLOCATE(intdos) 308 309 ! MG: It does not make sense to close the unit here since the routines 310 ! did not open the file here! 311 ! Close the DOS file 312 close(unitdos) 313 end if 314 315 ABI_DEALLOCATE(arg) 316 ABI_DEALLOCATE(derfun) 317 ABI_DEALLOCATE(ent) 318 ABI_DEALLOCATE(entfun) 319 ABI_DEALLOCATE(occfun) 320 ABI_DEALLOCATE(smdfun) 321 ABI_DEALLOCATE(xgrid) 322 323 DBG_EXIT("COLL") 324 325 end subroutine getnel