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