TABLE OF CONTENTS


ABINIT/linear_optics_paw [ Functions ]

[ Top ] [ 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