Line data Source code
1 :
2 : #if defined HAVE_CONFIG_H
3 : #include "config.h"
4 : #endif
5 :
6 : !!****m* ABINIT/m_Ctqmc
7 : !! NAME
8 : !! m_Ctqmc
9 : !!
10 : !! FUNCTION
11 : !! Manage and drive all the CTQMC
12 : !! Should not be used if you don't know what you do
13 : !! Please use CtqmcInterface
14 : !!
15 : !! COPYRIGHT
16 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
17 : !! This file is distributed under the terms of the
18 : !! GNU General Public License, see ~abinit/COPYING
19 : !! or http://www.gnu.org/copyleft/gpl.txt .
20 : !!
21 : !! NOTES
22 : !!
23 : !! SOURCE
24 :
25 : #include "defs.h"
26 :
27 : MODULE m_Ctqmc
28 :
29 : USE m_Global
30 : USE m_GreenHyb
31 : USE m_BathOperator
32 : USE m_ImpurityOperator
33 : USE m_Stat
34 : USE m_FFTHyb
35 : USE m_OurRng
36 : USE m_Vector
37 : USE m_errors
38 : use m_io_tools, only : open_file
39 : use defs_basis
40 :
41 : USE_MPI
42 :
43 : IMPLICIT NONE
44 :
45 : !!***
46 :
47 : PRIVATE
48 :
49 : INTEGER, PARAMETER :: CTQMC_SLICE1 = 100
50 : ! Coupe Sweeps en 100
51 : INTEGER, PARAMETER :: CTQMC_SLICE2 = 100
52 : ! Coupe modNoise1 en 100
53 : INTEGER, PARAMETER :: CTQMC_SEGME = 1
54 : INTEGER, PARAMETER :: CTQMC_ANTIS = -2
55 : INTEGER, PARAMETER :: CTQMC_ADDED = 3
56 : INTEGER, PARAMETER :: CTQMC_REMOV = 4
57 : INTEGER, PARAMETER :: CTQMC_DETSI = 5
58 :
59 : !!****t* m_Ctqmc/Ctqmc
60 : !! NAME
61 : !! Ctqmc
62 : !!
63 : !! FUNCTION
64 : !! This structured datatype contains the necessary data
65 : !!
66 : !! COPYRIGHT
67 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
68 : !! This file is distributed under the terms of the
69 : !! GNU General Public License, see ~abinit/COPYING
70 : !! or http://www.gnu.org/copyleft/gpl.txt .
71 : !!
72 : !! SOURCE
73 :
74 : TYPE, PUBLIC :: Ctqmc
75 :
76 : LOGICAL _PRIVATE :: init = .FALSE.
77 : ! Flag: is MC initialized
78 :
79 : LOGICAL _PRIVATE :: set = .FALSE.
80 : ! Flag: ??
81 :
82 : LOGICAL _PRIVATE :: setU = .FALSE.
83 : ! Flag: is U Set ?
84 :
85 : LOGICAL _PRIVATE :: inF = .FALSE.
86 : ! Flag: is hybridization fct in input ?
87 :
88 : LOGICAL :: done = .FALSE.
89 : ! Flag: is MC terminated ? (PUBLIC for parallel)
90 :
91 : LOGICAL :: para = .FALSE.
92 : ! Flag: do we have parameters in input (PUBLIC for parallel)
93 :
94 : LOGICAL _PRIVATE :: have_MPI = .FALSE.
95 : ! Flag:
96 :
97 : INTEGER :: opt_movie = 0
98 : ! (PUBLIC for parallel)
99 :
100 : INTEGER :: opt_analysis = 0
101 : ! correlations (PUBLIC for parallel)
102 :
103 : INTEGER :: opt_check = 0
104 : ! various check 0 (PUBLIC for parallel)
105 : ! various check 1 impurity
106 : ! various check 2 bath
107 : ! various check 3 both
108 :
109 : INTEGER :: opt_order = 0
110 : ! nb of segments max for analysis (PUBLIC for parallel)
111 :
112 : INTEGER :: opt_histo = 0
113 : ! Enable histograms (PUBLIC for parallel)
114 :
115 : INTEGER :: opt_noise = 0
116 : ! compute noise (PUBLIC for parallel)
117 :
118 : INTEGER :: opt_spectra = 0
119 : ! markov chain FT (correlation time) (PUBLIC for parallel)
120 :
121 : INTEGER :: opt_levels = 0
122 : ! do we have energy levels (PUBLIC for parallel)
123 :
124 : INTEGER :: flavors
125 : ! (PUBLIC for parallel)
126 : !
127 : INTEGER :: endDensity
128 : ! (PUBLIC for parallel)
129 : !
130 : INTEGER :: nspinor
131 : ! (PUBLIC for parallel)
132 : !
133 : INTEGER :: measurements
134 : ! nb of measure in the MC (PUBLIC for parallel)
135 :
136 : INTEGER :: samples
137 : ! nb of L points (PUBLIC for parallel)
138 :
139 : INTEGER(8) :: seed
140 : ! (PUBLIC for parallel - critical for thread-local RNG)
141 :
142 : INTEGER :: sweeps
143 : ! (PUBLIC for parallel)
144 :
145 : INTEGER :: thermalization
146 : ! (PUBLIC for parallel)
147 :
148 : INTEGER :: ostream
149 : ! output file (PUBLIC for parallel)
150 :
151 : INTEGER :: istream
152 : ! input file (PUBLIC for parallel)
153 :
154 : INTEGER :: modNoise1
155 : ! measure the noise each modNoise1 (PUBLIC for parallel)
156 :
157 : INTEGER :: modNoise2
158 : ! measure the noise each modNoise2 (PUBLIC for parallel)
159 :
160 : INTEGER _PRIVATE :: activeFlavor
161 : ! orbital on which one do sth now
162 :
163 : INTEGER, DIMENSION(1:2) :: modGlobalMove
164 : ! 1: gloabl move each modglobalmove(1) (PUBLIC for parallel)
165 : ! 2: we have done modglobalmove(2) for two different orbitals.
166 :
167 : INTEGER :: Wmax
168 : ! Max freq for FT (PUBLIC for parallel)
169 :
170 : DOUBLE PRECISION, DIMENSION(1:6) :: stats
171 : ! to now how many negative determinant, antisegments,seeme.e.twfs...j (PUBLIC for parallel)
172 :
173 : DOUBLE PRECISION :: swap
174 : ! nb of successfull GM (PUBLIC for parallel)
175 :
176 : INTEGER :: MY_COMM
177 : ! (PUBLIC for parallel)
178 :
179 : INTEGER :: tid = 1
180 : ! (PUBLIC for parallel)
181 :
182 : INTEGER :: rank
183 : ! (PUBLIC for parallel)
184 :
185 : INTEGER :: size
186 : ! size of MY_COMM (PUBLIC for parallel)
187 :
188 : DOUBLE PRECISION _PRIVATE :: runTime ! time for the run routine
189 : !
190 :
191 : DOUBLE PRECISION :: beta
192 : ! (PUBLIC for parallel)
193 :
194 : DOUBLE PRECISION :: U
195 : ! (PUBLIC for parallel)
196 :
197 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: mu
198 : ! levels (PUBLIC for parallel)
199 :
200 : TYPE(GreenHyb) , ALLOCATABLE, DIMENSION(: ) :: Greens
201 : ! (PUBLIC for parallel)
202 :
203 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,: ) :: measN
204 : ! measure of occupations (3or4,flavor) (PUBLIC for parallel)
205 :
206 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,: ) :: measDE
207 : ! (flavor,flavor) double occupancies (PUBLIC for parallel)
208 : ! (1,1): total energy of correlation.
209 :
210 : DOUBLE PRECISION _PRIVATE :: a_Noise
211 : ! Noise a exp (-bx) for the noise
212 :
213 : DOUBLE PRECISION _PRIVATE :: b_Noise
214 : ! Noise a exp (-bx) for the noise
215 :
216 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) _PRIVATE :: abNoiseG !(ab,tau,flavor)
217 : ! Noise but for G
218 :
219 : TYPE(Vector) , DIMENSION(1:2) _PRIVATE :: measNoise
220 : TYPE(Vector), ALLOCATABLE, DIMENSION(:,:,:) _PRIVATE :: measNoiseG !(tau,flavor,mod)
221 : ! accumulate each value relataed to measurenoise 1 2
222 :
223 : DOUBLE PRECISION _PRIVATE :: inv_dt
224 : ! 1/(beta/L)
225 :
226 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,: ) _PRIVATE :: measPerturbation
227 : ! opt_order,nflavor
228 :
229 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: occup_histo_time
230 : ! nflavor
231 :
232 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: occupconfig
233 : ! 2**nflavor
234 :
235 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: suscep
236 : ! samples
237 :
238 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: chi
239 : ! samples
240 :
241 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: chicharge
242 : ! samples
243 :
244 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: ntot
245 : ! occupation total, t2g, eg
246 :
247 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) _PRIVATE :: measCorrelation
248 : ! segment,antisegment,nflavor,nflavor
249 :
250 : DOUBLE PRECISION _PRIVATE :: errorImpurity
251 : ! check
252 :
253 : DOUBLE PRECISION _PRIVATE :: errorBath
254 : ! for check
255 :
256 : TYPE(BathOperator) :: Bath
257 : ! (PUBLIC for parallel - needed for deep copy)
258 :
259 : TYPE(ImpurityOperator) :: Impurity
260 : ! (PUBLIC for parallel - needed for deep copy)
261 :
262 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) _PRIVATE :: density
263 :
264 : END TYPE Ctqmc
265 : !!***
266 :
267 : PUBLIC :: Ctqmc_init
268 : PRIVATE :: Ctqmc_setParameters
269 : PUBLIC :: Ctqmc_setSweeps
270 : PUBLIC :: Ctqmc_setSeed
271 : PRIVATE :: Ctqmc_allocateAll
272 : PRIVATE :: Ctqmc_allocateOpt
273 : PUBLIC :: Ctqmc_setG0wFile
274 : PUBLIC :: Ctqmc_setG0wTab
275 : PUBLIC :: Ctqmc_setU
276 : PRIVATE :: Ctqmc_clear
277 : PUBLIC :: Ctqmc_reset
278 : PUBLIC :: Ctqmc_setMu
279 : PRIVATE :: Ctqmc_computeF
280 : PUBLIC :: Ctqmc_run
281 : PRIVATE :: Ctqmc_loop
282 : PRIVATE :: Ctqmc_tryAddRemove
283 : PRIVATE :: Ctqmc_trySwap
284 : PRIVATE :: Ctqmc_measN
285 : PRIVATE :: Ctqmc_measCorrelation
286 : PRIVATE :: Ctqmc_measPerturbation
287 : PUBLIC :: Ctqmc_getResult
288 : PUBLIC :: Ctqmc_symmetrizeGreen
289 : PUBLIC :: Ctqmc_getGreen
290 : PUBLIC :: Ctqmc_getD
291 : PUBLIC :: Ctqmc_getE
292 : PUBLIC :: Ctqmc_printAll
293 : PUBLIC :: Ctqmc_printQMC
294 : PUBLIC :: Ctqmc_printGreen
295 : PUBLIC :: Ctqmc_printD
296 : PUBLIC :: Ctqmc_printE
297 : PUBLIC :: Ctqmc_printPerturbation
298 : PUBLIC :: Ctqmc_printCorrelation
299 : PUBLIC :: Ctqmc_printSpectra
300 : PUBLIC :: Ctqmc_destroy
301 : PUBLIC :: Ctqmc_setMagmom
302 :
303 : CONTAINS
304 : !!***
305 :
306 : !!****f* ABINIT/m_Ctqmc/Ctqmc_init
307 : !! NAME
308 : !! Ctqmc_init
309 : !!
310 : !! FUNCTION
311 : !! Initialize the type Ctqmc
312 : !! Allocate all the non optional variables
313 : !!
314 : !! COPYRIGHT
315 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
316 : !! This file is distributed under the terms of the
317 : !! GNU General Public License, see ~abinit/COPYING
318 : !! or http://www.gnu.org/copyleft/gpl.txt .
319 : !!
320 : !! INPUTS
321 : !! this=ctqmc
322 : !! ostream=where to write
323 : !! istream=where to read the input parameters if so
324 : !! bFile=do we read in istream ?
325 : !! MY_COMM=mpi communicator for the CTQMC
326 : !! iBuffer=input parameters if bFile is false
327 : !!
328 : !! OUTPUT
329 : !!
330 : !! SIDE EFFECTS
331 : !!
332 : !! NOTES
333 : !!
334 : !! SOURCE
335 :
336 102 : SUBROUTINE Ctqmc_init(this, ostream, istream, bFile, MY_COMM, iBuffer)
337 :
338 :
339 : #ifdef HAVE_MPI1
340 : include 'mpif.h'
341 : #endif
342 : !Arguments ------------------------------------
343 : TYPE(Ctqmc), INTENT(INOUT) :: this
344 : INTEGER , INTENT(IN ) :: ostream
345 : INTEGER , INTENT(IN ) :: istream
346 : LOGICAL , INTENT(IN ) :: bFile
347 : DOUBLE PRECISION, DIMENSION(1:10), OPTIONAL, INTENT(IN) :: iBuffer
348 : INTEGER , OPTIONAL, INTENT(IN ) :: MY_COMM
349 : !Local variables ------------------------------
350 : #ifdef HAVE_MPI
351 : INTEGER :: ierr
352 : #endif
353 : INTEGER :: iflavor
354 : #ifdef __GFORTRAN__
355 : ! INTEGER :: pid
356 : ! CHARACTER(LEN=5) :: Cpid
357 : !
358 : #endif
359 : DOUBLE PRECISION, DIMENSION(1:10) :: buffer
360 :
361 102 : this%ostream = ostream
362 102 : this%istream = istream
363 102 : this%tid = 1
364 :
365 : ! --- RENICE ---
366 : !#ifdef __GFORTRAN__
367 : ! pid = GetPid()
368 : ! WRITE(Cpid,'(I5)') pid
369 : ! CALL SYSTEM('renice +19 '//TRIM(ADJUSTL(Cpid))//' > /dev/null')
370 : !#endif
371 : !! --- RENICE ---
372 :
373 102 : IF ( PRESENT(MY_COMM)) THEN
374 : #ifdef HAVE_MPI
375 102 : this%have_MPI = .TRUE.
376 102 : this%MY_COMM = MY_COMM
377 102 : CALL MPI_Comm_rank(this%MY_COMM, this%rank, ierr)
378 102 : CALL MPI_Comm_size(this%MY_COMM, this%size, ierr)
379 : #else
380 : CALL WARN("Ctqmc_init : MPI is not used ")
381 : this%have_MPI = .FALSE.
382 : this%MY_COMM = -1
383 : this%rank = 0
384 : this%size = 1
385 : #endif
386 : ELSE
387 0 : this%have_MPI = .FALSE.
388 0 : this%MY_COMM = -1
389 0 : this%rank = 0
390 0 : this%size = 1
391 : END IF
392 :
393 : !IF ( this%rank .EQ. 0 ) THEN
394 : ! WRITE(ostream,'(A20)') 'Job reniced with +19'
395 : !CALL FLUSH(ostream)
396 : !END IF
397 :
398 102 : IF ( bFile .EQV. .TRUE. ) THEN
399 0 : IF ( this%rank .EQ. 0 ) THEN
400 :
401 0 : READ(istream,*) buffer(1) !iseed
402 0 : READ(istream,*) buffer(2) !this%sweeps
403 0 : READ(istream,*) buffer(3) !this%thermalization
404 0 : READ(istream,*) buffer(4) !this%measurements
405 0 : READ(istream,*) buffer(5) !this%flavors
406 0 : READ(istream,*) buffer(6) !this%samples
407 0 : READ(istream,*) buffer(7) !this%beta
408 0 : READ(istream,*) buffer(8) !U
409 0 : READ(istream,*) buffer(9) !iTech
410 0 : READ(istream,*) buffer(10)!this%nspinor
411 : !READ(istream,*) buffer(9) !Wmax
412 : !#ifdef CTCtqmc_ANALYSIS
413 : !READ(istream,*) buffer(10) !order
414 : !#endif
415 : END IF
416 :
417 : #ifdef HAVE_MPI
418 0 : IF ( this%have_MPI .EQV. .TRUE. ) &
419 : CALL MPI_Bcast(buffer, 10, MPI_DOUBLE_PRECISION, 0, &
420 0 : this%MY_COMM, ierr)
421 : #endif
422 102 : ELSE IF ( PRESENT(iBuffer) ) THEN
423 102 : buffer(1:10) = iBuffer(1:10)
424 : ELSE
425 0 : CALL ERROR("Ctqmc_init : No input parameters ")
426 : END IF
427 :
428 102 : CALL Ctqmc_setParameters(this, buffer)
429 :
430 102 : CALL Ctqmc_allocateAll(this)
431 :
432 1178 : DO iflavor = 1, this%flavors
433 1178 : CALL GreenHyb_init(this%Greens(iflavor),this%samples, this%beta, iTech=INT(buffer(9)),MY_COMM=this%MY_COMM)
434 : END DO
435 :
436 :
437 : ! this%seg_added = 0.d0
438 : ! this%anti_added = 0.d0
439 : ! this%seg_removed = 0.d0
440 : ! this%anti_removed = 0.d0
441 : ! this%seg_sign = 0.d0
442 : ! this%anti_sign = 0.d0
443 714 : this%stats(:) = 0.d0
444 102 : this%swap = 0.d0
445 102 : this%runTime = 0.d0
446 :
447 102 : CALL Vector_init(this%measNoise(1),this%sweeps/this%modNoise1)
448 102 : CALL Vector_init(this%measNoise(2),(this%sweeps/this%modNoise1+1)*CTQMC_SLICE2)
449 : !CALL Vector_init(this%measNoise(3),101)
450 : !CALL Vector_init(this%measNoise(4),101)
451 :
452 102 : this%set = this%para .AND. this%inF
453 102 : this%done = .FALSE.
454 102 : this%init = .TRUE.
455 :
456 : !#ifdef CTCtqmc_CHECK
457 102 : this%errorImpurity = 0.d0
458 102 : this%errorBath = 0.d0
459 : !#endif
460 102 : END SUBROUTINE Ctqmc_init
461 : !!***
462 :
463 : !!****f* ABINIT/m_Ctqmc/Ctqmc_setParameters
464 : !! NAME
465 : !! Ctqmc_setParameters
466 : !!
467 : !! FUNCTION
468 : !! set all parameters and operators
469 : !!
470 : !! COPYRIGHT
471 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
472 : !! This file is distributed under the terms of the
473 : !! GNU General Public License, see ~abinit/COPYING
474 : !! or http://www.gnu.org/copyleft/gpl.txt .
475 : !!
476 : !! INPUTS
477 : !! this=ctqmc
478 : !! buffer=input parameters
479 : !!
480 : !! OUTPUT
481 : !!
482 : !! SIDE EFFECTS
483 : !!
484 : !! NOTES
485 : !!
486 : !! SOURCE
487 :
488 102 : SUBROUTINE Ctqmc_setParameters(this,buffer)
489 :
490 : !Arguments ------------------------------------
491 : TYPE(Ctqmc), INTENT(INOUT) :: this
492 : DOUBLE PRECISION, DIMENSION(1:10), INTENT(IN ) :: buffer
493 :
494 :
495 102 : this%thermalization = INT(buffer(3)) !this%thermalization
496 102 : CALL Ctqmc_setSeed(this,INT(buffer(1)))
497 102 : CALL Ctqmc_setSweeps(this,buffer(2))
498 :
499 102 : this%measurements = INT(buffer(4)) !this%measurements
500 102 : this%flavors = INT(buffer(5)) !this%flavors
501 102 : this%samples = INT(buffer(6)) !this%samples
502 102 : this%beta = buffer(7) !this%beta
503 102 : this%U = buffer(8) !U
504 102 : this%nspinor = INT(buffer(10))!this%nspinor
505 : ! this%mu = buffer(9) !this%mu
506 : !this%Wmax = INT(buffer(9)) !Freq
507 : !#ifdef CTCtqmc_ANALYSIS
508 : ! this%order = INT(buffer(10)) ! order
509 102 : this%inv_dt = this%samples / this%beta
510 : !#endif
511 :
512 : !CALL ImpurityOperator_init(this%Impurity,this%flavors,this%beta, this%samples)
513 102 : CALL ImpurityOperator_init(this%Impurity,this%flavors,this%beta)
514 102 : IF ( this%U .GE. 0.d0 ) THEN
515 102 : CALL ImpurityOperator_computeU(this%Impurity,this%U,0.d0)
516 102 : this%setU = .TRUE.
517 : END IF
518 : ! this%mu = this%mu + this%Impurity%shift_mu
519 :
520 102 : CALL BathOperator_init(this%Bath, this%flavors, this%samples, this%beta, INT(buffer(9)))
521 :
522 102 : this%para = .TRUE.
523 :
524 102 : END SUBROUTINE Ctqmc_setParameters
525 : !!***
526 :
527 : !!****f* ABINIT/m_Ctqmc/Ctqmc_setSweeps
528 : !! NAME
529 : !! Ctqmc_setSweeps
530 : !!
531 : !! FUNCTION
532 : !! set the number of sweeps
533 : !!
534 : !! COPYRIGHT
535 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
536 : !! This file is distributed under the terms of the
537 : !! GNU General Public License, see ~abinit/COPYING
538 : !! or http://www.gnu.org/copyleft/gpl.txt .
539 : !!
540 : !! INPUTS
541 : !! this=ctqmc
542 : !! sweeps=asked sweeps
543 : !!
544 : !! OUTPUT
545 : !!
546 : !! SIDE EFFECTS
547 : !!
548 : !! NOTES
549 : !!
550 : !! SOURCE
551 :
552 151 : SUBROUTINE Ctqmc_setSweeps(this,sweeps)
553 :
554 : !Arguments ------------------------------------
555 : TYPE(Ctqmc) , INTENT(INOUT) :: this
556 : DOUBLE PRECISION , INTENT(IN ) :: sweeps
557 :
558 151 : this%sweeps = NINT(sweeps / DBLE(this%size))
559 : ! write(6,*) this%sweeps,NINT(sweeps / DBLE(this%size)),ANINT(sweeps/DBLE(this%size))
560 151 : IF ( DBLE(this%sweeps) .NE. ANINT(sweeps/DBLE(this%size)) ) &
561 0 : CALL ERROR("Ctqmc_setSweeps : sweeps is negative or too big ")
562 151 : IF ( this%sweeps .LT. 2*CTQMC_SLICE1 ) THEN !202
563 0 : CALL WARNALL("Ctqmc_setSweeps : # sweeps automtically changed ")
564 0 : this%sweeps = 2*CTQMC_SLICE1
565 : ! ELSE IF ( this%sweeps .LT. this%thermalization ) THEN
566 : ! CALL WARNALL("Ctqmc_setSweeps : Thermalization > sweeps / cpu -> auto fix")
567 : ! this%sweeps = this%thermalization
568 : END IF
569 151 : IF ( DBLE(NINT(DBLE(this%sweeps)*DBLE(this%size)/DBLE(CTQMC_SLICE1))) .NE. &
570 151 : ANINT(DBLE(this%sweeps)*DBLE(this%size)/DBLE(CTQMC_SLICE1)) ) THEN
571 0 : this%modNoise1 = this%sweeps
572 : ELSE
573 151 : this%modNoise1 = MIN(this%sweeps,INT(DBLE(this%sweeps)*DBLE(this%size) / DBLE(CTQMC_SLICE1))) !101
574 : END IF
575 151 : this%modNoise2 = MAX(this%modNoise1 / CTQMC_SLICE2, 1) ! 100
576 : ! this%modGlobalMove(1) = this%thermalization / 10 + 1
577 : ! this%modGlobalMove(2) = 0
578 :
579 151 : END SUBROUTINE Ctqmc_setSweeps
580 : !!***
581 :
582 : !!****f* ABINIT/m_Ctqmc/Ctqmc_setSeed
583 : !! NAME
584 : !! Ctqmc_setSeed
585 : !!
586 : !! FUNCTION
587 : !! initialize random number generator
588 : !!
589 : !! COPYRIGHT
590 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
591 : !! This file is distributed under the terms of the
592 : !! GNU General Public License, see ~abinit/COPYING
593 : !! or http://www.gnu.org/copyleft/gpl.txt .
594 : !!
595 : !! INPUTS
596 : !! this=ctqmc
597 : !! iseed=seed from imput
598 : !!
599 : !! OUTPUT
600 : !!
601 : !! SIDE EFFECTS
602 : !!
603 : !! NOTES
604 : !!
605 : !! SOURCE
606 :
607 102 : SUBROUTINE Ctqmc_setSeed(this,iseed)
608 :
609 : !Arguments ------------------------------------
610 : TYPE(Ctqmc), INTENT(INOUT) :: this
611 : INTEGER , INTENT(IN ) :: iseed
612 : !Local variables ------------------------------
613 : !INTEGER :: n
614 : !INTEGER :: i
615 : !INTEGER, DIMENSION(:), ALLOCATABLE :: seed
616 :
617 :
618 : !CALL RANDOM_SEED(size = n)
619 : !MALLOC(seed,(n))
620 : !seed = iseed + (/ (i - 1, i = 1, n) /)
621 :
622 : !CALL RANDOM_SEED(PUT = seed+this%rank)
623 :
624 : !FREE(seed)
625 :
626 102 : this%seed=INT(iseed+this%rank,8)
627 :
628 0 : END SUBROUTINE Ctqmc_setSeed
629 : !!***
630 :
631 : !!****f* ABINIT/m_Ctqmc/Ctqmc_allocateAll
632 : !! NAME
633 : !! Ctqmc_allocateAll
634 : !!
635 : !! FUNCTION
636 : !! Allocate all non option varibales
637 : !!
638 : !! COPYRIGHT
639 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
640 : !! This file is distributed under the terms of the
641 : !! GNU General Public License, see ~abinit/COPYING
642 : !! or http://www.gnu.org/copyleft/gpl.txt .
643 : !!
644 : !! INPUTS
645 : !! this=ctqmc
646 : !!
647 : !! OUTPUT
648 : !!
649 : !! SIDE EFFECTS
650 : !!
651 : !! NOTES
652 : !!
653 : !! SOURCE
654 :
655 102 : SUBROUTINE Ctqmc_allocateAll(this)
656 :
657 : !Arguments ------------------------------------
658 : TYPE(Ctqmc), INTENT(INOUT) :: this
659 : !Local variables ------------------------------
660 : INTEGER :: flavors
661 :
662 102 : IF ( .NOT. this%para ) &
663 0 : CALL ERROR("Ctqmc_allocateAll : Ctqmc_setParameters never called ")
664 :
665 102 : flavors = this%flavors
666 :
667 102 : DT_FREEIF(this%Greens)
668 1382 : DT_MALLOC(this%Greens,(1:flavors))
669 :
670 102 : FREEIF(this%measN)
671 306 : MALLOC(this%measN,(1:4,1:flavors))
672 5482 : this%measN = 0.d0
673 :
674 102 : FREEIF(this%measDE)
675 408 : MALLOC(this%measDE,(1:flavors,1:flavors) )
676 13170 : this%measDE = 0.d0
677 :
678 102 : FREEIF(this%mu)
679 : #ifdef FC_LLVM
680 : ! LLVM 16 doesn't recognize this macro here
681 : MALLOC(this%mu, (1:flavors) )
682 : #else
683 306 : MALLOC(this%mu, (1:flavors))
684 : #endif
685 1178 : this%mu = 0.d0
686 102 : END SUBROUTINE Ctqmc_allocateAll
687 : !!***
688 :
689 : !#ifdef CTCtqmc_ANALYSIS
690 : !!****f* ABINIT/m_Ctqmc/Ctqmc_allocateOpt
691 : !! NAME
692 : !! Ctqmc_allocateOpt
693 : !!
694 : !! FUNCTION
695 : !! allocate all option variables
696 : !!
697 : !! COPYRIGHT
698 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
699 : !! This file is distributed under the terms of the
700 : !! GNU General Public License, see ~abinit/COPYING
701 : !! or http://www.gnu.org/copyleft/gpl.txt .
702 : !!
703 : !! INPUTS
704 : !! this=ctqmc
705 : !!
706 : !! OUTPUT
707 : !!
708 : !! SIDE EFFECTS
709 : !!
710 : !! NOTES
711 : !!
712 : !! SOURCE
713 :
714 49 : SUBROUTINE Ctqmc_allocateOpt(this)
715 :
716 : !Arguments ------------------------------------
717 : TYPE(Ctqmc), INTENT(INOUT) :: this
718 : !Local variables ------------------------------
719 : INTEGER :: i
720 : INTEGER :: j
721 : INTEGER :: k
722 :
723 49 : IF ( .NOT. this%para ) &
724 0 : CALL ERROR("Ctqmc_allocateOpt : Ctqmc_setParameters never called ")
725 :
726 49 : IF ( this%opt_analysis .EQ. 1 ) THEN
727 6 : FREEIF(this%measCorrelation)
728 24 : MALLOC(this%measCorrelation,(1:this%samples+1,1:3,1:this%flavors))
729 13194 : this%measCorrelation = 0.d0
730 : END IF
731 :
732 49 : IF ( this%opt_order .GT. 0 ) THEN
733 42 : FREEIF(this%measPerturbation)
734 168 : MALLOC(this%measPerturbation,(1:this%opt_order,1:this%flavors))
735 24318 : this%measPerturbation = 0.d0
736 : END IF
737 :
738 49 : IF ( this%opt_histo .GT. 0 ) THEN
739 3 : FREEIF(this%occup_histo_time)
740 9 : MALLOC(this%occup_histo_time,(1:this%flavors+1))
741 48 : this%occup_histo_time= 0.d0
742 3 : FREEIF(this%occupconfig)
743 9 : MALLOC(this%occupconfig,(1:2**this%flavors))
744 49155 : this%occupconfig= 0.d0
745 3 : FREEIF(this%suscep)
746 9 : MALLOC(this%suscep,(1:3,1:this%samples))
747 603 : this%suscep= 0.d0
748 3 : FREEIF(this%chi)
749 9 : MALLOC(this%chi,(1:3,1:this%samples))
750 603 : this%chi= 0.d0
751 3 : FREEIF(this%chicharge)
752 9 : MALLOC(this%chicharge,(1:3,1:this%samples))
753 603 : this%chicharge= 0.d0
754 3 : FREEIF(this%ntot)
755 3 : MALLOC(this%ntot,(1:3))
756 12 : this%ntot= 0.d0
757 : END IF
758 :
759 49 : IF ( this%opt_noise .EQ. 1 ) THEN
760 0 : IF ( ALLOCATED(this%measNoiseG) ) THEN
761 0 : DO i=1,2
762 0 : DO j = 1, this%flavors
763 0 : DO k= 1, this%samples+1
764 0 : CALL Vector_destroy(this%measNoiseG(k,j,i))
765 : END DO
766 : END DO
767 : END DO
768 0 : DT_FREE(this%measNoiseG)
769 : END IF
770 0 : DT_MALLOC(this%measNoiseG,(1:this%samples+1,1:this%flavors,1:2))
771 : !DO i=1,2
772 0 : DO j = 1, this%flavors
773 0 : DO k= 1, this%samples+1
774 0 : CALL Vector_init(this%measNoiseG(k,j,1),CTQMC_SLICE1)
775 : END DO
776 : END DO
777 0 : DO j = 1, this%flavors
778 0 : DO k= 1, this%samples+1
779 0 : CALL Vector_init(this%measNoiseG(k,j,2),CTQMC_SLICE1*CTQMC_SLICE2+1) ! +1 pour etre remplacer ceil
780 : END DO
781 : END DO
782 : !END DO
783 0 : FREEIF(this%abNoiseG)
784 0 : MALLOC(this%aBNoiseG,(1:2,1:this%samples+1,this%flavors))
785 0 : this%abNoiseG = 0.d0
786 : END IF
787 :
788 49 : IF (this%opt_spectra .GE. 1 ) THEN
789 0 : FREEIF(this%density)
790 0 : this%endDensity=-1
791 : !MALLOC(this%density,(1:this%thermalization,1:this%flavors))
792 0 : i = CEILING(DBLE(this%thermalization+this%sweeps)/DBLE(this%measurements*this%opt_spectra))
793 0 : MALLOC(this%density,(1:this%flavors+1,1:i))
794 0 : this%endDensity=i
795 0 : this%density = 0.d0
796 : END IF
797 : !#endif
798 49 : END SUBROUTINE Ctqmc_allocateOpt
799 : !!***
800 :
801 0 : SUBROUTINE Ctqmc_setG0wFile(this,istream,opt_fk)
802 :
803 :
804 : #ifdef HAVE_MPI1
805 : include 'mpif.h'
806 : #endif
807 : TYPE(Ctqmc) , INTENT(INOUT) :: this
808 : INTEGER , INTENT(IN ) :: istream
809 : INTEGER , INTENT(IN ) :: opt_fk
810 : COMPLEX(KIND=8), DIMENSION(:,:), ALLOCATABLE :: Gomega
811 : INTEGER :: flavors
812 : INTEGER :: iflavor
813 : INTEGER :: iomega
814 : #ifdef HAVE_MPI
815 : INTEGER :: ierr
816 : #endif
817 :
818 0 : IF ( .NOT. this%para ) &
819 0 : CALL ERROR("Ctqmc_setG0wFile : Ctqmc_setParameters never called ")
820 :
821 0 : flavors = this%flavors
822 :
823 0 : this%Wmax = this%samples
824 0 : MALLOC(Gomega,(1:this%Wmax,1:flavors))
825 :
826 0 : IF ( this%rank .EQ. 0 ) THEN
827 0 : DO iomega=1, this%Wmax
828 0 : READ(istream,*) (Gomega(iomega,iflavor),iflavor=1,flavors)
829 : END DO
830 : END IF
831 :
832 : #ifdef HAVE_MPI
833 : CALL MPI_Bcast(Gomega, this%Wmax*flavors, MPI_DOUBLE_COMPLEX, 0, &
834 0 : this%MY_COMM, ierr)
835 : #endif
836 :
837 0 : CALL Ctqmc_setG0wTab(this,Gomega,opt_fk)
838 0 : FREE(Gomega)
839 :
840 0 : END SUBROUTINE Ctqmc_setG0wFile
841 : !!***
842 :
843 : !!****f* ABINIT/m_Ctqmc/Ctqmc_setG0wTab
844 : !! NAME
845 : !! Ctqmc_setG0wTab
846 : !!
847 : !! FUNCTION
848 : !! Set Gow from input array
849 : !!
850 : !! COPYRIGHT
851 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
852 : !! This file is distributed under the terms of the
853 : !! GNU General Public License, see ~abinit/COPYING
854 : !! or http://www.gnu.org/copyleft/gpl.txt .
855 : !!
856 : !! INPUTS
857 : !! this=ctqmc
858 : !! Gomega=G0w
859 : !! opt_fk=F is already inversed with out iwn
860 : !!
861 : !! OUTPUT
862 : !!
863 : !! SIDE EFFECTS
864 : !!
865 : !! NOTES
866 : !!
867 : !! SOURCE
868 :
869 49 : SUBROUTINE Ctqmc_setG0wTab(this,Gomega,opt_fk)
870 :
871 : !Arguments ------------------------------------
872 : TYPE(Ctqmc), INTENT(INOUT) :: this
873 : COMPLEX(KIND=8), DIMENSION(:,:), INTENT(IN ) :: Gomega
874 : INTEGER , INTENT(IN ) :: opt_fk
875 : !Local variable -------------------------------
876 : DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: F
877 :
878 49 : IF ( .NOT. this%para ) &
879 0 : CALL ERROR("Ctqmc_setG0wTab : Ctqmc_setParameters never called ")
880 :
881 196 : MALLOC(F,(1:this%samples+1,1:this%flavors))
882 49 : CALL Ctqmc_computeF(this,Gomega, F, opt_fk) ! mu is changed
883 49 : CALL BathOperator_setF(this%Bath, F)
884 : !CALL BathOperator_printF(this%Bath)
885 49 : FREE(F)
886 :
887 49 : IF ( this%opt_levels .NE. 1 ) THEN ! We compute the mu by hand in computeF
888 0 : CALL ImpurityOperator_setMu(this%Impurity,this%mu)
889 : END IF
890 :
891 49 : this%inF = .TRUE.
892 49 : this%set = .TRUE.
893 :
894 49 : END SUBROUTINE Ctqmc_setG0wTab
895 : !!***
896 :
897 : !SUBROUTINE Ctqmc_setFwK(this,Gomega)
898 : ! COMPLEX*16 , DIMENSION(:,:), INTENT(IN ) :: Gomega
899 : ! DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: F
900 : !
901 : ! IF ( .NOT. this%para ) &
902 : ! CALL ERROR("Ctqmc_setG0wTab : Ctqmc_setParameters never called ")
903 : !
904 : ! MALLOC(F,(1:this%samples+1,1:this%flavors))
905 : ! CALL Ctqmc_computeFK(this,Gomega, this%Wmax, F) ! mu is changed
906 : ! CALL BathOperator_setF(this%Bath, F)
907 : ! CALL BathOperator_printF(this%Bath)
908 : ! FREE(F)
909 : !
910 : ! this%inF = .TRUE.
911 : ! this%set = .TRUE.
912 : !
913 : !END SUBROUTINE Ctqmc_setFwK
914 : !!***
915 :
916 : !SUBROUTINE Ctqmc_setBand(this, mu, U)
917 : ! DOUBLE PRECISION, INTENT(IN ) :: mu
918 : ! DOUBLE PRECISION, INTENT(IN ) :: U
919 : !
920 : ! IF ( .NOT. this%para ) &
921 : ! CALL ERROR("Ctqmc_setBand : Ctqmc_setParameters never called ")
922 : !
923 : ! CALL ImpurityOperator_setU(this%Impurity, U, 0.d0)
924 : ! !this%mu = mu + this%Impurity%shift_mu
925 : !END SUBROUTINE Ctqmc_setBand
926 : !!***
927 :
928 : !!****f* ABINIT/m_Ctqmc/Ctqmc_setU
929 : !! NAME
930 : !! Ctqmc_setU
931 : !!
932 : !! FUNCTION
933 : !! set the interaction this
934 : !!
935 : !! COPYRIGHT
936 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
937 : !! This file is distributed under the terms of the
938 : !! GNU General Public License, see ~abinit/COPYING
939 : !! or http://www.gnu.org/copyleft/gpl.txt .
940 : !!
941 : !! INPUTS
942 : !! this=ctqmc
943 : !! matU=interaction this
944 : !!
945 : !! OUTPUT
946 : !!
947 : !! SIDE EFFECTS
948 : !!
949 : !! NOTES
950 : !!
951 : !! SOURCE
952 :
953 49 : SUBROUTINE Ctqmc_setU(this,matU)
954 :
955 : !Arguments ------------------------------------
956 : TYPE(Ctqmc), INTENT(INOUT) :: this
957 : !Local variables ------------------------------
958 : DOUBLE PRECISION, DIMENSION(:,:), INTENT(IN) :: matU
959 :
960 147 : IF ( SIZE(matU) .NE. this%flavors*this%flavors ) &
961 0 : CALL ERROR("Ctqmc_setU : Wrong interaction this (size) ")
962 :
963 49 : CALL ImpurityOperator_setUmat(this%Impurity, matU)
964 49 : this%setU = .TRUE.
965 49 : END SUBROUTINE Ctqmc_setU
966 : !!***
967 :
968 : !!****f* ABINIT/m_Ctqmc/Ctqmc_clear
969 : !! NAME
970 : !! Ctqmc_clear
971 : !!
972 : !! FUNCTION
973 : !! clear a ctqmc run
974 : !!
975 : !! COPYRIGHT
976 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
977 : !! This file is distributed under the terms of the
978 : !! GNU General Public License, see ~abinit/COPYING
979 : !! or http://www.gnu.org/copyleft/gpl.txt .
980 : !!
981 : !! INPUTS
982 : !! this=ctqmc
983 : !!
984 : !! OUTPUT
985 : !!
986 : !! SIDE EFFECTS
987 : !!
988 : !! NOTES
989 : !!
990 : !! SOURCE
991 :
992 98 : SUBROUTINE Ctqmc_clear(this)
993 :
994 : !Arguments ------------------------------------
995 : TYPE(Ctqmc), INTENT(INOUT) :: this
996 : !Local variables ------------------------------
997 : INTEGER :: i
998 : INTEGER :: j
999 : INTEGER :: k
1000 :
1001 1134 : this%measN(1,:) = 0.d0
1002 1134 : this%measN(2,:) = 0.d0
1003 : !Do not set measN(3,:) to 0 to avoid erasing N between therm and ctqmc
1004 1134 : this%measN(4,:) = 0.d0
1005 12726 : this%measDE = 0.d0
1006 : ! this%seg_added = 0.d0
1007 : ! this%anti_added = 0.d0
1008 : ! this%seg_removed = 0.d0
1009 : ! this%anti_removed = 0.d0
1010 : ! this%seg_sign = 0.d0
1011 : ! this%anti_sign = 0.d0
1012 686 : this%stats(:) = 0.d0
1013 98 : this%swap = 0.d0
1014 98 : this%runTime = 0.d0
1015 98 : this%modGlobalMove(2) = 0
1016 98 : CALL Vector_clear(this%measNoise(1))
1017 98 : CALL Vector_clear(this%measNoise(2))
1018 : !#ifdef CTCtqmc_CHECK
1019 98 : this%errorImpurity = 0.d0
1020 98 : this%errorBath = 0.d0
1021 : !#endif
1022 1134 : DO j = 1, this%flavors
1023 1134 : CALL GreenHyb_clear(this%Greens(j))
1024 : END DO
1025 : !#ifdef CTCtqmc_ANALYSIS
1026 98 : IF ( this%opt_analysis .EQ. 1 .AND. ALLOCATED(this%measCorrelation) ) &
1027 13194 : this%measCorrelation = 0.d0
1028 98 : IF ( this%opt_order .GT. 0 .AND. ALLOCATED(this%measPerturbation) ) &
1029 24318 : this%measPerturbation = 0.d0
1030 98 : IF ( this%opt_noise .EQ. 1 .AND. ALLOCATED(this%measNoiseG) ) THEN
1031 0 : DO i=1,2
1032 0 : DO j = 1, this%flavors
1033 0 : DO k= 1, this%samples+1
1034 0 : CALL Vector_clear(this%measNoiseG(k,j,i))
1035 : END DO
1036 : END DO
1037 : END DO
1038 : END IF
1039 : !#endif
1040 98 : END SUBROUTINE Ctqmc_clear
1041 : !!***
1042 :
1043 : !!****f* ABINIT/m_Ctqmc/Ctqmc_reset
1044 : !! NAME
1045 : !! Ctqmc_reset
1046 : !!
1047 : !! FUNCTION
1048 : !! reset a ctqmc simulation
1049 : !!
1050 : !! COPYRIGHT
1051 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
1052 : !! This file is distributed under the terms of the
1053 : !! GNU General Public License, see ~abinit/COPYING
1054 : !! or http://www.gnu.org/copyleft/gpl.txt .
1055 : !!
1056 : !! INPUTS
1057 : !! this=ctqmc
1058 : !!
1059 : !! OUTPUT
1060 : !!
1061 : !! SIDE EFFECTS
1062 : !!
1063 : !! NOTES
1064 : !!
1065 : !! SOURCE
1066 :
1067 49 : SUBROUTINE Ctqmc_reset(this)
1068 :
1069 : !Arguments ------------------------------------
1070 : TYPE(Ctqmc), INTENT(INOUT) :: this
1071 : !Local variables ------------------------------
1072 : INTEGER :: iflavor
1073 : DOUBLE PRECISION :: sweeps
1074 :
1075 567 : DO iflavor = 1, this%flavors
1076 567 : CALL GreenHyb_reset(this%Greens(iflavor))
1077 : END DO
1078 49 : CALL Ctqmc_clear(this)
1079 49 : CALL ImpurityOperator_reset(this%Impurity)
1080 49 : CALL BathOperator_reset (this%Bath)
1081 567 : this%measN(3,:) = 0.d0
1082 : !complete restart -> measN=0
1083 49 : this%done = .FALSE.
1084 49 : this%set = .FALSE.
1085 49 : this%inF = .FALSE.
1086 49 : this%opt_movie = 0
1087 49 : this%opt_analysis = 0
1088 49 : this%opt_order = 0
1089 49 : this%opt_check = 0
1090 49 : this%opt_noise = 0
1091 49 : this%opt_spectra = 0
1092 49 : this%opt_levels = 0
1093 49 : sweeps = DBLE(this%sweeps)*DBLE(this%size)
1094 49 : CALL Ctqmc_setSweeps(this, sweeps)
1095 : !#ifdef HAVE_MPI
1096 : ! CALL MPI_BARRIER(this%MY_COMM,iflavor)
1097 : ! IF ( this%rank .EQ. 0 ) &
1098 : !#endif
1099 : ! WRITE(this%ostream,'(A9)') "QMC reset"
1100 : ! CALL FLUSH(this%ostream)
1101 49 : END SUBROUTINE Ctqmc_reset
1102 : !!***
1103 :
1104 : !!****f* ABINIT/m_Ctqmc/Ctqmc_setMu
1105 : !! NAME
1106 : !! Ctqmc_setMu
1107 : !!
1108 : !! FUNCTION
1109 : !! impose energy levels
1110 : !!
1111 : !! COPYRIGHT
1112 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
1113 : !! This file is distributed under the terms of the
1114 : !! GNU General Public License, see ~abinit/COPYING
1115 : !! or http://www.gnu.org/copyleft/gpl.txt .
1116 : !!
1117 : !! INPUTS
1118 : !! this=ctqmc
1119 : !! levels=energy levels vector
1120 : !!
1121 : !! OUTPUT
1122 : !! argout(sizeout)=description
1123 : !!
1124 : !! SIDE EFFECTS
1125 : !!
1126 : !! NOTES
1127 : !!
1128 : !! SOURCE
1129 :
1130 49 : SUBROUTINE Ctqmc_setMu(this, levels)
1131 :
1132 : !Arguments ------------------------------------
1133 : TYPE(Ctqmc) , INTENT(INOUT) :: this
1134 : DOUBLE PRECISION, DIMENSION(:), INTENT(IN ) :: levels
1135 :
1136 49 : IF ( this%flavors .NE. SIZE(levels,1) ) &
1137 0 : CALL WARNALL("Ctqmc_setMu : Taking energy levels from weiss G(iw)")
1138 :
1139 567 : this%mu(:)=-levels ! levels = \epsilon_j - \mu
1140 : !this%mu =\tilde{\mu} = \mu -\epsilon_j
1141 49 : CALL ImpurityOperator_setMu(this%Impurity,this%mu)
1142 49 : this%opt_levels = 1
1143 49 : END SUBROUTINE Ctqmc_setMu
1144 : !!***
1145 :
1146 : !!****f* ABINIT/m_Ctqmc/Ctqmc_computeF
1147 : !! NAME
1148 : !! Ctqmc_computeF
1149 : !!
1150 : !! FUNCTION
1151 : !! Compute the hybridization function
1152 : !!
1153 : !! COPYRIGHT
1154 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
1155 : !! This file is distributed under the terms of the
1156 : !! GNU General Public License, see ~abinit/COPYING
1157 : !! or http://www.gnu.org/copyleft/gpl.txt .
1158 : !!
1159 : !! INPUTS
1160 : !! this=ctqmc
1161 : !! Gomega=G0 to compute F
1162 : !! opt_fk=What is Gomega
1163 : !!
1164 : !! OUTPUT
1165 : !! F=hybridization function
1166 : !!
1167 : !! SIDE EFFECTS
1168 : !!
1169 : !! NOTES
1170 : !!
1171 : !! SOURCE
1172 :
1173 49 : SUBROUTINE Ctqmc_computeF(this, Gomega, F, opt_fk)
1174 :
1175 : !Arguments ------------------------------------
1176 : TYPE(Ctqmc) , INTENT(INOUT) :: this
1177 : COMPLEX(KIND=8), DIMENSION(:,:), INTENT(IN ) :: Gomega
1178 : !INTEGER , INTENT(IN ) :: Wmax
1179 : DOUBLE PRECISION, DIMENSION(:,:), INTENT(INOUT) :: F
1180 : INTEGER , INTENT(IN ) :: opt_fk
1181 : !Local variables ------------------------------
1182 : INTEGER :: flavors
1183 : INTEGER :: samples
1184 : INTEGER :: iflavor
1185 : INTEGER :: iomega
1186 : INTEGER :: itau
1187 : DOUBLE PRECISION :: pi_invBeta
1188 : DOUBLE PRECISION :: K
1189 : DOUBLE PRECISION :: re
1190 : DOUBLE PRECISION :: im
1191 49 : COMPLEX(KIND=8), DIMENSION(:,:), ALLOCATABLE :: F_omega
1192 49 : TYPE(GreenHyb) :: F_tmp
1193 :
1194 49 : flavors = this%flavors
1195 :
1196 49 : samples = this%samples
1197 49 : pi_invBeta = ACOS(-1.d0) / this%beta
1198 49 : this%Wmax=SIZE(Gomega,1)
1199 :
1200 49 : IF ( this%have_MPI .EQV. .TRUE. ) THEN
1201 49 : CALL GreenHyb_init(F_tmp,samples,this%beta,MY_COMM=this%MY_COMM)
1202 : ELSE
1203 0 : CALL GreenHyb_init(F_tmp,samples,this%beta)
1204 : END IF
1205 : ! K = this%mu
1206 :
1207 196 : MALLOC(F_omega,(1:this%Wmax,1:flavors))
1208 :
1209 : !IF ( this%rank .EQ. 0 ) &
1210 : !OPEN(UNIT=9876,FILE="K.dat",POSITION="APPEND")
1211 49 : IF ( opt_fk .EQ. 0 ) THEN
1212 0 : DO iflavor = 1, flavors
1213 0 : DO iomega=1,this%Wmax
1214 0 : re = REAL(Gomega(iomega,iflavor))
1215 0 : im = AIMAG(Gomega(iomega,iflavor))
1216 0 : F_omega(iomega,iflavor) = CMPLX(-re/(re*re+im*im),im/(re*re+im*im),8)
1217 : END DO
1218 : END DO
1219 : !F_omega = CMPLX(-1.d0,0,8)/Gomega
1220 : ELSE
1221 152516 : F_omega = Gomega
1222 : END IF
1223 :
1224 567 : DO iflavor = 1, flavors
1225 518 : IF ( this%opt_levels .EQ. 1 ) THEN
1226 518 : K = this%mu(iflavor)
1227 : ELSE
1228 0 : K = -REAL(F_omega(this%Wmax, iflavor))
1229 : ! this%mu = K
1230 0 : this%mu(iflavor) = K
1231 : END IF
1232 : !IF ( this%rank .EQ. 0 ) &
1233 : !WRITE(9876,'(I4,2E22.14)') iflavor, K, REAL(-F_omega(this%Wmax, iflavor))
1234 : !IF(this%rank .EQ.0) &
1235 : !WRITE(this%ostream,*) "CTQMC K, this%mu = ",K,this%mu(iflavor)
1236 : !WRITE(this%ostream,*) "CTQMC beta = ",this%beta
1237 518 : IF ( opt_fk .EQ. 0 ) THEN
1238 0 : DO iomega = 1, this%Wmax
1239 0 : re = REAL(F_omega(iomega,iflavor))
1240 0 : im = AIMAG(F_omega(iomega,iflavor))
1241 0 : F_omega(iomega,iflavor) = CMPLX(re + K, im + (2.d0*DBLE(iomega)-1.d0) * pi_invBeta, 8)
1242 : !if(iflavor==1.and.this%rank==0) then
1243 : !write(224,*) (2.d0*DBLE(iomega)-1.d0) * pi_invBeta, real(F_omega(iomega,iflavor)),imag(F_omega(iomega,iflavor))
1244 : !write(225,*) (2.d0*DBLE(iomega)-1.d0) * pi_invBeta, real(Gomega(iomega, iflavor)),imag(Gomega(iomega, iflavor))
1245 : !end if
1246 : END DO
1247 : ELSE
1248 152418 : DO iomega = 1, this%Wmax
1249 : F_omega(iomega,iflavor) = F_omega(iomega,iflavor) &
1250 152418 : + CMPLX(K, 0.d0, 8)
1251 : !if(iflavor==1.and.this%rank==0) then
1252 : !write(224,*) (2.d0*DBLE(iomega)-1.d0) * pi_invBeta, real(F_omega(iomega,iflavor)),imag(F_omega(iomega,iflavor))
1253 : !write(225,*) (2.d0*DBLE(iomega)-1.d0) * pi_invBeta, real(Gomega(iomega, iflavor)),imag(Gomega(iomega, iflavor))
1254 : !end if
1255 : END DO
1256 : END IF
1257 518 : K = REAL(CMPLX(0,(2.d0*DBLE(this%Wmax)-1.d0)*pi_invBeta,8)*F_omega(this%Wmax,iflavor))
1258 518 : CALL GreenHyb_setMuD1(this%Greens(iflavor),this%mu(iflavor),K)
1259 518 : CALL GreenHyb_setOperW(F_tmp,F_omega(:,iflavor))
1260 : !CALL GreenHyb_backFourier(F_tmp,F_omega(:,iflavor))
1261 518 : CALL GreenHyb_backFourier(F_tmp)
1262 305403 : F(1:samples+1,iflavor) = (/ (-F_tmp%oper(samples+1-itau),itau=0,samples) /)
1263 : END DO
1264 49 : IF ( this%rank .EQ. 0 ) THEN
1265 25 : open (unit=346,file='Hybridization.dat',status='unknown',form='formatted')
1266 279 : DO iflavor = 1, flavors
1267 254 : write(346,*) "#",iflavor
1268 58208 : do itau=1,this%samples+1
1269 58208 : write(346,*) itau,F(itau,iflavor)
1270 : enddo
1271 279 : write(346,*)
1272 : END DO
1273 25 : close(346)
1274 : ENDIF
1275 49 : FREE(F_omega)
1276 49 : CALL GreenHyb_destroy(F_tmp)
1277 49 : END SUBROUTINE Ctqmc_computeF
1278 : !!***
1279 :
1280 : !SUBROUTINE Ctqmc_computeFK(this, Gomega, Wmax, F)
1281 : ! COMPLEX*16 , DIMENSION(:,:), INTENT(IN ) :: Gomega
1282 : ! INTEGER , INTENT(IN ) :: Wmax
1283 : ! DOUBLE PRECISION, DIMENSION(:,:), INTENT(INOUT) :: F
1284 : ! INTEGER :: flavors
1285 : ! INTEGER :: samples
1286 : ! INTEGER :: iflavor
1287 : ! INTEGER :: iomega
1288 : ! INTEGER :: itau
1289 : ! DOUBLE PRECISION :: pi_invBeta
1290 : ! DOUBLE PRECISION :: K
1291 : ! COMPLEX*16 , DIMENSION(:,:), ALLOCATABLE :: F_omega
1292 : ! TYPE(GreenHyb) :: F_tmp
1293 : !
1294 : ! flavors = this%flavors
1295 : !
1296 : ! samples = this%samples
1297 : ! pi_invBeta = ACOS(-1.d0) / this%beta
1298 : !
1299 : ! IF ( this%have_MPI .EQV. .TRUE. ) THEN
1300 : ! CALL GreenHyb_init(F_tmp,samples,this%beta,this%MY_COMM)
1301 : ! ELSE
1302 : ! CALL GreenHyb_init(F_tmp,samples,this%beta)
1303 : ! END IF
1304 : !! K = this%mu
1305 : !
1306 : ! MALLOC(F_omega,(1:Wmax,1:flavors))
1307 : !
1308 : ! DO iflavor = 1, flavors
1309 : ! K = REAL(Gomega(Wmax, iflavor))
1310 : ! WRITE(this%ostream,*) "CTQMC K, this%mu = ",K,this%mu
1311 : ! WRITE(this%ostream,*) "CTQMC beta = ",this%beta
1312 : ! this%mu(iflavor) = K
1313 : ! DO iomega = 1, Wmax
1314 : ! F_omega(iomega,iflavor) = Gomega(iomega,iflavor) &
1315 : ! - CMPLX(K, 0.d0, 8)
1316 : ! !if(iflavor==1.and.this%rank==0) then
1317 : ! !write(224,*) (2.d0*DBLE(iomega)-1.d0) * pi_invBeta, real(F_omega(iomega,iflavor)),imag(F_omega(iomega,iflavor))
1318 : ! !write(225,*) (2.d0*DBLE(iomega)-1.d0) * pi_invBeta, real(Gomega(iomega, iflavor)),imag(Gomega(iomega, iflavor))
1319 : ! !end if
1320 : ! END DO
1321 : ! CALL GreenHyb_backFourier(F_tmp,F_omega(:,iflavor))
1322 : ! F(1:samples+1,iflavor) = (/ (-F_tmp%oper(samples+1-itau),itau=0,samples) /)
1323 : ! END DO
1324 : ! FREE(F_omega)
1325 : ! CALL GreenHyb_destroy(F_tmp)
1326 : !END SUBROUTINE Ctqmc_computeFK
1327 : !!***
1328 :
1329 :
1330 : !!****f* ABINIT/m_Ctqmc/Ctqmc_run
1331 : !! NAME
1332 : !! Ctqmc_run
1333 : !!
1334 : !! FUNCTION
1335 : !! set all options and run a simulation
1336 : !!
1337 : !! COPYRIGHT
1338 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
1339 : !! This file is distributed under the terms of the
1340 : !! GNU General Public License, see ~abinit/COPYING
1341 : !! or http://www.gnu.org/copyleft/gpl.txt .
1342 : !!
1343 : !! INPUTS
1344 : !! this=ctqmc
1345 : !! opt_order=maximal perturbation order to scope
1346 : !! opt_movie=draw a movie of the simulation
1347 : !! opt_analysis=compute correlation functions
1348 : !! opt_check=check fast calculations
1349 : !! opt_noise=compute noise for green function
1350 : !! opt_spectra=fourier transform of the time evolution of the number of electrons
1351 : !! opt_gMove=steps without global move
1352 : !!
1353 : !! OUTPUT
1354 : !!
1355 : !! SIDE EFFECTS
1356 : !!
1357 : !! NOTES
1358 : !!
1359 : !! SOURCE
1360 :
1361 49 : SUBROUTINE Ctqmc_run(this,opt_order,opt_histo,opt_movie,opt_analysis,opt_check,opt_noise,opt_spectra,opt_gMove)
1362 :
1363 :
1364 : #ifdef HAVE_MPI1
1365 : include 'mpif.h'
1366 : #endif
1367 : !Arguments ------------------------------------
1368 : TYPE(Ctqmc), INTENT(INOUT) :: this
1369 : INTEGER, OPTIONAL, INTENT(IN ) :: opt_order
1370 : INTEGER, OPTIONAL, INTENT(IN ) :: opt_histo
1371 : INTEGER, OPTIONAL, INTENT(IN ) :: opt_movie
1372 : INTEGER, OPTIONAL, INTENT(IN ) :: opt_analysis
1373 : INTEGER, OPTIONAL, INTENT(IN ) :: opt_check
1374 : INTEGER, OPTIONAL, INTENT(IN ) :: opt_noise
1375 : INTEGER, OPTIONAL, INTENT(IN ) :: opt_spectra
1376 : INTEGER, OPTIONAL, INTENT(IN ) :: opt_gMove
1377 : !Local variables ------------------------------
1378 : #ifdef HAVE_MPI
1379 : INTEGER :: ierr
1380 : #endif
1381 : !#ifdef CTCtqmc_MOVIE
1382 : INTEGER :: ilatex
1383 : CHARACTER(LEN=4) :: Cchar
1384 : !#endif
1385 : DOUBLE PRECISION :: estimatedTime(1)
1386 :
1387 49 : IF ( .NOT. this%set ) &
1388 0 : CALL ERROR("Ctqmc_run : QMC not set up ")
1389 49 : IF ( .NOT. this%setU ) &
1390 0 : CALL ERROR("Ctqmc_run : QMC does not have a U this ")
1391 :
1392 : #ifdef DEBUG_VERBOSE
1393 : write(this%ostream,'(a,i7,a,i3,a,i8,a,i15,a)') &
1394 : "Walker with rank ", this%rank, " and tid ", this%tid, " runs with seed ", this%seed, " for ", this%sweeps, " sweeps"
1395 : #endif
1396 :
1397 : ! OPTIONS of the run
1398 49 : IF ( PRESENT( opt_check ) ) THEN
1399 49 : this%opt_check = opt_check
1400 49 : CALL ImpurityOperator_doCheck(this%Impurity,opt_check)
1401 49 : CALL BathOperator_doCheck(this%Bath,opt_check)
1402 : END IF
1403 49 : IF ( PRESENT( opt_movie ) ) &
1404 49 : this%opt_movie = opt_movie
1405 49 : IF ( PRESENT( opt_analysis ) ) &
1406 49 : this%opt_analysis = opt_analysis
1407 49 : IF ( PRESENT ( opt_order ) ) &
1408 49 : this%opt_order = opt_order
1409 49 : IF ( PRESENT ( opt_histo ) ) &
1410 49 : this%opt_histo = opt_histo
1411 49 : IF ( PRESENT ( opt_noise ) ) THEN
1412 49 : this%opt_noise = opt_noise
1413 : END IF
1414 49 : IF ( PRESENT ( opt_spectra ) ) &
1415 49 : this%opt_spectra = opt_spectra
1416 :
1417 49 : this%modGlobalMove(1) = this%sweeps+1 ! No Global Move
1418 49 : this%modGlobalMove(2) = 0
1419 49 : IF ( PRESENT ( opt_gMove ) ) THEN
1420 49 : IF ( opt_gMove .LE. 0 .OR. opt_gMove .GT. this%sweeps ) THEN
1421 : this%modGlobalMove(1) = this%sweeps+1
1422 28 : CALL WARNALL("Ctqmc_run : global moves option is <= 0 or > sweeps/cpu -> No global Moves")
1423 : ELSE
1424 21 : this%modGlobalMove(1) = opt_gMove
1425 : END IF
1426 : END IF
1427 :
1428 49 : CALL Ctqmc_allocateOpt(this)
1429 :
1430 : !#ifdef CTCtqmc_MOVIE
1431 49 : ilatex = 0
1432 49 : IF ( this%opt_movie .EQ. 1 ) THEN
1433 0 : Cchar ="0000"
1434 0 : WRITE(Cchar,'(I4)') this%rank + (this%tid-1) * this%size
1435 0 : ilatex = 87+this%rank
1436 0 : OPEN(UNIT=ilatex, FILE="Movie_"//TRIM(ADJUSTL(Cchar))//".tex")
1437 0 : WRITE(ilatex,'(A)') "\documentclass{beamer}"
1438 0 : WRITE(ilatex,'(A)') "\usepackage{color}"
1439 0 : WRITE(ilatex,'(A)') "\setbeamersize{sidebar width left=0pt}"
1440 0 : WRITE(ilatex,'(A)') "\setbeamersize{sidebar width right=0pt}"
1441 : !FIXME Those lines don't work anymore, but output is still readable...
1442 : !WRITE(ilatex,'(A)') "\setbeamersize{text width left=0pt}"
1443 : !WRITE(ilatex,'(A)') "\setbeamersize{text width right=0pt}"
1444 0 : WRITE(ilatex,*)
1445 0 : WRITE(ilatex,'(A)') "\begin{document}"
1446 0 : WRITE(ilatex,*)
1447 : END IF
1448 : !#endif
1449 :
1450 49 : IF ( this%rank .EQ. 0 .AND. this%tid .EQ. 1 ) THEN
1451 19 : WRITE(this%ostream,'(A29)') "Starting QMC (Thermalization)"
1452 : END IF
1453 :
1454 : !=================================
1455 : ! STARTING THERMALIZATION
1456 : !=================================
1457 49 : CALL Ctqmc_loop(this,this%thermalization,ilatex)
1458 : !=================================
1459 : ! ENDING THERMALIZATION
1460 : !=================================
1461 :
1462 98 : estimatedTime = this%runTime
1463 : #ifdef HAVE_MPI
1464 49 : IF ( this%tid .EQ. 1 ) THEN
1465 : CALL MPI_REDUCE([this%runTime], estimatedTime, 1, MPI_DOUBLE_PRECISION, MPI_MAX, &
1466 68 : 0, this%MY_COMM, ierr)
1467 : END IF
1468 : #endif
1469 :
1470 49 : IF ( this%rank .EQ. 0 .AND. this%tid .EQ. 1 ) THEN
1471 38 : WRITE(this%ostream,'(A26,I6,A11)') "Thermalization done in ", CEILING(estimatedTime), " seconds"
1472 19 : WRITE(this%ostream,'(A25,I7,A15,I5,A5)') "The QMC should run in ", &
1473 19 : CEILING(estimatedTime(1)*DBLE(this%sweeps)/DBLE(this%thermalization)),&
1474 38 : " seconds on ", this%size, " CPUs"
1475 : END IF
1476 :
1477 : !=================================
1478 : ! CLEANING CTQMC
1479 : !=================================
1480 49 : CALL Ctqmc_clear(this)
1481 :
1482 : !=================================
1483 : ! STARTING CTQMC
1484 : !=================================
1485 49 : CALL Ctqmc_loop(this,this%sweeps,ilatex)
1486 : !=================================
1487 : ! ENDING CTQMC
1488 : !=================================
1489 :
1490 49 : IF ( this%opt_movie .EQ. 1 ) THEN
1491 0 : WRITE(ilatex,*) ""
1492 0 : WRITE(ilatex,'(A14)') "\end{document}"
1493 0 : CLOSE(ilatex)
1494 : END IF
1495 :
1496 49 : this%done = .TRUE.
1497 :
1498 49 : END SUBROUTINE Ctqmc_run
1499 : !!***
1500 :
1501 : !!****f* ABINIT/m_Ctqmc/Ctqmc_loop
1502 : !! NAME
1503 : !! Ctqmc_loop
1504 : !!
1505 : !! FUNCTION
1506 : !! Definition the main loop of the CT-QMC
1507 : !!
1508 : !! COPYRIGHT
1509 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
1510 : !! This file is distributed under the terms of the
1511 : !! GNU General Public License, see ~abinit/COPYING
1512 : !! or http://www.gnu.org/copyleft/gpl.txt .
1513 : !!
1514 : !! INPUTS
1515 : !! this=ctqmc
1516 : !! itotal=number of sweeps to perform : thermalization or sweeps
1517 : !! ilatex=unit of file to write movie if so
1518 : !!
1519 : !! OUTPUT
1520 : !!
1521 : !! SIDE EFFECTS
1522 : !!
1523 : !! NOTES
1524 : !!
1525 : !! SOURCE
1526 :
1527 98 : SUBROUTINE Ctqmc_loop(this,itotal,ilatex)
1528 :
1529 : !Arguments ------------------------------------
1530 : TYPE(Ctqmc), INTENT(INOUT) :: this
1531 : INTEGER , INTENT(IN ) :: itotal
1532 : INTEGER , INTENT(IN ) :: ilatex
1533 : !Local variables ------------------------------
1534 : LOGICAL :: updated
1535 : LOGICAL :: updated_seg
1536 98 : LOGICAL, DIMENSION(:), ALLOCATABLE :: updated_swap
1537 :
1538 : INTEGER :: flavors
1539 : INTEGER :: measurements
1540 : INTEGER :: modNoise1
1541 : INTEGER :: modNoise2
1542 : INTEGER :: modGlobalMove
1543 : INTEGER :: sp1
1544 : INTEGER :: itau
1545 : INTEGER :: ind
1546 : INTEGER :: endDensity
1547 : INTEGER :: indDensity
1548 : INTEGER :: swapUpdate1
1549 : INTEGER :: swapUpdate2
1550 : INTEGER :: old_percent
1551 : INTEGER :: new_percent
1552 : INTEGER :: ipercent
1553 : INTEGER :: iflavor
1554 : INTEGER :: isweep
1555 :
1556 : DOUBLE PRECISION :: cpu_time1
1557 : DOUBLE PRECISION :: cpu_time2
1558 : DOUBLE PRECISION :: NRJ_old1
1559 : DOUBLE PRECISION :: NRJ_old2
1560 : DOUBLE PRECISION :: NRJ_new
1561 : DOUBLE PRECISION :: total
1562 98 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: gtmp_old1
1563 98 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: gtmp_old2
1564 98 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: gtmp_new
1565 :
1566 98 : CALL CPU_TIME(cpu_time1)
1567 :
1568 98 : flavors = this%flavors
1569 98 : measurements = this%measurements
1570 98 : modNoise1 = this%modNoise1
1571 98 : modNoise2 = this%modNoise2
1572 98 : modGlobalMove = this%modGlobalMove(1)
1573 98 : sp1 = this%samples+1
1574 98 : IF ( this%opt_histo .GT. 0 ) THEN
1575 96 : this%occup_histo_time= 0.d0
1576 98310 : this%occupconfig= 0.d0
1577 1206 : this%suscep= 0.d0
1578 1206 : this%chi= 0.d0
1579 1206 : this%chicharge= 0.d0
1580 24 : this%ntot = 0.d0
1581 : END IF
1582 98 : old_percent = 0
1583 :
1584 294 : MALLOC(updated_swap,(1:flavors))
1585 1134 : updated_swap(:) = .FALSE.
1586 :
1587 98 : NRJ_old1 = 0.d0
1588 98 : NRJ_old2 = 0.d0
1589 98 : NRJ_new = 0.d0
1590 :
1591 98 : MALLOC(gtmp_new,(1,1))
1592 392 : gtmp_new = 0.d0
1593 98 : MALLOC(gtmp_old1,(1,1))
1594 392 : gtmp_old1 = 0.d0
1595 98 : MALLOC(gtmp_old2,(1,1))
1596 392 : gtmp_old2 = 0.d0
1597 :
1598 : !PROBLEM eos_gnu_13.2_openmpi . %endDensity was introduced throughout
1599 : ! endDensity = SIZE(this%density,2)
1600 98 : endDensity=this%endDensity
1601 : !ENDPROBLEM
1602 :
1603 98 : IF ( this%opt_noise .GT. 0 ) THEN
1604 0 : FREEIF(gtmp_new)
1605 0 : MALLOC(gtmp_new,(1:sp1,1:flavors))
1606 0 : FREEIF(gtmp_old1)
1607 0 : MALLOC(gtmp_old1,(1:sp1,1:flavors))
1608 0 : FREEIF(gtmp_old2)
1609 0 : MALLOC(gtmp_old2,(1:sp1,1:flavors))
1610 : END IF
1611 :
1612 98 : IF ( this%rank .EQ. 0 .AND. this%tid .EQ. 1 ) THEN
1613 : WRITE(this%ostream, '(1x,103A)') &
1614 38 : "|----------------------------------------------------------------------------------------------------|"
1615 38 : WRITE(this%ostream,'(1x,A)', ADVANCE="NO") "|"
1616 : END IF
1617 :
1618 98 : total = DBLE(itotal)
1619 :
1620 98 : indDensity = 1
1621 78610098 : DO isweep = 1, itotal
1622 1100990000 : DO iflavor = 1, flavors
1623 1022380000 : ImpurityOperator_QuickActivation(this%Impurity,iflavor)
1624 1022380000 : BathOperator_QuickActivation(this%Bath,iflavor)
1625 :
1626 1022380000 : CALL Ctqmc_tryAddRemove(this,updated_seg)
1627 1022380000 : updated = updated_seg .OR. updated_swap(iflavor)
1628 1022380000 : updated_swap(iflavor) = .FALSE.
1629 :
1630 1022380000 : CALL GreenHyb_measHybrid(this%Greens(iflavor), this%Bath%M(iflavor), this%Impurity%Particles(iflavor), updated)
1631 1022380000 : CALL Ctqmc_measN (this, iflavor, updated)
1632 1022380000 : IF ( this%opt_analysis .EQ. 1 ) &
1633 84840000 : CALL Ctqmc_measCorrelation (this, iflavor)
1634 1022380000 : IF ( this%opt_order .GT. 0 ) &
1635 2080950000 : CALL Ctqmc_measPerturbation(this, iflavor)
1636 : END DO
1637 :
1638 78610000 : IF ( MOD(isweep,modGlobalMove) .EQ. 0 ) THEN
1639 131386 : CALL Ctqmc_trySwap(this,swapUpdate1, swapUpdate2)
1640 131386 : IF ( swapUpdate1 .NE. 0 .AND. swapUpdate2 .NE. 0 ) THEN
1641 95797 : updated_swap(swapUpdate1) = .TRUE.
1642 95797 : updated_swap(swapUpdate2) = .TRUE.
1643 : END IF
1644 : END IF
1645 :
1646 78610000 : IF ( MOD(isweep,measurements) .EQ. 0 ) THEN
1647 78610000 : CALL ImpurityOperator_measDE(this%Impurity,this%measDE)
1648 78610000 : IF ( this%opt_spectra .GE. 1) THEN
1649 0 : IF ( measurements*this%opt_spectra /= 0) THEN
1650 0 : IF ( MOD(isweep,measurements*this%opt_spectra) .EQ. 0 ) THEN
1651 0 : this%density(1:flavors,indDensity) = this%measN(3,1:flavors)
1652 0 : indDensity = indDensity+1
1653 : END IF
1654 : END IF
1655 : END IF
1656 : END IF
1657 :
1658 78610000 : IF ( MOD(isweep,measurements) .EQ. 0 ) THEN
1659 78610000 : IF ( this%opt_histo .GT. 0 ) THEN
1660 : CALL ImpurityOperator_occup_histo_time(this%Impurity,this%occup_histo_time,this%occupconfig,this%suscep,this%samples,&
1661 3030000 : & this%chi,this%chicharge,this%ntot,this%opt_histo,this%nspinor)
1662 : ENDIF
1663 : ENDIF
1664 :
1665 78610000 : IF ( MOD(isweep, modNoise1) .EQ. 0 ) THEN
1666 : !modNext = isweep + modNoise2
1667 259963 : NRJ_new = (SUM(this%measDE(:,:))-this%measDE(1,1))*0.5d0 ! double occupation, avoid stat with 0 for U=J=0
1668 2049 : CALL Vector_pushBack(this%measNoise(1),NRJ_new - NRJ_old1)
1669 2049 : NRJ_old1 = NRJ_new
1670 :
1671 : !! Try to limit accumulation error
1672 2049 : CALL ImpurityOperator_cleanOverlaps(this%Impurity)
1673 :
1674 2049 : IF ( this%opt_noise .EQ. 1 ) THEN
1675 0 : DO iflavor = 1, flavors
1676 0 : DO ind = 1, this%Greens(iflavor)%this%tail
1677 0 : itau = this%Greens(iflavor)%this%listINT(ind)
1678 : gtmp_new(itau,iflavor) = this%Greens(iflavor)%oper(itau) &
1679 0 : +this%Greens(iflavor)%this%listDBLE(ind)*DBLE(this%Greens(iflavor)%factor)
1680 : END DO
1681 0 : DO itau = 1, sp1
1682 0 : CALL Vector_pushBack(this%measNoiseG(itau,iflavor,1), gtmp_new(itau,iflavor) - gtmp_old1(itau,iflavor))
1683 0 : gtmp_old1(itau,iflavor) = gtmp_new(itau,iflavor)
1684 : END DO
1685 : END DO
1686 : END IF
1687 : END IF
1688 :
1689 78610000 : IF ( MOD(isweep,modNoise2) .EQ. 0 ) THEN
1690 26045764 : NRJ_new = (SUM(this%measDE(:,:))-this%measDE(1,1))*0.5d0 ! double occupation, avoid stat with 0 for U=J=0
1691 205324 : CALL Vector_pushBack(this%measNoise(2),NRJ_new - NRJ_old2)
1692 205324 : NRJ_old2 = NRJ_new
1693 205324 : IF ( this%opt_noise .EQ. 1 ) THEN
1694 0 : DO iflavor = 1, flavors
1695 0 : DO ind = 1, this%Greens(iflavor)%this%tail
1696 0 : itau = this%Greens(iflavor)%this%listINT(ind)
1697 : gtmp_new(itau,iflavor) = this%Greens(iflavor)%oper(itau) &
1698 0 : +this%Greens(iflavor)%this%listDBLE(ind)*this%Greens(iflavor)%factor
1699 : END DO
1700 0 : DO itau = 1, sp1
1701 0 : CALL Vector_pushBack(this%measNoiseG(itau,iflavor,2), gtmp_new(itau,iflavor) - gtmp_old2(itau,iflavor))
1702 0 : gtmp_old2(itau,iflavor) = gtmp_new(itau,iflavor)
1703 : END DO
1704 : END DO
1705 : END IF
1706 :
1707 205324 : IF ( this%rank .EQ. 0 .AND. this%tid .EQ. 1 ) THEN
1708 146431 : new_percent = CEILING(DBLE(isweep)*100.d0/DBLE(itotal))
1709 150211 : DO ipercent = old_percent+1, new_percent
1710 150211 : WRITE(this%ostream,'(A)',ADVANCE="NO") "-"
1711 : END DO
1712 : old_percent = new_percent
1713 : END IF
1714 : END IF
1715 :
1716 78610098 : IF ( this%opt_movie .EQ. 1 ) THEN
1717 0 : WRITE(ilatex,'(A11,I9)') "%iteration ", isweep
1718 0 : CALL ImpurityOperator_printLatex(this%Impurity,ilatex,isweep)
1719 : END IF
1720 :
1721 : END DO
1722 :
1723 98 : IF ( this%rank .EQ. 0 .AND. this%tid .EQ. 1 ) THEN
1724 58 : DO ipercent = old_percent+1, 100
1725 58 : WRITE(this%ostream,'(A)',ADVANCE="NO") "-"
1726 : END DO
1727 38 : WRITE(this%ostream,'(A)') "|"
1728 : END IF
1729 :
1730 98 : FREE(gtmp_new)
1731 98 : FREE(gtmp_old1)
1732 98 : FREE(gtmp_old2)
1733 98 : FREE(updated_swap)
1734 :
1735 98 : IF ( this%opt_spectra .GE. 1 .AND. itotal .EQ. this%sweeps ) THEN
1736 0 : IF ( endDensity .NE. indDensity-1 ) THEN
1737 0 : this%density(:,endDensity) = -1.d0
1738 : END IF
1739 : END IF
1740 :
1741 98 : CALL CPU_TIME(cpu_time2)
1742 :
1743 98 : this%runTime = (cpu_time2 - cpu_time1)*1.05d0 ! facteur arbitraire de correction
1744 98 : END SUBROUTINE Ctqmc_loop
1745 : !!***
1746 :
1747 : !!****f* ABINIT/m_Ctqmc/Ctqmc_tryAddRemove
1748 : !! NAME
1749 : !! Ctqmc_tryAddRemove
1750 : !!
1751 : !! FUNCTION
1752 : !! Try to add or remove a segment and an anti-segment
1753 : !!
1754 : !! COPYRIGHT
1755 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
1756 : !! This file is distributed under the terms of the
1757 : !! GNU General Public License, see ~abinit/COPYING
1758 : !! or http://www.gnu.org/copyleft/gpl.txt .
1759 : !!
1760 : !! INPUTS
1761 : !! this=ctqmc
1762 : !!
1763 : !! OUTPUT
1764 : !! updated=something changed
1765 : !!
1766 : !! SIDE EFFECTS
1767 : !!
1768 : !! NOTES
1769 : !!
1770 : !! SOURCE
1771 :
1772 1022380000 : SUBROUTINE Ctqmc_tryAddRemove(this,updated)
1773 :
1774 : !Arguments ------------------------------------
1775 : TYPE(Ctqmc) , INTENT(INOUT) :: this
1776 : ! TYPE(BathOperator) , INTENT(INOUT) :: Bath
1777 : ! TYPE(ImpurityOperator), INTENT(INOUT) :: Impurity
1778 : LOGICAL , INTENT( OUT) :: updated
1779 : !Local variables ------------------------------
1780 : INTEGER :: position
1781 : INTEGER , DIMENSION(1:2) :: nature ! -2 for antiseg and 1 for seg
1782 : INTEGER :: i! -2 for antiseg and 1 for seg
1783 : DOUBLE PRECISION :: action
1784 : DOUBLE PRECISION :: beta
1785 : DOUBLE PRECISION :: time1
1786 : DOUBLE PRECISION :: time2
1787 : DOUBLE PRECISION :: time_avail
1788 : DOUBLE PRECISION :: det_ratio
1789 : DOUBLE PRECISION :: imp_trace
1790 : DOUBLE PRECISION :: signe
1791 : DOUBLE PRECISION :: tail
1792 : DOUBLE PRECISION, DIMENSION(1:2) :: CdagC_1
1793 :
1794 1022380000 : IF ( .NOT. this%set ) &
1795 0 : CALL ERROR("Ctqmc_trySegment : QMC not set ")
1796 :
1797 1022380000 : nature(1) = CTQMC_SEGME
1798 1022380000 : nature(2) = CTQMC_ANTIS
1799 1022380000 : beta = this%beta
1800 :
1801 1022380000 : updated = .FALSE.
1802 1022380000 : tail = DBLE(this%Impurity%particles(this%Impurity%activeFlavor)%tail)
1803 :
1804 :
1805 3067140000 : DO i = 1, 2
1806 2044760000 : signe = SIGN(1.d0,DBLE(nature(i)))
1807 :
1808 : !CALL RANDOM_NUMBER(action)
1809 2044760000 : CALL OurRng(this%seed,action)
1810 :
1811 3067140000 : IF ( action .LT. .5d0 ) THEN ! Ajout de segment
1812 : !CALL RANDOM_NUMBER(time1)
1813 1022389057 : CALL OurRng(this%seed,time1)
1814 1022389057 : time1 = time1 * beta
1815 1022389057 : time_avail = ImpurityOperator_getAvailableTime(this%Impurity,time1,position) * signe
1816 1022389057 : IF ( time_avail .GT. 0.d0 ) THEN
1817 : !CALL RANDOM_NUMBER(time2)
1818 511203414 : CALL OurRng(this%seed,time2)
1819 511203414 : IF ( time2 .EQ. 0.d0 ) THEN
1820 2 : CALL OurRng(this%seed,time2) ! Prevent null segment
1821 : END IF
1822 511203414 : time2 = time1 + time2 * time_avail
1823 511203414 : CdagC_1(Cdag_) = ((1.d0+signe)*time1+(1.d0-signe)*time2)*0.5d0
1824 511203414 : CdagC_1(C_ ) = ((1.d0+signe)*time2+(1.d0-signe)*time1)*0.5d0
1825 511203414 : det_ratio = BathOperator_getDetAdd(this%Bath,CdagC_1,position,this%Impurity%particles(this%Impurity%activeFlavor))
1826 511203414 : imp_trace = ImpurityOperator_getTraceAdd(this%Impurity,CdagC_1)
1827 : !CALL RANDOM_NUMBER(time1)
1828 511203414 : CALL OurRng(this%seed,time1)
1829 511203414 : IF ( det_ratio*imp_trace .LT. 0.d0 ) THEN
1830 197257 : this%stats(nature(i)+CTQMC_DETSI) = this%stats(nature(i)+CTQMC_DETSI) + 1.d0
1831 : END IF
1832 511203414 : IF ( (time1 * (tail + 1.d0 )) &
1833 : .LT. (beta * time_avail * det_ratio * imp_trace ) ) THEN
1834 85920850 : CALL ImpurityOperator_add(this%Impurity,CdagC_1,position)
1835 85920850 : CALL BathOperator_setMAdd(this%bath,this%Impurity%particles(this%Impurity%activeFlavor))
1836 85920850 : this%stats(nature(i)+CTQMC_ADDED) = this%stats(nature(i)+CTQMC_ADDED) + 1.d0
1837 85920850 : updated = .TRUE. .OR. updated
1838 85920850 : tail = tail + 1.d0
1839 : END IF
1840 : END IF
1841 :
1842 : ELSE ! Supprimer un segment
1843 1022370943 : IF ( tail .GT. 0.d0 ) THEN
1844 : !CALL RANDOM_NUMBER(time1)
1845 573945019 : CALL OurRng(this%seed,time1)
1846 573945019 : position = INT(((time1 * tail) + 1.d0) * signe )
1847 573945019 : time_avail = ImpurityOperator_getAvailedTime(this%Impurity,position)
1848 573945019 : det_ratio = BathOperator_getDetRemove(this%Bath,position)
1849 573945019 : imp_trace = ImpurityOperator_getTraceRemove(this%Impurity,position)
1850 : !CALL RANDOM_NUMBER(time1)
1851 573945019 : CALL OurRng(this%seed,time1)
1852 573945019 : IF ( det_ratio * imp_trace .LT. 0.d0 ) THEN
1853 15319 : this%stats(nature(i)+CTQMC_DETSI) = this%stats(nature(i)+CTQMC_DETSI) + 1.d0
1854 : END IF
1855 573945019 : IF ( (time1 * beta * time_avail ) &
1856 : .LT. (tail * det_ratio * imp_trace) ) THEN
1857 85918401 : CALL ImpurityOperator_remove(this%Impurity,position)
1858 85918401 : CALL BathOperator_setMRemove(this%Bath,this%Impurity%particles(this%Impurity%activeFlavor))
1859 85918401 : this%stats(nature(i)+CTQMC_REMOV) = this%stats(nature(i)+CTQMC_REMOV) + 1.d0
1860 85918401 : updated = .TRUE. .OR. updated
1861 85918401 : tail = tail -1.d0
1862 : END IF
1863 : END IF
1864 : END IF
1865 : END DO
1866 1022380000 : END SUBROUTINE Ctqmc_tryAddRemove
1867 : !!***
1868 :
1869 : !!****f* ABINIT/m_Ctqmc/Ctqmc_trySwap
1870 : !! NAME
1871 : !! Ctqmc_trySwap
1872 : !!
1873 : !! FUNCTION
1874 : !! try a global move (swap to flavors)
1875 : !!
1876 : !! COPYRIGHT
1877 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
1878 : !! This file is distributed under the terms of the
1879 : !! GNU General Public License, see ~abinit/COPYING
1880 : !! or http://www.gnu.org/copyleft/gpl.txt .
1881 : !!
1882 : !! INPUTS
1883 : !! this=ctqmc
1884 : !!
1885 : !! OUTPUT
1886 : !! flav_i=first flavor swaped
1887 : !! flav_j=second flavor swaped
1888 : !!
1889 : !! SIDE EFFECTS
1890 : !!
1891 : !! NOTES
1892 : !!
1893 : !! SOURCE
1894 :
1895 262772 : SUBROUTINE Ctqmc_trySwap(this,flav_i,flav_j)
1896 :
1897 : !Arguments ------------------------------------
1898 : TYPE(Ctqmc) , INTENT(INOUT) :: this
1899 : ! TYPE(BathOperator) , INTENT(INOUT) :: Bath
1900 : ! TYPE(ImpurityOperator), INTENT(INOUT) :: Impurity
1901 : INTEGER , INTENT( OUT) :: flav_i
1902 : INTEGER , INTENT( OUT) :: flav_j
1903 : !Local variables ------------------------------
1904 : INTEGER :: flavor_i
1905 : INTEGER :: flavor_j
1906 : DOUBLE PRECISION :: rnd
1907 : DOUBLE PRECISION :: lengthi
1908 : DOUBLE PRECISION :: lengthj
1909 : DOUBLE PRECISION :: overlapic1
1910 : DOUBLE PRECISION :: overlapjc1
1911 : DOUBLE PRECISION :: overlapic2
1912 : DOUBLE PRECISION :: overlapjc2
1913 : DOUBLE PRECISION :: detic1
1914 : DOUBLE PRECISION :: detjc1
1915 : DOUBLE PRECISION :: detic2
1916 : DOUBLE PRECISION :: detjc2
1917 : DOUBLE PRECISION :: det_ratio
1918 : DOUBLE PRECISION :: local_ratio
1919 :
1920 : !CALL RANDOM_NUMBER(rnd)
1921 131386 : CALL OurRng(this%seed,rnd)
1922 131386 : flavor_i = NINT(rnd*DBLE(this%flavors-1.d0))+1
1923 : !CALL RANDOM_NUMBER(rnd)
1924 131386 : CALL OurRng(this%seed,rnd)
1925 131386 : flavor_j = NINT(rnd*DBLE(this%flavors-1.d0))+1
1926 :
1927 131386 : flav_i = 0
1928 131386 : flav_j = 0
1929 :
1930 131386 : IF ( flavor_i .NE. flavor_j ) THEN
1931 : ! On tente d'intervertir i et j
1932 : ! Configuration actuelle :
1933 114290 : this%modGlobalMove(2) = this%modGlobalMove(2)+1
1934 114290 : detic1 = BathOperator_getDetF(this%Bath,flavor_i)
1935 114290 : detjc1 = BathOperator_getDetF(this%Bath,flavor_j)
1936 114290 : lengthi = ImpurityOperator_measN(this%Impurity,flavor_i)
1937 114290 : lengthj = ImpurityOperator_measN(this%Impurity,flavor_j)
1938 114290 : overlapic1 = ImpurityOperator_overlapFlavor(this%Impurity,flavor_i)
1939 114290 : overlapjc1 = ImpurityOperator_overlapFlavor(this%Impurity,flavor_j)
1940 : ! Configuration nouvelle :
1941 114290 : detic2 = BathOperator_getDetF(this%Bath,flavor_i,this%Impurity%particles(flavor_j))
1942 114290 : detjc2 = BathOperator_getDetF(this%Bath,flavor_j,this%Impurity%particles(flavor_i))
1943 : ! lengths unchanged
1944 114290 : overlapic2 = ImpurityOperator_overlapSwap(this%Impurity,flavor_i,flavor_j)
1945 114290 : overlapjc2 = ImpurityOperator_overlapSwap(this%Impurity,flavor_j,flavor_i)
1946 :
1947 : ! IF ( detic1*detjc1 .EQ. detic2*detjc2 ) THEN
1948 : ! det_ratio = 1.d0
1949 : ! ELSE IF ( detic1*detjc1 .EQ. 0.d0 ) THEN
1950 : ! det_ratio = detic2*detjc2 ! evite de diviser par 0 si pas de segment
1951 : ! ELSE
1952 114290 : det_ratio = detic2*detjc2/(detic1*detjc1)
1953 : ! END IF
1954 : local_ratio = DEXP(-overlapic2*overlapjc2+overlapic1*overlapjc1 &
1955 114290 : +(lengthj-lengthi)*(this%mu(flavor_i)-this%mu(flavor_j)))
1956 : ! Wloc = exp(muN-Uo)
1957 : !CALL RANDOM_NUMBER(rnd)
1958 114290 : CALL OurRng(this%seed,rnd)
1959 114290 : IF ( rnd .LT. local_ratio*det_ratio ) THEN ! swap accepted
1960 95797 : CALL ImpurityOperator_swap(this%Impurity, flavor_i,flavor_j)
1961 95797 : CALL BathOperator_swap (this%Bath , flavor_i,flavor_j)
1962 95797 : this%swap = this%swap + 1.d0
1963 95797 : flav_i = flavor_i
1964 95797 : flav_j = flavor_j
1965 : ! ELSE
1966 : ! CALL WARN("Swap refused")
1967 : ! WRITE(this%ostream,'(6E24.14)') local_ratio, det_ratio, detic1, detjc1, detic2, detjc2
1968 : END IF
1969 : END IF
1970 :
1971 131386 : END SUBROUTINE Ctqmc_trySwap
1972 : !!***
1973 :
1974 : !!****f* ABINIT/m_Ctqmc/Ctqmc_measN
1975 : !! NAME
1976 : !! Ctqmc_measN
1977 : !!
1978 : !! FUNCTION
1979 : !! measure the number of electron
1980 : !!
1981 : !! COPYRIGHT
1982 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
1983 : !! This file is distributed under the terms of the
1984 : !! GNU General Public License, see ~abinit/COPYING
1985 : !! or http://www.gnu.org/copyleft/gpl.txt .
1986 : !!
1987 : !! INPUTS
1988 : !! this=ctqmc
1989 : !! iflavor=which flavor to measure
1990 : !! updated=something has changed since last time
1991 : !!
1992 : !! OUTPUT
1993 : !!
1994 : !! SIDE EFFECTS
1995 : !!
1996 : !! NOTES
1997 : !!
1998 : !! SOURCE
1999 :
2000 1022380000 : SUBROUTINE Ctqmc_measN(this, iflavor, updated)
2001 :
2002 : !Arguments ------------------------------------
2003 : TYPE(Ctqmc) , INTENT(INOUT) :: this
2004 : !TYPE(ImpurityOperator), INTENT(IN ) :: impurity
2005 : INTEGER , INTENT(IN ) :: iflavor
2006 : LOGICAL , INTENT(IN ) :: updated
2007 :
2008 : ! IF ( .NOT. this%set ) &
2009 : ! CALL ERROR("Ctqmc_measN : QMC not set ")
2010 :
2011 :
2012 1022380000 : IF ( updated .EQV. .TRUE. ) THEN
2013 158969926 : this%measN(1,iflavor) = this%measN(1,iflavor) + this%measN(3,iflavor)*this%measN(4,iflavor)
2014 158969926 : this%measN(2,iflavor) = this%measN(2,iflavor) + this%measN(4,iflavor)
2015 158969926 : this%measN(3,iflavor) = ImpurityOperator_measN(this%impurity)
2016 158969926 : this%measN(4,iflavor) = 1.d0
2017 : ELSE
2018 863410074 : this%measN(4,iflavor) = this%measN(4,iflavor) + 1.d0
2019 : END IF
2020 1022380000 : END SUBROUTINE Ctqmc_measN
2021 : !!***
2022 :
2023 : !#ifdef CTCtqmc_ANALYSIS
2024 : !!****f* ABINIT/m_Ctqmc/Ctqmc_measCorrelation
2025 : !! NAME
2026 : !! Ctqmc_measCorrelation
2027 : !!
2028 : !! FUNCTION
2029 : !! measure all correlations in times for a flavor
2030 : !!
2031 : !! COPYRIGHT
2032 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
2033 : !! This file is distributed under the terms of the
2034 : !! GNU General Public License, see ~abinit/COPYING
2035 : !! or http://www.gnu.org/copyleft/gpl.txt .
2036 : !!
2037 : !! INPUTS
2038 : !! this=ctqmc
2039 : !! iflavor=the flavor to measure
2040 : !!
2041 : !! OUTPUT
2042 : !!
2043 : !! SIDE EFFECTS
2044 : !!
2045 : !! NOTES
2046 : !!
2047 : !! SOURCE
2048 :
2049 84840000 : SUBROUTINE Ctqmc_measCorrelation(this, iflavor)
2050 :
2051 : !Arguments ------------------------------------
2052 : TYPE(Ctqmc) , INTENT(INOUT) :: this
2053 : !TYPE(ImpurityOperator), INTENT(IN ) :: impurity
2054 : INTEGER , INTENT(IN ) :: iflavor
2055 : !Local variables ------------------------------
2056 : INTEGER :: iCdag
2057 : INTEGER :: iCdagBeta
2058 : INTEGER :: iC
2059 : INTEGER :: index
2060 : INTEGER :: size
2061 : DOUBLE PRECISION :: tC
2062 : DOUBLE PRECISION :: tCdag
2063 : !DOUBLE PRECISION :: time
2064 : DOUBLE PRECISION :: inv_dt
2065 : DOUBLE PRECISION :: beta
2066 :
2067 84840000 : IF ( .NOT. this%set ) &
2068 0 : CALL ERROR("Ctqmc_measCorrelation : QMC not set ")
2069 :
2070 84840000 : size = this%impurity%particles(this%impurity%activeFlavor)%tail
2071 84840000 : beta = this%beta
2072 :
2073 84840000 : IF ( size .EQ. 0 ) RETURN
2074 :
2075 35292468 : inv_dt = this%inv_dt
2076 :
2077 113957624 : DO iCdag = 1, size ! first segments
2078 78665156 : tCdag = this%impurity%particles(this%impurity%activeFlavor)%list(iCdag,Cdag_)
2079 78665156 : tC = this%impurity%particles(this%impurity%activeFlavor)%list(iCdag,C_ )
2080 78665156 : index = INT( ( (tC - tCdag) * inv_dt ) + .5d0 ) + 1
2081 78665156 : this%measCorrelation(index,1,iflavor) = this%measCorrelation(index,1,iflavor) + 1.d0
2082 78665156 : MODCYCLE(iCdag+1,size,iCdagBeta)
2083 : index = INT( ( ( &
2084 : this%impurity%particles(this%impurity%activeFlavor)%list(iCdagBeta,Cdag_) - tC &
2085 : + AINT(DBLE(iCdag)/DBLE(size))*beta &
2086 78665156 : ) * inv_dt ) + .5d0 ) + 1
2087 78665156 : IF ( index .LT. 1 .OR. index .GT. this%samples+1 ) THEN
2088 0 : CALL WARN("Ctqmc_measCorrelation : bad index line 1095 ")
2089 : ELSE
2090 78665156 : this%measCorrelation(index,2,iflavor) = this%measCorrelation(index,2,iflavor) + 1.d0
2091 : END IF
2092 : ! DO iC = 1, size
2093 : ! tC = impurity%particles(impurity%activeFlavor)%list(C_,iC)
2094 : ! time = tC - tCdag
2095 : ! IF ( time .LT. 0.d0 ) time = time + beta
2096 : ! index = INT( ( time * inv_dt ) + .5d0 ) + 1
2097 : ! this%measCorrelation(index,3,iflavor) = this%measCorrelation(index,3,iflavor) + 1.d0
2098 : ! END DO
2099 396182224 : DO iC = 1, size! this%Greens(iflavor)%index_old%tail
2100 : this%measCorrelation(this%Greens(iflavor)%this%listINT(iC+(iCdag-1)*size),3,iflavor) = &
2101 360889756 : this%measCorrelation(this%Greens(iflavor)%this%listINT(iC+(iCdag-1)*size),3,iflavor) + 1.d0
2102 : END DO
2103 : END DO
2104 :
2105 : END SUBROUTINE Ctqmc_measCorrelation
2106 : !!***
2107 :
2108 : !!****f* ABINIT/m_Ctqmc/Ctqmc_measPerturbation
2109 : !! NAME
2110 : !! Ctqmc_measPerturbation
2111 : !!
2112 : !! FUNCTION
2113 : !! measure perturbation order
2114 : !!
2115 : !! COPYRIGHT
2116 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
2117 : !! This file is distributed under the terms of the
2118 : !! GNU General Public License, see ~abinit/COPYING
2119 : !! or http://www.gnu.org/copyleft/gpl.txt .
2120 : !!
2121 : !! INPUTS
2122 : !! this=ctqmc
2123 : !! iflavor=the flavor to measure
2124 : !!
2125 : !! OUTPUT
2126 : !!
2127 : !! SIDE EFFECTS
2128 : !!
2129 : !! NOTES
2130 : !!
2131 : !! SOURCE
2132 :
2133 979960000 : SUBROUTINE Ctqmc_measPerturbation(this, iflavor)
2134 :
2135 : !Arguments ------------------------------------
2136 : TYPE(Ctqmc) , INTENT(INOUT) :: this
2137 : !TYPE(ImpurityOperator), INTENT(IN ) :: impurity
2138 : INTEGER , INTENT(IN ) :: iflavor
2139 : !Local variables ------------------------------
2140 : INTEGER :: index
2141 :
2142 979960000 : IF ( .NOT. this%set ) &
2143 0 : CALL ERROR("Ctqmc_measiPerturbation : QMC not set ")
2144 :
2145 979960000 : index = this%impurity%particles(this%impurity%activeFlavor)%tail + 1
2146 979960000 : IF ( index .LE. this%opt_order ) &
2147 979960000 : this%measPerturbation(index,iflavor) = this%measPerturbation(index,iflavor) + 1.d0
2148 :
2149 979960000 : END SUBROUTINE Ctqmc_measPerturbation
2150 : !!***
2151 : !#endif
2152 :
2153 : !!****f* ABINIT/m_Ctqmc/Ctqmc_getResult
2154 : !! NAME
2155 : !! Ctqmc_getResult
2156 : !!
2157 : !! FUNCTION
2158 : !! reduce everything to get the result of the simulation
2159 : !!
2160 : !! COPYRIGHT
2161 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder,F. Gendron)
2162 : !! This file is distributed under the terms of the
2163 : !! GNU General Public License, see ~abinit/COPYING
2164 : !! or http://www.gnu.org/copyleft/gpl.txt .
2165 : !!
2166 : !! INPUTS
2167 : !! this=ctqmc
2168 : !!
2169 : !! OUTPUT
2170 : !!
2171 : !! SIDE EFFECTS
2172 : !!
2173 : !! NOTES
2174 : !!
2175 : !! SOURCE
2176 :
2177 34 : SUBROUTINE Ctqmc_getResult(this,num_chains,chains,Iatom,fname)
2178 :
2179 :
2180 : #ifdef HAVE_MPI1
2181 : include 'mpif.h'
2182 : #endif
2183 : !Arguments ------------------------------------
2184 : TYPE(Ctqmc) , INTENT(INOUT) :: this
2185 : INTEGER, INTENT(IN ) :: Iatom, num_chains
2186 : TYPE(Ctqmc) , INTENT(INOUT) :: chains(num_chains)
2187 : character(len=fnlen), INTENT(IN) :: fname
2188 : !Local variables ------------------------------
2189 : INTEGER :: iflavor
2190 : INTEGER :: flavors
2191 : INTEGER :: itau
2192 : INTEGER :: endDensity
2193 : DOUBLE PRECISION :: inv_flavors
2194 : DOUBLE PRECISION :: a
2195 : DOUBLE PRECISION :: b
2196 : DOUBLE PRECISION :: r
2197 34 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: alpha
2198 34 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: beta
2199 : DOUBLE PRECISION, DIMENSION(1:2) :: TabX
2200 : DOUBLE PRECISION, DIMENSION(1:2) :: TabY
2201 34 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: freqs
2202 34 : INTEGER, ALLOCATABLE, DIMENSION(:) :: counts
2203 34 : INTEGER, ALLOCATABLE, DIMENSION(:) :: displs
2204 34 : INTEGER, ALLOCATABLE, DIMENSION(:) :: occtot
2205 34 : INTEGER, ALLOCATABLE, DIMENSION(:) :: spintot
2206 34 : INTEGER, ALLOCATABLE, DIMENSION(:,:) :: occ
2207 : INTEGER :: sp1,spinmax,spinmin,dspin,nelec,spin
2208 : INTEGER :: spAll
2209 : INTEGER :: last
2210 : INTEGER :: n1
2211 : INTEGER :: n2,n3,quotient,remainder,signe
2212 : INTEGER :: debut,i
2213 : ! INTEGER :: fin
2214 : character(len=2) :: atomnb
2215 : ! character(len=fnlen) :: tmpfile
2216 : ! INTEGER :: unt
2217 : #ifdef HAVE_MPI
2218 : INTEGER :: ierr
2219 : DOUBLE PRECISION, DIMENSION(1) :: rtime
2220 : #endif
2221 : DOUBLE PRECISION :: inv_size,sumh,sumtot
2222 34 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: buffer
2223 34 : TYPE(FFTHyb) :: FFTmrka
2224 :
2225 : #if defined HAVE_MPI && !defined HAVE_MPI2_INPLACE
2226 : DOUBLE PRECISION, ALLOCATABLE , DIMENSION(:) :: buffer1_out,freqs_buf
2227 : DOUBLE PRECISION, ALLOCATABLE , DIMENSION(:,:) :: buffer2_out
2228 : #endif
2229 :
2230 34 : IF ( .NOT. this%done ) &
2231 0 : CALL ERROR("Ctqmc_getResult : Simulation not run ")
2232 :
2233 34 : flavors = this%flavors
2234 34 : inv_flavors = 1.d0 / DBLE(flavors)
2235 :
2236 :
2237 34 : inv_size = 1.d0 / DBLE(this%size*num_chains)
2238 34 : sp1 = 0
2239 34 : spAll = 0
2240 :
2241 : !#ifdef CTCtqmc_CHECK
2242 34 : IF ( this%opt_check .GT. 0 ) THEN
2243 0 : do i=1, num_chains
2244 0 : chains(i)%errorImpurity = ImpurityOperator_getError(chains(i)%Impurity) * inv_flavors
2245 0 : chains(i)%errorBath = BathOperator_getError (chains(i)%Bath ) * inv_flavors
2246 : end do
2247 : END IF
2248 : !#endif
2249 :
2250 34 : MALLOC(alpha,(1,1))
2251 34 : MALLOC(beta,(1,1))
2252 34 : MALLOC(buffer,(1,1))
2253 34 : IF ( this%opt_noise .EQ. 1) THEN
2254 0 : FREEIF(alpha)
2255 0 : MALLOC(alpha,(1:this%samples+1,1:flavors))
2256 0 : FREEIF(beta)
2257 0 : MALLOC(beta,(1:this%samples+1,1:flavors))
2258 : END IF
2259 :
2260 34 : IF ( this%have_MPI .EQV. .TRUE.) THEN
2261 34 : sp1 = this%samples+1
2262 34 : spALL = sp1 + flavors + 6
2263 :
2264 : !#ifdef CTCtqmc_ANALYSIS
2265 34 : IF ( this%opt_analysis .EQ. 1 ) &
2266 6 : spAll = spAll + 3*sp1
2267 34 : IF ( this%opt_order .GT. 0 ) &
2268 27 : spAll = spAll + this%opt_order
2269 34 : IF ( this%opt_noise .EQ. 1 ) &
2270 0 : spAll = spAll + 2*(this%samples + 1)
2271 : !#endif
2272 :
2273 34 : FREEIF(buffer)
2274 136 : MALLOC(buffer,(1:spAll,1:MAX(2,flavors)))
2275 113326 : buffer = 0.0d0
2276 : END IF
2277 :
2278 : ! this%seg_added = this%seg_added * inv_flavors
2279 : ! this%seg_removed = this%seg_removed * inv_flavors
2280 : ! this%seg_sign = this%seg_sign * inv_flavors
2281 : ! this%anti_added = this%anti_added * inv_flavors
2282 : ! this%anti_removed = this%anti_removed * inv_flavors
2283 : ! this%anti_sign = this%anti_sign * inv_flavors
2284 :
2285 83 : DO i = 1, num_chains
2286 377 : chains(i)%stats(:) = chains(i)%stats(:) * inv_flavors
2287 : END DO
2288 83 : DO i = 1, num_chains
2289 567 : DO iflavor = 1, flavors
2290 518 : CALL GreenHyb_measHybrid(chains(i)%Greens(iflavor), chains(i)%Bath%M(iflavor), chains(i)%Impurity%Particles(iflavor), .TRUE.)
2291 518 : CALL GreenHyb_getHybrid(chains(i)%Greens(iflavor))
2292 : ! Accumule les dernieres mesure de N
2293 518 : chains(i)%measN(1,iflavor) = chains(i)%measN(1,iflavor) + chains(i)%measN(3,iflavor)*chains(i)%measN(4,iflavor)
2294 518 : chains(i)%measN(2,iflavor) = chains(i)%measN(2,iflavor) + chains(i)%measN(4,iflavor)
2295 : ! Reduction
2296 518 : chains(i)%measN(1,iflavor) = chains(i)%measN(1,iflavor) / ( chains(i)%measN(2,iflavor) * chains(i)%beta )
2297 : ! Correction
2298 518 : CALL GreenHyb_setN(chains(i)%Greens(iflavor), chains(i)%measN(1,iflavor))
2299 : !#ifdef CTCtqmc_ANALYSIS
2300 518 : IF ( chains(i)%opt_order .GT. 0 ) &
2301 : chains(i)%measPerturbation(: ,iflavor) = chains(i)%measPerturbation(:,iflavor) &
2302 48076 : / SUM(chains(i)%measPerturbation(:,iflavor))
2303 518 : IF ( chains(i)%opt_analysis .EQ. 1 ) THEN
2304 : chains(i)%measCorrelation (:,1,iflavor) = chains(i)%measCorrelation (:,1,iflavor) &
2305 : / SUM(chains(i)%measCorrelation (:,1,iflavor)) &
2306 8652 : * chains(i)%inv_dt
2307 : chains(i)%measCorrelation (:,2,iflavor) = chains(i)%measCorrelation (:,2,iflavor) &
2308 : / SUM(chains(i)%measCorrelation (:,2,iflavor)) &
2309 8652 : * chains(i)%inv_dt
2310 : chains(i)%measCorrelation (:,3,iflavor) = chains(i)%measCorrelation (:,3,iflavor) &
2311 : / SUM(chains(i)%measCorrelation (:,3,iflavor)) &
2312 8652 : * chains(i)%inv_dt
2313 : END IF
2314 : !#endif
2315 518 : IF ( chains(i)%opt_noise .EQ. 1 ) THEN
2316 0 : TabX(1) = DBLE(chains(i)%modNoise2)
2317 0 : TabX(2) = DBLE(chains(i)%modNoise1)
2318 0 : DO itau = 1, chains(i)%samples+1
2319 : chains(i)%measNoiseG(itau,iflavor,2)%vec = -chains(i)%measNoiseG(itau,iflavor,2)%vec*chains(i)%inv_dt &
2320 0 : /(chains(i)%beta*DBLE(chains(i)%modNoise2))
2321 : chains(i)%measNoiseG(itau,iflavor,1)%vec = -chains(i)%measNoiseG(itau,iflavor,1)%vec*chains(i)%inv_dt &
2322 0 : /(chains(i)%beta*DBLE(chains(i)%modNoise1))
2323 0 : n2 = chains(i)%measNoiseG(itau,iflavor,2)%tail
2324 0 : TabY(1) = Stat_deviation(chains(i)%measNoiseG(itau,iflavor,2)%vec(1:n2))!*SQRT(n2/(n2-1))
2325 0 : n1 = chains(i)%measNoiseG(itau,iflavor,1)%tail
2326 0 : TabY(2) = Stat_deviation(chains(i)%measNoiseG(itau,iflavor,1)%vec(1:n1))!*SQRT(n1/(n1-1))
2327 0 : CALL Stat_powerReg(TabX,SQRT(2.d0*LOG(2.d0))*TabY,alpha(itau,iflavor),beta(itau,iflavor),r)
2328 : ! ecart type -> 60%
2329 : ! largeur a mi-hauteur d'une gaussienne -> sqrt(2*ln(2))*sigma
2330 : END DO
2331 : END IF
2332 :
2333 567 : IF ( chains(i)%have_MPI .EQV. .TRUE. ) THEN
2334 152936 : buffer(1:sp1, iflavor) = buffer(1:sp1, iflavor) + chains(i)%Greens(iflavor)%oper(1:sp1)
2335 : END IF
2336 : END DO
2337 6397 : chains(i)%measDE(:,:) = chains(i)%measDE(:,:) * DBLE(chains(i)%measurements) /(DBLE(chains(i)%sweeps)*chains(i)%beta)
2338 : END DO
2339 34 : last = sp1
2340 :
2341 34 : IF ( this%opt_histo .GT. 0 ) THEN
2342 : !Note: this == chains(1), so scale here and sum other chains
2343 48 : this%occup_histo_time(:) = chains(1)%occup_histo_time(:) / INT(chains(1)%sweeps/chains(1)%measurements)
2344 49155 : this%occupconfig(:) = chains(1)%occupconfig(:) / INT(chains(1)%sweeps/chains(1)%measurements)
2345 603 : this%suscep(:,:) = chains(1)%suscep(:,:) / INT(chains(1)%sweeps/chains(1)%measurements)
2346 603 : this%chi(:,:) = chains(1)%chi(:,:) / INT(chains(1)%sweeps/chains(1)%measurements)
2347 603 : this%chicharge(:,:) = chains(1)%chicharge(:,:) / INT(chains(1)%sweeps/chains(1)%measurements)
2348 12 : this%ntot(:) = chains(1)%ntot(:) / INT(chains(1)%sweeps/chains(1)%measurements)
2349 3 : do i = 2, num_chains
2350 : this%occup_histo_time(:) = this%occup_histo_time(:) &
2351 0 : & + chains(i)%occup_histo_time(:) / INT(chains(i)%sweeps/chains(i)%measurements)
2352 0 : this%occupconfig(:) = this%occupconfig(:) + chains(i)%occupconfig(:) / INT(chains(i)%sweeps/chains(i)%measurements)
2353 0 : this%suscep(:,:) = this%suscep(:,:) + chains(i)%suscep(:,:) / INT(chains(i)%sweeps/chains(i)%measurements)
2354 0 : this%chi(:,:) = this%chi(:,:) + chains(i)%chi(:,:) / INT(chains(i)%sweeps/chains(i)%measurements)
2355 0 : this%chicharge(:,:) = this%chicharge(:,:) + chains(i)%chicharge(:,:) / INT(chains(i)%sweeps/chains(i)%measurements)
2356 3 : this%ntot(:) = this%ntot(:) + chains(i)%ntot(:) / INT(chains(i)%sweeps/chains(i)%measurements)
2357 : end do
2358 : END IF
2359 : ! write(6,*) "=== Histogram of occupations for complete simulation ====",INT(this%sweeps/this%measurements)
2360 : ! sumh=0
2361 : ! do n1=1,this%flavors+1
2362 : ! write(6,'(i4,f10.4)') n1-1, this%occup_histo_time(n1)
2363 : ! sumh=sumh+this%occup_histo_time(n1)
2364 : ! enddo
2365 : ! write(6,*) "=================================",sumh
2366 :
2367 34 : n1 = this%measNoise(1)%tail
2368 34 : n2 = this%measNoise(2)%tail
2369 34 : CALL Vector_setSize(this%measNoise(1),n1*num_chains)
2370 34 : CALL Vector_setSize(this%measNoise(2),n2*num_chains)
2371 49 : do i = 2, num_chains
2372 162 : this%measNoise(1)%vec(n1*(i-1)+1:n1*i) = chains(i)%measNoise(1)%vec(1:n1)
2373 15049 : this%measNoise(2)%vec(n2*(i-1)+1:n2*i) = chains(i)%measNoise(2)%vec(1:n2)
2374 : end do
2375 34 : n1 = n1 * num_chains
2376 34 : n2 = n2 * num_chains
2377 :
2378 : ! On utilise freqs comme tableau de regroupement
2379 : ! Gather de Noise1
2380 34 : IF ( this%have_MPI .EQV. .TRUE. ) THEN
2381 102 : MALLOC(counts,(1:this%size))
2382 68 : MALLOC(displs,(1:this%size))
2383 : FREEIF(freqs)
2384 102 : MALLOC(freqs,(1:this%size*n1))
2385 3418 : freqs = 0.d0
2386 128 : counts(:) = n1
2387 222 : displs(:) = (/ ( iflavor*n1, iflavor=0, this%size-1 ) /)
2388 : #ifdef HAVE_MPI
2389 : #if defined HAVE_MPI2_INPLACE
2390 1930 : freqs(n1*this%rank+1:n1*(this%rank+1)) = this%measNoise(1)%vec(1:n1)
2391 : CALL MPI_ALLGATHERV(MPI_IN_PLACE, 0, MPI_DOUBLE_PRECISION, &
2392 : freqs, counts, displs, &
2393 34 : MPI_DOUBLE_PRECISION, this%MY_COMM, ierr)
2394 : #else
2395 : MALLOC(freqs_buf,(n1))
2396 : freqs_buf(1:n1)=this%measNoise(1)%vec(1:n1)
2397 : CALL MPI_ALLGATHERV(freqs_buf, n1, MPI_DOUBLE_PRECISION, &
2398 : freqs, counts, displs, &
2399 : MPI_DOUBLE_PRECISION, this%MY_COMM, ierr)
2400 : FREE(freqs_buf)
2401 : #endif
2402 : #endif
2403 34 : n1 = this%size*n1
2404 34 : CALL Vector_setSize(this%measNoise(1),n1)
2405 3418 : this%measNoise(1)%vec(1:n1) = freqs(:)
2406 : ! Gather de Noise2
2407 34 : FREE(freqs)
2408 102 : MALLOC(freqs,(1:this%size*n2))
2409 340034 : freqs = 0.d0
2410 128 : counts(:) = n2
2411 222 : displs(:) = (/ ( iflavor*n2, iflavor=0, this%size-1 ) /)
2412 : #ifdef HAVE_MPI
2413 : #if defined HAVE_MPI2_INPLACE
2414 190034 : freqs(n2*this%rank+1:n2*(this%rank+1)) = this%measNoise(2)%vec(1:n2)
2415 : CALL MPI_ALLGATHERV(MPI_IN_PLACE, 0, MPI_DOUBLE_PRECISION, &
2416 : freqs, counts, displs, &
2417 34 : MPI_DOUBLE_PRECISION, this%MY_COMM, ierr)
2418 : #else
2419 : MALLOC(freqs_buf,(n2))
2420 : freqs_buf(1:n2)=this%measNoise(2)%vec(1:n2)
2421 : CALL MPI_ALLGATHERV(freqs_buf, n2, MPI_DOUBLE_PRECISION, &
2422 : freqs, counts, displs, &
2423 : MPI_DOUBLE_PRECISION, this%MY_COMM, ierr)
2424 : FREE(freqs_buf)
2425 : #endif
2426 : #endif
2427 34 : n2 = this%size*n2
2428 34 : CALL Vector_setSize(this%measNoise(2),n2)
2429 340034 : this%measNoise(2)%vec(1:n2) = freqs(:)
2430 34 : FREE(counts)
2431 34 : FREE(displs)
2432 34 : FREE(freqs)
2433 : END IF
2434 : !n1 = this%measNoise(1)%tail
2435 : !n2 = this%measNoise(2)%tail
2436 :
2437 : ! Transformation des paquets pour que ca fit a CTQMC_SLICE(1|2)
2438 34 : IF ( n1 .GT. CTQMC_SLICE1 ) THEN
2439 0 : itau = n1/CTQMC_SLICE1
2440 0 : MALLOC(freqs,(1:n1/itau))
2441 0 : DO debut=1, n1/itau
2442 0 : freqs(debut)=SUM(this%measNoise(1)%vec((debut-1)*itau+1:itau*debut))
2443 : END DO
2444 0 : freqs(:) = freqs(:)/DBLE(itau)
2445 0 : this%modNoise1 = this%modNoise1*itau
2446 0 : n1 = n1/itau
2447 0 : CALL Vector_setSize(this%measNoise(1),n1)
2448 0 : this%measNoise(1)%vec(1:n1) = freqs(:)
2449 0 : FREE(freqs)
2450 : END IF
2451 34 : IF ( n2 .GT. CTQMC_SLICE1*CTQMC_SLICE2 ) THEN
2452 0 : itau = n2/(CTQMC_SLICE1*CTQMC_SLICE2)
2453 0 : MALLOC(freqs,(1:n2/itau))
2454 0 : DO debut=1, n2/itau
2455 0 : freqs(debut)=SUM(this%measNoise(2)%vec((debut-1)*itau+1:itau*debut))
2456 : END DO
2457 0 : freqs(:) = freqs(:)/DBLE(itau)
2458 0 : this%modNoise2 = this%modNoise2*itau
2459 0 : n2 = n2/itau
2460 0 : CALL Vector_setSize(this%measNoise(2),n2)
2461 0 : this%measNoise(2)%vec(1:n2) = freqs(:)
2462 0 : FREE(freqs)
2463 : END IF
2464 : ! On peut s'amuser avec nos valeur d'energies
2465 : !MALLOC(TabX,(1:20))
2466 : !MALLOC(TabY,(1:20))
2467 :
2468 34 : TabX(1) = DBLE(this%modNoise2)
2469 34 : TabX(2) = DBLE(this%modNoise1)
2470 :
2471 : ! Il faut calculer pour chaque modulo 10 ecarts type sur les donnes acquises
2472 3418 : this%measNoise(1)%vec(1:n1) = this%measNoise(1)%vec(1:n1)/(this%beta*DBLE(this%modNoise1))*DBLE(this%measurements)
2473 340034 : this%measNoise(2)%vec(1:n2) = this%measNoise(2)%vec(1:n2)/(this%beta*DBLE(this%modNoise2))*DBLE(this%measurements)
2474 : ! CALL Vector_print(this%measNoise(1),this%rank+70)
2475 : ! CALL Vector_print(this%measNoise(2),this%rank+50)
2476 : ! DO iflavor=1,10
2477 : ! debut = (iflavor-1)*n2/10+1
2478 : ! fin = iflavor*n2/10
2479 : ! TabY(iflavor) = Stat_deviation(this%measNoise(2)%vec(debut:fin))
2480 : ! debut = (iflavor-1)*n1/10+1
2481 : ! fin = iflavor*n1/10
2482 : ! TabY(10+iflavor) = Stat_deviation(this%measNoise(1)%vec(debut:fin))
2483 : ! END DO
2484 : !! TabY(1:n) = (this%measNoise(2)%vec(1:n) &
2485 : !! )
2486 : !! !/(this%beta*DBLE(this%modNoise2))*DBLE(this%measurements) &
2487 : !! !- this%measDE(1,1))
2488 : !! TabY(this%measNoise(2)%tail+1:n+this%measNoise(2)%tail) = (this%measNoise(1)%vec(1:n) &
2489 : !! )
2490 : !! ! /(this%beta*DBLE(this%modNoise1))*DBLE(this%measurements) &
2491 : !! ! - this%measDE(1,1))
2492 : ! IF ( this%rank .EQ. 0 ) THEN
2493 : ! DO iflavor=1,20
2494 : ! write(45,*) TabX(iflavor), TabY(iflavor)
2495 : ! END DO
2496 : ! END IF
2497 : !
2498 :
2499 :
2500 34 : TabY(1) = Stat_deviation(this%measNoise(2)%vec(1:n2))!*SQRT(n2/(n2-1))
2501 : !! write(this%rank+10,*) TabX(2)
2502 : !! write(this%rank+40,*) TabX(1)
2503 : !! CALL Vector_print(this%measNoise(1),this%rank+10)
2504 : !! CALL Vector_print(this%measNoise(2),this%rank+40)
2505 : !! CLOSE(this%rank+10)
2506 : !! CLOSE(this%rank+40)
2507 34 : TabY(2) = Stat_deviation(this%measNoise(1)%vec(1:n1))!*SQRT(n1/(n1-1))
2508 : !! ! Ecart carre moyen ~ ecart type mais non biaise. Serait moins precis. Aucun
2509 : ! impact sur la pente, juste sur l'ordonnee a l'origine.
2510 :
2511 102 : CALL Stat_powerReg(TabX,SQRT(2.d0*LOG(2.d0))*TabY,a,b,r)
2512 : ! FREE(TabX)
2513 : ! FREE(TabY)
2514 : ! ecart type -> 60%
2515 : ! largeur a mi-hauteur d'une gaussienne -> sqrt(2*ln(2))*sigma
2516 :
2517 : !this%measDE(1,1) = SUM(this%measNoise(1)%vec(1:this%measNoise(1)%tail))/(DBLE(this%measNoise(1)%tail*this%modNoise1)*this%beta)
2518 : !this%measDE(2:flavors,1:flavors) = this%measDE(2:flavors,1:flavors) /(DBLE(this%sweeps)*this%beta)
2519 : do i = i, num_chains
2520 : CALL ImpurityOperator_getErrorOverlap(chains(i)%Impurity,chains(i)%measDE)
2521 : end do
2522 :
2523 : ! Add the difference between true calculation and quick calculation of the
2524 : ! last sweep overlap to measDE(2,2)
2525 : !chains(i)%measDE = chains(i)%measDE * DBLE(chains(i)%measurements)
2526 34 : IF ( this%have_MPI .EQV. .TRUE. ) THEN
2527 83 : do i = 1, num_chains
2528 49 : last = sp1
2529 49 : IF ( chains(i)%opt_analysis .EQ. 1 ) THEN
2530 4374 : buffer(last+1:last+sp1,:) = buffer(last+1:last+sp1,:) + chains(i)%measCorrelation(:,1,:)
2531 6 : last = last + sp1
2532 4374 : buffer(last+1:last+sp1,:) = buffer(last+1:last+sp1,:) + chains(i)%measCorrelation(:,2,:)
2533 6 : last = last + sp1
2534 4374 : buffer(last+1:last+sp1,:) = buffer(last+1:last+sp1,:) + chains(i)%measCorrelation(:,3,:)
2535 : last = last + sp1
2536 : END IF
2537 49 : IF ( chains(i)%opt_order .GT. 0 ) THEN
2538 24318 : buffer(last+1:last+chains(i)%opt_order, :) = buffer(last+1:last+chains(i)%opt_order, :) + chains(i)%measPerturbation(:,:)
2539 : last = last + chains(i)%opt_order
2540 : END IF
2541 49 : IF ( chains(i)%opt_noise .EQ. 1 ) THEN
2542 0 : buffer(last+1:last+chains(i)%samples+1,:) = buffer(last+1:last+chains(i)%samples+1,:) + alpha(:,:)
2543 0 : last = last + chains(i)%samples + 1
2544 0 : buffer(last+1:last+chains(i)%samples+1,:) = buffer(last+1:last+chains(i)%samples+1,:) + beta(:,:)
2545 49 : last = last + chains(i)%samples + 1
2546 : END IF
2547 : ! chains(i)%measDE(2,2) = a*EXP(b*LOG(DBLE(chains(i)%sweeps*chains(i)%size)))
2548 6363 : buffer(spall-(flavors+5):spAll-6,:) = buffer(spall-(flavors+5):spAll-6,:) + chains(i)%measDE(:,:)
2549 : ! buffer(spAll ,1) = chains(i)%seg_added
2550 : ! buffer(spAll-1,1) = chains(i)%seg_removed
2551 : ! buffer(spAll-2,1) = chains(i)%seg_sign
2552 : ! buffer(spAll ,2) = chains(i)%anti_added
2553 : ! buffer(spAll-1,2) = chains(i)%anti_removed
2554 : ! buffer(spAll-2,2) = chains(i)%anti_sign
2555 49 : buffer(spAll ,1) = buffer(spAll ,1) + chains(i)%stats(1)
2556 49 : buffer(spAll-1,1) = buffer(spAll-1,1) + chains(i)%stats(2)
2557 49 : buffer(spAll-2,1) = buffer(spAll-2,1) + chains(i)%stats(3)
2558 49 : buffer(spAll ,2) = buffer(spAll ,2) + chains(i)%stats(4)
2559 49 : buffer(spAll-1,2) = buffer(spAll-1,2) + chains(i)%stats(5)
2560 49 : buffer(spAll-2,2) = buffer(spAll-2,2) + chains(i)%stats(6)
2561 49 : buffer(spAll-3,1) = buffer(spAll-3,1) + chains(i)%swap
2562 49 : buffer(spAll-3,2) = buffer(spAll-3,2) + DBLE(chains(i)%modGlobalMove(2))
2563 49 : buffer(spAll-4,1) = buffer(spAll-4,1) + a
2564 49 : buffer(spAll-4,2) = buffer(spAll-4,2) + b
2565 : !#ifdef CTCtqmc_CHECK
2566 49 : buffer(spAll-5,1) = buffer(spAll-5,1) + chains(i)%errorImpurity
2567 83 : buffer(spAll-5,2) = buffer(spAll-5,2) + chains(i)%errorBath
2568 : !#endif
2569 : end do
2570 :
2571 : #ifdef HAVE_MPI
2572 68 : CALL MPI_ALLREDUCE([this%runTime], rtime, 1, MPI_DOUBLE_PRECISION, MPI_MAX, this%MY_COMM, ierr)
2573 34 : this%runTime=rtime(1)
2574 : #if defined HAVE_MPI2_INPLACE
2575 : CALL MPI_ALLREDUCE(MPI_IN_PLACE, buffer, spAll*flavors, &
2576 34 : MPI_DOUBLE_PRECISION, MPI_SUM, this%MY_COMM, ierr)
2577 34 : IF ( this%opt_histo .GT. 0 ) THEN
2578 : CALL MPI_ALLREDUCE(MPI_IN_PLACE, this%occup_histo_time, flavors+1, MPI_DOUBLE_PRECISION, MPI_SUM, &
2579 3 : this%MY_COMM, ierr)
2580 : CALL MPI_ALLREDUCE(MPI_IN_PLACE, this%occupconfig, 2**flavors, MPI_DOUBLE_PRECISION, MPI_SUM, &
2581 3 : this%MY_COMM, ierr)
2582 : CALL MPI_ALLREDUCE(MPI_IN_PLACE, this%suscep, 3*this%samples, MPI_DOUBLE_PRECISION, MPI_SUM, &
2583 3 : this%MY_COMM, ierr)
2584 : CALL MPI_ALLREDUCE(MPI_IN_PLACE, this%chi, 3*this%samples, MPI_DOUBLE_PRECISION, MPI_SUM, &
2585 3 : this%MY_COMM, ierr)
2586 : CALL MPI_ALLREDUCE(MPI_IN_PLACE, this%chicharge, 3*this%samples, MPI_DOUBLE_PRECISION, MPI_SUM, &
2587 3 : this%MY_COMM, ierr)
2588 : CALL MPI_ALLREDUCE(MPI_IN_PLACE, this%ntot, 3, MPI_DOUBLE_PRECISION, MPI_SUM, &
2589 3 : this%MY_COMM, ierr)
2590 : END IF
2591 : #else
2592 : MALLOC(buffer2_out,(spAll,flavors))
2593 : CALL MPI_ALLREDUCE(buffer, buffer2_out, spAll*flavors, &
2594 : MPI_DOUBLE_PRECISION, MPI_SUM, this%MY_COMM, ierr)
2595 : buffer(1:spAll,1:flavors)=buffer2_out(1:spAll,1:flavors)
2596 : FREE(buffer2_out)
2597 : IF ( this%opt_histo .GT. 0 ) THEN
2598 : MALLOC(buffer1_out,(flavors+1))
2599 : CALL MPI_ALLREDUCE(this%occup_histo_time, buffer1_out, flavors+1, MPI_DOUBLE_PRECISION, MPI_SUM, &
2600 : this%MY_COMM, ierr)
2601 : this%occup_histo_time(1:flavors+1)=buffer1_out(1:flavors+1)
2602 : FREE(buffer1_out)
2603 : END IF
2604 : #endif
2605 : #endif
2606 :
2607 :
2608 113326 : buffer = buffer * inv_size
2609 4698 : this%measDE(:,:) = buffer(spall-(flavors+5):spAll-6,:)
2610 : ! this%seg_added = buffer(spAll ,1)
2611 : ! this%seg_removed = buffer(spAll-1,1)
2612 : ! this%seg_sign = buffer(spAll-2,1)
2613 : ! this%anti_added = buffer(spAll ,2)
2614 : ! this%anti_removed = buffer(spAll-1,2)
2615 : ! this%anti_sign = buffer(spAll-2,2)
2616 34 : this%stats(1) = buffer(spAll ,1)
2617 34 : this%stats(2) = buffer(spAll-1,1)
2618 34 : this%stats(3) = buffer(spAll-2,1)
2619 34 : this%stats(4) = buffer(spAll ,2)
2620 34 : this%stats(5) = buffer(spAll-1,2)
2621 34 : this%stats(6) = buffer(spAll-2,2)
2622 34 : this%swap = buffer(spAll-3,1)
2623 34 : this%modGlobalMove(2) = NINT(buffer(spAll-3,2))
2624 34 : a = buffer(spAll-4,1)
2625 34 : b = buffer(spAll-4,2)
2626 : !!#ifdef CTCtqmc_CHECK
2627 34 : this%errorImpurity= buffer(spAll-5,1)
2628 34 : this%errorBath = buffer(spAll-5,2)
2629 : !#endif
2630 :
2631 :
2632 402 : DO iflavor = 1, flavors
2633 78038 : this%Greens(iflavor)%oper = buffer(1:sp1 , iflavor)
2634 : END DO
2635 34 : last = sp1
2636 34 : IF ( this%opt_analysis .EQ. 1 ) THEN
2637 4374 : this%measCorrelation(:,1,:) = buffer(last+1:last+sp1,:)
2638 6 : last = last + sp1
2639 4374 : this%measCorrelation(:,2,:) = buffer(last+1:last+sp1,:)
2640 6 : last = last + sp1
2641 4374 : this%measCorrelation(:,3,:) = buffer(last+1:last+sp1,:)
2642 : last = last + sp1
2643 : END IF
2644 34 : IF ( this%opt_order .GT. 0 ) THEN
2645 16653 : this%measPerturbation(:,:) = buffer(last+1:last+this%opt_order, :)
2646 : last = last + this%opt_order
2647 : END IF
2648 34 : IF ( this%opt_noise .EQ. 1 ) THEN
2649 0 : alpha(:,:) = buffer(last+1:last+this%samples+1,:)
2650 0 : last = last + this%samples + 1
2651 0 : beta(:,:) = buffer(last+1:last+this%samples+1,:)
2652 402 : last = last + this%samples + 1
2653 : END IF
2654 : END IF
2655 402 : DO iflavor = 1, flavors
2656 : ! complete DE this
2657 2366 : this%measDE(iflavor, iflavor+1:flavors) = this%measDE(iflavor+1:flavors,iflavor)
2658 : END DO
2659 34 : FREE(buffer)
2660 :
2661 :
2662 34 : IF ( this%opt_spectra .GE. 1 ) THEN
2663 : !PROBLEM eos_gnu_13.2_openmpi . %endDensity was introduced throughout
2664 : ! endDensity = SIZE(this%density,2)
2665 0 : endDensity=this%endDensity
2666 : !ENDPROBLEM
2667 0 : IF ( this%density(1,endDensity) .EQ. -1.d0 ) &
2668 0 : endDensity = endDensity - 1
2669 0 : CALL FFTHyb_init(FFTmrka,endDensity,DBLE(this%thermalization)/DBLE(this%measurements*this%opt_spectra))
2670 : ! Not very Beauty
2671 0 : MALLOC(freqs,(1:FFTmrka%size/2))
2672 0 : DO iflavor = 1, flavors
2673 : ! mean value is removed to supress the continue composent
2674 0 : CALL FFTHyb_setData(FFTmrka,this%density(iflavor,1:endDensity)/this%beta+this%Greens(iflavor)%oper(this%samples+1))
2675 0 : CALL FFTHyb_run(FFTmrka,1)
2676 0 : CALL FFTHyb_getData(FFTmrka,endDensity,this%density(iflavor,:),freqs)
2677 : END DO
2678 0 : this%density(flavors+1,:) = -1.d0
2679 0 : this%density(flavors+1,1:FFTmrka%size/2) = freqs
2680 0 : CALL FFTHyb_destroy(FFTmrka)
2681 0 : FREE(freqs)
2682 : END IF
2683 :
2684 34 : this%a_Noise = a
2685 34 : this%b_Noise = b
2686 34 : IF ( this%opt_noise .EQ. 1 ) THEN
2687 0 : this%abNoiseG(1,:,:) = alpha
2688 0 : this%abNoiseG(2,:,:) = beta
2689 : END IF
2690 34 : FREE(alpha)
2691 34 : FREE(beta)
2692 :
2693 34 : IF ( this%opt_histo .GT. 0 ) THEN
2694 3 : write(this%ostream,*) "=== Histogram of occupations for complete simulation ===="
2695 : ! write(6,*) "sumh over procs", sumh
2696 3 : sumh=0
2697 48 : do n1=1,this%flavors+1
2698 45 : write(this%ostream,'(i4,f10.4)') n1-1, this%occup_histo_time(n1)/float(this%size*num_chains)
2699 48 : sumh=sumh+this%occup_histo_time(n1)/float(this%size*num_chains)
2700 : enddo
2701 3 : write(this%ostream,'(a,f10.4)') " all" , sumh
2702 3 : write(this%ostream,*) "================================="
2703 :
2704 :
2705 12 : MALLOC(occ,(2**this%flavors,1:flavors))
2706 : #ifdef FC_LLVM
2707 : !FIXME I don't understand why LLVM fails here...
2708 : !I put preproc so others know extra spaces are on purpose
2709 : MALLOC(occtot,(2**this%flavors) )
2710 : MALLOC(spintot,(2**this%flavors) )
2711 : #else
2712 9 : MALLOC(occtot,(2**this%flavors))
2713 9 : MALLOC(spintot,(2**this%flavors))
2714 : #endif
2715 49155 : do n1=1,2**this%flavors
2716 : ! Compute occupations of individual Orbitals
2717 49152 : n3=n1-1
2718 49152 : occtot(n1)=0
2719 49152 : spintot(n1)=0
2720 49152 : signe=1
2721 737280 : do n2=1,this%flavors
2722 688128 : remainder=modulo(n3,2)
2723 688128 : quotient=(n3-remainder)/2
2724 688128 : occ(n1,n2)=remainder
2725 688128 : n3=quotient
2726 688128 : occtot(n1)=occtot(n1)+occ(n1,n2)
2727 688128 : if(n2>=6) signe =-1
2728 : !if(n2>=7) signe =0
2729 737280 : spintot(n1)=spintot(n1)+occ(n1,n2)*signe
2730 : enddo
2731 49155 : this%occupconfig(n1)=this%occupconfig(n1)/float(this%size*num_chains)
2732 : enddo
2733 :
2734 3 : write(this%ostream,*) "=== Histogram of occupations of configurations for complete simulation ===="
2735 3 : sumh=0
2736 3 : if(this%flavors==14) then
2737 49155 : do n1=1,2**this%flavors
2738 49152 : write(this%ostream,'(i4,14i2,f20.2)') n1, (occ(n1,n2),n2=1,14),this%occupconfig(n1)
2739 49155 : sumh=sumh+this%occupconfig(n1)
2740 : enddo
2741 0 : else if (this%flavors==10) then
2742 0 : do n1=1,2**this%flavors
2743 0 : write(this%ostream,'(i4,10i2,f20.2)') n1, (occ(n1,n2),n2=1,10),this%occupconfig(n1)
2744 0 : sumh=sumh+this%occupconfig(n1)
2745 : enddo
2746 : end if
2747 3 : write(this%ostream,'(a,f10.4)') " all" , sumh
2748 :
2749 3 : sumtot=0
2750 36 : do nelec=0,10
2751 33 : spinmin=modulo(nelec,2)
2752 33 : if(nelec<=5) spinmax=nelec
2753 33 : if(nelec>=6) spinmax=10-nelec
2754 33 : dspin=1
2755 33 : write(this%ostream,*) "=== Histogram of occupations of configurations for total number of electrons",nelec
2756 36 : do spin=spinmin,spinmax,dspin
2757 93 : sumh=0
2758 1523805 : do n1=1,2**this%flavors
2759 1523805 : if(occtot(n1)==nelec.and.abs(spintot(n1))==spin) then
2760 36672 : sumh=sumh+this%occupconfig(n1)
2761 36672 : if(this%flavors==10) then
2762 0 : write(this%ostream,'(i8,10i2,a,i2,i3,f10.4)') n1,(occ(n1,n2),n2=1,this%flavors)," ",occtot(n1),spintot(n1),&
2763 0 : &this%occupconfig(n1)
2764 36672 : else if(this%flavors==14) then
2765 36672 : write(this%ostream,'(i8,14i2,a,i2,i3,f10.4)') n1,(occ(n1,n2),n2=1,this%flavors)," ",occtot(n1),spintot(n1),&
2766 73344 : &this%occupconfig(n1)
2767 : end if
2768 : endif
2769 : enddo
2770 93 : write(this%ostream,'(a,i4,a,i4,a,f10.4)') " === Sum of weights for",nelec," electrons and spin",spin," is ",sumh
2771 93 : sumtot=sumtot+sumh
2772 : enddo
2773 : enddo
2774 3 : write(this%ostream,'(a,f10.4)') "Full sum is",sumtot
2775 3 : FREE(occ)
2776 3 : FREE(occtot)
2777 3 : FREE(spintot)
2778 :
2779 : !==============================
2780 : ! Print Susceptibilities
2781 : !==============================
2782 3 : if(Iatom .lt. 10) then
2783 3 : write(atomnb, '("0",i1)') Iatom
2784 : else
2785 0 : write(atomnb, '(i2)') Iatom
2786 : end if
2787 : !Local Magnetic Susceptibility
2788 3 : if(this%opt_histo .gt. 1) then
2789 : !Scalar
2790 0 : if(this%nspinor .eq. 1) then
2791 0 : open(unit=735,file=trim(fname)//'_LocalSpinSuscept_atom_'//atomnb//'.dat',status='unknown',form='formatted')
2792 0 : write(735,*) '#Tau Total t2g eg'
2793 0 : do n1=1,this%samples
2794 0 : this%suscep(:,n1)=this%suscep(:,n1)/float(this%size*num_chains)/float(this%samples)
2795 0 : write(735,'(1x,f14.8,2x,f12.8,2x,f12.8,2x,f12.8)') (n1-1)*this%beta/this%samples,(this%suscep(n2,n1),n2=1,3)
2796 : enddo
2797 : !add tau=beta
2798 0 : write(735,'(1x,f14.8,2x,f12.8,2x,f12.8,2x,f12.8)') (this%samples)*this%beta/this%samples,(this%suscep(n2,1),n2=1,3)
2799 :
2800 : else
2801 : !SOC
2802 0 : open (unit=735,file=trim(fname)//'_LocalMagnSuscept_atom_'//atomnb//'.dat',status='unknown',form='formatted')
2803 0 : write(735,*) '#Tau Total Orbital Spin'
2804 0 : do n1=1,this%samples
2805 0 : this%chi(:,n1)=this%chi(:,n1)/float(this%size*num_chains)/float(this%samples)
2806 0 : write(735,'(1x,f14.8,2x,f12.8,2x,f12.8,2x,f12.8)') (n1-1)*this%beta/this%samples,(this%chi(n2,n1),n2=1,3)
2807 : enddo
2808 : !add tau=beta
2809 0 : write(735,'(1x,f14.8,2x,f12.8,2x,f12.8,2x,f12.8)') (this%samples)*this%beta/this%samples,(this%chi(n2,1),n2=1,3)
2810 : endif
2811 0 : close(unit=735)
2812 : endif
2813 :
2814 : !Local Charge Susceptibility
2815 3 : if(this%opt_histo .gt. 2) then
2816 0 : this%ntot(:)=this%ntot(:)/float(this%size*num_chains)/float(this%samples)
2817 0 : open (unit=735,file=trim(fname)//'_LocalChargeSuscept_atom_'//atomnb//'.dat',status='unknown',form='formatted')
2818 0 : write(735,*) '#Tau Total <ntot> '
2819 0 : do n1=1,this%samples
2820 0 : this%chicharge(1,n1)=(this%chicharge(1,n1)/float(this%size*num_chains)/float(this%samples))-(this%ntot(1)*this%ntot(1))
2821 : !this%chicharge(2,n1)=(this%chicharge(2,n1)/float(this%size*num_chains)/float(this%samples))-(this%ntot(2)*this%ntot(2))
2822 : !this%chicharge(3,n1)=(this%chicharge(3,n1)/float(this%size*num_chains)/float(this%samples))-(this%ntot(3)*this%ntot(3))
2823 : !write(735 '(1x,f14.8,2x,f14.8,2x,f14.8,2x,f14.8,2x,f14.8)') (n1-1)*this%beta/this%samples,(this%chicharge(n2,n1),n2=1,3),this%ntot(1)
2824 0 : write(735, '(1x,f14.8,2x,f14.8,2x,f14.8)') (n1-1)*this%beta/this%samples,(this%chicharge(1,n1)),this%ntot(1)
2825 : enddo
2826 : !add tau=beta
2827 0 : write(735, '(1x,f14.8,2x,f14.8,2x,f14.8)') (this%samples)*this%beta/this%samples,(this%chicharge(1,1)),this%ntot(1)
2828 0 : close(unit=735)
2829 : endif
2830 :
2831 : #if defined HAVE_FC_FLUSH
2832 3 : call flush(735)
2833 : #elif defined HAVE_FC_FLUSH_
2834 : call flush(735)
2835 : #endif
2836 :
2837 :
2838 :
2839 :
2840 : ENDIF
2841 :
2842 34 : END SUBROUTINE Ctqmc_getResult
2843 : !!***
2844 :
2845 : !!****f* ABINIT/m_Ctqmc/Ctqmc_symmetrizeGreen
2846 : !! NAME
2847 : !! Ctqmc_symmetrizeGreen
2848 : !!
2849 : !! FUNCTION
2850 : !! optionnaly symmetrize the green functions
2851 : !!
2852 : !! COPYRIGHT
2853 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
2854 : !! This file is distributed under the terms of the
2855 : !! GNU General Public License, see ~abinit/COPYING
2856 : !! or http://www.gnu.org/copyleft/gpl.txt .
2857 : !!
2858 : !! INPUTS
2859 : !! this=ctqmc
2860 : !! syms=weight factors
2861 : !!
2862 : !! OUTPUT
2863 : !! argout(sizeout)=description
2864 : !!
2865 : !! SIDE EFFECTS
2866 : !!
2867 : !! NOTES
2868 : !!
2869 : !! SOURCE
2870 :
2871 0 : SUBROUTINE Ctqmc_symmetrizeGreen(this, syms)
2872 :
2873 : !Arguments ------------------------------------
2874 : TYPE(Ctqmc) , INTENT(INOUT) :: this
2875 : DOUBLE PRECISION, DIMENSION(:,:), INTENT(IN ) :: syms
2876 : !Local variables ------------------------------
2877 : INTEGER :: iflavor1
2878 : INTEGER :: iflavor2
2879 : INTEGER :: flavors
2880 0 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: green_tmp
2881 0 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(: ) :: n_tmp
2882 :
2883 0 : flavors = this%flavors
2884 0 : IF ( SIZE(syms,1) .NE. flavors .OR. SIZE(syms,2) .NE. flavors ) THEN
2885 0 : CALL WARNALL("Ctqmc_symmetrizeGreen : wrong opt_sym -> not symmetrizing")
2886 : RETURN
2887 : END IF
2888 :
2889 0 : MALLOC(green_tmp,(1:this%samples+1,flavors))
2890 0 : green_tmp(:,:) = 0.d0
2891 0 : MALLOC(n_tmp,(1:flavors))
2892 0 : n_tmp(:) = 0.d0
2893 0 : DO iflavor1=1, flavors
2894 0 : DO iflavor2=1,flavors
2895 : green_tmp(:,iflavor1) = green_tmp(:,iflavor1) &
2896 0 : + syms(iflavor2,iflavor1) * this%Greens(iflavor2)%oper(:)
2897 : n_tmp(iflavor1) = n_tmp(iflavor1) &
2898 0 : + syms(iflavor2,iflavor1) * this%measN(1,iflavor2)
2899 : END DO
2900 : END DO
2901 0 : DO iflavor1=1, flavors
2902 0 : this%Greens(iflavor1)%oper(:) = green_tmp(:,iflavor1)
2903 0 : this%measN(1,iflavor1) = n_tmp(iflavor1)
2904 : END DO
2905 0 : FREE(green_tmp)
2906 0 : FREE(n_tmp)
2907 : END SUBROUTINE Ctqmc_symmetrizeGreen
2908 : !!***
2909 :
2910 : !!****f* ABINIT/m_Ctqmc/Ctqmc_getGreen
2911 : !! NAME
2912 : !! Ctqmc_getGreen
2913 : !!
2914 : !! FUNCTION
2915 : !! Get the full green functions in time and/or frequency
2916 : !!
2917 : !! COPYRIGHT
2918 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
2919 : !! This file is distributed under the terms of the
2920 : !! GNU General Public License, see ~abinit/COPYING
2921 : !! or http://www.gnu.org/copyleft/gpl.txt .
2922 : !!
2923 : !! INPUTS
2924 : !! this=ctqmc
2925 : !!
2926 : !! OUTPUT
2927 : !! Gtau=green function in time
2928 : !! Gw=green function in frequency
2929 : !!
2930 : !! SIDE EFFECTS
2931 : !!
2932 : !! NOTES
2933 : !!
2934 : !! SOURCE
2935 :
2936 34 : SUBROUTINE Ctqmc_getGreen(this, Gtau, Gw)
2937 :
2938 : !Arguments ------------------------------------
2939 : TYPE(Ctqmc) , INTENT(INOUT) :: this
2940 : DOUBLE PRECISION, DIMENSION(:,:), OPTIONAL, INTENT(INOUT) :: Gtau
2941 : COMPLEX(KIND=8), DIMENSION(:,:), OPTIONAL, INTENT(INOUT) :: Gw
2942 : !Local variables ------------------------------
2943 : !INTEGER :: itime
2944 : INTEGER :: iflavor1
2945 : INTEGER :: iflavor2
2946 : INTEGER :: iflavor3
2947 : INTEGER :: flavors
2948 : DOUBLE PRECISION :: u1
2949 : DOUBLE PRECISION :: u2
2950 : DOUBLE PRECISION :: Un
2951 : DOUBLE PRECISION :: UUnn
2952 :
2953 34 : flavors = this%flavors
2954 402 : DO iflavor1 = 1, flavors
2955 368 : u1 = 0.d0
2956 368 : u2 = 0.d0
2957 4664 : DO iflavor2 = 1, flavors
2958 4296 : IF ( iflavor2 .EQ. iflavor1 ) CYCLE
2959 3928 : Un = this%Impurity%mat_U(iflavor2,iflavor1) * this%measN(1,iflavor2)
2960 3928 : u1 = u1 + Un
2961 3928 : u2 = u2 + Un*this%Impurity%mat_U(iflavor2,iflavor1)
2962 52928 : DO iflavor3 = 1, flavors
2963 48632 : IF ( iflavor3 .EQ. iflavor2 .OR. iflavor3 .EQ. iflavor1 ) CYCLE
2964 40776 : UUnn = (this%Impurity%mat_U(iflavor2,iflavor1)*this%Impurity%mat_U(iflavor3,iflavor1)) * this%measDE(iflavor2,iflavor3)
2965 52928 : u2 = u2 + UUnn
2966 : END DO
2967 : END DO
2968 :
2969 368 : CALL GreenHyb_setMoments(this%Greens(iflavor1),u1,u2)
2970 368 : IF ( PRESENT( Gtau ) ) THEN
2971 77268 : Gtau(1:this%samples,iflavor1) = this%Greens(iflavor1)%oper(1:this%samples)
2972 : END IF
2973 : !write(6,*) "present gw", present(gw)
2974 402 : IF ( PRESENT( Gw ) ) THEN
2975 : !write(6,*) "size gw",SIZE(Gw,DIM=2) ,flavors+1
2976 368 : IF ( SIZE(Gw,DIM=2) .EQ. flavors+1 ) THEN
2977 368 : CALL GreenHyb_forFourier(this%Greens(iflavor1), Gomega=Gw(:,iflavor1), omega=Gw(:,this%flavors+1))
2978 : !IF ( this%rank .EQ. 0 ) write(20,*) Gw(:,iflavor1)
2979 0 : ELSE IF ( SIZE(Gw,DIM=2) .EQ. flavors ) THEN
2980 0 : CALL GreenHyb_forFourier(this%Greens(iflavor1),Gomega=Gw(:,iflavor1))
2981 : ELSE
2982 0 : CALL WARNALL("Ctqmc_getGreen : Gw is not valid ")
2983 0 : CALL GreenHyb_forFourier(this%Greens(iflavor1),Wmax=this%Wmax)
2984 : END IF
2985 : ELSE
2986 0 : CALL GreenHyb_forFourier(this%Greens(iflavor1),Wmax=this%Wmax)
2987 : END IF
2988 : END DO
2989 34 : END SUBROUTINE Ctqmc_getGreen
2990 : !!***
2991 :
2992 : !!****f* ABINIT/m_Ctqmc/Ctqmc_getD
2993 : !! NAME
2994 : !! Ctqmc_getD
2995 : !!
2996 : !! FUNCTION
2997 : !! get double occupation
2998 : !!
2999 : !! COPYRIGHT
3000 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
3001 : !! This file is distributed under the terms of the
3002 : !! GNU General Public License, see ~abinit/COPYING
3003 : !! or http://www.gnu.org/copyleft/gpl.txt .
3004 : !!
3005 : !! INPUTS
3006 : !! this=ctqmc
3007 : !!
3008 : !! OUTPUT
3009 : !! D=full double occupation
3010 : !!
3011 : !! SIDE EFFECTS
3012 : !!
3013 : !! NOTES
3014 : !!
3015 : !! SOURCE
3016 :
3017 34 : SUBROUTINE Ctqmc_getD(this, D)
3018 :
3019 : !Arguments ------------------------------------
3020 : TYPE(Ctqmc) , INTENT(IN ) :: this
3021 : DOUBLE PRECISION, DIMENSION(:,:), INTENT(OUT) :: D
3022 : !Local variables ------------------------------
3023 : INTEGER :: iflavor1
3024 : INTEGER :: iflavor2
3025 :
3026 34 : IF ( SIZE(D,1) .LT. this%flavors .OR. SIZE(D,2) .LT. this%flavors ) &
3027 0 : CALL ERROR("Ctqmc_getD : Dimensions of array D are too small")
3028 :
3029 4698 : D = 0.d0
3030 :
3031 402 : DO iflavor1 = 1, this%flavors
3032 4664 : DO iflavor2 = 1, this%flavors
3033 4664 : D(iflavor2,iflavor1) = this%measDE(iflavor2,iflavor1)
3034 : END DO
3035 402 : D(iflavor1,iflavor1) = 0.d0
3036 : END DO
3037 :
3038 34 : END SUBROUTINE Ctqmc_getD
3039 : !!***
3040 :
3041 : !!****f* ABINIT/m_Ctqmc/Ctqmc_getE
3042 : !! NAME
3043 : !! Ctqmc_getE
3044 : !!
3045 : !! FUNCTION
3046 : !! get interaction energy and noise on it
3047 : !!
3048 : !! COPYRIGHT
3049 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
3050 : !! This file is distributed under the terms of the
3051 : !! GNU General Public License, see ~abinit/COPYING
3052 : !! or http://www.gnu.org/copyleft/gpl.txt .
3053 : !!
3054 : !! INPUTS
3055 : !! this=ctqmc
3056 : !!
3057 : !! OUTPUT
3058 : !! E=interaction energy
3059 : !! noise=noise on this value
3060 : !!
3061 : !! SIDE EFFECTS
3062 : !!
3063 : !! NOTES
3064 : !!
3065 : !! SOURCE
3066 :
3067 34 : SUBROUTINE Ctqmc_getE(this,E,noise)
3068 :
3069 : !Arguments ------------------------------------
3070 : TYPE(Ctqmc) , INTENT(IN ) :: this
3071 : DOUBLE PRECISION, OPTIONAL, INTENT(OUT) :: E
3072 : DOUBLE PRECISION, OPTIONAL, INTENT(OUT) :: Noise
3073 :
3074 34 : IF ( PRESENT(E) ) &
3075 34 : E = this%measDE(1,1)
3076 34 : IF ( PRESENT(Noise) ) &
3077 : Noise = SUM(this%Impurity%mat_U)/(this%flavors*(this%flavors-1)) &
3078 0 : * this%a_Noise*(DBLE(this%sweeps)*DBLE(this%size))**this%b_Noise
3079 34 : END SUBROUTINE Ctqmc_getE
3080 : !!***
3081 :
3082 : !!****f* ABINIT/m_Ctqmc/Ctqmc_printAll
3083 : !! NAME
3084 : !! Ctqmc_printAll
3085 : !!
3086 : !! FUNCTION
3087 : !! print different functions computed during the simulation
3088 : !!
3089 : !! COPYRIGHT
3090 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
3091 : !! This file is distributed under the terms of the
3092 : !! GNU General Public License, see ~abinit/COPYING
3093 : !! or http://www.gnu.org/copyleft/gpl.txt .
3094 : !!
3095 : !! INPUTS
3096 : !! this=ctqmc
3097 : !!
3098 : !! OUTPUT
3099 : !!
3100 : !! SIDE EFFECTS
3101 : !!
3102 : !! NOTES
3103 : !!
3104 : !! SOURCE
3105 :
3106 34 : SUBROUTINE Ctqmc_printAll(this, num_chains)
3107 :
3108 : !Arguments ------------------------------------
3109 : TYPE(Ctqmc), INTENT(INOUT) :: this
3110 : INTEGER, INTENT(IN) :: num_chains
3111 :
3112 34 : IF ( .NOT. this%done ) &
3113 0 : CALL WARNALL("Ctqmc_printAll : Simulation not run ")
3114 :
3115 34 : CALL Ctqmc_printQMC(this, num_chains)
3116 :
3117 34 : CALL Ctqmc_printGreen(this)
3118 :
3119 34 : CALL Ctqmc_printD(this)
3120 :
3121 : ! CALL Ctqmc_printE(this)
3122 :
3123 : !#ifdef CTCtqmc_ANALYSIS
3124 34 : CALL Ctqmc_printPerturbation(this)
3125 :
3126 34 : CALL Ctqmc_printCorrelation(this)
3127 : !#endif
3128 :
3129 34 : CALL Ctqmc_printSpectra(this)
3130 : !ABI_BUG("toto")
3131 :
3132 34 : END SUBROUTINE Ctqmc_printAll
3133 : !!***
3134 :
3135 : !!****f* ABINIT/m_Ctqmc/Ctqmc_printQMC
3136 : !! NAME
3137 : !! Ctqmc_printQMC
3138 : !!
3139 : !! FUNCTION
3140 : !! print ctqmc statistics
3141 : !!
3142 : !! COPYRIGHT
3143 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
3144 : !! This file is distributed under the terms of the
3145 : !! GNU General Public License, see ~abinit/COPYING
3146 : !! or http://www.gnu.org/copyleft/gpl.txt .
3147 : !!
3148 : !! INPUTS
3149 : !! this=ctqmc
3150 : !!
3151 : !! OUTPUT
3152 : !!
3153 : !! SIDE EFFECTS
3154 : !!
3155 : !! NOTES
3156 : !!
3157 : !! SOURCE
3158 :
3159 34 : SUBROUTINE Ctqmc_printQMC(this, num_chains)
3160 :
3161 : !Arguments ------------------------------------
3162 : TYPE(Ctqmc), INTENT(INOUT) :: this
3163 : INTEGER, INTENT(IN) :: num_chains
3164 : !Local variables ------------------------------
3165 : INTEGER :: ostream
3166 : INTEGER :: iflavor
3167 : DOUBLE PRECISION :: sweeps
3168 : DOUBLE PRECISION :: invSweeps
3169 : CHARACTER(LEN=2) :: a
3170 : CHARACTER(LEN=15) :: string
3171 :
3172 : !IF ( this%rank .NE. 0) RETURN
3173 34 : IF ( this%rank .NE. MOD(this%size,this%size)) RETURN
3174 :
3175 19 : ostream = this%ostream
3176 19 : sweeps = DBLE(this%sweeps)
3177 19 : invSweeps = 1.d0/sweeps
3178 :
3179 19 : WRITE(ostream,'(1x,F13.0,A11,F10.2,A12,I5,A5)') sweeps*DBLE(this%size)*num_chains, " sweeps in ", this%runTime, &
3180 38 : " seconds on ", this%size*num_chains, " CPUs"
3181 19 : WRITE(ostream,'(A28,F6.2)') "Segments added [%] : ", this%stats(CTQMC_SEGME+CTQMC_ADDED)*invSweeps*100.d0
3182 19 : WRITE(ostream,'(A28,F6.2)') "Segments removed [%] : ", this%stats(CTQMC_SEGME+CTQMC_REMOV)*invSweeps*100.d0
3183 19 : WRITE(ostream,'(A28,F6.2)') "Segments sign [%] : ", this%stats(CTQMC_SEGME+CTQMC_DETSI)*invSweeps*100.d0
3184 19 : WRITE(ostream,'(A28,F6.2)') "Anti-segments added [%] : ", this%stats(CTQMC_ANTIS+CTQMC_ADDED)*invSweeps*100.d0
3185 19 : WRITE(ostream,'(A28,F6.2)') "Anti-segments removed [%] : ", this%stats(CTQMC_ANTIS+CTQMC_REMOV)*invSweeps*100.d0
3186 19 : WRITE(ostream,'(A28,F6.2)') "Anti-segments sign [%] : ", this%stats(CTQMC_ANTIS+CTQMC_DETSI)*invSweeps*100.d0
3187 19 : IF ( this%modGlobalMove(1) .LT. this%sweeps + 1 ) THEN
3188 15 : WRITE(ostream,'(A28,F6.2)') "Global Move [%] : ", this%swap *invSweeps*100.d0*this%modGlobalMove(1)
3189 15 : WRITE(ostream,'(A28,F6.2)') "Global Move Reduced [%] : ", this%swap / DBLE(this%modGlobalMove(2))*100.d0
3190 : END IF
3191 : !#ifdef CTCtqmc_CHECK
3192 19 : IF ( this%opt_check .EQ. 1 .OR. this%opt_check .EQ. 3 ) &
3193 0 : WRITE(ostream,'(A28,E22.14)') "Impurity test [%] : ", this%errorImpurity*100.d0
3194 19 : IF ( this%opt_check .GE. 2 ) &
3195 0 : WRITE(ostream,'(A28,E22.14)') "Bath test [%] : ", this%errorBath *100.d0
3196 : !#endif
3197 19 : WRITE(ostream,'(A28,ES22.14,A5,ES21.14)') "<Epot> [U] : ", this%measDE(1,1), " +/- ",&
3198 : !#ifdef HAVE_MPI
3199 2452 : SUM(DBLE(this%Impurity%mat_U))/(this%flavors*(this%flavors-1)) * this%a_Noise*(sweeps*DBLE(this%size*num_chains))**this%b_Noise
3200 : !#else
3201 : ! this%a_Noise*(sweeps)**this%b_Noise
3202 : !#endif
3203 19 : WRITE(ostream,'(A28,F8.4,A3,F7.4)') "Noise [/U] : ", this%a_Noise, " x^", this%b_Noise
3204 19 : WRITE(ostream,'(A28,E10.2)') "Niquist puls. [/beta] : ", ACOS(-1.d0)*this%inv_dt
3205 19 : WRITE(ostream,'(A28,E22.14)') "Max Acc. Epot Error [U] : ", this%measDE(2,2)/(this%beta*this%modNoise1*2.d0)*sweeps
3206 :
3207 : !WRITE(ostream,'(A28,F7.4,A3,F7.4,A4,E20.14)') "Noise [G(tau)] : ", this%a_Noise(2), "x^", this%b_Noise(2), " -> ", &
3208 : !this%a_Noise(2)*(sweeps*DBLE(this%size))**this%b_Noise(2)
3209 19 : IF ( this%opt_order .GT. 0 ) THEN
3210 12 : WRITE(a,'(I2)') this%flavors
3211 12 : string = '(A28,'//TRIM(ADJUSTL(a))//'(1x,I3))'
3212 12 : WRITE(ostream,string) "Perturbation orders : ", &
3213 8080 : (/ (MAXLOC(this%measPerturbation(:, iflavor))-1, iflavor=1, this%flavors) /)
3214 : END IF
3215 : !CALL FLUSH(this%ostream)
3216 : IF ( ABS(((this%stats(CTQMC_SEGME+CTQMC_ADDED) *invSweeps*100.d0) / &
3217 : (this%stats(CTQMC_SEGME+CTQMC_REMOV) *invSweeps*100.d0) - 1.d0)) .GE. 0.02d0 &
3218 19 : .OR. ABS(((this%stats(CTQMC_ANTIS+CTQMC_ADDED)*invSweeps*100.d0) / &
3219 : (this%stats(CTQMC_ANTIS+CTQMC_REMOV)*invSweeps*100.d0) - 1.d0)) .GE. 0.02d0 ) &
3220 : THEN
3221 0 : CALL WARNALL("Ctqmc_printQMC : bad statistic according to moves. Increase sweeps")
3222 : END IF
3223 : ! Check sign problem for diagonal hybridization.
3224 19 : IF ( (this%stats(CTQMC_SEGME+CTQMC_DETSI) + this%stats(CTQMC_ANTIS+CTQMC_DETSI)) .GT. 1.d-10 ) THEN
3225 1 : CALL WARNALL("Ctqmc_printQMC : at least one negative sign occured. There might be a bug in the CT-QMC")
3226 : END IF
3227 :
3228 19 : IF ( ABS(this%b_Noise+0.5)/0.5d0 .GE. 0.05d0 ) &
3229 8 : CALL WARNALL("Ctqmc_printQMC : bad statistic according to Noise. Increase sweeps")
3230 : ! IF ( ISNAN(this%a_Noise) .OR. ISNAN(this%a_Noise) ) &
3231 : ! CALL WARNALL("Ctqmc_printQMC : NaN appeared. Increase sweeps ")
3232 :
3233 :
3234 : END SUBROUTINE Ctqmc_printQMC
3235 : !!***
3236 :
3237 : !!****f* ABINIT/m_Ctqmc/Ctqmc_printGreen
3238 : !! NAME
3239 : !! Ctqmc_printGreen
3240 : !!
3241 : !! FUNCTION
3242 : !! print green functions
3243 : !!
3244 : !! COPYRIGHT
3245 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
3246 : !! This file is distributed under the terms of the
3247 : !! GNU General Public License, see ~abinit/COPYING
3248 : !! or http://www.gnu.org/copyleft/gpl.txt .
3249 : !!
3250 : !! INPUTS
3251 : !! this=ctqmc
3252 : !! oFileIn=file stream
3253 : !!
3254 : !! OUTPUT
3255 : !!
3256 : !! SIDE EFFECTS
3257 : !!
3258 : !! NOTES
3259 : !!
3260 : !! SOURCE
3261 :
3262 34 : SUBROUTINE Ctqmc_printGreen(this, oFileIn)
3263 :
3264 : !Arguments ------------------------------------
3265 : TYPE(Ctqmc) , INTENT(IN) :: this
3266 : INTEGER , OPTIONAL, INTENT(IN) :: oFileIn
3267 : !Local variables ------------------------------
3268 : INTEGER :: oFile
3269 : INTEGER :: itime
3270 : INTEGER :: sp1
3271 : INTEGER :: iflavor
3272 : INTEGER :: flavors
3273 : CHARACTER(LEN=4) :: cflavors
3274 : CHARACTER(LEN=50) :: string
3275 : DOUBLE PRECISION :: dt
3276 : DOUBLE PRECISION :: sweeps
3277 :
3278 : !IF ( this%rank .NE. MOD(1,this%size)) RETURN
3279 34 : IF ( this%rank .NE. MOD(this%size+1,this%size)) RETURN
3280 :
3281 19 : oFile = 40
3282 19 : IF ( PRESENT(oFileIn) ) THEN
3283 0 : oFile = oFileIn
3284 : ELSE
3285 19 : OPEN(UNIT=oFile, FILE="Gtau.dat")
3286 : END IF
3287 :
3288 19 : sp1 = this%samples
3289 19 : dt = this%beta / DBLE(sp1)
3290 19 : sp1 = sp1 + 1
3291 19 : flavors = this%flavors
3292 19 : sweeps = DBLE(this%sweeps)*DBLE(this%size)
3293 :
3294 19 : IF ( this%opt_noise .EQ. 1) THEN
3295 0 : WRITE(cflavors,'(I4)') 2*flavors+1
3296 0 : string = '(1x,'//TRIM(ADJUSTL(cflavors))//'ES22.14)'
3297 0 : DO itime = 1, sp1
3298 0 : WRITE(oFile,string) DBLE(itime-1)*dt, &
3299 0 : (/ (this%Greens(iflavor)%oper(itime), iflavor=1, flavors) /), &
3300 0 : (/ (this%abNoiseG(1,itime,iflavor)*(sweeps)**this%abNoiseG(2,itime,iflavor), iflavor=1, flavors) /)
3301 : END DO
3302 : ELSE
3303 19 : WRITE(cflavors,'(I4)') flavors+1
3304 19 : string = '(1x,'//TRIM(ADJUSTL(cflavors))//'ES22.14)'
3305 2788 : DO itime = 1, sp1
3306 2769 : WRITE(oFile,string) DBLE(itime-1)*dt, &
3307 61345 : (/ (this%Greens(iflavor)%oper(itime), iflavor=1, flavors) /)
3308 : END DO
3309 : END IF
3310 :
3311 19 : IF ( .NOT. PRESENT(oFileIn) ) CLOSE(oFile)
3312 :
3313 : END SUBROUTINE Ctqmc_printGreen
3314 : !!***
3315 :
3316 : !!****f* ABINIT/m_Ctqmc/Ctqmc_printD
3317 : !! NAME
3318 : !! Ctqmc_printD
3319 : !!
3320 : !! FUNCTION
3321 : !! print individual double occupancy
3322 : !!
3323 : !! COPYRIGHT
3324 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
3325 : !! This file is distributed under the terms of the
3326 : !! GNU General Public License, see ~abinit/COPYING
3327 : !! or http://www.gnu.org/copyleft/gpl.txt .
3328 : !!
3329 : !! INPUTS
3330 : !! this=ctqmc
3331 : !! oFileIn=file stream
3332 : !!
3333 : !! OUTPUT
3334 : !!
3335 : !! SIDE EFFECTS
3336 : !!
3337 : !! NOTES
3338 : !!
3339 : !! SOURCE
3340 :
3341 34 : SUBROUTINE Ctqmc_printD(this,oFileIn)
3342 :
3343 : !Arguments ------------------------------------
3344 : TYPE(Ctqmc) , INTENT(IN) :: this
3345 : INTEGER , OPTIONAL, INTENT(IN) :: oFileIn
3346 : !Local variables ------------------------------
3347 : INTEGER :: oFile
3348 : INTEGER :: iflavor1
3349 : INTEGER :: iflavor2
3350 :
3351 : !IF ( this%rank .NE. MOD(2,this%size)) RETURN
3352 34 : IF ( this%rank .NE. MOD(this%size+2,this%size)) RETURN
3353 :
3354 19 : oFile = 41
3355 19 : IF ( PRESENT(oFileIn) ) THEN
3356 0 : oFile = oFileIn
3357 : ELSE
3358 19 : OPEN(UNIT=oFile, FILE="D.dat")
3359 : END IF
3360 :
3361 213 : DO iflavor1 = 1, this%flavors
3362 1226 : DO iflavor2 = iflavor1+1, this%flavors
3363 1207 : WRITE(oFile,'(1x,A8,I4,A1,I4,A3,ES21.14)') "Orbitals", iflavor1, "-", iflavor2, " : ", this%measDE(iflavor2,iflavor1)
3364 : END DO
3365 : END DO
3366 :
3367 19 : IF ( .NOT. PRESENT(oFileIn) ) CLOSE(oFile)
3368 :
3369 : END SUBROUTINE Ctqmc_printD
3370 : !!***
3371 :
3372 : !!****f* ABINIT/m_Ctqmc/Ctqmc_printE
3373 : !! NAME
3374 : !! Ctqmc_printE
3375 : !!
3376 : !! FUNCTION
3377 : !! print energy and noise
3378 : !!
3379 : !! COPYRIGHT
3380 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
3381 : !! This file is distributed under the terms of the
3382 : !! GNU General Public License, see ~abinit/COPYING
3383 : !! or http://www.gnu.org/copyleft/gpl.txt .
3384 : !!
3385 : !! INPUTS
3386 : !! this=ctqmc
3387 : !! oFileIn=file stream
3388 : !!
3389 : !! OUTPUT
3390 : !!
3391 : !! SIDE EFFECTS
3392 : !!
3393 : !! NOTES
3394 : !!
3395 : !! SOURCE
3396 :
3397 0 : SUBROUTINE Ctqmc_printE(this,oFileIn)
3398 :
3399 : !Arguments ------------------------------------
3400 : TYPE(Ctqmc) , INTENT(IN) :: this
3401 : INTEGER , OPTIONAL, INTENT(IN) :: oFileIn
3402 : !Local variables ------------------------------
3403 : INTEGER :: oFile
3404 : DOUBLE PRECISION :: E
3405 : DOUBLE PRECISION :: Noise
3406 :
3407 : !IF ( this%rank .NE. MOD(3,this%size)) RETURN
3408 0 : IF ( this%rank .NE. MOD(this%size+3,this%size)) RETURN
3409 :
3410 0 : oFile = 42
3411 0 : IF ( PRESENT(oFileIn) ) THEN
3412 0 : oFile = oFileIn
3413 : ELSE
3414 0 : OPEN(UNIT=oFile, FILE="BetaENoise.dat")
3415 : END IF
3416 :
3417 0 : CALL Ctqmc_getE(this,E,Noise)
3418 :
3419 0 : WRITE(oFile,'(1x,F5.2,A2,ES21.14,A2,ES21.14)') this%beta, " ", E, " ", Noise
3420 :
3421 0 : IF ( .NOT. PRESENT(oFileIn) ) CLOSE(oFile)
3422 :
3423 : END SUBROUTINE Ctqmc_printE
3424 : !!***
3425 :
3426 : !#ifdef CTCtqmc_ANALYSIS
3427 : !!****f* ABINIT/m_Ctqmc/Ctqmc_printPerturbation
3428 : !! NAME
3429 : !! Ctqmc_printPerturbation
3430 : !!
3431 : !! FUNCTION
3432 : !! print perturbation order
3433 : !!
3434 : !! COPYRIGHT
3435 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
3436 : !! This file is distributed under the terms of the
3437 : !! GNU General Public License, see ~abinit/COPYING
3438 : !! or http://www.gnu.org/copyleft/gpl.txt .
3439 : !!
3440 : !! INPUTS
3441 : !! this=ctqmc
3442 : !! oFileIn=file stream
3443 : !!
3444 : !! OUTPUT
3445 : !! argout(sizeout)=description
3446 : !!
3447 : !! SIDE EFFECTS
3448 : !!
3449 : !! NOTES
3450 : !!
3451 : !! SOURCE
3452 :
3453 34 : SUBROUTINE Ctqmc_printPerturbation(this, oFileIn)
3454 :
3455 : !Arguments ------------------------------------
3456 : TYPE(Ctqmc) , INTENT(IN) :: this
3457 : INTEGER , OPTIONAL, INTENT(IN) :: oFileIn
3458 : !Local variables-------------------------------
3459 : INTEGER :: oFile
3460 : INTEGER :: iorder
3461 : INTEGER :: order
3462 : INTEGER :: iflavor
3463 : INTEGER :: flavors
3464 : CHARACTER(LEN=2) :: a
3465 : CHARACTER(LEN=50) :: string
3466 :
3467 : !IF ( this%rank .NE. MOD(4,this%size)) RETURN
3468 41 : IF ( this%rank .NE. MOD(this%size+4,this%size)) RETURN
3469 19 : IF ( this%opt_order .LE. 0 ) RETURN
3470 :
3471 12 : oFile = 43
3472 12 : IF ( PRESENT(oFileIn) ) THEN
3473 0 : oFile = oFileIn
3474 : ELSE
3475 12 : OPEN(UNIT=oFile, FILE="Perturbation.dat")
3476 : END IF
3477 :
3478 12 : order = this%opt_order
3479 12 : flavors = this%flavors
3480 :
3481 12 : WRITE(a,'(I2)') flavors
3482 12 : string = '(I5,'//TRIM(ADJUSTL(a))//'F19.15)'
3483 612 : DO iorder = 1, order
3484 600 : WRITE(oFile,string) iorder-1, &
3485 16412 : (/ (this%measPerturbation(iorder, iflavor), iflavor=1, flavors) /)
3486 : END DO
3487 :
3488 12 : IF ( .NOT. PRESENT(oFileIn) ) CLOSE(oFile)
3489 : END SUBROUTINE Ctqmc_printPerturbation
3490 : !!***
3491 :
3492 : !!****f* ABINIT/m_Ctqmc/Ctqmc_printCorrelation
3493 : !! NAME
3494 : !! Ctqmc_printCorrelation
3495 : !!
3496 : !! FUNCTION
3497 : !! print correlation fonctions
3498 : !!
3499 : !! COPYRIGHT
3500 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
3501 : !! This file is distributed under the terms of the
3502 : !! GNU General Public License, see ~abinit/COPYING
3503 : !! or http://www.gnu.org/copyleft/gpl.txt .
3504 : !!
3505 : !! INPUTS
3506 : !! this=ctqmc
3507 : !! oFileIn=file stream
3508 : !!
3509 : !! OUTPUT
3510 : !!
3511 : !! SIDE EFFECTS
3512 : !!
3513 : !! NOTES
3514 : !!
3515 : !! SOURCE
3516 :
3517 34 : SUBROUTINE Ctqmc_printCorrelation(this, oFileIn)
3518 :
3519 : !Arguments ------------------------------------
3520 : TYPE(Ctqmc) , INTENT(IN) :: this
3521 : INTEGER , OPTIONAL, INTENT(IN) :: oFileIn
3522 : !Local variables ------------------------------
3523 : INTEGER :: oFile
3524 : INTEGER :: itime
3525 : INTEGER :: sp1
3526 : INTEGER :: iflavor
3527 : INTEGER :: i
3528 : INTEGER :: flavors
3529 : CHARACTER(LEN=2) :: a
3530 : CHARACTER(LEN=50) :: string
3531 : DOUBLE PRECISION :: dt
3532 :
3533 : !IF ( this%rank .NE. MOD(5,this%size)) RETURN
3534 47 : IF ( this%rank .NE. MOD(this%size+5,this%size)) RETURN
3535 19 : IF ( this%opt_analysis .NE. 1 ) RETURN
3536 :
3537 6 : oFile = 44
3538 6 : IF ( PRESENT(oFileIn) ) THEN
3539 0 : oFile = oFileIn
3540 : ELSE
3541 6 : OPEN(UNIT=oFile, FILE="Correlation.dat")
3542 : END IF
3543 :
3544 6 : sp1 = this%samples
3545 6 : dt = this%beta / sp1
3546 6 : sp1 = sp1 + 1
3547 6 : flavors = this%flavors
3548 :
3549 6 : i = 3*flavors + 1
3550 6 : WRITE(a,'(I2)') i
3551 6 : WRITE(oFile,*) "# time (/ (segement, antiseg, correl), i=1, flavor/)"
3552 6 : string = '(1x,'//TRIM(ADJUSTL(a))//'F19.15)'
3553 312 : DO itime = 1, sp1
3554 306 : WRITE(oFile,string) DBLE(itime-1)*dt, &
3555 : (/ ( &
3556 : (/ ( this%measCorrelation(itime, i, iflavor), i=1,3) /) &
3557 17754 : , iflavor=1, flavors) /)
3558 : END DO
3559 :
3560 6 : IF ( .NOT. PRESENT(oFileIn) ) CLOSE(oFile)
3561 :
3562 : END SUBROUTINE Ctqmc_printCorrelation
3563 : !!***
3564 : !#endif
3565 :
3566 : !!****f* ABINIT/m_Ctqmc/Ctqmc_printSpectra
3567 : !! NAME
3568 : !! Ctqmc_printSpectra
3569 : !!
3570 : !! FUNCTION
3571 : !! print fourier transform of time evolution of number of electrons
3572 : !!
3573 : !! COPYRIGHT
3574 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
3575 : !! This file is distributed under the terms of the
3576 : !! GNU General Public License, see ~abinit/COPYING
3577 : !! or http://www.gnu.org/copyleft/gpl.txt .
3578 : !!
3579 : !! INPUTS
3580 : !! this=ctqmc
3581 : !! oFileIn=file stream
3582 : !!
3583 : !! OUTPUT
3584 : !!
3585 : !! SIDE EFFECTS
3586 : !!
3587 : !! NOTES
3588 : !!
3589 : !! SOURCE
3590 :
3591 34 : SUBROUTINE Ctqmc_printSpectra(this, oFileIn)
3592 :
3593 : !Arguments ------------------------------------
3594 : TYPE(Ctqmc) , INTENT(IN) :: this
3595 : INTEGER , OPTIONAL, INTENT(IN) :: oFileIn
3596 : !Local variables ------------------------------
3597 : INTEGER :: oFile
3598 : INTEGER :: flavors
3599 : INTEGER :: indDensity
3600 : INTEGER :: endDensity
3601 : CHARACTER(LEN=4) :: a
3602 : CHARACTER(LEN=16) :: formatSpectra
3603 :
3604 : !IF ( this%rank .NE. MOD(6,this%size)) RETURN
3605 34 : IF ( this%opt_spectra .LT. 1 ) RETURN
3606 :
3607 0 : oFile = 45+this%rank
3608 0 : a ="0000"
3609 0 : WRITE(a,'(I4)') this%rank
3610 0 : IF ( PRESENT(oFileIn) ) THEN
3611 0 : oFile = oFileIn
3612 : ELSE
3613 0 : OPEN(UNIT=oFile, FILE="Markov_"//TRIM(ADJUSTL(a))//".dat")
3614 : END IF
3615 :
3616 0 : flavors = this%flavors
3617 0 : WRITE(a,'(I4)') flavors+1
3618 0 : formatSpectra ='(1x,'//TRIM(ADJUSTL(a))//'ES22.14)'
3619 0 : WRITE(oFile,*) "# freq[/hermalization] FFT"
3620 :
3621 : !PROBLEM eos_gnu_13.2_openmpi . %endDensity was introduced throughout
3622 : ! endDensity = SIZE(this%density,2)
3623 0 : endDensity=this%endDensity
3624 : !ENDPROBLEM
3625 0 : DO WHILE ( this%density(flavors+1,endDensity) .EQ. -1 )
3626 0 : endDensity = endDensity -1
3627 : END DO
3628 :
3629 0 : DO indDensity = 1, endDensity
3630 0 : WRITE(oFile,formatSpectra) this%density(flavors+1,indDensity), this%density(1:flavors,indDensity)
3631 : END DO
3632 :
3633 0 : IF ( .NOT. PRESENT(oFileIn) ) CLOSE(oFile)
3634 :
3635 : END SUBROUTINE Ctqmc_printSpectra
3636 : !!***
3637 :
3638 : !!****f* ABINIT/m_Ctqmc/Ctqmc_destroy
3639 : !! NAME
3640 : !! Ctqmc_destroy
3641 : !!
3642 : !! FUNCTION
3643 : !! destroy and deallocate all variables
3644 : !!
3645 : !! COPYRIGHT
3646 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
3647 : !! This file is distributed under the terms of the
3648 : !! GNU General Public License, see ~abinit/COPYING
3649 : !! or http://www.gnu.org/copyleft/gpl.txt .
3650 : !!
3651 : !! INPUTS
3652 : !! this=ctqmc
3653 : !!
3654 : !! OUTPUT
3655 : !!
3656 : !! SIDE EFFECTS
3657 : !!
3658 : !! NOTES
3659 : !!
3660 : !! SOURCE
3661 :
3662 102 : SUBROUTINE Ctqmc_destroy(this)
3663 :
3664 : !Arguments ------------------------------------
3665 : TYPE(Ctqmc), INTENT(INOUT) :: this
3666 : !Local variables ------------------------------
3667 : INTEGER :: iflavor
3668 : INTEGER :: flavors
3669 : INTEGER :: i
3670 : INTEGER :: j
3671 : INTEGER :: k
3672 :
3673 102 : if ( this%init .EQV. .FALSE. ) RETURN
3674 :
3675 102 : flavors = this%flavors
3676 :
3677 102 : CALL ImpurityOperator_destroy(this%Impurity)
3678 102 : CALL BathOperator_destroy(this%Bath)
3679 102 : CALL Vector_destroy(this%measNoise(1))
3680 102 : CALL Vector_destroy(this%measNoise(2))
3681 :
3682 102 : IF ( ALLOCATED(this%Greens) ) THEN
3683 1178 : DO iflavor = 1, flavors
3684 1178 : CALL GreenHyb_destroy(this%Greens(iflavor))
3685 : END DO
3686 1178 : DT_FREE( this%Greens )
3687 : END IF
3688 : !#ifdef CTCtqmc_ANALYSIS
3689 102 : FREEIF(this%measCorrelation)
3690 102 : FREEIF(this%measPerturbation)
3691 102 : IF ( this%opt_histo .GT. 0 ) THEN
3692 3 : FREEIF(this%occup_histo_time)
3693 3 : FREEIF(this%occupconfig)
3694 3 : FREEIF(this%suscep)
3695 3 : FREEIF(this%chi)
3696 3 : FREEIF(this%chicharge)
3697 3 : FREEIF(this%ntot)
3698 : ENDIF
3699 102 : FREEIF(this%measN)
3700 102 : FREEIF(this%measDE)
3701 102 : FREEIF(this%mu)
3702 102 : FREEIF(this%abNoiseG)
3703 102 : IF ( ALLOCATED(this%measNoiseG) ) THEN
3704 0 : DO i=1,2
3705 0 : DO j = 1, this%flavors
3706 0 : DO k= 1, this%samples+1
3707 0 : CALL Vector_destroy(this%measNoiseG(k,j,i))
3708 : END DO
3709 : END DO
3710 : END DO
3711 0 : DT_FREE(this%measNoiseG)
3712 : END IF
3713 102 : FREEIF(this%density)
3714 102 : this%endDensity=-1
3715 : !#endif
3716 102 : this%ostream = 0
3717 102 : this%istream = 0
3718 :
3719 102 : this%sweeps = 0
3720 102 : this%thermalization = 0
3721 102 : this%flavors = 0
3722 102 : this%samples = 0
3723 102 : this%beta = 0.d0
3724 : ! this%seg_added = 0.d0
3725 : ! this%anti_added = 0.d0
3726 : ! this%seg_removed = 0.d0
3727 : ! this%anti_removed = 0.d0
3728 : ! this%seg_sign = 0.d0
3729 : ! this%anti_sign = 0.d0
3730 714 : this%stats = 0.d0
3731 102 : this%swap = 0.d0
3732 :
3733 :
3734 102 : this%set = .FALSE.
3735 102 : this%done = .FALSE.
3736 102 : this%init = .FALSE.
3737 : END SUBROUTINE Ctqmc_destroy
3738 : !!***
3739 :
3740 : !!****f* ABINIT/m_Ctqmcoffdiag/Ctqmc_setMagmom
3741 : !! NAME
3742 : !! Ctqmcoffdiag_setMagmom
3743 : !!
3744 : !! FUNCTION
3745 : !! set the Magnetic moment matrix for susceptibility
3746 : !!
3747 : !! COPYRIGHT
3748 : !! Copyright (C) 2013-2026 ABINIT group (F. Gendron)
3749 : !! This file is distributed under the terms of the
3750 : !! GNU General Public License, see ~abinit/COPYING
3751 : !! or http://www.gnu.org/copyleft/gpl.txt .
3752 : !!
3753 : !! INPUTS
3754 : !!
3755 : !! OUTPUT
3756 : !!
3757 : !! SIDE EFFECTS
3758 : !!
3759 : !! SOURCE
3760 :
3761 49 : SUBROUTINE Ctqmc_setMagmom(this,Magmom_orb,Magmom_spin,Magmom_tot)
3762 :
3763 : !Arguments ------------------------------------
3764 : TYPE(Ctqmc), INTENT(INOUT) :: this
3765 : DOUBLE PRECISION, DIMENSION(:,:), INTENT(IN) :: Magmom_orb
3766 : DOUBLE PRECISION, DIMENSION(:,:), INTENT(IN) :: Magmom_spin
3767 : DOUBLE PRECISION, DIMENSION(:,:), INTENT(IN) :: Magmom_tot
3768 : !Local variables ------------------------------
3769 : ! INTEGER :: iflavor1,iflavor2
3770 :
3771 : ! do iflavor1=1,10
3772 : ! do iflavor2=1,10
3773 : ! if(iflavor1==iflavor2) THEN
3774 : ! write(6,*) iflavor1, iflavor2, Magmom(iflavor1,iflavor2)
3775 : ! end if
3776 : ! end do
3777 : ! end do
3778 :
3779 147 : IF ( SIZE(Magmom_orb) .NE. this%flavors*this%flavors ) &
3780 0 : CALL ERROR("Ctqmc_setMagmomm : Wrong Magnetic Moment matrix (size) ")
3781 :
3782 49 : CALL ImpurityOperator_setMagmommat(this%Impurity, Magmom_orb, Magmom_spin, Magmom_tot)
3783 :
3784 49 : END SUBROUTINE Ctqmc_setMagmom
3785 : !!***
3786 :
3787 :
3788 0 : END MODULE m_Ctqmc
3789 : !!***
|
