File: genwfpw.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 (223 lines) | stat: -rw-r--r-- 6,356 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
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223

! Copyright (C) 2012 J. K. Dewhurst, S. Sharma 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 genwfpw(vpl,ngp,igpig,vgpl,vgpc,gpc,tpgpc,sfacgp,nhp,vhpc,hpc, &
 tphpc,sfachp,wfpw)
use modmain
use modpw
implicit none
! arguments
real(8), intent(in) :: vpl(3)
integer, intent(in) :: ngp(nspnfv),igpig(ngkmax,nspnfv)
real(8), intent(in) :: vgpl(3,ngkmax,nspnfv),vgpc(3,ngkmax,nspnfv)
real(8), intent(in) :: gpc(ngkmax,nspnfv),tpgpc(2,ngkmax,nspnfv)
complex(8), intent(in) :: sfacgp(ngkmax,natmtot,nspnfv)
integer, intent(in) :: nhp(nspnfv)
real(8), intent(in) :: vhpc(3,nhkmax,nspnfv),hpc(nhkmax,nspnfv)
real(8), intent(in) :: tphpc(2,nhkmax,nspnfv)
complex(8), intent(in) :: sfachp(nhkmax,natmtot,nspnfv)
complex(8), intent(out) :: wfpw(nhkmax,nspinor,nstsv)
! local variables
integer ispn0,ispn1,ispn,jspn
integer ist,is,ia,ias
integer nrc,nrci,irco,irc
integer lmax,l,m,lm
integer npci,i,igp,ihp
real(8) t0,t1,t2
complex(8) zsum1,zsum2
complex(8) z1,z2,z3,z4
! automatic arrays
integer idx(nstsv)
real(8) fr1(nrcmtmax),fr2(nrcmtmax)
complex(8) ylm(lmmaxo)
! allocatable arrays
real(8), allocatable :: jl(:,:)
complex(8), allocatable :: apwalm(:,:,:,:,:)
complex(8), allocatable :: evecfv(:,:,:),evecsv(:,:)
complex(8), allocatable :: wfmt(:,:,:,:),wfir(:,:,:)
! external functions
real(8) splint
external splint
allocate(apwalm(ngkmax,apwordmax,lmmaxapw,natmtot,nspnfv))
allocate(evecfv(nmatmax,nstfv,nspnfv),evecsv(nstsv,nstsv))
allocate(wfmt(npcmtmax,natmtot,nspinor,nstsv))
allocate(wfir(ngkmax,nspinor,nstsv))
! get the eigenvectors from file
call getevecfv(filext,0,vpl,vgpl,evecfv)
call getevecsv(filext,0,vpl,evecsv)
! find the matching coefficients
do ispn=1,nspnfv
  call match(ngp(ispn),gpc(:,ispn),tpgpc(:,:,ispn),sfacgp(:,:,ispn), &
   apwalm(:,:,:,:,ispn))
end do
! index to all states
do ist=1,nstsv
  idx(ist)=ist
end do
! calculate the second-variational wavefunctions for all states
call genwfsv(.true.,.true.,nstsv,idx,ngridg,igfft,ngp,igpig,apwalm,evecfv, &
 evecsv,wfmt,ngkmax,wfir)
deallocate(apwalm,evecfv,evecsv)
! zero the plane wave coefficients
wfpw(:,:,:)=0.d0
!---------------------------!
!     interstitial part     !
!---------------------------!
do jspn=1,nspnfv
  if (spinsprl) then
    ispn0=jspn; ispn1=jspn
  else
    ispn0=1; ispn1=nspinor
  end if
  i=1
  do ihp=1,nhp(jspn)
    do igp=i,ngp(jspn)
      t1=abs(vhpc(1,ihp,jspn)-vgpc(1,igp,jspn)) &
        +abs(vhpc(2,ihp,jspn)-vgpc(2,igp,jspn)) &
        +abs(vhpc(3,ihp,jspn)-vgpc(3,igp,jspn))
      if (t1.lt.epslat) then
        do ist=1,nstsv
          do ispn=ispn0,ispn1
            wfpw(ihp,ispn,ist)=wfir(igp,ispn,ist)
          end do
        end do
        if (igp.eq.i) i=i+1
        exit
      end if
    end do
  end do
