File: dielectric_bse.f90

package info (click to toggle)
elkcode 5.4.24-2
  • links: PTS, VCS
  • area: main
  • in suites: buster, sid
  • size: 12,840 kB
  • sloc: f90: 48,415; fortran: 22,457; perl: 965; makefile: 384; sh: 369; python: 105; ansic: 67
file content (132 lines) | stat: -rw-r--r-- 3,430 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132

! Copyright (C) 2010 S. Sharma, J. K. Dewhurst and E. K. U. Gross.
! This file is distributed under the terms of the GNU General Public License.
! See the file COPYING for license details.

subroutine dielectric_bse
use modmain
use modtest
implicit none
! local variables
integer a1,a2,ik1,jk1
integer i1,j1,ist1,jst1
integer iw,i,j,l
integer iostat,nmbse_
real(8) e,eji,t1,t2
complex(8) eta,zv(3),z1
character(256) fname
! allocatable arrays
real(8), allocatable :: w(:)
complex(8), allocatable :: pmat(:,:,:),sigma(:,:,:)
! initialise global variables
call init0
call init1
! read Fermi energy from a file
call readfermi
! get the eigenvalues from file
do ik1=1,nkpt
  call getevalsv(filext,ik1,vkl(:,ik1),evalsv(:,ik1))
end do
! generate the BSE state index arrays
call genidxbse
! allocate global BSE arrays
if (allocated(evalbse)) deallocate(evalbse)
allocate(evalbse(nmbse))
if (allocated(hmlbse)) deallocate(hmlbse)
allocate(hmlbse(nmbse,nmbse))
! read in the BSE eigenvectors and eigenvalues
open(50,file='EVBSE.OUT',form='UNFORMATTED',status='OLD',iostat=iostat)
if (iostat.ne.0) then
  write(*,*)
  write(*,'("Error(dielectric_bse): error opening EVBSE.OUT")')
  write(*,*)
  stop
end if
read(50) nmbse_
if (nmbse.ne.nmbse_) then
  write(*,*)
  write(*,'("Error(dielectric_bse): differing nmbse")')
  write(*,'(" current   : ",I6)') nmbse
  write(*,'(" EVBSE.OUT : ",I6)') nmbse_
  stop
end if
read(50) evalbse
read(50) hmlbse
close(50)
! allocate local arrays
allocate(w(nwplot))
allocate(pmat(nstsv,nstsv,3))
allocate(sigma(3,3,nwplot))
! set up the frequency grid (starting from zero)
t1=wplot(2)/dble(nwplot)
do iw=1,nwplot
  w(iw)=t1*dble(iw-1)
end do
! i divided by the complex relaxation time
eta=cmplx(0.d0,swidth,8)
sigma(:,:,:)=0.d0
do a2=1,nmbse
  e=evalbse(a2)
  zv(:)=0.d0
! loop over non-reduced k-points
  do ik1=1,nkptnr
! equivalent reduced k-point
    jk1=ivkik(ivk(1,ik1),ivk(2,ik1),ivk(3,ik1))
! read the momentum matrix elements from file
    call getpmat(.false.,vkl(:,ik1),pmat)
    do i1=1,nvbse
      ist1=istbse(i1,ik1)
      do j1=1,ncbse
        jst1=jstbse(j1,ik1)
        a1=ijkbse(i1,j1,ik1)
        eji=evalsv(jst1,jk1)-evalsv(ist1,jk1)
        z1=(e/eji)*hmlbse(a1,a2)
        zv(:)=zv(:)+z1*pmat(ist1,jst1,:)
      end do
    end do
  end do
  if (abs(e).gt.1.d-8) then
    do i=1,3
      do j=1,3
        z1=zv(i)*conjg(zv(j))/e
        sigma(i,j,:)=sigma(i,j,:)+z1/(w(:)-e+eta)+conjg(z1)/(w(:)+e+eta)
      end do
    end do
  end if
end do
z1=zi*occmax*wkptnr/omega
sigma(:,:,:)=z1*sigma(:,:,:)
! loop over tensor components
do l=1,noptcomp
  i=optcomp(1,l)
  j=optcomp(2,l)
  t1=0.d0
  if (i.eq.j) t1=1.d0
  write(fname,'("EPSILON_BSE_",2I1,".OUT")') i,j
  open(50,file=trim(fname),form='FORMATTED')
  do iw=1,nwplot
    t2=t1-fourpi*aimag(sigma(i,j,iw)/(w(iw)+eta))
    write(50,'(2G18.10)') w(iw),t2
  end do
  write(50,*)
  do iw=1,nwplot
    t2=fourpi*dble(sigma(i,j,iw)/(w(iw)+eta))
    write(50,'(2G18.10)') w(iw),t2
  end do
  close(50)
end do
write(*,*)
write(*,'("Info(dielectric_bse):")')
write(*,'(" dielectric tensor written to EPSILON_BSE_ij.OUT")')
write(*,'(" for components")')
do l=1,noptcomp
  write(*,'("  i = ",I1,", j = ",I1)') optcomp(1:2,l)
end do
! write sigma to test file
call writetest(187,'BSE optical conductivity',nv=nwplot,tol=1.d-3,zva=sigma)
deallocate(w,pmat,sigma)
! deallocate global BSE arrays
deallocate(evalbse,hmlbse)
return
end subroutine