TABLE OF CONTENTS
ABINIT/exc_den [ Functions ]
NAME
exc_den
FUNCTION
This routines calculates the electron-hole excited state density.
COPYRIGHT
Copyright (C) 1992-2009 EXC group (L.Reining, V.Olevano, F.Sottile, S.Albrecht, G.Onida) Copyright (C) 2009-2018 ABINIT group (L.Reining, V.Olevano, F.Sottile, S.Albrecht, G.Onida, M.Giantomassi) 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
filbseig=Name of the file containing the excitonic eigenvectors and eigenvalues.
OUTPUT
PARENTS
bethe_salpeter
CHILDREN
wfd_get_ur,wrtout
SOURCE
30 #if defined HAVE_CONFIG_H 31 #include "config.h" 32 #endif 33 34 #include "abi_common.h" 35 36 subroutine exc_den(BSp,BS_files,ngfft,nfftot,Kmesh,ktabr,Wfd) 37 38 use defs_basis 39 use m_profiling_abi 40 use m_bs_defs 41 use m_errors 42 43 use m_io_tools, only : open_file 44 use m_bz_mesh, only : kmesh_t 45 use m_wfd, only : wfd_t, wfd_get_ur 46 47 !This section has been created automatically by the script Abilint (TD). 48 !Do not modify the following lines by hand. 49 #undef ABI_FUNC 50 #define ABI_FUNC 'exc_den' 51 use interfaces_14_hidewrite 52 !End of the abilint section 53 54 implicit none 55 56 !Arguments ------------------------------------ 57 !scalars 58 integer,intent(in) :: nfftot 59 type(excparam),intent(in) :: BSp 60 type(excfiles),intent(in) :: BS_files 61 type(kmesh_t),intent(in) :: Kmesh 62 type(wfd_t),intent(inout) :: Wfd 63 !arrays 64 integer,intent(in) :: ngfft(18) 65 integer,intent(in) :: ktabr(nfftot,BSp%nkbz) 66 67 !Local variables ------------------------------ 68 !scalars 69 integer :: nfft1,nfft2,nfft3,spin,ierr 70 integer :: it,itp,ik_ibz,ikp_bz,ik_bz,band,iv,ivp,ic,icp,ir,i1,i2,i3,iik 71 integer :: state,nstates,min_idx,eig_unt,den_unt,sden_unt,hsize,homo 72 real(dp) :: min_ene,cost 73 character(len=500) :: msg 74 character(len=fnlen) :: filbseig 75 !arrays 76 real(dp) :: n0(nfftot),rho_eh(nfftot),nexc(nfftot) 77 real(dp),allocatable :: exc_ene(:) 78 complex(dpc) :: rhor_h(nfftot),rhor_e(nfftot) 79 complex(gwpc),allocatable :: wfr(:,:,:), wfrk(:,:) 80 complex(dpc),allocatable :: exc_ene_cplx(:),exc_vec(:) 81 82 !************************************************************************ 83 84 ABI_CHECK(Wfd%nsppol==1,"nsppol==2 not coded") 85 86 if (BS_files%in_eig /= BSE_NOFILE) then 87 filbseig = BS_files%in_eig 88 else 89 filbseig = BS_files%out_eig 90 end if 91 92 call wrtout(std_out,' Calculating electron-hole, excited state density',"COLL") 93 call wrtout(std_out," Reading eigenstates from: "//TRIM(filbseig),"COLL") 94 MSG_ERROR("Not tested") 95 96 ! Prepare FFT tables to have u(r) on the ngfft_osc mesh. 97 !if ( ANY(ngfftf(1:3) /= Wfd%ngfft(1:3)) ) then 98 ! call wfd_change_ngfft(Wfd,Cryst,Psps,ngfftf) 99 !end if 100 101 nfft1 = ngfft(1) 102 nfft2 = ngfft(2) 103 nfft3 = ngfft(3) 104 ABI_CHECK(nfftot==PRODUCT(ngfft(1:3)),"Mismatch in FFT size") 105 106 !allocate and load wavefunctions in real space 107 ABI_STAT_MALLOC(wfr,(nfftot*Wfd%nspinor,BSp%nbnds,Wfd%nkibz), ierr) 108 ABI_CHECK(ierr==0, "out of memory: exc_den, wfr") 109 110 spin=1 ! SPIN support is missing. 111 112 do ik_ibz=1,Wfd%nkibz 113 do band=1,BSp%nbnds 114 call wfd_get_ur(Wfd,band,ik_ibz,spin,wfr(:,band,ik_ibz)) 115 end do 116 end do 117 118 if (open_file(filbseig,msg,newunit=eig_unt,form="unformatted", status="old", action="read") /= 0) then 119 MSG_ERROR(msg) 120 end if 121 122 read(eig_unt) hsize,nstates 123 124 if (nstates/=Bsp%nreh(1)) then 125 write(std_out,*) 'not resonant calculation' 126 close(eig_unt) 127 RETURN 128 end if 129 130 ABI_MALLOC(exc_ene_cplx,(nstates)) 131 ABI_MALLOC(exc_ene,(nstates)) 132 133 read(eig_unt) exc_ene_cplx(1:nstates) 134 135 ! Small imaginary part is always neglected. 136 exc_ene(:) = exc_ene_cplx(:) 137 ABI_FREE(exc_ene_cplx) 138 ! 139 ! Find the lowest non-negative eigenvalue. 140 min_ene = greatest_real 141 do state=1,nstates 142 if (exc_ene(state) < zero) cycle 143 if (exc_ene(state) < min_ene) then 144 min_ene = exc_ene(state) 145 min_idx = state 146 end if 147 end do 148 ABI_FREE(exc_ene) 149 150 write(std_out,*)' Lowest eigenvalue ', min_idx, min_ene*Ha_eV 151 ! 152 ! skip other eigenvectors 153 do state=1,min_idx-1 154 read(eig_unt) 155 end do 156 ! 157 ! read "lowest" eigenvector 158 ABI_MALLOC(exc_vec,(hsize)) 159 read(eig_unt) exc_vec 160 161 close(eig_unt) 162 163 ABI_STAT_MALLOC(wfrk,(nfftot,BSp%nbnds), ierr) 164 ABI_CHECK(ierr==0, 'out of memory: exc_den, wfrk') 165 166 !calculate ground state density 167 n0(:) = zero 168 spin = 1 169 homo = Bsp%homo_spin(spin) 170 do ik_bz = 1, BSp%nkbz 171 ik_ibz = Kmesh%tab(ik_bz) 172 iik = (3-Kmesh%tabi(ik_bz))/2 173 174 if (iik==1) then 175 do ir=1,nfftot 176 wfrk(ir,1:homo) = wfr(ktabr(ir,ik_bz),1:homo,ik_ibz) 177 end do 178 else 179 do ir=1,nfftot 180 wfrk(ir,1:homo) = conjg(wfr(ktabr(ir,ik_bz),1:homo,ik_ibz)) 181 end do 182 end if 183 184 do iv=1,homo 185 n0(:) = n0(:) + conjg(wfrk(:,band)) * wfrk(:,band) 186 end do 187 end do 188 ! 189 ! calculate electron and hole density 190 rhor_h=czero 191 rhor_e=czero 192 ABI_CHECK(BSp%nsppol==1,"nsppol=2 not coded") 193 194 do it=1,Bsp%nreh(1) 195 ik_bz = Bsp%Trans(it,1)%k 196 iv = Bsp%Trans(it,1)%v 197 ic = Bsp%Trans(it,1)%c 198 ik_ibz = Kmesh%tab(ik_bz) 199 iik = (3-Kmesh%tabi(ik_bz))/2 200 201 if (iik==1) then 202 do ir = 1, nfftot 203 wfrk(ir,:) = wfr(ktabr(ir,ik_bz),:,ik_ibz) 204 end do 205 else 206 do ir = 1, nfftot 207 wfrk(ir,:) = conjg(wfr(ktabr(ir,ik_bz),:,ik_ibz)) 208 end do 209 end if 210 ! 211 do itp=1,Bsp%nreh(1) 212 ikp_bz = Bsp%Trans(itp,1)%k 213 if (ikp_bz/=ik_bz) CYCLE 214 icp = Bsp%Trans(itp,1)%c 215 ivp = Bsp%Trans(itp,1)%v 216 if(icp==ic) then 217 rhor_h = rhor_h - conjg(exc_vec(it)) * exc_vec(itp) * wfrk(:,iv) * conjg(wfrk(:,ivp)) 218 end if 219 if(ivp==iv) then 220 rhor_e = rhor_e + conjg(exc_vec(it)) * exc_vec(itp) * conjg(wfrk(:,ic)) * wfrk(:,icp) 221 end if 222 end do 223 end do 224 225 ABI_FREE(exc_vec) 226 ! 227 ! calculate excited state density minus ground state density 228 ! n* - n0 = rhor_e + rhor_h 229 rho_eh(:) = rhor_e + rhor_h 230 ! 231 !calculate excited state density 232 ! n* = n0 + rhor_e + rhor_h 233 nexc(:) = n0(:) + rho_eh(:) 234 235 if (open_file('out.den',msg,newunit=den_unt) /= 0) then 236 MSG_ERROR(msg) 237 end if 238 239 ! here conversion to cartesian through a1, a2, a3 240 cost = zero 241 do i1=0,nfft1-1 242 do i2=0,nfft2-1 243 do i3=0,nfft3-1 244 ir=1 + i1 + i2*nfft1 + i3*nfft1*nfft2 245 write(den_unt,'(3i3,2x,5e11.4)')i1,i2,i3,n0(ir),nexc(ir),rho_eh(ir),real(rhor_e(ir)),real(rhor_h(ir)) 246 cost = cost + n0(ir) 247 end do 248 end do 249 end do 250 251 write(std_out,*) 'density normalization constant ', cost 252 close(den_unt) 253 254 if (open_file('out.sden',msg,newunit=sden_unt) /= 0) then 255 MSG_ERROR(msg) 256 end if 257 258 ! we are looking for the plane between (100) (111) 259 ! so it's the place where (111) [ (100) x v ] = 0 (mixed product) 260 ! finally v2 - v3 = 0 261 do i2=0, nfft2-1 262 do i3=0, nfft3-1 263 if(i2==i3) then 264 write(std_out,*) i2, i3 265 write(sden_unt,*) (rho_eh(1+i1+i2*nfft1+i3*nfft1*nfft2),i1=0,nfft1-1) 266 end if 267 end do 268 end do 269 270 close(sden_unt) 271 272 end subroutine exc_den