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
|
! 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 genvchi0(t3hw,ik,lock,vqpl,gclgq,jlgqr,ylmgq,sfacgq,nm,vchi0)
use modmain
use modomp
implicit none
! local variables
logical, intent(in) :: t3hw
integer, intent(in) :: ik
integer(omp_lock_kind), intent(inout) :: lock(nwrf)
real(8), intent(in) :: vqpl(3),gclgq(ngrf),jlgqr(njcmax,nspecies,ngrf)
complex(8), intent(in) :: ylmgq(lmmaxo,ngrf),sfacgq(ngrf,natmtot)
integer, intent(in) :: nm
complex(4), intent(inout) :: vchi0(nm,nm,nwrf)
! local variables
logical tq0
integer isym,jk,jkq,iw,nthd
integer nst,nstq,ist,jst,kst,lst
integer nm2,ig0,ig1,ig,jg,i,j
real(8) vkql(3),ei,ej,eij,t1
complex(8) a(3,3),z1
complex(4) c1
! automatic arrays
integer idx(nstsv),idxq(nstsv)
integer ngp(nspnfv),ngpq(nspnfv)
! allocatable arrays
integer, allocatable :: igpig(:,:),igpqig(:,:)
complex(4), allocatable :: wfmt(:,:,:,:),wfir(:,:,:)
complex(4), allocatable :: wfmtq(:,:,:,:),wfirq(:,:,:)
complex(4), allocatable :: crhomt(:,:),crhoir(:),cw(:),b(:,:)
complex(8), allocatable :: zrhoig(:),pmat(:,:,:)
! check if q = 0
tq0=.false.
if (sum(abs(vqpl(:))) < epslat) tq0=.true.
! k+q-vector in lattice coordinates
vkql(:)=vkl(:,ik)+vqpl(:)
! equivalent reduced k-points for k and k+q
jk=ivkik(ivk(1,ik),ivk(2,ik),ivk(3,ik))
call findkpt(vkql,isym,jkq)
! count and index states at k and k+q in energy window
nst=0
do ist=1,nstsv
if (abs(evalsv(ist,jk)-efermi) > emaxrf) cycle
nst=nst+1
idx(nst)=ist
end do
nstq=0
do ist=1,nstsv
if (abs(evalsv(ist,jkq)-efermi) > emaxrf) cycle
nstq=nstq+1
idxq(nstq)=ist
end do
! generate the wavefunctions for all states at k and k+q in energy window
allocate(igpig(ngkmax,nspnfv))
allocate(wfmt(npcmtmax,natmtot,nspinor,nst),wfir(ngtc,nspinor,nst))
call genwfsvp_sp(.false.,.false.,nst,idx,ngdgc,igfc,vkl(:,ik),ngp,igpig,wfmt, &
ngtc,wfir)
deallocate(igpig)
allocate(igpqig(ngkmax,nspnfv))
allocate(wfmtq(npcmtmax,natmtot,nspinor,nstq),wfirq(ngtc,nspinor,nstq))
call genwfsvp_sp(.false.,.false.,nstq,idxq,ngdgc,igfc,vkql,ngpq,igpqig,wfmtq, &
ngtc,wfirq)
deallocate(igpqig)
! read the momentum matrix elements from file for q = 0
if (tq0) then
allocate(pmat(nstsv,nstsv,3))
call getpmat(vkl(:,ik),pmat)
! divide by unit cell volume
t1=1.d0/omega
pmat(1:nstsv,1:nstsv,1:3)=t1*pmat(1:nstsv,1:nstsv,1:3)
end if
nm2=nm**2
call holdthd(nst,nthd)
!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP PRIVATE(crhomt,crhoir,zrhoig,cw,b) &
!$OMP PRIVATE(jst,kst,lst,ei,ej,eij,t1,iw) &
!$OMP PRIVATE(ig,jg,z1,c1,ig0,ig1,i,j,a) &
!$OMP NUM_THREADS(nthd)
allocate(crhomt(npcmtmax,natmtot),crhoir(ngtc))
allocate(zrhoig(ngrf),cw(nwrf))
if (tq0.and.t3hw) then
allocate(b(-1:ngrf,-1:ngrf))
else
allocate(b(ngrf,ngrf))
end if
!$OMP DO SCHEDULE(DYNAMIC)
do ist=1,nst
kst=idx(ist)
ei=evalsv(kst,jk)
do jst=1,nstq
lst=idxq(jst)
t1=wkptnr*omega*(occsv(kst,jk)-occsv(lst,jkq))
if (abs(t1) < 1.d-8) cycle
ej=evalsv(lst,jkq)
eij=ei-ej
! frequency-dependent part in response function formula for all frequencies
do iw=1,nwrf
cw(iw)=t1/(eij+wrf(iw))
end do
! compute the complex density in G+q-space
call gencrho(.true.,.true.,ngtc,wfmt(:,:,:,ist),wfir(:,:,ist), &
wfmtq(:,:,:,jst),wfirq(:,:,jst),crhomt,crhoir)
call zftcf(ngrf,jlgqr,ylmgq,ngrf,sfacgq,crhomt,crhoir,zrhoig)
! Hermitian part of body
do jg=1,ngrf
do ig=1,jg-1
b(ig,jg)=conjg(b(jg,ig))
end do
z1=gclgq(jg)*conjg(zrhoig(jg))
do ig=jg,ngrf
b(ig,jg)=gclgq(ig)*zrhoig(ig)*z1
end do
end do
! case of q = 0
if (tq0) then
if (t3hw) then
b(-1:1,-1:1)=0.e0
! calculate 3 x ngrf wings of matrix
t1=-sqrt(fourpi)/eij
do i=-1,1
z1=t1*pmat(kst,lst,i+2)
b(i,2:ngrf)=z1*conjg(zrhoig(2:ngrf))*gclgq(2:ngrf)
do j=2,ngrf
b(j,i)=conjg(b(i,j))
end do
end do
else
! use trace of 3 x 3 head of matrix
t1=sum(dble(pmat(kst,lst,1:3))**2+aimag(pmat(kst,lst,1:3))**2)/3.d0
b(1,1)=(fourpi/eij**2)*t1
! wings of matrix
t1=-sqrt(fourpi)/eij
z1=(t1/3.d0)*(pmat(kst,lst,1)+pmat(kst,lst,2)+pmat(kst,lst,3))
b(1,2:ngrf)=z1*conjg(zrhoig(2:ngrf))*gclgq(2:ngrf)
b(2:ngrf,1)=conjg(b(1,2:ngrf))
end if
end if
! add to body and wings of the response function
if (t3hw.or.(mbwgrf < 0)) then
! full matrix
do iw=1,nwrf
call omp_set_lock(lock(iw))
call caxpy(nm2,cw(iw),b,1,vchi0(1,1,iw),1)
call omp_unset_lock(lock(iw))
end do
else
! treat matrix as banded
do iw=1,nwrf
call omp_set_lock(lock(iw))
c1=cw(iw)
do jg=1,ngrf
ig0=max(jg-mbwgrf,1); ig1=min(jg+mbwgrf,ngrf)
vchi0(ig0:ig1,jg,iw)=vchi0(ig0:ig1,jg,iw)+c1*b(ig0:ig1,jg)
end do
call omp_unset_lock(lock(iw))
end do
end if
! calculate 3 x 3 head with alternative formula to improve numerical accuracy
if (tq0.and.t3hw) then
t1=-fourpi/eij
cw(1:nwrf)=cw(1:nwrf)/wrf(1:nwrf)
do j=1,3
do i=1,3
a(i,j)=t1*pmat(kst,lst,i)*conjg(pmat(kst,lst,j))
end do
end do
! add to the head of the response function
do iw=1,nwrf
call omp_set_lock(lock(iw))
vchi0(1:3,1:3,iw)=vchi0(1:3,1:3,iw)+cw(iw)*a(1:3,1:3)
call omp_unset_lock(lock(iw))
end do
end if
! end loop over jst
end do
! end loop over ist
end do
!$OMP END DO
deallocate(crhomt,crhoir,zrhoig,cw,b)
!$OMP END PARALLEL
call freethd(nthd)
deallocate(wfmt,wfir,wfmtq,wfirq)
if (tq0) deallocate(pmat)
end subroutine
|
