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