1
2! Copyright (C) 2018 J. K. Dewhurst, S. Sharma and E. K. U. Gross.
3! This file is distributed under the terms of the GNU General Public License.
4! See the file COPYING for license details.
5
6subroutine zfplot(np,vpl,zfmt,zfir,fp)
7use modmain
8use modomp
9implicit none
10! arguments
11integer, intent(in) :: np
12real(8), intent(in) :: vpl(3,np)
13complex(8), intent(in) :: zfmt(npcmtmax,natmtot),zfir(ngtot)
14complex(8), intent(out) :: fp(np)
15! local variables
16integer ias,is,ip,nthd
17! allocatable arrays
18complex(8), allocatable :: zfmt1(:,:,:),zfft(:)
19! unpack the muffin-tin function
20allocate(zfmt1(lmmaxo,nrcmtmax,natmtot))
21call holdthd(natmtot,nthd)
22!$OMP PARALLEL DO DEFAULT(SHARED) &
23!$OMP PRIVATE(is) &
24!$OMP NUM_THREADS(nthd)
25do ias=1,natmtot
26  is=idxis(ias)
27  call zfmtpack(.false.,nrcmt(is),nrcmti(is),zfmt(:,ias),zfmt1(:,:,ias))
28end do
29!$OMP END PARALLEL DO
30call freethd(nthd)
31! Fourier transform rfir to G-space
32allocate(zfft(ngtot))
33zfft(:)=zfir(:)
34call zfftifc(3,ngridg,-1,zfft)
35! begin loop over all points
36call holdthd(np,nthd)
37!$OMP PARALLEL DO DEFAULT(SHARED) &
38!$OMP NUM_THREADS(nthd)
39do ip=1,np
40  call zfip(ip)
41end do
42!$OMP END PARALLEL DO
43call freethd(nthd)
44deallocate(zfmt1,zfft)
45return
46
47contains
48
49subroutine zfip(ip)
50implicit none
51! arguments
52integer, intent(in) :: ip
53! local variables
54integer is,ia,ias,nrc,nrci
55integer irc0,irc,lmax,l,m,lm
56integer ig,ifg,i1,i2,i3,i,j
57real(8) rmt2,r,ya1(4),ya2(4),t1,t2
58real(8) v1(3),v2(3),v3(3),v4(3),v5(3)
59complex(8) z1
60! automatic arrays
61complex(8) ylm(lmmaxo)
62v2(:)=vpl(:,ip)
63call r3frac(epslat,v2)
64! convert point to Cartesian coordinates
65call r3mv(avec,v2,v1)
66! check if point is in a muffin-tin
67do is=1,nspecies
68  nrc=nrcmt(is)
69  nrci=nrcmti(is)
70  rmt2=rmt(is)**2
71  do ia=1,natoms(is)
72    ias=idxas(ia,is)
73    v2(:)=v1(:)-atposc(:,ia,is)
74    do i1=-1,1
75      v3(:)=v2(:)+dble(i1)*avec(:,1)
76      do i2=-1,1
77        v4(:)=v3(:)+dble(i2)*avec(:,2)
78        do i3=-1,1
79          v5(:)=v4(:)+dble(i3)*avec(:,3)
80          t1=v5(1)**2+v5(2)**2+v5(3)**2
81          if (t1.lt.rmt2) then
82            r=sqrt(t1)
83            call genylmv(lmaxo,v5,ylm)
84            do irc=1,nrc
85              if (rcmt(irc,is).ge.r) then
86                if (irc.le.3) then
87                  irc0=1
88                else if (irc.gt.nrc-2) then
89                  irc0=nrc-3
90                else
91                  irc0=irc-2
92                end if
93                r=max(r,rcmt(1,is))
94                if (irc0.le.nrci) then
95                  lmax=lmaxi
96                else
97                  lmax=lmaxo
98                end if
99                z1=0.d0
100                lm=0
101                do l=0,lmax
102                  do m=-l,l
103                    lm=lm+1
104                    do j=1,4
105                      i=irc0+j-1
106                      ya1(j)=dble(zfmt1(lm,i,ias))
107                      ya2(j)=aimag(zfmt1(lm,i,ias))
108                    end do
109                    t1=poly4(rcmt(irc0,is),ya1,r)
110                    t2=poly4(rcmt(irc0,is),ya2,r)
111                    z1=z1+cmplx(t1,t2,8)*ylm(lm)
112                  end do
113                end do
114                goto 10
115              end if
116            end do
117          end if
118        end do
119      end do
120    end do
121  end do
122end do
123! otherwise use direct Fourier transform of interstitial function
124z1=0.d0
125do ig=1,ngvec
126  ifg=igfft(ig)
127  t1=vgc(1,ig)*v1(1)+vgc(2,ig)*v1(2)+vgc(3,ig)*v1(3)
128  z1=z1+zfft(ifg)*cmplx(cos(t1),sin(t1),8)
129end do
13010 continue
131fp(ip)=z1
132end subroutine
133
134pure real(8) function poly4(xa,ya,x)
135implicit none
136! arguments
137real(8), intent(in) :: xa(4),ya(4),x
138! local variables
139real(8) x0,x1,x2,x3,y0,y1,y2,y3
140real(8) c1,c2,c3,t0,t1,t2,t3,t4,t5,t6
141! evaluate the polynomial coefficients
142x0=xa(1)
143x1=xa(2)-x0; x2=xa(3)-x0; x3=xa(4)-x0
144t4=x1-x2; t5=x1-x3; t6=x2-x3
145y0=ya(1)
146y1=ya(2)-y0; y2=ya(3)-y0; y3=ya(4)-y0
147t1=x1*x2*y3; t2=x2*x3*y1; t3=x1*x3
148t0=1.d0/(x2*t3*t4*t5*t6)
149t3=t3*y2
150c3=t1*t4+t2*t6-t3*t5
151t4=x1**2; t5=x2**2; t6=x3**2
152c2=t1*(t5-t4)+t2*(t6-t5)+t3*(t4-t6)
153c1=t1*(x2*t4-x1*t5)+t2*(x3*t5-x2*t6)+t3*(x1*t6-x3*t4)
154t1=x-x0
155! evaluate the polynomial
156poly4=y0+t0*t1*(c1+t1*(c2+c3*t1))
157end function
158
159end subroutine
160
161