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
|
! 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 hmldbsek(ik2)
use modmain
use modomp
implicit none
! arguments
integer, intent(in) :: ik2
! local variables
integer ik1,ist1,ist2,jst1,jst2
integer i1,i2,j1,j2,a1,a2,b1,b2
integer iv(3),iq,ig,jg,nthd
real(8) vl(3),vc(3),t0,t1,t2
complex(8) z1
! automatic arrays
integer ngp(nspnfv)
! allocatable arrays
integer, allocatable :: igpig(:,:)
real(8), allocatable :: vgqc(:,:),gqc(:),gclgq(:),jlgqr(:,:,:)
complex(8), allocatable :: ylmgq(:,:),sfacgq(:,:)
complex(4), allocatable :: wfmt1(:,:,:,:),wfir1(:,:,:)
complex(4), allocatable :: wfmt2(:,:,:,:),wfir2(:,:,:)
complex(4), allocatable :: crhomt(:,:),crhoir(:)
complex(8), allocatable :: zvv(:,:,:),zcc(:,:,:)
complex(8), allocatable :: zvc(:,:,:),zcv(:,:,:)
complex(8), allocatable :: epsi(:,:,:)
allocate(igpig(ngkmax,nspnfv))
allocate(vgqc(3,ngrf),gqc(ngrf),gclgq(ngrf))
allocate(jlgqr(njcmax,nspecies,ngrf))
allocate(ylmgq(lmmaxo,ngrf),sfacgq(ngrf,natmtot))
allocate(wfmt1(npcmtmax,natmtot,nspinor,nstsv),wfir1(ngtc,nspinor,nstsv))
allocate(wfmt2(npcmtmax,natmtot,nspinor,nstsv),wfir2(ngtc,nspinor,nstsv))
allocate(zvv(ngrf,nvbse,nvbse),zcc(ngrf,ncbse,ncbse))
allocate(epsi(ngrf,ngrf,nwrf))
if (bsefull) then
allocate(zvc(ngrf,nvbse,ncbse))
allocate(zcv(ngrf,ncbse,nvbse))
end if
! generate the wavefunctions for all states of k-point ik2
call genwfsvp_sp(.false.,.false.,nstsv,[0],ngdgc,igfc,vkl(:,ik2),ngp,igpig, &
wfmt2,ngtc,wfir2)
! begin loop over ik1
do ik1=1,nkptnr
! generate the wavefunctions for all states of k-point ik1
call genwfsvp_sp(.false.,.false.,nstsv,[0],ngdgc,igfc,vkl(:,ik1),ngp,igpig, &
wfmt1,ngtc,wfir1)
! determine equivalent q-vector in first Brillouin zone
iv(:)=ivk(:,ik1)-ivk(:,ik2)
iv(:)=modulo(iv(:),ngridk(:))
iq=ivqiq(iv(1),iv(2),iv(3))
! q-vector in lattice and Cartesian coordinates
vl(:)=vkl(:,ik1)-vkl(:,ik2)
vc(:)=vkc(:,ik1)-vkc(:,ik2)
! generate the G+q-vectors and related functions
call gengqf(ngrf,vc,vgqc,gqc,jlgqr,ylmgq,sfacgq)
! generate the regularised Coulomb Green's function in G+q-space
call gengclgq(.true.,iq,ngrf,gqc,gclgq)
! symmetrise the Coulomb Green's function
gclgq(:)=sqrt(gclgq(:))
! compute the <v|exp(-i(G+q).r)|v'> matrix elements
call holdthd(nvbse,nthd)
!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP PRIVATE(crhomt,crhoir,ist1,ist2,i2) &
!$OMP NUM_THREADS(nthd)
allocate(crhomt(npcmtmax,natmtot),crhoir(ngtc))
!$OMP DO SCHEDULE(DYNAMIC)
do i1=1,nvbse
ist1=istbse(i1,ik1)
do i2=1,nvbse
ist2=istbse(i2,ik2)
call gencrho(.true.,.true.,ngtc,wfmt2(:,:,:,ist2),wfir2(:,:,ist2), &
wfmt1(:,:,:,ist1),wfir1(:,:,ist1),crhomt,crhoir)
call zftcf(ngrf,jlgqr,ylmgq,ngrf,sfacgq,crhomt,crhoir,zvv(:,i1,i2))
end do
end do
!$OMP END DO
deallocate(crhomt,crhoir)
!$OMP END PARALLEL
call freethd(nthd)
! compute the <c|exp(-i(G+q).r)|c'> matrix elements
call holdthd(ncbse,nthd)
!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP PRIVATE(crhomt,crhoir,jst1,jst2,j2) &
!$OMP NUM_THREADS(nthd)
allocate(crhomt(npcmtmax,natmtot),crhoir(ngtc))
!$OMP DO SCHEDULE(DYNAMIC)
do j1=1,ncbse
jst1=jstbse(j1,ik1)
do j2=1,ncbse
jst2=jstbse(j2,ik2)
call gencrho(.true.,.true.,ngtc,wfmt2(:,:,:,jst2),wfir2(:,:,jst2), &
wfmt1(:,:,:,jst1),wfir1(:,:,jst1),crhomt,crhoir)
call zftcf(ngrf,jlgqr,ylmgq,ngrf,sfacgq,crhomt,crhoir,zcc(:,j1,j2))
end do
end do
!$OMP END DO
deallocate(crhomt,crhoir)
!$OMP END PARALLEL
call freethd(nthd)
! matrix elements for full BSE kernel if required
if (bsefull) then
! compute the <v|exp(-i(G+q).r)|c'> matrix elements
call holdthd(nvbse,nthd)
!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP PRIVATE(crhomt,crhoir,ist1,jst2,j2) &
!$OMP NUM_THREADS(nthd)
allocate(crhomt(npcmtmax,natmtot),crhoir(ngtc))
!$OMP DO SCHEDULE(DYNAMIC)
do i1=1,nvbse
ist1=istbse(i1,ik1)
do j2=1,ncbse
jst2=jstbse(j2,ik2)
call gencrho(.true.,.true.,ngtc,wfmt2(:,:,:,jst2),wfir2(:,:,jst2), &
wfmt1(:,:,:,ist1),wfir1(:,:,ist1),crhomt,crhoir)
call zftcf(ngrf,jlgqr,ylmgq,ngrf,sfacgq,crhomt,crhoir,zvc(:,i1,j2))
end do
end do
!$OMP END DO
deallocate(crhomt,crhoir)
!$OMP END PARALLEL
call freethd(nthd)
! compute the <c|exp(-i(G+q).r)|v'> matrix elements
call holdthd(ncbse,nthd)
!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP PRIVATE(crhomt,crhoir,jst1,ist2,i2) &
!$OMP NUM_THREADS(nthd)
allocate(crhomt(npcmtmax,natmtot),crhoir(ngtc))
!$OMP DO SCHEDULE(DYNAMIC)
do j1=1,ncbse
jst1=jstbse(j1,ik1)
do i2=1,nvbse
ist2=istbse(i2,ik2)
call gencrho(.true.,.true.,ngtc,wfmt2(:,:,:,ist2),wfir2(:,:,ist2), &
wfmt1(:,:,:,jst1),wfir1(:,:,jst1),crhomt,crhoir)
call zftcf(ngrf,jlgqr,ylmgq,ngrf,sfacgq,crhomt,crhoir,zcv(:,j1,i2))
end do
end do
!$OMP END DO
deallocate(crhomt,crhoir)
!$OMP END PARALLEL
call freethd(nthd)
end if
! get RPA inverse dielectric function from file
! this is the symmetric version: ϵ⁻¹ = 1 - v¹⸍² χ₀ v¹⸍²
call getcfgq('EPSINV.OUT',vl,ngrf,nwrf,epsi)
t0=wkptnr*omega
do i1=1,nvbse
do j1=1,ncbse
a1=ijkbse(i1,j1,ik1)
do i2=1,nvbse
do j2=1,ncbse
a2=ijkbse(i2,j2,ik2)
z1=0.d0
do ig=1,ngrf
t1=t0*gclgq(ig)
do jg=1,ngrf
t2=t1*gclgq(jg)
z1=z1+t2*epsi(ig,jg,1)*conjg(zcc(ig,j1,j2))*zvv(jg,i1,i2)
end do
end do
hmlbse(a1,a2)=hmlbse(a1,a2)-z1
! compute off-diagonal blocks if required
if (bsefull) then
b1=a1+nbbse
b2=a2+nbbse
hmlbse(b1,b2)=hmlbse(b1,b2)+conjg(z1)
z1=0.d0
do ig=1,ngrf
t1=t0*gclgq(ig)
do jg=1,ngrf
t2=t1*gclgq(jg)
z1=z1+t2*epsi(ig,jg,1)*conjg(zcv(ig,j1,i2))*zvc(jg,i1,j2)
end do
end do
hmlbse(a1,b2)=hmlbse(a1,b2)-z1
hmlbse(b1,a2)=hmlbse(b1,a2)+conjg(z1)
end if
! end loop over i2 and j2
end do
end do
! end loop over i1 and j1
end do
end do
! end loop over ik1
end do
deallocate(igpig,vgqc,gqc,gclgq,jlgqr)
deallocate(ylmgq,sfacgq)
deallocate(wfmt1,wfmt2,wfir1,wfir2)
deallocate(zvv,zcc,epsi)
if (bsefull) deallocate(zvc,zcv)
end subroutine
|
