Line data Source code
1 : MODULE m_kkintgr
2 :
3 : !------------------------------------------------------------------------------
4 : !
5 : ! MODULE: m_kkintgr
6 : !
7 : !> @author
8 : !> Henning Janßen
9 : !
10 : ! DESCRIPTION:
11 : !> Performs the Kramer-Kronig-Transformation to obtain the Green's function
12 : !> in the complex plane from the imaginary part calculated on the real axis
13 : !
14 : ! TODO: Look at FFT for Transformation
15 : ! How to do changing imaginary parts
16 : !------------------------------------------------------------------------------
17 :
18 : USE m_constants
19 : USE m_juDFT
20 : USE m_types_mat
21 : USE m_smooth
22 : USE m_lorentzian_smooth
23 :
24 : IMPLICIT NONE
25 :
26 : PRIVATE
27 :
28 : INTEGER, PARAMETER :: method_maclaurin = 1
29 : INTEGER, PARAMETER :: method_deriv = 2
30 : INTEGER, PARAMETER :: method_direct = 3
31 : INTEGER, PARAMETER :: method_fft = 4
32 :
33 : CHARACTER(len=10), PARAMETER :: smooth_method = 'lorentzian' !(lorentzian or gaussian)
34 : INTEGER, PARAMETER :: int_method(3) = [method_direct,method_direct,method_maclaurin]
35 :
36 : TYPE(t_mat),SAVE, ALLOCATABLE :: integration_weights(:)
37 : INTEGER,SAVE, ALLOCATABLE :: methods(:)
38 : REAL,SAVE, ALLOCATABLE :: sigma(:)
39 : REAL,SAVE, ALLOCATABLE :: energy_grid(:)
40 : REAL, SAVE :: del = -1.0
41 :
42 : INTERFACE kkintgr
43 : PROCEDURE :: kkintgr_single, kkintgr_mmpmat
44 : PROCEDURE :: kkintgr_mmpmat_extra1
45 : END INTERFACE
46 :
47 : PUBLIC kkintgr, kkintgr_init, kkintgr_free
48 :
49 : CONTAINS
50 :
51 517 : SUBROUTINE kkintgr_init(eMesh, contour, iContour, nContour, shape, additional_smearing)
52 :
53 : REAL, INTENT(IN) :: eMesh(:) !Energy grid on the real axis
54 : COMPLEX, INTENT(IN) :: contour(:) !Complex energy contour
55 : INTEGER, INTENT(IN) :: iContour,nContour,shape
56 : real, intent(in) :: additional_smearing
57 :
58 : INTEGER :: iz,izp,i,n1,n2
59 : REAL :: y
60 :
61 517 : IF(.NOT.ALLOCATED(integration_weights)) THEN
62 172 : ALLOCATE(integration_weights(nContour))
63 172 : ALLOCATE(methods(nContour), source=0)
64 172 : ALLOCATE(sigma(nContour), source=0.0)
65 271284 : energy_grid = eMesh
66 : ENDIF
67 :
68 517 : IF(integration_weights(iContour)%allocated()) RETURN
69 :
70 43 : IF(del < 0.0) del = eMesh(2)-eMesh(1)
71 43 : IF(ABS(del-(eMesh(2)-eMesh(1)))>1e-12) CALL juDFT_error("Inconsistent energy grids", calledby="kkintgr_init")
72 43 : methods(iContour) = int_method(shape)
73 :
74 43 : CALL integration_weights(iContour)%init(.false.,SIZE(eMesh),SIZE(contour))
75 :
76 43 : IF(int_method(shape) == method_direct) THEN
77 5463 : DO iz = 1, SIZE(contour)
78 34962021 : integration_weights(iContour)%data_c(:,iz) = 1.0/(contour(iz)-eMesh)
79 5421 : integration_weights(iContour)%data_c(1,iz) = integration_weights(iContour)%data_c(1,iz)/2.0
80 5463 : integration_weights(iContour)%data_c(SIZE(eMesh),iz) = integration_weights(iContour)%data_c(SIZE(eMesh),iz)/2.0
81 : ENDDO
82 42 : sigma(iContour) = additional_smearing
83 1 : ELSE IF(int_method(shape) == method_maclaurin) THEN
84 :
85 5401 : IF(ANY(ABS(AIMAG(contour)-AIMAG(contour(1)))>1e-12)) THEN
86 0 : CALL juDFT_error("Unsuitable contour for integration", calledby="kkintgr_init")
87 : ENDIF
88 :
89 1 : sigma(iContour) = AIMAG(contour(1))
90 :
91 1 : CALL integration_weights(iContour)%init(.false.,SIZE(eMesh),2*SIZE(contour))
92 :
93 5401 : DO iz = 1, SIZE(contour)
94 : !Next point to the left of the point
95 5400 : n1 = INT((REAL(contour(iz))-eMesh(1))/del) + 1
96 :
97 29165400 : integration_weights(iContour)%data_c(:,iz) = cmplx_0
98 : !Calculate the real part on the same energy points as the imaginary part
99 : !regardless of the contour
100 : !If i is odd skip the odd points and the other way around and use the trapezian method
101 8100 : DO i = MERGE(1,2,MOD(n1,2)==0), SIZE(eMesh), 2
102 14580000 : y = -1/pi_const * 2.0/REAL(n1-i)
103 : IF(i.EQ.1 .OR. i.EQ.2 .OR.&
104 14580000 : i.EQ.SIZE(eMesh) .OR. i.EQ.SIZE(eMesh)-1) y = y/2.0
105 14580000 : integration_weights(iContour)%data_c(i,iz) = y
106 : ENDDO
107 5400 : IF(n1>=1.AND.n1<= SIZE(eMesh)) THEN
108 5400 : integration_weights(iContour)%data_c(n1,iz) = ImagUnit
109 : ENDIF
110 :
111 29165401 : integration_weights(iContour)%data_c(:,iz) = integration_weights(iContour)%data_c(:,iz) * (1.0-(REAL(contour(iz))-(n1-1)*del-eMesh(1))/del)
112 :
113 : ENDDO
114 :
115 5401 : DO izp = SIZE(contour)+1, 2*SIZE(contour)
116 : !Next point to the right of the point
117 5400 : iz=izp-SIZE(contour)
118 5400 : n2 = INT((REAL(contour(iz))-eMesh(1))/del) + 2
119 :
120 29165400 : integration_weights(iContour)%data_c(:,izp) = cmplx_0
121 : !Calculate the real part on the same energy points as the imaginary part
122 : !regardless of the contour
123 : !If i is odd skip the odd points and the other way around and use the trapezian method
124 8100 : DO i = MERGE(1,2,MOD(n2,2)==0), SIZE(eMesh), 2
125 14580000 : y = -1/pi_const * 2.0/REAL(n2-i)
126 : IF(i.EQ.1 .OR. i.EQ.2 .OR.&
127 14580000 : i.EQ.SIZE(eMesh) .OR. i.EQ.SIZE(eMesh)-1) y = y/2.0
128 14580000 : integration_weights(iContour)%data_c(i,izp) = y
129 : ENDDO
130 5400 : IF(n2>=1.AND.n2<=SIZE(eMesh)) THEN
131 5399 : integration_weights(iContour)%data_c(n2,izp) = ImagUnit
132 : ENDIF
133 29165401 : integration_weights(iContour)%data_c(:,izp) = integration_weights(iContour)%data_c(:,izp) * (REAL(contour(iz))-(n2-2)*del-eMesh(1))/del
134 : ENDDO
135 :
136 : ELSE
137 0 : CALL juDFT_error("Not a supported method", calledby="kkintgr_init")
138 : ENDIF
139 :
140 : END SUBROUTINE kkintgr_init
141 :
142 42 : SUBROUTINE kkintgr_free()
143 :
144 : INTEGER :: iContour
145 :
146 42 : IF(.NOT.ALLOCATED(integration_weights)) RETURN
147 :
148 88 : DO iContour = 1, SIZE(integration_weights)
149 88 : IF(integration_weights(iContour)%allocated()) CALL integration_weights(iContour)%free()
150 : ENDDO
151 42 : del = -1.0
152 88 : DEALLOCATE(integration_weights,methods,energy_grid,sigma)
153 :
154 : END SUBROUTINE kkintgr_free
155 :
156 :
157 0 : SUBROUTINE kkintgr_single(im,l_conjg,g,iContour)
158 :
159 : !calculates the Kramer Kronig Transformation on the same contour where the imaginary part was calculated
160 : !Re(G(E+i * delta)) = -1/pi * int_bot^top dE' P(1/(E-E')) * Im(G(E'+i*delta))
161 :
162 : !The dominant source of error for this routine is a insufficiently dense energy mesh on the real axis
163 : !TODO: Some way to estimate the error (maybe search for the sharpest peak and estimate from width)
164 :
165 : COMPLEX, INTENT(IN) :: im(:) !Imaginary part of the green's function on the real axis
166 : LOGICAL, INTENT(IN) :: l_conjg !Switch determines wether we calculate g on the complex conjugate of the contour ez
167 : COMPLEX, INTENT(INOUT) :: g(:) !Green's function on the complex plane
168 : INTEGER, INTENT(IN) :: iContour !Which contour is used
169 :
170 : INTEGER :: ne,nz
171 : COMPLEX, ALLOCATABLE :: smoothed(:)
172 : CHARACTER(len=1) :: transA
173 :
174 0 : IF(.NOT.integration_weights(iContour)%allocated()) THEN
175 : CALL juDFT_error("Integration weights not initialized",&
176 : hint="This is a bug in FLEUR, please report",&
177 0 : calledby="kkintgr")
178 : ENDIF
179 :
180 0 : nz = integration_weights(iContour)%matsize2
181 0 : ne = integration_weights(iContour)%matsize1
182 :
183 0 : smoothed = im
184 0 : IF(ABS(sigma(iContour)).GT.1e-12) THEN
185 0 : CALL timestart('kkintgr: smoothing')
186 0 : SELECT CASE (TRIM(ADJUSTL(smooth_method)))
187 : CASE('lorentzian')
188 0 : CALL lorentzian_smooth(energy_grid,smoothed,sigma(iContour),ne)
189 : CASE('gaussian')
190 0 : CALL smooth(energy_grid,smoothed,sigma(iContour),ne)
191 : CASE DEFAULT
192 : CALL juDFT_error("No valid smooth_method set",&
193 : hint="This is a bug in FLEUR, please report",&
194 0 : calledby="kkintgr")
195 : END SELECT
196 :
197 0 : CALL timestop('kkintgr: smoothing')
198 : ENDIF
199 :
200 0 : transA = 'T'
201 0 : IF(l_conjg) transA = 'C'
202 :
203 0 : IF(methods(iContour)==method_direct) THEN
204 : CALL zgemm(transA,'N',&
205 : nz,1,ne,&
206 : CMPLX(-del/pi_const,0.0),&
207 : integration_weights(iContour)%data_c,ne,&
208 : smoothed,ne,&
209 : cmplx_0,&
210 0 : g,nz)
211 0 : ELSE IF(methods(iContour)==method_maclaurin) THEN
212 :
213 0 : nz = INT(integration_weights(iContour)%matsize2/2)
214 :
215 : CALL zgemm('T','N',&
216 : nz,1,ne,&
217 : cmplx_1,&
218 : integration_weights(iContour)%data_c(:,:nz),ne,&
219 : smoothed,ne,&
220 : cmplx_0,&
221 0 : g,nz)
222 : CALL zgemm('T','N',&
223 : nz,1,ne,&
224 : cmplx_1,&
225 : integration_weights(iContour)%data_c(:,nz+1:),ne,&
226 : smoothed,ne,&
227 : cmplx_1,&
228 0 : g,nz)
229 :
230 0 : IF(l_conjg) g = conjg(g)
231 :
232 : ENDIF
233 :
234 :
235 0 : END SUBROUTINE kkintgr_single
236 :
237 716 : SUBROUTINE kkintgr_mmpmat(im,l_conjg,g,iContour)
238 :
239 : !calculates the Kramer Kronig Transformation on the same contour where the imaginary part was calculated
240 : !Re(G(E+i * delta)) = -1/pi * int_bot^top dE' P(1/(E-E')) * Im(G(E'+i*delta))
241 :
242 : !The dominant source of error for this routine is a insufficiently dense energy mesh on the real axis
243 : !TODO: Some way to estimate the error (maybe search for the sharpest peak and estimate from width)
244 : USE m_smooth
245 : USE m_lorentzian_smooth
246 :
247 : COMPLEX, INTENT(IN) :: im(:,-lmaxU_const:,-lmaxU_const:) !Imaginary part of the green's function on the real axis
248 : LOGICAL, INTENT(IN) :: l_conjg !Switch determines wether we calculate g on the complex conjugate of the contour ez
249 : COMPLEX, INTENT(INOUT) :: g(:,-lmaxU_const:,-lmaxU_const:) !Green's function on the complex plane
250 : INTEGER, INTENT(IN) :: iContour !Which contour is used
251 :
252 : INTEGER :: ne,nz,m,mp
253 716 : COMPLEX, ALLOCATABLE :: smoothed(:,:,:)
254 : CHARACTER(len=1) :: transA
255 :
256 716 : IF(.NOT.integration_weights(iContour)%allocated()) THEN
257 : CALL juDFT_error("Integration weights not initialized",&
258 : hint="This is a bug in FLEUR, please report",&
259 0 : calledby="kkintgr")
260 : ENDIF
261 :
262 716 : nz = integration_weights(iContour)%matsize2
263 716 : ne = integration_weights(iContour)%matsize1
264 :
265 2864 : ALLOCATE(smoothed(SIZE(im,1),-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const))
266 237887528 : smoothed = im
267 716 : IF(ABS(sigma(iContour)).GT.1e-12) THEN
268 4 : CALL timestart('kkintgr: smoothing')
269 4 : SELECT CASE (TRIM(ADJUSTL(smooth_method)))
270 : CASE('lorentzian')
271 : !$OMP parallel do default(none) collapse(2) &
272 : !$OMP shared(energy_grid, smoothed,sigma,ne,iContour) &
273 4 : !$OMP private(m,mp)
274 : DO mp = -lmaxU_const, lmaxU_const
275 : DO m = -lmaxU_const, lmaxU_const
276 : CALL lorentzian_smooth(energy_grid,smoothed(:,m,mp),sigma(iContour),ne)
277 : ENDDO
278 : ENDDO
279 : !$OMP end parallel do
280 : CASE('gaussian')
281 : !$OMP parallel do default(none) collapse(2) &
282 : !$OMP shared(energy_grid, smoothed,sigma,ne,iContour) &
283 0 : !$OMP private(m,mp)
284 : DO mp = -lmaxU_const, lmaxU_const
285 : DO m = -lmaxU_const, lmaxU_const
286 : CALL smooth(energy_grid,smoothed(:,m,mp),sigma(iContour),ne)
287 : ENDDO
288 : ENDDO
289 : !$OMP end parallel do
290 : CASE DEFAULT
291 : CALL juDFT_error("No valid smooth_method set",&
292 : hint="This is a bug in FLEUR, please report",&
293 4 : calledby="kkintgr")
294 : END SELECT
295 :
296 4 : CALL timestop('kkintgr: smoothing')
297 : ENDIF
298 :
299 716 : transA = 'T'
300 716 : IF(l_conjg) transA = 'C'
301 :
302 716 : IF(methods(iContour)==method_direct) THEN
303 : CALL zgemm(transA,'N',&
304 : nz,(2*lmaxU_const+1)**2,ne,&
305 : CMPLX(-del/pi_const,0.0),&
306 : integration_weights(iContour)%data_c,ne,&
307 : smoothed,ne,&
308 : cmplx_0,&
309 712 : g,nz)
310 4 : ELSE IF(methods(iContour)==method_maclaurin) THEN
311 :
312 4 : nz = INT(integration_weights(iContour)%matsize2/2)
313 :
314 : CALL zgemm('T','N',&
315 : nz,(2*lmaxU_const+1)**2,ne,&
316 : cmplx_1,&
317 : integration_weights(iContour)%data_c(:,:nz),ne,&
318 : smoothed,ne,&
319 : cmplx_0,&
320 4 : g,nz)
321 : CALL zgemm('T','N',&
322 : nz,(2*lmaxU_const+1)**2,ne,&
323 : cmplx_1,&
324 : integration_weights(iContour)%data_c(:,nz+1:),ne,&
325 : smoothed,ne,&
326 : cmplx_1,&
327 4 : g,nz)
328 :
329 529316 : IF(l_conjg) g = conjg(g)
330 :
331 : ENDIF
332 :
333 :
334 716 : END SUBROUTINE kkintgr_mmpmat
335 :
336 0 : SUBROUTINE kkintgr_mmpmat_extra1(im,l_conjg,g,iContour)
337 :
338 : !calculates the Kramer Kronig Transformation on the same contour where the imaginary part was calculated
339 : !Re(G(E+i * delta)) = -1/pi * int_bot^top dE' P(1/(E-E')) * Im(G(E'+i*delta))
340 :
341 : !The dominant source of error for this routine is a insufficiently dense energy mesh on the real axis
342 : !TODO: Some way to estimate the error (maybe search for the sharpest peak and estimate from width)
343 : USE m_smooth
344 : USE m_lorentzian_smooth
345 :
346 : COMPLEX, INTENT(IN) :: im(:,-lmaxU_const:,-lmaxU_const:,:) !Imaginary part of the green's function on the real axis
347 : LOGICAL, INTENT(IN) :: l_conjg !Switch determines wether we calculate g on the complex conjugate of the contour ez
348 : COMPLEX, INTENT(INOUT) :: g(:,-lmaxU_const:,-lmaxU_const:,:) !Green's function on the complex plane
349 : INTEGER, INTENT(IN) :: iContour !Which contour is used
350 :
351 : INTEGER :: ne,nz,m,mp,i
352 0 : COMPLEX, ALLOCATABLE :: smoothed(:,:,:,:)
353 : CHARACTER(len=1) :: transA
354 :
355 0 : IF(.NOT.integration_weights(iContour)%allocated()) THEN
356 : CALL juDFT_error("Integration weights not initialized",&
357 : hint="This is a bug in FLEUR, please report",&
358 0 : calledby="kkintgr")
359 : ENDIF
360 :
361 0 : nz = integration_weights(iContour)%matsize2
362 0 : ne = integration_weights(iContour)%matsize1
363 :
364 0 : ALLOCATE(smoothed(SIZE(im,1),-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,SIZE(im,4)))
365 0 : smoothed = im
366 0 : IF(ABS(sigma(iContour)).GT.1e-12) THEN
367 0 : CALL timestart('kkintgr: smoothing')
368 0 : SELECT CASE (TRIM(ADJUSTL(smooth_method)))
369 : CASE('lorentzian')
370 : !$OMP parallel do default(none) collapse(3) &
371 : !$OMP shared(energy_grid, smoothed,sigma,ne,iContour) &
372 0 : !$OMP private(m,mp,i)
373 : DO i = 1, SIZE(smoothed,4)
374 : DO mp = -lmaxU_const, lmaxU_const
375 : DO m = -lmaxU_const, lmaxU_const
376 : CALL lorentzian_smooth(energy_grid,smoothed(:,m,mp,i),sigma(iContour),ne)
377 : ENDDO
378 : ENDDO
379 : ENDDO
380 : !$OMP end parallel do
381 : CASE('gaussian')
382 : !$OMP parallel do default(none) collapse(3) &
383 : !$OMP shared(energy_grid, smoothed,sigma,ne,iContour) &
384 0 : !$OMP private(m,mp)
385 : DO i = 1, SIZE(smoothed,4)
386 : DO mp = -lmaxU_const, lmaxU_const
387 : DO m = -lmaxU_const, lmaxU_const
388 : CALL smooth(energy_grid,smoothed(:,m,mp,i),sigma(iContour),ne)
389 : ENDDO
390 : ENDDO
391 : ENDDO
392 : !$OMP end parallel do
393 : CASE DEFAULT
394 : CALL juDFT_error("No valid smooth_method set",&
395 : hint="This is a bug in FLEUR, please report",&
396 0 : calledby="kkintgr")
397 : END SELECT
398 :
399 0 : CALL timestop('kkintgr: smoothing')
400 : ENDIF
401 :
402 0 : transA = 'T'
403 0 : IF(l_conjg) transA = 'C'
404 :
405 0 : IF(methods(iContour)==method_direct) THEN
406 : CALL zgemm(transA,'N',&
407 : nz,(2*lmaxU_const+1)**2*SIZE(im,4),ne,&
408 : CMPLX(-del/pi_const,0.0),&
409 : integration_weights(iContour)%data_c,ne,&
410 : smoothed,ne,&
411 : cmplx_0,&
412 0 : g,nz)
413 0 : ELSE IF(methods(iContour)==method_maclaurin) THEN
414 :
415 0 : nz = INT(integration_weights(iContour)%matsize2/2)
416 :
417 : CALL zgemm('T','N',&
418 : nz,(2*lmaxU_const+1)**2*SIZE(smoothed,4),ne,&
419 : cmplx_1,&
420 : integration_weights(iContour)%data_c(:,:nz),ne,&
421 : smoothed,ne,&
422 : cmplx_0,&
423 0 : g,nz)
424 : CALL zgemm('T','N',&
425 : nz,(2*lmaxU_const+1)**2*SIZE(smoothed,4),ne,&
426 : cmplx_1,&
427 : integration_weights(iContour)%data_c(:,nz+1:),ne,&
428 : smoothed,ne,&
429 : cmplx_1,&
430 0 : g,nz)
431 :
432 0 : IF(l_conjg) g = conjg(g)
433 :
434 : ENDIF
435 :
436 0 : END SUBROUTINE kkintgr_mmpmat_extra1
437 :
438 : END MODULE m_kkintgr
|
