Line data Source code
1 : !!****m* ABINIT/m_self
2 : !! NAME
3 : !! m_self
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_self
26 :
27 : use defs_basis
28 : use m_errors
29 : use m_abicore
30 :
31 : use m_datafordmft, only : compute_levels
32 : use m_fstrings, only : int2char4
33 : use m_hu, only : hu_type
34 : use m_io_tools, only : get_unit
35 : use m_matlu, only : add_matlu,copy_matlu,destroy_matlu,diag_matlu,fac_matlu, &
36 : & init_matlu,matlu_type,print_matlu,rotate_matlu,xmpi_matlu,zero_matlu
37 : use m_oper, only : destroy_oper,gather_oper,init_oper,oper_type,print_oper
38 : use m_paw_dmft, only : mpi_distrib_dmft_type,paw_dmft_type
39 : use m_paw_exactDC, only : compute_exactDC
40 : use m_pawtab, only : pawtab_type
41 : use m_xmpi, only : xmpi_bcast,xmpi_sum
42 :
43 : #ifdef HAVE_GPU_MARKERS
44 : use m_nvtx_data
45 : #endif
46 :
47 : implicit none
48 :
49 : private
50 :
51 : public :: alloc_self
52 : public :: initialize_self
53 : public :: destroy_self
54 : public :: print_self
55 : public :: rw_self
56 : public :: dc_self
57 : public :: new_self
58 : !public :: make_qmcshift_self
59 : public :: selfreal2imag_self
60 :
61 : !!***
62 :
63 : !!****t* m_self/self_type
64 : !! NAME
65 : !! self_type
66 : !!
67 : !! FUNCTION
68 : !! This structured datatype contains the necessary data
69 : !!
70 : !! SOURCE
71 :
72 : type, public :: self_type ! for each atom
73 :
74 : integer :: dmft_nwli
75 : ! Linear index of the last imaginary frequency
76 :
77 : integer :: dmft_nwlo
78 : ! Number of imaginary frequencies
79 :
80 : integer :: has_moments
81 : ! =1 if the high-frequency moments are computed
82 :
83 : integer :: iself_cv
84 : ! Integer for convergence of self-energy
85 :
86 : integer :: nmoments
87 : ! Number of high-frequency moments which will be computed
88 :
89 : integer :: nw
90 : ! Number of frequencies (equal to dmft_nwlo only if w_type="imag")
91 :
92 : character(len=4) :: w_type
93 : ! Type of frequencies used ("real" or "imag")
94 :
95 : !real(dp), allocatable :: qmc_shift(:)
96 : ! value of frequencies
97 :
98 : !real(dp), allocatable :: qmc_xmu(:)
99 : ! value of frequencies
100 :
101 : type(oper_type) :: hdc
102 : ! Operator for double counting
103 :
104 : type(oper_type), allocatable :: moments(:)
105 : ! High-frequency moments
106 :
107 : type(oper_type), allocatable :: oper(:)
108 : ! Operator for self-energy, for each frequency
109 :
110 : real(dp), ABI_CONTIGUOUS pointer :: omega(:) => null()
111 : ! Value of frequencies
112 :
113 : type(mpi_distrib_dmft_type), pointer :: distrib => null()
114 : ! Datastructure for MPI parallelization
115 :
116 : end type self_type
117 : !!***
118 :
119 : !----------------------------------------------------------------------
120 :
121 :
122 : CONTAINS !========================================================================================
123 : !!***
124 :
125 : !!****f* m_self/alloc_self
126 : !! NAME
127 : !! alloc_self
128 : !!
129 : !! FUNCTION
130 : !! Allocate variables used in type self_type.
131 : !!
132 : !! INPUTS
133 : !! self <type(self_type)>= variables related to self-energy
134 : !! paw_dmft <type(paw_dmft_type)> = variables related to self-consistent DFT+DMFT calculations.
135 : !! opt_oper = 1 Allocate only quantities in the KS basis.
136 : !! 2 Allocate only quantities in the local basis.
137 : !! 3 Allocate quantities in both the KS and local basis.
138 : !! wtype = "real" Self energy will be computed for real frequencies
139 : !! = "imag" (default) Self energy will be computed for imaginary frequencies
140 : !! opt_moments = 1 to allocate the high-frequency moments
141 : !!
142 : !! OUTPUTS
143 : !! self <type(self_type)>= variables related to self-energy
144 : !!
145 : !! SOURCE
146 :
147 211 : subroutine alloc_self(self,paw_dmft,opt_oper,wtype,opt_moments)
148 :
149 : !Arguments ------------------------------------
150 : type(self_type), intent(inout) :: self
151 : type(paw_dmft_type), target, intent(in) :: paw_dmft
152 : integer, optional, intent(in) :: opt_moments,opt_oper
153 : character(len=4), optional :: wtype
154 : !Local variables ------------------------------------
155 : integer :: i,ifreq,mkmem,optmoments,optoper,shift
156 : !************************************************************************
157 :
158 211 : optoper = 2
159 211 : self%w_type = "imag"
160 211 : self%nmoments = 0
161 211 : optmoments = 0
162 :
163 211 : if (present(opt_oper)) optoper = opt_oper
164 211 : if (present(wtype)) self%w_type = wtype
165 211 : if (present(opt_moments)) optmoments = opt_moments
166 :
167 211 : self%has_moments = optmoments
168 :
169 211 : if (self%w_type == "imag") then
170 207 : self%nw = paw_dmft%dmft_nwlo
171 207 : self%omega => paw_dmft%omega_lo(:)
172 207 : self%distrib => paw_dmft%distrib
173 4 : else if (self%w_type == "real") then
174 4 : self%nw = size(paw_dmft%omega_r(:))
175 4 : self%omega => paw_dmft%omega_r(:)
176 4 : self%distrib => paw_dmft%distrib_r
177 : end if ! w_type
178 :
179 211 : self%dmft_nwlo = paw_dmft%dmft_nwlo
180 211 : self%dmft_nwli = paw_dmft%dmft_nwli
181 211 : self%iself_cv = 0
182 :
183 211 : call init_oper(paw_dmft,self%hdc,opt_ksloc=optoper)
184 :
185 68311 : ABI_MALLOC(self%oper,(self%nw))
186 67889 : do ifreq=1,self%nw
187 67889 : call init_oper(paw_dmft,self%oper(ifreq),opt_ksloc=optoper)
188 : end do ! ifreq
189 :
190 211 : if (optmoments == 1) then
191 0 : self%nmoments = 4
192 0 : shift = self%distrib%shiftk
193 0 : mkmem = self%distrib%nkpt_mem(self%distrib%me_kpt+1)
194 0 : ABI_MALLOC(self%moments,(self%nmoments))
195 0 : call init_oper(paw_dmft,self%moments(1),nkpt=mkmem,shiftk=shift,opt_ksloc=2)
196 0 : do i=2,self%nmoments
197 0 : call init_oper(paw_dmft,self%moments(i),nkpt=mkmem,shiftk=shift,opt_ksloc=3)
198 : end do ! i
199 : end if ! optmoments
200 :
201 : !if (paw_dmft%dmft_solv == 4) then
202 : ! ABI_MALLOC(self%qmc_shift,(paw_dmft%natom))
203 : ! ABI_MALLOC(self%qmc_xmu,(paw_dmft%natom))
204 : ! self%qmc_shift(:) = zero
205 : ! self%qmc_xmu(:) = zero
206 : !end if ! dmft_solv=4
207 :
208 211 : end subroutine alloc_self
209 : !!***
210 :
211 : !!****f* m_self/initialize_self
212 : !! NAME
213 : !! initialize_self
214 : !!
215 : !! FUNCTION
216 : !! Initialize self-energy.
217 : !!
218 : !! INPUTS
219 : !! self <type(self_type)>= variables related to self-energy
220 : !! paw_dmft <type(paw_dmft_type)> = variables related to self-consistent DFT+DMFT calculations.
221 : !! wtype = "real" Self energy will be computed for real frequencies
222 : !! = "imag" (default) Self energy will be computed for imaginary frequencies
223 : !! opt_moments = 1 to compute the high frequency moments
224 : !!
225 : !! OUTPUTS
226 : !! self <type(self_type)>= variables related to self-energy
227 : !!
228 : !!
229 : !! SOURCE
230 :
231 211 : subroutine initialize_self(self,paw_dmft,wtype,opt_moments)
232 :
233 : !Arguments ------------------------------------
234 : type(self_type), intent(inout) :: self
235 : type(paw_dmft_type), intent(in) :: paw_dmft
236 : character(len=4), optional, intent(in) :: wtype
237 : integer, optional, intent(in) :: opt_moments
238 : !Local variables ------------------------------------
239 : ! character(len=500) :: message
240 : integer :: optmoments
241 : character(len=4) :: wtype2
242 : !************************************************************************
243 :
244 211 : optmoments = 0
245 211 : wtype2 = "imag"
246 211 : if (present(wtype)) wtype2 = wtype
247 211 : if (present(opt_moments)) optmoments = opt_moments
248 :
249 211 : call alloc_self(self,paw_dmft,opt_oper=2,wtype=wtype2,opt_moments=optmoments) ! opt_oper=1 is not useful and not implemented
250 : ! if(paw_dmft%dmft_rslf==1.and.opt_read==1) then
251 : ! call rw_self(cryst_struc,self,mpi_enreg,paw_dmft,pawtab,pawprtvol=2,opt_rw=1)
252 : ! endif
253 : ! write(message,'(a,a)') ch10," Self-energy for large frequency is"
254 : ! call wrtout(std_out,message,'COLL')
255 : ! call print_matlu(self%oper(paw_dmft%dmft_nwlo)%matlu, &
256 : !& paw_dmft%natom,3)
257 :
258 211 : end subroutine initialize_self
259 : !!***
260 :
261 : !!****f* m_self/destroy_self
262 : !! NAME
263 : !! destroy_self
264 : !!
265 : !! FUNCTION
266 : !! Deallocate self
267 : !!
268 : !! INPUTS
269 : !! self <type(self_type)>= variables related to self-energy
270 : !!
271 : !! OUTPUT
272 : !!
273 : !! SOURCE
274 :
275 211 : subroutine destroy_self(self)
276 :
277 : !Arguments ------------------------------------
278 : type(self_type), intent(inout) :: self
279 : !Local variables-------------------------------
280 : integer :: i,ifreq
281 : ! *********************************************************************
282 :
283 211 : if (allocated(self%oper)) then
284 67889 : do ifreq=1,self%nw
285 67889 : call destroy_oper(self%oper(ifreq))
286 : end do
287 67889 : ABI_FREE(self%oper)
288 : end if
289 :
290 211 : call destroy_oper(self%hdc)
291 :
292 211 : if (allocated(self%moments)) then
293 0 : do i=1,self%nmoments
294 0 : call destroy_oper(self%moments(i))
295 : end do ! i
296 0 : ABI_FREE(self%moments)
297 : end if
298 : !if (allocated(self%qmc_shift)) ABI_FREE(self%qmc_shift)
299 : !if (allocated(self%qmc_xmu)) ABI_FREE(self%qmc_xmu)
300 211 : self%distrib => null()
301 211 : self%omega => null()
302 :
303 211 : end subroutine destroy_self
304 : !!***
305 :
306 : !!****f* m_self/print_self
307 : !! NAME
308 : !! print_self
309 : !!
310 : !! FUNCTION
311 : !! Print self-energy
312 : !!
313 : !! INPUTS
314 : !! self <type(self_type)>= variables related to self-energy
315 : !! prtdc = print double counting if equal to "print_dc"
316 : !! paw_dmft <type(paw_dmft_type)> = variables related to self-consistent DFT+DMFT calculations.
317 : !! prtopt = integer which specifies the amount of printing in the subroutine called
318 : !!
319 : !! OUTPUT
320 : !! self <type(self_type)>= variables related to self-energy
321 : !!
322 : !! SOURCE
323 :
324 339 : subroutine print_self(self,prtdc,paw_dmft,prtopt)
325 :
326 : !Arguments ------------------------------------
327 : type(paw_dmft_type), intent(in) :: paw_dmft
328 : type(self_type), intent(in) :: self
329 : character(len=*), intent(in) :: prtdc
330 : integer, intent(in) :: prtopt
331 : !Local variables-------------------------------
332 : character(len=500) :: message
333 : ! *********************************************************************
334 :
335 339 : write(message,'(2a)') ch10," == The self-energy for smallest frequency is == "
336 339 : call wrtout(std_out,message,'COLL')
337 339 : call print_oper(self%oper(1),1,paw_dmft,prtopt)
338 : ! write(message,'(2a)') ch10," == The self-energy for small (3) frequency is == "
339 : ! call wrtout(std_out,message,'COLL')
340 : ! call print_oper(self%oper(3),1,paw_dmft,prtopt)
341 339 : write(message,'(2a)') ch10," == The self-energy for largest frequency is == "
342 339 : call wrtout(std_out,message,'COLL')
343 339 : call print_oper(self%oper(self%nw),1,paw_dmft,prtopt)
344 339 : if (prtdc == "print_dc") then
345 339 : write(message,'(2a)') ch10," == The double counting potential is == "
346 339 : call wrtout(std_out,message,'COLL')
347 339 : call print_matlu(self%hdc%matlu(:),paw_dmft%natom,prtopt)
348 : end if ! prtdc
349 :
350 339 : end subroutine print_self
351 : !!***
352 :
353 : !!****f* m_self/dc_self
354 : !! NAME
355 : !! dc_self
356 : !!
357 : !! FUNCTION
358 : !! Computes the double counting
359 : !!
360 : !! INPUTS
361 : !! charge_loc : local charge for each polarization and each atom
362 : !! hu <type(hu_type)>= U interaction
363 : !! paw_dmft <type(paw_dmft_type)> = variables related to self-consistent DFT+DMFT calculations.
364 : !! pawtab <type(pawtab_type)>=paw tabulated starting data
365 : !! occ_matlu : local occupation matrix
366 : !!
367 : !! OUTPUT
368 : !! hdc : double counting
369 : !!
370 : !! SOURCE
371 :
372 199 : subroutine dc_self(charge_loc,hdc,hu,paw_dmft,pawtab,occ_matlu)
373 :
374 : !Arguments ------------------------------------
375 : !type
376 : !type(crystal_t),intent(in) :: cryst_struc
377 : type(paw_dmft_type), intent(inout) :: paw_dmft
378 : type(matlu_type), intent(inout) :: hdc(paw_dmft%natom)
379 : real(dp), intent(in) :: charge_loc(paw_dmft%nsppol+1,paw_dmft%natom)
380 : type(hu_type), intent(in) :: hu(paw_dmft%ntypat)
381 : type(pawtab_type), intent(inout) :: pawtab(paw_dmft%ntypat)
382 : type(matlu_type), intent(in) :: occ_matlu(paw_dmft%natom)
383 : !type(hu_type),intent(inout) :: hu(cryst_struc%ntypat)
384 : !integer, intent(in) :: dmft_dc
385 : !Local variables-------------------------------
386 : integer :: dmft_dc,iatom,iatomc,ierr,im,ispinor,isppol
387 : integer :: itypat,lpawu,natom,ndim,nspinor,nsppol
388 : real(dp) :: dc,jpawu,ntot,upawu
389 : logical :: amf,fll,nmdc
390 : character(len=500) :: message
391 199 : complex(dp), allocatable :: occ(:,:),vdc(:,:)
392 : ! *********************************************************************
393 :
394 199 : dmft_dc = paw_dmft%dmft_dc
395 199 : natom = paw_dmft%natom
396 199 : nspinor = paw_dmft%nspinor
397 199 : nsppol = paw_dmft%nsppol
398 :
399 199 : amf = (dmft_dc == 2 .or. dmft_dc == 6) ! AMF double counting
400 199 : fll = (dmft_dc == 1 .or. dmft_dc == 4 .or. dmft_dc == 5) ! FLL double counting
401 199 : nmdc = (nspinor == 2 .or. (dmft_dc >= 5 .and. dmft_dc <= 8)) ! non-magnetic double counting
402 :
403 199 : if ((.not. fll) .and. (.not. amf) .and. dmft_dc /= 7 .and. dmft_dc /= 8) &
404 0 : & ABI_ERROR("not implemented")
405 199 : if (amf .and. (nspinor == 2)) then
406 0 : write(message,'(a,i4,i4,2x,e20.10)') " AMF Double counting is under test for SOC"
407 0 : ABI_WARNING(message)
408 : end if
409 :
410 199 : if (dmft_dc == 8) then
411 0 : paw_dmft%edc(:) = zero
412 0 : paw_dmft%edcdc(:) = zero
413 0 : call zero_matlu(hdc(:),natom)
414 : end if
415 :
416 : iatomc = -1
417 :
418 890 : do iatom=1,natom
419 691 : lpawu = paw_dmft%lpawu(iatom)
420 691 : if (lpawu == -1) cycle
421 235 : iatomc = iatomc + 1
422 10896 : hdc(iatom)%mat(:,:,:) = czero
423 235 : ntot = charge_loc(nsppol+1,iatom)
424 235 : ndim = 2*lpawu + 1
425 235 : itypat = paw_dmft%typat(iatom)
426 235 : upawu = hu(itypat)%upawu
427 235 : jpawu = merge(zero,hu(itypat)%jpawu,dmft_dc==4)
428 :
429 434 : if (dmft_dc == 8) then
430 0 : if (mod(iatomc,paw_dmft%nproc) /= paw_dmft%myproc) cycle
431 0 : ABI_MALLOC(occ,(ndim,ndim))
432 0 : ABI_MALLOC(vdc,(ndim,ndim))
433 0 : occ(:,:) = czero
434 0 : do isppol=1,nsppol
435 0 : do ispinor=1,nspinor
436 0 : occ(:,:) = occ(:,:) + occ_matlu(iatom)%mat(1+(ispinor-1)*ndim:ispinor*ndim,1+(ispinor-1)*ndim:ispinor*ndim,isppol)
437 : end do ! ispinor
438 : end do ! isppol
439 0 : if (nsppol == 1 .and. nspinor == 1) occ(:,:) = occ(:,:) * two
440 : call compute_exactDC(lpawu,pawtab(itypat),paw_dmft%radgrid(itypat),occ(:,:), &
441 0 : & vdc(:,:),paw_dmft%edc(iatom),paw_dmft%edcdc(iatom),paw_dmft%ixc)
442 0 : do isppol=1,nsppol
443 0 : do ispinor=1,nspinor
444 0 : hdc(iatom)%mat(1+(ispinor-1)*ndim:ispinor*ndim,1+(ispinor-1)*ndim:ispinor*ndim,isppol) = vdc(:,:)
445 : end do ! ispinor
446 : end do ! isppol
447 0 : ABI_FREE(occ)
448 0 : ABI_FREE(vdc)
449 : else
450 235 : if (nmdc) then ! Non-magnetic DC
451 29 : if (fll) dc = upawu*(ntot-half) - half*jpawu*(ntot-one)
452 29 : if (amf) dc = upawu*ntot*half + (upawu-jpawu)*ntot*half*dble(2*lpawu)/dble(2*lpawu+1)
453 29 : if (dmft_dc == 7) dc = upawu*(dble(paw_dmft%dmft_nominal(iatom))-half) - &
454 0 : & half*jpawu*(dble(paw_dmft%dmft_nominal(iatom))-one)
455 : end if ! nmdc
456 235 : ndim = nspinor * ndim
457 594 : do isppol=1,nsppol
458 359 : if (.not. nmdc) then ! Magnetic DC
459 304 : if (fll) dc = upawu*(ntot-half) - jpawu*(charge_loc(isppol,iatom)-half)
460 304 : if (amf) dc = upawu*charge_loc(min(3-isppol,nsppol),iatom) + &
461 0 : & (upawu-jpawu)*charge_loc(isppol,iatom)*dble(2*lpawu)/dble(2*lpawu+1)
462 : end if ! not nmdc
463 2728 : do im=1,ndim
464 2037 : hdc(iatom)%mat(im,im,isppol) = dc
465 : end do ! im
466 : end do ! isppol
467 : end if ! dc=8
468 : end do ! iatom
469 :
470 199 : if (dmft_dc == 8) then
471 0 : call xmpi_sum(paw_dmft%edc(:),paw_dmft%spacecomm,ierr)
472 0 : call xmpi_sum(paw_dmft%edcdc(:),paw_dmft%spacecomm,ierr)
473 0 : call xmpi_matlu(hdc(:),natom,paw_dmft%spacecomm)
474 : end if ! dmft_dc=8
475 :
476 199 : end subroutine dc_self
477 : !!***
478 :
479 : !!****f* m_self/rw_self
480 : !! NAME
481 : !! rw_self
482 : !!
483 : !! FUNCTION
484 : !! Read/write self-energy on file.
485 : !!
486 : !! INPUTS
487 : !! self <type(self_type)>= variables related to self-energy
488 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
489 : !! prtopt = flag for printing
490 : !! opt_rw = 0 (default) Set self-energy either to double counting or 0 depending on dmft_rslf
491 : !! 1 Read Self-Energy.
492 : !! 2 Write Self-Energy.
493 : !! 3 Impose Self-Energy.
494 : !! istep_iter = iteration step
495 : !! opt_char = char to add at the end of the filename
496 : !! opt_imagonly = only reads the imaginary part (useful when reading from Maxent)
497 : !! opt_selflimit = 0th order moment of the self-energy (useful when reading from Maxent)
498 : !! opt_hdc = double counting (useful when reading from Maxent)
499 : !! opt_stop = stop when encountering errors
500 : !! opt_maxent = if > 0, read/write Maxent files
501 : !!
502 : !! OUTPUT
503 : !!
504 : !! SOURCE
505 :
506 406 : subroutine rw_self(self,paw_dmft,prtopt,opt_rw,istep_iter,opt_char,opt_imagonly,&
507 406 : & opt_selflimit,opt_hdc,opt_stop,opt_maxent)
508 :
509 : !Arguments ------------------------------------
510 : !type
511 : type(self_type), intent(inout) :: self
512 : type(paw_dmft_type), intent(inout) :: paw_dmft
513 : integer, intent(in) :: prtopt
514 : integer, optional, intent(in) :: istep_iter,opt_imagonly,opt_maxent,opt_rw,opt_stop
515 : character(len=4), optional, intent(in) :: opt_char
516 : type(matlu_type), optional, intent(inout) :: opt_selflimit(paw_dmft%natom)
517 : type(matlu_type), optional, intent(in) :: opt_hdc(paw_dmft%natom)
518 : !local variables-------------------------------
519 : integer :: i,iall,iatom,iatu,icount,ier,iexist2,iexit,iflavor,ifreq,im,im1,ioerr
520 : integer :: ispinor,ispinor1,isppol,istep,istep_imp,istepiter,iter,iter_imp,lpawu,master
521 : integer :: myproc,natom,natom_read,ncount,ndim,ndim_read,nrecl,nspinor,nspinor_read
522 : integer :: nsppol,nsppol_read,nw_read,optmaxent,optrw,readimagonly,spacecomm,unitrot
523 : real(dp) :: fermie_read,x_r,x_i,xtemp
524 : logical :: lexist,lexist_rot,nondiaglevels,prtself
525 : character(len=30000) :: message ! Big buffer to avoid buffer overflow.
526 : character(len=fnlen) :: stringfile,tmpfil,tmpfil2,tmpfilrot,tmpmatrot
527 : character(len=1) :: tag_is
528 : character(len=3) :: self_iter
529 : character(len=4) :: chtemp
530 : character(len=5) :: tag_freq
531 : character(len=10) :: tag_at,tag_iflavor
532 : character(len=13) :: tag
533 : character(len=50) :: string_format
534 406 : type(oper_type) :: energy_level
535 406 : integer, allocatable :: unitselffunc_arr(:),unitselffunc_arr2(:),unitselfrot(:,:,:,:)
536 406 : real(dp), allocatable :: s_i(:,:),s_r(:,:) !,fermie_read2(:)
537 406 : complex(dp), allocatable :: buffer(:)
538 406 : type(matlu_type), allocatable :: eigvectmatlu(:),level_diag(:),selfmomrot(:,:)
539 406 : type(oper_type), allocatable :: selfrotmatlu(:)
540 : ! *********************************************************************
541 :
542 : ABI_NVTX_START_RANGE(NVTX_DMFT_RW_SELF)
543 406 : natom = paw_dmft%natom
544 406 : nspinor = paw_dmft%nspinor
545 406 : nsppol = paw_dmft%nsppol
546 : !mbandc = paw_dmft%mbandc
547 : !nkpt = paw_dmft%nkpt
548 :
549 : ! Initialise spaceComm, myproc, and nproc
550 406 : istep = 0
551 406 : iter = 0
552 406 : istep_imp = 0
553 406 : istepiter = 0
554 406 : iter_imp = 0
555 406 : optmaxent = 0
556 406 : optrw = 0
557 406 : prtself = (paw_dmft%dmft_prtself == 1)
558 406 : readimagonly = 0
559 :
560 406 : if (present(opt_rw)) optrw = opt_rw
561 406 : if (present(opt_maxent)) optmaxent = opt_maxent
562 :
563 406 : if (present(opt_imagonly)) then
564 4 : if (opt_imagonly == 1 .and. paw_dmft%dmft_solv >= 5) then
565 4 : readimagonly = opt_imagonly
566 4 : write(message,*)
567 4 : write(message,'(a,4x,2a)') ch10,"About to read imaginary part of Self energy"
568 4 : call wrtout(std_out,message,'COLL')
569 : else
570 0 : write(message,'(4x,2a)') "About to read both real and imaginary part of Self energy"
571 0 : call wrtout(std_out,message,'COLL')
572 : end if ! opt_imagonly
573 : end if ! present(opt_imagonly)
574 :
575 406 : if (present(istep_iter)) istepiter = istep_iter
576 :
577 406 : if (paw_dmft%use_fixed_self > 0) then
578 0 : istep = istepiter / 1000
579 0 : iter = istepiter - (istepiter/1000)*1000
580 0 : istep_imp = paw_dmft%use_fixed_self / 1000
581 0 : iter_imp = paw_dmft%use_fixed_self - (paw_dmft%use_fixed_self/1000)*1000
582 : end if !use_fixed_self
583 :
584 406 : if (paw_dmft%dmft_rslf <= 0 .and. optrw == 1) optrw = 0
585 :
586 406 : iexist2 = 1
587 406 : iexit = 0
588 406 : ioerr = 0
589 406 : lexist = .true.
590 406 : master = 0
591 406 : myproc = paw_dmft%myproc
592 406 : spacecomm = paw_dmft%spacecomm
593 : !nproc = paw_dmft%nproc
594 :
595 : ! write(std_out,*) "myproc,master",myproc,master
596 : !if(prtopt>200) then
597 : !endif
598 :
599 : ! - For the Tentative rotation of the self-energy file (begin init)
600 406 : if (optmaxent > 0) then
601 110 : if (optrw == 2) then
602 106 : write(message,'(a,2x,a)') ch10," == About to print self-energy for MAXENT code in basis that diagonalizes the atomic levels"
603 4 : else if (optrw == 1) then
604 4 : write(message,'(a,2x,a)') ch10," == About to read self-energy from MAXENT code"
605 : end if
606 110 : call wrtout(std_out,message,'COLL')
607 :
608 704 : ABI_MALLOC(eigvectmatlu,(natom))
609 110 : call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu(:),eigvectmatlu(:))
610 : end if ! optmaxent > 0
611 : ! - For the Tentative rotation of the self-energy file (end init)
612 :
613 : ! - For the Tentative rotation of the self-energy file (begin diag)
614 406 : if (optrw == 2 .and. optmaxent > 0) then
615 636 : ABI_MALLOC(unitselfrot,(2*paw_dmft%maxlpawu+1,nspinor,nsppol,natom)) ! 7 is the max ndim possible
616 688 : ABI_MALLOC(level_diag,(natom))
617 23407 : ABI_MALLOC(selfrotmatlu,(self%nw))
618 106 : call init_oper(paw_dmft,energy_level,opt_ksloc=2)
619 106 : call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu(:),level_diag(:))
620 106 : call compute_levels(energy_level,self%hdc,paw_dmft,nondiag=nondiaglevels)
621 106 : write(tag,'(f13.5)') paw_dmft%fermie
622 106 : if (self%nw >= 2) then
623 106 : write(message,'(a,2x,2a)') ch10," == Print non Diagonalized Self Energy for Fermi Level= ",adjustl(tag)
624 106 : call wrtout(std_out,message,'COLL')
625 106 : call print_matlu(self%oper(2)%matlu(:),natom,1,compl=1,opt_exp=1)
626 : end if
627 106 : call diag_matlu(energy_level%matlu(:),level_diag(:),natom,prtopt,eigvectmatlu(:),test=paw_dmft%dmft_solv)
628 106 : write(message,'(a,2x,2a)') ch10," == Print Diagonalized levels for Fermi Level= ",adjustl(tag)
629 106 : call wrtout(std_out,message,'COLL')
630 106 : call print_matlu(level_diag(:),natom,1,compl=1,opt_exp=1)
631 : ! Rotate self
632 23195 : do ifreq=1,self%nw
633 23089 : call init_oper(paw_dmft,selfrotmatlu(ifreq),opt_ksloc=2)
634 23089 : if (self%distrib%procf(ifreq) /= myproc) cycle
635 8032 : call copy_matlu(self%oper(ifreq)%matlu(:),selfrotmatlu(ifreq)%matlu(:),natom)
636 23195 : call rotate_matlu(selfrotmatlu(ifreq)%matlu(:),eigvectmatlu(:),natom,1)
637 : end do ! ifreq
638 106 : call gather_oper(selfrotmatlu(:),self%distrib,paw_dmft,2,master=master)
639 106 : if (self%has_moments == 1) then
640 0 : ABI_MALLOC(selfmomrot,(natom,self%nmoments))
641 0 : do i=1,self%nmoments
642 0 : call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu(:),selfmomrot(:,i))
643 0 : call copy_matlu(self%moments(i)%matlu(:),selfmomrot(:,i),natom)
644 0 : call rotate_matlu(selfmomrot(:,i),eigvectmatlu(:),natom,1)
645 : end do ! i
646 : end if
647 530 : do ifreq=1,min(3,self%nw)
648 318 : write(tag_freq,'(i5)') ifreq
649 424 : if (ifreq < 3 .and. myproc == master) then ! very important to call print_matlu only on master node
650 98 : write(message,'(a,2x,2a)') ch10," == Print non Rotated Self Energy for frequency ",adjustl(tag_freq)
651 98 : call wrtout(std_out,message,'COLL')
652 98 : call print_matlu(self%oper(ifreq)%matlu(:),natom,1,compl=1)
653 98 : write(message,'(a,2x,2a)') ch10," == Print Rotated Self Energy for frequency ",adjustl(tag_freq)
654 98 : call wrtout(std_out,message,'COLL')
655 98 : call print_matlu(selfrotmatlu(ifreq)%matlu(:),natom,1,compl=1)
656 220 : else if (ifreq == 3) then
657 106 : write(message,'(a,2x,a,i4)') ch10," (Other frequencies not printed)"
658 106 : call wrtout(std_out,message,'COLL')
659 : end if ! ifreq<3
660 : end do ! ifreq
661 : end if ! optrw=2 and optmaxent>0
662 : ! Create file for rotation
663 406 : if (optmaxent > 0 .and. myproc == master .and. (optrw == 1 .or. optrw == 2)) then
664 208 : do iatom=1,natom
665 158 : lpawu = paw_dmft%lpawu(iatom)
666 158 : if (lpawu == -1) cycle
667 62 : call int2char4(iatom,tag_at)
668 62 : ABI_CHECK((tag_at(1:1)/='#'),'Bug: string length too short!')
669 62 : if (optrw == 2) then
670 61 : tmpmatrot = trim(paw_dmft%filapp)//'_UnitaryMatrix_iatom'//trim(tag_at)
671 : else if (optrw == 1) then
672 1 : tmpmatrot = trim(paw_dmft%filnamei)//'_UnitaryMatrix_iatom'//trim(tag_at)
673 1 : inquire(file=trim(tmpmatrot),exist=lexist_rot)
674 1 : if (.not. lexist_rot) ABI_ERROR("File "//trim(tmpmatrot)//" does not exist !")
675 : end if ! optrw
676 62 : unitrot = 3189 + iatom
677 : #ifdef FC_NAG
678 : open(unit=unitrot,file=trim(tmpmatrot),status='unknown',form='formatted',recl=ABI_RECL)
679 : #else
680 62 : open(unit=unitrot,file=trim(tmpmatrot),status='unknown',form='formatted')
681 : #endif
682 62 : write(std_out,'(2a)') " Open file ",trim(tmpmatrot)
683 62 : rewind(unitrot)
684 62 : ndim = (2*lpawu+1) * nspinor
685 165 : do isppol=1,nsppol
686 674 : do im=1,ndim
687 3433 : do im1=1,ndim
688 3330 : if (optrw == 2) then
689 2796 : write(message,*) dble(eigvectmatlu(iatom)%mat(im,im1,isppol)),aimag(eigvectmatlu(iatom)%mat(im,im1,isppol))
690 2796 : call wrtout(unitrot,message,'COLL')
691 : else if (optrw == 1) then
692 25 : read(unitrot,*) x_r,x_i
693 25 : eigvectmatlu(iatom)%mat(im,im1,isppol) = cmplx(x_r,x_i,kind=dp)
694 : end if ! optrw
695 : end do ! im1
696 : end do ! im
697 : end do ! isppol
698 208 : close(unitrot)
699 : end do ! iatom
700 :
701 50 : if (optrw == 1) then
702 1 : write(message,'(2a)') ch10," == Print non-rotated high-frequency limit of the self-energy"
703 1 : call wrtout(std_out,message,'COLL')
704 1 : call print_matlu(opt_selflimit(:),natom,1,compl=1)
705 1 : call rotate_matlu(opt_selflimit(:),eigvectmatlu(:),natom,1)
706 1 : write(message,'(2a)') ch10," == Print rotated high-frequency limit of the self-energy"
707 1 : call wrtout(std_out,message,'COLL')
708 1 : call print_matlu(opt_selflimit(:),natom,1,compl=1)
709 : end if ! optrw=1
710 :
711 : end if ! optmaxent > 0
712 : ! - For the Tentative rotation of the self-energy file (end diag)
713 :
714 406 : if ((optrw == 2 .or. optrw == 1) .and. myproc == master) then
715 390 : ABI_MALLOC(unitselffunc_arr,(natom*nsppol))
716 130 : if (optrw == 2) then
717 182 : ABI_MALLOC(unitselffunc_arr2,(natom*nsppol))
718 91 : if (paw_dmft%idmftloop < 10) then
719 89 : write(self_iter,'("00",i1)') paw_dmft%idmftloop
720 2 : else if (paw_dmft%idmftloop >= 10 .and. paw_dmft%idmftloop < 100) then
721 2 : write(self_iter,'("0",i2)') paw_dmft%idmftloop
722 0 : else if (paw_dmft%idmftloop >= 100 .and. paw_dmft%idmftloop < 1000) then
723 0 : write(self_iter,'(i3)') paw_dmft%idmftloop
724 : else
725 0 : self_iter="xxx"
726 : end if ! idmftloop
727 : end if ! optrw=2
728 130 : iall = 0
729 532 : do iatom=1,natom
730 402 : lpawu = paw_dmft%lpawu(iatom)
731 402 : if (lpawu == -1) cycle
732 150 : ndim = 2*lpawu + 1
733 600 : ABI_MALLOC(s_r,(ndim*nspinor,ndim*nspinor))
734 450 : ABI_MALLOC(s_i,(ndim*nspinor,ndim*nspinor))
735 : ! write(std_out,*) "print_self",ndim
736 150 : call int2char4(iatom,tag_at)
737 150 : ABI_CHECK((tag_at(1:1)/='#'),'Bug: string length too short!')
738 390 : do isppol=1,nsppol
739 240 : write(tag_is,'(i1)') isppol
740 : ! do ispinor=1,nspinor
741 240 : iall = iall + 1
742 : ! write(tag_is2,'(i1)')ispinor
743 :
744 : ! ===========================
745 : ! == Create name for file
746 : ! ===========================
747 :
748 240 : if (self%w_type == "real") then
749 1 : tmpfil = trim(merge(paw_dmft%filnamei,paw_dmft%filapp,optrw==1))//'_Self_ra-omega_iatom'//trim(tag_at)//'_isppol'//tag_is
750 : else
751 239 : if (present(opt_char)) then
752 0 : tmpfil = trim(paw_dmft%filapp)//'Self_ra-omega_iatom'//trim(tag_at)//'_isppol'//tag_is//opt_char
753 : else
754 239 : stringfile = "_iatom" // trim(tag_at) // '_isppol' // tag_is
755 239 : if (optrw == 1 .and. paw_dmft%idmftloop == 0) then
756 32 : tmpfil = trim(paw_dmft%filselfin) // stringfile
757 32 : iexist2 = paw_dmft%ireadself
758 : else
759 207 : tmpfil = trim(paw_dmft%filapp) // '_Self-omega' // stringfile
760 : end if
761 : end if ! opt_char
762 : end if ! w_type
763 240 : if (optrw == 2) tmpfil2 = trim(tmpfil)//"_"//self_iter
764 :
765 240 : if (optrw == 1 .and. iexist2 == 1) then
766 30 : write(message,'(3a)') ch10," == Read self-energy and Fermi Level from file ",trim(tmpfil)
767 30 : call wrtout(std_out,message,'COLL')
768 210 : else if (optrw == 2) then
769 179 : write(message,'(3a)') ch10," == Write self-energy and Fermi Level on file ",trim(tmpfil)
770 179 : call wrtout(std_out,message,'COLL')
771 : end if
772 : !unitselffunc_arr(iall)=300+iall-1
773 240 : unitselffunc_arr(iall) = get_unit()
774 240 : ABI_CHECK(unitselffunc_arr(iall) > 0, "Cannot find free IO unit!")
775 :
776 : !- For the Tentative rotation of the self-energy file (create file)
777 240 : if (optrw == 2 .and. optmaxent > 0) then
778 102 : iflavor = 0
779 206 : do ispinor=1,nspinor
780 710 : do im=1,ndim
781 504 : iflavor = iflavor + 1
782 504 : call int2char4(iflavor,tag_iflavor)
783 504 : unitselfrot(im,ispinor,isppol,iatom) = 3000 + iflavor
784 : !ABI_CHECK(unitselfrot(im,ispinor,isppol,iatom) > 0, "Cannot find free IO unit for unitselfrot!")
785 : tmpfilrot = trim(paw_dmft%filapp)//'_Selfmxent'//&
786 504 : & trim(tag_at)//'_is'//tag_is//'_iflav'//trim(tag_iflavor)
787 504 : write(std_out,*) "Create file ",trim(tmpfilrot)," unit ",unitselfrot(im,ispinor,isppol,iatom)," for flavor",iflavor
788 : #ifdef FC_NAG
789 : open(unit=unitselfrot(im,ispinor,isppol,iatom),file=trim(tmpfilrot),status='unknown',form='formatted',recl=ABI_RECL)
790 : #else
791 504 : open(unit=unitselfrot(im,ispinor,isppol,iatom),file=trim(tmpfilrot),status='unknown',form='formatted')
792 : #endif
793 504 : rewind(unitselfrot(im,ispinor,isppol,iatom))
794 504 : write(unitselfrot(im,ispinor,isppol,iatom),'(3a)') "# Diagonal component of the self-energy, in the basis that diagonalizes the electronic levels.", &
795 1112 : & ch10,"# Frequency (Ha) Real part Imaginary part"
796 :
797 : end do ! im
798 : end do ! ispinor
799 : end if ! optrw=2 and optmaxent>0
800 : !- For the Tentative rotation of the self-energy file (create file)
801 :
802 : ! write(std_out,*) "1"
803 :
804 : ! ===========================
805 : ! == Read: check that the file exists
806 : ! ===========================
807 240 : if (optrw == 1) then
808 : ! write(std_out,*) "3"
809 61 : inquire(file=trim(tmpfil),exist=lexist,recl=nrecl)
810 61 : if (.not. lexist .and. (paw_dmft%ireadself == 1 .or. paw_dmft%idmftloop > 0)) then
811 : ! write(std_out,*) "4"
812 0 : iexist2 = 0
813 0 : write(message,'(4x,a,i5,3a)') "File number",unitselffunc_arr(iall),&
814 0 : & " called ",trim(tmpfil)," does not exist"
815 : ! write(std_out,*) lexist,nrecl
816 0 : call wrtout(std_out,message,'COLL')
817 : end if ! not lexist
818 : end if ! optrw=1
819 : !write(std_out,*) "2"
820 :
821 : ! ===========================
822 : ! == Open file
823 : ! ===========================
824 240 : if (optrw == 2 .or. (optrw == 1 .and. iexist2 == 1)) then
825 : !write(std_out,*) "5"
826 : #ifdef FC_NAG
827 : open(unit=unitselffunc_arr(iall),file=trim(tmpfil),status='unknown',form='formatted',recl=ABI_RECL)
828 : #else
829 209 : open(unit=unitselffunc_arr(iall),file=trim(tmpfil),status='unknown',form='formatted')
830 : #endif
831 209 : rewind(unitselffunc_arr(iall))
832 :
833 209 : if (optrw == 2 .and. prtself) then
834 0 : unitselffunc_arr2(iall) = get_unit()
835 0 : ABI_CHECK(unitselffunc_arr2(iall) > 0,"Cannot find free IO unit!")
836 : #ifdef FC_NAG
837 : open(unit=unitselffunc_arr2(iall),file=trim(tmpfil2),status="unknown",form="formatted",recl=ABI_RECL)
838 : #else
839 0 : open(unit=unitselffunc_arr2(iall),file=trim(tmpfil2),status="unknown",form="formatted")
840 : #endif
841 0 : rewind(unitselffunc_arr2(iall))
842 : end if ! optrw=2
843 :
844 : !write(std_out,*) "61",nrecl
845 209 : if (prtopt >= 3) then
846 2 : write(message,'(3a,i4)') ' Opened file : ',trim(tmpfil),' on unit ',unitselffunc_arr(iall)
847 2 : call wrtout(std_out,message,'COLL')
848 : end if ! prtopt>=3
849 : end if
850 : !write(std_out,*) "6",nrecl
851 :
852 : ! ===========================
853 : ! == Check Header
854 : ! ===========================
855 :
856 240 : if (optrw == 2) then
857 :
858 179 : write(message,'(11a,4i5,i6,2x,e25.17e3)') "# DFT+DMFT self-energy for each frequency",ch10, &
859 179 : & "# Columns are ordered this way:",ch10, &
860 179 : & "# Frequency (Ha) ((((Re(Sigma(im,im1,ispinor,ispinor1)) ", &
861 179 : & "Im(Sigma(im,im1,ispinor,ispinor1)),im=1,2*l+1),im1=1,2*l+1),", &
862 179 : & "ispinor=1,nspinor),ispinor1=1,nspinor) where the leftmost index varies first", &
863 179 : & ch10,"# natom,nsppol,nspinor,ndim,nw,fermilevel",ch10,&
864 358 : & "####",natom,nsppol,nspinor,ndim,self%nw,paw_dmft%fermie
865 179 : call wrtout(unitselffunc_arr(iall),message,'COLL')
866 179 : if (prtself) then
867 0 : call wrtout(unitselffunc_arr2(iall),message,'COLL')
868 : end if
869 61 : else if (optrw == 1 .and. iexist2 == 1 .and. readimagonly == 0) then
870 29 : read(unitselffunc_arr(iall),*)
871 29 : read(unitselffunc_arr(iall),*)
872 29 : read(unitselffunc_arr(iall),*)
873 29 : read(unitselffunc_arr(iall),*)
874 : read(unitselffunc_arr(iall),*,iostat=ioerr) &
875 29 : & chtemp,natom_read,nsppol_read,nspinor_read,ndim_read,nw_read,fermie_read
876 : !if(ioerr<0) then
877 : ! write(std_out,*)" HEADER IOERR"
878 : ! write(std_out,'(a4,2x,31(e15.8,2x))') chtemp,natom_read,nsppol_read,nspinor_read,ndim_read,nw_read,fermie_read
879 : !endif
880 29 : if (ioerr == 0) then
881 29 : write(message,'(a,3x,3a,i12,2a,i11,2a,i10,2a,i13,2a,i15,2a,e25.8)') ch10,"Data in Self Energy file corresponds to",&
882 29 : & ch10," natom",natom_read,&
883 29 : & ch10," nsppol",nsppol_read,&
884 29 : & ch10," nspinor",nspinor_read,&
885 29 : & ch10," ndim",ndim_read, &
886 29 : & ch10," nw",nw_read, &
887 58 : & ch10," Fermi level",fermie_read
888 29 : call wrtout(std_out,message,'COLL')
889 : if ((natom /= natom_read) .or. (nsppol_read /= nsppol) .or. &
890 29 : & (nspinor /= nspinor_read) .or. (nw_read /= self%nw)) then
891 0 : write(message,'(a,3x,3a,i12,2a,i11,2a,i10,2a,i13,2a,i15,2a,e25.8)') ch10,"Data required is ",&
892 0 : & ch10," natom",natom,&
893 0 : & ch10," nsppol",nsppol,&
894 0 : & ch10," nspinor",nspinor,&
895 0 : & ch10," ndim",ndim, &
896 0 : & ch10," nw",self%nw, &
897 0 : & ch10," Fermi level",paw_dmft%fermie
898 0 : call wrtout(std_out,message,'COLL')
899 0 : message = "Dimensions in self-energy file are not correct"
900 0 : if (present(opt_stop)) then
901 0 : ABI_ERROR(message)
902 : else
903 0 : ABI_WARNING(message)
904 : end if
905 0 : iexist2 = 2
906 : end if
907 : else
908 0 : ABI_WARNING("Self-energy file is empty")
909 : end if ! ioerr
910 : end if ! optrw
911 : !write(std_out,*) "7"
912 :
913 : ! ===========================
914 : ! == Write/Read self in the file
915 : ! ===========================
916 :
917 : !rewind(111)
918 46808 : do ifreq=1,self%nw
919 46808 : if (optrw == 2) then
920 : ! write(std_out,'(a,2x,31(e15.8,2x))') &
921 : !& "SETEST",self%omega(ifreq),&
922 : !& (self%oper(ifreq)%matlu(iatom)%mat(im,im,isppol,ispinor,ispinor)&
923 : !& ,im=1,ndim)
924 : ! write(std_out,*) self%omega(ifreq),&
925 : !& ((self%oper(ifreq)%matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor)&
926 : !& ,im=1,ndim),im1=1,ndim)
927 : ! write(message,'(2x,393(e25.17,2x))') self%omega(ifreq),&
928 : !& ((((self%oper(ifreq)%matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1)&
929 : !& ,im=1,ndim),im1=1,ndim),ispinor=1,nspinor),ispinor1=1,nspinor)
930 :
931 : !MGNAG Runtime Error: wrtout_cpp.f90, line 896: Buffer overflow on output
932 : !Is it possible to rewrite the code below to avoid such a long message
933 : !What about Netcdf binary files ?
934 28956 : string_format = merge('(2x,393(e25.17e3,2x))','(2x,393(e18.10e3,2x))',nspinor==1)
935 28956 : write(message,string_format) self%omega(ifreq),&
936 746272 : & ((((dble(self%oper(ifreq)%matlu(iatom)%mat(im+(ispinor-1)*ndim,im1+(ispinor1-1)*ndim,isppol)),&
937 892640 : & aimag(self%oper(ifreq)%matlu(iatom)%mat(im+(ispinor-1)*ndim,im1+(ispinor1-1)*ndim,isppol)),&
938 1009664 : & im=1,ndim),im1=1,ndim),ispinor=1,nspinor),ispinor1=1,nspinor)
939 28956 : call wrtout(unitselffunc_arr(iall),message,'COLL')
940 28956 : if (prtself) then
941 0 : call wrtout(unitselffunc_arr2(iall),message,'COLL')
942 : end if
943 :
944 : !- For the Tentative rotation of the self-energy file (begin rot)
945 : !----------------------------------------------------------------
946 28956 : if (optmaxent > 0) then
947 : !iflavor = 0
948 30024 : do ispinor=1,nspinor
949 103728 : do im=1,ndim
950 : !iflavor = iflavor + 1
951 : ! if(ifreq<5) then
952 : ! write(std_out,*) "Write in file unit",unitselfrot(iatom,isppol,ispinor,im),"for flavor",iflavor
953 : ! endif
954 73704 : write(message,'(2x,393(es24.16e3,2x))') self%omega(ifreq),&
955 73704 : & dble(selfrotmatlu(ifreq)%matlu(iatom)%mat(im+(ispinor-1)*ndim,im+(ispinor-1)*ndim,isppol)),&
956 147408 : & aimag(selfrotmatlu(ifreq)%matlu(iatom)%mat(im+(ispinor-1)*ndim,im+(ispinor-1)*ndim,isppol))
957 : ! write(6,'(2x,393(e18.10,2x))') self%omega(ifreq),&
958 : !& real(selfrotmatlu(iatom)%mat(im,im,isppol,ispinor,ispinor)),&
959 : !& aimag(selfrotmatlu(iatom)%mat(im,im,isppol,ispinor,ispinor))
960 : ! if(iflavor==1) then
961 : ! write(1024,*) iatom,isppol,ispinor,im,unitselfrot(iatom,isppol,ispinor,im)
962 : ! write(1024,'(2x,393(e18.10,2x))') self%omega(ifreq),&
963 : ! & real(selfrotmatlu(iatom)%mat(im,im,isppol,ispinor,ispinor)),&
964 : ! & aimag(selfrotmatlu(iatom)%mat(im,im,isppol,ispinor,ispinor))
965 : ! endif
966 88816 : call wrtout(unitselfrot(im,ispinor,isppol,iatom),message,'COLL')
967 : end do ! im
968 : end do ! ispinor
969 : end if ! optrw=2 and optmaxent>0
970 : !- For the Tentative rotation of the self-energy file (end rot)
971 :
972 : ! write(std_out,*) unitselffunc_arr(iall)
973 17612 : else if (optrw == 1 .and. iexist2 == 1 .and. ioerr == 0 .and. readimagonly == 0) then
974 : !write(std_out,*) "8"
975 : ! read(unitselffunc_arr(iall),'(2x,31(e15.8,2x))',iostat=ioerr) &
976 : !& xtemp,(s_r(im),s_i(im),im=1,ndim)
977 : !if (readimagonly == 0) then
978 3462 : read(unitselffunc_arr(iall),*,iostat=ioerr) xtemp,&
979 88928 : & ((((s_r(im+(ispinor-1)*ndim,im1+(ispinor1-1)*ndim),s_i(im+(ispinor-1)*ndim,im1+(ispinor1-1)*ndim), &
980 102776 : & im=1,ndim),im1=1,ndim),ispinor=1,nspinor),ispinor1=1,nspinor)
981 : ! if(ioerr<0) then
982 : ! write(std_out,*)" SELF IOERR<"
983 : ! else if(ioerr>0) then
984 : ! write(std_out,*)" SELF IOERR>"
985 : ! write(std_out,'(a4,2x,31(e15.8,2x))') xtemp,(s_r(im),s_i(im),im=1,ndim)
986 : ! endif
987 92390 : self%oper(ifreq)%matlu(iatom)%mat(:,:,isppol) = cmplx(s_r(:,:),s_i(:,:),kind=dp)
988 : end if ! optrw
989 : end do ! ifreq
990 :
991 : !- For the Tentative rotation of the self-energy file (begin close file)
992 240 : if (optrw == 2 .and. optmaxent > 0) then
993 206 : do ispinor=1,nspinor
994 710 : do im=1,ndim
995 504 : if (self%has_moments == 1) then
996 0 : do i=1,self%nmoments
997 0 : write(tag_is,'(i1)') i
998 0 : write(message,'(a,2x,2(es24.16e3,2x))') "#moments_"//tag_is, &
999 0 : & dble(selfmomrot(iatom,i)%mat(im+(ispinor-1)*ndim,im+(ispinor-1)*ndim,isppol)), &
1000 0 : & aimag(selfmomrot(iatom,i)%mat(im+(ispinor-1)*ndim,im+(ispinor-1)*ndim,isppol))
1001 0 : call wrtout(unitselfrot(im,ispinor,isppol,iatom),message,'COLL')
1002 : end do ! i
1003 : end if ! nmoments
1004 504 : close(unitselfrot(im,ispinor,isppol,iatom))
1005 608 : write(std_out,*) "Close file unit",unitselfrot(im,ispinor,isppol,iatom)
1006 : end do ! im
1007 : end do ! ispinor
1008 : end if ! optrw=2 and optmaxent>0
1009 :
1010 : !- For the Tentative rotation of the self-energy file (end close file)
1011 :
1012 240 : if (optrw == 1 .and. iexist2 == 1 .and. ioerr == 0 .and. readimagonly == 1) then
1013 :
1014 : ! Read self energy from Maxent (imag part) on the real axis
1015 : !----------------------------------------------------------
1016 3001 : do ifreq=1,self%nw
1017 96001 : self%oper(ifreq)%matlu(iatom)%mat(:,:,:) = czero
1018 : end do ! ifreq
1019 2 : do ispinor=1,nspinor
1020 7 : do im=1,ndim
1021 15006 : do ifreq=1,self%nw
1022 15000 : read(unitselffunc_arr(iall),*,iostat=ioerr) xtemp,s_i(im+(ispinor-1)*ndim,im+(ispinor-1)*ndim)
1023 : ! minus sign because - Im Sigma is the output of OmegaMaxent
1024 : self%oper(ifreq)%matlu(iatom)%mat(im+(ispinor-1)*ndim,im+(ispinor-1)*ndim,isppol) &
1025 15005 : & = cmplx(zero,-half*s_i(im+(ispinor-1)*ndim,im+(ispinor-1)*ndim),kind=dp)
1026 : end do ! ifreq
1027 : end do ! im
1028 : end do ! ispinor
1029 :
1030 1 : write(message,'(4x,2a)') " Read only diagonal self energy from Maxent"
1031 1 : call wrtout(std_out,message,'COLL')
1032 : !write(6,*) "opt_hdc",opt_hdc(1)%mat(1,1,1,1,1)
1033 :
1034 : end if ! optrw=1
1035 :
1036 : ! ===========================
1037 : ! == Write/Read hdc in the file
1038 : ! ===========================
1039 240 : if (optrw == 2) then
1040 : ! write(std_out,'(a,2x,31(e15.8,2x))') &
1041 : !& "SETEST #dc ",(self%hdc%matlu(iatom)%mat(im,im,isppol,ispinor,ispinor),im=1,ndim)
1042 : write(message,'(a,2x,500(e25.17e3,2x))') &
1043 5891 : & "#dc ",((((self%hdc%matlu(iatom)%mat(im+(ispinor-1)*ndim,im1+(ispinor1-1)*ndim,isppol),&
1044 6440 : & im=1,ndim),im1=1,ndim),ispinor=1,nspinor),ispinor1=1,nspinor)
1045 179 : call wrtout(unitselffunc_arr(iall),message,'COLL')
1046 179 : if (prtself) then
1047 0 : call wrtout(unitselffunc_arr2(iall),message,'COLL')
1048 : end if
1049 179 : if (self%has_moments == 1) then
1050 0 : do i=1,self%nmoments
1051 0 : write(tag_is,'(i1)') i
1052 0 : write(message,'(a,2x,500(e25.17e3,2x))') "#moments_"//trim(tag_is),&
1053 0 : & ((self%moments(i)%matlu(iatom)%mat(im,im1,isppol),im=1,nspinor*ndim),im1=1,nspinor*ndim)
1054 0 : call wrtout(unitselffunc_arr(iall),message,'COLL')
1055 0 : if (prtself) then
1056 0 : call wrtout(unitselffunc_arr2(iall),message,'COLL')
1057 : end if
1058 : end do ! i
1059 : end if ! moments
1060 61 : else if (optrw == 1 .and. iexist2 == 1 .and. ioerr == 0 .and. readimagonly == 0) then
1061 : !write(std_out,*) "8"
1062 : read(unitselffunc_arr(iall),*,iostat=ioerr) &
1063 775 : & chtemp,((((s_r(im+(ispinor-1)*ndim,im1+(ispinor1-1)*ndim),s_i(im+(ispinor-1)*ndim,im1+(ispinor1-1)*ndim),&
1064 862 : & im=1,ndim),im1=1,ndim),ispinor=1,nspinor),ispinor1=1,nspinor)
1065 : !if(ioerr<0) then
1066 : ! write(std_out,*)" HDC IOERR<",ioerr
1067 : !else if(ioerr>0) then
1068 : ! write(std_out,*)" HDC IOERR>",ioerr
1069 : !endif
1070 775 : self%hdc%matlu(iatom)%mat(:,:,isppol) = cmplx(s_r(:,:),s_i(:,:),kind=dp)
1071 :
1072 29 : if (self%has_moments == 1) then
1073 0 : do i=1,self%nmoments
1074 : read(unitselffunc_arr(iall),*,iostat=ioerr) &
1075 0 : & chtemp,((s_r(im,im1),s_i(im,im1),im=1,nspinor*ndim),im1=1,nspinor*ndim)
1076 0 : self%moments(i)%matlu(iatom)%mat(:,:,isppol) = cmplx(s_r(:,:),s_i(:,:),kind=dp)
1077 : end do ! i
1078 : end if ! moments
1079 : !write(6,*) "read selfhdc",self%hdc%matlu(1)%mat(1,1,1,1,1)
1080 32 : else if (readimagonly == 1 .and. (.not. present(opt_hdc))) then
1081 0 : self%hdc%matlu(iatom)%mat(:,:,isppol) = czero
1082 : !else
1083 : ! write(std_out,*) " self%hdc fixed in kramerskronig_self"
1084 : end if ! optrw
1085 240 : close(unitselffunc_arr(iall))
1086 390 : if (optrw == 2 .and. prtself) close(unitselffunc_arr2(iall))
1087 : ! enddo ! ispinor
1088 : end do ! isppol
1089 150 : ABI_FREE(s_r)
1090 532 : ABI_FREE(s_i)
1091 : end do ! iatom
1092 :
1093 130 : ABI_FREE(unitselffunc_arr)
1094 130 : ABI_SFREE(unitselffunc_arr2)
1095 : end if ! optrw==1 or 2 and myproc==master
1096 :
1097 : ! ===========================
1098 : ! == Error messages
1099 : ! ===========================
1100 289 : if (optrw == 1) then
1101 : ! call xmpi_barrier(spacecomm)
1102 : !write(std_out,*) ncount,maxval(pawtab(:)%lpawu)*2+1
1103 90 : call xmpi_bcast(iexist2,master,spacecomm,ier)
1104 90 : call xmpi_bcast(ioerr,master,spacecomm,ier)
1105 90 : if (iexist2 == 0 .or. ioerr /= 0) then
1106 32 : message = "Self-energy file does not exist or is incomplete"
1107 32 : if (readimagonly == 1 .or. present(opt_stop)) then
1108 0 : if (readimagonly == 1) message = "Self-energy file does not exist or is incomplete: check the number of self-energy data in file"
1109 0 : ABI_ERROR(message)
1110 32 : else if (paw_dmft%ireadself == 1 .or. paw_dmft%idmftloop > 0) then
1111 0 : ABI_WARNING(message)
1112 : end if ! readimagonly=1 or present(opt_stop)
1113 32 : if (iexist2 == 0 .and. (paw_dmft%ireadself == 1 .or. paw_dmft%idmftloop > 0)) then
1114 0 : write(message,'(4x,2a)') "File does not exist"
1115 0 : call wrtout(std_out,message,'COLL')
1116 : end if ! iexist2=0
1117 32 : if (ioerr < 0) then
1118 0 : write(message,'(4x,2a)') "End of file reached"
1119 0 : call wrtout(std_out,message,'COLL')
1120 : end if ! ioerr<0
1121 32 : if (ioerr > 0) then
1122 0 : write(message,'(4x,2a)') "Error during read statement"
1123 0 : call wrtout(std_out,message,'COLL')
1124 : end if ! ioerr>0
1125 : !if (paw_dmft%dmft_solv /= 4) then
1126 32 : write(message,'(4x,2a,5i5,2x,e14.7)') "-> Set Self-Energy Equal to double counting term"
1127 : !else if (paw_dmft%dmft_solv == 4) then
1128 : ! write(message,'(4x,a,a,5i5,2x,e14.7)') "-> Put Self-Energy Equal to dc term - shift"
1129 : ! call wrtout(std_out,message,'COLL')
1130 : ! write(message,'(4x,a,a,5i5,2x,e14.7)') " No self energy is given, change dmft_rslf"
1131 : ! ABI_ERROR(message)
1132 : !end if
1133 32 : call wrtout(std_out,message,'COLL')
1134 17513 : do ifreq=1,self%nw
1135 : ! write(std_out,*) "before",self%oper(1)%matlu(1)%mat(1,1,1,1,1)
1136 : ! write(std_out,*) "before",self%hdc%matlu(1)%mat(1,1,1,1,1)
1137 17481 : call copy_matlu(self%hdc%matlu(:),self%oper(ifreq)%matlu(:),natom)
1138 17513 : if (nspinor == 1 .and. nsppol == 2 .and. paw_dmft%dmft_dc >= 5) then
1139 32162 : do iatom=1,natom
1140 16081 : lpawu = paw_dmft%lpawu(iatom)
1141 16081 : if (lpawu == -1) cycle
1142 16081 : ndim = 2*lpawu + 1
1143 112567 : do im=1,ndim
1144 : self%oper(ifreq)%matlu(iatom)%mat(im,im,1) = self%oper(ifreq)%matlu(iatom)%mat(im,im,2) + &
1145 96486 : & paw_dmft%dmft_shiftself(iatom)
1146 : end do ! im
1147 : end do ! iatom
1148 : end if ! nspinor=1 and nsppol=2
1149 : ! write(std_out,*) "after",self%oper(1)%matlu(1)%mat(1,1,1,1,1)
1150 : ! write(std_out,*) "before",self%hdc%matlu(1)%mat(1,1,1,1,1)
1151 : !if (paw_dmft%dmft_solv == 4) then
1152 : ! if(ifreq==1) write(std_out,*) "shift",self%qmc_shift(1)
1153 : ! call shift_matlu(self%oper(ifreq)%matlu(:),natom,cmplx(self%qmc_shift(:),zero,kind=dp),1)
1154 : ! if(ifreq==1) write(std_out,*) "self after dc and shift",self%oper(ifreq)%matlu(1)%mat(1,1,1,1,1)
1155 : ! if(ifreq==1) write(std_out,*) "shift",self%qmc_shift(1)
1156 : !end if ! dmft_solv=4
1157 : end do ! ifreq
1158 32 : if (self%has_moments == 1) then
1159 0 : call copy_matlu(self%oper(1)%matlu(:),self%moments(1)%matlu(:),natom)
1160 : end if
1161 : else ! test read successfull
1162 : ! call xmpi_barrier(spacecomm)
1163 : !! lignes 924-928 semblent inutiles puisque la valeur de paw_dmft%fermie creee
1164 : !! en ligne 927 est ecrasee en ligne 992. BA+jmb
1165 : !! ABI_MALLOC(fermie_read2,(1))
1166 : !! fermie_read2(1)=fermie_read
1167 : !! call xmpi_sum(fermie_read2,spacecomm ,ier)
1168 : !! paw_dmft%fermie=fermie_read2(1)
1169 : !! ABI_FREE(fermie_read2)
1170 : !ncount = natom * nsppol * (nspinor**2) * (self%nw+1) *(maxval(paw_dmft%lpawu(:))*2+1)**2
1171 58 : ncount = 0
1172 279 : do iatom=1,natom
1173 221 : lpawu = paw_dmft%lpawu(iatom)
1174 221 : if (lpawu == -1) cycle
1175 59 : ndim = 2*lpawu + 1
1176 279 : ncount = ncount + ndim**2
1177 : end do ! iatom
1178 58 : ncount = ncount * (nspinor**2) * (self%nw+self%nmoments+1) * nsppol
1179 :
1180 : ! ===========================
1181 : ! bcast to other proc
1182 : ! ===========================
1183 174 : ABI_MALLOC(buffer,(ncount))
1184 : !! BA+jmb
1185 : !fermie_read2 = zero
1186 : !write(std_out,*) self%nw
1187 58 : if (myproc == master) then
1188 :
1189 : ! == Send read data to all process
1190 22 : if (readimagonly == 0) paw_dmft%fermie = fermie_read
1191 22 : icount = 0
1192 : ! Self energy-----------
1193 6020 : do ifreq=1,self%nw
1194 17058 : do iatom=1,natom
1195 11038 : lpawu = paw_dmft%lpawu(iatom)
1196 11038 : if (lpawu == -1) cycle
1197 6030 : ndim = (2*lpawu+1) * nspinor
1198 18490 : do isppol=1,nsppol
1199 47694 : do im1=1,ndim
1200 178928 : buffer(icount+1:icount+ndim) = self%oper(ifreq)%matlu(iatom)%mat(:,im1,isppol)
1201 36656 : icount = icount + ndim
1202 : !if (icount > ncount) then
1203 : ! write(message,'(2a,2i5)') ch10,"Error buffer",icount,ncount
1204 : ! iexit = 1
1205 : ! ABI_ERROR(message)
1206 : !end if ! icount > ncount
1207 : end do ! im1
1208 : end do ! isppol
1209 : end do ! iatom
1210 : end do ! ifreq
1211 : ! Double counting-------
1212 99 : do iatom=1,natom
1213 77 : lpawu = paw_dmft%lpawu(iatom)
1214 77 : if (lpawu == -1) cycle
1215 23 : ndim = (2*lpawu+1) * nspinor
1216 75 : do isppol=1,nsppol
1217 237 : do im1=1,ndim
1218 776 : buffer(icount+1:icount+ndim) = self%hdc%matlu(iatom)%mat(:,im1,isppol)
1219 160 : icount = icount + ndim
1220 : !if (icount > ncount) then
1221 : ! write(message,'(2a,2i5)') ch10,"Error buffer",icount,ncount
1222 : ! iexit = 1
1223 : ! ABI_ERROR(message)
1224 : !end if ! icount > ncount
1225 : end do ! im1
1226 : end do ! isppol
1227 : end do ! iatom
1228 22 : if (self%has_moments == 1) then
1229 0 : do i=1,self%nmoments
1230 0 : do iatom=1,natom
1231 0 : lpawu = paw_dmft%lpawu(iatom)
1232 0 : if (lpawu == -1) cycle
1233 0 : ndim = (2*lpawu+1) * nspinor
1234 0 : do isppol=1,nsppol
1235 0 : do im1=1,ndim
1236 0 : buffer(icount+1:icount+ndim) = self%moments(i)%matlu(iatom)%mat(:,im1,isppol)
1237 0 : icount = icount + ndim
1238 : end do ! im1
1239 : end do ! isppol
1240 : end do ! iatom
1241 : end do ! i
1242 : end if ! moments
1243 : end if ! proc=master
1244 58 : call xmpi_bcast(buffer(:),master,spacecomm,ier)
1245 : ! call xmpi_sum(iexit,spacecomm ,ier)
1246 : !!JB call xmpi_barrier(spacecomm)
1247 : !call xmpi_sum(buffer,spacecomm,ier)
1248 : !!JB call xmpi_barrier(spacecomm)
1249 :
1250 : ! bcast fermi level
1251 : !call xmpi_sum(fermie_read2,spacecomm,ier)
1252 58 : if (readimagonly == 0) then
1253 54 : call xmpi_bcast(paw_dmft%fermie,master,spacecomm,ier)
1254 : end if
1255 :
1256 58 : if (ier /= 0) then
1257 0 : message = "error in xmpi_sum in rw_self"
1258 0 : ABI_ERROR(message)
1259 : end if
1260 : !paw_dmft%fermie = fermie_read2(1)
1261 : ! write(std_out,*) "Fermi level",paw_dmft%fermie
1262 58 : icount = 0
1263 : ! Self ---------------
1264 21770 : do ifreq=1,self%nw
1265 59322 : do iatom=1,natom
1266 37552 : lpawu = paw_dmft%lpawu(iatom)
1267 37552 : if (lpawu == -1) cycle
1268 21744 : ndim = (2*lpawu+1) * nspinor
1269 65632 : do isppol=1,nsppol
1270 163632 : do im1=1,ndim
1271 608768 : self%oper(ifreq)%matlu(iatom)%mat(:,im1,isppol) = buffer(icount+1:icount+ndim)
1272 126080 : icount = icount + ndim
1273 : !write(6,*)'self procs', ifreq, self%oper(ifreq)%matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1)
1274 : end do ! im1
1275 : end do ! isppol
1276 : end do ! iatom
1277 : end do ! ifreq
1278 : ! hdc ---------------
1279 279 : do iatom=1,natom
1280 221 : lpawu = paw_dmft%lpawu(iatom)
1281 221 : if (lpawu == -1) cycle
1282 59 : ndim = (2*lpawu+1) * nspinor
1283 183 : do isppol=1,nsppol
1284 549 : do im1=1,ndim
1285 1448 : self%hdc%matlu(iatom)%mat(:,im1,isppol) = buffer(icount+1:icount+ndim)
1286 328 : icount = icount + ndim
1287 : end do ! im1
1288 : end do ! isppol
1289 : end do ! iatom
1290 58 : if (self%has_moments == 1) then
1291 0 : do i=1,self%nmoments
1292 0 : do iatom=1,natom
1293 0 : lpawu = paw_dmft%lpawu(iatom)
1294 0 : if (lpawu == -1) cycle
1295 0 : ndim = (2*lpawu+1) * nspinor
1296 0 : do isppol=1,nsppol
1297 0 : do im1=1,ndim
1298 0 : self%moments(i)%matlu(iatom)%mat(:,im1,isppol) = buffer(icount+1:icount+ndim)
1299 0 : icount = icount + ndim
1300 : end do ! im1
1301 : end do ! isppol
1302 : end do ! iatom
1303 : end do ! i
1304 : end if ! moments
1305 : !ABI_FREE(fermie_read2)
1306 58 : ABI_FREE(buffer)
1307 : end if ! test read successful
1308 : end if ! optrw==1
1309 :
1310 406 : if (optmaxent > 0 .and. optrw == 1 .and. iexist2 == 1 .and. ioerr == 0 .and. readimagonly == 1) then
1311 :
1312 4 : write(message,'(4x,2a)') " Rotate Back self-energy in the cubic basis"
1313 4 : call wrtout(std_out,message,'COLL')
1314 :
1315 : ! Kramers Kronig
1316 : !-----------------------------
1317 :
1318 4 : call kramerskronig_self(self,opt_selflimit(:),opt_hdc(:),paw_dmft,paw_dmft%filapp)
1319 :
1320 4 : call xmpi_matlu(eigvectmatlu(:),natom,spacecomm,master=master,option=2)
1321 :
1322 4 : call copy_matlu(opt_hdc(:),self%hdc%matlu(:),natom)
1323 :
1324 : ! Rotate back self
1325 : !-----------------------------
1326 :
1327 12004 : do ifreq=1,self%nw
1328 12000 : if (ifreq < 20) then
1329 76 : write(tag_freq,'(i5)') ifreq
1330 76 : write(message,'(a,2x,2a)') ch10," == Print Rotated Self Energy on real axis for frequency ",adjustl(tag_freq)
1331 76 : call wrtout(std_out,message,'COLL')
1332 76 : call print_matlu(self%oper(ifreq)%matlu(:),natom,1,compl=1)
1333 : end if ! ifreq<20
1334 12000 : if (self%distrib%procf(ifreq) /= myproc) cycle
1335 12004 : call rotate_matlu(self%oper(ifreq)%matlu(:),eigvectmatlu(:),natom,-1)
1336 : end do ! ifreq
1337 4 : call gather_oper(self%oper(:),self%distrib,paw_dmft,opt_ksloc=2)
1338 80 : do ifreq=1,min(19,self%nw)
1339 76 : write(tag_freq,'(i5)') ifreq
1340 76 : write(message,'(a,2x,2a)') ch10," == Print Self Energy rotated back in cubic basis on real axis for frequency ",adjustl(tag_freq)
1341 76 : call wrtout(std_out,message,'COLL')
1342 80 : call print_matlu(self%oper(ifreq)%matlu(:),natom,1,compl=1)
1343 : end do ! ifreq
1344 : end if ! rotate back self-energy
1345 :
1346 : ! call xmpi_barrier(spacecomm)
1347 : !write(std_out,*) "9"
1348 : ! - For the Tentative rotation of the self-energy file (begin destroy)
1349 406 : if (optmaxent > 0) then
1350 :
1351 110 : if (optrw == 2) then
1352 106 : call destroy_oper(energy_level)
1353 106 : call destroy_matlu(level_diag(:),natom)
1354 106 : if (self%has_moments == 1) then
1355 0 : do i=1,self%nmoments
1356 0 : call destroy_matlu(selfmomrot(:,i),natom)
1357 : end do
1358 0 : ABI_FREE(selfmomrot)
1359 : end if
1360 23195 : do ifreq=1,self%nw
1361 23195 : call destroy_oper(selfrotmatlu(ifreq))
1362 : end do ! ifreq
1363 476 : ABI_FREE(level_diag)
1364 23195 : ABI_FREE(selfrotmatlu)
1365 106 : ABI_FREE(unitselfrot) ! 7 is the max ndim possible
1366 : end if ! optrw=2
1367 :
1368 110 : call destroy_matlu(eigvectmatlu(:),natom)
1369 484 : ABI_FREE(eigvectmatlu)
1370 : end if ! optmaxent > 0
1371 : ! - For the Tentative rotation of the self-energy file (end destroy)
1372 :
1373 :
1374 : ! call flush_unit(std_out)
1375 : ! ABI_ERROR("Aboring now")
1376 406 : if (optrw == 0) then
1377 11 : if (paw_dmft%dmft_rslf == 0) then
1378 : !if (paw_dmft%dmft_solv /= 4) then
1379 11 : write(message,'(4x,a)') "-> Set Self-Energy Equal to double counting term"
1380 : !else if (paw_dmft%dmft_solv == 4) then
1381 : ! write(message,'(4x,a,a,5i5,2x,e14.7)') "-> Put Self-Energy Equal to dc term - shift"
1382 : !end if ! dmft_solv=4
1383 0 : else if (paw_dmft%dmft_rslf == -1) then
1384 0 : write(message,'(4x,a)') "-> Set Self-Energy Equal to zero"
1385 : end if ! dmft_rslf=0
1386 11 : call wrtout(std_out,message,'COLL')
1387 415 : do ifreq=1,self%nw
1388 415 : if (paw_dmft%dmft_rslf == 0) then
1389 404 : call copy_matlu(self%hdc%matlu(:),self%oper(ifreq)%matlu(:),natom)
1390 : ! if(ifreq==1) write(std_out,*) "self after dc",self%oper(ifreq)%matlu(1)%mat(1,1,1,1,1)
1391 : !if (paw_dmft%dmft_solv == 4) then
1392 : ! if(ifreq==1) write(std_out,*) "shift",self%qmc_shift(1)
1393 : ! call shift_matlu(self%oper(ifreq)%matlu(:),natom,cmplx(self%qmc_shift(:),zero,kind=dp),1)
1394 : ! if(ifreq==1) write(std_out,*) "self after dc and shift",self%oper(ifreq)%matlu(1)%mat(1,1,1,1,1)
1395 : ! if(ifreq==1) write(std_out,*) "shift",self%qmc_shift(1)
1396 : !end if ! dmft_solv=4
1397 0 : else if (paw_dmft%dmft_rslf == -1) then
1398 0 : call zero_matlu(self%oper(ifreq)%matlu(:),natom)
1399 : end if ! dmft_rslf
1400 : end do ! ifreq
1401 11 : if (paw_dmft%dmft_rslf == 0 .and. self%has_moments == 1) then
1402 0 : call copy_matlu(self%hdc%matlu(:),self%moments(1)%matlu(:),natom)
1403 : end if
1404 : end if ! optrw=0
1405 :
1406 : ! write(std_out,*) "optrw,use_fixed_self,istep,iter,istep_imp,iter_imp"
1407 : ! write(std_out,*) optrw,paw_dmft%use_fixed_self,istep,iter,istep_imp,iter_imp
1408 406 : if ((optrw == 1 .or. optrw == 3) .and. paw_dmft%use_fixed_self > 0 .and. istep <= istep_imp .and. iter <= iter_imp) then
1409 0 : write(message,'(4x,a)') "-> Set Self-Energy Equal to imposed self-energy"
1410 0 : call wrtout(std_out,message,'COLL')
1411 0 : do ifreq=1,self%nw
1412 0 : iatu = 0
1413 0 : do iatom=1,natom
1414 0 : lpawu = paw_dmft%lpawu(iatom)
1415 0 : if (lpawu == -1) cycle
1416 0 : self%oper(ifreq)%matlu(iatom)%mat(:,:,:) = czero
1417 0 : iatu = iatu + 1
1418 0 : ndim = 2*lpawu + 1
1419 0 : do isppol=1,nsppol
1420 0 : do ispinor=1,nspinor
1421 0 : do im=1,ndim
1422 0 : if (nspinor == 1) then
1423 0 : self%oper(ifreq)%matlu(iatom)%mat(im+(ispinor-1)*ndim,im+(ispinor-1)*ndim,isppol) = paw_dmft%fixed_self(im,im,isppol,iatu)
1424 : ! write(std_out,*) paw_dmft%fixed_self(im,im,isppol,iatu)
1425 : else
1426 0 : self%oper(ifreq)%matlu(iatom)%mat(im+(ispinor-1)*ndim,im+(ispinor-1)*ndim,isppol) = paw_dmft%fixed_self(im,im,ispinor,iatu)
1427 : ! write(message,'(a,i4,i4,2x,e20.10)') " Fixed self not implemented for nspinor==2"
1428 : ! call wrtout(std_out, message,'COLL')
1429 : ! ABI_ERROR("Aboring now")
1430 : end if ! nspinor
1431 : end do ! im
1432 : end do ! ispinor
1433 : end do ! isppol
1434 : end do ! iatom
1435 : end do ! ifreq
1436 0 : if (self%has_moments == 1) then
1437 0 : call copy_matlu(self%oper(1)%matlu(:),self%moments(1)%matlu(:),natom)
1438 0 : do i=2,self%nmoments
1439 0 : call zero_matlu(self%moments(i)%matlu(:),natom)
1440 : end do ! i
1441 : end if
1442 :
1443 : end if ! use_fixed_self
1444 : ABI_NVTX_END_RANGE()
1445 :
1446 406 : end subroutine rw_self
1447 : !!***
1448 :
1449 : !!****f* m_self/new_self
1450 : !! NAME
1451 : !! new_self
1452 : !!
1453 : !! FUNCTION
1454 : !!
1455 : !! Mix Old and New self_energy with the mixing coefficient dmft_mxsf
1456 : !!
1457 : !! INPUTS
1458 : !! self <type(self_type)>= variables related to self-energy
1459 : !! self_new <type(self_type)>= variables related to the new self-energy
1460 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
1461 : !!
1462 : !! OUTPUT
1463 : !! self <type(self_type)>= variables related to mixed self-energy
1464 : !!
1465 : !! SOURCE
1466 :
1467 106 : subroutine new_self(self,self_new,paw_dmft)
1468 :
1469 : !Arguments ------------------------------------
1470 : !type
1471 : ! type(crystal_t),intent(in) :: cryst_struc
1472 : type(self_type), intent(inout) :: self
1473 : type(self_type), intent(in) :: self_new
1474 : type(paw_dmft_type), intent(in) :: paw_dmft
1475 : !Local variables-------------------------------
1476 : integer :: i,iatom,icount,ifreq,im,im1,isppol,lpawu,natom,ndim,nspinor,nsppol
1477 : real(dp) :: alpha,diff_self,sum_self
1478 : character(len=500) :: message
1479 : ! *********************************************************************
1480 :
1481 106 : alpha = paw_dmft%dmft_mxsf
1482 106 : natom = paw_dmft%natom
1483 106 : nspinor = paw_dmft%nspinor
1484 106 : nsppol = paw_dmft%nsppol
1485 :
1486 23195 : do ifreq=1,self%nw
1487 67068 : do iatom=1,natom
1488 43873 : lpawu = paw_dmft%lpawu(iatom)
1489 43873 : if (lpawu == -1) cycle
1490 : self%oper(ifreq)%matlu(iatom)%mat(:,:,:) = (one-alpha)*self%oper(ifreq)%matlu(iatom)%mat(:,:,:) + &
1491 1355756 : & alpha*self_new%oper(ifreq)%matlu(iatom)%mat(:,:,:)
1492 : ! warning: self_new is the recent self-energy, which is mixed with self
1493 : ! to give self= mixed self energy. self_new is deallocated just after.
1494 : end do ! iatom
1495 : end do ! ifreq
1496 :
1497 106 : diff_self = zero
1498 106 : sum_self = zero
1499 106 : icount = 0
1500 :
1501 476 : do iatom=1,natom
1502 370 : lpawu = paw_dmft%lpawu(iatom)
1503 370 : if (lpawu == -1) cycle
1504 130 : ndim = nspinor * (2*lpawu+1)
1505 130 : icount = icount + nsppol*ndim
1506 440 : do isppol=1,nsppol
1507 1540 : do im1=1,ndim
1508 6180 : do im=1,ndim
1509 5010 : diff_self = diff_self + abs(self%hdc%matlu(iatom)%mat(im,im1,isppol)-self_new%hdc%matlu(iatom)%mat(im,im1,isppol))
1510 5010 : sum_self = sum_self + abs(self%hdc%matlu(iatom)%mat(im,im1,isppol))
1511 : self%hdc%matlu(iatom)%mat(im,im1,isppol) = (one-alpha)*self%hdc%matlu(iatom)%mat(im,im1,isppol) + &
1512 5976 : & alpha*self_new%hdc%matlu(iatom)%mat(im,im1,isppol)
1513 : end do ! im
1514 : end do ! im1
1515 : end do ! isppol
1516 : end do ! iatom
1517 :
1518 106 : if (self%has_moments == 1) then
1519 0 : do i=1,self%nmoments
1520 0 : do iatom=1,natom
1521 0 : lpawu = paw_dmft%lpawu(iatom)
1522 0 : if (lpawu == -1) cycle
1523 : self%moments(i)%matlu(iatom)%mat(:,:,:) = (one-alpha)*self%moments(i)%matlu(iatom)%mat(:,:,:) + &
1524 0 : & alpha*self_new%moments(i)%matlu(iatom)%mat(:,:,:)
1525 : end do ! iatom
1526 : end do ! i
1527 : end if ! moments
1528 :
1529 : !if(opt_mix==1) then
1530 : !endif
1531 :
1532 106 : diff_self = diff_self / dble(icount)
1533 :
1534 106 : write(message,'(8x,a,e12.5)') "DMFT Loop: Precision on self-energy is",diff_self
1535 106 : call wrtout(std_out,message,'COLL')
1536 106 : if (diff_self < paw_dmft%dmft_fermi_prec .and. sum_self > tol6 .and. paw_dmft%idmftloop >= 2) then
1537 0 : write(message,'(a,8x,a,e9.2,a,8x,a)') ch10, "Change of self =<", paw_dmft%dmft_fermi_prec,&
1538 0 : & ch10,"DMFT Loop: Self Energy is converged"
1539 0 : call wrtout(std_out,message,'COLL')
1540 0 : self%iself_cv = 1
1541 : else
1542 106 : write(message,'(a,8x,a)') ch10,"DMFT Loop: Self Energy is not converged"
1543 106 : call wrtout(std_out,message,'COLL')
1544 106 : self%iself_cv = 0
1545 : end if ! diff_self
1546 :
1547 106 : end subroutine new_self
1548 : !!***
1549 :
1550 : !!****f* m_self/make_qmcshift_self
1551 : !! NAME
1552 : !! make_qmcshift_hu
1553 : !!
1554 : !! FUNCTION
1555 : !!
1556 : !! INPUTS
1557 : !! hu <type(hu_type)> = U interaction
1558 : !! paw_dmft <type(paw_dmft_type)> = paw+dmft related data
1559 : !!
1560 : !! OUTPUT
1561 : !! self%qmc_shift in self <type(self_type)> = Self-energy
1562 : !!
1563 : !! SOURCE
1564 :
1565 : !subroutine make_qmcshift_self(cryst_struc,hu,self,apply)
1566 :
1567 : !Arguments ------------------------------------
1568 : !type
1569 : ! type(crystal_t),intent(in) :: cryst_struc
1570 : ! type(hu_type),intent(in) :: hu(cryst_struc%ntypat)
1571 : ! type(self_type),intent(inout) :: self
1572 : ! logical, optional :: apply
1573 :
1574 : !Local variables-------------------------------
1575 : ! integer :: im,iatom,ifreq,itypat,lpawu,tndim
1576 : ! real(dp) :: hu_shift2
1577 : ! character(len=500) :: message
1578 : ! *********************************************************************
1579 :
1580 : ! do iatom = 1 , cryst_struc%natom
1581 : ! lpawu=self%hdc%matlu(iatom)%lpawu
1582 : ! tndim=2*lpawu+1
1583 : ! itypat=cryst_struc%typat(iatom)
1584 : ! if(lpawu/=-1) then
1585 : ! self%qmc_shift(iatom) = zero
1586 : ! do im =1, 2*tndim-1
1587 : !! write(std_out,*)"make before",self%qmc_shift(iatom)
1588 : ! self%qmc_shift(iatom) = self%qmc_shift(iatom) + hu(itypat)%uqmc(im)
1589 : !! write(std_out,*)"make after",self%qmc_shift(iatom)
1590 : ! enddo
1591 : ! self%qmc_shift(iatom) = self%qmc_shift(iatom) / two
1592 : ! hu_shift2 = hu(itypat)%uqmc(1)
1593 : !
1594 : ! do im = 2*tndim, 2*tndim + 2*tndim -3
1595 : ! hu_shift2 = hu_shift2 + hu(itypat)%uqmc(im)
1596 : ! enddo
1597 :
1598 : ! hu_shift2 = hu_shift2 / two
1599 : ! write(message,'(2a,i4)') ch10,' -------> For Correlated atom',iatom
1600 : ! call wrtout(std_out, message,'COLL')
1601 :
1602 : ! if(abs(self%qmc_shift(iatom)-hu_shift2)>tol6) then
1603 : ! write(message,'(2a,2f16.7)') " Shift for QMC is not correctly"&
1604 : !& ," computed",self%qmc_shift(iatom),hu_shift2
1605 : ! ABI_ERROR(message)
1606 : ! endif ! shifts not equals
1607 : !
1608 : ! write(message,'(4x,a,f16.7)') &
1609 : !& " Shift for QMC (used to compute G(w)) is (in Ha) :",&
1610 : !& self%qmc_shift(iatom)
1611 : ! call wrtout(std_out, message,'COLL')
1612 :
1613 : ! self%qmc_xmu(iatom)=-self%qmc_shift(iatom)
1614 : ! self%qmc_xmu(iatom)=zero
1615 : ! write(message,'(4x,a,f16.7)') &
1616 : !& "Artificial Shift used in QMC AND to compute G is (in Ha) :",self%qmc_xmu(iatom)
1617 : ! ABI_WARNING(message)
1618 :
1619 : ! endif ! lpawu/=1
1620 : ! enddo ! natom
1621 :
1622 : ! if(present(apply)) then
1623 : ! if (apply) then
1624 : ! write(message,'(5x,a,f16.7,a)') " Shifts applied to self"
1625 : ! call wrtout(std_out, message,'COLL')
1626 : ! do ifreq=1,self%nw
1627 : ! call shift_matlu(self%oper(ifreq)%matlu,cryst_struc%natom,cmplx(self%qmc_shift,0.d0,kind=dp),1)
1628 : ! call shift_matlu(self%oper(ifreq)%matlu,cryst_struc%natom,cmplx(self%qmc_xmu,0.d0,kind=dp),1)
1629 : ! enddo
1630 : ! endif
1631 : ! endif
1632 :
1633 :
1634 : !end subroutine make_qmcshift_self
1635 : !!***
1636 :
1637 : !!****f* m_self/kramerskronig_self
1638 : !! NAME
1639 : !! kramerskronig_self
1640 : !!
1641 : !! FUNCTION
1642 : !! Compute the real part of the self-energy using Kramers-Kronig formula.
1643 : !!
1644 : !! INPUTS
1645 : !! self <type(self_type)>= variables related to self-energy
1646 : !! selflimit = 0th order moment of the self-energy for each atom
1647 : !! selfhdc = double counting
1648 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
1649 : !! filapp = name of the file
1650 : !!
1651 : !! OUTPUT
1652 : !!
1653 : !! SOURCE
1654 :
1655 4 : subroutine kramerskronig_self(self,selflimit,selfhdc,paw_dmft,filapp)
1656 :
1657 : !Arguments ------------------------------------
1658 : type(self_type), intent(inout) :: self
1659 : type(matlu_type), intent(in) :: selfhdc(self%hdc%natom),selflimit(self%hdc%natom)
1660 : type(paw_dmft_type), intent(in) :: paw_dmft
1661 : character(len=fnlen), intent(in) :: filapp
1662 : !Local variables-------------------------------
1663 : integer :: iatom,ifreq,im,im1,ispinor,ispinor1,isppol,jfreq
1664 : integer :: lpawu,myproc,natom,ndim,nspinor,nsppol,unt
1665 : character(len=2) :: tag_im,tag_im1
1666 : character(len=4) :: tag_at
1667 : character(len=50) :: tag_is
1668 : character(len=500) :: message
1669 : ! *********************************************************************
1670 :
1671 : !delta = 0.0000000
1672 :
1673 4 : myproc = paw_dmft%myproc
1674 4 : natom = self%hdc%natom
1675 4 : nspinor = self%hdc%nspinor
1676 4 : nsppol = self%hdc%nsppol
1677 :
1678 : !ABI_MALLOC(selftemp_imag,(self%nw))
1679 :
1680 4 : write(message,'(2a,i4)') ch10,' ------ High-frequency limit of the Self-Energy'
1681 4 : call wrtout(std_out,message,'COLL')
1682 :
1683 4 : call print_matlu(selflimit(:),natom,3)
1684 :
1685 : !print norms
1686 4 : write(message,'(2a,i4)') ch10,' ------ Double counting'
1687 4 : call wrtout(std_out,message,'COLL')
1688 :
1689 4 : call print_matlu(selfhdc(:),natom,3)
1690 :
1691 : ! Compute limit of Real Part and put in double counting energy.
1692 : ! call copy_matlu(selfhdc,self%hdc%matlu,natom)
1693 : !!write(6,*) "selfhdc kramerskronig",selfhdc(1)%mat(1,1,1,1,1)
1694 : !write(6,*) "selfr%hdc kramerskronig",self%hdc%matlu(1)%mat(1,1,1,1,1)
1695 :
1696 12004 : do ifreq=1,self%nw
1697 12000 : if (self%distrib%procf(ifreq) /= myproc) cycle
1698 9000000 : do jfreq=1,self%nw-1
1699 8997000 : if (jfreq == ifreq) cycle
1700 17991002 : do iatom=1,natom
1701 8994001 : lpawu = selfhdc(iatom)%lpawu
1702 8994001 : if (lpawu == -1) cycle
1703 : self%oper(ifreq)%matlu(iatom)%mat(:,:,:) = self%oper(ifreq)%matlu(iatom)%mat(:,:,:) - &
1704 : & cmplx(aimag(self%oper(jfreq)%matlu(iatom)%mat(:,:,:)) / (self%omega(ifreq)-self%omega(jfreq)) &
1705 296805032 : & * (self%omega(jfreq+1)-self%omega(jfreq)),zero,kind=dp)
1706 : end do ! iatom
1707 : end do ! jfreq
1708 6004 : do iatom=1,natom
1709 3000 : lpawu = selfhdc(iatom)%lpawu
1710 3000 : if (lpawu == -1) cycle
1711 : self%oper(ifreq)%matlu(iatom)%mat(:,:,:) = cmplx(dble(self%oper(ifreq)%matlu(iatom)%mat(:,:,:))/pi,&
1712 108000 : & aimag(self%oper(ifreq)%matlu(iatom)%mat(:,:,:)),kind=dp) + selflimit(iatom)%mat(:,:,:)
1713 : ! write(6,*) "TWO FACTOR IS PUT BECAUSE OF MAXENT CODE ??"
1714 : !self%oper(ifreq)%matlu(iatom)%mat(:,:,:) = self%oper(ifreq)%matlu(iatom)%mat(:,:,:) * half
1715 : ! & aimag(self%oper(ifreq)%matlu(iatom)%mat(:,:,:)),kind=dp) * half
1716 : end do ! iatom
1717 : end do ! ifreq
1718 :
1719 4 : call gather_oper(self%oper(:),self%distrib,paw_dmft,opt_ksloc=2)
1720 :
1721 4 : if (myproc == 0) then
1722 :
1723 1 : unt = get_unit()
1724 1 : open(unit=unt,file=trim(filapp)//"_DFTDMFT_Self_realaxis.dat",status='unknown',form='formatted')
1725 1 : rewind(unt)
1726 1 : write(unt,'(4a)') "# Self-energy on the real axis, with the real part computed using Kramers-Kronig relations.",ch10, &
1727 2 : "# Real Frequency (Ha.) Real part Imaginary part",ch10
1728 1 : tag_is = ""
1729 2 : do iatom=1,natom
1730 1 : lpawu = selfhdc(iatom)%lpawu
1731 1 : if (lpawu == -1) cycle
1732 1 : write(tag_at,'(i4)') iatom
1733 1 : ndim = 2*lpawu + 1
1734 3 : do isppol=1,nsppol
1735 1 : if (nsppol == 2) tag_is = trim(adjustl(" for spin " // merge("up ","down",isppol==1))) // "and"
1736 3 : do ispinor=1,nspinor
1737 3 : do ispinor1=1,nspinor
1738 7 : do im=1,ndim
1739 5 : write(tag_im,'(i2)') im + (ispinor-1)*ndim
1740 31 : do im1=1,ndim
1741 25 : write(tag_im1,'(i2)') im1 + (ispinor1-1)*ndim
1742 25 : write(unt,'(9a)') "## Sigma_{",trim(adjustl(tag_im)),",",trim(adjustl(tag_im1)),"}",trim(adjustl(tag_is))," for atom ",trim(adjustl(tag_at)),ch10
1743 75025 : do ifreq=1,self%nw
1744 75000 : write(unt,'(2x,393(es24.16e3,2x))') self%omega(ifreq),&
1745 75000 : & dble(self%oper(ifreq)%matlu(iatom)%mat(im+(ispinor-1)*ndim,im1+(ispinor1-1)*ndim,isppol)), &
1746 150025 : & aimag(self%oper(ifreq)%matlu(iatom)%mat(im+(ispinor-1)*ndim,im1+(ispinor1-1)*ndim,isppol))
1747 : end do ! ifreq
1748 30 : write(unt,*)
1749 : end do ! im1
1750 : end do ! im
1751 : end do ! ispinor1
1752 : end do ! ispinor
1753 : end do ! isppol
1754 : end do ! iatom
1755 :
1756 1 : close(unt)
1757 :
1758 : end if ! master node
1759 :
1760 : !do iatom=1,natom
1761 : ! lpawu = self%oper(1)%matlu(iatom)%lpawu
1762 : ! if (lpawu == -1) cycle
1763 : ! ndim = 2*lpawu + 1
1764 : ! do isppol=1,nsppol
1765 : ! do ispinor=1,nspinor
1766 : ! do ispinor1=1,nspinor
1767 : ! do im=1,ndim
1768 : ! do im1=1,ndim
1769 : !write(6,*)
1770 : !"realpart",real(selflimit(iatom)%mat(im,im1,isppol,ispinor,ispinor1))
1771 : ! do ifreq=1,self%nw
1772 : ! selftemp_re(ifreq)=zero
1773 : ! selftemp_imag(ifreq)=aimag(self%oper(ifreq)%matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1))
1774 : ! do jfreq=1,self%nw-1
1775 : ! if(jfreq==ifreq) cycle
1776 : ! selftemp_re(ifreq)=selftemp_re(ifreq) - &
1777 : ! &
1778 : ! aimag(self%oper(jfreq)%matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1)) &
1779 : ! & *(self%omega(ifreq)-self%omega(jfreq)) &
1780 : ! & /((self%omega(ifreq)-self%omega(jfreq))**2+delta**2)&
1781 : ! & *(self%omega(jfreq+1)-self%omega(jfreq))
1782 : ! selftemp_re(ifreq)=selftemp_re(ifreq) - &
1783 : !&
1784 : !aimag(self%oper(jfreq)%matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1)) &
1785 : !& /(self%omega(ifreq)-self%omega(jfreq)) *
1786 : !(self%omega(jfreq+1)-self%omega(jfreq))
1787 : ! enddo
1788 : ! selftemp_re(ifreq)=selftemp_re(ifreq)/pi
1789 : !write(671,*)
1790 : !self%omega(ifreq),selftemp_re(ifreq),selftemp_imag(ifreq)
1791 : ! enddo
1792 : ! TEST*************************
1793 : ! do ifreq=1,self%nw
1794 : ! selftemp_imag(ifreq)=zero
1795 : ! do jfreq=1,self%nw-1
1796 : ! if(jfreq==ifreq) cycle
1797 : !! selftemp_re(ifreq)=selftemp_re(ifreq) - &
1798 : !! &
1799 : !aimag(self%oper(jfreq)%matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1)) &
1800 : !! & *(self%omega(ifreq)-self%omega(jfreq)) &
1801 : !! & /((self%omega(ifreq)-self%omega(jfreq))**2+delta**2)&
1802 : !! & *(self%omega(jfreq+1)-self%omega(jfreq))
1803 : ! selftemp_imag(ifreq)=selftemp_imag(ifreq) + &
1804 : ! & selftemp_re(jfreq) &
1805 : ! & /(self%omega(ifreq)-self%omega(jfreq)) *
1806 : ! (self%omega(jfreq+1)-self%omega(jfreq))
1807 : ! enddo
1808 : ! selftemp_imag(ifreq)=selftemp_imag(ifreq)/pi
1809 : ! write(672,*)
1810 : ! self%omega(ifreq),selftemp_re(ifreq),selftemp_imag(ifreq)
1811 : ! enddo
1812 : ! TEST*************************
1813 : ! write(6,*) "TWO FACTOR IS PUT BECAUSE OF MAXENT CODE ??"
1814 : ! do ifreq=1,self%nw
1815 : ! write(68,*)
1816 : ! self%omega(ifreq),selftemp_re(ifreq),selftemp_imag(ifreq)
1817 : ! selftemp_re(ifreq)=selftemp_re(ifreq)+ &
1818 : !& real(selflimit(iatom)%mat(im,im1,isppol,ispinor,ispinor1)- &
1819 : !& selfhdc(iatom)%mat(im,im1,isppol,ispinor,ispinor1))
1820 : ! self%oper(ifreq)%matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1)&
1821 : ! &
1822 : ! =cmplx(selftemp_re(ifreq),selftemp_imag(ifreq),kind=dp)/two
1823 : !& =cmplx(0.d0,selftemp_imag(ifreq),kind=dp)/two
1824 : ! & =cmplx(selftemp_re(ifreq),0.d0,kind=dp)/two
1825 : ! self%oper(ifreq)%matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1)&
1826 : !& =cmplx(selftemp_re(ifreq),0.d0,kind=dp)/two
1827 : !& =cmplx(0.d0,0.d0,kind=dp)/two
1828 : ! The factor two is here to compensate for the factor two in OmegaMaxent..
1829 : ! & =cmplx(selftemp_re(ifreq),0.0,kind=dp)
1830 : ! write(67,*)
1831 : ! self%omega(ifreq),real(self%oper(ifreq)%matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1))&
1832 : ! ,aimag(self%oper(ifreq)%matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1))
1833 : ! enddo
1834 : ! write(67,*)
1835 : !write(68,*)
1836 : !!!!!!!!!! Z renormalization
1837 : ! i0=389
1838 : ! slope=(selftemp_re(i0+1)-selftemp_re(i0))/&
1839 : ! (self%omega(i0+1)-self%omega(i0))
1840 : ! y0= selftemp_re(i0)
1841 : ! do ifreq=1,self%nw
1842 : ! selftemp_re(ifreq)=slope * (self%omega(ifreq)-self%omega(i0))
1843 : ! + y0
1844 : ! selftemp_imag(ifreq)=zero
1845 : ! self%oper(ifreq)%matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1)&
1846 : ! &
1847 : ! =cmplx(selftemp_re(ifreq),selftemp_imag(ifreq),kind=dp)/two
1848 : ! write(6777,*)
1849 : ! self%omega(ifreq),real(self%oper(ifreq)%matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1)),aimag(self%oper(ifreq)%matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1))
1850 : ! enddo
1851 : !!!!!!!!!!
1852 : ! enddo
1853 : ! enddo
1854 : ! enddo
1855 : ! enddo
1856 : ! enddo ! isppol
1857 : ! end do ! iatom
1858 :
1859 4 : end subroutine kramerskronig_self
1860 : !!***
1861 :
1862 : !!****f* m_self/selfreal2imag_self
1863 : !! NAME
1864 : !! selfreal2imag_self
1865 : !!
1866 : !! FUNCTION
1867 : !! Print self on imaginary axis from real axis
1868 : !!
1869 : !! INPUTS
1870 : !! selfr = self on real axis
1871 : !! self = self on imaginary axis
1872 : !! filapp = name of the filename
1873 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
1874 : !!
1875 : !! OUTPUT
1876 : !! self%qmc_shift in self <type(self_type)> = Self-energy
1877 : !!
1878 : !! SOURCE
1879 :
1880 4 : subroutine selfreal2imag_self(selfr,self,filapp,paw_dmft)
1881 :
1882 : !Arguments ------------------------------------
1883 : !type
1884 : type(self_type), intent(in) :: selfr,self
1885 : character(len=fnlen), intent(in) :: filapp
1886 : type(paw_dmft_type), intent(in) :: paw_dmft
1887 : !Local variables-------------------------------
1888 : integer :: iatom,ifreq,im,im1,ispinor,ispinor1,isppol,jfreq
1889 : integer :: lpawu,myproc,natom,ndim,nspinor,nsppol,unt
1890 : logical :: triqs
1891 : complex(dp) :: omega
1892 : !real(dp) :: delta
1893 4 : type(self_type) :: selftempmatsub
1894 : character(len=2) :: tag_im,tag_im1
1895 : character(len=4) :: tag_at
1896 : character(len=50) :: tag_is
1897 4 : type(matlu_type), allocatable :: matlu_tmp(:)
1898 : ! *********************************************************************
1899 :
1900 : !delta=0.0000000
1901 4 : call initialize_self(selftempmatsub,paw_dmft)
1902 :
1903 4 : myproc = paw_dmft%myproc
1904 4 : natom = self%hdc%natom
1905 4 : nspinor = self%hdc%nspinor
1906 4 : nsppol = self%hdc%nsppol
1907 : ! Compute limit of Real Part and put in double counting energy.
1908 : ! call copy_matlu(selfhdc,self%hdc%matlu,natom)
1909 : !write(6,*) "self3",aimag(selfr%oper(489)%matlu(1)%mat(1,1,1,1,1))
1910 :
1911 4 : triqs = (paw_dmft%dmft_solv == 6 .or. paw_dmft%dmft_solv == 7)
1912 :
1913 4 : if (triqs) then
1914 0 : ABI_MALLOC(matlu_tmp,(natom))
1915 0 : call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu(:),matlu_tmp(:))
1916 : end if ! triqs
1917 :
1918 4996 : do ifreq=1,self%nw
1919 4992 : if (self%distrib%procf(ifreq) /= myproc) cycle
1920 1248 : omega = cmplx(zero,self%omega(ifreq),kind=dp)
1921 3744000 : do jfreq=1,selfr%nw-1
1922 7486752 : do iatom=1,natom
1923 3742752 : lpawu = paw_dmft%lpawu(iatom)
1924 3742752 : if (lpawu == -1) cycle
1925 : selftempmatsub%oper(ifreq)%matlu(iatom)%mat(:,:,:) = selftempmatsub%oper(ifreq)%matlu(iatom)%mat(:,:,:) - &
1926 : & aimag(selfr%oper(jfreq)%matlu(iatom)%mat(:,:,:)) / &
1927 123510816 : & (omega-selfr%omega(jfreq)) * (selfr%omega(jfreq+1)-selfr%omega(jfreq))
1928 : end do ! iatom
1929 : end do ! jfreq
1930 1248 : call fac_matlu(selftempmatsub%oper(ifreq)%matlu(:),natom,cone/pi)
1931 1252 : if (triqs) then ! add the missing high-frequency moment when using TRIQS
1932 0 : call add_matlu(selftempmatsub%oper(ifreq)%matlu(:),self%moments(1)%matlu(:),matlu_tmp(:),natom,1)
1933 0 : call copy_matlu(matlu_tmp(:),selftempmatsub%oper(ifreq)%matlu(:),natom)
1934 : end if ! triqs
1935 : end do ! ifreq
1936 :
1937 4 : if (triqs) then
1938 0 : call destroy_matlu(matlu_tmp(:),natom)
1939 0 : ABI_FREE(matlu_tmp)
1940 : end if ! triqs
1941 :
1942 4 : call gather_oper(selftempmatsub%oper(:),self%distrib,paw_dmft,opt_ksloc=2)
1943 :
1944 4 : if (myproc == 0) then
1945 1 : unt = get_unit()
1946 1 : open(unit=unt,file=trim(filapp)//"_DFTDMFT_Self_backtransform.dat",status='unknown',form='formatted')
1947 1 : write(unt,'(6a)') "# Hilbert transform of the analytically continued self-energy. To be compared with the actual",ch10, &
1948 2 : & "# self-energy on the imaginary axis in the cubic basis.",ch10,"# Matsubara Frequency (Ha.) Real part Imaginary part",ch10
1949 1 : tag_is = ""
1950 2 : do iatom=1,natom
1951 1 : lpawu = self%hdc%matlu(iatom)%lpawu
1952 1 : if (lpawu == -1) cycle
1953 1 : write(tag_at,'(i4)') iatom
1954 1 : ndim = 2*lpawu + 1
1955 3 : do isppol=1,nsppol
1956 1 : if (nsppol == 2) tag_is = trim(adjustl(" for spin " // merge("up ","down",isppol==1))) // "and"
1957 3 : do ispinor=1,nspinor
1958 3 : do ispinor1=1,nspinor
1959 7 : do im=1,ndim
1960 5 : write(tag_im,'(i2)') im + (ispinor-1)*ndim
1961 31 : do im1=1,ndim
1962 25 : write(tag_im1,'(i2)') im1 + (ispinor1-1)*ndim
1963 : !do jfreq=1,selfr%nw-1
1964 : ! write(6700,*) selfr%omega(jfreq),aimag(selfr%oper(jfreq)%matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1))
1965 : !enddo
1966 : ! write(6700,*)
1967 25 : write(unt,'(9a)') "## Sigma_{",trim(adjustl(tag_im)),",",trim(adjustl(tag_im1)),"}",trim(adjustl(tag_is))," for atom ",trim(adjustl(tag_at)),ch10
1968 31225 : do ifreq=1,self%nw
1969 31200 : write(unt,'(2x,393(es24.16e3,2x))') self%omega(ifreq),&
1970 31200 : & dble(selftempmatsub%oper(ifreq)%matlu(iatom)%mat(im+(ispinor-1)*ndim,im1+(ispinor1-1)*ndim,isppol)),&
1971 62425 : & aimag(selftempmatsub%oper(ifreq)%matlu(iatom)%mat(im+(ispinor-1)*ndim,im1+(ispinor1-1)*ndim,isppol))
1972 : end do ! ifreq
1973 30 : write(unt,*)
1974 : end do ! im1
1975 : end do ! im
1976 : end do ! ispinor1
1977 : end do ! ispinor
1978 : end do ! isppol
1979 : end do ! iatom
1980 1 : close(unt)
1981 : end if ! master node
1982 :
1983 4 : call destroy_self(selftempmatsub)
1984 :
1985 4 : end subroutine selfreal2imag_self
1986 : !!***
1987 :
1988 0 : END MODULE m_self
1989 : !!***
|
