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