TABLE OF CONTENTS
ABINIT/linear_optics_paw [ Functions ]
NAME
linear_optics_paw
FUNCTION
This program computes the elements of the optical frequency dependent linear susceptiblity using matrix elements <-i Nabla> obtained from a PAW ground state calculation. It uses formula 17 from Gadoc et al, Phys. Rev. B 73, 045112 (2006) together with a scissors correction. It uses a Kramers-Kronig transform to compute the real part from the imaginary part, and it will work on all types of unit cells. It outputs all tensor elements of both the real and imaginary parts.
COPYRIGHT
Copyright (C) 2002-2018 ABINIT group (VR, 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
filnam: base of file names to read data from mpi_enreg: mpi set up variable, not used in this code
OUTPUT
_real and _imag output files
NOTES
This routine is not tested
PARENTS
conducti
CHILDREN
destroy_mpi_enreg,hdr_free,hdr_io,hdr_read_from_fname,initmpi_seq kramerskronig,matrginv,metric,wffopen
SOURCE
41 #if defined HAVE_CONFIG_H 42 #include "config.h" 43 #endif 44 45 #include "abi_common.h" 46 47 subroutine linear_optics_paw(filnam,filnam_out,mpi_enreg_seq) 48 49 50 use defs_basis 51 use defs_abitypes 52 use m_xmpi 53 use m_errors 54 use m_wffile 55 use m_profiling_abi 56 use m_hdr 57 58 use m_geometry, only : metric 59 use m_io_tools, only : open_file, get_unit 60 use m_abilasi, only : matrginv 61 use m_mpinfo, only : destroy_mpi_enreg, nullify_mpi_enreg 62 use m_numeric_tools, only : kramerskronig 63 64 !This section has been created automatically by the script Abilint (TD). 65 !Do not modify the following lines by hand. 66 #undef ABI_FUNC 67 #define ABI_FUNC 'linear_optics_paw' 68 use interfaces_51_manage_mpi 69 !End of the abilint section 70 71 implicit none 72 73 !Arguments ----------------------------------- 74 !scalars 75 character(len=fnlen),intent(in) :: filnam,filnam_out 76 type(MPI_type),intent(inout) :: mpi_enreg_seq 77 78 !Local variables------------------------------- 79 integer,parameter :: master=0 80 integer :: iomode,bantot,bdtot_index,fform1,headform 81 integer :: iband,ierr,ii,ikpt,iom,iout,isppol,isym,jband,jj,me,mband 82 integer :: method,mom,nband_k,nkpt,nspinor,nsppol,nsym,occopt,only_check 83 integer :: rdwr,spaceComm,inpunt,reunt,imunt,wfunt 84 integer,allocatable :: nband(:),symrel(:,:,:) 85 real(dp) :: del,dom,fij,gdelta,omin,omax,paijpbij(2),mbpt_sciss,wij,ucvol 86 real(dp) :: diffwp, diffwm 87 real(dp) :: e2rot(3,3),gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3),rprimdinv(3,3),symd(3,3),symdinv(3,3) 88 real(dp),allocatable :: e1(:,:,:),e2(:,:,:,:),epsilon_tot(:,:,:,:),eigen0(:),eig0_k(:) 89 real(dp),allocatable :: kpts(:,:),occ(:),occ_k(:),oml1(:),wtk(:) 90 complex,allocatable :: eps_work(:) 91 character(len=fnlen) :: filnam1,filnam_gen 92 character(len=500) :: msg 93 type(hdr_type) :: hdr 94 type(wffile_type) :: wff1 95 !arrays 96 real(dp),allocatable :: psinablapsi(:,:,:,:) 97 98 ! ********************************************************************************* 99 100 DBG_ENTER("COLL") 101 102 !* Fake MPI_type for the sequential part. 103 !This routine should not be parallelized as communicating gbig and other 104 !tables takes more time than recalculating them in sequential. 105 call initmpi_seq(MPI_enreg_seq) 106 107 !write(std_out,'(a)')' Give the name of the output file ...' 108 !read(std_in, '(a)') filnam_out 109 !write(std_out,'(a)')' The name of the output file is :',filnam_out 110 111 !Read data file 112 if (open_file(filnam,msg,newunit=inpunt,form='formatted') /= 0 ) then 113 MSG_ERROR(msg) 114 end if 115 116 rewind(inpunt) 117 read(inpunt,*) 118 read(inpunt,'(a)')filnam_gen ! generic name for the files 119 filnam1=trim(filnam_gen)//'_OPT' ! nabla matrix elements file 120 121 !Open the Wavefunction and optic files 122 !These default values are typical of sequential use 123 iomode=IO_MODE_FORTRAN ; spaceComm=xmpi_comm_self; me=0 124 125 ! Read the header of the optic files 126 call hdr_read_from_fname(hdr, filnam1, fform1, spaceComm) 127 call hdr_free(hdr) 128 if (fform1 /= 610) then 129 MSG_ERROR("Abinit8 requires an OPT file with fform = 610") 130 end if 131 132 !Open the conducti optic files 133 wfunt = get_unit() 134 call WffOpen(iomode,spaceComm,filnam1,ierr,wff1,master,me,wfunt) 135 136 !Read the header from Ground state file 137 rdwr=1 138 call hdr_io(fform1,hdr,rdwr,wff1) 139 140 !Extract info from the header 141 headform=hdr%headform 142 bantot=hdr%bantot 143 nkpt=hdr%nkpt 144 ABI_ALLOCATE(kpts,(3,nkpt)) 145 ABI_ALLOCATE(wtk,(nkpt)) 146 kpts(:,:)=hdr%kptns(:,:) 147 wtk(:)=hdr%wtk(:) 148 nspinor=hdr%nspinor 149 nsppol=hdr%nsppol 150 occopt=hdr%occopt 151 rprimd(:,:)=hdr%rprimd(:,:) 152 rprimdinv(:,:) = rprimd(:,:) 153 call matrginv(rprimdinv,3,3) ! need the inverse of rprimd to symmetrize the tensors 154 ABI_ALLOCATE(nband,(nkpt*nsppol)) 155 ABI_ALLOCATE(occ,(bantot)) 156 occ(1:bantot)=hdr%occ(1:bantot) 157 nband(1:nkpt*nsppol)=hdr%nband(1:nkpt*nsppol) 158 nsym=hdr%nsym 159 ABI_ALLOCATE(symrel,(3,3,nsym)) 160 symrel(:,:,:)=hdr%symrel(:,:,:) 161 162 !Get mband, as the maximum value of nband(nkpt) 163 mband=maxval(nband(:)) 164 165 !get ucvol etc. 166 iout = -1 167 call metric(gmet,gprimd,iout,rmet,rprimd,ucvol) 168 169 write(std_out,*) 170 write(std_out,'(a,3f10.5,a)' )' rprimd(bohr) =',rprimd(1:3,1) 171 write(std_out,'(a,3f10.5,a)' )' ',rprimd(1:3,2) 172 write(std_out,'(a,3f10.5,a)' )' ',rprimd(1:3,3) 173 write(std_out,*) 174 write(std_out,'(a,3f10.5,a)' )' rprimdinv =',rprimdinv(1:3,1) 175 write(std_out,'(a,3f10.5,a)' )' ',rprimdinv(1:3,2) 176 write(std_out,'(a,3f10.5,a)' )' ',rprimdinv(1:3,3) 177 write(std_out,'(a,2i8)') ' nkpt,mband =',nkpt,mband 178 179 !get eigen0 180 ABI_ALLOCATE(eigen0,(mband*nkpt*nsppol)) 181 read(wfunt)(eigen0(iband),iband=1,mband*nkpt*nsppol) 182 183 read(inpunt,*)mbpt_sciss 184 read(inpunt,*)dom,omin,omax,mom 185 close(inpunt) 186 187 ABI_ALLOCATE(oml1,(mom)) 188 ABI_ALLOCATE(e1,(3,3,mom)) 189 ABI_ALLOCATE(e2,(2,3,3,mom)) 190 ABI_ALLOCATE(epsilon_tot,(2,3,3,mom)) 191 ABI_ALLOCATE(eps_work,(mom)) 192 del=(omax-omin)/(mom-1) 193 do iom=1,mom 194 oml1(iom)=omin+dble(iom-1)*del 195 end do 196 write(std_out,'(a,i8,4f10.5,a)')' npts,omin,omax,width,mbpt_sciss =',mom,omin,omax,dom,mbpt_sciss,' Ha' 197 198 ABI_ALLOCATE(psinablapsi,(2,3,mband,mband)) 199 200 !loop over spin components 201 do isppol=1,nsppol 202 bdtot_index = 0 203 ! loop over k points 204 do ikpt=1,nkpt 205 ! 206 ! number of bands for this k point 207 nband_k=nband(ikpt+(isppol-1)*nkpt) 208 ABI_ALLOCATE(eig0_k,(nband_k)) 209 ABI_ALLOCATE(occ_k,(nband_k)) 210 ! eigenvalues for this k-point 211 eig0_k(:)=eigen0(1+bdtot_index:nband_k+bdtot_index) 212 ! occupation numbers for this k-point 213 occ_k(:)=occ(1+bdtot_index:nband_k+bdtot_index) 214 ! values of -i*nabla matrix elements for this k point 215 psinablapsi=zero 216 read(wfunt)((psinablapsi(1:2,1,iband,jband),iband=1,nband_k),jband=1,nband_k) 217 read(wfunt)((psinablapsi(1:2,2,iband,jband),iband=1,nband_k),jband=1,nband_k) 218 read(wfunt)((psinablapsi(1:2,3,iband,jband),iband=1,nband_k),jband=1,nband_k) 219 220 ! occupation numbers for k-point 221 occ_k(:)=occ(1+bdtot_index:nband_k+bdtot_index) 222 ! accumulate e2 for this k point, Eq. 17 from PRB 73, 045112 (2006) 223 do iband = 1, nband_k 224 do jband = 1, nband_k 225 fij = occ_k(iband) - occ_k(jband) !occ number difference 226 wij = eig0_k(iband) - eig0_k(jband) !energy difference 227 if (abs(fij) > zero) then ! only consider states of differing occupation numbers 228 do ii = 1, 3 229 do jj = 1, 3 230 paijpbij(1) = psinablapsi(1,ii,iband,jband)*psinablapsi(1,jj,iband,jband) + & 231 & psinablapsi(2,ii,iband,jband)*psinablapsi(2,jj,iband,jband) 232 paijpbij(2) = psinablapsi(2,ii,iband,jband)*psinablapsi(1,jj,iband,jband) - & 233 & psinablapsi(1,ii,iband,jband)*psinablapsi(2,jj,iband,jband) 234 do iom = 1, mom 235 ! original version 236 ! diffw = wij + mbpt_sciss - oml1(iom) ! apply scissors term here 237 ! gdelta = exp(-diffw*diffw/(4.0*dom*dom))/(2.0*dom*sqrt(pi)) ! delta fnc resolved as Gaussian 238 ! e2(1,ii,jj,iom) = e2(1,ii,jj,iom) - (4.0*pi*pi/ucvol)*wtk(ikpt)*fij*paijpbij(1)*gdelta/(oml1(iom)*oml1(iom)) 239 ! e2(2,ii,jj,iom) = e2(2,ii,jj,iom) - (4.0*pi*pi/ucvol)*wtk(ikpt)*fij*paijpbij(2)*gdelta/(oml1(iom)*oml1(iom)) 240 diffwm = wij - mbpt_sciss + oml1(iom) ! apply scissors term here 241 diffwp = wij + mbpt_sciss - oml1(iom) ! apply scissors term here 242 gdelta = exp(-diffwp*diffwp/(4.0*dom*dom))/(2.0*dom*sqrt(pi)) 243 e2(1,ii,jj,iom) = e2(1,ii,jj,iom) - (4.0*pi*pi/ucvol)*wtk(ikpt)*fij*paijpbij(1)*gdelta/(wij*wij) 244 e2(2,ii,jj,iom) = e2(2,ii,jj,iom) - (4.0*pi*pi/ucvol)*wtk(ikpt)*fij*paijpbij(2)*gdelta/(wij*wij) 245 end do ! end loop over spectral points 246 end do ! end loop over jj = 1, 3 247 end do ! end loop over ii = 1, 3 248 end if ! end selection on fij /= 0 249 end do ! end loop over jband 250 end do ! end loop over iband 251 252 ABI_DEALLOCATE(eig0_k) 253 ABI_DEALLOCATE(occ_k) 254 bdtot_index=bdtot_index+nband_k 255 end do ! end loop over k points 256 end do ! end loop over spin polarizations 257 258 !here apply nsym symrel transformations to reconstruct full tensor from IBZ part 259 epsilon_tot(:,:,:,:) = zero 260 do isym = 1, nsym 261 symd(:,:)=matmul(rprimd(:,:),matmul(symrel(:,:,isym),rprimdinv(:,:))) 262 symdinv(:,:)=symd(:,:) 263 call matrginv(symdinv,3,3) 264 do iom = 1, mom 265 e2rot(:,:)=matmul(symdinv(:,:),matmul(e2(1,:,:,iom),symd(:,:))) 266 epsilon_tot(2,:,:,iom) = epsilon_tot(2,:,:,iom)+e2rot(:,:)/nsym 267 end do 268 end do 269 270 !generate e1 from e2 via KK transforma 271 method=0 ! use naive integration ( = 1 for simpson) 272 only_check=0 ! compute real part of eps in kk routine 273 do ii = 1, 3 274 do jj = 1, 3 275 eps_work(:) = cmplx(0.0,epsilon_tot(2,ii,jj,:)) 276 call kramerskronig(mom,oml1,eps_work,method,only_check) 277 epsilon_tot(1,ii,jj,:) = real(eps_work(:)) 278 if (ii /= jj) epsilon_tot(1,ii,jj,:) = epsilon_tot(1,ii,jj,:)- 1.0 279 end do ! end loop over jj 280 end do ! end loop over ii 281 282 if (open_file(trim(filnam_out)//'_imag',msg,newunit=reunt,form='formatted') /= 0) then 283 MSG_ERROR(msg) 284 end if 285 286 if (open_file(trim(filnam_out)//'_real',msg,unit=imunt,form='formatted') /= 0) then 287 MSG_ERROR(msg) 288 end if 289 290 write(reunt,'(a12,6a13)')' # Energy/Ha ','eps_2_xx','eps_2_yy','eps_2_zz',& 291 & 'eps_2_yz','eps_2_xz','eps_2_xy' 292 write(imunt,'(a12,6a13)')' # Energy/Ha ','eps_1_xx','eps_1_yy','eps_1_zz',& 293 & 'eps_1_yz','eps_1_xz','eps_1_xy' 294 295 do iom = 1, mom 296 write(reunt,'(ES12.4,a,ES12.4,a,ES12.4,a,ES12.4,a,ES12.4,a,ES12.4,a,ES12.4)') oml1(iom),' ',& 297 & epsilon_tot(2,1,1,iom),' ',epsilon_tot(2,2,2,iom),' ',epsilon_tot(2,3,3,iom),' ',& 298 & epsilon_tot(2,2,3,iom),' ',epsilon_tot(2,1,3,iom),' ',epsilon_tot(2,1,2,iom) 299 write(imunt,'(ES12.4,a,ES12.4,a,ES12.4,a,ES12.4,a,ES12.4,a,ES12.4,a,ES12.4)') oml1(iom),' ',& 300 & epsilon_tot(1,1,1,iom),' ',epsilon_tot(1,2,2,iom),' ',epsilon_tot(1,3,3,iom),' ',& 301 & epsilon_tot(1,2,3,iom),' ',epsilon_tot(1,1,3,iom),' ',epsilon_tot(1,1,2,iom) 302 end do 303 304 close(reunt) 305 close(imunt) 306 307 ABI_DEALLOCATE(nband) 308 ABI_DEALLOCATE(oml1) 309 ABI_DEALLOCATE(e2) 310 ABI_DEALLOCATE(e1) 311 ABI_DEALLOCATE(occ) 312 ABI_DEALLOCATE(psinablapsi) 313 ABI_DEALLOCATE(eigen0) 314 ABI_DEALLOCATE(wtk) 315 ABI_DEALLOCATE(kpts) 316 317 call hdr_free(hdr) 318 call destroy_mpi_enreg(MPI_enreg_seq) 319 320 DBG_EXIT("COLL") 321 322 end subroutine linear_optics_paw