TABLE OF CONTENTS
ABINIT/emispec_paw [ Functions ]
NAME
emispec_paw
FUNCTION
This program computes the elements of the emission spectra from the Kubo-Greenwood formula for PAW formalism largely inspired from the conducti_core_paw routine
COPYRIGHT
Copyright (C) 2002-2018 ABINIT group (SM, SVinko) 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/Infos/contributors .
INPUTS
(main routine)
OUTPUT
(main routine)
NOTES
bantot dom=frequency range eigen0(mband*nkpt_rbz*nsppol)=GS eigenvalues at k (hartree). ecut=kinetic energy planewave cutoff (hartree). fermie= fermi energy (Hartree) mom=number of frequency for conductivity computation mband=maximum number of bands. natom = number of atoms in the unit cell. nband(nkpt*nsppol)=number of bands at each RF k point for each spin. nkpt=number of k points in the IBZ for this perturbation ngfft(3)=integer fft box dimensions. nspinor=number of spinorial components of the wavefunctions. nsppol=1 for unpolarized, 2 for spin-polarized. ntypat = number of atom types. occ(mband*nkpt*nsppol)=occupation number for each band and k. occopt==option for occupancies psinablapsi2(2,3,mband,nphicor,natom)Matrix elements = <Phi_core|Nabla|Phi_i> rmet(3,3)=real space metric ($\textrm{bohr}^{2}$).sigx(mom,nphicor)) rprimd(3,3)=real space primitive translations. of primitive translations. ucvol=unit cell volume in ($\textrm{bohr}^{3}$). wind=frequency windows for computations of sigma wtk(nkpt)=weight assigned to each k point.
PARENTS
conducti
CHILDREN
hdr_free,hdr_io,hdr_read_from_fname,metric,wffclose,wffopen
SOURCE
57 #if defined HAVE_CONFIG_H 58 #include "config.h" 59 #endif 60 61 #include "abi_common.h" 62 63 subroutine emispec_paw(filnam,filnam_out,mpi_enreg) 64 65 use defs_basis 66 use defs_abitypes 67 use m_xmpi 68 use m_wffile 69 use m_profiling_abi 70 use m_errors 71 use m_hdr 72 73 use m_io_tools, only : open_file, get_unit 74 use m_geometry, only : metric 75 76 !This section has been created automatically by the script Abilint (TD). 77 !Do not modify the following lines by hand. 78 #undef ABI_FUNC 79 #define ABI_FUNC 'emispec_paw' 80 !End of the abilint section 81 82 implicit none 83 84 !Arguments ----------------------------------- 85 !scalars 86 character(len=fnlen) :: filnam,filnam_out 87 type(MPI_type),intent(in) :: mpi_enreg 88 89 !Local variables------------------------------- 90 !scalars 91 integer,parameter :: master=0 92 integer :: iomode,atnbr,bantot,bdtot_index 93 integer :: fform2,headform,iatom,iband,icor,ierr,ikpt 94 integer :: iom,isppol,l1,mband,me,mom 95 integer :: natom,nband_k,nkpt,nphicor,nspinor,nsppol,ntypat 96 integer :: occopt,rdwr,spaceComm,iunt,ems_unt,opt2_unt 97 real(dp) :: del,ecut,fermie 98 real(dp) :: omin,omax,dom,oml 99 real(dp) :: Tatm,tsmear,ucvol 100 character(len=fnlen) :: filnam2,filnam_gen 101 character(len=500) :: msg 102 type(hdr_type) :: hdr 103 type(wffile_type) :: wff2 104 !arrays 105 integer,allocatable :: nband(:),ncor(:),lcor(:) 106 real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3) 107 real(dp),allocatable :: dhdk2_g(:,:,:) 108 real(dp),allocatable :: eig0_k(:),eigen0(:) 109 real(dp),allocatable :: energy_cor(:) 110 real(dp),allocatable :: occ(:),occ_k(:),oml1(:) 111 real(dp),allocatable :: psinablapsi2(:,:,:,:,:) 112 real(dp),allocatable :: sigx(:,:,:),sigx_av(:,:),wtk(:) 113 114 ! ********************************************************************************* 115 ABI_UNUSED(mpi_enreg%paral_kgb) 116 117 !Read output file name 118 !write(std_out,'(a)')' emispec_paw : Please, give the name of the output file ...' 119 !read(std_in, '(a)') filnam_out 120 !write(std_out,'(a)')' The name of the output file is :',trim(filnam_out) 121 !Read data file 122 if (open_file(filnam,msg,newunit=iunt,form='formatted',action="read") /=0) then 123 MSG_ERROR(msg) 124 end if 125 rewind(iunt) 126 read(iunt,*) 127 read(iunt,'(a)')filnam_gen ! generic name for the files 128 filnam2=trim(filnam_gen)//'_OPT2' 129 !Read size of the frequency range 130 read(iunt,*) dom,omin,omax,mom,atnbr 131 close(iunt) 132 133 write(std_out,'(a,i8,3f10.5,a)')' npts,omin,omax,width =',mom,omin,omax,dom,' Ha' 134 write(std_out,'(a)')'--------------------------------------------' 135 write(std_out,'(a,i4)') 'selected atom for X ray emission',atnbr 136 write(std_out,'(a)')'--------------------------------------------' 137 138 !Open the Wavefunction and optic files 139 !These default values are typical of sequential use 140 iomode=IO_MODE_FORTRAN; spaceComm=xmpi_comm_self; me=0 141 142 ! Read the header of the OPT2 file. 143 call hdr_read_from_fname(hdr, filnam2, fform2, spaceComm) 144 call hdr_free(hdr) 145 146 if (fform2 /= 611) then 147 MSG_ERROR("Abinit8 requires an OPT2 file with fform = 611") 148 end if 149 150 !Open the optic files 151 opt2_unt = get_unit() 152 call WffOpen(iomode,spaceComm,filnam2,ierr,wff2,master,me,opt2_unt) 153 154 !Read the header 155 rdwr=1 156 call hdr_io(fform2,hdr,rdwr,wff2) 157 158 !Extract info from the header 159 headform=hdr%headform 160 bantot=hdr%bantot 161 ecut=hdr%ecut_eff 162 natom=hdr%natom 163 nkpt=hdr%nkpt 164 nspinor=hdr%nspinor 165 nsppol=hdr%nsppol 166 ntypat=hdr%ntypat 167 occopt=hdr%occopt 168 rprimd(:,:)=hdr%rprimd(:,:) 169 ABI_ALLOCATE(nband,(nkpt*nsppol)) 170 ABI_ALLOCATE(occ,(bantot)) 171 ABI_ALLOCATE(wtk,(nkpt)) 172 fermie=hdr%fermie 173 tsmear=hdr%tsmear 174 occ(1:bantot)=hdr%occ(1:bantot) 175 wtk(1:nkpt)=hdr%wtk(1:nkpt) 176 nband(1:nkpt*nsppol)=hdr%nband(1:nkpt*nsppol) 177 178 !Get mband, as the maximum value of nband(nkpt) 179 mband=maxval(nband(:)) 180 181 write(std_out,*) 182 write(std_out,'(a,3f10.5,a)' )' rprimd(bohr) =',rprimd(1,1:3) 183 write(std_out,'(a,3f10.5,a)' )' ',rprimd(2,1:3) 184 write(std_out,'(a,3f10.5,a)' )' ',rprimd(3,1:3) 185 write(std_out,'(a,i8)') ' natom =',natom 186 write(std_out,'(a,3i8)') ' nkpt,mband,nsppol =',nkpt,mband,nsppol 187 write(std_out, '(a, f10.5,a)' ) ' ecut =',ecut,' Ha' 188 write(std_out,'(a,f10.5,a,f10.5,a)' )' fermie =',fermie,' Ha',fermie*Ha_eV,' eV' 189 Tatm=tsmear*Ha_K 190 write(std_out,'(a,f12.5,a,f12.5,a)') ' Temp =',tsmear,' Ha ',Tatm,' Kelvin' 191 192 ABI_ALLOCATE(eigen0,(mband*nkpt*nsppol)) 193 read(opt2_unt)(eigen0(iband),iband=1,mband*nkpt*nsppol) 194 195 write(std_out,'(a)')'--------------------------------------------' 196 read(opt2_unt) nphicor 197 write(std_out,'(a,i4)') 'Number of core orbitals nc=',nphicor 198 ABI_ALLOCATE(ncor,(nphicor)) 199 ABI_ALLOCATE(lcor,(nphicor)) 200 ABI_ALLOCATE(energy_cor,(nphicor)) 201 do icor=1,nphicor 202 read(opt2_unt) ncor(icor),lcor(icor),energy_cor(icor) 203 write(std_out,'(a,2i4,f10.5)') ' n, l, Energy (Ha): ',ncor(icor),lcor(icor),energy_cor(icor) 204 end do 205 write(std_out,'(a)')'--------------------------------------------' 206 207 !--------------------------------------------------------------------------------- 208 !gmet inversion to get ucvol 209 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 210 211 !--------------------------------------------------------------------------------- 212 !frequency range fixed by the occupation of the valence states 213 bdtot_index=0 214 omax=0 215 do isppol=1,nsppol 216 do ikpt=1,nkpt 217 nband_k=nband(ikpt+(isppol-1)*nkpt) 218 do iband=1,nband_k 219 if(eigen0(1+bdtot_index)>=omax) omax=eigen0(1+bdtot_index) 220 bdtot_index=bdtot_index+1 221 end do 222 end do 223 end do 224 omin=eigen0(1) 225 del=(omax-omin)/(mom-1) 226 ABI_ALLOCATE(oml1,(mom)) 227 do iom=1,mom 228 oml1(iom)=omin+dble(iom-1)*del 229 end do 230 write(std_out,*) 'Valence state orbital energies: omin,omax',omin,omax 231 ABI_ALLOCATE(sigx,(natom,mom,nphicor)) 232 ABI_ALLOCATE(sigx_av,(mom,nphicor)) 233 !--------------------------------------------------------------------------------- 234 !emission X ------- 235 ! 236 ABI_ALLOCATE(psinablapsi2,(2,3,mband,nphicor,natom)) 237 sigx=zero 238 sigx_av=zero 239 bdtot_index = 0 240 241 !LOOP OVER SPINS 242 do isppol=1,nsppol 243 244 ! LOOP OVER K-POINTS 245 do ikpt=1,nkpt 246 nband_k=nband(ikpt+(isppol-1)*nkpt) 247 ABI_ALLOCATE(eig0_k,(nband_k)) 248 ABI_ALLOCATE(occ_k,(nband_k)) 249 ABI_ALLOCATE(dhdk2_g,(natom,nband_k,nphicor)) 250 251 dhdk2_g = zero 252 psinablapsi2=zero 253 254 ! eigenvalue for k-point 255 eig0_k(:)=eigen0(1+bdtot_index:nband_k+bdtot_index) 256 257 ! occupation numbers for k-point 258 occ_k(:)=occ(1+bdtot_index:nband_k+bdtot_index) 259 260 ! dipole matrix elements 261 do iatom=1,natom 262 read(opt2_unt) ((psinablapsi2(1:2,1,iband,icor,iatom),iband=1,nband_k),icor=1,nphicor) 263 read(opt2_unt) ((psinablapsi2(1:2,2,iband,icor,iatom),iband=1,nband_k),icor=1,nphicor) 264 read(opt2_unt) ((psinablapsi2(1:2,3,iband,icor,iatom),iband=1,nband_k),icor=1,nphicor) 265 end do 266 267 ! loops over atom, bands, core states 268 do iatom=1,natom 269 do iband=1,nband_k 270 do icor=1,nphicor 271 272 do l1=1,3 273 dhdk2_g(iatom,iband,icor)=dhdk2_g(iatom,iband,icor)+( & 274 & psinablapsi2(1,l1,iband,icor,iatom)*psinablapsi2(1,l1,iband,icor,iatom) & 275 & +psinablapsi2(2,l1,iband,icor,iatom)*psinablapsi2(2,l1,iband,icor,iatom)) 276 end do 277 278 ! emission for each omega 279 do iom=1,mom 280 oml=-energy_cor(icor)-oml1(iom) 281 sigx(iatom,iom,icor)=sigx(iatom,iom,icor)+ wtk(ikpt)*dhdk2_g(iatom,iband,icor)& 282 & *occ_k(iband)/oml*dexp(-((-energy_cor(icor)+eig0_k(iband)-oml)/dom)**2) 283 end do 284 end do !icor 285 end do !iband 286 end do !iatom 287 bdtot_index=bdtot_index+nband_k 288 ABI_DEALLOCATE(eig0_k) 289 ABI_DEALLOCATE(occ_k) 290 ABI_DEALLOCATE(dhdk2_g) 291 ! end loop over k 292 end do 293 ! end loop over spins 294 end do 295 ABI_DEALLOCATE(psinablapsi2) 296 297 do iatom=1,natom 298 do icor=1,nphicor 299 do iom=1,mom 300 if(sigx(iatom,iom,icor)<=tol16) sigx(iatom,iom,icor)=zero 301 end do 302 end do 303 end do ! iatom 304 305 sigx=sigx*two_pi*third*dble(natom)/(dom*ucvol)*half/dsqrt(pi) 306 307 do icor=1,nphicor 308 do iom=1,mom 309 do iatom=1,natom 310 sigx_av(iom,icor) =sigx_av(iom,icor)+sigx(iatom,iom,icor)/dble(natom) 311 end do 312 end do 313 end do 314 315 if (open_file(trim(filnam_out)//'_emisX',msg,newunit=ems_unt,form='formatted', action="write") /= 0) then 316 MSG_ERROR(msg) 317 end if 318 do iom=1,mom 319 write(ems_unt,'(9(1x,e15.8))') & 320 & ((-energy_cor(icor)+oml1(iom)),sigx_av(iom,icor),sigx(atnbr,iom,icor),icor=1,nphicor) 321 end do 322 close(ems_unt) 323 324 call WffClose(wff2,ierr) 325 326 ABI_DEALLOCATE(sigx) 327 ABI_DEALLOCATE(sigx_av) 328 ABI_DEALLOCATE(ncor) 329 ABI_DEALLOCATE(lcor) 330 ABI_DEALLOCATE(energy_cor) 331 ABI_DEALLOCATE(eigen0) 332 ABI_DEALLOCATE(nband) 333 ABI_DEALLOCATE(oml1) 334 ABI_DEALLOCATE(occ) 335 ABI_DEALLOCATE(wtk) 336 337 call hdr_free(hdr) 338 339 end subroutine emispec_paw