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