-
Notifications
You must be signed in to change notification settings - Fork 3
/
get_band_unfold.f90
202 lines (169 loc) · 7.54 KB
/
get_band_unfold.f90
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
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
#include "alias.inc"
subroutine get_band_unfold(WAVEC, PINPT)
use mpi_setup
use mykind
use parameters
use print_io
use utils
use do_math
implicit none
type(incar ) :: PINPT
type(eigen ) :: WAVEC
type(poscar ) :: PGEOM_PC, PGEOM_SC
if((.not. PINPT%flag_set_unfold .and. .not. PINPT%flag_unfold)) return
if(PINPT%flag_set_unfold) PINPT%flag_unfold = .FALSE.
if(PINPT%flag_set_unfold) then
call prepare_unfold_kpoint(PINPT, PGEOM_PC, PGEOM_SC, 1, PINPT%flag_reduce)
call write_unfold_kpoint(PINPT, PGEOM_SC)
return
endif
! calculate spectral weight P_Km(k)
call get_spectral_weight_band(WAVEC, PINPT, PGEOM_PC)
! calculate spectral function A(k,E)
call get_spectral_function(WAVEC, PINPT, PGEOM_PC)
! write result
call write_spectral_function(WAVEC, PINPT, PGEOM_PC)
return
endsubroutine
! compute spectral function
subroutine get_spectral_function(WAVEC, PINPT, PGEOM_PC)
use mpi_setup
use mykind
use parameters
use print_io
use utils
use do_math
implicit none
type(incar ) :: PINPT
type(eigen ) :: WAVEC
type(poscar ) :: PGEOM_PC
integer(kind=sp) ii
integer(kind=sp) nmin, nmax, nediv
real(kind=dp) sigma, erange(PINPT%nediv), de
integer(kind=sp) ik, ie, je, is
integer(kind=sp) ourjob(nprocs), ourjob_disp(0:nprocs-1)
real(kind=dp), allocatable :: SW(:,:,:)
! NOTE:
! spectral function SW -> A(E,k) = sum_m P_Km(k) * delta(Em - ef - E)
! ef : fermi level. It is already substituted when reading Em (see wavecar.f90:Enk_read)
nmin = PINPT%ie_init
nmax = PINPT%ie_fina ; if(nmax .gt. PINPT%nband) nmax = PINPT%nband
sigma = PINPT%sigma
nediv = PINPT%nediv
erange = PINPT%init_e + eta + &
dble((/(ii,ii=0,nediv-1)/))*(PINPT%fina_e-PINPT%init_e)/dble(nediv-1)
#ifdef MPI
call MPI_BARRIER(mpi_comm_earth, mpierr)
#endif
if(allocated(WAVEC%SW)) deallocate(WAVEC%SW)
allocate(WAVEC%SW(1,nediv,PINPT%ispin,PGEOM_PC%nkpts)) ; WAVEC%SW = 0d0
allocate( SW( nediv,PINPT%ispin,PGEOM_PC%nkpts)) ; SW = 0d0
call mpi_job_distribution_chain(PGEOM_PC%nkpts, nprocs, ourjob, ourjob_disp)
do ik = sum(ourjob(1:myid))+1, sum(ourjob(1:myid+1))
do ie = nmin, nmax
do is = 1, PINPT%ispin
SW(:,is,ik) = SW(:,is,ik) + &
florentz(sigma,WAVEC%E(ie,is,PGEOM_PC%ikpt(ik))-erange)*WAVEC%SW_BAND(ie,is,ik)
enddo
enddo
enddo
#ifdef MPI
call MPI_BARRIER(mpi_comm_earth, mpierr)
call MPI_ALLREDUCE(SW, WAVEC%SW(1,:,:,:), size(SW), MPI_REAL8, MPI_SUM, mpi_comm_earth, mpierr)
#else
WAVEC%SW(1,:,:,:) = SW
#endif
return
endsubroutine
! compute spectral weight in the primitive BZ for each band.
subroutine get_spectral_weight_band(WAVEC, PINPT, PGEOM_PC)
use mpi_setup
use mykind
use parameters
use print_io
use wavecar, only: get_Gid, Cnk, CSnk, CSid
use utils
implicit none
type(incar ) :: PINPT
type(eigen ) :: WAVEC
type(poscar ) :: PGEOM_PC, PGEOM_SC
integer(kind=sp) ik, ik_W, ispinor
integer(kind=sp) is, ie
integer(kind=sp), allocatable :: ikpt(:)
integer(kind=sp) Cid(WAVEC%nplw_max)
integer(kind=sp) nplw_screen
complex(kind=dp) CSik(WAVEC%nplw_max/PINPT%ispinor,PINPT%ispinor)
complex(kind=dp) Cik(WAVEC%nplw_max/PINPT%ispinor,PINPT%ispinor)
real(kind=dp), allocatable :: SW_BAND(:,:,:) ! spectral weight
integer(kind=sp) nmin, nmax
integer(kind=sp) ourjob(nprocs), ourjob_disp(0:nprocs-1)
! Get spectral weight for k0 [see eq(8) of PRB 85, 085201 (2012)]
!
! P_{Km}(k0) = \sum_n |<Km | kn>|^2
!
! --> adopting a plane-wave expansion, [ see eq(15) of PRB 85, 085201 (2012)]
!
! P_{Km}(k0) = \sum_{g} |C_{Km}(g + k0 - K)|^2
!
! --> k - K = G (shift G vector = PGEOM_SC%G)
! --> Summation runs over g vectors of primitive cell BZ.
! Therefore, one should find g vectors from the supercell G vectors.
! In "screen_PC_G_von_SC_G" routine, it finds g vectors from supercell G vectors
! and adding shifting G vectors (PGEOM_SC%G)
!
! g + k0 - K = g + PGEOM_SC%G => GS
!
! Set unique index for GS vectors using get_Gid function
! P_{Km}(k0) represents the probability of finding a set of PC states |k0,n>
! contributing to the SC state |K,m>, or, equivalently,
! the amount of Bloch character k0 preserved in |K,m> at the same energy En = Em.
nmin = PINPT%ie_init
nmax = PINPT%ie_fina ; if(nmax .gt. PINPT%nband) nmax = PINPT%nband
call prepare_unfold_kpoint(PINPT, PGEOM_PC, PGEOM_SC, 0, .false.)
write(message,'(A,I5,A)')' --> find ',PGEOM_SC%nkpts, ' k-points along K-path' ; write_msgi
allocate(SW_BAND(PINPT%nband, PINPT%ispin, PGEOM_PC%nkpts)) ; SW_BAND = 0d0
if( allocated(WAVEC%SW_BAND)) deallocate(WAVEC%SW_BAND)
allocate(WAVEC%SW_BAND(PINPT%nband, PINPT%ispin, PGEOM_PC%nkpts)) ; WAVEC%SW_BAND = 0d0
allocate(ikpt(PGEOM_PC%nkpts)) ; ikpt = 0
if( allocated(PGEOM_PC%ikpt)) deallocate(PGEOM_PC%ikpt)
allocate(PGEOM_PC%ikpt(PGEOM_PC%nkpts)) ; PGEOM_PC%ikpt = 0
write(message,'(A)')' ' ; write_msgi
write(message,'(A)')'#- START BAND UNFOLDING: ' ; write_msgi
call mpi_job_distribution_chain(PGEOM_PC%nkpts, nprocs, ourjob, ourjob_disp)
do ik = sum(ourjob(1:myid))+1, sum(ourjob(1:myid+1))
call K_index_von_WAVECAR(ik_W, WAVEC, PGEOM_SC%kpts(:,ik))
write(message,'(A,I4,A,3F9.4,A,3F9.4,A,3I4,A,I5)')' --> Processing ',ik,'-th k-point k0: ', PGEOM_PC%kpts(:,ik), &
' K0: ',WAVEC%kpts(:,ik_W),' G0: ',PGEOM_SC%G(:,ik),' ikp= ',ik_W ;call write_log(message,12,myid)
call screen_PC_G_von_SC_G(ik_W, WAVEC, PGEOM_SC%M, PGEOM_SC%G(:,ik), PINPT%ispinor)
call get_Gid(WAVEC%GSid(:,ik_W), WAVEC%GS(:,:,ik_W), WAVEC%nplw_screen(ik_W), WAVEC%nplw_max, WAVEC%nbmax)
ikpt(ik) = ik_W
nplw_screen = WAVEC%nplw_screen(ik_W)
Cid = CSid(WAVEC, PINPT, ik_W)
do is = 1, PINPT%ispin
do ie = nmin, nmax
CSik = CSnk(Cid, WAVEC, PINPT, ie, is, ik_W) ! obtain screened Coeff
do ispinor = 1, PINPT%ispinor
SW_BAND(ie,is,ik)=SW_BAND(ie,is,ik)+dble(dot_product(CSik(1:nplw_screen,ispinor),CSik(1:nplw_screen,ispinor)))
enddo
enddo
enddo
enddo
#ifdef MPI
call MPI_BARRIER(mpi_comm_earth, mpierr)
call MPI_ALLREDUCE(SW_BAND, WAVEC%SW_BAND, size(SW_BAND), MPI_REAL8, MPI_SUM, mpi_comm_earth, mpierr)
call MPI_ALLREDUCE(ikpt, PGEOM_PC%ikpt, size(ikpt), MPI_INTEGER4, MPI_SUM, mpi_comm_earth, mpierr)
#else
WAVEC%SW_BAND = SW_BAND
PGEOM_PC%ikpt = ikpt
#endif
call write_result_spectral_weight(WAVEC, PINPT, PGEOM_PC)
return
endsubroutine
! function integ_exp(init_z, fina_z, atten_z) result(f)
! use mykind
! implicit none
! real(kind=dp) init_z, fina_z, atten_z
! complex(kind=dp) f
!
! return
! endfunction