Line data Source code
1 :
2 : #if defined HAVE_CONFIG_H
3 : #include "config.h"
4 : #endif
5 : !!****m* ABINIT/m_BathOperator
6 : !! NAME
7 : !! m_BathOperator
8 : !!
9 : !! FUNCTION
10 : !! Manage all stuff related to the bath for the
11 : !! simgle Anderson Impurity Model
12 : !!
13 : !! COPYRIGHT
14 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
15 : !! This file is distributed under the terms of the
16 : !! GNU General Public License, see ~abinit/COPYING
17 : !! or http://www.gnu.org/copyleft/gpl.txt .
18 : !!
19 : !! NOTES
20 : !!
21 : !! SOURCE
22 :
23 : #include "defs.h"
24 : MODULE m_BathOperator
25 : USE m_MatrixHyb
26 : USE m_Vector
27 : USE m_VectorInt
28 : USE m_Global
29 : USE m_ListCdagC
30 :
31 : IMPLICIT NONE
32 :
33 : !!***
34 :
35 : PRIVATE
36 :
37 : !!****t* m_BathOperator/BathOperator
38 : !! NAME
39 : !! BathOperator
40 : !!
41 : !! FUNCTION
42 : !! This structured datatype contains the necessary data
43 : !!
44 : !! COPYRIGHT
45 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
46 : !! This file is distributed under the terms of the
47 : !! GNU General Public License, see ~abinit/COPYING
48 : !! or http://www.gnu.org/copyleft/gpl.txt .
49 : !!
50 : !! SOURCE
51 :
52 : TYPE, PUBLIC :: BathOperator
53 : LOGICAL _PRIVATE :: set = .FALSE.
54 : LOGICAL :: MAddFlag = .FALSE. ! Set to true if we can compute a new M (see updateDetXX)
55 : LOGICAL :: MRemoveFlag = .FALSE. ! Set to true if we can compute a new M (see updateDetXX)
56 : LOGICAL _PRIVATE :: antiShift = .FALSE. ! shift when M is updated with antiseg
57 : LOGICAL _PRIVATE :: doCheck = .FALSE.
58 : INTEGER _PRIVATE :: flavors
59 : INTEGER :: activeFlavor
60 : INTEGER _PRIVATE :: samples
61 : INTEGER _PRIVATE :: sizeHybrid
62 : INTEGER _PRIVATE :: updatePosRow
63 : INTEGER _PRIVATE :: updatePosCol
64 : INTEGER _PRIVATE :: iTech
65 : INTEGER _PRIVATE :: checkNumber
66 : DOUBLE PRECISION _PRIVATE :: beta
67 : DOUBLE PRECISION _PRIVATE :: dt
68 : DOUBLE PRECISION _PRIVATE :: inv_dt
69 : DOUBLE PRECISION _PRIVATE :: meanError
70 : DOUBLE PRECISION _PRIVATE :: S
71 : DOUBLE PRECISION _PRIVATE :: Stau
72 : DOUBLE PRECISION _PRIVATE :: Stilde
73 : TYPE(Vector) _PRIVATE :: R
74 : TYPE(Vector) _PRIVATE :: Q
75 : TYPE(Vector) _PRIVATE :: Rtau
76 : TYPE(Vector) _PRIVATE :: Qtau
77 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) _PRIVATE :: F ! sample,Flavors
78 : TYPE(MatrixHyb) , ALLOCATABLE, DIMENSION(:) :: M ! Flavors
79 : TYPE(MatrixHyb) , ALLOCATABLE, DIMENSION(:) _PRIVATE :: M_update ! Flavors
80 : END TYPE BathOperator
81 : !!***
82 :
83 : PUBLIC :: BathOperator_init
84 : PUBLIC :: BathOperator_reset
85 : PUBLIC :: BathOperator_activateParticle
86 : PRIVATE :: BathOperator_hybrid
87 : PUBLIC :: BathOperator_getDetAdd
88 : PUBLIC :: BathOperator_getDetRemove
89 : PUBLIC :: BathOperator_getDetF
90 : PUBLIC :: BathOperator_setMAdd
91 : PUBLIC :: BathOperator_setMRemove
92 : PUBLIC :: BathOperator_swap
93 : PUBLIC :: BathOperator_initF
94 : PUBLIC :: BathOperator_setF
95 : PUBLIC :: BathOperator_printF
96 : PUBLIC :: BathOperator_printM
97 : PUBLIC :: BathOperator_destroy
98 : PUBLIC :: BathOperator_doCheck
99 : PRIVATE :: BathOperator_checkM
100 : PUBLIC :: BathOperator_getError
101 :
102 : CONTAINS
103 : !!***
104 :
105 : !!****f* ABINIT/m_BathOperator/BathOperator_init
106 : !! NAME
107 : !! BathOperator_init
108 : !!
109 : !! FUNCTION
110 : !! Initialize and allocate data
111 : !!
112 : !! COPYRIGHT
113 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
114 : !! This file is distributed under the terms of the
115 : !! GNU General Public License, see ~abinit/COPYING
116 : !! or http://www.gnu.org/copyleft/gpl.txt .
117 : !!
118 : !! INPUTS
119 : !! this=bath object
120 : !! flavors=numbers of flavors we have (including spin)
121 : !! samples=Time slices in the input file
122 : !! beta=inverse temperature
123 : !! iTech=imaginary time or frequencies
124 : !! It is imposes to imaginary time
125 : !!
126 : !! OUTPUT
127 : !!
128 : !! SIDE EFFECTS
129 : !!
130 : !! NOTES
131 : !!
132 : !! SOURCE
133 :
134 102 : SUBROUTINE BathOperator_init(this, flavors, samples, beta, iTech)
135 :
136 : !Arguments ------------------------------------
137 : TYPE(BathOperator), INTENT(INOUT) :: this
138 : INTEGER , INTENT(IN ) :: flavors
139 : INTEGER , INTENT(IN ) :: samples
140 : DOUBLE PRECISION , INTENT(IN ) :: beta
141 : !Local variables ------------------------------
142 : INTEGER , INTENT(IN ) :: iTech
143 : INTEGER :: it
144 :
145 102 : this%MAddFlag = .FALSE.
146 102 : this%MRemoveFlag = .FALSE.
147 102 : this%flavors = flavors
148 102 : this%beta = beta
149 102 : this%samples = samples
150 102 : this%sizeHybrid = samples + 1
151 102 : this%dt = beta / DBLE(samples)
152 102 : this%inv_dt = DBLE(samples) / beta
153 102 : this%activeFlavor= 0
154 102 : this%updatePosRow = 0
155 102 : this%updatePosCol = 0
156 102 : this%iTech = iTech
157 : !#ifdef CTQMC_CHECK
158 102 : this%checkNumber = 0
159 102 : this%meanError = 0.d0
160 102 : this%doCheck = .FALSE.
161 : !#endif
162 :
163 102 : FREEIF(this%F)
164 408 : MALLOC(this%F,(1:this%sizeHybrid+1,1:flavors))
165 102 : DT_FREEIF(this%M)
166 1382 : DT_MALLOC(this%M,(1:flavors))
167 102 : DT_FREEIF(this%M_update)
168 1280 : DT_MALLOC(this%M_update,(1:flavors))
169 :
170 102 : CALL Vector_init(this%R,100)
171 102 : CALL Vector_init(this%Q,100)
172 102 : CALL Vector_init(this%Rtau,100)
173 102 : CALL Vector_init(this%Qtau,100)
174 :
175 1178 : DO it = 1, flavors
176 1076 : CALL MatrixHyb_init(this%M(it),this%iTech,size=Global_SIZE,Wmax=samples) !FIXME Should be consistent with ListCagC
177 1178 : CALL MatrixHyb_init(this%M_update(it),this%iTech,size=Global_SIZE,Wmax=samples) !FIXME Should be consistent with ListCagC
178 : END DO
179 327130 : this%F = 0.d0
180 102 : this%set = .TRUE.
181 :
182 102 : END SUBROUTINE BathOperator_init
183 : !!***
184 :
185 : !!****f* ABINIT/m_BathOperator/BathOperator_reset
186 : !! NAME
187 : !! BathOperator_reset
188 : !!
189 : !! FUNCTION
190 : !! Reset all internal variables
191 : !!
192 : !! COPYRIGHT
193 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
194 : !! This file is distributed under the terms of the
195 : !! GNU General Public License, see ~abinit/COPYING
196 : !! or http://www.gnu.org/copyleft/gpl.txt .
197 : !!
198 : !! INPUTS
199 : !! this=bath operator to reset
200 : !!
201 : !! OUTPUT
202 : !!
203 : !! SIDE EFFECTS
204 : !!
205 : !! NOTES
206 : !!
207 : !! SOURCE
208 :
209 49 : SUBROUTINE BathOperator_reset(this)
210 :
211 : !Arguments ------------------------------------
212 : TYPE(BathOperator), INTENT(INOUT) :: this
213 : !Local variables ------------------------------
214 : INTEGER :: it
215 49 : this%MAddFlag = .FALSE.
216 49 : this%MRemoveFlag = .FALSE.
217 49 : this%activeFlavor = 0
218 49 : this%updatePosRow = 0
219 49 : this%updatePosCol = 0
220 : !#ifdef CTQMC_CHECK
221 49 : this%checkNumber = 0
222 49 : this%meanError = 0.d0
223 : !#endif
224 49 : this%doCheck = .FALSE.
225 49 : CALL Vector_clear(this%R)
226 49 : CALL Vector_clear(this%Q)
227 49 : CALL Vector_clear(this%Rtau)
228 49 : CALL Vector_clear(this%Qtau)
229 :
230 567 : DO it = 1, this%flavors
231 567 : CALL MatrixHyb_clear(this%M(it)) !FIXME Should be consistent with ListCagC
232 : END DO
233 153503 : this%F = 0.d0
234 :
235 49 : END SUBROUTINE BathOperator_reset
236 : !!***
237 :
238 : !!****f* ABINIT/m_BathOperator/BathOperator_activateParticle
239 : !! NAME
240 : !! BathOperator_activateParticle
241 : !!
242 : !! FUNCTION
243 : !! Just save on wicht flavor we are working
244 : !! It is better to use the macro defined in defs.h
245 : !!
246 : !! COPYRIGHT
247 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
248 : !! This file is distributed under the terms of the
249 : !! GNU General Public License, see ~abinit/COPYING
250 : !! or http://www.gnu.org/copyleft/gpl.txt .
251 : !!
252 : !! INPUTS
253 : !! this=bath operator
254 : !! flavor=the flavor to activate
255 : !!
256 : !! OUTPUT
257 : !!
258 : !! SIDE EFFECTS
259 : !!
260 : !! NOTES
261 : !!
262 : !! SOURCE
263 :
264 0 : SUBROUTINE BathOperator_activateParticle(this,flavor)
265 :
266 : !Arguments ------------------------------------
267 : TYPE(BathOperator), INTENT(INOUT) :: this
268 : !Local variables ------------------------------
269 : INTEGER , INTENT(IN ) :: flavor
270 :
271 0 : IF ( flavor .GT. this%flavors ) &
272 0 : CALL ERROR("BathOperator_activateParticle : out of range ")
273 0 : IF ( this%set .EQV. .TRUE. .AND. ALLOCATED(this%M) ) THEN
274 0 : this%activeFlavor = flavor
275 0 : this%MAddFlag = .FALSE.
276 0 : this%MRemoveFlag = .FALSE.
277 : ELSE
278 0 : CALL ERROR("BathOperator_activateParticle : not allocated ")
279 : END IF
280 0 : END SUBROUTINE BathOperator_activateParticle
281 : !!***
282 :
283 : !!****f* ABINIT/m_BathOperator/BathOperator_hybrid
284 : !! NAME
285 : !! BathOperator_hybrid
286 : !!
287 : !! FUNCTION
288 : !! Compute the hybridization for the active flavor
289 : !! at time time
290 : !!
291 : !! COPYRIGHT
292 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
293 : !! This file is distributed under the terms of the
294 : !! GNU General Public License, see ~abinit/COPYING
295 : !! or http://www.gnu.org/copyleft/gpl.txt .
296 : !!
297 : !! INPUTS
298 : !! this=bath operator
299 : !! time=time F(time)
300 : !!
301 : !! OUTPUT
302 : !!
303 : !! SIDE EFFECTS
304 : !!
305 : !! NOTES
306 : !!
307 : !! SOURCE
308 :
309 : DOUBLE PRECISION FUNCTION BathOperator_hybrid(this,time)
310 :
311 : TYPE(BathOperator), INTENT(IN) :: this
312 : DOUBLE PRECISION , INTENT(IN) :: time
313 : #include "BathOperator_hybrid.h"
314 :
315 : IF ( this%activeFlavor .LE. 0 ) &
316 : CALL ERROR("BathOperator_hybrid : no active hybrid func ")
317 : #include "BathOperator_hybrid"
318 : BathOperator_hybrid = hybrid
319 :
320 : END FUNCTION BathOperator_hybrid
321 : !!***
322 :
323 : !!****f* ABINIT/m_BathOperator/BathOperator_getDetAdd
324 : !! NAME
325 : !! BathOperator_getDetAdd
326 : !!
327 : !! FUNCTION
328 : !! Compute the determinant ratio when a (anti)segment
329 : !! is trying to be added and store some array for setMadd
330 : !!
331 : !! COPYRIGHT
332 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
333 : !! This file is distributed under the terms of the
334 : !! GNU General Public License, see ~abinit/COPYING
335 : !! or http://www.gnu.org/copyleft/gpl.txt .
336 : !!
337 : !! INPUTS
338 : !! this=bath operator
339 : !! CdagC_1=segment to be added
340 : !! position=ordered position of the Cdag time
341 : !! particle=full list of CdagC for activeFlavor
342 : !!
343 : !! OUTPUT
344 : !! BathOperator_getDetAdd=the det
345 : !!
346 : !! SIDE EFFECTS
347 : !!
348 : !! NOTES
349 : !!
350 : !! SOURCE
351 511203414 : DOUBLE PRECISION FUNCTION BathOperator_getDetAdd(this,CdagC_1, position, particle)
352 :
353 : !Arguments ------------------------------------
354 : TYPE(BathOperator) , INTENT(INOUT) :: this
355 : DOUBLE PRECISION, DIMENSION(1:2), INTENT(IN ) :: CdagC_1
356 : INTEGER , INTENT(IN ) :: position
357 : TYPE(ListCdagC), INTENT(IN ) :: particle
358 : !Local variables-------------------------------
359 : INTEGER :: it1
360 : INTEGER :: it2
361 : INTEGER :: it3
362 : INTEGER :: tail
363 : INTEGER :: new_tail
364 : DOUBLE PRECISION :: C
365 : DOUBLE PRECISION :: Cbeta
366 : DOUBLE PRECISION :: Cibeta
367 : DOUBLE PRECISION :: Cdag
368 : DOUBLE PRECISION :: Cdagbeta
369 : DOUBLE PRECISION :: beta
370 : DOUBLE PRECISION :: ratio
371 : DOUBLE PRECISION :: time
372 : ! TYPE(CdagC) , POINTER, DIMENSION(:) :: list => NULL()
373 : #include "BathOperator_hybrid.h"
374 :
375 511203414 : this%antiShift = .FALSE.
376 511203414 : beta = this%beta
377 511203414 : C = CdagC_1(C_)
378 : ! Cbeta = C.MOD.beta
379 511203414 : MODCYCLE(C,beta,Cbeta)
380 511203414 : Cdag = CdagC_1(Cdag_)
381 : ! cdagbeta = Cdag.MOD.beta
382 511203414 : MODCYCLE(Cdag,beta,CdagBeta)
383 : ! IF ( Cdag .GE. beta ) &
384 : ! CALL ERROR("BathOperator_getDetAdd : bad case ... ")
385 511203414 : IF ( this%activeFlavor .LE. 0 ) &
386 0 : CALL ERROR("BathOperator_getDetAdd : no active hybrid function ")
387 :
388 511203414 : tail = particle%tail
389 511203414 : new_tail = tail+1
390 : ! list => particle%list
391 :
392 : IF ( ((C .GT. Cdag) .AND. (position .EQ. -1)) &
393 511203414 : .OR. ((C .LT. Cdag) .AND. (tail .EQ. 0))) THEN ! Possible only if it is a segment
394 88942177 : this%updatePosRow = tail + 1
395 88942177 : this%updatePosCol = tail + 1
396 : ELSE
397 422261237 : this%updatePosRow = ABS(position)
398 422261237 : this%updatePosCol = ABS(position)
399 : END IF
400 :
401 : ! If antisegment, the det ratio has to be by -1 ( sign of the signature of one
402 : ! permutation line in the this
403 511203414 : IF ( C .LT. Cdag .AND. tail .GT. 0) THEN ! if antiseg
404 : ! ratio = -ratio
405 45290271 : this%updatePosRow = (this%updatePosRow + 1) !position in [1;tail]
406 45290271 : IF ( CdagBeta .LT. particle%list(this%updatePosCol,Cdag_) ) this%antiShift = .TRUE.
407 : END IF
408 :
409 : ! CALL Vector_setSize(this%R,tail)
410 : ! CALL Vector_setSize(this%Q,tail)
411 511203414 : Vector_QuickResize(this%R,new_tail)
412 511203414 : Vector_QuickResize(this%Q,new_tail)
413 511203414 : Vector_QuickResize(this%Rtau,new_tail)
414 511203414 : Vector_QuickResize(this%Qtau,new_tail)
415 :
416 1239633524 : DO it1 = 1, tail
417 728430110 : it2 = it1 + ( 1+SIGN(1,it1-this%updatePosRow) )/2
418 728430110 : it3 = it1 + ( 1+SIGN(1,it1-this%updatePoscol) )/2
419 :
420 728430110 : this%Rtau%vec(it2)= C - particle%list(it1,Cdag_)
421 : !this%Rtau%vec(it1)= C - particle%list(it1,Cdag_)
422 728430110 : time = Cbeta - particle%list(it1,Cdag_)
423 : #include "BathOperator_hybrid"
424 728430110 : this%R%vec(it1) = hybrid
425 : ! this%R%vec(it) = BathOperator_hybrid(this, Cbeta - list(it)%Cdag)
426 : ! Cibeta = list(it)%C.MOD.beta
427 728430110 : MODCYCLE(particle%list(it1,C_),beta,Cibeta)
428 728430110 : time = Cibeta - Cdagbeta
429 728430110 : this%Qtau%vec(it3)= time
430 : !this%Qtau%vec(it1)= time
431 : #include "BathOperator_hybrid"
432 1239633524 : this%Q%vec(it1) = hybrid
433 : !this%Q%vec(it3) = hybrid
434 : ! Q(it) = BathOperator_hybrid(this, Cibeta - Cdagbeta)
435 : END DO
436 : ! Compute S
437 511203414 : this%Stau = C - Cdagbeta
438 511203414 : this%Rtau%vec(this%updatePosRow) = this%Stau
439 511203414 : this%Qtau%vec(this%updatePosCol) = this%Rtau%vec(this%updatePosRow)
440 :
441 511203414 : time = Cbeta-Cdagbeta
442 : #include "BathOperator_hybrid"
443 511203414 : this%S = hybrid
444 :
445 : !ratio = this%S - DOT_PRODUCT(MATMUL(this%R%vec(1:tail),this%M(this%activeFlavor)%mat(1:tail,1:tail)),this%Q%vec(1:tail))
446 511203414 : ratio = 0.d0
447 1239633524 : DO it1 = 1, tail
448 : time = 0.d0
449 4091562908 : DO it2 = 1, tail
450 4091562908 : time = time + this%R%vec(it2) * this%M(this%activeFlavor)%mat(it2,it1)
451 : END DO
452 1239633524 : ratio = ratio + this%Q%vec(it1) * time
453 : END DO
454 511203414 : ratio = this%S - ratio
455 :
456 511203414 : this%Stilde = 1.d0 / ratio
457 :
458 : ! This IF is the LAST "NON CORRECTION" in my opinion this should not appears.
459 : ! IF ( MAX(C,Cdag) .GT. this%beta ) THEN
460 : ! WRITE(*,*) this%Stilde
461 : ! this%Stilde = - ABS(this%Stilde)
462 : ! END IF
463 :
464 : ! If antisegment, the det ratio has to be by -1 ( sign of the signature of one
465 : ! permutation line in the this)
466 511203414 : IF ( C .LT. Cdag .AND. tail .GT. 0) THEN ! if antiseg
467 45290271 : ratio = -ratio
468 : ENDIF
469 :
470 511203414 : BathOperator_getDetAdd = ratio
471 511203414 : this%MAddFlag = .TRUE.
472 : !#ifdef CTQMC_CHECK
473 : ! this%ListCdagC = particle
474 : !!write(*,*) this%Stilde
475 : !!write(*,*) this%antishift
476 : !!write(*,*) this%updatePosRow
477 : !!write(*,*) this%updatePosCol
478 : !#endif
479 :
480 511203414 : END FUNCTION BathOperator_getDetAdd
481 : !!***
482 :
483 : !!****f* ABINIT/m_BathOperator/BathOperator_getDetRemove
484 : !! NAME
485 : !! BathOperator_getDetRemove
486 : !!
487 : !! FUNCTION
488 : !! Compute the determinant ratio when a (anti)segment
489 : !! is trying to be removed
490 : !!
491 : !! COPYRIGHT
492 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
493 : !! This file is distributed under the terms of the
494 : !! GNU General Public License, see ~abinit/COPYING
495 : !! or http://www.gnu.org/copyleft/gpl.txt .
496 : !!
497 : !! INPUTS
498 : !! this=bath operator
499 : !! position=position of segment to be removed
500 : !!
501 : !! OUTPUT
502 : !! BathOperator_getDetRemove=the det
503 : !!
504 : !! SIDE EFFECTS
505 : !!
506 : !! NOTES
507 : !!
508 : !! SOURCE
509 :
510 573945019 : DOUBLE PRECISION FUNCTION BathOperator_getDetRemove(this,position)
511 :
512 : !Arguments ------------------------------------
513 : TYPE(BathOperator), INTENT(INOUT) :: this
514 : !Local arguments-------------------------------
515 : INTEGER , INTENT(IN ) :: position
516 : INTEGER :: ABSposition
517 : INTEGER :: tail
518 :
519 573945019 : IF ( this%activeFlavor .LE. 0 ) &
520 0 : CALL ERROR("BathOperator_getDetRemove : no active hybrid fun ")
521 :
522 573945019 : this%antiShift = .FALSE.
523 573945019 : tail = this%M(this%activeFlavor)%tail
524 573945019 : ABSposition = ABS(position)
525 573945019 : IF ( ABSposition .GT. tail ) &
526 0 : CALL ERROR("BathOperator_getDetRemove : position > M size ")
527 573945019 : this%updatePosCol = ABSposition
528 573945019 : this%antiShift = .FALSE.
529 573945019 : IF ( position .GT. 0 ) THEN
530 286970080 : this%updatePosRow = ABSposition
531 : ELSE
532 286974939 : this%updatePosRow = ABSposition+1
533 286974939 : IF ( ABSposition .EQ. tail ) THEN
534 182719745 : this%antiShift = .TRUE.
535 182719745 : this%updatePosRow = 1 !ABSposition - 1
536 : ! this%updatePosRow = ABSposition
537 : ! IF ( this%updatePosCol .EQ. 0) this%updatePosCol = tail
538 : END IF
539 : ENDIF
540 573945019 : this%Stilde = this%M(this%activeflavor)%mat(this%updatePosRow,this%updatePosCol)
541 573945019 : this%MRemoveFlag = .TRUE.
542 573945019 : BathOperator_getDetRemove = this%Stilde
543 :
544 : ! If remove an antiseg , the det ratio has to be multiplied by -1
545 573945019 : IF ( position .LT. 0 .AND. tail .GT. 1 ) &
546 160672199 : BathOperator_getDetRemove = - BathOperator_getDetRemove
547 : !#ifdef CTQMC_CHECK
548 : ! this%ListCdagC = particle
549 : !!write(*,*) this%updatePosRow, this%updatePosCol, position
550 : !!CALL ListCdagC_print(particle)
551 : !#endif
552 :
553 573945019 : END FUNCTION BathOperator_getDetRemove
554 : !!***
555 :
556 : !!****f* ABINIT/m_BathOperator/BathOperator_getDetF
557 : !! NAME
558 : !! BathOperator_getDetF
559 : !!
560 : !! FUNCTION
561 : !! Compute the determinant of the F this
562 : !! using the hybridization of flavor and the
563 : !! segments of particle
564 : !!
565 : !! COPYRIGHT
566 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
567 : !! This file is distributed under the terms of the
568 : !! GNU General Public License, see ~abinit/COPYING
569 : !! or http://www.gnu.org/copyleft/gpl.txt .
570 : !!
571 : !! INPUTS
572 : !! this=bath operator
573 : !! flavor=hybridization function to take
574 : !! particles=segments to use
575 : !!
576 : !! OUTPUT
577 : !! BathOperator_getDetF=the det
578 : !!
579 : !! SIDE EFFECTS
580 : !!
581 : !! NOTES
582 : !!
583 : !! SOURCE
584 :
585 457160 : DOUBLE PRECISION FUNCTION BathOperator_getDetF(this,flavor,particle)
586 :
587 : !Arguments ------------------------------------
588 : TYPE(BathOperator) , INTENT(INOUT) :: this
589 : INTEGER , INTENT(IN ) :: flavor
590 : TYPE(ListCdagC), OPTIONAL, INTENT(IN ) :: particle
591 : !Local arguments-------------------------------
592 : INTEGER :: iCdag
593 : INTEGER :: iC
594 : INTEGER :: tail
595 : DOUBLE PRECISION :: time
596 : DOUBLE PRECISION :: tC
597 : DOUBLE PRECISION :: tCdag
598 : DOUBLE PRECISION :: beta
599 : DOUBLE PRECISION :: mbeta_two
600 : DOUBLE PRECISION :: signe
601 : DOUBLE PRECISION :: inv_dt
602 : #include "BathOperator_hybrid.h"
603 :
604 457160 : BathOperator_getDetF = 1.d0 ! pour eviter des divisions par 0
605 457160 : IF ( PRESENT( particle ) ) THEN
606 228580 : tail = particle%tail
607 228580 : activeF = flavor
608 228580 : beta = this%beta
609 228580 : mbeta_two = -beta*0.5d0
610 228580 : inv_dt = this%inv_dt
611 228580 : CALL MatrixHyb_setSize(this%M_update(flavor),tail)
612 499199 : DO iCdag = 1, tail
613 270619 : tCdag = particle%list(iCdag,Cdag_)
614 1192942 : DO iC = 1, tail
615 : !tC = particle%list(C_,iC).MOD.beta
616 693743 : MODCYCLE(particle%list(iC,C_),beta,tC)
617 693743 : time = tC - tCdag
618 : #include "BathOperator_hybrid"
619 964362 : this%M_update(flavor)%mat(iC,iCdag) = hybrid
620 : END DO
621 : END DO
622 : ! mat_tau needs to be transpose of ordered time mat (way of measuring
623 : ! G(tau))
624 499199 : DO iC = 1, tail
625 270619 : tC = particle%list(iC,C_)
626 1192942 : DO iCdag = 1, tail
627 693743 : tCdag = particle%list(iCdag,Cdag_)
628 693743 : time = tC - tCdag
629 693743 : signe = SIGN(1.d0,time)
630 693743 : time = time + (signe-1.d0)*mbeta_two
631 964362 : this%M_update(flavor)%mat_tau(iCdag,iC) = INT( ( time * inv_dt ) + 1.5d0 )
632 : END DO
633 : END DO
634 228580 : CALL MatrixHyb_inverse(this%M_update(flavor),BathOperator_getDetF) ! calcul le det de la matrice et l'inverse
635 : ELSE
636 228580 : CALL MatrixHyb_getDet(this%M(flavor),BathOperator_getDetF) ! det M = 1/detF !
637 228580 : BathOperator_getDetF = 1.d0 / BathOperator_getDetF
638 : ENDIF
639 457160 : END FUNCTION BathOperator_getDetF
640 : !!***
641 :
642 : !!****f* ABINIT/m_BathOperator/BathOperator_setMAdd
643 : !! NAME
644 : !! BathOperator_setMAdd
645 : !!
646 : !! FUNCTION
647 : !! Update de M this inserting a row and a column
648 : !!
649 : !! COPYRIGHT
650 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
651 : !! This file is distributed under the terms of the
652 : !! GNU General Public License, see ~abinit/COPYING
653 : !! or http://www.gnu.org/copyleft/gpl.txt .
654 : !!
655 : !! INPUTS
656 : !! this=bath operator
657 : !! particle=segments of active flavor
658 : !!
659 : !! OUTPUT
660 : !!
661 : !! SIDE EFFECTS
662 : !!
663 : !! NOTES
664 : !!
665 : !! SOURCE
666 :
667 85920850 : SUBROUTINE BathOperator_setMAdd(this,particle)
668 :
669 : !Arguments ------------------------------------
670 : TYPE(BathOperator), INTENT(INOUT) :: this
671 : TYPE(ListCdagC) , INTENT(IN ) :: particle
672 : !Local variables ------------------------------
673 : INTEGER :: tail
674 : INTEGER :: new_tail
675 : INTEGER :: col
676 : INTEGER :: col_move
677 : INTEGER :: row_move
678 : INTEGER :: row
679 : INTEGER :: positionRow
680 : INTEGER :: positionCol
681 : INTEGER :: aF
682 : DOUBLE PRECISION :: Stilde
683 : DOUBLE PRECISION :: time
684 : DOUBLE PRECISION :: mbeta_two
685 : DOUBLE PRECISION :: inv_dt
686 85920850 : TYPE(Vector) :: vec_tmp
687 85920850 : TYPE(VectorInt) :: vecI_tmp
688 : INTEGER :: m
689 : INTEGER :: count
690 : INTEGER :: i
691 : INTEGER :: j
692 : INTEGER :: p
693 :
694 85920850 : IF ( this%MAddFlag .EQV. .FALSE. ) &
695 0 : CALL ERROR("BathOperator_setMAdd : MAddFlag turn off ")
696 85920850 : af = this%activeFlavor
697 85920850 : IF ( aF .LE. 0 ) &
698 0 : CALL ERROR("BathOperator_setMAdd : no active hybrid function ")
699 85920850 : tail = this%M(aF)%tail
700 85920850 : new_tail = tail + 1
701 : !CALL this_print(M)
702 :
703 85920850 : positionRow = this%updatePosRow
704 85920850 : positionCol = this%updatePosCol
705 85920850 : Stilde = this%Stilde
706 : ! write(6,*) "before", positionRow, positionCol
707 : !CALL MatrixHyb_print(this%M(aF),opt_print=1)
708 85920850 : CALL MatrixHyb_setSize(this%M(aF),new_tail)
709 :
710 : ! Compute Qtilde with Q
711 : !this%Q%vec(1:tail) = (-1.d0) * MATMUL(this%M(aF)%mat(1:tail,1:tail),this%Q%vec(1:tail)) * Stilde
712 2632197880 : this%Q%vec(1:tail) = MATMUL(this%M(aF)%mat(1:tail,1:tail),this%Q%vec(1:tail))
713 : !this%Q%vec(PositionRow:new_tail) = EOSHIFT(this%Q%vec(PositionRow:new_tail), SHIFT=-1, BOUNDARY=-1.d0, DIM=1)
714 : ! this%Qtau%vec(PositionCol:new_tail) = EOSHIFT(this%Qtau%vec(PositionCol:new_tail), SHIFT=-1, BOUNDARY=1.d0, DIM=1)
715 : ! this%Qtau%vec(PositionCol) = this%Stau
716 :
717 : !Compute Rtilde with R and without multiplying by Stilde
718 : !this%R%vec(1:tail) = (-1.d0) * MATMUL(this%R%vec(1:tail),this%M(aF)%mat(1:tail,1:tail))
719 2632197880 : this%R%vec(1:tail) = MATMUL(this%R%vec(1:tail),this%M(aF)%mat(1:tail,1:tail))
720 : !this%R%vec(PositionCol:new_tail) = EOSHIFT(this%R%vec(PositionCol:new_tail), SHIFT=-1, BOUNDARY=-1.d0, DIM=1)
721 : ! this%Rtau%vec(PositionRow:new_tail) = EOSHIFT(this%Rtau%vec(PositionRow:new_tail), SHIFT=-1, BOUNDARY=1.d0, DIM=1)
722 : ! this%Rtau%vec(PositionRow) = this%Stau
723 :
724 : !Compute the new M this
725 : !this%M(aF)%mat(PositionRow:new_tail,1:new_tail) = &
726 : ! EOSHIFT(this%M(aF)%mat(PositionRow:new_tail,1:new_tail),SHIFT=-1, BOUNDARY=0.d0, DIM=1)
727 : !this%M(aF)%mat(1:new_tail,PositionCol:new_tail) = &
728 : ! EOSHIFT(this%M(aF)%mat(1:new_tail,PositionCol:new_tail),SHIFT=-1, BOUNDARY=0.d0, DIM=2)
729 : ! ! this%M(aF)%mat(1:new_tail,1:new_tail) = this%M(aF)%mat(1:new_tail,1:new_tail) + &
730 : ! ! Stilde * MATMUL(RESHAPE(this%Q%vec(1:new_tail),(/ new_tail,1 /)),RESHAPE(this%R%vec(1:new_tail),(/ 1,new_tail /)))
731 :
732 : !this%M(aF)%mat_tau(PositionRow:new_tail,1:new_tail) = &
733 : ! EOSHIFT(this%M(aF)%mat_tau(PositionRow:new_tail,1:new_tail),SHIFT=-1, BOUNDARY=0, DIM=1)
734 : !this%M(aF)%mat_tau(1:new_tail,PositionCol:new_tail) = &
735 : ! EOSHIFT(this%M(aF)%mat_tau(1:new_tail,PositionCol:new_tail),SHIFT=-1, BOUNDARY=0, DIM=2)
736 :
737 85920850 : mbeta_two = -this%beta*0.5d0
738 85920850 : inv_dt = this%inv_dt
739 : !Shift mat_tau
740 : !update old m
741 360855297 : DO col=tail,1,-1
742 274934447 : col_move = col + ( 1+SIGN(1,col-PositionCol) )/2
743 2082328986 : DO row=tail,1,-1
744 1721473689 : row_move = row + ( 1+SIGN(1,row-PositionRow) )/2
745 1721473689 : this%M(aF)%mat_tau(row_move,col_move) = this%M(aF)%mat_tau(row,col)
746 1996408136 : this%M(aF)%mat(row_move,col_move) = this%M(aF)%mat(row,col) + this%Q%vec(row)*this%R%vec(col) * Stilde
747 : END DO
748 : END DO
749 : ! Add new stuff for new row
750 360855297 : DO row = 1, tail
751 274934447 : row_move = row + ( 1+SIGN(1,row-PositionRow) )/2
752 274934447 : this%M(aF)%mat(row_move,PositionCol) = -this%Q%vec(row)*Stilde
753 274934447 : time = this%Rtau%vec(row)
754 274934447 : time = time + ( SIGN(1.d0,time) - 1.d0 )*mbeta_two
755 360855297 : this%M(aF)%mat_tau(row,PositionCol) = INT ( (time*inv_dt) +1.5d0 )
756 : END DO
757 : ! Add last time missing in the loops
758 85920850 : time = this%Rtau%vec(new_tail)
759 85920850 : time = time + ( SIGN(1.d0,time) - 1.d0 )*mbeta_two
760 85920850 : this%M(aF)%mat_tau(new_tail,PositionCol) = INT ( (time*inv_dt) +1.5d0 )
761 : ! Add new stuff for new col
762 360855297 : DO col = 1, tail
763 274934447 : col_move = col + ( 1+SIGN(1,col-PositionCol) )/2
764 274934447 : this%M(aF)%mat(PositionRow,col_move) = -this%R%vec(col)*Stilde
765 274934447 : time = this%Qtau%vec(col)
766 274934447 : time = time + ( SIGN(1.d0,time) - 1.d0 )*mbeta_two
767 360855297 : this%M(aF)%mat_tau(PositionRow,col) = INT ( (time*inv_dt) +1.5d0 )
768 : END DO
769 : ! Add last time missing in the loops
770 85920850 : time = this%Qtau%vec(new_tail)
771 85920850 : time = time + ( SIGN(1.d0,time) - 1.d0 )*mbeta_two
772 85920850 : this%M(aF)%mat_tau(PositionRow,new_tail) = INT ( (time*inv_dt) +1.5d0 )
773 :
774 85920850 : this%M(aF)%mat(PositionRow,PositionCol) = Stilde
775 :
776 : !CALL MatrixHyb_print(this%M(aF),opt_print=1)
777 :
778 : ! DO col = 1, new_tail
779 : ! time = this%Rtau%vec(col)
780 : ! time = time + ( SIGN(1.d0,time) - 1.d0 )*mbeta_two
781 : ! this%M(aF)%mat_tau(col,PositionCol) = INT ( (time*inv_dt) +1.5d0 )
782 : ! time = this%Qtau%vec(col)
783 : ! time = time + ( SIGN(1.d0,time) - 1.d0 )*mbeta_two
784 : ! this%M(aF)%mat_tau(PositionRow,Col) = INT ( (time*inv_dt) +1.5d0 )
785 : ! time = this%R%vec(col)*Stilde
786 : ! DO row = 1, new_tail
787 : ! this%M(aF)%mat(row,col) = this%M(aF)%mat(row,col) + this%Q%vec(row)*time
788 : ! END DO
789 : ! END DO
790 :
791 : !col_move = new_tail
792 : !col = tail
793 : !DO col_move = new_tail, 1, -1
794 : ! IF ( col_move .EQ. positionCol ) THEN
795 : ! ! on calcule rajoute Q tilde
796 : ! !row_move = new_tail
797 : ! row = tail
798 : ! DO row_move = new_tail, 1, -1
799 : ! ! calcul itau
800 : ! IF ( row_move .EQ. positionRow ) THEN
801 : ! this%M(aF)%mat(row_move,col_move) = Stilde
802 : ! !time = this%Stau
803 : ! ELSE
804 : ! this%M(aF)%mat(row_move,col_move) = -this%Q%vec(row)*Stilde
805 : ! !time = this%Rtau%vec(row_move)
806 : ! row = row - 1
807 : ! END IF
808 : ! !time = time + ( SIGN(1.d0,time) - 1.d0 )*mbeta_two
809 : ! !this%M(aF)%mat_tau(row_move,col_move) = INT ( (time*inv_dt) +1.5d0 )
810 : ! END DO
811 : ! ! realignement des indices
812 : ! ELSE
813 : ! ! on calcule Ptilde
814 : ! !row_move = new_tail
815 : ! row = tail
816 : ! DO row_move = new_tail, 1, -1
817 : ! IF ( row_move .EQ. positionRow ) THEN
818 : ! this%M(aF)%mat(row_move,col_move) = -this%R%vec(col) * Stilde
819 : ! ! calcul itau
820 : ! !time = this%Qtau%vec(col_move)
821 : ! !time = time + ( SIGN(1.d0,time) - 1.d0 )*mbeta_two
822 : ! !this%M(aF)%mat_tau(row_move,col_move) = INT ( (time*inv_dt) +1.5d0 )
823 : ! ELSE
824 : ! this%M(aF)%mat(row_move,col_move) = this%M(aF)%mat(row,col) + this%Q%vec(row)*this%R%vec(col)*Stilde
825 : ! ! copy itau
826 : ! !this%M(aF)%mat_tau(row_move,col_move) = this%M(aF)%mat_tau(row,col)
827 : ! row = row - 1
828 : ! END IF
829 : ! END DO
830 : ! col = col - 1
831 : ! END IF
832 : !END DO
833 : ! write(6,*) "after"
834 : ! CALL MatrixHyb_print(this%M(aF),opt_print=1)
835 : !CALL this_inverse(M)
836 : !CALL MatrixHyb_print(M)
837 : !CALL this_inverse(M)
838 :
839 85920850 : IF ( this%antiShift .EQV. .TRUE. ) THEN ! antisegment
840 2859166 : CALL Vector_init(vec_tmp,new_tail)
841 2859166 : CALL VectorInt_init(vecI_tmp,new_tail)
842 : ! Shift if necessary according to this%antishift
843 : ! shift DIM=2 (col)
844 2859166 : p = new_tail - 1
845 2859166 : m = 1
846 2859166 : count = 0
847 5718332 : DO WHILE ( count .NE. new_tail )
848 16457108 : vec_tmp%vec(1:new_tail) = this%M(aF)%mat(1:new_tail,m)
849 16457108 : vecI_tmp%vec(1:new_tail) = this%M(aF)%mat_tau(1:new_tail,m)
850 2859166 : i = m
851 : !j = m+p
852 2859166 : MODCYCLE(m+p, new_tail, j)
853 13597942 : DO WHILE (j .NE. m)
854 83824528 : this%M(aF)%mat(1:new_tail,i) = this%M(aF)%mat(1:new_tail,j)
855 83824528 : this%M(aF)%mat_tau(1:new_tail,i) = this%M(aF)%mat_tau(1:new_tail,j)
856 10738776 : i = j
857 10738776 : MODCYCLE(j+p, new_tail, j)
858 10738776 : count = count+1
859 : END DO
860 16457108 : this%M(aF)%mat(1:new_tail,i) = vec_tmp%vec(1:new_tail)
861 16457108 : this%M(aF)%mat_tau(1:new_tail,i) = vecI_tmp%vec(1:new_tail)
862 2859166 : count = count+1
863 2859166 : m = m+1
864 : END DO
865 : ! shift DIM=1 (row)
866 5718332 : p = new_tail - 1
867 : m = 1
868 : count = 0
869 5718332 : DO WHILE ( count .NE. new_tail)
870 16457108 : vec_tmp%vec(1:new_tail) = this%M(aF)%mat(m,1:new_tail)
871 16457108 : vecI_tmp%vec(1:new_tail) = this%M(aF)%mat_tau(m,1:new_tail)
872 2859166 : i = m
873 : !j = m+p
874 2859166 : MODCYCLE(m+p, new_tail, j)
875 13597942 : DO WHILE ( j .NE. m )
876 83824528 : this%M(aF)%mat(i,1:new_tail) = this%M(aF)%mat(j,1:new_tail)
877 83824528 : this%M(aF)%mat_tau(i,1:new_tail) = this%M(aF)%mat_tau(j,1:new_tail)
878 10738776 : i = j
879 10738776 : MODCYCLE(j+p, new_tail, j)
880 10738776 : count = count+1
881 : END DO
882 16457108 : this%M(aF)%mat(i,1:new_tail) = vec_tmp%vec(1:new_tail)
883 16457108 : this%M(aF)%mat_tau(i,1:new_tail) = vecI_tmp%vec(1:new_tail)
884 2859166 : count = count+1
885 2859166 : m = m+1
886 : END DO
887 2859166 : CALL Vector_destroy(vec_tmp)
888 2859166 : CALL VectorInt_destroy(vecI_tmp)
889 : !this%M(aF)%mat(1:new_tail,1:new_tail) = CSHIFT(this%M(aF)%mat(1:new_tail,1:new_tail), SHIFT=-1, DIM=1) ! Shift to the bottom
890 : !this%M(aF)%mat(1:new_tail,1:new_tail) = CSHIFT(this%M(aF)%mat(1:new_tail,1:new_tail), SHIFT=-1, DIM=2) ! Shift to the right
891 : !this%M(aF)%mat_tau(1:new_tail,1:new_tail) = CSHIFT(this%M(aF)%mat_tau(1:new_tail,1:new_tail), SHIFT=-1, DIM=1) ! Shift to the bottom
892 : !this%M(aF)%mat_tau(1:new_tail,1:new_tail) = CSHIFT(this%M(aF)%mat_tau(1:new_tail,1:new_tail), SHIFT=-1, DIM=2) ! Shift to the right
893 : !CALL this_print(M)
894 : END IF
895 :
896 85920850 : IF ( this%doCheck .EQV. .TRUE.) THEN
897 : !#ifdef CTQMC_CHECK
898 0 : CALL BathOperator_checkM(this,particle)
899 : !#endif
900 : END IF
901 :
902 85920850 : this%MAddFlag = .FALSE.
903 :
904 85920850 : END SUBROUTINE BathOperator_setMAdd
905 : !!***
906 :
907 : !!****f* ABINIT/m_BathOperator/BathOperator_setMRemove
908 : !! NAME
909 : !! BathOperator_setMRemove
910 : !!
911 : !! FUNCTION
912 : !! delete one row and one column of the M this
913 : !!
914 : !! COPYRIGHT
915 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
916 : !! This file is distributed under the terms of the
917 : !! GNU General Public License, see ~abinit/COPYING
918 : !! or http://www.gnu.org/copyleft/gpl.txt .
919 : !!
920 : !! INPUTS
921 : !! this=bath operator
922 : !! particle=segments of the active flavor
923 : !!
924 : !! OUTPUT
925 : !!
926 : !! SIDE EFFECTS
927 : !!
928 : !! NOTES
929 : !!
930 : !! SOURCE
931 :
932 85918401 : SUBROUTINE BathOperator_setMRemove(this,particle)
933 :
934 : !Arguments ------------------------------------
935 : TYPE(BathOperator), INTENT(INOUT) :: this
936 : TYPE(ListCdagC) , INTENT(IN ) :: particle
937 : !Local variables ------------------------------
938 : INTEGER :: tail
939 : INTEGER :: new_tail
940 : INTEGER :: col
941 : INTEGER :: col_move
942 : INTEGER :: row_move
943 : INTEGER :: row
944 : INTEGER :: positionCol
945 : INTEGER :: positionRow
946 : INTEGER :: aF
947 : INTEGER :: m
948 : INTEGER :: count
949 : INTEGER :: i
950 : INTEGER :: j
951 : INTEGER :: p
952 : DOUBLE PRECISION :: invStilde
953 : DOUBLE PRECISION :: invStilde2
954 85918401 : TYPE(VectorInt) :: vecI_tmp
955 85918401 : TYPE(Vector) :: vec_tmp
956 :
957 85918401 : IF ( this%MRemoveFlag .EQV. .FALSE. ) &
958 0 : CALL ERROR("BathOperator_setMRemove : MRemoveFlag turn off ")
959 85918401 : af = this%activeFlavor
960 85918401 : IF ( aF .LE. 0 ) &
961 0 : CALL ERROR("BathOperator_setMRemove : no active hybrid func ")
962 85918401 : tail = this%M(aF)%tail
963 85918401 : new_tail = tail - 1
964 85918401 : positionCol = this%updatePosCol
965 85918401 : positionRow = this%updatePosRow
966 85918401 : invStilde = 1.d0 / this%Stilde
967 :
968 : ! write(6,*) "before", positionRow, positionCol
969 : ! CALL MatrixHyb_print(this%M(aF),opt_print=1)
970 :
971 : ! IF ( new_tail .EQ. 0 ) THEN
972 : !! IF ( this%antiShift .EQV. .TRUE. ) THEN
973 : !! this%M(aF)%mat(1,1) = 1.d0/BathOperator_Hybrid(this, this%beta)
974 : !! this%MRemoveFlag = .FALSE.
975 : !! RETURN
976 : !! END IF
977 : ! CALL MatrixHyb_clear(this%M(aF))
978 : ! this%MRemoveFlag = .FALSE.
979 : ! RETURN
980 : ! END IF
981 :
982 : ! CALL Vector_setSize(this%Q,new_tail)
983 : ! CALL Vector_setSize(this%R,new_tail)
984 85918401 : Vector_QuickResize(this%Q,new_tail)
985 85918401 : Vector_QuickResize(this%R,new_tail)
986 :
987 : ! We use R and Q as this%R%vec and this%Q%vec
988 : ! this%R%vec => this%R
989 : ! this%Q%vec => this%Q
990 :
991 : !row = 1
992 : !row_move = 1
993 : !col = 1
994 : !col_move = 1
995 360844149 : DO row_move = 1, new_tail
996 : !IF ( row .EQ. positionRow ) row = row + 1
997 : !IF ( col .EQ. positionCol ) col = col + 1
998 274925748 : col = row_move + (1+SIGN(1,row_move-positionCol))/2
999 274925748 : row = row_move + (1+SIGN(1,row_move-positionRow))/2
1000 274925748 : this%R%vec(row_move) = this%M(aF)%mat(positionRow,col)
1001 360844149 : this%Q%vec(row_move) = this%M(aF)%mat(row,positionCol)
1002 : !row = row + 1
1003 : !col = col + 1
1004 : END DO
1005 : !! this%R%vec(1:positionCol-1) = this%M(aF)%mat(positionRow,1:positionCol-1)
1006 : !! this%R%vec(positionCol:new_tail) = this%M(aF)%mat(positionRow,positionCol+1:tail)
1007 : !! this%Q%vec(1:positionRow-1) = this%M(aF)%mat(1:positionRow-1,positionCol)
1008 : !! this%Q%vec(positionRow:new_tail) = this%M(aF)%mat(positionRow+1:tail,positionCol)
1009 : !write(*,*) positionRow, positionCol
1010 : !CALL MatrixHyb_print(M)
1011 : !CALL Vector_print(this%R)
1012 : !CALL Vector_print(this%Q)
1013 : !CALL ListCdagC_print(this%ListCdagC)
1014 :
1015 : !col = 1
1016 360844149 : DO col_move = 1, new_tail
1017 : !IF ( col_move .EQ. positionCol ) col = col + 1
1018 274925748 : col = col_move + (1+SIGN(1,col_move-positionCol))/2
1019 : !row = 1
1020 274925748 : invStilde2 = invStilde * this%R%vec(col_move)
1021 2082266517 : DO row_move = 1, new_tail
1022 : !IF ( row_move .EQ. positionRow ) row = row + 1
1023 1721422368 : row = row_move + (1+SIGN(1,row_move-positionRow))/2
1024 : this%M(aF)%mat(row_move,col_move) = this%M(aF)%mat(row,col) &
1025 1721422368 : - this%Q%vec(row_move)*invStilde2
1026 1996348116 : this%M(aF)%mat_tau(row_move,col_move) = this%M(aF)%mat_tau(row,col)
1027 : !row = row + 1
1028 : END DO
1029 : !col = col + 1
1030 : END DO
1031 85918401 : CALL MatrixHyb_setSize(this%M(aF),new_tail)
1032 :
1033 85918401 : IF ( this%antiShift .EQV. .TRUE. ) THEN ! antisegment
1034 : ! Shift if necessary according to this%antishift
1035 : ! shift DIM=2 (col)
1036 2896336 : CALL Vector_init(vec_tmp,new_tail)
1037 2896336 : CALL VectorInt_init(vecI_tmp,new_tail)
1038 2896336 : p = 1
1039 2896336 : m = 1
1040 2896336 : count = 0
1041 5756133 : DO WHILE ( count .NE. new_tail )
1042 13597195 : vec_tmp%vec(1:new_tail) = this%M(aF)%mat(1:new_tail,m)
1043 13597195 : vecI_tmp%vec(1:new_tail) = this%M(aF)%mat_tau(1:new_tail,m)
1044 2859797 : i = m
1045 : !j = m+p
1046 2859797 : MODCYCLE(m+p, new_tail, j)
1047 10737398 : DO WHILE (j .NE. m)
1048 59469081 : this%M(aF)%mat(1:new_tail,i) = this%M(aF)%mat(1:new_tail,j)
1049 59469081 : this%M(aF)%mat_tau(1:new_tail,i) = this%M(aF)%mat_tau(1:new_tail,j)
1050 7877601 : i = j
1051 7877601 : MODCYCLE(j+p, new_tail, j)
1052 7877601 : count = count+1
1053 : END DO
1054 13597195 : this%M(aF)%mat(1:new_tail,i) = vec_tmp%vec(1:new_tail)
1055 13597195 : this%M(aF)%mat_tau(1:new_tail,i) = vecI_tmp%vec(1:new_tail)
1056 2859797 : count = count+1
1057 2859797 : m = m+1
1058 : END DO
1059 2896336 : CALL Vector_destroy(vec_tmp)
1060 2896336 : CALL VectorInt_destroy(vecI_tmp)
1061 : !this%M(aF)%mat(1:new_tail,1:new_tail) = &
1062 : ! CSHIFT(this%M(aF)%mat(1:new_tail,1:new_tail), SHIFT=1, DIM=2) ! Shift to the top
1063 : !this%M(aF)%mat_tau(1:new_tail,1:new_tail) = &
1064 : ! CSHIFT(this%M(aF)%mat_tau(1:new_tail,1:new_tail), SHIFT=1, DIM=2) ! Shift to the top
1065 : END IF
1066 : ! write(6,*) "after "
1067 : ! CALL MatrixHyb_print(this%M(aF),opt_print=1)
1068 :
1069 85918401 : IF ( this%doCheck .EQV. .TRUE. ) THEN
1070 : !#ifdef CTQMC_CHECK
1071 0 : CALL BathOperator_checkM(this,particle)
1072 : !#endif
1073 : END IF
1074 :
1075 85918401 : this%MRemoveFlag = .FALSE.
1076 :
1077 85918401 : END SUBROUTINE BathOperator_setMRemove
1078 : !!***
1079 :
1080 : !!****f* ABINIT/m_BathOperator/BathOperator_swap
1081 : !! NAME
1082 : !! BathOperator_swap
1083 : !!
1084 : !! FUNCTION
1085 : !! Recompute 2 M this swaping the segments
1086 : !!
1087 : !! COPYRIGHT
1088 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
1089 : !! This file is distributed under the terms of the
1090 : !! GNU General Public License, see ~abinit/COPYING
1091 : !! or http://www.gnu.org/copyleft/gpl.txt .
1092 : !!
1093 : !! INPUTS
1094 : !! this=bath operator
1095 : !! iflavor1=flavor to swap with the next one
1096 : !! iflavor2=favor to swap with the previous one
1097 : !!
1098 : !! OUTPUT
1099 : !!
1100 : !! SIDE EFFECTS
1101 : !!
1102 : !! NOTES
1103 : !!
1104 : !! SOURCE
1105 :
1106 95797 : SUBROUTINE BathOperator_swap(this, flavor1, flavor2)
1107 :
1108 : !Arguments ------------------------------------
1109 : TYPE(BathOperator), INTENT(INOUT) :: this
1110 : INTEGER , INTENT(IN ) :: flavor1
1111 : INTEGER , INTENT(IN ) :: flavor2
1112 :
1113 : !CALL MatrixHyb_print(this%M(flavor1),234)
1114 95797 : this%M(flavor1) = this%M_update(flavor1)
1115 : !CALL MatrixHyb_print(this%M(flavor1),234)
1116 : !CALL MatrixHyb_print(this%M(flavor2),234)
1117 95797 : this%M(flavor2) = this%M_update(flavor2)
1118 : !CALL MatrixHyb_print(this%M(flavor2),234)
1119 :
1120 95797 : END SUBROUTINE BathOperator_swap
1121 : !!***
1122 :
1123 : !!****f* ABINIT/m_BathOperator/BathOperator_initF
1124 : !! NAME
1125 : !! BathOperator_initF
1126 : !!
1127 : !! FUNCTION
1128 : !! Copy input hybridization functions from a file
1129 : !!
1130 : !! COPYRIGHT
1131 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
1132 : !! This file is distributed under the terms of the
1133 : !! GNU General Public License, see ~abinit/COPYING
1134 : !! or http://www.gnu.org/copyleft/gpl.txt .
1135 : !!
1136 : !! INPUTS
1137 : !! this=bath operator
1138 : !! ifstream=file stream to read F
1139 : !!
1140 : !! OUTPUT
1141 : !! argout(sizeout)=description
1142 : !!
1143 : !! SIDE EFFECTS
1144 : !!
1145 : !! NOTES
1146 : !!
1147 : !! SOURCE
1148 :
1149 0 : SUBROUTINE BathOperator_initF(this,ifstream)
1150 :
1151 : !Arguments ----------------------
1152 : TYPE(BathOperator), INTENT(INOUT) :: this
1153 : INTEGER , INTENT(IN ) :: ifstream
1154 : !Local variables ----------------
1155 : INTEGER :: flavor
1156 : INTEGER :: sample
1157 :
1158 0 : IF ( this%set .EQV. .FALSE. ) &
1159 0 : CALL ERROR("BathOperator_initF : BathOperator not set ")
1160 :
1161 0 : DO flavor=1,this%flavors
1162 0 : DO sample = 1, this%sizeHybrid
1163 0 : READ(ifstream,*) this%F(sample,flavor)
1164 : END DO
1165 : END DO
1166 0 : END SUBROUTINE BathOperator_initF
1167 : !!***
1168 :
1169 : !!****f* ABINIT/m_BathOperator/BathOperator_setF
1170 : !! NAME
1171 : !! BathOperator_setF
1172 : !!
1173 : !! FUNCTION
1174 : !! Copy F from input array
1175 : !!
1176 : !! COPYRIGHT
1177 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
1178 : !! This file is distributed under the terms of the
1179 : !! GNU General Public License, see ~abinit/COPYING
1180 : !! or http://www.gnu.org/copyleft/gpl.txt .
1181 : !!
1182 : !! INPUTS
1183 : !! this=bath operator
1184 : !! F=array of the hybridization function
1185 : !!
1186 : !! OUTPUT
1187 : !!
1188 : !! SIDE EFFECTS
1189 : !!
1190 : !! NOTES
1191 : !!
1192 : !! SOURCE
1193 :
1194 49 : SUBROUTINE BathOperator_setF(this,F)
1195 :
1196 : !Arguments ------------------------------------
1197 : TYPE(BathOperator) , INTENT(INOUT) :: this
1198 : DOUBLE PRECISION, DIMENSION(:,:) , INTENT(IN ) :: F
1199 : !Arguments ------------------------------------
1200 : INTEGER :: flavor
1201 : INTEGER :: sample
1202 : INTEGER :: length
1203 :
1204 49 : IF ( this%set .EQV. .FALSE. ) &
1205 0 : CALL ERROR("BathOperator_setF : BathOperator not set ")
1206 :
1207 147 : length = SIZE(F)
1208 49 : IF ( length .NE. (this%flavors * this%sizeHybrid) ) &
1209 0 : CALL ERROR("BathOperator_setF : wrong input F ")
1210 :
1211 567 : DO flavor=1,this%flavors
1212 152985 : DO sample = 1, this%sizeHybrid
1213 152936 : this%F(sample,flavor) = F(sample,flavor)
1214 : END DO
1215 : END DO
1216 49 : END SUBROUTINE BathOperator_setF
1217 : !!***
1218 :
1219 : !!****f* ABINIT/m_BathOperator/BathOperator_printF
1220 : !! NAME
1221 : !! BathOperator_printF
1222 : !!
1223 : !! FUNCTION
1224 : !! print F function
1225 : !!
1226 : !! COPYRIGHT
1227 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
1228 : !! This file is distributed under the terms of the
1229 : !! GNU General Public License, see ~abinit/COPYING
1230 : !! or http://www.gnu.org/copyleft/gpl.txt .
1231 : !!
1232 : !! INPUTS
1233 : !! this=bath operator
1234 : !! ostream=file stream to write in
1235 : !!
1236 : !! OUTPUT
1237 : !!
1238 : !! SIDE EFFECTS
1239 : !!
1240 : !! NOTES
1241 : !!
1242 : !! SOURCE
1243 :
1244 0 : SUBROUTINE BathOperator_printF(this,ostream)
1245 :
1246 : !Arguments ------------------------------------
1247 : TYPE(BathOperator), INTENT(INOUT) :: this
1248 : INTEGER,OPTIONAL , INTENT(IN ) :: ostream
1249 : !Local variables ------------------------------
1250 : CHARACTER(LEN=4) :: aflavor
1251 : CHARACTER(LEN=50) :: string
1252 : INTEGER :: flavor
1253 : INTEGER :: sample
1254 : INTEGER :: ostream_val
1255 :
1256 0 : IF ( PRESENT(ostream) ) THEN
1257 0 : ostream_val = ostream
1258 : ELSE
1259 0 : ostream_val = 65
1260 0 : OPEN(UNIT=ostream_val, FILE="F.dat")
1261 : END IF
1262 :
1263 0 : WRITE(aflavor,'(I4)') this%flavors+1
1264 0 : string = '(1x,'//TRIM(ADJUSTL(aflavor))//'E22.14)'
1265 0 : DO sample = 1, this%sizeHybrid
1266 0 : WRITE(ostream_val,string) (sample-1)*this%dt, (this%F(sample,flavor), flavor=1,this%flavors)
1267 : END DO
1268 : !CALL FLUSH(ostream_val)
1269 :
1270 0 : IF ( .NOT. PRESENT(ostream) ) &
1271 0 : CLOSE(ostream_val)
1272 :
1273 0 : END SUBROUTINE BathOperator_printF
1274 : !!***
1275 :
1276 : !!****f* ABINIT/m_BathOperator/BathOperator_printM
1277 : !! NAME
1278 : !! BathOperator_printM
1279 : !!
1280 : !! FUNCTION
1281 : !! print M =F^{-1} this
1282 : !!
1283 : !! COPYRIGHT
1284 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
1285 : !! This file is distributed under the terms of the
1286 : !! GNU General Public License, see ~abinit/COPYING
1287 : !! or http://www.gnu.org/copyleft/gpl.txt .
1288 : !!
1289 : !! INPUTS
1290 : !! this=bath operator
1291 : !! ostream=file stream to write in
1292 : !!
1293 : !! OUTPUT
1294 : !! argout(sizeout)=description
1295 : !!
1296 : !! SIDE EFFECTS
1297 : !!
1298 : !! NOTES
1299 : !!
1300 : !! SOURCE
1301 :
1302 0 : SUBROUTINE BathOperator_printM(this,ostream)
1303 :
1304 : !Arguments ------------------------------------
1305 : TYPE(BathOperator), INTENT(IN) :: this
1306 : INTEGER, OPTIONAL , INTENT(IN) :: ostream
1307 : !Local variables ------------------------------
1308 : INTEGER :: ostream_val
1309 :
1310 0 : IF ( this%activeFlavor .LE. 0 ) &
1311 0 : CALL ERROR("BathOperator_printM : no active hybrid function ")
1312 0 : ostream_val = 6
1313 0 : IF ( PRESENT(ostream) ) ostream_val = ostream
1314 0 : CALL MatrixHyb_print(this%M(this%activeFlavor),ostream_val)
1315 0 : END SUBROUTINE BathOperator_printM
1316 : !!***
1317 :
1318 : !!****f* ABINIT/m_BathOperator/ BathOperator_destroy
1319 : !! NAME
1320 : !! BathOperator_destroy
1321 : !!
1322 : !! FUNCTION
1323 : !! Deallocate and reset every thing
1324 : !!
1325 : !! COPYRIGHT
1326 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
1327 : !! This file is distributed under the terms of the
1328 : !! GNU General Public License, see ~abinit/COPYING
1329 : !! or http://www.gnu.org/copyleft/gpl.txt .
1330 : !!
1331 : !! INPUTS
1332 : !! this=bath operator
1333 : !!
1334 : !! OUTPUT
1335 : !!
1336 : !! SIDE EFFECTS
1337 : !!
1338 : !! NOTES
1339 : !!
1340 : !! SOURCE
1341 :
1342 102 : SUBROUTINE BathOperator_destroy(this)
1343 :
1344 : TYPE(BathOperator), INTENT(INOUT) :: this
1345 : INTEGER :: it
1346 :
1347 1178 : DO it = 1, this%flavors
1348 1076 : CALL MatrixHyb_destroy(this%M(it))
1349 1178 : CALL MatrixHyb_destroy(this%M_update(it))
1350 : END DO
1351 :
1352 102 : CALL Vector_destroy(this%R)
1353 102 : CALL Vector_destroy(this%Q)
1354 102 : CALL Vector_destroy(this%Rtau)
1355 102 : CALL Vector_destroy(this%Qtau)
1356 102 : FREEIF(this%F)
1357 1178 : DT_FREEIF(this%M)
1358 1178 : DT_FREEIF(this%M_update)
1359 :
1360 102 : this%MAddFlag = .FALSE.
1361 102 : this%MRemoveFlag = .FALSE.
1362 102 : this%flavors = 0
1363 102 : this%beta = 0.d0
1364 102 : this%dt = 0.d0
1365 102 : this%inv_dt = 0.d0
1366 102 : this%samples = 0
1367 102 : this%sizeHybrid = 0
1368 102 : this%activeFlavor = 0
1369 102 : this%updatePosRow = 0
1370 102 : this%updatePosCol = 0
1371 :
1372 102 : END SUBROUTINE BathOperator_destroy
1373 : !!***
1374 :
1375 : !!****f* ABINIT/m_BathOperator/BathOperator_doCheck
1376 : !! NAME
1377 : !! BathOperator_doCheck
1378 : !!
1379 : !! FUNCTION
1380 : !! Just store if we perfom check for updates of M
1381 : !!
1382 : !! COPYRIGHT
1383 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
1384 : !! This file is distributed under the terms of the
1385 : !! GNU General Public License, see ~abinit/COPYING
1386 : !! or http://www.gnu.org/copyleft/gpl.txt .
1387 : !!
1388 : !! INPUTS
1389 : !! this=bath operator
1390 : !! opt_check=second bit should be one
1391 : !!
1392 : !! OUTPUT
1393 : !!
1394 : !! SIDE EFFECTS
1395 : !!
1396 : !! NOTES
1397 : !!
1398 : !! SOURCE
1399 :
1400 49 : SUBROUTINE BathOperator_doCheck(this,opt_check)
1401 :
1402 : !Arguments ------------------------------------
1403 : TYPE(BathOperator) , INTENT(INOUT) :: this
1404 : INTEGER , INTENT(IN ) :: opt_check
1405 :
1406 49 : IF ( opt_check .GE. 2 ) &
1407 0 : this%doCheck = .TRUE.
1408 49 : END SUBROUTINE BathOperator_doCheck
1409 : !!***
1410 :
1411 : !!****f* ABINIT/m_BathOperator/BathOperator_checkM
1412 : !! NAME
1413 : !! BathOperator_checkM
1414 : !!
1415 : !! FUNCTION
1416 : !! compute from scratch the M this and compar it
1417 : !! with the already computed M this
1418 : !!
1419 : !! COPYRIGHT
1420 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
1421 : !! This file is distributed under the terms of the
1422 : !! GNU General Public License, see ~abinit/COPYING
1423 : !! or http://www.gnu.org/copyleft/gpl.txt .
1424 : !!
1425 : !! INPUTS
1426 : !! this=bath operator
1427 : !! particle=list of all segments of the active flavor
1428 : !!
1429 : !! OUTPUT
1430 : !!
1431 : !! SIDE EFFECTS
1432 : !!
1433 : !! NOTES
1434 : !!
1435 : !! SOURCE
1436 :
1437 0 : SUBROUTINE BathOperator_checkM(this,particle)
1438 :
1439 : !Arguments ------------------------------------
1440 : TYPE(BathOperator) , INTENT(INOUT) :: this
1441 : TYPE(ListCdagC) , INTENT(IN ) :: particle
1442 : !Local variables ------------------------------
1443 : ! TYPE(MatrixHyb) :: checkMatrix
1444 : LOGICAL :: checkTau
1445 : INTEGER :: tail
1446 : INTEGER :: iC
1447 : INTEGER :: iCdag
1448 : INTEGER :: aF
1449 : CHARACTER(LEN=4) :: a
1450 : DOUBLE PRECISION :: time
1451 : DOUBLE PRECISION :: beta
1452 : DOUBLE PRECISION :: mbeta_two
1453 : DOUBLE PRECISION :: erreur
1454 : DOUBLE PRECISION :: tc
1455 : DOUBLE PRECISION :: tCdag
1456 : DOUBLE PRECISION :: sumMmat
1457 : DOUBLE PRECISION :: sumCheck
1458 : #include "BathOperator_hybrid.h"
1459 :
1460 0 : aF = this%activeFlavor
1461 : !Construction de la this
1462 0 : tail = particle%tail
1463 : ! CALL MatrixHyb_init(checkMatrix,this%iTech,size=tail,Wmax=this%samples)
1464 : ! CALL MatrixHyb_setSize(checkMatrix,tail)
1465 0 : CALL MatrixHyb_setSize(this%M_update(aF),tail)
1466 0 : beta = this%beta
1467 0 : mbeta_two = -beta*0.5d0
1468 0 : this%checkNumber = this%checkNumber + 1
1469 0 : IF ( tail .NE. this%M(aF)%tail ) THEN
1470 0 : CALL WARN("BathOperator_checkM : tails are different ")
1471 0 : RETURN
1472 : END IF
1473 :
1474 : !CALL ListCdagC_print(particle)
1475 0 : DO iCdag = 1, tail
1476 0 : tCdag = particle%list(iCdag,Cdag_)
1477 0 : DO iC = 1, tail
1478 : !tC = particle%list(C_,iC).MOD.beta
1479 0 : MODCYCLE(particle%list(iC,C_),beta,tC)
1480 0 : time = tC - tCdag
1481 : #include "BathOperator_hybrid"
1482 0 : this%M_update(aF)%mat(iC,iCdag) = hybrid
1483 :
1484 0 : time = time + ( SIGN(1.d0,time) - 1.d0 )*mbeta_two
1485 0 : this%M_update(aF)%mat_tau(iCdag,iC) = INT ( (time*this%inv_dt) +1.5d0 )
1486 : END DO
1487 : END DO
1488 :
1489 : ! CALL MatrixHyb_Print(checkMatrix)
1490 : !Inversion de la this
1491 0 : CALL MatrixHyb_inverse(this%M_update(aF))
1492 : ! CALL MatrixHyb_Print(checkMatrix)
1493 :
1494 : !Comparaison
1495 0 : sumMmat =0.d0
1496 0 : sumCheck=0.d0
1497 0 : erreur = 0.d0
1498 0 : checkTau = .FALSE.
1499 0 : DO iCdag = 1, tail
1500 0 : Do iC =1, tail
1501 0 : this%M_update(aF)%mat(iC,iCdag) = ABS((this%M_update(aF)%mat(iC, iCdag) - this%M(aF)%mat(iC,iCdag))/this%M(aF)%mat(iC,iCdag))
1502 : IF ( this%M_update(aF)%mat(iC,iCdag) .GT. erreur ) erreur = this%M_update(aF)%mat(ic,iCdag)
1503 0 : IF ( this%M_update(aF)%mat_tau(iC,iCdag) .NE. this%M(aF)%mat_tau(iC,iCdag) ) checkTau = .TRUE.
1504 : END DO
1505 : END DO
1506 :
1507 0 : IF ( checkTau .EQV. .TRUE. ) THEN
1508 0 : CALL WARN("BathOperator_checkM : mat_tau differs should be")
1509 0 : CALL MatrixHyb_print(this%M_update(aF),opt_print=1)
1510 0 : CALL WARN("BathOperator_checkM : whereas it is")
1511 0 : CALL MatrixHyb_print(this%M(aF),opt_print=1)
1512 : END IF
1513 0 : this%meanError = this%meanError + erreur
1514 0 : IF ( erreur .GT. 1.d0 ) THEN
1515 0 : WRITE(a,'(I4)') INT(erreur*100.d0)
1516 : ! CALL MatrixHyb_Print(this%M(aF)
1517 0 : CALL WARN("BathOperator_checkM : "//a//"% ")
1518 : END IF
1519 : ! CALL MatrixHyb_destroy(checkMatrix)
1520 : END SUBROUTINE BathOperator_checkM
1521 : !!***
1522 :
1523 : !!****f* ABINIT/m_BathOperator/BathOperator_getError
1524 : !! NAME
1525 : !! BathOperator_getError
1526 : !!
1527 : !! FUNCTION
1528 : !! compute a percentage error / checkM
1529 : !!
1530 : !! COPYRIGHT
1531 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
1532 : !! This file is distributed under the terms of the
1533 : !! GNU General Public License, see ~abinit/COPYING
1534 : !! or http://www.gnu.org/copyleft/gpl.txt .
1535 : !!
1536 : !! INPUTS
1537 : !! this=bath operator
1538 : !!
1539 : !! OUTPUT
1540 : !! BathOperator_getError=Error in percent
1541 : !!
1542 : !! SIDE EFFECTS
1543 : !!
1544 : !! NOTES
1545 : !!
1546 : !! SOURCE
1547 :
1548 0 : DOUBLE PRECISION FUNCTION BathOperator_getError(this)
1549 :
1550 : TYPE(BathOperator), INTENT(IN) :: this
1551 :
1552 0 : IF ( this%doCheck .EQV. .TRUE. ) THEN
1553 0 : BathOperator_getError = this%meanError / DBLE(this%checkNumber)
1554 : ELSE
1555 : BathOperator_getError = 0.d0
1556 : END IF
1557 0 : END FUNCTION BathOperator_getError
1558 : !!***
1559 : !#endif
1560 :
1561 85920850 : END MODULE m_BathOperator
1562 : !!***
|
