Line data Source code
1 :
2 : #if defined HAVE_CONFIG_H
3 : #include "config.h"
4 : #endif
5 : !!****m* ABINIT/m_CtqmcoffdiagComplex
6 : !! NAME
7 : !! m_CtqmcoffdiagComplex
8 : !!
9 : !! FUNCTION
10 : !! Manage and drive all the CTQMC
11 : !! Should not be used if you don't know what you do
12 : !! Please use CtqmcoffdiagComplexComplexInterface
13 : !!
14 : !! COPYRIGHT
15 : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder, B. Amadon, J. Denier)
16 : !! This file is distributed under the terms of the
17 : !! GNU General Public License, see ~abinit/COPYING
18 : !! or http://www.gnu.org/copyleft/gpl.txt .
19 : !!
20 : !! NOTES
21 : !!
22 : !! SOURCE
23 :
24 : #include "defs.h"
25 : #define CTQMC_SLICE1 100
26 : ! Coupe Sweeps en 100
27 : #define CTQMC_SLICE2 100
28 : ! Coupe modNoise1 en 2000
29 : #define CTQMC_SEGME 1
30 : #define CTQMC_ANTIS -2
31 : #define CTQMC_ADDED 3
32 : #define CTQMC_REMOV 4
33 : #define CTQMC_DETSI 5
34 : MODULE m_CtqmcoffdiagComplex
35 :
36 : USE m_Global
37 : USE m_GreenHyboffdiagComplex
38 : USE m_BathOperatoroffdiagComplex
39 : USE m_ImpurityOperator
40 : USE m_Stat
41 : USE m_FFTHyb
42 : USE m_OurRng
43 : use defs_basis
44 : use m_abicore
45 : use m_io_tools, only : open_file
46 :
47 : #ifdef HAVE_MPI2
48 : USE mpi
49 : #endif
50 : IMPLICIT NONE
51 :
52 : public :: CtqmcoffdiagComplex_init
53 : public :: CtqmcoffdiagComplex_setParameters
54 : public :: CtqmcoffdiagComplex_setSweeps
55 : public :: CtqmcoffdiagComplex_setSeed
56 : public :: CtqmcoffdiagComplex_allocateAll
57 : public :: CtqmcoffdiagComplex_allocateOpt
58 : public :: CtqmcoffdiagComplex_setG0wTab
59 : public :: CtqmcoffdiagComplex_setU
60 : public :: CtqmcoffdiagComplex_clear
61 : public :: CtqmcoffdiagComplex_reset
62 : public :: CtqmcoffdiagComplex_setMu
63 : public :: CtqmcoffdiagComplex_computeF
64 : public :: CtqmcoffdiagComplex_run
65 : public :: CtqmcoffdiagComplex_tryAddRemove
66 : public :: CtqmcoffdiagComplex_trySwap
67 : public :: CtqmcoffdiagComplex_measN
68 : public :: CtqmcoffdiagComplex_measCorrelation
69 : public :: CtqmcoffdiagComplex_measPerturbation
70 : public :: CtqmcoffdiagComplex_getResult
71 : public :: CtqmcoffdiagComplex_symmetrizeGreen
72 : public :: CtqmcoffdiagComplex_getGreen
73 : public :: CtqmcoffdiagComplex_getD
74 : public :: CtqmcoffdiagComplex_getE
75 : public :: CtqmcoffdiagComplex_printAll
76 : public :: CtqmcoffdiagComplex_printQMC
77 : public :: CtqmcoffdiagComplex_printGreen
78 : public :: CtqmcoffdiagComplex_printD
79 : public :: CtqmcoffdiagComplex_printE
80 : public :: CtqmcoffdiagComplex_printPerturbation
81 : public :: CtqmcoffdiagComplex_printCorrelation
82 : public :: CtqmcoffdiagComplex_printSpectra
83 : public :: CtqmcoffdiagComplex_destroy
84 : public :: CtqmcoffdiagComplex_setMagmom
85 :
86 : !!***
87 :
88 : !!****t* m_CtqmcoffdiagComplex/CtqmcoffdiagComplex
89 : !! NAME
90 : !! CtqmcoffdiagComplex
91 : !!
92 : !! FUNCTION
93 : !! This structured datatype contains the necessary data
94 : !!
95 : !! COPYRIGHT
96 : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
97 : !! This file is distributed under the terms of the
98 : !! GNU General Public License, see ~abinit/COPYING
99 : !! or http://www.gnu.org/copyleft/gpl.txt .
100 : !!
101 : !! SOURCE
102 :
103 : TYPE CtqmcoffdiagComplex
104 :
105 : LOGICAL :: init = .FALSE.
106 : ! Flag: is MC initialized
107 :
108 : LOGICAL :: set = .FALSE.
109 : ! Flag: ??
110 :
111 : LOGICAL :: setU = .FALSE.
112 : ! Flag: is U Set ?
113 :
114 : LOGICAL :: inF = .FALSE.
115 : ! Flag: is hybridization fct in input ?
116 :
117 : LOGICAL :: done = .FALSE.
118 : ! Flag: is MC terminated ?
119 :
120 : LOGICAL :: para = .FALSE.
121 : ! Flag: do we have parameters in input
122 :
123 : LOGICAL :: have_MPI = .FALSE.
124 : ! Flag:
125 :
126 : INTEGER :: opt_movie = 0
127 : !
128 :
129 : INTEGER :: opt_analysis = 0
130 : ! correlations
131 :
132 : INTEGER :: opt_check = 0
133 : ! various check 0
134 : ! various check 1 impurity
135 : ! various check 2 bath
136 : ! various check 3 both
137 :
138 : INTEGER :: opt_order = 0
139 : ! nb of segments max for analysis
140 :
141 : INTEGER :: opt_histo = 0
142 : ! Enable histo calc.
143 :
144 : INTEGER :: opt_noise = 0
145 : ! compute noise
146 :
147 : INTEGER :: opt_spectra = 0
148 : ! markov chain FT (correlation time)
149 :
150 : INTEGER :: opt_levels = 0
151 : ! do we have energy levels
152 :
153 : INTEGER :: opt_hybri_limit = 0
154 : ! do we have limit of hybridization (yes=1)
155 :
156 : INTEGER :: opt_nondiag = 0
157 : ! if opt_nondiag = 1 F is non diagonal.
158 :
159 : INTEGER :: prtopt = 1
160 : ! printing
161 :
162 : INTEGER :: flavors
163 : ! number of flavors
164 :
165 : INTEGER :: nspinor
166 : ! number of spinor
167 :
168 : INTEGER :: measurements
169 : ! The modulo used to measure the interaction energy and the number of electrons. Example : 2 means the measure is perform every two sweeps.
170 :
171 : INTEGER :: samples
172 : ! nb of L points (dmftqmc_l)
173 :
174 : ! INTEGER :: endDensity
175 : ! second dimension of density allocatable array
176 :
177 : INTEGER(8) :: seed
178 : !
179 :
180 : INTEGER :: sweeps
181 : !
182 :
183 : INTEGER :: thermalization
184 : !
185 :
186 : INTEGER :: ostream
187 : ! output file
188 :
189 : INTEGER :: istream
190 : ! input file
191 :
192 : INTEGER :: modNoise1
193 : ! measure the noise each modNoise1
194 :
195 : INTEGER :: modNoise2
196 : ! measure the noise each modNoise2
197 :
198 : INTEGER :: activeFlavor
199 : ! orbital on which one do sth now
200 :
201 : INTEGER, DIMENSION(1:2) :: modGlobalMove
202 : ! 1: global move each modglobalmove(1)
203 : ! 2: we have done modglobalmove(2) for two different orbitals.
204 :
205 : INTEGER :: Wmax
206 : ! Max freq for FT
207 :
208 : DOUBLE PRECISION, DIMENSION(1:6) :: stats
209 : ! to now how many negative determinant, antisegments,seeme.e.twfs...j
210 :
211 : DOUBLE PRECISION :: swap
212 : ! nb of successfull GM
213 :
214 : DOUBLE PRECISION :: signvalue
215 :
216 : COMPLEX(KIND=8) :: phasevalue
217 :
218 : INTEGER :: MY_COMM
219 : !
220 :
221 : INTEGER :: rank
222 : !
223 :
224 : INTEGER :: size
225 : ! size of MY_COMM
226 :
227 : DOUBLE PRECISION :: runTime ! time for the run routine
228 : !
229 :
230 : DOUBLE PRECISION :: beta
231 : !
232 :
233 : DOUBLE PRECISION :: U
234 :
235 : COMPLEX(KIND=8), ALLOCATABLE, DIMENSION(:) :: mu
236 : ! levels
237 :
238 : COMPLEX(KIND=8), ALLOCATABLE, DIMENSION(:,:) :: hybri_limit
239 : ! coeff A such that F=-A/(iwn)
240 :
241 : TYPE(GreenHyboffdiagComplex) :: Greens
242 : ! Green's function
243 :
244 : ! DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: measN
245 : COMPLEX(KIND=8), ALLOCATABLE, DIMENSION(:,:) :: measN
246 : ! measure of occupations (3or4,flavor)
247 :
248 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,: ) :: measDE
249 : ! (flavor,flavor) double occupancies
250 : ! (1,1): total energy of correlation.
251 :
252 : DOUBLE PRECISION :: a_Noise
253 : ! Noise a exp (-bx) for the noise
254 :
255 : DOUBLE PRECISION :: b_Noise
256 : ! Noise a exp (-bx) for the noise
257 :
258 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) :: abNoiseG !(ab,tau,flavor)
259 : ! Noise but for G
260 :
261 : TYPE(Vector) , DIMENSION(1:2) :: measNoise
262 : TYPE(Vector), ALLOCATABLE, DIMENSION(:,:,:) :: measNoiseG !(tau,flavor,mod)
263 : ! accumulate each value relataed to measurenoise 1 2
264 :
265 : !#ifdef CTCtqmcoffdiagComplex_ANALYSIS
266 : ! INTEGER :: order
267 : DOUBLE PRECISION :: inv_dt
268 : ! 1/(beta/L)
269 :
270 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,: ) :: measPerturbation
271 : ! opt_order,nflavor
272 :
273 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: occup_histo_time
274 : ! nflavor
275 :
276 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: occupconfig
277 : ! 2**nflavor
278 :
279 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: suscep
280 : ! samples
281 :
282 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: chi
283 : ! samples
284 :
285 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: chicharge
286 : ! samples
287 :
288 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: ntot
289 : ! occupation total, t2g, eg
290 :
291 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,: ) :: meas_fullemptylines
292 : ! opt_order,nflavor
293 :
294 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) :: measCorrelation
295 : ! segment,antisegment,nflavor,nflavor
296 :
297 : !#endif
298 : !#ifdef CTCtqmcoffdiagComplex_CHECK
299 : DOUBLE PRECISION :: errorImpurity
300 : ! check
301 :
302 : DOUBLE PRECISION :: errorBath
303 : ! for check
304 :
305 : !#endif
306 : TYPE(BathOperatoroffdiagComplex) :: Bath
307 :
308 :
309 : TYPE(ImpurityOperator) :: Impurity
310 :
311 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: density
312 :
313 : END TYPE CtqmcoffdiagComplex
314 : !!***
315 :
316 : !INTERFACE CtqmcoffdiagComplex_setG0w
317 : ! MODULE PROCEDURE CtqmcoffdiagComplex_setG0wFile, CtqmcoffdiagComplex_setG0wTab
318 : !END INTERFACE
319 :
320 : CONTAINS
321 : !!***
322 :
323 : !!****f* ABINIT/m_CtqmcoffdiagComplex/CtqmcoffdiagComplex_init
324 : !! NAME
325 : !! CtqmcoffdiagComplex_init
326 : !!
327 : !! FUNCTION
328 : !! Initialize the type CtqmcoffdiagComplex
329 : !! Allocate all the non optional variables
330 : !!
331 : !! COPYRIGHT
332 : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
333 : !! This file is distributed under the terms of the
334 : !! GNU General Public License, see ~abinit/COPYING
335 : !! or http://www.gnu.org/copyleft/gpl.txt .
336 : !!
337 : !! INPUTS
338 : !! op=ctqmc
339 : !! ostream=where to write
340 : !! istream=where to read the input parameters if so
341 : !! bFile=logical argument True if input is read from istream
342 : !! MY_COMM=mpi communicator for the CTQMC
343 : !! iBuffer=input parameters if bFile is false
344 : !!
345 : !! OUTPUT
346 : !!
347 : !! SIDE EFFECTS
348 : !!
349 : !! NOTES
350 : !!
351 : !! SOURCE
352 :
353 0 : SUBROUTINE CtqmcoffdiagComplex_init(op, ostream, istream, bFile, MY_COMM, iBuffer)
354 :
355 :
356 : #ifdef HAVE_MPI1
357 : include 'mpif.h'
358 : #endif
359 : !Arguments ------------------------------------
360 : TYPE(CtqmcoffdiagComplex), INTENT(INOUT) :: op
361 : INTEGER , INTENT(IN ) :: ostream
362 : INTEGER , INTENT(IN ) :: istream
363 : LOGICAL , INTENT(IN ) :: bFile
364 : DOUBLE PRECISION, DIMENSION(1:11), OPTIONAL, INTENT(IN) :: iBuffer
365 : INTEGER , OPTIONAL, INTENT(IN ) :: MY_COMM
366 : !Local variables ------------------------------
367 : #ifdef HAVE_MPI
368 : INTEGER :: ierr
369 : #endif
370 : !INTEGER :: iflavor
371 : #ifdef __GFORTRAN__
372 : ! INTEGER :: pid
373 : ! CHARACTER(LEN=5) :: Cpid
374 : !
375 : #endif
376 : DOUBLE PRECISION, DIMENSION(1:11) :: buffer
377 :
378 0 : op%ostream = ostream
379 0 : op%istream = istream
380 :
381 : ! --- RENICE ---
382 : !#ifdef __GFORTRAN__
383 : ! pid = GetPid()
384 : ! WRITE(Cpid,'(I5)') pid
385 : ! CALL SYSTEM('renice +19 '//TRIM(ADJUSTL(Cpid))//' > /dev/null')
386 : !#endif
387 : !! --- RENICE ---
388 :
389 0 : IF ( PRESENT(MY_COMM)) THEN
390 : #ifdef HAVE_MPI
391 0 : op%have_MPI = .TRUE.
392 0 : op%MY_COMM = MY_COMM
393 0 : CALL MPI_Comm_rank(op%MY_COMM, op%rank, ierr)
394 0 : CALL MPI_Comm_size(op%MY_COMM, op%size, ierr)
395 : #else
396 : CALL WARN("MPI is not used ")
397 : op%have_MPI = .FALSE.
398 : op%MY_COMM = -1
399 : op%rank = 0
400 : op%size = 1
401 : #endif
402 : ELSE
403 0 : op%have_MPI = .FALSE.
404 0 : op%MY_COMM = -1
405 0 : op%rank = 0
406 0 : op%size = 1
407 : END IF
408 :
409 : !IF ( op%rank .EQ. 0 ) THEN
410 : ! WRITE(ostream,'(A20)') 'Job reniced with +19'
411 : !CALL FLUSH(ostream)
412 : !END IF
413 :
414 0 : IF ( bFile .EQV. .TRUE. ) THEN
415 0 : IF ( op%rank .EQ. 0 ) THEN
416 :
417 0 : READ(istream,*) buffer(1) !iseed
418 0 : READ(istream,*) buffer(2) !op%sweeps
419 0 : READ(istream,*) buffer(3) !op%thermalization
420 0 : READ(istream,*) buffer(4) !op%measurements
421 0 : READ(istream,*) buffer(5) !op%flavors
422 0 : READ(istream,*) buffer(6) !op%samples
423 0 : READ(istream,*) buffer(7) !op%beta
424 0 : READ(istream,*) buffer(8) !U
425 0 : READ(istream,*) buffer(9) !iTech
426 0 : READ(istream,*) buffer(11)!op%nspinor
427 : !READ(istream,*) buffer(9) !Wmax
428 : !#ifdef CTCtqmcoffdiagComplex_ANALYSIS
429 : !READ(istream,*) buffer(10) !order
430 : !#endif
431 : END IF
432 :
433 : #ifdef HAVE_MPI
434 0 : IF ( op%have_MPI .EQV. .TRUE. ) &
435 : CALL MPI_Bcast(buffer, 11, MPI_DOUBLE_PRECISION, 0, &
436 0 : op%MY_COMM, ierr)
437 : #endif
438 0 : ELSE IF ( PRESENT(iBuffer) ) THEN
439 0 : buffer(1:11) = iBuffer(1:11)
440 : ELSE
441 0 : CALL ERROR("CtqmcoffdiagComplex_init : No input parameters ")
442 : END IF
443 :
444 0 : CALL CtqmcoffdiagComplex_setParameters(op, buffer)
445 :
446 0 : CALL CtqmcoffdiagComplex_allocateAll(op)
447 :
448 : CALL GreenHyboffdiagComplex_init(op%Greens,op%samples, op%beta,INT(buffer(5)), &
449 0 : iTech=INT(buffer(9)),MY_COMM=op%MY_COMM)
450 :
451 :
452 : ! op%seg_added = 0.d0
453 : ! op%anti_added = 0.d0
454 : ! op%seg_removed = 0.d0
455 : ! op%anti_removed = 0.d0
456 : ! op%seg_sign = 0.d0
457 : ! op%anti_sign = 0.d0
458 0 : op%stats(:) = 0.d0
459 : ! write(std_out,*) "op%stats",op%stats
460 0 : op%signvalue = 1.d0
461 0 : op%phasevalue = cmplx(1.d0,0.d0,kind=8)
462 : ! op%signvaluecurrent = 0.d0
463 : ! op%signvaluemeas = 0.d0
464 0 : op%swap = 0.d0
465 0 : op%runTime = 0.d0
466 :
467 0 : CALL Vector_init(op%measNoise(1),op%sweeps/op%modNoise1)
468 0 : CALL Vector_init(op%measNoise(2),(op%sweeps/op%modNoise1+1)*CTQMC_SLICE2)
469 : !CALL Vector_init(op%measNoise(3),101)
470 : !CALL Vector_init(op%measNoise(4),101)
471 :
472 0 : op%set = op%para .AND. op%inF
473 0 : op%done = .FALSE.
474 0 : op%init = .TRUE.
475 :
476 : !#ifdef CTCtqmcoffdiagComplex_CHECK
477 0 : op%errorImpurity = 0.d0
478 0 : op%errorBath = 0.d0
479 : !#endif
480 0 : END SUBROUTINE CtqmcoffdiagComplex_init
481 : !!***
482 :
483 : !!****f* ABINIT/m_CtqmcoffdiagComplex/CtqmcoffdiagComplex_setParameters
484 : !! NAME
485 : !! CtqmcoffdiagComplex_setParameters
486 : !!
487 : !! FUNCTION
488 : !! set all parameters and operators
489 : !!
490 : !! COPYRIGHT
491 : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
492 : !! This file is distributed under the terms of the
493 : !! GNU General Public License, see ~abinit/COPYING
494 : !! or http://www.gnu.org/copyleft/gpl.txt .
495 : !!
496 : !! INPUTS
497 : !! op=ctqmc
498 : !! buffer=input parameters
499 : !!
500 : !! OUTPUT
501 : !!
502 : !! SIDE EFFECTS
503 : !!
504 : !! NOTES
505 : !!
506 : !! SOURCE
507 :
508 0 : SUBROUTINE CtqmcoffdiagComplex_setParameters(op,buffer)
509 :
510 : !Arguments ------------------------------------
511 : TYPE(CtqmcoffdiagComplex), INTENT(INOUT) :: op
512 : DOUBLE PRECISION, DIMENSION(1:11), INTENT(IN ) :: buffer
513 :
514 :
515 0 : op%thermalization = INT(buffer(3)) !op%thermalization
516 0 : CALL CtqmcoffdiagComplex_setSeed(op,INT(buffer(1)))
517 0 : CALL CtqmcoffdiagComplex_setSweeps(op,buffer(2))
518 :
519 0 : op%measurements = INT(buffer(4)) !op%measurements
520 0 : op%flavors = INT(buffer(5))
521 0 : op%samples = INT(buffer(6)) !op%samples
522 0 : op%beta = buffer(7) !op%beta
523 0 : op%U = buffer(8) !U
524 0 : op%opt_nondiag = INT(buffer(10))
525 0 : op%nspinor = INT(buffer(11)) !op%nspinor
526 : ! op%mu = buffer(9) !op%mu
527 : !op%Wmax = INT(buffer(9)) !Freq
528 : !#ifdef CTCtqmcoffdiagComplex_ANALYSIS
529 : ! op%order = INT(buffer(10)) ! order
530 0 : op%inv_dt = op%samples / op%beta
531 : !#endif
532 :
533 : !CALL ImpurityOperator_init(op%Impurity,op%flavors,op%beta, op%samples)
534 0 : CALL ImpurityOperator_init(op%Impurity,op%flavors,op%beta)
535 0 : IF ( op%U .GE. 0.d0 ) THEN
536 0 : CALL ImpurityOperator_computeU(op%Impurity,op%U,0.d0)
537 0 : op%setU = .TRUE.
538 : END IF
539 : ! op%mu = op%mu + op%Impurity%shift_mu
540 : !sui!write(std_out,*) "op%opt_nondiag",op%opt_nondiag
541 0 : CALL BathOperatoroffdiagComplex_init(op%Bath, op%flavors, op%samples, op%beta, INT(buffer(9)), op%opt_nondiag)
542 :
543 0 : op%para = .TRUE.
544 :
545 0 : END SUBROUTINE CtqmcoffdiagComplex_setParameters
546 : !!***
547 :
548 : !!****f* ABINIT/m_CtqmcoffdiagComplex/CtqmcoffdiagComplex_setSweeps
549 : !! NAME
550 : !! CtqmcoffdiagComplex_setSweeps
551 : !!
552 : !! FUNCTION
553 : !! set the number of sweeps
554 : !!
555 : !! COPYRIGHT
556 : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
557 : !! This file is distributed under the terms of the
558 : !! GNU General Public License, see ~abinit/COPYING
559 : !! or http://www.gnu.org/copyleft/gpl.txt .
560 : !!
561 : !! INPUTS
562 : !! op=ctqmc
563 : !! sweeps=asked sweeps
564 : !!
565 : !! OUTPUT
566 : !!
567 : !! SIDE EFFECTS
568 : !!
569 : !! NOTES
570 : !!
571 : !! SOURCE
572 :
573 0 : SUBROUTINE CtqmcoffdiagComplex_setSweeps(op,sweeps)
574 :
575 : !Arguments ------------------------------------
576 : TYPE(CtqmcoffdiagComplex) , INTENT(INOUT) :: op
577 : DOUBLE PRECISION , INTENT(IN ) :: sweeps
578 :
579 0 : op%sweeps = NINT(sweeps / DBLE(op%size))
580 : ! !write(std_out,*) op%sweeps,NINT(sweeps / DBLE(op%size)),ANINT(sweeps/DBLE(op%size))
581 0 : IF ( DBLE(op%sweeps) .NE. ANINT(sweeps/DBLE(op%size)) ) &
582 0 : CALL ERROR("CtqmcoffdiagComplex_setSweeps : sweeps is negative or too big ")
583 0 : IF ( op%sweeps .LT. 2*CTQMC_SLICE1 ) THEN !202
584 0 : CALL WARNALL("CtqmcoffdiagComplex_setSweeps : # sweeps automtically changed ")
585 0 : op%sweeps = 2*CTQMC_SLICE1
586 : ! ELSE IF ( op%sweeps .LT. op%thermalization ) THEN
587 : ! CALL WARNALL("CtqmcoffdiagComplex_setSweeps : Thermalization > sweeps / cpu -> auto fix")
588 : ! op%sweeps = op%thermalization
589 : END IF
590 0 : IF ( DBLE(NINT(DBLE(op%sweeps)*DBLE(op%size)/DBLE(CTQMC_SLICE1))) .NE. &
591 0 : ANINT(DBLE(op%sweeps)*DBLE(op%size)/DBLE(CTQMC_SLICE1)) ) THEN
592 0 : op%modNoise1 = op%sweeps
593 : ELSE
594 0 : op%modNoise1 = MIN(op%sweeps,INT(DBLE(op%sweeps)*DBLE(op%size) / DBLE(CTQMC_SLICE1))) !101
595 : END IF
596 0 : op%modNoise2 = MAX(op%modNoise1 / CTQMC_SLICE2, 1) ! 100
597 : ! op%modGlobalMove(1) = op%thermalization / 10 + 1
598 : ! op%modGlobalMove(2) = 0
599 :
600 0 : END SUBROUTINE CtqmcoffdiagComplex_setSweeps
601 : !!***
602 :
603 : !!****f* ABINIT/m_CtqmcoffdiagComplex/CtqmcoffdiagComplex_setSeed
604 : !! NAME
605 : !! CtqmcoffdiagComplex_setSeed
606 : !!
607 : !! FUNCTION
608 : !! initialize random number generator
609 : !!
610 : !! COPYRIGHT
611 : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
612 : !! This file is distributed under the terms of the
613 : !! GNU General Public License, see ~abinit/COPYING
614 : !! or http://www.gnu.org/copyleft/gpl.txt .
615 : !!
616 : !! INPUTS
617 : !! op=ctqmc
618 : !! iseed=seed from imput
619 : !!
620 : !! OUTPUT
621 : !!
622 : !! SIDE EFFECTS
623 : !!
624 : !! NOTES
625 : !!
626 : !! SOURCE
627 :
628 0 : SUBROUTINE CtqmcoffdiagComplex_setSeed(op,iseed)
629 :
630 : !Arguments ------------------------------------
631 : TYPE(CtqmcoffdiagComplex), INTENT(INOUT) :: op
632 : INTEGER , INTENT(IN ) :: iseed
633 : !Local variables ------------------------------
634 : !INTEGER :: n
635 : !INTEGER :: i
636 : !INTEGER, DIMENSION(:), ALLOCATABLE :: seed
637 :
638 :
639 : !CALL RANDOM_SEED(size = n)
640 : !MALLOC(seed,(n))
641 : !seed = iseed + (/ (i - 1, i = 1, n) /)
642 :
643 : !CALL RANDOM_SEED(PUT = seed+op%rank)
644 :
645 : !FREE(seed)
646 :
647 0 : op%seed=INT(iseed+op%rank,8)
648 :
649 0 : END SUBROUTINE CtqmcoffdiagComplex_setSeed
650 : !!***
651 :
652 : !!****f* ABINIT/m_CtqmcoffdiagComplex/CtqmcoffdiagComplex_allocateAll
653 : !! NAME
654 : !! CtqmcoffdiagComplex_allocateAll
655 : !!
656 : !! FUNCTION
657 : !! Allocate all non option variables
658 : !!
659 : !! COPYRIGHT
660 : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
661 : !! This file is distributed under the terms of the
662 : !! GNU General Public License, see ~abinit/COPYING
663 : !! or http://www.gnu.org/copyleft/gpl.txt .
664 : !!
665 : !! INPUTS
666 : !! op=ctqmc
667 : !!
668 : !! OUTPUT
669 : !!
670 : !! SIDE EFFECTS
671 : !!
672 : !! NOTES
673 : !!
674 : !! SOURCE
675 :
676 0 : SUBROUTINE CtqmcoffdiagComplex_allocateAll(op)
677 :
678 : !Arguments ------------------------------------
679 : TYPE(CtqmcoffdiagComplex), INTENT(INOUT) :: op
680 : !Local variables ------------------------------
681 : INTEGER :: flavors
682 :
683 0 : IF ( .NOT. op%para ) &
684 0 : CALL ERROR("CtqmcoffdiagComplex_allocateAll : CtqmcoffdiagComplex_setParameters never called ")
685 :
686 0 : flavors = op%flavors
687 :
688 :
689 : ! number of electrons
690 0 : FREEIF(op%measN)
691 0 : MALLOC(op%measN,(1:4,1:flavors))
692 0 : op%measN = cmplx(0.d0,0.d0,kind=8)
693 :
694 : ! double occupancies
695 0 : FREEIF(op%measDE)
696 0 : MALLOC(op%measDE,(1:flavors,1:flavors) )
697 0 : op%measDE = real(cmplx(0.d0,0.d0,kind=8))
698 :
699 0 : FREEIF(op%mu)
700 : #ifdef FC_LLVM
701 : ! LLVM 16 doesn't recognize this macro here
702 : MALLOC(op%mu, (1:flavors) )
703 : #else
704 0 : MALLOC(op%mu, (1:flavors))
705 : #endif
706 0 : op%mu = 0.d0
707 0 : FREEIF(op%hybri_limit)
708 : #ifdef FC_LLVM
709 : ! LLVM 16 doesn't recognize this macro here
710 : MALLOC(op%hybri_limit, (flavors,flavors) )
711 : #else
712 0 : MALLOC(op%hybri_limit, (flavors,flavors))
713 : #endif
714 0 : op%hybri_limit = czero
715 0 : END SUBROUTINE CtqmcoffdiagComplex_allocateAll
716 : !!***
717 :
718 : !#ifdef CTCtqmcoffdiagComplex_ANALYSIS
719 : !!****f* ABINIT/m_CtqmcoffdiagComplex/CtqmcoffdiagComplex_allocateOpt
720 : !! NAME
721 : !! CtqmcoffdiagComplex_allocateOpt
722 : !!
723 : !! FUNCTION
724 : !! allocate all option variables
725 : !!
726 : !! COPYRIGHT
727 : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
728 : !! This file is distributed under the terms of the
729 : !! GNU General Public License, see ~abinit/COPYING
730 : !! or http://www.gnu.org/copyleft/gpl.txt .
731 : !!
732 : !! INPUTS
733 : !! op=ctqmc
734 : !!
735 : !! OUTPUT
736 : !!
737 : !! SIDE EFFECTS
738 : !!
739 : !! NOTES
740 : !!
741 : !! SOURCE
742 :
743 0 : SUBROUTINE CtqmcoffdiagComplex_allocateOpt(op)
744 :
745 : !Arguments ------------------------------------
746 : TYPE(CtqmcoffdiagComplex), INTENT(INOUT) :: op
747 : !Local variables ------------------------------
748 : INTEGER :: i
749 : INTEGER :: j
750 : INTEGER :: k
751 :
752 0 : IF ( .NOT. op%para ) &
753 0 : CALL ERROR("CtqmcoffdiagComplex_allocateOpt : CtqmcoffdiagComplex_setParameters never called ")
754 :
755 0 : IF ( op%opt_analysis .EQ. 1 ) THEN
756 0 : FREEIF(op%measCorrelation)
757 0 : MALLOC(op%measCorrelation,(1:op%samples+1,1:3,1:op%flavors))
758 0 : op%measCorrelation = 0.d0
759 : END IF
760 :
761 0 : IF ( op%opt_order .GT. 0 ) THEN
762 0 : FREEIF(op%measPerturbation)
763 0 : MALLOC(op%measPerturbation,(1:op%opt_order,1:op%flavors))
764 0 : op%measPerturbation = 0.d0
765 0 : FREEIF(op%meas_fullemptylines)
766 0 : MALLOC(op%meas_fullemptylines,(2,1:op%flavors))
767 0 : op%meas_fullemptylines = 0.d0
768 : END IF
769 :
770 0 : IF ( op%opt_histo .GT. 0 ) THEN
771 0 : FREEIF(op%occup_histo_time)
772 0 : MALLOC(op%occup_histo_time,(1:op%flavors+1))
773 0 : op%occup_histo_time= 0.d0
774 0 : FREEIF(op%occupconfig)
775 0 : MALLOC(op%occupconfig,(1:2**op%flavors))
776 0 : op%occupconfig= 0.d0
777 0 : FREEIF(op%suscep)
778 0 : MALLOC(op%suscep,(1:3,1:op%samples))
779 0 : op%suscep= 0.d0
780 0 : FREEIF(op%chi)
781 0 : MALLOC(op%chi,(1:3,1:op%samples))
782 0 : op%chi= 0.d0
783 0 : FREEIF(op%chicharge)
784 0 : MALLOC(op%chicharge,(1:3,1:op%samples))
785 0 : op%chicharge= 0.d0
786 0 : FREEIF(op%ntot)
787 0 : MALLOC(op%ntot,(1:3))
788 0 : op%ntot= 0.d0
789 : END IF
790 :
791 0 : IF ( op%opt_noise .EQ. 1 ) THEN
792 0 : IF ( ALLOCATED(op%measNoiseG) ) THEN
793 0 : DO i=1,2
794 0 : DO j = 1, op%flavors
795 0 : DO k= 1, op%samples+1
796 0 : CALL Vector_destroy(op%measNoiseG(k,j,i))
797 : END DO
798 : END DO
799 : END DO
800 0 : DT_FREE(op%measNoiseG)
801 : END IF
802 0 : DT_MALLOC(op%measNoiseG,(1:op%samples+1,1:op%flavors,1:2))
803 : !DO i=1,2
804 0 : DO j = 1, op%flavors
805 0 : DO k= 1, op%samples+1
806 0 : CALL Vector_init(op%measNoiseG(k,j,1),CTQMC_SLICE1)
807 : END DO
808 : END DO
809 0 : DO j = 1, op%flavors
810 0 : DO k= 1, op%samples+1
811 0 : CALL Vector_init(op%measNoiseG(k,j,2),CTQMC_SLICE1*CTQMC_SLICE2+1) ! +1 pour etre remplacer ceil
812 : END DO
813 : END DO
814 : !END DO
815 0 : FREEIF(op%abNoiseG)
816 0 : MALLOC(op%aBNoiseG,(1:2,1:op%samples+1,op%flavors))
817 0 : op%abNoiseG = 0.d0
818 : END IF
819 :
820 0 : IF (op%opt_spectra .GE. 1 ) THEN
821 0 : FREEIF(op%density)
822 : !MALLOC(op%density,(1:op%thermalization,1:op%flavors))
823 0 : i = CEILING(DBLE(op%thermalization+op%sweeps)/DBLE(op%measurements*op%opt_spectra))
824 0 : MALLOC(op%density,(1:op%flavors+1,1:i))
825 0 : op%density = 0.d0
826 : ! op%endDensity=i
827 : !ENDPROBLEM
828 :
829 : END IF
830 : !#endif
831 0 : END SUBROUTINE CtqmcoffdiagComplex_allocateOpt
832 : !!***
833 :
834 : !SUBROUTINE CtqmcoffdiagComplex_setG0wFile(op,istream,opt_fk)
835 : !include 'mpif.h'
836 : !#endif
837 : ! TYPE(CtqmcoffdiagComplex), INTENT(INOUT) :: op
838 : ! INTEGER , INTENT(IN ) :: istream
839 : ! INTEGER , INTENT(IN ) :: opt_fk
840 : !! DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: F
841 : ! COMPLEX*16 , DIMENSION(:,:), ALLOCATABLE :: Gomega
842 : ! INTEGER :: flavors
843 : ! INTEGER :: iflavor
844 : ! INTEGER :: iomega
845 : ! TYPE(GreenHyboffdiag) :: temp
846 : !#ifdef HAVE_MPI
847 : ! INTEGER :: ierr
848 : !#endif
849 : !
850 : ! IF ( .NOT. op%para ) &
851 : ! CALL ERROR("CtqmcoffdiagComplex_setG0wFile : CtqmcoffdiagComplex_setParameters never called ")
852 : !
853 : ! flavors = op%flavors
854 : !
855 : !! MALLOC(F,(1:op%samples+1,1:flavors))
856 : ! MALLOC(Gomega,(1:op%Wmax,1:flavors))
857 : ! IF ( op%have_MPI .EQV. .TRUE. ) THEN
858 : ! CALL GreenHyboffdiag_init(temp, op%samples, op%beta, op%MY_COMM)
859 : ! ELSE
860 : ! CALL GreenHyboffdiag_init(temp, op%samples, op%beta)
861 : ! END IF
862 : !
863 : ! DO iflavor = 1, flavors
864 : ! IF ( op%rank .EQ. 0 ) THEN
865 : ! DO iomega=1, op%Wmax
866 : ! READ(istream,*) Gomega(iomega,iflavor)
867 : ! END DO
868 : ! END IF
869 : ! END DO
870 : !
871 : !#ifdef HAVE_MPI
872 : ! CALL MPI_Bcast(Gomega, op%Wmax*flavors, MPI_DOUBLE_COMPLEX, 0, &
873 : ! op%MY_COMM, ierr)
874 : !#endif
875 : !
876 : !! CALL GreenHyboffdiag_backFourier(temp,Gomega(:,1))
877 : !! WRITE(10,*) Gomega(:,1)
878 : !! CALL GreenHyboffdiag_print(temp)
879 : !! CALL GreenHyboffdiag_forFourier (temp,Gomega(:,1))
880 : !! WRITE(11,*) Gomega(:,1)
881 : ! CALL CtqmcoffdiagComplex_setG0wTab(op,Gomega,opt_fk)
882 : !! CALL CtqmcoffdiagComplex_computeF(op,Gomega, op%Wmax, F) ! mu is changed
883 : ! FREE(Gomega)
884 : !! CALL BathOperatoroffdiag_setF(op%Bath, F)
885 : !! CALL BathOperatoroffdiag_printF(op%Bath)
886 : !! FREE(F)
887 : !
888 : !! stop
889 : !! op%inF = .TRUE.
890 : !! op%set = .TRUE.
891 : !
892 : !END SUBROUTINE CtqmcoffdiagComplex_setG0wFile
893 : !!***
894 :
895 : !!****f* ABINIT/m_CtqmcoffdiagComplex/CtqmcoffdiagComplex_setG0wTab
896 : !! NAME
897 : !! CtqmcoffdiagComplex_setG0wTab
898 : !!
899 : !! FUNCTION
900 : !! Set Gow from input array
901 : !!
902 : !! COPYRIGHT
903 : !! Copyright (C) 2013-2025 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=ctqmc
910 : !! Gomega=G0w
911 : !! opt_fk=F is already inversed with out iwn
912 : !!
913 : !! OUTPUT
914 : !!
915 : !! SIDE EFFECTS
916 : !!
917 : !! NOTES
918 : !!
919 : !! SOURCE
920 :
921 0 : SUBROUTINE CtqmcoffdiagComplex_setG0wTab(op,Gomega,opt_fk,Iatom,fname)
922 :
923 : !Arguments ------------------------------------
924 : TYPE(CtqmcoffdiagComplex), INTENT(INOUT) :: op
925 : COMPLEX(KIND=8), DIMENSION(:,:,:), INTENT(IN ) :: Gomega
926 : INTEGER , INTENT(IN ) :: opt_fk
927 : INTEGER, INTENT(IN) :: Iatom
928 : CHARACTER(LEN=fnlen), INTENT(INOUT) :: fname
929 : !Local variable -------------------------------
930 : COMPLEX(KIND=8), DIMENSION(:,:,:), ALLOCATABLE :: F
931 :
932 0 : IF ( .NOT. op%para ) &
933 0 : CALL ERROR("CtqmcoffdiagComplex_setG0wTab : CtqmcoffdiagComplex_setParameters never called ")
934 :
935 0 : MALLOC(F,(1:op%samples+1,1:op%flavors,1:op%flavors))
936 0 : CALL CtqmcoffdiagComplex_computeF(op,Gomega, F, opt_fk,Iatom,fname) ! mu is changed
937 : !write(6,*) "eee111"
938 0 : CALL BathOperatoroffdiagComplex_setF(op%Bath, F)
939 : ! CALL BathOperatoroffdiag_printF(op%Bath,333)
940 : !write(6,*) "eee"
941 0 : FREE(F)
942 :
943 0 : op%inF = .TRUE.
944 0 : op%set = .TRUE.
945 :
946 0 : END SUBROUTINE CtqmcoffdiagComplex_setG0wTab
947 : !!***
948 :
949 : !SUBROUTINE CtqmcoffdiagComplex_setFwK(op,Gomega)
950 : ! COMPLEX*16 , DIMENSION(:,:), INTENT(IN ) :: Gomega
951 : ! DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: F
952 : !
953 : ! IF ( .NOT. op%para ) &
954 : ! CALL ERROR("CtqmcoffdiagComplex_setG0wTab : CtqmcoffdiagComplex_setParameters never called ")
955 : !
956 : ! MALLOC(F,(1:op%samples+1,1:op%flavors))
957 : ! CALL CtqmcoffdiagComplex_computeFK(op,Gomega, op%Wmax, F) ! mu is changed
958 : ! CALL BathOperatoroffdiag_setF(op%Bath, F)
959 : ! CALL BathOperatoroffdiag_printF(op%Bath)
960 : ! FREE(F)
961 : !
962 : ! op%inF = .TRUE.
963 : ! op%set = .TRUE.
964 : !
965 : !END SUBROUTINE CtqmcoffdiagComplex_setFwK
966 : !!***
967 :
968 : !SUBROUTINE CtqmcoffdiagComplex_setBand(op, mu, U)
969 : ! DOUBLE PRECISION, INTENT(IN ) :: mu
970 : ! DOUBLE PRECISION, INTENT(IN ) :: U
971 : !
972 : ! IF ( .NOT. op%para ) &
973 : ! CALL ERROR("CtqmcoffdiagComplex_setBand : CtqmcoffdiagComplex_setParameters never called ")
974 : !
975 : ! CALL ImpurityOperator_setU(op%Impurity, U, 0.d0)
976 : ! !op%mu = mu + op%Impurity%shift_mu
977 : !END SUBROUTINE CtqmcoffdiagComplex_setBand
978 : !!***
979 :
980 : !!****f* ABINIT/m_CtqmcoffdiagComplex/CtqmcoffdiagComplex_setU
981 : !! NAME
982 : !! CtqmcoffdiagComplex_setU
983 : !!
984 : !! FUNCTION
985 : !! set the interaction matrix
986 : !!
987 : !! COPYRIGHT
988 : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
989 : !! This file is distributed under the terms of the
990 : !! GNU General Public License, see ~abinit/COPYING
991 : !! or http://www.gnu.org/copyleft/gpl.txt .
992 : !!
993 : !! INPUTS
994 : !! op=ctqmc
995 : !! matU=interaction matrix
996 : !!
997 : !! OUTPUT
998 : !!
999 : !! SIDE EFFECTS
1000 : !!
1001 : !! NOTES
1002 : !!
1003 : !! SOURCE
1004 :
1005 0 : SUBROUTINE CtqmcoffdiagComplex_setU(op,matU)
1006 :
1007 : !Arguments ------------------------------------
1008 : TYPE(CtqmcoffdiagComplex), INTENT(INOUT) ::op
1009 : !Local variables ------------------------------
1010 : COMPLEX(KIND=8), DIMENSION(:,:), INTENT(IN) :: matU
1011 :
1012 0 : IF ( SIZE(matU) .NE. op%flavors*op%flavors ) &
1013 0 : CALL ERROR("CtqmcoffdiagComplex_setU : Wrong interaction matrix (size) ")
1014 :
1015 0 : CALL ImpurityOperator_setUmatComplex(op%Impurity, matU)
1016 0 : op%setU = .TRUE.
1017 0 : END SUBROUTINE CtqmcoffdiagComplex_setU
1018 : !!***
1019 :
1020 : !!****f* ABINIT/m_CtqmcoffdiagComplex/CtqmcoffdiagComplex_clear
1021 : !! NAME
1022 : !! CtqmcoffdiagComplex_clear
1023 : !!
1024 : !! FUNCTION
1025 : !! clear a ctqmc run
1026 : !!
1027 : !! COPYRIGHT
1028 : !! Copyright (C) 2013-2025 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=ctqmc
1035 : !!
1036 : !! OUTPUT
1037 : !!
1038 : !! SIDE EFFECTS
1039 : !!
1040 : !! NOTES
1041 : !!
1042 : !! SOURCE
1043 :
1044 0 : SUBROUTINE CtqmcoffdiagComplex_clear(op)
1045 :
1046 : !Arguments ------------------------------------
1047 : TYPE(CtqmcoffdiagComplex), INTENT(INOUT) :: op
1048 : !Local variables ------------------------------
1049 : INTEGER :: i
1050 : INTEGER :: j
1051 : INTEGER :: k
1052 :
1053 0 : op%measN(1,:) = cmplx(0.d0,0.d0,kind=8)
1054 0 : op%measN(2,:) = cmplx(0.d0,0.d0,kind=8)
1055 : !Do not set measN(3,:) to 0 to avoid erasing N between therm and ctqmc
1056 0 : op%measN(4,:) = cmplx(0.d0,0.d0,kind=8)
1057 0 : op%measDE = real(cmplx(0.d0,0.d0,kind=8))
1058 : ! op%seg_added = 0.d0
1059 : ! op%anti_added = 0.d0
1060 : ! op%seg_removed = 0.d0
1061 : ! op%anti_removed = 0.d0
1062 : ! op%seg_sign = 0.d0
1063 : ! op%anti_sign = 0.d0
1064 0 : op%stats(:) = 0.d0
1065 : ! op%signvaluecurrent = 0.d0
1066 : ! op%signvaluemeas = 0.d0
1067 0 : op%swap = 0.d0
1068 0 : op%runTime = 0.d0
1069 0 : op%modGlobalMove(2) = 0
1070 0 : CALL Vector_clear(op%measNoise(1))
1071 0 : CALL Vector_clear(op%measNoise(2))
1072 : !#ifdef CTCtqmcoffdiagComplex_CHECK
1073 0 : op%errorImpurity = 0.d0
1074 0 : op%errorBath = 0.d0
1075 : !#endif
1076 0 : CALL GreenHyboffdiagComplex_clear(op%Greens)
1077 : !#ifdef CTCtqmcoffdiagComplex_ANALYSIS
1078 0 : IF ( op%opt_analysis .EQ. 1 .AND. ALLOCATED(op%measCorrelation) ) &
1079 0 : op%measCorrelation = 0.d0
1080 0 : IF ( op%opt_order .GT. 0 .AND. ALLOCATED(op%measPerturbation) ) &
1081 0 : op%measPerturbation = 0.d0
1082 0 : IF ( op%opt_order .GT. 0 .AND. ALLOCATED(op%meas_fullemptylines) ) &
1083 0 : op%meas_fullemptylines = 0.d0
1084 0 : IF ( op%opt_noise .EQ. 1 .AND. ALLOCATED(op%measNoiseG) ) THEN
1085 0 : DO i=1,2
1086 0 : DO j = 1, op%flavors
1087 0 : DO k= 1, op%samples+1
1088 0 : CALL Vector_clear(op%measNoiseG(k,j,i))
1089 : END DO
1090 : END DO
1091 : END DO
1092 : !DO j = 1, op%flavors
1093 : ! CALL GreenHyboffdiag_clear(op%Greens(j))
1094 : !END DO
1095 : END IF
1096 : !#endif
1097 0 : END SUBROUTINE CtqmcoffdiagComplex_clear
1098 : !!***
1099 :
1100 : !!****f* ABINIT/m_CtqmcoffdiagComplex/CtqmcoffdiagComplex_reset
1101 : !! NAME
1102 : !! CtqmcoffdiagComplex_reset
1103 : !!
1104 : !! FUNCTION
1105 : !! reset a ctqmc simulation
1106 : !!
1107 : !! COPYRIGHT
1108 : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
1109 : !! This file is distributed under the terms of the
1110 : !! GNU General Public License, see ~abinit/COPYING
1111 : !! or http://www.gnu.org/copyleft/gpl.txt .
1112 : !!
1113 : !! INPUTS
1114 : !! op=ctqmc
1115 : !!
1116 : !! OUTPUT
1117 : !!
1118 : !! SIDE EFFECTS
1119 : !!
1120 : !! NOTES
1121 : !!
1122 : !! SOURCE
1123 :
1124 0 : SUBROUTINE CtqmcoffdiagComplex_reset(op)
1125 :
1126 : !Arguments ------------------------------------
1127 : TYPE(CtqmcoffdiagComplex), INTENT(INOUT) :: op
1128 : !Local variables ------------------------------
1129 : !INTEGER :: iflavor
1130 : DOUBLE PRECISION :: sweeps
1131 :
1132 0 : CALL GreenHyboffdiagComplex_reset(op%Greens)
1133 0 : CALL CtqmcoffdiagComplex_clear(op)
1134 0 : CALL ImpurityOperator_reset(op%Impurity)
1135 0 : CALL BathOperatoroffdiagComplex_reset (op%Bath)
1136 0 : op%measN(3,:) = cmplx(0.d0,0.d0,kind=8)
1137 : !complete restart -> measN=0
1138 0 : op%done = .FALSE.
1139 0 : op%set = .FALSE.
1140 0 : op%inF = .FALSE.
1141 0 : op%opt_movie = 0
1142 0 : op%opt_analysis = 0
1143 0 : op%opt_order = 0
1144 0 : op%opt_check = 0
1145 0 : op%opt_noise = 0
1146 0 : op%opt_spectra = 0
1147 0 : op%opt_levels = 0
1148 0 : sweeps = DBLE(op%sweeps)*DBLE(op%size)
1149 0 : CALL CtqmcoffdiagComplex_setSweeps(op, sweeps)
1150 : !#ifdef HAVE_MPI
1151 : ! CALL MPI_BARRIER(op%MY_COMM,iflavor)
1152 : ! IF ( op%rank .EQ. 0 ) &
1153 : !#endif
1154 : ! WRITE(op%ostream,'(A9)') "QMC reset"
1155 : ! CALL FLUSH(op%ostream)
1156 0 : END SUBROUTINE CtqmcoffdiagComplex_reset
1157 : !!***
1158 :
1159 : !!****f* ABINIT/m_CtqmcoffdiagComplex/CtqmcoffdiagComplex_setMu
1160 : !! NAME
1161 : !! CtqmcoffdiagComplex_setMu
1162 : !!
1163 : !! FUNCTION
1164 : !! impose energy levels
1165 : !!
1166 : !! COPYRIGHT
1167 : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
1168 : !! This file is distributed under the terms of the
1169 : !! GNU General Public License, see ~abinit/COPYING
1170 : !! or http://www.gnu.org/copyleft/gpl.txt .
1171 : !!
1172 : !! INPUTS
1173 : !! op=ctqmc
1174 : !! levels=energy levels vector
1175 : !!
1176 : !! OUTPUT
1177 : !! argout(sizeout)=description
1178 : !!
1179 : !! SIDE EFFECTS
1180 : !!
1181 : !! NOTES
1182 : !!
1183 : !! SOURCE
1184 :
1185 0 : SUBROUTINE CtqmcoffdiagComplex_setMu(op, levels)
1186 :
1187 : !Arguments ------------------------------------
1188 : TYPE(CtqmcoffdiagComplex) , INTENT(INOUT) :: op
1189 : COMPLEX(KIND=8), DIMENSION(:), INTENT(IN ) :: levels
1190 :
1191 0 : IF ( op%flavors .NE. SIZE(levels,1) ) &
1192 0 : CALL WARNALL("CtqmcoffdiagComplex_setMu : Taking energy levels from weiss G(iw)")
1193 :
1194 0 : op%mu(:)=-levels(:) ! levels = \epsilon_j - \mu
1195 : !op%mu =\tilde{\mu} = \mu -\epsilon_j
1196 0 : op%opt_levels = 1
1197 0 : END SUBROUTINE CtqmcoffdiagComplex_setMu
1198 : !!***
1199 :
1200 : !!****f* ABINIT/m_CtqmcoffdiagComplex/CtqmcoffdiagComplex_sethybri_limit
1201 : !! NAME
1202 : !! CtqmcoffdiagComplex_sethybri_limit
1203 : !!
1204 : !! FUNCTION
1205 : !! use coefficient A such that F=-A/(iwn) given by DMFT code.
1206 : !!
1207 : !! COPYRIGHT
1208 : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
1209 : !! This file is distributed under the terms of the
1210 : !! GNU General Public License, see ~abinit/COPYING
1211 : !! or http://www.gnu.org/copyleft/gpl.txt .
1212 : !!
1213 : !! INPUTS
1214 : !! hybri_limit(nflavor,nflavor)=contains the limit for each couple of flavors
1215 : !!
1216 : !! OUTPUT
1217 : !! argout(sizeout)=description
1218 : !!
1219 : !! SIDE EFFECTS
1220 : !! op(CtqmcoffdiagComplex_type) = is the ctqmc main variable
1221 : !! op&limit is now filled
1222 : !! NOTES
1223 : !!
1224 : !! SOURCE
1225 :
1226 0 : SUBROUTINE CtqmcoffdiagComplex_sethybri_limit(op, hybri_limit,opthybri)
1227 :
1228 : !Arguments ------------------------------------
1229 : TYPE(CtqmcoffdiagComplex) , INTENT(INOUT) :: op
1230 : COMPLEX(KIND=8) , DIMENSION(:,:), INTENT(IN ) :: hybri_limit
1231 : INTEGER, INTENT(IN) :: opthybri
1232 : !----------------------------------------------
1233 0 : IF ( op%flavors .NE. SIZE(hybri_limit,1) ) &
1234 0 : CALL ERROR("Error in sethybri_limit")
1235 :
1236 0 : op%hybri_limit(:,:)=hybri_limit(:,:)
1237 0 : op%opt_hybri_limit = opthybri
1238 0 : END SUBROUTINE CtqmcoffdiagComplex_sethybri_limit
1239 : !!***
1240 : !!****f* ABINIT/m_CtqmcoffdiagComplex/CtqmcoffdiagComplex_computeF
1241 : !! NAME
1242 : !! CtqmcoffdiagComplex_computeF
1243 : !!
1244 : !! FUNCTION
1245 : !! Compute the hybridization function
1246 : !!
1247 : !! COPYRIGHT
1248 : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
1249 : !! This file is distributed under the terms of the
1250 : !! GNU General Public License, see ~abinit/COPYING
1251 : !! or http://www.gnu.org/copyleft/gpl.txt .
1252 : !!
1253 : !! INPUTS
1254 : !! op=ctqmc
1255 : !! Gomega=G0 to compute F
1256 : !! opt_fk=What is Gomega
1257 : !!
1258 : !! OUTPUT
1259 : !! F=hybridization function
1260 : !!
1261 : !! SIDE EFFECTS
1262 : !!
1263 : !! NOTES
1264 : !!
1265 : !! SOURCE
1266 :
1267 0 : SUBROUTINE CtqmcoffdiagComplex_computeF(op, Gomega, F, opt_fk,Iatom,fname)
1268 :
1269 : use m_hide_lapack, only : xginv
1270 : !Arguments ------------------------------------
1271 : TYPE(CtqmcoffdiagComplex) , INTENT(INOUT) :: op
1272 : COMPLEX(KIND=8), DIMENSION(:,:,:), INTENT(IN ) :: Gomega
1273 : !INTEGER , INTENT(IN ) :: Wmax
1274 : COMPLEX(KIND=8), DIMENSION(:,:,:), INTENT(INOUT) :: F
1275 : INTEGER , INTENT(IN ) :: opt_fk
1276 : INTEGER, INTENT(IN) :: Iatom
1277 : CHARACTER(LEN=fnlen), INTENT(INOUT) :: fname
1278 : !Local variables ------------------------------
1279 : INTEGER :: flavors
1280 : INTEGER :: samples
1281 : INTEGER :: iflavor
1282 : INTEGER :: iflavor2
1283 : INTEGER :: iomega
1284 : INTEGER :: itau,ioerr,debug
1285 : DOUBLE PRECISION :: pi_invBeta
1286 : COMPLEX(KIND=8) :: K
1287 : !DOUBLE PRECISION :: re
1288 : !DOUBLE PRECISION :: im
1289 : !DOUBLE PRECISION :: det
1290 0 : COMPLEX(KIND=8), DIMENSION(:,:,:), ALLOCATABLE :: F_omega
1291 0 : COMPLEX(KIND=8), DIMENSION(:,:), ALLOCATABLE :: F_omega_inv
1292 0 : COMPLEX(KIND=8), DIMENSION(:,:,:), ALLOCATABLE :: Gomega_tmp
1293 0 : TYPE(GreenHyboffdiagComplex) :: F_tmp
1294 : CHARACTER(LEN=200) :: message
1295 : character(len=2) :: atomnb
1296 0 : DOUBLE PRECISION, allocatable :: x_r(:,:,:), x_i(:,:,:)
1297 : CHARACTER(LEN=256) :: dummy
1298 : DOUBLE PRECISION :: tautemp
1299 : !character(len=30) :: tmpfil
1300 : !INTEGER :: unitnb
1301 :
1302 : ABI_UNUSED((/opt_fk/))
1303 :
1304 0 : flavors = op%flavors
1305 :
1306 0 : samples = op%samples
1307 0 : pi_invBeta = ACOS(-1.d0) / op%beta
1308 0 : op%Wmax=SIZE(Gomega,1)
1309 :
1310 0 : debug = 0 !set to one to build G0
1311 :
1312 : !=================================
1313 : ! --- Initialize F_tmp
1314 : !=================================
1315 0 : IF ( op%have_MPI .EQV. .TRUE. ) THEN
1316 0 : CALL GreenHyboffdiagComplex_init(F_tmp,samples,op%beta,flavors,MY_COMM=op%MY_COMM)
1317 : ELSE
1318 0 : CALL GreenHyboffdiagComplex_init(F_tmp,samples,op%beta,flavors)
1319 : END IF
1320 : ! K = op%mu
1321 :
1322 : !=================================
1323 : ! --- Allocate F_omega
1324 : !=================================
1325 0 : MALLOC(F_omega,(1:op%Wmax,1:flavors,1:flavors))
1326 0 : MALLOC(F_omega_inv,(1:flavors,1:flavors))
1327 0 : MALLOC(Gomega_tmp,(1:op%Wmax,1:flavors,1:flavors))
1328 0 : Gomega_tmp=Gomega
1329 :
1330 :
1331 : !This section of the code is not used / Only building G0 so just skip it
1332 : if(debug == 1) then
1333 : !=============================================================================================
1334 : ! --- Compute Bath Green's function from Hybridization function in imaginary time
1335 : !=============================================================================================
1336 : ! In the following we want to get G_0(iw_n) from F(iw_n).
1337 : ! Be careful : G_0 is called F and F is called G_0 (Why ????)
1338 : !-----------------------------------------------------------
1339 : ! G_0(iw_n)^-1 = iw_n + mu - F(iw_n)
1340 : !-----------------------------------
1341 : do iomega=1,op%Wmax
1342 : do iflavor=1,flavors
1343 : do iflavor2=1,flavors
1344 : if (iflavor==iflavor2) then
1345 : F_omega_inv(iflavor,iflavor2)= (cmplx(0.d0,(2.d0*DBLE(iomega)-1.d0) * pi_invBeta,kind=8) &
1346 : & + op%mu(iflavor)- Gomega_tmp(iomega,iflavor,iflavor2))
1347 : else
1348 : F_omega_inv(iflavor,iflavor2)= (- Gomega_tmp(iomega,iflavor,iflavor2))
1349 : endif
1350 : enddo
1351 : enddo
1352 :
1353 : !Inverse of G_0(iw_n)^-1 to get G_0(iw_n)
1354 : !--------------------
1355 : call xginv(F_omega_inv,flavors)
1356 :
1357 : do iflavor=1,flavors
1358 : do iflavor2=1,flavors
1359 : F_omega(iomega,iflavor,iflavor2) = F_omega_inv(iflavor,iflavor2)
1360 : enddo
1361 : enddo
1362 : enddo !iomega
1363 :
1364 : CALL GreenHyboffdiagComplex_setOperW(F_tmp,F_omega)
1365 :
1366 : CALL GreenHyboffdiagComplex_backFourier(F_tmp,func="green")
1367 :
1368 : !Put the result in F(tau) aka G_0(tau).
1369 : !--------------------------------------
1370 : DO iflavor = 1, flavors
1371 : DO iflavor2 = 1, flavors
1372 : DO itau=1,samples+1
1373 : ! This symetrization is general and valid even with SOC
1374 : ! Without SOC, it is useless ?
1375 : F(itau,iflavor,iflavor2) = ( F_tmp%oper(itau,iflavor,iflavor2) + &
1376 : & dconjg(F_tmp%oper(itau,iflavor2,iflavor)) )/2.d0
1377 : END DO
1378 : END DO
1379 : END DO
1380 :
1381 : !Print G_0(tau) in file
1382 : !------------------------
1383 : !Real part
1384 : open (unit=4367,file='G0tau_fromF_Re',status='unknown',form='formatted')
1385 : rewind(4367)
1386 : IF ( op%rank .EQ. 0 ) THEN
1387 : DO iflavor = 1, flavors
1388 : DO iflavor2 = 1, flavors
1389 : write(4367,*) "#",iflavor,iflavor2
1390 : do itau=1,op%samples+1
1391 : write(4367,*) (itau-1)*op%beta/(op%samples),real(F(itau,iflavor,iflavor2))
1392 : enddo
1393 : write(4367,*)
1394 : END DO
1395 : END DO
1396 : ENDIF
1397 : !call flush(4367) nag compiler problrem
1398 : close(4367)
1399 :
1400 : !Imag part
1401 : open (unit=4367,file='G0tau_fromF_Im',status='unknown',form='formatted')
1402 : rewind(4367)
1403 : IF ( op%rank .EQ. 0 ) THEN
1404 : DO iflavor = 1, flavors
1405 : DO iflavor2 = 1, flavors
1406 : write(4367,*) "#",iflavor,iflavor2
1407 : do itau=1,op%samples+1
1408 : write(4367,*) (itau-1)*op%beta/(op%samples),aimag(F(itau,iflavor,iflavor2))
1409 : enddo
1410 : write(4367,*)
1411 : END DO
1412 : END DO
1413 : ENDIF
1414 : !call flush(4367) !nag compiler problem
1415 : close(4367)
1416 :
1417 : endif ! end of building G0
1418 :
1419 0 : call xmpi_barrier(op%MY_COMM)
1420 :
1421 : !=============================================================================================
1422 : ! --- Restore Hybridization in F_omega
1423 : !=============================================================================================
1424 :
1425 : ! Restore Hybridization in F_omega for the following operations
1426 0 : F_omega = Gomega_tmp
1427 : ! From now F_omega is actually F(iw_n) calculated in m_forctqmc.f90
1428 : ! and corresponding to F(iw_n) = -[G_0(iw_n^-1 - iwn - levels)]
1429 :
1430 : !==================================================================
1431 : ! --- Full double loop on flavors to compute F (remove levels)
1432 : !==================================================================
1433 0 : DO iflavor = 1, flavors
1434 0 : DO iflavor2 = 1, flavors
1435 :
1436 : ! --- Compute or use the levels for the diagonal hybridization (else K=0)
1437 0 : IF(iflavor==iflavor2) THEN
1438 0 : IF ( op%opt_levels .EQ. 1 ) THEN
1439 : K = op%mu(iflavor)
1440 : ELSE
1441 0 : K = -(F_omega(op%Wmax,iflavor,iflavor))
1442 0 : op%mu(iflavor) = K
1443 : END IF
1444 : ELSE
1445 : K=cmplx(0.d0,0.d0,kind=8)
1446 : ENDIF
1447 :
1448 : ! --- compute residual K (?)
1449 : ! K corresponds to the first moment of hybridization (d1 in setMuD1)
1450 : ! ------------------------------------------------------------------
1451 0 : K = CMPLX(0,(2.d0*DBLE(op%Wmax)-1.d0)*pi_invBeta,8)*F_omega(op%Wmax,iflavor,iflavor2)
1452 0 : CALL GreenHyboffdiagComplex_setMuD1(op%Greens,iflavor,iflavor2,op%mu(iflavor),K)
1453 : END DO
1454 : END DO
1455 :
1456 : ! --- Creates F_tmp%oper_w
1457 0 : CALL GreenHyboffdiagComplex_setOperW(F_tmp,F_omega)
1458 :
1459 : ! For all iflavor and iflavor2, do the Fourier transformation to have F(tau)
1460 : ! ---------------------------------------------
1461 0 : if (op%opt_hybri_limit .eq. 0) then
1462 0 : write(message,'(a,a)') " == Not using the asymptotic limit of hybridization to enforce F(iw_n) -> -C_ij/iw_n (dmft_hybri_limit = 0)", ch10
1463 : else
1464 : ! Take into account asymptotic limit of hybridization function such that F(iw_n) -> -C_ij/iw_n
1465 : ! with C_ij calculated in m_forctqmc.f90
1466 : ! --------------------------------------
1467 0 : write(message,'(a,a)') " == Use asymptotic limit of hybridization function such that F(iw_n) -> -C_ij/iw_n (dmft_hybri_limit = 1)", ch10
1468 : endif
1469 0 : CALL wrtout(std_out,message,'COLL')
1470 0 : CALL GreenHyboffdiagComplex_backFourierComplex(F_tmp,hybri_limit=op%hybri_limit,opt_hybri_limit=op%opt_hybri_limit)
1471 :
1472 : ! --- Put the result in F
1473 0 : DO iflavor = 1, flavors
1474 0 : DO iflavor2 = 1, flavors
1475 0 : DO itau=1,samples+1
1476 0 : F(itau,iflavor,iflavor2) = -F_tmp%oper(samples+2-itau,iflavor,iflavor2)
1477 : END DO
1478 : END DO
1479 : END DO
1480 :
1481 0 : DO iflavor = 1, flavors
1482 0 : DO iflavor2 = 1, flavors
1483 0 : DO itau=1,samples+1
1484 : ! This symetrization is general and valid even with SOC
1485 : ! Without SOC, it leads to zero.
1486 : F(itau,iflavor,iflavor2) = -( (F_tmp%oper(samples+2-itau,iflavor,iflavor2)) + &
1487 0 : & dconjg(F_tmp%oper(samples+2-itau,iflavor2,iflavor)) )/2.d0
1488 : END DO
1489 : END DO
1490 : END DO
1491 :
1492 : !debug
1493 : if (3 .eq. 4) then
1494 : !== Read Hybridization function from file ==
1495 : !reset previous hybridization
1496 : write(message,'(a)') "Warning: About to read Hybridization from file!"
1497 : call wrtout(std_out,message,'COLL')
1498 : DO iflavor = 1, flavors
1499 : DO iflavor2 = 1, flavors
1500 : DO itau=1,samples+1
1501 : F(itau,iflavor,iflavor2) = czero
1502 : END DO
1503 : END DO
1504 : END DO
1505 :
1506 : if (op%rank .eq. 0 ) then
1507 : ABI_MALLOC(x_r,(op%samples+1,flavors,flavors))
1508 : ABI_MALLOC(x_i,(op%samples+1,flavors,flavors))
1509 : open(unit=735,file='Hybridization.dat',status='old',form='formatted',action='read',iostat=ioerr)
1510 :
1511 : read(735, '(a)') dummy
1512 :
1513 : do itau=1, op%samples+1
1514 : read(735, '(2x,393(e25.17e3,2x))') tautemp, &
1515 : & ((x_r(itau,iflavor,iflavor2), x_i(itau,iflavor,iflavor2), iflavor=1, flavors), iflavor2=1, flavors)
1516 : end do
1517 :
1518 : close(735)
1519 :
1520 : ! Now assign to F array (reverse the sign flip and time ordering)
1521 : do itau=1, op%samples+1
1522 : do iflavor=1, flavors
1523 : do iflavor2=1, flavors
1524 : F(op%samples+2-itau, iflavor, iflavor2) = &
1525 : CMPLX(-x_r(itau,iflavor,iflavor2), -x_i(itau,iflavor,iflavor2), kind=8)
1526 : enddo
1527 : enddo
1528 : enddo
1529 : ABI_FREE(x_r)
1530 : ABI_FREE(x_i)
1531 : endif
1532 : end if ! end reading F from file
1533 :
1534 : !== Write Hybridization function in file ==
1535 0 : If (Iatom .lt. 10) then
1536 0 : write(atomnb,'("0",i1)') Iatom
1537 : else
1538 0 : write(atomnb,'(i2)') Iatom
1539 : endif
1540 :
1541 0 : IF ( op%rank .EQ. 0 ) THEN
1542 0 : open (unit=735,file=trim(fname)//'_Hybridization_iatom_'//atomnb//'.dat',status='unknown',form='formatted')
1543 0 : write(735,'(6a)') "# Real and Imaginary part of the Hybridization function Delta(tau) in the CTQMC basis"
1544 0 : do itau=1,op%samples+1
1545 0 : write(735,'(2x,393(e25.17e3,2x))') DBLE(itau-1)*(op%beta/op%samples),&
1546 0 : &((-1*real(F(op%samples+2-itau,iflavor,iflavor2)),-1*AIMAG(F(op%samples+2-itau,iflavor,iflavor2)),iflavor=1,flavors),iflavor2=1,flavors)
1547 : enddo
1548 0 : close(735)
1549 : ENDIF
1550 :
1551 :
1552 : ! call xmpi_barrier(op%MY_COMM)
1553 : !write(6,*) "QQQQ3"
1554 0 : FREE(Gomega_tmp)
1555 0 : FREE(F_omega)
1556 0 : FREE(F_omega_inv)
1557 : !write(6,*) "QQQQ4"
1558 0 : CALL GreenHyboffdiagComplex_destroy(F_tmp)
1559 : !write(6,*) "QQQQ2"
1560 :
1561 :
1562 0 : END SUBROUTINE CtqmcoffdiagComplex_computeF
1563 : !!***
1564 :
1565 : !SUBROUTINE CtqmcoffdiagComplex_computeFK(op, Gomega, Wmax, F)
1566 : ! COMPLEX*16 , DIMENSION(:,:), INTENT(IN ) :: Gomega
1567 : ! INTEGER , INTENT(IN ) :: Wmax
1568 : ! DOUBLE PRECISION, DIMENSION(:,:), INTENT(INOUT) :: F
1569 : ! INTEGER :: flavors
1570 : ! INTEGER :: samples
1571 : ! INTEGER :: iflavor
1572 : ! INTEGER :: iomega
1573 : ! INTEGER :: itau
1574 : ! DOUBLE PRECISION :: pi_invBeta
1575 : ! DOUBLE PRECISION :: K
1576 : ! COMPLEX*16 , DIMENSION(:,:), ALLOCATABLE :: F_omega
1577 : ! TYPE(GreenHyboffdiag) :: F_tmp
1578 : !
1579 : ! flavors = op%flavors
1580 : !
1581 : ! samples = op%samples
1582 : ! pi_invBeta = ACOS(-1.d0) / op%beta
1583 : !
1584 : ! IF ( op%have_MPI .EQV. .TRUE. ) THEN
1585 : ! CALL GreenHyboffdiag_init(F_tmp,samples,op%beta,op%MY_COMM)
1586 : ! ELSE
1587 : ! CALL GreenHyboffdiag_init(F_tmp,samples,op%beta)
1588 : ! END IF
1589 : !! K = op%mu
1590 : !
1591 : ! MALLOC(F_omega,(1:Wmax,1:flavors))
1592 : !
1593 : ! DO iflavor = 1, flavors
1594 : ! K = REAL(Gomega(Wmax, iflavor))
1595 : ! WRITE(op%ostream,*) "CTQMC K, op%mu = ",K,op%mu
1596 : ! WRITE(op%ostream,*) "CTQMC beta = ",op%beta
1597 : ! op%mu(iflavor) = K
1598 : ! DO iomega = 1, Wmax
1599 : ! F_omega(iomega,iflavor) = Gomega(iomega,iflavor) &
1600 : ! - CMPLX(K, 0.d0, 8)
1601 : ! !if(iflavor==1.and.op%rank==0) then
1602 : ! !write(224,*) (2.d0*DBLE(iomega)-1.d0) * pi_invBeta, real(F_omega(iomega,iflavor)),imag(F_omega(iomega,iflavor))
1603 : ! !write(225,*) (2.d0*DBLE(iomega)-1.d0) * pi_invBeta, real(Gomega(iomega, iflavor)),imag(Gomega(iomega, iflavor))
1604 : ! !end if
1605 : ! END DO
1606 : ! CALL GreenHyboffdiag_backFourier(F_tmp,F_omega(:,iflavor))
1607 : ! F(1:samples+1,iflavor) = (/ (-F_tmp%oper(samples+1-itau),itau=0,samples) /)
1608 : ! END DO
1609 : ! FREE(F_omega)
1610 : ! CALL GreenHyboffdiag_destroy(F_tmp)
1611 : !END SUBROUTINE CtqmcoffdiagComplex_computeFK
1612 : !!***
1613 :
1614 : !!****f* ABINIT/m_CtqmcoffdiagComplex/CtqmcoffdiagComplex_run
1615 : !! NAME
1616 : !! CtqmcoffdiagComplex_run
1617 : !!
1618 : !! FUNCTION
1619 : !! set all options and run a simulation
1620 : !!
1621 : !! COPYRIGHT
1622 : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
1623 : !! This file is distributed under the terms of the
1624 : !! GNU General Public License, see ~abinit/COPYING
1625 : !! or http://www.gnu.org/copyleft/gpl.txt .
1626 : !!
1627 : !! INPUTS
1628 : !! op=ctqmc
1629 : !! opt_order=maximal perturbation order to scope
1630 : !! opt_movie=draw a movie of the simulation
1631 : !! opt_analysis=compute correlation functions
1632 : !! opt_check=check fast calculations
1633 : !! opt_noise=compute noise for green function
1634 : !! opt_spectra=fourier transform of the time evolution of the number of electrons
1635 : !! opt_gMove=steps without global move
1636 : !!
1637 : !! OUTPUT
1638 : !!
1639 : !! SIDE EFFECTS
1640 : !!
1641 : !! NOTES
1642 : !!
1643 : !! SOURCE
1644 :
1645 0 : SUBROUTINE CtqmcoffdiagComplex_run(op,opt_order,opt_histo,opt_movie,opt_analysis,opt_check,opt_noise,opt_spectra,opt_gMove)
1646 :
1647 :
1648 : #ifdef HAVE_MPI1
1649 : include 'mpif.h'
1650 : #endif
1651 : !Arguments ------------------------------------
1652 : TYPE(CtqmcoffdiagComplex), INTENT(INOUT) :: op
1653 : INTEGER, OPTIONAL, INTENT(IN ) :: opt_order
1654 : INTEGER, OPTIONAL, INTENT(IN ) :: opt_histo
1655 : INTEGER, OPTIONAL, INTENT(IN ) :: opt_movie
1656 : INTEGER, OPTIONAL, INTENT(IN ) :: opt_analysis
1657 : INTEGER, OPTIONAL, INTENT(IN ) :: opt_check
1658 : INTEGER, OPTIONAL, INTENT(IN ) :: opt_noise
1659 : INTEGER, OPTIONAL, INTENT(IN ) :: opt_spectra
1660 : INTEGER, OPTIONAL, INTENT(IN ) :: opt_gMove
1661 : !Local variables ------------------------------
1662 : #ifdef HAVE_MPI
1663 : INTEGER :: ierr
1664 : DOUBLE PRECISION :: rtime(1)
1665 : #endif
1666 : !#ifdef CTCtqmcoffdiagComplex_MOVIE
1667 : INTEGER :: ilatex
1668 : CHARACTER(LEN=4) :: Cchar
1669 : !#endif
1670 : DOUBLE PRECISION :: estimatedTime
1671 : CHARACTER(LEN=100) :: message
1672 :
1673 0 : IF ( .NOT. op%set ) &
1674 0 : CALL ERROR("CtqmcoffdiagComplex_run : QMC not set up ")
1675 0 : IF ( .NOT. op%setU ) &
1676 0 : CALL ERROR("CtqmcoffdiagComplex_run : QMC does not have a U matrix ")
1677 :
1678 :
1679 : ! OPTIONS of the run
1680 0 : IF ( PRESENT( opt_check ) ) THEN
1681 0 : op%opt_check = opt_check
1682 0 : CALL ImpurityOperator_doCheck(op%Impurity,opt_check)
1683 0 : CALL BathOperatoroffdiagComplex_doCheck(op%Bath,opt_check)
1684 : END IF
1685 0 : IF ( PRESENT( opt_movie ) ) &
1686 0 : op%opt_movie = opt_movie
1687 0 : IF ( PRESENT( opt_analysis ) ) &
1688 0 : op%opt_analysis = opt_analysis
1689 0 : IF ( PRESENT ( opt_order ) ) &
1690 0 : op%opt_order = opt_order
1691 0 : IF ( PRESENT ( opt_histo ) ) &
1692 0 : op%opt_histo = opt_histo
1693 0 : IF ( PRESENT ( opt_noise ) ) THEN
1694 0 : op%opt_noise = opt_noise
1695 : END IF
1696 0 : IF ( PRESENT ( opt_spectra ) ) &
1697 0 : op%opt_spectra = opt_spectra
1698 :
1699 0 : op%modGlobalMove(1) = max(op%sweeps,op%thermalization)+1 ! No Global Move
1700 : !!sui!write(std_out,*) "op%sweeps",op%thermalization,op%sweeps,opt_gMove
1701 0 : op%modGlobalMove(2) = 0
1702 0 : IF ( PRESENT ( opt_gMove ) ) THEN
1703 0 : IF ( opt_gMove .LE. 0 .OR. opt_gMove .GT. op%sweeps ) THEN
1704 : ! op%modGlobalMove(1) = op%sweeps+1
1705 : op%modGlobalMove(1) = max(op%sweeps,op%thermalization)+1 ! No Global Move
1706 : !write(std_out,*) "op%sweeps",op%sweeps, op%modGlobalMove(1)
1707 : !CALL WARNALL("CtqmcoffdiagComplex_run : global moves option is <= 0 or > sweeps/cpu -> No global Moves")
1708 0 : write(message,'(a,a)') " == No global moves are used in CT-QMC. Make sure it is correct", ch10
1709 0 : CALL wrtout(std_out,message,'COLL')
1710 : ELSE
1711 0 : op%modGlobalMove(1) = opt_gMove
1712 : END IF
1713 : END IF
1714 : !sui!write(std_out,*) "op%sweeps",op%thermalization,op%sweeps
1715 :
1716 0 : CALL CtqmcoffdiagComplex_allocateOpt(op)
1717 :
1718 : !#ifdef CTCtqmcoffdiagComplex_MOVIE
1719 0 : ilatex = 0
1720 0 : IF ( op%opt_movie .EQ. 1 ) THEN
1721 0 : Cchar ="0000"
1722 0 : WRITE(Cchar,'(I4)') op%rank
1723 0 : ilatex = 87+op%rank
1724 0 : OPEN(UNIT=ilatex, FILE="Movie_"//TRIM(ADJUSTL(Cchar))//".tex")
1725 0 : WRITE(ilatex,'(A)') "\documentclass{beamer}"
1726 0 : WRITE(ilatex,'(A)') "\usepackage{color}"
1727 0 : WRITE(ilatex,'(A)') "\setbeamersize{sidebar width left=0pt}"
1728 0 : WRITE(ilatex,'(A)') "\setbeamersize{sidebar width right=0pt}"
1729 0 : WRITE(ilatex,'(A)') "\setbeamersize{text width left=0pt}"
1730 0 : WRITE(ilatex,'(A)') "\setbeamersize{text width right=0pt}"
1731 0 : WRITE(ilatex,*)
1732 0 : WRITE(ilatex,'(A)') "\begin{document}"
1733 0 : WRITE(ilatex,*)
1734 : END IF
1735 : !#endif
1736 :
1737 0 : IF ( op%rank .EQ. 0 ) THEN
1738 0 : WRITE(op%ostream,'(A25)') "== Starting QMC Solver"
1739 : END IF
1740 :
1741 : !=================================
1742 : ! STARTING THERMALIZATION
1743 : !=================================
1744 : !write(std_out,*) "sweeps before thermalization",op%sweeps
1745 : !write(std_out,*) "op%stats",op%stats
1746 0 : CALL CtqmcoffdiagComplex_loop(op,op%thermalization,ilatex)
1747 : !=================================
1748 : ! ENDING THERMALIZATION
1749 : !=================================
1750 :
1751 0 : estimatedTime = op%runTime
1752 : #ifdef HAVE_MPI
1753 0 : CALL MPI_REDUCE([op%runTime], rtime, 1, MPI_DOUBLE_PRECISION, MPI_MAX, 0, op%MY_COMM, ierr)
1754 0 : estimatedTime=rtime(1)
1755 : #endif
1756 :
1757 0 : IF ( op%rank .EQ. 0 ) THEN
1758 0 : WRITE(op%ostream,'(A26,I6,A11)') "Thermalization done in ", CEILING(estimatedTime), " seconds"
1759 0 : WRITE(op%ostream,'(A25,I7,A15,I5,A5)') "The QMC should run in ", &
1760 0 : CEILING(estimatedTime*DBLE(op%sweeps)/DBLE(op%thermalization)),&
1761 0 : " seconds on ", op%size, " CPUs"
1762 : END IF
1763 :
1764 : !=================================
1765 : ! CLEANING CTQMC
1766 : !=================================
1767 0 : CALL CtqmcoffdiagComplex_clear(op)
1768 :
1769 : !=================================
1770 : ! STARTING CTQMC
1771 : !=================================
1772 : !write(std_out,*) "sweeps before loop",op%sweeps
1773 : !write(std_out,*) "op%stats",op%stats
1774 0 : CALL CtqmcoffdiagComplex_loop(op,op%sweeps,ilatex)
1775 : !=================================
1776 : ! ENDING CTQMC
1777 : !=================================
1778 :
1779 0 : IF ( op%opt_movie .EQ. 1 ) THEN
1780 0 : WRITE(ilatex,*) ""
1781 0 : WRITE(ilatex,'(A14)') "\end{document}"
1782 0 : CLOSE(ilatex)
1783 : END IF
1784 :
1785 0 : op%done = .TRUE.
1786 : !sui!write(std_out,*) "op%stats en of ctqmc_run",op%stats
1787 :
1788 0 : END SUBROUTINE CtqmcoffdiagComplex_run
1789 : !!***
1790 :
1791 : !!****f* ABINIT/m_CtqmcoffdiagComplex/CtqmcoffdiagComplex_loop
1792 : !! NAME
1793 : !! CtqmcoffdiagComplex_loop
1794 : !!
1795 : !! FUNCTION
1796 : !! Definition the main loop of the CT-QMC
1797 : !!
1798 : !! COPYRIGHT
1799 : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
1800 : !! This file is distributed under the terms of the
1801 : !! GNU General Public License, see ~abinit/COPYING
1802 : !! or http://www.gnu.org/copyleft/gpl.txt .
1803 : !!
1804 : !! INPUTS
1805 : !! op=ctqmc
1806 : !! itotal=number of sweeps to perform : thermalization or sweeps
1807 : !! ilatex=unit of file to write movie if so
1808 : !!
1809 : !! OUTPUT
1810 : !!
1811 : !! SIDE EFFECTS
1812 : !!
1813 : !! NOTES
1814 : !!
1815 : !! SOURCE
1816 :
1817 0 : SUBROUTINE CtqmcoffdiagComplex_loop(op,itotal,ilatex)
1818 :
1819 : !Arguments ------------------------------------
1820 : TYPE(CtqmcoffdiagComplex), INTENT(INOUT) :: op
1821 : INTEGER , INTENT(IN ) :: itotal
1822 : INTEGER , INTENT(IN ) :: ilatex
1823 : !Local variables ------------------------------
1824 : LOGICAL :: updated
1825 : LOGICAL :: updated_seg
1826 0 : LOGICAL, DIMENSION(:), ALLOCATABLE :: updated_swap
1827 :
1828 : INTEGER :: flavors
1829 : INTEGER :: measurements
1830 : INTEGER :: modNoise1
1831 : INTEGER :: modNoise2
1832 : INTEGER :: modGlobalMove
1833 : INTEGER :: sp1
1834 : INTEGER :: itau
1835 : INTEGER :: ind
1836 : INTEGER :: endDensity
1837 : INTEGER :: indDensity
1838 : INTEGER :: swapUpdate1
1839 : INTEGER :: swapUpdate2
1840 : INTEGER :: old_percent
1841 : INTEGER :: new_percent
1842 : INTEGER :: ipercent !,ii
1843 : INTEGER :: iflavor,ifl1,iflavor_d
1844 : INTEGER :: isweep
1845 :
1846 : DOUBLE PRECISION :: cpu_time1
1847 : DOUBLE PRECISION :: cpu_time2
1848 : DOUBLE PRECISION :: NRJ_old1
1849 : DOUBLE PRECISION :: NRJ_old2
1850 : DOUBLE PRECISION :: NRJ_new
1851 : DOUBLE PRECISION :: total
1852 0 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: gtmp_old1
1853 0 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: gtmp_old2
1854 0 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: gtmp_new
1855 :
1856 0 : CALL CPU_TIME(cpu_time1)
1857 :
1858 0 : flavors = op%flavors
1859 0 : measurements = op%measurements
1860 0 : modNoise1 = op%modNoise1
1861 0 : modNoise2 = op%modNoise2
1862 0 : modGlobalMove = op%modGlobalMove(1)
1863 0 : sp1 = op%samples+1
1864 0 : IF ( op%opt_histo .GT. 0 ) THEN
1865 0 : op%occup_histo_time= 0.d0
1866 0 : op%occupconfig= 0.d0
1867 0 : op%suscep= 0.d0
1868 0 : op%chi= 0.d0
1869 0 : op%chicharge= 0.d0
1870 0 : op%ntot= 0.d0
1871 : END IF
1872 :
1873 0 : old_percent = 0
1874 :
1875 0 : MALLOC(updated_swap,(1:flavors))
1876 0 : updated_swap(:) = .FALSE.
1877 :
1878 0 : NRJ_old1 = 0.d0
1879 0 : NRJ_old2 = 0.d0
1880 0 : NRJ_new = 0.d0
1881 :
1882 0 : MALLOC(gtmp_new,(1,1))
1883 0 : gtmp_new = 0.d0
1884 0 : MALLOC(gtmp_old1,(1,1))
1885 0 : gtmp_old1 = 0.d0
1886 0 : MALLOC(gtmp_old2,(1,1))
1887 0 : gtmp_old2 = 0.d0
1888 :
1889 : !PROBLEM eos_gnu_13.2_mpich . %endDensity was introduced throughout
1890 0 : endDensity = SIZE(op%density,2)
1891 : ! endDensity=op%endDensity
1892 : !ENDPROBLEM
1893 :
1894 :
1895 0 : IF ( op%opt_noise .GT. 0 ) THEN
1896 0 : FREEIF(gtmp_new)
1897 0 : MALLOC(gtmp_new,(1:sp1,1:flavors))
1898 0 : FREEIF(gtmp_old1)
1899 0 : MALLOC(gtmp_old1,(1:sp1,1:flavors))
1900 0 : FREEIF(gtmp_old2)
1901 0 : MALLOC(gtmp_old2,(1:sp1,1:flavors))
1902 : END IF
1903 :
1904 0 : IF ( op%rank .EQ. 0 ) THEN
1905 : WRITE(op%ostream, '(1x,103A)') &
1906 0 : "|----------------------------------------------------------------------------------------------------|"
1907 0 : WRITE(op%ostream,'(1x,A)', ADVANCE="NO") "|"
1908 : END IF
1909 :
1910 0 : total = DBLE(itotal)
1911 : !write(std_out,*) "itotal",itotal
1912 0 : indDensity = 1
1913 : !write(std_out,*) "op%stats",op%stats
1914 0 : DO isweep = 1, itotal
1915 : !ii if(op%prtopt==1) write(std_out,*) "======== Isweep = ",isweep
1916 : !updated_seg=.FALSE.
1917 0 : DO iflavor = 1, flavors
1918 : ! if(isweep==itotal) write(std_out,*) " Iflavor = ",iflavor,op%Impurity%Particles(iflavor)%tail
1919 : !ii if(op%prtopt==1) write(std_out,*) " ===Iflavor = ",iflavor
1920 0 : op%Impurity%activeFlavor=iflavor
1921 0 : op%Bath%activeFlavor=iflavor ; op%Bath%MAddFlag= .FALSE. ; op%Bath%MRemoveFlag = .FALSE.
1922 :
1923 : !write(std_out,*) "before tryaddremove"
1924 :
1925 : ! For iflavor, Try a move
1926 : !==========================
1927 0 : CALL CtqmcoffdiagComplex_tryAddRemove(op,updated_seg)
1928 : !sui!write(std_out,*) "after tryaddremove",updated_seg
1929 :
1930 0 : updated = updated_seg .OR. updated_swap(iflavor).OR.(isweep==1)
1931 0 : updated_swap(iflavor) = .FALSE.
1932 0 : if ( op%opt_nondiag >0 ) iflavor_d=0
1933 0 : if ( op%opt_nondiag==0 ) iflavor_d=iflavor
1934 0 : CALL GreenHyboffdiagComplex_measHybrid(op%Greens, op%Bath%M, op%Impurity%Particles, updated,op%signvalue,op%phasevalue,iflavor_d)
1935 :
1936 0 : CALL CtqmcoffdiagComplex_measN (op, iflavor, updated)
1937 0 : IF ( op%opt_analysis .EQ. 1 ) &
1938 0 : CALL CtqmcoffdiagComplex_measCorrelation (op, iflavor)
1939 0 : IF ( op%opt_order .GT. 0 ) &
1940 0 : CALL CtqmcoffdiagComplex_measPerturbation(op, iflavor)
1941 : END DO
1942 : !CALL GreenHyboffdiag_measHybrid(op%Greens, op%Bath%M, op%Impurity%Particles, updated,op%signvalue,iflavor_d)
1943 : !DO iflavor = 1,flavors
1944 : ! CALL CtqmcoffdiagComplex_measN (op, iflavor, updated)
1945 : !END DO
1946 :
1947 0 : IF ( MOD(isweep,modGlobalMove) .EQ. 0 ) THEN
1948 : ! !sui!write(std_out,*) "isweep,modGlobalMove,inside",isweep,modGlobalMove
1949 0 : CALL CtqmcoffdiagComplex_trySwap(op,swapUpdate1, swapUpdate2)
1950 : ! !write(std_out,*) "no global move yet for non diag hybridization"
1951 0 : IF ( swapUpdate1 .NE. 0 .AND. swapUpdate2 .NE. 0 ) THEN
1952 0 : updated_swap(swapUpdate1) = .TRUE.
1953 0 : updated_swap(swapUpdate2) = .TRUE.
1954 : END IF
1955 : END IF
1956 :
1957 0 : IF ( MOD(isweep,measurements) .EQ. 0 ) THEN ! default is always
1958 0 : CALL ImpurityOperator_measDE(op%Impurity,op%measDE)
1959 0 : IF ( op%opt_spectra .GE. 1 .AND. MOD(isweep,measurements*op%opt_spectra) .EQ. 0 ) THEN
1960 0 : op%density(1:flavors,indDensity) = real(op%measN(3,1:flavors))
1961 0 : indDensity = indDensity+1
1962 : END IF
1963 : END IF
1964 :
1965 0 : IF ( MOD(isweep,measurements) .EQ. 0 ) THEN
1966 0 : IF ( op%opt_histo .GT. 0 ) THEN
1967 : CALL ImpurityOperator_occup_histo_time(op%Impurity,op%occup_histo_time,op%occupconfig,op%suscep,op%samples,op%chi,&
1968 0 : & op%chicharge,op%ntot,op%opt_histo,op%nspinor)
1969 : END IF
1970 : ENDIF
1971 :
1972 0 : IF ( MOD(isweep, modNoise1) .EQ. 0 ) THEN
1973 : !modNext = isweep + modNoise2
1974 0 : NRJ_new = op%measDE(1,1)
1975 0 : CALL Vector_pushBack(op%measNoise(1),NRJ_new - NRJ_old1)
1976 0 : NRJ_old1 = NRJ_new
1977 :
1978 : !! Try to limit accumulation error
1979 0 : CALL ImpurityOperator_cleanOverlaps(op%Impurity)
1980 :
1981 0 : IF ( op%opt_noise .EQ. 1 ) THEN
1982 0 : DO ifl1 = 1, flavors
1983 0 : DO ind = 1, op%Greens%map(ifl1,ifl1)%tail
1984 0 : itau = op%Greens%map(ifl1,ifl1)%listINT(ind)
1985 : gtmp_new(itau,ifl1) = op%Greens%oper(itau,ifl1,ifl1) &
1986 0 : +op%Greens%map(ifl1,ifl1)%listDBLE(ind)*DBLE(op%Greens%factor)
1987 : END DO
1988 0 : DO itau = 1, sp1
1989 0 : CALL Vector_pushBack(op%measNoiseG(itau,ifl1,1), gtmp_new(itau,ifl1) - gtmp_old1(itau,ifl1))
1990 0 : gtmp_old1(itau,ifl1) = gtmp_new(itau,ifl1)
1991 : END DO
1992 : END DO
1993 : END IF
1994 : END IF
1995 :
1996 0 : IF ( MOD(isweep,modNoise2) .EQ. 0 ) THEN
1997 0 : NRJ_new = op%measDE(1,1)
1998 0 : CALL Vector_pushBack(op%measNoise(2),NRJ_new - NRJ_old2)
1999 0 : NRJ_old2 = NRJ_new
2000 0 : IF ( op%opt_noise .EQ. 1 ) THEN
2001 0 : DO ifl1 = 1, flavors
2002 0 : DO ind = 1, op%Greens%map(ifl1,ifl1)%tail
2003 0 : itau = op%Greens%map(ifl1,ifl1)%listINT(ind)
2004 : gtmp_new(itau,ifl1) = op%Greens%oper(itau,ifl1,ifl1) &
2005 0 : +op%Greens%map(ifl1,ifl1)%listDBLE(ind)*op%Greens%factor
2006 : END DO
2007 0 : DO itau = 1, sp1
2008 0 : CALL Vector_pushBack(op%measNoiseG(itau,ifl1,2), gtmp_new(itau,ifl1) - gtmp_old2(itau,ifl1))
2009 0 : gtmp_old2(itau,ifl1) = gtmp_new(itau,ifl1)
2010 : END DO
2011 : END DO
2012 : END IF
2013 :
2014 0 : IF ( op%rank .EQ. 0 ) THEN
2015 0 : new_percent = CEILING(DBLE(isweep)*100.d0/DBLE(itotal))
2016 0 : DO ipercent = old_percent+1, new_percent
2017 0 : WRITE(op%ostream,'(A)',ADVANCE="NO") "-"
2018 : END DO
2019 : old_percent = new_percent
2020 : END IF
2021 : END IF
2022 :
2023 0 : IF ( op%opt_movie .EQ. 1 ) THEN
2024 0 : WRITE(ilatex,'(A11,I9)') "%iteration ", isweep
2025 0 : CALL ImpurityOperator_printLatex(op%Impurity,ilatex,isweep)
2026 : END IF
2027 :
2028 : END DO
2029 :
2030 0 : IF ( op%rank .EQ. 0 ) THEN
2031 0 : DO ipercent = old_percent+1, 100
2032 0 : WRITE(op%ostream,'(A)',ADVANCE="NO") "-"
2033 : END DO
2034 0 : WRITE(op%ostream,'(A)') "|"
2035 : END IF
2036 :
2037 0 : FREE(gtmp_new)
2038 0 : FREE(gtmp_old1)
2039 0 : FREE(gtmp_old2)
2040 0 : FREE(updated_swap)
2041 :
2042 0 : IF ( op%opt_spectra .GE. 1 .AND. itotal .EQ. op%sweeps ) THEN
2043 0 : IF ( endDensity .NE. indDensity-1 ) THEN
2044 0 : op%density(:,endDensity) = -1.d0
2045 : END IF
2046 : END IF
2047 :
2048 0 : CALL CPU_TIME(cpu_time2)
2049 :
2050 0 : op%runTime = (cpu_time2 - cpu_time1)*1.05d0 ! facteur arbitraire de correction
2051 0 : END SUBROUTINE CtqmcoffdiagComplex_loop
2052 : !!***
2053 :
2054 :
2055 :
2056 : !!****f* ABINIT/m_CtqmcoffdiagComplex/CtqmcoffdiagComplex_tryAddRemove
2057 : !! NAME
2058 : !! CtqmcoffdiagComplex_tryAddRemove
2059 : !!
2060 : !! FUNCTION
2061 : !! Try to add or remove a segment and an anti-segment
2062 : !!
2063 : !! COPYRIGHT
2064 : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
2065 : !! This file is distributed under the terms of the
2066 : !! GNU General Public License, see ~abinit/COPYING
2067 : !! or http://www.gnu.org/copyleft/gpl.txt .
2068 : !!
2069 : !! INPUTS
2070 : !! op=ctqmc
2071 : !!
2072 : !! OUTPUT
2073 : !! updated=something changed
2074 : !!
2075 : !! SIDE EFFECTS
2076 : !!
2077 : !! NOTES
2078 : !!
2079 : !! SOURCE
2080 :
2081 0 : SUBROUTINE CtqmcoffdiagComplex_tryAddRemove(op,updated)
2082 :
2083 : !Arguments ------------------------------------
2084 : TYPE(CtqmcoffdiagComplex) , INTENT(INOUT) :: op
2085 : ! TYPE(BathOperatoroffdiag) , INTENT(INOUT) :: Bath
2086 : ! TYPE(ImpurityOperator), INTENT(INOUT) :: Impurity
2087 : LOGICAL , INTENT( OUT) :: updated
2088 : !Local variables ------------------------------
2089 : INTEGER :: position
2090 : INTEGER , DIMENSION(1:2) :: nature ! -2 for antiseg and 1 for seg
2091 : INTEGER :: i! -2 for antiseg and 1 for seg
2092 : !INTEGER :: it,it1 !ii,
2093 : DOUBLE PRECISION :: action
2094 : DOUBLE PRECISION :: beta
2095 : DOUBLE PRECISION :: time1
2096 : DOUBLE PRECISION :: time2
2097 : DOUBLE PRECISION :: time_avail
2098 : COMPLEX(KIND=8) :: det_ratio
2099 : DOUBLE PRECISION :: sign_det_ratio
2100 : DOUBLE PRECISION :: overlap
2101 : DOUBLE PRECISION :: length
2102 : DOUBLE PRECISION :: signe
2103 : DOUBLE PRECISION :: tail
2104 : INTEGER :: tailint
2105 : DOUBLE PRECISION :: signdet, signdetprev
2106 : DOUBLE PRECISION, DIMENSION(1:2) :: CdagC_1
2107 :
2108 0 : IF ( .NOT. op%set ) &
2109 0 : CALL ERROR("CtqmcoffdiagComplex_trySegment : QMC not set ")
2110 :
2111 : !write(std_out,*) " TryAddRemove start"
2112 0 : nature(1) = CTQMC_SEGME
2113 0 : nature(2) = CTQMC_ANTIS
2114 0 : beta = op%beta
2115 :
2116 0 : updated = .FALSE.
2117 : !tail is the number of segments for a given flavor
2118 0 : tailint = (op%Impurity%particles(op%Impurity%activeFlavor)%tail)
2119 0 : tail = DBLE(tailint)
2120 :
2121 :
2122 : !=====================================
2123 : ! First choose segment or antisegment
2124 : !=====================================
2125 0 : DO i = 1, 2
2126 0 : signe = SIGN(1.d0,DBLE(nature(i)))
2127 : ! ----- 1: segment signe= 1 ( CTQMC_SEGME = 1 )
2128 : ! ----- 2: antisegment signe=-1 ( CTQMC_ANTIS = -2 )
2129 :
2130 0 : tailint = (op%Impurity%particles(op%Impurity%activeFlavor)%tail)
2131 :
2132 : !CALL RANDOM_NUMBER(action)
2133 0 : CALL OurRng(op%seed,action)
2134 :
2135 : !==========================
2136 : ! Add segment/antisegment
2137 : !==========================
2138 0 : IF ( action .LT. .5d0 ) THEN ! Add a segment or antisegment
2139 :
2140 : ! Select time1 (>0) in [0,beta]
2141 : !==============================
2142 : !CALL RANDOM_NUMBER(time1)
2143 0 : CALL OurRng(op%seed,time1)
2144 0 : time1 = time1 * beta
2145 :
2146 : ! time_avail is the distance between between time1 and
2147 : ! - the next start of a segment for a segment addition
2148 : ! - the next end of a segment for an antisegment addition
2149 : ! ImpurityOperator_getAvailableTime > 0 for a segment (signe>0) -> time_avail>0
2150 : ! ImpurityOperator_getAvailableTime < 0 for an antisegment (signe<0) -> time_avail>0
2151 : !====================================================================
2152 0 : time_avail = ImpurityOperator_getAvailableTime(op%Impurity,time1,position) * signe
2153 :
2154 0 : IF ( time_avail .GT. 0.d0 ) THEN
2155 :
2156 : ! Time2 is the length of the proposed new (anti)segment
2157 : !=======================================================
2158 0 : CALL OurRng(op%seed,time2)
2159 0 : IF ( time2 .EQ. 0.d0 ) CALL OurRng(op%seed,time2) ! Prevent null segment
2160 :
2161 : ! Now time2 is the time at the end of the proposed new (anti) segment
2162 : ! time2 > time1
2163 : !====================================================================
2164 0 : time2 = time1 + time2 * time_avail
2165 :
2166 : ! CdagC_1 gives the stard/end times for the proposed new segment/antisegment
2167 : ! CdagC1(C_) can be above beta.
2168 : ! For a segment CdagC_1(Cdag_) = time1 < CdagC_1(C_) = time2, l=time2-time1 > 0
2169 : ! For a anti segment CdagC_1(Cdag_) = time2 > CdagC_1(C_) = time1, l=time1-time2 < 0
2170 : ! time2 can be above beta and thus for a segment CdagC_1(C_ ) > beta
2171 : ! time2 can be above beta and thus for an antisegment CdagC_1(Cdag_) > beta
2172 : ! length > 0 for segment
2173 : ! length < 0 for antisegment
2174 : !====================================================================================
2175 0 : CdagC_1(Cdag_) = ((1.d0+signe)*time1+(1.d0-signe)*time2)*0.5d0
2176 0 : CdagC_1(C_ ) = ((1.d0+signe)*time2+(1.d0-signe)*time1)*0.5d0
2177 0 : length = CdagC_1(C_ ) - CdagC_1(Cdag_)
2178 :
2179 : ! ----- Computes the determinant ratio
2180 0 : det_ratio = BathOperatoroffdiagComplex_getDetAdd(op%Bath,CdagC_1,position,op%Impurity%particles)
2181 :
2182 : ! ----- Computes the overlap
2183 0 : overlap = ImpurityOperator_getNewOverlap(op%Impurity,CdagC_1)
2184 0 : signdetprev = ImpurityOperator_getsign(op%Impurity, time2, i, action, position)
2185 : !write(std_out,*)"det_ratio",det_ratio,"signdetprev",signdetprev
2186 :
2187 : !Remove negative sign of det_ration
2188 0 : det_ratio=det_ratio*signdetprev
2189 :
2190 0 : IF ( real(det_ratio) .LT. 0.d0 ) THEN
2191 0 : det_ratio = - det_ratio
2192 0 : sign_det_ratio=-1
2193 0 : op%stats(nature(i)+CTQMC_DETSI) = op%stats(nature(i)+CTQMC_DETSI) + 1.d0
2194 : ELSE
2195 : sign_det_ratio=1
2196 : END IF
2197 :
2198 0 : CALL OurRng(op%seed,time1)
2199 : !write(std_out,*) " .................",(time1 * (tail + 1.d0 )),beta * time_avail * abs(det_ratio) * DEXP(real(op%mu(op%Impurity%activeFlavor))*length + overlap)
2200 : !write(std_out,*) " .................",beta , time_avail , op%mu(op%Impurity%activeFlavor),op%Impurity%activeFlavor
2201 :
2202 0 : IF ( (time1 * (tail + 1.d0 )) &
2203 : .LT. (beta * time_avail * abs(det_ratio) * DEXP(real(op%mu(op%Impurity%activeFlavor))*length + overlap) ) ) THEN
2204 0 : CALL ImpurityOperator_add(op%Impurity,CdagC_1,position)
2205 : ! write(*,*) "after "
2206 : ! CALL ListCdagC_print(op%Impurity%particles(op%Impurity%activeFlavor),6)
2207 0 : CALL BathOperatoroffdiagComplex_setMAdd(op%bath,op%Impurity%particles)
2208 0 : op%stats(nature(i)+CTQMC_ADDED) = op%stats(nature(i)+CTQMC_ADDED) + 1.d0
2209 0 : updated = .TRUE. .OR. updated
2210 0 : tail = tail + 1.d0
2211 0 : tailint = tailint + 1
2212 :
2213 0 : IF ( sign_det_ratio .LT. 0.d0 ) op%signvalue=-op%signvalue
2214 0 : op%phasevalue=op%phasevalue*det_ratio/abs(det_ratio)*sign_det_ratio
2215 : END IF
2216 :
2217 : END IF ! time_avail > 0
2218 :
2219 : !========================================
2220 : ! Remove segment/antisegment
2221 : !========================================
2222 : ELSE ! Remove a segment among the segment of the flavor activeflavor
2223 : !ii if(op%prtopt==1) write(6,*) " =try: Segment removed of type",i
2224 0 : IF ( tail .GT. 0.d0 ) THEN
2225 0 : CALL OurRng(op%seed,time1)
2226 0 : position = INT(((time1 * tail) + 1.d0) * signe )
2227 : !prt!if(op%prtopt==1) write(6,*) " position",position
2228 0 : time_avail = ImpurityOperator_getAvailedTime(op%Impurity,position)
2229 0 : det_ratio = BathOperatoroffdiagComplex_getDetRemove(op%Bath,position)
2230 : !write(6,*) " det_ratio", det_ratio
2231 0 : CdagC_1 = ImpurityOperator_getSegment(op%Impurity,position)
2232 : ! length = CdagC_length(CdagC_1)
2233 0 : length = CdagC_1(C_) - CdagC_1(Cdag_)
2234 : !write(6,*) " length ", length
2235 0 : overlap = ImpurityOperator_getNewOverlap(op%Impurity,CdagC_1)
2236 : !write(6,*) " overlap ", overlap
2237 0 : CALL OurRng(op%seed,time1)
2238 : !write(6,*) " Random ",time1
2239 0 : signdetprev = ImpurityOperator_getsign(op%Impurity, time2, i, action, position)
2240 0 : det_ratio=det_ratio*signdetprev
2241 0 : signdet=1.d0
2242 :
2243 0 : IF ( real(det_ratio) .LT. 0.d0 ) THEN
2244 : !sui!write(6,*) " NEGATIVE DET",det_ratio,signdetprev
2245 0 : det_ratio = -det_ratio
2246 0 : sign_det_ratio=-1
2247 0 : op%stats(nature(i)+CTQMC_DETSI) = op%stats(nature(i)+CTQMC_DETSI) + 1.d0
2248 0 : signdet=-1.d0
2249 : ELSE
2250 : sign_det_ratio=1
2251 : END IF
2252 :
2253 0 : IF ( (time1 * beta * time_avail * DEXP(real(op%mu(op%Impurity%activeFlavor))*length+overlap)) &
2254 : .LT. (tail * abs(det_ratio) ) ) THEN
2255 0 : CALL ImpurityOperator_remove(op%Impurity,position)
2256 0 : CALL BathOperatoroffdiagComplex_setMRemove(op%Bath,op%Impurity%particles)
2257 : !op%seg_removed = op%seg_removed + 1.d0
2258 0 : op%stats(nature(i)+CTQMC_REMOV) = op%stats(nature(i)+CTQMC_REMOV) + 1.d0
2259 0 : updated = .TRUE. .OR. updated
2260 0 : tail = tail -1.d0
2261 0 : tailint = tailint -1
2262 :
2263 0 : IF ( sign_det_ratio .LT. 0.d0 ) op%signvalue=-op%signvalue
2264 0 : op%phasevalue=op%phasevalue*det_ratio/abs(det_ratio)*sign_det_ratio
2265 : END IF
2266 : END IF !tail > O
2267 : END IF ! action
2268 : !========================================
2269 : ! End Add/Remove Antisegment
2270 : !========================================
2271 : END DO
2272 0 : END SUBROUTINE CtqmcoffdiagComplex_tryAddRemove
2273 : !!***
2274 :
2275 : !!SUBROUTINE CtqmcoffdiagComplex_trySegment(op,updated)
2276 : !!! TYPE(BathOperatoroffdiag) , INTENT(INOUT) :: Bath
2277 : !!! TYPE(ImpurityOperator), INTENT(INOUT) :: Impurity
2278 : !! LOGICAL , INTENT(INOUT) :: updated
2279 : !! INTEGER :: position
2280 : !! DOUBLE PRECISION :: action
2281 : !! DOUBLE PRECISION :: beta
2282 : !! DOUBLE PRECISION :: time1
2283 : !! DOUBLE PRECISION :: time2
2284 : !! DOUBLE PRECISION :: time_avail
2285 : !! DOUBLE PRECISION :: det_ratio
2286 : !! DOUBLE PRECISION :: overlap
2287 : !! DOUBLE PRECISION :: length
2288 : !! DOUBLE PRECISION :: tail
2289 : !! DOUBLE PRECISION, DIMENSION(1:2) :: CdagC_1
2290 : !!
2291 : !! IF ( .NOT. op%set ) &
2292 : !! CALL ERROR("CtqmcoffdiagComplex_trySegment : QMC not set ")
2293 : !!
2294 : !! beta = op%beta
2295 : !! tail = DBLE(op%Impurity%particles(op%Impurity%activeFlavor)%tail)
2296 : !!
2297 : !! CALL RANDOM_NUMBER(action)
2298 : !!
2299 : !! updated = .FALSE.
2300 : !!
2301 : !! IF ( action .LT. .5d0 ) THEN ! Ajout de segment
2302 : !! CALL RANDOM_NUMBER(time1)
2303 : !! time1 = time1 * beta
2304 : !! time_avail = ImpurityOperator_getAvailableTime(op%Impurity,time1,position)
2305 : !! IF ( time_avail .GT. 0.d0 ) THEN
2306 : !! CALL RANDOM_NUMBER(time2)
2307 : !! time2 = time1 + time2 * time_avail
2308 : !!! CALL CdagC_init(CdagC_1,time1,time2)
2309 : !! CdagC_1(Cdag_) = time1
2310 : !! CdagC_1(C_ ) = time2
2311 : !!! length = CdagC_length(CdagC_1)
2312 : !! length = time2 - time1
2313 : !! det_ratio = BathOperatoroffdiag_getDetAdd(op%Bath,CdagC_1,position,op%Impurity%particles(op%Impurity%activeFlavor))
2314 : !! overlap = ImpurityOperator_getNewOverlap(op%Impurity,CdagC_1)
2315 : !! CALL RANDOM_NUMBER(time1)
2316 : !! IF ( det_ratio .LT. 0.d0 ) THEN
2317 : !! det_ratio = -det_ratio
2318 : !! op%stats(CTQMC_SEGME+CTQMC_DETSI) = op%stats(CTQMC_SEGME+CTQMC_DETSI) + SIGN(1.d0,det_ratio)
2319 : !! END IF
2320 : !! IF ( (time1 * (tail + 1.d0 )) &
2321 : !! .LT. (beta * time_avail * det_ratio * DEXP(op%mu(op%Impurity%activeFlavor)*length + overlap) ) ) THEN
2322 : !! write(*,*) position
2323 : !! CALL ImpurityOperator_add(op%Impurity,CdagC_1,position)
2324 : !! CALL BathOperatoroffdiag_setMAdd(op%bath,op%Impurity%particles(op%Impurity%activeFlavor))
2325 : !! op%stats(CTQMC_SEGME+CTQMC_ADDED) = op%stats(CTQMC_SEGME+CTQMC_ADDED) + 1.d0
2326 : !! updated = .TRUE.
2327 : !! END IF
2328 : !! END IF
2329 : !!
2330 : !! ELSE ! Supprimer un segment
2331 : !! IF ( tail .GT. 0.d0 ) THEN
2332 : !! CALL RANDOM_NUMBER(time1)
2333 : !! position = INT(time1 * tail) + 1
2334 : !! time_avail = ImpurityOperator_getAvailedTime(op%Impurity,position)
2335 : !! det_ratio = BathOperatoroffdiag_getDetRemove(op%Bath,position)
2336 : !! CdagC_1 = ImpurityOperator_getSegment(op%Impurity,position)
2337 : !!! length = CdagC_length(CdagC_1)
2338 : !! length = CdagC_1(C_) - CdagC_1(Cdag_)
2339 : !! overlap = ImpurityOperator_getNewOverlap(op%Impurity,CdagC_1)
2340 : !! CALL RANDOM_NUMBER(time1)
2341 : !! IF ( det_ratio .LT. 0.d0 ) THEN
2342 : !! det_ratio = -det_ratio
2343 : !!! op%seg_sign = op%seg_sign + 1.d0
2344 : !! op%stats(CTQMC_SEGME+CTQMC_DETSI) = op%stats(CTQMC_SEGME+CTQMC_DETSI) + SIGN(1.d0,det_ratio)
2345 : !! END IF
2346 : !! IF ( (time1 * beta * time_avail * DEXP(op%mu(op%Impurity%activeFlavor)*length+overlap)) &
2347 : !! .LT. (tail * det_ratio ) ) THEN
2348 : !! write(*,*) position
2349 : !! CALL ImpurityOperator_remove(op%Impurity,position)
2350 : !! CALL BathOperatoroffdiag_setMRemove(op%Bath,op%Impurity%particles(op%Impurity%activeFlavor))
2351 : !! !op%seg_removed = op%seg_removed + 1.d0
2352 : !! op%stats(CTQMC_SEGME+CTQMC_REMOV) = op%stats(CTQMC_SEGME+CTQMC_REMOV) + 1.d0
2353 : !! updated = .TRUE.
2354 : !! END IF
2355 : !! END IF
2356 : !! END IF
2357 : !!END SUBROUTINE CtqmcoffdiagComplex_trySegment
2358 : !!***
2359 : !!
2360 : !!SUBROUTINE CtqmcoffdiagComplex_tryAntiSeg(op, updated)
2361 : !!! TYPE(BathOperatoroffdiag) , INTENT(INOUT) :: Bath
2362 : !!! TYPE(ImpurityOperator), INTENT(INOUT) :: Impurity
2363 : !! LOGICAL , INTENT(INOUT) :: updated
2364 : !! INTEGER :: position
2365 : !! DOUBLE PRECISION :: action
2366 : !! DOUBLE PRECISION :: beta
2367 : !! DOUBLE PRECISION :: time1
2368 : !! DOUBLE PRECISION :: time2
2369 : !! DOUBLE PRECISION :: time_avail
2370 : !! DOUBLE PRECISION :: det_ratio
2371 : !! DOUBLE PRECISION :: overlap
2372 : !! DOUBLE PRECISION :: length
2373 : !! DOUBLE PRECISION :: tail
2374 : !! DOUBLE PRECISION, DIMENSION(1:2) :: CdagC_1
2375 : !!
2376 : !! IF ( .NOT. op%set ) &
2377 : !! CALL ERROR("CtqmcoffdiagComplex_trySegment : QMC not set ")
2378 : !!
2379 : !! beta = op%beta
2380 : !! tail = DBLE(op%Impurity%particles(op%Impurity%activeFlavor)%tail)
2381 : !!
2382 : !! CALL RANDOM_NUMBER(action)
2383 : !!
2384 : !! updated = .FALSE.
2385 : !!
2386 : !! IF ( action .LT. .5d0 ) THEN ! Ajout d'un antiseg
2387 : !! CALL RANDOM_NUMBER(time1)
2388 : !! time1 = time1 * beta
2389 : !! time_avail = ImpurityOperator_getAvailableTime(op%Impurity,time1,position)
2390 : !! IF ( time_avail .LT. 0.d0 ) THEN
2391 : !! CALL RANDOM_NUMBER(time2)
2392 : !! time2 = time1 - time2 * time_avail
2393 : !!! CALL CdagC_init(CdagC_1,time2,time1)
2394 : !! CdagC_1(Cdag_) = time2
2395 : !! CdagC_1(C_ ) = time1
2396 : !!! length = CdagC_length(CdagC_1) ! /!\ length is negative
2397 : !! length = time1 - time2
2398 : !! det_ratio = BathOperatoroffdiag_getDetAdd(op%Bath,CdagC_1,position,op%Impurity%particles(op%Impurity%activeFlavor))
2399 : !! overlap = ImpurityOperator_getNewOverlap(op%Impurity,CdagC_1) ! OK
2400 : !! CALL RANDOM_NUMBER(time1)
2401 : !! IF ( det_ratio .LT. 0.d0 ) THEN
2402 : !! det_ratio = -det_ratio
2403 : !!! op%anti_sign = op%anti_sign + 1.d0
2404 : !! op%stats(CTQMC_ANTIS+CTQMC_DETSI) = op%stats(CTQMC_ANTIS+CTQMC_DETSI) + SIGN(1.d0,det_ratio)
2405 : !! END IF
2406 : !! IF ( (time1 * (tail + 1.d0 )) &
2407 : !! .LT. (beta * ABS(time_avail) * det_ratio * DEXP(op%mu(op%Impurity%activeFlavor)*length + overlap) ) ) THEN
2408 : !! CALL ImpurityOperator_add(op%Impurity,CdagC_1,position)
2409 : !! !write(*,*) position
2410 : !! CALL BathOperatoroffdiag_setMAdd(op%bath,op%Impurity%particles(op%Impurity%activeFlavor))
2411 : !! !op%anti_added = op%anti_added + 1.d0
2412 : !! op%stats(CTQMC_ANTIS+CTQMC_ADDED) = op%stats(CTQMC_ANTIS+CTQMC_ADDED) + 1.d0
2413 : !! updated = .TRUE.
2414 : !! END IF
2415 : !! END IF
2416 : !! ELSE ! Supprimer un antiseg
2417 : !! IF ( tail .GT. 0.d0 ) THEN
2418 : !! CALL RANDOM_NUMBER(time1)
2419 : !! position = -(INT(time1 * tail) + 1)
2420 : !! time_avail = ImpurityOperator_getAvailedTime(op%Impurity,position)!OK
2421 : !! det_ratio = BathOperatoroffdiag_getDetRemove(op%Bath,position)!OK
2422 : !! CdagC_1 = ImpurityOperator_getSegment(op%Impurity,position)!OK
2423 : !!! length = CdagC_length(CdagC_1) ! /!\ length is negative
2424 : !! length = CdagC_1(C_) - CdagC_1(Cdag_)
2425 : !! overlap = ImpurityOperator_getNewOverlap(op%Impurity,CdagC_1) !OK
2426 : !! CALL RANDOM_NUMBER(time1)
2427 : !! IF ( det_ratio .LT. 0.d0 ) THEN
2428 : !! det_ratio = -det_ratio
2429 : !!! op%anti_sign = op%anti_sign + 1.d0
2430 : !! op%stats(CTQMC_ANTIS+CTQMC_DETSI) = op%stats(CTQMC_ANTIS+CTQMC_DETSI) + SIGN(1.d0,det_ratio)
2431 : !! END IF
2432 : !! IF ( (time1 * beta * time_avail * DEXP(op%mu(op%Impurity%activeFlavor)*length+overlap)) &
2433 : !! .LT. (tail * det_ratio ) ) THEN
2434 : !! CALL ImpurityOperator_remove(op%Impurity,position)
2435 : !! !write(*,*) position
2436 : !! CALL BathOperatoroffdiag_setMRemove(op%Bath,op%Impurity%particles(op%Impurity%activeFlavor))
2437 : !! !op%anti_removed = op%anti_removed + 1.d0
2438 : !! op%stats(CTQMC_ANTIS+CTQMC_REMOV) = op%stats(CTQMC_ANTIS+CTQMC_REMOV) + 1.d0
2439 : !! updated = .TRUE.
2440 : !! END IF
2441 : !! END IF
2442 : !! END IF
2443 : !!END SUBROUTINE CtqmcoffdiagComplex_tryAntiSeg
2444 : !!***
2445 :
2446 : !!****f* ABINIT/m_CtqmcoffdiagComplex/CtqmcoffdiagComplex_trySwap
2447 : !! NAME
2448 : !! CtqmcoffdiagComplex_trySwap
2449 : !!
2450 : !! FUNCTION
2451 : !! try a global move (swap to flavors)
2452 : !!
2453 : !! COPYRIGHT
2454 : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
2455 : !! This file is distributed under the terms of the
2456 : !! GNU General Public License, see ~abinit/COPYING
2457 : !! or http://www.gnu.org/copyleft/gpl.txt .
2458 : !!
2459 : !! INPUTS
2460 : !! op=ctqmc
2461 : !!
2462 : !! OUTPUT
2463 : !! flav_i=first flavor swaped
2464 : !! flav_j=second flavor swaped
2465 : !!
2466 : !! SIDE EFFECTS
2467 : !!
2468 : !! NOTES
2469 : !!
2470 : !! SOURCE
2471 :
2472 0 : SUBROUTINE CtqmcoffdiagComplex_trySwap(op,flav_i,flav_j)
2473 :
2474 : !Arguments ------------------------------------
2475 : TYPE(CtqmcoffdiagComplex) , INTENT(INOUT) :: op
2476 : ! TYPE(BathOperatoroffdiag) , INTENT(INOUT) :: Bath
2477 : ! TYPE(ImpurityOperator), INTENT(INOUT) :: Impurity
2478 : INTEGER , INTENT( OUT) :: flav_i
2479 : INTEGER , INTENT( OUT) :: flav_j
2480 : !Local variables ------------------------------
2481 : INTEGER :: flavor_i
2482 : INTEGER :: flavor_j !,ii,it,it1 !,iflavor
2483 : DOUBLE PRECISION :: rnd
2484 : DOUBLE PRECISION :: lengthi
2485 : DOUBLE PRECISION :: lengthj
2486 : DOUBLE PRECISION :: overlapic1
2487 : DOUBLE PRECISION :: overlapjc1
2488 : DOUBLE PRECISION :: overlapic2
2489 : DOUBLE PRECISION :: overlapjc2
2490 : !DOUBLE PRECISION :: detic1
2491 : !DOUBLE PRECISION :: detjc1
2492 : !DOUBLE PRECISION :: detic2
2493 : !DOUBLE PRECISION :: detjc2
2494 : DOUBLE PRECISION :: det_ratio,detnew,detold
2495 : DOUBLE PRECISION :: local_ratio
2496 : ! TYPE(BathOperatoroffdiag) :: Bathnew
2497 :
2498 :
2499 : !CALL RANDOM_NUMBER(rnd)
2500 0 : CALL OurRng(op%seed,rnd)
2501 0 : flavor_i = NINT(rnd*DBLE(op%flavors-1.d0))+1
2502 : !CALL RANDOM_NUMBER(rnd)
2503 0 : CALL OurRng(op%seed,rnd)
2504 0 : flavor_j = NINT(rnd*DBLE(op%flavors-1.d0))+1
2505 : !ii write(6,'(a,2i4)') "--------------- new swap --------------------------------",flavor_i,flavor_j
2506 :
2507 0 : flav_i = 0
2508 0 : flav_j = 0
2509 : !ii do iflavor=1,op%flavors
2510 : !ii write(6,*) "BEFORE GMOVE For flavor", iflavor,"size is",op%Impurity%particles(iflavor)%tail," and Conf is :"
2511 : !ii do ii=1, op%Impurity%Particles(iflavor)%tail
2512 : !ii write(6,'(i4,100f12.3)') ii, op%Impurity%Particles(iflavor)%list(ii,1),&
2513 : !ii & op%Impurity%Particles(iflavor)%list(ii,2)
2514 : !ii enddo
2515 : !ii enddo
2516 : !ii write(6,*) " = M Matrix"
2517 : !ii write(6,'(a,2x,100(i12))') "Flavor=",((iflavor,it=1,op%Impurity%particles(iflavor)%tail),iflavor=1,op%flavors)
2518 : !ii write(6,'(i21,100i12)') ((it,it=1,op%Impurity%particles(iflavor)%tail),iflavor=1,op%flavors)
2519 : !ii do it=1,op%Bath%sumtails
2520 : !ii write(6,'(a,100f12.3)') " M before",(op%Bath%M%mat(it,it1),it1=1,op%Bath%sumtails)
2521 : !ii enddo
2522 :
2523 : ! todoba this part
2524 0 : IF ( flavor_i .NE. flavor_j ) THEN
2525 : !CALL BathOperatoroffdiag_init(Bathnew, op%flavors, op%samples, op%beta, 0, op%opt_nondiag)
2526 : ! On tente d'intervertir i et j
2527 : ! Configuration actuelle :
2528 :
2529 0 : op%modGlobalMove(2) = op%modGlobalMove(2)+1
2530 : ! ===========================================
2531 : ! First use M matrix to compute determinant
2532 : ! ===========================================
2533 0 : detold = BathOperatoroffdiagComplex_getDetF(op%Bath) ! use op%Bath%M
2534 :
2535 : ! ===========================================
2536 : ! Second build M_update matrix to compute determinant after.
2537 : ! ===========================================
2538 : !CALL ListCdagCoffdiag_print(particle)
2539 0 : call BathOperatoroffdiagComplex_recomputeM(op%Bath,op%impurity%particles,flavor_i,flavor_j) ! compute op%Bath%M_update
2540 0 : detnew = BathOperatoroffdiagComplex_getDetF(op%Bath,option=1) ! use op%Bath%M_update
2541 :
2542 0 : lengthi = real(ImpurityOperator_measN(op%Impurity,flavor_i))
2543 0 : lengthj = real(ImpurityOperator_measN(op%Impurity,flavor_j))
2544 0 : overlapic1 = ImpurityOperator_overlapFlavor(op%Impurity,flavor_i)
2545 0 : overlapjc1 = ImpurityOperator_overlapFlavor(op%Impurity,flavor_j)
2546 : ! lengths unchanged
2547 0 : overlapic2 = ImpurityOperator_overlapSwap(op%Impurity,flavor_i,flavor_j)
2548 0 : overlapjc2 = ImpurityOperator_overlapSwap(op%Impurity,flavor_j,flavor_i)
2549 :
2550 : ! IF ( detic1*detjc1 .EQ. detic2*detjc2 ) THEN
2551 : ! det_ratio = 1.d0
2552 : ! ELSE IF ( detic1*detjc1 .EQ. 0.d0 ) THEN
2553 : ! det_ratio = detic2*detjc2 ! evite de diviser par 0 si pas de segment
2554 : ! ELSE
2555 :
2556 0 : det_ratio = detnew/detold ! because the determinant is the determinant of F
2557 : !ii write(6,*) "det_ratio, detold,detnew",det_ratio, detold,detnew, detold/detnew
2558 :
2559 : ! END IF
2560 : local_ratio = DEXP(-overlapic2*overlapjc2+overlapic1*overlapjc1 &
2561 0 : +(lengthj-lengthi)*real((op%mu(flavor_i))-real(op%mu(flavor_j))))
2562 : !ii write(6,*) "local_ratio",local_ratio
2563 :
2564 : ! Wloc = exp(muN-Uo)
2565 : !CALL RANDOM_NUMBER(rnd)
2566 0 : CALL OurRng(op%seed,rnd)
2567 0 : IF ( rnd .LT. local_ratio*det_ratio ) THEN ! swap accepted
2568 : !ii write(6,*) " = M Matrix before swap"
2569 : !ii write(6,'(a,2x,100(i12))') "Flavor=",((iflavor,it=1,op%Impurity%particles(iflavor)%tail),iflavor=1,op%flavors)
2570 : !ii write(6,'(i21,100i12)') ((it,it=1,op%Impurity%particles(iflavor)%tail),iflavor=1,op%flavors)
2571 : !ii do it=1,op%Bath%sumtails
2572 : !ii write(6,'(a,100f12.3)') " M after ",(op%Bath%M%mat(it,it1),it1=1,op%Bath%sumtails)
2573 : !ii enddo
2574 : !ii do it=1,op%Bath%sumtails
2575 : !ii write(6,'(a,100f12.3)') " update M after ",(op%Bath%M_update%mat(it,it1),it1=1,op%Bath%sumtails)
2576 : !ii enddo
2577 : !ii write(6,*) "Gmove accepted",rnd,local_ratio*det_ratio
2578 0 : CALL ImpurityOperator_swap(op%Impurity, flavor_i,flavor_j)
2579 0 : CALL BathOperatoroffdiagComplex_swap (op%Bath , flavor_i,flavor_j) ! use op%Bath%M_update to built new op%Bath%M
2580 :
2581 0 : op%swap = op%swap + 1.d0
2582 0 : flav_i = flavor_i
2583 0 : flav_j = flavor_j
2584 : ELSE
2585 : !ii write(6,*) "Gmove refused",rnd,local_ratio*det_ratio
2586 : ! CALL WARN("Swap refused")
2587 : ! WRITE(op%ostream,'(6E24.14)') local_ratio, det_ratio, detic1, detjc1, detic2, detjc2
2588 : END IF
2589 : ! CALL BathOperatoroffdiag_destroy(Bathnew)
2590 : END IF
2591 : !ii do iflavor=1,op%flavors
2592 : !ii write(6,*) "AFTER GMOVE For flavor", iflavor,"size is",op%Impurity%particles(iflavor)%tail," and Conf is :"
2593 : !ii do ii=1, op%Impurity%Particles(iflavor)%tail
2594 : !ii write(6,'(15x,i4,100f12.3)') ii, op%Impurity%Particles(iflavor)%list(ii,1),&
2595 : !ii & op%Impurity%Particles(iflavor)%list(ii,2)
2596 : !ii enddo
2597 : !ii enddo
2598 : !ii write(6,*) " = M Matrix"
2599 : !ii write(6,'(a,2x,100(i12))') "Flavor=",((iflavor,it=1,op%Impurity%particles(iflavor)%tail),iflavor=1,op%flavors)
2600 : !ii write(6,'(i21,100i12)') ((it,it=1,op%Impurity%particles(iflavor)%tail),iflavor=1,op%flavors)
2601 : !ii do it=1,op%Bath%sumtails
2602 : !ii write(6,'(a,100f12.3)') " M after ",(op%Bath%M%mat(it,it1),it1=1,op%Bath%sumtails)
2603 : !ii enddo
2604 : !ii do it=1,op%Bath%sumtails
2605 : !ii write(6,'(a,100f12.3)') " update M after ",(op%Bath%M_update%mat(it,it1),it1=1,op%Bath%sumtails)
2606 : !ii enddo
2607 :
2608 0 : END SUBROUTINE CtqmcoffdiagComplex_trySwap
2609 : !!***
2610 :
2611 : !!****f* ABINIT/m_CtqmcoffdiagComplex/CtqmcoffdiagComplex_measN
2612 : !! NAME
2613 : !! CtqmcoffdiagComplex_measN
2614 : !!
2615 : !! FUNCTION
2616 : !! measures the number of electron
2617 : !! by taking into account the value for the move before before this one
2618 : !! with the correct weight.
2619 : !!
2620 : !! COPYRIGHT
2621 : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
2622 : !! This file is distributed under the terms of the
2623 : !! GNU General Public License, see ~abinit/COPYING
2624 : !! or http://www.gnu.org/copyleft/gpl.txt .
2625 : !!
2626 : !! INPUTS
2627 : !! op=ctqmc
2628 : !! iflavor=which flavor to measure
2629 : !! updated=something has changed since last time
2630 : !!
2631 : !! OUTPUT
2632 : !!
2633 : !! SIDE EFFECTS
2634 : !!
2635 : !! NOTES
2636 : !!
2637 : !! SOURCE
2638 :
2639 0 : SUBROUTINE CtqmcoffdiagComplex_measN(op, iflavor, updated)
2640 :
2641 : !Arguments ------------------------------------
2642 : TYPE(CtqmcoffdiagComplex) , INTENT(INOUT) :: op
2643 : !TYPE(ImpurityOperator), INTENT(IN ) :: impurity
2644 : INTEGER , INTENT(IN ) :: iflavor
2645 : LOGICAL , INTENT(IN ) :: updated
2646 :
2647 : ! IF ( .NOT. op%set ) &
2648 : ! CALL ERROR("CtqmcoffdiagComplex_measN : QMC not set ")
2649 :
2650 :
2651 0 : IF ( updated .EQV. .TRUE. ) THEN
2652 : ! --- accumulate occupations with values op%measN(3,iflavor) from the last measurements with the corresponding weight
2653 : ! --- op*measN(4,iflavor)
2654 0 : op%measN(1,iflavor) = op%measN(1,iflavor) + op%measN(3,iflavor)*op%measN(4,iflavor)
2655 : ! write(6,*) "Cllll42"
2656 :
2657 : ! --- Compute total number of new measurements
2658 0 : op%measN(2,iflavor) = op%measN(2,iflavor) + op%measN(4,iflavor)
2659 :
2660 : ! write(6,*) "Allll42"
2661 : ! --- Compute the occupation for this configuration (will be put in
2662 : ! --- op%measN(1,iflavor) at the next occurence of updated=.true.), with
2663 : ! --- the corresponding weight op%measN(4,iflavor) (we do not now it yet)
2664 0 : op%measN(3,iflavor) = ImpurityOperator_measN(op%impurity)
2665 :
2666 : ! --- set weight: as update=true, it is a new measurement , so put it to one
2667 0 : op%measN(4,iflavor) = 1.d0
2668 :
2669 : ELSE
2670 : ! --- increased the count so that at new move, we will be able to update measN(1) correctly.
2671 0 : op%measN(4,iflavor) = op%measN(4,iflavor) + 1.d0
2672 : ! write(6,*) "Bllll42"
2673 : END IF
2674 0 : END SUBROUTINE CtqmcoffdiagComplex_measN
2675 : !!***
2676 :
2677 : !#ifdef CTCtqmcoffdiagComplex_ANALYSIS
2678 : !!****f* ABINIT/m_CtqmcoffdiagComplex/CtqmcoffdiagComplex_measCorrelation
2679 : !! NAME
2680 : !! CtqmcoffdiagComplex_measCorrelation
2681 : !!
2682 : !! FUNCTION
2683 : !! measure all correlations in times for a flavor
2684 : !!
2685 : !! COPYRIGHT
2686 : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
2687 : !! This file is distributed under the terms of the
2688 : !! GNU General Public License, see ~abinit/COPYING
2689 : !! or http://www.gnu.org/copyleft/gpl.txt .
2690 : !!
2691 : !! INPUTS
2692 : !! op=ctqmc
2693 : !! iflavor=the flavor to measure
2694 : !!
2695 : !! OUTPUT
2696 : !!
2697 : !! SIDE EFFECTS
2698 : !!
2699 : !! NOTES
2700 : !!
2701 : !! SOURCE
2702 :
2703 0 : SUBROUTINE CtqmcoffdiagComplex_measCorrelation(op, iflavor)
2704 :
2705 : !Arguments ------------------------------------
2706 : TYPE(CtqmcoffdiagComplex) , INTENT(INOUT) :: op
2707 : !TYPE(ImpurityOperator), INTENT(IN ) :: impurity
2708 : INTEGER , INTENT(IN ) :: iflavor
2709 : !Local variables ------------------------------
2710 : INTEGER :: iCdag
2711 : INTEGER :: iCdagBeta
2712 : INTEGER :: iC
2713 : INTEGER :: index
2714 : INTEGER :: size
2715 : DOUBLE PRECISION :: tC
2716 : DOUBLE PRECISION :: tCdag
2717 : !DOUBLE PRECISION :: time
2718 : DOUBLE PRECISION :: inv_dt
2719 : DOUBLE PRECISION :: beta
2720 :
2721 0 : IF ( .NOT. op%set ) &
2722 0 : CALL ERROR("CtqmcoffdiagComplex_measCorrelation : QMC not set ")
2723 : !write(6,*) "not available"
2724 0 : stop
2725 :
2726 : size = op%impurity%particles(op%impurity%activeFlavor)%tail
2727 : beta = op%beta
2728 :
2729 : IF ( size .EQ. 0 ) RETURN
2730 :
2731 : inv_dt = op%inv_dt
2732 :
2733 : DO iCdag = 1, size ! first segments
2734 : tCdag = op%impurity%particles(op%impurity%activeFlavor)%list(iCdag,Cdag_)
2735 : tC = op%impurity%particles(op%impurity%activeFlavor)%list(iCdag,C_ )
2736 : index = INT( ( (tC - tCdag) * inv_dt ) + .5d0 ) + 1
2737 : op%measCorrelation(index,1,iflavor) = op%measCorrelation(index,1,iflavor) + 1.d0
2738 : MODCYCLE(iCdag+1,size,iCdagBeta)
2739 : index = INT( ( ( &
2740 : op%impurity%particles(op%impurity%activeFlavor)%list(iCdagBeta,Cdag_) - tC &
2741 : + AINT(DBLE(iCdag)/DBLE(size))*beta &
2742 : ) * inv_dt ) + .5d0 ) + 1
2743 : IF ( index .LT. 1 .OR. index .GT. op%samples+1 ) THEN
2744 : CALL WARN("CtqmcoffdiagComplex_measCorrelation : bad index line 1095 ")
2745 : ELSE
2746 : op%measCorrelation(index,2,iflavor) = op%measCorrelation(index,2,iflavor) + 1.d0
2747 : END IF
2748 : ! DO iC = 1, size
2749 : ! tC = impurity%particles(impurity%activeFlavor)%list(C_,iC)
2750 : ! time = tC - tCdag
2751 : ! IF ( time .LT. 0.d0 ) time = time + beta
2752 : ! index = INT( ( time * inv_dt ) + .5d0 ) + 1
2753 : ! op%measCorrelation(index,3,iflavor) = op%measCorrelation(index,3,iflavor) + 1.d0
2754 : ! END DO
2755 : DO iC = 1, size! op%Greens(iflavor)%index_old%tail
2756 : !todoba op%measCorrelation(op%Greens(iflavor)%map%listINT(iC+(iCdag-1)*size),3,iflavor) = &
2757 : !todoba op%measCorrelation(op%Greens(iflavor)%map%listINT(iC+(iCdag-1)*size),3,iflavor) + 1.d0
2758 : END DO
2759 : END DO
2760 :
2761 : END SUBROUTINE CtqmcoffdiagComplex_measCorrelation
2762 : !!***
2763 :
2764 : !!****f* ABINIT/m_CtqmcoffdiagComplex/CtqmcoffdiagComplex_measPerturbation
2765 : !! NAME
2766 : !! CtqmcoffdiagComplex_measPerturbation
2767 : !!
2768 : !! FUNCTION
2769 : !! measure perturbation order
2770 : !!
2771 : !! COPYRIGHT
2772 : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
2773 : !! This file is distributed under the terms of the
2774 : !! GNU General Public License, see ~abinit/COPYING
2775 : !! or http://www.gnu.org/copyleft/gpl.txt .
2776 : !!
2777 : !! INPUTS
2778 : !! op=ctqmc
2779 : !! iflavor=the flavor to measure
2780 : !!
2781 : !! OUTPUT
2782 : !!
2783 : !! SIDE EFFECTS
2784 : !!
2785 : !! NOTES
2786 : !!
2787 : !! SOURCE
2788 :
2789 0 : SUBROUTINE CtqmcoffdiagComplex_measPerturbation(op, iflavor)
2790 :
2791 : !Arguments ------------------------------------
2792 : TYPE(CtqmcoffdiagComplex) , INTENT(INOUT) :: op
2793 : !TYPE(ImpurityOperator), INTENT(IN ) :: impurity
2794 : INTEGER , INTENT(IN ) :: iflavor
2795 : !Local variables ------------------------------
2796 : INTEGER :: index
2797 :
2798 0 : IF ( .NOT. op%set ) &
2799 0 : CALL ERROR("CtqmcoffdiagComplex_measiPerturbation : QMC not set ")
2800 :
2801 0 : index = op%impurity%particles(op%impurity%activeFlavor)%tail + 1
2802 0 : IF ( index .LE. op%opt_order ) &
2803 0 : op%measPerturbation(index,iflavor) = op%measPerturbation(index,iflavor) + 1.d0
2804 0 : IF ( index == 1 ) THEN
2805 0 : IF (op%impurity%particles(iflavor)%list(0,C_) < op%impurity%particles(iflavor)%list(0,Cdag_) ) THEN
2806 0 : op%meas_fullemptylines(1,iflavor) = op%meas_fullemptylines(1,iflavor) + 1.d0
2807 : ELSE
2808 0 : op%meas_fullemptylines(2,iflavor) = op%meas_fullemptylines(2,iflavor) + 1.d0
2809 : ENDIF
2810 : ENDIF
2811 :
2812 0 : END SUBROUTINE CtqmcoffdiagComplex_measPerturbation
2813 : !!***
2814 : !#endif
2815 :
2816 : !!****f* ABINIT/m_CtqmcoffdiagComplex/CtqmcoffdiagComplex_getResult
2817 : !! NAME
2818 : !! CtqmcoffdiagComplex_getResult
2819 : !!
2820 : !! FUNCTION
2821 : !! reduce everything to get the result of the simulation
2822 : !!
2823 : !! COPYRIGHT
2824 : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder,F. Gendron)
2825 : !! This file is distributed under the terms of the
2826 : !! GNU General Public License, see ~abinit/COPYING
2827 : !! or http://www.gnu.org/copyleft/gpl.txt .
2828 : !!
2829 : !! INPUTS
2830 : !! op=ctqmc
2831 : !!
2832 : !! OUTPUT
2833 : !!
2834 : !! SIDE EFFECTS
2835 : !!
2836 : !! NOTES
2837 : !!
2838 : !! SOURCE
2839 :
2840 0 : SUBROUTINE CtqmcoffdiagComplex_getResult(op,Iatom,fname)
2841 :
2842 :
2843 : #ifdef HAVE_MPI1
2844 : include 'mpif.h'
2845 : #endif
2846 : !Arguments ------------------------------------
2847 : TYPE(CtqmcoffdiagComplex) , INTENT(INOUT) :: op
2848 : INTEGER, INTENT(IN ) :: Iatom
2849 : character(len=fnlen), INTENT(INOUT) :: fname
2850 : !Local variables ------------------------------
2851 : INTEGER :: iflavor
2852 : ! INTEGER :: iflavor1
2853 : ! INTEGER :: iflavor2
2854 : INTEGER :: flavors
2855 : INTEGER :: itau
2856 : INTEGER :: endDensity
2857 : character(len=2) :: atomnb
2858 : DOUBLE PRECISION :: inv_flavors
2859 : DOUBLE PRECISION :: a
2860 : DOUBLE PRECISION :: b
2861 : DOUBLE PRECISION :: r
2862 0 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: alpha
2863 0 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: beta
2864 0 : COMPLEX(kind=8), ALLOCATABLE, DIMENSION(:) :: measN_1
2865 : DOUBLE PRECISION, DIMENSION(1:2) :: TabX
2866 : DOUBLE PRECISION, DIMENSION(1:2) :: TabY
2867 0 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: freqs
2868 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: counts
2869 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: displs
2870 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: occtot
2871 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: spintot
2872 0 : INTEGER, ALLOCATABLE, DIMENSION(:,:) :: occ
2873 : INTEGER :: sp1,spinmax,spinmin,dspin,nelec,spin
2874 : INTEGER :: spAll
2875 : INTEGER :: last
2876 : INTEGER :: n1
2877 : INTEGER :: n2,n3,quotient,remainder,signe
2878 : INTEGER :: debut
2879 : DOUBLE PRECISION :: signvaluemeassum
2880 : COMPLEX(KIND=8) :: phasevaluemeassum
2881 : ! INTEGER :: fin
2882 : #ifdef HAVE_MPI
2883 : INTEGER :: ierr
2884 : DOUBLE PRECISION, DIMENSION(1) :: arr
2885 : #endif
2886 : INTEGER :: sizeoper,nbprocs,myrank
2887 : DOUBLE PRECISION :: inv_size,sumh,sumtot
2888 0 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: buffer
2889 0 : COMPLEX(KIND=8), ALLOCATABLE, DIMENSION(:,:,:) :: buffer2,buffer2s
2890 0 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: fullempty
2891 0 : TYPE(FFTHyb) :: FFTmrka
2892 :
2893 : #if defined HAVE_MPI && !defined HAVE_MPI2_INPLACE
2894 : DOUBLE PRECISION, ALLOCATABLE , DIMENSION(:) :: buffer1_out,freqs_buf
2895 : DOUBLE PRECISION, ALLOCATABLE , DIMENSION(:,:) :: buffer2_out
2896 : #endif
2897 :
2898 0 : IF ( .NOT. op%done ) &
2899 0 : CALL ERROR("CtqmcoffdiagComplex_getResult : Simulation not run ")
2900 :
2901 0 : flavors = op%flavors
2902 0 : inv_flavors = 1.d0 / DBLE(flavors)
2903 :
2904 :
2905 0 : inv_size = 1.d0 / DBLE(op%size)
2906 0 : sp1 = 0
2907 0 : spAll = 0
2908 :
2909 : !#ifdef CTCtqmcoffdiagComplex_CHECK
2910 0 : IF ( op%opt_check .GT. 0 ) THEN
2911 0 : op%errorImpurity = ImpurityOperator_getError(op%Impurity) * inv_flavors
2912 0 : op%errorBath = BathOperatoroffdiagComplex_getError(op%Bath) * inv_flavors
2913 : END IF
2914 : !#endif
2915 :
2916 0 : MALLOC(alpha,(1,1))
2917 0 : MALLOC(beta,(1,1))
2918 0 : MALLOC(buffer,(1,1))
2919 0 : IF ( op%opt_noise .EQ. 1) THEN
2920 0 : FREEIF(alpha)
2921 0 : MALLOC(alpha,(1:op%samples+1,1:flavors))
2922 0 : FREEIF(beta)
2923 0 : MALLOC(beta,(1:op%samples+1,1:flavors))
2924 : END IF
2925 :
2926 0 : IF ( op%have_MPI .EQV. .TRUE.) THEN
2927 0 : sp1 = 0
2928 0 : spAll = sp1 + flavors + 6
2929 :
2930 : !#ifdef CTCtqmcoffdiagComplex_ANALYSIS
2931 : IF ( op%opt_analysis .EQ. 1 ) &
2932 : spAll = spAll + 3*sp1
2933 0 : IF ( op%opt_order .GT. 0 ) &
2934 0 : spAll = spAll + op%opt_order
2935 0 : IF ( op%opt_noise .EQ. 1 ) &
2936 0 : spAll = spAll + 2*(op%samples + 1)
2937 : !#endif
2938 :
2939 0 : FREEIF(buffer)
2940 0 : MALLOC(buffer,(1:spAll,1:MAX(2,flavors)))
2941 : END IF
2942 :
2943 : ! op%seg_added = op%seg_added * inv_flavors
2944 : ! op%seg_removed = op%seg_removed * inv_flavors
2945 : ! op%seg_sign = op%seg_sign * inv_flavors
2946 : ! op%anti_added = op%anti_added * inv_flavors
2947 : ! op%anti_removed = op%anti_removed * inv_flavors
2948 : ! op%anti_sign = op%anti_sign * inv_flavors
2949 0 : op%stats(:) = op%stats(:) * inv_flavors
2950 :
2951 0 : DO iflavor = 1, flavors
2952 : ! Accumulate last values of N (see also ctqmc_measn)
2953 0 : op%measN(1,iflavor) = op%measN(1,iflavor) + op%measN(3,iflavor)*op%measN(4,iflavor)
2954 0 : op%measN(2,iflavor) = op%measN(2,iflavor) + op%measN(4,iflavor)
2955 : ! Reduction
2956 0 : op%measN(1,iflavor) = op%measN(1,iflavor) / ( op%measN(2,iflavor) * op%beta )
2957 : ! Correction
2958 : !#ifdef CTCtqmcoffdiagComplex_ANALYSIS
2959 0 : IF ( op%opt_order .GT. 0 ) &
2960 : op%measPerturbation(: ,iflavor) = op%measPerturbation(:,iflavor) &
2961 0 : / SUM(op%measPerturbation(:,iflavor))
2962 0 : IF ( op%opt_order .GT. 0 ) &
2963 : op%meas_fullemptylines(: ,iflavor) = op%meas_fullemptylines(:,iflavor) &
2964 0 : / SUM(op%meas_fullemptylines(:,iflavor))
2965 : !write(6,*) "sum fullempty",iflavor,op%meas_fullemptylines(:,iflavor)
2966 :
2967 0 : IF ( op%opt_analysis .EQ. 1 ) THEN
2968 : op%measCorrelation (:,1,iflavor) = op%measCorrelation (:,1,iflavor) &
2969 : / SUM(op%measCorrelation (:,1,iflavor)) &
2970 0 : * op%inv_dt
2971 : op%measCorrelation (:,2,iflavor) = op%measCorrelation (:,2,iflavor) &
2972 : / SUM(op%measCorrelation (:,2,iflavor)) &
2973 0 : * op%inv_dt
2974 : op%measCorrelation (:,3,iflavor) = op%measCorrelation (:,3,iflavor) &
2975 : / SUM(op%measCorrelation (:,3,iflavor)) &
2976 0 : * op%inv_dt
2977 : END IF
2978 : !#endif
2979 0 : IF ( op%opt_noise .EQ. 1 ) THEN
2980 0 : TabX(1) = DBLE(op%modNoise2)
2981 0 : TabX(2) = DBLE(op%modNoise1)
2982 0 : DO itau = 1, op%samples+1
2983 : op%measNoiseG(itau,iflavor,2)%vec = -op%measNoiseG(itau,iflavor,2)%vec*op%inv_dt &
2984 0 : /(op%beta*DBLE(op%modNoise2))
2985 : op%measNoiseG(itau,iflavor,1)%vec = -op%measNoiseG(itau,iflavor,1)%vec*op%inv_dt &
2986 0 : /(op%beta*DBLE(op%modNoise1))
2987 0 : n2 = op%measNoiseG(itau,iflavor,2)%tail
2988 0 : TabY(1) = Stat_deviation(op%measNoiseG(itau,iflavor,2)%vec(1:n2))!*SQRT(n2/(n2-1))
2989 0 : n1 = op%measNoiseG(itau,iflavor,1)%tail
2990 0 : TabY(2) = Stat_deviation(op%measNoiseG(itau,iflavor,1)%vec(1:n1))!*SQRT(n1/(n1-1))
2991 0 : CALL Stat_powerReg(TabX,SQRT(2.d0*LOG(2.d0))*TabY,alpha(itau,iflavor),beta(itau,iflavor),r)
2992 : ! ecart type -> 60%
2993 : ! largeur a mi-hauteur d'une gaussienne -> sqrt(2*ln(2))*sigma
2994 : END DO
2995 : END IF
2996 :
2997 : END DO
2998 : !sui!write(6,*) "getresults"
2999 0 : CALL GreenHyboffdiagComplex_measHybrid(op%Greens, op%Bath%M, op%Impurity%Particles, .TRUE.,op%signvalue,op%phasevalue)
3000 0 : CALL GreenHyboffdiagComplex_getHybrid(op%Greens)
3001 : !write(6,*) "op%measN",op%measN(1,:)
3002 0 : MALLOC(measN_1,(flavors))
3003 0 : do iflavor=1,flavors
3004 0 : measN_1(iflavor)=op%measN(1,iflavor)
3005 : enddo
3006 0 : CALL GreenHyboffdiagComplex_setN(op%Greens, measN_1(:))
3007 0 : FREE(measN_1)
3008 :
3009 : ! todoab case _nd and _d are not completely described.
3010 : FREEIF(buffer2)
3011 : FREEIF(buffer2s)
3012 0 : sizeoper=size(op%Greens%oper,1)
3013 : !write(6,*) "sss",size(op%Greens%oper,1),sizeoper
3014 : !write(6,*) "sss",size(op%Greens%oper,2),flavors
3015 : !write(6,*) "sss",size(op%Greens%oper,3),flavors
3016 0 : MALLOC(buffer2,(1:sizeoper,flavors,flavors))
3017 0 : MALLOC(buffer2s,(1:sizeoper,flavors,flavors))
3018 0 : MALLOC(fullempty,(2,flavors))
3019 : !sui!write(6,*) "greens1"
3020 0 : IF ( op%have_MPI .EQV. .TRUE. ) THEN
3021 : !sui!write(6,*) "greens2"
3022 0 : fullempty=0.d0
3023 0 : buffer2 = op%Greens%oper
3024 : !write(6,*) "buffer2",(op%Greens%oper(1,n1,n1),n1=1,flavors)
3025 0 : buffer2s= 0.d0
3026 : do iflavor=1,flavors
3027 : do itau=1,sizeoper
3028 : !sui!write(6,*) "greens",iflavor,itau,op%Greens%oper(itau,iflavor,iflavor)
3029 : enddo
3030 : enddo
3031 : !write(6,*) "beforempi",op%Greens%oper(1,1,1) ,buffer2(1,1,1)
3032 : #ifdef HAVE_MPI
3033 0 : CALL MPI_COMM_SIZE(op%MY_COMM,nbprocs,ierr)
3034 0 : CALL MPI_COMM_RANK(op%MY_COMM,myrank,ierr)
3035 : #endif
3036 : !write(6,*) "procs",nbprocs,myrank
3037 : END IF
3038 0 : last = sp1
3039 :
3040 0 : op%measDE(:,:) = op%measDE(:,:) * DBLE(op%measurements) /(DBLE(op%sweeps)*op%beta)
3041 :
3042 0 : IF ( op%opt_histo .GT. 0 ) THEN
3043 0 : op%occup_histo_time(:) = op%occup_histo_time(:) / INT(op%sweeps/op%measurements)
3044 0 : op%occupconfig(:) = op%occupconfig(:) / INT(op%sweeps/op%measurements)
3045 0 : op%suscep(:,:) = op%suscep(:,:) / INT(op%sweeps/op%measurements)
3046 0 : op%chi(:,:) = op%chi(:,:) / INT(op%sweeps/op%measurements)
3047 0 : op%chicharge(:,:) = op%chicharge(:,:) / INT(op%sweeps/op%measurements)
3048 0 : op%ntot(:) = op%ntot(:) / INT(op%sweeps/op%measurements)
3049 : END IF
3050 : ! HISTO before MPI_SUM
3051 : ! write(6,*) "=== Histogram of occupations for complete simulation 3 ====",INT(op%sweeps/op%measurements)
3052 : ! sumh=0
3053 : ! do n1=1,op%flavors+1
3054 : ! write(6,'(i4,f10.4)') n1-1, op%occup_histo_time(n1)
3055 : ! sumh=sumh+op%occup_histo_time(n1)
3056 : ! enddo
3057 : ! write(6,*) "=================================",sumh
3058 :
3059 0 : n1 = op%measNoise(1)%tail
3060 0 : n2 = op%measNoise(2)%tail
3061 :
3062 : ! On utilise freqs comme tableau de regroupement
3063 : ! Gather de Noise1
3064 0 : IF ( op%have_MPI .EQV. .TRUE. ) THEN
3065 0 : MALLOC(counts,(1:op%size))
3066 0 : MALLOC(displs,(1:op%size))
3067 : FREEIF(freqs)
3068 0 : MALLOC(freqs,(1:op%size*n1))
3069 0 : freqs = 0.d0
3070 0 : counts(:) = n1
3071 0 : displs(:) = (/ ( iflavor*n1, iflavor=0, op%size-1 ) /)
3072 : #ifdef HAVE_MPI
3073 : #if defined HAVE_MPI2_INPLACE
3074 0 : freqs(n1*op%rank+1:n1*(op%rank+1)) = op%measNoise(1)%vec(1:n1)
3075 : CALL MPI_ALLGATHERV(MPI_IN_PLACE, 0, MPI_DOUBLE_PRECISION, &
3076 : freqs, counts, displs, &
3077 0 : MPI_DOUBLE_PRECISION, op%MY_COMM, ierr)
3078 : #else
3079 : MALLOC(freqs_buf,(n1))
3080 : freqs_buf(1:n2)=op%measNoise(2)%vec(1:n1)
3081 : CALL MPI_ALLGATHERV(freqs_buf, n1, MPI_DOUBLE_PRECISION, &
3082 : freqs, counts, displs, &
3083 : MPI_DOUBLE_PRECISION, op%MY_COMM, ierr)
3084 : FREE(freqs_buf)
3085 : #endif
3086 : #endif
3087 0 : n1 = op%size*n1
3088 0 : CALL Vector_setSize(op%measNoise(1),n1)
3089 0 : op%measNoise(1)%vec(1:n1) = freqs(:)
3090 : ! Gather de Noise2
3091 0 : FREE(freqs)
3092 0 : MALLOC(freqs,(1:op%size*n2))
3093 0 : freqs = 0.d0
3094 0 : counts(:) = n2
3095 0 : displs(:) = (/ ( iflavor*n2, iflavor=0, op%size-1 ) /)
3096 : #ifdef HAVE_MPI
3097 : #if defined HAVE_MPI2_INPLACE
3098 0 : freqs(n2*op%rank+1:n2*(op%rank+1)) = op%measNoise(2)%vec(1:n2)
3099 : CALL MPI_ALLGATHERV(MPI_IN_PLACE, 0, MPI_DOUBLE_PRECISION, &
3100 : freqs, counts, displs, &
3101 0 : MPI_DOUBLE_PRECISION, op%MY_COMM, ierr)
3102 : #else
3103 : MALLOC(freqs_buf,(n2))
3104 : freqs_buf(1:n2)=op%measNoise(2)%vec(1:n2)
3105 : CALL MPI_ALLGATHERV(freqs_buf, n2, MPI_DOUBLE_PRECISION, &
3106 : freqs, counts, displs, &
3107 : MPI_DOUBLE_PRECISION, op%MY_COMM, ierr)
3108 : FREE(freqs_buf)
3109 : #endif
3110 : #endif
3111 0 : n2 = op%size*n2
3112 0 : CALL Vector_setSize(op%measNoise(2),n2)
3113 0 : op%measNoise(2)%vec(1:n2) = freqs(:)
3114 0 : FREE(counts)
3115 0 : FREE(displs)
3116 0 : FREE(freqs)
3117 : END IF
3118 : !n1 = op%measNoise(1)%tail
3119 : !n2 = op%measNoise(2)%tail
3120 :
3121 : ! Transformation des paquets pour que ca fit a CTQMC_SLICE(1|2)
3122 0 : IF ( n1 .GT. CTQMC_SLICE1 ) THEN
3123 0 : itau = n1/CTQMC_SLICE1
3124 0 : MALLOC(freqs,(1:n1/itau))
3125 0 : DO debut=1, n1/itau
3126 0 : freqs(debut)=SUM(op%measNoise(1)%vec((debut-1)*itau+1:itau*debut))
3127 : END DO
3128 0 : freqs(:) = freqs(:)/DBLE(itau)
3129 0 : op%modNoise1 = op%modNoise1*itau
3130 0 : n1 = n1/itau
3131 0 : CALL Vector_setSize(op%measNoise(1),n1)
3132 0 : op%measNoise(1)%vec(1:n1) = freqs(:)
3133 0 : FREE(freqs)
3134 : END IF
3135 0 : IF ( n2 .GT. CTQMC_SLICE1*CTQMC_SLICE2 ) THEN
3136 0 : itau = n2/(CTQMC_SLICE1*CTQMC_SLICE2)
3137 0 : MALLOC(freqs,(1:n2/itau))
3138 0 : DO debut=1, n2/itau
3139 0 : freqs(debut)=SUM(op%measNoise(2)%vec((debut-1)*itau+1:itau*debut))
3140 : END DO
3141 0 : freqs(:) = freqs(:)/DBLE(itau)
3142 0 : op%modNoise2 = op%modNoise2*itau
3143 0 : n2 = n2/itau
3144 0 : CALL Vector_setSize(op%measNoise(2),n2)
3145 0 : op%measNoise(2)%vec(1:n2) = freqs(:)
3146 0 : FREE(freqs)
3147 : END IF
3148 : ! On peut s'amuser avec nos valeur d'energies
3149 : !MALLOC(TabX,(1:20))
3150 : !MALLOC(TabY,(1:20))
3151 :
3152 0 : TabX(1) = DBLE(op%modNoise2)
3153 0 : TabX(2) = DBLE(op%modNoise1)
3154 :
3155 : ! Il faut calculer pour chaque modulo 10 ecarts type sur les donnes acquises
3156 0 : op%measNoise(1)%vec(1:n1) = op%measNoise(1)%vec(1:n1)/(op%beta*DBLE(op%modNoise1))*DBLE(op%measurements)
3157 0 : op%measNoise(2)%vec(1:n2) = op%measNoise(2)%vec(1:n2)/(op%beta*DBLE(op%modNoise2))*DBLE(op%measurements)
3158 : ! CALL Vector_print(op%measNoise(1),op%rank+70)
3159 : ! CALL Vector_print(op%measNoise(2),op%rank+50)
3160 : ! DO iflavor=1,10
3161 : ! debut = (iflavor-1)*n2/10+1
3162 : ! fin = iflavor*n2/10
3163 : ! TabY(iflavor) = Stat_deviation(op%measNoise(2)%vec(debut:fin))
3164 : ! debut = (iflavor-1)*n1/10+1
3165 : ! fin = iflavor*n1/10
3166 : ! TabY(10+iflavor) = Stat_deviation(op%measNoise(1)%vec(debut:fin))
3167 : ! END DO
3168 : !! TabY(1:n) = (op%measNoise(2)%vec(1:n) &
3169 : !! )
3170 : !! !/(op%beta*DBLE(op%modNoise2))*DBLE(op%measurements) &
3171 : !! !- op%measDE(1,1))
3172 : !! TabY(op%measNoise(2)%tail+1:n+op%measNoise(2)%tail) = (op%measNoise(1)%vec(1:n) &
3173 : !! )
3174 : !! ! /(op%beta*DBLE(op%modNoise1))*DBLE(op%measurements) &
3175 : !! ! - op%measDE(1,1))
3176 : ! IF ( op%rank .EQ. 0 ) THEN
3177 : ! DO iflavor=1,20
3178 : ! write(45,*) TabX(iflavor), TabY(iflavor)
3179 : ! END DO
3180 : ! END IF
3181 : !
3182 :
3183 :
3184 0 : TabY(1) = Stat_deviation(op%measNoise(2)%vec(1:n2))!*SQRT(n2/(n2-1))
3185 : !! write(op%rank+10,*) TabX(2)
3186 : !! write(op%rank+40,*) TabX(1)
3187 : !! CALL Vector_print(op%measNoise(1),op%rank+10)
3188 : !! CALL Vector_print(op%measNoise(2),op%rank+40)
3189 : !! CLOSE(op%rank+10)
3190 : !! CLOSE(op%rank+40)
3191 0 : TabY(2) = Stat_deviation(op%measNoise(1)%vec(1:n1))!*SQRT(n1/(n1-1))
3192 : !! ! Ecart carre moyen ~ ecart type mais non biaise. Serait moins precis. Aucun
3193 : ! impact sur la pente, juste sur l'ordonnee a l'origine.
3194 :
3195 0 : CALL Stat_powerReg(TabX,SQRT(2.d0*LOG(2.d0))*TabY,a,b,r)
3196 : ! FREE(TabX)
3197 : ! FREE(TabY)
3198 : ! ecart type -> 60%
3199 : ! largeur a mi-hauteur d'une gaussienne -> sqrt(2*ln(2))*sigma
3200 :
3201 : !op%measDE(1,1) = SUM(op%measNoise(1)%vec(1:op%measNoise(1)%tail))/(DBLE(op%measNoise(1)%tail*op%modNoise1)*op%beta)
3202 : !op%measDE(2:flavors,1:flavors) = op%measDE(2:flavors,1:flavors) /(DBLE(op%sweeps)*op%beta)
3203 0 : CALL ImpurityOperator_getErrorOverlap(op%Impurity,op%measDE)
3204 : ! Add the difference between true calculation and quick calculation of the
3205 : ! last sweep overlap to measDE(2,2)
3206 : !op%measDE = op%measDE * DBLE(op%measurements)
3207 0 : IF ( op%have_MPI .EQV. .TRUE. ) THEN
3208 0 : IF ( op%opt_analysis .EQ. 1 ) THEN
3209 0 : buffer(last+1:last+sp1,:) = op%measCorrelation(:,1,:)
3210 0 : last = last + sp1
3211 0 : buffer(last+1:last+sp1,:) = op%measCorrelation(:,2,:)
3212 0 : last = last + sp1
3213 0 : buffer(last+1:last+sp1,:) = op%measCorrelation(:,3,:)
3214 : last = last + sp1
3215 : END IF
3216 0 : IF ( op%opt_order .GT. 0 ) THEN
3217 0 : buffer(last+1:last+op%opt_order, :) = op%measPerturbation(:,:)
3218 : last = last + op%opt_order
3219 : END IF
3220 0 : IF ( op%opt_noise .EQ. 1 ) THEN
3221 0 : buffer(last+1:last+op%samples+1,:) = alpha(:,:)
3222 0 : last = last + op%samples + 1
3223 0 : buffer(last+1:last+op%samples+1,:) = beta(:,:)
3224 0 : last = last + op%samples + 1
3225 : END IF
3226 : ! op%measDE(2,2) = a*EXP(b*LOG(DBLE(op%sweeps*op%size)))
3227 0 : buffer(spall-(flavors+5):spAll-6,:) = op%measDE(:,:)
3228 : ! buffer(spAll ,1) = op%seg_added
3229 : ! buffer(spAll-1,1) = op%seg_removed
3230 : ! buffer(spAll-2,1) = op%seg_sign
3231 : ! buffer(spAll ,2) = op%anti_added
3232 : ! buffer(spAll-1,2) = op%anti_removed
3233 : ! buffer(spAll-2,2) = op%anti_sign
3234 0 : buffer(spAll ,1) = op%stats(1)
3235 0 : buffer(spAll-1,1) = op%stats(2)
3236 0 : buffer(spAll-2,1) = op%stats(3)
3237 0 : buffer(spAll ,2) = op%stats(4)
3238 0 : buffer(spAll-1,2) = op%stats(5)
3239 0 : buffer(spAll-2,2) = op%stats(6)
3240 0 : buffer(spAll-3,1) = op%swap
3241 0 : buffer(spAll-3,2) = DBLE(op%modGlobalMove(2))
3242 0 : buffer(spAll-4,1) = a
3243 0 : buffer(spAll-4,2) = b
3244 : !#ifdef CTCtqmcoffdiagComplex_CHECK
3245 0 : buffer(spAll-5,1) = op%errorImpurity
3246 0 : buffer(spAll-5,2) = op%errorBath
3247 0 : signvaluemeassum = 0
3248 0 : phasevaluemeassum = 0
3249 : !#endif
3250 :
3251 : #ifdef HAVE_MPI
3252 : #if defined HAVE_MPI2_INPLACE
3253 : CALL MPI_ALLREDUCE(MPI_IN_PLACE, buffer, spAll*flavors, &
3254 0 : MPI_DOUBLE_PRECISION, MPI_SUM, op%MY_COMM, ierr)
3255 : #else
3256 : MALLOC(buffer2_out,(spAll,flavors))
3257 : CALL MPI_ALLREDUCE(buffer, buffer2_out, spAll*flavors, &
3258 : MPI_DOUBLE_PRECISION, MPI_SUM, op%MY_COMM, ierr)
3259 : buffer(1:spAll,1:flavors)=buffer2_out(1:spAll,1:flavors)
3260 : FREE(buffer2_out)
3261 : #endif
3262 : CALL MPI_ALLREDUCE(buffer2, buffer2s, sizeoper*flavors*flavors, &
3263 0 : MPI_DOUBLE_COMPLEX, MPI_SUM, op%MY_COMM, ierr)
3264 0 : CALL MPI_ALLREDUCE([op%runtime], arr, 1, MPI_DOUBLE_PRECISION, MPI_MAX, op%MY_COMM, ierr)
3265 0 : op%runtime=arr(1)
3266 0 : CALL MPI_ALLREDUCE([op%Greens%signvaluemeas], arr, 1, MPI_DOUBLE_PRECISION, MPI_SUM, op%MY_COMM, ierr)
3267 0 : signvaluemeassum=arr(1)
3268 0 : CALL MPI_ALLREDUCE([op%Greens%phasevaluemeas], arr , 1, MPI_DOUBLE_COMPLEX, MPI_SUM, op%MY_COMM, ierr)
3269 0 : phasevaluemeassum=arr(1)
3270 : #if defined HAVE_MPI2_INPLACE
3271 0 : IF ( op%opt_histo .GT. 0 ) THEN
3272 : CALL MPI_ALLREDUCE(MPI_IN_PLACE,op%occup_histo_time, flavors+1, MPI_DOUBLE_PRECISION, MPI_SUM, &
3273 0 : op%MY_COMM, ierr)
3274 : CALL MPI_ALLREDUCE(MPI_IN_PLACE, op%occupconfig, 2**flavors, MPI_DOUBLE_PRECISION, MPI_SUM, &
3275 0 : op%MY_COMM, ierr)
3276 : CALL MPI_ALLREDUCE(MPI_IN_PLACE, op%suscep, 3*op%samples, MPI_DOUBLE_PRECISION, MPI_SUM, &
3277 0 : op%MY_COMM, ierr)
3278 : CALL MPI_ALLREDUCE(MPI_IN_PLACE, op%chi, 3*op%samples, MPI_DOUBLE_PRECISION, MPI_SUM, &
3279 0 : op%MY_COMM, ierr)
3280 : CALL MPI_ALLREDUCE(MPI_IN_PLACE, op%chicharge, 3*op%samples, MPI_DOUBLE_PRECISION, MPI_SUM, &
3281 0 : op%MY_COMM, ierr)
3282 : CALL MPI_ALLREDUCE(MPI_IN_PLACE, op%ntot, 3, MPI_DOUBLE_PRECISION, MPI_SUM, &
3283 0 : op%MY_COMM, ierr)
3284 : END IF
3285 : CALL MPI_ALLREDUCE(MPI_IN_PLACE, sumh, 1, MPI_DOUBLE_PRECISION, MPI_SUM, &
3286 0 : op%MY_COMM, ierr)
3287 : #else
3288 : MALLOC(buffer1_out,(flavors+1))
3289 : CALL MPI_ALLREDUCE(op%occup_histo_time, buffer1_out, flavors+1, MPI_DOUBLE_PRECISION, MPI_SUM, &
3290 : op%MY_COMM, ierr)
3291 : op%occup_histo_time(1:flavors+1) =buffer1_out(1:flavors+1)
3292 : FREE(buffer1_out)
3293 : #endif
3294 0 : CALL MPI_ALLREDUCE([sumh], arr, 1, MPI_DOUBLE_PRECISION, MPI_SUM, op%MY_COMM, ierr)
3295 0 : sumh=arr(1)
3296 0 : IF ( op%opt_order .GT. 0 ) THEN
3297 : CALL MPI_ALLREDUCE(op%meas_fullemptylines, fullempty, 2*flavors, MPI_DOUBLE_PRECISION, MPI_SUM, &
3298 0 : op%MY_COMM, ierr)
3299 : ENDIF
3300 : #endif
3301 :
3302 :
3303 0 : buffer = buffer * inv_size
3304 0 : op%measDE(:,:) = buffer(spall-(flavors+5):spAll-6,:)
3305 : ! op%seg_added = buffer(spAll ,1)
3306 : ! op%seg_removed = buffer(spAll-1,1)
3307 : ! op%seg_sign = buffer(spAll-2,1)
3308 : ! op%anti_added = buffer(spAll ,2)
3309 : ! op%anti_removed = buffer(spAll-1,2)
3310 : ! op%anti_sign = buffer(spAll-2,2)
3311 0 : op%stats(1) = buffer(spAll ,1)
3312 0 : op%stats(2) = buffer(spAll-1,1)
3313 0 : op%stats(3) = buffer(spAll-2,1)
3314 0 : op%stats(4) = buffer(spAll ,2)
3315 0 : op%stats(5) = buffer(spAll-1,2)
3316 0 : op%stats(6) = buffer(spAll-2,2)
3317 0 : op%swap = buffer(spAll-3,1)
3318 0 : op%modGlobalMove(2) = NINT(buffer(spAll-3,2))
3319 0 : a = buffer(spAll-4,1)
3320 0 : b = buffer(spAll-4,2)
3321 : !!#ifdef CTCtqmcoffdiagComplex_CHECK
3322 0 : op%errorImpurity= buffer(spAll-5,1)
3323 0 : op%errorBath = buffer(spAll-5,2)
3324 : !#endif
3325 :
3326 : ! DO iflavor = 1, flavors
3327 : ! op%Greens(iflavor)%oper = buffer(1:sp1 , iflavor)
3328 : ! END DO
3329 0 : op%Greens%oper = buffer2s/float(nbprocs)
3330 : ! write(6,*) "buffer2s",(op%Greens%oper(1,n1,n1),n1=1,flavors)
3331 0 : op%Greens%signvaluemeas = signvaluemeassum/float(nbprocs)
3332 0 : op%Greens%phasevaluemeas = phasevaluemeassum/float(nbprocs)
3333 : !sui!write(6,*) "nbprocs",nbprocs,op%Greens%signvaluemeas
3334 0 : op%Greens%oper = op%Greens%oper / op%Greens%phasevaluemeas
3335 : ! write(6,*) "buffer3s",(op%Greens%oper(1,n1,n1),n1=1,flavors)
3336 0 : IF ( op%opt_order .GT. 0 ) THEN
3337 0 : op%meas_fullemptylines= fullempty/float(nbprocs)
3338 : ENDIF
3339 0 : do iflavor=1,flavors
3340 : do itau=1,sizeoper
3341 : !sui!write(6,*) "greens_av",iflavor,itau,op%Greens%oper(itau,iflavor,iflavor)
3342 : enddo
3343 : enddo
3344 : !write(6,*) "aftermpi",op%Greens%oper(1,1,1) ,buffer2s(1,1,1)
3345 0 : last = sp1
3346 : IF ( op%opt_analysis .EQ. 1 ) THEN
3347 : op%measCorrelation(:,1,:) = buffer(last+1:last+sp1,:)
3348 : last = last + sp1
3349 : op%measCorrelation(:,2,:) = buffer(last+1:last+sp1,:)
3350 : last = last + sp1
3351 : op%measCorrelation(:,3,:) = buffer(last+1:last+sp1,:)
3352 : last = last + sp1
3353 : END IF
3354 0 : IF ( op%opt_order .GT. 0 ) THEN
3355 0 : op%measPerturbation(:,:) = buffer(last+1:last+op%opt_order, :)
3356 : last = last + op%opt_order
3357 : END IF
3358 0 : IF ( op%opt_noise .EQ. 1 ) THEN
3359 0 : alpha(:,:) = buffer(last+1:last+op%samples+1,:)
3360 0 : last = last + op%samples + 1
3361 0 : beta(:,:) = buffer(last+1:last+op%samples+1,:)
3362 0 : last = last + op%samples + 1
3363 : END IF
3364 : END IF
3365 0 : DO iflavor = 1, flavors
3366 : ! complete DE matrix
3367 0 : op%measDE(iflavor, iflavor+1:flavors) = op%measDE(iflavor+1:flavors,iflavor)
3368 : END DO
3369 0 : FREE(buffer)
3370 0 : FREE(buffer2)
3371 0 : FREE(buffer2s)
3372 0 : FREE(fullempty)
3373 :
3374 0 : IF ( op%opt_spectra .GE. 1 ) THEN
3375 : !PROBLEM eos_gnu_13.2_mpich . %endDensity was introduced throughout
3376 0 : endDensity = SIZE(op%density,2)
3377 : ! endDensity=op%endDensity
3378 : !ENDPROBLEM
3379 :
3380 0 : IF ( op%density(1,endDensity) .EQ. -1.d0 ) &
3381 0 : endDensity = endDensity - 1
3382 0 : CALL FFTHyb_init(FFTmrka,endDensity,DBLE(op%thermalization)/DBLE(op%measurements*op%opt_spectra))
3383 : ! Not very Beauty
3384 0 : MALLOC(freqs,(1:FFTmrka%size/2))
3385 0 : DO iflavor = 1, flavors
3386 : ! mean value is removed to supress the continue composent
3387 0 : CALL FFTHyb_setData(FFTmrka,op%density(iflavor,1:endDensity)/op%beta+real(op%Greens%oper(op%samples+1,iflavor,iflavor)))
3388 0 : CALL FFTHyb_run(FFTmrka,1)
3389 0 : CALL FFTHyb_getData(FFTmrka,endDensity,op%density(iflavor,:),freqs)
3390 : END DO
3391 0 : op%density(flavors+1,:) = -1.d0
3392 0 : op%density(flavors+1,1:FFTmrka%size/2) = freqs
3393 0 : CALL FFTHyb_destroy(FFTmrka)
3394 0 : FREE(freqs)
3395 : END IF
3396 :
3397 0 : op%a_Noise = a
3398 0 : op%b_Noise = b
3399 0 : IF ( op%opt_noise .EQ. 1 ) THEN
3400 0 : op%abNoiseG(1,:,:) = alpha
3401 0 : op%abNoiseG(2,:,:) = beta
3402 : END IF
3403 0 : FREE(alpha)
3404 0 : FREE(beta)
3405 0 : IF ( op%opt_histo .GT. 0 ) THEN
3406 0 : write(op%ostream,*) "=== Histogram of occupations for complete simulation ===="
3407 : ! write(6,*) "sumh over procs", sumh
3408 0 : sumh=0
3409 0 : do n1=1,op%flavors+1
3410 0 : write(op%ostream,'(i4,f10.4)') n1-1, op%occup_histo_time(n1)/float(nbprocs)
3411 0 : sumh=sumh+op%occup_histo_time(n1)/float(nbprocs)
3412 : enddo
3413 0 : write(op%ostream,'(a,f10.4)') " all" , sumh
3414 0 : write(op%ostream,*) "================================="
3415 :
3416 0 : MALLOC(occ,(2**op%flavors,1:flavors))
3417 : #ifdef FC_LLVM
3418 : MALLOC(occtot,(2**op%flavors) )
3419 : MALLOC(spintot,(2**op%flavors) )
3420 : #else
3421 0 : MALLOC(occtot,(2**op%flavors))
3422 0 : MALLOC(spintot,(2**op%flavors))
3423 : #endif
3424 0 : do n1=1,2**op%flavors
3425 : ! Compute occupations of individual Orbitals
3426 0 : n3=n1-1
3427 0 : occtot(n1)=0
3428 0 : spintot(n1)=0
3429 0 : signe=1
3430 0 : do n2=1,op%flavors
3431 0 : remainder=modulo(n3,2)
3432 0 : quotient=(n3-remainder)/2
3433 0 : occ(n1,n2)=remainder
3434 0 : n3=quotient
3435 0 : occtot(n1)=occtot(n1)+occ(n1,n2)
3436 0 : if(n2>=7) signe =-1
3437 : !if(n2>=7) signe =0
3438 0 : spintot(n1)=spintot(n1)+occ(n1,n2)*signe
3439 : enddo
3440 0 : op%occupconfig(n1)=op%occupconfig(n1)/float(nbprocs)
3441 : enddo
3442 :
3443 0 : write(op%ostream,*) "=== Histogram of occupations of configurations for complete simulation ===="
3444 0 : sumh=0
3445 0 : if(op%flavors==14) then
3446 0 : do n1=1,2**op%flavors
3447 0 : write(op%ostream,'(i4,14i2,f20.2)') n1, (occ(n1,n2),n2=1,op%flavors),op%occupconfig(n1)
3448 0 : sumh=sumh+op%occupconfig(n1)
3449 : enddo
3450 0 : else if(op%flavors==10) then
3451 0 : do n1=1,2**op%flavors
3452 0 : write(op%ostream,'(i4,10i2,f20.2)') n1, (occ(n1,n2),n2=1,op%flavors),op%occupconfig(n1)
3453 0 : sumh=sumh+op%occupconfig(n1)
3454 : enddo
3455 : end if
3456 0 : write(op%ostream,'(a,f10.4)') " all" , sumh
3457 :
3458 0 : sumtot=0
3459 0 : do nelec=0,10
3460 0 : spinmin=modulo(nelec,2)
3461 : if(nelec<=5) spinmax=nelec
3462 : if(nelec>=6) spinmax=10-nelec
3463 0 : spinmax=6
3464 0 : dspin=1
3465 0 : write(op%ostream,*) "=== Histogram of occupations of configurations for total number of electrons",nelec
3466 0 : do spin=spinmin,spinmax,dspin
3467 0 : sumh=0
3468 0 : do n1=1,2**op%flavors
3469 0 : if(occtot(n1)==nelec.and.abs(spintot(n1))==spin) then
3470 0 : sumh=sumh+op%occupconfig(n1)
3471 0 : if(op%flavors==14) then
3472 0 : write(op%ostream,'(i8,14i2,a,i2,i3,f10.4)') n1,(occ(n1,n2),n2=1,op%flavors)," ",occtot(n1),spintot(n1),&
3473 0 : & op%occupconfig(n1)
3474 0 : else if(op%flavors==10) then
3475 0 : write(op%ostream,'(i8,10i2,a,i2,i3,f10.4)') n1,(occ(n1,n2),n2=1,op%flavors)," ",occtot(n1),spintot(n1),&
3476 0 : op%occupconfig(n1)
3477 : end if
3478 : endif
3479 : enddo
3480 0 : write(op%ostream,'(a,i4,a,i4,a,f10.4)') " === Sum of weights for",nelec," electrons and spin",spin," is ",sumh
3481 0 : sumtot=sumtot+sumh
3482 : enddo
3483 : enddo
3484 0 : write(op%ostream,'(a,f10.4)') "Full sum is",sumtot
3485 :
3486 0 : FREE(occ)
3487 0 : FREE(occtot)
3488 0 : FREE(spintot)
3489 : !END IF
3490 :
3491 : !==============================
3492 : ! Write Susceptibilities
3493 : !==============================
3494 0 : if(Iatom .lt. 10) then
3495 0 : write(atomnb, '("0",i1)') Iatom
3496 : else
3497 0 : write(atomnb, '(i2)') Iatom
3498 : end if
3499 : ! Local Magnetic Susceptibility
3500 0 : if(op%opt_histo .gt. 1) then
3501 : ! Scalar
3502 0 : if(op%nspinor .eq. 1) then
3503 0 : open (unit=735,file=trim(fname)//'_LocalSpinSuscept_atom_'//atomnb//'.dat',status='unknown',form='formatted')
3504 0 : write(735,*) '#Tau Total t2g eg'
3505 0 : do n1=1,op%samples
3506 0 : op%suscep(:,n1)=op%suscep(:,n1)/float(nbprocs)/float(op%samples)
3507 0 : write(735,'(1x,f14.8,2x,f12.8,2x,f12.8,2x,f12.8)') (n1-1)*op%beta/op%samples,(op%suscep(n2,n1),n2=1,3)
3508 : enddo
3509 : !add tau=beta
3510 0 : write(735,'(1x,f14.8,2x,f12.8,2x,f12.8,2x,f12.8)') (op%samples)*op%beta/op%samples,(op%suscep(n2,1),n2=1,3)
3511 :
3512 : else
3513 : ! SOC
3514 0 : open (unit=735,file=trim(fname)//'_LocalMagnSuscept_atom_'//atomnb//'.dat',status='unknown',form='formatted')
3515 0 : write(735,*) '#Tau Total Orbital Spin'
3516 0 : do n1=1,op%samples
3517 0 : op%chi(:,n1) = op%chi(:,n1)/float(nbprocs)/float(op%samples)
3518 0 : write(735,'(1x,f14.8,2x,f12.8,2x,f12.8,2x,f12.8)') (n1-1)*op%beta/op%samples,(op%chi(n2,n1),n2=1,3)
3519 : end do
3520 : !add tau=beta
3521 0 : write(735,'(1x,f14.8,2x,f12.8,2x,f12.8,2x,f12.8)') (op%samples)*op%beta/op%samples,(op%chi(n2,1),n2=1,3)
3522 : endif
3523 0 : close(unit=735)
3524 : endif
3525 :
3526 : ! Local Charge Susceptiblity
3527 0 : if(op%opt_histo .gt. 2) then
3528 0 : op%ntot(:) = op%ntot(:)/float(nbprocs)/float(op%samples)
3529 0 : open (unit=735,file=trim(fname)//'_LocalChargeSuscept_atom_'//atomnb//'.dat',status='unknown',form='formatted')
3530 0 : write(735,*) '#Tau Total <ntot>'
3531 0 : do n1=1,op%samples
3532 0 : op%chicharge(1,n1)=(op%chicharge(1,n1)/float(nbprocs)/float(op%samples))-(op%ntot(1)*op%ntot(1))
3533 : !op%chicharge(2,n1)=(op%chicharge(2,n1)/float(nbprocs)/float(op%samples))-(op%ntot(2)*op%ntot(2))
3534 : !op%chicharge(3,n1)=(op%chicharge(3,n1)/float(nbprocs)/float(op%samples))-(op%ntot(3)*op%ntot(3))
3535 0 : write(735,'(1x,f14.8,2x,f12.8,2x,f12.8,2x,f12.8)') (n1-1)*op%beta/op%samples,(op%chicharge(1,n1)),op%ntot(1)
3536 : enddo
3537 : !add tau=beta
3538 0 : write(735,'(1x,f14.8,2x,f12.8,2x,f12.8,2x,f12.8)') (op%samples)*op%beta/op%samples,(op%chicharge(1,1)),op%ntot(1)
3539 0 : close(unit=735)
3540 : endif
3541 :
3542 : END IF
3543 :
3544 :
3545 0 : END SUBROUTINE CtqmcoffdiagComplex_getResult
3546 : !!***
3547 :
3548 : !!****f* ABINIT/m_CtqmcoffdiagComplex/CtqmcoffdiagComplex_symmetrizeGreen
3549 : !! NAME
3550 : !! CtqmcoffdiagComplex_symmetrizeGreen
3551 : !!
3552 : !! FUNCTION
3553 : !! optionnaly symmetrize the green functions
3554 : !!
3555 : !! COPYRIGHT
3556 : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
3557 : !! This file is distributed under the terms of the
3558 : !! GNU General Public License, see ~abinit/COPYING
3559 : !! or http://www.gnu.org/copyleft/gpl.txt .
3560 : !!
3561 : !! INPUTS
3562 : !! op=ctqmc
3563 : !! syms=weight factors
3564 : !!
3565 : !! OUTPUT
3566 : !! argout(sizeout)=description
3567 : !!
3568 : !! SIDE EFFECTS
3569 : !!
3570 : !! NOTES
3571 : !!
3572 : !! SOURCE
3573 :
3574 0 : SUBROUTINE CtqmcoffdiagComplex_symmetrizeGreen(op, syms)
3575 :
3576 : !Arguments ------------------------------------
3577 : TYPE(CtqmcoffdiagComplex) , INTENT(INOUT) :: op
3578 : DOUBLE PRECISION, DIMENSION(:,:), INTENT(IN ) :: syms
3579 : !Local variables ------------------------------
3580 : !INTEGER :: iflavor1
3581 : !INTEGER :: iflavor2
3582 : !INTEGER :: flavors
3583 : !DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: green_tmp
3584 : !DOUBLE PRECISION, ALLOCATABLE, DIMENSION(: ) :: n_tmp
3585 :
3586 : ABI_UNUSED((/syms(1,1), op%swap/))
3587 :
3588 : ! flavors = op%flavors
3589 : ! IF ( SIZE(syms,1) .NE. flavors .OR. SIZE(syms,2) .NE. flavors ) THEN
3590 : ! CALL WARNALL("CtqmcoffdiagComplex_symmetrizeGreen : wrong opt_sym -> not symmetrizing")
3591 : ! RETURN
3592 : ! END IF
3593 : !
3594 : ! MALLOC(green_tmp,(1:op%samples+1,flavors))
3595 : ! green_tmp(:,:) = 0.d0
3596 : ! MALLOC(n_tmp,(1:flavors))
3597 : ! n_tmp(:) = 0.d0
3598 : ! DO iflavor1=1, flavors
3599 : ! DO iflavor2=1,flavors
3600 : ! green_tmp(:,iflavor1) = green_tmp(:,iflavor1) &
3601 : ! + syms(iflavor2,iflavor1) * op%Greens(iflavor2)%oper(:)
3602 : ! n_tmp(iflavor1) = n_tmp(iflavor1) &
3603 : ! + syms(iflavor2,iflavor1) * op%measN(1,iflavor2)
3604 : ! END DO
3605 : ! END DO
3606 : ! DO iflavor1=1, flavors
3607 : ! op%Greens(iflavor1)%oper(:) = green_tmp(:,iflavor1)
3608 : ! op%measN(1,iflavor1) = n_tmp(iflavor1)
3609 : ! END DO
3610 : ! FREE(green_tmp)
3611 : ! FREE(n_tmp)
3612 0 : END SUBROUTINE CtqmcoffdiagComplex_symmetrizeGreen
3613 : !!***
3614 :
3615 : !!****f* ABINIT/m_CtqmcoffdiagComplex/CtqmcoffdiagComplex_getGreen
3616 : !! NAME
3617 : !! CtqmcoffdiagComplex_getGreen
3618 : !!
3619 : !! FUNCTION
3620 : !! Get the full green functions in time and/or frequency
3621 : !!
3622 : !! COPYRIGHT
3623 : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
3624 : !! This file is distributed under the terms of the
3625 : !! GNU General Public License, see ~abinit/COPYING
3626 : !! or http://www.gnu.org/copyleft/gpl.txt .
3627 : !!
3628 : !! INPUTS
3629 : !! op=ctqmc
3630 : !!
3631 : !! OUTPUT
3632 : !! Gtau=green function in time
3633 : !! Gw=green function in frequency
3634 : !!
3635 : !! SIDE EFFECTS
3636 : !!
3637 : !! NOTES
3638 : !!
3639 : !! SOURCE
3640 :
3641 0 : SUBROUTINE CtqmcoffdiagComplex_getGreen(op, Gtau, Gw)
3642 :
3643 : !Arguments ------------------------------------
3644 : USE m_GreenHyboffdiagComplex
3645 : TYPE(CtqmcoffdiagComplex) , INTENT(INOUT) :: op
3646 : COMPLEX(KIND=8), DIMENSION(:,:,:), OPTIONAL, INTENT(INOUT) :: Gtau
3647 : COMPLEX(KIND=8), DIMENSION(:,:,:), OPTIONAL, INTENT(INOUT) :: Gw
3648 : !Local variables ------------------------------
3649 : !INTEGER :: itime
3650 : INTEGER :: iflavor1
3651 : INTEGER :: iflavor1b !,iflavor,iflavorbis
3652 : INTEGER :: iflavor2
3653 : INTEGER :: iflavor3
3654 : INTEGER :: flavors,tail
3655 : INTEGER :: ifreq,itime
3656 : COMPLEX(KIND=8) :: u1
3657 : COMPLEX(KIND=8) :: u2
3658 : COMPLEX(KIND=8) :: u3
3659 : COMPLEX(KIND=8) :: Un
3660 : COMPLEX(KIND=8) :: UUnn
3661 : DOUBLE PRECISION :: iw !omega,
3662 : CHARACTER(LEN=4) :: cflavors
3663 : CHARACTER(LEN=50) :: string
3664 0 : TYPE(GreenHyboffdiagComplex) :: F_tmp
3665 :
3666 0 : flavors = op%flavors
3667 0 : DO iflavor1 = 1, flavors
3668 0 : u1 = cmplx(0.d0,0.d0,kind=8)
3669 0 : u2 = cmplx(0.d0,0.d0,kind=8)
3670 0 : u3 = cmplx(0.d0,0.d0,kind=8)
3671 0 : DO iflavor2 = 1, flavors
3672 0 : IF ( iflavor2 .EQ. iflavor1 ) CYCLE
3673 0 : Un = op%Impurity%mat_U(iflavor2,iflavor1) * op%measN(1,iflavor2)
3674 : ! Un = op%Impurity%mat_U(iflavor2,iflavor1) * (op%Greens%oper(1,iflavor2,iflavor2) + 1.d0)
3675 : !write(6,*) "forsetmoments",iflavor1,iflavor2,(op%Greens%oper(1,iflavor2,iflavor2) + 1.d0), Un
3676 0 : u1 = u1 + Un
3677 0 : u2 = u2 + Un*op%Impurity%mat_U(iflavor2,iflavor1)
3678 0 : DO iflavor3 = 1, flavors
3679 0 : IF ( iflavor3 .EQ. iflavor2 .OR. iflavor3 .EQ. iflavor1 ) CYCLE
3680 : UUnn = (op%Impurity%mat_U(iflavor2,iflavor1)*op%Impurity%mat_U(iflavor3,iflavor1)) * &
3681 0 : & op%measDE(iflavor2,iflavor3)
3682 0 : u2 = u2 + UUnn
3683 : END DO
3684 : END DO
3685 : ! write(6,*) "u1,u2",u1,u2
3686 :
3687 0 : DO iflavor1b = 1, flavors
3688 0 : u3 =-(op%Impurity%mat_U(iflavor1,iflavor1b))*op%Greens%oper(1,iflavor1,iflavor1b)
3689 : ! u3=U_{1,1b}*G_{1,1b}
3690 0 : CALL GreenHyboffdiagComplex_setMoments(op%Greens,iflavor1,iflavor1b,u1,u2,u3)
3691 : END DO ! iflavor1b
3692 :
3693 : END DO ! iflavor1
3694 :
3695 0 : IF ( PRESENT( Gtau ) ) THEN
3696 0 : DO iflavor1 = 1, flavors
3697 0 : DO iflavor2 = 1, flavors
3698 0 : Gtau(1:op%samples,iflavor1,iflavor2) = op%Greens%oper(1:op%samples,iflavor1,iflavor2)
3699 : END DO
3700 : END DO ! iflavor1
3701 : END IF
3702 : ! !--------- Write Occupation matrix before Gtau
3703 : ! write(ostream,'(17x,a)') "Occupation matrix"
3704 : ! write(ostream,'(17x,30i10)') (iflavorbis,iflavorbis=1,op%flavors)
3705 : ! write(ostream,'(18x,30a)') ("----------",iflavorbis=1,op%flavors)
3706 : ! do iflavor=1, op%flavors
3707 : ! write(ostream,'(7x,i10,a,30f10.4)') iflavor,"|",(-op%Greens%oper(op%samples,iflavor,iflavorbis),iflavorbis=1,op%flavors)
3708 : ! enddo
3709 : ! write(ostream,'(18x,30a)') ("----------",iflavorbis=1,op%flavors),ch10
3710 : ! !------------------------------------------------------------------------------------------
3711 : ! !--------- Write Occupation matrix Gtau
3712 : ! write(ostream,'(17x,a)') "Occupation matrix"
3713 : ! write(ostream,'(17x,30i10)') (iflavorbis,iflavorbis=1,op%flavors)
3714 : ! write(ostream,'(18x,30a)') ("----------",iflavorbis=1,op%flavors)
3715 : ! do iflavor=1, op%flavors
3716 : ! write(ostream,'(7x,i10,a,30f10.4)') iflavor,"|",(Gtau(op%samples,iflavor,iflavorbis),iflavorbis=1,op%flavors)
3717 : ! enddo
3718 : ! write(ostream,'(18x,30a)') ("----------",iflavorbis=1,op%flavors),ch10
3719 : ! !------------------------------------------------------------------------------------------
3720 :
3721 : !================================================
3722 : if(3==4) then
3723 : !================================================
3724 : DO iflavor1 = 1, flavors
3725 : DO iflavor1b = 1, flavors
3726 : !call nfourier3(op%Greens%oper(1:op%samples,iflavor1,iflavor1b),Gw(1:op%samples,iflavor1,iflavor1b),iflavor1==iflavor1b,op%Greens%samples,op%Greens%samples-1,op%Greens%beta,1.d0,op%Greens%Mk(iflavor1,iflavor1b,1),op%Greens%Mk(iflavor1,iflavor1b,2),op%Greens%Mk(iflavor1,iflavor1b,3))
3727 : END DO
3728 : END DO ! iflavor1
3729 : ! ============== write Gomega_nd.dat
3730 : if(op%rank==0) then
3731 : OPEN(UNIT=44, FILE="Gomega_nd_nfourier2.dat")
3732 : WRITE(cflavors,'(I4)') 2*(flavors*flavors+1)
3733 : string = '(1x,'//TRIM(ADJUSTL(cflavors))//'E15.5)'
3734 : !write(6,*) " op%Greens%Wmax", op%Greens%Wmax
3735 : do iflavor1=1, flavors
3736 : do iflavor1b=1, flavors
3737 : write(44,*) "#op%Greens%Mk(iflavor1,iflavor2,1",op%Greens%Mk(iflavor1,iflavor1b,:)
3738 : DO ifreq = 1, op%samples
3739 : ! !write(6,string) (DBLE(ifreq)*2-1)*3.1415/op%Greens%beta, &
3740 : ! (/ ((real(Gw(ifreq,iflavor1,iflavor1b)),imag(Gw(ifreq,iflavor1,iflavor1b)), iflavor1=1, flavors),iflavor1b=1,flavors) /)
3741 : ! WRITE(44,string) (DBLE(ifreq)*2.d0-1.d0)*3.1415926/op%Greens%beta, &
3742 : iw=aimag(Gw(ifreq,op%flavors,op%flavors+1))
3743 : WRITE(44,string) aimag(Gw(ifreq,op%flavors,op%flavors+1)),&
3744 : real(Gw(ifreq,iflavor1,iflavor1b)),aimag(Gw(ifreq,iflavor1,iflavor1b)),&
3745 : ( -op%Greens%Mk(iflavor1,iflavor1b,2) )/(iw*iw) , (op%Greens%Mk(iflavor1,iflavor1b,1))/iw!-op%Greens%Mk(iflavor1,iflavor1b,3)/(iw*iw))/iw
3746 : ! WRITE(102,*) aimag(Gw(ifreq,op%flavors,op%flavors+1)), (op%Greens%Mk(iflavor1,iflavor1b,1))/iw,op%Greens%Mk(iflavor1,iflavor1b,1),iw
3747 : END DO
3748 : WRITE(44,*)
3749 : END DO
3750 : END DO
3751 : close(44)
3752 : endif
3753 : !================================================
3754 : endif
3755 : !================================================
3756 : !!write(6,*) "present gw", present(gw)
3757 0 : IF ( PRESENT( Gw ) ) THEN
3758 : !!write(6,*) "size gw",SIZE(Gw,DIM=2) ,flavors+1
3759 0 : IF ( SIZE(Gw,DIM=3) .EQ. flavors+1 ) THEN
3760 : ! CALL GreenHyboffdiag_forFourier(op%Greens, Gomega=Gw, omega=Gw(:,op%flavors,op%flavors+1))
3761 0 : CALL GreenHyboffdiagComplex_forFourierComplex(op%Greens, Gomega=Gw, omega=Gw(:,op%flavors,op%flavors+1))
3762 : !write(6,*) "1"
3763 : !IF ( op%rank .EQ. 0 ) write(20,*) Gw(:,iflavor1)
3764 0 : ELSE IF ( SIZE(Gw,DIM=3) .EQ. flavors ) THEN
3765 0 : CALL GreenHyboffdiagComplex_forFourierComplex(op%Greens,Gomega=Gw)
3766 : !write(6,*) "2"
3767 : ELSE
3768 0 : CALL WARNALL("CtqmcoffdiagComplex_getGreen : Gw is not valid ")
3769 0 : CALL GreenHyboffdiagComplex_forFourierComplex(op%Greens,Wmax=op%Wmax)
3770 : !write(6,*) "3"
3771 : END IF
3772 : ELSE
3773 0 : CALL GreenHyboffdiagComplex_forFourierComplex(op%Greens,Wmax=op%Wmax)
3774 : END IF
3775 : ! ============== write Gomega_nd.dat
3776 : !================================================
3777 : ! if(3==4) then
3778 : !================================================
3779 : if(op%rank==0.and.3==4) then
3780 : OPEN(UNIT=44, FILE="Gomega_nd.dat")
3781 : WRITE(cflavors,'(I4)') 2*(flavors*flavors+1)
3782 : string = '(1x,'//TRIM(ADJUSTL(cflavors))//'E15.5)'
3783 : !write(6,*) " op%Greens%Wmax", op%Greens%Wmax
3784 : do iflavor1=1, flavors
3785 : do iflavor1b=1, flavors
3786 : write(44,*) "#op%Greens%Mk(iflavor1,iflavor2,1",op%Greens%Mk(iflavor1,iflavor1b,:)
3787 : DO ifreq = 1, SIZE(Gw,1)
3788 : ! !write(6,string) (DBLE(ifreq)*2-1)*3.1415/op%Greens%beta, &
3789 : ! (/ ((real(Gw(ifreq,iflavor1,iflavor1b)),imag(Gw(ifreq,iflavor1,iflavor1b)), iflavor1=1, flavors),iflavor1b=1,flavors) /)
3790 : ! WRITE(44,string) (DBLE(ifreq)*2.d0-1.d0)*3.1415926/op%Greens%beta, &
3791 : iw=aimag(Gw(ifreq,op%flavors,op%flavors+1))
3792 : WRITE(44,string) aimag(Gw(ifreq,op%flavors,op%flavors+1)),&
3793 : real(Gw(ifreq,iflavor1,iflavor1b)),aimag(Gw(ifreq,iflavor1,iflavor1b)),&
3794 : ( -op%Greens%Mk(iflavor1,iflavor1b,2) )/(iw*iw) , &
3795 : & (op%Greens%Mk(iflavor1,iflavor1b,1)-op%Greens%Mk(iflavor1,iflavor1b,3)/(iw*iw))/iw
3796 : END DO
3797 : WRITE(44,*)
3798 : END DO
3799 : END DO
3800 : endif
3801 : !================================================
3802 : ! endif
3803 : !================================================
3804 :
3805 :
3806 : ! ==============================
3807 : ! --- Initialize F_tmp
3808 : !write(6,*) "10"
3809 :
3810 0 : IF ( op%have_MPI .EQV. .TRUE. ) THEN
3811 : !CALL GreenHyboffdiag_init(F_tmp,op%samples,op%beta,op%flavors,MY_COMM=op%MY_COMM)
3812 0 : CALL GreenHyboffdiagComplex_init(F_tmp,op%samples,op%beta,flavors)
3813 : !write(6,*) "10a"
3814 : ELSE
3815 0 : CALL GreenHyboffdiagComplex_init(F_tmp,op%samples,op%beta,flavors)
3816 : !write(6,*) "10b"
3817 : END IF
3818 :
3819 : !write(6,*) "11"
3820 : ! CALL GreenHyboffdiag_setOperW(F_tmp,Gw)
3821 :
3822 0 : tail = op%samples
3823 0 : F_tmp%Wmax=op%samples ! backFourier only works for linear freq: calculation of A and etc..
3824 0 : MALLOC(F_tmp%oper_w,(1:tail,op%flavors,op%flavors))
3825 0 : F_tmp%oper_w(1:tail,1:F_tmp%nflavors,1:F_tmp%nflavors) = Gw(1:tail,1:F_tmp%nflavors,1:F_tmp%nflavors)
3826 : !write(6,*) "example",F_tmp%oper_w(1,1,1)
3827 : !write(6,*) "example",Gw(1,1,1)
3828 0 : F_tmp%setW = .TRUE.
3829 : !write(6,*) size(F_tmp%oper_w,1)
3830 : !write(6,*) size(F_tmp%oper_w,2)
3831 : !write(6,*) size(F_tmp%oper_w,3)
3832 : !write(6,*) size(Gw,1)
3833 : !write(6,*) size(Gw,2)
3834 : !write(6,*) size(Gw,3)
3835 :
3836 : !write(6,*) "eee", (2.d0*DBLE(ifreq)-1.d0) * 3.1415/op%beta,real(F_tmp%oper_w(1,1,1)),imag(F_tmp%oper_w(1,1,1))
3837 : !================================================
3838 : if(3==4) then
3839 : !================================================
3840 : OPEN(UNIT=3337, FILE="Gomega_nd2.dat")
3841 : do iflavor1=1, flavors
3842 : do iflavor1b=1, flavors
3843 : do ifreq=1, tail
3844 : ! write(3337,*) (2.d0*DBLE(ifreq)-1.d0) * 3.1415/op%beta,real(F_tmp%oper_w(ifreq,iflavor1,iflavor1b)),&
3845 : ! & imag(F_tmp%oper_w(ifreq,iflavor1,iflavor1b))
3846 : write(3337,*) aimag(Gw(ifreq,op%flavors,op%flavors+1)), real(F_tmp%oper_w(ifreq,iflavor1,iflavor1b)),&
3847 : & aimag(F_tmp%oper_w(ifreq,iflavor1,iflavor1b))
3848 :
3849 : ! omega=(2.d0*DBLE(ifreq)-1.d0) * 3.1415/op%beta
3850 : ! F_tmp%oper_w(ifreq,iflavor1,iflavor1b)=0.1**2/Gw(ifreq,op%flavors,op%flavors+1)
3851 : enddo
3852 : write(3337,*)
3853 : enddo
3854 : enddo
3855 : close(3337)
3856 : !================================================
3857 : endif
3858 : !================================================
3859 :
3860 : !write(6,*) "12",F_tmp%Wmax
3861 :
3862 : ! CALL GreenHyboffdiag_backFourier(F_tmp,func="green")
3863 :
3864 : !write(6,*) "13"
3865 :
3866 : !================================================
3867 : if(3==4) then
3868 : !================================================
3869 : OPEN(UNIT=48, FILE="Gtau_nd_2.dat")
3870 : ! --- Print full non diagonal Gtau in Gtau_nd.dat
3871 : WRITE(cflavors,'(I4)') flavors*flavors+1
3872 : string = '(1x,'//TRIM(ADJUSTL(cflavors))//'ES22.14)'
3873 : DO itime = 1, op%samples+1
3874 : WRITE(48,string) DBLE(itime-1)*op%beta/DBLE(op%samples), &
3875 : & ((F_tmp%oper(itime,iflavor1,iflavor1b), iflavor1=1, flavors),iflavor1b=1,flavors)
3876 : END DO
3877 : !================================================
3878 : endif
3879 : !================================================
3880 :
3881 0 : CALL GreenHyboffdiagComplex_destroy(F_tmp)
3882 :
3883 : !FREE(F_tmp%oper_w)
3884 : ! ==============================
3885 0 : END SUBROUTINE CtqmcoffdiagComplex_getGreen
3886 : !!***
3887 :
3888 : !!****f* ABINIT/m_CtqmcoffdiagComplex/CtqmcoffdiagComplex_getD
3889 : !! NAME
3890 : !! CtqmcoffdiagComplex_getD
3891 : !!
3892 : !! FUNCTION
3893 : !! get double occupation
3894 : !!
3895 : !! COPYRIGHT
3896 : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
3897 : !! This file is distributed under the terms of the
3898 : !! GNU General Public License, see ~abinit/COPYING
3899 : !! or http://www.gnu.org/copyleft/gpl.txt .
3900 : !!
3901 : !! INPUTS
3902 : !! op=ctqmc
3903 : !!
3904 : !! OUTPUT
3905 : !! D=full double occupation
3906 : !!
3907 : !! SIDE EFFECTS
3908 : !!
3909 : !! NOTES
3910 : !!
3911 : !! SOURCE
3912 :
3913 0 : SUBROUTINE CtqmcoffdiagComplex_getD(op, D)
3914 :
3915 : !Arguments ------------------------------------
3916 : TYPE(CtqmcoffdiagComplex) , INTENT(IN ) :: op
3917 : DOUBLE PRECISION, INTENT(OUT) :: D
3918 : !Local variables ------------------------------
3919 : INTEGER :: iflavor1
3920 : INTEGER :: iflavor2
3921 :
3922 0 : D = 0.d0
3923 :
3924 0 : DO iflavor1 = 1, op%flavors
3925 0 : DO iflavor2 = iflavor1+1, op%flavors
3926 0 : D = D + op%measDE(iflavor2,iflavor1)
3927 : END DO
3928 : END DO
3929 : !IF ( op%rank .EQ. 0 ) THEN
3930 : ! DO iflavor1 = 1, op%flavors
3931 : ! DO iflavor2 = iflavor1+1, op%flavors
3932 : ! write(4533,*) op%measDE(iflavor2,iflavor1)k
3933 : ! write(4534,*) op%Impurity%mat_U(iflavor2,iflavor1)k
3934 : ! END DO
3935 : ! END DO
3936 :
3937 : !ENDIF
3938 :
3939 0 : END SUBROUTINE CtqmcoffdiagComplex_getD
3940 : !!***
3941 :
3942 : !!****f* ABINIT/m_CtqmcoffdiagComplex/CtqmcoffdiagComplex_getE
3943 : !! NAME
3944 : !! CtqmcoffdiagComplex_getE
3945 : !!
3946 : !! FUNCTION
3947 : !! get interaction energy and noise on it
3948 : !!
3949 : !! COPYRIGHT
3950 : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
3951 : !! This file is distributed under the terms of the
3952 : !! GNU General Public License, see ~abinit/COPYING
3953 : !! or http://www.gnu.org/copyleft/gpl.txt .
3954 : !!
3955 : !! INPUTS
3956 : !! op=ctqmc
3957 : !!
3958 : !! OUTPUT
3959 : !! E=interaction energy
3960 : !! noise=noise on this value
3961 : !!
3962 : !! SIDE EFFECTS
3963 : !!
3964 : !! NOTES
3965 : !!
3966 : !! SOURCE
3967 :
3968 0 : SUBROUTINE CtqmcoffdiagComplex_getE(op,E,noise)
3969 :
3970 : !Arguments ------------------------------------
3971 : TYPE(CtqmcoffdiagComplex) , INTENT(IN ) :: op
3972 : DOUBLE PRECISION, INTENT(OUT) :: E
3973 : DOUBLE PRECISION, INTENT(OUT) :: Noise
3974 :
3975 0 : E = op%measDE(1,1)
3976 0 : Noise = op%a_Noise*(DBLE(op%sweeps)*DBLE(op%size))**op%b_Noise
3977 0 : END SUBROUTINE CtqmcoffdiagComplex_getE
3978 : !!***
3979 :
3980 : !!****f* ABINIT/m_CtqmcoffdiagComplex/CtqmcoffdiagComplex_printAll
3981 : !! NAME
3982 : !! CtqmcoffdiagComplex_printAll
3983 : !!
3984 : !! FUNCTION
3985 : !! print different functions computed during the simulation
3986 : !!
3987 : !! COPYRIGHT
3988 : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
3989 : !! This file is distributed under the terms of the
3990 : !! GNU General Public License, see ~abinit/COPYING
3991 : !! or http://www.gnu.org/copyleft/gpl.txt .
3992 : !!
3993 : !! INPUTS
3994 : !! op=ctqmc
3995 : !!
3996 : !! OUTPUT
3997 : !!
3998 : !! SIDE EFFECTS
3999 : !!
4000 : !! NOTES
4001 : !!
4002 : !! SOURCE
4003 :
4004 0 : SUBROUTINE CtqmcoffdiagComplex_printAll(op)
4005 :
4006 : !Arguments ------------------------------------
4007 : TYPE(CtqmcoffdiagComplex), INTENT(INOUT) :: op
4008 :
4009 0 : IF ( .NOT. op%done ) &
4010 0 : CALL WARNALL("CtqmcoffdiagComplex_printAll : Simulation not run ")
4011 :
4012 : !sui!write(6,*) "op%stats",op%stats
4013 0 : CALL CtqmcoffdiagComplex_printQMC(op)
4014 :
4015 0 : CALL CtqmcoffdiagComplex_printGreen(op)
4016 :
4017 0 : CALL CtqmcoffdiagComplex_printD(op)
4018 :
4019 : ! CALL CtqmcoffdiagComplex_printE(op)
4020 :
4021 : !#ifdef CTCtqmcoffdiagComplex_ANALYSIS
4022 : ! CALL CtqmcoffdiagComplex_printPerturbation(op)
4023 :
4024 0 : CALL CtqmcoffdiagComplex_printCorrelation(op)
4025 : !#endif
4026 :
4027 0 : CALL CtqmcoffdiagComplex_printSpectra(op)
4028 :
4029 0 : END SUBROUTINE CtqmcoffdiagComplex_printAll
4030 : !!***
4031 :
4032 : !!****f* ABINIT/m_CtqmcoffdiagComplex/CtqmcoffdiagComplex_printQMC
4033 : !! NAME
4034 : !! CtqmcoffdiagComplex_printQMC
4035 : !!
4036 : !! FUNCTION
4037 : !! print ctqmc statistics
4038 : !!
4039 : !! COPYRIGHT
4040 : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
4041 : !! This file is distributed under the terms of the
4042 : !! GNU General Public License, see ~abinit/COPYING
4043 : !! or http://www.gnu.org/copyleft/gpl.txt .
4044 : !!
4045 : !! INPUTS
4046 : !! op=ctqmc
4047 : !!
4048 : !! OUTPUT
4049 : !!
4050 : !! SIDE EFFECTS
4051 : !!
4052 : !! NOTES
4053 : !!
4054 : !! SOURCE
4055 :
4056 0 : SUBROUTINE CtqmcoffdiagComplex_printQMC(op)
4057 :
4058 : !Arguments ------------------------------------
4059 : TYPE(CtqmcoffdiagComplex), INTENT(INOUT) :: op
4060 : !Local variables ------------------------------
4061 : INTEGER :: ostream
4062 : INTEGER :: iflavor,iflavorbis,iorder
4063 : DOUBLE PRECISION :: sweeps
4064 : DOUBLE PRECISION :: invSweeps
4065 : CHARACTER(LEN=2) :: a
4066 : CHARACTER(LEN=15) :: string
4067 :
4068 : !IF ( op%rank .NE. 0) RETURN
4069 0 : IF ( op%rank .NE. MOD(op%size,op%size)) RETURN
4070 :
4071 0 : ostream = op%ostream
4072 0 : sweeps = DBLE(op%sweeps)
4073 0 : invSweeps = 1.d0/sweeps
4074 :
4075 0 : WRITE(ostream,'(a,F13.0,A11,F10.2,A12,I5,A5)') "The QMC Solver did ",sweeps*DBLE(op%size), " sweeps in ", op%runTime, &
4076 0 : " seconds on ", op%size, " CPUs"
4077 0 : WRITE(ostream,'(A28,F6.2)') "Segments added [%] : ", op%stats(4)*invSweeps*100.d0
4078 0 : WRITE(ostream,'(A28,F6.2)') "Segments removed [%] : ", op%stats(5)*invSweeps*100.d0
4079 0 : WRITE(ostream,'(A28,F6.2)') "Segments <0 sign [%] : ", op%stats(6)*invSweeps*100.d0
4080 : !WRITE(ostream,'(A28,F12.2)') "Number of meas [%] : ", op%stats(6)
4081 0 : WRITE(ostream,'(A28,F6.2)') "Anti-segments added [%] : ", op%stats(1)*invSweeps*100.d0
4082 0 : WRITE(ostream,'(A28,F6.2)') "Anti-segments removed [%] : ", op%stats(2)*invSweeps*100.d0
4083 0 : WRITE(ostream,'(A28,F6.2)') "Anti-segments <0 sign [%] : ", op%stats(3)*invSweeps*100.d0
4084 : !WRITE(ostream,'(A28,F12.2)') "Sum of sign [%] : ", op%stats(3)
4085 0 : WRITE(ostream,'(A28,F13.2)') "Signe value : ", op%Greens%signvaluemeas
4086 0 : WRITE(ostream,'(A28,F13.2,F13.2)') "Phase value : ", op%Greens%phasevaluemeas
4087 0 : IF ( op%modGlobalMove(1) .LT. op%sweeps + 1 ) THEN
4088 0 : WRITE(ostream,'(A28,F6.2)') "Global Move [%] : ", op%swap *invSweeps*100.d0*op%modGlobalMove(1)
4089 0 : WRITE(ostream,'(A28,F6.2)') "Global Move Reduced [%] : ", op%swap / DBLE(op%modGlobalMove(2))*100.d0
4090 : END IF
4091 : !#ifdef CTCtqmcoffdiagComplex_CHECK
4092 0 : IF ( op%opt_check .EQ. 1 .OR. op%opt_check .EQ. 3 ) &
4093 0 : WRITE(ostream,'(A28,E22.14)') "Impurity test [%] : ", op%errorImpurity*100.d0
4094 0 : IF ( op%opt_check .GE. 2 ) &
4095 0 : WRITE(ostream,'(A28,E22.14)') "Bath test [%] : ", op%errorBath *100.d0
4096 : !#endif
4097 0 : WRITE(ostream,'(A28,ES22.14,A5,ES21.14)') "<Epot> [U] : ", op%measDE(1,1), " +/- ",&
4098 : !#ifdef HAVE_MPI
4099 0 : op%a_Noise*(sweeps*DBLE(op%size))**op%b_Noise
4100 : !#else
4101 : ! op%a_Noise*(sweeps)**op%b_Noise
4102 : !#endif
4103 : !--------- Write double occupation between all pairs of orbitals --------------------------
4104 0 : write(ostream,'(17x,a)') "Double occupation between pairs of orbitals"
4105 0 : write(ostream,'(17x,30i10)') (iflavorbis,iflavorbis=1,op%flavors)
4106 0 : write(ostream,'(18x,30a)') ("----------",iflavorbis=1,op%flavors)
4107 0 : do iflavor=1, op%flavors
4108 0 : write(ostream,'(7x,i10,a,30f10.4)') iflavor,"|",(op%measDE(iflavor,iflavorbis),iflavorbis=1,op%flavors)
4109 : enddo
4110 0 : write(ostream,'(18x,30a)') ("----------",iflavorbis=1,op%flavors),ch10
4111 : !------------------------------------------------------------------------------------------
4112 :
4113 : !--------- Write number of segments for each orbitals
4114 : ! write(ostream,'(a)') "Number of segments for each orbitals"
4115 : ! write(ostream,'(17x,30i10)') (iflavorbis,iflavorbis=1,op%flavors)
4116 : ! write(ostream,'(17x,30a)') ("----------",iflavorbis=1,op%flavors)
4117 : ! do iflavor=1, op%flavors
4118 : ! write(ostream,'(i17,a,30f10.4)') iflavor,"|",(op%Impurity%particles(IT)%tail)
4119 : ! enddo
4120 : ! write(ostream,'(17x,30a)') ("----------",iflavorbis=1,op%flavors)
4121 : !------------------------------------------------------------------------------------------
4122 : !--------- Write G(L)
4123 0 : write(ostream,'(17x,a)') "G(tau = L)"
4124 0 : write(ostream,'(17x,30i10)') (iflavorbis,iflavorbis=1,op%flavors)
4125 0 : write(ostream,'(18x,30a)') ("----------",iflavorbis=1,op%flavors)
4126 0 : do iflavor=1, op%flavors
4127 0 : write(ostream,'(7x,i10,a,30f10.4)') iflavor,"|",(real(op%Greens%oper(op%samples,iflavor,iflavorbis)),iflavorbis=1,op%flavors)
4128 : enddo
4129 0 : write(ostream,'(18x,30a)') ("----------",iflavorbis=1,op%flavors),ch10
4130 : !------------------------------------------------------------------------------------------
4131 : !--------- Write G(1)
4132 0 : write(ostream,'(17x,a)') "G(tau = 1)"
4133 0 : write(ostream,'(17x,30i10)') (iflavorbis,iflavorbis=1,op%flavors)
4134 0 : write(ostream,'(18x,30a)') ("----------",iflavorbis=1,op%flavors)
4135 0 : do iflavor=1, op%flavors
4136 0 : write(ostream,'(7x,i10,a,30f10.4)') iflavor,"|",(real(op%Greens%oper(1,iflavor,iflavorbis)),iflavorbis=1,op%flavors)
4137 : enddo
4138 0 : write(ostream,'(18x,30a)') ("----------",iflavorbis=1,op%flavors),ch10
4139 : !------------------------------------------------------------------------------------------
4140 :
4141 0 : WRITE(ostream,'(A28,F8.4,A3,F7.4)') "Noise [U] : ", op%a_Noise, " x^", op%b_Noise
4142 0 : WRITE(ostream,'(A28,E10.2)') "Niquist puls. [/beta] : ", ACOS(-1.d0)*op%inv_dt
4143 0 : WRITE(ostream,'(A28,E22.14)') "Max Acc. Epot Error [U] : ", op%measDE(2,2)/(op%beta*op%modNoise1*2.d0)*sweeps
4144 :
4145 : !WRITE(ostream,'(A28,F7.4,A3,F7.4,A4,E20.14)') "Noise [G(tau)] : ", op%a_Noise(2), "x^", op%b_Noise(2), " -> ", &
4146 : !op%a_Noise(2)*(sweeps*DBLE(op%size))**op%b_Noise(2)
4147 : !----- PERTURBATION ORDER------------------------------------------------------------------
4148 0 : IF ( op%opt_order .GT. 0 ) THEN
4149 0 : write(ostream,*)
4150 0 : WRITE(a,'(I2)') op%flavors
4151 0 : string = '(A28,'//TRIM(ADJUSTL(a))//'(1x,I3))'
4152 0 : WRITE(ostream,string) "Perturbation orders : ",(/ (MAXLOC(op%measPerturbation(:, iflavor))-1, iflavor=1, op%flavors) /)
4153 0 : write(ostream,'(17x,a)') "order of Perturbation for flavors"
4154 0 : write(ostream,'(17x,30i10)') (iflavorbis,iflavorbis=1,op%flavors)
4155 0 : write(ostream,'(18x,30a)') ("----------",iflavorbis=1,op%flavors)
4156 0 : write(ostream,'(12x,a,30i10)') " max ",(/ (MAXLOC(op%measPerturbation(:, iflavor))-1, iflavor=1, op%flavors) /)
4157 0 : write(ostream,'(18x,30a)') ("----------",iflavorbis=1,op%flavors)
4158 0 : do iorder=0, op%opt_order-1
4159 0 : write(ostream,'(7x,i10,a,30f10.4)') iorder,"|",(op%measPerturbation(iorder+1,iflavor),iflavor=1,op%flavors)
4160 : enddo
4161 : END IF
4162 : !------------------------------------------------------------------------------------------
4163 : !----- PERTURBATION ORDER------------------------------------------------------------------
4164 0 : IF ( op%opt_order .GT. 0 ) THEN
4165 0 : write(ostream,*)
4166 0 : write(ostream,'(17x,a)') "Proportion of full and empty orbital for order 0"
4167 0 : write(ostream,'(17x,30i10)') (iflavorbis,iflavorbis=1,op%flavors)
4168 0 : write(ostream,'(18x,30a)') ("----------",iflavorbis=1,op%flavors)
4169 0 : write(ostream,'(2x,a,30f10.4)') " full orbital |",(op%meas_fullemptylines(1,iflavor),iflavor=1,op%flavors)
4170 0 : write(ostream,'(2x,a,30f10.4)') " empty orbital |",(op%meas_fullemptylines(2,iflavor),iflavor=1,op%flavors)
4171 : END IF
4172 : !------------------------------------------------------------------------------------------
4173 : !CALL FLUSH(op%ostream)
4174 : IF ( ABS(((op%stats(4) *invSweeps*100.d0) / (op%stats(5) *invSweeps*100.d0) - 1.d0)) .GE. 0.02d0 &
4175 0 : .OR. ABS(((op%stats(1)*invSweeps*100.d0) / (op%stats(2)*invSweeps*100.d0) - 1.d0)) .GE. 0.02d0 ) &
4176 : THEN
4177 0 : CALL WARNALL("CtqmcoffdiagComplex_printQMC : bad statistic according to moves. Increase sweeps")
4178 : END IF
4179 0 : IF ( ABS(op%b_Noise+0.5)/0.5d0 .GE. 0.05d0 ) &
4180 0 : CALL WARNALL("CtqmcoffdiagComplex_printQMC : bad statistic according to Noise. Increase sweeps")
4181 : ! IF ( ISNAN(op%a_Noise) .OR. ISNAN(op%a_Noise) ) &
4182 : ! CALL WARNALL("CtqmcoffdiagComplex_printQMC : NaN appeared. Increase sweeps ")
4183 :
4184 :
4185 : END SUBROUTINE CtqmcoffdiagComplex_printQMC
4186 : !!***
4187 :
4188 : !!****f* ABINIT/m_CtqmcoffdiagComplex/CtqmcoffdiagComplex_printGreen
4189 : !! NAME
4190 : !! CtqmcoffdiagComplex_printGreen
4191 : !!
4192 : !! FUNCTION
4193 : !! print green functions
4194 : !!
4195 : !! COPYRIGHT
4196 : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
4197 : !! This file is distributed under the terms of the
4198 : !! GNU General Public License, see ~abinit/COPYING
4199 : !! or http://www.gnu.org/copyleft/gpl.txt .
4200 : !!
4201 : !! INPUTS
4202 : !! op=ctqmc
4203 : !! oFileIn=file stream
4204 : !!
4205 : !! OUTPUT
4206 : !!
4207 : !! SIDE EFFECTS
4208 : !!
4209 : !! NOTES
4210 : !!
4211 : !! SOURCE
4212 :
4213 0 : SUBROUTINE CtqmcoffdiagComplex_printGreen(op, oFileIn)
4214 :
4215 : !Arguments ------------------------------------
4216 : use m_io_tools, only : flush_unit
4217 : TYPE(CtqmcoffdiagComplex) , INTENT(IN) :: op
4218 : INTEGER , OPTIONAL, INTENT(IN) :: oFileIn
4219 : !Local variables ------------------------------
4220 : INTEGER :: oFile
4221 : INTEGER :: itime
4222 : INTEGER :: sp1
4223 : INTEGER :: iflavor,iflavorb
4224 : INTEGER :: flavors !, iflavor2 !,iflavor1,
4225 : CHARACTER(LEN=4) :: cflavors,hcflavors
4226 : CHARACTER(LEN=50) :: string,hstring
4227 : DOUBLE PRECISION :: dt
4228 : DOUBLE PRECISION :: sweeps
4229 :
4230 : !IF ( op%rank .NE. MOD(1,op%size)) RETURN
4231 0 : IF ( op%rank .NE. MOD(op%size+1,op%size)) RETURN
4232 :
4233 0 : oFile = 40
4234 0 : IF ( PRESENT(oFileIn) ) THEN
4235 0 : oFile = oFileIn
4236 : ELSE
4237 0 : OPEN(UNIT=oFile, FILE="Gtau.dat")
4238 : END IF
4239 0 : OPEN(UNIT=43, FILE="Gtau_nd.dat")
4240 0 : rewind(43)
4241 0 : OPEN(UNIT=143, FILE="Gtau_nd_Re_ctqmc.dat")
4242 0 : rewind(143)
4243 0 : OPEN(UNIT=243, FILE="Gtau_nd_Im_ctqmc.dat")
4244 0 : rewind(243)
4245 0 : sp1 = op%samples
4246 0 : dt = op%beta / DBLE(sp1)
4247 0 : sp1 = sp1 + 1
4248 0 : flavors = op%flavors
4249 0 : sweeps = DBLE(op%sweeps)*DBLE(op%size)
4250 :
4251 0 : IF ( op%opt_noise .EQ. 1) THEN
4252 0 : WRITE(cflavors,'(I4)') (2*flavors+1)*2
4253 0 : string = '(1x,'//TRIM(ADJUSTL(cflavors))//'ES22.14)'
4254 0 : DO itime = 1, sp1
4255 0 : WRITE(oFile,string) DBLE(itime-1)*dt, &
4256 0 : (/ (op%Greens%oper(itime,iflavor,iflavor), iflavor=1, flavors) /), &
4257 0 : (/ (op%abNoiseG(1,itime,iflavor)*(sweeps)**op%abNoiseG(2,itime,iflavor), iflavor=1, flavors) /)
4258 : END DO
4259 : ELSE
4260 0 : WRITE(cflavors,'(I4)') (flavors+1)*2
4261 0 : string = '(1x,'//TRIM(ADJUSTL(cflavors))//'ES22.14)'
4262 0 : DO itime = 1, sp1
4263 : ! WRITE(45,string) DBLE(itime-1)*dt, &
4264 : ! (/ (op%Greens%oper(itime,iflavor,iflavor), iflavor=1, flavors) /)
4265 0 : WRITE(oFile,string) DBLE(itime-1)*dt, &
4266 0 : (/ (op%Greens%oper(itime,iflavor,iflavor), iflavor=1, flavors) /)
4267 : ! WRITE(46,*) DBLE(itime-1)*dt, &
4268 : ! & (/ ((op%Greens%oper(itime,iflavor,iflavorb), iflavor=1, flavors),iflavorb=1,flavors) /)
4269 : END DO
4270 : ! DO itime = 1, sp1
4271 : ! WRITE(47,*) DBLE(itime-1)*dt, &
4272 : ! & (/ ((op%Greens%oper(itime,iflavor,iflavorb), iflavor=1, flavors),iflavorb=1,flavors) /)
4273 : ! END DO
4274 : ! --- Print full non diagonal Gtau in Gtau_nd.dat
4275 0 : WRITE(cflavors,'(I4)') (flavors*flavors+1)
4276 0 : WRITE(hcflavors,'(I4)') (flavors*flavors+1)
4277 : ! write(47,*) "cflavors",cflavors
4278 0 : string = '(1x,'//TRIM(ADJUSTL(cflavors))//'ES22.14)'
4279 0 : hstring = '(1x,'//TRIM(ADJUSTL(hcflavors))//'ES22.14)'
4280 : ! write(47,*) string
4281 0 : DO itime = 1, sp1
4282 0 : WRITE(43,string) DBLE(itime-1)*dt, &
4283 0 : & (/ ((op%Greens%oper(itime,iflavor,iflavorb), iflavorb=1, flavors),iflavor=1,flavors) /)
4284 0 : WRITE(143,hstring) DBLE(itime-1)*dt, &
4285 0 : & (/ ((real(op%Greens%oper(itime,iflavor,iflavorb)), iflavorb=1, flavors),iflavor=1,flavors) /)
4286 0 : WRITE(243,hstring) DBLE(itime-1)*dt, &
4287 0 : & (/ ((aimag(op%Greens%oper(itime,iflavor,iflavorb)), iflavorb=1, flavors),iflavor=1,flavors) /)
4288 : ! WRITE(44,*) DBLE(itime-1)*dt, &
4289 : ! & (op%Greens%oper(itime,iflavor,iflavor), iflavor=1, flavors)
4290 : ! WRITE(44,string) DBLE(itime-1)*dt, &
4291 : ! & (op%Greens%oper(itime,iflavor,iflavor), iflavor=1, flavors)
4292 : END DO
4293 0 : WRITE(43,*)
4294 0 : WRITE(143,*)
4295 0 : WRITE(243,*)
4296 : END IF
4297 : ! DO iflavor = 1, flavors
4298 : ! DO iflavor2 = 1, flavors
4299 : ! write(4436,*) "#",iflavor,iflavor2
4300 : ! do itime=1,sp1
4301 : ! write(4436,*) DBLE(itime-1)*dt,real(op%Greens%oper(itime,iflavor,iflavor2))
4302 : ! enddo
4303 : ! write(4436,*)
4304 : ! END DO
4305 : ! END DO
4306 : ! close(4436)
4307 :
4308 0 : IF ( .NOT. PRESENT(oFileIn) ) CLOSE(oFile)
4309 0 : CLOSE(43)
4310 0 : CLOSE(143)
4311 0 : CLOSE(243)
4312 : ! CLOSE(46)
4313 : ! CLOSE(47)
4314 : !call flush_unit(43)
4315 :
4316 : END SUBROUTINE CtqmcoffdiagComplex_printGreen
4317 : !!***
4318 :
4319 : !!****f* ABINIT/m_CtqmcoffdiagComplex/CtqmcoffdiagComplex_printD
4320 : !! NAME
4321 : !! CtqmcoffdiagComplex_printD
4322 : !!
4323 : !! FUNCTION
4324 : !! print individual double occupancy
4325 : !!
4326 : !! COPYRIGHT
4327 : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
4328 : !! This file is distributed under the terms of the
4329 : !! GNU General Public License, see ~abinit/COPYING
4330 : !! or http://www.gnu.org/copyleft/gpl.txt .
4331 : !!
4332 : !! INPUTS
4333 : !! op=ctqmc
4334 : !! oFileIn=file stream
4335 : !!
4336 : !! OUTPUT
4337 : !!
4338 : !! SIDE EFFECTS
4339 : !!
4340 : !! NOTES
4341 : !!
4342 : !! SOURCE
4343 :
4344 0 : SUBROUTINE CtqmcoffdiagComplex_printD(op,oFileIn)
4345 :
4346 : !Arguments ------------------------------------
4347 : TYPE(CtqmcoffdiagComplex) , INTENT(IN) :: op
4348 : INTEGER , OPTIONAL, INTENT(IN) :: oFileIn
4349 : !Local variables ------------------------------
4350 : INTEGER :: oFile
4351 : INTEGER :: iflavor1
4352 : INTEGER :: iflavor2
4353 :
4354 : !IF ( op%rank .NE. MOD(2,op%size)) RETURN
4355 0 : IF ( op%rank .NE. MOD(op%size+2,op%size)) RETURN
4356 :
4357 0 : oFile = 41
4358 0 : IF ( PRESENT(oFileIn) ) THEN
4359 0 : oFile = oFileIn
4360 : ELSE
4361 0 : OPEN(UNIT=oFile, FILE="D.dat")
4362 : END IF
4363 :
4364 0 : DO iflavor1 = 1, op%flavors
4365 0 : DO iflavor2 = iflavor1+1, op%flavors
4366 0 : WRITE(oFile,'(1x,A8,I4,A1,I4,A3,ES21.14)') "Orbitals", iflavor1, "-", iflavor2, " : ", op%measDE(iflavor2,iflavor1)
4367 : END DO
4368 : END DO
4369 :
4370 0 : IF ( .NOT. PRESENT(oFileIn) ) CLOSE(oFile)
4371 :
4372 : END SUBROUTINE CtqmcoffdiagComplex_printD
4373 : !!***
4374 :
4375 : !!****f* ABINIT/m_CtqmcoffdiagComplex/CtqmcoffdiagComplex_printE
4376 : !! NAME
4377 : !! CtqmcoffdiagComplex_printE
4378 : !!
4379 : !! FUNCTION
4380 : !! print energy and noise
4381 : !!
4382 : !! COPYRIGHT
4383 : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
4384 : !! This file is distributed under the terms of the
4385 : !! GNU General Public License, see ~abinit/COPYING
4386 : !! or http://www.gnu.org/copyleft/gpl.txt .
4387 : !!
4388 : !! INPUTS
4389 : !! op=ctqmc
4390 : !! oFileIn=file stream
4391 : !!
4392 : !! OUTPUT
4393 : !!
4394 : !! SIDE EFFECTS
4395 : !!
4396 : !! NOTES
4397 : !!
4398 : !! SOURCE
4399 :
4400 0 : SUBROUTINE CtqmcoffdiagComplex_printE(op,oFileIn)
4401 :
4402 : !Arguments ------------------------------------
4403 : TYPE(CtqmcoffdiagComplex) , INTENT(IN) :: op
4404 : INTEGER , OPTIONAL, INTENT(IN) :: oFileIn
4405 : !Local variables ------------------------------
4406 : INTEGER :: oFile
4407 : DOUBLE PRECISION :: E
4408 : DOUBLE PRECISION :: Noise
4409 :
4410 : !IF ( op%rank .NE. MOD(3,op%size)) RETURN
4411 0 : IF ( op%rank .NE. MOD(op%size+3,op%size)) RETURN
4412 :
4413 0 : oFile = 42
4414 0 : IF ( PRESENT(oFileIn) ) THEN
4415 0 : oFile = oFileIn
4416 : ELSE
4417 0 : OPEN(UNIT=oFile, FILE="BetaENoise.dat")
4418 : END IF
4419 :
4420 0 : CALL CtqmcoffdiagComplex_getE(op,E,Noise)
4421 :
4422 0 : WRITE(oFile,'(1x,F3.2,A2,ES21.14,A2,ES21.14)') op%beta, " ", E, " ", Noise
4423 :
4424 0 : IF ( .NOT. PRESENT(oFileIn) ) CLOSE(oFile)
4425 :
4426 : END SUBROUTINE CtqmcoffdiagComplex_printE
4427 : !!***
4428 :
4429 : !#ifdef CTCtqmcoffdiagComplex_ANALYSIS
4430 : !!****f* ABINIT/m_CtqmcoffdiagComplex/CtqmcoffdiagComplex_printPerturbation
4431 : !! NAME
4432 : !! CtqmcoffdiagComplex_printPerturbation
4433 : !!
4434 : !! FUNCTION
4435 : !! print perturbation order
4436 : !!
4437 : !! COPYRIGHT
4438 : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
4439 : !! This file is distributed under the terms of the
4440 : !! GNU General Public License, see ~abinit/COPYING
4441 : !! or http://www.gnu.org/copyleft/gpl.txt .
4442 : !!
4443 : !! INPUTS
4444 : !! op=ctqmc
4445 : !! oFileIn=file stream
4446 : !!
4447 : !! OUTPUT
4448 : !! argout(sizeout)=description
4449 : !!
4450 : !! SIDE EFFECTS
4451 : !!
4452 : !! NOTES
4453 : !!
4454 : !! SOURCE
4455 :
4456 0 : SUBROUTINE CtqmcoffdiagComplex_printPerturbation(op, oFileIn)
4457 :
4458 : !Arguments ------------------------------------
4459 : TYPE(CtqmcoffdiagComplex) , INTENT(IN) :: op
4460 : INTEGER , OPTIONAL, INTENT(IN) :: oFileIn
4461 : !Local variables-------------------------------
4462 : INTEGER :: oFile
4463 : INTEGER :: iorder
4464 : INTEGER :: order
4465 : INTEGER :: iflavor
4466 : INTEGER :: flavors
4467 : CHARACTER(LEN=2) :: a
4468 : CHARACTER(LEN=50) :: string
4469 :
4470 : !IF ( op%rank .NE. MOD(4,op%size)) RETURN
4471 0 : IF ( op%rank .NE. MOD(op%size+4,op%size)) RETURN
4472 0 : IF ( op%opt_order .LE. 0 ) RETURN
4473 :
4474 0 : oFile = 43
4475 0 : IF ( PRESENT(oFileIn) ) THEN
4476 0 : oFile = oFileIn
4477 : ELSE
4478 0 : OPEN(UNIT=oFile, FILE="Perturbation.dat")
4479 : END IF
4480 :
4481 0 : order = op%opt_order
4482 0 : flavors = op%flavors
4483 :
4484 0 : WRITE(a,'(I2)') flavors
4485 0 : string = '(I5,'//TRIM(ADJUSTL(a))//'F19.15)'
4486 0 : DO iorder = 1, order
4487 0 : WRITE(oFile,string) iorder-1, &
4488 0 : (/ (op%measPerturbation(iorder, iflavor), iflavor=1, flavors) /)
4489 : END DO
4490 :
4491 0 : IF ( .NOT. PRESENT(oFileIn) ) CLOSE(oFile)
4492 : END SUBROUTINE CtqmcoffdiagComplex_printPerturbation
4493 : !!***
4494 :
4495 : !!****f* ABINIT/m_CtqmcoffdiagComplex/CtqmcoffdiagComplex_printCorrelation
4496 : !! NAME
4497 : !! CtqmcoffdiagComplex_printCorrelation
4498 : !!
4499 : !! FUNCTION
4500 : !! print correlation fonctions
4501 : !!
4502 : !! COPYRIGHT
4503 : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
4504 : !! This file is distributed under the terms of the
4505 : !! GNU General Public License, see ~abinit/COPYING
4506 : !! or http://www.gnu.org/copyleft/gpl.txt .
4507 : !!
4508 : !! INPUTS
4509 : !! op=ctqmc
4510 : !! oFileIn=file stream
4511 : !!
4512 : !! OUTPUT
4513 : !!
4514 : !! SIDE EFFECTS
4515 : !!
4516 : !! NOTES
4517 : !!
4518 : !! SOURCE
4519 :
4520 0 : SUBROUTINE CtqmcoffdiagComplex_printCorrelation(op, oFileIn)
4521 :
4522 : !Arguments ------------------------------------
4523 : TYPE(CtqmcoffdiagComplex) , INTENT(IN) :: op
4524 : INTEGER , OPTIONAL, INTENT(IN) :: oFileIn
4525 : !Local variables ------------------------------
4526 : INTEGER :: oFile
4527 : INTEGER :: itime
4528 : INTEGER :: sp1
4529 : INTEGER :: iflavor
4530 : INTEGER :: i
4531 : INTEGER :: flavors
4532 : CHARACTER(LEN=2) :: a
4533 : CHARACTER(LEN=50) :: string
4534 : DOUBLE PRECISION :: dt
4535 :
4536 : !IF ( op%rank .NE. MOD(5,op%size)) RETURN
4537 0 : IF ( op%rank .NE. MOD(op%size+5,op%size)) RETURN
4538 0 : IF ( op%opt_analysis .NE. 1 ) RETURN
4539 :
4540 0 : oFile = 44
4541 0 : IF ( PRESENT(oFileIn) ) THEN
4542 0 : oFile = oFileIn
4543 : ELSE
4544 0 : OPEN(UNIT=oFile, FILE="Correlation.dat")
4545 : END IF
4546 :
4547 0 : sp1 = op%samples
4548 0 : dt = op%beta / sp1
4549 0 : sp1 = sp1 + 1
4550 0 : flavors = op%flavors
4551 :
4552 0 : i = 3*flavors + 1
4553 0 : WRITE(a,'(I2)') i
4554 0 : WRITE(oFile,*) "# time (/ (segement, antiseg, correl), i=1, flavor/)"
4555 0 : string = '(1x,'//TRIM(ADJUSTL(a))//'F19.15)'
4556 0 : DO itime = 1, sp1
4557 0 : WRITE(oFile,string) DBLE(itime-1)*dt, &
4558 : (/ ( &
4559 : (/ ( op%measCorrelation(itime, i, iflavor), i=1,3) /) &
4560 0 : , iflavor=1, flavors) /)
4561 : END DO
4562 :
4563 0 : IF ( .NOT. PRESENT(oFileIn) ) CLOSE(oFile)
4564 :
4565 : END SUBROUTINE CtqmcoffdiagComplex_printCorrelation
4566 : !!***
4567 : !#endif
4568 :
4569 : !!****f* ABINIT/m_CtqmcoffdiagComplex/CtqmcoffdiagComplex_printSpectra
4570 : !! NAME
4571 : !! CtqmcoffdiagComplex_printSpectra
4572 : !!
4573 : !! FUNCTION
4574 : !! print fourier transform of time evolution of number of electrons
4575 : !!
4576 : !! COPYRIGHT
4577 : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
4578 : !! This file is distributed under the terms of the
4579 : !! GNU General Public License, see ~abinit/COPYING
4580 : !! or http://www.gnu.org/copyleft/gpl.txt .
4581 : !!
4582 : !! INPUTS
4583 : !! op=ctqmc
4584 : !! oFileIn=file stream
4585 : !!
4586 : !! OUTPUT
4587 : !!
4588 : !! SIDE EFFECTS
4589 : !!
4590 : !! NOTES
4591 : !!
4592 : !! SOURCE
4593 :
4594 0 : SUBROUTINE CtqmcoffdiagComplex_printSpectra(op, oFileIn)
4595 :
4596 : !Arguments ------------------------------------
4597 : TYPE(CtqmcoffdiagComplex) , INTENT(IN) :: op
4598 : INTEGER , OPTIONAL, INTENT(IN) :: oFileIn
4599 : !Local variables ------------------------------
4600 : INTEGER :: oFile
4601 : INTEGER :: flavors
4602 : INTEGER :: indDensity
4603 : INTEGER :: endDensity
4604 : CHARACTER(LEN=4) :: a
4605 : CHARACTER(LEN=16) :: formatSpectra
4606 :
4607 : !IF ( op%rank .NE. MOD(6,op%size)) RETURN
4608 0 : IF ( op%opt_spectra .LT. 1 ) RETURN
4609 :
4610 0 : oFile = 45+op%rank
4611 0 : a ="0000"
4612 0 : WRITE(a,'(I4)') op%rank
4613 0 : IF ( PRESENT(oFileIn) ) THEN
4614 0 : oFile = oFileIn
4615 : ELSE
4616 0 : OPEN(UNIT=oFile, FILE="Markov_"//TRIM(ADJUSTL(a))//".dat")
4617 : END IF
4618 :
4619 0 : flavors = op%flavors
4620 0 : WRITE(a,'(I4)') flavors+1
4621 0 : formatSpectra ='(1x,'//TRIM(ADJUSTL(a))//'ES22.14)'
4622 0 : WRITE(oFile,*) "# freq[/hermalization] FFT"
4623 :
4624 : !PROBLEM eos_gnu_13.2_mpich . %endDensity was introduced throughout
4625 0 : endDensity = SIZE(op%density,2)
4626 : ! endDensity=op%endDensity
4627 : !ENDPROBLEM
4628 0 : DO WHILE ( op%density(flavors+1,endDensity) .EQ. -1 )
4629 0 : endDensity = endDensity -1
4630 : END DO
4631 :
4632 0 : DO indDensity = 1, endDensity
4633 0 : WRITE(oFile,formatSpectra) op%density(flavors+1,indDensity), op%density(1:flavors,indDensity)
4634 : END DO
4635 :
4636 0 : IF ( .NOT. PRESENT(oFileIn) ) CLOSE(oFile)
4637 :
4638 : END SUBROUTINE CtqmcoffdiagComplex_printSpectra
4639 : !!***
4640 :
4641 : !!****f* ABINIT/m_CtqmcoffdiagComplex/CtqmcoffdiagComplex_destroy
4642 : !! NAME
4643 : !! CtqmcoffdiagComplex_destroy
4644 : !!
4645 : !! FUNCTION
4646 : !! destroy and deallocate all variables
4647 : !!
4648 : !! COPYRIGHT
4649 : !! Copyright (C) 2013-2025 ABINIT group (J. Bieder)
4650 : !! This file is distributed under the terms of the
4651 : !! GNU General Public License, see ~abinit/COPYING
4652 : !! or http://www.gnu.org/copyleft/gpl.txt .
4653 : !!
4654 : !! INPUTS
4655 : !! op=ctqmc
4656 : !!
4657 : !! OUTPUT
4658 : !!
4659 : !! SIDE EFFECTS
4660 : !!
4661 : !! NOTES
4662 : !!
4663 : !! SOURCE
4664 :
4665 0 : SUBROUTINE CtqmcoffdiagComplex_destroy(op)
4666 :
4667 : !Arguments ------------------------------------
4668 : TYPE(CtqmcoffdiagComplex), INTENT(INOUT) :: op
4669 : !Local variables ------------------------------
4670 : !INTEGER :: iflavor
4671 : INTEGER :: flavors
4672 : INTEGER :: i
4673 : INTEGER :: j
4674 : INTEGER :: k
4675 :
4676 0 : flavors = op%flavors
4677 :
4678 0 : CALL ImpurityOperator_destroy(op%Impurity)
4679 0 : CALL BathOperatoroffdiagComplex_destroy(op%Bath)
4680 0 : CALL Vector_destroy(op%measNoise(1))
4681 0 : CALL Vector_destroy(op%measNoise(2))
4682 :
4683 : !sui!write(6,*) "before greenhyb_destroy in ctmqc_destroy"
4684 0 : CALL GreenHyboffdiagComplex_destroy(op%Greens)
4685 : !#ifdef CTCtqmcoffdiagComplex_ANALYSIS
4686 0 : FREEIF(op%measCorrelation)
4687 0 : FREEIF(op%measPerturbation)
4688 0 : FREEIF(op%meas_fullemptylines)
4689 0 : FREEIF(op%measN)
4690 0 : IF ( op%opt_histo .GT. 0 ) THEN
4691 0 : FREEIF(op%occup_histo_time)
4692 0 : FREEIF(op%suscep)
4693 0 : FREEIF(op%occupconfig)
4694 0 : FREEIF(op%chi)
4695 0 : FREEIF(op%chicharge)
4696 0 : FREEIF(op%ntot)
4697 : END IF
4698 0 : FREEIF(op%measDE)
4699 0 : FREEIF(op%mu)
4700 0 : FREEIF(op%hybri_limit)
4701 0 : FREEIF(op%abNoiseG)
4702 0 : IF ( ALLOCATED(op%measNoiseG) ) THEN
4703 0 : DO i=1,2
4704 0 : DO j = 1, op%flavors
4705 0 : DO k= 1, op%samples+1
4706 0 : CALL Vector_destroy(op%measNoiseG(k,j,i))
4707 : END DO
4708 : END DO
4709 : END DO
4710 0 : DT_FREE(op%measNoiseG)
4711 : END IF
4712 0 : FREEIF(op%density)
4713 : !#endif
4714 0 : op%ostream = 0
4715 0 : op%istream = 0
4716 :
4717 0 : op%sweeps = 0
4718 0 : op%thermalization = 0
4719 0 : op%flavors = 0
4720 : ! op%endDensity = 0
4721 0 : op%samples = 0
4722 0 : op%beta = 0.d0
4723 : ! op%seg_added = 0.d0
4724 : ! op%anti_added = 0.d0
4725 : ! op%seg_removed = 0.d0
4726 : ! op%anti_removed = 0.d0
4727 : ! op%seg_sign = 0.d0
4728 : ! op%anti_sign = 0.d0
4729 0 : op%stats = 0.d0
4730 0 : op%swap = 0.d0
4731 :
4732 :
4733 0 : op%set = .FALSE.
4734 0 : op%done = .FALSE.
4735 0 : op%init = .FALSE.
4736 0 : END SUBROUTINE CtqmcoffdiagComplex_destroy
4737 : !!***
4738 :
4739 : !!****f* ABINIT/m_CtqmcoffdiagComplex/CtqmcoffdiagComplex_setMagmom
4740 : !! NAME
4741 : !! CtqmcoffdiagComplex_setMagmom
4742 : !!
4743 : !! FUNCTION
4744 : !! set the Magnetic moment matrix for susceptibility
4745 : !!
4746 : !! COPYRIGHT
4747 : !! Copyright (C) 2013-2025 ABINIT group (F. Gendron)
4748 : !! This file is distributed under the terms of the
4749 : !! GNU General Public License, see ~abinit/COPYING
4750 : !! or http://www.gnu.org/copyleft/gpl.txt .
4751 : !!
4752 : !! INPUTS
4753 : !!
4754 : !! OUTPUT
4755 : !!
4756 : !! SIDE EFFECTS
4757 : !!
4758 : !! NOTES
4759 : !!
4760 : !! PARENTS
4761 : !! Will be filled automatically by the parent script
4762 : !!
4763 : !! CHILDREN
4764 : !! Will be filled automatically by the parent script
4765 : !!
4766 : !! SOURCE
4767 :
4768 0 : SUBROUTINE CtqmcoffdiagComplex_setMagmom(op,Magmom_orb, Magmom_spin, Magmom_tot)
4769 :
4770 : !Arguments ------------------------------------
4771 : TYPE(CtqmcoffdiagComplex), INTENT(INOUT) ::op
4772 : DOUBLE PRECISION, DIMENSION(:,:), INTENT(IN) :: Magmom_orb
4773 : DOUBLE PRECISION, DIMENSION(:,:), INTENT(IN) :: Magmom_spin
4774 : DOUBLE PRECISION, DIMENSION(:,:), INTENT(IN) :: Magmom_tot
4775 : !Local variables ------------------------------
4776 : ! INTEGER :: iflavor1,iflavor2
4777 :
4778 : ! do iflavor1=1,10
4779 : ! do iflavor2=1,10
4780 : ! if(iflavor1==iflavor2) THEN
4781 : ! write(6,*) iflavor1, iflavor2, Magmom(iflavor1,iflavor2)
4782 : ! end if
4783 : ! end do
4784 : ! end do
4785 :
4786 0 : IF ( SIZE(Magmom_orb) .NE. op%flavors*op%flavors ) &
4787 0 : CALL ERROR("CtqmcoffdiagComplex_setU : Wrong Magnetic Moment matrix (size) ")
4788 :
4789 0 : CALL ImpurityOperator_setMagmommat(op%Impurity, Magmom_orb, Magmom_spin, Magmom_tot)
4790 :
4791 0 : END SUBROUTINE CtqmcoffdiagComplex_setMagmom
4792 : !!***
4793 :
4794 0 : END MODULE m_CtqmcoffdiagComplex
4795 : !!***
4796 :
4797 : #undef CTQMC_SLICE1
4798 : #undef CTQMC_SLICE2
4799 : #undef CTQMC_SEGME
4800 : #undef CTQMC_ANTIS
4801 : #undef CTQMC_ADDED
4802 : #undef CTQMC_REMOV
4803 : #undef CTQMC_DETSI
|
