Branch data 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 : : !!***
|
