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
|
! Copyright (C) 2018 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 zfpts(np,vrl,zfmt,zfir,fp)
use modmain
use modomp
implicit none
! arguments
integer, intent(in) :: np
real(8), intent(in) :: vrl(3,np)
complex(8), intent(in) :: zfmt(npcmtmax,natmtot),zfir(ngtot)
complex(8), intent(out) :: fp(np)
! local variables
integer ip,nthd
! allocatable arrays
complex(8), allocatable :: zfft(:)
! Fourier transform zfir to G-space
allocate(zfft(ngtot))
zfft(1:ngtot)=zfir(1:ngtot)
call zfftifc(3,ngridg,-1,zfft)
! begin loop over all points
call holdthd(np,nthd)
!$OMP PARALLEL DO DEFAULT(SHARED) &
!$OMP SCHEDULE(DYNAMIC) &
!$OMP NUM_THREADS(nthd)
do ip=1,np
call zfip(ip)
end do
!$OMP END PARALLEL DO
call freethd(nthd)
deallocate(zfft)
contains
subroutine zfip(ip)
implicit none
! arguments
integer, intent(in) :: ip
! local variables
integer is,ias,nrc,nrci
integer ir,irc,irc0,i0,j
integer lmax,lmmax,lm,ig
real(8) v(3),r,t1
complex(8) zsm,z1
! automatic arrays
complex(8) ya(4,lmmaxo),ylm(lmmaxo)
! check if point is in a muffin-tin
call findmtpt(vrl(:,ip),ias,ir,v,r)
if (ir > 0) then
! point is inside a muffin-tin
is=idxis(ias)
nrc=nrcmt(is)
nrci=nrcmti(is)
irc=(ir-1)/lradstp+1
if (irc <= 3) then
irc0=1
else if (irc > nrc-2) then
irc0=nrc-3
else
irc0=irc-2
end if
r=max(r,rcmt(1,is))
if (irc0 <= nrci) then
lmax=lmaxi
lmmax=lmmaxi
else
lmax=lmaxo
lmmax=lmmaxo
end if
do j=1,4
irc=irc0+j-1
if (irc <= nrci) then
i0=lmmaxi*(irc-1)
else
i0=lmmaxi*nrci+lmmaxo*(irc-nrci-1)
end if
ya(j,1:lmmax)=zfmt(i0+1:i0+lmmax,ias)
end do
! generate complex spherical harmonics
call genylmv(.false.,lmaxo,v,ylm)
zsm=0.d0
do lm=1,lmmax
z1=poly4(rcmt(irc0,is),ya(:,lm),r)
zsm=zsm+z1*ylm(lm)
end do
else
! otherwise use direct Fourier transform of interstitial function
v(1:3)=vrl(1,ip)*avec(1:3,1)+vrl(2,ip)*avec(1:3,2)+vrl(3,ip)*avec(1:3,3)
zsm=0.d0
!$OMP SIMD PRIVATE(t1) REDUCTION(+:zsm)
do ig=1,ngvec
t1=vgc(1,ig)*v(1)+vgc(2,ig)*v(2)+vgc(3,ig)*v(3)
zsm=zsm+zfft(igfft(ig))*cmplx(cos(t1),sin(t1),8)
end do
end if
fp(ip)=zsm
end subroutine
pure complex(8) function poly4(xa,ya,x)
implicit none
! arguments
real(8), intent(in) :: xa(4),x
complex(8), intent(in) :: ya(4)
! local variables
real(8) x0,x1,x2,x3,t0,t1,t2,t3,t4
complex(8) y0,y1,y2,y3,z1,z2,z3,c1,c2,c3
! evaluate the polynomial coefficients
x0=xa(1)
x1=xa(2)-x0; x2=xa(3)-x0; x3=xa(4)-x0
t1=x1-x2; t2=x1-x3; t3=x2-x3
y0=ya(1)
y1=ya(2)-y0; y2=ya(3)-y0; y3=ya(4)-y0
z1=x1*x2*y3; z2=x2*x3*y1; t4=x1*x3
t0=1.d0/(x2*t1*t2*t3*t4)
z3=t4*y2
c3=z1*t1+z2*t3-z3*t2
t1=x1**2; t2=x2**2; t3=x3**2
c2=z1*(t2-t1)+z2*(t3-t2)+z3*(t1-t3)
c1=z1*(x2*t1-x1*t2)+z2*(x3*t2-x2*t3)+z3*(x1*t3-x3*t1)
t1=x-x0
! evaluate the polynomial
poly4=y0+t0*t1*(c1+t1*(c2+c3*t1))
end function
end subroutine
|
