Line data Source code
1 : !!****m* ABINIT/m_green
2 : !! NAME
3 : !! m_green
4 : !!
5 : !! FUNCTION
6 : !!
7 : !! COPYRIGHT
8 : !! Copyright (C) 2006-2026 ABINIT group (BAmadon)
9 : !! This file is distributed under the terms of the
10 : !! GNU General Public License, see ~abinit/COPYING
11 : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : !!
13 : !! SOURCE
14 :
15 : #if defined HAVE_CONFIG_H
16 : #include "config.h"
17 : #endif
18 :
19 :
20 : #include "abi_common.h"
21 :
22 : ! nvtx related macro definition
23 : #include "nvtx_macros.h"
24 :
25 : MODULE m_green
26 :
27 : use, intrinsic :: iso_c_binding, only: c_size_t
28 : use defs_basis
29 : use m_abicore
30 : use m_errors
31 :
32 : use m_crystal, only : crystal_t
33 : use m_fstrings, only : int2char4
34 : use m_hide_lapack, only : xginv
35 : use m_io_tools, only : open_file
36 : use m_lib_four
37 : use m_matlu, only : add_matlu,copy_matlu,destroy_matlu,diff_matlu,init_matlu,matlu_type,print_matlu, &
38 : & prod_matlu,shift_matlu,sym_matlu,trace_matlu,trace_prod_matlu,xmpi_matlu,zero_matlu
39 : use m_oper, only : copy_oper,copy_oper_from_ndat,copy_oper_to_ndat,destroy_oper,downfold_oper,gather_oper, &
40 : & gather_oper_ks,init_oper,init_oper_ndat,inverse_oper,oper_type,print_oper,prod_oper, &
41 : & trace_oper,trace_prod_oper,upfold_oper
42 : use m_paw_dmft, only : construct_nwli_dmft,mpi_distrib_dmft_type,paw_dmft_type
43 : use m_self, only : self_type
44 : use m_splines
45 : use m_time, only : timab
46 : use m_xmpi, only : xmpi_barrier,xmpi_sum
47 : use m_gputk
48 : use m_abi_linalg
49 :
50 :
51 : #ifdef HAVE_GPU_MARKERS
52 : use m_nvtx_data
53 : #endif
54 :
55 : implicit none
56 :
57 : private
58 :
59 : public :: init_green
60 : public :: destroy_green
61 : public :: init_green_tau
62 : public :: destroy_green_tau
63 : public :: print_green
64 : public :: printocc_green
65 : public :: compute_green
66 : public :: integrate_green
67 : public :: icip_green
68 : public :: fourier_green
69 : public :: check_fourier_green
70 : public :: compa_occup_ks
71 : public :: copy_green
72 : public :: occup_green_tau
73 : public :: add_int_fct
74 : public :: int_fct
75 : public :: fourier_fct
76 : public :: spline_fct
77 : public :: distrib_paral
78 : public :: greendftcompute_green
79 : public :: fermi_green
80 : public :: newton
81 : public :: local_ks_green
82 : public :: compute_moments_ks
83 : public :: compute_trace_moments_ks
84 : public :: compute_moments_loc
85 : public :: occup_fd
86 : !!***
87 :
88 : !!****t* m_green/green_type
89 : !! NAME
90 : !! green_type
91 : !!
92 : !! FUNCTION
93 : !! This structured datatype contains the necessary data
94 : !!
95 : !! SOURCE
96 :
97 : type, public :: green_type ! for each atom
98 :
99 : integer :: dmft_nwli
100 : ! Linear index of the last imaginary frequency
101 :
102 : integer :: dmft_nwlo
103 : ! Number of imaginary frequencies
104 :
105 : integer :: dmftqmc_l
106 : ! Number of time slices for QMC
107 :
108 : integer :: fileprt_tau
109 : ! 1 if file created
110 :
111 : integer :: fileprt_w
112 : ! 1 if file created
113 :
114 : integer :: has_charge_matlu
115 : ! =2 if calculation of LOCAL CORRELATED occupations is done
116 :
117 : integer :: has_charge_matlu_prev
118 : ! =0 charge_matlu_prev not allocated
119 : ! =1 charge_matlu_prev is allocated
120 : ! =2 charge_matlu_prev is calculated (ie calculation of LOCAL CORRELATED occupations is done from
121 : ! solver green function)
122 :
123 : integer :: has_charge_matlu_solver
124 : ! =0 charge_matlu_solver not allocated
125 : ! =1 charge_matlu_solver is allocated
126 : ! =2 charge_matlu_solver is calculated (ie calculation of LOCAL CORRELATED occupations is done from
127 : ! solver green function)
128 :
129 : integer :: has_greenmatlu_xsum
130 : ! =1 green%oper%matlu xsum in compute_green
131 : ! =0 green%oper%matlu non xsumed in compute_green
132 : ! used in integrate_green to checked that green function was computed in
133 : ! integrate_green.
134 :
135 : integer :: has_moments
136 : ! =1 if the high-frequency moments are computed
137 :
138 : integer :: ichargeloc_cv
139 :
140 : integer :: ifermie_cv
141 :
142 : integer :: nmoments
143 : ! Number of high-frequency moments which will be computed
144 :
145 : integer :: nw
146 : ! Number of frequencies
147 :
148 : integer :: use_oper_tau_ks
149 : ! 0 do not use oper_tau_ks
150 : ! 1 use oper_tau_ks
151 :
152 : character(len=4) :: w_type
153 : ! type of frequencies used
154 :
155 : !character(len=12) :: whichgreen
156 : ! describe the type of green function computed (DFT, DMFT, ..)
157 :
158 : real(dp) :: charge_ks
159 : ! Total charge computed from ks orbitals
160 :
161 : real(dp) :: ekin_imp
162 : ! Kinetic energy of the impurity
163 :
164 : real(dp) :: fband_weiss
165 : ! Tr(log(G0))
166 :
167 : real(dp) :: integral
168 : ! Integral of the interaction energy divided by U
169 :
170 : real(dp) :: trace_log
171 : ! Tr(log(G)) in KS space
172 :
173 : !integer, allocatable :: procb(:,:)
174 :
175 : !integer, allocatable :: proct(:,:)
176 :
177 : type(oper_type) :: occup
178 : ! Occupation in different basis
179 :
180 : type(oper_type) :: occup_tau
181 : ! Occupation in different basis
182 :
183 : complex(dp) :: trace_fermie(12)
184 : ! Container to store useful quantities for a quick computation
185 : ! of the moments during the Fermi level search
186 :
187 : real(dp), allocatable :: charge_matlu(:,:)
188 : ! Total charge on correlated orbitals
189 : ! todo_ba name of charge_matlu is misleading: should be changed
190 :
191 : real(dp), allocatable :: charge_matlu_prev(:,:)
192 : ! Total charge on correlated orbitals from previous iteration
193 :
194 : real(dp), allocatable :: charge_matlu_solver(:,:)
195 : ! Total charge on correlated orbitals obtained from solver by
196 : ! integration over frequencies.
197 :
198 : real(dp), allocatable :: ecorr_qmc(:)
199 : ! Correlation energy for a given atom in qmc
200 :
201 : real(dp), allocatable :: tau(:)
202 : ! Value of time in imaginary space
203 :
204 : complex(dp), allocatable :: trace_moments_log_ks(:)
205 : ! Trace of the moments of log(G)+log(iw*Id) in KS space
206 :
207 : complex(dp), allocatable :: trace_moments_log_loc(:)
208 : ! Trace of the moments of log(G)+log(iw*Id) in local space
209 :
210 : type(oper_type), allocatable :: moments(:)
211 : ! High-frequency moments
212 :
213 : type(oper_type), allocatable :: oper(:)
214 : ! Green's function in different basis
215 :
216 : type(oper_type), allocatable :: oper_tau(:)
217 : ! Green's function in different basis
218 :
219 : real(dp), ABI_CONTIGUOUS pointer :: omega(:) => null()
220 : ! Value of frequencies
221 :
222 : type(mpi_distrib_dmft_type), pointer :: distrib => null()
223 : ! Datastructure for MPI parallelization
224 :
225 : end type green_type
226 :
227 : !----------------------------------------------------------------------
228 :
229 :
230 : CONTAINS
231 : !!***
232 :
233 : !!****f* m_green/init_green
234 : !! NAME
235 : !! init_green
236 : !!
237 : !! FUNCTION
238 : !! Allocate variables used in type green_type.
239 : !!
240 : !! INPUTS
241 : !! green <type(green_type)>= green function data
242 : !! energies_dmft = datastructure for dmft energy
243 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
244 : !! opt_oper_ksloc (optional) = option for init_oper
245 : !! wtype = "real" Green function will be computed for real frequencies
246 : !! = "imag" Green function will be computed for imaginary frequencies
247 : !! opt_moments = 1 : to allocate the high frequency moments
248 : !! opt_moments_ksloc = option to init green%moments
249 : !! opt_occup_ksloc = option to init green%occup
250 : !!
251 : !! OUTPUTS
252 : !! green = variable of type green_type
253 : !!
254 : !! SOURCE
255 :
256 551 : subroutine init_green(green,paw_dmft,opt_oper_ksloc,wtype,opt_moments,opt_moments_ksloc,opt_occup_ksloc)
257 :
258 : !Arguments ------------------------------------
259 : type(green_type), intent(inout) :: green
260 : type(paw_dmft_type), target, intent(in) :: paw_dmft
261 : integer, optional, intent(in) :: opt_moments,opt_moments_ksloc,opt_occup_ksloc,opt_oper_ksloc
262 : character(len=4), optional, intent(in) :: wtype
263 : !Local variables ------------------------------------
264 : integer :: i,ifreq,mkmem,natom,nsppol,nw,optmoments
265 : integer :: optmoments_ksloc,optoccup_ksloc,optoper_ksloc,shift
266 : !************************************************************************
267 :
268 551 : optoper_ksloc = 2
269 551 : if (present(opt_oper_ksloc)) optoper_ksloc = opt_oper_ksloc
270 :
271 551 : green%w_type = "imag"
272 551 : if (present(wtype)) green%w_type = wtype
273 :
274 551 : optmoments = 0
275 551 : if (present(opt_moments)) optmoments = opt_moments
276 :
277 551 : optoccup_ksloc = 3
278 551 : if (present(opt_occup_ksloc)) optoccup_ksloc = opt_occup_ksloc
279 :
280 551 : optmoments_ksloc = 3
281 551 : if (present(opt_moments_ksloc)) optmoments_ksloc = opt_moments_ksloc
282 :
283 551 : if (green%w_type == "imag") then
284 547 : nw = paw_dmft%dmft_nwlo
285 547 : green%omega => paw_dmft%omega_lo(:)
286 547 : green%distrib => paw_dmft%distrib
287 4 : else if (green%w_type == "real") then
288 4 : nw = size(paw_dmft%omega_r(:))
289 4 : green%omega => paw_dmft%omega_r(:)
290 4 : green%distrib => paw_dmft%distrib_r
291 : end if ! w_type
292 :
293 551 : natom = paw_dmft%natom
294 551 : nsppol = paw_dmft%nsppol
295 :
296 551 : green%dmft_nwlo = paw_dmft%dmft_nwlo
297 551 : green%dmft_nwli = paw_dmft%dmft_nwli
298 551 : green%charge_ks = zero
299 2204 : ABI_MALLOC(green%charge_matlu,(nsppol+1,natom))
300 7109 : green%charge_matlu(:,:) = zero
301 551 : green%has_charge_matlu = 1
302 551 : green%has_greenmatlu_xsum = 0
303 :
304 1653 : ABI_MALLOC(green%charge_matlu_solver,(nsppol+1,natom))
305 7109 : green%charge_matlu_solver(:,:) = zero
306 551 : green%has_charge_matlu_solver = 1
307 :
308 1653 : ABI_MALLOC(green%charge_matlu_prev,(nsppol+1,natom))
309 7109 : green%charge_matlu_prev(:,:) = zero
310 551 : green%has_charge_matlu_prev = 1
311 :
312 551 : call init_oper(paw_dmft,green%occup,opt_ksloc=optoccup_ksloc)
313 :
314 : ! build simple arrays to distribute the tasks in compute_green.
315 : !ABI_MALLOC(green%procb,(nw,paw_dmft%nkpt))
316 : !ABI_MALLOC(green%proct,(nw,0:paw_dmft%nproc-1))
317 :
318 : !call distrib_paral(paw_dmft%nkpt,paw_dmft%nproc,nw,nw_perproc,green%procb,green%proct)
319 551 : green%nw = nw
320 :
321 : ! need to distribute memory over frequencies
322 :
323 : !! begin of temporary modificatios
324 : ! ABI_MALLOC(green%oper,(green%nw_perproc))
325 : ! do ifreq=1,green%nw_perproc
326 : ! call init_oper(paw_dmft,green%oper(ifreq),opt_ksloc=optoper_ksloc)
327 : ! enddo
328 : !
329 : ! do ifreq=1,green%nw
330 : ! if(green%proct(ifreq,myproc)==1) then
331 : ! do ikpt = 1 , paw_dmft%nkpt
332 : ! if (green%procb(ifreq,ikpt)==myproc) then
333 : ! endif
334 : ! enddo ! ikpt
335 : ! endif ! parallelisation
336 : ! enddo ! ifreq
337 : !
338 : !
339 : !! end of temporary modificatios
340 132226 : ABI_MALLOC(green%oper,(nw))
341 131124 : do ifreq=1,nw
342 131124 : call init_oper(paw_dmft,green%oper(ifreq),opt_ksloc=optoper_ksloc)
343 : end do ! ifreq
344 551 : green%ichargeloc_cv = 0
345 551 : green%ifermie_cv = 0
346 :
347 551 : if (paw_dmft%dmft_solv >= 5) then
348 522 : ABI_MALLOC(green%ecorr_qmc,(natom))
349 488 : green%ecorr_qmc(:) = zero
350 : end if
351 :
352 551 : green%fileprt_tau = 0
353 551 : green%fileprt_w = 0
354 :
355 551 : green%has_moments = optmoments
356 :
357 551 : green%nmoments = 0
358 :
359 551 : if (green%has_moments == 1) then
360 0 : green%nmoments = 5
361 0 : shift = green%distrib%shiftk
362 0 : mkmem = green%distrib%nkpt_mem(green%distrib%me_kpt+1)
363 0 : ABI_MALLOC(green%moments,(green%nmoments))
364 0 : call init_oper(paw_dmft,green%moments(1),nkpt=mkmem,shiftk=shift,opt_ksloc=2)
365 0 : do i=2,green%nmoments
366 0 : call init_oper(paw_dmft,green%moments(i),nkpt=mkmem,shiftk=shift,opt_ksloc=optmoments_ksloc)
367 : end do ! i
368 0 : if (paw_dmft%dmft_triqs_entropy == 1) then
369 0 : ABI_MALLOC(green%trace_moments_log_ks,(green%nmoments-1))
370 0 : ABI_MALLOC(green%trace_moments_log_loc,(green%nmoments-1))
371 : end if ! entropy
372 : end if ! moments
373 :
374 551 : end subroutine init_green
375 : !!***
376 :
377 : !!****f* m_green/init_green_tau
378 : !! NAME
379 : !! init_green_tau
380 : !!
381 : !! FUNCTION
382 : !! Allocate variables used in type green_type.
383 : !!
384 : !! INPUTS
385 : !! green <type(green_type)>= green function data
386 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
387 : !!
388 : !! OUTPUTS
389 : !! green <type(green_type)>= green function data
390 : !!
391 : !! SOURCE
392 :
393 37 : subroutine init_green_tau(green,paw_dmft)
394 :
395 : !Arguments ------------------------------------
396 : type(green_type), intent(inout) :: green
397 : type(paw_dmft_type), intent(in) :: paw_dmft
398 : !Local variables ------------------------------------
399 : integer :: itau,optksloc
400 : !************************************************************************
401 :
402 37 : green%use_oper_tau_ks = 0
403 37 : optksloc = 2
404 :
405 37 : green%dmftqmc_l = paw_dmft%dmftqmc_l
406 111 : ABI_MALLOC(green%tau,(green%dmftqmc_l))
407 7779 : do itau=1,green%dmftqmc_l
408 7779 : green%tau(itau) = dble(itau-1) / dble(green%dmftqmc_l) / paw_dmft%temp
409 : end do
410 :
411 37 : call init_oper(paw_dmft,green%occup_tau,opt_ksloc=optksloc)
412 :
413 7853 : ABI_MALLOC(green%oper_tau,(paw_dmft%dmftqmc_l))
414 7779 : do itau=1,green%dmftqmc_l
415 7779 : call init_oper(paw_dmft,green%oper_tau(itau),opt_ksloc=optksloc)
416 : end do
417 :
418 37 : end subroutine init_green_tau
419 : !!***
420 :
421 : !!****f* m_green/destroy_green
422 : !! NAME
423 : !! destroy_green
424 : !!
425 : !! FUNCTION
426 : !! Deallocate green
427 : !!
428 : !! INPUTS
429 : !! green <type(green_type)>= green function data
430 : !!
431 : !! OUTPUT
432 : !!
433 : !! SOURCE
434 :
435 551 : subroutine destroy_green(green)
436 :
437 : !Arguments ------------------------------------
438 : type(green_type), intent(inout) :: green
439 : !Local variables-------------------------------
440 : integer :: i,ifreq
441 : ! *********************************************************************
442 :
443 551 : call destroy_oper(green%occup)
444 551 : if (allocated(green%oper)) then
445 131124 : do ifreq=1,green%nw
446 131124 : call destroy_oper(green%oper(ifreq))
447 : end do ! ifreq
448 131124 : ABI_FREE(green%oper)
449 : end if ! green%oper
450 :
451 551 : ABI_SFREE(green%charge_matlu)
452 551 : green%has_charge_matlu = 0
453 :
454 551 : ABI_SFREE(green%charge_matlu_prev)
455 551 : green%has_charge_matlu_prev = 0
456 :
457 551 : ABI_SFREE(green%charge_matlu_solver)
458 551 : green%has_charge_matlu_solver = 0
459 :
460 551 : ABI_SFREE(green%trace_moments_log_ks)
461 551 : ABI_SFREE(green%trace_moments_log_loc)
462 551 : ABI_SFREE(green%ecorr_qmc)
463 :
464 551 : if (allocated(green%moments)) then
465 0 : do i=1,green%nmoments
466 0 : call destroy_oper(green%moments(i))
467 : end do ! i
468 0 : ABI_FREE(green%moments)
469 : end if ! green%moments
470 :
471 : !ABI_SFREE(green%procb)
472 : !ABI_SFREE(green%proct)
473 551 : green%distrib => null()
474 551 : green%omega => null()
475 :
476 551 : end subroutine destroy_green
477 : !!***
478 :
479 : !!****f* m_green/destroy_green_tau
480 : !! NAME
481 : !! destroy_green_tau
482 : !!
483 : !! FUNCTION
484 : !! Deallocate green
485 : !!
486 : !! INPUTS
487 : !! green <type(green_type)>= green function data
488 : !!
489 : !! OUTPUT
490 : !!
491 : !! SOURCE
492 :
493 37 : subroutine destroy_green_tau(green)
494 :
495 : !Arguments ------------------------------------
496 : type(green_type), intent(inout) :: green
497 : ! integer, optional, intent(in) :: opt_ksloc
498 : !Local variables-------------------------------
499 : integer :: itau
500 : ! integer :: optksloc
501 : ! *********************************************************************
502 : ! if(present(opt_ksloc)) then
503 : ! optksloc=opt_ksloc
504 : ! else
505 : ! optksloc=3
506 : ! endif
507 :
508 37 : call destroy_oper(green%occup_tau)
509 37 : if (allocated(green%oper_tau)) then
510 7779 : do itau=1,green%dmftqmc_l
511 7779 : call destroy_oper(green%oper_tau(itau))
512 : end do ! itau
513 7779 : ABI_FREE(green%oper_tau)
514 : end if ! green%oper_tau
515 37 : ABI_SFREE(green%tau)
516 :
517 37 : end subroutine destroy_green_tau
518 : !!***
519 :
520 : !!****f* m_green/copy_green
521 : !! NAME
522 : !! copy_green
523 : !!
524 : !! FUNCTION
525 : !! Copy one data structure green1 into green2
526 : !!
527 : !! INPUTS
528 : !! green1 <type(green_type)>= green function data
529 : !! green2 <type(green_type)>= green function data
530 : !! opt_tw = option to specify which data to copy
531 : !! 1: copy only green%occup_tau and green%oper_tau data
532 : !! 2: copy only green%occup and green%oper data (frequency)
533 : !!
534 : !! OUTPUT
535 : !!
536 : !! SOURCE
537 :
538 36 : subroutine copy_green(green1,green2,opt_tw)
539 :
540 : !Arguments ------------------------------------
541 : type(green_type), intent(in) :: green1
542 : type(green_type), intent(inout) :: green2
543 : integer, intent(in) :: opt_tw
544 : !Local variables-------------------------------
545 : integer :: i,ifreq,itau
546 : ! *********************************************************************
547 :
548 36 : if (opt_tw == 2) then
549 36 : call copy_oper(green1%occup,green2%occup)
550 18029 : do ifreq=1,green1%nw
551 18029 : call copy_oper(green1%oper(ifreq),green2%oper(ifreq))
552 : end do
553 36 : if (green1%has_moments == 1) then
554 0 : do i=1,green1%nmoments
555 0 : call copy_oper(green1%moments(i),green2%moments(i))
556 : end do
557 : end if
558 : ! Indicate to integrate_green that xsum has been done
559 : ! for matlu in compute_green.
560 36 : if (green1%has_greenmatlu_xsum == 1) green2%has_greenmatlu_xsum = 1
561 0 : else if (opt_tw == 1) then
562 0 : call copy_oper(green1%occup_tau,green2%occup_tau)
563 0 : do itau=1,green1%dmftqmc_l
564 0 : call copy_oper(green1%oper_tau(itau),green2%oper_tau(itau))
565 : end do ! itau
566 : end if ! opt_tw
567 :
568 36 : end subroutine copy_green
569 : !!***
570 :
571 : !!****f* m_green/printocc_green
572 : !! NAME
573 : !! printocc_green
574 : !!
575 : !! FUNCTION
576 : !! Print occupations
577 : !!
578 : !! INPUTS
579 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
580 : !! green <type(green_type)>= green function data
581 : !! option= 1 :for G(w)
582 : !! 2 :for G(tau)
583 : !! 3 :for G(tau) and check % G(w)
584 : !! 4
585 : !! <5: write diagonal part of KS occupation matrix
586 : !! 5: for G(w)
587 : !! 6: for G(tau)
588 : !! 7 :for G(tau) and check % G(w)
589 : !! >8: write all elements of KS occup. matrix.
590 : !! 9: for G(w)
591 : !! pawprtvol: flag for print
592 : !! opt_weissgreen = 1 for Weiss field
593 : !! = 2 (default) for regular Green's function
594 : !! chtype: to specify the type of occupations being printed
595 : !!
596 : !! OUTPUT
597 : !!
598 : !! SOURCE
599 :
600 583 : subroutine printocc_green(green,option,paw_dmft,pawprtvol,opt_weissgreen,chtype)
601 :
602 : !Arguments ------------------------------------
603 : type(paw_dmft_type), intent(in) :: paw_dmft
604 : type(green_type), intent(in) :: green
605 : integer, intent(in) :: option,pawprtvol
606 : integer, optional, intent(in) :: opt_weissgreen
607 : character(len=*), optional, intent(in) :: chtype
608 : !Local variables-------------------------------
609 : character(len=500) :: message
610 : integer :: i_tau,optweissgreen
611 : ! *********************************************************************
612 :
613 583 : optweissgreen = 2
614 583 : if (present(opt_weissgreen)) optweissgreen = opt_weissgreen
615 :
616 583 : if (mod(option,4) == 1) then
617 582 : if (optweissgreen == 2) then
618 582 : if (present(chtype)) then
619 478 : write(message,'(4a)') ch10," == The ",trim(chtype)," occupations are == "
620 : else
621 104 : write(message,'(2a)') ch10," == The occupations (integral of the Green's function) are == "
622 : end if ! present(chtype)
623 0 : else if (optweissgreen == 1) then
624 0 : write(message,'(2a)') ch10," == The integrals of the Weiss function are == "
625 : end if ! optweissgreen
626 582 : call wrtout(std_out,message,'COLL')
627 582 : call print_oper(green%occup,option,paw_dmft,pawprtvol)
628 : end if ! mod(option,4)=1
629 :
630 583 : if (mod(option,4) >= 2) then
631 1 : if (optweissgreen == 2) then
632 1 : write(message,'(2a)') ch10," == The occupations (value of G(tau) for tau=0-) are == "
633 0 : else if (optweissgreen == 1) then
634 0 : write(message,'(2a)') ch10," == Values of G_0(tau) for tau=0- are == "
635 : end if ! optweissgreen
636 1 : call wrtout(std_out,message,'COLL')
637 1 : call print_oper(green%occup_tau,option,paw_dmft,pawprtvol)
638 : ! write(message,'(2a)') ch10," == check: occupations from Green functions are == "
639 : ! call wrtout(std_out,message,'COLL')
640 : ! call print_oper(green%occup,1,paw_dmft,pawprtvol)
641 1 : if (mod(option,4) >= 3) then
642 : call diff_matlu("Local occup from integral of G(iw) ","Local occup from G(tau=0-) ",&
643 0 : & green%occup%matlu(:),green%occup_tau%matlu(:),paw_dmft%natom,1,tol4)
644 0 : write(message,'(2a)') ch10,&
645 0 : & ' ***** => Calculations of occupations in omega and tau spaces are coherent ****'
646 0 : call wrtout(std_out,message,'COLL')
647 : end if ! mod(option,4)>=3
648 : end if ! mod(option,4)>=2
649 :
650 583 : if (present(chtype)) then
651 : if (paw_dmft%prtvol >= 4 .and. &
652 : & (chtype == "DFT+DMFT (end of DMFT loop)" .or. chtype == "converged DMFT") &
653 478 : & .and. green%occup%has_opermatlu == 1) then
654 87 : write(message,'(4a)') ch10," == The DFT+DMFT occupation matrix for correlated electrons is == "
655 87 : call wrtout(ab_out,message,'COLL')
656 87 : call print_matlu(green%occup%matlu(:),paw_dmft%natom,pawprtvol,opt_ab_out=1)
657 87 : write(message,'(a)') " "
658 87 : call wrtout(ab_out,message,'COLL')
659 : end if
660 : end if ! present(chtype)
661 :
662 583 : if (mod(option,4) >= 2) then
663 1 : i_tau = 1
664 1 : if (optweissgreen == 1) i_tau = -1
665 1 : call trace_matlu(green%occup_tau%matlu(:),paw_dmft%natom,itau=i_tau)
666 : end if ! mod(option,4)>=2
667 :
668 583 : end subroutine printocc_green
669 : !!***
670 :
671 : !!****f* m_green/print_green
672 : !! NAME
673 : !! print_green
674 : !!
675 : !! FUNCTION
676 : !! print green function
677 : !!
678 : !! INPUTS
679 : !! char1 = character which describes the type of green function
680 : !! green <type(green_type)>= green function data
681 : !! option=1 print local green function
682 : !! 2 print KS green function
683 : !! 3 print both local and KS green function
684 : !! 4 print spectral function is green%w_type="real"
685 : !! 5 print k-resolved spectral function is green%w_type="real"
686 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
687 : !! pawprtvol = printing option
688 : !! opt_wt=1 print green function as a function of frequency
689 : !! 2 print green function as a function of imaginary time
690 : !! opt_decim= if present, write more decimals
691 : !!
692 : !! OUTPUT
693 : !!
694 : !! SOURCE
695 :
696 1 : subroutine print_green(char1,green,option,paw_dmft,opt_wt,opt_decim)
697 :
698 : !Arguments ------------------------------------
699 : type(paw_dmft_type), intent(in) :: paw_dmft
700 : type(green_type), intent(inout) :: green
701 : integer, intent(in) :: option
702 : integer, optional, intent(in) :: opt_wt,opt_decim
703 : character(len=*), intent(in) :: char1
704 : !Local variables-------------------------------
705 : integer :: iall,iatom,ib,ifreq,ikpt,im,ispinor
706 : integer :: isppol,itau,lpawu,lsub,mbandc,natom
707 : integer :: ndim,nkpt,nspinor,nsppol,optwt,spf_unt
708 : integer :: spfkresolved_unt,spcorb_unt
709 : !real(dp) :: ima,re
710 : character(len=2000) :: message
711 : character(len=fnlen) :: tmpfil
712 : character(len=1) :: tag_is,tag_is2
713 : character(len=10) :: tag_at
714 : character(len=3) :: tag_ik
715 1 : integer, allocatable :: unitgreenfunc_arr(:),unitgreenloc_arr(:)
716 1 : complex(dp), allocatable :: sf(:,:),sf_corr(:),sf2(:)
717 : ! *********************************************************************
718 :
719 1 : optwt = 1
720 1 : if (present(opt_wt)) optwt = opt_wt
721 :
722 1 : mbandc = paw_dmft%mbandc
723 1 : natom = paw_dmft%natom
724 1 : nkpt = paw_dmft%nkpt
725 1 : nspinor = paw_dmft%nspinor
726 1 : nsppol = paw_dmft%nsppol
727 :
728 : ! == Print local Green Function
729 1 : if (option == 1 .or. option == 3) then
730 0 : ABI_MALLOC(unitgreenfunc_arr,(natom*nsppol*nspinor))
731 0 : iall = 0
732 0 : do iatom=1,natom
733 0 : lpawu = paw_dmft%lpawu(iatom)
734 0 : if (lpawu == -1) cycle
735 0 : ndim = 2*lpawu + 1
736 0 : call int2char4(iatom,tag_at)
737 0 : ABI_CHECK((tag_at(1:1)/='#'),'Bug: string length too short!')
738 0 : do isppol=1,nsppol
739 0 : write(tag_is,'(i1)') isppol
740 0 : do ispinor=1,nspinor
741 0 : iall = iall + 1
742 0 : write(tag_is2,'(i1)') ispinor
743 : ! == Create names
744 0 : if (optwt == 1) then
745 : tmpfil = trim(paw_dmft%filapp)//'Green-'//trim(char1)//'-omega_iatom'// &
746 0 : & trim(tag_at)//'_isppol'//tag_is//'_ispinor'//tag_is2
747 : else
748 : tmpfil = trim(paw_dmft%filapp)//'Green-'//trim(char1)//'-tau_iatom'// &
749 0 : & trim(tag_at)//'_isppol'//tag_is//'_ispinor'//tag_is2
750 : end if ! optwt
751 0 : if (iall <= 4) then
752 0 : write(message,'(3a)') ch10," == Print green function on file ",trim(tmpfil)
753 0 : call wrtout(std_out,message,'COLL')
754 0 : else if (iall == 5) then
755 0 : write(message,'(3a)') ch10," == following values are printed in files"
756 0 : call wrtout(std_out,message,'COLL')
757 : end if ! iall
758 0 : unitgreenfunc_arr(iall) = 300 + iall - 1
759 0 : if ((optwt == 1 .or. green%fileprt_tau == 0) .or. (optwt == 2 .and. green%fileprt_tau == 1)) &
760 0 : & open(unit=unitgreenfunc_arr(iall),file=trim(tmpfil),status='unknown',form='formatted',position='append')
761 :
762 : ! == Write in files
763 : ! rewind(unitgreenfunc_arr(iall))
764 : ! write(message,'(a,a,a,i4)') 'opened file : ', trim(tmpfil), ' unit', unitgreenfunc_arr(iall)
765 : ! call wrtout(std_out,message,'COLL')
766 0 : write(message,'(2a)') ch10,"# New record :"
767 0 : call wrtout(unitgreenfunc_arr(iall),message,'COLL')
768 0 : if (optwt == 1) then
769 0 : do ifreq=1,green%nw
770 0 : if (present(opt_decim)) then
771 0 : write(message,'(2x,30(e23.16,2x))') green%omega(ifreq), &
772 0 : & (green%oper(ifreq)%matlu(iatom)%mat(im+(ispinor-1)*ndim,im+(ispinor-1)*ndim,isppol),im=1,ndim)
773 : else
774 0 : write(message,'(2x,30(e10.3,2x))') green%omega(ifreq), &
775 0 : & (green%oper(ifreq)%matlu(iatom)%mat(im+(ispinor-1)*ndim,im+(ispinor-1)*ndim,isppol),im=1,ndim)
776 : end if ! opt_decim
777 0 : call wrtout(unitgreenfunc_arr(iall),message,'COLL')
778 : !re=real(green%oper(ifreq)%matlu(iatom)%mat(1,1,isppol,ispinor,ispinor))
779 : !ima=aimag(green%oper(ifreq)%matlu(iatom)%mat(1,1,isppol,ispinor,ispinor))
780 : ! write(228,*) green%omega(ifreq),re/(re**2+ima**2),ima/(re**2+ima**2)+green%omega(ifreq)
781 : end do ! ifreq
782 : else
783 0 : do itau=1,green%dmftqmc_l
784 0 : write(message,'(2x,30(e10.3,2x))') green%tau(itau), &
785 0 : & (green%oper_tau(itau)%matlu(iatom)%mat(im+(ispinor-1)*ndim,im+(ispinor-1)*ndim,isppol),im=1,ndim)
786 0 : call wrtout(unitgreenfunc_arr(iall),message,'COLL')
787 : end do ! itau
788 0 : write(message,'(2x,30(e10.3,2x))') one/paw_dmft%temp,&
789 0 : & (-green%oper_tau(1)%matlu(iatom)%mat(im+(ispinor-1)*ndim,im+(ispinor-1)*ndim,isppol)-one,im=1,ndim)
790 0 : call wrtout(unitgreenfunc_arr(iall),message,'COLL')
791 : end if ! optwt
792 0 : close(unitgreenfunc_arr(iall))
793 : end do ! ispinor
794 : end do ! isppol
795 : end do ! iatom
796 0 : ABI_FREE(unitgreenfunc_arr)
797 : end if ! option=1 or 3
798 :
799 : ! == Print ks green function
800 1 : if ((option == 2 .or. option == 3) .and. green%oper(1)%has_operks == 1) then
801 0 : ABI_MALLOC(unitgreenloc_arr,(nsppol*nkpt))
802 0 : iall = 0
803 0 : do isppol=1,nsppol
804 0 : write(tag_is,'(i1)') isppol
805 0 : do ikpt=1,nkpt
806 0 : write(tag_ik,'(i3)') ikpt
807 : ! do ib1 = 1, mbandc
808 0 : iall = iall + 1
809 : ! == Create names
810 0 : if (optwt == 1) then
811 0 : tmpfil = trim(paw_dmft%filapp)//'Green-'//trim(char1)//'-omega_isppol'//tag_is//'_ikpt'//trim(adjustl(tag_ik))
812 : else
813 0 : tmpfil = trim(paw_dmft%filapp)//'Green-'//trim(char1)//'-tau_isppol'//tag_is//'_ikpt'//trim(adjustl(tag_ik))
814 : end if ! optwt
815 0 : if (iall <= 4) then
816 0 : write(message,'(3a)') ch10," == Print green function on file ",trim(tmpfil)
817 0 : call wrtout(std_out,message,'COLL')
818 0 : else if (iall == 5) then
819 0 : write(message,'(3a)') ch10," == following values are printed in files"
820 0 : call wrtout(std_out,message,'COLL')
821 : end if ! iall
822 0 : unitgreenloc_arr(iall) = 400 + iall - 1
823 0 : open(unit=unitgreenloc_arr(iall),file=trim(tmpfil),status='unknown',form='formatted')
824 : ! rewind(unitgreenloc_arr(iall))
825 : ! write(message,'(a,a,a,i4)') 'opened file : ', trim(tmpfil), ' unit', unitgreenloc_arr(iall)
826 : ! call wrtout(std_out,message,'COLL')
827 :
828 : ! == Write in files
829 0 : write(message,'(2a)') ch10,"# New record : First 20 bands"
830 0 : call wrtout(unitgreenloc_arr(iall),message,'COLL')
831 : ! call flush(std_out)
832 0 : do lsub=1,mbandc/20+1
833 0 : if (optwt == 1) then
834 0 : do ifreq=1,green%nw
835 : ! call flush(std_out)
836 0 : write(message,'(2x,50(e10.3,2x))') green%omega(ifreq), &
837 0 : & (green%oper(ifreq)%ks(ib,ib,ikpt,isppol),ib=20*(lsub-1)+1,min(20*lsub,mbandc))
838 0 : call wrtout(unitgreenloc_arr(iall),message,'COLL')
839 : !re = dble(green%oper(ifreq)%ks(1,1,ikpt,isppol))
840 : !ima = aimag(green%oper(ifreq)%ks(1,1,ikpt,isppol))
841 : !if (ikpt == 1) write(229,*) green%omega(ifreq),re/(re**2+ima**2),ima/(re**2+ima**2)+green%omega(ifreq)
842 : end do ! ifreq
843 0 : else if (green%use_oper_tau_ks == 1) then
844 0 : do itau=1,green%dmftqmc_l
845 : ! call flush(std_out)
846 0 : write(message,'(2x,50(e10.3,2x))') green%tau(itau), &
847 0 : & (green%oper_tau(itau)%ks(ib,ib,ikpt,isppol),ib=20*(lsub-1)+1,min(20*lsub,mbandc))
848 0 : call wrtout(unitgreenloc_arr(iall),message,'COLL')
849 : end do ! itau
850 : end if ! optwt
851 0 : if (20*lsub < mbandc) write(message,'(2a,i5,a,i5)') ch10,"# Same record, Following bands : From ", &
852 0 : & 20*(lsub)," to ",min(20*(lsub+1),mbandc)
853 0 : call wrtout(unitgreenloc_arr(iall),message,'COLL')
854 : end do ! lsub
855 0 : close(unitgreenloc_arr(iall))
856 : end do ! ikpt
857 : end do ! isppol
858 0 : ABI_FREE(unitgreenloc_arr)
859 : end if ! option=2 or 3
860 :
861 1 : if ((green%w_type == "real" .and. option >= 4) .and. green%oper(1)%has_operks == 1) then
862 1 : write(message,'(2a)') ch10," == About to print spectral function"
863 1 : call wrtout(std_out,message,'COLL')
864 1 : if (option == 4) then
865 0 : tmpfil = trim(paw_dmft%filapp)//'_SpFunc-'//trim(char1)
866 0 : if (open_file(tmpfil,message,newunit=spf_unt,status='unknown',form='formatted') /= 0) &
867 0 : & ABI_ERROR(message)
868 0 : write(spf_unt,'(3a)') "# This is the total spectral function (DOS).",ch10, &
869 0 : & "# Real frequency (Ha) Spectral function"
870 :
871 : end if ! option=4
872 1 : if (option == 5) then
873 1 : tmpfil = trim(paw_dmft%filapp)//'_DFTDMFT_SpectralFunction_kres' !//trim(char1)
874 1 : if (open_file(tmpfil,message,newunit=spfkresolved_unt,status='unknown',form='formatted') /= 0) &
875 0 : & ABI_ERROR(message)
876 1 : write(spfkresolved_unt,'(3a)') "# This is the k-resolved spectral function.",ch10, &
877 2 : & "# Real frequency (eV) Spectral function (eV^-1) ikpt"
878 4 : ABI_MALLOC(sf,(nkpt,green%nw))
879 36001 : sf(:,:) = czero
880 3001 : do ifreq=1,green%nw
881 6001 : do isppol=1,nsppol
882 39000 : do ikpt=1,nkpt
883 333000 : do ib=1,mbandc
884 330000 : sf(ikpt,ifreq) = sf(ikpt,ifreq) + green%oper(ifreq)%ks(ib,ib,ikpt,isppol)
885 : end do ! ib
886 : end do ! ikpt
887 : end do ! isppol
888 : end do ! ifreq
889 12 : do ikpt=1,nkpt
890 33011 : do ifreq=1,green%nw
891 33000 : write(message,'(2x,2(es24.16e3,2x),i5)') green%omega(ifreq)*Ha_eV,(-aimag(sf(ikpt,ifreq)))/pi/Ha_eV,ikpt
892 33011 : call wrtout(spfkresolved_unt,message,'COLL')
893 : end do ! ifreq
894 11 : write(message,*)
895 12 : call wrtout(spfkresolved_unt,message,'COLL')
896 : end do ! ikpt
897 1 : write(message,*) ch10
898 1 : call wrtout(spfkresolved_unt,message,'COLL')
899 1 : ABI_FREE(sf)
900 1 : close(spfkresolved_unt)
901 : !
902 : ! do isppol = 1 , nsppol
903 : ! do ikpt = 1, nkpt
904 : ! do ib=1,mbandc
905 : ! sf=czero
906 : ! write(71,*)
907 : ! write(71,*) "#", ikpt, ib
908 : ! do ifreq=1,green%nw
909 : ! sf(ifreq)=sf(ifreq)+green%oper(ifreq)%ks(isppol,ikpt,ib,ib)
910 : ! write(71,*) green%omega(ifreq)*Ha_eV,(-aimag(sf(ifreq)))/pi/Ha_eV,ikpt
911 : ! enddo
912 : ! enddo
913 : ! enddo
914 : ! enddo
915 : end if ! option=5
916 :
917 1 : if (option == 4) then
918 0 : ABI_MALLOC(sf2,(green%nw))
919 0 : sf2(:) = czero
920 0 : do ifreq=1,green%nw
921 0 : do isppol=1,nsppol
922 0 : do ikpt=1,nkpt
923 0 : do ib=1,mbandc
924 0 : sf2(ifreq) = sf2(ifreq) + green%oper(ifreq)%ks(ib,ib,ikpt,isppol)*green%oper(1)%wtk(ikpt)
925 : end do ! ib
926 : end do ! ikpt
927 : end do ! isppol
928 : end do ! ifreq
929 0 : do ifreq=1,green%nw
930 0 : write(message,'(2x,2(es24.16e3,2x))') green%omega(ifreq),(-aimag(sf2(ifreq)))/pi
931 0 : call wrtout(spf_unt,message,'COLL')
932 : end do ! ifreq
933 0 : ABI_FREE(sf2)
934 0 : close(spf_unt)
935 : end if ! option=4
936 :
937 1 : if (paw_dmft%dmft_kspectralfunc == 1) then
938 3 : ABI_MALLOC(sf_corr,(green%nw))
939 2 : do iatom=1,natom
940 1 : lpawu = paw_dmft%lpawu(iatom)
941 1 : if (lpawu == -1) cycle
942 3001 : sf_corr(:) = czero
943 1 : call int2char4(iatom,tag_at)
944 1 : ABI_CHECK((tag_at(1:1)/='#'),'Bug: string length too short!')
945 1 : tmpfil = trim(paw_dmft%filapp)//'_DFTDMFT_SpFunloc_iatom'//trim(tag_at)
946 1 : if (open_file(tmpfil,message,newunit=spcorb_unt,status='unknown',form='formatted') /= 0) &
947 0 : & ABI_ERROR(message)
948 1 : ndim = 2*lpawu + 1
949 1 : write(message,'(3a,3(i3),i6,2a,i2,2a)') "# nspinor,nsppol,ndim,nw",ch10,"#",nspinor,nsppol,ndim,green%nw,ch10, &
950 2 : & "# lpawu",lpawu,ch10,"# Frequency (Ha.) Spectral function"
951 1 : call wrtout(spcorb_unt,message,'COLL')
952 1 : ndim = nspinor * ndim
953 3001 : do ifreq=1,green%nw
954 6001 : do isppol=1,nsppol
955 21000 : do im=1,ndim
956 18000 : sf_corr(ifreq) = sf_corr(ifreq) + green%oper(ifreq)%matlu(iatom)%mat(im,im,isppol)
957 : end do ! im
958 : end do ! isppol
959 : end do ! ifreq
960 3001 : do ifreq=1,green%nw
961 3000 : write(message,'(2x,2(es24.16e3,2x))') green%omega(ifreq),(-aimag(sf_corr(ifreq)))/pi
962 3001 : call wrtout(spcorb_unt,message,'COLL')
963 : end do ! ifreq
964 2 : close(spcorb_unt)
965 : end do ! iatom
966 1 : ABI_FREE(sf_corr)
967 : end if ! dmft_kspectralfunc=1
968 : end if ! (green%w_type == "real" .and. option >= 4) .and. green%oper(1)%has_operks == 1
969 :
970 1 : if (optwt == 2 .and. (option == 1 .or. option == 3)) green%fileprt_tau = 1 ! file for G(tau) has been created here
971 :
972 1 : end subroutine print_green
973 : !!***
974 :
975 : !!****f* m_green/compute_green_batched_core
976 : !! NAME
977 : !! compute_green_batched_core
978 : !!
979 : !! FUNCTION
980 : !! Variant for core loop of compute green function, useful for computation on GPU
981 : !!
982 : !! INPUTS
983 : !! green <type(green_type)>= green function data
984 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
985 : !! prtopt : option for printing
986 : !! self <type(self_type)>= variables related to self-energy
987 : !! opt_self = optional argument, if =1, upfold self-energy
988 : !! opt_log = if 1, compute Tr(log(G))
989 : !!
990 : !! OUTPUT
991 : !!
992 : !! SOURCE
993 0 : subroutine compute_green_batched_core(green,paw_dmft,self,optself,optlog)
994 :
995 : use m_abi_linalg, only : abi_xgemm
996 : use m_matlu, only : add_matlu,sym_matlu
997 : use m_oper, only : downfold_oper,inverse_oper,upfold_oper
998 : use m_time, only : timab
999 :
1000 : !Arguments ------------------------------------
1001 : type(green_type), target, intent(inout) :: green
1002 : type(paw_dmft_type), target, intent(in) :: paw_dmft
1003 : type(self_type), intent(inout) :: self
1004 : integer, intent(in) :: optlog,optself
1005 : !Local variables-------------------------------
1006 : logical :: oper_ndat_allocated
1007 : integer :: diag,ib,ib1,ifreq,ikpt,info,isppol,lwork,mbandc,ifreq_beg,ifreq_end,idat
1008 : integer :: me_kpt,mkmem,myproc,natom,nkpt,nmoments,nspinor,nsppol,gpu_option,ndat
1009 : integer :: option,shift,shift_green,spacecomm,optoper_ksloc
1010 : real(dp) :: fermilevel,wtk,temp
1011 : complex(dp) :: green_tmp,trace_tmp
1012 0 : real(dp), allocatable :: eig(:),rwork(:),fac(:)
1013 0 : complex(dp), allocatable :: mat_tmp(:,:),work(:),omega_current(:)
1014 0 : type(oper_type), target :: green_oper_ndat
1015 0 : real(dp), ABI_CONTIGUOUS pointer :: eigen_dft(:,:,:)
1016 0 : complex(dp), ABI_CONTIGUOUS pointer :: ks(:,:,:,:),occup_ks(:,:,:,:)
1017 : ! *********************************************************************
1018 :
1019 : ABI_NVTX_START_RANGE(NVTX_DMFT_COMPUTE_GREEN_BATCHED)
1020 :
1021 0 : diag = 1 - optself
1022 :
1023 : ! Initialise spaceComm, myproc, and nproc
1024 0 : me_kpt = green%distrib%me_kpt
1025 0 : myproc = paw_dmft%myproc
1026 0 : spacecomm = paw_dmft%spacecomm
1027 0 : oper_ndat_allocated = .false.
1028 :
1029 : ! Initialise integers
1030 0 : mbandc = paw_dmft%mbandc
1031 0 : natom = paw_dmft%natom
1032 0 : nkpt = paw_dmft%nkpt
1033 0 : nspinor = paw_dmft%nspinor
1034 0 : nsppol = paw_dmft%nsppol
1035 0 : temp = paw_dmft%temp
1036 0 : gpu_option = paw_dmft%gpu_option
1037 :
1038 0 : if (optlog == 1) then
1039 0 : ABI_MALLOC(eig,(mbandc))
1040 0 : ABI_MALLOC(work,(2*mbandc-1))
1041 0 : ABI_MALLOC(rwork,(3*mbandc-2))
1042 0 : ABI_MALLOC(mat_tmp,(mbandc,mbandc))
1043 0 : call zheev("n","u",mbandc,mat_tmp(:,:),mbandc,eig(:),work(:),-1,rwork(:),info)
1044 0 : lwork = int(work(1))
1045 0 : ABI_FREE(work)
1046 0 : ABI_MALLOC(work,(lwork))
1047 0 : green%trace_log = zero
1048 : end if ! optlog=1
1049 :
1050 0 : fermilevel = paw_dmft%fermie
1051 :
1052 : ! option for downfold_oper
1053 0 : option = 1
1054 0 : if (diag == 1) option = 3
1055 :
1056 : ! ================================
1057 : ! == Compute Green function G(k)
1058 : ! ================================
1059 0 : green%occup%ks(:,:,:,:) = czero
1060 :
1061 0 : nmoments = 1
1062 : if (green%has_moments == 1) then
1063 : nmoments = green%nmoments
1064 : end if ! moments
1065 :
1066 0 : shift = green%distrib%shiftk
1067 0 : mkmem = green%distrib%nkpt_mem(green%distrib%me_kpt+1)
1068 0 : shift_green = shift
1069 0 : ndat = green%distrib%nw_mem_kptparal(green%distrib%me_freq+1)
1070 0 : if (green%oper(1)%has_operks == 0) shift_green = 0
1071 :
1072 0 : if(green%oper(1)%has_opermatlu==0) optoper_ksloc=1
1073 0 : if(green%oper(1)%has_opermatlu==1) optoper_ksloc=3
1074 :
1075 0 : do ifreq=1,green%nw
1076 0 : if (green%distrib%proct(ifreq) == green%distrib%me_freq) then
1077 0 : ifreq_beg = ifreq; ifreq_end = ifreq_beg + ndat - 1
1078 0 : exit
1079 : end if
1080 : end do
1081 :
1082 0 : ABI_MALLOC(omega_current,(green%nw))
1083 0 : ABI_MALLOC(fac,(green%nw))
1084 : ! Initialise for compiler
1085 0 : omega_current = czero
1086 0 : do ifreq=1,green%nw
1087 : !if(present(iii)) write(6,*) ch10,'ifreq self', ifreq,self%oper(ifreq)%matlu(1)%mat(1,1,1,1,1)
1088 : ! ====================================================
1089 : ! First Upfold self-energy and double counting Self_imp -> self(k)
1090 : ! ====================================================
1091 : ! if(mod(ifreq-1,nproc)==myproc) then
1092 : ! write(6,*) "compute_green ifreq",ifreq, mod(ifreq-1,nproc)==myproc,proct(ifreq,myproc)==1
1093 0 : if (green%distrib%proct(ifreq) /= green%distrib%me_freq) cycle
1094 0 : if (green%w_type == "imag") then
1095 0 : omega_current(ifreq) = cmplx(zero,green%omega(ifreq),kind=dp)
1096 0 : fac(ifreq) = two * paw_dmft%wgt_wlo(ifreq) * paw_dmft%temp
1097 0 : else if (green%w_type == "real") then
1098 0 : omega_current(ifreq) = cmplx(green%omega(ifreq),paw_dmft%temp,kind=dp)
1099 : end if ! green%w_type
1100 :
1101 0 : if (green%oper(ifreq)%has_operks == 0) then
1102 0 : ABI_MALLOC(green%oper(ifreq)%ks,(mbandc,mbandc,mkmem,nsppol))
1103 0 : green%oper(ifreq)%nkpt = mkmem
1104 0 : green%oper(ifreq)%paral = 1
1105 0 : green%oper(ifreq)%shiftk = shift
1106 : end if
1107 0 : if(.not. oper_ndat_allocated) then
1108 0 : call init_oper_ndat(paw_dmft,green_oper_ndat,ndat,nkpt=green%oper(ifreq)%nkpt,opt_ksloc=optoper_ksloc,gpu_option=gpu_option)
1109 0 : if (green%oper(ifreq)%has_operks == 0) then
1110 0 : green_oper_ndat%paral = 1
1111 0 : green_oper_ndat%shiftk = shift
1112 : end if
1113 : oper_ndat_allocated=.true.
1114 : end if
1115 : end do ! ifreq
1116 :
1117 0 : if (optself == 0) then
1118 0 : if(green_oper_ndat%gpu_option==ABI_GPU_DISABLED) then
1119 0 : green_oper_ndat%ks(:,:,:,:) = czero
1120 0 : else if(green_oper_ndat%gpu_option==ABI_GPU_OPENMP) then
1121 0 : call gpu_set_to_zero_complex(green_oper_ndat%ks, int(nsppol,c_size_t)*ndat*mbandc*mbandc*mkmem)
1122 : end if
1123 0 : call copy_oper_to_ndat(green%oper,green_oper_ndat,ndat,green%nw,green%distrib%proct,green%distrib%me_freq,.false.)
1124 : else
1125 0 : do ifreq=ifreq_beg,ifreq_end
1126 0 : call add_matlu(self%hdc%matlu(:),self%oper(ifreq)%matlu(:),green%oper(ifreq)%matlu(:),natom,-1)
1127 : end do
1128 0 : call copy_oper_to_ndat(green%oper,green_oper_ndat,ndat,green%nw,green%distrib%proct,green%distrib%me_freq,.false.)
1129 0 : call upfold_oper(green_oper_ndat,paw_dmft,procb=green%distrib%procb(:),iproc=me_kpt,gpu_option=gpu_option)
1130 : end if ! optself
1131 :
1132 0 : ks => green_oper_ndat%ks
1133 0 : eigen_dft => paw_dmft%eigen_dft
1134 0 : if(gpu_option==ABI_GPU_DISABLED) then
1135 0 : do ifreq=ifreq_beg,ifreq_end
1136 0 : do isppol=1,nsppol
1137 0 : do ikpt=1,mkmem
1138 0 : do ib=1,mbandc
1139 0 : idat=ifreq-(ifreq_end-ndat)
1140 0 : green_tmp = omega_current(ifreq) + fermilevel - eigen_dft(ib,ikpt+shift,isppol)
1141 0 : if (optself == 0) then
1142 0 : ks(ib,ib+(idat-1)*mbandc,ikpt+shift_green,isppol) = cone / green_tmp
1143 : else
1144 : ks(ib,ib+(idat-1)*mbandc,ikpt+shift_green,isppol) = &
1145 0 : & ks(ib,ib+(idat-1)*mbandc,ikpt+shift_green,isppol) + green_tmp
1146 : end if
1147 : end do ! ib
1148 : end do ! ikpt
1149 : end do ! isppol
1150 : end do ! ifreq
1151 : else if(gpu_option==ABI_GPU_OPENMP) then
1152 : #ifdef HAVE_OPENMP_OFFLOAD
1153 : if (optself == 0) then
1154 : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:ks,eigen_dft) PRIVATE(ifreq,idat)
1155 : do ifreq=ifreq_beg,ifreq_end
1156 : idat=ifreq-(ifreq_end-ndat)
1157 : !$OMP PARALLEL DO COLLAPSE(3) PRIVATE(isppol,ikpt,ib,green_tmp)
1158 : do isppol=1,nsppol
1159 : do ikpt=1,mkmem
1160 : do ib=1,mbandc
1161 : green_tmp = omega_current(ifreq) + fermilevel - eigen_dft(ib,ikpt+shift,isppol)
1162 : ks(ib,ib+(idat-1)*mbandc,ikpt+shift_green,isppol) = cone / green_tmp
1163 : end do ! ib
1164 : end do ! ikpt
1165 : end do ! isppol
1166 : end do ! ifreq
1167 : else
1168 : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:ks,eigen_dft) PRIVATE(ifreq,idat)
1169 : do ifreq=ifreq_beg,ifreq_end
1170 : idat=ifreq-(ifreq_end-ndat)
1171 : !$OMP PARALLEL DO COLLAPSE(3) PRIVATE(isppol,ikpt,ib,green_tmp)
1172 : do isppol=1,nsppol
1173 : do ikpt=1,mkmem
1174 : do ib=1,mbandc
1175 : green_tmp = omega_current(ifreq) + fermilevel - eigen_dft(ib,ikpt+shift,isppol)
1176 : ks(ib,ib+(idat-1)*mbandc,ikpt+shift_green,isppol) = &
1177 : & ks(ib,ib+(idat-1)*mbandc,ikpt+shift_green,isppol) + green_tmp
1178 : end do ! ib
1179 : end do ! ikpt
1180 : end do ! isppol
1181 : end do ! ifreq
1182 : end if
1183 : #endif
1184 : end if
1185 :
1186 0 : if (optself /= 0) then
1187 0 : call inverse_oper(green_oper_ndat,1,procb=green%distrib%procb(:),iproc=me_kpt,gpu_option=gpu_option)
1188 : end if
1189 :
1190 0 : if (optlog == 1) then
1191 0 : call copy_oper_from_ndat(green_oper_ndat,green%oper,ndat,green%nw,green%distrib%proct,green%distrib%me_freq,.true.)
1192 0 : do ifreq=1,green%nw
1193 0 : if (green%distrib%proct(ifreq) /= green%distrib%me_freq) cycle
1194 : trace_tmp = czero
1195 0 : do isppol=1,nsppol
1196 0 : do ikpt=1,mkmem
1197 0 : wtk = green%oper(ifreq)%wtk(ikpt+shift)
1198 0 : if (optself == 0) then
1199 0 : do ib=1,mbandc
1200 0 : trace_tmp = trace_tmp + two*temp*wtk*log(green%oper(ifreq)%ks(ib,ib,ikpt+shift_green,isppol)*omega_current(ifreq))
1201 : end do ! ib
1202 : else
1203 : ! Use Tr(log(G(iw))) + Tr(log(G(-iw))) = Tr(log(G(iw)G(iw)^H)) so that we can use the much faster zheev instead of
1204 : ! zgeev. Even if G(iw) and G(iw)^H have no reason to commute, this equality still holds since we only care about the trace.
1205 : call abi_xgemm("n","c",mbandc,mbandc,mbandc,cone,green%oper(ifreq)%ks(:,:,ikpt+shift_green,isppol), &
1206 0 : & mbandc,green%oper(ifreq)%ks(:,:,ikpt+shift_green,isppol),mbandc,czero,mat_tmp(:,:),mbandc)
1207 0 : call zheev("n","u",mbandc,mat_tmp(:,:),mbandc,eig(:),work(:),lwork,rwork(:),info)
1208 : ! Do not use DOT_PRODUCT
1209 0 : trace_tmp = trace_tmp + sum(log(eig(:)*omega_current(ifreq)))*wtk*temp
1210 : end if ! optself
1211 : end do ! ikpt
1212 : end do ! isppol
1213 0 : if (nsppol == 1 .and. nspinor == 1) trace_tmp = trace_tmp * two
1214 0 : green%trace_log = green%trace_log + dble(trace_tmp)
1215 : end do ! ifreq
1216 0 : call copy_oper_to_ndat(green%oper,green_oper_ndat,ndat,green%nw,green%distrib%proct,green%distrib%me_freq,.true.)
1217 : end if ! optlog=1
1218 :
1219 0 : if (green%w_type /= "real") then
1220 : ABI_NVTX_START_RANGE(NVTX_DMFT_ADD_INT_FCT)
1221 0 : occup_ks => green%occup%ks
1222 0 : ks => green_oper_ndat%ks
1223 0 : if(gpu_option==ABI_GPU_DISABLED) then
1224 0 : do ifreq=ifreq_beg,ifreq_end
1225 :
1226 0 : do isppol=1,nsppol
1227 0 : do ikpt=1,mkmem
1228 0 : do ib1=1,mbandc
1229 0 : idat=ifreq-(ifreq_end-ndat)
1230 0 : if (diag == 1) then
1231 : green%occup%ks(ib1,ib1,ikpt+shift,isppol) = green%occup%ks(ib1,ib1,ikpt+shift,isppol) + &
1232 0 : & fac(ifreq)*green_oper_ndat%ks(ib1,ib1+(idat-1)*mbandc,ikpt+shift_green,isppol)
1233 : else
1234 0 : do ib=1,mbandc
1235 : green%occup%ks(ib,ib1,ikpt+shift,isppol) = green%occup%ks(ib,ib1,ikpt+shift,isppol) + &
1236 0 : & fac(ifreq)*green_oper_ndat%ks(ib,ib1+(idat-1)*mbandc,ikpt+shift_green,isppol)
1237 : end do ! ib
1238 : end if ! diag
1239 : end do ! ib1
1240 : end do ! ikpt
1241 : end do ! isppol
1242 :
1243 : end do ! ifreq
1244 : else if(gpu_option==ABI_GPU_OPENMP) then
1245 : #ifdef HAVE_OPENMP_OFFLOAD
1246 : if (diag == 1) then
1247 : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(2) &
1248 : !$OMP& MAP(tofrom:occup_ks) MAP(to:ks,fac) PRIVATE(isppol,ikpt)
1249 : do isppol=1,nsppol
1250 : do ikpt=1,mkmem
1251 : !$OMP PARALLEL DO PRIVATE(ib1,ifreq,idat)
1252 : do ib1=1,mbandc
1253 : do ifreq=ifreq_beg,ifreq_end
1254 : idat=ifreq-(ifreq_end-ndat)
1255 : occup_ks(ib1,ib1,ikpt+shift,isppol) = occup_ks(ib1,ib1,ikpt+shift,isppol) + &
1256 : & fac(ifreq)*ks(ib1,ib1+(idat-1)*mbandc,ikpt+shift_green,isppol)
1257 : end do ! ifreq
1258 : end do ! ib1
1259 : end do ! ikpt
1260 : end do ! isppol
1261 :
1262 : else
1263 :
1264 : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(3) &
1265 : !$OMP& MAP(tofrom:occup_ks) MAP(to:ks,fac) PRIVATE(isppol,ikpt,ib1)
1266 : do isppol=1,nsppol
1267 : do ikpt=1,mkmem
1268 : do ib1=1,mbandc
1269 : !$OMP PARALLEL DO PRIVATE(ib,ifreq,idat)
1270 : do ib=1,mbandc
1271 : do ifreq=ifreq_beg,ifreq_end
1272 : idat=ifreq-(ifreq_end-ndat)
1273 : occup_ks(ib,ib1,ikpt+shift,isppol) = occup_ks(ib,ib1,ikpt+shift,isppol) + &
1274 : & fac(ifreq)*ks(ib,ib1+(idat-1)*mbandc,ikpt+shift_green,isppol)
1275 : end do ! ib
1276 : end do ! ib1
1277 : end do ! ikpt
1278 : end do ! isppol
1279 : end do ! ifreq
1280 :
1281 : end if ! diag
1282 : #endif
1283 : end if
1284 : ABI_NVTX_END_RANGE()
1285 : end if ! w_type/="real"
1286 :
1287 0 : if (paw_dmft%lchipsiortho == 1 .or. optself == 1) then
1288 : call downfold_oper(green_oper_ndat,paw_dmft,&
1289 0 : & procb=green%distrib%procb(:),iproc=me_kpt,option=option,gpu_option=gpu_option)
1290 : end if ! lchipsiortho=1
1291 :
1292 : call copy_oper_from_ndat(green_oper_ndat,green%oper,ndat,green%nw,green%distrib%proct,&
1293 0 : & green%distrib%me_freq,green%oper(ifreq_beg)%has_operks /= 0)
1294 :
1295 0 : do ifreq=1,green%nw
1296 0 : if (green%distrib%proct(ifreq) /= green%distrib%me_freq) cycle
1297 0 : if (green%oper(ifreq)%has_operks == 0) then
1298 0 : ABI_FREE(green%oper(ifreq)%ks)
1299 : end if
1300 : end do ! ifreq
1301 :
1302 0 : call destroy_oper(green_oper_ndat)
1303 0 : ABI_FREE(omega_current)
1304 0 : ABI_FREE(fac)
1305 :
1306 : ABI_NVTX_END_RANGE()
1307 :
1308 0 : end subroutine compute_green_batched_core
1309 : !!***
1310 :
1311 : !!****f* m_green/compute_green
1312 : !! NAME
1313 : !! compute_green
1314 : !!
1315 : !! FUNCTION
1316 : !! compute green function from DFT and self-energy
1317 : !!
1318 : !! INPUTS
1319 : !! green <type(green_type)>= green function data
1320 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
1321 : !! prtopt : option for printing
1322 : !! self <type(self_type)>= variables related to self-energy
1323 : !! opt_self = optional argument, if =1, upfold self-energy
1324 : !! opt_nonxsum = 0 : do usual xsum after calculation of green(freq)%ks
1325 : !! = 1 : do not do xsum after calculation of green(freq)%ks: each proc as only a part
1326 : !! of the data: this is useful where only the total number of electron will be computed.
1327 : !! opt_nonxsum2 0 : green(ifreq)%matlu will be broadcasted
1328 : !! 1 : green(ifreq)%matlu will not be broadcasted in compute_green: calc
1329 : !! if occupations will not possible.
1330 : !! (a keyword: compute_local_green would be in fact equivalent and more clear)
1331 : !! opt_log = if 1, compute Tr(log(G))
1332 : !! opt_restart_moments = if 1, activate quick restart for calculation of the high frequency moments
1333 : !! (useful when compute_green is called after fermi_green)
1334 : !!
1335 : !! OUTPUT
1336 : !!
1337 : !! SOURCE
1338 :
1339 4884 : subroutine compute_green(green,paw_dmft,prtopt,self,opt_self,opt_nonxsum,opt_nonxsum2,opt_log,opt_restart_moments)
1340 :
1341 : !Arguments ------------------------------------
1342 : type(green_type), intent(inout) :: green
1343 : type(paw_dmft_type), intent(in) :: paw_dmft
1344 : !type(MPI_type), intent(in) :: mpi_enreg
1345 : type(self_type), intent(inout) :: self
1346 : integer, intent(in) :: prtopt
1347 : integer, optional, intent(in) :: opt_log,opt_nonxsum,opt_nonxsum2,opt_restart_moments,opt_self
1348 : !Local variables-------------------------------
1349 : !logical :: lintegrate
1350 : integer :: band_index,diag,i,ib,ierr,ifreq,ikpt,info,isppol,lwork,mbandc
1351 : integer :: me_kpt,mkmem,myproc,natom,nband_k,nkpt,nmoments,nspinor,nsppol
1352 : integer :: opt_quick_restart,option,optlog,optnonxsum,optnonxsum2,optself
1353 : integer :: shift,shift_green,spacecomm,gpu_option
1354 : real(dp) :: beta,correction,eigen,fac,fermilevel,freq2,temp,wtk
1355 : complex(dp) :: green_tmp,omega_current,trace_tmp
1356 : character(len=500) :: message
1357 : real(dp) :: tsec(2)
1358 4884 : real(dp), allocatable :: eig(:),rwork(:)
1359 4884 : complex(dp), allocatable :: mat_tmp(:,:),omega_fac(:),work(:)
1360 : #ifdef HAVE_OPENMP_OFFLOAD
1361 : integer :: ndat
1362 : type(oper_type), target :: green_oper_ndat
1363 : #endif
1364 : ! integer, allocatable :: procb(:,:),proct(:,:)
1365 : ! *********************************************************************
1366 :
1367 : ABI_NVTX_START_RANGE(NVTX_DMFT_COMPUTE_GREEN)
1368 : !lintegrate=.true.
1369 : !if(lintegrate.and.green%w_type=="real") then
1370 : !if(green%w_type=="real") then
1371 : ! message = 'integrate_green not implemented for real frequency'
1372 : ! ABI_BUG(message)
1373 : !endif
1374 4884 : call timab(624,1,tsec(:))
1375 4884 : optself = 0
1376 4884 : if (present(opt_self)) optself = opt_self
1377 4884 : optnonxsum = 0
1378 4884 : if (present(opt_nonxsum)) optnonxsum = opt_nonxsum
1379 4884 : optnonxsum2 = 0
1380 4884 : if (present(opt_nonxsum2)) optnonxsum2 = opt_nonxsum2
1381 4884 : optlog = 0
1382 4884 : if (present(opt_log)) optlog = opt_log
1383 4884 : opt_quick_restart = 0
1384 4884 : if (present(opt_restart_moments)) opt_quick_restart = opt_restart_moments
1385 :
1386 4884 : diag = 1 - optself
1387 :
1388 4884 : if (prtopt > 0) then
1389 434 : write(message,'(2a)') ch10," === Compute Green's function "
1390 434 : call wrtout(std_out,message,'COLL')
1391 : end if ! prtopt>0
1392 :
1393 4884 : if (self%nw /= green%nw) then
1394 0 : message = ' BUG: frequencies for green and self not coherent'
1395 0 : ABI_BUG(message)
1396 : end if
1397 :
1398 : ! Initialise spaceComm, myproc, and nproc
1399 4884 : me_kpt = green%distrib%me_kpt
1400 4884 : myproc = paw_dmft%myproc
1401 : !nproc = paw_dmft%nproc
1402 4884 : spacecomm = paw_dmft%spacecomm
1403 :
1404 : ! Initialise integers
1405 : !mband = paw_dmft%mband
1406 4884 : mbandc = paw_dmft%mbandc
1407 4884 : natom = paw_dmft%natom
1408 4884 : nkpt = paw_dmft%nkpt
1409 4884 : nspinor = paw_dmft%nspinor
1410 4884 : nsppol = paw_dmft%nsppol
1411 4884 : temp = paw_dmft%temp
1412 4884 : gpu_option = paw_dmft%gpu_option
1413 :
1414 4884 : if (optlog == 1) then
1415 0 : ABI_MALLOC(eig,(mbandc))
1416 0 : ABI_MALLOC(work,(2*mbandc-1))
1417 0 : ABI_MALLOC(rwork,(3*mbandc-2))
1418 0 : ABI_MALLOC(mat_tmp,(mbandc,mbandc))
1419 0 : call zheev("n","u",mbandc,mat_tmp(:,:),mbandc,eig(:),work(:),-1,rwork(:),info)
1420 0 : lwork = int(work(1))
1421 0 : ABI_FREE(work)
1422 0 : ABI_MALLOC(work,(lwork))
1423 0 : green%trace_log = zero
1424 : end if ! optlog=1
1425 :
1426 : !icomp_chloc = 0
1427 :
1428 : option = 1
1429 4884 : fermilevel = paw_dmft%fermie
1430 : if (option == 123) then
1431 : fermilevel = two
1432 : write(message,'(2a,e14.3,a)') ch10,' Warning (special case for check: fermi level=',fermilevel,')'
1433 : call wrtout(std_out,message,'COLL')
1434 : end if ! option=123
1435 :
1436 4884 : option = merge(3,1,diag==1) ! option for downfold_oper
1437 :
1438 : ! ====================================================
1439 : ! Upfold self-energy and double counting Self_imp -> self(k)
1440 : ! ====================================================
1441 : ! if(optself==1) then
1442 : ! do ifreq=1,green%nw
1443 : ! call upfold_oper(self%oper(ifreq),paw_dmft,1)
1444 : ! enddo ! ifreq
1445 : ! call upfold_oper(self%hdc,paw_dmft,1)
1446 : ! endif
1447 :
1448 : ! ================================
1449 : ! == Compute Green function G(k)
1450 : ! ================================
1451 7124010 : green%occup%ks(:,:,:,:) = czero
1452 :
1453 4884 : nmoments = 1
1454 4884 : if (green%has_moments == 1) then
1455 0 : nmoments = green%nmoments
1456 : call compute_moments_ks(green,self,paw_dmft,opt_self=optself,opt_log=optlog, &
1457 0 : & opt_quick_restart=opt_quick_restart)
1458 0 : if (paw_dmft%lchipsiortho == 1 .or. optself == 1) then
1459 0 : call downfold_oper(green%moments(1),paw_dmft,option=2)
1460 0 : do i=2,green%nmoments
1461 0 : call downfold_oper(green%moments(i),paw_dmft,option=option)
1462 : end do ! i
1463 0 : do i=1,green%nmoments
1464 0 : call xmpi_matlu(green%moments(i)%matlu(:),natom,green%distrib%comm_kpt)
1465 0 : call sym_matlu(green%moments(i)%matlu(:),paw_dmft)
1466 : end do ! i
1467 : end if ! lchipsiortho=1
1468 : end if ! moments
1469 :
1470 4884 : if (green%w_type /= "real") then
1471 14640 : ABI_MALLOC(omega_fac,(nmoments))
1472 9760 : do i=1,nmoments
1473 4880 : omega_fac(i) = czero
1474 991739 : do ifreq=green%nw,1,-1 ! NEVER change the summation order and DON'T use the intrinsic SUM
1475 991739 : omega_fac(i) = omega_fac(i) + paw_dmft%wgt_wlo(ifreq) / (paw_dmft%omega_lo(ifreq))**i
1476 : end do
1477 4880 : omega_fac(i) = - two * temp * omega_fac(i) / (j_dpc)**i
1478 4880 : if (i == 1) omega_fac(i) = omega_fac(i) + half
1479 4880 : if (i == 2) omega_fac(i) = omega_fac(i) - cone/(four*temp)
1480 9760 : if (i == 4) omega_fac(i) = omega_fac(i) + cone/(dble(48)*(temp**3))
1481 : end do ! i
1482 : end if ! w_type
1483 :
1484 4884 : shift = green%distrib%shiftk
1485 4884 : mkmem = green%distrib%nkpt_mem(green%distrib%me_kpt+1)
1486 4884 : shift_green = merge(0,shift,green%oper(1)%has_operks==0)
1487 :
1488 :
1489 4884 : if(mkmem>0 .and. gpu_option==ABI_GPU_OPENMP) then
1490 0 : call compute_green_batched_core(green,paw_dmft,self,optself,optlog)
1491 4884 : else if(mkmem>0) then
1492 : ABI_NVTX_START_RANGE(NVTX_DMFT_COMPUTE_GREEN_LOOP)
1493 : ! Initialise for compiler
1494 4763 : omega_current = czero
1495 999750 : do ifreq=1,green%nw
1496 : !if(present(iii)) write(6,*) ch10,'ifreq self', ifreq,self%oper(ifreq)%matlu(1)%mat(1,1,1,1,1)
1497 : ! ====================================================
1498 : ! First Upfold self-energy and double counting Self_imp -> self(k)
1499 : ! ====================================================
1500 : ! if(mod(ifreq-1,nproc)==myproc) then
1501 : ! write(6,*) "compute_green ifreq",ifreq, mod(ifreq-1,nproc)==myproc,proct(ifreq,myproc)==1
1502 994987 : if (green%distrib%proct(ifreq) /= green%distrib%me_freq) cycle
1503 994987 : if (green%w_type == "imag") then
1504 982987 : omega_current = cmplx(zero,green%omega(ifreq),kind=dp)
1505 982987 : fac = two * paw_dmft%wgt_wlo(ifreq) * temp
1506 12000 : else if (green%w_type == "real") then
1507 12000 : omega_current = cmplx(green%omega(ifreq),paw_dmft%temp,kind=dp)
1508 : end if ! green%w_type
1509 :
1510 994987 : if (green%oper(ifreq)%has_operks == 0) then
1511 5897922 : ABI_MALLOC(green%oper(ifreq)%ks,(mbandc,mbandc,mkmem,nsppol))
1512 982987 : green%oper(ifreq)%nkpt = mkmem
1513 982987 : green%oper(ifreq)%paral = 1
1514 982987 : green%oper(ifreq)%shiftk = shift
1515 : end if
1516 :
1517 994987 : if (optself == 0) then
1518 24492174 : green%oper(ifreq)%ks(:,:,:,:) = czero
1519 : else
1520 949905 : call add_matlu(self%hdc%matlu(:),self%oper(ifreq)%matlu(:),green%oper(ifreq)%matlu(:),natom,-1)
1521 : ! do iatom = 1 , natom
1522 : !write(6,*) 'self matlu', ifreq, self%oper(ifreq)%matlu(1)%mat(1,1,1,1,1)
1523 : !write(6,*) 'self hdc ', ifreq, self%hdc%matlu(1)%mat(1,1,1,1,1)
1524 : !write(6,*) 'self_minus_hdc_oper ', ifreq, self_minus_hdc_oper%matlu(1)%mat(1,1,1,1,1)
1525 : ! enddo ! natom
1526 : !if(paw_dmft%dmft_solv==4) then
1527 : ! call shift_matlu(self_minus_hdc_oper%matlu,paw_dmft%natom,cmplx(self%qmc_shift,0.d0,kind=dp),-1)
1528 : ! call shift_matlu(self_minus_hdc_oper%matlu,paw_dmft%natom,cmplx(self%qmc_xmu,0.d0,kind=dp),-1)
1529 : !endif
1530 949905 : call upfold_oper(green%oper(ifreq),paw_dmft,procb=green%distrib%procb(:),iproc=me_kpt)
1531 : end if ! optself
1532 :
1533 2670781 : do isppol=1,nsppol
1534 11483801 : do ikpt=1,mkmem
1535 68181838 : do ib=1,mbandc
1536 57693024 : green_tmp = omega_current + fermilevel - paw_dmft%eigen_dft(ib,ikpt+shift,isppol)
1537 66506044 : if (optself == 0) then
1538 2856368 : green%oper(ifreq)%ks(ib,ib,ikpt+shift_green,isppol) = cone / green_tmp
1539 : else
1540 : green%oper(ifreq)%ks(ib,ib,ikpt+shift_green,isppol) = &
1541 54836656 : & green%oper(ifreq)%ks(ib,ib,ikpt+shift_green,isppol) + green_tmp
1542 : end if
1543 : end do ! ib
1544 : end do ! ikpt
1545 : end do ! isppol
1546 :
1547 : ! do ib1 = 1 , paw_dmft%mbandc
1548 : ! do ib = 1 , paw_dmft%mbandc
1549 : ! do ikpt = 1 , paw_dmft%nkpt
1550 : ! do is = 1 , paw_dmft%nsppol
1551 : ! if (green%procb(ifreq,ikpt)==myproc) then
1552 : ! green%oper(ifreq)%ks(is,ikpt,ib,ib1)= &
1553 : ! green_temp%ks(is,ikpt,ib,ib1)= &
1554 : !& ( omega_current &
1555 : !& + fermilevel &
1556 : !& - paw_dmft%eigen_dft(is,ikpt,ib)) * Id(ib,ib1) &
1557 : !& - self_minus_hdc_oper%ks(is,ikpt,ib,ib1)
1558 : !if(ikpt==2.and.ib==ib1) then
1559 : ! write(6,*)
1560 : ! "self",ib1,ib,ikpt,is,ifreq,self_minus_hdc_oper%ks(is,ikpt,ib,ib1)
1561 : !endif
1562 : !& -
1563 : !(self%oper(ifreq)%ks(is,ikpt,ib,ib1)-self%hdc%ks(is,ikpt,ib,ib1))
1564 : ! if(prtopt>5) then
1565 : ! if(ikpt==1.and.(ifreq==1.or.ifreq==3).and.ib==16.and.ib1==16)
1566 : ! then
1567 : ! write(std_out,*) 'omega_current
1568 : ! ',omega_current
1569 : ! write(std_out,*) 'fermilevel
1570 : ! ',fermilevel
1571 : ! write(std_out,*) ' paw_dmft%eigen_dft(is,ikpt,ib) ',
1572 : ! paw_dmft%eigen_dft(is,ikpt,ib),Id(ib,ib1)
1573 : ! write(std_out,*)
1574 : ! 'self_minus_hdc_oper%ks(is,ikpt,ib,ib1)',self_minus_hdc_oper%ks(is,ikpt,ib,ib1)
1575 : ! write(std_out,*) 'green
1576 : ! ',green%oper(ifreq)%ks(is,ikpt,ib,ib1)
1577 : ! endif
1578 : ! if(ib==1.and.ib1==3) then
1579 : ! write(std_out,*) "ff compute",ikpt,ifreq,is,ikpt,ib,ib1
1580 : ! write(std_out,*) "ff compute",ikpt,ifreq,
1581 : ! green_temp%ks(is,ikpt,ib,ib1)
1582 : ! write(std_out,*) "ff details",paw_dmft%eigen_dft(is,ikpt,ib)
1583 : ! write(std_out,*) "ff details2",fermilevel
1584 : ! write(std_out,*) "ff details3",Id(ib,ib1)
1585 : ! ! write(std_out,*) "ff
1586 : ! details4",self_minus_hdc_oper%ks(is,ikpt,ib,ib1)
1587 : ! endif
1588 : ! endif
1589 : ! enddo ! is
1590 : ! enddo ! ikpt
1591 : !enddo ! ib
1592 : !enddo ! ib1
1593 :
1594 : ! call print_oper(green%oper(ifreq),9,paw_dmft,3)
1595 : ! write(std_out,*) 'after print_oper'
1596 : ! if(ifreq==1.or.ifreq==3) then
1597 : ! write(std_out,*) 'after print_oper', ifreq
1598 : ! write(std_out,*) 'green1 ifreq %ks(1,1,16,16)',ifreq,green%oper(ifreq)%ks(1,1,16,16)
1599 : ! endif
1600 : ! write(std_out,*) 'before inverse_oper'
1601 994987 : if (optself /= 0) then
1602 949905 : call inverse_oper(green%oper(ifreq),1,procb=green%distrib%procb(:),iproc=me_kpt)
1603 : end if
1604 : !if(ifreq==1) then
1605 : ! write(std_out,*) "1188",green_temp%ks(1,1,8,8)
1606 : ! write(std_out,*) "1189",green_temp%ks(1,1,8,9)
1607 : ! write(std_out,*) "1198",green_temp%ks(1,1,9,8)
1608 : ! write(std_out,*) "1199",green_temp%ks(1,1,9,9)
1609 : !endif
1610 994987 : if (optlog == 1) then
1611 0 : freq2 = paw_dmft%omega_lo(ifreq)**2
1612 0 : trace_tmp = czero
1613 0 : do isppol=1,nsppol
1614 0 : do ikpt=1,mkmem
1615 0 : wtk = green%oper(ifreq)%wtk(ikpt+shift)
1616 0 : if (optself == 0) then
1617 0 : do ib=1,mbandc
1618 0 : trace_tmp = trace_tmp + two*temp*wtk*log(green%oper(ifreq)%ks(ib,ib,ikpt+shift_green,isppol)*omega_current)
1619 : end do ! ib
1620 : else
1621 : ! Use Tr(log(G(iw))) + Tr(log(G(-iw))) = Tr(log(G(iw)G(iw)^H)) so that we can use the much faster zheev instead of
1622 : ! zgeev. Even if G(iw) and G(iw)^H have no reason to commute, this equality still holds since we only care about the trace.
1623 : call abi_xgemm("n","c",mbandc,mbandc,mbandc,cone,green%oper(ifreq)%ks(:,:,ikpt+shift_green,isppol), &
1624 0 : & mbandc,green%oper(ifreq)%ks(:,:,ikpt+shift_green,isppol),mbandc,czero,mat_tmp(:,:),mbandc)
1625 0 : call zheev("n","u",mbandc,mat_tmp(:,:),mbandc,eig(:),work(:),lwork,rwork(:),info)
1626 : ! Do not use DOT_PRODUCT
1627 0 : trace_tmp = trace_tmp + sum(log(eig(:)*freq2))*wtk*temp
1628 : end if ! optself
1629 : end do ! ikpt
1630 : end do ! isppol
1631 0 : if (nsppol == 1 .and. nspinor == 1) trace_tmp = trace_tmp * two
1632 0 : green%trace_log = green%trace_log + dble(trace_tmp)
1633 : end if ! optlog=1
1634 :
1635 : !if(lintegrate) then
1636 : ! accumulate integration
1637 994987 : if (green%w_type /= "real") then
1638 :
1639 : !do isppol=1,nsppol
1640 : ! do ikpt=1,nkpt
1641 : ! if (green%distrib%procb(ikpt+(isppol-1)*nkpt) /= me_spkpt) cycle
1642 :
1643 : !do ib1=1,mbandc
1644 : ! do ib=1,mbandc
1645 : ! if (green%procb(ifreq,ikpt)==myproc) then
1646 : ! call
1647 : ! add_int_fct(ifreq,green%oper(ifreq)%ks(ib,ib1,ikpt,isppol),ib==ib1,
1648 : ! &
1649 : ! & omega_current,2,green%occup%ks(ib,ib1,ikpt,isppol), &
1650 : ! & paw_dmft%temp,paw_dmft%wgt_wlo(ifreq),paw_dmft%dmft_nwlo)
1651 :
1652 : !endif
1653 : ! end do ! ib
1654 : !end do ! ib1
1655 : ! end do ! ikpt
1656 : !end do ! isppol
1657 :
1658 982987 : if (diag == 1) then
1659 402510 : do ib=1,mbandc
1660 : green%occup%ks(ib,ib,1+shift:mkmem+shift,:) = green%occup%ks(ib,ib,1+shift:mkmem+shift,:) + &
1661 3835334 : & fac*green%oper(ifreq)%ks(ib,ib,1+shift_green:mkmem+shift_green,:)
1662 : end do ! ib
1663 : else
1664 : green%occup%ks(:,:,1+shift:mkmem+shift,:) = green%occup%ks(:,:,1+shift:mkmem+shift,:) + &
1665 492153235 : & fac*green%oper(ifreq)%ks(:,:,1+shift_green:mkmem+shift_green,:)
1666 : end if ! diag
1667 :
1668 : end if ! green%wtype
1669 :
1670 : ! write(std_out,*) 'after inverse_oper'
1671 : ! if(ikpt==1.and.is==1.and.ib==1.and.ib1==1) then
1672 : ! write(6,*) 'occup(is,ikpt,ib,ib1)',ifreq,green%occup%ks(1,1,1,1),green_temp%ks(1,1,1,1)
1673 : ! endif
1674 : ! write(std_out,*) 'green1afterinversion %ks(1,1,16,16)',ifreq,green%oper(ifreq)%ks(1,1,16,16)
1675 : ! endif
1676 : ! write(std_out,*) 'before flush'
1677 : ! call flush(std_out)
1678 : !endif
1679 : ! ================================
1680 : ! == Compute Local Green function
1681 : ! ================================
1682 : !write(message,'(2a)') ch10,' loc'
1683 : !call wrtout(std_out,message,'COLL')
1684 : !call flush(std_out)
1685 994987 : if (paw_dmft%lchipsiortho == 1 .or. optself == 1) then
1686 972414 : call downfold_oper(green%oper(ifreq),paw_dmft,procb=green%distrib%procb(:),iproc=me_kpt,option=option)
1687 : !if(ifreq==1) then
1688 : ! write(std_out,*) "4411",green_temp%matlu(1)%mat(4,4,1,1,1)
1689 : ! write(std_out,*) "4512",green_temp%matlu(1)%mat(4,5,1,1,2)
1690 : ! write(std_out,*) "5421",green_temp%matlu(1)%mat(5,4,1,2,1)
1691 : ! write(std_out,*) "5522",green_temp%matlu(1)%mat(5,5,1,2,2)
1692 : ! write(std_out,*) "(5512)",green_temp%matlu(1)%mat(5,5,1,1,2)
1693 : !endif
1694 : ! write(std_out,*) ifreq,nproc,'before if after loc_oper'
1695 : ! if(ifreq==1.or.ifreq==11) then
1696 : ! write(std_out,*) ifreq,nproc,'before sym'
1697 : ! call print_matlu(green%oper(ifreq)%matlu,green%oper(ifreq)%natom,2,-1,0)
1698 : ! ! ok
1699 : ! endif
1700 : ! call flush(std_out)
1701 : end if ! lchipsiortho=1
1702 : !call copy_matlu(green_temp%matlu,green%oper(ifreq)%matlu,natom)
1703 : ! if(ifreq==1.and.ifreq==11) then
1704 : ! write(std_out,*) ifreq,nproc,'after sym'
1705 : ! call print_matlu(green%oper(ifreq)%matlu,green%oper(ifreq)%natom,2,-1,0)
1706 : ! ! ok
1707 : ! endif
1708 999750 : if (green%oper(ifreq)%has_operks == 0) then
1709 982987 : ABI_FREE(green%oper(ifreq)%ks)
1710 : end if
1711 : ! call flush(std_out)
1712 : end do ! ifreq
1713 : ABI_NVTX_END_RANGE()
1714 : end if
1715 :
1716 :
1717 4884 : if (optlog == 1) then
1718 :
1719 0 : ABI_FREE(eig)
1720 0 : ABI_FREE(work)
1721 0 : ABI_FREE(rwork)
1722 0 : ABI_FREE(mat_tmp)
1723 :
1724 0 : call xmpi_sum(green%trace_log,spacecomm,ierr)
1725 :
1726 0 : correction = - log(two) * mbandc * nsppol * temp
1727 :
1728 0 : band_index = 0
1729 0 : beta = one / paw_dmft%temp
1730 0 : do isppol=1,nsppol
1731 0 : do ikpt=1,nkpt
1732 0 : nband_k = paw_dmft%nband(ikpt+(isppol-1)*nkpt)
1733 0 : do ib=1,nband_k
1734 0 : if (paw_dmft%band_in(ib)) cycle
1735 0 : eigen = paw_dmft%eigen(ib+band_index)
1736 0 : if (eigen-fermilevel >= zero) then
1737 0 : correction = correction - temp*paw_dmft%wtk(ikpt)*log(one+exp(-beta*(eigen-fermilevel)))
1738 : else
1739 0 : correction = correction - temp*paw_dmft%wtk(ikpt)*(log(one+exp(beta*(eigen-fermilevel)))-beta*(eigen-fermilevel))
1740 : end if ! eigen-fermilevel>=0
1741 : end do ! ib
1742 0 : band_index = band_index + nband_k
1743 : end do ! ikpt
1744 : end do ! isppol
1745 0 : if (nsppol == 1 .and. nspinor == 1) correction = correction * two
1746 :
1747 : ! Do not use DOT_PRODUCT
1748 : green%trace_log = green%trace_log + correction + fermilevel*paw_dmft%nelectval + &
1749 0 : & dble(sum(green%trace_moments_log_ks(1:nmoments-1)*omega_fac(1:nmoments-1)))
1750 :
1751 : end if ! optlog=1
1752 :
1753 : ! =============================================
1754 : ! Build total green function (sum over procs).
1755 : ! =============================================
1756 : !call xmpi_barrier(spacecomm)
1757 : ! call xmpi_sum(green%oper(ifreq)%ks,spacecomm,ierr)
1758 : ! print *, "myproc, proct, ifreq ------------------- ", myproc, ifreq
1759 : ! do ikpt=1,paw_dmft%nkpt
1760 : ! call xmpi_bcast(green%oper(ifreq)%ks(:,ikpt,:,:),procb(ifreq,ikpt),spacecomm,ierr)
1761 : ! enddo
1762 : !! or
1763 4884 : if (optnonxsum == 0 .and. green%oper(1)%has_operks == 1) then
1764 4 : call gather_oper(green%oper(:),green%distrib,paw_dmft,opt_ksloc=1,opt_diag=diag)
1765 0 : else if (optnonxsum == 0 .and. green%oper(1)%has_operks == 0) then
1766 0 : message = 'optnonxsum=0 and green%oper(1)%has_operks=0: not compatible'
1767 0 : ABI_BUG(message)
1768 : end if ! optnonxsum
1769 : ! endif
1770 : ! enddo ! ifreq
1771 : ! print *,"myproc", myproc
1772 4884 : if (optnonxsum2 == 0) then
1773 681 : if (paw_dmft%lchipsiortho == 1 .or. optself == 1) then
1774 587 : call gather_oper(green%oper(:),green%distrib,paw_dmft,opt_ksloc=2,opt_commkpt=1)
1775 :
1776 587 : if(gpu_option==ABI_GPU_DISABLED) then
1777 149153 : do ifreq=1,green%nw
1778 148566 : if (green%distrib%procf(ifreq) /= myproc) cycle
1779 149153 : call sym_matlu(green%oper(ifreq)%matlu(:),paw_dmft)
1780 : end do
1781 : else if(gpu_option==ABI_GPU_OPENMP) then
1782 : #ifdef HAVE_OPENMP_OFFLOAD
1783 : !FIXME: Remove those CPU-GPU transfers once gather_oper has been ported to GPU
1784 : ! 1) Init green_oper_ndat
1785 : ndat = green%distrib%nw_mem_kptparal(green%distrib%me_freq+1)
1786 : do ifreq=1,green%nw
1787 : if (green%distrib%procf(ifreq) /= myproc) cycle
1788 : call init_oper_ndat(paw_dmft,green_oper_ndat,ndat,nkpt=green%oper(ifreq)%nkpt,opt_ksloc=2,gpu_option=gpu_option)
1789 : if (green%oper(ifreq)%has_operks == 0) then
1790 : green_oper_ndat%paral = 1
1791 : green_oper_ndat%shiftk = shift
1792 : end if
1793 : exit
1794 : end do
1795 : ! 2) Copy green%oper(:)%matlu into green_oper_ndat (CPU->GPU transfer)
1796 : call copy_oper_to_ndat(green%oper,green_oper_ndat,ndat,green%nw,green%distrib%proct,green%distrib%me_freq,.false.)
1797 :
1798 : ! 3) Perform sym_matlu on green_oper_ndat (GPU enabled)
1799 : call sym_matlu(green_oper_ndat%matlu(:),paw_dmft)
1800 :
1801 : ! 4) Copy back green%oper(:)%matlu from green_oper_ndat (GPU->CPU transfer)
1802 : call copy_oper_from_ndat(green_oper_ndat,green%oper,ndat,green%nw,green%distrib%proct,&
1803 : & green%distrib%me_freq,.false.)
1804 : ! 5) Destroy green_oper_ndat
1805 : call destroy_oper(green_oper_ndat)
1806 : #endif
1807 : end if
1808 :
1809 587 : call gather_oper(green%oper(:),green%distrib,paw_dmft,opt_ksloc=2)
1810 : end if
1811 681 : green%has_greenmatlu_xsum = 1
1812 4203 : else if (optnonxsum2 == 1) then
1813 4203 : green%has_greenmatlu_xsum = 0
1814 : end if ! optnonxsum2
1815 : ! if(ifreq==1.or.ifreq==11) then
1816 : ! write(std_out,*) ifreq,nproc,'after xsum'
1817 : ! call print_matlu(green%oper(ifreq)%matlu,green%oper(ifreq)%natom,2,-1,0)
1818 : ! ! ok
1819 : ! endif
1820 : !end do ! ifreq
1821 : !call xmpi_barrier(spacecomm)
1822 : ! write(std_out,*) 'afterxsum sym %matlu(1)%mat(2,5,1,1,1) 1',green%oper(1)%matlu(1)%mat(2,5,1,1,1)
1823 :
1824 4884 : if (green%w_type /= "real") then
1825 :
1826 4880 : if (green%distrib%me_freq == 0) then
1827 :
1828 68179 : do ib=1,mbandc
1829 355207 : green%occup%ks(ib,ib,1+shift:mkmem+shift,:) = green%occup%ks(ib,ib,1+shift:mkmem+shift,:) + omega_fac(1)
1830 : end do ! ib
1831 :
1832 4759 : do i=2,nmoments
1833 4759 : if (diag == 1) then
1834 0 : do ib=1,mbandc
1835 : green%occup%ks(ib,ib,1+shift:mkmem+shift,:) = green%occup%ks(ib,ib,1+shift:mkmem+shift,:) + &
1836 0 : & omega_fac(i)*green%moments(i)%ks(ib,ib,:,:)
1837 : end do ! ib
1838 : else
1839 : green%occup%ks(:,:,1+shift:mkmem+shift,:) = green%occup%ks(:,:,1+shift:mkmem+shift,:) + &
1840 0 : & omega_fac(i)*green%moments(i)%ks(:,:,:,:)
1841 : end if ! diag
1842 : end do ! i
1843 :
1844 : end if ! me_freq = 0
1845 :
1846 4880 : call gather_oper_ks(green%occup,green%distrib,paw_dmft,opt_diag=diag)
1847 :
1848 : end if ! w_type/="real"
1849 :
1850 4884 : ABI_SFREE(omega_fac)
1851 :
1852 4884 : if (prtopt /= 0 .and. prtopt > -100) then
1853 470 : write(message,'(2a)') ch10," === Green's function is computed"
1854 470 : call wrtout(std_out,message,'COLL')
1855 : end if
1856 :
1857 4884 : if (prtopt /= 0 .and. prtopt > -100 .and. (paw_dmft%lchipsiortho == 1 .or. optself == 1)) then
1858 376 : write(message,'(2a)') ch10,&
1859 752 : & " === Local Green's function has been computed and projected on local orbitals"
1860 376 : call wrtout(std_out,message,'COLL')
1861 : end if
1862 : ! useless test
1863 4884 : if (abs(prtopt) >= 4 .and. prtopt > -100) then
1864 0 : write(message,'(2a)') ch10," == Green's function is now printed for first frequency"
1865 0 : call wrtout(std_out,message,'COLL')
1866 0 : call print_oper(green%oper(1),9,paw_dmft,3)
1867 0 : write(message,'(2a)') ch10," == Green's function is now printed for second frequency"
1868 0 : call wrtout(std_out,message,'COLL')
1869 0 : call print_oper(green%oper(2),9,paw_dmft,3)
1870 0 : if (paw_dmft%dmft_nwlo >= 11) then
1871 0 : write(message,'(2a)') ch10," == Green's function is now printed for 11th frequency"
1872 0 : call wrtout(std_out,message,'COLL')
1873 0 : call print_oper(green%oper(11),9,paw_dmft,3)
1874 : end if
1875 : end if
1876 : ! call flush(std_out)
1877 :
1878 : ABI_NVTX_END_RANGE()
1879 4884 : call timab(624,2,tsec(:))
1880 :
1881 4884 : end subroutine compute_green
1882 : !!***
1883 :
1884 : !!****f* m_green/integrate_green
1885 : !! NAME
1886 : !! integrate_green
1887 : !!
1888 : !! FUNCTION
1889 : !! Integrate green function
1890 : !!
1891 : !! INPUTS
1892 : !! green <type(green_type)>=green function (green%oper(:))
1893 : !! paw_dmft <type(m_paw_dmft)>= paw+dmft data
1894 : !! prtopt : flag for printing
1895 : !! opt_ksloc= 1: only integrate on the KS basis
1896 : !! 2: do the integration on the local basis
1897 : !! (This can work only if chipsi are renormalized!!!)
1898 : !! 3: do both calculations and test the consistency of it
1899 : !! -1: do the integration on the KS basis, but only
1900 : !! compute diagonal part of the band-band density matrix
1901 : !! in order to compute the total charge for fermi_green
1902 : !! opt_after_solver= to be activated after a call to impurity_solve
1903 : !! opt_diff= check the change in the local charge compared to the previous iteration
1904 : !! opt_fill_occnd= 0 (default) : do not fill paw_dmft%occnd with the Green's function occupations
1905 : !! = 1 : fill paw_dmft%occnd with the Green's function occupations
1906 : !! opt_self= 0 : assume the green%ks is diagonal
1907 : !! = 1 (default) : assume there is a self-energy and green%ks is non-diagonal
1908 : !!
1909 : !! OUTPUT
1910 : !! green%occup = occupations
1911 : !!
1912 : !! SOURCE
1913 :
1914 4984 : subroutine integrate_green(green,paw_dmft,prtopt,opt_ksloc,opt_after_solver,opt_diff,opt_fill_occnd,opt_self)
1915 :
1916 : !Arguments ------------------------------------
1917 : type(green_type), intent(inout) :: green
1918 : !type(MPI_type), intent(in) :: mpi_enreg
1919 : type(paw_dmft_type), intent(inout) :: paw_dmft
1920 : integer, intent(in) :: prtopt
1921 : integer, optional, intent(in) :: opt_after_solver,opt_diff,opt_fill_occnd,opt_ksloc,opt_self
1922 : !Local variables-------------------------------
1923 : integer :: band_index,i,iatom,ib,ib1,icomp_chloc,ifreq,ikpt,isppol
1924 : integer :: lpawu,mbandc,myproc,natom,nband_k,nkpt,nmoments,nspinor
1925 : integer :: nsppol,optaftsolv,optdiff,optfilloccnd,option,optksloc,optself,optiondiff
1926 : real(dp) :: correction,diff_chloc,fac,temp
1927 : character(len=12) :: tag
1928 : character(len=500) :: message
1929 : real(dp) :: tsec(2)
1930 4984 : complex(dp), allocatable :: omega_fac(:),shift(:)
1931 4984 : type(matlu_type), allocatable :: matlu_temp(:)
1932 : ! real(dp), allocatable :: charge_loc_old(:,:)
1933 : ! type(oper_type) :: oper_c
1934 : ! *********************************************************************
1935 :
1936 : DBG_ENTER("COLL")
1937 4984 : call timab(625,1,tsec(:))
1938 : ABI_NVTX_START_RANGE(NVTX_DMFT_INTEGRATE_GREEN)
1939 :
1940 4984 : if (prtopt > 0) then
1941 582 : write(message,'(2a,i3,13x,a)') ch10," === Integrate Green's function"
1942 582 : call wrtout(std_out,message,'COLL')
1943 : end if
1944 4984 : if (green%w_type == "real") then
1945 0 : message = 'integrate_green not implemented for real frequency'
1946 0 : ABI_BUG(message)
1947 : end if
1948 :
1949 4984 : optfilloccnd = 0
1950 4984 : if (present(opt_fill_occnd)) optfilloccnd = opt_fill_occnd
1951 :
1952 4984 : optself = 1
1953 4984 : if (present(opt_self)) optself = opt_self
1954 279 : option = merge(3,1,optself==0) ! option for downfold_oper
1955 :
1956 : ! Initialize spaceComm, myproc, and master
1957 : !spacecomm=paw_dmft%spacecomm
1958 4984 : myproc = paw_dmft%myproc
1959 : !nproc=paw_dmft%nproc
1960 :
1961 : ! Initialise integers
1962 : !mband = paw_dmft%mband
1963 4984 : mbandc = paw_dmft%mbandc
1964 4984 : natom = paw_dmft%natom
1965 4984 : nkpt = paw_dmft%nkpt
1966 4984 : nspinor = paw_dmft%nspinor
1967 4984 : nsppol = paw_dmft%nsppol
1968 4984 : temp = paw_dmft%temp
1969 :
1970 4984 : nmoments = merge(green%nmoments,1,green%has_moments==1)
1971 :
1972 : ! Initialize green%oper before calculation (important for xmpi_sum)
1973 : ! allocate(charge_loc_old(paw_dmft%natom,paw_dmft%nsppol+1))
1974 : ! if(.not.present(opt_diff)) then ! if integrate_green is called in m_dmft after calculation of self
1975 : ! charge_loc_old=green%charge_matlu
1976 : ! endif
1977 4984 : icomp_chloc = 0
1978 4984 : optdiff = 0
1979 4984 : if (present(opt_diff)) optdiff = opt_diff
1980 :
1981 : ! Choose what to compute
1982 4984 : optksloc = 3
1983 4984 : if (present(opt_ksloc)) optksloc = opt_ksloc
1984 4984 : optaftsolv = 0
1985 4984 : if (present(opt_after_solver)) optaftsolv = opt_after_solver
1986 :
1987 103 : if (optaftsolv == 1 .and. abs(optksloc) /= 2) then
1988 0 : message = "integration of ks green function should not be done after call to solver : it has not been computed"
1989 0 : ABI_BUG(message)
1990 : end if
1991 :
1992 4984 : if (abs(optksloc) >= 2 .and. green%has_greenmatlu_xsum == 0) then
1993 0 : write(message,'(4a)') ch10,&
1994 0 : & "BUG: integrate_green is asked to integrate local green function",ch10,&
1995 0 : & " and local green function was non broadcasted in compute_green"
1996 0 : ABI_BUG(message)
1997 : end if
1998 :
1999 : ! Allocations
2000 33684 : ABI_MALLOC(matlu_temp,(natom))
2001 4984 : call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu(:),matlu_temp(:))
2002 :
2003 : ! =================================================
2004 : ! == Integrate Local Green function ===============
2005 4984 : if (abs(optksloc)/2 == 1) then ! optksloc=2 or 3
2006 : ! =================================================
2007 :
2008 : ! == Calculation of \int{G_{LL'}{\sigma\sigma',s}(R)(i\omega_n)}
2009 781 : if (paw_dmft%lchipsiortho == 1) then
2010 : ! - Calculation of frequency sum over positive frequency
2011 : !if (nspinor==1) option=1
2012 : !if (nspinor==2) option=2
2013 :
2014 : !do atom=1, natom
2015 : ! ndim=2*green%oper(1)%matlu(iatom)%lpawu+1
2016 : ! if(green%oper(1)%matlu(iatom)%lpawu.ne.-1) then
2017 : ! do ispinor1 = 1, nspinor
2018 : ! do ispinor = 1, nspinor
2019 : ! do is = 1 , nsppol
2020 : ! do im=1,ndim
2021 : ! do im1=1,ndim
2022 : ! do ifreq=1,green%nw
2023 : ! ff(ifreq)= &
2024 : !& green%oper(ifreq)%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1)
2025 : ! write(std_out,*) green%omega(ifreq),ff(ifreq)," integrate
2026 : ! green fw_lo"
2027 : !if(present(iii).and.im==1.and.im1==1) write(std_out_default,*)
2028 : !ch10,ifreq,ff(ifreq),"#ff"
2029 : ! enddo
2030 : ! call int_fct(ff,(im==im1).and.(ispinor==ispinor1),&
2031 : !& option,paw_dmft,integral)
2032 : ! call int_fct(ff,(im==im1).and.(ispinor==ispinor1),&
2033 : !& 2,paw_dmft,integral) ! test_1
2034 : ! green%occup%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1)=integral
2035 : !if(present(iii).and.im==1.and.im1==1) write(std_out_default,*)
2036 : !ch10,'integral',im,im1,ifreq,integral
2037 : ! if(im==2.and.im1==5.and.is==1.and.iatom==1) then
2038 : ! write(std_out,*) " occup
2039 : ! %matlu(1)%mat(2,5,1,1,1)",green%occup%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1)
2040 : ! endif
2041 : ! enddo
2042 : ! enddo
2043 : ! enddo ! ispinor1
2044 : ! enddo ! ispinor
2045 : ! enddo ! is
2046 : ! matlu_temp(iatom)%mat=green%occup%matlu(iatom)%mat
2047 : ! endif ! lpawu=/-1
2048 : ! enddo ! iatom
2049 :
2050 682 : call zero_matlu(green%occup%matlu(:),natom)
2051 :
2052 2046 : ABI_MALLOC(omega_fac,(nmoments))
2053 :
2054 1364 : do i=1,nmoments
2055 682 : omega_fac(i) = czero
2056 160113 : do ifreq=green%nw,1,-1 ! NEVER change the summation order and DON'T use the intrinsic SUM
2057 160113 : omega_fac(i) = omega_fac(i) + paw_dmft%wgt_wlo(ifreq) / (paw_dmft%omega_lo(ifreq))**i
2058 : end do
2059 682 : omega_fac(i) = - two * temp * omega_fac(i) / (j_dpc)**i
2060 682 : if (i == 1) omega_fac(i) = omega_fac(i) + half
2061 682 : if (i == 2) omega_fac(i) = omega_fac(i) - cone/(four*temp)
2062 1364 : if (i == 4) omega_fac(i) = omega_fac(i) + cone/(dble(48)*(temp**3))
2063 : end do ! i
2064 :
2065 160113 : do ifreq=1,green%nw
2066 :
2067 159431 : if (green%distrib%procf(ifreq) /= myproc) cycle
2068 54800 : fac = two * temp * paw_dmft%wgt_wlo(ifreq)
2069 161146 : do iatom=1,natom
2070 105664 : lpawu = paw_dmft%lpawu(iatom)
2071 105664 : if (lpawu == -1) cycle
2072 : green%occup%matlu(iatom)%mat(:,:,:) = green%occup%matlu(iatom)%mat(:,:,:) + &
2073 3437655 : & fac*green%oper(ifreq)%matlu(iatom)%mat(:,:,:)
2074 : end do ! iatom
2075 :
2076 : end do ! ifreq
2077 :
2078 682 : call xmpi_matlu(green%occup%matlu(:),natom,paw_dmft%spacecomm)
2079 :
2080 682 : if (green%has_moments == 1) then
2081 0 : do i=1,nmoments
2082 0 : do iatom=1,natom
2083 0 : lpawu = paw_dmft%lpawu(iatom)
2084 0 : if (lpawu == -1) cycle
2085 : green%occup%matlu(iatom)%mat(:,:,:) = green%occup%matlu(iatom)%mat(:,:,:) + &
2086 0 : omega_fac(i)*green%moments(i)%matlu(iatom)%mat(:,:,:)
2087 : end do ! iatom
2088 : end do ! i
2089 : else
2090 2046 : ABI_MALLOC(shift,(natom))
2091 3060 : shift(:) = omega_fac(1)
2092 682 : call shift_matlu(green%occup%matlu(:),natom,shift(:),signe=-1)
2093 682 : ABI_FREE(shift)
2094 : end if
2095 :
2096 682 : ABI_FREE(omega_fac)
2097 :
2098 : ! Print density matrix if prtopt high
2099 682 : if (abs(prtopt) > 2) then
2100 105 : write(message,'(2a,i10,a)') ch10," = green%occup%matlu from int(gloc(w))"
2101 105 : call wrtout(std_out,message,'COLL')
2102 105 : call print_matlu(green%occup%matlu(:),natom,prtopt=3,opt_diag=-1)
2103 : end if ! abs(prtopt)>2
2104 :
2105 : ! - Symmetrize: continue sum over k-point: Full BZ
2106 682 : call sym_matlu(green%occup%matlu(:),paw_dmft)
2107 682 : if (abs(prtopt) > 2) then
2108 105 : write(message,'(2a,i10,a)') ch10, &
2109 210 : & " = green%occup%matlu from int(gloc(w)) with symmetrization"
2110 105 : call wrtout(std_out,message,'COLL')
2111 105 : call print_matlu(green%occup%matlu(:),natom,prtopt=3,opt_diag=-1)
2112 : end if ! abs(prtopt)>2
2113 :
2114 : ! - Post-treatment for summation over negative and positive frequencies:
2115 : ! necessary in the case of nspinor==2 AND nspinor==1, but valid anywhere
2116 : ! N(ll'sigmasigma')= (N(ll'sigmasigma')+ N*(l'lsigma'sigma))/2
2117 : ! because [G_{LL'}^{sigma,sigma'}(iomega_n)]*= G_{L'L}^{sigma',sigma}(-iomega_n)
2118 3060 : do iatom=1,natom
2119 2378 : lpawu = paw_dmft%lpawu(iatom)
2120 2378 : if (lpawu == -1) cycle
2121 2018 : do isppol=1,nsppol
2122 : green%occup%matlu(iatom)%mat(:,:,isppol) = (green%occup%matlu(iatom)%mat(:,:,isppol) + &
2123 71522 : & conjg(transpose(green%occup%matlu(iatom)%mat(:,:,isppol)))) * half
2124 : end do ! isppol
2125 39028 : matlu_temp(iatom)%mat(:,:,:) = green%occup%matlu(iatom)%mat(:,:,:)
2126 : end do ! iatom
2127 682 : if (abs(prtopt) > 2) then
2128 105 : write(message,'(2a,i10,a)') ch10, &
2129 210 : & " = green%occup%matlu from int(gloc(w)) symmetrized with post-treatment"
2130 105 : call wrtout(std_out,message,'COLL')
2131 105 : call print_matlu(green%occup%matlu(:),natom,prtopt=3,opt_diag=-1)
2132 : end if ! abs(prtopt)>2
2133 :
2134 682 : if (optaftsolv == 0) then
2135 579 : call trace_oper(green%occup,green%charge_ks,green%charge_matlu(:,:),2)
2136 579 : green%has_charge_matlu = 2
2137 579 : green%has_charge_matlu_solver = 0
2138 579 : icomp_chloc = 1
2139 103 : else if (optaftsolv == 1) then ! .and. paw_dmft%dmft_solv /= 4) then
2140 : ! else if(optaftsolv==1) then
2141 : !if(paw_dmft%dmft_solv==4) then
2142 : ! write(message,'(a,a,a)') ch10,&
2143 : !& " = WARNING: Double Counting will be computed with solver charge.." ,&
2144 : !& "might be problematic with Hirsch Fye QMC"
2145 : ! call wrtout(std_out,message,'COLL')
2146 : ! endif
2147 : ! This is done only when called from impurity_solver with solver
2148 : ! green function. (And if QMC is NOT used).
2149 103 : if (paw_dmft%dmft_solv >= 5 .and. green%has_charge_matlu_solver /= 2) then
2150 0 : write(message,'(2a,i3)') ch10,&
2151 0 : & " = BUG : has_charge_matlu_solver should be 2 and is",green%has_charge_matlu_solver
2152 0 : ABI_BUG(message)
2153 : end if
2154 103 : if (paw_dmft%dmft_solv <= 4) then
2155 69 : call trace_oper(green%occup,green%charge_ks,green%charge_matlu_solver(:,:),2)
2156 69 : green%has_charge_matlu_solver = 2
2157 : end if
2158 : end if ! optaftsolv
2159 : else
2160 99 : write(message,'(a,4x,3a,4x,a)') ch10,&
2161 99 : & " Local basis is not (yet) orthonormal:",&
2162 99 : & " local Green's function is thus not integrated",ch10,&
2163 198 : & " Local occupations are computed by downfolding the Kohn-Sham occupations instead"
2164 99 : call wrtout(std_out,message,'COLL')
2165 : end if ! lchipsiortho=1
2166 :
2167 : end if ! optksloc
2168 : ! =================================================
2169 :
2170 : ! =================================================
2171 : ! == Integrate Kohn Sham Green function ===========
2172 4984 : if (mod(abs(optksloc),2) == 1) then ! optksloc=1 or 3 or -1
2173 : ! green%occup%ks=czero ! important for xmpi_sum
2174 : !! =================================================
2175 : !! == Calculation of \int{G_{\nu\nu'}{k,s}(i\omega_n)}
2176 : ! ff=czero
2177 : ! do is = 1 , nsppol
2178 : ! do ikpt = 1, nkpt
2179 : !!020212 if(mod(ikpt-1,nproc)==myproc) then
2180 : !!! print'(A,3I4)', "P ",myproc,is,ikpt
2181 : ! do ib = 1, mbandc
2182 : ! do ib1 = 1, mbandc
2183 : ! if(optksloc==-1.and.(ib/=ib1)) cycle
2184 : ! ff(:)=czero
2185 : ! do ifreq=1,green%nw
2186 : ! !! the following line is a quick and dirty tricks and should be removed when the data will
2187 : ! !! be correctly distributed.
2188 : !!! if((optnonxsum==1).and.(green%procb(ifreq,ikpt)==myproc)).or.(optnonxsum==0)) then
2189 : ! if(green%procb(ifreq,ikpt)==myproc) then
2190 : ! ff(ifreq)=green%oper(ifreq)%ks(is,ikpt,ib,ib1)
2191 : !!! print'(A,5I4,3(2E15.5,3x))', "P1",myproc,is,ikpt,ib,ib1,ff(:)
2192 : ! endif
2193 : ! !endif
2194 : ! enddo
2195 : !! call int_fct(ff,ib==ib1,nspinor,paw_dmft,integral) ! here, option==1 even if nspinor==2
2196 : !! green%occup%ks(is,ikpt,ib,ib1)=integral
2197 : !!! print'(A,5I4,3(2E15.5,3x))', "PP",myproc,is,ikpt,ib,ib1,ff(:)
2198 : ! call int_fct(ff,ib==ib1,2,paw_dmft,integral,green%procb(:,ikpt),myproc) ! here, option==1 even if nspinor==2
2199 : !!020212 call int_fct(ff,ib==ib1,2,paw_dmft,integral) ! here, option==1 even if nspinor==2
2200 : ! green%occup%ks(is,ikpt,ib,ib1)=integral
2201 : !!! print'(A,5I4,2E15.8)', "P2",myproc,is,ikpt,ib,ib1,integral
2202 : ! ! write(std_out,*) "integral",ikpt,green%occup%ks(is,ikpt,ib,ib1)
2203 : ! ! endif
2204 : !! write(std_out,'(a,4i6,e14.5,e14.5)') "ks",is,ikpt,ib,ib1,integral
2205 : ! enddo ! ib1
2206 : ! enddo ! ib
2207 : !!020212 endif
2208 : ! enddo ! ikpt
2209 : ! enddo ! isppol
2210 :
2211 :
2212 : ! call xmpi_barrier(spacecomm)
2213 : ! call xmpi_sum(green%occup%ks,spacecomm,ierr)
2214 :
2215 :
2216 : !! do is = 1 , nsppol
2217 : !! do ikpt = 1, nkpt
2218 : !! do ib = 1, mbandc
2219 : !! do ib1 = 1, mbandc
2220 : !! write(6,'(A,5I4,2E15.8)') "AAAFTERXSUM",myproc,is,ikpt,ib,ib1,green%occup%ks(is,ikpt,ib,ib1)
2221 : !! print '(A,5I4,2E15.8)', "AFTERXSUM P",myproc,is,ikpt,ib,ib1,green%occup%ks(is,ikpt,ib,ib1)
2222 : !! enddo ! ib1
2223 : !! enddo ! ib
2224 : !! enddo ! ikpt
2225 : !! enddo ! isppol
2226 : !write(std_out,*) "occup%ks ik1",green%occup%ks(1,1,1,3)
2227 : !write(std_out,*) "occup%ks ik2",green%occup%ks(1,2,1,3)
2228 : ! - Post-treatment for summation over negative and positive frequencies:
2229 : ! necessary in the case of nspinor==2, but valid everywhere
2230 : ! N(k,n_1,n_2)= (N(k,n_1,n_2)+ N*(k,n_2,n_1))/2
2231 : ! because [G_{k}^{n_1,n_2}(iomega_n)]*= G_{k}^{n_2,n_1}(-iomega_n)
2232 11564 : do isppol=1,nsppol
2233 60270 : do ikpt=1,nkpt
2234 : green%occup%ks(:,:,ikpt,isppol) = (green%occup%ks(:,:,ikpt,isppol)+ &
2235 14174846 : & conjg(transpose(green%occup%ks(:,:,ikpt,isppol)))) * half
2236 : end do ! ikpt
2237 : end do ! isppol
2238 : !write(std_out,*) "occup%ks ik1 BB",green%occup%ks(1,1,1,3)
2239 : !write(std_out,*) "occup%ks ik2 AA",green%occup%ks(1,2,1,3)
2240 4880 : if (optfilloccnd == 1) then
2241 199 : band_index = 0
2242 199 : fac = merge(two,one,nsppol==1.and.nspinor==1)
2243 486 : do isppol=1,nsppol
2244 2708 : do ikpt=1,nkpt
2245 23638 : do ib1=1,mbandc
2246 319390 : do ib=1,mbandc
2247 : paw_dmft%occnd(1,paw_dmft%include_bands(ib),paw_dmft%include_bands(ib1),&
2248 295752 : & ikpt,isppol) = fac * dble(green%occup%ks(ib,ib1,ikpt,isppol))
2249 : paw_dmft%occnd(2,paw_dmft%include_bands(ib),paw_dmft%include_bands(ib1),&
2250 317168 : & ikpt,isppol) = fac * aimag(green%occup%ks(ib,ib1,ikpt,isppol))
2251 : end do ! ib
2252 : end do ! ib1
2253 2509 : if (paw_dmft%dmft_solv == 6 .or. paw_dmft%dmft_solv == 7) then
2254 0 : nband_k = paw_dmft%nband(ikpt+(isppol-1)*nkpt)
2255 0 : do ib=1,nband_k
2256 0 : if (paw_dmft%band_in(ib)) cycle
2257 0 : paw_dmft%occnd(:,:,ib,ikpt,isppol) = zero
2258 : paw_dmft%occnd(1,ib,ib,ikpt,isppol) = fac * &
2259 0 : & occup_fd(paw_dmft%eigen(ib+band_index),paw_dmft%fermie,paw_dmft%temp)
2260 : end do ! ib
2261 0 : band_index = band_index + nband_k
2262 : end if ! use_all_bands
2263 : end do ! ikpt
2264 : end do ! isppol
2265 : end if ! optfilloccnd
2266 :
2267 4880 : if (optksloc > 0) then
2268 : ! - Compute local occupations
2269 : call downfold_oper(green%occup,paw_dmft,procb=green%distrib%procb(:), &
2270 677 : & iproc=green%distrib%me_kpt,option=option)
2271 677 : call xmpi_matlu(green%occup%matlu(:),natom,green%distrib%comm_kpt)
2272 677 : if (abs(prtopt) > 2) then
2273 106 : write(message,'(2a,i10,a)') ch10,&
2274 212 : & " = green%occup%matlu from projection of int(gks(w)) without symmetrization"
2275 106 : call wrtout(std_out,message,'COLL')
2276 106 : call print_matlu(green%occup%matlu(:),natom,prtopt=3,opt_diag=-1)
2277 : end if ! abs(prtopt)>2
2278 :
2279 : ! - Symmetrize: continue sum over k-point: Full BZ
2280 677 : call sym_matlu(green%occup%matlu(:),paw_dmft)
2281 677 : if (abs(prtopt) >= 2) then
2282 : ! write(message,'(a,a,i10,a)') ch10,&
2283 : !& " = green%occup%matlu from projection of int(gks(w)) with symetrization"
2284 478 : write(message,'(2a,i10,a)') ch10," == Occupation matrix from downfolded Kohn-Sham occupations"
2285 478 : call wrtout(std_out,message,'COLL')
2286 478 : call print_matlu(green%occup%matlu(:),natom,prtopt=3,opt_diag=0)
2287 : end if ! abs(prtopt)>=2
2288 :
2289 : ! - If the trace of occup matrix in the LOCAL basis was not done
2290 : ! before because lchipsiortho/=1 , do it now
2291 677 : if (paw_dmft%lchipsiortho /= 1 .and. abs(prtopt) > 0) then
2292 97 : call trace_oper(green%occup,green%charge_ks,green%charge_matlu(:,:),2)
2293 97 : green%has_charge_matlu = 2
2294 97 : icomp_chloc = 1
2295 : end if ! lchipsiortho/=1 and abs(prtopt)>0
2296 : end if ! optksloc>0: only diagonal elements of G\nunu' are computed
2297 :
2298 : ! - Compute trace over ks density matrix
2299 4880 : call trace_oper(green%occup,green%charge_ks,green%charge_matlu(:,:),1)
2300 4880 : if (paw_dmft%dmft_solv == 6 .or. paw_dmft%dmft_solv == 7) then
2301 :
2302 : band_index = 0
2303 : correction = zero
2304 0 : do isppol=1,nsppol
2305 0 : do ikpt=1,nkpt
2306 0 : nband_k = paw_dmft%nband(ikpt+(isppol-1)*nkpt)
2307 0 : do ib=1,nband_k
2308 0 : if (paw_dmft%band_in(ib)) cycle
2309 : correction = correction + &
2310 0 : & occup_fd(paw_dmft%eigen(ib+band_index),paw_dmft%fermie,paw_dmft%temp)*green%occup%wtk(ikpt)
2311 : end do ! ib
2312 0 : band_index = band_index + nband_k
2313 : end do ! ikpt
2314 : end do ! isppol
2315 0 : if (nsppol == 1 .and. nspinor == 1) correction = correction * two
2316 0 : green%charge_ks = green%charge_ks + correction
2317 : end if ! use_all_bands
2318 4880 : if (abs(prtopt) > 0) then
2319 478 : write(tag,'(f12.6)') green%charge_ks
2320 478 : write(message,'(3a)') ch10,&
2321 956 : & " == Total number of electrons from integration of Kohn-Sham Green's function is : ",adjustl(tag)
2322 478 : call wrtout(std_out,message,'COLL')
2323 478 : write(tag,'(f12.6)') paw_dmft%nelectval
2324 478 : write(message,'(8x,3a)') " (should be ",trim(adjustl(tag)),")"
2325 478 : call wrtout(std_out,message,'COLL')
2326 : end if ! abs(prtopt)>0
2327 : end if ! optksloc
2328 : ! =================================================
2329 :
2330 : ! =================================================
2331 : ! Tests and compute precision on local charge
2332 : ! =================================================
2333 : ! - Check that if, renormalized psichi are used, occupations matrices
2334 : ! obtained directly from local green function or, through kohn sham
2335 : ! occupations are the same.
2336 4984 : if ((abs(optksloc) == 3) .and. (paw_dmft%lchipsiortho == 1)) then ! optksloc= 3
2337 :
2338 578 : if(nspinor == 2) then
2339 8 : optiondiff = 2
2340 : else
2341 570 : optiondiff = 1
2342 : endif
2343 : call diff_matlu("Local projection of Kohn-Sham occupations ",&
2344 : & "Integration of local Green's function ",&
2345 578 : & green%occup%matlu(:),matlu_temp(:),natom,optiondiff,tol4)
2346 578 : write(message,'(2a)') ch10,&
2347 1156 : & " ***** => Calculations of Green's function in Kohn-Sham and local spaces are coherent ****"
2348 578 : call wrtout(std_out,message,'COLL')
2349 : end if ! abs(optksloc)=3
2350 :
2351 : ! == Precision on charge_matlu (done only if local charge was computed ie not for optksloc=-1)
2352 4984 : if (icomp_chloc == 1 .and. paw_dmft%idmftloop >= 1 .and. optdiff == 1) then ! if the computation was done here.
2353 106 : if (green%has_charge_matlu_prev == 2) then
2354 : diff_chloc = zero
2355 62 : do iatom=1,natom
2356 49 : lpawu = paw_dmft%lpawu(iatom)
2357 49 : if (lpawu == -1) cycle
2358 87 : do isppol=1,nsppol
2359 : diff_chloc = diff_chloc + &
2360 98 : & (green%charge_matlu_prev(isppol,iatom)-green%charge_matlu(isppol,iatom))**2
2361 : end do ! isppol
2362 : end do ! iatom
2363 13 : if (sqrt(diff_chloc) < paw_dmft%dmft_lcpr) then
2364 0 : green%ichargeloc_cv = 1
2365 0 : write(message,'(a,8x,a,e9.2,a,8x,a)') ch10,"Change of correlated number of electron =<",&
2366 0 : & paw_dmft%dmft_lcpr,ch10,"DMFT Loop: Local charge is converged"
2367 0 : call wrtout(std_out,message,'COLL')
2368 : else
2369 13 : green%ichargeloc_cv = 0
2370 13 : write(message,'(a,8x,a,e9.2,a,8x,a)') ch10,"Change of correlated number of electron >",&
2371 26 : & paw_dmft%dmft_lcpr,ch10,"DMFT Loop: Local charge is not converged"
2372 13 : call wrtout(std_out,message,'COLL')
2373 : end if ! sqrt(diff_chloc)<dmft_lcpr
2374 : end if ! green%has_charge_matlu_prev=2
2375 1370 : green%charge_matlu_prev(:,:) = green%charge_matlu(:,:)
2376 106 : green%has_charge_matlu_prev = 2
2377 : end if ! icomp_chloc=1 and present(opt_diff) and idmftloop>=1
2378 :
2379 4984 : call destroy_matlu(matlu_temp(:),natom)
2380 23716 : ABI_FREE(matlu_temp)
2381 : ! deallocate(charge_loc_old)
2382 : ABI_NVTX_END_RANGE()
2383 4984 : call timab(625,2,tsec(:))
2384 : DBG_EXIT("COLL")
2385 :
2386 4984 : end subroutine integrate_green
2387 : !!***
2388 :
2389 : !!****f* m_green/icip_green
2390 : !! NAME
2391 : !! icip_green
2392 : !!
2393 : !! FUNCTION
2394 : !! init, compute, integrate and print lda green function
2395 : !!
2396 : !! INPUTS
2397 : !! char1 = character which precises the type of green function computed.
2398 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
2399 : !! green <type(green_type)>= green function data
2400 : !! pawprtvol = option for printing
2401 : !! self <type(self_type)>= variables related to self-energy
2402 : !! opt_self = optional argument, if =1, upfold self-energy
2403 : !! opt_moments = 1 to activate the computation of high-frequency moments
2404 : !! opt_log = 1 to compute Tr(log(G))
2405 : !!
2406 : !! OUTPUT
2407 : !!
2408 : !! SOURCE
2409 :
2410 279 : subroutine icip_green(char1,green,paw_dmft,pawprtvol,self,opt_self,opt_moments,opt_log)
2411 :
2412 : !Arguments ------------------------------------
2413 : !type(MPI_type), intent(in) :: mpi_enreg
2414 : type(green_type), intent(inout) :: green
2415 : type(paw_dmft_type), intent(inout) :: paw_dmft
2416 : type(self_type), intent(inout) :: self
2417 : integer, intent(in) :: pawprtvol
2418 : character(len=*), intent(in) :: char1
2419 : integer, optional, intent(in) :: opt_log,opt_moments,opt_self
2420 : !Local variables ------------------------------------
2421 : integer :: option,optlog,optmoments,optself,optlocks,prtopt_for_integrate_green,opt_nonxsum_icip
2422 : character(len=500) :: message
2423 : ! *********************************************************************
2424 :
2425 : !green%whichgreen="DFT"
2426 279 : prtopt_for_integrate_green = 2
2427 :
2428 279 : optself = 0
2429 279 : if (present(opt_self)) optself = opt_self
2430 279 : opt_nonxsum_icip = 1
2431 : !if(paw_dmft%dmftcheck==1) then ! for fourier_green
2432 : ! opt_nonxsum_icip=0
2433 : !endif
2434 :
2435 279 : optlog = 0
2436 279 : if (present(opt_log)) optlog = opt_log
2437 :
2438 279 : optmoments = 0
2439 279 : if (present(opt_moments)) optmoments = opt_moments
2440 :
2441 279 : call init_green(green,paw_dmft,opt_moments=optmoments)
2442 :
2443 : ! useless test ok
2444 : ! call printocc_green(green,1,paw_dmft,2)
2445 : ! write(std_out,*)" printocc_green zero finished "
2446 :
2447 : !== Compute green%oper(:)%ks
2448 : !== Deduce green%oper(:)%matlu(:)%mat
2449 279 : call compute_green(green,paw_dmft,pawprtvol,self,opt_self=optself,opt_nonxsum=opt_nonxsum_icip,opt_log=optlog)
2450 279 : if (paw_dmft%dmft_prgn >= 1 .and. paw_dmft%dmft_prgn <= 2) then
2451 0 : optlocks = paw_dmft%dmft_prgn*2 + 1 ! if dmft_prgn==2 => do not print
2452 0 : if (paw_dmft%lchipsiortho == 1 .and. pawprtvol > -100) then
2453 0 : call print_green(char1,green,optlocks,paw_dmft)
2454 : ! call print_green('inicip',green,1,paw_dmft,pawprtvol=1,opt_wt=1)
2455 : end if
2456 : end if
2457 :
2458 : !== Integrate green%oper(:)%ks
2459 : !== Integrate green%oper(:)%matlu(:)%mat
2460 279 : call integrate_green(green,paw_dmft,prtopt_for_integrate_green,opt_ksloc=3,opt_self=optself)!,opt_nonxsum=opt_nonxsum_icip)
2461 : !== Print green%oper(:)%ks
2462 : !== Print green%oper(:)%matlu(:)%mat
2463 279 : if (char1 == "DFT") then
2464 93 : option = 1
2465 93 : if (self%oper(1)%matlu(1)%lpawu /= -1) then
2466 93 : if (abs(dble(self%oper(1)%matlu(1)%mat(1,1,1))) > tol7) then
2467 : ! todo_ab: generalise this
2468 0 : write(message,'(2a)') ch10,&
2469 0 : & "Warning: a DFT calculation is carried out and self is not zero"
2470 0 : ABI_WARNING(message)
2471 : ! call abi_abort('COLL')
2472 : end if
2473 : end if
2474 : else
2475 186 : option = 5
2476 : end if ! char1
2477 :
2478 279 : if (pawprtvol > -100) then
2479 279 : call printocc_green(green,option,paw_dmft,3,chtype=char1)
2480 : end if
2481 :
2482 279 : end subroutine icip_green
2483 : !!***
2484 :
2485 : !!****f* m_green/fourier_green
2486 : !! NAME
2487 : !! fourier_green
2488 : !!
2489 : !! FUNCTION
2490 : !! integrate green function in the band index basis
2491 : !!
2492 : !! INPUTS
2493 : !! cryst_struc <type(crystal_t)>= crystal structure data.
2494 : !! green <type(green_type)>= green function data
2495 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
2496 : !!
2497 : !! OUTPUT
2498 : !!
2499 : !! SOURCE
2500 :
2501 2 : subroutine fourier_green(cryst_struc,green,paw_dmft,opt_ksloc,opt_tw)
2502 :
2503 : !Arguments ------------------------------------
2504 : !type
2505 : type(crystal_t),intent(in) :: cryst_struc
2506 : type(green_type),intent(inout) :: green
2507 : !type(MPI_type), intent(in) :: mpi_enreg
2508 : type(paw_dmft_type), intent(in) :: paw_dmft
2509 : integer,intent(in) :: opt_ksloc,opt_tw ! fourier on ks or local
2510 :
2511 : !local variables-------------------------------
2512 : integer :: iatom,ib,ib1,ierr,ifreq,ikpt,im,im1,iparal,is,ispinor,ispinor1,itau
2513 : integer :: mband,mbandc,myproc,natom,ndim,nkpt,nproc,nspinor,nsppol,spacecomm!,opt_four
2514 : character(len=500) :: message
2515 : ! complex(dp):: ybcbeg,ybcend
2516 : ! arrays
2517 2 : complex(dp), allocatable :: fw(:)
2518 2 : complex(dp), allocatable :: ft(:)
2519 2 : type(green_type) :: green_temp
2520 : ! *********************************************************************
2521 : ! ybcbeg=czero
2522 : ! ybcend=czero
2523 :
2524 : ! define spaceComm, myproc, and nproc from world communicator
2525 : ! and mpi_enreg
2526 2 : spacecomm=paw_dmft%spacecomm
2527 2 : myproc=paw_dmft%myproc
2528 2 : nproc=paw_dmft%nproc
2529 :
2530 : ! Initialise integers
2531 2 : mband = paw_dmft%mband
2532 2 : mbandc = paw_dmft%mbandc
2533 2 : nkpt = paw_dmft%nkpt
2534 2 : nspinor = paw_dmft%nspinor
2535 2 : nsppol = paw_dmft%nsppol
2536 2 : natom = cryst_struc%natom
2537 :
2538 : ! Only imaginary frequencies here
2539 2 : if(green%w_type=="real") then
2540 0 : message = 'fourier_green not implemented for real frequency'
2541 0 : ABI_BUG(message)
2542 : endif
2543 :
2544 : ! Initialise temporary green function
2545 2 : call init_green(green_temp,paw_dmft)
2546 2 : call init_green_tau(green_temp,paw_dmft)
2547 :
2548 : !green%oper(:)%matlu(1)%mat(1,1,1,1,1)
2549 6 : ABI_MALLOC(fw,(green%nw))
2550 6 : ABI_MALLOC(ft,(green%dmftqmc_l))
2551 :
2552 : !==============================================
2553 : ! == Inverse fourier transformation ===========
2554 : !==============================================
2555 :
2556 2 : if(opt_tw==-1) then
2557 :
2558 : ! = For Kohn Sham green function
2559 : !==============================================
2560 1 : if(opt_ksloc==1.and.green%use_oper_tau_ks==1) then
2561 :
2562 0 : do is = 1 , nsppol
2563 0 : do ikpt = 1, nkpt
2564 0 : do ib = 1, mbandc
2565 0 : do ib1 = 1, mbandc
2566 0 : do ifreq=1,green%nw
2567 0 : fw(ifreq)=green%oper(ifreq)%ks(ib,ib1,ikpt,is)
2568 : enddo
2569 0 : call fourier_fct(fw,ft,ib==ib1,green%dmftqmc_l,-1,paw_dmft) ! inverse fourier
2570 0 : do itau=1,green%dmftqmc_l
2571 0 : green_temp%oper_tau(itau)%ks(ib,ib1,ikpt,is)=ft(itau)
2572 : enddo
2573 0 : if(ib==ib1) then
2574 0 : green%occup_tau%ks(ib,ib1,ikpt,is)=ft(1)+one
2575 : else
2576 0 : green%occup_tau%ks(ib,ib1,ikpt,is)=ft(1)
2577 : endif
2578 : enddo ! ib1
2579 : enddo ! ib
2580 : enddo ! ikpt
2581 : enddo ! isppol
2582 :
2583 : ! = Post-treatment: necessary in the case of nspinor==2, but valid anywhere
2584 : ! because G(tau)=G_{nu,nu'}(tau)+[G_{nu,nu'}(tau)]*
2585 :
2586 0 : do is = 1 , nsppol
2587 0 : do ikpt = 1, nkpt
2588 0 : do ib = 1, mbandc
2589 0 : do ib1 = 1, mbandc
2590 0 : do itau=1,green%dmftqmc_l
2591 : green%oper_tau(itau)%ks(ib,ib1,ikpt,is)=&
2592 : & (green_temp%oper_tau(itau)%ks(ib,ib1,ikpt,is)+ &
2593 0 : & conjg(green_temp%oper_tau(itau)%ks(ib1,ib,ikpt,is)))/two
2594 0 : if(ib==ib1) then
2595 0 : green%occup_tau%ks(ib,ib1,ikpt,is)=green%oper_tau(1)%ks(ib,ib1,ikpt,is)+one
2596 : else
2597 0 : green%occup_tau%ks(ib,ib1,ikpt,is)=green%oper_tau(1)%ks(ib,ib1,ikpt,is)
2598 : endif
2599 : enddo
2600 : enddo ! ib1
2601 : enddo ! ib
2602 : enddo ! ikpt
2603 : enddo ! isppol
2604 0 : call downfold_oper(green%occup_tau,paw_dmft)
2605 0 : call sym_matlu(green%occup_tau%matlu,paw_dmft)
2606 0 : write(message,'(a,a,i10,a)') ch10," green%occup_tau%matlu from green_occup_tau%ks"
2607 0 : call wrtout(std_out,message,'COLL')
2608 0 : call print_matlu(green%occup_tau%matlu,natom,prtopt=3)
2609 : endif
2610 :
2611 : ! = For local green function
2612 : !==============================================
2613 1 : if(opt_ksloc ==2) then
2614 :
2615 : iparal=0
2616 3 : do iatom=1, natom
2617 3 : if(green%oper(1)%matlu(iatom)%lpawu.ne.-1) then
2618 1 : ndim=2*green%oper(1)%matlu(iatom)%lpawu+1
2619 3 : do is = 1 , nsppol
2620 5 : do ispinor = 1, nspinor
2621 6 : do ispinor1 = 1, nspinor
2622 14 : do im=1,ndim
2623 62 : do im1=1,ndim
2624 50 : iparal=iparal+1
2625 60 : if(mod(iparal-1,nproc)==myproc) then
2626 1650 : do ifreq=1,green%nw
2627 1650 : fw(ifreq)=green%oper(ifreq)%matlu(iatom)%mat(im+(ispinor-1)*ndim,im1+(ispinor1-1)*ndim,is)
2628 : enddo
2629 : ! inverse fourier
2630 : ! write(std_out,'(a)') "fourierbeforeposttreatement,ispinor,ispinor1,is,im,im1"
2631 : ! write(std_out,'(a,5i4,f12.5,f12.5)') "fourier",ispinor,ispinor1,is,im,im1
2632 : ! write(std_out,'(a,e12.5,e12.5)')&
2633 : ! &"green%oper(4)%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1)"&
2634 : !& ,green%oper(4)%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1)
2635 50 : call fourier_fct(fw,ft,(im==im1).and.(ispinor==ispinor1),green%dmftqmc_l,-1,paw_dmft)
2636 3250 : do itau=1,green%dmftqmc_l
2637 3250 : green_temp%oper_tau(itau)%matlu(iatom)%mat(im+(ispinor-1)*ndim,im1+(ispinor1-1)*ndim,is)=ft(itau)
2638 : ! write(std_out,*) itau,green_temp%oper_tau(itau)%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1)
2639 : enddo
2640 : ! if((im==im1).and.(ispinor==ispinor1)) then
2641 : ! green%occup_tau%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1)=ft(1)+one
2642 : ! else
2643 : ! green%occup_tau%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1)=ft(1)
2644 : ! endif
2645 : endif ! iparal
2646 : enddo
2647 : enddo
2648 : enddo ! ispinor1
2649 : enddo ! ispinor
2650 : enddo ! isppol
2651 : endif ! lpawu.ne.-1
2652 : enddo ! iatom
2653 :
2654 : ! Parallelisation must be finished here, because in the post
2655 : ! treatment, value from different proc will be mixed.
2656 1 : call xmpi_barrier(spacecomm)
2657 3 : do iatom=1,natom
2658 2 : if (green%oper(1)%matlu(iatom)%lpawu==-1) cycle
2659 66 : do itau=1,green%dmftqmc_l
2660 66 : call xmpi_sum(green_temp%oper_tau(itau)%matlu(iatom)%mat,spacecomm,ierr)
2661 : enddo
2662 : enddo
2663 1 : call xmpi_barrier(spacecomm)
2664 :
2665 : ! = Post-treatment: necessary in the case of nspinor==2, but valid anywhere
2666 : ! because G(tau)=G_{LL'}^{sigma,sigma'}(tau)+[G_{L'L}^{sigma',sigma}(tau)]*
2667 :
2668 1 : if(nspinor>0) Then
2669 3 : do iatom=1, natom
2670 3 : if(green%oper(1)%matlu(iatom)%lpawu.ne.-1) then
2671 1 : ndim=2*green%oper(1)%matlu(iatom)%lpawu+1
2672 3 : do is = 1 , nsppol
2673 5 : do ispinor = 1, nspinor
2674 6 : do ispinor1 = 1, nspinor
2675 14 : do im=1,ndim
2676 62 : do im1=1,ndim
2677 : ! write(std_out,'(a,5i4,f12.5,f12.5)') "fourier -1",ispinor,ispinor1,is,im,im1
2678 3260 : do itau=1,green%dmftqmc_l
2679 : green%oper_tau(itau)%matlu(iatom)%mat(im+(ispinor-1)*ndim,im1+(ispinor1-1)*ndim,is)=&
2680 : & ((green_temp%oper_tau(itau)%matlu(iatom)%mat(im+(ispinor-1)*ndim,im1+(ispinor1-1)*ndim,is)+ &
2681 3200 : & conjg(green_temp%oper_tau(itau)%matlu(iatom)%mat(im1+(ispinor1-1)*ndim,im+(ispinor-1)*ndim,is))))/two
2682 : ! write(std_out,*) itau,green%oper_tau(itau)%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1)
2683 3250 : if((im==im1).and.(ispinor==ispinor1)) then
2684 : green%occup_tau%matlu(iatom)%mat(im+(ispinor-1)*ndim,im1+(ispinor1-1)*ndim,is)=&
2685 640 : & green%oper_tau(1)%matlu(iatom)%mat(im+(ispinor-1)*ndim,im1+(ispinor1-1)*ndim,is)+one
2686 : else
2687 : green%occup_tau%matlu(iatom)%mat(im+(ispinor-1)*ndim,im1+(ispinor1-1)*ndim,is)=&
2688 2560 : & green%oper_tau(1)%matlu(iatom)%mat(im+(ispinor-1)*ndim,im1+(ispinor1-1)*ndim,is)
2689 : endif ! test diag
2690 : enddo ! itau
2691 : enddo ! im1
2692 : enddo ! im
2693 : enddo ! ispinor1
2694 : enddo ! ispinor
2695 : enddo ! isppol
2696 : endif
2697 : enddo ! iatom
2698 : endif ! nspinor>0
2699 : endif ! opt_ksloc=2
2700 : endif ! opt_tw==-1
2701 :
2702 :
2703 : !==============================================
2704 : ! == Direct fourier transformation
2705 : !==============================================
2706 : ! todo_ba ft useful only for diagonal elements ...
2707 :
2708 2 : if(opt_tw==1) then
2709 :
2710 : ! = For local green function
2711 : !==============================================
2712 1 : if(opt_ksloc ==2) then
2713 :
2714 : iparal=0
2715 3 : do iatom=1, natom
2716 : ! put to zero (usefull at least for parallelism)
2717 3 : if(green%oper(1)%matlu(iatom)%lpawu.ne.-1) then
2718 33 : do ifreq=1,green%nw
2719 2017 : green%oper(ifreq)%matlu(iatom)%mat=czero
2720 : enddo
2721 1 : ndim=2*green%oper(1)%matlu(iatom)%lpawu+1
2722 3 : do is = 1 , nsppol
2723 5 : do ispinor = 1, nspinor
2724 6 : do ispinor1 = 1, nspinor
2725 14 : do im=1,ndim
2726 62 : do im1=1,ndim
2727 50 : iparal=iparal+1
2728 60 : if(mod(iparal-1,nproc)==myproc) then
2729 3250 : do itau=1,green%dmftqmc_l
2730 3250 : ft(itau)=green%oper_tau(itau)%matlu(iatom)%mat(im+(ispinor-1)*ndim,im1+(ispinor1-1)*ndim,is)
2731 : enddo
2732 : ! if(im==1.and.im1==1) write(std_out,*) "ft(itau=1)",is,ft(1) ! debug
2733 50 : call fourier_fct(fw,ft,(im==im1).and.(ispinor==ispinor1),green%dmftqmc_l,1,paw_dmft)
2734 1650 : do ifreq=1,green%nw
2735 1650 : green%oper(ifreq)%matlu(iatom)%mat(im+(ispinor-1)*ndim,im1+(ispinor1-1)*ndim,is)=fw(ifreq)
2736 : enddo
2737 : ! if(im==1.and.im1==1) write(std_out,*) "fw(ifreq=1)",is,fw(1) ! debug
2738 : endif
2739 : enddo
2740 : enddo
2741 : enddo ! ispinor1
2742 : enddo ! ispinor
2743 : enddo ! isppol
2744 : endif ! lpawu=/-1
2745 : enddo ! iatom
2746 1 : call xmpi_barrier(spacecomm)
2747 3 : do iatom=1,natom
2748 2 : if (green%oper(1)%matlu(iatom)%lpawu==-1) cycle
2749 34 : do ifreq=1,green%nw
2750 34 : call xmpi_sum(green%oper(ifreq)%matlu(iatom)%mat,spacecomm,ierr)
2751 : enddo
2752 : enddo
2753 : endif ! opt_ksloc=2
2754 : endif ! opt_tw==-1
2755 2 : ABI_FREE(fw)
2756 2 : ABI_FREE(ft)
2757 2 : call destroy_green_tau(green_temp)
2758 2 : call destroy_green(green_temp)
2759 :
2760 2 : end subroutine fourier_green
2761 : !!***
2762 :
2763 :
2764 : !!****f* m_green/check_fourier_green
2765 : !! NAME
2766 : !! check_fourier_green
2767 : !!
2768 : !! FUNCTION
2769 : !! Check fourier transformations
2770 : !!
2771 : !! INPUTS
2772 : !! cryst_struc <type(crystal_t)>= crystal structure data.
2773 : !! green <type(green_type)>= green function data
2774 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
2775 : !! pawang <type(pawang)>=paw angular mesh and related data
2776 : !!
2777 : !! OUTPUT
2778 : !!
2779 : !! SOURCE
2780 :
2781 1 : subroutine check_fourier_green(cryst_struc,green,paw_dmft)
2782 :
2783 : !Arguments ------------------------------------
2784 : !type
2785 : type(crystal_t),intent(in) :: cryst_struc
2786 : type(green_type),intent(inout) :: green
2787 : !type(MPI_type), intent(in) :: mpi_enreg
2788 : type(paw_dmft_type), intent(inout) :: paw_dmft
2789 :
2790 : !local variables-------------------------------
2791 1 : type(green_type) :: green_check
2792 : character(len=500) :: message
2793 : ! *********************************************************************
2794 : ! Only imaginary frequencies here
2795 1 : if(green%w_type=="real") then
2796 0 : message = 'check_fourier_green not implemented for real frequency'
2797 0 : ABI_BUG(message)
2798 : endif
2799 :
2800 1 : call init_green(green_check,paw_dmft)
2801 1 : call init_green_tau(green_check,paw_dmft)
2802 1 : call copy_green(green,green_check,opt_tw=2)
2803 :
2804 1 : write(message,'(2a,i3,13x,a)') ch10,' === Inverse Fourier Transform w->t of Weiss Field'
2805 1 : call wrtout(std_out,message,'COLL')
2806 : call fourier_green(cryst_struc,green_check,paw_dmft,&
2807 1 : & opt_ksloc=2,opt_tw=-1)
2808 :
2809 1 : write(message,'(3a)') ch10,' === Print (for check by user) of occupation matrix'&
2810 2 : & ,' after fourier transform with respect to initial one'
2811 1 : call wrtout(std_out,message,'COLL')
2812 1 : call printocc_green(green_check,6,paw_dmft,3)
2813 :
2814 1 : write(message,'(2a,i3,13x,a)') ch10,' === Direct Fourier Transform t->w of Weiss Field'
2815 1 : call wrtout(std_out,message,'COLL')
2816 : call fourier_green(cryst_struc,green_check,paw_dmft,&
2817 1 : & opt_ksloc=2,opt_tw=1)
2818 : ! call print_matlu(green%oper(1)%matlu,paw_dmft%natom,1) ! debug
2819 :
2820 : call integrate_green(green_check,paw_dmft,&
2821 1 : & prtopt=2,opt_ksloc=2)
2822 :
2823 1 : write(message,'(3a)') ch10,' === Print (for check by user) of occupation matrix'&
2824 2 : & ,' after double fourier transform with respect to initial one'
2825 1 : call wrtout(std_out,message,'COLL')
2826 1 : call printocc_green(green_check,5,paw_dmft,3)
2827 :
2828 1 : call destroy_green_tau(green_check)
2829 1 : call destroy_green(green_check)
2830 :
2831 1 : end subroutine check_fourier_green
2832 : !!***
2833 :
2834 : !!****f* m_green/compa_occup_ks
2835 : !! NAME
2836 : !! compa_occup_ks
2837 : !!
2838 : !! FUNCTION
2839 : !! Compare occupations to Fermi Dirac Occupations
2840 : !!
2841 : !! INPUTS
2842 : !! green <type(green_type)>= green function data
2843 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
2844 : !!
2845 : !! OUTPUT
2846 : !!
2847 : !! SOURCE
2848 :
2849 0 : subroutine compa_occup_ks(green,paw_dmft)
2850 :
2851 : !Arguments ------------------------------------
2852 : !type
2853 : type(green_type),intent(inout) :: green
2854 : type(paw_dmft_type), intent(inout) :: paw_dmft
2855 :
2856 : !local variables-------------------------------
2857 : character(len=500) :: message
2858 : integer :: ib,ikpt,isppol
2859 : integer :: ib_,ikpt_,isppol_
2860 : integer :: ib__,ikpt__,isppol__
2861 : real(dp) :: diffmax,occ1,occ2,occ_a,occ_b,diffrel,occ_aa,occ_bb
2862 : ! *********************************************************************
2863 0 : diffmax=zero
2864 0 : diffrel=zero
2865 0 : do isppol=1,paw_dmft%nsppol
2866 0 : do ikpt=1,paw_dmft%nkpt
2867 0 : do ib=1,paw_dmft%mbandc
2868 0 : occ1=occup_fd(paw_dmft%eigen_dft(ib,ikpt,isppol),paw_dmft%fermie,paw_dmft%temp)
2869 0 : occ2=real(green%occup%ks(ib,ib,ikpt,isppol))
2870 0 : if(abs(occ1-occ2)>diffmax) then
2871 0 : diffmax=abs(occ1-occ2)
2872 0 : occ_a=occ1
2873 0 : occ_b=occ2
2874 0 : ib_=ib;isppol_=isppol;ikpt_=ikpt
2875 : endif
2876 0 : if(abs(two*(occ1-occ2)/(occ1+occ2))>diffrel) then
2877 0 : diffrel=abs(two*(occ1-occ2)/(occ1+occ2))
2878 0 : occ_aa=occ1
2879 0 : occ_bb=occ2
2880 0 : ib__=ib;isppol__=isppol;ikpt__=ikpt
2881 : endif
2882 : enddo
2883 : enddo
2884 : enddo
2885 :
2886 :
2887 0 : write(message,'(2a)') ch10,' === Compare green function occupations and Fermi Dirac occupations'
2888 0 : call wrtout(std_out,message,'COLL')
2889 0 : write(message,'(2a,f12.5,2a,f12.5,2a,f12.5,2a,3i5)') ch10,' = Max difference is',diffmax,&
2890 0 : & ch10,' Corresponding Occupation from green function is',occ_b,&
2891 0 : & ch10,' Corresponding Occupation Fermi Dirac weight is',occ_a,&
2892 0 : & ch10,' (For polarization, k-point and band index) ',isppol_,ikpt_,ib_
2893 0 : call wrtout(std_out,message,'COLL')
2894 0 : write(message,'(2a,f12.5,2a,f12.5,2a,f12.5,2a,3i5)') ch10,' = Max relative difference is',diffrel,&
2895 0 : & ch10,' Corresponding Occupation from green function is',occ_bb,&
2896 0 : & ch10,' Corresponding Occupation Fermi Dirac weight is',occ_aa,&
2897 0 : & ch10,' (For polarization, k-point and band index) ',isppol__,ikpt__,ib__
2898 0 : call wrtout(std_out,message,'COLL')
2899 :
2900 0 : end subroutine compa_occup_ks
2901 : !!***
2902 :
2903 : !!****f* m_green/add_int_fct
2904 : !! NAME
2905 : !! add_int_fct
2906 : !!
2907 : !! FUNCTION
2908 : !! Do integration in matsubara space
2909 : !!
2910 : !! COPYRIGHT
2911 : !! Copyright (C) 2006-2026 ABINIT group (BAmadon)
2912 : !! This file is distributed under the terms of the
2913 : !! GNU General Public License, see ~abinit/COPYING
2914 : !! or http://www.gnu.org/copyleft/gpl.txt .
2915 : !!
2916 : !! INPUTS
2917 : !! ff= function is frequency space
2918 : !! ldiag = option according to diagonal or non-diagonal elements
2919 : !! option = nspinor
2920 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
2921 : !! proct= for parallelism
2922 : !!
2923 : !! SIDE EFFECTS
2924 : !! * integral = integral of ff over matsubara frequencies (there is an accumulation in the present routine, so intent(inout))
2925 : !! * ft= function is time space
2926 : !!
2927 : !! SOURCE
2928 0 : subroutine add_int_fct(ifreq,ff,ldiag,omega_current,option,integral,temp,wgt_wlo,dmft_nwlo)
2929 :
2930 : !Arguments ------------------------------------
2931 : !type
2932 : integer,intent(in) :: ifreq
2933 : logical,intent(in) :: ldiag
2934 : integer,intent(in) :: option,dmft_nwlo
2935 : complex(dp),intent(inout) :: integral
2936 : complex(dp), intent(in) :: ff
2937 : complex(dp), intent(in) :: omega_current
2938 : real(dp), intent(in) :: temp, wgt_wlo
2939 :
2940 : !local variables-------------------------------
2941 : real(dp) :: omega
2942 : ! *********************************************************************
2943 0 : omega=aimag(omega_current)
2944 0 : if(ldiag) then
2945 :
2946 0 : if(option==1) then ! nspinor==1
2947 : ! write(500,*) paw_dmft%omega_lo(ifreq),real(ff(ifreq)),imag(ff(ifreq))
2948 : integral=integral+2.d0*temp * &
2949 : & real( ff-one / ( j_dpc*omega ) ) * &
2950 0 : & wgt_wlo
2951 0 : if(ifreq==dmft_nwlo) integral=integral+half
2952 : ! integral=integral+half
2953 : ! the if is here, to count only one time this correction
2954 : endif
2955 :
2956 0 : if(option==2) then ! nspinor==2
2957 : integral=integral+2.d0*temp * &
2958 : & ( ff-one / ( j_dpc*omega ) ) * &
2959 0 : & wgt_wlo
2960 0 : if(ifreq==dmft_nwlo) integral=integral+half
2961 : ! integral=integral+half
2962 : ! the if is here, to count only one time this correction
2963 : endif
2964 :
2965 :
2966 : else ! ldiag
2967 :
2968 : ! write(std_out,*) "nondiag"
2969 0 : if(option==1) then
2970 : integral=integral+2.d0*temp * &
2971 : & real( ff ) * &
2972 0 : & wgt_wlo
2973 : endif
2974 0 : if(option==2) then
2975 : integral=integral+2.d0*temp * &
2976 : & ff * &
2977 0 : & wgt_wlo
2978 : endif
2979 : endif ! ldiag
2980 :
2981 0 : end subroutine add_int_fct
2982 : !!***
2983 :
2984 : !!****m* m_green/int_fct
2985 : !! NAME
2986 : !! int_fct
2987 : !!
2988 : !! FUNCTION
2989 : !! Do integration in matsubara space
2990 : !!
2991 : !! COPYRIGHT
2992 : !! Copyright (C) 2006-2026 ABINIT group (BAmadon)
2993 : !! This file is distributed under the terms of the
2994 : !! GNU General Public License, see ~abinit/COPYING
2995 : !! or http://www.gnu.org/copyleft/gpl.txt .
2996 : !!
2997 : !! INPUTS
2998 : !! ff= function is frequency space
2999 : !! ldiag = option according to diagonal or non-diagonal elements
3000 : !! option = nspinor
3001 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
3002 : !! proct= for parallelism
3003 : !!
3004 : !! OUTPUT
3005 : !! integral = integral of ff over matsubara frequencies
3006 : !!
3007 : !! SIDE EFFECTS
3008 : !! ft= function is time space
3009 : !!
3010 : !! SOURCE
3011 0 : subroutine int_fct(ff,ldiag,option,paw_dmft,integral,procb,myproc)
3012 :
3013 : !Arguments ------------------------------------
3014 : !type
3015 : logical,intent(in) :: ldiag
3016 : integer,intent(in) :: option
3017 : complex(dp),intent(out) :: integral
3018 : type(paw_dmft_type), intent(in) :: paw_dmft
3019 : complex(dp), intent(in) :: ff(paw_dmft%dmft_nwlo)
3020 : integer, optional, intent(in) :: procb(paw_dmft%dmft_nwlo)
3021 : integer, optional, intent(in) :: myproc
3022 :
3023 : !local variables-------------------------------
3024 0 : logical, allocatable :: procb2(:)
3025 : character(len=500) :: message
3026 : integer :: ifreq
3027 : ! *********************************************************************
3028 0 : ABI_MALLOC(procb2,(paw_dmft%dmft_nwlo))
3029 0 : if(present(procb).and.present(myproc)) then
3030 0 : do ifreq=1,paw_dmft%dmft_nwlo
3031 0 : procb2(ifreq)=(procb(ifreq)==myproc)
3032 : enddo
3033 0 : else if(present(procb).and..not.present(myproc)) then
3034 0 : write(message,'(a,a,2(e15.4))') ch10,&
3035 0 : & "BUG: procb is present and not myproc in int_fct"
3036 0 : ABI_BUG(message)
3037 0 : else if(.not.present(procb).and.present(myproc)) then
3038 0 : write(message,'(a,a,2(e15.4))') ch10,&
3039 0 : & "BUG: procb is not present and myproc is in int_fct"
3040 0 : ABI_BUG(message)
3041 : else
3042 0 : do ifreq=1,paw_dmft%dmft_nwlo
3043 0 : procb2(ifreq)=(1==1)
3044 : enddo
3045 : endif
3046 :
3047 0 : integral=czero
3048 0 : if(ldiag) then
3049 :
3050 0 : if(option==1) then ! nspinor==1
3051 0 : do ifreq=1,paw_dmft%dmft_nwlo
3052 0 : if(procb2(ifreq)) then
3053 : ! write(500,*) paw_dmft%omega_lo(ifreq),real(ff(ifreq)),imag(ff(ifreq))
3054 : integral=integral+2.d0*paw_dmft%temp * &
3055 : & real( ff(ifreq)-one / ( j_dpc*paw_dmft%omega_lo(ifreq) ) ) * &
3056 0 : & paw_dmft%wgt_wlo(ifreq)
3057 : endif
3058 : enddo
3059 0 : if(procb2(paw_dmft%dmft_nwlo)) integral=integral+half
3060 : ! integral=integral+half
3061 : ! the if is here, to count only one time this correction
3062 : endif
3063 :
3064 0 : if(option==2) then ! nspinor==2
3065 0 : do ifreq=1,paw_dmft%dmft_nwlo
3066 0 : if(procb2(ifreq)) then
3067 : integral=integral+2.d0*paw_dmft%temp * &
3068 : & ( ff(ifreq)-one / ( j_dpc*paw_dmft%omega_lo(ifreq) ) ) * &
3069 0 : & paw_dmft%wgt_wlo(ifreq)
3070 : endif
3071 : enddo
3072 0 : if(procb2(paw_dmft%dmft_nwlo)) integral=integral+half
3073 : ! integral=integral+half
3074 : ! the if is here, to count only one time this correction
3075 : endif
3076 :
3077 :
3078 : else ! ldiag
3079 :
3080 : ! write(std_out,*) "nondiag"
3081 0 : if(option==1) then
3082 0 : do ifreq=1,paw_dmft%dmft_nwlo
3083 0 : if(procb2(ifreq)) then
3084 : integral=integral+2.d0*paw_dmft%temp * &
3085 : & real( ff(ifreq) ) * &
3086 0 : & paw_dmft%wgt_wlo(ifreq)
3087 : endif
3088 : enddo
3089 : endif
3090 0 : if(option==2) then
3091 0 : do ifreq=1,paw_dmft%dmft_nwlo
3092 0 : if(procb2(ifreq)) then
3093 : integral=integral+2.d0*paw_dmft%temp * &
3094 : & ff(ifreq) * &
3095 0 : & paw_dmft%wgt_wlo(ifreq)
3096 : endif
3097 : enddo
3098 : endif
3099 : endif ! ldiag
3100 0 : ABI_FREE(procb2)
3101 :
3102 0 : end subroutine int_fct
3103 : !!***
3104 :
3105 : !!****f* m_green/fourier_fct
3106 : !! NAME
3107 : !! fourier_fct
3108 : !!
3109 : !! FUNCTION
3110 : !! Do fourier transformation from matsubara space to imaginary time
3111 : !! (A spline is performed )
3112 : !!
3113 : !! COPYRIGHT
3114 : !! Copyright (C) 2006-2026 ABINIT group (BAmadon)
3115 : !! This file is distributed under the terms of the
3116 : !! GNU General Public License, see ~abinit/COPYING
3117 : !! or http://www.gnu.org/copyleft/gpl.txt .
3118 : !!
3119 : !! INPUTS
3120 : !! ldiag = option according to diagonal or non-diagonal elements
3121 : !! opt_four = option for direct (1) or inverse (-1) fourier transform
3122 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
3123 : !!
3124 : !! OUTPUT
3125 : !!
3126 : !! SIDE EFFECTS
3127 : !! fw= function is frequency space
3128 : !! ft= function is time space
3129 : !!
3130 : !! SOURCE
3131 100 : subroutine fourier_fct(fw,ft,ldiag,ltau,opt_four,paw_dmft)
3132 :
3133 : !Arguments ------------------------------------
3134 : !type
3135 : logical,intent(in) :: ldiag
3136 : integer,intent(in) :: ltau,opt_four
3137 : type(paw_dmft_type), intent(in) :: paw_dmft
3138 : complex(dp), intent(inout) :: fw(paw_dmft%dmft_nwlo)
3139 : complex(dp), intent(inout) :: ft(ltau)
3140 :
3141 : !local variables-------------------------------
3142 100 : complex(dp), allocatable :: splined_li(:)
3143 100 : complex(dp), allocatable :: tospline_li(:)
3144 : ! complex(dp), allocatable :: fw1(:)
3145 100 : real(dp), allocatable :: ftr(:)
3146 : real(dp) :: beta
3147 : complex(dp) :: xsto
3148 : integer :: iflag,ifreq,itau,iwarn,log_direct
3149 : character(len=500) :: message
3150 100 : real(dp), allocatable :: omega_li(:)
3151 : ! *********************************************************************
3152 100 : beta=one/paw_dmft%temp
3153 100 : iflag=0
3154 100 : log_direct=1
3155 :
3156 100 : if(ldiag) iflag=1
3157 :
3158 : ! == inverse fourier transform
3159 100 : if(opt_four==-1) then
3160 :
3161 150 : ABI_MALLOC(splined_li,(paw_dmft%dmft_nwli))
3162 : ! allocate(fw1(0:paw_dmft%dmft_nwlo-1))
3163 50 : if(paw_dmft%dmft_log_freq==1) then
3164 150 : ABI_MALLOC(omega_li,(1:paw_dmft%dmft_nwli))
3165 50 : call construct_nwli_dmft(paw_dmft,paw_dmft%dmft_nwli,omega_li)
3166 : call spline_c(paw_dmft%dmft_nwlo,paw_dmft%dmft_nwli,paw_dmft%omega_lo,&
3167 50 : & omega_li,splined_li,fw)
3168 50 : ABI_FREE(omega_li)
3169 : else
3170 0 : splined_li=fw
3171 : endif
3172 50 : call invfourier(splined_li,ft,paw_dmft%dmft_nwli,ltau,iflag,beta)
3173 : ! deallocate(fw1)
3174 50 : ABI_FREE(splined_li)
3175 :
3176 : ! == direct fourier transform
3177 50 : else if(opt_four==1) then
3178 :
3179 150 : ABI_MALLOC(ftr,(ltau))
3180 :
3181 3250 : iwarn=0
3182 3250 : do itau=1,ltau
3183 3250 : if(abs(aimag(ft(itau)))>tol12) then
3184 0 : if(ldiag) then
3185 0 : write(message,'(a,a,2(e15.4))') ch10,&
3186 0 : & "green function is not real in imaginary time space",ft(itau)
3187 0 : ABI_ERROR(message)
3188 : else
3189 0 : iwarn=iwarn+1
3190 0 : ftr(itau)=real(ft(itau))
3191 0 : xsto=ft(itau)
3192 : endif
3193 : else
3194 3200 : ftr(itau)=real(ft(itau))
3195 : endif
3196 : ! write(std_out,*) itau,ftr(itau)
3197 : enddo
3198 :
3199 150 : ABI_MALLOC(tospline_li,(paw_dmft%dmft_nwli))
3200 : if (log_direct==1) then
3201 : ! do not have a physical meaning..because the log frequency is not
3202 : ! one of the linear frequency.
3203 : ! if log frequencies are also one of linear frequency, it should be good
3204 1650 : do ifreq=1, paw_dmft%dmft_nwlo
3205 1650 : call nfourier2(ftr,fw(ifreq),iflag,paw_dmft%omega_lo(ifreq),ltau,beta)
3206 : enddo
3207 : else
3208 : call nfourier(ftr,tospline_li,iflag,paw_dmft%dmft_nwli-1,ltau,beta)
3209 : do ifreq=1,paw_dmft%dmft_nwli
3210 : write(1112,*) paw_dmft%temp*pi*real(2*ifreq-1,kind=dp),real(tospline_li(ifreq),kind=dp),aimag(tospline_li(ifreq))
3211 : !write(1112,*) paw_dmft%omega_li(ifreq),real(tospline_li(ifreq)),aimag(tospline_li(ifreq))
3212 : enddo
3213 : if(paw_dmft%dmft_log_freq==1) then
3214 : ABI_MALLOC(omega_li,(1:paw_dmft%dmft_nwli))
3215 : call construct_nwli_dmft(paw_dmft,paw_dmft%dmft_nwli,omega_li)
3216 : call spline_c(paw_dmft%dmft_nwli,paw_dmft%dmft_nwlo,omega_li,&
3217 : & paw_dmft%omega_lo,fw,tospline_li)
3218 : ABI_FREE(omega_li)
3219 : else
3220 : fw=tospline_li
3221 : endif
3222 : endif
3223 :
3224 50 : ABI_FREE(tospline_li)
3225 :
3226 50 : ABI_FREE(ftr)
3227 50 : if(iwarn>0) then
3228 0 : write(message,'(a,a,2(e15.4))') ch10,&
3229 0 : & "WARNING: off-diag green function is not real in imaginary time space",xsto
3230 0 : call wrtout(std_out,message,'COLL')
3231 : endif
3232 :
3233 : endif
3234 :
3235 100 : end subroutine fourier_fct
3236 : !!***
3237 :
3238 : !!****f* m_green/spline_fct
3239 : !! NAME
3240 : !! spline_fct
3241 : !!
3242 : !! FUNCTION
3243 : !! Do fourier transformation from matsubara space to imaginary time
3244 : !! (A spline is performed )
3245 : !!
3246 : !! COPYRIGHT
3247 : !! Copyright (C) 2006-2026 ABINIT group (BAmadon)
3248 : !! This file is distributed under the terms of the
3249 : !! GNU General Public License, see ~abinit/COPYING
3250 : !! or http://www.gnu.org/copyleft/gpl.txt .
3251 : !!
3252 : !! INPUTS
3253 : !! opt_spline = option for the spline
3254 : !! -1 log to linear.
3255 : !! 1 linear to log.
3256 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
3257 : !!
3258 : !! OUTPUT
3259 : !!
3260 : !! SIDE EFFECTS
3261 : !! fw= function is frequency space
3262 : !! ft= function is time space
3263 : !!
3264 : !! SOURCE
3265 :
3266 0 : subroutine spline_fct(fw1,fw2,opt_spline,paw_dmft)
3267 :
3268 : !Arguments ------------------------------------
3269 : !type
3270 : integer,intent(in) :: opt_spline
3271 : type(paw_dmft_type), intent(in) :: paw_dmft
3272 : complex(dp), intent(inout) :: fw1(:)
3273 : complex(dp), intent(inout) :: fw2(:)
3274 : integer :: size_fw1
3275 : integer :: size_fw2
3276 0 : real(dp), allocatable :: omega_li(:)
3277 : ! *********************************************************************
3278 :
3279 0 : size_fw1 = size(fw1)
3280 0 : size_fw2 = size(fw2)
3281 :
3282 : ! == inverse fourier transform
3283 0 : if(opt_spline==-1) then
3284 :
3285 : ! allocate(fw1(0:paw_dmft%dmft_nwlo-1))
3286 0 : if(paw_dmft%dmft_log_freq==1) then
3287 0 : ABI_MALLOC(omega_li,(1:size_fw2))
3288 0 : call construct_nwli_dmft(paw_dmft,size_fw2,omega_li)
3289 : call spline_c(size_fw1,size_fw2,paw_dmft%omega_lo(1:size_fw1),&
3290 0 : & omega_li(1:size_fw2),fw2,fw1)
3291 0 : ABI_FREE(omega_li)
3292 : else
3293 0 : fw2=fw1
3294 : endif
3295 :
3296 : ! == direct fourier transform
3297 0 : else if(opt_spline==1) then
3298 :
3299 :
3300 0 : if(paw_dmft%dmft_log_freq==1) then
3301 0 : ABI_MALLOC(omega_li,(1:size_fw2))
3302 0 : call construct_nwli_dmft(paw_dmft,size_fw2,omega_li)
3303 : call spline_c(size_fw2,size_fw1,omega_li(1:size_fw2),&
3304 0 : & paw_dmft%omega_lo(1:size_fw1),fw1,fw2)
3305 0 : ABI_FREE(omega_li)
3306 : else
3307 0 : fw1=fw2
3308 : endif
3309 :
3310 : endif
3311 :
3312 0 : end subroutine spline_fct
3313 : !!***
3314 :
3315 : !!****f* m_green/occup_green_tau
3316 : !! NAME
3317 : !! occup_green_tau
3318 : !!
3319 : !! FUNCTION
3320 : !! Compute occup_tau from green%oper_tau
3321 : !!
3322 : !! INPUTS
3323 : !! green <type(green_type)>= green function data
3324 : !!
3325 : !! OUTPUT
3326 : !!
3327 : !! SOURCE
3328 :
3329 68 : subroutine occup_green_tau(green)
3330 :
3331 : !Arguments ------------------------------------
3332 : type(green_type), intent(inout) :: green
3333 : !Local variables-------------------------------
3334 : integer :: natom
3335 68 : complex(dp), allocatable :: shift(:)
3336 : ! *********************************************************************
3337 :
3338 68 : natom = green%oper_tau(1)%natom
3339 :
3340 204 : ABI_MALLOC(shift,(natom))
3341 :
3342 192 : shift(:) = - cone
3343 :
3344 68 : call copy_matlu(green%oper_tau(1)%matlu(:),green%occup_tau%matlu(:),natom)
3345 68 : call shift_matlu(green%occup_tau%matlu(:),natom,shift(:))
3346 :
3347 68 : ABI_FREE(shift)
3348 :
3349 68 : end subroutine occup_green_tau
3350 : !!***
3351 :
3352 : !!****f* m_green/occupfd
3353 : !! NAME
3354 : !! occupfd
3355 : !!
3356 : !! FUNCTION
3357 : !!
3358 : !! INPUTS
3359 : !!
3360 : !! OUTPUT
3361 : !!
3362 : !!
3363 : !! SOURCE
3364 :
3365 : !function occupfd(eig,fermie,temp)
3366 :
3367 :
3368 : !Arguments ------------------------------------
3369 : !type
3370 : ! Integrate analytic tail 1/(iw-mu)
3371 : !real(dp),intent(in) :: eig,fermie,temp
3372 : !real(dp) :: occupfd
3373 : !Local variables-------------------------------
3374 : ! *********************************************************************
3375 :
3376 : !if ((eig-fermie) > zero) then
3377 : ! occupfd = exp(-(eig-fermie)/temp)/(one+exp(-(eig-fermie)/temp))
3378 : !else
3379 : ! occupfd = one / (one+exp((eig-fermie)/temp))
3380 : !end if
3381 :
3382 : !end function occupfd
3383 : !!***
3384 :
3385 : !!****f* m_green/distrib_paral
3386 : !! NAME
3387 : !! distrib_paral
3388 : !!
3389 : !! FUNCTION
3390 : !!
3391 : !! INPUTS
3392 : !! nw : number of frequencies
3393 : !! nkpt : number of k-points
3394 : !! nproc : number of procs
3395 : !!
3396 : !! OUTPUT
3397 : !! procb(iw,ikpt) : number of the proc that is in charge of the combination {iw,ikpt}.
3398 : !! proct(iw,iproc) : 1 if the frequency "iw" should be computed by the proc "iproc"
3399 : !! 0 if iw should not " "
3400 : !! careful: if proct=1, it does not mean that each combination {iw,ikpt} is treated by this proc.
3401 : !! SOURCE
3402 :
3403 0 : subroutine distrib_paral(nkpt,nproc,nw,nw_perproc,procb,proct)
3404 :
3405 : !Arguments ------------------------------------
3406 : !type
3407 : ! Integrate analytic tail 1/(iw-mu)
3408 : integer, intent(in) :: nw,nkpt,nproc
3409 : integer, intent(out) :: nw_perproc
3410 : integer, intent(out):: procb(nw,nkpt),proct(nw,0:nproc-1)
3411 : !Local variables-------------------------------
3412 : integer :: ikpt,iw,ir,ratio,m,iproc
3413 0 : integer, allocatable :: proca(:,:)
3414 : ! *********************************************************************
3415 :
3416 :
3417 0 : proct(:,:)=0
3418 0 : procb(:,:)=-1
3419 :
3420 0 : if(nproc.ge.2*nw) then
3421 : !write(6,*) "AA"
3422 0 : ratio=nproc/nw
3423 0 : ABI_MALLOC(proca,(nw,nproc/nw))
3424 0 : do ir=1,ratio
3425 0 : do iw=1,nw
3426 0 : proca(iw,ir)=iw+(ir-1)*nw-1
3427 0 : proct(iw,iw+(ir-1)*nw-1)=1
3428 : enddo
3429 : enddo
3430 : ! do iw=1,nw
3431 : ! write(6,*) " freq, procs"
3432 : ! write(6,*) iw,proca(iw,:)
3433 : ! enddo
3434 0 : do iw=1,nw
3435 0 : do ikpt=1,nkpt
3436 0 : if(ikpt.le.ratio) procb(iw,ikpt)=proca(iw,ikpt)
3437 0 : if(ikpt.ge.ratio) then
3438 0 : m=mod(ikpt,ratio)
3439 0 : if(m==0) then
3440 0 : procb(iw,ikpt)=proca(iw,ratio)
3441 : else
3442 0 : procb(iw,ikpt)=proca(iw,m)
3443 : endif
3444 : endif
3445 : ! write(6,*) " freq, k-point, proc"
3446 : ! write(6,*) iw,ikpt,procb(iw,ikpt)
3447 : enddo
3448 : enddo
3449 0 : nw_perproc=1
3450 0 : ABI_FREE(proca)
3451 :
3452 0 : else if (nproc.ge.nw) then
3453 : !write(6,*) "BB"
3454 0 : do iw=1,nw
3455 0 : procb(iw,:)= iw
3456 0 : proct(iw,iw)=1
3457 : enddo
3458 : ! do iw=1,nw
3459 : ! write(6,*) " freq, proc"
3460 : ! write(6,*) iw, procb(iw,1)
3461 : ! enddo
3462 0 : nw_perproc=1
3463 :
3464 : else if (nproc.le.nw) then
3465 : !write(6,*) "CC"
3466 0 : ratio=nw/nproc
3467 : !write(6,*) "ratio", ratio
3468 0 : do iproc=0,nproc-1
3469 0 : do iw=1,nw
3470 0 : if (mod(iw-1,nproc)==iproc) then
3471 0 : procb(iw,:)=iproc
3472 0 : proct(iw,iproc)=1
3473 : endif
3474 : enddo
3475 : enddo
3476 : ! do iw=1,nw
3477 : ! write(6,*) "iw, iproc", iw, procb(iw,1)
3478 : ! enddo
3479 0 : nw_perproc=ratio+1
3480 : ! some procs will compute a number of frequency which is ratio and some
3481 : ! other will compute a number of frequency which is ratio+1.
3482 :
3483 : ! do iw=1,nw
3484 : ! write(6,*) " freq, proc"
3485 : ! write(6,*) iw, procb(iw,1)
3486 : ! enddo
3487 : endif
3488 :
3489 : ! do iw=1,nw
3490 : ! do iproc=0,nproc-1
3491 : ! write(6,*) " freq, procs,?"
3492 : ! write(6,*) iw,iproc,proct(iw,iproc)
3493 : ! enddo
3494 : ! enddo
3495 :
3496 :
3497 :
3498 0 : end subroutine distrib_paral
3499 : !!***
3500 :
3501 : !!****f* ABINIT/greendftcompute_green
3502 : !! NAME
3503 : !! greendftcompute_green
3504 : !!
3505 : !! FUNCTION
3506 : !! Compute levels for ctqmc
3507 : !!
3508 : !! COPYRIGHT
3509 : !! Copyright (C) 1999-2026 ABINIT group (BAmadon)
3510 : !! This file is distributed under the terms of the
3511 : !! GNU General Public License, see ~abinit/COPYING
3512 : !! or http://www.gnu.org/copyleft/gpl.txt .
3513 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
3514 : !!
3515 : !! INPUTS
3516 : !!
3517 : !! OUTPUT
3518 : !!
3519 : !! SOURCE
3520 :
3521 0 : subroutine greendftcompute_green(green,paw_dmft)
3522 :
3523 : !Arguments ------------------------------------
3524 : !scalars
3525 : type(paw_dmft_type), intent(in) :: paw_dmft
3526 : type(green_type),intent(inout) :: green
3527 :
3528 : !Local variables ------------------------------
3529 : ! scalars
3530 : integer :: ifreq,iband,ikpt,isppol
3531 : integer :: mbandc,natom,nspinor,nsppol,nkpt
3532 : character(len=500) :: message
3533 : ! arrays
3534 : !************************************************************************
3535 :
3536 0 : mbandc=paw_dmft%mbandc
3537 0 : nkpt=paw_dmft%nkpt
3538 0 : nsppol=paw_dmft%nsppol
3539 0 : nspinor=paw_dmft%nspinor
3540 0 : natom=paw_dmft%natom
3541 :
3542 : ! write(6,*) green%oper(1)%ks(1,1,1,1)
3543 0 : if(green%oper(1)%has_operks==0) then
3544 0 : ABI_ERROR("greendft%oper(1)%ks not allocated")
3545 : endif
3546 :
3547 : !======================================
3548 : !Get Green's function G=1/(iw_n+mu-e_nk)
3549 : !======================================
3550 0 : do ifreq=1,green%nw
3551 0 : do iband=1,mbandc
3552 0 : do ikpt=1,nkpt
3553 0 : do isppol=1,nsppol
3554 : ! write(6,*) green%oper(ifreq)%ks(isppol,ikpt,iband,iband)
3555 : ! write(6,*) ifreq,iband,ikpt,isppol
3556 : ! write(6,*) cmplx(0.d0,paw_dmft%omega_lo(ifreq),kind=dp)
3557 : ! write(6,*) paw_dmft%fermie
3558 : ! write(6,*) paw_dmft%eigen_dft(isppol,ikpt,iband)
3559 : green%oper(ifreq)%ks(iband,iband,ikpt,isppol)=&
3560 0 : cone/(cmplx(0.d0,paw_dmft%omega_lo(ifreq),kind=dp)+paw_dmft%fermie-paw_dmft%eigen_dft(iband,ikpt,isppol))
3561 : end do
3562 : end do
3563 : end do
3564 :
3565 :
3566 : !======================================================================
3567 : ! Compute local Green's function
3568 : !======================================================================
3569 0 : call downfold_oper(green%oper(ifreq),paw_dmft)
3570 :
3571 : !======================================================================
3572 : ! Symetrize
3573 : !======================================================================
3574 0 : call sym_matlu(green%oper(ifreq)%matlu,paw_dmft)
3575 : enddo
3576 0 : write(message,'(a,2x,a,f13.5)') ch10," == Print DFT Green's function for last frequency"
3577 0 : call wrtout(std_out,message,'COLL')
3578 0 : call print_matlu(green%oper(paw_dmft%dmft_nwlo)%matlu,natom,1)
3579 :
3580 0 : end subroutine greendftcompute_green
3581 : !!***
3582 :
3583 : !!****f* m_green/fermi_green
3584 : !! NAME
3585 : !! fermi_green
3586 : !!
3587 : !! FUNCTION
3588 : !! Compute Fermi level for DMFT or DFT.
3589 : !!
3590 : !! COPYRIGHT
3591 : !! Copyright (C) 2006-2026 ABINIT group (BAmadon)
3592 : !! This file is distributed under the terms of the
3593 : !! GNU General Public License, see ~abinit/COPYING
3594 : !! or http://www.gnu.org/copyleft/gpl.txt .
3595 : !!
3596 : !! INPUTS
3597 : !! green <type(green_type)>= green function data
3598 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
3599 : !! self <type(self_type)>= variables related to self-energy
3600 : !!
3601 : !! OUTPUT
3602 : !!
3603 : !! SOURCE
3604 :
3605 199 : subroutine fermi_green(green,paw_dmft,self)
3606 :
3607 : !Arguments ------------------------------------
3608 : type(green_type), intent(inout) :: green
3609 : type(paw_dmft_type), intent(inout) :: paw_dmft
3610 : !type(MPI_type), intent(in) :: mpi_enreg
3611 : type(self_type), intent(inout) :: self
3612 : !Local variables-------------------------------
3613 : integer :: ierr_hh,max_iter
3614 : real(dp) :: f_precision,fermi_old,x_precision
3615 : ! real(dp) :: hx
3616 : character(len=13) :: tag
3617 : character(len=500) :: message
3618 : !************************************************************************
3619 : !
3620 : ABI_NVTX_START_RANGE(NVTX_DMFT_FERMI_GREEN)
3621 199 : write(message,'(a,8x,a)') ch10," == Compute Fermi level"
3622 199 : call wrtout(std_out,message,'COLL')
3623 :
3624 : !=============
3625 : !headers
3626 : !=============
3627 199 : write(message,'(2a)') ch10," |--- Newton method to search for Fermi level ------------|"
3628 199 : call wrtout(std_out,message,'COLL')
3629 199 : write(tag,'(f13.6)') paw_dmft%fermie
3630 199 : write(message,'(3a)') ch10," |--- Initial value for Fermi level ",adjustl(tag)
3631 199 : call wrtout(std_out,message,'COLL')
3632 :
3633 : !========================================
3634 : !Define precision and nb of iterations
3635 : !=========================================
3636 199 : fermi_old = paw_dmft%fermie
3637 199 : ierr_hh = 0
3638 199 : f_precision = paw_dmft%dmft_charge_prec
3639 : !f_precision=0.01
3640 199 : x_precision = tol5
3641 : !if(option==1) then
3642 : !f_precision=(erreursurlacharge)/100d0
3643 : !else
3644 : !f_precision=tol11
3645 : !endif
3646 : !max_iter=1 ! for tests only
3647 : !write(6,*) "for tests max_iter=1"
3648 199 : max_iter = 50
3649 : !opt_noninter = 4
3650 :
3651 : ! Precompute some useful quantities so that the update of the moments can be done in constant time at each iteration
3652 199 : green%trace_fermie(1) = cmplx(paw_dmft%nsppol*paw_dmft%mbandc,zero,kind=dp)
3653 199 : if (paw_dmft%nsppol == 1 .and. paw_dmft%nspinor == 1) green%trace_fermie(1) = &
3654 108 : & green%trace_fermie(1) * two
3655 199 : if (green%has_moments == 1) then
3656 0 : call compute_trace_moments_ks(green,self,paw_dmft)
3657 : end if
3658 :
3659 : !=====================
3660 : !Call newton method
3661 : !=====================
3662 199 : write(message,'(a,4x,a,e13.6)') ch10," Precision required :",f_precision
3663 199 : call wrtout(std_out,message,'COLL')
3664 199 : if (f_precision < ten) then
3665 199 : call newton(green,self,paw_dmft,paw_dmft%fermie,x_precision,max_iter,f_precision,ierr_hh)
3666 : end if
3667 :
3668 : !===========================
3669 : !Deals with errors signals
3670 : !===========================
3671 199 : if (ierr_hh == -314) then
3672 0 : write(message,'(a)') "Warning, check Fermi level"
3673 0 : call wrtout(std_out,message,'COLL')
3674 : ! call abi_abort('COLL')
3675 0 : write(message,'(2a,f13.6)') ch10," |--- Final value for Fermi level (check)",paw_dmft%fermie
3676 0 : call wrtout(std_out,message,'COLL')
3677 199 : else if (ierr_hh == -123) then
3678 0 : write(message,'(a,f13.6)') " Fermi level is set to",fermi_old
3679 0 : paw_dmft%fermie = fermi_old
3680 0 : call wrtout(std_out,message,'COLL')
3681 :
3682 : ! =====================================
3683 : ! If fermi level search was successful
3684 : ! =====================================
3685 : else
3686 199 : write(tag,'(e13.6)') x_precision
3687 199 : write(message,'(a,4x,2a)') ch10," Precision achieved on Fermi Level : ",adjustl(tag)
3688 199 : call wrtout(std_out,message,'COLL')
3689 199 : write(tag,'(e13.6)') f_precision
3690 199 : write(message,'(4x,2a)') " Precision achieved on number of electrons : ",adjustl(tag)
3691 199 : call wrtout(std_out,message,'COLL')
3692 199 : write(tag,'(f13.6)') paw_dmft%fermie
3693 199 : write(message,'(3a)') ch10," |--- Final value for Fermi level ",adjustl(tag)
3694 199 : call wrtout(std_out,message,'COLL')
3695 : end if ! ierr_hh
3696 :
3697 : !========================================================
3698 : ! Check convergence of fermi level during DMFT iterations
3699 : !========================================================
3700 199 : if (paw_dmft%idmftloop >= 2) then
3701 110 : if (abs(paw_dmft%fermie-fermi_old) <= paw_dmft%dmft_fermi_prec) then
3702 : ! write(message,'(a,8x,a,e9.2,a,8x,a,e12.5)') ch10,"|fermie(n)-fermie(n-1)|=<",paw_dmft%dmft_fermi_prec,ch10,&
3703 1 : write(message,'(a,8x,a,e9.2,a,e9.2,a,8x,a,e12.5)') ch10,"|fermie(n)-fermie(n-1)|=",&
3704 1 : & abs(paw_dmft%fermie-fermi_old),"<",paw_dmft%dmft_fermi_prec,ch10,&
3705 2 : & "=> DMFT Loop: Fermi level is converged to:",paw_dmft%fermie
3706 1 : call wrtout(std_out,message,'COLL')
3707 1 : green%ifermie_cv = 1
3708 : else
3709 109 : write(tag,'(f12.5)') paw_dmft%fermie
3710 109 : write(message,'(a,8x,2a)') ch10,"DMFT Loop: Fermi level is not converged: ",adjustl(tag)
3711 109 : call wrtout(std_out,message,'COLL')
3712 109 : green%ifermie_cv = 0
3713 : end if ! convergence of Fermi level
3714 : end if ! idmftloop>=2
3715 199 : write(message,'(2a)') ch10, " |---------------------------------------------------|"
3716 199 : call wrtout(std_out,message,'COLL')
3717 : !
3718 : !==========================================================
3719 : !Recompute full green function including non diag elements
3720 : !==========================================================
3721 : !call compute_green(green,paw_dmft,0,self,opt_self=1,opt_nonxsum=1)
3722 : !call integrate_green(green,paw_dmft,prtopt=0,opt_ksloc=3) !,opt_nonxsum=1)
3723 :
3724 : ABI_NVTX_END_RANGE()
3725 : !return
3726 199 : end subroutine fermi_green
3727 : !!***
3728 :
3729 : !!****f* m_green/newton
3730 : !! NAME
3731 : !! newton
3732 : !!
3733 : !! FUNCTION
3734 : !! Compute root of a function with newton methods (newton/halley)
3735 : !!
3736 : !! COPYRIGHT
3737 : !! Copyright (C) 2006-2026 ABINIT group (BAmadon)
3738 : !! This file is distributed under the terms of the
3739 : !! GNU General Public License, see ~abinit/COPYING
3740 : !! or http://www.gnu.org/copyleft/gpl.txt .
3741 : !!
3742 : !! INPUTS
3743 : !! green <type(green_type)>= green function data
3744 : !! self <type(self_type)>= variables related to self-energy
3745 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
3746 : !! x_input : input value for x
3747 : !! max_iter : maximum number of iterations
3748 : !! f_precision : required precision on function F
3749 : !! opt_algo : 1 for Newton, 2 for Halley
3750 : !!
3751 : !! OUTPUT
3752 : !! x_precision : output precision on x
3753 : !! ierr_hh : different from zero if an error occurs
3754 : !!
3755 : !! SOURCE
3756 :
3757 199 : subroutine newton(green,self,paw_dmft,x_input,x_precision,max_iter,&
3758 : & f_precision,ierr_hh,opt_algo)
3759 :
3760 : !Arguments ------------------------------------
3761 : type(green_type), intent(inout) :: green
3762 : type(self_type), intent(inout) :: self
3763 : !type(MPI_type), intent(in) :: mpi_enreg
3764 : type(paw_dmft_type), intent(inout) :: paw_dmft
3765 : integer, intent(in) :: max_iter
3766 : integer, intent(out) :: ierr_hh
3767 : real(dp), intent(in) :: f_precision
3768 : real(dp), intent(inout) :: x_input,x_precision
3769 : integer, optional, intent(in) :: opt_algo
3770 : !Local variables-------------------------------
3771 : integer :: iter,option
3772 : logical :: dmft_optim,l_minus,l_plus
3773 : real(dp) :: Fx,Fxdouble,Fxoptimum,Fxprime,nb_elec_x,step
3774 : real(dp) :: x_minus,x_optimum,x_plus,xold
3775 : character(len=500) :: message
3776 : ! *********************************************************************
3777 :
3778 199 : x_minus = - dble(10)
3779 199 : x_plus = - dble(11)
3780 199 : xold = - dble(12)
3781 :
3782 199 : ierr_hh = 0
3783 199 : option = 1
3784 199 : if (present(opt_algo)) option = opt_algo
3785 199 : step = paw_dmft%dmft_fermi_step
3786 :
3787 : !write(std_out,*) "ldaprint",opt_noninter
3788 :
3789 : !--- Start of iterations
3790 199 : write(message,'(a,3a)') " Fermi level Charge Difference"
3791 199 : call wrtout(std_out,message,'COLL')
3792 : !do iter=1, 40
3793 : !x_input=float(iter)/100_dp
3794 : !call function_and_deriv(cryst_struc,f_precision,green,iter,mpi_enreg,paw_dmft,pawang,self,&
3795 : !& x_input,x_before,x_precision,Fx,Fxprime,Fxdouble,opt_noninter,option)
3796 : !write(std_out,*) x_input,Fx
3797 : !enddo
3798 : !call abi_abort('COLL')
3799 :
3800 199 : l_minus = .false.
3801 199 : l_plus = .false.
3802 199 : Fxoptimum = one
3803 199 : x_optimum = zero
3804 :
3805 199 : dmft_optim = (paw_dmft%dmft_solv == 6 .or. paw_dmft%dmft_solv == 7)
3806 :
3807 199 : if (dmft_optim) xold = x_input
3808 :
3809 : !========================================
3810 : ! Start iteration to find fermi level
3811 : !========================================
3812 2394 : do iter=1,max_iter
3813 :
3814 : ! ========================================
3815 : ! If zero is located between two values: apply newton method or dichotomy
3816 : ! ========================================
3817 2394 : if ((l_minus .and. l_plus) .or. dmft_optim) then
3818 :
3819 : ! ==============================================
3820 : ! Compute the function and derivatives for newton
3821 : ! ==============================================
3822 : call function_and_deriv(green,self,paw_dmft,x_input,x_precision, &
3823 1975 : & f_precision,Fx,Fxprime,Fxdouble,option)
3824 :
3825 : ! Apply stop criterion on Fx
3826 1975 : if (abs(Fx) < f_precision) then
3827 : ! write(message,'(a,2f12.6)') "Fx,f_precision",Fx,f_precision
3828 : ! call wrtout(std_out,message,'COLL')
3829 199 : x_precision = x_input - xold
3830 199 : return
3831 : end if ! abs(Fx)<f_precision
3832 1776 : if (iter == max_iter) then
3833 0 : write(message,'(a,2f12.6)') " Fermi level could not be found"
3834 0 : call wrtout(std_out,message,'COLL')
3835 0 : x_input = x_optimum
3836 0 : ierr_hh = -123
3837 0 : return
3838 : end if ! iter=max_iter
3839 :
3840 : ! Cannot divide by Fxprime if too small
3841 1776 : if (abs(Fxprime) <= tol15) then
3842 0 : ierr_hh = -314
3843 0 : write(message,'(a,f12.7)') "Fxprime=",Fxprime
3844 0 : call wrtout(std_out,message,'COLL')
3845 0 : return
3846 : end if ! abs(Fxprime)<=tol15
3847 :
3848 1776 : x_precision = x_input - xold
3849 :
3850 : ! ==============================================
3851 : ! Newton/Halley's formula for next iteration
3852 : ! ==============================================
3853 1776 : xold = x_input
3854 1776 : if (option == 1) then
3855 1776 : if (dmft_optim) then
3856 0 : x_input = x_input - sign(one,Fx)*merge(step,min(step,abs(Fx/Fxprime)),Fxprime<0)
3857 : else
3858 1776 : x_input = x_input - Fx/Fxprime
3859 : end if ! dmft_optim
3860 : end if ! option=1
3861 1776 : if (option == 2) x_input = x_input - two*Fx*Fxprime/(two*(Fxprime**2)-Fx*Fxdouble)
3862 :
3863 : ! ==============================================
3864 : ! If newton does not work well, use dichotomy.
3865 : ! ==============================================
3866 :
3867 1776 : if (dmft_optim) then
3868 0 : if (Fx < 0) then
3869 0 : l_minus = .true.
3870 0 : x_minus = xold
3871 : end if
3872 0 : if (Fx > 0) then
3873 0 : l_plus = .true.
3874 0 : x_plus = xold
3875 : end if
3876 : end if ! dmft_optim
3877 :
3878 1776 : if ((x_input < x_minus .or. x_input > x_plus) .and. (l_minus .and. l_plus)) then
3879 :
3880 33 : if (.not. dmft_optim) then
3881 33 : call compute_nb_elec(green,self,paw_dmft,Fx,nb_elec_x,xold)
3882 : end if
3883 :
3884 33 : write(message,'(a,3f12.6)') " ---",x_input,Fx+paw_dmft%nelectval,Fx
3885 33 : call wrtout(std_out,message,'COLL')
3886 33 : if (.not. dmft_optim) then
3887 33 : if (Fx > 0) then
3888 : x_plus = xold
3889 17 : else if (Fx < 0) then
3890 33 : x_minus = xold
3891 : end if ! Fx>0
3892 : end if
3893 33 : x_input = (x_plus+x_minus) * half
3894 :
3895 : end if ! x_input<x_minus or x_input>x_plus
3896 : ! write(std_out,'(a,2f12.6)') " Q(xold) and dQ/dx=",Fx,Fxprime
3897 : ! write(std_out,'(a,f12.6)') " => new Fermi level",x_input
3898 : ! ========================================
3899 : ! Locate zero between two values
3900 : ! ========================================
3901 : else
3902 419 : call compute_nb_elec(green,self,paw_dmft,Fx,nb_elec_x,x_input)
3903 419 : write(message,'(a,3f12.6)') " --",x_input,nb_elec_x,Fx
3904 : ! Possible improvement for large systems, removed temporarely for
3905 : ! automatic tests: more study is necessary: might worsen the convergency
3906 : ! if(iter==1) then
3907 : ! f_precision=max(abs(Fx/50),f_precision)
3908 : ! write(message,'(a,4x,a,e12.6)') ch10," Precision required changed to:",f_precision
3909 : ! call wrtout(std_out,message,'COLL')
3910 : ! endif
3911 419 : call wrtout(std_out,message,'COLL')
3912 419 : if (Fx > 0) then
3913 217 : l_plus = .true.
3914 217 : x_plus = x_input
3915 217 : x_input = x_input - step
3916 202 : else if (Fx < 0) then
3917 202 : l_minus = .true.
3918 202 : x_minus = x_input
3919 202 : x_input = x_input + step
3920 : end if ! Fx>0
3921 :
3922 : end if ! l_minus and l_plus
3923 :
3924 2195 : if (abs(Fx) < abs(Fxoptimum) .or. (iter == 1 .and. dmft_optim)) then
3925 1169 : Fxoptimum = Fx
3926 1169 : x_optimum = merge(xold,x_input,dmft_optim)
3927 : end if ! abs(Fx)<abs(Fxoptimum)
3928 :
3929 :
3930 :
3931 : ! if(myid==master) then
3932 : ! write(std_out,'(a,i4,3f12.6)') "i,xnew,F,Fprime",i,x_input,Fx,Fxprime
3933 : ! endif
3934 :
3935 :
3936 : ! ! Apply stop criterion on x
3937 : ! if(abs(x_input-xold) .le. x_input*x_precision) then
3938 : ! ! write(std_out,'(a,4f12.6)') "x_input-xold, x_precision*x_input "&
3939 : ! !& ,x_input-xold,x_precision*x_input,x_precision
3940 : ! f_precision=Fx
3941 : ! return
3942 : ! endif
3943 :
3944 : end do ! iter
3945 : !--- End of iterations
3946 :
3947 :
3948 0 : ierr_hh = 1
3949 : !return
3950 :
3951 : CONTAINS !========================================================================================
3952 : !-----------------------------------------------------------------------
3953 : !!***
3954 :
3955 : !!****f* newton/function_and_deriv
3956 : !! NAME
3957 : !! function_and_deriv
3958 : !!
3959 : !! FUNCTION
3960 : !! Compute value of a function and its numerical derivatives
3961 : !!
3962 : !! INPUTS
3963 : !! green <type(green_type)>= green function data
3964 : !! self <type(self_type)>= variables related to self-energy
3965 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
3966 : !! x_input = input value for x
3967 : !! option = 1 to only compute the first derivative
3968 : !! = 2 to compute the two first derivatives.
3969 : !!
3970 : !! OUTPUTS
3971 : !! Fx : Value of F(x)
3972 : !! Fxprime : Value of F'(x)
3973 : !! Fxdouble : Value of F''(x)
3974 : !!
3975 : !! SOURCE
3976 :
3977 1975 : subroutine function_and_deriv(green,self,paw_dmft,x_input,x_precision, &
3978 : & f_precision,Fx,Fxprime,Fxdouble,option)
3979 :
3980 : !Arguments ------------------------------------
3981 : type(green_type), intent(inout) :: green
3982 : type(self_type), intent(inout) :: self
3983 : !type(MPI_type), intent(in) :: mpi_enreg
3984 : type(paw_dmft_type), intent(inout) :: paw_dmft
3985 : integer, intent(in) :: option
3986 : real(dp), intent(in) :: f_precision,x_input,x_precision
3987 : real(dp), intent(out) :: Fx,Fxprime,Fxdouble
3988 : !Local variables-------------------------------
3989 : logical :: dmft_optim
3990 : real(dp) :: deltax,Fxminus,Fxplus,nb_elec_x,xminus,x0,xplus
3991 : character(len=500) :: message
3992 : ! *********************************************************************
3993 :
3994 1975 : dmft_optim = (paw_dmft%dmft_solv == 6 .or. paw_dmft%dmft_solv == 7)
3995 :
3996 1975 : if (.not. dmft_optim) then
3997 :
3998 : ! Choose deltax: for numeric evaluation of derivative
3999 : !if(iter==1) then
4000 : ! deltax=0.02
4001 : !end if
4002 : ! deltax=max((x_input-x_old)/10.d0,min(0.00001_dp,x_precision/100_dp))
4003 1975 : deltax = min(tol5,x_precision/dble(100)) ! small but efficient
4004 : ! endif
4005 : ! write(std_out,*) "iter,x_input,deltax",iter,x_input,deltax
4006 1975 : x0 = x_input
4007 1975 : xminus = x0 - deltax
4008 1975 : xplus = x0 + deltax
4009 :
4010 1975 : call compute_nb_elec(green,self,paw_dmft,Fx,nb_elec_x,x0)
4011 :
4012 1975 : write(message,'(a,3f12.6)') " - ",x0,nb_elec_x,Fx
4013 1975 : call wrtout(std_out,message,'COLL')
4014 :
4015 : ! write(std_out,*) "Fx", Fx
4016 1975 : if (abs(Fx) < f_precision) return
4017 :
4018 1776 : call compute_nb_elec(green,self,paw_dmft,Fxplus,nb_elec_x,xplus)
4019 :
4020 1776 : write(message,'(a,3f12.6)') " - ",xplus,nb_elec_x,Fxplus
4021 1776 : call wrtout(std_out,message,'COLL')
4022 :
4023 1776 : if (option == 2) then
4024 0 : call compute_nb_elec(green,self,paw_dmft,Fxminus,nb_elec_x,xminus)
4025 :
4026 0 : write(message,'(a,3f12.6)') " - ",xminus,nb_elec_x,Fxminus
4027 0 : call wrtout(std_out,message,'COLL')
4028 : end if ! option=2
4029 :
4030 1776 : if (option == 1) then
4031 1776 : Fxprime = (Fxplus-Fx) / deltax
4032 0 : else if (option == 2) then
4033 0 : Fxprime = (Fxplus-Fxminus) / (two*deltax)
4034 0 : Fxdouble = (Fxplus+Fxminus-two*Fx) / (deltax**2)
4035 : end if ! option
4036 : ! write(std_out,*) "after computation of Fxprime",myid
4037 :
4038 : else
4039 :
4040 0 : call compute_nb_elec(green,self,paw_dmft,Fx,nb_elec_x,x_input,Fxprime=Fxprime)
4041 :
4042 0 : write(message,'(a,3f12.6)') " - ",x_input,nb_elec_x,Fx
4043 0 : call wrtout(std_out,message,'COLL')
4044 :
4045 : end if ! dmft_optim
4046 :
4047 1776 : if (Fxprime < zero) then
4048 0 : write(message,'(a,f12.6)') " Warning: slope of charge versus fermi level is negative !",Fxprime
4049 0 : call wrtout(std_out,message,'COLL')
4050 : end if
4051 :
4052 : end subroutine function_and_deriv
4053 : !!***
4054 :
4055 : !!****f* newton/compute_nb_elec
4056 : !! NAME
4057 : !! compute_nb_elec
4058 : !!
4059 : !! FUNCTION
4060 : !! Compute nb of electrons as a function of Fermi level
4061 : !!
4062 : !! INPUTS
4063 : !! fermie : input of energy
4064 : !! opt_noninter
4065 : !!
4066 : !! OUTPUTS
4067 : !! Fx : Value of F(x).
4068 : !! nb_elec_x : Number of electrons for the value of x
4069 : !!
4070 : !! SOURCE
4071 :
4072 4203 : subroutine compute_nb_elec(green,self,paw_dmft,Fx,nb_elec_x,fermie,Fxprime)
4073 :
4074 : !Arguments ------------------------------------
4075 : type(green_type), intent(inout) :: green
4076 : type(self_type), intent(inout) :: self
4077 : !type(MPI_type), intent(in) :: mpi_enreg
4078 : type(paw_dmft_type), intent(inout) :: paw_dmft
4079 : !integer,intent(in) :: opt_noninter
4080 : real(dp), intent(in) :: fermie
4081 : real(dp), intent(out) :: Fx,nb_elec_x
4082 : real(dp), optional, intent(out) :: Fxprime
4083 : !Local variables-------------------------------
4084 : integer :: band_index,i,ib,ierr,ifreq,ikpt,isppol,mbandc
4085 : integer :: mkmem,nband_k,nkpt,nmoments,nspinor,nsppol,shift
4086 : logical :: dmft_optim
4087 : real(dp) :: correction,correction_prime,eig
4088 : real(dp) :: fac,occ_prime,temp,wtk
4089 : complex(dp) :: omega
4090 4203 : type(oper_type) :: oper_tmp
4091 4203 : complex(dp), allocatable :: omega_fac(:),trace_moments(:),trace_moments_prime(:)
4092 : ! *********************************************************************
4093 :
4094 : ABI_NVTX_START_RANGE(NVTX_DMFT_COMPUTE_NB_ELEC)
4095 4203 : dmft_optim = (paw_dmft%dmft_solv == 6 .or. paw_dmft%dmft_solv == 7)
4096 :
4097 4203 : if (.not. dmft_optim) then
4098 :
4099 4203 : paw_dmft%fermie = fermie
4100 4203 : call compute_green(green,paw_dmft,0,self,opt_self=1,opt_nonxsum=1,opt_nonxsum2=1)
4101 4203 : call integrate_green(green,paw_dmft,0,opt_ksloc=-1) !,opt_nonxsum=1)
4102 : ! opt_ksloc=-1, compute total charge
4103 :
4104 : else
4105 :
4106 0 : green%charge_ks = zero
4107 :
4108 0 : mkmem = green%distrib%nkpt_mem(green%distrib%me_kpt+1)
4109 0 : shift = green%distrib%shiftk
4110 0 : nmoments = merge(green%nmoments,1,green%has_moments==1)
4111 :
4112 0 : mbandc = paw_dmft%mbandc
4113 0 : nkpt = paw_dmft%nkpt
4114 0 : nspinor = paw_dmft%nspinor
4115 0 : nsppol = paw_dmft%nsppol
4116 0 : temp = paw_dmft%temp
4117 :
4118 0 : ABI_MALLOC(trace_moments,(nmoments))
4119 0 : ABI_MALLOC(trace_moments_prime,(nmoments))
4120 :
4121 0 : trace_moments(1) = green%trace_fermie(1)
4122 0 : trace_moments_prime(:) = czero
4123 0 : Fxprime = zero
4124 :
4125 0 : if (green%has_moments == 1) then
4126 0 : call compute_trace_moments(fermie,green%trace_fermie(:),trace_moments(:),trace_moments_prime(:))
4127 : end if
4128 :
4129 0 : call init_oper(paw_dmft,oper_tmp,nkpt=mkmem,shiftk=shift)
4130 :
4131 0 : do ifreq=1,green%nw
4132 0 : if (green%distrib%proct(ifreq) /= green%distrib%me_freq) cycle
4133 0 : omega = cmplx(zero,green%omega(ifreq),kind=dp)
4134 0 : fac = two * paw_dmft%wgt_wlo(ifreq) * merge(two*temp,temp,nsppol==1.and.nspinor==1)
4135 :
4136 0 : call add_matlu(self%hdc%matlu(:),self%oper(ifreq)%matlu(:),oper_tmp%matlu(:),paw_dmft%natom,-1)
4137 0 : call upfold_oper(oper_tmp,paw_dmft)
4138 :
4139 0 : do isppol=1,nsppol
4140 0 : do ikpt=1,mkmem
4141 0 : wtk = paw_dmft%wtk(ikpt+shift)
4142 0 : do ib=1,mbandc
4143 : oper_tmp%ks(ib,ib,ikpt,isppol) = oper_tmp%ks(ib,ib,ikpt,isppol) + omega + &
4144 0 : & fermie - paw_dmft%eigen_dft(ib,ikpt+shift,isppol)
4145 : end do ! ib
4146 0 : call xginv(oper_tmp%ks(:,:,ikpt,isppol),mbandc)
4147 :
4148 0 : Fxprime = Fxprime - dble(sum(oper_tmp%ks(:,:,ikpt,isppol)*transpose(oper_tmp%ks(:,:,ikpt,isppol))))*wtk*fac
4149 0 : do ib=1,mbandc
4150 0 : green%charge_ks = green%charge_ks + dble(oper_tmp%ks(ib,ib,ikpt,isppol))*wtk*fac
4151 : end do ! ib
4152 : end do ! ikpt
4153 : end do ! isppol
4154 :
4155 : end do ! ifreq
4156 :
4157 0 : call xmpi_sum(green%charge_ks,paw_dmft%spacecomm,ierr)
4158 0 : call xmpi_sum(Fxprime,paw_dmft%spacecomm,ierr)
4159 :
4160 0 : call destroy_oper(oper_tmp)
4161 :
4162 0 : ABI_MALLOC(omega_fac,(nmoments))
4163 :
4164 0 : do i=1,nmoments
4165 0 : omega_fac(i) = czero
4166 0 : do ifreq=green%nw,1,-1 ! NEVER change the summation order and DON'T use the intrinsic SUM
4167 0 : omega_fac(i) = omega_fac(i) + paw_dmft%wgt_wlo(ifreq) / (paw_dmft%omega_lo(ifreq))**i
4168 : end do
4169 0 : omega_fac(i) = - two * temp * omega_fac(i) / (j_dpc)**i
4170 0 : if (i == 1) omega_fac(i) = omega_fac(i) + half
4171 0 : if (i == 2) omega_fac(i) = omega_fac(i) - cone/(four*temp)
4172 0 : if (i == 4) omega_fac(i) = omega_fac(i) + cone/(dble(48)*(temp**3))
4173 : end do ! i
4174 :
4175 : ! Do not use DOT_PRODUCT
4176 0 : green%charge_ks = green%charge_ks + dble(sum(trace_moments(1:nmoments)*omega_fac(1:nmoments)))
4177 : ! Do not use DOT_PRODUCT
4178 0 : Fxprime = Fxprime + dble(sum(trace_moments_prime(1:nmoments)*omega_fac(1:nmoments)))
4179 :
4180 0 : ABI_FREE(trace_moments)
4181 0 : ABI_FREE(trace_moments_prime)
4182 0 : ABI_FREE(omega_fac)
4183 :
4184 0 : if (paw_dmft%dmft_solv == 6 .or. paw_dmft%dmft_solv == 7) then
4185 : band_index = 0
4186 : correction = zero
4187 : correction_prime = zero
4188 0 : do isppol=1,nsppol
4189 0 : do ikpt=1,nkpt
4190 0 : nband_k = paw_dmft%nband(ikpt+(isppol-1)*nkpt)
4191 0 : wtk = paw_dmft%wtk(ikpt)
4192 0 : do ib=1,nband_k
4193 0 : if (paw_dmft%band_in(ib)) cycle
4194 0 : eig = paw_dmft%eigen(ib+band_index)
4195 0 : correction = correction + occup_fd(eig,fermie,temp)*wtk
4196 0 : if ((eig-fermie) > zero) then
4197 0 : occ_prime = exp(-(eig-fermie)/temp) / (one+exp(-(eig-fermie)/temp))**2 / temp
4198 : else
4199 0 : occ_prime = exp((eig-fermie)/temp) / (one+exp((eig-fermie)/temp))**2 / temp
4200 : end if
4201 0 : correction_prime = correction_prime + occ_prime*wtk
4202 : end do ! ib
4203 0 : band_index = band_index + nband_k
4204 : end do ! ikpt
4205 : end do ! isppol
4206 0 : if (nsppol == 1 .and. nspinor == 1) correction = correction * two
4207 0 : green%charge_ks = green%charge_ks + correction
4208 0 : if (nsppol == 1 .and. nspinor == 1) correction_prime = correction_prime * two
4209 0 : Fxprime = Fxprime + correction_prime
4210 : end if ! use_all_bands
4211 :
4212 : end if ! dmft_optim
4213 :
4214 4203 : nb_elec_x = green%charge_ks
4215 4203 : Fx = green%charge_ks - paw_dmft%nelectval
4216 :
4217 : ABI_NVTX_END_RANGE()
4218 4203 : end subroutine compute_nb_elec
4219 : !!***
4220 :
4221 0 : subroutine compute_trace_moments(fermie,trace_fermie,trace_moments,trace_moments_prime)
4222 :
4223 : !Arguments ------------------------------------
4224 : real(dp), intent(in) :: fermie
4225 : complex(dp), intent(in) :: trace_fermie(:)
4226 : complex(dp), intent(inout) :: trace_moments(:),trace_moments_prime(:)
4227 : !Local variables-------------------------------
4228 : real(dp) :: fermie2,fermie3,fermie4
4229 : ! *********************************************
4230 :
4231 0 : fermie2 = fermie * fermie
4232 0 : fermie3 = fermie2 * fermie
4233 0 : fermie4 = fermie3 * fermie
4234 :
4235 0 : trace_moments(2) = trace_fermie(2) - fermie*trace_fermie(1)
4236 : trace_moments(3) = trace_fermie(3) + trace_fermie(4) - two*fermie*trace_fermie(2) + &
4237 0 : & fermie2*trace_fermie(1)
4238 : trace_moments(4) = trace_fermie(5) + two*(trace_fermie(6)-fermie*trace_fermie(3)) + &
4239 : & trace_fermie(7) - three*fermie*trace_fermie(4) + &
4240 0 : & three*fermie2*trace_fermie(2) - fermie3*trace_fermie(1)
4241 : trace_moments(5) = trace_fermie(8) + two*(trace_fermie(9)-fermie*trace_fermie(5)) + &
4242 : & trace_fermie(10) + three*(trace_fermie(11)-two*fermie*trace_fermie(6)+ &
4243 : & fermie2*trace_fermie(3)) + trace_fermie(12) - four*fermie*trace_fermie(7) + &
4244 : & six*fermie2*trace_fermie(4) - four*fermie3*trace_fermie(2) + &
4245 0 : & fermie4*trace_fermie(1)
4246 :
4247 0 : trace_moments_prime(2) = - trace_fermie(1)
4248 0 : trace_moments_prime(3) = two * (fermie*trace_fermie(1)-trace_fermie(2))
4249 : trace_moments_prime(4) = - two*trace_fermie(3) - three*trace_fermie(4) + &
4250 0 : & six*fermie*trace_fermie(2) - three*fermie2*trace_fermie(1)
4251 : trace_moments_prime(5) = - two*trace_fermie(5) + three*(-two*trace_fermie(6)+ &
4252 : & two*fermie*trace_fermie(3)) - four*trace_fermie(7) + &
4253 : & dble(12)*fermie*trace_fermie(4) - dble(12)*fermie2*trace_fermie(2) + &
4254 0 : & four*fermie3*trace_fermie(1)
4255 :
4256 0 : end subroutine compute_trace_moments
4257 : !!***
4258 :
4259 : end subroutine newton
4260 : !!***
4261 :
4262 : !!****f* m_green/local_ks_green
4263 : !! NAME
4264 : !! local_ks_green
4265 : !!
4266 : !! FUNCTION
4267 : !! Compute the sum over k-point of ks green function.
4268 : !! do the fourier transformation and print it
4269 : !!
4270 : !! COPYRIGHT
4271 : !! Copyright (C) 1999-2026 ABINIT group (BAmadon)
4272 : !! This file is distributed under the terms of the
4273 : !! GNU General Public License, see ~abinit/COPYING
4274 : !! or http://www.gnu.org/copyleft/gpl.txt .
4275 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
4276 : !!
4277 : !! INPUTS
4278 : !! cryst_struc
4279 : !! istep = step of iteration for DFT.
4280 : !! dft_occup
4281 : !! mpi_enreg=information about MPI parallelization
4282 : !! paw_dmft = data for self-consistent DFT+DMFT calculations.
4283 : !!
4284 : !! OUTPUT
4285 : !! paw_dmft = data for self-consistent DFT+DMFT calculations.
4286 : !!
4287 : !! SOURCE
4288 :
4289 0 : subroutine local_ks_green(green,paw_dmft,prtopt)
4290 :
4291 : !Arguments ------------------------------------
4292 : !scalars
4293 : type(green_type), intent(in) :: green
4294 : type(paw_dmft_type), intent(in) :: paw_dmft
4295 : integer, intent(in) :: prtopt
4296 :
4297 : !Local variables ------------------------------
4298 : character(len=500) :: message
4299 : integer :: iband,ifreq,ikpt,isppol,itau,lsub,ltau,mbandc,nkpt,nsppol
4300 : character(len=1) :: tag_is
4301 : character(len=fnlen) :: tmpfil
4302 0 : integer,allocatable :: unitgreenlocks_arr(:)
4303 : real(dp) :: beta
4304 0 : real(dp), allocatable :: tau(:)
4305 0 : complex(dp), allocatable :: loc_ks(:,:,:)
4306 0 : complex(dp), allocatable :: loc_ks_tau(:,:,:),fw(:),ft(:)
4307 : !************************************************************************
4308 0 : mbandc=paw_dmft%mbandc
4309 0 : nkpt=paw_dmft%nkpt
4310 0 : nsppol=paw_dmft%nsppol
4311 0 : ltau=128
4312 0 : ABI_MALLOC(tau,(ltau))
4313 0 : do itau=1,ltau
4314 0 : tau(itau)=float(itau-1)/float(ltau)/paw_dmft%temp
4315 : end do
4316 0 : beta=one/paw_dmft%temp
4317 :
4318 : !Only imaginary frequencies here
4319 0 : if(green%w_type=="real") then
4320 0 : message = ' compute_energy not implemented for real frequency'
4321 0 : ABI_BUG(message)
4322 : end if
4323 :
4324 : !=========================================
4325 : !Compute local band ks green function
4326 : ! should be computed in compute_green: it would be less costly in memory.
4327 : !=========================================
4328 0 : ABI_MALLOC(loc_ks,(nsppol,mbandc,paw_dmft%dmft_nwlo))
4329 0 : if(green%oper(1)%has_operks==1) then
4330 0 : loc_ks(:,:,:)=czero
4331 0 : do isppol=1,nsppol
4332 0 : do iband=1,mbandc
4333 0 : do ifreq=1,paw_dmft%dmft_nwlo
4334 0 : do ikpt=1,nkpt
4335 : loc_ks(isppol,iband,ifreq)=loc_ks(isppol,iband,ifreq)+ &
4336 0 : & green%oper(ifreq)%ks(iband,iband,ikpt,isppol)*paw_dmft%wtk(ikpt)
4337 : end do
4338 : end do
4339 : end do
4340 : end do
4341 : else
4342 0 : message = ' green fct is not computed in ks space'
4343 0 : ABI_BUG(message)
4344 : end if
4345 :
4346 : !=========================================
4347 : !Compute fourier transformation
4348 : !=========================================
4349 :
4350 0 : ABI_MALLOC(loc_ks_tau,(nsppol,mbandc,ltau))
4351 0 : ABI_MALLOC(fw,(paw_dmft%dmft_nwlo))
4352 0 : ABI_MALLOC(ft,(ltau))
4353 0 : loc_ks_tau(:,:,:)=czero
4354 0 : do isppol=1,nsppol
4355 0 : do iband=1,mbandc
4356 0 : do ifreq=1,paw_dmft%dmft_nwlo
4357 0 : fw(ifreq)=loc_ks(isppol,iband,ifreq)
4358 : end do
4359 0 : call fourier_fct(fw,ft,.true.,ltau,-1,paw_dmft) ! inverse fourier
4360 0 : do itau=1,ltau
4361 0 : loc_ks_tau(isppol,iband,itau)=ft(itau)
4362 : end do
4363 : end do
4364 : end do
4365 0 : ABI_FREE(fw)
4366 0 : ABI_FREE(ft)
4367 0 : do isppol=1,nsppol
4368 0 : do iband=1,mbandc
4369 0 : do itau=1,ltau
4370 0 : loc_ks_tau(isppol,iband,itau)=(loc_ks_tau(isppol,iband,itau)+conjg(loc_ks_tau(isppol,iband,itau)))/two
4371 : end do
4372 : end do
4373 : end do
4374 :
4375 : !=========================================
4376 : !Print out ksloc green function
4377 : !=========================================
4378 0 : if(abs(prtopt)==1) then
4379 0 : ABI_MALLOC(unitgreenlocks_arr,(nsppol))
4380 0 : do isppol=1,nsppol
4381 0 : write(tag_is,'(i1)')isppol
4382 0 : tmpfil = trim(paw_dmft%filapp)//'Gtau_locks_isppol'//tag_is
4383 0 : write(message,'(3a)') ch10," == Print green function on file ",trim(tmpfil)
4384 0 : call wrtout(std_out,message,'COLL')
4385 0 : unitgreenlocks_arr(isppol)=500+isppol-1
4386 0 : open (unit=unitgreenlocks_arr(isppol),file=trim(tmpfil),status='unknown',form='formatted')
4387 0 : rewind(unitgreenlocks_arr(isppol))
4388 0 : write(message,'(a,a,a,i4)') 'opened file : ', trim(tmpfil), ' unit', unitgreenlocks_arr(isppol)
4389 0 : write(message,'(a,a)') ch10,"# New record : First 40 bands"
4390 0 : call wrtout(unitgreenlocks_arr(isppol),message,'COLL')
4391 0 : do lsub=1,mbandc/40+1
4392 0 : do itau=1, ltau
4393 0 : write(message,'(2x,50(e10.3,2x))') tau(itau), &
4394 0 : & (real(loc_ks_tau(isppol,iband,itau)),iband=40*(lsub-1)+1,min(40*lsub,mbandc))
4395 0 : call wrtout(unitgreenlocks_arr(isppol),message,'COLL')
4396 : end do
4397 0 : write(message,'(2x,50(e10.3,2x))') beta, &
4398 0 : & ((-one-real(loc_ks_tau(isppol,iband,1))),iband=40*(lsub-1)+1,min(40*lsub,mbandc))
4399 0 : call wrtout(unitgreenlocks_arr(isppol),message,'COLL')
4400 0 : if(40*lsub<mbandc) then
4401 : write(message,'(a,a,i5,a,i5)') &
4402 0 : & ch10,"# Same record, Following bands : From ", &
4403 0 : & 40*(lsub)," to ",min(40*(lsub+1),mbandc)
4404 0 : call wrtout(unitgreenlocks_arr(isppol),message,'COLL')
4405 : end if
4406 : end do
4407 : ! call flush(unitgreenlocks_arr(isppol))
4408 : end do
4409 0 : ABI_FREE(unitgreenlocks_arr)
4410 : end if
4411 :
4412 : !Deallocations
4413 0 : ABI_FREE(loc_ks)
4414 0 : ABI_FREE(loc_ks_tau)
4415 0 : ABI_FREE(tau)
4416 :
4417 0 : end subroutine local_ks_green
4418 : !!***
4419 :
4420 : !!****f* m_green/compute_moments_ks
4421 : !! NAME
4422 : !! compute_moments_ks
4423 : !!
4424 : !! FUNCTION
4425 : !! Compute the high-frequency moments of the Green's function in KS space
4426 : !!
4427 : !! COPYRIGHT
4428 : !! Copyright (C) 1999-2026 ABINIT group (BAmadon)
4429 : !! This file is distributed under the terms of the
4430 : !! GNU General Public License, see ~abinit/COPYING
4431 : !! or http://www.gnu.org/copyleft/gpl.txt .
4432 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
4433 : !!
4434 : !! INPUTS
4435 : !! green <type(green_type)>= green function data
4436 : !! self <type(self_type)>= variables related to self-energy
4437 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
4438 : !! opt_self = 0 : do not use the self-energy
4439 : !! = 1 (default) : use the self-energy, by upfolding its moments
4440 : !! opt_log = 0 (default): the moments of G are computed in green%moments using self%moments
4441 : !! = 1 : the trace of the moments of log(G)+log(iw*Id)=-log(Id+(mu-sigma)/iw) is computed in green%trace_moments_log_ks
4442 : !! opt_quick_restart = 0 (default) : default behavior, recompute everything from scratch
4443 : !! = 1 : use precomputed quantities from fermi_green : assume that self%moments(i>=2) have
4444 : !! already been upfolded, that self%moments(1)-self%hdc has already been upfolded in
4445 : !! green%moments(2), and that (self%moments(1)-self%hdc+eigen_dft)**2 has been computed
4446 : !! in green%moments(3)
4447 : !! The conventions are: green%moments(i) -> ith order moment of the Green's function in the usual case
4448 : !! (the 1st one is not computed, since this is simply the identity)
4449 : !! self%moments(i) -> (i-1)th order moment of the self-energy (the 1st one is not upfolded, since we directly
4450 : !! upfold self%moments(1) - self%hdc
4451 : !! where the ith order moment corresponds to the factor in front of 1/(jw_n)^i in the asymptotic expansion
4452 : !!
4453 : !! OUTPUT
4454 : !!
4455 : !! SOURCE
4456 :
4457 0 : subroutine compute_moments_ks(green,self,paw_dmft,opt_self,opt_log,opt_quick_restart)
4458 :
4459 : !Arguments ------------------------------------
4460 : !scalars
4461 : type(green_type), intent(inout) :: green
4462 : type(self_type), intent(inout) :: self
4463 : type(paw_dmft_type), intent(in) :: paw_dmft
4464 : integer, optional, intent(in) :: opt_log,opt_quick_restart,opt_self
4465 : !Local variables ------------------------------
4466 : integer :: diag,i,ib,ierr,mkmem,natom,nsppol
4467 : integer :: optlog,optquickrestart,optself,shift
4468 : real(dp) :: dum,mu,mu2,mu3,mu4
4469 : complex(dp) :: trace_tmp
4470 0 : type(oper_type) :: oper(2)
4471 0 : real(dp), allocatable :: trace_loc(:,:)
4472 : character(len=500) :: message
4473 : !************************************************************************
4474 :
4475 0 : optlog = 0
4476 0 : if (present(opt_log)) optlog = opt_log
4477 :
4478 0 : optself = 1
4479 0 : if (present(opt_self)) optself = opt_self
4480 :
4481 0 : optquickrestart = 0
4482 0 : if (present(opt_quick_restart)) optquickrestart = opt_quick_restart
4483 :
4484 0 : if (optself == 0 .and. optquickrestart == 1) then
4485 0 : message = "Case optself=0 and optquickrestart=1 is not supposed to happen"
4486 0 : ABI_ERROR(message)
4487 : end if
4488 :
4489 0 : diag = 1 - optself
4490 0 : mkmem = green%moments(1)%nkpt
4491 0 : mu = paw_dmft%fermie
4492 0 : mu2 = mu * mu
4493 0 : mu3 = mu2 * mu
4494 0 : mu4 = mu3 * mu
4495 0 : natom = paw_dmft%natom
4496 0 : nsppol = paw_dmft%nsppol
4497 0 : shift = green%moments(1)%shiftk
4498 :
4499 0 : if (optlog == 1 .and. optquickrestart == 0) then
4500 0 : ABI_MALLOC(trace_loc,(nsppol+1,natom))
4501 : end if
4502 :
4503 0 : if (optself == 1 .and. optquickrestart == 0) then
4504 0 : call add_matlu(self%moments(1)%matlu(:),self%hdc%matlu(:),green%moments(2)%matlu(:),natom,-1)
4505 0 : call upfold_oper(green%moments(2),paw_dmft)
4506 0 : do i=2,self%nmoments
4507 0 : call upfold_oper(self%moments(i),paw_dmft)
4508 : end do ! i
4509 : end if ! optself=1
4510 :
4511 0 : do ib=1,paw_dmft%mbandc
4512 0 : if (optself == 1) then
4513 : green%moments(2)%ks(ib,ib,:,:) = green%moments(2)%ks(ib,ib,:,:) + &
4514 0 : & paw_dmft%eigen_dft(ib,1+shift:mkmem+shift,:) - mu
4515 : else
4516 0 : green%moments(2)%ks(ib,ib,:,:) = paw_dmft%eigen_dft(ib,1+shift:mkmem+shift,:) - mu
4517 : end if
4518 : ! Careful, we need to substract mu**2 below, not add it (think about it)
4519 0 : if (optquickrestart == 1) green%moments(3)%ks(ib,ib,:,:) = green%moments(3)%ks(ib,ib,:,:) - mu2
4520 : end do ! ib
4521 :
4522 0 : if (optquickrestart == 1) then
4523 0 : green%moments(3)%ks(:,:,:,:) = green%moments(3)%ks(:,:,:,:) - two*mu*green%moments(2)%ks(:,:,:,:)
4524 : else
4525 0 : call prod_oper(green%moments(2),green%moments(2),green%moments(3),1,opt_diag=diag)
4526 : end if ! optquickrestart
4527 0 : do i=3,green%nmoments-1
4528 0 : call prod_oper(green%moments(2),green%moments(i),green%moments(i+1),1,opt_diag=diag)
4529 : end do ! i
4530 :
4531 0 : if (optlog == 1 .and. optquickrestart == 0) then
4532 0 : do i=1,green%nmoments-1
4533 0 : call trace_oper(green%moments(i+1),dum,trace_loc(:,:),1,trace_ks_cmplx=trace_tmp)
4534 0 : green%trace_moments_log_ks(i) = trace_tmp / dble(i)
4535 : end do ! i
4536 0 : if (optself == 1) then
4537 0 : do i=2,self%nmoments
4538 0 : call trace_oper(self%moments(i),dum,trace_loc(:,:),1,trace_ks_cmplx=trace_tmp)
4539 0 : green%trace_moments_log_ks(i) = green%trace_moments_log_ks(i) + trace_tmp
4540 : end do ! i
4541 : end if ! optself=1
4542 : end if ! optlog=1 and optquickrestart=0
4543 :
4544 0 : if (optlog == 1 .and. optquickrestart == 1) then
4545 0 : green%trace_moments_log_ks(1) = green%trace_fermie(2) - mu*green%trace_fermie(1)
4546 : green%trace_moments_log_ks(2) = green%trace_fermie(3) + &
4547 0 : & half*(green%trace_fermie(4)-two*mu*green%trace_fermie(2)+mu2*green%trace_fermie(1))
4548 : green%trace_moments_log_ks(3) = green%trace_fermie(5) + green%trace_fermie(6) - &
4549 : & mu*green%trace_fermie(3) + third*(green%trace_fermie(7)+ &
4550 : & three*mu2*green%trace_fermie(2)-three*mu*green%trace_fermie(4)- &
4551 0 : & mu3*green%trace_fermie(1))
4552 : green%trace_moments_log_ks(4) = green%trace_fermie(8) + green%trace_fermie(9) - &
4553 : & mu*green%trace_fermie(5) + half*green%trace_fermie(10) + &
4554 : & green%trace_fermie(11) - two*mu*green%trace_fermie(6) + mu2*green%trace_fermie(3) + &
4555 : & quarter*(green%trace_fermie(12)-four*mu3*green%trace_fermie(2)+ &
4556 : & six*mu2*green%trace_fermie(4)-four*mu*green%trace_fermie(7)+ &
4557 0 : & mu4*green%trace_fermie(1))
4558 : end if ! optlog=1 and optquickrestart=1
4559 :
4560 0 : if (optself == 1) then
4561 0 : do i=1,2
4562 0 : call init_oper(paw_dmft,oper(i),nkpt=mkmem,shiftk=shift,opt_ksloc=1)
4563 : end do
4564 0 : call prod_oper(green%moments(2),self%moments(2),oper(1),1)
4565 0 : if (optlog == 1 .and. optquickrestart == 0) then
4566 0 : call trace_oper(oper(1),dum,trace_loc(:,:),1,trace_ks_cmplx=trace_tmp)
4567 0 : green%trace_moments_log_ks(3) = green%trace_moments_log_ks(3) + trace_tmp
4568 : end if
4569 0 : call prod_oper(self%moments(2),green%moments(2),oper(2),1)
4570 0 : oper(2)%ks(:,:,:,:) = oper(2)%ks(:,:,:,:) + oper(1)%ks(:,:,:,:)
4571 0 : green%moments(4)%ks(:,:,:,:) = oper(2)%ks(:,:,:,:) + green%moments(4)%ks(:,:,:,:)
4572 0 : call prod_oper(oper(2),green%moments(2),oper(1),1)
4573 0 : call prod_oper(green%moments(3),self%moments(2),oper(2),1)
4574 0 : if (optlog == 1 .and. optquickrestart == 0) then
4575 0 : call trace_oper(oper(2),dum,trace_loc(:,:),1,trace_ks_cmplx=trace_tmp)
4576 0 : green%trace_moments_log_ks(4) = green%trace_moments_log_ks(4) + trace_tmp
4577 : end if
4578 0 : oper(2)%ks(:,:,:,:) = oper(2)%ks(:,:,:,:) + oper(1)%ks(:,:,:,:)
4579 0 : call prod_oper(green%moments(2),self%moments(3),oper(1),1)
4580 0 : if (optlog == 1 .and. optquickrestart == 0) then
4581 0 : call trace_oper(oper(1),dum,trace_loc(:,:),1,trace_ks_cmplx=trace_tmp)
4582 0 : green%trace_moments_log_ks(4) = green%trace_moments_log_ks(4) + trace_tmp
4583 : end if
4584 0 : oper(2)%ks(:,:,:,:) = oper(2)%ks(:,:,:,:) + oper(1)%ks(:,:,:,:)
4585 0 : call prod_oper(self%moments(3),green%moments(2),oper(1),1)
4586 0 : oper(2)%ks(:,:,:,:) = oper(2)%ks(:,:,:,:) + oper(1)%ks(:,:,:,:)
4587 0 : call prod_oper(self%moments(2),self%moments(2),oper(1),1)
4588 0 : if (optlog == 1 .and. optquickrestart == 0) then
4589 0 : call trace_oper(oper(1),dum,trace_loc(:,:),1,trace_ks_cmplx=trace_tmp)
4590 0 : green%trace_moments_log_ks(4) = green%trace_moments_log_ks(4) + trace_tmp*half
4591 : end if
4592 0 : oper(2)%ks(:,:,:,:) = oper(2)%ks(:,:,:,:) + oper(1)%ks(:,:,:,:)
4593 0 : green%moments(5)%ks(:,:,:,:) = oper(2)%ks(:,:,:,:) + green%moments(5)%ks(:,:,:,:)
4594 0 : do i=2,self%nmoments
4595 0 : green%moments(i+1)%ks(:,:,:,:) = green%moments(i+1)%ks(:,:,:,:) + self%moments(i)%ks(:,:,:,:)
4596 : end do ! i
4597 0 : do i=1,2
4598 0 : call destroy_oper(oper(i))
4599 : end do
4600 : end if ! optself>=1
4601 :
4602 0 : if (optlog == 1 .and. optquickrestart == 0) then
4603 0 : call xmpi_sum(green%trace_moments_log_ks(1:green%nmoments-1),green%distrib%comm_kpt,ierr)
4604 : end if
4605 :
4606 0 : ABI_SFREE(trace_loc)
4607 :
4608 0 : end subroutine compute_moments_ks
4609 : !!***
4610 :
4611 : !!****f* m_green/compute_trace_moments_ks
4612 : !! NAME
4613 : !! compute_trace_moments_ks
4614 : !!
4615 : !! FUNCTION
4616 : !! Precompute some useful quantities for the calculation of the trace of the moments (for
4617 : !! Fermi level search)
4618 : !!
4619 : !! COPYRIGHT
4620 : !! Copyright (C) 1999-2026 ABINIT group (BAmadon)
4621 : !! This file is distributed under the terms of the
4622 : !! GNU General Public License, see ~abinit/COPYING
4623 : !! or http://www.gnu.org/copyleft/gpl.txt .
4624 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
4625 : !!
4626 : !! INPUTS
4627 : !! green <type(green_type)>= green function data
4628 : !! self <type(self_type)>= variables related to self-energy
4629 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
4630 : !! The conventions are:
4631 : !! trace_fermie(1) = nsppol*mbandc (=tr(Id))
4632 : !! trace_fermie(2) = tr(m0) (with m0=self%moments(1)-self%hdc+eigen_dft*Id
4633 : !! trace_fermie(3) = tr(self%moments(2))
4634 : !! trace_fermie(4) = tr(m0**2)
4635 : !! trace_fermie(5) = tr(self%moments(3))
4636 : !! trace_fermie(6) = tr(m0*self%moments(2))
4637 : !! trace_fermie(7) = tr(m0**3)
4638 : !! trace_fermie(8) = tr(self%moments(4))
4639 : !! trace_fermie(9) = tr(m0*self%moments(3))
4640 : !! trace_fermie(10) = tr(self%moments(2)**2)
4641 : !! trace_fermie(11) = tr(m0**2*self%moments(2))
4642 : !! trace_fermie(12) = tr(m0**4)
4643 : !!
4644 : !! OUTPUT
4645 : !!
4646 : !! SOURCE
4647 :
4648 0 : subroutine compute_trace_moments_ks(green,self,paw_dmft)
4649 :
4650 : !Arguments ------------------------------------
4651 : !scalars
4652 : type(green_type), intent(inout) :: green
4653 : type(self_type), intent(inout) :: self
4654 : type(paw_dmft_type), intent(in) :: paw_dmft
4655 : !Local variables ------------------------------
4656 : integer :: i,ib,ierr,mkmem,natom,nsppol,shift
4657 : real(dp) :: dum
4658 0 : type(oper_type) :: oper_tmp
4659 0 : real(dp), allocatable :: trace_loc(:,:)
4660 : !************************************************************************
4661 :
4662 0 : mkmem = green%distrib%nkpt_mem(green%distrib%me_kpt+1)
4663 0 : natom = paw_dmft%natom
4664 0 : nsppol = paw_dmft%nsppol
4665 0 : shift = green%distrib%shiftk
4666 :
4667 0 : ABI_MALLOC(trace_loc,(nsppol+1,natom))
4668 :
4669 0 : call init_oper(paw_dmft,oper_tmp,nkpt=mkmem,shiftk=shift,opt_ksloc=1)
4670 :
4671 0 : call add_matlu(self%moments(1)%matlu(:),self%hdc%matlu(:),green%moments(2)%matlu(:),natom,-1)
4672 0 : call upfold_oper(green%moments(2),paw_dmft)
4673 0 : do i=2,self%nmoments
4674 0 : call upfold_oper(self%moments(i),paw_dmft)
4675 : end do ! i
4676 :
4677 0 : oper_tmp%ks(:,:,:,:) = green%moments(2)%ks(:,:,:,:)
4678 :
4679 0 : do ib=1,paw_dmft%mbandc
4680 0 : oper_tmp%ks(ib,ib,:,:) = oper_tmp%ks(ib,ib,:,:) + paw_dmft%eigen_dft(ib,1+shift:mkmem+shift,:)
4681 : end do ! ib
4682 :
4683 0 : call prod_oper(oper_tmp,oper_tmp,green%moments(3),1)
4684 0 : call trace_oper(oper_tmp,dum,trace_loc(:,:),1,trace_ks_cmplx=green%trace_fermie(2))
4685 0 : call trace_oper(self%moments(2),dum,trace_loc(:,:),1,trace_ks_cmplx=green%trace_fermie(3))
4686 0 : call trace_oper(green%moments(3),dum,trace_loc(:,:),1,trace_ks_cmplx=green%trace_fermie(4))
4687 0 : call trace_oper(self%moments(3),dum,trace_loc(:,:),1,trace_ks_cmplx=green%trace_fermie(5))
4688 0 : call trace_prod_oper(oper_tmp,self%moments(2),green%trace_fermie(6))
4689 0 : call trace_prod_oper(oper_tmp,green%moments(3),green%trace_fermie(7))
4690 0 : call trace_oper(self%moments(4),dum,trace_loc(:,:),1,trace_ks_cmplx=green%trace_fermie(8))
4691 0 : call trace_prod_oper(oper_tmp,self%moments(3),green%trace_fermie(9))
4692 0 : call trace_prod_oper(self%moments(2),self%moments(2),green%trace_fermie(10))
4693 0 : call trace_prod_oper(green%moments(3),self%moments(2),green%trace_fermie(11))
4694 0 : call trace_prod_oper(green%moments(3),green%moments(3),green%trace_fermie(12))
4695 :
4696 0 : call destroy_oper(oper_tmp)
4697 0 : ABI_FREE(trace_loc)
4698 :
4699 0 : call xmpi_sum(green%trace_fermie(2:12),green%distrib%comm_kpt,ierr)
4700 :
4701 0 : end subroutine compute_trace_moments_ks
4702 : !!***
4703 :
4704 : !!****f* m_green/compute_moments_loc
4705 : !! NAME
4706 : !! compute_moments_loc
4707 : !!
4708 : !! FUNCTION
4709 : !! Computes the high-frequency moments in local space.
4710 : !!
4711 : !! INPUTS
4712 : !! green <type(green_type)>= green function data
4713 : !! self <type(self_type)>= variables related to self-energy
4714 : !! energy_level = local energy levels
4715 : !! weiss <type(green_type)>= weiss field
4716 : !! option = 0 : computes the hybridization moments in weiss%moments up to 4th order (including
4717 : !! the spurious 0th order moment, which is removed from the weiss field later),
4718 : !! using the moments of the self-energy and Green's function as inputs
4719 : !! = 1 : compute the self-energy moments in self%moments using the moments of the hybridization
4720 : !! and the Green's function (assuming the spurious 0th order moment of the hybridization
4721 : !! has been removed)
4722 : !! opt_log = if set to 1, also computes the trace of the moments of log(G)+log(iw*Id) and
4723 : !! log(G0)+log(iw*Id) in green%trace_moments_log_loc and weiss%trace_moments_log_loc
4724 : !! shift_mu = shift to apply to the chemical potential for the moments of G0
4725 : !!
4726 : !! OUTPUTS
4727 : !!
4728 : !! SOURCE
4729 :
4730 0 : subroutine compute_moments_loc(green,self,energy_level,weiss,option,opt_log,shift_mu)
4731 :
4732 : !Arguments ------------------------------------
4733 : type(green_type), intent(inout) :: green,weiss
4734 : type(self_type), intent(inout) :: self
4735 : type(oper_type), target, intent(in) :: energy_level
4736 : integer, intent(in) :: option
4737 : integer, optional, intent(in) :: opt_log
4738 : real(dp), optional, intent(in) :: shift_mu
4739 : !Local variables ------------------------------
4740 : integer :: i,natom,nspinor,nsppol,optlog
4741 : complex(dp) :: trace
4742 0 : type(matlu_type), target, allocatable :: level_shift(:)
4743 0 : integer, allocatable :: lpawu(:)
4744 0 : complex(dp), allocatable :: shift(:),trace_loc(:)
4745 0 : type(matlu_type), allocatable :: matlu(:,:)
4746 : !************************************************************************
4747 :
4748 0 : optlog = 0
4749 0 : if (present(opt_log)) optlog = opt_log
4750 :
4751 0 : natom = energy_level%natom
4752 0 : nspinor = energy_level%nspinor
4753 0 : nsppol = energy_level%nsppol
4754 :
4755 0 : ABI_MALLOC(trace_loc,(natom))
4756 :
4757 0 : ABI_MALLOC(lpawu,(natom))
4758 0 : ABI_MALLOC(matlu,(natom,6))
4759 :
4760 0 : lpawu(:) = energy_level%matlu(:)%lpawu ! to prevent the creation of a temporary array when calling init_matlu
4761 :
4762 0 : do i=1,6
4763 0 : call init_matlu(natom,nspinor,nsppol,lpawu(:),matlu(:,i))
4764 : end do
4765 :
4766 0 : if (option == 0) then
4767 0 : call add_matlu(green%moments(2)%matlu(:),energy_level%matlu(:),matlu(:,1),natom,-1)
4768 0 : call add_matlu(matlu(:,1),self%moments(1)%matlu(:),weiss%moments(1)%matlu(:),natom,-1)
4769 : else
4770 0 : call add_matlu(green%moments(2)%matlu(:),energy_level%matlu(:),self%moments(1)%matlu(:),natom,-1)
4771 : end if ! option
4772 :
4773 0 : if (optlog > 0) then
4774 :
4775 0 : ABI_MALLOC(level_shift,(natom))
4776 0 : ABI_MALLOC(shift,(natom))
4777 0 : call init_matlu(natom,nspinor,nsppol,lpawu(:),level_shift(:))
4778 0 : call copy_matlu(energy_level%matlu(:),level_shift(:),natom)
4779 :
4780 0 : shift(:) = cmplx(shift_mu,zero,kind=dp)
4781 0 : call shift_matlu(level_shift(:),natom,shift(:),signe=1)
4782 :
4783 0 : call trace_matlu(level_shift(:),natom,itau=0,trace=weiss%trace_moments_log_loc(1))
4784 0 : do i=2,green%nmoments-1
4785 0 : call trace_matlu(weiss%moments(i)%matlu(:),natom,itau=0,trace=weiss%trace_moments_log_loc(i))
4786 : end do ! i
4787 0 : call prod_matlu(level_shift(:),level_shift(:),matlu(:,1),natom)
4788 0 : call trace_matlu(matlu(:,1),natom,itau=0,trace=trace)
4789 0 : weiss%trace_moments_log_loc(2) = weiss%trace_moments_log_loc(2) + trace*half
4790 0 : call trace_prod_matlu(level_shift(:),weiss%moments(2)%matlu(:),natom,trace_loc(:),trace_tot=trace)
4791 0 : weiss%trace_moments_log_loc(3) = weiss%trace_moments_log_loc(3) + trace
4792 0 : call trace_prod_matlu(level_shift(:),matlu(:,1),natom,trace_loc(:),trace_tot=trace)
4793 0 : weiss%trace_moments_log_loc(3) = weiss%trace_moments_log_loc(3) + trace*third
4794 0 : call trace_prod_matlu(level_shift(:),weiss%moments(3)%matlu(:),natom,trace_loc(:),trace_tot=trace)
4795 0 : weiss%trace_moments_log_loc(4) = weiss%trace_moments_log_loc(4) + trace
4796 0 : call trace_prod_matlu(weiss%moments(2)%matlu(:),weiss%moments(2)%matlu(:),natom,trace_loc(:),trace_tot=trace)
4797 0 : weiss%trace_moments_log_loc(4) = weiss%trace_moments_log_loc(4) + trace*half
4798 0 : call trace_prod_matlu(weiss%moments(2)%matlu(:),matlu(:,1),natom,trace_loc(:),trace_tot=trace)
4799 0 : weiss%trace_moments_log_loc(4) = weiss%trace_moments_log_loc(4) + trace
4800 0 : call trace_prod_matlu(matlu(:,1),matlu(:,1),natom,trace_loc(:),trace_tot=trace)
4801 0 : weiss%trace_moments_log_loc(4) = weiss%trace_moments_log_loc(4) + trace*quarter
4802 :
4803 0 : call trace_matlu(green%moments(2)%matlu(:),natom,itau=0,trace=green%trace_moments_log_loc(1))
4804 :
4805 0 : call destroy_matlu(level_shift(:),natom)
4806 0 : ABI_FREE(level_shift)
4807 0 : ABI_FREE(shift)
4808 :
4809 : end if ! optlog
4810 :
4811 0 : call prod_matlu(green%moments(2)%matlu(:),green%moments(2)%matlu(:),matlu(:,1),natom) ! matlu=m0m0
4812 :
4813 0 : call add_matlu(green%moments(3)%matlu(:),matlu(:,1),matlu(:,2),natom,-1) ! matlu2=m1
4814 :
4815 0 : if (optlog > 0) then
4816 0 : call trace_matlu(matlu(:,2),natom,itau=0,trace=green%trace_moments_log_loc(2))
4817 0 : call trace_matlu(matlu(:,1),natom,itau=0,trace=trace)
4818 0 : green%trace_moments_log_loc(2) = green%trace_moments_log_loc(2) + half*trace
4819 : end if ! optlog
4820 :
4821 0 : if (option == 0) then
4822 0 : call add_matlu(matlu(:,2),self%moments(2)%matlu(:),weiss%moments(2)%matlu(:),natom,-1)
4823 0 : else if (option == 1) then
4824 0 : call add_matlu(matlu(:,2),weiss%moments(2)%matlu(:),self%moments(2)%matlu(:),natom,-1)
4825 : end if ! option
4826 :
4827 0 : call prod_matlu(green%moments(2)%matlu(:),matlu(:,1),matlu(:,3),natom) ! matlu3=m0m0m0
4828 0 : call add_matlu(green%moments(4)%matlu(:),matlu(:,3),matlu(:,1),natom,-1) ! matlu=green%moments(4)-m0m0m0
4829 :
4830 0 : if (optlog > 0) then
4831 0 : call trace_matlu(matlu(:,3),natom,itau=0,trace=trace)
4832 0 : green%trace_moments_log_loc(3) = trace * third
4833 : end if ! optlog
4834 :
4835 0 : call prod_matlu(green%moments(2)%matlu(:),matlu(:,2),matlu(:,4),natom) ! matlu4=m0m1
4836 0 : call prod_matlu(matlu(:,2),green%moments(2)%matlu(:),matlu(:,5),natom) ! matlu5=m1m0
4837 :
4838 0 : call add_matlu(matlu(:,4),matlu(:,5),matlu(:,6),natom,1) ! matlu6=m0m1+m1m0
4839 :
4840 0 : if (optlog > 0) then
4841 0 : call trace_matlu(matlu(:,4),natom,itau=0,trace=trace)
4842 0 : green%trace_moments_log_loc(3) = trace + green%trace_moments_log_loc(3)
4843 : end if ! optlog
4844 :
4845 0 : call add_matlu(matlu(:,1),matlu(:,6),matlu(:,5),natom,-1) ! matlu5=m2
4846 :
4847 0 : if (optlog > 0) then
4848 0 : call trace_matlu(matlu(:,5),natom,itau=0,trace=trace)
4849 0 : green%trace_moments_log_loc(3) = trace + green%trace_moments_log_loc(3)
4850 : end if ! optlog
4851 :
4852 0 : if (option == 0) then
4853 0 : call add_matlu(matlu(:,5),self%moments(3)%matlu(:),weiss%moments(3)%matlu(:),natom,-1)
4854 0 : else if (option == 1) then
4855 0 : call add_matlu(matlu(:,5),weiss%moments(3)%matlu(:),self%moments(3)%matlu(:),natom,-1)
4856 : end if ! option
4857 :
4858 0 : call prod_matlu(green%moments(2)%matlu(:),matlu(:,3),matlu(:,1),natom) ! matlu=m0m0m0m0
4859 0 : call prod_matlu(matlu(:,6),green%moments(2)%matlu(:),matlu(:,3),natom) ! matlu3 = m0m1m0+m1m0m0
4860 0 : call prod_matlu(green%moments(2)%matlu(:),matlu(:,4),matlu(:,6),natom) ! matlu6 = m0m0m1
4861 0 : call add_matlu(matlu(:,3),matlu(:,6),matlu(:,4),natom,1) ! matlu4 = m0m1m0+m1m0m0+m0m0m1
4862 :
4863 0 : if (optlog > 0) then
4864 0 : call trace_matlu(matlu(:,4),natom,itau=0,trace=trace)
4865 0 : green%trace_moments_log_loc(4) = trace * third
4866 0 : call trace_matlu(matlu(:,1),natom,itau=0,trace=trace)
4867 0 : green%trace_moments_log_loc(4) = green%trace_moments_log_loc(4) + trace*quarter
4868 : end if ! optlog
4869 :
4870 0 : call add_matlu(matlu(:,1),matlu(:,4),matlu(:,3),natom,1) ! matlu3 = m0m1m0+m1m0m0+m0m0m1+m0m0m0m0
4871 :
4872 0 : call prod_matlu(green%moments(2)%matlu(:),matlu(:,5),matlu(:,1),natom) ! matlu=m0m2
4873 0 : call prod_matlu(matlu(:,5),green%moments(2)%matlu(:),matlu(:,4),natom) ! matlu4=m2m0
4874 0 : call add_matlu(matlu(:,1),matlu(:,4),matlu(:,5),natom,1) ! matlu5=m0m2+m2m0
4875 0 : call prod_matlu(matlu(:,2),matlu(:,2),matlu(:,1),natom) ! matlu=m1m1
4876 0 : call add_matlu(matlu(:,1),matlu(:,5),matlu(:,2),natom,1) ! matlu2=m1m1+m0m2+m2m0
4877 :
4878 0 : if (optlog > 0) then
4879 0 : call trace_matlu(matlu(:,2),natom,itau=0,trace=trace)
4880 0 : green%trace_moments_log_loc(4) = green%trace_moments_log_loc(4) + trace*half
4881 : end if ! optlog
4882 :
4883 0 : call add_matlu(green%moments(5)%matlu(:),matlu(:,2),matlu(:,1),natom,-1)
4884 0 : call add_matlu(matlu(:,1),matlu(:,3),matlu(:,2),natom,-1) ! matlu2=m3
4885 :
4886 0 : if (optlog > 0) then
4887 0 : call trace_matlu(matlu(:,2),natom,itau=0,trace=trace)
4888 0 : green%trace_moments_log_loc(4) = green%trace_moments_log_loc(4) + trace
4889 : end if ! optlog
4890 :
4891 0 : if (option == 0) then
4892 0 : call add_matlu(matlu(:,2),self%moments(4)%matlu(:),weiss%moments(4)%matlu(:),natom,-1)
4893 0 : else if (option == 1) then
4894 0 : call add_matlu(matlu(:,2),weiss%moments(4)%matlu(:),self%moments(4)%matlu(:),natom,-1)
4895 : end if ! option
4896 :
4897 0 : do i=1,6
4898 0 : call destroy_matlu(matlu(:,i),natom)
4899 : end do
4900 :
4901 0 : ABI_FREE(lpawu)
4902 0 : ABI_FREE(matlu)
4903 0 : ABI_FREE(trace_loc)
4904 :
4905 0 : end subroutine compute_moments_loc
4906 : !!***
4907 :
4908 : !!****f* m_green/occup_fd
4909 : !! NAME
4910 : !! occup_fd
4911 : !!
4912 : !! FUNCTION
4913 : !! Computes the Fermi-Dirac occupation
4914 : !!
4915 : !! INPUTS
4916 : !! eig = eigenvalues
4917 : !! fermie = Fermi level
4918 : !! temp = temperature
4919 : !!
4920 : !! OUTPUT
4921 : !!
4922 : !!
4923 : !! SOURCE
4924 :
4925 9852 : function occup_fd(eig,fermie,temp)
4926 :
4927 : !Arguments ------------------------------------
4928 : !type
4929 : ! Integrate analytic tail 1/(iw-mu)
4930 : real(dp), intent(in) :: eig,fermie,temp
4931 : real(dp) :: occup_fd
4932 : !Local variables-------------------------------
4933 : ! *********************************************************************
4934 :
4935 9852 : if ((eig-fermie) > zero) then
4936 5067 : occup_fd = exp(-(eig-fermie)/temp) / (one+exp(-(eig-fermie)/temp))
4937 : else
4938 4785 : occup_fd = one / (one+exp((eig-fermie)/temp))
4939 : end if ! eig-fermie>0
4940 :
4941 9852 : end function occup_fd
4942 :
4943 0 : END MODULE m_green
4944 : !!***
|
