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