Branch data Line data Source code
1 : :
2 : : #if defined HAVE_CONFIG_H
3 : : #include "config.h"
4 : : #endif
5 : : !!****m* ABINIT/m_GreenHyboffdiagComplex
6 : : !! NAME
7 : : !! m_GreenHyboffdiagComplex
8 : : !!
9 : : !! FUNCTION
10 : : !! Manage a green function for one orbital
11 : : !!
12 : : !! COPYRIGHT
13 : : !! Copyright (C) 2013-2025 ABINIT group (B.Amadon, J. Denier and 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_GreenHyboffdiagComplex
24 : :
25 : : USE m_global
26 : : USE m_MatrixHybComplex
27 : : USE m_VectorComplex
28 : : USE m_VectorInt
29 : : USE m_ListCdagC
30 : : USE m_MapHybComplex
31 : : #ifdef HAVE_MPI2
32 : : USE mpi
33 : : #endif
34 : :
35 : : IMPLICIT NONE
36 : :
37 : : public :: GreenHyboffdiagComplex_init
38 : : public :: GreenHyboffdiagComplex_reset
39 : : public :: GreenHyboffdiagComplex_clear
40 : : public :: GreenHyboffdiagComplex_setOperW
41 : : public :: GreenHyboffdiagComplex_measHybrid
42 : : public :: GreenHyboffdiagComplex_getHybrid
43 : : public :: GreenHyboffdiagComplex_setN
44 : : public :: GreenHyboffdiagComplex_setMuD1
45 : : public :: GreenHyboffdiagComplex_setMoments
46 : : public :: GreenHyboffdiagComplex_backFourier
47 : : public :: GreenHyboffdiagComplex_backFourierComplex
48 : : public :: GreenHyboffdiagComplex_forFourier
49 : : public :: GreenHyboffdiagComplex_forFourierComplex
50 : : public :: GreenHyboffdiagComplex_print
51 : : public :: GreenHyboffdiagComplex_destroy
52 : : public :: nfourier3
53 : :
54 : : !!***
55 : :
56 : : !!****t* m_GreenHyboffdiagComplex/GreenHyboffdiagComplex
57 : : !! NAME
58 : : !! GreenHyboffdiagComplex
59 : : !!
60 : : !! FUNCTION
61 : : !! This structured datatype contains the necessary data
62 : : !!
63 : : !! COPYRIGHT
64 : : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
65 : : !! This file is distributed under the terms of the
66 : : !! GNU General Public License, see ~abinit/COPYING
67 : : !! or http://www.gnu.org/copyleft/gpl.txt .
68 : : !!
69 : : !! SOURCE
70 : :
71 : : TYPE GreenHyboffdiagComplex
72 : :
73 : : LOGICAL :: set = .FALSE.
74 : : ! True if variable of type GreenHyboffdiagComplex is initialized
75 : :
76 : : LOGICAL :: setT = .FALSE.
77 : : ! True if variable oper contains data
78 : :
79 : : LOGICAL :: setW = .FALSE.
80 : : ! True if variable oper_w contains data
81 : :
82 : : LOGICAL :: have_MPI = .FALSE.
83 : : ! True if MPI is used.
84 : :
85 : : INTEGER :: setMk = 0
86 : : ! setMk=0 is moments for Fourier transform are not computed
87 : :
88 : : INTEGER :: samples
89 : : ! samples=imaginary time slices (dmftqmc_l+1)
90 : :
91 : : INTEGER :: measurements
92 : : ! number of measurements for the Green's function
93 : :
94 : : INTEGER :: factor
95 : : ! if the move is not accepted, the statistic weight has to be
96 : : ! increased for the current configuration.
97 : :
98 : : INTEGER :: MY_COMM
99 : : ! MPI Communicator
100 : :
101 : : INTEGER :: size
102 : : ! size=1
103 : :
104 : : INTEGER :: rank
105 : : ! rank=0
106 : :
107 : : INTEGER :: Wmax
108 : : ! samples-1 if frequency Green's function
109 : :
110 : : INTEGER :: iTech
111 : : ! Precise if Frequency Green's function is computed or not
112 : :
113 : : INTEGER :: nflavors
114 : : ! Number of flavors
115 : :
116 : : DOUBLE PRECISION :: beta
117 : : ! Inverse of temperature
118 : :
119 : : DOUBLE PRECISION :: inv_beta
120 : : ! Temperature
121 : :
122 : : DOUBLE PRECISION :: delta_t
123 : : ! 1/inv_dt
124 : :
125 : : DOUBLE PRECISION :: inv_dt
126 : : ! (samples-1)/beta
127 : : DOUBLE PRECISION :: signvaluemeas
128 : :
129 : : DOUBLE PRECISION :: signvalueold
130 : :
131 : : COMPLEX(KIND=8) :: phasevaluemeas
132 : :
133 : : COMPLEX(KIND=8) :: phasevalueold
134 : :
135 : : COMPLEX(KIND=8), ALLOCATABLE, DIMENSION(:,:,:) :: oper
136 : : ! oper(samples)
137 : :
138 : : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: omega
139 : : ! omega(Wmax)
140 : :
141 : : COMPLEX(KIND=8), ALLOCATABLE, DIMENSION(:,:,:) :: Mk
142 : : ! Moments for FT
143 : :
144 : : COMPLEX(KIND=8) , ALLOCATABLE, DIMENSION(:,:,:) :: oper_w
145 : : ! Frequency Green's function
146 : :
147 : : COMPLEX(KIND=8) , ALLOCATABLE, DIMENSION(:) :: oper_w_old
148 : : ! Old frequency Green's function (not used)
149 : :
150 : : TYPE(VectorComplex) :: oper_old
151 : : ! useless data
152 : :
153 : : TYPE(VectorInt) :: index_old
154 : : ! useless data
155 : :
156 : : TYPE(MapHybComplex), ALLOCATABLE, DIMENSION(:,:) :: map
157 : : ! value of time and Green's functions computed in GreenHyboffdiagComplex_measHybrid
158 : : ! These values are used to fill op%oper in the same routine.
159 : :
160 : : END TYPE GreenHyboffdiagComplex
161 : : !!***
162 : :
163 : : CONTAINS
164 : : !!***
165 : :
166 : : !!****f* ABINIT/m_GreenHyboffdiagComplex/GreenHyboffdiagComplex_init
167 : : !! NAME
168 : : !! GreenHyboffdiagComplex_init
169 : : !!
170 : : !! FUNCTION
171 : : !! Initialize and allocate
172 : : !!
173 : : !! COPYRIGHT
174 : : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
175 : : !! This file is distributed under the terms of the
176 : : !! GNU General Public License, see ~abinit/COPYING
177 : : !! or http://www.gnu.org/copyleft/gpl.txt .
178 : : !!
179 : : !! INPUTS
180 : : !! op=Green
181 : : !! samples=imaginary time slices
182 : : !! beta=inverse temperature
183 : : !! iTech=SHOULD NOT BE USED => BUGGY
184 : : !! MY_COMM=mpi_communicator
185 : : !!
186 : : !! OUTPUT
187 : : !!
188 : : !! SIDE EFFECTS
189 : : !!
190 : : !! NOTES
191 : : !!
192 : : !! SOURCE
193 : :
194 : 0 : SUBROUTINE GreenHyboffdiagComplex_init(op, samples, beta,nflavors,iTech,MY_COMM)
195 : :
196 : :
197 : : #ifdef HAVE_MPI1
198 : : include 'mpif.h'
199 : : #endif
200 : : !Arguments ------------------------------------
201 : : TYPE(GreenHyboffdiagComplex) , INTENT(INOUT) :: op
202 : : INTEGER , INTENT(IN ) :: samples
203 : : DOUBLE PRECISION, INTENT(IN ) :: beta
204 : : INTEGER , INTENT(IN ) :: nflavors
205 : : !INTEGER , INTENT(IN ) :: Wmax
206 : : INTEGER, OPTIONAL, INTENT(IN ) :: iTech
207 : : INTEGER, OPTIONAL, INTENT(IN ) :: MY_COMM
208 : : !Local variables ------------------------------
209 : : INTEGER :: iflavor,iflavorbis,sp1
210 : : DOUBLE PRECISION :: dt
211 : : #ifdef HAVE_MPI
212 : : INTEGER :: ierr
213 : : #endif
214 : :
215 [ # # ]: 0 : IF ( PRESENT(MY_COMM)) THEN
216 : : #ifdef HAVE_MPI
217 : 0 : op%have_MPI = .TRUE.
218 : 0 : op%MY_COMM = MY_COMM
219 : 0 : CALL MPI_Comm_rank(op%MY_COMM, op%rank, ierr)
220 : 0 : CALL MPI_Comm_size(op%MY_COMM, op%size, ierr)
221 : : #else
222 : : CALL WARN("GreenHyboffdiagComplex_init : MPI is not used ")
223 : : op%have_MPI = .FALSE.
224 : : op%MY_COMM = -1
225 : : op%rank = 0
226 : : op%size = 1
227 : : #endif
228 : : ELSE
229 : 0 : op%have_MPI = .FALSE.
230 : 0 : op%MY_COMM = -1
231 : 0 : op%rank = 0
232 : 0 : op%size = 1
233 : : END IF
234 : :
235 : 0 : sp1 = samples + 1
236 : 0 : op%samples = sp1
237 : 0 : op%measurements = 0
238 : 0 : op%nflavors = nflavors
239 : 0 : op%beta = beta
240 : 0 : op%inv_beta = 1.d0 / beta
241 : 0 : op%inv_dt = DBLE(samples) * op%inv_beta
242 : 0 : dt = 1.d0 / op%inv_dt
243 : 0 : op%delta_t = dt
244 : : !op%Wmax = Wmax
245 : 0 : op%Wmax = -1
246 [ # # ]: 0 : FREEIF(op%Mk)
247 [ # # # # : 0 : MALLOC(op%Mk,(nflavors,nflavors,3))
# # # # #
# # # # #
# # # # #
# ]
248 [ # # ]: 0 : FREEIF(op%oper)
249 [ # # # # : 0 : MALLOC(op%oper,(sp1,nflavors,nflavors))
# # # # #
# # # # #
# # # # #
# # # ]
250 : : ! If we want to measure in frequences
251 : : ! let assume we first have "samples" frequences
252 [ # # ]: 0 : IF ( PRESENT(iTech) ) THEN
253 : 0 : op%iTech = iTech
254 : 0 : SELECT CASE (op%iTech)
255 : : CASE (GREENHYB_TAU) ! omega
256 : 0 : op%iTech = GREENHYB_TAU
257 : : CASE (GREENHYB_OMEGA) ! omega
258 : 0 : op%Wmax = samples
259 [ # # ]: 0 : FREEIF(op%oper_w)
260 [ # # # # : 0 : MALLOC(op%oper_w,(1:op%Wmax,nflavors,nflavors))
# # # # #
# # # # #
# # # # #
# # # ]
261 [ # # ]: 0 : FREEIF(op%oper_w_old)
262 [ # # # # : 0 : MALLOC(op%oper_w_old,(1:op%Wmax))
# # # # #
# # # #
# ]
263 [ # # # # : 0 : op%oper_w = CMPLX(0.d0,0.d0,8)
# # ]
264 [ # # ]: 0 : op%oper_w_old = CMPLX(0.d0,0.d0,8)
265 [ # # ]: 0 : FREEIF(op%omega)
266 [ # # # # : 0 : MALLOC(op%omega,(1:op%Wmax))
# # # # #
# # # #
# ]
267 [ # # # # : 0 : op%omega = (/ ((2.d0 * DBLE(sp1) - 1.d0)*ACOS(-1.d0)*op%inv_beta, sp1=1, op%Wmax) /)
# # # # #
# # # # #
# # # #
# ]
268 : : END SELECT
269 : : ELSE
270 : 0 : op%iTech = GREENHYB_TAU
271 : : END IF
272 : : ! end if
273 : 0 : CALL VectorComplex_init(op%oper_old,10000)
274 : 0 : CALL VectorInt_init(op%index_old,10000)
275 [ # # # # : 0 : DT_FREEIF(op%map)
# # # # #
# ]
276 [ # # # # : 0 : MALLOC(op%map,(nflavors,nflavors))
# # # # #
# # # # #
# # # # #
# # # #
# ]
277 [ # # ]: 0 : do iflavor=1,nflavors
278 [ # # ]: 0 : do iflavorbis=1,nflavors
279 : 0 : CALL MapHybComplex_init(op%map(iflavor,iflavorbis),10000)
280 : : enddo
281 : : enddo
282 : :
283 [ # # # # : 0 : op%oper = cmplx(0.d0,0.d0,kind=8)
# # ]
284 : 0 : op%signvaluemeas = 0.d0
285 : 0 : op%signvalueold = 0.d0
286 : 0 : op%phasevaluemeas = cmplx(0.d0,0.d0,kind=8)
287 : 0 : op%phasevalueold = cmplx(0.d0,0.d0,kind=8)
288 : 0 : op%set = .TRUE.
289 : 0 : op%factor = 1
290 : 0 : op%setMk = 0
291 [ # # # # : 0 : op%Mk = cmplx(0.d0,0.d0,kind=8)
# # ]
292 : 0 : END SUBROUTINE GreenHyboffdiagComplex_init
293 : : !!***
294 : :
295 : : !!****f* ABINIT/m_GreenHyboffdiagComplex/GreenHyboffdiagComplex_reset
296 : : !! NAME
297 : : !! GreenHyboffdiagComplex_reset
298 : : !!
299 : : !! FUNCTION
300 : : !! reset green function
301 : : !!
302 : : !! COPYRIGHT
303 : : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
304 : : !! This file is distributed under the terms of the
305 : : !! GNU General Public License, see ~abinit/COPYING
306 : : !! or http://www.gnu.org/copyleft/gpl.txt .
307 : : !!
308 : : !! INPUTS
309 : : !! op=Green
310 : : !!
311 : : !! OUTPUT
312 : : !!
313 : : !! SIDE EFFECTS
314 : : !!
315 : : !! NOTES
316 : : !!
317 : : !! SOURCE
318 : :
319 : 0 : SUBROUTINE GreenHyboffdiagComplex_reset(op)
320 : :
321 : : !Arguments ------------------------------------
322 : : TYPE(GreenHyboffdiagComplex) , INTENT(INOUT) :: op
323 : :
324 : 0 : CALL GreenHyboffdiagComplex_clear(op)
325 : 0 : op%setMk = 0
326 [ # # # # : 0 : op%Mk = cmplx(0.d0,0.d0,kind=8)
# # ]
327 : 0 : op%setT = .FALSE.
328 : 0 : op%setW = .FALSE.
329 : 0 : END SUBROUTINE GreenHyboffdiagComplex_reset
330 : : !!***
331 : :
332 : :
333 : : !!****f* ABINIT/m_GreenHyboffdiagComplex/GreenHyboffdiagComplex_clear
334 : : !! NAME
335 : : !! GreenHyboffdiagComplex_clear
336 : : !!
337 : : !! FUNCTION
338 : : !! clear green function
339 : : !!
340 : : !! COPYRIGHT
341 : : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
342 : : !! This file is distributed under the terms of the
343 : : !! GNU General Public License, see ~abinit/COPYING
344 : : !! or http://www.gnu.org/copyleft/gpl.txt .
345 : : !!
346 : : !! INPUTS
347 : : !! op=Green
348 : : !!
349 : : !! OUTPUT
350 : : !!
351 : : !! SIDE EFFECTS
352 : : !!
353 : : !! NOTES
354 : : !!
355 : : !! SOURCE
356 : :
357 : 0 : SUBROUTINE GreenHyboffdiagComplex_clear(op)
358 : :
359 : : !Arguments ------------------------------------
360 : : TYPE(GreenHyboffdiagComplex) , INTENT(INOUT) :: op
361 : : INTEGER :: iflavor,iflavorbis
362 : :
363 : : !CALL VectorComplex_clear(op%oper_old)
364 : : !CALL VectorInt_clear(op%index_old)
365 [ # # ]: 0 : do iflavor=1,op%nflavors
366 [ # # ]: 0 : do iflavorbis=1,op%nflavors
367 : 0 : CALL MapHybComplex_clear(op%map(iflavor,iflavorbis))
368 : : enddo
369 : : enddo
370 : 0 : op%measurements = 0
371 [ # # ]: 0 : IF ( ALLOCATED(op%oper) ) &
372 [ # # # # : 0 : op%oper = cmplx(0.d0,0.d0,kind=8)
# # ]
373 : 0 : op%signvaluemeas = 0.d0
374 : 0 : op%signvalueold = 1.d0
375 : 0 : op%phasevaluemeas = cmplx(0.d0,0.d0,kind=8)
376 : 0 : op%phasevalueold = cmplx(1.d0,0.d0,kind=8)
377 [ # # ]: 0 : IF ( op%iTech .EQ. GREENHYB_OMEGA ) THEN
378 [ # # ]: 0 : IF ( ALLOCATED(op%oper_w) ) &
379 [ # # # # : 0 : op%oper_w = CMPLX(0.d0,0.d0,kind=8)
# # ]
380 [ # # ]: 0 : IF ( ALLOCATED(op%oper_w_old) ) &
381 [ # # ]: 0 : op%oper_w_old = CMPLX(0.d0,0.d0,kind=8)
382 : : END IF
383 : 0 : op%factor = 0
384 : 0 : END SUBROUTINE GreenHyboffdiagComplex_clear
385 : : !!***
386 : :
387 : : !!****f* ABINIT/m_GreenHyboffdiagComplex/GreenHyboffdiagComplex_setOperW
388 : : !! NAME
389 : : !! GreenHyboffdiagComplex_setOperW
390 : : !!
391 : : !! FUNCTION
392 : : !! set Green function in frequencies
393 : : !!
394 : : !! COPYRIGHT
395 : : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
396 : : !! This file is distributed under the terms of the
397 : : !! GNU General Public License, see ~abinit/COPYING
398 : : !! or http://www.gnu.org/copyleft/gpl.txt .
399 : : !!
400 : : !! INPUTS
401 : : !! op=Green
402 : : !! Gomega=Input values
403 : : !!
404 : : !! OUTPUT
405 : : !!
406 : : !! SIDE EFFECTS
407 : : !!
408 : : !! NOTES
409 : : !!
410 : : !! SOURCE
411 : :
412 [ # # ]: 0 : SUBROUTINE GreenHyboffdiagComplex_setOperW(op, Gomega)
413 : :
414 : : !Arguments ------------------------------------
415 : : TYPE(GreenHyboffdiagComplex) , INTENT(INOUT) :: op
416 : : COMPLEX(KIND=8), DIMENSION(:,:,:), INTENT(IN ) :: Gomega
417 : : !Loval variables ------------------------------
418 : : INTEGER :: tail
419 : :
420 : 0 : tail = SIZE(Gomega,1)
421 [ # # ]: 0 : IF ( .NOT. op%set ) &
422 : 0 : CALL ERROR("GreenHyboffdiagComplex_setOperW : Uninitialized GreenHyboffdiagComplex structure")
423 [ # # ]: 0 : IF ( ALLOCATED(op%oper_w) ) THEN
424 [ # # # # ]: 0 : IF ( SIZE(op%oper_w) .NE. tail ) THEN
425 : 0 : FREE(op%oper_w)
426 [ # # # # : 0 : MALLOC(op%oper_w,(1:tail,op%nflavors,op%nflavors))
# # # # #
# # # # #
# # # # #
# ]
427 : : END IF
428 : : ELSE
429 [ # # # # : 0 : MALLOC(op%oper_w,(1:tail,op%nflavors,op%nflavors))
# # # # #
# # # # #
# # # # #
# # # ]
430 : : END IF
431 [ # # # # : 0 : op%oper_w(:,:,:) = Gomega(:,:,:)
# # ]
432 : 0 : op%Wmax = tail
433 : 0 : op%setW = .TRUE.
434 : 0 : END SUBROUTINE GreenHyboffdiagComplex_setOperW
435 : : !!***
436 : :
437 : : !!****f* ABINIT/m_GreenHyboffdiagComplex/GreenHyboffdiagComplex_measHybrid
438 : : !! NAME
439 : : !! GreenHyboffdiagComplex_measHybrid
440 : : !!
441 : : !! FUNCTION
442 : : !! Measure Green's function
443 : : !!
444 : : !! COPYRIGHT
445 : : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
446 : : !! This file is distributed under the terms of the
447 : : !! GNU General Public License, see ~abinit/COPYING
448 : : !! or http://www.gnu.org/copyleft/gpl.txt .
449 : : !!
450 : : !! INPUTS
451 : : !! op=Green
452 : : !! Mmatrix=M matrix for the current flavor
453 : : !! ListCdagC_1=list of all creator and annhilator operators
454 : : !! updated=should we accumulate or not
455 : : !!
456 : : !! OUTPUT
457 : : !!
458 : : !! SIDE EFFECTS
459 : : !!
460 : : !! NOTES
461 : : !!
462 : : !! SOURCE
463 : :
464 : 0 : SUBROUTINE GreenHyboffdiagComplex_measHybrid(op, Mmatrix, ListCdagC_1, updated,signvalue,phasevalue,activeflavor)
465 : :
466 : : !Arguments ------------------------------------
467 : : TYPE(GreenHyboffdiagComplex) , INTENT(INOUT) :: op
468 : : TYPE(MatrixHybComplex) , INTENT(IN ) :: Mmatrix
469 : : TYPE(ListCdagC) , INTENT(IN ) :: ListCdagC_1(op%nflavors)
470 : : DOUBLE PRECISION , INTENT(IN ) :: signvalue
471 : : COMPLEX(KIND=8) , INTENT(IN) :: phasevalue
472 : : LOGICAL , INTENT(IN ) :: updated
473 : : INTEGER, OPTIONAL , INTENT(IN ) :: activeflavor
474 : : !Local variables ------------------------------
475 : : INTEGER :: iC
476 : : INTEGER :: iCdag
477 : : INTEGER :: tail
478 : : INTEGER :: tailbis
479 : : !INTEGER :: index
480 : : INTEGER :: idx_old
481 : : INTEGER :: old_size
482 : : ! INTEGER :: omegaSamples
483 : : ! INTEGER :: iomega
484 : : INTEGER :: iflavor
485 : : INTEGER :: iflavorbis
486 : : INTEGER :: iC_m,iC_m_add
487 : : INTEGER :: iCdag_m,iCdag_m_add
488 : : INTEGER :: stail !,ii
489 : : ! DOUBLE PRECISION :: pi_invBeta
490 : : DOUBLE PRECISION :: mbeta_two
491 : : DOUBLE PRECISION :: beta
492 : : DOUBLE PRECISION :: beta_tc
493 : : DOUBLE PRECISION :: tcbeta_tc
494 : : DOUBLE PRECISION :: inv_dt
495 : : DOUBLE PRECISION :: tC,tc_phys
496 : : DOUBLE PRECISION :: tCdag
497 : : DOUBLE PRECISION :: time
498 : : DOUBLE PRECISION :: signe,signe2
499 : : DOUBLE PRECISION :: argument
500 : : COMPLEX(KIND=8) :: argumentc
501 : : INTEGER :: iflavorbegin,iflavorend,prtopt
502 : : !DOUBLE PRECISION :: taupi_invbeta
503 : : !COMPLEX(KIND=8) :: cargument
504 : : !COMPLEX(2*8) :: base_exp
505 : : !COMPLEX(2*8) :: increm_exp
506 : : !write(6,*) "measHybrid"
507 : 0 : prtopt=0
508 [ # # ]: 0 : IF ( op%set .EQV. .FALSE. ) &
509 : 0 : CALL ERROR("GreenHyboffdiagComplex_measHybrid : green operator not set ")
510 : 0 : stail=0
511 [ # # ]: 0 : do iflavor=1,op%nflavors
512 : 0 : stail=stail + ListCdagC_1(iflavor)%tail
513 : : enddo
514 : 0 : iflavorbegin = 1
515 : 0 : iflavorend = op%nflavors
516 : :
517 : : if(present(activeflavor)) then
518 : : if(activeflavor.ne.0) then
519 : : !sui!write(6,*) "measHybrid activeflavor",activeflavor
520 : : iflavorbegin = activeflavor
521 : : iflavorend = activeflavor
522 : : endif
523 : : endif
524 : :
525 [ # # ]: 0 : IF ( stail .NE. Mmatrix%tail ) &
526 : 0 : CALL ERROR("GreenHyboffdiagComplex_measHybrid : ListCdagC & M unconsistent ")
527 : :
528 [ # # ]: 0 : IF ( updated .EQV. .TRUE. ) THEN ! NEW change in the configuration
529 : : ! FIXME SHOULD be much more faster
530 : :
531 : :
532 : : ! write(6,*) "LKLLL2b"
533 : 0 : SELECT CASE(op%iTech)
534 : : CASE (GREENHYB_TAU)
535 : 0 : argument = DBLE(op%factor)
536 : : ! At the beginning old_size=0, then it increases
537 : : ! until
538 : : ! for all values of iC, increment green%oper with the value of the
539 : : ! Green's function in listDBLE(iC) obtained from previous iteration
540 : : ! (below)
541 : : ! ===============================================================
542 : : ! An update has been done. So the Green's function will change
543 : : ! It is thus the good moment to store the previous Green's
544 : : ! function with argument, the number of times this Green's
545 : : ! function has been constant
546 : : ! ===============================================================
547 [ # # ]: 0 : DO iflavor=1, op%nflavors
548 [ # # ]: 0 : DO iflavorbis=1, op%nflavors
549 : 0 : old_size = op%map(iflavor,iflavorbis)%tail
550 : : !write(6,*) "size listDBLE",size(op%map(iflavor,iflavorbis)%listDBLE)
551 : : !write(6,*) " measHybrid",old_size,iflavor,iflavorbis
552 [ # # ]: 0 : DO iC = 1, old_size
553 : : !write(6,*) " measHybrid iflavor,iflavorbis,iC listINT ",iflavor,iflavorbis,iC,op%map(iflavor,iflavorbis)%listINT(iC)
554 : : !write(6,*) " measHybrid listDBLE ",iflavor,iflavorbis,iC,op%map(iflavor,iflavorbis)%listDBLE(iC),argument
555 : : ! write(6,*) "measHybrid ic iflavor1 iflavor2 op%map(iflav1,iflav2) oper",ic,iflavor,iflavorbis,op%map(iflavor,iflavorbis)%listDBLE(iC),op%map(iflavor,iflavorbis)%listDBLE(iC) * op%signvalueold * argument
556 [ # # ]: 0 : if (iflavor .eq. iflavorbis) then
557 : : op%oper(op%map(iflavor,iflavorbis)%listINT(iC),iflavor,iflavorbis) = &
558 : : & cmplx(dble(op%oper(op%map(iflavor,iflavorbis)%listINT(iC),iflavor,iflavorbis) &
559 : 0 : & + op%map(iflavor,iflavorbis)%listDBLE(iC) * op%phasevalueold * argument),0.d0,kind=8)
560 : : else
561 : : op%oper(op%map(iflavor,iflavorbis)%listINT(iC),iflavor,iflavorbis) = &
562 : : op%oper(op%map(iflavor,iflavorbis)%listINT(iC),iflavor,iflavorbis) &
563 : 0 : + op%map(iflavor,iflavorbis)%listDBLE(iC) * op%phasevalueold * argument
564 : : endif
565 : : !if(op%map(iflavor,iflavorbis)%listINT(iC)==1.and.iflavor==iflavorbis) then
566 : : ! if(iflavor==iflavorbis) then
567 : : ! !sui!write(6,*) "G(0)", op%map(iflavor,iflavorbis)%listINT(iC),op%map(iflavor,iflavorbis)%listDBLE(iC) * op%signvalueold,op%oper(op%map(iflavor,iflavorbis)%listINT(iC),iflavor,iflavorbis),iflavor
568 : : ! endif
569 : : ! if(iflavor==1.and.iflavorbis==6.and.op%map(iflavor,iflavorbis)%listINT(iC)==1) then
570 : : ! !prt!if(prtopt==1) write(6,*) "G16(0)", op%map(iflavor,iflavorbis)%listDBLE(iC) * op%signvalueold*argument,op%oper(op%map(iflavor,iflavorbis)%listINT(iC),iflavor,iflavorbis),op%signvalueold
571 : : ! endif
572 : : ! if(iflavor==6.and.iflavorbis==1.and.op%map(iflavor,iflavorbis)%listINT(iC)==1) then
573 : : ! !prt!if(prtopt==1) write(6,*) "G61(0)", op%map(iflavor,iflavorbis)%listDBLE(iC) * op%signvalueold*argument,op%oper(op%map(iflavor,iflavorbis)%listINT(iC),iflavor,iflavorbis),op%signvalueold
574 : : ! endif
575 : : END DO
576 : : ! tail**2 is the number of possible t-t'
577 : : ! MapHybComplex_setSize with resize map tail*tail will thus be the new
578 : : ! op%map%tail
579 : : ! update size of map and map%tail
580 : : CALL MapHybComplex_setSize(op%map(iflavor,iflavorbis),&
581 : 0 : & ListCdagC_1(iflavor)%tail*ListCdagC_1(iflavorbis)%tail)
582 : : END DO
583 : : END DO
584 : 0 : op%signvaluemeas = op%signvaluemeas + op%signvalueold * argument
585 : 0 : op%measurements = op%measurements + op%factor
586 : 0 : op%phasevaluemeas = op%phasevaluemeas + op%phasevalueold * argument
587 : :
588 : : !sui!write(6,*) " measurements", op%measurements
589 : : !sui! write(6,*) " signvaluemeas",op%signvaluemeas,op%signvalueold*argument
590 : : !sui! write(6,*) " signvaluemeas/measurements",op%signvaluemeas/op%measurements
591 : :
592 : : ! This is new measurement, thus op%factor should be put to one
593 : 0 : op%factor = 1
594 : : ! write(6,*) "LKLLL2C"
595 : :
596 : :
597 : : ! initialized index idx_old for the doubles loops over flavors and segments.
598 : :
599 : : ! setup usefull quantities
600 : 0 : beta = op%beta
601 : 0 : mbeta_two = -(beta*0.5d0)
602 : 0 : inv_dt = op%inv_dt
603 : :
604 : : ! WARNING time is not the time but just a temporary variable.
605 : : ! Index Time has been calculated previously and is in mat_tau
606 : :
607 : : ! initialized index for each annihilation time of a segment for a given flavor
608 : :
609 : : ! initialized index for each creation time of a segment for another flavor
610 : :
611 : 0 : iC_m=0
612 : 0 : iC_m_add=0
613 [ # # ]: 0 : DO iflavor=1,op%nflavors
614 : 0 : tail=ListCdagC_1(iflavor)%tail
615 : : !write(6,*) " measHybrid iflavor",iflavor,tail
616 : :
617 : 0 : iCdag_m=0
618 : 0 : iCdag_m_add=0
619 [ # # ]: 0 : DO iflavorbis=1,op%nflavors
620 : 0 : tailbis=ListCdagC_1(iflavorbis)%tail
621 : : !write(6,*) " measHybrid iflavorbis",iflavorbis,tailbis
622 : 0 : idx_old = 0
623 : :
624 [ # # ]: 0 : DO iC = 1, tail
625 : : ! tC is the annihilation (C_) time for segment iC and flavor iflavor
626 : : !-------------------------------------------------------------------
627 : 0 : tC = ListCdagC_1(iflavor)%list(iC,C_)
628 : :
629 : : !iC_m=iC_m+1 ! for Mmatrix%mat
630 : : ! For each flavor iflavor, iC start at \sum_{iflavor1<iflavor} tail(iflavor1)
631 : : ! It thus explains the presence of iC_m_add (same below for iCdag_m_add)
632 : : ! ---------------------------------------------------------------------------------
633 : 0 : iC_m=iC_m_add+iC
634 : 0 : beta_tc = beta - tC
635 : 0 : tcbeta_tc = tC * beta_tc
636 : :
637 : : !write(6,*) " measHybrid iC_m",iC_m
638 : : !write(6,*) " measHybrid tailbis",tailbis
639 [ # # ]: 0 : DO iCdag = 1, tailbis
640 : : !iCdag_m=iCdag_m+1
641 : 0 : iCdag_m=iCdag_m_add+iCdag
642 : : !write(6,*) " measHybrid iCdag_m",iCdag_m
643 : :
644 : : ! tCdag is the creation time for segment iCdag and flavor iflavorbis
645 : 0 : tCdag = ListCdagC_1(iflavorbis)%list(iCdag,Cdag_)
646 : :
647 : : ! --- time is equivalent to time=(tc-tcdag)*(beta-tc) and is only
648 : : ! --- useful for signe
649 : 0 : time = tcbeta_tc - tCdag*beta_tc
650 : :
651 : : !signe = SIGN(1.d0,time)
652 : : !time = time + (signe-1.d0)*mbeta_two
653 : : !signe = signe * SIGN(1.d0,beta-tC)
654 : : !signe = SIGN(1.d0,time) * SIGN(1.d0,beta-tC)
655 : 0 : tc_phys=tc
656 [ # # ]: 0 : if(tc>beta) tc_phys=tc-beta
657 : 0 : signe2=SIGN(1.d0,tc_phys-tcdag)
658 : :
659 [ # # ]: 0 : if(iflavor==iflavorbis) signe = SIGN(1.d0,time)
660 [ # # ]: 0 : if(iflavor/=iflavorbis) signe = signe2
661 : : ! signe = SIGN(1.d0,tc-tcdag)
662 : : ! --- tc>tcdag and beta>tc signe=1 ! segment in the middle or antisegment at the edge
663 : : ! ! tc-tcdag > 0
664 : : ! --- tc<tcdag and beta<tc signe=1 ! never
665 : : ! --- tc>tcdag and beta<tc signe=-1 ! segment at the edges
666 : : ! ! tc'-tcdag < 0 (with tc'=tc-beta) -> signe < 0
667 : : ! --- tc<tcdag and beta>tc signe=-1 ! antisegment in the middle
668 : : ! ! tc-tcdag < 0 (with tc'=tc-beta) -> signe < 0
669 : : ! 22/09/14:
670 : : ! ListCdagC_1 is the list of segment, so we are dealing
671 : : ! only with segment here. However all combination of Cdag
672 : : ! and C are taken, this it is possible that tc<tcdag
673 : :
674 : : ! 21/10/14: Wagt is important are the true times (between
675 : : ! 0 and beta). If tauC>tauCdag signe=+1
676 : : ! If tauC<tauCdag signe=-1
677 : : ! if(tc<tcdag.and.(iflavor==iflavorbis)) then
678 : : ! write(6,*) ListCdagC_1(iflavorbis)%tail
679 : : ! do ii=1, ListCdagC_1(iflavorbis)%tail
680 : : ! write(6,*) ListCdagC_1(iflavorbis)%list(ii,1), ListCdagC_1(iflavorbis)%list(ii,2)
681 : : ! enddo
682 : : ! write(6,*) "tc<tcdag", tc,tcdag,beta,iflavor,iflavorbis
683 : : ! stop
684 : : ! endif
685 : :
686 : : if(tc-tcdag>beta) then
687 : : ! write(6,*) " tc-tcdag > beta ", tcdag-tc,beta
688 : : endif
689 : : !if(tc>beta) then
690 : : ! write(6,*) " TC>BETA"
691 : : ! write(6,*) " iflavor,iflavorbis",iflavor,iflavorbis
692 : : ! write(6,*) " ic,icdag ",ic,icdag
693 : : ! write(6,*) " signe ",signe
694 : : ! write(6,*) " tc,tcdag ",tc,tcdag
695 : : ! write(6,*) " Mmatrix%mat ",Mmatrix%mat(iCdag_m,iC_m)
696 : : ! write(6,*) " Mmatrix%mat_tau ",Mmatrix%mat_tau(iCdag_m,iC_m)
697 : : !endif
698 : : ! Si iflavor/=iflavorbis, tc-tcdag can be negative..so in
699 : : ! this case, on should add beta to tc-tcdag with the minus
700 : : ! sign. NOT DONE HERE??
701 : :
702 : : ! ----- Compute the Green's function as the value of the matrix M for times iCdag and iC.
703 : 0 : argumentc = signe*Mmatrix%mat(iCdag_m,iC_m)
704 : :
705 : : !index = INT( ( time * inv_dt ) + 1.5d0 )
706 : : !IF (index .NE. Mmatrix%mat_tau(iCdag,iC)) THEN
707 : : ! WRITE(*,*) index, Mmatrix%mat_tau(iCdag,iC)
708 : : !! CALL ERROR("Plantage")
709 : : !END IF
710 : :
711 : 0 : idx_old = idx_old + 1
712 : :
713 : : ! --- define the value of listDBLE as a function of idx_old
714 : 0 : op%map(iflavor,iflavorbis)%listDBLE(idx_old) = argumentc
715 : : !write(6,*) " measHybrid listDBLE2 ",iflavor,iflavorbis,idx_old,argumentc
716 : : !op%map%listINT(idx_old) = index
717 : :
718 : : ! --- define the new corresponding value of listINT(idx_old) from mat_tau (integers)
719 : : ! --- idx_old has no meaning but listINT(idx_old) has.
720 : 0 : op%map(iflavor,iflavorbis)%listINT(idx_old) = Mmatrix%mat_tau(iCdag_m,iC_m)
721 : : !write(6,*) " measHybrid idx_old listINT ",idx_old,op%map(iflavor,iflavorbis)%listINT(idx_old)
722 : : !write(6,*) " measHybrid iCdag_m, iC_m, mat_tau",iCdag_m,iC_m,Mmatrix%mat_tau(iCdag_m,iC_m)
723 : : ! if(iflavor==1.and.iflavorbis==2.and.op%map(iflavor,iflavorbis)%listINT(idx_old)==1) then
724 : : ! !prt!if(prtopt==1) write(6,*) "---------------------------"
725 : : ! !prt!if(prtopt==1) write(6,*) "GG12(0)", op%map(iflavor,iflavorbis)%listINT(idx_old),op%map(iflavor,iflavorbis)%listDBLE(idx_old),tcdag,tc,signe,signe2
726 : : ! !prt!if(prtopt==1) write(6,*) " ", tc-tcdag,tc_phys-tcdag
727 : : ! do ii=1, tail
728 : : ! !prt!if(prtopt==1) write(6,*) ii, ListCdagC_1(iflavor)%list(ii,1), ListCdagC_1(iflavor)%list(ii,2)
729 : : ! enddo
730 : : ! do ii=1, tailbis
731 : : ! !prt!if(prtopt==1) write(6,*) ii, ListCdagC_1(iflavorbis)%list(ii,1), ListCdagC_1(iflavorbis)%list(ii,2)
732 : : ! enddo
733 : : ! !prt!if(prtopt==1) write(6,*) "---------------------------"
734 : : ! endif
735 : : ! if(iflavor==1.and.iflavorbis==2.and.op%map(iflavor,iflavorbis)%listINT(idx_old)==999) then
736 : : ! !prt!if(prtopt==1) write(66,*) "---------------------------"
737 : : ! !prt!if(prtopt==1) write(66,*) "GG12(0)", op%map(iflavor,iflavorbis)%listINT(idx_old),op%map(iflavor,iflavorbis)%listDBLE(idx_old),tcdag,tc,signe,signe2
738 : : ! !prt!if(prtopt==1) write(66,*) " ", tc-tcdag,tc_phys-tcdag
739 : : ! do ii=1, tail
740 : : ! !prt!if(prtopt==1) write(66,*) ii, ListCdagC_1(iflavor)%list(ii,1), ListCdagC_1(iflavor)%list(ii,2)
741 : : ! enddo
742 : : ! do ii=1, tailbis
743 : : ! !prt!if(prtopt==1) write(66,*) ii, ListCdagC_1(iflavorbis)%list(ii,1), ListCdagC_1(iflavorbis)%list(ii,2)
744 : : ! enddo
745 : : ! !prt!if(prtopt==1) write(66,*) "---------------------------"
746 : : ! endif
747 : : ! !if(iflavor==2.and.iflavorbis==1.and.op%map(iflavor,iflavorbis)%listINT(idx_old)==1) then
748 : : ! !prt!if(prtopt==1) write(6,*) "GG21(0)", op%map(iflavor,iflavorbis)%listINT(idx_old),op%map(iflavor,iflavorbis)%listDBLE(idx_old),tcdag,tc,signe
749 : : !endif
750 : :
751 : :
752 : : END DO
753 : : END DO
754 : : ! do ii=1,tail*tailbis
755 : : ! !write(6,*) " measHybrid ii,op%map(iflavor,iflavorbis)%listINT(ii)", ii,op%map(iflavor,iflavorbis)%listINT(ii)
756 : : ! enddo
757 : 0 : iCdag_m_add=iCdag_m_add+tailbis
758 : : END DO ! iflavorbis
759 : 0 : iC_m_add=iC_m_add+tail
760 : : END DO ! iflavor
761 : 0 : op%signvalueold = signvalue
762 [ # # ]: 0 : op%phasevalueold = phasevalue
763 : : ! write(6,*) "LKLLL2D"
764 : : CASE (GREENHYB_OMEGA)
765 : : ! argument = DBLE(op%factor)
766 : : ! DO iomega = 1, omegaSamples
767 : : ! op%oper_w(iomega) = op%oper_w(iomega) + op%oper_w_old(iomega) * argument
768 : : ! END DO
769 : : ! op%measurements = op%measurements + op%factor
770 : :
771 : : ! op%factor = 1
772 : : ! beta = op%beta
773 : : ! mbeta_two = -(beta*0.5d0)
774 : : ! pi_invBeta = ACOS(-1.d0)/beta
775 : : ! omegaSamples = op%samples-1
776 : : ! DO iC = 1, tail
777 : : ! tC = ListCdagC_1%list(iC,C_)
778 : : ! DO iCdag = 1, tail
779 : : ! tCdag = ListCdagC_1%list(iCdag,Cdag_)
780 : : ! time = tC - tCdag
781 : :
782 : : ! signe = SIGN(1.d0,time)
783 : : ! time = time + (signe-1.d0)*mbeta_two
784 : : ! signe = signe * SIGN(1.d0,beta-tC)
785 : : ! argument = signe*Mmatrix%mat(iCdag,iC)
786 : :
787 : : ! DO iomega = 1, omegaSamples
788 : : ! !op%oper_w_old(iomega) = Mmatrix%mat_tau(iCdag,iC)*CMPLX(0.d0,argument)
789 : : ! op%oper_w_old(iomega) = EXP(CMPLX(0.d0,op%omega(iomega)*time))*CMPLX(0.d0,argument)
790 : : ! END DO
791 : : ! END DO
792 : : ! END DO
793 : : END SELECT
794 : : ELSE
795 : 0 : op%factor = op%factor + 1
796 : : END IF
797 : 0 : END SUBROUTINE GreenHyboffdiagComplex_measHybrid
798 : : !!***
799 : :
800 : : !!****f* ABINIT/m_GreenHyboffdiagComplex/GreenHyboffdiagComplex_getHybrid
801 : : !! NAME
802 : : !! GreenHyboffdiagComplex_getHybrid
803 : : !!
804 : : !! FUNCTION
805 : : !! reduce green function
806 : : !!
807 : : !! COPYRIGHT
808 : : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
809 : : !! This file is distributed under the terms of the
810 : : !! GNU General Public License, see ~abinit/COPYING
811 : : !! or http://www.gnu.org/copyleft/gpl.txt .
812 : : !!
813 : : !! INPUTS
814 : : !! op=Green
815 : : !!
816 : : !! OUTPUT
817 : : !!
818 : : !! SIDE EFFECTS
819 : : !!
820 : : !! NOTES
821 : : !!
822 : : !! SOURCE
823 : :
824 : 0 : SUBROUTINE GreenHyboffdiagComplex_getHybrid(op)
825 : :
826 : : !Arguments ------------------------------------
827 : : TYPE(GreenHyboffdiagComplex), INTENT(INOUT) :: op
828 : :
829 [ # # ]: 0 : IF ( op%set .EQV. .FALSE. ) &
830 : 0 : CALL ERROR("GreenHyboffdiagComplex_getHybrid : green operator not set ")
831 : :
832 : 0 : SELECT CASE(op%iTech)
833 : : CASE (GREENHYB_TAU)
834 [ # # # # : 0 : op%oper = -(op%oper * op%inv_beta) / (DBLE(op%measurements) * op%delta_t)
# # ]
835 : : !write(6,*) "measurements",op%measurements,op%delta_t,op%inv_beta,op%oper(1,1,1)
836 : : !sui!write(6,*) "signevaluemeas meas",op%signvaluemeas,op%measurements
837 : 0 : op%signvaluemeas = op%signvaluemeas / DBLE(op%measurements)
838 : 0 : op%phasevaluemeas = op%phasevaluemeas / DBLE(op%measurements)
839 : : ! print*, "op%oper",op%oper(1,1,1)
840 : : !sui!write(6,*) "signevaluemeas/meas",op%signvaluemeas
841 : : ! print*, "signevaluemeas/meas",op%signvaluemeas
842 : 0 : op%setT = .TRUE.
843 : : CASE (GREENHYB_OMEGA)
844 [ # # # # : 0 : op%oper_w = -(op%oper_w * op%inv_beta) / (DBLE(op%measurements) * op%delta_t)
# # ]
845 : 0 : op%setW = .TRUE.
846 [ # # # ]: 0 : CALL GreenHyboffdiagComplex_backFourier(op)
847 : : END SELECT
848 : :
849 : 0 : END SUBROUTINE GreenHyboffdiagComplex_getHybrid
850 : : !!***
851 : :
852 : : !!****f* ABINIT/m_GreenHyboffdiagComplex/GreenHyboffdiagComplex_setN
853 : : !! NAME
854 : : !! GreenHyboffdiagComplex_setN
855 : : !!
856 : : !! FUNCTION
857 : : !! impose number of electrons for this flavor
858 : : !!
859 : : !! COPYRIGHT
860 : : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
861 : : !! This file is distributed under the terms of the
862 : : !! GNU General Public License, see ~abinit/COPYING
863 : : !! or http://www.gnu.org/copyleft/gpl.txt .
864 : : !!
865 : : !! INPUTS
866 : : !! op=Green
867 : : !! N=number of electrons
868 : : !!
869 : : !! OUTPUT
870 : : !!
871 : : !! SIDE EFFECTS
872 : : !!
873 : : !! NOTES
874 : : !!
875 : : !! SOURCE
876 : :
877 : 0 : SUBROUTINE GreenHyboffdiagComplex_setN(op,N)
878 : :
879 : : !Arguments ------------------------------------
880 : : TYPE(GreenHyboffdiagComplex) , INTENT(INOUT) :: op
881 : : COMPLEX(KIND=8) , INTENT(IN ) :: N(op%nflavors)
882 : : INTEGER :: iflavor,iflavor2
883 : : !COMPLEX(KIND=8) :: tmpoper
884 : :
885 [ # # ]: 0 : IF ( op%set .EQV. .FALSE. ) &
886 : 0 : CALL ERROR("GreenHyboffdiagComplex_setN: green op%operator not set ")
887 [ # # ]: 0 : DO iflavor=1, op%nflavors
888 : : ! write(6,*) "iflavor",-N(iflavor)*op%signvaluemeas ,2*op%oper(op%samples,iflavor,iflavor),(N(iflavor)-1.d0)*op%signvaluemeas
889 : : ! the mulplication by signvaluemeas is necessary because N is
890 : : ! exactly the number of electrons in the flavor iflavor whereas
891 : : ! op%oper is not exact, because it still has to be divided by
892 : : ! signvaluemeas after the MPIREDUCE
893 : 0 : op%oper(1,iflavor,iflavor) = cmplx(dble((N(iflavor) - 1.d0)*op%phasevaluemeas),0.d0,kind=8)
894 : 0 : op%oper(op%samples,iflavor,iflavor) = - N(iflavor)*op%phasevaluemeas
895 : : !op%oper(op%samples,iflavor,iflavor) = 2*op%oper(op%samples,iflavor,iflavor)
896 : : !op%oper(1,iflavor,iflavor) = 2*op%oper(1,iflavor,iflavor)
897 [ # # ]: 0 : DO iflavor2=1, op%nflavors
898 [ # # ]: 0 : if(iflavor/=iflavor2) then
899 : : ! UNEXPLAINED but MANDATORY to have exact results for U=0 nspinor=4 with pawspnorb=0
900 : : ! Correction: The fact 2 is necessary for edge points because the points are at the
901 : : ! edges.
902 : : ! It is of course necessary to fulfill exact results (U=0).
903 : : !tmpoper=(op%oper(op%samples,iflavor,iflavor2)-op%oper(1,iflavor,iflavor2))
904 : : !op%oper(op%samples,iflavor,iflavor2) = tmpoper
905 : : !op%oper(1,iflavor,iflavor2) = -tmpoper
906 : : !op%oper(op%samples,iflavor,iflavor2) = 2*op%oper(op%samples,iflavor,iflavor2)
907 : : !op%oper(1,iflavor,iflavor2) = 2*op%oper(1,iflavor,iflavor2)
908 : 0 : op%oper(op%samples,iflavor,iflavor2) = 2*op%oper(op%samples,iflavor,iflavor2)
909 : 0 : op%oper(1,iflavor,iflavor2) = 2*op%oper(1,iflavor,iflavor2)
910 : : !op%oper(op%samples,iflavor,iflavor2) = op%oper(op%samples,iflavor,iflavor2)+op%oper(op%samples,iflavor2,iflavor)
911 : : !op%oper(1,iflavor,iflavor2) = op%oper(1,iflavor,iflavor2)+op%oper(1,iflavor2,iflavor)
912 : : endif
913 : : ENDDO
914 : : ENDDO
915 : 0 : END SUBROUTINE GreenHyboffdiagComplex_setN
916 : : !!***
917 : :
918 : : !!****f* ABINIT/m_GreenHyboffdiagComplex/GreenHyboffdiagComplex_setMuD1
919 : : !! NAME
920 : : !! GreenHyboffdiagComplex_setMuD1
921 : : !!
922 : : !! FUNCTION
923 : : !! Set first moments for G
924 : : !!
925 : : !! COPYRIGHT
926 : : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
927 : : !! This file is distributed under the terms of the
928 : : !! GNU General Public License, see ~abinit/COPYING
929 : : !! or http://www.gnu.org/copyleft/gpl.txt .
930 : : !!
931 : : !! INPUTS
932 : : !! op=Green
933 : : !! mu=energy level (irrespectige with fermi level)
934 : : !! d1=first moment of hybridization function ("K")
935 : : !!
936 : : !! OUTPUT
937 : : !! argout(sizeout)=description
938 : : !!
939 : : !! SIDE EFFECTS
940 : : !!
941 : : !! NOTES
942 : : !!
943 : : !! SOURCE
944 : :
945 : 0 : SUBROUTINE GreenHyboffdiagComplex_setMuD1(op,iflavor,iflavor2,mu,d1)
946 : :
947 : : !Arguments ------------------------------------
948 : : !Arguments ------------------------------------
949 : : !scalars
950 : : COMPLEX(KIND=8), INTENT(IN ) :: mu
951 : : COMPLEX(KIND=8), INTENT(IN ) :: d1
952 : : INTEGER , INTENT(IN ) :: iflavor
953 : : INTEGER , INTENT(IN ) :: iflavor2
954 : : !type
955 : : TYPE(GreenHyboffdiagComplex) , INTENT(INOUT) :: op
956 : : !Local variables ------------------------------------
957 : : DOUBLE PRECISION :: mu2
958 : : !*********************************************************************
959 : :
960 : : ! ABI_UNUSED((/d1/))
961 : :
962 : 0 : mu2=0
963 : 0 : if(iflavor==iflavor2) mu2=mu
964 : :
965 [ # # ]: 0 : if(iflavor==iflavor2) then
966 : 0 : op%Mk(iflavor,iflavor2,3) = -d1-(mu*mu)
967 : : !op%Mk(iflavor,iflavor2,3) = -(mu*mu)
968 : 0 : op%Mk(iflavor,iflavor2,2) = -mu
969 : : !sui!write(6,*) "setmud1",iflavor,iflavor2, op%Mk(iflavor,iflavor2,2), op%Mk(iflavor,iflavor2,3)
970 : : else
971 : 0 : op%Mk(iflavor,iflavor2,3) = cmplx(0.d0,0.d0,kind=8)
972 : 0 : op%Mk(iflavor,iflavor2,2) = cmplx(0.d0,0.d0,kind=8)
973 : : endif
974 : 0 : op%setMk = op%setMk + 1
975 : : !write(6,*) "mom1",op%Mk(iflavor,iflavor2,:)
976 : 0 : END SUBROUTINE GreenHyboffdiagComplex_setMuD1
977 : : !!***
978 : :
979 : : !!****f* ABINIT/m_GreenHyboffdiagComplex/GreenHyboffdiagComplex_setMoments
980 : : !! NAME
981 : : !! GreenHyboffdiagComplex_setMoments
982 : : !!
983 : : !! FUNCTION
984 : : !! Compute full moments
985 : : !!
986 : : !! COPYRIGHT
987 : : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
988 : : !! This file is distributed under the terms of the
989 : : !! GNU General Public License, see ~abinit/COPYING
990 : : !! or http://www.gnu.org/copyleft/gpl.txt .
991 : : !!
992 : : !! INPUTS
993 : : !! op=Greenb
994 : : !! u1_iflavor1=\sum_{iflavor2} U_{iflavor2,iflavor1} N_iflavor2
995 : : !! (useful for first moment)
996 : : !! u2=\sum_{iflavor1,iflavor2,iflavor3} U_{iflavor1,iflavor2} N_iflavor2
997 : : !!
998 : : !! OUTPUT
999 : : !!
1000 : : !! SIDE EFFECTS
1001 : : !!
1002 : : !! NOTES
1003 : : !!
1004 : : !! CHI
1005 : : !! Will be filled automatically by the parent script
1006 : : !!
1007 : : !! SOURCE
1008 : :
1009 : 0 : SUBROUTINE GreenHyboffdiagComplex_setMoments(op,iflavor1,iflavor1b,u1,u2,u3)
1010 : :
1011 : : !Arguments ------------------------------------
1012 : : TYPE(GreenHyboffdiagComplex) , INTENT(INOUT) :: op
1013 : : COMPLEX(KIND=8), INTENT(IN ) :: u1
1014 : : COMPLEX(KIND=8), INTENT(IN ) :: u2
1015 : : COMPLEX(KIND=8), INTENT(IN ) :: u3
1016 : : INTEGER , INTENT(IN ) :: iflavor1
1017 : : INTEGER , INTENT(IN ) :: iflavor1b
1018 : :
1019 [ # # ]: 0 : if(iflavor1==iflavor1b) then
1020 : 0 : op%Mk(iflavor1,iflavor1b,1) = cmplx(-1.d0,0.d0,kind=8)
1021 : : ! c_a(3)=-d1-mu*mu-2(-mu)(\sum_{b.ne.a} Uab nb)
1022 : 0 : op%Mk(iflavor1,iflavor1b,3) = op%Mk(iflavor1,iflavor1b,3) - cmplx(2.d0,0.d0,kind=8)*(op%Mk(iflavor1,iflavor1b,2)*u1)
1023 : :
1024 : : ! c_a(2)=-mu+\sum_{b.ne.a} Uab n_b
1025 : 0 : op%Mk(iflavor1,iflavor1b,2) = op%Mk(iflavor1,iflavor1b,2) + u1
1026 : : !sui!write(6,*) "setmiments",iflavor1,iflavor1b,u1
1027 : :
1028 : : ! c_a(3)=c_a(3) + \sum Uab^2 nb + \sum Uba Uca <nbnc>
1029 : : ! ie c_a(3)=-d1+mu*mu-2mu*\sumb Uab nb + \sum Uab^2 nb + \sum Uba Uca <nbnc>
1030 : 0 : op%Mk(iflavor1,iflavor1b,3) = op%Mk(iflavor1,iflavor1b,3) - u2
1031 : : else
1032 : 0 : op%Mk(iflavor1,iflavor1b,1) = cmplx(0.d0,0.d0,kind=8)
1033 : 0 : op%Mk(iflavor1,iflavor1b,2) = u3
1034 : 0 : op%Mk(iflavor1,iflavor1b,3) = cmplx(0.d0,0.d0,kind=8)
1035 : : endif
1036 : : !write(6,*) "mom",iflavor1,iflavor1b, op%Mk(iflavor1,iflavor1b,:)
1037 : :
1038 : 0 : op%setMk = op%setMk + 1
1039 : :
1040 : 0 : END SUBROUTINE GreenHyboffdiagComplex_setMoments
1041 : : !!***
1042 : :
1043 : : !!****f* ABINIT/m_GreenHyboffdiagComplex/GreenHyboffdiagComplex_backFourier
1044 : : !! NAME
1045 : : !! GreenHyboffdiagComplex_backFourier
1046 : : !!
1047 : : !! FUNCTION
1048 : : !! perform back fourier transform
1049 : : !!
1050 : : !! COPYRIGHT
1051 : : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
1052 : : !! This file is distributed under the terms of the
1053 : : !! GNU General Public License, see ~abinit/COPYING
1054 : : !! or http://www.gnu.org/copyleft/gpl.txt .
1055 : : !!
1056 : : !! INPUTS
1057 : : !! op=Green
1058 : : !! dvgc=divergence parameter
1059 : : !!
1060 : : !! OUTPUT
1061 : : !!
1062 : : !! SIDE EFFECTS
1063 : : !!
1064 : : !! NOTES
1065 : : !!
1066 : : !! SOURCE
1067 : :
1068 : 0 : SUBROUTINE GreenHyboffdiagComplex_backFourier(op,dvgc,func,hybri_limit,opt_hybri_limit)
1069 : :
1070 : : use m_fstrings, only : int2char4
1071 : :
1072 : : #ifdef HAVE_MPI1
1073 : : include 'mpif.h'
1074 : : #endif
1075 : : !Arguments ------------------------------------
1076 : : TYPE(GreenHyboffdiagComplex) , INTENT(INOUT) :: op
1077 : : DOUBLE PRECISION, OPTIONAL, INTENT(IN ) :: dvgc
1078 : : CHARACTER(len=5) ,OPTIONAL, INTENT(IN) :: func
1079 : : COMPLEX(KIND=8), DIMENSION(op%nflavors,op%nflavors), OPTIONAL, INTENT(IN) :: hybri_limit
1080 : : INTEGER, OPTIONAL, INTENT(IN) :: opt_hybri_limit
1081 : : !Local variables ------------------------------
1082 : : INTEGER :: itau
1083 : : INTEGER :: iomega
1084 : : INTEGER :: omegaSamples
1085 : : INTEGER :: tauSamples
1086 : : INTEGER :: tauBegin
1087 : : INTEGER :: tauEnd
1088 : : INTEGER :: delta
1089 : : INTEGER :: residu
1090 : : INTEGER :: iflavor1
1091 : : INTEGER :: iflavor2,unitnb !,unitnb1
1092 : 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: counts
1093 : 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: displs
1094 : : DOUBLE PRECISION :: A,AA ! Correction factor
1095 : : COMPLEX(KIND=8) :: B,BB ! Correction factor
1096 : : COMPLEX(KIND=8) :: C !,CC ! Correction factor
1097 : : DOUBLE PRECISION :: inv_beta
1098 : : DOUBLE PRECISION :: pi_invBeta
1099 : : DOUBLE PRECISION :: two_invBeta
1100 : : DOUBLE PRECISION :: minusDt
1101 : : DOUBLE PRECISION :: minusOmegaTau
1102 : : DOUBLE PRECISION :: omegaa
1103 : : DOUBLE PRECISION :: minusTau
1104 : : COMPLEX(KIND=8) :: sumTerm
1105 : : DOUBLE PRECISION :: pi
1106 : : DOUBLE PRECISION :: twoPi
1107 : : COMPLEX(KIND=8) :: correction
1108 : 0 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: Domega
1109 : 0 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: A_omega
1110 : 0 : COMPLEX(KIND=8) , ALLOCATABLE, DIMENSION(:) :: C_omega
1111 : 0 : COMPLEX(KIND=8), ALLOCATABLE, DIMENSION(:) :: opertau
1112 : : CHARACTER(len=5) :: funct
1113 : : character(len=4) :: tag_proc
1114 : : character(len=30) :: tmpfil
1115 : :
1116 : : #if defined HAVE_MPI && !defined HAVE_MPI2_INPLACE
1117 : : INTEGER :: my_count
1118 : : COMPLEX(KIND=8), ALLOCATABLE , DIMENSION(:) :: opertau_buf
1119 : : #endif
1120 : :
1121 [ # # ]: 0 : IF ( op%set .EQV. .FALSE. ) &
1122 : 0 : CALL ERROR("GreenHyboffdiagComplex_backFourier : Uninitialized GreenHyboffdiagComplex structure")
1123 [ # # ]: 0 : IF ( op%setW .EQV. .FALSE. ) &
1124 : 0 : CALL ERROR("GreenHyboffdiagComplex_backFourier : no G(iw)")
1125 : :
1126 : 0 : funct="hybri"
1127 [ # # ]: 0 : if(present(func)) funct=func
1128 : : !sui!write(6,*) funct
1129 : 0 : inv_beta = op%inv_beta
1130 : 0 : two_invBeta = 2.d0 * inv_beta
1131 : 0 : minusDt = - op%delta_t
1132 : 0 : omegaSamples = op%Wmax
1133 : 0 : tauSamples = op%samples-1
1134 : 0 : pi = ACOS(-1.d0)
1135 : 0 : twoPi = 2.d0 * pi
1136 : 0 : pi_invBeta = pi * inv_beta
1137 : : !sui!write(6,*) "omegaSamples",omegaSamples
1138 [ # # # # : 0 : MALLOC(Domega,(1:omegaSamples))
# # # # #
# # # ]
1139 [ # # # # : 0 : MALLOC(A_omega,(1:omegaSamples))
# # # # ]
1140 [ # # # # : 0 : MALLOC(C_omega,(1:omegaSamples))
# # # # #
# # # ]
1141 : : IF ( op%rank .EQ. 0 ) THEN
1142 : : !DO iflavor1 = 1, op%nflavors
1143 : : ! DO iflavor2 = 1, op%nflavors
1144 : : ! write(22236,*) "#",iflavor1,iflavor2
1145 : : ! do iomega=1,op%Wmax
1146 : : ! write(22236,*) (2.d0*DBLE(iomega)-1.d0) * pi_invBeta,real(op%oper_w(iomega,iflavor1,iflavor2))
1147 : : ! enddo
1148 : : ! write(22236,*)
1149 : : ! ENDDO
1150 : : !ENDDO
1151 : : ENDIF
1152 : :
1153 [ # # # # : 0 : op%oper = 0.d0
# # ]
1154 : :
1155 [ # # ]: 0 : DO iflavor1 = 1, op%nflavors
1156 [ # # ]: 0 : DO iflavor2 = 1, op%nflavors
1157 : : ! -- compute limit of function G*(i\omega_n)
1158 [ # # ]: 0 : if(funct=="hybri") then
1159 [ # # ]: 0 : IF ( PRESENT(dvgc) ) THEN
1160 : 0 : A = dvgc
1161 : : ELSE
1162 : : A = AIMAG(op%oper_w(omegaSamples,iflavor1,iflavor2))&! A = \lim_\infty Imag(G(iwn))*(wn)
1163 : 0 : *(2.d0*DBLE(omegaSamples)-1.d0) * pi_invBeta
1164 : : AA = AIMAG(op%oper_w(omegaSamples-10,iflavor1,iflavor2))&! A = \lim_\infty Imag(G(iwn))*(wn)
1165 : 0 : *(2.d0*DBLE(omegaSamples-10)-1.d0) * pi_invBeta
1166 : : B = op%oper_w(omegaSamples,iflavor1,iflavor2)&! A = \lim_\infty (G(iwn))*(iwn)
1167 : 0 : *(2.d0*DBLE(omegaSamples)-1.d0) *pi_invBeta*cmplx(0.d0,1.d0,kind=8)
1168 : : BB = op%oper_w(omegaSamples-10,iflavor1,iflavor2)&! A = \lim_\infty (G(iwn))*(iwn)
1169 : 0 : *(2.d0*DBLE(omegaSamples-10)-1.d0) *pi_invBeta*cmplx(0.d0,1.d0,kind=8)
1170 : : !sui!write(6,*) "B=",iflavor1,iflavor2,B,BB
1171 : : END IF
1172 [ # # # # ]: 0 : else if(iflavor1==iflavor2.and.funct=="green") then
1173 : : A = -1.d0
1174 [ # # # # ]: 0 : else if(iflavor1/=iflavor2.and.funct=="green") then
1175 : 0 : A = 0.d0
1176 : : endif ! funct
1177 : :
1178 : 0 : C=cmplx(-A,0.d0,kind=8)
1179 [ # # ]: 0 : if(present(hybri_limit)) then
1180 [ # # ]: 0 : if(present(opt_hybri_limit)) then
1181 [ # # ]: 0 : if(opt_hybri_limit==1) then
1182 : 0 : C = (hybri_limit(iflavor1,iflavor2))
1183 : : !write(6,*) "Hello C= ",C
1184 : : endif
1185 : : endif
1186 : : endif
1187 : :
1188 : : ! -- correction on G(tau=0) is thus
1189 : 0 : correction = -C*cmplx(0.5d0,0.d0,kind=8)
1190 : :
1191 : : ! -- built frequency mesh
1192 [ # # # # : 0 : Domega = (/ ((2.d0 * DBLE(iomega) - 1.d0)*pi_invbeta, iomega=1, omegaSamples) /)
# # # # #
# # # # #
# # ]
1193 : :
1194 : : ! -- built asymptotic function (-C_ij / iw_n)
1195 [ # # # # : 0 : C_omega = C / (Domega*cmplx(0.d0,1.d0,kind=8))
# # # # ]
1196 : :
1197 : : ! -- built time mesh
1198 [ # # ]: 0 : IF (op%have_MPI .EQV. .TRUE.) THEN
1199 : 0 : delta = tauSamples / op%size
1200 : 0 : residu = tauSamples - op%size*delta
1201 [ # # ]: 0 : IF ( op%rank .LT. op%size - residu ) THEN
1202 : 0 : tauBegin = 1 + op%rank*delta
1203 : 0 : tauEnd = (op%rank + 1)*delta
1204 : : ELSE
1205 : : ! tauBegin = (op%size-residu)*delta + 1 + (op%rank-op%size+residu)*(delta+1)
1206 : 0 : tauBegin = 1 + op%rank*(delta + 1) -op%size + residu
1207 : 0 : tauEnd = tauBegin + delta
1208 : : END IF
1209 [ # # # # : 0 : MALLOC(counts,(1:op%size))
# # # # #
# # # #
# ]
1210 [ # # # # : 0 : MALLOC(displs,(1:op%size))
# # # # #
# ]
1211 : : counts = (/ (delta, iTau=1, op%size-residu), &
1212 [ # # # # : 0 : (delta+1, iTau=op%size-residu+1, op%size) /)
# # # # #
# # # # #
# # # # #
# ]
1213 : 0 : displs(1)=0
1214 [ # # ]: 0 : DO iTau = 2, op%size
1215 : 0 : displs(iTau) = displs(iTau-1) + counts (iTau-1)
1216 : : END DO
1217 : : ELSE
1218 : : tauBegin = 1
1219 : : tauEnd = tauSamples
1220 : : END IF
1221 [ # # # # : 0 : MALLOC(opertau,(1:tauSamples+1))
# # # # #
# # # ]
1222 : : !do iomega=1,omegaSamples
1223 : : ! write(6,*) iomega, imag(op%oper_w(iomega,iflavor1,iflavor2)), A_omega(iomega) ,"#diff"
1224 : : !enddo
1225 : 0 : unitnb=70000+op%rank
1226 : 0 : call int2char4(op%rank,tag_proc)
1227 : 0 : tmpfil = 'counts'//tag_proc
1228 : : ! open (unit=unitnb,file=trim(tmpfil),status='unknown',form='formatted')
1229 : : ! write(unitnb,*) "#",iflavor1,iflavor2
1230 : : ! do itau=1,op%size
1231 : : ! write(unitnb,*) itau,counts(itau),displs(itau)
1232 : : ! enddo
1233 : : ! write(unitnb,*)
1234 : :
1235 : 0 : unitnb=10000+op%rank
1236 : 0 : call int2char4(op%rank,tag_proc)
1237 : 0 : tmpfil = 'oper_w'//tag_proc
1238 : : ! open (unit=unitnb,file=trim(tmpfil),status='unknown',form='formatted')
1239 : : ! write(unitnb,*) "#",iflavor1,iflavor2,C
1240 : : ! ! C_omega et oper_w differents Domega identique. Est ce du a des
1241 : : !! ! diago differentes pour chaque procs dans qmc_prep_ctqmc
1242 : : !! do iomega=1,op%Wmax
1243 : : !! write(unitnb,*) (2.d0*DBLE(iomega)-1.d0) * pi_invBeta,real(op%oper_w(iomega,iflavor1,iflavor2)),C_omega(iomega),Domega(iomega)
1244 : : ! enddo
1245 : : ! write(unitnb,*)
1246 : :
1247 : : ! unitnb=40000+op%rank
1248 : : ! unitnb1=50000+op%rank
1249 : : ! call int2char4(op%rank,tag_proc)
1250 : : ! tmpfil = 'tauend'//tag_proc
1251 : : ! open (unit=unitnb,file=trim(tmpfil),status='unknown',form='formatted')
1252 : : ! tmpfil = 'taubegin'//tag_proc
1253 : : ! open (unit=unitnb1,file=trim(tmpfil),status='unknown',form='formatted')
1254 : : ! write(unitnb,*) "#",iflavor1,iflavor2
1255 : : ! write(unitnb1,*) "#",iflavor1,iflavor2
1256 : :
1257 : : ! -- compute Fourier transformation
1258 [ # # ]: 0 : opertau=0.d0
1259 [ # # ]: 0 : DO itau = tauBegin, tauEnd
1260 : : !DO itau = max(tauBegin-1,1), tauEnd
1261 : 0 : minusTau = DBLE(itau -1) * minusDt
1262 [ # # ]: 0 : DO iomega = 1, omegaSamples
1263 : 0 : omegaa = Domega(iomega)
1264 : 0 : minusOmegaTau = MOD(omegaa*minusTau, TwoPi)
1265 : : sumTerm = ( op%oper_w(iomega,iflavor1,iflavor2) &
1266 : : - C_omega(iomega) ) &
1267 : : !- CMPLX(0.d0, A_omega(iomega),8) ) &
1268 : 0 : * EXP( CMPLX(0.d0, minusOmegaTau, 8))
1269 : 0 : opertau(itau) = opertau(itau) + sumTerm
1270 : :
1271 : : END DO
1272 : 0 : opertau(itau) = correction + two_invBeta*opertau(itau)
1273 : : !write(*,*) "itau opertau(itau)",opertau(itau)
1274 : : END DO
1275 : :
1276 : : ! -- Gather
1277 [ # # ]: 0 : IF ( op%have_MPI .EQV. .TRUE. ) THEN
1278 : : ! rassembler les resultats
1279 : : #ifdef HAVE_MPI
1280 : : #if defined HAVE_MPI2_INPLACE
1281 : : CALL MPI_ALLGATHERV(MPI_IN_PLACE, 0, MPI_DOUBLE_COMPLEX, &
1282 : : opertau, counts, displs, &
1283 : 0 : MPI_DOUBLE_COMPLEX, op%MY_COMM, residu)
1284 : : #else
1285 : : my_count=tauBegin-tauEnd+1
1286 : : MALLOC(opertau_buf,(my_count))
1287 : : opertau_buf(1:my_count)=opertau(tauBegin:tauEnd)
1288 : : CALL MPI_ALLGATHERV(opertau_buf, my_count, MPI_DOUBLE_COMPLEX, &
1289 : : opertau, counts, displs, &
1290 : : MPI_DOUBLE_COMPLEX, op%MY_COMM, residu)
1291 : : FREE(opertau_buf)
1292 : : #endif
1293 : : #endif
1294 [ # # ]: 0 : FREE(counts)
1295 [ # # ]: 0 : FREE(displs)
1296 : : END IF
1297 : : ! unitnb=30000+op%rank
1298 : : ! call int2char4(op%rank,tag_proc)
1299 : : ! tmpfil = 'opertau_MPI_'//tag_proc
1300 : : ! open (unit=unitnb,file=trim(tmpfil),status='unknown',form='formatted')
1301 : : ! write(unitnb,*) "#",iflavor1,iflavor2
1302 : : ! do itau=tauBegin, tauEnd
1303 : : ! write(unitnb,*) itau,opertau(itau)
1304 : : ! enddo
1305 : : ! write(unitnb,*)
1306 : :
1307 : : ! -- Add correction for discontinuity.
1308 : : ! if(iflavor1==iflavor2) then
1309 : : !G(0+)-G(0-)=G(0+)+G(beta-)=A
1310 : 0 : opertau(tauSamples+1) = -C - opertau(1)
1311 : : !sui!write(6,*) "BackFourier",opertau(tauSamples+1),opertau(1),real(C)
1312 : :
1313 : 0 : op%setT = .TRUE.
1314 : : ! endif
1315 [ # # ]: 0 : op%oper(:,iflavor1,iflavor2)=opertau(:)
1316 : 0 : FREE(opertau)
1317 : : END DO ! iflavor2
1318 : : END DO ! iflavor1
1319 : : ! -- End loop over flavors.
1320 : :
1321 [ # # ]: 0 : FREE(Domega)
1322 [ # # ]: 0 : FREE(A_omega)
1323 [ # # ]: 0 : FREE(C_omega)
1324 : 0 : close(236)
1325 : 0 : close(237)
1326 : :
1327 [ # # # # : 0 : END SUBROUTINE GreenHyboffdiagComplex_backFourier
# # ]
1328 : : !!***
1329 : :
1330 : : !****f* ABINIT/m_GreenHyboffdiagComplex/GreenHyboffdiagComplex_backFourierComplex
1331 : : !! NAME
1332 : : !! GreenHyboffdiagComplex_backFourierComplex
1333 : : !!
1334 : : !! FUNCTION
1335 : : !! perform back fourier transform for complex green's function
1336 : : !!
1337 : : !! COPYRIGHT
1338 : : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
1339 : : !! This file is distributed under the terms of the
1340 : : !! GNU General Public License, see ~abinit/COPYING
1341 : : !! or http://www.gnu.org/copyleft/gpl.txt .
1342 : : !!
1343 : : !! INPUTS
1344 : : !! op=Green
1345 : : !! dvgc=divergence parameter
1346 : : !!
1347 : : !! OUTPUT
1348 : : !!
1349 : : !! SIDE EFFECTS
1350 : : !!
1351 : : !! NOTES
1352 : : !!
1353 : : !! SOURCE
1354 : :
1355 : 0 : SUBROUTINE GreenHyboffdiagComplex_backFourierComplex(op,dvgc,func,hybri_limit,opt_hybri_limit)
1356 : :
1357 : : use m_fstrings, only : int2char4
1358 : :
1359 : : #ifdef HAVE_MPI1
1360 : : include 'mpif.h'
1361 : : #endif
1362 : : !Arguments ------------------------------------
1363 : : TYPE(GreenHyboffdiagComplex) , INTENT(INOUT) :: op
1364 : : DOUBLE PRECISION, OPTIONAL, INTENT(IN ) :: dvgc
1365 : : CHARACTER(len=5) ,OPTIONAL, INTENT(IN) :: func
1366 : : COMPLEX(KIND=8), DIMENSION(op%nflavors,op%nflavors), OPTIONAL, INTENT(IN) :: hybri_limit
1367 : : INTEGER, OPTIONAL, INTENT(IN) :: opt_hybri_limit
1368 : : !Local variables ------------------------------
1369 : : INTEGER :: itau
1370 : : INTEGER :: iomega
1371 : : INTEGER :: omegaSamples
1372 : : INTEGER :: tauSamples
1373 : : INTEGER :: tauBegin
1374 : : INTEGER :: tauEnd
1375 : : INTEGER :: delta
1376 : : INTEGER :: residu
1377 : : INTEGER :: iflavor1
1378 : : INTEGER :: iflavor2,unitnb !,unitnb1
1379 : 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: counts
1380 : 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: displs
1381 : : DOUBLE PRECISION :: A,AA ! Correction factor
1382 : : COMPLEX(KIND=8) :: B,BB ! Correction factor
1383 : : COMPLEX(KIND=8) :: C !,CC ! Correction factor
1384 : : DOUBLE PRECISION :: inv_beta
1385 : : DOUBLE PRECISION :: pi_invBeta
1386 : : DOUBLE PRECISION :: two_invBeta
1387 : : DOUBLE PRECISION :: minusDt
1388 : : DOUBLE PRECISION :: minusOmegaTau
1389 : : DOUBLE PRECISION :: omegaa
1390 : : DOUBLE PRECISION :: minusTau
1391 : : COMPLEX(KIND=8) :: sumTerm_ab,sumTerm_ba
1392 : : DOUBLE PRECISION :: pi
1393 : : DOUBLE PRECISION :: twoPi
1394 : : COMPLEX(KIND=8) :: correction
1395 : 0 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: Domega
1396 : 0 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: A_omega
1397 : 0 : COMPLEX(KIND=8) , ALLOCATABLE, DIMENSION(:) :: C_omega
1398 : 0 : COMPLEX(KIND=8), ALLOCATABLE, DIMENSION(:) :: opertau
1399 : : CHARACTER(len=5) :: funct
1400 : : character(len=4) :: tag_proc
1401 : : character(len=30) :: tmpfil
1402 : :
1403 : : #if defined HAVE_MPI && !defined HAVE_MPI2_INPLACE
1404 : : INTEGER :: my_count
1405 : : COMPLEX(KIND=8), ALLOCATABLE , DIMENSION(:) :: opertau_buf,opertau_bufab,opertau_bufba
1406 : : #endif
1407 : :
1408 [ # # ]: 0 : IF ( op%set .EQV. .FALSE. ) &
1409 : 0 : CALL ERROR("GreenHyboffdiagComplex_backFourier : Uninitialized GreenHyboffdiagComplex structure")
1410 [ # # ]: 0 : IF ( op%setW .EQV. .FALSE. ) &
1411 : 0 : CALL ERROR("GreenHyboffdiagComplex_backFourier : no G(iw)")
1412 : :
1413 : 0 : funct="hybri"
1414 [ # # ]: 0 : if(present(func)) funct=func
1415 : 0 : inv_beta = op%inv_beta
1416 : 0 : two_invBeta = 2.d0 * inv_beta
1417 : 0 : minusDt = - op%delta_t
1418 : 0 : omegaSamples = op%Wmax
1419 : 0 : tauSamples = op%samples-1
1420 : 0 : pi = ACOS(-1.d0)
1421 : 0 : twoPi = 2.d0 * pi
1422 : 0 : pi_invBeta = pi * inv_beta
1423 [ # # # # : 0 : MALLOC(Domega,(1:omegaSamples))
# # # # #
# # # ]
1424 [ # # # # : 0 : MALLOC(A_omega,(1:omegaSamples))
# # # # ]
1425 [ # # # # : 0 : MALLOC(C_omega,(1:omegaSamples))
# # # # #
# # # ]
1426 : : !IF ( op%rank .EQ. 0 ) THEN
1427 : : !DO iflavor1 = 1, op%nflavors
1428 : : ! DO iflavor2 = 1, op%nflavors
1429 : : ! write(22236,*) "#",iflavor1,iflavor2
1430 : : ! do iomega=1,op%Wmax
1431 : : ! write(22236,*) (2.d0*DBLE(iomega)-1.d0) * pi_invBeta,real(op%oper_w(iomega,iflavor1,iflavor2))
1432 : : ! enddo
1433 : : ! write(22236,*)
1434 : : ! ENDDO
1435 : : !ENDDO
1436 : : !ENDIF
1437 : :
1438 [ # # # # : 0 : op%oper = 0.d0
# # ]
1439 : :
1440 [ # # ]: 0 : DO iflavor1 = 1, op%nflavors
1441 [ # # ]: 0 : DO iflavor2 = 1, op%nflavors
1442 : : ! -- compute limit of function G*(i\omega_n)
1443 [ # # ]: 0 : if(funct=="hybri") then
1444 [ # # ]: 0 : IF ( PRESENT(dvgc) ) THEN
1445 : 0 : A = dvgc
1446 : : ELSE
1447 : : A = AIMAG(op%oper_w(omegaSamples,iflavor1,iflavor2))&! A = \lim_\infty Imag(G(iwn))*(wn)
1448 : 0 : *(2.d0*DBLE(omegaSamples)-1.d0) * pi_invBeta
1449 : : AA = AIMAG(op%oper_w(omegaSamples-10,iflavor1,iflavor2))&! A = \lim_\infty Imag(G(iwn))*(wn)
1450 : 0 : *(2.d0*DBLE(omegaSamples-10)-1.d0) * pi_invBeta
1451 : : B = op%oper_w(omegaSamples,iflavor1,iflavor2)&! A = \lim_\infty (G(iwn))*(iwn)
1452 : 0 : *(2.d0*DBLE(omegaSamples)-1.d0) *pi_invBeta*cmplx(0.d0,1.d0,kind=8)
1453 : : BB = op%oper_w(omegaSamples-10,iflavor1,iflavor2)&! A = \lim_\infty (G(iwn))*(iwn)
1454 : 0 : *(2.d0*DBLE(omegaSamples-10)-1.d0) *pi_invBeta*cmplx(0.d0,1.d0,kind=8)
1455 : : !sui!write(6,*) "B=",iflavor1,iflavor2,B,BB
1456 : : END IF
1457 [ # # # # ]: 0 : else if(iflavor1==iflavor2.and.funct=="green") then
1458 : : A = -1.d0
1459 [ # # # # ]: 0 : else if(iflavor1/=iflavor2.and.funct=="green") then
1460 : 0 : A = 0.d0
1461 : : endif ! funct
1462 : :
1463 : 0 : C=cmplx(-A,0.d0,kind=8)
1464 : : !C=B
1465 : : !write(*,*) "Hello C before hybri_limit", C
1466 [ # # ]: 0 : if(present(hybri_limit)) then
1467 [ # # ]: 0 : if(present(opt_hybri_limit)) then
1468 [ # # ]: 0 : if(opt_hybri_limit==1) then
1469 : 0 : C= (hybri_limit(iflavor1,iflavor2))
1470 : : !write(6,*) "C= ",C
1471 : : endif
1472 : : endif
1473 : : endif
1474 : :
1475 : : ! -- correction on G(tau=0) is thus
1476 : 0 : correction = -C*cmplx(0.5d0,0.d0,kind=8)
1477 : :
1478 : : ! -- built frequency mesh
1479 [ # # # # : 0 : Domega = (/ ((2.d0 * DBLE(iomega) - 1.d0)*pi_invbeta, iomega=1, omegaSamples) /)
# # # # #
# # # # #
# # ]
1480 : :
1481 : : ! -- built asymptotic function (-C_ij / iw_n)
1482 [ # # # # : 0 : C_omega = C / (Domega*cmplx(0.d0,1.d0,kind=8))
# # # # ]
1483 : :
1484 : : ! -- built time mesh
1485 [ # # ]: 0 : IF (op%have_MPI .EQV. .TRUE.) THEN
1486 : 0 : delta = tauSamples / op%size
1487 : 0 : residu = tauSamples - op%size*delta
1488 [ # # ]: 0 : IF ( op%rank .LT. op%size - residu ) THEN
1489 : 0 : tauBegin = 1 + op%rank*delta
1490 : 0 : tauEnd = (op%rank + 1)*delta
1491 : : ELSE
1492 : : ! tauBegin = (op%size-residu)*delta + 1 + (op%rank-op%size+residu)*(delta+1)
1493 : 0 : tauBegin = 1 + op%rank*(delta + 1) -op%size + residu
1494 : 0 : tauEnd = tauBegin + delta
1495 : : END IF
1496 [ # # # # : 0 : MALLOC(counts,(1:op%size))
# # # # #
# # # #
# ]
1497 [ # # # # : 0 : MALLOC(displs,(1:op%size))
# # # # #
# ]
1498 : : counts = (/ (delta, iTau=1, op%size-residu), &
1499 [ # # # # : 0 : (delta+1, iTau=op%size-residu+1, op%size) /)
# # # # #
# # # # #
# # # # #
# ]
1500 : 0 : displs(1)=0
1501 [ # # ]: 0 : DO iTau = 2, op%size
1502 : 0 : displs(iTau) = displs(iTau-1) + counts (iTau-1)
1503 : : END DO
1504 : : ELSE
1505 : : tauBegin = 1
1506 : : tauEnd = tauSamples
1507 : : END IF
1508 [ # # # # : 0 : MALLOC(opertau,(1:tauSamples+1))
# # # # #
# # # ]
1509 : : ! MALLOC(opertau_ab,(1:tauSamples+1))
1510 : : ! MALLOC(opertau_ba,(1:tauSamples+1))
1511 : : !do iomega=1,omegaSamples
1512 : : ! write(6,*) iomega, imag(op%oper_w(iomega,iflavor1,iflavor2)), A_omega(iomega) ,"#diff"
1513 : : !enddo
1514 : 0 : unitnb=70000+op%rank
1515 : 0 : call int2char4(op%rank,tag_proc)
1516 : 0 : tmpfil = 'counts'//tag_proc
1517 : : ! open (unit=unitnb,file=trim(tmpfil),status='unknown',form='formatted')
1518 : : ! write(unitnb,*) "#",iflavor1,iflavor2
1519 : : ! do itau=1,op%size
1520 : : ! write(unitnb,*) itau,counts(itau),displs(itau)
1521 : : ! enddo
1522 : : ! write(unitnb,*)
1523 : :
1524 : 0 : unitnb=10000+op%rank
1525 : 0 : call int2char4(op%rank,tag_proc)
1526 : 0 : tmpfil = 'oper_w'//tag_proc
1527 : :
1528 : :
1529 : : ! -- compute Fourier transformation
1530 [ # # ]: 0 : opertau=cmplx(0.d0,0.d0,kind=8)
1531 : :
1532 [ # # ]: 0 : DO itau = tauBegin, tauEnd
1533 : 0 : minusTau = DBLE(itau -1) * minusDt
1534 [ # # ]: 0 : DO iomega = 1, omegaSamples
1535 : 0 : omegaa = Domega(iomega)
1536 : 0 : minusOmegaTau = MOD(omegaa*minusTau, TwoPi)
1537 : : !== Original code from Jordan ==
1538 : : !sumTerm = REAL((op%oper_w(iomega,iflavor1,iflavor2) - C_omega(iomega) ) * EXP( CMPLX(0.d0, minusOmegaTau, 8))
1539 : : !opertau(itau) = opertau(itau) + sumTerm
1540 : : !!
1541 : : !== Complex version (related to shinaoka) ==
1542 : 0 : sumTerm_ab = (op%oper_w(iomega,iflavor1,iflavor2) - C_omega(iomega) ) * EXP( CMPLX(0.d0, minusOmegaTau, 8))
1543 : 0 : sumTerm_ba = (conjg( op%oper_w(iomega,iflavor2,iflavor1)) + C_omega(iomega) ) * EXP( -1*CMPLX(0.d0, minusOmegaTau, 8))
1544 : 0 : opertau(itau) = opertau(itau) + sumTerm_ab + sumTerm_ba
1545 : : END DO
1546 : : !== Jordan ==
1547 : : !opertau(itau) = correction + two_invbeta*(opertau(itau))
1548 : : !== Complex version ==
1549 : 0 : opertau(itau) = correction + inv_beta*(opertau(itau))
1550 : : END DO
1551 : :
1552 : : ! -- Gather
1553 [ # # ]: 0 : IF ( op%have_MPI .EQV. .TRUE. ) THEN
1554 : : ! rassembler les resultats
1555 : : #ifdef HAVE_MPI
1556 : : #if defined HAVE_MPI2_INPLACE
1557 : : CALL MPI_ALLGATHERV(MPI_IN_PLACE, 0, MPI_DOUBLE_COMPLEX, &
1558 : : opertau, counts, displs, &
1559 : 0 : MPI_DOUBLE_COMPLEX, op%MY_COMM, residu)
1560 : : ! CALL MPI_ALLGATHERV(MPI_IN_PLACE, 0, MPI_DOUBLE_COMPLEX, &
1561 : : ! opertau_ab, counts, displs, &
1562 : : ! MPI_DOUBLE_COMPLEX, op%MY_COMM, residu)
1563 : : ! CALL MPI_ALLGATHERV(MPI_IN_PLACE, 0, MPI_DOUBLE_COMPLEX, &
1564 : : ! opertau_ab, counts, displs, &
1565 : : ! MPI_DOUBLE_COMPLEX, op%MY_COMM, residu)
1566 : : #else
1567 : : my_count=tauBegin-tauEnd+1
1568 : : MALLOC(opertau_buf,(my_count))
1569 : : ! MALLOC(opertau_bufab,(my_count))
1570 : : ! MALLOC(opertau_bufba,(my_count))
1571 : : opertau_buf(1:my_count)=opertau(tauBegin:tauEnd)
1572 : : ! opertau_bufab(1:my_count)=opertau_ab(tauBegin:tauEnd)
1573 : : ! opertau_bufba(1:my_count)=opertau_ba(tauBegin:tauEnd)
1574 : : CALL MPI_ALLGATHERV(opertau_buf, my_count, MPI_DOUBLE_COMPLEX, &
1575 : : opertau, counts, displs, &
1576 : : MPI_DOUBLE_COMPLEX, op%MY_COMM, residu)
1577 : : ! CALL MPI_ALLGATHERV(opertau_bufab, my_count, MPI_DOUBLE_COMPLEX, &
1578 : : ! opertau_ab, counts, displs, &
1579 : : ! MPI_DOUBLE_COMPLEX, op%MY_COMM, residu)
1580 : : ! CALL MPI_ALLGATHERV(opertau_bufba, my_count, MPI_DOUBLE_COMPLEX, &
1581 : : ! opertau_ba, counts, displs, &
1582 : : ! MPI_DOUBLE_COMPLEX, op%MY_COMM, residu)
1583 : :
1584 : : FREE(opertau_buf)
1585 : : ! FREE(opertau_bufab)
1586 : : ! FREE(opertau_bufba)
1587 : : #endif
1588 : : #endif
1589 [ # # ]: 0 : FREE(counts)
1590 [ # # ]: 0 : FREE(displs)
1591 : : END IF
1592 : : ! -- Add correction for discontinuity.
1593 : : ! if(iflavor1==iflavor2) then
1594 : : !G(0+)-G(0-)=G(0+)+G(beta-)=A
1595 : 0 : opertau(tauSamples+1) = -C - opertau(1)
1596 : : ! opertau_ab(tauSamples+1) = -C - opertau_ab(1)
1597 : : ! opertau_ba(tauSamples+1) = -C - opertau_ba(1)
1598 : : !sui!write(6,*) "BackFourier",opertau(tauSamples+1),opertau(1),real(C)
1599 : :
1600 : 0 : op%setT = .TRUE.
1601 : : ! endif
1602 [ # # ]: 0 : op%oper(:,iflavor1,iflavor2)=opertau(:)
1603 : 0 : FREE(opertau)
1604 : : ! FREE(opertau_ab)
1605 : : ! FREE(opertau_ba)
1606 : : END DO ! iflavor2
1607 : : END DO ! iflavor1
1608 : : ! -- End loop over flavors.
1609 : :
1610 [ # # ]: 0 : FREE(Domega)
1611 [ # # ]: 0 : FREE(A_omega)
1612 [ # # ]: 0 : FREE(C_omega)
1613 : 0 : close(236)
1614 : 0 : close(237)
1615 : :
1616 [ # # # # : 0 : END SUBROUTINE GreenHyboffdiagComplex_backFourierComplex
# # ]
1617 : : !!***
1618 : :
1619 : : !!****f* ABINIT/m_GreenHyboffdiagComplex/GreenHyboffdiagComplex_forFourier
1620 : : !! NAME
1621 : : !! GreenHyboffdiagComplex_forFourier
1622 : : !!
1623 : : !! FUNCTION
1624 : : !! perform forward fourier transform without imaginary part of Gtau
1625 : : !!
1626 : : !! COPYRIGHT
1627 : : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
1628 : : !! This file is distributed under the terms of the
1629 : : !! GNU General Public License, see ~abinit/COPYING
1630 : : !! or http://www.gnu.org/copyleft/gpl.txt .
1631 : : !!
1632 : : !! INPUTS
1633 : : !! op=Green
1634 : : !! Wmax=linear maximum frequency
1635 : : !!
1636 : : !! OUTPUT
1637 : : !! Gomega=Results for omega frequencies
1638 : : !! omega=ask frequencies
1639 : : !!
1640 : : !! SIDE EFFECTS
1641 : : !!
1642 : : !! NOTES
1643 : : !!
1644 : : !! SOURCE
1645 : :
1646 [ # # # # : 0 : SUBROUTINE GreenHyboffdiagComplex_forFourier(op, Gomega, omega, Wmax)
# # # # #
# # # ]
1647 : : !Arguments ------------------------------------
1648 : :
1649 : : #ifdef HAVE_MPI1
1650 : : include 'mpif.h'
1651 : : #endif
1652 : : TYPE(GreenHyboffdiagComplex) , INTENT(INOUT) :: op
1653 : : COMPLEX(KIND=8), DIMENSION(:,:,:), OPTIONAL, INTENT(INOUT) :: Gomega ! INOUT for MPI
1654 : : COMPLEX(KIND=8), DIMENSION(:), OPTIONAL, INTENT(IN ) :: omega
1655 : : INTEGER , OPTIONAL, INTENT(IN ) :: Wmax
1656 : : INTEGER :: i
1657 : : INTEGER :: j
1658 : : INTEGER :: iflavor1
1659 : : INTEGER :: iflavor2
1660 : : INTEGER :: nflavors
1661 : : INTEGER :: L
1662 : : INTEGER :: Lspline
1663 : : INTEGER :: Nom
1664 : : INTEGER :: omegaBegin
1665 : : INTEGER :: omegaEnd
1666 : : INTEGER :: deltaw
1667 : : INTEGER :: residu
1668 : 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: counts
1669 : 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: displs
1670 : : DOUBLE PRECISION :: beta
1671 : : DOUBLE PRECISION :: tau
1672 : : DOUBLE PRECISION :: delta
1673 : : DOUBLE PRECISION :: deltabis
1674 : : DOUBLE PRECISION :: inv_delta
1675 : : DOUBLE PRECISION :: inv_delta2
1676 : : DOUBLE PRECISION :: omdeltabis
1677 : : DOUBLE PRECISION :: tmp
1678 : : DOUBLE PRECISION :: xpi
1679 : 0 : DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: diag
1680 : 0 : DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: diagL
1681 : 0 : DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: lastR
1682 : 0 : DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: lastC
1683 : 0 : DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: XM
1684 : 0 : DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: X2
1685 : : DOUBLE PRECISION :: iw
1686 : : COMPLEX(KIND=8) :: iwtau
1687 : 0 : COMPLEX(KIND=8), ALLOCATABLE, DIMENSION(:) :: Gwtmp
1688 : 0 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: omegatmp
1689 : : #if defined HAVE_MPI && !defined HAVE_MPI2_INPLACE
1690 : : INTEGER :: my_count
1691 : : COMPLEX(KIND=8), ALLOCATABLE , DIMENSION(:) :: Gwtmp_buf
1692 : : #endif
1693 : :
1694 : 0 : nflavors=op%nflavors
1695 : :
1696 : : !sui!write(6,*) " Fourier transformation begin"
1697 : :
1698 [ # # ]: 0 : IF ( op%set .EQV. .FALSE. ) &
1699 : 0 : CALL ERROR("GreenHyboffdiag_forFourier : Uninitialized GreenHyboffdiag structure")
1700 [ # # ]: 0 : IF ( op%setT .EQV. .FALSE. ) &
1701 : 0 : CALL ERROR("GreenHyboffdiag_forFourier : no G(tau)")
1702 : : !write(6,*) "op%setMk=", op%setMk
1703 [ # # ]: 0 : IF ( op%setMk .NE. 2*nflavors*nflavors ) &
1704 : 0 : CALL WARNALL("GreenHyboffdiag_forFourier : green does not have moments ")
1705 : :
1706 : 0 : L = op%samples
1707 : :
1708 : 0 : xpi=acos(-1.d0) !!! XPI=PI
1709 : 0 : beta = op%beta
1710 : 0 : Nom = op%Wmax
1711 [ # # ]: 0 : IF ( PRESENT(Gomega) ) THEN
1712 : 0 : Nom = SIZE(Gomega,1)
1713 : : !IF ( op%rank .EQ. 0 ) &
1714 : : !!write(6,*) "size Gomega", Nom
1715 : : END IF
1716 [ # # ]: 0 : IF ( PRESENT(omega) ) THEN
1717 [ # # # # ]: 0 : IF ( PRESENT(Gomega) .AND. SIZE(omega) .NE. Nom ) THEN
1718 : 0 : CALL ERROR("GreenHyboffdiag_forFourier : sizes mismatch ")
1719 : : !ELSE
1720 : : !Nom = SIZE(omega)
1721 : : END IF
1722 : : END IF
1723 [ # # ]: 0 : IF ( .NOT. PRESENT(Gomega) .AND. .NOT. PRESENT(omega) ) THEN
1724 [ # # ]: 0 : IF ( PRESENT(Wmax) ) THEN
1725 : 0 : Nom=Wmax
1726 : : ELSE
1727 : 0 : CALL ERROR("GreenHyboffdiag_forFourier : Missing argument Wmax")
1728 : : END IF
1729 : : END IF
1730 : :
1731 : : !!IF ( ALLOCATED(op%oper_w) ) THEN
1732 : : !! IF ( SIZE(op%oper_w,1) .NE. Nom ) THEN
1733 : : !! FREE(op%oper_w)
1734 : : !! MALLOC(op%oper_w,(1:Nom,nflavors,nflavors))
1735 : : !! END IF
1736 : : !!ELSE
1737 : : !! MALLOC(op%oper_w,(1:Nom,nflavors,nflavors))
1738 : : !!END IF
1739 : :
1740 : : !!write(6,*) "PRESENT(GOMEGA)", PRESENT(GOMEGA)
1741 : : !!write(6,*) "PRESENT(OMEGA)", PRESENT(OMEGA)
1742 : : !call flush(6)
1743 : :
1744 : 0 : delta=op%delta_t
1745 : 0 : inv_delta = op%inv_dt
1746 : 0 : inv_delta2 = inv_delta*inv_delta
1747 : :
1748 [ # # # # : 0 : MALLOC(diagL,(L-1))
# # # # #
# # # #
# ]
1749 [ # # # # : 0 : MALLOC(lastR,(L-1))
# # # # ]
1750 [ # # # # : 0 : MALLOC(diag,(L))
# # # # #
# # # ]
1751 [ # # # # : 0 : MALLOC(lastC,(L-1))
# # # # ]
1752 : :
1753 : : !(cf Stoer) for the spline interpolation :
1754 : : ! second derivatives XM solution of A*XM=B.
1755 : : !A=(2.4.2.11) of Stoer&Bulirsch + 2 limit conditions
1756 : : !The LU decomposition of A is known explicitly;
1757 : :
1758 : 0 : diag (1) = 4.d0 ! 1.d0 *4.d0 factor 4 added for conditionning
1759 : 0 : diagL(1) = 0.25d0 !1.d0/4.d0
1760 : 0 : lastR(1) = -0.5d0 ! -2.d0/4.d0
1761 : 0 : lastC(1) = 4.d0 ! 1.d0*4.d0
1762 : 0 : diag (2) = 4.d0
1763 : 0 : diagL(2) = 0.25d0
1764 : 0 : lastR(2) = -0.25d0
1765 : 0 : lastC(2) = -1.d0
1766 : :
1767 [ # # ]: 0 : DO i = 3, L-2
1768 : 0 : tmp = 4.d0 - diagL(i-1)
1769 : 0 : diagL(i) = 1.d0 / tmp
1770 : : END DO
1771 [ # # ]: 0 : DO i = 3, L-2
1772 : 0 : diag (i) = 1.d0 / diagL(i)
1773 : 0 : lastR(i) = -(lastR(i-1)*diagL(i))
1774 : 0 : lastC(i) = -(lastC(i-1)*diagL(i-1))
1775 : : END DO
1776 : :
1777 : 0 : tmp = 1.d0/diag(L-2)
1778 : 0 : diag (L-1) = 4.d0 - tmp
1779 : 0 : lastR(L-1) = (1.d0 - lastR(L-2))/ diag(L-1)
1780 : : !diagL(L-1) = lastR(L-1)
1781 : 0 : diagL(L-1) = 0.d0 ! for the Lq=B resolution
1782 : : !lastC(L-1) = 1.d0 - lastC(L-2)*diagL(L-1) ! equivalent to the next line
1783 : 0 : lastC(L-1) = 1.d0 - (lastC(L-2)*lastR(L-1)) ! True value
1784 : 0 : diag (L ) = 2.d0! - DOT_PRODUCT( lastR , lastC )
1785 : 0 : tmp = 0.d0
1786 [ # # ]: 0 : DO i = 1, L-1
1787 : 0 : tmp = tmp + lastR(i)*lastC(i)
1788 : : END DO
1789 : 0 : diag (L ) = diag (L ) - tmp
1790 : 0 : lastC(L-1) = lastC(L-1)-1.d0 ! 1 is removed for the u.XM=q resolution
1791 : :
1792 [ # # # # : 0 : MALLOC(XM,(L))
# # # # ]
1793 [ # # # # : 0 : MALLOC(Gwtmp,(1:Nom))
# # # # #
# # # ]
1794 : :
1795 : 0 : Lspline = L-1
1796 [ # # # # : 0 : MALLOC(X2,(1:Lspline+1)) ! We impose L = Nom
# # # # #
# # # ]
1797 : :
1798 [ # # ]: 0 : IF ( op%have_MPI .EQV. .TRUE. ) THEN
1799 : 0 : deltaw = Nom / op%size
1800 : 0 : residu = Nom - op%size*deltaw
1801 [ # # ]: 0 : IF ( op%rank .LT. op%size - residu ) THEN
1802 : 0 : omegaBegin = 1 + op%rank*deltaw
1803 : 0 : omegaEnd = (op%rank + 1)*deltaw
1804 : : ELSE
1805 : : ! tauBegin = (op%size-residu)*deltaw + 1 + (op%rank-op%size+residu)*(deltaw+1)
1806 : 0 : omegaBegin = 1 + op%rank*(deltaw + 1) -op%size + residu
1807 : 0 : omegaEnd = omegaBegin + deltaw
1808 : : END IF
1809 [ # # # # : 0 : MALLOC(counts,(1:op%size))
# # # # #
# # # ]
1810 [ # # # # : 0 : MALLOC(displs,(1:op%size))
# # # # ]
1811 : : counts = (/ (deltaw, i=1, op%size-residu), &
1812 [ # # # # : 0 : (deltaw+1, i=op%size-residu+1, op%size) /)
# # # # #
# # # # #
# # # # #
# ]
1813 : 0 : displs(1)=0
1814 [ # # ]: 0 : DO i = 2, op%size
1815 : 0 : displs(i) = displs(i-1) + counts (i-1)
1816 : : END DO
1817 : : ELSE
1818 : : omegaBegin = 1
1819 : : omegaEnd = Nom
1820 : : END IF
1821 : :
1822 : : ! op%Mk(iflavor1,iflavor2,1) = 0.d0
1823 : : ! DO iflavor1 = 1, nflavors
1824 : : ! op%Mk(iflavor1,iflavor1,1) = -1.d0
1825 : : ! ENDDO
1826 : : ! op%Mk(:,:,3) = 0.d0
1827 : :
1828 [ # # # # : 0 : MALLOC(omegatmp,(omegaBegin:omegaEnd))
# # # # #
# # # ]
1829 [ # # ]: 0 : IF ( PRESENT(omega) ) THEN
1830 [ # # # # : 0 : omegatmp(omegaBegin:omegaEnd) = (/ (AIMAG(omega(i)),i=omegaBegin,omegaEnd) /)
# # ]
1831 : : ELSE
1832 [ # # # # : 0 : omegatmp(omegaBegin:omegaEnd) = (/ ((((2.d0*DBLE(i)-1.d0)*xpi)/Beta), i=omegaBegin,omegaEnd) /)
# # ]
1833 : : END IF
1834 : :
1835 [ # # ]: 0 : DO iflavor1 = 1, nflavors
1836 [ # # ]: 0 : DO iflavor2 = 1, nflavors
1837 : : ! write(6,*) " Moments:",op%Mk(iflavor1,iflavor2,:),iflavor1,iflavor2
1838 : :
1839 : : ! construct the B vector from A.Xm=B
1840 : 0 : XM(1) = 4.d0*op%Mk(iflavor1,iflavor2,3)
1841 : : XM(L) = (6.d0 * inv_delta) * ( op%Mk(iflavor1,iflavor2,2) - ( &
1842 : : (op%oper(2,iflavor1,iflavor2)-op%oper(1,iflavor1,iflavor2)) + &
1843 : 0 : (op%oper(L,iflavor1,iflavor2)-op%oper(L-1,iflavor1,iflavor2)) ) * inv_delta )
1844 : : ! built generic second derivative of oper
1845 : : !sui!write(6,*) "XM 1 L",XM(1),XM(L),op%Mk(iflavor1,iflavor2,2),op%Mk(iflavor1,iflavor2,3)
1846 [ # # ]: 0 : DO i = 2, L-1
1847 : : XM(i) = (6.d0 * inv_delta2) * ( (op%oper(i+1,iflavor1,iflavor2) &
1848 : : - 2.d0 * op%oper(i,iflavor1,iflavor2)) &
1849 : 0 : + op%oper(i-1,iflavor1,iflavor2) )
1850 : : !sui!write(6,*) "XM",i,XM(i),op%oper(i,iflavor1,iflavor2)
1851 : : END DO
1852 : :
1853 : : ! Find second derivatives XM: Solve the system
1854 : : ! SOLVING Lq= XM
1855 : : ! q = XM
1856 [ # # ]: 0 : do j=1,L-1
1857 : 0 : XM(j+1)=XM(j+1)-(diagL(j)*XM(j))
1858 : 0 : XM(L) =XM(L) -(lastR(j)*XM(j))
1859 : : end do
1860 : :
1861 : :
1862 : : ! SOLVING U.XM=q
1863 : : ! XM = q
1864 [ # # ]: 0 : do j=L-1,2,-1
1865 : 0 : XM(j+1) = XM(j+1) / diag(j+1)
1866 : 0 : XM(j)= (XM(j)-(XM(L)*lastC(j)))-XM(j+1)
1867 : : end do
1868 : 0 : XM(2) = XM(2) / diag(2)
1869 : 0 : XM(1) = (XM(1)-XM(L)*lastC(1)) / diag(1)
1870 : :
1871 : :
1872 : :
1873 : : !Construct L2 second derivative from known derivatives XM
1874 : 0 : deltabis = beta / DBLE(Lspline)
1875 [ # # ]: 0 : DO i = 1, Lspline
1876 : 0 : tau = deltabis * DBLE(i-1)
1877 : 0 : j = ((L-1)*(i-1))/Lspline + 1!INT(tau * inv_delta) + 1
1878 : 0 : X2(i) = inv_delta * ( XM(j)*(DBLE(j)*delta - tau ) + XM(j+1)*(tau - DBLE(j-1)*delta) )
1879 : : END DO
1880 : 0 : X2(Lspline+1) = XM(L)
1881 : :
1882 : :
1883 [ # # ]: 0 : DO i = omegaBegin, omegaEnd
1884 : 0 : iw = omegatmp(i)
1885 : 0 : omdeltabis = iw*deltabis
1886 : 0 : Gwtmp(i)=CMPLX(0.d0,0.d0,8)
1887 [ # # ]: 0 : DO j=2, Lspline ! We impose L+1 = Nom
1888 : 0 : iwtau = CMPLX(0.d0,omdeltabis*DBLE(j-1),8)
1889 : 0 : Gwtmp(i) = Gwtmp(i) + EXP(iwtau) * CMPLX((X2(j+1) + X2(j-1))-2.d0*X2(j),0.d0,8)
1890 : : !write(6,*) "ww",i,j,Gwtmp(i),X2(j),iwtau
1891 : : END DO
1892 : : Gwtmp(i) = Gwtmp(i)/CMPLX(((iw*iw)*(iw*iw)*deltabis),0.d0,8) &
1893 : : + CMPLX( ( ((X2(2)-X2(1))+(X2(Lspline+1)-X2(Lspline)))/((iw*iw)*deltabis) -real(op%Mk(iflavor1,iflavor2,2)) ) &
1894 : 0 : /(iw*iw) , real(op%Mk(iflavor1,iflavor2,1)-op%Mk(iflavor1,iflavor2,3)/(iw*iw))/iw , 8)
1895 : : !+ CMPLX( (X2(2)-X2(1))+(X2(Lspline+1)-X2(Lspline)), 0.d0, 8 ) ) &
1896 : : ! / (((iw*iw)*(iw*iw))*CMPLX(deltabis,0.d0,8)) &
1897 : : !- CMPLX(op%Mk(1),0.d0,8)/iw &
1898 : : !+ CMPLX(op%Mk(2),0.d0,8)/(iw*iw) &
1899 : : !- CMPLX(op%Mk(3),0.d0,8)/((iw*iw)*iw)
1900 : : !IF ( op%rank .EQ. 0 ) write(12819,*) iw,gwtmp(i)
1901 : : END DO
1902 : : !call flush(12819)
1903 [ # # ]: 0 : IF ( op%have_MPI .EQV. .TRUE. ) THEN
1904 : : #ifdef HAVE_MPI
1905 : : #if defined HAVE_MPI2_INPLACE
1906 : : CALL MPI_ALLGATHERV(MPI_IN_PLACE, 0, MPI_DOUBLE_COMPLEX, &
1907 : : Gwtmp , counts, displs, &
1908 : 0 : MPI_DOUBLE_COMPLEX, op%MY_COMM, residu)
1909 : : #else
1910 : : my_count=omegaBegin-omegaEnd+1
1911 : : MALLOC(Gwtmp_buf,(my_count))
1912 : : Gwtmp_buf(1:my_count)=Gwtmp(omegaBegin:omegaEnd)
1913 : : CALL MPI_ALLGATHERV(Gwtmp_buf, my_count, MPI_DOUBLE_COMPLEX, &
1914 : : Gwtmp , counts, displs, &
1915 : : MPI_DOUBLE_COMPLEX, op%MY_COMM, residu)
1916 : : FREE(Gwtmp_buf)
1917 : : #endif
1918 : : #endif
1919 : : END IF
1920 [ # # ]: 0 : IF ( PRESENT(Gomega) ) THEN
1921 [ # # ]: 0 : Gomega(:,iflavor1,iflavor2) = Gwtmp(:)
1922 : : END IF
1923 : 0 : op%setW = .TRUE.
1924 : : ENDDO ! iflavor1
1925 : : ENDDO ! iflavor2
1926 : : !!op%oper_w=Gomega
1927 : : do iflavor1=1,nflavors
1928 : : !sui!write(6,*) iflavor1
1929 : : do i=1,Nom
1930 : : !write(6,*) "w",i,op%oper_w(i,iflavor1,iflavor1)
1931 : : enddo
1932 : : enddo
1933 : :
1934 [ # # ]: 0 : FREE(Gwtmp)
1935 : 0 : FREE(diagL)
1936 : 0 : FREE(lastR)
1937 : 0 : FREE(diag)
1938 : 0 : FREE(lastC)
1939 [ # # ]: 0 : FREE(XM)
1940 [ # # ]: 0 : FREE(omegatmp)
1941 [ # # ]: 0 : FREE(X2)
1942 [ # # ]: 0 : FREE(counts)
1943 [ # # ]: 0 : FREE(displs)
1944 : :
1945 : 0 : END SUBROUTINE GreenHyboffdiagComplex_forFourier
1946 : : !!***
1947 : :
1948 : : !!****f* ABINIT/m_GreenHyboffdiagComplex/GreenHyboffdiagComplex_forFourierComplex
1949 : : !! NAME
1950 : : !! GreenHyboffdiagComplex_forFourier
1951 : : !!
1952 : : !! FUNCTION
1953 : : !! perform forward fourier transform
1954 : : !!
1955 : : !! COPYRIGHT
1956 : : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
1957 : : !! This file is distributed under the terms of the
1958 : : !! GNU General Public License, see ~abinit/COPYING
1959 : : !! or http://www.gnu.org/copyleft/gpl.txt .
1960 : : !!
1961 : : !! INPUTS
1962 : : !! op=Green
1963 : : !! Wmax=linear maximum frequency
1964 : : !!
1965 : : !! OUTPUT
1966 : : !! Gomega=Results for omega frequencies
1967 : : !! omega=ask frequencies
1968 : : !!
1969 : : !! SIDE EFFECTS
1970 : : !!
1971 : : !! NOTES
1972 : : !!
1973 : : !! SOURCE
1974 : :
1975 [ # # # # : 0 : SUBROUTINE GreenHyboffdiagComplex_forFourierComplex(op, Gomega, omega, Wmax)
# # # # #
# # # ]
1976 : : !Arguments ------------------------------------
1977 : :
1978 : : #ifdef HAVE_MPI1
1979 : : include 'mpif.h'
1980 : : #endif
1981 : : TYPE(GreenHyboffdiagComplex) , INTENT(INOUT) :: op
1982 : : COMPLEX(KIND=8), DIMENSION(:,:,:), OPTIONAL, INTENT(INOUT) :: Gomega ! INOUT for MPI
1983 : : COMPLEX(KIND=8), DIMENSION(:), OPTIONAL, INTENT(IN ) :: omega
1984 : : INTEGER , OPTIONAL, INTENT(IN ) :: Wmax
1985 : : INTEGER :: i
1986 : : INTEGER :: j
1987 : : INTEGER :: iflavor1
1988 : : INTEGER :: iflavor2
1989 : : INTEGER :: nflavors
1990 : : INTEGER :: L
1991 : : INTEGER :: Lspline
1992 : : INTEGER :: Nom
1993 : : INTEGER :: omegaBegin
1994 : : INTEGER :: omegaEnd
1995 : : INTEGER :: deltaw
1996 : : INTEGER :: residu
1997 : 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: counts
1998 : 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: displs
1999 : : DOUBLE PRECISION :: beta
2000 : : DOUBLE PRECISION :: tau
2001 : : DOUBLE PRECISION :: delta
2002 : : DOUBLE PRECISION :: deltabis
2003 : : DOUBLE PRECISION :: inv_delta
2004 : : DOUBLE PRECISION :: inv_delta2
2005 : : DOUBLE PRECISION :: omdeltabis
2006 : : DOUBLE PRECISION :: tmp
2007 : : DOUBLE PRECISION :: xpi
2008 : 0 : DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: diag
2009 : 0 : DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: diagL
2010 : 0 : DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: lastR
2011 : 0 : DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: lastC
2012 : 0 : DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: XM
2013 : 0 : COMPLEX(KIND=8), DIMENSION(:), ALLOCATABLE :: XMC
2014 : 0 : DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: X2
2015 : 0 : COMPLEX(KIND=8), DIMENSION(:), ALLOCATABLE :: X2C
2016 : : DOUBLE PRECISION :: iw
2017 : : COMPLEX(KIND=8) :: iwtau
2018 : 0 : COMPLEX(KIND=8), ALLOCATABLE, DIMENSION(:) :: Gwtmp
2019 : 0 : COMPLEX(KIND=8), ALLOCATABLE, DIMENSION(:) :: Gwtmpc
2020 : 0 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: omegatmp
2021 : :
2022 : : #if defined HAVE_MPI && !defined HAVE_MPI2_INPLACE
2023 : : INTEGER :: my_count
2024 : : COMPLEX(KIND=8), ALLOCATABLE , DIMENSION(:) :: Gwtmp_buf
2025 : : #endif
2026 : :
2027 : 0 : nflavors=op%nflavors
2028 : :
2029 : : !sui!write(6,*) " Fourier transformation begin"
2030 : :
2031 [ # # ]: 0 : IF ( op%set .EQV. .FALSE. ) &
2032 : 0 : CALL ERROR("GreenHyboffdiagComplex_forFourier : Uninitialized GreenHyboffdiagComplex structure")
2033 [ # # ]: 0 : IF ( op%setT .EQV. .FALSE. ) &
2034 : 0 : CALL ERROR("GreenHyboffdiagComplex_forFourier : no G(tau)")
2035 : : !write(6,*) "op%setMk=", op%setMk
2036 [ # # ]: 0 : IF ( op%setMk .NE. 2*nflavors*nflavors ) &
2037 : 0 : CALL WARNALL("GreenHyboffdiagComplex_forFourier : green does not have moments ")
2038 : :
2039 : 0 : L = op%samples
2040 : :
2041 : 0 : xpi=acos(-1.d0) !!! XPI=PI
2042 : 0 : beta = op%beta
2043 : 0 : Nom = op%Wmax
2044 : :
2045 [ # # ]: 0 : IF ( PRESENT(Gomega) ) THEN
2046 : 0 : Nom = SIZE(Gomega,1)
2047 : : !IF ( op%rank .EQ. 0 ) &
2048 : : !!write(6,*) "size Gomega", Nom
2049 : : END IF
2050 [ # # ]: 0 : IF ( PRESENT(omega) ) THEN
2051 [ # # # # ]: 0 : IF ( PRESENT(Gomega) .AND. SIZE(omega) .NE. Nom ) THEN
2052 : 0 : CALL ERROR("GreenHyboffdiagComplex_forFourier : sizes mismatch ")
2053 : : !ELSE
2054 : : !Nom = SIZE(omega)
2055 : : END IF
2056 : : END IF
2057 [ # # ]: 0 : IF ( .NOT. PRESENT(Gomega) .AND. .NOT. PRESENT(omega) ) THEN
2058 [ # # ]: 0 : IF ( PRESENT(Wmax) ) THEN
2059 : 0 : Nom=Wmax
2060 : : ELSE
2061 : 0 : CALL ERROR("GreenHyboffdiagComplex_forFourier : Missing argument Wmax")
2062 : : END IF
2063 : : END IF
2064 : :
2065 : : !!IF ( ALLOCATED(op%oper_w) ) THEN
2066 : : !! IF ( SIZE(op%oper_w,1) .NE. Nom ) THEN
2067 : : !! FREE(op%oper_w)
2068 : : !! MALLOC(op%oper_w,(1:Nom,nflavors,nflavors))
2069 : : !! END IF
2070 : : !!ELSE
2071 : : !! MALLOC(op%oper_w,(1:Nom,nflavors,nflavors))
2072 : : !!END IF
2073 : :
2074 : : !!write(6,*) "PRESENT(GOMEGA)", PRESENT(GOMEGA)
2075 : : !!write(6,*) "PRESENT(OMEGA)", PRESENT(OMEGA)
2076 : : !call flush(6)
2077 : :
2078 : 0 : delta=op%delta_t
2079 : 0 : inv_delta = op%inv_dt
2080 : 0 : inv_delta2 = inv_delta*inv_delta
2081 : :
2082 [ # # # # : 0 : MALLOC(diagL,(L-1))
# # # # #
# # # #
# ]
2083 [ # # # # : 0 : MALLOC(lastR,(L-1))
# # # # ]
2084 [ # # # # : 0 : MALLOC(diag,(L))
# # # # #
# # # ]
2085 [ # # # # : 0 : MALLOC(lastC,(L-1))
# # # # ]
2086 : :
2087 : : !(cf Stoer) for the spline interpolation :
2088 : : ! second derivatives XM solution of A*XM=B.
2089 : : !A=(2.4.2.11) of Stoer&Bulirsch + 2 limit conditions
2090 : : !The LU decomposition of A is known explicitly;
2091 : :
2092 : : !Construction of LU matrix elements (see Appendix Thesis J. Bieder)
2093 : 0 : diag (1) = 4.d0 ! 1.d0 *4.d0 factor 4 added for conditionning
2094 : 0 : diagL(1) = 0.25d0 !1.d0/4.d0
2095 : 0 : lastR(1) = -0.5d0 ! -2.d0/4.d0
2096 : 0 : lastC(1) = 4.d0 ! 1.d0*4.d0
2097 : 0 : diag (2) = 4.d0
2098 : 0 : diagL(2) = 0.25d0
2099 : 0 : lastR(2) = -0.25d0
2100 : 0 : lastC(2) = -1.d0
2101 : :
2102 : : ! -- Fill matrix (LU decomposition)
2103 [ # # ]: 0 : DO i = 3, L-2
2104 : 0 : tmp = 4.d0 - diagL(i-1)
2105 : 0 : diagL(i) = 1.d0 / tmp
2106 : : END DO
2107 [ # # ]: 0 : DO i = 3, L-2
2108 : 0 : diag (i) = 1.d0 / diagL(i)
2109 : 0 : lastR(i) = -(lastR(i-1)*diagL(i))
2110 : 0 : lastC(i) = -(lastC(i-1)*diagL(i-1))
2111 : : END DO
2112 : :
2113 : : ! -- Boundary conditions
2114 : 0 : tmp = 1.d0/diag(L-2)
2115 : 0 : diag (L-1) = 4.d0 - tmp
2116 : 0 : lastR(L-1) = (1.d0 - lastR(L-2))/ diag(L-1)
2117 : : !diagL(L-1) = lastR(L-1)
2118 : 0 : diagL(L-1) = 0.d0 ! for the Lq=B resolution
2119 : : !lastC(L-1) = 1.d0 - lastC(L-2)*diagL(L-1) ! equivalent to the next line
2120 : 0 : lastC(L-1) = 1.d0 - (lastC(L-2)*lastR(L-1)) ! True value
2121 : 0 : diag (L ) = 2.d0! - DOT_PRODUCT( lastR , lastC )
2122 : 0 : tmp = 0.d0
2123 [ # # ]: 0 : DO i = 1, L-1
2124 : 0 : tmp = tmp + lastR(i)*lastC(i)
2125 : : END DO
2126 : 0 : diag (L ) = diag (L ) - tmp
2127 : 0 : lastC(L-1) = lastC(L-1)-1.d0 ! 1 is removed for the u.XM=q resolution
2128 : : !end construction of LU matrix elements
2129 : :
2130 [ # # # # : 0 : MALLOC(XM,(L))
# # # # ]
2131 [ # # # # : 0 : MALLOC(XMC,(L))
# # # # #
# # # ]
2132 [ # # # # : 0 : MALLOC(Gwtmp,(1:Nom))
# # # # #
# # # ]
2133 [ # # # # : 0 : MALLOC(Gwtmpc,(1:Nom))
# # # # ]
2134 : :
2135 : 0 : Lspline = L-1
2136 [ # # # # : 0 : MALLOC(X2,(1:Lspline+1)) ! We impose L = Nom
# # # # #
# # # ]
2137 [ # # # # : 0 : MALLOC(X2C,(1:Lspline+1))
# # # # #
# # # ]
2138 : :
2139 [ # # ]: 0 : IF ( op%have_MPI .EQV. .TRUE. ) THEN
2140 : 0 : deltaw = Nom / op%size
2141 : 0 : residu = Nom - op%size*deltaw
2142 [ # # ]: 0 : IF ( op%rank .LT. op%size - residu ) THEN
2143 : 0 : omegaBegin = 1 + op%rank*deltaw
2144 : 0 : omegaEnd = (op%rank + 1)*deltaw
2145 : : ELSE
2146 : : ! tauBegin = (op%size-residu)*deltaw + 1 + (op%rank-op%size+residu)*(deltaw+1)
2147 : 0 : omegaBegin = 1 + op%rank*(deltaw + 1) -op%size + residu
2148 : 0 : omegaEnd = omegaBegin + deltaw
2149 : : END IF
2150 [ # # # # : 0 : MALLOC(counts,(1:op%size))
# # # # #
# # # ]
2151 [ # # # # : 0 : MALLOC(displs,(1:op%size))
# # # # ]
2152 : : counts = (/ (deltaw, i=1, op%size-residu), &
2153 [ # # # # : 0 : (deltaw+1, i=op%size-residu+1, op%size) /)
# # # # #
# # # # #
# # # # #
# ]
2154 : 0 : displs(1)=0
2155 [ # # ]: 0 : DO i = 2, op%size
2156 : 0 : displs(i) = displs(i-1) + counts (i-1)
2157 : : END DO
2158 : : ELSE
2159 : : omegaBegin = 1
2160 : : omegaEnd = Nom
2161 : : END IF
2162 : :
2163 [ # # # # : 0 : MALLOC(omegatmp,(omegaBegin:omegaEnd))
# # # # #
# # # ]
2164 [ # # ]: 0 : IF ( PRESENT(omega) ) THEN
2165 [ # # # # : 0 : omegatmp(omegaBegin:omegaEnd) = (/ (AIMAG(omega(i)),i=omegaBegin,omegaEnd) /)
# # ]
2166 : : ELSE
2167 [ # # # # : 0 : omegatmp(omegaBegin:omegaEnd) = (/ ((((2.d0*DBLE(i)-1.d0)*xpi)/Beta), i=omegaBegin,omegaEnd) /)
# # ]
2168 : : END IF
2169 : :
2170 : : ! -- Main loop over flavors --
2171 [ # # ]: 0 : DO iflavor1 = 1, nflavors
2172 [ # # ]: 0 : DO iflavor2 = 1, nflavors
2173 : : ! write(6,*) " Moments:",op%Mk(iflavor1,iflavor2,:),iflavor1,iflavor2
2174 : : ! construct the B vector from A.Xm=B
2175 : 0 : XM(1) = 4.d0*real(op%Mk(iflavor1,iflavor2,3))
2176 : : XM(L) = (6.d0 * inv_delta) * ( real(op%Mk(iflavor1,iflavor2,2)) - ( &
2177 : : (real(op%oper(2,iflavor1,iflavor2))-real(op%oper(1,iflavor1,iflavor2))) + &
2178 : 0 : (real(op%oper(L,iflavor1,iflavor2))-real(op%oper(L-1,iflavor1,iflavor2))) ) * inv_delta )
2179 : : ! build d_ni terms of B
2180 [ # # ]: 0 : DO i = 2, L-1
2181 : : XM(i) = (6.d0 * inv_delta2) * ( (real(op%oper(i+1,iflavor1,iflavor2)) &
2182 : : - 2.d0 * real(op%oper(i,iflavor1,iflavor2))) &
2183 : 0 : + real(op%oper(i-1,iflavor1,iflavor2)) )
2184 : : END DO
2185 : : !complex version
2186 : 0 : XMC(1) = 4.d0*op%Mk(iflavor1,iflavor2,3)
2187 : : XMC(L) = (6.d0 * inv_delta) * ( op%Mk(iflavor1,iflavor2,2) - ( &
2188 : : (op%oper(2,iflavor1,iflavor2)-op%oper(1,iflavor1,iflavor2)) + &
2189 : 0 : (op%oper(L,iflavor1,iflavor2)-op%oper(L-1,iflavor1,iflavor2)) ) * inv_delta)
2190 : : ! build d_ni terms of B
2191 [ # # ]: 0 : DO i = 2, L-1
2192 : : XMC(i) = (6.d0 * inv_delta2) * ( (op%oper(i+1,iflavor1,iflavor2) &
2193 : : - 2.d0 * op%oper(i,iflavor1,iflavor2)) &
2194 : 0 : + op%oper(i-1,iflavor1,iflavor2))
2195 : : END DO
2196 : :
2197 : : ! Find second derivatives XM: Solve the system
2198 : : ! SOLVING Lq= XM
2199 : : ! q = XM
2200 [ # # ]: 0 : do j=1,L-1
2201 : 0 : XM(j+1)=XM(j+1)-(diagL(j)*XM(j))
2202 : 0 : XM(L) =XM(L) -(lastR(j)*XM(j))
2203 : : !complex version
2204 : : !write(*,*)"diagL(j) lastR(j)",diagL(j),lastR(j)
2205 : 0 : XMC(j+1)=XMC(j+1)-(cmplx(diagL(j),0.d0,kind=8)*XMC(j))
2206 : 0 : XMC(L) =XMC(L) -(cmplx(lastR(j),0.d0,kind=8)*XMC(j))
2207 : : !write(*,*)"XM XMC",XM(j),XMC(j)
2208 : : end do
2209 : :
2210 : :
2211 : : ! SOLVING U.XM=q
2212 : : ! XM = q
2213 [ # # ]: 0 : do j=L-1,2,-1
2214 : 0 : XM(j+1) = XM(j+1) / diag(j+1)
2215 : 0 : XM(j)= (XM(j)-(XM(L)*lastC(j)))-XM(j+1)
2216 : : !complex version
2217 : 0 : XMC(j+1) = XMC(j+1) / cmplx(diag(j+1),0.d0,kind=8)
2218 : 0 : XMC(j)= (XMC(j)-(XMC(L)*cmplx(lastC(j),0.d0,kind=8)))-XMC(j+1)
2219 : : end do
2220 : 0 : XM(2) = XM(2) / diag(2)
2221 : 0 : XM(1) = (XM(1)-XM(L)*lastC(1)) / diag(1)
2222 : : !version complex
2223 : 0 : XMC(2) = XMC(2) / cmplx(diag(2),0.d0,kind=8)
2224 : 0 : XMC(1) = (XMC(1)-XMC(L)*cmplx(lastC(1),0.d0,kind=8)) / cmplx(diag(1),0.d0,kind=8)
2225 : : !write(*,*) "XM(1) XMC(1) XM(2) XMC(2)",XM(1),XMC(1),XM(2),XMC(2)
2226 : :
2227 : : !Construct L2 second derivative from known derivatives XM
2228 : 0 : deltabis = beta / DBLE(Lspline)
2229 [ # # ]: 0 : DO i = 1, Lspline
2230 : 0 : tau = deltabis * DBLE(i-1)
2231 : 0 : j = ((L-1)*(i-1))/Lspline + 1!INT(tau * inv_delta) + 1
2232 : 0 : X2(i) = inv_delta * ( XM(j)*(DBLE(j)*delta - tau ) + XM(j+1)*(tau - DBLE(j-1)*delta) )
2233 : : !version complex
2234 : 0 : X2C(i) = cmplx(inv_delta,0.d0,kind=8) * ( XMC(j)*(DBLE(j)*delta - tau ) + XMC(j+1)*(tau - DBLE(j-1)*delta) )
2235 : : END DO
2236 : 0 : X2(Lspline+1) = XM(L)
2237 : 0 : X2C(Lspline+1) = XMC(L)
2238 : : !version complex
2239 : : !write(*,*) "X2(1) X2C(1) X2(2) X2C(2)",X2(1),X2C(1),X2(2),X2C(2)
2240 : :
2241 : : ! -- Fourier Transform --
2242 [ # # ]: 0 : DO i = omegaBegin, omegaEnd
2243 : 0 : iw = omegatmp(i)
2244 : 0 : omdeltabis = iw*deltabis
2245 : : !Gwtmp(i)=CMPLX(0.d0,0.d0,8)
2246 : 0 : Gwtmpc(i)=cmplx(0.d0,0.d0,kind=8)
2247 [ # # ]: 0 : DO j=2, Lspline ! We impose L+1 = Nom
2248 : 0 : iwtau = CMPLX(0.d0,omdeltabis*DBLE(j-1),8)
2249 : : !Gwtmp(i) = Gwtmp(i) + EXP(iwtau) * CMPLX((X2(j+1) + X2(j-1))-2.d0*X2(j),0.d0,kind=8)
2250 : : !write(6,*) "ww",i,j,Gwtmp(i),X2(j),iwtau
2251 : : !version complex
2252 : 0 : Gwtmpc(i) = Gwtmpc(i) + EXP(iwtau) * (X2C(j+1) + X2C(j-1)-2.d0*X2C(j))
2253 : : !write(*,*)"Gwtmp(l) Gwtmpc(l)",Gwtmp(i),Gwtmpc(i)
2254 : : END DO
2255 : : !Gwtmp(i) = Gwtmp(i)/CMPLX(((iw*iw)*(iw*iw)*deltabis),0.d0,8) &
2256 : : ! + CMPLX( ( ((X2(2)-X2(1))+(X2(Lspline+1)-X2(Lspline)))/((iw*iw*iw*iw)*deltabis)),0.d0,kind=8) & !C4/iomega^4
2257 : : !-cmplx(real(op%Mk(iflavor1,iflavor2,2))/(iw*iw),0.d0,kind=8) & !C2/iomega^2
2258 : : !& +cmplx(0.d0,(real(op%Mk(iflavor1,iflavor2,1))-real(op%Mk(iflavor1,iflavor2,3))/(iw*iw))/iw,kind=8) !C3/iomega^3
2259 : : !
2260 : : ! Complex version
2261 : : !
2262 : : Gwtmpc(i) = Gwtmpc(i)/((iw*iw)*(iw*iw)*deltabis) & ! int_O^\beta G^(4)/iomega^4
2263 : : & + ((X2C(2)-X2C(1))+(X2C(Lspline+1)-X2C(Lspline)))/((iw*iw*iw*iw)*deltabis) & ! C4/iomega^4
2264 : : & - op%Mk(iflavor1,iflavor2,2)/(iw*iw) & ! C2/iomega^2
2265 : : & - op%Mk(iflavor1,iflavor2,1)/cmplx(0.d0,iw,kind=8) & ! C1/iomega (careful with sign)
2266 : 0 : & - op%Mk(iflavor1,iflavor2,3)/(cmplx(0.d0,iw,kind=8)*cmplx(0.d0,iw,kind=8)*cmplx(0.d0,iw,kind=8)) ! C3/iomega^3
2267 : : !
2268 : : !IF ( op%rank .EQ. 0 ) write(*,*)"iw,Gwtmp(iw),Gwtmpc(iw)",iw,Gwtmp(i),Gwtmpc(i)
2269 : : END DO
2270 : : !call flush(12819)
2271 [ # # ]: 0 : IF ( op%have_MPI .EQV. .TRUE. ) THEN
2272 : : #ifdef HAVE_MPI
2273 : : #if defined HAVE_MPI2_INPLACE
2274 : : CALL MPI_ALLGATHERV(MPI_IN_PLACE, 0, MPI_DOUBLE_COMPLEX, &
2275 : : Gwtmpc , counts, displs, &
2276 : 0 : MPI_DOUBLE_COMPLEX, op%MY_COMM, residu)
2277 : : #else
2278 : : my_count=omegaBegin-omegaEnd+1
2279 : : MALLOC(Gwtmp_buf,(my_count))
2280 : : Gwtmp_buf(1:my_count)=Gwtmpc(omegaBegin:omegaEnd)
2281 : : CALL MPI_ALLGATHERV(Gwtmp_buf, my_count, MPI_DOUBLE_COMPLEX, &
2282 : : Gwtmpc , counts, displs, &
2283 : : MPI_DOUBLE_COMPLEX, op%MY_COMM, residu)
2284 : : FREE(Gwtmp_buf)
2285 : : #endif
2286 : : #endif
2287 : : END IF
2288 [ # # ]: 0 : IF ( PRESENT(Gomega) ) THEN
2289 [ # # ]: 0 : Gomega(:,iflavor1,iflavor2) = Gwtmpc(:)
2290 : : END IF
2291 : 0 : op%setW = .TRUE.
2292 : : ENDDO ! iflavor1
2293 : : ENDDO ! iflavor2
2294 : :
2295 [ # # ]: 0 : FREE(Gwtmp)
2296 [ # # ]: 0 : FREE(Gwtmpc)
2297 : 0 : FREE(diagL)
2298 : 0 : FREE(lastR)
2299 : 0 : FREE(diag)
2300 : 0 : FREE(lastC)
2301 [ # # ]: 0 : FREE(XM)
2302 [ # # ]: 0 : FREE(XMC)
2303 [ # # ]: 0 : FREE(omegatmp)
2304 [ # # ]: 0 : FREE(X2)
2305 [ # # ]: 0 : FREE(X2C)
2306 [ # # ]: 0 : FREE(counts)
2307 [ # # ]: 0 : FREE(displs)
2308 : :
2309 : 0 : END SUBROUTINE GreenHyboffdiagComplex_forFourierComplex
2310 : : !!***
2311 : :
2312 : : !!****f* ABINIT/m_GreenHyboffdiagComplex/GreenHyboffdiagComplex_print
2313 : : !! NAME
2314 : : !! GreenHyboffdiagComplex_print
2315 : : !!
2316 : : !! FUNCTION
2317 : : !! print Green function
2318 : : !!
2319 : : !! COPYRIGHT
2320 : : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
2321 : : !! This file is distributed under the terms of the
2322 : : !! GNU General Public License, see ~abinit/COPYING
2323 : : !! or http://www.gnu.org/copyleft/gpl.txt .
2324 : : !!
2325 : : !! INPUTS
2326 : : !! op=Green
2327 : : !! ostream=file stream
2328 : : !!
2329 : : !! OUTPUT
2330 : : !!
2331 : : !! SIDE EFFECTS
2332 : : !!
2333 : : !! NOTES
2334 : : !!
2335 : : !! SOURCE
2336 : :
2337 : 0 : SUBROUTINE GreenHyboffdiagComplex_print(op, ostream)
2338 : :
2339 : : !Arguments ------------------------------------
2340 : : TYPE(GreenHyboffdiagComplex), INTENT(IN) :: op
2341 : : INTEGER, OPTIONAL , INTENT(IN) :: ostream
2342 : : !Local variables ------------------------------
2343 : : INTEGER :: ostream_val
2344 : : INTEGER :: isample
2345 : : INTEGER :: samples
2346 : : INTEGER :: iflavor1
2347 : : INTEGER :: iflavor2
2348 : :
2349 : :
2350 [ # # ]: 0 : IF ( op%set .EQV. .FALSE. ) &
2351 : 0 : CALL ERROR("GreenHyboffdiagComplex_print : green op%operator not set ")
2352 : :
2353 [ # # ]: 0 : IF ( PRESENT(ostream) ) THEN
2354 : 0 : ostream_val = ostream
2355 : : ELSE
2356 : 0 : ostream_val = 66
2357 : 0 : OPEN(UNIT=ostream_val,FILE="Green.dat")
2358 : : END IF
2359 : :
2360 : 0 : samples = op%samples
2361 : :
2362 [ # # ]: 0 : DO iflavor1=1,op%nflavors
2363 [ # # ]: 0 : DO iflavor2=1,op%nflavors
2364 : 0 : WRITE(ostream_val,'(a,i3,a,i3,a)') "## (iflavor1,iflavor2)= (", iflavor1,",",iflavor2,")"
2365 [ # # ]: 0 : DO isample = 1, samples
2366 : 0 : WRITE(ostream_val,*) DBLE(isample-1)*op%delta_t, op%oper(isample,iflavor1,iflavor2)
2367 : : END DO
2368 : 0 : WRITE(ostream_val,*)
2369 : : END DO
2370 : : END DO
2371 : :
2372 [ # # ]: 0 : IF ( .NOT. PRESENT(ostream) ) &
2373 : 0 : CLOSE(ostream_val)
2374 : 0 : END SUBROUTINE GreenHyboffdiagComplex_print
2375 : : !!***
2376 : :
2377 : : !!****f* ABINIT/m_GreenHyboffdiagComplex/GreenHyboffdiagComplex_destroy
2378 : : !! NAME
2379 : : !! GreenHyboffdiagComplex_destroy
2380 : : !!
2381 : : !! FUNCTION
2382 : : !! destroy green function
2383 : : !!
2384 : : !! COPYRIGHT
2385 : : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
2386 : : !! This file is distributed under the terms of the
2387 : : !! GNU General Public License, see ~abinit/COPYING
2388 : : !! or http://www.gnu.org/copyleft/gpl.txt .
2389 : : !!
2390 : : !! INPUTS
2391 : : !! op=Green
2392 : : !!
2393 : : !! OUTPUT
2394 : : !!
2395 : : !! SIDE EFFECTS
2396 : : !!
2397 : : !! NOTES
2398 : : !!
2399 : : !! SOURCE
2400 : :
2401 : 0 : SUBROUTINE GreenHyboffdiagComplex_destroy(op)
2402 : :
2403 : : !Arguments ------------------------------------
2404 : : TYPE(GreenHyboffdiagComplex), INTENT(INOUT) :: op
2405 : : INTEGER :: iflavor,iflavorbis
2406 : :
2407 : 0 : op%set = .FALSE.
2408 : 0 : op%setT = .FALSE.
2409 : 0 : op%setW = .FALSE.
2410 : 0 : op%samples = 0
2411 : 0 : op%measurements = 0
2412 : 0 : op%beta = 0.d0
2413 : 0 : op%inv_beta = 0.d0
2414 : 0 : op%inv_dt = 0.d0
2415 : 0 : op%delta_t = 0.d0
2416 : 0 : CALL VectorInt_destroy(op%index_old)
2417 : 0 : CALL VectorComplex_destroy(op%oper_old)
2418 [ # # ]: 0 : do iflavor=1,op%nflavors
2419 [ # # ]: 0 : do iflavorbis=1,op%nflavors
2420 : : !sui!write(6,*) "test",iflavor,iflavorbis
2421 : 0 : CALL MapHybComplex_destroy(op%map(iflavor,iflavorbis))
2422 : : enddo
2423 : : enddo
2424 [ # # # # : 0 : DT_FREEIF(op%map)
# # # # #
# ]
2425 [ # # ]: 0 : FREEIF(op%oper)
2426 [ # # ]: 0 : FREEIF(op%Mk)
2427 [ # # ]: 0 : FREEIF(op%oper_w)
2428 [ # # ]: 0 : FREEIF(op%oper_w_old)
2429 [ # # ]: 0 : FREEIF(op%omega)
2430 : 0 : END SUBROUTINE GreenHyboffdiagComplex_destroy
2431 : :
2432 : : !!***
2433 : : ! This routine contains direct and inverse fourier transformation
2434 : : ! It is a modification of a routine of the GNU GPL
2435 : : ! code available on http://dmft.rutgers.edu/ and
2436 : : ! described in the RMP 2006 paper written by
2437 : : ! G.Kotliar, S.Y.Savrasov, K.Haule, V.S.Oudovenko, O.Parcollet, C.A.Marianetti
2438 : : !=======+=========+=========+=========+=========+=========+=========+=$
2439 : : ! TYPE : SUBROUTINE
2440 : : ! PROGRAM: nfourier3
2441 : : ! PURPOSE: fourier-transform the natural-spline interpolation
2442 : : ! of function Green(tau)
2443 : : ! calculate function Green(omega)
2444 : : ! I/O :
2445 : : ! VERSION: 2-16-92
2446 : : ! 29-Nov-95 removal of minimal bug concerning
2447 : : ! DIMENSION of rindata
2448 : : ! COMMENT: cf J. Stoer R. Bulirsch, Introduction to numerical
2449 : : ! analysis (Springer, New York, 1980)
2450 : : !=======+=========+=========+=========+=========+=========+=========+=$
2451 : : !
2452 : 0 : SUBROUTINE nfourier3(rindata,coutdata,lflag,Iwmax,L,Beta,AA,c1,c2,c3)
2453 : :
2454 : : ! use m_profiling
2455 : : ! include 'param.dat'
2456 : : ! use defs_basis
2457 : : !Arguments ------------------------------------
2458 : : integer, intent(in) :: Iwmax,L
2459 : : logical, intent(in) :: lflag
2460 : : real*8, intent(in) :: beta,AA,c1,c2,c3
2461 : : real*8, intent(in) :: rindata(L)
2462 : : complex*16, intent(out) :: coutdata(Iwmax+1)
2463 : : !Local variables ------------------------------
2464 : : integer :: i,j,k,p
2465 : 0 : real*8 :: rincopy(L+1),a(L),b(L),c(L),d(L),u(L+1), q(L+1),XM(L+1)
2466 : : complex*16 :: cdummy,explus,ex,j_dpc
2467 : : real*8 :: one,two,zero,three,six,tau,xpi,delta,om !wn,
2468 : : complex*16 :: czero
2469 : : !***********************************************
2470 : :
2471 : : ABI_UNUSED((/aa, c3/))
2472 : : ABI_UNUSED((/lflag/))
2473 : 0 : czero=cmplx(0.d0,0.d0,kind=8)
2474 : 0 : zero=0.d0
2475 : 0 : one=1.d0
2476 : 0 : two=2.d0
2477 : 0 : three=3.d0
2478 : 0 : six=6.d0
2479 : 0 : j_dpc=cmplx(0.d0,1.d0,kind=8)
2480 : 0 : xpi = ACOS(-One)
2481 : 0 : delta = beta/float(L)
2482 : : ! c2 devrait etre nul en symetrie particule trou.
2483 [ # # ]: 0 : DO i = 1,L
2484 : 0 : tau=beta/dble(L)*dble(i-1)
2485 : : !worksrincopy(i) = rindata(i)-c1/two!-c2/4.d0*(-Beta+2*tau)+c3/4.d0*(beta*tau-tau*tau)
2486 : 0 : rincopy(i) = rindata(i)-c1/two-c2/4.d0*(-Beta+2.d0*tau)!+c3/4.d0*(beta*tau-tau*tau)
2487 : : ! rincopy(i) = -c2/4.d0*(-Beta+2*tau)
2488 : : ! write(99,*) i,rindata(i)
2489 : : ! write(98,*) i,tau,rincopy(i),rindata(i),rindata(i)-c1/two,-c1/two,-c2/4.d0*(-Beta+2.d0*tau)
2490 : : ! write(97,*) i,(-Beta+two*tau),c2/4.d0,2.d0*tau,-c2/4.d0*(-Beta+2.d0*tau)
2491 : : ENDDO
2492 : : ! write(99,*)
2493 : : ! write(98,*)
2494 : : ! if(lflag) then
2495 : : ! rincopy(L+1) = AA-rindata(1)
2496 : : ! else
2497 : 0 : rincopy(L+1) = -rindata(1)
2498 : : ! endif
2499 : : !DO i = 1,L+1
2500 : : ! write(999,*) i,rincopy(i)
2501 : : !ENDDO
2502 : : !write(6,*) lflag,Iwmax,L,Beta,delta
2503 : : ! Three = Two+One
2504 : : ! six = Two*Three
2505 : :
2506 : : !c
2507 : : !c spline interpolation: the spline is given by
2508 : : !c G(tau) = a(i) + b(i) (tau-tau_i) + c(i) ( )^2 + d(i) ( )^3
2509 : : !c The following formulas are taken directly from Stoer and
2510 : : !c Bulirsch p. 102
2511 : : !c
2512 : 0 : q(1) = Zero
2513 : 0 : u(1) = Zero
2514 [ # # ]: 0 : DO k = 2,L
2515 : 0 : p = q(k-1)/Two+Two
2516 : 0 : q(k)=-One/Two/p
2517 : : ! this is equation 2.4.2.10 or Bulirsch for dn. here uk=dn
2518 : 0 : u(k)=Three/delta**2*(rincopy(k+1)+rincopy(k-1)-Two*rincopy(k))
2519 : 0 : u(k)=(u(k)-u(k-1)/Two)/p
2520 : : ENDDO
2521 : 0 : XM(L+1) = 0
2522 [ # # ]: 0 : DO k = L,1,-1
2523 : 0 : XM(k) = q(k)*XM(k+1)+u(k)
2524 : : ENDDO
2525 : : !c
2526 : : !c The following formulas are taken directly from Stoer and
2527 : : !c Bulirsch p. 98 second edition.
2528 : : !c a b c d are the spline coefficients.
2529 : : !c XM(j) is the second derivative at node j
2530 : : !c
2531 : :
2532 [ # # ]: 0 : DO j = 1, L
2533 : 0 : a(j) = rincopy(j)
2534 : 0 : c(j) = XM(j)/Two
2535 : : b(j) = (rincopy(j+1)-rincopy(j))/delta - &
2536 : 0 : & (Two*XM(j)+XM(j+1))*delta/6.
2537 : 0 : d(j) = (XM(j+1)-XM(j))/(6.*delta)
2538 : : ENDDO
2539 : :
2540 : : !c
2541 : : !c The Spline multiplied by the exponential can now be exlicitely
2542 : : !c integrated. The following formulas were obtained using
2543 : : !c MATHEMATICA
2544 : : !c
2545 [ # # ]: 0 : DO i = 0,Iwmax
2546 : 0 : om = (Two*(i)+One)*xpi/Beta
2547 : 0 : coutdata(i+1) = czero
2548 [ # # ]: 0 : DO j = 1,L
2549 : 0 : cdummy = j_dpc*om*delta*j
2550 : 0 : explus = exp(cdummy)
2551 : 0 : cdummy = j_dpc*om*delta*(j-1)
2552 : 0 : ex = exp(cdummy)
2553 : : coutdata(i+1) = coutdata(i+1) + explus*(&
2554 : : & ( -six* d(j) )/om**4 + &
2555 : : & ( Two*j_dpc*c(j) + six*delta*j_dpc*d(j) )/om**3 +&
2556 : : & ( b(j)+ Two*delta*c(j)+ three*delta**2*d(j) )/om**2 +&
2557 : : & (- j_dpc*a(j) - delta*j_dpc*b(j) - delta**2*j_dpc*c(j) -&
2558 : 0 : & delta**3*j_dpc*d(j))/om)
2559 : :
2560 : : coutdata(i+1) = coutdata(i+1) + ex*(&
2561 : : & six*d(j)/om**4 - Two*j_dpc*c(j)/om**3 &
2562 : 0 : & -b(j)/om**2 + j_dpc*a(j)/om)
2563 : : ENDDO
2564 : : !write(100,*) i,real(coutdata(i+1)),imag(coutdata(i+1))
2565 : : ENDDO
2566 : : ! DO i = 0,Iwmax
2567 : : ! wn=3.1415926/beta*dble(2*i+1)
2568 : : ! write(101,*) wn,real(coutdata(i+1)),imag(coutdata(i+1)),c1/wn,-c2/(wn*wn),-c3/(wn*wn*wn)
2569 : : ! coutdata(i+1)=coutdata(i+1)+cmplx(-c2/(wn*wn),c1/(wn))!-c3/(wn*wn*wn))
2570 : : ! !works coutdata(i+1)=coutdata(i+1)+cmplx(0.d0,c1/(wn))
2571 : : ! write(103,*) wn,real(coutdata(i+1)),imag(coutdata(i+1)),c1/wn,-c2/(wn*wn),-c3/(wn*wn*wn)
2572 : : ! ENDDO
2573 : : ! write(101,*)
2574 : : ! write(103,*)
2575 : 0 : end subroutine nfourier3
2576 : :
2577 [ # # # # : 0 : END MODULE m_GreenHyboffdiagComplex
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
2578 : : !!***
|
