File: genvchi0.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 (185 lines) | stat: -rw-r--r-- 5,394 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
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

! 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,scsr,vqpl,gclgq,jlgqr,ylmgq,sfacgq,nm,vchi0)
use modmain
use modomp
implicit none
! local variables
logical, intent(in) :: t3hw
integer, intent(in) :: ik
integer(8), intent(in) :: lock(nwrf)
real(8), intent(in) :: scsr,vqpl(3),gclgq(ngrf)
real(8), intent(in) :: jlgqr(njcmax,nspecies,ngrf)
complex(8), intent(in) :: ylmgq(lmmaxo,ngrf)
complex(8), intent(in) :: sfacgq(ngrf,natmtot)
integer, intent(in) :: nm
complex(8), intent(inout) :: vchi0(nm,nm,nwrf)
! local variables
logical tq0
integer isym,jk,jkq,iw
integer nst,nstq,ist,jst,kst,lst
integer nm2,ig,jg,i,j,nthd
real(8) vkql(3),eij,t1,t2
complex(8) a(3,3),z1
! automatic arrays
integer idx(nstsv),idxq(nstsv)
integer ngp(nspnfv),ngpq(nspnfv)
! allocatable arrays
integer, allocatable :: igpig(:,:),igpqig(:,:)
complex(8), allocatable :: wfmt(:,:,:,:),wfir(:,:,:)
complex(8), allocatable :: wfmtq(:,:,:,:),wfirq(:,:,:)
complex(8), allocatable :: zrhomt(:,:),zrhoir(:),zrhoig(:)
complex(8), allocatable :: pmat(:,:,:),zw(:),b(:,:)
! check if q=0
tq0=.false.
if (sum(abs(vqpl(:))).lt.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).lt.emaxrf) then
    nst=nst+1
    idx(nst)=ist
  end if
end do
nstq=0
do jst=1,nstsv
  if (abs(evalsv(jst,jkq)-efermi).lt.emaxrf) then
    nstq=nstq+1
    idxq(nstq)=jst
  end if
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(.false.,.false.,nst,idx,ngdc,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(.false.,.false.,nstq,idxq,ngdc,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(.false.,vkl(:,ik),pmat)
! divide by unit cell volume
  t1=1.d0/omega
  pmat(:,:,:)=t1*pmat(:,:,:)
end if
nm2=nm**2
call omp_hold(nst,nthd)
!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP PRIVATE(zrhomt,zrhoir,zrhoig,zw,b) &
!$OMP PRIVATE(jst,kst,lst,t1,t2,eij) &
!$OMP PRIVATE(iw,ig,jg,z1,i,j,a) &
!$OMP NUM_THREADS(nthd)
!$OMP DO
do ist=1,nst
  allocate(zrhomt(npcmtmax,natmtot),zrhoir(ngtc))
  allocate(zrhoig(ngrf),zw(nwrf))
  if (tq0.and.t3hw) then
    allocate(b(-1:ngrf,-1:ngrf))
  else
    allocate(b(ngrf,ngrf))
  end if
  kst=idx(ist)
  do jst=1,nstq
    lst=idxq(jst)
    t1=wkptnr*omega*(occsv(kst,jk)-occsv(lst,jkq))
    if (abs(t1).lt.1.d-8) cycle
    eij=evalsv(kst,jk)-evalsv(lst,jkq)
! scissor operator
    if (abs(scsr).gt.1.d-8) then
      t2=eij
      if (eij.gt.0.d0) then
        eij=eij+scsr
      else
        eij=eij-scsr
      end if
      t2=eij/t2
! scale the momentum matrix elements for q=0
      if (tq0) pmat(kst,lst,:)=t2*pmat(kst,lst,:)
    end if
! frequency-dependent part in response function formula for all frequencies
    do iw=1,nwrf
      zw(iw)=t1/(eij+wrf(iw))
    end do
! compute the complex density in G+q-space
    call genzrho(.true.,.true.,ngtc,wfmt(:,:,:,ist),wfir(:,:,ist), &
     wfmtq(:,:,:,jst),wfirq(:,:,jst),zrhomt,zrhoir)
    call zftzf(ngrf,jlgqr,ylmgq,ngrf,sfacgq,zrhomt,zrhoir,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.d0
! 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:)=z1*conjg(zrhoig(2:))*gclgq(2:)
          b(2:,i)=conjg(b(i,2:))
        end do
      else
! use trace of 3 x 3 head of matrix
        t1=sum(dble(pmat(kst,lst,:))**2+aimag(pmat(kst,lst,:))**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:)=z1*conjg(zrhoig(2:))*gclgq(2:)
        b(2:,1)=conjg(b(1,2:))
      end if
    end if
! add to body of response function
    do iw=1,nwrf
      call omp_set_lock(lock(iw))
      call zaxpy(nm2,zw(iw),b,1,vchi0(1,1,iw),1)
      call omp_unset_lock(lock(iw))
    end do
! calculate 3 x 3 head
    if (tq0.and.t3hw) then
      t1=-fourpi/eij
      zw(:)=zw(:)/wrf(:)
      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
      do iw=1,nwrf
        call omp_set_lock(lock(iw))
        vchi0(1:3,1:3,iw)=vchi0(1:3,1:3,iw)+a(1:3,1:3)*zw(iw)
        call omp_unset_lock(lock(iw))
      end do
    end if
! end loop over jst
  end do
  deallocate(zrhomt,zrhoir,zrhoig,zw,b)
! end loop over ist
end do
!$OMP END DO
!$OMP END PARALLEL
call omp_free(nthd)
deallocate(wfmt,wfir,wfmtq,wfirq)
if (tq0) deallocate(pmat)
return
end subroutine