Branch data Line data Source code
1 : :
2 : : #if defined HAVE_CONFIG_H
3 : : #include "config.h"
4 : : #endif
5 : : !!****m* ABINIT/m_GreenHyboffdiag
6 : : !! NAME
7 : : !! m_GreenHyboffdiag
8 : : !!
9 : : !! FUNCTION
10 : : !! Manage a green function for one orbital
11 : : !!
12 : : !! COPYRIGHT
13 : : !! Copyright (C) 2013-2026 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_GreenHyboffdiag
24 : :
25 : : USE m_global
26 : : USE m_MatrixHyb
27 : : USE m_Vector
28 : : USE m_VectorInt
29 : : USE m_ListCdagC
30 : : USE m_MapHyb
31 : : USE_MPI
32 : :
33 : : IMPLICIT NONE
34 : :
35 : : public :: GreenHyboffdiag_init
36 : : public :: GreenHyboffdiag_reset
37 : : public :: GreenHyboffdiag_clear
38 : : public :: GreenHyboffdiag_setOperW
39 : : public :: GreenHyboffdiag_measHybrid
40 : : public :: GreenHyboffdiag_getHybrid
41 : : public :: GreenHyboffdiag_setN
42 : : public :: GreenHyboffdiag_setMuD1
43 : : public :: GreenHyboffdiag_setMoments
44 : : public :: GreenHyboffdiag_backFourier
45 : : public :: GreenHyboffdiag_forFourier
46 : : public :: GreenHyboffdiag_print
47 : : public :: GreenHyboffdiag_destroy
48 : : public :: nfourier3
49 : :
50 : : !!***
51 : :
52 : : !!****t* m_GreenHyboffdiag/GreenHyboffdiag
53 : : !! NAME
54 : : !! GreenHyboffdiag
55 : : !!
56 : : !! FUNCTION
57 : : !! This structured datatype contains the necessary data
58 : : !!
59 : : !! COPYRIGHT
60 : : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
61 : : !! This file is distributed under the terms of the
62 : : !! GNU General Public License, see ~abinit/COPYING
63 : : !! or http://www.gnu.org/copyleft/gpl.txt .
64 : : !!
65 : : !! SOURCE
66 : :
67 : : TYPE GreenHyboffdiag
68 : :
69 : : LOGICAL :: set = .FALSE.
70 : : ! True if variable of type GreenHyboffdiag is initialized
71 : :
72 : : LOGICAL :: setT = .FALSE.
73 : : ! True if variable oper contains data
74 : :
75 : : LOGICAL :: setW = .FALSE.
76 : : ! True if variable oper_w contains data
77 : :
78 : : LOGICAL :: have_MPI = .FALSE.
79 : : ! True if MPI is used.
80 : :
81 : : INTEGER :: setMk = 0
82 : : ! setMk=0 is moments for Fourier transform are not computed
83 : :
84 : : INTEGER :: samples
85 : : ! samples=imaginary time slices (dmftqmc_l+1)
86 : :
87 : : INTEGER :: measurements
88 : : ! number of measurements for the Green's function
89 : :
90 : : INTEGER :: factor
91 : : ! if the move is not accepted, the statistic weight has to be
92 : : ! increased for the current configuration.
93 : :
94 : : INTEGER :: MY_COMM
95 : : ! MPI Communicator
96 : :
97 : : INTEGER :: size
98 : : ! size=1
99 : :
100 : : INTEGER :: rank
101 : : ! rank=0
102 : :
103 : : INTEGER :: Wmax
104 : : ! samples-1 if frequency Green's function
105 : :
106 : : INTEGER :: iTech
107 : : ! Precise if Frequency Green's function is computed or not
108 : :
109 : : INTEGER :: nflavors
110 : : ! Number of flavors
111 : :
112 : : DOUBLE PRECISION :: beta
113 : : ! Inverse of temperature
114 : :
115 : : DOUBLE PRECISION :: inv_beta
116 : : ! Temperature
117 : :
118 : : DOUBLE PRECISION :: delta_t
119 : : ! 1/inv_dt
120 : :
121 : : DOUBLE PRECISION :: inv_dt
122 : : ! (samples-1)/beta
123 : : DOUBLE PRECISION :: signvaluemeas
124 : :
125 : : DOUBLE PRECISION :: signvalueold
126 : :
127 : : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) :: oper
128 : : ! oper(samples)
129 : :
130 : : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: omega
131 : : ! omega(Wmax)
132 : :
133 : : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) :: Mk
134 : : ! Moments for FT
135 : :
136 : : COMPLEX(KIND=8) , ALLOCATABLE, DIMENSION(:,:,:) :: oper_w
137 : : ! Frequency Green's function
138 : :
139 : : COMPLEX(KIND=8) , ALLOCATABLE, DIMENSION(:) :: oper_w_old
140 : : ! Old frequency Green's function (not used)
141 : :
142 : : TYPE(Vector) :: oper_old
143 : : ! useless data
144 : :
145 : : TYPE(VectorInt) :: index_old
146 : : ! useless data
147 : :
148 : : TYPE(MapHyb), ALLOCATABLE, DIMENSION(:,:) :: map
149 : : ! value of time and Green's functions computed in GreenHyboffdiag_measHybrid
150 : : ! These values are used to fill op%oper in the same routine.
151 : :
152 : : END TYPE GreenHyboffdiag
153 : : !!***
154 : :
155 : : CONTAINS
156 : : !!***
157 : :
158 : : !!****f* ABINIT/m_GreenHyboffdiag/GreenHyboffdiag_init
159 : : !! NAME
160 : : !! GreenHyboffdiag_init
161 : : !!
162 : : !! FUNCTION
163 : : !! Initialize and allocate
164 : : !!
165 : : !! COPYRIGHT
166 : : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
167 : : !! This file is distributed under the terms of the
168 : : !! GNU General Public License, see ~abinit/COPYING
169 : : !! or http://www.gnu.org/copyleft/gpl.txt .
170 : : !!
171 : : !! INPUTS
172 : : !! op=Green
173 : : !! samples=imaginary time slices
174 : : !! beta=inverse temperature
175 : : !! iTech=SHOULD NOT BE USED => BUGGY
176 : : !! MY_COMM=mpi_communicator
177 : : !!
178 : : !! OUTPUT
179 : : !!
180 : : !! SIDE EFFECTS
181 : : !!
182 : : !! NOTES
183 : : !!
184 : : !! SOURCE
185 : :
186 : 0 : SUBROUTINE GreenHyboffdiag_init(op, samples, beta,nflavors,iTech,MY_COMM)
187 : :
188 : :
189 : : #ifdef HAVE_MPI1
190 : : include 'mpif.h'
191 : : #endif
192 : : !Arguments ------------------------------------
193 : : TYPE(GreenHyboffdiag) , INTENT(INOUT) :: op
194 : : INTEGER , INTENT(IN ) :: samples
195 : : DOUBLE PRECISION, INTENT(IN ) :: beta
196 : : INTEGER , INTENT(IN ) :: nflavors
197 : : !INTEGER , INTENT(IN ) :: Wmax
198 : : INTEGER, OPTIONAL, INTENT(IN ) :: iTech
199 : : INTEGER, OPTIONAL, INTENT(IN ) :: MY_COMM
200 : : !Local variables ------------------------------
201 : : INTEGER :: iflavor,iflavorbis,sp1
202 : : DOUBLE PRECISION :: dt
203 : : #ifdef HAVE_MPI
204 : : INTEGER :: ierr
205 : : #endif
206 : :
207 [ # # ]: 0 : IF ( PRESENT(MY_COMM)) THEN
208 : : #ifdef HAVE_MPI
209 : 0 : op%have_MPI = .TRUE.
210 : 0 : op%MY_COMM = MY_COMM
211 : 0 : CALL MPI_Comm_rank(op%MY_COMM, op%rank, ierr)
212 : 0 : CALL MPI_Comm_size(op%MY_COMM, op%size, ierr)
213 : : #else
214 : : CALL WARN("GreenHyboffdiag_init : MPI is not used ")
215 : : op%have_MPI = .FALSE.
216 : : op%MY_COMM = -1
217 : : op%rank = 0
218 : : op%size = 1
219 : : #endif
220 : : ELSE
221 : 0 : op%have_MPI = .FALSE.
222 : 0 : op%MY_COMM = -1
223 : 0 : op%rank = 0
224 : 0 : op%size = 1
225 : : END IF
226 : :
227 : 0 : sp1 = samples + 1
228 : 0 : op%samples = sp1
229 : 0 : op%measurements = 0
230 : 0 : op%nflavors = nflavors
231 : 0 : op%beta = beta
232 : 0 : op%inv_beta = 1.d0 / beta
233 : 0 : op%inv_dt = DBLE(samples) * op%inv_beta
234 : 0 : dt = 1.d0 / op%inv_dt
235 : 0 : op%delta_t = dt
236 : : !op%Wmax = Wmax
237 : 0 : op%Wmax = -1
238 [ # # ]: 0 : FREEIF(op%Mk)
239 [ # # # # : 0 : MALLOC(op%Mk,(nflavors,nflavors,3))
# # # # #
# # # # #
# # # # #
# ]
240 [ # # ]: 0 : FREEIF(op%oper)
241 [ # # # # : 0 : MALLOC(op%oper,(sp1,nflavors,nflavors))
# # # # #
# # # # #
# # # # #
# # # ]
242 : : ! If we want to measure in frequences
243 : : ! let assume we first have "samples" frequences
244 [ # # ]: 0 : IF ( PRESENT(iTech) ) THEN
245 : 0 : op%iTech = iTech
246 : 0 : SELECT CASE (op%iTech)
247 : : CASE (GREENHYB_TAU) ! omega
248 : 0 : op%iTech = GREENHYB_TAU
249 : : CASE (GREENHYB_OMEGA) ! omega
250 : 0 : op%Wmax = samples
251 [ # # ]: 0 : FREEIF(op%oper_w)
252 [ # # # # : 0 : MALLOC(op%oper_w,(1:op%Wmax,nflavors,nflavors))
# # # # #
# # # # #
# # # # #
# # # ]
253 [ # # ]: 0 : FREEIF(op%oper_w_old)
254 [ # # # # : 0 : MALLOC(op%oper_w_old,(1:op%Wmax))
# # # # #
# # # #
# ]
255 [ # # # # : 0 : op%oper_w = CMPLX(0.d0,0.d0,8)
# # ]
256 [ # # ]: 0 : op%oper_w_old = CMPLX(0.d0,0.d0,8)
257 [ # # ]: 0 : FREEIF(op%omega)
258 [ # # # # : 0 : MALLOC(op%omega,(1:op%Wmax))
# # # # #
# # # #
# ]
259 [ # # # # : 0 : op%omega = (/ ((2.d0 * DBLE(sp1) - 1.d0)*ACOS(-1.d0)*op%inv_beta, sp1=1, op%Wmax) /)
# # # # #
# # # # #
# # # #
# ]
260 : : END SELECT
261 : : ELSE
262 : 0 : op%iTech = GREENHYB_TAU
263 : : END IF
264 : : ! end if
265 : 0 : CALL Vector_init(op%oper_old,10000)
266 : 0 : CALL VectorInt_init(op%index_old,10000)
267 [ # # # # : 0 : DT_FREEIF(op%map)
# # # # #
# ]
268 [ # # # # : 0 : MALLOC(op%map,(nflavors,nflavors))
# # # # #
# # # # #
# # # # #
# # # #
# ]
269 [ # # ]: 0 : do iflavor=1,nflavors
270 [ # # ]: 0 : do iflavorbis=1,nflavors
271 : 0 : CALL MapHyb_init(op%map(iflavor,iflavorbis),10000)
272 : : enddo
273 : : enddo
274 : :
275 [ # # # # : 0 : op%oper = 0.d0
# # ]
276 : 0 : op%signvaluemeas = 0.d0
277 : 0 : op%signvalueold = 0.d0
278 : 0 : op%set = .TRUE.
279 : 0 : op%factor = 1
280 : 0 : op%setMk = 0
281 [ # # # # : 0 : op%Mk = 0.d0
# # ]
282 : 0 : END SUBROUTINE GreenHyboffdiag_init
283 : : !!***
284 : :
285 : : !!****f* ABINIT/m_GreenHyboffdiag/GreenHyboffdiag_reset
286 : : !! NAME
287 : : !! GreenHyboffdiag_reset
288 : : !!
289 : : !! FUNCTION
290 : : !! reset green function
291 : : !!
292 : : !! COPYRIGHT
293 : : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
294 : : !! This file is distributed under the terms of the
295 : : !! GNU General Public License, see ~abinit/COPYING
296 : : !! or http://www.gnu.org/copyleft/gpl.txt .
297 : : !!
298 : : !! INPUTS
299 : : !! op=Green
300 : : !!
301 : : !! OUTPUT
302 : : !!
303 : : !! SIDE EFFECTS
304 : : !!
305 : : !! NOTES
306 : : !!
307 : : !! SOURCE
308 : :
309 : 0 : SUBROUTINE GreenHyboffdiag_reset(op)
310 : :
311 : : !Arguments ------------------------------------
312 : : TYPE(GreenHyboffdiag) , INTENT(INOUT) :: op
313 : :
314 : 0 : CALL GreenHyboffdiag_clear(op)
315 : 0 : op%setMk = 0
316 [ # # # # : 0 : op%Mk = 0.d0
# # ]
317 : 0 : op%setT = .FALSE.
318 : 0 : op%setW = .FALSE.
319 : 0 : END SUBROUTINE GreenHyboffdiag_reset
320 : : !!***
321 : :
322 : :
323 : : !!****f* ABINIT/m_GreenHyboffdiag/GreenHyboffdiag_clear
324 : : !! NAME
325 : : !! GreenHyboffdiag_clear
326 : : !!
327 : : !! FUNCTION
328 : : !! clear green function
329 : : !!
330 : : !! COPYRIGHT
331 : : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
332 : : !! This file is distributed under the terms of the
333 : : !! GNU General Public License, see ~abinit/COPYING
334 : : !! or http://www.gnu.org/copyleft/gpl.txt .
335 : : !!
336 : : !! INPUTS
337 : : !! op=Green
338 : : !!
339 : : !! OUTPUT
340 : : !!
341 : : !! SIDE EFFECTS
342 : : !!
343 : : !! NOTES
344 : : !!
345 : : !! SOURCE
346 : :
347 : 0 : SUBROUTINE GreenHyboffdiag_clear(op)
348 : :
349 : : !Arguments ------------------------------------
350 : : TYPE(GreenHyboffdiag) , INTENT(INOUT) :: op
351 : : INTEGER :: iflavor,iflavorbis
352 : :
353 : : !CALL Vector_clear(op%oper_old)
354 : : !CALL VectorInt_clear(op%index_old)
355 [ # # ]: 0 : do iflavor=1,op%nflavors
356 [ # # ]: 0 : do iflavorbis=1,op%nflavors
357 : 0 : CALL MapHyb_clear(op%map(iflavor,iflavorbis))
358 : : enddo
359 : : enddo
360 : 0 : op%measurements = 0
361 [ # # ]: 0 : IF ( ALLOCATED(op%oper) ) &
362 [ # # # # : 0 : op%oper = 0.d0
# # ]
363 : 0 : op%signvaluemeas = 0.d0
364 : 0 : op%signvalueold = 1.d0
365 [ # # ]: 0 : IF ( op%iTech .EQ. GREENHYB_OMEGA ) THEN
366 [ # # ]: 0 : IF ( ALLOCATED(op%oper_w) ) &
367 [ # # # # : 0 : op%oper_w = CMPLX(0.d0,0.d0,8)
# # ]
368 [ # # ]: 0 : IF ( ALLOCATED(op%oper_w_old) ) &
369 [ # # ]: 0 : op%oper_w_old = CMPLX(0.d0,0.d0,8)
370 : : END IF
371 : 0 : op%factor = 0
372 : 0 : END SUBROUTINE GreenHyboffdiag_clear
373 : : !!***
374 : :
375 : : !!****f* ABINIT/m_GreenHyboffdiag/GreenHyboffdiag_setOperW
376 : : !! NAME
377 : : !! GreenHyboffdiag_setOperW
378 : : !!
379 : : !! FUNCTION
380 : : !! set Green function in frequencies
381 : : !!
382 : : !! COPYRIGHT
383 : : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
384 : : !! This file is distributed under the terms of the
385 : : !! GNU General Public License, see ~abinit/COPYING
386 : : !! or http://www.gnu.org/copyleft/gpl.txt .
387 : : !!
388 : : !! INPUTS
389 : : !! op=Green
390 : : !! Gomega=Input values
391 : : !!
392 : : !! OUTPUT
393 : : !!
394 : : !! SIDE EFFECTS
395 : : !!
396 : : !! NOTES
397 : : !!
398 : : !! SOURCE
399 : :
400 [ # # ]: 0 : SUBROUTINE GreenHyboffdiag_setOperW(op, Gomega)
401 : :
402 : : !Arguments ------------------------------------
403 : : TYPE(GreenHyboffdiag) , INTENT(INOUT) :: op
404 : : COMPLEX(KIND=8), DIMENSION(:,:,:), INTENT(IN ) :: Gomega
405 : : !Loval variables ------------------------------
406 : : INTEGER :: tail
407 : :
408 : 0 : tail = SIZE(Gomega,1)
409 [ # # ]: 0 : IF ( .NOT. op%set ) &
410 : 0 : CALL ERROR("GreenHyboffdiag_setOperW : Uninitialized GreenHyboffdiag structure")
411 [ # # ]: 0 : IF ( ALLOCATED(op%oper_w) ) THEN
412 [ # # # # ]: 0 : IF ( SIZE(op%oper_w) .NE. tail ) THEN
413 : 0 : FREE(op%oper_w)
414 [ # # # # : 0 : MALLOC(op%oper_w,(1:tail,op%nflavors,op%nflavors))
# # # # #
# # # # #
# # # # #
# ]
415 : : END IF
416 : : ELSE
417 [ # # # # : 0 : MALLOC(op%oper_w,(1:tail,op%nflavors,op%nflavors))
# # # # #
# # # # #
# # # # #
# # # ]
418 : : END IF
419 [ # # # # : 0 : op%oper_w(:,:,:) = Gomega(:,:,:)
# # ]
420 : 0 : op%Wmax = tail
421 : 0 : op%setW = .TRUE.
422 : 0 : END SUBROUTINE GreenHyboffdiag_setOperW
423 : : !!***
424 : :
425 : : !!****f* ABINIT/m_GreenHyboffdiag/GreenHyboffdiag_measHybrid
426 : : !! NAME
427 : : !! GreenHyboffdiag_measHybrid
428 : : !!
429 : : !! FUNCTION
430 : : !! Measure Green's function
431 : : !!
432 : : !! COPYRIGHT
433 : : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
434 : : !! This file is distributed under the terms of the
435 : : !! GNU General Public License, see ~abinit/COPYING
436 : : !! or http://www.gnu.org/copyleft/gpl.txt .
437 : : !!
438 : : !! INPUTS
439 : : !! op=Green
440 : : !! Mmatrix=M matrix for the current flavor
441 : : !! ListCdagC_1=list of all creator and annhilator operators
442 : : !! updated=should we accumulate or not
443 : : !!
444 : : !! OUTPUT
445 : : !!
446 : : !! SIDE EFFECTS
447 : : !!
448 : : !! NOTES
449 : : !!
450 : : !! SOURCE
451 : :
452 : 0 : SUBROUTINE GreenHyboffdiag_measHybrid(op, Mmatrix, ListCdagC_1, updated,signvalue,activeflavor)
453 : :
454 : : !Arguments ------------------------------------
455 : : TYPE(GreenHyboffdiag) , INTENT(INOUT) :: op
456 : : TYPE(MatrixHyb) , INTENT(IN ) :: Mmatrix
457 : : TYPE(ListCdagC) , INTENT(IN ) :: ListCdagC_1(op%nflavors)
458 : : DOUBLE PRECISION , INTENT(IN ) :: signvalue
459 : : LOGICAL , INTENT(IN ) :: updated
460 : : INTEGER, OPTIONAL , INTENT(IN ) :: activeflavor
461 : : !Local variables ------------------------------
462 : : INTEGER :: iC
463 : : INTEGER :: iCdag
464 : : INTEGER :: tail
465 : : INTEGER :: tailbis
466 : : !INTEGER :: index
467 : : INTEGER :: idx_old
468 : : INTEGER :: old_size
469 : : ! INTEGER :: omegaSamples
470 : : ! INTEGER :: iomega
471 : : INTEGER :: iflavor
472 : : INTEGER :: iflavorbis
473 : : INTEGER :: iC_m,iC_m_add
474 : : INTEGER :: iCdag_m,iCdag_m_add
475 : : INTEGER :: stail !,ii
476 : : ! DOUBLE PRECISION :: pi_invBeta
477 : : DOUBLE PRECISION :: mbeta_two
478 : : DOUBLE PRECISION :: beta
479 : : DOUBLE PRECISION :: beta_tc
480 : : DOUBLE PRECISION :: tcbeta_tc
481 : : DOUBLE PRECISION :: inv_dt
482 : : DOUBLE PRECISION :: tC,tc_phys
483 : : DOUBLE PRECISION :: tCdag
484 : : DOUBLE PRECISION :: time
485 : : DOUBLE PRECISION :: signe,signe2
486 : : DOUBLE PRECISION :: argument
487 : : INTEGER :: iflavorbegin,iflavorend,prtopt
488 : : !DOUBLE PRECISION :: taupi_invbeta
489 : : !COMPLEX(KIND=8) :: cargument
490 : : !COMPLEX(2*8) :: base_exp
491 : : !COMPLEX(2*8) :: increm_exp
492 : : !write(6,*) "measHybrid"
493 : 0 : prtopt=0
494 [ # # ]: 0 : IF ( op%set .EQV. .FALSE. ) &
495 : 0 : CALL ERROR("GreenHyboffdiag_measHybrid : green operator not set ")
496 : 0 : stail=0
497 [ # # ]: 0 : do iflavor=1,op%nflavors
498 : 0 : stail=stail + ListCdagC_1(iflavor)%tail
499 : : enddo
500 : 0 : iflavorbegin = 1
501 : 0 : iflavorend = op%nflavors
502 : :
503 : : if(present(activeflavor)) then
504 : : if(activeflavor.ne.0) then
505 : : !sui!write(6,*) "measHybrid activeflavor",activeflavor
506 : : iflavorbegin = activeflavor
507 : : iflavorend = activeflavor
508 : : endif
509 : : endif
510 : :
511 [ # # ]: 0 : IF ( stail .NE. Mmatrix%tail ) &
512 : 0 : CALL ERROR("GreenHyboffdiag_measHybrid : ListCdagC & M unconsistent ")
513 : :
514 [ # # ]: 0 : IF ( updated .EQV. .TRUE. ) THEN ! NEW change in the configuration
515 : : ! FIXME SHOULD be much more faster
516 : :
517 : :
518 : : ! write(6,*) "LKLLL2b"
519 : 0 : SELECT CASE(op%iTech)
520 : : CASE (GREENHYB_TAU)
521 : 0 : argument = DBLE(op%factor)
522 : : ! At the beginning old_size=0, then it increases
523 : : ! until
524 : : ! for all values of iC, increment green%oper with the value of the
525 : : ! Green's function in listDBLE(iC) obtained from previous iteration
526 : : ! (below)
527 : : ! ===============================================================
528 : : ! An update has been done. So the Green's function will change
529 : : ! It is thus the good moment to store the previous Green's
530 : : ! function with argument, the number of times this Green's
531 : : ! function has been constant
532 : : ! ===============================================================
533 [ # # ]: 0 : DO iflavor=1, op%nflavors
534 [ # # ]: 0 : DO iflavorbis=1, op%nflavors
535 : 0 : old_size = op%map(iflavor,iflavorbis)%tail
536 : : !write(6,*) "size listDBLE",size(op%map(iflavor,iflavorbis)%listDBLE)
537 : : !write(6,*) " measHybrid",old_size,iflavor,iflavorbis
538 [ # # ]: 0 : DO iC = 1, old_size
539 : : if(op%map(iflavor,iflavorbis)%listINT(iC)==0) then
540 : : !write(6,*) "listINT(iC)=",iC,op%map(iflavor,iflavorbis)%listINT(iC)
541 : : endif
542 : : !write(6,*) " measHybrid iflavor,iflavorbis,iC listINT ",iflavor,iflavorbis,iC,op%map(iflavor,iflavorbis)%listINT(iC)
543 : : !write(6,*) " measHybrid listDBLE ",iflavor,iflavorbis,iC,op%map(iflavor,iflavorbis)%listDBLE(iC),argument
544 : : !write(6,*) "measHybrid ic iflavor1 iflavor2 oper",ic,iflavor,iflavorbis,op%map(iflavor,iflavorbis)%listDBLE(iC) * op%signvalueold * argument
545 : : op%oper(op%map(iflavor,iflavorbis)%listINT(iC),iflavor,iflavorbis) = &
546 : : op%oper(op%map(iflavor,iflavorbis)%listINT(iC),iflavor,iflavorbis) &
547 : 0 : + op%map(iflavor,iflavorbis)%listDBLE(iC) * op%signvalueold * argument
548 : : !if(op%map(iflavor,iflavorbis)%listINT(iC)==1.and.iflavor==iflavorbis) then
549 : : ! if(iflavor==iflavorbis) then
550 : : ! !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
551 : : ! endif
552 : : ! if(iflavor==1.and.iflavorbis==6.and.op%map(iflavor,iflavorbis)%listINT(iC)==1) then
553 : : ! !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
554 : : ! endif
555 : : ! if(iflavor==6.and.iflavorbis==1.and.op%map(iflavor,iflavorbis)%listINT(iC)==1) then
556 : : ! !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
557 : : ! endif
558 : : END DO
559 : : ! tail**2 is the number of possible t-t'
560 : : ! MapHyb_setSize with resize map tail*tail will thus be the new
561 : : ! op%map%tail
562 : : ! update size of map and map%tail
563 : : CALL MapHyb_setSize(op%map(iflavor,iflavorbis),&
564 : 0 : & ListCdagC_1(iflavor)%tail*ListCdagC_1(iflavorbis)%tail)
565 : : END DO
566 : : END DO
567 : 0 : op%signvaluemeas = op%signvaluemeas + op%signvalueold * argument
568 : 0 : op%measurements = op%measurements + op%factor
569 : : !sui!write(6,*) " measurements", op%measurements
570 : : !sui! write(6,*) " signvaluemeas",op%signvaluemeas,op%signvalueold*argument
571 : : !sui! write(6,*) " signvaluemeas/measurements",op%signvaluemeas/op%measurements
572 : :
573 : : ! This is new measurement, thus op%factor should be put to one
574 : 0 : op%factor = 1
575 : : ! write(6,*) "LKLLL2C"
576 : :
577 : :
578 : : ! initialized index idx_old for the doubles loops over flavors and segments.
579 : :
580 : : ! setup usefull quantities
581 : 0 : beta = op%beta
582 : 0 : mbeta_two = -(beta*0.5d0)
583 : 0 : inv_dt = op%inv_dt
584 : :
585 : : ! WARNING time is not the time but just a temporary variable.
586 : : ! Index Time has been calculated previously and is in mat_tau
587 : :
588 : : ! initialized index for each annihilation time of a segment for a given flavor
589 : :
590 : : ! initialized index for each creation time of a segment for another flavor
591 : :
592 : 0 : iC_m=0
593 : 0 : iC_m_add=0
594 [ # # ]: 0 : DO iflavor=1,op%nflavors
595 : 0 : tail=ListCdagC_1(iflavor)%tail
596 : : !write(6,*) " measHybrid iflavor",iflavor,tail
597 : :
598 : 0 : iCdag_m=0
599 : 0 : iCdag_m_add=0
600 [ # # ]: 0 : DO iflavorbis=1,op%nflavors
601 : 0 : tailbis=ListCdagC_1(iflavorbis)%tail
602 : : !write(6,*) " measHybrid iflavorbis",iflavorbis,tailbis
603 : 0 : idx_old = 0
604 : :
605 [ # # ]: 0 : DO iC = 1, tail
606 : : ! tC is the annihilation (C_) time for segment iC and flavor iflavor
607 : : !-------------------------------------------------------------------
608 : 0 : tC = ListCdagC_1(iflavor)%list(iC,C_)
609 : :
610 : : !iC_m=iC_m+1 ! for Mmatrix%mat
611 : : ! For each flavor iflavor, iC start at \sum_{iflavor1<iflavor} tail(iflavor1)
612 : : ! It thus explains the presence of iC_m_add (same below for iCdag_m_add)
613 : : ! ---------------------------------------------------------------------------------
614 : 0 : iC_m=iC_m_add+iC
615 : 0 : beta_tc = beta - tC
616 : 0 : tcbeta_tc = tC * beta_tc
617 : :
618 : : !write(6,*) " measHybrid iC_m",iC_m
619 : : !write(6,*) " measHybrid tailbis",tailbis
620 [ # # ]: 0 : DO iCdag = 1, tailbis
621 : : !iCdag_m=iCdag_m+1
622 : 0 : iCdag_m=iCdag_m_add+iCdag
623 : : !write(6,*) " measHybrid iCdag_m",iCdag_m
624 : :
625 : : ! tCdag is the creation time for segment iCdag and flavor iflavorbis
626 : 0 : tCdag = ListCdagC_1(iflavorbis)%list(iCdag,Cdag_)
627 : :
628 : : ! --- time is equivalent to time=(tc-tcdag)*(beta-tc) and is only
629 : : ! --- useful for signe
630 : 0 : time = tcbeta_tc - tCdag*beta_tc
631 : :
632 : : !signe = SIGN(1.d0,time)
633 : : !time = time + (signe-1.d0)*mbeta_two
634 : : !signe = signe * SIGN(1.d0,beta-tC)
635 : : !signe = SIGN(1.d0,time) * SIGN(1.d0,beta-tC)
636 : 0 : tc_phys=tc
637 [ # # ]: 0 : if(tc>beta) tc_phys=tc-beta
638 : 0 : signe2=SIGN(1.d0,tc_phys-tcdag)
639 : :
640 [ # # ]: 0 : if(iflavor==iflavorbis) signe = SIGN(1.d0,time)
641 [ # # ]: 0 : if(iflavor/=iflavorbis) signe = signe2
642 : : ! signe = SIGN(1.d0,tc-tcdag)
643 : : ! --- tc>tcdag and beta>tc signe=1 ! segment in the middle or antisegment at the edge
644 : : ! ! tc-tcdag > 0
645 : : ! --- tc<tcdag and beta<tc signe=1 ! never
646 : : ! --- tc>tcdag and beta<tc signe=-1 ! segment at the edges
647 : : ! ! tc'-tcdag < 0 (with tc'=tc-beta) -> signe < 0
648 : : ! --- tc<tcdag and beta>tc signe=-1 ! antisegment in the middle
649 : : ! ! tc-tcdag < 0 (with tc'=tc-beta) -> signe < 0
650 : : ! 22/09/14:
651 : : ! ListCdagC_1 is the list of segment, so we are dealing
652 : : ! only with segment here. However all combination of Cdag
653 : : ! and C are taken, this it is possible that tc<tcdag
654 : :
655 : : ! 21/10/14: Wagt is important are the true times (between
656 : : ! 0 and beta). If tauC>tauCdag signe=+1
657 : : ! If tauC<tauCdag signe=-1
658 : : ! if(tc<tcdag.and.(iflavor==iflavorbis)) then
659 : : ! write(6,*) ListCdagC_1(iflavorbis)%tail
660 : : ! do ii=1, ListCdagC_1(iflavorbis)%tail
661 : : ! write(6,*) ListCdagC_1(iflavorbis)%list(ii,1), ListCdagC_1(iflavorbis)%list(ii,2)
662 : : ! enddo
663 : : ! write(6,*) "tc<tcdag", tc,tcdag,beta,iflavor,iflavorbis
664 : : ! stop
665 : : ! endif
666 : :
667 : : if(tc-tcdag>beta) then
668 : : ! write(6,*) " tc-tcdag > beta ", tcdag-tc,beta
669 : : endif
670 : : !if(tc>beta) then
671 : : ! write(6,*) " TC>BETA"
672 : : ! write(6,*) " iflavor,iflavorbis",iflavor,iflavorbis
673 : : ! write(6,*) " ic,icdag ",ic,icdag
674 : : ! write(6,*) " signe ",signe
675 : : ! write(6,*) " tc,tcdag ",tc,tcdag
676 : : ! write(6,*) " Mmatrix%mat ",Mmatrix%mat(iCdag_m,iC_m)
677 : : ! write(6,*) " Mmatrix%mat_tau ",Mmatrix%mat_tau(iCdag_m,iC_m)
678 : : !endif
679 : : ! Si iflavor/=iflavorbis, tc-tcdag can be negative..so in
680 : : ! this case, on should add beta to tc-tcdag with the minus
681 : : ! sign. NOT DONE HERE??
682 : :
683 : : ! ----- Compute the Green's function as the value of the matrix M for times iCdag and iC.
684 : 0 : argument = signe*Mmatrix%mat(iCdag_m,iC_m)
685 : :
686 : : !index = INT( ( time * inv_dt ) + 1.5d0 )
687 : : !IF (index .NE. Mmatrix%mat_tau(iCdag,iC)) THEN
688 : : ! WRITE(*,*) index, Mmatrix%mat_tau(iCdag,iC)
689 : : !! CALL ERROR("Plantage")
690 : : !END IF
691 : :
692 : 0 : idx_old = idx_old + 1
693 : :
694 : : ! --- define the value of listDBLE as a function of idx_old
695 : 0 : op%map(iflavor,iflavorbis)%listDBLE(idx_old) = argument
696 : : !write(6,*) " measHybrid listDBLE2 ",iflavor,iflavorbis,idx_old,argument
697 : : !op%map%listINT(idx_old) = index
698 : :
699 : : ! --- define the new corresponding value of listINT(idx_old) from mat_tau (integers)
700 : : ! --- idx_old has no meaning but listINT(idx_old) has.
701 : 0 : op%map(iflavor,iflavorbis)%listINT(idx_old) = Mmatrix%mat_tau(iCdag_m,iC_m)
702 : : !write(6,*) " measHybrid idx_old listINT ",idx_old,op%map(iflavor,iflavorbis)%listINT(idx_old)
703 : : !write(6,*) " measHybrid iCdag_m, iC_m, mat_tau",iCdag_m,iC_m,Mmatrix%mat_tau(iCdag_m,iC_m)
704 : : ! if(iflavor==1.and.iflavorbis==2.and.op%map(iflavor,iflavorbis)%listINT(idx_old)==1) then
705 : : ! !prt!if(prtopt==1) write(6,*) "---------------------------"
706 : : ! !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
707 : : ! !prt!if(prtopt==1) write(6,*) " ", tc-tcdag,tc_phys-tcdag
708 : : ! do ii=1, tail
709 : : ! !prt!if(prtopt==1) write(6,*) ii, ListCdagC_1(iflavor)%list(ii,1), ListCdagC_1(iflavor)%list(ii,2)
710 : : ! enddo
711 : : ! do ii=1, tailbis
712 : : ! !prt!if(prtopt==1) write(6,*) ii, ListCdagC_1(iflavorbis)%list(ii,1), ListCdagC_1(iflavorbis)%list(ii,2)
713 : : ! enddo
714 : : ! !prt!if(prtopt==1) write(6,*) "---------------------------"
715 : : ! endif
716 : : ! if(iflavor==1.and.iflavorbis==2.and.op%map(iflavor,iflavorbis)%listINT(idx_old)==999) then
717 : : ! !prt!if(prtopt==1) write(66,*) "---------------------------"
718 : : ! !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
719 : : ! !prt!if(prtopt==1) write(66,*) " ", tc-tcdag,tc_phys-tcdag
720 : : ! do ii=1, tail
721 : : ! !prt!if(prtopt==1) write(66,*) ii, ListCdagC_1(iflavor)%list(ii,1), ListCdagC_1(iflavor)%list(ii,2)
722 : : ! enddo
723 : : ! do ii=1, tailbis
724 : : ! !prt!if(prtopt==1) write(66,*) ii, ListCdagC_1(iflavorbis)%list(ii,1), ListCdagC_1(iflavorbis)%list(ii,2)
725 : : ! enddo
726 : : ! !prt!if(prtopt==1) write(66,*) "---------------------------"
727 : : ! endif
728 : : ! !if(iflavor==2.and.iflavorbis==1.and.op%map(iflavor,iflavorbis)%listINT(idx_old)==1) then
729 : : ! !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
730 : : !endif
731 : :
732 : :
733 : : END DO
734 : : END DO
735 : : ! do ii=1,tail*tailbis
736 : : ! !write(6,*) " measHybrid ii,op%map(iflavor,iflavorbis)%listINT(ii)", ii,op%map(iflavor,iflavorbis)%listINT(ii)
737 : : ! enddo
738 : 0 : iCdag_m_add=iCdag_m_add+tailbis
739 : : END DO ! iflavorbis
740 : 0 : iC_m_add=iC_m_add+tail
741 : : END DO ! iflavor
742 [ # # ]: 0 : op%signvalueold = signvalue
743 : : ! write(6,*) "LKLLL2D"
744 : : CASE (GREENHYB_OMEGA)
745 : : ! argument = DBLE(op%factor)
746 : : ! DO iomega = 1, omegaSamples
747 : : ! op%oper_w(iomega) = op%oper_w(iomega) + op%oper_w_old(iomega) * argument
748 : : ! END DO
749 : : ! op%measurements = op%measurements + op%factor
750 : :
751 : : ! op%factor = 1
752 : : ! beta = op%beta
753 : : ! mbeta_two = -(beta*0.5d0)
754 : : ! pi_invBeta = ACOS(-1.d0)/beta
755 : : ! omegaSamples = op%samples-1
756 : : ! DO iC = 1, tail
757 : : ! tC = ListCdagC_1%list(iC,C_)
758 : : ! DO iCdag = 1, tail
759 : : ! tCdag = ListCdagC_1%list(iCdag,Cdag_)
760 : : ! time = tC - tCdag
761 : :
762 : : ! signe = SIGN(1.d0,time)
763 : : ! time = time + (signe-1.d0)*mbeta_two
764 : : ! signe = signe * SIGN(1.d0,beta-tC)
765 : : ! argument = signe*Mmatrix%mat(iCdag,iC)
766 : :
767 : : ! DO iomega = 1, omegaSamples
768 : : ! !op%oper_w_old(iomega) = Mmatrix%mat_tau(iCdag,iC)*CMPLX(0.d0,argument)
769 : : ! op%oper_w_old(iomega) = EXP(CMPLX(0.d0,op%omega(iomega)*time))*CMPLX(0.d0,argument)
770 : : ! END DO
771 : : ! END DO
772 : : ! END DO
773 : : END SELECT
774 : : ELSE
775 : 0 : op%factor = op%factor + 1
776 : : END IF
777 : 0 : END SUBROUTINE GreenHyboffdiag_measHybrid
778 : : !!***
779 : :
780 : : !!****f* ABINIT/m_GreenHyboffdiag/GreenHyboffdiag_getHybrid
781 : : !! NAME
782 : : !! GreenHyboffdiag_getHybrid
783 : : !!
784 : : !! FUNCTION
785 : : !! reduce green function
786 : : !!
787 : : !! COPYRIGHT
788 : : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
789 : : !! This file is distributed under the terms of the
790 : : !! GNU General Public License, see ~abinit/COPYING
791 : : !! or http://www.gnu.org/copyleft/gpl.txt .
792 : : !!
793 : : !! INPUTS
794 : : !! op=Green
795 : : !!
796 : : !! OUTPUT
797 : : !!
798 : : !! SIDE EFFECTS
799 : : !!
800 : : !! NOTES
801 : : !!
802 : : !! SOURCE
803 : :
804 : 0 : SUBROUTINE GreenHyboffdiag_getHybrid(op)
805 : :
806 : : !Arguments ------------------------------------
807 : : TYPE(GreenHyboffdiag), INTENT(INOUT) :: op
808 : :
809 [ # # ]: 0 : IF ( op%set .EQV. .FALSE. ) &
810 : 0 : CALL ERROR("GreenHyboffdiag_getHybrid : green operator not set ")
811 : :
812 : 0 : SELECT CASE(op%iTech)
813 : : CASE (GREENHYB_TAU)
814 [ # # # # : 0 : op%oper = -(op%oper * op%inv_beta) / (DBLE(op%measurements) * op%delta_t)
# # ]
815 : : !sui!write(6,*) "measurements",op%measurements,op%delta_t,op%inv_beta
816 : : !sui!write(6,*) "signevaluemeas meas",op%signvaluemeas,op%measurements
817 : 0 : op%signvaluemeas = op%signvaluemeas / DBLE(op%measurements)
818 : : ! print*, "op%oper",op%oper(1,1,1)
819 : : !sui!write(6,*) "signevaluemeas/meas",op%signvaluemeas
820 : : ! print*, "signevaluemeas/meas",op%signvaluemeas
821 : 0 : op%setT = .TRUE.
822 : : CASE (GREENHYB_OMEGA)
823 [ # # # # : 0 : op%oper_w = -(op%oper_w * op%inv_beta) / (DBLE(op%measurements) * op%delta_t)
# # ]
824 : 0 : op%setW = .TRUE.
825 [ # # # ]: 0 : CALL GreenHyboffdiag_backFourier(op)
826 : : END SELECT
827 : :
828 : 0 : END SUBROUTINE GreenHyboffdiag_getHybrid
829 : : !!***
830 : :
831 : : !!****f* ABINIT/m_GreenHyboffdiag/GreenHyboffdiag_setN
832 : : !! NAME
833 : : !! GreenHyboffdiag_setN
834 : : !!
835 : : !! FUNCTION
836 : : !! impose number of electrons for this flavor
837 : : !!
838 : : !! COPYRIGHT
839 : : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
840 : : !! This file is distributed under the terms of the
841 : : !! GNU General Public License, see ~abinit/COPYING
842 : : !! or http://www.gnu.org/copyleft/gpl.txt .
843 : : !!
844 : : !! INPUTS
845 : : !! op=Green
846 : : !! N=number of electrons
847 : : !!
848 : : !! OUTPUT
849 : : !!
850 : : !! SIDE EFFECTS
851 : : !!
852 : : !! NOTES
853 : : !!
854 : : !! SOURCE
855 : :
856 : 0 : SUBROUTINE GreenHyboffdiag_setN(op,N)
857 : :
858 : : !Arguments ------------------------------------
859 : : TYPE(GreenHyboffdiag) , INTENT(INOUT) :: op
860 : : DOUBLE PRECISION , INTENT(IN ) :: N(op%nflavors)
861 : : INTEGER :: iflavor,iflavor2
862 : : !COMPLEX(KIND=8) :: tmpoper
863 : :
864 [ # # ]: 0 : IF ( op%set .EQV. .FALSE. ) &
865 : 0 : CALL ERROR("GreenHyboffdiag_setN: green op%operator not set ")
866 [ # # ]: 0 : DO iflavor=1, op%nflavors
867 : : ! write(6,*) "iflavor",-N(iflavor)*op%signvaluemeas ,2*op%oper(op%samples,iflavor,iflavor),(N(iflavor)-1.d0)*op%signvaluemeas
868 : : ! the mulplication by signvaluemeas is necessary because N is
869 : : ! exactly the number of electrons in the flavor iflavor whereas
870 : : ! op%oper is not exact, because it still has to be divided by
871 : : ! signvaluemeas after the MPIREDUCE
872 : 0 : op%oper(1,iflavor,iflavor) = (N(iflavor) - 1.d0)*op%signvaluemeas
873 : 0 : op%oper(op%samples,iflavor,iflavor) = - N(iflavor)*op%signvaluemeas
874 : : !op%oper(op%samples,iflavor,iflavor) = 2*op%oper(op%samples,iflavor,iflavor)
875 : : !op%oper(1,iflavor,iflavor) = 2*op%oper(1,iflavor,iflavor)
876 [ # # ]: 0 : DO iflavor2=1, op%nflavors
877 [ # # ]: 0 : if(iflavor/=iflavor2) then
878 : : ! UNEXPLAINED but MANDATORY to have exact results for U=0 nspinor=4 with pawspnorb=0
879 : : ! Correction: The fact 2 is necessary for edge points because the points are at the
880 : : ! edges.
881 : : ! It is of course necessary to fulfill exact results (U=0).
882 : : !tmpoper=(op%oper(op%samples,iflavor,iflavor2)-op%oper(1,iflavor,iflavor2))
883 : : !op%oper(op%samples,iflavor,iflavor2) = tmpoper
884 : : !op%oper(1,iflavor,iflavor2) = -tmpoper
885 : : !op%oper(op%samples,iflavor,iflavor2) = 2*op%oper(op%samples,iflavor,iflavor2)
886 : : !op%oper(1,iflavor,iflavor2) = 2*op%oper(1,iflavor,iflavor2)
887 : 0 : op%oper(op%samples,iflavor,iflavor2) = 2*op%oper(op%samples,iflavor,iflavor2)
888 : 0 : op%oper(1,iflavor,iflavor2) = 2*op%oper(1,iflavor,iflavor2)
889 : : endif
890 : : ENDDO
891 : : ENDDO
892 : 0 : END SUBROUTINE GreenHyboffdiag_setN
893 : : !!***
894 : :
895 : : !!****f* ABINIT/m_GreenHyboffdiag/GreenHyboffdiag_setMuD1
896 : : !! NAME
897 : : !! GreenHyboffdiag_setMuD1
898 : : !!
899 : : !! FUNCTION
900 : : !! Set first moments for G
901 : : !!
902 : : !! COPYRIGHT
903 : : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
904 : : !! This file is distributed under the terms of the
905 : : !! GNU General Public License, see ~abinit/COPYING
906 : : !! or http://www.gnu.org/copyleft/gpl.txt .
907 : : !!
908 : : !! INPUTS
909 : : !! op=Green
910 : : !! mu=energy level (irrespectige with fermi level)
911 : : !! d1=first moment of hybridization function ("K")
912 : : !!
913 : : !! OUTPUT
914 : : !! argout(sizeout)=description
915 : : !!
916 : : !! SIDE EFFECTS
917 : : !!
918 : : !! NOTES
919 : : !!
920 : : !! SOURCE
921 : :
922 : 0 : SUBROUTINE GreenHyboffdiag_setMuD1(op,iflavor,iflavor2,mu,d1)
923 : :
924 : : !Arguments ------------------------------------
925 : : !Arguments ------------------------------------
926 : : !scalars
927 : : DOUBLE PRECISION, INTENT(IN ) :: mu
928 : : DOUBLE PRECISION, INTENT(IN ) :: d1
929 : : INTEGER , INTENT(IN ) :: iflavor
930 : : INTEGER , INTENT(IN ) :: iflavor2
931 : : !type
932 : : TYPE(GreenHyboffdiag) , INTENT(INOUT) :: op
933 : : !Local variables ------------------------------------
934 : : DOUBLE PRECISION :: mu2
935 : : !*********************************************************************
936 : :
937 : : ABI_UNUSED((/d1/))
938 : :
939 : 0 : mu2=0
940 : 0 : if(iflavor==iflavor2) mu2=mu
941 : :
942 [ # # ]: 0 : if(iflavor==iflavor2) then
943 : 0 : op%Mk(iflavor,iflavor2,3) = -d1-(mu*mu)
944 : : !op%Mk(iflavor,iflavor2,3) = -(mu*mu)
945 : 0 : op%Mk(iflavor,iflavor2,2) = -mu
946 : : !sui!write(6,*) "setmud1",iflavor,iflavor2, op%Mk(iflavor,iflavor2,2), op%Mk(iflavor,iflavor2,3)
947 : : else
948 : 0 : op%Mk(iflavor,iflavor2,3) = 0.d0
949 : 0 : op%Mk(iflavor,iflavor2,2) = 0.d0
950 : : endif
951 : 0 : op%setMk = op%setMk + 1
952 : : !write(6,*) "mom1",op%Mk(iflavor,iflavor2,:)
953 : 0 : END SUBROUTINE GreenHyboffdiag_setMuD1
954 : : !!***
955 : :
956 : : !!****f* ABINIT/m_GreenHyboffdiag/GreenHyboffdiag_setMoments
957 : : !! NAME
958 : : !! GreenHyboffdiag_setMoments
959 : : !!
960 : : !! FUNCTION
961 : : !! Compute full moments
962 : : !!
963 : : !! COPYRIGHT
964 : : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
965 : : !! This file is distributed under the terms of the
966 : : !! GNU General Public License, see ~abinit/COPYING
967 : : !! or http://www.gnu.org/copyleft/gpl.txt .
968 : : !!
969 : : !! INPUTS
970 : : !! op=Greenb
971 : : !! u1_iflavor1=\sum_{iflavor2} U_{iflavor2,iflavor1} N_iflavor2
972 : : !! (useful for first moment)
973 : : !! u2=\sum_{iflavor1,iflavor2,iflavor3} U_{iflavor1,iflavor2} N_iflavor2
974 : : !!
975 : : !! OUTPUT
976 : : !!
977 : : !! SIDE EFFECTS
978 : : !!
979 : : !! NOTES
980 : : !!
981 : : !! CHI
982 : : !! Will be filled automatically by the parent script
983 : : !!
984 : : !! SOURCE
985 : :
986 : 0 : SUBROUTINE GreenHyboffdiag_setMoments(op,iflavor1,iflavor1b,u1,u2,u3)
987 : :
988 : : !Arguments ------------------------------------
989 : : TYPE(GreenHyboffdiag) , INTENT(INOUT) :: op
990 : : DOUBLE PRECISION, INTENT(IN ) :: u1
991 : : DOUBLE PRECISION, INTENT(IN ) :: u2
992 : : DOUBLE PRECISION, INTENT(IN ) :: u3
993 : : INTEGER , INTENT(IN ) :: iflavor1
994 : : INTEGER , INTENT(IN ) :: iflavor1b
995 : :
996 [ # # ]: 0 : if(iflavor1==iflavor1b) then
997 : 0 : op%Mk(iflavor1,iflavor1b,1) = -1.d0
998 : : ! c_a(3)=-d1-mu*mu-2(-mu)(\sum_{b.ne.a} Uab nb)
999 : 0 : op%Mk(iflavor1,iflavor1b,3) = op%Mk(iflavor1,iflavor1b,3) - 2.d0*(op%Mk(iflavor1,iflavor1b,2)*u1)
1000 : :
1001 : : ! c_a(2)=-mu+\sum_{b.ne.a} Uab n_b
1002 : 0 : op%Mk(iflavor1,iflavor1b,2) = op%Mk(iflavor1,iflavor1b,2) + u1
1003 : : !sui!write(6,*) "setmiments",iflavor1,iflavor1b,u1
1004 : :
1005 : : ! c_a(3)=c_a(3) + \sum Uab^2 nb + \sum Uba Uca <nbnc>
1006 : : ! ie c_a(3)=-d1+mu*mu-2mu*\sumb Uab nb + \sum Uab^2 nb + \sum Uba Uca <nbnc>
1007 : 0 : op%Mk(iflavor1,iflavor1b,3) = op%Mk(iflavor1,iflavor1b,3) - u2
1008 : : else
1009 : 0 : op%Mk(iflavor1,iflavor1b,1) = 0.d0
1010 : 0 : op%Mk(iflavor1,iflavor1b,2) = u3
1011 : 0 : op%Mk(iflavor1,iflavor1b,3) = 0.d0
1012 : : endif
1013 : : !write(6,*) "mom",iflavor1,iflavor1b, op%Mk(iflavor1,iflavor1b,:)
1014 : :
1015 : 0 : op%setMk = op%setMk + 1
1016 : :
1017 : 0 : END SUBROUTINE GreenHyboffdiag_setMoments
1018 : : !!***
1019 : :
1020 : : !!****f* ABINIT/m_GreenHyboffdiag/GreenHyboffdiag_backFourier
1021 : : !! NAME
1022 : : !! GreenHyboffdiag_backFourier
1023 : : !!
1024 : : !! FUNCTION
1025 : : !! perform back fourier transform
1026 : : !!
1027 : : !! COPYRIGHT
1028 : : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
1029 : : !! This file is distributed under the terms of the
1030 : : !! GNU General Public License, see ~abinit/COPYING
1031 : : !! or http://www.gnu.org/copyleft/gpl.txt .
1032 : : !!
1033 : : !! INPUTS
1034 : : !! op=Green
1035 : : !! dvgc=divergence parameter
1036 : : !!
1037 : : !! OUTPUT
1038 : : !!
1039 : : !! SIDE EFFECTS
1040 : : !!
1041 : : !! NOTES
1042 : : !!
1043 : : !! SOURCE
1044 : :
1045 : 0 : SUBROUTINE GreenHyboffdiag_backFourier(op,dvgc,func,hybri_limit,opt_hybri_limit)
1046 : :
1047 : : use m_fstrings, only : int2char4
1048 : :
1049 : : #ifdef HAVE_MPI1
1050 : : include 'mpif.h'
1051 : : #endif
1052 : : !Arguments ------------------------------------
1053 : : TYPE(GreenHyboffdiag) , INTENT(INOUT) :: op
1054 : : DOUBLE PRECISION, OPTIONAL, INTENT(IN ) :: dvgc
1055 : : CHARACTER(len=5) ,OPTIONAL, INTENT(IN) :: func
1056 : : COMPLEX(KIND=8), DIMENSION(op%nflavors,op%nflavors), OPTIONAL, INTENT(IN) :: hybri_limit
1057 : : INTEGER, OPTIONAL, INTENT(IN) :: opt_hybri_limit
1058 : : !Local variables ------------------------------
1059 : : INTEGER :: itau
1060 : : INTEGER :: iomega
1061 : : INTEGER :: omegaSamples
1062 : : INTEGER :: tauSamples
1063 : : INTEGER :: tauBegin
1064 : : INTEGER :: tauEnd
1065 : : INTEGER :: delta
1066 : : INTEGER :: residu
1067 : : INTEGER :: iflavor1
1068 : : INTEGER :: iflavor2,unitnb !,unitnb1
1069 : 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: counts
1070 : 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: displs
1071 : : DOUBLE PRECISION :: A,AA ! Correction factor
1072 : : COMPLEX(KIND=8) :: B,BB ! Correction factor
1073 : : COMPLEX(KIND=8) :: C !,CC ! Correction factor
1074 : : DOUBLE PRECISION :: inv_beta
1075 : : DOUBLE PRECISION :: pi_invBeta
1076 : : DOUBLE PRECISION :: two_invBeta
1077 : : DOUBLE PRECISION :: minusDt
1078 : : DOUBLE PRECISION :: minusOmegaTau
1079 : : DOUBLE PRECISION :: omegaa
1080 : : DOUBLE PRECISION :: minusTau
1081 : : DOUBLE PRECISION :: sumTerm
1082 : : DOUBLE PRECISION :: pi
1083 : : DOUBLE PRECISION :: twoPi
1084 : : DOUBLE PRECISION :: correction
1085 : 0 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: Domega
1086 : 0 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: A_omega
1087 : 0 : COMPLEX(KIND=8) , ALLOCATABLE, DIMENSION(:) :: C_omega
1088 : 0 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: opertau
1089 : : CHARACTER(len=5) :: funct
1090 : : character(len=4) :: tag_proc
1091 : : character(len=30) :: tmpfil
1092 : :
1093 : : #if defined HAVE_MPI && !defined HAVE_MPI2_INPLACE
1094 : : INTEGER :: my_count
1095 : : DOUBLE PRECISION, ALLOCATABLE , DIMENSION(:) :: opertau_buf
1096 : : #endif
1097 : :
1098 [ # # ]: 0 : IF ( op%set .EQV. .FALSE. ) &
1099 : 0 : CALL ERROR("GreenHyboffdiag_backFourier : Uninitialized GreenHyboffdiag structure")
1100 [ # # ]: 0 : IF ( op%setW .EQV. .FALSE. ) &
1101 : 0 : CALL ERROR("GreenHyboffdiag_backFourier : no G(iw)")
1102 : :
1103 : 0 : funct="hybri"
1104 [ # # ]: 0 : if(present(func)) funct=func
1105 : : !sui!write(6,*) funct
1106 : 0 : inv_beta = op%inv_beta
1107 : 0 : two_invBeta = 2.d0 * inv_beta
1108 : 0 : minusDt = - op%delta_t
1109 : 0 : omegaSamples = op%Wmax
1110 : 0 : tauSamples = op%samples-1
1111 : 0 : pi = ACOS(-1.d0)
1112 : 0 : twoPi = 2.d0 * pi
1113 : 0 : pi_invBeta = pi * inv_beta
1114 : : !sui!write(6,*) "omegaSamples",omegaSamples
1115 [ # # # # : 0 : MALLOC(Domega,(1:omegaSamples))
# # # # #
# # # ]
1116 [ # # # # : 0 : MALLOC(A_omega,(1:omegaSamples))
# # # # ]
1117 [ # # # # : 0 : MALLOC(C_omega,(1:omegaSamples))
# # # # #
# # # ]
1118 : : IF ( op%rank .EQ. 0 ) THEN
1119 : : !DO iflavor1 = 1, op%nflavors
1120 : : ! DO iflavor2 = 1, op%nflavors
1121 : : ! write(22236,*) "#",iflavor1,iflavor2
1122 : : ! do iomega=1,op%Wmax
1123 : : ! write(22236,*) (2.d0*DBLE(iomega)-1.d0) * pi_invBeta,real(op%oper_w(iomega,iflavor1,iflavor2))
1124 : : ! enddo
1125 : : ! write(22236,*)
1126 : : ! ENDDO
1127 : : !ENDDO
1128 : : ENDIF
1129 : :
1130 [ # # # # : 0 : op%oper = 0.d0
# # ]
1131 : :
1132 [ # # ]: 0 : DO iflavor1 = 1, op%nflavors
1133 [ # # ]: 0 : DO iflavor2 = 1, op%nflavors
1134 : : ! -- compute limit of function G*(i\omega_n)
1135 [ # # ]: 0 : if(funct=="hybri") then
1136 [ # # ]: 0 : IF ( PRESENT(dvgc) ) THEN
1137 : 0 : A = dvgc
1138 : : ELSE
1139 : : A = AIMAG(op%oper_w(omegaSamples,iflavor1,iflavor2))&! A = \lim_\infty Imag(G(iwn))*(wn)
1140 : 0 : *(2.d0*DBLE(omegaSamples)-1.d0) * pi_invBeta
1141 : : AA = AIMAG(op%oper_w(omegaSamples-10,iflavor1,iflavor2))&! A = \lim_\infty Imag(G(iwn))*(wn)
1142 : 0 : *(2.d0*DBLE(omegaSamples-10)-1.d0) * pi_invBeta
1143 : : B = op%oper_w(omegaSamples,iflavor1,iflavor2)&! A = \lim_\infty (G(iwn))*(iwn)
1144 : 0 : *(2.d0*DBLE(omegaSamples)-1.d0) *pi_invBeta*cmplx(0.d0,1.d0,kind=8)
1145 : : BB = op%oper_w(omegaSamples-10,iflavor1,iflavor2)&! A = \lim_\infty (G(iwn))*(iwn)
1146 : 0 : *(2.d0*DBLE(omegaSamples-10)-1.d0) *pi_invBeta*cmplx(0.d0,1.d0,kind=8)
1147 : : !sui!write(6,*) "B=",iflavor1,iflavor2,B,BB
1148 : : END IF
1149 [ # # # # ]: 0 : else if(iflavor1==iflavor2.and.funct=="green") then
1150 : : A = -1.d0
1151 [ # # # # ]: 0 : else if(iflavor1/=iflavor2.and.funct=="green") then
1152 : 0 : A = 0.d0
1153 : : endif
1154 : : !sui!write(6,*) "A=",iflavor1,iflavor2,A,AA
1155 : 0 : C=cmplx(-A,0.d0,kind=8)
1156 [ # # ]: 0 : if(present(hybri_limit)) then
1157 [ # # ]: 0 : if(present(opt_hybri_limit)) then
1158 [ # # ]: 0 : if(opt_hybri_limit==1) C= (hybri_limit(iflavor1,iflavor2))
1159 : : !sui!write(6,*) "C= ",C
1160 : : endif
1161 : : endif
1162 : :
1163 : :
1164 : :
1165 : : ! -- correction on G(tau=0) is thus
1166 : 0 : correction = -C*0.5d0
1167 : :
1168 : : ! -- built frequency mesh
1169 [ # # # # : 0 : Domega = (/ ((2.d0 * DBLE(iomega) - 1.d0)*pi_invbeta, iomega=1, omegaSamples) /)
# # # # #
# # # # #
# # ]
1170 : :
1171 : : ! -- built asymptotic function
1172 : : !if(present(hybri_limit)) then
1173 [ # # # # : 0 : C_omega = C / (Domega*cmplx(0.d0,1.d0,kind=8))
# # # # ]
1174 : : !else
1175 : : ! A_omega = A / Domega
1176 : : ! C_omega=cmplx(0.d0,A_omega,kind=8)
1177 : : !endif
1178 : : !write(6,*) "AC 1",A_omega(2),C_omega(2)
1179 : :
1180 : : IF ( op%rank .EQ. 0 ) THEN
1181 : : ! write(236,*) "#",iflavor1,iflavor2
1182 : : ! write(237,*) "#",iflavor1,iflavor2
1183 : : ! write(2236,*) "#",iflavor1,iflavor2
1184 : : ! write(2237,*) "#",iflavor1,iflavor2
1185 : : ! write(238,*) "#",iflavor1,iflavor2
1186 : : ! do iomega=1,op%Wmax
1187 : : ! write(2236,*) (2.d0*DBLE(iomega)-1.d0) * pi_invBeta,real(op%oper_w(iomega,iflavor1,iflavor2))
1188 : : ! write(2237,*) (2.d0*DBLE(iomega)-1.d0) * pi_invBeta,imag(op%oper_w(iomega,iflavor1,iflavor2))
1189 : : ! write(236,*) (2.d0*DBLE(iomega)-1.d0) * pi_invBeta,real(op%oper_w(iomega,iflavor1,iflavor2)-C_omega(iomega))
1190 : : ! write(237,*) (2.d0*DBLE(iomega)-1.d0) * pi_invBeta,imag(op%oper_w(iomega,iflavor1,iflavor2)-C_omega(iomega))
1191 : : ! write(238,*) (2.d0*DBLE(iomega)-1.d0) * pi_invBeta,C_omega(iomega)
1192 : : ! enddo
1193 : : ! write(236,*)
1194 : : ! write(237,*)
1195 : : ! write(2236,*)
1196 : : ! write(2237,*)
1197 : : ! write(238,*)
1198 : : ENDIF
1199 : : IF ( op%rank .EQ. 1 ) THEN
1200 : : ! write(22360,*) "#",iflavor1,iflavor2
1201 : : ! do iomega=1,op%Wmax
1202 : : ! write(22360,*) (2.d0*DBLE(iomega)-1.d0) * pi_invBeta,real(op%oper_w(iomega,iflavor1,iflavor2))
1203 : : ! enddo
1204 : : ! write(22360,*)
1205 : : ENDIF
1206 : : !open (unit=unitnb,file=trim(tmpfil),status='unknown',form='formatted')
1207 : : !write(unitnb,*) "#",iflavor1,iflavor2
1208 : : !do iomega=1,op%Wmax
1209 : : ! write(unitnb,*) (2.d0*DBLE(iomega)-1.d0) * pi_invBeta,real(op%oper_w(iomega,iflavor1,iflavor2))
1210 : : !enddo
1211 : : !write(unitnb,*)
1212 : : ! -- built time mesh
1213 [ # # ]: 0 : IF (op%have_MPI .EQV. .TRUE.) THEN
1214 : 0 : delta = tauSamples / op%size
1215 : 0 : residu = tauSamples - op%size*delta
1216 [ # # ]: 0 : IF ( op%rank .LT. op%size - residu ) THEN
1217 : 0 : tauBegin = 1 + op%rank*delta
1218 : 0 : tauEnd = (op%rank + 1)*delta
1219 : : ELSE
1220 : : ! tauBegin = (op%size-residu)*delta + 1 + (op%rank-op%size+residu)*(delta+1)
1221 : 0 : tauBegin = 1 + op%rank*(delta + 1) -op%size + residu
1222 : 0 : tauEnd = tauBegin + delta
1223 : : END IF
1224 [ # # # # : 0 : MALLOC(counts,(1:op%size))
# # # # #
# # # #
# ]
1225 [ # # # # : 0 : MALLOC(displs,(1:op%size))
# # # # #
# ]
1226 : : counts = (/ (delta, iTau=1, op%size-residu), &
1227 [ # # # # : 0 : (delta+1, iTau=op%size-residu+1, op%size) /)
# # # # #
# # # # #
# # # # #
# ]
1228 : 0 : displs(1)=0
1229 [ # # ]: 0 : DO iTau = 2, op%size
1230 : 0 : displs(iTau) = displs(iTau-1) + counts (iTau-1)
1231 : : END DO
1232 : : ELSE
1233 : : tauBegin = 1
1234 : : tauEnd = tauSamples
1235 : : END IF
1236 [ # # # # : 0 : MALLOC(opertau,(1:tauSamples+1))
# # # # #
# # # ]
1237 : 0 : do iomega=1,omegaSamples
1238 : : ! write(6,*) iomega, imag(op%oper_w(iomega,iflavor1,iflavor2)), A_omega(iomega) ,"#diff"
1239 : : enddo
1240 : 0 : unitnb=70000+op%rank
1241 : 0 : call int2char4(op%rank,tag_proc)
1242 : 0 : tmpfil = 'counts'//tag_proc
1243 : : ! open (unit=unitnb,file=trim(tmpfil),status='unknown',form='formatted')
1244 : : ! write(unitnb,*) "#",iflavor1,iflavor2
1245 : : ! do itau=1,op%size
1246 : : ! write(unitnb,*) itau,counts(itau),displs(itau)
1247 : : ! enddo
1248 : : ! write(unitnb,*)
1249 : :
1250 : 0 : unitnb=10000+op%rank
1251 : 0 : call int2char4(op%rank,tag_proc)
1252 : 0 : tmpfil = 'oper_w'//tag_proc
1253 : : ! open (unit=unitnb,file=trim(tmpfil),status='unknown',form='formatted')
1254 : : ! write(unitnb,*) "#",iflavor1,iflavor2,C
1255 : : ! ! C_omega et oper_w differents Domega identique. Est ce du a des
1256 : : !! ! diago differentes pour chaque procs dans qmc_prep_ctqmc
1257 : : !! do iomega=1,op%Wmax
1258 : : !! write(unitnb,*) (2.d0*DBLE(iomega)-1.d0) * pi_invBeta,real(op%oper_w(iomega,iflavor1,iflavor2)),C_omega(iomega),Domega(iomega)
1259 : : ! enddo
1260 : : ! write(unitnb,*)
1261 : :
1262 : : ! unitnb=40000+op%rank
1263 : : ! unitnb1=50000+op%rank
1264 : : ! call int2char4(op%rank,tag_proc)
1265 : : ! tmpfil = 'tauend'//tag_proc
1266 : : ! open (unit=unitnb,file=trim(tmpfil),status='unknown',form='formatted')
1267 : : ! tmpfil = 'taubegin'//tag_proc
1268 : : ! open (unit=unitnb1,file=trim(tmpfil),status='unknown',form='formatted')
1269 : : ! write(unitnb,*) "#",iflavor1,iflavor2
1270 : : ! write(unitnb1,*) "#",iflavor1,iflavor2
1271 : :
1272 : : ! -- compute Fourier transformation
1273 [ # # ]: 0 : opertau=0.d0
1274 [ # # ]: 0 : DO itau = tauBegin, tauEnd
1275 : : !DO itau = max(tauBegin-1,1), tauEnd
1276 : 0 : minusTau = DBLE(itau -1) * minusDt
1277 [ # # ]: 0 : DO iomega = 1, omegaSamples
1278 : 0 : omegaa = Domega(iomega)
1279 : 0 : minusOmegaTau = MOD(omegaa*minusTau, TwoPi)
1280 : : sumTerm = REAL(( op%oper_w(iomega,iflavor1,iflavor2) &
1281 : : - C_omega(iomega) ) &
1282 : : !- CMPLX(0.d0, A_omega(iomega),8) ) &
1283 : 0 : * EXP( CMPLX(0.d0, minusOmegaTau, 8)))
1284 : 0 : opertau(itau) = opertau(itau) + sumTerm
1285 : : ! Domega et minusomegatau identique MAIS oper_w different
1286 : : !write(unitnb,*) iomega,Domega(iomega),real(C_omega(iomega)),imag(C_omega(iomega))
1287 : : !if(itau==tauend) then
1288 : : ! write(unitnb,*) iomega, sumTerm,opertau(itau),minusOmegaTau,op%oper_w(iomega,iflavor1,iflavor2),Domega(iomega)
1289 : : !endif
1290 : : !if(itau==max(tauBegin-1,1)) then
1291 : : ! write(unitnb1,*) iomega, sumTerm,opertau(itau),minusOmegaTau,op%oper_w(iomega,iflavor1,iflavor2),Domega(iomega)
1292 : : !endif
1293 : :
1294 : : END DO
1295 : : ! if(itau==tauEnd) write(unitnb,*)
1296 : : ! if(itau==max(tauBegin-1,1)) write(unitnb1,*)
1297 : : ! if(iflavor1==iflavor2) then
1298 : 0 : opertau(itau) = correction + two_invBeta*opertau(itau)
1299 : : ! if(itau==tauend) then
1300 : : ! write(unitnb,*) "final", opertau(itau),correction
1301 : : ! endif
1302 : : ! if(itau==max(tauBegin-1,1)) then
1303 : : ! write(unitnb1,*) "final",opertau(itau),correction
1304 : : ! endif
1305 : : !write(66666,*) itau, opertau(itau),correction
1306 : : ! else
1307 : : ! opertau(itau) = &
1308 : : ! two_invBeta*opertau(itau)
1309 : : ! endif
1310 : : END DO
1311 : : !unitnb=20000+op%rank
1312 : : !call int2char4(op%rank,tag_proc)
1313 : : !tmpfil = 'opertau'//tag_proc
1314 : : !open (unit=unitnb,file=trim(tmpfil),status='unknown',form='formatted')
1315 : : !write(unitnb,*) "#",iflavor1,iflavor2,tauBegin,tauEnd
1316 : : !!do itau=tauBegin, tauEnd
1317 : : !do itau=1,tauSamples
1318 : : ! write(unitnb,*) itau,opertau(itau)
1319 : : !enddo
1320 : : !write(unitnb,*)
1321 : : !opertau(tauBegin-1)=0.d0
1322 : : !opertau(tauEnd+1)=0.d0
1323 : : !write(66666,*)
1324 : :
1325 : : ! -- Gather
1326 [ # # ]: 0 : IF ( op%have_MPI .EQV. .TRUE. ) THEN
1327 : : ! rassembler les resultats
1328 : : #ifdef HAVE_MPI
1329 : : #if defined HAVE_MPI2_INPLACE
1330 : : CALL MPI_ALLGATHERV(MPI_IN_PLACE, 0, MPI_DOUBLE_PRECISION, &
1331 : : opertau, counts, displs, &
1332 : 0 : MPI_DOUBLE_PRECISION, op%MY_COMM, residu)
1333 : : #else
1334 : : my_count=tauBegin-tauEnd+1
1335 : : MALLOC(opertau_buf,(my_count))
1336 : : opertau_buf(1:my_count)=opertau(tauBegin:tauEnd)
1337 : : CALL MPI_ALLGATHERV(opertau_buf, my_count, MPI_DOUBLE_PRECISION, &
1338 : : opertau, counts, displs, &
1339 : : MPI_DOUBLE_PRECISION, op%MY_COMM, residu)
1340 : : FREE(opertau_buf)
1341 : : #endif
1342 : : #endif
1343 [ # # ]: 0 : FREE(counts)
1344 [ # # ]: 0 : FREE(displs)
1345 : : END IF
1346 : : ! unitnb=30000+op%rank
1347 : : ! call int2char4(op%rank,tag_proc)
1348 : : ! tmpfil = 'opertau_MPI_'//tag_proc
1349 : : ! open (unit=unitnb,file=trim(tmpfil),status='unknown',form='formatted')
1350 : : ! write(unitnb,*) "#",iflavor1,iflavor2
1351 : : ! do itau=tauBegin, tauEnd
1352 : : ! write(unitnb,*) itau,opertau(itau)
1353 : : ! enddo
1354 : : ! write(unitnb,*)
1355 : :
1356 : : ! -- Add correction for discontinuity.
1357 : : ! if(iflavor1==iflavor2) then
1358 : : !G(0+)-G(0-)=G(0+)+G(beta-)=A
1359 : 0 : opertau(tauSamples+1) = -real(C) - opertau(1)
1360 : : !sui!write(6,*) "BackFourier",opertau(tauSamples+1),opertau(1),real(C)
1361 : :
1362 : 0 : op%setT = .TRUE.
1363 : : ! endif
1364 [ # # ]: 0 : op%oper(:,iflavor1,iflavor2)=opertau(:)
1365 : 0 : FREE(opertau)
1366 : : END DO ! iflavor2
1367 : : END DO ! iflavor1
1368 : : ! -- End loop over flavors.
1369 : :
1370 [ # # ]: 0 : FREE(Domega)
1371 [ # # ]: 0 : FREE(A_omega)
1372 [ # # ]: 0 : FREE(C_omega)
1373 : 0 : close(236)
1374 : 0 : close(237)
1375 : :
1376 [ # # # # : 0 : END SUBROUTINE GreenHyboffdiag_backFourier
# # ]
1377 : : !!***
1378 : :
1379 : : !!****f* ABINIT/m_GreenHyboffdiag/GreenHyboffdiag_forFourier
1380 : : !! NAME
1381 : : !! GreenHyboffdiag_forFourier
1382 : : !!
1383 : : !! FUNCTION
1384 : : !! perform forward fourier transform
1385 : : !!
1386 : : !! COPYRIGHT
1387 : : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
1388 : : !! This file is distributed under the terms of the
1389 : : !! GNU General Public License, see ~abinit/COPYING
1390 : : !! or http://www.gnu.org/copyleft/gpl.txt .
1391 : : !!
1392 : : !! INPUTS
1393 : : !! op=Green
1394 : : !! Wmax=linear maximum frequency
1395 : : !!
1396 : : !! OUTPUT
1397 : : !! Gomega=Results for omega frequencies
1398 : : !! omega=ask frequencies
1399 : : !!
1400 : : !! SIDE EFFECTS
1401 : : !!
1402 : : !! NOTES
1403 : : !!
1404 : : !! SOURCE
1405 : :
1406 [ # # # # : 0 : SUBROUTINE GreenHyboffdiag_forFourier(op, Gomega, omega, Wmax)
# # # # #
# # # ]
1407 : : !Arguments ------------------------------------
1408 : :
1409 : : #ifdef HAVE_MPI1
1410 : : include 'mpif.h'
1411 : : #endif
1412 : : TYPE(GreenHyboffdiag) , INTENT(INOUT) :: op
1413 : : COMPLEX(KIND=8), DIMENSION(:,:,:), OPTIONAL, INTENT(INOUT) :: Gomega ! INOUT for MPI
1414 : : COMPLEX(KIND=8), DIMENSION(:), OPTIONAL, INTENT(IN ) :: omega
1415 : : INTEGER , OPTIONAL, INTENT(IN ) :: Wmax
1416 : : INTEGER :: i
1417 : : INTEGER :: j
1418 : : INTEGER :: iflavor1
1419 : : INTEGER :: iflavor2
1420 : : INTEGER :: nflavors
1421 : : INTEGER :: L
1422 : : INTEGER :: Lspline
1423 : : INTEGER :: Nom
1424 : : INTEGER :: omegaBegin
1425 : : INTEGER :: omegaEnd
1426 : : INTEGER :: deltaw
1427 : : INTEGER :: residu
1428 : 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: counts
1429 : 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: displs
1430 : : DOUBLE PRECISION :: beta
1431 : : DOUBLE PRECISION :: tau
1432 : : DOUBLE PRECISION :: delta
1433 : : DOUBLE PRECISION :: deltabis
1434 : : DOUBLE PRECISION :: inv_delta
1435 : : DOUBLE PRECISION :: inv_delta2
1436 : : DOUBLE PRECISION :: omdeltabis
1437 : : DOUBLE PRECISION :: tmp
1438 : : DOUBLE PRECISION :: xpi
1439 : 0 : DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: diag
1440 : 0 : DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: diagL
1441 : 0 : DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: lastR
1442 : 0 : DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: lastC
1443 : 0 : DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: XM
1444 : 0 : DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: X2
1445 : : DOUBLE PRECISION :: iw
1446 : : COMPLEX(KIND=8) :: iwtau
1447 : 0 : COMPLEX(KIND=8), ALLOCATABLE, DIMENSION(:) :: Gwtmp
1448 : 0 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: omegatmp
1449 : : #if defined HAVE_MPI && !defined HAVE_MPI2_INPLACE
1450 : : INTEGER :: my_count
1451 : : COMPLEX(KIND=8), ALLOCATABLE , DIMENSION(:) :: Gwtmp_buf
1452 : : #endif
1453 : :
1454 : 0 : nflavors=op%nflavors
1455 : :
1456 : : !sui!write(6,*) " Fourier transformation begin"
1457 : :
1458 [ # # ]: 0 : IF ( op%set .EQV. .FALSE. ) &
1459 : 0 : CALL ERROR("GreenHyboffdiag_forFourier : Uninitialized GreenHyboffdiag structure")
1460 [ # # ]: 0 : IF ( op%setT .EQV. .FALSE. ) &
1461 : 0 : CALL ERROR("GreenHyboffdiag_forFourier : no G(tau)")
1462 : : !write(6,*) "op%setMk=", op%setMk
1463 [ # # ]: 0 : IF ( op%setMk .NE. 2*nflavors*nflavors ) &
1464 : 0 : CALL WARNALL("GreenHyboffdiag_forFourier : green does not have moments ")
1465 : :
1466 : 0 : L = op%samples
1467 : :
1468 : 0 : xpi=acos(-1.d0) !!! XPI=PI
1469 : 0 : beta = op%beta
1470 : 0 : Nom = op%Wmax
1471 [ # # ]: 0 : IF ( PRESENT(Gomega) ) THEN
1472 : 0 : Nom = SIZE(Gomega,1)
1473 : : !IF ( op%rank .EQ. 0 ) &
1474 : : !!write(6,*) "size Gomega", Nom
1475 : : END IF
1476 [ # # ]: 0 : IF ( PRESENT(omega) ) THEN
1477 [ # # # # ]: 0 : IF ( PRESENT(Gomega) .AND. SIZE(omega) .NE. Nom ) THEN
1478 : 0 : CALL ERROR("GreenHyboffdiag_forFourier : sizes mismatch ")
1479 : : !ELSE
1480 : : !Nom = SIZE(omega)
1481 : : END IF
1482 : : END IF
1483 [ # # ]: 0 : IF ( .NOT. PRESENT(Gomega) .AND. .NOT. PRESENT(omega) ) THEN
1484 [ # # ]: 0 : IF ( PRESENT(Wmax) ) THEN
1485 : 0 : Nom=Wmax
1486 : : ELSE
1487 : 0 : CALL ERROR("GreenHyboffdiag_forFourier : Missing argument Wmax")
1488 : : END IF
1489 : : END IF
1490 : :
1491 : : !!IF ( ALLOCATED(op%oper_w) ) THEN
1492 : : !! IF ( SIZE(op%oper_w,1) .NE. Nom ) THEN
1493 : : !! FREE(op%oper_w)
1494 : : !! MALLOC(op%oper_w,(1:Nom,nflavors,nflavors))
1495 : : !! END IF
1496 : : !!ELSE
1497 : : !! MALLOC(op%oper_w,(1:Nom,nflavors,nflavors))
1498 : : !!END IF
1499 : :
1500 : : !!write(6,*) "PRESENT(GOMEGA)", PRESENT(GOMEGA)
1501 : : !!write(6,*) "PRESENT(OMEGA)", PRESENT(OMEGA)
1502 : : !call flush(6)
1503 : :
1504 : 0 : delta=op%delta_t
1505 : 0 : inv_delta = op%inv_dt
1506 : 0 : inv_delta2 = inv_delta*inv_delta
1507 : :
1508 [ # # # # : 0 : MALLOC(diagL,(L-1))
# # # # #
# # # #
# ]
1509 [ # # # # : 0 : MALLOC(lastR,(L-1))
# # # # ]
1510 [ # # # # : 0 : MALLOC(diag,(L))
# # # # #
# # # ]
1511 [ # # # # : 0 : MALLOC(lastC,(L-1))
# # # # ]
1512 : :
1513 : : !(cf Stoer) for the spline interpolation :
1514 : : ! second derivatives XM solution of A*XM=B.
1515 : : !A=(2.4.2.11) of Stoer&Bulirsch + 2 limit conditions
1516 : : !The LU decomposition of A is known explicitly;
1517 : :
1518 : 0 : diag (1) = 4.d0 ! 1.d0 *4.d0 factor 4 added for conditionning
1519 : 0 : diagL(1) = 0.25d0 !1.d0/4.d0
1520 : 0 : lastR(1) = -0.5d0 ! -2.d0/4.d0
1521 : 0 : lastC(1) = 4.d0 ! 1.d0*4.d0
1522 : 0 : diag (2) = 4.d0
1523 : 0 : diagL(2) = 0.25d0
1524 : 0 : lastR(2) = -0.25d0
1525 : 0 : lastC(2) = -1.d0
1526 : :
1527 [ # # ]: 0 : DO i = 3, L-2
1528 : 0 : tmp = 4.d0 - diagL(i-1)
1529 : 0 : diagL(i) = 1.d0 / tmp
1530 : : END DO
1531 [ # # ]: 0 : DO i = 3, L-2
1532 : 0 : diag (i) = 1.d0 / diagL(i)
1533 : 0 : lastR(i) = -(lastR(i-1)*diagL(i))
1534 : 0 : lastC(i) = -(lastC(i-1)*diagL(i-1))
1535 : : END DO
1536 : :
1537 : 0 : tmp = 1.d0/diag(L-2)
1538 : 0 : diag (L-1) = 4.d0 - tmp
1539 : 0 : lastR(L-1) = (1.d0 - lastR(L-2))/ diag(L-1)
1540 : : !diagL(L-1) = lastR(L-1)
1541 : 0 : diagL(L-1) = 0.d0 ! for the Lq=B resolution
1542 : : !lastC(L-1) = 1.d0 - lastC(L-2)*diagL(L-1) ! equivalent to the next line
1543 : 0 : lastC(L-1) = 1.d0 - (lastC(L-2)*lastR(L-1)) ! True value
1544 : 0 : diag (L ) = 2.d0! - DOT_PRODUCT( lastR , lastC )
1545 : 0 : tmp = 0.d0
1546 [ # # ]: 0 : DO i = 1, L-1
1547 : 0 : tmp = tmp + lastR(i)*lastC(i)
1548 : : END DO
1549 : 0 : diag (L ) = diag (L ) - tmp
1550 : 0 : lastC(L-1) = lastC(L-1)-1.d0 ! 1 is removed for the u.XM=q resolution
1551 : :
1552 [ # # # # : 0 : MALLOC(XM,(L))
# # # # ]
1553 [ # # # # : 0 : MALLOC(Gwtmp,(1:Nom))
# # # # #
# # # ]
1554 : :
1555 : 0 : Lspline = L-1
1556 [ # # # # : 0 : MALLOC(X2,(1:Lspline+1)) ! We impose L = Nom
# # # # #
# # # ]
1557 : :
1558 [ # # ]: 0 : IF ( op%have_MPI .EQV. .TRUE. ) THEN
1559 : 0 : deltaw = Nom / op%size
1560 : 0 : residu = Nom - op%size*deltaw
1561 [ # # ]: 0 : IF ( op%rank .LT. op%size - residu ) THEN
1562 : 0 : omegaBegin = 1 + op%rank*deltaw
1563 : 0 : omegaEnd = (op%rank + 1)*deltaw
1564 : : ELSE
1565 : : ! tauBegin = (op%size-residu)*deltaw + 1 + (op%rank-op%size+residu)*(deltaw+1)
1566 : 0 : omegaBegin = 1 + op%rank*(deltaw + 1) -op%size + residu
1567 : 0 : omegaEnd = omegaBegin + deltaw
1568 : : END IF
1569 [ # # # # : 0 : MALLOC(counts,(1:op%size))
# # # # #
# # # ]
1570 [ # # # # : 0 : MALLOC(displs,(1:op%size))
# # # # ]
1571 : : counts = (/ (deltaw, i=1, op%size-residu), &
1572 [ # # # # : 0 : (deltaw+1, i=op%size-residu+1, op%size) /)
# # # # #
# # # # #
# # # # #
# ]
1573 : 0 : displs(1)=0
1574 [ # # ]: 0 : DO i = 2, op%size
1575 : 0 : displs(i) = displs(i-1) + counts (i-1)
1576 : : END DO
1577 : : ELSE
1578 : : omegaBegin = 1
1579 : : omegaEnd = Nom
1580 : : END IF
1581 : :
1582 : : ! op%Mk(iflavor1,iflavor2,1) = 0.d0
1583 : : ! DO iflavor1 = 1, nflavors
1584 : : ! op%Mk(iflavor1,iflavor1,1) = -1.d0
1585 : : ! ENDDO
1586 : : ! op%Mk(:,:,3) = 0.d0
1587 : :
1588 [ # # # # : 0 : MALLOC(omegatmp,(omegaBegin:omegaEnd))
# # # # #
# # # ]
1589 [ # # ]: 0 : IF ( PRESENT(omega) ) THEN
1590 [ # # # # : 0 : omegatmp(omegaBegin:omegaEnd) = (/ (AIMAG(omega(i)),i=omegaBegin,omegaEnd) /)
# # ]
1591 : : ELSE
1592 [ # # # # : 0 : omegatmp(omegaBegin:omegaEnd) = (/ ((((2.d0*DBLE(i)-1.d0)*xpi)/Beta), i=omegaBegin,omegaEnd) /)
# # ]
1593 : : END IF
1594 : :
1595 [ # # ]: 0 : DO iflavor1 = 1, nflavors
1596 [ # # ]: 0 : DO iflavor2 = 1, nflavors
1597 : : ! write(6,*) " Moments:",op%Mk(iflavor1,iflavor2,:),iflavor1,iflavor2
1598 : :
1599 : : ! construct the B vector from A.Xm=B
1600 : 0 : XM(1) = 4.d0*op%Mk(iflavor1,iflavor2,3)
1601 : : XM(L) = (6.d0 * inv_delta) * ( op%Mk(iflavor1,iflavor2,2) - ( &
1602 : : (op%oper(2,iflavor1,iflavor2)-op%oper(1,iflavor1,iflavor2)) + &
1603 : 0 : (op%oper(L,iflavor1,iflavor2)-op%oper(L-1,iflavor1,iflavor2)) ) * inv_delta )
1604 : : ! built generic second derivative of oper
1605 : : !sui!write(6,*) "XM 1 L",XM(1),XM(L),op%Mk(iflavor1,iflavor2,2),op%Mk(iflavor1,iflavor2,3)
1606 [ # # ]: 0 : DO i = 2, L-1
1607 : : XM(i) = (6.d0 * inv_delta2) * ( (op%oper(i+1,iflavor1,iflavor2) &
1608 : : - 2.d0 * op%oper(i,iflavor1,iflavor2)) &
1609 : 0 : + op%oper(i-1,iflavor1,iflavor2) )
1610 : : !sui!write(6,*) "XM",i,XM(i),op%oper(i,iflavor1,iflavor2)
1611 : : END DO
1612 : :
1613 : : ! Find second derivatives XM: Solve the system
1614 : : ! SOLVING Lq= XM
1615 : : ! q = XM
1616 [ # # ]: 0 : do j=1,L-1
1617 : 0 : XM(j+1)=XM(j+1)-(diagL(j)*XM(j))
1618 : 0 : XM(L) =XM(L) -(lastR(j)*XM(j))
1619 : : end do
1620 : :
1621 : :
1622 : : ! SOLVING U.XM=q
1623 : : ! XM = q
1624 [ # # ]: 0 : do j=L-1,2,-1
1625 : 0 : XM(j+1) = XM(j+1) / diag(j+1)
1626 : 0 : XM(j)= (XM(j)-(XM(L)*lastC(j)))-XM(j+1)
1627 : : end do
1628 : 0 : XM(2) = XM(2) / diag(2)
1629 : 0 : XM(1) = (XM(1)-XM(L)*lastC(1)) / diag(1)
1630 : :
1631 : :
1632 : :
1633 : : !Construct L2 second derivative from known derivatives XM
1634 : 0 : deltabis = beta / DBLE(Lspline)
1635 [ # # ]: 0 : DO i = 1, Lspline
1636 : 0 : tau = deltabis * DBLE(i-1)
1637 : 0 : j = ((L-1)*(i-1))/Lspline + 1!INT(tau * inv_delta) + 1
1638 : 0 : X2(i) = inv_delta * ( XM(j)*(DBLE(j)*delta - tau ) + XM(j+1)*(tau - DBLE(j-1)*delta) )
1639 : : END DO
1640 : 0 : X2(Lspline+1) = XM(L)
1641 : :
1642 : :
1643 [ # # ]: 0 : DO i = omegaBegin, omegaEnd
1644 : 0 : iw = omegatmp(i)
1645 : 0 : omdeltabis = iw*deltabis
1646 : 0 : Gwtmp(i)=CMPLX(0.d0,0.d0,8)
1647 [ # # ]: 0 : DO j=2, Lspline ! We impose L+1 = Nom
1648 : 0 : iwtau = CMPLX(0.d0,omdeltabis*DBLE(j-1),8)
1649 : 0 : Gwtmp(i) = Gwtmp(i) + EXP(iwtau) * CMPLX((X2(j+1) + X2(j-1))-2.d0*X2(j),0.d0,8)
1650 : : !write(6,*) "ww",i,j,Gwtmp(i),X2(j),iwtau
1651 : : END DO
1652 : : Gwtmp(i) = Gwtmp(i)/CMPLX(((iw*iw)*(iw*iw)*deltabis),0.d0,8) &
1653 : : + CMPLX( ( ((X2(2)-X2(1))+(X2(Lspline+1)-X2(Lspline)))/((iw*iw)*deltabis) -op%Mk(iflavor1,iflavor2,2) ) &
1654 : 0 : /(iw*iw) , (op%Mk(iflavor1,iflavor2,1)-op%Mk(iflavor1,iflavor2,3)/(iw*iw))/iw , 8)
1655 : : !+ CMPLX( (X2(2)-X2(1))+(X2(Lspline+1)-X2(Lspline)), 0.d0, 8 ) ) &
1656 : : ! / (((iw*iw)*(iw*iw))*CMPLX(deltabis,0.d0,8)) &
1657 : : !- CMPLX(op%Mk(1),0.d0,8)/iw &
1658 : : !+ CMPLX(op%Mk(2),0.d0,8)/(iw*iw) &
1659 : : !- CMPLX(op%Mk(3),0.d0,8)/((iw*iw)*iw)
1660 : : !IF ( op%rank .EQ. 0 ) write(12819,*) iw,gwtmp(i)
1661 : : END DO
1662 : : !call flush(12819)
1663 [ # # ]: 0 : IF ( op%have_MPI .EQV. .TRUE. ) THEN
1664 : : #ifdef HAVE_MPI
1665 : : #if defined HAVE_MPI2_INPLACE
1666 : : CALL MPI_ALLGATHERV(MPI_IN_PLACE, 0, MPI_DOUBLE_COMPLEX, &
1667 : : Gwtmp , counts, displs, &
1668 : 0 : MPI_DOUBLE_COMPLEX, op%MY_COMM, residu)
1669 : : #else
1670 : : my_count=omegaBegin-omegaEnd+1
1671 : : MALLOC(Gwtmp_buf,(my_count))
1672 : : Gwtmp_buf(1:my_count)=Gwtmp(omegaBegin:omegaEnd)
1673 : : CALL MPI_ALLGATHERV(Gwtmp_buf, my_count, MPI_DOUBLE_COMPLEX, &
1674 : : Gwtmp , counts, displs, &
1675 : : MPI_DOUBLE_COMPLEX, op%MY_COMM, residu)
1676 : : FREE(Gwtmp_buf)
1677 : : #endif
1678 : : #endif
1679 : : END IF
1680 [ # # ]: 0 : IF ( PRESENT(Gomega) ) THEN
1681 [ # # ]: 0 : Gomega(:,iflavor1,iflavor2) = Gwtmp(:)
1682 : : END IF
1683 : 0 : op%setW = .TRUE.
1684 : : ENDDO ! iflavor1
1685 : : ENDDO ! iflavor2
1686 : : !!op%oper_w=Gomega
1687 : : do iflavor1=1,nflavors
1688 : : !sui!write(6,*) iflavor1
1689 : : do i=1,Nom
1690 : : !write(6,*) "w",i,op%oper_w(i,iflavor1,iflavor1)
1691 : : enddo
1692 : : enddo
1693 : :
1694 [ # # ]: 0 : FREE(Gwtmp)
1695 : 0 : FREE(diagL)
1696 : 0 : FREE(lastR)
1697 : 0 : FREE(diag)
1698 : 0 : FREE(lastC)
1699 [ # # ]: 0 : FREE(XM)
1700 [ # # ]: 0 : FREE(omegatmp)
1701 [ # # ]: 0 : FREE(X2)
1702 [ # # ]: 0 : FREE(counts)
1703 [ # # ]: 0 : FREE(displs)
1704 : :
1705 : 0 : END SUBROUTINE GreenHyboffdiag_forFourier
1706 : : !!***
1707 : :
1708 : : !!****f* ABINIT/m_GreenHyboffdiag/GreenHyboffdiag_print
1709 : : !! NAME
1710 : : !! GreenHyboffdiag_print
1711 : : !!
1712 : : !! FUNCTION
1713 : : !! print Green function
1714 : : !!
1715 : : !! COPYRIGHT
1716 : : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
1717 : : !! This file is distributed under the terms of the
1718 : : !! GNU General Public License, see ~abinit/COPYING
1719 : : !! or http://www.gnu.org/copyleft/gpl.txt .
1720 : : !!
1721 : : !! INPUTS
1722 : : !! op=Green
1723 : : !! ostream=file stream
1724 : : !!
1725 : : !! OUTPUT
1726 : : !!
1727 : : !! SIDE EFFECTS
1728 : : !!
1729 : : !! NOTES
1730 : : !!
1731 : : !! SOURCE
1732 : :
1733 : 0 : SUBROUTINE GreenHyboffdiag_print(op, ostream)
1734 : :
1735 : : !Arguments ------------------------------------
1736 : : TYPE(GreenHyboffdiag), INTENT(IN) :: op
1737 : : INTEGER, OPTIONAL , INTENT(IN) :: ostream
1738 : : !Local variables ------------------------------
1739 : : INTEGER :: ostream_val
1740 : : INTEGER :: isample
1741 : : INTEGER :: samples
1742 : : INTEGER :: iflavor1
1743 : : INTEGER :: iflavor2
1744 : :
1745 : :
1746 [ # # ]: 0 : IF ( op%set .EQV. .FALSE. ) &
1747 : 0 : CALL ERROR("GreenHyboffdiag_print : green op%operator not set ")
1748 : :
1749 [ # # ]: 0 : IF ( PRESENT(ostream) ) THEN
1750 : 0 : ostream_val = ostream
1751 : : ELSE
1752 : 0 : ostream_val = 66
1753 : 0 : OPEN(UNIT=ostream_val,FILE="Green.dat")
1754 : : END IF
1755 : :
1756 : 0 : samples = op%samples
1757 : :
1758 [ # # ]: 0 : DO iflavor1=1,op%nflavors
1759 [ # # ]: 0 : DO iflavor2=1,op%nflavors
1760 : 0 : WRITE(ostream_val,'(a,i3,a,i3,a)') "## (iflavor1,iflavor2)= (", iflavor1,",",iflavor2,")"
1761 [ # # ]: 0 : DO isample = 1, samples
1762 : 0 : WRITE(ostream_val,*) DBLE(isample-1)*op%delta_t, op%oper(isample,iflavor1,iflavor2)
1763 : : END DO
1764 : 0 : WRITE(ostream_val,*)
1765 : : END DO
1766 : : END DO
1767 : :
1768 [ # # ]: 0 : IF ( .NOT. PRESENT(ostream) ) &
1769 : 0 : CLOSE(ostream_val)
1770 : 0 : END SUBROUTINE GreenHyboffdiag_print
1771 : : !!***
1772 : :
1773 : : !!****f* ABINIT/m_GreenHyboffdiag/GreenHyboffdiag_destroy
1774 : : !! NAME
1775 : : !! GreenHyboffdiag_destroy
1776 : : !!
1777 : : !! FUNCTION
1778 : : !! destroy green function
1779 : : !!
1780 : : !! COPYRIGHT
1781 : : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
1782 : : !! This file is distributed under the terms of the
1783 : : !! GNU General Public License, see ~abinit/COPYING
1784 : : !! or http://www.gnu.org/copyleft/gpl.txt .
1785 : : !!
1786 : : !! INPUTS
1787 : : !! op=Green
1788 : : !!
1789 : : !! OUTPUT
1790 : : !!
1791 : : !! SIDE EFFECTS
1792 : : !!
1793 : : !! NOTES
1794 : : !!
1795 : : !! SOURCE
1796 : :
1797 : 0 : SUBROUTINE GreenHyboffdiag_destroy(op)
1798 : :
1799 : : !Arguments ------------------------------------
1800 : : TYPE(GreenHyboffdiag), INTENT(INOUT) :: op
1801 : : INTEGER :: iflavor,iflavorbis
1802 : :
1803 : 0 : op%set = .FALSE.
1804 : 0 : op%setT = .FALSE.
1805 : 0 : op%setW = .FALSE.
1806 : 0 : op%samples = 0
1807 : 0 : op%measurements = 0
1808 : 0 : op%beta = 0.d0
1809 : 0 : op%inv_beta = 0.d0
1810 : 0 : op%inv_dt = 0.d0
1811 : 0 : op%delta_t = 0.d0
1812 : 0 : CALL VectorInt_destroy(op%index_old)
1813 : 0 : CALL Vector_destroy(op%oper_old)
1814 [ # # ]: 0 : do iflavor=1,op%nflavors
1815 [ # # ]: 0 : do iflavorbis=1,op%nflavors
1816 : : !sui!write(6,*) "test",iflavor,iflavorbis
1817 : 0 : CALL MapHyb_destroy(op%map(iflavor,iflavorbis))
1818 : : enddo
1819 : : enddo
1820 [ # # # # : 0 : DT_FREEIF(op%map)
# # # # #
# ]
1821 [ # # ]: 0 : FREEIF(op%oper)
1822 [ # # ]: 0 : FREEIF(op%Mk)
1823 [ # # ]: 0 : FREEIF(op%oper_w)
1824 [ # # ]: 0 : FREEIF(op%oper_w_old)
1825 [ # # ]: 0 : FREEIF(op%omega)
1826 : 0 : END SUBROUTINE GreenHyboffdiag_destroy
1827 : :
1828 : : !!***
1829 : : ! This routine contains direct and inverse fourier transformation
1830 : : ! It is a modification of a routine of the GNU GPL
1831 : : ! code available on http://dmft.rutgers.edu/ and
1832 : : ! described in the RMP 2006 paper written by
1833 : : ! G.Kotliar, S.Y.Savrasov, K.Haule, V.S.Oudovenko, O.Parcollet, C.A.Marianetti
1834 : : !=======+=========+=========+=========+=========+=========+=========+=$
1835 : : ! TYPE : SUBROUTINE
1836 : : ! PROGRAM: nfourier3
1837 : : ! PURPOSE: fourier-transform the natural-spline interpolation
1838 : : ! of function Green(tau)
1839 : : ! calculate function Green(omega)
1840 : : ! I/O :
1841 : : ! VERSION: 2-16-92
1842 : : ! 29-Nov-95 removal of minimal bug concerning
1843 : : ! DIMENSION of rindata
1844 : : ! COMMENT: cf J. Stoer R. Bulirsch, Introduction to numerical
1845 : : ! analysis (Springer, New York, 1980)
1846 : : !=======+=========+=========+=========+=========+=========+=========+=$
1847 : : !
1848 : 0 : SUBROUTINE nfourier3(rindata,coutdata,lflag,Iwmax,L,Beta,AA,c1,c2,c3)
1849 : :
1850 : : ! use m_profiling
1851 : : ! include 'param.dat'
1852 : : ! use defs_basis
1853 : : !Arguments ------------------------------------
1854 : : integer, intent(in) :: Iwmax,L
1855 : : logical, intent(in) :: lflag
1856 : : real*8, intent(in) :: beta,AA,c1,c2,c3
1857 : : real*8, intent(in) :: rindata(L)
1858 : : complex*16, intent(out) :: coutdata(Iwmax+1)
1859 : : !Local variables ------------------------------
1860 : : integer :: i,j,k,p
1861 : 0 : real*8 :: rincopy(L+1),a(L),b(L),c(L),d(L),u(L+1), q(L+1),XM(L+1)
1862 : : complex*16 :: cdummy,explus,ex,j_dpc
1863 : : real*8 :: one,two,zero,three,six,tau,xpi,delta,om !wn,
1864 : : complex*16 :: czero
1865 : : !***********************************************
1866 : :
1867 : : ABI_UNUSED((/aa, c3/))
1868 : : ABI_UNUSED((/lflag/))
1869 : 0 : czero=cmplx(0.d0,0.d0)
1870 : 0 : zero=0.d0
1871 : 0 : one=1.d0
1872 : 0 : two=2.d0
1873 : 0 : three=3.d0
1874 : 0 : six=6.d0
1875 : 0 : j_dpc=dcmplx(0.d0,1.d0)
1876 : 0 : xpi = ACOS(-One)
1877 : 0 : delta = beta/float(L)
1878 : : ! c2 devrait etre nul en symetrie particule trou.
1879 [ # # ]: 0 : DO i = 1,L
1880 : 0 : tau=beta/dble(L)*dble(i-1)
1881 : : !worksrincopy(i) = rindata(i)-c1/two!-c2/4.d0*(-Beta+2*tau)+c3/4.d0*(beta*tau-tau*tau)
1882 : 0 : rincopy(i) = rindata(i)-c1/two-c2/4.d0*(-Beta+2.d0*tau)!+c3/4.d0*(beta*tau-tau*tau)
1883 : : ! rincopy(i) = -c2/4.d0*(-Beta+2*tau)
1884 : : ! write(99,*) i,rindata(i)
1885 : : ! write(98,*) i,tau,rincopy(i),rindata(i),rindata(i)-c1/two,-c1/two,-c2/4.d0*(-Beta+2.d0*tau)
1886 : : ! write(97,*) i,(-Beta+two*tau),c2/4.d0,2.d0*tau,-c2/4.d0*(-Beta+2.d0*tau)
1887 : : ENDDO
1888 : : ! write(99,*)
1889 : : ! write(98,*)
1890 : : ! if(lflag) then
1891 : : ! rincopy(L+1) = AA-rindata(1)
1892 : : ! else
1893 : 0 : rincopy(L+1) = -rindata(1)
1894 : : ! endif
1895 : : !DO i = 1,L+1
1896 : : ! write(999,*) i,rincopy(i)
1897 : : !ENDDO
1898 : : !write(6,*) lflag,Iwmax,L,Beta,delta
1899 : : ! Three = Two+One
1900 : : ! six = Two*Three
1901 : :
1902 : : !c
1903 : : !c spline interpolation: the spline is given by
1904 : : !c G(tau) = a(i) + b(i) (tau-tau_i) + c(i) ( )^2 + d(i) ( )^3
1905 : : !c The following formulas are taken directly from Stoer and
1906 : : !c Bulirsch p. 102
1907 : : !c
1908 : 0 : q(1) = Zero
1909 : 0 : u(1) = Zero
1910 [ # # ]: 0 : DO k = 2,L
1911 : 0 : p = q(k-1)/Two+Two
1912 : 0 : q(k)=-One/Two/p
1913 : : ! this is equation 2.4.2.10 or Bulirsch for dn. here uk=dn
1914 : 0 : u(k)=Three/delta**2*(rincopy(k+1)+rincopy(k-1)-Two*rincopy(k))
1915 : 0 : u(k)=(u(k)-u(k-1)/Two)/p
1916 : : ENDDO
1917 : 0 : XM(L+1) = 0
1918 [ # # ]: 0 : DO k = L,1,-1
1919 : 0 : XM(k) = q(k)*XM(k+1)+u(k)
1920 : : ENDDO
1921 : : !c
1922 : : !c The following formulas are taken directly from Stoer and
1923 : : !c Bulirsch p. 98 second edition.
1924 : : !c a b c d are the spline coefficients.
1925 : : !c XM(j) is the second derivative at node j
1926 : : !c
1927 : :
1928 [ # # ]: 0 : DO j = 1, L
1929 : 0 : a(j) = rincopy(j)
1930 : 0 : c(j) = XM(j)/Two
1931 : : b(j) = (rincopy(j+1)-rincopy(j))/delta - &
1932 : 0 : & (Two*XM(j)+XM(j+1))*delta/6.
1933 : 0 : d(j) = (XM(j+1)-XM(j))/(6.*delta)
1934 : : ENDDO
1935 : :
1936 : : !c
1937 : : !c The Spline multiplied by the exponential can now be exlicitely
1938 : : !c integrated. The following formulas were obtained using
1939 : : !c MATHEMATICA
1940 : : !c
1941 [ # # ]: 0 : DO i = 0,Iwmax
1942 : 0 : om = (Two*(i)+One)*xpi/Beta
1943 : 0 : coutdata(i+1) = czero
1944 [ # # ]: 0 : DO j = 1,L
1945 : 0 : cdummy = j_dpc*om*delta*j
1946 : 0 : explus = exp(cdummy)
1947 : 0 : cdummy = j_dpc*om*delta*(j-1)
1948 : 0 : ex = exp(cdummy)
1949 : : coutdata(i+1) = coutdata(i+1) + explus*(&
1950 : : & ( -six* d(j) )/om**4 + &
1951 : : & ( Two*j_dpc*c(j) + six*delta*j_dpc*d(j) )/om**3 +&
1952 : : & ( b(j)+ Two*delta*c(j)+ three*delta**2*d(j) )/om**2 +&
1953 : : & (- j_dpc*a(j) - delta*j_dpc*b(j) - delta**2*j_dpc*c(j) -&
1954 : 0 : & delta**3*j_dpc*d(j))/om)
1955 : :
1956 : : coutdata(i+1) = coutdata(i+1) + ex*(&
1957 : : & six*d(j)/om**4 - Two*j_dpc*c(j)/om**3 &
1958 : 0 : & -b(j)/om**2 + j_dpc*a(j)/om)
1959 : : ENDDO
1960 : : !write(100,*) i,real(coutdata(i+1)),imag(coutdata(i+1))
1961 : : ENDDO
1962 : : ! DO i = 0,Iwmax
1963 : : ! wn=3.1415926/beta*dble(2*i+1)
1964 : : ! write(101,*) wn,real(coutdata(i+1)),imag(coutdata(i+1)),c1/wn,-c2/(wn*wn),-c3/(wn*wn*wn)
1965 : : ! coutdata(i+1)=coutdata(i+1)+cmplx(-c2/(wn*wn),c1/(wn))!-c3/(wn*wn*wn))
1966 : : ! !works coutdata(i+1)=coutdata(i+1)+cmplx(0.d0,c1/(wn))
1967 : : ! write(103,*) wn,real(coutdata(i+1)),imag(coutdata(i+1)),c1/wn,-c2/(wn*wn),-c3/(wn*wn*wn)
1968 : : ! ENDDO
1969 : : ! write(101,*)
1970 : : ! write(103,*)
1971 : 0 : end subroutine nfourier3
1972 : :
1973 [ # # # # : 0 : END MODULE m_GreenHyboffdiag
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
1974 : : !!***
|
