TABLE OF CONTENTS


ABINIT/conducti_paw_core [ Functions ]

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