Line data Source code
1 : !!****f* ABINIT/nfourier
2 : !! NAME
3 : !! nfourier
4 : !!
5 : !! FUNCTION
6 : !!
7 : !! COPYRIGHT
8 : !! Copyright (C) 2014-2026 ABINIT group (XG)
9 : !! This file is distributed under the terms of the
10 : !! GNU General Public License, see ~abinit/COPYING
11 : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : !!
13 : !! SOURCE
14 :
15 : #if defined HAVE_CONFIG_H
16 : #include "config.h"
17 : #endif
18 :
19 : #include "abi_common.h"
20 :
21 : MODULE m_lib_four
22 :
23 : use defs_basis
24 : use m_errors
25 : use m_abicore
26 :
27 : contains
28 :
29 : ! This routine contains direct and inverse fourier transformation
30 : ! It is a modification of a routine of the GNU GPL
31 : ! code available on http://dmft.rutgers.edu/ and
32 : ! described in the RMP 2006 paper written by
33 : ! G.Kotliar, S.Y.Savrasov, K.Haule, V.S.Oudovenko, O.Parcollet, C.A.Marianetti
34 : !=======+=========+=========+=========+=========+=========+=========+=$
35 : ! TYPE : SUBROUTINE
36 : ! PROGRAM: nfourier
37 : ! PURPOSE: fourier-transform the natural-spline interpolation
38 : ! of function Green(tau)
39 : ! calculate function Green(omega)
40 : ! I/O :
41 : ! VERSION: 2-16-92
42 : ! 29-Nov-95 removal of minimal bug concerning
43 : ! DIMENSION of rindata
44 : ! COMMENT: cf J. Stoer R. Bulirsch, Introduction to numerical
45 : ! analysis (Springer, New York, 1980)
46 : !=======+=========+=========+=========+=========+=========+=========+=$
47 : !
48 0 : SUBROUTINE nfourier(rindata,coutdata,iflag,Iwmax,L,Beta)
49 :
50 : !include 'param.dat'
51 : IMPLICIT DOUBLE PRECISION (A-H,O-Z)
52 : IMPLICIT INTEGER(I-N)
53 : integer Iwmax,L,iflag
54 : DIMENSION rindata(L)
55 0 : DIMENSION rincopy(L+1),a(L),b(L),c(L),d(L),u(L+1), q(L+1),XM(L+1)
56 : complex*16 :: coutdata(Iwmax+1)
57 : complex*16 cdummy,explus,ex
58 0 : xpi = ACOS(-One)
59 0 : delta = Beta/L
60 0 : DO i = 1,L
61 0 : rincopy(i) = rindata(i)
62 : ENDDO
63 0 : if(iflag==1) then
64 0 : rincopy(L+1) = -1-rindata(1)
65 : else
66 0 : rincopy(L+1) = -rindata(1)
67 : endif
68 : ! Three = Two+One
69 : ! six = Two*Three
70 :
71 : !c
72 : !c spline interpolation: the spline is given by
73 : !c G(tau) = a(i) + b(i) (tau-tau_i) + c(i) ( )^2 + d(i) ( )^3
74 : !c The following formulas are taken directly from Stoer and
75 : !c Bulirsch p. 102
76 : !c
77 0 : q(1) = Zero
78 0 : u(1) = Zero
79 0 : DO k = 2,L
80 0 : p = q(k-1)/Two+Two
81 0 : q(k)=-One/Two/p
82 0 : u(k)=Three/delta**2*(rincopy(k+1)+rincopy(k-1)-Two*rincopy(k))
83 0 : u(k)=(u(k)-u(k-1)/Two)/p
84 : ENDDO
85 0 : XM(L+1) = 0
86 0 : DO k = L,1,-1
87 0 : XM(k) = q(k)*XM(k+1)+u(k)
88 : ENDDO
89 : !c
90 : !c The following formulas are taken directly from Stoer and
91 : !c Bulirsch p. 98
92 : !c
93 0 : DO j = 1, L
94 0 : a(j) = rincopy(j)
95 0 : c(j) = XM(j)/Two
96 : b(j) = (rincopy(j+1)-rincopy(j))/delta - &
97 0 : & (Two*XM(j)+XM(j+1))*delta/6.
98 0 : d(j) = (XM(j+1)-XM(j))/(6.*delta)
99 : ENDDO
100 :
101 : !c
102 : !c The Spline multiplied by the exponential can now be exlicitely
103 : !c integrated. The following formulas were obtained using
104 : !c MATHEMATICA
105 : !c
106 0 : DO i = 0,Iwmax
107 0 : om = (Two*(i)+One)*xpi/Beta
108 0 : coutdata(i+1) = czero
109 0 : DO j = 1,L
110 0 : cdummy = j_dpc*om*delta*j
111 0 : explus = exp(cdummy)
112 0 : cdummy = j_dpc*om*delta*(j-1)
113 0 : ex = exp(cdummy)
114 : coutdata(i+1) = coutdata(i+1) + explus*(&
115 : & ( -six* d(j) )/om**4 + &
116 : & ( Two*j_dpc*c(j) + six*delta*j_dpc*d(j) )/om**3 +&
117 : & ( b(j)+ Two*delta*c(j)+ three*delta**2*d(j) )/om**2 +&
118 : & (- j_dpc*a(j) - delta*j_dpc*b(j) - delta**2*j_dpc*c(j) -&
119 0 : & delta**3*j_dpc*d(j))/om)
120 :
121 : coutdata(i+1) = coutdata(i+1) + ex*(&
122 : & six*d(j)/om**4 - Two*j_dpc*c(j)/om**3 &
123 0 : & -b(j)/om**2 + j_dpc*a(j)/om)
124 : ENDDO
125 : ENDDO
126 0 : end subroutine nfourier
127 : !!***
128 :
129 : ! This routine contains direct and inverse fourier transformation
130 : ! It is a modification of a routine of the GNU GPL
131 : ! code available on http://dmft.rutgers.edu/ and
132 : ! described in the RMP 2006 paper written by
133 : ! G.Kotliar, S.Y.Savrasov, K.Haule, V.S.Oudovenko, O.Parcollet, C.A.Marianetti
134 : !=======+=========+=========+=========+=========+=========+=========+=$
135 : ! TYPE : SUBROUTINE
136 : ! PROGRAM: nfourier
137 : ! PURPOSE: fourier-transform the natural-spline interpolation
138 : ! of function Green(tau)
139 : ! calculate function Green(omega)
140 : ! I/O :
141 : ! VERSION: 2-16-92
142 : ! 29-Nov-95 removal of minimal bug concerning
143 : ! DIMENSION of rindata
144 : ! COMMENT: cf J. Stoer R. Bulirsch, Introduction to numerical
145 : ! analysis (Springer, New York, 1980)
146 : !=======+=========+=========+=========+=========+=========+=========+=$
147 : !
148 1600 : SUBROUTINE nfourier2(rindata,coutdata,iflag,om,L,Beta)
149 : ! include 'param.dat'
150 : IMPLICIT DOUBLE PRECISION (A-H,O-Z)
151 : IMPLICIT INTEGER(I-N)
152 : integer L,iflag
153 : DIMENSION rindata(L)
154 3200 : DIMENSION rincopy(L+1),a(L),b(L),c(L),d(L),u(L+1), q(L+1),XM(L+1)
155 : complex*16 :: coutdata
156 : real*8 :: om
157 : complex*16 cdummy,explus,ex
158 1600 : xpi = ACOS(-One)
159 1600 : delta = Beta/L
160 104000 : DO i = 1,L
161 104000 : rincopy(i) = rindata(i)
162 : ENDDO
163 1600 : if(iflag==1 .and. L.ge.1) then
164 320 : rincopy(L+1) = -1-rindata(1)
165 1280 : elseif(iflag==0 .and. L.ge.1) then
166 1280 : rincopy(L+1) = -rindata(1)
167 : else
168 0 : write(std_out,*) "Warning : Check nfourier2"
169 : endif
170 : ! Three = Two+One
171 : ! six = Two*Three
172 :
173 : !c
174 : !c spline interpolation: the spline is given by
175 : !c G(tau) = a(i) + b(i) (tau-tau_i) + c(i) ( )^2 + d(i) ( )^3
176 : !c The following formulas are taken directly from Stoer and
177 : !c Bulirsch p. 102
178 : !c
179 1600 : q(1) = Zero
180 1600 : u(1) = Zero
181 102400 : DO k = 2,L
182 100800 : p = q(k-1)/Two+Two
183 100800 : q(k)=-One/Two/p
184 100800 : u(k)=Three/delta**2*(rincopy(k+1)+rincopy(k-1)-Two*rincopy(k))
185 102400 : u(k)=(u(k)-u(k-1)/Two)/p
186 : ENDDO
187 1600 : XM(L+1) = 0
188 104000 : DO k = L,1,-1
189 104000 : XM(k) = q(k)*XM(k+1)+u(k)
190 : ENDDO
191 : !c
192 : !c The following formulas are taken directly from Stoer and
193 : !c Bulirsch p. 98
194 : !c
195 104000 : DO j = 1, L
196 102400 : a(j) = rincopy(j)
197 102400 : c(j) = XM(j)/Two
198 : b(j) = (rincopy(j+1)-rincopy(j))/delta - &
199 102400 : & (Two*XM(j)+XM(j+1))*delta/6.
200 104000 : d(j) = (XM(j+1)-XM(j))/(6.*delta)
201 : ENDDO
202 :
203 : !c
204 : !c The Spline multiplied by the exponential can now be exlicitely
205 : !c integrated. The following formulas were obtained using
206 : !c MATHEMATICA
207 : !c
208 1600 : coutdata = czero
209 104000 : DO j = 1,L
210 102400 : cdummy = j_dpc*om*delta*j
211 102400 : explus = exp(cdummy)
212 102400 : cdummy = j_dpc*om*delta*(j-1)
213 102400 : ex = exp(cdummy)
214 : coutdata = coutdata + explus*(&
215 : & ( -six* d(j) )/om**4 + &
216 : & ( Two*j_dpc*c(j) + six*delta*j_dpc*d(j) )/om**3 +&
217 : & ( b(j)+ Two*delta*c(j)+ three*delta**2*d(j) )/om**2 +&
218 : & (- j_dpc*a(j) - delta*j_dpc*b(j) - delta**2*j_dpc*c(j) -&
219 102400 : & delta**3*j_dpc*d(j))/om)
220 :
221 : coutdata = coutdata + ex*(&
222 : & six*d(j)/om**4 - Two*j_dpc*c(j)/om**3 &
223 104000 : & -b(j)/om**2 + j_dpc*a(j)/om)
224 : ENDDO
225 1600 : end subroutine nfourier2
226 : !=======+=========+=========+=========+=========+=========+=========+=$
227 : ! TYPE : SUBROUTINE
228 : ! PROGRAM: invfourier
229 : ! PURPOSE: inverse fourier transform
230 : ! Greent, Greenw use physical definition
231 : ! Greent(i) = G((i-1)*deltau) for i = 1,...,L
232 : ! Greenw(n) = G(i w_n), for n = 0,L/2-1
233 : ! w_n = (2*n+1)pi/beta
234 : ! Symmetry property:
235 : ! G(iw_(-n) = G(iw_(n-1))*
236 : ! coupled to the impurity
237 : ! I/O :
238 : ! VERSION: 6-16-92
239 : ! COMMENT:
240 : !=======+=========+=========+=========+=========+=========+=========+=$
241 : !
242 50 : SUBROUTINE invfourier(cindata,routdata,Iwmax,L,iflag,beta)
243 :
244 : ! include 'param.dat'
245 : implicit none
246 : integer, intent(in) :: Iwmax
247 : complex*16, intent(in) :: cindata(1:Iwmax) !vz_d
248 : integer, intent(in) :: L
249 : complex*16, intent(inout) :: routdata(1:L) !vz_d
250 : integer, intent(in) :: iflag
251 : double precision, intent(in) :: beta
252 :
253 : double precision :: xpi
254 : double precision :: tau
255 : double precision :: om
256 : complex*16 :: cdummy,dummy
257 : integer :: i,j
258 :
259 50 : xpi = ACOS(-One)
260 3250 : DO 1 i = 1,L
261 3200 : routdata(i) = Zero
262 3200 : tau = (i-1)*beta/real(L)
263 64003200 : DO 2 j = 1,Iwmax
264 : ! om = mod((2*(j)+One)*xpi/Beta*tau,2*xpi)
265 64000000 : om = ((2*(j)-One)*xpi/Beta*tau)
266 64000000 : cdummy = CMPLX(Zero,om)
267 64000000 : dummy = cindata(j)*exp(-cdummy)
268 64000000 : routdata(i) = routdata(i)+Two/beta*dummy
269 3200 : 2 CONTINUE
270 : ! write(std_out,*) "FT",i,routdata(i)
271 50 : 1 CONTINUE
272 : !c
273 : !c special treatment for tau = 0
274 : !c
275 50 : if(iflag==1 .and. L.ge.1 ) then
276 10 : routdata(1) = -One/Two+routdata(1)
277 : endif
278 50 : END SUBROUTINE invfourier
279 :
280 : END MODULE m_lib_four
281 : !!***
|
