Branch data Line data Source code
1 : :
2 : : #if defined HAVE_CONFIG_H
3 : : #include "config.h"
4 : : #endif
5 : : !!****m* ABINIT/m_GreenHyb
6 : : !! NAME
7 : : !! m_GreenHyb
8 : : !!
9 : : !! FUNCTION
10 : : !! Manage a green function for one orbital
11 : : !!
12 : : !! COPYRIGHT
13 : : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
14 : : !! This file is distributed under the terms of the
15 : : !! GNU General Public License, see ~abinit/COPYING
16 : : !! or http://www.gnu.org/copyleft/gpl.txt .
17 : : !!
18 : : !! NOTES
19 : : !!
20 : : !! SOURCE
21 : :
22 : : #include "defs.h"
23 : : MODULE m_GreenHyb
24 : : USE m_Global
25 : : USE m_MatrixHyb
26 : : USE m_Vector
27 : : USE m_VectorInt
28 : : USE m_ListCdagC
29 : : USE m_MapHyb
30 : : USE_MPI
31 : :
32 : : IMPLICIT NONE
33 : :
34 : : !!***
35 : :
36 : : PRIVATE
37 : :
38 : : !!****t* m_GreenHyb/GreenHyb
39 : : !! NAME
40 : : !! GreenHyb
41 : : !!
42 : : !! FUNCTION
43 : : !! This structured datatype contains the necessary data
44 : : !!
45 : : !! COPYRIGHT
46 : : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
47 : : !! This file is distributed under the terms of the
48 : : !! GNU General Public License, see ~abinit/COPYING
49 : : !! or http://www.gnu.org/copyleft/gpl.txt .
50 : : !!
51 : : !! SOURCE
52 : :
53 : : TYPE, PUBLIC :: GreenHyb
54 : : LOGICAL _PRIVATE :: set = .FALSE.
55 : : LOGICAL _PRIVATE :: setT = .FALSE.
56 : : LOGICAL _PRIVATE :: setW = .FALSE.
57 : : LOGICAL _PRIVATE :: have_MPI = .FALSE.
58 : : INTEGER _PRIVATE :: setMk = 0
59 : : INTEGER _PRIVATE :: samples
60 : : INTEGER _PRIVATE :: measurements
61 : : INTEGER :: factor
62 : : INTEGER _PRIVATE :: MY_COMM
63 : : INTEGER _PRIVATE :: size
64 : : INTEGER _PRIVATE :: rank
65 : : INTEGER _PRIVATE :: Wmax
66 : : INTEGER _PRIVATE :: iTech
67 : : DOUBLE PRECISION _PRIVATE :: beta
68 : : DOUBLE PRECISION _PRIVATE :: inv_beta
69 : : DOUBLE PRECISION _PRIVATE :: delta_t
70 : : DOUBLE PRECISION _PRIVATE :: inv_dt
71 : : DOUBLE PRECISION, ALLOCATABLE , DIMENSION(:) :: oper
72 : : DOUBLE PRECISION, ALLOCATABLE , DIMENSION(:) _PRIVATE :: omega
73 : : DOUBLE PRECISION , DIMENSION(1:3) _PRIVATE :: Mk
74 : : COMPLEX(KIND=8) , ALLOCATABLE, DIMENSION(:) _PRIVATE :: oper_w
75 : : COMPLEX(KIND=8) , ALLOCATABLE, DIMENSION(:) _PRIVATE :: oper_w_old
76 : : TYPE(MapHyb) :: this
77 : : END TYPE GreenHyb
78 : : !!***
79 : :
80 : : PUBLIC :: GreenHyb_init
81 : : PUBLIC :: GreenHyb_clear
82 : : PUBLIC :: GreenHyb_reset
83 : : PUBLIC :: GreenHyb_setOperW
84 : : PUBLIC :: GreenHyb_measHybrid
85 : : PUBLIC :: GreenHyb_getHybrid
86 : : PUBLIC :: GreenHyb_setN
87 : : PUBLIC :: GreenHyb_setMuD1
88 : : PUBLIC :: GreenHyb_setMoments
89 : : PUBLIC :: GreenHyb_backFourier
90 : : PUBLIC :: GreenHyb_forFourier
91 : : PUBLIC :: GreenHyb_print
92 : : PUBLIC :: GreenHyb_destroy
93 : :
94 : : CONTAINS
95 : : !!***
96 : :
97 : : !!****f* ABINIT/m_GreenHyb/GreenHyb_init
98 : : !! NAME
99 : : !! GreenHyb_init
100 : : !!
101 : : !! FUNCTION
102 : : !! Initialize and allocate
103 : : !!
104 : : !! COPYRIGHT
105 : : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
106 : : !! This file is distributed under the terms of the
107 : : !! GNU General Public License, see ~abinit/COPYING
108 : : !! or http://www.gnu.org/copyleft/gpl.txt .
109 : : !!
110 : : !! INPUTS
111 : : !! this=Green
112 : : !! samples=imaginary time slices
113 : : !! beta=inverse temperature
114 : : !! iTech=SHOULD NOT BE USED => BUGGY
115 : : !! MY_COMM=mpi_communicator
116 : : !!
117 : : !! OUTPUT
118 : : !!
119 : : !! SIDE EFFECTS
120 : : !!
121 : : !! NOTES
122 : : !!
123 : : !! SOURCE
124 : :
125 : 1125 : SUBROUTINE GreenHyb_init(this, samples, beta,iTech,MY_COMM)
126 : :
127 : :
128 : : #ifdef HAVE_MPI1
129 : : include 'mpif.h'
130 : : #endif
131 : : !Arguments ------------------------------------
132 : : TYPE(GreenHyb) , INTENT(INOUT) :: this
133 : : INTEGER , INTENT(IN ) :: samples
134 : : DOUBLE PRECISION, INTENT(IN ) :: beta
135 : : !INTEGER , INTENT(IN ) :: Wmax
136 : : INTEGER, OPTIONAL, INTENT(IN ) :: iTech
137 : : INTEGER, OPTIONAL, INTENT(IN ) :: MY_COMM
138 : : !Local variables ------------------------------
139 : : INTEGER :: sp1
140 : : DOUBLE PRECISION :: dt
141 : : #ifdef HAVE_MPI
142 : : INTEGER :: ierr
143 : : #endif
144 : :
145 [ + - ]: 1125 : IF ( PRESENT(MY_COMM)) THEN
146 : : #ifdef HAVE_MPI
147 : 1125 : this%have_MPI = .TRUE.
148 : 1125 : this%MY_COMM = MY_COMM
149 : 1125 : CALL MPI_Comm_rank(this%MY_COMM, this%rank, ierr)
150 : 1125 : CALL MPI_Comm_size(this%MY_COMM, this%size, ierr)
151 : : #else
152 : : CALL WARN("GreenHyb_init : MPI is not used ")
153 : : this%have_MPI = .FALSE.
154 : : this%MY_COMM = -1
155 : : this%rank = 0
156 : : this%size = 1
157 : : #endif
158 : : ELSE
159 : 0 : this%have_MPI = .FALSE.
160 : 0 : this%MY_COMM = -1
161 : 0 : this%rank = 0
162 : 0 : this%size = 1
163 : : END IF
164 : :
165 : 1125 : sp1 = samples + 1
166 : 1125 : this%samples = sp1
167 : 1125 : this%measurements = 0
168 : 1125 : this%beta = beta
169 : 1125 : this%inv_beta = 1.d0 / beta
170 : 1125 : this%inv_dt = DBLE(samples) * this%inv_beta
171 : 1125 : dt = 1.d0 / this%inv_dt
172 : 1125 : this%delta_t = dt
173 : : !this%Wmax = Wmax
174 : 1125 : this%Wmax = -1
175 [ - + ]: 1125 : FREEIF(this%oper)
176 [ - + + - : 3375 : MALLOC(this%oper,(1:sp1))
+ - + - -
+ - + -
+ ]
177 : : ! If we want to measure in frequences
178 : : ! let assume we first have "samples" frequences
179 [ + + ]: 1125 : IF ( PRESENT(iTech) ) THEN
180 : 1076 : this%iTech = iTech
181 : 1076 : SELECT CASE (this%iTech)
182 : : CASE (GREENHYB_TAU) ! omega
183 : 1076 : this%iTech = GREENHYB_TAU
184 : : CASE (GREENHYB_OMEGA) ! omega
185 : 0 : this%Wmax = samples
186 [ # # ]: 0 : FREEIF(this%oper_w)
187 [ # # # # : 0 : MALLOC(this%oper_w,(1:this%Wmax))
# # # # #
# # # #
# ]
188 [ # # ]: 0 : FREEIF(this%oper_w_old)
189 [ # # # # : 0 : MALLOC(this%oper_w_old,(1:this%Wmax))
# # # # #
# # # #
# ]
190 [ # # ]: 0 : this%oper_w = CMPLX(0.d0,0.d0,8)
191 [ # # ]: 0 : this%oper_w_old = CMPLX(0.d0,0.d0,8)
192 [ # # ]: 0 : FREEIF(this%omega)
193 [ # # # # : 0 : MALLOC(this%omega,(1:this%Wmax))
# # # # #
# # # #
# ]
194 [ + - - - : 1076 : this%omega = (/ ((2.d0 * DBLE(sp1) - 1.d0)*ACOS(-1.d0)*this%inv_beta, sp1=1, this%Wmax) /)
- - - - -
- - - - -
- - - -
- ]
195 : : END SELECT
196 : : ELSE
197 : 49 : this%iTech = GREENHYB_TAU
198 : : END IF
199 : : ! end if
200 : : !CALL Vector_init(this%oper_old,10000)
201 : : !CALL VectorInt_init(this%index_old,10000)
202 : 1125 : CALL MapHyb_init(this%this,10000)
203 : :
204 [ + + ]: 341100 : this%oper = 0.d0
205 : 1125 : this%set = .TRUE.
206 : 1125 : this%factor = 1
207 : 1125 : this%setMk = 0
208 [ + + ]: 4500 : this%Mk = 0.d0
209 : 1125 : END SUBROUTINE GreenHyb_init
210 : : !!***
211 : :
212 : : !!****f* ABINIT/m_GreenHyb/GreenHyb_clear
213 : : !! NAME
214 : : !! GreenHyb_clear
215 : : !!
216 : : !! FUNCTION
217 : : !! clear green function
218 : : !!
219 : : !! COPYRIGHT
220 : : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
221 : : !! This file is distributed under the terms of the
222 : : !! GNU General Public License, see ~abinit/COPYING
223 : : !! or http://www.gnu.org/copyleft/gpl.txt .
224 : : !!
225 : : !! INPUTS
226 : : !! this=Green
227 : : !!
228 : : !! OUTPUT
229 : : !!
230 : : !! SIDE EFFECTS
231 : : !!
232 : : !! NOTES
233 : : !!
234 : : !! SOURCE
235 : :
236 : 1554 : SUBROUTINE GreenHyb_clear(this)
237 : :
238 : : !Arguments ------------------------------------
239 : : TYPE(GreenHyb) , INTENT(INOUT) :: this
240 : :
241 : : !CALL Vector_clear(this%oper_old)
242 : : !CALL VectorInt_clear(this%index_old)
243 : 1554 : CALL MapHyb_clear(this%this)
244 : 1554 : this%measurements = 0
245 [ + - ]: 1554 : IF ( ALLOCATED(this%oper) ) &
246 [ + + ]: 458808 : this%oper = 0.d0
247 [ - + ]: 1554 : IF ( this%iTech .EQ. GREENHYB_OMEGA ) THEN
248 [ # # ]: 0 : IF ( ALLOCATED(this%oper_w) ) &
249 [ # # ]: 0 : this%oper_w = CMPLX(0.d0,0.d0,8)
250 [ # # ]: 0 : IF ( ALLOCATED(this%oper_w_old) ) &
251 [ # # ]: 0 : this%oper_w_old = CMPLX(0.d0,0.d0,8)
252 : : END IF
253 : 1554 : this%factor = 0
254 : 1554 : END SUBROUTINE GreenHyb_clear
255 : : !!***
256 : :
257 : : !!****f* ABINIT/m_GreenHyb/GreenHyb_reset
258 : : !! NAME
259 : : !! GreenHyb_reset
260 : : !!
261 : : !! FUNCTION
262 : : !! reset green function
263 : : !!
264 : : !! COPYRIGHT
265 : : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
266 : : !! This file is distributed under the terms of the
267 : : !! GNU General Public License, see ~abinit/COPYING
268 : : !! or http://www.gnu.org/copyleft/gpl.txt .
269 : : !!
270 : : !! INPUTS
271 : : !! this=Green
272 : : !!
273 : : !! OUTPUT
274 : : !!
275 : : !! SIDE EFFECTS
276 : : !!
277 : : !! NOTES
278 : : !!
279 : : !! SOURCE
280 : :
281 : 518 : SUBROUTINE GreenHyb_reset(this)
282 : :
283 : : !Arguments ------------------------------------
284 : : TYPE(GreenHyb) , INTENT(INOUT) :: this
285 : :
286 : 518 : CALL GreenHyb_clear(this)
287 : 518 : this%setMk = 0
288 [ + + ]: 2072 : this%Mk = 0.d0
289 : 518 : this%setT = .FALSE.
290 : 518 : this%setW = .FALSE.
291 : 518 : END SUBROUTINE GreenHyb_reset
292 : : !!***
293 : :
294 : : !!****f* ABINIT/m_GreenHyb/GreenHyb_setOperW
295 : : !! NAME
296 : : !! GreenHyb_setOperW
297 : : !!
298 : : !! FUNCTION
299 : : !! set Green function in frequencies
300 : : !!
301 : : !! COPYRIGHT
302 : : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
303 : : !! This file is distributed under the terms of the
304 : : !! GNU General Public License, see ~abinit/COPYING
305 : : !! or http://www.gnu.org/copyleft/gpl.txt .
306 : : !!
307 : : !! INPUTS
308 : : !! this=Green
309 : : !! Gomega=Input values
310 : : !!
311 : : !! OUTPUT
312 : : !!
313 : : !! SIDE EFFECTS
314 : : !!
315 : : !! NOTES
316 : : !!
317 : : !! SOURCE
318 : :
319 [ - + ]: 518 : SUBROUTINE GreenHyb_setOperW(this, Gomega)
320 : :
321 : : !Arguments ------------------------------------
322 : : TYPE(GreenHyb) , INTENT(INOUT) :: this
323 : : COMPLEX(KIND=8), DIMENSION(:), INTENT(IN ) :: Gomega
324 : : !Loval variables ------------------------------
325 : : INTEGER :: tail
326 : :
327 : 518 : tail = SIZE(Gomega)
328 [ - + ]: 518 : IF ( .NOT. this%set ) &
329 : 0 : CALL ERROR("GreenHyb_setOperW : Uninitialized GreenHyb structure")
330 [ + + ]: 518 : IF ( ALLOCATED(this%oper_w) ) THEN
331 [ - + ]: 469 : IF ( SIZE(this%oper_w) .NE. tail ) THEN
332 : 0 : FREE(this%oper_w)
333 [ # # # # : 0 : MALLOC(this%oper_w,(1:tail))
# # # # #
# # # ]
334 : : END IF
335 : : ELSE
336 [ - + + - : 147 : MALLOC(this%oper_w,(1:tail))
+ - + - -
+ - + -
+ ]
337 : : END IF
338 [ + + ]: 152418 : this%oper_w(:) = Gomega(:)
339 : 518 : this%Wmax = tail
340 : 518 : this%setW = .TRUE.
341 : 518 : END SUBROUTINE GreenHyb_setOperW
342 : : !!***
343 : :
344 : : !!****f* ABINIT/m_GreenHyb/GreenHyb_measHybrid
345 : : !! NAME
346 : : !! GreenHyb_measHybrid
347 : : !!
348 : : !! FUNCTION
349 : : !! Measure Green function
350 : : !!
351 : : !! COPYRIGHT
352 : : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
353 : : !! This file is distributed under the terms of the
354 : : !! GNU General Public License, see ~abinit/COPYING
355 : : !! or http://www.gnu.org/copyleft/gpl.txt .
356 : : !!
357 : : !! INPUTS
358 : : !! this=Green
359 : : !! Mthis=M this for the current flavor
360 : : !! ListCdagC_1=list of all creator and annhilator operators
361 : : !! updated=should we accumulate or not
362 : : !!
363 : : !! OUTPUT
364 : : !!
365 : : !! SIDE EFFECTS
366 : : !!
367 : : !! NOTES
368 : : !!
369 : : !! SOURCE
370 : :
371 : 1022380518 : SUBROUTINE GreenHyb_measHybrid(this, Mthis, ListCdagC_1, updated)
372 : :
373 : : !Arguments ------------------------------------
374 : : TYPE(GreenHyb) , INTENT(INOUT) :: this
375 : : TYPE(MatrixHyb) , INTENT(IN ) :: Mthis
376 : : TYPE(ListcdagC), INTENT(IN ) :: ListCdagC_1
377 : : LOGICAL , INTENT(IN ) :: updated
378 : : !Local variables ------------------------------
379 : : INTEGER :: iC
380 : : INTEGER :: iCdag
381 : : INTEGER :: tail
382 : : !INTEGER :: index
383 : : INTEGER :: idx_old
384 : : INTEGER :: old_size
385 : : INTEGER :: omegaSamples
386 : : INTEGER :: iomega
387 : : DOUBLE PRECISION :: pi_invBeta
388 : : DOUBLE PRECISION :: mbeta_two
389 : : DOUBLE PRECISION :: beta
390 : : DOUBLE PRECISION :: beta_tc
391 : : DOUBLE PRECISION :: tcbeta_tc
392 : : DOUBLE PRECISION :: inv_dt
393 : : DOUBLE PRECISION :: tC
394 : : DOUBLE PRECISION :: tCdag
395 : : DOUBLE PRECISION :: time
396 : : DOUBLE PRECISION :: signe
397 : : DOUBLE PRECISION :: argument
398 : : !DOUBLE PRECISION :: taupi_invbeta
399 : : !COMPLEX(KIND=8) :: cargument
400 : : !COMPLEX(2*8) :: base_exp
401 : : !COMPLEX(2*8) :: increm_exp
402 : :
403 [ - + ]: 1022380518 : IF ( this%set .EQV. .FALSE. ) &
404 : 0 : CALL ERROR("GreenHyb_measHybrid : green operator not set ")
405 : :
406 : 1022380518 : tail = ListCdagC_1%tail
407 [ - + ]: 1022380518 : IF ( tail .NE. Mthis%tail ) &
408 : 0 : CALL ERROR("GreenHyb_measHybrid : ListCdagC & M unconsistent ")
409 : :
410 [ + + ]: 1022380518 : IF ( updated .EQV. .TRUE. ) THEN ! NEW change in the configuration
411 : : ! FIXME SHOULD be much more faster
412 : :
413 : 158970444 : old_size = this%this%tail
414 : 317940888 : SELECT CASE(this%iTech)
415 : : CASE (GREENHYB_TAU)
416 : 158970444 : argument = DBLE(this%factor)
417 [ + + ]: 3591438449 : DO iC = 1, old_size
418 : 3591438449 : this%oper(this%this%listINT(iC)) = this%oper(this%this%listINT(iC)) + this%this%listDBLE(iC) * argument
419 : : END DO
420 : 158970444 : this%measurements = this%measurements + this%factor
421 : :
422 : 158970444 : CALL MapHyb_setSize(this%this,tail*tail)
423 : 158970444 : this%factor = 1
424 : 158970444 : idx_old = 0
425 : 158970444 : beta = this%beta
426 : 158970444 : mbeta_two = -(beta*0.5d0)
427 : 158970444 : inv_dt = this%inv_dt
428 : : ! WARNING time is not the time but just a temporary variable.
429 : : ! Index Time has been calculated previously and is in mat_tau
430 [ + + ]: 715187494 : DO iC = 1, tail
431 : 556217050 : tC = ListCdagC_1%list(iC,C_)
432 : 556217050 : beta_tc = beta - tC
433 : 556217050 : tcbeta_tc = tC * beta_tc
434 [ + + ]: 4147695554 : DO iCdag = 1, tail
435 : 3432508060 : tCdag = ListCdagC_1%list(iCdag,Cdag_)
436 : 3432508060 : time = tcbeta_tc - tCdag*beta_tc
437 : :
438 : : !signe = SIGN(1.d0,time)
439 : : !time = time + (signe-1.d0)*mbeta_two
440 : :
441 : : !signe = signe * SIGN(1.d0,beta-tC)
442 : :
443 : : !signe = SIGN(1.d0,time) * SIGN(1.d0,beta-tC)
444 : 3432508060 : signe = SIGN(1.d0,time)
445 : :
446 : 3432508060 : argument = signe*Mthis%mat(iCdag,iC)
447 : :
448 : : !index = INT( ( time * inv_dt ) + 1.5d0 )
449 : : !IF (index .NE. Mthis%mat_tau(iCdag,iC)) THEN
450 : : ! WRITE(*,*) index, Mthis%mat_tau(iCdag,iC)
451 : : !! CALL ERROR("Plantage")
452 : : !END IF
453 : :
454 : 3432508060 : idx_old = idx_old + 1
455 : 3432508060 : this%this%listDBLE(idx_old) = argument
456 : : !this%this%listINT(idx_old) = index
457 : 3988725110 : this%this%listINT(idx_old) = Mthis%mat_tau(iCdag,iC)
458 : : END DO
459 : : END DO
460 : : CASE (GREENHYB_OMEGA)
461 : 0 : omegaSamples = this%Wmax
462 : 0 : argument = DBLE(this%factor)
463 [ # # ]: 0 : DO iomega = 1, omegaSamples
464 : 0 : this%oper_w(iomega) = this%oper_w(iomega) + this%oper_w_old(iomega) * argument
465 : : END DO
466 : 0 : this%measurements = this%measurements + this%factor
467 : :
468 : 0 : this%factor = 1
469 : 0 : beta = this%beta
470 : 0 : mbeta_two = -(beta*0.5d0)
471 : 0 : pi_invBeta = ACOS(-1.d0)/beta
472 [ + - - - : 158970444 : DO iC = 1, tail
- ]
473 : 0 : tC = ListCdagC_1%list(iC,C_)
474 [ # # ]: 0 : DO iCdag = 1, tail
475 : 0 : tCdag = ListCdagC_1%list(iCdag,Cdag_)
476 : 0 : time = tC - tCdag
477 : :
478 : 0 : signe = SIGN(1.d0,time)
479 : 0 : time = time + (signe-1.d0)*mbeta_two
480 : 0 : signe = signe * SIGN(1.d0,beta-tC)
481 : 0 : argument = signe*Mthis%mat(iCdag,iC)
482 : :
483 [ # # ]: 0 : DO iomega = 1, omegaSamples
484 : : !this%oper_w_old(iomega) = Mthis%mat_tau(iCdag,iC)*CMPLX(0.d0,argument)
485 : 0 : this%oper_w_old(iomega) = EXP(CMPLX(0.d0,this%omega(iomega)*time,8))*CMPLX(0.d0,argument,8)
486 : : END DO
487 : : END DO
488 : : END DO
489 : : END SELECT
490 : : ELSE
491 : 863410074 : this%factor = this%factor + 1
492 : : END IF
493 : 1022380518 : END SUBROUTINE GreenHyb_measHybrid
494 : : !!***
495 : :
496 : : !!****f* ABINIT/m_GreenHyb/GreenHyb_getHybrid
497 : : !! NAME
498 : : !! GreenHyb_getHybrid
499 : : !!
500 : : !! FUNCTION
501 : : !! reduce green function
502 : : !!
503 : : !! COPYRIGHT
504 : : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
505 : : !! This file is distributed under the terms of the
506 : : !! GNU General Public License, see ~abinit/COPYING
507 : : !! or http://www.gnu.org/copyleft/gpl.txt .
508 : : !!
509 : : !! INPUTS
510 : : !! this=Green
511 : : !!
512 : : !! OUTPUT
513 : : !!
514 : : !! SIDE EFFECTS
515 : : !!
516 : : !! NOTES
517 : : !!
518 : : !! SOURCE
519 : :
520 : 518 : SUBROUTINE GreenHyb_getHybrid(this)
521 : :
522 : : !Arguments ------------------------------------
523 : : TYPE(GreenHyb), INTENT(INOUT) :: this
524 : :
525 [ - + ]: 518 : IF ( this%set .EQV. .FALSE. ) &
526 : 0 : CALL ERROR("GreenHyb_getHybrid : green operator not set ")
527 : :
528 : 1036 : SELECT CASE(this%iTech)
529 : : CASE (GREENHYB_TAU)
530 [ + + ]: 152936 : this%oper = -(this%oper * this%inv_beta) / (DBLE(this%measurements) * this%delta_t)
531 : 518 : this%setT = .TRUE.
532 : : CASE (GREENHYB_OMEGA)
533 [ # # ]: 0 : this%oper_w = -(this%oper_w * this%inv_beta) / (DBLE(this%measurements) * this%delta_t)
534 : 0 : this%setW = .TRUE.
535 [ + - - ]: 518 : CALL GreenHyb_backFourier(this)
536 : : END SELECT
537 : :
538 : 518 : END SUBROUTINE GreenHyb_getHybrid
539 : : !!***
540 : :
541 : : !!****f* ABINIT/m_GreenHyb/GreenHyb_setN
542 : : !! NAME
543 : : !! GreenHyb_setN
544 : : !!
545 : : !! FUNCTION
546 : : !! impose number of electrons for this flavor
547 : : !!
548 : : !! COPYRIGHT
549 : : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
550 : : !! This file is distributed under the terms of the
551 : : !! GNU General Public License, see ~abinit/COPYING
552 : : !! or http://www.gnu.org/copyleft/gpl.txt .
553 : : !!
554 : : !! INPUTS
555 : : !! this=Green
556 : : !! N=number of electrons
557 : : !!
558 : : !! OUTPUT
559 : : !!
560 : : !! SIDE EFFECTS
561 : : !!
562 : : !! NOTES
563 : : !!
564 : : !! SOURCE
565 : :
566 : 518 : SUBROUTINE GreenHyb_setN(this,N)
567 : :
568 : : !Arguments ------------------------------------
569 : : TYPE(GreenHyb) , INTENT(INOUT) :: this
570 : : DOUBLE PRECISION, INTENT(IN ) :: N
571 : :
572 [ - + ]: 518 : IF ( this%set .EQV. .FALSE. ) &
573 : 0 : CALL ERROR("GreenHyb_setN: green this%operator not set ")
574 : 518 : this%oper(1) = N - 1.d0
575 : 518 : this%oper(this%samples) = - N
576 : 518 : END SUBROUTINE GreenHyb_setN
577 : : !!***
578 : :
579 : : !!****f* ABINIT/m_GreenHyb/GreenHyb_setMuD1
580 : : !! NAME
581 : : !! GreenHyb_setMuD1
582 : : !!
583 : : !! FUNCTION
584 : : !! Set first moments for G
585 : : !!
586 : : !! COPYRIGHT
587 : : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
588 : : !! This file is distributed under the terms of the
589 : : !! GNU General Public License, see ~abinit/COPYING
590 : : !! or http://www.gnu.org/copyleft/gpl.txt .
591 : : !!
592 : : !! INPUTS
593 : : !! this=Green
594 : : !! mu=energy level (irrespectige with fermi level)
595 : : !! d1=firt moment of hybridization function
596 : : !!
597 : : !! OUTPUT
598 : : !! argout(sizeout)=description
599 : : !!
600 : : !! SIDE EFFECTS
601 : : !!
602 : : !! NOTES
603 : : !!
604 : : !! SOURCE
605 : :
606 : 518 : SUBROUTINE GreenHyb_setMuD1(this,mu,d1)
607 : :
608 : : !Arguments ------------------------------------
609 : : TYPE(GreenHyb) , INTENT(INOUT) :: this
610 : : DOUBLE PRECISION, INTENT(IN ) :: mu
611 : : DOUBLE PRECISION, INTENT(IN ) :: d1
612 : :
613 : 518 : this%Mk(3) = -d1-(mu*mu)
614 : 518 : this%Mk(2) = -mu
615 : 518 : this%setMk = this%setMk + 1
616 : 518 : END SUBROUTINE GreenHyb_setMuD1
617 : : !!***
618 : :
619 : : !!****f* ABINIT/m_GreenHyb/GreenHyb_setMoments
620 : : !! NAME
621 : : !! GreenHyb_setMoments
622 : : !!
623 : : !! FUNCTION
624 : : !! Compute full moments
625 : : !!
626 : : !! COPYRIGHT
627 : : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
628 : : !! This file is distributed under the terms of the
629 : : !! GNU General Public License, see ~abinit/COPYING
630 : : !! or http://www.gnu.org/copyleft/gpl.txt .
631 : : !!
632 : : !! INPUTS
633 : : !! this=Green
634 : : !! u1=interaction energi like
635 : : !! u2=idem order2
636 : : !!
637 : : !! OUTPUT
638 : : !!
639 : : !! SIDE EFFECTS
640 : : !!
641 : : !! NOTES
642 : : !!
643 : : !! SOURCE
644 : :
645 : 368 : SUBROUTINE GreenHyb_setMoments(this,u1,u2)
646 : :
647 : : !Arguments ------------------------------------
648 : : TYPE(GreenHyb) , INTENT(INOUT) :: this
649 : : DOUBLE PRECISION, INTENT(IN ) :: u1
650 : : DOUBLE PRECISION, INTENT(IN ) :: u2
651 : :
652 : 368 : this%Mk(1) = -1.d0
653 : 368 : this%Mk(3) = this%Mk(3) - 2.d0*(this%Mk(2)*u1)
654 : 368 : this%Mk(2) = this%Mk(2) + u1
655 : 368 : this%Mk(3) = this%Mk(3) - u2
656 : :
657 : 368 : this%setMk = this%setMk + 1
658 : :
659 : 368 : END SUBROUTINE GreenHyb_setMoments
660 : : !!***
661 : :
662 : : !!****f* ABINIT/m_GreenHyb/GreenHyb_backFourier
663 : : !! NAME
664 : : !! GreenHyb_backFourier
665 : : !!
666 : : !! FUNCTION
667 : : !! perform back fourier transform
668 : : !!
669 : : !! COPYRIGHT
670 : : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
671 : : !! This file is distributed under the terms of the
672 : : !! GNU General Public License, see ~abinit/COPYING
673 : : !! or http://www.gnu.org/copyleft/gpl.txt .
674 : : !!
675 : : !! INPUTS
676 : : !! this=Green
677 : : !! dvgc=divergence parameter
678 : : !!
679 : : !! OUTPUT
680 : : !!
681 : : !! SIDE EFFECTS
682 : : !!
683 : : !! NOTES
684 : : !!
685 : : !! SOURCE
686 : :
687 : 518 : SUBROUTINE GreenHyb_backFourier(this,dvgc)
688 : :
689 : :
690 : : #ifdef HAVE_MPI1
691 : : include 'mpif.h'
692 : : #endif
693 : : !Arguments ------------------------------------
694 : : TYPE(GreenHyb) , INTENT(INOUT) :: this
695 : : DOUBLE PRECISION, OPTIONAL, INTENT(IN ) :: dvgc
696 : : !Local variables ------------------------------
697 : : INTEGER :: itau
698 : : INTEGER :: iomega
699 : : INTEGER :: omegaSamples
700 : : INTEGER :: tauSamples
701 : : INTEGER :: tauBegin
702 : : INTEGER :: tauEnd
703 : : INTEGER :: delta
704 : : INTEGER :: residu
705 : 518 : INTEGER, ALLOCATABLE, DIMENSION(:) :: counts
706 : 518 : INTEGER, ALLOCATABLE, DIMENSION(:) :: displs
707 : : DOUBLE PRECISION :: A ! Correction factor
708 : : DOUBLE PRECISION :: inv_beta
709 : : DOUBLE PRECISION :: pi_invBeta
710 : : DOUBLE PRECISION :: two_invBeta
711 : : DOUBLE PRECISION :: minusDt
712 : : DOUBLE PRECISION :: minusOmegaTau
713 : : DOUBLE PRECISION :: omega
714 : : DOUBLE PRECISION :: minusTau
715 : : DOUBLE PRECISION :: sumTerm
716 : : DOUBLE PRECISION :: pi
717 : : DOUBLE PRECISION :: twoPi
718 : : DOUBLE PRECISION :: correction
719 : 518 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: Domega
720 : 518 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: A_omega
721 : :
722 : : #if defined HAVE_MPI && !defined HAVE_MPI2_INPLACE
723 : : INTEGER :: my_count
724 : : DOUBLE PRECISION, ALLOCATABLE , DIMENSION(:) :: oper_buf
725 : : #endif
726 : :
727 [ - + ]: 518 : IF ( this%set .EQV. .FALSE. ) &
728 : 0 : CALL ERROR("GreenHyb_backFourier : Uninitialized GreenHyb structure")
729 [ - + ]: 518 : IF ( this%setW .EQV. .FALSE. ) &
730 : 0 : CALL ERROR("GreenHyb_backFourier : no G(iw)")
731 : :
732 : 518 : inv_beta = this%inv_beta
733 : 518 : two_invBeta = 2.d0 * inv_beta
734 : 518 : minusDt = - this%delta_t
735 : 518 : omegaSamples = this%Wmax
736 : 518 : tauSamples = this%samples-1
737 : :
738 [ + + ]: 152936 : this%oper = 0.d0
739 : :
740 : 518 : pi = ACOS(-1.d0)
741 : 518 : twoPi = 2.d0 * pi
742 : 518 : pi_invBeta = pi * inv_beta
743 : :
744 [ - + ]: 518 : IF ( PRESENT(dvgc) ) THEN
745 : 0 : A = dvgc
746 : : ELSE
747 : : A = AIMAG(this%oper_w(omegaSamples)) & ! A = \lim_\infty G(w)
748 : 518 : *(2.d0*DBLE(omegaSamples)-1.d0) * pi_invBeta
749 : : END IF
750 : :
751 : 518 : correction = A*0.5d0
752 : :
753 [ - + + - : 1554 : MALLOC(Domega,(1:omegaSamples))
+ - + - -
+ - + ]
754 [ - + + - : 1036 : MALLOC(A_omega,(1:omegaSamples))
- + - + ]
755 [ - + + + : 456736 : Domega = (/ ((2.d0 * DBLE(iomega) - 1.d0)*pi_invbeta, iomega=1, omegaSamples) /)
+ - + + +
- - + - -
+ + ]
756 [ + - - + : 152936 : A_omega = A / Domega
- - + + ]
757 [ + - ]: 518 : IF (this%have_MPI .EQV. .TRUE.) THEN
758 : 518 : delta = tauSamples / this%size
759 : 518 : residu = tauSamples - this%size*delta
760 [ + + ]: 518 : IF ( this%rank .LT. this%size - residu ) THEN
761 : 462 : tauBegin = 1 + this%rank*delta
762 : 462 : tauEnd = (this%rank + 1)*delta
763 : : ELSE
764 : : ! tauBegin = (this%size-residu)*delta + 1 + (this%rank-this%size+residu)*(delta+1)
765 : 56 : tauBegin = 1 + this%rank*(delta + 1) -this%size + residu
766 : 56 : tauEnd = tauBegin + delta
767 : : END IF
768 [ - + + - : 1554 : MALLOC(counts,(1:this%size))
+ - + - -
+ - + ]
769 [ - + + - : 1036 : MALLOC(displs,(1:this%size))
- + - + ]
770 : : counts = (/ (delta, iTau=1, this%size-residu), &
771 [ - + + + : 5758 : (delta+1, iTau=this%size-residu+1, this%size) /)
- + + + +
- + + + -
- + - - +
+ ]
772 : 518 : displs(1)=0
773 [ + + ]: 1574 : DO iTau = 2, this%size
774 : 1574 : displs(iTau) = displs(iTau-1) + counts (iTau-1)
775 : : END DO
776 : : ELSE
777 : : tauBegin = 1
778 : : tauEnd = tauSamples
779 : : END IF
780 [ + + ]: 58218 : DO iTau = tauBegin, tauEnd
781 : 57700 : minusTau = DBLE(itau -1) * minusDt
782 [ + + ]: 25442700 : DO iomega = 1, omegaSamples
783 : 25385000 : omega = Domega(iomega)
784 : 25385000 : minusOmegaTau = MOD(omega*minusTau, TwoPi)
785 : : sumTerm = REAL(( this%oper_w(iomega) - CMPLX(0.d0, A_omega(iomega),8) ) &
786 : 25385000 : * EXP( CMPLX(0.d0, minusOmegaTau, 8)))
787 : 25442700 : this%oper(itau) = this%oper(itau) + sumTerm
788 : : END DO
789 : 58218 : this%oper(itau) = correction + two_invBeta*this%oper(itau)
790 : : END DO
791 [ + - ]: 518 : IF ( this%have_MPI .EQV. .TRUE. ) THEN
792 : : ! rassembler les resultats
793 : : #ifdef HAVE_MPI
794 : : #if defined HAVE_MPI2_INPLACE
795 : : CALL MPI_ALLGATHERV(MPI_IN_PLACE, 0, MPI_DOUBLE_PRECISION, &
796 : : this%oper, counts, displs, &
797 : 518 : MPI_DOUBLE_PRECISION, this%MY_COMM, residu)
798 : : #else
799 : : my_count=tauBegin-tauEnd+1
800 : : MALLOC(oper_buf,(my_count))
801 : : oper_buf(1:my_count)=this%oper(tauBegin:tauEnd)
802 : : CALL MPI_ALLGATHERV(oper_buf, my_count, MPI_DOUBLE_PRECISION, &
803 : : this%oper, counts, displs, &
804 : : MPI_DOUBLE_PRECISION, this%MY_COMM, residu)
805 : : FREE(oper_buf)
806 : : #endif
807 : : #endif
808 [ - + ]: 518 : FREE(counts)
809 [ - + ]: 518 : FREE(displs)
810 : : END IF
811 : 518 : this%oper(tauSamples+1) = A - this%oper(1) !G(0+)-G(0-)=G(0+)+G(beta-)=A
812 : 518 : this%setT = .TRUE.
813 [ - + ]: 518 : FREE(Domega)
814 [ - + ]: 518 : FREE(A_omega)
815 : :
816 [ - + - + ]: 518 : END SUBROUTINE GreenHyb_backFourier
817 : : !!***
818 : :
819 : : !!****f* ABINIT/m_GreenHyb/GreenHyb_forFourier
820 : : !! NAME
821 : : !! GreenHyb_forFourier
822 : : !!
823 : : !! FUNCTION
824 : : !! perform forward fourier transform
825 : : !!
826 : : !! COPYRIGHT
827 : : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
828 : : !! This file is distributed under the terms of the
829 : : !! GNU General Public License, see ~abinit/COPYING
830 : : !! or http://www.gnu.org/copyleft/gpl.txt .
831 : : !!
832 : : !! INPUTS
833 : : !! this=Green
834 : : !! Wmax=linear maximum frequency
835 : : !!
836 : : !! OUTPUT
837 : : !! Gomega=Results for omega frequencies
838 : : !! omega=ask frequencies
839 : : !!
840 : : !! SIDE EFFECTS
841 : : !!
842 : : !! NOTES
843 : : !!
844 : : !! SOURCE
845 : :
846 [ + - + - : 368 : SUBROUTINE GreenHyb_forFourier(this, Gomega, omega, Wmax)
- + + - +
- - + ]
847 : : !Arguments ------------------------------------
848 : :
849 : : #ifdef HAVE_MPI1
850 : : include 'mpif.h'
851 : : #endif
852 : : TYPE(GreenHyb) , INTENT(INOUT) :: this
853 : : COMPLEX(KIND=8), DIMENSION(:), OPTIONAL, INTENT(INOUT) :: Gomega ! INOUT for MPI
854 : : COMPLEX(KIND=8), DIMENSION(:), OPTIONAL, INTENT(IN ) :: omega
855 : : INTEGER , OPTIONAL, INTENT(IN ) :: Wmax
856 : : INTEGER :: i
857 : : INTEGER :: j
858 : : INTEGER :: L
859 : : INTEGER :: Lspline
860 : : INTEGER :: Nom
861 : : INTEGER :: omegaBegin
862 : : INTEGER :: omegaEnd
863 : : INTEGER :: deltaw
864 : : INTEGER :: residu
865 : 368 : INTEGER, ALLOCATABLE, DIMENSION(:) :: counts
866 : 368 : INTEGER, ALLOCATABLE, DIMENSION(:) :: displs
867 : : DOUBLE PRECISION :: beta
868 : : DOUBLE PRECISION :: tau
869 : : DOUBLE PRECISION :: delta
870 : : DOUBLE PRECISION :: deltabis
871 : : DOUBLE PRECISION :: inv_delta
872 : : DOUBLE PRECISION :: inv_delta2
873 : : DOUBLE PRECISION :: omdeltabis
874 : : DOUBLE PRECISION :: tmp
875 : : DOUBLE PRECISION :: xpi
876 : 368 : DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: diag
877 : 368 : DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: diagL
878 : 368 : DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: lastR
879 : 368 : DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: lastC
880 : 368 : DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: XM
881 : 368 : DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: X2
882 : : DOUBLE PRECISION :: iw
883 : : COMPLEX(KIND=8) :: iwtau
884 : 368 : COMPLEX(KIND=8), ALLOCATABLE, DIMENSION(:) :: Gwtmp
885 : 368 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: omegatmp
886 : :
887 : : #if defined HAVE_MPI && !defined HAVE_MPI2_INPLACE
888 : : INTEGER :: my_count
889 : : COMPLEX(KIND=8), ALLOCATABLE , DIMENSION(:) :: Gwtmp_buf
890 : : #endif
891 : :
892 [ - + ]: 368 : IF ( this%set .EQV. .FALSE. ) &
893 : 0 : CALL ERROR("GreenHyb_forFourier : Uninitialized GreenHyb structure")
894 [ - + ]: 368 : IF ( this%setT .EQV. .FALSE. ) &
895 : 0 : CALL ERROR("GreenHyb_forFourier : no G(tau)")
896 [ - + ]: 368 : IF ( this%setMk .NE. 2 ) &
897 : 0 : CALL WARNALL("GreenHyb_forFourier : green does not have moments ")
898 : :
899 : 368 : L = this%samples
900 : :
901 : 368 : xpi=acos(-1.d0) !!! XPI=PI
902 : 368 : beta = this%beta
903 : 368 : Nom = this%Wmax
904 [ + - ]: 368 : IF ( PRESENT(Gomega) ) THEN
905 : 368 : Nom = SIZE(Gomega)
906 : : !IF ( this%rank .EQ. 0 ) &
907 : : !write(6,*) "size Gomega", Nom
908 : : END IF
909 [ + - ]: 368 : IF ( PRESENT(omega) ) THEN
910 [ + - - + ]: 368 : IF ( PRESENT(Gomega) .AND. SIZE(omega) .NE. Nom ) THEN
911 : 0 : CALL ERROR("GreenHyb_forFourier : sizes mismatch ")
912 : : !ELSE
913 : : !Nom = SIZE(omega)
914 : : END IF
915 : : END IF
916 [ + - ]: 368 : IF ( .NOT. PRESENT(Gomega) .AND. .NOT. PRESENT(omega) ) THEN
917 [ # # ]: 0 : IF ( PRESENT(Wmax) ) THEN
918 : 0 : Nom=Wmax
919 : : ELSE
920 : 0 : CALL ERROR("GreenHyb_forFourier : Missing argument Wmax")
921 : : END IF
922 : : END IF
923 : :
924 [ - + ]: 368 : IF ( ALLOCATED(this%oper_w) ) THEN
925 [ # # ]: 0 : IF ( SIZE(this%oper_w) .NE. Nom ) THEN
926 : 0 : FREE(this%oper_w)
927 [ # # # # : 0 : MALLOC(this%oper_w,(1:Nom))
# # # # #
# # # ]
928 : : END IF
929 : : ELSE
930 [ - + + - : 1104 : MALLOC(this%oper_w,(1:Nom))
+ - + - -
+ - + -
+ ]
931 : : END IF
932 : :
933 : : !write(6,*) "PRESENT(GOMEGA)", PRESENT(GOMEGA)
934 : : !write(6,*) "PRESENT(OMEGA)", PRESENT(OMEGA)
935 : : !call flush(6)
936 : :
937 : 368 : delta=this%delta_t
938 : 368 : inv_delta = this%inv_dt
939 : 368 : inv_delta2 = inv_delta*inv_delta
940 : :
941 [ + - - + : 1104 : MALLOC(diagL,(L-1))
+ - + - +
- - + -
+ ]
942 [ - + + - : 736 : MALLOC(lastR,(L-1))
- + - + ]
943 [ - + + - : 1104 : MALLOC(diag,(L))
+ - + - -
+ - + ]
944 [ - + + - : 736 : MALLOC(lastC,(L-1))
- + - + ]
945 : :
946 : : !(cf Stoer) for the spline interpolation :
947 : : ! second derivatives XM solution of A*XM=B.
948 : : !A=(2.4.2.11) of Stoer&Bulirsch + 2 limit conditions
949 : : !The LU decomposition of A is known explicitly;
950 : :
951 : 368 : diag (1) = 4.d0 ! 1.d0 *4.d0 factor 4 added for conditionning
952 : 368 : diagL(1) = 0.25d0 !1.d0/4.d0
953 : 368 : lastR(1) = -0.5d0 ! -2.d0/4.d0
954 : 368 : lastC(1) = 4.d0 ! 1.d0*4.d0
955 : 368 : diag (2) = 4.d0
956 : 368 : diagL(2) = 0.25d0
957 : 368 : lastR(2) = -0.25d0
958 : 368 : lastC(2) = -1.d0
959 : :
960 [ + + ]: 76164 : DO i = 3, L-2
961 : 75796 : tmp = 4.d0 - diagL(i-1)
962 : 76164 : diagL(i) = 1.d0 / tmp
963 : : END DO
964 [ + + ]: 76164 : DO i = 3, L-2
965 : 75796 : diag (i) = 1.d0 / diagL(i)
966 : 75796 : lastR(i) = -(lastR(i-1)*diagL(i))
967 : 76164 : lastC(i) = -(lastC(i-1)*diagL(i-1))
968 : : END DO
969 : :
970 : 368 : tmp = 1.d0/diag(L-2)
971 : 368 : diag (L-1) = 4.d0 - tmp
972 : 368 : lastR(L-1) = (1.d0 - lastR(L-2))/ diag(L-1)
973 : : !diagL(L-1) = lastR(L-1)
974 : 368 : diagL(L-1) = 0.d0 ! for the Lq=B resolution
975 : : !lastC(L-1) = 1.d0 - lastC(L-2)*diagL(L-1) ! equivalent to the next line
976 : 368 : lastC(L-1) = 1.d0 - (lastC(L-2)*lastR(L-1)) ! True value
977 : 368 : diag (L ) = 2.d0! - DOT_PRODUCT( lastR , lastC )
978 : 368 : tmp = 0.d0
979 [ + + ]: 77268 : DO i = 1, L-1
980 : 77268 : tmp = tmp + lastR(i)*lastC(i)
981 : : END DO
982 : 368 : diag (L ) = diag (L ) - tmp
983 : 368 : lastC(L-1) = lastC(L-1)-1.d0 ! 1 is removed for the u.XM=q resolution
984 : :
985 : : ! construct the B vector from A.Xm=B
986 [ - + + - : 736 : MALLOC(XM,(L))
- + - + ]
987 : 368 : XM(1) = 4.d0*this%Mk(3)
988 : : XM(L) = (6.d0 * inv_delta) * ( this%Mk(2) - ( &
989 : : (this%oper(2)-this%oper(1)) + &
990 : 368 : (this%oper(L)-this%oper(L-1)) ) * inv_delta )
991 [ + + ]: 76900 : DO i = 2, L-1
992 : : XM(i) = (6.d0 * inv_delta2) * ( (this%oper(i+1) &
993 : : - 2.d0 * this%oper(i)) &
994 : 76900 : + this%oper(i-1) )
995 : : END DO
996 : :
997 : : ! Find second derivatives XM: Solve the system
998 : : ! SOLVING Lq= XM
999 : : ! q = XM
1000 [ + + ]: 77268 : do j=1,L-1
1001 : 76900 : XM(j+1)=XM(j+1)-(diagL(j)*XM(j))
1002 : 77268 : XM(L) =XM(L) -(lastR(j)*XM(j))
1003 : : end do
1004 : 368 : FREE(diagL)
1005 : 368 : FREE(lastR)
1006 : :
1007 : : ! SOLVING U.XM=q
1008 : : ! XM = q
1009 [ + + ]: 76900 : do j=L-1,2,-1
1010 : 76532 : XM(j+1) = XM(j+1) / diag(j+1)
1011 : 76900 : XM(j)= (XM(j)-(XM(L)*lastC(j)))-XM(j+1)
1012 : : end do
1013 : 368 : XM(2) = XM(2) / diag(2)
1014 : 368 : XM(1) = (XM(1)-XM(L)*lastC(1)) / diag(1)
1015 : :
1016 : 368 : FREE(diag)
1017 : 368 : FREE(lastC)
1018 : :
1019 : 368 : Lspline = L-1
1020 [ - + + - : 1104 : MALLOC(X2,(1:Lspline+1)) ! We impose L = Nom
+ - + - -
+ - + ]
1021 : : !Construct L2 second derivative from known derivatives XM
1022 : 368 : deltabis = beta / DBLE(Lspline)
1023 [ + + ]: 77268 : DO i = 1, Lspline
1024 : 76900 : tau = deltabis * DBLE(i-1)
1025 : 76900 : j = ((L-1)*(i-1))/Lspline + 1!INT(tau * inv_delta) + 1
1026 : 77268 : X2(i) = inv_delta * ( XM(j)*(DBLE(j)*delta - tau ) + XM(j+1)*(tau - DBLE(j-1)*delta) )
1027 : : END DO
1028 : 368 : X2(Lspline+1) = XM(L)
1029 : 368 : FREE(XM)
1030 : :
1031 [ + - ]: 368 : IF ( this%have_MPI .EQV. .TRUE. ) THEN
1032 : 368 : deltaw = Nom / this%size
1033 : 368 : residu = Nom - this%size*deltaw
1034 [ + + ]: 368 : IF ( this%rank .LT. this%size - residu ) THEN
1035 : 282 : omegaBegin = 1 + this%rank*deltaw
1036 : 282 : omegaEnd = (this%rank + 1)*deltaw
1037 : : ELSE
1038 : : ! tauBegin = (this%size-residu)*deltaw + 1 + (this%rank-this%size+residu)*(deltaw+1)
1039 : 86 : omegaBegin = 1 + this%rank*(deltaw + 1) -this%size + residu
1040 : 86 : omegaEnd = omegaBegin + deltaw
1041 : : END IF
1042 [ - + + - : 1104 : MALLOC(counts,(1:this%size))
+ - + - -
+ - + ]
1043 [ - + + - : 736 : MALLOC(displs,(1:this%size))
- + - + ]
1044 : : counts = (/ (deltaw, i=1, this%size-residu), &
1045 [ - + + + : 3928 : (deltaw+1, i=this%size-residu+1, this%size) /)
- + + + +
- + + + -
- + - - +
+ ]
1046 : 368 : displs(1)=0
1047 [ + + ]: 1064 : DO i = 2, this%size
1048 : 1064 : displs(i) = displs(i-1) + counts (i-1)
1049 : : END DO
1050 : : ELSE
1051 : : omegaBegin = 1
1052 : : omegaEnd = Nom
1053 : : END IF
1054 : :
1055 : 368 : this%Mk(1) = -1.d0
1056 [ - + + - : 1104 : MALLOC(omegatmp,(omegaBegin:omegaEnd))
+ - + - -
+ - + ]
1057 [ + - ]: 368 : IF ( PRESENT(omega) ) THEN
1058 [ + - + + : 126544 : omegatmp(omegaBegin:omegaEnd) = (/ (AIMAG(omega(i)),i=omegaBegin,omegaEnd) /)
+ + ]
1059 : : ELSE
1060 [ # # # # : 0 : omegatmp(omegaBegin:omegaEnd) = (/ ((((2.d0*DBLE(i)-1.d0)*xpi)/Beta), i=omegaBegin,omegaEnd) /)
# # ]
1061 : : END IF
1062 [ - + + - : 1104 : MALLOC(Gwtmp,(1:Nom))
+ - + - -
+ - + ]
1063 [ + + ]: 63456 : DO i = omegaBegin, omegaEnd
1064 : 63088 : iw = omegatmp(i)
1065 : 63088 : omdeltabis = iw*deltabis
1066 : 63088 : Gwtmp(i)=CMPLX(0.d0,0.d0,8)
1067 [ + + ]: 25420400 : DO j=2, Lspline ! We impose L+1 = Nom
1068 : 25357312 : iwtau = CMPLX(0.d0,omdeltabis*DBLE(j-1),8)
1069 : 25420400 : Gwtmp(i) = Gwtmp(i) + EXP(iwtau) * CMPLX((X2(j+1) + X2(j-1))-2.d0*X2(j),0.d0,8)
1070 : : END DO
1071 : : Gwtmp(i) = Gwtmp(i)/CMPLX(((iw*iw)*(iw*iw)*deltabis),0.d0,8) &
1072 : :
1073 : : + CMPLX( &
1074 : : ( ((X2(2)-X2(1))+(X2(Lspline+1)-X2(Lspline)))/((iw*iw)*deltabis) -this%Mk(2) ) &
1075 : : /(iw*iw) &
1076 : : , &
1077 : : (this%Mk(1)-this%Mk(3)/(iw*iw))/iw &
1078 : 63456 : , 8)
1079 : : !+ CMPLX( (X2(2)-X2(1))+(X2(Lspline+1)-X2(Lspline)), 0.d0, 8 ) ) &
1080 : : ! / (((iw*iw)*(iw*iw))*CMPLX(deltabis,0.d0,8)) &
1081 : : !- CMPLX(this%Mk(1),0.d0,8)/iw &
1082 : : !+ CMPLX(this%Mk(2),0.d0,8)/(iw*iw) &
1083 : : !- CMPLX(this%Mk(3),0.d0,8)/((iw*iw)*iw)
1084 : : !IF ( this%rank .EQ. 0 ) write(12819,*) iw,gwtmp(i)
1085 : : END DO
1086 [ - + ]: 368 : FREE(omegatmp)
1087 : : !call flush(12819)
1088 : 368 : FREE(X2)
1089 [ + - ]: 368 : IF ( this%have_MPI .EQV. .TRUE. ) THEN
1090 : : #ifdef HAVE_MPI
1091 : : #if defined HAVE_MPI2_INPLACE
1092 : : CALL MPI_ALLGATHERV(MPI_IN_PLACE, 0, MPI_DOUBLE_COMPLEX, &
1093 : : Gwtmp , counts, displs, &
1094 : 368 : MPI_DOUBLE_COMPLEX, this%MY_COMM, residu)
1095 : : #else
1096 : : my_count=omegaBegin-omegaEnd+1
1097 : : MALLOC(Gwtmp_buf,(my_count))
1098 : : Gwtmp_buf(1:my_count)=Gwtmp(omegaBegin:omegaEnd)
1099 : : CALL MPI_ALLGATHERV(Gwtmp_buf, my_count, MPI_DOUBLE_COMPLEX, &
1100 : : Gwtmp , counts, displs, &
1101 : : MPI_DOUBLE_COMPLEX, this%MY_COMM, residu)
1102 : : FREE(Gwtmp_buf)
1103 : : #endif
1104 : : #endif
1105 [ - + ]: 368 : FREE(counts)
1106 [ - + ]: 368 : FREE(displs)
1107 : : END IF
1108 [ + - ]: 368 : IF ( PRESENT(Gomega) ) THEN
1109 [ + + ]: 182346 : Gomega = Gwtmp
1110 : : END IF
1111 [ + - - + : 182714 : this%oper_w = Gwtmp
- - + + ]
1112 : 368 : this%setW = .TRUE.
1113 [ - + ]: 368 : FREE(Gwtmp)
1114 [ - + - + ]: 368 : END SUBROUTINE GreenHyb_forFourier
1115 : : !!***
1116 : :
1117 : : !!****f* ABINIT/m_GreenHyb/GreenHyb_print
1118 : : !! NAME
1119 : : !! GreenHyb_print
1120 : : !!
1121 : : !! FUNCTION
1122 : : !! print Green function
1123 : : !!
1124 : : !! COPYRIGHT
1125 : : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
1126 : : !! This file is distributed under the terms of the
1127 : : !! GNU General Public License, see ~abinit/COPYING
1128 : : !! or http://www.gnu.org/copyleft/gpl.txt .
1129 : : !!
1130 : : !! INPUTS
1131 : : !! this=Green
1132 : : !! ostream=file stream
1133 : : !!
1134 : : !! OUTPUT
1135 : : !!
1136 : : !! SIDE EFFECTS
1137 : : !!
1138 : : !! NOTES
1139 : : !!
1140 : : !! SOURCE
1141 : :
1142 : 0 : SUBROUTINE GreenHyb_print(this, ostream)
1143 : :
1144 : : !Arguments ------------------------------------
1145 : : TYPE(GreenHyb), INTENT(IN) :: this
1146 : : INTEGER, OPTIONAL , INTENT(IN) :: ostream
1147 : : !Local variables ------------------------------
1148 : : INTEGER :: ostream_val
1149 : : INTEGER :: i
1150 : : INTEGER :: samples
1151 : :
1152 : :
1153 [ # # ]: 0 : IF ( this%set .EQV. .FALSE. ) &
1154 : 0 : CALL ERROR("GreenHyb_print : green this%operator not set ")
1155 : :
1156 [ # # ]: 0 : IF ( PRESENT(ostream) ) THEN
1157 : 0 : ostream_val = ostream
1158 : : ELSE
1159 : 0 : ostream_val = 66
1160 : 0 : OPEN(UNIT=ostream_val,FILE="Green.dat")
1161 : : END IF
1162 : :
1163 : 0 : samples = this%samples
1164 : :
1165 [ # # ]: 0 : DO i = 1, samples
1166 : 0 : WRITE(ostream_val,*) DBLE(i-1)*this%delta_t, this%oper(i)
1167 : : END DO
1168 : :
1169 [ # # ]: 0 : IF ( .NOT. PRESENT(ostream) ) &
1170 : 0 : CLOSE(ostream_val)
1171 : 0 : END SUBROUTINE GreenHyb_print
1172 : : !!***
1173 : :
1174 : : !!****f* ABINIT/m_GreenHyb/GreenHyb_destroy
1175 : : !! NAME
1176 : : !! GreenHyb_destroy
1177 : : !!
1178 : : !! FUNCTION
1179 : : !! destroy green function
1180 : : !!
1181 : : !! COPYRIGHT
1182 : : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
1183 : : !! This file is distributed under the terms of the
1184 : : !! GNU General Public License, see ~abinit/COPYING
1185 : : !! or http://www.gnu.org/copyleft/gpl.txt .
1186 : : !!
1187 : : !! INPUTS
1188 : : !! this=Green
1189 : : !!
1190 : : !! OUTPUT
1191 : : !!
1192 : : !! SIDE EFFECTS
1193 : : !!
1194 : : !! NOTES
1195 : : !!
1196 : : !! SOURCE
1197 : :
1198 : 1125 : SUBROUTINE GreenHyb_destroy(this)
1199 : :
1200 : : !Arguments ------------------------------------
1201 : : TYPE(GreenHyb), INTENT(INOUT) :: this
1202 : :
1203 : 1125 : this%set = .FALSE.
1204 : 1125 : this%setT = .FALSE.
1205 : 1125 : this%setW = .FALSE.
1206 : 1125 : this%samples = 0
1207 : 1125 : this%measurements = 0
1208 : 1125 : this%beta = 0.d0
1209 : 1125 : this%inv_beta = 0.d0
1210 : 1125 : this%inv_dt = 0.d0
1211 : 1125 : this%delta_t = 0.d0
1212 : : !CALL VectorInt_destroy(this%index_old)
1213 : : !CALL Vector_destroy(this%oper_old)
1214 : 1125 : CALL MapHyb_destroy(this%this)
1215 [ + - ]: 1125 : FREEIF(this%oper)
1216 [ + + ]: 1125 : FREEIF(this%oper_w)
1217 [ - + ]: 1125 : FREEIF(this%oper_w_old)
1218 [ - + ]: 1125 : FREEIF(this%omega)
1219 : 1125 : END SUBROUTINE GreenHyb_destroy
1220 : : !!***
1221 : :
1222 [ # # # # : 0 : END MODULE m_GreenHyb
# # # # #
# # # # #
# # # # #
# # # #
# ]
1223 : : !!***
|
