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