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