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
|
! Copyright (C) 2002-2005 J. K. Dewhurst, S. Sharma and C. Ambrosch-Draxl.
! This file is distributed under the terms of the GNU Lesser General Public
! License. See the file COPYING for license details.
!BOP
! !ROUTINE: brzint
! !INTERFACE:
subroutine brzint(nsm,ngridk,nsk,ivkik,nw,wint,n,ld,e,f,g)
! !USES:
use modomp
! !INPUT/OUTPUT PARAMETERS:
! nsm : level of smoothing for output function (in,integer)
! ngridk : k-point grid size (in,integer(3))
! nsk : k-point subdivision grid size (in,integer(3))
! ivkik : map from (i1,i2,i3) to k-point index
! (in,integer(0:ngridk(1)-1,0:ngridk(2)-1,0:ngridk(3)-1))
! nw : number of energy divisions (in,integer)
! wint : energy interval (in,real(2))
! n : number of functions to integrate (in,integer)
! ld : leading dimension (in,integer)
! e : array of energies as a function of k-points (in,real(ld,*))
! f : array of weights as a function of k-points (in,real(ld,*))
! g : output function (out,real(nw))
! !DESCRIPTION:
! Given energy and weight functions, $e$ and $f$, on the Brillouin zone and a
! set of equidistant energies $\omega_i$, this routine computes the integrals
! $$ g(\omega_i)=\frac{\Omega}{(2\pi)^3}\int_{\rm BZ} f({\bf k})
! \delta(\omega_i-e({\bf k}))d{\bf k}, $$
! where $\Omega$ is the unit cell volume. This is done by first interpolating
! $e$ and $f$ on a finer $k$-point grid using the trilinear method. Then for
! each $e({\bf k})$ on the finer grid the nearest $\omega_i$ is found and
! $f({\bf k})$ is accumulated in $g(\omega_i)$. If the output function is
! noisy then either {\tt nsk} should be increased or {\tt nw} decreased.
! Alternatively, the output function can be artificially smoothed up to a
! level given by {\tt nsm}. See routine {\tt fsmooth}.
!
! !REVISION HISTORY:
! Created October 2003 (JKD)
! Improved efficiency, May 2007 (Sebastian Lebegue)
! Added parallelism, March 2020 (JKD)
!EOP
!BOC
implicit none
! arguments
integer, intent(in) :: nsm,ngridk(3),nsk(3)
integer, intent(in) :: ivkik(0:ngridk(1)-1,0:ngridk(2)-1,0:ngridk(3)-1)
integer, intent(in) :: nw
real(8), intent(in) :: wint(2)
integer, intent(in) :: n,ld
real(8), intent(in) :: e(ld,*),f(ld,*)
real(8), intent(out) :: g(nw)
! local variables
integer nk,i1,i2,i3,j1,j2,j3,k1,k2,k3,i,iw,nthd
integer i000,i001,i010,i011,i100,i101,i110,i111
real(8) wd,dw,dwi,w1,t1,t2
! automatic arrays
real(8) f0(n),f1(n),e0(n),e1(n)
real(8) f00(n),f01(n),f10(n),f11(n)
real(8) e00(n),e01(n),e10(n),e11(n)
if ((ngridk(1) < 1).or.(ngridk(2) < 1).or.(ngridk(3) < 1)) then
write(*,*)
write(*,'("Error(brzint): ngridk < 1 :",3(X,I0))') ngridk
write(*,*)
stop
end if
if ((nsk(1) < 1).or.(nsk(2) < 1).or.(nsk(3) < 1)) then
write(*,*)
write(*,'("Error(brzint): nsk < 1 :",3(X,I0))') nsk
write(*,*)
stop
end if
! total number of k-points
nk=ngridk(1)*ngridk(2)*ngridk(3)
! length of interval
wd=wint(2)-wint(1)
! energy step size
dw=wd/dble(nw)
dwi=1.d0/dw
w1=wint(1)*dwi
g(1:nw)=0.d0
call holdthd(nk,nthd)
!$OMP PARALLEL DO DEFAULT(SHARED) &
!$OMP PRIVATE(f0,f1,e0,e1) &
!$OMP PRIVATE(f00,f01,f10,f11) &
!$OMP PRIVATE(e00,e01,e10,e11) &
!$OMP PRIVATE(k1,k2,k3,i1,i2,i3) &
!$OMP PRIVATE(i000,i001,i010,i011) &
!$OMP PRIVATE(i100,i101,i110,i111) &
!$OMP PRIVATE(t1,t2,i,iw) &
!$OMP REDUCTION(+:g) &
!$OMP SCHEDULE(DYNAMIC) &
!$OMP NUM_THREADS(nthd) &
!$OMP COLLAPSE(3)
do j1=0,ngridk(1)-1
do j2=0,ngridk(2)-1
do j3=0,ngridk(3)-1
k1=mod(j1+1,ngridk(1))
k2=mod(j2+1,ngridk(2))
k3=mod(j3+1,ngridk(3))
i000=ivkik(j1,j2,j3); i001=ivkik(j1,j2,k3)
i010=ivkik(j1,k2,j3); i011=ivkik(j1,k2,k3)
i100=ivkik(k1,j2,j3); i101=ivkik(k1,j2,k3)
i110=ivkik(k1,k2,j3); i111=ivkik(k1,k2,k3)
do i1=0,nsk(1)-1
t2=dble(i1)/dble(nsk(1))
t1=1.d0-t2
f00(1:n)=f(1:n,i000)*t1+f(1:n,i100)*t2
f01(1:n)=f(1:n,i001)*t1+f(1:n,i101)*t2
f10(1:n)=f(1:n,i010)*t1+f(1:n,i110)*t2
f11(1:n)=f(1:n,i011)*t1+f(1:n,i111)*t2
t1=t1*dwi
t2=t2*dwi
e00(1:n)=e(1:n,i000)*t1+e(1:n,i100)*t2-w1
e01(1:n)=e(1:n,i001)*t1+e(1:n,i101)*t2-w1
e10(1:n)=e(1:n,i010)*t1+e(1:n,i110)*t2-w1
e11(1:n)=e(1:n,i011)*t1+e(1:n,i111)*t2-w1
do i2=0,nsk(2)-1
t2=dble(i2)/dble(nsk(2))
t1=1.d0-t2
f0(1:n)=f00(1:n)*t1+f10(1:n)*t2
f1(1:n)=f01(1:n)*t1+f11(1:n)*t2
e0(1:n)=e00(1:n)*t1+e10(1:n)*t2
e1(1:n)=e01(1:n)*t1+e11(1:n)*t2
do i3=0,nsk(3)-1
t2=dble(i3)/dble(nsk(3))
t1=1.d0-t2
do i=1,n
iw=nint(e0(i)*t1+e1(i)*t2)+1
if ((iw >= 1).and.(iw <= nw)) g(iw)=g(iw)+f0(i)*t1+f1(i)*t2
end do
end do
end do
end do
end do
end do
end do
!$OMP END PARALLEL DO
call freethd(nthd)
! normalise function
t1=dw*dble(nk)*dble(nsk(1)*nsk(2)*nsk(3))
t1=1.d0/t1
g(1:nw)=t1*g(1:nw)
! smooth output function if required
if (nsm > 0) call fsmooth(nsm,nw,g)
end subroutine
!EOC
|
