TABLE OF CONTENTS
ABINIT/conducti_paw_core [ Functions ]
NAME
conducti_paw_core
FUNCTION
This program computes the elements of the optical frequency dependent conductivity tensor and the conductivity along the three principal axes from the Kubo-Greenwood formula for PAW formalism
COPYRIGHT
Copyright (C) 2002-2018 ABINIT group (VRecoules, PGhosh) 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 conducti_paw_core(filnam,filnam_out,mpi_enreg) 64 65 use defs_basis 66 use defs_abitypes 67 use m_wffile 68 use m_xmpi 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 'conducti_paw_core' 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,opt2_unt,sigx_unt 97 real(dp) :: del,diff_occ,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)')' conducti_paw_core : Please, give the name of the output file ...' 119 !read(5, '(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", status="old") /= 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 spectro X',atnbr 136 write(std_out,'(a)')'--------------------------------------------' 137 138 !These default values are typical of sequential use 139 iomode=IO_MODE_FORTRAN; spaceComm=xmpi_comm_self; me=0 140 141 ! Read the header of the OPT2 file. 142 call hdr_read_from_fname(hdr, filnam2, fform2, spaceComm) 143 call hdr_free(hdr) 144 145 if (fform2 /= 611) then 146 MSG_ERROR("Abinit8 requires an OPT2 file with fform = 611") 147 end if 148 149 !Open the optic files 150 opt2_unt = get_unit() 151 call WffOpen(iomode,spaceComm,filnam2,ierr,wff2,master,me,opt2_unt) 152 153 !Read the header 154 rdwr=1 155 call hdr_io(fform2,hdr,rdwr,wff2) 156 157 !Extract info from the header 158 headform=hdr%headform 159 bantot=hdr%bantot 160 ecut=hdr%ecut_eff 161 natom=hdr%natom 162 nkpt=hdr%nkpt 163 nspinor=hdr%nspinor 164 nsppol=hdr%nsppol 165 ntypat=hdr%ntypat 166 occopt=hdr%occopt 167 rprimd(:,:)=hdr%rprimd(:,:) 168 ABI_ALLOCATE(nband,(nkpt*nsppol)) 169 ABI_ALLOCATE(occ,(bantot)) 170 ABI_ALLOCATE(wtk,(nkpt)) 171 fermie=hdr%fermie 172 tsmear=hdr%tsmear 173 occ(1:bantot)=hdr%occ(1:bantot) 174 wtk(1:nkpt)=hdr%wtk(1:nkpt) 175 nband(1:nkpt*nsppol)=hdr%nband(1:nkpt*nsppol) 176 177 !Get mband, as the maximum value of nband(nkpt) 178 mband=maxval(nband(:)) 179 180 write(std_out,*) 181 write(std_out,'(a,3f10.5,a)' )' rprimd(bohr) =',rprimd(1,1:3) 182 write(std_out,'(a,3f10.5,a)' )' ',rprimd(2,1:3) 183 write(std_out,'(a,3f10.5,a)' )' ',rprimd(3,1:3) 184 write(std_out,'(a,i8)') ' natom =',natom 185 write(std_out,'(a,3i8)') ' nkpt,mband,nsppol =',nkpt,mband,nsppol 186 write(std_out, '(a, f10.5,a)' ) ' ecut =',ecut,' Ha' 187 write(std_out,'(a,f10.5,a,f10.5,a)' )' fermie =',fermie,' Ha',fermie*Ha_eV,' eV' 188 Tatm=tsmear*Ha_K 189 write(std_out,'(a,f12.5,a,f12.5,a)') ' Temp =',tsmear,' Ha ',Tatm,' Kelvin' 190 191 ABI_ALLOCATE(eigen0,(mband*nkpt*nsppol)) 192 read(opt2_unt)(eigen0(iband),iband=1,mband*nkpt*nsppol) 193 ! 194 write(std_out,'(a)')'--------------------------------------------' 195 read(opt2_unt) nphicor 196 write(std_out,'(a,i4)') 'Number of core orbitals nc=',nphicor 197 ABI_ALLOCATE(ncor,(nphicor)) 198 ABI_ALLOCATE(lcor,(nphicor)) 199 ABI_ALLOCATE(energy_cor,(nphicor)) 200 do icor=1,nphicor 201 read(opt2_unt) ncor(icor),lcor(icor),energy_cor(icor) 202 write(std_out,'(a,2i4,f10.5)') ' n, l, Energy (Ha): ',ncor(icor),lcor(icor),energy_cor(icor) 203 end do 204 write(std_out,'(a)')'--------------------------------------------' 205 206 !--------------------------------------------------------------------------------- 207 !gmet inversion to get ucvol 208 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 209 210 !--------------------------------------------------------------------------------- 211 !size of the frequency range 212 del=(omax-omin)/(mom-1) 213 ABI_ALLOCATE(oml1,(mom)) 214 do iom=1,mom 215 oml1(iom)=omin+dble(iom-1)*del 216 end do 217 218 ABI_ALLOCATE(sigx,(natom,mom,nphicor)) 219 ABI_ALLOCATE(sigx_av,(mom,nphicor)) 220 !--------------------------------------------------------------------------------- 221 !SpectroX ------- 222 ! 223 ABI_ALLOCATE(psinablapsi2,(2,3,mband,nphicor,natom)) 224 sigx=zero 225 sigx_av=zero 226 bdtot_index = 0 227 228 !LOOP OVER SPINS 229 do isppol=1,nsppol 230 231 ! LOOP OVER K-POINTS 232 do ikpt=1,nkpt 233 nband_k=nband(ikpt+(isppol-1)*nkpt) 234 ABI_ALLOCATE(eig0_k,(nband_k)) 235 ABI_ALLOCATE(occ_k,(nband_k)) 236 ABI_ALLOCATE(dhdk2_g,(natom,nband_k,nphicor)) 237 238 dhdk2_g = zero 239 psinablapsi2=zero 240 241 ! eigenvalue for k-point 242 eig0_k(:)=eigen0(1+bdtot_index:nband_k+bdtot_index) 243 244 ! occupation numbers for k-point 245 occ_k(:)=occ(1+bdtot_index:nband_k+bdtot_index) 246 247 do iatom=1,natom 248 ! first derivative eigenvalues for k-point 249 read(opt2_unt) ((psinablapsi2(1:2,1,iband,icor,iatom),iband=1,nband_k),icor=1,nphicor) 250 read(opt2_unt) ((psinablapsi2(1:2,2,iband,icor,iatom),iband=1,nband_k),icor=1,nphicor) 251 read(opt2_unt) ((psinablapsi2(1:2,3,iband,icor,iatom),iband=1,nband_k),icor=1,nphicor) 252 end do 253 254 ! LOOP OVER ATOMS/BANDS 255 do iatom=1,natom 256 do iband=1,nband_k 257 do icor=1,nphicor 258 259 do l1=1,3 260 dhdk2_g(iatom,iband,icor)=dhdk2_g(iatom,iband,icor)+( & 261 & psinablapsi2(1,l1,iband,icor,iatom)*psinablapsi2(1,l1,iband,icor,iatom) & 262 & +psinablapsi2(2,l1,iband,icor,iatom)*psinablapsi2(2,l1,iband,icor,iatom)) 263 end do 264 265 diff_occ = (two/dble(nsppol))-occ_k(iband) 266 ! Spectro for each omega 267 omin = -1.0 268 do iom=1,mom 269 oml=-energy_cor(icor)+oml1(iom)+omin 270 sigx(iatom,iom,icor)=sigx(iatom,iom,icor)+ wtk(ikpt)*dhdk2_g(iatom,iband,icor)& 271 & *(diff_occ)/oml*dexp(-((-energy_cor(icor)+eig0_k(iband)-oml)/dom)**2) 272 end do 273 end do !icor 274 end do !iband 275 end do !iatom 276 bdtot_index=bdtot_index+nband_k 277 ABI_DEALLOCATE(eig0_k) 278 ABI_DEALLOCATE(occ_k) 279 ABI_DEALLOCATE(dhdk2_g) 280 ! end loop over k 281 end do 282 ! end loop over spins 283 end do 284 ABI_DEALLOCATE(psinablapsi2) 285 286 do iatom=1,natom 287 do icor=1,nphicor 288 do iom=1,mom 289 if(sigx(iatom,iom,icor)<=tol16) sigx(iatom,iom,icor)=zero 290 end do 291 end do 292 end do ! iatom 293 294 sigx=sigx*two_pi*third*dble(natom)/(dom*ucvol)*half/dsqrt(pi) 295 296 do icor=1,nphicor 297 do iom=1,mom 298 do iatom=1,natom 299 sigx_av(iom,icor) =sigx_av(iom,icor)+sigx(iatom,iom,icor)/dble(natom) 300 end do 301 end do 302 end do 303 304 if (open_file(trim(filnam_out)//'_sigX', msg, newunit=sigx_unt, form='formatted', action="write") /= 0) then 305 MSG_ERROR(msg) 306 end if 307 do iom=1,mom 308 write(sigx_unt,'(9(1x,e14.8))') & 309 & ((-energy_cor(icor)+oml1(iom)+omin),sigx_av(iom,icor),sigx(atnbr,iom,icor),icor=1,nphicor) 310 end do 311 close(sigx_unt) 312 313 call WffClose(wff2,ierr) 314 315 ABI_DEALLOCATE(sigx) 316 ABI_DEALLOCATE(sigx_av) 317 ABI_DEALLOCATE(ncor) 318 ABI_DEALLOCATE(lcor) 319 ABI_DEALLOCATE(energy_cor) 320 ABI_DEALLOCATE(eigen0) 321 ABI_DEALLOCATE(nband) 322 ABI_DEALLOCATE(oml1) 323 ABI_DEALLOCATE(occ) 324 ABI_DEALLOCATE(wtk) 325 326 call hdr_free(hdr) 327 328 end subroutine conducti_paw_core