end do
!-------------------------!
!     muffin-tin part     !
!-------------------------!
allocate(jl(0:lmaxo,nrcmtmax))
t0=fourpi/sqrt(omega)
! remove continuation of interstitial function into muffin-tin
do jspn=1,nspnfv
  if (spinsprl) then
    ispn0=jspn; ispn1=jspn
  else
    ispn0=1; ispn1=nspinor
  end if
! loop over G+p-vectors
  do igp=1,ngp(jspn)
! generate the conjugate spherical harmonics Y_lm*(G+p)
    call genylm(lmaxo,tpgpc(:,igp,jspn),ylm)
    ylm(:)=conjg(ylm(:))
    do is=1,nspecies
      nrc=nrcmt(is)
      nrci=nrcmti(is)
      irco=nrci+1
      npci=npcmti(is)
! generate spherical Bessel functions
      lmax=lmaxi
      do irc=1,nrc
        t1=gpc(igp,jspn)*rcmt(irc,is)
        call sbessel(lmax,t1,jl(:,irc))
        if (irc.eq.nrci) lmax=lmaxo
      end do
! loop over atoms
      do ia=1,natoms(is)
        ias=idxas(ia,is)
        z1=t0*sfacgp(igp,ias,jspn)
        do ist=1,nstsv
          do ispn=ispn0,ispn1
            z2=z1*wfir(igp,ispn,ist)
            lm=0
            do l=0,lmaxi
              z3=z2*zil(l)
              do m=-l,l
                lm=lm+1
                z4=z3*ylm(lm)
                i=lm
                do irc=1,nrci
                  wfmt(i,ias,ispn,ist)=wfmt(i,ias,ispn,ist)-z4*jl(l,irc)
                  i=i+lmmaxi
                end do
              end do
            end do
            lm=0
            do l=0,lmaxo
              z3=z2*zil(l)
              do m=-l,l
                lm=lm+1
                z4=z3*ylm(lm)
                i=npci+lm
                do irc=irco,nrc
                  wfmt(i,ias,ispn,ist)=wfmt(i,ias,ispn,ist)-z4*jl(l,irc)
                  i=i+lmmaxo
                end do
              end do
            end do
          end do
        end do
      end do
    end do
  end do
end do
! Fourier transform the muffin-tin wavefunctions
do jspn=1,nspnfv
  if (spinsprl) then
    ispn0=jspn; ispn1=jspn
  else
    ispn0=1; ispn1=nspinor
  end if
! loop over H+p-vectors
  do ihp=1,nhp(jspn)
! generate the spherical harmonics Y_lm(H+p)
    call genylm(lmaxo,tphpc(:,ihp,jspn),ylm)
    do is=1,nspecies
      nrc=nrcmt(is)
      nrci=nrcmti(is)
! generate spherical Bessel functions
      lmax=lmaxi
      do irc=1,nrc
        t1=hpc(ihp,jspn)*rcmt(irc,is)
        call sbessel(lmax,t1,jl(:,irc))
        if (irc.eq.nrci) lmax=lmaxo
      end do
      do ia=1,natoms(is)
        ias=idxas(ia,is)
! conjugate structure factor
        z1=t0*conjg(sfachp(ihp,ias,jspn))
! loop over states
        do ist=1,nstsv
          do ispn=ispn0,ispn1
            lmax=lmaxi
            i=0
            do irc=1,nrc
              i=i+1
              zsum1=jl(0,irc)*wfmt(i,ias,ispn,ist)*ylm(1)
              lm=1
              do l=1,lmax
                lm=lm+1
                i=i+1
                zsum2=wfmt(i,ias,ispn,ist)*ylm(lm)
                do m=1-l,l
                  lm=lm+1
                  i=i+1
                  zsum2=zsum2+wfmt(i,ias,ispn,ist)*ylm(lm)
                end do
                zsum1=zsum1+jl(l,irc)*zilc(l)*zsum2
              end do
              zsum1=zsum1*r2cmt(irc,is)
              fr1(irc)=dble(zsum1)
              fr2(irc)=aimag(zsum1)
              if (irc.eq.nrci) lmax=lmaxo
            end do
            t1=splint(nrc,rcmt(:,is),fr1)
            t2=splint(nrc,rcmt(:,is),fr2)
! add to the H+p wavefunction
            wfpw(ihp,ispn,ist)=wfpw(ihp,ispn,ist)+z1*cmplx(t1,t2,8)
          end do
        end do
      end do
    end do
  end do
end do
deallocate(jl,wfmt,wfir)
return
end subroutine