Line data Source code
1 :
2 : #if defined HAVE_CONFIG_H
3 : #include "config.h"
4 : #endif
5 :
6 : !!****m* ABINIT/m_BathOperatoroffdiag
7 : !! NAME
8 : !! m_BathOperatoroffdiag
9 : !!
10 : !! FUNCTION
11 : !! Manage all stuff related to the bath for the
12 : !! simgle Anderson Impurity Model
13 : !!
14 : !! COPYRIGHT
15 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder, J. Denier, B. Amadon)
16 : !! This file is distributed under the terms of the
17 : !! GNU General Public License, see ~abinit/COPYING
18 : !! or http://www.gnu.org/copyleft/gpl.txt .
19 : !!
20 : !! NOTES
21 : !!
22 : !! SOURCE
23 :
24 : #include "defs.h"
25 : MODULE m_BathOperatoroffdiag
26 : USE m_MatrixHyb
27 : USE m_Vector
28 : USE m_VectorInt
29 : USE m_Global
30 : USE m_ListCdagC
31 : IMPLICIT NONE
32 :
33 : ! subroutines
34 : public :: BathOperatoroffdiag_init
35 : public :: BathOperatoroffdiag_reset
36 : public :: BathOperatoroffdiag_activateParticle
37 : public :: BathOperatoroffdiag_setMAdd
38 : public :: BathOperatoroffdiag_setMRemove
39 : public :: BathOperatoroffdiag_swap
40 : public :: BathOperatoroffdiag_initF
41 : public :: BathOperatoroffdiag_setF
42 : public :: BathOperatoroffdiag_printF
43 : public :: BathOperatoroffdiag_printM
44 : public :: BathOperatoroffdiag_destroy
45 : public :: BathOperatoroffdiag_doCheck
46 : public :: BathOperatoroffdiag_checkM
47 :
48 : ! functions
49 : ! public :: BathOperatoroffdiag_hybrid
50 : public :: BathOperatoroffdiag_getDetAdd
51 : public :: BathOperatoroffdiag_getDetRemove
52 : public :: BathOperatoroffdiag_getDetF
53 : public :: BathOperatoroffdiag_getError
54 : !!***
55 :
56 : !!****t* m_BathOperatoroffdiag/BathOperatoroffdiag
57 : !! NAME
58 : !! BathOperatoroffdiag
59 : !!
60 : !! FUNCTION
61 : !! This structured datatype contains the necessary data
62 : !!
63 : !! COPYRIGHT
64 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
65 : !! This file is distributed under the terms of the
66 : !! GNU General Public License, see ~abinit/COPYING
67 : !! or http://www.gnu.org/copyleft/gpl.txt .
68 : !!
69 : !! SOURCE
70 :
71 : TYPE BathOperatoroffdiag
72 : LOGICAL :: set = .FALSE.
73 : ! True if the BathOperatoroffdiag is initialized in BathOperatoroffdiag_init
74 :
75 : LOGICAL :: MAddFlag = .FALSE.
76 : ! Set to true if we can compute a new M (see updateDetXX) (ie in
77 : ! BathOperatoroffdiag_getDetAdd)
78 :
79 : LOGICAL :: MRemoveFlag = .FALSE.
80 : ! Set to true if we can compute a new M (see updateDetXX) (ie in
81 : ! BathOperatoroffdiag_getDetRemove)
82 :
83 : LOGICAL :: antiShift = .FALSE.
84 : ! shift when M is updated with antiseg
85 :
86 : LOGICAL :: doCheck = .FALSE.
87 : ! TRUE is checks are activated
88 :
89 : INTEGER :: opt_nondiag = 0
90 : ! if opt_nondiag = 1 F is non diagonal.
91 :
92 : INTEGER :: flavors
93 : ! number of flavors
94 : ! if opt_nondiag = 0 , flavors= number of flavor
95 : ! if opt_nondiag = 1 , flavors= 1
96 :
97 : INTEGER :: activeFlavor
98 : ! Active flavor on which a segment is added/suppressed...
99 :
100 : INTEGER :: samples
101 : ! Number of time slices (given in the input file)
102 :
103 : INTEGER :: sizeHybrid
104 : ! Number of time slices (given in the input file) + 1 (=qmc_l+1)
105 :
106 : INTEGER :: updatePosRow
107 : ! Gives the position of new Row to add
108 : ! Modified in BathOperatoroffdiag_getDetAdd and BathOperatoroffdiag_getDetRemove
109 : ! could be the Row in the Full matrix for the non diag implementation
110 :
111 : INTEGER :: updatePosCol
112 : ! Gives the position of new Col to add
113 : ! Modified in BathOperatoroffdiag_getDetAdd and BathOperatoroffdiag_getDetRemove
114 :
115 : INTEGER :: iTech
116 : ! iTech is an integer which precise the technics used to compute the
117 : ! Green's function (in time or frequency)
118 :
119 : INTEGER :: sumtails
120 : ! size of the full F matrix (sums of tails(iflavor) over iflavor)
121 :
122 : INTEGER, ALLOCATABLE, DIMENSION(:) :: tails
123 : ! tails(iflavor) is the current number of segments for the flavor iflavor
124 :
125 : INTEGER, ALLOCATABLE, DIMENSION(:) :: Fshift
126 : ! Fshift(iflavor) is the sum of number of segments for all flavors iflavor2
127 : ! such that iflavor< iflavor
128 : ! It is thus the shift in the F matrix to have the first segment of the flavor
129 : ! iflavor
130 : ! Fshift(nflavor+1) is the total nb of tails (=sumtails)
131 :
132 : DOUBLE PRECISION :: beta
133 : ! Inverse of Temperature
134 : !
135 :
136 : DOUBLE PRECISION :: dt
137 : ! dt=beta/samples
138 :
139 : DOUBLE PRECISION :: inv_dt
140 : ! inv_dt=1/dt
141 :
142 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) :: F ! qmc_l+2,Flavors
143 : ! Hybridization function F(1:op%sizeHybrid+1,1:flavors,1:flavors)
144 :
145 : DOUBLE PRECISION :: S
146 : ! Sherman Morrison notations
147 :
148 : DOUBLE PRECISION :: Stau
149 : ! Sherman Morrison notations
150 :
151 : DOUBLE PRECISION :: Stilde
152 : ! Sherman Morrison notations
153 :
154 : TYPE(Vector) :: R
155 : ! Sherman Morrison notations R%vec(size).
156 : ! computed for each flavor (As matrices are made of Blocks for each
157 : ! flavor because the code is restricted to diagonal F matrices)
158 :
159 : TYPE(Vector) :: Q
160 : ! Sherman Morrison notations
161 : ! computed for each flavor (As matrices are made of Blocks for each
162 : ! flavor because the code is restricted to diagonal F matrices)
163 :
164 : TYPE(Vector) :: Rtau
165 : ! Sherman Morrison notations
166 : ! Rtau gives the time length for each elements of R
167 : ! computed for each flavor (As matrices are made of Blocks for each
168 : ! flavor because the code is restricted to diagonal F matrices)
169 :
170 : TYPE(Vector) :: Qtau
171 : ! Sherman Morrison notations
172 : ! Qtau gives the time length for each elements of Q
173 : ! computed for each flavor (As matrices are made of Blocks for each
174 : ! flavor because the code is restricted to diagonal F matrices)
175 :
176 : TYPE(MatrixHyb) :: M ! Flavors
177 : ! inverse of Hybridization matrix M%mat(global_size,global_size)
178 : ! contains the value of the hybridization for all flavor and segments times,
179 : ! the times (mat_tau), and possibly the
180 : ! frequency
181 :
182 : TYPE(MatrixHyb) :: M_update ! Flavors
183 : ! used in BathOperatoroffdiag_getdetF and in BathOperatoroffdiag_checkM
184 : ! for checks
185 :
186 : !#ifdef CTQMC_CHECK
187 : INTEGER :: checkNumber
188 : DOUBLE PRECISION :: meanError
189 : ! TYPE(ListCdagC) :: ListCdagC
190 : !#endif
191 : END TYPE BathOperatoroffdiag
192 : !!***
193 :
194 : CONTAINS
195 : !!***
196 :
197 : !!****f* ABINIT/m_BathOperatoroffdiag/BathOperatoroffdiag_init
198 : !! NAME
199 : !! BathOperatoroffdiag_init
200 : !!
201 : !! FUNCTION
202 : !! Initialize and allocate data
203 : !!
204 : !! COPYRIGHT
205 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
206 : !! This file is distributed under the terms of the
207 : !! GNU General Public License, see ~abinit/COPYING
208 : !! or http://www.gnu.org/copyleft/gpl.txt .
209 : !!
210 : !! INPUTS
211 : !! op=bath object
212 : !! flavors=numbers of flavors we have (including spin)
213 : !! samples=Time slices in the input file
214 : !! beta=inverse temperature
215 : !! iTech=imaginary time or frequencies
216 : !! It is imposes to imaginary time
217 : !!
218 : !! OUTPUT
219 : !!
220 : !! SIDE EFFECTS
221 : !!
222 : !! NOTES
223 : !!
224 : !! SOURCE
225 :
226 0 : SUBROUTINE BathOperatoroffdiag_init(op, flavors, samples, beta, iTech,opt_nondiag)
227 :
228 : !Arguments ------------------------------------
229 : TYPE(BathOperatoroffdiag), INTENT(INOUT) :: op
230 : INTEGER , INTENT(IN ) :: flavors
231 : INTEGER , INTENT(IN ) :: samples
232 : INTEGER , INTENT(IN ) :: opt_nondiag
233 : DOUBLE PRECISION , INTENT(IN ) :: beta
234 : !Local variables ------------------------------
235 : INTEGER , INTENT(IN ) :: iTech
236 : !INTEGER :: it
237 :
238 0 : op%MAddFlag = .FALSE.
239 0 : op%MRemoveFlag = .FALSE.
240 0 : op%flavors = flavors
241 0 : op%opt_nondiag = opt_nondiag
242 0 : op%beta = beta
243 0 : op%samples = samples
244 0 : op%sizeHybrid = samples + 1
245 0 : op%dt = beta / DBLE(samples)
246 0 : op%inv_dt = DBLE(samples) / beta
247 0 : op%activeFlavor= 0
248 0 : op%updatePosRow = 0
249 0 : op%updatePosCol = 0
250 0 : op%iTech = iTech
251 : !#ifdef CTQMC_CHECK
252 0 : op%checkNumber = 0
253 0 : op%meanError = 0.d0
254 0 : op%doCheck = .FALSE.
255 : !#endif
256 :
257 0 : FREEIF(op%F)
258 0 : MALLOC(op%F,(1:op%sizeHybrid+1,1:flavors,1:flavors))
259 0 : DT_FREEIF(op%tails)
260 : #ifdef FC_LLVM
261 : ! LLVM 16 doesn't recognize this macro here
262 : DT_MALLOC(op%tails, (1:op%flavors))
263 : #else
264 0 : DT_MALLOC(op%tails,(1:op%flavors))
265 : #endif
266 0 : op%tails=0
267 0 : DT_FREEIF(op%Fshift)
268 : #ifdef FC_LLVM
269 : ! LLVM 16 doesn't recognize this macro here
270 : DT_MALLOC(op%Fshift, (1:op%flavors+1))
271 : #else
272 0 : DT_MALLOC(op%Fshift,(1:op%flavors+1))
273 : #endif
274 0 : op%Fshift=0
275 :
276 0 : CALL Vector_init(op%R,100*op%flavors)
277 0 : CALL Vector_init(op%Q,100*op%flavors)
278 0 : CALL Vector_init(op%Rtau,100*op%flavors)
279 0 : CALL Vector_init(op%Qtau,100*op%flavors)
280 :
281 0 : CALL MatrixHyb_init(op%M,op%iTech,size=Global_SIZE*op%flavors,Wmax=samples) !FIXME Should be consistent with ListCagC
282 0 : CALL MatrixHyb_init(op%M_update,op%iTech,size=Global_SIZE*op%flavors,Wmax=samples) !FIXME Should be consistent with ListCagC
283 0 : op%F = 0.d0
284 0 : op%set = .TRUE.
285 :
286 0 : END SUBROUTINE BathOperatoroffdiag_init
287 : !!***
288 :
289 : !!****f* ABINIT/m_BathOperatoroffdiag/BathOperatoroffdiag_reset
290 : !! NAME
291 : !! BathOperatoroffdiag_reset
292 : !!
293 : !! FUNCTION
294 : !! Reset all internal variables
295 : !!
296 : !! COPYRIGHT
297 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
298 : !! This file is distributed under the terms of the
299 : !! GNU General Public License, see ~abinit/COPYING
300 : !! or http://www.gnu.org/copyleft/gpl.txt .
301 : !!
302 : !! INPUTS
303 : !! op=bath operator to reset
304 : !!
305 : !! OUTPUT
306 : !!
307 : !! SIDE EFFECTS
308 : !!
309 : !! NOTES
310 : !!
311 : !! SOURCE
312 :
313 0 : SUBROUTINE BathOperatoroffdiag_reset(op)
314 :
315 : !Arguments ------------------------------------
316 : TYPE(BathOperatoroffdiag), INTENT(INOUT) :: op
317 : !Local variables ------------------------------
318 : INTEGER :: iflavor
319 0 : op%MAddFlag = .FALSE.
320 0 : op%MRemoveFlag = .FALSE.
321 0 : op%activeFlavor = 0
322 0 : op%updatePosRow = 0
323 0 : op%updatePosCol = 0
324 : !#ifdef CTQMC_CHECK
325 0 : op%checkNumber = 0
326 0 : op%meanError = 0.d0
327 0 : op%sumtails = 0
328 : !#endif
329 0 : op%doCheck = .FALSE.
330 0 : CALL Vector_clear(op%R)
331 0 : CALL Vector_clear(op%Q)
332 0 : CALL Vector_clear(op%Rtau)
333 0 : CALL Vector_clear(op%Qtau)
334 :
335 0 : CALL MatrixHyb_clear(op%M) !FIXME Should be consistent with ListCagC
336 0 : op%F = 0.d0
337 0 : do iflavor=1,op%flavors
338 0 : op%tails(iflavor)=0
339 0 : op%Fshift(iflavor)=0
340 : enddo
341 :
342 0 : END SUBROUTINE BathOperatoroffdiag_reset
343 : !!***
344 :
345 : !!****f* ABINIT/m_BathOperatoroffdiag/BathOperatoroffdiag_activateParticle
346 : !! NAME
347 : !! BathOperatoroffdiag_activateParticle
348 : !!
349 : !! FUNCTION
350 : !! Just save on wicht flavor we are working
351 : !! It is better to use the macro defined in defs.h
352 : !!
353 : !! COPYRIGHT
354 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
355 : !! This file is distributed under the terms of the
356 : !! GNU General Public License, see ~abinit/COPYING
357 : !! or http://www.gnu.org/copyleft/gpl.txt .
358 : !!
359 : !! INPUTS
360 : !! op=bath operator
361 : !! flavor=the flavor to activate
362 : !!
363 : !! OUTPUT
364 : !!
365 : !! SIDE EFFECTS
366 : !!
367 : !! NOTES
368 : !!
369 : !! SOURCE
370 :
371 0 : SUBROUTINE BathOperatoroffdiag_activateParticle(op,flavor)
372 :
373 : !Arguments ------------------------------------
374 : TYPE(BathOperatoroffdiag), INTENT(INOUT) :: op
375 : !Local variables ------------------------------
376 : INTEGER , INTENT(IN ) :: flavor
377 :
378 0 : IF ( flavor .GT. op%flavors ) &
379 0 : CALL ERROR("BathOperatoroffdiag_activateParticle : out of range ")
380 0 : IF ( op%set .EQV. .TRUE. ) THEN
381 0 : op%activeFlavor = flavor
382 0 : op%MAddFlag = .FALSE.
383 0 : op%MRemoveFlag = .FALSE.
384 : ELSE
385 0 : CALL ERROR("BathOperatoroffdiag_activateParticle : not allocated ")
386 : END IF
387 0 : END SUBROUTINE BathOperatoroffdiag_activateParticle
388 : !!***
389 :
390 : !!****f* ABINIT/m_BathOperatoroffdiag/BathOperatoroffdiag_getDetAdd
391 : !! NAME
392 : !! BathOperatoroffdiag_getDetAdd
393 : !!
394 : !! FUNCTION
395 : !! Compute the determinant ratio when a (anti)segment
396 : !! is trying to be added and store some array for setMAdd
397 : !!
398 : !! COPYRIGHT
399 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
400 : !! This file is distributed under the terms of the
401 : !! GNU General Public License, see ~abinit/COPYING
402 : !! or http://www.gnu.org/copyleft/gpl.txt .
403 : !!
404 : !! INPUTS
405 : !! op=bath operator
406 : !! CdagC_1=segment to be added
407 : !! position=ordered position of the Cdag time
408 : !! particle=full list of CdagC for activeFlavor
409 : !!
410 : !! OUTPUT
411 : !! BathOperatoroffdiag_getDetAdd=the det
412 : !!
413 : !! SIDE EFFECTS
414 : !!
415 : !! NOTES
416 : !!
417 : !! SOURCE
418 0 : DOUBLE PRECISION FUNCTION BathOperatoroffdiag_getDetAdd(op,CdagC_1, position, particle)
419 :
420 : !Arguments ------------------------------------
421 : TYPE(BathOperatoroffdiag) , INTENT(INOUT) :: op
422 : DOUBLE PRECISION, DIMENSION(1:2), INTENT(IN ) :: CdagC_1
423 : INTEGER , INTENT(IN ) :: position
424 : TYPE(ListCdagC), INTENT(IN ) :: particle(:)
425 : !Local variables-------------------------------
426 : INTEGER :: it1
427 : INTEGER :: it2
428 : INTEGER :: it3,iflavor,iflavora,iflavorb
429 : INTEGER :: iflavorbegin,iflavorend
430 : INTEGER :: tail,tailbegin,tailend
431 : INTEGER :: tcheck
432 : INTEGER :: new_tail
433 : DOUBLE PRECISION :: C
434 : DOUBLE PRECISION :: Cbeta
435 : DOUBLE PRECISION :: Cibeta
436 : DOUBLE PRECISION :: Cdag
437 : DOUBLE PRECISION :: Cdagbeta
438 : DOUBLE PRECISION :: beta
439 : DOUBLE PRECISION :: ratio
440 : DOUBLE PRECISION :: time
441 : ! TYPE(CdagC) , POINTER, DIMENSION(:) :: list => NULL()
442 : #include "BathOperatoroffdiag_hybrid.h"
443 :
444 0 : op%antiShift = .FALSE.
445 0 : beta = op%beta
446 0 : C = CdagC_1(C_)
447 : ! Cbeta = C.MOD.beta
448 0 : MODCYCLE(C,beta,Cbeta)
449 0 : Cdag = CdagC_1(Cdag_)
450 : ! cdagbeta = Cdag.MOD.beta
451 0 : MODCYCLE(Cdag,beta,Cdagbeta)
452 : ! IF ( Cdag .GE. beta ) &
453 : ! CALL ERROR("BathOperatoroffdiag_getDetAdd : bad case ... ")
454 0 : IF ( op%activeFlavor .LE. 0 ) &
455 0 : CALL ERROR("BathOperatoroffdiag_getDetAdd : no active hybrid function ")
456 :
457 0 : IF ( size(particle)/=op%flavors ) &
458 0 : CALL ERROR("BathOperatoroffdiag_getDetAdd : size of particle is erroneous ")
459 :
460 : ! tail is now the complete size of the F matrix Fshift(nflavors+1)
461 0 : tail = op%sumtails
462 0 : new_tail = tail+1
463 : ! list => particle%list
464 :
465 0 : if(op%opt_nondiag==1) then
466 0 : iflavorbegin = 1
467 0 : iflavorend = op%flavors
468 0 : tailbegin = 1
469 0 : tailend = tail
470 : else
471 : !sui!write(6,*) "Bathoperator opt_nondiag=0"
472 0 : iflavorbegin = op%activeflavor
473 0 : iflavorend = op%activeflavor
474 0 : tailbegin = op%Fshift(op%activeflavor)+1
475 0 : tailend = op%Fshift(op%activeflavor)+op%tails(op%activeflavor)
476 : endif
477 :
478 : IF ( ((C .GT. Cdag) .AND. (position .EQ. -1)) & ! Segment added at the end of the segment
479 0 : .OR. ((C .LT. Cdag) .AND. (tail .EQ. 0))) THEN ! empty orbital case: only adding a segment is possible
480 : ! If ones add a segment to an empty orbital or a segment at the end
481 : ! of a segment, then:
482 0 : op%updatePosRow = op%tails(op%activeFlavor) + 1
483 0 : op%updatePosCol = op%tails(op%activeFlavor) + 1
484 : ELSE
485 : ! For all the other cases, ABS(position) is the true position.
486 0 : op%updatePosRow = ABS(position)
487 0 : op%updatePosCol = ABS(position)
488 : END IF
489 : !write(6,*) " BathOperatoroffdiag_getDetAdd : op%updatePosRow",op%updatePosRow
490 : !write(6,*) " BathOperatoroffdiag_getDetAdd : op%updatePosCol",op%updatePosCol
491 : !write(6,*) " BathOperatoroffdiag_getDetAdd : C,Cdag",C,Cdag
492 :
493 0 : IF ( C .LT. Cdag .AND. op%tails(op%activeFlavor) .GT. 0) THEN ! only if an antisegment is added
494 : ! ratio = -ratio
495 0 : op%updatePosRow = (op%updatePosRow + 1) !position in [1;tail]
496 : ! If the antisegment created is such that a segment with tcdagger> tc
497 : ! is suppressed
498 : !write(6,*) " BathOperatoroffdiag_getDetAdd : op%updatePosRow",op%updatePosRow
499 : !write(6,*) " BathOperatoroffdiag_getDetAdd : op%updatePosCol",op%updatePosCol
500 0 : IF ( Cdagbeta .LT. particle(op%activeFlavor)%list(op%updatePosCol,Cdag_) ) op%antiShift = .TRUE.
501 : END IF
502 :
503 : ! CALL Vector_setSize(op%R,tail)
504 : ! CALL Vector_setSize(op%Q,tail)
505 0 : Vector_QuickResize(op%R,new_tail)
506 0 : Vector_QuickResize(op%Q,new_tail)
507 0 : Vector_QuickResize(op%Rtau,new_tail)
508 0 : Vector_QuickResize(op%Qtau,new_tail)
509 :
510 : ! This loop compute all Row and Col except op%updatePosRow
511 0 : tcheck=0
512 0 : DO iflavor = iflavorbegin,iflavorend
513 : !write(6,*) " BathOperatoroffdiag_getDetAdd : tails(iflavor)",iflavor,op%tails(iflavor)
514 0 : DO it1 = 1, op%tails(iflavor)
515 0 : tcheck=tcheck+1
516 0 : it2 = it1
517 0 : it3 = it1
518 0 : IF ( iflavor .GE. op%activeFlavor ) THEN
519 0 : it2 = it1 + 1
520 0 : it3 = it1 + 1
521 0 : IF ( iflavor .EQ. op%activeFlavor .AND. it1 .LT. op%updatePosRow ) it2 = it1
522 0 : IF ( iflavor .EQ. op%activeFlavor .AND. it1 .LT. op%updatePosCol ) it3 = it1
523 : !it2 = it1 + ( 1+SIGN(1,it1-op%updatePosRow) )/2
524 : !it3 = it1 + ( 1+SIGN(1,it1-op%updatePoscol) )/2
525 : ! if it1>=op%updatePosRow and iflavor> activeflavor, then it2=it1+1
526 : ! if it1< op%updatePosRow and iflavor> activeflavor, then it2=it1
527 : END IF
528 :
529 : !!write(6,*) size(op%Rtau%vec)
530 : !!write(6,*) size(particle(iflavor)%list,1)
531 : !!write(6,*) size(particle(iflavor)%list,2)
532 : !!write(6,*) size(op%Fshift)
533 : !!write(6,*) it1,Cdag_,op%Fshift(iflavor)+it2
534 0 : op%Rtau%vec(op%Fshift(iflavor)+it2)= C - particle(iflavor)%list(it1,Cdag_)
535 : ! the following line happend only for nondiag case
536 0 : IF(op%Rtau%vec(op%Fshift(iflavor)+it2) .GT. beta) op%Rtau%vec(op%Fshift(iflavor)+it2)=op%Rtau%vec(op%Fshift(iflavor)+it2)-beta
537 : !op%Rtau%vec(it1)= C - particle%list(it1,Cdag_)
538 0 : time = Cbeta - particle(iflavor)%list(it1,Cdag_)
539 0 : if(op%Rtau%vec(op%Fshift(iflavor)+it2)>beta) then
540 : !write(6,*) "Rtau sup beta",op%Rtau%vec(op%Fshift(iflavor)+it2),C,particle(iflavor)%list(it1,Cdag_)
541 : !write(6,*) time
542 0 : stop
543 : endif
544 :
545 : ! "BathOperatoroffdiag_hybrid" interpolates between known values of F for the
546 : ! selected time.
547 0 : iflavora=iflavor
548 0 : iflavorb=op%activeFlavor
549 : #include "BathOperatoroffdiag_hybrid"
550 :
551 0 : op%R%vec(op%Fshift(iflavor)+it1) = hybrid
552 : ! op%R%vec(it) = BathOperatoroffdiag_hybrid(op, Cbeta - list(it)%Cdag)
553 : ! Cibeta = list(it)%C.MOD.beta
554 0 : MODCYCLE(particle(iflavor)%list(it1,C_),beta,Cibeta)
555 0 : time = Cibeta - Cdagbeta
556 0 : op%Qtau%vec(op%Fshift(iflavor)+it3)= time
557 : !op%Qtau%vec(it1)= time
558 :
559 0 : iflavora=op%activeFlavor
560 0 : iflavorb=iflavor
561 : #include "BathOperatoroffdiag_hybrid"
562 0 : op%Q%vec(op%Fshift(iflavor)+it1) = hybrid
563 :
564 : !op%Q%vec(it3) = hybrid
565 : ! Q(it) = BathOperatoroffdiag_hybrid(op, Cibeta - Cdagbeta)
566 : END DO
567 : END DO
568 0 : if(tcheck.ne.tail) then
569 : !write(6,*) " PRB in the loop tail tcheck",tail,tcheck
570 0 : stop
571 : endif
572 :
573 : ! Compute S
574 0 : op%Stau = C - Cdagbeta
575 0 : op%Rtau%vec(op%Fshift(op%activeFlavor)+op%updatePosRow) = op%Stau
576 0 : if(op%Rtau%vec(op%Fshift(op%activeFlavor)+op%updatePosRow)>beta) then
577 : !write(6,*) "Rtau sup beta", op%Stau,C,Cdagbeta
578 0 : stop
579 : endif
580 0 : op%Qtau%vec(op%Fshift(op%activeFlavor)+op%updatePosCol) = op%Rtau%vec(op%Fshift(op%activeFlavor)+op%updatePosRow)
581 : !write(6,*) " getdetAdd op%Stau",op%Stau
582 :
583 0 : time = Cbeta-Cdagbeta
584 : !write(6,*) " getdetAdd time",time
585 0 : iflavora=op%activeFlavor
586 0 : iflavorb=op%activeFlavor
587 : !write(6,*) "time",time
588 : #include "BathOperatoroffdiag_hybrid"
589 : !write(6,*) "hybrid",hybrid
590 0 : op%S = hybrid
591 : !write(6,*) " getdetAdd hybrid=op%S",hybrid
592 :
593 : !ratio = op%S - DOT_PRODUCT(MATMUL(op%R%vec(1:tail),op%M(op%activeFlavor)%mat(1:tail,1:tail)),op%Q%vec(1:tail))
594 :
595 : ! product of matrix R and M(k) is computed now:
596 0 : ratio = 0.d0
597 0 : DO it1 = tailbegin, tailend
598 : time = 0.d0
599 0 : DO it2 = tailbegin, tailend
600 0 : time = time + op%R%vec(it2) * op%M%mat(it2,it1)
601 : END DO
602 0 : ratio = ratio + op%Q%vec(it1) * time
603 : END DO
604 : !sui!write(6,*) " = R Matrix",tail
605 : !sui!write(6,*) " R ",(op%R%vec(it1),it1=1,tail)
606 : !sui!write(6,*) " = Q Matrix",tail
607 : !sui!do it1=1,tail
608 : !sui!write(6,*) " Q ",op%Q%vec(it1)
609 : !sui!enddo
610 : !sui!write(6,*) " = M Matrix",tail
611 : !sui!do it2=1,tail
612 : !sui!write(6,*) " M ",(op%M%mat(it2,it1),it1=1,tail)
613 : !sui!enddo
614 : !sui!write(6,*) " RMQ =", ratio
615 : !sui!write(6,*) " S =", op%S
616 0 : ratio = op%S - ratio
617 : !sui!write(6,*) " S-RMQ =", ratio
618 : !sui!write(6,*) " getdetAdd ratio",ratio
619 :
620 0 : op%Stilde = 1.d0 / ratio
621 : ! If antisegment, the det ratio has to be multiplied by -1 ( sign of the signature of one
622 : ! permutation line in the matrix)
623 0 : IF ( C .LT. Cdag .AND. op%tails(op%activeFlavor) .GT. 0) THEN ! only if an antisegment is added
624 0 : ratio=-ratio
625 : ENDIF
626 :
627 : ! This IF is the LAST "NON CORRECTION" in my opinion this should not appears.
628 : ! IF ( MAX(C,Cdag) .GT. op%beta ) THEN
629 : ! WRITE(*,*) op%Stilde
630 : ! op%Stilde = - ABS(op%Stilde)
631 : ! END IF
632 0 : BathOperatoroffdiag_getDetAdd = ratio
633 : !write(6,*) " getdetAdd",ratio,BathOperatoroffdiag_getDetAdd
634 0 : op%MAddFlag = .TRUE.
635 : !#ifdef CTQMC_CHECK
636 : ! op%ListCdagC = particle
637 : !!write(*,*) op%Stilde
638 : !!write(*,*) op%antishift
639 : !!write(*,*) op%updatePosRow
640 : !!write(*,*) op%updatePosCol
641 : !#endif
642 :
643 0 : END FUNCTION BathOperatoroffdiag_getDetAdd
644 : !!***
645 :
646 : !!****f* ABINIT/m_BathOperatoroffdiag/BathOperatoroffdiag_getDetRemove
647 : !! NAME
648 : !! BathOperatoroffdiag_getDetRemove
649 : !!
650 : !! FUNCTION
651 : !! Compute the determinant ratio when a (anti)segment
652 : !! is trying to be removed
653 : !!
654 : !! COPYRIGHT
655 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
656 : !! This file is distributed under the terms of the
657 : !! GNU General Public License, see ~abinit/COPYING
658 : !! or http://www.gnu.org/copyleft/gpl.txt .
659 : !!
660 : !! INPUTS
661 : !! op=bath operator
662 : !! position=position of segment to be removed
663 : !!
664 : !! OUTPUT
665 : !! BathOperatoroffdiag_getDetRemove=the det
666 : !!
667 : !! SIDE EFFECTS
668 : !!
669 : !! NOTES
670 : !!
671 : !! SOURCE
672 :
673 0 : DOUBLE PRECISION FUNCTION BathOperatoroffdiag_getDetRemove(op,position)
674 :
675 : !Arguments ------------------------------------
676 : TYPE(BathOperatoroffdiag), INTENT(INOUT) :: op
677 : !Local arguments-------------------------------
678 : INTEGER , INTENT(IN ) :: position
679 : INTEGER :: ABSposition
680 : INTEGER :: tail !,it,it1
681 :
682 0 : IF ( op%activeFlavor .LE. 0 ) &
683 0 : CALL ERROR("BathOperatoroffdiag_getDetRemove : no active hybrid fun ")
684 :
685 0 : op%antiShift = .FALSE.
686 0 : tail = op%sumtails
687 0 : ABSposition = ABS(position)
688 0 : IF ( ABSposition .GT. op%tails(op%activeFlavor) ) &
689 0 : CALL ERROR("BathOperatoroffdiag_getDetRemove : position > M size ")
690 0 : op%updatePosCol = ABSposition
691 0 : op%antiShift = .FALSE.
692 0 : IF ( position .GT. 0 ) THEN
693 0 : op%updatePosRow = ABSposition
694 : ELSE
695 0 : op%updatePosRow = ABSposition+1
696 0 : IF ( ABSposition .EQ. op%tails(op%activeFlavor) ) THEN
697 0 : op%antiShift = .TRUE.
698 0 : op%updatePosRow = 1 !ABSposition - 1
699 : ! op%updatePosRow = ABSposition
700 : ! IF ( op%updatePosCol .EQ. 0) op%updatePosCol = tail
701 : END IF
702 : ENDIF
703 : op%Stilde = op%M%mat(op%Fshift(op%activeFlavor)+&
704 0 : & op%updatePosRow,op%Fshift(op%activeFlavor)+op%updatePosCol)
705 : !sui!write(6,*) "Fshift",op%Fshift(op%activeFlavor)
706 : !sui!write(6,*) "updatepos",op%updatePosRow,op%updatePosCol
707 :
708 :
709 0 : op%MRemoveFlag = .TRUE.
710 : !write(6,*) " getdetRemove",op%Stilde
711 0 : BathOperatoroffdiag_getDetRemove = op%Stilde
712 0 : if(position<0.and.op%tails(op%activeFlavor)>1) then
713 0 : BathOperatoroffdiag_getDetRemove = -op%Stilde
714 : endif
715 : !do it=1,op%sumtails
716 : !!sui!write(6,*) " getdetRemove M",(op%M%mat(it,it1),it1=1,op%sumtails)
717 : !enddo
718 : !#ifdef CTQMC_CHECK
719 : ! op%ListCdagC = particle
720 : !!write(*,*) op%updatePosRow, op%updatePosCol, position
721 : !!CALL ListCdagC_print(particle)
722 : !#endif
723 :
724 0 : END FUNCTION BathOperatoroffdiag_getDetRemove
725 : !!***
726 :
727 : !!****f* ABINIT/m_BathOperatoroffdiag/BathOperatoroffdiag_getDetF
728 : !! NAME
729 : !! BathOperatoroffdiag_getDetF
730 : !!
731 : !! FUNCTION
732 : !! Compute the determinant of the F matrix
733 : !! using the hybridization of flavor and the
734 : !! segments of particle
735 : !! used for Gloval moves only
736 : !!
737 : !! COPYRIGHT
738 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
739 : !! This file is distributed under the terms of the
740 : !! GNU General Public License, see ~abinit/COPYING
741 : !! or http://www.gnu.org/copyleft/gpl.txt .
742 : !!
743 : !! INPUTS
744 : !! op=bath operator
745 : !! flavor=hybridization function to take
746 : !! particles=segments to use
747 : !!
748 : !! OUTPUT
749 : !! BathOperatoroffdiag_getDetF=the det
750 : !!
751 : !! SIDE EFFECTS
752 : !!
753 : !! NOTES
754 : !!
755 : !! SOURCE
756 :
757 0 : DOUBLE PRECISION FUNCTION BathOperatoroffdiag_getDetF(op,particle,option)
758 :
759 : !Arguments ------------------------------------
760 : TYPE(BathOperatoroffdiag) , INTENT(INOUT) :: op
761 : TYPE(ListCdagC), OPTIONAL, INTENT(IN ) :: particle(:)
762 : INTEGER , optional :: option
763 : !Local arguments-------------------------------
764 : INTEGER :: iCdag
765 : INTEGER :: iC
766 : INTEGER :: tail
767 : DOUBLE PRECISION :: time
768 : DOUBLE PRECISION :: tC
769 : DOUBLE PRECISION :: tCdag
770 : DOUBLE PRECISION :: beta
771 : DOUBLE PRECISION :: mbeta_two
772 : DOUBLE PRECISION :: signe
773 : DOUBLE PRECISION :: inv_dt
774 : INTEGER :: iflavor,iflavora
775 : INTEGER :: iflavordag,iflavorb
776 : #include "BathOperatoroffdiag_hybrid.h"
777 :
778 0 : BathOperatoroffdiag_getDetF = 1.d0 ! pour eviter des divisions par 0
779 0 : IF ( PRESENT( particle ) ) THEN
780 0 : tail = op%sumtails
781 0 : beta = op%beta
782 0 : mbeta_two = -beta*0.5d0
783 0 : inv_dt = op%inv_dt
784 0 : CALL MatrixHyb_setSize(op%M_update,tail)
785 0 : DO iflavordag=1,op%flavors
786 0 : DO iCdag = 1, op%tails(iflavordag)
787 0 : tCdag = particle(iflavordag)%list(iCdag,Cdag_)
788 0 : DO iflavor=1,op%flavors
789 0 : DO iC = 1, op%tails(iflavor)
790 : !tC = particle%list(C_,iC).MOD.beta
791 0 : MODCYCLE(particle(iflavor)%list(iC,C_),beta,tC)
792 0 : time = tC - tCdag
793 0 : iflavora=iflavordag
794 0 : iflavorb=iflavor
795 : #include "BathOperatoroffdiag_hybrid"
796 0 : op%M_update%mat(op%Fshift(iflavor)+iC,op%Fshift(iflavordag)+iCdag) = hybrid
797 : END DO
798 : END DO
799 : END DO
800 : END DO
801 : ! mat_tau needs to be transpose of ordered time mat (way of measuring
802 : ! G(tau))
803 0 : DO iflavor=1,op%flavors
804 0 : DO iC = 1, tail
805 0 : tC = particle(iflavor)%list(iC,C_)
806 0 : DO iflavordag=1,op%flavors
807 0 : DO iCdag = 1, tail
808 : !sui!write(6,*) iCdag,Cdag_,size(particle(iflavordag)%list,1)
809 : !stop
810 0 : tCdag = particle(iflavordag)%list(iCdag,Cdag_)
811 0 : time = tC - tCdag
812 0 : signe = SIGN(1.d0,time)
813 0 : time = time + (signe-1.d0)*mbeta_two
814 0 : op%M_update%mat_tau(op%Fshift(iflavordag)+iCdag,op%Fshift(iflavor)+iC) = INT( ( time * inv_dt ) + 1.5d0 )
815 : END DO
816 : END DO
817 : END DO
818 : END DO
819 0 : CALL MatrixHyb_inverse(op%M_update,BathOperatoroffdiag_getDetF) ! calcul le det de la matrice et l'inverse
820 : ELSE
821 0 : if(present(option)) then
822 0 : CALL MatrixHyb_getDet(op%M_update,BathOperatoroffdiag_getDetF) ! det M = 1/detF !
823 : else
824 0 : CALL MatrixHyb_getDet(op%M,BathOperatoroffdiag_getDetF) ! det M = 1/detF !
825 : endif
826 0 : BathOperatoroffdiag_getDetF = 1.d0 / BathOperatoroffdiag_getDetF
827 : ENDIF
828 0 : END FUNCTION BathOperatoroffdiag_getDetF
829 : !!***
830 :
831 : !!****f* ABINIT/m_BathOperatoroffdiag/BathOperatoroffdiag_setMAdd
832 : !! NAME
833 : !! BathOperatoroffdiag_setMAdd
834 : !!
835 : !! FUNCTION
836 : !! Update de M matrix inserting a row and a column
837 : !!
838 : !! COPYRIGHT
839 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
840 : !! This file is distributed under the terms of the
841 : !! GNU General Public License, see ~abinit/COPYING
842 : !! or http://www.gnu.org/copyleft/gpl.txt .
843 : !!
844 : !! INPUTS
845 : !! op=bath operator
846 : !! particle=segments of active flavor
847 : !!
848 : !! OUTPUT
849 : !!
850 : !! SIDE EFFECTS
851 : !!
852 : !! NOTES
853 : !!
854 : !! SOURCE
855 :
856 0 : SUBROUTINE BathOperatoroffdiag_setMAdd(op,particle)
857 :
858 : !Arguments ------------------------------------
859 : TYPE(BathOperatoroffdiag), INTENT(INOUT) :: op
860 : TYPE(ListCdagC) , INTENT(IN ) :: particle(:)
861 : !Local variables ------------------------------
862 : INTEGER :: tail
863 : INTEGER :: new_tail
864 : INTEGER :: col
865 : INTEGER :: col_move
866 : INTEGER :: row_move
867 : INTEGER :: row
868 : INTEGER :: positionRow
869 : INTEGER :: positionCol
870 : INTEGER :: aF,indice
871 : INTEGER :: tailb,taile
872 : DOUBLE PRECISION :: Stilde
873 : DOUBLE PRECISION :: time
874 : DOUBLE PRECISION :: mbeta_two
875 : DOUBLE PRECISION :: inv_dt
876 0 : TYPE(Vector) :: vec_tmp
877 0 : TYPE(VectorInt) :: vecI_tmp
878 : INTEGER :: m
879 : INTEGER :: count
880 : INTEGER :: i
881 : INTEGER :: j
882 : INTEGER :: p,it !,it1
883 :
884 :
885 :
886 : ! --- op%MAddFlag is put to .TRUE. in BathOperatoroffdiag_getDetAdd.
887 0 : IF ( op%MAddFlag .EQV. .FALSE. ) &
888 0 : CALL ERROR("BathOperatoroffdiag_setMAdd : MAddFlag turn off ")
889 :
890 : ! --- op%activeFlavor is put in ctqmc_loop
891 0 : aF = op%activeFlavor
892 0 : IF ( aF .LE. 0 ) &
893 0 : CALL ERROR("BathOperatoroffdiag_setMAdd : no active hybrid function ")
894 :
895 : !! do it=1,op%sumtails
896 : !write(6,*) " setMAdd begin M",(op%M%mat(it,it1),it1=1,op%sumtails)
897 : !! enddo
898 : !! do it=1,op%sumtails
899 : !write(6,*) " setMAdd begin M%mat_tau",(op%M%mat_tau(it,it1),it1=1,op%sumtails)
900 : !! enddo
901 : ! old tail
902 : !write(6,*) " BathOperatoroffdiag_setMAdd op%sumtails",op%sumtails
903 0 : tail = op%sumtails
904 0 : new_tail = tail + 1
905 0 : op%tails(aF)= op%tails(aF) + 1
906 0 : DO indice = aF +1, op%flavors+1
907 0 : op%Fshift(indice) = op%Fshift(indice) + 1
908 : END DO
909 0 : op%sumtails = op%Fshift(op%flavors) + op%tails(op%flavors) !last slot of Fshift is the tail of full matrix
910 : !write(6,*) " BathOperatoroffdiag_setMAdd op%sumtails",op%sumtails
911 : !write(6,*) " setMAdd actualized Fshift",(op%Fshift(it),it=1,op%flavors+1)
912 : !write(6,*) " setMAdd actualized tails",(op%tails(it),it=1,op%flavors)
913 : !CALL matrix_print(M)
914 :
915 0 : if(op%opt_nondiag==1) then
916 : tailb = 1
917 : taile = tail
918 : else
919 : !sui!write(6,*) "Bathoperator a opt_nondiag=0"
920 0 : tailb = op%Fshift(aF)+1
921 0 : taile = op%Fshift(aF)+op%tails(aF)
922 : endif
923 :
924 : ! --- data obtained from BathOperatoroffdiag_getDetAdd
925 0 : PositionRow = op%updatePosRow + op%Fshift(aF) ! position in the full matrix
926 0 : PositionCol = op%updatePosCol + op%Fshift(aF) ! position in the full matrix
927 0 : Stilde = op%Stilde
928 :
929 : ! !write(6,*) "before", positionRow, positionCol
930 : !CALL MatrixHyb_print(op%M(aF),opt_print=1)
931 : ! --- MatrixHyb_setSize
932 : !write(6,*) " BathOperatoroffdiag_setMAdd before setsize",size(op%M%mat,1)
933 0 : CALL MatrixHyb_setSize(op%M,new_tail)
934 : !write(6,*) " BathOperatoroffdiag_setMAdd after setsize",size(op%M%mat,1)
935 :
936 : ! Compute Qtilde with Q
937 : !op%Q%vec(1:tail) = (-1.d0) * MATMUL(op%M(aF)%mat(1:tail,1:tail),op%Q%vec(1:tail)) * Stilde
938 :
939 : ! --- M*Q => Q
940 0 : op%Q%vec(tailb:taile) = MATMUL(op%M%mat(tailb:taile,tailb:taile),op%Q%vec(tailb:taile))
941 :
942 : !op%Q%vec(PositionRow:new_tail) = EOSHIFT(op%Q%vec(PositionRow:new_tail), SHIFT=-1, BOUNDARY=-1.d0, DIM=1)
943 : ! op%Qtau%vec(PositionCol:new_tail) = EOSHIFT(op%Qtau%vec(PositionCol:new_tail), SHIFT=-1, BOUNDARY=1.d0, DIM=1)
944 : ! op%Qtau%vec(PositionCol) = op%Stau
945 :
946 : !Compute Rtilde with R and without multiplying by Stilde
947 : !op%R%vec(1:tail) = (-1.d0) * MATMUL(op%R%vec(1:tail),op%M(aF)%mat(1:tail,1:tail))
948 :
949 : ! --- R*M => R
950 0 : op%R%vec(tailb:taile) = MATMUL(op%R%vec(tailb:taile),op%M%mat(tailb:taile,tailb:taile))
951 :
952 : !op%R%vec(PositionCol:new_tail) = EOSHIFT(op%R%vec(PositionCol:new_tail), SHIFT=-1, BOUNDARY=-1.d0, DIM=1)
953 : ! op%Rtau%vec(PositionRow:new_tail) = EOSHIFT(op%Rtau%vec(PositionRow:new_tail), SHIFT=-1, BOUNDARY=1.d0, DIM=1)
954 : ! op%Rtau%vec(PositionRow) = op%Stau
955 :
956 : !Compute the new M matrix
957 : !op%M(aF)%mat(PositionRow:new_tail,1:new_tail) = &
958 : ! EOSHIFT(op%M(aF)%mat(PositionRow:new_tail,1:new_tail),SHIFT=-1, BOUNDARY=0.d0, DIM=1)
959 : !op%M(aF)%mat(1:n12 characters (ABI_MALLOC) instead of 4 (FREE)ew_tail,PositionCol:new_tail) = &
960 : ! EOSHIFT(op%M(aF)%mat(1:new_tail,PositionCol:new_tail),SHIFT=-1, BOUNDARY=0.d0, DIM=2)
961 : ! ! op%M(aF)%mat(1:new_tail,1:new_tail) = op%M(aF)%mat(1:new_tail,1:new_tail) + &
962 : ! ! Stilde * MATMUL(RESHAPE(op%Q%vec(1:new_tail),(/ new_tail,1 /)),RESHAPE(op%R%vec(1:new_tail),(/ 1,new_tail /)))
963 :
964 : !op%M(aF)%mat_tau(PositionRow:new_tail,1:new_tail) = &
965 : ! EOSHIFT(op%M(aF)%mat_tau(PositionRow:new_tail,1:new_tail),SHIFT=-1, BOUNDARY=0, DIM=1)
966 : !op%M(aF)%mat_tau(1:new_tail,PositionCol:new_tail) = &
967 : ! EOSHIFT(op%M(aF)%mat_tau(1:new_tail,PositionCol:new_tail),SHIFT=-1, BOUNDARY=0, DIM=2)
968 :
969 0 : mbeta_two = -op%beta*0.5d0
970 0 : inv_dt = op%inv_dt
971 :
972 : ! ------ Shift mat_tau and update old M=Ptilde
973 0 : DO col=tail,1,-1 ! decreasing order to avoid overwrite of data
974 0 : col_move = col + ( 1+SIGN(1,col-PositionCol) )/2
975 : ! if col>= PositionCol col_move=col+1
976 : ! if col< PositionCol col_move=col
977 0 : DO row=tail,1,-1
978 0 : row_move = row + ( 1+SIGN(1,row-PositionRow) )/2
979 : ! --- times for Ptilde are kept unchanged. But we have to copy it at the right place
980 : op%M%mat_tau(row_move,col_move) = &
981 0 : op%M%mat_tau(row,col)
982 : ! --- Update Ptilde with the same indices as mat_tau
983 : ! --- M + M*Q Stilde R*M => Ptilde => M
984 : !if(row>=tailb.and.row<=taile.and.col>=tailb.and.col<=taile) then
985 : op%M%mat(row_move,col_move) = &
986 0 : op%M%mat(row,col) + op%Q%vec(row)*op%R%vec(col) * Stilde
987 : !else
988 : ! op%M%mat(row_move,col_move) = op%M%mat(row,col)
989 : !endif
990 : END DO
991 : END DO
992 :
993 : ! ------ Add new stuff for new row
994 0 : DO row = 1, tail
995 0 : row_move = row + ( 1+SIGN(1,row-PositionRow) )/2
996 : ! --- M*Q Stilde => Qtilde => M with the good indices
997 : !if(row>=tailb.and.row<=taile) then
998 0 : op%M%mat(row_move,PositionCol) = -op%Q%vec(row)*Stilde
999 : !else
1000 : ! op%M%mat(row_move,PositionCol) = op%M%mat(row,PositionCol)
1001 : !endif
1002 :
1003 0 : time = op%Rtau%vec(row) ! pourquoi Rtau et pas Qtau ici ?
1004 0 : time = time + ( SIGN(1.d0,time) - 1.d0 )*mbeta_two
1005 : ! if time>=0 time=time
1006 : ! if time< 0 time=time + beta
1007 : ! --- mat_tau=int(time*L/beta+1.5)
1008 0 : op%M%mat_tau(row,PositionCol) = INT ( (time*inv_dt) +1.5d0 )
1009 : !write(6,*) " setMadd new row", op%Rtau%vec(row),op%M%mat_tau(row,PositionCol)
1010 : ! if(op%M%mat_tau(row,PositionCol)>301) then
1011 : ! !write(6,*) ">301 a", time,inv_dt, op%M%mat_tau(row,PositionCol)
1012 : ! time = op%Rtau%vec(row) ! pourquoi Rtau et pas Qtau ici ?
1013 : ! !write(6,*) time,mbeta_two
1014 : ! time = time + ( SIGN(1.d0,time) - 1.d0 )*mbeta_two
1015 : ! !write(6,*) time
1016 : ! !write(6,*) INT ( (time*inv_dt) +1.5d0 )
1017 : ! stop
1018 : ! endif
1019 : END DO
1020 : ! Add last time missing in the loops
1021 0 : time = op%Rtau%vec(new_tail)
1022 0 : time = time + ( SIGN(1.d0,time) - 1.d0 )*mbeta_two
1023 0 : op%M%mat_tau(new_tail,PositionCol) = INT ( (time*inv_dt) +1.5d0 )
1024 : !write(6,*) " setMadd last time", op%Rtau%vec(new_tail),op%M%mat_tau(new_tail,PositionCol)
1025 : ! if(op%M%mat_tau(new_tail,PositionCol)>301) then
1026 : ! !write(6,*) ">301 b", time,inv_dt, op%M%mat_tau(new_tail,PositionCol)
1027 : ! time = op%Rtau%vec(new_tail) ! pourquoi Rtau et pas Qtau ici ?
1028 : ! !write(6,*) time,mbeta_two
1029 : ! time = time + ( SIGN(1.d0,time) - 1.d0 )*mbeta_two
1030 : ! !write(6,*) time
1031 : ! !write(6,*) INT ( (time*inv_dt) +1.5d0 )
1032 : ! stop
1033 : ! endif
1034 :
1035 : ! Add new stuff for new col
1036 0 : DO col = 1, tail
1037 0 : col_move = col + ( 1+SIGN(1,col-PositionCol) )/2
1038 : ! --- Stilde RN => Rtilde => M
1039 : !if(col>=tailb.and.col<=taile) then
1040 0 : op%M%mat(PositionRow,col_move) = -op%R%vec(col)*Stilde
1041 : !else
1042 : ! op%M%mat(PositionRow,col_move) = op%M%mat(PositionRow,col)
1043 : !endif
1044 0 : time = op%Qtau%vec(col)
1045 0 : time = time + ( SIGN(1.d0,time) - 1.d0 )*mbeta_two
1046 0 : op%M%mat_tau(PositionRow,col) = INT ( (time*inv_dt) +1.5d0 )
1047 : !write(6,*) " setMadd new col", op%Qtau%vec(col),op%M%mat_tau(PositionRow,col)
1048 : ! if(op%M%mat_tau(PositionRow,col)>301) then
1049 : ! !write(6,*) ">301 c", time,inv_dt, op%M%mat_tau(PositionRow,col)
1050 : ! time = op%Qtau%vec(col) ! pourquoi Rtau et pas Qtau ici ?
1051 : ! !write(6,*) time,mbeta_two
1052 : ! time = time + ( SIGN(1.d0,time) - 1.d0 )*mbeta_two
1053 : ! !write(6,*) time
1054 : ! !write(6,*) INT ( (time*inv_dt) +1.5d0 )
1055 : ! stop
1056 : ! endif
1057 : END DO
1058 : ! Add last time missing in the loops
1059 0 : time = op%Qtau%vec(new_tail)
1060 0 : time = time + ( SIGN(1.d0,time) - 1.d0 )*mbeta_two
1061 0 : op%M%mat_tau(PositionRow,new_tail) = INT ( (time*inv_dt) +1.5d0 )
1062 : !write(6,*) " setMadd last time", op%Qtau%vec(new_tail),op%M%mat_tau(PositionRow,new_tail)
1063 : ! if(op%M%mat_tau(PositionRow,new_tail)>301) then
1064 : ! !write(6,*) ">301 d", time,inv_dt, op%M%mat_tau(PositionRow,new_tail)
1065 : ! time = op%Qtau%vec(new_tail) ! pourquoi Rtau et pas Qtau ici ?
1066 : ! !write(6,*) time,mbeta_two
1067 : ! time = time + ( SIGN(1.d0,time) - 1.d0 )*mbeta_two
1068 : ! !write(6,*) time
1069 : ! !write(6,*) INT ( (time*inv_dt) +1.5d0 )
1070 : ! stop
1071 : ! endif
1072 :
1073 0 : op%M%mat(PositionRow,PositionCol) = Stilde
1074 :
1075 : !CALL MatrixHyb_print(op%M,opt_print=1)
1076 :
1077 : ! DO col = 1, new_tail
1078 : ! time = op%Rtau%vec(col)
1079 : ! time = time + ( SIGN(1.d0,time) - 1.d0 )*mbeta_two
1080 : ! op%M(aF)%mat_tau(col,PositionCol) = INT ( (time*inv_dt) +1.5d0 )
1081 : ! time = op%Qtau%vec(col)
1082 : ! time = time + ( SIGN(1.d0,time) - 1.d0 )*mbeta_two
1083 : ! op%M(aF)%mat_tau(PositionRow,Col) = INT ( (time*inv_dt) +1.5d0 )
1084 : ! time = op%R%vec(col)*Stilde
1085 : ! DO row = 1, new_tail
1086 : ! op%M(aF)%mat(row,col) = op%M(aF)%mat(row,col) + op%Q%vec(row)*time
1087 : ! END DO
1088 : ! END DO
1089 :
1090 : !col_move = new_tail
1091 : !col = tail
1092 : !DO col_move = new_tail, 1, -1
1093 : ! IF ( col_move .EQ. positionCol ) THEN
1094 : ! ! on calcule rajoute Q tilde
1095 : ! !row_move = new_tail
1096 : ! row = tail
1097 : ! DO row_move = new_tail, 1, -1
1098 : ! ! calcul itau
1099 : ! IF ( row_move .EQ. positionRow ) THEN
1100 : ! op%M(aF)%mat(row_move,col_move) = Stilde
1101 : ! !time = op%Stau
1102 : ! ELSE
1103 : ! op%M(aF)%mat(row_move,col_move) = -op%Q%vec(row)*Stilde
1104 : ! !time = op%Rtau%vec(row_move)
1105 : ! row = row - 1
1106 : ! END IF
1107 : ! !time = time + ( SIGN(1.d0,time) - 1.d0 )*mbeta_two
1108 : ! !op%M(aF)%mat_tau(row_move,col_move) = INT ( (time*inv_dt) +1.5d0 )
1109 : ! END DO
1110 : ! ! realignement des indices
1111 : ! ELSE
1112 : ! ! on calcule Ptilde
1113 : ! !row_move = new_tail
1114 : ! row = tail
1115 : ! DO row_move = new_tail, 1, -1
1116 : ! IF ( row_move .EQ. positionRow ) THEN
1117 : ! op%M(aF)%mat(row_move,col_move) = -op%R%vec(col) * Stilde
1118 : ! ! calcul itau
1119 : ! !time = op%Qtau%vec(col_move)
1120 : ! !time = time + ( SIGN(1.d0,time) - 1.d0 )*mbeta_two
1121 : ! !op%M(aF)%mat_tau(row_move,col_move) = INT ( (time*inv_dt) +1.5d0 )
1122 : ! ELSE
1123 : ! op%M(aF)%mat(row_move,col_move) = op%M(aF)%mat(row,col) + op%Q%vec(row)*op%R%vec(col)*Stilde
1124 : ! ! copy itau
1125 : ! !op%M(aF)%mat_tau(row_move,col_move) = op%M(aF)%mat_tau(row,col)
1126 : ! row = row - 1
1127 : ! END IF
1128 : ! END DO
1129 : ! col = col - 1
1130 : ! END IF
1131 : !END DO
1132 : ! !write(6,*) "after"
1133 : ! CALL MatrixHyb_print(op%M(aF),opt_print=1)
1134 : !CALL matrix_inverse(M)
1135 : !CALL MatrixHyb_print(M)
1136 : !CALL matrix_inverse(M)
1137 :
1138 0 : IF ( op%antiShift .EQV. .TRUE. ) THEN ! antisegment
1139 : if(3==4) then
1140 : CALL Vector_init(vec_tmp,new_tail)
1141 : CALL VectorInt_init(vecI_tmp,new_tail)
1142 : ! Shift if necessary according to op%antishift
1143 : ! shift DIM=2 (col)
1144 :
1145 : ! For new_tail=4, the following lines transform
1146 : ! M=(a,b,c,d) vith a,b,c,d column vectors into
1147 : ! M=(d,a,b,c)
1148 : p = new_tail - 1 ! = tail
1149 : m = 1
1150 : ! count increases in the loop from 0 to new_tail-1
1151 : count = 0
1152 : DO WHILE ( count .NE. new_tail )
1153 : ! put column b in vec_tmp
1154 : vec_tmp%vec(1:new_tail) = op%M%mat(1:new_tail,m)
1155 : vecI_tmp%vec(1:new_tail) = op%M%mat_tau(1:new_tail,m)
1156 : i = m
1157 : !j = m+p
1158 : MODCYCLE(m+p, new_tail, j) ! j=m+p modulo new_tail
1159 : DO WHILE (j .NE. m)
1160 : op%M%mat(1:new_tail,i) = op%M%mat(1:new_tail,j)
1161 : op%M%mat_tau(1:new_tail,i) = op%M%mat_tau(1:new_tail,j)
1162 : i = j
1163 : MODCYCLE(j+p, new_tail, j)
1164 : count = count+1
1165 : END DO
1166 : op%M%mat(1:new_tail,i) = vec_tmp%vec(1:new_tail)
1167 : op%M%mat_tau(1:new_tail,i) = vecI_tmp%vec(1:new_tail)
1168 : count = count+1
1169 : m = m+1
1170 : END DO
1171 : ! shift DIM=1 (row)
1172 :
1173 : ! below is similar to above but for rows instead of columns.
1174 : p = new_tail - 1
1175 : m = 1
1176 : count = 0
1177 : DO WHILE ( count .NE. new_tail)
1178 : vec_tmp%vec(1:new_tail) = op%M%mat(m,1:new_tail)
1179 : vecI_tmp%vec(1:new_tail) = op%M%mat_tau(m,1:new_tail)
1180 : i = m
1181 : !j = m+p
1182 : MODCYCLE(m+p, new_tail, j)
1183 : DO WHILE ( j .NE. m )
1184 : op%M%mat(i,1:new_tail) = op%M%mat(j,1:new_tail)
1185 : op%M%mat_tau(i,1:new_tail) = op%M%mat_tau(j,1:new_tail)
1186 : i = j
1187 : MODCYCLE(j+p, new_tail, j)
1188 : count = count+1
1189 : END DO
1190 : op%M%mat(i,1:new_tail) = vec_tmp%vec(1:new_tail)
1191 : op%M%mat_tau(i,1:new_tail) = vecI_tmp%vec(1:new_tail)
1192 : count = count+1
1193 : m = m+1
1194 : END DO
1195 : CALL Vector_destroy(vec_tmp)
1196 : CALL VectorInt_destroy(vecI_tmp)
1197 : endif
1198 : !op%M(aF)%mat(1:new_tail,1:new_tail) = CSHIFT(op%M(aF)%mat(1:new_tail,1:new_tail), SHIFT=-1, DIM=1) ! Shift to the bottom
1199 : !op%M(aF)%mat(1:new_tail,1:new_tail) = CSHIFT(op%M(aF)%mat(1:new_tail,1:new_tail), SHIFT=-1, DIM=2) ! Shift to the right
1200 : !op%M(aF)%mat_tau(1:new_tail,1:new_tail) = CSHIFT(op%M(aF)%mat_tau(1:new_tail,1:new_tail), SHIFT=-1, DIM=1) ! Shift to the bottom
1201 : !op%M(aF)%mat_tau(1:new_tail,1:new_tail) = CSHIFT(op%M(aF)%mat_tau(1:new_tail,1:new_tail), SHIFT=-1, DIM=2) ! Shift to the right
1202 : !write(6,*) " setMAdd size M%mat",size(op%M%mat,1),size(op%M%mat,2),new_tail
1203 : !write(6,*) " setMAdd arguement M%mat",aF,op%Fshift(aF)
1204 : !write(6,*) " setMAdd arguement M%mat",aF,op%Fshift(aF),op%Fshift(aF+1)
1205 0 : do it=1,op%sumtails
1206 : !write(6,*) " setMAdd before antishift M%mat_tau",(op%M%mat_tau(it,it1),it1=1,op%sumtails)
1207 : enddo
1208 0 : if (new_tail>0.and.op%Fshift(aF+1)>op%Fshift(aF)) then
1209 : op%M%mat(op%Fshift(aF)+1:op%Fshift(aF+1) , 1:new_tail) = &
1210 0 : CSHIFT( op%M%mat(op%Fshift(aF)+1:op%Fshift(aF+1) , 1:new_tail) , SHIFT=-1 , DIM=1) ! Shift to the bottom
1211 :
1212 : op%M%mat(1:new_tail , op%Fshift(aF)+1:op%Fshift(aF+1)) = &
1213 0 : CSHIFT( op%M%mat(1:new_tail , op%Fshift(aF)+1:op%Fshift(aF+1)) , SHIFT=-1 , DIM=2) ! Shift to the right
1214 :
1215 : op%M%mat_tau(op%Fshift(aF)+1:op%Fshift(aF+1) , 1:new_tail) = &
1216 0 : CSHIFT( op%M%mat_tau(op%Fshift(aF)+1:op%Fshift(aF+1) , 1:new_tail) , SHIFT=-1 , DIM=1) ! Shift to the bottom
1217 :
1218 : op%M%mat_tau(1:new_tail , op%Fshift(aF)+1:op%Fshift(aF+1)) = &
1219 0 : CSHIFT( op%M%mat_tau(1:new_tail , op%Fshift(aF)+1:op%Fshift(aF+1)) , SHIFT=-1 , DIM=2) ! Shift to the right
1220 : end if
1221 : !CALL matrix_print(M)
1222 : END IF
1223 :
1224 : !! do it=1,op%sumtails
1225 : !write(6,*) " setMAdd end M",(op%M%mat(it,it1),it1=1,op%sumtails)
1226 : !! enddo
1227 : !! do it=1,op%sumtails
1228 : !! !write(6,*) " setMAdd end M%mat_tau",(op%M%mat_tau(it,it1),it1=1,op%sumtails)
1229 : !! enddo
1230 0 : IF ( op%doCheck .EQV. .TRUE.) THEN
1231 : !#ifdef CTQMC_CHECK
1232 0 : CALL BathOperatoroffdiag_checkM(op,particle)
1233 : !#endif
1234 : END IF
1235 :
1236 0 : op%MAddFlag = .FALSE.
1237 :
1238 0 : END SUBROUTINE BathOperatoroffdiag_setMAdd
1239 : !!***
1240 :
1241 : !!****f* ABINIT/m_BathOperatoroffdiag/BathOperatoroffdiag_setMRemove
1242 : !! NAME
1243 : !! BathOperatoroffdiag_setMRemove
1244 : !!
1245 : !! FUNCTION
1246 : !! delete one row and one column of the M matrix
1247 : !!
1248 : !! COPYRIGHT
1249 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
1250 : !! This file is distributed under the terms of the
1251 : !! GNU General Public License, see ~abinit/COPYING
1252 : !! or http://www.gnu.org/copyleft/gpl.txt .
1253 : !!
1254 : !! INPUTS
1255 : !! op=bath operator
1256 : !! particle=segments of the active flavor
1257 : !!
1258 : !! OUTPUT
1259 : !!
1260 : !! SIDE EFFECTS
1261 : !!
1262 : !! NOTES
1263 : !!
1264 : !! SOURCE
1265 :
1266 0 : SUBROUTINE BathOperatoroffdiag_setMRemove(op,particle)
1267 :
1268 : !Arguments ------------------------------------
1269 : TYPE(BathOperatoroffdiag), INTENT(INOUT) :: op
1270 : TYPE(ListCdagC) , INTENT(IN ) :: particle(:)
1271 : !Local variables ------------------------------
1272 : INTEGER :: tail,tailb,taile
1273 : INTEGER :: new_tail
1274 : INTEGER :: col
1275 : INTEGER :: col_move
1276 : INTEGER :: row_move
1277 : INTEGER :: row
1278 : INTEGER :: positionCol
1279 : INTEGER :: positionRow
1280 : INTEGER :: aF,iaf
1281 : INTEGER :: m
1282 : INTEGER :: count
1283 : INTEGER :: i
1284 : INTEGER :: j,it !,it1
1285 : INTEGER :: p
1286 : DOUBLE PRECISION :: invStilde
1287 : DOUBLE PRECISION :: invStilde2
1288 0 : TYPE(VectorInt) :: vecI_tmp
1289 0 : TYPE(Vector) :: vec_tmp
1290 :
1291 0 : IF ( op%MRemoveFlag .EQV. .FALSE. ) &
1292 0 : CALL ERROR("BathOperatoroffdiag_setMRemove : MRemoveFlag turn off ")
1293 0 : aF = op%activeFlavor
1294 0 : IF ( aF .LE. 0 ) &
1295 0 : CALL ERROR("BathOperatoroffdiag_setMRemove : no active hybrid func ")
1296 0 : do it=1,op%sumtails
1297 : !write(6,*) " setMRemove begin M",(op%M%mat(it,it1),it1=1,op%sumtails)
1298 : enddo
1299 0 : tail = op%sumtails
1300 0 : new_tail = tail - 1
1301 0 : op%tails(af)= op%tails(af) - 1
1302 0 : DO iaf=af+1 , op%flavors+1
1303 0 : op%Fshift(iaf) = op%Fshift(iaf) - 1
1304 : END DO
1305 0 : op%sumtails = op%Fshift(op%flavors) + op%tails(op%flavors)
1306 0 : positionCol = op%updatePosCol + op%Fshift(af)
1307 0 : positionRow = op%updatePosRow + op%Fshift(af)
1308 0 : invStilde = 1.d0 / op%Stilde
1309 : if(op%opt_nondiag==1) then
1310 : tailb = 1
1311 : taile = new_tail
1312 : else
1313 : !sui!write(6,*) "Bathoperator c opt_nondiag=0"
1314 : tailb = op%Fshift(aF)+1
1315 : taile = op%Fshift(aF)+op%tails(aF)
1316 : endif
1317 :
1318 : ! !write(6,*) "before", positionRow, positionCol
1319 : ! CALL MatrixHyb_print(op%M(aF),opt_print=1)
1320 :
1321 : ! IF ( new_tail .EQ. 0 ) THEN
1322 : !! IF ( op%antiShift .EQV. .TRUE. ) THEN
1323 : !! op%M(aF)%mat(1,1) = 1.d0/BathOperatoroffdiag_Hybrid(op, op%beta)
1324 : !! op%MRemoveFlag = .FALSE.
1325 : !! RETURN
1326 : !! END IF
1327 : ! CALL MatrixHyb_clear(op%M(aF))
1328 : ! op%MRemoveFlag = .FALSE.
1329 : ! RETURN
1330 : ! END IF
1331 :
1332 : ! CALL Vector_setSize(op%Q,new_tail)
1333 : ! CALL Vector_setSize(op%R,new_tail)
1334 0 : Vector_QuickResize(op%Q,new_tail)
1335 0 : Vector_QuickResize(op%R,new_tail)
1336 :
1337 : ! We use R and Q as op%R%vec and op%Q%vec
1338 : ! op%R%vec => op%R
1339 : ! op%Q%vec => op%Q
1340 :
1341 0 : row = 1
1342 : !row_move = 1
1343 0 : col = 1
1344 : !col_move = 1
1345 0 : DO row_move = 1, new_tail
1346 0 : IF ( row .EQ. positionRow ) row = row + 1
1347 0 : IF ( col .EQ. positionCol ) col = col + 1
1348 : !col = row_move + (1+SIGN(1,row_move-positionCol))/2
1349 : !row = row_move + (1+SIGN(1,row_move-positionRow))/2
1350 0 : op%R%vec(row_move) = op%M%mat(positionRow,col)
1351 0 : op%Q%vec(row_move) = op%M%mat(row,positionCol)
1352 0 : row = row + 1
1353 0 : col = col + 1
1354 : END DO
1355 : !! op%R%vec(1:positionCol-1) = op%M(aF)%mat(positionRow,1:positionCol-1)
1356 : !! op%R%vec(positionCol:new_tail) = op%M(aF)%mat(positionRow,positionCol+1:tail)
1357 : !! op%Q%vec(1:positionRow-1) = op%M(aF)%mat(1:positionRow-1,positionCol)
1358 : !! op%Q%vec(positionRow:new_tail) = op%M(aF)%mat(positionRow+1:tail,positionCol)
1359 : !write(*,*) positionRow, positionCol
1360 : !CALL MatrixHyb_print(M)
1361 : !CALL Vector_print(op%R)
1362 : !CALL Vector_print(op%Q)
1363 : !CALL ListCdagC_print(op%ListCdagC)
1364 :
1365 : col = 1
1366 0 : DO col_move = 1, new_tail
1367 0 : IF ( col_move .EQ. positionCol ) col = col + 1
1368 : !col = col_move + (1+SIGN(1,col_move-positionCol))/2
1369 0 : row = 1
1370 0 : invStilde2 = invStilde * op%R%vec(col_move)
1371 0 : DO row_move = 1, new_tail
1372 0 : IF ( row_move .EQ. positionRow ) row = row + 1
1373 : !row = row_move + (1+SIGN(1,row_move-positionRow))/2
1374 : ! Compute for all rows and cols M <= M - Q 1/S R
1375 : !if(row_move>=tailb.and.row_move<=taile.and.col_move>=tailb.and.col_move<=taile) then
1376 : op%M%mat(row_move,col_move) = op%M%mat(row,col) &
1377 0 : - op%Q%vec(row_move)*invStilde2
1378 : !else
1379 : ! op%M%mat(row_move,col_move) = op%M%mat(row,col)
1380 : !endif
1381 0 : op%M%mat_tau(row_move,col_move) = op%M%mat_tau(row,col)
1382 0 : row = row + 1
1383 : END DO
1384 0 : col = col + 1
1385 : END DO
1386 0 : CALL MatrixHyb_setSize(op%M,new_tail)
1387 :
1388 0 : IF ( op%antiShift .EQV. .TRUE. ) THEN ! antisegment
1389 : if(3==4) then
1390 : ! Shift if necessary according to op%antishift
1391 : ! shift DIM=2 (col)
1392 : CALL Vector_init(vec_tmp,new_tail)
1393 : CALL VectorInt_init(vecI_tmp,new_tail)
1394 : p = 1
1395 : m = 1
1396 : count = 0
1397 : DO WHILE ( count .NE. new_tail )
1398 : vec_tmp%vec(1:new_tail) = op%M%mat(1:new_tail,m)
1399 : vecI_tmp%vec(1:new_tail) = op%M%mat_tau(1:new_tail,m)
1400 : i = m
1401 : !j = m+p
1402 : MODCYCLE(m+p, new_tail, j)
1403 : DO WHILE (j .NE. m)
1404 : op%M%mat(1:new_tail,i) = op%M%mat(1:new_tail,j)
1405 : op%M%mat_tau(1:new_tail,i) = op%M%mat_tau(1:new_tail,j)
1406 : i = j
1407 : MODCYCLE(j+p, new_tail, j)
1408 : count = count+1
1409 : END DO
1410 : op%M%mat(1:new_tail,i) = vec_tmp%vec(1:new_tail)
1411 : op%M%mat_tau(1:new_tail,i) = vecI_tmp%vec(1:new_tail)
1412 : count = count+1
1413 : m = m+1
1414 : END DO
1415 : CALL Vector_destroy(vec_tmp)
1416 : CALL VectorInt_destroy(vecI_tmp)
1417 : !op%M(aF)%mat(1:new_tail,1:new_tail) = &
1418 : ! CSHIFT(op%M(aF)%mat(1:new_tail,1:new_tail), SHIFT=1, DIM=2) ! Shift to the top
1419 : !op%M(aF)%mat_tau(1:new_tail,1:new_tail) = &
1420 : ! CSHIFT(op%M(aF)%mat_tau(1:new_tail,1:new_tail), SHIFT=1, DIM=2) ! Shift to the top
1421 : endif
1422 0 : if (new_tail>0.and.op%Fshift(af+1)>op%Fshift(af)) then
1423 : op%M%mat(1:new_tail,op%Fshift(af)+1:op%Fshift(af+1)) = &
1424 0 : CSHIFT(op%M%mat(1:new_tail,op%Fshift(af)+1:op%Fshift(af+1)), SHIFT=1, DIM=2) ! Shift to the top
1425 : op%M%mat_tau(1:new_tail,op%Fshift(af)+1:op%Fshift(af+1)) = &
1426 0 : CSHIFT(op%M%mat_tau(1:new_tail,op%Fshift(af)+1:op%Fshift(af+1)), SHIFT=1, DIM=2) ! Shift to the top
1427 : end if
1428 : END IF
1429 : ! !write(6,*) "after "
1430 : ! CALL MatrixHyb_print(op%M(aF),opt_print=1)
1431 :
1432 0 : IF ( op%doCheck .EQV. .TRUE. ) THEN
1433 : !#ifdef CTQMC_CHECK
1434 0 : CALL BathOperatoroffdiag_checkM(op,particle)
1435 : !#endif
1436 : END IF
1437 0 : do it=1,op%sumtails
1438 : !write(6,*) " setMRemove end M",(op%M%mat(it,it1),it1=1,op%sumtails)
1439 : enddo
1440 :
1441 0 : op%MRemoveFlag = .FALSE.
1442 :
1443 0 : END SUBROUTINE BathOperatoroffdiag_setMRemove
1444 : !!***
1445 :
1446 : !!****f* ABINIT/m_BathOperatoroffdiag/BathOperatoroffdiag_swap
1447 : !! NAME
1448 : !! BathOperatoroffdiag_swap
1449 : !!
1450 : !! FUNCTION
1451 : !! Recompute 2 M matrix swaping the segments (used for Global moves)
1452 : !!
1453 : !! COPYRIGHT
1454 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
1455 : !! This file is distributed under the terms of the
1456 : !! GNU General Public License, see ~abinit/COPYING
1457 : !! or http://www.gnu.org/copyleft/gpl.txt .
1458 : !!
1459 : !! INPUTS
1460 : !! op=bath operator
1461 : !! iflavor1=flavor to swap with the next one
1462 : !! iflavor2=favor to swap with the previous one
1463 : !!
1464 : !! OUTPUT
1465 : !!
1466 : !! SIDE EFFECTS
1467 : !!
1468 : !! NOTES
1469 : !!
1470 : !! SOURCE
1471 :
1472 0 : SUBROUTINE BathOperatoroffdiag_swap(op, flavor1, flavor2)
1473 :
1474 : !Arguments ------------------------------------
1475 : TYPE(BathOperatoroffdiag), INTENT(INOUT) :: op
1476 : INTEGER , INTENT(IN ) :: flavor1
1477 : INTEGER , INTENT(IN ) :: flavor2
1478 : INTEGER :: ii,iflavort,itmptail,flavora,flavorb
1479 0 : DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: mat_temp
1480 0 : INTEGER , ALLOCATABLE, DIMENSION(:,:) :: mat_tau_temp
1481 :
1482 0 : if(flavor1>flavor2) then
1483 : flavora=flavor2
1484 : flavorb=flavor1
1485 : else
1486 0 : flavora=flavor1
1487 0 : flavorb=flavor2
1488 : endif
1489 0 : MALLOC(mat_temp,(1:op%sumtails,1:op%sumtails))
1490 0 : MALLOC(mat_tau_temp,(1:op%sumtails,1:op%sumtails))
1491 : !mat_temp= op%M%mat
1492 : !mat_tau_temp= op%M%mat_tau
1493 : !it1=0
1494 : !do iflav1=1,op%flavors
1495 : ! do ii1=1,op%tails(iflav1)
1496 : ! it1=it1+1
1497 : ! it2=0
1498 : ! do iflav2=1,op%flavors
1499 : ! do ii2=1,op%tails(iflav1)
1500 : ! it2=it2+1
1501 : ! op%M%mat(it1,it2)=
1502 : ! enddo
1503 : ! enddo
1504 : ! enddo
1505 : !enddo
1506 : if(3==3) then
1507 0 : op%M=op%M_update
1508 0 : if (op%sumtails>0) then
1509 : ! shift block flavorb at the place of flavora (column)
1510 0 : if (op%Fshift(flavorb+1)>op%Fshift(flavora)) then
1511 0 : do ii=1, op%tails(flavorb)
1512 : op%M%mat(op%Fshift(flavora)+1:op%Fshift(flavorb+1) , 1:op%sumtails) = &
1513 0 : CSHIFT( op%M%mat(op%Fshift(flavora)+1:op%Fshift(flavorb+1) , 1:op%sumtails) , SHIFT=-1 , DIM=1)
1514 : op%M%mat_tau(op%Fshift(flavora)+1:op%Fshift(flavorb+1) , 1:op%sumtails) = &
1515 0 : CSHIFT( op%M%mat_tau(op%Fshift(flavora)+1:op%Fshift(flavorb+1) , 1:op%sumtails) , SHIFT=-1 , DIM=1)
1516 : enddo
1517 : end if
1518 :
1519 : ! shift block flavora at the place of flavorb (column)
1520 0 : if (op%Fshift(flavorb)>op%Fshift(flavora)) then
1521 0 : do ii=1, op%tails(flavora)
1522 : op%M%mat(op%Fshift(flavora)+op%tails(flavorb)+&
1523 : & 1:op%Fshift(flavorb)+op%tails(flavorb) , 1:op%sumtails) = &
1524 : CSHIFT( op%M%mat( op%Fshift(flavora)+op%tails(flavorb)&
1525 0 : & +1:op%Fshift(flavorb)+op%tails(flavorb) , 1:op%sumtails) , SHIFT=1 , DIM=1)
1526 : op%M%mat_tau(op%Fshift(flavora)+op%tails(flavorb)+1:op%Fshift(flavorb)+&
1527 : & op%tails(flavorb) , 1:op%sumtails) = &
1528 : CSHIFT( op%M%mat_tau( op%Fshift(flavora)+op%tails(flavorb)+&
1529 0 : & 1:op%Fshift(flavorb)+op%tails(flavorb) , 1:op%sumtails) , SHIFT=1 , DIM=1)
1530 : enddo
1531 : end if
1532 :
1533 : ! shift block flavorb at the place of flavora (row)
1534 0 : if (op%Fshift(flavorb+1)>op%Fshift(flavora)) then
1535 0 : do ii=1, op%tails(flavorb)
1536 : op%M%mat(1:op%sumtails , op%Fshift(flavora)+1:op%Fshift(flavorb+1)) = &
1537 0 : CSHIFT( op%M%mat(1:op%sumtails , op%Fshift(flavora)+1:op%Fshift(flavorb+1)) , SHIFT=-1 , DIM=2)
1538 : op%M%mat_tau(1:op%sumtails , op%Fshift(flavora)+1:op%Fshift(flavorb+1)) = &
1539 0 : CSHIFT( op%M%mat_tau(1:op%sumtails , op%Fshift(flavora)+1:op%Fshift(flavorb+1)) , SHIFT=-1 , DIM=2)
1540 : enddo
1541 : end if
1542 :
1543 : ! shift block flavora at the place of flavorb (row)
1544 0 : if (op%Fshift(flavorb)>op%Fshift(flavora)) then
1545 0 : do ii=1, op%tails(flavora)
1546 : op%M%mat(1:op%sumtails , op%Fshift(flavora)+op%tails(flavorb)+1:op%Fshift(flavorb)+op%tails(flavorb)) = &
1547 : CSHIFT( op%M%mat(1:op%sumtails ,op%Fshift(flavora)+op%tails(flavorb)&
1548 0 : & +1:op%Fshift(flavorb)+op%tails(flavorb)) , SHIFT=1 , DIM=2)
1549 : op%M%mat_tau(1:op%sumtails ,op%Fshift(flavora)+op%tails(flavorb)+&
1550 : & 1:op%Fshift(flavorb)+op%tails(flavorb) ) = &
1551 : CSHIFT( op%M%mat_tau(1:op%sumtails ,op%Fshift(flavora)+&
1552 0 : & op%tails(flavorb)+1:op%Fshift(flavorb)+op%tails(flavorb) ) , SHIFT=1 , DIM=2)
1553 : enddo
1554 : end if
1555 : end if
1556 : endif
1557 : if(3==4) then
1558 : op%M=op%M_update
1559 : endif
1560 :
1561 :
1562 0 : do iflavort=flavora+1,flavorb
1563 0 : op%Fshift(iflavort)=op%Fshift(iflavort)+op%tails(flavorb)-op%tails(flavora)
1564 : enddo
1565 :
1566 0 : itmptail=op%tails(flavora)
1567 0 : op%tails(flavora)=op%tails(flavorb)
1568 0 : op%tails(flavorb)=itmptail
1569 0 : FREE(mat_temp)
1570 0 : FREE(mat_tau_temp)
1571 :
1572 0 : END SUBROUTINE BathOperatoroffdiag_swap
1573 : !!***
1574 :
1575 : !!****f* ABINIT/m_BathOperatoroffdiag/BathOperatoroffdiag_initF
1576 : !! NAME
1577 : !! BathOperatoroffdiag_initF
1578 : !!
1579 : !! FUNCTION
1580 : !! Copy input hybridization functions from a file
1581 : !!
1582 : !! COPYRIGHT
1583 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
1584 : !! This file is distributed under the terms of the
1585 : !! GNU General Public License, see ~abinit/COPYING
1586 : !! or http://www.gnu.org/copyleft/gpl.txt .
1587 : !!
1588 : !! INPUTS
1589 : !! op=bath operator
1590 : !! ifstream=file stream to read F
1591 : !!
1592 : !! OUTPUT
1593 : !! argout(sizeout)=description
1594 : !!
1595 : !! SIDE EFFECTS
1596 : !!
1597 : !! NOTES
1598 : !!
1599 : !! SOURCE
1600 :
1601 0 : SUBROUTINE BathOperatoroffdiag_initF(op,ifstream)
1602 :
1603 : !Arguments ----------------------
1604 : TYPE(BathOperatoroffdiag), INTENT(INOUT) :: op
1605 : INTEGER , INTENT(IN ) :: ifstream
1606 : !Local variables ----------------
1607 : INTEGER :: iflavor1
1608 : INTEGER :: iflavor2
1609 : INTEGER :: sample
1610 :
1611 0 : IF ( op%set .EQV. .FALSE. ) &
1612 0 : CALL ERROR("BathOperatoroffdiag_initF : BathOperatoroffdiag not set ")
1613 :
1614 0 : DO iflavor1=1,op%flavors
1615 0 : DO iflavor2=2,op%flavors
1616 0 : DO sample = 1, op%sizeHybrid
1617 0 : READ(ifstream,*) op%F(sample,iflavor1,iflavor2)
1618 : END DO
1619 : END DO
1620 : END DO
1621 0 : END SUBROUTINE BathOperatoroffdiag_initF
1622 : !!***
1623 :
1624 : !!****f* ABINIT/m_BathOperatoroffdiag/BathOperatoroffdiag_setF
1625 : !! NAME
1626 : !! BathOperatoroffdiag_setF
1627 : !!
1628 : !! FUNCTION
1629 : !! Copy F from input array
1630 : !!
1631 : !! COPYRIGHT
1632 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
1633 : !! This file is distributed under the terms of the
1634 : !! GNU General Public License, see ~abinit/COPYING
1635 : !! or http://www.gnu.org/copyleft/gpl.txt .
1636 : !!
1637 : !! INPUTS
1638 : !! op=bath operator
1639 : !! F=array of the hybridization function
1640 : !!
1641 : !! OUTPUT
1642 : !!
1643 : !! SIDE EFFECTS
1644 : !!
1645 : !! NOTES
1646 : !!
1647 : !! SOURCE
1648 :
1649 0 : SUBROUTINE BathOperatoroffdiag_setF(op,F)
1650 :
1651 : !Arguments ------------------------------------
1652 : TYPE(BathOperatoroffdiag) , INTENT(INOUT) :: op
1653 : DOUBLE PRECISION, DIMENSION(:,:,:) , INTENT(IN ) :: F
1654 : !Arguments ------------------------------------
1655 : INTEGER :: iflavor1
1656 : INTEGER :: iflavor2
1657 : INTEGER :: sample
1658 : INTEGER :: length
1659 :
1660 0 : IF ( op%set .EQV. .FALSE. ) &
1661 0 : CALL ERROR("BathOperatoroffdiag_setF : BathOperatoroffdiag not set ")
1662 :
1663 0 : length = SIZE(F)
1664 0 : IF ( length .NE. (op%flavors * op%flavors * op%sizeHybrid) ) &
1665 0 : CALL ERROR("BathOperatoroffdiag_setF : wrong input F ")
1666 :
1667 0 : DO iflavor1=1,op%flavors
1668 0 : DO iflavor2=1,op%flavors
1669 0 : DO sample = 1, op%sizeHybrid
1670 0 : op%F(sample,iflavor1,iflavor2) = F(sample,iflavor1,iflavor2)
1671 : END DO
1672 : END DO
1673 : END DO
1674 0 : END SUBROUTINE BathOperatoroffdiag_setF
1675 : !!***
1676 :
1677 : !!****f* ABINIT/m_BathOperatoroffdiag/BathOperatoroffdiag_printF
1678 : !! NAME
1679 : !! BathOperatoroffdiag_printF
1680 : !!
1681 : !! FUNCTION
1682 : !! print F function
1683 : !!
1684 : !! COPYRIGHT
1685 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
1686 : !! This file is distributed under the terms of the
1687 : !! GNU General Public License, see ~abinit/COPYING
1688 : !! or http://www.gnu.org/copyleft/gpl.txt .
1689 : !!
1690 : !! INPUTS
1691 : !! op=bath operator
1692 : !! ostream=file stream to write in
1693 : !!
1694 : !! OUTPUT
1695 : !!
1696 : !! SIDE EFFECTS
1697 : !!
1698 : !! NOTES
1699 : !!
1700 : !! SOURCE
1701 :
1702 0 : SUBROUTINE BathOperatoroffdiag_printF(op,ostream)
1703 :
1704 : !Arguments ------------------------------------
1705 : TYPE(BathOperatoroffdiag), INTENT(INOUT) :: op
1706 : INTEGER,OPTIONAL , INTENT(IN ) :: ostream
1707 : !Local variables ------------------------------
1708 : CHARACTER(LEN=4) :: aflavor
1709 : CHARACTER(LEN=50) :: string
1710 : INTEGER :: iflavor1
1711 : INTEGER :: iflavor2
1712 : INTEGER :: sample
1713 : INTEGER :: ostream_val
1714 :
1715 0 : IF ( PRESENT(ostream) ) THEN
1716 0 : ostream_val = ostream
1717 : ELSE
1718 0 : ostream_val = 65
1719 0 : OPEN(UNIT=ostream_val, FILE="F.dat")
1720 : END IF
1721 :
1722 0 : WRITE(aflavor,'(I4)') (op%flavors*op%flavors+1)
1723 0 : string = '(1x,'//TRIM(ADJUSTL(aflavor))//'E22.14)'
1724 0 : DO sample = 1, op%sizeHybrid
1725 0 : WRITE(ostream_val,string) (sample-1)*op%dt, ((op%F(sample,iflavor1,iflavor2),&
1726 0 : iflavor1=1,op%flavors),iflavor2=1,op%flavors)
1727 : END DO
1728 : !CALL FLUSH(ostream_val)
1729 :
1730 0 : IF ( .NOT. PRESENT(ostream) ) &
1731 0 : CLOSE(ostream_val)
1732 :
1733 0 : END SUBROUTINE BathOperatoroffdiag_printF
1734 : !!***
1735 :
1736 : !!****f* ABINIT/m_BathOperatoroffdiag/BathOperatoroffdiag_printM
1737 : !! NAME
1738 : !! BathOperatoroffdiag_printM
1739 : !!
1740 : !! FUNCTION
1741 : !! print M =F^{-1} matrix
1742 : !!
1743 : !! COPYRIGHT
1744 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
1745 : !! This file is distributed under the terms of the
1746 : !! GNU General Public License, see ~abinit/COPYING
1747 : !! or http://www.gnu.org/copyleft/gpl.txt .
1748 : !!
1749 : !! INPUTS
1750 : !! op=bath operator
1751 : !! ostream=file stream to write in
1752 : !!
1753 : !! OUTPUT
1754 : !! argout(sizeout)=description
1755 : !!
1756 : !! SIDE EFFECTS
1757 : !!
1758 : !! NOTES
1759 : !!
1760 : !! SOURCE
1761 :
1762 0 : SUBROUTINE BathOperatoroffdiag_printM(op,ostream)
1763 :
1764 : !Arguments ------------------------------------
1765 : TYPE(BathOperatoroffdiag), INTENT(IN) :: op
1766 : INTEGER, OPTIONAL , INTENT(IN) :: ostream
1767 : !Local variables ------------------------------
1768 : INTEGER :: ostream_val
1769 :
1770 0 : IF ( op%activeFlavor .LE. 0 ) &
1771 0 : CALL ERROR("BathOperatoroffdiag_printM : no active hybrid function ")
1772 0 : ostream_val = 6
1773 0 : IF ( PRESENT(ostream) ) ostream_val = ostream
1774 0 : CALL MatrixHyb_print(op%M,ostream_val)
1775 0 : END SUBROUTINE BathOperatoroffdiag_printM
1776 : !!***
1777 :
1778 : !!****f* ABINIT/m_BathOperatoroffdiag/BathOperatoroffdiag_printM_matrix
1779 : !! NAME
1780 : !! BathOperatoroffdiag_printM_matrix
1781 : !!
1782 : !! FUNCTION
1783 : !! print M =F^{-1} matrix
1784 : !!
1785 : !! COPYRIGHT
1786 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
1787 : !! This file is distributed under the terms of the
1788 : !! GNU General Public License, see ~abinit/COPYING
1789 : !! or http://www.gnu.org/copyleft/gpl.txt .
1790 : !!
1791 : !! INPUTS
1792 : !! op=bath operator
1793 : !! ostream=file stream to write in
1794 : !!
1795 : !! OUTPUT
1796 : !! argout(sizeout)=description
1797 : !!
1798 : !! SIDE EFFECTS
1799 : !!
1800 : !! NOTES
1801 : !!
1802 : !! SOURCE
1803 :
1804 0 : SUBROUTINE BathOperatoroffdiag_printM_matrix(op,ostream)
1805 :
1806 : !Arguments ------------------------------------
1807 : TYPE(BathOperatoroffdiag), INTENT(IN) :: op
1808 : INTEGER, OPTIONAL , INTENT(IN) :: ostream
1809 : !Local variables ------------------------------
1810 : INTEGER :: iflavor1
1811 : INTEGER :: i1,it1,it2
1812 : CHARACTER(LEN=22) :: string
1813 : CHARACTER(LEN=22) :: string2
1814 : CHARACTER(LEN=4 ) :: size
1815 :
1816 : ABI_UNUSED(ostream)
1817 :
1818 0 : WRITE(size,'(I4)') op%sumtails
1819 0 : string ='(i2,x,i3,a,'//TRIM(ADJUSTL(size))//'(E5.2,1x))'
1820 0 : string2 ='(6x,'//TRIM(ADJUSTL(size))//'(i6))'
1821 0 : open(unit=222, file="M_matrix.dat")
1822 0 : open(unit=223, file="M_matrix_tau.dat")
1823 0 : it1=0
1824 0 : write(222,string2) ((i1,i1=1,op%tails(iflavor1)),iflavor1=1,op%flavors)
1825 0 : do iflavor1=1, op%flavors
1826 0 : do i1=1, op%tails(iflavor1)
1827 0 : it1=it1+1
1828 0 : write(222,string) iflavor1,i1,'|',(op%M%mat(it1,it2),it2=1,op%sumtails)
1829 : enddo
1830 : enddo
1831 :
1832 :
1833 0 : END SUBROUTINE BathOperatoroffdiag_printM_matrix
1834 : !!***
1835 :
1836 : !!****f* ABINIT/m_BathOperatoroffdiag/ BathOperatoroffdiag_destroy
1837 : !! NAME
1838 : !! BathOperatoroffdiag_destroy
1839 : !!
1840 : !! FUNCTION
1841 : !! Deallocate and reset every thing
1842 : !!
1843 : !! COPYRIGHT
1844 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
1845 : !! This file is distributed under the terms of the
1846 : !! GNU General Public License, see ~abinit/COPYING
1847 : !! or http://www.gnu.org/copyleft/gpl.txt .
1848 : !!
1849 : !! INPUTS
1850 : !! op=bath operator
1851 : !!
1852 : !! OUTPUT
1853 : !!
1854 : !! SIDE EFFECTS
1855 : !!
1856 : !! NOTES
1857 : !!
1858 : !! SOURCE
1859 :
1860 0 : SUBROUTINE BathOperatoroffdiag_destroy(op)
1861 :
1862 : TYPE(BathOperatoroffdiag), INTENT(INOUT) :: op
1863 :
1864 0 : CALL MatrixHyb_destroy(op%M)
1865 0 : CALL MatrixHyb_destroy(op%M_update)
1866 :
1867 0 : CALL Vector_destroy(op%R)
1868 0 : CALL Vector_destroy(op%Q)
1869 0 : CALL Vector_destroy(op%Rtau)
1870 0 : CALL Vector_destroy(op%Qtau)
1871 0 : FREEIF(op%F)
1872 0 : FREEIF(op%Fshift)
1873 0 : FREEIF(op%tails)
1874 :
1875 0 : op%MAddFlag = .FALSE.
1876 0 : op%MRemoveFlag = .FALSE.
1877 0 : op%flavors = 0
1878 0 : op%beta = 0.d0
1879 0 : op%dt = 0.d0
1880 0 : op%inv_dt = 0.d0
1881 0 : op%samples = 0
1882 0 : op%sizeHybrid = 0
1883 0 : op%activeFlavor = 0
1884 0 : op%updatePosRow = 0
1885 0 : op%updatePosCol = 0
1886 :
1887 0 : END SUBROUTINE BathOperatoroffdiag_destroy
1888 : !!***
1889 :
1890 : !!****f* ABINIT/m_BathOperatoroffdiag/BathOperatoroffdiag_doCheck
1891 : !! NAME
1892 : !! BathOperatoroffdiag_doCheck
1893 : !!
1894 : !! FUNCTION
1895 : !! Just store if we perfom check for updates of M
1896 : !!
1897 : !! COPYRIGHT
1898 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
1899 : !! This file is distributed under the terms of the
1900 : !! GNU General Public License, see ~abinit/COPYING
1901 : !! or http://www.gnu.org/copyleft/gpl.txt .
1902 : !!
1903 : !! INPUTS
1904 : !! op=bath operator
1905 : !! opt_check=second bit should be one
1906 : !!
1907 : !! OUTPUT
1908 : !!
1909 : !! SIDE EFFECTS
1910 : !!
1911 : !! NOTES
1912 : !!
1913 : !! SOURCE
1914 :
1915 0 : SUBROUTINE BathOperatoroffdiag_doCheck(op,opt_check)
1916 :
1917 : !Arguments ------------------------------------
1918 : TYPE(BathOperatoroffdiag) , INTENT(INOUT) :: op
1919 : INTEGER , INTENT(IN ) :: opt_check
1920 :
1921 0 : IF ( opt_check .GE. 2 ) &
1922 0 : op%doCheck = .TRUE.
1923 0 : END SUBROUTINE BathOperatoroffdiag_doCheck
1924 : !!***
1925 :
1926 : !!****f* ABINIT/m_BathOperatoroffdiag/BathOperatoroffdiag_checkM
1927 : !! NAME
1928 : !! BathOperatoroffdiag_checkM
1929 : !!
1930 : !! FUNCTION
1931 : !! compute from scratch the M matrix and compar it
1932 : !! with the already computed M matrix
1933 : !!
1934 : !! COPYRIGHT
1935 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
1936 : !! This file is distributed under the terms of the
1937 : !! GNU General Public License, see ~abinit/COPYING
1938 : !! or http://www.gnu.org/copyleft/gpl.txt .
1939 : !!
1940 : !! INPUTS
1941 : !! op=bath operator
1942 : !! particle=list of all segments of the active flavor
1943 : !!
1944 : !! OUTPUT
1945 : !!
1946 : !! SIDE EFFECTS
1947 : !!
1948 : !! NOTES
1949 : !!
1950 : !! SOURCE
1951 :
1952 0 : SUBROUTINE BathOperatoroffdiag_checkM(op,particle)
1953 :
1954 : !Arguments ------------------------------------
1955 : TYPE(BathOperatoroffdiag) , INTENT(INOUT) :: op
1956 : TYPE(ListCdagC) , INTENT(IN ) :: particle(:)
1957 : !Local variables ------------------------------
1958 : ! TYPE(MatrixHyb) :: checkMatrix
1959 : LOGICAL :: checkTau
1960 : INTEGER :: tail
1961 : INTEGER :: iC
1962 : INTEGER :: iCdag
1963 : INTEGER :: aF
1964 : INTEGER :: iflavora
1965 : INTEGER :: iflavorb,it !,it1
1966 : CHARACTER(LEN=6) :: a
1967 : DOUBLE PRECISION :: time
1968 : DOUBLE PRECISION :: beta
1969 : DOUBLE PRECISION :: mbeta_two
1970 : DOUBLE PRECISION :: errorabs
1971 : DOUBLE PRECISION :: errormax
1972 : DOUBLE PRECISION :: error1
1973 : DOUBLE PRECISION :: errorrel
1974 : DOUBLE PRECISION :: tc
1975 : DOUBLE PRECISION :: tCdag
1976 : DOUBLE PRECISION :: sumMmat
1977 : DOUBLE PRECISION :: sumCheck
1978 : #include "BathOperatoroffdiag_hybrid.h"
1979 :
1980 0 : aF = op%activeFlavor
1981 : !Construction de la matrix
1982 0 : tail = op%sumtails
1983 : ! CALL MatrixHyb_init(checkMatrix,op%iTech,size=tail,Wmax=op%samples)
1984 : ! CALL MatrixHyb_setSize(checkMatrix,tail)
1985 :
1986 : ! --- set size of the matrix
1987 0 : CALL MatrixHyb_setSize(op%M_update,tail)
1988 :
1989 : ! --- compute useful quantities
1990 0 : beta = op%beta
1991 0 : mbeta_two = -beta*0.5d0
1992 0 : op%checkNumber = op%checkNumber + 1
1993 0 : IF ( tail .NE. op%M%tail ) THEN
1994 0 : CALL WARN("BathOperatoroffdiag_checkM : tails are different ")
1995 0 : RETURN
1996 : END IF
1997 :
1998 0 : do it=1,op%sumtails
1999 : !write(6,*) " checkM begin M_update%mat_tau",(op%M_update%mat_tau(it,it1),it1=1,op%sumtails)
2000 : enddo
2001 : ! --- build matrix
2002 : !CALL ListCdagC_print(particle)
2003 0 : DO iflavora = 1, op%flavors
2004 0 : DO iCdag = 1, op%tails(iflavora)
2005 0 : tCdag = particle(iflavora)%list(iCdag,Cdag_)
2006 : !write(6,*) " checkM a",iflavora,tCdag
2007 0 : DO iflavorb = 1, op%flavors
2008 0 : DO iC = 1, op%tails(iflavorb)
2009 : !tC = particle%list(C_,iC).MOD.beta
2010 0 : MODCYCLE(particle(iflavorb)%list(iC,C_),beta,tC) ! tC is tC, or Tc-Beta if tc>beta
2011 : !write(6,*) " checkM b",iflavorb,tC
2012 0 : time = tC - tCdag ! time is positive or negative but lower than beta
2013 : !write(6,*) " checkM time",time
2014 :
2015 : #include "BathOperatoroffdiag_hybrid"
2016 :
2017 0 : op%M_update%mat(op%Fshift(iflavorb)+iC,op%Fshift(iflavora)+iCdag) = hybrid
2018 :
2019 0 : time = time + ( SIGN(1.d0,time) - 1.d0 )*mbeta_two
2020 0 : op%M_update%mat_tau(op%Fshift(iflavora)+iCdag,op%Fshift(iflavorb)+iC) = INT ( (time*op%inv_dt) +1.5d0 )
2021 : !write(6,*) " checkM mat_tau",INT ( (time*op%inv_dt) +1.5d0 )
2022 : !write(6,*) " checkM shifts",op%Fshift(iflavorb),iCdag,op%Fshift(iflavora),iC
2023 : END DO ! iC
2024 : END DO ! iflavorb
2025 : END DO ! iCdag
2026 : END DO ! iflavora
2027 :
2028 : ! CALL MatrixHyb_Print(checkMatrix)
2029 : ! --- Inverse matrix
2030 0 : CALL MatrixHyb_inverse(op%M_update)
2031 :
2032 : ! CALL MatrixHyb_Print(checkMatrix)
2033 0 : do it=1,op%sumtails
2034 : !write(6,*) " checkM end M_update%mat_tau",(op%M_update%mat_tau(it,it1),it1=1,op%sumtails)
2035 : enddo
2036 : do it=1,op%sumtails
2037 : !write(6,*) " checkM end M_update",(op%M%mat(it,it1),it1=1,op%sumtails)
2038 : enddo
2039 :
2040 : ! --- Compare M_update and M to check if calculation of M is correct
2041 : sumMmat =0.d0
2042 : sumCheck=0.d0
2043 : error1 = 0.d0
2044 : errormax = 0.d0
2045 : checkTau = .FALSE.
2046 0 : DO iCdag = 1, tail
2047 0 : Do iC =1, tail
2048 : errorrel= ABS((op%M_update%mat(iC, iCdag) - &
2049 0 : op%M%mat(iC,iCdag))/op%M_update%mat(iC,iCdag))
2050 : errorabs= ABS(op%M_update%mat(iC, iCdag) - &
2051 0 : op%M%mat(iC,iCdag))
2052 0 : IF ( errorrel .gt. errormax .and. errorabs .gt. 0.001d0 ) errormax = errorrel
2053 : ! write(6,*) " checkM ", errorrel,errorabs
2054 0 : IF ( op%M_update%mat_tau(iC,iCdag) .NE. op%M%mat_tau(iC,iCdag) ) then
2055 0 : checkTau = .TRUE.
2056 : !write(6,*) "op%M_update%mat_tau(iC,iCdag), op%M%mat_tau(iC,iCdag)",op%M_update%mat_tau(iC,iCdag), op%M%mat_tau(iC,iCdag)
2057 : !call flush(6)
2058 0 : CALL ERROR("BathOperatoroffdiag_checkM : "//a//"% ")
2059 : ENDIF
2060 :
2061 : END DO
2062 : END DO
2063 :
2064 0 : IF ( checkTau .EQV. .TRUE. ) THEN
2065 0 : CALL WARN("BathOperatoroffdiag_checkM : mat_tau differs should be")
2066 0 : CALL MatrixHyb_print(op%M_update,opt_print=1)
2067 0 : CALL WARN("BathOperatoroffdiag_checkM : whereas it is")
2068 0 : CALL MatrixHyb_print(op%M,opt_print=1)
2069 : END IF
2070 0 : op%meanError = op%meanError + errormax
2071 0 : IF ( errormax .GT. 1.d0 ) THEN
2072 0 : WRITE(a,'(I4)') INT(error1*100.d0)
2073 : !write(6,'(I4)') INT(error1*100.d0)
2074 : ! CALL MatrixHyb_Print(op%M)
2075 0 : CALL WARN("BathOperatoroffdiag_checkM")
2076 : END IF
2077 : ! CALL MatrixHyb_destroy(checkMatrix)
2078 : END SUBROUTINE BathOperatoroffdiag_checkM
2079 : !!***
2080 :
2081 : !!****f* ABINIT/m_BathOperatoroffdiag/BathOperatoroffdiag_recomputeM
2082 : !! NAME
2083 : !! BathOperatoroffdiag_recomputeM
2084 : !!
2085 : !! FUNCTION
2086 : !! compute from scratch the M matrix
2087 : !!
2088 : !! COPYRIGHT
2089 : !! Copyright (C) 2013-2026 ABINIT group (B. Amadon)
2090 : !! This file is distributed under the terms of the
2091 : !! GNU General Public License, see ~abinit/COPYING
2092 : !! or http://www.gnu.org/copyleft/gpl.txt .
2093 : !!
2094 : !! INPUTS
2095 : !! op=bath operator
2096 : !! particle=list of all segments of the active flavor
2097 : !!
2098 : !! OUTPUT
2099 : !!
2100 : !! SIDE EFFECTS
2101 : !!
2102 : !! NOTES
2103 : !!
2104 : !! SOURCE
2105 :
2106 0 : SUBROUTINE BathOperatoroffdiag_recomputeM(op,particle,flav_i,flav_j)
2107 :
2108 : !Arguments ------------------------------------
2109 : TYPE(BathOperatoroffdiag) , INTENT(INOUT) :: op
2110 : TYPE(ListCdagC) , INTENT(IN ) :: particle(:)
2111 : INTEGER :: flav_i,flav_j
2112 : !Local variables ------------------------------
2113 : ! TYPE(MatrixHyb) :: checkMatrix
2114 : INTEGER :: tail
2115 : INTEGER :: iC
2116 : INTEGER :: iCdag
2117 : INTEGER :: aF
2118 : INTEGER :: iflavora
2119 : INTEGER :: iflavorb,it !,it1
2120 : INTEGER :: iflavora_imp
2121 : INTEGER :: iflavorb_imp
2122 : !CHARACTER(LEN=6) :: a
2123 : DOUBLE PRECISION :: time
2124 : DOUBLE PRECISION :: beta
2125 : DOUBLE PRECISION :: mbeta_two
2126 : DOUBLE PRECISION :: tc
2127 : DOUBLE PRECISION :: tCdag
2128 : !DOUBLE PRECISION :: sumMmat
2129 : !DOUBLE PRECISION :: sumCheck
2130 : #include "BathOperatoroffdiag_hybrid.h"
2131 :
2132 0 : aF = op%activeFlavor
2133 : !Construction de la matrix
2134 0 : tail = op%sumtails
2135 : ! CALL MatrixHyb_init(checkMatrix,op%iTech,size=tail,Wmax=op%samples)
2136 : ! CALL MatrixHyb_setSize(checkMatrix,tail)
2137 :
2138 : ! --- set size of the matrix
2139 0 : CALL MatrixHyb_setSize(op%M_update,tail)
2140 :
2141 : ! --- compute useful quantities
2142 0 : beta = op%beta
2143 0 : mbeta_two = -beta*0.5d0
2144 0 : op%checkNumber = op%checkNumber + 1
2145 0 : IF ( tail .NE. op%M%tail ) THEN
2146 0 : CALL WARN("BathOperatoroffdiag_checkM : tails are different ")
2147 0 : RETURN
2148 : END IF
2149 :
2150 0 : do it=1,op%sumtails
2151 : !write(6,*) " checkM begin M_update%mat_tau",(op%M_update%mat_tau(it,it1),it1=1,op%sumtails)
2152 : enddo
2153 : ! --- build matrix
2154 : !CALL ListCdagC_print(particle)
2155 0 : DO iflavora = 1, op%flavors
2156 0 : iflavora_imp=iflavora
2157 0 : if(iflavora==flav_i) iflavora_imp=flav_j
2158 0 : if(iflavora==flav_j) iflavora_imp=flav_i
2159 0 : DO iCdag = 1, op%tails(iflavora_imp)
2160 0 : tCdag = particle(iflavora_imp)%list(iCdag,Cdag_)
2161 : !write(6,*) " checkM a",iflavora,tCdag
2162 0 : DO iflavorb = 1, op%flavors
2163 0 : iflavorb_imp=iflavorb
2164 0 : if(iflavorb==flav_j) iflavorb_imp=flav_i
2165 0 : if(iflavorb==flav_i) iflavorb_imp=flav_j
2166 0 : DO iC = 1, op%tails(iflavorb_imp)
2167 : !tC = particle%list(C_,iC).MOD.beta
2168 0 : MODCYCLE(particle(iflavorb_imp)%list(iC,C_),beta,tC) ! tC is tC, or Tc-Beta if tc>beta
2169 : !write(6,*) " checkM b",iflavorb,tC
2170 0 : time = tC - tCdag ! time is positive or negative but lower than beta
2171 : !write(6,*) " checkM time",time
2172 :
2173 : #include "BathOperatoroffdiag_hybrid"
2174 :
2175 0 : op%M_update%mat(op%Fshift(iflavorb_imp)+iC,op%Fshift(iflavora_imp)+iCdag) = hybrid
2176 :
2177 0 : time = time + ( SIGN(1.d0,time) - 1.d0 )*mbeta_two
2178 0 : op%M_update%mat_tau(op%Fshift(iflavora_imp)+iCdag,op%Fshift(iflavorb_imp)+iC) = INT ( (time*op%inv_dt) +1.5d0 )
2179 : !write(6,*) " checkM mat_tau",INT ( (time*op%inv_dt) +1.5d0 )
2180 : !write(6,*) " checkM shifts",op%Fshift(iflavorb),iCdag,op%Fshift(iflavora),iC
2181 : END DO ! iC
2182 : END DO ! iflavorb
2183 : END DO ! iCdag
2184 : END DO ! iflavora
2185 :
2186 : ! CALL MatrixHyb_Print(checkMatrix)
2187 : ! --- Inverse matrix
2188 0 : CALL MatrixHyb_inverse(op%M_update)
2189 :
2190 : ! CALL MatrixHyb_Print(checkMatrix)
2191 0 : do it=1,op%sumtails
2192 : !write(6,*) " checkM end M_update%mat_tau",(op%M_update%mat_tau(it,it1),it1=1,op%sumtails)
2193 : enddo
2194 : do it=1,op%sumtails
2195 : !write(6,*) " checkM end M_update",(op%M%mat(it,it1),it1=1,op%sumtails)
2196 : enddo
2197 :
2198 : ! --- Compare M_update and M to check if calculation of M is correct
2199 : END SUBROUTINE BathOperatoroffdiag_recomputeM
2200 : !!***
2201 :
2202 : !!****f* ABINIT/m_BathOperatoroffdiag/BathOperatoroffdiag_getError
2203 : !! NAME
2204 : !! BathOperatoroffdiag_getError
2205 : !!
2206 : !! FUNCTION
2207 : !! compute a percentage error / checkM
2208 : !!
2209 : !! COPYRIGHT
2210 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
2211 : !! This file is distributed under the terms of the
2212 : !! GNU General Public License, see ~abinit/COPYING
2213 : !! or http://www.gnu.org/copyleft/gpl.txt .
2214 : !!
2215 : !! INPUTS
2216 : !! op=bath operator
2217 : !!
2218 : !! OUTPUT
2219 : !! BathOperatoroffdiag_getError=Error in percent
2220 : !!
2221 : !! SIDE EFFECTS
2222 : !!
2223 : !! NOTES
2224 : !!
2225 : !! SOURCE
2226 :
2227 0 : DOUBLE PRECISION FUNCTION BathOperatoroffdiag_getError(op)
2228 :
2229 : TYPE(BathOperatoroffdiag), INTENT(IN) :: op
2230 :
2231 0 : IF ( op%doCheck .EQV. .TRUE. ) THEN
2232 0 : BathOperatoroffdiag_getError = op%meanError / DBLE(op%checkNumber)
2233 : ELSE
2234 : BathOperatoroffdiag_getError = 0.d0
2235 : END IF
2236 0 : END FUNCTION BathOperatoroffdiag_getError
2237 : !!***
2238 : !#endif
2239 :
2240 0 : END MODULE m_BathOperatoroffdiag
2241 : !!***
|
